diff --git a/CHANGELOG b/CHANGELOG
index 7c8a717dc4e192716d85de30d7ad9f79a8ce07e9..6326f4208a86129047313511dc69371eb08351ea 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -1,3 +1,33 @@
+Differences Between DuMuX 2.7 and DuMuX 2.8
+===================================================
+
+* IMPORTANT NOTES:
+
+* IMPROVEMENTS and ENHANCEMENTS:
+ - The zero equation turbulence models (zeroeq) have been added as new models
+ to the freeflow folder.
+ - Tests for coupling a turbulent free flow using zeroeq turbulence models
+ with flow in a porous medium have been added.
+
+* IMMEDIATE INTERFACE CHANGES not allowing/requiring a deprecation period:
+
+* Deprecated PROPERTY and PARAMETER NAMES, to be removed after 2.7: BEWARE: The
+ compiler will not print any warning if a deprecated property or parameter name
+ is used. However, a run-time warning should appear in the summary lines after
+ the corresponding run.
+
+* Deprecated CLASSES/FILES, to be removed after 2.7:
+
+* Deprecated MEMBER FUNCTIONS, to be removed after 2.7:
+
+* Deprecated protected MEMBER VARIABLES, to be removed after 2.7: BEWARE: Older
+ compilers will not print any warning if a deprecated protected member variable
+ is used.
+
+* DELETED classes/files, property names, constants/enums,
+ member functions, which have been deprecated in DuMuX 2.6:
+ Everything listed as deprecated below has been removed.
+
Differences Between DuMuX 2.6 and DuMuX 2.7
===================================================
diff --git a/configure.ac b/configure.ac
index 4302bd7ca4615c382311bf9d6dfb0571f826b09a..df7af336f0fbb53e5b8fd7571a673670496e5f44 100644
--- a/configure.ac
+++ b/configure.ac
@@ -137,6 +137,8 @@ AC_CONFIG_FILES([dumux.pc
test/multidomain/Makefile
test/multidomain/2cstokes2p2c/Makefile
test/multidomain/2cnistokes2p2cni/Makefile
+ test/multidomain/2czeroeq2p2c/Makefile
+ test/multidomain/2cnizeroeq2p2cni/Makefile
test/references/Makefile
tutorial/Makefile
])
diff --git a/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cniproblem.hh b/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cniproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ef4ed40017ca5904f0a6e6fc1017d8018d5c1a1f
--- /dev/null
+++ b/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cniproblem.hh
@@ -0,0 +1,237 @@
+// -*- 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 The problem class for the coupling of a non-isothermal two-component ZeroEq
+ * and a non-isothermal two-phase two-component Darcy model.
+ */
+#ifndef DUMUX_TWOCNIZEROEQTWOPTWOCNIPROBLEM_HH
+#define DUMUX_TWOCNIZEROEQTWOPTWOCNIPROBLEM_HH
+
+#include <dune/grid/multidomaingrid.hh>
+#include <dune/grid/common/gridinfo.hh>
+#include <dune/grid/io/file/dgfparser.hh>
+
+#include <dumux/material/fluidsystems/h2oairfluidsystem.hh>
+#include <dumux/multidomain/common/multidomainproblem.hh>
+#include <dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh>
+#include <dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnipropertydefaults.hh>
+
+#include "2cnizeroeq2p2cnispatialparameters.hh"
+#include "zeroeq2cnisubproblem.hh"
+#include "2p2cnisubproblem.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class TwoCNIZeroEqTwoPTwoCNIProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(TwoCNIZeroEqTwoPTwoCNIProblem, INHERITS_FROM(TwoCNIStokesTwoPTwoCNI));
+
+// Set the grid type
+#if HAVE_UG
+SET_TYPE_PROP(TwoCNIZeroEqTwoPTwoCNIProblem, Grid, Dune::UGGrid<2>);
+#elif HAVE_ALUGRID
+SET_TYPE_PROP(TwoCNIZeroEqTwoPTwoCNIProblem, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>);
+#else
+#error Requires UG or ALUGrid.
+#endif
+
+// Set the global problem
+SET_TYPE_PROP(TwoCNIZeroEqTwoPTwoCNIProblem, Problem, TwoCNIZeroEqTwoPTwoCNIProblem<TypeTag>);
+
+// Set the local coupling operator
+SET_TYPE_PROP(TwoCNIZeroEqTwoPTwoCNIProblem, MultiDomainCouplingLocalOperator,
+ Dumux::TwoCNIStokesTwoPTwoCNILocalOperator<TypeTag>);
+
+// Set the two sub-problems of the global problem
+SET_TYPE_PROP(TwoCNIZeroEqTwoPTwoCNIProblem, SubDomain1TypeTag, TTAG(ZeroEq2cniSubProblem));
+SET_TYPE_PROP(TwoCNIZeroEqTwoPTwoCNIProblem, SubDomain2TypeTag, TTAG(TwoPTwoCNISubProblem));
+
+// Set the global problem in the context of the two sub-problems
+SET_TYPE_PROP(ZeroEq2cniSubProblem, MultiDomainTypeTag, TTAG(TwoCNIZeroEqTwoPTwoCNIProblem));
+SET_TYPE_PROP(TwoPTwoCNISubProblem, MultiDomainTypeTag, TTAG(TwoCNIZeroEqTwoPTwoCNIProblem));
+
+// Set the other sub-problem for each of the two sub-problems
+SET_TYPE_PROP(ZeroEq2cniSubProblem, OtherSubDomainTypeTag, TTAG(TwoPTwoCNISubProblem));
+SET_TYPE_PROP(TwoPTwoCNISubProblem, OtherSubDomainTypeTag, TTAG(ZeroEq2cniSubProblem));
+
+// Set the same spatial parameters for both sub-problems
+SET_TYPE_PROP(ZeroEq2cniSubProblem, SpatialParams, Dumux::TwoCNIZeroEqTwoPTwoCNISpatialParams<TypeTag>);
+SET_TYPE_PROP(TwoPTwoCNISubProblem, SpatialParams, Dumux::TwoCNIZeroEqTwoPTwoCNISpatialParams<TypeTag>);
+
+// Set the fluid system
+SET_PROP(TwoCNIZeroEqTwoPTwoCNIProblem, FluidSystem)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef Dumux::FluidSystems::H2OAir<Scalar,
+ Dumux::H2O<Scalar>,
+ /*useComplexrelations=*/true> type;
+};
+
+// If SuperLU is not available, the UMFPack solver is used:
+#ifdef HAVE_SUPERLU
+SET_TYPE_PROP(TwoCNIZeroEqTwoPTwoCNIProblem, LinearSolver, SuperLUBackend<TypeTag>);
+#else
+SET_TYPE_PROP(TwoCNIZeroEqTwoPTwoCNIProblem, LinearSolver, UMFPackBackend<TypeTag>);
+#endif
+}
+
+/*!
+ * \brief The problem class for the coupling of a non-isothermal two-component ZeroEq (zeroeq2cni)
+ * and a non-isothermal two-phase two-component Darcy model (2p2cni).
+ *
+ * The problem class for the coupling of a non-isothermal two-component ZeroEq (zeroeq2cni)
+ * and a non-isothermal two-phase two-component Darcy model (2p2cni).
+ * It uses the 2p2cniCoupling model and the ZeroEq2cnicoupling model and provides
+ * the problem specifications for common parameters of the two submodels.
+ * The initial and boundary conditions of the submodels are specified in the two subproblems,
+ * 2p2cnisubproblem.hh and zeroeq2cnisubproblem.hh, which are accessible via the coupled problem.
+ */
+template <class TypeTag = TTAG(TwoCNIZeroEqTwoPTwoCNIProblem) >
+class TwoCNIZeroEqTwoPTwoCNIProblem : public MultiDomainProblem<TypeTag>
+{
+ typedef MultiDomainProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGrid) MDGrid;
+ typedef typename MDGrid::LeafGridView MDGridView;
+ enum { dim = MDGridView::dimension };
+ typedef Dune::FieldVector<Scalar, dim> GlobalPosition;
+
+ typedef typename MDGrid::template Codim<0>::LeafIterator ElementIterator;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+public:
+ /*!
+ * \brief docme
+ *
+ * \param hostGrid docme
+ * \param timeManager The TimeManager which is used by the simulation
+ *
+ */
+ TwoCNIZeroEqTwoPTwoCNIProblem(MDGrid &mdGrid,
+ TimeManager &timeManager)
+ : ParentType(mdGrid, timeManager)
+ {
+ dtInit_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
+ episodeLength_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, EpisodeLength);
+
+ // define location of the interface
+ interfacePos_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
+ noDarcyX1_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, NoDarcyX1);
+ noDarcyX2_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, NoDarcyX2);
+
+ // define output options
+ freqRestart_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqRestart);
+ freqOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqOutput);
+
+ zeroeq2cni_ = this->sdID1();
+ twoPtwoCNI_ = this->sdID2();
+
+ initializeGrid();
+
+ // initialize the tables of the fluid system
+ FluidSystem::init();
+// FluidSystem::init(/*tempMin=*/273.15, /*tempMax=*/373.15, /*numTemp=*/200,
+// /*pMin=*/1e4, /*pMax=*/2e5, /*numP=*/200);
+
+ this->timeManager().startNextEpisode(episodeLength_);
+ }
+
+ /*!
+ * \brief Initialization of the grids
+ *
+ * This function splits the multidomain grid in the two
+ * individual subdomain grids and takes care of parallelization.
+ */
+ void initializeGrid()
+ {
+ MDGrid& mdGrid = this->mdGrid();
+ mdGrid.startSubDomainMarking();
+
+ // subdivide grid in two subdomains
+ ElementIterator eendit = mdGrid.template leafend<0>();
+ for (ElementIterator elementIt = mdGrid.template leafbegin<0>();
+ elementIt != eendit; ++elementIt)
+ {
+ // this is required for parallelization checks if element is within a partition
+ if (elementIt->partitionType() != Dune::InteriorEntity)
+ continue;
+
+ GlobalPosition globalPos = elementIt->geometry().center();
+
+ if (globalPos[1] > interfacePos_)
+ mdGrid.addToSubDomain(zeroeq2cni_,*elementIt);
+ else
+ if(globalPos[0] > noDarcyX1_ && globalPos[0] < noDarcyX2_)
+ mdGrid.addToSubDomain(twoPtwoCNI_,*elementIt);
+ }
+ mdGrid.preUpdateSubDomains();
+ mdGrid.updateSubDomains();
+ mdGrid.postUpdateSubDomains();
+
+ gridinfo(this->sdGrid1());
+ gridinfo(this->sdGrid2());
+ }
+
+ //! \copydoc Dumux::CoupledProblem::episodeEnd()
+ void episodeEnd()
+ { this->timeManager().startNextEpisode(episodeLength_); }
+
+ //! \copydoc Dumux::CoupledProblem::shouldWriteRestartFile()
+ bool shouldWriteRestartFile() const
+ {
+ return ( ((this->timeManager().timeStepIndex() > 0)
+ && (this->timeManager().timeStepIndex() % freqRestart_ == 0))
+ // also write a restart file at the end of each episode
+ || this->timeManager().episodeWillBeOver());
+ }
+
+ //! \copydoc Dumux::CoupledProblem::shouldWriteOutput()
+ bool shouldWriteOutput() const
+ {
+ return ( ((this->timeManager().timeStepIndex() > 0)
+ && (this->timeManager().timeStepIndex() % freqOutput_ == 0))
+ // also write a restart file at the end of each episode
+ || this->timeManager().episodeWillBeOver());
+ }
+
+private:
+ typename MDGrid::SubDomainType zeroeq2cni_;
+ typename MDGrid::SubDomainType twoPtwoCNI_;
+
+ unsigned freqRestart_;
+ unsigned freqOutput_;
+
+ Scalar interfacePos_;
+ Scalar noDarcyX1_;
+ Scalar noDarcyX2_;
+ Scalar episodeLength_;
+ Scalar dtInit_;
+};
+
+} //end namespace
+
+#endif // DUMUX_TWOCNIZEROEQTWOPTWOCNIPROBLEM_HH
diff --git a/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cnispatialparameters.hh b/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cnispatialparameters.hh
new file mode 100644
index 0000000000000000000000000000000000000000..20e8c5e408dc3403000c2317432a11ea27a0bdec
--- /dev/null
+++ b/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cnispatialparameters.hh
@@ -0,0 +1,206 @@
+// -*- 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 Spatial parameters for the
+ * coupling of a non-isothermal two-component ZeroEq
+ * and a non-isothermal two-phase two-component Darcy model.
+ */
+#ifndef DUMUX_TWOCNIZEROEQTWOPTWOCNISPATIALPARAMS_HH
+#define DUMUX_TWOCNIZEROEQTWOPTWOCNISPATIALPARAMS_HH
+
+#include <dune/grid/io/file/vtk/common.hh>
+
+#include <dumux/material/spatialparams/implicitspatialparams.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+namespace Dumux
+{
+//forward declaration
+template<class TypeTag>
+class TwoCNIZeroEqTwoPTwoCNISpatialParams;
+
+namespace Properties
+{
+// The spatial parameters TypeTag
+NEW_TYPE_TAG(TwoCNIZeroEqTwoPTwoCNISpatialParams);
+
+// Set the spatial parameters
+SET_TYPE_PROP(TwoCNIZeroEqTwoPTwoCNISpatialParams, SpatialParams,
+ TwoCNIZeroEqTwoPTwoCNISpatialParams<TypeTag>);
+
+// Set the material law parameterized by absolute saturations
+SET_TYPE_PROP(TwoCNIZeroEqTwoPTwoCNISpatialParams,
+ MaterialLaw,
+ EffToAbsLaw<RegularizedVanGenuchten<typename GET_PROP_TYPE(TypeTag, Scalar)>>);
+// EffToAbsLaw<RegularizedBrooksCorey<typename GET_PROP_TYPE(TypeTag, Scalar)> >);
+}
+
+
+/*!
+ * \ingroup TwoPTwoCNIModel
+ * \ingroup ZeroEqTwoCNIModel
+ * \ingroup ImplicitTestProblems
+ * \brief Definition of the spatial parameters for
+ * the coupling of a non-isothermal two-component ZeroEq
+ * and a non-isothermal two-phase two-component Darcy model.
+ */
+template<class TypeTag>
+class TwoCNIZeroEqTwoPTwoCNISpatialParams : public ImplicitSpatialParams<TypeTag>
+{
+ typedef ImplicitSpatialParams<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GridView::ctype CoordScalar;
+
+ enum {
+ dim=GridView::dimension,
+ dimWorld=GridView::dimensionworld
+ };
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+public:
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ /*!
+ * \brief Spatial parameters for the
+ * coupling of a non-isothermal two-component ZeroEq
+ * and a non-isothermal two-phase two-component Darcy model.
+ *
+ * \param gridView The GridView which is used by the problem
+ */
+ TwoCNIZeroEqTwoPTwoCNISpatialParams(const GridView& gridView)
+ : ParentType(gridView)
+ {
+ permeability_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Permeability);
+ porosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Porosity);
+ thermalConductivitySolid_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, ThermalConductivitySolid);
+
+ spatialParams_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Swr));
+ spatialParams_.setSnr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Snr));
+ spatialParams_.setVgAlpha(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, VgAlpha));
+ spatialParams_.setVgn(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, VgN));
+ }
+
+ /*!
+ * \brief Returns the intrinsic permeability tensor \f$[m^2]\f$
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param scvIdx The local index of the sub-control volume
+ */
+ Scalar intrinsicPermeability(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ return permeability_;
+ }
+
+ /*!
+ * \brief Returns the porosity \f$[-]\f$
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param scvIdx The local index of the sub-control volume
+ */
+ Scalar porosity(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ return porosity_;
+ }
+
+ /*!
+ * \brief Returns the parameter object for the material law
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param scvIdx The local index of the sub-control volume
+ */
+ const MaterialLawParams& materialLawParams(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ return spatialParams_;
+ }
+
+ /*!
+ * \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume
+ */
+ Scalar solidHeatCapacity(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ return 790.0;
+ }
+
+ /*!
+ * \brief Returns the density of the solid material (not the bulk density) \f$[kg/m^3]\f$
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param scvIdx The local index of the sub-control volume
+ */
+ Scalar solidDensity(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ return 2650.0;
+ }
+
+ /*!
+ * \brief Returns the thermal conductivity \f$[W/(m*K)]\f$ of the solid
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param scvIdx The local index of the sub-control volume
+ */
+ Scalar solidThermalConductivity(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ return thermalConductivitySolid_;
+ }
+
+private:
+ Scalar permeability_;
+ Scalar porosity_;
+ Scalar thermalConductivitySolid_;
+ MaterialLawParams spatialParams_;
+};
+} // end namespace
+
+#endif // DUMUX_TWOCNIZEROEQTWOPTWOCNISPATIALPARAMS_HH
diff --git a/test/multidomain/2cnizeroeq2p2cni/2p2cnisubproblem.hh b/test/multidomain/2cnizeroeq2p2cni/2p2cnisubproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8412a9258097c0f51b85e4c007c2d88b8d0e7978
--- /dev/null
+++ b/test/multidomain/2cnizeroeq2p2cni/2p2cnisubproblem.hh
@@ -0,0 +1,414 @@
+// -*- 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 Non-isothermal two-phase two-component porous-medium subproblem
+ * with coupling at the top boundary.
+ */
+#ifndef DUMUX_2P2CNISUB_PROBLEM_HH
+#define DUMUX_2P2CNISUB_PROBLEM_HH
+
+#include <dumux/implicit/2p2c/2p2cindices.hh>
+#include <dumux/implicit/common/implicitporousmediaproblem.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivityjohansen.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
+#include <dumux/multidomain/common/subdomainpropertydefaults.hh>
+#include <dumux/multidomain/common/multidomainlocaloperator.hh>
+#include <dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh>
+
+#include "2cnizeroeq2p2cnispatialparameters.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class TwoPTwoCNISubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(TwoPTwoCNISubProblem,
+ INHERITS_FROM(BoxTwoPTwoCNI, SubDomain, TwoCNIZeroEqTwoPTwoCNISpatialParams));
+
+// Set the problem property
+SET_TYPE_PROP(TwoPTwoCNISubProblem, Problem, TwoPTwoCNISubProblem<TTAG(TwoPTwoCNISubProblem)>);
+
+// Use the 2p2cni local jacobian operator for the 2p2cniCoupling model
+SET_TYPE_PROP(TwoPTwoCNISubProblem, LocalResidual, TwoPTwoCNICouplingLocalResidual<TypeTag>);
+
+// Choose pn and Sw as primary variables
+SET_INT_PROP(TwoPTwoCNISubProblem, Formulation, TwoPTwoCFormulation::pnsw);
+
+// The gas component balance (air) is replaced by the total mass balance
+SET_PROP(TwoPTwoCNISubProblem, ReplaceCompEqIdx)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ static const int value = Indices::contiNEqIdx;
+};
+
+// Use the fluid system from the coupled problem
+SET_TYPE_PROP(TwoPTwoCNISubProblem,
+ FluidSystem,
+ typename GET_PROP_TYPE(typename GET_PROP_TYPE(TypeTag, MultiDomainTypeTag), FluidSystem));
+
+// Johanson is used as model to compute the effective thermal heat conductivity
+SET_TYPE_PROP(TwoPTwoCNISubProblem, ThermalConductivityModel,
+ ThermalConductivityJohansen<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Use formulation based on mass fractions
+SET_BOOL_PROP(TwoPTwoCNISubProblem, UseMoles, false);
+
+// Enable/disable velocity output
+SET_BOOL_PROP(TwoPTwoCNISubProblem, VtkAddVelocity, true);
+
+// Enable gravity
+SET_BOOL_PROP(TwoPTwoCNISubProblem, ProblemEnableGravity, true);
+}
+
+/*!
+ * \ingroup ImplicitTestProblems
+ * \ingroup MultidomainProblems
+ * \brief Non-isothermal two-phase two-component porous-medium subproblem
+ * with coupling at the top boundary.
+ *
+ * \todo update description
+ *
+ * This sub problem uses the \ref TwoPTwoCModel. It is part of the 2p2cni model and
+ * is combined with the zeroeq2cnisubproblem for the free flow domain.
+ */
+template <class TypeTag = TTAG(TwoPTwoCNISubProblem) >
+class TwoPTwoCNISubProblem : public ImplicitPorousMediaProblem<TypeTag>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GridView::Grid Grid;
+
+ typedef TwoPTwoCNISubProblem<TypeTag> ThisType;
+ typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
+
+ // copy some indices for convenience
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ // the type tag of the coupled problem
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainTypeTag) CoupledTypeTag;
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+ enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+ enum { // the equation indices
+ contiTotalMassIdx = Indices::contiNEqIdx,
+ contiWEqIdx = Indices::contiWEqIdx,
+ energyEqIdx = Indices::energyEqIdx
+ };
+ enum { // the indices of the primary variables
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx,
+ temperatureIdx = Indices::temperatureIdx
+ };
+ enum {
+ wPhaseOnly = Indices::wPhaseOnly,
+ nPhaseOnly = Indices::nPhaseOnly,
+ bothPhases = Indices::bothPhases
+ };
+ enum {
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx
+ };
+ enum { // grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ /*!
+ * \brief The sub-problem for the porous-medium subdomain
+ *
+ * \param timeManager The TimeManager which is used by the simulation
+ * \param gridView The simulation's idea about physical space
+ */
+ TwoPTwoCNISubProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView)
+ {
+ Scalar noDarcyX1 = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, NoDarcyX1);
+ Scalar noDarcyX2 = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, NoDarcyX2);
+ Scalar xMin = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMin);
+ Scalar xMax = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
+
+ bBoxMin_[0] = std::max(xMin,noDarcyX1);
+ bBoxMax_[0] = std::min(xMax,noDarcyX2);
+ bBoxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMin);
+ bBoxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
+ runUpDistanceX1_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, RunUpDistanceX1); // first part of the interface without coupling
+ runUpDistanceX2_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, RunUpDistanceX2); // second part of the interface without coupling
+
+ refTemperature_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, PorousMedium, RefTemperaturePM);
+ refPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, PorousMedium, RefPressurePM);
+ refSw_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, PorousMedium, RefSw);
+
+ freqMassOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqMassOutput);
+
+ storageLastTimestep_ = Scalar(0);
+ lastMassOutputTime_ = Scalar(0);
+
+ outfile.open("storage.out");
+ outfile << "Time[s]" << ";"
+ << "TotalMassChange[kg/(s*mDepth)]" << ";"
+ << "WaterMassChange[kg/(s*mDepth))]" << ";"
+ << "IntEnergyChange[J/(m^3*s*mDepth)]" << ";"
+ << "WaterMass[kg/mDepth]" << ";"
+ << "WaterMassLoss[kg/mDepth]" << ";"
+ << "EvaporationRate[mm/s]"
+ << std::endl;
+ }
+
+ //! \brief The destructor
+ ~TwoPTwoCNISubProblem()
+ {
+ outfile.close();
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ //! \copydoc Dumux::ImplicitProblem::name()
+ const std::string &name() const
+ { return GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Output, NamePM); }
+
+ //! \copydoc Dumux::ImplicitProblem::init()
+ void init()
+ {
+ ParentType::init();
+ this->model().globalStorage(storageLastTimestep_);
+ }
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ //! \copydoc Dumux::ImplicitProblem::boundaryTypesAtPos()
+ void boundaryTypesAtPos(BoundaryTypes &values,
+ const GlobalPosition &globalPos) const
+ {
+ values.setAllNeumann();
+
+ if (onLowerBoundary_(globalPos))
+ {
+ values.setDirichlet(temperatureIdx, energyEqIdx);
+ }
+
+ if (onUpperBoundary_(globalPos))
+ {
+ if (globalPos[0] > runUpDistanceX1_ - eps_
+ && globalPos[0] < runUpDistanceX2_ + eps_)
+ {
+ values.setAllCouplingInflow();
+ }
+ else
+ values.setAllNeumann();
+ }
+ }
+
+ //! \copydoc Dumux::ImplicitProblem::dirichletAtPos()
+ void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ {
+ initial_(values, globalPos);
+ }
+
+ //! \copydoc Dumux::ImplicitProblem::neumannAtPos()
+ void neumannAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ {
+ values = 0.;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ //! \copydoc Dumux::ImplicitProblem::sourceAtPos()
+ void sourceAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ {
+ values = 0.;
+ }
+
+ //! \copydoc Dumux::ImplicitProblem::initialAtPos()
+ void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ {
+ initial_(values, globalPos);
+ }
+
+ /*!
+ * \brief Return the initial phase state inside a control volume.
+ *
+ * \param vertex The vertex
+ * \param globalIdx The index of the global vertex
+ * \param globalPos The global position
+ */
+ int initialPhasePresence(const Vertex &vertex,
+ const int &globalIdx,
+ const GlobalPosition &globalPos) const
+ {
+ return bothPhases;
+ }
+
+ /*!
+ * \brief Called by the time manager after the time integration to
+ * do some post processing on the solution.
+ */
+ void postTimeStep()
+ {
+ // Calculate masses
+ PrimaryVariables storage;
+
+ this->model().globalStorage(storage);
+ const Scalar time = this->timeManager().time() + this->timeManager().timeStepSize();
+
+ static Scalar initialWaterContent = 0.0; ;
+ if (this->timeManager().time() < this->timeManager().timeStepSize() + 1e-10)
+ initialWaterContent = storage[contiWEqIdx];
+
+ // Write mass balance information for rank 0
+ if (this->gridView().comm().rank() == 0)
+ {
+ if (this->timeManager().timeStepIndex() % freqMassOutput_ == 0
+ || this->timeManager().episodeWillBeOver())
+ {
+ PrimaryVariables storageChange(0.);
+ storageChange = storageLastTimestep_ - storage;
+
+ assert(time - lastMassOutputTime_ != 0);
+ storageChange /= (time - lastMassOutputTime_);
+
+ std::cout << "Time[s]: " << time
+ << " TotalMass[kg]: " << storage[contiTotalMassIdx]
+ << " WaterMass[kg]: " << storage[contiWEqIdx]
+ << " IntEnergy[J/m^3]: " << storage[energyEqIdx]
+ << " WaterMassChange[kg/s]: " << storageChange[contiWEqIdx]
+ << std::endl;
+ if (this->timeManager().time() != 0.)
+ outfile << time << ";"
+ << storageChange[contiTotalMassIdx] << ";"
+ << storageChange[contiWEqIdx] << ";"
+ << storageChange[energyEqIdx] << ";"
+ << storage[contiWEqIdx] << ";"
+ << initialWaterContent - storage[contiWEqIdx] << ";"
+ << storageChange[contiWEqIdx] / (bBoxMax_[0]-bBoxMin_[0])
+ << std::endl;
+
+ storageLastTimestep_ = storage;
+ lastMassOutputTime_ = time;
+ }
+ }
+ }
+
+ /*!
+ * \brief Determine if we are on a corner of the grid
+ *
+ * \param globalPos The global position
+ */
+ bool isCornerPoint(const GlobalPosition &globalPos)
+ {
+ return ((onLeftBoundary_(globalPos) && onLowerBoundary_(globalPos))
+ || (onLeftBoundary_(globalPos) && onUpperBoundary_(globalPos))
+ || (onRightBoundary_(globalPos) && onLowerBoundary_(globalPos))
+ || (onRightBoundary_(globalPos) && onUpperBoundary_(globalPos)));
+ }
+
+ /*!
+ * \brief Returns whether the position is an interface corner point
+ *
+ * This function is required in case of mortar coupling otherwise it should return false
+ *
+ * \param globalPos The global position
+ */
+ bool isInterfaceCornerPoint(const GlobalPosition &globalPos) const
+ { return false; }
+
+ // \}
+
+private:
+ /*!
+ * \brief Internal method for the initial condition
+ * (reused for the dirichlet conditions!)
+ */
+ void initial_(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ values[pressureIdx] = refPressure_
+ + 1000. * this->gravity()[1] * (globalPos[1] - bBoxMax_[1]);
+ values[switchIdx] = refSw_;
+ values[temperatureIdx] = refTemperature_;
+ }
+
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < bBoxMin_[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > bBoxMax_[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < bBoxMin_[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > bBoxMax_[1] - eps_; }
+
+ bool onBoundary_(const GlobalPosition &globalPos) const
+ {
+ return (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)
+ || onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos));
+ }
+
+ static constexpr Scalar eps_ = 1e-8;
+ GlobalPosition bBoxMin_;
+ GlobalPosition bBoxMax_;
+
+ int freqMassOutput_;
+
+ PrimaryVariables storageLastTimestep_;
+ Scalar lastMassOutputTime_;
+
+ Scalar refTemperature_;
+ Scalar refPressure_;
+ Scalar refSw_;
+
+ Scalar runUpDistanceX1_;
+ Scalar runUpDistanceX2_;
+ std::ofstream outfile;
+};
+} //end namespace Dumux
+
+#endif // DUMUX_TWOPTWOCNI_SUBPROBLEM_HH
diff --git a/test/multidomain/2cnizeroeq2p2cni/CMakeLists.txt b/test/multidomain/2cnizeroeq2p2cni/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f6ce5add56ce64d6923f22f9d9f25a4d6521c0d0
--- /dev/null
+++ b/test/multidomain/2cnizeroeq2p2cni/CMakeLists.txt
@@ -0,0 +1,21 @@
+file(COPY evaporationRates.gp DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+
+# test: test_2cnizeroeq2p2cni_1ff
+add_dumux_test(test_2cnizeroeq2p2cni_1ff test_2cnizeroeq2p2cni test_2cnizeroeq2p2cni.cc
+ ${CMAKE_SOURCE_DIR}/bin/runTest.sh
+ ${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
+ ${CMAKE_SOURCE_DIR}/test/references/2cnizeroeq2p2cni-ff-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/zeroeq2cni-00008.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_2cnizeroeq2p2cni
+ -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/test_2cnizeroeq2p2cni_reference.input
+ )
+
+# test: test_2cnizeroeq2p2cni_2pm
+add_dumux_test(test_2cnizeroeq2p2cni_2pm test_2cnizeroeq2p2cni test_2cnizeroeq2p2cni.cc
+ ${CMAKE_SOURCE_DIR}/bin/runTest.sh
+ ${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
+ ${CMAKE_SOURCE_DIR}/test/references/2cnizeroeq2p2cni-pm-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/2p2cni-00008.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_2cnizeroeq2p2cni
+ -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/test_2cnizeroeq2p2cni_reference.input
+ )
diff --git a/test/multidomain/2cnizeroeq2p2cni/Makefile.am b/test/multidomain/2cnizeroeq2p2cni/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..04e515bc8ffad8a74dfdf0ae087c6eb8068d000b
--- /dev/null
+++ b/test/multidomain/2cnizeroeq2p2cni/Makefile.am
@@ -0,0 +1,6 @@
+# tests where program to build and program to run are equal
+check_PROGRAMS = test_2cnizeroeq2p2cni
+
+test_2cnizeroeq2p2cni_SOURCES = test_2cnizeroeq2p2cni.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/multidomain/2cnizeroeq2p2cni/evaporationRates.gp b/test/multidomain/2cnizeroeq2p2cni/evaporationRates.gp
new file mode 100644
index 0000000000000000000000000000000000000000..3b50921549ee3a0bf0932867c513b922b87c3462
--- /dev/null
+++ b/test/multidomain/2cnizeroeq2p2cni/evaporationRates.gp
@@ -0,0 +1,13 @@
+reset
+set datafile separator ';'
+
+set xlabel 'Time [d]'
+set ylabel 'Evaporation rate [mm/d]'
+set xrange [0:5]
+set yrange [0:5]
+plot \
+'storage.out' u ($1/86400):($3*86400) w l lw 2 t 'current'
+
+set terminal pngcairo size 1200,900
+set output 'evaporationRates.png'
+replot
\ No newline at end of file
diff --git a/test/multidomain/2cnizeroeq2p2cni/test_2cnizeroeq2p2cni.cc b/test/multidomain/2cnizeroeq2p2cni/test_2cnizeroeq2p2cni.cc
new file mode 100644
index 0000000000000000000000000000000000000000..82557555a92df481132a441028a3abeddfa5cdce
--- /dev/null
+++ b/test/multidomain/2cnizeroeq2p2cni/test_2cnizeroeq2p2cni.cc
@@ -0,0 +1,342 @@
+// -*- 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 for the coupled non-isothermal two-component ZeroEq and
+ * non-isothermal two-phase two-component Darcy model
+ */
+
+#include "config.h"
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dumux/common/start.hh>
+#include <dumux/io/interfacemeshcreator.hh>
+
+#include "2cnizeroeq2p2cniproblem.hh"
+
+/*!
+ * \brief Print a usage string for simulations.
+ *
+ * \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 printUsage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\nThe list of mandatory options for this program is:\n"
+ "[Grid]\n"
+ "XMin Minumum x-coordinate [m]\n"
+ "XMax Maximum x-coordinate [m]\n"
+ "YMin Minumum y-coordinate [m]\n"
+ "YMax Maximum y-coordinate [m]\n"
+ "CellsX Number of cells in x-direction\n"
+ "CellsY Number of cells in y-direction\n"
+ "GradingY Vertical grading of the cells\n"
+ "RefineTop Specifies whethter the top of the free flow will be refined\n"
+ "InterfacePos Vertical position of the interface [m]\n"
+ "NoDarcyX1 Horizontal position where the porous medium starts [m]\n"
+ "NoDarcyX2 Horizontal position where the porous medium ends [m]\n"
+ "RunUpDistanceX1 Horizontal position where the coupling starts [m]\n"
+ "RunUpDistanceX2 Horizontal position where the coupling ends [m]\n"
+ "\n"
+ "[SpatialParams]\n"
+ "AlphaBJ Beavers-Joseph coefficient [-]\n"
+ "Permeability Hydraulic conductivity [m^2]\n"
+ "Porosity Porosity [-]\n"
+ "Swr Residual water saturation [-]\n"
+ "Snr Residual gas saturation [-]\n"
+ "VgAlpha Van-Genuchten parameter [1/Pa]\n"
+ "VgN Van-Genuchten parameter [-]\n"
+ "ThermalConductivitySolid Thermal conductivity of the solid material [W/(m*K)]\n"
+ "\n"
+ "[FreeFlow]\n"
+ "RefVelocity Inflow velocity [m/s]\n"
+ "RefPressure Reference pressure [Pa]\n"
+ "RefMassfrac Inflow water mass fraction [-]\n"
+ "RefTemperature Inflow temperature [K]\n"
+ "\n"
+ "[PorousMedium]\n"
+ "RefSw Initial water saturation [-]\n"
+ "RefPressurePM Initial pressure [Pa]\n"
+ "RefTemperaturePM Initial temperature [K]\n"
+ "\n"
+ "[Output]\n"
+ "NameFF Name free flow .vtu files\n"
+ "NamePM Name porous medium .vtu files\n"
+ "FreqRestart Frequency of writting restart information\n"
+ "FreqOutput Frequency of writting vtu output\n"
+ "FreqMassOutput Frequency of writting storage output\n"
+ "FreqFluxOutput Frequency of writting flux output\n"
+ "FreqVaporFluxOutput Frequency of writting vapor flux output\n"
+ "\n"
+ "[TimeManager]\n"
+ "EpisodeLength Length of one episode [s]\n"
+ "\n"
+ "[BoundaryLayer]\n"
+ "Model Enable use of boundary layer models (discouraged)\n"
+ "\n"
+ "[MassTransfer]\n"
+ "Model Enable use of mass transfer models (discouraged)\n"
+ "\n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+template <class TypeTag>
+int startLocal_(int argc, char **argv,
+ void (*usage)(const char *, const std::string &))
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGrid) MDGrid;
+ typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+ // print dumux start message
+ Dumux::dumuxMessage_(true);
+
+ ////////////////////////////////////////////////////////////
+ // Load the input parameters
+ ////////////////////////////////////////////////////////////
+
+ typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
+ Dune::ParameterTreeParser::readOptions(argc, argv, ParameterTree::tree());
+
+ if (ParameterTree::tree().hasKey("ParameterFile") || argc==1)
+ {
+ // read input file, but do not overwrite options specified
+ // on the command line, since the latter have precedence.
+ std::string inputFileName ;
+ if(argc==1) // if there are no arguments given (and there is a file ./<programname>.input) we use it as input file
+ {
+ std::cout<< "\nNo parameter file given. \n"
+ << "Defaulting to '"
+ << argv[0]
+ << ".input' for input file.\n";
+ inputFileName = argv[0];
+ inputFileName += ".input";
+ }
+ else
+ inputFileName = GET_RUNTIME_PARAM(TypeTag, std::string, ParameterFile); // otherwise we read from the command line
+
+ std::ifstream parameterFile;
+
+ // check whether the parameter file exists.
+ parameterFile.open(inputFileName.c_str());
+ if (not parameterFile.is_open()){
+ std::cout<< "\n\t -> Could not open file"
+ << inputFileName
+ << ". <- \n\n\n\n";
+ printUsage(argv[0], Dumux::usageTextBlock());
+ return 1;
+ }
+ parameterFile.close();
+
+ Dune::ParameterTreeParser::readINITree(inputFileName,
+ ParameterTree::tree(),
+ /*overwrite=*/false);
+ }
+
+ // initialize MPI, finalize is done automatically on exit
+ static Dune::MPIHelper& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // define the problem dimensions
+ const int dim=2;
+
+ std::string dgfFileName = "temp.dgf";
+ Scalar dt, tEnd;
+ Dune::FieldVector<int, dim> nElements;
+ Scalar interfacePos, gradingFactor;
+ bool refineTop;
+
+ try
+ {
+ dt = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
+ tEnd = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, TEnd);
+
+ nElements[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsX);
+ nElements[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsY);
+ interfacePos = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
+ gradingFactor = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, GradingY);
+ refineTop = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Grid, RefineTop);
+ }
+ catch (Dumux::ParameterException &e) {
+ std::cerr << e << ". Abort!\n";
+ exit(1) ;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ exit(1);
+ }
+ std::cout << "Starting with timestep size = " << dt << "s, simulation end = " << tEnd << "s\n";
+
+ try
+ {
+ Dune::FieldVector<Scalar, dim> min;
+ Dune::FieldVector<Scalar, dim> max;
+ Dune::FieldVector<int, dim> cells(1);
+ min[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMin);
+ min[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMin);
+ max[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
+ max[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMax);
+
+ std::ofstream outfile;
+ outfile.open(dgfFileName);
+ outfile << "DGF" << std::endl;
+ outfile << "Interval" << std::endl;
+ outfile << min[0] << " " << min[1] << " % first corner" << std::endl;
+ outfile << max[0] << " " << max[1] << " % second corner" << std::endl;
+ outfile << cells[0] << " " << cells[1] << " % cells in x and y direction" << std::endl;
+ outfile << "#" << std::endl;
+ outfile << "" << std::endl;
+ outfile << "Cube" << std::endl;
+ outfile << "0 1 2 3" << std::endl;
+ outfile << "#" << std::endl;
+ outfile << "" << std::endl;
+ outfile << "BOUNDARYDOMAIN" << std::endl;
+ outfile << "default 1 % all boundaries have id 1" << std::endl;
+ outfile << "BOUNDARYDOMAIN" << std::endl;
+ outfile << "# unitcube.dgf" << std::endl;
+ outfile.close();
+ }
+ catch (Dumux::ParameterException &e) {
+ std::cerr << e << ". Abort!\n";
+ exit(1) ;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ exit(1);
+ }
+
+ Dumux::InterfaceMeshCreator<Grid> interfaceMeshCreator;
+ GridCreator::gridPtr() = interfaceMeshCreator.create(dgfFileName, nElements, interfacePos, gradingFactor, refineTop);
+
+ if (mpiHelper.size() > 1) {
+ if (!Dune::Capabilities::isParallel<Grid>::v) {
+ std::cerr << "WARNING: THE PROGRAM IS STARTED USING MPI, BUT THE GRID IMPLEMENTATION\n"
+ << " YOU HAVE CHOSEN IS NOT PARALLEL!\n";
+ }
+ GridCreator::loadBalance();
+ }
+
+ // Instantiate the time manager
+ TimeManager timeManager;
+
+ // instantiate grid
+ Dune::shared_ptr<MDGrid> mdGrid_ = Dune::make_shared<MDGrid> (GridCreator::grid());
+
+ // instantiate coupled problem
+ Problem problem(*mdGrid_,
+ timeManager);
+
+ // print all properties and properties
+ bool printProps = false;
+ if (ParameterTree::tree().hasKey("PrintProperties")
+ || ParameterTree::tree().hasKey("TimeManager.PrintProperties"))
+ printProps = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, TimeManager, PrintProperties);
+ if (printProps)
+ {
+ Dumux::Properties::print<TypeTag>();
+ Dumux::Parameters::print<TypeTag>();
+ }
+
+ // deal with the restart stuff
+ bool restart = false;
+ Scalar restartTime = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
+ if (ParameterTree::tree().hasKey("Restart")
+ || ParameterTree::tree().hasKey("TimeManager.Restart")) {
+ restart = true;
+ restartTime = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, Restart);
+ }
+
+ // run the simulation
+ timeManager.init(problem,
+ restartTime, // initial time
+ dt, // initial time step
+ tEnd, // final time
+ restart);
+
+ timeManager.run();
+
+ // print all parameters
+ bool printParams = true;
+ if (ParameterTree::tree().hasKey("PrintParameters")
+ || ParameterTree::tree().hasKey("TimeManager.PrintParameters"))
+ printParams = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, TimeManager, PrintParameters);
+ if (printParams)
+ Dumux::Parameters::print<TypeTag>();
+
+ // print dumux end message
+ Dumux::dumuxMessage_(false);
+
+ return 0;
+}
+
+/*!
+ * \brief Provides a main function which reads in parameters from the
+ * command line and a parameter file.
+ *
+ * In this function only the differentiation between debugger
+ * or not is made.
+ *
+ * \tparam TypeTag The type tag of the problem which needs to be solved
+ *
+ * \param argc The number of command line arguments of the program
+ * \param argv The contents of the command line arguments of the program
+ */
+template <class TypeTag>
+int startLocal(int argc, char **argv,
+ void (*printUsage)(const char *, const std::string &))
+{
+ try {
+ return startLocal_<TypeTag>(argc, argv, printUsage);
+ }
+ catch (Dumux::ParameterException &e)
+ {
+ std::cerr << e << ". Abort!\n";
+ printUsage(argv[0], Dumux::usageTextBlock());
+ return 1;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ return 2;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ return 3;
+ }
+}
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(TwoCNIZeroEqTwoPTwoCNIProblem) ProblemTypeTag;
+ return startLocal<ProblemTypeTag>(argc, argv, printUsage);
+}
diff --git a/test/multidomain/2cnizeroeq2p2cni/test_2cnizeroeq2p2cni.input b/test/multidomain/2cnizeroeq2p2cni/test_2cnizeroeq2p2cni.input
new file mode 100644
index 0000000000000000000000000000000000000000..6ad03a64e9a9b91b52512f1fa65774a54c469fca
--- /dev/null
+++ b/test/multidomain/2cnizeroeq2p2cni/test_2cnizeroeq2p2cni.input
@@ -0,0 +1,90 @@
+[TimeManager]
+DtInitial = 5e-5 # [s]
+MaxTimeStepSize = 360 # [s]
+TEnd = 3600 # [s]
+EpisodeLength = 14400 # [s] # 14400s = 4h
+
+[Grid]
+XMin = 0.0
+XMax = 0.75
+YMin = 0.0
+YMax = 0.75
+# Number of elements in x-, y-direction
+CellsX = 15
+CellsY = 30
+# Grading and refinement of the mesh
+GradingY = 1.16
+RefineTop = true
+# Position information
+NoDarcyX1 = 0.25 # [m] # Beginning of PM below
+NoDarcyX2 = 0.5 # [m] # End of PM below
+RunUpDistanceX1 = 0.251 # [m] # Beginning of without Coupling to PM (x-coordinate)
+RunUpDistanceX2 = 0.499 # [m] # End of without Coupling to PM (x-coordinate)
+InterfacePos = 0.25 # [m] # Vertical position of coupling interface
+
+[Output]
+NameFF = zeroeq2cni
+NamePM = 2p2cni
+# Frequency of restart file, flux and VTK output
+FreqRestart = 5 # how often restart files are written out
+FreqOutput = 5 # 10 # frequency of VTK output
+FreqMassOutput = 5 # 20 # frequency of mass and evaporation rate output (Darcy)
+FreqFluxOutput = 100 # 100 # frequency of detailed flux output
+FreqVaporFluxOutput = 3 # 5 # frequency of summarized flux output
+
+[FreeFlow]
+RefVelocity = 1.0 # [m/s]
+RefPressure = 1e5 # [Pa]
+RefMassfrac = 0.008 # [-]
+RefTemperature = 298.15 # [K]
+
+[BoundaryLayer]
+Model = 0 # disable boundary layer models
+
+[MassTransfer]
+Model = 0 # disable mass transfer models
+
+[PorousMedium]
+RefPressurePM = 1e5 # [Pa]
+RefTemperaturePM = 298.15 # [K]
+RefSw = 0.98 # [-]
+
+[SpatialParams]
+AlphaBJ = 1.0 # [-]
+Permeability = 2.65e-10 # [m^2]
+Porosity = 0.41 # [-]
+Swr = 0.005 # [-]
+Snr = 0.01 # [-]
+VgAlpha = 6.371e-4 # [1/Pa]
+VgN = 8.0 # [-]
+ThermalConductivitySolid = 5.26 # [W/(m*K)]
+
+[Newton]
+MaxRelativeShift = 1e-5
+TargetSteps = 8
+MaxSteps = 12
+WriteConvergence = false
+
+[LinearSolver]
+ResidualReduction = 1e-9
+Verbosity = 0
+MaxIterations = 100
+
+[ZeroEq]
+WriteAllSCVData = -1.0
+# Eddy Viscosity Models
+# 0 = none
+# 1 = Prandtl
+# 2 = modified Van Driest
+# 3 = Baldwin Lomax
+EddyViscosityModel = 2
+# Eddy Diffusivity and Eddy Conductivity Models
+# 0 = none
+# 1 = Reynolds analogy
+# 2 = modified Van Driest
+# 3 = Deissler
+# 4 = Meier and Rotta
+EddyDiffusivityModel = 3
+EddyConductivityModel = 3
+BBoxMinSandGrainRoughness = 0.0 # [m]
+BBoxMaxSandGrainRoughness = 0.0 # [m]
diff --git a/test/multidomain/2cnizeroeq2p2cni/test_2cnizeroeq2p2cni_reference.input b/test/multidomain/2cnizeroeq2p2cni/test_2cnizeroeq2p2cni_reference.input
new file mode 100644
index 0000000000000000000000000000000000000000..2a8cea7706e0d61bc5cab745d92e5c270f3e4765
--- /dev/null
+++ b/test/multidomain/2cnizeroeq2p2cni/test_2cnizeroeq2p2cni_reference.input
@@ -0,0 +1,90 @@
+[TimeManager]
+DtInitial = 1e-4 # [s]
+MaxTimeStepSize = 360 # [s]
+TEnd = 3600 # [s]
+EpisodeLength = 14400 # [s] # 14400s = 4h
+
+[Grid]
+XMin = 0.0
+XMax = 0.75
+YMin = 0.0
+YMax = 0.75
+# Number of elements in x-, y-direction
+CellsX = 15
+CellsY = 15
+# Grading and refinement of the mesh
+GradingY = 1.75
+RefineTop = true
+# Position information
+NoDarcyX1 = 0.25 # [m] # Beginning of PM below
+NoDarcyX2 = 0.5 # [m] # End of PM below
+RunUpDistanceX1 = 0.251 # [m] # Beginning of without Coupling to PM (x-coordinate)
+RunUpDistanceX2 = 0.499 # [m] # End of without Coupling to PM (x-coordinate)
+InterfacePos = 0.25 # [m] # Vertical position of coupling interface
+
+[Output]
+NameFF = zeroeq2cni
+NamePM = 2p2cni
+# Frequency of restart file, flux and VTK output
+FreqRestart = 1000 # how often restart files are written out
+FreqOutput = 5 # 10 # frequency of VTK output
+FreqMassOutput = 1 # 20 # frequency of mass and evaporation rate output (Darcy)
+FreqFluxOutput = 100 # 100 # frequency of detailed flux output
+FreqVaporFluxOutput = 2 # 5 # frequency of summarized flux output
+
+[FreeFlow]
+RefVelocity = 3.5 # [m/s]
+RefPressure = 1e5 # [Pa]
+RefMassfrac = 0.008 # [-]
+RefTemperature = 298.15 # [K]
+
+[BoundaryLayer]
+Model = 0 # disable boundary layer models
+
+[MassTransfer]
+Model = 0 # disable mass transfer models
+
+[PorousMedium]
+RefPressurePM = 1e5 # [Pa]
+RefTemperaturePM = 298.15 # [K]
+RefSw = 0.28 # [-]
+
+[SpatialParams]
+AlphaBJ = 1.0 # [-]
+Permeability = 2.65e-10 # [m^2]
+Porosity = 0.41 # [-]
+Swr = 0.005 # [-]
+Snr = 0.01 # [-]
+VgAlpha = 6.371e-4 # [1/Pa]
+VgN = 8.0 # [-]
+ThermalConductivitySolid = 5.26 # [W/(m*K)]
+
+[Newton]
+MaxRelativeShift = 1e-5
+TargetSteps = 8
+MaxSteps = 12
+WriteConvergence = false
+
+[LinearSolver]
+ResidualReduction = 1e-9
+Verbosity = 0
+MaxIterations = 100
+
+[ZeroEq]
+WriteAllSCVData = -1.0
+# Eddy Viscosity Models
+# 0 = none
+# 1 = Prandtl
+# 2 = modified Van Driest
+# 3 = Baldwin Lomax
+EddyViscosityModel = 2
+# Eddy Diffusivity and Eddy Conductivity Models
+# 0 = none
+# 1 = Reynolds analogy
+# 2 = modified Van Driest
+# 3 = Deissler
+# 4 = Meier and Rotta
+EddyDiffusivityModel = 3
+EddyConductivityModel = 3
+BBoxMinSandGrainRoughness = 0.0 # [m]
+BBoxMaxSandGrainRoughness = 0.0 # [m]
diff --git a/test/multidomain/2cnizeroeq2p2cni/zeroeq2cnisubproblem.hh b/test/multidomain/2cnizeroeq2p2cni/zeroeq2cnisubproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c045774b3ce16c2169a992b25d07aed9da0c8536
--- /dev/null
+++ b/test/multidomain/2cnizeroeq2p2cni/zeroeq2cnisubproblem.hh
@@ -0,0 +1,422 @@
+// -*- 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 Non-isothermal two-component ZeroEq subproblem with air flowing
+ * from the left to the right and coupling at the bottom.
+ */
+#ifndef DUMUX_ZEROEQTWOCNI_SUBPROBLEM_HH
+#define DUMUX_ZEROEQTWOCNI_SUBPROBLEM_HH
+
+#include <dumux/freeflow/zeroeqncni/zeroeqncnimodel.hh>
+#include <dumux/multidomain/common/subdomainpropertydefaults.hh>
+#include <dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh>
+
+#include "2cnizeroeq2p2cnispatialparameters.hh"
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class ZeroEq2cniSubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(ZeroEq2cniSubProblem,
+ INHERITS_FROM(BoxZeroEqncni, SubDomain, TwoCNIZeroEqTwoPTwoCNISpatialParams));
+
+// Set the problem property
+SET_TYPE_PROP(ZeroEq2cniSubProblem, Problem, Dumux::ZeroEq2cniSubProblem<TypeTag>);
+
+// Use the StokesncniCouplingLocalResidual for the computation of the local residual in the ZeroEq domain
+SET_TYPE_PROP(ZeroEq2cniSubProblem, LocalResidual, StokesncniCouplingLocalResidual<TypeTag>);
+
+// Set the property for the material parameters by extracting it from the material law.
+SET_TYPE_PROP(ZeroEq2cniSubProblem,
+ MaterialLawParams,
+ typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params);
+
+// Use the fluid system from the coupled problem
+SET_TYPE_PROP(ZeroEq2cniSubProblem,
+ FluidSystem,
+ typename GET_PROP_TYPE(typename GET_PROP_TYPE(TypeTag, MultiDomainTypeTag), FluidSystem));
+
+// Disable use of mole formulation
+SET_BOOL_PROP(ZeroEq2cniSubProblem, UseMoles, false);
+
+// Disable gravity
+SET_BOOL_PROP(ZeroEq2cniSubProblem, ProblemEnableGravity, false);
+
+// Enable Navier-Stokes
+SET_BOOL_PROP(ZeroEq2cniSubProblem, EnableNavierStokes, true);
+
+// Set the properties for variable inflow BC
+NEW_PROP_TAG(FreeFlowSinusVelocityAmplitude);
+NEW_PROP_TAG(FreeFlowSinusVelocityPeriod);
+SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusVelocityAmplitude, 0.0);
+SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusVelocityPeriod, 3600.0);
+NEW_PROP_TAG(FreeFlowSinusPressureAmplitude);
+NEW_PROP_TAG(FreeFlowSinusPressurePeriod);
+SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusPressureAmplitude, 0.0);
+SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusPressurePeriod, 3600.0);
+NEW_PROP_TAG(FreeFlowSinusConcentrationAmplitude);
+NEW_PROP_TAG(FreeFlowSinusConcentrationPeriod);
+SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusConcentrationAmplitude, 0.0);
+SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusConcentrationPeriod, 3600.0);
+NEW_PROP_TAG(FreeFlowSinusTemperatureAmplitude);
+NEW_PROP_TAG(FreeFlowSinusTemperaturePeriod);
+SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusTemperatureAmplitude, 0.0);
+SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusTemperaturePeriod, 3600.0);
+}
+
+/*!
+ * \ingroup ImplicitTestProblems
+ * \ingroup MultidomainProblems
+ * \brief ZeroEq2cni problem with air flowing from the left to the right.
+ *
+ * \todo update test description
+ * This sub problem uses the \ref ZeroEq2cniModel. It is part of the 2cnizeroeq2p2cni model and
+ * is combined with the 2p2csubproblem for the Darcy domain.
+ */
+template <class TypeTag>
+class ZeroEq2cniSubProblem : public ZeroEqProblem<TypeTag>
+{
+ typedef ZeroEq2cniSubProblem<TypeTag> ThisType;
+ typedef ZeroEqProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+ typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ enum {
+ dim = GridView::dimension
+ };
+ enum { // equation indices
+ massBalanceIdx = Indices::massBalanceIdx,
+ momentumXIdx = Indices::momentumXIdx, // Index of the x-component of the momentum balance
+ momentumYIdx = Indices::momentumYIdx, // Index of the y-component of the momentum balance
+ momentumZIdx = Indices::momentumZIdx, // Index of the z-component of the momentum balance
+ transportEqIdx = Indices::transportEqIdx, // Index of the transport equation (massfraction)
+ energyEqIdx = Indices::energyEqIdx // Index of the energy equation (temperature)
+ };
+ enum { // primary variable indices
+ pressureIdx = Indices::pressureIdx,
+ velocityXIdx = Indices::velocityXIdx,
+ velocityYIdx = Indices::velocityYIdx,
+ velocityZIdx = Indices::velocityZIdx,
+ massOrMoleFracIdx = Indices::massOrMoleFracIdx,
+ temperatureIdx = Indices::temperatureIdx
+ };
+ enum {
+ transportCompIdx = Indices::transportCompIdx, // water component index
+ phaseCompIdx = Indices::phaseCompIdx // air component index
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::ctype CoordScalar;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef Dune::FieldVector<CoordScalar, dim> GlobalPosition;
+
+
+public:
+ /*!
+ * \brief The sub-problem for the ZeroEq subdomain
+ *
+ * \param timeManager The TimeManager which is used by the simulation
+ * \param gridView The simulation's idea about physical space
+ */
+ ZeroEq2cniSubProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView),
+ spatialParams_(gridView)
+ {
+ bBoxMin_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMin);
+ bBoxMax_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
+ bBoxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
+ bBoxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMax);
+ runUpDistanceX1_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, RunUpDistanceX1); // first part of the interface without coupling
+ runUpDistanceX2_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, RunUpDistanceX2); // second part of the interface without coupling
+
+ refVelocity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefVelocity);
+ refPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefPressure);
+ refMassfrac_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefMassfrac);
+ refTemperature_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefTemperature);
+
+ alphaBJ_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, AlphaBJ);
+ }
+
+ // functions have to be overwritten, otherwise they remain uninitialised
+ //! \copydoc BoxProblem::&bBoxMin()
+ const GlobalPosition &bBoxMin() const
+ { return bBoxMin_; }
+
+ //! \copydoc BoxProblem::&bBoxMax()
+ const GlobalPosition &bBoxMax() const
+ { return bBoxMax_; }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ //! \copydoc Dumux::ImplicitProblem::name()
+ const std::string &name() const
+ { return GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Output, NameFF); }
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ //! \copydoc Dumux::ImplicitProblem::boundaryTypesAtPos()
+ void boundaryTypesAtPos(BoundaryTypes &values,
+ const GlobalPosition &globalPos) const
+ {
+ values.setAllDirichlet();
+
+
+ if (onUpperBoundary_(globalPos))
+ {
+ values.setNeumann(transportEqIdx);
+ values.setDirichlet(temperatureIdx, energyEqIdx);
+ }
+
+ if (onLowerBoundary_(globalPos))
+ {
+ values.setNeumann(transportEqIdx);
+ values.setDirichlet(temperatureIdx, energyEqIdx);
+
+ if (globalPos[0] > runUpDistanceX1_ - eps_
+ && globalPos[0] < runUpDistanceX2_ + eps_)
+ {
+ values.setAllCouplingOutflow();
+ }
+ }
+
+ if (onRightBoundary_(globalPos))
+ {
+ values.setAllOutflow();
+
+ if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos)) // corner points
+ values.setAllDirichlet();
+ }
+
+ if (onLeftBoundary_(globalPos))
+ {
+ values.setAllDirichlet();
+ }
+
+ // the mass balance has to be of type outflow
+ // it does not get a coupling condition, since pn is a condition for stokes
+ values.setOutflow(massBalanceIdx);
+
+ if (onRightBoundary_(globalPos))
+ {
+ values.setAllOutflow();
+ values.setDirichlet(pressureIdx, massBalanceIdx);
+ }
+ }
+
+ //! \copydoc Dumux::ImplicitProblem::dirichletAtPos()
+ void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ {
+ values = 0.0;
+
+ values[velocityXIdx] = xVelocity_(globalPos);
+ values[velocityYIdx] = 0.0;
+ values[pressureIdx] = refPressure()
+ + 1.189 * this->gravity()[1] * (globalPos[1] - bBoxMin_[1]);
+ values[massOrMoleFracIdx] = refMassfrac();
+ values[temperatureIdx] = refTemperature();
+ }
+
+ //! \copydoc Dumux::ImplicitProblem::neumannAtPos()
+ void neumannAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ {
+ values = 0.;
+ }
+
+ /*!
+ * \brief Evaluate the Beavers-Joseph coefficient at given position
+ *
+ * \param globalPos The global position
+ */
+ Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
+ {
+ return alphaBJ_;
+ }
+
+ /*!
+ * \brief Evaluate the intrinsic permeability at the corner of a given element
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param scvIdx The local subcontrolvolume index
+ */
+ Scalar permeability(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ return spatialParams_.intrinsicPermeability(element,
+ fvGeometry,
+ scvIdx);
+ }
+
+ //! \copydoc Dumux::ImplicitProblem::sourceAtPos()
+ void sourceAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ values = 0.0;
+ }
+
+ //! \copydoc Dumux::ImplicitProblem::initialAtPos()
+ void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ {
+ initial_(values, globalPos);
+ }
+
+ // \}
+
+ /*!
+ * \brief Determine if we are on a corner of the grid
+ *
+ * \param globalPos The global position
+ */
+ bool isCornerPoint(const GlobalPosition &globalPos)
+ {
+ return ((onLeftBoundary_(globalPos) && onLowerBoundary_(globalPos))
+ || (onLeftBoundary_(globalPos) && onUpperBoundary_(globalPos))
+ || (onRightBoundary_(globalPos) && onLowerBoundary_(globalPos))
+ || (onRightBoundary_(globalPos) && onUpperBoundary_(globalPos)));
+ }
+
+ /*!
+ * \brief Auxiliary function used for the mortar coupling, if mortar coupling,
+ * this should return true
+ *
+ * \param globalPos The global position
+ */
+ bool isInterfaceCornerPoint(const GlobalPosition &globalPos) const
+ {
+ return false;
+ }
+
+ //! \brief Returns the velocity at the inflow.
+ const Scalar refVelocity() const
+ {
+ return refVelocity_ + variation_(GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusVelocityAmplitude),
+ GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusVelocityPeriod));
+ }
+
+ //! \brief Returns the pressure at the inflow.
+ const Scalar refPressure() const
+ {
+ return refPressure_ + variation_(GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusPressureAmplitude),
+ GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusPressurePeriod));
+ }
+
+ //! \brief Returns the mass fraction at the inflow.
+ const Scalar refMassfrac() const
+ {
+ return refMassfrac_ + variation_(GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusConcentrationAmplitude),
+ GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusConcentrationPeriod));
+ }
+
+ //! \brief Returns the temperature at the inflow.
+ const Scalar refTemperature() const
+ {
+ return refTemperature_ + variation_(GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusTemperatureAmplitude),
+ GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusTemperaturePeriod));
+ }
+
+private:
+ // Internal method for the initial and Dirichlet conditions
+ void initial_(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ values[velocityXIdx] = xVelocity_(globalPos);
+ values[velocityYIdx] = 0.;
+
+ values[pressureIdx] = refPressure() + 1.189 * this->gravity()[1] * (globalPos[1] - bBoxMin_[1]);
+ values[massOrMoleFracIdx] = refMassfrac();
+ values[temperatureIdx] = refTemperature();
+ }
+
+ // returns the inflow velocity profile
+ const Scalar xVelocity_(const GlobalPosition &globalPos) const
+ {
+ if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
+ return 0.0;
+ return refVelocity();
+ }
+
+ // can be used for the variation of a boundary condition
+ const Scalar variation_(const Scalar amplitude, const Scalar period) const
+ { return sin(2*M_PI*this->timeManager().time()/period) * amplitude; }
+
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < bBoxMin_[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > bBoxMax_[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < bBoxMin_[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > bBoxMax_[1] - eps_; }
+
+ bool onBoundary_(const GlobalPosition &globalPos) const
+ {
+ return (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)
+ || onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos));
+ }
+
+ SpatialParams spatialParams_;
+
+ static constexpr Scalar eps_ = 1e-8;
+ GlobalPosition bBoxMin_;
+ GlobalPosition bBoxMax_;
+
+ Scalar refVelocity_;
+ Scalar refPressure_;
+ Scalar refMassfrac_;
+ Scalar refTemperature_;
+
+ Scalar alphaBJ_;
+
+ Scalar runUpDistanceX1_;
+ Scalar runUpDistanceX2_;
+};
+} //end namespace
+
+#endif // DUMUX_ZEROEQTWOCNI_SUBPROBLEM_HH
diff --git a/test/multidomain/2czeroeq2p2c/2czeroeq2p2cproblem.hh b/test/multidomain/2czeroeq2p2c/2czeroeq2p2cproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3c34ff8fd98799c5b4f96cfb1b4b06007946d43a
--- /dev/null
+++ b/test/multidomain/2czeroeq2p2c/2czeroeq2p2cproblem.hh
@@ -0,0 +1,239 @@
+// -*- 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 The problem class for the coupling of an isothermal two-component ZeroEq
+ * and an isothermal two-phase two-component Darcy model.
+ */
+#ifndef DUMUX_TWOCZEROEQTWOPTWOCPROBLEM_HH
+#define DUMUX_TWOCZEROEQTWOPTWOCPROBLEM_HH
+
+#include <dune/grid/multidomaingrid.hh>
+#include <dune/grid/common/gridinfo.hh>
+#include <dune/grid/io/file/dgfparser.hh>
+
+#include <dumux/material/fluidsystems/h2oairfluidsystem.hh>
+#include <dumux/multidomain/common/multidomainproblem.hh>
+#include <dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh>
+#include <dumux/multidomain/2cstokes2p2c/2cstokes2p2cpropertydefaults.hh>
+
+#include "2czeroeq2p2cspatialparameters.hh"
+#include "zeroeq2csubproblem.hh"
+#include "2p2csubproblem.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class TwoCZeroEqTwoPTwoCProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(TwoCZeroEqTwoPTwoCProblem, INHERITS_FROM(TwoCStokesTwoPTwoC));
+
+// Set the grid type
+#if HAVE_UG
+SET_TYPE_PROP(TwoCZeroEqTwoPTwoCProblem, Grid, Dune::UGGrid<2>);
+#elif HAVE_ALUGRID
+SET_TYPE_PROP(TwoCZeroEqTwoPTwoCProblem, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>);
+#else
+#error Requires UG or ALUGrid.
+#endif
+
+// Set the global problem
+SET_TYPE_PROP(TwoCZeroEqTwoPTwoCProblem, Problem, TwoCZeroEqTwoPTwoCProblem<TypeTag>);
+
+// Set the two sub-problems of the global problem
+SET_TYPE_PROP(TwoCZeroEqTwoPTwoCProblem, SubDomain1TypeTag, TTAG(ZeroEq2cSubProblem));
+SET_TYPE_PROP(TwoCZeroEqTwoPTwoCProblem, SubDomain2TypeTag, TTAG(TwoPTwoCSubProblem));
+
+// Set the local coupling operator
+SET_TYPE_PROP(TwoCZeroEqTwoPTwoCProblem, MultiDomainCouplingLocalOperator,
+ Dumux::TwoCStokesTwoPTwoCLocalOperator<TypeTag>);
+
+// Set the global problem in the context of the two sub-problems
+SET_TYPE_PROP(ZeroEq2cSubProblem, MultiDomainTypeTag, TTAG(TwoCZeroEqTwoPTwoCProblem));
+SET_TYPE_PROP(TwoPTwoCSubProblem, MultiDomainTypeTag, TTAG(TwoCZeroEqTwoPTwoCProblem));
+
+// Set the other sub-problem for each of the two sub-problems
+SET_TYPE_PROP(ZeroEq2cSubProblem, OtherSubDomainTypeTag, TTAG(TwoPTwoCSubProblem));
+SET_TYPE_PROP(TwoPTwoCSubProblem, OtherSubDomainTypeTag, TTAG(ZeroEq2cSubProblem));
+
+// Set the same spatial parameters for both sub-problems
+SET_TYPE_PROP(ZeroEq2cSubProblem, SpatialParams, Dumux::TwoCZeroEqTwoPTwoCSpatialParams<TypeTag>);
+SET_TYPE_PROP(TwoPTwoCSubProblem, SpatialParams, Dumux::TwoCZeroEqTwoPTwoCSpatialParams<TypeTag>);
+
+// Set the fluid system
+SET_PROP(TwoCZeroEqTwoPTwoCProblem, FluidSystem)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef Dumux::FluidSystems::H2OAir<Scalar,
+ Dumux::H2O<Scalar>,
+ /*useComplexrelations=*/false> type;
+};
+
+// If SuperLU is not available, the UMFPack solver is used:
+#ifdef HAVE_SUPERLU
+SET_TYPE_PROP(TwoCZeroEqTwoPTwoCProblem, LinearSolver, SuperLUBackend<TypeTag>);
+#else
+SET_TYPE_PROP(TwoCZeroEqTwoPTwoCProblem, LinearSolver, UMFPackBackend<TypeTag>);
+#endif
+}
+
+/*!
+ * \brief The problem class for the coupling of an isothermal two-component ZeroEq (zeroeq2c)
+ * and an isothermal two-phase two-component Darcy model (2p2c).
+ *
+ * The problem class for the coupling of a non-isothermal two-component ZeroEq (zeroeq2c)
+ * and an isothermal two-phase two-component Darcy model (2p2c).
+ * It uses the 2p2cCoupling model and the ZeroEq2ccoupling model and provides
+ * the problem specifications for common parameters of the two submodels.
+ * The initial and boundary conditions of the submodels are specified in the two subproblems,
+ * 2p2csubproblem.hh and zeroeq2csubproblem.hh, which are accessible via the coupled problem.
+ */
+template <class TypeTag = TTAG(TwoCZeroEqTwoPTwoCProblem) >
+class TwoCZeroEqTwoPTwoCProblem : public MultiDomainProblem<TypeTag>
+{
+ typedef MultiDomainProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGrid) MDGrid;
+ typedef typename MDGrid::LeafGridView MDGridView;
+ enum { dim = MDGridView::dimension };
+ typedef Dune::FieldVector<Scalar, dim> GlobalPosition;
+
+ typedef typename MDGrid::template Codim<0>::LeafIterator ElementIterator;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+public:
+ /*!
+ * \brief The problem for the coupling of ZeroEq and Darcy flow
+ *
+ * \param mdGrid The multidomain grid
+ * \param timeManager The time manager
+ */
+ TwoCZeroEqTwoPTwoCProblem(MDGrid &mdGrid,
+ TimeManager &timeManager)
+ : ParentType(mdGrid, timeManager)
+ {
+ dtInit_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
+ episodeLength_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, EpisodeLength);
+
+ // define location of the interface
+ interfacePos_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
+ noDarcyX_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, NoDarcyX);
+
+ // define output options
+ freqRestart_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqRestart);
+ freqOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqOutput);
+
+ zeroeq2c_ = this->sdID1();
+ twoPtwoC_ = this->sdID2();
+
+ initializeGrid();
+
+ // initialize the tables of the fluid system
+ FluidSystem::init(/*tempMin=*/273.15, /*tempMax=*/373.15, /*numTemp=*/200,
+ /*pMin=*/1e4, /*pMax=*/2e5, /*numP=*/200);
+
+ this->timeManager().startNextEpisode(episodeLength_);
+ }
+
+ /*!
+ * \brief Initialization of the grids
+ *
+ * This function splits the multidomain grid in the two
+ * individual subdomain grids and takes care of parallelization.
+ */
+ void initializeGrid()
+ {
+ MDGrid& mdGrid = this->mdGrid();
+ mdGrid.startSubDomainMarking();
+
+ // subdivide grid in two subdomains
+ ElementIterator eendit = mdGrid.template leafend<0>();
+ for (ElementIterator elementIt = mdGrid.template leafbegin<0>();
+ elementIt != eendit; ++elementIt)
+ {
+ // this is required for parallelization
+ // checks if element is within a partition
+ if (elementIt->partitionType() != Dune::InteriorEntity)
+ continue;
+
+ GlobalPosition globalPos = elementIt->geometry().center();
+
+ if (globalPos[1] > interfacePos_)
+ mdGrid.addToSubDomain(zeroeq2c_,*elementIt);
+ else
+ if(globalPos[0] > noDarcyX_)
+ mdGrid.addToSubDomain(twoPtwoC_,*elementIt);
+ }
+ mdGrid.preUpdateSubDomains();
+ mdGrid.updateSubDomains();
+ mdGrid.postUpdateSubDomains();
+
+ gridinfo(this->sdGrid1());
+ gridinfo(this->sdGrid2());
+ }
+
+ /*!
+ * \brief Called when the end of an simulation episode is reached.
+ *
+ * Typically a new episode should be started in this method.
+ */
+ void episodeEnd()
+ { this->timeManager().startNextEpisode(episodeLength_); }
+
+ //! \copydoc Dumux::CoupledProblem::shouldWriteRestartFile()
+ bool shouldWriteRestartFile() const
+ {
+ return ( ((this->timeManager().timeStepIndex() > 0)
+ && (this->timeManager().timeStepIndex() % freqRestart_ == 0))
+ // also write a restart file at the end of each episode
+ || this->timeManager().episodeWillBeOver());
+ }
+
+ //! \copydoc Dumux::CoupledProblem::shouldWriteOutput()
+ bool shouldWriteOutput() const
+ {
+ return ( ((this->timeManager().timeStepIndex() > 0)
+ && (this->timeManager().timeStepIndex() % freqOutput_ == 0))
+ // also write a restart file at the end of each episode
+ || this->timeManager().episodeWillBeOver());
+ }
+
+private:
+ typename MDGrid::SubDomainType zeroeq2c_;
+ typename MDGrid::SubDomainType twoPtwoC_;
+
+ unsigned freqRestart_;
+ unsigned freqOutput_;
+
+ Scalar interfacePos_;
+ Scalar noDarcyX_;
+ Scalar episodeLength_;
+ Scalar initializationTime_;
+ Scalar dtInit_;
+};
+
+} // namespace Dumux
+
+#endif // DUMUX_TWOCZEROEQTWOPTWOCPROBLEM_HH
diff --git a/test/multidomain/2czeroeq2p2c/2czeroeq2p2cspatialparameters.hh b/test/multidomain/2czeroeq2p2c/2czeroeq2p2cspatialparameters.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cc3f382933737f0ba45abeeda9ee459e8cab6546
--- /dev/null
+++ b/test/multidomain/2czeroeq2p2c/2czeroeq2p2cspatialparameters.hh
@@ -0,0 +1,160 @@
+// -*- 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 Spatial parameters for the
+ * coupling of an isothermal two-component ZeroEq
+ * and an isothermal two-phase two-component Darcy model.
+ */
+
+#ifndef DUMUX_TWOCZEROEQTWOPTWOCSPATIALPARAMS_HH
+#define DUMUX_TWOCZEROEQTWOPTWOCSPATIALPARAMS_HH
+
+#include <dune/grid/io/file/vtk/common.hh>
+
+#include <dumux/material/spatialparams/implicitspatialparams.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+namespace Dumux
+{
+
+//forward declaration
+template<class TypeTag>
+class TwoCZeroEqTwoPTwoCSpatialParams;
+
+namespace Properties
+{
+// The spatial parameters TypeTag
+NEW_TYPE_TAG(TwoCZeroEqTwoPTwoCSpatialParams);
+
+// Set the spatial parameters
+SET_TYPE_PROP(TwoCZeroEqTwoPTwoCSpatialParams, SpatialParams, TwoCZeroEqTwoPTwoCSpatialParams<TypeTag>);
+
+// Set the material law parameterized by absolute saturations
+SET_TYPE_PROP(TwoCZeroEqTwoPTwoCSpatialParams,
+ MaterialLaw,
+ EffToAbsLaw<RegularizedVanGenuchten<typename GET_PROP_TYPE(TypeTag, Scalar)>>);
+// EffToAbsLaw<RegularizedBrooksCorey<typename GET_PROP_TYPE(TypeTag, Scalar)> >);
+}
+
+
+/*!
+ * \ingroup TwoPTwoCModel
+ * \ingroup ZeroEqModel
+ * \ingroup ImplicitTestProblems
+ * \brief Definition of the spatial parameters for
+ * the coupling of an isothermal two-component ZeroEq
+ * and an isothermal two-phase two-component Darcy model.
+ */
+template<class TypeTag>
+class TwoCZeroEqTwoPTwoCSpatialParams : public ImplicitSpatialParams<TypeTag>
+{
+ typedef ImplicitSpatialParams<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GridView::ctype CoordScalar;
+
+ enum {
+ dim=GridView::dimension,
+ dimWorld=GridView::dimensionworld
+ };
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+public:
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ /*!
+ * \brief Spatial parameters for the
+ * coupling of an isothermal two-component ZeroEq
+ * and an isothermal two-phase two-component Darcy model.
+ *
+ * \param gridView The GridView which is used by the problem
+ */
+ TwoCZeroEqTwoPTwoCSpatialParams(const GridView& gridView)
+ : ParentType(gridView)
+ {
+ permeability_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Permeability);
+ porosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Porosity);
+
+ spatialParams_.setSwr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Swr));
+ spatialParams_.setSnr(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, Snr));
+ spatialParams_.setVgAlpha(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, VgAlpha));
+ spatialParams_.setVgn(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, VgN));
+ }
+
+ /*!
+ * \brief Returns the intrinsic permeability tensor \f$[m^2]\f$
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param scvIdx The local index of the sub-control volume
+ */
+ const Scalar intrinsicPermeability(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ return permeability_;
+ }
+
+ /*!
+ * \brief Returns the porosity \f$[-]\f$
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param scvIdx The local index of the sub-control volume
+ */
+ Scalar porosity(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ return porosity_;
+ }
+
+ /*!
+ * \brief Returns the parameter object for the material law
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param scvIdx The local index of the sub-control volume
+ */
+ const MaterialLawParams& materialLawParams(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ return spatialParams_;
+ }
+
+private:
+ Scalar permeability_;
+ Scalar porosity_;
+ MaterialLawParams spatialParams_;
+};
+
+} // end namespace Dumux
+
+#endif // DUMUX_TWOCZEROEQTWOPTWOCSPATIALPARAMS_HH
diff --git a/test/multidomain/2czeroeq2p2c/2p2csubproblem.hh b/test/multidomain/2czeroeq2p2c/2p2csubproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7c0ec9394df437ab1a2e56fb7612b21e042f5493
--- /dev/null
+++ b/test/multidomain/2czeroeq2p2c/2p2csubproblem.hh
@@ -0,0 +1,399 @@
+// -*- 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 Isothermal two-phase two-component porous-medium subproblem
+ * with coupling at the top boundary.
+ */
+#ifndef DUMUX_TWOPTWOC_SUBPROBLEM_HH
+#define DUMUX_TWOPTWOC_SUBPROBLEM_HH
+
+#include <dumux/implicit/2p2c/2p2cindices.hh>
+#include <dumux/implicit/common/implicitporousmediaproblem.hh>
+#include <dumux/multidomain/common/subdomainpropertydefaults.hh>
+#include <dumux/multidomain/common/multidomainlocaloperator.hh>
+#include <dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh>
+
+#include "2czeroeq2p2cspatialparameters.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class TwoPTwoCSubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(TwoPTwoCSubProblem,
+ INHERITS_FROM(BoxTwoPTwoC, SubDomain, TwoCZeroEqTwoPTwoCSpatialParams));
+
+// Set the problem property
+SET_TYPE_PROP(TwoPTwoCSubProblem, Problem, TwoPTwoCSubProblem<TTAG(TwoPTwoCSubProblem)>);
+
+// Set the local residual extended for the coupling
+SET_TYPE_PROP(TwoPTwoCSubProblem, LocalResidual, TwoPTwoCCouplingLocalResidual<TypeTag>);
+
+// Set pn and Sw as primary variables
+SET_INT_PROP(TwoPTwoCSubProblem, Formulation, TwoPTwoCFormulation::pnsw);
+
+// Set the gas component balance (air) to be replaced by the total mass balance
+SET_PROP(TwoPTwoCSubProblem, ReplaceCompEqIdx)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ static const int value = Indices::contiNEqIdx;
+};
+
+// Used the fluid system from the coupled problem
+SET_TYPE_PROP(TwoPTwoCSubProblem,
+ FluidSystem,
+ typename GET_PROP_TYPE(typename GET_PROP_TYPE(TypeTag, MultiDomainTypeTag), FluidSystem));
+
+// Disable use of mole formulation
+SET_BOOL_PROP(TwoPTwoCSubProblem, UseMoles, false);
+
+// Enable velocity output
+SET_BOOL_PROP(TwoPTwoCSubProblem, VtkAddVelocity, true);
+
+// Enable gravity
+SET_BOOL_PROP(TwoPTwoCSubProblem, ProblemEnableGravity, true);
+}
+
+/*!
+ * \ingroup ImplicitTestProblems
+ * \ingroup MultidomainProblems
+ * \brief Isothermal two-phase two-component porous-medium subproblem
+ * with coupling at the top boundary.
+ *
+ * \todo update description
+ *
+ * This sub problem uses the \ref TwoPTwoCModel. It is part of the 2p2c model and
+ * is combined with the zeroeq2csubproblem for the free flow domain.
+ */
+template <class TypeTag = TTAG(TwoPTwoCSubProblem) >
+class TwoPTwoCSubProblem : public ImplicitPorousMediaProblem<TypeTag>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GridView::Grid Grid;
+
+ typedef TwoPTwoCSubProblem<TypeTag> ThisType;
+ typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
+
+ // copy some indices for convenience
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ // the type tag of the coupled problem
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainTypeTag) CoupledTypeTag;
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq)};
+ enum { // the equation indices
+ contiTotalMassIdx = Indices::contiNEqIdx,
+ contiWEqIdx = Indices::contiWEqIdx,
+ };
+ enum { // the indices of the primary variables
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx,
+ };
+ enum { // the indices for the phase presence
+ wPhaseOnly = Indices::wPhaseOnly,
+ nPhaseOnly = Indices::nPhaseOnly,
+ bothPhases = Indices::bothPhases
+ };
+ enum { // grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ /*!
+ * \brief The sub-problem for the porous-medium subdomain
+ *
+ * \param timeManager The TimeManager which is used by the simulation
+ * \param gridView The simulation's idea about physical space
+ */
+ TwoPTwoCSubProblem(TimeManager &timeManager, const GridView gridView)
+ : ParentType(timeManager, gridView)
+ {
+ Scalar noDarcyX = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, NoDarcyX);
+ Scalar xMin = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMin);
+
+ bBoxMin_[0] = std::max(xMin,noDarcyX);
+ bBoxMax_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
+ bBoxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMin);
+ bBoxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
+ runUpDistanceX_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, RunUpDistanceX); // first part of the interface without coupling
+
+ refTemperature_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, PorousMedium, RefTemperaturePM);
+ refPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, PorousMedium, RefPressurePM);
+ refSw_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, PorousMedium, RefSw);
+
+ freqMassOutput_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Output, FreqMassOutput);
+
+ storageLastTimestep_ = Scalar(0);
+ lastMassOutputTime_ = Scalar(0);
+
+ outfile.open("storage.out");
+ outfile << "Time;"
+ << "TotalMassChange;"
+ << "WaterMassChange;"
+ << "WaterMass"
+ << std::endl;
+ }
+
+ //! \brief The destructor
+ ~TwoPTwoCSubProblem()
+ {
+ outfile.close();
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Returns the problem name
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const std::string &name() const
+ { return GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Output, NamePM); }
+
+ /*!
+ * \brief Called by the Dumux::TimeManager in order to
+ * initialize the problem.
+ *
+ * If you overload this method don't forget to call
+ * ParentType::init()
+ */
+ void init()
+ {
+ ParentType::init();
+ this->model().globalStorage(storageLastTimestep_);
+ }
+
+ /*!
+ * \brief Returns the temperature \f$ K \f$
+ *
+ * \param globalPos The global position
+ */
+ Scalar temperatureAtPos(const GlobalPosition &globalPos) const
+ {
+ return refTemperature_;
+ }
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ //! \copydoc ImplicitProblem::boundaryTypesAtPos()
+ void boundaryTypesAtPos(BoundaryTypes &values,
+ const GlobalPosition &globalPos) const
+ {
+ values.setAllNeumann();
+
+ if (onUpperBoundary_(globalPos))
+ {
+ if (globalPos[0] > runUpDistanceX_ - eps_)
+ values.setAllCouplingInflow();
+ else
+ values.setAllNeumann();
+ }
+ }
+
+ //! \copydoc ImplicitProblem::dirichletAtPos()
+ void dirichletAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ initial_(values, globalPos);
+ }
+
+ //! \copydoc ImplicitProblem::neumannAtPos()
+ void neumannAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ values = 0.;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ //! \copydoc ImplicitProblem::sourceAtPos()
+ void sourceAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ values = Scalar(0);
+ }
+
+ //! \copydoc ImplicitProblem::initialAtPos()
+ void initialAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ initial_(values, globalPos);
+ }
+
+ // \}
+
+ /*!
+ * \brief Return the initial phase state inside a control volume.
+ *
+ * \param vertex The vertex
+ * \param globalIdx The index of the global vertex
+ * \param globalPos The global position
+ */
+ int initialPhasePresence(const Vertex &vertex,
+ const int &globalIdx,
+ const GlobalPosition &globalPos) const
+ {
+ return bothPhases;
+ }
+
+ /*!
+ * \brief Called by the time manager after the time integration to
+ * do some post processing on the solution.
+ */
+ void postTimeStep()
+ {
+ // Calculate masses
+ PrimaryVariables storage;
+
+ this->model().globalStorage(storage);
+ const Scalar time = this->timeManager().time() + this->timeManager().timeStepSize();
+
+ // Write mass balance information for rank 0
+ if (this->gridView().comm().rank() == 0)
+ {
+ if (this->timeManager().timeStepIndex() % freqMassOutput_ == 0
+ || this->timeManager().episodeWillBeOver())
+ {
+ PrimaryVariables storageChange(0.);
+ storageChange = storageLastTimestep_ - storage;
+
+ assert(time - lastMassOutputTime_ != 0);
+ storageChange /= (time - lastMassOutputTime_);
+ // 2d: interface length has to be accounted for
+ // in order to obtain kg/m²s
+ storageChange /= (bBoxMax_[0]-bBoxMin_[0]);
+
+ std::cout << "Time: " << time
+ << " TotalMass: " << storage[contiTotalMassIdx]
+ << " WaterMass: " << storage[contiWEqIdx]
+ << " WaterMassChange: " << storageChange[contiWEqIdx]
+ << std::endl;
+ if (this->timeManager().time() != 0.)
+ outfile << time << ";"
+ << storageChange[contiTotalMassIdx] << ";"
+ << storageChange[contiWEqIdx] << ";"
+ << storage[contiWEqIdx]
+ << std::endl;
+
+ storageLastTimestep_ = storage;
+ lastMassOutputTime_ = time;
+ }
+ }
+ }
+
+ /*!
+ * \brief Determines if globalPos is a corner of the grid
+ *
+ * \param globalPos The global position
+ */
+ bool isCornerPoint(const GlobalPosition &globalPos)
+ {
+ return ((onLeftBoundary_(globalPos) && onLowerBoundary_(globalPos))
+ || (onLeftBoundary_(globalPos) && onUpperBoundary_(globalPos))
+ || (onRightBoundary_(globalPos) && onLowerBoundary_(globalPos))
+ || (onRightBoundary_(globalPos) && onUpperBoundary_(globalPos)));
+ }
+
+ /*!
+ * \brief Returns whether the position is an interface corner point
+ *
+ * This function is required in case of mortar coupling otherwise it should return false
+ *
+ * \param globalPos The global position
+ */
+ bool isInterfaceCornerPoint(const GlobalPosition &globalPos) const
+ { return false; }
+
+private:
+ // Internal method for the initial condition (reused for the dirichlet conditions!)
+ void initial_(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ values[pressureIdx] = refPressure_
+ + 1000. * this->gravity()[1] * (globalPos[1] - bBoxMax_[1]);
+ values[switchIdx] = refSw_;
+ }
+
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < bBoxMin_[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > bBoxMax_[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < bBoxMin_[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > bBoxMax_[1] - eps_; }
+
+ bool onBoundary_(const GlobalPosition &globalPos) const
+ {
+ return (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)
+ || onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos));
+ }
+
+ static constexpr Scalar eps_ = 1e-8;
+ GlobalPosition bBoxMin_;
+ GlobalPosition bBoxMax_;
+ Scalar runUpDistanceX_;
+
+ Scalar refTemperature_;
+ Scalar refPressure_;
+ Scalar refSw_;
+
+ PrimaryVariables storageLastTimestep_;
+ Scalar lastMassOutputTime_;
+ int freqMassOutput_;
+ std::ofstream outfile;
+};
+} //end namespace Dumux
+
+#endif // DUMUX_TWOPTWOC_SUBPROBLEM_HH
diff --git a/test/multidomain/2czeroeq2p2c/CMakeLists.txt b/test/multidomain/2czeroeq2p2c/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..144e6cf8966711c1ba4155784dfb36a411fe0be7
--- /dev/null
+++ b/test/multidomain/2czeroeq2p2c/CMakeLists.txt
@@ -0,0 +1,19 @@
+# test: test_2czeroeq2p2c_1ff
+add_dumux_test(test_2czeroeq2p2c_1ff test_2czeroeq2p2c test_2czeroeq2p2c.cc
+ ${CMAKE_SOURCE_DIR}/bin/runTest.sh
+ ${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
+ ${CMAKE_SOURCE_DIR}/test/references/2czeroeq2p2c-ff-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/zeroeq2c-00026.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_2czeroeq2p2c
+ -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/test_2czeroeq2p2c_reference.input
+ )
+
+# test: test_2czeroeq2p2c_2pm
+add_dumux_test(test_2czeroeq2p2c_2pm test_2czeroeq2p2c test_2czeroeq2p2c.cc
+ ${CMAKE_SOURCE_DIR}/bin/runTest.sh
+ ${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
+ ${CMAKE_SOURCE_DIR}/test/references/2czeroeq2p2c-pm-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/2p2c-00026.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_2czeroeq2p2c
+ -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/test_2czeroeq2p2c_reference.input
+ )
diff --git a/test/multidomain/2czeroeq2p2c/Makefile.am b/test/multidomain/2czeroeq2p2c/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..3b459d4ee7a361e5cbd16377a79dfd5e9d6870f7
--- /dev/null
+++ b/test/multidomain/2czeroeq2p2c/Makefile.am
@@ -0,0 +1,6 @@
+# tests where program to build and program to run are equal
+check_PROGRAMS = test_2czeroeq2p2c
+
+test_2czeroeq2p2c_SOURCES = test_2czeroeq2p2c.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/multidomain/2czeroeq2p2c/test_2czeroeq2p2c.cc b/test/multidomain/2czeroeq2p2c/test_2czeroeq2p2c.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d67b86b6b3eb94789b088de0ad4ecde3b524b297
--- /dev/null
+++ b/test/multidomain/2czeroeq2p2c/test_2czeroeq2p2c.cc
@@ -0,0 +1,330 @@
+// -*- 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 for the coupled isothermal two-component ZeroEq and
+ * isothermal two-phase two-component Darcy model
+ */
+
+#include "config.h"
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dumux/common/start.hh>
+#include <dumux/io/interfacemeshcreator.hh>
+
+#include "2czeroeq2p2cproblem.hh"
+
+/*!
+ * \brief Print a usage string for simulations.
+ *
+ * \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 printUsage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\nThe list of mandatory options for this program is:\n"
+ "[Grid]\n"
+ "XMin Minumum x-coordinate [m]\n"
+ "XMax Maximum x-coordinate [m]\n"
+ "YMin Minumum y-coordinate [m]\n"
+ "YMax Maximum y-coordinate [m]\n"
+ "CellsX Number of cells in x-direction\n"
+ "CellsY Number of cells in y-direction\n"
+ "GradingY Vertical grading of the cells\n"
+ "RefineTop Specifies whethter the top of the free flow will be refined\n"
+ "InterfacePos Vertical position of the interface [m]\n"
+ "NoDarcyX Horizontal position where the porous medium starts [m]\n"
+ "\n"
+ "[SpatialParams]\n"
+ "AlphaBJ Beavers-Joseph coefficient [-]\n"
+ "Permeability Hydraulic conductivity [m^2]\n"
+ "Porosity Porosity [-]\n"
+ "Swr Residual water saturation [-]\n"
+ "Snr Residual gas saturation [-]\n"
+ "VgAlpha Van-Genuchten parameter [1/Pa]\n"
+ "VgN Van-Genuchten parameter [-]\n"
+ "\n"
+ "[FreeFlow]\n"
+ "RefVelocity Inflow velocity [m/s]\n"
+ "RefPressure Reference pressure [Pa]\n"
+ "RefMassfrac Inflow water mass fraction [-]\n"
+ "RefTemperature Inflow temperature [K]\n"
+ "\n"
+ "[PorousMedium]\n"
+ "RefSw Initial water saturation [-]\n"
+ "RefPressurePM Initial pressure [Pa]\n"
+ "RefTemperaturePM Initial temperature [K]\n"
+ "\n"
+ "[Output]\n"
+ "NameFF Name free flow .vtu files\n"
+ "NamePM Name porous medium .vtu files\n"
+ "FreqRestart Frequency of writting restart information\n"
+ "FreqOutput Frequency of writting vtu output\n"
+ "FreqMassOutput Frequency of writting storage output\n"
+ "FreqFluxOutput Frequency of writting flux output\n"
+ "FreqVaporFluxOutput Frequency of writting vapor flux output\n"
+ "\n"
+ "[TimeManager]\n"
+ "EpisodeLength Length of one episode [s]\n"
+ "\n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+template <class TypeTag>
+int startLocal_(int argc, char **argv,
+ void (*usage)(const char *, const std::string &))
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGrid) MDGrid;
+ typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+ // print dumux start message
+ Dumux::dumuxMessage_(true);
+
+ ////////////////////////////////////////////////////////////
+ // Load the input parameters
+ ////////////////////////////////////////////////////////////
+
+ typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
+ Dune::ParameterTreeParser::readOptions(argc, argv, ParameterTree::tree());
+
+ if (ParameterTree::tree().hasKey("ParameterFile") || argc==1)
+ {
+ // read input file, but do not overwrite options specified
+ // on the command line, since the latter have precedence.
+ std::string inputFileName ;
+ if(argc==1) // if there are no arguments given (and there is a file ./<programname>.input) we use it as input file
+ {
+ std::cout<< "\nNo parameter file given. \n"
+ << "Defaulting to '"
+ << argv[0]
+ << ".input' for input file.\n";
+ inputFileName = argv[0];
+ inputFileName += ".input";
+ }
+ else
+ inputFileName = GET_RUNTIME_PARAM(TypeTag, std::string, ParameterFile); // otherwise we read from the command line
+
+ std::ifstream parameterFile;
+
+ // check whether the parameter file exists.
+ parameterFile.open(inputFileName.c_str());
+ if (not parameterFile.is_open()){
+ std::cout<< "\n\t -> Could not open file"
+ << inputFileName
+ << ". <- \n\n\n\n";
+ printUsage(argv[0], Dumux::usageTextBlock());
+ return 1;
+ }
+ parameterFile.close();
+
+ Dune::ParameterTreeParser::readINITree(inputFileName,
+ ParameterTree::tree(),
+ /*overwrite=*/false);
+ }
+
+ // initialize MPI, finalize is done automatically on exit
+ static Dune::MPIHelper& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // define the problem dimensions
+ const int dim=2;
+
+ std::string dgfFileName = "temp.dgf";
+ Scalar dt, tEnd;
+ Dune::FieldVector<int, dim> nElements;
+ Scalar interfacePos, gradingFactor;
+ bool refineTop;
+
+ try {
+ dt = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
+ tEnd = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, TEnd);
+
+ nElements[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsX);
+ nElements[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsY);
+ interfacePos = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
+ gradingFactor = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, GradingY);
+ refineTop = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Grid, RefineTop);
+ }
+ catch (Dumux::ParameterException &e) {
+ std::cerr << e << ". Abort!\n";
+ exit(1) ;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ exit(1);
+ }
+ std::cout << "Starting with timestep size = " << dt << "s, simulation end = " << tEnd << "s\n";
+
+ try
+ {
+ Dune::FieldVector<Scalar, dim> min;
+ Dune::FieldVector<Scalar, dim> max;
+ Dune::FieldVector<int, dim> cells(1);
+ min[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMin);
+ min[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMin);
+ max[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
+ max[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMax);
+
+ std::ofstream outfile;
+ outfile.open(dgfFileName);
+ outfile << "DGF" << std::endl;
+ outfile << "Interval" << std::endl;
+ outfile << min[0] << " " << min[1] << " % first corner" << std::endl;
+ outfile << max[0] << " " << max[1] << " % second corner" << std::endl;
+ outfile << cells[0] << " " << cells[1] << " % cells in x and y direction" << std::endl;
+ outfile << "#" << std::endl;
+ outfile << "" << std::endl;
+ outfile << "Cube" << std::endl;
+ outfile << "0 1 2 3" << std::endl;
+ outfile << "#" << std::endl;
+ outfile << "" << std::endl;
+ outfile << "BOUNDARYDOMAIN" << std::endl;
+ outfile << "default 1 % all boundaries have id 1" << std::endl;
+ outfile << "BOUNDARYDOMAIN" << std::endl;
+ outfile << "# unitcube.dgf" << std::endl;
+ outfile.close();
+ }
+ catch (Dumux::ParameterException &e) {
+ std::cerr << e << ". Abort!\n";
+ exit(1) ;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ exit(1);
+ }
+
+ Dumux::InterfaceMeshCreator<Grid> interfaceMeshCreator;
+ GridCreator::gridPtr() = interfaceMeshCreator.create(dgfFileName, nElements, interfacePos, gradingFactor, refineTop);
+
+ if (mpiHelper.size() > 1) {
+ if (!Dune::Capabilities::isParallel<Grid>::v) {
+ std::cerr << "WARNING: THE PROGRAM IS STARTED USING MPI, BUT THE GRID IMPLEMENTATION\n"
+ << " YOU HAVE CHOSEN IS NOT PARALLEL!\n";
+ }
+ GridCreator::loadBalance();
+ }
+
+ // Instantiate the time manager
+ TimeManager timeManager;
+
+ // instantiate grid
+ Dune::shared_ptr<MDGrid> mdGrid_ = Dune::make_shared<MDGrid> (GridCreator::grid());
+
+ // instantiate coupled problem
+ Problem problem(*mdGrid_, timeManager);
+
+ // print all properties and properties
+ bool printProps = false;
+ if (ParameterTree::tree().hasKey("PrintProperties")
+ || ParameterTree::tree().hasKey("TimeManager.PrintProperties"))
+ printProps = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, TimeManager, PrintProperties);
+ if (printProps)
+ {
+ Dumux::Properties::print<TypeTag>();
+ Dumux::Parameters::print<TypeTag>();
+ }
+
+ // deal with the restart stuff
+ bool restart = false;
+ Scalar restartTime = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
+ if (ParameterTree::tree().hasKey("Restart")
+ || ParameterTree::tree().hasKey("TimeManager.Restart")) {
+ restart = true;
+ restartTime = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, Restart);
+ }
+
+ // run the simulation
+ timeManager.init(problem,
+ restartTime, // initial time
+ dt, // initial time step
+ tEnd, // final time
+ restart);
+
+ timeManager.run();
+
+ // print all parameters
+ bool printParams = true;
+ if (ParameterTree::tree().hasKey("PrintParameters")
+ || ParameterTree::tree().hasKey("TimeManager.PrintParameters"))
+ printParams = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, TimeManager, PrintParameters);
+ if (printParams)
+ Dumux::Parameters::print<TypeTag>();
+
+ // print dumux end message
+ Dumux::dumuxMessage_(false);
+
+ return 0;
+}
+
+/*!
+ * \brief Provides a main function which reads in parameters from the
+ * command line and a parameter file.
+ *
+ * In this function only the differentiation between debugger
+ * or not is made.
+ *
+ * \tparam TypeTag The type tag of the problem which needs to be solved
+ *
+ * \param argc The number of command line arguments of the program
+ * \param argv The contents of the command line arguments of the program
+ */
+template <class TypeTag>
+int startLocal(int argc, char **argv,
+ void (*printUsage)(const char *, const std::string &))
+{
+ try {
+ return startLocal_<TypeTag>(argc, argv, printUsage);
+ }
+ catch (Dumux::ParameterException &e)
+ {
+ std::cerr << e << ". Abort!\n";
+ printUsage(argv[0], Dumux::usageTextBlock());
+ return 1;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ return 2;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ return 3;
+ }
+}
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(TwoCZeroEqTwoPTwoCProblem) ProblemTypeTag;
+ return startLocal<ProblemTypeTag>(argc, argv, printUsage);
+}
diff --git a/test/multidomain/2czeroeq2p2c/test_2czeroeq2p2c.input b/test/multidomain/2czeroeq2p2c/test_2czeroeq2p2c.input
new file mode 100644
index 0000000000000000000000000000000000000000..d896a91fd1e268490eed9b6ad4dcf3deb860d698
--- /dev/null
+++ b/test/multidomain/2czeroeq2p2c/test_2czeroeq2p2c.input
@@ -0,0 +1,83 @@
+[TimeManager]
+DtInitial = 5e-5 # [s]
+MaxTimeStepSize = 900 # [s]
+TEnd = 518400 # [s] # 518400 = 6 days
+EpisodeLength = 8e5 # [s] # 43200
+
+[Grid]
+XMin = 0.0 # [m]
+XMax = 0.5 # [m]
+YMin = 0.0 # [m]
+YMax = 0.75 # [m]
+# Number of elements in x-, y-direction
+CellsX = 20
+CellsY = 40
+# Grading and refinement of the mesh
+GradingY = 1.16
+Refinement = 0
+RefineTop = false
+# Position information
+NoDarcyX = 0.25 # [m] # Horizontal position without PM below
+RunUpDistanceX = 0.26 # [m] # Horizontal position without Coupling to PM
+InterfacePos = 0.25 # [m] # Vertical position of coupling interface
+
+[Output]
+NameFF = zeroeq2c
+NamePM = 2p2c
+# Frequency of restart file, flux and VTK output
+FreqRestart = 50 # 500 # how often restart files are written out
+FreqOutput = 5 # 10 # frequency of VTK output
+FreqMassOutput = 5 # 20 # frequency of mass and evaporation rate output (Darcy)
+FreqFluxOutput = 3 # 100 # frequency of detailed flux output
+FreqVaporFluxOutput = 3 # 5 # frequency of summarized flux output
+
+[FreeFlow]
+RefVelocity = 3.5 # [m/s]
+RefPressure = 1e5 # [Pa]
+RefMassfrac = 0.008 # [-]
+RefTemperature = 298.15 # [K]
+SlipVelocity = 0.0 # [m/s]
+
+[PorousMedium]
+RefPressurePM = 1e5 # [Pa]
+RefTemperaturePM = 298.15 # [K]
+RefSw = 0.9 # [-]
+
+[SpatialParams]
+AlphaBJ = 1.0 # [-]
+Permeability = 2.65e-10 # [m^2]
+Porosity = 0.41 # [-]
+Swr = 0.005 # [-]
+Snr = 0.01 # [-]
+VgAlpha = 6.371e-4 # [1/Pa]
+VgN = 6.9 # [-]
+
+# optional parameters
+[Newton]
+MaxRelativeShift = 1e-5
+TargetSteps = 10
+MaxSteps = 15
+WriteConvergence = false
+
+[LinearSolver]
+ResidualReduction = 1e-8
+Verbosity = 0
+MaxIterations = 100
+
+[ZeroEq]
+WriteAllSCVData = -1.0
+# Eddy Viscosity Models
+# 0 = none
+# 1 = Prandtl
+# 2 = modified Van Driest
+# 3 = Baldwin Lomax
+EddyViscosityModel = 3
+# Eddy Diffusivity Models
+# 0 = none
+# 1 = Reynolds analogy
+# 2 = modified Van Driest
+# 3 = Deissler
+# 4 = Meier and Rotta
+EddyDiffusivityModel = 3
+BBoxMinSandGrainRoughness = 0.0 # [m] # 0.6e-3
+BBoxMaxSandGrainRoughness = 0.0 # [m] # 0
diff --git a/test/multidomain/2czeroeq2p2c/test_2czeroeq2p2c_reference.input b/test/multidomain/2czeroeq2p2c/test_2czeroeq2p2c_reference.input
new file mode 100644
index 0000000000000000000000000000000000000000..93bfb438070bed44a522b19d88380d887493010a
--- /dev/null
+++ b/test/multidomain/2czeroeq2p2c/test_2czeroeq2p2c_reference.input
@@ -0,0 +1,82 @@
+[TimeManager]
+DtInitial = 5e-5 # [s]
+MaxTimeStepSize = 900 # [s]
+TEnd = 86400 # [s]
+EpisodeLength = 8e5 # [s]
+
+[Grid]
+XMin = 0.0
+XMax = 0.5
+YMin = 0.0
+YMax = 0.75
+# Number of elements in x-, y-direction
+CellsX = 10
+CellsY = 25
+# Grading and refinement of the mesh
+GradingY = 1.1
+RefineTop = false
+# Position information
+NoDarcyX = 0.25 # [m] # Horizontal position without PM below
+RunUpDistanceX = 0.26 # [m] # Horizontal position without Coupling to PM
+InterfacePos = 0.25 # [m] # Vertical position of coupling interface
+
+[Output]
+NameFF = zeroeq2c
+NamePM = 2p2c
+# Frequency of restart file, flux and VTK output
+FreqRestart = 100 # 500 # how often restart files are written out
+FreqOutput = 5 # 10 # frequency of VTK output
+FreqMassOutput = 5 # 20 # frequency of mass and evaporation rate output (Darcy)
+FreqFluxOutput = 3 # 100 # frequency of detailed flux output
+FreqVaporFluxOutput = 3 # 5 # frequency of summarized flux output
+
+[FreeFlow]
+RefVelocity = 3.5 # [m/s]
+RefPressure = 1e5 # [Pa]
+RefMassfrac = 0.008 # [-]
+RefTemperature = 298.15 # [K]
+
+[PorousMedium]
+RefPressurePM = 1e5 # [Pa]
+RefTemperaturePM = 298.15 # [K]
+RefSw = 0.9 # [-]
+
+[SpatialParams]
+AlphaBJ = 1.0 # [-]
+Permeability = 2.65e-10 # [m^2]
+Porosity = 0.41 # [-]
+Swr = 0.005 # [-]
+Snr = 0.01 # [-]
+VgAlpha = 6.371e-4 # [1/Pa]
+VgN = 6.9 # [-]
+
+# optional parameters
+[Newton]
+MaxRelativeShift = 1e-5
+TargetSteps = 10
+MaxSteps = 15
+WriteConvergence = false
+MaxTimeStepDivisions = 20
+
+[LinearSolver]
+ResidualReduction = 1e-8
+Verbosity = 0
+MaxIterations = 100
+
+[ZeroEq]
+WriteAllSCVData = -1.0
+# Eddy Viscosity Models
+# 0 = none
+# 1 = Prandtl
+# 2 = modified Van Driest
+# 3 = Baldwin Lomax
+EddyViscosityModel = 3
+# Eddy Diffusivity Models
+# 0 = none
+# 1 = Reynolds analogy
+# 2 = modified Van Driest
+# 3 = Deissler
+# 4 = Meier and Rotta
+EddyDiffusivityModel = 3
+BBoxMinSandGrainRoughness = 0.0 # [m] # 0.6e-3
+BBoxMaxSandGrainRoughness = 0.0 # [m] # 0
diff --git a/test/multidomain/2czeroeq2p2c/zeroeq2csubproblem.hh b/test/multidomain/2czeroeq2p2c/zeroeq2csubproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..47471c7a993a79fff91879876fa8bc340d702cea
--- /dev/null
+++ b/test/multidomain/2czeroeq2p2c/zeroeq2csubproblem.hh
@@ -0,0 +1,466 @@
+// -*- 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 Isothermal two-component ZeroEq subproblem with air flowing
+ * from the left to the right and coupling at the bottom.
+ */
+#ifndef DUMUX_ZEROEQTWOCSUBPROBLEM_HH
+#define DUMUX_ZEROEQTWOCSUBPROBLEM_HH
+
+#include <dumux/freeflow/zeroeqnc/zeroeqncmodel.hh>
+#include <dumux/multidomain/common/subdomainpropertydefaults.hh>
+#include <dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh>
+
+#include "2czeroeq2p2cspatialparameters.hh"
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class ZeroEq2cSubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(ZeroEq2cSubProblem,
+ INHERITS_FROM(BoxZeroEqnc, SubDomain, TwoCZeroEqTwoPTwoCSpatialParams));
+
+// Set the problem property
+SET_TYPE_PROP(ZeroEq2cSubProblem, Problem, Dumux::ZeroEq2cSubProblem<TypeTag>);
+
+// Set the property for the material parameters by extracting it from the material law.
+SET_TYPE_PROP(ZeroEq2cSubProblem,
+ MaterialLawParams,
+ typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params);
+
+// Use the StokencCouplingLocalResidual for the computation of the local residual in the ZeroEq domain
+SET_TYPE_PROP(ZeroEq2cSubProblem, LocalResidual,
+ StokesncCouplingLocalResidual<TypeTag>);
+
+// Used the fluid system from the coupled problem
+SET_TYPE_PROP(ZeroEq2cSubProblem,
+ FluidSystem,
+ typename GET_PROP_TYPE(typename GET_PROP_TYPE(TypeTag, MultiDomainTypeTag), FluidSystem));
+
+// Disable use of mole formulation
+SET_BOOL_PROP(ZeroEq2cSubProblem, UseMoles, false);
+
+// Disable gravity
+SET_BOOL_PROP(ZeroEq2cSubProblem, ProblemEnableGravity, false);
+
+// Enable Navier-Stokes
+SET_BOOL_PROP(ZeroEq2cSubProblem, EnableNavierStokes, true);
+
+// Set the properties for variable inflow BC
+NEW_PROP_TAG(FreeFlowSinusVelocityAmplitude);
+NEW_PROP_TAG(FreeFlowSinusVelocityPeriod);
+SET_SCALAR_PROP(ZeroEq2cSubProblem, FreeFlowSinusVelocityAmplitude, 0.0);
+SET_SCALAR_PROP(ZeroEq2cSubProblem, FreeFlowSinusVelocityPeriod, 3600.0);
+NEW_PROP_TAG(FreeFlowSinusPressureAmplitude);
+NEW_PROP_TAG(FreeFlowSinusPressurePeriod);
+SET_SCALAR_PROP(ZeroEq2cSubProblem, FreeFlowSinusPressureAmplitude, 0.0);
+SET_SCALAR_PROP(ZeroEq2cSubProblem, FreeFlowSinusPressurePeriod, 3600.0);
+NEW_PROP_TAG(FreeFlowSinusConcentrationAmplitude);
+NEW_PROP_TAG(FreeFlowSinusConcentrationPeriod);
+SET_SCALAR_PROP(ZeroEq2cSubProblem, FreeFlowSinusConcentrationAmplitude, 0.0);
+SET_SCALAR_PROP(ZeroEq2cSubProblem, FreeFlowSinusConcentrationPeriod, 3600.0);
+NEW_PROP_TAG(FreeFlowSinusTemperatureAmplitude);
+NEW_PROP_TAG(FreeFlowSinusTemperaturePeriod);
+SET_SCALAR_PROP(ZeroEq2cSubProblem, FreeFlowSinusTemperatureAmplitude, 0.0);
+SET_SCALAR_PROP(ZeroEq2cSubProblem, FreeFlowSinusTemperaturePeriod, 3600.0);
+}
+
+/*!
+ * \ingroup ImplicitTestProblems
+ * \ingroup MultidomainProblems
+ * \brief ZeroEq2c problem with air flowing from the left to the right.
+ *
+ * \todo update test description
+ * The zeroeq subdomain is sized 1m times 1m. The boundary conditions for the momentum balances
+ * are all set to Dirichlet. The mass balance receives
+ * outflow bcs, which are replaced in the localresidual by the sum
+ * of the two momentum balances. In the middle of the right boundary,
+ * one vertex receives Dirichlet bcs, to set the pressure level.
+ *
+ * This sub problem uses the \ref ZeroEqTwoCModel. It is part of the 2czeroeq2p2c model and
+ * is combined with the 2p2csubproblem for the Darcy domain.
+ */
+template <class TypeTag>
+class ZeroEq2cSubProblem : public ZeroEqProblem<TypeTag>
+{
+ typedef ZeroEq2cSubProblem<TypeTag> ThisType;
+ typedef ZeroEqProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+ typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ enum {
+ dim = GridView::dimension
+ };
+ enum { // equation indices
+ massBalanceIdx = Indices::massBalanceIdx,
+ momentumXIdx = Indices::momentumXIdx, // Index of the x-component of the momentum balance
+ momentumYIdx = Indices::momentumYIdx, // Index of the y-component of the momentum balance
+ momentumZIdx = Indices::momentumZIdx, // Index of the z-component of the momentum balance
+ transportEqIdx = Indices::transportEqIdx // Index of the transport equation (massfraction)
+ };
+ enum { // primary variable indices
+ pressureIdx = Indices::pressureIdx,
+ velocityXIdx = Indices::velocityXIdx,
+ velocityYIdx = Indices::velocityYIdx,
+ velocityZIdx = Indices::velocityZIdx,
+ massOrMoleFracIdx = Indices::massOrMoleFracIdx
+ };
+ enum { phaseIdx = Indices::phaseIdx };
+ enum { numComponents = Indices::numComponents };
+ enum {
+ transportCompIdx = Indices::transportCompIdx, // water component index
+ phaseCompIdx = Indices::phaseCompIdx // air component index
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::ctype CoordScalar;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef Dune::FieldVector<CoordScalar, dim> GlobalPosition;
+
+
+public:
+ /*!
+ * \brief The sub-problem for the ZeroEq subdomain
+ *
+ * \param timeManager The TimeManager which is used by the simulation
+ * \param gridView The simulation's idea about physical space
+ */
+ ZeroEq2cSubProblem(TimeManager &timeManager, const GridView gridView)
+ : ParentType(timeManager, gridView),
+ spatialParams_(gridView)
+ {
+ bBoxMin_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMin);
+ bBoxMax_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, XMax);
+ bBoxMin_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
+ bBoxMax_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, YMax);
+ runUpDistanceX_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, RunUpDistanceX); // first part of the interface without coupling
+
+ refVelocity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefVelocity);
+ refPressure_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefPressure);
+ refMassfrac_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefMassfrac);
+ refTemperature_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefTemperature);
+ alphaBJ_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, AlphaBJ);
+ }
+
+ // functions have to be overwritten, otherwise they remain uninitialised
+ //! \copydoc BoxProblem::&bBoxMin()
+ const GlobalPosition &bBoxMin() const
+ { return bBoxMin_; }
+
+ //! \copydoc BoxProblem::&bBoxMax()
+ const GlobalPosition &bBoxMax() const
+ { return bBoxMax_; }
+
+ /*!
+ * \brief Returns the problem name
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const std::string &name() const
+ { return GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Output, NameFF); }
+
+ /*!
+ * \brief Returns the temperature \f$ K \f$
+ *
+ * \param globalPos The global position
+ */
+ Scalar temperatureAtPos(const GlobalPosition &globalPos) const
+ {
+ return refTemperature_;
+ };
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ //! \copydoc ImplicitProblem::boundaryTypesAtPos()
+ void boundaryTypesAtPos(BoundaryTypes &values,
+ const GlobalPosition &globalPos) const
+ {
+ values.setAllDirichlet();
+
+ if (onUpperBoundary_(globalPos))
+ {
+ values.setNeumann(transportEqIdx);
+ }
+
+ if (onRightBoundary_(globalPos))
+ {
+ values.setAllOutflow();
+
+ if (onUpperBoundary_(globalPos)) // corner point
+ values.setAllDirichlet();
+ }
+
+ if (onLowerBoundary_(globalPos))
+ {
+ values.setAllDirichlet();
+ values.setNeumann(transportEqIdx);
+
+ if (globalPos[0] > runUpDistanceX_-eps_)
+ values.setAllCouplingOutflow();
+ }
+ if (onLeftBoundary_(globalPos))
+ {
+ // Left inflow boundaries should be Neumann, otherwise the
+ // evaporative fluxes are much more grid dependent
+ values.setNeumann(transportEqIdx);
+
+ if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos)) // corner point
+ values.setAllDirichlet();
+ }
+
+ // the mass balance has to be of type outflow
+ // it does not get a coupling condition, since pn is a condition for ZeroEq
+ values.setOutflow(massBalanceIdx);
+
+ if (onRightBoundary_(globalPos))
+ values.setDirichlet(pressureIdx, massBalanceIdx);
+ }
+
+ //! \copydoc ImplicitProblem::dirichletAtPos()
+ void dirichletAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ initial_(values, globalPos);
+ }
+
+ //! \copydoc ImplicitProblem::neumannAtPos()
+ void neumannAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ values = 0.0;
+
+ FluidState fluidState;
+ updateFluidStateForBC_(fluidState);
+ const Scalar density = FluidSystem::density(fluidState, phaseIdx);
+
+ // rho*v*X at inflow
+ if (onLeftBoundary_(globalPos)
+ && globalPos[1] > bBoxMin_[1] && globalPos[1] < bBoxMax_[1])
+ {
+ values[transportEqIdx] = -xVelocity_(globalPos) * density * refMassfrac();
+ }
+ }
+
+ /*!
+ * \brief Evaluate the Beavers-Joseph coefficient at given position
+ *
+ * \param globalPos The global position
+ *
+ * \return Beavers-Joseph coefficient
+ */
+ Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
+ {
+ return alphaBJ_;
+ }
+
+ /*!
+ * \brief Returns the intrinsic permeability tensor \f$[m^2]\f$
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param scvIdx The local index of the sub-control volume
+ */
+ Scalar permeability(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ return spatialParams_.intrinsicPermeability(element,
+ fvGeometry,
+ scvIdx);
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ //! \copydoc ImplicitProblem::sourceAtPos()
+ void sourceAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ values = Scalar(0);
+ }
+
+ //! \copydoc ImplicitProblem::initialAtPos()
+ void initialAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ initial_(values, globalPos);
+ }
+
+ // \}
+
+ /*!
+ * \brief Determines if globalPos is a corner of the grid
+ *
+ * \param globalPos The global position
+ */
+ bool isCornerPoint(const GlobalPosition &globalPos)
+ {
+ return ((onLeftBoundary_(globalPos) && onLowerBoundary_(globalPos))
+ || (onLeftBoundary_(globalPos) && onUpperBoundary_(globalPos))
+ || (onRightBoundary_(globalPos) && onLowerBoundary_(globalPos))
+ || (onRightBoundary_(globalPos) && onUpperBoundary_(globalPos)));
+ }
+
+ /*!
+ * \brief Auxiliary function used for the mortar coupling, if mortar coupling,
+ * this should return true
+ *
+ * \param globalPos The global position
+ */
+ bool isInterfaceCornerPoint(const GlobalPosition &globalPos) const
+ { return false; }
+
+ //! \brief Returns the velocity at the inflow.
+ const Scalar refVelocity() const
+ {
+ return refVelocity_ + variation_(GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusVelocityAmplitude),
+ GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusVelocityPeriod));
+ }
+
+ //! \brief Returns the pressure at the inflow.
+ const Scalar refPressure() const
+ {
+ return refPressure_ + variation_(GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusPressureAmplitude),
+ GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusPressurePeriod));
+ }
+
+ //! \brief Returns the mass fraction at the inflow.
+ const Scalar refMassfrac() const
+ {
+ return refMassfrac_ + variation_(GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusConcentrationAmplitude),
+ GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusConcentrationPeriod));
+ }
+
+ //! \brief Returns the temperature at the inflow.
+ const Scalar refTemperature() const
+ {
+ return refTemperature_ + variation_(GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusTemperatureAmplitude),
+ GET_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, SinusTemperaturePeriod));
+ }
+
+private:
+ // Internal method for the initial condition (reused for the dirichlet conditions!)
+ void initial_(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+
+ FluidState fluidState;
+ updateFluidStateForBC_(fluidState);
+ const Scalar density = FluidSystem::density(fluidState, phaseIdx);
+
+ values[velocityXIdx] = xVelocity_(globalPos);
+ values[velocityYIdx] = 0.0;
+ values[pressureIdx] = refPressure()
+ + density * this->gravity()[1] * (globalPos[1] - bBoxMin_[1]);
+ values[massOrMoleFracIdx] = refMassfrac();
+ }
+
+ // Set the profile of the inflow velocity (horizontal direction)
+ const Scalar xVelocity_(const GlobalPosition &globalPos) const
+ {
+ if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
+ return 0.0;
+ return refVelocity();
+ }
+
+ // Updates the fluid state to obtain required quantities for IC/BC
+ void updateFluidStateForBC_(FluidState& fluidState) const
+ {
+ fluidState.setTemperature(refTemperature());
+ fluidState.setPressure(phaseIdx, refPressure());
+
+ Scalar massFraction[numComponents];
+ massFraction[transportCompIdx] = refMassfrac();
+ massFraction[phaseCompIdx] = 1 - massFraction[transportCompIdx];
+
+ // calculate average molar mass of the gas phase
+ Scalar M1 = FluidSystem::molarMass(transportCompIdx);
+ Scalar M2 = FluidSystem::molarMass(phaseCompIdx);
+ Scalar X2 = massFraction[phaseCompIdx];
+ Scalar massToMoleDenominator = M2 + X2*(M1 - M2);
+
+ fluidState.setMoleFraction(phaseIdx, transportCompIdx, massFraction[transportCompIdx]*M2/massToMoleDenominator);
+ fluidState.setMoleFraction(phaseIdx, phaseCompIdx, massFraction[phaseCompIdx]*M1/massToMoleDenominator);
+ }
+
+ // can be used for the variation of a boundary condition
+ const Scalar variation_(const Scalar amplitude, const Scalar period) const
+ { return sin(2*M_PI*this->timeManager().time()/period) * amplitude; }
+
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < bBoxMin_[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > bBoxMax_[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < bBoxMin_[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > bBoxMax_[1] - eps_; }
+
+ bool onBoundary_(const GlobalPosition &globalPos) const
+ {
+ return (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)
+ || onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos));
+ }
+
+ // spatial parameters
+ SpatialParams spatialParams_;
+
+ static constexpr Scalar eps_ = 1e-8;
+ GlobalPosition bBoxMin_;
+ GlobalPosition bBoxMax_;
+ Scalar runUpDistanceX_;
+
+ Scalar refVelocity_;
+ Scalar refPressure_;
+ Scalar refMassfrac_;
+ Scalar refTemperature_;
+ Scalar alphaBJ_;
+};
+} //end namespace Dumux
+
+#endif // DUMUX_ZEROEQ2C_SUBPROBLEM_HH
diff --git a/test/multidomain/CMakeLists.txt b/test/multidomain/CMakeLists.txt
index 3cc98ca7e0978f812ae16ac6c6fd49c0f2aebcb0..f562a68846f0beecb7604f8c34aeca444438a0bc 100644
--- a/test/multidomain/CMakeLists.txt
+++ b/test/multidomain/CMakeLists.txt
@@ -1,2 +1,4 @@
add_subdirectory("2cstokes2p2c")
add_subdirectory("2cnistokes2p2cni")
+add_subdirectory("2czeroeq2p2c")
+add_subdirectory("2cnizeroeq2p2cni")
diff --git a/test/multidomain/Makefile.am b/test/multidomain/Makefile.am
index ca63840bcb1e0757359ef33b91ca7061f2d79482..ec367d39e78df4fdd84c2c43ad227a87ccbadd99 100644
--- a/test/multidomain/Makefile.am
+++ b/test/multidomain/Makefile.am
@@ -1,8 +1,9 @@
SUBDIRS = \
2cstokes2p2c \
- 2cnistokes2p2cni
+ 2cnistokes2p2cni \
+ 2czeroeq2p2c \
+ 2cnizeroeq2p2cni
EXTRA_DIST = CMakeLists.txt installRequiredModules.sh
include $(top_srcdir)/am/global-rules
-
diff --git a/test/references/2cnizeroeq2p2cni-ff-reference.vtu b/test/references/2cnizeroeq2p2cni-ff-reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..182c10dfa1d72b26e8cbaed59ff898a6b8361055
--- /dev/null
+++ b/test/references/2cnizeroeq2p2cni-ff-reference.vtu
@@ -0,0 +1,289 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+ <UnstructuredGrid>
+ <Piece NumberOfCells="150" NumberOfPoints="176">
+ <PointData Scalars="P" Vectors="v">
+ <DataArray type="Float32" Name="P" NumberOfComponents="1" format="ascii">
+ 100015 99988.6 100015 99988.6 100014 100014 99989.6 99989.6 100012 100012 99991.3 99991.2
+ 100010 100010 99993.5 99993 100008 100008 99995 99994.8 100006 100006 99996.5 99996.5
+ 100004 100004 99998.2 99998.2 100002 100002 100000 100000 100018 99987.6 100014 99989.5
+ 100012 99991.3 100010 99993.3 100008 99994.9 100006 99996.5 100004 99998.2 100002 100000
+ 100017 99987.4 100014 99989.2 100012 99991.1 100010 99992.9 100008 99994.7 100006 99996.5
+ 100004 99998.2 100003 100000 100017 99987.1 100015 99989.1 100012 99991 100010 99993
+ 100008 99994.7 100006 99996.5 100005 99998.2 100003 100000 100017 99987.1 100015 99989
+ 100012 99990.9 100010 99992.7 100009 99994.5 100007 99996.4 100005 99998.2 100003 100000
+ 100017 99987.1 100015 99989.1 100012 99991 100010 99992.8 100008 99994.6 100006 99996.4
+ 100005 99998.2 100003 100000 100017 99987.4 100014 99989.2 100012 99991 100010 99992.8
+ 100008 99994.6 100006 99996.4 100004 99998.2 100003 100000 100018 99987.6 100014 99989.4
+ 100012 99991.1 100010 99992.9 100008 99994.7 100006 99996.5 100004 99998.2 100002 100000
+ 100015 99988.6 100014 99989.6 100012 99991.2 100010 99992.9 100008 99994.7 100006 99996.5
+ 100004 99998.2 100002 100000 100015 99988.6 100014 99989.5 100012 99991.2 100010 99992.9
+ 100008 99994.7 100006 99996.5 100004 99998.2 100002 100000
+ </DataArray>
+ <DataArray type="Float32" Name="delP" NumberOfComponents="1" format="ascii">
+ 14.9338 -11.4362 14.9567 -11.4133 13.6823 13.6848 -10.4036 -10.4002 11.9229 11.9373 -8.69525 -8.76971
+ 9.67047 10.1029 -6.48061 -6.96785 7.95471 8.18343 -5.02132 -5.19894 6.31562 6.35073 -3.54476 -3.54993
+ 4.40077 4.40192 -1.76312 -1.76231 2.37932 2.36774 0 0 17.6635 -12.3723 14.1585 -10.544
+ 12.0087 -8.71708 9.84568 -6.65659 8.0536 -5.06219 6.30137 -3.49891 4.39634 -1.76843 2.43834 0
+ 16.9772 -12.621 14.4644 -10.7638 12.2541 -8.92822 10.2815 -7.12957 8.31592 -5.29699 6.38596 -3.5383
+ 4.45793 -1.7844 2.56562 0 16.617 -12.863 14.5095 -10.8791 12.3664 -8.9555 10.2771 -7.04887
+ 8.32222 -5.26678 6.41372 -3.52887 4.52993 -1.78052 2.70329 0 16.5325 -12.9067 14.5336 -10.9932
+ 12.5005 -9.12334 10.4926 -7.29367 8.50707 -5.46085 6.55654 -3.6378 4.63928 -1.81918 2.76135 0
+ 16.6077 -12.8631 14.5293 -10.9078 12.4203 -9.01611 10.3719 -7.15972 8.39268 -5.35021 6.457 -3.56974
+ 4.55989 -1.79681 2.72104 0 16.9641 -12.6195 14.4832 -10.7907 12.3003 -8.96944 10.2949 -7.15503
+ 8.34436 -5.34937 6.41146 -3.56058 4.48226 -1.79292 2.57921 0 17.6386 -12.3734 14.1919 -10.6009
+ 12.1215 -8.85314 10.2014 -7.08082 8.28975 -5.30248 6.36978 -3.52894 4.4347 -1.77498 2.45624 0
+ 14.9327 -11.4144 13.7137 -10.4456 12.0166 -8.83874 10.175 -7.08921 8.28138 -5.3101 6.36411 -3.53001
+ 4.42257 -1.76134 2.38053 0 14.9098 -11.4373 13.7117 -10.4515 12.0161 -8.83968 10.175 -7.08905
+ 8.28152 -5.30993 6.36418 -3.53011 4.42294 -1.7621 2.39274 0
+ </DataArray>
+ <DataArray type="Float32" Name="X^H2O" NumberOfComponents="1" format="ascii">
+ 0.008 0.00622658 0.008 0.00792512 0.00598599 0.00788434 0.00583605 0.00784301 0.00591219 0.00781546 0.00621079 0.00779806
+ 0.00754476 0.00779776 0.00939182 0.00782734 0.0107062 0.00786771 0.0109197 0.00792815 0.0115109 0.00799964 0.0109309 0.0080463
+ 0.0104995 0.00807898 0.0103333 0.00810682 0.00775501 0.00809998 0.0059721 0.00766512 0.008 0.00799624 0.00799233 0.00798685
+ 0.00798182 0.00797669 0.00797301 0.00796895 0.00796663 0.00796574 0.00796705 0.00796983 0.00797317 0.00797692 0.00797625 0.00771151
+ 0.008 0.00799992 0.00799975 0.00799942 0.00799902 0.00799848 0.00799796 0.00799728 0.00799661 0.00799587 0.0079952 0.00799462
+ 0.00799418 0.00799353 0.00798616 0.00775718 0.008 0.008 0.008 0.00799999 0.00799999 0.00799997 0.00799995 0.00799991
+ 0.00799987 0.00799982 0.00799976 0.00799968 0.00799957 0.0079988 0.00798684 0.00775569 0.008 0.008 0.008 0.008
+ 0.00799999 0.00799999 0.00799998 0.00799998 0.00799997 0.00799996 0.00799995 0.00799994 0.00799988 0.00799908 0.00798441 0.00775215
+ 0.008 0.00799999 0.00799996 0.00799993 0.00799988 0.00799983 0.00799976 0.00799967 0.00799958 0.00799948 0.00799937 0.00799925
+ 0.00799913 0.00799881 0.0079906 0.00776135 0.008 0.00799971 0.0079993 0.0079987 0.00799806 0.00799731 0.00799654 0.00799569
+ 0.00799484 0.00799394 0.00799308 0.00799218 0.00799137 0.00799059 0.00798595 0.00775814 0.008 0.00799441 0.00798944 0.00798328
+ 0.00797817 0.00797299 0.00796856 0.00796425 0.00796052 0.00795694 0.00795389 0.00795102 0.00794875 0.00794656 0.00794231 0.0076764
+ 0.008 0.00791749 0.0078772 0.00783769 0.0078172 0.0078026 0.00779499 0.00779114 0.00779009 0.00779107 0.00779335 0.00779664
+ 0.00780174 0.00780312 0.00779054 0.00737342 0.008 0.00624493 0.00601681 0.00585947 0.00602312 0.006231 0.00647328 0.00672531
+ 0.00693001 0.00711547 0.00724865 0.00735113 0.00745775 0.00772897 0.006229 0.00494063
+ </DataArray>
+ <DataArray type="Float32" Name="rho" NumberOfComponents="1" format="ascii">
+ 1.16276 1.1637 1.16276 1.16251 1.16416 1.16283 1.16399 1.16257 1.16419 1.16285 1.16374 1.16262
+ 1.17425 1.16286 1.19112 1.16282 1.20091 1.16305 1.19934 1.16311 1.1602 1.16321 1.16049 1.16306
+ 1.16089 1.16313 1.16093 1.16304 1.16279 1.16307 1.18744 1.16835 1.16279 1.16245 1.16276 1.16247
+ 1.16274 1.1625 1.16272 1.16253 1.16271 1.16256 1.16269 1.16258 1.16267 1.1626 1.16265 1.166
+ 1.16278 1.16244 1.16276 1.16246 1.16273 1.16249 1.16271 1.16251 1.16269 1.16253 1.16267 1.16255
+ 1.16264 1.16257 1.16263 1.16549 1.16278 1.16244 1.16276 1.16246 1.16273 1.16248 1.16271 1.16251
+ 1.16268 1.16253 1.16266 1.16255 1.16264 1.16257 1.16263 1.16552 1.16278 1.16244 1.16276 1.16246
+ 1.16273 1.16248 1.16271 1.1625 1.16269 1.16252 1.16266 1.16255 1.16264 1.16257 1.16263 1.16556
+ 1.16278 1.16244 1.16276 1.16246 1.16273 1.16248 1.16271 1.1625 1.16269 1.16253 1.16266 1.16255
+ 1.16264 1.16257 1.16263 1.16553 1.16278 1.16244 1.16276 1.16246 1.16273 1.16249 1.16271 1.16251
+ 1.16269 1.16253 1.16267 1.16255 1.16265 1.16257 1.16263 1.16551 1.16279 1.16245 1.16276 1.16248
+ 1.16274 1.1625 1.16273 1.16253 1.16271 1.16255 1.16269 1.16258 1.16267 1.1626 1.16266 1.16602
+ 1.16276 1.16251 1.16283 1.16258 1.16285 1.16262 1.16284 1.16265 1.16283 1.16267 1.1628 1.16268
+ 1.16277 1.1627 1.16277 1.1681 1.16276 1.16369 1.16414 1.16397 1.16412 1.16373 1.16378 1.1634
+ 1.16343 1.16315 1.16319 1.163 1.16302 1.16276 1.16386 1.18809
+ </DataArray>
+ <DataArray type="Float32" Name="mu" NumberOfComponents="1" format="ascii">
+ 1.81222e-05 1.81479e-05 1.81222e-05 1.81233e-05 1.81514e-05 1.81239e-05 1.81536e-05 1.81245e-05 1.81525e-05 1.8125e-05 1.81482e-05 1.81252e-05
+ 1.79936e-05 1.8125e-05 1.77501e-05 1.81224e-05 1.76091e-05 1.81208e-05 1.76206e-05 1.81168e-05 1.80716e-05 1.81156e-05 1.80799e-05 1.81149e-05
+ 1.80861e-05 1.81144e-05 1.80885e-05 1.81141e-05 1.81258e-05 1.81144e-05 1.78719e-05 1.80603e-05 1.81222e-05 1.81223e-05 1.81223e-05 1.81224e-05
+ 1.81225e-05 1.81226e-05 1.81226e-05 1.81227e-05 1.81227e-05 1.81227e-05 1.81226e-05 1.81225e-05 1.81224e-05 1.81224e-05 1.81224e-05 1.80876e-05
+ 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81223e-05 1.81223e-05 1.81223e-05 1.81223e-05 1.81223e-05 1.81223e-05 1.81223e-05
+ 1.81223e-05 1.81223e-05 1.81224e-05 1.80927e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05
+ 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81224e-05 1.80924e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05
+ 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81224e-05 1.8092e-05
+ 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05
+ 1.81222e-05 1.81222e-05 1.81223e-05 1.80921e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81222e-05 1.81223e-05 1.81223e-05 1.81223e-05 1.81223e-05
+ 1.81223e-05 1.81223e-05 1.81223e-05 1.81223e-05 1.81224e-05 1.81224e-05 1.81224e-05 1.80924e-05 1.81222e-05 1.81223e-05 1.81224e-05 1.81225e-05
+ 1.81226e-05 1.81226e-05 1.81227e-05 1.81228e-05 1.81228e-05 1.81229e-05 1.81229e-05 1.8123e-05 1.8123e-05 1.8123e-05 1.81231e-05 1.80881e-05
+ 1.81222e-05 1.81234e-05 1.8124e-05 1.81246e-05 1.81249e-05 1.81251e-05 1.81253e-05 1.81253e-05 1.81253e-05 1.81253e-05 1.81253e-05 1.81252e-05
+ 1.81252e-05 1.81251e-05 1.81252e-05 1.807e-05 1.81222e-05 1.81477e-05 1.8151e-05 1.81533e-05 1.81509e-05 1.81479e-05 1.81444e-05 1.81407e-05
+ 1.81377e-05 1.8135e-05 1.81331e-05 1.81316e-05 1.81301e-05 1.81262e-05 1.81479e-05 1.78877e-05
+ </DataArray>
+ <DataArray type="Float32" Name="v" NumberOfComponents="3" format="ascii">
+ 0 0 0 0 0 0 3.5 0 0 3.42882 0.102831 0
+ 0 0 0 3.12033 0.00888965 0 0 0 0 3.23492 0.0250088 0
+ 0 0 0 2.91421 0.00796738 0 0 0 0 3.11557 0.00537881 0
+ 0.0265927 -0.00442212 0 2.79992 -0.0368189 0 0.260362 0.00671295 0 3.025 0.0125029 0
+ 0.0135791 -0.00557923 0 2.73487 -0.00563452 0 0.25174 0.00373404 0 2.96154 0.00574155 0
+ 0 0 0 2.67089 0.0155556 0 0 0 0 2.92057 0.00184655 0
+ 0 0 0 2.6197 0.00197139 0 0 0 0 2.8943 0.00118524 0
+ 0 0 0 2.62412 -0.0366102 0 1.00648 -0.00725145 0 2.87305 -0.047548 0
+ 3.5 0 0 3.83764 0.1139 0 3.60914 0.0365758 0 3.81297 0.0620856 0
+ 3.56647 0.0264748 0 3.76552 0.035235 0 3.49958 0.0121173 0 3.67212 0.0273948 0
+ 3.44497 0.0127998 0 3.61618 0.0231502 0 3.39317 0.0201092 0 3.58326 0.0158747 0
+ 3.34994 0.00755705 0 3.55237 0.00380736 0 3.3234 -0.0268679 0 3.50477 -0.0323767 0
+ 3.5 0 0 3.7965 0.0530881 0 3.58192 0.0348213 0 3.83174 0.0529557 0
+ 3.62172 0.030784 0 3.84598 0.0396834 0 3.63792 0.0129947 0 3.85077 0.036085 0
+ 3.64583 0.0221042 0 3.8514 0.0338985 0 3.6521 0.0276761 0 3.85637 0.0266335 0
+ 3.65489 0.0158968 0 3.85515 0.0109235 0 3.64941 -0.011512 0 3.83131 -0.015253 0
+ 3.5 0 0 3.74876 0.0173781 0 3.52306 0.0184163 0 3.76854 0.0272178 0
+ 3.55177 0.0188408 0 3.78773 0.0246647 0 3.57395 0.0173603 0 3.79758 0.0221225 0
+ 3.59405 0.0156268 0 3.81228 0.0206952 0 3.61376 0.0146725 0 3.82848 0.01665 0
+ 3.6294 0.00902373 0 3.83785 0.00764436 0 3.63502 -0.00146545 0 3.83273 -0.00293924 0
+ 3.5 0 0 3.73744 -6.56461e-05 0 3.51322 -0.000644705 0 3.75131 -1.30716e-05 0
+ 3.5357 -0.00150479 0 3.76758 -0.000135594 0 3.55741 -0.00264728 0 3.78369 0.000528393 0
+ 3.57659 -0.000911225 0 3.79717 0.00153294 0 3.59346 0.00068055 0 3.80858 0.00175858 0
+ 3.60655 0.000986566 0 3.81564 0.00129446 0 3.6126 0.000723964 0 3.81597 0.000595173 0
+ 3.5 0 0 3.74616 -0.0175151 0 3.52292 -0.0192653 0 3.76553 -0.0273057 0
+ 3.55131 -0.0206185 0 3.7844 -0.0249855 0 3.57481 -0.017902 0 3.80017 -0.0217802 0
+ 3.59481 -0.0160212 0 3.81513 -0.0189575 0 3.61262 -0.0135217 0 3.82864 -0.0148049 0
+ 3.62705 -0.00807985 0 3.83695 -0.0065884 0 3.63233 0.00194918 0 3.83119 0.00330661 0
+ 3.5 0 0 3.79274 -0.0530063 0 3.58152 -0.0354653 0 3.82772 -0.0522742 0
+ 3.62022 -0.0336115 0 3.84364 -0.0406984 0 3.63772 -0.0298172 0 3.85132 -0.032684 0
+ 3.64621 -0.0259332 0 3.85445 -0.0273553 0 3.65001 -0.0213722 0 3.85506 -0.0218203 0
+ 3.65176 -0.0122009 0 3.85357 -0.00826005 0 3.64669 0.0139402 0 3.82972 0.0175634 0
+ 3.5 0 0 3.83175 -0.113618 0 3.60887 -0.0371117 0 3.80533 -0.0610134 0
+ 3.56638 -0.0285282 0 3.74568 -0.0359244 0 3.50107 -0.0209872 0 3.68242 -0.0246889 0
+ 3.43861 -0.0165681 0 3.62575 -0.0187842 0 3.3855 -0.0126929 0 3.57974 -0.0142077 0
+ 3.3455 -0.00622179 0 3.55098 -0.00351178 0 3.32176 0.0274053 0 3.50484 0.0330586 0
+ 3.5 0 0 3.42268 -0.102815 0 3.12092 -0.00817394 0 3.22822 -0.025862 0
+ 2.91668 -0.00359975 0 3.10551 -0.00508346 0 2.80001 -0.000452951 0 3.02456 5.83809e-05 0
+ 2.722 0.000217902 0 2.96428 0.000234495 0 2.66352 -1.98939e-05 0 2.91884 -0.000836122 0
+ 2.62082 0.00107814 0 2.89698 -9.37616e-05 0 2.62635 0.0384585 0 2.87645 0.0493796 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 1.01298 0.00752403 0
+ </DataArray>
+ <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
+ 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.151 298.15 298.151 298.15 298.151
+ 295.299 298.147 290.745 298.1 288.19 298.079 288.493 298.014 298.15 298.009 298.15 298.009
+ 298.15 298.009 298.15 298.01 298.15 298.015 292.268 296.739 298.15 298.15 298.15 298.15
+ 298.15 298.15 298.15 298.15 298.15 298.149 298.147 298.147 298.146 298.146 298.145 297.33
+ 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15
+ 298.15 298.15 298.149 297.451 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15
+ 298.15 298.15 298.15 298.15 298.15 298.15 298.149 297.444 298.15 298.15 298.15 298.15
+ 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.149 297.435
+ 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15
+ 298.15 298.15 298.149 297.44 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15
+ 298.15 298.15 298.15 298.15 298.15 298.15 298.15 297.447 298.15 298.15 298.15 298.15
+ 298.15 298.15 298.15 298.15 298.15 298.15 298.151 298.151 298.151 298.151 298.15 297.331
+ 298.15 298.15 298.151 298.151 298.151 298.151 298.151 298.151 298.151 298.151 298.151 298.151
+ 298.151 298.151 298.149 296.856 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15
+ 298.15 298.15 298.15 298.15 298.15 298.15 298.15 292.289
+ </DataArray>
+ </PointData>
+ <Points>
+ <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+ 0 0.25 0 0.05 0.25 0 0 0.262165 0 0.05 0.262165 0
+ 0.1 0.25 0 0.1 0.262165 0 0.15 0.25 0 0.15 0.262165 0
+ 0.2 0.25 0 0.2 0.262165 0 0.25 0.25 0 0.25 0.262165 0
+ 0.3 0.25 0 0.3 0.262165 0 0.35 0.25 0 0.35 0.262165 0
+ 0.4 0.25 0 0.4 0.262165 0 0.45 0.25 0 0.45 0.262165 0
+ 0.5 0.25 0 0.5 0.262165 0 0.55 0.25 0 0.55 0.262165 0
+ 0.6 0.25 0 0.6 0.262165 0 0.65 0.25 0 0.65 0.262165 0
+ 0.7 0.25 0 0.7 0.262165 0 0.75 0.25 0 0.75 0.262165 0
+ 0 0.283454 0 0.05 0.283454 0 0.1 0.283454 0 0.15 0.283454 0
+ 0.2 0.283454 0 0.25 0.283454 0 0.3 0.283454 0 0.35 0.283454 0
+ 0.4 0.283454 0 0.45 0.283454 0 0.5 0.283454 0 0.55 0.283454 0
+ 0.6 0.283454 0 0.65 0.283454 0 0.7 0.283454 0 0.75 0.283454 0
+ 0 0.320709 0 0.05 0.320709 0 0.1 0.320709 0 0.15 0.320709 0
+ 0.2 0.320709 0 0.25 0.320709 0 0.3 0.320709 0 0.35 0.320709 0
+ 0.4 0.320709 0 0.45 0.320709 0 0.5 0.320709 0 0.55 0.320709 0
+ 0.6 0.320709 0 0.65 0.320709 0 0.7 0.320709 0 0.75 0.320709 0
+ 0 0.385906 0 0.05 0.385906 0 0.1 0.385906 0 0.15 0.385906 0
+ 0.2 0.385906 0 0.25 0.385906 0 0.3 0.385906 0 0.35 0.385906 0
+ 0.4 0.385906 0 0.45 0.385906 0 0.5 0.385906 0 0.55 0.385906 0
+ 0.6 0.385906 0 0.65 0.385906 0 0.7 0.385906 0 0.75 0.385906 0
+ 0 0.5 0 0.05 0.5 0 0.1 0.5 0 0.15 0.5 0
+ 0.2 0.5 0 0.25 0.5 0 0.3 0.5 0 0.35 0.5 0
+ 0.4 0.5 0 0.45 0.5 0 0.5 0.5 0 0.55 0.5 0
+ 0.6 0.5 0 0.65 0.5 0 0.7 0.5 0 0.75 0.5 0
+ 0 0.614094 0 0.05 0.614094 0 0.1 0.614094 0 0.15 0.614094 0
+ 0.2 0.614094 0 0.25 0.614094 0 0.3 0.614094 0 0.35 0.614094 0
+ 0.4 0.614094 0 0.45 0.614094 0 0.5 0.614094 0 0.55 0.614094 0
+ 0.6 0.614094 0 0.65 0.614094 0 0.7 0.614094 0 0.75 0.614094 0
+ 0 0.679291 0 0.05 0.679291 0 0.1 0.679291 0 0.15 0.679291 0
+ 0.2 0.679291 0 0.25 0.679291 0 0.3 0.679291 0 0.35 0.679291 0
+ 0.4 0.679291 0 0.45 0.679291 0 0.5 0.679291 0 0.55 0.679291 0
+ 0.6 0.679291 0 0.65 0.679291 0 0.7 0.679291 0 0.75 0.679291 0
+ 0 0.716546 0 0.05 0.716546 0 0.1 0.716546 0 0.15 0.716546 0
+ 0.2 0.716546 0 0.25 0.716546 0 0.3 0.716546 0 0.35 0.716546 0
+ 0.4 0.716546 0 0.45 0.716546 0 0.5 0.716546 0 0.55 0.716546 0
+ 0.6 0.716546 0 0.65 0.716546 0 0.7 0.716546 0 0.75 0.716546 0
+ 0 0.737835 0 0.05 0.737835 0 0.1 0.737835 0 0.15 0.737835 0
+ 0.2 0.737835 0 0.25 0.737835 0 0.3 0.737835 0 0.35 0.737835 0
+ 0.4 0.737835 0 0.45 0.737835 0 0.5 0.737835 0 0.55 0.737835 0
+ 0.6 0.737835 0 0.65 0.737835 0 0.7 0.737835 0 0.75 0.737835 0
+ 0 0.75 0 0.05 0.75 0 0.1 0.75 0 0.15 0.75 0
+ 0.2 0.75 0 0.25 0.75 0 0.3 0.75 0 0.35 0.75 0
+ 0.4 0.75 0 0.45 0.75 0 0.5 0.75 0 0.55 0.75 0
+ 0.6 0.75 0 0.65 0.75 0 0.7 0.75 0 0.75 0.75 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 20 22 23 21 22 24 25 23
+ 24 26 27 25 26 28 29 27 28 30 31 29
+ 2 3 33 32 3 5 34 33 5 7 35 34
+ 7 9 36 35 9 11 37 36 11 13 38 37
+ 13 15 39 38 15 17 40 39 17 19 41 40
+ 19 21 42 41 21 23 43 42 23 25 44 43
+ 25 27 45 44 27 29 46 45 29 31 47 46
+ 32 33 49 48 33 34 50 49 34 35 51 50
+ 35 36 52 51 36 37 53 52 37 38 54 53
+ 38 39 55 54 39 40 56 55 40 41 57 56
+ 41 42 58 57 42 43 59 58 43 44 60 59
+ 44 45 61 60 45 46 62 61 46 47 63 62
+ 48 49 65 64 49 50 66 65 50 51 67 66
+ 51 52 68 67 52 53 69 68 53 54 70 69
+ 54 55 71 70 55 56 72 71 56 57 73 72
+ 57 58 74 73 58 59 75 74 59 60 76 75
+ 60 61 77 76 61 62 78 77 62 63 79 78
+ 64 65 81 80 65 66 82 81 66 67 83 82
+ 67 68 84 83 68 69 85 84 69 70 86 85
+ 70 71 87 86 71 72 88 87 72 73 89 88
+ 73 74 90 89 74 75 91 90 75 76 92 91
+ 76 77 93 92 77 78 94 93 78 79 95 94
+ 80 81 97 96 81 82 98 97 82 83 99 98
+ 83 84 100 99 84 85 101 100 85 86 102 101
+ 86 87 103 102 87 88 104 103 88 89 105 104
+ 89 90 106 105 90 91 107 106 91 92 108 107
+ 92 93 109 108 93 94 110 109 94 95 111 110
+ 96 97 113 112 97 98 114 113 98 99 115 114
+ 99 100 116 115 100 101 117 116 101 102 118 117
+ 102 103 119 118 103 104 120 119 104 105 121 120
+ 105 106 122 121 106 107 123 122 107 108 124 123
+ 108 109 125 124 109 110 126 125 110 111 127 126
+ 112 113 129 128 113 114 130 129 114 115 131 130
+ 115 116 132 131 116 117 133 132 117 118 134 133
+ 118 119 135 134 119 120 136 135 120 121 137 136
+ 121 122 138 137 122 123 139 138 123 124 140 139
+ 124 125 141 140 125 126 142 141 126 127 143 142
+ 128 129 145 144 129 130 146 145 130 131 147 146
+ 131 132 148 147 132 133 149 148 133 134 150 149
+ 134 135 151 150 135 136 152 151 136 137 153 152
+ 137 138 154 153 138 139 155 154 139 140 156 155
+ 140 141 157 156 141 142 158 157 142 143 159 158
+ 144 145 161 160 145 146 162 161 146 147 163 162
+ 147 148 164 163 148 149 165 164 149 150 166 165
+ 150 151 167 166 151 152 168 167 152 153 169 168
+ 153 154 170 169 154 155 171 170 155 156 172 171
+ 156 157 173 172 157 158 174 173 158 159 175 174
+ </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
+ </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
+ </DataArray>
+ </Cells>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>
diff --git a/test/references/2cnizeroeq2p2cni-pm-reference.vtu b/test/references/2cnizeroeq2p2cni-pm-reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..f44c037b81f9177ed9e6d11d026ff7c3294c0588
--- /dev/null
+++ b/test/references/2cnizeroeq2p2cni-pm-reference.vtu
@@ -0,0 +1,174 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+ <UnstructuredGrid>
+ <Piece NumberOfCells="25" NumberOfPoints="36">
+ <PointData Scalars="Sn" Vectors="velocityW">
+ <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
+ 0.0750879 0.0750906 0.875107 0.874938 0.0750629 0.8746 0.0750438 0.874156 0.0750472 0.873755 0.0750519 0.873594
+ 0.938627 0.9388 0.937982 0.937143 0.936472 0.93608 0.957838 0.96014 0.958001 0.955501 0.954155 0.953244
+ 0.97054 0.979307 0.976986 0.971314 0.966526 0.963334 0.965619 1 1 0.984275 0.976549 0.96919
+ </DataArray>
+ <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
+ 0.924912 0.924909 0.124893 0.125062 0.924937 0.1254 0.924956 0.125844 0.924953 0.126245 0.924948 0.126406
+ 0.0613725 0.0611999 0.0620181 0.0628566 0.0635281 0.0639198 0.0421615 0.0398596 0.0419992 0.044499 0.0458449 0.046756
+ 0.0294595 0.0206932 0.0230136 0.0286856 0.0334739 0.036666 0.0343812 0 0 0.0157249 0.0234506 0.0308099
+ </DataArray>
+ <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
+ 100007 100007 100004 100004 100006 100004 100006 100003 100006 100002 100006 100002
+ 100004 100004 100003 100002 100001 100001 100005 100005 100001 100002 99999.6 99999.2
+ 100005 100007 99997.8 100005 99997.5 99998.7 100005 100010 99993.7 100008 99995 99999.1
+ </DataArray>
+ <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
+ 98859.7 98859.7 97908.5 97908.8 98859.5 97909.2 98859.3 97909.8 98859.1 97910.5 98859 97910.7
+ 97653.9 97652.3 97656.1 97660.6 97663.6 97665.5 97505.9 97481.9 97500.1 97525.3 97534.6 97542.2
+ 97349.1 97174.5 97221.3 97335.9 97399.2 97440.2 97418.1 96325.4 96309.5 97016.2 97228.2 97363.4
+ </DataArray>
+ <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
+ 1146.91 1146.92 2095.7 2095.24 1146.85 2094.32 1146.8 2093.11 1146.81 2092.02 1146.82 2091.58
+ 2350.55 2351.62 2346.57 2341.49 2337.47 2335.15 2499.36 2522.84 2500.96 2477.14 2464.99 2457.03
+ 2656.29 2832.23 2776.46 2668.7 2598.33 2558.45 2586.5 3684.13 3684.13 2991.69 2766.87 2635.66
+ </DataArray>
+ <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
+ 997.055 997.055 997.114 997.117 997.055 997.121 997.055 997.12 997.055 997.112 997.055 997.108
+ 997.323 997.376 997.435 997.452 997.382 997.319 997.534 997.788 997.973 998.118 997.924 997.598
+ 997.554 998.08 998.351 998.719 998.52 997.669 997.668 997.756 998.682 999.103 999.057 997.504
+ </DataArray>
+ <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
+ 1.15435 1.15435 1.15541 1.15546 1.15435 1.15554 1.15435 1.15551 1.15434 1.15537 1.15434 1.15529
+ 1.15935 1.16036 1.16148 1.16179 1.16042 1.15922 1.16345 1.16855 1.17241 1.17557 1.17134 1.16465
+ 1.16383 1.1748 1.18081 1.18999 1.18487 1.16605 1.16611 1.17425 1.19112 1.20091 1.19934 1.16279
+ </DataArray>
+ <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
+ 907.079 907.066 2.45157 2.46061 907.198 2.47902 907.289 2.50503 907.273 2.53001 907.251 2.54044
+ 0.291851 0.287916 0.298081 0.309989 0.3223 0.330315 0.0894122 0.0729038 0.0843765 0.0996158 0.11178 0.122983
+ 0.0278011 0.00763081 0.0108525 0.0221623 0.0380559 0.0564579 0.0458204 0 0 0.00229257 0.0104817 0.0324517
+ </DataArray>
+ <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
+ 242.962 242.983 45242 45227.3 242.768 45197.6 242.621 45156.2 242.647 45117.1 242.684 45101
+ 51931.4 51970.2 51900.7 51812.5 51710.4 51643.1 53991.8 54304.6 54197.3 54033.4 53812.4 53583.1
+ 55016.8 55761.1 55786.9 55666.3 55200.9 54524.3 54706.7 55575.2 56337.7 56602.8 56223.1 54903.3
+ </DataArray>
+ <DataArray type="Float32" Name="X_liquid^H2O" NumberOfComponents="1" format="ascii">
+ 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978
+ 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999977 0.999977 0.999977 0.999977 0.999978
+ 0.999978 0.999977 0.999976 0.999975 0.999976 0.999977 0.999977 0.452322 0.745811 0.999974 0.999974 0.999978
+ </DataArray>
+ <DataArray type="Float32" Name="X_liquid^Air" NumberOfComponents="1" format="ascii">
+ 2.15483e-05 2.15483e-05 2.16354e-05 2.16394e-05 2.15483e-05 2.16463e-05 2.15482e-05 2.1644e-05 2.15482e-05 2.1633e-05 2.15482e-05 2.16266e-05
+ 2.19593e-05 2.20444e-05 2.21408e-05 2.21677e-05 2.20522e-05 2.19511e-05 2.23077e-05 2.27596e-05 2.31174e-05 2.3416e-05 2.30191e-05 2.24169e-05
+ 2.23408e-05 2.3338e-05 2.39345e-05 2.48871e-05 2.43496e-05 2.25407e-05 2.25412e-05 2.30241e-05 2.49001e-05 2.61252e-05 2.59543e-05 2.22553e-05
+ </DataArray>
+ <DataArray type="Float32" Name="X_gas^H2O" NumberOfComponents="1" format="ascii">
+ 0.0199555 0.0199556 0.0196794 0.0196667 0.0199556 0.0196451 0.0199556 0.0196518 0.0199557 0.019686 0.0199557 0.0197058
+ 0.0187047 0.0184602 0.018189 0.0181141 0.0184367 0.0187266 0.0177352 0.0165855 0.015752 0.0151054 0.0159738 0.0174441
+ 0.017647 0.015272 0.0140715 0.0124361 0.0133201 0.0171252 0.0171265 0.00754476 0.00939182 0.0107062 0.0109197 0.0178733
+ </DataArray>
+ <DataArray type="Float32" Name="X_gas^Air" NumberOfComponents="1" format="ascii">
+ 0.980044 0.980044 0.980321 0.980333 0.980044 0.980355 0.980044 0.980348 0.980044 0.980314 0.980044 0.980294
+ 0.981295 0.98154 0.981811 0.981886 0.981563 0.981273 0.982265 0.983414 0.984248 0.984895 0.984026 0.982556
+ 0.982353 0.984728 0.985929 0.987564 0.98668 0.982875 0.982873 0.992455 0.990608 0.989294 0.98908 0.982127
+ </DataArray>
+ <DataArray type="Float32" Name="x_liquid^H2O" NumberOfComponents="1" format="ascii">
+ 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987
+ 0.999986 0.999986 0.999986 0.999986 0.999986 0.999986 0.999986 0.999986 0.999986 0.999985 0.999986 0.999986
+ 0.999986 0.999985 0.999985 0.999985 0.999985 0.999986 0.999986 0.45233 0.745821 0.999984 0.999984 0.999986
+ </DataArray>
+ <DataArray type="Float32" Name="x_liquid^Air" NumberOfComponents="1" format="ascii">
+ 1.34047e-05 1.34047e-05 1.34589e-05 1.34614e-05 1.34047e-05 1.34656e-05 1.34046e-05 1.34642e-05 1.34046e-05 1.34574e-05 1.34046e-05 1.34534e-05
+ 1.36603e-05 1.37133e-05 1.37733e-05 1.379e-05 1.37182e-05 1.36553e-05 1.38771e-05 1.41582e-05 1.43808e-05 1.45665e-05 1.43196e-05 1.3945e-05
+ 1.38977e-05 1.4518e-05 1.48891e-05 1.54817e-05 1.51473e-05 1.40221e-05 1.40224e-05 1.43229e-05 1.54899e-05 1.62519e-05 1.61456e-05 1.38445e-05
+ </DataArray>
+ <DataArray type="Float32" Name="x_gas^H2O" NumberOfComponents="1" format="ascii">
+ 0.0316949 0.031695 0.0312615 0.0312416 0.031695 0.0312076 0.0316951 0.0312182 0.0316952 0.0312718 0.0316952 0.031303
+ 0.0297306 0.0293463 0.0289198 0.028802 0.0293094 0.0297651 0.028206 0.0263958 0.0250819 0.0240615 0.0254316 0.0277479
+ 0.0280673 0.0243246 0.0224286 0.0198415 0.0212407 0.0272459 0.0272479 0.0120731 0.015012 0.0170994 0.0174381 0.0284233
+ </DataArray>
+ <DataArray type="Float32" Name="x_gas^Air" NumberOfComponents="1" format="ascii">
+ 0.968305 0.968305 0.968739 0.968758 0.968305 0.968792 0.968305 0.968782 0.968305 0.968728 0.968305 0.968697
+ 0.970269 0.970654 0.97108 0.971198 0.970691 0.970235 0.971794 0.973604 0.974918 0.975938 0.974568 0.972252
+ 0.971933 0.975675 0.977571 0.980159 0.978759 0.972754 0.972752 0.987927 0.984988 0.982901 0.982562 0.971577
+ </DataArray>
+ <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
+ 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41
+ 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41
+ 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41
+ </DataArray>
+ <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
+ 298.15 298.15 297.919 297.908 298.15 297.89 298.15 297.895 298.15 297.924 298.15 297.941
+ 297.081 296.864 296.621 296.553 296.843 297.099 296.208 295.116 294.281 293.607 294.507 295.937
+ 296.126 293.784 292.473 290.52 291.602 295.636 295.638 295.299 290.745 288.19 288.493 296.334
+ </DataArray>
+ <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
+ 3 3 3 3 3 3 3 3 3 3 3 3
+ 3 3 3 3 3 3 3 3 3 3 3 3
+ 3 3 3 3 3 3 3 2 2 3 3 3
+ </DataArray>
+ <DataArray type="Float32" Name="velocityW" NumberOfComponents="3" format="ascii">
+ -1.13827e-08 -9.38867e-07 0 -2.02936e-08 -9.43886e-07 0 3.67637e-10 -6.9645e-07 0 -4.89952e-09 -6.95475e-07 0
+ -3.32727e-08 -9.54168e-07 0 -1.1732e-08 -7.10208e-07 0 -3.44622e-08 -9.68821e-07 0 -1.22953e-08 -7.29101e-07 0
+ -2.19662e-08 -9.83091e-07 0 -8.49334e-09 -7.47566e-07 0 -1.23489e-08 -9.89044e-07 0 -5.69332e-09 -7.57778e-07 0
+ 6.46854e-09 -2.94085e-07 0 -8.5928e-10 -2.75075e-07 0 -9.47256e-09 -2.95962e-07 0 -8.4042e-09 -3.25321e-07 0
+ -5.12068e-09 -3.49662e-07 0 -4.19098e-09 -3.6815e-07 0 1.89538e-08 -6.94712e-08 0 5.01846e-09 -1.6798e-08 0
+ -1.33765e-08 -2.99025e-08 0 -1.34456e-08 -8.18932e-08 0 -8.12791e-09 -1.10338e-07 0 -7.20068e-09 -1.39335e-07 0
+ 3.98571e-08 1.01552e-08 0 1.89117e-08 8.65335e-08 0 -1.14392e-08 1.16956e-07 0 -1.79959e-08 6.05178e-08 0
+ -1.48476e-08 7.61521e-09 0 -1.4548e-08 -4.10402e-08 0 2.09619e-07 2.51157e-08 0 1.04809e-07 1.03575e-07 0
+ -3.39414e-09 1.68041e-07 0 -8.24859e-09 1.23561e-07 0 -1.44586e-08 4.8624e-08 0 -1.92083e-08 -1.31981e-08 0
+ </DataArray>
+ <DataArray type="Float32" Name="velocityN" NumberOfComponents="3" format="ascii">
+ 1.39233e-07 6.47251e-07 0 2.4638e-07 7.19854e-07 0 5.75125e-05 -9.52283e-05 0 9.57313e-05 -6.36067e-05 0
+ 3.88318e-07 8.84205e-07 0 0.00013264 -1.45967e-06 0 3.78512e-07 1.06416e-06 0 0.000122537 1.91172e-05 0
+ 2.30399e-07 1.19748e-06 0 8.0125e-05 5.83971e-05 0 1.26884e-07 1.25052e-06 0 4.65057e-05 7.91642e-05 0
+ 0.000138223 -0.000329088 0 0.000268815 -0.000246146 0 0.000254242 0.000119787 0 0.000208719 -4.34298e-05 0
+ 0.000208455 0.000207451 0 0.000108549 0.000265968 0 7.57763e-05 -0.000522431 0 0.00062318 -0.000653991 0
+ 0.000326743 0.000691165 0 0.000188181 -0.000491347 0 0.000475768 0.000521091 0 5.80777e-05 0.000417074 0
+ -0.000471595 -0.000435573 0 0.00112488 -0.00143926 0 0.000312129 0.00189921 0 2.83532e-05 -0.00151907 0
+ 0.000886834 0.00113801 0 -0.000380138 0.000346331 0 -0.00119016 -0.000293355 0 0.00149089 -0.0019349 0
+ 0.000224842 0.00265946 0 -0.000144805 -0.00217949 0 0.0012347 0.00153314 0 -0.000963249 0.000233374 0
+ </DataArray>
+ </PointData>
+ <CellData Scalars="process rank">
+ <DataArray type="Float32" Name="process rank" 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
+ </DataArray>
+ </CellData>
+ <Points>
+ <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+ 0.25 0 0 0.3 0 0 0.25 0.114094 0 0.3 0.114094 0
+ 0.35 0 0 0.35 0.114094 0 0.4 0 0 0.4 0.114094 0
+ 0.45 0 0 0.45 0.114094 0 0.5 0 0 0.5 0.114094 0
+ 0.25 0.179291 0 0.3 0.179291 0 0.35 0.179291 0 0.4 0.179291 0
+ 0.45 0.179291 0 0.5 0.179291 0 0.25 0.216546 0 0.3 0.216546 0
+ 0.35 0.216546 0 0.4 0.216546 0 0.45 0.216546 0 0.5 0.216546 0
+ 0.25 0.237835 0 0.3 0.237835 0 0.35 0.237835 0 0.4 0.237835 0
+ 0.45 0.237835 0 0.5 0.237835 0 0.25 0.25 0 0.3 0.25 0
+ 0.35 0.25 0 0.4 0.25 0 0.45 0.25 0 0.5 0.25 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 2 3 13 12
+ 3 5 14 13 5 7 15 14 7 9 16 15
+ 9 11 17 16 12 13 19 18 13 14 20 19
+ 14 15 21 20 15 16 22 21 16 17 23 22
+ 18 19 25 24 19 20 26 25 20 21 27 26
+ 21 22 28 27 22 23 29 28 24 25 31 30
+ 25 26 32 31 26 27 33 32 27 28 34 33
+ 28 29 35 34
+ </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
+ </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
+ </DataArray>
+ </Cells>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>
diff --git a/test/references/2czeroeq2p2c-ff-reference.vtu b/test/references/2czeroeq2p2c-ff-reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..eff1e3a54112e29b1c66108dc71d98fb5ff11d2b
--- /dev/null
+++ b/test/references/2czeroeq2p2c-ff-reference.vtu
@@ -0,0 +1,258 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+ <UnstructuredGrid>
+ <Piece NumberOfCells="140" NumberOfPoints="165">
+ <PointData Scalars="P" Vectors="v">
+ <DataArray type="Float32" Name="P" NumberOfComponents="1" format="ascii">
+ 100013 99987.5 100013 99987.5 100012 100012 99989.5 99989.5 100009 100009 99991.8 99991.4
+ 100004 100006 99995.2 99993.9 100002 100003 99997.4 99996.6 100000 100000 100017 99986.2
+ 100012 99989.1 100009 99991.6 100004 99994.8 100002 99997.1 100000 100016 99986.1 100012
+ 99988.9 100009 99991.3 100006 99993.7 100003 99996.4 100000 100016 99985.8 100012 99988.7
+ 100009 99991.4 100005 99994.3 100003 99996.8 100000 100016 99985.9 100013 99988.6 100009
+ 99991.2 100006 99993.6 100003 99996.2 100000 100016 99985.6 100012 99988.5 100009 99991.2
+ 100006 99994 100003 99996.5 100000 100016 99985.8 100013 99988.5 100010 99991.1 100007
+ 99993.6 100003 99996.1 100000 100015 99985.6 100013 99988.4 100009 99991.1 100006 99993.8
+ 100003 99996.3 100000 100016 99985.9 100013 99988.6 100010 99991.2 100006 99993.7 100003
+ 99996.2 100000 100016 99985.6 100013 99988.5 100009 99991.1 100006 99993.7 100003 99996.3
+ 100000 100016 99986.2 100013 99988.9 100009 99991.5 100006 99993.9 100003 99996.4 100000
+ 100017 99985.9 100013 99988.9 100009 99991.5 100006 99994 100003 99996.4 100000 100014
+ 99988.8 100011 99990.3 100008 99992.3 100005 99994.4 100003 99996.8 100000 100014 99988.8
+ 100011 99990.2 100008 99992.2 100005 99994.4 100003 99996.7 100000
+ </DataArray>
+ <DataArray type="Float32" Name="delP" NumberOfComponents="1" format="ascii">
+ 12.7076 -12.4945 12.7147 -12.4873 11.5107 11.5158 -10.5338 -10.5425 8.73663 8.82835 -8.18645 -8.55845
+ 3.71356 6.01525 -4.74737 -6.10381 2.06491 3.03852 -2.64549 -3.44214 0 0 17.0914 -13.8383
+ 12.1097 -10.8542 8.76432 -8.37632 4.4706 -5.18672 2.29655 -2.85157 0 16.3375 -13.8691 12.3292
+ -11.076 9.14511 -8.69825 6.32756 -6.2847 3.16635 -3.64969 0 15.9681 -14.1798 12.4226 -11.2617
+ 8.97967 -8.61866 5.23928 -5.70375 2.62092 -3.25249 0 15.8914 -14.1375 12.5676 -11.3627 9.39897
+ -8.82097 6.47907 -6.38522 3.25645 -3.81512 0 15.6062 -14.3606 12.4863 -11.4896 9.17212 -8.80719
+ 5.734 -6.04666 2.861 -3.53865 0 15.7283 -14.2093 12.6345 -11.4624 9.5261 -8.87624 6.51486
+ -6.40337 3.28666 -3.87781 0 15.4858 -14.4277 12.5174 -11.5889 9.30705 -8.90985 6.02325 -6.23511
+ 3.01126 -3.70012 0 15.8087 -14.1284 12.6719 -11.3985 9.53782 -8.81664 6.4479 -6.32766 3.24639
+ -3.82553 0 15.5736 -14.3812 12.6133 -11.5458 9.36205 -8.87881 6.12033 -6.26005 3.04599 -3.72436
+ 0 16.2371 -13.7486 12.6555 -11.0583 9.35481 -8.53976 6.20458 -6.09135 3.08707 -3.63627 0
+ 16.8246 -14.0623 12.6584 -11.1259 9.10147 -8.52942 5.87792 -6.01484 2.8939 -3.61811 0 14.0211
+ -11.1779 11.285 -9.69386 8.39836 -7.73738 5.49217 -5.58777 2.66224 -3.24038 0 13.9475 -11.2514
+ 11.2436 -9.77022 8.35548 -7.79895 5.45528 -5.63374 2.63353 -3.28055 0
+ </DataArray>
+ <DataArray type="Float32" Name="X^H2O" NumberOfComponents="1" format="ascii">
+ 0.008 0.00639858 0.00884764 0.00866148 0.00568373 0.00850738 0.00475249 0.00828663 0.00388155 0.00802916 0.00354285 0.00780191
+ 0.00561738 0.00776443 0.00701838 0.007743 0.00815446 0.00774103 0.00940333 0.00779163 0.0199569 0.00888273 0.00758332 0.00769322
+ 0.00774766 0.00779033 0.0078032 0.00779355 0.0077811 0.00776858 0.00775815 0.00775621 0.00819061 0.00752683 0.00754476 0.0075634
+ 0.0075875 0.00760789 0.00762438 0.00763296 0.00764089 0.00764682 0.00765974 0.00799899 0.00756122 0.00756113 0.00756274 0.0075666
+ 0.00757175 0.00757794 0.00758306 0.00758864 0.00759415 0.00760953 0.00789157 0.00758727 0.00758609 0.00758514 0.00758443 0.00758434
+ 0.00758484 0.00758576 0.00758707 0.00758924 0.00760374 0.00790229 0.00760326 0.00760237 0.00760149 0.00760055 0.00759976 0.00759909
+ 0.0075987 0.00759846 0.00759931 0.00761493 0.00791603 0.00761129 0.00761075 0.00761018 0.00760955 0.00760896 0.00760838 0.0076079
+ 0.00760756 0.00760848 0.0076266 0.00793841 0.00761523 0.00761435 0.0076136 0.00761294 0.00761243 0.00761204 0.00761182 0.00761183
+ 0.00761333 0.00763434 0.00794521 0.00761071 0.00761049 0.00761036 0.00761041 0.00761072 0.00761134 0.00761227 0.00761359 0.00761624
+ 0.00763889 0.00794919 0.00760424 0.00760393 0.00760459 0.00760631 0.00760889 0.00761242 0.00761642 0.00762111 0.00762659 0.0076476
+ 0.00795234 0.00757301 0.00758345 0.0075951 0.00760946 0.00762337 0.00763827 0.00765176 0.00766519 0.00767758 0.00770162 0.00800484
+ 0.00761882 0.00767841 0.00772312 0.00776918 0.00780216 0.00783126 0.00785091 0.00786636 0.00787653 0.00789583 0.00820428 0.00891555
+ 0.00884878 0.00878432 0.00870188 0.00862544 0.00853867 0.00846084 0.00837967 0.00830662 0.00825266 0.00872767 0.008 0.00733974
+ 0.0068826 0.00631456 0.0059062 0.00548498 0.00517694 0.00489343 0.00467724 0.0044644 0.008
+ </DataArray>
+ <DataArray type="Float32" Name="rho" NumberOfComponents="1" format="ascii">
+ 1.16273 1.16356 1.16213 1.16197 1.16435 1.16236 1.16475 1.16226 1.16559 1.16266 1.16563 1.16262
+ 1.1643 1.16282 1.16322 1.16269 1.1625 1.1628 1.16157 1.16269 1.15424 1.16196 1.16307 1.16264
+ 1.1629 1.1626 1.16282 1.16263 1.16279 1.16268 1.16278 1.16272 1.16245 1.1631 1.16274 1.16303
+ 1.16274 1.16296 1.16274 1.16291 1.16276 1.16287 1.16278 1.16258 1.16308 1.16273 1.16303 1.16276
+ 1.16299 1.16278 1.16294 1.1628 1.1629 1.16282 1.16266 1.16306 1.16271 1.16302 1.16274 1.16298
+ 1.16277 1.16295 1.1628 1.16291 1.16282 1.16265 1.16304 1.16269 1.16301 1.16273 1.16297 1.16276
+ 1.16293 1.16279 1.1629 1.16281 1.16264 1.16304 1.16269 1.163 1.16272 1.16297 1.16275 1.16293
+ 1.16278 1.16289 1.1628 1.16262 1.16303 1.16268 1.163 1.16272 1.16296 1.16275 1.16292 1.16278
+ 1.16289 1.1628 1.16262 1.16304 1.16269 1.163 1.16272 1.16297 1.16275 1.16293 1.16278 1.16289
+ 1.16279 1.16262 1.16304 1.16269 1.16301 1.16272 1.16297 1.16275 1.16292 1.16278 1.16288 1.16279
+ 1.16261 1.16307 1.16271 1.16301 1.16273 1.16296 1.16274 1.1629 1.16275 1.16284 1.16275 1.16258
+ 1.16305 1.16264 1.16292 1.16261 1.16283 1.1626 1.16275 1.16261 1.1627 1.16261 1.16244 1.1621
+ 1.16186 1.16216 1.16198 1.16224 1.16211 1.16232 1.16225 1.1624 1.16237 1.16207 1.16274 1.16291
+ 1.1635 1.16365 1.16415 1.16426 1.16463 1.1647 1.16495 1.16503 1.16258
+ </DataArray>
+ <DataArray type="Float32" Name="mu" NumberOfComponents="1" format="ascii">
+ 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05
+ 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05
+ 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05
+ 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05
+ 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05
+ 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05
+ 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05
+ 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05
+ 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05
+ 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05
+ 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05
+ 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05
+ 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05
+ 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05 1.82387e-05
+ </DataArray>
+ <DataArray type="Float32" Name="v" NumberOfComponents="3" format="ascii">
+ 0 0 0 0 0 0 3.5 0 0 3.52261 0.133439 0
+ 0 0 0 3.10318 0.0397147 0 0 0 0 3.32717 0.0969562 0
+ 0 0 0 2.78192 0.0977522 0 0 0 0 3.06062 0.102388 0
+ -0.115849 -0.00175612 0 2.55383 -0.0717298 0 0.32921 0.00308935 0 2.90739 0.0240212 0
+ -0.112088 -0.0019973 0 2.42871 -0.0148873 0 0.383762 0.0016687 0 2.78055 0.0462522 0
+ -0.416799 -0.00103634 0 2.21883 0.147873 0 3.5 0 0 3.97611 0.135411 0
+ 3.63171 0.0331452 0 3.96465 0.105193 0 3.62018 0.0667938 0 3.99838 0.0771861 0
+ 3.62136 0.0251646 0 3.85319 0.0462281 0 3.53092 0.029934 0 3.75459 0.0704872 0
+ 3.46606 0.162208 0 3.5 0 0 3.96827 0.0929385 0 3.64746 0.0569344 0
+ 4.0091 0.115685 0 3.73122 0.0926953 0 4.0703 0.103988 0 3.78528 -0.00884531 0
+ 4.09312 0.0518312 0 3.80518 0.0221414 0 4.10276 0.0843086 0 3.86436 0.101274 0
+ 3.5 0 0 3.91554 0.0540075 0 3.58704 0.0297174 0 3.97234 0.0668087 0
+ 3.68263 0.0384946 0 4.07347 0.0504149 0 3.76906 0.0303677 0 4.06877 0.0404479 0
+ 3.80092 0.032808 0 4.09855 0.0586242 0 3.8703 0.0664113 0 3.5 0 0
+ 3.88719 0.0369148 0 3.56405 0.0374395 0 3.93426 0.0664135 0 3.6435 0.0470753 0
+ 3.99244 0.0581781 0 3.71145 -0.00316295 0 4.05398 0.0303077 0 3.77579 0.0104698 0
+ 4.11531 0.0473784 0 3.85955 0.0566774 0 3.5 0 0 3.86659 0.0129515 0
+ 3.53694 0.00779407 0 3.9143 0.0182623 0 3.61177 0.00450286 0 3.9875 0.00911529 0
+ 3.68695 0.00574065 0 4.017 0.00687795 0 3.73889 0.00457519 0 4.06504 0.0128807 0
+ 3.81561 0.015201 0 3.5 0 0 3.85981 0.0108036 0 3.53731 0.0102988 0
+ 3.90102 0.0195087 0 3.60194 0.00621148 0 3.95428 0.00999858 0 3.66873 -0.0207267 0
+ 4.02149 -0.00651645 0 3.74162 -0.0182141 0 4.09276 -0.00166231 0 3.82756 0.0015045 0
+ 3.5 0 0 3.85133 -0.0120295 0 3.52373 -0.0193841 0 3.89305 -0.0279505 0
+ 3.58843 -0.0318146 0 3.95588 -0.0377894 0 3.65991 -0.0330108 0 4.00349 -0.0399422 0
+ 3.7231 -0.0381885 0 4.06268 -0.0436643 0 3.80329 -0.0466272 0 3.5 0 0
+ 3.86121 -0.00870985 0 3.54081 -0.0211788 0 3.9056 -0.0295636 0 3.60759 -0.0375571 0
+ 3.96355 -0.0447559 0 3.6795 -0.0556996 0 4.03222 -0.056937 0 3.75518 -0.0612776 0
+ 4.1044 -0.061733 0 3.83627 -0.0653025 0 3.5 0 0 3.85782 -0.0444244 0
+ 3.53418 -0.0548395 0 3.90939 -0.0852571 0 3.61127 -0.0779007 0 3.97607 -0.0981892 0
+ 3.69005 -0.083874 0 4.03606 -0.100417 0 3.7624 -0.0928305 0 4.10025 -0.107962 0
+ 3.83348 -0.114762 0 3.5 0 0 3.91481 -0.0504797 0 3.60303 -0.0663328 0
+ 3.97994 -0.102633 0 3.69545 -0.0925367 0 4.04625 -0.11732 0 3.77283 -0.106078 0
+ 4.10656 -0.122519 0 3.83404 -0.115238 0 4.15444 -0.123512 0 3.88891 -0.132938 0
+ 3.5 0 0 3.93779 -0.160023 0 3.62328 -0.102881 0 3.99174 -0.179873 0
+ 3.69479 -0.132112 0 4.02409 -0.172299 0 3.72755 -0.132414 0 4.02713 -0.154781 0
+ 3.7226 -0.129947 0 4.0036 -0.151507 0 3.69493 -0.169734 0 3.5 0 0
+ 3.66153 -0.162261 0 3.33193 -0.0838218 0 3.51551 -0.151727 0 3.13615 -0.0987867 0
+ 3.34681 -0.129778 0 2.94834 -0.0922446 0 3.18652 -0.103511 0 2.7813 -0.0786421 0
+ 3.0359 -0.096603 0 2.61508 -0.0181357 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>
+ </PointData>
+ <Points>
+ <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+ 0 0.25 0 0.05 0.25 0 0 0.267873 0 0.05 0.267873 0
+ 0.1 0.25 0 0.1 0.267873 0 0.15 0.25 0 0.15 0.267873 0
+ 0.2 0.25 0 0.2 0.267873 0 0.25 0.25 0 0.25 0.267873 0
+ 0.3 0.25 0 0.3 0.267873 0 0.35 0.25 0 0.35 0.267873 0
+ 0.4 0.25 0 0.4 0.267873 0 0.45 0.25 0 0.45 0.267873 0
+ 0.5 0.25 0 0.5 0.267873 0 0 0.287534 0 0.05 0.287534 0
+ 0.1 0.287534 0 0.15 0.287534 0 0.2 0.287534 0 0.25 0.287534 0
+ 0.3 0.287534 0 0.35 0.287534 0 0.4 0.287534 0 0.45 0.287534 0
+ 0.5 0.287534 0 0 0.30916 0 0.05 0.30916 0 0.1 0.30916 0
+ 0.15 0.30916 0 0.2 0.30916 0 0.25 0.30916 0 0.3 0.30916 0
+ 0.35 0.30916 0 0.4 0.30916 0 0.45 0.30916 0 0.5 0.30916 0
+ 0 0.332949 0 0.05 0.332949 0 0.1 0.332949 0 0.15 0.332949 0
+ 0.2 0.332949 0 0.25 0.332949 0 0.3 0.332949 0 0.35 0.332949 0
+ 0.4 0.332949 0 0.45 0.332949 0 0.5 0.332949 0 0 0.359117 0
+ 0.05 0.359117 0 0.1 0.359117 0 0.15 0.359117 0 0.2 0.359117 0
+ 0.25 0.359117 0 0.3 0.359117 0 0.35 0.359117 0 0.4 0.359117 0
+ 0.45 0.359117 0 0.5 0.359117 0 0 0.387902 0 0.05 0.387902 0
+ 0.1 0.387902 0 0.15 0.387902 0 0.2 0.387902 0 0.25 0.387902 0
+ 0.3 0.387902 0 0.35 0.387902 0 0.4 0.387902 0 0.45 0.387902 0
+ 0.5 0.387902 0 0 0.419565 0 0.05 0.419565 0 0.1 0.419565 0
+ 0.15 0.419565 0 0.2 0.419565 0 0.25 0.419565 0 0.3 0.419565 0
+ 0.35 0.419565 0 0.4 0.419565 0 0.45 0.419565 0 0.5 0.419565 0
+ 0 0.454395 0 0.05 0.454395 0 0.1 0.454395 0 0.15 0.454395 0
+ 0.2 0.454395 0 0.25 0.454395 0 0.3 0.454395 0 0.35 0.454395 0
+ 0.4 0.454395 0 0.45 0.454395 0 0.5 0.454395 0 0 0.492708 0
+ 0.05 0.492708 0 0.1 0.492708 0 0.15 0.492708 0 0.2 0.492708 0
+ 0.25 0.492708 0 0.3 0.492708 0 0.35 0.492708 0 0.4 0.492708 0
+ 0.45 0.492708 0 0.5 0.492708 0 0 0.534851 0 0.05 0.534851 0
+ 0.1 0.534851 0 0.15 0.534851 0 0.2 0.534851 0 0.25 0.534851 0
+ 0.3 0.534851 0 0.35 0.534851 0 0.4 0.534851 0 0.45 0.534851 0
+ 0.5 0.534851 0 0 0.58121 0 0.05 0.58121 0 0.1 0.58121 0
+ 0.15 0.58121 0 0.2 0.58121 0 0.25 0.58121 0 0.3 0.58121 0
+ 0.35 0.58121 0 0.4 0.58121 0 0.45 0.58121 0 0.5 0.58121 0
+ 0 0.632204 0 0.05 0.632204 0 0.1 0.632204 0 0.15 0.632204 0
+ 0.2 0.632204 0 0.25 0.632204 0 0.3 0.632204 0 0.35 0.632204 0
+ 0.4 0.632204 0 0.45 0.632204 0 0.5 0.632204 0 0 0.688297 0
+ 0.05 0.688297 0 0.1 0.688297 0 0.15 0.688297 0 0.2 0.688297 0
+ 0.25 0.688297 0 0.3 0.688297 0 0.35 0.688297 0 0.4 0.688297 0
+ 0.45 0.688297 0 0.5 0.688297 0 0 0.75 0 0.05 0.75 0
+ 0.1 0.75 0 0.15 0.75 0 0.2 0.75 0 0.25 0.75 0
+ 0.3 0.75 0 0.35 0.75 0 0.4 0.75 0 0.45 0.75 0
+ 0.5 0.75 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 38 49 48 38 39 50 49
+ 39 40 51 50 40 41 52 51 41 42 53 52
+ 42 43 54 53 44 45 56 55 45 46 57 56
+ 46 47 58 57 47 48 59 58 48 49 60 59
+ 49 50 61 60 50 51 62 61 51 52 63 62
+ 52 53 64 63 53 54 65 64 55 56 67 66
+ 56 57 68 67 57 58 69 68 58 59 70 69
+ 59 60 71 70 60 61 72 71 61 62 73 72
+ 62 63 74 73 63 64 75 74 64 65 76 75
+ 66 67 78 77 67 68 79 78 68 69 80 79
+ 69 70 81 80 70 71 82 81 71 72 83 82
+ 72 73 84 83 73 74 85 84 74 75 86 85
+ 75 76 87 86 77 78 89 88 78 79 90 89
+ 79 80 91 90 80 81 92 91 81 82 93 92
+ 82 83 94 93 83 84 95 94 84 85 96 95
+ 85 86 97 96 86 87 98 97 88 89 100 99
+ 89 90 101 100 90 91 102 101 91 92 103 102
+ 92 93 104 103 93 94 105 104 94 95 106 105
+ 95 96 107 106 96 97 108 107 97 98 109 108
+ 99 100 111 110 100 101 112 111 101 102 113 112
+ 102 103 114 113 103 104 115 114 104 105 116 115
+ 105 106 117 116 106 107 118 117 107 108 119 118
+ 108 109 120 119 110 111 122 121 111 112 123 122
+ 112 113 124 123 113 114 125 124 114 115 126 125
+ 115 116 127 126 116 117 128 127 117 118 129 128
+ 118 119 130 129 119 120 131 130 121 122 133 132
+ 122 123 134 133 123 124 135 134 124 125 136 135
+ 125 126 137 136 126 127 138 137 127 128 139 138
+ 128 129 140 139 129 130 141 140 130 131 142 141
+ 132 133 144 143 133 134 145 144 134 135 146 145
+ 135 136 147 146 136 137 148 147 137 138 149 148
+ 138 139 150 149 139 140 151 150 140 141 152 151
+ 141 142 153 152 143 144 155 154 144 145 156 155
+ 145 146 157 156 146 147 158 157 147 148 159 158
+ 148 149 160 159 149 150 161 160 150 151 162 161
+ 151 152 163 162 152 153 164 163
+ </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
+ </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
+ </DataArray>
+ </Cells>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>
diff --git a/test/references/2czeroeq2p2c-pm-reference.vtu b/test/references/2czeroeq2p2c-pm-reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..23a98f9b68b350a6e99a3e3dfba4cb728f043df1
--- /dev/null
+++ b/test/references/2czeroeq2p2c-pm-reference.vtu
@@ -0,0 +1,277 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+ <UnstructuredGrid>
+ <Piece NumberOfCells="55" NumberOfPoints="72">
+ <PointData Scalars="Sn" Vectors="velocityW">
+ <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
+ 0.0140853 0.0140839 0.0243754 0.0243607 0.0140812 0.024332 0.0140789 0.024308 0.0140775 0.0242924 0.014077 0.0242871
+ 0.115755 0.115683 0.115543 0.115426 0.115349 0.115323 0.34387 0.343736 0.343462 0.343232 0.343084 0.343031
+ 0.598441 0.598316 0.598049 0.597839 0.597697 0.597641 0.768045 0.767974 0.767786 0.767656 0.767543 0.76748
+ 0.861092 0.861128 0.861028 0.860967 0.860838 0.860708 0.911419 0.91173 0.911748 0.911738 0.911514 0.911146
+ 0.940176 0.941222 0.941415 0.941453 0.941026 0.93998 0.957922 0.960868 0.961251 0.961356 0.96063 0.957807
+ 0.969141 0.97723 0.977527 0.977785 0.976824 0.969236 0.965501 1 1 1 1 0.965964
+ </DataArray>
+ <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
+ 0.985915 0.985916 0.975625 0.975639 0.985919 0.975668 0.985921 0.975692 0.985923 0.975708 0.985923 0.975713
+ 0.884245 0.884317 0.884457 0.884574 0.884651 0.884677 0.65613 0.656264 0.656538 0.656768 0.656916 0.656969
+ 0.401559 0.401684 0.401951 0.402161 0.402303 0.402359 0.231955 0.232026 0.232214 0.232344 0.232457 0.23252
+ 0.138908 0.138872 0.138972 0.139033 0.139162 0.139292 0.088581 0.0882699 0.0882516 0.0882616 0.0884865 0.0888542
+ 0.0598236 0.0587778 0.0585848 0.0585468 0.0589737 0.0600202 0.0420784 0.0391316 0.0387488 0.038644 0.0393696 0.0421927
+ 0.0308594 0.0227695 0.0224725 0.0222153 0.0231765 0.0307642 0.0344987 0 0 0 0 0.0340357
+ </DataArray>
+ <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
+ 100003 100003 100003 100002 100003 100002 100002 100002 100002 100002 100002 100002
+ 100002 100002 100002 100002 100002 100001 100002 100002 100002 100001 100001 100001
+ 100002 100002 100001 100001 100001 100001 100001 100001 100001 100001 100001 100000
+ 100001 100001 100001 100000 100000 100000 100001 100001 100000 100000 100000 99999.9
+ 100001 100001 99999.7 100000 99999.7 99999.7 100001 100001 99999 100000 99999.3 99999.6
+ 100001 100002 99997.7 100001 99998.7 99999.6 100001 100003 99995.6 100002 99997.6 100000
+ </DataArray>
+ <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
+ 99472.7 99472.7 99130.5 99130.5 99472.7 99130.5 99472.7 99130.5 99472.7 99130.5 99472.7 99130.5
+ 98819.3 98819.3 98819.3 98819.3 98819.3 98819.3 98536.5 98536.5 98536.5 98536.5 98536.5 98536.5
+ 98279.3 98279.3 98279.3 98279.3 98279.3 98279.3 98045.4 98045.4 98045.3 98045.3 98045.4 98045.4
+ 97832.4 97832 97831.8 97831.8 97832 97832.4 97636.7 97634.8 97634 97634 97634.9 97636.7
+ 97453.2 97444.2 97441.3 97441.3 97444.5 97453.1 97273.5 97233.9 97226.3 97226 97235.3 97273.1
+ 97098.4 96906 96893 96888.1 96914.9 97094.9 97163 95700.5 95693.1 95699.4 95695.1 97154.4
+ </DataArray>
+ <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
+ 530.282 530.154 872.149 872.018 529.901 871.761 529.685 871.546 529.546 871.406 529.499 871.358
+ 1182.93 1182.8 1182.54 1182.33 1182.19 1182.14 1465.47 1465.34 1465.07 1464.84 1464.69 1464.64
+ 1722.37 1722.23 1721.93 1721.69 1721.53 1721.47 1956.04 1955.91 1955.58 1955.35 1955.15 1955.04
+ 2169.01 2169.12 2168.81 2168.63 2168.24 2167.84 2364.59 2366.17 2366.27 2366.21 2365.07 2363.2
+ 2548.18 2556.8 2558.41 2558.73 2555.17 2546.58 2727.9 2767.27 2772.67 2774.16 2763.95 2726.44
+ 2902.88 3095.76 3104.7 3112.58 3083.79 2904.72 2837.83 4302.49 4302.49 4302.49 4302.49 2845.58
+ </DataArray>
+ <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
+ 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047
+ 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047
+ 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047
+ 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047 997.047
+ 997.047 997.047 997.047 997.047 997.047 997.047 997.046 997.046 997.046 997.046 997.046 997.046
+ 997.046 997.046 997.046 997.046 997.046 997.046 997.046 997.046 997.046 997.046 997.046 997.046
+ </DataArray>
+ <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
+ 1.15427 1.15427 1.15427 1.15427 1.15427 1.15426 1.15427 1.15426 1.15426 1.15426 1.15426 1.15426
+ 1.15426 1.15426 1.15426 1.15426 1.15426 1.15425 1.15426 1.15426 1.15426 1.15425 1.15425 1.15425
+ 1.15426 1.15426 1.15425 1.15425 1.15425 1.15425 1.15425 1.15425 1.15425 1.15425 1.15424 1.15424
+ 1.15425 1.15425 1.15424 1.15424 1.15424 1.15424 1.15425 1.15425 1.15424 1.15424 1.15424 1.15424
+ 1.15425 1.15425 1.15423 1.15424 1.15423 1.15423 1.15425 1.15425 1.15423 1.15424 1.15423 1.15423
+ 1.15425 1.15426 1.15421 1.15424 1.15422 1.15423 1.15425 1.16429 1.16322 1.16249 1.16157 1.15424
+ </DataArray>
+ <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
+ 1122.18 1122.18 1107.45 1107.48 1122.19 1107.54 1122.19 1107.59 1122.19 1107.63 1122.19 1107.64
+ 734.377 734.578 734.969 735.298 735.513 735.587 285.574 285.753 286.118 286.423 286.622 286.692
+ 65.7264 65.7869 65.9163 66.0183 66.0869 66.114 13.0872 13.0989 13.1302 13.1517 13.1705 13.1811
+ 2.88784 2.88562 2.89178 2.89557 2.90355 2.91159 0.752847 0.744891 0.744426 0.744681 0.750423 0.759879
+ 0.226654 0.214569 0.212385 0.211958 0.216799 0.228973 0.074536 0.0589032 0.0570449 0.0565427 0.0600782 0.0751911
+ 0.0267663 0.0092188 0.0087878 0.00842538 0.00983122 0.0264872 0.0389117 0 0 0 0 0.0372013
+ </DataArray>
+ <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
+ 0.970566 0.969906 12.6741 12.6477 0.968603 12.5962 0.967496 12.5531 0.966782 12.5253 0.966536 12.5157
+ 737.705 736.691 734.716 733.06 731.98 731.604 7430.79 7424.87 7412.82 7402.72 7396.18 7393.86
+ 22392.5 22383.5 22364.2 22349 22338.8 22334.8 35810.4 35804.4 35788.4 35777.5 35767.9 35762.5
+ 43865.9 43869 43860.3 43854.9 43843.6 43832.2 48404 48436.1 48437.9 48436.9 48413.7 48375.9
+ 51436.7 51542.9 51562.4 51566.2 51523 51416.6 53099.9 53338.8 53368.9 53377.1 53320 53090.4
+ 53933.3 54387.3 54401.3 54413.2 54367.8 53939.4 53686.5 54828.6 54828.6 54828.6 54828.6 53719.3
+ </DataArray>
+ <DataArray type="Float32" Name="X_liquid^H2O" NumberOfComponents="1" format="ascii">
+ 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978
+ 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978
+ 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978
+ 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978
+ 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978
+ 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.283916 0.354401 0.411513 0.47416 0.999978
+ </DataArray>
+ <DataArray type="Float32" Name="X_liquid^Air" NumberOfComponents="1" format="ascii">
+ 2.15475e-05 2.15475e-05 2.15474e-05 2.15474e-05 2.15474e-05 2.15473e-05 2.15474e-05 2.15473e-05 2.15474e-05 2.15473e-05 2.15473e-05 2.15473e-05
+ 2.15474e-05 2.15473e-05 2.15473e-05 2.15472e-05 2.15472e-05 2.15472e-05 2.15473e-05 2.15473e-05 2.15472e-05 2.15471e-05 2.15471e-05 2.15471e-05
+ 2.15472e-05 2.15472e-05 2.15471e-05 2.15471e-05 2.1547e-05 2.1547e-05 2.15472e-05 2.15471e-05 2.15471e-05 2.1547e-05 2.1547e-05 2.1547e-05
+ 2.15472e-05 2.15471e-05 2.1547e-05 2.15469e-05 2.15469e-05 2.15469e-05 2.15471e-05 2.15471e-05 2.15469e-05 2.15469e-05 2.15468e-05 2.15468e-05
+ 2.15471e-05 2.15471e-05 2.15468e-05 2.15469e-05 2.15468e-05 2.15468e-05 2.15472e-05 2.15471e-05 2.15466e-05 2.15469e-05 2.15467e-05 2.15468e-05
+ 2.15471e-05 2.15472e-05 2.15463e-05 2.1547e-05 2.15466e-05 2.15468e-05 2.1547e-05 2.20521e-05 2.20009e-05 2.19621e-05 2.1917e-05 2.15469e-05
+ </DataArray>
+ <DataArray type="Float32" Name="X_gas^H2O" NumberOfComponents="1" format="ascii">
+ 0.0199563 0.0199563 0.0199564 0.0199564 0.0199564 0.0199564 0.0199564 0.0199565 0.0199564 0.0199565 0.0199564 0.0199565
+ 0.0199564 0.0199565 0.0199565 0.0199565 0.0199566 0.0199566 0.0199565 0.0199565 0.0199566 0.0199566 0.0199566 0.0199567
+ 0.0199565 0.0199566 0.0199566 0.0199567 0.0199567 0.0199567 0.0199566 0.0199566 0.0199567 0.0199567 0.0199568 0.0199568
+ 0.0199566 0.0199567 0.0199568 0.0199568 0.0199568 0.0199568 0.0199566 0.0199567 0.0199568 0.0199568 0.0199569 0.0199569
+ 0.0199566 0.0199567 0.0199569 0.0199569 0.019957 0.0199569 0.0199566 0.0199566 0.0199571 0.0199569 0.019957 0.019957
+ 0.0199566 0.0199565 0.0199574 0.0199568 0.0199572 0.019957 0.0199567 0.00561738 0.00701838 0.00815446 0.00940333 0.0199569
+ </DataArray>
+ <DataArray type="Float32" Name="X_gas^Air" NumberOfComponents="1" format="ascii">
+ 0.980044 0.980044 0.980044 0.980044 0.980044 0.980044 0.980044 0.980044 0.980044 0.980043 0.980044 0.980043
+ 0.980044 0.980044 0.980043 0.980043 0.980043 0.980043 0.980044 0.980043 0.980043 0.980043 0.980043 0.980043
+ 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043
+ 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043
+ 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043
+ 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.980043 0.994383 0.992982 0.991846 0.990597 0.980043
+ </DataArray>
+ <DataArray type="Float32" Name="x_liquid^H2O" NumberOfComponents="1" format="ascii">
+ 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987
+ 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987
+ 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987
+ 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987
+ 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987
+ 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.283924 0.354409 0.411521 0.474168 0.999987
+ </DataArray>
+ <DataArray type="Float32" Name="x_liquid^Air" NumberOfComponents="1" format="ascii">
+ 1.34042e-05 1.34042e-05 1.34041e-05 1.34041e-05 1.34041e-05 1.34041e-05 1.34041e-05 1.34041e-05 1.34041e-05 1.3404e-05 1.34041e-05 1.3404e-05
+ 1.34041e-05 1.34041e-05 1.3404e-05 1.3404e-05 1.3404e-05 1.3404e-05 1.34041e-05 1.3404e-05 1.3404e-05 1.3404e-05 1.34039e-05 1.34039e-05
+ 1.3404e-05 1.3404e-05 1.3404e-05 1.34039e-05 1.34039e-05 1.34039e-05 1.3404e-05 1.3404e-05 1.34039e-05 1.34039e-05 1.34039e-05 1.34038e-05
+ 1.3404e-05 1.34039e-05 1.34039e-05 1.34038e-05 1.34038e-05 1.34038e-05 1.3404e-05 1.34039e-05 1.34038e-05 1.34038e-05 1.34038e-05 1.34038e-05
+ 1.3404e-05 1.34039e-05 1.34037e-05 1.34038e-05 1.34037e-05 1.34037e-05 1.3404e-05 1.3404e-05 1.34036e-05 1.34038e-05 1.34037e-05 1.34037e-05
+ 1.3404e-05 1.3404e-05 1.34035e-05 1.34039e-05 1.34036e-05 1.34037e-05 1.34039e-05 1.37184e-05 1.36864e-05 1.36622e-05 1.36342e-05 1.34038e-05
+ </DataArray>
+ <DataArray type="Float32" Name="x_gas^H2O" NumberOfComponents="1" format="ascii">
+ 0.0316961 0.0316961 0.0316962 0.0316963 0.0316962 0.0316963 0.0316963 0.0316964 0.0316963 0.0316964 0.0316963 0.0316965
+ 0.0316963 0.0316964 0.0316964 0.0316965 0.0316966 0.0316966 0.0316964 0.0316965 0.0316966 0.0316966 0.0316967 0.0316967
+ 0.0316965 0.0316966 0.0316967 0.0316967 0.0316968 0.0316968 0.0316966 0.0316966 0.0316968 0.0316968 0.0316969 0.0316969
+ 0.0316966 0.0316967 0.0316969 0.0316969 0.031697 0.031697 0.0316966 0.0316967 0.031697 0.031697 0.0316971 0.0316971
+ 0.0316966 0.0316967 0.0316971 0.031697 0.0316972 0.0316971 0.0316966 0.0316967 0.0316974 0.031697 0.0316973 0.0316972
+ 0.0316966 0.0316965 0.0316978 0.0316968 0.0316975 0.0316972 0.0316968 0.00899941 0.0112344 0.0130439 0.0150303 0.031697
+ </DataArray>
+ <DataArray type="Float32" Name="x_gas^Air" NumberOfComponents="1" format="ascii">
+ 0.968304 0.968304 0.968304 0.968304 0.968304 0.968304 0.968304 0.968304 0.968304 0.968304 0.968304 0.968304
+ 0.968304 0.968304 0.968304 0.968304 0.968303 0.968303 0.968304 0.968304 0.968303 0.968303 0.968303 0.968303
+ 0.968304 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303
+ 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303
+ 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303 0.968303
+ 0.968303 0.968304 0.968302 0.968303 0.968303 0.968303 0.968303 0.991001 0.988766 0.986956 0.98497 0.968303
+ </DataArray>
+ <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
+ 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41
+ 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41
+ 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41
+ 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41
+ 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41
+ 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41
+ </DataArray>
+ <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
+ 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15
+ 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15
+ 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15
+ 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15
+ 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15
+ 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15 298.15
+ </DataArray>
+ <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
+ 3 3 3 3 3 3 3 3 3 3 3 3
+ 3 3 3 3 3 3 3 3 3 3 3 3
+ 3 3 3 3 3 3 3 3 3 3 3 3
+ 3 3 3 3 3 3 3 3 3 3 3 3
+ 3 3 3 3 3 3 3 3 3 3 3 3
+ 3 3 3 3 3 3 3 2 2 2 2 3
+ </DataArray>
+ <DataArray type="Float32" Name="velocityW" NumberOfComponents="3" format="ascii">
+ 3.10381e-09 -1.86032e-09 0 4.06777e-09 -3.62789e-10 0 4.52264e-09 -1.74433e-09 0 5.8113e-09 1.88955e-09 0
+ 2.46815e-09 2.10586e-09 0 3.49318e-09 7.60247e-09 0 -2.62968e-09 2.08528e-09 0 -3.6771e-09 7.54959e-09 0
+ -4.15043e-09 -3.97595e-10 0 -5.8919e-09 1.8046e-09 0 -3.13693e-09 -1.88365e-09 0 -4.54319e-09 -1.79452e-09 0
+ 6.94243e-09 8.05697e-09 0 8.27558e-09 1.71883e-08 0 4.73953e-09 2.98709e-08 0 -4.93102e-09 2.97373e-08 0
+ -8.33281e-09 1.69921e-08 0 -6.93327e-09 7.95837e-09 0 7.32583e-09 3.7112e-08 0 8.08261e-09 5.31766e-08 0
+ 4.35641e-09 7.24513e-08 0 -4.53139e-09 7.2228e-08 0 -8.10286e-09 5.28827e-08 0 -7.26953e-09 3.70301e-08 0
+ 6.33295e-09 7.22187e-08 0 6.19503e-09 9.50427e-08 0 2.97302e-09 1.18659e-07 0 -3.12602e-09 1.1837e-07 0
+ -6.19717e-09 9.47028e-08 0 -6.25338e-09 7.22544e-08 0 6.43265e-09 9.59252e-08 0 5.56601e-09 1.25196e-07 0
+ 2.29565e-09 1.50992e-07 0 -2.45353e-09 1.50661e-07 0 -5.56479e-09 1.24827e-07 0 -6.3306e-09 9.61195e-08 0
+ 7.78875e-09 1.06514e-07 0 6.05569e-09 1.42907e-07 0 2.1079e-09 1.69619e-07 0 -2.28851e-09 1.69259e-07 0
+ -6.05839e-09 1.42476e-07 0 -7.64658e-09 1.06866e-07 0 1.02872e-08 1.08487e-07 0 7.22141e-09 1.53668e-07 0
+ 2.03625e-09 1.80249e-07 0 -2.24648e-09 1.79897e-07 0 -7.23882e-09 1.53115e-07 0 -1.00678e-08 1.0901e-07 0
+ 1.42747e-08 1.04834e-07 0 8.98346e-09 1.61494e-07 0 1.84977e-09 1.86577e-07 0 -2.08683e-09 1.86327e-07 0
+ -9.03024e-09 1.60718e-07 0 -1.38795e-08 1.0559e-07 0 2.16512e-08 9.62593e-08 0 1.20328e-08 1.68955e-07 0
+ 1.33199e-09 1.90267e-07 0 -1.58444e-09 1.90333e-07 0 -1.20573e-08 1.67734e-07 0 -2.06963e-08 9.74483e-08 0
+ 4.71047e-08 7.71452e-08 0 2.38369e-08 1.692e-07 0 3.53107e-10 1.87898e-07 0 -4.71711e-10 1.83816e-07 0
+ -2.31524e-08 1.74893e-07 0 -4.52241e-08 7.94672e-08 0 2.36137e-07 6.35992e-08 0 1.18069e-07 1.65533e-07 0
+ 0 1.84231e-07 0 0 1.75768e-07 0 -1.12336e-07 1.7844e-07 0 -2.24672e-07 6.67604e-08 0
+ </DataArray>
+ <DataArray type="Float32" Name="velocityN" NumberOfComponents="3" format="ascii">
+ 6.65527e-10 -7.77578e-10 0 9.87679e-10 -5.3726e-10 0 8.84014e-09 -9.79259e-09 0 1.30401e-08 -6.45179e-09 0
+ 1.20768e-09 -2.40279e-10 0 1.5815e-08 -2.39377e-09 0 9.07379e-10 -1.73045e-10 0 1.18024e-08 -1.54563e-09 0
+ 4.7598e-10 -1.48442e-10 0 6.18097e-09 -1.22053e-09 0 2.42732e-10 -1.41791e-10 0 3.14712e-09 -1.13734e-09 0
+ 5.22306e-07 -3.49776e-07 0 7.70695e-07 -1.76963e-07 0 9.32476e-07 3.41563e-08 0 6.93205e-07 7.58981e-08 0
+ 3.62552e-07 9.19486e-08 0 1.84562e-07 9.61057e-08 0 5.47481e-06 -3.72401e-06 0 8.23758e-06 -1.88796e-06 0
+ 9.97164e-06 6.15181e-07 0 7.26592e-06 1.0563e-06 0 3.74706e-06 1.17943e-06 0 1.90519e-06 1.2225e-06 0
+ 1.85232e-05 -1.60647e-05 0 2.81583e-05 -8.27875e-06 0 3.27752e-05 3.63809e-06 0 2.24806e-05 4.58818e-06 0
+ 1.15512e-05 4.99031e-06 0 5.89821e-06 5.19116e-06 0 3.62742e-05 -4.18775e-05 0 5.55921e-05 -2.19169e-05 0
+ 5.87063e-05 1.34327e-05 0 3.55434e-05 1.04024e-05 0 1.90989e-05 1.21194e-05 0 9.61364e-06 1.25074e-05 0
+ 5.68274e-05 -8.19032e-05 0 9.12839e-05 -4.49265e-05 0 8.24469e-05 3.85259e-05 0 3.95387e-05 1.35007e-05 0
+ 2.54479e-05 2.24009e-05 0 1.0972e-05 2.14403e-05 0 7.6703e-05 -0.000134297 0 0.000143492 -8.41305e-05 0
+ 0.000105261 9.55282e-05 0 3.14239e-05 1.8791e-06 0 3.42354e-05 3.8966e-05 0 5.86438e-06 2.82625e-05 0
+ 8.05915e-05 -0.000190896 0 0.000224798 -0.000158102 0 0.000128241 0.000216531 0 5.34087e-06 -4.95329e-05 0
+ 5.20097e-05 7.3689e-05 0 -1.91851e-05 2.41316e-05 0 2.52377e-05 -0.000226901 0 0.000346918 -0.000306389 0
+ 0.000147862 0.000459181 0 -4.94446e-05 -0.000189107 0 8.81045e-05 0.000155098 0 -9.77756e-05 -1.22691e-05 0
+ -0.000177105 -0.000185125 0 0.000526485 -0.000590764 0 0.000160483 0.000924177 0 -0.000147491 -0.000492576 0
+ 0.00015734 0.000343024 0 -0.000299449 -0.000125361 0 -0.00049421 -0.000135058 0 0.000707074 -0.000775835 0
+ 0.000162631 0.0012283 0 -0.000257662 -0.000696346 0 0.000241856 0.000473668 0 -0.000584061 -0.000208449 0
+ </DataArray>
+ </PointData>
+ <CellData Scalars="process rank">
+ <DataArray type="Float32" Name="process rank" 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
+ </DataArray>
+ </CellData>
+ <Points>
+ <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+ 0.25 0 0 0.3 0 0 0.25 0.0349916 0 0.3 0.0349916 0
+ 0.35 0 0 0.35 0.0349916 0 0.4 0 0 0.4 0.0349916 0
+ 0.45 0 0 0.45 0.0349916 0 0.5 0 0 0.5 0.0349916 0
+ 0.25 0.0668022 0 0.3 0.0668022 0 0.35 0.0668022 0 0.4 0.0668022 0
+ 0.45 0.0668022 0 0.5 0.0668022 0 0.25 0.0957209 0 0.3 0.0957209 0
+ 0.35 0.0957209 0 0.4 0.0957209 0 0.45 0.0957209 0 0.5 0.0957209 0
+ 0.25 0.122011 0 0.3 0.122011 0 0.35 0.122011 0 0.4 0.122011 0
+ 0.45 0.122011 0 0.5 0.122011 0 0.25 0.14591 0 0.3 0.14591 0
+ 0.35 0.14591 0 0.4 0.14591 0 0.45 0.14591 0 0.5 0.14591 0
+ 0.25 0.167637 0 0.3 0.167637 0 0.35 0.167637 0 0.4 0.167637 0
+ 0.45 0.167637 0 0.5 0.167637 0 0.25 0.187389 0 0.3 0.187389 0
+ 0.35 0.187389 0 0.4 0.187389 0 0.45 0.187389 0 0.5 0.187389 0
+ 0.25 0.205345 0 0.3 0.205345 0 0.35 0.205345 0 0.4 0.205345 0
+ 0.45 0.205345 0 0.5 0.205345 0 0.25 0.221669 0 0.3 0.221669 0
+ 0.35 0.221669 0 0.4 0.221669 0 0.45 0.221669 0 0.5 0.221669 0
+ 0.25 0.236509 0 0.3 0.236509 0 0.35 0.236509 0 0.4 0.236509 0
+ 0.45 0.236509 0 0.5 0.236509 0 0.25 0.25 0 0.3 0.25 0
+ 0.35 0.25 0 0.4 0.25 0 0.45 0.25 0 0.5 0.25 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 2 3 13 12
+ 3 5 14 13 5 7 15 14 7 9 16 15
+ 9 11 17 16 12 13 19 18 13 14 20 19
+ 14 15 21 20 15 16 22 21 16 17 23 22
+ 18 19 25 24 19 20 26 25 20 21 27 26
+ 21 22 28 27 22 23 29 28 24 25 31 30
+ 25 26 32 31 26 27 33 32 27 28 34 33
+ 28 29 35 34 30 31 37 36 31 32 38 37
+ 32 33 39 38 33 34 40 39 34 35 41 40
+ 36 37 43 42 37 38 44 43 38 39 45 44
+ 39 40 46 45 40 41 47 46 42 43 49 48
+ 43 44 50 49 44 45 51 50 45 46 52 51
+ 46 47 53 52 48 49 55 54 49 50 56 55
+ 50 51 57 56 51 52 58 57 52 53 59 58
+ 54 55 61 60 55 56 62 61 56 57 63 62
+ 57 58 64 63 58 59 65 64 60 61 67 66
+ 61 62 68 67 62 63 69 68 63 64 70 69
+ 64 65 71 70
+ </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
+ </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
+ </DataArray>
+ </Cells>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>