From aeba4e08159f0e16bc24c671130ca88d7f1e2e11 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Fri, 16 Dec 2011 13:23:31 +0000
Subject: [PATCH] Implemented completeFluidState for 1p, 1p2c, richards. Moved
convenience function to the base box model.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7065 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/boxmodels/1p/1pvolumevariables.hh | 65 ++++++++------
dumux/boxmodels/1p2c/1p2cvolumevariables.hh | 88 +++++++++++--------
dumux/boxmodels/2p/2pmodel.hh | 12 ---
dumux/boxmodels/2p/2pvolumevariables.hh | 2 +-
dumux/boxmodels/common/boxmodel.hh | 11 +++
.../richards/richardsvolumevariables.hh | 77 ++++++++--------
6 files changed, 140 insertions(+), 115 deletions(-)
diff --git a/dumux/boxmodels/1p/1pvolumevariables.hh b/dumux/boxmodels/1p/1pvolumevariables.hh
index 356d29d5dc..0f59889221 100644
--- a/dumux/boxmodels/1p/1pvolumevariables.hh
+++ b/dumux/boxmodels/1p/1pvolumevariables.hh
@@ -89,36 +89,43 @@ public:
{
ParentType::update(priVars, problem, element, elemGeom, scvIdx, isOldSol);
- // set the phase temperature and pressure
- asImp_().updateTemperature_(priVars,
- element,
- elemGeom,
- scvIdx,
- problem);
- fluidState_.setPressure(/*phaseIdx=*/0, priVars[Indices::pressureIdx]);
+ completeFluidState(priVars, problem, element, elemGeom, scvIdx, fluidState_);
+ // porosity
+ porosity_ = problem.spatialParameters().porosity(element,
+ elemGeom,
+ scvIdx);
+
+ // energy related quantities not contained in the fluid state
+ asImp_().updateEnergy_(priVars, problem, element, elemGeom, scvIdx, isOldSol);
+ };
+
+ static void completeFluidState(const PrimaryVariables& priVars,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& elemGeom,
+ int scvIdx,
+ FluidState& fluidState)
+ {
+ Scalar t = Implementation::temperature_(priVars, problem, element,
+ elemGeom, scvIdx);
+ fluidState.setTemperature(t);
+
+ fluidState.setPressure(/*phaseIdx=*/0, priVars[Indices::pressureIdx]);
// saturation in a single phase is always 1 and thus redundant
// to set. But since we use the fluid state shared by the
// immiscible multi-phase models, so we have to set it here...
- fluidState_.setSaturation(/*phaseIdx=*/0, 1.0);
+ fluidState.setSaturation(/*phaseIdx=*/0, 1.0);
typename FluidSystem::ParameterCache paramCache;
- paramCache.updatePhase(fluidState_, /*phaseIdx=*/0);
+ paramCache.updatePhase(fluidState, /*phaseIdx=*/0);
- Scalar value = FluidSystem::density(fluidState_, paramCache, /*phaseIdx=*/0);
- fluidState_.setDensity(/*phaseIdx=*/0, value);
+ Scalar value = FluidSystem::density(fluidState, paramCache, /*phaseIdx=*/0);
+ fluidState.setDensity(/*phaseIdx=*/0, value);
- value = FluidSystem::viscosity(fluidState_, paramCache, /*phaseIdx=*/0);
- fluidState_.setViscosity(/*phaseIdx=*/0, value);
-
- // porosity
- porosity_ = problem.spatialParameters().porosity(element,
- elemGeom,
- scvIdx);
-
- // energy related quantities
- asImp_().updateEnergy_(paramCache, priVars, problem, element, elemGeom, scvIdx, isOldSol);
- };
+ value = FluidSystem::viscosity(fluidState, paramCache, /*phaseIdx=*/0);
+ fluidState.setViscosity(/*phaseIdx=*/0, value);
+ }
/*!
* \brief Returns temperature inside the sub-control volume.
@@ -164,19 +171,19 @@ public:
{ return fluidState_; }
protected:
- void updateTemperature_(const PrimaryVariables &priVars,
+ static Scalar temperature_(const PrimaryVariables &priVars,
+ const Problem& problem,
const Element &element,
const FVElementGeometry &elemGeom,
- int scvIdx,
- const Problem &problem)
- { fluidState_.setTemperature(problem.boxTemperature(element, elemGeom, scvIdx)); }
+ int scvIdx)
+ {
+ return problem.boxTemperature(element, elemGeom, scvIdx);
+ }
/*!
* \brief Called by update() to compute the energy related quantities
*/
- template <class ParameterCache>
- void updateEnergy_(ParameterCache ¶mCache,
- const PrimaryVariables &sol,
+ void updateEnergy_(const PrimaryVariables &sol,
const Problem &problem,
const Element &element,
const FVElementGeometry &elemGeom,
diff --git a/dumux/boxmodels/1p2c/1p2cvolumevariables.hh b/dumux/boxmodels/1p2c/1p2cvolumevariables.hh
index 4d41328539..d7324bb9b4 100644
--- a/dumux/boxmodels/1p2c/1p2cvolumevariables.hh
+++ b/dumux/boxmodels/1p2c/1p2cvolumevariables.hh
@@ -99,13 +99,45 @@ public:
{
ParentType::update(priVars, problem, element, elemGeom, scvIdx, isOldSol);
- asImp_().updateTemperature_(priVars,
- element,
- elemGeom,
- scvIdx,
- problem);
+ completeFluidState(priVars, problem, element, elemGeom, scvIdx, fluidState_);
- fluidState_.setPressure(/*phaseIdx=*/0, priVars[pressureIdx]);
+ porosity_ = problem.spatialParameters().porosity(element, elemGeom, scvIdx);
+ tortuosity_ = problem.spatialParameters().tortuosity(element, elemGeom, scvIdx);
+ dispersivity_ = problem.spatialParameters().dispersivity(element, elemGeom, scvIdx);
+
+ // Second instance of a parameter cache.
+ // Could be avoided if diffusion coefficients also
+ // became part of the fluid state.
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState_, /*phaseIdx=*/0);
+
+ diffCoeff_ = FluidSystem::binaryDiffusionCoefficient(fluidState_,
+ paramCache,
+ phaseIdx,
+ comp0Idx,
+ comp1Idx);
+
+ Valgrind::CheckDefined(porosity_);
+ Valgrind::CheckDefined(tortuosity_);
+ Valgrind::CheckDefined(dispersivity_);
+ Valgrind::CheckDefined(diffCoeff_);
+
+ // energy related quantities not contained in the fluid state
+ asImp_().updateEnergy_(priVars, problem, element, elemGeom, scvIdx, isOldSol);
+ }
+
+ static void completeFluidState(const PrimaryVariables& priVars,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& elemGeom,
+ int scvIdx,
+ FluidState& fluidState)
+ {
+ Scalar t = Implementation::temperature_(priVars, problem, element,
+ elemGeom, scvIdx);
+ fluidState.setTemperature(t);
+
+ fluidState.setPressure(/*phaseIdx=*/0, priVars[pressureIdx]);
Scalar x1 = priVars[x1Idx]; //mole or mass fraction of component 1
if(!useMoles) //mass-fraction formulation
@@ -118,35 +150,17 @@ public:
x1 *= meanMolarMass/M1;
}
- fluidState_.setMoleFraction(/*phaseIdx=*/0, /*compIdx=*/0, 1 - x1);
- fluidState_.setMoleFraction(/*phaseIdx=*/0, /*compIdx=*/1, x1);
+ fluidState.setMoleFraction(/*phaseIdx=*/0, /*compIdx=*/0, 1 - x1);
+ fluidState.setMoleFraction(/*phaseIdx=*/0, /*compIdx=*/1, x1);
typename FluidSystem::ParameterCache paramCache;
- paramCache.updatePhase(fluidState_, /*phaseIdx=*/0);
+ paramCache.updatePhase(fluidState, /*phaseIdx=*/0);
Scalar value;
- value = FluidSystem::density(fluidState_, paramCache, /*phaseIdx=*/0);
- fluidState_.setDensity(/*phaseIdx=*/0, value);
- value = FluidSystem::viscosity(fluidState_, paramCache, /*phaseIdx=*/0);
- fluidState_.setViscosity(/*phaseIdx=*/0, value);
-
- porosity_ = problem.spatialParameters().porosity(element, elemGeom, scvIdx);
- tortuosity_ = problem.spatialParameters().tortuosity(element, elemGeom, scvIdx);
- dispersivity_ = problem.spatialParameters().dispersivity(element, elemGeom, scvIdx);
-
- diffCoeff_ = FluidSystem::binaryDiffusionCoefficient(fluidState_,
- paramCache,
- phaseIdx,
- comp0Idx,
- comp1Idx);
-
- Valgrind::CheckDefined(porosity_);
- Valgrind::CheckDefined(tortuosity_);
- Valgrind::CheckDefined(dispersivity_);
- Valgrind::CheckDefined(diffCoeff_);
-
- // energy related quantities
- asImp_().updateEnergy_(paramCache, priVars, problem, element, elemGeom, scvIdx, isOldSol);
+ value = FluidSystem::density(fluidState, paramCache, /*phaseIdx=*/0);
+ fluidState.setDensity(/*phaseIdx=*/0, value);
+ value = FluidSystem::viscosity(fluidState, paramCache, /*phaseIdx=*/0);
+ fluidState.setViscosity(/*phaseIdx=*/0, value);
}
/*!
@@ -236,21 +250,19 @@ public:
{ return porosity_; }
protected:
- void updateTemperature_(const PrimaryVariables &priVars,
+ static Scalar temperature_(const PrimaryVariables &priVars,
+ const Problem& problem,
const Element &element,
const FVElementGeometry &elemGeom,
- int scvIdx,
- const Problem &problem)
+ int scvIdx)
{
- fluidState_.setTemperature(problem.boxTemperature(element, elemGeom, scvIdx));
+ return problem.boxTemperature(element, elemGeom, scvIdx);
}
/*!
* \brief Called by update() to compute the energy related quantities
*/
- template <class ParameterCache>
- void updateEnergy_(ParameterCache ¶mCache,
- const PrimaryVariables &sol,
+ void updateEnergy_(const PrimaryVariables &sol,
const Problem &problem,
const Element &element,
const FVElementGeometry &elemGeom,
diff --git a/dumux/boxmodels/2p/2pmodel.hh b/dumux/boxmodels/2p/2pmodel.hh
index 9401466e43..4aba693c60 100644
--- a/dumux/boxmodels/2p/2pmodel.hh
+++ b/dumux/boxmodels/2p/2pmodel.hh
@@ -125,18 +125,6 @@ public:
return 1;
}
- template <class PrimaryVariables, class Problem, class Element, class FluidState>
- static void completeFluidState(const PrimaryVariables& primaryVariables,
- const Problem& problem,
- const Element& element,
- const FVElementGeometry& elementGeometry,
- int scvIdx,
- FluidState& fluidState)
- {
- VolumeVariables::completeFluidState(primaryVariables, problem, element,
- elementGeometry, scvIdx, fluidState);
- }
-
/*!
* \brief Append all quantities of interest which can be derived
* from the solution of the current time step to the VTK
diff --git a/dumux/boxmodels/2p/2pvolumevariables.hh b/dumux/boxmodels/2p/2pvolumevariables.hh
index e3b80aa2a0..49ab450ac9 100644
--- a/dumux/boxmodels/2p/2pvolumevariables.hh
+++ b/dumux/boxmodels/2p/2pvolumevariables.hh
@@ -102,7 +102,7 @@ public:
scvIdx,
isOldSol);
- completeFluidState(priVars, problem, element, elemGeom, scvIdx, fluidState_);
+ completeFluidState(priVars, problem, element, elemGeom, scvIdx, fluidState_);
const MaterialLawParams &materialParams =
problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
diff --git a/dumux/boxmodels/common/boxmodel.hh b/dumux/boxmodels/common/boxmodel.hh
index 70f5dd8ee4..a4bd9528fd 100644
--- a/dumux/boxmodels/common/boxmodel.hh
+++ b/dumux/boxmodels/common/boxmodel.hh
@@ -706,6 +706,17 @@ public:
const GridView &gridView() const
{ return problem_().gridView(); }
+ template <class FluidState>
+ static void completeFluidState(const PrimaryVariables& primaryVariables,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& elementGeometry,
+ int scvIdx,
+ FluidState& fluidState)
+ {
+ VolumeVariables::completeFluidState(primaryVariables, problem, element,
+ elementGeometry, scvIdx, fluidState);
+ }
protected:
/*!
* \brief A reference to the problem on which the model is applied.
diff --git a/dumux/boxmodels/richards/richardsvolumevariables.hh b/dumux/boxmodels/richards/richardsvolumevariables.hh
index 9336df54a3..558092781e 100644
--- a/dumux/boxmodels/richards/richardsvolumevariables.hh
+++ b/dumux/boxmodels/richards/richardsvolumevariables.hh
@@ -102,49 +102,58 @@ public:
scvIdx,
isOldSol);
- asImp_().updateTemperature_(priVars,
- element,
- elemGeom,
- scvIdx,
- problem);
-
+ completeFluidState(priVars, problem, element, elemGeom, scvIdx, fluidState_);
+
//////////
- // specify the fluid state
+ // specify the other parameters
//////////
+ const MaterialLawParams &matParams =
+ problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
+ relativePermeabilityWetting_ = MaterialLaw::krw(matParams,
+ fluidState_.saturation(wPhaseIdx));
+
+ porosity_ = problem.spatialParameters().porosity(element, elemGeom, scvIdx);
+
+ // energy related quantities not contained in the fluid state
+ asImp_().updateEnergy_(priVars, problem, element, elemGeom, scvIdx, isOldSol);
+ }
+
+ static void completeFluidState(const PrimaryVariables& priVars,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& elemGeom,
+ int scvIdx,
+ FluidState& fluidState)
+ {
+ // temperature
+ Scalar t = Implementation::temperature_(priVars, problem, element,
+ elemGeom, scvIdx);
+ fluidState.setTemperature(t);
// pressures
Scalar pnRef = problem.referencePressure(element, elemGeom, scvIdx);
const MaterialLawParams &matParams =
problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
Scalar minPc = MaterialLaw::pC(matParams, 1.0);
- fluidState_.setPressure(wPhaseIdx, priVars[pwIdx]);
- fluidState_.setPressure(nPhaseIdx, std::max(pnRef, priVars[pwIdx] + minPc));
+ fluidState.setPressure(wPhaseIdx, priVars[pwIdx]);
+ fluidState.setPressure(nPhaseIdx, std::max(pnRef, priVars[pwIdx] + minPc));
// saturations
- Scalar Sw = MaterialLaw::Sw(matParams, capillaryPressure());
- fluidState_.setSaturation(wPhaseIdx, Sw);
- fluidState_.setSaturation(nPhaseIdx, 1 - Sw);
+ Scalar Sw = MaterialLaw::Sw(matParams, fluidState.pressure(nPhaseIdx)
+ - fluidState.pressure(wPhaseIdx));
+ fluidState.setSaturation(wPhaseIdx, Sw);
+ fluidState.setSaturation(nPhaseIdx, 1 - Sw);
// density and viscosity
typename FluidSystem::ParameterCache paramCache;
- paramCache.updateAll(fluidState_);
- fluidState_.setDensity(wPhaseIdx, FluidSystem::density(fluidState_, paramCache, wPhaseIdx));
- fluidState_.setDensity(nPhaseIdx, 1e-10);
-
- fluidState_.setViscosity(wPhaseIdx, FluidSystem::viscosity(fluidState_, paramCache, wPhaseIdx));
- fluidState_.setViscosity(nPhaseIdx, 1e-10);
-
- //////////
- // specify the other parameters
- //////////
- relativePermeabilityWetting_ = MaterialLaw::krw(matParams, Sw);
+ paramCache.updateAll(fluidState);
+ fluidState.setDensity(wPhaseIdx, FluidSystem::density(fluidState, paramCache, wPhaseIdx));
+ fluidState.setDensity(nPhaseIdx, 1e-10);
- porosity_ = problem.spatialParameters().porosity(element, elemGeom, scvIdx);
-
- // energy related quantities
- asImp_().updateEnergy_(paramCache, priVars, problem, element, elemGeom, scvIdx, isOldSol);
+ fluidState.setViscosity(wPhaseIdx, FluidSystem::viscosity(fluidState, paramCache, wPhaseIdx));
+ fluidState.setViscosity(nPhaseIdx, 1e-10);
}
-
+
/*!
* \brief Return the fluid configuration at the given primary
* variables
@@ -248,21 +257,19 @@ public:
{ return fluidState_.pressure(nPhaseIdx) - fluidState_.pressure(wPhaseIdx); }
protected:
- void updateTemperature_(const PrimaryVariables &priVars,
+ static Scalar temperature_(const PrimaryVariables &priVars,
+ const Problem& problem,
const Element &element,
const FVElementGeometry &elemGeom,
- int scvIdx,
- const Problem &problem)
+ int scvIdx)
{
- fluidState_.setTemperature(problem.boxTemperature(element, elemGeom, scvIdx));
+ return problem.boxTemperature(element, elemGeom, scvIdx);
}
/*!
* \brief Called by update() to compute the energy related quantities
*/
- template <class ParameterCache>
- void updateEnergy_(ParameterCache ¶mCache,
- const PrimaryVariables &sol,
+ void updateEnergy_(const PrimaryVariables &sol,
const Problem &problem,
const Element &element,
const FVElementGeometry &elemGeom,
--
GitLab