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