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