3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
2pnc/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_2PNC_PRIMARY_VARIABLE_SWITCH_HH
26#define DUMUX_2PNC_PRIMARY_VARIABLE_SWITCH_HH
27
28#include <iostream>
29
32
33namespace Dumux {
34
40: public PrimaryVariableSwitch<TwoPNCPrimaryVariableSwitch>
41{
43 friend ParentType;
44public:
45 using ParentType::ParentType;
46
47protected:
48 // perform variable switch at a degree of freedom location
49 template<class VolumeVariables, class IndexType, class GlobalPosition>
50 bool update_(typename VolumeVariables::PrimaryVariables& priVars,
51 const VolumeVariables& volVars,
52 IndexType dofIdxGlobal,
53 const GlobalPosition& globalPos)
54 {
55 using Scalar = typename VolumeVariables::PrimaryVariables::value_type;
56
57 using FluidSystem = typename VolumeVariables::FluidSystem;
58 static constexpr int phase0Idx = FluidSystem::phase0Idx;
59 static constexpr int phase1Idx = FluidSystem::phase1Idx;
60 static constexpr int comp0Idx = FluidSystem::comp0Idx;
61 static constexpr int comp1Idx = FluidSystem::comp1Idx;
62
63 static constexpr auto numComponents = VolumeVariables::numFluidComponents();
64 static constexpr bool useMoles = VolumeVariables::useMoles();
65 static_assert(useMoles || numComponents < 3, "!useMoles is only implemented for numComponents < 3.");
66 static constexpr auto numMajorComponents = VolumeVariables::numFluidPhases();
67 static constexpr auto formulation = VolumeVariables::priVarFormulation();
68 static_assert( (formulation == TwoPFormulation::p0s1 || formulation == TwoPFormulation::p1s0),
69 "Chosen TwoPFormulation not supported!");
70
71 using Indices = typename VolumeVariables::Indices;
72 static constexpr int switchIdx = Indices::switchIdx;
73
74 // evaluate primary variable switch
75 bool wouldSwitch = false;
76 int phasePresence = priVars.state();
77 int newPhasePresence = phasePresence;
78
79 //check if a primary variable switch is necessary
80 if (phasePresence == Indices::bothPhases)
81 {
82 Scalar Smin = 0; //saturation threshold
83 if (this->wasSwitched_[dofIdxGlobal])
84 Smin = -0.01;
85
86 // if saturation of first phase is smaller 0: switch
87 if (volVars.saturation(phase0Idx) <= Smin)
88 {
89 wouldSwitch = true;
90 // first phase has to disappear
91 if (this->verbosity() > 1)
92 std::cout << "First phase (" << FluidSystem::phaseName(phase0Idx) << ")"
93 << " disappears at dof " << dofIdxGlobal
94 << ", coordinates: " << globalPos
95 << ", S_" << FluidSystem::phaseName(phase0Idx) << ": " << volVars.saturation(phase0Idx)
96 << std::endl;
97 newPhasePresence = Indices::secondPhaseOnly;
98
99 // switch not depending on formulation, switch "S0" to "x10"
100 if(useMoles) // mole-fraction formulation
101 priVars[switchIdx] = volVars.moleFraction(phase1Idx, comp0Idx);
102 else // mass-fraction formulation
103 priVars[switchIdx] = volVars.massFraction(phase1Idx, comp0Idx);
104
105 // switch all secondary components to mole fraction in nonwetting phase
106 if(useMoles) // mole-fraction formulation
107 for (int compIdx = numMajorComponents; compIdx < numComponents; ++compIdx)
108 priVars[compIdx] = volVars.moleFraction(phase1Idx, compIdx);
109 else // mass-fraction formulation
110 for (int compIdx = numMajorComponents; compIdx < numComponents; ++compIdx)
111 priVars[compIdx] = volVars.massFraction(phase1Idx, compIdx);
112 }
113
114 // if saturation of second phase is smaller than 0: switch
115 else if (volVars.saturation(phase1Idx) <= Smin)
116 {
117 wouldSwitch = true;
118 // second phase has to disappear
119 if (this->verbosity() > 1)
120 std::cout << "Second phase (" << FluidSystem::phaseName(phase1Idx) << ")"
121 << " disappears at dof " << dofIdxGlobal
122 << ", coordinates: " << globalPos
123 << ", S_" << FluidSystem::phaseName(phase1Idx) << ": " << volVars.saturation(phase1Idx)
124 << std::endl;
125 newPhasePresence = Indices::firstPhaseOnly;
126
127 // switch "S1" to "x01"
128 if(useMoles) // mole-fraction formulation
129 priVars[switchIdx] = volVars.moleFraction(phase0Idx, comp1Idx);
130 else // mass-fraction formulation
131 priVars[switchIdx] = volVars.massFraction(phase0Idx, comp1Idx);
132 }
133 }
134 else if (phasePresence == Indices::secondPhaseOnly)
135 {
136 Scalar x0Max = 1;
137 Scalar x0Sum = 0;
138 // Calculate sum of mole fractions in the hypothetical first phase
139 for (int compIdx = 0; compIdx < numComponents; compIdx++)
140 x0Sum += volVars.moleFraction(phase0Idx, compIdx);
141
142 if (x0Sum > x0Max)
143 wouldSwitch = true;
144 if (this->wasSwitched_[dofIdxGlobal])
145 x0Max *= 1.02;
146
147 // first phase appears if sum is larger than one
148 if (x0Sum/*sum of mole fractions*/ > x0Max/*1*/)
149 {
150 if (this->verbosity() > 1)
151 std::cout << "Second phase (" << FluidSystem::phaseName(phase0Idx) << ")"
152 << " appears at dof " << dofIdxGlobal
153 << ", coordinates: " << globalPos
154 << ", sum x^i_" << FluidSystem::phaseName(phase0Idx) << ": " << x0Sum
155 << std::endl;
156 newPhasePresence = Indices::bothPhases;
157
158 // saturation of the first phase set to 0.0001 (if formulation TwoPFormulation::p1s0 and vice versa)
159 if (formulation == TwoPFormulation::p1s0)
160 priVars[switchIdx] = 0.0001;
161 else
162 priVars[switchIdx] = 0.9999;
163
164 // switch all secondary components back to first component mole fraction
165 for (int compIdx = numMajorComponents; compIdx < numComponents; ++compIdx)
166 priVars[compIdx] = volVars.moleFraction(phase0Idx,compIdx);
167 }
168 }
169 else if (phasePresence == Indices::firstPhaseOnly)
170 {
171 Scalar x1Max = 1;
172 Scalar x1Sum = 0;
173
174 // Calculate sum of mole fractions in the hypothetical wetting phase
175 for (int compIdx = 0; compIdx < numComponents; compIdx++)
176 x1Sum += volVars.moleFraction(phase1Idx, compIdx);
177
178 if (x1Sum > x1Max)
179 wouldSwitch = true;
180 if (this->wasSwitched_[dofIdxGlobal])
181 x1Max *= 1.02;
182
183 // wetting phase appears if sum is larger than one
184 if (x1Sum > x1Max)
185 {
186 if (this->verbosity() > 1)
187 std::cout << "Second phase (" << FluidSystem::phaseName(phase1Idx) << ")"
188 << " appears at dof " << dofIdxGlobal
189 << ", coordinates: " << globalPos
190 << ", sum x^i_" << FluidSystem::phaseName(phase1Idx) << ": " << x1Sum
191 << std::endl;
192 newPhasePresence = Indices::bothPhases;
193 //saturation of the wetting phase set to 0.9999 (if formulation TwoPFormulation::pnsw and vice versa)
194 if (formulation == TwoPFormulation::p1s0)
195 priVars[switchIdx] = 0.9999;
196 else
197 priVars[switchIdx] = 0.0001;
198 }
199 }
200
201 priVars.setState(newPhasePresence);
202 this->wasSwitched_[dofIdxGlobal] = wouldSwitch;
203 return phasePresence != newPhasePresence;
204 }
205};
206
207} // end namespace Dumux
208
209#endif
Defines an enumeration for the formulations accepted by the two-phase model.
@ p1s0
first phase saturation and second phase pressure as primary variables
@ p0s1
first phase pressure and second phase saturation as primary variables
Definition: adapt.hh:29
std::string phasePresence() noexcept
I/O name of phase presence.
Definition: name.hh:147
The primary variable switch controlling the phase presence state variable.
Definition: 2pnc/primaryvariableswitch.hh:41
bool update_(typename VolumeVariables::PrimaryVariables &priVars, const VolumeVariables &volVars, IndexType dofIdxGlobal, const GlobalPosition &globalPos)
Definition: 2pnc/primaryvariableswitch.hh:50
The primary variable switch controlling the phase presence state variable.
Definition: compositional/primaryvariableswitch.hh:60
std::vector< bool > wasSwitched_
Definition: compositional/primaryvariableswitch.hh:463
int verbosity() const
The verbosity level.
Definition: compositional/primaryvariableswitch.hh:301
The primary variable switch base class for compositional models.