Commit 3d2871fd authored by Andreas Lauser's avatar Andreas Lauser
Browse files

3p3c: use central differences to approximate the partial derivatives

this seems to solve the convergence problems, but why??

also, use tabulated Iapws water for the H2O-Air-Mesitylene fluid
system. this increases performance by an order of magnitute.

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7581 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 1d60b678
......@@ -268,13 +268,12 @@ public:
// floating point values, the resolution is about 10^-16 and
// the base epsilon is thus approximately 10^-8.
static const Scalar baseEps
= Dumux::geometricMean<Scalar>(std::numeric_limits<Scalar>::epsilon(),
1.0);
= Dumux::geometricMean<Scalar>(std::numeric_limits<Scalar>::epsilon(), 1.0);
// the epsilon value used for the numeric differentiation is
// now scaled by the absolute value of the primary variable...
Scalar pv = this->curVolVars_[scvIdx].primaryVar(pvIdx);
return baseEps*(std::abs(pv) + 1);
return baseEps*(std::abs(pv) + 1.0);
}
protected:
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2011 by Holger Class *
* Copyright (C) 2011-2012 by Holger Class *
* Copyright (C) 2009-2010 by Andreas Lauser *
* Institute of Hydraulic Engineering *
* University of Stuttgart, Germany *
......@@ -32,6 +32,8 @@
#include <dumux/material/idealgas.hh>
#include <dumux/material/components/air.hh>
#include <dumux/material/components/h2o.hh>
#include <dumux/material/components/tabulatedcomponent.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/components/mesitylene.hh>
#include <dumux/material/components/tabulatedcomponent.hh>
......@@ -57,10 +59,14 @@ class H2OAirMesitylene
typedef H2OAirMesitylene<Scalar> ThisType;
typedef BaseFluidSystem<Scalar, ThisType> Base;
typedef Dumux::H2O<Scalar> IapwsH2O;
typedef Dumux::TabulatedComponent<Scalar, IapwsH2O> TabulatedH2O;
typedef Dumux::SimpleH2O<Scalar> SimpleH2O;
public:
typedef Dumux::H2O<Scalar> H2O;
typedef Dumux::Mesitylene<Scalar> NAPL;
typedef Dumux::Air<Scalar> Air;
typedef TabulatedH2O H2O;
static const int numPhases = 3;
static const int numComponents = 3;
......@@ -73,8 +79,46 @@ public:
static const int NAPLIdx = 1;
static const int airIdx = 2;
/*!
* \brief Initialize the fluid system's static parameters generically
*
* If a tabulated H2O component is used, we do our best to create
* tables that always work.
*/
static void init()
{ }
{
init(/*tempMin=*/273.15,
/*tempMax=*/623.15,
/*numTemp=*/100,
/*pMin=*/0.0,
/*pMax=*/20e6,
/*numP=*/200);
}
/*!
* \brief Initialize the fluid system's static parameters using
* problem specific temperature and pressure ranges
*
* \param tempMin The minimum temperature used for tabulation of water [K]
* \param tempMax The maximum temperature used for tabulation of water [K]
* \param nTemp The number of ticks on the temperature axis of the table of water
* \param pressMin The minimum pressure used for tabulation of water [Pa]
* \param pressMax The maximum pressure used for tabulation of water [Pa]
* \param nPress The number of ticks on the pressure axis of the table of water
*/
static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
Scalar pressMin, Scalar pressMax, unsigned nPress)
{
if (H2O::isTabulated) {
std::cout << "Initializing tables for the H2O fluid properties ("
<< nTemp*nPress
<< " entries).\n";
TabulatedH2O::init(tempMin, tempMax, nTemp,
pressMin, pressMax, nPress);
}
}
/*!
* \brief Return whether a phase is liquid
......
......@@ -66,6 +66,9 @@ SET_BOOL_PROP(InfiltrationProblem, NewtonWriteConvergence, false);
// Maximum tolerated relative error in the Newton method
SET_SCALAR_PROP(InfiltrationProblem, NewtonRelTolerance, 1e-8);
// -1 backward differences, 0: central differences, +1: forward differences
SET_INT_PROP(InfiltrationProblem, NumericDifferenceMethod, 0);
}
/*!
......@@ -86,7 +89,6 @@ class InfiltrationProblem : public ThreePThreeCProblem<TypeTag>
// copy some indices for convenience
typedef typename GET_PROP_TYPE(TypeTag, ThreePThreeCIndices) Indices;
enum {
pressureIdx = Indices::pressureIdx,
switch1Idx = Indices::switch1Idx,
switch2Idx = Indices::switch2Idx,
......@@ -124,7 +126,13 @@ public:
: ParentType(timeManager, gridView)
, eps_(1e-6)
{
FluidSystem::init();
temperature_ = 273.15 + 10.0; // -> 10 degrees Celsius
FluidSystem::init(/*tempMin=*/temperature_ - 1,
/*tempMax=*/temperature_ + 1,
/*nTemp=*/3,
/*pressMin=*/0.8*1e5,
/*pressMax=*/3*1e5,
/*nPress=*/200);
}
/*!
......@@ -153,7 +161,7 @@ public:
const FVElementGeometry &fvElemGeom,
int scvIdx) const
{
return (273.15 + 10.0); // -> Temperatur 10°C
return temperature_;
};
void sourceAtPos(PrimaryVariables &values,
......@@ -362,6 +370,7 @@ private:
return(Sw);
}
Scalar temperature_;
Scalar eps_;
};
} //end namespace
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment