[poroelastic] add test

test/geomechanics/poroelastic/CMakeLists.txt

+dune_symlink_to_source_files(FILES "poroelastic.input")
+
+# using box and numeric differentiation
+dune_add_test(NAME test_poroelastic_box
+              SOURCES test_poroelastic.cc
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS  --script fuzzy
+                        --files ${CMAKE_SOURCE_DIR}/test/references/poroelasticbox-reference.vtu
+                                ${CMAKE_CURRENT_BINARY_DIR}/poroelastic-00001.vtu
+                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_poroelastic_box poroelastic.input")
+
+set(CMAKE_BUILD_TYPE Release)
+
+install(FILES
+problem.hh
+spatialparams.hh
+test_poroelastic.cc
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/test/geomechanics/poroelastic)

test/geomechanics/poroelastic/poroelastic.input

+[Grid]
+UpperRight = 1 1
+Cells = 10 10
+
+[Problem]
+Name = poroelastic
+EnableGravity = false
+
+[Assembly.NumericDifference]
+PriVarMagnitude = 1e5 1e5
+
+[Component]
+SolidDensity = 2700

test/geomechanics/poroelastic/problem.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 Definition of a test problem for the poro-elastic model
+ */
+#ifndef DUMUX_POROELASTIC_PROBLEM_HH
+#define DUMUX_POROELASTIC_PROBLEM_HH
+
+#include <dune/common/fmatrix.hh>
+
+#include <dumux/discretization/box/properties.hh>
+#include <dumux/geomechanics/poroelastic/model.hh>
+#include <dumux/geomechanics/fvproblem.hh>
+
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/material/components/constant.hh>
+
+#include "spatialparams.hh"
+
+namespace Dumux {
+
+template <class TypeTag>
+class PoroElasticProblem;
+
+namespace Properties {
+NEW_TYPE_TAG(PoroElasticTypeTag, INHERITS_FROM(BoxModel, PoroElastic));
+// Set the grid type
+SET_TYPE_PROP(PoroElasticTypeTag, Grid, Dune::YaspGrid<2>);
+// Set the problem property
+SET_TYPE_PROP(PoroElasticTypeTag, Problem, Dumux::PoroElasticProblem<TypeTag>);
+// The fluid phase consists of one constant component
+SET_TYPE_PROP(PoroElasticTypeTag,
+              FluidSystem,
+              Dumux::FluidSystems::OnePLiquid< typename GET_PROP_TYPE(TypeTag, Scalar),
+                                               Dumux::Components::Constant<0, typename GET_PROP_TYPE(TypeTag, Scalar)> >);
+// The spatial parameters property
+SET_TYPE_PROP(PoroElasticTypeTag, SpatialParams, PoroElasticSpatialParams< typename GET_PROP_TYPE(TypeTag, Scalar),
+                                                                           typename GET_PROP_TYPE(TypeTag, FVGridGeometry) >);
+}
+
+/*!
+ * \ingroup Geomechanics
+ * \ingroup PoroElastic
+ *
+ * \brief Problem definition for the deformation of a poro-elastic body
+ */
+template<class TypeTag>
+class PoroElasticProblem : public GeomechanicsFVProblem<TypeTag>
+{
+    using ParentType = GeomechanicsFVProblem<TypeTag>;
+
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using FVElementGeometry = typename FVGridGeometry::LocalView;
+    using SubControlVolume = typename FVGridGeometry::SubControlVolume;
+    using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
+
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+    static constexpr int dim = GridView::dimension;
+    static constexpr int dimWorld = GridView::dimensionworld;
+    using GradU = Dune::FieldMatrix<Scalar, dim, dimWorld>;
+
+public:
+    //! The constructor
+    PoroElasticProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry) {}
+
+    //! The temperature in the domain
+    static constexpr Scalar temperature() { return 273.15; }
+    //! Evaluate the initial value for a control volume.
+    PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const { return PrimaryVariables(0.0); }
+    //! Evaluate the boundary conditions for a Dirichlet boundary segment.
+    PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const { return PrimaryVariables(0.0); }
+    //! Evaluate the boundary conditions for a Neumannboundary segment.
+    PrimaryVariables neumannAtPos(const GlobalPosition& globalPos) const { return PrimaryVariables(0.0); }
+
+    /*!
+     * \brief Returns the effective fluid density
+     * \param globalPos The global position
+     */
+    Scalar effectiveFluidDensityAtPos(const GlobalPosition& globalPos) const
+    {
+        // This test uses the constant component, obtain density only once
+        using FS = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+        static const Scalar rho = FS::density( effectivePorePressureAtPos(globalPos), temperature() );
+        return rho;
+    }
+
+    /*!
+     * \brief Returns the effective pore pressure
+     * \note We use the x-displacement as pressure solution. The shift to
+     *       higher values is done to see a mor pronounced effect in stresses.
+     *
+     * \param globalPos The global position
+     */
+    Scalar effectivePorePressureAtPos(const GlobalPosition& globalPos) const
+    { return exactSolution(globalPos)[0] + 10; }
+
+    /*!
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary segment.
+     * \param globalPos The global position
+     */
+    BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
+    {
+        BoundaryTypes values;
+        values.setAllDirichlet();
+        return values;
+    }
+
+    /*!
+     * \brief Evaluate the source term for all phases within a given
+     *        sub-control-volume.
+     */
+    PrimaryVariables source(const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const SubControlVolume& scv) const
+    {
+        using std::sin;
+        using std::cos;
+
+        static const Scalar pi = 3.14159265358979323846;
+        const auto ipGlobal = scv.center();
+        const auto x = ipGlobal[0];
+        const auto y = ipGlobal[1];
+
+        // the lame parameters (we know they only depend on position here)
+        const auto& lameParams = this->spatialParams().lameParamsAtPos(scv.center());
+        const auto lambda = lameParams.lambda();
+        const auto mu = lameParams.mu();
+
+        // precalculated products
+        const Scalar pi_2 = 2.0*pi;
+        const Scalar pi_2_square = pi_2*pi_2;
+        const Scalar cos_2pix = cos(pi_2*x);
+        const Scalar sin_2pix = sin(pi_2*x);
+        const Scalar cos_2piy = cos(pi_2*y);
+        const Scalar sin_2piy = sin(pi_2*y);
+
+        const Scalar dE11_dx = -2.0*sin_2piy;
+        const Scalar dE22_dx = pi_2_square*cos_2pix*cos_2piy;
+        const Scalar dE11_dy = pi_2*(1.0-2.0*x)*cos_2piy;
+        const Scalar dE22_dy = -1.0*pi_2_square*sin_2pix*sin_2piy;
+        const Scalar dE12_dy = 0.5*pi_2_square*(cos_2pix*cos_2piy - (x-x*x)*sin_2piy);
+        const Scalar dE21_dx = 0.5*((1.0-2*x)*pi_2*cos_2piy - pi_2_square*sin_2pix*sin_2piy);
+
+        // compute exact divergence of sigma
+        PrimaryVariables divSigma(0.0);
+        divSigma[Indices::momentum(/*x-dir*/0)] = lambda*(dE11_dx + dE22_dx) + 2*mu*(dE11_dx + dE12_dy);
+        divSigma[Indices::momentum(/*y-dir*/1)] = lambda*(dE11_dy + dE22_dy) + 2*mu*(dE21_dx + dE22_dy);
+        return divSigma;
+    }
+
+    /*!
+     * \brief Evaluate the exact displacement to this problem at a given position.
+     */
+    PrimaryVariables exactSolution(const GlobalPosition& globalPos) const
+    {
+        using std::sin;
+
+        static const Scalar pi = 3.14159265358979323846;
+        const auto x = globalPos[0];
+        const auto y = globalPos[1];
+
+        PrimaryVariables exact(0.0);
+        exact[Indices::momentum(/*x-dir*/0)] = (x-x*x)*sin(2*pi*y);
+        exact[Indices::momentum(/*y-dir*/1)] = sin(2*pi*x)*sin(2*pi*y);
+        return exact;
+    }
+
+    /*!
+     * \brief Evaluate the exact displacement gradient to this problem at a given position.
+     */
+    GradU exactGradient(const GlobalPosition& globalPos) const
+    {
+        using std::sin;
+        using std::cos;
+
+        static const Scalar pi = 3.14159265358979323846;
+        const auto x = globalPos[0];
+        const auto y = globalPos[1];
+
+        static constexpr int xIdx = Indices::momentum(/*x-dir*/0);
+        static constexpr int yIdx = Indices::momentum(/*y-dir*/1);
+
+        GradU exactGrad(0.0);
+        exactGrad[xIdx][xIdx] = (1-2*x)*sin(2*pi*y);
+        exactGrad[xIdx][yIdx] = (x - x*x)*2*pi*cos(2*pi*y);
+        exactGrad[yIdx][xIdx] = 2*pi*cos(2*pi*x)*sin(2*pi*y);
+        exactGrad[yIdx][yIdx] = 2*pi*sin(2*pi*x)*cos(2*pi*y);
+        return exactGrad;
+    }
+
+private:
+    static constexpr Scalar eps_ = 3e-6;
+};
+
+} //end namespace
+
+#endif

test/geomechanics/poroelastic/spatialparams.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 Definition of the spatial parameters for the poro-elastic problem
+ */
+#ifndef DUMUX_POROELASTIC_SPATIAL_PARAMS_HH
+#define DUMUX_POROELASTIC_SPATIAL_PARAMS_HH
+
+#include <dumux/geomechanics/lameparams.hh>
+#include <dumux/material/spatialparams/fvporoelastic.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Geomechanics
+ * \ingroup PoroElastic
+ * \brief Definition of the spatial parameters for the poro-elastic problem
+ */
+template<class Scalar, class FVGridGeometry>
+class PoroElasticSpatialParams : public FVSpatialParamsPoroElastic< Scalar,
+                                                                    FVGridGeometry,
+                                                                    PoroElasticSpatialParams<Scalar, FVGridGeometry> >
+{
+    using ThisType = PoroElasticSpatialParams<Scalar, FVGridGeometry>;
+    using ParentType = FVSpatialParamsPoroElastic<Scalar, FVGridGeometry, ThisType>;
+
+    using SubControlVolume = typename FVGridGeometry::SubControlVolume;
+    using GridView = typename FVGridGeometry::GridView;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+    //! export the type of the lame parameters
+    using LameParams = Dumux::LameParams<Scalar>;
+
+    //! The constructor
+    PoroElasticSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry)
+    {
+        lameParams_.setLambda(2);
+        lameParams_.setMu(2);
+    }
+
+    //! Define the Lame parameters
+    const LameParams& lameParamsAtPos(const GlobalPosition& globalPos) const { return lameParams_; }
+    //! Return the porosity of the porous medium
+    Scalar porosityAtPos(const GlobalPosition& globalPos) const { return 0.3; }
+    //! Return the biot coefficient of the porous medium
+    Scalar biotCoefficientAtPos(const GlobalPosition& globalPos) const { return 1.0; }
+
+private:
+    LameParams lameParams_;
+};
+} // end namespace Dumux
+#endif

test/geomechanics/poroelastic/test_poroelastic.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/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief test for the one-phase CC model
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/common/fvector.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include "problem.hh"
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/elementsolution.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+
+// function to evaluate the element stresses
+template< class StressType,
+          class Problem,
+          class GridVariables,
+          class SolutionVector,
+          class SigmaStorage >
+void assembleElementStresses(SigmaStorage& sigmaStorage,
+                             SigmaStorage& effSigmaStorage,
+                             const Problem& problem,
+                             const typename GridVariables::GridGeometry& fvGridGeometry,
+                             const GridVariables& gridVariables,
+                             const SolutionVector& x)
+{
+    for (const auto& element : elements(fvGridGeometry.gridView()))
+    {
+        auto fvGeometry = localView(fvGridGeometry);
+        auto elemVolVars = localView(gridVariables.curGridVolVars());
+
+        fvGeometry.bind(element);
+        elemVolVars.bind(element, fvGeometry, x);
+
+        // evaluate flux variables cache at cell center
+        using FluxVarsCache = typename GridVariables::GridFluxVariablesCache::FluxVariablesCache;
+        FluxVarsCache fluxVarsCache;
+        fluxVarsCache.update(problem, element, fvGeometry, elemVolVars, element.geometry().center());
+
+        // get lame parameters, the pressure and compute stress tensor
+        const auto sigma = StressType::stressTensor(problem, element, fvGeometry, elemVolVars, fluxVarsCache);
+        const auto effSigma = StressType::effectiveStressTensor(problem, element, fvGeometry, elemVolVars, fluxVarsCache);
+
+        // pass values into storage container
+        using FVGridGeometry = typename GridVariables::GridGeometry;
+        for (int dir = 0; dir < FVGridGeometry::GridView::dimension; ++dir)
+        {
+            const auto eIdx = fvGridGeometry.elementMapper().index(element);
+            sigmaStorage[dir][eIdx] = sigma[dir];
+            effSigmaStorage[dir][eIdx] = effSigma[dir];
+        }
+    }
+}
+
+// main function
+int main(int argc, char** argv) try
+{
+    using namespace Dumux;
+
+    // define the type tag for this problem
+    using TypeTag = TTAG(PoroElasticTypeTag);
+
+    // stop time for the entire computation
+    Dune::Timer timer;
+
+    // initialize MPI, finalize is done automatically on exit
+    const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+    // print dumux start message
+    if (mpiHelper.rank() == 0)
+        DumuxMessage::print(/*firstCall=*/true);
+
+    // parse command line arguments and input file
+    Parameters::init(argc, argv);
+
+    // try to create a grid (from the given grid file or the input file)
+    using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+    GridCreator::makeGrid();
+    GridCreator::loadBalance();
+
+    ////////////////////////////////////////////////////////////
+    // run non-linear problem on this grid
+    ////////////////////////////////////////////////////////////
+
+    // we compute on the leaf grid view
+    const auto& leafGridView = GridCreator::grid().leafGridView();
+
+    // create the finite volume grid geometry
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+    fvGridGeometry->update();
+
+    // the problem (initial and boundary conditions)
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+    // the solution vector
+    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+    SolutionVector x(fvGridGeometry->numDofs());
+    problem->applyInitialSolution(x);
+
+    // the grid variables
+    using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+    auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+    gridVariables->init(x);
+
+    // intialize the vtk output module and output fields
+    using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+    using VtkOutputModule = VtkOutputModule<TypeTag>;
+    VtkOutputModule vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+    VtkOutputFields::init(vtkWriter);
+
+    // also, add displacement and exact solution to the output
+    vtkWriter.addField(x, "u");
+    SolutionVector xExact(fvGridGeometry->numDofs());
+    for (const auto& v : vertices(leafGridView))
+        xExact[ fvGridGeometry->vertexMapper().index(v) ] = problem->exactSolution(v.geometry().center());
+    vtkWriter.addField(xExact, "u_exact");
+
+    // Furthermore, write out element stress tensors
+    static constexpr int dim = FVGridGeometry::GridView::dimension;
+    static constexpr int dimWorld = FVGridGeometry::GridView::dimensionworld;
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using ForceVector = Dune::FieldVector< Scalar, dimWorld >;
+
+    // containers to store sigma/effective sigma
+    std::array< std::vector<ForceVector>, dim > sigmaStorage;
+    std::array< std::vector<ForceVector>, dim > effSigmaStorage;
+
+    const auto numCells = fvGridGeometry->gridView().size(0);
+    std::for_each(sigmaStorage.begin(), sigmaStorage.end(), [numCells] (auto& sigma) { sigma.resize(numCells); });
+    std::for_each(effSigmaStorage.begin(), effSigmaStorage.end(), [numCells] (auto& effSigma) { effSigma.resize(numCells); });
+
+    for (int dir = 0; dir < dim; ++dir)
+    {
+        vtkWriter.addField(sigmaStorage[dir], "sigma_" + std::to_string(dir), VtkOutputModule::FieldType::element);
+        vtkWriter.addField(effSigmaStorage[dir], "effSigma_" + std::to_string(dir), VtkOutputModule::FieldType::element);
+    }
+
+    // use convenience function to compute stresses
+    using StressType = GET_PROP_TYPE(TypeTag, StressType);
+    assembleElementStresses<StressType>(sigmaStorage, effSigmaStorage, *problem, *fvGridGeometry, *gridVariables, x);
+
+    // write initial solution
+    vtkWriter.write(0.0);
+
+    // the assembler with time loop for instationary problem
+    using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+    auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables);
+
+    // the linear solver
+    using LinearSolver = AMGBackend<TypeTag>;
+    auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
+
+    // the non-linear solver
+    using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+    NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+    // linearize & solve
+    nonLinearSolver.solve(x);
+
+    // the grid variables need to be up to date for subsequent output
+    gridVariables->update(x);
+
+    // write vtk output
+    assembleElementStresses<StressType>(sigmaStorage, effSigmaStorage, *problem, *fvGridGeometry, *gridVariables, x);
+    vtkWriter.write(1.0);
+
+    // print time and say goodbye
+    const auto& comm = Dune::MPIHelper::getCollectiveCommunication();
+    if (mpiHelper.rank() == 0)
+        std::cout << "Simulation took " << timer.elapsed() << " seconds on "
+                  << comm.size() << " processes.\n"
+                  << "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
+
+    // print parameters
+    if (mpiHelper.rank() == 0)
+        Parameters::print();
+
+    return 0;
+} // end main
+catch (Dumux::ParameterException& e)
+{
+    std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+    return 1;
+}
+catch (Dune::DGFException& e)
+{
+    std::cerr << "DGF exception thrown (" << e <<
+                 "). Most likely, the DGF file name is wrong "
+                 "or the DGF file is corrupted, "
+                 "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+                 << " ---> Abort!" << std::endl;
+    return 2;
+}
+catch (Dune::Exception& e)
+{
+    std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+    return 3;
+}
+catch (...)
+{
+    std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+    return 4;
+}