mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
cosmetic cleanups
- remove stray AC_MSG_RESULT from DUMUX_CHECK_QUAD - slight reformulations of the usage message of test_2p - replace tabulators with four spaces - remove trailing whitespaceformating
This commit is contained in:
committed by
Andreas Lauser
parent
8732eed885
commit
388080231e
@@ -78,7 +78,7 @@ public:
|
||||
const Scalar pg0=1.e5; /* reference pressure */
|
||||
const Scalar T0=273.15; /* reference temperature */
|
||||
Scalar Dgaw;
|
||||
|
||||
|
||||
Dgaw=Daw*(pg0/pressure)*pow((temperature/T0),Theta);
|
||||
|
||||
return Dgaw;
|
||||
|
||||
@@ -6,8 +6,8 @@
|
||||
* Technical University of Berlin *
|
||||
* email:lstadler@wahyd.tu-berlin.de *
|
||||
* *
|
||||
* Copyright (C) 2011 by Benjamin Faigle *
|
||||
* Copyright (C) 2011 by Holger Class *
|
||||
* Copyright (C) 2011 by Benjamin Faigle *
|
||||
* Copyright (C) 2011 by Holger Class *
|
||||
* Institute for Modelling Hydraulic and Environmental Systems *
|
||||
* University of Stuttgart, Germany *
|
||||
* email: <givenname>.<name>@iws.uni-stuttgart.de *
|
||||
@@ -79,7 +79,7 @@ public:
|
||||
*/
|
||||
static Scalar criticalPressure()
|
||||
{ return 37.86e5; /* [Pa] */ }
|
||||
|
||||
|
||||
/*!
|
||||
* \brief The density of \f$AIR\f$ at a given pressure and temperature [kg/m^3].
|
||||
*
|
||||
@@ -206,7 +206,7 @@ public:
|
||||
Scalar pressure)
|
||||
{
|
||||
return
|
||||
gasEnthalpy(temperature, pressure)
|
||||
gasEnthalpy(temperature, pressure)
|
||||
-
|
||||
IdealGas::R * temperature // = pressure * molar volume for an ideal gas
|
||||
/ molarMass(); // conversion from [J/(mol K)] to [J/(kg K)]
|
||||
|
||||
@@ -233,7 +233,7 @@ public:
|
||||
Scalar pressure)
|
||||
{
|
||||
return
|
||||
gasEnthalpy(temperature, pressure) -
|
||||
gasEnthalpy(temperature, pressure) -
|
||||
pressure/gasDensity(temperature, pressure);
|
||||
}
|
||||
|
||||
|
||||
@@ -46,7 +46,7 @@ class Component
|
||||
{
|
||||
public:
|
||||
static const bool isTabulated = false;
|
||||
|
||||
|
||||
/*!
|
||||
* \brief A default routine for initialization, not needed for components and must not be called.
|
||||
*
|
||||
|
||||
@@ -289,7 +289,7 @@ public:
|
||||
result +=
|
||||
n(i) *
|
||||
std::pow(7.1 - pi_, I(i)) *
|
||||
J(i) *
|
||||
J(i) *
|
||||
(J(i) - 1) *
|
||||
std::pow(tau_ - 1.222, J(i) - 2);
|
||||
}
|
||||
|
||||
@@ -182,7 +182,7 @@ public:
|
||||
|
||||
Scalar Fp0 = 1.0;
|
||||
Scalar xi = 0.00474;
|
||||
Scalar eta_xi =
|
||||
Scalar eta_xi =
|
||||
Fp0*(0.807*std::pow(Tr,0.618)
|
||||
- 0.357*std::exp(-0.449*Tr)
|
||||
+ 0.34*std::exp(-4.058*Tr)
|
||||
@@ -277,7 +277,7 @@ protected:
|
||||
}
|
||||
else {
|
||||
assert(temperature>=348.0);
|
||||
|
||||
|
||||
/* take care: extrapolation for Temperature>373 */
|
||||
H = 16.7+2.1*(temperature-348.0)/25.; /* leads probably to underestimates */
|
||||
CH3 = 45.8+2.5*(temperature-348.0)/25.;
|
||||
|
||||
@@ -142,7 +142,7 @@ public:
|
||||
*/
|
||||
static bool gasIsCompressible()
|
||||
{ return true; }
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Returns true iff the gas phase is assumed to be ideal
|
||||
*/
|
||||
@@ -216,7 +216,7 @@ public:
|
||||
return
|
||||
gasEnthalpy(temperature, pressure) -
|
||||
1/molarMass()* // conversion from [J/(mol K)] to [J/(kg K)]
|
||||
IdealGas::R*temperature; // = pressure * spec. volume for an ideal gas
|
||||
IdealGas::R*temperature; // = pressure * spec. volume for an ideal gas
|
||||
}
|
||||
|
||||
/*!
|
||||
|
||||
@@ -122,13 +122,13 @@ public:
|
||||
// reference phase
|
||||
paramCache.updatePhase(fluidState, refPhaseIdx);
|
||||
fluidState.setDensity(refPhaseIdx,
|
||||
FluidSystem::density(fluidState,
|
||||
FluidSystem::density(fluidState,
|
||||
paramCache,
|
||||
refPhaseIdx));
|
||||
|
||||
if (setEnthalpy)
|
||||
fluidState.setEnthalpy(refPhaseIdx,
|
||||
FluidSystem::enthalpy(fluidState,
|
||||
FluidSystem::enthalpy(fluidState,
|
||||
paramCache,
|
||||
refPhaseIdx));
|
||||
|
||||
@@ -140,7 +140,7 @@ public:
|
||||
|
||||
// compute the fugacities of all components in the reference phase
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
fluidState.setFugacityCoefficient(refPhaseIdx,
|
||||
fluidState.setFugacityCoefficient(refPhaseIdx,
|
||||
compIdx,
|
||||
FluidSystem::fugacityCoefficient(fluidState,
|
||||
paramCache,
|
||||
@@ -159,10 +159,10 @@ public:
|
||||
|
||||
if (setViscosity)
|
||||
fluidState.setViscosity(phaseIdx,
|
||||
FluidSystem::viscosity(fluidState,
|
||||
paramCache,
|
||||
FluidSystem::viscosity(fluidState,
|
||||
paramCache,
|
||||
phaseIdx));
|
||||
|
||||
|
||||
if (setEnthalpy)
|
||||
fluidState.setEnthalpy(phaseIdx,
|
||||
FluidSystem::enthalpy(fluidState,
|
||||
|
||||
@@ -67,7 +67,7 @@ class ImmiscibleFlash
|
||||
{
|
||||
static constexpr int numPhases = FluidSystem::numPhases;
|
||||
static constexpr int numComponents = FluidSystem::numComponents;
|
||||
static_assert(numPhases == numComponents,
|
||||
static_assert(numPhases == numComponents,
|
||||
"Immiscibility assumes that the number of phases"
|
||||
" is equal to the number of components");
|
||||
|
||||
@@ -94,7 +94,7 @@ public:
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
|
||||
sumMoles += globalMolarities[compIdx];
|
||||
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
|
||||
// pressure. use atmospheric pressure as initial guess
|
||||
fluidState.setPressure(phaseIdx, 2e5);
|
||||
|
||||
@@ -127,14 +127,14 @@ public:
|
||||
break;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
if (allIncompressible) {
|
||||
paramCache.updateAll(fluidState);
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
|
||||
fluidState.setDensity(phaseIdx, rho);
|
||||
|
||||
Scalar saturation =
|
||||
|
||||
Scalar saturation =
|
||||
globalMolarities[/*compIdx=*/phaseIdx]
|
||||
/ fluidState.molarDensity(phaseIdx);
|
||||
fluidState.setSaturation(phaseIdx, saturation);
|
||||
@@ -144,7 +144,7 @@ public:
|
||||
/////////////////////////
|
||||
// Newton method
|
||||
/////////////////////////
|
||||
|
||||
|
||||
// Jacobian matrix
|
||||
Matrix J;
|
||||
// solution, i.e. phase composition
|
||||
@@ -154,7 +154,7 @@ public:
|
||||
|
||||
Valgrind::SetUndefined(J);
|
||||
Valgrind::SetUndefined(deltaX);
|
||||
Valgrind::SetUndefined(b);
|
||||
Valgrind::SetUndefined(b);
|
||||
|
||||
completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
|
||||
|
||||
@@ -164,18 +164,18 @@ public:
|
||||
linearize_<MaterialLaw>(J, b, fluidState, paramCache, matParams, globalMolarities);
|
||||
Valgrind::CheckDefined(J);
|
||||
Valgrind::CheckDefined(b);
|
||||
|
||||
|
||||
// Solve J*x = b
|
||||
deltaX = 0;
|
||||
|
||||
try { J.solve(deltaX, b); }
|
||||
catch (Dune::FMatrixError e)
|
||||
{
|
||||
{
|
||||
/*
|
||||
std::cout << "error: " << e << "\n";
|
||||
std::cout << "b: " << b << "\n";
|
||||
*/
|
||||
|
||||
|
||||
throw Dumux::NumericalProblem(e.what());
|
||||
}
|
||||
Valgrind::CheckDefined(deltaX);
|
||||
@@ -187,7 +187,7 @@ public:
|
||||
for (int j = 0; j < numEq; ++j) {
|
||||
std::ostringstream os;
|
||||
os << J[i][j];
|
||||
|
||||
|
||||
std::string s(os.str());
|
||||
do {
|
||||
s += " ";
|
||||
@@ -203,7 +203,7 @@ public:
|
||||
|
||||
// update the fluid quantities.
|
||||
Scalar relError = update_<MaterialLaw>(fluidState, paramCache, matParams, deltaX);
|
||||
|
||||
|
||||
if (relError < 1e-9)
|
||||
return;
|
||||
}
|
||||
@@ -218,7 +218,7 @@ public:
|
||||
|
||||
DUNE_THROW(NumericalProblem,
|
||||
"Flash calculation failed."
|
||||
" {c_alpha^kappa} = {" << globalMolarities << "}, T = "
|
||||
" {c_alpha^kappa} = {" << globalMolarities << "}, T = "
|
||||
<< fluidState.temperature(/*phaseIdx=*/0));
|
||||
}
|
||||
|
||||
@@ -231,7 +231,7 @@ protected:
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
|
||||
std::cout << fs.saturation(phaseIdx) << " ";
|
||||
std::cout << "\n";
|
||||
|
||||
|
||||
std::cout << "pressures: ";
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
|
||||
std::cout << fs.pressure(phaseIdx) << " ";
|
||||
@@ -263,7 +263,7 @@ protected:
|
||||
{
|
||||
FluidState origFluidState(fluidState);
|
||||
ParameterCache origParamCache(paramCache);
|
||||
|
||||
|
||||
Vector tmp;
|
||||
|
||||
// reset jacobian
|
||||
@@ -286,7 +286,7 @@ protected:
|
||||
const Scalar eps = 1e-10/quantityWeight_(fluidState, pvIdx);
|
||||
setQuantity_<MaterialLaw>(fluidState, paramCache, matParams, pvIdx, x_i + eps);
|
||||
assert(getQuantity_(fluidState, pvIdx) == x_i + eps);
|
||||
|
||||
|
||||
// compute derivative of the defect
|
||||
calculateDefect_(tmp, origFluidState, fluidState, globalMolarities);
|
||||
tmp -= b;
|
||||
@@ -312,10 +312,10 @@ protected:
|
||||
{
|
||||
// global molarities are given
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
b[phaseIdx] =
|
||||
b[phaseIdx] =
|
||||
fluidState.saturation(phaseIdx)
|
||||
* fluidState.molarity(phaseIdx, /*compIdx=*/phaseIdx);
|
||||
|
||||
|
||||
b[phaseIdx] -= globalMolarities[/*compIdx=*/phaseIdx];
|
||||
}
|
||||
}
|
||||
@@ -330,9 +330,9 @@ protected:
|
||||
for (int pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
|
||||
Scalar tmp = getQuantity_(fluidState, pvIdx);
|
||||
Scalar delta = deltaX[pvIdx];
|
||||
|
||||
|
||||
relError = std::max(relError, std::abs(delta)*quantityWeight_(fluidState, pvIdx));
|
||||
|
||||
|
||||
if (isSaturationIdx_(pvIdx)) {
|
||||
// dampen to at most 20% change in saturation per
|
||||
// iteration
|
||||
@@ -343,10 +343,10 @@ protected:
|
||||
// iteration
|
||||
delta = std::min(0.30*fluidState.pressure(0), std::max(-0.30*fluidState.pressure(0), delta));
|
||||
};
|
||||
|
||||
|
||||
setQuantityRaw_(fluidState, pvIdx, tmp - delta);
|
||||
}
|
||||
|
||||
|
||||
completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
|
||||
|
||||
return relError;
|
||||
@@ -376,14 +376,14 @@ protected:
|
||||
|
||||
// set the last saturation
|
||||
fluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);
|
||||
|
||||
|
||||
// update the pressures using the material law (saturations
|
||||
// and first pressure are already set because it is implicitly
|
||||
// solved for.)
|
||||
ComponentVector pC;
|
||||
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
|
||||
for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
|
||||
fluidState.setPressure(phaseIdx,
|
||||
fluidState.setPressure(phaseIdx,
|
||||
fluidState.pressure(0)
|
||||
+ (pC[phaseIdx] - pC[0]));
|
||||
|
||||
@@ -415,7 +415,7 @@ protected:
|
||||
return fs.pressure(phaseIdx);
|
||||
}
|
||||
// saturations
|
||||
else {
|
||||
else {
|
||||
assert(pvIdx < numPhases);
|
||||
int phaseIdx = pvIdx - 1;
|
||||
return fs.saturation(phaseIdx);
|
||||
@@ -457,7 +457,7 @@ protected:
|
||||
|
||||
// set last saturation (-> minus the change of the
|
||||
// saturation of the other phases)
|
||||
fs.setSaturation(/*phaseIdx=*/numPhases - 1,
|
||||
fs.setSaturation(/*phaseIdx=*/numPhases - 1,
|
||||
fs.saturation(numPhases - 1) - delta);
|
||||
|
||||
// update all fluid pressures using the capillary pressure
|
||||
@@ -465,7 +465,7 @@ protected:
|
||||
ComponentVector pC;
|
||||
MaterialLaw::capillaryPressures(pC, matParams, fs);
|
||||
for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
|
||||
fs.setPressure(phaseIdx,
|
||||
fs.setPressure(phaseIdx,
|
||||
fs.pressure(0)
|
||||
+ (pC[phaseIdx] - pC[0]));
|
||||
paramCache.updateAllPressures(fs);
|
||||
@@ -483,7 +483,7 @@ protected:
|
||||
static void setQuantityRaw_(FluidState &fs, int pvIdx, Scalar value)
|
||||
{
|
||||
assert(pvIdx < numEq);
|
||||
|
||||
|
||||
// first pressure
|
||||
if (pvIdx < 1) {
|
||||
int phaseIdx = 0;
|
||||
|
||||
@@ -46,7 +46,7 @@ namespace Dumux {
|
||||
* - temperatures of *all* phases
|
||||
* - saturations of *all* phases
|
||||
* - pressures of *all* phases
|
||||
*
|
||||
*
|
||||
* It also assumes that the mole/mass fractions of all phases sum up
|
||||
* to 1. After calling the solve() method the following quantities
|
||||
* are calculated in addition:
|
||||
@@ -64,7 +64,7 @@ class MiscibleMultiPhaseComposition
|
||||
{
|
||||
static constexpr int numPhases = FluidSystem::numPhases;
|
||||
static constexpr int numComponents = FluidSystem::numComponents;
|
||||
|
||||
|
||||
static_assert(numComponents == numPhases,
|
||||
"This solver requires that the number fluid phases is equal "
|
||||
"to the number of components");
|
||||
@@ -81,7 +81,7 @@ public:
|
||||
* - temperatures of *all* phases
|
||||
* - saturations of *all* phases
|
||||
* - pressures of *all* phases
|
||||
*
|
||||
*
|
||||
* It also assumes that the mole/mass fractions of all phases sum up
|
||||
* to 1. After calling the solve() method the following quantities
|
||||
* are calculated in addition:
|
||||
@@ -122,14 +122,14 @@ public:
|
||||
fluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
// create the linear system of equations which defines the
|
||||
// mole fractions
|
||||
Dune::FieldMatrix<Scalar, numComponents*numPhases, numComponents*numPhases> M(0.0);
|
||||
Dune::FieldVector<Scalar, numComponents*numPhases> x(0.0);
|
||||
Dune::FieldVector<Scalar, numComponents*numPhases> b(0.0);
|
||||
|
||||
|
||||
// assemble the equations expressing the assumption that the
|
||||
// sum of all mole fractions in each phase must be 1
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
@@ -142,11 +142,11 @@ public:
|
||||
M[rowIdx][colIdx] = 1.0;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// assemble the equations expressing the fact that the
|
||||
// fugacities of each component is equal in all phases
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
Scalar entryCol1 =
|
||||
Scalar entryCol1 =
|
||||
fluidState.fugacityCoefficient(/*phaseIdx=*/0, compIdx)
|
||||
* fluidState.pressure(/*phaseIdx=*/0);
|
||||
int col1Idx = compIdx;
|
||||
@@ -155,7 +155,7 @@ public:
|
||||
int rowIdx = (phaseIdx - 1)*numComponents + compIdx;
|
||||
int col2Idx = phaseIdx*numComponents + compIdx;
|
||||
|
||||
Scalar entryCol2 =
|
||||
Scalar entryCol2 =
|
||||
fluidState.fugacityCoefficient(phaseIdx, compIdx)
|
||||
* fluidState.pressure(phaseIdx);
|
||||
|
||||
@@ -163,7 +163,7 @@ public:
|
||||
M[rowIdx][col2Idx] = -entryCol2;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// solve for all mole fractions
|
||||
M.solve(x, b);
|
||||
|
||||
@@ -175,7 +175,7 @@ public:
|
||||
fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]);
|
||||
}
|
||||
paramCache.updateComposition(fluidState, phaseIdx);
|
||||
|
||||
|
||||
Scalar value = FluidSystem::density(fluidState, paramCache, phaseIdx);
|
||||
fluidState.setDensity(phaseIdx, value);
|
||||
|
||||
@@ -183,7 +183,7 @@ public:
|
||||
value = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
|
||||
fluidState.setViscosity(phaseIdx, value);
|
||||
}
|
||||
|
||||
|
||||
if (setInternalEnergy) {
|
||||
value = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
|
||||
fluidState.setEnthalpy(phaseIdx, value);
|
||||
|
||||
@@ -110,14 +110,14 @@ public:
|
||||
fluidState.setMoleFraction(phaseIdx,
|
||||
compIdx,
|
||||
globalMolarities[compIdx]/sumMoles);
|
||||
|
||||
|
||||
// pressure. use atmospheric pressure as initial guess
|
||||
fluidState.setPressure(phaseIdx, 2e5);
|
||||
|
||||
// saturation. assume all fluids to be equally distributed
|
||||
fluidState.setSaturation(phaseIdx, 1.0/numPhases);
|
||||
}
|
||||
|
||||
|
||||
// set the fugacity coefficients of all components in all phases
|
||||
paramCache.updateAll(fluidState);
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
|
||||
@@ -156,7 +156,7 @@ public:
|
||||
Valgrind::SetUndefined(J);
|
||||
Valgrind::SetUndefined(deltaX);
|
||||
Valgrind::SetUndefined(b);
|
||||
|
||||
|
||||
// make the fluid state consistent with the fluid system.
|
||||
completeFluidState_<MaterialLaw>(fluidState,
|
||||
paramCache,
|
||||
@@ -176,13 +176,13 @@ public:
|
||||
linearize_<MaterialLaw>(J, b, fluidState, paramCache, matParams, globalMolarities);
|
||||
Valgrind::CheckDefined(J);
|
||||
Valgrind::CheckDefined(b);
|
||||
|
||||
|
||||
// Solve J*x = b
|
||||
deltaX = 0;
|
||||
|
||||
try { J.solve(deltaX, b); }
|
||||
catch (Dune::FMatrixError e)
|
||||
{
|
||||
{
|
||||
/*
|
||||
printFluidState_(fluidState);
|
||||
std::cout << "error: " << e << "\n";
|
||||
@@ -201,7 +201,7 @@ public:
|
||||
for (int j = 0; j < numEq; ++j) {
|
||||
std::ostringstream os;
|
||||
os << J[i][j];
|
||||
|
||||
|
||||
std::string s(os.str());
|
||||
do {
|
||||
s += " ";
|
||||
@@ -217,7 +217,7 @@ public:
|
||||
|
||||
// update the fluid quantities.
|
||||
Scalar relError = update_<MaterialLaw>(fluidState, paramCache, matParams, deltaX);
|
||||
|
||||
|
||||
if (relError < 1e-9)
|
||||
return;
|
||||
}
|
||||
@@ -232,7 +232,7 @@ public:
|
||||
|
||||
DUNE_THROW(NumericalProblem,
|
||||
"Flash calculation failed."
|
||||
" {c_alpha^kappa} = {" << globalMolarities << "}, T = "
|
||||
" {c_alpha^kappa} = {" << globalMolarities << "}, T = "
|
||||
<< fluidState.temperature(/*phaseIdx=*/0));
|
||||
}
|
||||
|
||||
@@ -245,7 +245,7 @@ protected:
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
|
||||
std::cout << fs.saturation(phaseIdx) << " ";
|
||||
std::cout << "\n";
|
||||
|
||||
|
||||
std::cout << "pressures: ";
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
|
||||
std::cout << fs.pressure(phaseIdx) << " ";
|
||||
@@ -293,7 +293,7 @@ protected:
|
||||
{
|
||||
FluidState origFluidState(fluidState);
|
||||
ParameterCache origParamCache(paramCache);
|
||||
|
||||
|
||||
Vector tmp;
|
||||
|
||||
// reset jacobian
|
||||
@@ -316,7 +316,7 @@ protected:
|
||||
const Scalar eps = 1e-8/quantityWeight_(fluidState, pvIdx);
|
||||
setQuantity_<MaterialLaw>(fluidState, paramCache, matParams, pvIdx, x_i + eps);
|
||||
assert(getQuantity_(fluidState, pvIdx) == x_i + eps);
|
||||
|
||||
|
||||
// compute derivative of the defect
|
||||
calculateDefect_(tmp, origFluidState, fluidState, globalMolarities);
|
||||
tmp -= b;
|
||||
@@ -346,32 +346,32 @@ protected:
|
||||
// fugacity of any component must be equal in all phases
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
|
||||
b[eqIdx] =
|
||||
fluidState.fugacity(/*phaseIdx=*/0, compIdx) -
|
||||
b[eqIdx] =
|
||||
fluidState.fugacity(/*phaseIdx=*/0, compIdx) -
|
||||
fluidState.fugacity(phaseIdx, compIdx);
|
||||
++eqIdx;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
assert(eqIdx == numComponents*(numPhases - 1));
|
||||
|
||||
|
||||
// the fact saturations must sum up to 1 is included explicitly!
|
||||
|
||||
|
||||
// capillary pressures are explicitly included!
|
||||
|
||||
|
||||
// global molarities are given
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
b[eqIdx] = 0.0;
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
b[eqIdx] +=
|
||||
b[eqIdx] +=
|
||||
fluidState.saturation(phaseIdx)
|
||||
* fluidState.molarity(phaseIdx, compIdx);
|
||||
}
|
||||
|
||||
|
||||
b[eqIdx] -= globalMolarities[compIdx];
|
||||
++eqIdx;
|
||||
}
|
||||
|
||||
|
||||
// model assumptions (-> non-linear complementarity functions)
|
||||
// must be adhered
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
@@ -403,9 +403,9 @@ protected:
|
||||
for (int pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
|
||||
Scalar tmp = getQuantity_(fluidState, pvIdx);
|
||||
Scalar delta = deltaX[pvIdx];
|
||||
|
||||
|
||||
relError = std::max(relError, std::abs(delta)*quantityWeight_(fluidState, pvIdx));
|
||||
|
||||
|
||||
if (isSaturationIdx_(pvIdx)) {
|
||||
// dampen to at most 20% change in saturation per
|
||||
// iteration
|
||||
@@ -421,7 +421,7 @@ protected:
|
||||
// iteration
|
||||
delta = std::min(0.15*fluidState.pressure(0), std::max(-0.15*fluidState.pressure(0), delta));
|
||||
};
|
||||
|
||||
|
||||
setQuantityRaw_(fluidState, pvIdx, tmp - delta);
|
||||
}
|
||||
|
||||
@@ -451,9 +451,9 @@ protected:
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
Scalar value = fluidState.saturation(phaseIdx)/(0.95*sumSat);
|
||||
fluidState.setSaturation(phaseIdx, value);
|
||||
};
|
||||
};
|
||||
};
|
||||
|
||||
|
||||
completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
|
||||
|
||||
return relError;
|
||||
@@ -470,14 +470,14 @@ protected:
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
|
||||
sumSat += fluidState.saturation(phaseIdx);
|
||||
fluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);
|
||||
|
||||
|
||||
// update the pressures using the material law (saturations
|
||||
// and first pressure are already set because it is implicitly
|
||||
// solved for.)
|
||||
ComponentVector pC;
|
||||
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
|
||||
for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
|
||||
fluidState.setPressure(phaseIdx,
|
||||
fluidState.setPressure(phaseIdx,
|
||||
fluidState.pressure(0)
|
||||
+ (pC[phaseIdx] - pC[0]));
|
||||
|
||||
@@ -488,7 +488,7 @@ protected:
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
|
||||
Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
|
||||
fluidState.setDensity(phaseIdx, rho);
|
||||
|
||||
|
||||
for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
|
||||
Scalar phi = FluidSystem::fugacityCoefficient( fluidState, paramCache, phaseIdx, compIdx);
|
||||
fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
|
||||
@@ -522,7 +522,7 @@ protected:
|
||||
return fs.saturation(phaseIdx);
|
||||
}
|
||||
// mole fractions
|
||||
else // if (pvIdx < numPhases + numPhases*numComponents)
|
||||
else // if (pvIdx < numPhases + numPhases*numComponents)
|
||||
{
|
||||
int phaseIdx = (pvIdx - numPhases)/numComponents;
|
||||
int compIdx = (pvIdx - numPhases)%numComponents;
|
||||
@@ -553,7 +553,7 @@ protected:
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
|
||||
fs.setDensity(phaseIdx, rho);
|
||||
|
||||
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
Scalar phi = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx);
|
||||
fs.setFugacityCoefficient(phaseIdx, compIdx, phi);
|
||||
@@ -566,7 +566,7 @@ protected:
|
||||
|
||||
// set last saturation (-> minus the change of the
|
||||
// satuation of the other phase)
|
||||
fs.setSaturation(/*phaseIdx=*/numPhases - 1,
|
||||
fs.setSaturation(/*phaseIdx=*/numPhases - 1,
|
||||
fs.saturation(numPhases - 1) - delta);
|
||||
|
||||
// update all fluid pressures using the capillary pressure
|
||||
@@ -574,7 +574,7 @@ protected:
|
||||
ComponentVector pC;
|
||||
MaterialLaw::capillaryPressures(pC, matParams, fs);
|
||||
for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
|
||||
fs.setPressure(phaseIdx,
|
||||
fs.setPressure(phaseIdx,
|
||||
fs.pressure(0)
|
||||
+ (pC[phaseIdx] - pC[0]));
|
||||
paramCache.updateAllPressures(fs);
|
||||
@@ -583,7 +583,7 @@ protected:
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
|
||||
fs.setDensity(phaseIdx, rho);
|
||||
|
||||
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
Scalar phi = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx);
|
||||
fs.setFugacityCoefficient(phaseIdx, compIdx, phi);
|
||||
@@ -621,7 +621,7 @@ protected:
|
||||
static void setQuantityRaw_(FluidState &fs, int pvIdx, Scalar value)
|
||||
{
|
||||
assert(pvIdx < numEq);
|
||||
|
||||
|
||||
// first pressure
|
||||
if (pvIdx < 1) {
|
||||
int phaseIdx = 0;
|
||||
@@ -633,7 +633,7 @@ protected:
|
||||
fs.setSaturation(phaseIdx, value);
|
||||
}
|
||||
// mole fractions
|
||||
else // if (pvIdx < numPhases + numPhases*numComponents)
|
||||
else // if (pvIdx < numPhases + numPhases*numComponents)
|
||||
{
|
||||
int phaseIdx = (pvIdx - numPhases)/numComponents;
|
||||
int compIdx = (pvIdx - numPhases)%numComponents;
|
||||
@@ -651,7 +651,7 @@ protected:
|
||||
else if (pvIdx < numPhases)
|
||||
return 1.0;
|
||||
// mole fractions
|
||||
else // if (pvIdx < numPhases + numPhases*numComponents)
|
||||
else // if (pvIdx < numPhases + numPhases*numComponents)
|
||||
return 1.0;
|
||||
}
|
||||
};
|
||||
|
||||
@@ -93,7 +93,7 @@ public:
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Predicts the vapor pressure for the temperature given in
|
||||
* setTP().
|
||||
@@ -201,7 +201,7 @@ public:
|
||||
else {
|
||||
// find the extrema (if they are present)
|
||||
Scalar Vmin, Vmax, pmin, pmax;
|
||||
if (findExtrema_(Vmin, Vmax,
|
||||
if (findExtrema_(Vmin, Vmax,
|
||||
pmin, pmax,
|
||||
a, b, T))
|
||||
{
|
||||
@@ -214,7 +214,7 @@ public:
|
||||
Vm = VmCubic;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Valgrind::CheckDefined(Vm);
|
||||
assert(std::isfinite(Vm) && Vm > 0);
|
||||
return Vm;
|
||||
@@ -287,7 +287,7 @@ protected:
|
||||
Scalar &Vcrit,
|
||||
Scalar a,
|
||||
Scalar b)
|
||||
{
|
||||
{
|
||||
Scalar minVm;
|
||||
Scalar maxVm;
|
||||
|
||||
@@ -328,7 +328,7 @@ protected:
|
||||
bool DUMUX_UNUSED hasExtrema = findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps);
|
||||
assert(hasExtrema);
|
||||
Scalar fStar = maxVm - minVm;
|
||||
|
||||
|
||||
// derivative of the difference between the maximum's
|
||||
// molar volume and the minimum's molar volume regarding
|
||||
// temperature
|
||||
|
||||
@@ -114,14 +114,14 @@ public:
|
||||
|
||||
// calculate delta_i (see: Reid, p. 145)
|
||||
Scalar sumMoleFractions = 0.0;
|
||||
for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx)
|
||||
for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx)
|
||||
sumMoleFractions += fs.moleFraction(phaseIdx, compJIdx);
|
||||
Scalar deltai = 2*std::sqrt(params.aPure(phaseIdx, compIdx))/params.a(phaseIdx);
|
||||
Scalar tmp = 0;
|
||||
for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
|
||||
tmp +=
|
||||
fs.moleFraction(phaseIdx, compJIdx)
|
||||
/ sumMoleFractions
|
||||
fs.moleFraction(phaseIdx, compJIdx)
|
||||
/ sumMoleFractions
|
||||
* std::sqrt(params.aPure(phaseIdx, compJIdx))
|
||||
* (1.0 - StaticParameters::interactionCoefficient(compIdx, compJIdx));
|
||||
};
|
||||
|
||||
@@ -107,8 +107,8 @@ public:
|
||||
Scalar pc = FluidSystem::criticalPressure(i);
|
||||
Scalar omega = FluidSystem::acentricFactor(i);
|
||||
Scalar Tr = temperature/FluidSystem::criticalTemperature(i);
|
||||
Scalar RTc = R*FluidSystem::criticalTemperature(i);
|
||||
|
||||
Scalar RTc = R*FluidSystem::criticalTemperature(i);
|
||||
|
||||
Scalar f_omega;
|
||||
|
||||
if (useSpe5Relations) {
|
||||
@@ -151,7 +151,7 @@ public:
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
|
||||
sumx += fs.moleFraction(phaseIdx, compIdx);
|
||||
sumx = std::max(1e-10, sumx);
|
||||
|
||||
|
||||
// Calculate the Peng-Robinson parameters of the mixture
|
||||
//
|
||||
// See: R. Reid, et al.: The Properties of Gases and Liquids,
|
||||
@@ -161,14 +161,14 @@ public:
|
||||
for (int compIIdx = 0; compIIdx < numComponents; ++compIIdx) {
|
||||
Scalar xi = fs.moleFraction(phaseIdx, compIIdx) / sumx;
|
||||
Valgrind::CheckDefined(xi);
|
||||
|
||||
|
||||
for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
|
||||
Scalar xj = fs.moleFraction(phaseIdx, compJIdx) / sumx;
|
||||
Valgrind::CheckDefined(xj);
|
||||
|
||||
|
||||
// mixing rule from Reid, page 82
|
||||
a += xi * xj * aCache_[compIIdx][compJIdx];
|
||||
|
||||
|
||||
assert(std::isfinite(a));
|
||||
}
|
||||
|
||||
@@ -176,7 +176,7 @@ public:
|
||||
b += xi * this->pureParams_[compIIdx].b();
|
||||
assert(std::isfinite(b));
|
||||
}
|
||||
|
||||
|
||||
this->setA(a);
|
||||
this->setB(b);
|
||||
Valgrind::CheckDefined(this->a());
|
||||
@@ -193,7 +193,7 @@ public:
|
||||
* this method!
|
||||
*/
|
||||
template <class FluidState>
|
||||
void updateSingleMoleFraction(const FluidState &fs,
|
||||
void updateSingleMoleFraction(const FluidState &fs,
|
||||
int compIdx)
|
||||
{
|
||||
updateMix(fs);
|
||||
@@ -234,20 +234,20 @@ protected:
|
||||
|
||||
private:
|
||||
void updateACache_()
|
||||
{
|
||||
{
|
||||
for (int compIIdx = 0; compIIdx < numComponents; ++ compIIdx) {
|
||||
for (int compJIdx = 0; compJIdx < numComponents; ++ compJIdx) {
|
||||
// interaction coefficient as given in SPE5
|
||||
Scalar Psi = FluidSystem::interactionCoefficient(compIIdx, compJIdx);
|
||||
|
||||
aCache_[compIIdx][compJIdx] =
|
||||
aCache_[compIIdx][compJIdx] =
|
||||
std::sqrt(this->pureParams_[compIIdx].a()
|
||||
* this->pureParams_[compJIdx].a())
|
||||
* (1 - Psi);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Scalar aCache_[numComponents][numComponents];
|
||||
};
|
||||
|
||||
|
||||
@@ -45,7 +45,7 @@ namespace Dumux
|
||||
*
|
||||
* \sa VanGenuchten, VanGenuchtenThreephase
|
||||
*/
|
||||
template <class ScalarT, class ParamsT = ParkerVanGen3PParams<ScalarT> >
|
||||
template <class ScalarT, class ParamsT = ParkerVanGen3PParams<ScalarT> >
|
||||
class ParkerVanGen3P
|
||||
{
|
||||
|
||||
@@ -65,9 +65,9 @@ public:
|
||||
static Scalar pCGW(const Params ¶ms, Scalar Sw)
|
||||
{
|
||||
/*
|
||||
Sw = wetting phase saturation, or,
|
||||
sum of wetting phase saturations
|
||||
alpha : VanGenuchten-alpha
|
||||
Sw = wetting phase saturation, or,
|
||||
sum of wetting phase saturations
|
||||
alpha : VanGenuchten-alpha
|
||||
this function is just copied from MUFTE/pml/constrel3p3cni.c
|
||||
that is why variable names do not yet fulfill Dumux rules, TODO Change */
|
||||
|
||||
@@ -109,9 +109,9 @@ public:
|
||||
static Scalar pCNW(const Params ¶ms, Scalar Sw)
|
||||
{
|
||||
/*
|
||||
Sw = wetting phase saturation, or,
|
||||
sum of wetting phase saturations
|
||||
alpha : VanGenuchten-alpha
|
||||
Sw = wetting phase saturation, or,
|
||||
sum of wetting phase saturations
|
||||
alpha : VanGenuchten-alpha
|
||||
this function is just copied from MUFTE/pml/constrel3p3cni.c
|
||||
that is why variable names do not yet fulfill Dumux rules, TODO Change */
|
||||
|
||||
@@ -153,8 +153,8 @@ public:
|
||||
static Scalar pCGN(const Params ¶ms, Scalar St)
|
||||
{
|
||||
/*
|
||||
St = sum of wetting (liquid) phase saturations
|
||||
alpha : VanGenuchten-alpha
|
||||
St = sum of wetting (liquid) phase saturations
|
||||
alpha : VanGenuchten-alpha
|
||||
this function is just copied from MUFTE/pml/constrel3p3cni.c
|
||||
that is why variable names do not yet fulfill Dumux rules, TODO Change */
|
||||
|
||||
@@ -253,7 +253,7 @@ public:
|
||||
*/
|
||||
static Scalar krw(const Params ¶ms, Scalar saturation, Scalar Sn, Scalar Sg)
|
||||
{
|
||||
|
||||
|
||||
//transformation to effective saturation
|
||||
Scalar Se = (saturation - params.Swr()) / (1-params.Swr());
|
||||
|
||||
@@ -365,7 +365,7 @@ public:
|
||||
};
|
||||
|
||||
/*
|
||||
* \brief the basis for calculating adsorbed NAPL in storage term
|
||||
* \brief the basis for calculating adsorbed NAPL in storage term
|
||||
* \param bulk density of porous medium, adsorption coefficient
|
||||
*/
|
||||
static Scalar bulkDensTimesAdsorpCoeff (const Params ¶ms)
|
||||
|
||||
@@ -208,7 +208,7 @@ public:
|
||||
*/
|
||||
Scalar KdNAPL() const
|
||||
{ return KdNAPL_; }
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Set the adsorption coefficient
|
||||
*/
|
||||
@@ -224,7 +224,7 @@ private:
|
||||
Scalar Snr_;
|
||||
Scalar Sgr_;
|
||||
Scalar Swrx_; /* (Sw+Sn)_r */
|
||||
|
||||
|
||||
Scalar KdNAPL_;
|
||||
Scalar rhoBulk_;
|
||||
|
||||
|
||||
@@ -52,11 +52,11 @@ public:
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
|
||||
moleFraction_[phaseIdx][compIdx] = 0;
|
||||
|
||||
|
||||
averageMolarMass_[phaseIdx] = 0;
|
||||
sumMoleFractions_[phaseIdx] = 0;
|
||||
}
|
||||
|
||||
|
||||
// make everything undefined so that valgrind will complain
|
||||
Valgrind::SetUndefined(*this);
|
||||
}
|
||||
@@ -255,7 +255,7 @@ public:
|
||||
* \brief Set the mole fraction of a component in a phase []
|
||||
*/
|
||||
void setMoleFraction(int phaseIdx, int compIdx, Scalar value)
|
||||
{
|
||||
{
|
||||
Valgrind::CheckDefined(value);
|
||||
Valgrind::SetDefined(sumMoleFractions_[phaseIdx]);
|
||||
Valgrind::SetDefined(averageMolarMass_[phaseIdx]);
|
||||
@@ -265,13 +265,13 @@ public:
|
||||
Scalar delta = value - moleFraction_[phaseIdx][compIdx];
|
||||
|
||||
moleFraction_[phaseIdx][compIdx] = value;
|
||||
|
||||
|
||||
sumMoleFractions_[phaseIdx] += delta;
|
||||
averageMolarMass_[phaseIdx] += delta*FluidSystem::molarMass(compIdx);
|
||||
}
|
||||
else {
|
||||
else {
|
||||
moleFraction_[phaseIdx][compIdx] = value;
|
||||
|
||||
|
||||
// re-calculate the mean molar mass
|
||||
sumMoleFractions_[phaseIdx] = 0.0;
|
||||
averageMolarMass_[phaseIdx] = 0.0;
|
||||
|
||||
@@ -49,16 +49,16 @@ public:
|
||||
enum { numComponents = FluidSystem::numComponents };
|
||||
|
||||
NonEquilibriumFluidState()
|
||||
{
|
||||
{
|
||||
// set the composition to 0
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
|
||||
moleFraction_[phaseIdx][compIdx] = 0;
|
||||
|
||||
|
||||
averageMolarMass_[phaseIdx] = 0;
|
||||
sumMoleFractions_[phaseIdx] = 0;
|
||||
}
|
||||
|
||||
|
||||
// make everything undefined so that valgrind will complain
|
||||
Valgrind::SetUndefined(*this);
|
||||
}
|
||||
@@ -205,7 +205,7 @@ public:
|
||||
* \brief Set the mole fraction of a component in a phase []
|
||||
*/
|
||||
void setMoleFraction(int phaseIdx, int compIdx, Scalar value)
|
||||
{
|
||||
{
|
||||
Valgrind::CheckDefined(value);
|
||||
Valgrind::SetDefined(sumMoleFractions_[phaseIdx]);
|
||||
Valgrind::SetDefined(averageMolarMass_[phaseIdx]);
|
||||
@@ -213,15 +213,15 @@ public:
|
||||
|
||||
if (std::isfinite(averageMolarMass_[phaseIdx])) {
|
||||
Scalar delta = value - moleFraction_[phaseIdx][compIdx];
|
||||
|
||||
|
||||
moleFraction_[phaseIdx][compIdx] = value;
|
||||
|
||||
|
||||
sumMoleFractions_[phaseIdx] += delta;
|
||||
averageMolarMass_[phaseIdx] += delta*FluidSystem::molarMass(compIdx);
|
||||
}
|
||||
else {
|
||||
else {
|
||||
moleFraction_[phaseIdx][compIdx] = value;
|
||||
|
||||
|
||||
// re-calculate the mean molar mass
|
||||
sumMoleFractions_[phaseIdx] = 0.0;
|
||||
averageMolarMass_[phaseIdx] = 0.0;
|
||||
|
||||
@@ -46,7 +46,7 @@ class SaturationOverlayFluidState
|
||||
public:
|
||||
enum { numPhases = FluidState::numPhases };
|
||||
enum { numComponents = FluidState::numComponents };
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Constructor
|
||||
*
|
||||
@@ -54,7 +54,7 @@ public:
|
||||
* argument, so it initially behaves exactly like the underlying
|
||||
* fluid state.
|
||||
*/
|
||||
SaturationOverlayFluidState(const FluidState &fs)
|
||||
SaturationOverlayFluidState(const FluidState &fs)
|
||||
: fs_(&fs)
|
||||
{
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
|
||||
|
||||
@@ -322,7 +322,7 @@ public:
|
||||
{
|
||||
DUNE_THROW(Dune::InvalidStateException, "Not applicable: Diffusion coefficients");
|
||||
}
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Given a phase's composition, temperature and pressure,
|
||||
* return the binary diffusion coefficient for components
|
||||
@@ -343,7 +343,7 @@ public:
|
||||
{
|
||||
DUNE_THROW(Dune::InvalidStateException, "Not applicable: Binary diffusion coefficients");
|
||||
}
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Given a phase's composition, temperature, pressure and
|
||||
* density, calculate its specific enthalpy [J/kg].
|
||||
@@ -403,7 +403,7 @@ public:
|
||||
Scalar pressure = fluidState.pressure(phaseIdx);
|
||||
return Fluid::heatCapacity(temperature, pressure);
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
|
||||
} // end namepace
|
||||
|
||||
@@ -376,7 +376,7 @@ public:
|
||||
p = - 1e100;
|
||||
Valgrind::SetUndefined(p);
|
||||
}
|
||||
|
||||
|
||||
|
||||
Scalar sumMoleFrac = 0;
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
|
||||
@@ -385,7 +385,7 @@ public:
|
||||
if (phaseIdx == lPhaseIdx)
|
||||
{
|
||||
//return 1000.;
|
||||
if (!useComplexRelations)
|
||||
if (!useComplexRelations)
|
||||
// assume pure water
|
||||
return H2O::liquidDensity(T, p);
|
||||
else
|
||||
@@ -404,7 +404,7 @@ public:
|
||||
}
|
||||
else if (phaseIdx == gPhaseIdx)
|
||||
{
|
||||
//return 1.189;
|
||||
//return 1.189;
|
||||
if (!useComplexRelations)
|
||||
// for the gas phase assume an ideal gas
|
||||
return
|
||||
@@ -446,7 +446,7 @@ public:
|
||||
|
||||
if (phaseIdx == lPhaseIdx)
|
||||
{
|
||||
return 0.0013;
|
||||
return 0.0013;
|
||||
// assume pure water for the liquid phase
|
||||
// TODO: viscosity of mixture
|
||||
// couldn't find a way to solve the mixture problem
|
||||
@@ -454,7 +454,7 @@ public:
|
||||
}
|
||||
else if (phaseIdx == gPhaseIdx)
|
||||
{
|
||||
return 1.8e-5;
|
||||
return 1.8e-5;
|
||||
if(!useComplexRelations){
|
||||
return Air::gasViscosity(T, p);
|
||||
}
|
||||
|
||||
@@ -86,7 +86,7 @@ public:
|
||||
* tables that always work.
|
||||
*/
|
||||
static void init()
|
||||
{
|
||||
{
|
||||
init(/*tempMin=*/273.15,
|
||||
/*tempMax=*/623.15,
|
||||
/*numTemp=*/100,
|
||||
@@ -133,7 +133,7 @@ public:
|
||||
|
||||
static bool isIdealGas(int phaseIdx)
|
||||
{ return phaseIdx == gPhaseIdx && H2O::gasIsIdeal() && Air::gasIsIdeal() && NAPL::gasIsIdeal(); }
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Returns true if and only if a fluid phase is assumed to
|
||||
* be an ideal mixture.
|
||||
|
||||
@@ -89,7 +89,7 @@ public:
|
||||
|
||||
static bool isIdealGas(int phaseIdx)
|
||||
{ return phaseIdx == gPhaseIdx && H2O::gasIsIdeal() && Air::gasIsIdeal() && NAPL::gasIsIdeal(); }
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Returns true if and only if a fluid phase is assumed to
|
||||
* be an ideal mixture.
|
||||
|
||||
@@ -360,7 +360,7 @@ public:
|
||||
Scalar sumMoleFrac = 0;
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
|
||||
sumMoleFrac += fluidState.moleFraction(phaseIdx, compIdx);
|
||||
|
||||
|
||||
// liquid phase
|
||||
if (phaseIdx == lPhaseIdx) {
|
||||
if (!useComplexRelations)
|
||||
@@ -375,7 +375,7 @@ public:
|
||||
// this assumes each nitrogen molecule displaces exactly one
|
||||
// water molecule in the liquid
|
||||
return
|
||||
clH2O
|
||||
clH2O
|
||||
* (H2O::molarMass()*fluidState.moleFraction(lPhaseIdx, H2OIdx)
|
||||
+
|
||||
N2::molarMass()*fluidState.moleFraction(lPhaseIdx, N2Idx))
|
||||
@@ -419,7 +419,7 @@ public:
|
||||
// assume pure water for the liquid phase
|
||||
return H2O::liquidViscosity(T, p);
|
||||
}
|
||||
|
||||
|
||||
// gas phase
|
||||
if (!useComplexRelations)
|
||||
{
|
||||
@@ -444,7 +444,7 @@ public:
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
|
||||
sumx += fluidState.moleFraction(phaseIdx, compIdx);
|
||||
sumx = std::max(1e-10, sumx);
|
||||
|
||||
|
||||
for (int i = 0; i < numComponents; ++i) {
|
||||
Scalar divisor = 0;
|
||||
for (int j = 0; j < numComponents; ++j) {
|
||||
@@ -469,9 +469,9 @@ public:
|
||||
* fraction \f$x^\kappa_\alpha\f$ by means of the relation
|
||||
*
|
||||
* \f[
|
||||
f^\kappa_\alpha = \phi^\kappa_\alpha\;x^\kappa_\alpha\;p_\alpha
|
||||
f^\kappa_\alpha = \phi^\kappa_\alpha\;x^\kappa_\alpha\;p_\alpha
|
||||
\f]
|
||||
* where \f$p_\alpha\f$ is the pressure of the fluid phase.
|
||||
* where \f$p_\alpha\f$ is the pressure of the fluid phase.
|
||||
*
|
||||
* The quantity "fugacity" itself is just an other way to express
|
||||
* the chemical potential \f$\zeta^\kappa_\alpha\f$ of the
|
||||
@@ -507,7 +507,7 @@ public:
|
||||
|
||||
// for the gas phase, assume an ideal gas when it comes to
|
||||
// fugacity (-> fugacity == partial pressure)
|
||||
return 1.0;
|
||||
return 1.0;
|
||||
}
|
||||
|
||||
|
||||
@@ -629,7 +629,7 @@ public:
|
||||
// "see" the molecules of the other component, which means
|
||||
// that the total specific enthalpy is the sum of the
|
||||
// "partial specific enthalpies" of the components.
|
||||
Scalar hH2O =
|
||||
Scalar hH2O =
|
||||
fluidState.massFraction(gPhaseIdx, H2OIdx)
|
||||
* H2O::gasEnthalpy(T,
|
||||
p*fluidState.moleFraction(gPhaseIdx, H2OIdx));
|
||||
@@ -683,7 +683,7 @@ public:
|
||||
Scalar pressure = fluidState.pressure(phaseIdx);
|
||||
Scalar partialPressure = pressure * xH2O;
|
||||
|
||||
Scalar lambdaH2O =
|
||||
Scalar lambdaH2O =
|
||||
xH2O
|
||||
* H2O::gasThermalConductivity(temperature, partialPressure);
|
||||
return lambdaN2 + lambdaH2O;
|
||||
@@ -715,13 +715,13 @@ public:
|
||||
// component
|
||||
|
||||
Scalar c_pN2;
|
||||
Scalar c_pH2O;
|
||||
Scalar c_pH2O;
|
||||
// let the water and nitrogen components do things their own way
|
||||
if (useComplexRelations) {
|
||||
c_pN2 = N2::gasHeatCapacity(fluidState.temperature(phaseIdx),
|
||||
fluidState.pressure(phaseIdx)
|
||||
* fluidState.moleFraction(phaseIdx, N2Idx));
|
||||
|
||||
|
||||
c_pH2O = H2O::gasHeatCapacity(fluidState.temperature(phaseIdx),
|
||||
fluidState.pressure(phaseIdx)
|
||||
* fluidState.moleFraction(phaseIdx, H2OIdx));
|
||||
|
||||
@@ -38,11 +38,11 @@ class ParameterCacheBase
|
||||
public:
|
||||
enum ExceptQuantities {
|
||||
None = 0,
|
||||
Temperature = 1,
|
||||
Temperature = 1,
|
||||
Pressure = 2,
|
||||
Composition = 4
|
||||
};
|
||||
|
||||
|
||||
ParameterCacheBase()
|
||||
{}
|
||||
|
||||
@@ -84,7 +84,7 @@ public:
|
||||
* changed between two update*() calls. If more changed, call
|
||||
* updatePhase()!
|
||||
*/
|
||||
template <class FluidState>
|
||||
template <class FluidState>
|
||||
void updateTemperature(const FluidState &fs, int phaseIdx)
|
||||
{
|
||||
asImp_().updatePhase(fs, phaseIdx);
|
||||
|
||||
@@ -81,10 +81,10 @@ public:
|
||||
static const int wPhaseIdx = 1;
|
||||
//! Index of the oil phase
|
||||
static const int oPhaseIdx = 2;
|
||||
|
||||
|
||||
//! The component for pure water to be used
|
||||
typedef Dumux::H2O<Scalar> H2O;
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Return the human readable name of a fluid phase
|
||||
*/
|
||||
@@ -158,7 +158,7 @@ public:
|
||||
|
||||
//! Number of components in the fluid system
|
||||
static const int numComponents = 7;
|
||||
|
||||
|
||||
static const int H2OIdx = 0;
|
||||
static const int C1Idx = 1;
|
||||
static const int C3Idx = 2;
|
||||
@@ -185,7 +185,7 @@ public:
|
||||
assert(0 <= compIdx && compIdx < numComponents);
|
||||
return name[compIdx];
|
||||
}
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Return the molar mass of a component in [kg/mol].
|
||||
*/
|
||||
@@ -338,11 +338,11 @@ public:
|
||||
int phaseIdx)
|
||||
{
|
||||
assert(0 <= phaseIdx && phaseIdx <= numPhases);
|
||||
|
||||
|
||||
if (phaseIdx == gPhaseIdx) {
|
||||
// given by SPE-5 in table on page 64. we use a constant
|
||||
// viscosity, though...
|
||||
return 0.0170e-2 * 0.1;
|
||||
return 0.0170e-2 * 0.1;
|
||||
}
|
||||
else if (phaseIdx == wPhaseIdx)
|
||||
// given by SPE-5: 0.7 centi-Poise = 0.0007 Pa s
|
||||
@@ -351,7 +351,7 @@ public:
|
||||
assert(phaseIdx == oPhaseIdx);
|
||||
// given by SPE-5 in table on page 64. we use a constant
|
||||
// viscosity, though...
|
||||
return 0.208e-2 * 0.1;
|
||||
return 0.208e-2 * 0.1;
|
||||
}
|
||||
}
|
||||
|
||||
@@ -369,7 +369,7 @@ public:
|
||||
template <class FluidState>
|
||||
static Scalar fugacityCoefficient(const FluidState &fs,
|
||||
const ParameterCache ¶mCache,
|
||||
int phaseIdx,
|
||||
int phaseIdx,
|
||||
int compIdx)
|
||||
{
|
||||
assert(0 <= phaseIdx && phaseIdx <= numPhases);
|
||||
@@ -378,16 +378,16 @@ public:
|
||||
if (phaseIdx == oPhaseIdx || phaseIdx == gPhaseIdx)
|
||||
return PengRobinsonMixture::computeFugacityCoefficient(fs,
|
||||
paramCache,
|
||||
phaseIdx,
|
||||
phaseIdx,
|
||||
compIdx);
|
||||
else {
|
||||
assert(phaseIdx == wPhaseIdx);
|
||||
return
|
||||
return
|
||||
henryCoeffWater_(compIdx, fs.temperature(wPhaseIdx))
|
||||
/ fs.pressure(wPhaseIdx);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Calculate the binary molecular diffusion coefficient for
|
||||
@@ -396,7 +396,7 @@ public:
|
||||
* Molecular diffusion of a compoent \f$\kappa\f$ is caused by a
|
||||
* gradient of the chemical potential and follows the law
|
||||
*
|
||||
* \f[ J = - D \grad mu_\kappa \f]
|
||||
* \f[ J = - D \grad mu_\kappa \f]
|
||||
*
|
||||
* where \f$\mu_\kappa\\f$ is the component's chemical potential,
|
||||
* \f$D\f$ is the diffusion coefficient and \f$J\f$ is the
|
||||
@@ -456,14 +456,14 @@ public:
|
||||
{ DUNE_THROW(Dune::NotImplemented, "Heat capacities"); }
|
||||
|
||||
|
||||
private:
|
||||
private:
|
||||
static Scalar henryCoeffWater_(int compIdx, Scalar temperature)
|
||||
{
|
||||
// use henry's law for the solutes and the vapor pressure for
|
||||
// the solvent.
|
||||
switch (compIdx) {
|
||||
case H2OIdx: return H2O::vaporPressure(temperature);
|
||||
|
||||
|
||||
// the values of the Henry constant for the solutes have
|
||||
// are faked so far...
|
||||
case C1Idx: return 5.57601e+09;
|
||||
|
||||
@@ -44,11 +44,11 @@ namespace Dumux
|
||||
*/
|
||||
template <class Scalar, class FluidSystem>
|
||||
class Spe5ParameterCache
|
||||
: public Dumux::ParameterCacheBase<Spe5ParameterCache<Scalar, FluidSystem> >
|
||||
: public Dumux::ParameterCacheBase<Spe5ParameterCache<Scalar, FluidSystem> >
|
||||
{
|
||||
typedef Spe5ParameterCache<Scalar, FluidSystem> ThisType;
|
||||
typedef Dumux::ParameterCacheBase<ThisType> ParentType;
|
||||
|
||||
|
||||
typedef Dumux::PengRobinson<Scalar> PengRobinson;
|
||||
|
||||
enum { numPhases = FluidSystem::numPhases };
|
||||
@@ -78,10 +78,10 @@ public:
|
||||
* calculate some quantities for the phase.
|
||||
*/
|
||||
template <class FluidState>
|
||||
void updatePhase(const FluidState &fs,
|
||||
int phaseIdx,
|
||||
void updatePhase(const FluidState &fs,
|
||||
int phaseIdx,
|
||||
int except = ParentType::None)
|
||||
{
|
||||
{
|
||||
updateEosParams(fs, phaseIdx, except);
|
||||
|
||||
// if we don't need to recalculate the molar volume, we exit
|
||||
@@ -92,7 +92,7 @@ public:
|
||||
// update the phase's molar volume
|
||||
updateMolarVolume_(fs, phaseIdx);
|
||||
}
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Update all cached parameters of a specific fluid phase
|
||||
* which depend on the mole fraction of a single component
|
||||
@@ -111,7 +111,7 @@ public:
|
||||
oilPhaseParams_.updateSingleMoleFraction(fs, compIdx);
|
||||
else if (phaseIdx == gPhaseIdx)
|
||||
gasPhaseParams_.updateSingleMoleFraction(fs, compIdx);
|
||||
|
||||
|
||||
// update the phase's molar volume
|
||||
updateMolarVolume_(fs, phaseIdx);
|
||||
}
|
||||
@@ -120,13 +120,13 @@ public:
|
||||
* \brief The Peng-Robinson attractive parameter for a phase.
|
||||
*/
|
||||
Scalar a(int phaseIdx) const
|
||||
{
|
||||
{
|
||||
switch (phaseIdx)
|
||||
{
|
||||
case oPhaseIdx: return oilPhaseParams_.a();
|
||||
case gPhaseIdx: return gasPhaseParams_.a();
|
||||
default:
|
||||
DUNE_THROW(Dune::InvalidStateException,
|
||||
default:
|
||||
DUNE_THROW(Dune::InvalidStateException,
|
||||
"The a() parameter is only defined for "
|
||||
"oil and gas phases");
|
||||
};
|
||||
@@ -136,13 +136,13 @@ public:
|
||||
* \brief The Peng-Robinson covolume for a phase.
|
||||
*/
|
||||
Scalar b(int phaseIdx) const
|
||||
{
|
||||
{
|
||||
switch (phaseIdx)
|
||||
{
|
||||
case oPhaseIdx: return oilPhaseParams_.b();
|
||||
case gPhaseIdx: return gasPhaseParams_.b();
|
||||
default:
|
||||
DUNE_THROW(Dune::InvalidStateException,
|
||||
default:
|
||||
DUNE_THROW(Dune::InvalidStateException,
|
||||
"The b() parameter is only defined for "
|
||||
"oil and gas phases");
|
||||
};
|
||||
@@ -154,13 +154,13 @@ public:
|
||||
* phase.
|
||||
*/
|
||||
Scalar aPure(int phaseIdx, int compIdx) const
|
||||
{
|
||||
{
|
||||
switch (phaseIdx)
|
||||
{
|
||||
case oPhaseIdx: return oilPhaseParams_.pureParams(compIdx).a();
|
||||
case gPhaseIdx: return gasPhaseParams_.pureParams(compIdx).a();
|
||||
default:
|
||||
DUNE_THROW(Dune::InvalidStateException,
|
||||
default:
|
||||
DUNE_THROW(Dune::InvalidStateException,
|
||||
"The a() parameter is only defined for "
|
||||
"oil and gas phases");
|
||||
};
|
||||
@@ -171,13 +171,13 @@ public:
|
||||
* the same temperature and pressure of the phase.
|
||||
*/
|
||||
Scalar bPure(int phaseIdx, int compIdx) const
|
||||
{
|
||||
{
|
||||
switch (phaseIdx)
|
||||
{
|
||||
case oPhaseIdx: return oilPhaseParams_.pureParams(compIdx).b();
|
||||
case gPhaseIdx: return gasPhaseParams_.pureParams(compIdx).b();
|
||||
default:
|
||||
DUNE_THROW(Dune::InvalidStateException,
|
||||
default:
|
||||
DUNE_THROW(Dune::InvalidStateException,
|
||||
"The b() parameter is only defined for "
|
||||
"oil and gas phases");
|
||||
};
|
||||
@@ -188,8 +188,8 @@ public:
|
||||
*/
|
||||
Scalar molarVolume(int phaseIdx) const
|
||||
{ assert(VmUpToDate_[phaseIdx]); return Vm_[phaseIdx]; }
|
||||
|
||||
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Returns the Peng-Robinson mixture parameters for the oil
|
||||
* phase.
|
||||
@@ -203,16 +203,16 @@ public:
|
||||
*/
|
||||
const GasPhaseParams &gasPhaseParams() const
|
||||
{ return gasPhaseParams_; }
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Update all parameters required by the equation of state to
|
||||
* calculate some quantities for the phase.
|
||||
*/
|
||||
template <class FluidState>
|
||||
void updateEosParams(const FluidState &fs,
|
||||
void updateEosParams(const FluidState &fs,
|
||||
int phaseIdx,
|
||||
int exceptQuantities = ParentType::None)
|
||||
{
|
||||
{
|
||||
if (!(exceptQuantities & ParentType::Temperature))
|
||||
{
|
||||
updatePure_(fs, phaseIdx);
|
||||
@@ -231,7 +231,7 @@ public:
|
||||
|
||||
protected:
|
||||
/*!
|
||||
* \brief Update all parameters of a phase which only depend on
|
||||
* \brief Update all parameters of a phase which only depend on
|
||||
* temperature and/or pressure.
|
||||
*
|
||||
* This usually means the parameters for the pure components.
|
||||
@@ -263,11 +263,11 @@ protected:
|
||||
Valgrind::CheckDefined(fs.averageMolarMass(phaseIdx));
|
||||
switch (phaseIdx)
|
||||
{
|
||||
case oPhaseIdx:
|
||||
case oPhaseIdx:
|
||||
oilPhaseParams_.updateMix(fs);
|
||||
break;
|
||||
case gPhaseIdx:
|
||||
gasPhaseParams_.updateMix(fs);
|
||||
gasPhaseParams_.updateMix(fs);
|
||||
break;
|
||||
case wPhaseIdx:
|
||||
break;
|
||||
@@ -275,11 +275,11 @@ protected:
|
||||
}
|
||||
|
||||
template <class FluidState>
|
||||
void updateMolarVolume_(const FluidState &fs,
|
||||
void updateMolarVolume_(const FluidState &fs,
|
||||
int phaseIdx)
|
||||
{
|
||||
VmUpToDate_[phaseIdx] = true;
|
||||
|
||||
|
||||
// calculate molar volume of the phase (we will need this for the
|
||||
// fugacity coefficients and the density anyway)
|
||||
switch (phaseIdx) {
|
||||
@@ -289,7 +289,7 @@ protected:
|
||||
// molar volume appears in basically every quantity the fluid
|
||||
// system can get queried, so it is okay to calculate it
|
||||
// here...
|
||||
Vm_[gPhaseIdx] =
|
||||
Vm_[gPhaseIdx] =
|
||||
PengRobinson::computeMolarVolume(fs,
|
||||
*this,
|
||||
phaseIdx,
|
||||
@@ -301,7 +301,7 @@ protected:
|
||||
// molar volume appears in basically every quantity the fluid
|
||||
// system can get queried, so it is okay to calculate it
|
||||
// here...
|
||||
Vm_[oPhaseIdx] =
|
||||
Vm_[oPhaseIdx] =
|
||||
PengRobinson::computeMolarVolume(fs,
|
||||
*this,
|
||||
phaseIdx,
|
||||
@@ -316,7 +316,7 @@ protected:
|
||||
// Water compressibility is specified as 3.3e-6 per psi
|
||||
// overpressure, where 1 psi = 6894.7573 Pa
|
||||
Scalar overPressure = fs.pressure(wPhaseIdx) - 1.013e5; // [Pa]
|
||||
Scalar waterDensity =
|
||||
Scalar waterDensity =
|
||||
stockTankWaterDensity * (1 + 3.3e-6*overPressure/6894.7573);
|
||||
|
||||
// convert water density [kg/m^3] to molar volume [m^3/mol]
|
||||
|
||||
@@ -23,7 +23,7 @@
|
||||
* \file
|
||||
*
|
||||
* \brief This file provides the actual code for the fluid systems
|
||||
* test.
|
||||
* test.
|
||||
*
|
||||
* It is not directly in test_fluidsystems.cc so that external modules
|
||||
* like dumux-devel can use it easily
|
||||
@@ -53,16 +53,16 @@
|
||||
|
||||
// this is a fluid state which makes sure that only the quantities
|
||||
// allowed are accessed
|
||||
template <class Scalar,
|
||||
template <class Scalar,
|
||||
class FluidSystem,
|
||||
class BaseFluidState = Dumux::CompositionalFluidState<Scalar, FluidSystem> >
|
||||
class HairSplittingFluidState
|
||||
class HairSplittingFluidState
|
||||
: protected BaseFluidState
|
||||
{
|
||||
public:
|
||||
enum { numPhases = FluidSystem::numPhases };
|
||||
enum { numComponents = FluidSystem::numComponents };
|
||||
|
||||
|
||||
HairSplittingFluidState()
|
||||
{
|
||||
// set some fake values
|
||||
@@ -73,16 +73,16 @@ public:
|
||||
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
BaseFluidState::setMoleFraction(phaseIdx, compIdx, 1.0 / numComponents);
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// initially, do not allow anything
|
||||
allowTemperature(false);
|
||||
allowPressure(false);
|
||||
allowComposition(false);
|
||||
allowDensity(false);
|
||||
|
||||
|
||||
// do not allow accessing any phase
|
||||
restrictToPhase(1000);
|
||||
}
|
||||
@@ -95,11 +95,11 @@ public:
|
||||
|
||||
void allowComposition(bool yesno)
|
||||
{ allowComposition_ = yesno; }
|
||||
|
||||
void allowDensity(bool yesno)
|
||||
|
||||
void allowDensity(bool yesno)
|
||||
{ allowDensity_ = yesno; }
|
||||
|
||||
void restrictToPhase(int phaseIdx)
|
||||
void restrictToPhase(int phaseIdx)
|
||||
{ restrictPhaseIdx_ = phaseIdx; }
|
||||
|
||||
Scalar temperature(int phaseIdx) const
|
||||
@@ -108,31 +108,31 @@ public:
|
||||
assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx);
|
||||
return BaseFluidState::temperature(phaseIdx);
|
||||
}
|
||||
|
||||
|
||||
Scalar pressure(int phaseIdx) const
|
||||
{
|
||||
{
|
||||
assert(allowPressure_);
|
||||
assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx);
|
||||
return BaseFluidState::pressure(phaseIdx);
|
||||
}
|
||||
|
||||
Scalar moleFraction(int phaseIdx, int compIdx) const
|
||||
{
|
||||
assert(allowComposition_);
|
||||
{
|
||||
assert(allowComposition_);
|
||||
assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx);
|
||||
return BaseFluidState::moleFraction(phaseIdx, compIdx);
|
||||
}
|
||||
|
||||
Scalar massFraction(int phaseIdx, int compIdx) const
|
||||
{
|
||||
assert(allowComposition_);
|
||||
{
|
||||
assert(allowComposition_);
|
||||
assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx);
|
||||
return BaseFluidState::massFraction(phaseIdx, compIdx);
|
||||
}
|
||||
|
||||
Scalar averageMolarMass(int phaseIdx) const
|
||||
{
|
||||
assert(allowComposition_);
|
||||
assert(allowComposition_);
|
||||
assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx);
|
||||
return BaseFluidState::averageMolarMass(phaseIdx);
|
||||
}
|
||||
@@ -153,7 +153,7 @@ public:
|
||||
|
||||
Scalar molarVolume(int phaseIdx) const
|
||||
{
|
||||
assert(allowDensity_);
|
||||
assert(allowDensity_);
|
||||
assert(restrictPhaseIdx_ < 0 || restrictPhaseIdx_ == phaseIdx);
|
||||
return BaseFluidState::molarVolume(phaseIdx);
|
||||
}
|
||||
@@ -166,7 +166,7 @@ public:
|
||||
}
|
||||
|
||||
Scalar saturation(int phaseIdx) const
|
||||
{
|
||||
{
|
||||
assert(false);
|
||||
return BaseFluidState::saturation(phaseIdx);
|
||||
}
|
||||
@@ -196,7 +196,7 @@ public:
|
||||
}
|
||||
|
||||
Scalar viscosity(int phaseIdx) const
|
||||
{
|
||||
{
|
||||
assert(false);
|
||||
return BaseFluidState::viscosity(phaseIdx);
|
||||
}
|
||||
@@ -219,7 +219,7 @@ void checkFluidState(const BaseFluidState &fs)
|
||||
|
||||
// a fluid state must provide a checkDefined() method
|
||||
fs.checkDefined();
|
||||
|
||||
|
||||
// make sure the fluid state provides all mandatory methods
|
||||
while (false) {
|
||||
Scalar DUMUX_UNUSED val;
|
||||
@@ -296,7 +296,7 @@ void checkFluidSystem()
|
||||
try { val = FluidSystem::enthalpy(fs, paramCache, phaseIdx); } catch (...) {};
|
||||
try { val = FluidSystem::heatCapacity(fs, paramCache, phaseIdx); } catch (...) {};
|
||||
try { val = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx); } catch (...) {};
|
||||
|
||||
|
||||
for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
|
||||
fs.allowComposition(!FluidSystem::isIdealMixture(phaseIdx));
|
||||
try { val = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx); } catch (...) {};
|
||||
@@ -314,7 +314,7 @@ void checkFluidSystem()
|
||||
bool DUMUX_UNUSED bVal = FluidSystem::isLiquid(phaseIdx);
|
||||
bVal = FluidSystem::isIdealGas(phaseIdx);
|
||||
}
|
||||
|
||||
|
||||
// test for componentName()
|
||||
for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
|
||||
val = FluidSystem::molarMass(compIdx);
|
||||
|
||||
@@ -58,11 +58,11 @@ int main()
|
||||
|
||||
typedef Dumux::LiquidPhase<Scalar, H2O> Liquid;
|
||||
typedef Dumux::GasPhase<Scalar, N2> Gas;
|
||||
|
||||
|
||||
// check all fluid states
|
||||
{
|
||||
typedef Dumux::FluidSystems::H2ON2<Scalar, /*enableComplexRelations=*/false> FluidSystem;
|
||||
|
||||
|
||||
// CompositionalFluidState
|
||||
{ Dumux::CompositionalFluidState<Scalar, FluidSystem> fs;
|
||||
checkFluidState<Scalar>(fs); }
|
||||
|
||||
@@ -67,7 +67,7 @@ void checkSame(const FluidState &fsRef, const FluidState &fsFlash)
|
||||
// check the saturations
|
||||
error = fsRef.saturation(phaseIdx) - fsFlash.saturation(phaseIdx);
|
||||
if (std::abs(error) > 1e-6)
|
||||
std::cout << "saturation error phase " << phaseIdx << ": "
|
||||
std::cout << "saturation error phase " << phaseIdx << ": "
|
||||
<< fsFlash.saturation(phaseIdx) << " flash vs "
|
||||
<< fsRef.saturation(phaseIdx) << " reference"
|
||||
<< " error=" << error << "\n";
|
||||
@@ -76,7 +76,7 @@ void checkSame(const FluidState &fsRef, const FluidState &fsFlash)
|
||||
for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
|
||||
error = fsRef.moleFraction(phaseIdx, compIdx) - fsFlash.moleFraction(phaseIdx, compIdx);
|
||||
if (std::abs(error) > 1e-6)
|
||||
std::cout << "composition error phase " << phaseIdx << ", component " << compIdx << ": "
|
||||
std::cout << "composition error phase " << phaseIdx << ", component " << compIdx << ": "
|
||||
<< fsFlash.moleFraction(phaseIdx, compIdx) << " flash vs "
|
||||
<< fsRef.moleFraction(phaseIdx, compIdx) << " reference"
|
||||
<< " error=" << error << "\n";
|
||||
@@ -85,7 +85,7 @@ void checkSame(const FluidState &fsRef, const FluidState &fsFlash)
|
||||
}
|
||||
|
||||
template <class Scalar, class FluidSystem, class MaterialLaw, class FluidState>
|
||||
void checkImmiscibleFlash(const FluidState &fsRef,
|
||||
void checkImmiscibleFlash(const FluidState &fsRef,
|
||||
typename MaterialLaw::Params &matParams)
|
||||
{
|
||||
enum { numPhases = FluidSystem::numPhases };
|
||||
@@ -119,7 +119,7 @@ void checkImmiscibleFlash(const FluidState &fsRef,
|
||||
|
||||
|
||||
template <class Scalar, class FluidSystem, class MaterialLaw, class FluidState>
|
||||
void completeReferenceFluidState(FluidState &fs,
|
||||
void completeReferenceFluidState(FluidState &fs,
|
||||
typename MaterialLaw::Params &matParams,
|
||||
int refPhaseIdx)
|
||||
{
|
||||
@@ -127,7 +127,7 @@ void completeReferenceFluidState(FluidState &fs,
|
||||
|
||||
typedef Dumux::ComputeFromReferencePhase<Scalar, FluidSystem> ComputeFromReferencePhase;
|
||||
typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
|
||||
|
||||
|
||||
int otherPhaseIdx = 1 - refPhaseIdx;
|
||||
|
||||
// calculate the other saturation
|
||||
@@ -136,7 +136,7 @@ void completeReferenceFluidState(FluidState &fs,
|
||||
// calulate the capillary pressure
|
||||
PhaseVector pC;
|
||||
MaterialLaw::capillaryPressures(pC, matParams, fs);
|
||||
fs.setPressure(otherPhaseIdx,
|
||||
fs.setPressure(otherPhaseIdx,
|
||||
fs.pressure(refPhaseIdx)
|
||||
+ (pC[otherPhaseIdx] - pC[refPhaseIdx]));
|
||||
|
||||
@@ -170,12 +170,12 @@ int main()
|
||||
typedef MaterialLaw::Params MaterialLawParams;
|
||||
|
||||
Scalar T = 273.15 + 25;
|
||||
|
||||
|
||||
// initialize the tables of the fluid system
|
||||
Scalar Tmin = T - 1.0;
|
||||
Scalar Tmax = T + 1.0;
|
||||
int nT = 3;
|
||||
|
||||
|
||||
Scalar pmin = 0.0;
|
||||
Scalar pmax = 1.25 * 2e6;
|
||||
int np = 100;
|
||||
@@ -190,12 +190,12 @@ int main()
|
||||
matParams.setLambda(2.0);
|
||||
|
||||
ImmiscibleFluidState fsRef;
|
||||
|
||||
|
||||
// create an fluid state which is consistent
|
||||
|
||||
// set the fluid temperatures
|
||||
fsRef.setTemperature(T);
|
||||
|
||||
|
||||
////////////////
|
||||
// only liquid
|
||||
////////////////
|
||||
@@ -204,10 +204,10 @@ int main()
|
||||
// set liquid saturation and pressure
|
||||
fsRef.setSaturation(lPhaseIdx, 1.0);
|
||||
fsRef.setPressure(lPhaseIdx, 1e6);
|
||||
|
||||
|
||||
// set the remaining parameters of the reference fluid state
|
||||
completeReferenceFluidState<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams, lPhaseIdx);
|
||||
|
||||
|
||||
// check the flash calculation
|
||||
checkImmiscibleFlash<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams);
|
||||
|
||||
@@ -222,7 +222,7 @@ int main()
|
||||
|
||||
// set the remaining parameters of the reference fluid state
|
||||
completeReferenceFluidState<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams, gPhaseIdx);
|
||||
|
||||
|
||||
// check the flash calculation
|
||||
checkImmiscibleFlash<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams);
|
||||
|
||||
@@ -237,7 +237,7 @@ int main()
|
||||
|
||||
// set the remaining parameters of the reference fluid state
|
||||
completeReferenceFluidState<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams, lPhaseIdx);
|
||||
|
||||
|
||||
// check the flash calculation
|
||||
checkImmiscibleFlash<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams);
|
||||
|
||||
@@ -254,13 +254,13 @@ int main()
|
||||
|
||||
// set liquid saturation
|
||||
fsRef.setSaturation(lPhaseIdx, 0.5);
|
||||
|
||||
|
||||
// set pressure of the liquid phase
|
||||
fsRef.setPressure(lPhaseIdx, 1e6);
|
||||
|
||||
// set the remaining parameters of the reference fluid state
|
||||
completeReferenceFluidState<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams2, lPhaseIdx);
|
||||
|
||||
|
||||
// check the flash calculation
|
||||
checkImmiscibleFlash<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams2);
|
||||
|
||||
|
||||
@@ -67,7 +67,7 @@ void checkSame(const FluidState &fsRef, const FluidState &fsFlash)
|
||||
// check the saturations
|
||||
error = fsRef.saturation(phaseIdx) - fsFlash.saturation(phaseIdx);
|
||||
if (std::abs(error) > 1e-6)
|
||||
std::cout << "saturation error phase " << phaseIdx << ": "
|
||||
std::cout << "saturation error phase " << phaseIdx << ": "
|
||||
<< fsFlash.saturation(phaseIdx) << " flash vs "
|
||||
<< fsRef.saturation(phaseIdx) << " reference"
|
||||
<< " error=" << error << "\n";
|
||||
@@ -76,7 +76,7 @@ void checkSame(const FluidState &fsRef, const FluidState &fsFlash)
|
||||
for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
|
||||
error = fsRef.moleFraction(phaseIdx, compIdx) - fsFlash.moleFraction(phaseIdx, compIdx);
|
||||
if (std::abs(error) > 1e-6)
|
||||
std::cout << "composition error phase " << phaseIdx << ", component " << compIdx << ": "
|
||||
std::cout << "composition error phase " << phaseIdx << ", component " << compIdx << ": "
|
||||
<< fsFlash.moleFraction(phaseIdx, compIdx) << " flash vs "
|
||||
<< fsRef.moleFraction(phaseIdx, compIdx) << " reference"
|
||||
<< " error=" << error << "\n";
|
||||
@@ -85,7 +85,7 @@ void checkSame(const FluidState &fsRef, const FluidState &fsFlash)
|
||||
}
|
||||
|
||||
template <class Scalar, class FluidSystem, class MaterialLaw, class FluidState>
|
||||
void checkNcpFlash(const FluidState &fsRef,
|
||||
void checkNcpFlash(const FluidState &fsRef,
|
||||
typename MaterialLaw::Params &matParams)
|
||||
{
|
||||
enum { numPhases = FluidSystem::numPhases };
|
||||
@@ -119,7 +119,7 @@ void checkNcpFlash(const FluidState &fsRef,
|
||||
|
||||
|
||||
template <class Scalar, class FluidSystem, class MaterialLaw, class FluidState>
|
||||
void completeReferenceFluidState(FluidState &fs,
|
||||
void completeReferenceFluidState(FluidState &fs,
|
||||
typename MaterialLaw::Params &matParams,
|
||||
int refPhaseIdx)
|
||||
{
|
||||
@@ -127,7 +127,7 @@ void completeReferenceFluidState(FluidState &fs,
|
||||
|
||||
typedef Dumux::ComputeFromReferencePhase<Scalar, FluidSystem> ComputeFromReferencePhase;
|
||||
typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
|
||||
|
||||
|
||||
int otherPhaseIdx = 1 - refPhaseIdx;
|
||||
|
||||
// calculate the other saturation
|
||||
@@ -136,15 +136,15 @@ void completeReferenceFluidState(FluidState &fs,
|
||||
// calulate the capillary pressure
|
||||
PhaseVector pC;
|
||||
MaterialLaw::capillaryPressures(pC, matParams, fs);
|
||||
fs.setPressure(otherPhaseIdx,
|
||||
fs.setPressure(otherPhaseIdx,
|
||||
fs.pressure(refPhaseIdx)
|
||||
+ (pC[otherPhaseIdx] - pC[refPhaseIdx]));
|
||||
|
||||
|
||||
// make the fluid state consistent with local thermodynamic
|
||||
// equilibrium
|
||||
typename FluidSystem::ParameterCache paramCache;
|
||||
ComputeFromReferencePhase::solve(fs,
|
||||
paramCache,
|
||||
paramCache,
|
||||
refPhaseIdx,
|
||||
/*setViscosity=*/false,
|
||||
/*setEnthalpy=*/false);
|
||||
@@ -171,12 +171,12 @@ int main()
|
||||
typedef MaterialLaw::Params MaterialLawParams;
|
||||
|
||||
Scalar T = 273.15 + 25;
|
||||
|
||||
|
||||
// initialize the tables of the fluid system
|
||||
Scalar Tmin = T - 1.0;
|
||||
Scalar Tmax = T + 1.0;
|
||||
int nT = 3;
|
||||
|
||||
|
||||
Scalar pmin = 0.0;
|
||||
Scalar pmax = 1.25 * 2e6;
|
||||
int np = 100;
|
||||
@@ -191,12 +191,12 @@ int main()
|
||||
matParams.setLambda(2.0);
|
||||
|
||||
CompositionalFluidState fsRef;
|
||||
|
||||
|
||||
// create an fluid state which is consistent
|
||||
|
||||
// set the fluid temperatures
|
||||
fsRef.setTemperature(T);
|
||||
|
||||
|
||||
////////////////
|
||||
// only liquid
|
||||
////////////////
|
||||
@@ -204,14 +204,14 @@ int main()
|
||||
|
||||
// set liquid saturation
|
||||
fsRef.setSaturation(lPhaseIdx, 1.0);
|
||||
|
||||
|
||||
// set pressure of the liquid phase
|
||||
fsRef.setPressure(lPhaseIdx, 2e5);
|
||||
|
||||
|
||||
// set the liquid composition to pure water
|
||||
fsRef.setMoleFraction(lPhaseIdx, N2Idx, 0.0);
|
||||
fsRef.setMoleFraction(lPhaseIdx, H2OIdx, 1.0);
|
||||
|
||||
|
||||
// "complete" the fluid state
|
||||
completeReferenceFluidState<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams, lPhaseIdx);
|
||||
|
||||
@@ -224,14 +224,14 @@ int main()
|
||||
std::cout << "testing single-phase gas\n";
|
||||
// set gas saturation
|
||||
fsRef.setSaturation(gPhaseIdx, 1.0);
|
||||
|
||||
|
||||
// set pressure of the gas phase
|
||||
fsRef.setPressure(gPhaseIdx, 1e6);
|
||||
|
||||
|
||||
// set the gas composition to 99.9% nitrogen and 0.1% water
|
||||
fsRef.setMoleFraction(gPhaseIdx, N2Idx, 0.999);
|
||||
fsRef.setMoleFraction(gPhaseIdx, H2OIdx, 0.001);
|
||||
|
||||
|
||||
// "complete" the fluid state
|
||||
completeReferenceFluidState<Scalar, FluidSystem, MaterialLaw>(fsRef, matParams, gPhaseIdx);
|
||||
|
||||
@@ -246,11 +246,11 @@ int main()
|
||||
// set saturations
|
||||
fsRef.setSaturation(lPhaseIdx, 0.5);
|
||||
fsRef.setSaturation(gPhaseIdx, 0.5);
|
||||
|
||||
|
||||
// set pressures
|
||||
fsRef.setPressure(lPhaseIdx, 1e6);
|
||||
fsRef.setPressure(gPhaseIdx, 1e6);
|
||||
|
||||
|
||||
FluidSystem::ParameterCache paramCache;
|
||||
typedef Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem> MiscibleMultiPhaseComposition;
|
||||
MiscibleMultiPhaseComposition::solve(fsRef, paramCache,
|
||||
@@ -273,7 +273,7 @@ int main()
|
||||
// set gas saturation
|
||||
fsRef.setSaturation(gPhaseIdx, 0.5);
|
||||
fsRef.setSaturation(lPhaseIdx, 0.5);
|
||||
|
||||
|
||||
// set pressure of the liquid phase
|
||||
fsRef.setPressure(lPhaseIdx, 1e6);
|
||||
|
||||
@@ -281,10 +281,10 @@ int main()
|
||||
typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
|
||||
PhaseVector pC;
|
||||
MaterialLaw::capillaryPressures(pC, matParams2, fsRef);
|
||||
fsRef.setPressure(gPhaseIdx,
|
||||
fsRef.setPressure(gPhaseIdx,
|
||||
fsRef.pressure(lPhaseIdx)
|
||||
+ (pC[gPhaseIdx] - pC[lPhaseIdx]));
|
||||
|
||||
|
||||
typedef Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem> MiscibleMultiPhaseComposition;
|
||||
MiscibleMultiPhaseComposition::solve(fsRef, paramCache,
|
||||
/*setViscosity=*/false,
|
||||
|
||||
@@ -65,7 +65,7 @@ int main(int argc, char** argv)
|
||||
////////////
|
||||
// Initialize the fluid system and create the capillary pressure
|
||||
// parameters
|
||||
////////////
|
||||
////////////
|
||||
FluidSystem::init();
|
||||
|
||||
// set the parameters for the capillary pressure law
|
||||
@@ -77,7 +77,7 @@ int main(int argc, char** argv)
|
||||
|
||||
////////////
|
||||
// Create a fluid state
|
||||
////////////
|
||||
////////////
|
||||
FluidState fluidState;
|
||||
ParameterCache paramCache;
|
||||
|
||||
@@ -127,7 +127,7 @@ int main(int argc, char** argv)
|
||||
for (int i = 0; i < n; ++i) {
|
||||
Scalar minAlpha = 0.98;
|
||||
Scalar maxAlpha = 5.0;
|
||||
|
||||
|
||||
// ratio between the original and the current volume
|
||||
Scalar alpha = minAlpha + (maxAlpha - minAlpha)*i/(n - 1);
|
||||
|
||||
@@ -138,7 +138,7 @@ int main(int argc, char** argv)
|
||||
// "flash" the modified reservoir oil
|
||||
Flash::solve<MaterialLaw>(flashFluidState, paramCache, matParams, curMolarities);
|
||||
|
||||
/*
|
||||
/*
|
||||
// bring the gas and the oil to the surface
|
||||
surfaceFluidState.setPressure(gPhaseIdx, 1.013e5);
|
||||
surfaceFluidState.setPressure(oPhaseIdx, 1.013e5);
|
||||
@@ -147,7 +147,7 @@ int main(int argc, char** argv)
|
||||
paramCache.updatePhase(surfaceFluidState, oPhaseIdx);
|
||||
Scalar rhoGSurface = FluidSystem::density(surfaceFluidState, paramCache, gPhaseIdx);
|
||||
Scalar rhoOSurface = FluidSystem::density(surfaceFluidState, paramCache, gPhaseIdx);
|
||||
*/
|
||||
*/
|
||||
std::cout << alpha << " "
|
||||
<< flashFluidState.pressure(oPhaseIdx) << " "
|
||||
<< flashFluidState.saturation(gPhaseIdx) << " "
|
||||
@@ -155,6 +155,6 @@ int main(int argc, char** argv)
|
||||
<< flashFluidState.density(gPhaseIdx) << " "
|
||||
<< "\n";
|
||||
}
|
||||
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user