diff --git a/test/freeflow/navierstokes/angeli/CMakeLists.txt b/test/freeflow/navierstokes/angeli/CMakeLists.txt
index cb411f6e9a0d0f26db0f7422b50e5a280dddcba9..ff3f34615c093302af3a102bccb5fcc3b91c02ea 100644
--- a/test/freeflow/navierstokes/angeli/CMakeLists.txt
+++ b/test/freeflow/navierstokes/angeli/CMakeLists.txt
@@ -20,7 +20,6 @@ dumux_add_test(NAME test_ff_navierstokes_angeli_averaged
--command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_angeli params.input
-TimeLoop.TEnd 1e-6 -TimeLoop.DtInitial 1e-6
-Problem.Name test_ff_navierstokes_angeli_averaged
- -Problem.UseVelocityAveragingForInitial true
- -Problem.UseVelocityAveragingForDirichlet true")
+ -Problem.InterpolateExactVelocity true")
dune_symlink_to_source_files(FILES "params.input")
diff --git a/test/freeflow/navierstokes/angeli/main.cc b/test/freeflow/navierstokes/angeli/main.cc
index 7da33e3f6e4921118e1009fd798146ed734e0a43..1911467ee65194bfaafc25d33a8063bfec24b1f9 100644
--- a/test/freeflow/navierstokes/angeli/main.cc
+++ b/test/freeflow/navierstokes/angeli/main.cc
@@ -32,18 +32,20 @@
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
-#include <dune/grid/io/file/vtk.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/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager_yasp.hh>
+
#include <dumux/linear/seqsolverbackend.hh>
-#include <dumux/nonlinear/newtonsolver.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/traits.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+
+#include <dumux/freeflow/navierstokes/velocityoutput.hh>
#include <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
#include <test/freeflow/navierstokes/errors.hh>
@@ -53,8 +55,10 @@ int main(int argc, char** argv)
{
using namespace Dumux;
- // define the type tag for this problem
- using TypeTag = Properties::TTag::AngeliTest;
+ // define the type tag for this problem
+ using MomentumTypeTag = Properties::TTag::AngeliTestMomentum;
+ using MassTypeTag = Properties::TTag::AngeliTestMass;
+
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
@@ -67,7 +71,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();
////////////////////////////////////////////////////////////
@@ -78,11 +82,27 @@ 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 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);
// get some time loop parameters
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using Scalar = typename Traits::Scalar;
const auto tStart = getParam<Scalar>("TimeLoop.TStart", 0.0);
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
@@ -92,62 +112,72 @@ int main(int argc, char** argv)
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(tStart, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
- // the problem (initial and boundary conditions)
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- auto problem = std::make_shared<Problem>(gridGeometry);
- problem->setTime(timeLoop->time());
-
// 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());
- problem->applyInitialSolution(x);
+ momentumProblem->applyInitialSolution(x[momentumIdx]);
+ massProblem->applyInitialSolution(x[massIdx]);
auto xOld = x;
// 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, xOld);
- // initialize the vtk output module
- StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
- using IOFields = GetPropType<TypeTag, Properties::IOFields>;
+ massGridVariables->init(x[massIdx]);
+ momentumGridVariables->init(x[momentumIdx]);
+
+ // intialize 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
- NavierStokesAnalyticalSolutionVectors analyticalSolVectors(problem, tStart);
- vtkWriter.addField(analyticalSolVectors.getAnalyticalPressureSolution(), "pressureExact");
- vtkWriter.addField(analyticalSolVectors.getAnalyticalVelocitySolution(), "velocityExact");
- vtkWriter.addFaceField(analyticalSolVectors.getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
+
+ vtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<MassGridVariables>>());
+ NavierStokesTest::AnalyticalSolutionVectors analyticalSolVectors(momentumProblem, massProblem);
+ vtkWriter.addField(analyticalSolVectors.analyticalPressureSolution(), "pressureExact");
+ vtkWriter.addField(analyticalSolVectors.analyticalVelocitySolution(), "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, timeLoop, xOld);
+ 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::UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
- using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
- NewtonSolver nonLinearSolver(assembler, linearSolver);
+ // the non-linear solver
+ using NewtonSolver = Dumux::MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
// the discrete L2 and Linfity errors
const bool printErrors = getParam<bool>("Problem.PrintErrors", false);
- NavierStokesErrors errors(problem, x, 0.0);
- NavierStokesErrorCSVWriter errorCSVWriter(problem);
+ NavierStokesTest::Errors errors(momentumProblem, massProblem, x);
+ NavierStokesTest::ErrorCSVWriter errorCSVWriter(momentumProblem, massProblem);
// time loop
timeLoop->start(); do
{
- problem->setTime(timeLoop->time() + timeLoop->timeStepSize());
+ massProblem->updateTime(timeLoop->time() + timeLoop->timeStepSize());
+ momentumProblem->updateTime(timeLoop->time() + timeLoop->timeStepSize());
// solve the non-linear system with time step control
nonLinearSolver.solve(x, *timeLoop);
// make the new solution the old solution
xOld = x;
- gridVariables->advanceTimeStep();
+ momentumGridVariables->advanceTimeStep();
+ massGridVariables->advanceTimeStep();
// print discrete L2 and Linfity errors
if (printErrors)
@@ -161,6 +191,7 @@ int main(int argc, char** argv)
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
+ analyticalSolVectors.update(timeLoop->time());
// write vtk output
vtkWriter.write(timeLoop->time());
@@ -175,6 +206,7 @@ int main(int argc, char** argv)
timeLoop->finalize(leafGridView.comm());
+
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
diff --git a/test/freeflow/navierstokes/angeli/params.input b/test/freeflow/navierstokes/angeli/params.input
index f6d7188543b2c1bdc414bfe99dfd741f58570503..d5b8163b7957832dfb816942d7f98ed23a3839e1 100644
--- a/test/freeflow/navierstokes/angeli/params.input
+++ b/test/freeflow/navierstokes/angeli/params.input
@@ -12,8 +12,6 @@ Name = test_angeli # name passed to the output routines
EnableGravity = false
PrintErrors = false
EnableInertiaTerms = true
-UseVelocityAveragingForInitial = false
-UseVelocityAveragingForDirichlet = false
[Component]
LiquidDensity = 1
diff --git a/test/freeflow/navierstokes/angeli/problem.hh b/test/freeflow/navierstokes/angeli/problem.hh
index fd37c83dec48e7badb868bfb260e275d94e683b9..546a89f2b04e13a3f7c338d0e9ed25100f7a8721 100644
--- a/test/freeflow/navierstokes/angeli/problem.hh
+++ b/test/freeflow/navierstokes/angeli/problem.hh
@@ -51,30 +51,36 @@ class AngeliTestProblem : 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 GridView = typename GridGeometry::GridView;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
- using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
- using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
+ using NumEqVector = typename ParentType::NumEqVector;
+ using PrimaryVariables = typename ParentType::PrimaryVariables;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
-
- using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
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 FVElementGeometry::Element;
+
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+
public:
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
- AngeliTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
- : ParentType(gridGeometry)
+ AngeliTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(gridGeometry, couplingManager)
{
kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
- rho_ = getParam<Scalar>("Component.LiquidDensity", 1.0);
- useVelocityAveragingForDirichlet_ = getParam<bool>("Problem.UseVelocityAveragingForDirichlet", false);
- useVelocityAveragingForInitial_ = getParam<bool>("Problem.UseVelocityAveragingForInitial", false);
+ useNeumann_ = getParam<bool>("Problem.UseNeumann", false);
+ interpolateExactVelocity_ = getParam<bool>("Problem.InterpolateExactVelocity", false);
}
/*!
@@ -99,34 +105,63 @@ public:
{
BoundaryTypes values;
- // set Dirichlet values for the velocity everywhere
- values.setDirichlet(Indices::velocityXIdx);
- values.setDirichlet(Indices::velocityYIdx);
+ if constexpr (ParentType::isMomentumProblem())
+ {
+ static constexpr Scalar eps = 1e-8;
+ if (useNeumann_ && (globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps ||
+ globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps))
+ values.setAllNeumann();
+ else
+ values.setAllDirichlet();
+ }
+ else
+ values.setAllNeumann();
return values;
+
}
+ //! Enable internal Dirichlet constraints
+ static constexpr bool enableInternalDirichletConstraints()
+ { return !ParentType::isMomentumProblem(); }
+
+
/*!
- * \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
+ * \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 fvGeometry The finite-volume geometry
- * \param scv The sub control volume
- * \param pvIdx The primary variable index in the solution vector
+ * \param scv The sub-control volume
*/
- bool isDirichletCell(const Element& element,
- const typename GridGeometry::LocalView& fvGeometry,
- const typename GridGeometry::SubControlVolume& scv,
- int pvIdx) const
+ std::bitset<PrimaryVariables::dimension> hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const
{
- // set a fixed pressure in all cells at the boundary
- for (const auto& scvf : scvfs(fvGeometry))
- if (scvf.boundary())
- return true;
+ std::bitset<PrimaryVariables::dimension> values;
+
+ // We don't need internal Dirichlet conditions if a Neumann BC is set for the momentum balance (which accounts for the pressure).
+ // If only Dirichlet BCs are set for the momentum balance, fix the pressure at some cells such that the solution is fully defined.
+ if (!useNeumann_)
+ {
+ const auto fvGeometry = localView(this->gridGeometry()).bindElement(element);
+ if (fvGeometry.hasBoundaryScvf())
+ values.set(Indices::pressureIdx);
+ }
- return false;
+ 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 analyticalSolution(scv.center(), time_); }
+
+
+
/*!
* \brief Evaluate the boundary conditions for a dirichlet
* control volume face (velocities)
@@ -138,26 +173,58 @@ public:
*/
PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
{
- if (useVelocityAveragingForDirichlet_)
- return averagedVelocity_(scvf, time_);
+ if (ParentType::isMomentumProblem() && interpolateExactVelocity_)
+ return velocityDirichlet_(scvf);
else
return analyticalSolution(scvf.center(), time_);
}
/*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * control volume face (pressure)
+ * \brief Evaluates the boundary conditions for a Neumann control volume.
*
- * \param element The finite element
- * \param scv the sub 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
*/
- PrimaryVariables dirichlet(const Element& element, const SubControlVolume& scv) const
+ template<class ElementVolumeVariables, class ElementFluxVariablesCache>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
{
- PrimaryVariables priVars(0.0);
- priVars[Indices::pressureIdx] = analyticalSolution(scv.center(), time_)[Indices::pressureIdx];
- return priVars;
+ NumEqVector values(0.0);
+
+ if constexpr (!ParentType::isMomentumProblem())
+ {
+ const auto insideDensity = elemVolVars[scvf.insideScvIdx()].density();
+ values[Indices::conti0EqIdx] = this->faceVelocity(element, fvGeometry, scvf) * insideDensity * scvf.unitOuterNormal();
+ }
+ else
+ {
+ using std::sin;
+ using std::cos;
+ Dune::FieldMatrix<Scalar, 2, 2> momentumFlux(0.0);
+ const auto x = scvf.ipGlobal()[0];
+ const auto y = scvf.ipGlobal()[1];
+ const Scalar mu = kinematicViscosity_;
+ const Scalar t = time_;
+ momentumFlux[0][0] = M_PI*M_PI *(-4.0*mu*exp(10.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y) + (4.0*sin(2*M_PI*y)*sin(2*M_PI*y)*cos(M_PI*x)*cos(M_PI*x) - 1.0*cos(2*M_PI*x) - 0.25*cos(4*M_PI*y))*exp(5.0*M_PI*M_PI*mu*t))*exp(-15.0*M_PI*M_PI*mu*t);
+ momentumFlux[0][1] = M_PI*M_PI *( 3.0*mu*exp(10.0*M_PI*M_PI*mu*t) - 2.0*exp(5.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y))*exp(-15.0*M_PI*M_PI*mu*t)*cos(M_PI*x)*cos(2*M_PI*y);
+ momentumFlux[1][0] = momentumFlux[0][1];
+ momentumFlux[1][1] = M_PI*M_PI *( 4.0*mu*exp(10.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y) + (sin(M_PI*x)*sin(M_PI*x)*cos(2*M_PI*y)*cos(2*M_PI*y) - 1.0*cos(2*M_PI*x) - 0.25*cos(4*M_PI*y))*exp(5.0*M_PI*M_PI*mu*t))*exp(-15.0*M_PI*M_PI*mu*t);
+
+ const auto normal = scvf.unitOuterNormal();
+ momentumFlux.mv(normal, values);
+ }
+
+ return values;
}
+
+
/*!
* \brief Returns the analytical solution of the problem at a given time and position.
* \param globalPos The global position
@@ -169,10 +236,17 @@ public:
const Scalar y = globalPos[1];
PrimaryVariables values;
+ using std::exp;
+ using std::sin;
+ using std::cos;
- values[Indices::pressureIdx] = - 0.25 * std::exp(-10.0 * kinematicViscosity_ * M_PI * M_PI * time) * M_PI * M_PI * (4.0 * std::cos(2.0 * M_PI * x) + std::cos(4.0 * M_PI * y))*rho_;
- values[Indices::velocityXIdx] = - 2.0 * M_PI * std::exp(- 5.0 * kinematicViscosity_ * M_PI * M_PI * time) * std::cos(M_PI * x) * std::sin(2.0 * M_PI * y);
- values[Indices::velocityYIdx] = M_PI * std::exp(- 5.0 * kinematicViscosity_ * M_PI * M_PI * time) * std::sin(M_PI * x) * std::cos(2.0 * M_PI * y);
+ if constexpr (ParentType::isMomentumProblem())
+ {
+ values[Indices::velocityXIdx] = - 2.0 * M_PI * exp(- 5.0 * kinematicViscosity_ * M_PI * M_PI * time) * cos(M_PI * x) * sin(2.0 * M_PI * y);
+ values[Indices::velocityYIdx] = M_PI * exp(- 5.0 * kinematicViscosity_ * M_PI * M_PI * time) * sin(M_PI * x) * cos(2.0 * M_PI * y);
+ }
+ else
+ values[Indices::pressureIdx] = - 0.25 * exp(-10.0 * kinematicViscosity_ * M_PI * M_PI * time) * M_PI * M_PI * (4.0 * cos(2.0 * M_PI * x) + cos(4.0 * M_PI * y));
return values;
}
@@ -185,31 +259,42 @@ public:
// \{
/*!
- * \brief Evaluates the initial value for a control volume (pressure)
+ * \brief Evaluates the initial value for a control volume (velocity)
*/
PrimaryVariables initial(const SubControlVolume& scv) const
{
- PrimaryVariables priVars(0.0);
- priVars[Indices::pressureIdx] = analyticalSolution(scv.center(), time_)[Indices::pressureIdx];
- return priVars;
+ // We can either evaluate the analytical solution point-wise at the
+ // momentum balance integration points or use an averaged velocity.
+
+ // Not using the quadrature to integrate and average the velocity
+ // will yield a spurious pressure solution in the first time step
+ // because the discrete mass balance equation is not fulfilled when just taking
+ // point-wise values of the analytical solution.
+
+ if (!interpolateExactVelocity_)
+ return analyticalSolution(scv.dofPosition(), 0.0);
+
+ // Get the element intersection/facet corresponding corresponding to the dual grid scv
+ // (where the velocity DOFs are located) and use a quadrature to get the average velocity
+ // on that facet.
+ const auto& element = this->gridGeometry().element(scv.elementIndex());
+
+ for (const auto& intersection : intersections(this->gridGeometry().gridView(), element))
+ {
+ if (intersection.indexInInside() == scv.indexInElement())
+ return velocityDirichlet_(intersection);
+
+ }
+ DUNE_THROW(Dune::InvalidStateException, "No intersection found");
}
+
/*!
- * \brief Evaluates the initial value for a sub control volume face (velocities)
- *
- * Simply assigning the value of the analytical solution at the face center
- * gives a discrete solution that is not divergence-free. For small initial
- * time steps, this has a negative impact on the pressure solution
- * after the first time step. The flag UseVelocityAveragingForInitial triggers the
- * function averagedVelocity_ which uses a higher order quadrature formula to
- * bring the discrete solution sufficiently close to being divergence-free.
+ * \brief Evaluates the initial value for an element (pressure)
*/
- PrimaryVariables initial(const SubControlVolumeFace& scvf) const
+ PrimaryVariables initial(const Element& element) const
{
- if (useVelocityAveragingForInitial_)
- return averagedVelocity_(scvf, time_);
- else
- return analyticalSolution(scvf.center(), time_);
+ return analyticalSolution(element.geometry().center(), 0.0);
}
// \}
@@ -217,34 +302,33 @@ public:
/*!
* \brief Updates the time information
*/
- void setTime(Scalar t)
+ void updateTime(Scalar t)
{
time_ = t;
}
private:
- PrimaryVariables averagedVelocity_(const SubControlVolumeFace& scvf, Scalar t) const
+ template<class Entity>
+ PrimaryVariables velocityDirichlet_(const Entity& entity) const
{
PrimaryVariables priVars(0.0);
- const auto geo = scvf.geometry();
- const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension-1>::rule(geo.type(), 3);
+ const auto geo = entity.geometry();
+ const auto& quad = Dune::QuadratureRules<Scalar, decltype(geo)::mydimension>::rule(geo.type(), 3);
for (auto&& qp : quad)
{
const auto w = qp.weight()*geo.integrationElement(qp.position());
const auto globalPos = geo.global(qp.position());
- const auto sol = analyticalSolution(globalPos, t);
- priVars[Indices::velocityXIdx] += sol[Indices::velocityXIdx]*w;
- priVars[Indices::velocityYIdx] += sol[Indices::velocityYIdx]*w;
+ const auto sol = analyticalSolution(globalPos, time_);
+ priVars += sol*w;
}
- priVars[Indices::velocityXIdx] /= scvf.area();
- priVars[Indices::velocityYIdx] /= scvf.area();
+ priVars /= geo.volume();
return priVars;
}
- Scalar kinematicViscosity_, rho_;
- Scalar time_ = 0;
- bool useVelocityAveragingForDirichlet_;
- bool useVelocityAveragingForInitial_;
+ Scalar kinematicViscosity_;
+ Scalar time_ = 0.0;
+ bool useNeumann_;
+ bool interpolateExactVelocity_;
};
} // end namespace Dumux
diff --git a/test/freeflow/navierstokes/angeli/properties.hh b/test/freeflow/navierstokes/angeli/properties.hh
index 51217da6891b4103e0158cfea590a62d6e9e74db..976b49397d9776922d8f2dc4ad0b9fb5d9f28b5f 100644
--- a/test/freeflow/navierstokes/angeli/properties.hh
+++ b/test/freeflow/navierstokes/angeli/properties.hh
@@ -26,19 +26,27 @@
#define DUMUX_ANGELI_TEST_PROPERTIES_HH
#include <dune/grid/yaspgrid.hh>
-#include <dumux/discretization/staggered/freeflow/properties.hh>
-#include <dumux/freeflow/navierstokes/model.hh>
+#include <dumux/freeflow/navierstokes/momentum/model.hh>
+#include <dumux/freeflow/navierstokes/mass/1p/model.hh>
+
+#include <dumux/discretization/fcstaggered.hh>
+#include <dumux/discretization/cctpfa.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 AngeliTest { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
+struct AngeliTest {};
+struct AngeliTestMomentum { using InheritsFrom = std::tuple<AngeliTest, NavierStokesMomentum, FaceCenteredStaggeredModel>; };
+struct AngeliTestMass { using InheritsFrom = std::tuple<AngeliTest, NavierStokesMassOneP, CCTpfaModel>; };
} // end namespace TTag
// the fluid system
@@ -66,6 +74,15 @@ struct EnableGridFluxVariablesCache<TypeTag, TTag::AngeliTest> { static constexp
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::AngeliTest> { static constexpr bool value = true; };
-} // end namespace Dumux::Properties
+// Set the problem property
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::AngeliTest>
+{
+private:
+ using Traits = MultiDomainTraits<TTag::AngeliTestMomentum, TTag::AngeliTestMass>;
+public:
+ using type = StaggeredFreeFlowCouplingManager<Traits>;
+};
+} // end namespace Dumux::Properties
#endif