diff --git a/configure.ac b/configure.ac
index 8929970f41bd685cc591108a8fae11ed79098b8a..f92dc7366a16b4c0eccd3e567787b09c1f7738d9 100644
--- a/configure.ac
+++ b/configure.ac
@@ -113,7 +113,7 @@ AC_CONFIG_FILES([dumux.pc
dumux/decoupled/2p/diffusion/mimetic/Makefile
dumux/decoupled/2p/impes/Makefile
dumux/decoupled/2p/transport/Makefile
- dumux/decoupled/2p/transport/fv/Makefile
+ dumux/decoupled/2p/transport/fv/Makefile
dumux/decoupled/common/Makefile
dumux/io/Makefile
dumux/material/Makefile
@@ -140,6 +140,7 @@ AC_CONFIG_FILES([dumux.pc
test/decoupled/Makefile
test/decoupled/1p/Makefile
test/decoupled/2p/Makefile
+ test/decoupled/2p2c/Makefile
tutorial/Makefile
])
diff --git a/dumux/decoupled/2p2c/2p2cproblem.hh b/dumux/decoupled/2p2c/2p2cproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0a75294ee741cf911b558854b895df5c3a7e8ecd
--- /dev/null
+++ b/dumux/decoupled/2p2c/2p2cproblem.hh
@@ -0,0 +1,147 @@
+// $Id:
+/*****************************************************************************
+ * Copyright (C) 2010 by Benjamin Faigle, Markus Wolff *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for sequential 2p2c compositional problems
+ */
+#ifndef DUMUX_IMPETPROBLEM_2P2C_HH
+#define DUMUX_IMPETPROBLEM_2P2C_HH
+
+#include "boundaryconditions2p2c.hh"
+#include <dumux/decoupled/common/impet.hh>
+#include <dumux/decoupled/common/impetproblem.hh>
+#include <dumux/decoupled/2p2c/variableclass2p2c.hh>
+#include <dumux/decoupled/2p2c/2p2cproperties.hh>
+
+
+namespace Dumux
+{
+/*!
+ * \ingroup IMPEC
+ * \defgroup IMPECproblems IMPEC problems
+ */
+/*!
+ * \ingroup IMPECproblems
+ * \brief Base class for all compositional 2-phase problems which use an impet algorithm
+ *
+ * Differs from .../2p/impes/impesproblem2p.hh only in the includes:
+ * Usage of the compositional properties and variableclass.
+ */
+template<class TypeTag, class Implementation>
+class IMPETProblem2P2C : public IMPETProblem<TypeTag, Implementation>
+{
+ typedef IMPETProblem<TypeTag, Implementation> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GridView::Grid Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+ // material properties
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ 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 gridView The grid view
+ * \param verbose Output flag for the time manager.
+ */
+ IMPETProblem2P2C(const GridView &gridView, bool verbose = true)
+ : ParentType(gridView, verbose),
+ gravity_(0)
+ {
+ newSpatialParams_ = true;
+ spatialParameters_ = new SpatialParameters(gridView);
+
+ gravity_ = 0;
+ if (GET_PROP_VALUE(TypeTag, PTAG(EnableGravity)))
+ gravity_[dim - 1] = - 9.81;
+ }
+ /*!
+ * \brief The constructor
+ *
+ * \param gridView The grid view
+ * \param SpatialParameters SpatialParameters instantiation
+ * \param verbose Output flag for the time manager.
+ */
+ IMPETProblem2P2C(const GridView &gridView, SpatialParameters &spatialParameters, bool verbose = true)
+ : ParentType(gridView, verbose),
+ gravity_(0),spatialParameters_(&spatialParameters)
+ {
+ newSpatialParams_ = false;
+ gravity_ = 0;
+ if (GET_PROP_VALUE(TypeTag, PTAG(EnableGravity)))
+ gravity_[dim - 1] = - 9.81;
+ }
+
+ virtual ~IMPETProblem2P2C()
+ {
+ if (newSpatialParams_)
+ {
+ delete spatialParameters_;
+ }
+ }
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \copydoc Dumux::IMPESProblem2P::temperature()
+ */
+ Scalar temperature() const
+ { return this->asImp_()->temperature(); };
+
+ /*!
+ * \copydoc Dumux::IMPESProblem2P::gravity()
+ */
+ const GlobalPosition &gravity() const
+ { return gravity_; }
+
+ /*!
+ * \copydoc Dumux::IMPESProblem2P::spatialParameters()
+ */
+ SpatialParameters &spatialParameters()
+ { return *spatialParameters_; }
+
+ /*!
+ * \copydoc Dumux::IMPESProblem2P::spatialParameters()
+ */
+ const SpatialParameters &spatialParameters() const
+ { return *spatialParameters_; }
+
+ // \}
+
+private:
+ GlobalPosition gravity_;
+
+ // fluids and material properties
+ SpatialParameters* spatialParameters_;
+ bool newSpatialParams_;
+};
+
+}
+
+#endif
diff --git a/dumux/decoupled/2p2c/2p2cproperties.hh b/dumux/decoupled/2p2c/2p2cproperties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b1b5b42e7bd03ae69c348d767cfb13eaa76391d0
--- /dev/null
+++ b/dumux/decoupled/2p2c/2p2cproperties.hh
@@ -0,0 +1,264 @@
+// $Id:
+/*****************************************************************************
+ * Copyright (C) 20010 by Benjamin Faigle *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+
+/*!
+ * \ingroup IMPET
+ * \defgroup IMPEC Miscible IMPEC
+ */
+/*!
+ * \ingroup IMPEC
+ * \defgroup multiphase Multiphase Compositional Models
+ */
+/*!
+ * \ingroup IMPEC
+ * \defgroup multiphysics Multiphysics Compositional Models
+ */
+/*!
+ * \ingroup IMPEC
+ * \file
+ *
+ * \brief Defines the properties required for the decoupled 2p2c models.
+ */
+
+#ifndef DUMUX_2P2CPROPERTIES_HH
+#define DUMUX_2P2CPROPERTIES_HH
+
+//Dune-includes
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
+
+//DUMUX includes
+#include <dumux/decoupled/common/impetproperties.hh>
+
+namespace Dumux
+{
+
+////////////////////////////////
+// forward declarations
+////////////////////////////////
+
+
+template<class TypeTag>
+class VariableClass2P2C;
+
+template<class TypeTag>
+class DecTwoPTwoCFluidState;
+
+template <class TypeTag>
+struct TwoPCommonIndices;
+
+////////////////////////////////
+// properties
+////////////////////////////////
+namespace Properties
+{
+
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the two-phase problems
+NEW_TYPE_TAG(DecoupledTwoPTwoC, INHERITS_FROM(IMPET));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG ( TwoPIndices );
+NEW_PROP_TAG( SpatialParameters ); //!< The type of the soil properties object
+NEW_PROP_TAG( EnableGravity); //!< Returns whether gravity is considered in the problem
+NEW_PROP_TAG( PressureFormulation); //!< The formulation of the model
+NEW_PROP_TAG( SaturationFormulation); //!< The formulation of the model
+NEW_PROP_TAG( VelocityFormulation); //!< The formulation of the model
+NEW_PROP_TAG( EnableCompressibility);// ! Returns whether compressibility is allowed
+NEW_PROP_TAG(EnableCapillarity); //!< Returns whether capillarity is regarded
+NEW_PROP_TAG( BoundaryMobility );
+NEW_PROP_TAG( FluidSystem );
+NEW_PROP_TAG( FluidState );
+
+//Properties for linear solvers
+NEW_PROP_TAG(PressureCoefficientMatrix);
+NEW_PROP_TAG(PressureRHSVector);
+NEW_PROP_TAG( PressurePreconditioner );
+NEW_PROP_TAG( PressureSolver );
+NEW_PROP_TAG( SolverParameters );
+NEW_PROP_TAG(ReductionSolver);
+NEW_PROP_TAG(MaxIterationNumberSolver);
+NEW_PROP_TAG(IterationNumberPreconditioner);
+NEW_PROP_TAG(VerboseLevelSolver);
+NEW_PROP_TAG(RelaxationPreconditioner);
+
+//////////////////////////////////////////////////////////////////
+// Properties
+//////////////////////////////////////////////////////////////////
+
+SET_PROP(DecoupledTwoPTwoC, TwoPIndices)
+{
+ typedef TwoPCommonIndices<TypeTag> type;
+};
+
+// set fluid/component information
+SET_PROP(DecoupledTwoPTwoC, NumPhases) //!< The number of phases in the 2p model is 2
+{
+ // the property is created in decoupledproperties.hh
+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 2p2c model!");
+};
+
+SET_PROP(DecoupledTwoPTwoC, NumComponents) //!< The number of components in the 2p2c model is 2
+{
+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 2p2c model!");
+};
+
+//! Set the default formulation
+SET_INT_PROP(DecoupledTwoPTwoC,
+ PressureFormulation,
+ TwoPCommonIndices<TypeTag>::pressureW);
+
+SET_INT_PROP(DecoupledTwoPTwoC,
+ SaturationFormulation,
+ TwoPCommonIndices<TypeTag>::saturationW);
+
+SET_INT_PROP(DecoupledTwoPTwoC,
+ VelocityFormulation,
+ TwoPCommonIndices<TypeTag>::velocityTotal);
+
+SET_PROP(DecoupledTwoPTwoC, TransportSolutionType)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename Dune::BlockVector<Dune::BlockVector<Dune::FieldVector<Scalar,1> > > type;//!<type for vector of vector (of scalars)
+
+};
+
+SET_BOOL_PROP(DecoupledTwoPTwoC, EnableCompressibility, true);
+
+SET_BOOL_PROP(DecoupledTwoPTwoC, EnableCapillarity, false);
+
+
+SET_PROP_DEFAULT(BoundaryMobility)
+{
+ static const int value = TwoPCommonIndices<TypeTag>::satDependent;
+};
+
+
+SET_TYPE_PROP(DecoupledTwoPTwoC, Variables, VariableClass2P2C<TypeTag>);
+
+SET_TYPE_PROP(DecoupledTwoPTwoC, FluidState, DecTwoPTwoCFluidState<TypeTag>);
+
+// solver stuff
+SET_PROP(DecoupledTwoPTwoC, PressureCoefficientMatrix)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef Dune::FieldMatrix<Scalar, 1, 1> MB;
+
+public:
+ typedef Dune::BCRSMatrix<MB> type;
+};
+SET_PROP(DecoupledTwoPTwoC, PressureRHSVector)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+public:
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > type;
+};
+
+SET_PROP(DecoupledTwoPTwoC, PressurePreconditioner)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureCoefficientMatrix)) Matrix;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureRHSVector)) Vector;
+public:
+ typedef Dune::SeqILUn<Matrix, Vector, Vector> type;
+};
+SET_PROP(DecoupledTwoPTwoC, PressureSolver)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureCoefficientMatrix)) Matrix;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureRHSVector)) Vector;
+public:
+ typedef Dune::BiCGSTABSolver<Vector> type;
+};
+//SET_INT_PROP(DecoupledTwoPTwoC, PressurePreconditioner, SolverIndices::seqILU0);
+//SET_INT_PROP(DecoupledTwoPTwoC, PressureSolver, SolverIndices::biCGSTAB);
+SET_SCALAR_PROP(DecoupledTwoPTwoC, ReductionSolver, 1E-12);
+SET_INT_PROP(DecoupledTwoPTwoC, MaxIterationNumberSolver, 10000);
+SET_INT_PROP(DecoupledTwoPTwoC, IterationNumberPreconditioner, 0);
+SET_INT_PROP(DecoupledTwoPTwoC, VerboseLevelSolver, 1);
+SET_SCALAR_PROP(DecoupledTwoPTwoC, RelaxationPreconditioner, 1.0);
+SET_PROP(DecoupledTwoPTwoC, SolverParameters)
+{
+public:
+ //solver parameters
+ static const double reductionSolver = GET_PROP_VALUE(TypeTag, PTAG(ReductionSolver));
+ static const int maxIterationNumberSolver = GET_PROP_VALUE(TypeTag, PTAG(MaxIterationNumberSolver));
+ static const int iterationNumberPreconditioner = GET_PROP_VALUE(TypeTag, PTAG(IterationNumberPreconditioner));
+ static const int verboseLevelSolver = GET_PROP_VALUE(TypeTag, PTAG(VerboseLevelSolver));
+ static const double relaxationPreconditioner = GET_PROP_VALUE(TypeTag, PTAG(RelaxationPreconditioner));
+};
+
+}
+
+/*!
+ * \brief The common indices for the 2p2c are the same as for the two-phase model.
+ */
+template <class TypeTag>
+struct TwoPCommonIndices
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ // Formulations
+ //saturation flags
+ static const int saturationW = 0;
+ static const int saturationNW = 1;
+ //pressure flags
+ static const int pressureW = 0;
+ static const int pressureNW = 1;
+ static const int pressureGlobal = 2;
+ //velocity flags
+ static const int velocityW = 0;
+ static const int velocityNW = 1;
+ static const int velocityTotal = 2;
+
+ // Phase indices
+ static const int wPhaseIdx = FluidSystem::wPhaseIdx; //!< Index of the wetting phase in a phase vector
+ static const int nPhaseIdx = FluidSystem::nPhaseIdx; //!< Index of the non-wetting phase in a phase vector
+
+ // BoundaryCondition flags
+ static const int satDependent = 0;
+ static const int permDependent = 1;
+};
+
+// \}
+
+}
+
+#endif
diff --git a/dumux/decoupled/2p2c/boundaryconditions2p2c.hh b/dumux/decoupled/2p2c/boundaryconditions2p2c.hh
new file mode 100644
index 0000000000000000000000000000000000000000..34affed587933f51d2858a6196816c94215c0bb9
--- /dev/null
+++ b/dumux/decoupled/2p2c/boundaryconditions2p2c.hh
@@ -0,0 +1,45 @@
+// $Id: boundaryconditions2p2c.hh 3357 2010-03-25 13:02:05Z lauser $
+/*****************************************************************************
+ * Copyright (C) 2007-2008 by Jochen Fritz *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+
+#ifndef DUMUX_BOUNDARYCONDITIONS2P2C_HH
+#define DUMUX_BOUNDARYCONDITIONS2P2C_HH
+
+
+namespace Dumux
+{
+/** \ingroup IMPEC
+ * \brief Defines type of boundary conditions for 2p2c processes
+ *
+ * This is to distinguish BC types for 2p2c processes similar to
+ * the class Dumux::BoundaryConditions which distinguishes between
+ * dirichlet, process and neumann.
+ * BoundaryConditions for compositional transport can either be
+ * defined as saturations or as total concentrations (decoupled method).
+ * Each leads via pressure and constitutive relationships to the other
+ * and to the phase concentrations.
+ */
+struct BoundaryConditions2p2c
+{
+ enum Flags
+ {
+ saturation=1,
+ concentration=2,
+ };
+};
+
+}
+
+#endif
diff --git a/dumux/decoupled/2p2c/dec2p2cfluidstate.hh b/dumux/decoupled/2p2c/dec2p2cfluidstate.hh
new file mode 100644
index 0000000000000000000000000000000000000000..65498170bf3f8e65862d7a6e3671036e56334f96
--- /dev/null
+++ b/dumux/decoupled/2p2c/dec2p2cfluidstate.hh
@@ -0,0 +1,379 @@
+/*****************************************************************************
+ * Copyright (C) 2009 by Benjamin Faigle *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Calculates the 2p2c phase state for compositional models.
+ */
+#ifndef DUMUX_DEC2P2C_FLUID_STATE_HH
+#define DUMUX_DEC2P2C_FLUID_STATE_HH
+
+#include <dumux/material/fluidstate.hh>
+#include <dumux/decoupled/2p2c/2p2cproperties.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup multiphysics multiphase
+ * \brief Calculates the phase state from the primary variables in the
+ * sequential 2p2c model.
+ *
+ * This boils down to so-called "flash calculation", in this case isothermal and isobaric.
+ *
+ * \tparam TypeTag The property Type Tag
+ */
+template <class TypeTag>
+class DecTwoPTwoCFluidState : public FluidState<typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)),
+ DecTwoPTwoCFluidState<TypeTag> >
+{
+ typedef DecTwoPTwoCFluidState<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(TwoPIndices)) Indices;
+ static const int pressureType = GET_PROP_VALUE(TypeTag, PTAG(PressureFormulation)); //!< gives kind of pressure used (\f$ 0 = p_w\f$, \f$ 1 = p_n\f$, \f$ 2 = p_{global}\f$)
+
+
+ enum {
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+
+ wCompIdx = Indices::wPhaseIdx,
+ nCompIdx = Indices::nPhaseIdx,
+ };
+
+public:
+ enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
+ numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),};
+
+public:
+ /*!
+ * \name flash calculation routines
+ * Routines to determine the phase composition after the transport step.
+ */
+ //@{
+
+ //! the update routine equals the 2p2c - concentration flash for constant p & t.
+ /*!
+ * Routine goes as follows:
+ * - determination of the equilibrium constants from the fluid system
+ * - determination of maximum solubilities (mole fractions) according to phase pressures
+ * - comparison with Z1 to determine phase presence => phase mass fractions
+ * - round off fluid properties
+ * \param Z1 Feed mass fraction [-]
+ * \param phasePressure Vector holding the pressure of the phases [Pa]
+ * \param poro Porosity [-]
+ * \param temperature Temperature [K]
+ */
+ void update(Scalar Z1, Scalar pw, Scalar poro, Scalar temperature)
+ {
+ if (pressureType == Indices::pressureGlobal)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Pressure type not supported in fluidState!");
+ }
+ else
+ phasePressure_[wPhaseIdx] = pw;
+ phasePressure_[nPhaseIdx] = pw; // as long as capillary pressure is neglected
+ temperature_=temperature;
+
+
+ //mole equilibrium ratios K for in case wPhase is reference phase
+ double k1 = FluidSystem::activityCoeff(wPhaseIdx, wCompIdx, temperature_, phasePressure_[nPhaseIdx], *this)
+ / phasePressure_[nPhaseIdx]; // = p^wComp_vap / p_n
+ double k2 = FluidSystem::activityCoeff(wPhaseIdx, nCompIdx, temperature_, phasePressure_[nPhaseIdx], *this)
+ / phasePressure_[nPhaseIdx]; // = H^nComp_w / p_n
+
+ // get mole fraction from equilibrium konstants
+ Scalar xw1 = (1. - k2) / (k1 -k2);
+ Scalar xn1 = xw1 * k1;
+
+
+ // transform mole to mass fractions
+ massfrac_[wPhaseIdx][wCompIdx] = xw1 * FluidSystem::molarMass(wCompIdx)
+ / ( xw1 * FluidSystem::molarMass(wCompIdx) + (1.-xw1) * FluidSystem::molarMass(nCompIdx) );
+ massfrac_[nPhaseIdx][wCompIdx] = xn1 * FluidSystem::molarMass(wCompIdx)
+ / ( xn1 * FluidSystem::molarMass(wCompIdx) + (1.-xn1) * FluidSystem::molarMass(nCompIdx) );
+
+
+ //mass equilibrium ratios
+ equilRatio_[nPhaseIdx][wCompIdx] = massfrac_[nPhaseIdx][wCompIdx] / massfrac_[wPhaseIdx][wCompIdx]; // = Xn1 / Xw1 = K1
+ equilRatio_[nPhaseIdx][nCompIdx] = (1.-massfrac_[nPhaseIdx][wCompIdx])/ (1.-massfrac_[wPhaseIdx][wCompIdx]); // =(1.-Xn1) / (1.-Xw1) = K2
+ equilRatio_[wPhaseIdx][nCompIdx] = equilRatio_[wPhaseIdx][wCompIdx] = 1.;
+
+ // phase fraction of nPhase [mass/totalmass]
+ nu_[nPhaseIdx] = 0;
+
+ // check if there is enough of component 1 to form a phase
+ if (Z1 > massfrac_[nPhaseIdx][wCompIdx] && Z1 < massfrac_[wPhaseIdx][wCompIdx])
+ nu_[nPhaseIdx] = -((equilRatio_[nPhaseIdx][wCompIdx]-1)*Z1 + (equilRatio_[nPhaseIdx][nCompIdx]-1)*(1-Z1)) / (equilRatio_[nPhaseIdx][wCompIdx]-1) / (equilRatio_[nPhaseIdx][nCompIdx] -1);
+ else if (Z1 <= massfrac_[nPhaseIdx][wCompIdx]) // too little wComp to form a phase
+ {
+ nu_[nPhaseIdx] = 1; // only nPhase
+ massfrac_[nPhaseIdx][wCompIdx] = Z1; // hence, assign complete mass soluted into nPhase
+ }
+ else // (Z1 >= Xw1) => no nPhase
+ {
+ nu_[nPhaseIdx] = 0; // no second phase
+ massfrac_[wPhaseIdx][wCompIdx] = Z1;
+ }
+
+ // complete array of mass fractions
+ massfrac_[wPhaseIdx][nCompIdx] = 1. - massfrac_[wPhaseIdx][wCompIdx];
+ massfrac_[nPhaseIdx][nCompIdx] = 1. - massfrac_[nPhaseIdx][wCompIdx];
+
+ // complete phase mass fractions
+ nu_[wPhaseIdx] = 1. - nu_[nPhaseIdx];
+
+// // let FluidSystem compute partial pressures
+// FluidSystem::computePartialPressures(temperature_, phasePressure_[nPhaseIdx], *this);
+
+ // get densities with correct composition
+ density_[wPhaseIdx] = FluidSystem::phaseDensity(wPhaseIdx,
+ temperature,
+ phasePressure_[wPhaseIdx],
+ *this);
+ density_[nPhaseIdx] = FluidSystem::phaseDensity(nPhaseIdx,
+ temperature,
+ phasePressure_[nPhaseIdx],
+ *this);
+
+ Sw_ = (nu_[wPhaseIdx]) / density_[wPhaseIdx];
+ Sw_ /= (nu_[wPhaseIdx]/density_[wPhaseIdx] + nu_[nPhaseIdx]/density_[nPhaseIdx]);
+
+ massConcentration_[wCompIdx] =
+ poro * (massfrac_[wPhaseIdx][wCompIdx] * Sw_ * density_[wPhaseIdx]
+ + massfrac_[nPhaseIdx][wCompIdx] * (1.-Sw_) * density_[nPhaseIdx]);
+ massConcentration_[nCompIdx] =
+ poro * (massfrac_[wPhaseIdx][nCompIdx] * Sw_ * density_[wPhaseIdx]
+ + massfrac_[nPhaseIdx][nCompIdx] * (1-Sw_) * density_[nPhaseIdx]);
+ }
+
+ //! a flash routine for 2p2c if the saturation instead of total concentration is known.
+ /*!
+ * Routine goes as follows:
+ * - determination of the equilibrium constants from the fluid system
+ * - determination of maximum solubilities (mole fractions) according to phase pressures
+ * - round off fluid properties
+ * \param sat saturation of phase 1 [-]
+ * \param phasePressure Vector holding the pressure of the phases [Pa]
+ * \param poro Porosity [-]
+ * \param temperature Temperature [K]
+ */
+ void satFlash(Scalar sat, Scalar pw, Scalar poro, Scalar temperature)
+ {
+ if (pressureType == Indices::pressureGlobal)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Pressure type not supported in fluidState!");
+ }
+// else if (sat <= 0 || sat >= 1)
+// DUNE_THROW(Dune::RangeError,
+// "Decoupled2p2c :: saturation initial and boundary conditions may not equal zero or one!");
+
+ // assign values
+ Sw_ = sat;
+ phasePressure_[wPhaseIdx] = pw;
+ phasePressure_[nPhaseIdx] = pw; // as long as capillary pressure is neglected
+ temperature_=temperature;
+ nu_[nPhaseIdx] = nu_[wPhaseIdx] = NAN; //in contrast to the standard update() method, satflash() does not calculate nu.
+
+
+ //mole equilibrium ratios K for in case wPhase is reference phase
+ double k1 = FluidSystem::activityCoeff(wPhaseIdx, wCompIdx, temperature_, phasePressure_[nPhaseIdx], *this)
+ / phasePressure_[nPhaseIdx];
+ double k2 = FluidSystem::activityCoeff(wPhaseIdx, nCompIdx, temperature_, phasePressure_[nPhaseIdx], *this)
+ / phasePressure_[nPhaseIdx];
+
+ if (Sw_==1.) //only wPhase present
+ k1 = 1.;
+ if (Sw_==0.) //only nPhase present
+ k2 = 0.;
+
+ // get mole fraction from equilibrium konstants
+ Scalar xw1 = (1. - k2) / (k1 -k2);
+ Scalar xn1 = xw1 * k1;
+
+
+ // transform mole to mass fractions
+ massfrac_[wPhaseIdx][wCompIdx] = xw1 * FluidSystem::molarMass(wCompIdx)
+ / ( xw1 * FluidSystem::molarMass(wCompIdx) + (1.-xw1) * FluidSystem::molarMass(nCompIdx) );
+ massfrac_[nPhaseIdx][wCompIdx] = xn1 * FluidSystem::molarMass(wCompIdx)
+ / ( xn1 * FluidSystem::molarMass(wCompIdx) + (1.-xn1) * FluidSystem::molarMass(nCompIdx) );
+
+ // complete array of mass fractions
+ massfrac_[wPhaseIdx][nCompIdx] = 1. - massfrac_[wPhaseIdx][wCompIdx];
+ massfrac_[nPhaseIdx][nCompIdx] = 1. - massfrac_[nPhaseIdx][wCompIdx];
+
+ //mass equilibrium ratios
+ equilRatio_[nPhaseIdx][wCompIdx] = massfrac_[nPhaseIdx][wCompIdx] / massfrac_[wPhaseIdx][wCompIdx]; // = Xn1 / Xw1 = K1
+ equilRatio_[nPhaseIdx][nCompIdx] = (1.-massfrac_[nPhaseIdx][wCompIdx])/ (1.-massfrac_[wPhaseIdx][wCompIdx]); // =(1.-Xn1) / (1.-Xw1) = K2
+ equilRatio_[wPhaseIdx][nCompIdx] = equilRatio_[wPhaseIdx][wCompIdx] = 1.;
+
+ // get densities with correct composition
+ density_[wPhaseIdx] = FluidSystem::phaseDensity(wPhaseIdx,
+ temperature,
+ phasePressure_[wPhaseIdx],
+ *this);
+ density_[nPhaseIdx] = FluidSystem::phaseDensity(nPhaseIdx,
+ temperature,
+ phasePressure_[nPhaseIdx],
+ *this);
+
+ massConcentration_[wCompIdx] =
+ poro * (massfrac_[wPhaseIdx][wCompIdx] * Sw_ * density_[wPhaseIdx]
+ + massfrac_[nPhaseIdx][wCompIdx] * (1.-Sw_) * density_[nPhaseIdx]);
+ massConcentration_[nCompIdx] =
+ poro * (massfrac_[wPhaseIdx][nCompIdx] * Sw_ * density_[wPhaseIdx]
+ + massfrac_[nPhaseIdx][nCompIdx] * (1-Sw_) * density_[nPhaseIdx]);
+ }
+ //@}
+ /*!
+ * \name acess functions
+ */
+ //@{
+ /*!
+ * \brief Returns the saturation of a phase.
+ *
+ * \param phaseIdx the index of the phase
+ */
+ Scalar saturation(int phaseIdx) const
+ {
+ if (phaseIdx == wPhaseIdx)
+ return Sw_;
+ else
+ return Scalar(1.0) - Sw_;
+ };
+
+ /*!
+ * \brief Returns the mass fraction of a component in a phase.
+ *
+ * \param phaseIdx the index of the phase
+ * \param compIdx the index of the component
+ */
+ Scalar massFrac(int phaseIdx, int compIdx) const
+ {
+ return massfrac_[phaseIdx][compIdx];
+
+ }
+
+ /*!
+ * \brief Returns the molar fraction of a component in a fluid phase.
+ *
+ * \param phaseIdx the index of the phase
+ * \param compIdx the index of the component
+ */
+ Scalar moleFrac(int phaseIdx, int compIdx) const
+ {
+ // as the moass fractions are calculated, it is used to determine the mole fractions
+ double moleFrac_ = ( massfrac_[phaseIdx][compIdx] / FluidSystem::molarMass(compIdx) ); // = moles of compIdx
+ moleFrac_ /= ( massfrac_[phaseIdx][wCompIdx] / FluidSystem::molarMass(wCompIdx)
+ + massfrac_[phaseIdx][nCompIdx] / FluidSystem::molarMass(nCompIdx) ); // /= total moles in phase
+ return moleFrac_;
+ }
+
+ /*!
+ * \brief Returns the total mass concentration of a component [kg / m^3].
+ *
+ * This is equivalent to the sum of the component concentrations for all
+ * phases multiplied with the phase density.
+ *
+ * \param compIdx the index of the component
+ */
+ Scalar massConcentration(int compIdx) const
+ {
+ return massConcentration_[compIdx];
+ };
+
+
+ /*!
+ * \brief Returns the density of a phase [kg / m^3].
+ *
+ * \param phaseIdx the index of the phase
+ */
+ Scalar density(int phaseIdx) const
+ { return density_[phaseIdx]; }
+
+ /*!
+ * \brief Return the partial pressure of a component in the gas phase.
+ *
+ * For an ideal gas, this means R*T*c.
+ * Unit: [Pa] = [N/m^2]
+ *
+ * \param compIdx the index of the component
+ */
+ Scalar partialPressure(int componentIdx) const
+ {
+ if(componentIdx==nCompIdx)
+ return phasePressure_[nPhaseIdx]*moleFrac(nPhaseIdx, nCompIdx);
+ if(componentIdx == wCompIdx)
+ return phasePressure_[nPhaseIdx]*moleFrac(nPhaseIdx, wCompIdx);
+ else
+ DUNE_THROW(Dune::NotImplemented, "component not found in fluidState!");
+ }
+
+ /*!
+ * \brief Returns the pressure of a fluid phase [Pa].
+ *
+ * \param phaseIdx the index of the phase
+ */
+ Scalar phasePressure(int phaseIdx) const
+ { return phasePressure_[phaseIdx]; }
+
+ /*!
+ * \brief Returns the capillary pressure [Pa]
+ */
+ Scalar capillaryPressure() const
+ { return phasePressure_[nPhaseIdx] - phasePressure_[wPhaseIdx]; }
+
+ /*!
+ * \brief Returns the temperature of the fluids [K].
+ *
+ * Note that we assume thermodynamic equilibrium, so all fluids
+ * and the rock matrix exhibit the same temperature.
+ */
+ Scalar temperature() const
+ { return temperature_; };
+
+ /*!
+ * \brief Returns the phase mass fraction. phase mass per total mass [kg/kg].
+ *
+ * \param phaseIdx the index of the phase
+ */
+ Scalar phaseMassFraction(int phaseIdx)
+ {
+ if (std::isnan(nu_[phaseIdx])) //in contrast to the standard update() method, satflash() does not calculate nu.
+ {
+ nu_[wPhaseIdx] = Sw_ * density_[wPhaseIdx] / (Sw_ * density_[wPhaseIdx] + (1. - Sw_) * density_[nPhaseIdx]);
+ nu_[nPhaseIdx] = 1. - nu_[wPhaseIdx];
+ return nu_[phaseIdx];
+ }
+ else
+ return nu_[phaseIdx];
+ }
+ //@}
+
+private:
+ Scalar density_[numPhases];
+ Scalar massConcentration_[numComponents];
+ Scalar massfrac_[numPhases][numComponents];
+ Scalar equilRatio_[numPhases][numComponents];
+ Scalar phasePressure_[numPhases];
+ Scalar temperature_;
+ Scalar Sw_;
+ Scalar nu_[numPhases]; //phase mass fraction
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/decoupled/2p2c/fvpressure2p2c.hh b/dumux/decoupled/2p2c/fvpressure2p2c.hh
new file mode 100644
index 0000000000000000000000000000000000000000..427dbd449cccfda84257c72e37344de115212dff
--- /dev/null
+++ b/dumux/decoupled/2p2c/fvpressure2p2c.hh
@@ -0,0 +1,1350 @@
+// $Id:
+/*****************************************************************************
+ * Copyright (C) 2007-2009 by Bernd Flemisch *
+ * Copyright (C) 2007-2009 by Jochen Fritz *
+ * Copyright (C) 2008-2009 by Markus Wolff *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_FVPRESSURE2P2C_HH
+#define DUMUX_FVPRESSURE2P2C_HH
+
+// dune environent:
+#include <dune/istl/bvector.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
+
+// dumux environment
+#include "dumux/common/pardiso.hh"
+#include "dumux/common/math.hh"
+#include <dumux/decoupled/2p2c/2p2cproperties.hh>
+
+/**
+ * @file
+ * @brief Finite Volume Diffusion Model
+ * @author Bernd Flemisch, Jochen Fritz, Markus Wolff
+ */
+
+namespace Dumux
+{
+//! The finite volume model for the solution of the compositional pressure equation
+/*! \ingroup multiphase
+ * Provides a Finite Volume implementation for the pressure equation of a gas-liquid
+ * system with two components. An IMPES-like method is used for the sequential
+ * solution of the problem. Capillary forces and diffusion are
+ * neglected. Isothermal conditions and local thermodynamic
+ * equilibrium are assumed. Gravity is included.
+ * \f[
+ c_{total}\frac{\partial p}{\partial t} + \sum_{\kappa} \frac{\partial v_{total}}{\partial C^{\kappa}} \nabla \cdot \left( \sum_{\alpha} X^{\kappa}_{\alpha} \varrho_{alpha} \bf{v}_{\alpha}\right)
+ = \sum_{\kappa} \frac{\partial v_{total}}{\partial C^{\kappa}} q^{\kappa},
+ * \f]
+ * where \f$\bf{v}_{\alpha} = - \lambda_{\alpha} \bf{K} \left(\nabla p_{\alpha} + \rho_{\alpha} \bf{g} \right)\f$.
+ * \f$ c_{total} \f$ represents the total compressibility, for constant porosity this yields \f$ - \frac{\partial V_{total}}{\partial p_{\alpha}}\f$,
+ * \f$p_{\alpha} \f$ denotes the phase pressure, \f$ \bf{K} \f$ the absolute permeability, \f$ \lambda_{\alpha} \f$ the phase mobility,
+ * \f$ \rho_{\alpha} \f$ the phase density and \f$ \bf{g} \f$ the gravity constant and \f$ C^{\kappa} \f$ the total Component concentration.
+ * See paper SPE 99619 or "Analysis of a Compositional Model for Fluid
+ * Flow in Porous Media" by Chen, Qin and Ewing for derivation.
+ *
+ * The partial derivatives of the actual fluid volume \f$ v_{total} \f$ are gained by using a secant method.
+ *
+ * \tparam TypeTag The Type Tag
+ */
+template<class TypeTag> class FVPressure2P2C
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) ReferenceElements;
+ typedef typename ReferenceElements::Container ReferenceElementContainer;
+ typedef typename ReferenceElements::ContainerFaces ReferenceElementFaceContainer;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+ typedef typename SpatialParameters::MaterialLaw MaterialLaw;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+ enum
+ {
+ pw = Indices::pressureW,
+ pn = Indices::pressureNW,
+ pglobal = Indices::pressureGlobal,
+ Sw = Indices::saturationW,
+ Sn = Indices::saturationNW,
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
+ wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx
+ };
+
+ // typedefs to abbreviate several dune classes...
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::Grid Grid;
+ typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ // convenience shortcuts for Vectors/Matrices
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix;
+ typedef Dune::FieldVector<Scalar, 2> PhaseVector;
+
+ // the typenames used for the stiffness matrix and solution vector
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureCoefficientMatrix)) Matrix;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureRHSVector)) RHSVector;
+
+ //initializes the matrix to store the system of equations
+ void initializeMatrix();
+
+ //function which assembles the system of equations to be solved
+ void assemble(bool first);
+
+ //solves the system of equations to get the spatial distribution of the pressure
+ void solve();
+
+protected:
+ Problem& problem()
+ {
+ return problem_;
+ }
+ const Problem& problem() const
+ {
+ return problem_;
+ }
+
+public:
+ //the variables object is initialized, non-compositional before and compositional after first pressure calculation
+ void initialMaterialLaws(bool compositional);
+
+ //constitutive functions are initialized and stored in the variables object
+ void updateMaterialLaws();
+
+ //initialization routine to prepare first timestep
+ void initialize(bool solveTwice = false);
+
+ //pressure solution routine: update estimate for secants, assemble, solve.
+ void pressure(bool solveTwice = true)
+ {
+ //pre-transport to estimate update vector
+ Scalar dt_estimate = 0.;
+ Dune::dinfo << "secant guess"<< std::endl;
+ problem_.transportModel().update(-1, dt_estimate, problem_.variables().updateEstimate(), false);
+ //last argument false in update() makes shure that this is estimate and no "real" transport step
+ problem_.variables().updateEstimate() *= problem_.timeManager().timeStepSize();
+
+ assemble(false); Dune::dinfo << "pressure calculation"<< std::endl;
+ solve();
+
+ return;
+ }
+
+ void calculateVelocity()
+ {
+ // TODO: do this if we use a total velocity description. else, its just a vaste of time.
+ return;
+ }
+
+ //numerical volume derivatives wrt changes in mass, pressure
+ void volumeDerivatives(GlobalPosition globalPos, ElementPointer ep, Scalar& dv_dC1, Scalar& dv_dC2, Scalar& dV_dp);
+
+ /*! \name general methods for serialization, output */
+ //@{
+ // serialization methods
+ template<class Restarter>
+ void serialize(Restarter &res)
+ {
+ return;
+ }
+
+ template<class Restarter>
+ void deserialize(Restarter &res)
+ {
+ return;
+ }
+
+ //! \brief Write data files
+ /* \param name file name */
+ template<class MultiWriter>
+ void addOutputVtkFields(MultiWriter &writer)
+ {
+ problem().variables().addOutputVtkFields(writer);
+
+#if DUNE_MINIMAL_DEBUG_LEVEL <= 3
+ // add debug stuff
+ Dune::BlockVector<Dune::FieldVector<double,1> > *errorCorrPtr = writer.template createField<double, 1> (dV_dp.size());
+ *errorCorrPtr = errorCorrection;
+// int size = subdomainPtr.size();
+ writer.addCellData(errorCorrPtr, "Error Correction");
+
+ Dune::BlockVector<Dune::FieldVector<double,1> > *dV_dpPtr = writer.template createField<double, 1> (dV_dp.size());
+ *dV_dpPtr = dV_dp;
+ writer.addCellData(dV_dpPtr, "dV_dP");
+ Dune::BlockVector<Dune::FieldVector<double,1> > *dV_dC1Ptr = writer.template createField<double, 1> (dV_dp.size());
+ Dune::BlockVector<Dune::FieldVector<double,1> > *dV_dC2Ptr = writer.template createField<double, 1> (dV_dp.size());
+ *dV_dC1Ptr = dV_[0];
+ *dV_dC2Ptr = dV_[1];
+ writer.addCellData(dV_dC1Ptr, "dV_dC1");
+ writer.addCellData(dV_dC2Ptr, "dV_dC2");
+
+ Dune::BlockVector<Dune::FieldVector<double,1> > *updEstimate1 = writer.template createField<double, 1> (dV_dp.size());
+ Dune::BlockVector<Dune::FieldVector<double,1> > *updEstimate2 = writer.template createField<double, 1> (dV_dp.size());
+ *updEstimate1 = problem_.variables().updateEstimate()[0];
+ *updEstimate2 = problem_.variables().updateEstimate()[1];
+ writer.addCellData(updEstimate1, "updEstimate comp 1");
+ writer.addCellData(updEstimate2, "updEstimate comp 2");
+#endif
+ return;
+ }
+
+ //! \brief Write additional debug info in a special writer
+ void debugOutput(double pseudoTS = 0.)
+ {
+ std::cout << "Writing debug for current time step\n";
+ debugWriter_.beginTimestep(problem_.timeManager().time() + pseudoTS,
+ problem_.gridView());
+ problem().variables().addOutputVtkFields(debugWriter_);
+
+ debugWriter_.endTimestep();
+ return;
+ }
+ //@}
+
+ //! Constructs a FVPressure2P2C object
+ /**
+ * \param problem a problem class object
+ */
+ FVPressure2P2C(Problem& problem) :
+ problem_(problem), A_(problem.variables().gridSize(), problem.variables().gridSize(), (2 * dim + 1)
+ * problem.variables().gridSize(), Matrix::random), f_(problem.variables().gridSize()),
+ debugWriter_("debugOutput2p2c"), gravity(problem.gravity())
+ {
+ if (pressureType != pw && pressureType != pn && pressureType != pglobal)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
+ }
+ if (saturationType != Sw && saturationType != Sn)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
+ }
+
+ //rezise block vectors
+ dV_[wPhaseIdx].resize(problem_.variables().gridSize());
+ dV_[nPhaseIdx].resize(problem_.variables().gridSize());
+ dV_dp.resize(problem_.variables().gridSize());
+
+ errorCorrection.resize(problem_.variables().gridSize());
+ //prepare stiffness Matrix and Vectors
+ initializeMatrix();
+ }
+
+private:
+ Problem& problem_;
+ Matrix A_;
+ RHSVector f_;
+ std::string solverName_;
+ std::string preconditionerName_;
+
+ //vectors for partial derivatives
+ typename SolutionTypes::PhaseProperty dV_;
+ typename SolutionTypes::ScalarSolution dV_dp;
+
+ // debug
+ Dumux::VtkMultiWriter<GridView> debugWriter_;
+ typename SolutionTypes::ScalarSolution errorCorrection; //debug output
+protected:
+ const Dune::FieldVector<Scalar, dimWorld>& gravity; //!< vector including the gravity constant
+ static const Scalar cFLFactor_ = GET_PROP_VALUE(TypeTag, PTAG(CFLFactor)); //!< determines the CFLfactor
+ static const int pressureType = GET_PROP_VALUE(TypeTag, PTAG(PressureFormulation)); //!< gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
+ static const int saturationType = GET_PROP_VALUE(TypeTag, PTAG(SaturationFormulation)); //!< gives kind of saturation used (\f$ 0 = S_w \f$, \f$ 1 = S_n \f$)
+};
+
+//! initializes the matrix to store the system of equations
+template<class TypeTag>
+void FVPressure2P2C<TypeTag>::initializeMatrix()
+{
+ // determine matrix row sizes
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // cell index
+ int globalIdxI = problem_.variables().index(*eIt);
+
+ // initialize row size
+ int rowSize = 1;
+
+ // run through all intersections with neighbors
+ IntersectionIterator isItEnd = problem_.gridView().template iend(*eIt);
+ for (IntersectionIterator isIt = problem_.gridView().template ibegin(*eIt); isIt != isItEnd; ++isIt)
+ if (isIt->neighbor())
+ rowSize++;
+ A_.setrowsize(globalIdxI, rowSize);
+ }
+ A_.endrowsizes();
+
+ // determine position of matrix entries
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // cell index
+ int globalIdxI = problem_.variables().index(*eIt);
+
+ // add diagonal index
+ A_.addindex(globalIdxI, globalIdxI);
+
+ // run through all intersections with neighbors
+ IntersectionIterator isItEnd = problem_.gridView().template iend(*eIt);
+ for (IntersectionIterator isIt = problem_.gridView().template ibegin(*eIt); isIt != isItEnd; ++isIt)
+ if (isIt->neighbor())
+ {
+ // access neighbor
+ ElementPointer outside = isIt->outside();
+ int globalIdxJ = problem_.variables().index(*outside);
+
+ // add off diagonal index
+ A_.addindex(globalIdxI, globalIdxJ);
+ }
+ }
+ A_.endindices();
+
+ return;
+}
+
+//! initializes the simulation run
+/*!
+ * Initializes the simulation to gain the initial pressure field.
+ *
+ * \param solveTwice flag to determine possible iterations of the initialization process
+ */
+template<class TypeTag>
+void FVPressure2P2C<TypeTag>::initialize(bool solveTwice)
+{
+ // initialguess: set saturations, determine visco and mobility for initial pressure equation
+ // at this moment, the pressure is unknown. Hence, dont regard compositional effects.
+ initialMaterialLaws(false); Dune::dinfo << "first saturation guess"<<std::endl; //=J: initialGuess()
+ #if DUNE_MINIMAL_DEBUG_LEVEL <= 3
+ debugOutput();
+ #endif
+ assemble(true); Dune::dinfo << "first pressure guess"<<std::endl;
+ solve();
+ #if DUNE_MINIMAL_DEBUG_LEVEL <= 3
+ debugOutput(1e-6);
+ #endif
+ // update the compositional variables (hence true)
+ initialMaterialLaws(true); Dune::dinfo << "first guess for mass fractions"<<std::endl;//= J: transportInitial()
+
+ // perform concentration update to determine secants
+ Scalar dt_estimate = 0.;
+ problem_.transportModel().update(0., dt_estimate, problem_.variables().updateEstimate(), false); Dune::dinfo << "secant guess"<< std::endl;
+ dt_estimate = std::min ( problem_.timeManager().timeStepSize(), dt_estimate);
+ problem_.variables().updateEstimate() *= dt_estimate;
+ #if DUNE_MINIMAL_DEBUG_LEVEL <= 3
+ debugOutput(2e-6);
+ #endif
+ // pressure calculation
+ assemble(false); Dune::dinfo << "second pressure guess"<<std::endl;
+ solve();
+ // update the compositional variables
+ initialMaterialLaws(true);
+
+ if (solveTwice)
+ {
+ Dune::BlockVector<Dune::FieldVector<Scalar, 1> > pressureOld(this->problem().variables().pressure());
+ Dune::BlockVector<Dune::FieldVector<Scalar, 1> > pressureDiff;
+ Scalar pressureNorm = 1.; //dummy initialization to perform at least 1 iteration
+ int numIter = 1;
+
+ while (pressureNorm > 1e-5 && numIter < 10)
+ {
+ Scalar dt_dummy=0.; // without this dummy, it will never converge!
+ // update for secants
+ problem_.transportModel().update(0., dt_dummy, problem_.variables().updateEstimate(), false); Dune::dinfo << "secant guess"<< std::endl;
+ problem_.variables().updateEstimate() *= dt_estimate;
+ // pressure calculation
+ assemble(false); Dune::dinfo << "pressure guess number "<< numIter <<std::endl;
+ solve();
+ // update the compositional variables
+ initialMaterialLaws(true);
+
+ pressureDiff = pressureOld;
+ pressureDiff -= this->problem().variables().pressure();
+ pressureNorm = pressureDiff.infinity_norm();
+ pressureOld = this->problem().variables().pressure();
+ pressureNorm /= pressureOld.infinity_norm();
+
+ numIter++;
+ }
+ }
+ return;
+}
+
+//! function which assembles the system of equations to be solved
+/** for first == true, this function assembles the matrix and right hand side for
+ * the solution of the pressure field in the same way as in the class FVPressure2P.
+ * for first == false, the approach is changed to \f[-\frac{\partial V}{\partial p}
+ * \frac{\partial p}{\partial t}+\sum_{\kappa}\frac{\partial V}{\partial m^{\kappa}}\nabla\cdot
+ * \left(\sum_{\alpha}C_{\alpha}^{\kappa}\mathbf{v}_{\alpha}\right)
+ * =\sum_{\kappa}\frac{\partial V}{\partial m^{\kappa}}q^{\kappa} \f]. See Paper SPE 99619.
+ * This is done to account for the volume effects which appear when gas and liquid are dissolved in each other.
+ * \param first Flag if pressure field is unknown
+ */
+template<class TypeTag>
+void FVPressure2P2C<TypeTag>::assemble(bool first)
+{
+ // initialization: set matrix A_ to zero
+ A_ = 0;
+ f_ = 0;
+
+ // initialization: set the fluid volume derivatives to zero
+ Scalar dv_dC1(0.), dv_dC2(0.);
+ for (int i = 0; i < problem_.variables().gridSize(); i++)
+ {
+ dV_[wPhaseIdx][i] = 0; // dv_dC1 = dV/dm1
+ dV_[nPhaseIdx][i] = 0; // dv / dC2
+ dV_dp[i] = 0; // dv / dp
+ }
+
+ // determine maximum error to scale error-term
+ Scalar timestep_ = problem_.timeManager().timeStepSize();
+ Scalar maxErr = fabs(problem_.variables().volErr().infinity_norm()) / timestep_;
+
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // cell geometry type
+ Dune::GeometryType gt = eIt->geometry().type();
+
+ // number of Faces of current element
+ int numberOfFaces = ReferenceElementContainer::general(gt).size(1);
+
+
+ // get global coordinate of cell center
+ const GlobalPosition& globalPos = eIt->geometry().center();
+
+ // cell index
+ int globalIdxI = problem_.variables().index(*eIt);
+
+ // cell volume, assume linear map here
+ Scalar volume = eIt->geometry().volume();
+
+ Scalar densityWI = problem_.variables().densityWetting(globalIdxI);
+ Scalar densityNWI = problem_.variables().densityNonwetting(globalIdxI);
+
+ /****************implement sources and sinks************************/
+ Dune::FieldVector<Scalar,2> source(problem_.source(globalPos, *eIt));
+ if(first)
+ {
+ source[wPhaseIdx] /= densityWI;
+ source[nPhaseIdx] /= densityNWI;
+ }
+ else
+ {
+ // derivatives of the fluid volume with respect to concentration of components, or pressure
+ if (dV_[0][globalIdxI] == 0)
+ volumeDerivatives(globalPos, *eIt, dV_[wPhaseIdx][globalIdxI][0], dV_[nPhaseIdx][globalIdxI][0], dV_dp[globalIdxI][0]);
+
+ source[wPhaseIdx] *= dV_[wPhaseIdx][globalIdxI]; // note: dV_[i][1] = dv_dC1 = dV/dm1
+ source[nPhaseIdx] *= dV_[nPhaseIdx][globalIdxI];
+ }
+ f_[globalIdxI] = volume * (source[wPhaseIdx] + source[nPhaseIdx]);
+ /***********************************/
+
+ // get absolute permeability
+ FieldMatrix permeabilityI(problem_.spatialParameters().intrinsicPermeability(globalPos, *eIt));
+
+ // get mobilities and fractional flow factors
+ Scalar lambdaWI = problem_.variables().mobilityWetting(globalIdxI);
+ Scalar lambdaNWI = problem_.variables().mobilityNonwetting(globalIdxI);
+ Scalar fractionalWI=0, fractionalNWI=0;
+ if (first)
+ {
+ fractionalWI = lambdaWI / (lambdaWI+ lambdaNWI);
+ fractionalNWI = lambdaNWI / (lambdaWI+ lambdaNWI);
+ }
+
+ Scalar pcI = problem_.variables().capillaryPressure(globalIdxI);
+
+ // prepare meatrix entries
+ Scalar entry=0.;
+ Scalar rightEntry=0.;
+
+ // iterate over all faces of the cell
+ IntersectionIterator isItEnd = problem_.gridView().template iend(*eIt);
+ for (IntersectionIterator isIt = problem_.gridView().template ibegin(*eIt); isIt != isItEnd; ++isIt)
+ {
+ // get geometry type of face
+ Dune::GeometryType faceGT = isIt->geometryInInside().type();
+
+ // center in face's reference element
+ const Dune::FieldVector<Scalar, dim - 1>& faceLocal = ReferenceElementFaceContainer::general(faceGT).position(0,0);
+
+ int isIndex = isIt->indexInInside();
+
+ // center of face inside volume reference element
+ const LocalPosition localPosFace(0);
+
+ // get normal vector
+ Dune::FieldVector<Scalar, dimWorld> unitOuterNormal = isIt->unitOuterNormal(faceLocal);
+
+ // get face volume
+ Scalar faceArea = isIt->geometry().volume();
+
+ /************* handle interior face *****************/
+ if (isIt->neighbor())
+ {
+ // access neighbor
+ ElementPointer neighborPointer = isIt->outside();
+ int globalIdxJ = problem_.variables().index(*neighborPointer);
+
+ // gemotry info of neighbor
+ Dune::GeometryType neighborGT = neighborPointer->geometry().type();
+ const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
+
+ // distance vector between barycenters
+ Dune::FieldVector<Scalar, dimWorld> distVec = globalPosNeighbor - globalPos;
+
+ // compute distance between cell centers
+ Scalar dist = distVec.two_norm();
+
+ Dune::FieldVector<Scalar, dimWorld> unitDistVec(distVec);
+ unitDistVec /= dist;
+
+ FieldMatrix permeabilityJ
+ = problem_.spatialParameters().intrinsicPermeability(globalPosNeighbor,
+ *neighborPointer);
+
+ // compute vectorized permeabilities
+ FieldMatrix meanPermeability(0);
+ Dumux::harmonicMeanMatrix(meanPermeability, permeabilityI, permeabilityJ);
+
+ Dune::FieldVector<Scalar, dim> permeability(0);
+ meanPermeability.mv(unitDistVec, permeability);
+
+ // get mobilities in neighbor
+ Scalar lambdaWJ = problem_.variables().mobilityWetting(globalIdxJ);
+ Scalar lambdaNWJ = problem_.variables().mobilityNonwetting(globalIdxJ);
+
+ // phase densities in cell in neighbor
+ Scalar densityWJ = problem_.variables().densityWetting(globalIdxJ);
+ Scalar densityNWJ = problem_.variables().densityNonwetting(globalIdxJ);
+
+ Scalar pcJ = problem_.variables().capillaryPressure(globalIdxJ);
+
+ Scalar rhoMeanW = 0.5 * (densityWI + densityWJ);
+ Scalar rhoMeanNW = 0.5 * (densityNWI + densityNWJ);
+
+ // reset potential gradients, density
+ Scalar potentialW = 0;
+ Scalar potentialNW = 0;
+ Scalar densityW = 0;
+ Scalar densityNW = 0;
+
+ if (first) // if we are at the very first iteration we can't calculate phase potentials
+ {
+ // get fractional flow factors in neigbor
+ Scalar fractionalWJ = lambdaWJ / (lambdaWJ+ lambdaNWJ);
+ Scalar fractionalNWJ = lambdaNWJ / (lambdaWJ+ lambdaNWJ);
+
+ // perform central weighting
+ Scalar lambda = (lambdaWI + lambdaWJ) * 0.5 + (lambdaNWI + lambdaNWJ) * 0.5;
+
+ entry = fabs(lambda*faceArea*(permeability*unitOuterNormal)/(dist));
+
+ Scalar factor = (fractionalWI + fractionalWJ) * (rhoMeanW) * 0.5 + (fractionalNWI + fractionalNWJ) * (rhoMeanNW) * 0.5;
+ rightEntry = factor * lambda * faceArea * (permeability * gravity) * (unitOuterNormal * unitDistVec);
+ }
+ else
+ {
+ // determine volume derivatives
+ if (dV_[0][globalIdxJ] == 0)
+ volumeDerivatives(globalPosNeighbor, *neighborPointer, dV_[wPhaseIdx][globalIdxJ][0], dV_[nPhaseIdx][globalIdxJ][0], dV_dp[globalIdxJ][0]);
+ dv_dC1 = (dV_[wPhaseIdx][globalIdxI] + dV_[wPhaseIdx][globalIdxJ]) / 2; // dV/dm1= dV/dC^1
+ dv_dC2 = (dV_[nPhaseIdx][globalIdxI] + dV_[nPhaseIdx][globalIdxJ]) / 2;
+
+ Scalar graddv_dC1 = (dV_[wPhaseIdx][globalIdxJ] - dV_[wPhaseIdx][globalIdxI]) / dist;
+ Scalar graddv_dC2 = (dV_[nPhaseIdx][globalIdxJ] - dV_[nPhaseIdx][globalIdxI]) / dist;
+
+
+// potentialW = problem_.variables().potentialWetting(globalIdxI, isIndex);
+// potentialNW = problem_.variables().potentialNonwetting(globalIdxI, isIndex);
+//
+// densityW = (potentialW > 0.) ? densityWI : densityWJ;
+// densityNW = (potentialNW > 0.) ? densityNWI : densityNWJ;
+//
+// densityW = (potentialW == 0.) ? rhoMeanW : densityW;
+// densityNW = (potentialNW == 0.) ? rhoMeanNW : densityNW;
+ //jochen: central weighting for gravity term
+ densityW = rhoMeanW; densityNW = rhoMeanNW;
+
+ switch (pressureType) //Markus: hab (unitOuterNormal * distVec)/dist hinzugefuegt
+ {
+ case pw:
+ {
+ potentialW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI]
+ - problem_.variables().pressure()[globalIdxJ]) / (dist*dist);
+ potentialNW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI]
+ - problem_.variables().pressure()[globalIdxJ] + pcI - pcJ) / (dist*dist);
+ break;
+ }
+ case pn:
+ {
+ potentialW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI]
+ - problem_.variables().pressure()[globalIdxJ] - pcI + pcJ) / (dist*dist);
+ potentialNW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI]
+ - problem_.variables().pressure()[globalIdxJ]) / (dist*dist);
+ break;
+ }
+ case pglobal:
+ {
+ DUNE_THROW(Dune::NotImplemented, "Global pressure not yet implemented for 2p2c");
+// potentialW = (problem_.variables().pressure()[globalIdxI]
+// - problem_.variables().pressure()[globalIdxJ] - fMeanNW * (pcI - pcJ)) / dist;
+// potentialNW = (problem_.variables().pressure()[globalIdxI]
+// - problem_.variables().pressure()[globalIdxJ] + fMeanW * (pcI - pcJ)) / dist;
+ break;
+ }
+ }
+
+ potentialW += densityW * (unitDistVec * gravity);
+ potentialNW += densityNW * (unitDistVec * gravity);
+
+ //store potentials for further calculations (velocity, saturation, ...)
+ problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialW;
+ problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
+
+ // initialize convenience shortcuts
+ Scalar lambdaW, lambdaN;
+ Scalar dV_w(0.), dV_n(0.); // dV_a = \sum_k \rho_a * dv/dC^k * X^k_a
+ Scalar gV_w(0.), gV_n(0.); // multipaper eq(3.3) zeile 3 analogon dV_w
+
+
+ //do the upwinding of the mobility depending on the phase potentials
+ if (potentialW >= 0.)
+ {
+ dV_w = (dv_dC1 * problem_.variables().wet_X1(globalIdxI) + dv_dC2 * (1. - problem_.variables().wet_X1(globalIdxI)));
+ lambdaW = problem_.variables().mobilityWetting(globalIdxI);
+ gV_w = (graddv_dC1 * problem_.variables().wet_X1(globalIdxI) + graddv_dC2 * (1. - problem_.variables().wet_X1(globalIdxI)));
+ dV_w *= densityWI; gV_w *= densityWI;
+ }
+ else
+ {
+ dV_w = (dv_dC1 * problem_.variables().wet_X1(globalIdxJ) + dv_dC2 * (1. - problem_.variables().wet_X1(globalIdxJ)));
+ lambdaW = problem_.variables().mobilityWetting(globalIdxJ);
+ gV_w = (graddv_dC1 * problem_.variables().wet_X1(globalIdxJ) + graddv_dC2 * (1. - problem_.variables().wet_X1(globalIdxJ)));
+ dV_w *= densityWJ; gV_w *= densityWJ;
+ }
+ if (potentialNW >= 0.)
+ {
+ dV_n = (dv_dC1 * problem_.variables().nonwet_X1(globalIdxI) + dv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxI)));
+ lambdaN = problem_.variables().mobilityNonwetting(globalIdxI);
+ gV_n = (graddv_dC1 * problem_.variables().nonwet_X1(globalIdxI) + graddv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxI)));
+ dV_n *= densityNWI; gV_n *= densityNWI;
+ }
+ else
+ {
+ dV_n = (dv_dC1 * problem_.variables().nonwet_X1(globalIdxJ) + dv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxJ)));
+ lambdaN = problem_.variables().mobilityNonwetting(globalIdxJ);
+ gV_n = (graddv_dC1 * problem_.variables().nonwet_X1(globalIdxJ) + graddv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxJ)));
+ dV_n *= densityNWJ; gV_n *= densityNWJ;
+ }
+
+ //calculate current matrix entry
+ entry = faceArea * (lambdaW * dV_w + lambdaN * dV_n)
+ - volume / numberOfFaces * (lambdaW * gV_w + lambdaN * gV_n); // randintegral - gebietsintegral
+ entry *= fabs((permeability*unitOuterNormal)/(dist)); // TODO: markus nimmt hier statt fabs() ein unitDistVec * unitDistVec
+
+ //calculate right hand side
+ rightEntry = faceArea * (unitOuterNormal * unitDistVec) * (densityW * lambdaW * dV_w + densityNW * lambdaN * dV_n);
+ rightEntry -= volume / numberOfFaces * (densityW * lambdaW * gV_w + densityNW * lambdaN * gV_n);
+ rightEntry *= (permeability * gravity); // = multipaper eq(3.3) zeile 2+3 ohne (p-p_k)
+
+ } // end !first
+
+ //set right hand side
+ f_[globalIdxI] -= rightEntry;
+ // set diagonal entry
+ A_[globalIdxI][globalIdxI] += entry;
+ // set off-diagonal entry
+ A_[globalIdxI][globalIdxJ] = -entry;
+ } // end neighbor
+ /****************************************************/
+
+
+ /************* boundary face ************************/
+ else
+ {
+ // get volume derivatives inside the cell
+ dv_dC1 = dV_[wPhaseIdx][globalIdxI];
+ dv_dC2 = dV_[nPhaseIdx][globalIdxI];
+
+ // center of face in global coordinates
+ const GlobalPosition& globalPosFace = isIt->geometry().global(faceLocal);
+
+ // geometrical information
+ Dune::FieldVector<Scalar, dimWorld> distVec(globalPosFace - globalPos);
+ Scalar dist = distVec.two_norm();
+ Dune::FieldVector<Scalar, dimWorld> unitDistVec(distVec);
+ unitDistVec /= dist;
+
+ //get boundary condition for boundary face center
+ BoundaryConditions::Flags bctype = problem_.bcTypePress(globalPosFace, *isIt);
+
+ /********** Dirichlet Boundary *************/
+ if (bctype == BoundaryConditions::dirichlet)
+ {
+ //permeability vector at boundary
+ Dune::FieldVector<Scalar, dim> permeability(0);
+ permeabilityI.mv(unitDistVec, permeability);
+
+ // create a fluid state for the boundary
+ FluidState BCfluidState;
+
+ Scalar temperatureBC = problem_.temperature(globalPosFace, *eIt);
+
+ //get dirichlet pressure boundary condition
+ Scalar pressBound = problem_.dirichletPress(globalPosFace, *isIt);
+
+ Scalar pressBC = 0.;
+ Scalar pcBound = 0.;
+ if(pressureType==pw)
+ {
+ pressBC = pressBound;
+ }
+ else if(pressureType==pn)
+ {
+ //TODO: take pC from variables or from MaterialLaw?
+ // if the latter, one needs Sw
+ pcBound = problem_.variables().capillaryPressure(globalIdxI);
+ pressBC = pressBound - pcBound;
+ }
+
+ if (first)
+ {
+ Scalar lambda = lambdaWI+lambdaNWI;
+ A_[globalIdxI][globalIdxI] += lambda * faceArea * (permeability * unitOuterNormal) / (dist);
+ Scalar pressBC = problem_.dirichletPress(globalPosFace, *isIt);
+ f_[globalIdxI] += lambda * faceArea * pressBC * (permeability * unitOuterNormal) / (dist);
+ Scalar rightentry = (fractionalWI * densityWI
+ + fractionalNWI * densityNWI)
+ * lambda * faceArea * (unitOuterNormal * unitDistVec) *(permeability*gravity);
+ f_[globalIdxI] -= rightentry;
+ }
+ else //not first
+ {
+ //get boundary condition type for compositional transport
+ BoundaryConditions2p2c::Flags bcform = problem_.bcFormulation(globalPosFace, *isIt);
+ if (bcform == BoundaryConditions2p2c::saturation) // saturation given
+ {
+ Scalar satBound = problem_.dirichletTransport(globalPosFace, *isIt);
+ BCfluidState.satFlash(satBound, pressBC, problem_.spatialParameters().porosity(globalPos, *eIt), temperatureBC);
+ }
+ else if (bcform == Dumux::BoundaryConditions2p2c::concentration) // mass fraction given
+ {
+ Scalar Z1Bound = problem_.dirichletTransport(globalPosFace, *isIt);
+ BCfluidState.update(Z1Bound, pressBC, problem_.spatialParameters().porosity(globalPos, *eIt), temperatureBC);
+ }
+ else // nothing declared at boundary
+ {
+ Scalar satBound = problem_.variables().saturation()[globalIdxI];
+ BCfluidState.satFlash(satBound, pressBC, problem_.spatialParameters().porosity(globalPos, *eIt), temperatureBC);
+ Dune::dwarn << "no boundary saturation/concentration specified on boundary pos " << globalPosFace << std::endl;
+ }
+
+
+ // determine fluid properties at the boundary
+ Scalar lambdaWBound = 0.;
+ Scalar lambdaNWBound = 0.;
+
+ Scalar densityWBound =
+ FluidSystem::phaseDensity(wPhaseIdx, temperatureBC, pressBC, BCfluidState);
+ Scalar densityNWBound =
+ FluidSystem::phaseDensity(nPhaseIdx, temperatureBC, pressBC+pcBound, BCfluidState);
+ Scalar viscosityWBound =
+ FluidSystem::phaseViscosity(wPhaseIdx, temperatureBC, pressBC, BCfluidState);
+ Scalar viscosityNWBound =
+ FluidSystem::phaseViscosity(nPhaseIdx, temperatureBC, pressBC+pcBound, BCfluidState);
+
+ // mobility at the boundary
+ switch (GET_PROP_VALUE(TypeTag, PTAG(BoundaryMobility)))
+ {
+ case Indices::satDependent:
+ {
+ lambdaWBound = BCfluidState.saturation(wPhaseIdx)
+ / viscosityWBound;
+ lambdaNWBound = BCfluidState.saturation(nPhaseIdx)
+ / viscosityNWBound;
+ break;
+ }
+ case Indices::permDependent:
+ {
+ lambdaWBound = MaterialLaw::krw(
+ problem_.spatialParameters().materialLawParams(globalPos, *eIt), BCfluidState.saturation(wPhaseIdx))
+ / viscosityWBound;
+ lambdaNWBound = MaterialLaw::krn(
+ problem_.spatialParameters().materialLawParams(globalPos, *eIt), BCfluidState.saturation(wPhaseIdx))
+ / viscosityNWBound;
+ }
+ }
+
+ Scalar rhoMeanW = 0.5 * (densityWI + densityWBound);
+ Scalar rhoMeanNW = 0.5 * (densityNWI + densityNWBound);
+
+ Scalar potentialW = 0;
+ Scalar potentialNW = 0;
+
+ Scalar densityW = 0;
+ Scalar densityNW = 0;
+
+ if (!first)
+ {
+// potentialW = problem_.variables().potentialWetting(globalIdxI, isIndex);
+// potentialNW = problem_.variables().potentialNonwetting(globalIdxI, isIndex);
+//
+// // do potential upwinding according to last potGradient vs Jochen: central weighting
+// densityW = (potentialW > 0.) ? densityWI : densityWBound;
+// densityNW = (potentialNW > 0.) ? densityNWI : densityNWBound;
+//
+// densityW = (potentialW == 0.) ? rhoMeanW : densityW;
+// densityNW = (potentialNW == 0.) ? rhoMeanNW : densityNW;
+ densityW=rhoMeanW; densityNW=rhoMeanNW;
+
+ //calculate potential gradient
+ switch (pressureType)
+ {
+ case pw:
+ {
+ potentialW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI] - pressBound) / (dist * dist);
+ potentialNW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI] + pcI - pressBound - pcBound)
+ / (dist * dist);
+ break;
+ }
+ case pn:
+ {
+ potentialW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI] - pcI - pressBound + pcBound)
+ / (dist * dist);
+ potentialNW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI] - pressBound) / (dist * dist);
+ break;
+ }
+ case pglobal:
+ {
+ Scalar fractionalWBound = lambdaWBound / (lambdaWBound + lambdaNWBound);
+ Scalar fractionalNWBound = lambdaNWBound / (lambdaWBound + lambdaNWBound);
+ Scalar fMeanW = 0.5 * (fractionalWI + fractionalWBound);
+ Scalar fMeanNW = 0.5 * (fractionalNWI + fractionalNWBound);
+
+ potentialW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI] - pressBound - fMeanNW * (pcI
+ - pcBound)) / (dist * dist);
+ potentialNW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI] - pressBound + fMeanW * (pcI
+ - pcBound)) / (dist * dist);
+ break;
+ }
+ }
+ potentialW += densityW * (unitDistVec * gravity);
+ potentialNW += densityNW * (unitDistVec * gravity);
+
+ //store potential gradients for further calculations
+ problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialW;
+ problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
+ } //end !first
+
+ //do the upwinding of the mobility depending on the phase potentials
+ Scalar lambdaW, lambdaNW;
+ Scalar dV_w, dV_n; // gV_a weglassen, da dV/dc am Rand ortsunabhängig angenommen -> am rand nicht bestimmbar -> nur Randintegral ohne Gebietsintegral
+
+ if (potentialW >= 0.)
+ {
+ densityW = (potentialW == 0) ? rhoMeanW : densityWI;
+ dV_w = (dv_dC1 * problem_.variables().wet_X1(globalIdxI)
+ + dv_dC2 * (1. - problem_.variables().wet_X1(globalIdxI)));
+ dV_w *= densityW;
+ lambdaW = (potentialW == 0) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaWI;
+ }
+ else
+ {
+ densityW = densityWBound;
+ dV_w = (dv_dC1 * BCfluidState.massFrac(wPhaseIdx, wCompIdx)
+ + dv_dC2 * BCfluidState.massFrac(wPhaseIdx, nCompIdx));
+ dV_w *= densityW;
+ lambdaW = lambdaWBound;
+ }
+ if (potentialNW >= 0.)
+ {
+ densityNW = (potentialNW == 0) ? rhoMeanNW : densityNWI;
+ dV_n = (dv_dC1 * problem_.variables().nonwet_X1(globalIdxI) + dv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxI)));
+ dV_n *= densityNW;
+ lambdaNW = (potentialNW == 0) ? 0.5 * (lambdaNWI + lambdaNWBound) : lambdaNWI;
+ }
+ else
+ {
+ densityNW = densityNWBound;
+ dV_n = (dv_dC1 * BCfluidState.massFrac(nPhaseIdx, wCompIdx)
+ + dv_dC2 * BCfluidState.massFrac(nPhaseIdx, nCompIdx));
+ dV_n *= densityNW;
+ lambdaNW = lambdaNWBound;
+ }
+
+ //calculate current matrix entry
+ Scalar entry = (lambdaW * dV_w + lambdaNW * dV_n) * ((permeability * unitDistVec) / dist) * faceArea
+ * (unitOuterNormal * unitDistVec);
+
+ //calculate right hand side
+ Scalar rightEntry = (lambdaW * densityW * dV_w + lambdaNW * densityNW * dV_n) * (permeability * gravity)
+ * faceArea ;
+
+
+
+ // set diagonal entry and right hand side entry
+ A_[globalIdxI][globalIdxI] += entry;
+ f_[globalIdxI] += entry * pressBound;
+ f_[globalIdxI] -= rightEntry * (unitOuterNormal * unitDistVec);
+ } //end of if(first) ... else{...
+ } // end dirichlet
+
+ /**********************************
+ * set neumann boundary condition
+ **********************************/
+ else
+ {
+ Dune::FieldVector<Scalar,2> J = problem_.neumann(globalPosFace, *isIt);
+ if (first)
+ {
+ J[wPhaseIdx] /= densityWI;
+ J[nPhaseIdx] /= densityNWI;
+ }
+ else
+ {
+ J[wPhaseIdx] *= dv_dC1;
+ J[nPhaseIdx] *= dv_dC2;
+ }
+
+ f_[globalIdxI] -= (J[wPhaseIdx] + J[nPhaseIdx]) * faceArea;
+
+ //Assumes that the phases flow in the same direction at the neumann boundary, which is the direction of the total flux!!!
+ //needed to determine the upwind direction in the saturation equation
+ problem_.variables().potentialWetting(globalIdxI, isIndex) = J[wPhaseIdx];
+ problem_.variables().potentialNonwetting(globalIdxI, isIndex) = J[nPhaseIdx];
+ }
+ }
+ } // end all intersections
+
+ // compressibility term
+ if (!first && timestep_ != 0)
+ {
+ Scalar compress_term = dV_dp[globalIdxI] / timestep_;
+
+ A_[globalIdxI][globalIdxI] -= compress_term*volume;
+ f_[globalIdxI] -= problem_.variables().pressure()[globalIdxI] * compress_term * volume;
+
+ if (isnan(compress_term) || isinf(compress_term))
+ DUNE_THROW(Dune::MathError, "Compressibility term leads to NAN matrix entry at index " << globalIdxI);
+
+ if(!GET_PROP_VALUE(TypeTag, PTAG(EnableCompressibility)))
+ DUNE_THROW(Dune::NotImplemented, "Compressibility is switched off???");
+ }
+
+ // error reduction routine: volumetric error is damped and inserted to right hand side
+ // if damping is not done, the solution method gets unstable!
+ problem_.variables().volErr()[globalIdxI] /= timestep_;
+ Scalar erri = fabs(problem_.variables().volErr()[globalIdxI]);
+ Scalar x_lo = 0.6;
+ Scalar x_mi = 0.9;
+ Scalar fac = 0.05;
+ Scalar lofac = 0.;
+ Scalar hifac = 0.;
+ hifac /= fac;
+
+ if (erri*timestep_ > 5e-5)
+ if (erri > x_lo * maxErr)
+ {
+ if (erri <= x_mi * maxErr)
+ f_[globalIdxI] += errorCorrection[globalIdxI]= fac* (1-x_mi*(lofac/fac-1)/(x_lo-x_mi) + (lofac/fac-1)/(x_lo-x_mi)*erri/maxErr)
+ * problem_.variables().volErr()[globalIdxI] * volume;
+ else
+ f_[globalIdxI] += errorCorrection[globalIdxI]= fac * (1 + x_mi - hifac*x_mi/(1-x_mi) + (hifac/(1-x_mi)-1)*erri/maxErr)
+ * problem_.variables().volErr()[globalIdxI] * volume;
+ }
+ } // end grid traversal
+ return;
+}
+
+//! solves the system of equations to get the spatial distribution of the pressure
+template<class TypeTag>
+void FVPressure2P2C<TypeTag>::solve()
+{
+ typedef typename GET_PROP(TypeTag, PTAG(SolverParameters)) SolverParameters;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressurePreconditioner)) Preconditioner;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureSolver)) Solver;
+
+// printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
+// printvector(std::cout, f_, "right hand side", "row", 200, 1, 3);
+
+ Dune::MatrixAdapter<Matrix, RHSVector, RHSVector> op(A_); // make linear operator from A_
+ Dune::InverseOperatorResult result;
+ Scalar reduction = SolverParameters::reductionSolver;
+ int maxItSolver = SolverParameters::maxIterationNumberSolver;
+ int iterPreconditioner = SolverParameters::iterationNumberPreconditioner;
+ int verboseLevelSolver = SolverParameters::verboseLevelSolver;
+ Scalar relaxation = SolverParameters::relaxationPreconditioner;
+
+ if (verboseLevelSolver)
+ std::cout << "FVPressure2P: solve for pressure" << std::endl;
+
+ Preconditioner preconditioner(A_, iterPreconditioner, relaxation);
+ Solver solver(op, preconditioner, reduction, maxItSolver, verboseLevelSolver);
+ solver.apply(problem_.variables().pressure(), f_, result);
+
+ return;
+}
+
+/*!
+ * \brief initializes the fluid distribution and hereby the variables container
+ *
+ * This function equals the method initialguess() and transportInitial() in the old model.
+ * It differs from updateMaterialLaws because there are two possible initial conditions:
+ * saturations and concentration.
+ * \param compositional flag that determines if compositional effects are regarded, i.e.
+ * a reasonable pressure field is known.
+ */
+template<class TypeTag>
+void FVPressure2P2C<TypeTag>::initialMaterialLaws(bool compositional)
+{
+ // initialize the fluid system
+ FluidState fluidState;
+
+ // iterate through leaf grid an evaluate c0 at cell center
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // get geometry type
+ Dune::GeometryType gt = eIt->geometry().type();
+
+ // get global coordinate of cell center
+ GlobalPosition globalPos = eIt->geometry().center();
+
+ // assign an Index for convenience
+ int globalIdx = problem_.variables().index(*eIt);
+
+ // get the temperature
+ Scalar temperature_ = problem_.temperature(globalPos, *eIt);
+
+ // initial conditions
+ Scalar pressW = 0;
+ Scalar pressNW = 0;
+ Scalar sat_0=0.;
+
+ BoundaryConditions2p2c::Flags ictype = problem_.initFormulation(globalPos, *eIt); // get type of initial condition
+
+ if(!compositional) //means that we do the first approximate guess without compositions
+ {
+ // phase pressures are unknown, so start with an exemplary
+ Scalar exemplaryPressure = problem_.referencePressure(globalPos, *eIt);
+ pressW = pressNW = problem_.variables().pressure()[globalIdx] = exemplaryPressure;
+
+ if (ictype == BoundaryConditions2p2c::saturation) // saturation initial condition
+ {
+ sat_0 = problem_.initSat(globalPos, *eIt);
+ fluidState.satFlash(sat_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ else if (ictype == BoundaryConditions2p2c::concentration) // concentration initial condition
+ {
+ Scalar Z1_0 = problem_.initConcentration(globalPos, *eIt);
+ fluidState.update(Z1_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ }
+ else if(compositional) //means we regard compositional effects since we know an estimate pressure field
+ {
+ //determine phase pressures from primary pressure variable
+ switch (pressureType)
+ {
+ case pw:
+ {
+ pressW = problem_.variables().pressure()[globalIdx];
+ pressNW = problem_.variables().pressure()[globalIdx] + problem_.variables().capillaryPressure(globalIdx);
+ break;
+ }
+ case pn:
+ {
+ pressW = problem_.variables().pressure()[globalIdx] - problem_.variables().capillaryPressure(globalIdx);
+ pressNW = problem_.variables().pressure()[globalIdx];
+ break;
+ }
+ }
+
+ if (ictype == Dumux::BoundaryConditions2p2c::saturation) // saturation initial condition
+ {
+ sat_0 = problem_.initSat(globalPos, *eIt);
+ fluidState.satFlash(sat_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ else if (ictype == Dumux::BoundaryConditions2p2c::concentration) // concentration initial condition
+ {
+ Scalar Z1_0 = problem_.initConcentration(globalPos, *eIt);
+ fluidState.update(Z1_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ }
+
+ // initialize densities
+ problem_.variables().densityWetting(globalIdx) = FluidSystem::phaseDensity(wPhaseIdx, temperature_, pressW, fluidState);
+ problem_.variables().densityNonwetting(globalIdx) = FluidSystem::phaseDensity(nPhaseIdx, temperature_, pressNW, fluidState);
+
+ // initialize mass fractions
+ problem_.variables().wet_X1(globalIdx) = fluidState.massFrac(wPhaseIdx, wCompIdx);
+ problem_.variables().nonwet_X1(globalIdx) = fluidState.massFrac(nPhaseIdx, wCompIdx);
+
+
+ // initialize viscosities
+ problem_.variables().viscosityWetting(globalIdx) = FluidSystem::phaseViscosity(wPhaseIdx, temperature_, pressW, fluidState);
+ problem_.variables().viscosityNonwetting(globalIdx) = FluidSystem::phaseViscosity(nPhaseIdx, temperature_, pressNW, fluidState);
+
+ // initialize mobilities
+ problem_.variables().mobilityWetting(globalIdx) = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(globalPos, *eIt), fluidState.saturation(wPhaseIdx))
+ / problem_.variables().viscosityWetting(globalIdx);
+ problem_.variables().mobilityNonwetting(globalIdx) = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(globalPos, *eIt), fluidState.saturation(wPhaseIdx))
+ / problem_.variables().viscosityNonwetting(globalIdx);
+
+ // initialize cell concentration
+ problem_.variables().totalConcentration(globalIdx, wCompIdx) = fluidState.massConcentration(wCompIdx);
+ problem_.variables().totalConcentration(globalIdx, nCompIdx) = fluidState.massConcentration(nCompIdx);
+ problem_.variables().saturation()[globalIdx] = fluidState.saturation(wPhaseIdx);
+ }
+ return;
+}
+
+//! constitutive functions are updated once if new concentrations are calculated and stored in the variables container
+template<class TypeTag>
+void FVPressure2P2C<TypeTag>::updateMaterialLaws()
+{
+ // this method only completes the variables: = old postprocessupdate()
+
+ // instantiate a brandnew fluid state object
+ FluidState fluidState;
+
+ //get timestep for error term
+ Scalar dt = problem_.timeManager().timeStepSize();
+
+ // iterate through leaf grid an evaluate c0 at cell center
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // get geometry type
+ Dune::GeometryType gt = eIt->geometry().type();
+
+ // get global coordinate of cell center
+ GlobalPosition globalPos = eIt->geometry().center();
+
+ int globalIdx = problem_.variables().index(*eIt);
+
+ Scalar temperature_ = problem_.temperature(globalPos, *eIt);
+
+
+ // get the overall mass of component 1 Z1 = C^k / (C^1+C^2) [-]
+ Scalar Z1 = problem_.variables().totalConcentration(globalIdx, wCompIdx) / (problem_.variables().totalConcentration(globalIdx, wCompIdx) + problem_.variables().totalConcentration(globalIdx, nCompIdx));
+
+ //determine phase pressures from primary pressure variable
+ Scalar pressW(0.), pressNW(0.);
+ switch (pressureType)
+ {
+ case pw:
+ {
+ pressW = problem_.variables().pressure()[globalIdx];
+
+ pressNW = problem_.variables().pressure()[globalIdx];
+ break;
+ }
+ case pn:
+ {
+ //todo: check this case for consistency throughout the model!
+ pressNW = problem_.variables().pressure()[globalIdx];
+
+ pressW = problem_.variables().pressure()[globalIdx];
+
+ break;
+ }
+ }
+
+ //complete fluid state
+ fluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+
+ /*************************************
+ * update variables in variableclass
+ *************************************/
+ // initialize saturation
+ problem_.variables().saturation(globalIdx) = fluidState.saturation(wPhaseIdx);
+
+ // initialize pC todo: remove this dummy implementation
+ problem_.variables().capillaryPressure(globalIdx) = 0.0;
+
+
+ // initialize viscosities
+ problem_.variables().viscosityWetting(globalIdx)
+ = FluidSystem::phaseViscosity(wPhaseIdx, temperature_, pressW, fluidState);
+ problem_.variables().viscosityNonwetting(globalIdx)
+ = FluidSystem::phaseViscosity(nPhaseIdx, temperature_, pressNW, fluidState);
+
+ // initialize mobilities
+ problem_.variables().mobilityWetting(globalIdx) =
+ MaterialLaw::krw(problem_.spatialParameters().materialLawParams(globalPos, *eIt), fluidState.saturation(wPhaseIdx))
+ / problem_.variables().viscosityWetting(globalIdx);
+ problem_.variables().mobilityNonwetting(globalIdx) =
+ MaterialLaw::krn(problem_.spatialParameters().materialLawParams(globalPos, *eIt), fluidState.saturation(wPhaseIdx))
+ / problem_.variables().viscosityNonwetting(globalIdx);
+
+ // initialize mass fractions
+ problem_.variables().wet_X1(globalIdx) = fluidState.massFrac(wPhaseIdx, wCompIdx);
+ problem_.variables().nonwet_X1(globalIdx) = fluidState.massFrac(nPhaseIdx, wCompIdx);
+
+ // initialize densities
+ problem_.variables().densityWetting(globalIdx)
+ = FluidSystem::phaseDensity(wPhaseIdx, temperature_, pressW, fluidState);
+ problem_.variables().densityNonwetting(globalIdx)
+ = FluidSystem::phaseDensity(nPhaseIdx, temperature_, pressNW, fluidState);
+
+ problem_.spatialParameters().update(fluidState.saturation(wPhaseIdx), *eIt);
+
+ // determine volume mismatch between actual fluid volume and pore volume
+ Scalar sumConc = (problem_.variables().totalConcentration(globalIdx, wCompIdx)
+ + problem_.variables().totalConcentration(globalIdx, nCompIdx));
+ Scalar massw = sumConc * fluidState.phaseMassFraction(wPhaseIdx);
+ Scalar massn = sumConc * fluidState.phaseMassFraction(nPhaseIdx);
+
+ if ((problem_.variables().densityWetting(globalIdx)*problem_.variables().densityNonwetting(globalIdx)) == 0)
+ DUNE_THROW(Dune::MathError, "Decoupled2p2c::postProcessUpdate: try to divide by 0 density");
+ Scalar vol = massw / problem_.variables().densityWetting(globalIdx) + massn / problem_.variables().densityNonwetting(globalIdx);
+ if (dt != 0)
+ {
+ problem_.variables().volErr()[globalIdx] = (vol - problem_.spatialParameters().porosity(globalPos, *eIt));
+
+ Scalar volErrI = problem_.variables().volErr(globalIdx);
+ if (std::isnan(volErrI))
+ {
+ DUNE_THROW(Dune::MathError, "Decoupled2p2c::postProcessUpdate:\n"
+ << "volErr[" << globalIdx << "] isnan: vol = " << vol
+ << ", massw = " << massw << ", rho_l = " << problem_.variables().densityWetting(globalIdx)
+ << ", massn = " << massn << ", rho_g = " << problem_.variables().densityNonwetting(globalIdx)
+ << ", poro = " << problem_.spatialParameters().porosity(globalPos, *eIt) << ", dt = " << dt);
+ }
+ }
+ else
+ {
+ problem_.variables().volErr()[globalIdx] = 0;
+ }
+ }
+ return;
+}
+
+//! partial derivatives of the volumes w.r.t. changes in total concentration and pressure
+/*!
+ * This method calculates the volume derivatives via a secant method, where the
+ * secants are gained in a pre-computational step via the transport equation and
+ * the last TS size.
+ * The partial derivatives w.r.t. mass are defined as
+ * \f$ \frac{\partial v}{\partial C^{\kappa}} = \frac{\partial V}{\partial m^{\kappa}}\f$
+ *
+ * \param globalPos The global position of the current element
+ * \param ep A pointer to the current element
+ * \param[out] dv_dC1 partial derivative of fluid volume w.r.t. mass of component 1 [m^3/kg]
+ * \param[out] dv_dC2 partial derivative of fluid volume w.r.t. mass of component 2 [m^3/kg]
+ * \param[out] dv_dp partial derivative of fluid volume w.r.t. pressure [1/Pa]
+ */
+template<class TypeTag>
+void FVPressure2P2C<TypeTag>::volumeDerivatives(GlobalPosition globalPos, ElementPointer ep, Scalar& dv_dC1, Scalar& dv_dC2, Scalar& dV_dp)
+{
+ // cell index
+ int globalIdx = problem_.variables().index(*ep);
+
+ // get cell temperature
+ Scalar temperature_ = problem_.temperature(globalPos, *ep);
+
+ // initialize an Fluid state for the update
+ FluidState updFluidState;
+
+ /**********************************
+ * a) get necessary variables
+ **********************************/
+ //determine phase pressures from primary pressure variable
+ Scalar pressW=0.;
+ switch (pressureType)
+ {
+ case pw:
+ {
+ pressW = problem_.variables().pressure()[globalIdx];
+ break;
+ }
+ case pn:
+ {
+ pressW = problem_.variables().pressure()[globalIdx] - problem_.variables().capillaryPressure(globalIdx);
+ break;
+ }
+ }
+
+ Scalar v_w = 1. / problem_.variables().densityWetting(globalIdx);
+ Scalar v_g = 1. / problem_.variables().densityNonwetting(globalIdx);
+ Scalar sati = problem_.variables().saturation()[globalIdx];
+ // mass of components inside the cell
+ Scalar m1 = problem_.variables().totalConcentration(globalIdx, wCompIdx);
+ Scalar m2 = problem_.variables().totalConcentration(globalIdx, nCompIdx);
+ // mass fraction of wetting phase
+ Scalar nuw1 = sati/ v_w / (sati/v_w + (1-sati)/v_g);
+ // actual fluid volume
+ Scalar volalt = (m1+m2) * (nuw1 * v_w + (1-nuw1) * v_g);
+
+ /**********************************
+ * b) define increments
+ **********************************/
+ // increments for numerical derivatives
+ Scalar inc1 = (fabs(problem_.variables().updateEstimate(globalIdx, wCompIdx)) > 1e-8 / v_w) ? problem_.variables().updateEstimate(globalIdx,wCompIdx) : 1e-8/v_w;
+ Scalar inc2 =(fabs(problem_.variables().updateEstimate(globalIdx, nCompIdx)) > 1e-8 / v_g) ? problem_.variables().updateEstimate(globalIdx,nCompIdx) : 1e-8 / v_g;
+ Scalar incp = 1e-2;
+
+
+ /**********************************
+ * c) Secant method for derivatives
+ **********************************/
+
+ // numerical derivative of fluid volume with respect to pressure
+ Scalar p_ = pressW + incp;
+ Scalar Z1 = m1 / (m1 + m2);
+ updFluidState.update(Z1,
+ p_, problem_.spatialParameters().porosity(globalPos, *ep), temperature_);
+ Scalar v_w_ = 1. / FluidSystem::phaseDensity(wPhaseIdx, temperature_, p_, updFluidState);
+ Scalar v_g_ = 1. / FluidSystem::phaseDensity(nPhaseIdx, temperature_, p_, updFluidState);
+ dV_dp = ((m1+m2) * (nuw1 * v_w_ + (1-nuw1) * v_g_) - volalt) /incp;
+
+ // numerical derivative of fluid volume with respect to mass of component 1
+ m1 += inc1;
+ Z1 = m1 / (m1 + m2);
+ updFluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *ep), temperature_);
+ Scalar satt = updFluidState.saturation(wPhaseIdx);
+ Scalar nuw = satt / v_w / (satt/v_w + (1-satt)/v_g);
+ dv_dC1 = ((m1+m2) * (nuw * v_w + (1-nuw) * v_g) - volalt) /inc1;
+ m1 -= inc1;
+
+ // numerical derivative of fluid volume with respect to mass of component 2
+ m2 += inc2;
+ Z1 = m1 / (m1 + m2);
+ updFluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *ep), temperature_);
+ satt = updFluidState.saturation(wPhaseIdx);
+ nuw = satt / v_w / (satt/v_w + (1-satt)/v_g);
+ dv_dC2 = ((m1+m2) * (nuw * v_w + (1-nuw) * v_g) - volalt)/ inc2;
+ m2 -= inc2;
+}
+
+
+}//end namespace Dumux
+#endif
diff --git a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
new file mode 100644
index 0000000000000000000000000000000000000000..102e3b08cffcab8bf1ddddd4c9f25feac78754fa
--- /dev/null
+++ b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
@@ -0,0 +1,1594 @@
+// $Id: fvpressure2p.hh 3357 2010-03-25 13:02:05Z lauser $
+/*****************************************************************************
+ * Copyright (C) 2007-2009 by Bernd Flemisch *
+ * Copyright (C) 2007-2009 by Jochen Fritz *
+ * Copyright (C) 2008-2009 by Markus Wolff *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_FVPRESSURE2P2CMULTIPHYSICS_HH
+#define DUMUX_FVPRESSURE2P2CMULTIPHYSICS_HH
+
+// dune environent:
+#include <dune/istl/bvector.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
+
+// dumux environment
+#include "dumux/common/pardiso.hh"
+#include "dumux/common/math.hh"
+#include <dumux/decoupled/2p2c/2p2cproperties.hh>
+#include <dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh>
+
+/**
+ * @file
+ * @brief Finite Volume Diffusion Model
+ * @author Bernd Flemisch, Jochen Fritz, Markus Wolff
+ */
+
+namespace Dumux
+{
+//! The finite volume model for the solution of the compositional pressure equation
+/*! \ingroup multiphysics
+ * Provides a Finite Volume implementation for the pressure equation of a gas-liquid
+ * system with two components. An IMPES-like method is used for the sequential
+ * solution of the problem. Capillary forces and diffusion are
+ * neglected. Isothermal conditions and local thermodynamic
+ * equilibrium are assumed. Gravity is included.
+ * \f[
+ c_{total}\frac{\partial p}{\partial t} + \sum_{\kappa} \frac{\partial v_{total}}{\partial C^{\kappa}} \nabla \cdot \left( \sum_{\alpha} X^{\kappa}_{\alpha} \varrho_{alpha} \bf{v}_{\alpha}\right)
+ = \sum_{\kappa} \frac{\partial v_{total}}{\partial C^{\kappa}} q^{\kappa},
+ * \f]
+ * where \f$\bf{v}_{\alpha} = - \lambda_{\alpha} \bf{K} \left(\nabla p_{\alpha} + \rho_{\alpha} \bf{g} \right) \f$.
+ * \f$ c_{total} \f$ represents the total compressibility, for constant porosity this yields \f$ - \frac{\partial V_{total}}{\partial p_{\alpha}} \f$,
+ * \f$p_{\alpha} \f$ denotes the phase pressure, \f$ \bf{K} \f$ the absolute permeability, \f$ \lambda_{\alpha} \f$ the phase mobility,
+ * \f$ \rho_{\alpha} \f$ the phase density and \f$ \bf{g} \f$ the gravity constant and \f$ C^{\kappa} \f$ the total Component concentration.
+ * See paper SPE 99619 or "Analysis of a Compositional Model for Fluid
+ * Flow in Porous Media" by Chen, Qin and Ewing for derivation.
+ *
+ * The partial derivatives of the actual fluid volume \f$ v_{total} \f$ are gained by using a secant method.
+ *
+ * The model domain is automatically divided
+ * in a single-phase and a two-phase domain. The full 2p2c model is only evaluated within the
+ * two-phase subdomain, whereas a single-phase transport model is computed in the rest of the
+ * domain.
+ *
+ * \tparam TypeTag The Type Tag
+ */
+template<class TypeTag> class FVPressure2P2CMultiPhysics
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) ReferenceElements;
+ typedef typename ReferenceElements::Container ReferenceElementContainer;
+ typedef typename ReferenceElements::ContainerFaces ReferenceElementFaceContainer;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+ typedef typename SpatialParameters::MaterialLaw MaterialLaw;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+ enum
+ {
+ pw = Indices::pressureW,
+ pn = Indices::pressureNW,
+ pglobal = Indices::pressureGlobal,
+ Sw = Indices::saturationW,
+ Sn = Indices::saturationNW,
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
+ wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx
+ };
+
+ // typedefs to abbreviate several dune classes...
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::Grid Grid;
+ typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ // convenience shortcuts for Vectors/Matrices
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix;
+
+ // the typenames used for the stiffness matrix and solution vector
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureCoefficientMatrix)) Matrix;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureRHSVector)) RHSVector;
+
+ //initializes the matrix to store the system of equations
+ void initializeMatrix();
+
+ //function which assembles the system of equations to be solved
+ void assemble(bool first);
+
+ //solves the system of equations to get the spatial distribution of the pressure
+ void solve();
+
+protected:
+ Problem& problem()
+ {
+ return problem_;
+ }
+ const Problem& problem() const
+ {
+ return problem_;
+ }
+
+public:
+ //the variables object is initialized, non-compositional before and compositional after first pressure calculation
+ void initialMaterialLaws(bool compositional);
+
+ //constitutive functions are initialized and stored in the variables object
+ void updateMaterialLaws();
+
+ //initialization routine to prepare first timestep
+ void initialize(bool solveTwice = false);
+
+ //pressure solution routine: update estimate for secants, assemble, solve.
+ void pressure(bool solveTwice = true)
+ {
+ //pre-transport to estimate update vector
+ Scalar dt_estimate = 0.;
+ Dune::dinfo << "secant guess"<< std::endl;
+ problem_.transportModel().update(-1, dt_estimate, problem_.variables().updateEstimate(), false);
+ //last argument false in update() makes shure that this is estimate and no "real" transport step
+ problem_.variables().updateEstimate() *= problem_.timeManager().timeStepSize();
+
+ assemble(false); Dune::dinfo << "pressure calculation"<< std::endl;
+ solve();
+
+ return;
+ }
+
+ void calculateVelocity()
+ {
+ // TODO: do this if we use a total velocity description. else, its just a vaste of time.
+ return;
+ }
+
+ //numerical volume derivatives wrt changes in mass, pressure
+ void volumeDerivatives(GlobalPosition globalPos, ElementPointer ep, Scalar& dv_dC1, Scalar& dv_dC2, Scalar& dV_dp);
+
+ /*! \name general methods for serialization, output */
+ //@{
+ // serialization methods
+ template<class Restarter>
+ void serialize(Restarter &res)
+ {
+ return;
+ }
+
+ template<class Restarter>
+ void deserialize(Restarter &res)
+ {
+ return;
+ }
+
+ //! \brief Write data files
+ /* \param name file name */
+ template<class MultiWriter>
+ void addOutputVtkFields(MultiWriter &writer)
+ {
+ problem().variables().addOutputVtkFields(writer);
+
+ // add multiphysics stuff
+ Dune::BlockVector<Dune::FieldVector<int,1> > *subdomainPtr = writer.template createField<int, 1> (dV_dp.size());
+ *subdomainPtr = problem_.variables().subdomain();
+ writer.addCellData(subdomainPtr, "subdomain");
+
+#if DUNE_MINIMAL_DEBUG_LEVEL <= 3
+ // add debug stuff
+ Dune::BlockVector<Dune::FieldVector<double,1> > *errorCorrPtr = writer.template createField<double, 1> (dV_dp.size());
+ *errorCorrPtr = errorCorrection;
+ writer.addCellData(errorCorrPtr, "Error Correction");
+ // add debug stuff
+ Dune::BlockVector<Dune::FieldVector<double,1> > *dV_dpPtr = writer.template createField<double, 1> (dV_dp.size());
+ *dV_dpPtr = dV_dp;
+ writer.addCellData(dV_dpPtr, "dV_dP");
+ Dune::BlockVector<Dune::FieldVector<double,1> > *dV_dC1Ptr = writer.template createField<double, 1> (dV_dp.size());
+ Dune::BlockVector<Dune::FieldVector<double,1> > *dV_dC2Ptr = writer.template createField<double, 1> (dV_dp.size());
+ *dV_dC1Ptr = dV_[0];
+ *dV_dC2Ptr = dV_[1];
+ writer.addCellData(dV_dC1Ptr, "dV_dC1");
+ writer.addCellData(dV_dC2Ptr, "dV_dC2");
+#endif
+
+ return;
+ }
+
+ //! Constructs a FVPressure2P2C object
+ /**
+ * \param problem a problem class object
+ */
+ FVPressure2P2CMultiPhysics(Problem& problem) :
+ problem_(problem), A_(problem.variables().gridSize(), problem.variables().gridSize(), (2 * dim + 1)
+ * problem.variables().gridSize(), Matrix::random), f_(problem.variables().gridSize()),
+ gravity(problem.gravity())
+ {
+ if (pressureType != pw && pressureType != pn && pressureType != pglobal)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
+ }
+ if (saturationType != Sw && saturationType != Sn)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
+ }
+
+ //rezise block vectors
+ dV_[wPhaseIdx].resize(problem_.variables().gridSize());
+ dV_[nPhaseIdx].resize(problem_.variables().gridSize());
+ dV_dp.resize(problem_.variables().gridSize());
+ problem_.variables().subdomain().resize(problem_.variables().gridSize());
+ nextSubdomain.resize(problem_.variables().gridSize());
+
+ errorCorrection.resize(problem_.variables().gridSize());
+
+ //prepare stiffness Matrix and Vectors
+ initializeMatrix();
+ }
+
+private:
+ Problem& problem_;
+ Matrix A_;
+ RHSVector f_;
+ std::string solverName_;
+ std::string preconditionerName_;
+
+ //vectors for partial derivatives
+ typename SolutionTypes::PhaseProperty dV_;
+ typename SolutionTypes::ScalarSolution dV_dp;
+
+ // subdomain map
+ Dune::BlockVector<Dune::FieldVector<int,1> > nextSubdomain; //! vector holding next subdomain
+
+ // debug
+ typename SolutionTypes::ScalarSolution errorCorrection; //debug output
+
+protected:
+ const Dune::FieldVector<Scalar, dimWorld>& gravity; //!< vector including the gravity constant
+ static const Scalar cFLFactor_ = GET_PROP_VALUE(TypeTag, PTAG(CFLFactor));
+ static const int pressureType = GET_PROP_VALUE(TypeTag, PTAG(PressureFormulation)); //!< gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
+ static const int saturationType = GET_PROP_VALUE(TypeTag, PTAG(SaturationFormulation)); //!< gives kind of saturation used (\f$ 0 = S_w \f$, \f$ 1 = S_n \f$)
+};
+
+//! initializes the matrix to store the system of equations
+template<class TypeTag>
+void FVPressure2P2CMultiPhysics<TypeTag>::initializeMatrix()
+{
+ // determine matrix row sizes
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // cell index
+ int globalIdxI = problem_.variables().index(*eIt);
+
+ // initialize row size
+ int rowSize = 1;
+
+ // run through all intersections with neighbors
+ IntersectionIterator isItEnd = problem_.gridView().template iend(*eIt);
+ for (IntersectionIterator isIt = problem_.gridView().template ibegin(*eIt); isIt != isItEnd; ++isIt)
+ if (isIt->neighbor())
+ rowSize++;
+ A_.setrowsize(globalIdxI, rowSize);
+ }
+ A_.endrowsizes();
+
+ // determine position of matrix entries
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // cell index
+ int globalIdxI = problem_.variables().index(*eIt);
+
+ // add diagonal index
+ A_.addindex(globalIdxI, globalIdxI);
+
+ // run through all intersections with neighbors
+ IntersectionIterator isItEnd = problem_.gridView().template iend(*eIt);
+ for (IntersectionIterator isIt = problem_.gridView().template ibegin(*eIt); isIt != isItEnd; ++isIt)
+ if (isIt->neighbor())
+ {
+ // access neighbor
+ ElementPointer outside = isIt->outside();
+ int globalIdxJ = problem_.variables().index(*outside);
+
+ // add off diagonal index
+ A_.addindex(globalIdxI, globalIdxJ);
+ }
+ }
+ A_.endindices();
+
+ return;
+}
+
+//! initializes the simulation run
+/*!
+ * Initializes the simulation to gain the initial pressure field.
+ *
+ * \param solveTwice flag to determine possible iterations of the initialization process
+ */
+template<class TypeTag>
+void FVPressure2P2CMultiPhysics<TypeTag>::initialize(bool solveTwice)
+{
+ // assign whole domain to most complex subdomain, => 2p
+ problem_.variables().subdomain() = 2;
+ // initialguess: set saturations, determine visco and mobility for initial pressure equation
+ // at this moment, the pressure is unknown. Hence, dont regard compositional effects.
+ initialMaterialLaws(false); Dune::dinfo << "first saturation guess"<<std::endl; //=J: initialGuess()
+
+ assemble(true); Dune::dinfo << "first pressure guess"<<std::endl;
+ solve();
+
+ // update the compositional variables (hence true)
+ initialMaterialLaws(true); Dune::dinfo << "first guess for mass fractions"<<std::endl;//= J: transportInitial()
+
+ // perform concentration update to determine secants
+ Scalar dt_estimate = 0.;
+ problem_.transportModel().update(0., dt_estimate, problem_.variables().updateEstimate(), false); Dune::dinfo << "secant guess"<< std::endl;
+ dt_estimate = std::min ( problem_.timeManager().timeStepSize(), dt_estimate);
+ problem_.variables().updateEstimate() *= dt_estimate;
+
+ // pressure calculation
+ assemble(false); Dune::dinfo << "second pressure guess"<<std::endl;
+ solve();
+ // update the compositional variables
+ initialMaterialLaws(true);
+
+ if (solveTwice)
+ {
+ Dune::BlockVector<Dune::FieldVector<Scalar, 1> > pressureOld(this->problem().variables().pressure());
+ Dune::BlockVector<Dune::FieldVector<Scalar, 1> > pressureDiff;
+ Scalar pressureNorm = 1.; //dummy initialization to perform at least 1 iteration
+ int numIter = 1;
+
+ while (pressureNorm > 1e-5 && numIter < 10)
+ {
+ Scalar dt_dummy=0.; // without this dummy, it will never converge!
+ // update for secants
+ problem_.transportModel().update(-1, dt_dummy, problem_.variables().updateEstimate(), false); Dune::dinfo << "secant guess"<< std::endl;
+ problem_.variables().updateEstimate() *= dt_estimate;
+ // pressure calculation
+ assemble(false); Dune::dinfo << "pressure guess number "<< numIter <<std::endl;
+ solve();
+ // update the compositional variables
+ initialMaterialLaws(true);
+
+ pressureDiff = pressureOld;
+ pressureDiff -= this->problem().variables().pressure();
+ pressureNorm = pressureDiff.infinity_norm();
+ pressureOld = this->problem().variables().pressure();
+ pressureNorm /= pressureOld.infinity_norm();
+
+ numIter++;
+ }
+ }
+ return;
+}
+
+//! function which assembles the system of equations to be solved
+/** for first == true, this function assembles the matrix and right hand side for
+ * the solution of the pressure field in the same way as in the class FVPressure2P.
+ * for first == false, the approach is changed to \f[-\frac{\partial V}{\partial p}
+ * \frac{\partial p}{\partial t}+\sum_{\kappa}\frac{\partial V}{\partial m^{\kappa}}\nabla\cdot
+ * \left(\sum_{\alpha}C_{\alpha}^{\kappa}\mathbf{v}_{\alpha}\right)
+ * =\sum_{\kappa}\frac{\partial V}{\partial m^{\kappa}}q^{\kappa} \f]. See Paper SPE 99619.
+ * This is done to account for the volume effects which appear when gas and liquid are dissolved in each other.
+ * \param first Flag if pressure field is unknown
+ */
+template<class TypeTag>
+void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
+{
+ // initialization: set matrix A_ to zero
+ A_ = 0;
+ f_ = 0;
+
+ // initialization: set the fluid volume derivatives to zero
+ Scalar dv_dC1(0.), dv_dC2(0.);
+ for (int i = 0; i < problem_.variables().gridSize(); i++)
+ {
+ dV_[wPhaseIdx][i] = 0; // dv_dC1 = dV/dm1
+ dV_[nPhaseIdx][i] = 0; // dv / dC2
+ dV_dp[i] = 0; // dv / dp
+ }
+
+ // determine maximum error to scale error-term
+ Scalar timestep_ = problem_.timeManager().timeStepSize();
+ Scalar maxErr = fabs(problem_.variables().volErr().infinity_norm()) / timestep_;
+
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // cell geometry type
+ Dune::GeometryType gt = eIt->geometry().type();
+
+ // number of Faces of current element
+ int numberOfFaces = ReferenceElementContainer::general(gt).size(1);
+
+
+ // get global coordinate of cell center
+ const GlobalPosition& globalPos = eIt->geometry().center();
+
+ // cell index
+ int globalIdxI = problem_.variables().index(*eIt);
+
+ // cell volume, assume linear map here
+ Scalar volume = eIt->geometry().volume();
+
+ Scalar densityWI = problem_.variables().densityWetting(globalIdxI);
+ Scalar densityNWI = problem_.variables().densityNonwetting(globalIdxI);
+
+ /****************implement sources and sinks************************/
+ Dune::FieldVector<Scalar,2> source(problem_.source(globalPos, *eIt));
+ if(first || problem_.variables().subdomain(globalIdxI) == 1)
+ {
+ source[wPhaseIdx] /= densityWI;
+ source[nPhaseIdx] /= densityNWI;
+ }
+ else
+ {
+ // derivatives of the fluid volume with respect to concentration of components, or pressure
+ if (dV_[0][globalIdxI] == 0)
+ volumeDerivatives(globalPos, *eIt, dV_[wPhaseIdx][globalIdxI][0], dV_[nPhaseIdx][globalIdxI][0], dV_dp[globalIdxI][0]);
+
+ source[wPhaseIdx] *= dV_[wPhaseIdx][globalIdxI]; // note: dV_[i][1] = dv_dC1 = dV/dm1
+ source[nPhaseIdx] *= dV_[nPhaseIdx][globalIdxI];
+ }
+ f_[globalIdxI] = volume * (source[wPhaseIdx] + source[nPhaseIdx]);
+ /********************************************************************/
+
+ // get absolute permeability
+ FieldMatrix permeabilityI(problem_.spatialParameters().intrinsicPermeability(globalPos, *eIt));
+
+ // get mobilities and fractional flow factors
+ Scalar lambdaWI = problem_.variables().mobilityWetting(globalIdxI);
+ Scalar lambdaNWI = problem_.variables().mobilityNonwetting(globalIdxI);
+ Scalar fractionalWI=0, fractionalNWI=0;
+ if (first)
+ {
+ fractionalWI = lambdaWI / (lambdaWI+ lambdaNWI);
+ fractionalNWI = lambdaNWI / (lambdaWI+ lambdaNWI);
+ }
+
+ Scalar pcI = problem_.variables().capillaryPressure(globalIdxI);
+
+ // prepare meatrix entries
+ Scalar entry=0.;
+ Scalar rightEntry=0.;
+
+ // iterate over all faces of the cell
+ IntersectionIterator isItEnd = problem_.gridView().template iend(*eIt);
+ for (IntersectionIterator isIt = problem_.gridView().template ibegin(*eIt); isIt != isItEnd; ++isIt)
+ {
+ // get geometry type of face
+ Dune::GeometryType faceGT = isIt->geometryInInside().type();
+
+ // center in face's reference element
+ const Dune::FieldVector<Scalar, dim - 1>& faceLocal = ReferenceElementFaceContainer::general(faceGT).position(0,0);
+
+ int isIndex = isIt->indexInInside();
+
+ // center of face inside volume reference element
+ const LocalPosition localPosFace(0);
+
+ // get normal vector
+ Dune::FieldVector<Scalar, dimWorld> unitOuterNormal = isIt->unitOuterNormal(faceLocal);
+
+ // get face volume
+ Scalar faceArea = isIt->geometry().volume();
+
+ /************* handle interior face *****************/
+ if (isIt->neighbor())
+ {
+ // access neighbor
+ ElementPointer neighborPointer = isIt->outside();
+ int globalIdxJ = problem_.variables().index(*neighborPointer);
+
+ // gemotry info of neighbor
+ Dune::GeometryType neighborGT = neighborPointer->geometry().type();
+ const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
+
+ // distance vector between barycenters
+ Dune::FieldVector<Scalar, dimWorld> distVec = globalPosNeighbor - globalPos;
+
+ // compute distance between cell centers
+ Scalar dist = distVec.two_norm();
+
+ Dune::FieldVector<Scalar, dimWorld> unitDistVec(distVec);
+ unitDistVec /= dist;
+
+ FieldMatrix permeabilityJ
+ = problem_.spatialParameters().intrinsicPermeability(globalPosNeighbor,
+ *neighborPointer);
+
+ // compute vectorized permeabilities
+ FieldMatrix meanPermeability(0);
+ Dumux::harmonicMeanMatrix(meanPermeability, permeabilityI, permeabilityJ);
+
+ Dune::FieldVector<Scalar, dim> permeability(0);
+ meanPermeability.mv(unitDistVec, permeability);
+
+ // get mobilities in neighbor
+ Scalar lambdaWJ = problem_.variables().mobilityWetting(globalIdxJ);
+ Scalar lambdaNWJ = problem_.variables().mobilityNonwetting(globalIdxJ);
+
+ // phase densities in cell in neighbor
+ Scalar densityWJ = problem_.variables().densityWetting(globalIdxJ);
+ Scalar densityNWJ = problem_.variables().densityNonwetting(globalIdxJ);
+
+ Scalar pcJ = problem_.variables().capillaryPressure(globalIdxJ);
+
+ Scalar rhoMeanW = 0.5 * (densityWI + densityWJ);
+ Scalar rhoMeanNW = 0.5 * (densityNWI + densityNWJ);
+
+ // reset potential gradients, density
+ Scalar potentialW = 0;
+ Scalar potentialNW = 0;
+ Scalar densityW = 0;
+ Scalar densityNW = 0;
+
+ if (first) // if we are at the very first iteration we can't calculate phase potentials
+ {
+ // get fractional flow factors in neigbor
+ Scalar fractionalWJ = lambdaWJ / (lambdaWJ+ lambdaNWJ);
+ Scalar fractionalNWJ = lambdaNWJ / (lambdaWJ+ lambdaNWJ);
+
+ // perform central weighting
+ Scalar lambda = (lambdaWI + lambdaWJ) * 0.5 + (lambdaNWI + lambdaNWJ) * 0.5;
+
+ entry = fabs(lambda*faceArea*(permeability*unitOuterNormal)/(dist));
+
+ Scalar factor = (fractionalWI + fractionalWJ) * (rhoMeanW) * 0.5 + (fractionalNWI + fractionalNWJ) * (rhoMeanNW) * 0.5;
+ rightEntry = factor * lambda * faceArea * (permeability * gravity) * (unitOuterNormal * unitDistVec);
+ }
+ else if (problem_.variables().subdomain(globalIdxI) == 1) //the current cell in the 1p domain
+ {
+ // 1p => no pC => only 1 pressure, potential
+ Scalar potential = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI]
+ - problem_.variables().pressure()[globalIdxJ]) / (dist*dist);
+ potentialW = potentialNW = potential;
+
+ potentialW += densityW * (unitDistVec * gravity);
+ potentialNW += densityNW * (unitDistVec * gravity);
+
+ double lambdaW, lambdaNW;
+
+ if (potentialW >= 0.)
+ {
+ lambdaW = lambdaWI;
+ densityW = densityWI;
+ }
+ else
+ {
+ lambdaW = lambdaWJ;
+ densityW = densityWJ;
+ }
+
+ if (potentialNW >= 0.)
+ {
+ lambdaNW = lambdaNWI;
+ densityNW = densityNWI;
+ }
+ else
+ {
+ lambdaNW = lambdaNWJ;
+ densityNW = densityNWJ;
+ }
+
+ entry = (lambdaW + lambdaNW) * faceArea * (permeability * unitOuterNormal) / (dist);
+ rightEntry = (densityW * lambdaW + densityNW * lambdaNW) * faceArea * (permeability * gravity) * (unitOuterNormal * unitDistVec);
+ }
+ else //current cell is in the 2p domain
+ {
+ Scalar graddv_dC1(0.), graddv_dC2(0.);
+
+ //check neighboring cell
+ if (problem_.variables().subdomain(globalIdxJ)== 2)
+ {
+ // determine volume derivatives
+ if (dV_[0][globalIdxJ] == 0)
+ volumeDerivatives(globalPosNeighbor, *neighborPointer, dV_[wPhaseIdx][globalIdxJ][0], dV_[nPhaseIdx][globalIdxJ][0], dV_dp[globalIdxJ][0]);
+ dv_dC1 = (dV_[wPhaseIdx][globalIdxI] + dV_[wPhaseIdx][globalIdxJ]) / 2; // dV/dm1= dV/dC^1
+ dv_dC2 = (dV_[nPhaseIdx][globalIdxI] + dV_[nPhaseIdx][globalIdxJ]) / 2;
+
+ graddv_dC1 = (dV_[wPhaseIdx][globalIdxJ] - dV_[wPhaseIdx][globalIdxI]) / dist;
+ graddv_dC2 = (dV_[nPhaseIdx][globalIdxJ] - dV_[nPhaseIdx][globalIdxI]) / dist;
+ }
+ else
+ {
+ dv_dC1 = dV_[wPhaseIdx][globalIdxI]; //only regard volume changes in 2p area => no averageing
+ dv_dC2 = dV_[nPhaseIdx][globalIdxI];
+
+ }
+
+// potentialW = problem_.variables().potentialWetting(globalIdxI, isIndex);
+// potentialNW = problem_.variables().potentialNonwetting(globalIdxI, isIndex);
+//
+// densityW = (potentialW > 0.) ? densityWI : densityWJ;
+// densityNW = (potentialNW > 0.) ? densityNWI : densityNWJ;
+//
+// densityW = (potentialW == 0.) ? rhoMeanW : densityW;
+// densityNW = (potentialNW == 0.) ? rhoMeanNW : densityNW;
+ //jochen: central weighting for gravity term
+ densityW = rhoMeanW; densityNW = rhoMeanNW;
+
+ switch (pressureType) //Markus: hab (unitOuterNormal * distVec)/dist hinzugefuegt
+ {
+ case pw:
+ {
+ potentialW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI]
+ - problem_.variables().pressure()[globalIdxJ]) / (dist*dist);
+ potentialNW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI]
+ - problem_.variables().pressure()[globalIdxJ] + pcI - pcJ) / (dist*dist);
+ break;
+ }
+ case pn:
+ {
+ potentialW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI]
+ - problem_.variables().pressure()[globalIdxJ] - pcI + pcJ) / (dist*dist);
+ potentialNW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI]
+ - problem_.variables().pressure()[globalIdxJ]) / (dist*dist);
+ break;
+ }
+ case pglobal:
+ {
+ DUNE_THROW(Dune::NotImplemented, "Global pressure not yet implemented for 2p2c");
+// potentialW = (problem_.variables().pressure()[globalIdxI]
+// - problem_.variables().pressure()[globalIdxJ] - fMeanNW * (pcI - pcJ)) / dist;
+// potentialNW = (problem_.variables().pressure()[globalIdxI]
+// - problem_.variables().pressure()[globalIdxJ] + fMeanW * (pcI - pcJ)) / dist;
+ break;
+ }
+ }
+ potentialW += densityW * (unitDistVec * gravity);
+ potentialNW += densityNW * (unitDistVec * gravity);
+
+ //store potentials for further calculations (velocity, saturation, ...)
+ problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialW;
+ problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
+
+ // initialize convenience shortcuts
+ Scalar lambdaW, lambdaN;
+ Scalar dV_w(0.), dV_n(0.); // dV_a = \sum_k \rho_a * dv/dC^k * X^k_a
+ Scalar gV_w(0.), gV_n(0.); // multipaper eq(3.3) zeile 3 analogon dV_w
+
+
+ //do the upwinding of the mobility & density depending on the phase potentials
+ if (potentialW >= 0.)
+ {
+ dV_w = (dv_dC1 * problem_.variables().wet_X1(globalIdxI) + dv_dC2 * (1. - problem_.variables().wet_X1(globalIdxI)));
+ dV_w *= densityWI;
+ lambdaW = problem_.variables().mobilityWetting(globalIdxI);
+
+ if (problem_.variables().subdomain(globalIdxJ)== 2)
+ {
+ gV_w = (graddv_dC1 * problem_.variables().wet_X1(globalIdxI) + graddv_dC2 * (1. - problem_.variables().wet_X1(globalIdxI)));
+ gV_w *= densityWI;
+ }
+ }
+ else
+ {
+ dV_w = (dv_dC1 * problem_.variables().wet_X1(globalIdxJ) + dv_dC2 * (1. - problem_.variables().wet_X1(globalIdxJ)));
+ dV_w *= densityWJ;
+ lambdaW = problem_.variables().mobilityWetting(globalIdxJ);
+ if (problem_.variables().subdomain(globalIdxJ)== 2)
+ {
+ gV_w = (graddv_dC1 * problem_.variables().wet_X1(globalIdxJ) + graddv_dC2 * (1. - problem_.variables().wet_X1(globalIdxJ)));
+ gV_w *= densityWJ;
+ }
+ }
+ if (potentialNW >= 0.)
+ {
+ dV_n = (dv_dC1 * problem_.variables().nonwet_X1(globalIdxI) + dv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxI)));
+ dV_n *= densityNWI;
+ lambdaN = problem_.variables().mobilityNonwetting(globalIdxI);
+ if (problem_.variables().subdomain(globalIdxJ)== 2)
+ {
+ gV_n = (graddv_dC1 * problem_.variables().nonwet_X1(globalIdxI) + graddv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxI)));
+ gV_n *= densityNWI;
+ }
+ }
+ else
+ {
+ dV_n = (dv_dC1 * problem_.variables().nonwet_X1(globalIdxJ) + dv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxJ)));
+ dV_n *= densityNWJ;
+ lambdaN = problem_.variables().mobilityNonwetting(globalIdxJ);
+ if (problem_.variables().subdomain(globalIdxJ)== 2)
+ {
+ gV_n = (graddv_dC1 * problem_.variables().nonwet_X1(globalIdxJ) + graddv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxJ)));
+ gV_n *= densityNWJ;
+ }
+ }
+
+ //calculate current matrix entry
+ entry = faceArea * (lambdaW * dV_w + lambdaN * dV_n);
+
+ //calculate right hand side
+ rightEntry = faceArea * (unitOuterNormal * unitDistVec) * (densityW * lambdaW * dV_w + densityNW * lambdaN * dV_n);
+ if (problem_.variables().subdomain(globalIdxJ)!= 1) // complex 2p subdomain
+ {
+ //subtract area integral
+ entry -= volume / numberOfFaces * (lambdaW * gV_w + lambdaN * gV_n);
+ rightEntry -= volume / numberOfFaces * (densityW * lambdaW * gV_w + densityNW * lambdaN * gV_n);
+ }
+ entry *= fabs((permeability*unitOuterNormal)/(dist)); // TODO: markus nimmt hier statt fabs() ein unitDistVec * unitDistVec
+ rightEntry *= (permeability * gravity);
+
+ // TODO: implement a proper capillary pressure
+// switch (pressureType)
+// {
+// case pw:
+// {
+// // calculate capillary pressure gradient
+// Dune::FieldVector<Scalar, dim> pCGradient = unitDistVec;
+// pCGradient *= (pcI - pcJ) / dist;
+//
+// //add capillary pressure term to right hand side
+// rightEntry += 0.5 * (lambdaNWI + lambdaNWJ) * (permeability * pCGradient) * faceArea;
+// break;
+// }
+// case pn:
+// {
+// // calculate capillary pressure gradient
+// Dune::FieldVector<Scalar, dim> pCGradient = unitDistVec;
+// pCGradient *= (pcI - pcJ) / dist;
+//
+// //add capillary pressure term to right hand side
+// rightEntry -= 0.5 * (lambdaWI + lambdaWJ) * (permeability * pCGradient) * faceArea;
+// break;
+// }
+// }
+ } // end !first
+
+ //set right hand side
+ f_[globalIdxI] -= rightEntry;
+ // set diagonal entry
+ A_[globalIdxI][globalIdxI] += entry;
+ // set off-diagonal entry
+ A_[globalIdxI][globalIdxJ] = -entry;
+ } // end neighbor
+ /****************************************************/
+
+
+ /************* boundary face ************************/
+ else
+ {
+ // get volume derivatives inside the cell
+ dv_dC1 = dV_[wPhaseIdx][globalIdxI];
+ dv_dC2 = dV_[nPhaseIdx][globalIdxI];
+
+ // center of face in global coordinates
+ const GlobalPosition& globalPosFace = isIt->geometry().global(faceLocal);
+
+ // geometrical information
+ Dune::FieldVector<Scalar, dimWorld> distVec(globalPosFace - globalPos);
+ Scalar dist = distVec.two_norm();
+ Dune::FieldVector<Scalar, dimWorld> unitDistVec(distVec);
+ unitDistVec /= dist;
+
+ //get boundary condition for boundary face center
+ BoundaryConditions::Flags bctype = problem_.bcTypePress(globalPosFace, *isIt);
+// BoundaryConditions::Flags bcTypeSat = problem_.bctypeSat(globalPosFace, *isIt);
+
+ /********** Dirichlet Boundary *************/
+ if (bctype == BoundaryConditions::dirichlet)
+ {
+ //permeability vector at boundary
+ Dune::FieldVector<Scalar, dim> permeability(0);
+ permeabilityI.mv(unitDistVec, permeability);
+
+ // create a fluid state for the boundary
+ FluidState BCfluidState;
+
+ Scalar temperatureBC = problem_.temperature(globalPosFace, *eIt);
+
+ //get dirichlet pressure boundary condition
+ Scalar pressBound = problem_.dirichletPress(globalPosFace, *isIt);
+
+ Scalar pressBC = 0.;
+ Scalar pcBound = 0.;
+ if(pressureType==pw)
+ {
+ pressBC = pressBound;
+ }
+ else if(pressureType==pn)
+ {
+ //TODO: take pC from variables or from MaterialLaw?
+ // if the latter, one needs Sw
+ pcBound = problem_.variables().capillaryPressure(globalIdxI);
+ pressBC = pressBound - pcBound;
+ }
+
+ if (first)
+ {
+ Scalar lambda = lambdaWI+lambdaNWI;
+ A_[globalIdxI][globalIdxI] += lambda * faceArea * (permeability * unitOuterNormal) / (dist);
+ pressBC = problem_.dirichletPress(globalPosFace, *isIt);
+ f_[globalIdxI] += lambda * faceArea * pressBC * (permeability * unitOuterNormal) / (dist);
+ rightEntry = (fractionalWI * densityWI
+ + fractionalNWI * densityNWI)
+ * lambda * faceArea * (unitOuterNormal * unitDistVec) *(permeability*gravity);
+ f_[globalIdxI] -= rightEntry;
+ }
+ else //not first
+ {
+ //get boundary condition type for compositional transport
+ BoundaryConditions2p2c::Flags bcform = problem_.bcFormulation(globalPosFace, *isIt);
+ if (bcform == BoundaryConditions2p2c::saturation) // saturation given
+ {
+ Scalar satBound = problem_.dirichletTransport(globalPosFace, *isIt);
+ BCfluidState.satFlash(satBound, pressBC, problem_.spatialParameters().porosity(globalPos, *eIt), temperatureBC);
+ }
+ else if (bcform == Dumux::BoundaryConditions2p2c::concentration) // mass fraction given
+ {
+ Scalar Z1Bound = problem_.dirichletTransport(globalPosFace, *isIt);
+ BCfluidState.update(Z1Bound, pressBC, problem_.spatialParameters().porosity(globalPos, *eIt), temperatureBC);
+ }
+ else // nothing declared at boundary
+ {
+ Scalar satBound = problem_.variables().saturation()[globalIdxI];
+ BCfluidState.satFlash(satBound, pressBC, problem_.spatialParameters().porosity(globalPos, *eIt), temperatureBC);
+ Dune::dwarn << "no boundary saturation/concentration specified on boundary pos " << globalPosFace << std::endl;
+ }
+
+ // determine fluid properties at the boundary
+ Scalar lambdaWBound = 0.;
+ Scalar lambdaNWBound = 0.;
+
+ Scalar densityWBound =
+ FluidSystem::phaseDensity(wPhaseIdx, temperatureBC, pressBC, BCfluidState);
+ Scalar densityNWBound =
+ FluidSystem::phaseDensity(nPhaseIdx, temperatureBC, pressBC+pcBound, BCfluidState);
+ Scalar viscosityWBound =
+ FluidSystem::phaseViscosity(wPhaseIdx, temperatureBC, pressBC, BCfluidState);
+ Scalar viscosityNWBound =
+ FluidSystem::phaseViscosity(nPhaseIdx, temperatureBC, pressBC+pcBound, BCfluidState);
+
+ // mobility at the boundary
+ switch (GET_PROP_VALUE(TypeTag, PTAG(BoundaryMobility)))
+ {
+ case Indices::satDependent:
+ {
+ lambdaWBound = BCfluidState.saturation(wPhaseIdx)
+ / viscosityWBound;
+ lambdaNWBound = BCfluidState.saturation(nPhaseIdx)
+ / viscosityNWBound;
+ break;
+ }
+ case Indices::permDependent:
+ {
+ lambdaWBound = MaterialLaw::krw(
+ problem_.spatialParameters().materialLawParams(globalPos, *eIt), BCfluidState.saturation(wPhaseIdx))
+ / viscosityWBound;
+ lambdaNWBound = MaterialLaw::krn(
+ problem_.spatialParameters().materialLawParams(globalPos, *eIt), BCfluidState.saturation(wPhaseIdx))
+ / viscosityNWBound;
+ }
+ }
+
+ Scalar rhoMeanW = 0.5 * (densityWI + densityWBound);
+ Scalar rhoMeanNW = 0.5 * (densityNWI + densityNWBound);
+
+ Scalar potentialW = 0, potentialNW = 0;
+
+// potentialW = problem_.variables().potentialWetting(globalIdxI, isIndex);
+// potentialNW = problem_.variables().potentialNonwetting(globalIdxI, isIndex);
+//
+// // do potential upwinding according to last potGradient vs Jochen: central weighting
+// densityW = (potentialW > 0.) ? densityWI : densityWBound;
+// densityNW = (potentialNW > 0.) ? densityNWI : densityNWBound;
+//
+// densityW = (potentialW == 0.) ? rhoMeanW : densityW;
+// densityNW = (potentialNW == 0.) ? rhoMeanNW : densityNW;
+ Scalar densityW=rhoMeanW;
+ Scalar densityNW=rhoMeanNW;
+
+ //calculate potential gradient
+ switch (pressureType)
+ {
+ case pw:
+ {
+ potentialW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI] - pressBound) / (dist * dist);
+ potentialNW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI] + pcI - pressBound - pcBound)
+ / (dist * dist);
+ break;
+ }
+ case pn:
+ {
+ potentialW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI] - pcI - pressBound + pcBound)
+ / (dist * dist);
+ potentialNW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI] - pressBound) / (dist * dist);
+ break;
+ }
+ case pglobal:
+ {
+ Scalar fractionalWBound = lambdaWBound / (lambdaWBound + lambdaNWBound);
+ Scalar fractionalNWBound = lambdaNWBound / (lambdaWBound + lambdaNWBound);
+ Scalar fMeanW = 0.5 * (fractionalWI + fractionalWBound);
+ Scalar fMeanNW = 0.5 * (fractionalNWI + fractionalNWBound);
+
+ potentialW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI] - pressBound - fMeanNW * (pcI
+ - pcBound)) / (dist * dist);
+ potentialNW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI] - pressBound + fMeanW * (pcI
+ - pcBound)) / (dist * dist);
+ break;
+ }
+ }
+
+ potentialW += densityW * (unitDistVec * gravity);
+ potentialNW += densityNW * (unitDistVec * gravity);
+
+ //store potential gradients for further calculations
+ problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialW;
+ problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
+
+ //do the upwinding of the mobility depending on the phase potentials
+ Scalar lambdaW, lambdaNW;
+ Scalar dV_w, dV_n; // gV_a weglassen, da dV/dc am Rand ortsunabhängig angenommen -> am rand nicht bestimmbar -> nur Randintegral ohne Gebietsintegral
+
+ if(problem_.variables().subdomain(globalIdxI)==1) // easy 1p subdomain
+ {
+ if (potentialW >= 0.)
+ {
+ densityW = (potentialW == 0) ? rhoMeanW : densityWI;
+ lambdaW = (potentialW == 0) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaWI;
+ }
+ else
+ {
+ densityW = densityWBound;
+ lambdaW = lambdaWBound;
+ }
+ if (potentialNW >= 0.)
+ {
+ densityNW = (potentialNW == 0) ? rhoMeanNW : densityNWI;
+ lambdaNW = (potentialNW == 0) ? 0.5 * (lambdaNWI + lambdaNWBound) : lambdaNWI;
+ }
+ else
+ {
+ densityNW = densityNWBound;
+ lambdaNW = lambdaNWBound;
+ }
+
+ //calculate current matrix entry
+ entry = (lambdaW + lambdaNW) * ((permeability * unitDistVec) / dist) * faceArea
+ * (unitOuterNormal * unitDistVec);
+
+ //calculate right hand side
+ rightEntry = (lambdaW * densityW + lambdaNW * densityNW) * (permeability * gravity)
+ * faceArea ;
+ }
+ else // 2p subdomain
+ {
+ if (potentialW >= 0.)
+ {
+ densityW = (potentialW == 0) ? rhoMeanW : densityWI;
+ dV_w = (dv_dC1 * problem_.variables().wet_X1(globalIdxI)
+ + dv_dC2 * (1. - problem_.variables().wet_X1(globalIdxI)));
+ dV_w *= densityW;
+ lambdaW = (potentialW == 0) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaWI;
+ }
+ else
+ {
+ densityW = densityWBound;
+ dV_w = (dv_dC1 * BCfluidState.massFrac(wPhaseIdx, wCompIdx) + dv_dC2 * BCfluidState.massFrac(wPhaseIdx, nCompIdx));
+ dV_w *= densityW;
+ lambdaW = lambdaWBound;
+ }
+ if (potentialNW >= 0.)
+ {
+ densityNW = (potentialNW == 0) ? rhoMeanNW : densityNWI;
+ dV_n = (dv_dC1 * problem_.variables().nonwet_X1(globalIdxI)
+ + dv_dC2 * (1. - problem_.variables().nonwet_X1(globalIdxI)));
+ dV_n *= densityNW;
+ lambdaNW = (potentialNW == 0) ? 0.5 * (lambdaNWI + lambdaNWBound) : lambdaNWI;
+ }
+ else
+ {
+ densityNW = densityNWBound;
+ dV_n = (dv_dC1 * BCfluidState.massFrac(nPhaseIdx, wCompIdx) + dv_dC2 * BCfluidState.massFrac(nPhaseIdx, nCompIdx));
+ dV_n *= densityNW;
+ lambdaNW = lambdaNWBound;
+ }
+
+ //calculate current matrix entry
+ entry = (lambdaW * dV_w + lambdaNW * dV_n) * ((permeability * unitDistVec) / dist) * faceArea
+ * (unitOuterNormal * unitDistVec);
+
+ //calculate right hand side
+ rightEntry = (lambdaW * densityW * dV_w + lambdaNW * densityNW * dV_n) * (permeability * gravity)
+ * faceArea ;
+
+ // switch (pressureType)
+ // {
+ // case pw:
+ // {
+ // // calculate capillary pressure gradient
+ // Dune::FieldVector<Scalar, dim> pCGradient = unitDistVec;
+ // pCGradient *= (pcI - pcBound) / dist;
+ //
+ // //add capillary pressure term to right hand side
+ // rightEntry += 0.5 * (lambdaNWI + lambdaNWBound) * (permeability * pCGradient) * faceArea;
+ // break;
+ // }
+ // case pn:
+ // {
+ // // calculate capillary pressure gradient
+ // Dune::FieldVector<Scalar, dim> pCGradient = unitDistVec;
+ // pCGradient *= (pcI - pcBound) / dist;
+ //
+ // //add capillary pressure term to right hand side
+ // rightEntry -= 0.5 * (lambdaWI + lambdaWBound) * (permeability * pCGradient) * faceArea;
+ // break;
+ // }
+ // }
+ } //end 2p subdomain
+
+
+ // set diagonal entry and right hand side entry
+ A_[globalIdxI][globalIdxI] += entry;
+ f_[globalIdxI] += entry * pressBound;
+ f_[globalIdxI] -= rightEntry * (unitOuterNormal * unitDistVec);
+ } //end of if(first) ... else{...
+ } // end dirichlet
+
+ /**********************************
+ * set neumann boundary condition
+ **********************************/
+ else
+ {
+ Dune::FieldVector<Scalar,2> J = problem_.neumann(globalPosFace, *isIt);
+ if (first || problem_.variables().subdomain(globalIdxI)==1)
+ {
+ J[wPhaseIdx] /= densityWI;
+ J[nPhaseIdx] /= densityNWI;
+ }
+ else
+ {
+ J[wPhaseIdx] *= dv_dC1;
+ J[nPhaseIdx] *= dv_dC2;
+ }
+
+ f_[globalIdxI] -= (J[wPhaseIdx] + J[nPhaseIdx]) * faceArea;
+
+ //Assumes that the phases flow in the same direction at the neumann boundary, which is the direction of the total flux!!!
+ //needed to determine the upwind direction in the saturation equation
+ problem_.variables().potentialWetting(globalIdxI, isIndex) = J[wPhaseIdx];
+ problem_.variables().potentialNonwetting(globalIdxI, isIndex) = J[nPhaseIdx];
+ }
+ /*************************************************/
+ }
+ } // end all intersections
+
+ // compressibility term
+ if (!first && timestep_ != 0.)
+ {
+ if (dV_dp[globalIdxI] == 0.)
+ {
+ assert(problem_.variables().subdomain(globalIdxI)==1);
+
+ // numerical derivative of fluid volume with respect to pressure
+ Scalar p_ = problem_.variables().pressure()[globalIdxI] + 1e-2;
+ Scalar Z1 = problem_.variables().totalConcentration(globalIdxI, wCompIdx)
+ / (problem_.variables().totalConcentration(globalIdxI, wCompIdx)
+ + problem_.variables().totalConcentration(globalIdxI, nCompIdx));
+ PseudoOnePTwoCFluidState<TypeTag> pseudoFluidState;
+ pseudoFluidState.update(Z1, p_, problem_.variables().saturation(globalIdxI));
+ if(problem_.variables().saturation(globalIdxI)==1) //only w-phase
+ {
+ Scalar v_w_ = 1. / FluidSystem::phaseDensity(wPhaseIdx,
+ problem_.temperature(globalPos, *eIt),
+ p_, pseudoFluidState);
+ dV_dp[globalIdxI] = (problem_.variables().totalConcentration(globalIdxI, wCompIdx)
+ + problem_.variables().totalConcentration(globalIdxI, nCompIdx))
+ * ( v_w_ - 1./densityWI) / 1e-2;
+ }
+ if(problem_.variables().saturation(globalIdxI)==0.) //only nw-phase
+ {
+ Scalar v_n_ = 1. / FluidSystem::phaseDensity(nPhaseIdx,
+ problem_.temperature(globalPos, *eIt),
+ p_, pseudoFluidState);
+ dV_dp[globalIdxI] = (problem_.variables().totalConcentration(globalIdxI, wCompIdx)
+ + problem_.variables().totalConcentration(globalIdxI, nCompIdx))
+ * ( v_n_ - 1./densityNWI) / 1e-2;
+ }
+ } //end calculation of 1p compress_term
+ Scalar compress_term = dV_dp[globalIdxI] / timestep_;
+
+ A_[globalIdxI][globalIdxI] -= compress_term*volume;
+ f_[globalIdxI] -= problem_.variables().pressure()[globalIdxI] * compress_term * volume;
+
+ if (isnan(compress_term) || isinf(compress_term))
+ DUNE_THROW(Dune::MathError, "Compressibility term leads to NAN matrix entry at index " << globalIdxI);
+
+ if(!GET_PROP_VALUE(TypeTag, PTAG(EnableCompressibility)))
+ DUNE_THROW(Dune::NotImplemented, "Compressibility is switched off???");
+ }
+
+ // error reduction routine: volumetric error is damped and inserted to right hand side
+ // if damping is not done, the solution method gets unstable!
+ problem_.variables().volErr()[globalIdxI] /= timestep_;
+ Scalar erri = fabs(problem_.variables().volErr()[globalIdxI]);
+ Scalar x_lo = 0.6;
+ Scalar x_mi = 0.9;
+ Scalar fac = 0.05;
+ Scalar lofac = 0.;
+ Scalar hifac = 0.;
+ hifac /= fac;
+
+ if ((erri*timestep_ > 5e-5) && (erri > x_lo * maxErr))
+ {
+ if (erri <= x_mi * maxErr)
+ f_[globalIdxI] += errorCorrection[globalIdxI] = fac* (1-x_mi*(lofac/fac-1)/(x_lo-x_mi) + (lofac/fac-1)/(x_lo-x_mi)*erri/maxErr)
+ * problem_.variables().volErr()[globalIdxI] * volume;
+ else
+ f_[globalIdxI] += errorCorrection[globalIdxI] = fac * (1 + x_mi - hifac*x_mi/(1-x_mi) + (hifac/(1-x_mi)-1)*erri/maxErr)
+ * problem_.variables().volErr()[globalIdxI] * volume;
+ }
+ } // end grid traversal
+ return;
+}
+
+//! solves the system of equations to get the spatial distribution of the pressure
+template<class TypeTag>
+void FVPressure2P2CMultiPhysics<TypeTag>::solve()
+{
+ typedef typename GET_PROP(TypeTag, PTAG(SolverParameters)) SolverParameters;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressurePreconditioner)) Preconditioner;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureSolver)) Solver;
+
+// printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
+// printvector(std::cout, f_, "right hand side", "row", 200, 1, 3);
+
+ Dune::MatrixAdapter<Matrix, RHSVector, RHSVector> op(A_); // make linear operator from A_
+ Dune::InverseOperatorResult result;
+ Scalar reduction = SolverParameters::reductionSolver;
+ int maxItSolver = SolverParameters::maxIterationNumberSolver;
+ int iterPreconditioner = SolverParameters::iterationNumberPreconditioner;
+ int verboseLevelSolver = SolverParameters::verboseLevelSolver;
+ Scalar relaxation = SolverParameters::relaxationPreconditioner;
+
+ if (verboseLevelSolver)
+ std::cout << "FVPressure2P: solve for pressure" << std::endl;
+
+ Preconditioner preconditioner(A_, iterPreconditioner, relaxation);
+ Solver solver(op, preconditioner, reduction, maxItSolver, verboseLevelSolver);
+ solver.apply(problem_.variables().pressure(), f_, result);
+
+ return;
+}
+
+/*!
+ * \brief initializes the fluid distribution and hereby the variables container
+ *
+ * This function equals the method initialguess() and transportInitial() in the old model.
+ * It differs from updateMaterialLaws because there are two possible initial conditions:
+ * saturations and concentration.
+ * \param compositional flag that determines if compositional effects are regarded, i.e.
+ * a reasonable pressure field is known.
+ */
+template<class TypeTag>
+void FVPressure2P2CMultiPhysics<TypeTag>::initialMaterialLaws(bool compositional)
+{
+ // initialize the fluid system
+ FluidState fluidState;
+
+ // iterate through leaf grid an evaluate c0 at cell center
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // get geometry type
+ Dune::GeometryType gt = eIt->geometry().type();
+
+ // get global coordinate of cell center
+ GlobalPosition globalPos = eIt->geometry().center();
+
+ // assign an Index for convenience
+ int globalIdx = problem_.variables().index(*eIt);
+
+ // get the temperature
+ Scalar temperature_ = problem_.temperature(globalPos, *eIt);
+
+ // initial conditions
+ Scalar pressW = 0;
+ Scalar pressNW = 0;
+ Scalar sat_0=0.;
+
+ BoundaryConditions2p2c::Flags ictype = problem_.initFormulation(globalPos, *eIt); // get type of initial condition
+
+ if(!compositional) //means that we do the first approximate guess without compositions
+ {
+ // phase pressures are unknown, so start with an exemplary
+ Scalar exemplaryPressure = problem_.referencePressure(globalPos, *eIt);
+ pressW = pressNW = problem_.variables().pressure()[globalIdx] = exemplaryPressure;
+
+ if (ictype == BoundaryConditions2p2c::saturation) // saturation initial condition
+ {
+ sat_0 = problem_.initSat(globalPos, *eIt);
+ fluidState.satFlash(sat_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ else if (ictype == BoundaryConditions2p2c::concentration) // concentration initial condition
+ {
+ Scalar Z1_0 = problem_.initConcentration(globalPos, *eIt);
+ fluidState.update(Z1_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ }
+ else if(compositional) //means we regard compositional effects since we know an estimate pressure field
+ {
+ //determine phase pressures from primary pressure variable
+ switch (pressureType)
+ {
+ case pw:
+ {
+ pressW = problem_.variables().pressure()[globalIdx];
+ pressNW = problem_.variables().pressure()[globalIdx] + problem_.variables().capillaryPressure(globalIdx);
+ break;
+ }
+ case pn:
+ {
+ pressW = problem_.variables().pressure()[globalIdx] - problem_.variables().capillaryPressure(globalIdx);
+ pressNW = problem_.variables().pressure()[globalIdx];
+ break;
+ }
+ }
+
+ if (ictype == Dumux::BoundaryConditions2p2c::saturation) // saturation initial condition
+ {
+ sat_0 = problem_.initSat(globalPos, *eIt);
+ fluidState.satFlash(sat_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ else if (ictype == Dumux::BoundaryConditions2p2c::concentration) // concentration initial condition
+ {
+ Scalar Z1_0 = problem_.initConcentration(globalPos, *eIt);
+ fluidState.update(Z1_0, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+ }
+ }
+
+ // initialize densities
+ problem_.variables().densityWetting(globalIdx) = FluidSystem::phaseDensity(wPhaseIdx, temperature_, pressW, fluidState);
+ problem_.variables().densityNonwetting(globalIdx) = FluidSystem::phaseDensity(nPhaseIdx, temperature_, pressNW, fluidState);
+
+ // initialize mass fractions
+ problem_.variables().wet_X1(globalIdx) = fluidState.massFrac(wPhaseIdx, wCompIdx);
+ problem_.variables().nonwet_X1(globalIdx) = fluidState.massFrac(nPhaseIdx, wCompIdx);
+
+
+ // initialize viscosities
+ problem_.variables().viscosityWetting(globalIdx) = FluidSystem::phaseViscosity(wPhaseIdx, temperature_, pressW, fluidState);
+ problem_.variables().viscosityNonwetting(globalIdx) = FluidSystem::phaseViscosity(nPhaseIdx, temperature_, pressNW, fluidState);
+
+ // initialize mobilities
+ problem_.variables().mobilityWetting(globalIdx) = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(globalPos, *eIt), fluidState.saturation(wPhaseIdx))
+ / problem_.variables().viscosityWetting(globalIdx);
+ problem_.variables().mobilityNonwetting(globalIdx) = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(globalPos, *eIt), fluidState.saturation(wPhaseIdx))
+ / problem_.variables().viscosityNonwetting(globalIdx);
+
+ // initialize cell concentration
+ problem_.variables().totalConcentration(globalIdx, wCompIdx) = fluidState.massConcentration(wCompIdx);
+ problem_.variables().totalConcentration(globalIdx, nCompIdx) = fluidState.massConcentration(nCompIdx);
+ problem_.variables().saturation()[globalIdx] = fluidState.saturation(wPhaseIdx);
+
+ }
+ return;
+}
+
+//! constitutive functions are updated once if new concentrations are calculated and stored in the variables container
+/*!
+ * In contrast to the standard sequential 2p2c model, this method also holds routines
+ * to adapt the subdomain. The subdomain indicates weather we are in 1p domain (value = 1)
+ * or in the two phase subdomain (value = 2).
+ */
+template<class TypeTag>
+void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
+{
+ // this method only completes the variables: = old postprocessupdate()
+
+ // instantiate standard 2p2c and pseudo1p fluid state objects
+ FluidState fluidState;
+ PseudoOnePTwoCFluidState<TypeTag> pseudoFluidState;
+
+ //get timestep for error term
+ Scalar dt = problem_.timeManager().timeStepSize();
+
+ // next subdomain map
+ if (problem_.timeManager().time() == 0.)
+ nextSubdomain = 2; // start with complicated sub in initialization
+ else
+ nextSubdomain = 1; // reduce complexity after first TS
+
+ // iterate through leaf grid an evaluate c0 at cell center
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // get geometry type
+ Dune::GeometryType gt = eIt->geometry().type();
+
+ // get global coordinate of cell center
+ GlobalPosition globalPos = eIt->geometry().center();
+
+ int globalIdx = problem_.variables().index(*eIt);
+
+ Scalar temperature_ = problem_.temperature(globalPos, *eIt);
+
+ // get the overall mass of component 1: Z1 = C^k / (C^1+C^2) [-]
+ Scalar Z1 = problem_.variables().totalConcentration(globalIdx, wCompIdx)
+ / (problem_.variables().totalConcentration(globalIdx, wCompIdx)
+ + problem_.variables().totalConcentration(globalIdx, nCompIdx));
+
+ if (problem_.variables().subdomain(globalIdx)==2) //=> 2p domain
+ {
+ //determine phase pressures from primary pressure variable
+ Scalar pressW(0.), pressNW(0.);
+ switch (pressureType)
+ {
+ case pw:
+ {
+ pressW = problem_.variables().pressure()[globalIdx];
+ Scalar oldSatW = problem_.variables().saturation(globalIdx);
+ pressNW = problem_.variables().pressure()[globalIdx]
+ + MaterialLaw::pC(problem_.spatialParameters().materialLawParams(globalPos, *eIt), oldSatW);
+ break;
+ }
+ case pn:
+ {
+ //todo: check this case for consistency throughout the model!
+ pressNW = problem_.variables().pressure()[globalIdx];
+ Scalar oldSatW = problem_.variables().saturation(globalIdx);
+ pressW = problem_.variables().pressure()[globalIdx]
+ - MaterialLaw::pC(problem_.spatialParameters().materialLawParams(globalPos, *eIt), oldSatW);
+ break;
+ }
+ }
+ //complete fluid state
+ fluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
+
+ /******** update variables in variableclass **********/
+ // initialize saturation
+ problem_.variables().saturation(globalIdx) = fluidState.saturation(wPhaseIdx);
+
+ // initialize pC todo: remove this dummy implementation
+ problem_.variables().capillaryPressure(globalIdx) = 0.0;
+
+ // initialize viscosities
+ problem_.variables().viscosityWetting(globalIdx)
+ = FluidSystem::phaseViscosity(wPhaseIdx, temperature_, pressW, fluidState);
+ problem_.variables().viscosityNonwetting(globalIdx)
+ = FluidSystem::phaseViscosity(nPhaseIdx, temperature_, pressNW, fluidState);
+
+ // initialize mobilities
+ problem_.variables().mobilityWetting(globalIdx) =
+ MaterialLaw::krw(problem_.spatialParameters().materialLawParams(globalPos, *eIt), fluidState.saturation(wPhaseIdx))
+ / problem_.variables().viscosityWetting(globalIdx);
+ problem_.variables().mobilityNonwetting(globalIdx) =
+ MaterialLaw::krn(problem_.spatialParameters().materialLawParams(globalPos, *eIt), fluidState.saturation(wPhaseIdx))
+ / problem_.variables().viscosityNonwetting(globalIdx);
+
+ // initialize mass fractions
+ problem_.variables().wet_X1(globalIdx) = fluidState.massFrac(wPhaseIdx, wCompIdx);
+ problem_.variables().nonwet_X1(globalIdx) = fluidState.massFrac(nPhaseIdx, wCompIdx);
+
+ // initialize densities
+ problem_.variables().densityWetting(globalIdx)
+ = FluidSystem::phaseDensity(wPhaseIdx, temperature_, pressW, fluidState);
+ problem_.variables().densityNonwetting(globalIdx)
+ = FluidSystem::phaseDensity(nPhaseIdx, temperature_, pressNW, fluidState);
+
+ problem_.spatialParameters().update(fluidState.saturation(wPhaseIdx), *eIt);
+
+ // determine volume mismatch between actual fluid volume and pore volume
+ Scalar sumConc = (problem_.variables().totalConcentration(globalIdx, wCompIdx)
+ + problem_.variables().totalConcentration(globalIdx, nCompIdx));
+ Scalar massw = sumConc * fluidState.phaseMassFraction(wPhaseIdx);
+ Scalar massn = sumConc * fluidState.phaseMassFraction(nPhaseIdx);
+ Scalar vol = massw / problem_.variables().densityWetting(globalIdx)
+ + massn / problem_.variables().densityNonwetting(globalIdx);
+ if (dt != 0)
+ {
+ problem_.variables().volErr()[globalIdx] = (vol - problem_.spatialParameters().porosity(globalPos, *eIt));
+
+ Scalar volErrI = problem_.variables().volErr(globalIdx);
+ if (std::isnan(volErrI))
+ {
+ DUNE_THROW(Dune::MathError, "Decoupled2p2c::postProcessUpdate:\n"
+ << "volErr[" << globalIdx << "] isnan: vol = " << vol
+ << ", massw = " << massw << ", rho_l = " << problem_.variables().densityWetting(globalIdx)
+ << ", massn = " << massn << ", rho_g = " << problem_.variables().densityNonwetting(globalIdx)
+ << ", poro = " << problem_.spatialParameters().porosity(globalPos, *eIt) << ", dt = " << dt);
+ }
+ }
+ else
+ {
+ problem_.variables().volErr()[globalIdx] = 0;
+ }
+
+ // check subdomain consistency
+ if (problem_.variables().saturation(globalIdx) != (1. || 0.)) // cell still 2p
+ {
+ // mark this element
+ nextSubdomain[globalIdx] = 2;
+
+ // mark neighbors
+ IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
+ for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt!=isItEnd; ++isIt)
+ {
+ if (isIt->neighbor())
+ {
+ int globalIdxJ = problem_.variables().index(*(isIt->outside()));
+ // mark neighbor Element
+ nextSubdomain[globalIdxJ] = 2;
+ }
+ }
+ }
+ // end subdomain check
+ }
+ else // simple
+ {
+ // Todo: update only variables of present phase!!
+ Scalar press1p = problem_.variables().pressure()[globalIdx];
+ pseudoFluidState.update(Z1, press1p, problem_.variables().saturation(globalIdx));
+
+ // initialize viscosities
+ problem_.variables().viscosityWetting(globalIdx)
+ = FluidSystem::phaseViscosity(wPhaseIdx, temperature_, press1p, pseudoFluidState);
+ problem_.variables().viscosityNonwetting(globalIdx)
+ = FluidSystem::phaseViscosity(nPhaseIdx, temperature_, press1p, pseudoFluidState);
+
+ // initialize mobilities
+ problem_.variables().mobilityWetting(globalIdx) =
+ MaterialLaw::krw(problem_.spatialParameters().materialLawParams(globalPos, *eIt), pseudoFluidState.saturation(wPhaseIdx))
+ / problem_.variables().viscosityWetting(globalIdx);
+ problem_.variables().mobilityNonwetting(globalIdx) =
+ MaterialLaw::krn(problem_.spatialParameters().materialLawParams(globalPos, *eIt), pseudoFluidState.saturation(wPhaseIdx))
+ / problem_.variables().viscosityNonwetting(globalIdx);
+
+ // initialize mass fractions
+ problem_.variables().wet_X1(globalIdx) = pseudoFluidState.massFrac(wPhaseIdx, wCompIdx);
+ problem_.variables().nonwet_X1(globalIdx) = pseudoFluidState.massFrac(nPhaseIdx, wCompIdx);
+
+ // initialize densities
+ problem_.variables().densityWetting(globalIdx)
+ = FluidSystem::phaseDensity(wPhaseIdx, temperature_, press1p, pseudoFluidState);
+ problem_.variables().densityNonwetting(globalIdx)
+ = FluidSystem::phaseDensity(nPhaseIdx, temperature_, press1p, pseudoFluidState);
+
+ // error term handling
+ Scalar sumConc = (problem_.variables().totalConcentration(globalIdx, wCompIdx)
+ + problem_.variables().totalConcentration(globalIdx, nCompIdx));
+ Scalar vol(0.);
+ if(problem_.variables().saturation(globalIdx) == 1.) //only w-phase
+ vol = sumConc / problem_.variables().densityWetting(globalIdx);
+ else //only nw-phase
+ vol = sumConc / problem_.variables().densityNonwetting(globalIdx);
+
+ if (dt != 0)
+ problem_.variables().volErr()[globalIdx] = (vol - problem_.spatialParameters().porosity(globalPos, *eIt));
+
+
+ }
+ }// end grid traversal
+
+ problem_.variables().subdomain() = nextSubdomain;
+
+ return;
+}
+
+//! partial derivatives of the volumes w.r.t. changes in total concentration and pressure
+/*!
+ * This method calculates the volume derivatives via a secant method, where the
+ * secants are gained in a pre-computational step via the transport equation and
+ * the last TS size.
+ * The partial derivatives w.r.t. mass are defined as
+ * \f$ \frac{\partial v}{\partial C^{\kappa}} = \frac{\partial V}{\partial m^{\kappa}}\f$
+ *
+ * \param globalPos The global position of the current element
+ * \param ep A pointer to the current element
+ * \param[out] dv_dC1 partial derivative of fluid volume w.r.t. mass of component 1 [m^3/kg]
+ * \param[out] dv_dC2 partial derivative of fluid volume w.r.t. mass of component 2 [m^3/kg]
+ * \param[out] dv_dp partial derivative of fluid volume w.r.t. pressure [1/Pa]
+ */
+template<class TypeTag>
+void FVPressure2P2CMultiPhysics<TypeTag>::volumeDerivatives(GlobalPosition globalPos, ElementPointer ep, Scalar& dv_dC1, Scalar& dv_dC2, Scalar& dV_dp)
+{
+ // cell index
+ int globalIdx = problem_.variables().index(*ep);
+
+ // get cell temperature
+ Scalar temperature_ = problem_.temperature(globalPos, *ep);
+
+ // initialize an Fluid state for the update
+ FluidState updFluidState;
+
+ /**********************************
+ * a) get necessary variables
+ **********************************/
+ //determine phase pressures from primary pressure variable
+ Scalar pressW=0.;
+ switch (pressureType)
+ {
+ case pw:
+ {
+ pressW = problem_.variables().pressure()[globalIdx];
+ break;
+ }
+ case pn:
+ {
+ pressW = problem_.variables().pressure()[globalIdx] - problem_.variables().capillaryPressure(globalIdx);
+ break;
+ }
+ }
+
+ Scalar v_w = 1. / problem_.variables().densityWetting(globalIdx);
+ Scalar v_g = 1. / problem_.variables().densityNonwetting(globalIdx);
+ Scalar sati = problem_.variables().saturation()[globalIdx];
+ // mass of components inside the cell
+ Scalar m1 = problem_.variables().totalConcentration(globalIdx, wCompIdx);
+ Scalar m2 = problem_.variables().totalConcentration(globalIdx, nCompIdx);
+ // mass fraction of wetting phase
+ Scalar nuw1 = sati/ v_w / (sati/v_w + (1-sati)/v_g);
+ // actual fluid volume
+ Scalar volalt = (m1+m2) * (nuw1 * v_w + (1-nuw1) * v_g);
+
+ /**********************************
+ * b) define increments
+ **********************************/
+ // increments for numerical derivatives
+ Scalar inc1 = (fabs(problem_.variables().updateEstimate(globalIdx, wCompIdx)) > 1e-8 / v_w) ? problem_.variables().updateEstimate(globalIdx,wCompIdx) : 1e-8/v_w;
+ Scalar inc2 =(fabs(problem_.variables().updateEstimate(globalIdx, nCompIdx)) > 1e-8 / v_g) ? problem_.variables().updateEstimate(globalIdx,nCompIdx) : 1e-8 / v_g;
+ Scalar incp = 1e-2;
+
+
+ /**********************************
+ * c) Secant method for derivatives
+ **********************************/
+
+ // numerical derivative of fluid volume with respect to pressure
+ Scalar p_ = pressW + incp;
+ Scalar Z1 = m1 / (m1 + m2);
+ updFluidState.update(Z1,
+ p_, problem_.spatialParameters().porosity(globalPos, *ep), temperature_);
+ Scalar v_w_ = 1. / FluidSystem::phaseDensity(wPhaseIdx, temperature_, p_, updFluidState);
+ Scalar v_g_ = 1. / FluidSystem::phaseDensity(nPhaseIdx, temperature_, p_, updFluidState);
+ dV_dp = ((m1+m2) * (nuw1 * v_w_ + (1-nuw1) * v_g_) - volalt) /incp;
+
+ // numerical derivative of fluid volume with respect to mass of component 1
+ m1 += inc1;
+ Z1 = m1 / (m1 + m2);
+ updFluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *ep), temperature_);
+ Scalar satt = updFluidState.saturation(wPhaseIdx);
+ Scalar nuw = satt / v_w / (satt/v_w + (1-satt)/v_g);
+ dv_dC1 = ((m1+m2) * (nuw * v_w + (1-nuw) * v_g) - volalt) /inc1;
+ m1 -= inc1;
+
+ // numerical derivative of fluid volume with respect to mass of component 2
+ m2 += inc2;
+ Z1 = m1 / (m1 + m2);
+ updFluidState.update(Z1, pressW, problem_.spatialParameters().porosity(globalPos, *ep), temperature_);
+ satt = updFluidState.saturation(wPhaseIdx);
+ nuw = satt / v_w / (satt/v_w + (1-satt)/v_g);
+ dv_dC2 = ((m1+m2) * (nuw * v_w + (1-nuw) * v_g) - volalt)/ inc2;
+ m2 -= inc2;
+}
+
+
+}//end namespace Dumux
+#endif
diff --git a/dumux/decoupled/2p2c/fvtransport2p2c.hh b/dumux/decoupled/2p2c/fvtransport2p2c.hh
new file mode 100644
index 0000000000000000000000000000000000000000..362dec0a40469a04a594e04637a578b6023f2391
--- /dev/null
+++ b/dumux/decoupled/2p2c/fvtransport2p2c.hh
@@ -0,0 +1,571 @@
+// $Id:
+/*****************************************************************************
+ * Copyright (C) 2010 by Markus Wolff, Benjamin Faigle *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_FVTRANSPORT2P2C_HH
+#define DUMUX_FVTRANSPORT2P2C_HH
+
+#include <dumux/decoupled/2p2c/2p2cproperties.hh>
+#include <dumux/common/math.hh>
+
+/**
+ * @file
+ * @brief Finite Volume discretization of the component transport equation
+ * @author Markus Wolff, Jochen Fritz, Benjamin Faigle
+ */
+
+namespace Dumux
+{
+//! Miscible Transport step in a Finite Volume discretization
+/*! \ingroup multiphase
+ * The finite volume model for the solution of the transport equation for compositional
+ * two-phase flow.
+ * \f[
+ \frac{\partial C^\kappa}{\partial t} = - \nabla \cdot \left( \sum_{\alpha} X^{\kappa}_{\alpha} \varrho_{alpha} \bf{v}_{\alpha}\right) + q^{\kappa},
+ * \f]
+ * where \f$ \bf{v}_{\alpha} = - \lambda_{\alpha} \bf{K} \left(\nabla p_{\alpha} + \rho_{\alpha} \bf{g} \right) \f$.
+ * \f$ p_{\alpha} \f$ denotes the phase pressure, \f$ \bf{K} \f$ the absolute permeability, \f$ \lambda_{\alpha} \f$ the phase mobility,
+ * \f$ \rho_{\alpha} \f$ the phase density and \f$ \bf{g} \f$ the gravity constant and \f$ C^{\kappa} \f$ the total Component concentration.
+ *
+ * \tparam TypeTag The Type Tag
+ */
+template<class TypeTag>
+class FVTransport2P2C
+{
+ 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(Problem)) Problem;
+ typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) ReferenceElements;
+ typedef typename ReferenceElements::Container ReferenceElementContainer;
+ typedef typename ReferenceElements::ContainerFaces ReferenceElementFaceContainer;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+ typedef typename SpatialParameters::MaterialLaw MaterialLaw;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TransportSolutionType)) TransportSolutionType;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+ enum
+ {
+ pw = Indices::pressureW,
+ pn = Indices::pressureNW,
+ pglobal = Indices::pressureGlobal,
+ vw = Indices::velocityW,
+ vn = Indices::velocityNW,
+ vt = Indices::velocityTotal,
+ Sw = Indices::saturationW,
+ Sn = Indices::saturationNW,
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
+ wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx
+ };
+
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::Grid Grid;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, 2> PhaseVector;
+
+public:
+ virtual void update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec, bool impes);
+
+ //! Set the initial values before the first pressure equation
+ /*!
+ * This method is called before first pressure equation is solved from Dumux::IMPET.
+ */
+ void initialize()
+ {};
+
+ void evalBoundary(GlobalPosition globalPosFace, IntersectionIterator& isIt,
+ ElementIterator& eIt, FluidState &BCfluidState, PhaseVector &pressure);
+
+ //! \brief Write data files
+ /* \param writer applied VTK-writer */
+ template<class MultiWriter>
+ void addOutputVtkFields(MultiWriter &writer)
+ {
+#if DUNE_MINIMAL_DEBUG_LEVEL <= 3
+ // add debug stuff
+ Dune::BlockVector<Dune::FieldVector<double,1> > *potentialDifferenceW = writer.template createField<double, 1> (problem_.variables().gridSize());
+ *potentialDifferenceW = potential_[wPhaseIdx];
+ writer.addCellData(potentialDifferenceW, "potentialDifferenceW pre/post pressure");
+ Dune::BlockVector<Dune::FieldVector<double,1> > *potentialDifferenceNW = writer.template createField<double, 1> (problem_.variables().gridSize());
+ *potentialDifferenceNW = potential_[nPhaseIdx];
+ writer.addCellData(potentialDifferenceNW, "potentialDifferenceNW pre/post pressure");
+#endif
+
+ return;
+ }
+
+
+ //! Constructs a FVTransport2P2C object
+ /*!
+ * Currently, the miscible transport scheme can not be applied with a global pressure / total velocity
+ * formulation.
+ *
+ * \param problem a problem class object
+ */
+ FVTransport2P2C(Problem& problem) :
+ problem_(problem), switchNormals(false)
+ {
+ const int velocityType = GET_PROP_VALUE(TypeTag, PTAG(VelocityFormulation));
+ if (GET_PROP_VALUE(TypeTag, PTAG(EnableCompressibility)) && velocityType == vt)
+ {
+ DUNE_THROW(Dune::NotImplemented,
+ "Total velocity - global pressure - model cannot be used with compressible fluids!");
+ }
+ const int saturationType = GET_PROP_VALUE(TypeTag, PTAG(SaturationFormulation));
+ if (saturationType != Sw && saturationType != Sn)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
+ }
+ if (pressureType != pw && pressureType != pn && pressureType != pglobal)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
+ }
+ if (velocityType != vw && velocityType != vn && velocityType != vt)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Velocity type not supported!");
+ }
+
+ potential_.resize(2);
+ potential_[0].resize(problem_.variables().gridSize());
+ potential_[1].resize(problem_.variables().gridSize());
+ }
+
+ ~FVTransport2P2C()
+ {
+ }
+
+private:
+ Problem& problem_;
+
+ TransportSolutionType potential_;
+
+ static const int pressureType = GET_PROP_VALUE(TypeTag, PTAG(PressureFormulation)); //!< gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
+ bool switchNormals;
+};
+
+//! \brief Calculate the update vector and determine timestep size
+/*!
+ * This method calculates the update vector \f$ u \f$ of the discretized equation
+ * \f[
+ C^{\kappa , new} = C^{\kappa , old} + u,
+ * \f]
+ * where \f$ u = \sum_{element faces} \boldsymbol{v}_{\alpha} * \varrho_{\alpha} * X^{\kappa}_{\alpha} * \boldsymbol{n} * A_{element face} \f$,
+ * \f$ \boldsymbol{n} \f$ is the face normal and \f$ A_{element face} \f$ is the face area.
+ *
+ * In addition to the \a update vector, the recommended time step size \a dt is calculated
+ * employing a CFL condition.
+ * This method = old concentrationUpdate()
+ *
+ * \param t Current simulation time [s]
+ * \param[out] dt Time step size [s]
+ * \param[out] updateVec Update vector, or update estimate for secants, resp. Here in [kg/m^3]
+ * \param impes Flag that determines if it is a real impes step or an update estimate for volume derivatives
+ */
+template<class TypeTag>
+void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec, bool impes = false)
+{
+ // initialize dt very large
+ dt = 1E100;
+
+ // set update vector to zero
+ updateVec[wCompIdx] = 0;
+ updateVec[nCompIdx] = 0;
+
+ // Cell which restricts time step size
+ int restrictingCell = -1;
+
+ // some phase properties
+ Dune::FieldVector<Scalar, dimWorld> gravity_ = problem_.gravity();
+
+ // compute update vector
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // cell index
+ int globalIdxI = problem_.variables().index(*eIt);
+
+ // cell geometry type
+ Dune::GeometryType gt = eIt->geometry().type();
+
+ // get position
+ GlobalPosition globalPos = eIt->geometry().center();
+
+ // cell volume, assume linear map here
+ Scalar volume = eIt->geometry().volume();
+
+ // get values of cell I
+ Scalar pressI = problem_.variables().pressure()[globalIdxI];
+
+ Dune::FieldMatrix<Scalar,dim,dim> K_I(problem_.spatialParameters().intrinsicPermeability(globalPos, *eIt));
+
+ Scalar SwmobI = std::max((problem_.variables().saturation(globalIdxI)
+ - problem_.spatialParameters().materialLawParams(globalPos, *eIt).Swr())
+ , 1e-2);
+ Scalar SnmobI = std::max((1. - problem_.variables().saturation(globalIdxI)
+ - problem_.spatialParameters().materialLawParams(globalPos, *eIt).Snr())
+ , 1e-2);
+
+
+ double Xw1_I = problem_.variables().wet_X1(globalIdxI);
+ double Xn1_I = problem_.variables().nonwet_X1(globalIdxI);
+
+ Scalar densityWI = problem_.variables().densityWetting(globalIdxI);
+ Scalar densityNWI = problem_.variables().densityNonwetting(globalIdxI);
+
+ // some variables for time step calculation
+ double sumfactorin = 0;
+ double sumfactorout = 0;
+
+ // run through all intersections with neighbors and boundary
+ IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
+ for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItEnd; ++isIt)
+ {
+ // get geometry type of face
+ Dune::GeometryType faceGT = isIt->geometryInInside().type();
+
+ // center in face's reference element
+ const Dune::FieldVector<Scalar, dim - 1>& faceLocal =
+ ReferenceElementFaceContainer::general(faceGT).position(0, 0);
+
+ // center of face inside volume reference element
+ const LocalPosition localPosFace(0);
+
+ Dune::FieldVector<Scalar, dimWorld> unitOuterNormal = isIt->unitOuterNormal(faceLocal);
+ if (switchNormals)
+ unitOuterNormal *= -1.0;
+
+ Scalar faceArea = isIt->geometry().volume();
+
+ // create vector for timestep and for update
+ Dune::FieldVector<Scalar, 2> factor (0.);
+ Dune::FieldVector<Scalar, 2> updFactor (0.);
+
+ // handle interior face
+ if (isIt->neighbor())
+ {
+ // access neighbor
+ ElementPointer neighborPointer = isIt->outside();
+ int globalIdxJ = problem_.variables().index(*neighborPointer);
+
+ // compute factor in neighbor
+ Dune::GeometryType neighborGT = neighborPointer->geometry().type();
+
+ // cell center in global coordinates
+ const GlobalPosition& globalPos = eIt->geometry().center();
+
+ // neighbor cell center in global coordinates
+ const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
+
+ // distance vector between barycenters
+ Dune::FieldVector<Scalar, dimWorld> distVec = globalPosNeighbor - globalPos;
+ // compute distance between cell centers
+ Scalar dist = distVec.two_norm();
+
+ Dune::FieldVector<Scalar, dimWorld> unitDistVec(distVec);
+ unitDistVec /= dist;
+
+ // get saturation and concentration value at neighbor cell center
+ double Xw1_J = problem_.variables().wet_X1(globalIdxJ);
+ double Xn1_J = problem_.variables().nonwet_X1(globalIdxJ);
+
+ // phase densities in neighbor
+ Scalar densityWJ = problem_.variables().densityWetting(globalIdxJ);
+ Scalar densityNWJ = problem_.variables().densityNonwetting(globalIdxJ);
+
+ // average phase densities with central weighting
+ double densityW_mean = (densityWI + densityWJ) * 0.5;
+ double densityNW_mean = (densityNWI + densityNWJ) * 0.5;
+
+ double pressJ = problem_.variables().pressure()[globalIdxJ];
+
+ // compute mean permeability
+ Dune::FieldMatrix<Scalar,dim,dim> meanK_(0.);
+ Dumux::harmonicMeanMatrix(meanK_,
+ K_I,
+ problem_.spatialParameters().intrinsicPermeability(globalPosNeighbor, *neighborPointer));
+ Dune::FieldVector<Scalar,dim> K(0);
+ meanK_.umv(unitDistVec,K);
+
+ // velocities
+ double potentialW = (unitOuterNormal * distVec) * (pressI - pressJ) / (dist * dist);
+ double potentialNW = potentialW + densityNW_mean * (unitDistVec * gravity_);
+ potentialW += densityW_mean * (unitDistVec * gravity_);
+
+ // upwind mobility
+ double lambdaW, lambdaN;
+ if (potentialW >= 0.)
+ lambdaW = problem_.variables().mobilityWetting(globalIdxI);
+ else
+ lambdaW = problem_.variables().mobilityWetting(globalIdxJ);
+
+ if (potentialNW >= 0.)
+ lambdaN = problem_.variables().mobilityNonwetting(globalIdxI);
+ else
+ lambdaN = problem_.variables().mobilityNonwetting(globalIdxJ);
+
+ double velocityW = lambdaW * ((K * unitOuterNormal) * (pressI - pressJ) / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean);
+ double velocityN = lambdaN * ((K * unitOuterNormal) * (pressI - pressJ) / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean);
+
+ // standardized velocity
+ double velocityJIw = std::max(-velocityW * faceArea / volume, 0.0);
+ double velocityIJw = std::max( velocityW * faceArea / volume, 0.0);
+ double velocityJIn = std::max(-velocityN * faceArea / volume, 0.0);
+ double velocityIJn = std::max( velocityN * faceArea / volume, 0.0);
+
+ // for timestep control
+ factor[0] = velocityJIw + velocityJIn;
+
+ double foutw = velocityIJw/SwmobI;
+ double foutn = velocityIJn/SnmobI;
+ if (std::isnan(foutw) || std::isinf(foutw) || foutw < 0) foutw = 0;
+ if (std::isnan(foutn) || std::isinf(foutn) || foutn < 0) foutn = 0;
+ factor[1] = foutw + foutn;
+
+ updFactor[wCompIdx] =
+ velocityJIw * Xw1_J * densityWJ
+ - velocityIJw * Xw1_I * densityWI
+ + velocityJIn * Xn1_J * densityNWJ
+ - velocityIJn * Xn1_I * densityNWI;
+ updFactor[nCompIdx] =
+ velocityJIw * (1. - Xw1_J) * densityWJ
+ - velocityIJw * (1. - Xw1_I) * densityWI
+ + velocityJIn * (1. - Xn1_J) * densityNWJ
+ - velocityIJn * (1. - Xn1_I) * densityNWI;
+ }
+
+ /******************************************
+ * Boundary Face
+ ******************************************/
+ if (isIt->boundary())
+ {
+ // center of face in global coordinates
+ const GlobalPosition& globalPosFace = isIt->geometry().global(faceLocal);
+
+ // cell center in global coordinates
+ const GlobalPosition& globalPos = eIt->geometry().center();
+
+ // distance vector between barycenters
+ Dune::FieldVector<Scalar, dimWorld> distVec = globalPosFace - globalPos;
+ // compute distance between cell centers
+ Scalar dist = distVec.two_norm();
+
+ Dune::FieldVector<Scalar, dimWorld> unitDistVec(distVec);
+ unitDistVec /= dist;
+
+ //instantiate a fluid state
+ FluidState BCfluidState;
+
+ //get boundary type
+ BoundaryConditions::Flags bcTypeTransport_ = problem_.bcTypeTransport(globalPosFace, *isIt);
+
+ /********** Dirichlet Boundary *************/
+ if (bcTypeTransport_ == BoundaryConditions::dirichlet)
+ {
+ //get dirichlet pressure boundary condition
+ Scalar pressBound = 0.;
+ //TODO: take pC from variables or from MaterialLaw?
+ // if the latter, one needs Sw
+ Scalar pcBound = problem_.variables().capillaryPressure(globalIdxI);
+ switch (pressureType)
+ {
+ case pw:
+ {
+ pressBound = problem_.dirichletPress(globalPosFace, *isIt);
+ break;
+ }
+ case pn:
+ {
+ pressBound = problem_.dirichletPress(globalPosFace, *isIt) - pcBound;
+ break;
+ }
+ }
+
+ // read boundary values
+ BoundaryConditions2p2c::Flags bctype = problem_.bcFormulation(globalPosFace, *isIt);
+ if (bctype == BoundaryConditions2p2c::saturation)
+ {
+ Scalar satBound = problem_.dirichletTransport(globalPosFace, *isIt);
+ BCfluidState.satFlash(satBound, pressBound,
+ problem_.spatialParameters().porosity(globalPos, *eIt),
+ problem_.temperature(globalPosFace, *eIt));
+
+ }
+ if (bctype == BoundaryConditions2p2c::concentration)
+ {
+ // saturation and hence pc and hence corresponding pressure unknown
+ Scalar Z1Bound = problem_.dirichletTransport(globalPosFace, *isIt);
+ BCfluidState.update(Z1Bound, pressBound,
+ problem_.spatialParameters().porosity(globalPos, *eIt),
+ problem_.temperature(globalPosFace, *eIt));
+ }
+ // determine fluid properties at the boundary
+ Scalar Xw1Bound = BCfluidState.massFrac(wPhaseIdx, wCompIdx);
+ Scalar Xn1Bound = BCfluidState.massFrac(nPhaseIdx, wCompIdx);
+ Scalar densityWBound = BCfluidState.density(wPhaseIdx);
+ Scalar densityNWBound = BCfluidState.density(nPhaseIdx);
+ Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx,
+ problem_.temperature(globalPosFace, *eIt),
+ pressBound, BCfluidState);
+ Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx,
+ problem_.temperature(globalPosFace, *eIt),
+ pressBound+pcBound, BCfluidState);
+
+ // average
+ double densityW_mean = (densityWI + densityWBound) / 2;
+ double densityNW_mean = (densityNWI + densityNWBound) / 2;
+
+ // velocities
+ double potentialW = (unitOuterNormal * distVec) * (pressI - pressBound) / (dist * dist);
+ double potentialNW = potentialW + densityNW_mean * (unitDistVec * gravity_);
+ potentialW += densityW_mean * (unitDistVec * gravity_);
+ potentialNW += densityNW_mean * (unitDistVec * gravity_);
+
+ // do upwinding for lambdas
+ double lambdaW, lambdaN;
+ if (potentialW >= 0.)
+ lambdaW = problem_.variables().mobilityWetting(globalIdxI);
+ else
+ {
+ if(GET_PROP_VALUE(TypeTag, PTAG(BoundaryMobility))==Indices::satDependent)
+ lambdaW = BCfluidState.saturation(wPhaseIdx) / viscosityWBound;
+ else
+ lambdaW = MaterialLaw::krw(
+ problem_.spatialParameters().materialLawParams(globalPos, *eIt), BCfluidState.saturation(wPhaseIdx))
+ / viscosityWBound;
+ }
+ if (potentialNW >= 0.)
+ lambdaN = problem_.variables().mobilityNonwetting(globalIdxI);
+ else
+ {
+ if(GET_PROP_VALUE(TypeTag, PTAG(BoundaryMobility))==Indices::satDependent)
+ lambdaN = BCfluidState.saturation(nPhaseIdx) / viscosityNWBound;
+ else
+ lambdaN = MaterialLaw::krn(
+ problem_.spatialParameters().materialLawParams(globalPos, *eIt), BCfluidState.saturation(nPhaseIdx))
+ / viscosityNWBound;
+ }
+
+ // prepare K
+ Dune::FieldVector<Scalar,dim> K(0);
+ K_I.umv(unitDistVec,K);
+
+ double velocityW = lambdaW * ((K * unitOuterNormal) * (pressI - pressBound)
+ / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean);
+ double velocityN = lambdaN * ((K * unitOuterNormal) * (pressI - pressBound)
+ / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean);
+
+ // standardized velocity
+ double velocityJIw = std::max(-velocityW * faceArea / volume, 0.0);
+ double velocityIJw = std::max( velocityW * faceArea / volume, 0.0);
+ double velocityJIn = std::max(-velocityN * faceArea / volume, 0.0);
+ double velocityIJn = std::max( velocityN* faceArea / volume, 0.0);
+
+ // for timestep control
+ factor[0] = velocityJIw + velocityJIn;
+
+ double foutw = velocityIJw/SwmobI;
+ double foutn = velocityIJn/SnmobI;
+ if (std::isnan(foutw) || std::isinf(foutw) || foutw < 0) foutw = 0;
+ if (std::isnan(foutn) || std::isinf(foutn) || foutn < 0) foutn = 0;
+ factor[1] = foutw + foutn;
+
+ updFactor[wCompIdx] =
+ + velocityJIw * Xw1Bound * densityWBound
+ - velocityIJw * Xw1_I * densityWI
+ + velocityJIn * Xn1Bound * densityNWBound
+ - velocityIJn * Xn1_I * densityNWI ;
+ updFactor[nCompIdx] =
+ velocityJIw * (1. - Xw1Bound) * densityWBound
+ - velocityIJw * (1. - Xw1_I) * densityWI
+ + velocityJIn * (1. - Xn1Bound) * densityNWBound
+ - velocityIJn * (1. - Xn1_I) * densityNWI ;
+ }//end dirichlet boundary
+
+ if (bcTypeTransport_ == BoundaryConditions::neumann)
+ {
+ Dune::FieldVector<Scalar,2> J = problem_.neumann(globalPosFace, *isIt);
+ updFactor[wCompIdx] = J[0] * faceArea / volume;
+ updFactor[nCompIdx] = J[1] * faceArea / volume;
+
+ // for timestep control
+ #define cflIgnoresNeumann
+ #ifdef cflIgnoresNeumann
+ factor[0] = 0;
+ factor[1] = 0;
+ #else
+ double inflow = updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW;
+ if (inflow>0)
+ {
+ factor[0] = updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW; // =factor in
+ factor[1] = -(updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI); // =factor out
+ }
+ else
+ {
+ factor[0] = -(updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW); // =factor in
+ factor[1] = updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI; // =factor out
+ }
+ #endif
+ }//end neumann boundary
+ }//end boundary
+ // add to update vector
+ updateVec[wCompIdx][globalIdxI] += updFactor[wCompIdx];
+ updateVec[nCompIdx][globalIdxI] += updFactor[nCompIdx];
+
+ // for time step calculation
+ sumfactorin += factor[0];
+ sumfactorout += factor[1];
+
+ }// end all intersections
+
+ /************************************
+ * Handle source term
+ ***********************************/
+ Dune::FieldVector<double,2> q = problem_.source(globalPos, *eIt);
+ updateVec[wCompIdx][globalIdxI] += q[wCompIdx];
+ updateVec[nCompIdx][globalIdxI] += q[nCompIdx];
+
+ // account for porosity
+ sumfactorin = std::max(sumfactorin,sumfactorout)
+ / problem_.spatialParameters().porosity(globalPos, *eIt);
+
+ if ( 1./sumfactorin < dt)
+ {
+ dt = 1./sumfactorin;
+ restrictingCell= globalIdxI;
+ }
+ } // end grid traversal
+ if(impes)
+ Dune::dinfo << "Timestep restricted by CellIdx " << restrictingCell << " leads to dt = "<<dt * GET_PROP_VALUE(TypeTag, PTAG(CFLFactor))<< std::endl;
+
+ return;
+ }
+}
+#endif
diff --git a/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh b/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4cd0b17b33dcd8013e0196d073119269ed02ff58
--- /dev/null
+++ b/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
@@ -0,0 +1,560 @@
+// $Id:
+/*****************************************************************************
+ * Copyright (C) 2010 by Markus Wolff, Benjamin Faigle *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_FVTRANSPORT2P2CMULTIPHYSICS_HH
+#define DUMUX_FVTRANSPORT2P2CMULTIPHYSICS_HH
+
+#include <dumux/decoupled/2p2c/2p2cproperties.hh>
+#include <dumux/common/math.hh>
+
+/**
+ * @file
+ * @brief Finite Volume discretization of the component transport equation
+ * @author Markus Wolff, Jochen Fritz, Benjamin Faigle
+ */
+
+namespace Dumux
+{
+//! Miscible Transport step in a Finite Volume discretization
+/*!
+ * \ingroup multiphysics
+ * The finite volume model for the solution of the transport equation for compositional
+ * two-phase flow.
+ * \f[
+ \frac{\partial C^\kappa}{\partial t} = - \nabla \cdot \left( \sum_{\alpha} X^{\kappa}_{\alpha} \varrho_{alpha} \bf{v}_{\alpha}\right) + q^{\kappa},
+ * \f]
+ * where \f$ \bf{v}_{\alpha} = - \lambda_{\alpha} \bf{K} \left(\nabla p_{\alpha} + \rho_{\alpha} \bf{g} \right) \f$.
+ * \f$ p_{\alpha} \f$ denotes the phase pressure, \f$ \bf{K} \f$ the absolute permeability, \f$ \lambda_{\alpha} \f$ the phase mobility,
+ * \f$ \rho_{\alpha} \f$ the phase density and \f$ \bf{g} \f$ the gravity constant and \f$ C^{\kappa} \f$ the total Component concentration.
+ *
+ * The model domain is automatically divided
+ * in a single-phase and a two-phase domain. The full 2p2c model is only evaluated within the
+ * two-phase subdomain, whereas a single-phase transport model is computed in the rest of the
+ * domain.
+ *
+ * \tparam TypeTag The Type Tag
+ */
+template<class TypeTag>
+class FVTransport2P2CMultiPhysics
+{
+ 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(Problem)) Problem;
+ typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) ReferenceElements;
+ typedef typename ReferenceElements::Container ReferenceElementContainer;
+ typedef typename ReferenceElements::ContainerFaces ReferenceElementFaceContainer;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+ typedef typename SpatialParameters::MaterialLaw MaterialLaw;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TransportSolutionType)) TransportSolutionType;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+ enum
+ {
+ pw = Indices::pressureW,
+ pn = Indices::pressureNW,
+ pglobal = Indices::pressureGlobal,
+ vw = Indices::velocityW,
+ vn = Indices::velocityNW,
+ vt = Indices::velocityTotal,
+ Sw = Indices::saturationW,
+ Sn = Indices::saturationNW,
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
+ wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx
+ };
+
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::Grid Grid;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ virtual void update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec, bool impes);
+
+ //! Set the initial values before the first pressure equation
+ /*!
+ * This method is called before first pressure equation is solved from Dumux::IMPET.
+ */
+ void initialize()
+ {};
+
+ //! \brief Write data files
+ /* \param writer applied VTK-writer */
+ template<class MultiWriter>
+ void addOutputVtkFields(MultiWriter &writer)
+ {
+ return;
+ }
+
+
+ //! Constructs a FVTransport2P2CMultiPhysics object
+ /**
+ * \param problem a problem class object
+ */
+
+ FVTransport2P2CMultiPhysics(Problem& problem) :
+ problem_(problem), switchNormals(false)
+ {
+ const int velocityType = GET_PROP_VALUE(TypeTag, PTAG(VelocityFormulation));
+ if (GET_PROP_VALUE(TypeTag, PTAG(EnableCompressibility)) && velocityType == vt)
+ {
+ DUNE_THROW(Dune::NotImplemented,
+ "Total velocity - global pressure - model cannot be used with compressible fluids!");
+ }
+ const int saturationType = GET_PROP_VALUE(TypeTag, PTAG(SaturationFormulation));
+ if (saturationType != Sw && saturationType != Sn)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
+ }
+ if (pressureType != pw && pressureType != pn && pressureType != pglobal)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
+ }
+ if (velocityType != vw && velocityType != vn && velocityType != vt)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Velocity type not supported!");
+ }
+ }
+
+ ~FVTransport2P2CMultiPhysics()
+ {
+ }
+
+private:
+ Problem& problem_;
+
+ static const int pressureType = GET_PROP_VALUE(TypeTag, PTAG(PressureFormulation)); //!< gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
+ bool switchNormals;
+};
+
+//! \brief Calculate the update vector and determine timestep size
+/*!
+ * This method calculates the update vector \f$ u \f$ of the discretized equation
+ * \f[
+ C^{\kappa , new} = C^{\kappa , old} + u,
+ * \f]
+ * where \f$ u = \sum_{element faces} \boldsymbol{v}_{\alpha} * \varrho_{\alpha} * X^{\kappa}_{\alpha} * \boldsymbol{n} * A_{element face} \f$,
+ * \f$ \boldsymbol{n} \f$ is the face normal and \f$ A_{element face} \f$ is the face area.
+ *
+ * In addition to the \a update vector, the recommended time step size \a dt is calculated
+ * employing a CFL condition.
+ * This method = old concentrationUpdate()
+ *
+ * \param t Current simulation time [s]
+ * \param[out] dt Time step size [s]
+ * \param[out] updateVec Update vector, or update estimate for secants, resp. Here in [kg/m^3]
+ * \param impes Flag that determines if it is a real impes step or an update estimate for volume derivatives
+ */
+template<class TypeTag>
+void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec, bool impet = false)
+{
+ // initialize dt very large
+ dt = 1E100;
+
+ // set update vector to zero
+ updateVec[wCompIdx] = 0;
+ updateVec[nCompIdx] = 0;
+
+ // Cell which restricts time step size
+ int restrictingCell = -1;
+
+ // some phase properties
+ Dune::FieldVector<Scalar, dimWorld> gravity_ = problem_.gravity();
+
+ // compute update vector
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // cell index
+ int globalIdxI = problem_.variables().index(*eIt);
+
+ if(impet || (problem_.variables().subdomain(globalIdxI)==2)) // estimate only necessary in subdomain
+ {
+ // cell geometry type
+ Dune::GeometryType gt = eIt->geometry().type();
+
+ // get position
+ GlobalPosition globalPos = eIt->geometry().center();
+
+ // cell volume, assume linear map here
+ Scalar volume = eIt->geometry().volume();
+
+ // get values of cell I
+ Scalar pressI = problem_.variables().pressure()[globalIdxI];
+ // Scalar pcI = problem_.variables().capillaryPressure(globalIdxI);
+ Dune::FieldMatrix<Scalar,dim,dim> K_I(problem_.spatialParameters().intrinsicPermeability(globalPos, *eIt));
+
+ Scalar SwmobI = std::max((problem_.variables().saturation(globalIdxI)
+ - problem_.spatialParameters().materialLawParams(globalPos, *eIt).Swr())
+ , 1e-2);
+ Scalar SnmobI = std::max((1. - problem_.variables().saturation(globalIdxI)
+ - problem_.spatialParameters().materialLawParams(globalPos, *eIt).Snr())
+ , 1e-2);
+
+ double Xw1_I = problem_.variables().wet_X1(globalIdxI);
+ double Xn1_I = problem_.variables().nonwet_X1(globalIdxI);
+
+ Scalar densityWI = problem_.variables().densityWetting(globalIdxI);
+ Scalar densityNWI = problem_.variables().densityNonwetting(globalIdxI);
+
+ // some variables for time step calculation
+ double sumfactorin = 0;
+ double sumfactorout = 0;
+
+ // run through all intersections with neighbors and boundary
+ IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
+ for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItEnd; ++isIt)
+ {
+ // get geometry type of face
+ Dune::GeometryType faceGT = isIt->geometryInInside().type();
+
+ // center in face's reference element
+ const Dune::FieldVector<Scalar, dim - 1>& faceLocal =
+ ReferenceElementFaceContainer::general(faceGT).position(0, 0);
+
+ // center of face inside volume reference element
+ const LocalPosition localPosFace(0);
+
+ Dune::FieldVector<Scalar, dimWorld> unitOuterNormal = isIt->unitOuterNormal(faceLocal);
+ if (switchNormals) // TODO: whats this??
+ unitOuterNormal *= -1.0;
+
+ Scalar faceArea = isIt->geometry().volume();
+
+ // create vector for timestep and for update
+ Dune::FieldVector<Scalar, 2> factor (0.);
+ Dune::FieldVector<Scalar, 2> updFactor (0.);
+
+ // handle interior face
+ if (isIt->neighbor())
+ {
+ // access neighbor
+ ElementPointer neighborPointer = isIt->outside();
+ int globalIdxJ = problem_.variables().index(*neighborPointer);
+
+ // compute factor in neighbor
+ Dune::GeometryType neighborGT = neighborPointer->geometry().type();
+
+ // cell center in global coordinates
+ const GlobalPosition& globalPos = eIt->geometry().center();
+
+ // neighbor cell center in global coordinates
+ const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
+
+ // distance vector between barycenters
+ Dune::FieldVector<Scalar, dimWorld> distVec = globalPosNeighbor - globalPos;
+ // compute distance between cell centers
+ Scalar dist = distVec.two_norm();
+
+ Dune::FieldVector<Scalar, dimWorld> unitDistVec(distVec);
+ unitDistVec /= dist;
+
+ // get saturation and concentration value at neighbor cell center
+ double Xw1_J = problem_.variables().wet_X1(globalIdxJ);
+ double Xn1_J = problem_.variables().nonwet_X1(globalIdxJ);
+
+ // phase densities in neighbor
+ Scalar densityWJ = problem_.variables().densityWetting(globalIdxJ);
+ Scalar densityNWJ = problem_.variables().densityNonwetting(globalIdxJ);
+
+ // average phase densities with central weighting
+ double densityW_mean = (densityWI + densityWJ) * 0.5;
+ double densityNW_mean = (densityNWI + densityNWJ) * 0.5;
+
+ double pressJ = problem_.variables().pressure()[globalIdxJ];
+
+ // compute mean permeability
+ Dune::FieldMatrix<Scalar,dim,dim> meanK_(0.);
+ Dumux::harmonicMeanMatrix(meanK_,
+ K_I,
+ problem_.spatialParameters().intrinsicPermeability(globalPosNeighbor, *neighborPointer));
+ Dune::FieldVector<Scalar,dim> K(0);
+ meanK_.umv(unitDistVec,K);
+
+ // velocities
+ double potentialW = (unitOuterNormal * distVec) * (pressI - pressJ) / (dist * dist);
+ double potentialNW = potentialW + densityNW_mean * (unitDistVec * gravity_);
+ potentialW += densityW_mean * (unitDistVec * gravity_);
+ // lambda = 2 * lambdaI * lambdaJ / (lambdaI + lambdaJ);
+
+ double lambdaW, lambdaN;
+
+ if (potentialW >= 0.)
+ lambdaW = problem_.variables().mobilityWetting(globalIdxI);
+ else
+ lambdaW = problem_.variables().mobilityWetting(globalIdxJ);
+
+ if (potentialNW >= 0.)
+ lambdaN = problem_.variables().mobilityNonwetting(globalIdxI);
+ else
+ lambdaN = problem_.variables().mobilityNonwetting(globalIdxJ);
+
+ double velocityW = lambdaW * ((K * unitOuterNormal) * (pressI - pressJ) / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean);
+ double velocityN = lambdaN * ((K * unitOuterNormal) * (pressI - pressJ) / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean);
+
+ // standardized velocity
+ double velocityJIw = std::max(-velocityW * faceArea / volume, 0.0);
+ double velocityIJw = std::max( velocityW * faceArea / volume, 0.0);
+ double velocityJIn = std::max(-velocityN * faceArea / volume, 0.0);
+ double velocityIJn = std::max( velocityN * faceArea / volume, 0.0);
+
+ // for timestep control
+ factor[0] = velocityJIw + velocityJIn;
+
+ double foutw = velocityIJw/SwmobI;
+ double foutn = velocityIJn/SnmobI;
+ if (std::isnan(foutw) || std::isinf(foutw) || foutw < 0) foutw = 0;
+ if (std::isnan(foutn) || std::isinf(foutn) || foutn < 0) foutn = 0;
+ factor[1] = foutw + foutn;
+
+ updFactor[wCompIdx] =
+ velocityJIw * Xw1_J * densityWJ
+ - velocityIJw * Xw1_I * densityWI
+ + velocityJIn * Xn1_J * densityNWJ
+ - velocityIJn * Xn1_I * densityNWI;
+ updFactor[nCompIdx] =
+ velocityJIw * (1. - Xw1_J) * densityWJ
+ - velocityIJw * (1. - Xw1_I) * densityWI
+ + velocityJIn * (1. - Xn1_J) * densityNWJ
+ - velocityIJn * (1. - Xn1_I) * densityNWI;
+ }
+
+ /******************************************
+ * Boundary Face
+ ******************************************/
+ if (isIt->boundary())
+ {
+ // center of face in global coordinates
+ GlobalPosition globalPosFace = isIt->geometry().global(faceLocal);
+
+ // cell center in global coordinates
+ GlobalPosition globalPos = eIt->geometry().center();
+
+ // distance vector between barycenters
+ Dune::FieldVector<Scalar, dimWorld> distVec = globalPosFace - globalPos;
+ // compute distance between cell centers
+ Scalar dist = distVec.two_norm();
+
+ Dune::FieldVector<Scalar, dimWorld> unitDistVec(distVec);
+ unitDistVec /= dist;
+
+ //instantiate a fluid state
+ FluidState fluidState;
+
+ //get boundary type
+ BoundaryConditions::Flags bcTypeTransport_ = problem_.bcTypeTransport(globalPosFace, *isIt);
+
+ if (bcTypeTransport_ == BoundaryConditions::dirichlet)
+ {
+ //get dirichlet pressure boundary condition
+ Scalar pressBound = 0.;
+ //TODO: take pC from variables or from MaterialLaw?
+ // if the latter, one needs Sw
+ Scalar pcBound = problem_.variables().capillaryPressure(globalIdxI);
+ switch (pressureType)
+ {
+ case pw:
+ {
+ pressBound = problem_.dirichletPress(globalPosFace, *isIt);
+ break;
+ }
+ case pn:
+ {
+ pressBound = problem_.dirichletPress(globalPosFace, *isIt) - pcBound;
+ break;
+ }
+ }
+
+ // read boundary values
+ BoundaryConditions2p2c::Flags bctype = problem_.bcFormulation(globalPosFace, *isIt);
+ if (bctype == BoundaryConditions2p2c::saturation)
+ {
+ Scalar satBound = problem_.dirichletTransport(globalPosFace, *isIt);
+ fluidState.satFlash(satBound, pressBound,
+ problem_.spatialParameters().porosity(globalPos, *eIt),
+ problem_.temperature(globalPosFace, *eIt));
+
+ }
+ if (bctype == BoundaryConditions2p2c::concentration)
+ {
+ Scalar Z1Bound = problem_.dirichletTransport(globalPosFace, *isIt);
+ fluidState.update(Z1Bound, pressBound,
+ problem_.spatialParameters().porosity(globalPos, *eIt),
+ problem_.temperature(globalPosFace, *eIt));
+ }
+ // determine fluid properties at the boundary
+ Scalar Xw1Bound = fluidState.massFrac(wPhaseIdx, wCompIdx);
+ Scalar Xn1Bound = fluidState.massFrac(nPhaseIdx, wCompIdx);
+ Scalar densityWBound = FluidSystem::phaseDensity(wPhaseIdx,
+ problem_.temperature(globalPosFace, *eIt),
+ pressBound, fluidState);
+ Scalar densityNWBound = FluidSystem::phaseDensity(nPhaseIdx,
+ problem_.temperature(globalPosFace, *eIt),
+ pressBound+pcBound, fluidState);
+ Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx,
+ problem_.temperature(globalPosFace, *eIt),
+ pressBound, fluidState);
+ Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx,
+ problem_.temperature(globalPosFace, *eIt),
+ pressBound+pcBound, fluidState);
+
+ // average
+ double densityW_mean = (densityWI + densityWBound) / 2;
+ double densityNW_mean = (densityNWI + densityNWBound) / 2;
+
+ // velocities
+ double potentialW = (unitOuterNormal * distVec) * (pressI - pressBound) / (dist * dist);
+ double potentialNW = potentialW + densityNW_mean * (unitDistVec * gravity_);
+ potentialW += densityW_mean * (unitDistVec * gravity_);
+
+ double lambdaW, lambdaN;
+
+ if (potentialW >= 0.)
+ lambdaW = problem_.variables().mobilityWetting(globalIdxI);
+ else
+ {
+ if(GET_PROP_VALUE(TypeTag, PTAG(BoundaryMobility))==Indices::satDependent)
+ lambdaW = fluidState.saturation(wPhaseIdx) / viscosityWBound;
+ else
+ lambdaW = MaterialLaw::krw(
+ problem_.spatialParameters().materialLawParams(globalPos, *eIt), fluidState.saturation(wPhaseIdx))
+ / viscosityWBound;
+ }
+ if (potentialNW >= 0.)
+ lambdaN = problem_.variables().mobilityNonwetting(globalIdxI);
+ else
+ {
+ if(GET_PROP_VALUE(TypeTag, PTAG(BoundaryMobility))==Indices::satDependent)
+ lambdaN = fluidState.saturation(nPhaseIdx) / viscosityNWBound;
+ else
+ lambdaN = MaterialLaw::krn(
+ problem_.spatialParameters().materialLawParams(globalPos, *eIt), fluidState.saturation(nPhaseIdx))
+ / viscosityNWBound;
+ }
+
+ // prepare K
+ Dune::FieldVector<Scalar,dim> K(0);
+ K_I.umv(unitDistVec,K);
+
+ double velocityW = lambdaW * ((K * unitOuterNormal) * (pressI - pressBound)
+ / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean);
+ double velocityN = lambdaN * ((K * unitOuterNormal) * (pressI - pressBound)
+ / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean);
+
+ // standardized velocity
+ double velocityJIw = std::max(-velocityW * faceArea / volume, 0.0);
+ double velocityIJw = std::max( velocityW * faceArea / volume, 0.0);
+ double velocityJIn = std::max(-velocityN * faceArea / volume, 0.0);
+ double velocityIJn = std::max( velocityN* faceArea / volume, 0.0);
+
+ // for timestep control
+ factor[0] = velocityJIw + velocityJIn;
+
+ double foutw = velocityIJw/SwmobI;
+ double foutn = velocityIJn/SnmobI;
+ if (std::isnan(foutw) || std::isinf(foutw) || foutw < 0) foutw = 0;
+ if (std::isnan(foutn) || std::isinf(foutn) || foutn < 0) foutn = 0;
+ factor[1] = foutw + foutn;
+
+ updFactor[wCompIdx] =
+ + velocityJIw * Xw1Bound * densityWBound
+ - velocityIJw * Xw1_I * densityWI
+ + velocityJIn * Xn1Bound * densityNWBound
+ - velocityIJn * Xn1_I * densityNWI ;
+ updFactor[nCompIdx] =
+ velocityJIw * (1. - Xw1Bound) * densityWBound
+ - velocityIJw * (1. - Xw1_I) * densityWI
+ + velocityJIn * (1. - Xn1Bound) * densityNWBound
+ - velocityIJn * (1. - Xn1_I) * densityNWI ;
+ }//end dirichlet boundary
+
+ if (bcTypeTransport_ == BoundaryConditions::neumann)
+ {
+ Dune::FieldVector<Scalar,2> J = problem_.neumann(globalPosFace, *isIt);
+ updFactor[wCompIdx] = J[0] * faceArea / volume;
+ updFactor[nCompIdx] = J[1] * faceArea / volume;
+
+ // for timestep control
+ #define cflIgnoresNeumann
+ #ifdef cflIgnoresNeumann
+ factor[0] = 0;
+ factor[1] = 0;
+ #else
+ double inflow = updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW;
+ if (inflow>0)
+ {
+ factor[0] = updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW; // =factor in
+ factor[1] = -(updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI); // =factor out
+ }
+ else
+ {
+ factor[0] = -(updFactor[wCompIdx] / densityW + updFactor[nCompIdx] / densityNW); // =factor in
+ factor[1] = updFactor[wCompIdx] / densityW /SwmobI + updFactor[nCompIdx] / densityNW / SnmobI; // =factor out
+ }
+ #endif
+ }//end neumann boundary
+ }//end boundary
+ // add to update vector
+ updateVec[wCompIdx][globalIdxI] += updFactor[wCompIdx];
+ updateVec[nCompIdx][globalIdxI] += updFactor[nCompIdx];
+
+ // for time step calculation
+ sumfactorin += factor[0];
+ sumfactorout += factor[1];
+
+ }// end all intersections
+
+ /************************************
+ * Handle source term
+ ***********************************/
+ Dune::FieldVector<double,2> q = problem_.source(globalPos, *eIt);
+ updateVec[wCompIdx][globalIdxI] += q[wCompIdx];
+ updateVec[nCompIdx][globalIdxI] += q[nCompIdx];
+
+ // account for porosity
+ sumfactorin = std::max(sumfactorin,sumfactorout)
+ / problem_.spatialParameters().porosity(globalPos, *eIt);
+
+ if ( 1./sumfactorin < dt)
+ {
+ dt = 1./sumfactorin;
+ restrictingCell= globalIdxI;
+ }
+ }
+ } // end grid traversal
+ if(impet)
+ Dune::dinfo << "Timestep restricted by CellIdx " << restrictingCell << " leads to dt = "<<dt * GET_PROP_VALUE(TypeTag, PTAG(CFLFactor))<< std::endl;
+
+ return;
+}
+
+}
+#endif
diff --git a/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh b/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9da2c2f82d2ff383333e56d0d74962bc0001e2d5
--- /dev/null
+++ b/dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh
@@ -0,0 +1,155 @@
+/*****************************************************************************
+ * Copyright (C) 2009 by Benjamin Faigle *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Calculates phase state for a single phase but two-component state.
+ */
+#ifndef DUMUX_PSEUDO1P2C_FLUID_STATE_HH
+#define DUMUX_PSEUDO1P2C_FLUID_STATE_HH
+
+#include <dumux/material/fluidstate.hh>
+#include <dumux/decoupled/2p2c/2p2cproperties.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup multiphysics
+ * \brief Container for compositional variables in a 1p2c situation
+ *
+ * This class represents a pseudo flash calculation when only single-phase situations
+ * occur. This is used in case of a multiphysics approach.
+ * \tparam TypeTag The property Type Tag
+ */
+template <class TypeTag>
+class PseudoOnePTwoCFluidState : public FluidState<typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)),
+ PseudoOnePTwoCFluidState<TypeTag> >
+{
+ typedef DecTwoPTwoCFluidState<TypeTag> ThisType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+
+public:
+ enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
+ numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),
+ wPhaseIdx = 0,
+ nPhaseIdx = 1,};
+
+public:
+ /*!
+ * \name flash calculation routines
+ * Routine to determine the phase composition after the transport step.
+ */
+ //@{
+ //! The simplest possible update routine for 1p2c "flash" calculations
+ /*!
+ * Routine goes as follows:
+ * - Check if we are in single phase condition
+ * - Assign total concentration to the present phase
+ *
+ * \param Z1 Feed mass fraction [-]
+ * \param press1p Pressure value for present phase [Pa]
+ * \param satW Saturation of the wetting phase [-]
+ */
+ void update(Scalar Z1,Scalar press1p, Scalar satW)
+ {
+ Sw_ = satW;
+ pressure1p_=press1p;
+
+ if (satW == 1.)
+ {
+ massfracX1_[wPhaseIdx] = Z1;
+ massfracX1_[nPhaseIdx] = 0.;
+ }
+ else if (satW == 0.)
+ {
+ massfracX1_[wPhaseIdx] = 0.;
+ massfracX1_[nPhaseIdx] = Z1;
+ }
+ else
+ Dune::dgrave << "Twophase conditions in simple-phase flash! Saturation is " << satW << std::endl;
+
+ return;
+ }
+ //@}
+ /*! \name Acess functions */
+ //@{
+ /*! \brief Returns the saturation of a phase.
+ * \param phaseIdx Index of the phase
+ */
+ Scalar saturation(int phaseIdx) const
+ {
+ if (phaseIdx == wPhaseIdx)
+ return Sw_;
+ else
+ return Scalar(1.0) - Sw_;
+ };
+
+ /*! \brief Returns the pressure of a fluid phase [Pa].
+ * \param phaseIdx Index of the phase
+ */
+ Scalar phasePressure(int phaseIdx) const
+ { return pressure1p_; }
+
+
+ /*!
+ * \brief Returns the mass fraction of a component in a phase.
+ * \param phaseIdx Index of the phase
+ * \param compIdx Index of the component
+ */
+ Scalar massFrac(int phaseIdx, int compIdx) const
+ {
+ if (compIdx == wPhaseIdx)
+ return massfracX1_[phaseIdx];
+ else
+ return 1.-massfracX1_[phaseIdx];
+
+ }
+
+ /*!
+ * \brief Returns the molar fraction of a component in a fluid phase.
+ * \param phaseIdx Index of the phase
+ * \param compIdx Index of the component
+ */
+ Scalar moleFrac(int phaseIdx, int compIdx) const
+ {
+ double moleFrac_(0.);
+ // as the mass fractions are calculated, these are used to determine the mole fractions
+ if (compIdx == wPhaseIdx) //only interested in wComp => X1
+ {
+ moleFrac_ = ( massfracX1_[phaseIdx] / FluidSystem::molarMass(compIdx) ); // = moles of compIdx
+ moleFrac_ /= ( massfracX1_[phaseIdx] / FluidSystem::molarMass(0)
+ + (1.-massfracX1_[phaseIdx]) / FluidSystem::molarMass(1) ); // /= total moles in phase
+ }
+ else // interested in nComp => 1-X1
+ {
+ moleFrac_ = ( (1.-massfracX1_[phaseIdx]) / FluidSystem::molarMass(compIdx) ); // = moles of compIdx
+ moleFrac_ /= ((1.- massfracX1_[phaseIdx] )/ FluidSystem::molarMass(0)
+ + massfracX1_[phaseIdx] / FluidSystem::molarMass(1) ); // /= total moles in phase
+ }
+ return moleFrac_;
+ }
+ //@}
+
+public:
+ Scalar massfracX1_[numPhases];
+ Scalar Sw_;
+ Scalar pressure1p_;
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/decoupled/2p2c/variableclass2p2c.hh b/dumux/decoupled/2p2c/variableclass2p2c.hh
new file mode 100644
index 0000000000000000000000000000000000000000..81498d80a7e21ddf8138b5d0b7386a5cf3c6bcfc
--- /dev/null
+++ b/dumux/decoupled/2p2c/variableclass2p2c.hh
@@ -0,0 +1,488 @@
+// $Id: variableclass2p.hh 3357 2010-03-25 13:02:05Z lauser $
+/*****************************************************************************
+ * Copyright (C) 2009 by Markus Wolff *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_VARIABLECLASS2P2C_HH
+#define DUMUX_VARIABLECLASS2P2C_HH
+
+#include <dumux/decoupled/common/variableclass.hh>
+#include <dumux/decoupled/2p2c/dec2p2cfluidstate.hh>
+
+/**
+ * @file
+ * @brief Class including the variables and data of discretized data of the constitutive relations
+ * @author Markus Wolff, Benjamin Faigle
+ */
+
+namespace Dumux
+{
+/*!
+ * \ingroup IMPEC
+ */
+//! Class including the variables and data of discretized data of the constitutive relations.
+/*! The variables of compositional two-phase flow, which are one pressure and one saturation are stored in this class.
+ *
+ * \tparam TypeTag The Type Tag
+ */
+template<class TypeTag>
+class VariableClass2P2C: public VariableClass<TypeTag>
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TransportSolutionType)) TransportSolutionType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ typedef VariableClass<TypeTag> ParentClass;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+ enum
+ {
+ pw = Indices::pressureW, pn = Indices::pressureNW, pglobal = Indices::pressureGlobal,
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+ };
+
+ static const int pressureType_ = GET_PROP_VALUE(TypeTag, PTAG(PressureFormulation));
+
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef typename GridView::IndexSet IndexSet;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+public:
+ typedef typename SolutionTypes::ScalarSolution ScalarSolutionType;//!<type for vector of scalars
+ typedef typename SolutionTypes::PhaseProperty PhasePropertyType;//!<type for vector of phase properties = FluidPropertyType
+ typedef typename SolutionTypes::FluidProperty FluidPropertyType;//!<type for vector of fluid properties = PhasePropertyType -> decoupledproperties.hh
+ typedef typename SolutionTypes::PhasePropertyElemFace PhasePropertyElemFaceType;//!<type for vector of vectors (of size 2 x dimension) of scalars
+ typedef typename SolutionTypes::DimVecElemFace DimVecElemFaceType;//!<type for vector of vectors (of size 2 x dimension) of vector (of size dimension) of scalars
+
+private:
+ const int codim_;
+
+ ScalarSolutionType saturation_;
+ PhasePropertyType mobility_; //store lambda for efficiency reasons
+ ScalarSolutionType capillaryPressure_;
+ DimVecElemFaceType velocitySecondPhase_; //necessary??
+
+ FluidPropertyType density_;
+ FluidPropertyType viscosity_;
+
+ ScalarSolutionType volErr_;
+
+ TransportSolutionType totalConcentration_;
+ PhasePropertyType massfrac_;
+
+ TransportSolutionType updateEstimate_;
+ Dune::BlockVector<Dune::FieldVector<int,1> > subdomain_;
+
+public:
+
+ //! Constructs a VariableClass object
+ /**
+ * @param gridView a DUNE gridview object corresponding to diffusion and transport equation
+ * @param initialSat initial value for the saturation (only necessary if only diffusion part is solved)
+ * @param initialVel initial value for the velocity (only necessary if only transport part is solved)
+ */
+
+ VariableClass2P2C(const GridView& gridView, Scalar& initialSat = *(new Scalar(1)),
+ Dune::FieldVector<Scalar, dim>& initialVel = *(new Dune::FieldVector<Scalar, dim>(0))) :
+ VariableClass<TypeTag> (gridView, initialVel), codim_(0)
+ {
+ initialize2p2cVariables(initialVel, initialSat);
+ }
+
+ //! Constructs a VariableClass object
+ /**
+ * @param gridView a DUNE gridview object corresponding to diffusion and transport equation
+ * @param codim codimension of the entity of which data has to be strored
+ * @param initialSat initial value for the saturation (only necessary if only diffusion part is solved)
+ * @param initialVel initial value for the velocity (only necessary if only transport part is solved)
+ */
+ VariableClass2P2C(const GridView& gridView, int codim, Scalar& initialSat = *(new Scalar(1)), Dune::FieldVector<
+ Scalar, dim>& initialVel = *(new Dune::FieldVector<Scalar, dim>(0))) :
+ VariableClass<TypeTag> (gridView, codim, initialVel), codim_(codim)
+ {
+ initialize2p2cVariables(initialVel, initialSat);
+ }
+
+ //! serialization methods -- same as Dumux::VariableClass2P
+ //@{
+ template<class Restarter>
+ void serialize(Restarter &res)
+ {
+ res.template serializeEntities<0> (*this, this->gridView());
+ }
+ template<class Restarter>
+ void deserialize(Restarter &res)
+ {
+ res.template deserializeEntities<0> (*this, this->gridView());
+ }
+
+ void serializeEntity(std::ostream &outstream, const Element &element)
+ {
+ Dune::dwarn << "here, rather total concentrations should be used!" << std::endl;
+ int globalIdx = this->elementMapper().map(element);
+ outstream << this->pressure()[globalIdx] << " " << saturation_[globalIdx];
+ }
+ void deserializeEntity(std::istream &instream, const Element &element)
+ {
+ int globalIdx = this->elementMapper().map(element);
+ instream >> this->pressure()[globalIdx] >> saturation_[globalIdx];
+ }
+ //@}
+
+
+private:
+ //! initializes the miscible 2p variables
+ /*! Internal method that prepares the vectors holding the 2p2c variables for the current
+ * simulation.
+ *
+ *\param initialVel Vector containing the initial velocity
+ *\param initialSat Initial value for the saturation
+ */
+ void initialize2p2cVariables(Dune::FieldVector<Scalar, dim>& initialVel, int initialSat)
+ {
+ //resize to grid size
+ int size_ = this->gridSize();
+ // a) global variables
+ velocitySecondPhase_.resize(size_);//depends on pressure
+ velocitySecondPhase_ = initialVel;
+ for (int i=0; i<2; i++) //for both phases
+ {
+ density_[i].resize(size_);//depends on pressure
+ viscosity_[i].resize(size_);//depends on pressure
+ };
+ // b) transport variables
+ saturation_.resize(size_);
+ saturation_ = initialSat;
+ capillaryPressure_.resize(size_);//depends on saturation
+ mobility_[0].resize(size_);//lambda is dependent on saturation! ->choose same size
+ mobility_[1].resize(size_);
+
+ // c) compositional stuff
+ totalConcentration_.resize(2); // resize solution vector to 2 primary transport variables
+ updateEstimate_.resize(2);
+ for (int i=0; i<2; i++) //for both phases
+ {
+ totalConcentration_[i].resize(size_); // resize vector to number of cells
+ updateEstimate_[i].resize(size_);
+ massfrac_[i].resize(size_);
+ }
+ volErr_.resize(size_);
+ }
+
+
+public:
+
+ //! Writes output into file
+ /*!
+ * Creates output for standard 2p2c models.
+ * \param writer Applied VTK-writer
+ */
+ template<class MultiWriter>
+ void addOutputVtkFields(MultiWriter &writer)
+ {
+ int size_ = this->gridSize();
+ if (codim_ == 0)
+ {
+ // pressure & saturation
+ ScalarSolutionType *pressure = writer.template createField<Scalar, 1> (size_);
+ ScalarSolutionType *saturation = writer.template createField<Scalar, 1> (size_);
+ ScalarSolutionType *pc = writer.template createField<Scalar, 1> (size_);
+ *pressure = this->pressure();
+ *saturation = saturation_;
+ *pc = capillaryPressure_;
+ writer.addCellData(pressure, "pressure");
+ writer.addCellData(saturation, "water saturation");
+ writer.addCellData(pc, "capillary pressure");
+
+ //total Concentration
+ ScalarSolutionType *totalConcentration1 = writer.template createField<Scalar, 1> (size_);
+ ScalarSolutionType *totalConcentration2 = writer.template createField<Scalar, 1> (size_);
+ *totalConcentration1 = totalConcentration_[0];
+ *totalConcentration2 = totalConcentration_[1];
+ writer.addCellData(totalConcentration1, "totalConcentration w-Comp");
+ writer.addCellData(totalConcentration2, "totalConcentration n-Comp");
+
+ // numerical stuff
+ ScalarSolutionType *volErr = writer.template createField<Scalar, 1> (size_);
+ *volErr = volErr_;
+ writer.addCellData(volErr, "volume Error");
+
+ // composition
+ ScalarSolutionType *massfraction1W = writer.template createField<Scalar, 1> (size_);
+ ScalarSolutionType *massfraction1NW = writer.template createField<Scalar, 1> (size_);
+ *massfraction1W = massfrac_[0];
+ *massfraction1NW = massfrac_[1];
+ writer.addCellData(massfraction1W, "massfraction1 in w-phase");
+ writer.addCellData(massfraction1NW, "massfraction1NW nw-phase");
+
+ // phase properties
+ ScalarSolutionType *mobilityW = writer.template createField<Scalar, 1> (size_);
+ ScalarSolutionType *mobilityNW = writer.template createField<Scalar, 1> (size_);
+ ScalarSolutionType *densityW = writer.template createField<Scalar, 1> (size_);
+ ScalarSolutionType *densityNW = writer.template createField<Scalar, 1> (size_);
+ *mobilityW = mobility_[0];
+ *mobilityNW = mobility_[1];
+ *densityW = density_[0];
+ *densityNW = density_[1];
+ writer.addCellData(mobilityW, "mobility w-phase");
+ writer.addCellData(mobilityNW, "mobility nw-phase");
+ writer.addCellData(densityW, "density w-phase");
+ writer.addCellData(densityNW, "density nw-phase");
+ }
+ if (codim_ == dim)
+ {
+ ScalarSolutionType *pressure = writer.template createField<Scalar, 1> (size_);
+ ScalarSolutionType *saturation = writer.template createField<Scalar, 1> (size_);
+
+ *pressure = this->pressure();
+ *saturation = saturation_;
+
+ writer.addVertexData(pressure, "pressure");
+ writer.addVertexData(saturation, "saturation");
+ }
+
+ return;
+ }
+
+ //! \name Access functions
+ //@{
+ //! Return the vector of the transported quantity
+ /*! For an immiscible IMPES scheme, this is the saturation. For Miscible simulations, however,
+ * the total concentration of all components is transported.
+ */
+ TransportSolutionType& transportedQuantity()
+ {
+ return totalConcentration_;
+ }
+
+ //! Returs a reference to the total concentration vector
+ const Scalar& totalConcentration() const
+ {
+ return totalConcentration_;
+ }
+
+ //! Returs a reference to the total concentration
+ /*!\param Idx Element index
+ * \param compIdx Index of the component */
+ Scalar& totalConcentration(int Idx, int compIdx)
+ {
+ return totalConcentration_[compIdx][Idx][0];
+ }
+ //! \copydoc Dumux::VariableClass2P2C::totalConcentration()
+ const Scalar& totalConcentration(int Idx, int compIdx) const
+ {
+ return totalConcentration_[compIdx][Idx][0];
+ }
+
+ //! Return saturation vector
+ ScalarSolutionType& saturation()
+ {
+ return saturation_;
+ }
+ //! Return saturation vector
+ /*! \param Idx Element index
+ * \param compIdx Index of the component */
+ const ScalarSolutionType& saturation() const
+ {
+ return saturation_;
+ }
+ //! \copydoc Dumux::VariableClass2P2C::saturation()
+ Scalar& saturation(int Idx)
+ {
+ return saturation_[Idx][0];
+ }
+
+ //! Return vector of wetting phase potential gradients
+ /*! \param Idx Element index
+ * \param isIdx Intersection index */
+ Scalar& potentialWetting(int Idx1, int isIdx)
+ {
+ return this->potential(Idx1, isIdx)[wPhaseIdx];
+ }
+ //! \copydoc Dumux::VariableClass2P2C::potentialWetting()
+ const Scalar& potentialWetting(int Idx1, int isIdx) const
+ {
+ return this->potential(Idx1, isIdx)[wPhaseIdx];
+ }
+ //! Return vector of nonwetting phase potential gradients
+ /*! \param Idx Element index
+ * \param isIdx Intersection index */
+ Scalar& potentialNonwetting(int Idx1, int isIdx)
+ {
+ return this->potential(Idx1, isIdx)[nPhaseIdx];
+ }
+ //! \copydoc Dumux::VariableClass2P2C::potentialNonwetting()
+ const Scalar& potentialNonwetting(int Idx1, int isIdx) const
+ {
+ return this->potential(Idx1, isIdx)[nPhaseIdx];
+ }
+
+ //! Return vector of wetting phase mobilities
+ /*! \param Idx Element index*/
+ Scalar& mobilityWetting(int Idx)
+ {
+ return mobility_[wPhaseIdx][Idx][0];
+ }
+ //! \copydoc Dumux::VariableClass2P2C::mobilityWetting()
+ const Scalar& mobilityWetting(int Idx) const
+ {
+ return mobility_[wPhaseIdx][Idx][0];
+ }
+
+ //! Return vector of non-wetting phase mobilities
+ /*! \param Idx Element index*/
+ Scalar& mobilityNonwetting(int Idx)
+ {
+ return mobility_[nPhaseIdx][Idx][0];
+ }
+ //! \copydoc Dumux::VariableClass2P2C::mobilityNonwetting()
+ const Scalar& mobilityNonwetting(int Idx) const
+ {
+ return mobility_[nPhaseIdx][Idx][0];
+ }
+
+ //! Return capillary pressure vector
+ /*! \param Idx Element index*/
+ Scalar& capillaryPressure(int Idx)
+ {
+ return capillaryPressure_[Idx][0];
+ }
+ //! \copydoc Dumux::VariableClass2P2C::capillaryPressure()
+ const Scalar& capillaryPressure(int Idx) const
+ {
+ return capillaryPressure_[Idx][0];
+ }
+
+ //! Return wetting phase density
+ /*! \param Idx Element index*/
+ Scalar& densityWetting(int Idx)
+ {
+ return density_[wPhaseIdx][Idx][0];
+ }
+ //! \copydoc Dumux::VariableClass2P2C::densityWetting()
+ const Scalar& densityWetting(int Idx) const
+ {
+ return density_[wPhaseIdx][Idx][0];
+ }
+
+ //! Return nonwetting phase density
+ /*! \param Idx Element index*/
+ Scalar& densityNonwetting(int Idx)
+ {
+ return density_[nPhaseIdx][Idx][0];
+ }
+ //! \copydoc Dumux::VariableClass2P2C::densityNonwetting()
+ const Scalar& densityNonwetting(int Idx) const
+ {
+ return density_[nPhaseIdx][Idx][0];
+ }
+
+ //! Return viscosity of the wetting phase
+ /*! \param Idx Element index*/
+ Scalar& viscosityWetting(int Idx)
+ {
+ return viscosity_[wPhaseIdx][Idx][0];
+ }
+ //! \copydoc Dumux::VariableClass2P2C::viscosityWetting()
+ const Scalar& viscosityWetting(int Idx) const
+ {
+ return viscosity_[wPhaseIdx][Idx][0];
+ }
+
+ //! Return viscosity of the nonwetting phase
+ /*! \param Idx Element index*/
+ Scalar& viscosityNonwetting(int Idx)
+ {
+ return viscosity_[nPhaseIdx][Idx][0];
+ }
+ //! \copydoc Dumux::VariableClass2P2C::viscosityNonwetting()
+ const Scalar& viscosityNonwetting(int Idx) const
+ {
+ return viscosity_[nPhaseIdx][Idx][0];
+ }
+
+ //! Returns a reference to the vector holding the volume error
+ ScalarSolutionType& volErr()
+ {
+ return volErr_;
+ }
+ //! Returns a reference to the volume error
+ /*! \param Idx Element index*/
+ const Scalar& volErr(int Idx) const
+ {
+ return volErr_[Idx][0];
+ }
+
+
+ //! Returns a reference to mass fraction of first component in wetting phase
+ /*! \param Idx Element index*/
+ Scalar& wet_X1(int Idx)
+ {
+ return massfrac_[wPhaseIdx][Idx][0];
+ }
+ //! Returns a reference to mass fraction of first component in nonwetting phase
+ /*! \param Idx Element index*/
+ Scalar& nonwet_X1(int Idx)
+ {
+ return massfrac_[nPhaseIdx][Idx][0];
+ }
+ //! Returns a reference to mass fractions of first component
+ /*!\param Idx Element index
+ * \param phaseIdx Index of the phase */
+ const Scalar& massFracComp1(int Idx, int phaseIdx) const
+ {
+ return massfrac_[phaseIdx][Idx][0];
+ }
+
+ //! Return the estimated update vector for pressure calculation
+ TransportSolutionType& updateEstimate()
+ {
+ return updateEstimate_;
+ }
+ //! Return the estimate of the next update
+ /*!\param Idx Element index
+ * \param compIdx Index of the component */
+ Scalar& updateEstimate(int Idx, int compIdx)
+ {
+ return updateEstimate_[compIdx][Idx][0];
+ }
+
+ //! Return the vector holding subdomain information
+ Dune::BlockVector<Dune::FieldVector<int,1> >& subdomain()
+ {
+ return subdomain_;
+ }
+ //! Return specific subdomain information
+ /*!\param Idx Element index */
+ int& subdomain(int Idx)
+ {
+ return subdomain_[Idx][0];
+ }
+ //@}
+
+};
+}
+#endif
diff --git a/dumux/decoupled/Makefile.am b/dumux/decoupled/Makefile.am
index f704a24b393a9bad2c25084ba52a9e72f6bcc512..4c408a528e397312fcbe433a426e76a8b844a2e8 100644
--- a/dumux/decoupled/Makefile.am
+++ b/dumux/decoupled/Makefile.am
@@ -1,4 +1,4 @@
-SUBDIRS = 2p common
+SUBDIRS = 1p 2p 2p2c common
decoupleddir = $(includedir)/dumux/decoupled
diff --git a/test/decoupled/2p2c/Makefile.am b/test/decoupled/2p2c/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..4b55af671edbf9158febb3c52ccbc89153e8e3e6
--- /dev/null
+++ b/test/decoupled/2p2c/Makefile.am
@@ -0,0 +1,11 @@
+# tests where program to build and program to run are equal
+NORMALTESTS = test_dec2p2c test_multiphysics2p2c
+
+# programs just to build when "make check" is used
+check_PROGRAMS = $(NORMALTESTS)
+
+test_dec2p2c_SOURCES = test_dec2p2c.cc
+
+test_multiphysics2p2c_SOURCES = test_multiphysics2p2c.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/decoupled/2p2c/test_dec2p2c.cc b/test/decoupled/2p2c/test_dec2p2c.cc
new file mode 100644
index 0000000000000000000000000000000000000000..77add08fe4f397e9468213f5428a2a5586117456
--- /dev/null
+++ b/test/decoupled/2p2c/test_dec2p2c.cc
@@ -0,0 +1,119 @@
+// $Id: test_2pinjection.cc 3151 2010-02-15 11:41:23Z mwolff $
+/*****************************************************************************
+ * Copyright (C) 20010 by Markus Wolff *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+
+#include "test_dec2p2cproblem.hh"
+
+#include <dune/grid/common/gridinfo.hh>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/mpihelper.hh>
+
+#include <iostream>
+#include <boost/format.hpp>
+
+////////////////////////
+// the main function
+////////////////////////
+void usage(const char *progname)
+{
+ std::cout << boost::format("usage: %s [--restart restartTime] tEnd firstDt\n")%progname;
+ exit(1);
+}
+
+int main(int argc, char** argv)
+{
+ try {
+ typedef TTAG(TestTwoPTwoCProblem) TypeTag;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef Dune::FieldVector<Scalar, Grid::dimensionworld> GlobalPosition;
+
+ static const int dim = Grid::dimension;
+
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper::instance(argc, argv);
+
+ ////////////////////////////////////////////////////////////
+ // parse the command line arguments
+ ////////////////////////////////////////////////////////////
+ // deal with the restart stuff
+ int argPos = 1;
+ bool restart = false;
+ double restartTime = 0;
+ // deal with start parameters
+ double tEnd= 3e3;
+ double firstDt = 2e3;
+ if (argc != 1)
+ {
+ // deal with the restart stuff
+ if (std::string("--restart") == argv[argPos]) {
+ restart = true;
+ ++argPos;
+
+ std::istringstream(argv[argPos++]) >> restartTime;
+ }
+ if (argc - argPos == 2)
+ {
+ // read the initial time step and the end time
+ std::istringstream(argv[argPos++]) >> tEnd;
+ std::istringstream(argv[argPos++]) >> firstDt;
+ }
+ else
+ usage(argv[0]);
+ }
+ else
+ {
+ Dune::dwarn << "simulation started with predefs" << std::endl;
+ }
+
+ ////////////////////////////////////////////////////////////
+ // create the grid
+ ////////////////////////////////////////////////////////////
+ Dune::FieldVector<int,dim> N(10);
+ Dune::FieldVector<double ,dim> L(0);
+ Dune::FieldVector<double,dim> H(10);
+ Grid grid(N,L,H);
+
+ ////////////////////////////////////////////////////////////
+ // instantiate and run the concrete problem
+ ////////////////////////////////////////////////////////////
+
+ Problem problem(grid.leafView(), L, H);
+
+ // load restart file if necessarry
+ if (restart)
+ problem.deserialize(restartTime);
+
+ // initialize the simulation
+ problem.timeManager().init(problem, 0, firstDt, tEnd, !restart);
+ // run the simulation
+ problem.timeManager().run();
+ return 0;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ throw;
+ }
+
+ return 3;
+}
diff --git a/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh b/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f08bed0e564a660cee01894ff2433d413b1df373
--- /dev/null
+++ b/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh
@@ -0,0 +1,107 @@
+// $Id: test_2p_spatialparamsinjection.hh 3456 2010-04-09 12:11:51Z mwolff $
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Markus Wolff *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef TEST_2P2C_SPATIALPARAMETERS_HH
+#define TEST_2P2C_SPATIALPARAMETERS_HH
+
+
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+//#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+namespace Dumux
+{
+
+/** \todo Please doc me! */
+
+template<class TypeTag>
+class Test2P2CSpatialParams
+{
+ 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=Grid::dimension, dimWorld=Grid::dimensionworld, numEq=1};
+ typedef typename Grid::Traits::template Codim<0>::Entity Element;
+
+ typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<CoordScalar, dim> LocalPosition;
+ typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix;
+
+// typedef RegularizedBrooksCorey<Scalar> RawMaterialLaw;
+ typedef LinearMaterial<Scalar> RawMaterialLaw;
+public:
+ typedef EffToAbsLaw<RawMaterialLaw> MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ void update (Scalar saturationW, const Element& element)
+ {
+
+ }
+
+ const FieldMatrix& intrinsicPermeability (const GlobalPosition& globalPos, const Element& element) const
+ {
+ return constPermeability_;
+ }
+
+ double porosity(const GlobalPosition& globalPos, const Element& element) const
+ {
+ return 0.2;
+ }
+
+
+ // return the brooks-corey context depending on the position
+ const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
+ {
+ return materialLawParams_;
+ }
+
+
+ Test2P2CSpatialParams(const GridView& gridView)
+ : constPermeability_(0)
+ {
+ // residual saturations
+ materialLawParams_.setSwr(0);
+ materialLawParams_.setSnr(0);
+
+// // parameters for the Brooks-Corey Law
+// // entry pressures
+// materialLawParams_.setPe(10000);
+//
+// // Brooks-Corey shape parameters
+// materialLawParams_.setAlpha(2);
+
+ // parameters for the linear
+ // entry pressures function
+
+ materialLawParams_.setEntryPC(0);
+ materialLawParams_.setMaxPC(0);
+
+ for(int i = 0; i < dim; i++)
+ {
+ constPermeability_[i][i] = 1e-12;
+ }
+ }
+
+private:
+ MaterialLawParams materialLawParams_;
+ FieldMatrix constPermeability_;
+
+};
+
+} // end namespace
+#endif
diff --git a/test/decoupled/2p2c/test_dec2p2cproblem.hh b/test/decoupled/2p2c/test_dec2p2cproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..017e0bc194b63455344593db08964f79e88d45f8
--- /dev/null
+++ b/test/decoupled/2p2c/test_dec2p2cproblem.hh
@@ -0,0 +1,269 @@
+// $Id: test_2p_injectionproblem.hh 3456 2010-04-09 12:11:51Z mwolff $
+/*****************************************************************************
+ * Copyright (C) 2007-2008 by Klaus Mosthaf *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_TEST_2P2C_PROBLEM_HH
+#define DUMUX_TEST_2P2C_PROBLEM_HH
+
+#if HAVE_UG
+#include <dune/grid/uggrid.hh>
+#endif
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/sgrid.hh>
+
+// fluid properties
+//#include <dumux/material/fluidsystems/simple_h2o_n2_system.hh>
+#include <dumux/material/fluidsystems/h2o_n2_system.hh>
+
+#include <dumux/decoupled/2p2c/2p2cproblem.hh>
+#include <dumux/decoupled/2p2c/fvpressure2p2c.hh>
+#include <dumux/decoupled/2p2c/fvtransport2p2c.hh>
+
+#include "test_dec2p2c_spatialparams.hh"
+
+namespace Dumux
+{
+
+template<class TypeTag>
+class TestTwoPTwoCProblem;
+
+//////////
+// Specify the properties for the lens problem
+//////////
+namespace Properties
+{
+NEW_TYPE_TAG(TestTwoPTwoCProblem, INHERITS_FROM(DecoupledTwoPTwoC/*, Transport*/));
+
+// Set the grid type
+SET_PROP(TestTwoPTwoCProblem, Grid)
+{
+ // typedef Dune::YaspGrid<2> type;
+ typedef Dune::SGrid<3, 3> type;
+};
+
+// Set the problem property
+SET_PROP(TestTwoPTwoCProblem, Problem)
+{
+ typedef Dumux::TestTwoPTwoCProblem<TTAG(TestTwoPTwoCProblem)> type;
+};
+
+// Set the model properties
+SET_PROP(TestTwoPTwoCProblem, TransportModel)
+{
+ typedef Dumux::FVTransport2P2C<TTAG(TestTwoPTwoCProblem)> type;
+};
+
+SET_PROP(TestTwoPTwoCProblem, PressureModel)
+{
+ typedef Dumux::FVPressure2P2C<TTAG(TestTwoPTwoCProblem)> type;
+};
+
+SET_INT_PROP(TestTwoPTwoCProblem, VelocityFormulation,
+ GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::velocityW);
+
+SET_INT_PROP(TestTwoPTwoCProblem, PressureFormulation,
+ GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::pressureW);
+
+
+// Select fluid system
+SET_PROP(TestTwoPTwoCProblem, FluidSystem)
+{
+ typedef Dumux::H2O_N2_System<TypeTag> type;
+};
+
+// Select water formulation
+SET_PROP(TestTwoPTwoCProblem, Components) : public GET_PROP(TypeTag, PTAG(DefaultComponents))
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+// typedef Dumux::TabulatedComponent<Scalar, typename Dumux::H2O<Scalar> > H20;
+ typedef Dumux::SimpleH2O<Scalar> H2O;
+};
+
+// Set the soil properties
+SET_PROP(TestTwoPTwoCProblem, SpatialParameters)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+public:
+ typedef Dumux::Test2P2CSpatialParams<TypeTag> type;
+};
+
+// Enable gravity
+SET_BOOL_PROP(TestTwoPTwoCProblem, EnableGravity, true);
+SET_INT_PROP(DecoupledTwoPTwoC,
+ BoundaryMobility,
+ TwoPCommonIndices<TypeTag>::satDependent);
+SET_SCALAR_PROP(TestTwoPTwoCProblem, CFLFactor, 0.8);
+}
+
+/*!
+ * \ingroup DecoupledProblems
+ */
+template<class TypeTag = TTAG(TestTwoPTwoCProblem)>
+class TestTwoPTwoCProblem: public IMPETProblem2P2C<TypeTag, TestTwoPTwoCProblem<TypeTag> >
+{
+typedef TestTwoPTwoCProblem<TypeTag> ThisType;
+typedef IMPETProblem2P2C<TypeTag, ThisType> ParentType;
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+
+enum
+{
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+};
+
+enum
+{
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+};
+
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+typedef typename GridView::Traits::template Codim<0>::Entity Element;
+typedef typename GridView::Intersection Intersection;
+typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+
+public:
+TestTwoPTwoCProblem(const GridView &gridView, const GlobalPosition lowerLeft = 0, const GlobalPosition upperRight = 0) :
+ParentType(gridView), lowerLeft_(lowerLeft), upperRight_(upperRight)
+{
+ // initialize the tables of the fluid system
+// WaterFormulation::init(273.15, 623.15, 100,
+// -10, 20e6, 200);
+ FluidSystem::init();
+}
+
+/*!
+ * \name Problem parameters
+ */
+// \{
+
+/*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+const char *name() const
+{
+ return "test_dec2p2c";
+}
+
+bool shouldWriteRestartFile() const
+{
+ return false;
+}
+
+/*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+Scalar temperature(const GlobalPosition& globalPos, const Element& element) const
+{
+ return 273.15 + 10; // -> 10°C
+}
+
+// \}
+
+Scalar referencePressure(const GlobalPosition& globalPos, const Element& element) const
+{
+ return 1e6; // TODO: proper description!!
+}
+
+typename BoundaryConditions::Flags bcTypePress(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ if (globalPos[0] > 10-1E-6 || globalPos[0] < 1e-6)
+ return BoundaryConditions::dirichlet;
+ // all other boundaries
+ return BoundaryConditions::neumann;
+}
+
+typename BoundaryConditions::Flags bcTypeTransport(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ return bcTypePress(globalPos, intersection);
+}
+
+BoundaryConditions2p2c::Flags bcFormulation(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ return BoundaryConditions2p2c::concentration;
+}
+
+Scalar dirichletPress(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ const Element& element = *(intersection.inside());
+
+ Scalar pRef = referencePressure(globalPos, element);
+ Scalar temp = temperature(globalPos, element);
+
+ // all other boundaries
+ return (globalPos[0] < 1e-6) ? (2.5e5 - FluidSystem::H2O::liquidDensity(temp, pRef) * this->gravity()[dim-1])
+ : (2e5 - FluidSystem::H2O::liquidDensity(temp, pRef) * this->gravity()[dim-1]);
+}
+
+Scalar dirichletTransport(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ return 1.;
+}
+
+Dune::FieldVector<Scalar,2> neumann(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ Dune::FieldVector<Scalar,2> neumannFlux(0.0);
+// if (globalPos[1] < 15 && globalPos[1]> 5)
+// {
+// neumannFlux[nPhaseIdx] = -0.015;
+// }
+ return neumannFlux;
+}
+
+Dune::FieldVector<Scalar,2> source(const GlobalPosition& globalPos, const Element& element)
+{
+ Dune::FieldVector<Scalar,2> q_(0);
+ if (fabs(globalPos[0] - 4.5) < 1 && fabs(globalPos[1] - 4.5) < 1) q_[1] = 0.0001;
+ return q_;
+}
+
+const BoundaryConditions2p2c::Flags initFormulation (const GlobalPosition& globalPos, const Element& element) const
+{
+ return BoundaryConditions2p2c::concentration;
+}
+
+Scalar initSat(const GlobalPosition& globalPos, const Element& element) const
+{
+ return 1;
+}
+Scalar initConcentration(const GlobalPosition& globalPos, const Element& element) const
+{
+ return 1;
+}
+
+private:
+GlobalPosition lowerLeft_;
+GlobalPosition upperRight_;
+
+static const Scalar eps_ = 1e-6;
+static const Scalar depthBOR_ = 1000;
+};
+} //end namespace
+
+#endif
diff --git a/test/decoupled/2p2c/test_multiphysics2p2c.cc b/test/decoupled/2p2c/test_multiphysics2p2c.cc
new file mode 100644
index 0000000000000000000000000000000000000000..44e5950e42ac4029a8772a7870c489f8303027f7
--- /dev/null
+++ b/test/decoupled/2p2c/test_multiphysics2p2c.cc
@@ -0,0 +1,119 @@
+// $Id:
+/*****************************************************************************
+ * Copyright (C) 20010 by Benjamin Faigle *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+
+#include "test_multiphysics2p2cproblem.hh"
+
+#include <dune/grid/common/gridinfo.hh>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/mpihelper.hh>
+
+#include <iostream>
+#include <boost/format.hpp>
+
+////////////////////////
+// the main function
+////////////////////////
+void usage(const char *progname)
+{
+ std::cout << boost::format("usage: %s [--restart restartTime] tEnd firstDt\n")%progname;
+ exit(1);
+}
+
+int main(int argc, char** argv)
+{
+ try {
+ typedef TTAG(TestTwoPTwoCProblem) TypeTag;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef Dune::FieldVector<Scalar, Grid::dimensionworld> GlobalPosition;
+
+ static const int dim = Grid::dimension;
+
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper::instance(argc, argv);
+
+ ////////////////////////////////////////////////////////////
+ // parse the command line arguments
+ ////////////////////////////////////////////////////////////
+ // deal with the restart stuff
+ int argPos = 1;
+ bool restart = false;
+ double restartTime = 0;
+ // deal with start parameters
+ double tEnd= 3e3;
+ double firstDt = 2e3;
+ if (argc != 1)
+ {
+ // deal with the restart stuff
+ if (std::string("--restart") == argv[argPos]) {
+ restart = true;
+ ++argPos;
+
+ std::istringstream(argv[argPos++]) >> restartTime;
+ }
+ if (argc - argPos == 2)
+ {
+ // read the initial time step and the end time
+ std::istringstream(argv[argPos++]) >> tEnd;
+ std::istringstream(argv[argPos++]) >> firstDt;
+ }
+ else
+ usage(argv[0]);
+ }
+ else
+ {
+ Dune::dwarn << "simulation started with predefs" << std::endl;
+ }
+
+ ////////////////////////////////////////////////////////////
+ // create the grid
+ ////////////////////////////////////////////////////////////
+ Dune::FieldVector<int,dim> N(10);
+ Dune::FieldVector<double ,dim> L(0);
+ Dune::FieldVector<double,dim> H(10);
+ Grid grid(N,L,H);
+
+ ////////////////////////////////////////////////////////////
+ // instantiate and run the concrete problem
+ ////////////////////////////////////////////////////////////
+
+ Problem problem(grid.leafView(), L, H);
+
+ // load restart file if necessarry
+ if (restart)
+ problem.deserialize(restartTime);
+
+ // initialize the simulation
+ problem.timeManager().init(problem, 0, firstDt, tEnd, !restart);
+ // run the simulation
+ problem.timeManager().run();
+ return 0;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ throw;
+ }
+
+ return 3;
+}
diff --git a/test/decoupled/2p2c/test_multiphysics2p2cproblem.hh b/test/decoupled/2p2c/test_multiphysics2p2cproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..74dabd4284ecf0b0e211e8f4e328f11e0b5c27ec
--- /dev/null
+++ b/test/decoupled/2p2c/test_multiphysics2p2cproblem.hh
@@ -0,0 +1,270 @@
+// $Id: test_2p_injectionproblem.hh 3456 2010-04-09 12:11:51Z mwolff $
+/*****************************************************************************
+ * Copyright (C) 2007-2008 by Klaus Mosthaf *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_TEST_2P2C_PROBLEM_HH
+#define DUMUX_TEST_2P2C_PROBLEM_HH
+
+#if HAVE_UG
+#include <dune/grid/uggrid.hh>
+#endif
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/sgrid.hh>
+
+// fluid properties
+//#include <dumux/material/fluidsystems/simple_h2o_n2_system.hh>
+#include <dumux/material/fluidsystems/h2o_n2_system.hh>
+
+#include <dumux/decoupled/2p2c/2p2cproblem.hh>
+#include <dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh>
+#include <dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh>
+
+#include "test_dec2p2c_spatialparams.hh"
+
+namespace Dumux
+{
+
+template<class TypeTag>
+class TestTwoPTwoCProblem;
+
+//////////
+// Specify the properties for the lens problem
+//////////
+namespace Properties
+{
+NEW_TYPE_TAG(TestTwoPTwoCProblem, INHERITS_FROM(DecoupledTwoPTwoC/*, Transport*/));
+
+// Set the grid type
+SET_PROP(TestTwoPTwoCProblem, Grid)
+{
+ // typedef Dune::YaspGrid<2> type;
+ typedef Dune::SGrid<3, 3> type;
+};
+
+// Set the problem property
+SET_PROP(TestTwoPTwoCProblem, Problem)
+{
+ typedef Dumux::TestTwoPTwoCProblem<TTAG(TestTwoPTwoCProblem)> type;
+};
+
+// Set the model properties
+SET_PROP(TestTwoPTwoCProblem, TransportModel)
+{
+ typedef Dumux::FVTransport2P2CMultiPhysics<TTAG(TestTwoPTwoCProblem)> type;
+};
+
+SET_PROP(TestTwoPTwoCProblem, PressureModel)
+{
+ typedef Dumux::FVPressure2P2CMultiPhysics<TTAG(TestTwoPTwoCProblem)> type;
+};
+
+SET_INT_PROP(TestTwoPTwoCProblem, VelocityFormulation,
+ GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::velocityW);
+
+SET_INT_PROP(TestTwoPTwoCProblem, PressureFormulation,
+ GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::pressureW);
+
+
+// Select fluid system
+SET_PROP(TestTwoPTwoCProblem, FluidSystem)
+{
+ typedef Dumux::H2O_N2_System<TypeTag> type;
+};
+
+// Select water formulation
+SET_PROP(TestTwoPTwoCProblem, Components) : public GET_PROP(TypeTag, PTAG(DefaultComponents))
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+// typedef Dumux::TabulatedComponent<Scalar, typename Dumux::H2O<Scalar> > type;
+ typedef Dumux::H2O<Scalar> H2O;
+// static void init(){};
+};
+
+// Set the soil properties
+SET_PROP(TestTwoPTwoCProblem, SpatialParameters)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+public:
+ typedef Dumux::Test2P2CSpatialParams<TypeTag> type;
+};
+
+// Enable gravity
+SET_BOOL_PROP(TestTwoPTwoCProblem, EnableGravity, true);
+SET_INT_PROP(DecoupledTwoPTwoC,
+ BoundaryMobility,
+ TwoPCommonIndices<TypeTag>::satDependent);
+SET_SCALAR_PROP(TestTwoPTwoCProblem, CFLFactor, 0.8);
+}
+
+/*!
+ * \ingroup DecoupledProblems
+ */
+template<class TypeTag = TTAG(TestTwoPTwoCProblem)>
+class TestTwoPTwoCProblem: public IMPETProblem2P2C<TypeTag, TestTwoPTwoCProblem<TypeTag> >
+{
+typedef TestTwoPTwoCProblem<TypeTag> ThisType;
+typedef IMPETProblem2P2C<TypeTag, ThisType> ParentType;
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+
+enum
+{
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+};
+
+enum
+{
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+};
+
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+typedef typename GridView::Traits::template Codim<0>::Entity Element;
+typedef typename GridView::Intersection Intersection;
+typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+
+public:
+TestTwoPTwoCProblem(const GridView &gridView, const GlobalPosition lowerLeft = 0, const GlobalPosition upperRight = 0) :
+ParentType(gridView), lowerLeft_(lowerLeft), upperRight_(upperRight)
+{
+ // initialize the tables of the fluid system
+// WaterFormulation::init(273.15, 623.15, 100,
+// -10, 20e6, 200);
+ FluidSystem::init();
+}
+
+/*!
+ * \name Problem parameters
+ */
+// \{
+
+/*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+const char *name() const
+{
+ return "test_multiphysics2p2c";
+}
+
+bool shouldWriteRestartFile() const
+{
+ return false;
+}
+
+/*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+Scalar temperature(const GlobalPosition& globalPos, const Element& element) const
+{
+ return 273.15 + 10; // -> 10°C
+}
+
+// \}
+
+Scalar referencePressure(const GlobalPosition& globalPos, const Element& element) const
+{
+ return 1e6; // TODO: proper description!!
+}
+
+typename BoundaryConditions::Flags bcTypePress(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ if (globalPos[0] > 10-1E-6 || globalPos[0] < 1e-6)
+ return BoundaryConditions::dirichlet;
+ // all other boundaries
+ return BoundaryConditions::neumann;
+}
+
+typename BoundaryConditions::Flags bcTypeTransport(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ return bcTypePress(globalPos, intersection);
+}
+
+BoundaryConditions2p2c::Flags bcFormulation(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ return BoundaryConditions2p2c::concentration;
+}
+
+Scalar dirichletPress(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ const Element& element = *(intersection.inside());
+
+ Scalar pRef = referencePressure(globalPos, element);
+ Scalar temp = temperature(globalPos, element);
+
+ // all other boundaries
+ return (globalPos[0] < 1e-6) ? (2.5e5 - FluidSystem::H2O::liquidDensity(temp, pRef) * this->gravity()[dim-1])
+ : (2e5 - FluidSystem::H2O::liquidDensity(temp, pRef) * this->gravity()[dim-1]);
+}
+
+Scalar dirichletTransport(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ return 1.;
+}
+
+Dune::FieldVector<Scalar,2> neumann(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ Dune::FieldVector<Scalar,2> neumannFlux(0.0);
+// if (globalPos[1] < 15 && globalPos[1]> 5)
+// {
+// neumannFlux[nPhaseIdx] = -0.015;
+// }
+ return neumannFlux;
+}
+
+Dune::FieldVector<Scalar,2> source(const GlobalPosition& globalPos, const Element& element)
+{
+ Dune::FieldVector<Scalar,2> q_(0);
+ if (fabs(globalPos[0] - 4.5) < 1 && fabs(globalPos[1] - 4.5) < 1) q_[1] = 0.0001;
+ return q_;
+}
+
+const BoundaryConditions2p2c::Flags initFormulation (const GlobalPosition& globalPos, const Element& element) const
+{
+ return BoundaryConditions2p2c::concentration;
+}
+
+Scalar initSat(const GlobalPosition& globalPos, const Element& element) const
+{
+ return 1;
+}
+Scalar initConcentration(const GlobalPosition& globalPos, const Element& element) const
+{
+ return 1;
+}
+
+private:
+GlobalPosition lowerLeft_;
+GlobalPosition upperRight_;
+
+static const Scalar eps_ = 1e-6;
+static const Scalar depthBOR_ = 1000;
+};
+} //end namespace
+
+#endif
diff --git a/test/decoupled/Makefile.am b/test/decoupled/Makefile.am
index 09685656b90442673ad795c4f28deb09d09ae023..97a27f90a80a7ea1f105f37976ab7e94e0c1fe57 100644
--- a/test/decoupled/Makefile.am
+++ b/test/decoupled/Makefile.am
@@ -1,5 +1,6 @@
SUBDIRS = 1p \
- 2p
+ 2p \
+ 2p2c
include $(top_srcdir)/am/global-rules