From 239af28c82549a30d1324822d806bd7611de45cc Mon Sep 17 00:00:00 2001
From: Edward 'Ned' Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Wed, 26 Apr 2023 12:06:12 +0200
Subject: [PATCH] [test][navierstokesnc] Port channel test
---
.../navierstokesnc/channel/CMakeLists.txt | 15 +-
test/freeflow/navierstokesnc/channel/main.cc | 83 ++++--
.../navierstokesnc/channel/problem.hh | 278 ++++++++++--------
.../navierstokesnc/channel/properties.hh | 17 +-
.../navierstokesnc/maxwellstefan/main.cc | 4 +-
.../navierstokesnc/maxwellstefan/problem.hh | 10 +-
.../maxwellstefan/properties.hh | 33 ++-
7 files changed, 262 insertions(+), 178 deletions(-)
diff --git a/test/freeflow/navierstokesnc/channel/CMakeLists.txt b/test/freeflow/navierstokesnc/channel/CMakeLists.txt
index 0e68ff1588..34318acecc 100644
--- a/test/freeflow/navierstokesnc/channel/CMakeLists.txt
+++ b/test/freeflow/navierstokesnc/channel/CMakeLists.txt
@@ -4,11 +4,11 @@
dune_symlink_to_source_files(FILES "params_advection.input" "params_diffusion.input" "params_advectionni.input" "params_diffusionni.input")
add_executable(test_ff_stokes2c_mass EXCLUDE_FROM_ALL main.cc)
-target_compile_definitions(test_ff_stokes2c_mass PUBLIC "USE_MASS=1")
+target_compile_definitions(test_ff_stokes2c_mass PUBLIC "USE_MOLES=0")
dumux_add_test(NAME test_ff_stokes2c_diffusion_mass
TARGET test_ff_stokes2c_mass
- LABELS freeflow navierstokes
+ LABELS freeflow navierstokes navierstokesnc
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
@@ -19,10 +19,11 @@ dumux_add_test(NAME test_ff_stokes2c_diffusion_mass
--zeroThreshold {"velocity_liq \(m/s\)":1e-22})
add_executable(test_ff_stokes2c EXCLUDE_FROM_ALL main.cc)
+target_compile_definitions(test_ff_stokes2c PUBLIC "USE_MOLES=1")
dumux_add_test(NAME test_ff_stokes2c_diffusion_mole
TARGET test_ff_stokes2c
- LABELS freeflow navierstokes
+ LABELS freeflow navierstokes navierstokesnc
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
@@ -34,7 +35,7 @@ dumux_add_test(NAME test_ff_stokes2c_diffusion_mole
dumux_add_test(NAME test_ff_stokes2c_advection
TARGET test_ff_stokes2c
- LABELS freeflow navierstokes
+ LABELS freeflow navierstokes navierstokesnc
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
@@ -45,7 +46,7 @@ dumux_add_test(NAME test_ff_stokes2c_advection
dumux_add_test(NAME test_ff_stokes2c_advection_nocaching
SOURCES main.cc
- LABELS freeflow navierstokes
+ LABELS freeflow navierstokes navierstokesnc
CMAKE_GUARD HAVE_UMFPACK
COMPILE_DEFINITIONS ENABLECACHING=0
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
@@ -60,7 +61,7 @@ target_compile_definitions(test_ff_stokes2cni PUBLIC "NONISOTHERMAL=1")
dumux_add_test(NAME test_ff_stokes2cni_advection
TARGET test_ff_stokes2cni
- LABELS freeflow navierstokes
+ LABELS freeflow navierstokes navierstokesnc
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
@@ -71,7 +72,7 @@ dumux_add_test(NAME test_ff_stokes2cni_advection
dumux_add_test(NAME test_ff_stokes2cni_diffusion
TARGET test_ff_stokes2cni
- LABELS freeflow navierstokes
+ LABELS freeflow navierstokes navierstokesnc
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
diff --git a/test/freeflow/navierstokesnc/channel/main.cc b/test/freeflow/navierstokesnc/channel/main.cc
index 066e4e4a15..f313120b59 100644
--- a/test/freeflow/navierstokesnc/channel/main.cc
+++ b/test/freeflow/navierstokesnc/channel/main.cc
@@ -17,11 +17,12 @@
#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_yasp.hh>
#include <dumux/linear/istlsolvers.hh>
#include <dumux/linear/linearsolvertraits.hh>
@@ -34,6 +35,11 @@
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/freeflow/navierstokes/velocityoutput.hh>
+
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/traits.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+
#include "properties.hh"
int main(int argc, char** argv)
@@ -41,7 +47,8 @@ int main(int argc, char** argv)
using namespace Dumux;
// define the type tag for this problem
- using TypeTag = Properties::TTag::ChannelNCTest;
+ using MomentumTypeTag = Properties::TTag::ChannelNCTestMomentum;
+ using MassTypeTag = Properties::TTag::ChannelNCTestMass;
// maybe initialize MPI and/or multithreading backend
initialize(argc, argv);
@@ -55,7 +62,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<MassTypeTag, Properties::Grid>> gridManager;
+ GridManager<GetPropType<MomentumTypeTag, Properties::Grid>> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
@@ -66,8 +73,10 @@ 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 CouplingManager = GetPropType<MomentumTypeTag, Properties::CouplingManager>;
@@ -80,15 +89,13 @@ int main(int argc, char** argv)
auto massProblem = std::make_shared<MassProblem>(massGridGeometry, couplingManager);
// get some time loop parameters
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using Scalar = GetPropType<MomentumTypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
- // instantiate time loop
- auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, dt, tEnd);
- timeLoop->setMaxTimeStepSize(maxDt);
- problem->setTimeLoop(timeLoop);
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = getParam<Scalar>("Restart.Time", 0);
// the solution vector
constexpr auto momentumIdx = CouplingManager::freeFlowMomentumIndex;
@@ -96,34 +103,51 @@ int main(int argc, char** argv)
using Traits = MultiDomainTraits<MomentumTypeTag, MassTypeTag>;
using SolutionVector = typename Traits::SolutionVector;
SolutionVector x;
- x[GridGeometry::cellCenterIdx()].resize(gridGeometry->numCellCenterDofs());
- x[GridGeometry::faceIdx()].resize(gridGeometry->numFaceDofs());
- problem->applyInitialSolution(x);
+ x[momentumIdx].resize(momentumGridGeometry->numDofs());
+ x[massIdx].resize(massGridGeometry->numDofs());
+ momentumProblem->applyInitialSolution(x[momentumIdx]);
+ massProblem->applyInitialSolution(x[massIdx]);
auto xOld = x;
+ // instantiate time loop
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ massProblem->setTimeLoop(timeLoop);
+ momentumProblem->setTimeLoop(timeLoop);
+
// 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);
- // initialize the vtk output module
- using IOFields = GetPropType<TypeTag, Properties::IOFields>;
- StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
- IOFields::initOutputModule(vtkWriter); // Add model specific output fields
- vtkWriter.addField(problem->getDeltaP(), "deltaP");
- vtkWriter.write(0.0);
+ using MassGridVariables = GetPropType<MassTypeTag, Properties::GridVariables>;
+ auto massGridVariables = std::make_shared<MassGridVariables>(massProblem, massGridGeometry);
- // the assembler with time loop for instationary problem
- using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
- auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
+ couplingManager->init(momentumProblem, massProblem, std::make_tuple(momentumGridVariables, massGridVariables), x, xOld);
+ momentumGridVariables->init(x[momentumIdx]);
+ massGridVariables->init(x[massIdx]);
+ // 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.addVolumeVariable([](const auto& volVars){ return volVars.pressure() - 1.1e5; }, "deltaP");
+ vtkWriter.write(0);
+
+ 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,
+ timeLoop, xOld);
// the linear solver
using LinearSolver = Dumux::UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
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);
// time loop
timeLoop->start(); do
@@ -133,13 +157,12 @@ int main(int argc, char** argv)
// make the new solution the old solution
xOld = x;
- gridVariables->advanceTimeStep();
+ momentumGridVariables->advanceTimeStep();
+ massGridVariables->advanceTimeStep();
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
- problem->calculateDeltaP(*gridVariables, x);
-
// write vtk output
vtkWriter.write(timeLoop->time());
diff --git a/test/freeflow/navierstokesnc/channel/problem.hh b/test/freeflow/navierstokesnc/channel/problem.hh
index 06c4699b3f..53638e84bb 100644
--- a/test/freeflow/navierstokesnc/channel/problem.hh
+++ b/test/freeflow/navierstokesnc/channel/problem.hh
@@ -13,8 +13,8 @@
#ifndef DUMUX_CHANNEL_NC_TEST_PROBLEM_HH
#define DUMUX_CHANNEL_NC_TEST_PROBLEM_HH
-#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
+#include <dumux/common/properties.hh>
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
@@ -32,202 +32,250 @@ namespace Dumux {
* a fixed value on the right boundary. The top and bottom boundaries
* represent solid walls with no-slip/no-flow conditions.
*/
-template <class TypeTag>
-class ChannelNCTestProblem : public NavierStokesStaggeredProblem<TypeTag>
+template <class TypeTag, class BaseProblem>
+class ChannelNCTestProblem : public BaseProblem
{
- using ParentType = NavierStokesStaggeredProblem<TypeTag>;
-
- using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
- using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using ParentType = BaseProblem;
+ 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 PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
- using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+ using InitialValues = typename ParentType::InitialValues;
+ using Sources = typename ParentType::Sources;
+ using DirichletValues = typename ParentType::DirichletValues;
+ using BoundaryFluxes = typename ParentType::BoundaryFluxes;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ 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>;
using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
+ static constexpr bool useMoles = ModelTraits::useMoles();
+ static constexpr bool isNonisothermal = ModelTraits::enableEnergyBalance();
static constexpr auto compIdx = 1;
- static constexpr auto transportCompIdx = Indices::conti0EqIdx + compIdx;
- static constexpr auto transportEqIdx = Indices::conti0EqIdx + compIdx;
public:
- ChannelNCTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
- : ParentType(gridGeometry), eps_(1e-6)
+ ChannelNCTestProblem(std::shared_ptr<const GridGeometry> gridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(gridGeometry, couplingManager)
+ , eps_(1e-6)
{
- inletVelocity_ = getParam<Scalar>("Problem.InletVelocity");
FluidSystem::init();
- deltaP_.resize(this->gridGeometry().numCellCenterDofs());
+ inletVelocity_ = getParam<Scalar>("Problem.InletVelocity");
}
- /*!
- * \name Boundary conditions
+ /*!
+ * \brief Returns a reference pressure at a given sub control volume face.
+ * This pressure is subtracted from the actual pressure for the momentum balance
+ * which potentially helps to improve numerical accuracy by avoiding issues related do floating point arithmetic.
*/
- // \{
+ Scalar referencePressure(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf) const
+ { return 1.1e5; }
- /*!
+ /*!
* \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 boundaryTypesAtPos(const GlobalPosition& globalPos) const
{
BoundaryTypes values;
- if(isInlet_(globalPos))
+ if constexpr (ParentType::isMomentumProblem())
{
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
- values.setDirichlet(transportCompIdx);
-#if NONISOTHERMAL
- values.setDirichlet(Indices::temperatureIdx);
-#endif
- }
- else if(isOutlet_(globalPos))
- {
- values.setDirichlet(Indices::pressureIdx);
- values.setOutflow(transportEqIdx);
-#if NONISOTHERMAL
- values.setOutflow(Indices::energyEqIdx);
-#endif
+
+ if (isOutlet_(globalPos))
+ values.setAllNeumann();
}
else
{
- // set Dirichlet values for the velocity everywhere
- values.setDirichlet(Indices::velocityXIdx);
- values.setDirichlet(Indices::velocityYIdx);
- values.setNeumann(Indices::conti0EqIdx);
- values.setNeumann(transportEqIdx);
-#if NONISOTHERMAL
- values.setNeumann(Indices::energyEqIdx);
-#endif
+ values.setAllNeumann();
+
+ if (isInlet_(globalPos))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ values.setDirichlet(Indices::conti0EqIdx + compIdx);
+ if constexpr (isNonisothermal)
+ values.setDirichlet(Indices::temperatureIdx);
+ }
}
return values;
}
- /*!
+ /*!
* \brief Evaluates the boundary conditions for a Dirichlet control volume.
*
* \param globalPos The center of the finite volume which ought to be set.
*/
- PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+ DirichletValues dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
{
- PrimaryVariables values = initialAtPos(globalPos);
+ const auto& globalPos = scvf.ipGlobal();
+ DirichletValues values = initialAtPos(globalPos);
- // give the system some time so that the pressure can equilibrate, then start the injection of the tracer
- if(isInlet_(globalPos))
+ if constexpr (!ParentType::isMomentumProblem())
{
- if(time() >= 10.0 || inletVelocity_ < eps_)
+ // give the system some time so that the pressure can equilibrate, then start the injection of the tracer
+ if (isInlet_(globalPos))
{
- Scalar moleFracTransportedComp = 1e-3;
-#if USE_MASS
- Scalar averageMolarMassPhase = moleFracTransportedComp * FluidSystem::molarMass(compIdx)
- + (1. - moleFracTransportedComp) * FluidSystem::molarMass(1-compIdx);
- values[transportCompIdx] = moleFracTransportedComp * FluidSystem::molarMass(compIdx)
- / averageMolarMassPhase;
-#else
- values[transportCompIdx] = moleFracTransportedComp;
-#endif
-#if NONISOTHERMAL
- values[Indices::temperatureIdx] = 293.15;
-#endif
+ values[Indices::pressureIdx] = this->couplingManager().cellPressure(element, scvf);
+
+ if (time() >= 10.0 || inletVelocity_ < eps_)
+ {
+ Scalar moleFracTransportedComp = 1e-3;
+ if constexpr (useMoles)
+ values[Indices::conti0EqIdx+compIdx] = moleFracTransportedComp;
+ else
+ {
+ Scalar averageMolarMassPhase = moleFracTransportedComp * FluidSystem::molarMass(compIdx)
+ + (1. - moleFracTransportedComp) * FluidSystem::molarMass(1-compIdx);
+ values[Indices::conti0EqIdx+compIdx] = moleFracTransportedComp * FluidSystem::molarMass(compIdx)
+ / averageMolarMassPhase;
+ }
+
+ if constexpr (isNonisothermal)
+ values[Indices::temperatureIdx] = 293.15;
+ }
}
}
return values;
}
- // \}
-
- /*!
- * \name Volume terms
- */
- // \{
-
- /*!
- * \brief Evaluates the initial value for a control volume.
+ /*!
+ * \brief Evaluates the boundary conditions for a Neumann control volume.
*
- * \param globalPos The global position
+ * \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
*/
- PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ template<class ElementVolumeVariables, class ElementFluxVariablesCache>
+ BoundaryFluxes neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
{
- PrimaryVariables values;
- values[Indices::pressureIdx] = 1.1e+5;
- values[transportCompIdx] = 0.0;
-#if NONISOTHERMAL
- values[Indices::temperatureIdx] = 283.15;
-#endif
+ BoundaryFluxes values(0.0);
- // parabolic velocity profile
- values[Indices::velocityXIdx] = inletVelocity_*(globalPos[1] - this->gridGeometry().bBoxMin()[1])*(this->gridGeometry().bBoxMax()[1] - globalPos[1])
- / (0.25*(this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1])*(this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1]));
+ if constexpr (ParentType::isMomentumProblem())
+ {
+ using FluxHelper = NavierStokesMomentumBoundaryFluxHelper;
+ values = FluxHelper::fixedPressureMomentumFlux(*this, fvGeometry, scvf,
+ elemVolVars, elemFluxVarsCache,
+ referencePressure(element, fvGeometry, scvf), true /*zeroNormalVelocityGradient*/);
+ }
+ else
+ {
+ const auto insideDensity = elemVolVars[scvf.insideScvIdx()].density();
+ values[Indices::conti0EqIdx] = this->faceVelocity(element, fvGeometry, scvf) * insideDensity * scvf.unitOuterNormal();
- values[Indices::velocityYIdx] = 0.0;
+ using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>;
+ if (isOutlet_(scvf.ipGlobal()))
+ values = FluxHelper::scalarOutflowFlux( *this, element, fvGeometry, scvf, elemVolVars);
+ }
return values;
}
- /*!
- * \brief Adds additional VTK output data to the VTKWriter.
- *
- * Function is called by the output module on every write.
+ /*!
+ * \brief Evaluates the initial value for a control volume.
*
- * \param gridVariables The grid variables
- * \param sol The solution vector
+ * \param globalPos The global position
*/
- template<class GridVariables, class SolutionVector>
- void calculateDeltaP(const GridVariables& gridVariables, const SolutionVector& sol)
+ InitialValues initialAtPos(const GlobalPosition& globalPos) const
{
- auto fvGeometry = localView(this->gridGeometry());
- auto elemVolVars = localView(gridVariables.curGridVolVars());
- for (const auto& element : elements(this->gridGeometry().gridView()))
+ InitialValues values(0.0);
+
+ if constexpr (ParentType::isMomentumProblem())
{
- fvGeometry.bindElement(element);
- elemVolVars.bindElement(element, fvGeometry, sol);
- for (auto&& scv : scvs(fvGeometry))
- {
- const auto ccDofIdx = scv.dofIndex();
- deltaP_[ccDofIdx] = elemVolVars[scv].pressure() - 1.1e5;
- }
+ // parabolic velocity profile
+ values[Indices::velocityXIdx] = parabolicProfile(globalPos[1], inletVelocity_);
+ }
+ else
+ {
+ values[Indices::pressureIdx] = 1.1e+5;
+ if constexpr (isNonisothermal)
+ values[Indices::temperatureIdx] = 283.15;
}
- }
-
- auto& getDeltaP() const
- { return deltaP_; }
- // \}
+ return values;
+ }
- void setTimeLoop(TimeLoopPtr timeLoop)
+ /*!
+ * \brief A parabolic velocity profile.
+ *
+ * \param y The position where the velocity is evaluated.
+ * \param vMax The profile's maximum velocity.
+ */
+ Scalar parabolicProfile(const Scalar y, const Scalar vMax) const
{
- timeLoop_ = timeLoop;
+ const Scalar yMin = this->gridGeometry().bBoxMin()[1];
+ const Scalar yMax = this->gridGeometry().bBoxMax()[1];
+ return vMax * (y - yMin)*(yMax - y) / (0.25*(yMax - yMin)*(yMax - yMin));
}
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ { timeLoop_ = timeLoop; }
+
Scalar time() const
+ { return timeLoop_->time(); }
+
+ //! 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<DirichletValues::dimension> hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const
{
- return timeLoop_->time();
+ std::bitset<DirichletValues::dimension> values;
+ return values;
}
-private:
+ /*!
+ * \brief Define the values of internal Dirichlet constraints for a degree of freedom.
+ * \param element The finite element
+ * \param scv The sub-control volume
+ */
+ DirichletValues internalDirichlet(const Element& element, const SubControlVolume& scv) const
+ { return DirichletValues(1.1e5); }
+private:
bool isInlet_(const GlobalPosition& globalPos) const
- {
- return globalPos[0] < eps_;
- }
+ { return globalPos[0] < eps_; }
bool isOutlet_(const GlobalPosition& globalPos) const
- {
- return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_;
- }
+ { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
const Scalar eps_;
Scalar inletVelocity_;
TimeLoopPtr timeLoop_;
- std::vector<Scalar> deltaP_;
};
} // end namespace Dumux
diff --git a/test/freeflow/navierstokesnc/channel/properties.hh b/test/freeflow/navierstokesnc/channel/properties.hh
index 68f02f994e..60e537702e 100644
--- a/test/freeflow/navierstokesnc/channel/properties.hh
+++ b/test/freeflow/navierstokesnc/channel/properties.hh
@@ -16,6 +16,10 @@
#define ENABLECACHING 1
#endif
+#ifndef USEMOLES
+#define USEMOLES 1
+#endif
+
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/cctpfa.hh>
@@ -40,10 +44,12 @@ namespace Dumux::Properties {
// Create new type tags
namespace TTag {
+struct ChannelNCTest {};
+struct ChannelNCTestMomentum { using InheritsFrom = std::tuple<ChannelNCTest, NavierStokesMomentum, FaceCenteredStaggeredModel>; };
#if !NONISOTHERMAL
-struct ChannelNCTest { using InheritsFrom = std::tuple<NavierStokesNC, StaggeredFreeFlowModel>; };
+struct ChannelNCTestMass { using InheritsFrom = std::tuple<ChannelNCTest, NavierStokesMassOnePNC, CCTpfaModel>; };
#else
-struct ChannelNCTest { using InheritsFrom = std::tuple<NavierStokesNCNI, StaggeredFreeFlowModel>; };
+struct ChannelNCTestMass { using InheritsFrom = std::tuple<ChannelNCTest, NavierStokesMassOnePNCNI, CCTpfaModel>; };
#endif
} // end namespace TTag
@@ -84,13 +90,8 @@ template<class TypeTag>
struct EnableGridFaceVariablesCache<TypeTag, TTag::ChannelNCTest> { static constexpr bool value = ENABLECACHING; };
// Use mole fraction formulation
-#if USE_MASS
template<class TypeTag>
-struct UseMoles<TypeTag, TTag::ChannelNCTest> { static constexpr bool value = false; };
-#else
-template<class TypeTag>
-struct UseMoles<TypeTag, TTag::ChannelNCTest> { static constexpr bool value = true; };
-#endif
+struct UseMoles<TypeTag, TTag::ChannelNCTest> { static constexpr bool value = USEMOLES; };
// Set the problem property
template<class TypeTag>
diff --git a/test/freeflow/navierstokesnc/maxwellstefan/main.cc b/test/freeflow/navierstokesnc/maxwellstefan/main.cc
index 21ceac603c..cbcbfc8e27 100644
--- a/test/freeflow/navierstokesnc/maxwellstefan/main.cc
+++ b/test/freeflow/navierstokesnc/maxwellstefan/main.cc
@@ -199,8 +199,7 @@ int main(int argc, char** argv)
auto massGridGeometry = std::make_shared<MassGridGeometry>(leafGridView);
// the coupling manager
- using Traits = MultiDomainTraits<MomentumTypeTag, MassTypeTag>;
- using CouplingManager = StaggeredFreeFlowCouplingManager<Traits>;
+ using CouplingManager = GetPropType<MomentumTypeTag, Properties::CouplingManager>;
auto couplingManager = std::make_shared<CouplingManager>();
// the problem (boundary conditions)
@@ -210,6 +209,7 @@ int main(int argc, char** argv)
auto massProblem = std::make_shared<MassProblem>(massGridGeometry, couplingManager);
// get some time loop parameters
+ using Traits = MultiDomainTraits<MomentumTypeTag, MassTypeTag>;
using Scalar = typename Traits::Scalar;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
diff --git a/test/freeflow/navierstokesnc/maxwellstefan/problem.hh b/test/freeflow/navierstokesnc/maxwellstefan/problem.hh
index f1bdfd7f54..8b2f320472 100644
--- a/test/freeflow/navierstokesnc/maxwellstefan/problem.hh
+++ b/test/freeflow/navierstokesnc/maxwellstefan/problem.hh
@@ -16,10 +16,10 @@
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
+#include <dumux/freeflow/navierstokes/boundarytypes.hh>
+
#include <dumux/freeflow/navierstokes/momentum/fluxhelper.hh>
#include <dumux/freeflow/navierstokes/scalarfluxhelper.hh>
-#include <dumux/freeflow/navierstokes/boundarytypes.hh>
-#include <dumux/freeflow/navierstokes/staggered/problem.hh>
namespace Dumux {
@@ -27,10 +27,10 @@ namespace Dumux {
* \ingroup NavierStokesNCTests
* \brief Test problem for the Maxwell-Stefan model
*/
-template <class TypeTag>
-class MaxwellStefanNCTestProblem : public NavierStokesStaggeredProblem<TypeTag>
+template <class TypeTag, class BaseProblem>
+class MaxwellStefanNCTestProblem : public BaseProblem
{
- using ParentType = NavierStokesStaggeredProblem<TypeTag>;
+ using ParentType = BaseProblem;
using BoundaryTypes = typename ParentType::BoundaryTypes;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
diff --git a/test/freeflow/navierstokesnc/maxwellstefan/properties.hh b/test/freeflow/navierstokesnc/maxwellstefan/properties.hh
index 35b5f5a734..55469cf76c 100644
--- a/test/freeflow/navierstokesnc/maxwellstefan/properties.hh
+++ b/test/freeflow/navierstokesnc/maxwellstefan/properties.hh
@@ -14,8 +14,11 @@
#include <dune/grid/yaspgrid.hh>
-#include <dumux/freeflow/navierstokes/momentum/model.hh>
#include <dumux/freeflow/navierstokes/mass/1pnc/model.hh>
+#include <dumux/freeflow/navierstokes/mass/problem.hh>
+
+#include <dumux/freeflow/navierstokes/momentum/problem.hh>
+#include <dumux/freeflow/navierstokes/momentum/model.hh>
#include <dumux/flux/maxwellstefanslaw.hh>
#include <dumux/material/fluidsystems/base.hh>
@@ -23,7 +26,8 @@
#include <dumux/discretization/fcstaggered.hh>
#include <dumux/discretization/cctpfa.hh>
-#include <dumux/multidomain/staggeredfreeflow/couplingmanager.hh>
+#include <dumux/multidomain/traits.hh>
+#include <dumux/multidomain/freeflow/couplingmanager.hh>
#include "problem.hh"
@@ -41,11 +45,17 @@ struct ReplaceCompEqIdx<TypeTag, TTag::MaxwellStefanNCTest> { static constexpr i
// Set the grid type
template<class TypeTag>
-struct Grid<TypeTag, TTag::MaxwellStefanNCTest> { using type = Dune::YaspGrid<2>; };
+struct Grid<TypeTag, TTag::MaxwellStefanNCTest>
+{ using type = Dune::YaspGrid<2>; };
// Set the problem property
template<class TypeTag>
-struct Problem<TypeTag, TTag::MaxwellStefanNCTest> { using type = Dumux::MaxwellStefanNCTestProblem<TypeTag> ; };
+struct Problem<TypeTag, TTag::MaxwellStefanTestMomentum>
+{ using type = MaxwellStefanNCTestProblem<TypeTag, Dumux::NavierStokesMomentumProblem<TypeTag>>; };
+
+template<class TypeTag>
+struct Problem<TypeTag, TTag::MaxwellStefanTestMass>
+{ using type = MaxwellStefanNCTestProblem<TypeTag, Dumux::NavierStokesMassProblem<TypeTag>>; };
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::MaxwellStefanNCTest> { static constexpr bool value = true; };
@@ -61,13 +71,6 @@ struct UseMoles<TypeTag, TTag::MaxwellStefanNCTest> { static constexpr bool valu
template<class TypeTag>
struct MolecularDiffusionType<TypeTag, TTag::MaxwellStefanNCTest> { using type = MaxwellStefansLaw<TypeTag>; };
-template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::MaxwellStefanNCTest>
-{
- using Traits = MultiDomainTraits<TTag::MaxwellStefanTestMomentum, TTag::MaxwellStefanTestMass>;
- using type = StaggeredFreeFlowCouplingManager<Traits>;
-};
-
/*!
* \ingroup NavierStokesNCTests
* \brief A simple fluid system with three components for testing the multi-component diffusion with the Maxwell-Stefan formulation.
@@ -194,6 +197,14 @@ public:
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::MaxwellStefanNCTest> { using type = MaxwellStefanFluidSystem<TypeTag>; };
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::MaxwellStefanNCTest>
+{
+ using Traits = MultiDomainTraits<TTag::MaxwellStefanTestMomentum, TTag::MaxwellStefanTestMass>;
+ using type = FreeFlowCouplingManager<Traits>;
+};
+
+
} // end namespace Dumux::Properties
#endif
--
GitLab