version 3.10-dev
3p3c/primaryvariableswitch.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3//
4// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5// SPDX-License-Identifier: GPL-3.0-or-later
6//
13#ifndef DUMUX_3P3C_PRIMARY_VARIABLE_SWITCH_HH
14#define DUMUX_3P3C_PRIMARY_VARIABLE_SWITCH_HH
15
16#include <iostream>
17
18#include <dumux/porousmediumflow/compositional/primaryvariableswitch.hh>
19
20namespace Dumux {
21
26class ThreePThreeCPrimaryVariableSwitch
27: public PrimaryVariableSwitch<ThreePThreeCPrimaryVariableSwitch>
28{
29 using ParentType = PrimaryVariableSwitch<ThreePThreeCPrimaryVariableSwitch>;
30 friend ParentType;
31
32public:
33 using ParentType::ParentType;
34
35protected:
36
37 // perform variable switch at a degree of freedom location
38 template<class VolumeVariables, class GlobalPosition>
39 bool update_(typename VolumeVariables::PrimaryVariables& priVars,
40 const VolumeVariables& volVars,
41 std::size_t dofIdxGlobal,
42 const GlobalPosition& globalPos)
43 {
44 using PrimaryVariables = typename VolumeVariables::PrimaryVariables;
45 using Scalar = typename PrimaryVariables::value_type;
46 using Indices = typename VolumeVariables::Indices;
47 using FluidSystem = typename VolumeVariables::FluidSystem;
48
49 // evaluate if the primary variable switch would switch
50 bool wouldSwitch = false;
51 auto phasePresence = priVars.state();
52 auto newPhasePresence = phasePresence;
53
54 // check if a primary variable switch is necessary
55 if (phasePresence == Indices::threePhases)
56 {
57 Scalar Smin = 0;
58 if (this->wasSwitched_[dofIdxGlobal])
59 Smin = -0.01;
60
61 if (volVars.saturation(FluidSystem::gPhaseIdx) <= Smin)
62 {
63 wouldSwitch = true;
64 // gas phase disappears
65 if (this->verbosity() > 1)
66 std::cout << "Gas phase disappears at dof " << dofIdxGlobal
67 << ", coordinates: " << globalPos << ", sg: "
68 << volVars.saturation(FluidSystem::gPhaseIdx) << std::endl;
69 newPhasePresence = Indices::wnPhaseOnly;
70
71 priVars[Indices::switch1Idx] = volVars.moleFraction(FluidSystem::wPhaseIdx, FluidSystem::gCompIdx);
72 }
73 else if (volVars.saturation(FluidSystem::wPhaseIdx) <= Smin)
74 {
75 wouldSwitch = true;
76 // water phase disappears
77 if (this->verbosity() > 1)
78 std::cout << "Water phase disappears at dof " << dofIdxGlobal
79 << ", coordinates: " << globalPos << ", sw: "
80 << volVars.saturation(FluidSystem::wPhaseIdx) << std::endl;
81 newPhasePresence = Indices::gnPhaseOnly;
82
83 priVars[Indices::switch1Idx] = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::wCompIdx);
84 }
85 else if (volVars.saturation(FluidSystem::nPhaseIdx) <= Smin)
86 {
87 wouldSwitch = true;
88 // NAPL phase disappears
89 if (this->verbosity() > 1)
90 std::cout << "NAPL phase disappears at dof " << dofIdxGlobal
91 << ", coordinates: " << globalPos << ", sn: "
92 << volVars.saturation(FluidSystem::nPhaseIdx) << std::endl;
93 newPhasePresence = Indices::wgPhaseOnly;
94
95 priVars[Indices::switch2Idx] = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::nCompIdx);
96 }
97 }
98 else if (phasePresence == Indices::wPhaseOnly)
99 {
100 bool gasPresent = false;
101 bool nonwettingPresent = false;
102 // calculate fractions of the partial pressures in the
103 // hypothetical gas phase
104 Scalar xwg = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::wCompIdx);
105 Scalar xgg = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::gCompIdx);
106 Scalar xng = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::nCompIdx);
107 /* take care:
108 for xgg in case Indices::wPhaseOnly we compute xgg=henry_air*x2w
109 for xwg in case Indices::wPhaseOnly we compute xwg=pwsat
110 for xng in case Indices::wPhaseOnly we compute xng=henry_NAPL*x1w
111 */
112
113 Scalar xgMax = 1.0;
114 if (xwg + xgg + xng > xgMax)
115 wouldSwitch = true;
116 if (this->wasSwitched_[dofIdxGlobal])
117 xgMax *= 1.02;
118
119 // if the sum of the mole fractions would be larger than
120 // 100%, gas phase appears
121 if (xwg + xgg + xng > xgMax)
122 {
123 // gas phase appears
124 if (this->verbosity() > 1)
125 std::cout << "gas phase appears at dof " << dofIdxGlobal
126 << ", coordinates: " << globalPos << ", xwg + xgg + xng: "
127 << xwg + xgg + xng << std::endl;
128 gasPresent = true;
129 }
130
131 // calculate fractions in the hypothetical NAPL phase
132 Scalar xnn = volVars.moleFraction(FluidSystem::nPhaseIdx, FluidSystem::nCompIdx);
133 /* take care:
134 for xnn in case Indices::wPhaseOnly we compute xnn=henry_mesitylene*x1w,
135 where a hypothetical gas pressure is assumed for the Henry
136 x0n is set to NULL (all NAPL phase is dirty)
137 x2n is set to NULL (all NAPL phase is dirty)
138 */
139
140 Scalar xnMax = 1.0;
141 if (xnn > xnMax)
142 wouldSwitch = true;
143 if (this->wasSwitched_[dofIdxGlobal])
144 xnMax *= 1.02;
145
146 // if the sum of the hypothetical mole fractions would be larger than
147 // 100%, NAPL phase appears
148 if (xnn > xnMax)
149 {
150 // NAPL phase appears
151 if (this->verbosity() > 1)
152 std::cout << "NAPL phase appears at dof " << dofIdxGlobal
153 << ", coordinates: " << globalPos << ", xnn: "
154 << xnn << std::endl;
155 nonwettingPresent = true;
156 }
157
158 if (gasPresent && !nonwettingPresent)
159 {
160 newPhasePresence = Indices::wgPhaseOnly;
161 priVars[Indices::switch1Idx] = 0.9999;
162 priVars[Indices::switch2Idx] = 0.0001;
163 }
164 else if (gasPresent && nonwettingPresent)
165 {
166 newPhasePresence = Indices::threePhases;
167 priVars[Indices::switch1Idx] = 0.9999;
168 priVars[Indices::switch2Idx] = 0.0001;
169 }
170 else if (!gasPresent && nonwettingPresent)
171 {
172 newPhasePresence = Indices::wnPhaseOnly;
173 priVars[Indices::switch1Idx] = volVars.moleFraction(FluidSystem::wPhaseIdx, FluidSystem::gCompIdx);
174 priVars[Indices::switch2Idx] = 0.0001;
175 }
176 }
177 else if (phasePresence == Indices::gnPhaseOnly)
178 {
179 bool nonwettingPresent = false;
180 bool wettingPresent = false;
181
182 Scalar Smin = 0.0;
183 if (this->wasSwitched_[dofIdxGlobal])
184 Smin = -0.01;
185
186 if (volVars.saturation(FluidSystem::nPhaseIdx) <= Smin)
187 {
188 wouldSwitch = true;
189 // NAPL phase disappears
190 if (this->verbosity() > 1)
191 std::cout << "NAPL phase disappears at dof " << dofIdxGlobal
192 << ", coordinates: " << globalPos << ", sn: "
193 << volVars.saturation(FluidSystem::nPhaseIdx) << std::endl;
194 nonwettingPresent = true;
195 }
196
197
198 // calculate fractions of the hypothetical water phase
199 Scalar xww = volVars.moleFraction(FluidSystem::wPhaseIdx, FluidSystem::wCompIdx);
200 /*
201 take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
202 If this is larger than 1, then water appears
203 */
204 Scalar xwMax = 1.0;
205 if (xww > xwMax)
206 wouldSwitch = true;
207 if (this->wasSwitched_[dofIdxGlobal])
208 xwMax *= 1.02;
209
210 // if the sum of the mole fractions would be larger than
211 // 100%, gas phase appears
212 if (xww > xwMax)
213 {
214 // water phase appears
215 if (this->verbosity() > 1)
216 std::cout << "water phase appears at dof " << dofIdxGlobal
217 << ", coordinates: " << globalPos << ", xww=xwg*pg/pwsat : "
218 << xww << std::endl;
219 wettingPresent = true;
220 }
221
222 if (wettingPresent && !nonwettingPresent)
223 {
224 newPhasePresence = Indices::threePhases;
225 priVars[Indices::switch1Idx] = 0.0001;
226 priVars[Indices::switch2Idx] = volVars.saturation(FluidSystem::nPhaseIdx);
227 }
228 else if (wettingPresent && nonwettingPresent)
229 {
230 newPhasePresence = Indices::wgPhaseOnly;
231 priVars[Indices::switch1Idx] = 0.0001;
232 priVars[Indices::switch2Idx] = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::nCompIdx);
233 }
234 else if (!wettingPresent && nonwettingPresent)
235 {
236 newPhasePresence = Indices::gPhaseOnly;
237 priVars[Indices::switch1Idx] = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::wCompIdx);
238 priVars[Indices::switch2Idx] = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::nCompIdx);
239 }
240 }
241 else if (phasePresence == Indices::wnPhaseOnly)
242 {
243 bool nonwettingPresent = false;
244 bool gasPresent = false;
245
246 Scalar Smin = 0.0;
247 if (this->wasSwitched_[dofIdxGlobal])
248 Smin = -0.01;
249
250 if (volVars.saturation(FluidSystem::nPhaseIdx) <= Smin)
251 {
252 wouldSwitch = true;
253 // NAPL phase disappears
254 if (this->verbosity() > 1)
255 std::cout << "NAPL phase disappears at dof " << dofIdxGlobal
256 << ", coordinates: " << globalPos << ", sn: "
257 << volVars.saturation(FluidSystem::nPhaseIdx) << std::endl;
258 nonwettingPresent = true;
259 }
260
261 // calculate fractions of the partial pressures in the
262 // hypothetical gas phase
263 Scalar xwg = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::wCompIdx);
264 Scalar xgg = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::gCompIdx);
265 Scalar xng = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::nCompIdx);
266 /* take care:
267 for xgg in case Indices::wPhaseOnly we compute xgg=henry_air*x2w
268 for xwg in case Indices::wPhaseOnly we compute xwg=pwsat
269 for xng in case Indices::wPhaseOnly we compute xng=henry_NAPL*x1w
270 */
271 Scalar xgMax = 1.0;
272 if (xwg + xgg + xng > xgMax)
273 wouldSwitch = true;
274 if (this->wasSwitched_[dofIdxGlobal])
275 xgMax *= 1.02;
276
277 // if the sum of the mole fractions would be larger than
278 // 100%, gas phase appears
279 if (xwg + xgg + xng > xgMax)
280 {
281 // gas phase appears
282 if (this->verbosity() > 1)
283 std::cout << "gas phase appears at dof " << dofIdxGlobal
284 << ", coordinates: " << globalPos << ", xwg + xgg + xng: "
285 << xwg + xgg + xng << std::endl;
286 gasPresent = true;
287 }
288
289 if (gasPresent && !nonwettingPresent)
290 {
291 newPhasePresence = Indices::threePhases;
292 priVars[Indices::switch1Idx] = volVars.saturation(FluidSystem::wPhaseIdx);
293 priVars[Indices::switch2Idx] = volVars.saturation(FluidSystem::nPhaseIdx);
294 }
295 else if (gasPresent && nonwettingPresent)
296 {
297 newPhasePresence = Indices::wgPhaseOnly;
298 priVars[Indices::switch1Idx] = volVars.saturation(FluidSystem::wPhaseIdx);
299 priVars[Indices::switch2Idx] = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::nCompIdx);
300 }
301 else if (!gasPresent && nonwettingPresent)
302 {
303 newPhasePresence = Indices::wPhaseOnly;
304 priVars[Indices::switch1Idx] = volVars.moleFraction(FluidSystem::wPhaseIdx, FluidSystem::gCompIdx);
305 priVars[Indices::switch2Idx] = volVars.moleFraction(FluidSystem::wPhaseIdx, FluidSystem::nCompIdx);
306 }
307 }
308 else if (phasePresence == Indices::gPhaseOnly)
309 {
310 bool nonwettingPresent = false;
311 bool wettingPresent = false;
312
313 // calculate fractions in the hypothetical NAPL phase
314 Scalar xnn = volVars.moleFraction(FluidSystem::nPhaseIdx, FluidSystem::nCompIdx);
315 /*
316 take care:, xnn, if no NAPL phase is there, take xnn=xng*pg/pcsat
317 if this is larger than 1, then NAPL appears
318 */
319
320 Scalar xnMax = 1.0;
321 if (xnn > xnMax)
322 wouldSwitch = true;
323 if (this->wasSwitched_[dofIdxGlobal])
324 xnMax *= 1.02;
325
326 // if the sum of the hypothetical mole fraction would be larger than
327 // 100%, NAPL phase appears
328 if (xnn > xnMax)
329 {
330 // NAPL phase appears
331 if (this->verbosity() > 1)
332 std::cout << "NAPL phase appears at dof " << dofIdxGlobal
333 << ", coordinates: " << globalPos << ", xnn: "
334 << xnn << std::endl;
335 nonwettingPresent = true;
336 }
337 // calculate fractions of the hypothetical water phase
338 Scalar xww = volVars.moleFraction(FluidSystem::wPhaseIdx, FluidSystem::wCompIdx);
339 /*
340 take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
341 If this is larger than 1, then water appears
342 */
343 Scalar xwMax = 1.0;
344 if (xww > xwMax)
345 wouldSwitch = true;
346 if (this->wasSwitched_[dofIdxGlobal])
347 xwMax *= 1.02;
348
349 // if the sum of the mole fractions would be larger than
350 // 100%, gas phase appears
351 if (xww > xwMax)
352 {
353 // water phase appears
354 if (this->verbosity() > 1)
355 std::cout << "water phase appears at dof " << dofIdxGlobal
356 << ", coordinates: " << globalPos << ", xww=xwg*pg/pwsat : "
357 << xww << std::endl;
358 wettingPresent = true;
359 }
360 if (wettingPresent && !nonwettingPresent)
361 {
362 newPhasePresence = Indices::wgPhaseOnly;
363 priVars[Indices::switch1Idx] = 0.0001;
364 priVars[Indices::switch2Idx] = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::nCompIdx);
365 }
366 else if (wettingPresent && nonwettingPresent)
367 {
368 newPhasePresence = Indices::threePhases;
369 priVars[Indices::switch1Idx] = 0.0001;
370 priVars[Indices::switch2Idx] = 0.0001;
371 }
372 else if (!wettingPresent && nonwettingPresent)
373 {
374 newPhasePresence = Indices::gnPhaseOnly;
375 priVars[Indices::switch1Idx] = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::wCompIdx);
376 priVars[Indices::switch2Idx] = 0.0001;
377 }
378 }
379 else if (phasePresence == Indices::wgPhaseOnly)
380 {
381 bool nonwettingPresent = false;
382 bool gasPresent = false;
383 bool wettingPresent = false;
384
385 // get the fractions in the hypothetical NAPL phase
386 Scalar xnn = volVars.moleFraction(FluidSystem::nPhaseIdx, FluidSystem::nCompIdx);
387
388 // take care: if the NAPL phase is not present, take
389 // xnn=xng*pg/pcsat if this is larger than 1, then NAPL
390 // appears
391 Scalar xnMax = 1.0;
392 if (xnn > xnMax)
393 wouldSwitch = true;
394 if (this->wasSwitched_[dofIdxGlobal])
395 xnMax *= 1.02;
396
397 // if the sum of the hypothetical mole fraction would be larger than
398 // 100%, NAPL phase appears
399 if (xnn > xnMax)
400 {
401 // NAPL phase appears
402 if (this->verbosity() > 1)
403 std::cout << "NAPL phase appears at dof " << dofIdxGlobal
404 << ", coordinates: " << globalPos << ", xnn: "
405 << xnn << std::endl;
406 nonwettingPresent = true;
407 }
408
409 Scalar Smin = -1.e-6;
410 if (this->wasSwitched_[dofIdxGlobal])
411 Smin = -0.01;
412
413 if (volVars.saturation(FluidSystem::gPhaseIdx) <= Smin)
414 {
415 wouldSwitch = true;
416 // gas phase disappears
417 if (this->verbosity() > 1)
418 std::cout << "Gas phase disappears at dof " << dofIdxGlobal
419 << ", coordinates: " << globalPos << ", sg: "
420 << volVars.saturation(FluidSystem::gPhaseIdx) << std::endl;
421 gasPresent = true;
422 }
423
424 Smin = 0.0;
425 if (this->wasSwitched_[dofIdxGlobal])
426 Smin = -0.01;
427
428 if (volVars.saturation(FluidSystem::wPhaseIdx) <= Smin)
429 {
430 wouldSwitch = true;
431 // gas phase disappears
432 if (this->verbosity() > 1)
433 std::cout << "Water phase disappears at dof " << dofIdxGlobal
434 << ", coordinates: " << globalPos << ", sw: "
435 << volVars.saturation(FluidSystem::wPhaseIdx) << std::endl;
436 wettingPresent = true;
437 }
438
439 if (!gasPresent && nonwettingPresent && wettingPresent)
440 {
441 newPhasePresence = Indices::gnPhaseOnly;
442 priVars[Indices::switch1Idx] = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::wCompIdx);
443 priVars[Indices::switch2Idx] = 0.0001;
444 }
445 else if (!gasPresent && nonwettingPresent && !wettingPresent)
446 {
447 newPhasePresence = Indices::threePhases;
448 priVars[Indices::switch1Idx] = volVars.saturation(FluidSystem::wPhaseIdx);
449 priVars[Indices::switch2Idx] = 0.0;
450 }
451 else if (gasPresent && !nonwettingPresent && !wettingPresent)
452 {
453 newPhasePresence = Indices::wPhaseOnly;
454 priVars[Indices::switch1Idx] = volVars.moleFraction(FluidSystem::wPhaseIdx, FluidSystem::gCompIdx);
455 priVars[Indices::switch2Idx] = volVars.moleFraction(FluidSystem::wPhaseIdx, FluidSystem::nCompIdx);
456 }
457 else if (!gasPresent && !nonwettingPresent && wettingPresent)
458 {
459 newPhasePresence = Indices::gPhaseOnly;
460 priVars[Indices::switch1Idx] = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::wCompIdx);
461 priVars[Indices::switch2Idx] = volVars.moleFraction(FluidSystem::gPhaseIdx, FluidSystem::nCompIdx);
462 }
463 }
464
465 priVars.setState(newPhasePresence);
466 this->wasSwitched_[dofIdxGlobal] = wouldSwitch;
467 return phasePresence != newPhasePresence;
468 }
469};
470
471} // end namespace dumux
472
473#endif
Dumux::PrimaryVariableSwitch
The primary variable switch controlling the phase presence state variable.
Definition: compositional/primaryvariableswitch.hh:49
Dumux::PrimaryVariableSwitch< ThreePThreeCPrimaryVariableSwitch >::wasSwitched_
std::vector< bool > wasSwitched_
Definition: compositional/primaryvariableswitch.hh:433
Dumux::PrimaryVariableSwitch< ThreePThreeCPrimaryVariableSwitch >::verbosity
int verbosity() const
The verbosity level.
Definition: compositional/primaryvariableswitch.hh:271
Dumux::ThreePThreeCPrimaryVariableSwitch
The primary variable switch controlling the phase presence state variable.
Definition: 3p3c/primaryvariableswitch.hh:28
Dumux::ThreePThreeCPrimaryVariableSwitch::update_
bool update_(typename VolumeVariables::PrimaryVariables &priVars, const VolumeVariables &volVars, std::size_t dofIdxGlobal, const GlobalPosition &globalPos)
Definition: 3p3c/primaryvariableswitch.hh:39
primaryvariableswitch.hh
The primary variable switch base class for compositional models.
Dumux::IOName::phasePresence
std::string phasePresence() noexcept
I/O name of phase presence.
Definition: name.hh:135
Dumux
Definition: adapt.hh:17