| Line | Branch | Exec | Source |
|---|---|---|---|
| 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-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup TwoPNCModel | ||
| 10 | * \brief The primary variable switch for the 2pnc model. | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_2PNC_PRIMARY_VARIABLE_SWITCH_HH | ||
| 14 | #define DUMUX_2PNC_PRIMARY_VARIABLE_SWITCH_HH | ||
| 15 | |||
| 16 | #include <iostream> | ||
| 17 | |||
| 18 | #include <dumux/porousmediumflow/compositional/primaryvariableswitch.hh> | ||
| 19 | #include <dumux/porousmediumflow/2p/formulation.hh> | ||
| 20 | |||
| 21 | namespace Dumux { | ||
| 22 | |||
| 23 | /*! | ||
| 24 | * \ingroup TwoPNCModel | ||
| 25 | * \brief The primary variable switch controlling the phase presence state variable. | ||
| 26 | */ | ||
| 27 | 68 | class TwoPNCPrimaryVariableSwitch | |
| 28 | : public PrimaryVariableSwitch<TwoPNCPrimaryVariableSwitch> | ||
| 29 | { | ||
| 30 | using ParentType = PrimaryVariableSwitch<TwoPNCPrimaryVariableSwitch>; | ||
| 31 | friend ParentType; | ||
| 32 | public: | ||
| 33 | 34 | using ParentType::ParentType; | |
| 34 | |||
| 35 | protected: | ||
| 36 | // perform variable switch at a degree of freedom location | ||
| 37 | template<class VolumeVariables, class IndexType, class GlobalPosition> | ||
| 38 | 1912281 | bool update_(typename VolumeVariables::PrimaryVariables& priVars, | |
| 39 | const VolumeVariables& volVars, | ||
| 40 | IndexType dofIdxGlobal, | ||
| 41 | const GlobalPosition& globalPos) | ||
| 42 | { | ||
| 43 | using Scalar = typename VolumeVariables::PrimaryVariables::value_type; | ||
| 44 | |||
| 45 | using FluidSystem = typename VolumeVariables::FluidSystem; | ||
| 46 | static constexpr int phase0Idx = FluidSystem::phase0Idx; | ||
| 47 | static constexpr int phase1Idx = FluidSystem::phase1Idx; | ||
| 48 | static constexpr int comp0Idx = FluidSystem::comp0Idx; | ||
| 49 | static constexpr int comp1Idx = FluidSystem::comp1Idx; | ||
| 50 | |||
| 51 | static constexpr auto numComponents = VolumeVariables::numFluidComponents(); | ||
| 52 | static constexpr bool useMoles = VolumeVariables::useMoles(); | ||
| 53 | static_assert(useMoles || numComponents < 3, "!useMoles is only implemented for numComponents < 3."); | ||
| 54 | static constexpr auto numMajorComponents = VolumeVariables::numFluidPhases(); | ||
| 55 | static constexpr auto formulation = VolumeVariables::priVarFormulation(); | ||
| 56 | static_assert( (formulation == TwoPFormulation::p0s1 || formulation == TwoPFormulation::p1s0), | ||
| 57 | "Chosen TwoPFormulation not supported!"); | ||
| 58 | |||
| 59 | using Indices = typename VolumeVariables::Indices; | ||
| 60 | static constexpr int switchIdx = Indices::switchIdx; | ||
| 61 | |||
| 62 | // evaluate primary variable switch | ||
| 63 | 1912281 | bool wouldSwitch = false; | |
| 64 | 1912281 | int phasePresence = priVars.state(); | |
| 65 | 1912281 | int newPhasePresence = phasePresence; | |
| 66 | |||
| 67 | //check if a primary variable switch is necessary | ||
| 68 |
2/2✓ Branch 0 taken 825807 times.
✓ Branch 1 taken 1086474 times.
|
1912281 | if (phasePresence == Indices::bothPhases) |
| 69 | { | ||
| 70 | 825807 | Scalar Smin = 0; //saturation threshold | |
| 71 |
2/2✓ Branch 0 taken 377 times.
✓ Branch 1 taken 825430 times.
|
825807 | if (this->wasSwitched_[dofIdxGlobal]) |
| 72 | 377 | Smin = -0.01; | |
| 73 | |||
| 74 | // if saturation of first phase is smaller 0: switch | ||
| 75 |
2/2✓ Branch 0 taken 95 times.
✓ Branch 1 taken 825712 times.
|
825807 | if (volVars.saturation(phase0Idx) <= Smin) |
| 76 | { | ||
| 77 | 95 | wouldSwitch = true; | |
| 78 | // first phase has to disappear | ||
| 79 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 95 times.
|
95 | if (this->verbosity() > 1) |
| 80 | std::cout << "First phase (" << FluidSystem::phaseName(phase0Idx) << ")" | ||
| 81 | ✗ | << " disappears at dof " << dofIdxGlobal | |
| 82 | ✗ | << ", coordinates: " << globalPos | |
| 83 | ✗ | << ", S_" << FluidSystem::phaseName(phase0Idx) << ": " << volVars.saturation(phase0Idx) | |
| 84 | ✗ | << std::endl; | |
| 85 | 95 | newPhasePresence = Indices::secondPhaseOnly; | |
| 86 | |||
| 87 | // switch not depending on formulation, switch "S0" to "x10" | ||
| 88 | if(useMoles) // mole-fraction formulation | ||
| 89 | 95 | priVars[switchIdx] = volVars.moleFraction(phase1Idx, comp0Idx); | |
| 90 | else // mass-fraction formulation | ||
| 91 | priVars[switchIdx] = volVars.massFraction(phase1Idx, comp0Idx); | ||
| 92 | |||
| 93 | // switch all secondary components to mole fraction in nonwetting phase | ||
| 94 | if(useMoles) // mole-fraction formulation | ||
| 95 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
95 | for (int compIdx = numMajorComponents; compIdx < numComponents; ++compIdx) |
| 96 | ✗ | priVars[compIdx] = volVars.moleFraction(phase1Idx, compIdx); | |
| 97 | else // mass-fraction formulation | ||
| 98 | for (int compIdx = numMajorComponents; compIdx < numComponents; ++compIdx) | ||
| 99 | priVars[compIdx] = volVars.massFraction(phase1Idx, compIdx); | ||
| 100 | } | ||
| 101 | |||
| 102 | // if saturation of second phase is smaller than 0: switch | ||
| 103 |
2/2✓ Branch 0 taken 315 times.
✓ Branch 1 taken 825397 times.
|
825712 | else if (volVars.saturation(phase1Idx) <= Smin) |
| 104 | { | ||
| 105 | 315 | wouldSwitch = true; | |
| 106 | // second phase has to disappear | ||
| 107 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 315 times.
|
315 | if (this->verbosity() > 1) |
| 108 | std::cout << "Second phase (" << FluidSystem::phaseName(phase1Idx) << ")" | ||
| 109 | ✗ | << " disappears at dof " << dofIdxGlobal | |
| 110 | ✗ | << ", coordinates: " << globalPos | |
| 111 | ✗ | << ", S_" << FluidSystem::phaseName(phase1Idx) << ": " << volVars.saturation(phase1Idx) | |
| 112 | ✗ | << std::endl; | |
| 113 | 315 | newPhasePresence = Indices::firstPhaseOnly; | |
| 114 | |||
| 115 | // switch "S1" to "x01" | ||
| 116 | if(useMoles) // mole-fraction formulation | ||
| 117 | 315 | priVars[switchIdx] = volVars.moleFraction(phase0Idx, comp1Idx); | |
| 118 | else // mass-fraction formulation | ||
| 119 | priVars[switchIdx] = volVars.massFraction(phase0Idx, comp1Idx); | ||
| 120 | } | ||
| 121 | } | ||
| 122 |
2/2✓ Branch 0 taken 1047 times.
✓ Branch 1 taken 1085427 times.
|
1086474 | else if (phasePresence == Indices::secondPhaseOnly) |
| 123 | { | ||
| 124 | 3141 | Scalar x0Max = 1; | |
| 125 | Scalar x0Sum = 0; | ||
| 126 | // Calculate sum of mole fractions in the hypothetical first phase | ||
| 127 |
2/2✓ Branch 0 taken 2094 times.
✓ Branch 1 taken 1047 times.
|
3141 | for (int compIdx = 0; compIdx < numComponents; compIdx++) |
| 128 | 2094 | x0Sum += volVars.moleFraction(phase0Idx, compIdx); | |
| 129 | |||
| 130 |
2/2✓ Branch 0 taken 40 times.
✓ Branch 1 taken 1007 times.
|
1047 | if (x0Sum > x0Max) |
| 131 | 40 | wouldSwitch = true; | |
| 132 |
2/2✓ Branch 0 taken 68 times.
✓ Branch 1 taken 979 times.
|
1047 | if (this->wasSwitched_[dofIdxGlobal]) |
| 133 | 68 | x0Max *= 1.02; | |
| 134 | |||
| 135 | // first phase appears if sum is larger than one | ||
| 136 |
2/2✓ Branch 0 taken 40 times.
✓ Branch 1 taken 1007 times.
|
1047 | if (x0Sum/*sum of mole fractions*/ > x0Max/*1*/) |
| 137 | { | ||
| 138 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
|
40 | if (this->verbosity() > 1) |
| 139 | std::cout << "Second phase (" << FluidSystem::phaseName(phase0Idx) << ")" | ||
| 140 | ✗ | << " appears at dof " << dofIdxGlobal | |
| 141 | ✗ | << ", coordinates: " << globalPos | |
| 142 | ✗ | << ", sum x^i_" << FluidSystem::phaseName(phase0Idx) << ": " << x0Sum | |
| 143 | ✗ | << std::endl; | |
| 144 | 40 | newPhasePresence = Indices::bothPhases; | |
| 145 | |||
| 146 | // saturation of the first phase set to 0.0001 (if formulation TwoPFormulation::p1s0 and vice versa) | ||
| 147 | if (formulation == TwoPFormulation::p1s0) | ||
| 148 | ✗ | priVars[switchIdx] = 0.0001; | |
| 149 | else | ||
| 150 | 40 | priVars[switchIdx] = 0.9999; | |
| 151 | |||
| 152 | // switch all secondary components back to first component mole fraction | ||
| 153 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
40 | for (int compIdx = numMajorComponents; compIdx < numComponents; ++compIdx) |
| 154 | ✗ | priVars[compIdx] = volVars.moleFraction(phase0Idx,compIdx); | |
| 155 | } | ||
| 156 | } | ||
| 157 |
1/2✓ Branch 0 taken 1085427 times.
✗ Branch 1 not taken.
|
1085427 | else if (phasePresence == Indices::firstPhaseOnly) |
| 158 | { | ||
| 159 | 3899717 | Scalar x1Max = 1; | |
| 160 | Scalar x1Sum = 0; | ||
| 161 | |||
| 162 | // Calculate sum of mole fractions in the hypothetical wetting phase | ||
| 163 |
2/2✓ Branch 0 taken 2814290 times.
✓ Branch 1 taken 1085427 times.
|
3899717 | for (int compIdx = 0; compIdx < numComponents; compIdx++) |
| 164 | 2814290 | x1Sum += volVars.moleFraction(phase1Idx, compIdx); | |
| 165 | |||
| 166 |
2/2✓ Branch 0 taken 347 times.
✓ Branch 1 taken 1085080 times.
|
1085427 | if (x1Sum > x1Max) |
| 167 | 347 | wouldSwitch = true; | |
| 168 |
2/2✓ Branch 0 taken 315 times.
✓ Branch 1 taken 1085112 times.
|
1085427 | if (this->wasSwitched_[dofIdxGlobal]) |
| 169 | 315 | x1Max *= 1.02; | |
| 170 | |||
| 171 | // wetting phase appears if sum is larger than one | ||
| 172 |
2/2✓ Branch 0 taken 344 times.
✓ Branch 1 taken 1085083 times.
|
1085427 | if (x1Sum > x1Max) |
| 173 | { | ||
| 174 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 344 times.
|
344 | if (this->verbosity() > 1) |
| 175 | std::cout << "Second phase (" << FluidSystem::phaseName(phase1Idx) << ")" | ||
| 176 | ✗ | << " appears at dof " << dofIdxGlobal | |
| 177 | ✗ | << ", coordinates: " << globalPos | |
| 178 | ✗ | << ", sum x^i_" << FluidSystem::phaseName(phase1Idx) << ": " << x1Sum | |
| 179 | ✗ | << std::endl; | |
| 180 | 344 | newPhasePresence = Indices::bothPhases; | |
| 181 | //saturation of the wetting phase set to 0.9999 (if formulation TwoPFormulation::pnsw and vice versa) | ||
| 182 | if (formulation == TwoPFormulation::p1s0) | ||
| 183 | 41 | priVars[switchIdx] = 0.9999; | |
| 184 | else | ||
| 185 | 303 | priVars[switchIdx] = 0.0001; | |
| 186 | } | ||
| 187 | } | ||
| 188 | |||
| 189 | 1912281 | priVars.setState(newPhasePresence); | |
| 190 |
2/2✓ Branch 0 taken 797 times.
✓ Branch 1 taken 1911484 times.
|
1912281 | this->wasSwitched_[dofIdxGlobal] = wouldSwitch; |
| 191 | 1912281 | return phasePresence != newPhasePresence; | |
| 192 | } | ||
| 193 | }; | ||
| 194 | |||
| 195 | } // end namespace Dumux | ||
| 196 | |||
| 197 | #endif | ||
| 198 |