diff --git a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
index c906d0f6e94d4edf1ca4790bb0da3276786971f0..4050e7872263546fbc33299a2bfb3d48d380d1ad 100644
--- a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
+++ b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
@@ -90,7 +90,7 @@ class FreeFlowStaggeredGeometryHelper
static constexpr int codimIntersection = 1;
static constexpr int codimCommonEntity = 2;
static constexpr int numFacetSubEntities = (dim == 2) ? 2 : 4;
- static constexpr int numPairs = (dimWorld == 2) ? 2 : 4;
+ static constexpr int numPairs = 2 * (dimWorld - 1);
public:
@@ -376,6 +376,8 @@ private:
};
Indices indices;
+ if (dim == 1)
+ return indices;
switch(localFacetIdx)
{
diff --git a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
index 5b9cafec19c37ed952fda7baba0fd2ad1fa8941e..e43857366c85f141f8c5f1bf1a0e5d3e8edb2a56 100644
--- a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
@@ -55,7 +55,7 @@ class FreeFlowStaggeredSubControlVolumeFace
static const int dim = Geometry::mydimension;
static const int dimworld = Geometry::coorddimension;
- static constexpr int numPairs = (dimworld == 2) ? 2 : 4;
+ static constexpr int numPairs = 2 * (dimworld - 1);
public:
using GlobalPosition = typename T::GlobalPosition;
diff --git a/test/freeflow/navierstokes/CMakeLists.txt b/test/freeflow/navierstokes/CMakeLists.txt
index 6808d3f93ce7efb8ab4ec0e84d2f5dda1fda474e..6056f26f402b57379fdd8219829da90708993108 100644
--- a/test/freeflow/navierstokes/CMakeLists.txt
+++ b/test/freeflow/navierstokes/CMakeLists.txt
@@ -1,4 +1,5 @@
add_input_file_links()
+dune_symlink_to_source_files(FILES convergence.sh)
set(CMAKE_BUILD_TYPE Debug)
add_executable(test_closedsystem EXCLUDE_FROM_ALL test_closedsystem.cc)
@@ -80,6 +81,15 @@ dune_add_test(NAME test_stokes_donea
${CMAKE_CURRENT_BINARY_DIR}/test_donea-00001.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_stokes_donea")
+dune_add_test(NAME test_navierstokes_1d
+ SOURCES test_navierstokes_1d.cc
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_navierstokes_1d.vtp
+ ${CMAKE_CURRENT_BINARY_DIR}/test_navierstokes_1d-00001.vtp
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_navierstokes_1d")
+
dune_add_test(NAME test_navierstokes_kovasznay
SOURCES test_kovasznay.cc
CMAKE_GUARD HAVE_UMFPACK
diff --git a/test/freeflow/navierstokes/convergence.sh b/test/freeflow/navierstokes/convergence.sh
new file mode 100755
index 0000000000000000000000000000000000000000..945249cbff00987e0d99236f5cbd5878d25f5ae2
--- /dev/null
+++ b/test/freeflow/navierstokes/convergence.sh
@@ -0,0 +1,44 @@
+#! /bin/bash
+
+if [ $# -ne 3 ]; then
+ echo "usage: $0 EXECUTABLE_NAME DIMENSION NUM_REFINEMENT_STEPS"
+ exit
+fi
+
+make $1
+
+LOGFILE=logfile_$1.out
+rm $LOGFILE
+L2ERRORFILE=l2error_$1.txt
+
+if [ $2 == 1 ]; then
+GRIDCELLS="4"
+elif [ $2 == 2 ]; then
+GRIDCELLS="4 4"
+else
+GRIDCELLS="4 4 4"
+fi
+
+MAX=$3
+for (( i=0; i <= $MAX; ++i )); do
+ printf "refinement $i / $MAX "
+ ./$1 $1.input -Grid.Cells "$GRIDCELLS" -Grid.Refinement $i -Problem.PrintL2Error true &>> $LOGFILE
+ echo "done."
+done
+
+grep "L2 error (abs/rel) for" $LOGFILE | tee $L2ERRORFILE
+echo "reset; \
+set log x; \
+set log y; \
+set arrow from graph 0,1 to graph 1,0 nohead lc rgb 'gray'; \
+set arrow from graph 0,1 to graph 1,0.5 nohead lc rgb 'gray'" > $1.gp
+
+PLOT="plot '$L2ERRORFILE' u 6:17 w lp t 'pressure', '$L2ERRORFILE' u 6:23 w lp t 'velocity'"
+if [ $2 == 2 ]; then
+PLOT=$PLOT", '$L2ERRORFILE' u 6:29 w lp t 'velocity'"
+elif [ $2 == 3 ]; then
+PLOT=$PLOT", '$L2ERRORFILE' u 6:35 w lp t 'velocity'"
+fi
+
+echo $PLOT >> $1.gp
+gnuplot --persist $1.gp
diff --git a/test/freeflow/navierstokes/doneatestproblem.hh b/test/freeflow/navierstokes/doneatestproblem.hh
index bbd58ba1bdeb53432cdadba4f4ac5a8ef3798dc4..db0020715bc6b0ef1f30912f87a5b635b76b05d4 100644
--- a/test/freeflow/navierstokes/doneatestproblem.hh
+++ b/test/freeflow/navierstokes/doneatestproblem.hh
@@ -121,8 +121,8 @@ public:
<< 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(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx]
+ << " , L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx]
<< std::endl;
}
}
diff --git a/test/freeflow/navierstokes/kovasznaytestproblem.hh b/test/freeflow/navierstokes/kovasznaytestproblem.hh
index 8927693ffae2bb769b38b5eb6d2325ba43a1a7e1..a3ce7237c0c141b9996bc00c5dcbf33bb2847061 100644
--- a/test/freeflow/navierstokes/kovasznaytestproblem.hh
+++ b/test/freeflow/navierstokes/kovasznaytestproblem.hh
@@ -127,8 +127,8 @@ public:
<< 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(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx]
+ << " , L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx]
<< std::endl;
}
}
diff --git a/test/freeflow/navierstokes/navierstokesanalyticproblem.hh b/test/freeflow/navierstokes/navierstokesanalyticproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..fa6ca53d978591da699bce9885f0105578dd99d7
--- /dev/null
+++ b/test/freeflow/navierstokes/navierstokesanalyticproblem.hh
@@ -0,0 +1,374 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup NavierStokesTests
+ * \brief Test for the 1-D Navier-Stokes model with an analytical solution
+ *
+ * \copydoc NavierStokesAnalyticProblem
+ */
+#ifndef DUMUX_DONEA_TEST_PROBLEM_HH
+#define DUMUX_DONEA_TEST_PROBLEM_HH
+
+#include <dumux/material/components/constant.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/freeflow/navierstokes/model.hh>
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include "l2error.hh"
+
+
+namespace Dumux
+{
+template <class TypeTag>
+class NavierStokesAnalyticProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(NavierStokesAnalyticTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokes));
+
+// the fluid system
+SET_PROP(NavierStokesAnalyticTypeTag, FluidSystem)
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
+};
+
+// Set the grid type
+SET_TYPE_PROP(NavierStokesAnalyticTypeTag, Grid, Dune::YaspGrid<1>);
+
+// Set the problem property
+SET_TYPE_PROP(NavierStokesAnalyticTypeTag, Problem, Dumux::NavierStokesAnalyticProblem<TypeTag> );
+
+SET_BOOL_PROP(NavierStokesAnalyticTypeTag, EnableFVGridGeometryCache, true);
+
+SET_BOOL_PROP(NavierStokesAnalyticTypeTag, EnableGridFluxVariablesCache, true);
+SET_BOOL_PROP(NavierStokesAnalyticTypeTag, EnableGridVolumeVariablesCache, true);
+
+SET_BOOL_PROP(NavierStokesAnalyticTypeTag, EnableInertiaTerms, true);
+SET_BOOL_PROP(NavierStokesAnalyticTypeTag, NormalizePressure, false);
+}
+
+/*!
+ * \ingroup NavierStokesTests
+ * \brief Test for the 1-D Navier-Stokes model with an analytical solution
+ *
+ * The 1-D analytic solution is given by
+ * \f[ p = 2 - 2 \cdot x \f]
+ * \f[ v_\text{x} = 2 \cdot x^3 \f]
+ */
+template <class TypeTag>
+class NavierStokesAnalyticProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+
+ static constexpr auto dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ using DimVector = Dune::FieldVector<Scalar, dimWorld>;
+ using DimMatrix = Dune::FieldVector<Dune::FieldVector<Scalar, dimWorld>, dimWorld>;
+
+public:
+ NavierStokesAnalyticProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry), eps_(1e-6)
+ {
+ printL2Error_ = getParam<bool>("Problem.PrintL2Error");
+ density_ = getParam<Scalar>("Problem.LiquidDensity");
+ kinematicViscosity_ = getParam<Scalar>("Problem.LiquidKinematicViscosity");
+ createAnalyticalSolution_();
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ bool shouldWriteRestartFile() const
+ {
+ return false;
+ }
+
+ void postTimeStep(const SolutionVector& curSol) const
+ {
+ if(printL2Error_)
+ {
+ using L2Error = NavierStokesTestL2Error<Scalar, ModelTraits, PrimaryVariables>;
+ const auto l2error = L2Error::calculateL2Error(*this, curSol);
+ const int numCellCenterDofs = this->fvGridGeometry().numCellCenterDofs();
+ const int numFaceDofs = this->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]
+ << std::endl;
+ }
+ }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature() const
+ { return 298.0; }
+
+ /*!
+ * \brief Return the sources within the domain.
+ *
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+ {
+ NumEqVector source(0.0);
+
+ // mass balance - term div(rho*v)
+ for (unsigned int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
+ {
+ source[Indices::conti0EqIdx] += dvdx(globalPos)[dimIdx][dimIdx];
+ }
+ source[Indices::conti0EqIdx] *= density_;
+
+ // momentum balance
+ for (unsigned int velIdx = 0; velIdx < dimWorld; ++velIdx)
+ {
+ for (unsigned int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
+ {
+ // inertia term
+ if (GET_PROP_VALUE(TypeTag, EnableInertiaTerms))
+ source[Indices::velocity(velIdx)] += density_ * dv2dx(globalPos)[velIdx][dimIdx];
+
+ // viscous term (molecular)
+ source[Indices::velocity(velIdx)] -= density_ * kinematicViscosity_* dvdx2(globalPos)[velIdx][dimIdx];
+ static bool enableUnsymmetrizedVelocityGradient = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup),
+ "FreeFlow.EnableUnsymmetrizedVelocityGradient", false);
+ if (!enableUnsymmetrizedVelocityGradient)
+ source[Indices::velocity(velIdx)] -= density_ * kinematicViscosity_* dvdx2(globalPos)[dimIdx][velIdx];
+ }
+ // pressure term
+ source[Indices::velocity(velIdx)] += dpdx(globalPos)[velIdx];
+
+ // gravity term
+ static bool enableGravity = getParam<bool>("Problem.EnableGravity");
+ if (enableGravity)
+ {
+ source[Indices::velocity(velIdx)] -= density_ * this->gravity()[velIdx];
+ }
+ }
+
+ 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 and pressure everywhere
+ values.setDirichletCell(Indices::conti0EqIdx);
+ values.setDirichlet(Indices::momentumXBalanceIdx);
+
+ return values;
+ }
+
+ /*!
+ * \brief Return dirichlet boundary values at a given position
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+ {
+ // use the values of the analytical solution
+ return analyticalSolution(globalPos);
+ }
+
+ /*!
+ * \brief Return the analytical solution of the problem at a given position
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables values;
+ values[Indices::pressureIdx] = p(globalPos);
+ values[Indices::velocityXIdx] = v(globalPos);
+ return values;
+ }
+
+ //! \brief The velocity
+ const DimVector v(const DimVector& globalPos) const
+ {
+ DimVector v(0.0);
+ v[0] = 2.0 * globalPos[0] * globalPos[0] * globalPos[0];
+ return v;
+ }
+
+ //! \brief The velocity gradient
+ const DimMatrix dvdx(const DimVector& globalPos) const
+ {
+ DimMatrix dvdx(0.0);
+ dvdx[0][0] = 6.0 * globalPos[0] * globalPos[0];
+ return dvdx;
+ }
+
+ //! \brief The gradient of the velocity squared (using product rule -> nothing to do here)
+ const DimMatrix dv2dx(const DimVector& globalPos) const
+ {
+ DimMatrix dv2dx;
+ for (unsigned int velIdx = 0; velIdx < dimWorld; ++velIdx)
+ {
+ for (unsigned int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
+ {
+ dv2dx[velIdx][dimIdx] = dvdx(globalPos)[velIdx][dimIdx] * v(globalPos)[dimIdx]
+ + dvdx(globalPos)[dimIdx][dimIdx] * v(globalPos)[velIdx];
+ }
+ }
+ return dv2dx;
+ }
+
+ //! \brief The gradient of the velocity gradient
+ const DimMatrix dvdx2(const DimVector& globalPos) const
+ {
+ DimMatrix dvdx2(0.0);
+ dvdx2[0][0] = 12.0 * globalPos[0];
+ return dvdx2;
+ }
+
+ //! \brief The pressure
+ const Scalar p(const DimVector& globalPos) const
+ { return 2.0 - 2.0 * globalPos[0]; }
+
+ //! \brief The pressure gradient
+ const DimVector dpdx(const DimVector& globalPos) const
+ {
+ DimVector dpdx(0.0);
+ dpdx[0] = -2.0;
+ return dpdx;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
+ {
+ return analyticalSolution(globalPos);
+ }
+
+ /*!
+ * \brief Returns the analytical solution for the pressure
+ */
+ auto& getAnalyticalPressureSolution() const
+ {
+ return analyticalPressure_;
+ }
+
+ /*!
+ * \brief Returns the analytical solution for the velocity
+ */
+ auto& getAnalyticalVelocitySolution() const
+ {
+ return analyticalVelocity_;
+ }
+
+ /*!
+ * \brief Returns the analytical solution for the velocity at the faces
+ */
+ auto& getAnalyticalVelocitySolutionOnFace() const
+ {
+ return analyticalVelocityOnFace_;
+ }
+
+private:
+
+ /*!
+ * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
+ */
+ void createAnalyticalSolution_()
+ {
+ analyticalPressure_.resize(this->fvGridGeometry().numCellCenterDofs());
+ analyticalVelocity_.resize(this->fvGridGeometry().numCellCenterDofs());
+ analyticalVelocityOnFace_.resize(this->fvGridGeometry().numFaceDofs());
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ auto ccDofIdx = scv.dofIndex();
+ auto ccDofPosition = scv.dofPosition();
+ auto analyticalSolutionAtCc = analyticalSolution(ccDofPosition);
+
+ // velocities on faces
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ const auto faceDofIdx = scvf.dofIndex();
+ const auto faceDofPosition = scvf.center();
+ const auto dirIdx = scvf.directionIndex();
+ const auto analyticalSolutionAtFace = analyticalSolution(faceDofPosition);
+ analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
+ }
+
+ analyticalPressure_[ccDofIdx] = analyticalSolutionAtCc[Indices::pressureIdx];
+
+ for(int dirIdx = 0; dirIdx < ModelTraits::dim(); ++dirIdx)
+ analyticalVelocity_[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)];
+ }
+ }
+ }
+
+ Scalar eps_;
+ bool printL2Error_;
+ Scalar density_;
+ Scalar kinematicViscosity_;
+ std::vector<Scalar> analyticalPressure_;
+ std::vector<GlobalPosition> analyticalVelocity_;
+ std::vector<GlobalPosition> analyticalVelocityOnFace_;
+};
+} //end namespace
+
+#endif
diff --git a/test/freeflow/navierstokes/test_navierstokes_1d.cc b/test/freeflow/navierstokes/test_navierstokes_1d.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4d8c4598ea6f7bec49abb8153b098e27f709c9eb
--- /dev/null
+++ b/test/freeflow/navierstokes/test_navierstokes_1d.cc
@@ -0,0 +1,197 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Test for a 1-D staggered grid Navier-Stokes model
+ */
+
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+#include <dune/istl/io.hh>
+
+#include <dumux/common/defaultusagemessage.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+
+#include <dumux/assembly/staggeredfvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+
+#include "navierstokesanalyticproblem.hh"
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n"
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n";
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(NavierStokesAnalyticTypeTag);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv, usage);
+
+ // try to create a grid (from the given grid file or the input file)
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ GridCreator::makeGrid();
+ GridCreator::loadBalance();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = GridCreator::grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ const auto numDofsCellCenter = leafGridView.size(0);
+ const auto numDofsFace = leafGridView.size(1);
+ SolutionVector x;
+ x[FVGridGeometry::cellCenterIdx()].resize(numDofsCellCenter);
+ x[FVGridGeometry::faceIdx()].resize(numDofsFace);
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x);
+
+ // intialize the vtk output module
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ StaggeredVtkOutputModule<TypeTag, GET_PROP_VALUE(TypeTag, PhaseIdx)> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
+ vtkWriter.addField(problem->getAnalyticalPressureSolution(), "pressureExact");
+ vtkWriter.addField(problem->getAnalyticalVelocitySolution(), "velocityExact");
+ vtkWriter.addFaceField(problem->getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
+ vtkWriter.write(0.0);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables);
+
+ // 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);
+
+ // linearize & solve
+ Dune::Timer timer;
+ nonLinearSolver.solve(x);
+
+ // write vtk output
+ problem->postTimeStep(x);
+ 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;
+}
diff --git a/test/freeflow/navierstokes/test_navierstokes_1d.input b/test/freeflow/navierstokes/test_navierstokes_1d.input
new file mode 100644
index 0000000000000000000000000000000000000000..2215feafa6a633f0d67f4421f21c26cfdda8926b
--- /dev/null
+++ b/test/freeflow/navierstokes/test_navierstokes_1d.input
@@ -0,0 +1,18 @@
+[Grid]
+UpperRight = 1
+Cells = 32
+
+[Problem]
+Name = test_navierstokes_1d
+EnableGravity = false
+PrintL2Error = true
+LiquidKinematicViscosity = 1.
+LiquidDensity = 1.
+
+[Newton]
+MaxSteps = 10
+MaxRelativeShift = 1e-9
+
+[Vtk]
+WriteFaceData = false
+AddVelocity = true
diff --git a/test/references/test_navierstokes_1d.vtp b/test/references/test_navierstokes_1d.vtp
new file mode 100644
index 0000000000000000000000000000000000000000..55ed9f74ae3516156328c29cd479a76e795d127c
--- /dev/null
+++ b/test/references/test_navierstokes_1d.vtp
@@ -0,0 +1,72 @@
+<?xml version="1.0"?>
+<VTKFile type="PolyData" version="0.1" byte_order="LittleEndian">
+ <PolyData>
+ <Piece NumberOfLines="32" NumberOfPoints="33">
+ <CellData Scalars="p">
+ <DataArray type="Float32" Name="p" NumberOfComponents="1" format="ascii">
+ 1.96875 1.81039 1.74788 1.68538 1.62288 1.56038 1.49789 1.43542 1.37299 1.31061 1.2483 1.18611
+ 1.12406 1.06219 1.00058 0.939263 0.878327 0.817848 0.757919 0.698646 0.640142 0.582539 0.525978 0.470617
+ 0.416626 0.364194 0.313523 0.264835 0.218367 0.174376 0.133136 0.03125
+ </DataArray>
+ <DataArray type="Float32" Name="rhoMolar" NumberOfComponents="1" format="ascii">
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1
+ </DataArray>
+ <DataArray type="Float32" Name="rho" NumberOfComponents="1" format="ascii">
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1
+ </DataArray>
+ <DataArray type="Float32" Name="velocity_component (m/s)" NumberOfComponents="1" format="ascii">
+ 0.000771803 0.0017496 0.0025278 0.00422152 0.00719699 0.0118204 0.018458 0.0274759 0.0392404 0.0541178 0.0724741 0.0946756
+ 0.121089 0.152079 0.188014 0.229258 0.276179 0.329142 0.388514 0.454661 0.527949 0.608744 0.697413 0.794322
+ 0.899836 1.01432 1.13815 1.27168 1.41528 1.56932 1.73416 1.90966
+ </DataArray>
+ <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0
+ </DataArray>
+ <DataArray type="Float32" Name="pressureExact" NumberOfComponents="1" format="ascii">
+ 1.96875 1.90625 1.84375 1.78125 1.71875 1.65625 1.59375 1.53125 1.46875 1.40625 1.34375 1.28125
+ 1.21875 1.15625 1.09375 1.03125 0.96875 0.90625 0.84375 0.78125 0.71875 0.65625 0.59375 0.53125
+ 0.46875 0.40625 0.34375 0.28125 0.21875 0.15625 0.09375 0.03125
+ </DataArray>
+ <DataArray type="Float32" Name="velocityExact" NumberOfComponents="1" format="ascii">
+ 7.62939e-06 0.000205994 0.000953674 0.00261688 0.00556183 0.0101547 0.0167618 0.0257492 0.0374832 0.05233 0.0706558 0.0928268
+ 0.119209 0.150169 0.186073 0.227287 0.274178 0.32711 0.386452 0.452568 0.525826 0.60659 0.695229 0.792107
+ 0.897591 1.01205 1.13584 1.26934 1.41291 1.56692 1.73173 1.90771
+ </DataArray>
+ </CellData>
+ <Points>
+ <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+ 0 0 0 0.03125 0 0 0.0625 0 0 0.09375 0 0
+ 0.125 0 0 0.15625 0 0 0.1875 0 0 0.21875 0 0
+ 0.25 0 0 0.28125 0 0 0.3125 0 0 0.34375 0 0
+ 0.375 0 0 0.40625 0 0 0.4375 0 0 0.46875 0 0
+ 0.5 0 0 0.53125 0 0 0.5625 0 0 0.59375 0 0
+ 0.625 0 0 0.65625 0 0 0.6875 0 0 0.71875 0 0
+ 0.75 0 0 0.78125 0 0 0.8125 0 0 0.84375 0 0
+ 0.875 0 0 0.90625 0 0 0.9375 0 0 0.96875 0 0
+ 1 0 0
+ </DataArray>
+ </Points>
+ <Lines>
+ <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+ 0 1 1 2 2 3 3 4 4 5 5 6
+ 6 7 7 8 8 9 9 10 10 11 11 12
+ 12 13 13 14 14 15 15 16 16 17 17 18
+ 18 19 19 20 20 21 21 22 22 23 23 24
+ 24 25 25 26 26 27 27 28 28 29 29 30
+ 30 31 31 32
+ </DataArray>
+ <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+ 2 4 6 8 10 12 14 16 18 20 22 24
+ 26 28 30 32 34 36 38 40 42 44 46 48
+ 50 52 54 56 58 60 62 64
+ </DataArray>
+ </Lines>
+ </Piece>
+ </PolyData>
+</VTKFile>