Commit a133773f authored by Roman Winter's avatar Roman Winter Committed by Timo Koch
Browse files

[port] channel/1d to new staggered - failing test

parent 3c8c8b8e
Pipeline #9436 passed with stages
......@@ -24,24 +24,25 @@
#include <config.h>
#include <ctime>
#include <iostream>
#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.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/traits.hh>
#include <dumux/multidomain/staggeredfreeflow/couplingmanager.hh>
#include <dumux/multidomain/newtonsolver.hh>
#include <dumux/freeflow/navierstokes/velocityoutput.hh>
#include <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
#include <test/freeflow/navierstokes/errors.hh>
......@@ -52,7 +53,8 @@ int main(int argc, char** argv)
using namespace Dumux;
// define the type tag for this problem
using TypeTag = Properties::TTag::NavierStokesAnalytic;
using MomentumTypeTag = Properties::TTag::NavierStokesAnalyticMomentum;
using MassTypeTag = Properties::TTag::NavierStokesAnalyticMass;
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
......@@ -65,7 +67,7 @@ int main(int argc, char** argv)
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<GetPropType<MomentumTypeTag, Properties::Grid>> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
......@@ -76,65 +78,84 @@ 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);
// the problem (boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
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 (initial and 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
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]);
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);
// 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");
NavierStokesTest::AnalyticalSolutionVectors analyticalSolVectors(momentumProblem, massProblem);
vtkWriter.addField(analyticalSolVectors.analyticalPressureSolution(), "pressureExact");
vtkWriter.addField(analyticalSolVectors.analyticalVelocitySolution(), "velocityExact");
vtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<MassGridVariables>>());
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);
using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
// the discrete L2 and Linfity errors
const bool printErrors = getParam<bool>("Problem.PrintErrors", false);
const bool printConvergenceTestFile = getParam<bool>("Problem.PrintConvergenceTestFile", false);
// linearize & solve
Dune::Timer timer;
nonLinearSolver.solve(x);
// print discrete L2 and Linfity errors
const bool printErrors = getParam<bool>("Problem.PrintErrors", false);
const bool printConvergenceTestFile = getParam<bool>("Problem.PrintConvergenceTestFile", false);
if (printErrors || printConvergenceTestFile)
{
NavierStokesErrors errors(problem, x);
NavierStokesErrorCSVWriter(
problem, std::to_string(x[GridGeometry::cellCenterIdx()].size())
).printErrors(errors);
NavierStokesTest::Errors errors(momentumProblem, massProblem, x);
NavierStokesTest::ErrorCSVWriter errorCSVWriter(momentumProblem, massProblem);
errors.update(x);
errorCSVWriter.printErrors(errors);
if (printConvergenceTestFile)
convergenceTestAppendErrors(problem, errors);
convergenceTestAppendErrors(momentumProblem, massProblem, errors);
}
// write vtk output
......
......@@ -56,18 +56,22 @@ class NavierStokesAnalyticProblem : public NavierStokesProblem<TypeTag>
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using FVElementGeometry = typename GridGeometry::LocalView;
static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
using DimVector = GlobalPosition;
using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
public:
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
NavierStokesAnalyticProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
NavierStokesAnalyticProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
: ParentType(gridGeometry, couplingManager)
{
density_ = getParam<Scalar>("Component.LiquidDensity");
kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity");
......@@ -90,36 +94,41 @@ public:
{
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)
if constexpr (!ParentType::isMomentumProblem())
{
// mass balance - term div(rho*v)
for (unsigned int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
{
// inertia term
if (this->enableInertiaTerms())
source[Indices::velocity(velIdx)] += density_ * dv2dx(globalPos)[velIdx][dimIdx];
// viscous term (molecular)
source[Indices::velocity(velIdx)] -= density_ * kinematicViscosity_* dvdx2(globalPos)[velIdx][dimIdx];
static const bool enableUnsymmetrizedVelocityGradient = getParam<bool>("FreeFlow.EnableUnsymmetrizedVelocityGradient", false);
if (!enableUnsymmetrizedVelocityGradient)
source[Indices::velocity(velIdx)] -= density_ * kinematicViscosity_* dvdx2(globalPos)[dimIdx][velIdx];
source[Indices::conti0EqIdx] += dvdx(globalPos)[dimIdx][dimIdx];
}
// pressure term
source[Indices::velocity(velIdx)] += dpdx(globalPos)[velIdx];
// gravity term
static const bool enableGravity = getParam<bool>("Problem.EnableGravity");
if (enableGravity)
source[Indices::conti0EqIdx] *= density_;
}
else
{
// momentum balance
for (unsigned int velIdx = 0; velIdx < dimWorld; ++velIdx)
{
source[Indices::velocity(velIdx)] -= density_ * this->gravity()[velIdx];
for (unsigned int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
{
// inertia term
if (this->enableInertiaTerms())
source[Indices::velocity(velIdx)] += density_ * dv2dx(globalPos)[velIdx][dimIdx];
// viscous term (molecular)
source[Indices::velocity(velIdx)] -= density_ * kinematicViscosity_* dvdx2(globalPos)[velIdx][dimIdx];
static const bool enableUnsymmetrizedVelocityGradient = getParam<bool>("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 const bool enableGravity = getParam<bool>("Problem.EnableGravity");
if (enableGravity)
{
source[Indices::velocity(velIdx)] -= density_ * this->gravity()[velIdx];
}
}
}
......@@ -141,33 +150,17 @@ public:
{
BoundaryTypes values;
// set Dirichlet values for the velocity everywhere
values.setDirichlet(Indices::momentumXBalanceIdx);
if constexpr (ParentType::isMomentumProblem())
{
// set Dirichlet values for the velocity everywhere
values.setDirichlet(Indices::momentumXBalanceIdx);
}
else
values.setNeumann(Indices::pressureIdx);
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 GridGeometry::LocalView& fvGeometry,
const typename GridGeometry::SubControlVolume& scv,
int pvIdx) const
{
// set a fixed pressure in all cells at the boundary
for (const auto& scvf : scvfs(fvGeometry))
if (scvf.boundary())
return true;
return false;
}
/*!
* \brief Returns Dirichlet boundary values at a given position
*
......@@ -188,8 +181,12 @@ public:
PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
{
PrimaryVariables values;
values[Indices::pressureIdx] = p(globalPos);
values[Indices::velocityXIdx] = v(globalPos);
if constexpr (ParentType::isMomentumProblem())
values[Indices::velocityXIdx] = v(globalPos);
else
values[Indices::pressureIdx] = p(globalPos);
return values;
}
......@@ -244,6 +241,51 @@ public:
return dpdx;
}
//! 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;
auto fvGeometry = localView(this->gridGeometry());
fvGeometry.bindElement(element);
bool onBoundary = false;
for (const auto& scvf : scvfs(fvGeometry))
onBoundary = std::max(onBoundary, scvf.boundary());
if (onBoundary)
values.set(0);
// TODO: only use one cell or pass fvGeometry to hasInternalDirichletConstraint
// if (scv.dofIndex() == 0)
// values.set(0);
// 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]); }
// \}
/*!
......
......@@ -25,19 +25,26 @@
#define DUMUX_DONEA_TEST_PROPERTIES_HH
#include <dune/grid/yaspgrid.hh>
#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/freeflow/navierstokes/momentum/model.hh>
#include <dumux/freeflow/navierstokes/mass/1p/model.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/multidomain/staggeredfreeflow/couplingmanager.hh>
#include "problem.hh"
namespace Dumux::Properties {
// Create new type tags
namespace TTag {
struct NavierStokesAnalytic { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
struct NavierStokesAnalytic {};
struct NavierStokesAnalyticMomentum { using InheritsFrom = std::tuple<NavierStokesAnalytic, NavierStokesMomentum, FaceCenteredStaggeredModel>; };
struct NavierStokesAnalyticMass { using InheritsFrom = std::tuple<NavierStokesAnalytic, NavierStokesMassOneP, CCTpfaModel>; };
} // end namespace TTag
// the fluid system
......@@ -66,6 +73,13 @@ struct EnableGridVolumeVariablesCache<TypeTag, TTag::NavierStokesAnalytic> { sta
template<class TypeTag>
struct NormalizePressure<TypeTag, TTag::NavierStokesAnalytic> { static constexpr bool value = false; };
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::NavierStokesAnalytic>
{
using Traits = MultiDomainTraits<TTag::NavierStokesAnalyticMomentum, TTag::NavierStokesAnalyticMass>;
using type = StaggeredFreeFlowCouplingManager<Traits>;
};
} // end namespace Dumux::Properties
#endif
Supports Markdown
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