From 17cd1b63e0e13699eba85dbe8e960b351539c2c0 Mon Sep 17 00:00:00 2001
From: Benjamin Faigle <benjamin.faigle@posteo.de>
Date: Fri, 28 Sep 2012 14:48:55 +0000
Subject: [PATCH] Adding compositional adaptive module with mpfa on hanging
nodes (as announced in last Dumux meeting). A test and reference is also
provided. The module is stable, but doxygen docu is still needing some
improvement.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@9145 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/decoupled/2p2c/cellData2p2cadaptive.hh | 276 ++++
.../decoupled/2p2c/fvpressure2p2cadaptive.hh | 1277 +++++++++++++++++
.../decoupled/2p2c/fvtransport2p2cadaptive.hh | 526 +++++++
.../2p2c/variableclass2p2cadaptive.hh | 301 ++++
test/decoupled/2p2c/Makefile.am | 4 +-
.../2p2c/test_adaptive2p2c-reference.vtu | 457 ++++++
test/decoupled/2p2c/test_adaptive2p2c.cc | 57 +
test/decoupled/2p2c/test_adaptive2p2c.input | 43 +
.../2p2c/test_adaptive2p2cproblem.hh | 294 ++++
9 files changed, 3233 insertions(+), 2 deletions(-)
create mode 100644 dumux/decoupled/2p2c/cellData2p2cadaptive.hh
create mode 100644 dumux/decoupled/2p2c/fvpressure2p2cadaptive.hh
create mode 100644 dumux/decoupled/2p2c/fvtransport2p2cadaptive.hh
create mode 100644 dumux/decoupled/2p2c/variableclass2p2cadaptive.hh
create mode 100644 test/decoupled/2p2c/test_adaptive2p2c-reference.vtu
create mode 100644 test/decoupled/2p2c/test_adaptive2p2c.cc
create mode 100644 test/decoupled/2p2c/test_adaptive2p2c.input
create mode 100644 test/decoupled/2p2c/test_adaptive2p2cproblem.hh
diff --git a/dumux/decoupled/2p2c/cellData2p2cadaptive.hh b/dumux/decoupled/2p2c/cellData2p2cadaptive.hh
new file mode 100644
index 0000000000..d289fc1664
--- /dev/null
+++ b/dumux/decoupled/2p2c/cellData2p2cadaptive.hh
@@ -0,0 +1,276 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Markus Wolff *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+#ifndef DUMUX_ELEMENTDATA2P2C_ADAPTIVE_HH
+#define DUMUX_ELEMENTDATA2P2C_ADAPTIVE_HH
+
+#include <dumux/decoupled/2p2c/cellData2p2c.hh>
+#include <dumux/decoupled/2p2c/cellData2p2cmultiphysics.hh>
+/**
+ * @file
+ * @brief Class including the variables and data of discretized data of the constitutive relations for one element
+ * @author Benjamin Faigle, Markus Wolff
+ */
+namespace Dumux
+{
+/*!
+ * \ingroup IMPET
+ */
+//! Class including the data of a grid cell needed if an adaptive grid is used.
+/*! The class provides model-specific functions needed to adapt the stored cell data to a new (adapted) grid.
+ * Additionally, it provides the storage-infrastructure for explicit front tracking.
+ *
+ * @tparam TypeTag The Type Tag
+ */
+template<class TypeTag>
+class CellData2P2CAdaptive: public CellData2P2Cmultiphysics<TypeTag>
+{
+private:
+ typedef CellData2P2CAdaptive<TypeTag> ThisType;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GridView::Grid Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
+ wCompIdx = Indices::wCompIdx, nCompIdx = Indices::nCompIdx
+ };
+ enum
+ {
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
+ };
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ static constexpr int pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation); //!< gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
+ int upwindError_[numPhases];
+
+public:
+ //! A container for all necessary variables to map an old solution to a new grid
+ /*
+ * If the primary variables (pressure, total concentrations) are mapped to a new grid,
+ * the secondary variables can be calulated. For the mapping between sons and father, it
+ * is in addition necessary to know about how many sons live in each father ("count").
+ * While only one phase pressure is PV, the method updateMaterialLaws() that
+ * recalculates the secondary variables needs both phase pressures (initiated via the
+ * capillary pressure of the last time step) to start the iteration to determine
+ * appropriate phase composition and pc. Hence, both phase pressures are mapped
+ * to the new solution.
+ */
+ struct AdaptedValues
+ {
+ Dune::FieldVector<Scalar,2> totalConcentration_;
+ Dune::FieldVector<Scalar,2> pressure_;
+ Dune::FieldVector<Scalar,3> volumeDerivatives_;
+ Scalar cellVolume;
+ FluxData2P2C<TypeTag> fluxData_;
+ int subdomain;
+ int count;
+ int isRefined;
+ AdaptedValues()
+ {
+ totalConcentration_=0.;
+ pressure_ = 0.;
+ volumeDerivatives_ =0.;
+ cellVolume = 0.;
+ subdomain=-1;
+ count = 0;
+ isRefined = false;
+ }
+ };
+
+
+ //! Constructs an adaptive CellData object
+ CellData2P2CAdaptive() : CellData2P2Cmultiphysics<TypeTag>()
+ {
+ for (int i = 0; i < numPhases;i++)
+ upwindError_[i] = 0;
+ }
+
+ //! stores leaf cell primary variables to transfer to new indexing
+ /**
+ * Stores values to be adapted from the current CellData objects into
+ * the adaptation container in order to be mapped on a new grid.
+ *
+ * @param adaptedValues Container for model-specific values to be adapted
+ * @param element The element to be stored
+ */
+ void storeAdaptionValues(AdaptedValues& adaptedValues, const Element& element)
+ {
+ adaptedValues.totalConcentration_[wCompIdx]= this->massConcentration(wCompIdx);
+ adaptedValues.totalConcentration_[nCompIdx]= this->massConcentration(nCompIdx);
+ adaptedValues.pressure_[wPhaseIdx] = this->pressure(wPhaseIdx);
+ adaptedValues.pressure_[nPhaseIdx] = this->pressure(nPhaseIdx);
+ adaptedValues.volumeDerivatives_[wCompIdx] = this->dv(wPhaseIdx);
+ adaptedValues.volumeDerivatives_[nCompIdx] = this->dv(nPhaseIdx);
+ adaptedValues.volumeDerivatives_[2] = this->dv_dp();
+ adaptedValues.cellVolume = element.geometry().volume();
+ adaptedValues.subdomain = this->subdomain();
+ adaptedValues.fluxData_=this->fluxData();
+ }
+
+ //! adds cell information to father element for possible averageing / coarsening
+ /**
+ * Sum up the adaptedValues (sons values) into father element. We store from leaf
+ * upwards, so sons are stored first, then cells on the next leaf (=fathers)
+ * can be averaged.
+ *
+ * @param adaptedValues Container for model-specific values to be adapted
+ * @param adaptedValuesFather Values to be adapted of father cell
+ * @param fatherElement The element of the father
+ */
+ static void storeAdaptionValues(AdaptedValues& adaptedValues,
+ AdaptedValues& adaptedValuesFather,
+ const Element& fatherElement)
+ {
+ adaptedValuesFather.totalConcentration_[wCompIdx]
+ += adaptedValues.cellVolume* adaptedValues.totalConcentration_[wCompIdx];
+ adaptedValuesFather.totalConcentration_[nCompIdx]
+ += adaptedValues.cellVolume* adaptedValues.totalConcentration_[nCompIdx];
+ // if all cells are summed up, re-convert mass into total concentrations
+ Scalar fatherVolume = fatherElement.geometry().volume();
+ if(adaptedValuesFather.count == (pow(2,dim)))
+ {
+ adaptedValuesFather.totalConcentration_[wCompIdx] /= fatherVolume;
+ adaptedValuesFather.totalConcentration_[nCompIdx] /= fatherVolume;
+ }
+ adaptedValuesFather.cellVolume = fatherVolume;
+
+ adaptedValuesFather.pressure_[wPhaseIdx] += adaptedValues.pressure_[wPhaseIdx] / adaptedValues.count;
+ adaptedValuesFather.pressure_[nPhaseIdx] += adaptedValues.pressure_[nPhaseIdx] / adaptedValues.count;
+ // apply maximum complexity for new cell
+ adaptedValuesFather.subdomain = std::max(adaptedValuesFather.subdomain, adaptedValues.subdomain);
+ }
+ //! Set adapted values in CellData
+ /**
+ * This methods stores reconstructed values into the cellData object, by
+ * this setting a newly mapped solution to the storage container of the
+ * decoupled models.
+ * In new cells, update estimate does not give meaningful results. We therefore
+ * copy volume derivatives from old time step, and indicate that those are already availabe.
+ *
+ * @param adaptedValues Container for model-specific values to be adapted
+ * @param element The element where things are stored.
+ */
+ void setAdaptionValues(AdaptedValues& adaptedValues, const Element& element)
+ {
+ // in new cells, there is a cellData object, but not yet a fluidstate.
+ // To write the adapted values in this fluidstate, we have to enshure that
+ // it was created. We therefore specify the type of FluidState that should be stored.
+ this->setSubdomainAndFluidStateType(adaptedValues.subdomain);
+ this->setMassConcentration(wCompIdx,
+ adaptedValues.totalConcentration_[wCompIdx]);
+ this->setMassConcentration(nCompIdx,
+ adaptedValues.totalConcentration_[nCompIdx]);
+ this->setPressure(wPhaseIdx,
+ adaptedValues.pressure_[wPhaseIdx] / adaptedValues.count);
+ this->setPressure(nPhaseIdx,
+ adaptedValues.pressure_[nPhaseIdx] / adaptedValues.count);
+
+ //copy flux directions
+ this->fluxData()=adaptedValues.fluxData_;
+
+ //indicate if cell was just refined in this TS
+ this->wasRefined()= adaptedValues.isRefined;
+ if(adaptedValues.isRefined)
+ {
+ this->dv(wPhaseIdx) = adaptedValues.volumeDerivatives_[wCompIdx];
+ this->dv(nPhaseIdx) = adaptedValues.volumeDerivatives_[nCompIdx];
+ this->dv_dp() = adaptedValues.volumeDerivatives_[2];
+ this->volumeDerivativesAvailable(true);
+ }
+ else
+ this->volumeDerivativesAvailable(false); // recalculate volume derivatives
+ }
+ //! Reconstructs sons entries from data of father cell
+ /**
+ * Reconstructs an new solution from a father cell into for a newly
+ * generated son cell. The new cell is stored into the global
+ * adaptationMap.
+ *
+ * @param adaptationMap Global map storing all values to be adapted
+ * @param father Entity Pointer to the father cell
+ * @param son Entity Pointer to the newly created son cell
+ * @param problem The problem
+ */
+ static void reconstructAdaptionValues(Dune::PersistentContainer<Grid, AdaptedValues>& adaptationMap,
+ const Element& father, const Element& son, const Problem& problem)
+ {
+ // acess father and son
+ AdaptedValues& adaptedValues = adaptationMap[son];
+ AdaptedValues& adaptedValuesFather = adaptationMap[father];
+
+ adaptedValues.subdomain = adaptedValuesFather.subdomain;
+ adaptedValues.totalConcentration_[wCompIdx]
+ = adaptedValuesFather.totalConcentration_[wCompIdx] / adaptedValuesFather.count;
+ adaptedValues.totalConcentration_[nCompIdx]
+ = adaptedValuesFather.totalConcentration_[nCompIdx] / adaptedValuesFather.count;
+
+ // Introduce a hydrostatic pressure distribution inside the father cell
+ Scalar gTimesHeight = problem.gravity()
+ * (son.geometry().center() - father.geometry().center());
+ gTimesHeight *= (adaptedValuesFather.totalConcentration_[wCompIdx]+adaptedValuesFather.totalConcentration_[nCompIdx])/ problem.spatialParams().porosity(son);
+// int globalIdxSon = problem.variables().index(son);
+
+
+ // recalculate new Primary variable and store pc (i.e. other pressure)
+ if(pressureType == wPhaseIdx)
+ {
+ // recompute pressure and pc
+ Scalar pressure = adaptedValuesFather.pressure_[wPhaseIdx] / adaptedValuesFather.count;
+ Scalar pc = adaptedValuesFather.pressure_[nPhaseIdx] / adaptedValuesFather.count
+ - pressure;
+
+ adaptedValues.pressure_[wPhaseIdx]
+ = pressure + gTimesHeight ;
+ adaptedValues.pressure_[nPhaseIdx]
+ = pc + adaptedValues.pressure_[wPhaseIdx];
+ }
+ else
+ {
+ // recompute pressure and pc
+ Scalar pressure = adaptedValuesFather.pressure_[nPhaseIdx] / adaptedValuesFather.count;
+ Scalar pc = pressure
+ - adaptedValuesFather.pressure_[wPhaseIdx] / adaptedValuesFather.count;
+
+ adaptedValues.pressure_[nPhaseIdx]
+ = pressure + gTimesHeight ;
+ adaptedValues.pressure_[wPhaseIdx]
+ = adaptedValues.pressure_[nPhaseIdx] - pc;
+ }
+ // copy volume derivatives and compressibility, and indicate that cell was refined
+ adaptedValues.volumeDerivatives_[wCompIdx] = adaptedValuesFather.volumeDerivatives_[wCompIdx];
+ adaptedValues.volumeDerivatives_[nCompIdx] = adaptedValuesFather.volumeDerivatives_[nCompIdx];
+ adaptedValues.volumeDerivatives_[2] = adaptedValuesFather.volumeDerivatives_[2];
+ adaptedValues.isRefined = true;
+ }
+};
+
+}
+#endif
diff --git a/dumux/decoupled/2p2c/fvpressure2p2cadaptive.hh b/dumux/decoupled/2p2c/fvpressure2p2cadaptive.hh
new file mode 100644
index 0000000000..256f7fb257
--- /dev/null
+++ b/dumux/decoupled/2p2c/fvpressure2p2cadaptive.hh
@@ -0,0 +1,1277 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2010 by Benjamin Faigle *
+ * Copyright (C) 2007-2009 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Markus Wolff *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+#ifndef DUMUX_FVPRESSURE2P2C_ADAPTIVE_HH
+#define DUMUX_FVPRESSURE2P2C_ADAPTIVE_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/decoupled/2p2c/fvpressure2p2c.hh>
+#include <dumux/common/math.hh>
+#include <dumux/io/vtkmultiwriter.hh>
+#include <dumux/decoupled/2p2c/2p2cadaptiveproperties.hh>
+
+// include pressure model from Markus
+#include <dumux/decoupled/common/fv/mpfa/fvmpfaproperties.hh>
+#include <dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2pfaboundpressure2padaptive.hh>
+
+/**
+ * @file
+ * @brief Finite Volume Diffusion Model
+ * @author Benjamin Faigle, 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 compressible
+ * system with two components. An IMPES-like method is used for the sequential
+ * solution of the problem. Diffusion is neglected, capillarity can be regarded.
+ * 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 pressure base class FVPressure assembles the matrix and right-hand-side vector and solves for the pressure vector,
+ * whereas this class provides the actual entries for the matrix and RHS vector.
+ * 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 FVPressure2P2CAdaptive
+: public FVPressure2P2C<TypeTag>
+{
+ //the model implementation
+ typedef typename GET_PROP_TYPE(TypeTag, PressureModel) Implementation;
+ typedef FVPressure2P2C<TypeTag> ParentType;
+ typedef FVPressure<TypeTag> BaseType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+
+ typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+ typedef typename SpatialParams::MaterialLaw MaterialLaw;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+
+ typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
+ 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,
+ contiWEqIdx = Indices::contiWEqIdx, contiNEqIdx = Indices::contiNEqIdx
+ };
+ enum
+ {
+ rhs = BaseType::rhs, matrix = BaseType::matrix,
+ };
+
+ // 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::Intersection Intersection;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ // convenience shortcuts for Vectors/Matrices
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar,dim+1> TransmissivityMatrix;
+ typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
+ typedef Dune::FieldVector<Scalar, GET_PROP_VALUE(TypeTag, NumPhases)> PhaseVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+
+ // the typenames used for the stiffness matrix and solution vector
+ typedef typename GET_PROP_TYPE(TypeTag, PressureCoefficientMatrix) Matrix;
+ typedef typename GET_PROP_TYPE(TypeTag, PressureRHSVector) RHSVector;
+
+protected:
+ Problem& problem()
+ {
+ return this->problem_;
+ }
+ const Problem& problem() const
+ {
+ return this->problem_;
+ }
+
+public:
+ //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);
+ void getMpfaFlux(const IntersectionIterator&, const CellData&);
+
+ // mpfa transmissibilities
+ int computeTransmissibilities(const IntersectionIterator&,
+ IntersectionIterator&,
+ TransmissivityMatrix&,
+ TransmissivityMatrix&,
+ GlobalPosition&,
+ int&);
+
+ //! Adapt primary variables vector after adapting the grid
+ void adaptPressure()
+ {
+ int gridSize = problem().gridView().size(0);
+ this->pressure().resize(gridSize);
+
+ for(int i=0; i< gridSize; i++)
+ {
+ this->pressure()[i]
+ = problem().variables().cellData(i).pressure(this->pressureType);
+ }
+ }
+
+ //! Constructs a FVPressure2P2CAdaptive object
+ /**
+ * \param problem a problem class object
+ */
+ FVPressure2P2CAdaptive(Problem& problem) : FVPressure2P2C<TypeTag>(problem),
+ problem_(problem), pressureModelAdaptive2p_(0)
+ {
+ enableVolumeIntegral = GET_PARAM_FROM_GROUP(TypeTag,bool, Impet, EnableVolumeIntegral);
+ enableMPFA= GET_PARAM_FROM_GROUP(TypeTag,bool, GridAdapt, EnableMultiPointFluxApproximation);
+
+ enableSecondHalfEdge = GET_PARAM_FROM_GROUP(TypeTag,bool, GridAdapt, EnableSecondHalfEdge);
+ // prepare rotation matrix: R_ is initialized as zero matrix
+ R_=0.;
+ // evaluate matrix R
+ if (dim == 2)
+ {
+ R_[0][1] = 1;
+ R_[1][0] = -1;
+ }
+ }
+
+private:
+ Problem& problem_;
+ //! Returns the implementation of the problem (i.e. static polymorphism)
+ Implementation &asImp_()
+ { return *static_cast<Implementation *>(this);}
+
+ //! \copydoc Dumux::IMPETProblem::asImp_()
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation *>(this);}
+
+protected:
+ // mpfa transmissibilities from mpfal2pfa
+ int transmissibilityAdapter_(const IntersectionIterator&,
+ const IntersectionIterator&,
+ IntersectionIterator&,
+ GlobalPosition&,
+ TransmissivityMatrix&);
+
+ // Matrix for vector rotation of mpfa
+ // introduce matrix R for vector rotation and R is initialized as zero matrix
+ DimMatrix R_;
+ bool enableVolumeIntegral;
+ bool enableMPFA;
+ bool enableSecondHalfEdge;
+ FVMPFAL2PFABoundPressure2PAdaptive<TypeTag>* pressureModelAdaptive2p_;
+};
+
+
+//! initializes the matrix to store the system of equations
+template<class TypeTag>
+void FVPressure2P2CAdaptive<TypeTag>::initializeMatrix()
+{
+ int gridSize_ = problem().gridView().size(0);
+ // update RHS vector, matrix
+ this->A_.setSize (gridSize_,gridSize_); //
+ this->f_.resize(gridSize_);
+
+ // 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);
+ CellData& cellDataI = problem().variables().cellData(globalIdxI);
+
+ // initialize row size
+ int rowSize = 1;
+
+ // set perimeter to zero
+ cellDataI.perimeter() = 0;
+
+ // prepare storage for all found 3rd cells
+ std::vector<int> foundAdditionals;
+
+ int numberOfIntersections = 0;
+ // run through all intersections with neighbors
+ IntersectionIterator isItEnd = problem().gridView().template iend(*eIt);
+ for (IntersectionIterator isIt = problem().gridView().template ibegin(*eIt); isIt != isItEnd; ++isIt)
+ {
+ cellDataI.perimeter()
+ += isIt->geometry().volume();
+ numberOfIntersections++;
+ if (isIt->neighbor())
+ {
+ rowSize++;
+
+ // if mpfa is used, more entries might be needed if both halfedges are regarded
+ if (enableMPFA && enableSecondHalfEdge && isIt->outside()->level()!=eIt.level())
+ {
+ GlobalPosition globalPos3(0.);
+ int globalIdx3=-1;
+ TransmissivityMatrix T(0.);
+ IntersectionIterator additionalIsIt = isIt;
+ TransmissivityMatrix additionalT(0.);
+ // compute Transmissibilities: also examines subcontrolvolume information
+ int halfedgesStored
+ = problem().variables().getMpfaData(*isIt, additionalIsIt, T, additionalT, globalPos3, globalIdx3);
+ if (halfedgesStored == 0)
+ halfedgesStored = problem().pressureModel().computeTransmissibilities(isIt,additionalIsIt, T,additionalT, globalPos3, globalIdx3 );
+
+ if(halfedgesStored == 2)
+ {
+ bool increaseRowSize = true;
+ //check if additional cell is ordinary neighbor of eIt
+ for (IntersectionIterator isIt2 = problem().gridView().template ibegin(*eIt); isIt2 != isItEnd; ++isIt2)
+ {
+ if(!isIt2->neighbor())
+ continue;
+ if(additionalIsIt->outside() == isIt2->outside() )
+ increaseRowSize = false;
+ }
+ //also check if additional cell was already used for another interaction triangle
+ for(int i =0; i<foundAdditionals.size(); i++)
+ if(foundAdditionals[i] == problem().variables().index(*additionalIsIt->outside()))
+ increaseRowSize = false;
+
+ if (increaseRowSize)
+ {
+ rowSize++;
+ foundAdditionals.push_back(problem().variables().index(*additionalIsIt->outside()));
+ }
+ }
+ }
+ }
+ }
+ cellDataI.fluxData().resize(numberOfIntersections);
+ this->A_.setrowsize(globalIdxI, rowSize);
+ }
+ this->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
+ this->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
+ this->A_.addindex(globalIdxI, globalIdxJ);
+
+ // if mpfa is used, more entries might be needed if both halfedges are regarded
+ if (enableMPFA && enableSecondHalfEdge && isIt->outside()->level()!=eIt.level())
+ {
+ GlobalPosition globalPos3(0.);
+ int globalIdx3=-1;
+ TransmissivityMatrix T(0.);
+ IntersectionIterator additionalIsIt = isIt;
+ TransmissivityMatrix additionalT(0.);
+ // compute Transmissibilities: also examines subcontrolvolume information
+ int halfedgesStored
+ = problem().variables().getMpfaData(*isIt, additionalIsIt, T, additionalT, globalPos3, globalIdx3);
+ if (halfedgesStored == 0)
+ halfedgesStored = problem().pressureModel().computeTransmissibilities(isIt,additionalIsIt, T,additionalT, globalPos3, globalIdx3 );
+ // add off diagonal index if 2 half-edges regarded
+ if(halfedgesStored == 2)
+ this->A_.addindex(globalIdxI, problem().variables().index(*additionalIsIt->outside()));
+ }
+ }
+ }
+ this->A_.endindices();
+
+ return;
+}
+
+//! function which assembles the system of equations to be solved
+/* This function assembles the Matrix and the RHS vectors to solve for
+ * a pressure field with an Finite-Volume Discretization in an implicit
+ * fasion. In the implementations to this base class, the methods
+ * getSource(), getStorage(), getFlux() and getFluxOnBoundary() have to
+ * be provided.
+ * \param first Flag if pressure field is unknown
+ */
+template<class TypeTag>
+void FVPressure2P2CAdaptive<TypeTag>::assemble(bool first)
+{
+ if(first)
+ {
+ BaseType::assemble(true);
+ return;
+ }
+
+ initializeMatrix();
+ // initialization: set matrix A_ to zero
+ this->A_ = 0;
+ this->f_ = 0;
+
+ ElementIterator eItEnd = problem().gridView().template end<0>();
+ for (ElementIterator eIt = problem().gridView().template begin<0>(); eIt != eItEnd; ++eIt)
+ {
+ // cell information
+ int globalIdxI = problem().variables().index(*eIt);
+ CellData& cellDataI = problem().variables().cellData(globalIdxI);
+
+ Dune::FieldVector<Scalar, 2> entries(0.);
+
+ /***** source term ***********/
+ problem().pressureModel().getSource(entries, *eIt, cellDataI, first);
+ this->f_[globalIdxI] += entries[rhs];
+
+ /***** flux term ***********/
+ // 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)
+ {
+ /************* handle interior face *****************/
+ if (isIt->neighbor())
+ {
+ ElementPointer elementNeighbor = isIt->outside();
+
+ int globalIdxJ = problem().variables().index(*elementNeighbor);
+
+ //check for hanging nodes
+ //take a hanging node never from the element with smaller level!
+ bool haveSameLevel = (eIt->level() == elementNeighbor->level());
+ //calculate only from one side, but add matrix entries for both sides
+ if (GET_PROP_VALUE(TypeTag, VisitFacesOnlyOnce) && (globalIdxI > globalIdxJ) && haveSameLevel)
+ continue;
+
+ entries = 0;
+ //check for hanging nodes
+ if(!haveSameLevel && enableMPFA)
+ {
+ problem().pressureModel().getMpfaFlux(isIt, cellDataI);
+ }
+ else
+ {
+ problem().pressureModel().getFlux(entries, *isIt, cellDataI, first);
+
+ //set right hand side
+ this->f_[globalIdxI] -= entries[rhs];
+
+ // set diagonal entry
+ this->A_[globalIdxI][globalIdxI] += entries[matrix];
+
+ // set off-diagonal entry
+ this->A_[globalIdxI][globalIdxJ] -= entries[matrix];
+
+ if (GET_PROP_VALUE(TypeTag, VisitFacesOnlyOnce))
+ {
+ this->f_[globalIdxJ] += entries[rhs];
+ this->A_[globalIdxJ][globalIdxJ] += entries[matrix];
+ this->A_[globalIdxJ][globalIdxI] -= entries[matrix];
+ }
+ }
+
+ } // end neighbor
+
+ /************* boundary face ************************/
+ else
+ {
+ entries = 0;
+ problem().pressureModel().getFluxOnBoundary(entries, *isIt, cellDataI, first);
+
+ //set right hand side
+ this->f_[globalIdxI] += entries[rhs];
+ // set diagonal entry
+ this->A_[globalIdxI][globalIdxI] += entries[matrix];
+ }
+ } //end interfaces loop
+// printmatrix(std::cout, this->A_, "global stiffness matrix", "row", 11, 3);
+
+ /***** storage term ***********/
+ entries = 0;
+ problem().pressureModel().getStorage(entries, *eIt, cellDataI, first);
+ this->f_[globalIdxI] += entries[rhs];
+// set diagonal entry
+ this->A_[globalIdxI][globalIdxI] += entries[matrix];
+ } // end grid traversal
+// printmatrix(std::cout, this->A_, "global stiffness matrix after assempling", "row", 11,3);
+// printvector(std::cout, this->f_, "right hand side", "row", 10);
+ return;
+}
+
+//! Get flux at an interface between two cells
+/** for first == true, the flux is calculated in traditional fractional-flow forn as in FVPressure2P.
+ * for first == false, the flux thorugh \f$ \gamma{ij} \f$ is calculated via a volume balance formulation
+ * \f[ - A_{\gamma_{ij}} \cdot \mathbf{K} \cdot \mathbf{u} \cdot (\mathbf{n}_{\gamma_{ij}} \cdot \mathbf{u})
+ \sum_{\alpha} \varrho_{\alpha} \lambda_{\alpha} \sum_{\kappa} \frac{\partial v_{t}}{\partial C^{\kappa}} X^{\kappa}_{\alpha}
+ \left( \frac{p_{\alpha,j}^t - p^{t}_{\alpha,i}}{\Delta x} + \varrho_{\alpha} \mathbf{g}\right)
+ + V_i \frac{A_{\gamma_{ij}}}{U_i} \mathbf{K}
+ \sum_{\alpha} \varrho_{\alpha} \lambda_{\alpha} \sum_{\kappa}
+\frac{\frac{\partial v_{t,j}}{\partial C^{\kappa}_j}-\frac{\partial v_{t,i}}{\partial C^{\kappa}_i}}{\Delta x} X^{\kappa}_{\alpha}
+ \left( \frac{p_{\alpha,j}^t - p^{t}_{\alpha,i}}{\Delta x} + \varrho_{\alpha} \mathbf{g}\right) \f]
+ * This includes a boundary integral and a volume integral, because
+ * \f$ \frac{\partial v_{t,i}}{\partial C^{\kappa}_i} \f$ is not constant.
+ * Here, \f$ \mathbf{u} \f$ is the normalized vector connecting the cell centers, and \f$ \mathbf{n}_{\gamma_{ij}} \f$
+ * represents the normal of the face \f$ \gamma{ij} \f$.
+ * \param entries The Matrix and RHS entries
+ * \param intersection Intersection between cell I and J
+ * \param cellDataI Data of cell I
+ * \param first Flag if pressure field is unknown
+ */
+template<class TypeTag>
+void FVPressure2P2CAdaptive<TypeTag>::getMpfaFlux(const IntersectionIterator& intersectionIterator,
+ const CellData& cellDataI)
+{
+ // acess Cell I
+ ElementPointer elementPointerI = intersectionIterator->inside();
+ int globalIdxI = problem().variables().index(*elementPointerI);
+
+ // get global coordinate of cell center
+ const GlobalPosition& globalPos = elementPointerI->geometry().center();
+
+ // cell volume & perimeter, assume linear map here
+ Scalar volume = elementPointerI->geometry().volume();
+ Scalar perimeter = cellDataI.perimeter();
+
+ const GlobalPosition& gravity_ = problem().gravity();
+
+ // get absolute permeability
+ DimMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(*elementPointerI));
+
+ // access neighbor
+ ElementPointer neighborPointer = intersectionIterator->outside();
+ int globalIdxJ = problem().variables().index(*neighborPointer);
+ CellData& cellDataJ = problem().variables().cellData(globalIdxJ);
+
+ // gemotry info of neighbor
+ const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
+
+ // distance vector between barycenters
+ GlobalPosition distVec = globalPosNeighbor - globalPos;
+
+ // compute distance between cell centers
+ Scalar dist = distVec.two_norm();
+
+ GlobalPosition unitDistVec(distVec);
+ unitDistVec /= dist;
+
+ DimMatrix permeabilityJ
+ = problem().spatialParams().intrinsicPermeability(*neighborPointer);
+
+ // compute vectorized permeabilities
+ DimMatrix meanPermeability(0);
+ Dumux::harmonicMeanMatrix(meanPermeability, permeabilityI, permeabilityJ);
+
+ Dune::FieldVector<Scalar, dim> permeability(0);
+ meanPermeability.mv(unitDistVec, permeability);
+
+ // get average density for gravity flux
+ Scalar rhoMeanW = 0.5 * (cellDataI.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx));
+ Scalar rhoMeanNW = 0.5 * (cellDataI.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx));
+
+ // reset potential gradients
+ Scalar potentialW = 0;
+ Scalar potentialNW = 0;
+
+ // determine volume derivatives in neighbor
+ if (!cellDataJ.hasVolumeDerivatives())
+ this->volumeDerivatives(globalPosNeighbor, *neighborPointer);
+
+ Scalar dv_dC1 = (cellDataJ.dv(wPhaseIdx)
+ + cellDataI.dv(wPhaseIdx)) / 2; // dV/dm1= dv/dC^1
+ Scalar dv_dC2 = (cellDataJ.dv(nPhaseIdx)
+ + cellDataI.dv(nPhaseIdx)) / 2;
+
+ Scalar graddv_dC1 = (cellDataJ.dv(wPhaseIdx)
+ - cellDataI.dv(wPhaseIdx)) / dist;
+ Scalar graddv_dC2 = (cellDataJ.dv(nPhaseIdx)
+ - cellDataI.dv(nPhaseIdx)) / 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
+ Scalar densityW = rhoMeanW;
+ Scalar densityNW = rhoMeanNW;
+
+ // Prepare MPFA
+ /** get geometrical Info, transmissibility matrix */
+ GlobalPosition globalPos3(0.);
+ int globalIdx3=-1;
+ TransmissivityMatrix T(0.);
+ // prepare second half-edge
+ IntersectionIterator additionalIsIt = intersectionIterator;
+ TransmissivityMatrix additionalT(0.);
+
+ int halfedgesStored = problem().variables().getMpfaData(*intersectionIterator,
+ additionalIsIt, T, additionalT,
+ globalPos3, globalIdx3 );
+ if (halfedgesStored == 0)
+ halfedgesStored = problem().pressureModel().computeTransmissibilities(intersectionIterator,
+ additionalIsIt, T, additionalT,
+ globalPos3, globalIdx3 );
+
+ if(!halfedgesStored)
+ Dune::dgrave << "something went wrong getting mpfa data on cell " << globalIdxI << std::endl;
+
+ // shortcurts mpfa case
+ CellData& cellData3 = problem().variables().cellData(globalIdx3);
+ Scalar temp1 = globalPos * gravity_;
+ Scalar temp2 = globalPosNeighbor * gravity_;
+ Scalar temp3 = globalPos3 * gravity_;
+
+ potentialW = (cellDataI.pressure(wPhaseIdx)-temp1*densityW) * T[2]
+ + (cellDataJ.pressure(wPhaseIdx)-temp2*densityW) * T[0]
+ + (cellData3.pressure(wPhaseIdx)-temp3*densityW) * T[1];
+ potentialNW = (cellDataI.pressure(nPhaseIdx)-temp1*densityNW) * T[2]
+ + (cellDataJ.pressure(nPhaseIdx)-temp2*densityNW) * T[0]
+ + (cellData3.pressure(nPhaseIdx)-temp3*densityNW) * T[1];
+ // regard second half edge, if there is one
+ if(halfedgesStored == 2)
+ {
+ int AdditionalIdx = problem().variables().index(*(additionalIsIt->outside()));
+ CellData& cellDataAdditional = problem().variables().cellData(AdditionalIdx);
+ potentialW += (cellDataI.pressure(wPhaseIdx)-temp1*densityW) * additionalT[2]
+ + (cellDataJ.pressure(wPhaseIdx)-temp2*densityW) * additionalT[0]
+ + (cellDataAdditional.pressure(wPhaseIdx)
+ -(additionalIsIt->outside()->geometry().center()*gravity_)
+ *densityW) * additionalT[1];
+ potentialNW += (cellDataI.pressure(nPhaseIdx)-temp1*densityNW) * additionalT[2]
+ + (cellDataJ.pressure(nPhaseIdx)-temp2*densityNW) * additionalT[0]
+ + (cellDataAdditional.pressure(nPhaseIdx)
+ -(additionalIsIt->outside()->geometry().center()*gravity_)
+ *densityNW) * additionalT[1];
+ }
+
+
+ // initialize convenience shortcuts
+ Scalar lambdaW(0.), lambdaN(0.);
+ 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) line 3 analogon dV_w
+
+
+ //do the upwinding of the mobility depending on the phase potentials
+ if (potentialW >= 0.)
+ {
+ dV_w = (dv_dC1 * cellDataI.massFraction(wPhaseIdx, wCompIdx)
+ + dv_dC2 * cellDataI.massFraction(wPhaseIdx, nCompIdx));
+ lambdaW = cellDataI.mobility(wPhaseIdx);
+ gV_w = (graddv_dC1 * cellDataI.massFraction(wPhaseIdx, wCompIdx)
+ + graddv_dC2 * cellDataI.massFraction(wPhaseIdx, nCompIdx));
+ dV_w *= cellDataI.density(wPhaseIdx);
+ gV_w *= cellDataI.density(wPhaseIdx);
+ }
+ else
+ {
+ dV_w = (dv_dC1 * cellDataJ.massFraction(wPhaseIdx, wCompIdx)
+ + dv_dC2 * cellDataJ.massFraction(wPhaseIdx, nCompIdx));
+ lambdaW = cellDataJ.mobility(wPhaseIdx);
+ gV_w = (graddv_dC1 * cellDataJ.massFraction(wPhaseIdx, wCompIdx)
+ + graddv_dC2 * cellDataJ.massFraction(wPhaseIdx, nCompIdx));
+ dV_w *= cellDataJ.density(wPhaseIdx);
+ gV_w *= cellDataJ.density(wPhaseIdx);
+ }
+ if (potentialNW >= 0.)
+ {
+ dV_n = (dv_dC1 * cellDataI.massFraction(nPhaseIdx, wCompIdx)
+ + dv_dC2 * cellDataI.massFraction(nPhaseIdx, nCompIdx));
+ lambdaN = cellDataI.mobility(nPhaseIdx);
+ gV_n = (graddv_dC1 * cellDataI.massFraction(nPhaseIdx, wCompIdx)
+ + graddv_dC2 * cellDataI.massFraction(nPhaseIdx, nCompIdx));
+ dV_n *= cellDataI.density(nPhaseIdx);
+ gV_n *= cellDataI.density(nPhaseIdx);
+ }
+ else
+ {
+ dV_n = (dv_dC1 * cellDataJ.massFraction(nPhaseIdx, wCompIdx)
+ + dv_dC2 * cellDataJ.massFraction(nPhaseIdx, nCompIdx));
+ lambdaN = cellDataJ.mobility(nPhaseIdx);
+ gV_n = (graddv_dC1 * cellDataJ.massFraction(nPhaseIdx, wCompIdx)
+ + graddv_dC2 * cellDataJ.massFraction(nPhaseIdx, nCompIdx));
+ dV_n *= cellDataJ.density(nPhaseIdx);
+ gV_n *= cellDataJ.density(nPhaseIdx);
+ }
+
+ /** compute matrix entry: advective fluxes */
+ /* extend T with other matrix entries and assemble to A_ */
+ this->A_[globalIdxI][globalIdxJ] += (lambdaW * dV_w + lambdaN * dV_n) * T[0];
+ this->A_[globalIdxI][globalIdx3] += (lambdaW * dV_w + lambdaN * dV_n) * T[1];
+ this->A_[globalIdxI][globalIdxI] += (lambdaW * dV_w + lambdaN * dV_n) * T[2];
+
+ // add gravity to RHS vector
+ this->f_[globalIdxI] += (densityW * lambdaW * dV_w + densityNW * lambdaN * dV_n) * temp2 * T[0];
+ this->f_[globalIdxI] += (densityW * lambdaW * dV_w + densityNW * lambdaN * dV_n) * temp3 * T[1];
+ this->f_[globalIdxI] += (densityW * lambdaW * dV_w + densityNW * lambdaN * dV_n) * temp1 * T[2];
+
+ // weithing accounts for the fraction of the subcontrol volume
+ Scalar weightingFactor = volume / perimeter; // transforms flux through area A -> V * A/perimeter
+ if(enableVolumeIntegral) // switch off volume integral for mpfa case
+ {
+ // correct for area integral
+ this->A_[globalIdxI][globalIdxJ] -= weightingFactor * (lambdaW * gV_w + lambdaN * gV_n)*T[0];
+ this->A_[globalIdxI][globalIdx3] -= weightingFactor * (lambdaW * gV_w + lambdaN * gV_n)*T[1];
+ this->A_[globalIdxI][globalIdxI] -= weightingFactor * (lambdaW * gV_w + lambdaN * gV_n)*T[2];
+
+ // add gravity to RHS vector
+ this->f_[globalIdxI] -= weightingFactor * (densityW * lambdaW * gV_w + densityNW * lambdaN * gV_n) * temp2 * T[0];
+ this->f_[globalIdxI] -= weightingFactor * (densityW * lambdaW * gV_w + densityNW * lambdaN * gV_n) * temp3 * T[1];
+ this->f_[globalIdxI] -= weightingFactor * (densityW * lambdaW * gV_w + densityNW * lambdaN * gV_n) * temp1 * T[2];
+ }
+
+ // capillary pressure flux
+ Scalar pcGradient = cellDataI.capillaryPressure() * T[2]
+ + cellDataJ.capillaryPressure() * T[0]
+ + cellData3.capillaryPressure() * T[1];
+
+ if (this->pressureType == pw)
+ pcGradient *= + lambdaN * dV_n - enableVolumeIntegral * weightingFactor * lambdaN * gV_n;
+ else if (this->pressureType == pn)
+ pcGradient *= - lambdaW * dV_w + enableVolumeIntegral * weightingFactor * lambdaW * gV_w;
+
+ this->f_[globalIdxI] += pcGradient;
+
+ // include second half-edge
+ if(halfedgesStored == 2)
+ {
+ int AdditionalIdx = problem().variables().index(*(additionalIsIt->outside()));
+ CellData& cellDataAdditional = problem().variables().cellData(AdditionalIdx);
+
+ /* extend T with other matrix entries and assemble to A_ */
+ this->A_[globalIdxI][globalIdxJ] += (lambdaW * dV_w + lambdaN * dV_n) * additionalT[0];
+ this->A_[globalIdxI][AdditionalIdx] += (lambdaW * dV_w + lambdaN * dV_n) * additionalT[1];
+ this->A_[globalIdxI][globalIdxI] += (lambdaW * dV_w + lambdaN * dV_n) * additionalT[2];
+
+ // add gravity to RHS vector
+ this->f_[globalIdxI] += (densityW * lambdaW * dV_w + densityNW * lambdaN * dV_n) * temp2 * additionalT[0];
+ this->f_[globalIdxI] += (densityW * lambdaW * dV_w + densityNW * lambdaN * dV_n)
+ * (additionalIsIt->outside()->geometry().center()*gravity_) * additionalT[1];
+ this->f_[globalIdxI] += (densityW * lambdaW * dV_w + densityNW * lambdaN * dV_n) * temp1 * additionalT[2];
+
+ if(enableVolumeIntegral) // switch off volume integral for mpfa case
+ {
+ // correct for area integral
+ this->A_[globalIdxI][globalIdxJ] -= weightingFactor * (lambdaW * gV_w + lambdaN * gV_n)*additionalT[0];
+ this->A_[globalIdxI][AdditionalIdx] -= weightingFactor * (lambdaW * gV_w + lambdaN * gV_n)*additionalT[1];
+ this->A_[globalIdxI][globalIdxI] -= weightingFactor * (lambdaW * gV_w + lambdaN * gV_n)*additionalT[2];
+
+ // add gravity to RHS vector
+ this->f_[globalIdxI] -= weightingFactor * (densityW * lambdaW * gV_w + densityNW * lambdaN * gV_n) * temp2 * additionalT[0];
+ this->f_[globalIdxI] -= weightingFactor * (densityW * lambdaW * gV_w + densityNW * lambdaN * gV_n)
+ * (additionalIsIt->outside()->geometry().center()*gravity_) * additionalT[1];
+ this->f_[globalIdxI] -= weightingFactor * (densityW * lambdaW * gV_w + densityNW * lambdaN * gV_n) * temp1 * additionalT[2];
+ }
+
+ // capillary pressure flux
+ pcGradient = cellDataI.capillaryPressure() * additionalT[2]
+ + cellDataJ.capillaryPressure() * additionalT[0]
+ + cellDataAdditional.capillaryPressure() * additionalT[1];
+ if (this->pressureType == pw)
+ pcGradient *= + lambdaN * dV_n - enableVolumeIntegral * weightingFactor * lambdaN * gV_n;
+ else if (this->pressureType == pn)
+ pcGradient *= - lambdaW * dV_w + enableVolumeIntegral * weightingFactor * lambdaW * gV_w;
+
+ this->f_[globalIdxI] += pcGradient;
+ }
+}
+
+//! Computes the transmissibility coefficients for the MPFA-l method
+/* // Indices used in a interaction volume of the MPFA-l method
+ ___________________________________________________
+ | nux: cell geometry | nx: face normal |
+ | | =integrationOuterNormal|
+ | | |
+ | nextIt________________| |
+ | | nu4 3 |
+ | | \;.´ | |
+ | | .´ <--|nu3 |
+ | elementnumber | .´ ^n1 | |
+ | 1...........x_____|_____|____________|
+ | `. |nu5 !`. _nu7 | |
+ | `.`'´ ! `./` <--| |
+ lastIt | `._ ! `. nu2| |
+ \ | nu6/`. !nu1^ ` . | |
+ \ | `!---|------`2 |
+ `| | |
+ |\ isIt__________________|->n12 |
+ | \ | |
+ |_face14_________________|________________________|
+
+
+ _
+ /` is a normal directing upwards to the right,
+ while a normal directing downwards to the right looks like \;
+*
+*/
+template <class TypeTag>
+int FVPressure2P2CAdaptive<TypeTag>::computeTransmissibilities(const IntersectionIterator& isIt,
+ IntersectionIterator& additionalIntersectionIt,
+ TransmissivityMatrix& T,
+ TransmissivityMatrix& additionalT,
+ GlobalPosition& globalPos3,
+ int& globalIdx3)
+{
+ // get geometry information of cellI = cell1, cellJ = cell2
+ ElementPointer eIt = isIt->inside();
+ ElementPointer neighborPointer = isIt->outside();
+ GlobalPosition globalPos1 = eIt->geometry().center();
+ GlobalPosition globalPos2 = neighborPointer->geometry().center();
+ DimMatrix K1(problem().spatialParams().intrinsicPermeability(*eIt));
+ DimMatrix K2(problem().spatialParams().intrinsicPermeability(*neighborPointer));
+// int debugIdx = problem().variables().index(*isIt->outside());
+
+ /* 1. get geometrical information of interaction triangle */
+ // geometry and Data of face IJ in nomenclature of mpfa
+ GlobalPosition globalPosFace12 = isIt->geometry().center();
+ GlobalPosition integrationOuterNormaln12 = isIt->centerUnitOuterNormal();
+ integrationOuterNormaln12 *= isIt->geometry().volume() / 2.0; // TODO: 2.0 only in 2D
+
+
+ // nextIs points to next intersection
+ IntersectionIterator nextIs = isIt;
+ ++nextIs;
+ if (nextIs== problem().gridView().template iend(*eIt))
+ nextIs = problem().gridView().template ibegin(*eIt);
+
+ // get last intersection : --intersection does not exist
+ // paceingIt loops one IS bevore prevIs
+ IntersectionIterator prevIs = problem().gridView().template ibegin(*eIt);
+ IntersectionIterator paceingIt = prevIs;
+ for (++paceingIt; paceingIt != problem().gridView().template iend(*eIt); ++paceingIt)
+ {
+// if((paceingIt->neighbor()))
+// debugIdx = problem().variables().index(*paceingIt->outside());
+// if((prevIs->neighbor()))
+// debugIdx = problem().variables().index(*prevIs->outside());
+
+ if (!paceingIt->neighbor()) // continue if no neighbor found
+ ++prevIs; // we investigate next paceingIt -> prevIs is also increased
+ else if (paceingIt->outside() == isIt->outside()) // we already found prevIs
+ break;
+ else if (paceingIt == problem().gridView().template iend(*eIt))
+ prevIs = paceingIt; // this could only happen if isIt is begin, so prevIs has to be last.
+ else
+ ++prevIs; // we investigate next paceingIt -> prevIs is also increased
+ }
+
+// GlobalPosition prevIsCoords = prevIs->geometry().center();
+// GlobalPosition nextIsCoords = nextIs->geometry().center();
+
+ // search for face13, face23
+ IntersectionIterator face13 = isIt; // as long as isIt13=isIt, it is still not found!
+ IntersectionIterator face23 = isIt; // as long as face23=isIt, it is still not found!
+ // store other intersection for the other interaction region for the other half-edge
+
+ IntersectionIterator isIt23End = problem().gridView().iend(*neighborPointer);
+ for (IntersectionIterator isIt23 = problem().gridView().ibegin(*neighborPointer); isIt23 != isIt23End; ++isIt23)
+ {
+// if((isIt23->neighbor()))
+// debugIdx = problem().variables().index(*isIt23->outside());
+ // stop search if found
+ if( (face13->outside() != isIt->outside()))
+ break;
+
+ if(!(isIt23->neighbor()))
+ continue;
+
+ // either prevIs or nextIs is face13, it is that with common interface with cell2
+ // investigate if prevIs points to cell 3
+ if (prevIs->neighbor())
+ {
+// if((prevIs->neighbor()))
+// debugIdx = problem().variables().index(*prevIs->outside());
+
+ if (prevIs->outside() == isIt23->outside())
+ {
+ face23 = isIt23;
+ face13 = prevIs;
+ additionalIntersectionIt = nextIs;
+ }
+ }
+ // investigate if nextIs points to cell 3
+ if (nextIs->neighbor())
+ {
+// if((nextIs->neighbor()))
+// debugIdx = problem().variables().index(*nextIs->outside());
+
+ if (nextIs->outside() == isIt23->outside())
+ {
+ face23 = isIt23;
+ face13 = nextIs;
+ additionalIntersectionIt = prevIs;
+ }
+ }
+ }
+ if(face13->outside() == isIt->outside()) //means isIt13 not found yet
+ Dune::dgrave << "is 13 not found!!!" << std::endl;
+
+ // get information of cell3
+ globalPos3 = face13->outside()->geometry().center();
+ globalIdx3 = problem().variables().index(*(face13->outside()));
+ // get absolute permeability of neighbor cell 3
+ DimMatrix K3(problem().spatialParams().intrinsicPermeability(*(face13->outside())));
+
+
+ // get the intersection node /bar^{x_3} between 'isIt' and 'isIt13', denoted as 'corner123'
+ // initialization of corner123
+ GlobalPosition corner123(0.);
+
+ // get the global coordinate of corner123
+ for (int i = 0; i < isIt->geometry().corners(); ++i)
+ {
+ for (int j = 0; j < face13->geometry().corners(); ++j)
+ {
+ if (face13->geometry().corner(j) == isIt->geometry().corner(i))
+ {
+ corner123 = face13->geometry().corner(j);
+
+ // stop outer (i) and inner (j) loop
+ i = isIt->geometry().corners();
+ break;
+ }
+ }
+ }
+
+ /** 3. Calculate omega, chi for matrices **/
+ // center of face in global coordinates, i.e., the midpoint of edge 'isIt24'
+ GlobalPosition globalPosFace23 = face23->geometry().center();
+
+ // get face volume
+ Scalar face23vol = face23->geometry().volume();
+
+ // get outer normal vector scaled with half volume of face 'isIt24'
+ Dune::FieldVector<Scalar,dimWorld> integrationOuterNormaln23 = face23->centerUnitOuterNormal();
+ integrationOuterNormaln23 *= face23vol / 2.0; //TODO: 2.0 only for 2D
+
+ // compute normal vectors nu1-nu7 in triangle R for first half edge
+ GlobalPosition nu1(0);
+ R_.umv(globalPosFace12-globalPos2 ,nu1); //globalPos2 = globalPosNeighbor
+
+ GlobalPosition nu2(0);
+ R_.umv(globalPos2-globalPosFace23, nu2);
+
+ GlobalPosition nu3(0);
+ R_.umv(globalPosFace23-globalPos3, nu3);
+
+ GlobalPosition nu4(0);
+ R_.umv(globalPos3-corner123, nu4);
+
+ GlobalPosition nu5(0);
+ R_.umv(corner123-globalPos1, nu5); //globalPos1 = globalPos
+
+ GlobalPosition nu6(0);
+ R_.umv(globalPos1-globalPosFace12, nu6);
+
+ GlobalPosition nu7(0);
+ R_.umv(corner123-globalPos2, nu7);
+
+ // compute T, i.e., the area of quadrilateral made by normal vectors 'nu'
+ GlobalPosition Rnu2(0);
+ R_.umv(nu2, Rnu2);
+ Scalar T1 = nu1 * Rnu2;
+
+ GlobalPosition Rnu4(0);
+ R_.umv(nu4, Rnu4);
+ Scalar T2 = nu3 * Rnu4;
+
+ GlobalPosition Rnu6(0);
+ R_.umv(nu6, Rnu6);
+ Scalar T3 = nu5 * Rnu6;
+
+ // compute components needed for flux calculation, denoted as 'omega' and 'chi'
+ GlobalPosition K2nu1(0);
+ K2.umv(nu1, K2nu1); // permeabilityJ = K2 = perm of cell 2 around x_1
+ GlobalPosition K2nu2(0);
+ K2.umv(nu2, K2nu2);
+ GlobalPosition K3nu3(0);
+ K3.umv(nu3, K3nu3);
+ GlobalPosition K3nu4(0);
+ K3.umv(nu4, K3nu4);
+ GlobalPosition K1nu5(0);
+ K1.umv(nu5, K1nu5); // permeabilityI = K1 = perm of cell 1 around x_3
+ GlobalPosition K1nu6(0);
+ K1.umv(nu6, K1nu6);
+ GlobalPosition Rnu1(0);
+ R_.umv(nu1, Rnu1);
+
+ double omega111 = /*lambda2 */ (integrationOuterNormaln23 * K2nu1)/T1;
+ double omega112 = /*lambda2 */ (integrationOuterNormaln23 * K2nu2)/T1;
+ double omega211 = /*lambda2 */ (integrationOuterNormaln12 * K2nu1)/T1;
+ double omega212 = /*lambda2 */ (integrationOuterNormaln12 * K2nu2)/T1;
+ double omega123 = /*lambda3 */ (integrationOuterNormaln23 * K3nu3)/T2;
+ double omega124 = /*lambda3 */ (integrationOuterNormaln23 * K3nu4)/T2;
+ double omega235 = /*lambda1 */ (integrationOuterNormaln12 * K1nu5)/T3;
+ double omega236 = /*lambda1 */ (integrationOuterNormaln12 * K1nu6)/T3;
+ double chi711 = (nu7 * Rnu1)/T1;
+ double chi712 = (nu7 * Rnu2)/T1;
+
+ /** 4. Calculate A, B, C, D and solve for T **/
+ // compute transmissibility matrix T = CA^{-1}B+D
+ DimMatrix C(0), A(0);
+ Dune::FieldMatrix<Scalar,dim,2*dim-dim+1> D(0), B(0);
+
+ // evaluate matrix C, D, A, B
+ C[0][0] = -omega111;
+ C[0][1] = -omega112;
+ C[1][0] = -omega211;
+ C[1][1] = -omega212;
+
+ D[0][0] = omega111 + omega112;
+ D[1][0] = omega211 + omega212;
+
+ A[0][0] = omega111 - omega124 - omega123*chi711;
+ A[0][1] = omega112 - omega123*chi712;
+ A[1][0] = omega211 - omega236*chi711;
+ A[1][1] = omega212 - omega235 - omega236*chi712;
+
+ B[0][0] = omega111 + omega112 + omega123*(1.0 - chi711 - chi712);
+ B[0][1] = -omega123 - omega124;
+ B[1][0] = omega211 + omega212 + omega236*(1.0 - chi711 - chi712);
+ B[1][2] = -omega235 - omega236;
+
+ // compute T
+ A.invert();
+ D += B.leftmultiply(C.rightmultiply(A));
+ T = D[1];
+ if(!enableSecondHalfEdge )//or fabs(isIt->centerUnitOuterNormal()[0])<0.5) // [0]<0.5 => switch off vertical 2hes
+ {
+ T *= 2;
+ // set your map entry
+ problem().variables().storeMpfaData(*isIt, T, globalPos3, globalIdx3);
+ return 1; // indicates that only 1 halfedge was regarded
+ }
+ else
+ {
+ // derive additional T for second half edge
+
+ // get the intersection node /bar^{x_3} between 'isIt' and 'aditionalIntersection', denoted as 'corner1245'
+ // initialization of corner1245
+ GlobalPosition corner1245(0.);
+
+ // get the global coordinate of corner1245, which is not connected with face13
+ IntersectionIterator tempIntersection = face13;
+ bool corner1245found = false;
+ // ensure iterator increases over local end
+ if (tempIntersection== problem().gridView().template iend(*eIt))
+ tempIntersection = problem().gridView().template ibegin(*eIt);
+ while (!corner1245found)
+ {
+ ++tempIntersection;
+
+ // ensure iterator increases over local end
+ if (tempIntersection== problem().gridView().template iend(*eIt))
+ tempIntersection = problem().gridView().template ibegin(*eIt);
+ // enshure we do not arrive at isIt
+ if (tempIntersection == isIt)
+ continue;
+
+ // loop over both corners of is
+ for (int i = 0; i < isIt->geometry().corners(); ++i)
+ {
+ // test if a corner of additionalIntersectionIt also lies on is
+ for (int j = 0; j < tempIntersection->geometry().corners(); ++j)
+ {
+ if (tempIntersection->geometry().corner(j) == isIt->geometry().corner(i))
+ {
+ corner1245 = tempIntersection->geometry().corner(j);
+ additionalIntersectionIt = tempIntersection;
+
+ // stop outer (i) and inner (j) loop
+ i = isIt->geometry().corners();
+ // stop also Intersection loop
+ corner1245found = true;
+ break;
+ }
+ }
+ }
+ }
+ if (!additionalIntersectionIt->neighbor())// or (additionalIntersectionIt->outside().level() == isIt->inside().level()))
+ {
+ T *= 2;
+ // set your map entry
+ problem().variables().storeMpfaData(*isIt, T, globalPos3, globalIdx3);
+ return 1;
+ }
+
+ // use Markus mpfa-l implementation
+ this->transmissibilityAdapter_(isIt, face23, additionalIntersectionIt,
+ corner1245, additionalT);
+
+ // store transmissivity data for second half edge
+ problem().variables().storeMpfaData(*isIt, additionalIntersectionIt, T,additionalT, globalPos3, globalIdx3);
+ // return that information about two interaction volumes have been stored
+ return 2;
+ }
+}
+
+// mpfa transmissibilities from mpfal2pfa
+// Indices used in a interaction volume of the MPFA-o method
+// | | |
+// | 4-----------3-----------3 |
+// | | --> nu43 | nu34 <-- | |
+// | | |nu41 1|--> n43 ||nu32 |
+// | | v ^ |0 ^ v| |
+// |____________4__0__|n14__|__n23_|_1__2____________|
+// | | 1 0 | 0 | |
+// | | ^ |1 nu23 ^ | |
+// | | |nu14 0|--> n12 | | |
+// | | -->nu12 | nu21<-- | |
+// | J = 1-----------1-----------2 = I |
+// | | |
+template<class TypeTag>
+int FVPressure2P2CAdaptive<TypeTag>::transmissibilityAdapter_(const IntersectionIterator& isIt,
+ const IntersectionIterator& face23_2p2cnaming,
+ IntersectionIterator& additionalIntersectionIt,
+ GlobalPosition& corner1234,
+ TransmissivityMatrix& additionalT)
+{
+ /****** find all 4 faces *********/
+ // store additionalIntersection, which is face23 in 2p mpfa-l naming scheme
+ IntersectionIterator isIt23 = additionalIntersectionIt;
+ IntersectionIterator isIt14 = isIt; //has to be initialized with something
+
+ // search for 4th intersection that connects to corner1245, which is needed
+ // for markus implementation of the mpfa
+
+ // reuse corner1245found bool
+ bool face14found = false;
+ // repeat search procedure with isIt and face23 (in 2p2c naming)
+ IntersectionIterator tempIntersection = face23_2p2cnaming;
+
+ while (!face14found)
+ {
+ ++tempIntersection;
+
+ // ensure iterator increases over local end of neighbor J
+ if (tempIntersection== problem().gridView().template iend(*isIt->outside()))
+ tempIntersection = problem().gridView().template ibegin(*isIt->outside());
+
+ if(!tempIntersection->neighbor())
+ continue;
+
+ // enshure we hace not arrived at isIt (but seen from other side)
+ if (tempIntersection->outside() == isIt->inside())
+ continue; // restart loop to recheck after next increase for local end
+
+ // loop over both corners of is
+ for (int i = 0; i < isIt->geometry().corners(); ++i)
+ {
+ // test if a corner of additionalIntersectionIt also lies on is
+ for (int j = 0; j < tempIntersection->geometry().corners(); ++j)
+ {
+ if (tempIntersection->geometry().corner(j) == isIt->geometry().corner(i))
+ {
+// // test if this tempIntersection is better than additionalIntersectionIt
+// if (tempIntersection->outside().level() > additionalIntersectionIt->outside().level())
+// {
+// // select this as element for second half edge
+// additionalIntersectionIt = tempIntersection;
+// corner1245 = additionalIntersectionIt->geometry().corner(j);
+// }
+ isIt14 = tempIntersection;
+ // stop outer (i) and inner (j) loop
+ i = isIt->geometry().corners();
+ // stop also Intersection loop
+ face14found = true;
+ break;
+ }
+ }
+ }
+ }
+ /**** end find 4 faces **/
+
+
+ if(!pressureModelAdaptive2p_)
+ pressureModelAdaptive2p_= new FVMPFAL2PFABoundPressure2PAdaptive<TypeTag>(problem()) ;
+ // create Interaction Volume object
+ Dumux::FVMPFALInteractionVolume<TypeTag> interactionVolume;
+
+ interactionVolume.setCenterPosition(corner1234);
+
+ //*************** store pointer 1
+ interactionVolume.setSubVolumeElement(*isIt->outside(), 0);
+ interactionVolume.setIndexOnElement(isIt->indexInOutside(), 0, 0);
+ interactionVolume.setIndexOnElement(isIt14->indexInInside(), 0, 1);
+
+ // center of face in global coordinates, i.e., the midpoint of edge 'isIt12'
+ const GlobalPosition& globalPosFace12 = isIt->geometry().center();
+
+ // get face volume
+ Scalar faceVol12 = isIt->geometry().volume() / 2.0;
+
+ // get outer normal vector scaled with half volume of face 'isIt12'
+ Dune::FieldVector<Scalar, dimWorld> unitOuterNormal12 = isIt->centerUnitOuterNormal();
+ unitOuterNormal12 *=-1;
+
+ // center of face in global coordinates, i.e., the midpoint of edge 'isIt14'
+ const GlobalPosition& globalPosFace41 = isIt14->geometry().center();
+
+ // get face volume
+ Scalar faceVol41 = isIt14->geometry().volume() / 2.0;
+
+ // get outer normal vector scaled with half volume of face 'isIt14': for numbering of n see Aavatsmark, Eigestad
+ Dune::FieldVector<Scalar, dimWorld> unitOuterNormal14 = isIt14->centerUnitOuterNormal();
+
+ interactionVolume.setNormal(unitOuterNormal12, 0, 0);
+ interactionVolume.setNormal(unitOuterNormal14, 0, 1);
+ //get the normals of from cell 2 and 4
+ unitOuterNormal14 *= -1;
+ unitOuterNormal12 *= -1;
+ interactionVolume.setFaceArea(faceVol12, 0, 0);
+ interactionVolume.setFaceArea(faceVol41, 0, 1);
+ interactionVolume.setFacePosition(globalPosFace12, 0, 0);
+ interactionVolume.setFacePosition(globalPosFace41, 0, 1);
+
+ // access neighbor cell 2 of 'isIt12'
+ ElementPointer elementPointer2 = isIt->inside();
+
+ //**************** store pointer 2
+ interactionVolume.setSubVolumeElement(elementPointer2, 1);
+// interactionVolume.setIndexOnElement(isIt->indexInInside(), 1, 1);
+ interactionVolume.setNormal(unitOuterNormal12, 1, 1);
+ interactionVolume.setFaceArea(faceVol12, 1, 1);
+ interactionVolume.setFacePosition(globalPosFace12, 1, 1);
+
+
+ //**************** data for cell 4
+ ElementPointer elementPointer4 = isIt14->outside();
+
+ //store pointer 4
+ interactionVolume.setSubVolumeElement(elementPointer4, 3);
+// interactionVolume.setIndexOnElement(globalIdx4, 3, 0);
+
+ interactionVolume.setNormal(unitOuterNormal14, 3, 0);
+ interactionVolume.setFaceArea(faceVol41, 3, 0);
+ interactionVolume.setFacePosition(globalPosFace41, 3, 0);
+
+ //**************** data for cell 3
+ ElementPointer elementPointer3 = isIt23->outside();
+ //store pointer 3
+ interactionVolume.setSubVolumeElement(elementPointer3, 2);
+
+ GlobalPosition globalPosFace23 = isIt23->geometry().center();
+ Scalar faceVol23 = isIt23->geometry().volume() / 2.0;
+ // get outer normal vector scaled with half volume of face : for numbering of n see Aavatsmark, Eigestad
+ GlobalPosition unitOuterNormal23 = isIt23->centerUnitOuterNormal();
+
+ interactionVolume.setNormal(unitOuterNormal23, 1, 0);
+ unitOuterNormal23 *= -1;
+ interactionVolume.setNormal(unitOuterNormal23, 2, 1);
+ interactionVolume.setFaceArea(faceVol23, 1, 0);
+ interactionVolume.setFaceArea(faceVol23, 2, 1);
+ Scalar dummy=NAN;
+ interactionVolume.setFaceArea(dummy, 2, 0);
+ interactionVolume.setFaceArea(dummy, 3, 1);
+ interactionVolume.setFacePosition(globalPosFace23, 1, 0);
+ interactionVolume.setFacePosition(globalPosFace23, 2, 1);
+
+// int globalIdx1 = problem().variables().index(*isIt->outside());
+// int globalIdx2 = problem().variables().index(*elementPointer2);
+// int globalIdx3 = problem().variables().index(*elementPointer3);
+// int globalIdx4 = problem().variables().index(*isIt14->outside());
+
+ Dune::FieldVector<Scalar, dim> unity(1.);
+ std::vector<Dune::FieldVector<Scalar, dim> > lambda = {unity, unity, unity, unity};
+
+ Dune::FieldMatrix<Scalar,dim,2*dim-dim+1> T;
+ int triangleType = pressureModelAdaptive2p_->calculateTransmissibility(
+ T, interactionVolume, lambda,
+ 0, 1, 2, 3);
+
+ // 3.decide which triangle (which transmissibility coefficients) to use
+ if (triangleType == pressureModelAdaptive2p_->rightTriangle)
+ {
+ additionalIntersectionIt = isIt23;
+ // Translate flux from 2p mpfa-l local indexing to
+ // 2p2c naming with i and j (reverse direction)
+ // and adjust for the fact that it is right Triangle
+ additionalT[0] = -T[1][2];
+ additionalT[1] = -T[1][1];
+ additionalT[2] = -T[1][0];
+ }
+ else if (triangleType == pressureModelAdaptive2p_->leftTriangle)
+ {
+ additionalIntersectionIt = isIt14;
+ // Translate flux from 2p mpfa-l local indexing to
+ // 2p2c naming with i and j (reverse direction)
+ // and adjust for the fact that it is left Triangle
+ additionalT[0] = -T[1][0];
+ additionalT[1] = -T[1][1];
+ additionalT[2] = -T[1][2];
+
+ }
+ else
+ DUNE_THROW(Dune::MathError, "No transmissivity for second half edge found!");
+
+ return triangleType;
+}
+
+
+}//end namespace Dumux
+#endif
diff --git a/dumux/decoupled/2p2c/fvtransport2p2cadaptive.hh b/dumux/decoupled/2p2c/fvtransport2p2cadaptive.hh
new file mode 100644
index 0000000000..2753f9d195
--- /dev/null
+++ b/dumux/decoupled/2p2c/fvtransport2p2cadaptive.hh
@@ -0,0 +1,526 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2010 by Markus Wolff, Benjamin Faigle *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+#ifndef DUMUX_FVTRANSPORT2P2C_ADAPTIVE_HH
+#define DUMUX_FVTRANSPORT2P2C_ADAPTIVE_HH
+
+#include <dune/grid/common/gridenums.hh>
+#include <dumux/decoupled/2p2c/2p2cadaptiveproperties.hh>
+#include <dumux/decoupled/2p2c/fvtransport2p2c.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 FVTransport2P2CAdaptive : public FVTransport2P2C<TypeTag>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+
+ typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+ typedef typename SpatialParams::MaterialLaw MaterialLaw;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+
+ typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
+
+ typedef typename GET_PROP_TYPE(TypeTag, TransportSolutionType) TransportSolutionType;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld,
+ NumPhases = GET_PROP_VALUE(TypeTag, NumPhases)
+ };
+ enum
+ {
+ pw = Indices::pressureW,
+ pn = Indices::pressureNW,
+ Sw = Indices::saturationW,
+ Sn = Indices::saturationNW
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
+ wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx,
+ contiWEqIdx=Indices::contiWEqIdx, contiNEqIdx=Indices::contiNEqIdx
+ };
+
+ 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 typename GridView::Intersection Intersection;
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldMatrix<Scalar,dim,dim> DimMatrix;
+ typedef Dune::FieldVector<Scalar,dim+1> TransmissivityMatrix;
+ typedef Dune::FieldVector<Scalar, NumPhases> PhaseVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+
+ //! Acess function for the current problem
+ Problem& problem()
+ {return problem_;};
+
+public:
+ virtual void update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec, bool impes);
+
+ void getMpfaFlux(Dune::FieldVector<Scalar, 2>&, Dune::FieldVector<Scalar, 2>&,
+ const IntersectionIterator&, CellData&);
+
+ //! 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
+ */
+ FVTransport2P2CAdaptive(Problem& problem) : FVTransport2P2C<TypeTag>(problem),
+ problem_(problem)
+ {
+ enableMPFA = GET_PARAM_FROM_GROUP(TypeTag,bool,GridAdapt, EnableMultiPointFluxApproximation);
+ }
+
+ virtual ~FVTransport2P2CAdaptive()
+ {
+ }
+
+protected:
+ Problem& problem_;
+
+ bool enableMPFA;
+ static const int pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation); //!< gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
+};
+
+//! \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 \f$\mathrm{[s]}\f$
+ * \param[out] dt Time step size \f$\mathrm{[s]}\f$
+ * \param[out] updateVec Update vector, or update estimate for secants, resp. Here in \f$\mathrm{[kg/m^3]}\f$
+ * \param impet Flag that determines if it is a real impet step or an update estimate for volume derivatives
+ */
+template<class TypeTag>
+void FVTransport2P2CAdaptive<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec, bool impet = false)
+{
+ this->impet_ = impet;
+ // initialize dt very large
+ dt = 1E100;
+ this->averagedFaces_ = 0;
+
+ // resize update vector and set to zero
+ int size_ = problem_.gridView().size(0);
+ updateVec.resize(GET_PROP_VALUE(TypeTag, NumComponents));
+ updateVec[wCompIdx].resize(size_);
+ updateVec[nCompIdx].resize(size_);
+ updateVec[wCompIdx] = 0;
+ updateVec[nCompIdx] = 0;
+ //also resize PV vector if necessary
+ if(this->totalConcentration_.size() != size_)
+ {
+ this->totalConcentration_[wCompIdx].resize(size_);
+ this->totalConcentration_[nCompIdx].resize(size_);
+ // copy data //TODO: remove this, remove PM Transport Vector!!
+ // loop through all elements
+ for (int i = 0; i< problem().gridView().size(0); i++)
+ {
+ CellData& cellDataI = problem().variables().cellData(i);
+ for(int compIdx = 0; compIdx < GET_PROP_VALUE(TypeTag, NumComponents); compIdx++)
+ {
+ this->totalConcentration_[compIdx][i]
+ = cellDataI.totalConcentration(compIdx);
+ }
+ }
+ }
+
+ // Cell which restricts time step size
+ int restrictingCell = -1;
+
+ Dune::FieldVector<Scalar, 2> entries(0.), timestepFlux(0.);
+ // compute update vector
+ ElementIterator eItEnd = problem().gridView().template end<0> ();
+ for (ElementIterator eIt = problem().gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // get cell infos
+ int globalIdxI = problem().variables().index(*eIt);
+// const GlobalPosition globalPos = eIt->geometry().center();
+ CellData& cellDataI = problem().variables().cellData(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)
+ {
+ // handle interior face
+ if (isIt->neighbor())
+ {
+ if (enableMPFA && isIt->outside()->level()!=eIt.level())
+ getMpfaFlux(entries, timestepFlux, isIt, cellDataI);
+ else
+ getFlux(entries, timestepFlux, *isIt, cellDataI);
+ }
+
+ // Boundary Face
+ if (isIt->boundary())
+ {
+ getFluxOnBoundary(entries, timestepFlux, *isIt, cellDataI);
+ }
+
+ // add to update vector
+ updateVec[wCompIdx][globalIdxI] += entries[wCompIdx];
+ updateVec[nCompIdx][globalIdxI] += entries[nCompIdx];
+
+ // for time step calculation
+ sumfactorin += timestepFlux[0];
+ sumfactorout += timestepFlux[1];
+
+ }// end all intersections
+
+ /*********** Handle source term ***************/
+ PrimaryVariables q(NAN);
+ problem().source(q, *eIt);
+ updateVec[wCompIdx][globalIdxI] += q[contiWEqIdx];
+ updateVec[nCompIdx][globalIdxI] += q[contiNEqIdx];
+
+ // account for porosity
+ sumfactorin = std::max(sumfactorin,sumfactorout)
+ / problem().spatialParams().porosity(*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_PARAM_FROM_GROUP(TypeTag, Scalar,Impet, CFLFactor)<< std::endl;
+ if(this->averagedFaces_ != 0)
+ Dune::dinfo << " Averageing done for " << this->averagedFaces_ << " faces. "<< std::endl;
+ }
+ return;
+}
+
+template<class TypeTag>
+void FVTransport2P2CAdaptive<TypeTag>::getMpfaFlux(Dune::FieldVector<Scalar, 2>& fluxEntries, Dune::FieldVector<Scalar, 2>& timestepFlux,
+ const IntersectionIterator& intersectionIterator, CellData& cellDataI)
+{
+ fluxEntries = 0.;
+ timestepFlux = 0.;
+ // cell information
+ ElementPointer elementI= intersectionIterator->inside();
+ int globalIdxI = problem().variables().index(*elementI);
+
+ // get position
+ const GlobalPosition globalPos = elementI->geometry().center();
+ const GlobalPosition& gravity_ = problem().gravity();
+ // cell volume, assume linear map here
+ Scalar volume = elementI->geometry().volume();
+
+ // get values of cell I
+ Scalar pressI = problem().pressureModel().pressure(globalIdxI);
+ Scalar pcI = cellDataI.capillaryPressure();
+ DimMatrix K_I(problem().spatialParams().intrinsicPermeability(*elementI));
+
+ PhaseVector SmobI(0.);
+ SmobI[wPhaseIdx] = std::max((cellDataI.saturation(wPhaseIdx)
+ - problem().spatialParams().materialLawParams(*elementI).Swr())
+ , 1e-2);
+ SmobI[nPhaseIdx] = std::max((cellDataI.saturation(nPhaseIdx)
+ - problem().spatialParams().materialLawParams(*elementI).Snr())
+ , 1e-2);
+
+ Scalar densityWI (0.), densityNWI(0.);
+ densityWI= cellDataI.density(wPhaseIdx);
+ densityNWI = cellDataI.density(nPhaseIdx);
+
+ PhaseVector potential(0.);
+
+ // access neighbor
+ ElementPointer neighborPointer = intersectionIterator->outside();
+ int globalIdxJ = problem().variables().index(*neighborPointer);
+ CellData& cellDataJ = problem().variables().cellData(globalIdxJ);
+
+ // neighbor cell center in global coordinates
+ const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
+
+ // distance vector between barycenters
+ GlobalPosition distVec = globalPosNeighbor - globalPos;
+ // compute distance between cell centers
+ Scalar dist = distVec.two_norm();
+
+ GlobalPosition unitDistVec(distVec);
+ unitDistVec /= dist;
+
+ // phase densities in neighbor
+ Scalar densityWJ (0.), densityNWJ(0.);
+ densityWJ = cellDataJ.density(wPhaseIdx);
+ densityNWJ = cellDataJ.density(nPhaseIdx);
+
+ // average phase densities with central weighting
+ double densityW_mean = (densityWI + densityWJ) * 0.5;
+ double densityNW_mean = (densityNWI + densityNWJ) * 0.5;
+
+ double pressJ = problem().pressureModel().pressure(globalIdxJ);
+ Scalar pcJ = cellDataJ.capillaryPressure();
+
+ // determine potentials for upwind
+ /** get geometrical Info, transmissibility matrix */
+ GlobalPosition globalPos3(0.);
+ int globalIdx3=-1;
+ TransmissivityMatrix T(0.);
+ IntersectionIterator additionalIsIt = intersectionIterator;
+ TransmissivityMatrix additionalT(0.);
+
+ int halfedgesStored
+ = problem().variables().getMpfaData(*intersectionIterator, additionalIsIt, T, additionalT, globalPos3, globalIdx3);
+ if (halfedgesStored == 0)
+ halfedgesStored = problem().pressureModel().computeTransmissibilities(intersectionIterator,additionalIsIt, T,additionalT, globalPos3, globalIdx3 );
+
+ // acess cell 3 and prepare mpfa
+ Scalar press3 = problem().pressureModel().pressure(globalIdx3);
+ CellData& cellData3 = problem().variables().cellData(globalIdx3);
+ Scalar pc3 = cellData3.capillaryPressure();
+ Scalar temp1 = globalPos * gravity_;
+ Scalar temp2 = globalPosNeighbor * gravity_;
+ Scalar temp3 = globalPos3 * gravity_;
+ if(pressureType==pw)
+ {
+ potential[wPhaseIdx] += (pressI-temp1*densityW_mean) * T[2]
+ +(pressJ-temp2*densityW_mean) * T[0]
+ +(press3- temp3*densityW_mean) * T[1];
+ potential[nPhaseIdx] += (pressI+pcI-temp1*densityNW_mean) * T[2]
+ +(pressJ+pcJ-temp2*densityNW_mean) * T[0]
+ +(press3+pc3- temp3*densityNW_mean) * T[1];
+ // second half edge, if there is one
+ if(halfedgesStored == 2)
+ {
+ int AdditionalIdx = problem().variables().index(*(additionalIsIt->outside()));
+ CellData& cellDataAdditional = problem().variables().cellData(AdditionalIdx);
+ potential[wPhaseIdx] += (pressI-temp1*densityW_mean) * additionalT[2]
+ +(pressJ-temp2*densityW_mean) * additionalT[0]
+ +(problem().pressureModel().pressure(AdditionalIdx)
+ -(additionalIsIt->outside()->geometry().center()*gravity_*densityW_mean)
+ ) * additionalT[1];
+ potential[nPhaseIdx] += (pressI+pcI-temp1*densityNW_mean) * additionalT[2]
+ +(pressJ+pcJ-temp2*densityNW_mean) * additionalT[0]
+ +(problem().pressureModel().pressure(AdditionalIdx)
+ + cellDataAdditional.capillaryPressure()
+ -(additionalIsIt->outside()->geometry().center()*gravity_*densityNW_mean)
+ ) * additionalT[1];
+ }
+ }
+ else if(pressureType==pn)
+ {
+ potential[wPhaseIdx] += (pressI-pcI-temp1*densityW_mean) * T[2]
+ + (pressJ-pcJ-temp2*densityW_mean) * T[0]
+ + (press3-pc3- temp3*densityW_mean) * T[1];
+ potential[nPhaseIdx] += (pressI-temp1*densityNW_mean) * T[2]
+ + (pressJ-temp2*densityNW_mean) * T[0]
+ + (press3-temp3*densityNW_mean) * T[1];
+ // second half edge, if there is one
+ if(halfedgesStored == 2)
+ {
+ int AdditionalIdx = problem().variables().index(*(additionalIsIt->outside()));
+ CellData& cellDataAdditional = problem().variables().cellData(AdditionalIdx);
+
+ potential[wPhaseIdx] += (pressI-pcI-temp1*densityW_mean) * additionalT[2]
+ +(pressJ-pcJ-temp2*densityW_mean) * additionalT[0]
+ +(problem().pressureModel().pressure(AdditionalIdx)
+ - cellDataAdditional.capillaryPressure()
+ -(additionalIsIt->outside()->geometry().center()*gravity_*densityW_mean)
+ ) * additionalT[1];
+ potential[nPhaseIdx] += (pressI-temp1*densityNW_mean) * additionalT[2]
+ +(pressJ-temp2*densityNW_mean) * additionalT[0]
+ +(problem().pressureModel().pressure(AdditionalIdx)
+ -(additionalIsIt->outside()->geometry().center()*gravity_*densityNW_mean))
+ * additionalT[1];
+ }
+ }
+ //end of mpfa specific stuff
+
+ // determine upwinding direction, perform upwinding if possible
+ Dune::FieldVector<bool, NumPhases> doUpwinding(true);
+ PhaseVector lambda(0.);
+ for(int phaseIdx = 0; phaseIdx < NumPhases; phaseIdx++)
+ {
+ int contiEqIdx = 0;
+ if(phaseIdx == wPhaseIdx)
+ contiEqIdx = contiWEqIdx;
+ else
+ contiEqIdx = contiNEqIdx;
+
+ if(!this->impet_ or !this->restrictFluxInTransport_) // perform a strict uwpind scheme
+ {
+ if(potential[phaseIdx] > 0.)
+ {
+ lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
+ cellDataI.setUpwindCell(intersectionIterator->indexInInside(), contiEqIdx, true);
+ }
+ else if(potential[phaseIdx] < 0.)
+ {
+ lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
+ cellDataI.setUpwindCell(intersectionIterator->indexInInside(), contiEqIdx, false);
+ }
+ else
+ {
+ doUpwinding[phaseIdx] = false;
+ cellDataI.setUpwindCell(intersectionIterator->indexInInside(), contiEqIdx, false);
+ cellDataJ.setUpwindCell(intersectionIterator->indexInOutside(), contiEqIdx, false);
+ }
+ }
+ else // Transport after PE with check on flow direction
+ {
+ bool cellIwasUpwindCell;
+ //get the information from smaller (higher level) cell, as its IS is unique
+ if(elementI->level()>neighborPointer->level())
+ cellIwasUpwindCell = cellDataI.isUpwindCell(intersectionIterator->indexInInside(), contiEqIdx);
+ else // reverse neighbors information gathered
+ cellIwasUpwindCell = !cellDataJ.isUpwindCell(intersectionIterator->indexInOutside(), contiEqIdx);
+
+ if (potential[phaseIdx] > 0. && cellIwasUpwindCell)
+ lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
+ else if (potential[phaseIdx] < 0. && !cellIwasUpwindCell)
+ lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
+ // potential direction does not coincide with that of P.E.
+ else if(this->restrictFluxInTransport_ == 2) // perform central averageing for all direction changes
+ doUpwinding[phaseIdx] = false;
+ else // i.e. restrictFluxInTransport == 1
+ {
+ //check if harmonic weithing is necessary
+ if (!cellIwasUpwindCell && cellDataJ.mobility(phaseIdx) != 0.) // check if outflow induce neglected (i.e. mob=0) phase flux
+ lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
+ else if (cellIwasUpwindCell && cellDataI.mobility(phaseIdx) != 0.) // check if inflow induce neglected phase flux
+ lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
+ else
+ doUpwinding[phaseIdx] = false;
+ }
+
+ // do not perform upwinding if so desired
+ if(!doUpwinding[phaseIdx])
+ {
+ //a) no flux if there wont be any flux regardless how to average/upwind
+ if(cellDataI.mobility(phaseIdx)+cellDataJ.mobility(phaseIdx)==0.)
+ {
+ potential[phaseIdx] = 0;
+ continue;
+ }
+
+ //b) perform harmonic averageing
+ fluxEntries[wCompIdx] -= potential[phaseIdx] / volume
+ * harmonicMean(cellDataI.massFraction(phaseIdx, wCompIdx) * cellDataI.mobility(phaseIdx) * cellDataI.density(phaseIdx),
+ cellDataJ.massFraction(phaseIdx, wCompIdx) * cellDataJ.mobility(phaseIdx) * cellDataJ.density(phaseIdx));
+ fluxEntries[nCompIdx] -= potential[phaseIdx] / volume
+ * harmonicMean(cellDataI.massFraction(phaseIdx, nCompIdx) * cellDataI.mobility(phaseIdx) * cellDataI.density(phaseIdx),
+ cellDataJ.massFraction(phaseIdx, nCompIdx) * cellDataJ.mobility(phaseIdx) * cellDataJ.density(phaseIdx));
+ // c) timestep control
+ // for timestep control : influx
+ timestepFlux[0] += std::max(0.,
+ - potential[phaseIdx] / volume
+ * harmonicMean(cellDataI.mobility(phaseIdx),cellDataJ.mobility(phaseIdx)));
+ // outflux
+ timestepFlux[1] += std::max(0.,
+ potential[phaseIdx] / volume
+ * harmonicMean(cellDataI.mobility(phaseIdx),cellDataJ.mobility(phaseIdx))/SmobI[phaseIdx]);
+
+ //d) output (only for one side)
+ this->averagedFaces_++;
+ #if DUNE_MINIMAL_DEBUG_LEVEL < 3
+ // verbose (only for one side)
+ if(globalIdxI > globalIdxJ)
+ Dune::dinfo << "harmonicMean flux of phase" << phaseIdx <<" used from cell" << globalIdxI<< " into " << globalIdxJ
+ << " ; TE upwind I = "<< cellDataI.isUpwindCell(intersectionIterator->indexInInside(), contiEqIdx) << " but pot = "<< potential[phaseIdx] << " \n";
+ #endif
+
+ //e) stop further standard calculations
+ potential[phaseIdx] = 0;
+ }
+ }
+ }
+
+ // calculate and standardized velocity
+ double velocityJIw = std::max((-lambda[wPhaseIdx] * potential[wPhaseIdx]) / volume, 0.0);
+ double velocityIJw = std::max(( lambda[wPhaseIdx] * potential[wPhaseIdx]) / volume, 0.0);
+ double velocityJIn = std::max((-lambda[nPhaseIdx] * potential[nPhaseIdx]) / volume, 0.0);
+ double velocityIJn = std::max(( lambda[nPhaseIdx] * potential[nPhaseIdx]) / volume, 0.0);
+
+ // for timestep control
+ timestepFlux[0] += velocityJIw + velocityJIn;
+
+ double foutw = velocityIJw/SmobI[wPhaseIdx];
+ double foutn = velocityIJn/SmobI[nPhaseIdx];
+ if (std::isnan(foutw) || std::isinf(foutw) || foutw < 0) foutw = 0;
+ if (std::isnan(foutn) || std::isinf(foutn) || foutn < 0) foutn = 0;
+ timestepFlux[1] += foutw + foutn;
+
+ fluxEntries[wCompIdx] +=
+ velocityJIw * cellDataJ.massFraction(wPhaseIdx, wCompIdx) * densityWJ
+ - velocityIJw * cellDataI.massFraction(wPhaseIdx, wCompIdx) * densityWI
+ + velocityJIn * cellDataJ.massFraction(nPhaseIdx, wCompIdx) * densityNWJ
+ - velocityIJn * cellDataI.massFraction(nPhaseIdx, wCompIdx) * densityNWI;
+ fluxEntries[nCompIdx] +=
+ velocityJIw * cellDataJ.massFraction(wPhaseIdx, nCompIdx) * densityWJ
+ - velocityIJw * cellDataI.massFraction(wPhaseIdx, nCompIdx) * densityWI
+ + velocityJIn * cellDataJ.massFraction(nPhaseIdx, nCompIdx) * densityNWJ
+ - velocityIJn * cellDataI.massFraction(nPhaseIdx, nCompIdx) * densityNWI;
+ return;
+}
+
+}
+#endif
diff --git a/dumux/decoupled/2p2c/variableclass2p2cadaptive.hh b/dumux/decoupled/2p2c/variableclass2p2cadaptive.hh
new file mode 100644
index 0000000000..94a9594a79
--- /dev/null
+++ b/dumux/decoupled/2p2c/variableclass2p2cadaptive.hh
@@ -0,0 +1,301 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2012 by Benjamin Faigle *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * T1his 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. *
+ * *
+ * T1his program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+#ifndef DUMUX_VARIABLECLASS2P2C_ADAPTIVE_HH
+#define DUMUX_VARIABLECLASS2P2C_ADAPTIVE_HH
+
+
+// for parallelization
+#include <dumux/decoupled/common/variableclassadaptive.hh>
+
+/**
+ * @file
+ * @brief Base class holding the variables for sequential models.
+ * @author Markus Wolff
+ */
+
+namespace Dumux
+{
+/*!
+ * \ingroup IMPET
+ */
+//! Class holding additionally mpfa data of adaptive compositional models.
+/*!
+ * This class provides the possibility to store and load the transmissibilities (and associated infos)
+ * of the mpfa method per irregular face. While this class provides the storage container for both 2D
+ * and 3D implementation, only the 2D storage methods are included here.
+ * Note that according to the number of half-edges (sub-faces) regarded, one ore two transmissibility
+ * can be stored and loaded.
+ *
+ * @tparam TypeTag The Type Tag
+ */
+template<class TypeTag>
+class VariableClass2P2CAdaptive: public VariableClassAdaptive<TypeTag>
+{
+private:
+ typedef VariableClassAdaptive<TypeTag> ParentType;
+ typedef VariableClass<TypeTag> BaseType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
+ typedef typename CellData::AdaptedValues AdaptedValues;
+
+ typedef typename GridView::Grid Grid;
+ typedef typename Grid::LevelGridView LevelGridView;
+ typedef typename Grid::LocalIdSet::IdType IdType;
+ typedef typename LevelGridView::template Codim<0>::Iterator LevelIterator;
+ typedef typename GridView::Traits::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename GridView::Intersection Intersection;
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+ enum // for first and second half edge 2D / subvolume face 3D
+ {
+ first = 0, second = 1
+ };
+ // convenience shortcuts for Vectors/Matrices
+ typedef Dune::FieldVector<Scalar, GridView::dimensionworld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar,dim+1> T1ransmissivityMatrix;
+
+protected:
+ // in the 2D case, we need to store 1 additional cell. In 3D, we store
+ // dim-1=2 cells for one interaction volume, and 4 if two interaction volumes
+ // are regarded.
+ const static int storageRequirement = (dim-1)*(dim-1);
+ //! Storage object for data related to the MPFA method
+ /*
+ * This Struct stores the transmissibility coefficients
+ * for the two half-eges of an irregular interface (one
+ * near a hanging node) in an h-adaptive simulation.
+ */
+ struct mpfaData
+ {
+ T1ransmissivityMatrix T1_[2];
+ int globalIdx3_[storageRequirement];
+ GlobalPosition globalPos3_[storageRequirement];
+ std::vector<IntersectionIterator> secondHalfEdgeIntersection_;
+ bool hasSecondHalfEdge;
+
+ mpfaData()
+ {
+ hasSecondHalfEdge = false;
+ }
+
+ void setIntersection(IntersectionIterator& is23)
+ {
+ secondHalfEdgeIntersection_.push_back(is23);
+ hasSecondHalfEdge = true;
+ };
+ const IntersectionIterator& getIntersection()
+ {
+ return secondHalfEdgeIntersection_[0];
+ };
+ };
+ std::map<IdType, mpfaData> irregularInterfaceMap_;
+
+ const Grid& grid_;
+
+public:
+ //! 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 initialVel initial value for the velocity (only necessary if only transport part is solved)
+ */
+ VariableClass2P2CAdaptive(const GridView& gridView) :
+ ParentType(gridView), grid_(gridView.grid())
+ {}
+
+ //! Resizes decoupled variable vectors
+ /*! Method that change the size of the vectors for h-adaptive simulations, and clears recently
+ * stored transmissibility matrices for the newly adapted grid.
+ *
+ *\param size Size of the current (refined and coarsened) grid
+ */
+ void adaptVariableSize(int size)
+ {
+ BaseType::adaptVariableSize(size);
+ // clear mapper
+ irregularInterfaceMap_.clear();
+ }
+ /*!
+ * Reconstruct missing primary variables (where elements are created/deleted)
+ *
+ * To reconstruct the solution in father elements, problem properties might
+ * need to be accessed.
+ *
+ * @param problem The current problem
+ */
+ void reconstructPrimVars(Problem& problem)
+ {
+ ParentType::reconstructPrimVars(problem);
+
+ problem.pressureModel().adaptPressure();
+ }
+
+ void storeMpfaData(const typename GridView::Intersection & irregularIs,
+ const T1ransmissivityMatrix& T1,
+ const GlobalPosition& globalPos3,
+ const int& globalIdx3)
+ {
+ IdType intersectionID
+ = grid_.localIdSet().subId(*irregularIs.inside(),
+ irregularIs.indexInInside(), 1);
+ // mapping is only unique from smaller cell (if *inside and not *outside)
+ if (irregularIs.inside().level() < irregularIs.outside().level())
+ {
+ // IS is regarded from larger cell: get the unique number as seen from smaller
+ intersectionID
+ = grid_.localIdSet().subId(*irregularIs.outside(),
+ irregularIs.indexInOutside(), 1);
+
+ // store as if it was seen from smaller: change i & j
+ irregularInterfaceMap_[intersectionID].T1_[first][2] = - T1[0];
+ irregularInterfaceMap_[intersectionID].T1_[first][1] = - T1[1];
+ irregularInterfaceMap_[intersectionID].T1_[first][0] = - T1[2];
+ }
+ else // we look from smaller cell = unique interface
+ // proceed with numbering according to Aavatsmark, seen from cell i
+ irregularInterfaceMap_[intersectionID].T1_[first] = T1;
+
+ irregularInterfaceMap_[intersectionID].globalPos3_[0] = globalPos3;
+ irregularInterfaceMap_[intersectionID].globalIdx3_[0] = globalIdx3;
+ return;
+ }
+
+ void storeMpfaData(const typename GridView::Intersection & irregularIs,
+ IntersectionIterator& secondHalfEdgeIntersectionIt,
+ const T1ransmissivityMatrix& T1,
+ const T1ransmissivityMatrix& T11_secondHalfEdge,
+ const GlobalPosition& globalPos3,
+ const int& globalIdx3)
+ {
+ IdType intersectionID
+ = grid_.localIdSet().subId(*irregularIs.inside(),
+ irregularIs.indexInInside(), 1);
+
+ // mapping is only unique from smaller cell (if *inside and not *outside)
+ if (irregularIs.inside().level() < irregularIs.outside().level())
+ {
+ // IS is regarded from larger cell: get the unique number as seen from smaller
+ intersectionID
+ = grid_.localIdSet().subId(*irregularIs.outside(),
+ irregularIs.indexInOutside(), 1);
+
+ // store as if it was seen from smaller: change i & j
+ irregularInterfaceMap_[intersectionID].T1_[first][2] = - T1[0];
+ irregularInterfaceMap_[intersectionID].T1_[first][1] = - T1[1];
+ irregularInterfaceMap_[intersectionID].T1_[first][0] = - T1[2];
+
+ irregularInterfaceMap_[intersectionID].T1_[second][2] = - T11_secondHalfEdge[0];
+ irregularInterfaceMap_[intersectionID].T1_[second][1] = - T11_secondHalfEdge[1];
+ irregularInterfaceMap_[intersectionID].T1_[second][0] = - T11_secondHalfEdge[2];
+
+ }
+ else // we look from smaller cell = unique interface
+ {
+ // proceed with numbering according to Aavatsmark, seen from cell i
+ irregularInterfaceMap_[intersectionID].T1_[first] = T1;
+ irregularInterfaceMap_[intersectionID].T1_[second] = T11_secondHalfEdge;
+ }
+
+ irregularInterfaceMap_[intersectionID].globalPos3_[0] = globalPos3;
+ irregularInterfaceMap_[intersectionID].globalIdx3_[0] = globalIdx3;
+ // second half edge
+ irregularInterfaceMap_[intersectionID].setIntersection(secondHalfEdgeIntersectionIt);
+
+
+ return;
+ }
+
+ int getMpfaData(const Intersection& irregularIs,
+ IntersectionIterator& secondHalfEdgeIntersectionIt,
+ T1ransmissivityMatrix& T1,
+ T1ransmissivityMatrix& T1_secondHalfEdge,
+ GlobalPosition& globalPos3,
+ int& globalIdx3)
+ {
+ IdType intersectionID
+ = grid_.localIdSet().subId(*irregularIs.inside(),
+ irregularIs.indexInInside(), 1);
+ // mapping is only unique from smaller cell (if *inside and not *outside)
+ if (irregularIs.inside().level() < irregularIs.outside().level())
+ {
+ // IS is regarded from larger cell: get the unique number as seen from smaller
+ intersectionID
+ = grid_.localIdSet().subId(*irregularIs.outside(),
+ irregularIs.indexInOutside(), 1);
+
+ // check if T1ransmissibility matrix was stored for that IF
+ if (irregularInterfaceMap_.find(intersectionID) == irregularInterfaceMap_.end())
+ return 0; // no stored data!
+
+ // If data is stored, it is so as if the IF is regarded from smaller cell.
+ // since we are looking from larger cell, cells i & j have to be changed
+ // Additionally, flux points in opposite direction: - sign
+ T1[0] = -irregularInterfaceMap_[intersectionID].T1_[first][2];
+ T1[1] = -irregularInterfaceMap_[intersectionID].T1_[first][1];
+ T1[2] = -irregularInterfaceMap_[intersectionID].T1_[first][0];
+ globalPos3 = irregularInterfaceMap_[intersectionID].globalPos3_[0];
+ globalIdx3 = irregularInterfaceMap_[intersectionID].globalIdx3_[0];
+ //second half edge
+ if(irregularInterfaceMap_[intersectionID].hasSecondHalfEdge)
+ {
+ secondHalfEdgeIntersectionIt = irregularInterfaceMap_[intersectionID].getIntersection();
+ T1_secondHalfEdge[0] = -irregularInterfaceMap_[intersectionID].T1_[second][2];
+ T1_secondHalfEdge[1] = -irregularInterfaceMap_[intersectionID].T1_[second][1];
+ T1_secondHalfEdge[2] = -irregularInterfaceMap_[intersectionID].T1_[second][0];
+ // Dune::dinfo << "mpfa Info retrieved for isID " << intersectionID
+ // << "at coordinate " << irregularIs.geometry().center() << " from GlobalIdx " << this->index(*irregularIs.inside())<<std::endl;
+ return 2;
+ }
+ return 1;
+ }
+ // check if T1ransmissibility matrix was stored for that IF
+ if (irregularInterfaceMap_.find(intersectionID) == irregularInterfaceMap_.end())
+ return 0; // no stored data!
+
+ T1[0] = irregularInterfaceMap_[intersectionID].T1_[first][0];
+ T1[1] = irregularInterfaceMap_[intersectionID].T1_[first][1];
+ T1[2] = irregularInterfaceMap_[intersectionID].T1_[first][2];
+ globalPos3 = irregularInterfaceMap_[intersectionID].globalPos3_[0];
+ globalIdx3 = irregularInterfaceMap_[intersectionID].globalIdx3_[0];
+ //second half edge
+ if(irregularInterfaceMap_[intersectionID].hasSecondHalfEdge)
+ {
+ secondHalfEdgeIntersectionIt = irregularInterfaceMap_[intersectionID].getIntersection();
+ T1_secondHalfEdge = irregularInterfaceMap_[intersectionID].T1_[second];
+
+ // Dune::dinfo << "mpfa Info retrieved for isID " << intersectionID
+ // << "at coordinate " << irregularIs.geometry().center() << " from GlobalIdx " << this->index(*irregularIs.inside())<<std::endl;
+ return 2;
+ }
+ return 1;
+ }
+ //@}
+
+};
+}
+#endif
diff --git a/test/decoupled/2p2c/Makefile.am b/test/decoupled/2p2c/Makefile.am
index 24ef9a813d..f3cd2db21a 100644
--- a/test/decoupled/2p2c/Makefile.am
+++ b/test/decoupled/2p2c/Makefile.am
@@ -1,10 +1,10 @@
# programs just to build when "make check" is used
-check_PROGRAMS = test_dec2p2c test_multiphysics2p2c
+check_PROGRAMS = test_adaptive2p2c test_dec2p2c test_multiphysics2p2c
noinst_HEADERS = *.hh
+test_adaptive2p2c_SOURCES = test_adaptive2p2c.cc
test_dec2p2c_SOURCES = test_dec2p2c.cc
-
test_multiphysics2p2c_SOURCES = test_multiphysics2p2c.cc
EXTRA_DIST=*reference.vtu *.input CMakeLists.txt
diff --git a/test/decoupled/2p2c/test_adaptive2p2c-reference.vtu b/test/decoupled/2p2c/test_adaptive2p2c-reference.vtu
new file mode 100644
index 0000000000..6b2dcfc0d4
--- /dev/null
+++ b/test/decoupled/2p2c/test_adaptive2p2c-reference.vtu
@@ -0,0 +1,457 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+ <UnstructuredGrid>
+ <Piece NumberOfCells="163" NumberOfPoints="194">
+ <CellData Scalars="wetting pressure">
+ <DataArray type="Float32" Name="wetting pressure" NumberOfComponents="1" format="ascii">
+ 265462 270504 271922 270840 267760 262867 256119 247265 235777 220565 261915 264320
+ 264615 263029 259766 254916 248417 240093 229693 217053 260169 260246 259189 256896
+ 253361 248612 242540 234998 225848 215386 259060 257307 255000 252005 250427 248358
+ 245990 248042 246062 243529 241370 243782 240696 237590 235794 238714 231512 223317
+ 214405 258202 254922 251498 249784 247922 245956 247848 245929 243941 242236 244051
+ 241812 239524 238048 240228 237046 234414 233185 235692 231723 227956 226991 230527
+ 221504 213702 257407 252680 248163 245971 243997 241952 243999 242097 240378 238337
+ 240029 238536 236500 234711 236605 234287 231922 230483 232655 229376 225793 224045
+ 227752 219522 212978 256529 250229 244524 241931 239773 237380 239602 237816 236110
+ 233690 235461 234438 232661 230368 232036 230733 228662 225778 228492 225439 222119
+ 219780 222707 216830 212043 255388 247183 240129 234374 232722 230981 227368 229264
+ 229325 227255 223116 225405 224890 222290 218594 220880 217040 213369 210787 253607
+ 242987 234435 227673 222377 217907 213872 210643 208822 208901 250036 236722 226952
+ 219507 213705 209039 205314 202841 202376 205273
+ </DataArray>
+ <DataArray type="Float32" Name="nonwetting pressure" NumberOfComponents="1" format="ascii">
+ 265462 270504 271922 270840 267760 262867 256119 247265 235777 220565 261915 264320
+ 264615 263029 259766 254916 248417 240093 229693 217053 260169 260246 259189 256896
+ 253361 248612 242540 234998 225848 215386 259060 257307 255000 252005 250427 248358
+ 245990 248042 246062 243529 241370 243782 240696 237590 235794 238714 231512 223317
+ 214405 258202 254922 251498 249784 247922 245956 247848 246345 244520 242989 244624
+ 242172 239751 238419 240754 237196 234515 233376 235959 231745 227956 226991 230633
+ 221504 213702 257407 252680 248163 245971 243997 241952 243999 242435 240884 238679
+ 240226 238979 236854 235008 236948 234561 232125 230648 232895 229483 225793 224045
+ 227752 219522 212978 256529 250229 244524 241931 239773 237380 239602 237931 236342
+ 233839 235522 234696 232897 230529 232220 230925 228752 225778 228587 225439 222119
+ 219780 222707 216830 212043 255388 247183 240129 234374 232722 231049 227368 229264
+ 229412 227255 223116 225405 224890 222290 218594 220880 217040 213369 210787 253607
+ 242987 234435 227673 222377 217907 213872 210643 208822 208901 250036 236722 226952
+ 219507 213705 209039 205314 202841 202376 205273
+ </DataArray>
+ <DataArray type="Float32" Name="wetting saturation" NumberOfComponents="1" format="ascii">
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 0.958403 0.94218 0.924758 0.942679
+ 0.964063 0.977305 0.962871 0.94743 0.985015 0.989905 0.980925 0.97332 0.99773 1 1 0.989424
+ 1 1 1 1 1 1 1 1 1 0.966258 0.949408 0.96573
+ 0.980345 0.955664 0.96463 0.970305 0.965664 0.972637 0.979666 0.983491 0.975984 0.989357 1 1
+ 1 1 1 1 1 1 1 1 1 1 0.988525 0.976799
+ 0.985079 0.993889 0.974166 0.976355 0.983841 0.98158 0.980871 0.990961 1 0.990513 1 1
+ 1 1 1 1 1 1 1 1 1 0.993177 1 1
+ 0.99124 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1
+ </DataArray>
+ <DataArray type="Float32" Name="C^w from cellData" NumberOfComponents="1" format="ascii">
+ 199.956 199.956 199.957 199.956 199.956 199.956 199.955 199.954 199.953 199.952 199.956 199.956
+ 199.956 199.956 199.955 199.955 199.954 199.954 199.953 199.951 199.956 199.955 199.955 199.955
+ 199.955 199.954 199.954 199.953 199.952 199.951 199.955 199.955 199.955 199.954 199.954 199.954
+ 199.78 199.839 199.952 199.952 199.944 199.963 199.953 199.952 199.948 199.97 199.952 199.952
+ 199.951 199.955 199.955 199.955 199.954 199.759 199.76 199.954 191.819 188.626 185.175 188.711
+ 192.782 195.315 192.509 189.564 196.777 197.414 195.902 194.491 198.998 199.921 199.367 197.171
+ 199.952 199.951 199.955 199.955 199.954 199.954 199.899 199.897 199.954 193.15 189.941 193.064
+ 195.867 191.115 192.81 193.887 193.015 194.327 195.587 195.899 194.941 197.063 199.621 199.951
+ 198.989 199.952 199.951 199.955 199.955 199.954 199.953 199.909 199.932 199.954 197.425 195.17
+ 196.67 198.224 194.626 195.018 195.938 195.944 195.676 197.649 199.852 197.531 199.918 199.948
+ 199.948 199.948 199.951 199.951 199.955 199.954 199.954 199.953 199.515 198.076 199.749 199.953
+ 197.538 199.019 199.942 199.75 199.908 199.944 199.944 199.944 199.951 199.951 199.951 199.955
+ 199.954 199.953 199.952 199.952 199.951 199.951 199.951 199.951 199.951 199.955 199.953 199.952
+ 199.952 199.951 199.951 199.95 199.95 199.95 199.95
+ </DataArray>
+ <DataArray type="Float32" Name="C^n from cellData" NumberOfComponents="1" format="ascii">
+ 0 0 0 7.92093e-11 3.24337e-10 2.11506e-10 7.22924e-11 5.06063e-12 1.64029e-13 4.65438e-15 0 0
+ 0 1.88533e-08 8.24478e-08 4.3008e-08 1.85189e-08 1.99666e-09 4.60804e-11 8.19964e-13 0 0 0 1.73815e-06
+ 7.97736e-06 3.37589e-06 1.7112e-06 4.81668e-07 7.25419e-09 8.47836e-11 0 0 0 7.64542e-05 0.000283102 0.000449254
+ 0.0049682 0.0034181 0.00015551 0.000118114 0.00161755 0.00199534 8.04903e-05 6.32126e-05 0.000872258 0.00114886 8.24128e-05 7.38825e-07
+ 5.91603e-09 0 0 0 3.86124e-08 0.00288273 0.00422177 4.29922e-07 0.0390057 0.0482056 0.0580305 0.0479313
+ 0.0350357 0.0270904 0.0351612 0.0444121 0.0224289 0.0194069 0.024334 0.028898 0.014863 0.000738621 0.00734964 0.0193333
+ 2.25224e-05 1.47795e-07 0 0 0 3.87075e-06 4.1301e-06 3.80509e-06 3.53094e-06 0.0337913 0.043294 0.0335691
+ 0.0253952 0.039364 0.0339264 0.0304702 0.0333549 0.0291033 0.0248916 0.0225784 0.0270297 0.0192639 0.00401202 0.000138304
+ 0.0133138 4.04253e-06 1.1637e-08 0 0 0 0 0 0 0 0.0205024 0.0269822
+ 0.0220611 0.0172593 0.0282639 0.0268239 0.0223712 0.0238393 0.0240806 0.0183549 0.00126675 0.0185827 0.000437198 8.65817e-05
+ 8.65262e-05 8.70207e-05 1.00488e-07 9.02684e-11 0 0 0 0 0.00508785 0.0173203 0.00222419 0
+ 0.0182402 0.0131763 0.000121281 0.00223952 0.000507655 0.000102718 0.000102686 0.000102702 5.57498e-08 2.49248e-11 1.48737e-14 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0
+ </DataArray>
+ <DataArray type="Float32" Name="capillary pressure" NumberOfComponents="1" format="ascii">
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 415.975 578.198 752.415 573.214
+ 359.373 226.953 371.287 525.696 149.847 100.949 190.748 266.803 22.6961 0 0 105.758
+ 0 0 0 0 0 0 0 0 0 337.422 505.916 342.704
+ 196.554 443.359 353.7 296.95 343.355 273.626 203.344 165.092 240.159 106.432 0 0
+ 0 0 0 0 0 0 0 0 0 0 114.752 232.013
+ 149.215 61.1145 258.339 236.444 161.588 184.195 191.294 90.3875 0 94.8702 0 0
+ 0 0 0 0 0 0 0 0 0 68.229 0 0
+ 87.5962 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0
+ </DataArray>
+ <DataArray type="Float32" Name="wetting density" NumberOfComponents="1" format="ascii">
+ 999.78 999.783 999.783 999.783 999.781 999.779 999.776 999.771 999.766 999.759 999.779 999.78
+ 999.78 999.779 999.777 999.775 999.772 999.768 999.763 999.757 999.778 999.778 999.777 999.776
+ 999.774 999.772 999.769 999.766 999.761 999.756 999.777 999.776 999.775 999.774 999.774 999.773
+ 999.78 999.778 999.771 999.77 999.772 999.774 999.768 999.767 999.768 999.77 999.764 999.76
+ 999.756 999.777 999.775 999.773 999.773 999.777 999.779 999.772 999.798 999.797 999.796 999.797
+ 999.796 999.794 999.794 999.795 999.793 999.791 999.791 999.792 999.79 999.764 999.776 999.789
+ 999.759 999.755 999.776 999.774 999.772 999.771 999.77 999.769 999.77 999.796 999.795 999.794
+ 999.795 999.794 999.793 999.792 999.793 999.791 999.79 999.789 999.79 999.788 999.769 999.761
+ 999.787 999.758 999.755 999.776 999.773 999.77 999.769 999.768 999.767 999.768 999.793 999.792
+ 999.791 999.792 999.791 999.79 999.789 999.79 999.789 999.788 999.764 999.788 999.762 999.76
+ 999.759 999.76 999.757 999.755 999.775 999.771 999.768 999.765 999.774 999.789 999.766 999.763
+ 999.788 999.787 999.76 999.765 999.762 999.76 999.758 999.759 999.757 999.755 999.754 999.775
+ 999.769 999.765 999.762 999.76 999.757 999.756 999.754 999.753 999.753 999.773 999.766 999.762
+ 999.758 999.755 999.753 999.751 999.75 999.75 999.751
+ </DataArray>
+ <DataArray type="Float32" Name="nonwetting density" NumberOfComponents="1" format="ascii">
+ 3.2598 3.32182 3.33927 3.32595 3.28807 3.22788 3.14487 3.03596 2.89464 2.70751 3.21616 3.24575
+ 3.24938 3.22987 3.18974 3.13006 3.05012 2.94773 2.81979 2.66431 3.19468 3.19564 3.18264 3.15443
+ 3.11094 3.05253 2.97783 2.88506 2.7725 2.6438 3.18105 3.15948 3.1311 3.09427 3.07485 3.0494
+ 3.02026 3.04552 3.02116 2.99 2.96344 2.99311 2.95514 2.91694 2.89485 2.93077 2.84218 2.74137
+ 2.63174 3.1705 3.13014 3.08803 3.06694 3.04403 3.01985 3.04312 3.02462 3.00217 2.98333 3.00345
+ 2.97329 2.94352 2.92713 2.95584 2.91209 2.87911 2.8651 2.89687 2.84504 2.79842 2.78656 2.83135
+ 2.71906 2.62309 3.16071 3.10257 3.047 3.02003 2.99575 2.9706 2.99578 2.97653 2.95744 2.93033
+ 2.94936 2.93402 2.90787 2.88517 2.90904 2.87967 2.84971 2.83154 2.85918 2.81721 2.77182 2.75032
+ 2.79592 2.69468 2.61418 3.14992 3.07242 3.00224 2.97034 2.94379 2.91436 2.94169 2.92113 2.90158
+ 2.8708 2.89151 2.88133 2.85921 2.83008 2.85088 2.83494 2.80822 2.77164 2.80619 2.76746 2.72662
+ 2.69785 2.73386 2.66156 2.60269 3.13588 3.03495 2.94817 2.87738 2.85705 2.83647 2.7912 2.81452
+ 2.81634 2.7898 2.73889 2.76705 2.76072 2.72873 2.68326 2.71139 2.66415 2.619 2.58723 3.11396
+ 2.98333 2.87813 2.79495 2.72981 2.67481 2.62518 2.58546 2.56306 2.56403 3.07004 2.90626 2.78608
+ 2.6945 2.62312 2.56572 2.51991 2.48949 2.48376 2.5194
+ </DataArray>
+ <DataArray type="Float32" Name="mobility w_phase" NumberOfComponents="1" format="ascii">
+ 765.842 765.844 765.845 765.845 765.843 765.84 765.837 765.832 765.826 765.818 765.84 765.841
+ 765.841 765.841 765.839 765.836 765.833 765.828 765.823 765.816 765.839 765.839 765.839 765.837
+ 765.835 765.833 765.83 765.826 765.821 765.815 765.838 765.838 765.836 765.835 765.834 765.833
+ 765.831 765.833 765.832 765.83 765.829 765.83 765.829 765.827 765.826 765.828 765.824 765.82
+ 765.815 765.838 765.836 765.834 765.833 765.833 765.831 765.832 733.975 721.55 708.207 721.932
+ 738.307 748.447 737.393 725.569 754.351 758.094 751.217 745.394 764.086 765.822 765.821 757.724
+ 765.819 765.814 765.838 765.835 765.833 765.831 765.83 765.829 765.83 739.989 727.084 739.582
+ 750.776 731.874 738.739 743.084 739.531 744.87 750.251 753.18 747.432 757.672 765.821 765.82
+ 765.822 765.818 765.814 765.837 765.834 765.831 765.829 765.828 765.827 765.828 757.039 748.058
+ 754.398 761.146 746.041 747.717 753.448 751.718 751.174 758.9 765.821 758.557 765.821 765.819
+ 765.818 765.819 765.816 765.814 765.836 765.832 765.828 765.825 765.824 760.598 765.822 765.823
+ 759.114 765.822 765.819 765.821 765.82 765.819 765.817 765.818 765.816 765.814 765.813 765.836
+ 765.83 765.825 765.822 765.819 765.817 765.815 765.813 765.812 765.812 765.834 765.827 765.821
+ 765.818 765.814 765.812 765.81 765.809 765.808 765.81
+ </DataArray>
+ <DataArray type="Float32" Name="mobility nw_phase" NumberOfComponents="1" format="ascii">
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 2380.54 3308.98 4306.1 3280.46
+ 2056.71 1298.9 2124.99 3008.64 857.628 577.781 1091.77 1527.03 129.926 0 0 605.337
+ 0 0 0 0 0 0 0 0 0 1931.08 2895.43 1961.39
+ 1124.91 2537.46 2024.37 1699.6 1965.16 1566.11 1163.88 944.951 1374.59 609.208 0 0
+ 0 0 0 0 0 0 0 0 0 0 656.757 1327.91
+ 854.044 349.784 1478.61 1353.32 924.896 1054.28 1094.93 517.39 0 543.048 0 0
+ 0 0 0 0 0 0 0 0 0 390.531 0 0
+ 501.399 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0
+ </DataArray>
+ <DataArray type="Float32" Name="mass fraction H2O in liquid-phase" NumberOfComponents="1" format="ascii">
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 0.999999 0.999998
+ 0.999975 0.999983 0.999999 0.999999 0.999992 0.99999 1 1 0.999996 0.999994 1 1
+ 1 1 1 1 1 0.999986 0.999979 1 0.999928 0.999928 0.999929 0.999928
+ 0.999929 0.999929 0.99993 0.999929 0.99993 0.999931 0.999931 0.999931 0.999932 0.999996 0.999963 0.999932
+ 1 1 1 1 1 1 1 1 1 0.999929 0.999929 0.99993
+ 0.999929 0.99993 0.99993 0.999931 0.99993 0.999931 0.999932 0.999932 0.999932 0.999932 0.99998 0.999999
+ 0.999933 1 1 1 1 1 1 1 1 1 0.99993 0.99993
+ 0.999931 0.999931 0.999931 0.999932 0.999932 0.999932 0.999932 0.999933 0.999994 0.999933 0.999998 1
+ 1 1 1 1 1 1 1 1 0.999974 0.999932 0.999989 1
+ 0.999932 0.999934 0.999999 0.999989 0.999997 0.999999 0.999999 0.999999 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1
+ </DataArray>
+ <DataArray type="Float32" Name="mass fraction H2O in gas-phase" NumberOfComponents="1" format="ascii">
+ 0.00288297 0.00282914 0.00281436 0.00282563 0.00285818 0.00291148 0.00298834 0.00309554 0.00324668 0.00347108 0.00292209 0.00289545
+ 0.00289221 0.00290968 0.0029463 0.00300247 0.00308117 0.0031882 0.00333285 0.00352736 0.00294173 0.00294086 0.00295287 0.00297928
+ 0.00302092 0.00307873 0.00315597 0.00325746 0.00338971 0.00355472 0.00295434 0.00297451 0.00300147 0.0030372 0.00305639 0.0030819
+ 0.00311163 0.00308582 0.0031107 0.00314312 0.0031713 0.00313986 0.0031802 0.00322185 0.00324644 0.00320665 0.00330661 0.00342821
+ 0.00357102 0.00296418 0.00300239 0.00304334 0.00306427 0.00308733 0.00311205 0.00308825 0.0031124 0.0031378 0.00315993 0.00313639
+ 0.00316549 0.00319578 0.00321565 0.0031864 0.00322926 0.00326559 0.00328283 0.00324785 0.0033036 0.00335831 0.0033726 0.00332077
+ 0.00345633 0.0035828 0.00297335 0.00302908 0.00308432 0.00311186 0.00313709 0.00316365 0.00313705 0.00316175 0.00318441 0.00321174
+ 0.00318905 0.00320905 0.00323673 0.00326145 0.00323529 0.00326736 0.00330075 0.0033214 0.00329033 0.00333746 0.00339054 0.00341705
+ 0.00336131 0.00348761 0.00359501 0.00298354 0.0030588 0.00313031 0.00316392 0.00319246 0.00322471 0.00319475 0.00321879 0.00324209
+ 0.00327573 0.00325104 0.00326526 0.00329025 0.00332307 0.00329913 0.00331779 0.00334791 0.00339077 0.0033504 0.00339588 0.00344675
+ 0.0034835 0.00343762 0.003531 0.00361088 0.0029969 0.00309657 0.00318772 0.00326615 0.00328939 0.00331423 0.003367 0.0033391
+ 0.00333822 0.00336868 0.0034313 0.00339639 0.00340418 0.00344408 0.00350245 0.00346611 0.00352758 0.0035884 0.00363245 0.00301799
+ 0.00315016 0.0032653 0.00336248 0.00344273 0.00351351 0.00357995 0.00363494 0.00366671 0.00366533 0.00306117 0.00323369 0.00337319
+ 0.00348784 0.00358275 0.0036629 0.00372951 0.00377508 0.00378378 0.00373026
+ </DataArray>
+ <DataArray type="Float32" Name="volume Error" NumberOfComponents="1" format="ascii">
+ -2.06043e-08 -6.30168e-08 -1.07765e-07 -1.34875e-07 -1.44058e-07 -1.32446e-07 -9.93583e-08 -4.56964e-08 3.30681e-08 3.7297e-08 -3.16524e-08 -8.05503e-08
+ -1.23608e-07 -1.54229e-07 -1.66849e-07 -1.56086e-07 -1.21156e-07 -6.86887e-08 -3.7375e-09 1.23142e-08 -3.2728e-08 -9.04513e-08 -1.41352e-07 -1.84024e-07
+ -2.03549e-07 -1.94603e-07 -1.51116e-07 -9.3437e-08 -3.07356e-08 -4.4461e-10 -3.34329e-08 -9.76669e-08 -1.61881e-07 -2.28732e-07 -4.87404e-07 -4.80496e-07
+ -0.000171057 -0.000113132 -2.3303e-06 -2.27194e-06 -8.65653e-06 1.03443e-05 -1.02469e-06 -1.06813e-06 -5.09582e-06 1.75149e-05 -9.97759e-07 -5.0181e-08
+ -9.69074e-09 -3.28179e-08 -1.00344e-07 -1.76566e-07 -4.7453e-07 -0.000193902 -0.000191828 -2.99158e-08 0.000199412 0.00025602 0.000296497 0.000240859
+ 2.36815e-05 -9.42712e-05 -1.29871e-05 0.000137294 -0.000173932 -0.000517111 -0.000233351 -0.000122325 -0.000493446 -3.07238e-05 -0.000580904 -0.00066566
+ -6.48427e-08 -1.62292e-08 -2.92768e-08 -9.33132e-08 -1.70544e-07 -6.36117e-07 -5.53033e-05 -5.68674e-05 -1.92721e-07 -5.04253e-05 0.000117187 -2.9785e-05
+ -0.000150674 3.60848e-05 -6.49446e-05 -0.000123579 -6.68367e-05 -0.000150944 -0.000297419 -0.000756702 -0.000206197 -0.000761879 -0.000329041 -7.86535e-07
+ -0.000955115 -6.93698e-08 -2.04984e-08 -2.21075e-08 -7.50255e-08 -1.39324e-07 -4.34793e-07 -4.5046e-05 -2.09748e-05 3.83093e-08 -0.000227929 -0.000138876
+ -0.00029499 -0.00050176 -0.000156979 -0.000204074 -0.000788184 -0.000323943 -0.000451927 -0.000492895 -9.9739e-05 -0.000520946 -3.442e-05 -4.19408e-06
+ -3.72476e-06 -3.72482e-06 -6.86633e-08 -2.33817e-08 -1.21327e-08 -5.06069e-08 -1.03849e-07 -1.45336e-07 -0.000435082 -0.000507778 -0.000202047 7.94434e-08
+ -0.000661198 -0.000925088 -9.49768e-06 -0.0002008 -4.37327e-05 -7.85378e-06 -7.38488e-06 -7.38494e-06 -1.02385e-07 -6.84013e-08 -2.59389e-08 5.05843e-10
+ -2.14527e-08 -7.36224e-08 -1.10385e-07 -1.27868e-07 -1.2785e-07 -1.1446e-07 -9.35935e-08 -6.52531e-08 -2.81496e-08 2.42134e-08 1.44817e-08 -4.95027e-08
+ -8.86586e-08 -1.08728e-07 -1.12056e-07 -1.02352e-07 -8.31326e-08 -5.30539e-08 -2.00632e-08
+ </DataArray>
+ <DataArray type="Float32" Name="Error Correction" NumberOfComponents="1" format="ascii">
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0
+ </DataArray>
+ <DataArray type="Float32" Name="dv_dp" NumberOfComponents="1" format="ascii">
+ -9.55458e-11 -9.55402e-11 -9.55458e-11 -9.55486e-11 -9.55347e-11 -9.55569e-11 -9.55458e-11 -9.5543e-11 -9.55347e-11 -9.5543e-11 -9.55319e-11 -9.5543e-11
+ -9.55402e-11 -9.55375e-11 -9.55375e-11 -9.55402e-11 -9.5543e-11 -9.55458e-11 -9.5543e-11 -9.55375e-11 -9.5543e-11 -9.55486e-11 -9.55375e-11 -9.55541e-11
+ -9.55375e-11 -9.55375e-11 -9.5543e-11 -9.55458e-11 -9.55541e-11 -9.55513e-11 -9.55319e-11 -9.55402e-11 -9.55486e-11 -9.55486e-11 -9.55513e-11 -9.5543e-11
+ -9.54681e-11 -9.54958e-11 -9.55486e-11 -9.55319e-11 -9.5543e-11 -9.55541e-11 -9.55513e-11 -9.55458e-11 -9.55458e-11 -9.55569e-11 -9.5543e-11 -9.55458e-11
+ -9.55541e-11 -9.55486e-11 -9.55347e-11 -9.55541e-11 -9.5543e-11 -9.54514e-11 -9.54486e-11 -9.55264e-11 -5.29354e-08 -6.63459e-08 -8.07696e-08 -6.59242e-08
+ -4.91753e-08 -3.88138e-08 -5.08388e-08 -6.29582e-08 -3.28557e-08 -2.90967e-08 -3.67622e-08 -4.26823e-08 -2.25612e-08 -9.55375e-11 -9.5271e-11 -2.99535e-08
+ -9.55402e-11 -9.55486e-11 -9.55375e-11 -9.55402e-11 -9.55486e-11 -9.55347e-11 -9.55097e-11 -9.55208e-11 -9.55458e-11 -4.73408e-08 -6.13218e-08 -4.8452e-08
+ -3.62578e-08 -5.66143e-08 -4.96806e-08 -4.53295e-08 -4.88084e-08 -4.34705e-08 -3.79843e-08 -3.49718e-08 -4.09594e-08 -3.00655e-08 -9.54153e-11 -9.5543e-11
+ -9.51406e-11 -9.55513e-11 -9.55597e-11 -9.55402e-11 -9.55375e-11 -9.55402e-11 -9.55402e-11 -9.55291e-11 -9.55486e-11 -9.55513e-11 -2.98626e-08 -3.97382e-08
+ -3.31996e-08 -2.56703e-08 -4.21765e-08 -4.06503e-08 -3.47149e-08 -3.63489e-08 -3.71098e-08 -2.86384e-08 -9.55291e-11 -2.90684e-08 -9.55569e-11 -9.55569e-11
+ -9.55569e-11 -9.55569e-11 -9.55513e-11 -9.55347e-11 -9.5543e-11 -9.5543e-11 -9.55458e-11 -9.55458e-11 -9.53432e-11 -2.6664e-08 -9.54736e-11 -9.5543e-11
+ -2.84276e-08 -9.51461e-11 -9.55569e-11 -9.54653e-11 -9.55569e-11 -9.55569e-11 -9.55569e-11 -9.55569e-11 -9.55597e-11 -9.55458e-11 -9.55652e-11 -9.55402e-11
+ -9.55458e-11 -9.55402e-11 -9.55486e-11 -9.5543e-11 -9.55513e-11 -9.55486e-11 -9.55513e-11 -9.55458e-11 -9.55486e-11 -9.55347e-11 -9.55402e-11 -9.55486e-11
+ -9.55597e-11 -9.55735e-11 -9.55486e-11 -9.55597e-11 -9.55569e-11 -9.55569e-11 -9.55541e-11
+ </DataArray>
+ <DataArray type="Float32" Name="dV_dC1" NumberOfComponents="1" format="ascii">
+ 0.00100022 0.00100022 0.00100022 0.00100022 0.00100022 0.00100022 0.00100022 0.00100023 0.00100023 0.00100024 0.00100022 0.00100022
+ 0.00100022 0.00100022 0.00100022 0.00100022 0.00100023 0.00100023 0.00100024 0.00100024 0.00100022 0.00100022 0.00100022 0.00100022
+ 0.00100023 0.00100023 0.00100023 0.00100023 0.00100024 0.00100024 0.00100022 0.00100022 0.00100022 0.00100023 0.00100023 0.00100023
+ 0.00100023 0.00100023 0.00100023 0.00100023 0.00100023 0.00100023 0.00100023 0.00100023 0.00100023 0.00100023 0.00100024 0.00100024
+ 0.00100024 0.00100022 0.00100022 0.00100023 0.00100023 0.00100023 0.00100023 0.00100023 0.000976249 0.000976266 0.000976284 0.000976266
+ 0.000976246 0.000976233 0.000976248 0.000976263 0.000976226 0.000976222 0.000976232 0.000976239 0.000976215 0.00100024 0.00100024 0.000976224
+ 0.00100024 0.00100024 0.00100022 0.00100023 0.00100023 0.00100023 0.00100023 0.00100023 0.00100023 0.000976243 0.000976261 0.000976245
+ 0.00097623 0.000976255 0.000976247 0.000976242 0.000976246 0.00097624 0.000976234 0.00097623 0.000976237 0.000976225 0.00100024 0.00100024
+ 0.00100024 0.00100024 0.00100024 0.00100022 0.00100023 0.00100023 0.00100023 0.00100023 0.00100023 0.00100023 0.000976223 0.000976235
+ 0.000976228 0.000976218 0.000976238 0.000976237 0.00097623 0.000976232 0.000976233 0.000976223 0.00100024 0.000976223 0.00100024 0.00100024
+ 0.00100024 0.00100024 0.00100024 0.00100025 0.00100022 0.00100023 0.00100023 0.00100023 0.00100024 0.00097622 0.00100024 0.00100024
+ 0.000976223 0.00100024 0.00100024 0.00100024 0.00100024 0.00100024 0.00100024 0.00100024 0.00100024 0.00100024 0.00100025 0.00100023
+ 0.00100023 0.00100023 0.00100024 0.00100024 0.00100024 0.00100024 0.00100025 0.00100025 0.00100025 0.00100023 0.00100023 0.00100024
+ 0.00100024 0.00100024 0.00100025 0.00100025 0.00100025 0.00100025 0.00100025
+ </DataArray>
+ <DataArray type="Float32" Name="dV_dC2" NumberOfComponents="1" format="ascii">
+ 0.000622207 0.000622206 0.000622205 0.000622207 0.000622207 0.00062221 0.00062221 0.000622211 0.000622217 0.000622221 0.000622209 0.000622207
+ 0.000622207 0.000622209 0.000622209 0.000622212 0.000622215 0.000622215 0.000622217 0.000622222 0.00062221 0.000622209 0.00062221 0.00062221
+ 0.000622211 0.000622211 0.000622216 0.000622216 0.000622218 0.000622222 0.00062221 0.000622209 0.000622213 0.000622212 0.000622212 0.000622213
+ 0.000622214 0.000622213 0.000622214 0.000622214 0.000622215 0.000622214 0.000622215 0.000622216 0.000622217 0.000622216 0.000622218 0.00062222
+ 0.000622225 0.000622211 0.000622211 0.000622212 0.000622212 0.000622213 0.000622214 0.000622214 0.331724 0.334196 0.336298 0.334054
+ 0.337469 0.34094 0.342852 0.339458 0.344721 0.34885 0.350538 0.346533 0.353425 0.000622219 0.000622219 0.354955
+ 0.000622221 0.000622224 0.000622211 0.00062221 0.000622214 0.000622214 0.000622216 0.000622215 0.000622216 0.337096 0.339269 0.342453
+ 0.340241 0.342015 0.345161 0.347948 0.345007 0.348653 0.352494 0.354943 0.351241 0.356761 0.00062222 0.00062222
+ 0.0473524 0.000622222 0.000622222 0.000622211 0.000622214 0.000622215 0.000622214 0.000622215 0.000622217 0.000622216 0.34359 0.345906
+ 0.34969 0.347198 0.348405 0.351216 0.355072 0.352244 0.354323 0.357683 0.00062222 0.357935 0.000622221 0.000622221
+ 0.000622221 0.000622221 0.000622222 0.000622225 0.000622213 0.000622214 0.000622216 0.000622218 0.000622218 0.353954 0.000622219 0.000622216
+ 0.356606 0.000622219 0.000622221 0.00062222 0.000622221 0.000622221 0.000622221 0.000622221 0.000622222 0.000622222 0.000622225 0.000622212
+ 0.000622215 0.000622218 0.00062222 0.000622221 0.000622222 0.000622221 0.000622225 0.000622226 0.000622226 0.000622212 0.000622216 0.000622219
+ 0.000622219 0.000622225 0.000622224 0.000622227 0.000622228 0.000622226 0.000622227
+ </DataArray>
+ <DataArray type="Float32" Name="updEstimate comp 1" NumberOfComponents="1" format="ascii">
+ 1.96191e-07 5.65419e-07 8.67442e-07 1.16772e-06 1.41852e-06 1.59637e-06 1.67433e-06 1.62181e-06 1.46218e-06 6.73928e-07 1.49355e-07 5.20991e-07
+ 8.77534e-07 1.20391e-06 1.47312e-06 1.67827e-06 1.79388e-06 1.74211e-06 1.4253e-06 5.89113e-07 1.8251e-07 5.84287e-07 9.84557e-07 1.33043e-06
+ 1.58795e-06 1.88716e-06 2.14072e-06 2.1893e-06 1.68525e-06 6.60889e-07 2.32804e-07 7.08527e-07 1.16909e-06 1.26094e-06 4.26079e-07 -1.46302e-07
+ 0.000162797 0.000402969 1.04153e-06 1.40354e-06 3.42569e-05 5.36551e-05 1.87408e-06 2.29902e-06 -0.000354575 4.31895e-05 -0.000337097 2.27123e-06
+ 8.51136e-07 2.97167e-07 8.79149e-07 1.41239e-06 1.69466e-06 -0.000209079 -0.000254725 1.89647e-06 0.0159753 0.00741613 -0.0149873 0.0098538
+ 0.00445898 0.00284899 -0.0175628 -0.0183546 0.00249114 -0.00161958 -0.0171747 -0.0164226 -0.0385423 -0.00173165 -0.032402 -0.0491562
+ 1.94549e-06 1.14106e-06 3.76499e-07 1.09647e-06 1.74428e-06 2.07606e-06 -0.000485905 -0.00050629 2.38583e-06 0.00852323 -0.0141509 -0.0113201
+ 0.00721519 -0.0242673 -0.0251288 -0.0254044 -0.0225277 -0.0246707 -0.0247783 -0.0461592 -0.0265096 -0.0990562 -0.0331152 0.168581
+ -0.0656533 3.93268e-06 1.50981e-06 4.71055e-07 1.36107e-06 2.17406e-06 2.65833e-06 -0.000521456 -0.000701362 2.97164e-06 0.00774117 -0.00934717
+ -0.0105364 -0.00796375 -0.0200449 -0.0247371 -0.0434241 -0.0179707 -0.0167422 -0.195948 -0.0425728 -0.209916 0.0903865 -0.566319
+ -0.0215949 0.524278 0.0262807 1.7748e-06 5.75842e-07 1.65176e-06 2.6669e-06 3.79176e-06 -0.0157004 -0.056696 -0.0210985 5.1503e-06
+ -0.111661 -0.0823378 -0.0919392 -0.0239204 0.158124 -0.400725 -0.0856641 0.321697 9.91621e-05 5.17468e-06 1.66925e-06 6.84764e-07
+ 1.93948e-06 3.08007e-06 4.47094e-06 6.32372e-06 8.25817e-06 -0.0111811 7.12758e-06 4.40992e-06 1.44529e-06 8.38604e-07 2.19935e-06 3.29397e-06
+ 4.62327e-06 6.11108e-06 7.30653e-06 7.46325e-06 6.18586e-06 3.98411e-06 1.34396e-06
+ </DataArray>
+ <DataArray type="Float32" Name="updEstimate comp 2" NumberOfComponents="1" format="ascii">
+ 0 0 0 2.56695e-12 1.01682e-11 5.4929e-12 2.7076e-12 3.77806e-13 1.45777e-14 4.65785e-16 0 0
+ 0 4.35318e-10 1.84521e-09 8.20345e-10 4.67878e-10 1.55233e-10 4.05339e-12 7.98731e-14 0 0 0 2.53675e-08
+ 1.13977e-07 4.62554e-08 2.67523e-08 3.47027e-08 6.24337e-10 8.08183e-12 0 0 0 5.24956e-07 2.11376e-06 3.13585e-06
+ 6.27237e-06 6.55836e-06 1.54356e-06 1.1776e-06 1.08285e-05 1.169e-05 8.95384e-07 7.03186e-07 1.27638e-05 1.1862e-05 3.86232e-06 5.52475e-08
+ 5.02818e-10 0 0 0 3.05196e-10 -4.04385e-07 -5.74498e-06 2.70426e-09 -2.00359e-05 -1.84233e-05 9.57014e-06 -2.09697e-05
+ -6.73763e-06 -1.68319e-06 4.90065e-05 4.46177e-05 -1.29666e-06 8.81908e-06 5.03791e-05 4.54805e-05 0.000257584 0.000104715 0.000461455 7.4356e-05
+ 2.03391e-06 1.07982e-08 0 0 0 -2.53452e-09 -4.11013e-10 -3.05844e-09 -4.93963e-09 -1.51408e-05 3.53213e-05 3.26247e-05
+ -1.26256e-05 6.64858e-05 7.37487e-05 7.74501e-05 6.62767e-05 7.20134e-05 7.72354e-05 9.9888e-05 7.83046e-05 8.2677e-05 0.000471595 1.24921e-05
+ 0.000869579 9.23757e-07 1.36303e-09 0 0 0 0 0 0 0 -1.16614e-05 2.57312e-05
+ 3.52921e-05 -2.41436e-06 6.02181e-05 7.40684e-05 9.28649e-05 6.87393e-05 0.000137388 4.94424e-05 0.000697467 1.23993e-06 0.000428892 -2.45228e-07
+ -7.48738e-08 1.02576e-06 4.42541e-08 2.16156e-11 0 0 0 0 0.000203801 4.69961e-05 0.000301059 0
+ 7.06372e-05 0.00090251 1.06256e-06 0.000418631 0.000517888 -2.05866e-07 -7.8177e-08 1.33097e-07 4.65682e-08 5.60091e-12 3.26627e-15 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0
+ </DataArray>
+ </CellData>
+ <Points>
+ <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+ 0 0 0 1 0 0 0 1 0 1 1 0
+ 2 0 0 2 1 0 3 0 0 3 1 0
+ 4 0 0 4 1 0 5 0 0 5 1 0
+ 6 0 0 6 1 0 7 0 0 7 1 0
+ 8 0 0 8 1 0 9 0 0 9 1 0
+ 10 0 0 10 1 0 0 2 0 1 2 0
+ 2 2 0 3 2 0 4 2 0 5 2 0
+ 6 2 0 7 2 0 8 2 0 9 2 0
+ 10 2 0 0 3 0 1 3 0 2 3 0
+ 3 3 0 4 3 0 5 3 0 6 3 0
+ 7 3 0 8 3 0 9 3 0 10 3 0
+ 0 4 0 1 4 0 2 4 0 3 4 0
+ 4 4 0 4.5 3 0 4 3.5 0 4.5 3.5 0
+ 5 3.5 0 4.5 4 0 5 4 0 5.5 3 0
+ 5.5 3.5 0 6 3.5 0 5.5 4 0 6 4 0
+ 6.5 3 0 6.5 3.5 0 7 3.5 0 6.5 4 0
+ 7 4 0 8 4 0 9 4 0 10 4 0
+ 0 5 0 1 5 0 2 5 0 3 5 0
+ 3.5 4 0 3 4.5 0 3.5 4.5 0 4 4.5 0
+ 3.5 5 0 4 5 0 4.5 4.5 0 5 4.5 0
+ 4.5 5 0 5 5 0 5.5 4.5 0 6 4.5 0
+ 5.5 5 0 6 5 0 6.5 4.5 0 7 4.5 0
+ 6.5 5 0 7 5 0 7.5 4 0 7.5 4.5 0
+ 8 4.5 0 7.5 5 0 8 5 0 9 5 0
+ 10 5 0 0 6 0 1 6 0 2 6 0
+ 3 6 0 3 5.5 0 3.5 5.5 0 4 5.5 0
+ 3.5 6 0 4 6 0 4.5 5.5 0 5 5.5 0
+ 4.5 6 0 5 6 0 5.5 5.5 0 6 5.5 0
+ 5.5 6 0 6 6 0 6.5 5.5 0 7 5.5 0
+ 6.5 6 0 7 6 0 7.5 5.5 0 8 5.5 0
+ 7.5 6 0 8 6 0 9 6 0 10 6 0
+ 0 7 0 1 7 0 2 7 0 3 7 0
+ 3 6.5 0 3.5 6.5 0 4 6.5 0 3.5 7 0
+ 4 7 0 4.5 6.5 0 5 6.5 0 4.5 7 0
+ 5 7 0 5.5 6.5 0 6 6.5 0 5.5 7 0
+ 6 7 0 6.5 6.5 0 7 6.5 0 6.5 7 0
+ 7 7 0 7.5 6.5 0 8 6.5 0 7.5 7 0
+ 8 7 0 9 7 0 10 7 0 0 8 0
+ 1 8 0 2 8 0 3 8 0 4 8 0
+ 4 7.5 0 4.5 7.5 0 5 7.5 0 4.5 8 0
+ 5 8 0 5.5 7.5 0 6 7.5 0 5.5 8 0
+ 6 8 0 6.5 7.5 0 7 7.5 0 6.5 8 0
+ 7 8 0 8 8 0 9 8 0 10 8 0
+ 0 9 0 1 9 0 2 9 0 3 9 0
+ 4 9 0 5 9 0 6 9 0 7 9 0
+ 8 9 0 9 9 0 10 9 0 0 10 0
+ 1 10 0 2 10 0 3 10 0 4 10 0
+ 5 10 0 6 10 0 7 10 0 8 10 0
+ 9 10 0 10 10 0
+ </DataArray>
+ </Points>
+ <Cells>
+ <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+ 0 1 3 2 1 4 5 3 4 6 7 5
+ 6 8 9 7 8 10 11 9 10 12 13 11
+ 12 14 15 13 14 16 17 15 16 18 19 17
+ 18 20 21 19 2 3 23 22 3 5 24 23
+ 5 7 25 24 7 9 26 25 9 11 27 26
+ 11 13 28 27 13 15 29 28 15 17 30 29
+ 17 19 31 30 19 21 32 31 22 23 34 33
+ 23 24 35 34 24 25 36 35 25 26 37 36
+ 26 27 38 37 27 28 39 38 28 29 40 39
+ 29 30 41 40 30 31 42 41 31 32 43 42
+ 33 34 45 44 34 35 46 45 35 36 47 46
+ 36 37 48 47 37 49 51 50 49 38 52 51
+ 51 52 54 53 50 51 53 48 38 55 56 52
+ 55 39 57 56 56 57 59 58 52 56 58 54
+ 39 60 61 57 60 40 62 61 61 62 64 63
+ 57 61 63 59 40 41 65 64 41 42 66 65
+ 42 43 67 66 44 45 69 68 45 46 70 69
+ 46 47 71 70 47 72 74 73 72 48 75 74
+ 74 75 77 76 73 74 76 71 48 53 78 75
+ 53 54 79 78 78 79 81 80 75 78 80 77
+ 54 58 82 79 58 59 83 82 82 83 85 84
+ 79 82 84 81 59 63 86 83 63 64 87 86
+ 86 87 89 88 83 86 88 85 64 90 91 87
+ 90 65 92 91 91 92 94 93 87 91 93 89
+ 65 66 95 94 66 67 96 95 68 69 98 97
+ 69 70 99 98 70 71 100 99 71 76 102 101
+ 76 77 103 102 102 103 105 104 101 102 104 100
+ 77 80 106 103 80 81 107 106 106 107 109 108
+ 103 106 108 105 81 84 110 107 84 85 111 110
+ 110 111 113 112 107 110 112 109 85 88 114 111
+ 88 89 115 114 114 115 117 116 111 114 116 113
+ 89 93 118 115 93 94 119 118 118 119 121 120
+ 115 118 120 117 94 95 122 121 95 96 123 122
+ 97 98 125 124 98 99 126 125 99 100 127 126
+ 100 104 129 128 104 105 130 129 129 130 132 131
+ 128 129 131 127 105 108 133 130 108 109 134 133
+ 133 134 136 135 130 133 135 132 109 112 137 134
+ 112 113 138 137 137 138 140 139 134 137 139 136
+ 113 116 141 138 116 117 142 141 141 142 144 143
+ 138 141 143 140 117 120 145 142 120 121 146 145
+ 145 146 148 147 142 145 147 144 121 122 149 148
+ 122 123 150 149 124 125 152 151 125 126 153 152
+ 126 127 154 153 127 132 155 154 132 135 157 156
+ 135 136 158 157 157 158 160 159 156 157 159 155
+ 136 139 161 158 139 140 162 161 161 162 164 163
+ 158 161 163 160 140 143 165 162 143 144 166 165
+ 165 166 168 167 162 165 167 164 144 148 169 168
+ 148 149 170 169 149 150 171 170 151 152 173 172
+ 152 153 174 173 153 154 175 174 154 155 176 175
+ 155 160 177 176 160 164 178 177 164 168 179 178
+ 168 169 180 179 169 170 181 180 170 171 182 181
+ 172 173 184 183 173 174 185 184 174 175 186 185
+ 175 176 187 186 176 177 188 187 177 178 189 188
+ 178 179 190 189 179 180 191 190 180 181 192 191
+ 181 182 193 192
+ </DataArray>
+ <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+ 4 8 12 16 20 24 28 32 36 40 44 48
+ 52 56 60 64 68 72 76 80 84 88 92 96
+ 100 104 108 112 116 120 124 128 132 136 140 144
+ 148 152 156 160 164 168 172 176 180 184 188 192
+ 196 200 204 208 212 216 220 224 228 232 236 240
+ 244 248 252 256 260 264 268 272 276 280 284 288
+ 292 296 300 304 308 312 316 320 324 328 332 336
+ 340 344 348 352 356 360 364 368 372 376 380 384
+ 388 392 396 400 404 408 412 416 420 424 428 432
+ 436 440 444 448 452 456 460 464 468 472 476 480
+ 484 488 492 496 500 504 508 512 516 520 524 528
+ 532 536 540 544 548 552 556 560 564 568 572 576
+ 580 584 588 592 596 600 604 608 612 616 620 624
+ 628 632 636 640 644 648 652
+ </DataArray>
+ <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9
+ </DataArray>
+ </Cells>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>
diff --git a/test/decoupled/2p2c/test_adaptive2p2c.cc b/test/decoupled/2p2c/test_adaptive2p2c.cc
new file mode 100644
index 0000000000..a0ea70c22d
--- /dev/null
+++ b/test/decoupled/2p2c/test_adaptive2p2c.cc
@@ -0,0 +1,57 @@
+/*****************************************************************************
+ * Copyright (C) 2010 by Benjamin Faigle *
+ * Copyright (C) 20010 by Markus Wolff *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+
+#include "test_adaptive2p2cproblem.hh"
+#include <dumux/common/start.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+ "\t-Grid.File Name of the file containing the grid \n"
+ "\t definition in DGF format\n"
+ "\t-SpatialParams.LensLowerLeftX x-coordinate of the lower left corner of the lens [m] \n"
+ "\t-SpatialParams.LensLowerLeftY y-coordinate of the lower left corner of the lens [m] \n"
+ "\t-SpatialParams.LensUpperRightX x-coordinate of the upper right corner of the lens [m] \n"
+ "\t-SpatialParams.LensUpperRightY y-coordinate of the upper right corner of the lens [m] \n"
+ "\n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+// The main function using the standard start procedure
+int main(int argc, char** argv)
+{
+ typedef TTAG(Adaptive2p2c) TypeTag;
+ return Dumux::start<TypeTag>(argc, argv, usage);
+}
diff --git a/test/decoupled/2p2c/test_adaptive2p2c.input b/test/decoupled/2p2c/test_adaptive2p2c.input
new file mode 100644
index 0000000000..88294e507f
--- /dev/null
+++ b/test/decoupled/2p2c/test_adaptive2p2c.input
@@ -0,0 +1,43 @@
+#############################################################
+#Configuration file for test_adaptive2p2c.cc
+#############################################################
+[Impet]
+CFLFactor = 0.8
+ErrorTermLowerBound = 0.2
+ErrorTermUpperBound = 0.9
+ErrorTermFactor = 0.5
+EnableVolumeIntegral = 1
+RestrictFluxInTransport = 1
+
+#############################################################
+#Refinement control
+#############################################################
+[GridAdapt]
+MinLevel = 0
+MaxLevel = 1
+# CoarsenTolerance = 0.1
+# RefineTolerance = 0.15
+EnableMultiPointFluxApproximation = 1
+EnableSecondHalfEdge = 1
+
+#############################################################
+#Problem parameters
+#############################################################
+[Problem]
+
+SimulationName = test_adaptive2p2c # name of the output files
+OutputInterval = 5
+
+[TimeManager]
+#Simulated time
+TEnd= 3e3
+#Initial time step
+DtInitial = 200
+
+[Grid]
+NumberOfCellsX = 10# [-] level 0 resolution in x-direction
+NumberOfCellsY = 10# [-] level 0 resolution in y-direction
+
+UpperRightX = 10# [m] dimension of the grid
+UpperRightY = 10# [m] dimension of the grid # 2D
+
diff --git a/test/decoupled/2p2c/test_adaptive2p2cproblem.hh b/test/decoupled/2p2c/test_adaptive2p2cproblem.hh
new file mode 100644
index 0000000000..3aea9de28f
--- /dev/null
+++ b/test/decoupled/2p2c/test_adaptive2p2cproblem.hh
@@ -0,0 +1,294 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief test problem for the sequential 2p2c model
+ */
+#ifndef DUMUX_TEST_ADAPTIVE_2P2C_PROBLEM_HH
+#define DUMUX_TEST_ADAPTIVE_2P2C_PROBLEM_HH
+
+#if HAVE_ALUGRID
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/sgrid.hh>
+#include <dune/grid/alugrid/2d/alugrid.hh>
+
+#include <dumux/common/cubegridcreator.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/decoupled/2p2c/2p2cadaptiveproperties.hh>
+#include <dumux/decoupled/2p2c/2p2cproblem.hh>
+#include <dumux/decoupled/2p/impes/gridadaptionindicator2p.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+
+#include <dumux/material/fluidsystems/h2oairfluidsystem.hh>
+#include <dumux/material/fluidsystems/h2on2fluidsystem.hh>
+
+#include "test_dec2p2c_spatialparams.hh"
+
+namespace Dumux
+{
+
+template<class TypeTag>
+class Adaptive2p2c;
+
+namespace Properties
+{
+NEW_TYPE_TAG(Adaptive2p2c, INHERITS_FROM(DecoupledTwoPTwoCAdaptive,Test2P2CSpatialParams));
+
+// Set the grid type
+SET_PROP(Adaptive2p2c, Grid)
+{
+ typedef Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming> type;
+// typedef Dune::UGGrid<2> type;
+};
+// set the GridCreator property
+SET_TYPE_PROP(Adaptive2p2c, GridCreator, CubeGridCreator<TypeTag>);
+
+// Set the problem property
+SET_PROP(Adaptive2p2c, Problem)
+{
+ typedef Dumux::Adaptive2p2c<TTAG(Adaptive2p2c)> type;
+};
+
+// Select fluid system
+SET_PROP(Adaptive2p2c, FluidSystem)
+{
+ typedef Dumux::H2OAirFluidSystem<TypeTag> type;
+// typedef Dumux::H2ON2FluidSystem<TypeTag> type;
+// typedef Dumux::Brine_CO2_System<TypeTag, Dumux::Benchmark3::CO2Tables> type;
+};
+
+SET_BOOL_PROP(Adaptive2p2c, EnableComplicatedFluidSystem, true);
+
+// Select water formulation
+SET_PROP(Adaptive2p2c, Components) : public GET_PROP(TypeTag, DefaultComponents)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+// typedef Dumux::TabulatedComponent<Scalar, typename Dumux::H2O<Scalar> > H20;
+ typedef Dumux::H2O<Scalar> H2O;
+};
+// Specify indicator
+SET_TYPE_PROP(Adaptive2p2c, AdaptionIndicator, GridAdaptionIndicator2P<TypeTag>);
+
+// Enable gravity
+SET_BOOL_PROP(Adaptive2p2c, EnableGravity, true);
+SET_INT_PROP(Adaptive2p2c,
+ BoundaryMobility,
+ GET_PROP_TYPE(TypeTag, Indices)::permDependent);
+SET_BOOL_PROP(Adaptive2p2c, EnableCapillarity, true);
+SET_INT_PROP(Adaptive2p2c, PressureFormulation,
+ GET_PROP_TYPE(TypeTag, Indices)::pressureNW);
+
+}
+
+/*!
+ * \ingroup Adaptive2p2cs
+ */
+template<class TypeTag = TTAG(Adaptive2p2c)>
+class Adaptive2p2c: public IMPETProblem2P2C<TypeTag>
+{
+typedef IMPETProblem2P2C<TypeTag> ParentType;
+typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;
+typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+
+
+// boundary typedefs
+typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+
+enum
+{
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+};
+
+enum
+{
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+};
+
+typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+typedef typename GridView::Traits::template Codim<0>::Entity Element;
+typedef typename Grid::Traits::template Codim<0>::EntityPointer ElementPointer;
+typedef typename GridView::Intersection Intersection;
+typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+Adaptive2p2c(TimeManager &timeManager, const GridView& gridView) :
+ ParentType(timeManager, gridView),
+ debugWriter_(gridView, "gridAfterAdapt")
+{
+ this->setGrid(GridCreator::grid());
+ std::string s = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, SimulationName);
+ this->setName(s.c_str());
+ this->setOutputInterval(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Problem, OutputInterval));
+ // initialize the tables of the fluid system
+// WaterFormulation::init(273.15, 623.15, 100,
+// -10, 20e6, 200);
+// FluidSystem::init();
+}
+
+//void preTimeStep()
+//{
+// ParentType::preTimeStep();
+// // use second writer
+// debugWriter_.gridChanged();
+// // write
+// debugWriter_.beginWrite(this->timeManager().time());
+// //write stuff out
+// typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes))::ScalarSolution ScalarSolutionType;
+// typedef typename GET_PROP_TYPE(TypeTag, PTAG(CellData)) CellData;
+// int size = this->gridView().size(0);
+// ScalarSolutionType *pressureW = debugWriter_.allocateManagedBuffer (size);
+// ScalarSolutionType *pressureN = debugWriter_.allocateManagedBuffer (size);
+// ScalarSolutionType *totalConcentration1 = debugWriter_.allocateManagedBuffer (size);
+// ScalarSolutionType *totalConcentration2 = debugWriter_.allocateManagedBuffer (size);
+// for (int i = 0; i < size; i++)
+// {
+// CellData& cellData = this->variables().cellData(i);
+// (*pressureW)[i] = cellData.pressure(wPhaseIdx);
+// (*pressureN)[i] = cellData.pressure(nPhaseIdx);
+// (*totalConcentration1)[i] = cellData.massConcentration(wPhaseIdx);
+// (*totalConcentration2)[i] = cellData.massConcentration(nPhaseIdx);
+// }
+// debugWriter_.attachCellData(*pressureW, "wetting pressure");
+// debugWriter_.attachCellData(*pressureN, "nonwetting pressure");
+// debugWriter_.attachCellData(*totalConcentration1, "C^w from cellData");
+// debugWriter_.attachCellData(*totalConcentration2, "C^n from cellData");
+// debugWriter_.endWrite();
+// return;
+//}
+
+/*!
+ * \name Problem parameters
+ */
+// \{
+
+bool shouldWriteRestartFile() const
+{
+ return false;
+}
+
+//! Returns the temperature within the domain.
+/*! This problem assumes a temperature of 10 degrees Celsius.
+ * \param globalPos The global Position
+ */
+Scalar temperatureAtPos(const GlobalPosition& globalPos) const
+{
+ return 273.15 + 10; // -> 10°C
+}
+
+// \}
+/*!
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::referencePressureAtPos()
+ */
+Scalar referencePressureAtPos(const GlobalPosition& globalPos) const
+{
+ return 1e6;
+}
+/*!
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::boundaryTypesAtPos()
+ */
+void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition& globalPos) const
+{
+ if (globalPos[0] > this->bboxMax()[0]-1E-6 || globalPos[0] < 1e-6)
+ bcTypes.setAllDirichlet();
+ else
+ // all other boundaries
+ bcTypes.setAllNeumann();
+}
+
+/*!
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::boundaryFormulation()
+ */
+const void boundaryFormulation(typename Indices::BoundaryFormulation &bcFormulation, const Intersection& intersection) const
+{
+ bcFormulation = Indices::BoundaryFormulation::concentration;
+}
+/*!
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::dirichletAtPos()
+ */
+void dirichletAtPos(PrimaryVariables &bcValues, const GlobalPosition& globalPos) const
+{
+ Scalar pRef = referencePressureAtPos(globalPos);
+ Scalar temp = temperatureAtPos(globalPos);
+
+ // Dirichlet for pressure equation
+ bcValues[Indices::pressureEqIdx] = (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]);
+
+ // Dirichlet values for transport equations
+ bcValues[Indices::contiWEqIdx] = 1.;
+ bcValues[Indices::contiNEqIdx] = 1.- bcValues[Indices::contiWEqIdx];
+
+}
+/*!
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::neumannAtPos()
+ */
+void neumannAtPos(PrimaryVariables &neumannValues, const GlobalPosition& globalPos) const
+{
+ setZero(neumannValues, Indices::contiWEqIdx);
+}
+/*!
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::sourceAtPos()
+ */
+void source(PrimaryVariables &values, const Element &element)
+{
+ setZero(values, Indices::contiWEqIdx);
+ ElementPointer father(element);
+ // access level 1 entity
+ while (father->level() != this->gridAdapt().getMinLevel())
+ {
+ father = father->father();
+ }
+ GlobalPosition globalPos = father->geometry().center();
+ if (fabs(globalPos[0] - 4.8) < 0.5 && fabs(globalPos[1] - 4.8) < 0.5)
+ values[Indices::contiNEqIdx] = 0.0001;
+}
+/*!
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::initialFormulation()
+ */
+const void initialFormulation(typename Indices::BoundaryFormulation &initialFormulation, const Element& element) const
+{
+ initialFormulation = Indices::concentration;
+}
+/*!
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::initConcentrationAtPos()
+ */
+Scalar initConcentrationAtPos(const GlobalPosition& globalPos) const
+{
+ return 1.0;
+}
+
+private:
+Grid grid_;
+Dumux::VtkMultiWriter<GridView> debugWriter_;
+};
+} //end namespace
+#endif //have alu
+
+#endif
--
GitLab