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