[implicit] add a test that uses the general setDirichlet method

The new test employs the function
setDirichlet(pvIdx, eqIdx)
with pvIdx != eqIdx.

On the right boundary, a Dirichlet value for the nonwetting-phase
pressure should be set. At the same time, nonwetting phase should
be allowed to leave the domain by means of an outflow boundary
condition. In order to achieve this the general 'setDirichlet'
method is used in the function 'boundaryTypesAtPos' which allows
to set the value of a primary variable by replacing an arbitrary
balance equation. In this case, the index of the primary variable
is Indices::pnIdx, while the index of the replaced equation is
Indices::contiWEqIdx. This allows to set the outflow condition for
equation index Indices::contiNEqIdx.

This completes FS#270.

Reviewed by Timo.



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@15167 2fb0f335-1f38-0410-981e-8018bf24f1b0

test/implicit/2p/CMakeLists.txt

@@ -29,6 +29,13 @@ add_dumux_test(test_ccadaptive2p test_ccadaptive2p test_ccadaptive2p.cc
                          ${CMAKE_CURRENT_BINARY_DIR}/lensccadaptive-00018.vtu
                  --command "${CMAKE_CURRENT_BINARY_DIR}/test_ccadaptive2p")
 
+add_dumux_test(test_generalizeddirichlet test_generalizeddirichlet test_generalizeddirichlet.cc
+               python ${CMAKE_SOURCE_DIR}/bin/runtest.py
+                 --script fuzzy
+                 --files ${CMAKE_SOURCE_DIR}/test/references/generalizeddirichlet-reference.vtp
+                         ${CMAKE_CURRENT_BINARY_DIR}/generalizeddirichlet-00035.vtp
+                 --command "${CMAKE_CURRENT_BINARY_DIR}/test_generalizeddirichlet")
+
 # non-isothermal tests
 if(HAVE_UG)
 add_dumux_test(test_box2pni test_box2pni test_box2pni.cc

test/implicit/2p/generalizeddirichletproblem.hh

+// -*- 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 Problem that uses a generalized Dirichlet boundary condition.
+ */
+#ifndef DUMUX_GENERALIZED_DIRICHLET_PROBLEM_HH
+#define DUMUX_GENERALIZED_DIRICHLET_PROBLEM_HH
+
+// The numerical model
+#include <dumux/implicit/2p/2pmodel.hh>
+
+// The base porous media box problem
+#include <dumux/implicit/common/implicitporousmediaproblem.hh>
+
+// The DUNE grid used
+#include <dune/grid/yaspgrid.hh>
+
+// Spatially dependent parameters
+#include "generalizeddirichletspatialparams.hh"
+
+// The components that are used
+#include <dumux/material/components/h2o.hh>
+#include <dumux/material/components/lnapl.hh>
+#include <dumux/io/cubegridcreator.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+
+namespace Dumux{
+// Forward declaration of the problem class
+template <class TypeTag>
+class GeneralizedDirichletProblem;
+
+namespace Properties {
+// Create a new type tag for the problem
+NEW_TYPE_TAG(GeneralizedDirichletProblem, INHERITS_FROM(BoxTwoP, GeneralizedDirichletSpatialParams));
+
+// Set the "Problem" property
+SET_PROP(GeneralizedDirichletProblem, Problem)
+{ typedef Dumux::GeneralizedDirichletProblem<TypeTag> type;};
+
+// Set grid and the grid creator to be used
+SET_TYPE_PROP(GeneralizedDirichletProblem, Grid, Dune::YaspGrid<1>);
+SET_TYPE_PROP(GeneralizedDirichletProblem, GridCreator, Dumux::CubeGridCreator<TypeTag>);
+
+// Set the wetting phase
+SET_PROP(GeneralizedDirichletProblem, WettingPhase)
+{
+private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+public: typedef Dumux::LiquidPhase<Scalar, Dumux::H2O<Scalar> > type;
+};
+
+// Set the non-wetting phase
+SET_PROP(GeneralizedDirichletProblem, NonwettingPhase)
+{
+private: typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+public: typedef Dumux::LiquidPhase<Scalar, Dumux::LNAPL<Scalar> > type;
+};
+
+SET_INT_PROP(GeneralizedDirichletProblem, Formulation, TwoPFormulation::pnsw);
+
+SET_BOOL_PROP(GeneralizedDirichletProblem, ProblemEnableGravity, false);
+}
+
+/*!
+ * \ingroup TwoPBoxModel
+ *
+ * \brief  Problem that uses a generalized Dirichlet boundary condition.
+ *
+ * On the right boundary, a Dirichlet value for the nonwetting-phase pressure
+ * should be set. At the same time, nonwetting phase should be allowed to leave
+ * the domain by means of an outflow boundary condition. In order to achieve
+ * this the general 'setDirichlet' method is used in the function
+ * 'boundaryTypesAtPos' which allows to set the value of a primary variable by
+ * replacing an arbitrary balance equation. In this case, the index of the
+ * primary variable is Indices::pnIdx, while the index of the replaced equation
+ * is Indices::contiWEqIdx. This allows to set the outflow condition for
+ * equation index Indices::contiNEqIdx.
+ */
+template <class TypeTag>
+class GeneralizedDirichletProblem : public ImplicitPorousMediaProblem<TypeTag>
+{
+    typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+    // Grid dimension
+    enum { dim = GridView::dimension,
+           dimWorld = GridView::dimensionworld
+    };
+
+    // Types from DUNE-Grid
+    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+    // Dumux specific types
+    typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+
+public:
+    GeneralizedDirichletProblem(TimeManager &timeManager,
+                                const GridView &gridView)
+    : ParentType(timeManager, gridView)
+    , eps_(3e-6)
+    {
+        name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
+    }
+
+    //! Specifies the problem name. This is used as a prefix for files
+    //! generated by the simulation.
+    const char *name() const
+    { return name_.c_str(); }
+
+    //! Returns the temperature within a finite volume. We use constant
+    //! 10 degrees Celsius.
+    Scalar temperature() const
+    { return 283.15; };
+
+    bool shouldWriteRestartFile() const
+    {
+        return false;
+    }
+
+    //! Specifies which kind of boundary condition should be used for
+    //! which equation for a finite volume on the boundary.
+    void boundaryTypesAtPos(BoundaryTypes &bcTypes,
+                            const GlobalPosition &globalPos) const
+    {
+        if (globalPos[0] > this->bBoxMax()[0] - eps_)
+        {
+           bcTypes.setDirichlet(Indices::pnIdx, Indices::contiWEqIdx);
+           bcTypes.setOutflow(Indices::contiNEqIdx);
+        }
+        else // Neumann for the remaining boundaries
+           bcTypes.setAllNeumann();
+
+    }
+
+    //! Evaluates the Dirichlet boundary conditions for a finite volume
+    //! on the grid boundary. Here, the 'values' parameter stores
+    //! primary variables.
+    void dirichletAtPos(PrimaryVariables &values,
+                        const GlobalPosition &globalPos) const
+    {
+        values = 0;
+        if (globalPos[0] > this->bBoxMax()[0] - eps_)
+            values[Indices::pnIdx] = 200e3; // 200 kPa = 2 bar oil pressure on right boundary
+    }
+
+    //! Evaluates the boundary conditions for a Neumann boundary
+    //! segment. Here, the 'values' parameter stores the mass flux in
+    //! [kg/(m^2 * s)] in normal direction of each phase. Negative
+    //! values mean influx.
+    void neumannAtPos(PrimaryVariables &values,
+                      const GlobalPosition &globalPos) const
+    {
+        // water influx of 0.05 kg / (m.s) on the left boundary
+        if (globalPos[0] < eps_) {
+            values[Indices::contiWEqIdx] = -5e-2;
+            values[Indices::contiNEqIdx] = 0;
+        }
+    }
+
+    //! Evaluates the initial value for a control volume. For this
+    //! method, the 'values' parameter stores primary variables.
+    void initialAtPos(PrimaryVariables &values,
+                      const GlobalPosition &globalPos) const
+    {
+        values[Indices::pnIdx] = 200e3; // 200 kPa = 2 bar
+        values[Indices::swIdx] = 0;
+    }
+
+    //! Evaluates the source term for all phases within a given
+    //! sub-control-volume. In this case, the 'values' parameter
+    //! stores the rate mass generated or annihilated per volume unit
+    //! in [kg / (m^3 * s)]. Positive values mean that mass is created.
+    void sourceAtPos(PrimaryVariables &values,
+                     const GlobalPosition &globalPos) const
+    {
+        values = 0;
+    }
+
+private:
+    Scalar eps_;
+    std::string name_;
+};
+}
+
+#endif

test/implicit/2p/generalizeddirichletspatialparams.hh

+// -*- 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 spatial parameters for the problem that uses a generalized
+ *        Dirichlet boundary condition.
+ */
+#ifndef DUMUX_GENERALIZED_DIRICHLET_SPATIAL_PARAMS_COUPLED_HH
+#define DUMUX_GENERALIZED_DIRICHLET_SPATIAL_PARAMS_COUPLED_HH
+
+// include parent spatialparameters
+#include <dumux/material/spatialparams/implicitspatialparams.hh>
+
+// include material laws
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+
+namespace Dumux {
+//forward declaration
+template<class TypeTag>
+class GeneralizedDirichletSpatialParams;
+
+namespace Properties
+{
+// The spatial parameters TypeTag
+NEW_TYPE_TAG(GeneralizedDirichletSpatialParams);
+
+// Set the spatial parameters
+SET_TYPE_PROP(GeneralizedDirichletSpatialParams, SpatialParams,
+        Dumux::GeneralizedDirichletSpatialParams<TypeTag>);
+
+// Set the material law
+SET_PROP(GeneralizedDirichletSpatialParams, MaterialLaw)
+{
+private:
+    // material law typedefs
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    // select material law to be used
+    typedef RegularizedBrooksCorey<Scalar> RawMaterialLaw;
+public:
+    // adapter for absolute law
+    typedef EffToAbsLaw<RawMaterialLaw> type;
+};
+}
+
+/*!
+ * \ingroup TwoPBoxModel
+ *
+ * \brief The spatial parameters for the problem that uses a generalized
+ *        Dirichlet boundary condition.
+ */
+template<class TypeTag>
+class GeneralizedDirichletSpatialParams: public ImplicitSpatialParams<TypeTag>
+{
+    // Get informations for current implementation via property system
+    typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    enum
+    {
+        dim = Grid::dimension,
+        dimWorld = GridView::dimensionworld
+    };
+    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+    typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
+
+public:
+    // get material law from property system
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+    // determine appropriate parameters depending on selected materialLaw
+    typedef typename MaterialLaw::Params MaterialLawParams;
+
+    /*! Intrinsic permeability tensor K \f$[m^2]\f$ depending
+     *  on the position in the domain
+     */
+    const DimWorldMatrix& intrinsicPermeabilityAtPos(const GlobalPosition& globalPos) const
+    { return K_; }
+
+    /*! Defines the porosity \f$[-]\f$ of the porous medium depending
+     * on the position in the domain
+     */
+    Scalar porosityAtPos(const GlobalPosition& globalPos) const
+    { return 0.2; }
+
+    /*! Returns the parameter object for the material law (i.e. Brooks-Corey)
+     *  depending on the position in the domain
+     */
+    const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
+    {
+        return materialParams_;
+    }
+
+    // constructor
+    GeneralizedDirichletSpatialParams(const GridView& gridView) :
+        ImplicitSpatialParams<TypeTag>(gridView),
+        K_(0)
+    {
+        //set main diagonal entries of the permeability tensor to a value
+        //setting to one value means: isotropic, homogeneous
+        for (int i = 0; i < dim; i++)
+            K_[i][i] = 1e-7;
+
+        //set residual saturations
+        materialParams_.setSwr(0.0);
+        materialParams_.setSnr(0.0);
+
+        //parameters of Brooks & Corey Law
+        materialParams_.setPe(500.0);
+        materialParams_.setLambda(2);
+    }
+
+private:
+    DimWorldMatrix K_;
+    // Object that holds the values/parameters of the selected material law.
+    MaterialLawParams materialParams_;
+};
+} // end namespace
+#endif

test/implicit/2p/test_generalizeddirichlet.cc

+// -*- 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/>.   *
+ *****************************************************************************/
+#include "config.h"
+#include "generalizeddirichletproblem.hh"
+#include <dumux/common/start.hh>
+
+void usage(const char *progName, const std::string &errorMsg)
+{
+    std::cout
+        <<  "\nUsage: " << progName << " [options]\n";
+    if (errorMsg.size() > 0)
+        std::cout << errorMsg << "\n";
+    std::cout
+        << "\n"
+        << "The list of mandatory arguments for this program is:\n"
+        << "\t-TEnd                The end of the simulation [s]\n"
+        << "\t-DtInitial           The initial timestep size [s]\n"
+        << "\t-Grid.UpperRightX    The coordinate of the right corner [m]\n"
+        << "\t-Grid.NumberOfCellsX The grid's resolution\n"
+        << "\n";
+}
+
+int main(int argc, char** argv)
+{
+    typedef TTAG(GeneralizedDirichletProblem) TypeTag;
+    return Dumux::start<TypeTag>(argc, argv, usage);
+}

test/implicit/2p/test_generalizeddirichlet.input

+[TimeManager]
+TEnd = 500000 # duration of the simulation [s]
+DtInitial = 10 # initial time step size [s]
+
+[Problem]
+Name = generalizeddirichlet # name passed to the output routines
+
+[Grid]
+UpperRightX = 300 # coordinate of the right corner of the grid [m]
+NumberOfCellsX = 100 # resolution of the grid

test/references/generalizeddirichlet-reference.vtp

+<?xml version="1.0"?>
+<VTKFile type="PolyData" version="0.1" byte_order="LittleEndian">
+  <PolyData>
+    <Piece NumberOfLines="100" NumberOfPoints="101">
+      <PointData Scalars="Sn">
+        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
+          0.21989 0.232707 0.243682 0.253351 0.262038 0.269954 0.277249 0.284031 0.290382 0.296366 0.302033 0.307425
+          0.312575 0.317513 0.322262 0.326844 0.331275 0.335573 0.33975 0.343819 0.347793 0.35168 0.35549 0.359232
+          0.362914 0.366543 0.370127 0.373673 0.377187 0.380677 0.384147 0.387606 0.391059 0.394514 0.397976 0.401453
+          0.404953 0.408483 0.412051 0.415667 0.419339 0.423079 0.426897 0.430806 0.434819 0.438953 0.443224 0.44765
+          0.452254 0.457059 0.462093 0.467387 0.472977 0.478905 0.485219 0.491973 0.499233 0.507075 0.515591 0.524885
+          0.535087 0.546345 0.558842 0.57279 0.588442 0.606094 0.626088 0.648806 0.674668 0.704121 0.737637 0.775705
+          0.818757 0.866848 0.918452 0.966678 0.996098 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1
+        </DataArray>
+        <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
+          0.78011 0.767293 0.756318 0.746649 0.737962 0.730046 0.722751 0.715969 0.709618 0.703634 0.697967 0.692575
+          0.687425 0.682487 0.677738 0.673156 0.668725 0.664427 0.66025 0.656181 0.652207 0.64832 0.64451 0.640768
+          0.637086 0.633457 0.629873 0.626327 0.622813 0.619323 0.615853 0.612394 0.608941 0.605486 0.602024 0.598547
+          0.595047 0.591517 0.587949 0.584333 0.580661 0.576921 0.573103 0.569194 0.565181 0.561047 0.556776 0.55235
+          0.547746 0.542941 0.537907 0.532613 0.527023 0.521095 0.514781 0.508027 0.500767 0.492925 0.484409 0.475115
+          0.464914 0.453655 0.441158 0.42721 0.411558 0.393906 0.373912 0.351194 0.325332 0.295879 0.262363 0.224295
+          0.181243 0.133152 0.0815484 0.0333222 0.00390185 1.90716e-07 1.32837e-26 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
+          201197 201196 201195 201192 201190 201186 201182 201178 201173 201168 201163 201157
+          201151 201145 201138 201131 201124 201116 201109 201101 201093 201084 201075 201067
+          201057 201048 201038 201029 201019 201008 200998 200987 200976 200965 200954 200942
+          200930 200918 200906 200893 200881 200868 200855 200841 200828 200814 200800 200785
+          200770 200756 200740 200725 200709 200693 200677 200660 200643 200626 200609 200591
+          200573 200554 200535 200517 200497 200478 200459 200439 200420 200401 200382 200364
+          200347 200330 200315 200301 200288 200276 200264 200252 200240 200228 200216 200204
+          200192 200180 200168 200156 200144 200132 200120 200108 200096 200084 200072 200060
+          200048 200036 200024 200012 200000
+        </DataArray>
+        <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
+          200631 200625 200620 200614 200608 200601 200594 200587 200580 200572 200564 200556
+          200548 200539 200531 200522 200512 200503 200493 200484 200473 200463 200453 200442
+          200431 200420 200408 200397 200385 200373 200361 200348 200335 200322 200309 200296
+          200282 200268 200254 200239 200225 200210 200194 200178 200162 200146 200129 200112
+          200095 200077 200059 200040 200020 200001 199980 199959 199937 199914 199890 199865
+          199839 199812 199783 199752 199718 199681 199641 199595 199543 199481 199406 199308
+          199172 198960 198564 197562 193764 192776 192764 192752 192740 192728 192716 192704
+          192692 192680 192668 192656 192644 192632 192620 192608 192596 192584 192572 192560
+          192548 192536 192524 192512 192500
+        </DataArray>
+        <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
+          566.099 570.807 574.934 578.645 582.04 585.187 588.133 590.912 593.551 596.069 598.484 600.809
+          603.056 605.233 607.35 609.413 611.429 613.404 615.341 617.246 619.123 620.977 622.809 624.625
+          626.428 628.22 630.004 631.785 633.565 635.347 637.135 638.932 640.741 642.566 644.411 646.28
+          648.178 650.109 652.079 654.093 656.158 658.282 660.471 662.735 665.084 667.529 670.085 672.764
+          675.585 678.568 681.736 685.116 688.74 692.647 696.881 701.498 706.565 712.164 718.396 725.389
+          733.304 742.348 752.788 764.978 779.389 796.662 817.684 843.716 876.61 919.207 976.154 1055.75
+          1174.46 1370.24 1750.9 2739.07 6524.54 7499.95 7500 7500 7500 7500 7500 7500
+          7500 7500 7500 7500 7500 7500 7500 7500 7500 7500 7500 7500
+          7500 7500 7500 7500 7500
+        </DataArray>
+        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
+          999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749
+          999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749
+          999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749
+          999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749
+          999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749
+          999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749 999.749
+          999.749 999.748 999.748 999.748 999.746 999.745 999.745 999.745 999.745 999.745 999.745 999.745
+          999.745 999.745 999.745 999.745 999.745 999.745 999.745 999.745 999.745 999.745 999.745 999.745
+          999.745 999.745 999.745 999.745 999.745
+        </DataArray>
+        <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
+          890 890 890 890 890 890 890 890 890 890 890 890
+          890 890 890 890 890 890 890 890 890 890 890 890
+          890 890 890 890 890 890 890 890 890 890 890 890
+          890 890 890 890 890 890 890 890 890 890 890 890
+          890 890 890 890 890 890 890 890 890 890 890 890
+          890 890 890 890 890 890 890 890 890 890 890 890
+          890 890 890 890 890 890 890 890 890 890 890 890
+          890 890 890 890 890 890 890 890 890 890 890 890
+          890 890 890 890 890
+        </DataArray>
+        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
+          283.624 265.439 250.574 238.004 227.12 217.53 208.965 201.231 194.186 187.719 181.744 176.193
+          171.009 166.149 161.572 157.247 153.147 149.248 145.53 141.975 138.568 135.294 132.141 129.099
+          126.157 123.307 120.54 117.848 115.226 112.665 110.161 107.707 105.298 102.929 100.595 98.2905
+          96.0118 93.7538 91.512 89.2817 87.0583 84.8372 82.6135 80.3825 78.1391 75.8781 73.5941 71.2815
+          68.9347 66.5474 64.1134 61.6263 59.0794 56.4658 53.7788 51.0114 48.1573 45.2107 42.1667 39.0223
+          35.7773 32.4355 29.0066 25.5086 21.9707 18.4369 14.9692 11.6495 8.57882 5.86913 3.62855 1.9382
+          0.826347 0.240717 0.0338673 0.000944182 1.775e-07 1.01312e-24 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
+          2.36578 2.78387 3.17677 3.55047 3.90879 4.2544 4.58927 4.91493 5.23261 5.54332 5.84792 6.14716
+          6.44167 6.73204 7.0188 7.30242 7.58336 7.86202 8.13882 8.41413 8.68832 8.96176 9.23481 9.50783
+          9.78119 10.0552 10.3304 10.607 10.8855 11.1664 11.45 11.7369 12.0276 12.3226 12.6226 12.9283
+          13.2403 13.5595 13.8867 14.223 14.5695 14.9274 15.298 15.6831 16.0843 16.5037 16.9436 17.4067
+          17.896 18.4152 18.9683 19.5602 20.1965 20.8841 21.6308 22.4462 23.3417 24.3313 25.4319 26.6642
+          28.0539 29.6328 31.4404 33.5262 35.9517 38.794 42.1478 46.1288 50.875 56.5479 63.3318 71.4308
+          81.0428 92.263 104.743 116.679 124.025 125 125 125 125 125 125 125
+          125 125 125 125 125 125 125 125 125 125 125 125
+          125 125 125 125 125
+        </DataArray>
+        <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
+          0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
+          0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
+          0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
+          0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
+          0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
+          0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
+          0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
+          0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
+          0.2 0.2 0.2 0.2 0.2
+        </DataArray>
+        <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
+          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
+          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
+          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
+          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
+          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
+          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
+          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
+          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
+          283.15 283.15 283.15 283.15 283.15
+        </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 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 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 0 0 3 0 0 6 0 0 9 0 0
+          12 0 0 15 0 0 18 0 0 21 0 0
+          24 0 0 27 0 0 30 0 0 33 0 0
+          36 0 0 39 0 0 42 0 0 45 0 0
+          48 0 0 51 0 0 54 0 0 57 0 0
+          60 0 0 63 0 0 66 0 0 69 0 0
+          72 0 0 75 0 0 78 0 0 81 0 0
+          84 0 0 87 0 0 90 0 0 93 0 0
+          96 0 0 99 0 0 102 0 0 105 0 0
+          108 0 0 111 0 0 114 0 0 117 0 0
+          120 0 0 123 0 0 126 0 0 129 0 0
+          132 0 0 135 0 0 138 0 0 141 0 0
+          144 0 0 147 0 0 150 0 0 153 0 0
+          156 0 0 159 0 0 162 0 0 165 0 0
+          168 0 0 171 0 0 174 0 0 177 0 0
+          180 0 0 183 0 0 186 0 0 189 0 0
+          192 0 0 195 0 0 198 0 0 201 0 0
+          204 0 0 207 0 0 210 0 0 213 0 0
+          216 0 0 219 0 0 222 0 0 225 0 0
+          228 0 0 231 0 0 234 0 0 237 0 0
+          240 0 0 243 0 0 246 0 0 249 0 0
+          252 0 0 255 0 0 258 0 0 261 0 0
+          264 0 0 267 0 0 270 0 0 273 0 0
+          276 0 0 279 0 0 282 0 0 285 0 0
+          288 0 0 291 0 0 294 0 0 297 0 0
+          300 0 0
+        </DataArray>
+      </Points>
+      <Lines>
+        <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+          0 1 1 2 2 3 3 4 4 5 5 6
+          6 7 7 8 8 9 9 10 10 11 11 12
+          12 13 13 14 14 15 15 16 16 17 17 18
+          18 19 19 20 20 21 21 22 22 23 23 24
+          24 25 25 26 26 27 27 28 28 29 29 30
+          30 31 31 32 32 33 33 34 34 35 35 36
+          36 37 37 38 38 39 39 40 40 41 41 42
+          42 43 43 44 44 45 45 46 46 47 47 48
+          48 49 49 50 50 51 51 52 52 53 53 54
+          54 55 55 56 56 57 57 58 58 59 59 60
+          60 61 61 62 62 63 63 64 64 65 65 66
+          66 67 67 68 68 69 69 70 70 71 71 72
+          72 73 73 74 74 75 75 76 76 77 77 78
+          78 79 79 80 80 81 81 82 82 83 83 84
+          84 85 85 86 86 87 87 88 88 89 89 90
+          90 91 91 92 92 93 93 94 94 95 95 96
+          96 97 97 98 98 99 99 100
+        </DataArray>
+        <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+          2 4 6 8 10 12 14 16 18 20 22 24
+          26 28 30 32 34 36 38 40 42 44 46 48
+          50 52 54 56 58 60 62 64 66 68 70 72
+          74 76 78 80 82 84 86 88 90 92 94 96
+          98 100 102 104 106 108 110 112 114 116 118 120
+          122 124 126 128 130 132 134 136 138 140 142 144
+          146 148 150 152 154 156 158 160 162 164 166 168
+          170 172 174 176 178 180 182 184 186 188 190 192
+          194 196 198 200
+        </DataArray>
+      </Lines>
+    </Piece>
+  </PolyData>
+</VTKFile>