[NavierStokes][test] Add transient test

test/freeflow/staggered/CMakeLists.txt

@@ -79,3 +79,12 @@ dune_add_test(NAME test_navierstokes_kovasznay
                              --files ${CMAKE_SOURCE_DIR}/test/references/test_kovasznay-reference.vtu
                                      ${CMAKE_CURRENT_BINARY_DIR}/test_kovasznay-00002.vtu
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_navierstokes_kovasznay")
+
+dune_add_test(NAME test_navierstokes_angeli
+              SOURCES test_angeli.cc
+              CMAKE_GUARD HAVE_UMFPACK
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS       --script fuzzy
+                             --files ${CMAKE_SOURCE_DIR}/test/references/test_angeli-reference.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_angeli-00045.vtu
+                             --command "${CMAKE_CURRENT_BINARY_DIR}/test_navierstokes_angeli")

test/freeflow/staggered/angelitestproblem.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 Test for the instationary staggered grid Navier-Stokes model with analytical solution (Angeli et al., 2017)
+ */
+#ifndef DUMUX_ANGELI_TEST_PROBLEM_HH
+#define DUMUX_ANGELI_TEST_PROBLEM_HH
+
+#include <dumux/freeflow/staggered/problem.hh>
+#include <dumux/discretization/staggered/properties.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
+#include <dumux/material/components/constant.hh>
+
+#include <dumux/discretization/staggered/properties.hh>
+#include <dumux/freeflow/staggered/properties.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class AngeliTestProblem;
+
+namespace Capabilities
+{
+    template<class TypeTag>
+    struct isStationary<AngeliTestProblem<TypeTag>>
+    { static const bool value = false; };
+}
+
+namespace Properties
+{
+NEW_TYPE_TAG(AngeliTestProblem, INHERITS_FROM(StaggeredModel, NavierStokes));
+
+SET_PROP(AngeliTestProblem, Fluid)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+public:
+    typedef FluidSystems::LiquidPhase<Scalar, Dumux::Components::Constant<1, Scalar> > type;
+};
+
+// Set the grid type
+SET_TYPE_PROP(AngeliTestProblem, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+
+// Set the problem property
+SET_TYPE_PROP(AngeliTestProblem, Problem, Dumux::AngeliTestProblem<TypeTag> );
+
+SET_BOOL_PROP(AngeliTestProblem, EnableFVGridGeometryCache, true);
+
+SET_BOOL_PROP(AngeliTestProblem, EnableGlobalFluxVariablesCache, true);
+SET_BOOL_PROP(AngeliTestProblem, EnableGlobalVolumeVariablesCache, true);
+
+SET_BOOL_PROP(AngeliTestProblem, EnableInertiaTerms, true);
+}
+
+/*!
+ * \ingroup ImplicitTestProblems
+ * \brief  Test problem for the staggered grid (Angeli 1947)
+ * \todo doc me!
+ */
+template <class TypeTag>
+class AngeliTestProblem : public NavierStokesProblem<TypeTag>
+{
+    using ParentType = NavierStokesProblem<TypeTag>;
+
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+
+    // copy some indices for convenience
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+    enum {
+        // Grid and world dimension
+        dim = GridView::dimension,
+        dimWorld = GridView::dimensionworld
+    };
+    enum {
+        massBalanceIdx = Indices::massBalanceIdx,
+        momentumBalanceIdx = Indices::momentumBalanceIdx,
+        momentumXBalanceIdx = Indices::momentumXBalanceIdx,
+        momentumYBalanceIdx = Indices::momentumYBalanceIdx,
+        pressureIdx = Indices::pressureIdx,
+        velocityXIdx = Indices::velocityXIdx,
+        velocityYIdx = Indices::velocityYIdx
+    };
+
+    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+    using Element = typename GridView::template Codim<0>::Entity;
+    using Intersection = typename GridView::Intersection;
+
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+
+    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+
+    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+    using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+    using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+    using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+    typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+    typename DofTypeIndices::FaceIdx faceIdx;
+
+public:
+    AngeliTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry), eps_(1e-6)
+    {
+        printL2Error_ = getParam<bool>("Problem.PrintL2Error");
+
+        kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
+
+        using CellArray = std::array<unsigned int, dimWorld>;
+        const auto numCells = getParam<CellArray>("Grid.Cells");
+
+        cellSizeX_ = this->fvGridGeometry().bBoxMax()[0] / numCells[0];
+    }
+
+    /*!
+     * \name Problem parameters
+     */
+    // \{
+
+    bool shouldWriteRestartFile() const
+    {
+        return false;
+    }
+
+    void postTimeStep(const SolutionVector& curSol) const
+    {
+        if(printL2Error_)
+        {
+            const auto l2error = calculateL2Error(curSol);
+            const int numCellCenterDofs = this->fvGridGeometry().gridView().size(0);
+            const int numFaceDofs = this->fvGridGeometry().gridView().size(1);
+            std::cout << std::setprecision(8) << "** L2 error (abs/rel) for "
+                    << std::setw(6) << numCellCenterDofs << " cc dofs and " << numFaceDofs << " face dofs (total: " << numCellCenterDofs + numFaceDofs << "): "
+                    << std::scientific
+                    << "L2(p) = " << l2error.first[pressureIdx] << " / " << l2error.second[pressureIdx]
+                    << ", L2(vx) = " << l2error.first[velocityXIdx] << " / " << l2error.second[velocityXIdx]
+                    << ", L2(vy) = " << l2error.first[velocityYIdx] << " / " << l2error.second[velocityYIdx]
+                    << std::endl;
+        }
+    }
+
+    /*!
+     * \brief Return the temperature within the domain in [K].
+     *
+     * This problem assumes a temperature of 10 degrees Celsius.
+     */
+    Scalar temperature() const
+    { return 298.0; }
+
+
+    /*!
+     * \brief Return the sources within the domain.
+     *
+     * \param globalPos The global position
+     */
+    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
+    {
+        return SourceValues(0.0);
+    }
+
+    // \}
+    /*!
+     * \name Boundary conditions
+     */
+    // \{
+
+    /*!
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary control volume.
+     *
+     * \param globalPos The position of the center of the finite volume
+     */
+    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+    {
+        BoundaryTypes values;
+
+        // set Dirichlet values for the velocity everywhere
+        values.setDirichlet(momentumBalanceIdx);
+
+        // set a fixed pressure in one cell
+        values.setDirichletCell(massBalanceIdx);
+
+        return values;
+    }
+
+    /*!
+     * \brief Return dirichlet boundary values at a given position
+     *
+     * \param globalPos The global position
+     */
+    PrimaryVariables dirichletAtPos(const GlobalPosition & globalPos) const
+    {
+        // use the values of the analytical solution
+        return analyticalSolution(globalPos, time());
+    }
+
+     /*!
+     * \brief Return the analytical solution of the problem at a given position
+     *
+     * \param globalPos The global position
+     */
+    PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, const Scalar time) const
+    {
+        Scalar x = globalPos[0];
+        Scalar y = globalPos[1];
+
+        Scalar t = time + timeLoop_->timeStepSize();
+
+        PrimaryVariables values;
+
+        values[pressureIdx] = - 0.25 * std::exp(-10.0 * kinematicViscosity_ * M_PI * M_PI * t) * M_PI * M_PI * (4.0 * std::cos(2.0 * M_PI * x) + std::cos(4.0 * M_PI * y));
+        values[velocityXIdx] = - 2.0 * M_PI * std::exp(- 5.0 * kinematicViscosity_ * M_PI * M_PI * t) * std::cos(M_PI * x) * std::sin(2.0 * M_PI * y);
+        values[velocityYIdx] = M_PI * std::exp(- 5.0 * kinematicViscosity_ * M_PI * M_PI * t) * std::sin(M_PI * x) * std::cos(2.0 * M_PI * y);
+
+        return values;
+    }
+
+    // \}
+
+    /*!
+     * \name Volume terms
+     */
+    // \{
+
+    /*!
+     * \brief Evaluate the initial value for a control volume.
+     *
+     * \param globalPos The global position
+     */
+    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+    {
+        return analyticalSolution(globalPos, -timeLoop_->timeStepSize());
+    }
+
+
+    /*!
+     * \brief Calculate the L2 error between the analytical solution and the numerical approximation.
+     *
+     */
+    auto calculateL2Error(const SolutionVector& curSol) const
+    {
+        PrimaryVariables sumError(0.0), sumReference(0.0), l2NormAbs(0.0), l2NormRel(0.0);
+
+        const int numFaceDofs = this->fvGridGeometry().gridView().size(1);
+
+        std::vector<Scalar> staggeredVolume(numFaceDofs);
+        std::vector<Scalar> errorVelocity(numFaceDofs);
+        std::vector<Scalar> velocityReference(numFaceDofs);
+        std::vector<int> directionIndex(numFaceDofs);
+
+        Scalar totalVolume = 0.0;
+
+        for (const auto& element : elements(this->fvGridGeometry().gridView()))
+        {
+            auto fvGeometry = localView(this->fvGridGeometry());
+            fvGeometry.bindElement(element);
+
+            for (auto&& scv : scvs(fvGeometry))
+            {
+                // treat cell-center dofs
+                const auto dofIdxCellCenter = scv.dofIndex();
+                const auto& posCellCenter = scv.dofPosition();
+                const auto analyticalSolutionCellCenter = dirichletAtPos(posCellCenter)[cellCenterIdx];
+                const auto numericalSolutionCellCenter = curSol[cellCenterIdx][dofIdxCellCenter];
+                sumError[cellCenterIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * scv.volume();
+                sumReference[cellCenterIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * scv.volume();
+                totalVolume += scv.volume();
+
+                // treat face dofs
+                for (auto&& scvf : scvfs(fvGeometry))
+                {
+                    const int dofIdxFace = scvf.dofIndex();
+                    const int dirIdx = scvf.directionIndex();
+                    const auto analyticalSolutionFace = dirichletAtPos(scvf.center())[faceIdx][dirIdx];
+                    const auto numericalSolutionFace = curSol[faceIdx][dofIdxFace][momentumBalanceIdx];
+                    directionIndex[dofIdxFace] = dirIdx;
+                    errorVelocity[dofIdxFace] = squaredDiff_(analyticalSolutionFace, numericalSolutionFace);
+                    velocityReference[dofIdxFace] = squaredDiff_(analyticalSolutionFace, 0.0);
+                    const Scalar staggeredHalfVolume = 0.5 * scv.volume();
+                    staggeredVolume[dofIdxFace] = staggeredVolume[dofIdxFace] + staggeredHalfVolume;
+                }
+            }
+        }
+
+        // get the absolute and relative discrete L2-error for cell-center dofs
+        l2NormAbs[cellCenterIdx] = std::sqrt(sumError[cellCenterIdx] / totalVolume);
+        l2NormRel[cellCenterIdx] = std::sqrt(sumError[cellCenterIdx] / sumReference[cellCenterIdx]);
+
+        // get the absolute and relative discrete L2-error for face dofs
+        for(int i = 0; i < numFaceDofs; ++i)
+        {
+            const int dirIdx = directionIndex[i];
+            const auto error = errorVelocity[i];
+            const auto ref = velocityReference[i];
+            const auto volume = staggeredVolume[i];
+            sumError[faceIdx][dirIdx] += error * volume;
+            sumReference[faceIdx][dirIdx] += ref * volume;
+        }
+
+        for(int dirIdx = 0; dirIdx < dimWorld; ++dirIdx)
+        {
+            l2NormAbs[faceIdx][dirIdx] = std::sqrt(sumError[faceIdx][dirIdx] / totalVolume);
+            l2NormRel[faceIdx][dirIdx] = std::sqrt(sumError[faceIdx][dirIdx] / sumReference[faceIdx][dirIdx]);
+        }
+        return std::make_pair(l2NormAbs, l2NormRel);
+    }
+
+    /*!
+     * \brief Returns the analytical solution for the pressure
+     */
+    auto& getAnalyticalPressureSolution() const
+    {
+        return analyticalPressure_;
+    }
+
+    /*!
+     * \brief Returns the analytical solution for the velocity
+     */
+    auto& getAnalyticalVelocitySolution() const
+    {
+        return analyticalVelocity_;
+    }
+
+    /*!
+     * \brief Returns the analytical solution for the velocity at the faces
+     */
+    auto& getAnalyticalVelocitySolutionOnFace() const
+    {
+        return analyticalVelocityOnFace_;
+    }
+
+    void setTimeLoop(TimeLoopPtr timeLoop)
+    {
+        timeLoop_ = timeLoop;
+        createAnalyticalSolution();
+    }
+
+    Scalar time() const
+    {
+        return timeLoop_->time();
+    }
+
+    /*!
+     * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
+     */
+    void createAnalyticalSolution()
+    {
+        analyticalPressure_.resize(this->fvGridGeometry().numCellCenterDofs());
+        analyticalVelocity_.resize(this->fvGridGeometry().numCellCenterDofs());
+        analyticalVelocityOnFace_.resize(this->fvGridGeometry().numFaceDofs());
+
+        for (const auto& element : elements(this->fvGridGeometry().gridView()))
+        {
+            auto fvGeometry = localView(this->fvGridGeometry());
+            fvGeometry.bindElement(element);
+            for (auto&& scv : scvs(fvGeometry))
+            {
+                auto ccDofIdx = scv.dofIndex();
+                auto ccDofPosition = scv.dofPosition();
+                auto analyticalSolutionAtCc = analyticalSolution(ccDofPosition, time());
+
+                // velocities on faces
+                for (auto&& scvf : scvfs(fvGeometry))
+                {
+                    const auto faceDofIdx = scvf.dofIndex();
+                    const auto faceDofPosition = scvf.center();
+                    const auto dirIdx = scvf.directionIndex();
+                    const auto analyticalSolutionAtFace = analyticalSolution(faceDofPosition, time());
+                    analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[faceIdx][dirIdx];
+                }
+
+                analyticalPressure_[ccDofIdx] = analyticalSolutionAtCc[pressureIdx];
+                analyticalVelocity_[ccDofIdx] = analyticalSolutionAtCc[faceIdx];
+            }
+        }
+    }
+
+    private:
+
+    template<class T>
+    T squaredDiff_(const T& a, const T& b) const
+    {
+        return (a-b)*(a-b);
+    }
+
+    bool isLowerLeftCell_(const GlobalPosition& globalPos) const
+    {
+        return globalPos[0] < (this->fvGridGeometry().bBoxMin()[0] + 0.5*cellSizeX_ + eps_);
+    }
+
+    Scalar eps_;
+    Scalar cellSizeX_;
+
+    Scalar kinematicViscosity_;
+    TimeLoopPtr timeLoop_;
+    bool printL2Error_;
+    std::vector<Scalar> analyticalPressure_;
+    std::vector<GlobalPosition> analyticalVelocity_;
+    std::vector<GlobalPosition> analyticalVelocityOnFace_;
+};
+} //end namespace
+
+#endif

test/freeflow/staggered/test_angeli.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 instationary staggered grid Navier-Stokes model with analytical solution (Angeli et al., 2017)
+ */
+ #include <config.h>
+
+ #include <ctime>
+ #include <iostream>
+
+ #include <dune/common/parallel/mpihelper.hh>
+ #include <dune/common/timer.hh>
+ #include <dune/grid/io/file/dgfparser/dgfexception.hh>
+ #include <dune/grid/io/file/vtk.hh>
+ #include <dune/istl/io.hh>
+
+
+#include "angelitestproblem.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/seqsolverbackend.hh>
+#include <dumux/nonlinear/newtonmethod.hh>
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+#include <dumux/assembly/staggeredfvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/implicit/staggered/newtoncontroller.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ *        reading in parameters.
+ *
+ * \param progName  The name of the program, that was tried to be started.
+ * \param errorMsg  The error message that was issued by the start function.
+ *                  Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+    if (errorMsg.size() > 0) {
+        std::string errorMessageOut = "\nUsage: ";
+                    errorMessageOut += progName;
+                    errorMessageOut += " [options]\n";
+                    errorMessageOut += errorMsg;
+                    errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n"
+                                        "\t-TimeManager.TEnd               End of the simulation [s] \n"
+                                        "\t-TimeManager.DtInitial          Initial timestep size [s] \n"
+                                        "\t-Grid.File                      Name of the file containing the grid \n"
+                                        "\t                                definition in DGF format\n"
+                                        "\t-SpatialParams.LensLowerLeftX   x-coordinate of the lower left corner of the lens [m] \n"
+                                        "\t-SpatialParams.LensLowerLeftY   y-coordinate of the lower left corner of the lens [m] \n"
+                                        "\t-SpatialParams.LensUpperRightX  x-coordinate of the upper right corner of the lens [m] \n"
+                                        "\t-SpatialParams.LensUpperRightY  y-coordinate of the upper right corner of the lens [m] \n"
+                                        "\t-SpatialParams.Permeability     Permeability of the domain [m^2] \n"
+                                        "\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
+
+        std::cout << errorMessageOut
+                  << "\n";
+    }
+}
+
+int main(int argc, char** argv) try
+{
+    using namespace Dumux;
+
+    // define the type tag for this problem
+    using TypeTag = TTAG(AngeliTestProblem);
+
+    // 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, usage);
+
+    // 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 instationary 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();
+
+    // get some time loop parameters
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+    const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
+    const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+     auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+    // check if we are about to restart a previously interrupted simulation
+    Scalar restartTime = 0;
+    if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+        restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+    // instantiate time loop
+    auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+    timeLoop->setMaxTimeStepSize(maxDt);
+
+    // the problem (initial and boundary conditions)
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    auto problem = std::make_shared<Problem>(fvGridGeometry);
+    problem->setTimeLoop(timeLoop);
+
+    // the solution vector
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+    using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+    typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+    typename DofTypeIndices::FaceIdx faceIdx;
+    const auto numDofsCellCenter = leafGridView.size(0);
+    const auto numDofsFace = leafGridView.size(1);
+    SolutionVector x;
+    x[cellCenterIdx].resize(numDofsCellCenter);
+    x[faceIdx].resize(numDofsFace);
+    problem->applyInitialSolution(x);
+    auto xOld = x;
+
+    // the grid variables
+    using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+    auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+    gridVariables->init(x, xOld);
+
+    // intialize the vtk output module
+    using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+    StaggeredVtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+    VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+    vtkWriter.addField(problem->getAnalyticalPressureSolution(), "pressureExact", 1);
+    vtkWriter.addField(problem->getAnalyticalVelocitySolution(), "velocityExact", GridView::dimensionworld);
+    vtkWriter.addFaceField(problem->getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
+    vtkWriter.write(0.0);
+
+    // the assembler with time loop for instationary problem
+    using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
+    auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+    // the linear solver
+    using LinearSolver = Dumux::UMFPackBackend<TypeTag>;
+    auto linearSolver = std::make_shared<LinearSolver>();
+
+    // the non-linear solver
+    using NewtonController = StaggeredNewtonController<TypeTag>;
+    using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
+    auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
+    NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
+
+    // time loop
+    timeLoop->start(); do
+    {
+        // set previous solution for storage evaluations
+        assembler->setPreviousSolution(xOld);
+
+        // try solving the non-linear system
+        for (int i = 0; i < maxDivisions; ++i)
+        {
+            // linearize & solve
+            auto converged = nonLinearSolver.solve(x);
+
+            if (converged)
+                break;
+
+            if (!converged && i == maxDivisions-1)
+                DUNE_THROW(Dune::MathError,
+                           "Newton solver didn't converge after "
+                           << maxDivisions
+                           << " time-step divisions. dt="
+                           << timeLoop->timeStepSize()
+                           << ".\nThe solutions of the current and the previous time steps "
+                           << "have been saved to restart files.");
+        }
+
+        // make the new solution the old solution
+        xOld = x;
+        gridVariables->advanceTimeStep();
+
+        problem->createAnalyticalSolution();
+
+        // advance to the time loop to the next step
+        timeLoop->advanceTimeStep();
+
+        problem->postTimeStep(x);
+
+        // write vtk output
+        vtkWriter.write(timeLoop->time());
+
+        // report statistics of this time step
+        timeLoop->reportTimeStep();
+
+        // set new dt as suggested by newton controller
+        timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+
+    } while (!timeLoop->finished());
+
+    timeLoop->finalize(leafGridView.comm());
+
+    ////////////////////////////////////////////////////////////
+    // finalize, print dumux message to say goodbye
+    ////////////////////////////////////////////////////////////
+
+    // print dumux end message
+    if (mpiHelper.rank() == 0)
+    {
+        Parameters::print();
+        DumuxMessage::print(/*firstCall=*/false);
+    }
+
+    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;
+}

test/freeflow/staggered/test_navierstokes_angeli.input

+[TimeLoop]
+DtInitial = 1e-12 # [s]
+TEnd = 1e-2 # [s]
+
+[Grid]
+LowerLeft = 0 0
+UpperRight = 1 1
+Cells = 50 50
+
+[Problem]
+Name = test_angeli # name passed to the output routines
+EnableGravity = false
+PrintL2Error = false
+
+[Component]
+LiquidDensity = 1
+LiquidKinematicViscosity = 0.1
+
+[ Newton ]
+MaxSteps = 10
+MaxRelativeShift = 1e-5
+
+[Vtk]
+WriteFaceData = false
+AddVelocity = true