From 96855f6c07bb35b173edf72cc924e258b420b41e Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Thu, 17 Dec 2020 16:09:13 +0100
Subject: [PATCH] [test][navierstokesnc] Port 1pnc MaxwellStefan test
---
.../navierstokes/mass/1pnc/CMakeLists.txt | 3 +
.../navierstokes/mass/1pnc/localresidual.hh | 3 +-
.../freeflow/navierstokes/mass/1pnc/model.hh | 11 +-
test/freeflow/navierstokes/angeli/main.cc | 1 -
test/freeflow/navierstokesnc/channel/main.cc | 32 ++-
.../navierstokesnc/channel/problem.hh | 6 +-
.../navierstokesnc/channel/properties.hh | 40 ++-
.../maxwellstefan/CMakeLists.txt | 10 +-
.../navierstokesnc/maxwellstefan/main.cc | 215 +++++++++++---
.../navierstokesnc/maxwellstefan/params.input | 3 +-
.../navierstokesnc/maxwellstefan/problem.hh | 265 +++++++-----------
.../maxwellstefan/properties.hh | 28 +-
12 files changed, 378 insertions(+), 239 deletions(-)
create mode 100644 dumux/freeflow/navierstokes/mass/1pnc/CMakeLists.txt
diff --git a/dumux/freeflow/navierstokes/mass/1pnc/CMakeLists.txt b/dumux/freeflow/navierstokes/mass/1pnc/CMakeLists.txt
new file mode 100644
index 0000000000..80a3545121
--- /dev/null
+++ b/dumux/freeflow/navierstokes/mass/1pnc/CMakeLists.txt
@@ -0,0 +1,3 @@
+file(GLOB DUMUX_FREEFLOW_NAVIERSTOKES_MASS_1PNC_HEADERS *.hh *.inc)
+install(FILES ${DUMUX_FREEFLOW_NAVIERSTOKES_MASS_1PNC_HEADERS}
+ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/navierstokes/mass/1pnc)
diff --git a/dumux/freeflow/navierstokes/mass/1pnc/localresidual.hh b/dumux/freeflow/navierstokes/mass/1pnc/localresidual.hh
index 5878df9c16..eea71d29d8 100644
--- a/dumux/freeflow/navierstokes/mass/1pnc/localresidual.hh
+++ b/dumux/freeflow/navierstokes/mass/1pnc/localresidual.hh
@@ -26,6 +26,7 @@
#include <dumux/common/numeqvector.hh>
#include <dumux/common/properties.hh>
+#include <dumux/discretization/method.hh>
namespace Dumux {
@@ -58,7 +59,7 @@ class NavierStokesMassOnePNCLocalResidual : public GetPropType<TypeTag, Properti
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
- static_assert(GridGeometry::discMethod == DiscretizationMethod::cctpfa);
+ static_assert(GridGeometry::discMethod == DiscretizationMethods::cctpfa, "");
static constexpr bool useMoles = ModelTraits::useMoles();
static constexpr auto numComponents = ModelTraits::numFluidComponents();
diff --git a/dumux/freeflow/navierstokes/mass/1pnc/model.hh b/dumux/freeflow/navierstokes/mass/1pnc/model.hh
index 73fb1d51a9..34eb915401 100644
--- a/dumux/freeflow/navierstokes/mass/1pnc/model.hh
+++ b/dumux/freeflow/navierstokes/mass/1pnc/model.hh
@@ -53,17 +53,8 @@
#include "volumevariables.hh"
#include "fluxvariables.hh"
#include "indices.hh"
-#include "./../../iofields.hh"
#include "iofields.hh"
-#include <dumux/flux/cctpfa/fourierslaw.hh>
-
-#include <dumux/discretization/method.hh>
-#include <dumux/freeflow/navierstokes/energy/model.hh>
-#include <dumux/flux/fickslaw.hh>
-#include <dumux/freeflow/navierstokes/scalarfluxvariablescachefiller.hh>
-
-
namespace Dumux {
/*!
@@ -360,7 +351,7 @@ struct ThermalConductivityModel<TypeTag, TTag::NavierStokesMassOnePNCNI>
template<class TypeTag>
struct HeatConductionType<TypeTag, TTag::NavierStokesMassOnePNCNI>
-{ using type = FouriersLawImplementation<TypeTag, DiscretizationMethod::cctpfa>; };
+{ using type = FouriersLaw<TypeTag>; };
//! The flux variables
template<class TypeTag>
diff --git a/test/freeflow/navierstokes/angeli/main.cc b/test/freeflow/navierstokes/angeli/main.cc
index 02013dc541..81208a3d97 100644
--- a/test/freeflow/navierstokes/angeli/main.cc
+++ b/test/freeflow/navierstokes/angeli/main.cc
@@ -147,7 +147,6 @@ int main(int argc, char** argv)
using LinearSolver = Dumux::UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
auto linearSolver = std::make_shared<LinearSolver>();
- // the non-linear solver
// the non-linear solver
using NewtonSolver = Dumux::MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
diff --git a/test/freeflow/navierstokesnc/channel/main.cc b/test/freeflow/navierstokesnc/channel/main.cc
index f51db75508..066e4e4a15 100644
--- a/test/freeflow/navierstokesnc/channel/main.cc
+++ b/test/freeflow/navierstokesnc/channel/main.cc
@@ -17,21 +17,22 @@
#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/initialize.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
-#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/linear/istlsolvers.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <dumux/linear/linearalgebratraits.hh>
-#include <dumux/nonlinear/newtonsolver.hh>
+
+#include <dumux/multidomain/newtonsolver.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/traits.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/freeflow/navierstokes/velocityoutput.hh>
#include "properties.hh"
@@ -54,7 +55,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<MassTypeTag, Properties::Grid>> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
@@ -68,9 +69,15 @@ int main(int argc, char** argv)
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
- // the problem (initial and boundary conditions)
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- auto problem = std::make_shared<Problem>(gridGeometry);
+ // the coupling manager
+ using CouplingManager = GetPropType<MomentumTypeTag, Properties::CouplingManager>;
+ 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>;
@@ -84,7 +91,10 @@ int main(int argc, char** argv)
problem->setTimeLoop(timeLoop);
// the solution vector
- using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+ constexpr auto momentumIdx = CouplingManager::freeFlowMomentumIndex;
+ constexpr auto massIdx = CouplingManager::freeFlowMassIndex;
+ using Traits = MultiDomainTraits<MomentumTypeTag, MassTypeTag>;
+ using SolutionVector = typename Traits::SolutionVector;
SolutionVector x;
x[GridGeometry::cellCenterIdx()].resize(gridGeometry->numCellCenterDofs());
x[GridGeometry::faceIdx()].resize(gridGeometry->numFaceDofs());
diff --git a/test/freeflow/navierstokesnc/channel/problem.hh b/test/freeflow/navierstokesnc/channel/problem.hh
index 851a90ea4e..06c4699b3f 100644
--- a/test/freeflow/navierstokesnc/channel/problem.hh
+++ b/test/freeflow/navierstokesnc/channel/problem.hh
@@ -15,11 +15,11 @@
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
-#include <dumux/common/timeloop.hh>
-#include <dumux/common/numeqvector.hh>
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
-#include <dumux/freeflow/navierstokes/staggered/problem.hh>
+
+#include <dumux/freeflow/navierstokes/momentum/fluxhelper.hh>
+#include <dumux/freeflow/navierstokes/scalarfluxhelper.hh>
namespace Dumux {
diff --git a/test/freeflow/navierstokesnc/channel/properties.hh b/test/freeflow/navierstokesnc/channel/properties.hh
index a2657253a5..68f02f994e 100644
--- a/test/freeflow/navierstokesnc/channel/properties.hh
+++ b/test/freeflow/navierstokesnc/channel/properties.hh
@@ -18,13 +18,22 @@
#include <dune/grid/yaspgrid.hh>
-#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/discretization/cctpfa.hh>
+#include <dumux/discretization/fcstaggered.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/freeflow/compositional/navierstokesncmodel.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/1padapter.hh>
#include <dumux/material/fluidsystems/h2oair.hh>
+#include <dumux/multidomain/traits.hh>
+#include <dumux/multidomain/freeflow/couplingmanager.hh>
+
#include "problem.hh"
namespace Dumux::Properties {
@@ -38,6 +47,15 @@ struct ChannelNCTest { using InheritsFrom = std::tuple<NavierStokesNCNI, Stagger
#endif
} // end namespace TTag
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::ChannelNCTestMomentum>
+{ using type = ChannelNCTestProblem<TypeTag, Dumux::NavierStokesMomentumProblem<TypeTag>>; };
+
+template<class TypeTag>
+struct Problem<TypeTag, TTag::ChannelNCTestMass>
+{ using type = ChannelNCTestProblem<TypeTag, Dumux::NavierStokesMassProblem<TypeTag>>; };
+
// Select the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::ChannelNCTest>
@@ -48,15 +66,13 @@ struct FluidSystem<TypeTag, TTag::ChannelNCTest>
};
template<class TypeTag>
-struct ReplaceCompEqIdx<TypeTag, TTag::ChannelNCTest> { static constexpr int value = 0; };
+struct ReplaceCompEqIdx<TypeTag, TTag::ChannelNCTest>
+{ static constexpr int value = 0; };
// Set the grid type
template<class TypeTag>
-struct Grid<TypeTag, TTag::ChannelNCTest> { using type = Dune::YaspGrid<2>; };
-
-// Set the problem property
-template<class TypeTag>
-struct Problem<TypeTag, TTag::ChannelNCTest> { using type = Dumux::ChannelNCTestProblem<TypeTag> ; };
+struct Grid<TypeTag, TTag::ChannelNCTest>
+{ using type = Dune::YaspGrid<2>; };
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::ChannelNCTest> { static constexpr bool value = ENABLECACHING; };
@@ -76,6 +92,14 @@ template<class TypeTag>
struct UseMoles<TypeTag, TTag::ChannelNCTest> { static constexpr bool value = true; };
#endif
+// Set the problem property
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::ChannelNCTest>
+{
+ using Traits = MultiDomainTraits<TTag::ChannelNCTestMomentum, TTag::ChannelNCTestMass>;
+ using type = FreeFlowCouplingManager<Traits>;
+};
+
} // end namespace Dumux::Properties
#endif
diff --git a/test/freeflow/navierstokesnc/maxwellstefan/CMakeLists.txt b/test/freeflow/navierstokesnc/maxwellstefan/CMakeLists.txt
index da04883876..23a5ef924b 100644
--- a/test/freeflow/navierstokesnc/maxwellstefan/CMakeLists.txt
+++ b/test/freeflow/navierstokesnc/maxwellstefan/CMakeLists.txt
@@ -1,8 +1,11 @@
+
# SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
# SPDX-License-Identifier: GPL-3.0-or-later
+dune_symlink_to_source_files(FILES "params.input")
+
dumux_add_test(NAME test_ff_stokes2c_maxwellstefan
- LABELS freeflow navierstokes
+ LABELS freeflow navierstokesnc navierstokes
SOURCES main.cc
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
@@ -10,6 +13,5 @@ dumux_add_test(NAME test_ff_stokes2c_maxwellstefan
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_stokes2c_maxwellstefan-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_ff_stokes2c_maxwellstefan-00005.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_stokes2c_maxwellstefan params.input
- -Problem.Name test_ff_stokes2c_maxwellstefan")
-
-dune_symlink_to_source_files(FILES "params.input")
+ -Problem.Name test_ff_stokes2c_maxwellstefan
+ -Problem.PlotOutput false")
diff --git a/test/freeflow/navierstokesnc/maxwellstefan/main.cc b/test/freeflow/navierstokesnc/maxwellstefan/main.cc
index a0fc14d160..21ceac603c 100644
--- a/test/freeflow/navierstokesnc/maxwellstefan/main.cc
+++ b/test/freeflow/navierstokesnc/maxwellstefan/main.cc
@@ -14,33 +14,161 @@
#include <ctime>
#include <iostream>
+#include <type_traits>
+#include <tuple>
#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/initialize.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
+
#include <dumux/io/grid/gridmanager_yasp.hh>
-#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/gnuplotinterface.hh>
+#include <dumux/freeflow/navierstokes/velocityoutput.hh>
+
#include <dumux/linear/istlsolvers.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <dumux/linear/linearalgebratraits.hh>
#include <dumux/nonlinear/newtonsolver.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/traits.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+
#include "properties.hh"
+namespace Dumux {
+
+template<class Scalar>
+class PlotConcentration
+{
+public:
+ /*!
+ * \brief Writes out the diffusion rates from left to right
+ *
+ * Called after every time step
+ *
+ * \param curSol Vector containing the current solution
+ * \param gridVariables The grid variables
+ * \param time The time
+ */
+ template<class SolutionVector, class GridVariables>
+ void plotComponentsOverTime(const SolutionVector& curSol,
+ const GridVariables& gridVariables,
+ const Scalar time)
+ {
+ using FluidSystem = typename GridVariables::VolumeVariables::FluidSystem;
+ const auto& gridGeometry = gridVariables.curGridVolVars().problem().gridGeometry();
+ Scalar x_co2_left = 0.0;
+ Scalar x_n2_left = 0.0;
+ Scalar x_co2_right = 0.0;
+ Scalar x_n2_right = 0.0;
+ Scalar x_h2_left = 0.0;
+ Scalar x_h2_right = 0.0;
+ Scalar i = 0.0;
+ Scalar j = 0.0;
+ for (const auto& element : elements(gridGeometry.gridView()))
+ {
+ auto fvGeometry = localView(gridGeometry);
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(gridVariables.curGridVolVars());
+ elemVolVars.bind(element, fvGeometry, curSol);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto globalPos = scv.dofPosition();
+
+ if (globalPos[0] < 0.5)
+ {
+ x_co2_left += elemVolVars[scv].moleFraction(FluidSystem::CO2Idx);
+ x_n2_left += elemVolVars[scv].moleFraction(FluidSystem::N2Idx);
+ x_h2_left += elemVolVars[scv].moleFraction(FluidSystem::H2Idx);
+ i +=1;
+ }
+ else
+ {
+ x_co2_right += elemVolVars[scv].moleFraction(FluidSystem::CO2Idx);
+ x_n2_right += elemVolVars[scv].moleFraction(FluidSystem::N2Idx);
+ x_h2_right += elemVolVars[scv].moleFraction(FluidSystem::H2Idx);
+ j +=1;
+ }
+ }
+ }
+ x_co2_left /= i;
+ x_n2_left /= i;
+ x_h2_left /= i;
+ x_co2_right /= j;
+ x_n2_right /= j;
+ x_h2_right /= j;
+
+ //do a gnuplot
+ x_.push_back(time); // in seconds
+ y1_.push_back(x_n2_left);
+ y2_.push_back(x_n2_right);
+ y3_.push_back(x_co2_left);
+ y4_.push_back(x_co2_right);
+ y5_.push_back(x_h2_left);
+ y6_.push_back(x_h2_right);
+
+ gnuplot_.resetPlot();
+ gnuplot_.setXRange(0, std::max(time, 72000.0));
+ gnuplot_.setYRange(0.4, 0.6);
+ gnuplot_.setXlabel("time [s]");
+ gnuplot_.setYlabel("mole fraction mol/mol");
+ gnuplot_.addDataSetToPlot(x_, y1_, "N2_left.dat", "w l t 'N_2 left'");
+ gnuplot_.addDataSetToPlot(x_, y2_, "N2_right.dat", "w l t 'N_2 right'");
+ gnuplot_.plot("mole_fraction_N2");
+
+ gnuplot2_.resetPlot();
+ gnuplot2_.setXRange(0, std::max(time, 72000.0));
+ gnuplot2_.setYRange(0.0, 0.6);
+ gnuplot2_.setXlabel("time [s]");
+ gnuplot2_.setYlabel("mole fraction mol/mol");
+ gnuplot2_.addDataSetToPlot(x_, y3_, "CO2_left.dat", "w l t 'CO_2 left'");
+ gnuplot2_.addDataSetToPlot(x_, y4_, "CO2_right.dat", "w l t 'CO_2 right'");
+ gnuplot2_.plot("mole_fraction_C02");
+
+ gnuplot3_.resetPlot();
+ gnuplot3_.setXRange(0, std::max(time, 72000.0));
+ gnuplot3_.setYRange(0.0, 0.6);
+ gnuplot3_.setXlabel("time [s]");
+ gnuplot3_.setYlabel("mole fraction mol/mol");
+ gnuplot3_.addDataSetToPlot(x_, y5_, "H2_left.dat", "w l t 'H_2 left'");
+ gnuplot3_.addDataSetToPlot(x_, y6_, "H2_right.dat", "w l t 'H_2 right'");
+ gnuplot3_.plot("mole_fraction_H2");
+ }
+
+private:
+ GnuplotInterface<Scalar> gnuplot_;
+ GnuplotInterface<Scalar> gnuplot2_;
+ GnuplotInterface<Scalar> gnuplot3_;
+
+ std::vector<Scalar> x_;
+ std::vector<Scalar> y1_;
+ std::vector<Scalar> y2_;
+ std::vector<Scalar> y3_;
+ std::vector<Scalar> y4_;
+ std::vector<Scalar> y5_;
+ std::vector<Scalar> y6_;
+};
+
+}
+
+
+
+
int main(int argc, char** argv)
{
using namespace Dumux;
// define the type tag for this problem
- using TypeTag = Properties::TTag::MaxwellStefanNCTest;
+ using MomentumTypeTag = Properties::TTag::MaxwellStefanTestMomentum;
+ using MassTypeTag = Properties::TTag::MaxwellStefanTestMass;
// maybe initialize MPI and/or multithreading backend
initialize(argc, argv);
@@ -54,7 +182,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();
////////////////////////////////////////////////////////////
@@ -65,15 +193,24 @@ 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 (initial and 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 (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 tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
@@ -83,39 +220,49 @@ int main(int argc, char** argv)
timeLoop->setMaxTimeStepSize(maxDt);
// 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);
+ momentumGridVariables->init(x[momentumIdx]);
+ massGridVariables->init(x[massIdx]);
// 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
+ 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 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::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 = Dumux::MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
//! set some check points for the time loop
timeLoop->setPeriodicCheckPoint(tEnd/5.0);
+ PlotConcentration<Scalar> plotConcentration;
+
// time loop
timeLoop->start(); do
{
@@ -124,8 +271,12 @@ int main(int argc, char** argv)
// make the new solution the old solution
xOld = x;
- gridVariables->advanceTimeStep();
- problem->plotComponentsOverTime(x, *gridVariables,timeLoop->time()+timeLoop->timeStepSize());
+ momentumGridVariables->advanceTimeStep();
+ massGridVariables->advanceTimeStep();
+
+ static const bool plotOutput = getParam<bool>("Problem.PlotOutput");
+ if (plotOutput)
+ plotConcentration.plotComponentsOverTime(x[massIdx], *massGridVariables, timeLoop->time()+timeLoop->timeStepSize());
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
diff --git a/test/freeflow/navierstokesnc/maxwellstefan/params.input b/test/freeflow/navierstokesnc/maxwellstefan/params.input
index 3bc552f490..f406f1b138 100644
--- a/test/freeflow/navierstokesnc/maxwellstefan/params.input
+++ b/test/freeflow/navierstokesnc/maxwellstefan/params.input
@@ -11,9 +11,10 @@ Cells = 30 30
Name = test_msfreeflow
EnableGravity = false
EnableInertiaTerms = false
+PlotOutput = true
[Assembly]
-NumericDifference.BaseEpsilon = 1e-8
+NumericDifference.BaseEpsilon = 1e-4
[Vtk]
WriteFaceData = false
diff --git a/test/freeflow/navierstokesnc/maxwellstefan/problem.hh b/test/freeflow/navierstokesnc/maxwellstefan/problem.hh
index ef12f3bebf..f1bdfd7f54 100644
--- a/test/freeflow/navierstokesnc/maxwellstefan/problem.hh
+++ b/test/freeflow/navierstokesnc/maxwellstefan/problem.hh
@@ -15,13 +15,12 @@
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
-#include <dumux/common/numeqvector.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>
-#include <dumux/io/gnuplotinterface.hh>
-
namespace Dumux {
/*!
@@ -33,211 +32,151 @@ class MaxwellStefanNCTestProblem : public NavierStokesStaggeredProblem<TypeTag>
{
using ParentType = NavierStokesStaggeredProblem<TypeTag>;
- using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
- using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using BoundaryTypes = typename ParentType::BoundaryTypes;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
- using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
- using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
+ 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 SubControlVolume = typename GridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using FVElementGeometry = typename GridGeometry::LocalView;
- using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
+ static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
+ using Element = typename FVElementGeometry::Element;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
- static constexpr auto compOneIdx = Indices::conti0EqIdx;
- static constexpr auto compTwoIdx = Indices::conti0EqIdx + FluidSystem::N2Idx;
- static constexpr auto compThreeIdx = Indices::conti0EqIdx + FluidSystem::CO2Idx;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
-public:
- MaxwellStefanNCTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
- : ParentType(gridGeometry)
- {
- plotOutput_ = getParam<bool>("Problem.PlotOutput", false);
- }
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+ using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
- /*!
- * \name Problem parameters
- */
- // \{
+public:
+ MaxwellStefanNCTestProblem(std::shared_ptr<const GridGeometry> gridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(gridGeometry, couplingManager)
+ {}
/*!
- * \brief Writes out the diffusion rates from left to right
- *
- * Called after every time step
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary control volume.
*
- * \param curSol Vector containing the current solution
- * \param gridVariables The grid variables
- * \param time The time
+ * \param globalPos The position of the center of the finite volume
*/
- template<class SolutionVector, class GridVariables>
- void plotComponentsOverTime(const SolutionVector& curSol,
- const GridVariables& gridVariables,
- const Scalar time)
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
- if (plotOutput_)
+ BoundaryTypes values;
+
+ if constexpr (ParentType::isMomentumProblem())
{
- Scalar x_co2_left = 0.0;
- Scalar x_n2_left = 0.0;
- Scalar x_co2_right = 0.0;
- Scalar x_n2_right = 0.0;
- Scalar x_h2_left = 0.0;
- Scalar x_h2_right = 0.0;
- Scalar i = 0.0;
- Scalar j = 0.0;
- auto fvGeometry = localView(this->gridGeometry());
- auto elemVolVars = localView(gridVariables.curGridVolVars());
- for (const auto& element : elements(this->gridGeometry().gridView()))
- {
- fvGeometry.bindElement(element);
- elemVolVars.bind(element, fvGeometry, curSol);
- for (auto&& scv : scvs(fvGeometry))
- {
- const auto globalPos = scv.dofPosition();
-
- if (globalPos[0] < 0.5)
- {
- x_co2_left += elemVolVars[scv].moleFraction(FluidSystem::CO2Idx);
- x_n2_left += elemVolVars[scv].moleFraction(FluidSystem::N2Idx);
- x_h2_left += elemVolVars[scv].moleFraction(FluidSystem::H2Idx);
- i +=1;
- }
- else
- {
- x_co2_right += elemVolVars[scv].moleFraction(FluidSystem::CO2Idx);
- x_n2_right += elemVolVars[scv].moleFraction(FluidSystem::N2Idx);
- x_h2_right += elemVolVars[scv].moleFraction(FluidSystem::H2Idx);
- j +=1;
- }
- }
- }
- x_co2_left /= i;
- x_n2_left /= i;
- x_h2_left /= i;
- x_co2_right /= j;
- x_n2_right /= j;
- x_h2_right /= j;
-
- //do a gnuplot
- x_.push_back(time); // in seconds
- y1_.push_back(x_n2_left);
- y2_.push_back(x_n2_right);
- y3_.push_back(x_co2_left);
- y4_.push_back(x_co2_right);
- y5_.push_back(x_h2_left);
- y6_.push_back(x_h2_right);
-
- gnuplot_.resetPlot();
- gnuplot_.setXRange(0, std::max(time, 72000.0));
- gnuplot_.setYRange(0.4, 0.6);
- gnuplot_.setXlabel("time [s]");
- gnuplot_.setYlabel("mole fraction mol/mol");
- gnuplot_.addDataSetToPlot(x_, y1_, "N2_left.dat", "w l t 'N_2 left'");
- gnuplot_.addDataSetToPlot(x_, y2_, "N2_right.dat", "w l t 'N_2 right'");
- gnuplot_.plot("mole_fraction_N2");
-
- gnuplot2_.resetPlot();
- gnuplot2_.setXRange(0, std::max(time, 72000.0));
- gnuplot2_.setYRange(0.0, 0.6);
- gnuplot2_.setXlabel("time [s]");
- gnuplot2_.setYlabel("mole fraction mol/mol");
- gnuplot2_.addDataSetToPlot(x_, y3_, "CO2_left.dat", "w l t 'CO_2 left'");
- gnuplot2_.addDataSetToPlot(x_, y4_, "CO2_right.dat", "w l t 'CO_2 right'");
- gnuplot2_.plot("mole_fraction_C02");
-
- gnuplot3_.resetPlot();
- gnuplot3_.setXRange(0, std::max(time, 72000.0));
- gnuplot3_.setYRange(0.0, 0.6);
- gnuplot3_.setXlabel("time [s]");
- gnuplot3_.setYlabel("mole fraction mol/mol");
- gnuplot3_.addDataSetToPlot(x_, y5_, "H2_left.dat", "w l t 'H_2 left'");
- gnuplot3_.addDataSetToPlot(x_, y6_, "H2_right.dat", "w l t 'H_2 right'");
- gnuplot3_.plot("mole_fraction_H2");
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
}
+ else
+ values.setAllNeumann();
+
+ return values;
}
- // \}
+ //! Enable internal Dirichlet constraints
+ static constexpr bool enableInternalDirichletConstraints()
+ { return !ParentType::isMomentumProblem(); }
/*!
- * \name Boundary conditions
- */
- // \{
-
- /*!
- * \brief Specifies which kind of boundary condition should be
- * used for which equation on a given boundary control volume.
+ * \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 globalPos The position of the center of the finite volume
+ * \param element The finite element
+ * \param scv The sub-control volume
*/
- BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+ std::bitset<DirichletValues::dimension> hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const
{
- BoundaryTypes values;
- // Set no-flow/no-slip everywhere. Neumann defaults to zero.
- values.setDirichlet(Indices::velocityXIdx);
- values.setDirichlet(Indices::velocityYIdx);
- values.setNeumann(compOneIdx);
- values.setNeumann(compTwoIdx);
- values.setNeumann(compThreeIdx);
+ std::bitset<DirichletValues::dimension> values;
+
+ // the pure Neumann problem is only defined up to a constant
+ // we create a well-posed problem by fixing the pressure at one dof
+ if (scv.dofIndex() == 0)
+ values.set(0);
+
return values;
}
- /*!
- * \brief Evaluates the boundary conditions for a Dirichlet control volume.
+ /*!
+ * \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); }
+
+ /*!
+ * \brief Returns Dirichlet boundary values at a given position.
*
- * \param globalPos The center of the finite volume which ought to be set.
+ * \param globalPos The global position
*/
- PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
- {
- return initialAtPos(globalPos);
- }
+ DirichletValues dirichletAtPos(const GlobalPosition &globalPos) const
+ { return initialAtPos(globalPos); }
- /*!
- * \name Volume terms
+ /*!
+ * \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>
+ BoundaryFluxes neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
+ { return BoundaryFluxes(0.0); }
- /*!
+ /*!
* \brief Evaluates the initial value for a control volume.
*
* \param globalPos The global position
*/
- PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ InitialValues initialAtPos(const GlobalPosition &globalPos) const
{
- PrimaryVariables initialValues(0.0);
- if (globalPos[0] < 0.5)
+ InitialValues initialValues(0.0);
+
+ if constexpr (!ParentType::isMomentumProblem())
{
- initialValues[compTwoIdx] = 0.50086;
- initialValues[compThreeIdx] = 0.49914;
+ static constexpr auto compTwoIdx = Indices::conti0EqIdx + FluidSystem::N2Idx;
+ static constexpr auto compThreeIdx = Indices::conti0EqIdx + FluidSystem::CO2Idx;
+
+ initialValues[Indices::pressureIdx] = 1.1e+5;
+ if (globalPos[0] < 0.5)
+ {
+ initialValues[compTwoIdx] = 0.50086;
+ initialValues[compThreeIdx] = 0.49914;
+ }
+ else
+ {
+ initialValues[compTwoIdx] = 0.49879;
+ initialValues[compThreeIdx] = 0.0;
+ }
}
else
{
- initialValues[compTwoIdx] = 0.49879;
- initialValues[compThreeIdx] = 0.0;
+ initialValues[Indices::velocityXIdx] = 0.0;
+ initialValues[Indices::velocityYIdx] = 0.0;
}
- initialValues[Indices::pressureIdx] = 1.1e+5;
- initialValues[Indices::velocityXIdx] = 0.0;
- initialValues[Indices::velocityYIdx] = 0.0;
-
return initialValues;
}
-private:
-
- bool plotOutput_;
+private:
const Scalar eps_{1e-6};
-
- Dumux::GnuplotInterface<Scalar> gnuplot_;
- Dumux::GnuplotInterface<Scalar> gnuplot2_;
- Dumux::GnuplotInterface<Scalar> gnuplot3_;
-
- std::vector<Scalar> x_;
- std::vector<Scalar> y1_;
- std::vector<Scalar> y2_;
- std::vector<Scalar> y3_;
- std::vector<Scalar> y4_;
- std::vector<Scalar> y5_;
- std::vector<Scalar> y6_;
};
} // end namespace Dumux
diff --git a/test/freeflow/navierstokesnc/maxwellstefan/properties.hh b/test/freeflow/navierstokesnc/maxwellstefan/properties.hh
index 63980d3202..35b5f5a734 100644
--- a/test/freeflow/navierstokesnc/maxwellstefan/properties.hh
+++ b/test/freeflow/navierstokesnc/maxwellstefan/properties.hh
@@ -14,13 +14,16 @@
#include <dune/grid/yaspgrid.hh>
-#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/freeflow/navierstokes/momentum/model.hh>
+#include <dumux/freeflow/navierstokes/mass/1pnc/model.hh>
#include <dumux/flux/maxwellstefanslaw.hh>
+#include <dumux/material/fluidsystems/base.hh>
-#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
-#include <dumux/material/components/simpleh2o.hh>
-#include <dumux/material/fluidsystems/h2oair.hh>
+#include <dumux/discretization/fcstaggered.hh>
+#include <dumux/discretization/cctpfa.hh>
+
+#include <dumux/multidomain/staggeredfreeflow/couplingmanager.hh>
#include "problem.hh"
@@ -28,7 +31,9 @@ namespace Dumux::Properties {
// Create new type tags
namespace TTag {
-struct MaxwellStefanNCTest { using InheritsFrom = std::tuple<NavierStokesNC, StaggeredFreeFlowModel>; };
+struct MaxwellStefanNCTest {};
+struct MaxwellStefanTestMomentum { using InheritsFrom = std::tuple<MaxwellStefanNCTest, NavierStokesMomentum, FaceCenteredStaggeredModel>; };
+struct MaxwellStefanTestMass { using InheritsFrom = std::tuple<MaxwellStefanNCTest, NavierStokesMassOnePNC, CCTpfaModel>; };
} // end namespace TTag
template<class TypeTag>
@@ -56,6 +61,12 @@ 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
@@ -91,6 +102,13 @@ public:
static Scalar molarMass(unsigned int compIdx)
{ return 0.02896; }
+ //! Returns whether the fluids are miscible
+ static constexpr bool isMiscible()
+ { return false; }
+
+ //! Returns whether the fluids are compressible
+ static constexpr bool isCompressible(int phaseIdx = 0)
+ { return false; }
using Base::binaryDiffusionCoefficient;
/*!
--
GitLab