Commit df4b7b3f authored by Timo Koch's avatar Timo Koch
Browse files

[test][ff][diamond] Add DGF benchmark test flow around cylinder

parent 93822f60
Pipeline #20015 waiting for manual action with stages
in 6 seconds
......@@ -4,5 +4,6 @@ add_subdirectory(donea)
add_subdirectory(kovasznay)
add_subdirectory(periodic)
add_subdirectory(sincos)
add_subdirectory(unstructured)
dune_symlink_to_source_files(FILES convergence.sh)
dune_symlink_to_source_files(FILES params.input cylinder_channel.msh cylinder_channel_quad.msh)
# Navier-Stokes version of the test (Re=20)
dumux_add_test(NAME test_ff_navierstokes_dfg_benchmark_stationary
SOURCES main.cc
LABELS freeflow navierstokes
CMAKE_GUARD "( HAVE_UMFPACK AND dune-uggrid_FOUND )"
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_navierstokes_dfg_benchmark_stationary.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_dfg_benchmark_stationary-00001.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_dfg_benchmark_stationary params.input
-Problem.Name test_ff_navierstokes_dfg_benchmark_stationary -Problem.EnableInertiaTerms true")
# Stokes version of the test (no inertia)
dumux_add_test(NAME test_ff_stokes_dfg_benchmark_stationary
TARGET test_ff_navierstokes_dfg_benchmark_stationary
LABELS freeflow navierstokes
CMAKE_GUARD "( HAVE_UMFPACK AND dune-uggrid_FOUND )"
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_stokes_dfg_benchmark_stationary.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_ff_stokes_dfg_benchmark_stationary-00001.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_dfg_benchmark_stationary params.input
-Problem.Name test_ff_stokes_dfg_benchmark_stationary -Problem.EnableInertiaTerms false")
// channel measurements according to
// the DFG benchmark
//
// domain mesh size
cl_ = 0.01;
// cylinder surface mesh size
cl2_ = 0.002;
//+
Point(1) = {0.0, 0.0, 0, cl_};
//+
Point(2) = {0.0, 0.0, 0, cl_};
//+
Point(3) = {2.2, 0, 0, cl_};
//+
Point(4) = {0.0, 0.41, 0.0, cl_};
//+
Point(5) = {2.2, 0.41, 0.0, cl_};
//+
Point(6) = {0.15, 0.2, 0.0, cl2_};
//+
Point(7) = {0.25, 0.2, 0.0, cl2_};
//+
Point(8) = {0.2, 0.15, 0.0, cl2_};
//+
Point(9) = {0.2, 0.25, 0.0, cl2_};
//+
Point(10) = {0.2, 0.2, 0, cl2_};
//+
Line(1) = {4, 5};
//+
Line(2) = {5, 3};
//+
Line(3) = {3, 1};
//+
Line(4) = {1, 4};
//+
Circle(5) = {6, 10, 9};
//+
Circle(6) = {9, 10, 7};
//+
Circle(7) = {7, 10, 8};
//+
Circle(8) = {8, 10, 6};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Curve Loop(2) = {6, 7, 8, 5};
//+
Plane Surface(1) = {1, 2};
This diff is collapsed.
// channel measurements according to
// the DFG benchmark
// using quad mesh
Include "cylinder_channel.geo";
Recombine Surface {1};
// -*- 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 <random>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dumux/common/initialize.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager_ug.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <dumux/linear/incompressiblestokessolver.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/traits.hh>
#include <dumux/multidomain/newtonsolver.hh>
#include <dumux/freeflow/navierstokes/velocityoutput.hh>
#include "properties.hh"
int main(int argc, char** argv)
{
using namespace Dumux;
// define the type tag for this problem
using MomentumTypeTag = Properties::TTag::DFGChannelTestMomentum;
using MassTypeTag = Properties::TTag::DFGChannelTestMass;
// initialize MPI, finalize is done automatically on exit
initialize(argc, argv);
const auto& mpiHelper = Dune::MPIHelper::instance();
// print dumux start message
if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/true);
// parse command line arguments and input file
Parameters::init(argc, argv);
// create a grid
using Grid = GetPropType<MassTypeTag, Properties::Grid>;
Dumux::GridManager<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 MomentumGridGeometry = GetPropType<MomentumTypeTag, Properties::GridGeometry>;
auto momentumGridGeometry = std::make_shared<MomentumGridGeometry>(leafGridView);
using MassGridGeometry = GetPropType<MassTypeTag, Properties::GridGeometry>;
auto massGridGeometry = std::make_shared<MassGridGeometry>(leafGridView);
// the coupling manager
using Traits = MultiDomainTraits<MomentumTypeTag, MassTypeTag>;
using CouplingManager = StaggeredFreeFlowCouplingManager<Traits>;
auto couplingManager = std::make_shared<CouplingManager>();
// the problems (boundary conditions)
using MomentumProblem = GetPropType<MomentumTypeTag, Properties::Problem>;
auto momentumProblem = std::make_shared<MomentumProblem>(momentumGridGeometry, couplingManager);
using MassProblem = GetPropType<MassTypeTag, Properties::Problem>;
auto massProblem = std::make_shared<MassProblem>(massGridGeometry, couplingManager);
// the solution vector
constexpr auto momentumIdx = CouplingManager::freeFlowMomentumIndex;
constexpr auto massIdx = CouplingManager::freeFlowMassIndex;
using SolutionVector = typename Traits::SolutionVector;
SolutionVector x;
x[momentumIdx].resize(momentumGridGeometry->numDofs());
x[massIdx].resize(massGridGeometry->numDofs());
std::cout << "Total number of dofs: "
<< massGridGeometry->numDofs() + momentumGridGeometry->numDofs()*MomentumGridGeometry::GridView::dimension
<< std::endl;
// the grid variables
using MomentumGridVariables = GetPropType<MomentumTypeTag, Properties::GridVariables>;
auto momentumGridVariables = std::make_shared<MomentumGridVariables>(momentumProblem, momentumGridGeometry);
using MassGridVariables = GetPropType<MassTypeTag, Properties::GridVariables>;
auto massGridVariables = std::make_shared<MassGridVariables>(massProblem, massGridGeometry);
// compute coupling stencil and afterwards initialize grid variables (need coupling information)
couplingManager->init(momentumProblem, massProblem, std::make_tuple(momentumGridVariables, massGridVariables), x);
massGridVariables->init(x[massIdx]);
momentumGridVariables->init(x[momentumIdx]);
using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(std::make_tuple(momentumProblem, massProblem),
std::make_tuple(momentumGridGeometry, massGridGeometry),
std::make_tuple(momentumGridVariables, massGridVariables),
couplingManager);
// initialize the vtk output module
using IOFields = GetPropType<MassTypeTag, Properties::IOFields>;
VtkOutputModule vtkWriter(*massGridVariables, x[massIdx], massProblem->name());
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
vtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<MassGridVariables>>());
vtkWriter.write(0.0);
// the linear solver
using M = typename Assembler::JacobianMatrix; using V = typename Assembler::ResidualType;
using LinearSolver = IncompressibleStokesSolver<M, V, MassGridGeometry, LinearSolverTraits<MomentumGridGeometry>>;
//using LinearSolver = UMFPackBackend;
V dirichletDofs;
dirichletDofs[momentumIdx].resize(momentumGridGeometry->numDofs());
dirichletDofs[massIdx].resize(massGridGeometry->numDofs());
for (const auto& element : elements(leafGridView))
{
const auto fvGeometry = localView(*momentumGridGeometry).bind(element);
for (const auto& scvf : scvfs(fvGeometry))
{
if (scvf.boundary())
{
const auto bcTypes = momentumProblem->boundaryTypes(element, scvf);
for (int i = 0; i < bcTypes.size(); ++i)
if (bcTypes.isDirichlet(i))
dirichletDofs[momentumIdx][fvGeometry.scv(scvf.insideScvIdx()).dofIndex()][i] = 1.0;
}
}
}
auto linearSolver = std::make_shared<LinearSolver>(massGridGeometry, dirichletDofs);
// the non-linear solver
using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
// linearize & solve
Dune::Timer timer;
nonLinearSolver.solve(x);
// write vtk output
vtkWriter.write(1.0);
// update coupling manager for output
couplingManager->updateSolution(x);
// evaluate benchmark indicators (only valid with inertia term at Re=20)
if (momentumProblem->enableInertiaTerms())
{
static constexpr double cDrafReference = 5.57953523384;
static constexpr double cLiftReference = 0.010618948146;
static constexpr double pDiffReference = 0.11752016697;
const auto [cDrag, cLift] = momentumProblem->evalDragAndLiftCoefficient(*momentumGridVariables, x[momentumIdx]);
std::cout << "cDrag: " << cDrag
<< " (reference: " << cDrafReference << ")"
<< "\n"
<< "cLift: " << cLift
<< " (reference: " << cLiftReference << ")"
<< std::endl;
const auto pDiff = massProblem->evalPressureDifference(x[massIdx]);
std::cout << "pDiff: " << pDiff
<< " (reference: " << pDiffReference << ")"
<< std::endl;
if (getParam<bool>("Problem.CheckIndicators", true))
{
using std::abs;
if (abs(cDrag - cDrafReference) > 0.002)
DUNE_THROW(Dune::Exception,
"Deviation from reference too large for drag coefficient: "
<< cDrag << " (ref: " << cDrafReference << ")"
);
if (abs(cLift - cLiftReference) > 0.002)
DUNE_THROW(Dune::Exception,
"Deviation from reference too large for lift coefficient: "
<< cLift << " (ref: " << cLiftReference << ")"
);
if (abs(pDiff - pDiffReference) > 0.003)
DUNE_THROW(Dune::Exception,
"Deviation from reference too large for pressure difference: "
<< pDiff << " (ref: " << pDiffReference << ")"
);
}
}
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;
}
[Grid]
#File = cylinder_channel.msh
LowerLeft = 0 0
UpperRight = 1 1
Cells = 10 10
[Problem]
Name = benchmark_stationary
EnableGravity = false
EnableInertiaTerms = false
Alpha = 0.1
Beta = 2.0
[FreeFlow]
EnableUnsymmetrizedVelocityGradient = true
[Flux]
UpwindWeight = 0.5
[Component]
LiquidDensity = 1
LiquidDynamicViscosity = 1e-3
[Mass.Assembly.NumericDifference]
PriVarMagnitude = 1e-2
BaseEpsilon = 100
[Momentum.Assembly.NumericDifference]
PriVarMagnitude = 0.2 0.2
BaseEpsilon = 100
[LinearSolver]
Type = bicgstab
Verbosity = 2
MaxIterations = 500
ResidualReduction = 1e-8
GMResRestart = 30
SymmetrizeDirichlet = true
DirectSolverForVelocity = false
AmgForPressure = true
[LinearSolver.Preconditioner]
Verbosity = 0
AmgDefaultAggregationDimension = 2
AmgMinAggregateSize = 2
AmgMaxAggregateSize = 4
AmgPreSmoothingSteps = 2
AmgPostSmoothingSteps = 2
AmgAdditive = false
AmgGamma = 1 # 1: V-cycle 2: W-cycle
AmgCriterionSymmetric = true
MassMatrixWeight = 2.0
[Newton]
MinSteps = 1
EnableAbsoluteResidualCriterion = true
MaxAbsoluteResidual = 1e-7
// -*- 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_NAVIERSTOKES_3D_CHANNEL_PROBLEM_HH
#define DUMUX_TEST_FREEFLOW_NAVIERSTOKES_3D_CHANNEL_PROBLEM_HH
#include <cmath>
#include <numeric>
#include <dune/common/float_cmp.hh>
#include <dune/common/fmatrix.hh>
#include <dumux/common/math.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/io/vtk/function.hh>
#include <dumux/io/grid/griddata.hh>
#include <dumux/freeflow/navierstokes/problem.hh>
namespace Dumux {
/*!
* \brief Test problem for the (Navier-) Stokes model in a 3D channel
*
* Benchmark case from
* Turek, Schaefer et a (1996) Benchmark computations of laminar flow around cylinder.
* Flow Simulation with High-Performance Computers II,
* Notes on Numerical Fluid Mechanics 52, 547-566, Vieweg 1996
* https://doi.org/10.1007/978-3-322-89849-4_39
*/
template <class TypeTag>
class DFGChannelTestProblem : public NavierStokesProblem<TypeTag>
{
using ParentType = NavierStokesProblem<TypeTag>;
using BoundaryTypes = typename ParentType::BoundaryTypes;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using BoundaryFluxes = typename ParentType::BoundaryFluxes;
using DirichletValues = typename ParentType::DirichletValues;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
static constexpr int dim = GridGeometry::GridView::dimension;
static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
public:
DFGChannelTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
: ParentType(gridGeometry, couplingManager, ParentType::isMomentumProblem() ? "Momentum" : "Mass")
{
density_ = getParam<Scalar>("Component.LiquidDensity");
viscosity_ = getParam<Scalar>("Component.LiquidDynamicViscosity");
for (int i = 0; i < dimWorld; ++i)
domainSize_[i] = gridGeometry->bBoxMax()[i] - gridGeometry->bBoxMin()[i];
}
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
* \param element The finite element
* \param scvf The sub control volume face
*/
BoundaryTypes boundaryTypes(const Element& element,
const SubControlVolumeFace& scvf) const
{
BoundaryTypes values;
if constexpr (ParentType::isMomentumProblem())
{
if (isOutlet_(scvf))
values.setAllNeumann();
else
values.setAllDirichlet();
}
else
values.setNeumann(Indices::conti0EqIdx);
return values;
}
/*!
* \brief Evaluates the boundary conditions for a Dirichlet boundary face
*/
DirichletValues dirichletAtPos(const GlobalPosition& globalPos) const
{
if (globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_)
return inflowVelocity_(globalPos[1]);
else
return DirichletValues(0.0);
}
/*!
* \brief Evaluates the boundary conditions for a Neumann boundary face
*
* \param element The element for which the Neumann boundary condition is set
* \param fvGeometry The fvGeometry
* \param elemVolVars The element volume variables
* \param elemFaceVars The element face variables
* \param scvf The boundary sub control volume face
*/
template<class ElementVolumeVariables, class ElementFluxVariablesCache>
BoundaryFluxes neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
BoundaryFluxes values(0.0);
if constexpr (ParentType::isMomentumProblem())
{
if (this->enableInertiaTerms())
{
if (isOutlet_(scvf))
{
// advective term: vv*n
const auto v = elemVolVars[scvf.insideScvIdx()].velocity();
values.axpy(density_*(v*scvf.unitOuterNormal()), v);
}
}
}
else
{
values[Indices::conti0EqIdx]
= this->faceVelocity(element, fvGeometry, scvf)
* density_ * scvf.unitOuterNormal();
}
return values;
}
template<class ScvOrScvf>
Scalar density(const Element&,
const FVElementGeometry&,
const ScvOrScvf&,
const bool isPreviousTimeStep = false) const
{
return density_;
}
template<class ScvOrScvf>
Scalar effectiveViscosity(const Element&,
const FVElementGeometry&,
const ScvOrScvf&) const
{
return viscosity_;
}
//! Computes pressure difference benchmark indicator
//! (only works for correct domain size 2.2 x 0.41 and boundary conditions, Re=20)
template<class SolutionVector, class P = ParentType, typename std::enable_if_t<!P::isMomentumProblem(), int> = 0>
Scalar evalPressureDifference(const SolutionVector& p) const
{
const auto& tree = this->gridGeometry().boundingBoxTree();
GlobalPosition evalPoint1({0.15, 0.2});
GlobalPosition evalPoint2({0.25, 0.2});
const Scalar stepSize = 0.001;
auto entities1 = intersectingEntities(evalPoint1, tree);
while (entities1.empty())
{
evalPoint1[0] -= stepSize;
entities1 = intersectingEntities(evalPoint1, tree);
}
auto entities2 = intersectingEntities(evalPoint2, tree);
while (entities2.empty())
{
evalPoint2[0] += stepSize;
entities2 = intersectingEntities(evalPoint2, tree);
}
const auto p1 = p[entities1[0]][0];
const auto p2 = p[entities2[0]][0];
return p1 - p2;
}
//! Computes drag and lift coefficient benchmark indicator
//! (only works for correct domain size 2.2 x 0.41 and boundary conditions, Re=20)
template<class SolutionVector, class GridVariables, class P = ParentType, typename std::enable_if_t<P::isMomentumProblem(), int> = 0>
auto evalDragAndLiftCoefficient(const GridVariables& gridVariables, const SolutionVector& v) const
{
auto fvGeometry = localView(this->gridGeometry());
auto elemVolVars = localView(gridVariables.curGridVolVars());
auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());