Skip to content

GitLab

Commit 717c9e2a authored by Timo Koch's avatar Timo Koch
Browse files

[volvars][2p...] Store wetting phase index in fluid state

parent 3e7ac93a
Hide whitespace changes
Inline Side-by-side
dumux/porousmediumflow/2p/volumevariables.hh
View file @ 717c9e2a
......@@ -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_);
}
}
......@@ -288,14 +290,13 @@ public:
* \brief Returns the wetting phase index
*/
int wettingPhase() const
{ return wPhaseIdx_; }
{ return fluidState_.wettingPhase(); }
protected:
FluidState fluidState_;
SolidState solidState_;
private:
int wPhaseIdx_;
Scalar pc_;
Scalar porosity_;
PermeabilityType permeability_;
......
dumux/porousmediumflow/2p1c/volumevariables.hh
View file @ 717c9e2a
......@@ -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.");
......@@ -402,14 +403,13 @@ public:
* \brief Returns the wetting phase index
*/
int wettingPhase() const
{ return wPhaseIdx_; }
{ return fluidState_.wettingPhase(); }
protected:
FluidState fluidState_;
SolidState solidState_;
private:
int wPhaseIdx_;
Scalar pc_; // The capillary pressure
PermeabilityType permeability_; // Effective permeability within the control volume
......
dumux/porousmediumflow/2p2c/volumevariables.hh
View file @ 717c9e2a
......@@ -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_);
}
}
......@@ -407,14 +408,12 @@ public:
* \brief Returns the wetting phase index
*/
int wettingPhase() const
{ return wPhaseIdx_; }
{ return fluidState_.wettingPhase(); }
private:
FluidState fluidState_;
SolidState solidState_;
int wPhaseIdx_;
Scalar pc_; // The capillary pressure
PermeabilityType permeability_; // Effective permeability within the control volume
......
dumux/porousmediumflow/2pnc/volumevariables.hh
View file @ 717c9e2a
......@@ -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_);
}
......@@ -496,7 +497,7 @@ public:
* \brief Returns the wetting phase index
*/
int wettingPhase() const
{ return wPhaseIdx_; }
{ 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_;
};
......
dumux/porousmediumflow/co2/volumevariables.hh
View file @ 717c9e2a
......@@ -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_);
}
......@@ -492,10 +493,9 @@ public:
* \brief Returns the wetting phase index
*/
int wettingPhase() const
{ return wPhaseIdx_; }
{ return fluidState_.wettingPhase(); }
private:
int wPhaseIdx_;
FluidState fluidState_;
SolidState solidState_;
Scalar pc_; // The capillary pressure
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment