Commit 5ec5aaa0 authored by Andreas Lauser's avatar Andreas Lauser
Browse files

new fluid systems: implement the changes discussed last friday

If you have any major objections, speak at meeting on friday or
forever hold your peace.  (Where "forever" means "until the next major
fluid framework refactoring". [Which won't be initiated by me.])

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@6807 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 9cc71304
......@@ -67,7 +67,7 @@ NEW_PROP_TAG(MPNCVtkAddVarPressures);
NEW_PROP_TAG(MPNCVtkAddVelocities);
NEW_PROP_TAG(MPNCVtkAddDensities);
NEW_PROP_TAG(MPNCVtkAddMobilities);
NEW_PROP_TAG(MPNCVtkAddMeanMolarMass);
NEW_PROP_TAG(MPNCVtkAddAverageMolarMass);
NEW_PROP_TAG(MPNCVtkAddMassFractions);
NEW_PROP_TAG(MPNCVtkAddMoleFractions);
NEW_PROP_TAG(MPNCVtkAddMolarities);
......
......@@ -214,7 +214,7 @@ SET_BOOL_PROP(BoxMPNC, MPNCVtkAddVarPressures, false);
SET_BOOL_PROP(BoxMPNC, MPNCVtkAddVelocities, false);
SET_BOOL_PROP(BoxMPNC, MPNCVtkAddDensities, true);
SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMobilities, true);
SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMeanMolarMass, false);
SET_BOOL_PROP(BoxMPNC, MPNCVtkAddAverageMolarMass, false);
SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMassFractions, false);
SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMoleFractions, true);
SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMolarities, false);
......
......@@ -90,7 +90,6 @@ class MPNCVolumeVariables
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
public:
//! The return type of the fluidState() method
typedef typename MassVolumeVariables::FluidState FluidState;
......@@ -124,8 +123,9 @@ public:
elemGeom,
scvIdx,
isOldSol);
ParentType::checkDefined();
typename FluidSystem::MutableParameters mutParams;
typename FluidSystem::ParameterCache paramCache;
/////////////
// set the phase saturations
......@@ -133,24 +133,21 @@ public:
Scalar sumSat = 0;
for (int i = 0; i < numPhases - 1; ++i) {
sumSat += priVars[S0Idx + i];
mutParams.setSaturation(i, priVars[S0Idx + i]);
fluidState_.setSaturation(i, priVars[S0Idx + i]);
}
Valgrind::CheckDefined(sumSat);
mutParams.setSaturation(numPhases - 1, 1.0 - sumSat);
fluidState_.setSaturation(numPhases - 1, 1.0 - sumSat);
/////////////
// set the fluid phase temperatures
/////////////
// update the temperature part of the energy module
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx){
Scalar T = EnergyVolumeVariables::getTemperature(priVars,
element,
elemGeom,
scvIdx,
problem,
phaseIdx);
mutParams.setTemperature(phaseIdx, T);
}
EnergyVolumeVariables::updateTemperatures(fluidState_,
paramCache,
priVars,
element,
elemGeom,
scvIdx,
problem);
/////////////
// set the phase pressures
......@@ -161,16 +158,17 @@ public:
problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
// capillary pressures
Scalar capPress[numPhases];
MaterialLaw::capillaryPressures(capPress, materialParams, mutParams);
MaterialLaw::capillaryPressures(capPress, materialParams, fluidState_);
// add to the pressure of the first fluid phase
Scalar p0 = priVars[p0Idx];
for (int i = 0; i < numPhases; ++ i)
mutParams.setPressure(i, p0 - capPress[0] + capPress[i]);
fluidState_.setPressure(i, p0 - capPress[0] + capPress[i]);
/////////////
// set the fluid compositions
/////////////
MassVolumeVariables::update(mutParams,
MassVolumeVariables::update(fluidState_,
paramCache,
priVars,
hint_,
problem,
......@@ -178,7 +176,6 @@ public:
elemGeom,
scvIdx);
MassVolumeVariables::checkDefined();
ParentType::checkDefined();
/////////////
// Porosity
......@@ -189,7 +186,6 @@ public:
elemGeom,
scvIdx);
Valgrind::CheckDefined(porosity_);
ParentType::checkDefined();
/////////////
// Phase mobilities
......@@ -198,31 +194,30 @@ public:
// relative permeabilities
MaterialLaw::relativePermeabilities(relativePermeability_,
materialParams,
mutParams);
fluidState_);
// dynamic viscosities
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
// viscosities
Scalar mu = FluidSystem::computeViscosity(mutParams, phaseIdx);
mutParams.setViscosity(phaseIdx, mu);
Scalar mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
fluidState_.setViscosity(phaseIdx, mu);
}
ParentType::checkDefined();
/////////////
// diffusion
/////////////
// update the diffusion part of the volume data
DiffusionVolumeVariables::update(mutParams, *this, problem);
DiffusionVolumeVariables::update(fluidState_, paramCache, *this, problem);
DiffusionVolumeVariables::checkDefined();
ParentType::checkDefined();
/////////////
// energy
/////////////
// update the remaining parts of the energy module
EnergyVolumeVariables::update(mutParams,
EnergyVolumeVariables::update(fluidState_,
paramCache,
priVars,
element,
elemGeom,
......@@ -230,18 +225,16 @@ public:
problem);
EnergyVolumeVariables::checkDefined();
ParentType::checkDefined();
// assign the fluid state to be stored (needed in the next update)
fluidState_.assign(mutParams);
// make sure the quantities in the fluid state are well-defined
fluidState_.checkDefined();
// specific interfacial area,
// well also all the dimensionless numbers :-)
// well, also the mass transfer rate
IAVolumeVariables::update(*this,
mutParams,
fluidState_,
paramCache,
priVars,
problem,
element,
......@@ -249,11 +242,7 @@ public:
scvIdx);
IAVolumeVariables::checkDefined();
ParentType::checkDefined();
checkDefined();
ParentType::checkDefined();
}
/*!
......@@ -330,7 +319,7 @@ public:
// difference of sum of mole fractions in the phase from 100%
Scalar a = 1;
for (int i = 0; i < numComponents; ++i)
a -= fluidState.moleFrac(phaseIdx, i);
a -= fluidState.moleFraction(phaseIdx, i);
return a;
}
......
......@@ -61,9 +61,10 @@ public:
/*!
* \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
*/
template <class MutableParams>
template <class FluidState, class ParameterCache>
void update(const VolumeVariables & volVars,
const MutableParams & mutParams,
const FluidState &fluidState,
const ParameterCache &paramCache,
const PrimaryVariables &priVars,
const Problem &problem,
const Element & element,
......
......@@ -85,7 +85,7 @@ public:
velocityOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddVelocities);
densityOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddDensities);
mobilityOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMobilities);
meanMolarMassOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMeanMolarMass);
averageMolarMassOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddAverageMolarMass);
massFracOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMassFractions);
moleFracOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMoleFractions);
molarityOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMolarities);
......@@ -116,7 +116,7 @@ public:
if (pressureOutput_) this->resizePhaseBuffer_(pressure_);
if (densityOutput_) this->resizePhaseBuffer_(density_);
if (mobilityOutput_) this->resizePhaseBuffer_(mobility_);
if (meanMolarMassOutput_) this->resizePhaseBuffer_(meanMolarMass_);
if (averageMolarMassOutput_) this->resizePhaseBuffer_(averageMolarMass_);
if (moleFracOutput_) this->resizePhaseComponentBuffer_(moleFrac_);
if (massFracOutput_) this->resizePhaseComponentBuffer_(massFrac_);
......@@ -154,10 +154,10 @@ public:
if (pressureOutput_) pressure_[phaseIdx][I] = volVars.fluidState().pressure(phaseIdx);
if (densityOutput_) density_[phaseIdx][I] = volVars.fluidState().density(phaseIdx);
if (mobilityOutput_) mobility_[phaseIdx][I] = volVars.mobility(phaseIdx);
if (meanMolarMassOutput_) meanMolarMass_[phaseIdx][I] = volVars.fluidState().meanMolarMass(phaseIdx);
if (averageMolarMassOutput_) averageMolarMass_[phaseIdx][I] = volVars.fluidState().averageMolarMass(phaseIdx);
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
if (moleFracOutput_) moleFrac_[phaseIdx][compIdx][I] = volVars.fluidState().moleFrac(phaseIdx, compIdx);
if (massFracOutput_) massFrac_[phaseIdx][compIdx][I] = volVars.fluidState().massFrac(phaseIdx, compIdx);
if (moleFracOutput_) moleFrac_[phaseIdx][compIdx][I] = volVars.fluidState().moleFraction(phaseIdx, compIdx);
if (massFracOutput_) massFrac_[phaseIdx][compIdx][I] = volVars.fluidState().massFraction(phaseIdx, compIdx);
if (molarityOutput_) molarity_[phaseIdx][compIdx][I] = volVars.fluidState().molarity(phaseIdx, compIdx);
}
}
......@@ -213,8 +213,8 @@ public:
if (densityOutput_)
this->commitPhaseBuffer_(writer, "rho_%s", density_);
if (meanMolarMassOutput_)
this->commitPhaseBuffer_(writer, "M_%s", meanMolarMass_);
if (averageMolarMassOutput_)
this->commitPhaseBuffer_(writer, "M_%s", averageMolarMass_);
if (mobilityOutput_)
this->commitPhaseBuffer_(writer, "lambda_%s", mobility_);
......@@ -263,13 +263,13 @@ private:
bool massFracOutput_;
bool moleFracOutput_;
bool molarityOutput_;
bool meanMolarMassOutput_;
bool averageMolarMassOutput_;
PhaseBuffer saturation_;
PhaseBuffer pressure_;
PhaseBuffer density_;
PhaseBuffer mobility_;
PhaseBuffer meanMolarMass_;
PhaseBuffer averageMolarMass_;
PhaseVelocityField velocity_;
ScalarBuffer boxSurface_;
......
......@@ -129,18 +129,17 @@ class MPNCVolumeVariablesDiffusion<TypeTag, false>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename FluidSystem::MutableParameters MutableParameters;
public:
MPNCVolumeVariablesDiffusion()
{}
void update(MutableParameters &mutParams,
const VolumeVariables &volVars,
template <class FluidState, class ParameterCache>
void update(FluidState &fluidState,
ParameterCache &paramCache,
const VolumeVariables &volVars,
const Problem &problem)
{
};
{ };
Scalar diffCoeffL(int compIdx) const
{ return 0; };
......
......@@ -27,6 +27,8 @@
#ifndef DUMUX_MPNC_ENERGY_VOLUME_VARIABLES_HH
#define DUMUX_MPNC_ENERGY_VOLUME_VARIABLES_HH
#include <dumux/material/MpNcfluidstates/equilibriumfluidstate.hh>
namespace Dumux
{
/*!
......@@ -57,19 +59,24 @@ class MPNCVolumeVariablesEnergy
enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)) };
//typedef typename GET_PROP_TYPE(TypeTag, PTAG(MPNCEnergyIndices)) EnergyIndices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename FluidSystem::ParameterCache ParameterCache;
typedef Dumux::EquilibriumFluidState<Scalar, FluidSystem> FluidState;
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 int temperatureIdx) const
void updateTemperatures(FluidState &fs,
ParameterCache &paramCache,
const PrimaryVariables &sol,
const Element &element,
const FVElementGeometry &elemGeom,
int scvIdx,
const Problem &problem) const
{
return problem.boxTemperature(element, elemGeom, scvIdx);
Scalar T = problem.boxTemperature(element, elemGeom, scvIdx);
fs.setTemperature(T);
}
......@@ -79,39 +86,23 @@ public:
*
* Since we are isothermal, we don't need to do anything!
*/
template <class MutableParams>
void update(MutableParams &mutParams,
void update(FluidState &fs,
ParameterCache &paramCache,
const PrimaryVariables &sol,
const Element &element,
const FVElementGeometry &elemGeom,
int scvIdx,
const Problem &problem)
{
temperature_ = problem.boxTemperature(element, elemGeom, scvIdx);
// set the fluid temperatures
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
mutParams.setTemperature(phaseIdx, temperature_);
}
/*!
* \brief Return the temperature of any given phase [K]
*/
Scalar temperature(int phaseIdx = 0) const
{ return temperature_; }
/*!
* \brief If running under valgrind this produces an error message
* if some of the object's attributes is undefined.
*/
void checkDefined() const
{
Valgrind::CheckDefined(temperature_);
}
protected:
Scalar temperature_;
};
/*!
......@@ -130,7 +121,6 @@ class MPNCVolumeVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*kineticEnergyT
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MPNCIndices)) Indices;
enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)) };
......@@ -139,28 +129,33 @@ class MPNCVolumeVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*kineticEnergyT
enum { numEnergyEqs = Indices::NumPrimaryEnergyVars};
enum { temperature0Idx = Indices::temperatureIdx };
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename FluidSystem::ParameterCache ParameterCache;
typedef Dumux::EquilibriumFluidState<Scalar, FluidSystem> FluidState;
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 int dummy) const
void updateTemperatures(FluidState &fs,
ParameterCache &paramCache,
const PrimaryVariables &sol,
const Element &element,
const FVElementGeometry &elemGeom,
int scvIdx,
const Problem &problem) const
{
// retrieve temperature from solution vector
return sol[temperatureIdx];
Scalar T = sol[temperatureIdx];
fs.setTemperature(T);
}
/*!
* \brief Update the enthalpy and the internal energy for a given
* control volume.
*/
template <class MutableParams>
void update(MutableParams &mutParams,
const PrimaryVariables &sol,
void update(FluidState &fs,
ParameterCache &paramCache,
const Element &element,
const FVElementGeometry &elemGeom,
int scvIdx,
......@@ -168,24 +163,20 @@ public:
{
Valgrind::SetUndefined(*this);
// set the fluid temperatures
temperature_ = sol[temperature0Idx];
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
mutParams.setTemperature(phaseIdx, temperature_);
// heat capacities of the fluids plus the porous medium
heatCapacity_ =
problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
Valgrind::CheckDefined(heatCapacity_);
soilDensity_ =
problem.spatialParameters().soilDensity(element, elemGeom, scvIdx);
problem.spatialParameters().soilDensity(element, elemGeom, scvIdx);
Valgrind::CheckDefined(soilDensity_);
// set the enthalpies
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
Scalar h =
FluidSystem::computeEnthalpy(mutParams, phaseIdx);
mutParams.setEnthalpy(phaseIdx, h);
Scalar u =
FluidSystem::internalEnergy(fs, paramCache, phaseIdx);
fs.setInternalEnergy(phaseIdx, u);
}
}
......@@ -196,13 +187,6 @@ public:
Scalar heatCapacity() const
{ return heatCapacity_; };
/*!
* \brief Returns the temperature in fluid / solid phase(s)
* the sub-control volume.
*/
Scalar temperature(int phaseIdx = 0) const
{ return temperature_; }
/*!
* \brief Returns the total density of the given soil [kg / m^3] in
* the sub-control volume.
......@@ -218,13 +202,11 @@ public:
{
Valgrind::CheckDefined(heatCapacity_);
Valgrind::CheckDefined(soilDensity_);
Valgrind::CheckDefined(temperature_);
};
protected:
Scalar heatCapacity_;
Scalar soilDensity_;
Scalar temperature_;
};
} // end namepace
......
......@@ -60,7 +60,8 @@ class MPNCVolumeVariablesMass
enum { fug0Idx = Indices::fug0Idx };
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
typedef typename FluidSystem::ParameterCache ParameterCache;
public:
/*!
* \brief The fluid state which is used by the volume variables to
......@@ -75,11 +76,10 @@ public:
* \brief Update composition of all phases in the mutable
* parameters from the primary variables.
*/
template <class MutableParams>
void update(MutableParams &mutParams,
void update(FluidState &fs,
ParameterCache &paramCache,
const PrimaryVariables &priVars,
const VolumeVariables *hint,
const Problem &problem,
const Element &element,
const FVElementGeometry &elemGeom,
......@@ -97,10 +97,10 @@ public:
Scalar x_ij = 1.0/numComponents;
if (hint)
// use the hint for the initial mole fraction!
x_ij = hint->fluidState().moleFrac(phaseIdx, compIdx);
x_ij = hint->fluidState().moleFraction(phaseIdx, compIdx);
// set initial guess of the component's mole fraction
mutParams.setMoleFrac(phaseIdx,
fs.setMoleFraction(phaseIdx,
compIdx,
x_ij);
......@@ -111,18 +111,18 @@ public:
if (!hint)
// if no hint was given, we ask the constraint solver
// to make an initial guess
CompositionFromFugacitiesSolver::guessInitial(mutParams, phaseIdx, fug);
CompositionFromFugacitiesSolver::solve(mutParams, phaseIdx, fug);
CompositionFromFugacitiesSolver::guessInitial(fs, paramCache, phaseIdx, fug);
CompositionFromFugacitiesSolver::solve(fs, paramCache, phaseIdx, fug);
/*
std::cout << "final composition: " << FluidSystem::phaseName(phaseIdx) << "="
<< mutParams.moleFrac(phaseIdx, 0) << " "
<< mutParams.moleFrac(phaseIdx, 1) << " "
<< mutParams.moleFrac(phaseIdx, 2) << " "
<< mutParams.moleFrac(phaseIdx, 3) << " "
<< mutParams.moleFrac(phaseIdx, 4) << " "
<< mutParams.moleFrac(phaseIdx, 5) << " "
<< mutParams.moleFrac(phaseIdx, 6) << "\n";
<< fs.moleFrac(phaseIdx, 0) << " "
<< fs.moleFrac(phaseIdx, 1) << " "
<< fs.moleFrac(phaseIdx, 2) << " "
<< fs.moleFrac(phaseIdx, 3) << " "
<< fs.moleFrac(phaseIdx, 4) << " "
<< fs.moleFrac(phaseIdx, 5) << " "
<< fs.moleFrac(phaseIdx, 6) << "\n";
*/
}
......
......@@ -38,18 +38,22 @@ namespace Dumux {
template <class Scalar, class FluidSystem>
class CompositionFromFugacities
{
typedef typename FluidSystem::MutableParameters MutableParameters;
enum { numPhases = FluidSystem::numPhases };
enum { numComponents = FluidSystem::numComponents };
typedef typename FluidSystem::ParameterCache ParameterCache;
public:
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
/*!
* \brief Guess an initial value for the composition of the phase.
*/
static void guessInitial(MutableParameters &mutParams, int phaseIdx, const ComponentVector &fugVec)
template <class FluidState>
static void guessInitial(FluidState &fluidState,
ParameterCache &paramCache,
int phaseIdx,
const ComponentVector &fugVec)
{
if (FluidSystem::isIdealMixture(phaseIdx))
return;
......@@ -57,9 +61,9 @@ public:
// Pure component fugacities
for (int i = 0; i < numComponents; ++ i) {
//std::cout << f << " -> " << mutParams.fugacity(phaseIdx, i)/f << "\n";
mutParams.setMoleFrac(phaseIdx,
i,
1.0/numComponents);
fluidState.setMoleFraction(phaseIdx,
i,
1.0/numComponents);
}
};
......@@ -69,14 +73,16 @@ public:
*
* The phase's fugacities must already be set.
*/
static void solve(MutableParameters &mutParams,
template <class FluidState>
static void solve(FluidState &fluidState,
ParameterCache &paramCache,
int phaseIdx,
const ComponentVector &targetFug)
{
// use a much more efficient method in case the phase is an
// ideal mixture
if (FluidSystem::isIdealMixture(phaseIdx)) {
solveIdealMix_(mutParams, phaseIdx, targetFug);
solveIdealMix_(fluidState, paramCache, phaseIdx, targetFug);