Commit fc783410 authored by Yue Wang's avatar Yue Wang Committed by Timo Koch
Browse files

[test][ff][navierstokes] Port Kovazsny test to new staggered

parent 071b1a8e
Pipeline #9308 failed with stages
......@@ -29,24 +29,26 @@
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/istl/io.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/io/grid/gridmanager_sub.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
#include <test/freeflow/navierstokes/errors.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/traits.hh>
#include <dumux/multidomain/newtonsolver.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/freeflow/navierstokes/velocityoutput.hh>
#include "properties.hh"
int main(int argc, char** argv)
......@@ -54,7 +56,8 @@ int main(int argc, char** argv)
using namespace Dumux;
// define the type tag for this problem
using TypeTag = Properties::TTag::KovasznayTest;
using MomentumTypeTag = Properties::TTag::KovasznayTestMomentum;
using MassTypeTag = Properties::TTag::KovasznayTestMass;
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
......@@ -67,7 +70,7 @@ int main(int argc, char** argv)
Parameters::init(argc, argv);
// create a grid
using Grid = GetPropType<TypeTag, Properties::Grid>;
using Grid = GetPropType<MomentumTypeTag, Properties::Grid>;
Dumux::GridManager<Grid> gridManager;
#if HAVE_DUNE_SUBGRID
......@@ -106,47 +109,71 @@ int main(int argc, char** argv)
const auto& leafGridView = gridManager.grid().leafGridView();
// create the finite volume grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
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 problem (boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
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
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
constexpr auto momentumIdx = CouplingManager::freeFlowMomentumIndex;
constexpr auto massIdx = CouplingManager::freeFlowMassIndex;
using SolutionVector = typename Traits::SolutionVector;
SolutionVector x;
x[GridGeometry::cellCenterIdx()].resize(gridGeometry->numCellCenterDofs());
x[GridGeometry::faceIdx()].resize(gridGeometry->numFaceDofs());
x[momentumIdx].resize(momentumGridGeometry->numDofs());
x[massIdx].resize(massGridGeometry->numDofs());
// the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
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);
couplingManager->init(momentumProblem, massProblem, std::make_tuple(momentumGridVariables, massGridVariables), x);
massGridVariables->init(x[massIdx]);
momentumGridVariables->init(x[momentumIdx]);
// intialize the vtk output module
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
using IOFields = GetPropType<MassTypeTag, Properties::IOFields>;
VtkOutputModule vtkWriter(*massGridVariables, x[massIdx], massProblem->name());
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
NavierStokesAnalyticalSolutionVectors analyticalSolVectors(problem);
vtkWriter.addField(analyticalSolVectors.getAnalyticalPressureSolution(), "pressureExact");
vtkWriter.addField(analyticalSolVectors.getAnalyticalVelocitySolution(), "velocityExact");
vtkWriter.addFaceField(analyticalSolVectors.getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
using AnalyticalSolutionVectors = NavierStokesTest::AnalyticalSolutionVectors<MomentumProblem,MassProblem>;
AnalyticalSolutionVectors analyticalSolutionVectors(momentumProblem, massProblem);
vtkWriter.write(0.0);
vtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<MassGridVariables>>());
const auto exactPressure = analyticalSolutionVectors.analyticalPressureSolution();
const auto exactVelocity = analyticalSolutionVectors.analyticalVelocitySolution();
vtkWriter.addField(exactPressure, "pressureExact");
vtkWriter.addField(exactVelocity, "velocityExact");
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);
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);
// the linear solver
using LinearSolver = Dumux::UzawaBiCGSTABBackend<LinearSolverTraits<GridGeometry>>; // if rho \neq 1, use UMFPack rather than BIGCStab
using LinearSolver = Dumux::UzawaBiCGSTABBackend<LinearSolverTraits<MassGridGeometry>>; // if rho \neq 1, use UMFPack rather than BIGCStab
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
// linearize & solve
Dune::Timer timer;
......@@ -158,13 +185,15 @@ int main(int argc, char** argv)
if (printErrors || printConvergenceTestFile)
{
NavierStokesErrors errors(problem, x);
NavierStokesErrorCSVWriter(
problem, std::to_string(x[GridGeometry::cellCenterIdx()].size())
NavierStokesTest::Errors errors(
momentumProblem, massProblem, x
);
NavierStokesTest::ErrorCSVWriter(
momentumProblem, massProblem
).printErrors(errors);
if (printConvergenceTestFile)
convergenceTestAppendErrors(problem, errors);
convergenceTestAppendErrors(momentumProblem, massProblem, errors);
}
// write vtk output
......
......@@ -47,12 +47,14 @@ class KovasznayTestProblem : public NavierStokesProblem<TypeTag>
{
using ParentType = NavierStokesProblem<TypeTag>;
using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
using BoundaryTypes = typename ParentType::BoundaryTypes;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using PrimaryVariables = typename ParentType::PrimaryVariables;
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
......@@ -60,15 +62,18 @@ class KovasznayTestProblem : public NavierStokesProblem<TypeTag>
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
public:
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
KovasznayTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
KovasznayTestProblem(std::shared_ptr<const GridGeometry> gridGeometry,
std::shared_ptr<CouplingManager> couplingManager)
: ParentType(gridGeometry, couplingManager)
{
std::cout<< "upwindSchemeOrder is: " << GridGeometry::upwindStencilOrder() << "\n";
std::cout<< "upwindSchemeOrder is: " << upwindSchemeOrder << "\n";
rho_ = getParam<Scalar>("Component.LiquidDensity", 1.0);
kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
Scalar reynoldsNumber = 1.0 / kinematicViscosity_;
......@@ -99,41 +104,17 @@ public:
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
// set Dirichlet values for the velocity everywhere
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
if constexpr (ParentType::isMomentumProblem())
{
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
}
else
values.setAllNeumann();
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 FVElementGeometry& fvGeometry,
const SubControlVolume& scv,
int pvIdx) const
{
// set fixed pressure in all cells at the left boundary
auto isAtLeftBoundary = [&](const FVElementGeometry& fvGeometry)
{
if (fvGeometry.hasBoundaryScvf())
{
for (const auto& scvf : scvfs(fvGeometry))
if (scvf.boundary() && scvf.center()[0] < this->gridGeometry().bBoxMin()[0] + eps_)
return true;
}
return false;
};
return (isAtLeftBoundary(fvGeometry) && pvIdx == Indices::pressureIdx);
}
/*!
* \brief Returns Dirichlet boundary values at a given position.
*
......@@ -144,7 +125,30 @@ public:
// use the values of the analytical solution
return analyticalSolution(globalPos);
}
/*!
* \brief Evaluates the boundary conditions for a Neumann control volume.
*
* \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>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
if constexpr (!ParentType::isMomentumProblem())
{
const auto insideDensity = elemVolVars[scvf.insideScvIdx()].density();
values[Indices::conti0EqIdx] = insideDensity * (this->faceVelocity(element, fvGeometry, scvf) * scvf.unitOuterNormal());
}
return values;
}
/*!
* \brief Returns the analytical solution of the problem at a given position.
*
......@@ -153,13 +157,19 @@ public:
*/
PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
{
Scalar x = globalPos[0];
Scalar y = globalPos[1];
const Scalar x = globalPos[0];
PrimaryVariables values;
values[Indices::pressureIdx] = rho_ * 0.5 * (1.0 - std::exp(2.0 * lambda_ * x));
values[Indices::velocityXIdx] = 1.0 - std::exp(lambda_ * x) * std::cos(2.0 * M_PI * y);
values[Indices::velocityYIdx] = 0.5 * lambda_ / M_PI * std::exp(lambda_ * x) * std::sin(2.0 * M_PI * y);
if constexpr (ParentType::isMomentumProblem())
{
const Scalar y = globalPos[1];
values[Indices::velocityXIdx] = 1.0 - std::exp(lambda_ * x) * std::cos(2.0 * M_PI * y);
values[Indices::velocityYIdx] = 0.5 * lambda_ / M_PI * std::exp(lambda_ * x) * std::sin(2.0 * M_PI * y);
}
else
{ values[Indices::pressureIdx] = rho_ *0.5 * (1.0 - std::exp(2.0 * lambda_ * x));
}
return values;
}
......@@ -178,14 +188,53 @@ public:
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values;
values[Indices::pressureIdx] = 0.0;
values[Indices::velocityXIdx] = 0.0;
values[Indices::velocityYIdx] = 0.0;
return PrimaryVariables(0.0);
}
//! Enable internal Dirichlet constraints
static constexpr bool enableInternalDirichletConstraints()
{ return !ParentType::isMomentumProblem(); }
/*!
* \brief Tag a degree of freedom to carry internal Dirichlet constraints.
* If true is returned for a dof, the equation for this dof is replaced
* by the constraint that its primary variable values must match the
* user-defined values obtained from the function internalDirichlet(),
* which must be defined in the problem.
*
* \param element The finite element
* \param scv The sub-control volume
*/
std::bitset<PrimaryVariables::dimension> hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const
{
std::bitset<PrimaryVariables::dimension> values;
const auto fvGeometry = localView(this->gridGeometry()).bindElement(element);
const bool isAtLeftBoundary = [&]{
if (fvGeometry.hasBoundaryScvf())
for (const auto& scvf : scvfs(fvGeometry))
if (scvf.boundary() && scvf.center()[0] < this->gridGeometry().bBoxMin()[0] + eps_)
return true;
return false;
}();
if (isAtLeftBoundary)
values.set(0);
// TODO: only use one cell or pass fvGeometry to hasInternalDirichletConstraint
// the pure Neumann problem is only defined up to a constant
// we create a well-posed problem by fixing the pressure at one dof
return values;
}
/*!
* \brief Define the values of internal Dirichlet constraints for a degree of freedom.
* \param element The finite element
* \param scv The sub-control volume
*/
PrimaryVariables internalDirichlet(const Element& element, const SubControlVolume& scv) const
{ return PrimaryVariables(analyticalSolution(scv.center())[Indices::pressureIdx]); }
private:
static constexpr Scalar eps_=1e-6;
......
......@@ -34,18 +34,25 @@
#include <dune/subgrid/subgrid.hh>
#endif
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/freeflow/navierstokes/model.hh>
#include <dumux/discretization/fcstaggered.hh>
#include <dumux/discretization/cctpfa.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/freeflow/navierstokes/momentum/model.hh>
#include <dumux/freeflow/navierstokes/mass/1p/model.hh>
#include <dumux/multidomain/staggeredfreeflow/couplingmanager.hh>
#include "problem.hh"
namespace Dumux::Properties {
// Create new type tags
namespace TTag {
struct KovasznayTest { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
struct KovasznayTest {};
struct KovasznayTestMomentum { using InheritsFrom = std::tuple<KovasznayTest, NavierStokesMomentum, FaceCenteredStaggeredModel>; };
struct KovasznayTestMass { using InheritsFrom = std::tuple<KovasznayTest, NavierStokesMassOneP, CCTpfaModel>; };
} // end namespace TTag
// the fluid system
......@@ -69,6 +76,16 @@ struct Grid<TypeTag, TTag::KovasznayTest>
#endif
};
// set the coupling manager property
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::KovasznayTest>
{
private:
using Traits = MultiDomainTraits<TTag::KovasznayTestMomentum, TTag::KovasznayTestMass>;
public:
using type = StaggeredFreeFlowCouplingManager<Traits>;
};
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::KovasznayTest> { using type = Dumux::KovasznayTestProblem<TypeTag> ; };
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment