Skip to content

GitLab

Commit 7d3239ae authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'feature/wetting-phase' into 'master' 

Cleanup/wetting-phase-index

Closes #834

See merge request !1918
parents 8a2fd7d7 717c9e2a
Hide whitespace changes
Inline Side-by-side
dumux/io/plotthermalconductivitymodel.hh
View file @ 7d3239ae
......@@ -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_; }
......
dumux/material/fluidmatrixinteractions/2p/thermalconductivityjohansen.hh
View file @ 7d3239ae
......@@ -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();
......
dumux/material/fluidstates/compositional.hh
View file @ 7d3239ae
......@@ -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};
};
......
dumux/material/fluidstates/immiscible.hh
View file @ 7d3239ae
......@@ -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
......
dumux/material/fluidstates/isothermalimmiscible.hh
View file @ 7d3239ae
......@@ -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
......
dumux/porousmediumflow/2p/volumevariables.hh
View file @ 7d3239ae
......@@ -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_;
......
dumux/porousmediumflow/2p1c/volumevariables.hh
View file @ 7d3239ae
......@@ -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
......
dumux/porousmediumflow/2p2c/volumevariables.hh
View file @ 7d3239ae
......@@ -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
......
dumux/porousmediumflow/2pnc/volumevariables.hh
View file @ 7d3239ae
......@@ -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_;
};
......
dumux/porousmediumflow/co2/volumevariables.hh
View file @ 7d3239ae
......@@ -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
......
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