diff --git a/dumux/boxmodels/2p2c/2p2cfluidstate.hh b/dumux/boxmodels/2p2c/2p2cfluidstate.hh
index 1233d106dc92a514783ad546f9b1e81412978087..5e689939b7961a1a3d657ea207bb21c665f2785d 100644
--- a/dumux/boxmodels/2p2c/2p2cfluidstate.hh
+++ b/dumux/boxmodels/2p2c/2p2cfluidstate.hh
@@ -26,7 +26,6 @@
* \brief Calculates the phase state from the primary variables in the
* 2p2c model.
*/
-#if 0
#ifndef DUMUX_2P2C_PHASE_STATE_HH
#define DUMUX_2P2C_PHASE_STATE_HH
@@ -368,4 +367,3 @@ public:
} // end namepace
#endif
-#endif
diff --git a/dumux/boxmodels/2p2c/2p2cfluxvariables.hh b/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
index 569179e8042be7ce068bf2b709931f72888bbb25..c48f64207ad25ae0c889f327baea85cb2122b0ad 100644
--- a/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
+++ b/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
@@ -244,7 +244,6 @@ private:
const Element &element,
const ElementVolumeVariables &elemDat)
{
-#if 0
const VolumeVariables &vDat_i = elemDat[face().i];
const VolumeVariables &vDat_j = elemDat[face().j];
@@ -261,6 +260,7 @@ private:
// calculate tortuosity at the nodes i and j needed
// for porous media diffusion coefficient
+
Scalar tau_i =
1.0/(vDat_i.porosity() * vDat_i.porosity()) *
pow(vDat_i.porosity() * vDat_i.saturation(phaseIdx), 7.0/3);
@@ -272,8 +272,8 @@ private:
// -> harmonic mean
porousDiffCoeff_[phaseIdx] = harmonicMean(vDat_i.porosity() * vDat_i.saturation(phaseIdx) * tau_i * vDat_i.diffCoeff(phaseIdx),
vDat_j.porosity() * vDat_j.saturation(phaseIdx) * tau_j * vDat_j.diffCoeff(phaseIdx));
+
}
-#endif
}
public:
@@ -306,13 +306,11 @@ public:
int downstreamIdx(int phaseIdx) const
{ return downstreamIdx_[phaseIdx]; }
-#if 0
/*!
* \brief The binary diffusion coefficient for each fluid phase.
*/
Scalar porousDiffCoeff(int phaseIdx) const
{ return porousDiffCoeff_[phaseIdx]; };
-#endif
/*!
* \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase at the integration
@@ -364,10 +362,8 @@ protected:
// local index of the downwind vertex for each phase
int downstreamIdx_[numPhases];
-/*
// the diffusion coefficient for the porous medium
Scalar porousDiffCoeff_[numPhases];
-*/
};
} // end namepace
diff --git a/dumux/boxmodels/2p2c/2p2clocalresidual.hh b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
index e1473a362a51b830641621f3b905b53ffc3ff3ea..64f5bf73a1af49605b77ffef5f7802f719f7c105 100644
--- a/dumux/boxmodels/2p2c/2p2clocalresidual.hh
+++ b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
@@ -74,6 +74,8 @@ protected:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
+ typedef TwoPTwoCFluidState<TypeTag> FluidState;
+
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
enum
@@ -205,9 +207,7 @@ public:
flux = 0;
asImp_()->computeAdvectiveFlux(flux, vars);
- Valgrind::CheckDefined(flux);
asImp_()->computeDiffusiveFlux(flux, vars);
- Valgrind::CheckDefined(flux);
}
/*!
@@ -248,13 +248,6 @@ public:
- mobilityUpwindAlpha) * (dn.density(phaseIdx)
* dn.mobility(phaseIdx) * dn.fluidState().massFrac(
phaseIdx, compIdx));
- Valgrind::CheckDefined(vars.KmvpNormal(phaseIdx));
- Valgrind::CheckDefined(up.density(phaseIdx));
- Valgrind::CheckDefined(up.mobility(phaseIdx));
- Valgrind::CheckDefined(up.fluidState().massFrac(phaseIdx, compIdx));
- Valgrind::CheckDefined(dn.density(phaseIdx));
- Valgrind::CheckDefined(dn.mobility(phaseIdx));
- Valgrind::CheckDefined(dn.fluidState().massFrac(phaseIdx, compIdx));
}
}
}
@@ -268,7 +261,6 @@ public:
*/
void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &vars) const
{
-#if 0
// add diffusive flux of gas component in liquid phase
Scalar tmp =
- vars.porousDiffCoeff(lPhaseIdx) *
@@ -284,7 +276,6 @@ public:
(vars.molarConcGrad(gPhaseIdx) * vars.face().normal);
flux[contiLEqIdx] += tmp * FluidSystem::molarMass(lCompIdx);
flux[contiGEqIdx] -= tmp * FluidSystem::molarMass(gCompIdx);
-#endif
}
/*!
diff --git a/dumux/boxmodels/2p2c/2p2cmodel.hh b/dumux/boxmodels/2p2c/2p2cmodel.hh
index 6cfd2ba56a98c727f86f4b63de18e46feec19f4f..4869d604ae0c62a721143c2acb43dc0febecf2c5 100644
--- a/dumux/boxmodels/2p2c/2p2cmodel.hh
+++ b/dumux/boxmodels/2p2c/2p2cmodel.hh
@@ -142,6 +142,8 @@ class TwoPTwoCModel: public BoxModel<TypeTag>
formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation))
};
+ typedef TwoPTwoCFluidState<TypeTag> FluidState;
+
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
@@ -368,7 +370,7 @@ public:
{
(*massFrac[phaseIdx][compIdx])[globalIdx]
= volVars.fluidState().massFrac(phaseIdx,
- compIdx);
+ compIdx);
Valgrind::CheckDefined(
(*massFrac[phaseIdx][compIdx])[globalIdx][0]);
diff --git a/dumux/boxmodels/2p2c/2p2cvolumevariables.hh b/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
index 801b394a2b4fbec697082cd794499462f008632e..d6da6e6e8398a5e0b2edfbfc905d2ac9796ff657 100644
--- a/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
+++ b/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
@@ -33,15 +33,11 @@
#include <dumux/common/math.hh>
#include <dune/common/collectivecommunication.hh>
-
-#include <dumux/material/fluidstates/equilibriumfluidstate.hh>
-#include <dumux/material/constraintsolvers/compositionfromfugacities.hh>
-
#include <vector>
#include <iostream>
#include "2p2cproperties.hh"
-#include "2p2cindices.hh"
+#include "2p2cfluidstate.hh"
namespace Dumux
{
@@ -73,34 +69,17 @@ class TwoPTwoCVolumeVariables : public BoxVolumeVariables<TypeTag>
numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation)),
- // component indices
lCompIdx = Indices::lCompIdx,
gCompIdx = Indices::gCompIdx,
- // phase indices
lPhaseIdx = Indices::lPhaseIdx,
- gPhaseIdx = Indices::gPhaseIdx,
-
- // primary variable indices
- pressureIdx = Indices::pressureIdx,
- switchIdx = Indices::switchIdx,
-
- // phase states
- lPhaseOnly = Indices::lPhaseOnly,
- gPhaseOnly = Indices::gPhaseOnly,
- bothPhases = Indices::bothPhases,
-
- // formulations
- plSg = TwoPTwoCFormulation::plSg,
- pgSl = TwoPTwoCFormulation::pgSl,
+ gPhaseIdx = Indices::gPhaseIdx
};
typedef typename GridView::template Codim<0>::Entity Element;
+ typedef TwoPTwoCFluidState<TypeTag> FluidState;
public:
- //! The return type of the fluidState() method
- typedef Dumux::EquilibriumFluidState<Scalar, FluidSystem> FluidState;
-
/*!
* \brief Update all quantities for a given control volume.
*
@@ -125,227 +104,62 @@ public:
scvIdx,
isOldSol);
- typename FluidSystem::MutableParameters mutParams;
- typename FluidSystem::MutableParameters::FluidState &fs
- = mutParams.fluidState();
-
- int globalVertIdx = problem.model().vertexMapper().map(element, scvIdx, dim);
- int phasePresence = problem.model().phasePresence(globalVertIdx, isOldSol);
-
- /////////////
- // set the phase saturations
- /////////////
- if (phasePresence == gPhaseOnly) {
- fs.setSaturation(lPhaseIdx, 0.0);
- fs.setSaturation(gPhaseIdx, 1.0);
- }
- else if (phasePresence == lPhaseOnly) {
- fs.setSaturation(lPhaseIdx, 1.0);
- fs.setSaturation(gPhaseIdx, 0.0);
- }
- else if (phasePresence == bothPhases) {
- Scalar Sg;
- if (formulation == plSg)
- Sg = priVars[switchIdx];
- else if (formulation == pgSl)
- Sg = 1.0 - priVars[switchIdx];
-
- fs.setSaturation(lPhaseIdx, 1 - Sg);
- fs.setSaturation(gPhaseIdx, Sg);
- }
-
- /////////////
- // set the phase temperatures
- /////////////
- // update the temperature part of the energy module
- Scalar T = asImp_().getTemperature(priVars,
- element,
- elemGeom,
- scvIdx,
- problem);
- Valgrind::CheckDefined(T);
- for (int i = 0; i < numPhases; ++i)
- fs.setTemperature(i, T);
-
- /////////////
- // set the phase pressures
- /////////////
+ asImp().updateTemperature_(priVars,
+ element,
+ elemGeom,
+ scvIdx,
+ problem);
// capillary pressure parameters
const MaterialLawParams &materialParams =
problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
- Scalar pC = MaterialLaw::pC(materialParams, fs.saturation(lPhaseIdx));
- if (formulation == plSg) {
- fs.setPressure(lPhaseIdx, priVars[pressureIdx]);
- fs.setPressure(gPhaseIdx, priVars[pressureIdx] + pC);
- }
- else if (formulation == pgSl){
- fs.setPressure(lPhaseIdx, priVars[pressureIdx] - pC);
- fs.setPressure(gPhaseIdx, priVars[pressureIdx]);
- }
- Valgrind::CheckDefined(fs.pressure(lPhaseIdx));
- Valgrind::CheckDefined(fs.pressure(gPhaseIdx));
-
- // update the mutable parameters for the pure components
- mutParams.updatePure(lPhaseIdx);
- mutParams.updatePure(gPhaseIdx);
-
- /////////////
- // set the phase compositions.
- /////////////
- if (phasePresence == gPhaseOnly) {
- // mass fractions
- Scalar Xg1 = priVars[switchIdx];
- Scalar Xg2 = 1 - Xg1;
-
- // molar masses
- Scalar M1 = FluidSystem::molarMass(lCompIdx);
- Scalar M2 = FluidSystem::molarMass(gCompIdx);
-
- // convert mass to mole fractions
- Scalar meanM = M1*M2/(M2 + Xg2*(M1 - M2));
- fs.setMoleFrac(gPhaseIdx, lCompIdx, Xg1 * M1/meanM);
- fs.setMoleFrac(gPhaseIdx, gCompIdx, Xg2 * M2/meanM);
- mutParams.updateMix(gPhaseIdx);
-
- // calculate component fugacities in gas phase
- Scalar fug1 = FluidSystem::computeFugacity(mutParams, gPhaseIdx, lCompIdx);
- Scalar fug2 = FluidSystem::computeFugacity(mutParams, gPhaseIdx, gCompIdx);
- fs.setFugacity(gPhaseIdx, lCompIdx, fug1);
- fs.setFugacity(gPhaseIdx, gCompIdx, fug2);
-
- // use same fugacities in liquid phase
- fs.setFugacity(lPhaseIdx, lCompIdx, fug1);
- fs.setFugacity(lPhaseIdx, gCompIdx, fug2);
-
- // initial guess of liquid composition
- fs.setMoleFrac(lPhaseIdx, lCompIdx, 0.98);
- fs.setMoleFrac(lPhaseIdx, gCompIdx, 0.02);
-
- // calculate liquid composition from fugacities
- typedef Dumux::CompositionFromFugacities<Scalar, FluidSystem> CompFromFug;
- CompFromFug::run(mutParams, lPhaseIdx);
-
- Valgrind::CheckDefined(fs.moleFrac(gPhaseIdx, lCompIdx));
- Valgrind::CheckDefined(fs.moleFrac(gPhaseIdx, gCompIdx));
- Valgrind::CheckDefined(fs.moleFrac(lPhaseIdx, lCompIdx));
- Valgrind::CheckDefined(fs.moleFrac(lPhaseIdx, gCompIdx));
-
- // calculate molar volume of gas phase
- Scalar Vmg = FluidSystem::computeMolarVolume(mutParams, gPhaseIdx);
- fs.setMolarVolume(gPhaseIdx, Vmg);
- }
- else if (phasePresence == lPhaseOnly) {
- // mass fractions
- Scalar Xl2 = priVars[switchIdx];
- Scalar Xl1 = 1 - Xl2;
-
- // molar masses
- Scalar M1 = FluidSystem::molarMass(lCompIdx);
- Scalar M2 = FluidSystem::molarMass(gCompIdx);
-
- // convert mass to mole fractions
- Scalar meanM = M1*M2/(M2 + Xl2*(M1 - M2));
- fs.setMoleFrac(lPhaseIdx, lCompIdx, Xl1 * M1/meanM);
- fs.setMoleFrac(lPhaseIdx, gCompIdx, Xl2 * M2/meanM);
- mutParams.updateMix(lPhaseIdx);
-
- // calculate component fugacities in liquid phase
- Scalar fug1 = FluidSystem::computeFugacity(mutParams, lPhaseIdx, lCompIdx);
- Scalar fug2 = FluidSystem::computeFugacity(mutParams, lPhaseIdx, gCompIdx);
- fs.setFugacity(lPhaseIdx, lCompIdx, fug1);
- fs.setFugacity(lPhaseIdx, gCompIdx, fug2);
-
- // use same fugacities in gas phase
- fs.setFugacity(gPhaseIdx, lCompIdx, fug1);
- fs.setFugacity(gPhaseIdx, gCompIdx, fug2);
-
- // initial guess of gas composition
- fs.setMoleFrac(gPhaseIdx, lCompIdx, 0.05);
- fs.setMoleFrac(gPhaseIdx, gCompIdx, 0.95);
-
- // calculate liquid composition from fugacities
- typedef Dumux::CompositionFromFugacities<Scalar, FluidSystem> CompFromFug;
- CompFromFug::run(mutParams, gPhaseIdx);
- Valgrind::CheckDefined(fs.moleFrac(gPhaseIdx, lCompIdx));
- Valgrind::CheckDefined(fs.moleFrac(gPhaseIdx, gCompIdx));
- Valgrind::CheckDefined(fs.moleFrac(lPhaseIdx, lCompIdx));
- Valgrind::CheckDefined(fs.moleFrac(lPhaseIdx, gCompIdx));
-
- // calculate molar volume of liquid phase
- Scalar Vml = FluidSystem::computeMolarVolume(mutParams, lPhaseIdx);
- fs.setMolarVolume(lPhaseIdx, Vml);
- }
- else if (phasePresence == bothPhases) {
- // HACK: assume both phases to be an ideal mixture,
- // i.e. the fugacity coefficents do not depend on the
- // composition
- Scalar phi_l1 = FluidSystem::computeFugacityCoeff(mutParams, lPhaseIdx, lCompIdx);
- Scalar phi_l2 = FluidSystem::computeFugacityCoeff(mutParams, lPhaseIdx, gCompIdx);
- Scalar phi_g1 = FluidSystem::computeFugacityCoeff(mutParams, gPhaseIdx, lCompIdx);
- Scalar phi_g2 = FluidSystem::computeFugacityCoeff(mutParams, gPhaseIdx, gCompIdx);
- Scalar pl = fs.pressure(lPhaseIdx);
- Scalar pg = fs.pressure(gPhaseIdx);
- Valgrind::CheckDefined(phi_l1);
- Valgrind::CheckDefined(phi_l2);
- Valgrind::CheckDefined(phi_g1);
- Valgrind::CheckDefined(phi_g2);
-
- Scalar xg2 = (phi_g2/phi_l1 - pl/pg) / (phi_g1/phi_l1 - phi_g2/phi_l2);
- Scalar xg1 = 1 - xg2;
- Scalar xl2 = (xg2*pg*phi_g2)/(pl*phi_l2);
- Scalar xl1 = 1 - xl2;
-
- fs.setMoleFrac(lPhaseIdx, lCompIdx, xl1);
- fs.setMoleFrac(lPhaseIdx, gCompIdx, xl2);
- fs.setMoleFrac(gPhaseIdx, lCompIdx, xg1);
- fs.setMoleFrac(gPhaseIdx, gCompIdx, xg2);
-
- mutParams.updateMix(lPhaseIdx);
- mutParams.updateMix(gPhaseIdx);
-
- Scalar Vml = FluidSystem::computeMolarVolume(mutParams, lPhaseIdx);
- Scalar Vmg = FluidSystem::computeMolarVolume(mutParams, gPhaseIdx);
- fs.setMolarVolume(lPhaseIdx, Vml);
- fs.setMolarVolume(gPhaseIdx, Vmg);
+
+ int globalVertIdx = problem.model().dofMapper().map(element, scvIdx, dim);
+ int phasePresence = problem.model().phasePresence(globalVertIdx, isOldSol);
+
+ // calculate phase state
+ fluidState_.update(priVars, materialParams, temperature(), phasePresence);
+ Valgrind::CheckDefined(fluidState_);
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ if (this->fluidState().saturation(phaseIdx) != 0.0) {
+ // Mobilities
+ const Scalar mu =
+ FluidSystem::phaseViscosity(phaseIdx,
+ fluidState().temperature(),
+ fluidState().phasePressure(lPhaseIdx),
+ fluidState());
+ 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));
+ mobility_[phaseIdx] = kr / mu;
+ Valgrind::CheckDefined(mobility_[phaseIdx]);
+
+ // binary diffusion coefficents
+ diffCoeff_[phaseIdx] =
+ FluidSystem::diffCoeff(phaseIdx,
+ lCompIdx,
+ gCompIdx,
+ fluidState_.temperature(),
+ fluidState_.phasePressure(phaseIdx),
+ fluidState_);
+ Valgrind::CheckDefined(diffCoeff_[phaseIdx]);
+ }
+ else {
+ mobility_[phaseIdx] = 0;
+ diffCoeff_[phaseIdx] = 0;
+ }
}
-
-
- // Mobilities
- Scalar muL = FluidSystem::computeViscosity(mutParams, lPhaseIdx);
- Scalar krL = MaterialLaw::krw(materialParams, fs.saturation(lPhaseIdx));
- mobility_[lPhaseIdx] = krL / muL;
- Valgrind::CheckDefined(mobility_[lPhaseIdx]);
-
- // ATTENTION: krn requires the liquid saturation as parameter!
- Scalar muG = FluidSystem::computeViscosity(mutParams, gPhaseIdx);
- Scalar krG = MaterialLaw::krn(materialParams, fs.saturation(lPhaseIdx));
- mobility_[gPhaseIdx] = krG / muG;
- Valgrind::CheckDefined(mobility_[gPhaseIdx]);
-
-#if 0
- // binary diffusion coefficents
- diffCoeff_[phaseIdx] =
- FluidSystem::diffCoeff(phaseIdx,
- lCompIdx,
- gCompIdx,
- fluidState_.temperature(),
- fluidState_.phasePressure(phaseIdx),
- fluidState_);
- Valgrind::CheckDefined(diffCoeff_[phaseIdx]);
-#endif
-
+
// porosity
porosity_ = problem.spatialParameters().porosity(element,
elemGeom,
scvIdx);
Valgrind::CheckDefined(porosity_);
-
- asImp_().updateEnergy(mutParams, priVars, element, elemGeom, scvIdx, problem);
-
- // assign the equilibrium fluid state from the generic one
- fluidState_.assign(fs);
- }
+ }
/*!
* \brief Returns the phase state for the control-volume.
@@ -378,7 +192,7 @@ public:
* \param phaseIdx The phase index
*/
Scalar molarDensity(int phaseIdx) const
- { return fluidState_.molarDensity(phaseIdx); }
+ { return fluidState_.density(phaseIdx) / fluidState_.averageMolarMass(phaseIdx); }
/*!
* \brief Returns the effective pressure of a given phase within
@@ -387,7 +201,7 @@ public:
* \param phaseIdx The phase index
*/
Scalar pressure(int phaseIdx) const
- { return fluidState_.pressure(phaseIdx); }
+ { return fluidState_.phasePressure(phaseIdx); }
/*!
* \brief Returns temperature inside the sub-control volume.
@@ -397,7 +211,7 @@ public:
* identical.
*/
Scalar temperature() const
- { return fluidState_.temperature(); }
+ { return temperature_; }
/*!
* \brief Returns the effective mobility of a given phase within
@@ -414,7 +228,7 @@ public:
* \brief Returns the effective capillary pressure within the control volume.
*/
Scalar capillaryPressure() const
- { return fluidState_.pressure(lPhaseIdx) - fluidState_.pressure(gPhaseIdx); }
+ { return fluidState_.capillaryPressure(); }
/*!
* \brief Returns the average porosity within the control volume.
@@ -422,35 +236,35 @@ public:
Scalar porosity() const
{ return porosity_; }
-#if 0
/*!
* \brief Returns the binary diffusion coefficients for a phase
*/
Scalar diffCoeff(int phaseIdx) const
{ return diffCoeff_[phaseIdx]; }
-#endif
+
protected:
- Scalar getTemperature_(const PrimaryVariables &priVars,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int scvIdx,
- const Problem &problem)
+ void updateTemperature_(const PrimaryVariables &priVars,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ const Problem &problem)
{
- return problem.temperature(element, elemGeom, scvIdx);
+ temperature_ = problem.temperature(element, elemGeom, scvIdx);
}
+ Scalar temperature_; //!< Temperature within the control volume
Scalar porosity_; //!< Effective porosity within the control volume
Scalar mobility_[numPhases]; //!< Effective mobility within the control volume
-// Scalar diffCoeff_[numPhases]; //!< Binary diffusion coefficients for the phases
+ Scalar diffCoeff_[numPhases]; //!< Binary diffusion coefficients for the phases
FluidState fluidState_;
private:
- Implementation &asImp_()
+ Implementation &asImp()
{ return *static_cast<Implementation*>(this); }
- const Implementation &asImp_() const
+ const Implementation &asImp() const
{ return *static_cast<const Implementation*>(this); }
};
diff --git a/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh b/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh
index 50d1d0b69860db45ddb32b209bbd2015393df3c3..e9abac3a8a5eb63af8c11a40827322e833f9ae46 100644
--- a/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh
+++ b/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh
@@ -112,11 +112,13 @@ public:
// compute the energy storage
result[energyEqIdx] =
vertDat.porosity()*(vertDat.density(lPhaseIdx) *
- vertDat.fluidState().internalEnergy(lPhaseIdx) *
+ vertDat.internalEnergy(lPhaseIdx) *
+ //vertDat.enthalpy(lPhaseIdx) *
vertDat.saturation(lPhaseIdx)
+
vertDat.density(gPhaseIdx) *
- vertDat.fluidState().internalEnergy(gPhaseIdx) *
+ vertDat.internalEnergy(gPhaseIdx) *
+ //vertDat.enthalpy(gPhaseIdx) *
vertDat.saturation(gPhaseIdx))
+
vertDat.temperature()*vertDat.heatCapacity();
@@ -150,12 +152,12 @@ public:
mobilityUpwindAlpha * // upstream vertex
( up.density(phase) *
up.mobility(phase) *
- up.fluidState().enthalpy(phase))
+ up.enthalpy(phase))
+
(1-mobilityUpwindAlpha) * // downstream vertex
( dn.density(phase) *
dn.mobility(phase) *
- dn.fluidState().enthalpy(phase)) );
+ dn.enthalpy(phase)) );
}
}
diff --git a/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh b/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
index 2a12e12e13fd60640b86057b9f4be3a96bfb5d31..4046591746d4545b50ab49ba71b4c523943fef34 100644
--- a/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
+++ b/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
@@ -70,19 +70,6 @@ class TwoPTwoCNIVolumeVariables : public TwoPTwoCVolumeVariables<TypeTag>
//! \endcond
public:
- /*!
- * \brief Update the temperature of the sub-control volume.
- */
- Scalar getTemperature(const PrimaryVariables &sol,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int scvIdx,
- const Problem &problem) const
- {
- // retrieve temperature from solution vector
- return sol[temperatureIdx];
- }
-
/*!
* \brief Update all quantities for a given control volume.
*
@@ -93,26 +80,61 @@ public:
* \param vertIdx The local index of the SCV (sub-control volume)
* \param isOldSol Evaluate function with solution of current or previous time step
*/
- template <class MutableParams>
- void updateEnergy(MutableParams &mutParams,
- const PrimaryVariables &sol,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int scvIdx,
- const Problem &problem)
+ void update(const PrimaryVariables &sol,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int vertIdx,
+ bool isOldSol)
{
- heatCapacity_ =
- problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
- Valgrind::CheckDefined(heatCapacity_);
+ // vertex update data for the mass balance
+ ParentType::update(sol,
+ problem,
+ element,
+ elemGeom,
+ vertIdx,
+ isOldSol);
// the internal energies and the enthalpies
- typename MutableParams::FluidState &fs = mutParams.fluidState();
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
- Scalar h =
- FluidSystem::computeEnthalpy(mutParams, phaseIdx);
- fs.setEnthalpy(phaseIdx, h);
+ if (this->fluidState().saturation(phaseIdx) != 0.0) {
+ enthalpy_[phaseIdx] =
+ FluidSystem::phaseEnthalpy(phaseIdx,
+ this->fluidState().temperature(),
+ this->fluidState().phasePressure(phaseIdx),
+ this->fluidState());
+ internalEnergy_[phaseIdx] =
+ FluidSystem::phaseInternalEnergy(phaseIdx,
+ this->fluidState().temperature(),
+ this->fluidState().phasePressure(phaseIdx),
+ this->fluidState());
+ }
+ else {
+ enthalpy_[phaseIdx] = 0;
+ internalEnergy_[phaseIdx] = 0;
+ }
}
- }
+ Valgrind::CheckDefined(internalEnergy_);
+ Valgrind::CheckDefined(enthalpy_);
+ };
+
+ /*!
+ * \brief Returns the total internal energy of a phase in the
+ * sub-control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar internalEnergy(int phaseIdx) const
+ { return internalEnergy_[phaseIdx]; };
+
+ /*!
+ * \brief Returns the total enthalpy of a phase in the sub-control
+ * volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar enthalpy(int phaseIdx) const
+ { return enthalpy_[phaseIdx]; };
/*!
* \brief Returns the total heat capacity \f$\mathrm{[J/(K*m^3]}\f$ of the rock matrix in
@@ -122,6 +144,27 @@ public:
{ return heatCapacity_; };
protected:
+ // this method gets called by the parent class. since this method
+ // is protected, we are friends with our parent..
+ friend class TwoPTwoCVolumeVariables<TypeTag>;
+ void updateTemperature_(const PrimaryVariables &sol,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ const Problem &problem)
+ {
+ // retrieve temperature from solution vector
+ this->temperature_ = sol[temperatureIdx];
+
+ heatCapacity_ =
+ problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
+
+ Valgrind::CheckDefined(this->temperature_);
+ Valgrind::CheckDefined(heatCapacity_);
+ }
+
+ Scalar internalEnergy_[numPhases];
+ Scalar enthalpy_[numPhases];
Scalar heatCapacity_;
};
diff --git a/dumux/boxmodels/common/boxassembler.hh b/dumux/boxmodels/common/boxassembler.hh
index 5e10a8d60819673c0e8e167e7707cf17f0565cea..2ebc31f1312279be302ae208d03b9996f3d8e7b6 100644
--- a/dumux/boxmodels/common/boxassembler.hh
+++ b/dumux/boxmodels/common/boxassembler.hh
@@ -135,10 +135,11 @@ public:
problemPtr_ = 0;
matrix_ = 0;
- // set reassemble accuracy to 0, so that if partial reassembly
- // of the jacobian matrix is disabled, the reassemble accuracy
- // is always smaller than the current relative tolerance
- reassembleAccuracy_ = 0.0;
+ // set reassemble tolerance to 0, so that if partial
+ // reassembly of the jacobian matrix is disabled, the
+ // reassemble tolerance is always smaller than the current
+ // relative tolerance
+ reassembleTolerance_ = 0.0;
}
~BoxAssembler()
@@ -238,13 +239,12 @@ public:
if (enablePartialReassemble) {
greenElems_ = gridView_().comm().sum(greenElems_);
- reassembleAccuracy_ = nextReassembleAccuracy_;
+ reassembleTolerance_ = nextReassembleTolerance_;
// print some information at the end of the iteration
problem_().newtonController().endIterMsg()
<< ", reassembled "
<< totalElems_ - greenElems_ << "/" << totalElems_
- << " (" << 100*Scalar(totalElems_ - greenElems_)/totalElems_ << "%) elems @accuracy="
- << reassembleAccuracy_;
+ << " (" << 100*Scalar(totalElems_ - greenElems_)/totalElems_ << "%) elems";
}
return;
@@ -270,7 +270,7 @@ public:
*/
void reassembleAll()
{
- nextReassembleAccuracy_ = 0.0;
+ nextReassembleTolerance_ = 0.0;
if (enablePartialReassemble) {
std::fill(vertexColor_.begin(),
@@ -286,17 +286,15 @@ public:
}
/*!
- * \brief Returns the largest relative error of a "green" vertex
- * for the most recent call of the assemble() method.
- *
- * This only has an effect if partial Jacobian reassembly is
- * enabled. If it is disabled, then this method always returns 0.
+ * \brief Returns the relative error below which a vertex is
+ * considered to be "green" if partial Jacobian reassembly
+ * is enabled.
*
* This returns the _actual_ relative computed seen by
* computeColors(), not the tolerance which it was given.
*/
- Scalar reassembleAccuracy() const
- { return reassembleAccuracy_; }
+ Scalar reassembleTolerance() const
+ { return reassembleTolerance_; }
/*!
* \brief Update the distance where the non-linear system was
@@ -360,16 +358,16 @@ public:
// mark the red vertices and update the tolerance of the
// linearization which actually will get achieved
- nextReassembleAccuracy_ = 0;
+ nextReassembleTolerance_ = 0;
for (int i = 0; i < vertexColor_.size(); ++i) {
vertexColor_[i] = Green;
- if (vertexDelta_[i] > relTol)
+ if (vertexDelta_[i] > relTol) {
// mark vertex as red if discrepancy is larger than
// the relative tolerance
vertexColor_[i] = Red;
- else
- nextReassembleAccuracy_ =
- std::max(nextReassembleAccuracy_, vertexDelta_[i]);
+ }
+ nextReassembleTolerance_ =
+ std::max(nextReassembleTolerance_, vertexDelta_[i]);
};
// Mark all red elements
@@ -732,8 +730,8 @@ private:
int totalElems_;
int greenElems_;
- Scalar nextReassembleAccuracy_;
- Scalar reassembleAccuracy_;
+ Scalar nextReassembleTolerance_;
+ Scalar reassembleTolerance_;
#if HAVE_DUNE_PDELAB
// PDELab stuff
diff --git a/dumux/boxmodels/common/boxlocalresidual.hh b/dumux/boxmodels/common/boxlocalresidual.hh
index 511eeae021d4e4c72336ddad11fc8df31004edb6..c9e8815d75edd2b0059c52c3e3f1855610c436b8 100644
--- a/dumux/boxmodels/common/boxlocalresidual.hh
+++ b/dumux/boxmodels/common/boxlocalresidual.hh
@@ -274,7 +274,7 @@ public:
Valgrind::CheckDefined(prevVolVars);
Valgrind::CheckDefined(curVolVars);
-#if 0 // HAVE_VALGRIND
+#if HAVE_VALGRIND
for (int i=0; i < fvGeom.numVertices; i++) {
Valgrind::CheckDefined(prevVolVars[i]);
Valgrind::CheckDefined(curVolVars[i]);
diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
index 4dfe3a9e04bea33d71a7db118b13a7a1d6da271e..4d85fab244026c37db0317aa8dd859f9916b0a09 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -512,10 +512,8 @@ public:
// compute the vertex and element colors for partial
// reassembly
if (enablePartialReassemble) {
- Scalar minReasmTol = 0.1*tolerance_;
- Scalar tmp = Dumux::geometricMean(error_, minReasmTol);
- Scalar reassembleTol = Dumux::geometricMean(error_, tmp);
- reassembleTol = std::max(reassembleTol, minReasmTol);
+ Scalar reassembleTol = Dumux::geometricMean(error_, 0.1*tolerance_);
+ reassembleTol = std::max(reassembleTol, 0.1*tolerance_);
this->model_().jacobianAssembler().updateDiscrepancy(uLastIter, deltaU);
this->model_().jacobianAssembler().computeColors(reassembleTol);
}
diff --git a/test/boxmodels/2p2cni/waterairproblem.hh b/test/boxmodels/2p2cni/waterairproblem.hh
index bd50e90580e0aaed4ebb357321f7c97f156e0887..d2b2faabaa3aef3140d5a72706346d574444cb2b 100644
--- a/test/boxmodels/2p2cni/waterairproblem.hh
+++ b/test/boxmodels/2p2cni/waterairproblem.hh
@@ -37,7 +37,6 @@
#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
#include <dumux/material/fluidsystems/h2o_n2_system.hh>
-#include <dumux/material/new_fluidsystems/h2on2fluidsystem.hh>
#include <dumux/boxmodels/2p2cni/2p2cnimodel.hh>
@@ -76,15 +75,11 @@ public:
// Set the problem property
SET_PROP(WaterAirProblem, Problem)
{
- typedef Dumux::WaterAirProblem<TypeTag> type;
+ typedef Dumux::WaterAirProblem<TTAG(WaterAirProblem)> type;
};
-// Set fluid configuration
-SET_PROP(WaterAirProblem, FluidSystem)
-{
-private: typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
-public: typedef Dumux::H2ON2FluidSystem<Scalar> type;
-};
+// Set the wetting phase
+SET_TYPE_PROP(WaterAirProblem, FluidSystem, Dumux::H2O_N2_System<TypeTag>);
// Set the spatial parameters
SET_TYPE_PROP(WaterAirProblem,
@@ -135,7 +130,7 @@ SET_BOOL_PROP(WaterAirProblem, NewtonWriteConvergence, false);
* To run the simulation execute the following line in shell:
* <tt>./test_2p2cni ./grids/test_2p2cni.dgf 300000 1000</tt>
* */
-template <class TypeTag>
+template <class TypeTag = TTAG(WaterAirProblem) >
class WaterAirProblem : public TwoPTwoCNIProblem<TypeTag>
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;