From 5ab08062d80bb7d31817fefc5addaab1d901b511 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Tue, 4 Jan 2011 15:49:24 +0000
Subject: [PATCH] convert the 2p2c(ni) models to the new fluid systems
TODO: move code to calculate composition in twophase case to
constraint solver for composition from pressure + phase presence
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4949 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/boxmodels/2p2c/2p2cfluidstate.hh | 2 +
dumux/boxmodels/2p2c/2p2cfluxvariables.hh | 8 +-
dumux/boxmodels/2p2c/2p2clocalresidual.hh | 13 +-
dumux/boxmodels/2p2c/2p2cmodel.hh | 4 +-
dumux/boxmodels/2p2c/2p2cvolumevariables.hh | 314 ++++++++++++++----
dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh | 10 +-
.../boxmodels/2p2cni/2p2cnivolumevariables.hh | 99 ++----
dumux/boxmodels/common/boxassembler.hh | 42 +--
dumux/boxmodels/common/boxlocalresidual.hh | 2 +-
dumux/nonlinear/newtoncontroller.hh | 6 +-
test/boxmodels/2p2cni/waterairproblem.hh | 13 +-
11 files changed, 338 insertions(+), 175 deletions(-)
diff --git a/dumux/boxmodels/2p2c/2p2cfluidstate.hh b/dumux/boxmodels/2p2c/2p2cfluidstate.hh
index ae77441be8..28626ca8ac 100644
--- a/dumux/boxmodels/2p2c/2p2cfluidstate.hh
+++ b/dumux/boxmodels/2p2c/2p2cfluidstate.hh
@@ -21,6 +21,7 @@
* \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
@@ -362,3 +363,4 @@ public:
} // end namepace
#endif
+#endif
diff --git a/dumux/boxmodels/2p2c/2p2cfluxvariables.hh b/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
index 472854f0e9..780fcf8a25 100644
--- a/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
+++ b/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
@@ -239,6 +239,7 @@ private:
const Element &element,
const ElementVolumeVariables &elemDat)
{
+#if 0
const VolumeVariables &vDat_i = elemDat[face().i];
const VolumeVariables &vDat_j = elemDat[face().j];
@@ -255,7 +256,6 @@ 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);
@@ -267,8 +267,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:
@@ -301,11 +301,13 @@ 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
@@ -357,8 +359,10 @@ 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 b17bb1c219..2a21508b44 100644
--- a/dumux/boxmodels/2p2c/2p2clocalresidual.hh
+++ b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
@@ -69,8 +69,6 @@ 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
@@ -202,7 +200,9 @@ public:
flux = 0;
asImp_()->computeAdvectiveFlux(flux, vars);
+ Valgrind::CheckDefined(flux);
asImp_()->computeDiffusiveFlux(flux, vars);
+ Valgrind::CheckDefined(flux);
}
/*!
@@ -243,6 +243,13 @@ 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));
}
}
}
@@ -256,6 +263,7 @@ 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) *
@@ -271,6 +279,7 @@ 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 cf82487655..3ed9b171fa 100644
--- a/dumux/boxmodels/2p2c/2p2cmodel.hh
+++ b/dumux/boxmodels/2p2c/2p2cmodel.hh
@@ -137,8 +137,6 @@ 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;
@@ -365,7 +363,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 15532f15ad..64287cd154 100644
--- a/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
+++ b/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
@@ -28,11 +28,15 @@
#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 "2p2cfluidstate.hh"
+#include "2p2cindices.hh"
namespace Dumux
{
@@ -64,17 +68,34 @@ 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
+ 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,
};
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.
*
@@ -99,62 +120,227 @@ public:
scvIdx,
isOldSol);
- asImp().updateTemperature_(priVars,
- element,
- elemGeom,
- scvIdx,
- problem);
-
- // capillary pressure parameters
- const MaterialLawParams &materialParams =
- problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
+ typename FluidSystem::MutableParameters mutParams;
+ typename FluidSystem::MutableParameters::FluidState &fs
+ = mutParams.fluidState();
- int globalVertIdx = problem.model().dofMapper().map(element, scvIdx, dim);
+ int globalVertIdx = problem.model().vertexMapper().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;
- }
+ /////////////
+ // 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
+ /////////////
+
+ // 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);
+ }
+
+
+ // 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.
@@ -187,7 +373,7 @@ public:
* \param phaseIdx The phase index
*/
Scalar molarDensity(int phaseIdx) const
- { return fluidState_.density(phaseIdx) / fluidState_.averageMolarMass(phaseIdx); }
+ { return fluidState_.molarDensity(phaseIdx); }
/*!
* \brief Returns the effective pressure of a given phase within
@@ -196,7 +382,7 @@ public:
* \param phaseIdx The phase index
*/
Scalar pressure(int phaseIdx) const
- { return fluidState_.phasePressure(phaseIdx); }
+ { return fluidState_.pressure(phaseIdx); }
/*!
* \brief Returns temperature inside the sub-control volume.
@@ -206,7 +392,7 @@ public:
* identical.
*/
Scalar temperature() const
- { return temperature_; }
+ { return fluidState_.temperature(); }
/*!
* \brief Returns the effective mobility of a given phase within
@@ -223,7 +409,7 @@ public:
* \brief Returns the effective capillary pressure within the control volume.
*/
Scalar capillaryPressure() const
- { return fluidState_.capillaryPressure(); }
+ { return fluidState_.pressure(lPhaseIdx) - fluidState_.pressure(gPhaseIdx); }
/*!
* \brief Returns the average porosity within the control volume.
@@ -231,35 +417,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:
- void updateTemperature_(const PrimaryVariables &priVars,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int scvIdx,
- const Problem &problem)
+ Scalar getTemperature_(const PrimaryVariables &priVars,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ const Problem &problem)
{
- temperature_ = problem.temperature(element, elemGeom, scvIdx);
+ return 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 7e63800945..1f9b75d568 100644
--- a/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh
+++ b/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh
@@ -107,13 +107,11 @@ public:
// compute the energy storage
result[energyEqIdx] =
vertDat.porosity()*(vertDat.density(lPhaseIdx) *
- vertDat.internalEnergy(lPhaseIdx) *
- //vertDat.enthalpy(lPhaseIdx) *
+ vertDat.fluidState().internalEnergy(lPhaseIdx) *
vertDat.saturation(lPhaseIdx)
+
vertDat.density(gPhaseIdx) *
- vertDat.internalEnergy(gPhaseIdx) *
- //vertDat.enthalpy(gPhaseIdx) *
+ vertDat.fluidState().internalEnergy(gPhaseIdx) *
vertDat.saturation(gPhaseIdx))
+
vertDat.temperature()*vertDat.heatCapacity();
@@ -147,12 +145,12 @@ public:
mobilityUpwindAlpha * // upstream vertex
( up.density(phase) *
up.mobility(phase) *
- up.enthalpy(phase))
+ up.fluidState().enthalpy(phase))
+
(1-mobilityUpwindAlpha) * // downstream vertex
( dn.density(phase) *
dn.mobility(phase) *
- dn.enthalpy(phase)) );
+ dn.fluidState().enthalpy(phase)) );
}
}
diff --git a/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh b/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
index 7e222fedb8..a6c158ca5f 100644
--- a/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
+++ b/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
@@ -65,6 +65,19 @@ 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.
*
@@ -75,61 +88,26 @@ public:
* \param vertIdx The local index of the SCV (sub-control volume)
* \param isOldSol Evaluate function with solution of current or previous time step
*/
- void update(const PrimaryVariables &sol,
- const Problem &problem,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int vertIdx,
- bool isOldSol)
+ template <class MutableParams>
+ void updateEnergy(MutableParams &mutParams,
+ const PrimaryVariables &sol,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ const Problem &problem)
{
- // vertex update data for the mass balance
- ParentType::update(sol,
- problem,
- element,
- elemGeom,
- vertIdx,
- isOldSol);
+ heatCapacity_ =
+ problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
+ Valgrind::CheckDefined(heatCapacity_);
// the internal energies and the enthalpies
+ typename MutableParams::FluidState &fs = mutParams.fluidState();
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
- 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;
- }
+ Scalar h =
+ FluidSystem::computeEnthalpy(mutParams, phaseIdx);
+ fs.setEnthalpy(phaseIdx, h);
}
- 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
@@ -139,27 +117,6 @@ 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 a64c827693..29d57529ef 100644
--- a/dumux/boxmodels/common/boxassembler.hh
+++ b/dumux/boxmodels/common/boxassembler.hh
@@ -130,11 +130,10 @@ public:
problemPtr_ = 0;
matrix_ = 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;
+ // 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;
}
~BoxAssembler()
@@ -234,12 +233,13 @@ public:
if (enablePartialReassemble) {
greenElems_ = gridView_().comm().sum(greenElems_);
- reassembleTolerance_ = nextReassembleTolerance_;
+ reassembleAccuracy_ = nextReassembleAccuracy_;
// print some information at the end of the iteration
problem_().newtonController().endIterMsg()
<< ", reassembled "
<< totalElems_ - greenElems_ << "/" << totalElems_
- << " (" << 100*Scalar(totalElems_ - greenElems_)/totalElems_ << "%) elems";
+ << " (" << 100*Scalar(totalElems_ - greenElems_)/totalElems_ << "%) elems @accuracy="
+ << reassembleAccuracy_;
}
return;
@@ -265,7 +265,7 @@ public:
*/
void reassembleAll()
{
- nextReassembleTolerance_ = 0.0;
+ nextReassembleAccuracy_ = 0.0;
if (enablePartialReassemble) {
std::fill(vertexColor_.begin(),
@@ -281,15 +281,17 @@ public:
}
/*!
- * \brief Returns the relative error below which a vertex is
- * considered to be "green" if partial Jacobian reassembly
- * is enabled.
+ * \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.
*
* This returns the _actual_ relative computed seen by
* computeColors(), not the tolerance which it was given.
*/
- Scalar reassembleTolerance() const
- { return reassembleTolerance_; }
+ Scalar reassembleAccuracy() const
+ { return reassembleAccuracy_; }
/*!
* \brief Update the distance where the non-linear system was
@@ -353,16 +355,16 @@ public:
// mark the red vertices and update the tolerance of the
// linearization which actually will get achieved
- nextReassembleTolerance_ = 0;
+ nextReassembleAccuracy_ = 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;
- }
- nextReassembleTolerance_ =
- std::max(nextReassembleTolerance_, vertexDelta_[i]);
+ else
+ nextReassembleAccuracy_ =
+ std::max(nextReassembleAccuracy_, vertexDelta_[i]);
};
// Mark all red elements
@@ -725,8 +727,8 @@ private:
int totalElems_;
int greenElems_;
- Scalar nextReassembleTolerance_;
- Scalar reassembleTolerance_;
+ Scalar nextReassembleAccuracy_;
+ Scalar reassembleAccuracy_;
#if HAVE_DUNE_PDELAB
// PDELab stuff
diff --git a/dumux/boxmodels/common/boxlocalresidual.hh b/dumux/boxmodels/common/boxlocalresidual.hh
index 6892fb7b68..71302b8805 100644
--- a/dumux/boxmodels/common/boxlocalresidual.hh
+++ b/dumux/boxmodels/common/boxlocalresidual.hh
@@ -269,7 +269,7 @@ public:
Valgrind::CheckDefined(prevVolVars);
Valgrind::CheckDefined(curVolVars);
-#if HAVE_VALGRIND
+#if 0 // 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 843030e552..ccf77d366a 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -507,8 +507,10 @@ public:
// compute the vertex and element colors for partial
// reassembly
if (enablePartialReassemble) {
- Scalar reassembleTol = Dumux::geometricMean(error_, 0.1*tolerance_);
- reassembleTol = std::max(reassembleTol, 0.1*tolerance_);
+ Scalar minReasmTol = 0.1*tolerance_;
+ Scalar tmp = Dumux::geometricMean(error_, minReasmTol);
+ Scalar reassembleTol = Dumux::geometricMean(error_, tmp);
+ reassembleTol = std::max(reassembleTol, minReasmTol);
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 9a3a1bacf0..63f4b26519 100644
--- a/test/boxmodels/2p2cni/waterairproblem.hh
+++ b/test/boxmodels/2p2cni/waterairproblem.hh
@@ -32,6 +32,7 @@
#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>
@@ -70,11 +71,15 @@ public:
// Set the problem property
SET_PROP(WaterAirProblem, Problem)
{
- typedef Dumux::WaterAirProblem<TTAG(WaterAirProblem)> type;
+ typedef Dumux::WaterAirProblem<TypeTag> type;
};
-// Set the wetting phase
-SET_TYPE_PROP(WaterAirProblem, FluidSystem, Dumux::H2O_N2_System<TypeTag>);
+// 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 spatial parameters
SET_TYPE_PROP(WaterAirProblem,
@@ -125,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 = TTAG(WaterAirProblem) >
+template <class TypeTag>
class WaterAirProblem : public TwoPTwoCNIProblem<TypeTag>
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
--
GitLab