[freeflow][test] Added test sincos steady

test/freeflow/navierstokes/sincos/CMakeLists.txt

+add_subdirectory(steady)
 add_subdirectory(unsteady)

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/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 problem the problem for which to evaluate the analytical solution
+*/
+template<class Problem>
+auto createAnalyticalSolution(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 Scalar = double;
+    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);
+
+            // 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);
+                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::SincosSteadyTest;
+
+    // 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();
+
+    // the problem (initial and boundary conditions)
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+    // the solution vector
+    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+    SolutionVector x;
+    x[FVGridGeometry::cellCenterIdx()].resize(fvGridGeometry->numCellCenterDofs());
+    x[FVGridGeometry::faceIdx()].resize(fvGridGeometry->numFaceDofs());
+
+    // 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(*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);
+
+    // 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");
+
+    // linearize & solve
+    Dune::Timer timer;
+    nonLinearSolver.solve(x);
+
+    if (printL2Error)
+    {
+        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;
+    }
+
+    // write vtk output
+    analyticalSolution = createAnalyticalSolution(*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";
+
+    ////////////////////////////////////////////////////////////
+    // 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/navierstokes/sincos/steady/params.input

+[Grid]
+UpperRight = 1 1
+Cells = 50 50
+
+[Problem]
+Name = test_sincos_steady # name passed to the output routines
+EnableGravity = false
+PrintL2Error = true
+EnableInertiaTerms = true
+
+[Component]
+LiquidDensity = 1.0
+LiquidKinematicViscosity = 0.001
+
+[Newton]
+MaxSteps = 10
+MaxRelativeShift = 1e-5
+
+[Vtk]
+WriteFaceData = false
+AddProcessRank = false

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 ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+    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 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