Merge branch 'feature/simplify-time-dependency-angeli' into 'master'

Improve angeli test

See merge request !2607

test/freeflow/navierstokes/angeli/CMakeLists.txt

@@ -9,4 +9,18 @@ dumux_add_test(NAME test_ff_navierstokes_angeli
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_angeli params.input
                              -Problem.Name test_ff_navierstokes_angeli")
 
+dumux_add_test(NAME test_ff_navierstokes_angeli_averaged
+              LABELS freeflow navierstokes
+              TARGET test_ff_navierstokes_angeli
+              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_angeli_averaged-reference.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_angeli_averaged-00001.vtu
+                             --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_angeli params.input
+                             -TimeLoop.TEnd 1e-6 -TimeLoop.DtInitial 1e-6
+                             -Problem.Name test_ff_navierstokes_angeli_averaged
+                             -Problem.UseVelocityAveragingForInitial true
+                             -Problem.UseVelocityAveragingForDirichlet true")
+
 dune_symlink_to_source_files(FILES "params.input")

test/freeflow/navierstokes/angeli/main.cc

@@ -49,11 +49,10 @@
 /*!
 * \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)
+template<class Problem, class Scalar = double>
+auto createAnalyticalSolution(const Problem& problem)
 {
     const auto& gridGeometry = problem.gridGeometry();
     using GridView = typename std::decay_t<decltype(gridGeometry)>::GridView;
@@ -80,7 +79,7 @@ auto createAnalyticalSolution(const Scalar time, const Problem& problem)
         {
             auto ccDofIdx = scv.dofIndex();
             auto ccDofPosition = scv.dofPosition();
-            auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition, time);
+            auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition);
 
             // velocities on faces
             for (auto&& scvf : scvfs(fvGeometry))
@@ -88,7 +87,7 @@ auto createAnalyticalSolution(const Scalar time, const Problem& problem)
                 const auto faceDofIdx = scvf.dofIndex();
                 const auto faceDofPosition = scvf.center();
                 const auto dirIdx = scvf.directionIndex();
-                const auto analyticalSolutionAtFace = problem.analyticalSolution(faceDofPosition, time);
+                const auto analyticalSolutionAtFace = problem.analyticalSolution(faceDofPosition);
                 analyticalVelocityOnFace[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
             }
 
@@ -137,18 +136,19 @@ int main(int argc, char** argv)
 
     // get some time loop parameters
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    const auto tStart = getParam<Scalar>("TimeLoop.TStart", 0.0);
     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, dt, tEnd);
+    auto timeLoop = std::make_shared<TimeLoop<Scalar>>(tStart, dt, tEnd);
     timeLoop->setMaxTimeStepSize(maxDt);
 
     // the problem (initial and boundary conditions)
     using Problem = GetPropType<TypeTag, Properties::Problem>;
     auto problem = std::make_shared<Problem>(gridGeometry);
-    problem->updateTimeStepSize(timeLoop->timeStepSize());
+    problem->setTime(timeLoop->time());
 
     // the solution vector
     using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
@@ -168,7 +168,7 @@ int main(int argc, char** argv)
     using IOFields = GetPropType<TypeTag, Properties::IOFields>;
     IOFields::initOutputModule(vtkWriter); // Add model specific output fields
 
-    auto analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
+    auto analyticalSolution = createAnalyticalSolution(*problem);
     vtkWriter.addField(std::get<0>(analyticalSolution), "pressureExact");
     vtkWriter.addField(std::get<1>(analyticalSolution), "velocityExact");
     vtkWriter.addFaceField(std::get<2>(analyticalSolution), "faceVelocityExact");
@@ -190,6 +190,8 @@ int main(int argc, char** argv)
     const bool printL2Error = getParam<bool>("Problem.PrintL2Error");
     timeLoop->start(); do
     {
+        problem->setTime(timeLoop->time() + timeLoop->timeStepSize());
+
         // solve the non-linear system with time step control
         nonLinearSolver.solve(x, *timeLoop);
 
@@ -204,22 +206,23 @@ int main(int argc, char** argv)
             using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
 
             using L2Error = NavierStokesTestL2Error<Scalar, ModelTraits, PrimaryVariables>;
-            const auto l2error = L2Error::calculateL2Error(*problem, x);
+            const auto [l2errorAbs, l2errorRel] = L2Error::calculateL2Error(*problem, x);
             const int numCellCenterDofs = gridGeometry->numCellCenterDofs();
             const int numFaceDofs = gridGeometry->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]
+                    << "L2(p) = " << l2errorAbs[Indices::pressureIdx] << " / " << l2errorRel[Indices::pressureIdx]
+                    << ", L2(vx) = " << l2errorAbs[Indices::velocityXIdx] << " / " << l2errorRel[Indices::velocityXIdx]
+                    << ", L2(vy) = " << l2errorAbs[Indices::velocityYIdx] << " / " << l2errorRel[Indices::velocityYIdx]
                     << std::endl;
         }
 
+        // update the analytical solution
+        analyticalSolution = createAnalyticalSolution(*problem);
+
         // 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());
@@ -229,7 +232,6 @@ int main(int argc, char** argv)
 
         // set new dt as suggested by newton solver
         timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
-        problem->updateTimeStepSize(timeLoop->timeStepSize());
 
     } while (!timeLoop->finished());
 

test/freeflow/navierstokes/angeli/params.input

@@ -12,6 +12,8 @@ Name = test_angeli # name passed to the output routines
 EnableGravity = false
 PrintL2Error = false
 EnableInertiaTerms = true
+UseVelocityAveragingForInitial = false
+UseVelocityAveragingForDirichlet = false
 
 [Component]
 LiquidDensity = 1

test/freeflow/navierstokes/angeli/problem.hh

@@ -26,6 +26,8 @@
 #ifndef DUMUX_ANGELI_TEST_PROBLEM_HH
 #define DUMUX_ANGELI_TEST_PROBLEM_HH
 
+#include <dune/geometry/quadraturerules.hh>
+
 #include <dumux/common/parameters.hh>
 #include <dumux/common/properties.hh>
 #include <dumux/common/numeqvector.hh>
@@ -53,6 +55,7 @@ class AngeliTestProblem : public NavierStokesProblem<TypeTag>
 
     using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
     using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+    using GridView = typename GridGeometry::GridView;
     using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
     using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
     using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
@@ -61,6 +64,8 @@ class AngeliTestProblem : public NavierStokesProblem<TypeTag>
 
     static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
     using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
+    using SubControlVolume = typename GridGeometry::SubControlVolume;
+    using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
     using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
 
@@ -72,23 +77,23 @@ public:
     {
         kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
         rho_ = getParam<Scalar>("Component.LiquidDensity", 1.0);
+        useVelocityAveragingForDirichlet_ = getParam<bool>("Problem.UseVelocityAveragingForDirichlet", false);
+        useVelocityAveragingForInitial_ = getParam<bool>("Problem.UseVelocityAveragingForInitial", false);
     }
 
-   /*!
+    /*!
      * \brief Returns the temperature within the domain in [K].
-     *
-     * This problem assumes a temperature of 10 degrees Celsius.
+     * This problem assumes a temperature of 20 degrees Celsius (unused)
      */
     Scalar temperature() const
-    { return 298.0; }
+    { return 293.15; }
 
-    // \}
-   /*!
+    /*!
      * \name Boundary conditions
      */
     // \{
 
-   /*!
+    /*!
      * \brief Specifies which kind of boundary condition should be
      *        used for which equation on a given boundary control volume.
      *
@@ -126,29 +131,46 @@ public:
         return false;
     }
 
-   /*!
-     * \brief Returns Dirichlet boundary values at a given position.
+    /*!
+     * \brief Evaluate the boundary conditions for a dirichlet
+     *        control volume face (velocities)
      *
-     * \param globalPos The global position
+     * \param element The finite element
+     * \param scvf the sub control volume face
+     *
+     * Concerning the usage of averagedVelocity_, see the explanation of the initial function.
      */
-    PrimaryVariables dirichletAtPos(const GlobalPosition & globalPos) const
+    PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
     {
-        // use the values of the analytical solution
-        return analyticalSolution(globalPos, time_+timeStepSize_);
+        if (useVelocityAveragingForDirichlet_)
+            return averagedVelocity_(scvf);
+        else
+            return analyticalSolution(scvf.center());
     }
 
     /*!
-     * \brief Returns the analytical solution of the problem at a given time and position.
+     * \brief Evaluate the boundary conditions for a dirichlet
+     *        control volume face (pressure)
      *
+     * \param element The finite element
+     * \param scv the sub control volume
+     */
+    PrimaryVariables dirichlet(const Element& element, const SubControlVolume& scv) const
+    {
+        PrimaryVariables priVars(0.0);
+        priVars[Indices::pressureIdx] = analyticalSolution(scv.center())[Indices::pressureIdx];
+        return priVars;
+    }
+
+    /*!
+     * \brief Returns the analytical solution of the problem at a given time and position.
      * \param globalPos The global position
-     * \param time The current simulation time
      */
-    PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, const Scalar time) const
+    PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
     {
         const Scalar x = globalPos[0];
         const Scalar y = globalPos[1];
-
-        const Scalar t = time;
+        const Scalar t = time_;
 
         PrimaryVariables values;
 
@@ -159,54 +181,74 @@ public:
         return values;
     }
 
-    /*!
-     * \brief Returns the analytical solution of the problem at a given position.
-     *
-     * \param globalPos The global position
-     */
-    PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
-    {
-        return analyticalSolution(globalPos, time_+timeStepSize_);
-    }
-
     // \}
 
-   /*!
+    /*!
      * \name Volume terms
      */
     // \{
 
-   /*!
-     * \brief Evaluates the initial value for a control volume.
-     *
-     * \param globalPos The global position
+    /*!
+     * \brief Evaluates the initial value for a control volume (pressure)
      */
-    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+    PrimaryVariables initial(const SubControlVolume& scv) const
     {
-        return analyticalSolution(globalPos, 0);
+        PrimaryVariables priVars(0.0);
+        priVars[Indices::pressureIdx] = analyticalSolution(scv.center())[Indices::pressureIdx];
+        return priVars;
     }
 
     /*!
-     * \brief Updates the time
+     * \brief Evaluates the initial value for a sub control volume face (velocities)
+     *
+     * Simply assigning the value of the analytical solution at the face center
+     * gives a discrete solution that is not divergence-free. For small initial
+     * time steps, this has a negative impact on the pressure solution
+     * after the first time step. The flag UseVelocityAveragingForInitial triggers the
+     * function averagedVelocity_ which uses a higher order quadrature formula to
+     * bring the discrete solution sufficiently close to being divergence-free.
      */
-    void updateTime(const Scalar time)
+    PrimaryVariables initial(const SubControlVolumeFace& scvf) const
     {
-        time_ = time;
+        if (useVelocityAveragingForInitial_)
+            return averagedVelocity_(scvf);
+        else
+            return analyticalSolution(scvf.center());
     }
 
+    // \}
+
     /*!
-     * \brief Updates the time step size
+     * \brief Updates the time information
      */
-    void updateTimeStepSize(const Scalar timeStepSize)
+    void setTime(Scalar t)
     {
-        timeStepSize_ = timeStepSize;
+        time_ = t;
     }
 
 private:
-    Scalar kinematicViscosity_;
-    Scalar rho_;
+    PrimaryVariables averagedVelocity_(const SubControlVolumeFace& scvf) const
+    {
+        PrimaryVariables priVars(0.0);
+        const auto geo = scvf.geometry();
+        const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension-1>::rule(geo.type(), 3);
+        for (auto&& qp : quad)
+        {
+            const auto w = qp.weight()*geo.integrationElement(qp.position());
+            const auto globalPos = geo.global(qp.position());
+            const auto sol = analyticalSolution(globalPos);
+            priVars[Indices::velocityXIdx] += sol[Indices::velocityXIdx]*w;
+            priVars[Indices::velocityYIdx] += sol[Indices::velocityYIdx]*w;
+        }
+        priVars[Indices::velocityXIdx] /= scvf.area();
+        priVars[Indices::velocityYIdx] /= scvf.area();
+        return priVars;
+    }
+
+    Scalar kinematicViscosity_, rho_;
     Scalar time_ = 0;
-    Scalar timeStepSize_ = 0;
+    bool useVelocityAveragingForDirichlet_;
+    bool useVelocityAveragingForInitial_;
 };
 } // end namespace Dumux