OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

make the constraint solvers usable with local-AD

Browse Source
This commit is contained in:
Andreas Lauser 2015-05-21 15:33:21 +02:00
parent 8aaf2bfdab
commit 662531fee8
4 changed files with 193 additions and 209 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

154
opm/material/constraintsolvers/CompositionFromFugacities.hpp
View File

@ -24,6 +24,8 @@
#ifndef OPM_COMPOSITION_FROM_FUGACITIES_HPP
#define OPM_COMPOSITION_FROM_FUGACITIES_HPP
#include <opm/material/common/MathToolbox.hpp>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
@ -39,7 +41,7 @@ namespace Opm {
* \brief Calculates the chemical equilibrium from the component
* fugacities in a phase.
*/
template <class Scalar, class FluidSystem>
template <class Scalar, class FluidSystem, class Evaluation = Scalar>
class CompositionFromFugacities
{
enum { numComponents = FluidSystem::numComponents };
@ -47,7 +49,7 @@ class CompositionFromFugacities
typedef typename FluidSystem::ParameterCache ParameterCache;
public:
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
/*!
* \brief Guess an initial value for the composition of the phase.
@ -65,8 +67,8 @@ public:
for (int i = 0; i < numComponents; ++ i) {
//std::cout << f << " -> " << mutParams.fugacity(phaseIdx, i)/f << "\n";
fluidState.setMoleFraction(phaseIdx,
i,
1.0/numComponents);
i,
1.0/numComponents);
}
}
@ -82,6 +84,8 @@ public:
int phaseIdx,
const ComponentVector &targetFug)
{
typedef MathToolbox<Evaluation> Toolbox;
// use a much more efficient method in case the phase is an
// ideal mixture
if (FluidSystem::isIdealMixture(phaseIdx)) {
@ -89,10 +93,10 @@ public:
return;
}
Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-25);
//Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-25);
// save initial composition in case something goes wrong
Dune::FieldVector<Scalar, numComponents> xInit;
Dune::FieldVector<Evaluation, numComponents> xInit;
for (int i = 0; i < numComponents; ++i) {
xInit[i] = fluidState.moleFraction(phaseIdx, i);
}
@ -102,11 +106,11 @@ public:
/////////////////////////
// Jacobian matrix
Dune::FieldMatrix<Scalar, numComponents, numComponents> J;
Dune::FieldMatrix<Evaluation, numComponents, numComponents> J;
// solution, i.e. phase composition
Dune::FieldVector<Scalar, numComponents> x;
Dune::FieldVector<Evaluation, numComponents> x;
// right hand side
Dune::FieldVector<Scalar, numComponents> b;
Dune::FieldVector<Evaluation, numComponents> b;
paramCache.updatePhase(fluidState, phaseIdx);
@ -130,7 +134,7 @@ public:
*/
// Solve J*x = b
x = 0;
x = Toolbox::createConstant(0.0);
try { J.solve(x, b); }
catch (Dune::FMatrixError e)
{ throw Opm::NumericalIssue(e.what()); }
@ -159,7 +163,7 @@ public:
Scalar relError = update_(fluidState, paramCache, x, b, phaseIdx, targetFug);
if (relError < 1e-9) {
Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
fluidState.setDensity(phaseIdx, rho);
//std::cout << "num iterations: " << nIdx << "\n";
@ -187,11 +191,11 @@ protected:
const ComponentVector &fugacities)
{
for (int i = 0; i < numComponents; ++ i) {
Scalar phi = FluidSystem::fugacityCoefficient(fluidState,
paramCache,
phaseIdx,
i);
Scalar gamma = phi * fluidState.pressure(phaseIdx);
const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
paramCache,
phaseIdx,
i);
const Evaluation& gamma = phi * fluidState.pressure(phaseIdx);
Valgrind::CheckDefined(phi);
Valgrind::CheckDefined(gamma);
Valgrind::CheckDefined(fugacities[i]);
@ -201,19 +205,21 @@ protected:
paramCache.updatePhase(fluidState, phaseIdx);
Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
fluidState.setDensity(phaseIdx, rho);
return;
}
template <class FluidState>
static Scalar linearize_(Dune::FieldMatrix<Scalar, numComponents, numComponents> &J,
Dune::FieldVector<Scalar, numComponents> &defect,
static Scalar linearize_(Dune::FieldMatrix<Evaluation, numComponents, numComponents> &J,
Dune::FieldVector<Evaluation, numComponents> &defect,
FluidState &fluidState,
ParameterCache &paramCache,
int phaseIdx,
const ComponentVector &targetFug)
{
typedef MathToolbox<Evaluation> Toolbox;
// reset jacobian
J = 0;
@ -221,15 +227,15 @@ protected:
// calculate the defect (deviation of the current fugacities
// from the target fugacities)
for (int i = 0; i < numComponents; ++ i) {
Scalar phi = FluidSystem::fugacityCoefficient(fluidState,
const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
paramCache,
phaseIdx,
i);
Scalar f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i);
const Evaluation& f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i);
fluidState.setFugacityCoefficient(phaseIdx, i, phi);
defect[i] = targetFug[i] - f;
absError = std::max(absError, std::abs(defect[i]));
absError = std::max(absError, std::abs(Toolbox::value(defect[i])));
}
// assemble jacobian matrix of the constraints for the composition
@ -242,7 +248,7 @@ protected:
// forward differences
// deviate the mole fraction of the i-th component
Scalar xI = fluidState.moleFraction(phaseIdx, i);
Evaluation xI = fluidState.moleFraction(phaseIdx, i);
fluidState.setMoleFraction(phaseIdx, i, xI + eps);
paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
@ -250,17 +256,17 @@ protected:
// fugacities
for (int j = 0; j < numComponents; ++j) {
// compute the j-th component's fugacity coefficient ...
Scalar phi = FluidSystem::fugacityCoefficient(fluidState,
paramCache,
phaseIdx,
j);
const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
paramCache,
phaseIdx,
j);
// ... and its fugacity ...
Scalar f =
const Evaluation& f =
phi *
fluidState.pressure(phaseIdx) *
fluidState.moleFraction(phaseIdx, j);
// as well as the defect for this fugacity
Scalar defJPlusEps = targetFug[j] - f;
const Evaluation& defJPlusEps = targetFug[j] - f;
// use forward differences to calculate the defect's
// derivative
@ -281,86 +287,48 @@ protected:
template <class FluidState>
static Scalar update_(FluidState &fluidState,
ParameterCache &paramCache,
Dune::FieldVector<Scalar, numComponents> &x,
Dune::FieldVector<Scalar, numComponents> &b,
Dune::FieldVector<Evaluation, numComponents> &x,
Dune::FieldVector<Evaluation, numComponents> &b,
int phaseIdx,
const Dune::FieldVector<Scalar, numComponents> &targetFug)
const Dune::FieldVector<Evaluation, numComponents> &targetFug)
{
typedef MathToolbox<Evaluation> Toolbox;
// store original composition and calculate relative error
Dune::FieldVector<Scalar, numComponents> origComp;
Dune::FieldVector<Evaluation, numComponents> origComp;
Scalar relError = 0;
Scalar sumDelta = 0;
Scalar sumx = 0;
Evaluation sumDelta = Toolbox::createConstant(0.0);
Evaluation sumx = Toolbox::createConstant(0.0);
for (int i = 0; i < numComponents; ++i) {
origComp[i] = fluidState.moleFraction(phaseIdx, i);
relError = std::max(relError, std::abs(x[i]));
relError = std::max(relError, std::abs(Toolbox::value(x[i])));
sumx += std::abs(fluidState.moleFraction(phaseIdx, i));
sumDelta += std::abs(x[i]);
sumx += Toolbox::abs(fluidState.moleFraction(phaseIdx, i));
sumDelta += Toolbox::abs(x[i]);
}
#if 1
// chop update to at most 20% change in composition
const Scalar maxDelta = 0.2;
if (sumDelta > maxDelta)
x /= (sumDelta/maxDelta);
#endif
//Scalar curDefect = calculateDefect_(fluidState, phaseIdx, targetFug);
//Scalar nextDefect;
//Scalar sumMoleFrac = 0.0;
//ComponentVector newB(1e100);
//for (int numTries = 0; numTries < 1; ++numTries) {
// change composition
for (int i = 0; i < numComponents; ++i) {
Scalar newx = origComp[i] - x[i];
#if 1
// only allow negative mole fractions if the target fugacity is negative
if (targetFug[i] > 0)
newx = std::max(0.0, newx);
// only allow positive mole fractions if the target fugacity is positive
else if (targetFug[i] < 0)
newx = std::min(0.0, newx);
// if the target fugacity is zero, the mole fraction must also be zero
else
newx = 0;
#endif
fluidState.setMoleFraction(phaseIdx, i, newx);
//sumMoleFrac += std::abs(newx);
}
// change composition
for (int i = 0; i < numComponents; ++i) {
Evaluation newx = origComp[i] - x[i];
// only allow negative mole fractions if the target fugacity is negative
if (targetFug[i] > 0)
newx = Toolbox::max(0.0, newx);
// only allow positive mole fractions if the target fugacity is positive
else if (targetFug[i] < 0)
newx = Toolbox::min(0.0, newx);
// if the target fugacity is zero, the mole fraction must also be zero
else
newx = 0;
paramCache.updateComposition(fluidState, phaseIdx);
/*
// if the sum of the mole fractions gets 0, we take the
// original composition divided by 100
if (sumMoleFrac < 1e-10) {
for (int i = 0; i < numComponents; ++i) {
fluidState.setMoleFraction(phaseIdx, i, origComp[i]/100);
}
return relError;
}
*/
/*
// calculate new residual
for (int i = 0; i < numComponents; ++i) {
Scalar phi = FluidSystem::computeFugacityCoeff(fluidState,
phaseIdx,
i);
fluidState.setFugacityCoeff(phaseIdx, i, phi);
}
nextDefect = calculateDefect_(fluidState, phaseIdx, targetFug);
//std::cout << "try delta: " << x << "\n";
//std::cout << "defect: old=" << curDefect << " new=" << nextDefect << "\n";
if (nextDefect <= curDefect)
break;
// divide delta vector
x /= 2;
fluidState.setMoleFraction(phaseIdx, i, newx);
}
*/
paramCache.updateComposition(fluidState, phaseIdx);
return relError;
}

13
opm/material/constraintsolvers/ComputeFromReferencePhase.hpp
View File

@ -58,13 +58,13 @@ namespace Opm {
* - 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
*/
template <class Scalar, class FluidSystem>
template <class Scalar, class FluidSystem, class Evaluation = Scalar>
class ComputeFromReferencePhase
{
enum { numPhases = FluidSystem::numPhases };
enum { numComponents = FluidSystem::numComponents };
typedef Opm::CompositionFromFugacities<Scalar, FluidSystem> CompositionFromFugacities;
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
typedef Opm::CompositionFromFugacities<Scalar, FluidSystem, Evaluation> CompositionFromFugacities;
typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
public:
/*!
@ -108,6 +108,7 @@ public:
bool setViscosity,
bool setEnthalpy)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
// compute the density and enthalpy of the
// reference phase
@ -145,8 +146,10 @@ public:
continue; // reference phase is already calculated
ComponentVector fugVec;
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
fugVec[compIdx] = fluidState.fugacity(refPhaseIdx, compIdx);
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
const auto& fug = fluidState.fugacity(refPhaseIdx, compIdx);
fugVec[compIdx] = FsToolbox::template toLhs<Evaluation>(fug);
}
CompositionFromFugacities::solve(fluidState, paramCache, phaseIdx, fugVec);

45
opm/material/constraintsolvers/MiscibleMultiPhaseComposition.hpp
View File

@ -25,6 +25,8 @@
#ifndef OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HPP
#define OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HPP
#include <opm/material/common/MathToolbox.hpp>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
@ -114,12 +116,14 @@ private:
* - if the setViscosity parameter is true, also dynamic viscosities of *all* phases
* - if the setInternalEnergy parameter is true, also specific enthalpies and internal energies of *all* phases
*/
template <class Scalar, class FluidSystem>
template <class Scalar, class FluidSystem, class Evaluation = Scalar>
class MiscibleMultiPhaseComposition
{
static const int numPhases = FluidSystem::numPhases;
static const int numComponents = FluidSystem::numComponents;
typedef MathToolbox<Evaluation> Toolbox;
static_assert(numPhases <= numComponents,
"This solver requires that the number fluid phases is smaller or equal "
"to the number of components");
@ -153,11 +157,15 @@ public:
static void solve(FluidState &fluidState,
ParameterCache &paramCache,
int phasePresence,
const MMPCAuxConstraint<Scalar> *auxConstraints,
const MMPCAuxConstraint<Evaluation> *auxConstraints,
int numAuxConstraints,
bool setViscosity,
bool setInternalEnergy)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
static_assert(std::is_same<typename FluidState::Scalar, Evaluation>::value,
"The scalar type of the fluid state must be 'Evaluation'");
#ifndef NDEBUG
// currently this solver can only handle fluid systems which
// assume ideal mixtures of all fluids. TODO: relax this
@ -176,7 +184,8 @@ public:
// coefficients of the components cannot depend on
// composition, i.e. the parameters in the cache are valid
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
Scalar fugCoeff = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
Evaluation fugCoeff = FsToolbox::template toLhs<Evaluation>(
FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
fluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
}
}
@ -184,25 +193,25 @@ public:
// create the linear system of equations which defines the
// mole fractions
static const int numEq = numComponents*numPhases;
Dune::FieldMatrix<Scalar, numEq, numEq> M(0.0);
Dune::FieldVector<Scalar, numEq> x(0.0);
Dune::FieldVector<Scalar, numEq> b(0.0);
Dune::FieldMatrix<Evaluation, numEq, numEq> M(Toolbox::createConstant(0.0));
Dune::FieldVector<Evaluation, numEq> x(Toolbox::createConstant(0.0));
Dune::FieldVector<Evaluation, numEq> b(Toolbox::createConstant(0.0));
// assemble the equations expressing the fact that the
// fugacities of each component are equal in all phases
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
Scalar entryCol1 =
const Evaluation& entryCol1 =
fluidState.fugacityCoefficient(/*phaseIdx=*/0, compIdx)
* fluidState.pressure(/*phaseIdx=*/0);
*fluidState.pressure(/*phaseIdx=*/0);
int col1Idx = compIdx;
for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
int rowIdx = (phaseIdx - 1)*numComponents + compIdx;
int col2Idx = phaseIdx*numComponents + compIdx;
Scalar entryCol2 =
const Evaluation& entryCol2 =
fluidState.fugacityCoefficient(phaseIdx, compIdx)
* fluidState.pressure(phaseIdx);
*fluidState.pressure(phaseIdx);
M[rowIdx][col1Idx] = entryCol1;
M[rowIdx][col2Idx] = -entryCol2;
@ -220,11 +229,11 @@ public:
int rowIdx = numComponents*(numPhases - 1) + presentPhases;
presentPhases += 1;
b[rowIdx] = 1.0;
b[rowIdx] = Toolbox::createConstant(1.0);
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
int colIdx = phaseIdx*numComponents + compIdx;
M[rowIdx][colIdx] = 1.0;
M[rowIdx][colIdx] = Toolbox::createConstant(1.0);
}
}
@ -262,17 +271,17 @@ public:
}
paramCache.updateComposition(fluidState, phaseIdx);
Scalar value = FluidSystem::density(fluidState, paramCache, phaseIdx);
fluidState.setDensity(phaseIdx, value);
const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
fluidState.setDensity(phaseIdx, rho);
if (setViscosity) {
value = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
fluidState.setViscosity(phaseIdx, value);
const Evaluation& mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
fluidState.setViscosity(phaseIdx, mu);
}
if (setInternalEnergy) {
value = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
fluidState.setEnthalpy(phaseIdx, value);
const Evaluation& h = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
fluidState.setEnthalpy(phaseIdx, h);
}
}
}

190
opm/material/constraintsolvers/NcpFlash.hpp
View File

@ -28,6 +28,7 @@
#include <opm/material/fluidmatrixinteractions/NullMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/common/ErrorMacros.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <opm/material/common/Means.hpp>
@ -87,22 +88,17 @@ class NcpFlash
static const int numEq = numPhases*(numComponents + 1);
typedef Dune::FieldMatrix<Scalar, numEq, numEq> Matrix;
typedef Dune::FieldVector<Scalar, numEq> Vector;
public:
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
/*!
* \brief Guess initial values for all quantities.
*/
template <class FluidState>
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static void guessInitial(FluidState &fluidState,
ParameterCache &paramCache,
const ComponentVector &globalMolarities)
const Dune::FieldVector<Evaluation, numComponents>& globalMolarities)
{
// the sum of all molarities
Scalar sumMoles = 0;
Evaluation sumMoles = 0;
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
sumMoles += globalMolarities[compIdx];
@ -124,7 +120,8 @@ public:
paramCache.updateAll(fluidState);
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
Scalar phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
const typename FluidState::Scalar phi =
FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
}
}
@ -140,15 +137,19 @@ public:
static void solve(FluidState &fluidState,
ParameterCache &paramCache,
const typename MaterialLaw::Params &matParams,
const ComponentVector &globalMolarities,
const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
Scalar tolerance = 0.0)
{
typedef typename FluidState::Scalar Evaluation;
typedef Dune::FieldMatrix<Evaluation, numEq, numEq> Matrix;
typedef Dune::FieldVector<Evaluation, numEq> Vector;
Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-35);
if (tolerance <= 0.0) {
tolerance = std::min(1e-10,
Opm::geometricMean(Scalar(1.0),
std::numeric_limits<Scalar>::epsilon()));
tolerance = std::min<Scalar>(1e-10,
Opm::geometricMean(Scalar(1.0),
std::numeric_limits<Scalar>::epsilon()));
}
/////////////////////////
@ -257,7 +258,7 @@ public:
* This is a convenience method which assumes that the capillary pressure is
* zero...
*/
template <class FluidState>
template <class FluidState, class ComponentVector>
static void solve(FluidState &fluidState,
const ComponentVector &globalMolarities,
Scalar tolerance = 0.0)
@ -320,14 +321,20 @@ protected:
std::cout << "\n";
}
template <class MaterialLaw, class FluidState>
static Scalar linearize_(Matrix &J,
Vector &b,
FluidState &fluidState,
ParameterCache &paramCache,
const typename MaterialLaw::Params &matParams,
const ComponentVector &globalMolarities)
template <class MaterialLaw,
class FluidState,
class Matrix,
class Vector,
class ComponentVector>
static void linearize_(Matrix &J,
Vector &b,
FluidState &fluidState,
ParameterCache &paramCache,
const typename MaterialLaw::Params &matParams,
const ComponentVector &globalMolarities)
{
typedef typename FluidState::Scalar Evaluation;
FluidState origFluidState(fluidState);
ParameterCache origParamCache(paramCache);
@ -340,28 +347,6 @@ protected:
calculateDefect_(b, fluidState, fluidState, globalMolarities);
Valgrind::CheckDefined(b);
///////
// calculate the absolute error
///////
Scalar absError = 0;
#if 0
// fugacities are equal
int eqIdx = 0;
for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
absError = std::max(std::abs(b[eqIdx + compIdx]), absError);
eqIdx += numComponents;
// sum of concentrations are given
for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
absError = std::max(std::abs(b[eqIdx + compIdx]), absError);
eqIdx += numComponents;
// NCP model assumptions
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
absError = std::max(std::abs(b[eqIdx + phaseIdx]), absError);
eqIdx += numPhases;
#endif
///////
// assemble jacobian matrix
///////
@ -373,17 +358,16 @@ protected:
// forward differences
// deviate the mole fraction of the i-th component
Scalar x_i = getQuantity_(fluidState, pvIdx);
const Evaluation& x_i = getQuantity_(fluidState, pvIdx);
const Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e7/(quantityWeight_(fluidState, pvIdx));
setQuantity_<MaterialLaw>(fluidState, paramCache, matParams, pvIdx, x_i + eps);
assert(std::abs(getQuantity_(fluidState, pvIdx) - (x_i + eps))
<= std::max(1.0, std::abs(x_i))*std::numeric_limits<Scalar>::epsilon()*100);
// compute derivative of the defect
calculateDefect_(tmp, origFluidState, fluidState, globalMolarities);
tmp -= b;
tmp /= eps;
for (int i = 0; i < numEq; ++ i)
tmp[i] /= eps;
// store derivative in jacobian matrix
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
@ -396,16 +380,16 @@ protected:
// end forward differences
////////
}
return absError;
}
template <class FluidState>
template <class FluidState, class Vector, class ComponentVector>
static void calculateDefect_(Vector &b,
const FluidState &fluidStateEval,
const FluidState &fluidState,
const ComponentVector &globalMolarities)
{
typedef typename FluidState::Scalar Evaluation;
int eqIdx = 0;
// fugacity of any component must be equal in all phases
@ -420,9 +404,8 @@ protected:
assert(eqIdx == numComponents*(numPhases - 1));
// the fact saturations must sum up to 1 is included explicitly!
// capillary pressures are explicitly included!
// the fact saturations must sum up to 1 is included implicitly and also,
// capillary pressures are treated implicitly!
// global molarities are given
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
@ -437,10 +420,9 @@ protected:
++eqIdx;
}
// model assumptions (-> non-linear complementarity functions)
// must be adhered
// model assumptions (-> non-linear complementarity functions) must be adhered
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
Scalar sumMoleFracEval = 0.0;
Evaluation sumMoleFracEval = 0.0;
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
sumMoleFracEval += fluidStateEval.moleFraction(phaseIdx, compIdx);
@ -448,7 +430,7 @@ protected:
b[eqIdx] = fluidState.saturation(phaseIdx);
}
else {
Scalar sumMoleFrac = 0.0;
Evaluation sumMoleFrac = 0.0;
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
sumMoleFrac += fluidState.moleFraction(phaseIdx, compIdx);
b[eqIdx] = 1.0 - sumMoleFrac;
@ -458,36 +440,43 @@ protected:
}
}
template <class MaterialLaw, class FluidState>
template <class MaterialLaw, class FluidState, class Vector>
static Scalar update_(FluidState &fluidState,
ParameterCache &paramCache,
const typename MaterialLaw::Params &matParams,
const Vector &deltaX)
{
typedef typename FluidState::Scalar Evaluation;
typedef Opm::MathToolbox<Evaluation> Toolbox;
// make sure we don't swallow non-finite update vectors
assert(std::isfinite(deltaX.two_norm()));
#ifndef NDEBUG
assert(deltaX.dimension == numEq);
for (int i = 0; i < numEq; ++i)
assert(std::isfinite(Toolbox::value(deltaX[i])));
#endif
Scalar relError = 0;
for (int pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
Scalar tmp = getQuantity_(fluidState, pvIdx);
Scalar delta = deltaX[pvIdx];
const Evaluation& tmp = getQuantity_(fluidState, pvIdx);
Evaluation delta = deltaX[pvIdx];
relError = std::max(relError, std::abs(delta)*quantityWeight_(fluidState, pvIdx));
relError = std::max<Scalar>(relError,
std::abs(Toolbox::value(delta))
* quantityWeight_(fluidState, pvIdx));
if (isSaturationIdx_(pvIdx)) {
// dampen to at most 25% change in saturation per
// iteration
delta = std::min(0.25, std::max(-0.25, delta));
// dampen to at most 25% change in saturation per iteration
delta = Toolbox::min(0.25, Toolbox::max(-0.25, delta));
}
else if (isMoleFracIdx_(pvIdx)) {
// dampen to at most 20% change in mole fraction per
// iteration
delta = std::min(0.20, std::max(-0.20, delta));
// dampen to at most 20% change in mole fraction per iteration
delta = Toolbox::min(0.20, Toolbox::max(-0.20, delta));
}
else if (isPressureIdx_(pvIdx)) {
// dampen to at most 50% change in pressure per
// iteration
delta = std::min(0.5*fluidState.pressure(0), std::max(-0.5*fluidState.pressure(0), delta));
// dampen to at most 50% change in pressure per iteration
delta = Toolbox::min(0.5*fluidState.pressure(0),
Toolbox::max(-0.5*fluidState.pressure(0), delta));
}
setQuantityRaw_(fluidState, pvIdx, tmp - delta);
@ -534,9 +523,11 @@ protected:
ParameterCache &paramCache,
const typename MaterialLaw::Params &matParams)
{
typedef typename FluidState::Scalar Evaluation;
// calculate the saturation of the last phase as a function of
// the other saturations
Scalar sumSat = 0.0;
Evaluation sumSat = 0.0;
for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
sumSat += fluidState.saturation(phaseIdx);
fluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);
@ -544,7 +535,7 @@ protected:
// update the pressures using the material law (saturations
// and first pressure are already set because it is implicitly
// solved for.)
ComponentVector pC;
Dune::FieldVector<Scalar, numPhases> pC;
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
fluidState.setPressure(phaseIdx,
@ -556,11 +547,11 @@ protected:
// update all densities and fugacity coefficients
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
const Evaluation& 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);
const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
}
}
@ -577,7 +568,7 @@ protected:
// retrieves a quantity from the fluid state
template <class FluidState>
static Scalar getQuantity_(const FluidState &fluidState, int pvIdx)
static const typename FluidState::Scalar& getQuantity_(const FluidState &fluidState, int pvIdx)
{
assert(pvIdx < numEq);
@ -606,12 +597,14 @@ protected:
ParameterCache &paramCache,
const typename MaterialLaw::Params &matParams,
int pvIdx,
Scalar value)
const typename FluidState::Scalar& value)
{
typedef typename FluidState::Scalar Evaluation;
assert(0 <= pvIdx && pvIdx < numEq);
if (pvIdx < 1) { // <- first pressure
Scalar delta = value - fluidState.pressure(0);
Evaluation delta = value - fluidState.pressure(0);
// set all pressures. here we assume that the capillary
// pressure does not depend on absolute pressure.
@ -621,41 +614,46 @@ protected:
// update all densities and fugacity coefficients
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
Valgrind::CheckDefined(rho);
fluidState.setDensity(phaseIdx, rho);
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
Scalar phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
Valgrind::CheckDefined(phi);
fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
}
}
}
else if (pvIdx < numPhases) { // <- first M - 1 saturations
Scalar delta = value - fluidState.saturation(/*phaseIdx=*/pvIdx - 1);
const Evaluation& delta = value - fluidState.saturation(/*phaseIdx=*/pvIdx - 1);
Valgrind::CheckDefined(delta);
fluidState.setSaturation(/*phaseIdx=*/pvIdx - 1, value);
// set last saturation (-> minus the change of the
// satuation of the other phase)
// set last saturation (-> minus the change of the saturation of the other
// phase)
fluidState.setSaturation(/*phaseIdx=*/numPhases - 1,
fluidState.saturation(numPhases - 1) - delta);
// update all fluid pressures using the capillary pressure
// law
ComponentVector pC;
// update all fluid pressures using the capillary pressure law
Dune::FieldVector<Evaluation, numPhases> pC;
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
Valgrind::CheckDefined(pC);
for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
fluidState.setPressure(phaseIdx,
fluidState.pressure(0)
+ (pC[phaseIdx] - pC[0]));
fluidState.pressure(0)
+ (pC[phaseIdx] - pC[0]));
paramCache.updateAllPressures(fluidState);
// update all densities and fugacity coefficients
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
Valgrind::CheckDefined(rho);
fluidState.setDensity(phaseIdx, rho);
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
Scalar phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
Valgrind::CheckDefined(phi);
fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
}
}
@ -665,18 +663,21 @@ protected:
int phaseIdx = (pvIdx - numPhases)/numComponents;
int compIdx = (pvIdx - numPhases)%numComponents;
Valgrind::CheckDefined(value);
fluidState.setMoleFraction(phaseIdx, compIdx, value);
paramCache.updateSingleMoleFraction(fluidState, phaseIdx, compIdx);
// update the density of the phase
Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
Valgrind::CheckDefined(rho);
fluidState.setDensity(phaseIdx, rho);
// if the phase's fugacity coefficients are composition
// dependent, update them as well.
if (!FluidSystem::isIdealMixture(phaseIdx)) {
for (int fugCompIdx = 0; fugCompIdx < numComponents; ++fugCompIdx) {
Scalar phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, fugCompIdx);
const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, fugCompIdx);
Valgrind::CheckDefined(phi);
fluidState.setFugacityCoefficient(phaseIdx, fugCompIdx, phi);
}
}
@ -688,10 +689,13 @@ protected:
// set a quantity in the fluid state
template <class FluidState>
static void setQuantityRaw_(FluidState &fluidState, int pvIdx, Scalar value)
static void setQuantityRaw_(FluidState &fluidState,
int pvIdx,
const typename FluidState::Scalar& value)
{
assert(pvIdx < numEq);
Valgrind::CheckDefined(value);
// first pressure
if (pvIdx < 1) {
int phaseIdx = 0;