diff --git a/dumux/io/plotthermalconductivitymodel.hh b/dumux/io/plotthermalconductivitymodel.hh index d7be67324457084db66a698e06463c9db3dcb1b1..c89c8b0e3739bc0d304cca584e83e836357ba056 100644 --- a/dumux/io/plotthermalconductivitymodel.hh +++ b/dumux/io/plotthermalconductivitymodel.hh @@ -122,9 +122,9 @@ private: using FluidSystem = typename PlotThermalConductivityModel::FluidSystem; - Scalar saturation(int phaseIdx) const + Scalar saturation(const int phaseIdx) const { - if (phaseIdx == wettingPhaseIdx()) + if (phaseIdx == wettingPhase()) return saturation_; else return 1.0 - saturation_; @@ -132,14 +132,14 @@ private: Scalar fluidThermalConductivity(const int phaseIdx) const { - if (phaseIdx == wettingPhaseIdx()) + if (phaseIdx == wettingPhase()) return lambdaW_; else return lambdaN_; } - Scalar wettingPhaseIdx() const - { return phase0Idx;} + int wettingPhase() const + { return phase0Idx; } Scalar porosity() const { return porosity_; } diff --git a/dumux/material/fluidmatrixinteractions/2p/thermalconductivityjohansen.hh b/dumux/material/fluidmatrixinteractions/2p/thermalconductivityjohansen.hh index 20f8f0fe8453a2058c7fe7623be5a4f0bf09cfb2..2bee93d2ff5e5d003b10628feaf21f7b3f574a6d 100644 --- a/dumux/material/fluidmatrixinteractions/2p/thermalconductivityjohansen.hh +++ b/dumux/material/fluidmatrixinteractions/2p/thermalconductivityjohansen.hh @@ -82,8 +82,8 @@ public: // TODO: there should be an assertion that the indices are correct and 0 is actually the wetting phase! const Scalar sw = volVars.saturation(volVars.wettingPhaseIdx()); - const Scalar lambdaW = volVars.fluidThermalConductivity(volVars.wettingPhaseIdx()); - const Scalar lambdaN = volVars.fluidThermalConductivity(1-volVars.wettingPhaseIdx()); + const Scalar lambdaW = volVars.fluidThermalConductivity(volVars.wettingPhase()); + const Scalar lambdaN = volVars.fluidThermalConductivity(1-volVars.wettingPhase()); const Scalar lambdaSolid = volVars.solidThermalConductivity(); const Scalar porosity = volVars.porosity(); const Scalar rhoSolid = volVars.solidDensity(); diff --git a/dumux/material/fluidstates/compositional.hh b/dumux/material/fluidstates/compositional.hh index 91bf96b5c7500465ead6f2817045b2deb2575ffa..5cb672441e0cf70eeed60ddb7155d7f62f3ebb95 100644 --- a/dumux/material/fluidstates/compositional.hh +++ b/dumux/material/fluidstates/compositional.hh @@ -71,7 +71,7 @@ public: * on thermodynamic equilibrium required) *****************************************************/ /*! - * \brief Returns the index of the wetting phase in the + * \brief Returns the index of the most wetting phase in the * fluid-solid configuration (for porous medium systems). */ int wettingPhase() const { return wPhaseIdx_; } @@ -441,7 +441,7 @@ public: { viscosity_[phaseIdx] = value; } /*! - * \brief Set the index of the wetting phase + * \brief Set the index of the most wetting phase */ void setWettingPhase(int phaseIdx) { wPhaseIdx_ = phaseIdx; } @@ -460,8 +460,6 @@ protected: std::array<Scalar, numPhases> viscosity_ = {}; std::array<Scalar, numPhases> temperature_ = {}; - // For porous medium flow models, here we ... the index of the wetting - // phase (needed for vapor pressure evaluation if kelvin equation is used) int wPhaseIdx_{0}; }; diff --git a/dumux/material/fluidstates/immiscible.hh b/dumux/material/fluidstates/immiscible.hh index c84cb2ea6d11c5feba0c814a7be30870d93fb165..451b054ef0a8a3d66326471feb3a6fe5115aa9be 100644 --- a/dumux/material/fluidstates/immiscible.hh +++ b/dumux/material/fluidstates/immiscible.hh @@ -66,6 +66,12 @@ public: * on thermodynamic equilibrium required) *****************************************************/ + /*! + * \brief Returns the index of the most wetting phase in the + * fluid-solid configuration (for porous medium systems). + */ + int wettingPhase() const { return wPhaseIdx_; } + /*! * \brief Returns the saturation \f$S_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[-]}\f$. * @@ -333,6 +339,12 @@ public: */ void setViscosity(int phaseIdx, Scalar value) { viscosity_[phaseIdx] = value; } + + /*! + * \brief Set the index of the most wetting phase + */ + void setWettingPhase(int phaseIdx) + { wPhaseIdx_ = phaseIdx; } protected: //! zero-initialize all data members with braces syntax Scalar pressure_[numPhases] = {}; @@ -342,6 +354,8 @@ protected: Scalar enthalpy_[numPhases] = {}; Scalar viscosity_[numPhases] = {}; Scalar temperature_[numPhases] = {}; + + int wPhaseIdx_{0}; }; } // end namespace Dumux diff --git a/dumux/material/fluidstates/isothermalimmiscible.hh b/dumux/material/fluidstates/isothermalimmiscible.hh index bfa2425b570026be8f1c62fc9155b4ea3075527e..88e7feaec891a765161928dbdb5bd3085004bfd4 100644 --- a/dumux/material/fluidstates/isothermalimmiscible.hh +++ b/dumux/material/fluidstates/isothermalimmiscible.hh @@ -64,6 +64,13 @@ public: /***************************************************** * Generic access to fluid properties *****************************************************/ + + /*! + * \brief Returns the index of the most wetting phase in the + * fluid-solid configuration (for porous medium systems). + */ + int wettingPhase() const { return wPhaseIdx_; } + /*! * \brief Returns the saturation \f$S_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[-]}\f$. * @@ -276,6 +283,11 @@ public: void setViscosity(int phaseIdx, Scalar value) { viscosity_[phaseIdx] = value; } + /*! + * \brief Set the index of the most wetting phase + */ + void setWettingPhase(int phaseIdx) + { wPhaseIdx_ = phaseIdx; } protected: Scalar pressure_[numPhases] = {}; Scalar saturation_[numPhases] = {}; @@ -283,6 +295,8 @@ protected: Scalar molarDensity_[numPhases] = {}; Scalar viscosity_[numPhases] = {}; Scalar temperature_ = 0.0; + + int wPhaseIdx_{0}; }; } // end namespace Dumux diff --git a/dumux/porousmediumflow/2p/volumevariables.hh b/dumux/porousmediumflow/2p/volumevariables.hh index 4508ef25c9f88f693a6c8da543f7868c711ffaba..1812fcf8141eb068c84b5f2cfe1c3a7a84f56d2f 100644 --- a/dumux/porousmediumflow/2p/volumevariables.hh +++ b/dumux/porousmediumflow/2p/volumevariables.hh @@ -94,14 +94,15 @@ public: using MaterialLaw = typename Problem::SpatialParams::MaterialLaw; const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol); - const int nPhaseIdx = 1 - wPhaseIdx_; + const int wPhaseIdx = fluidState_.wettingPhase(); + const int nPhaseIdx = 1 - wPhaseIdx; - mobility_[wPhaseIdx_] = - MaterialLaw::krw(materialParams, fluidState_.saturation(wPhaseIdx_)) - / fluidState_.viscosity(wPhaseIdx_); + mobility_[wPhaseIdx] = + MaterialLaw::krw(materialParams, fluidState_.saturation(wPhaseIdx)) + / fluidState_.viscosity(wPhaseIdx); mobility_[nPhaseIdx] = - MaterialLaw::krn(materialParams, fluidState_.saturation(wPhaseIdx_)) + MaterialLaw::krn(materialParams, fluidState_.saturation(wPhaseIdx)) / fluidState_.viscosity(nPhaseIdx); // porosity calculation over inert volumefraction @@ -137,15 +138,16 @@ public: const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol); const auto& priVars = elemSol[scv.localDofIndex()]; - wPhaseIdx_ = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol); + const auto wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol); + fluidState.setWettingPhase(wPhaseIdx); if (formulation == TwoPFormulation::p0s1) { fluidState.setPressure(phase0Idx, priVars[pressureIdx]); - if (wPhaseIdx_ == phase1Idx) + if (fluidState.wettingPhase() == phase1Idx) { fluidState.setSaturation(phase1Idx, priVars[saturationIdx]); fluidState.setSaturation(phase0Idx, 1 - priVars[saturationIdx]); - pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx_)); + pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx)); fluidState.setPressure(phase1Idx, priVars[pressureIdx] - pc_); } else @@ -154,27 +156,27 @@ public: scv, elemSol, priVars[saturationIdx]); fluidState.setSaturation(phase1Idx, Sn); fluidState.setSaturation(phase0Idx, 1 - Sn); - pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx_)); + pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx)); fluidState.setPressure(phase1Idx, priVars[pressureIdx] + pc_); } } else if (formulation == TwoPFormulation::p1s0) { fluidState.setPressure(phase1Idx, priVars[pressureIdx]); - if (wPhaseIdx_ == phase1Idx) + if (wPhaseIdx == phase1Idx) { const auto Sn = Traits::SaturationReconstruction::reconstructSn(problem.spatialParams(), element, scv, elemSol, priVars[saturationIdx]); fluidState.setSaturation(phase0Idx, Sn); fluidState.setSaturation(phase1Idx, 1 - Sn); - pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx_)); + pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx)); fluidState.setPressure(phase0Idx, priVars[pressureIdx] + pc_); } else { fluidState.setSaturation(phase0Idx, priVars[saturationIdx]); fluidState.setSaturation(phase1Idx, 1 - priVars[saturationIdx]); - pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx_)); + pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx)); fluidState.setPressure(phase0Idx, priVars[pressureIdx] - pc_); } } @@ -287,15 +289,14 @@ public: /*! * \brief Returns the wetting phase index */ - const int wettingPhaseIdx() const - { return wPhaseIdx_; } + int wettingPhase() const + { return fluidState_.wettingPhase(); } protected: FluidState fluidState_; SolidState solidState_; private: - int wPhaseIdx_; Scalar pc_; Scalar porosity_; PermeabilityType permeability_; diff --git a/dumux/porousmediumflow/2p1c/volumevariables.hh b/dumux/porousmediumflow/2p1c/volumevariables.hh index f72c02c061813f5a790ee37158be15c79d069555..9438ad3594d9837a0898eae953bb9eba2d39d6a8 100644 --- a/dumux/porousmediumflow/2p1c/volumevariables.hh +++ b/dumux/porousmediumflow/2p1c/volumevariables.hh @@ -135,15 +135,16 @@ public: // became part of the fluid state. typename FluidSystem::ParameterCache paramCache; paramCache.updateAll(fluidState_); + const int wPhaseIdx = fluidState_.wettingPhase(); for (int phaseIdx = 0; phaseIdx < numFluidPhases; ++phaseIdx) { // relative permeabilities Scalar kr; - if (phaseIdx == wPhaseIdx_) - kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx_)); + if (phaseIdx == wPhaseIdx) + kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx)); else // ATTENTION: krn requires the wetting phase saturation // as parameter! - kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx_)); + kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx)); relativePermeability_[phaseIdx] = kr; Valgrind::CheckDefined(relativePermeability_[phaseIdx]); } @@ -178,8 +179,8 @@ public: // capillary pressure parameters const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol); - wPhaseIdx_ = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol); - fluidState.setWettingPhase(wPhaseIdx_); + const auto wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol); + fluidState.setWettingPhase(wPhaseIdx); const auto& priVars = elemSol[scv.localDofIndex()]; const auto phasePresence = priVars.state(); @@ -213,17 +214,17 @@ public: // set pressures of the fluid phases using MaterialLaw = typename Problem::SpatialParams::MaterialLaw; - pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx_)); + pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx)); if (formulation == TwoPFormulation::p0s1) { fluidState.setPressure(liquidPhaseIdx, priVars[pressureIdx]); - fluidState.setPressure(gasPhaseIdx, (wPhaseIdx_ == liquidPhaseIdx) ? priVars[pressureIdx] + pc_ + fluidState.setPressure(gasPhaseIdx, (wPhaseIdx == liquidPhaseIdx) ? priVars[pressureIdx] + pc_ : priVars[pressureIdx] - pc_); } else { fluidState.setPressure(gasPhaseIdx, priVars[pressureIdx]); - fluidState.setPressure(liquidPhaseIdx, (wPhaseIdx_ == liquidPhaseIdx) ? priVars[pressureIdx] - pc_ + fluidState.setPressure(liquidPhaseIdx, (wPhaseIdx == liquidPhaseIdx) ? priVars[pressureIdx] - pc_ : priVars[pressureIdx] + pc_); } @@ -276,7 +277,7 @@ public: if (phasePresence == liquidPhaseOnly || phasePresence == gasPhaseOnly) fluidTemperature = priVars[switchIdx]; else if (phasePresence == twoPhases) - fluidTemperature = FluidSystem::vaporTemperature(fluidState, wPhaseIdx_); + fluidTemperature = FluidSystem::vaporTemperature(fluidState, fluidState_.wettingPhase()); else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid."); @@ -401,15 +402,14 @@ public: /*! * \brief Returns the wetting phase index */ - const int wettingPhaseIdx() const - { return wPhaseIdx_; } + int wettingPhase() const + { return fluidState_.wettingPhase(); } protected: FluidState fluidState_; SolidState solidState_; private: - int wPhaseIdx_; Scalar pc_; // The capillary pressure PermeabilityType permeability_; // Effective permeability within the control volume diff --git a/dumux/porousmediumflow/2p2c/volumevariables.hh b/dumux/porousmediumflow/2p2c/volumevariables.hh index fc5e2ebfe4c295452aa4d77727106094166a696f..d0aa354de2b91a670d767aa776349da5ded65aa6 100644 --- a/dumux/porousmediumflow/2p2c/volumevariables.hh +++ b/dumux/porousmediumflow/2p2c/volumevariables.hh @@ -147,11 +147,12 @@ public: using MaterialLaw = typename Problem::SpatialParams::MaterialLaw; const auto& matParams = problem.spatialParams().materialLawParams(element, scv, elemSol); - const int nPhaseIdx = 1 - wPhaseIdx_; + const int wPhaseIdx = fluidState_.wettingPhase(); + const int nPhaseIdx = 1 - wPhaseIdx; // relative permeabilities -> require wetting phase saturation as parameter! - relativePermeability_[wPhaseIdx_] = MaterialLaw::krw(matParams, saturation(wPhaseIdx_)); - relativePermeability_[nPhaseIdx] = MaterialLaw::krn(matParams, saturation(wPhaseIdx_)); + relativePermeability_[wPhaseIdx] = MaterialLaw::krw(matParams, saturation(wPhaseIdx)); + relativePermeability_[nPhaseIdx] = MaterialLaw::krn(matParams, saturation(wPhaseIdx)); // porosity & permeabilty updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps); @@ -206,8 +207,8 @@ public: using MaterialLaw = typename Problem::SpatialParams::MaterialLaw; const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol); - wPhaseIdx_ = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol); - fluidState.setWettingPhase(wPhaseIdx_); + const auto wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol); + fluidState.setWettingPhase(wPhaseIdx); // set the saturations if (phasePresence == firstPhaseOnly) @@ -237,17 +238,17 @@ public: DUNE_THROW(Dune::InvalidStateException, "Invalid phase presence."); // set pressures of the fluid phases - pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx_)); + pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx)); if (formulation == TwoPFormulation::p0s1) { fluidState.setPressure(phase0Idx, priVars[pressureIdx]); - fluidState.setPressure(phase1Idx, (wPhaseIdx_ == phase0Idx) ? priVars[pressureIdx] + pc_ + fluidState.setPressure(phase1Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] + pc_ : priVars[pressureIdx] - pc_); } else { fluidState.setPressure(phase1Idx, priVars[pressureIdx]); - fluidState.setPressure(phase0Idx, (wPhaseIdx_ == phase0Idx) ? priVars[pressureIdx] - pc_ + fluidState.setPressure(phase0Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] - pc_ : priVars[pressureIdx] + pc_); } } @@ -406,15 +407,13 @@ public: /*! * \brief Returns the wetting phase index */ - const int wettingPhaseIdx() const - { return wPhaseIdx_; } + int wettingPhase() const + { return fluidState_.wettingPhase(); } private: FluidState fluidState_; SolidState solidState_; - int wPhaseIdx_; - Scalar pc_; // The capillary pressure PermeabilityType permeability_; // Effective permeability within the control volume diff --git a/dumux/porousmediumflow/2pnc/volumevariables.hh b/dumux/porousmediumflow/2pnc/volumevariables.hh index 3efdd790b746aa63c0410c2015f91d29b453bfb4..fc19b96e245c2218424805c3b77adb3483935aaa 100644 --- a/dumux/porousmediumflow/2pnc/volumevariables.hh +++ b/dumux/porousmediumflow/2pnc/volumevariables.hh @@ -145,11 +145,12 @@ public: using MaterialLaw = typename Problem::SpatialParams::MaterialLaw; const auto& matParams = problem.spatialParams().materialLawParams(element, scv, elemSol); - const int nPhaseIdx = 1 - wPhaseIdx_; + const int wPhaseIdx = fluidState_.wettingPhase(); + const int nPhaseIdx = 1 - wPhaseIdx; // mobilities -> require wetting phase saturation as parameter! - mobility_[wPhaseIdx_] = MaterialLaw::krw(matParams, saturation(wPhaseIdx_))/fluidState_.viscosity(wPhaseIdx_); - mobility_[nPhaseIdx] = MaterialLaw::krn(matParams, saturation(wPhaseIdx_))/fluidState_.viscosity(nPhaseIdx); + mobility_[wPhaseIdx] = MaterialLaw::krw(matParams, saturation(wPhaseIdx))/fluidState_.viscosity(wPhaseIdx); + mobility_[nPhaseIdx] = MaterialLaw::krn(matParams, saturation(wPhaseIdx))/fluidState_.viscosity(nPhaseIdx); //update porosity before calculating the effective properties depending on it updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps); @@ -204,8 +205,8 @@ public: using MaterialLaw = typename Problem::SpatialParams::MaterialLaw; const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol); - wPhaseIdx_ = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol); - fluidState.setWettingPhase(wPhaseIdx_); + const auto wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol); + fluidState.setWettingPhase(wPhaseIdx); // set the saturations if (phasePresence == secondPhaseOnly) @@ -235,17 +236,17 @@ public: DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid."); // set pressures of the fluid phases - pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx_)); + pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx)); if (formulation == TwoPFormulation::p0s1) { fluidState.setPressure(phase0Idx, priVars[pressureIdx]); - fluidState.setPressure(phase1Idx, (wPhaseIdx_ == phase0Idx) ? priVars[pressureIdx] + pc_ + fluidState.setPressure(phase1Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] + pc_ : priVars[pressureIdx] - pc_); } else { fluidState.setPressure(phase1Idx, priVars[pressureIdx]); - fluidState.setPressure(phase0Idx, (wPhaseIdx_ == phase0Idx) ? priVars[pressureIdx] - pc_ + fluidState.setPressure(phase0Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] - pc_ : priVars[pressureIdx] + pc_); } @@ -495,8 +496,8 @@ public: /*! * \brief Returns the wetting phase index */ - const int wettingPhaseIdx() const - { return wPhaseIdx_; } + int wettingPhase() const + { return fluidState_.wettingPhase(); } protected: FluidState fluidState_; @@ -510,7 +511,6 @@ private: Scalar mobility_[ModelTraits::numFluidPhases()]; // Effective mobility within the control volume DiffusionCoefficients diffCoeff_; DiffusionCoefficients effectiveDiffCoeff_; - int wPhaseIdx_; }; diff --git a/dumux/porousmediumflow/co2/volumevariables.hh b/dumux/porousmediumflow/co2/volumevariables.hh index ea08710c4ff29488f7aa6d03b42775b88f505d27..4c6613cba7c9b05d4bc83decae5f24469f4c115b 100644 --- a/dumux/porousmediumflow/co2/volumevariables.hh +++ b/dumux/porousmediumflow/co2/volumevariables.hh @@ -136,11 +136,12 @@ public: using MaterialLaw = typename Problem::SpatialParams::MaterialLaw; const auto& matParams = problem.spatialParams().materialLawParams(element, scv, elemSol); - const int nPhaseIdx = 1 - wPhaseIdx_; + const int wPhaseIdx = fluidState_.wettingPhase(); + const int nPhaseIdx = 1 - wPhaseIdx; // relative permeabilities -> require wetting phase saturation as parameter! - relativePermeability_[wPhaseIdx_] = MaterialLaw::krw(matParams, saturation(wPhaseIdx_)); - relativePermeability_[nPhaseIdx] = MaterialLaw::krn(matParams, saturation(wPhaseIdx_)); + relativePermeability_[wPhaseIdx] = MaterialLaw::krw(matParams, saturation(wPhaseIdx)); + relativePermeability_[nPhaseIdx] = MaterialLaw::krn(matParams, saturation(wPhaseIdx)); // porosity & permeabilty updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps); @@ -197,8 +198,8 @@ public: using MaterialLaw = typename Problem::SpatialParams::MaterialLaw; const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol); - wPhaseIdx_ = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol); - fluidState.setWettingPhase(wPhaseIdx_); + const auto wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol); + fluidState.setWettingPhase(wPhaseIdx); // set the saturations if (phasePresence == secondPhaseOnly) @@ -228,17 +229,17 @@ public: DUNE_THROW(Dune::InvalidStateException, "Invalid phase presence."); // set pressures of the fluid phases - pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx_)); + pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx)); if (formulation == TwoPFormulation::p0s1) { fluidState.setPressure(phase0Idx, priVars[pressureIdx]); - fluidState.setPressure(phase1Idx, (wPhaseIdx_ == phase0Idx) ? priVars[pressureIdx] + pc_ + fluidState.setPressure(phase1Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] + pc_ : priVars[pressureIdx] - pc_); } else { fluidState.setPressure(phase1Idx, priVars[pressureIdx]); - fluidState.setPressure(phase0Idx, (wPhaseIdx_ == phase0Idx) ? priVars[pressureIdx] - pc_ + fluidState.setPressure(phase0Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] - pc_ : priVars[pressureIdx] + pc_); } @@ -491,12 +492,10 @@ public: /*! * \brief Returns the wetting phase index */ - const int wettingPhaseIdx() const - { return wPhaseIdx_; } - + int wettingPhase() const + { return fluidState_.wettingPhase(); } private: - int wPhaseIdx_; FluidState fluidState_; SolidState solidState_; Scalar pc_; // The capillary pressure