From af8308dd35106a17bafc5692918ad2053946fed5 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Mon, 18 Dec 2017 17:11:15 +0100
Subject: [PATCH] [navierstokes][test] Adapt test_kovasznay to changes

* make stationary
---
 test/freeflow/staggered/CMakeLists.txt        |  2 +-
 .../staggered/kovasznaytestproblem.hh         | 43 ++++------
 test/freeflow/staggered/test_kovasznay.cc     | 80 ++++---------------
 .../test_navierstokes_kovasznay.input         |  4 -
 4 files changed, 31 insertions(+), 98 deletions(-)

diff --git a/test/freeflow/staggered/CMakeLists.txt b/test/freeflow/staggered/CMakeLists.txt
index 37f17e8202..e59bd2a8c4 100644
--- a/test/freeflow/staggered/CMakeLists.txt
+++ b/test/freeflow/staggered/CMakeLists.txt
@@ -77,7 +77,7 @@ dune_add_test(NAME test_navierstokes_kovasznay
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS       --script fuzzy
                              --files ${CMAKE_SOURCE_DIR}/test/references/test_kovasznay-reference.vtu
-                                     ${CMAKE_CURRENT_BINARY_DIR}/test_kovasznay-00002.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_kovasznay-00001.vtu
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_navierstokes_kovasznay")
 
 dune_add_test(NAME test_navierstokes_angeli
diff --git a/test/freeflow/staggered/kovasznaytestproblem.hh b/test/freeflow/staggered/kovasznaytestproblem.hh
index 49415b0fa1..3b04ff4342 100644
--- a/test/freeflow/staggered/kovasznaytestproblem.hh
+++ b/test/freeflow/staggered/kovasznaytestproblem.hh
@@ -32,22 +32,11 @@
 #include <dumux/discretization/staggered/freeflow/properties.hh>
 #include <dumux/freeflow/navierstokes/model.hh>
 
-// solve Navier-Stokes equations
-#define ENABLE_NAVIERSTOKES 1
-
-
 namespace Dumux
 {
 template <class TypeTag>
 class KovasznayTestProblem;
 
-namespace Capabilities
-{
-    template<class TypeTag>
-    struct isStationary<KovasznayTestProblem<TypeTag>>
-    { static const bool value = true; };
-}
-
 namespace Properties
 {
 NEW_TYPE_TAG(KovasznayTestProblem, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokes));
@@ -70,11 +59,7 @@ SET_BOOL_PROP(KovasznayTestProblem, EnableFVGridGeometryCache, true);
 SET_BOOL_PROP(KovasznayTestProblem, EnableGridFluxVariablesCache, true);
 SET_BOOL_PROP(KovasznayTestProblem, EnableGridVolumeVariablesCache, true);
 
-#if ENABLE_NAVIERSTOKES
 SET_BOOL_PROP(KovasznayTestProblem, EnableInertiaTerms, true);
-#else
-SET_BOOL_PROP(KovasznayTestProblem, EnableInertiaTerms, false);
-#endif
 }
 
 /*!
@@ -299,10 +284,10 @@ public:
                 // 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();
+                const auto analyticalSolutionCellCenter = dirichletAtPos(posCellCenter)[pressureIdx];
+                const auto numericalSolutionCellCenter = curSol[cellCenterIdx][dofIdxCellCenter][pressureIdx];
+                sumError[pressureIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * scv.volume();
+                sumReference[pressureIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * scv.volume();
                 totalVolume += scv.volume();
 
                 // treat face dofs
@@ -310,7 +295,7 @@ public:
                 {
                     const int dofIdxFace = scvf.dofIndex();
                     const int dirIdx = scvf.directionIndex();
-                    const auto analyticalSolutionFace = dirichletAtPos(scvf.center())[faceIdx][dirIdx];
+                    const auto analyticalSolutionFace = dirichletAtPos(scvf.center())[Indices::velocity(dirIdx)];
                     const auto numericalSolutionFace = curSol[faceIdx][dofIdxFace][momentumBalanceIdx];
                     directionIndex[dofIdxFace] = dirIdx;
                     errorVelocity[dofIdxFace] = squaredDiff_(analyticalSolutionFace, numericalSolutionFace);
@@ -322,8 +307,8 @@ public:
         }
 
         // 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]);
+        l2NormAbs[pressureIdx] = std::sqrt(sumError[pressureIdx] / totalVolume);
+        l2NormRel[pressureIdx] = std::sqrt(sumError[pressureIdx] / sumReference[pressureIdx]);
 
         // get the absolute and relative discrete L2-error for face dofs
         for(int i = 0; i < numFaceDofs; ++i)
@@ -332,14 +317,14 @@ public:
             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;
+            sumError[Indices::velocity(dirIdx)] += error * volume;
+            sumReference[Indices::velocity(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]);
+            l2NormAbs[Indices::velocity(dirIdx)] = std::sqrt(sumError[Indices::velocity(dirIdx)] / totalVolume);
+            l2NormRel[Indices::velocity(dirIdx)] = std::sqrt(sumError[Indices::velocity(dirIdx)] / sumReference[Indices::velocity(dirIdx)]);
         }
         return std::make_pair(l2NormAbs, l2NormRel);
     }
@@ -396,11 +381,13 @@ private:
                     const auto faceDofPosition = scvf.center();
                     const auto dirIdx = scvf.directionIndex();
                     const auto analyticalSolutionAtFace = analyticalSolution(faceDofPosition);
-                    analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[faceIdx][dirIdx];
+                    analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
                 }
 
                 analyticalPressure_[ccDofIdx] = analyticalSolutionAtCc[pressureIdx];
-                analyticalVelocity_[ccDofIdx] = analyticalSolutionAtCc[faceIdx];
+
+                for(int dirIdx = 0; dirIdx < dim; ++dirIdx)
+                    analyticalVelocity_[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)];
             }
         }
     }
diff --git a/test/freeflow/staggered/test_kovasznay.cc b/test/freeflow/staggered/test_kovasznay.cc
index b1c855a136..50afa4e0a8 100644
--- a/test/freeflow/staggered/test_kovasznay.cc
+++ b/test/freeflow/staggered/test_kovasznay.cc
@@ -19,7 +19,7 @@
 /*!
  * \file
  *
- * \brief Test for the staggered grid Stokes model (Donea et al., 2003)
+ * \brief Stationary test for the staggered grid Navier-Stokes model (Kovasznay 1947)
  */
  #include <config.h>
 
@@ -119,7 +119,7 @@ int main(int argc, char** argv) try
     auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
     fvGridGeometry->update();
 
-    // the problem (initial and boundary conditions)
+    // the problem (boundary conditions)
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     auto problem = std::make_shared<Problem>(fvGridGeometry);
 
@@ -133,25 +133,11 @@ int main(int argc, char** argv) try
     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);
-
-    // 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");
+    gridVariables->init(x);
 
     // intialize the vtk output module
     using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
@@ -162,13 +148,9 @@ int main(int argc, char** argv) try
     vtkWriter.addFaceField(problem->getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
     vtkWriter.write(0.0);
 
-    // instantiate time loop
-    auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
-    timeLoop->setMaxTimeStepSize(maxDt);
-
     // the assembler with time loop for instationary problem
     using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
-    auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+    auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables);
 
     // the linear solver
     using LinearSolver = Dumux::UMFPackBackend<TypeTag>;
@@ -177,55 +159,23 @@ int main(int argc, char** argv) try
     // the non-linear solver
     using NewtonController = StaggeredNewtonController<TypeTag>;
     using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
-    auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
+    auto newtonController = std::make_shared<NewtonController>(leafGridView.comm());
     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();
-
-        // advance to the time loop to the next step
-        timeLoop->advanceTimeStep();
-
-        problem->postTimeStep(x);
-
-        // write vtk output
-        vtkWriter.write(timeLoop->time());
+    // linearize & solve
+    Dune::Timer timer;
+    nonLinearSolver.solve(x);
 
-        // report statistics of this time step
-        timeLoop->reportTimeStep();
+    // write vtk output
+    vtkWriter.write(1.0);
 
-        // set new dt as suggested by newton controller
-        timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+    timer.stop();
 
-    } while (!timeLoop->finished());
+    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";
 
-    timeLoop->finalize(leafGridView.comm());
 
     ////////////////////////////////////////////////////////////
     // finalize, print dumux message to say goodbye
diff --git a/test/freeflow/staggered/test_navierstokes_kovasznay.input b/test/freeflow/staggered/test_navierstokes_kovasznay.input
index 1bcc53d179..03029c6f46 100644
--- a/test/freeflow/staggered/test_navierstokes_kovasznay.input
+++ b/test/freeflow/staggered/test_navierstokes_kovasznay.input
@@ -1,7 +1,3 @@
-[TimeLoop]
-DtInitial = 1 # [s]
-TEnd = 2 # [s]
-
 [Grid]
 LowerLeft = -0.5 -0.5
 UpperRight = 2 1.5
-- 
GitLab