| 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-FileCopyrightInfo: 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 ThreePWaterOilModel | ||
| 10 | * \brief The primary variable switch for the 3p3c model. | ||
| 11 | */ | ||
| 12 | |||
| 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 | |||
| 19 | namespace Dumux { | ||
| 20 | |||
| 21 | /*! | ||
| 22 | * \ingroup ThreePWaterOilModel | ||
| 23 | * \brief The primary variable switch controlling the phase presence state variable. | ||
| 24 | */ | ||
| 25 | ✗ | class ThreePWaterOilPrimaryVariableSwitch | |
| 26 | : public PrimaryVariableSwitch<ThreePWaterOilPrimaryVariableSwitch> | ||
| 27 | { | ||
| 28 | using ParentType = PrimaryVariableSwitch<ThreePWaterOilPrimaryVariableSwitch>; | ||
| 29 | friend ParentType; | ||
| 30 | |||
| 31 | public: | ||
| 32 | ✗ | using ParentType::ParentType; | |
| 33 | |||
| 34 | protected: | ||
| 35 | |||
| 36 | // perform variable switch at a degree of freedom location | ||
| 37 | template<class VolumeVariables, class GlobalPosition> | ||
| 38 | 153726 | 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 | 153726 | bool wouldSwitch = false; | |
| 50 | 153726 | auto phasePresence = priVars.state(); | |
| 51 | 153726 | int newPhasePresence = phasePresence; | |
| 52 | |||
| 53 | if(VolumeVariables::onlyGasPhaseCanDisappear()) | ||
| 54 | { | ||
| 55 | // check if a primary var switch is necessary | ||
| 56 |
2/2✓ Branch 0 taken 28 times.
✓ Branch 1 taken 153698 times.
|
153726 | if (phasePresence == Indices::threePhases) |
| 57 | { | ||
| 58 | 28 | Scalar Smin = 0; | |
| 59 |
6/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 24 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 24 times.
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 24 times.
|
84 | if (this->wasSwitched_[dofIdxGlobal]) |
| 60 | 4 | Smin = -0.01; | |
| 61 | |||
| 62 |
4/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 26 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 26 times.
|
56 | if (volVars.saturation(FluidSystem::gPhaseIdx) <= Smin) |
| 63 | { | ||
| 64 | 2 | wouldSwitch = true; | |
| 65 | // gas phase disappears | ||
| 66 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | 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 | 2 | newPhasePresence = Indices::wnPhaseOnly; | |
| 71 | |||
| 72 | 6 | priVars[Indices::pressureIdx] = volVars.fluidState().pressure(FluidSystem::wPhaseIdx); | |
| 73 | 8 | priVars[Indices::switch1Idx] = volVars.fluidState().temperature(); | |
| 74 | } | ||
| 75 | } | ||
| 76 |
1/2✓ Branch 0 taken 153698 times.
✗ Branch 1 not taken.
|
153698 | else if (phasePresence == Indices::wnPhaseOnly) |
| 77 | { | ||
| 78 | 153698 | bool gasFlag = 0; | |
| 79 | |||
| 80 | // calculate fractions of the partial pressures in the | ||
| 81 | // hypothetical gas phase | ||
| 82 |
4/4✓ Branch 0 taken 6 times.
✓ Branch 1 taken 153692 times.
✓ Branch 2 taken 6 times.
✓ Branch 3 taken 153692 times.
|
307396 | Scalar xwg = volVars.fluidState().moleFraction(FluidSystem::gPhaseIdx, FluidSystem::wCompIdx); |
| 83 |
4/4✓ Branch 0 taken 6 times.
✓ Branch 1 taken 153692 times.
✓ Branch 2 taken 6 times.
✓ Branch 3 taken 153692 times.
|
307396 | Scalar xng = volVars.fluidState().moleFraction(FluidSystem::gPhaseIdx, FluidSystem::nCompIdx); |
| 84 | |||
| 85 | 153698 | Scalar xgMax = 1.0; | |
| 86 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 153692 times.
|
153698 | if (xwg + xng > xgMax) |
| 87 | 6 | wouldSwitch = true; | |
| 88 |
6/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 153696 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 153696 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 153696 times.
|
461094 | if (this->wasSwitched_[dofIdxGlobal]) |
| 89 | 2 | xgMax *= 1.02; | |
| 90 | |||
| 91 | // if the sum of the mole fractions would be larger than | ||
| 92 | // 100%, gas phase appears | ||
| 93 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 153692 times.
|
153698 | if (xwg + xng > xgMax) |
| 94 | { | ||
| 95 | // gas phase appears | ||
| 96 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | if (this->verbosity() > 1) |
| 97 | ✗ | std::cout << "gas phase appears at dof " << dofIdxGlobal | |
| 98 | ✗ | << ", coordinates: " << globalPos << ", xwg + xng: " | |
| 99 | ✗ | << xwg + xng << std::endl; | |
| 100 | 6 | gasFlag = 1; | |
| 101 | } | ||
| 102 | |||
| 103 | if (gasFlag == 1) | ||
| 104 | { | ||
| 105 | 6 | newPhasePresence = Indices::threePhases; | |
| 106 | 12 | priVars[Indices::pressureIdx] = volVars.pressure(FluidSystem::gPhaseIdx); | |
| 107 | 18 | 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 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 153718 times.
|
153726 | priVars.setState(newPhasePresence); |
| 521 |
4/4✓ Branch 0 taken 8 times.
✓ Branch 1 taken 153718 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 153718 times.
|
307452 | this->wasSwitched_[dofIdxGlobal] = wouldSwitch; |
| 522 | 153726 | return phasePresence != newPhasePresence; | |
| 523 | } | ||
| 524 | |||
| 525 | }; | ||
| 526 | |||
| 527 | } // end namespace dumux | ||
| 528 | |||
| 529 | #endif | ||
| 530 |