Commit b027b376 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

Merge branch 'cleanup/remove-enthalpy-and-viscosity-from-constraintsolvers' into 'master'

Cleanup/remove enthalpy and viscosity from constraintsolvers

Closes #519

See merge request !1136
parents 49d40c4b e863f0f0
......@@ -71,10 +71,6 @@ namespace Dumux {
* - mean molar masses of *all* phases \f$M_\alpha\f$, \f$M_\beta\f$
* - fugacity coefficients of *all* components in *all* phases
* \f$\Phi^\kappa_\alpha\f$, \f$\Phi^\kappa_\beta\f$
* - if the setViscosity parameter is true, also dynamic viscosities of *all* phases
* \f$\mu_\alpha\f$, \f$\mu_\beta\f$
* - if the setEnthalpy parameter is true, also specific enthalpies and internal energies of *all* phases
* \f$h_\alpha\f$, \f$h_\beta\f$, \f$u_\alpha\f$, \f$u_\beta\f$
*/
template <class Scalar, class FluidSystem>
class ComputeFromReferencePhase
......@@ -107,24 +103,15 @@ public:
* - composition in mole and mass fractions and molaries of *all* phases
* - mean molar masses of *all* phases
* - fugacity coefficients of *all* components in *all* phases
* - if the setViscosity parameter is true, also dynamic viscosities of *all* phases
* - if the setEnthalpy parameter is true, also specific enthalpies and internal energies of *all* phases
*
* \param fluidState Thermodynamic state of the fluids
* \param paramCache Container for cache parameters
* \param refPhaseIdx The phase index of the reference phase
* \param setViscosity Specify whether the dynamic viscosity of
* each phase should also be set.
* \param setEnthalpy Specify whether the specific
* enthalpy/internal energy of each phase
* should also be set.
*/
template <class FluidState, class ParameterCache>
static void solve(FluidState &fluidState,
ParameterCache &paramCache,
int refPhaseIdx,
bool setViscosity,
bool setEnthalpy)
int refPhaseIdx)
{
ComponentVector fugVec;
......@@ -140,18 +127,6 @@ public:
paramCache,
refPhaseIdx));
if (setEnthalpy)
fluidState.setEnthalpy(refPhaseIdx,
FluidSystem::enthalpy(fluidState,
paramCache,
refPhaseIdx));
if (setViscosity)
fluidState.setViscosity(refPhaseIdx,
FluidSystem::viscosity(fluidState,
paramCache,
refPhaseIdx));
// compute the fugacities of all components in the reference phase
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
......@@ -172,18 +147,6 @@ public:
CompositionFromFugacities::guessInitial(fluidState, paramCache, phaseIdx, fugVec);
CompositionFromFugacities::solve(fluidState, paramCache, phaseIdx, fugVec);
if (setViscosity)
fluidState.setViscosity(phaseIdx,
FluidSystem::viscosity(fluidState,
paramCache,
phaseIdx));
if (setEnthalpy)
fluidState.setEnthalpy(phaseIdx,
FluidSystem::enthalpy(fluidState,
paramCache,
phaseIdx));
}
}
};
......
......@@ -54,8 +54,6 @@ namespace Dumux {
* - composition in mole and mass fractions and molarities of *all* phases
* - mean molar masses of *all* phases
* - fugacity coefficients of *all* components in *all* phases
* - if the setViscosity parameter is true, also dynamic viscosities of *all* phases
* - if the setEnthalpy parameter is true, also specific enthalpies of *all* phases
*/
template <class Scalar, class FluidSystem>
class MiscibleMultiPhaseComposition
......@@ -78,14 +76,10 @@ public:
*
* \param fluidState A container with the current (physical) state of the fluid
* \param paramCache A container for iterative calculation of fluid composition
* \param setViscosity Should the viscosity be set in the fluidstate?
* \param setEnthalpy Should the enthalpy be set in the fluidstate?
*/
template <class FluidState, class ParameterCache>
static void solve(FluidState &fluidState,
ParameterCache &paramCache,
bool setViscosity,
bool setEnthalpy,
int knownPhaseIdx = 0)
{
#ifndef NDEBUG
......@@ -205,16 +199,6 @@ public:
value = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
fluidState.setMolarDensity(phaseIdx, value);
if (setViscosity) {
value = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
fluidState.setViscosity(phaseIdx, value);
}
if (setEnthalpy) {
value = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
fluidState.setEnthalpy(phaseIdx, value);
}
}
}
};
......
......@@ -250,9 +250,7 @@ public:
// of the "MiscibleMultiPhaseComposition" constraint solver
if(useConstraintSolver)
MiscibleMultiPhaseComposition::solve(fluidState,
paramCache,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
paramCache);
// ... or calculated explicitly this way ...
else
{
......@@ -299,9 +297,7 @@ public:
if (useConstraintSolver)
ComputeFromReferencePhase::solve(fluidState,
paramCache,
phase1Idx,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
phase1Idx);
// ... or calculated explicitly this way ...
else
{
......@@ -349,9 +345,7 @@ public:
if (useConstraintSolver)
ComputeFromReferencePhase::solve(fluidState,
paramCache,
phase0Idx,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
phase0Idx);
// ... or calculated explicitly this way ...
else
{
......@@ -387,11 +381,11 @@ public:
fluidState.setDensity(phaseIdx, rho);
Scalar rhoMolar = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
fluidState.setMolarDensity(phaseIdx, rhoMolar);
const Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
fluidState.setViscosity(phaseIdx,mu);
}
// compute and set the enthalpy
const Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
fluidState.setViscosity(phaseIdx,mu);
Scalar h = EnergyVolVars::enthalpy(fluidState, paramCache, phaseIdx);
fluidState.setEnthalpy(phaseIdx, h);
}
......
......@@ -262,8 +262,6 @@ public:
MiscibleMultiPhaseComposition::solve(fluidState,
paramCache,
/*setViscosity=*/true,
/*setEnthalpy=*/false,
knownPhaseIdx);
}
else if (phasePresence == secondPhaseOnly)
......@@ -291,9 +289,7 @@ public:
// of the "ComputeFromReferencePhase" constraint solver
ComputeFromReferencePhase::solve(fluidState,
paramCache,
phase1Idx,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
phase1Idx);
}
else if (phasePresence == firstPhaseOnly)
{
......@@ -321,9 +317,7 @@ public:
// of the "ComputeFromReferencePhase" constraint solver
ComputeFromReferencePhase::solve(fluidState,
paramCache,
phase0Idx,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
phase0Idx);
}
paramCache.updateAll(fluidState);
for (int phaseIdx = 0; phaseIdx < ModelTraits::numPhases(); ++phaseIdx)
......
......@@ -209,9 +209,7 @@ public:
// constraint solver ...
if (useConstraintSolver) {
MiscibleMultiPhaseComposition::solve(fluidState_,
paramCache,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
paramCache);
}
// ... or calculated explicitly this way ...
// please note that we experienced some problems with un-regularized
......@@ -293,9 +291,7 @@ public:
{
ComputeFromReferencePhase::solve(fluidState_,
paramCache,
wPhaseIdx,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
wPhaseIdx);
}
// ... or calculated explicitly this way ...
else {
......@@ -368,9 +364,7 @@ public:
// of the "ComputeFromReferencePhase" constraint solver
ComputeFromReferencePhase::solve(fluidState_,
paramCache,
gPhaseIdx,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
gPhaseIdx);
}
else if (phasePresence == wnPhaseOnly) {
// only water and NAPL phases are present
......@@ -392,9 +386,7 @@ public:
// of the "ComputeFromReferencePhase" constraint solver
ComputeFromReferencePhase::solve(fluidState_,
paramCache,
wPhaseIdx,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
wPhaseIdx);
}
else if (phasePresence == gPhaseOnly) {
// only the gas phase is present, gas phase composition is
......@@ -416,9 +408,7 @@ public:
{
ComputeFromReferencePhase::solve(fluidState_,
paramCache,
gPhaseIdx,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
gPhaseIdx);
}
// ... or calculated explicitly this way ...
else {
......@@ -484,9 +474,7 @@ public:
{
ComputeFromReferencePhase::solve(fluidState_,
paramCache,
gPhaseIdx,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
gPhaseIdx);
}
// ... or calculated explicitly this way ...
else {
......
......@@ -34,10 +34,11 @@
#include <dumux/material/constants.hh>
#include <dumux/material/fluidstates/compositional.hh>
#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
#include <dumux/material/solidstates/updatesolidvolumefractions.hh>
#include <dumux/common/valgrind.hh>
#include <dumux/common/exceptions.hh>
namespace Dumux {
/*!
......
......@@ -713,41 +713,39 @@ public:
// // For using the ... other way of calculating equilibrium
// THIS IS ONLY FOR silencing Valgrind but is not used in this model
for(int phaseIdx=0; phaseIdx<numPhases(); ++phaseIdx)
for (int compIdx = 0; compIdx < numFluidComps; ++compIdx) {
const Scalar phi = FluidSystem::fugacityCoefficient(actualFluidState,
for(int phaseIdx=0; phaseIdx<numPhases(); ++phaseIdx)
for (int compIdx = 0; compIdx < numFluidComps; ++compIdx) {
const Scalar phi = FluidSystem::fugacityCoefficient(actualFluidState,
paramCache,
phaseIdx,
compIdx);
actualFluidState.setFugacityCoefficient(phaseIdx,
actualFluidState.setFugacityCoefficient(phaseIdx,
compIdx,
phi);
}
FluidState equilFluidState; // the fluidState *on the interface* i.e. chemical equilibrium
equilFluidState.assign(actualFluidState) ;
ConstraintSolver::solve(equilFluidState,
paramCache,
/*setViscosity=*/false,
/*setEnthalpy=*/false) ;
// Setting the equilibrium composition (in a kinetic model not necessarily the same as the actual mole fraction)
for(int phaseIdx=0; phaseIdx<numPhases(); ++phaseIdx){
for (int compIdx=0; compIdx< numFluidComps; ++ compIdx){
xEquil_[phaseIdx][compIdx] = equilFluidState.moleFraction(phaseIdx, compIdx);
}
}
FluidState equilFluidState; // the fluidState *on the interface* i.e. chemical equilibrium
equilFluidState.assign(actualFluidState) ;
ConstraintSolver::solve(equilFluidState,
paramCache) ;
// compute densities of all phases
for(int phaseIdx=0; phaseIdx<numPhases(); ++phaseIdx){
const Scalar rho = FluidSystem::density(actualFluidState, paramCache, phaseIdx);
actualFluidState.setDensity(phaseIdx, rho);
const Scalar rhoMolar = FluidSystem::molarDensity(actualFluidState, paramCache, phaseIdx);
actualFluidState.setMolarDensity(phaseIdx, rhoMolar);
// Setting the equilibrium composition (in a kinetic model not necessarily the same as the actual mole fraction)
for(int phaseIdx=0; phaseIdx<numPhases(); ++phaseIdx){
for (int compIdx=0; compIdx< numFluidComps; ++ compIdx){
xEquil_[phaseIdx][compIdx] = equilFluidState.moleFraction(phaseIdx, compIdx);
}
}
// compute densities of all phases
for(int phaseIdx=0; phaseIdx<numPhases(); ++phaseIdx){
const Scalar rho = FluidSystem::density(actualFluidState, paramCache, phaseIdx);
actualFluidState.setDensity(phaseIdx, rho);
const Scalar rhoMolar = FluidSystem::molarDensity(actualFluidState, paramCache, phaseIdx);
actualFluidState.setMolarDensity(phaseIdx, rhoMolar);
}
}
/*!
* \brief The mole fraction we would have in the case of chemical equilibrium /
* on the interface.
......
......@@ -143,9 +143,7 @@ void completeReferenceFluidState(FluidState &fs,
typename FluidSystem::ParameterCache paramCache;
ComputeFromReferencePhase::solve(fs,
paramCache,
refPhaseIdx,
/*setViscosity=*/false,
/*setEnthalpy=*/false);
refPhaseIdx);
}
......@@ -251,9 +249,7 @@ int main()
FluidSystem::ParameterCache paramCache;
using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
MiscibleMultiPhaseComposition::solve(fsRef, paramCache,
/*setViscosity=*/false,
/*setEnthalpy=*/false);
MiscibleMultiPhaseComposition::solve(fsRef, paramCache);
// check the flash calculation
checkNcpFlash<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams);
......@@ -284,9 +280,7 @@ int main()
+ (pc[gasPhaseIdx] - pc[liquidPhaseIdx]));
using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
MiscibleMultiPhaseComposition::solve(fsRef, paramCache,
/*setViscosity=*/false,
/*setEnthalpy=*/false);
MiscibleMultiPhaseComposition::solve(fsRef, paramCache);
// check the flash calculation
......
......@@ -275,10 +275,8 @@ private:
using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
ParameterCache paramCache;
MiscibleMultiPhaseComposition::solve(fs,
paramCache,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
MiscibleMultiPhaseComposition::solve(fs, paramCache);
///////////
// assign the primary variables
///////////
......
......@@ -465,9 +465,7 @@ private:
ParameterCache paramCache;
ComputeFromReferencePhase::solve(fluidState,
paramCache,
refPhaseIdx,
/*setViscosity=*/false,
/*setEnthalpy=*/false);
refPhaseIdx);
//////////////////////////////////////
// Set fugacities
......
......@@ -323,9 +323,7 @@ public:
// This solves the system of equations defining x=x(p,T)
ConstraintSolver::solve(fluidState,
dummyCache,
/*setViscosity=*/false,
/*setEnthalpy=*/false) ;
dummyCache) ;
// Now let's make the air phase less than fully saturated with water
fluidState.setMoleFraction(gasPhaseIdx, wCompIdx, fluidState.moleFraction(gasPhaseIdx, wCompIdx)*percentOfEquil_ ) ;
......@@ -480,9 +478,7 @@ private:
// This solves the system of equations defining x=x(p,T)
ParameterCache dummyCache;
ConstraintSolver::solve(equilibriumFluidState,
dummyCache,
/*setViscosity=*/false,
/*setEnthalpy=*/false) ;
dummyCache) ;
FluidState dryFluidState(equilibriumFluidState);
// Now let's make the air phase less than fully saturated with vapor
......
......@@ -358,9 +358,7 @@ private:
ParameterCache paramCache;
ComputeFromReferencePhase::solve(fs,
paramCache,
refPhaseIdx,
/*setViscosity=*/false,
/*setEnthalpy=*/false);
refPhaseIdx);
///////////
// assign the primary variables
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment