diff --git a/test/multidomain/boundary/CMakeLists.txt b/test/multidomain/boundary/CMakeLists.txt
index addd460feb4d2b502e6a92a0630ad445b11020d9..57ed3d76834b0045aa3e70443493256df325bd9b 100644
--- a/test/multidomain/boundary/CMakeLists.txt
+++ b/test/multidomain/boundary/CMakeLists.txt
@@ -1,2 +1,3 @@
add_subdirectory(darcydarcy)
add_subdirectory(stokesdarcy)
+add_subdirectory(ransdarcy)
diff --git a/test/multidomain/boundary/ransdarcy/1p2c_2p2c/CMakeLists.txt b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..061b603d5233c61e3484b93508b058fae38718ad
--- /dev/null
+++ b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/CMakeLists.txt
@@ -0,0 +1,72 @@
+add_input_file_links()
+dune_symlink_to_source_files(FILES scripts)
+
+dune_add_test(NAME test_md_boundary_darcy2p2cni_rans1p2cnikepsilon
+ SOURCES main.cc
+ COMPILE_DEFINITIONS RANSTYPETAG=RANSSubProblemKEpsilonNCNI
+ LABELS freeflow rans multidomain
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy2p2cni_rans1p2cnikepsilon_rans-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnikepsilon_rans-00044.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy2p2cni_rans1p2cnikepsilon_darcy-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnikepsilon_darcy-00044.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnikepsilon params.input
+ -Vtk.OutputName test_md_boundary_darcy2p2cni_rans1p2cnikepsilon")
+
+dune_add_test(NAME test_md_boundary_darcy2p2cni_rans1p2cnikomega
+ SOURCES main.cc
+ COMPILE_DEFINITIONS RANSTYPETAG=RANSSubProblemKOmegaNCNI
+ LABELS freeflow rans multidomain
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy2p2cni_rans1p2cnikomega_rans-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnikomega_rans-00044.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy2p2cni_rans1p2cnikomega_darcy-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnikomega_darcy-00044.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnikomega params.input
+ -Vtk.OutputName test_md_boundary_darcy2p2cni_rans1p2cnikomega")
+
+dune_add_test(NAME test_md_boundary_darcy2p2cni_rans1p2cnilowrekepsilon
+ SOURCES main.cc
+ COMPILE_DEFINITIONS RANSTYPETAG=RANSSubProblemLowReKEpsilonNCNI
+ LABELS freeflow rans multidomain
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy2p2cni_rans1p2cnilowrekepsilon_rans-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnilowrekepsilon_rans-00041.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy2p2cni_rans1p2cnilowrekepsilon_darcy-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnilowrekepsilon_darcy-00041.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnilowrekepsilon params.input
+ -Vtk.OutputName test_md_boundary_darcy2p2cni_rans1p2cnilowrekepsilon")
+
+dune_add_test(NAME test_md_boundary_darcy2p2cni_rans1p2cnioneeq
+ SOURCES main.cc
+ COMPILE_DEFINITIONS RANSTYPETAG=RANSSubProblemOneEqNCNI
+ LABELS freeflow rans multidomain
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy2p2cni_rans1p2cnioneeq_rans-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnioneeq_rans-00041.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy2p2cni_rans1p2cnioneeq_darcy-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnioneeq_darcy-00041.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnioneeq params.input
+ -Vtk.OutputName test_md_boundary_darcy2p2cni_rans1p2cnioneeq")
+
+dune_add_test(NAME test_md_boundary_darcy2p2cni_rans1p2cnizeroeq
+ SOURCES main.cc
+ COMPILE_DEFINITIONS RANSTYPETAG=RANSSubProblemZeroEqNCNI
+ LABELS freeflow rans multidomain
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy2p2cni_rans1p2cnizeroeq_rans-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnizeroeq_rans-00040.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy2p2cni_rans1p2cnizeroeq_darcy-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnizeroeq_darcy-00040.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_rans1p2cnizeroeq params.input
+ -Vtk.OutputName test_md_boundary_darcy2p2cni_rans1p2cnizeroeq")
diff --git a/test/multidomain/boundary/ransdarcy/1p2c_2p2c/main.cc b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/main.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c9417fdd43b679c7f0c9eb2b8e1a0f4009250a9e
--- /dev/null
+++ b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/main.cc
@@ -0,0 +1,288 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A test problem for the nonisothermal compositional coupled RANS/Darcy problem (1p2cni/2p2cni)
+ */
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+
+#include "problem_darcy.hh"
+#include "problem_rans.hh"
+
+namespace Dumux {
+namespace Properties {
+
+template
+struct CouplingManager
+{
+ using Traits = StaggeredMultiDomainTraits;
+ using type = Dumux::StokesDarcyCouplingManager;
+};
+
+template
+struct CouplingManager
+{
+ using Traits = StaggeredMultiDomainTraits;
+ using type = Dumux::StokesDarcyCouplingManager;
+};
+
+} // end namespace Properties
+} // end namespace Dumux
+
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // Define the sub problem type tags
+ using RANSTypeTag = Properties::TTag::RANSTYPETAG;
+ using DarcyTypeTag = Properties::TTag::DarcyTwoPTwoCNI;
+
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using RANSGridManager = Dumux::GridManager>;
+ RANSGridManager ransGridManager;
+ ransGridManager.init("RANS"); // pass parameter group
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridManager.grid().leafGridView();
+ const auto& ransGridView = ransGridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using RANSFVGridGeometry = GetPropType;
+ auto ransFvGridGeometry = std::make_shared(ransGridView);
+ ransFvGridGeometry->update();
+ using DarcyFVGridGeometry = GetPropType;
+ auto darcyFvGridGeometry = std::make_shared(darcyGridView);
+ darcyFvGridGeometry->update();
+
+ using Traits = StaggeredMultiDomainTraits;
+
+ // the coupling manager
+ using CouplingManager = StokesDarcyCouplingManager;
+ auto couplingManager = std::make_shared(ransFvGridGeometry, darcyFvGridGeometry);
+
+ // the indices
+ constexpr auto ransCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+ constexpr auto ransFaceIdx = CouplingManager::stokesFaceIdx;
+ constexpr auto darcyIdx = CouplingManager::darcyIdx;
+
+ // initialize the fluidsystem (tabulation)
+ GetPropType::init();
+
+ // the problem (initial and boundary conditions)
+ using RANSProblem = GetPropType;
+ auto ransProblem = std::make_shared(ransFvGridGeometry, couplingManager);
+ using DarcyProblem = GetPropType;
+ auto darcyProblem = std::make_shared(darcyFvGridGeometry, couplingManager);
+
+ // get some time loop parameters
+ using Scalar = GetPropType;
+ const auto tEnd = getParam("TimeLoop.TEnd");
+ const auto maxDt = getParam("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam("TimeLoop.DtInitial");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared>(0.0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // set timeloop for the subproblems, needed for boundary value variations
+ ransProblem->setTimeLoop(timeLoop);
+ darcyProblem->setTimeLoop(timeLoop);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[ransCellCenterIdx].resize(ransFvGridGeometry->numCellCenterDofs());
+ sol[ransFaceIdx].resize(ransFvGridGeometry->numFaceDofs());
+ sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+
+ // get a solution vector storing references to the two Stokes solution vectors
+ auto ransSol = partial(sol, ransCellCenterIdx, ransFaceIdx);
+
+ // apply initial solution for instationary problems
+ ransProblem->applyInitialSolution(ransSol);
+ ransProblem->updateStaticWallProperties();
+ ransProblem->updateDynamicWallProperties(ransSol);
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+
+ auto solOld = sol;
+
+ couplingManager->init(ransProblem, darcyProblem, sol);
+
+ // the grid variables
+ using RANSGridVariables = GetPropType;
+ auto ransGridVariables = std::make_shared(ransProblem, ransFvGridGeometry);
+ ransGridVariables->init(ransSol);
+ using DarcyGridVariables = GetPropType;
+ auto darcyGridVariables = std::make_shared(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx]);
+ darcyProblem->init(sol[darcyIdx], *darcyGridVariables, ransProblem->ransName());
+
+ // intialize the vtk output module
+ StaggeredVtkOutputModule ransVtkWriter(*ransGridVariables, ransSol, ransProblem->name());
+ GetPropType::initOutputModule(ransVtkWriter);
+ ransVtkWriter.write(0.0);
+
+ VtkOutputModule> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyProblem->name());
+ using DarcyVelocityOutput = GetPropType;
+ darcyVtkWriter.addVelocityOutput(std::make_shared(*darcyGridVariables));
+ GetPropType::initOutputModule(darcyVtkWriter);
+ darcyVtkWriter.write(0.0);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = MultiDomainFVAssembler;
+ auto assembler = std::make_shared(std::make_tuple(ransProblem, ransProblem, darcyProblem),
+ std::make_tuple(ransFvGridGeometry->cellCenterFVGridGeometryPtr(),
+ ransFvGridGeometry->faceFVGridGeometryPtr(),
+ darcyFvGridGeometry),
+ std::make_tuple(ransGridVariables->cellCenterGridVariablesPtr(),
+ ransGridVariables->faceGridVariablesPtr(),
+ darcyGridVariables),
+ couplingManager,
+ timeLoop);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared();
+
+ // the non-linear solver
+ using NewtonSolver = MultiDomainNewtonSolver;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(solOld);
+
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(sol, *timeLoop);
+
+ // make the new solution the old solution
+ solOld = sol;
+
+ // post time step treatment of Darcy problem
+ darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables);
+
+ // advance grid variables to the next time step
+ ransGridVariables->advanceTimeStep();
+ darcyGridVariables->advanceTimeStep();
+
+ ransSol[ransCellCenterIdx] = sol[ransCellCenterIdx];
+ ransSol[ransFaceIdx] = sol[ransFaceIdx];
+
+ // update the auxiliary free flow solution vector
+ ransProblem->updateDynamicWallProperties(ransSol);
+ assembler->updateGridVariables(sol);
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+
+ // write vtk output
+ ransVtkWriter.write(timeLoop->time());
+ darcyVtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(ransGridView.comm());
+ timeLoop->finalize(darcyGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/test/multidomain/boundary/ransdarcy/1p2c_2p2c/params.input b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/params.input
new file mode 100644
index 0000000000000000000000000000000000000000..8e8068d869a11f334090cfb84000d3d5b6acde7a
--- /dev/null
+++ b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/params.input
@@ -0,0 +1,79 @@
+[TimeLoop]
+DtInitial = 1e-3 # [s]
+MaxTimeStepSize = 20 # [s]
+TEnd = 200 # [s]
+
+[RANS.Grid]
+Verbosity = true
+Positions0 = -1.0 0.0 1.0 2.0
+Positions1 = 0.5 1.0
+Grading0 = -1.0 1.0 1.0
+Grading1 = 1.3
+Cells0 = 5 10 5
+Cells1 = 15
+
+[Darcy.Grid]
+Verbosity = true
+Positions0 = 0.0 1.0
+Positions1 = 0.0 0.5
+Cells0 = 10
+Cells1 = 15
+Grading0 = 1.0
+Grading1 = -1.3
+
+[RANS.Problem]
+Name = rans
+InletVelocity = 2.5 # [m/s]
+InletPressure = 1e5 # [Pa]
+InletMassFrac = 0.008 # [-]
+Temperature = 298.15 # [K]
+
+[Darcy.Problem]
+Name = darcy
+Pressure = 1e5 # [Pa]
+Saturation = 0.8 # initial Sw
+Temperature = 298.15 # [K]
+InitPhasePresence = 3 # bothPhases
+InitialSaturation = 0.85 # [-]
+
+[Darcy.SpatialParams]
+Porosity = 0.41
+Permeability = 2.65e-10
+AlphaBeaversJoseph = 1.0
+Swr = 0.005
+Snr = 0.01
+VgAlpha = 6.37e-4
+VgN = 8.0
+
+[SpatialParams]
+ForchCoeff = 0.55 #
+
+[Problem]
+EnableGravity = true
+InterfaceDiffusionCoefficientAvg = FreeFlowOnly # Harmonic
+
+[RANS.RANS]
+UseStoredEddyViscosity = false
+EddyViscosityModel = "baldwinLomax"
+
+[RANS.KEpsilon]
+YPlusThreshold = 40. # should be small (10-30) for coarse grids large (30-60) for fine grids
+
+[Vtk]
+AddVelocity = true
+WriteFaceData = false
+OutputName = rectangle_interface
+
+[Component]
+SolidDensity = 2700
+SolidThermalConductivity = 2.8
+SolidHeatCapacity = 790
+
+[Newton]
+TargetSteps = 8
+MaxSteps = 10
+MaxRelativeShift = 1e-8
+
+[Assembly]
+NumericDifferenceMethod = 0
+NumericDifference.BaseEpsilon = 1e-8
diff --git a/test/multidomain/boundary/ransdarcy/1p2c_2p2c/problem_darcy.hh b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/problem_darcy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ce3d20c2eb527e205ff3d6871c0c53cfb6b5557b
--- /dev/null
+++ b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/problem_darcy.hh
@@ -0,0 +1,408 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief A Darcy Subproblem (cell-centered finite volume method).
+ */
+
+#ifndef DUMUX_DARCY2P2C_SUBPROBLEM_HH
+#define DUMUX_DARCY2P2C_SUBPROBLEM_HH
+
+#include
+
+#include
+#include
+
+#include
+#include
+
+#include
+#include
+
+#include "spatialparams.hh"
+
+namespace Dumux {
+template
+class DarcySubProblem;
+
+namespace Properties {
+// Create new type tags
+namespace TTag {
+struct DarcyTwoPTwoCNI { using InheritsFrom = std::tuple; };
+} // end namespace TTag
+
+// Set the problem property
+template
+struct Problem { using type = Dumux::DarcySubProblem; };
+
+// the fluid system
+template
+struct FluidSystem { using type = FluidSystems::H2OAir>; };
+
+template
+struct ReplaceCompEqIdx { static constexpr int value = 10; };
+
+// Set the grid type
+template
+struct Grid
+{ using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates, 2> >; };
+
+template
+struct SpatialParams
+{
+ using FVGridGeometry = GetPropType;
+ using Scalar = GetPropType;
+ using type = TwoPTwoCSpatialParams;
+};
+
+template
+struct UseMoles { static constexpr bool value = true; };
+
+//! Set the default formulation to pw-Sn: This can be over written in the problem.
+template
+struct Formulation
+{ static constexpr auto value = TwoPFormulation::p1s0; };
+
+} // end namespace Properties
+
+template
+class DarcySubProblem : public PorousMediumFlowProblem
+{
+ using ParentType = PorousMediumFlowProblem;
+ using GridView = GetPropType;
+ using Scalar = GetPropType;
+ using PrimaryVariables = GetPropType;
+ using NumEqVector = GetPropType;
+ using BoundaryTypes = GetPropType;
+ using VolumeVariables = GetPropType;
+ using FVElementGeometry = typename GetPropType::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = GetPropType;
+
+ using ElementVolumeVariables = typename GetPropType::LocalView;
+
+ using FluidSystem = GetPropType;
+ using FluidState = GetPropType;
+
+ using Indices = typename GetPropType::Indices;
+
+ enum {
+ // grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ conti0EqIdx = Indices::conti0EqIdx,
+ contiWEqIdx = Indices::conti0EqIdx + FluidSystem::H2OIdx,
+ contiNEqIdx = Indices::conti0EqIdx + FluidSystem::AirIdx,
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx
+ };
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using CouplingManager = GetPropType;
+ using TimeLoopPtr = std::shared_ptr>;
+
+ using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
+
+
+public:
+ DarcySubProblem(std::shared_ptr fvGridGeometry,
+ std::shared_ptr couplingManager)
+ : ParentType(fvGridGeometry, "Darcy")
+ , couplingManager_(couplingManager)
+ , eps_(1e-7)
+
+ {
+ pressure_ = getParamFromGroup(this->paramGroup(), "Problem.Pressure");
+ initialSaturation_ = getParamFromGroup(this->paramGroup(), "Problem.InitialSaturation");
+ temperature_ = getParamFromGroup(this->paramGroup(), "Problem.Temperature");
+ initialPhasePresence_ = getParamFromGroup(this->paramGroup(), "Problem.InitPhasePresence");
+
+ diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{},
+ getParamFromGroup(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg"));
+
+ problemName_ = getParam("Vtk.OutputName") + "_" + getParamFromGroup(this->paramGroup(), "Problem.Name");
+ }
+
+
+ /*!
+ * \name Mass Exchange Evaluations
+ */
+ // \{
+
+ /*!
+ * \brief Initialize the problem.
+ */
+ template
+ void init(const SolutionVector& curSol,
+ const GridVariables& gridVariables,
+ const std::string ransName)
+ {
+ initialWaterContent_ = evaluateWaterMassStorageTerm(curSol, gridVariables);
+ lastWaterMass_ = initialWaterContent_;
+
+ storageFileName_ = "storage_" + problemName_ + "_" + ransName + ".csv";
+ storageFile_.open(storageFileName_, std::ios::app);
+ storageFile_ << "#Time[s]" << ";"
+ << "WaterMass[kg]" << ";"
+ << "WaterMassLoss[kg]" << ";"
+ << "EvaporationRate[mm/d]"
+ << std::endl;
+ storageFile_.close();
+ }
+
+ /*!
+ * \brief Evaluate the exchange processes after each timestep.
+ */
+ template
+ void postTimeStep(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+
+ {
+ evaluateWaterMassStorageTerm(curSol, gridVariables);
+ }
+
+ /*!
+ * \brief Calculate the total water mass storage exchange.
+ */
+ template
+ Scalar evaluateWaterMassStorageTerm(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+ {
+ // compute the mass in the entire domain
+ Scalar waterMass = 0.0;
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(gridVariables.curGridVolVars());
+ elemVolVars.bindElement(element, fvGeometry, curSol);
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ for(int phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
+ {
+ // insert calculation of the water mass here
+ waterMass += volVars.massFraction(phaseIdx, FluidSystem::H2OIdx) * volVars.density(phaseIdx)
+ * volVars.saturation(phaseIdx) * volVars.porosity()
+ * scv.volume() * volVars.extrusionFactor();
+ }
+ }
+ }
+ Scalar cumMassLoss = initialWaterContent_ - waterMass;
+ Scalar evaporationRate = (lastWaterMass_ - waterMass) * 86400
+ / (this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0])
+ / timeLoop_->timeStepSize();
+ lastWaterMass_ = waterMass;
+ storageFile_.open(this->storageFileName_, std::ios::app);
+ storageFile_ << timeLoop_->time() << ";"
+ << waterMass << ";"
+ << cumMassLoss << ";"
+ << evaporationRate
+ << "\n";
+ storageFile_.close();
+
+ std::cout << "Mass of water is: " << waterMass << "\n";
+ std::cout << "Evaporation rate is: " << evaporationRate << "\n";
+ std::cout << "\n";
+ return waterMass;
+ }
+
+ // \}
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ */
+ const std::string& name() const
+ {
+ return problemName_;
+ }
+
+
+ /*!
+ * \brief Returns the temperature within the domain in [K].
+ *
+ */
+ Scalar temperature() const
+ { return temperature_; }
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary control volume.
+ *
+ * \param element The element
+ * \param scvf The boundary sub control volume face
+ */
+ BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ BoundaryTypes values;
+ values.setAllNeumann();
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values.setAllCouplingNeumann();
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate 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 scvf The boundary sub control volume face
+ *
+ */
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ {
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
+
+ for(int i = 0; i< massFlux.size(); ++i)
+ values[i] = massFlux[i];
+
+ values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(element, fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
+ }
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * \param element The element for which the source term is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scv The subcontrolvolume
+ */
+ template
+ NumEqVector source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ { return NumEqVector(0.0); }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * For this method, the \a priVars parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ values.setState(Indices::bothPhases);
+
+ FluidState fluidState;
+ fluidState.setPressure(0, pressure_);
+ fluidState.setTemperature(temperature());
+ fluidState.setMoleFraction(0, 0, 1);
+ Scalar density = FluidSystem::density(fluidState, 0);
+
+ values[Indices::pressureIdx] = pressure_ - density * this->gravity()[dimWorld-1]* (this->fvGridGeometry().bBoxMax()[1] - globalPos[1]);
+ values[switchIdx] = initialSaturation_;
+ values[Indices::temperatureIdx] = temperature();
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \brief Set the coupling manager
+ */
+ void setCouplingManager(std::shared_ptr cm)
+ { couplingManager_ = cm; }
+
+ /*!
+ * \brief Get the coupling manager
+ */
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ { timeLoop_ = timeLoop; }
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ Scalar pressure_;
+ Scalar temperature_;
+ int initialPhasePresence_;
+ Scalar initialSaturation_;
+
+ TimeLoopPtr timeLoop_;
+ std::shared_ptr couplingManager_;
+ DiffusionCoefficientAveragingType diffCoeffAvgType_;
+ Scalar eps_;
+
+ std::string problemName_;
+ std::string storageFileName_;
+ std::ofstream storageFile_;
+
+ bool plotFluxes_;
+ bool plotStorage_;
+ Scalar initialWaterContent_ = 0.0;
+ Scalar lastWaterMass_ = 0.0;
+
+};
+} // end namespace Dumux
+
+#endif //DUMUX_DARCY_SUBPROBLEM_HH
diff --git a/test/multidomain/boundary/ransdarcy/1p2c_2p2c/problem_rans.hh b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/problem_rans.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5af2bc743e9fa9cfb441f113c46241b7cf570b5c
--- /dev/null
+++ b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/problem_rans.hh
@@ -0,0 +1,668 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief The free-flow RANS subproblem (staggered grid model).
+ */
+
+#ifndef DUMUX_RANS_SUBPROBLEM_HH
+#define DUMUX_RANS_SUBPROBLEM_HH
+
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+
+#include
+
+#include
+#include
+
+#include
+#include
+
+#include
+#include
+
+#include
+#include
+
+#include
+
+namespace Dumux {
+template
+class RANSSubProblem;
+
+namespace Properties {
+
+// Create new type tags
+namespace TTag {
+
+// Base Typetag
+struct RANSSubModel { using InheritsFrom = std::tuple; };
+
+// Typetags
+struct RANSSubProblemZeroEqNCNI { using InheritsFrom = std::tuple; };
+struct RANSSubProblemOneEqNCNI { using InheritsFrom = std::tuple; };
+struct RANSSubProblemKOmegaNCNI { using InheritsFrom = std::tuple; };
+struct RANSSubProblemLowReKEpsilonNCNI { using InheritsFrom = std::tuple; };
+struct RANSSubProblemKEpsilonNCNI { using InheritsFrom = std::tuple; };
+
+} // end namespace TTag
+
+// The fluid system
+template
+struct FluidSystem
+{
+ using H2OAir = FluidSystems::H2OAir>;
+ static constexpr auto phaseIdx = H2OAir::gasPhaseIdx;
+ using type = FluidSystems::OnePAdapter;
+};
+
+// Set the grid type
+template
+struct Grid
+{ using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates, 2> >; };
+
+template
+struct ReplaceCompEqIdx { static constexpr int value = 10; };
+
+// Use formulation based on mole fractions
+template
+struct UseMoles { static constexpr bool value = true; };
+
+// Set the problem property
+template
+struct Problem { using type = Dumux::RANSSubProblem ; };
+
+template
+struct EnableFVGridGeometryCache { static constexpr bool value = true; };
+template
+struct EnableGridFluxVariablesCache { static constexpr bool value = true; };
+template
+struct EnableGridVolumeVariablesCache { static constexpr bool value = true; };
+} // end namespace Properties
+
+/*!
+ * \ingroup BoundaryTests
+ * \brief Test problem for the one-phase Reynolds Averaged (Navier-) Stokes problem.
+ *
+ * Horizontal flow from left to right with a turbulent velocity profile.
+ */
+template
+class RANSSubProblem : public RANSProblem
+{
+ using ParentType = RANSProblem;
+ using CouplingManager = GetPropType;
+
+ using BoundaryTypes = GetPropType;
+ using FluidSystem = GetPropType;
+ using FluidState = GetPropType;
+ using FVGridGeometry = GetPropType;
+ using Indices = typename GetPropType::Indices;
+ using NumEqVector = GetPropType;
+ using PrimaryVariables = GetPropType;
+ using Scalar = GetPropType;
+
+ using ModelTraits = GetPropType;
+ using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using FVElementGeometry = typename GetPropType::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using GridView = GetPropType;
+
+ using TimeLoopPtr = std::shared_ptr>;
+ using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
+
+ static constexpr auto transportEqIdx = Indices::conti0EqIdx + 1;
+ static constexpr auto transportCompIdx = Indices::conti0EqIdx + 1;
+ static constexpr auto dimWorld = FVGridGeometry::GridView::dimensionworld;
+
+public:
+ RANSSubProblem(std::shared_ptr fvGridGeometry, std::shared_ptr couplingManager)
+ : ParentType(fvGridGeometry, "RANS"), eps_(1e-6), couplingManager_(couplingManager)
+ {
+ pressure_ = getParamFromGroup(this->paramGroup(), "Problem.InletPressure");
+ inletMassFrac_ = getParamFromGroup(this->paramGroup(), "Problem.InletMassFrac");
+ temperature_ = getParamFromGroup(this->paramGroup(), "Problem.Temperature");
+ inletVelocity_ = getParamFromGroup(this->paramGroup(), "Problem.InletVelocity");
+
+ Dumux::TurbulenceProperties turbulenceProperties;
+ FluidState fluidState;
+ const auto phaseIdx = 0;
+ fluidState.setPressure(phaseIdx, 1e5);
+ fluidState.setTemperature(temperature());
+ fluidState.setMassFraction(phaseIdx, phaseIdx, inletMassFrac());
+ Scalar density = FluidSystem::density(fluidState, 0);
+ Scalar kinematicViscosity = FluidSystem::viscosity(fluidState, 0) / density;
+ Scalar charLength = this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1];
+ diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{},
+ getParam("Problem.InterfaceDiffusionCoefficientAvg"));
+
+ std::cout << "The reynolds Number (l_char = channel diameter) is: " << inletVelocity_ * charLength / kinematicViscosity << "\n";
+
+
+ // ideally the viscosityTilde parameter as inflow for the Spalart-Allmaras model should be zero
+ viscosityTilde_ = 1e-3 * turbulenceProperties.viscosityTilde(inletVelocity_, charLength, kinematicViscosity);
+ turbulentKineticEnergy_ = turbulenceProperties.turbulentKineticEnergy(inletVelocity_, charLength, kinematicViscosity);
+
+ if (ModelTraits::turbulenceModel() == TurbulenceModel::komega)
+ dissipation_ = turbulenceProperties.dissipationRate(inletVelocity_, charLength, kinematicViscosity);
+ else
+ dissipation_ = turbulenceProperties.dissipation(inletVelocity_, charLength, kinematicViscosity);
+ std::cout << std::endl;
+
+ problemName_ = getParam("Vtk.OutputName") + "_" + getParamFromGroup(this->paramGroup(), "Problem.Name");
+ turbulenceModelName_ = turbulenceModelToString(ModelTraits::turbulenceModel());
+ std::cout << "Using the "<< turbulenceModelName_ << " Turbulence Model. \n";
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ //! \brief The problem name.
+ const std::string& name() const
+ { return problemName_; }
+
+ //! \brief The name of the RANS model.
+ const std::string& ransName() const
+ { return turbulenceModelName_; }
+
+ //! \brief Return the temperature within the domain in [K].
+ const Scalar temperature() const
+ { return temperature_; }
+
+ //! \brief Returns the inlet velocity.
+ const Scalar inletVelocity() const
+ { return inletVelocity_ ;}
+
+ //! \brief Returns the inlet pressure.
+ const Scalar inletPressure() const
+ { return pressure_; }
+
+ //! \brief Returns the inlet mass fraction.
+ const Scalar inletMassFrac() const
+ { return inletMassFrac_; }
+
+ // \}
+
+ /*!
+ * \name Boundary Locations
+ */
+ // \{
+
+ /*!
+ * \brief Returns if the location is on the wall.
+ *
+ * \param scvf The sub control volume face
+ */
+ bool isOnWall(const SubControlVolumeFace &scvf) const
+ {
+ GlobalPosition globalPos = scvf.ipGlobal();
+ return isOnWallAtPos(globalPos);
+ }
+
+ /*!
+ * \brief Returns if the location is on the wall.
+ *
+ * \param globalPos The global position
+ */
+ bool isOnWallAtPos(const GlobalPosition &globalPos) const
+ {
+ return onLowerBoundary_(globalPos);
+ }
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ *
+ * \param element The finite element
+ * \param scvf The sub control volume face
+ */
+ BoundaryTypes boundaryTypes(const Element& element,
+ const SubControlVolumeFace& scvf) const
+ {
+ BoundaryTypes bTypes;
+ GlobalPosition globalPos = scvf.ipGlobal();
+
+ // turbulence model-specific boundary types
+ static constexpr auto numEq = numTurbulenceEq(ModelTraits::turbulenceModel());
+ setBcTypes_(bTypes, globalPos, Dune::index_constant{});
+
+ // common boundary types for all turbulence models
+ if(isInlet_(globalPos))
+ {
+ bTypes.setDirichlet(Indices::velocityXIdx);
+ bTypes.setDirichlet(Indices::velocityYIdx);
+ bTypes.setDirichlet(Indices::temperatureIdx);
+ bTypes.setDirichlet(transportCompIdx);
+ }
+ else if (isOutlet_(globalPos))
+ {
+ bTypes.setDirichlet(Indices::pressureIdx);
+ bTypes.setOutflow(transportEqIdx);
+ bTypes.setOutflow(Indices::energyEqIdx);
+ }
+ else if (isOnWall(scvf))
+ {
+ bTypes.setDirichlet(Indices::velocityXIdx);
+ bTypes.setDirichlet(Indices::velocityYIdx);
+ bTypes.setNeumann(Indices::energyEqIdx);
+ bTypes.setNeumann(Indices::conti0EqIdx + 1);
+ }
+ else
+ {
+ bTypes.setAllSymmetry();
+ }
+
+ if (isOnCouplingWall_(scvf))
+ {
+ bTypes.setCouplingNeumann(Indices::conti0EqIdx);
+ bTypes.setCouplingNeumann(Indices::conti0EqIdx + 1);
+ bTypes.setCouplingNeumann(Indices::energyEqIdx);
+ bTypes.setCouplingNeumann(scvf.directionIndex());
+ bTypes.setBJS(1 - scvf.directionIndex());
+ }
+
+
+ return bTypes;
+ }
+
+ /*!
+ * \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param scv The sub control volume
+ */
+ template
+ bool isDirichletCell(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolume& scv,
+ int pvIdx) const
+ {
+ using IsKOmegaKEpsilon = std::integral_constant;
+ return isDirichletCell_(element, fvGeometry, scv, pvIdx, IsKOmegaKEpsilon{});
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for fixed values at cell centers
+ *
+ * \param element The finite element
+ * \param scv the sub control volume
+ * \note used for cell-centered discretization schemes
+ */
+ template = 0>
+ PrimaryVariables dirichlet(const Element& element, const SubControlVolume& scv) const
+ {
+ const auto globalPos = scv.center();
+ PrimaryVariables values(initialAtPos(globalPos));
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for fixed values at cell centers
+ *
+ * \param element The finite element
+ * \param scv the sub control volume
+ * \note used for cell-centered discretization schemes
+ */
+ template = 0>
+ PrimaryVariables dirichlet(const Element& element, const SubControlVolume& scv) const
+ {
+ using SetDirichletCellForBothTurbEq = std::integral_constant;
+
+ return dirichletTurbulentTwoEq_(element, scv, SetDirichletCellForBothTurbEq{});
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet values at the boundary.
+ *
+ * \param element The finite element
+ * \param scvf the sub control volume face
+ */
+ PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
+ {
+ const auto globalPos = scvf.ipGlobal();
+ PrimaryVariables values(initialAtPos(globalPos));
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFaceVars The element face variables
+ * \param scvf The boundary sub control volume face
+ */
+ template
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ PrimaryVariables values(0.0);
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_);
+ values[Indices::conti0EqIdx] = massFlux[0];
+ values[Indices::conti0EqIdx + 1] = massFlux[1];
+
+ values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_);
+
+ }
+ return values;
+ }
+
+
+ // \}
+
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+
+ values[Indices::velocityXIdx] = inletVelocity();
+
+ if (isOnWallAtPos(globalPos))
+ values[Indices::velocityXIdx] = 0.0;
+
+ FluidState fluidState;
+ fluidState.setPressure(0, pressure_);
+ fluidState.setTemperature(temperature());
+ fluidState.setMoleFraction(0, 0, 1);
+ Scalar density = FluidSystem::density(fluidState, 0);
+ values[Indices::pressureIdx] = pressure_ - density * this->gravity()[dimWorld-1]* (this->fvGridGeometry().bBoxMax()[1] - globalPos[1]);
+
+ values[transportCompIdx] = inletMassFrac();
+ values[Indices::temperatureIdx] = temperature();
+
+ // turbulence model-specific initial conditions
+ static constexpr auto numEq = numTurbulenceEq(ModelTraits::turbulenceModel());
+ setInitialAtPos_(values, globalPos, Dune::index_constant{});
+
+ return values;
+ }
+
+
+ /*!
+ * \brief Returns the intrinsic permeability of required as input parameter
+ for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(element, scvf);
+ }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the
+ Beavers-Joseph-Saffman boundary condition.
+ */
+ Scalar alphaBJ(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
+ }
+
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ {
+ timeLoop_ = timeLoop;
+ }
+
+ Scalar time() const
+ {
+ return timeLoop_->time();
+ }
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ // \}
+
+private:
+ bool isInlet_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool isOutlet_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ bool isOnCouplingWall_(const SubControlVolumeFace& scvf) const
+ { return couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf);}
+
+
+ /*!
+ * \name Turbulent Boundary and Initial Conditions
+ */
+ // \{
+
+ //! Initial conditions for the zero-eq turbulence model (none)
+ void setInitialAtPos_(PrimaryVariables& values, const GlobalPosition &globalPos, Dune::index_constant<0>) const {}
+
+ //! Initial conditions for the one-eq turbulence model
+ void setInitialAtPos_(PrimaryVariables& values, const GlobalPosition &globalPos, Dune::index_constant<1>) const
+ {
+ values[Indices::viscosityTildeIdx] = viscosityTilde_;
+ if (isOnWallAtPos(globalPos))
+ values[Indices::viscosityTildeIdx] = 0.0;
+ }
+
+ //! Initial conditions for the komega, kepsilon and lowrekepsilon turbulence models
+ void setInitialAtPos_(PrimaryVariables& values, const GlobalPosition &globalPos, Dune::index_constant<2>) const
+ {
+ values[Indices::turbulentKineticEnergyIdx] = turbulentKineticEnergy_;
+ values[Indices::dissipationIdx] = dissipation_;
+ if (isOnWallAtPos(globalPos))
+ {
+ values[Indices::turbulentKineticEnergyIdx] = 0.0;
+ values[Indices::dissipationIdx] = 0.0;
+ }
+ }
+
+ //! Boundary condition types for the zero-eq turbulence model (none)
+ void setBcTypes_(BoundaryTypes& values, const GlobalPosition& pos, Dune::index_constant<0>) const {}
+
+ //! Boundary condition types for the one-eq turbulence model
+ void setBcTypes_(BoundaryTypes& values, const GlobalPosition& pos, Dune::index_constant<1>) const
+ {
+ if(isOutlet_(pos))
+ values.setOutflow(Indices::viscosityTildeIdx);
+ else // walls and inflow
+ values.setDirichlet(Indices::viscosityTildeIdx);
+ }
+
+ //! Boundary condition types for the komega, kepsilon and lowrekepsilon turbulence models
+ void setBcTypes_(BoundaryTypes& values,const GlobalPosition& pos, Dune::index_constant<2>) const
+ {
+ if(isOutlet_(pos))
+ {
+ values.setOutflow(Indices::turbulentKineticEnergyEqIdx);
+ values.setOutflow(Indices::dissipationEqIdx);
+ }
+ else
+ {
+ // walls and inflow
+ values.setDirichlet(Indices::turbulentKineticEnergyIdx);
+ values.setDirichlet(Indices::dissipationIdx);
+ }
+ }
+
+ //! Forward to ParentType
+ template
+ bool isDirichletCell_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolume& scv,
+ int pvIdx,
+ std::false_type) const
+ {
+ return ParentType::isDirichletCell(element, fvGeometry, scv, pvIdx);
+ }
+
+ //! Specialization for the KOmega and KEpsilon Models
+ template
+ bool isDirichletCell_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolume& scv,
+ int pvIdx,
+ std::true_type) const
+ {
+ using SetDirichletCellForBothTurbEq = std::integral_constant;
+ return isDirichletCellTurbulentTwoEq_(element, fvGeometry, scv, pvIdx, SetDirichletCellForBothTurbEq{});
+ }
+
+ //! Specialization for the KEpsilon Model
+ template
+ bool isDirichletCellTurbulentTwoEq_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolume& scv,
+ int pvIdx,
+ std::true_type) const
+ {
+ const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
+
+ // set a fixed turbulent kinetic energy and dissipation near the wall
+ if (this->inNearWallRegion(eIdx))
+ return pvIdx == Indices::turbulentKineticEnergyEqIdx || pvIdx == Indices::dissipationEqIdx;
+
+ // set a fixed dissipation at the matching point
+ if (this->isMatchingPoint(eIdx))
+ return pvIdx == Indices::dissipationEqIdx;// set a fixed dissipation (omega) for all cells at the wall
+
+ return false;
+ }
+
+ //! Specialization for the KOmega Model
+ template
+ bool isDirichletCellTurbulentTwoEq_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolume& scv,
+ int pvIdx,
+ std::false_type) const
+ {
+ // set a fixed dissipation (omega) for all cells at the wall
+ for (const auto& scvf : scvfs(fvGeometry))
+ if (isOnWallAtPos(scvf.center()) && pvIdx == Indices::dissipationIdx)
+ return true;
+
+ return false;
+
+ }
+
+ //! Specialization for the kepsilon
+ template
+ PrimaryVariables dirichletTurbulentTwoEq_(const Element& element,
+ const SubControlVolume& scv,
+ std::true_type) const
+ {
+ const auto globalPos = scv.center();
+ PrimaryVariables values(initialAtPos(globalPos));
+ unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element);
+
+ // fixed value for the turbulent kinetic energy
+ values[Indices::turbulentKineticEnergyEqIdx] = this->turbulentKineticEnergyWallFunction(elementIdx);
+
+ // fixed value for the dissipation
+ values[Indices::dissipationEqIdx] = this->dissipationWallFunction(elementIdx);
+
+ return values;
+ }
+
+ //! Specialization for the KOmega
+ template
+ PrimaryVariables dirichletTurbulentTwoEq_(const Element& element,
+ const SubControlVolume& scv,
+ std::false_type) const
+ {
+ const auto globalPos = scv.center();
+ PrimaryVariables values(initialAtPos(globalPos));
+ unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element);
+
+ const auto wallDistance = ParentType::wallDistance_[elementIdx];
+ using std::pow;
+ values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementIdx]
+ / (ParentType::betaOmega() * pow(wallDistance, 2));
+ return values;
+ }
+
+
+ // \}
+
+ Scalar inletVelocity_;
+ Scalar pressure_;
+ Scalar inletMassFrac_;
+ Scalar temperature_;
+
+ Scalar viscosityTilde_;
+ Scalar turbulentKineticEnergy_;
+ Scalar dissipation_;
+
+ DiffusionCoefficientAveragingType diffCoeffAvgType_;
+ Scalar eps_;
+ std::string problemName_;
+ std::string turbulenceModelName_;
+ TimeLoopPtr timeLoop_;
+ std::shared_ptr couplingManager_;
+};
+} // end namespace Dumux
+
+#endif // DUMUX_RANS_SUBPROBLEM_HH
diff --git a/test/multidomain/boundary/ransdarcy/1p2c_2p2c/scripts/params_evaporation.input b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/scripts/params_evaporation.input
new file mode 100644
index 0000000000000000000000000000000000000000..bac077ef932717f6f7a30589c90a352a0dd1cf45
--- /dev/null
+++ b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/scripts/params_evaporation.input
@@ -0,0 +1,77 @@
+[TimeLoop]
+DtInitial = 1e-4 # [s]
+MaxTimeStepSize = 350 # [s]
+TEnd = 691200 # [s] 8 days
+
+[RANS.Grid]
+Verbosity = true
+Positions0 = -1.0 0.0 1.0 2.0
+Positions1 = 0.5 1.0
+Grading0 = -1.1 1.0 1.1
+Grading1 = 1.3
+Cells0 = 15 30 15
+Cells1 = 20
+
+[Darcy.Grid]
+Verbosity = true
+Positions0 = 0.0 1.0
+Positions1 = 0.0 0.5
+Cells0 = 30
+Cells1 = 20
+Grading0 = 1.0
+Grading1 = -1.3
+
+[RANS.Problem]
+Name = rans
+InletVelocity = 2.5 # [m/s]
+InletPressure = 1e5 # [Pa]
+InletMassFrac = 0.008 # [-]
+Temperature = 298.15 # [K]
+
+[Darcy.Problem]
+Name = darcy
+Pressure = 1e5 # [Pa]
+Temperature = 298.15 # [K]
+InitPhasePresence = 3 # bothPhases
+InitialSaturation = 0.90 # [-]
+
+[Darcy.SpatialParams]
+Porosity = 0.41
+Permeability = 2.65e-10
+AlphaBeaversJoseph = 1.0
+Swr = 0.005
+Snr = 0.01
+VgAlpha = 6.37e-4
+VgN = 8.0
+
+[Problem]
+EnableGravity = true
+InterfaceDiffusionCoefficientAvg = FreeFlowOnly # Harmonic
+
+[RANS.RANS]
+UseStoredEddyViscosity = false
+EddyViscosityModel = "baldwinLomax"
+TurbulentSchmidtNumber = 0.7
+TurbulentPrandtlNumber = 0.85
+
+[RANS.KEpsilon]
+YPlusThreshold = 50. # should be small (10-30) for coarse grids large (30-60) for fine grids
+
+[Vtk]
+AddVelocity = true
+WriteFaceData = false
+OutputName = flat_interface
+
+[Component]
+SolidDensity = 2700
+SolidThermalConductivity = 2.8
+SolidHeatCapacity = 790
+
+[Newton]
+TargetSteps = 6
+MaxSteps = 8
+MaxRelativeShift = 1e-7
+
+[Assembly]
+NumericDifferenceMethod = 0
+NumericDifference.BaseEpsilon = 1e-8
diff --git a/test/multidomain/boundary/ransdarcy/1p2c_2p2c/scripts/plotevaporationrates.gp b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/scripts/plotevaporationrates.gp
new file mode 100644
index 0000000000000000000000000000000000000000..57e3925154bc7d5788305d2090b3b6b8f067a6d3
--- /dev/null
+++ b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/scripts/plotevaporationrates.gp
@@ -0,0 +1,16 @@
+reset
+set term pngcairo size 800,600 solid
+set output "./EvaporationRate.png"
+set datafile separator ';'
+DATA='./'
+set ylabel "Evaporation Rate [mm/d]"
+set xlabel "time [days]"
+set yrange [0:7]
+set xrange [0:8]
+set key right top samplen 1
+
+plot DATA.'storage_KEpsilon.csv' u ($1/86400):4 w l lc rgb "red" t 'KEpsilon',\
+ DATA.'storage_LowReKEpsilon.csv' u ($1/86400):4 w l lc rgb "orange" t 'Low Re KEpsilon',\
+ DATA.'storage_KOmega.csv' u ($1/86400):4 w l lc rgb "green" t 'KOmega',\
+ DATA.'storage_OneEq.csv' u ($1/86400):4 w l lc rgb "cyan" t 'OneEq',\
+ DATA.'storage_ZeroEq.csv' u ($1/86400):4 w l lc rgb "purple" t 'ZeroEq'
diff --git a/test/multidomain/boundary/ransdarcy/1p2c_2p2c/scripts/runevaporationturbulencemodels.sh b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/scripts/runevaporationturbulencemodels.sh
new file mode 100644
index 0000000000000000000000000000000000000000..0f996980a66281d8ec56386a3d90662fd8de73fe
--- /dev/null
+++ b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/scripts/runevaporationturbulencemodels.sh
@@ -0,0 +1,18 @@
+runSim () {
+./$1 $INPUT -Vtk.OutputName $2 -RANS.Problem.InletVelocity $3| tee -a logfile.out
+}
+
+### Evaporation Tests (Drying)
+INPUT=scripts/params_evaporation.input
+
+runSim test_md_boundary_darcy2p2cni_rans1p2cnikepsilon evaporation_kepsilon_3 3.0
+runSim test_md_boundary_darcy2p2cni_rans1p2cnikomega evaporation_komega_3 3.0
+runSim test_md_boundary_darcy2p2cni_rans1p2cnilowrekepsilon evaporation_lowrekepsilon_3 3.0
+runSim test_md_boundary_darcy2p2cni_rans1p2cnioneeq evaporation_oneeq_3 3.0
+runSim test_md_boundary_darcy2p2cni_rans1p2cnizeroeq evaporation_zeroeq_3 3.0
+runSim test_md_boundary_darcy2p2cni_rans1p2cnikomega evaporation_komega_0p3 0.3
+runSim test_md_boundary_darcy2p2cni_rans1p2cnikomega evaporation_komega_15 15.0
+runSim test_md_boundary_darcy2p2cni_rans1p2cnikomega evaporation_komega_30 30.0
+
+
+gnuplot scripts/plotevaporationrates.gp
diff --git a/test/multidomain/boundary/ransdarcy/1p2c_2p2c/spatialparams.hh b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/spatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a841a2e6017ebd6db3466849899b4bd2821aa955
--- /dev/null
+++ b/test/multidomain/boundary/ransdarcy/1p2c_2p2c/spatialparams.hh
@@ -0,0 +1,137 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief The spatial parameters class for the test problem using the 2pnc cc model
+ */
+#ifndef DUMUX_CONSERVATION_SPATIAL_PARAMS_HH
+#define DUMUX_CONSERVATION_SPATIAL_PARAMS_HH
+
+#include
+#include
+#include
+#include
+
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPNCModel
+ * \ingroup ImplicitTestProblems
+ *
+ * \brief The spatial parameters class for the test problem using the
+ * 1p cc model
+ */
+template
+class TwoPTwoCSpatialParams
+: public FVSpatialParams>
+{
+ using GridView = typename FVGridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using ParentType = FVSpatialParams>;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using EffectiveLaw = RegularizedVanGenuchten;
+
+public:
+ using MaterialLaw = EffToAbsLaw;
+ using MaterialLawParams = typename MaterialLaw::Params;
+ using PermeabilityType = Scalar;
+
+ TwoPTwoCSpatialParams(std::shared_ptr fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ permeability_ = getParam("Darcy.SpatialParams.Permeability");
+ alphaBJ_ = getParam("Darcy.SpatialParams.AlphaBeaversJoseph");
+ porosity_ = getParam("Darcy.SpatialParams.Porosity");
+
+ // residual saturations
+ params_.setSwr(getParam("Darcy.SpatialParams.Swr"));
+ params_.setSnr(getParam("Darcy.SpatialParams.Snr"));
+ // parameters for the vanGenuchten law
+ params_.setVgAlpha(getParam("Darcy.SpatialParams.VgAlpha"));
+ params_.setVgn(getParam("Darcy.SpatialParams.VgN"));
+ Scalar threshold = 0.01 * (1.0 - params_.swr() - params_.snr());
+ params_.setPcLowSw(getParam("Darcy.SpatialParams.PcLowSw", threshold));
+ params_.setPcHighSw(getParam("Darcy.SpatialParams.PcHighSw", 1.0-threshold));
+ }
+
+ /*!
+ * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
+ *
+ * \param globalPos The global position
+ * \return the intrinsic permeability
+ */
+ PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
+ { return permeability_; }
+
+ /*! \brief Defines the porosity in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return porosity_; }
+
+ /*! \brief Defines the Beavers-Joseph coefficient in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
+ { return alphaBJ_; }
+
+ /*!
+ * \brief Returns the parameter object for the Brooks-Corey material law.
+ * In this test, we use element-wise distributed material parameters.
+ *
+ * \param element The current element
+ * \param scv The sub-control volume inside the element.
+ * \param elemSol The solution at the dofs connected to the element.
+ * \return the material parameters object
+ */
+ template
+ const MaterialLawParams& materialLawParams(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ { return params_; }
+
+ /*!
+ * \brief Function for defining which phase is to be considered as the wetting phase.
+ *
+ * \return the wetting phase index
+ * \param globalPos The global position
+ */
+ template
+ int wettingPhaseAtPos(const GlobalPosition& globalPos) const
+ { return FluidSystem::phase0Idx; }
+
+private:
+ Scalar permeability_;
+ Scalar alphaBJ_;
+ Scalar porosity_;
+ MaterialLawParams params_;
+ static constexpr Scalar eps_ = 1.0e-7;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/boundary/ransdarcy/1p_1p/CMakeLists.txt b/test/multidomain/boundary/ransdarcy/1p_1p/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..73a7c38eb506825d7c411a18eaba2754dee3fda6
--- /dev/null
+++ b/test/multidomain/boundary/ransdarcy/1p_1p/CMakeLists.txt
@@ -0,0 +1,85 @@
+add_input_file_links()
+
+dune_add_test(NAME test_md_boundary_darcy1p_rans1pkepsilon
+ SOURCES main.cc
+ COMPILE_DEFINITIONS RANSTYPETAG=RANSSubProblemKEpsilon
+ CMAKE_GUARD HAVE_UMFPACK
+ LABELS freeflow rans multidomain
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_rans1pkepsilon_darcy-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1pkepsilon_darcy-00031.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_rans1pkepsilon_rans-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1pkepsilon_rans-00031.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1pkepsilon params.input
+ -Vtk.OutputName test_md_boundary_darcy1p_rans1pkepsilon")
+
+dune_add_test(NAME test_md_boundary_darcy1p_rans1pkomega
+ SOURCES main.cc
+ COMPILE_DEFINITIONS RANSTYPETAG=RANSSubProblemKOmega
+ LABELS freeflow rans multidomain
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_rans1pkomega_darcy-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1pkomega_darcy-00033.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_rans1pkomega_rans-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1pkomega_rans-00033.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1pkomega params.input
+ -Vtk.OutputName test_md_boundary_darcy1p_rans1pkomega")
+
+dune_add_test(NAME test_md_boundary_darcy1p_rans1plowrekepsilon
+ SOURCES main.cc
+ COMPILE_DEFINITIONS RANSTYPETAG=RANSSubProblemLowReKEpsilon
+ LABELS freeflow rans multidomain
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_rans1plowrekepsilon_darcy-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1plowrekepsilon_darcy-00035.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_rans1plowrekepsilon_rans-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1plowrekepsilon_rans-00035.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1plowrekepsilon params.input
+ -Vtk.OutputName test_md_boundary_darcy1p_rans1plowrekepsilon")
+
+dune_add_test(NAME test_md_boundary_darcy1p_rans1poneeq
+ SOURCES main.cc
+ COMPILE_DEFINITIONS RANSTYPETAG=RANSSubProblemOneEq
+ LABELS freeflow rans multidomain
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_rans1poneeq_darcy-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1poneeq_darcy-00032.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_rans1poneeq_rans-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1poneeq_rans-00032.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1poneeq params.input
+ -Vtk.OutputName test_md_boundary_darcy1p_rans1poneeq")
+
+dune_add_test(NAME test_md_boundary_darcy1p_rans1pzeroeq
+ SOURCES main.cc
+ COMPILE_DEFINITIONS RANSTYPETAG=RANSSubProblemZeroEq
+ LABELS freeflow rans multidomain
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_rans1pzeroeq_darcy-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1pzeroeq_darcy-00029.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_rans1pzeroeq_rans-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1pzeroeq_rans-00029.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_rans1pzeroeq params.input
+ -Vtk.OutputName test_md_boundary_darcy1p_rans1pzeroeq")
+
+dune_add_test(NAME test_md_boundary_forchheimer1p_rans1pzeroeq
+ SOURCES main.cc
+ COMPILE_DEFINITIONS RANSTYPETAG=RANSSubProblemZeroEq FORCHHEIMER=1
+ LABELS freeflow rans multidomain
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_forchheimer1p_rans1pzeroeq_darcy-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_forchheimer1p_rans1pzeroeq_darcy-00029.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_forchheimer1p_rans1pzeroeq_rans-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_forchheimer1p_rans1pzeroeq_rans-00029.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_forchheimer1p_rans1pzeroeq params.input
+ -Vtk.OutputName test_md_boundary_forchheimer1p_rans1pzeroeq")
diff --git a/test/multidomain/boundary/ransdarcy/1p_1p/main.cc b/test/multidomain/boundary/ransdarcy/1p_1p/main.cc
new file mode 100644
index 0000000000000000000000000000000000000000..60d78c707b9d1c37ee0a80e45935ffa59e7e0b0d
--- /dev/null
+++ b/test/multidomain/boundary/ransdarcy/1p_1p/main.cc
@@ -0,0 +1,284 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A test problem for the isothermal coupled RANS/Darcy problem (1p/1p)
+ */
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+
+#include "problem_darcy.hh"
+#include "problem_rans.hh"
+
+namespace Dumux {
+namespace Properties {
+
+template
+struct CouplingManager