From 3b7ffc22180a19ff600fd1eca49f53c1273024a7 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Wed, 2 Mar 2011 16:42:52 +0000
Subject: [PATCH] make the modified 2p2c(ni) models run and rename them to
new_2p2c(ni)
also, restore the 2p2c(ni) models to their previous state.in order to
compile the new_2p2c(ni) tests, you have to set symbolic links to the
devel repository for the following the direcories
dumux/material/fluidstates
dumux/material/new_fluidsystems
dumux/material/eos
dumux/material/constraintsolvers
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@5335 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
configure.ac | 4 +
dumux/boxmodels/2p2c/2p2cfluidstate.hh | 2 -
dumux/boxmodels/2p2c/2p2cfluxvariables.hh | 8 +-
dumux/boxmodels/2p2c/2p2clocalresidual.hh | 13 +-
dumux/boxmodels/2p2c/2p2cmodel.hh | 24 +-
dumux/boxmodels/2p2c/2p2cnewtoncontroller.hh | 2 +-
dumux/boxmodels/2p2c/2p2cvolumevariables.hh | 314 ++------
dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh | 10 +-
.../boxmodels/2p2cni/2p2cnivolumevariables.hh | 93 ++-
dumux/boxmodels/Makefile.am | 2 +-
.../new_2p2c/2p2cboundaryvariables.hh | 252 ++++++
dumux/boxmodels/new_2p2c/2p2cfluidstate.hh | 371 +++++++++
dumux/boxmodels/new_2p2c/2p2cfluxvariables.hh | 382 +++++++++
dumux/boxmodels/new_2p2c/2p2cindices.hh | 128 +++
dumux/boxmodels/new_2p2c/2p2clocalresidual.hh | 342 ++++++++
dumux/boxmodels/new_2p2c/2p2cmodel.hh | 760 ++++++++++++++++++
.../new_2p2c/2p2cnewtoncontroller.hh | 181 +++++
dumux/boxmodels/new_2p2c/2p2cproblem.hh | 129 +++
dumux/boxmodels/new_2p2c/2p2cproperties.hh | 68 ++
.../new_2p2c/2p2cpropertydefaults.hh | 165 ++++
.../boxmodels/new_2p2c/2p2cvolumevariables.hh | 456 +++++++++++
dumux/boxmodels/new_2p2c/Makefile.am | 4 +
.../new_2p2cni/2p2cniboundaryvariables.hh | 272 +++++++
.../new_2p2cni/2p2cnifluxvariables.hh | 131 +++
dumux/boxmodels/new_2p2cni/2p2cniindices.hh | 57 ++
.../new_2p2cni/2p2cnilocalresidual.hh | 183 +++++
dumux/boxmodels/new_2p2cni/2p2cnimodel.hh | 108 +++
dumux/boxmodels/new_2p2cni/2p2cniproblem.hh | 80 ++
.../boxmodels/new_2p2cni/2p2cniproperties.hh | 54 ++
.../new_2p2cni/2p2cnipropertydefaults.hh | 85 ++
.../new_2p2cni/2p2cnivolumevariables.hh | 129 +++
dumux/boxmodels/new_2p2cni/Makefile.am | 4 +
test/boxmodels/2p2cni/waterairproblem.hh | 11 +-
test/boxmodels/Makefile.am | 4 +-
test/boxmodels/new_2p2c/CMakeLists.txt | 16 +
test/boxmodels/new_2p2c/Makefile.am | 8 +
test/boxmodels/new_2p2c/grids/test_2p2c.dgf | 10 +
test/boxmodels/new_2p2c/injectionproblem.hh | 379 +++++++++
.../new_2p2c/injectionspatialparameters.hh | 268 ++++++
test/boxmodels/new_2p2c/test_2p2c.cc | 34 +
test/boxmodels/new_2p2cni/CMakeLists.txt | 16 +
test/boxmodels/new_2p2cni/Makefile.am | 8 +
.../new_2p2cni/grids/test_2p2cni.dgf | 15 +
test/boxmodels/new_2p2cni/test_2p2cni.cc | 37 +
test/boxmodels/new_2p2cni/waterairproblem.hh | 397 +++++++++
.../new_2p2cni/waterairspatialparameters.hh | 269 +++++++
46 files changed, 5951 insertions(+), 334 deletions(-)
create mode 100644 dumux/boxmodels/new_2p2c/2p2cboundaryvariables.hh
create mode 100644 dumux/boxmodels/new_2p2c/2p2cfluidstate.hh
create mode 100644 dumux/boxmodels/new_2p2c/2p2cfluxvariables.hh
create mode 100644 dumux/boxmodels/new_2p2c/2p2cindices.hh
create mode 100644 dumux/boxmodels/new_2p2c/2p2clocalresidual.hh
create mode 100644 dumux/boxmodels/new_2p2c/2p2cmodel.hh
create mode 100644 dumux/boxmodels/new_2p2c/2p2cnewtoncontroller.hh
create mode 100644 dumux/boxmodels/new_2p2c/2p2cproblem.hh
create mode 100644 dumux/boxmodels/new_2p2c/2p2cproperties.hh
create mode 100644 dumux/boxmodels/new_2p2c/2p2cpropertydefaults.hh
create mode 100644 dumux/boxmodels/new_2p2c/2p2cvolumevariables.hh
create mode 100644 dumux/boxmodels/new_2p2c/Makefile.am
create mode 100644 dumux/boxmodels/new_2p2cni/2p2cniboundaryvariables.hh
create mode 100644 dumux/boxmodels/new_2p2cni/2p2cnifluxvariables.hh
create mode 100644 dumux/boxmodels/new_2p2cni/2p2cniindices.hh
create mode 100644 dumux/boxmodels/new_2p2cni/2p2cnilocalresidual.hh
create mode 100644 dumux/boxmodels/new_2p2cni/2p2cnimodel.hh
create mode 100644 dumux/boxmodels/new_2p2cni/2p2cniproblem.hh
create mode 100644 dumux/boxmodels/new_2p2cni/2p2cniproperties.hh
create mode 100644 dumux/boxmodels/new_2p2cni/2p2cnipropertydefaults.hh
create mode 100644 dumux/boxmodels/new_2p2cni/2p2cnivolumevariables.hh
create mode 100644 dumux/boxmodels/new_2p2cni/Makefile.am
create mode 100644 test/boxmodels/new_2p2c/CMakeLists.txt
create mode 100644 test/boxmodels/new_2p2c/Makefile.am
create mode 100644 test/boxmodels/new_2p2c/grids/test_2p2c.dgf
create mode 100644 test/boxmodels/new_2p2c/injectionproblem.hh
create mode 100644 test/boxmodels/new_2p2c/injectionspatialparameters.hh
create mode 100644 test/boxmodels/new_2p2c/test_2p2c.cc
create mode 100644 test/boxmodels/new_2p2cni/CMakeLists.txt
create mode 100644 test/boxmodels/new_2p2cni/Makefile.am
create mode 100644 test/boxmodels/new_2p2cni/grids/test_2p2cni.dgf
create mode 100644 test/boxmodels/new_2p2cni/test_2p2cni.cc
create mode 100644 test/boxmodels/new_2p2cni/waterairproblem.hh
create mode 100644 test/boxmodels/new_2p2cni/waterairspatialparameters.hh
diff --git a/configure.ac b/configure.ac
index f449b1fdc3..c73c682c54 100644
--- a/configure.ac
+++ b/configure.ac
@@ -25,6 +25,8 @@ AC_CONFIG_FILES([dumux.pc
dumux/boxmodels/2p/Makefile
dumux/boxmodels/2p2c/Makefile
dumux/boxmodels/2p2cni/Makefile
+ dumux/boxmodels/new_2p2c/Makefile
+ dumux/boxmodels/new_2p2cni/Makefile
dumux/boxmodels/2pni/Makefile
dumux/boxmodels/common/Makefile
dumux/boxmodels/richards/Makefile
@@ -62,6 +64,8 @@ AC_CONFIG_FILES([dumux.pc
test/boxmodels/1p2c/Makefile
test/boxmodels/2p2c/Makefile
test/boxmodels/2p2cni/Makefile
+ test/boxmodels/new_2p2c/Makefile
+ test/boxmodels/new_2p2cni/Makefile
test/boxmodels/richards/Makefile
test/common/Makefile
test/common/pardiso/Makefile
diff --git a/dumux/boxmodels/2p2c/2p2cfluidstate.hh b/dumux/boxmodels/2p2c/2p2cfluidstate.hh
index 1233d106dc..5e689939b7 100644
--- a/dumux/boxmodels/2p2c/2p2cfluidstate.hh
+++ b/dumux/boxmodels/2p2c/2p2cfluidstate.hh
@@ -26,7 +26,6 @@
* \brief Calculates the phase state from the primary variables in the
* 2p2c model.
*/
-#if 0
#ifndef DUMUX_2P2C_PHASE_STATE_HH
#define DUMUX_2P2C_PHASE_STATE_HH
@@ -368,4 +367,3 @@ public:
} // end namepace
#endif
-#endif
diff --git a/dumux/boxmodels/2p2c/2p2cfluxvariables.hh b/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
index ccebd24833..828268f48e 100644
--- a/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
+++ b/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
@@ -243,7 +243,6 @@ private:
const Element &element,
const ElementVolumeVariables &elemDat)
{
-#if 0
const VolumeVariables &vDat_i = elemDat[face().i];
const VolumeVariables &vDat_j = elemDat[face().j];
@@ -260,6 +259,7 @@ private:
// calculate tortuosity at the nodes i and j needed
// for porous media diffusion coefficient
+
Scalar tau_i =
1.0/(vDat_i.porosity() * vDat_i.porosity()) *
pow(vDat_i.porosity() * vDat_i.saturation(phaseIdx), 7.0/3);
@@ -271,8 +271,8 @@ private:
// -> harmonic mean
porousDiffCoeff_[phaseIdx] = harmonicMean(vDat_i.porosity() * vDat_i.saturation(phaseIdx) * tau_i * vDat_i.diffCoeff(phaseIdx),
vDat_j.porosity() * vDat_j.saturation(phaseIdx) * tau_j * vDat_j.diffCoeff(phaseIdx));
+
}
-#endif
}
public:
@@ -312,13 +312,11 @@ public:
int downstreamIdx(int phaseIdx) const
{ return downstreamIdx_[phaseIdx]; }
-#if 0
/*!
* \brief The binary diffusion coefficient for each fluid phase.
*/
Scalar porousDiffCoeff(int phaseIdx) const
{ return porousDiffCoeff_[phaseIdx]; };
-#endif
/*!
* \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase at the integration
@@ -371,10 +369,8 @@ protected:
// local index of the downwind vertex for each phase
int downstreamIdx_[numPhases];
-/*
// the diffusion coefficient for the porous medium
Scalar porousDiffCoeff_[numPhases];
-*/
};
} // end namepace
diff --git a/dumux/boxmodels/2p2c/2p2clocalresidual.hh b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
index 73d7c658e9..937e8c635b 100644
--- a/dumux/boxmodels/2p2c/2p2clocalresidual.hh
+++ b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
@@ -74,6 +74,8 @@ protected:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
+ typedef TwoPTwoCFluidState<TypeTag> FluidState;
+
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
enum
@@ -205,9 +207,7 @@ public:
flux = 0;
asImp_()->computeAdvectiveFlux(flux, vars);
- Valgrind::CheckDefined(flux);
asImp_()->computeDiffusiveFlux(flux, vars);
- Valgrind::CheckDefined(flux);
}
/*!
@@ -248,13 +248,6 @@ public:
- mobilityUpwindAlpha) * (dn.density(phaseIdx)
* dn.mobility(phaseIdx) * dn.fluidState().massFrac(
phaseIdx, compIdx));
- Valgrind::CheckDefined(vars.KmvpNormal(phaseIdx));
- Valgrind::CheckDefined(up.density(phaseIdx));
- Valgrind::CheckDefined(up.mobility(phaseIdx));
- Valgrind::CheckDefined(up.fluidState().massFrac(phaseIdx, compIdx));
- Valgrind::CheckDefined(dn.density(phaseIdx));
- Valgrind::CheckDefined(dn.mobility(phaseIdx));
- Valgrind::CheckDefined(dn.fluidState().massFrac(phaseIdx, compIdx));
}
}
}
@@ -268,7 +261,6 @@ public:
*/
void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &vars) const
{
-#if 0
// add diffusive flux of gas component in liquid phase
Scalar tmp =
- vars.porousDiffCoeff(lPhaseIdx) *
@@ -284,7 +276,6 @@ public:
(vars.molarConcGrad(gPhaseIdx) * vars.face().normal);
flux[contiLEqIdx] += tmp * FluidSystem::molarMass(lCompIdx);
flux[contiGEqIdx] -= tmp * FluidSystem::molarMass(gCompIdx);
-#endif
}
/*!
diff --git a/dumux/boxmodels/2p2c/2p2cmodel.hh b/dumux/boxmodels/2p2c/2p2cmodel.hh
index c59755e211..899ee39df6 100644
--- a/dumux/boxmodels/2p2c/2p2cmodel.hh
+++ b/dumux/boxmodels/2p2c/2p2cmodel.hh
@@ -143,6 +143,8 @@ class TwoPTwoCModel: public BoxModel<TypeTag>
formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation))
};
+ typedef TwoPTwoCFluidState<TypeTag> FluidState;
+
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
@@ -220,6 +222,9 @@ public:
setSwitched_(false);
resetPhasePresence_();
+ /*this->localJacobian().updateStaticData(this->curSolFunction(),
+ this->prevSolFunction());
+ */
};
/*!
@@ -370,7 +375,7 @@ public:
{
(*massFrac[phaseIdx][compIdx])[globalIdx]
= volVars.fluidState().massFrac(phaseIdx,
- compIdx);
+ compIdx);
Valgrind::CheckDefined(
(*massFrac[phaseIdx][compIdx])[globalIdx][0]);
@@ -638,8 +643,7 @@ protected:
// perform variable switch at a vertex; Returns true if a
// variable switch was performed.
bool primaryVarSwitch_(SolutionVector &globalSol,
- const VolumeVariables &volVars,
- int globalIdx,
+ const VolumeVariables &volVars, int globalIdx,
const GlobalPosition &globalPos)
{
// evaluate primary variable switch
@@ -665,16 +669,14 @@ protected:
if (xll + xlg > xlMax)
{
// liquid phase appears
- /*
- std::cout << "Liquid phase appears at vertex " << globalIdx
+ std::cout << "liquid phase appears at vertex " << globalIdx
<< ", coordinates: " << globalPos << ", xll + xlg: "
<< xll + xlg << std::endl;
- */
newPhasePresence = bothPhases;
if (formulation == pgSl)
- globalSol[globalIdx][switchIdx] = 0.001;
+ globalSol[globalIdx][switchIdx] = 0.0;
else if (formulation == plSg)
- globalSol[globalIdx][switchIdx] = 0.999;
+ globalSol[globalIdx][switchIdx] = 1.0;
};
}
else if (phasePresence == lPhaseOnly)
@@ -695,11 +697,9 @@ protected:
if (xgl + xgg > xgMax)
{
// gas phase appears
- /*
std::cout << "gas phase appears at vertex " << globalIdx
<< ", coordinates: " << globalPos << ", x_gl + x_gg: "
<< xgl + xgg << std::endl;
- */
newPhasePresence = bothPhases;
if (formulation == pgSl)
globalSol[globalIdx][switchIdx] = 0.999;
@@ -717,11 +717,9 @@ protected:
{
wouldSwitch = true;
// gas phase disappears
- /*
std::cout << "Gas phase disappears at vertex " << globalIdx
<< ", coordinates: " << globalPos << ", Sg: "
<< volVars.saturation(gPhaseIdx) << std::endl;
- */
newPhasePresence = lPhaseOnly;
globalSol[globalIdx][switchIdx]
@@ -731,11 +729,9 @@ protected:
{
wouldSwitch = true;
// liquid phase disappears
- /*
std::cout << "Liquid phase disappears at vertex " << globalIdx
<< ", coordinates: " << globalPos << ", Sl: "
<< volVars.saturation(lPhaseIdx) << std::endl;
- */
newPhasePresence = gPhaseOnly;
globalSol[globalIdx][switchIdx]
diff --git a/dumux/boxmodels/2p2c/2p2cnewtoncontroller.hh b/dumux/boxmodels/2p2c/2p2cnewtoncontroller.hh
index 1b4ad9376f..0001582c22 100644
--- a/dumux/boxmodels/2p2c/2p2cnewtoncontroller.hh
+++ b/dumux/boxmodels/2p2c/2p2cnewtoncontroller.hh
@@ -68,7 +68,7 @@ public:
TwoPTwoCNewtonController()
{
this->setRelTolerance(1e-7);
- this->setTargetSteps(10);
+ this->setTargetSteps(9);
this->setMaxSteps(18);
};
diff --git a/dumux/boxmodels/2p2c/2p2cvolumevariables.hh b/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
index 04734d91a2..004226f89a 100644
--- a/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
+++ b/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
@@ -33,15 +33,11 @@
#include <dumux/common/math.hh>
#include <dune/common/collectivecommunication.hh>
-
-#include <dumux/material/fluidstates/equilibriumfluidstate.hh>
-#include <dumux/material/constraintsolvers/compositionfromfugacities.hh>
-
#include <vector>
#include <iostream>
#include "2p2cproperties.hh"
-#include "2p2cindices.hh"
+#include "2p2cfluidstate.hh"
namespace Dumux
{
@@ -73,34 +69,17 @@ class TwoPTwoCVolumeVariables : public BoxVolumeVariables<TypeTag>
numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation)),
- // component indices
lCompIdx = Indices::lCompIdx,
gCompIdx = Indices::gCompIdx,
- // phase indices
lPhaseIdx = Indices::lPhaseIdx,
- gPhaseIdx = Indices::gPhaseIdx,
-
- // primary variable indices
- pressureIdx = Indices::pressureIdx,
- switchIdx = Indices::switchIdx,
-
- // phase states
- lPhaseOnly = Indices::lPhaseOnly,
- gPhaseOnly = Indices::gPhaseOnly,
- bothPhases = Indices::bothPhases,
-
- // formulations
- plSg = TwoPTwoCFormulation::plSg,
- pgSl = TwoPTwoCFormulation::pgSl,
+ gPhaseIdx = Indices::gPhaseIdx
};
typedef typename GridView::template Codim<0>::Entity Element;
+ typedef TwoPTwoCFluidState<TypeTag> FluidState;
public:
- //! The return type of the fluidState() method
- typedef Dumux::EquilibriumFluidState<Scalar, FluidSystem> FluidState;
-
/*!
* \brief Update all quantities for a given control volume.
*
@@ -125,227 +104,56 @@ public:
scvIdx,
isOldSol);
- typename FluidSystem::MutableParameters mutParams;
- typename FluidSystem::MutableParameters::FluidState &fs
- = mutParams.fluidState();
-
- int globalVertIdx = problem.model().vertexMapper().map(element, scvIdx, dim);
- int phasePresence = problem.model().phasePresence(globalVertIdx, isOldSol);
-
- /////////////
- // set the phase saturations
- /////////////
- if (phasePresence == gPhaseOnly) {
- fs.setSaturation(lPhaseIdx, 0.0);
- fs.setSaturation(gPhaseIdx, 1.0);
- }
- else if (phasePresence == lPhaseOnly) {
- fs.setSaturation(lPhaseIdx, 1.0);
- fs.setSaturation(gPhaseIdx, 0.0);
- }
- else if (phasePresence == bothPhases) {
- Scalar Sg;
- if (formulation == plSg)
- Sg = priVars[switchIdx];
- else if (formulation == pgSl)
- Sg = 1.0 - priVars[switchIdx];
-
- fs.setSaturation(lPhaseIdx, 1 - Sg);
- fs.setSaturation(gPhaseIdx, Sg);
- }
-
- /////////////
- // set the phase temperatures
- /////////////
- // update the temperature part of the energy module
- Scalar T = asImp_().getTemperature(priVars,
- element,
- elemGeom,
- scvIdx,
- problem);
- Valgrind::CheckDefined(T);
- for (int i = 0; i < numPhases; ++i)
- fs.setTemperature(i, T);
-
- /////////////
- // set the phase pressures
- /////////////
+ asImp().updateTemperature_(priVars,
+ element,
+ elemGeom,
+ scvIdx,
+ problem);
// capillary pressure parameters
const MaterialLawParams &materialParams =
problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
- Scalar pC = MaterialLaw::pC(materialParams, fs.saturation(lPhaseIdx));
- if (formulation == plSg) {
- fs.setPressure(lPhaseIdx, priVars[pressureIdx]);
- fs.setPressure(gPhaseIdx, priVars[pressureIdx] + pC);
- }
- else if (formulation == pgSl){
- fs.setPressure(lPhaseIdx, priVars[pressureIdx] - pC);
- fs.setPressure(gPhaseIdx, priVars[pressureIdx]);
- }
- Valgrind::CheckDefined(fs.pressure(lPhaseIdx));
- Valgrind::CheckDefined(fs.pressure(gPhaseIdx));
-
- // update the mutable parameters for the pure components
- mutParams.updatePure(lPhaseIdx);
- mutParams.updatePure(gPhaseIdx);
-
- /////////////
- // set the phase compositions.
- /////////////
- if (phasePresence == gPhaseOnly) {
- // mass fractions
- Scalar Xg1 = priVars[switchIdx];
- Scalar Xg2 = 1 - Xg1;
-
- // molar masses
- Scalar M1 = FluidSystem::molarMass(lCompIdx);
- Scalar M2 = FluidSystem::molarMass(gCompIdx);
-
- // convert mass to mole fractions
- Scalar meanM = M1*M2/(M2 + Xg2*(M1 - M2));
- fs.setMoleFrac(gPhaseIdx, lCompIdx, Xg1 * M1/meanM);
- fs.setMoleFrac(gPhaseIdx, gCompIdx, Xg2 * M2/meanM);
- mutParams.updateMix(gPhaseIdx);
-
- // calculate component fugacities in gas phase
- Scalar fug1 = FluidSystem::computeFugacity(mutParams, gPhaseIdx, lCompIdx);
- Scalar fug2 = FluidSystem::computeFugacity(mutParams, gPhaseIdx, gCompIdx);
- fs.setFugacity(gPhaseIdx, lCompIdx, fug1);
- fs.setFugacity(gPhaseIdx, gCompIdx, fug2);
-
- // use same fugacities in liquid phase
- fs.setFugacity(lPhaseIdx, lCompIdx, fug1);
- fs.setFugacity(lPhaseIdx, gCompIdx, fug2);
-
- // initial guess of liquid composition
- fs.setMoleFrac(lPhaseIdx, lCompIdx, 0.98);
- fs.setMoleFrac(lPhaseIdx, gCompIdx, 0.02);
-
- // calculate liquid composition from fugacities
- typedef Dumux::CompositionFromFugacities<Scalar, FluidSystem> CompFromFug;
- CompFromFug::run(mutParams, lPhaseIdx);
-
- Valgrind::CheckDefined(fs.moleFrac(gPhaseIdx, lCompIdx));
- Valgrind::CheckDefined(fs.moleFrac(gPhaseIdx, gCompIdx));
- Valgrind::CheckDefined(fs.moleFrac(lPhaseIdx, lCompIdx));
- Valgrind::CheckDefined(fs.moleFrac(lPhaseIdx, gCompIdx));
-
- // calculate molar volume of gas phase
- Scalar Vmg = FluidSystem::computeMolarVolume(mutParams, gPhaseIdx);
- fs.setMolarVolume(gPhaseIdx, Vmg);
- }
- else if (phasePresence == lPhaseOnly) {
- // mass fractions
- Scalar Xl2 = priVars[switchIdx];
- Scalar Xl1 = 1 - Xl2;
-
- // molar masses
- Scalar M1 = FluidSystem::molarMass(lCompIdx);
- Scalar M2 = FluidSystem::molarMass(gCompIdx);
-
- // convert mass to mole fractions
- Scalar meanM = M1*M2/(M2 + Xl2*(M1 - M2));
- fs.setMoleFrac(lPhaseIdx, lCompIdx, Xl1 * M1/meanM);
- fs.setMoleFrac(lPhaseIdx, gCompIdx, Xl2 * M2/meanM);
- mutParams.updateMix(lPhaseIdx);
+
+ int globalVertIdx = problem.model().dofMapper().map(element, scvIdx, dim);
+ int phasePresence = problem.model().phasePresence(globalVertIdx, isOldSol);
+
+ // calculate phase state
+ fluidState_.update(priVars, materialParams, temperature(), phasePresence);
+ Valgrind::CheckDefined(fluidState_);
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ // Mobilities
+ const Scalar mu =
+ FluidSystem::phaseViscosity(phaseIdx,
+ fluidState().temperature(),
+ fluidState().phasePressure(lPhaseIdx),
+ fluidState());
+ Scalar kr;
+ if (phaseIdx == lPhaseIdx)
+ kr = MaterialLaw::krw(materialParams, saturation(lPhaseIdx));
+ else // ATTENTION: krn requires the liquid saturation
+ // as parameter!
+ kr = MaterialLaw::krn(materialParams, saturation(lPhaseIdx));
+ mobility_[phaseIdx] = kr / mu;
+ Valgrind::CheckDefined(mobility_[phaseIdx]);
- // calculate component fugacities in liquid phase
- Scalar fug1 = FluidSystem::computeFugacity(mutParams, lPhaseIdx, lCompIdx);
- Scalar fug2 = FluidSystem::computeFugacity(mutParams, lPhaseIdx, gCompIdx);
- fs.setFugacity(lPhaseIdx, lCompIdx, fug1);
- fs.setFugacity(lPhaseIdx, gCompIdx, fug2);
-
- // use same fugacities in gas phase
- fs.setFugacity(gPhaseIdx, lCompIdx, fug1);
- fs.setFugacity(gPhaseIdx, gCompIdx, fug2);
-
- // initial guess of gas composition
- fs.setMoleFrac(gPhaseIdx, lCompIdx, 0.05);
- fs.setMoleFrac(gPhaseIdx, gCompIdx, 0.95);
-
- // calculate liquid composition from fugacities
- typedef Dumux::CompositionFromFugacities<Scalar, FluidSystem> CompFromFug;
- CompFromFug::run(mutParams, gPhaseIdx);
- Valgrind::CheckDefined(fs.moleFrac(gPhaseIdx, lCompIdx));
- Valgrind::CheckDefined(fs.moleFrac(gPhaseIdx, gCompIdx));
- Valgrind::CheckDefined(fs.moleFrac(lPhaseIdx, lCompIdx));
- Valgrind::CheckDefined(fs.moleFrac(lPhaseIdx, gCompIdx));
-
- // calculate molar volume of liquid phase
- Scalar Vml = FluidSystem::computeMolarVolume(mutParams, lPhaseIdx);
- fs.setMolarVolume(lPhaseIdx, Vml);
- }
- else if (phasePresence == bothPhases) {
- // HACK: assume both phases to be an ideal mixture,
- // i.e. the fugacity coefficents do not depend on the
- // composition
- Scalar phi_l1 = FluidSystem::computeFugacityCoeff(mutParams, lPhaseIdx, lCompIdx);
- Scalar phi_l2 = FluidSystem::computeFugacityCoeff(mutParams, lPhaseIdx, gCompIdx);
- Scalar phi_g1 = FluidSystem::computeFugacityCoeff(mutParams, gPhaseIdx, lCompIdx);
- Scalar phi_g2 = FluidSystem::computeFugacityCoeff(mutParams, gPhaseIdx, gCompIdx);
- Scalar pl = fs.pressure(lPhaseIdx);
- Scalar pg = fs.pressure(gPhaseIdx);
- Valgrind::CheckDefined(phi_l1);
- Valgrind::CheckDefined(phi_l2);
- Valgrind::CheckDefined(phi_g1);
- Valgrind::CheckDefined(phi_g2);
-
- Scalar xg2 = (phi_g2/phi_l1 - pl/pg) / (phi_g1/phi_l1 - phi_g2/phi_l2);
- Scalar xg1 = 1 - xg2;
- Scalar xl2 = (xg2*pg*phi_g2)/(pl*phi_l2);
- Scalar xl1 = 1 - xl2;
-
- fs.setMoleFrac(lPhaseIdx, lCompIdx, xl1);
- fs.setMoleFrac(lPhaseIdx, gCompIdx, xl2);
- fs.setMoleFrac(gPhaseIdx, lCompIdx, xg1);
- fs.setMoleFrac(gPhaseIdx, gCompIdx, xg2);
-
- mutParams.updateMix(lPhaseIdx);
- mutParams.updateMix(gPhaseIdx);
-
- Scalar Vml = FluidSystem::computeMolarVolume(mutParams, lPhaseIdx);
- Scalar Vmg = FluidSystem::computeMolarVolume(mutParams, gPhaseIdx);
- fs.setMolarVolume(lPhaseIdx, Vml);
- fs.setMolarVolume(gPhaseIdx, Vmg);
+ // binary diffusion coefficents
+ diffCoeff_[phaseIdx] =
+ FluidSystem::diffCoeff(phaseIdx,
+ lCompIdx,
+ gCompIdx,
+ fluidState_.temperature(),
+ fluidState_.phasePressure(phaseIdx),
+ fluidState_);
+ Valgrind::CheckDefined(diffCoeff_[phaseIdx]);
}
-
- // Mobilities
- Scalar muL = FluidSystem::computeViscosity(mutParams, lPhaseIdx);
- Scalar krL = MaterialLaw::krw(materialParams, fs.saturation(lPhaseIdx));
- mobility_[lPhaseIdx] = krL / muL;
- Valgrind::CheckDefined(mobility_[lPhaseIdx]);
-
- // ATTENTION: krn requires the liquid saturation as parameter!
- Scalar muG = FluidSystem::computeViscosity(mutParams, gPhaseIdx);
- Scalar krG = MaterialLaw::krn(materialParams, fs.saturation(lPhaseIdx));
- mobility_[gPhaseIdx] = krG / muG;
- Valgrind::CheckDefined(mobility_[gPhaseIdx]);
-
-#if 0
- // binary diffusion coefficents
- diffCoeff_[phaseIdx] =
- FluidSystem::diffCoeff(phaseIdx,
- lCompIdx,
- gCompIdx,
- fluidState_.temperature(),
- fluidState_.phasePressure(phaseIdx),
- fluidState_);
- Valgrind::CheckDefined(diffCoeff_[phaseIdx]);
-#endif
-
// porosity
porosity_ = problem.spatialParameters().porosity(element,
elemGeom,
scvIdx);
Valgrind::CheckDefined(porosity_);
-
- asImp_().updateEnergy(mutParams, priVars, element, elemGeom, scvIdx, problem);
-
- // assign the equilibrium fluid state from the generic one
- fluidState_.assign(fs);
- }
+ }
/*!
* \brief Returns the phase state for the control-volume.
@@ -378,7 +186,7 @@ public:
* \param phaseIdx The phase index
*/
Scalar molarDensity(int phaseIdx) const
- { return fluidState_.molarDensity(phaseIdx); }
+ { return fluidState_.density(phaseIdx) / fluidState_.averageMolarMass(phaseIdx); }
/*!
* \brief Returns the effective pressure of a given phase within
@@ -387,7 +195,7 @@ public:
* \param phaseIdx The phase index
*/
Scalar pressure(int phaseIdx) const
- { return fluidState_.pressure(phaseIdx); }
+ { return fluidState_.phasePressure(phaseIdx); }
/*!
* \brief Returns temperature inside the sub-control volume.
@@ -397,7 +205,7 @@ public:
* identical.
*/
Scalar temperature() const
- { return fluidState_.temperature(); }
+ { return temperature_; }
/*!
* \brief Returns the effective mobility of a given phase within
@@ -414,7 +222,7 @@ public:
* \brief Returns the effective capillary pressure within the control volume.
*/
Scalar capillaryPressure() const
- { return fluidState_.pressure(gPhaseIdx) - fluidState_.pressure(lPhaseIdx); }
+ { return fluidState_.capillaryPressure(); }
/*!
* \brief Returns the average porosity within the control volume.
@@ -422,43 +230,35 @@ public:
Scalar porosity() const
{ return porosity_; }
-#if 0
/*!
* \brief Returns the binary diffusion coefficients for a phase
*/
Scalar diffCoeff(int phaseIdx) const
{ return diffCoeff_[phaseIdx]; }
-#endif
+
protected:
- Scalar getTemperature(const PrimaryVariables &priVars,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int scvIdx,
- const Problem &problem)
+
+ void updateTemperature_(const PrimaryVariables &priVars,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ const Problem &problem)
{
- return problem.temperature(element, elemGeom, scvIdx);
+ temperature_ = problem.temperature(element, elemGeom, scvIdx);
}
-
- template <class MutableParameters>
- void updateEnergy(MutableParameters &mutParams,
- const PrimaryVariables &priVars,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int scvIdx,
- const Problem &problem)
- {};
+ Scalar temperature_; //!< Temperature within the control volume
Scalar porosity_; //!< Effective porosity within the control volume
Scalar mobility_[numPhases]; //!< Effective mobility within the control volume
-// Scalar diffCoeff_[numPhases]; //!< Binary diffusion coefficients for the phases
+ Scalar diffCoeff_[numPhases]; //!< Binary diffusion coefficients for the phases
FluidState fluidState_;
private:
- Implementation &asImp_()
+ Implementation &asImp()
{ return *static_cast<Implementation*>(this); }
- const Implementation &asImp_() const
+ const Implementation &asImp() const
{ return *static_cast<const Implementation*>(this); }
};
diff --git a/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh b/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh
index 50d1d0b698..e9abac3a8a 100644
--- a/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh
+++ b/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh
@@ -112,11 +112,13 @@ public:
// compute the energy storage
result[energyEqIdx] =
vertDat.porosity()*(vertDat.density(lPhaseIdx) *
- vertDat.fluidState().internalEnergy(lPhaseIdx) *
+ vertDat.internalEnergy(lPhaseIdx) *
+ //vertDat.enthalpy(lPhaseIdx) *
vertDat.saturation(lPhaseIdx)
+
vertDat.density(gPhaseIdx) *
- vertDat.fluidState().internalEnergy(gPhaseIdx) *
+ vertDat.internalEnergy(gPhaseIdx) *
+ //vertDat.enthalpy(gPhaseIdx) *
vertDat.saturation(gPhaseIdx))
+
vertDat.temperature()*vertDat.heatCapacity();
@@ -150,12 +152,12 @@ public:
mobilityUpwindAlpha * // upstream vertex
( up.density(phase) *
up.mobility(phase) *
- up.fluidState().enthalpy(phase))
+ up.enthalpy(phase))
+
(1-mobilityUpwindAlpha) * // downstream vertex
( dn.density(phase) *
dn.mobility(phase) *
- dn.fluidState().enthalpy(phase)) );
+ dn.enthalpy(phase)) );
}
}
diff --git a/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh b/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
index 2a12e12e13..a2187dc62e 100644
--- a/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
+++ b/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
@@ -70,19 +70,6 @@ class TwoPTwoCNIVolumeVariables : public TwoPTwoCVolumeVariables<TypeTag>
//! \endcond
public:
- /*!
- * \brief Update the temperature of the sub-control volume.
- */
- Scalar getTemperature(const PrimaryVariables &sol,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int scvIdx,
- const Problem &problem) const
- {
- // retrieve temperature from solution vector
- return sol[temperatureIdx];
- }
-
/*!
* \brief Update all quantities for a given control volume.
*
@@ -93,26 +80,55 @@ public:
* \param vertIdx The local index of the SCV (sub-control volume)
* \param isOldSol Evaluate function with solution of current or previous time step
*/
- template <class MutableParams>
- void updateEnergy(MutableParams &mutParams,
- const PrimaryVariables &sol,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int scvIdx,
- const Problem &problem)
+ void update(const PrimaryVariables &sol,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int vertIdx,
+ bool isOldSol)
{
- heatCapacity_ =
- problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
- Valgrind::CheckDefined(heatCapacity_);
+ // vertex update data for the mass balance
+ ParentType::update(sol,
+ problem,
+ element,
+ elemGeom,
+ vertIdx,
+ isOldSol);
// the internal energies and the enthalpies
- typename MutableParams::FluidState &fs = mutParams.fluidState();
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
- Scalar h =
- FluidSystem::computeEnthalpy(mutParams, phaseIdx);
- fs.setEnthalpy(phaseIdx, h);
+ enthalpy_[phaseIdx] =
+ FluidSystem::phaseEnthalpy(phaseIdx,
+ this->fluidState().temperature(),
+ this->fluidState().phasePressure(phaseIdx),
+ this->fluidState());
+ internalEnergy_[phaseIdx] =
+ FluidSystem::phaseInternalEnergy(phaseIdx,
+ this->fluidState().temperature(),
+ this->fluidState().phasePressure(phaseIdx),
+ this->fluidState());
}
- }
+ Valgrind::CheckDefined(internalEnergy_);
+ Valgrind::CheckDefined(enthalpy_);
+ };
+
+ /*!
+ * \brief Returns the total internal energy of a phase in the
+ * sub-control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar internalEnergy(int phaseIdx) const
+ { return internalEnergy_[phaseIdx]; };
+
+ /*!
+ * \brief Returns the total enthalpy of a phase in the sub-control
+ * volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar enthalpy(int phaseIdx) const
+ { return enthalpy_[phaseIdx]; };
/*!
* \brief Returns the total heat capacity \f$\mathrm{[J/(K*m^3]}\f$ of the rock matrix in
@@ -122,6 +138,27 @@ public:
{ return heatCapacity_; };
protected:
+ // this method gets called by the parent class. since this method
+ // is protected, we are friends with our parent..
+ friend class TwoPTwoCVolumeVariables<TypeTag>;
+ void updateTemperature_(const PrimaryVariables &sol,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ const Problem &problem)
+ {
+ // retrieve temperature from solution vector
+ this->temperature_ = sol[temperatureIdx];
+
+ heatCapacity_ =
+ problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
+
+ Valgrind::CheckDefined(this->temperature_);
+ Valgrind::CheckDefined(heatCapacity_);
+ }
+
+ Scalar internalEnergy_[numPhases];
+ Scalar enthalpy_[numPhases];
Scalar heatCapacity_;
};
diff --git a/dumux/boxmodels/Makefile.am b/dumux/boxmodels/Makefile.am
index 0fd56c7e68..9d175728d1 100644
--- a/dumux/boxmodels/Makefile.am
+++ b/dumux/boxmodels/Makefile.am
@@ -1,4 +1,4 @@
-SUBDIRS = common 1p 1p2c 2p 2p2c 2p2cni 2pni richards
+SUBDIRS = common 1p 1p2c 2p 2p2c 2p2cni new_2p2c new_2p2cni 2pni richards
boxmodelsdir = $(includedir)/dumux/boxmodels
diff --git a/dumux/boxmodels/new_2p2c/2p2cboundaryvariables.hh b/dumux/boxmodels/new_2p2c/2p2cboundaryvariables.hh
new file mode 100644
index 0000000000..58e4cb7d5c
--- /dev/null
+++ b/dumux/boxmodels/new_2p2c/2p2cboundaryvariables.hh
@@ -0,0 +1,252 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2010 by Katherina Baber, Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch, Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief This file contains the data which is required to calculate
+ * all fluxes of the fluid phase over the boundary of a finite volume.
+ *
+ * This means pressure and temperature gradients, phase densities at
+ * the integration point of the boundary, etc.
+ */
+#ifndef DUMUX_2P2C_BOUNDARY_VARIABLES_HH
+#define DUMUX_2P2C_BOUNDARY_VARIABLES_HH
+
+#include <dumux/common/math.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPTwoCModel
+ * \brief This template class contains the data which is required to
+ * calculate the fluxes of the fluid phases over the boundary of a
+ * finite volume for the 2p2c model.
+ *
+ * This means pressure and velocity gradients, phase density and viscosity at
+ * the integration point of the boundary, etc.
+ */
+template <class TypeTag>
+class TwoPTwoCBoundaryVariables
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+
+ enum {
+ dim = GridView::dimension,
+ };
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef Dune::FieldVector<Scalar, dim> VelocityVector;
+ typedef Dune::FieldVector<Scalar, dim> ScalarGradient;
+ typedef Dune::FieldMatrix<Scalar, dim, dim> VectorGradient;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename FVElementGeometry::SubControlVolume SCV;
+ typedef typename FVElementGeometry::BoundaryFace BoundaryFace;
+
+ enum {
+ lPhaseIdx = Indices::lPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+
+ lCompIdx = Indices::lCompIdx,
+ gCompIdx = Indices::gCompIdx,
+
+ numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases))
+ };
+
+public:
+ TwoPTwoCBoundaryVariables(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int boundaryFaceIdx,
+ const ElementVolumeVariables &elemDat,
+ int scvIdx)
+ : fvElemGeom_(elemGeom), scvIdx_(scvIdx)
+ {
+ boundaryFace_ = &fvElemGeom_.boundaryFace[boundaryFaceIdx];
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ densityAtIP_[phaseIdx] = Scalar(0);
+ pressureAtIP_[phaseIdx] = Scalar(0);
+ molarDensityAtIP_[phaseIdx] = Scalar(0);
+ potentialGrad_[phaseIdx] = Scalar(0);
+ concentrationGrad_[phaseIdx] = Scalar(0);
+ molarConcGrad_[phaseIdx] = Scalar(0);
+ }
+
+ calculateBoundaryValues_(problem, element, elemDat);
+ };
+
+ Scalar KmvpNormal(int phaseIdx) const
+ { return KmvpNormal_[phaseIdx]; }
+
+ /*!
+ * \brief The binary diffusion coefficient for each fluid phase.
+ */
+ Scalar porousDiffCoeff(int phaseIdx) const
+ { return porousDiffCoeff_[phaseIdx]; };
+
+ /*!
+ * \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase at the integration
+ * point.
+ */
+ Scalar densityAtIP(int phaseIdx) const
+ { return densityAtIP_[phaseIdx]; }
+
+ /*!
+ * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ of a phase at the integration
+ * point.
+ */
+ Scalar molarDensityAtIP(int phaseIdx) const
+ { return molarDensityAtIP_[phaseIdx]; }
+
+ /*!
+ * \brief The concentration gradient of a component in a phase.
+ */
+ const ScalarGradient &concentrationGrad(int phaseIdx) const
+ { return concentrationGrad_[phaseIdx]; };
+
+ /*!
+ * \brief The molar concentration gradient of a component in a phase.
+ */
+ const ScalarGradient &molarConcGrad(int phaseIdx) const
+ { return molarConcGrad_[phaseIdx]; };
+
+ const FVElementGeometry &fvElemGeom() const
+ { return fvElemGeom_; }
+
+ const BoundaryFace& boundaryFace() const
+ { return *boundaryFace_; }
+
+protected:
+ void calculateBoundaryValues_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemDat)
+ {
+ ScalarGradient tmp(0.0);
+
+ // calculate gradients and secondary variables at IPs of the boundary
+ for (int idx = 0;
+ idx < fvElemGeom_.numVertices;
+ idx++) // loop over adjacent vertices
+ {
+ // FE gradient at vertex idx
+ const ScalarGradient& feGrad = boundaryFace_->grad[idx];
+
+ // compute sum of pressure gradients for each phase
+ for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
+ {
+ // the pressure gradient
+ tmp = feGrad;
+ tmp *= elemDat[idx].pressure(phaseIdx);
+ potentialGrad_[phaseIdx] += tmp;
+ }
+
+ // the concentration gradient of the non-wetting
+ // component in the wetting phase
+ tmp = feGrad;
+ tmp *= elemDat[idx].fluidState().massFrac(lPhaseIdx, gCompIdx);
+ concentrationGrad_[lPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemDat[idx].fluidState().moleFrac(lPhaseIdx, gCompIdx);
+ molarConcGrad_[lPhaseIdx] += tmp;
+
+ // the concentration gradient of the wetting component
+ // in the non-wetting phase
+ tmp = feGrad;
+ tmp *= elemDat[idx].fluidState().massFrac(gPhaseIdx, lCompIdx);
+ concentrationGrad_[gPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemDat[idx].fluidState().moleFrac(gPhaseIdx, lCompIdx);
+ molarConcGrad_[gPhaseIdx] += tmp;
+
+ for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
+ {
+ pressureAtIP_[phaseIdx] += elemDat[idx].pressure(phaseIdx)*boundaryFace_->shapeValue[idx];
+ densityAtIP_[phaseIdx] += elemDat[idx].density(phaseIdx)*boundaryFace_->shapeValue[idx];
+ molarDensityAtIP_[phaseIdx] += elemDat[idx].molarDensity(phaseIdx)*boundaryFace_->shapeValue[idx];
+ }
+ }
+
+ for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
+ {
+ tmp = problem.gravity();
+ tmp *= densityAtIP_[phaseIdx];
+
+ potentialGrad_[phaseIdx] -= tmp;
+
+ Scalar k = problem.spatialParameters().intrinsicPermeability(element, fvElemGeom_, scvIdx_);
+ VectorGradient K(0);
+ K[0][0] = K[1][1] = k;
+ ScalarGradient Kmvp;
+ K.mv(potentialGrad_[phaseIdx], Kmvp);
+ KmvpNormal_[phaseIdx] = - (Kmvp*boundaryFace_->normal);
+
+ const VolumeVariables &vertDat = elemDat[scvIdx_];
+
+ if (vertDat.saturation(phaseIdx) <= 0)
+ porousDiffCoeff_[phaseIdx] = 0.0;
+ else
+ {
+ Scalar tau = 1.0/(vertDat.porosity()*vertDat.porosity())*
+ pow(vertDat.porosity()*vertDat.saturation(phaseIdx), 7.0/3);
+
+ porousDiffCoeff_[phaseIdx] = vertDat.porosity()*vertDat.saturation(phaseIdx)*tau*vertDat.diffCoeff(phaseIdx);
+ }
+ }
+ }
+
+ const FVElementGeometry &fvElemGeom_;
+ const BoundaryFace *boundaryFace_;
+
+ // gradients
+ ScalarGradient potentialGrad_[numPhases];
+ ScalarGradient concentrationGrad_[numPhases];
+ ScalarGradient molarConcGrad_[numPhases];
+
+ // density of each face at the integration point
+ Scalar pressureAtIP_[numPhases];
+ Scalar densityAtIP_[numPhases];
+ Scalar molarDensityAtIP_[numPhases];
+
+ // intrinsic permeability times pressure potential gradient
+ // projected on the face normal
+ Scalar KmvpNormal_[numPhases];
+
+ // the diffusion coefficient for the porous medium
+ Scalar porousDiffCoeff_[numPhases];
+
+ int scvIdx_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/boxmodels/new_2p2c/2p2cfluidstate.hh b/dumux/boxmodels/new_2p2c/2p2cfluidstate.hh
new file mode 100644
index 0000000000..1233d106dc
--- /dev/null
+++ b/dumux/boxmodels/new_2p2c/2p2cfluidstate.hh
@@ -0,0 +1,371 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2010 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Calculates the phase state from the primary variables in the
+ * 2p2c model.
+ */
+#if 0
+#ifndef DUMUX_2P2C_PHASE_STATE_HH
+#define DUMUX_2P2C_PHASE_STATE_HH
+
+#include "2p2cproperties.hh"
+#include "2p2cindices.hh"
+
+#include <dumux/material/fluidstate.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPTwoCModel
+ * \brief Calculates the phase state from the primary variables in the
+ * 2p2c model.
+ */
+template <class TypeTag>
+class TwoPTwoCFluidState : public FluidState<typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)),
+ TwoPTwoCFluidState<TypeTag> >
+{
+ typedef TwoPTwoCFluidState<TypeTag> ThisType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ enum {
+ lPhaseIdx = Indices::lPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+
+ lCompIdx = Indices::lCompIdx,
+ gCompIdx = Indices::gCompIdx,
+
+ switchIdx = Indices::switchIdx,
+ pressureIdx = Indices::pressureIdx,
+ };
+
+ // present phases
+ enum {
+ lPhaseOnly = Indices::lPhaseOnly,
+ gPhaseOnly = Indices::gPhaseOnly,
+ bothPhases = Indices::bothPhases
+ };
+
+ // formulations
+ enum {
+ formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation)),
+ plSg = TwoPTwoCFormulation::plSg,
+ pgSl = TwoPTwoCFormulation::pgSl
+ };
+
+public:
+ enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)) };
+ enum { numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)) };
+ enum { numSolvents = 1 };
+
+ /*!
+ * \brief Update the phase state from the primary variables.
+ *
+ * \param primaryVars The primary variables
+ * \param pcParams The parameters for the material law
+ * \param temperature The temperature
+ * \param phasePresence Stands either for nonwetting phase, wetting phase or both phases
+ */
+ void update(const PrimaryVariables &primaryVars,
+ const MaterialLawParams &pcParams,
+ Scalar temperature,
+ int phasePresence)
+ {
+ Valgrind::CheckDefined(primaryVars);
+
+ Valgrind::SetUndefined(concentration_);
+ Valgrind::SetUndefined(phaseConcentration_);
+ Valgrind::SetUndefined(avgMolarMass_);
+ Valgrind::SetUndefined(phasePressure_);
+ Valgrind::SetUndefined(temperature_);
+ Valgrind::SetUndefined(Sg_);
+
+ temperature_ = temperature;
+ Valgrind::CheckDefined(temperature_);
+
+ // extract non-wetting phase pressure
+ if (phasePresence == gPhaseOnly)
+ Sg_ = 1.0;
+ else if (phasePresence == lPhaseOnly) {
+ Sg_ = 0.0;
+ }
+ else if (phasePresence == bothPhases) {
+ if (formulation == plSg)
+ Sg_ = primaryVars[switchIdx];
+ else if (formulation == pgSl)
+ Sg_ = 1.0 - primaryVars[switchIdx];
+ else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
+ }
+ else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
+ Valgrind::CheckDefined(Sg_);
+
+ // calculate capillary pressure
+ Scalar pC = MaterialLaw::pC(pcParams, 1 - Sg_);
+
+ // extract the pressures
+ if (formulation == plSg) {
+ phasePressure_[lPhaseIdx] = primaryVars[pressureIdx];
+ phasePressure_[gPhaseIdx] = phasePressure_[lPhaseIdx] + pC;
+ }
+ else if (formulation == pgSl) {
+ phasePressure_[gPhaseIdx] = primaryVars[pressureIdx];
+ phasePressure_[lPhaseIdx] = phasePressure_[gPhaseIdx] - pC;
+ }
+ else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
+ Valgrind::CheckDefined(phasePressure_);
+
+ // now comes the tricky part: calculate phase composition
+ if (phasePresence == bothPhases) {
+ // both phases are present, phase composition results from
+ // the gas <-> liquid equilibrium.
+ FluidSystem::computeEquilibrium(*this, -1);
+ }
+ else if (phasePresence == gPhaseOnly) {
+ // only the gas phase is present, gas phase composition is
+ // stored explicitly.
+
+ // extract _mass_ (!) fractions in the gas phase, misuse
+ // the concentration_ array for temporary storage
+ concentration_[gPhaseIdx][lCompIdx] = primaryVars[switchIdx];
+ concentration_[gPhaseIdx][gCompIdx] = 1 - concentration_[gPhaseIdx][lCompIdx];
+
+ // calculate average molar mass of the gas phase
+ Scalar M1 = FluidSystem::molarMass(lCompIdx);
+ Scalar M2 = FluidSystem::molarMass(gCompIdx);
+ Scalar X2 = concentration_[gPhaseIdx][gCompIdx]; // mass fraction of solvent in gas
+ avgMolarMass_[gPhaseIdx] = M1*M2/(M2 + X2*(M1 - M2));
+
+ // convert mass to mole fractions
+ concentration_[gPhaseIdx][lCompIdx] *= avgMolarMass_[gPhaseIdx]/M1;
+ concentration_[gPhaseIdx][gCompIdx] *= avgMolarMass_[gPhaseIdx]/M2;
+
+ // convert to real concentrations
+ phaseConcentration_[gPhaseIdx] = 1.0;
+ Scalar rhog = FluidSystem::phaseDensity(gPhaseIdx,
+ temperature_,
+ phasePressure_[gPhaseIdx],
+ *this);
+ phaseConcentration_[gPhaseIdx] = rhog / avgMolarMass_[gPhaseIdx];
+ concentration_[gPhaseIdx][lCompIdx] *= phaseConcentration_[gPhaseIdx];
+ concentration_[gPhaseIdx][gCompIdx] *= phaseConcentration_[gPhaseIdx];
+
+ // tell the fluid system to calculate the composition of
+ // the remaining phases
+ FluidSystem::computeEquilibrium(*this, gPhaseIdx);
+ }
+ else if (phasePresence == lPhaseOnly) {
+ // only the gas phase is present, liquid phase composition is
+ // stored explicitly.
+
+ // extract _mass_ (!) fractions in the liquid phase, misuse
+ // the concentration_ array for temporary storage
+ concentration_[lPhaseIdx][gCompIdx] = primaryVars[switchIdx];
+ concentration_[lPhaseIdx][lCompIdx] = 1 - concentration_[lPhaseIdx][gCompIdx];
+
+ Scalar M1 = FluidSystem::molarMass(lCompIdx);
+ Scalar M2 = FluidSystem::molarMass(gCompIdx);
+ Scalar X2 = concentration_[lPhaseIdx][gCompIdx]; // mass fraction of solvent in gas
+ avgMolarMass_[lPhaseIdx] = M1*M2/(M2 + X2*(M1 - M2));
+
+ // convert mass to mole fractions
+ concentration_[lPhaseIdx][lCompIdx] *= avgMolarMass_[lPhaseIdx]/M1;
+ concentration_[lPhaseIdx][gCompIdx] *= avgMolarMass_[lPhaseIdx]/M2;
+
+ // convert to real concentrations
+ phaseConcentration_[lPhaseIdx] = 1.0;
+ Scalar rhol = FluidSystem::phaseDensity(lPhaseIdx,
+ temperature_,
+ phasePressure_[lPhaseIdx],
+ *this);
+ phaseConcentration_[lPhaseIdx] = rhol / avgMolarMass_[lPhaseIdx];
+ concentration_[lPhaseIdx][lCompIdx] *= phaseConcentration_[lPhaseIdx];
+ concentration_[lPhaseIdx][gCompIdx] *= phaseConcentration_[lPhaseIdx];
+
+ // tell the fluid system to calculate the composition of
+ // the remaining phases
+ FluidSystem::computeEquilibrium(*this, lPhaseIdx);
+ }
+
+ Valgrind::CheckDefined(concentration_);
+ Valgrind::CheckDefined(phaseConcentration_);
+ Valgrind::CheckDefined(avgMolarMass_);
+ Valgrind::CheckDefined(phasePressure_);
+ Valgrind::CheckDefined(temperature_);
+ Valgrind::CheckDefined(Sg_);
+ }
+
+public:
+ /*!
+ * \brief Retrieves the phase composition and pressure from a
+ * phase composition class.
+ *
+ * \param phaseIdx The phase index
+ * \param compo The index of the component
+ *
+ * This method is called by the fluid system's
+ * computeEquilibrium()
+ */
+ template <class PhaseCompo>
+ void assignPhase(int phaseIdx, const PhaseCompo &compo)
+ {
+ avgMolarMass_[phaseIdx] = compo.meanMolarMass();
+ phasePressure_[phaseIdx] = compo.pressure();
+
+ phaseConcentration_[phaseIdx] = compo.phaseConcentration();
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ concentration_[phaseIdx][compIdx] = compo.concentration(compIdx);
+ };
+
+ /*!
+ * \brief Returns the saturation of a phase.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar saturation(int phaseIdx) const
+ {
+ if (phaseIdx == lPhaseIdx)
+ return Scalar(1.0) - Sg_;
+ else
+ return Sg_;
+ };
+
+ /*!
+ * \brief Returns the mole fraction of a component in a fluid phase.
+ *
+ * \param phaseIdx The phase index
+ * \param compIdx The index of the component
+ */
+ Scalar moleFrac(int phaseIdx, int compIdx) const
+ {
+ return
+ concentration_[phaseIdx][compIdx]/
+ phaseConcentration_[phaseIdx];
+ }
+
+ /*!
+ * \brief Returns the molar density of a phase \f$\mathrm{[mol/m^3]}\f$.
+ *
+ * \param phaseIdx The phase index
+ *
+ * This is equivalent to the sum of all molar component concentrations.
+ */
+ Scalar phaseConcentration(int phaseIdx) const
+ { return phaseConcentration_[phaseIdx]; };
+
+ /*!
+ * \brief Returns the molar concentration of a component in a phase \f$\mathrm{[mol/m^3]}\f$.
+ *
+ * \param phaseIdx The phase index
+ * \param compIdx The index of the component
+ *
+ */
+ Scalar concentration(int phaseIdx, int compIdx) const
+ { return concentration_[phaseIdx][compIdx]; };
+
+ /*!
+ * \brief Returns the mass fraction of a component in a phase.
+ *
+ * \param phaseIdx The phase index
+ * \param compIdx The index of the component
+ *
+ */
+ Scalar massFrac(int phaseIdx, int compIdx) const
+ {
+ return
+ moleFrac(phaseIdx, compIdx) *
+ FluidSystem::molarMass(compIdx)/avgMolarMass_[phaseIdx];
+ }
+
+ /*!
+ * \brief Returns the mass density of a phase \f$\mathrm{[kg/m^3]}\f$.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar density(int phaseIdx) const
+ { return phaseConcentration_[phaseIdx]*avgMolarMass_[phaseIdx]; }
+
+ /*!
+ * \brief Returns mean molar mass of a phase \f$\mathrm{[kg/mol]}\f$.
+ *
+ * \param phaseIdx The phase index
+ *
+ * This is equivalent to the sum of all component molar masses
+ * weighted by their respective mole fraction.
+ */
+ Scalar averageMolarMass(int phaseIdx) const
+ { return avgMolarMass_[phaseIdx]; };
+
+ /*!
+ * \brief Returns the pressure of a fluid phase \f$\mathrm{[Pa]}\f$.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar phasePressure(int phaseIdx) const
+ { return phasePressure_[phaseIdx]; }
+
+ /*!
+ * \brief Return the fugacity of a component \f$\mathrm{[Pa]}\f$.
+ *
+ * \param compIdx The index of the component
+ */
+ Scalar fugacity(int compIdx) const
+ { return moleFrac(gPhaseIdx, compIdx)*phasePressure(gPhaseIdx); };
+
+ /*!
+ * \brief Returns the capillary pressure \f$\mathrm{[Pa]}\f$
+ */
+ Scalar capillaryPressure() const
+ { return phasePressure_[gPhaseIdx] - phasePressure_[lPhaseIdx]; }
+
+ /*!
+ * \brief Returns the temperature of the fluids \f$\mathrm{[K]}\f$.
+ *
+ * Note that we assume thermodynamic equilibrium, so all fluids
+ * and the rock matrix exhibit the same temperature.
+ */
+ Scalar temperature() const
+ { return temperature_; };
+
+public:
+ Scalar concentration_[numPhases][numComponents];
+ Scalar phaseConcentration_[numPhases];
+ Scalar avgMolarMass_[numPhases];
+ Scalar phasePressure_[numPhases];
+ Scalar temperature_;
+ Scalar Sg_;
+};
+
+} // end namepace
+
+#endif
+#endif
diff --git a/dumux/boxmodels/new_2p2c/2p2cfluxvariables.hh b/dumux/boxmodels/new_2p2c/2p2cfluxvariables.hh
new file mode 100644
index 0000000000..ccebd24833
--- /dev/null
+++ b/dumux/boxmodels/new_2p2c/2p2cfluxvariables.hh
@@ -0,0 +1,382 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief This file contains the data which is required to calculate
+ * all fluxes of components over a face of a finite volume for
+ * the two-phase, two-component model.
+ */
+/*!
+ * \ingroup TwoPTwoCModel
+ */
+#ifndef DUMUX_2P2C_FLUX_VARIABLES_HH
+#define DUMUX_2P2C_FLUX_VARIABLES_HH
+
+#include <dumux/common/math.hh>
+#include <dumux/common/spline.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \brief This template class contains the data which is required to
+ * calculate all fluxes of components over a face of a finite
+ * volume for the two-phase, two-component model.
+ *
+ * This means pressure and concentration gradients, phase densities at
+ * the integration point, etc.
+ */
+template <class TypeTag>
+class TwoPTwoCFluxVariables
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
+
+ typedef typename GridView::ctype CoordScalar;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+ numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases))
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+ typedef typename FVElementGeometry::SubControlVolume SCV;
+ typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
+
+ typedef Dune::FieldVector<CoordScalar, dimWorld> Vector;
+ typedef Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld> Tensor;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+ enum {
+ lPhaseIdx = Indices::lPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+
+ lCompIdx = Indices::lCompIdx,
+ gCompIdx = Indices::gCompIdx,
+ };
+
+public:
+ /*
+ * \brief The constructor
+ *
+ * \param problem The problem
+ * \param element The finite element
+ * \param elemGeom The finite-volume geometry in the box scheme
+ * \param faceIdx The local index of the SCV (sub-control-volume) face
+ * \param elemDat The volume variables of the current element
+ */
+ TwoPTwoCFluxVariables(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int faceIdx,
+ const ElementVolumeVariables &elemDat)
+ : fvElemGeom_(elemGeom)
+ {
+ scvfIdx_ = faceIdx;
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ densityAtIP_[phaseIdx] = Scalar(0);
+ molarDensityAtIP_[phaseIdx] = Scalar(0);
+ potentialGrad_[phaseIdx] = Scalar(0);
+ concentrationGrad_[phaseIdx] = Scalar(0);
+ molarConcGrad_[phaseIdx] = Scalar(0);
+ }
+
+ calculateGradients_(problem, element, elemDat);
+ calculateVelocities_(problem, element, elemDat);
+ calculateDiffCoeffPM_(problem, element, elemDat);
+ };
+
+private:
+ void calculateGradients_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemDat)
+ {
+ // calculate gradients
+ Vector tmp(0.0);
+ for (int idx = 0;
+ idx < fvElemGeom_.numVertices;
+ idx++) // loop over adjacent vertices
+ {
+ // FE gradient at vertex idx
+ const Vector &feGrad = face().grad[idx];
+
+ // compute sum of pressure gradients for each phase
+ for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
+ {
+ // the pressure gradient
+ tmp = feGrad;
+ tmp *= elemDat[idx].pressure(phaseIdx);
+ potentialGrad_[phaseIdx] += tmp;
+ }
+
+ // the concentration gradient of the non-wetting
+ // component in the wetting phase
+ tmp = feGrad;
+ tmp *= elemDat[idx].fluidState().massFrac(lPhaseIdx, gCompIdx);
+ concentrationGrad_[lPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemDat[idx].fluidState().moleFrac(lPhaseIdx, gCompIdx);
+ molarConcGrad_[lPhaseIdx] += tmp;
+
+// // the concentration gradient of the wetting component
+// // in the non-wetting phase
+ tmp = feGrad;
+ tmp *= elemDat[idx].fluidState().massFrac(gPhaseIdx, lCompIdx);
+ concentrationGrad_[gPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemDat[idx].fluidState().moleFrac(gPhaseIdx, lCompIdx);
+ molarConcGrad_[gPhaseIdx] += tmp;
+ }
+
+ // correct the pressure gradients by the hydrostatic
+ // pressure due to gravity
+ for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
+ {
+ int i = face().i;
+ int j = face().j;
+ Scalar fI = rhoFactor_(phaseIdx, i, elemDat);
+ Scalar fJ = rhoFactor_(phaseIdx, j, elemDat);
+ if (fI + fJ <= 0)
+ fI = fJ = 0.5; // doesn't matter because no phase is
+ // present in both cells!
+ densityAtIP_[phaseIdx] =
+ (fI*elemDat[i].density(phaseIdx) +
+ fJ*elemDat[j].density(phaseIdx))
+ /
+ (fI + fJ);
+ // phase density
+ molarDensityAtIP_[phaseIdx]
+ =
+ (fI*elemDat[i].molarDensity(phaseIdx) +
+ fJ*elemDat[j].molarDensity(phaseIdx))
+ /
+ (fI + fJ);
+
+ tmp = problem.gravity();
+ tmp *= densityAtIP_[phaseIdx];
+
+ potentialGrad_[phaseIdx] -= tmp;
+ }
+ }
+
+ Scalar rhoFactor_(int phaseIdx, int scvIdx, const ElementVolumeVariables &vDat)
+ {
+ static const Scalar eps = 1e-2;
+ const Scalar sat = vDat[scvIdx].density(phaseIdx);
+ if (sat > eps)
+ return 0.5;
+ if (sat <= 0)
+ return 0;
+
+ static const Dumux::Spline<Scalar> sp(0, eps, // x0, x1
+ 0, 0.5, // y0, y1
+ 0, 0); // m0, m1
+ return sp.eval(sat);
+ }
+
+ void calculateVelocities_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemDat)
+ {
+ const SpatialParameters &spatialParams = problem.spatialParameters();
+ // multiply the pressure potential with the intrinsic
+ // permeability
+ Tensor K;
+ for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
+ {
+ spatialParams.meanK(K,
+ spatialParams.intrinsicPermeability(element,
+ fvElemGeom_,
+ face().i),
+ spatialParams.intrinsicPermeability(element,
+ fvElemGeom_,
+ face().j));
+ K.mv(potentialGrad_[phaseIdx], Kmvp_[phaseIdx]);
+ KmvpNormal_[phaseIdx] = - (Kmvp_[phaseIdx] * face().normal);
+ }
+
+ // set the upstream and downstream vertices
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ upstreamIdx_[phaseIdx] = face().i;
+ downstreamIdx_[phaseIdx] = face().j;
+
+ if (KmvpNormal_[phaseIdx] < 0) {
+ std::swap(upstreamIdx_[phaseIdx],
+ downstreamIdx_[phaseIdx]);
+ }
+ }
+ }
+
+ void calculateDiffCoeffPM_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemDat)
+ {
+#if 0
+ const VolumeVariables &vDat_i = elemDat[face().i];
+ const VolumeVariables &vDat_j = elemDat[face().j];
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ // make sure to only calculate diffusion coefficents
+ // for phases which exist in both finite volumes
+ if (vDat_i.saturation(phaseIdx) <= 0 ||
+ vDat_j.saturation(phaseIdx) <= 0)
+ {
+ porousDiffCoeff_[phaseIdx] = 0.0;
+ continue;
+ }
+
+ // calculate tortuosity at the nodes i and j needed
+ // for porous media diffusion coefficient
+ Scalar tau_i =
+ 1.0/(vDat_i.porosity() * vDat_i.porosity()) *
+ pow(vDat_i.porosity() * vDat_i.saturation(phaseIdx), 7.0/3);
+ Scalar tau_j =
+ 1.0/(vDat_j.porosity() * vDat_j.porosity()) *
+ pow(vDat_j.porosity() * vDat_j.saturation(phaseIdx), 7.0/3);
+ // Diffusion coefficient in the porous medium
+
+ // -> harmonic mean
+ porousDiffCoeff_[phaseIdx] = harmonicMean(vDat_i.porosity() * vDat_i.saturation(phaseIdx) * tau_i * vDat_i.diffCoeff(phaseIdx),
+ vDat_j.porosity() * vDat_j.saturation(phaseIdx) * tau_j * vDat_j.diffCoeff(phaseIdx));
+ }
+#endif
+ }
+
+public:
+ /*!
+ * \brief Return the pressure potential multiplied with the
+ * intrinsic permeability which goes from vertex i to
+ * vertex j.
+ *
+ * Note that the length of the face's normal is the area of the
+ * phase, so this is not the actual velocity by the integral of
+ * the velocity over the face's area. Also note that the phase
+ * mobility is not yet included here since this would require a
+ * decision on the upwinding approach (which is done in the
+ * actual model).
+ */
+ Scalar KmvpNormal(int phaseIdx) const
+ { return KmvpNormal_[phaseIdx]; }
+
+ /*!
+ * \brief Return the pressure potential multiplied with the
+ * intrinsic permeability as vector (for velocity output)
+ */
+ Vector Kmvp(int phaseIdx) const
+ { return Kmvp_[phaseIdx]; }
+
+ /*!
+ * \brief Return the local index of the upstream control volume
+ * for a given phase.
+ */
+ int upstreamIdx(int phaseIdx) const
+ { return upstreamIdx_[phaseIdx]; }
+
+ /*!
+ * \brief Return the local index of the downstream control volume
+ * for a given phase.
+ */
+ int downstreamIdx(int phaseIdx) const
+ { return downstreamIdx_[phaseIdx]; }
+
+#if 0
+ /*!
+ * \brief The binary diffusion coefficient for each fluid phase.
+ */
+ Scalar porousDiffCoeff(int phaseIdx) const
+ { return porousDiffCoeff_[phaseIdx]; };
+#endif
+
+ /*!
+ * \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase at the integration
+ * point.
+ */
+ Scalar densityAtIP(int phaseIdx) const
+ { return densityAtIP_[phaseIdx]; }
+
+ /*!
+ * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ of a phase at the integration
+ * point.
+ */
+ Scalar molarDensityAtIP(int phaseIdx) const
+ { return molarDensityAtIP_[phaseIdx]; }
+
+ /*!
+ * \brief The concentration gradient of a component in a phase.
+ */
+ const Vector &concentrationGrad(int phaseIdx) const
+ { return concentrationGrad_[phaseIdx]; };
+
+ /*!
+ * \brief The molar concentration gradient of a component in a phase.
+ */
+ const Vector &molarConcGrad(int phaseIdx) const
+ { return molarConcGrad_[phaseIdx]; };
+
+ const SCVFace &face() const
+ { return fvElemGeom_.subContVolFace[scvfIdx_]; }
+
+protected:
+ const FVElementGeometry &fvElemGeom_;
+ int scvfIdx_;
+
+ // gradients
+ Vector potentialGrad_[numPhases];
+ Vector concentrationGrad_[numPhases];
+ Vector molarConcGrad_[numPhases];
+
+ // density of each face at the integration point
+ Scalar densityAtIP_[numPhases], molarDensityAtIP_[numPhases];
+
+ // intrinsic permeability times pressure potential gradient
+ Vector Kmvp_[numPhases];
+ // projected on the face normal
+ Scalar KmvpNormal_[numPhases];
+
+ // local index of the upwind vertex for each phase
+ int upstreamIdx_[numPhases];
+ // local index of the downwind vertex for each phase
+ int downstreamIdx_[numPhases];
+
+/*
+ // the diffusion coefficient for the porous medium
+ Scalar porousDiffCoeff_[numPhases];
+*/
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/boxmodels/new_2p2c/2p2cindices.hh b/dumux/boxmodels/new_2p2c/2p2cindices.hh
new file mode 100644
index 0000000000..8862e7f9e8
--- /dev/null
+++ b/dumux/boxmodels/new_2p2c/2p2cindices.hh
@@ -0,0 +1,128 @@
+// $Id: 2p2cproperties.hh 3784 2010-06-24 13:43:57Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2008 by Klaus Mosthaf, Andreas Lauser, Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+
+/*!
+ * \file
+ *
+ * \brief Defines the indices required for the 2p2c BOX model.
+ */
+#ifndef DUMUX_2P2C_INDICES_HH
+#define DUMUX_2P2C_INDICES_HH
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPTwoCModel
+ */
+// \{
+
+/*!
+ * \brief Enumerates the formulations which the 2p2c model accepts.
+ */
+struct TwoPTwoCFormulation
+{
+ enum {
+ plSg,
+ pgSl
+ };
+};
+
+/*!
+ * \brief The indices for the isothermal TwoPTwoC model.
+ *
+ * \tparam formulation The formulation, either pwSn or pnSw.
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag,
+ int formulation = TwoPTwoCFormulation::plSg,
+ int PVOffset = 0>
+class TwoPTwoCIndices
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ // Phase indices
+ static const int lPhaseIdx = FluidSystem::lPhaseIdx; //!< Index of the liquid phase
+ static const int gPhaseIdx = FluidSystem::gPhaseIdx; //!< Index of the gas phase
+
+ // Component indices
+ static const int lCompIdx = 0; //!< Index of the liquid's primary component
+ static const int gCompIdx = 1; //!< Index of the gas' primary component
+
+ // present phases (-> 'pseudo' primary variable)
+ static const int lPhaseOnly = 1; //!< Only the non-wetting phase is present
+ static const int gPhaseOnly = 0; //!< Only the wetting phase is present
+ static const int bothPhases = 2; //!< Both phases are present
+
+ // Primary variable indices
+ static const int pressureIdx = PVOffset + 0; //!< Index for wetting/non-wetting phase pressure (depending on formulation) in a solution vector
+ static const int switchIdx = PVOffset + 1; //!< Index of the either the saturation or the mass fraction of the non-wetting/wetting phase
+
+ static const int plIdx = pressureIdx; //!< Index for liquid phase pressure in a solution vector
+ static const int SgOrXIdx = switchIdx; //!< Index of the either the saturation of the gas phase or the mass fraction secondary component in the only phase
+
+ // equation indices
+ static const int contiLEqIdx = PVOffset + lCompIdx; //!< Index of the mass conservation equation for the liquid's primary component
+ static const int contiGEqIdx = PVOffset + gCompIdx; //!< Index of the mass conservation equation for the gas' primary component
+};
+
+/*!
+ * \brief The indices for the isothermal TwoPTwoC model in the pg-Sl
+ * formulation.
+ *
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag, int PVOffset>
+class TwoPTwoCIndices<TypeTag, TwoPTwoCFormulation::pgSl, PVOffset>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ // Phase indices
+ static const int lPhaseIdx = FluidSystem::lPhaseIdx; //!< Index of the liquid phase
+ static const int gPhaseIdx = FluidSystem::gPhaseIdx; //!< Index of the gas phase
+
+ // Component indices
+ static const int lCompIdx = 0; //!< Index of the liquid's primary component
+ static const int gCompIdx = 1; //!< Index of the gas' primary component
+
+ // present phases (-> 'pseudo' primary variable)
+ static const int lPhaseOnly = 1; //!< Only the non-wetting phase is present
+ static const int gPhaseOnly = 2; //!< Only the wetting phase is present
+ static const int bothPhases = 3; //!< Both phases are present
+
+ // Primary variable indices
+ static const int pressureIdx = PVOffset + 0; //!< Index for wetting/non-wetting phase pressure (depending on formulation) in a solution vector
+ static const int switchIdx = PVOffset + 1; //!< Index of the either the saturation or the mass fraction of the non-wetting/wetting phase
+
+ static const int pgIdx = pressureIdx; //!< Index for gas phase pressure in a solution vector
+ static const int SlOrXIdx = switchIdx; //!< Index of the either the saturation of the liquid phase or the mass fraction secondary component in the only phase
+
+ // Equation indices
+ static const int contiLEqIdx = PVOffset + lCompIdx; //!< Index of the mass conservation equation for the liquid's primary component
+ static const int contiGEqIdx = PVOffset + gCompIdx; //!< Index of the mass conservation equation for the gas' primary component
+};
+
+// \}
+
+}
+
+#endif
diff --git a/dumux/boxmodels/new_2p2c/2p2clocalresidual.hh b/dumux/boxmodels/new_2p2c/2p2clocalresidual.hh
new file mode 100644
index 0000000000..73d7c658e9
--- /dev/null
+++ b/dumux/boxmodels/new_2p2c/2p2clocalresidual.hh
@@ -0,0 +1,342 @@
+// $Id: 2p2clocalresidual.hh 3795 2010-06-25 16:08:04Z melanie $
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * Copyright (C) 2009-2010 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the two-phase two-component box model.
+ */
+
+#ifndef DUMUX_NEW_2P2C_LOCAL_RESIDUAL_BASE_HH
+#define DUMUX_NEW_2P2C_LOCAL_RESIDUAL_BASE_HH
+
+#include <dumux/boxmodels/common/boxmodel.hh>
+#include <dumux/common/math.hh>
+
+#include "2p2cproperties.hh"
+
+#include "2p2cvolumevariables.hh"
+
+#include "2p2cfluxvariables.hh"
+
+#include "2p2cnewtoncontroller.hh"
+
+#include <iostream>
+#include <vector>
+
+//#define VELOCITY_OUTPUT 1 // uncomment this line if an output of the velocity is needed
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPTwoCModel
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the two-phase two-component box model.
+ *
+ * This class is used to fill the gaps in BoxLocalResidual for the 2P-2C flow.
+ */
+template<class TypeTag>
+class TwoPTwoCLocalResidual: public BoxLocalResidual<TypeTag>
+{
+protected:
+ typedef TwoPTwoCLocalResidual<TypeTag> ThisType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) Implementation;
+ typedef BoxLocalResidual<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+
+ enum
+ {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
+ numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
+ numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),
+
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx,
+
+ contiLEqIdx = Indices::contiLEqIdx,
+ contiGEqIdx = Indices::contiGEqIdx,
+
+ lPhaseIdx = Indices::lPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+
+ lCompIdx = Indices::lCompIdx,
+ gCompIdx = Indices::gCompIdx,
+
+ lPhaseOnly = Indices::lPhaseOnly,
+ gPhaseOnly = Indices::gPhaseOnly,
+ bothPhases = Indices::bothPhases,
+
+ plSg = TwoPTwoCFormulation::plSg,
+ pgSl = TwoPTwoCFormulation::pgSl,
+ formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation))
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::ctype CoordScalar;
+
+ typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
+
+ static const Scalar mobilityUpwindAlpha =
+ GET_PROP_VALUE(TypeTag, PTAG(MobilityUpwindAlpha));
+
+public:
+ /*!
+ * \brief Evaluate the storage term of the current solution in a
+ * single phase.
+ *
+ * \param element The element
+ * \param phaseIdx The index of the fluid phase
+ */
+ void evalPhaseStorage(const Element &element, int phaseIdx)
+ {
+ FVElementGeometry fvGeom;
+ fvGeom.update(this->gridView_(), element);
+ ElementBoundaryTypes bcTypes;
+ bcTypes.update(this->problem_(), element, fvGeom);
+ ElementVolumeVariables volVars;
+ volVars.update(this->problem_(), element, fvGeom, false);
+
+ this->residual_.resize(fvGeom.numVertices);
+ this->residual_ = 0;
+
+ this->elemPtr_ = &element;
+ this->fvElemGeomPtr_ = &fvGeom;
+ this->bcTypesPtr_ = &bcTypes;
+ this->prevVolVarsPtr_ = 0;
+ this->curVolVarsPtr_ = &volVars;
+ evalPhaseStorage_(phaseIdx);
+ }
+
+ /*!
+ * \brief Evaluate the amount all conservation quantities
+ * (e.g. phase mass) within a sub-control volume.
+ *
+ * The result should be averaged over the volume (e.g. phase mass
+ * inside a sub control volume divided by the volume)
+ *
+ * \param result The mass of the component within the sub-control volume
+ * \param scvIdx The SCV (sub-control-volume) index
+ * \param usePrevSol Evaluate function with solution of current or previous time step
+ */
+ void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const
+ {
+ // if flag usePrevSol is set, the solution from the previous
+ // time step is used, otherwise the current solution is
+ // used. The secondary variables are used accordingly. This
+ // is required to compute the derivative of the storage term
+ // using the implicit euler method.
+ const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_()
+ : this->curVolVars_();
+ const VolumeVariables &volVars = elemVolVars[scvIdx];
+
+ // compute storage term of all components within all phases
+ result = 0;
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ int eqIdx = (compIdx == lCompIdx) ? contiLEqIdx : contiGEqIdx;
+ result[eqIdx] += volVars.density(phaseIdx)
+ * volVars.saturation(phaseIdx)
+ * volVars.fluidState().massFrac(phaseIdx, compIdx);
+ }
+ }
+ result *= volVars.porosity();
+ }
+
+ /*!
+ * \brief Evaluates the total flux of all conservation quantities
+ * over a face of a sub-control volume.
+ *
+ * \param flux The flux over the SCV (sub-control-volume) face for each component
+ * \param faceIdx The index of the SCV face
+ */
+ void computeFlux(PrimaryVariables &flux, int faceIdx) const
+ {
+ FluxVariables vars(this->problem_(),
+ this->elem_(),
+ this->fvElemGeom_(),
+ faceIdx,
+ this->curVolVars_());
+
+ flux = 0;
+ asImp_()->computeAdvectiveFlux(flux, vars);
+ Valgrind::CheckDefined(flux);
+ asImp_()->computeDiffusiveFlux(flux, vars);
+ Valgrind::CheckDefined(flux);
+ }
+
+ /*!
+ * \brief Evaluates the advective mass flux of all components over
+ * a face of a subcontrol volume.
+ *
+ * \param flux The advective flux over the sub-control-volume face for each component
+ * \param vars The flux variables at the current SCV
+ */
+ void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &vars) const
+ {
+ ////////
+ // advective fluxes of all components in all phases
+ ////////
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ // data attached to upstream and the downstream vertices
+ // of the current phase
+ const VolumeVariables &up =
+ this->curVolVars_(vars.upstreamIdx(phaseIdx));
+ const VolumeVariables &dn =
+ this->curVolVars_(vars.downstreamIdx(phaseIdx));
+
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ int eqIdx = (compIdx == lCompIdx) ? contiLEqIdx : contiGEqIdx;
+ // add advective flux of current component in current
+ // phase
+ if (mobilityUpwindAlpha > 0.0)
+ // upstream vertex
+ flux[eqIdx] += vars.KmvpNormal(phaseIdx)
+ * mobilityUpwindAlpha * (up.density(phaseIdx)
+ * up.mobility(phaseIdx) * up.fluidState().massFrac(
+ phaseIdx, compIdx));
+ if (mobilityUpwindAlpha < 1.0)
+ // downstream vertex
+ flux[eqIdx] += vars.KmvpNormal(phaseIdx) * (1
+ - mobilityUpwindAlpha) * (dn.density(phaseIdx)
+ * dn.mobility(phaseIdx) * dn.fluidState().massFrac(
+ phaseIdx, compIdx));
+ Valgrind::CheckDefined(vars.KmvpNormal(phaseIdx));
+ Valgrind::CheckDefined(up.density(phaseIdx));
+ Valgrind::CheckDefined(up.mobility(phaseIdx));
+ Valgrind::CheckDefined(up.fluidState().massFrac(phaseIdx, compIdx));
+ Valgrind::CheckDefined(dn.density(phaseIdx));
+ Valgrind::CheckDefined(dn.mobility(phaseIdx));
+ Valgrind::CheckDefined(dn.fluidState().massFrac(phaseIdx, compIdx));
+ }
+ }
+ }
+
+ /*!
+ * \brief Adds the diffusive mass flux of all components over
+ * a face of a subcontrol volume.
+ *
+ * \param flux The diffusive flux over the sub-control-volume face for each component
+ * \param vars The flux variables at the current SCV
+ */
+ void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &vars) const
+ {
+#if 0
+ // add diffusive flux of gas component in liquid phase
+ Scalar tmp =
+ - vars.porousDiffCoeff(lPhaseIdx) *
+ vars.molarDensityAtIP(lPhaseIdx) *
+ (vars.molarConcGrad(lPhaseIdx) * vars.face().normal);
+ flux[contiGEqIdx] += tmp * FluidSystem::molarMass(gCompIdx);
+ flux[contiLEqIdx] -= tmp * FluidSystem::molarMass(lCompIdx);
+
+ // add diffusive flux of liquid component in gas phase
+ tmp =
+ - vars.porousDiffCoeff(gPhaseIdx) *
+ vars.molarDensityAtIP(gPhaseIdx) *
+ (vars.molarConcGrad(gPhaseIdx) * vars.face().normal);
+ flux[contiLEqIdx] += tmp * FluidSystem::molarMass(lCompIdx);
+ flux[contiGEqIdx] -= tmp * FluidSystem::molarMass(gCompIdx);
+#endif
+ }
+
+ /*!
+ * \brief Calculate the source term of the equation
+ *
+ * \param q The source/sink in the SCV for each component
+ * \param localVertexIdx The index of the SCV
+ */
+ void computeSource(PrimaryVariables &q, int localVertexIdx)
+ {
+ this->problem_().source(q,
+ this->elem_(),
+ this->fvElemGeom_(),
+ localVertexIdx);
+ }
+
+protected:
+ void evalPhaseStorage_(int phaseIdx)
+ {
+ // evaluate the storage terms of a single phase
+ for (int i=0; i < this->fvElemGeom_().numVertices; i++) {
+ PrimaryVariables &result = this->residual_[i];
+ const ElementVolumeVariables &elemVolVars = this->curVolVars_();
+ const VolumeVariables &volVars = elemVolVars[i];
+
+ // compute storage term of all components within all phases
+ result = 0;
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ int eqIdx = (compIdx == lCompIdx) ? contiLEqIdx : contiGEqIdx;
+ result[eqIdx] += volVars.density(phaseIdx)
+ * volVars.saturation(phaseIdx)
+ * volVars.fluidState().massFrac(phaseIdx, compIdx);
+ }
+
+ result *= volVars.porosity();
+ result *= this->fvElemGeom_().subContVol[i].volume;
+ }
+ }
+
+
+ Implementation *asImp_()
+ {
+ return static_cast<Implementation *> (this);
+ }
+ const Implementation *asImp_() const
+ {
+ return static_cast<const Implementation *> (this);
+ }
+
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/boxmodels/new_2p2c/2p2cmodel.hh b/dumux/boxmodels/new_2p2c/2p2cmodel.hh
new file mode 100644
index 0000000000..c59755e211
--- /dev/null
+++ b/dumux/boxmodels/new_2p2c/2p2cmodel.hh
@@ -0,0 +1,760 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008 by Klaus Mosthaf, Andreas Lauser, Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+
+/*!
+* \file
+*
+* \brief Adaption of the BOX scheme to the two-phase two-component flow model.
+*/
+
+#ifndef DUMUX_2P2C_MODEL_HH
+#define DUMUX_2P2C_MODEL_HH
+
+#include "2p2cproperties.hh"
+#include "2p2clocalresidual.hh"
+#include "2p2cproblem.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup BoxModels
+ * \defgroup TwoPTwoCModel Two-phase two-component box model
+ */
+
+/*!
+ * \ingroup TwoPTwoCModel
+ * \brief Adaption of the BOX scheme to the two-phase two-component flow model.
+ *
+ * This model implements two-phase two-component flow of two compressible and
+ * partially miscible fluids \f$\alpha \in \{ w, n \}\f$ composed of the two components
+ * \f$\kappa \in \{ w, a \}\f$. The standard multiphase Darcy
+ * approach is used as the equation for the conservation of momentum:
+ * \f[
+ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K}
+ \left(\text{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right)
+ * \f]
+ *
+ * By inserting this into the equations for the conservation of the
+ * components, one gets one transport equation for each component
+ * \f{eqnarray}
+ && \phi \frac{\partial (\sum_\alpha \varrho_\alpha X_\alpha^\kappa S_\alpha )}
+ {\partial t}
+ - \sum_\alpha \text{div} \left\{ \varrho_\alpha X_\alpha^\kappa
+ \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K}
+ (\text{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g}) \right\}
+ \nonumber \\ \nonumber \\
+ &-& \sum_\alpha \text{div} \left\{{\bf D_{\alpha, pm}^\kappa} \varrho_{\alpha} \text{grad}\, X^\kappa_{\alpha} \right\}
+ - \sum_\alpha q_\alpha^\kappa = 0 \qquad \kappa \in \{w, a\} \, ,
+ \alpha \in \{w, g\}
+ \f}
+ *
+ * This is discretized using a fully-coupled vertex
+ * centered finite volume (box) scheme as spatial and
+ * the implicit Euler method as temporal discretization.
+ *
+ * By using constitutive relations for the capillary pressure \f$p_c =
+ * p_n - p_w\f$ and relative permeability \f$k_{r\alpha}\f$ and taking
+ * advantage of the fact that \f$S_w + S_n = 1\f$ and \f$X^\kappa_w + X^\kappa_n = 1\f$, the number of
+ * unknowns can be reduced to two.
+ * The used primary variables are, like in the two-phase model, either \f$p_w\f$ and \f$S_n\f$
+ * or \f$p_n\f$ and \f$S_w\f$. The formulation which ought to be used can be
+ * specified by setting the <tt>Formulation</tt> property to either
+ * TwoPTwoCIndices::pWsN or TwoPTwoCIndices::pNsW. By
+ * default, the model uses \f$p_w\f$ and \f$S_n\f$.
+ * Moreover, the second primary variable depends on the phase state, since a
+ * primary variable switch is included. The phase state is stored for all nodes
+ * of the system. Following cases can be distinguished:
+ * <ul>
+ * <li> Both phases are present: The saturation is used (either \f$S_n\f$ or \f$S_w\f$, dependent on the chosen <tt>Formulation</tt>),
+ * as long as \f$ 0 < S_\alpha < 1\f$</li>.
+ * <li> Only wetting phase is present: The mass fraction of, e.g., air in the wetting phase \f$X^a_w\f$ is used,
+ * as long as the maximum mass fraction is not exceeded (\f$X^a_w<X^a_{w,max}\f$)</li>
+ * <li> Only non-wetting phase is present: The mass fraction of, e.g., water in the non-wetting phase, \f$X^w_n\f$, is used,
+ * as long as the maximum mass fraction is not exceeded (\f$X^w_n<X^w_{n,max}\f$)</li>
+ * </ul>
+ */
+
+template<class TypeTag>
+class TwoPTwoCModel: public BoxModel<TypeTag>
+{
+ typedef TwoPTwoCModel<TypeTag> ThisType;
+ typedef BoxModel<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementMapper)) ElementMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
+ numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
+ numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),
+
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx,
+
+ contiLEqIdx = Indices::contiLEqIdx,
+ contiGEqIdx = Indices::contiGEqIdx,
+
+ lPhaseIdx = Indices::lPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+
+ lCompIdx = Indices::lCompIdx,
+ gCompIdx = Indices::gCompIdx,
+
+ lPhaseOnly = Indices::lPhaseOnly,
+ gPhaseOnly = Indices::gPhaseOnly,
+ bothPhases = Indices::bothPhases,
+
+ plSg = TwoPTwoCFormulation::plSg,
+ pgSl = TwoPTwoCFormulation::pgSl,
+ formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation))
+ };
+
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ static const Scalar mobilityUpwindAlpha =
+ GET_PROP_VALUE(TypeTag, PTAG(MobilityUpwindAlpha));
+
+public:
+ /*!
+ * \brief Initialize the static data with the initial solution.
+ *
+ * \param problem The problem to be solved
+ */
+ void init(Problem &problem)
+ {
+ ParentType::init(problem);
+
+ staticVertexDat_.resize(this->gridView_().size(dim));
+
+ setSwitched_(false);
+
+ VertexIterator it = this->gridView_().template begin<dim> ();
+ const VertexIterator &endit = this->gridView_().template end<dim> ();
+ for (; it != endit; ++it)
+ {
+ int globalIdx = this->dofMapper().map(*it);
+ const GlobalPosition &globalPos = it->geometry().corner(0);
+
+ // initialize phase presence
+ staticVertexDat_[globalIdx].phasePresence
+ = this->problem_().initialPhasePresence(*it, globalIdx,
+ globalPos);
+ staticVertexDat_[globalIdx].wasSwitched = false;
+
+ staticVertexDat_[globalIdx].oldPhasePresence
+ = staticVertexDat_[globalIdx].phasePresence;
+ }
+ }
+
+ /*!
+ * \brief Compute the total storage inside one phase of all
+ * conservation quantities.
+ *
+ * \param dest Contains the storage of each component for one phase
+ * \param phaseIdx The phase index
+ */
+ void globalPhaseStorage(PrimaryVariables &dest, int phaseIdx)
+ {
+ dest = 0;
+
+ ElementIterator elemIt = this->gridView_().template begin<0>();
+ const ElementIterator elemEndIt = this->gridView_().template end<0>();
+ for (; elemIt != elemEndIt; ++elemIt) {
+ this->localResidual().evalPhaseStorage(*elemIt, phaseIdx);
+
+ for (int i = 0; i < elemIt->template count<dim>(); ++i)
+ dest += this->localResidual().residual(i);
+ };
+
+ if (this->gridView_().comm().size() > 1)
+ dest = this->gridView_().comm().sum(dest);
+ }
+
+ /*!
+ * \brief Called by the update() method if applying the newton
+ * method was unsuccessful.
+ */
+ void updateFailed()
+ {
+ ParentType::updateFailed();
+
+ setSwitched_(false);
+ resetPhasePresence_();
+ };
+
+ /*!
+ * \brief Returns the relative weight of a primary variable for
+ * calculating relative errors.
+ *
+ * \param globalVertexIdx The global vertex index
+ * \param pvIdx The primary variable index
+ */
+ Scalar primaryVarWeight(int globalVertexIdx, int pvIdx) const
+ {
+ if (Indices::pressureIdx == pvIdx)
+ return std::min(1.0/this->prevSol()[globalVertexIdx][pvIdx], 1.0);
+ return 1;
+ }
+
+ /*!
+ * \brief Called by the problem if a time integration was
+ * successful, post processing of the solution is done and the
+ * result has been written to disk.
+ *
+ * This should prepare the model for the next time integration.
+ */
+ void advanceTimeLevel()
+ {
+ ParentType::advanceTimeLevel();
+
+ // update the phase state
+ updateOldPhasePresence_();
+ setSwitched_(false);
+ }
+
+ /*!
+ * \brief Return true if the primary variables were switched for
+ * at least one vertex after the last timestep.
+ */
+ bool switched() const
+ {
+ return switchFlag_;
+ }
+
+ /*!
+ * \brief Returns the phase presence of the current or the old solution of a vertex.
+ *
+ * \param globalVertexIdx The global vertex index
+ * \param oldSol Evaluate function with solution of current or previous time step
+ */
+ int phasePresence(int globalVertexIdx, bool oldSol) const
+ {
+ return oldSol ? staticVertexDat_[globalVertexIdx].oldPhasePresence
+ : staticVertexDat_[globalVertexIdx].phasePresence;
+ }
+
+ /*!
+ * \brief Append all quantities of interest which can be derived
+ * from the solution of the current time step to the VTK
+ * writer.
+ *
+ * \param sol The solution vector
+ * \param writer The writer for multi-file VTK datasets
+ */
+ template<class MultiWriter>
+ void addOutputVtkFields(const SolutionVector &sol,
+ MultiWriter &writer)
+ {
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
+
+ // create the required scalar fields
+ unsigned numVertices = this->problem_().gridView().size(dim);
+
+ ScalarField *Sg = writer.template createField<Scalar, 1> (numVertices);
+ ScalarField *Sl = writer.template createField<Scalar, 1> (numVertices);
+ ScalarField *pg = writer.template createField<Scalar, 1> (numVertices);
+ ScalarField *pl = writer.template createField<Scalar, 1> (numVertices);
+ ScalarField *pc = writer.template createField<Scalar, 1> (numVertices);
+ ScalarField *rhoL =
+ writer.template createField<Scalar, 1> (numVertices);
+ ScalarField *rhoG =
+ writer.template createField<Scalar, 1> (numVertices);
+ ScalarField *mobL =
+ writer.template createField<Scalar, 1> (numVertices);
+ ScalarField *mobG =
+ writer.template createField<Scalar, 1> (numVertices);
+ ScalarField *phasePresence = writer.template createField<Scalar, 1> (
+ numVertices);
+ ScalarField *massFrac[numPhases][numComponents];
+ for (int i = 0; i < numPhases; ++i)
+ for (int j = 0; j < numComponents; ++j)
+ massFrac[i][j] = writer.template createField<Scalar, 1>(numVertices);
+ ScalarField *temperature = writer.template createField<Scalar, 1>(numVertices);
+ ScalarField *poro = writer.template createField<Scalar, 1>(numVertices);
+
+#ifdef VELOCITY_OUTPUT // check if velocity output is demanded
+ ScalarField *velocityX = writer.template createField<Scalar, 1>(numVertices);
+ ScalarField *velocityY = writer.template createField<Scalar, 1>(numVertices);
+ ScalarField *velocityZ = writer.template createField<Scalar, 1>(numVertices);
+// Scalar maxV=0.; // variable to store the maximum face velocity
+
+ // initialize velocity fields
+ Scalar boxSurface[numVertices];
+ for (int i = 0; i < numVertices; ++i)
+ {
+ (*velocityX)[i] = 0;
+ if (dim > 1)
+ (*velocityY)[i] = 0;
+ if (dim > 2)
+ (*velocityZ)[i] = 0;
+ boxSurface[i] = 0.0; // initialize the boundary surface of the fv-boxes
+ }
+#endif
+
+ unsigned numElements = this->gridView_().size(0);
+ ScalarField *rank =
+ writer.template createField<Scalar, 1> (numElements);
+
+ FVElementGeometry fvElemGeom;
+ VolumeVariables volVars;
+
+ ElementIterator elemIt = this->gridView_().template begin<0>();
+ ElementIterator elemEndIt = this->gridView_().template end<0>();
+ for (; elemIt != elemEndIt; ++elemIt)
+ {
+ int idx = this->problem_().elementMapper().map(*elemIt);
+ (*rank)[idx] = this->gridView_().comm().rank();
+ fvElemGeom.update(this->gridView_(), *elemIt);
+
+ int numVerts = elemIt->template count<dim> ();
+ for (int i = 0; i < numVerts; ++i)
+ {
+ int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ volVars.update(sol[globalIdx],
+ this->problem_(),
+ *elemIt,
+ fvElemGeom,
+ i,
+ false);
+ (*Sg)[globalIdx] = volVars.saturation(gPhaseIdx);
+ (*Sl)[globalIdx] = volVars.saturation(lPhaseIdx);
+ (*pg)[globalIdx] = volVars.pressure(gPhaseIdx);
+ (*pl)[globalIdx] = volVars.pressure(lPhaseIdx);
+ (*pc)[globalIdx] = volVars.capillaryPressure();
+ (*rhoL)[globalIdx] = volVars.fluidState().density(lPhaseIdx);
+ (*rhoG)[globalIdx] = volVars.fluidState().density(gPhaseIdx);
+ (*mobL)[globalIdx] = volVars.mobility(lPhaseIdx);
+ (*mobG)[globalIdx] = volVars.mobility(gPhaseIdx);
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ (*massFrac[phaseIdx][compIdx])[globalIdx]
+ = volVars.fluidState().massFrac(phaseIdx,
+ compIdx);
+
+ Valgrind::CheckDefined(
+ (*massFrac[phaseIdx][compIdx])[globalIdx][0]);
+ }
+ (*poro)[globalIdx] = volVars.porosity();
+ (*temperature)[globalIdx] = volVars.temperature();
+ (*phasePresence)[globalIdx]
+ = staticVertexDat_[globalIdx].phasePresence;
+ };
+
+#ifdef VELOCITY_OUTPUT // check if velocity output is demanded
+ // In the box method, the velocity is evaluated on the FE-Grid. However, to get an
+ // average apparent velocity at the vertex, all contributing velocities have to be interpolated.
+ GlobalPosition velocity(0.);
+ ElementVolumeVariables elemVolVars;
+
+ elemVolVars.update(this->problem_(),
+ *elemIt,
+ fvElemGeom,
+ false /* isOldSol? */);
+
+ // loop over the phases
+ for (int faceIdx = 0; faceIdx< fvElemGeom.numEdges; faceIdx++)
+ {
+ //prepare the flux calculations (set up and prepare geometry, FE gradients)
+ FluxVariables fluxDat(this->problem_(),
+ *elemIt,
+ fvElemGeom,
+ faceIdx,
+ elemVolVars);
+
+ // choose phase of interest. Alternatively, a loop over all phases would be possible.
+ int phaseIdx = gPhaseIdx;
+
+ // get darcy velocity
+ velocity = fluxDat.Kmvp(phaseIdx); // mind the sign: vDarcy = kf grad p
+
+ // up+downstream mobility
+ const VolumeVariables &up = elemVolVars[fluxDat.upstreamIdx(phaseIdx)];
+ const VolumeVariables &down = elemVolVars[fluxDat.downstreamIdx(phaseIdx)];
+ Scalar scvfArea = fluxDat.face().normal.two_norm(); //get surface area to weight velocity at the IP with the surface area
+ velocity *= (mobilityUpwindAlpha*up.mobility(phaseIdx) + (1-mobilityUpwindAlpha)*down.mobility(phaseIdx))* scvfArea;
+
+ int vertIIdx = this->problem_().vertexMapper().map(*elemIt,
+ fluxDat.face().i,
+ dim);
+ int vertJIdx = this->problem_().vertexMapper().map(*elemIt,
+ fluxDat.face().j,
+ dim);
+ // add surface area for weighting purposes
+ boxSurface[vertIIdx] += scvfArea;
+ boxSurface[vertJIdx] += scvfArea;
+
+ // Add velocity to upstream and downstream vertex.
+ // Beware: velocity has to be substracted because of the (wrong) sign of vDarcy
+ (*velocityX)[vertIIdx] -= velocity[0];
+ (*velocityX)[vertJIdx] -= velocity[0];
+ if (dim >= 2)
+ {
+ (*velocityY)[vertIIdx] -= velocity[1];
+ (*velocityY)[vertJIdx] -= velocity[1];
+ }
+ if (dim == 3)
+ {
+ (*velocityZ)[vertIIdx] -= velocity[2];
+ (*velocityZ)[vertJIdx] -= velocity[2];
+ }
+ }
+#endif
+ }
+
+#ifdef VELOCITY_OUTPUT // check if velocity output is demanded
+ // normalize the velocities at the vertices
+ for (int i = 0; i < numVertices; ++i)
+ {
+ (*velocityX)[i] /= boxSurface[i];
+ if (dim >= 2)
+ (*velocityY)[i] /= boxSurface[i];
+ if (dim == 3)
+ (*velocityZ)[i] /= boxSurface[i];
+ }
+#endif
+
+ writer.addVertexData(Sg, "Sg");
+ writer.addVertexData(Sl, "Sl");
+ writer.addVertexData(pg, "pg");
+ writer.addVertexData(pl, "pl");
+ writer.addVertexData(pc, "pc");
+ writer.addVertexData(rhoL, "rhoL");
+ writer.addVertexData(rhoG, "rhoG");
+ writer.addVertexData(mobL, "mobL");
+ writer.addVertexData(mobG, "mobG");
+ for (int i = 0; i < numPhases; ++i)
+ {
+ for (int j = 0; j < numComponents; ++j)
+ {
+ std::string name = (boost::format("X_%s%s")
+ % ((i == lPhaseIdx) ? "l" : "g")
+ % FluidSystem::componentName(j)).str();
+ writer.addVertexData(massFrac[i][j], name.c_str());
+ }
+ }
+ writer.addVertexData(poro, "porosity");
+ writer.addVertexData(temperature, "temperature");
+ writer.addVertexData(phasePresence, "phase presence");
+
+#ifdef VELOCITY_OUTPUT // check if velocity output is demanded
+ writer.addVertexData(velocityX, "Vx");
+ if (dim >= 2)
+ writer.addVertexData(velocityY, "Vy");
+ if (dim == 3)
+ writer.addVertexData(velocityZ, "Vz");
+#endif
+ writer.addCellData(rank, "process rank");
+ }
+
+ /*!
+ * \brief Write the current solution to a restart file.
+ *
+ * \param outStream The output stream of one vertex for the restart file
+ * \param vert The vertex
+ */
+ void serializeEntity(std::ostream &outStream, const Vertex &vert)
+ {
+ // write primary variables
+ ParentType::serializeEntity(outStream, vert);
+
+ int vertIdx = this->dofMapper().map(vert);
+ if (!outStream.good())
+ DUNE_THROW(Dune::IOError, "Could not serialize vertex " << vertIdx);
+
+ outStream << staticVertexDat_[vertIdx].phasePresence << " ";
+ }
+
+ /*!
+ * \brief Reads the current solution for a vertex from a restart
+ * file.
+ *
+ * \param inStream The input stream of one vertex from the restart file
+ * \param vert The vertex
+ */
+ void deserializeEntity(std::istream &inStream, const Vertex &vert)
+ {
+ // read primary variables
+ ParentType::deserializeEntity(inStream, vert);
+
+ // read phase presence
+ int vertIdx = this->dofMapper().map(vert);
+ if (!inStream.good())
+ DUNE_THROW(Dune::IOError,
+ "Could not deserialize vertex " << vertIdx);
+
+ inStream >> staticVertexDat_[vertIdx].phasePresence;
+ staticVertexDat_[vertIdx].oldPhasePresence
+ = staticVertexDat_[vertIdx].phasePresence;
+
+ }
+
+ /*!
+ * \brief Update the static data of all vertices in the grid.
+ *
+ * \param curGlobalSol The current global solution
+ * \param oldGlobalSol The previous global solution
+ */
+ void updateStaticData(SolutionVector &curGlobalSol,
+ const SolutionVector &oldGlobalSol)
+ {
+ bool wasSwitched = false;
+
+ for (unsigned i = 0; i < staticVertexDat_.size(); ++i)
+ staticVertexDat_[i].visited = false;
+
+ FVElementGeometry fvElemGeom;
+ static VolumeVariables volVars;
+ ElementIterator it = this->gridView_().template begin<0> ();
+ const ElementIterator &endit = this->gridView_().template end<0> ();
+ for (; it != endit; ++it)
+ {
+ fvElemGeom.update(this->gridView_(), *it);
+ for (int i = 0; i < fvElemGeom.numVertices; ++i)
+ {
+ int globalIdx = this->vertexMapper().map(*it, i, dim);
+
+ if (staticVertexDat_[globalIdx].visited)
+ continue;
+
+ staticVertexDat_[globalIdx].visited = true;
+ volVars.update(curGlobalSol[globalIdx],
+ this->problem_(),
+ *it,
+ fvElemGeom,
+ i,
+ false);
+ const GlobalPosition &global = it->geometry().corner(i);
+ if (primaryVarSwitch_(curGlobalSol,
+ volVars,
+ globalIdx,
+ global))
+ wasSwitched = true;
+ }
+ }
+
+ // make sure that if there was a variable switch in an
+ // other partition we will also set the switch flag
+ // for our partition.
+ if (this->gridView_().comm().size() > 1)
+ wasSwitched = this->gridView_().comm().max(wasSwitched);
+
+ setSwitched_(wasSwitched);
+ }
+
+protected:
+ /*!
+ * \brief Data which is attached to each vertex and is not only
+ * stored locally.
+ */
+ struct StaticVars
+ {
+ int phasePresence;
+ bool wasSwitched;
+
+ int oldPhasePresence;
+ bool visited;
+ };
+
+ /*!
+ * \brief Reset the current phase presence of all vertices to the old one.
+ *
+ * This is done after an update failed.
+ */
+ void resetPhasePresence_()
+ {
+ int numVertices = this->gridView_().size(dim);
+ for (int i = 0; i < numVertices; ++i)
+ {
+ staticVertexDat_[i].phasePresence
+ = staticVertexDat_[i].oldPhasePresence;
+ staticVertexDat_[i].wasSwitched = false;
+ }
+ }
+
+ /*!
+ * \brief Set the old phase of all verts state to the current one.
+ */
+ void updateOldPhasePresence_()
+ {
+ int numVertices = this->gridView_().size(dim);
+ for (int i = 0; i < numVertices; ++i)
+ {
+ staticVertexDat_[i].oldPhasePresence
+ = staticVertexDat_[i].phasePresence;
+ staticVertexDat_[i].wasSwitched = false;
+ }
+ }
+
+ /*!
+ * \brief Set whether there was a primary variable switch after in
+ * the last timestep.
+ */
+ void setSwitched_(bool yesno)
+ {
+ switchFlag_ = yesno;
+ }
+
+ // perform variable switch at a vertex; Returns true if a
+ // variable switch was performed.
+ bool primaryVarSwitch_(SolutionVector &globalSol,
+ const VolumeVariables &volVars,
+ int globalIdx,
+ const GlobalPosition &globalPos)
+ {
+ // evaluate primary variable switch
+ bool wouldSwitch = false;
+ int phasePresence = staticVertexDat_[globalIdx].phasePresence;
+ int newPhasePresence = phasePresence;
+
+ // check if a primary var switch is necessary
+ if (phasePresence == gPhaseOnly)
+ {
+ // calculate mole fraction in the hypothetic liquid phase
+ Scalar xll = volVars.fluidState().moleFrac(lPhaseIdx, lCompIdx);
+ Scalar xlg = volVars.fluidState().moleFrac(lPhaseIdx, gCompIdx);
+
+ Scalar xlMax = 1.0;
+ if (xll + xlg > xlMax)
+ wouldSwitch = true;
+ if (staticVertexDat_[globalIdx].wasSwitched)
+ xlMax *= 1.02;
+
+ // if the sum of the mole fractions would be larger than
+ // 100%, liquid phase appears
+ if (xll + xlg > xlMax)
+ {
+ // liquid phase appears
+ /*
+ std::cout << "Liquid phase appears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", xll + xlg: "
+ << xll + xlg << std::endl;
+ */
+ newPhasePresence = bothPhases;
+ if (formulation == pgSl)
+ globalSol[globalIdx][switchIdx] = 0.001;
+ else if (formulation == plSg)
+ globalSol[globalIdx][switchIdx] = 0.999;
+ };
+ }
+ else if (phasePresence == lPhaseOnly)
+ {
+ // calculate fractions of the partial pressures in the
+ // hypothetic gas phase
+ Scalar xgl = volVars.fluidState().moleFrac(gPhaseIdx, lCompIdx);
+ Scalar xgg = volVars.fluidState().moleFrac(gPhaseIdx, gCompIdx);
+
+ Scalar xgMax = 1.0;
+ if (xgl + xgg > xgMax)
+ wouldSwitch = true;
+ if (staticVertexDat_[globalIdx].wasSwitched)
+ xgMax *= 1.02;
+
+ // if the sum of the mole fractions would be larger than
+ // 100%, gas phase appears
+ if (xgl + xgg > xgMax)
+ {
+ // gas phase appears
+ /*
+ std::cout << "gas phase appears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", x_gl + x_gg: "
+ << xgl + xgg << std::endl;
+ */
+ newPhasePresence = bothPhases;
+ if (formulation == pgSl)
+ globalSol[globalIdx][switchIdx] = 0.999;
+ else if (formulation == plSg)
+ globalSol[globalIdx][switchIdx] = 0.001;
+ }
+ }
+ else if (phasePresence == bothPhases)
+ {
+ Scalar Smin = 0.0;
+ if (staticVertexDat_[globalIdx].wasSwitched)
+ Smin = -0.01;
+
+ if (volVars.saturation(gPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // gas phase disappears
+ /*
+ std::cout << "Gas phase disappears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", Sg: "
+ << volVars.saturation(gPhaseIdx) << std::endl;
+ */
+ newPhasePresence = lPhaseOnly;
+
+ globalSol[globalIdx][switchIdx]
+ = volVars.fluidState().massFrac(lPhaseIdx, gCompIdx);
+ }
+ else if (volVars.saturation(lPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // liquid phase disappears
+ /*
+ std::cout << "Liquid phase disappears at vertex " << globalIdx
+ << ", coordinates: " << globalPos << ", Sl: "
+ << volVars.saturation(lPhaseIdx) << std::endl;
+ */
+ newPhasePresence = gPhaseOnly;
+
+ globalSol[globalIdx][switchIdx]
+ = volVars.fluidState().massFrac(gPhaseIdx, lCompIdx);
+ }
+ }
+
+ staticVertexDat_[globalIdx].phasePresence = newPhasePresence;
+ staticVertexDat_[globalIdx].wasSwitched = wouldSwitch;
+ return phasePresence != newPhasePresence;
+ }
+
+ // parameters given in constructor
+ std::vector<StaticVars> staticVertexDat_;
+ bool switchFlag_;
+};
+
+}
+
+#include "2p2cpropertydefaults.hh"
+
+#endif
diff --git a/dumux/boxmodels/new_2p2c/2p2cnewtoncontroller.hh b/dumux/boxmodels/new_2p2c/2p2cnewtoncontroller.hh
new file mode 100644
index 0000000000..1b4ad9376f
--- /dev/null
+++ b/dumux/boxmodels/new_2p2c/2p2cnewtoncontroller.hh
@@ -0,0 +1,181 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008 by Bernd Flemisch, Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief A 2p2c specific controller for the newton solver.
+ *
+ * This controller 'knows' what a 'physically meaningful' solution is
+ * which allows the newton method to abort quicker if the solution is
+ * way out of bounds.
+ */
+#ifndef DUMUX_2P2C_NEWTON_CONTROLLER_HH
+#define DUMUX_2P2C_NEWTON_CONTROLLER_HH
+
+#include "2p2cproperties.hh"
+
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+namespace Dumux {
+/*!
+ * \ingroup TwoPTwoCModel
+ * \brief A 2p2c specific controller for the newton solver.
+ *
+ * This controller 'knows' what a 'physically meaningful' solution is
+ * which allows the newton method to abort quicker if the solution is
+ * way out of bounds.
+ */
+template <class TypeTag>
+class TwoPTwoCNewtonController : public NewtonController<TypeTag>
+{
+ typedef NewtonController<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
+
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+
+ enum {
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx
+ };
+
+ enum { enablePartialReassemble = GET_PROP_VALUE(TypeTag, PTAG(EnablePartialReassemble)) };
+
+public:
+ TwoPTwoCNewtonController()
+ {
+ this->setRelTolerance(1e-7);
+ this->setTargetSteps(10);
+ this->setMaxSteps(18);
+ };
+
+
+ /*!
+ * \brief
+ * Suggest a new time step size based either on the number of newton
+ * iterations required or on the variable switch
+ *
+ * \param u The current global solution vector
+ * \param uLastIter The previous global solution vector
+ *
+ */
+ void newtonEndStep(SolutionVector &uCurrentIter,
+ const SolutionVector &uLastIter)
+ {
+ // call the method of the base class
+ this->method().model().updateStaticData(uCurrentIter, uLastIter);
+ ParentType::newtonEndStep(uCurrentIter, uLastIter);
+ }
+
+ /*!
+ * \brief Update the current solution function with a delta vector.
+ *
+ * The error estimates required for the newtonConverged() and
+ * newtonProceed() methods should be updated here.
+ *
+ * Different update strategies, such as line search and chopped
+ * updates can be implemented. The default behaviour is just to
+ * subtract deltaU from uLastIter.
+ *
+ * \param uCurrentIter The solution after the current Newton iteration
+ * \param uLastIter The solution after the last Newton iteration
+ * \param deltaU The delta as calculated from solving the linear
+ * system of equations. This parameter also stores
+ * the updated solution.
+ */
+ void newtonUpdate(SolutionVector &uCurrentIter,
+ const SolutionVector &uLastIter,
+ const SolutionVector &deltaU)
+ {
+ this->writeConvergence_(uLastIter, deltaU);
+
+ this->newtonUpdateRelError(uLastIter, deltaU);
+
+ // compute the vertex and element colors for partial
+ // reassembly
+ if (enablePartialReassemble) {
+ Scalar reassembleTol = Dumux::geometricMean(this->error_, 0.1*this->tolerance_);
+ reassembleTol = std::max(reassembleTol, 0.1*this->tolerance_);
+ this->model_().jacobianAssembler().updateDiscrepancy(uLastIter, deltaU);
+ this->model_().jacobianAssembler().computeColors(reassembleTol);
+ }
+
+ if (GET_PROP_VALUE(TypeTag, PTAG(NewtonUseLineSearch)))
+ lineSearchUpdate_(uCurrentIter, uLastIter, deltaU);
+ else {
+ uCurrentIter = uLastIter;
+ uCurrentIter -= deltaU;
+ }
+ }
+
+
+ /*!
+ * \brief
+ * Returns true if the current solution can be considered to
+ * be accurate enough
+ */
+ bool newtonConverged()
+ {
+ if (this->method().model().switched())
+ return false;
+
+ return ParentType::newtonConverged();
+ };
+
+private:
+ void lineSearchUpdate_(SolutionVector &uCurrentIter,
+ const SolutionVector &uLastIter,
+ const SolutionVector &deltaU)
+ {
+ Scalar lambda = 1.0;
+ Scalar globDef;
+
+ // calculate the residual of the current solution
+ SolutionVector tmp(uLastIter);
+ Scalar oldGlobDef = this->method().model().globalResidual(tmp, uLastIter);
+
+ while (true) {
+ uCurrentIter = deltaU;
+ uCurrentIter *= -lambda;
+ uCurrentIter += uLastIter;
+
+ // calculate the residual of the current solution
+ globDef = this->method().model().globalResidual(tmp, uCurrentIter);
+
+ if (globDef < oldGlobDef || lambda <= 1.0/8) {
+ this->endIterMsg() << ", defect " << oldGlobDef << "->" << globDef << "@lambda=" << lambda;
+ return;
+ }
+
+ // try with a smaller update
+ lambda /= 2;
+ }
+ };
+
+};
+}
+
+#endif
diff --git a/dumux/boxmodels/new_2p2c/2p2cproblem.hh b/dumux/boxmodels/new_2p2c/2p2cproblem.hh
new file mode 100644
index 0000000000..53c91d2820
--- /dev/null
+++ b/dumux/boxmodels/new_2p2c/2p2cproblem.hh
@@ -0,0 +1,129 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Base class for all problems which use the two-phase,
+ * two-component box model
+ */
+#ifndef DUMUX_2P2C_PROBLEM_HH
+#define DUMUX_2P2C_PROBLEM_HH
+
+#include <dumux/boxmodels/common/boxproblem.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPTwoCModel
+ * \brief Base class for all problems which use the two-phase, two-component box model
+ *
+ */
+template<class TypeTag>
+class TwoPTwoCProblem : public BoxProblem<TypeTag>
+{
+ typedef BoxProblem<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GridView::Grid Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
+
+ // material properties
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+
+ enum {
+ dim = Grid::dimension,
+ dimWorld = Grid::dimensionworld
+ };
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param timeManager The time manager
+ * \param gridView The grid view
+ */
+ TwoPTwoCProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView),
+ gravity_(0),
+ spatialParams_(gridView)
+ {
+ if (GET_PROP_VALUE(TypeTag, PTAG(EnableGravity)))
+ gravity_[dim-1] = -9.81;
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This method MUST be overwritten by the actual problem.
+ */
+ Scalar temperature() const
+ { DUNE_THROW(Dune::NotImplemented, "The Problem must implement a temperature() method for isothermal problems!"); };
+
+ /*!
+ * \brief Returns the acceleration due to gravity.
+ *
+ * If the <tt>EnableGravity</tt> property is true, this means
+ * \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$, else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$
+ */
+ const GlobalPosition &gravity() const
+ { return gravity_; }
+
+ /*!
+ * \brief Returns the spatial parameters object.
+ */
+ SpatialParameters &spatialParameters()
+ { return spatialParams_; }
+
+ /*!
+ * \brief Returns the spatial parameters object.
+ */
+ const SpatialParameters &spatialParameters() const
+ { return spatialParams_; }
+
+ // \}
+
+private:
+ //! Returns the implementation of the problem (i.e. static polymorphism)
+ Implementation *asImp_()
+ { return static_cast<Implementation *>(this); }
+
+ //! \copydoc asImp_()
+ const Implementation *asImp_() const
+ { return static_cast<const Implementation *>(this); }
+
+ GlobalPosition gravity_;
+
+ // spatial parameters
+ SpatialParameters spatialParams_;
+};
+
+}
+
+#endif
diff --git a/dumux/boxmodels/new_2p2c/2p2cproperties.hh b/dumux/boxmodels/new_2p2c/2p2cproperties.hh
new file mode 100644
index 0000000000..5204442601
--- /dev/null
+++ b/dumux/boxmodels/new_2p2c/2p2cproperties.hh
@@ -0,0 +1,68 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Melanie Darcis *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup TwoPTwoCModel
+ */
+/*!
+ * \file
+ *
+ * \brief Defines the properties required for the 2p2c BOX model.
+ */
+#ifndef DUMUX_2P2C_PROPERTIES_HH
+#define DUMUX_2P2C_PROPERTIES_HH
+
+#include <dumux/boxmodels/common/boxproperties.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the isothermal single phase problems
+NEW_TYPE_TAG(BoxTwoPTwoC, INHERITS_FROM(BoxModel));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system
+NEW_PROP_TAG(NumComponents); //!< Number of fluid components in the system
+NEW_PROP_TAG(TwoPTwoCIndices); //!< Enumerations for the 2p2c models
+NEW_PROP_TAG(Formulation); //!< The formulation of the model
+NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters
+NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
+
+NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (extracted from the spatial parameters)
+NEW_PROP_TAG(MaterialLawParams); //!< The parameters of the material law (extracted from the spatial parameters)
+
+NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
+NEW_PROP_TAG(MobilityUpwindAlpha); //!< The value of the upwind parameter for the mobility
+}
+}
+
+#endif
diff --git a/dumux/boxmodels/new_2p2c/2p2cpropertydefaults.hh b/dumux/boxmodels/new_2p2c/2p2cpropertydefaults.hh
new file mode 100644
index 0000000000..a499883722
--- /dev/null
+++ b/dumux/boxmodels/new_2p2c/2p2cpropertydefaults.hh
@@ -0,0 +1,165 @@
+// $Id: 2p2cproperties.hh 3784 2010-06-24 13:43:57Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2008 by Klaus Mosthaf, Andreas Lauser, Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup TwoPTwoCModel
+ */
+/*!
+ * \file
+ *
+ * \brief Defines default values for most properties required by the
+ * 2p2c box model.
+ */
+#ifndef DUMUX_2P2C_PROPERTY_DEFAULTS_HH
+#define DUMUX_2P2C_PROPERTY_DEFAULTS_HH
+
+#include "2p2cindices.hh"
+
+#include "2p2cmodel.hh"
+#include "2p2cproblem.hh"
+#include "2p2cindices.hh"
+#include "2p2cfluxvariables.hh"
+#include "2p2cvolumevariables.hh"
+#include "2p2cfluidstate.hh"
+#include "2p2cproperties.hh"
+#include "2p2cnewtoncontroller.hh"
+#include "2p2cboundaryvariables.hh"
+
+namespace Dumux
+{
+
+namespace Properties {
+//////////////////////////////////////////////////////////////////
+// Property values
+//////////////////////////////////////////////////////////////////
+
+/*!
+ * \brief Set the property for the number of components.
+ *
+ * We just forward the number from the fluid system and use an static
+ * assert to make sure it is 2.
+ */
+SET_PROP(BoxTwoPTwoC, NumComponents)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ static const int value = FluidSystem::numComponents;
+
+ static_assert(value == 2,
+ "Only fluid systems with 2 components are supported by the 2p-2c model!");
+};
+
+/*!
+ * \brief Set the property for the number of fluid phases.
+ *
+ * We just forward the number from the fluid system and use an static
+ * assert to make sure it is 2.
+ */
+SET_PROP(BoxTwoPTwoC, NumPhases)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ static const int value = FluidSystem::numPhases;
+ static_assert(value == 2,
+ "Only fluid systems with 2 phases are supported by the 2p-2c model!");
+};
+
+SET_INT_PROP(BoxTwoPTwoC, NumEq, 2); //!< set the number of equations to 2
+
+//! Set the default formulation to pl-Sg
+SET_INT_PROP(BoxTwoPTwoC,
+ Formulation,
+ TwoPTwoCFormulation::plSg);
+
+/*!
+ * \brief Set the property for the material law by retrieving it from
+ * the spatial parameters.
+ */
+SET_PROP(BoxTwoPTwoC, MaterialLaw)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+
+public:
+ typedef typename SpatialParameters::MaterialLaw type;
+};
+
+/*!
+ * \brief Set the property for the material parameters by extracting
+ * it from the material law.
+ */
+SET_PROP(BoxTwoPTwoC, MaterialLawParams)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
+
+public:
+ typedef typename MaterialLaw::Params type;
+};
+
+//! Use the 2p2c local jacobian operator for the 2p2c model
+SET_TYPE_PROP(BoxTwoPTwoC,
+ LocalResidual,
+ TwoPTwoCLocalResidual<TypeTag>);
+
+//! Use the 2p2c specific newton controller for the 2p2c model
+SET_TYPE_PROP(BoxTwoPTwoC, NewtonController, TwoPTwoCNewtonController<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(BoxTwoPTwoC, Model, TwoPTwoCModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(BoxTwoPTwoC, VolumeVariables, TwoPTwoCVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+SET_TYPE_PROP(BoxTwoPTwoC, FluxVariables, TwoPTwoCFluxVariables<TypeTag>);
+
+//! the BoundaryVariables property
+SET_TYPE_PROP(BoxTwoPTwoC, BoundaryVariables, TwoPTwoCBoundaryVariables<TypeTag>);
+
+//! the upwind factor for the mobility.
+SET_SCALAR_PROP(BoxTwoPTwoC, MobilityUpwindAlpha, 1.0);
+
+//! The indices required by the isothermal 2p2c model
+SET_PROP(BoxTwoPTwoC,
+ TwoPTwoCIndices)
+{ private:
+ enum { Formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation)) };
+public:
+ typedef TwoPTwoCIndices<TypeTag, Formulation, 0> type;
+};
+
+// enable jacobian matrix recycling by default
+SET_BOOL_PROP(BoxTwoPTwoC, EnableJacobianRecycling, true);
+// enable partial reassembling by default
+SET_BOOL_PROP(BoxTwoPTwoC, EnablePartialReassemble, true);
+// enable time-step ramp up by default
+SET_BOOL_PROP(BoxTwoPTwoC, EnableTimeStepRampUp, false);
+
+//
+}
+
+}
+
+#endif
diff --git a/dumux/boxmodels/new_2p2c/2p2cvolumevariables.hh b/dumux/boxmodels/new_2p2c/2p2cvolumevariables.hh
new file mode 100644
index 0000000000..239e5ad450
--- /dev/null
+++ b/dumux/boxmodels/new_2p2c/2p2cvolumevariables.hh
@@ -0,0 +1,456 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008,2009 by Klaus Mosthaf, *
+ * Andreas Lauser, *
+ * Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Contains the quantities which are constant within a
+ * finite volume in the two-phase, two-component model.
+ */
+#ifndef DUMUX_2P2C_VOLUME_VARIABLES_HH
+#define DUMUX_2P2C_VOLUME_VARIABLES_HH
+
+#include <dumux/boxmodels/common/boxmodel.hh>
+#include <dumux/common/math.hh>
+
+#include <dune/common/collectivecommunication.hh>
+
+#include <dumux/material/fluidstates/equilibriumfluidstate.hh>
+#include <dumux/material/constraintsolvers/compositionfromfugacities.hh>
+
+#include <vector>
+#include <iostream>
+
+#include "2p2cproperties.hh"
+#include "2p2cindices.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPTwoCModel
+ * \brief Contains the quantities which are are constant within a
+ * finite volume in the two-phase, two-component model.
+ */
+template <class TypeTag>
+class TwoPTwoCVolumeVariables : public BoxVolumeVariables<TypeTag>
+{
+ typedef BoxVolumeVariables<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLawParams)) MaterialLawParams;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+ enum {
+ dim = GridView::dimension,
+
+ numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
+ formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation)),
+
+ // component indices
+ lCompIdx = Indices::lCompIdx,
+ gCompIdx = Indices::gCompIdx,
+
+ // phase indices
+ lPhaseIdx = Indices::lPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+
+ // primary variable indices
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx,
+
+ // phase states
+ lPhaseOnly = Indices::lPhaseOnly,
+ gPhaseOnly = Indices::gPhaseOnly,
+ bothPhases = Indices::bothPhases,
+
+ // formulations
+ plSg = TwoPTwoCFormulation::plSg,
+ pgSl = TwoPTwoCFormulation::pgSl,
+ };
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+public:
+ //! The return type of the fluidState() method
+ typedef Dumux::EquilibriumFluidState<Scalar, FluidSystem> FluidState;
+
+ /*!
+ * \brief Update all quantities for a given control volume.
+ *
+ * \param priVars The primary variables
+ * \param problem The problem
+ * \param element The element
+ * \param elemGeom The finite-volume geometry in the box scheme
+ * \param scvIdx The local index of the SCV (sub-control volume)
+ * \param isOldSol Evaluate function with solution of current or previous time step
+ */
+ void update(const PrimaryVariables &priVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ bool isOldSol)
+ {
+ ParentType::update(priVars,
+ problem,
+ element,
+ elemGeom,
+ scvIdx,
+ isOldSol);
+
+ typename FluidSystem::MutableParameters mutParams;
+
+ int globalVertIdx = problem.model().vertexMapper().map(element, scvIdx, dim);
+ int phasePresence = problem.model().phasePresence(globalVertIdx, isOldSol);
+
+ /////////////
+ // set the phase saturations
+ /////////////
+ if (phasePresence == gPhaseOnly) {
+ mutParams.setSaturation(lPhaseIdx, 0.0);
+ mutParams.setSaturation(gPhaseIdx, 1.0);
+ }
+ else if (phasePresence == lPhaseOnly) {
+ mutParams.setSaturation(lPhaseIdx, 1.0);
+ mutParams.setSaturation(gPhaseIdx, 0.0);
+ }
+ else if (phasePresence == bothPhases) {
+ Scalar Sg;
+ if (formulation == plSg)
+ Sg = priVars[switchIdx];
+ else if (formulation == pgSl)
+ Sg = 1.0 - priVars[switchIdx];
+
+ mutParams.setSaturation(lPhaseIdx, 1 - Sg);
+ mutParams.setSaturation(gPhaseIdx, Sg);
+ }
+
+ /////////////
+ // set the phase temperatures
+ /////////////
+ // update the temperature part of the energy module
+ Scalar T = asImp_().getTemperature(priVars,
+ element,
+ elemGeom,
+ scvIdx,
+ problem);
+ Valgrind::CheckDefined(T);
+ for (int i = 0; i < numPhases; ++i)
+ mutParams.setTemperature(i, T);
+
+ /////////////
+ // set the phase pressures
+ /////////////
+
+ // capillary pressure parameters
+ const MaterialLawParams &materialParams =
+ problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
+ Scalar pC = MaterialLaw::pC(materialParams, mutParams.saturation(lPhaseIdx));
+ if (formulation == plSg) {
+ mutParams.setPressure(lPhaseIdx, priVars[pressureIdx]);
+ mutParams.setPressure(gPhaseIdx, priVars[pressureIdx] + pC);
+ }
+ else if (formulation == pgSl){
+ mutParams.setPressure(lPhaseIdx, priVars[pressureIdx] - pC);
+ mutParams.setPressure(gPhaseIdx, priVars[pressureIdx]);
+ }
+ Valgrind::CheckDefined(mutParams.pressure(lPhaseIdx));
+ Valgrind::CheckDefined(mutParams.pressure(gPhaseIdx));
+
+ /////////////
+ // set the phase compositions.
+ /////////////
+ if (phasePresence == gPhaseOnly) {
+ // mass fractions
+ Scalar Xg1 = priVars[switchIdx];
+ Scalar Xg2 = 1 - Xg1;
+
+ // molar masses
+ Scalar M1 = FluidSystem::molarMass(lCompIdx);
+ Scalar M2 = FluidSystem::molarMass(gCompIdx);
+
+ // convert mass to mole fractions
+ Scalar meanM = M1*M2/(M2 + Xg2*(M1 - M2));
+ mutParams.setMoleFrac(gPhaseIdx, lCompIdx, Xg1 * M1/meanM);
+ mutParams.setMoleFrac(gPhaseIdx, gCompIdx, Xg2 * M2/meanM);
+
+ // calculate component fugacities in gas phase
+ Scalar fug1 = FluidSystem::computeFugacity(mutParams, gPhaseIdx, lCompIdx);
+ Scalar fug2 = FluidSystem::computeFugacity(mutParams, gPhaseIdx, gCompIdx);
+ mutParams.setFugacity(gPhaseIdx, lCompIdx, fug1);
+ mutParams.setFugacity(gPhaseIdx, gCompIdx, fug2);
+
+ // use same fugacities in liquid phase
+ mutParams.setFugacity(lPhaseIdx, lCompIdx, fug1);
+ mutParams.setFugacity(lPhaseIdx, gCompIdx, fug2);
+
+ // initial guess of liquid composition
+ mutParams.setMoleFrac(lPhaseIdx, lCompIdx, 0.98);
+ mutParams.setMoleFrac(lPhaseIdx, gCompIdx, 0.02);
+
+ // calculate liquid composition from fugacities
+ typedef Dumux::CompositionFromFugacities<Scalar, FluidSystem> CompFromFug;
+ CompFromFug::solve(mutParams, lPhaseIdx);
+
+ Valgrind::CheckDefined(mutParams.moleFrac(gPhaseIdx, lCompIdx));
+ Valgrind::CheckDefined(mutParams.moleFrac(gPhaseIdx, gCompIdx));
+ Valgrind::CheckDefined(mutParams.moleFrac(lPhaseIdx, lCompIdx));
+ Valgrind::CheckDefined(mutParams.moleFrac(lPhaseIdx, gCompIdx));
+
+ // calculate molar volume of gas phase
+ Scalar Vmg = FluidSystem::computeMolarVolume(mutParams, gPhaseIdx);
+ mutParams.setMolarVolume(gPhaseIdx, Vmg);
+ }
+ else if (phasePresence == lPhaseOnly) {
+ // mass fractions
+ Scalar Xl2 = priVars[switchIdx];
+ Scalar Xl1 = 1 - Xl2;
+
+ // molar masses
+ Scalar M1 = FluidSystem::molarMass(lCompIdx);
+ Scalar M2 = FluidSystem::molarMass(gCompIdx);
+
+ // convert mass to mole fractions
+ Scalar meanM = M1*M2/(M2 + Xl2*(M1 - M2));
+ mutParams.setMoleFrac(lPhaseIdx, lCompIdx, Xl1 * M1/meanM);
+ mutParams.setMoleFrac(lPhaseIdx, gCompIdx, Xl2 * M2/meanM);
+
+ // calculate component fugacities in liquid phase
+ Scalar fug1 = FluidSystem::computeFugacity(mutParams, lPhaseIdx, lCompIdx);
+ Scalar fug2 = FluidSystem::computeFugacity(mutParams, lPhaseIdx, gCompIdx);
+ mutParams.setFugacity(lPhaseIdx, lCompIdx, fug1);
+ mutParams.setFugacity(lPhaseIdx, gCompIdx, fug2);
+
+ // use same fugacities in gas phase
+ mutParams.setFugacity(gPhaseIdx, lCompIdx, fug1);
+ mutParams.setFugacity(gPhaseIdx, gCompIdx, fug2);
+
+ // initial guess of gas composition
+ mutParams.setMoleFrac(gPhaseIdx, lCompIdx, 0.05);
+ mutParams.setMoleFrac(gPhaseIdx, gCompIdx, 0.95);
+
+ // calculate liquid composition from fugacities
+ typedef Dumux::CompositionFromFugacities<Scalar, FluidSystem> CompFromFug;
+ CompFromFug::solve(mutParams, gPhaseIdx);
+ Valgrind::CheckDefined(mutParams.moleFrac(gPhaseIdx, lCompIdx));
+ Valgrind::CheckDefined(mutParams.moleFrac(gPhaseIdx, gCompIdx));
+ Valgrind::CheckDefined(mutParams.moleFrac(lPhaseIdx, lCompIdx));
+ Valgrind::CheckDefined(mutParams.moleFrac(lPhaseIdx, gCompIdx));
+
+ // calculate molar volume of liquid phase
+ Scalar Vml = FluidSystem::computeMolarVolume(mutParams, lPhaseIdx);
+ mutParams.setMolarVolume(lPhaseIdx, Vml);
+ }
+ else if (phasePresence == bothPhases) {
+ // HACK: assume both phases to be an ideal mixture,
+ // i.e. the fugacity coefficents do not depend on the
+ // composition
+ Scalar phi_l1 = FluidSystem::computeFugacityCoeff(mutParams, lPhaseIdx, lCompIdx);
+ Scalar phi_l2 = FluidSystem::computeFugacityCoeff(mutParams, lPhaseIdx, gCompIdx);
+ Scalar phi_g1 = FluidSystem::computeFugacityCoeff(mutParams, gPhaseIdx, lCompIdx);
+ Scalar phi_g2 = FluidSystem::computeFugacityCoeff(mutParams, gPhaseIdx, gCompIdx);
+ Scalar pl = mutParams.pressure(lPhaseIdx);
+ Scalar pg = mutParams.pressure(gPhaseIdx);
+ Valgrind::CheckDefined(phi_l1);
+ Valgrind::CheckDefined(phi_l2);
+ Valgrind::CheckDefined(phi_g1);
+ Valgrind::CheckDefined(phi_g2);
+
+ Scalar xg2 = (phi_g2/phi_l1 - pl/pg) / (phi_g1/phi_l1 - phi_g2/phi_l2);
+ Scalar xg1 = 1 - xg2;
+ Scalar xl2 = (xg2*pg*phi_g2)/(pl*phi_l2);
+ Scalar xl1 = 1 - xl2;
+
+ mutParams.setMoleFrac(lPhaseIdx, lCompIdx, xl1);
+ mutParams.setMoleFrac(lPhaseIdx, gCompIdx, xl2);
+ mutParams.setMoleFrac(gPhaseIdx, lCompIdx, xg1);
+ mutParams.setMoleFrac(gPhaseIdx, gCompIdx, xg2);
+
+ Scalar Vml = FluidSystem::computeMolarVolume(mutParams, lPhaseIdx);
+ Scalar Vmg = FluidSystem::computeMolarVolume(mutParams, gPhaseIdx);
+ mutParams.setMolarVolume(lPhaseIdx, Vml);
+ mutParams.setMolarVolume(gPhaseIdx, Vmg);
+ }
+
+
+ // Mobilities
+ Scalar muL = FluidSystem::computeViscosity(mutParams, lPhaseIdx);
+ Scalar krL = MaterialLaw::krw(materialParams, mutParams.saturation(lPhaseIdx));
+ mobility_[lPhaseIdx] = krL / muL;
+ Valgrind::CheckDefined(mobility_[lPhaseIdx]);
+
+ // ATTENTION: krn requires the liquid saturation as parameter!
+ Scalar muG = FluidSystem::computeViscosity(mutParams, gPhaseIdx);
+ Scalar krG = MaterialLaw::krn(materialParams, mutParams.saturation(lPhaseIdx));
+ mobility_[gPhaseIdx] = krG / muG;
+ Valgrind::CheckDefined(mobility_[gPhaseIdx]);
+
+#if 0
+ // binary diffusion coefficents
+ diffCoeff_[phaseIdx] =
+ FluidSystem::diffCoeff(phaseIdx,
+ lCompIdx,
+ gCompIdx,
+ fluidState_.temperature(),
+ fluidState_.phasePressure(phaseIdx),
+ fluidState_);
+ Valgrind::CheckDefined(diffCoeff_[phaseIdx]);
+#endif
+
+ // porosity
+ porosity_ = problem.spatialParameters().porosity(element,
+ elemGeom,
+ scvIdx);
+ Valgrind::CheckDefined(porosity_);
+
+ asImp_().updateEnergy(mutParams, priVars, element, elemGeom, scvIdx, problem);
+
+ // assign the equilibrium fluid state from the generic one
+ fluidState_.assign(mutParams);
+ }
+
+ /*!
+ * \brief Returns the phase state for the control-volume.
+ */
+ const FluidState &fluidState() const
+ { return fluidState_; }
+
+ /*!
+ * \brief Returns the effective saturation of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar saturation(int phaseIdx) const
+ { return fluidState_.saturation(phaseIdx); }
+
+ /*!
+ * \brief Returns the mass density of a given phase within the
+ * control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar density(int phaseIdx) const
+ { return fluidState_.density(phaseIdx); }
+
+ /*!
+ * \brief Returns the mass density of a given phase within the
+ * control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar molarDensity(int phaseIdx) const
+ { return fluidState_.molarDensity(phaseIdx); }
+
+ /*!
+ * \brief Returns the effective pressure of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar pressure(int phaseIdx) const
+ { return fluidState_.pressure(phaseIdx); }
+
+ /*!
+ * \brief Returns temperature inside the sub-control volume.
+ *
+ * Note that we assume thermodynamic equilibrium, i.e. the
+ * temperature of the rock matrix and of all fluid phases are
+ * identical.
+ */
+ Scalar temperature() const
+ { return fluidState_.temperature(); }
+
+ /*!
+ * \brief Returns the effective mobility of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar mobility(int phaseIdx) const
+ {
+ return mobility_[phaseIdx];
+ }
+
+ /*!
+ * \brief Returns the effective capillary pressure within the control volume.
+ */
+ Scalar capillaryPressure() const
+ { return fluidState_.pressure(gPhaseIdx) - fluidState_.pressure(lPhaseIdx); }
+
+ /*!
+ * \brief Returns the average porosity within the control volume.
+ */
+ Scalar porosity() const
+ { return porosity_; }
+
+#if 0
+ /*!
+ * \brief Returns the binary diffusion coefficients for a phase
+ */
+ Scalar diffCoeff(int phaseIdx) const
+ { return diffCoeff_[phaseIdx]; }
+#endif
+
+protected:
+ Scalar getTemperature(const PrimaryVariables &priVars,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ const Problem &problem)
+ {
+ return problem.temperature(element, elemGeom, scvIdx);
+ }
+
+ template <class MutableParameters>
+ void updateEnergy(MutableParameters &mutParams,
+ const PrimaryVariables &priVars,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ const Problem &problem)
+ {};
+
+ Scalar porosity_; //!< Effective porosity within the control volume
+ Scalar mobility_[numPhases]; //!< Effective mobility within the control volume
+// Scalar diffCoeff_[numPhases]; //!< Binary diffusion coefficients for the phases
+ FluidState fluidState_;
+
+private:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/boxmodels/new_2p2c/Makefile.am b/dumux/boxmodels/new_2p2c/Makefile.am
new file mode 100644
index 0000000000..f15c521752
--- /dev/null
+++ b/dumux/boxmodels/new_2p2c/Makefile.am
@@ -0,0 +1,4 @@
+2p2cdir = $(includedir)/dumux/boxmodels/2p2c
+2p2c_HEADERS = *.hh
+
+include $(top_srcdir)/am/global-rules
diff --git a/dumux/boxmodels/new_2p2cni/2p2cniboundaryvariables.hh b/dumux/boxmodels/new_2p2cni/2p2cniboundaryvariables.hh
new file mode 100644
index 0000000000..6856d2068c
--- /dev/null
+++ b/dumux/boxmodels/new_2p2cni/2p2cniboundaryvariables.hh
@@ -0,0 +1,272 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2010 by Katherina Baber, Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch, Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief This file contains the data which is required to calculate
+ * all fluxes of the fluid phase over the boundary of a finite volume.
+ *
+ * This means pressure and temperature gradients, phase densities at
+ * the integration point of the boundary, etc.
+ */
+#ifndef DUMUX_2P2CNI_BOUNDARY_VARIABLES_HH
+#define DUMUX_2P2CNI_BOUNDARY_VARIABLES_HH
+
+#include <dumux/common/math.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPTwoCNIModel
+ * \brief This template class contains the data which is required to
+ * calculate the fluxes of the fluid phases over the boundary of a
+ * finite volume for the 2p2cni model.
+ *
+ * This means pressure and velocity gradients, phase density and viscosity at
+ * the integration point of the boundary, etc.
+ */
+template <class TypeTag>
+class TwoPTwoCNIBoundaryVariables
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+
+ enum {
+ dim = GridView::dimension,
+ };
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef Dune::FieldVector<Scalar, dim> VelocityVector;
+ typedef Dune::FieldVector<Scalar, dim> ScalarGradient;
+ typedef Dune::FieldMatrix<Scalar, dim, dim> VectorGradient;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename FVElementGeometry::SubControlVolume SCV;
+ typedef typename FVElementGeometry::BoundaryFace BoundaryFace;
+
+ enum {
+ lPhaseIdx = Indices::lPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+
+ lCompIdx = Indices::lCompIdx,
+ gCompIdx = Indices::gCompIdx,
+
+ numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases))
+ };
+
+public:
+ TwoPTwoCNIBoundaryVariables(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int boundaryFaceIdx,
+ const ElementVolumeVariables &elemDat,
+ int scvIdx)
+ : fvElemGeom_(elemGeom), scvIdx_(scvIdx)
+ {
+ boundaryFace_ = &fvElemGeom_.boundaryFace[boundaryFaceIdx];
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ densityAtIP_[phaseIdx] = Scalar(0);
+ molarDensityAtIP_[phaseIdx] = Scalar(0);
+ potentialGrad_[phaseIdx] = Scalar(0);
+ concentrationGrad_[phaseIdx] = Scalar(0);
+ molarConcGrad_[phaseIdx] = Scalar(0);
+ }
+
+ calculateBoundaryValues_(problem, element, elemDat);
+ };
+
+ Scalar KmvpNormal(int phaseIdx) const
+ { return KmvpNormal_[phaseIdx]; }
+
+ /*!
+ * \brief The binary diffusion coefficient for each fluid phase.
+ */
+ Scalar porousDiffCoeff(int phaseIdx) const
+ { return porousDiffCoeff_[phaseIdx]; };
+
+ /*!
+ * \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase at the integration
+ * point.
+ */
+ Scalar densityAtIP(int phaseIdx) const
+ { return densityAtIP_[phaseIdx]; }
+
+ /*!
+ * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ of a phase at the integration
+ * point.
+ */
+ Scalar molarDensityAtIP(int phaseIdx) const
+ { return molarDensityAtIP_[phaseIdx]; }
+
+ /*!
+ * \brief The concentration gradient of a component in a phase.
+ */
+ const ScalarGradient &concentrationGrad(int phaseIdx) const
+ { return concentrationGrad_[phaseIdx]; };
+
+ /*!
+ * \brief The molar concentration gradient of a component in a phase.
+ */
+ const ScalarGradient &molarConcGrad(int phaseIdx) const
+ { return molarConcGrad_[phaseIdx]; };
+
+ /*!
+ * \brief The total heat flux \f$\mathrm{[J/s]}\f$ due to heat conduction
+ * of the rock matrix over the sub-control volume's face in
+ * direction of the face normal.
+ */
+ Scalar normalMatrixHeatFlux() const
+ { return normalMatrixHeatFlux_; }
+
+ const BoundaryFace& boundaryFace() const
+ { return *boundaryFace_; }
+
+protected:
+ void calculateBoundaryValues_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemDat)
+ {
+ ScalarGradient tmp(0.0);
+
+ // calculate gradients and secondary variables at IPs of the boundary
+ for (int idx = 0;
+ idx < fvElemGeom_.numVertices;
+ idx++) // loop over adjacent vertices
+ {
+ // FE gradient at vertex idx
+ const ScalarGradient& feGrad = boundaryFace_->grad[idx];
+
+ // compute sum of pressure gradients for each phase
+ for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
+ {
+ // the pressure gradient
+ tmp = feGrad;
+ tmp *= elemDat[idx].pressure(phaseIdx);
+ potentialGrad_[phaseIdx] += tmp;
+ }
+
+ // the concentration gradient of the non-wetting
+ // component in the wetting phase
+ tmp = feGrad;
+ tmp *= elemDat[idx].fluidState().massFrac(lPhaseIdx, gCompIdx);
+ concentrationGrad_[lPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemDat[idx].fluidState().moleFrac(lPhaseIdx, gCompIdx);
+ molarConcGrad_[lPhaseIdx] += tmp;
+
+ // the concentration gradient of the wetting component
+ // in the non-wetting phase
+ tmp = feGrad;
+ tmp *= elemDat[idx].fluidState().massFrac(gPhaseIdx, lCompIdx);
+ concentrationGrad_[gPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemDat[idx].fluidState().moleFrac(gPhaseIdx, lCompIdx);
+ molarConcGrad_[gPhaseIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemDat[idx].temperature();
+ temperatureGrad_ += tmp;
+
+ for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
+ {
+ densityAtIP_[phaseIdx] += elemDat[idx].density(phaseIdx)*boundaryFace_->shapeValue[idx];
+ molarDensityAtIP_[phaseIdx] += elemDat[idx].molarDensity(phaseIdx)*boundaryFace_->shapeValue[idx];
+ }
+ }
+
+ for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
+ {
+ tmp = problem.gravity();
+ tmp *= densityAtIP_[phaseIdx];
+
+ potentialGrad_[phaseIdx] -= tmp;
+
+ Scalar k = problem.spatialParameters().intrinsicPermeability(element, fvElemGeom_, scvIdx_);
+ VectorGradient K(0);
+ K[0][0] = K[1][1] = k;
+ ScalarGradient Kmvp;
+ K.mv(potentialGrad_[phaseIdx], Kmvp);
+ KmvpNormal_[phaseIdx] = - (Kmvp*boundaryFace_->normal);
+
+ const VolumeVariables &vertDat = elemDat[scvIdx_];
+
+ if (vertDat.saturation(phaseIdx) <= 0)
+ porousDiffCoeff_[phaseIdx] = 0.0;
+ else
+ {
+ Scalar tau = 1.0/(vertDat.porosity()*vertDat.porosity())*
+ pow(vertDat.porosity()*vertDat.saturation(phaseIdx), 7.0/3);
+
+ porousDiffCoeff_[phaseIdx] = vertDat.porosity()*vertDat.saturation(phaseIdx)*tau*vertDat.diffCoeff(phaseIdx);
+ }
+ }
+
+ // The spatial parameters calculates the actual heat flux vector
+ problem.spatialParameters().boundaryMatrixHeatFlux(tmp,
+ elemDat,
+ temperatureGrad_,
+ element,
+ fvElemGeom_,
+ scvIdx_);
+ // project the heat flux vector on the face's normal vector
+ normalMatrixHeatFlux_ = tmp*boundaryFace_->normal;
+
+ }
+
+ const FVElementGeometry &fvElemGeom_;
+ const BoundaryFace *boundaryFace_;
+
+ // gradients
+ ScalarGradient potentialGrad_[numPhases];
+ ScalarGradient concentrationGrad_[numPhases];
+ ScalarGradient molarConcGrad_[numPhases];
+
+ // density of each face at the integration point
+ Scalar densityAtIP_[numPhases];
+ Scalar molarDensityAtIP_[numPhases];
+
+ // intrinsic permeability times pressure potential gradient
+ // projected on the face normal
+ Scalar KmvpNormal_[numPhases];
+
+ // the diffusion coefficient for the porous medium
+ Scalar porousDiffCoeff_[numPhases];
+
+ ScalarGradient temperatureGrad_;
+ Scalar normalMatrixHeatFlux_;
+
+ int scvIdx_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/boxmodels/new_2p2cni/2p2cnifluxvariables.hh b/dumux/boxmodels/new_2p2cni/2p2cnifluxvariables.hh
new file mode 100644
index 0000000000..df52333574
--- /dev/null
+++ b/dumux/boxmodels/new_2p2cni/2p2cnifluxvariables.hh
@@ -0,0 +1,131 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Melanie Darcis *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief This file contains the data which is required to calculate
+ * all fluxes (mass of components and energy) over a face of a finite volume.
+ *
+ * This means pressure, concentration and temperature gradients, phase
+ * densities at the integration point, etc.
+ */
+#ifndef DUMUX_2P2CNI_FLUX_VARIABLES_HH
+#define DUMUX_2P2CNI_FLUX_VARIABLES_HH
+
+#include <dumux/common/math.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPTwoCNIModel
+ * \brief This template class contains the data which is required to
+ * calculate all fluxes (mass of components and energy) over a face of a finite
+ * volume for the non-isothermal two-phase, two-component model.
+ *
+ * This means pressure and concentration gradients, phase densities at
+ * the integration point, etc.
+ */
+template <class TypeTag>
+class TwoPTwoCNIFluxVariables : public TwoPTwoCFluxVariables<TypeTag>
+{
+ typedef TwoPTwoCFluxVariables<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
+
+ typedef typename GridView::ctype CoordScalar;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename FVElementGeometry::SubControlVolume SCV;
+ typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
+
+ typedef Dune::FieldVector<CoordScalar, dimWorld> Vector;
+
+public:
+ /*
+ * \brief The constructor
+ *
+ * \param problem The problem
+ * \param element The finite element
+ * \param elemGeom The finite-volume geometry in the box scheme
+ * \param faceIdx The local index of the SCV (sub-control-volume) face
+ * \param elemDat The volume variables of the current element
+ */
+ TwoPTwoCNIFluxVariables(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvfIdx,
+ const ElementVolumeVariables &elemDat)
+ : ParentType(problem, element, elemGeom, scvfIdx, elemDat)
+ {
+ // calculate temperature gradient using finite element
+ // gradients
+ Vector temperatureGrad(0);
+ Vector tmp(0.0);
+ for (int vertIdx = 0; vertIdx < elemGeom.numVertices; vertIdx++)
+ {
+ tmp = elemGeom.subContVolFace[scvfIdx].grad[vertIdx];
+ tmp *= elemDat[vertIdx].temperature();
+ temperatureGrad += tmp;
+ }
+
+ // The spatial parameters calculates the actual heat flux vector
+ problem.spatialParameters().matrixHeatFlux(tmp,
+ *this,
+ elemDat,
+ temperatureGrad,
+ element,
+ elemGeom,
+ scvfIdx);
+ // project the heat flux vector on the face's normal vector
+ normalMatrixHeatFlux_ = tmp*elemGeom.subContVolFace[scvfIdx].normal;
+ }
+
+ /*!
+ * \brief The total heat flux \f$\mathrm{[J/s]}\f$ due to heat conduction
+ * of the rock matrix over the sub-control volume's face in
+ * direction of the face normal.
+ */
+ Scalar normalMatrixHeatFlux() const
+ { return normalMatrixHeatFlux_; }
+
+private:
+ Scalar normalMatrixHeatFlux_;
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/boxmodels/new_2p2cni/2p2cniindices.hh b/dumux/boxmodels/new_2p2cni/2p2cniindices.hh
new file mode 100644
index 0000000000..e912f6503d
--- /dev/null
+++ b/dumux/boxmodels/new_2p2cni/2p2cniindices.hh
@@ -0,0 +1,57 @@
+// $Id: 2p2cniproperties.hh 3784 2010-06-24 13:43:57Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Melanie Darcis *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+
+/*!
+ * \file
+ *
+ * \brief Defines the indices used by the 2p2cni box model
+ */
+#ifndef DUMUX_2P2CNI_INDICES_HH
+#define DUMUX_2P2CNI_INDICES_HH
+
+#include <dumux/boxmodels/new_2p2c/2p2cindices.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPTwoCNIModel
+ */
+// \{
+
+/*!
+ * \brief Enumerations for the non-isothermal 2-phase 2-component model
+ *
+ * \tparam formulation The formulation, either pwSn or pnSw.
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag, int formulation, int PVOffset>
+class TwoPTwoCNIIndices : public TwoPTwoCIndices<TypeTag, formulation, PVOffset>
+{
+public:
+ static const int temperatureIdx = PVOffset + 2; //! The index for temperature in primary variable vectors.
+ static const int energyEqIdx = PVOffset + 2; //! The index for energy in equation vectors.
+};
+
+}
+#endif
diff --git a/dumux/boxmodels/new_2p2cni/2p2cnilocalresidual.hh b/dumux/boxmodels/new_2p2cni/2p2cnilocalresidual.hh
new file mode 100644
index 0000000000..7bed80375e
--- /dev/null
+++ b/dumux/boxmodels/new_2p2cni/2p2cnilocalresidual.hh
@@ -0,0 +1,183 @@
+// $Id: 2p2cnilocalresidual.hh 3784 2010-06-24 13:43:57Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Melanie Darcis *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the non-isothermal two-phase two-component box model.
+ *
+ */
+#ifndef DUMUX_NEW_2P2CNI_LOCAL_RESIDUAL_HH
+#define DUMUX_NEW_2P2CNI_LOCAL_RESIDUAL_HH
+
+#include <dumux/boxmodels/new_2p2c/2p2clocalresidual.hh>
+
+#include "2p2cnivolumevariables.hh"
+#include "2p2cnifluxvariables.hh"
+#include "2p2cniproperties.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPTwoCNIModel
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the two-phase two-component box model.
+ */
+template<class TypeTag>
+class TwoPTwoCNILocalResidual : public TwoPTwoCLocalResidual<TypeTag>
+{
+ typedef TwoPTwoCNILocalResidual<TypeTag> ThisType;
+ typedef TwoPTwoCLocalResidual<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
+
+ energyEqIdx = Indices::energyEqIdx,
+ temperatureIdx = Indices::temperatureIdx,
+
+ lPhaseIdx = Indices::lPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx
+ };
+
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ static const Scalar mobilityUpwindAlpha = GET_PROP_VALUE(TypeTag, PTAG(MobilityUpwindAlpha));
+
+public:
+ /*!
+ * \brief Evaluate the amount all conservation quantities
+ * (e.g. phase mass) within a sub-control volume.
+ *
+ * The result should be averaged over the volume (e.g. phase mass
+ * inside a sub control volume divided by the volume)
+ *
+ * \param result The storage of the conservation quantitiy (mass or energy) within the sub-control volume
+ * \param scvIdx The SCV (sub-control-volume) index
+ * \param usePrevSol Evaluate function with solution of current or previous time step
+ */
+ void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const
+ {
+ // compute the storage term for phase mass
+ ParentType::computeStorage(result, scvIdx, usePrevSol);
+
+ // if flag usePrevSol is set, the solution from the previous
+ // time step is used, otherwise the current solution is
+ // used. The secondary variables are used accordingly. This
+ // is required to compute the derivative of the storage term
+ // using the implicit euler method.
+ const ElementVolumeVariables &elemDat = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
+ const VolumeVariables &vertDat = elemDat[scvIdx];
+
+ // compute the energy storage
+ result[energyEqIdx] =
+ vertDat.porosity()*(vertDat.density(lPhaseIdx) *
+ vertDat.fluidState().internalEnergy(lPhaseIdx) *
+ vertDat.saturation(lPhaseIdx)
+ +
+ vertDat.density(gPhaseIdx) *
+ vertDat.fluidState().internalEnergy(gPhaseIdx) *
+ vertDat.saturation(gPhaseIdx))
+ +
+ vertDat.temperature()*vertDat.heatCapacity();
+ }
+
+ /*!
+ * \brief Evaluates the advective mass flux and the heat flux
+ * over a face of a subcontrol volume and writes the result in
+ * the flux vector.
+ *
+ * \param flux The advective flux over the SCV (sub-control-volume) face for each component
+ * \param fluxData The flux variables at the current SCV face
+ *
+ * This method is called by compute flux (base class)
+ */
+ void computeAdvectiveFlux(PrimaryVariables &flux,
+ const FluxVariables &fluxData) const
+ {
+ // advective mass flux
+ ParentType::computeAdvectiveFlux(flux, fluxData);
+
+ // advective heat flux in all phases
+ flux[energyEqIdx] = 0;
+ for (int phase = 0; phase < numPhases; ++phase) {
+ // vertex data of the upstream and the downstream vertices
+ const VolumeVariables &up = this->curVolVars_(fluxData.upstreamIdx(phase));
+ const VolumeVariables &dn = this->curVolVars_(fluxData.downstreamIdx(phase));
+
+ flux[energyEqIdx] +=
+ fluxData.KmvpNormal(phase) * (
+ mobilityUpwindAlpha * // upstream vertex
+ ( up.density(phase) *
+ up.mobility(phase) *
+ up.fluidState().enthalpy(phase))
+ +
+ (1-mobilityUpwindAlpha) * // downstream vertex
+ ( dn.density(phase) *
+ dn.mobility(phase) *
+ dn.fluidState().enthalpy(phase)) );
+ }
+ }
+
+ /*!
+ * \brief Adds the diffusive heat flux to the flux vector over
+ * the face of a sub-control volume.
+ *
+ * \param flux The diffusive flux over the SCV (sub-control-volume) face for each conservation quantity (mass, energy)
+ * \param fluxData The flux variables at the current SCV face
+ *
+ * This method is called by compute flux (base class)
+ */
+ void computeDiffusiveFlux(PrimaryVariables &flux,
+ const FluxVariables &fluxData) const
+ {
+ // diffusive mass flux
+ ParentType::computeDiffusiveFlux(flux, fluxData);
+
+ // diffusive heat flux
+ flux[temperatureIdx] +=
+ fluxData.normalMatrixHeatFlux();
+ }
+};
+
+}
+
+#endif
diff --git a/dumux/boxmodels/new_2p2cni/2p2cnimodel.hh b/dumux/boxmodels/new_2p2cni/2p2cnimodel.hh
new file mode 100644
index 0000000000..73defa1c98
--- /dev/null
+++ b/dumux/boxmodels/new_2p2cni/2p2cnimodel.hh
@@ -0,0 +1,108 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Melanie Darcis *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Adaption of the BOX scheme to the non-isothermal two-phase two-component flow model.
+ */
+#ifndef DUMUX_NEW_2P2CNI_MODEL_HH
+#define DUMUX_NEW_2P2CNI_MODEL_HH
+
+#include <dumux/boxmodels/new_2p2c/2p2cmodel.hh>
+
+namespace Dumux {
+/*!
+ * \ingroup BoxModels
+ * \defgroup TwoPTwoCNIModel Non-isothermal two-phase two-component box model
+ */
+
+/*!
+ * \ingroup TwoPTwoCNIModel
+ * \brief Adaption of the BOX scheme to the non-isothermal two-phase two-component flow model.
+ *
+ * This model implements a non-isothermal two-phase flow of two compressible and partly miscible fluids
+ * \f$\alpha \in \{ w, n \}\f$. Thus each component \f$\kappa \{ w, a \}\f$ can be present in
+ * each phase.
+ * Using the standard multiphase Darcy approach a mass balance equation is
+ * solved:
+ * \f{eqnarray*}
+ && \phi \frac{\partial (\sum_\alpha \varrho_\alpha X_\alpha^\kappa S_\alpha )}{\partial t}
+ - \sum_\alpha \text{div} \left\{ \varrho_\alpha X_\alpha^\kappa
+ \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K}
+ (\text{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g}) \right\}\\
+ &-& \sum_\alpha \text{div} \left\{{\bf D_{\alpha, pm}^\kappa} \varrho_{\alpha} \text{grad}\, X^\kappa_{\alpha} \right\}
+ - \sum_\alpha q_\alpha^\kappa = 0 \qquad \kappa \in \{w, a\} \, ,
+ \alpha \in \{w, n\}
+ * \f}
+ * For the energy balance, local thermal equilibrium is assumed which results in one
+ * energy conservation equation for the porous solid matrix and the fluids:
+ * \f{eqnarray*}
+ && \phi \frac{\partial \left( \sum_\alpha \varrho_\alpha u_\alpha S_\alpha \right)}{\partial t}
+ + \left( 1 - \phi \right) \frac{\partial (\varrho_s c_s T)}{\partial t}
+ - \sum_\alpha \text{div} \left\{ \varrho_\alpha h_\alpha
+ \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K} \left( \text{grad}\,
+ p_\alpha
+ - \varrho_\alpha \mathbf{g} \right) \right\} \\
+ &-& \text{div} \left( \lambda_{pm} \text{grad} \, T \right)
+ - q^h = 0 \qquad \alpha \in \{w, n\}
+\f}
+ *
+ * This is discretized using a fully-coupled vertex
+ * centered finite volume (box) scheme as spatial and
+ * the implicit Euler method as temporal discretization.
+ *
+ * By using constitutive relations for the capillary pressure \f$p_c =
+ * p_n - p_w\f$ and relative permeability \f$k_{r\alpha}\f$ and taking
+ * advantage of the fact that \f$S_w + S_n = 1\f$ and \f$X^\kappa_w + X^\kappa_n = 1\f$, the number of
+ * unknowns can be reduced to two.
+ * If both phases are present the primary variables are, like in the nonisothermal two-phase model, either \f$p_w\f$, \f$S_n\f$ and
+ * temperature or \f$p_n\f$, \f$S_w\f$ and temperature. The formulation which ought to be used can be
+ * specified by setting the <tt>Formulation</tt> property to either
+ * <tt>TwoPTwoIndices::pWsN</tt> or <tt>TwoPTwoCIndices::pNsW</tt>. By
+ * default, the model uses \f$p_w\f$ and \f$S_n\f$.
+ * In case that only one phase (nonwetting or wetting phase) is present the second primary
+ * variable represents a mass fraction. The correct assignment of the second
+ * primary variable is performed by a phase state dependent primary variable switch. The phase state is stored for all nodes of the system. The following cases can be distinguished:
+ * <ul>
+ * <li>
+ * Both phases are present: The saturation is used (either\f$S_n\f$ or \f$S_w\f$, dependent on the chosen formulation).
+ * </li>
+ * <li>
+ * Only wetting phase is present: The mass fraction of air in the wetting phase \f$X^a_w\f$ is used.
+ * </li>
+ * <li>
+ * Only non-wetting phase is present: The mass fraction of water in the non-wetting phase, \f$X^w_n\f$, is used.
+ * </li>
+ * </ul>
+ */
+template<class TypeTag>
+class TwoPTwoCNIModel : public TwoPTwoCModel<TypeTag>
+{
+};
+
+}
+
+#include "2p2cnipropertydefaults.hh"
+
+#endif
diff --git a/dumux/boxmodels/new_2p2cni/2p2cniproblem.hh b/dumux/boxmodels/new_2p2cni/2p2cniproblem.hh
new file mode 100644
index 0000000000..703f3b6184
--- /dev/null
+++ b/dumux/boxmodels/new_2p2cni/2p2cniproblem.hh
@@ -0,0 +1,80 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Base class for all problems which use the non-isothermal two-phase,
+ * two-component box model
+ */
+#ifndef DUMUX_2P2CNI_PROBLEM_HH
+#define DUMUX_2P2CNI_PROBLEM_HH
+
+#include <dumux/boxmodels/new_2p2c/2p2cproblem.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPTwoCNIModel
+ * \brief Base class for all problems which use the non-isothermal
+ * two-phase, two-component box model.
+ */
+template<class TypeTag>
+class TwoPTwoCNIProblem : public TwoPTwoCProblem<TypeTag>
+{
+ typedef TwoPTwoCProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param timeManager The time manager
+ * \param gridView The grid view
+ */
+ TwoPTwoCNIProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView)
+ {
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This is a non-isothermal model, so this method just throws an
+ * exception. This method MUST NOT be overwritten by the actual
+ * problem.
+ */
+ Scalar temperature() const
+ { DUNE_THROW(Dune::Exception, "temperature() method called for a 2p2cni problem"); };
+
+ // \}
+};
+
+}
+
+#endif
diff --git a/dumux/boxmodels/new_2p2cni/2p2cniproperties.hh b/dumux/boxmodels/new_2p2cni/2p2cniproperties.hh
new file mode 100644
index 0000000000..3bff17423b
--- /dev/null
+++ b/dumux/boxmodels/new_2p2cni/2p2cniproperties.hh
@@ -0,0 +1,54 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008 by Klaus Mosthaf, Andreas Lauser, Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup TwoPTwoCNIModel
+ */
+/*!
+ * \file
+ *
+ * \brief Defines the properties required for the non-isothermal two-phase,
+ * two-component BOX model.
+ */
+#ifndef DUMUX_2P2CNI_PROPERTIES_HH
+#define DUMUX_2P2CNI_PROPERTIES_HH
+
+#include <dumux/boxmodels/new_2p2c/2p2cproperties.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the non-isothermal two-phase, two-component problems
+NEW_TYPE_TAG(BoxTwoPTwoCNI, INHERITS_FROM(BoxTwoPTwoC));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+NEW_PROP_TAG(TwoPTwoCNIIndices); //!< Enumerations for the 2p2cni models
+}
+}
+
+#endif
diff --git a/dumux/boxmodels/new_2p2cni/2p2cnipropertydefaults.hh b/dumux/boxmodels/new_2p2cni/2p2cnipropertydefaults.hh
new file mode 100644
index 0000000000..9b02d09a3d
--- /dev/null
+++ b/dumux/boxmodels/new_2p2cni/2p2cnipropertydefaults.hh
@@ -0,0 +1,85 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Melanie Darcis *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup TwoPTwoCNIModel
+ */
+/*!
+ * \file
+ *
+ * \brief Defines default values for most properties required by the 2p2cni
+ * box model.
+ */
+#ifndef DUMUX_2P2CNI_PROPERTY_DEFAULTS_HH
+#define DUMUX_2P2CNI_PROPERTY_DEFAULTS_HH
+
+#include <dumux/boxmodels/new_2p2c/2p2cpropertydefaults.hh>
+
+#include "2p2cnimodel.hh"
+#include "2p2cniproblem.hh"
+#include "2p2cniindices.hh"
+#include "2p2cnilocalresidual.hh"
+#include "2p2cnivolumevariables.hh"
+#include "2p2cnifluxvariables.hh"
+
+namespace Dumux
+{
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Property values
+//////////////////////////////////////////////////////////////////
+
+SET_INT_PROP(BoxTwoPTwoCNI, NumEq, 3); //!< set the number of equations to 3
+
+//! Use the 2p2cni local jacobian operator for the 2p2cni model
+SET_TYPE_PROP(BoxTwoPTwoCNI,
+ LocalResidual,
+ TwoPTwoCNILocalResidual<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(BoxTwoPTwoCNI, Model, TwoPTwoCNIModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(BoxTwoPTwoCNI, VolumeVariables, TwoPTwoCNIVolumeVariables<TypeTag>);
+
+
+//! the FluxVariables property
+SET_TYPE_PROP(BoxTwoPTwoCNI, FluxVariables, TwoPTwoCNIFluxVariables<TypeTag>);
+
+//! The indices required by the non-isothermal 2p2c model
+SET_PROP(BoxTwoPTwoCNI, TwoPTwoCIndices)
+{ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCNIIndices)) type; };
+
+SET_PROP(BoxTwoPTwoCNI, TwoPTwoCNIIndices)
+{ private:
+ enum { formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation)) };
+public:
+ typedef TwoPTwoCNIIndices<TypeTag, formulation, 0> type;
+};
+
+}
+
+}
+#endif
diff --git a/dumux/boxmodels/new_2p2cni/2p2cnivolumevariables.hh b/dumux/boxmodels/new_2p2cni/2p2cnivolumevariables.hh
new file mode 100644
index 0000000000..b5ff5240ec
--- /dev/null
+++ b/dumux/boxmodels/new_2p2cni/2p2cnivolumevariables.hh
@@ -0,0 +1,129 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Melanie Darcis *
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Contains the quantities which are constant within a
+ * finite volume in the non-isothermal two-phase, two-component
+ * model.
+ */
+#ifndef DUMUX_2P2CNI_VOLUME_VARIABLES_HH
+#define DUMUX_2P2CNI_VOLUME_VARIABLES_HH
+
+#include <dumux/boxmodels/new_2p2c/2p2cvolumevariables.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPTwoCNIModel
+ * \brief Contains the quantities which are are constant within a
+ * finite volume in the non-isothermal two-phase, two-component
+ * model.
+ */
+template <class TypeTag>
+class TwoPTwoCNIVolumeVariables : public TwoPTwoCVolumeVariables<TypeTag>
+{
+ //! \cond 0
+ typedef TwoPTwoCVolumeVariables<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+ enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)) };
+ enum { numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)) };
+ enum { temperatureIdx = Indices::temperatureIdx };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
+ //! \endcond
+
+public:
+ /*!
+ * \brief Update the temperature of the sub-control volume.
+ */
+ Scalar getTemperature(const PrimaryVariables &sol,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ const Problem &problem) const
+ {
+ // retrieve temperature from solution vector
+ return sol[temperatureIdx];
+ }
+
+ /*!
+ * \brief Update all quantities for a given control volume.
+ *
+ * \param sol The solution primary variables
+ * \param problem The problem
+ * \param element The element
+ * \param elemGeom Evaluate function with solution of current or previous time step
+ * \param vertIdx The local index of the SCV (sub-control volume)
+ * \param isOldSol Evaluate function with solution of current or previous time step
+ */
+ template <class MutableParams>
+ void updateEnergy(MutableParams &mutParams,
+ const PrimaryVariables &sol,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ const Problem &problem)
+ {
+ heatCapacity_ =
+ problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
+ Valgrind::CheckDefined(heatCapacity_);
+
+ // the internal energies and the enthalpies
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ Scalar h =
+ FluidSystem::computeEnthalpy(mutParams, phaseIdx);
+ mutParams.setEnthalpy(phaseIdx, h);
+ }
+ }
+
+ /*!
+ * \brief Returns the total heat capacity \f$\mathrm{[J/(K*m^3]}\f$ of the rock matrix in
+ * the sub-control volume.
+ */
+ Scalar heatCapacity() const
+ { return heatCapacity_; };
+
+protected:
+ Scalar heatCapacity_;
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/boxmodels/new_2p2cni/Makefile.am b/dumux/boxmodels/new_2p2cni/Makefile.am
new file mode 100644
index 0000000000..79135402f3
--- /dev/null
+++ b/dumux/boxmodels/new_2p2cni/Makefile.am
@@ -0,0 +1,4 @@
+2p2cnidir = $(includedir)/dumux/boxmodels/2p2cni
+2p2cni_HEADERS = *.hh
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/boxmodels/2p2cni/waterairproblem.hh b/test/boxmodels/2p2cni/waterairproblem.hh
index f329e1888e..8ccc478661 100644
--- a/test/boxmodels/2p2cni/waterairproblem.hh
+++ b/test/boxmodels/2p2cni/waterairproblem.hh
@@ -33,7 +33,6 @@
#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
#include <dumux/material/fluidsystems/h2o_n2_system.hh>
-#include <dumux/material/new_fluidsystems/h2on2fluidsystem.hh>
#include <dumux/boxmodels/2p2cni/2p2cnimodel.hh>
@@ -75,12 +74,8 @@ SET_PROP(WaterAirProblem, Problem)
typedef Dumux::WaterAirProblem<TypeTag> type;
};
-// Set fluid configuration
-SET_PROP(WaterAirProblem, FluidSystem)
-{
-private: typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
-public: typedef Dumux::H2ON2FluidSystem<Scalar> type;
-};
+// Set the wetting phase
+SET_TYPE_PROP(WaterAirProblem, FluidSystem, Dumux::H2O_N2_System<TypeTag>);
// Set the spatial parameters
SET_TYPE_PROP(WaterAirProblem,
@@ -131,7 +126,7 @@ SET_BOOL_PROP(WaterAirProblem, NewtonWriteConvergence, false);
* To run the simulation execute the following line in shell:
* <tt>./test_2p2cni ./grids/test_2p2cni.dgf 300000 1000</tt>
* */
-template <class TypeTag>
+template <class TypeTag >
class WaterAirProblem : public TwoPTwoCNIProblem<TypeTag>
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
diff --git a/test/boxmodels/Makefile.am b/test/boxmodels/Makefile.am
index d64fbf4578..29b296e1ec 100644
--- a/test/boxmodels/Makefile.am
+++ b/test/boxmodels/Makefile.am
@@ -1,5 +1,5 @@
-SUBDIRS = 1p 1p2c 2p 2p2c 2p2cni 2pni richards
-
+SUBDIRS = 1p 1p2c 2p 2p2c 2p2cni new_2p2c new_2p2cni 2pni richards
+
EXTRA_DIST = CMakeLists.txt
include $(top_srcdir)/am/global-rules
diff --git a/test/boxmodels/new_2p2c/CMakeLists.txt b/test/boxmodels/new_2p2c/CMakeLists.txt
new file mode 100644
index 0000000000..9faa5f7802
--- /dev/null
+++ b/test/boxmodels/new_2p2c/CMakeLists.txt
@@ -0,0 +1,16 @@
+# build target for the twophase-twocomponent test problem
+ADD_EXECUTABLE("test_2p2c" test_2p2c.cc)
+TARGET_LINK_LIBRARIES("test_2p2c" ${DumuxLinkLibraries})
+
+# add required libraries and includes to the build flags
+LINK_DIRECTORIES(${DumuxLinkDirectories})
+INCLUDE_DIRECTORIES(${DumuxIncludeDirectories})
+
+# make sure the grids are present in the build directory
+add_custom_command(TARGET "test_2p2c"
+ POST_BUILD
+ COMMAND ${CMAKE_COMMAND} -E
+ copy_directory
+ "${CMAKE_CURRENT_SOURCE_DIR}/grids"
+ "${CMAKE_CURRENT_BINARY_DIR}/grids")
+
diff --git a/test/boxmodels/new_2p2c/Makefile.am b/test/boxmodels/new_2p2c/Makefile.am
new file mode 100644
index 0000000000..9d441002d4
--- /dev/null
+++ b/test/boxmodels/new_2p2c/Makefile.am
@@ -0,0 +1,8 @@
+check_PROGRAMS = test_2p2c
+
+noinst_HEADERS = *.hh
+EXTRA_DIST = grids/*.dgf CMakeLists.txt
+
+test_2p2c_SOURCES = test_2p2c.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/boxmodels/new_2p2c/grids/test_2p2c.dgf b/test/boxmodels/new_2p2c/grids/test_2p2c.dgf
new file mode 100644
index 0000000000..765233c04a
--- /dev/null
+++ b/test/boxmodels/new_2p2c/grids/test_2p2c.dgf
@@ -0,0 +1,10 @@
+DGF
+Interval
+0 0 % first corner
+60 40 % second corner
+24 16 % 24 cells in x and 16 in y direction
+#
+BOUNDARYDOMAIN
+default 1 % all boundaries have id 1
+#BOUNDARYDOMAIN
+# unitcube.dgf
diff --git a/test/boxmodels/new_2p2c/injectionproblem.hh b/test/boxmodels/new_2p2c/injectionproblem.hh
new file mode 100644
index 0000000000..b77da478e0
--- /dev/null
+++ b/test/boxmodels/new_2p2c/injectionproblem.hh
@@ -0,0 +1,379 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Definition of a problem, where air is injected under a low permeable layer.
+ */
+
+#ifndef DUMUX_INJECTIONPROBLEM_HH
+#define DUMUX_INJECTIONPROBLEM_HH
+
+#include <dune/grid/io/file/dgfparser/dgfug.hh>
+#include <dune/grid/io/file/dgfparser/dgfs.hh>
+#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
+
+#include <dumux/boxmodels/new_2p2c/2p2cmodel.hh>
+#include <dumux/material/new_fluidsystems/h2on2fluidsystem.hh>
+
+#include "injectionspatialparameters.hh"
+
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class InjectionProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(InjectionProblem, INHERITS_FROM(BoxTwoPTwoC));
+
+// Set the grid type
+SET_PROP(InjectionProblem, Grid)
+{
+ typedef Dune::SGrid<2,2> type;
+};
+
+
+#if HAVE_DUNE_PDELAB
+SET_PROP(InjectionProblem, LocalFEMSpace)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum{dim = GridView::dimension};
+
+public:
+ typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for cubes
+// typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices
+};
+#endif // HAVE_DUNE_PDELAB
+
+// Set the problem property
+SET_PROP(InjectionProblem, Problem)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+ typedef Dumux::InjectionProblem<TTAG(InjectionProblem)> type;
+};
+
+// Set fluid configuration
+SET_PROP(InjectionProblem, FluidSystem)
+{
+private: typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public: typedef Dumux::H2ON2FluidSystem<Scalar> type;
+};
+
+// Set the spatial parameters
+SET_TYPE_PROP(InjectionProblem,
+ SpatialParameters,
+ Dumux::InjectionSpatialParameters<TypeTag>);
+
+// Enable gravity
+SET_BOOL_PROP(InjectionProblem, EnableGravity, true);
+
+// Enable gravity
+SET_INT_PROP(InjectionProblem, NewtonLinearSolverVerbosity, 0);
+
+SET_BOOL_PROP(InjectionProblem, EnableJacobianRecycling, false);
+SET_BOOL_PROP(InjectionProblem, EnablePartialReassemble, false);
+}
+
+
+/*!
+ * \ingroup TwoPTwoCModel
+ * \brief Problem where air is injected under a low permeable layer in a depth of 800m.
+ *
+ * The domain is sized 60m times 40m and consists of two layers, a moderately
+ * permeable spatial parameters (\f$ K=10e-12\f$) for \f$ y>22m\f$ and one with a lower permeablility (\f$ K=10e-13\f$)
+ * in the rest of the domain.
+ *
+ * Air enters a water-filled aquifer, which is situated 800m below sea level, at the right boundary
+ * (\f$ 5m<y<15m\f$) and migrates upwards due to buoyancy. It accumulates and
+ * partially enters the lower permeable aquitard.
+ * This problem uses the \ref TwoPTwoCModel.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./test_2p2c grids/test_2p2c.dgf 1e6 1e4 </tt>
+ */
+template <class TypeTag = TTAG(InjectionProblem) >
+class InjectionProblem : public TwoPTwoCProblem<TypeTag>
+{
+ typedef InjectionProblem<TypeTag> ThisType;
+ typedef TwoPTwoCProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+ enum {
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+ };
+
+ // copy some indices for convenience
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+ enum {
+ lPhaseIdx = Indices::lPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+
+ lCompIdx = Indices::lCompIdx,
+ gCompIdx = Indices::gCompIdx,
+
+ contiLEqIdx = Indices::contiLEqIdx,
+ contiGEqIdx = Indices::contiGEqIdx,
+ };
+
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param timeManager The time manager
+ * \param gridView The grid view
+ */
+ InjectionProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView)
+ {
+ Scalar Tmin = 273.0 + 35;
+ Scalar Tmax = 273.0 + 45;
+ Scalar nT = 5;
+
+ Scalar pmin = 25e6;
+ Scalar pmax = 30e6;
+ Scalar np = 250;
+ FluidSystem::init(Tmin, Tmax, nT,
+ pmin, pmax, np);
+ }
+
+ /*!
+ * \brief Called directly after the time integration.
+ */
+ void postTimeStep()
+ {
+ // Calculate storage terms
+ PrimaryVariables storageL, storageG;
+ this->model().globalPhaseStorage(storageL, lPhaseIdx);
+ this->model().globalPhaseStorage(storageG, gPhaseIdx);
+
+ // Write mass balance information for rank 0
+ if (this->gridView().comm().rank() == 0) {
+ std::cout<<"Storage: liquid=[" << storageL << "]"
+ << " gas=[" << storageG << "]\n";
+ }
+ }
+
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ { return "injection"; }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * \param element The element
+ * \param fvElemGeom The finite-volume geometry in the box scheme
+ * \param scvIdx The local vertex index (SCV index)
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return temperature_;
+ };
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ *
+ * \param values The boundary types for the conservation equations
+ * \param vertex The vertex for which the boundary type is set
+ */
+ void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
+
+ if (globalPos[0] < eps_)
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * boundary segment.
+ *
+ * \param values The dirichlet values for the primary variables
+ * \param vertex The vertex for which the boundary type is set
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
+
+ initial_(values, globalPos);
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * \param values The neumann values for the conservation equations
+ * \param element The finite element
+ * \param fvElemGeom The finite-volume geometry in the box scheme
+ * \param is The intersection between element and boundary
+ * \param scvIdx The local vertex index
+ * \param boundaryFaceIdx The index of the boundary face
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each phase. Negative values mean influx.
+ */
+ void neumann(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ values = 0;
+ if (globalPos[1] < 15 && globalPos[1] > 5) {
+ values[contiGEqIdx] = -1e-3; // kg/(s*m^2)
+ }
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * \param values The source and sink values for the conservation equations
+ * \param element The finite element
+ * \param fvElemGeom The finite-volume geometry in the box scheme
+ * \param scvIdx The local vertex index
+ *
+ * For this method, the \a values parameter stores the rate mass
+ * generated or annihilate per volume unit. Positive values mean
+ * that mass is created, negative ones mean that it vanishes.
+ */
+ void source(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ values = Scalar(0.0);
+ }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param values The initial values for the primary variables
+ * \param element The finite element
+ * \param fvElemGeom The finite-volume geometry in the box scheme
+ * \param scvIdx The local vertex index
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ void initial(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+ initial_(values, globalPos);
+ }
+
+ /*!
+ * \brief Return the initial phase state inside a control volume.
+ *
+ * \param vert The vertex
+ * \param globalIdx The index of the global vertex
+ * \param globalPos The global position
+ */
+ int initialPhasePresence(const Vertex &vert,
+ int &globalIdx,
+ const GlobalPosition &globalPos) const
+ { return Indices::lPhaseOnly; }
+
+ // \}
+
+private:
+ // the internal method for the initial condition
+ void initial_(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ Scalar densityW = FluidSystem::H2O::liquidDensity(temperature_, 1e5);
+
+ values[Indices::plIdx] = 1e5 - densityW*this->gravity()[1]*(depthBOR_ - globalPos[1]);
+ values[Indices::SgOrXIdx] =
+ values[Indices::plIdx]*0.95/
+ BinaryCoeff::H2O_N2::henry(temperature_);
+ }
+
+ static const Scalar temperature_ = 273.15 + 40; // [K]
+ static const Scalar depthBOR_ = 2700.0; // [m]
+ static const Scalar eps_ = 1e-6;
+};
+} //end namespace
+
+#endif
diff --git a/test/boxmodels/new_2p2c/injectionspatialparameters.hh b/test/boxmodels/new_2p2c/injectionspatialparameters.hh
new file mode 100644
index 0000000000..54715db3d6
--- /dev/null
+++ b/test/boxmodels/new_2p2c/injectionspatialparameters.hh
@@ -0,0 +1,268 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Definition of the spatial parameters for the injection
+ * problem which uses the isothermal 2p2c box model
+ */
+
+#ifndef DUMUX_INJECTION_SPATIAL_PARAMETERS_HH
+#define DUMUX_INJECTION_SPATIAL_PARAMETERS_HH
+
+#include <dumux/material/spatialparameters/boxspatialparameters.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPTwoCModel
+ *
+ * \brief Definition of the spatial parameters for the injection
+ * problem which uses the isothermal 2p2c box model
+ */
+template<class TypeTag>
+class InjectionSpatialParameters : public BoxSpatialParameters<TypeTag>
+{
+ typedef BoxSpatialParameters<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename Grid::ctype CoordScalar;
+ enum {
+ dim=GridView::dimension,
+ dimWorld=GridView::dimensionworld,
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ enum {
+ lPhaseIdx = FluidSystem::lPhaseIdx,
+ gPhaseIdx = FluidSystem::gPhaseIdx,
+ };
+
+ typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
+ typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<CoordScalar,dimWorld> Vector;
+
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ typedef RegularizedBrooksCorey<Scalar> EffMaterialLaw;
+ //typedef LinearMaterial<Scalar> EffMaterialLaw;
+public:
+ typedef EffToAbsLaw<EffMaterialLaw> MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ /*!
+ * \brief The constructor
+ *
+ * \param gv The grid view
+ */
+ InjectionSpatialParameters(const GridView &gv)
+ : ParentType(gv)
+ {
+ layerBottom_ = 22.0;
+
+ // intrinsic permeabilities
+ fineK_ = 1e-13;
+ coarseK_ = 1e-12;
+
+ // porosities
+ finePorosity_ = 0.3;
+ coarsePorosity_ = 0.3;
+
+ // residual saturations
+ fineMaterialParams_.setSwr(0.2);
+ fineMaterialParams_.setSnr(0.0);
+ coarseMaterialParams_.setSwr(0.2);
+ coarseMaterialParams_.setSnr(0.0);
+
+ // parameters for the Brooks-Corey law
+ fineMaterialParams_.setPe(1e4);
+ coarseMaterialParams_.setPe(1e4);
+ fineMaterialParams_.setAlpha(2.0);
+ coarseMaterialParams_.setAlpha(2.0);
+ }
+
+ ~InjectionSpatialParameters()
+ {}
+
+
+ /*!
+ * \brief Update the spatial parameters with the flow solution
+ * after a timestep.
+ *
+ * \param globalSolution The global solution vector
+ */
+ void update(const SolutionVector &globalSolution)
+ {
+ };
+
+ /*!
+ * \brief Apply the intrinsic permeability tensor to a pressure
+ * potential gradient.
+ *
+ * \param element The current finite element
+ * \param fvElemGeom The current finite volume geometry of the element
+ * \param scvIdx The index of the sub-control volume
+ */
+ const Scalar intrinsicPermeability(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
+ if (isFineMaterial_(pos))
+ return fineK_;
+ return coarseK_;
+ }
+
+ /*!
+ * \brief Define the porosity \f$[-]\f$ of the spatial parameters
+ *
+ * \param element The finite element
+ * \param fvElemGeom The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume where
+ * the porosity needs to be defined
+ */
+ double porosity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
+ if (isFineMaterial_(pos))
+ return finePorosity_;
+ return coarsePorosity_;
+ }
+
+
+ /*!
+ * \brief return the parameter object for the Brooks-Corey material law which depends on the position
+ *
+ * \param element The current finite element
+ * \param fvElemGeom The current finite volume geometry of the element
+ * \param scvIdx The index of the sub-control volume
+ */
+ const MaterialLawParams& materialLawParams(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
+ if (isFineMaterial_(pos))
+ return fineMaterialParams_;
+ return coarseMaterialParams_;
+ }
+
+ /*!
+ * \brief Returns the heat capacity \f$[J/m^3 K]\f$ of the rock matrix.
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param element The finite element
+ * \param fvElemGeom The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume where
+ * the heat capacity needs to be defined
+ */
+ double heatCapacity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return
+ 790 // specific heat capacity of granite [J / (kg K)]
+ * 2700 // density of granite [kg/m^3]
+ * (1 - porosity(element, fvElemGeom, scvIdx));
+ }
+
+ /*!
+ * \brief Calculate the heat flux \f$[W/m^2]\f$ through the
+ * rock matrix based on the temperature gradient \f$[K / m]\f$
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param heatFlux The resulting heat flux vector
+ * \param fluxDat The flux variables
+ * \param vDat The volume variables
+ * \param tempGrad The temperature gradient
+ * \param element The current finite element
+ * \param fvElemGeom The finite volume geometry of the current element
+ * \param scvfIdx The local index of the sub-control volume face where
+ * the matrix heat flux should be calculated
+ */
+ void matrixHeatFlux(Vector &heatFlux,
+ const FluxVariables &fluxDat,
+ const ElementVolumeVariables &vDat,
+ const Vector &tempGrad,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvfIdx) const
+ {
+ static const Scalar lWater = 0.6;
+ static const Scalar lGranite = 2.8;
+
+ // arithmetic mean of the liquid saturation and the porosity
+ const int i = fvElemGeom.subContVolFace[scvfIdx].i;
+ const int j = fvElemGeom.subContVolFace[scvfIdx].j;
+ Scalar Sl = std::max(0.0, (vDat[i].saturation(lPhaseIdx) +
+ vDat[j].saturation(lPhaseIdx)) / 2);
+ Scalar poro = (porosity(element, fvElemGeom, i) +
+ porosity(element, fvElemGeom, j)) / 2;
+
+ Scalar lsat = pow(lGranite, (1-poro)) * pow(lWater, poro);
+ Scalar ldry = pow(lGranite, (1-poro));
+
+ // the heat conductivity of the matrix. in general this is a
+ // tensorial value, but we assume isotropic heat conductivity.
+ Scalar heatCond = ldry + sqrt(Sl) * (ldry - lsat);
+
+ // the matrix heat flux is the negative temperature gradient
+ // times the heat conductivity.
+ heatFlux = tempGrad;
+ heatFlux *= -heatCond;
+ }
+
+private:
+ bool isFineMaterial_(const GlobalPosition &pos) const
+ { return pos[dim-1] > layerBottom_; };
+
+ Scalar fineK_;
+ Scalar coarseK_;
+ Scalar layerBottom_;
+
+ Scalar finePorosity_;
+ Scalar coarsePorosity_;
+
+ MaterialLawParams fineMaterialParams_;
+ MaterialLawParams coarseMaterialParams_;
+};
+
+}
+
+#endif
diff --git a/test/boxmodels/new_2p2c/test_2p2c.cc b/test/boxmodels/new_2p2c/test_2p2c.cc
new file mode 100644
index 0000000000..05d1c1ccdd
--- /dev/null
+++ b/test/boxmodels/new_2p2c/test_2p2c.cc
@@ -0,0 +1,34 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008 by Klaus Mosthaf *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief test for the 2p2c box model
+ */
+#include "config.h"
+#include "injectionproblem.hh"
+#include <dumux/common/start.hh>
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(InjectionProblem) ProblemTypeTag;
+ return Dumux::startFromDGF<ProblemTypeTag>(argc, argv);
+}
diff --git a/test/boxmodels/new_2p2cni/CMakeLists.txt b/test/boxmodels/new_2p2cni/CMakeLists.txt
new file mode 100644
index 0000000000..b88b7db7f3
--- /dev/null
+++ b/test/boxmodels/new_2p2cni/CMakeLists.txt
@@ -0,0 +1,16 @@
+# build target for the 2p2cni test problem
+ADD_EXECUTABLE("test_2p2cni" test_2p2cni.cc)
+TARGET_LINK_LIBRARIES("test_2p2cni" ${DumuxLinkLibraries})
+
+# add required libraries and includes to the build flags
+LINK_DIRECTORIES(${DumuxLinkDirectories})
+INCLUDE_DIRECTORIES(${DumuxIncludeDirectories})
+
+# make sure the grids are present in the build directory
+add_custom_command(TARGET "test_2p2cni"
+ POST_BUILD
+ COMMAND ${CMAKE_COMMAND} -E
+ copy_directory
+ "${CMAKE_CURRENT_SOURCE_DIR}/grids"
+ "${CMAKE_CURRENT_BINARY_DIR}/grids")
+
diff --git a/test/boxmodels/new_2p2cni/Makefile.am b/test/boxmodels/new_2p2cni/Makefile.am
new file mode 100644
index 0000000000..5b3d1ff7e0
--- /dev/null
+++ b/test/boxmodels/new_2p2cni/Makefile.am
@@ -0,0 +1,8 @@
+check_PROGRAMS = test_2p2cni
+
+noinst_HEADERS = *.hh
+EXTRA_DIST = grids/*.dgf CMakeLists.txt
+
+test_2p2cni_SOURCES = test_2p2cni.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/boxmodels/new_2p2cni/grids/test_2p2cni.dgf b/test/boxmodels/new_2p2cni/grids/test_2p2cni.dgf
new file mode 100644
index 0000000000..a3f2629ae7
--- /dev/null
+++ b/test/boxmodels/new_2p2cni/grids/test_2p2cni.dgf
@@ -0,0 +1,15 @@
+DGF
+Interval
+0 0 % first corner
+40 40 % second corner
+64 64 % 24 cells in x and 16 in y direction
+#
+GridParameter
+overlap 0
+periodic
+closure green
+#GridParameter
+BOUNDARYDOMAIN
+default 1 % all boundaries have id 1
+#BOUNDARYDOMAIN
+# unitcube.dgf
diff --git a/test/boxmodels/new_2p2cni/test_2p2cni.cc b/test/boxmodels/new_2p2cni/test_2p2cni.cc
new file mode 100644
index 0000000000..76b5ffb3e8
--- /dev/null
+++ b/test/boxmodels/new_2p2cni/test_2p2cni.cc
@@ -0,0 +1,37 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2007-2008 by Klaus Mosthaf *
+ * Copyright (C) 2007-2008 by Melanie Darcis *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief test for the 2p2cni box model
+ */
+#include "config.h"
+#include "waterairproblem.hh"
+#include <dumux/common/start.hh>
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(WaterAirProblem) ProblemTypeTag;
+ return Dumux::startFromDGF<ProblemTypeTag>(argc, argv);
+}
diff --git a/test/boxmodels/new_2p2cni/waterairproblem.hh b/test/boxmodels/new_2p2cni/waterairproblem.hh
new file mode 100644
index 0000000000..3582d954d8
--- /dev/null
+++ b/test/boxmodels/new_2p2cni/waterairproblem.hh
@@ -0,0 +1,397 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2009 by Klaus Mosthaf *
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Non-isothermal gas injection problem where a gas (e.g. air)
+ * is injected into a fully water saturated medium.
+ */
+#ifndef DUMUX_WATERAIRPROBLEM_HH
+#define DUMUX_WATERAIRPROBLEM_HH
+
+#include <dune/grid/io/file/dgfparser/dgfug.hh>
+#include <dune/grid/io/file/dgfparser/dgfs.hh>
+#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
+
+//#include <dumux/material/fluidsystems/h2o_n2_system.hh>
+#include <dumux/material/new_fluidsystems/h2on2fluidsystem.hh>
+
+#include <dumux/boxmodels/new_2p2cni/2p2cnimodel.hh>
+
+#include "waterairspatialparameters.hh"
+
+#define ISOTHERMAL 0
+
+namespace Dumux
+{
+template <class TypeTag>
+class WaterAirProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(WaterAirProblem, INHERITS_FROM(BoxTwoPTwoCNI));
+
+// Set the grid type
+SET_PROP(WaterAirProblem, Grid)
+{
+ typedef Dune::YaspGrid<2> type;
+};
+
+#if HAVE_DUNE_PDELAB
+SET_PROP(WaterAirProblem, LocalFEMSpace)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum{dim = GridView::dimension};
+
+public:
+ typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for cubes
+// typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices
+};
+#endif // HAVE_DUNE_PDELAB
+
+// Set the problem property
+SET_PROP(WaterAirProblem, Problem)
+{
+ typedef Dumux::WaterAirProblem<TypeTag> type;
+};
+
+// Set fluid configuration
+SET_PROP(WaterAirProblem, FluidSystem)
+{
+private: typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public: typedef Dumux::H2ON2FluidSystem<Scalar> type;
+};
+
+// Set the spatial parameters
+SET_TYPE_PROP(WaterAirProblem,
+ SpatialParameters,
+ Dumux::WaterAirSpatialParameters<TypeTag>);
+
+// Enable gravity
+SET_BOOL_PROP(WaterAirProblem, EnableGravity, true);
+
+// Use forward differences instead of central differences
+SET_INT_PROP(WaterAirProblem, NumericDifferenceMethod, +1);
+
+// Write newton convergence
+SET_BOOL_PROP(WaterAirProblem, NewtonWriteConvergence, false);
+}
+
+
+/*!
+ * \ingroup TwoPTwoCNIBoxModel
+ *
+ * \brief Non-isothermal gas injection problem where a gas (e.g. air)
+ * is injected into a fully water saturated medium. During
+ * buoyancy driven upward migration the gas passes a high
+ * temperature area.
+ *
+ * The domain is sized 40 m times 40 m. The rectangular area with the
+ * increased temperature (380 K) starts at (20 m, 5 m) and ends at (30
+ * m, 35 m).
+ *
+ * For the mass conservation equation neumann boundary conditions are used on
+ * the top and on the bottom of the domain, while dirichlet conditions
+ * apply on the left and the right boundary.
+ * For the energy conservation equation dirichlet boundary conditions are applied
+ * on all boundaries.
+ *
+ * Gas is injected at the bottom boundary from 15 m to 25 m at a rate of
+ * 0.001 kg/(s m), the remaining neumann boundaries are no-flow
+ * boundaries.
+ *
+ * At the dirichlet boundaries a hydrostatic pressure, a gas saturation of zero and
+ * a geothermal temperature gradient of 0.03 K/m are applied.
+ *
+ * This problem uses the \ref TwoPTwoCNIModel.
+ *
+ * This problem should typically be simulated for 300000 s.
+ * A good choice for the initial time step size is 1000 s.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./test_2p2cni ./grids/test_2p2cni.dgf 300000 1000</tt>
+ * */
+template <class TypeTag>
+class WaterAirProblem : public TwoPTwoCNIProblem<TypeTag>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
+ typedef typename GridView::Grid Grid;
+
+ typedef WaterAirProblem<TypeTag> ThisType;
+ typedef TwoPTwoCNIProblem<TypeTag> ParentType;
+
+ // copy some indices for convenience
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+ enum {
+ numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
+
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx,
+#if !ISOTHERMAL
+ temperatureIdx = Indices::temperatureIdx,
+ energyEqIdx = Indices::energyEqIdx,
+#endif
+
+ // Phase State
+ lPhaseOnly = Indices::lPhaseOnly,
+ gPhaseOnly = Indices::gPhaseOnly,
+ bothPhases = Indices::bothPhases,
+
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+ };
+
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param timeManager The time manager
+ * \param gridView The grid view
+ */
+ WaterAirProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView)
+ {
+ Scalar Tmin = 280.0;
+ Scalar Tmax = 390.0;
+ Scalar nT = 50;
+
+ Scalar pmin = 8e6;
+ Scalar pmax = 12e6;
+ Scalar np = 250;
+ FluidSystem::init(Tmin, Tmax, nT,
+ pmin, pmax, np);
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ { return "waterair"; }
+
+#if ISOTHERMAL
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * \param element The element
+ * \param fvElemGeom The finite-volume geometry in the box scheme
+ * \param scvIdx The local vertex index (SCV index)
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return 273.15 + 10; // -> 10°C
+ };
+#endif
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ *
+ * \param values The boundary types for the conservation equations
+ * \param vertex The vertex for which the boundary type is set
+ */
+ void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
+
+ if(globalPos[0] > 40 - eps_ || globalPos[0] < eps_)
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+
+#if !ISOTHERMAL
+ values.setDirichlet(temperatureIdx, energyEqIdx);
+#endif
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * boundary segment.
+ *
+ * \param values The dirichlet values for the primary variables
+ * \param vertex The vertex for which the boundary type is set
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+ {
+ const GlobalPosition globalPos = vertex.geometry().center();
+
+ initial_(values, globalPos);
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * \param values The neumann values for the conservation equations
+ * \param element The finite element
+ * \param fvElemGeom The finite-volume geometry in the box scheme
+ * \param is The intersection between element and boundary
+ * \param scvIdx The local vertex index
+ * \param boundaryFaceIdx The index of the boundary face
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each phase. Negative values mean influx.
+ */
+ void neumann(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+ values = 0;
+
+ // negative values for injection
+ if (globalPos[0] > 15 && globalPos[0] < 25 &&
+ globalPos[1] < eps_)
+ {
+ values[Indices::contiGEqIdx] = -1e-3;
+ }
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * \param values The source and sink values for the conservation equations
+ * \param element The finite element
+ * \param fvElemGeom The finite-volume geometry in the box scheme
+ * \param scvIdx The local vertex index
+ *
+ * For this method, the \a values parameter stores the rate mass
+ * generated or annihilate per volume unit. Positive values mean
+ * that mass is created, negative ones mean that it vanishes.
+ */
+ void source(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ values = Scalar(0.0);
+ }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param values The initial values for the primary variables
+ * \param element The finite element
+ * \param fvElemGeom The finite-volume geometry in the box scheme
+ * \param scvIdx The local vertex index
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ void initial(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ initial_(values, globalPos);
+
+#if !ISOTHERMAL
+ if (globalPos[0] > 20 && globalPos[0] < 30 && globalPos[1] < 30)
+ values[temperatureIdx] = 380;
+#endif
+ }
+
+ /*!
+ * \brief Return the initial phase state inside a control volume.
+ *
+ * \param vert The vertex
+ * \param globalIdx The index of the global vertex
+ * \param globalPos The global position
+ */
+ int initialPhasePresence(const Vertex &vert,
+ int &globalIdx,
+ const GlobalPosition &globalPos) const
+ {
+ return lPhaseOnly;
+ }
+
+private:
+ // internal method for the initial condition (reused for the
+ // dirichlet conditions!)
+ void initial_(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ Scalar densityW = 1000.0;
+ values[pressureIdx] = 1e5 + (depthBOR_ - globalPos[1])*densityW*9.81;
+ values[switchIdx] = 0.0;
+#if !ISOTHERMAL
+ values[temperatureIdx] = 283.0 + (depthBOR_ - globalPos[1])*0.03;
+#endif
+ }
+
+ static const Scalar depthBOR_ = 1000.0; // [m]
+ static const Scalar eps_ = 1e-6;
+};
+} //end namespace
+
+#endif
diff --git a/test/boxmodels/new_2p2cni/waterairspatialparameters.hh b/test/boxmodels/new_2p2cni/waterairspatialparameters.hh
new file mode 100644
index 0000000000..b525b1e639
--- /dev/null
+++ b/test/boxmodels/new_2p2cni/waterairspatialparameters.hh
@@ -0,0 +1,269 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Definition of the spatial parameters for the water-air problem.
+ */
+#ifndef DUMUX_WATER_AIR_SPATIAL_PARAMETERS_HH
+#define DUMUX_WATER_AIR_SPATIAL_PARAMETERS_HH
+
+#include <dumux/material/spatialparameters/boxspatialparameters.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPTwoCNIModel
+ *
+ * \brief Definition of the spatial parameters for the water-air problem
+ */
+template<class TypeTag>
+class WaterAirSpatialParameters : public BoxSpatialParameters<TypeTag>
+{
+ typedef BoxSpatialParameters<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename Grid::ctype CoordScalar;
+ enum {
+ dim=GridView::dimension,
+ dimWorld=GridView::dimensionworld,
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+ enum {
+ lPhaseIdx = Indices::lPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+ };
+
+ typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
+ typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<CoordScalar,dimWorld> Vector;
+
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+ typedef RegularizedBrooksCorey<Scalar> EffMaterialLaw;
+public:
+ typedef EffToAbsLaw<EffMaterialLaw> MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ /*!
+ * \brief The constructor
+ *
+ * \param gv The grid view
+ */
+ WaterAirSpatialParameters(const GridView &gv)
+ : ParentType(gv)
+ {
+ layerBottom_ = 22.0;
+
+ // intrinsic permeabilities
+ fineK_ = 1e-13;
+ coarseK_ = 1e-12;
+
+ // porosities
+ finePorosity_ = 0.3;
+ coarsePorosity_ = 0.3;
+
+ // residual saturations
+ fineMaterialParams_.setSwr(0.2);
+ fineMaterialParams_.setSnr(0.0);
+ coarseMaterialParams_.setSwr(0.2);
+ coarseMaterialParams_.setSnr(0.0);
+
+ // parameters for the Brooks-Corey law
+ fineMaterialParams_.setPe(1e4);
+ coarseMaterialParams_.setPe(1e4);
+ fineMaterialParams_.setAlpha(2.0);
+ coarseMaterialParams_.setAlpha(2.0);
+ }
+
+ ~WaterAirSpatialParameters()
+ {}
+
+
+ /*!
+ * \brief Update the spatial parameters with the flow solution
+ * after a timestep.
+ *
+ * \param globalSolution The global solution vector
+ */
+ void update(const SolutionVector &globalSolution)
+ {
+ };
+
+ /*!
+ * \brief Apply the intrinsic permeability tensor to a pressure
+ * potential gradient.
+ *
+ * \param element The current finite element
+ * \param fvElemGeom The current finite volume geometry of the element
+ * \param scvIdx The index of the sub-control volume
+ */
+ const Scalar intrinsicPermeability(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
+ if (isFineMaterial_(pos))
+ return fineK_;
+ return coarseK_;
+ }
+
+ /*!
+ * \brief Define the porosity \f$[-]\f$ of the spatial parameters
+ *
+ * \param element The finite element
+ * \param fvElemGeom The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume where
+ * the porosity needs to be defined
+ */
+ double porosity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
+ if (isFineMaterial_(pos))
+ return finePorosity_;
+ else
+ return coarsePorosity_;
+ }
+
+
+ /*!
+ * \brief return the parameter object for the Brooks-Corey material law which depends on the position
+ *
+ * \param element The current finite element
+ * \param fvElemGeom The current finite volume geometry of the element
+ * \param scvIdx The index of the sub-control volume
+ */
+ const MaterialLawParams& materialLawParams(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
+ if (isFineMaterial_(pos))
+ return fineMaterialParams_;
+ else
+ return coarseMaterialParams_;
+ }
+
+ /*!
+ * \brief Returns the heat capacity \f$[J/m^3 K]\f$ of the rock matrix.
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param element The finite element
+ * \param fvElemGeom The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume where
+ * the heat capacity needs to be defined
+ */
+ double heatCapacity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return
+ 790 // specific heat capacity of granite [J / (kg K)]
+ * 2700 // density of granite [kg/m^3]
+ * (1 - porosity(element, fvElemGeom, scvIdx));
+ }
+
+ /*!
+ * \brief Calculate the heat flux \f$[W/m^2]\f$ through the
+ * rock matrix based on the temperature gradient \f$[K / m]\f$
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param heatFlux The resulting heat flux vector
+ * \param fluxDat The flux variables
+ * \param vDat The volume variables
+ * \param tempGrad The temperature gradient
+ * \param element The current finite element
+ * \param fvElemGeom The finite volume geometry of the current element
+ * \param scvfIdx The local index of the sub-control volume face where
+ * the matrix heat flux should be calculated
+ */
+ void matrixHeatFlux(Vector &heatFlux,
+ const FluxVariables &fluxDat,
+ const ElementVolumeVariables &vDat,
+ const Vector &tempGrad,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvfIdx) const
+ {
+ static const Scalar lWater = 0.6;
+ static const Scalar lGranite = 2.8;
+
+ // arithmetic mean of the liquid saturation and the porosity
+ const int i = fvElemGeom.subContVolFace[scvfIdx].i;
+ const int j = fvElemGeom.subContVolFace[scvfIdx].j;
+ Scalar Sl = std::max(0.0, (vDat[i].saturation(lPhaseIdx) +
+ vDat[j].saturation(lPhaseIdx)) / 2);
+ Scalar poro = (porosity(element, fvElemGeom, i) +
+ porosity(element, fvElemGeom, j)) / 2;
+
+ Scalar lsat = pow(lGranite, (1-poro)) * pow(lWater, poro);
+ Scalar ldry = pow(lGranite, (1-poro));
+
+ // the heat conductivity of the matrix. in general this is a
+ // tensorial value, but we assume isotropic heat conductivity.
+ Scalar heatCond = ldry + sqrt(Sl) * (ldry - lsat);
+
+ // the matrix heat flux is the negative temperature gradient
+ // times the heat conductivity.
+ heatFlux = tempGrad;
+ heatFlux *= -heatCond;
+ }
+
+private:
+ bool isFineMaterial_(const GlobalPosition &pos) const
+ { return pos[dim-1] > layerBottom_; };
+
+ Scalar fineK_;
+ Scalar coarseK_;
+ Scalar layerBottom_;
+
+ Scalar finePorosity_;
+ Scalar coarsePorosity_;
+
+ MaterialLawParams fineMaterialParams_;
+ MaterialLawParams coarseMaterialParams_;
+};
+
+}
+
+#endif
--
GitLab