From 902a81689ec546730eac3d3a6be4dc5909bc88da Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Wed, 10 Jun 2020 18:57:18 +0200
Subject: [PATCH] [ff][test] Add rotational-symmetric pipe flow test
---
.../navierstokes/channel/CMakeLists.txt | 1 +
.../navierstokes/channel/pipe/CMakeLists.txt | 8 +
.../channel/pipe/convergencetest.py | 81 +++++++++
.../navierstokes/channel/pipe/main.cc | 159 ++++++++++++++++++
.../navierstokes/channel/pipe/params.input | 22 +++
.../navierstokes/channel/pipe/problem.hh | 150 +++++++++++++++++
.../navierstokes/channel/pipe/properties.hh | 83 +++++++++
test/freeflow/navierstokes/l2error.hh | 12 +-
8 files changed, 512 insertions(+), 4 deletions(-)
create mode 100644 test/freeflow/navierstokes/channel/pipe/CMakeLists.txt
create mode 100755 test/freeflow/navierstokes/channel/pipe/convergencetest.py
create mode 100644 test/freeflow/navierstokes/channel/pipe/main.cc
create mode 100644 test/freeflow/navierstokes/channel/pipe/params.input
create mode 100644 test/freeflow/navierstokes/channel/pipe/problem.hh
create mode 100644 test/freeflow/navierstokes/channel/pipe/properties.hh
diff --git a/test/freeflow/navierstokes/channel/CMakeLists.txt b/test/freeflow/navierstokes/channel/CMakeLists.txt
index 642e957dc2..31ded5a8b0 100644
--- a/test/freeflow/navierstokes/channel/CMakeLists.txt
+++ b/test/freeflow/navierstokes/channel/CMakeLists.txt
@@ -1,3 +1,4 @@
add_subdirectory(1d)
add_subdirectory(2d)
add_subdirectory(3d)
+add_subdirectory(pipe)
diff --git a/test/freeflow/navierstokes/channel/pipe/CMakeLists.txt b/test/freeflow/navierstokes/channel/pipe/CMakeLists.txt
new file mode 100644
index 0000000000..b6df8b1811
--- /dev/null
+++ b/test/freeflow/navierstokes/channel/pipe/CMakeLists.txt
@@ -0,0 +1,8 @@
+dune_symlink_to_source_files(FILES params.input convergencetest.py)
+dumux_add_test(NAME test_ff_stokes_channel_pipe
+ LABELS freeflow navierstokes
+ SOURCES main.cc
+ LABELS freeflow
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ./convergencetest.py
+ CMD_ARGS test_ff_stokes_channel_pipe)
diff --git a/test/freeflow/navierstokes/channel/pipe/convergencetest.py b/test/freeflow/navierstokes/channel/pipe/convergencetest.py
new file mode 100755
index 0000000000..24750a8962
--- /dev/null
+++ b/test/freeflow/navierstokes/channel/pipe/convergencetest.py
@@ -0,0 +1,81 @@
+#!/usr/bin/env python3
+
+from math import *
+import subprocess
+import sys
+
+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, 2]:
+ 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'])
+
+def mean(numbers):
+ return float(sum(numbers)) / len(numbers)
+
+# check the rates, we expect rates around 2
+if mean(resultsP) < 2.1 and mean(resultsP) < 1.8:
+ 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 mean(resultsVx) < 2.05 and mean(resultsVx) < 1.9:
+ 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 mean(resultsVy) < 2.05 and mean(resultsVy) < 1.9:
+ 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)
diff --git a/test/freeflow/navierstokes/channel/pipe/main.cc b/test/freeflow/navierstokes/channel/pipe/main.cc
new file mode 100644
index 0000000000..4ca8e828d6
--- /dev/null
+++ b/test/freeflow/navierstokes/channel/pipe/main.cc
@@ -0,0 +1,159 @@
+// -*- 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/>. *
+ *****************************************************************************/
+
+#include <config.h>
+
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/partial.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager_yasp.hh>
+
+#include <dumux/assembly/staggeredfvassembler.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+
+#include "properties.hh"
+#include "../../l2error.hh"
+
+template<class Problem, class SolutionVector, class GridGeometry>
+void printL2Error(const Problem& problem, const SolutionVector& x, const GridGeometry& gridGeometry)
+{
+ using namespace Dumux;
+ using Scalar = double;
+ using TypeTag = Properties::TTag::PipeFlow;
+ 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 [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) = " << l2errorAbs[Indices::pressureIdx] << " / " << l2errorRel[Indices::pressureIdx]
+ << " , L2(vx) = " << l2errorAbs[Indices::velocityXIdx] << " / " << l2errorRel[Indices::velocityXIdx]
+ << " , L2(vy) = " << l2errorAbs[Indices::velocityYIdx] << " / " << l2errorRel[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) = " << l2errorAbs[Indices::pressureIdx] << " L2(vx) = " << l2errorAbs[Indices::velocityXIdx] << " L2(vy) = " << l2errorAbs[Indices::velocityYIdx] << std::endl;
+ logFile.close();
+}
+
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // Define the sub problem type tags
+ using TypeTag = Properties::TTag::PipeFlow;
+
+ // try to create a grid (from the given grid file or the input file)
+ Dumux::GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
+ gridManager.init();
+
+ // we compute on the leaf grid view
+ const auto& gridView = gridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ auto gridGeometry = std::make_shared<GridGeometry>(gridView);
+ gridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ auto problem = std::make_shared<Problem>(gridGeometry);
+
+ // the solution vector
+ using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+ SolutionVector sol;
+ sol[GridGeometry::cellCenterIdx()].resize(gridGeometry->numCellCenterDofs());
+ sol[GridGeometry::faceIdx()].resize(gridGeometry->numFaceDofs());
+ problem->applyInitialSolution(sol);
+
+ // the grid variables
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
+ gridVariables->init(sol);
+
+ // intialize the vtk output module
+ using IOFields = GetPropType<TypeTag, Properties::IOFields>;
+ StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, sol, problem->name());
+ IOFields::initOutputModule(vtkWriter); // Add model specific output fields
+ 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, gridGeometry, 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);
+
+ nonLinearSolver.solve(sol);
+ vtkWriter.write(1.0);
+
+ printL2Error(problem, sol, gridGeometry);
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ Parameters::print();
+
+ return 0;
+
+} // end main
+catch (const Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (const 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 (const Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
diff --git a/test/freeflow/navierstokes/channel/pipe/params.input b/test/freeflow/navierstokes/channel/pipe/params.input
new file mode 100644
index 0000000000..268e5b9bf9
--- /dev/null
+++ b/test/freeflow/navierstokes/channel/pipe/params.input
@@ -0,0 +1,22 @@
+[Grid]
+Positions0 = 0 0.1e-3 # [m]
+Positions1 = 0 1.0e-3 # [m]
+Cells0 = 10 # [-]
+Cells1 = 100 # [-]
+Grading0 = 1
+Grading1 = 1
+
+[Problem]
+Name = test_ff_stokes_channel_pipe
+MeanInletVelocity = 0.1 # [m/s]
+
+[Problem]
+EnableGravity = false
+EnableInertiaTerms = false
+
+[Component]
+LiquidKinematicViscosity = 1e-6
+LiquidDensity = 1000
+
+[Vtk]
+AddVelocity = 1
diff --git a/test/freeflow/navierstokes/channel/pipe/problem.hh b/test/freeflow/navierstokes/channel/pipe/problem.hh
new file mode 100644
index 0000000000..cf5cc443b9
--- /dev/null
+++ b/test/freeflow/navierstokes/channel/pipe/problem.hh
@@ -0,0 +1,150 @@
+// -*- 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/>. *
+ *****************************************************************************/
+
+#ifndef DUMUX_TEST_FREEFLOW_PIPE_PROBLEM_HH
+#define DUMUX_TEST_FREEFLOW_PIPE_PROBLEM_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/freeflow/navierstokes/boundarytypes.hh>
+
+namespace Dumux {
+
+/*!
+ * \brief Freeflow problem for pipe flow
+ * Simulation of a radially-symmetric pipe flow with circular cross-section
+ */
+template <class TypeTag>
+class FreeFlowPipeProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using GridView = typename GridGeometry::GridView;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using BoundaryTypes = NavierStokesBoundaryTypes<PrimaryVariables::size()>;
+
+public:
+ FreeFlowPipeProblem(std::shared_ptr<const GridGeometry> gridGeometry)
+ : ParentType(gridGeometry)
+
+ {
+ name_ = getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
+ meanInletVelocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MeanInletVelocity");
+ mu_ = getParam<Scalar>("Component.LiquidKinematicViscosity")*getParam<Scalar>("Component.LiquidDensity");
+
+ pipeRadius_ = this->gridGeometry().bBoxMax()[0] - this->gridGeometry().bBoxMin()[0];
+ pipeLength_ = this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1];
+ eps_ = 1e-7*pipeRadius_;
+
+ std::cout << "-- Reynolds number: " << 2*pipeRadius_*meanInletVelocity_/getParam<Scalar>("Component.LiquidKinematicViscosity") << std::endl;
+ }
+
+ const std::string& name() const
+ {
+ return name_;
+ }
+
+ Scalar temperature() const
+ { return 293.15; }
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ */
+ BoundaryTypes boundaryTypes(const Element& element,
+ const SubControlVolumeFace& scvf) const
+ {
+ BoundaryTypes values;
+ const auto& globalPos = scvf.dofPosition();
+
+ // inlet
+ if (onLowerBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+ // outlet
+ else if (onUpperBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ }
+ // pipe centerline
+ else if (onInnerBoundary_(globalPos))
+ {
+ values.setAllSymmetry();
+ }
+ // pipe wall
+ else if (onOuterBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+
+ return values;
+ }
+
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+ { return analyticalSolution(globalPos); }
+
+ PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
+ { return analyticalSolution(globalPos); }
+
+ PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables values(0.0);
+
+ // paraboloid velocity profile
+ const auto r = globalPos[0] - this->gridGeometry().bBoxMin()[0];
+ const auto y = globalPos[1] - this->gridGeometry().bBoxMin()[1];
+ values[Indices::velocityXIdx] = 0.0;
+ values[Indices::velocityYIdx] = 2.0*meanInletVelocity_*(1.0 - r*r/(pipeRadius_*pipeRadius_));
+ values[Indices::pressureIdx] = 1e5 + (pipeLength_-y)*meanInletVelocity_*8.0*mu_/(pipeRadius_*pipeRadius_);
+
+ return values;
+ }
+
+private:
+ bool onInnerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onOuterBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
+
+ std::string name_;
+ Scalar meanInletVelocity_;
+ Scalar mu_;
+ Scalar pipeRadius_, pipeLength_;
+ Scalar eps_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/freeflow/navierstokes/channel/pipe/properties.hh b/test/freeflow/navierstokes/channel/pipe/properties.hh
new file mode 100644
index 0000000000..9de9fd669a
--- /dev/null
+++ b/test/freeflow/navierstokes/channel/pipe/properties.hh
@@ -0,0 +1,83 @@
+// -*- 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/>. *
+ *****************************************************************************/
+
+#ifndef DUMUX_TEST_FREEFLOW_PIPE_PROPERTIES_HH
+#define DUMUX_TEST_FREEFLOW_PIPE_PROPERTIES_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/common/properties.hh>
+
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/discretization/extrusion.hh>
+#include <dumux/freeflow/navierstokes/model.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/material/components/constant.hh>
+
+#include "problem.hh"
+
+namespace Dumux::Properties {
+
+// Create new type tags
+namespace TTag {
+struct PipeFlow { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
+} // end namespace TTag
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::PipeFlow>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::Constant<1, Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::PipeFlow>
+{ using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::PipeFlow>
+{ using type = FreeFlowPipeProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableGridGeometryCache<TypeTag, TTag::PipeFlow> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::PipeFlow> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::PipeFlow> { static constexpr bool value = true; };
+
+// rotation-symmetric grid geometry forming a cylinder channel
+template<class TypeTag>
+struct GridGeometry<TypeTag, TTag::PipeFlow>
+{
+ static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
+ static constexpr bool enableCache = getPropValue<TypeTag, Properties::EnableGridGeometryCache>();
+ using GridView = typename GetPropType<TypeTag, Properties::Grid>::LeafGridView;
+
+ struct GGTraits : public StaggeredFreeFlowDefaultFVGridGeometryTraits<GridView, upwindSchemeOrder>
+ { using Extrusion = RotationalExtrusion<0>; };
+
+ using type = StaggeredFVGridGeometry<GridView, enableCache, GGTraits>;
+};
+
+} // end namespace Dumux::Properties
+
+#endif
diff --git a/test/freeflow/navierstokes/l2error.hh b/test/freeflow/navierstokes/l2error.hh
index 8c001783b3..ec418b0ab8 100644
--- a/test/freeflow/navierstokes/l2error.hh
+++ b/test/freeflow/navierstokes/l2error.hh
@@ -26,6 +26,7 @@
#include <vector>
#include <cmath>
+#include <dumux/discretization/extrusion.hh>
namespace Dumux {
@@ -50,6 +51,7 @@ public:
static auto calculateL2Error(const Problem& problem, const SolutionVector& curSol)
{
using GridGeometry = std::decay_t<decltype(problem.gridGeometry())>;
+ using Extrusion = Extrusion_t<GridGeometry>;
PrimaryVariables sumError(0.0), sumReference(0.0), l2NormAbs(0.0), l2NormRel(0.0);
const int numFaceDofs = problem.gridGeometry().numFaceDofs();
@@ -73,9 +75,9 @@ public:
const auto& posCellCenter = scv.dofPosition();
const auto analyticalSolutionCellCenter = problem.analyticalSolution(posCellCenter)[Indices::pressureIdx];
const auto numericalSolutionCellCenter = curSol[GridGeometry::cellCenterIdx()][dofIdxCellCenter][Indices::pressureIdx - ModelTraits::dim()];
- sumError[Indices::pressureIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * scv.volume();
- sumReference[Indices::pressureIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * scv.volume();
- totalVolume += scv.volume();
+ sumError[Indices::pressureIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * Extrusion::volume(scv);
+ sumReference[Indices::pressureIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * Extrusion::volume(scv);
+ totalVolume += Extrusion::volume(scv);
// treat face dofs
for (auto&& scvf : scvfs(fvGeometry))
@@ -87,7 +89,9 @@ public:
directionIndex[dofIdxFace] = dirIdx;
errorVelocity[dofIdxFace] = squaredDiff_(analyticalSolutionFace, numericalSolutionFace);
velocityReference[dofIdxFace] = squaredDiff_(analyticalSolutionFace, 0.0);
- const Scalar staggeredHalfVolume = 0.5 * scv.volume();
+ auto faceScvCenter = scv.center() + scvf.center(); faceScvCenter *= 0.5;
+ typename GridGeometry::Traits::FaceSubControlVolume faceScv(faceScvCenter, 0.5*scv.volume());
+ const Scalar staggeredHalfVolume = Extrusion::volume(faceScv);
staggeredVolume[dofIdxFace] = staggeredVolume[dofIdxFace] + staggeredHalfVolume;
}
}
--
GitLab