[freeflow][test] Add convergence sincos test. Unite steady and unsteady sincos...

[freeflow][test] Add convergence sincos test. Unite steady and unsteady sincos test in one main file, probelem file and executable.

test/freeflow/navierstokes/sincos/CMakeLists.txt

-add_subdirectory(steady)
-add_subdirectory(unsteady)
+dune_symlink_to_source_files(FILES "convergencetest.py")
+
+add_executable(test_ff_navierstokes_sincos EXCLUDE_FROM_ALL main.cc)
+
+dune_add_test(NAME test_ff_navierstokes_sincos_unsteady
+              TARGET test_ff_navierstokes_sincos
+              LABELS freeflow
+              CMAKE_GUARD HAVE_UMFPACK
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS       --script fuzzy
+                             --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_navierstokes_sincos_unsteady-reference.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_sincos_unsteady-00017.vtu
+                             --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_sincos params_unsteady.input
+                             -Problem.Name test_ff_navierstokes_sincos_unsteady")
+
+dune_add_test(NAME test_ff_navierstokes_sincos_steady
+              TARGET test_ff_navierstokes_sincos
+              LABELS freeflow
+              CMAKE_GUARD HAVE_UMFPACK
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS       --script fuzzy
+                             --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_navierstokes_sincos_steady-reference.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_sincos_steady-00001.vtu
+                             --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_sincos params_steady.input
+                             -Problem.Name test_ff_navierstokes_sincos_steady")
+
+dune_add_test(NAME test_ff_navierstokes_sincos_steady_convergence
+              TARGET test_ff_navierstokes_sincos
+              LABELS freeflow
+              CMAKE_GUARD HAVE_UMFPACK
+              COMMAND ./convergencetest.py
+              CMD_ARGS       test_ff_navierstokes_sincos params_steady.input  -Grid.Cells "150 150"
+                             -Problem.Name test_ff_navierstokes_sincos_steady_convergence)
+
+dune_symlink_to_source_files(FILES "params_steady.input")
+dune_symlink_to_source_files(FILES "params_unsteady.input")

test/freeflow/navierstokes/sincos/convergencetest.py

+#!/usr/bin/env python
+
+from math import *
+import subprocess
+import sys
+import numpy
+
+if len(sys.argv) < 2:
+    sys.stderr.write('Please provide a single argument <testname> to the script\n')
+    sys.exit(1)
+
+testname = str(sys.argv[1])
+testargs = [str(i) for i in sys.argv][2:]
+
+# remove the old log file
+subprocess.call(['rm', testname + '.log'])
+print("Removed old log file ({})!".format(testname + '.log'))
+
+# do the runs with different refinement
+for i in [0, 1]:
+    subprocess.call(['./' + testname] + testargs + ['-Grid.Refinement', str(i),
+                                      '-Problem.Name', testname])
+
+# check the rates and append them to the log file
+logfile = open(testname + '.log', "r+")
+
+errorP = []
+errorVx = []
+errorVy = []
+for line in logfile:
+    line = line.strip("\n")
+    line = line.strip("\[ConvergenceTest\]")
+    line = line.split()
+    errorP.append(float(line[2]))
+    errorVx.append(float(line[5]))
+    errorVy.append(float(line[8]))
+
+resultsP = []
+resultsVx = []
+resultsVy = []
+logfile.truncate(0)
+logfile.write("n\terrorP\t\trateP\t\terrorVx\t\trateVx\t\terrorVy\t\trateVy\n")
+logfile.write("-"*50 + "\n")
+for i in range(len(errorP)-1):
+    if isnan(errorP[i]) or isinf(errorP[i]):
+        continue
+    if not ((errorP[i] < 1e-12 or errorP[i+1] < 1e-12) and (errorVx[i] < 1e-12 or errorVx[i+1] < 1e-12) and (errorVy[i] < 1e-12 or errorVy[i+1] < 1e-12)):
+        rateP = (log(errorP[i])-log(errorP[i+1]))/log(2)
+        rateVx = (log(errorVx[i])-log(errorVx[i+1]))/log(2)
+        rateVy = (log(errorVy[i])-log(errorVy[i+1]))/log(2)
+        message = "{}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\n".format(i, errorP[i], rateP,  errorVx[i], rateVx, errorVy[i], rateVy)
+        logfile.write(message)
+        resultsP.append(rateP)
+        resultsVx.append(rateVx)
+        resultsVy.append(rateVy)
+    else:
+        logfile.write("error: exact solution!?")
+i = len(errorP)-1
+message = "{}\t{:0.4e}\t\t{}\t{:0.4e}\t\t{}\t{:0.4e}\t\t{}\n".format(i, errorP[i], "",  errorVx[i], "", errorVy[i], "")
+logfile.write(message)
+
+logfile.close()
+print("\nComputed the following convergence rates for {}:\n".format(testname))
+subprocess.call(['cat', testname + '.log'])
+
+# check the rates, we expect rates around 2
+if int(numpy.mean(resultsP) < 2.05 and numpy.mean(resultsP) < 1.95):
+    sys.stderr.write("*"*70 + "\n" + "The convergence rates for pressure were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
+    sys.exit(1)
+
+if int(numpy.mean(resultsVx) < 2.05 and numpy.mean(resultsVx) < 1.95):
+    sys.stderr.write("*"*70 + "\n" + "The convergence rates for x-velocity were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
+    sys.exit(1)
+
+if int(numpy.mean(resultsVy) < 2.05 and numpy.mean(resultsVy) < 1.95):
+    sys.stderr.write("*"*70 + "\n" + "The convergence rates for y-velocity were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
+    sys.exit(1)
+
+sys.exit(0)

test/freeflow/navierstokes/sincos/steady/main.cc → test/freeflow/navierstokes/sincos/main.cc

@@ -22,6 +22,7 @@
  * \brief Test for the instationary staggered grid Navier-Stokes model
  *        with analytical solution.
  */
+
 #include <config.h>
 
 #include <ctime>
@@ -50,10 +51,11 @@
 /*!
 * \brief Creates analytical solution.
 * Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces
+* \param time the time at which to evaluate the analytical solution
 * \param problem the problem for which to evaluate the analytical solution
 */
-template<class Problem>
-auto createAnalyticalSolution(const Problem& problem)
+template<class Scalar, class Problem>
+auto createAnalyticalSolution(const Scalar time, const Problem& problem)
 {
     const auto& fvGridGeometry = problem.fvGridGeometry();
     using GridView = typename std::decay_t<decltype(fvGridGeometry)>::GridView;
@@ -61,7 +63,6 @@ auto createAnalyticalSolution(const Problem& problem)
     static constexpr auto dim = GridView::dimension;
     static constexpr auto dimWorld = GridView::dimensionworld;
 
-    using Scalar = double;
     using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
 
     std::vector<Scalar> analyticalPressure;
@@ -81,7 +82,7 @@ auto createAnalyticalSolution(const Problem& problem)
         {
             auto ccDofIdx = scv.dofIndex();
             auto ccDofPosition = scv.dofPosition();
-            auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition);
+            auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition, time);
 
             // velocities on faces
             for (auto&& scvf : scvfs(fvGeometry))
@@ -89,7 +90,7 @@ auto createAnalyticalSolution(const Problem& problem)
                 const auto faceDofIdx = scvf.dofIndex();
                 const auto faceDofPosition = scvf.center();
                 const auto dirIdx = scvf.directionIndex();
-                const auto analyticalSolutionAtFace = problem.analyticalSolution(faceDofPosition);
+                const auto analyticalSolutionAtFace = problem.analyticalSolution(faceDofPosition, time);
                 analyticalVelocityOnFace[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
             }
 
@@ -136,12 +137,41 @@ auto createSource(const Problem& problem)
     return source;
 }
 
+template<class Problem, class SolutionVector, class FVGridGeometry>
+void printL2Error(const Problem& problem, const SolutionVector& x, const FVGridGeometry& fvGridGeometry)
+{
+    using namespace Dumux;
+    using Scalar = double;
+    using TypeTag = Properties::TTag::SincosTest;
+    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+
+    using L2Error = NavierStokesTestL2Error<Scalar, ModelTraits, PrimaryVariables>;
+    const auto l2error = L2Error::calculateL2Error(*problem, x);
+    const int numCellCenterDofs = fvGridGeometry->numCellCenterDofs();
+    const int numFaceDofs = fvGridGeometry->numFaceDofs();
+    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[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx]
+                << " , L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx]
+                << " , L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx]
+                << std::endl;
+
+    // write the norm into a log file
+    std::ofstream logFile;
+    logFile.open(problem->name() + ".log", std::ios::app);
+    logFile << "[ConvergenceTest] L2(p) = " << l2error.first[Indices::pressureIdx] << " L2(vx) = " << l2error.first[Indices::velocityXIdx] << " L2(vy) = " << l2error.first[Indices::velocityYIdx] << std::endl;
+    logFile.close();
+}
+
 int main(int argc, char** argv) try
 {
     using namespace Dumux;
 
     // define the type tag for this problem
-    using TypeTag = Properties::TTag::SincosSteadyTest;
+    using TypeTag = Properties::TTag::SincosTest;
 
     // initialize MPI, finalize is done automatically on exit
     const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
@@ -169,15 +199,28 @@ int main(int argc, char** argv) try
     auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
     fvGridGeometry->update();
 
+    // get some time loop parameters
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+    const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+    auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+    // instantiate time loop
+    auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
+    timeLoop->setMaxTimeStepSize(maxDt);
+
     // the problem (initial and boundary conditions)
     using Problem = GetPropType<TypeTag, Properties::Problem>;
     auto problem = std::make_shared<Problem>(fvGridGeometry);
+    problem->updateTimeStepSize(timeLoop->timeStepSize());
 
     // the solution vector
     using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
     SolutionVector x;
     x[FVGridGeometry::cellCenterIdx()].resize(fvGridGeometry->numCellCenterDofs());
     x[FVGridGeometry::faceIdx()].resize(fvGridGeometry->numFaceDofs());
+    problem->applyInitialSolution(x);
+    auto xOld = x;
 
     // the grid variables
     using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
@@ -195,16 +238,17 @@ int main(int argc, char** argv) try
     vtkWriter.addField(sourceX, "sourceX");
     vtkWriter.addField(sourceY, "sourceY");
 
-    auto analyticalSolution = createAnalyticalSolution(*problem);
+    auto analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
     vtkWriter.addField(std::get<0>(analyticalSolution), "pressureExact");
     vtkWriter.addField(std::get<1>(analyticalSolution), "velocityExact");
     vtkWriter.addFaceField(std::get<2>(analyticalSolution), "faceVelocityExact");
-
     vtkWriter.write(0.0);
 
+    const bool isStationary = getParam<bool>("Problem.IsStationary");
+
     // the assembler with time loop for instationary problem
     using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
-    auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables);
+    auto assembler = isStationary ? std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables) : std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
 
     // the linear solver
     using LinearSolver = Dumux::UMFPackBackend;
@@ -214,42 +258,69 @@ int main(int argc, char** argv) try
     using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
-    const bool printL2Error = getParam<bool>("Problem.PrintL2Error");
+    const bool shouldPrintL2Error = getParam<bool>("Problem.PrintL2Error");
 
-    // linearize & solve
-    Dune::Timer timer;
-    nonLinearSolver.solve(x);
-
-    if (printL2Error)
+    if (isStationary)
     {
-        using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
-        using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
-        using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
-
-        using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-        using L2Error = NavierStokesTestL2Error<Scalar, ModelTraits, PrimaryVariables>;
-        const auto l2error = L2Error::calculateL2Error(*problem, x);
-        const int numCellCenterDofs = fvGridGeometry->numCellCenterDofs();
-        const int numFaceDofs = fvGridGeometry->numFaceDofs();
-        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[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx]
-                  << " , L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx]
-                  << " , L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx]
-                  << std::endl;
+        // linearize & solve
+        Dune::Timer timer;
+        nonLinearSolver.solve(x);
+
+        if (shouldPrintL2Error)
+        {
+            printL2Error(problem, x, fvGridGeometry);
+        }
+
+        // write vtk output
+        analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
+        vtkWriter.write(1.0);
+
+        timer.stop();
+
+        const auto& comm = Dune::MPIHelper::getCollectiveCommunication();
+        std::cout << "Simulation took " << timer.elapsed() << " seconds on "
+                << comm.size() << " processes.\n"
+                << "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
     }
+    else
+    {
+        // time loop
+        timeLoop->start(); do
+        {
+            // set previous solution for storage evaluations
+            assembler->setPreviousSolution(xOld);
 
-    // write vtk output
-    analyticalSolution = createAnalyticalSolution(*problem);
-    vtkWriter.write(1.0);
+            // solve the non-linear system with time step control
+            nonLinearSolver.solve(x, *timeLoop);
 
-    timer.stop();
+            // make the new solution the old solution
+            xOld = x;
+            gridVariables->advanceTimeStep();
 
-    const auto& comm = Dune::MPIHelper::getCollectiveCommunication();
-    std::cout << "Simulation took " << timer.elapsed() << " seconds on "
-              << comm.size() << " processes.\n"
-              << "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
+            if (shouldPrintL2Error)
+            {
+                printL2Error(problem, x, fvGridGeometry);
+            }
+
+            // advance to the time loop to the next step
+            timeLoop->advanceTimeStep();
+            problem->updateTime(timeLoop->time());
+            analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
+
+            // write vtk output
+            vtkWriter.write(timeLoop->time());
+
+            // report statistics of this time step
+            timeLoop->reportTimeStep();
+
+            // set new dt as suggested by newton solver
+            timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+            problem->updateTimeStepSize(timeLoop->timeStepSize());
+
+        } while (!timeLoop->finished());
+
+        timeLoop->finalize(leafGridView.comm());
+    }
 
     ////////////////////////////////////////////////////////////
     // finalize, print dumux message to say goodbye

test/freeflow/navierstokes/sincos/steady/params.input → test/freeflow/navierstokes/sincos/params_steady.input

+[TimeLoop]
+DtInitial = 0.01 # [s]
+TEnd = 1.0 # [s]
+
 [Grid]
 UpperRight = 6.28 6.28
 Cells = 50 50
@@ -7,6 +11,7 @@ Name = test_sincos_steady # name passed to the output routines
 EnableGravity = false
 PrintL2Error = true
 EnableInertiaTerms = true
+IsStationary = true
 
 [Component]
 LiquidDensity = 1.0

test/freeflow/navierstokes/sincos/unsteady/params.input → test/freeflow/navierstokes/sincos/params_unsteady.input

@@ -12,6 +12,7 @@ Name = test_sincos_unsteady # name passed to the output routines
 EnableGravity = false
 PrintL2Error = true
 EnableInertiaTerms = true
+IsStationary = false
 
 [Component]
 LiquidDensity = 1.0

test/freeflow/navierstokes/sincos/unsteady/problem.hh → test/freeflow/navierstokes/sincos/problem.hh

@@ -19,11 +19,11 @@
 /*!
  * \file
  * \ingroup NavierStokesTests
- * \brief Test for the instationary staggered grid Navier-Stokes model
+ * \brief Test for the staggered grid Navier-Stokes model
  *        with analytical solution.
  */
-#ifndef DUMUX_SINCOS_UNSTEADY_TEST_PROBLEM_HH
-#define DUMUX_SINCOS_UNSTEADY_TEST_PROBLEM_HH
+#ifndef DUMUX_SINCOS_STEADY_TEST_PROBLEM_HH
+#define DUMUX_SINCOS_STEADY_TEST_PROBLEM_HH
 
 #include <dune/grid/yaspgrid.hh>
 
@@ -33,21 +33,21 @@
 #include <dumux/freeflow/navierstokes/problem.hh>
 #include <dumux/discretization/staggered/freeflow/properties.hh>
 #include <dumux/freeflow/navierstokes/model.hh>
-#include "../../l2error.hh"
+#include "../l2error.hh"
 
 namespace Dumux {
 template <class TypeTag>
-class SincosUnsteadyTestProblem;
+class SincosTestProblem;
 
 namespace Properties {
 // Create new type tags
 namespace TTag {
-struct SincosUnsteadyTest { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
+struct SincosTest { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
 } // end namespace TTag
 
 // the fluid system
 template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::SincosUnsteadyTest>
+struct FluidSystem<TypeTag, TTag::SincosTest>
 {
 private:
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -57,37 +57,36 @@ public:
 
 // Set the grid type
 template<class TypeTag>
-struct Grid<TypeTag, TTag::SincosUnsteadyTest> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+struct Grid<TypeTag, TTag::SincosTest> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
 
 // Set the problem property
 template<class TypeTag>
-struct Problem<TypeTag, TTag::SincosUnsteadyTest> { using type = Dumux::SincosUnsteadyTestProblem<TypeTag> ; };
+struct Problem<TypeTag, TTag::SincosTest> { using type = Dumux::SincosTestProblem<TypeTag> ; };
 
 template<class TypeTag>
-struct EnableFVGridGeometryCache<TypeTag, TTag::SincosUnsteadyTest> { static constexpr bool value = true; };
+struct EnableFVGridGeometryCache<TypeTag, TTag::SincosTest> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridFluxVariablesCache<TypeTag, TTag::SincosUnsteadyTest> { static constexpr bool value = true; };
+struct EnableGridFluxVariablesCache<TypeTag, TTag::SincosTest> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridVolumeVariablesCache<TypeTag, TTag::SincosUnsteadyTest> { static constexpr bool value = true; };
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::SincosTest> { static constexpr bool value = true; };
 } // end namespace Properties
 
 /*!
  * \ingroup NavierStokesTests
  * \brief  Test problem for the staggered grid.
  *
- * The unsteady, 2D, incompressible Navier-Stokes equations for zero gravity and a Newtonian
+ * The 2D, incompressible Navier-Stokes equations for zero gravity and a Newtonian
  * flow is solved and compared to an analytical solution (sums/products of trigonometric functions).
- * The velocities and pressures are periodical in time. The Dirichlet boundary conditions are
- * time-dependent and consistent with the analytical solution.
+ * For the unsteady case, the velocities and pressures are periodical in time. The Dirichlet boundary conditions are
+ * consistent with the analytical solution and in the unsteady case time-dependent.
  */
 template <class TypeTag>
-class SincosUnsteadyTestProblem : public NavierStokesProblem<TypeTag>
+class SincosTestProblem : public NavierStokesProblem<TypeTag>
 {
     using ParentType = NavierStokesProblem<TypeTag>;
 
     using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
     using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
-    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
     using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
     using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -102,11 +101,14 @@ class SincosUnsteadyTestProblem : public NavierStokesProblem<TypeTag>
     using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
 
 public:
+    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
-    SincosUnsteadyTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    SincosTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
     : ParentType(fvGridGeometry), time_(0.0), timeStepSize_(0.0)
     {
+        isStationary_ = getParam<bool>("Problem.IsStationary");
+        enableInertiaTerms_ = getParam<bool>("Problem.EnableInertiaTerms");
         kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
 
         using CellArray = std::array<unsigned int, dimWorld>;
@@ -146,8 +148,22 @@ public:
         using std::cos;
         using std::sin;
 
-        source[Indices::momentumXBalanceIdx] = -2.0 * cos(x) * sin(y) * (cos(2.0 * t) + sin(2.0 * t) * kinematicViscosity_);
-        source[Indices::momentumYBalanceIdx] =  2.0 * sin(x) * cos(y) * (cos(2.0 * t) + sin(2.0 * t) * kinematicViscosity_);
+        if (isStationary_)
+        {
+            source[Indices::momentumXBalanceIdx] = -2.0 * kinematicViscosity_ * cos(x) * sin(y);
+            source[Indices::momentumYBalanceIdx] =  2.0 * kinematicViscosity_ * cos(y) * sin(x);
+
+            if (!enableInertiaTerms_)
+            {
+                source[Indices::momentumXBalanceIdx] += 0.5 * sin(2.0 * x);
+                source[Indices::momentumYBalanceIdx] += 0.5 * sin(2.0 * y);
+            }
+        }
+        else
+        {
+            source[Indices::momentumXBalanceIdx] = -2.0 * cos(x) * sin(y) * (cos(2.0 * t) + sin(2.0 * t) * kinematicViscosity_);
+            source[Indices::momentumYBalanceIdx] =  2.0 * sin(x) * cos(y) * (cos(2.0 * t) + sin(2.0 * t) * kinematicViscosity_);
+        }
 
         return source;
     }
@@ -221,9 +237,16 @@ public:
         using std::sin;
         using std::cos;
 
-        values[Indices::pressureIdx] = -0.25 * (cos(2.0 * x) + cos(2.0 * y)) * sin(2.0 * t) * sin(2.0 * t);
-        values[Indices::velocityXIdx] = -1.0 * cos(x) * sin(y) * sin(2.0 * t);
-        values[Indices::velocityYIdx] = sin(x) * cos(y) * sin(2.0 * t);
+        values[Indices::pressureIdx] = -0.25 * (cos(2.0 * x) + cos(2.0 * y));
+        values[Indices::velocityXIdx] = -1.0 * cos(x) * sin(y);
+        values[Indices::velocityYIdx] = sin(x) * cos(y);
+
+        if (!isStationary_)
+        {
+            values[Indices::pressureIdx] *= sin(2.0 * t) * sin(2.0 * t);
+            values[Indices::velocityXIdx] *= sin(2.0 * t);
+            values[Indices::velocityYIdx] *= sin(2.0 * t);
+        }
 
         return values;
     }
@@ -242,7 +265,19 @@ public:
      */
     PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
     {
-        return analyticalSolution(globalPos, 0.0);
+        if (isStationary_)
+        {
+            PrimaryVariables values;
+            values[Indices::pressureIdx] = 0.0;
+            values[Indices::velocityXIdx] = 0.0;
+            values[Indices::velocityYIdx] = 0.0;
+
+            return values;
+        }
+        else
+        {
+            return analyticalSolution(globalPos, 0.0);
+        }
     }
 
     /*!
@@ -268,10 +303,13 @@ private:
     Scalar cellSizeY_;
 
     Scalar kinematicViscosity_;
+    bool enableInertiaTerms_;
     Scalar time_;
     Scalar timeStepSize_;
+
+    bool isStationary_;
 };
 
 } // end namespace Dumux
 
-#endif // DUMUX_SINCOS_UNSTEADY_TEST_PROBLEM_HH
+#endif // DUMUX_SINCOS_TEST_PROBLEM_HH

test/freeflow/navierstokes/sincos/steady/CMakeLists.txt

-dune_add_test(NAME test_ff_navierstokes_sincos_steady
-              SOURCES main.cc
-              LABELS freeflow
-              CMAKE_GUARD HAVE_UMFPACK
-              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
-              CMD_ARGS       --script fuzzy
-                             --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_navierstokes_sincos_steady-reference.vtu
-                                     ${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_sincos_steady-00001.vtu
-                             --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_sincos_steady params.input
-                             -Problem.Name test_ff_navierstokes_sincos_steady")
-
-dune_symlink_to_source_files(FILES "params.input")

test/freeflow/navierstokes/sincos/steady/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 3 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
- * \ingroup NavierStokesTests
- * \brief Test for the stationary staggered grid Navier-Stokes model
- *        with analytical solution.
- */
-#ifndef DUMUX_SINCOS_STEADY_TEST_PROBLEM_HH
-#define DUMUX_SINCOS_STEADY_TEST_PROBLEM_HH
-
-#include <dune/grid/yaspgrid.hh>
-
-#include <dumux/material/fluidsystems/1pliquid.hh>
-#include <dumux/material/components/constant.hh>
-
-#include <dumux/freeflow/navierstokes/problem.hh>
-#include <dumux/discretization/staggered/freeflow/properties.hh>
-#include <dumux/freeflow/navierstokes/model.hh>
-#include "../../l2error.hh"
-
-namespace Dumux {
-template <class TypeTag>
-class SincosSteadyTestProblem;
-
-namespace Properties {
-// Create new type tags
-namespace TTag {
-struct SincosSteadyTest { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
-} // end namespace TTag
-
-// the fluid system
-template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::SincosSteadyTest>
-{
-private:
-    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-public:
-    using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
-};
-
-// Set the grid type
-template<class TypeTag>
-struct Grid<TypeTag, TTag::SincosSteadyTest> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
-
-// Set the problem property
-template<class TypeTag>
-struct Problem<TypeTag, TTag::SincosSteadyTest> { using type = Dumux::SincosSteadyTestProblem<TypeTag> ; };
-
-template<class TypeTag>
-struct EnableFVGridGeometryCache<TypeTag, TTag::SincosSteadyTest> { static constexpr bool value = true; };
-template<class TypeTag>
-struct EnableGridFluxVariablesCache<TypeTag, TTag::SincosSteadyTest> { static constexpr bool value = true; };
-template<class TypeTag>
-struct EnableGridVolumeVariablesCache<TypeTag, TTag::SincosSteadyTest> { static constexpr bool value = true; };
-} // end namespace Properties
-
-/*!
- * \ingroup NavierStokesTests
- * \brief  Test problem for the staggered grid.
- *
- * The steady, 2D, incompressible Navier-Stokes equations for zero gravity and a Newtonian
- * flow is solved and compared to an analytical solution (sums/products of trigonometric functions).
- * The Dirichlet boundary conditions are consistent with the analytical solution.
- */
-template <class TypeTag>
-class SincosSteadyTestProblem : public NavierStokesProblem<TypeTag>
-{
-    using ParentType = NavierStokesProblem<TypeTag>;
-
-    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
-    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
-    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
-
-    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
-    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
-    using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
-    using SubControlVolume = typename FVGridGeometry::SubControlVolume;
-    using FVElementGeometry = typename FVGridGeometry::LocalView;
-
-    static constexpr auto dimWorld = GetPropType<TypeTag, Properties::GridView>::dimensionworld;
-    using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
-    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
-    using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
-
-public:
-    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
-    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
-
-    SincosSteadyTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
-    : ParentType(fvGridGeometry)
-    {
-        kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
-        enableInertiaTerms_ = getParam<bool>("Problem.EnableInertiaTerms");
-
-        using CellArray = std::array<unsigned int, dimWorld>;
-        CellArray numCells = getParam<CellArray>("Grid.Cells");
-
-        const unsigned int refinement = getParam<unsigned int>("Grid.Refinement", 0);
-        for(unsigned int i = 0; i < refinement; i++)
-        {
-            numCells[0] *= 2;
-            numCells[1] *= 2;
-        }
-
-        cellSizeX_ = (this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0]) / numCells[0];
-        cellSizeY_ = (this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]) / numCells[1];
-    }
-
-   /*!
-     * \brief Returns 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
-     */
-    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
-    {
-        NumEqVector source(0.0);
-        const Scalar x = globalPos[0];
-        const Scalar y = globalPos[1];
-
-        using std::cos;
-        using std::sin;
-
-        source[Indices::momentumXBalanceIdx] = -2.0 * kinematicViscosity_ * cos(x) * sin(y);
-        source[Indices::momentumYBalanceIdx] =  2.0 * kinematicViscosity_ * cos(y) * sin(x);
-
-        if (!enableInertiaTerms_)
-        {
-            source[Indices::momentumXBalanceIdx] += 0.5 * sin(2.0 * x);
-            source[Indices::momentumYBalanceIdx] += 0.5 * sin(2.0 * y);
-        }
-
-        return source;
-    }
-
-    // \}
-   /*!
-     * \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(Indices::velocityXIdx);
-        values.setDirichlet(Indices::velocityYIdx);
-
-        return values;
-    }
-
-    /*!
-     * \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
-     *
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry
-     * \param scv The sub control volume
-     * \param pvIdx The primary variable index in the solution vector
-     */
-    bool isDirichletCell(const Element& element,
-                         const typename FVGridGeometry::LocalView& fvGeometry,
-                         const typename FVGridGeometry::SubControlVolume& scv,
-                         int pvIdx) const
-    {
-        // set fixed pressure in one cell
-        return (scv.dofIndex() == 0) && pvIdx == Indices::pressureIdx;
-    }
-
-   /*!
-     * \brief Returns 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);
-    }
-
-    /*!
-     * \brief Returns the analytical solution of the problem at a given position.
-     *
-     * \param globalPos The global position
-     */
-    PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
-    {
-        const Scalar x = globalPos[0];
-        const Scalar y = globalPos[1];
-
-        PrimaryVariables values;
-
-        using std::sin;
-        using std::cos;
-
-        values[Indices::pressureIdx] = -0.25 * (cos(2.0 * x) + cos(2.0 * y));
-        values[Indices::velocityXIdx] = -1.0 * cos(x) * sin(y);
-        values[Indices::velocityYIdx] = sin(x) * cos(y);
-
-        return values;
-    }
-
-    // \}
-
-   /*!
-     * \name Volume terms
-     */
-    // \{
-
-   /*!
-     * \brief Evaluates the initial value for a control volume.
-     *
-     * \param globalPos The global position
-     */
-    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
-    {
-        PrimaryVariables values;
-        values[Indices::pressureIdx] = 0.0;
-        values[Indices::velocityXIdx] = 0.0;
-        values[Indices::velocityYIdx] = 0.0;
-
-        return values;
-    }
-
-private:
-    static constexpr Scalar eps_ = 1e-6;
-
-    Scalar cellSizeX_;
-    Scalar cellSizeY_;
-
-    Scalar kinematicViscosity_;
-    bool enableInertiaTerms_;
-};
-
-} // end namespace Dumux
-
-#endif // DUMUX_SINCOS_STEADY_TEST_PROBLEM_HH

test/freeflow/navierstokes/sincos/unsteady/CMakeLists.txt

-dune_add_test(NAME test_ff_navierstokes_sincos_unsteady
-              SOURCES main.cc
-              LABELS freeflow
-              CMAKE_GUARD HAVE_UMFPACK
-              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
-              CMD_ARGS       --script fuzzy
-                             --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_navierstokes_sincos_unsteady-reference.vtu
-                                     ${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_sincos_unsteady-00017.vtu
-                             --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_sincos_unsteady params.input
-                             -Problem.Name test_ff_navierstokes_sincos_unsteady")
-
-dune_symlink_to_source_files(FILES "params.input")

test/freeflow/navierstokes/sincos/unsteady/main.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 3 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
- * \ingroup NavierStokesTests
- * \brief Test for the instationary staggered grid Navier-Stokes model
- *        with analytical solution.
- */
-
-#include <config.h>
-
-#include <ctime>
-#include <iostream>
-#include <type_traits>
-#include <tuple>
-
-#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 <dumux/assembly/staggeredfvassembler.hh>
-#include <dumux/assembly/diffmethod.hh>
-#include <dumux/common/dumuxmessage.hh>
-#include <dumux/common/parameters.hh>
-#include <dumux/common/properties.hh>
-#include <dumux/common/valgrind.hh>
-#include <dumux/io/grid/gridmanager.hh>
-#include <dumux/io/staggeredvtkoutputmodule.hh>
-#include <dumux/linear/seqsolverbackend.hh>
-#include <dumux/nonlinear/newtonsolver.hh>
-
-#include "problem.hh"
-
-/*!
-* \brief Creates analytical solution.
-* Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces
-* \param time the time at which to evaluate the analytical solution
-* \param problem the problem for which to evaluate the analytical solution
-*/
-template<class Scalar, class Problem>
-auto createAnalyticalSolution(const Scalar time, const Problem& problem)
-{
-    const auto& fvGridGeometry = problem.fvGridGeometry();
-    using GridView = typename std::decay_t<decltype(fvGridGeometry)>::GridView;
-
-    static constexpr auto dim = GridView::dimension;
-    static constexpr auto dimWorld = GridView::dimensionworld;
-
-    using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
-
-    std::vector<Scalar> analyticalPressure;
-    std::vector<VelocityVector> analyticalVelocity;
-    std::vector<VelocityVector> analyticalVelocityOnFace;
-
-    analyticalPressure.resize(fvGridGeometry.numCellCenterDofs());
-    analyticalVelocity.resize(fvGridGeometry.numCellCenterDofs());
-    analyticalVelocityOnFace.resize(fvGridGeometry.numFaceDofs());
-
-    using Indices = typename Problem::Indices;
-    for (const auto& element : elements(fvGridGeometry.gridView()))
-    {
-        auto fvGeometry = localView(fvGridGeometry);
-        fvGeometry.bindElement(element);
-        for (auto&& scv : scvs(fvGeometry))
-        {
-            auto ccDofIdx = scv.dofIndex();
-            auto ccDofPosition = scv.dofPosition();
-            auto analyticalSolutionAtCc = problem.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 = problem.analyticalSolution(faceDofPosition, time);
-                analyticalVelocityOnFace[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
-            }
-
-            analyticalPressure[ccDofIdx] = analyticalSolutionAtCc[Indices::pressureIdx];
-
-            for(int dirIdx = 0; dirIdx < dim; ++dirIdx)
-                analyticalVelocity[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)];
-        }
-    }
-
-    return std::make_tuple(analyticalPressure, analyticalVelocity, analyticalVelocityOnFace);
-}
-
-int main(int argc, char** argv) try
-{
-    using namespace Dumux;
-
-    // define the type tag for this problem
-    using TypeTag = Properties::TTag::SincosUnsteadyTest;
-
-    // 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)
-    GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
-    gridManager.init();
-
-    ////////////////////////////////////////////////////////////
-    // run instationary non-linear problem on this grid
-    ////////////////////////////////////////////////////////////
-
-    // we compute on the leaf grid view
-    const auto& leafGridView = gridManager.grid().leafGridView();
-
-    // create the finite volume grid geometry
-    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
-    auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
-    fvGridGeometry->update();
-
-    // get some time loop parameters
-    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-    const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
-    const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
-    auto dt = getParam<Scalar>("TimeLoop.DtInitial");
-
-    // instantiate time loop
-    auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
-    timeLoop->setMaxTimeStepSize(maxDt);
-
-    // the problem (initial and boundary conditions)
-    using Problem = GetPropType<TypeTag, Properties::Problem>;
-    auto problem = std::make_shared<Problem>(fvGridGeometry);
-    problem->updateTimeStepSize(timeLoop->timeStepSize());
-
-    // the solution vector
-    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
-    SolutionVector x;
-    x[FVGridGeometry::cellCenterIdx()].resize(fvGridGeometry->numCellCenterDofs());
-    x[FVGridGeometry::faceIdx()].resize(fvGridGeometry->numFaceDofs());
-    problem->applyInitialSolution(x);
-    auto xOld = x;
-
-    // the grid variables
-    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
-    auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
-    gridVariables->init(x);
-
-    // initialize the vtk output module
-    StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
-    using IOFields = GetPropType<TypeTag, Properties::IOFields>;
-    IOFields::initOutputModule(vtkWriter); // Add model specific output fields
-
-    auto analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
-    vtkWriter.addField(std::get<0>(analyticalSolution), "pressureExact");
-    vtkWriter.addField(std::get<1>(analyticalSolution), "velocityExact");
-    vtkWriter.addFaceField(std::get<2>(analyticalSolution), "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;
-    auto linearSolver = std::make_shared<LinearSolver>();
-
-    // the non-linear solver
-    using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
-    NewtonSolver nonLinearSolver(assembler, linearSolver);
-
-    const bool printL2Error = getParam<bool>("Problem.PrintL2Error");
-
-    // time loop
-    timeLoop->start(); do
-    {
-        // set previous solution for storage evaluations
-        assembler->setPreviousSolution(xOld);
-
-        // solve the non-linear system with time step control
-        nonLinearSolver.solve(x, *timeLoop);
-
-        // make the new solution the old solution
-        xOld = x;
-        gridVariables->advanceTimeStep();
-
-        if (printL2Error)
-        {
-            using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
-            using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
-            using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
-
-            using L2Error = NavierStokesTestL2Error<Scalar, ModelTraits, PrimaryVariables>;
-            const auto l2error = L2Error::calculateL2Error(*problem, x);
-            const int numCellCenterDofs = fvGridGeometry->numCellCenterDofs();
-            const int numFaceDofs = fvGridGeometry->numFaceDofs();
-            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[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx]
-                      << " , L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx]
-                      << " , L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx]
-                      << std::endl;
-        }
-
-        // advance to the time loop to the next step
-        timeLoop->advanceTimeStep();
-        problem->updateTime(timeLoop->time());
-        analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
-
-        // write vtk output
-        vtkWriter.write(timeLoop->time());
-
-        // report statistics of this time step
-        timeLoop->reportTimeStep();
-
-        // set new dt as suggested by newton solver
-        timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
-        problem->updateTimeStepSize(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;
-}