diff --git a/dumux/boxmodels/2p2c/2p2cvolumevariables.hh b/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
index 675772a839bd65301d56ac09f820534d9e2f98f2..f6f6b7583b5482aa6a56e261d33cabd899be7104 100644
--- a/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
+++ b/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
@@ -133,11 +133,61 @@ public:
scvIdx,
isOldSol);
- asImp_().updateTemperature_(priVars,
- element,
- elemGeom,
- scvIdx,
- problem);
+ completeFluidState(priVars, problem, element, elemGeom, scvIdx, fluidState_, isOldSol);
+
+ /////////////
+ // calculate the remaining quantities
+ /////////////
+ const MaterialLawParams &materialParams =
+ problem.spatialParameters().materialLawParams(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.updateAll(fluidState_);
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ // relative permeabilities
+ Scalar kr;
+ if (phaseIdx == lPhaseIdx)
+ kr = MaterialLaw::krw(materialParams, saturation(lPhaseIdx));
+ else // ATTENTION: krn requires the liquid saturation
+ // as parameter!
+ kr = MaterialLaw::krn(materialParams, saturation(lPhaseIdx));
+ relativePermeability_[phaseIdx] = kr;
+ Valgrind::CheckDefined(relativePermeability_[phaseIdx]);
+
+ // binary diffusion coefficents
+ diffCoeff_[phaseIdx] =
+ FluidSystem::binaryDiffusionCoefficient(fluidState_,
+ paramCache,
+ phaseIdx,
+ lCompIdx,
+ gCompIdx);
+ Valgrind::CheckDefined(diffCoeff_[phaseIdx]);
+ }
+
+ // porosity
+ porosity_ = problem.spatialParameters().porosity(element,
+ elemGeom,
+ scvIdx);
+ Valgrind::CheckDefined(porosity_);
+
+ // 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,
+ bool isOldSol = false)
+ {
+ Scalar t = Implementation::temperature_(priVars, problem, element,
+ elemGeom, scvIdx);
+ fluidState.setTemperature(t);
int globalVertIdx = problem.model().dofMapper().map(element, scvIdx, dim);
int phasePresence = problem.model().phasePresence(globalVertIdx, isOldSol);
@@ -159,8 +209,8 @@ public:
else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
}
else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
- fluidState_.setSaturation(lPhaseIdx, 1 - Sg);
- fluidState_.setSaturation(gPhaseIdx, Sg);
+ fluidState.setSaturation(lPhaseIdx, 1 - Sg);
+ fluidState.setSaturation(gPhaseIdx, Sg);
/////////////
// set the pressures of the fluid phases
@@ -172,12 +222,12 @@ public:
Scalar pC = MaterialLaw::pC(materialParams, 1 - Sg);
if (formulation == plSg) {
- fluidState_.setPressure(lPhaseIdx, priVars[pressureIdx]);
- fluidState_.setPressure(gPhaseIdx, priVars[pressureIdx] + pC);
+ fluidState.setPressure(lPhaseIdx, priVars[pressureIdx]);
+ fluidState.setPressure(gPhaseIdx, priVars[pressureIdx] + pC);
}
else if (formulation == pgSl) {
- fluidState_.setPressure(gPhaseIdx, priVars[pressureIdx]);
- fluidState_.setPressure(lPhaseIdx, priVars[pressureIdx] - pC);
+ fluidState.setPressure(gPhaseIdx, priVars[pressureIdx]);
+ fluidState.setPressure(lPhaseIdx, priVars[pressureIdx] - pC);
}
else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
@@ -192,7 +242,7 @@ public:
// result of the the gas <-> liquid equilibrium. This is
// the job of the "MiscibleMultiPhaseComposition"
// constraint solver
- MiscibleMultiPhaseComposition::solve(fluidState_,
+ MiscibleMultiPhaseComposition::solve(fluidState,
paramCache,
/*setViscosity=*/true,
/*setInternalEnergy=*/false);
@@ -214,13 +264,13 @@ public:
Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
// convert mass to mole fractions and set the fluid state
- fluidState_.setMoleFraction(gPhaseIdx, lCompIdx, massFractionG[lCompIdx]*avgMolarMass/M1);
- fluidState_.setMoleFraction(gPhaseIdx, gCompIdx, massFractionG[gCompIdx]*avgMolarMass/M2);
+ fluidState.setMoleFraction(gPhaseIdx, lCompIdx, massFractionG[lCompIdx]*avgMolarMass/M1);
+ fluidState.setMoleFraction(gPhaseIdx, gCompIdx, massFractionG[gCompIdx]*avgMolarMass/M2);
// calculate the composition of the remaining phases (as
// well as the densities of all phases). this is the job
// of the "ComputeFromReferencePhase" constraint solver
- ComputeFromReferencePhase::solve(fluidState_,
+ ComputeFromReferencePhase::solve(fluidState,
paramCache,
gPhaseIdx,
/*setViscosity=*/true,
@@ -242,51 +292,24 @@ public:
Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
// convert mass to mole fractions and set the fluid state
- fluidState_.setMoleFraction(lPhaseIdx, lCompIdx, massFractionL[lCompIdx]*avgMolarMass/M1);
- fluidState_.setMoleFraction(lPhaseIdx, gCompIdx, massFractionL[gCompIdx]*avgMolarMass/M2);
+ fluidState.setMoleFraction(lPhaseIdx, lCompIdx, massFractionL[lCompIdx]*avgMolarMass/M1);
+ fluidState.setMoleFraction(lPhaseIdx, gCompIdx, massFractionL[gCompIdx]*avgMolarMass/M2);
// calculate the composition of the remaining phases (as
// well as the densities of all phases). this is the job
// of the "ComputeFromReferencePhase" constraint solver
- ComputeFromReferencePhase::solve(fluidState_,
+ ComputeFromReferencePhase::solve(fluidState,
paramCache,
lPhaseIdx,
/*setViscosity=*/true,
/*setInternalEnergy=*/false);
}
-
- /////////////
- // calculate the remaining quantities
- /////////////
+
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
- // relative permeabilities
- Scalar kr;
- if (phaseIdx == lPhaseIdx)
- kr = MaterialLaw::krw(materialParams, saturation(lPhaseIdx));
- else // ATTENTION: krn requires the liquid saturation
- // as parameter!
- kr = MaterialLaw::krn(materialParams, saturation(lPhaseIdx));
- relativePermeability_[phaseIdx] = kr;
- Valgrind::CheckDefined(relativePermeability_[phaseIdx]);
-
- // binary diffusion coefficents
- diffCoeff_[phaseIdx] =
- FluidSystem::binaryDiffusionCoefficient(fluidState_,
- paramCache,
- phaseIdx,
- lCompIdx,
- gCompIdx);
- Valgrind::CheckDefined(diffCoeff_[phaseIdx]);
+ // compute and set the enthalpy
+ Scalar h = Implementation::enthalpy_(fluidState, paramCache, phaseIdx);
+ fluidState.setEnthalpy(phaseIdx, h);
}
-
- // porosity
- porosity_ = problem.spatialParameters().porosity(element,
- elemGeom,
- scvIdx);
- Valgrind::CheckDefined(porosity_);
-
- // energy related quantities
- asImp_().updateEnergy_(paramCache, priVars, problem, element, elemGeom, scvIdx, isOldSol);
}
/*!
@@ -383,21 +406,27 @@ public:
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)
+ {
+ return problem.boxTemperature(element, elemGeom, scvIdx);
+ }
+
+ template<class ParameterCache>
+ static Scalar enthalpy_(const FluidState& fluidState,
+ const ParameterCache& paramCache,
+ int phaseIdx)
{
- fluidState_.setTemperature(problem.boxTemperature(element, elemGeom, scvIdx));
+ return 0;
}
/*!
* \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/2p2cni/2p2cnivolumevariables.hh b/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
index b3d040276c35b3be0cbeebb952c55d6d72d890d5..7473eb5fd86ce0af5a293fef93fd81a18551c75d 100644
--- a/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
+++ b/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
@@ -60,6 +60,7 @@ class TwoPTwoCNIVolumeVariables : public TwoPTwoCVolumeVariables<TypeTag>
};
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename ParentType::FluidState FluidState;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)) };
enum { numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)) };
@@ -100,15 +101,21 @@ protected:
// is protected, we are friends with our parent..
friend class TwoPTwoCVolumeVariables<TypeTag>;
- // set the fluid state's temperature
- void updateTemperature_(const PrimaryVariables &sol,
+ static Scalar temperature_(const PrimaryVariables &priVars,
+ const Problem& problem,
const Element &element,
const FVElementGeometry &elemGeom,
- int scvIdx,
- const Problem &problem)
+ int scvIdx)
{
- // retrieve temperature from solution vector
- this->fluidState_.setTemperature(sol[temperatureIdx]);
+ return priVars[temperatureIdx];
+ }
+
+ template<class ParameterCache>
+ static Scalar enthalpy_(const FluidState& fluidState,
+ const ParameterCache& paramCache,
+ int phaseIdx)
+ {
+ return FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
}
/*!
@@ -121,22 +128,13 @@ protected:
* \param scvIdx The local index of the SCV (sub-control volume)
* \param isOldSol Evaluate function with solution of current or previous time step
*/
- template <class ParameterCache>
- void updateEnergy_(ParameterCache ¶mCache,
- const PrimaryVariables &sol,
+ void updateEnergy_(const PrimaryVariables &sol,
const Problem &problem,
const Element &element,
const FVElementGeometry &elemGeom,
int scvIdx,
bool isOldSol)
{
- // copmute and set the internal energies of the fluid phases
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
- Scalar h = FluidSystem::enthalpy(this->fluidState_, paramCache, phaseIdx);
-
- this->fluidState_.setEnthalpy(phaseIdx, h);
- }
-
// copmute and set the heat capacity of the solid phase
heatCapacity_ = problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
Valgrind::CheckDefined(heatCapacity_);