diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt
index 9bc324aad612457f366991aaa6da3c965ec4237b..45acb89ecd3ed8e3f2b58c3260ba52309dd9ee7f 100644
--- a/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt
+++ b/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt
@@ -1,3 +1,5 @@
+add_subdirectory(convergencetest)
+
add_input_file_links()
add_executable(test_md_boundary_darcy1p_stokes1p EXCLUDE_FROM_ALL main.cc)
diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..69e2b6d8ed08292a1949a37457adee6e5b3615c2
--- /dev/null
+++ b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/CMakeLists.txt
@@ -0,0 +1,10 @@
+add_input_file_links()
+dune_symlink_to_source_files(FILES "convergencetest.py")
+
+dune_add_test(NAME test_md_boundary_darcy1p_stokes1p_convtest
+ SOURCES main.cc
+ LABELS multidomain multidomain_boundary stokesdarcy
+ TIMEOUT 1000
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ./convergencetest.py
+ CMD_ARGS test_md_boundary_darcy1p_stokes1p_convtest params.input)
diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/convergencetest.py b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/convergencetest.py
new file mode 100755
index 0000000000000000000000000000000000000000..6548362670d941409bb996c1ff91cbface471381
--- /dev/null
+++ b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/convergencetest.py
@@ -0,0 +1,132 @@
+#!/usr/bin/env python
+
+from math import *
+import subprocess
+import sys
+
+if len(sys.argv) < 2:
+ sys.stderr.write('Please provide a single argument <testname> to the script\n')
+ sys.exit(1)
+
+testname = str(sys.argv[1])
+testargs = [str(i) for i in sys.argv][2:]
+
+# remove the old log files
+subprocess.call(['rm', testname + '_stokes.log'])
+print("Removed old log file ({})!".format(testname + '_stokes.log'))
+subprocess.call(['rm', testname + '_darcy.log'])
+print("Removed old log file ({})!".format(testname + '_darcy.log'))
+
+# do the runs with different refinement
+for i in [0, 1, 2]:
+ subprocess.call(['./' + testname] + testargs + ['-Grid.Refinement', str(i),
+ '-Vtk.OutputName', testname])
+
+def checkRatesStokes():
+ # check the rates and append them to the log file
+ logfile = open(testname + '_stokes.log', "r+")
+
+ errorP = []
+ errorVx = []
+ errorVy = []
+ for line in logfile:
+ line = line.strip("\n")
+ line = line.strip("\[ConvergenceTest\]")
+ line = line.split()
+ errorP.append(float(line[2]))
+ errorVx.append(float(line[5]))
+ errorVy.append(float(line[8]))
+
+ resultsP = []
+ resultsVx = []
+ resultsVy = []
+ logfile.truncate(0)
+ logfile.write("n\terrorP\t\trateP\t\terrorVx\t\trateVx\t\terrorVy\t\trateVy\n")
+ logfile.write("-"*50 + "\n")
+ for i in range(len(errorP)-1):
+ if isnan(errorP[i]) or isinf(errorP[i]):
+ continue
+ if not ((errorP[i] < 1e-12 or errorP[i+1] < 1e-12) and (errorVx[i] < 1e-12 or errorVx[i+1] < 1e-12) and (errorVy[i] < 1e-12 or errorVy[i+1] < 1e-12)):
+ rateP = (log(errorP[i])-log(errorP[i+1]))/log(2)
+ rateVx = (log(errorVx[i])-log(errorVx[i+1]))/log(2)
+ rateVy = (log(errorVy[i])-log(errorVy[i+1]))/log(2)
+ message = "{}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\n".format(i, errorP[i], rateP, errorVx[i], rateVx, errorVy[i], rateVy)
+ logfile.write(message)
+ resultsP.append(rateP)
+ resultsVx.append(rateVx)
+ resultsVy.append(rateVy)
+ else:
+ logfile.write("error: exact solution!?")
+ i = len(errorP)-1
+ message = "{}\t{:0.4e}\t\t{}\t{:0.4e}\t\t{}\t{:0.4e}\t\t{}\n".format(i, errorP[i], "", errorVx[i], "", errorVy[i], "")
+ logfile.write(message)
+
+ logfile.close()
+ print("\nComputed the following convergence rates for {}:\n".format(testname))
+
+ subprocess.call(['cat', testname + '_stokes.log'])
+
+ return {"p" : resultsP, "v_x" : resultsVx, "v_y" : resultsVy}
+
+def checkRatesDarcy():
+ # check the rates and append them to the log file
+ logfile = open(testname + '_darcy.log', "r+")
+
+ errorP = []
+ for line in logfile:
+ line = line.strip("\n")
+ line = line.strip("\[ConvergenceTest\]")
+ line = line.split()
+ errorP.append(float(line[2]))
+
+ resultsP = []
+ logfile.truncate(0)
+ logfile.write("n\terrorP\t\trateP\n")
+ logfile.write("-"*50 + "\n")
+ for i in range(len(errorP)-1):
+ if isnan(errorP[i]) or isinf(errorP[i]):
+ continue
+ if not ((errorP[i] < 1e-12 or errorP[i+1] < 1e-12)):
+ rateP = (log(errorP[i])-log(errorP[i+1]))/log(2)
+ message = "{}\t{:0.4e}\t{:0.4e}\n".format(i, errorP[i], rateP)
+ logfile.write(message)
+ resultsP.append(rateP)
+ else:
+ logfile.write("error: exact solution!?")
+ i = len(errorP)-1
+ message = "{}\t{:0.4e}\n".format(i, errorP[i], "")
+ logfile.write(message)
+
+ logfile.close()
+ print("\nComputed the following convergence rates for {}:\n".format(testname))
+
+ subprocess.call(['cat', testname + '_darcy.log'])
+
+ return {"p" : resultsP}
+
+def checkRatesStokesAndDarcy():
+ resultsStokes = checkRatesStokes()
+ resultsDarcy = checkRatesDarcy()
+
+ def mean(numbers):
+ return float(sum(numbers)) / len(numbers)
+
+ # check the rates, we expect rates around 2
+ if mean(resultsStokes["p"]) < 2.05 and mean(resultsStokes["p"]) < 1.84:
+ sys.stderr.write("*"*70 + "\n" + "The convergence rates for pressure were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
+ sys.exit(1)
+
+ if mean(resultsStokes["v_x"]) < 2.05 and mean(resultsStokes["v_x"]) < 1.95:
+ sys.stderr.write("*"*70 + "\n" + "The convergence rates for x-velocity were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
+ sys.exit(1)
+
+ if mean(resultsStokes["v_y"]) < 2.05 and mean(resultsStokes["v_y"]) < 1.95:
+ sys.stderr.write("*"*70 + "\n" + "The convergence rates for y-velocity were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
+ sys.exit(1)
+
+ if mean(resultsDarcy["p"]) < 2.05 and mean(resultsDarcy["p"]) < 1.95:
+ sys.stderr.write("*"*70 + "\n" + "The convergence rates for pressure were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
+ sys.exit(1)
+
+
+checkRatesStokesAndDarcy()
diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/main.cc b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/main.cc
new file mode 100644
index 0000000000000000000000000000000000000000..331fc2cdc2c728113d8c39cc3286c04f7206b596
--- /dev/null
+++ b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/main.cc
@@ -0,0 +1,401 @@
+// -*- 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 <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief A test problem for the coupled Stokes/Darcy problem (1p).
+ * The analytical solution is given in Shiue et al., 2018:
+ * "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
+ */
+
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/partial.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager_yasp.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+#include <test/freeflow/navierstokes/l2error.hh>
+
+#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
+
+#include "problem_darcy.hh"
+#include "problem_stokes.hh"
+
+namespace Dumux {
+namespace Properties {
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::StokesOneP>
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOneP>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyOneP>
+{
+ using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesOneP, Properties::TTag::StokesOneP, TypeTag>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+} // end namespace Properties
+} // end namespace Dumux
+
+/*!
+* \brief Creates analytical solution.
+* Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces
+* \param problem the problem for which to evaluate the analytical solution
+*/
+template<class Scalar, class Problem>
+auto createStokesAnalyticalSolution(const Problem& problem)
+{
+ const auto& gridGeometry = problem.gridGeometry();
+ using GridView = typename std::decay_t<decltype(gridGeometry)>::GridView;
+
+ static constexpr auto dim = GridView::dimension;
+ static constexpr auto dimWorld = GridView::dimensionworld;
+
+ using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
+
+ std::vector<Scalar> analyticalPressure;
+ std::vector<VelocityVector> analyticalVelocity;
+ std::vector<Scalar> analyticalVelocityOnFace;
+
+ analyticalPressure.resize(gridGeometry.numCellCenterDofs());
+ analyticalVelocity.resize(gridGeometry.numCellCenterDofs());
+ analyticalVelocityOnFace.resize(gridGeometry.numFaceDofs());
+
+ using Indices = typename Problem::Indices;
+ for (const auto& element : elements(gridGeometry.gridView()))
+ {
+ auto fvGeometry = localView(gridGeometry);
+ fvGeometry.bindElement(element);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ auto ccDofIdx = scv.dofIndex();
+ auto ccDofPosition = scv.dofPosition();
+ auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition);
+
+ // velocities on faces
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ const auto faceDofIdx = scvf.dofIndex();
+ const auto faceDofPosition = scvf.center();
+ const auto dirIdx = scvf.directionIndex();
+ const auto analyticalSolutionAtFace = problem.analyticalSolution(faceDofPosition);
+ analyticalVelocityOnFace[faceDofIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
+ }
+
+ analyticalPressure[ccDofIdx] = analyticalSolutionAtCc[Indices::pressureIdx];
+
+ for (int dirIdx = 0; dirIdx < dim; ++dirIdx)
+ analyticalVelocity[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)];
+ }
+ }
+
+ return std::make_tuple(analyticalPressure, analyticalVelocity, analyticalVelocityOnFace);
+}
+
+/*!
+* \brief Creates analytical solution.
+* Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces
+* \param problem the problem for which to evaluate the analytical solution
+*/
+template<class Scalar, class Problem>
+auto createDarcyAnalyticalSolution(const Problem& problem)
+{
+ const auto& gridGeometry = problem.gridGeometry();
+ using GridView = typename std::decay_t<decltype(gridGeometry)>::GridView;
+
+ static constexpr auto dim = GridView::dimension;
+ static constexpr auto dimWorld = GridView::dimensionworld;
+
+ using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
+
+ std::vector<Scalar> analyticalPressure;
+ std::vector<VelocityVector> analyticalVelocity;
+
+ analyticalPressure.resize(gridGeometry.numDofs());
+ analyticalVelocity.resize(gridGeometry.numDofs());
+
+ for (const auto& element : elements(gridGeometry.gridView()))
+ {
+ auto fvGeometry = localView(gridGeometry);
+ fvGeometry.bindElement(element);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto ccDofIdx = scv.dofIndex();
+ const auto ccDofPosition = scv.dofPosition();
+ const auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition);
+ analyticalPressure[ccDofIdx] = analyticalSolutionAtCc[dim];
+
+ for (int dirIdx = 0; dirIdx < dim; ++dirIdx)
+ analyticalVelocity[ccDofIdx][dirIdx] = analyticalSolutionAtCc[dirIdx];
+ }
+ }
+
+ return std::make_tuple(analyticalPressure, analyticalVelocity);
+}
+
+template<class Problem, class SolutionVector>
+void printStokesL2Error(const Problem& problem, const SolutionVector& x)
+{
+ using namespace Dumux;
+ using Scalar = double;
+ using TypeTag = Properties::TTag::StokesOneP;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+
+ using L2Error = NavierStokesTestL2Error<Scalar, ModelTraits, PrimaryVariables>;
+ const auto l2error = L2Error::calculateL2Error(problem, x);
+ const int numCellCenterDofs = problem.gridGeometry().numCellCenterDofs();
+ const int numFaceDofs = problem.gridGeometry().numFaceDofs();
+ std::ostream tmpOutputObject(std::cout.rdbuf()); // create temporary output with fixed formatting without affecting std::cout
+ tmpOutputObject << std::setprecision(8) << "** L2 error (abs/rel) for "
+ << std::setw(6) << numCellCenterDofs << " cc dofs and " << numFaceDofs << " face dofs (total: " << numCellCenterDofs + numFaceDofs << "): "
+ << std::scientific
+ << "L2(p) = " << l2error.first[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx]
+ << " , L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx]
+ << " , L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx]
+ << std::endl;
+
+ // write the norm into a log file
+ std::ofstream logFile;
+ logFile.open(problem.name() + ".log", std::ios::app);
+ logFile << "[ConvergenceTest] L2(p) = " << l2error.first[Indices::pressureIdx] << " L2(vx) = " << l2error.first[Indices::velocityXIdx] << " L2(vy) = " << l2error.first[Indices::velocityYIdx] << std::endl;
+ logFile.close();
+}
+
+template<class Problem, class SolutionVector>
+void printDarcyL2Error(const Problem& problem, const SolutionVector& x)
+{
+ using namespace Dumux;
+ using Scalar = double;
+
+ Scalar l2error = 0.0;
+
+ for (const auto& element : elements(problem.gridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(problem.gridGeometry());
+ fvGeometry.bindElement(element);
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto dofIdx = scv.dofIndex();
+ const Scalar delta = x[dofIdx] - problem.analyticalSolution(scv.center())[2/*pressureIdx*/];
+ l2error += scv.volume()*(delta*delta);
+ }
+ }
+ using std::sqrt;
+ l2error = sqrt(l2error);
+
+ const auto numDofs = problem.gridGeometry().numDofs();
+ std::ostream tmp(std::cout.rdbuf());
+ tmp << std::setprecision(8) << "** L2 error (abs) for "
+ << std::setw(6) << numDofs << " cc dofs "
+ << std::scientific
+ << "L2 error = " << l2error
+ << std::endl;
+
+ // write the norm into a log file
+ std::ofstream logFile;
+ logFile.open(problem.name() + ".log", std::ios::app);
+ logFile << "[ConvergenceTest] L2(p) = " << l2error << std::endl;
+ logFile.close();
+}
+
+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 StokesTypeTag = Properties::TTag::StokesOneP;
+ using DarcyTypeTag = Properties::TTag::DarcyOneP;
+
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
+ StokesGridManager stokesGridManager;
+ stokesGridManager.init("Stokes"); // pass parameter group
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridManager.grid().leafGridView();
+ const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using StokesGridGeometry = GetPropType<StokesTypeTag, Properties::GridGeometry>;
+ auto stokesGridGeometry = std::make_shared<StokesGridGeometry>(stokesGridView);
+ stokesGridGeometry->update();
+ using DarcyGridGeometry = GetPropType<DarcyTypeTag, Properties::GridGeometry>;
+ auto darcyGridGeometry = std::make_shared<DarcyGridGeometry>(darcyGridView);
+ darcyGridGeometry->update();
+
+ using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+
+ // the coupling manager
+ using CouplingManager = StokesDarcyCouplingManager<Traits>;
+ auto couplingManager = std::make_shared<CouplingManager>(stokesGridGeometry, darcyGridGeometry);
+
+ // the indices
+ constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+ constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
+ constexpr auto darcyIdx = CouplingManager::darcyIdx;
+
+ // the problem (initial and boundary conditions)
+ using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesGridGeometry, couplingManager);
+ using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyGridGeometry, couplingManager);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[stokesCellCenterIdx].resize(stokesGridGeometry->numCellCenterDofs());
+ sol[stokesFaceIdx].resize(stokesGridGeometry->numFaceDofs());
+ sol[darcyIdx].resize(darcyGridGeometry->numDofs());
+
+ // get a solution vector storing references to the two Stokes solution vectors
+ auto stokesSol = partial(sol, stokesFaceIdx, stokesCellCenterIdx);
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesGridGeometry);
+ stokesGridVariables->init(stokesSol);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx]);
+
+ // intialize the vtk output module
+ using Scalar = typename Traits::Scalar;
+ StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name());
+ GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
+ const auto stokesAnalyticalSolution = createStokesAnalyticalSolution<Scalar>(*stokesProblem);
+ stokesVtkWriter.addField(std::get<0>(stokesAnalyticalSolution), "pressureExact");
+ stokesVtkWriter.addField(std::get<1>(stokesAnalyticalSolution), "velocityExact");
+ stokesVtkWriter.addFaceField(std::get<2>(stokesAnalyticalSolution), "faceVelocityExact");
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyProblem->name());
+ using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
+ darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
+ GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
+ const auto darcyAnalyticalSolution = createDarcyAnalyticalSolution<Scalar>(*darcyProblem);
+ darcyVtkWriter.addField(std::get<0>(darcyAnalyticalSolution), "pressureExact");
+ darcyVtkWriter.addField(std::get<1>(darcyAnalyticalSolution), "velocityExact");
+ darcyVtkWriter.write(0.0);
+
+ // the assembler for a stationary problem
+ using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
+ std::make_tuple(stokesGridGeometry->faceFVGridGeometryPtr(),
+ stokesGridGeometry->cellCenterFVGridGeometryPtr(),
+ darcyGridGeometry),
+ std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
+ stokesGridVariables->cellCenterGridVariablesPtr(),
+ darcyGridVariables),
+ couplingManager);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ // solve the non-linear system
+ nonLinearSolver.solve(sol);
+
+ // write vtk output
+ stokesVtkWriter.write(1.0);
+ darcyVtkWriter.write(1.0);
+
+ printStokesL2Error(*stokesProblem, stokesSol);
+ printDarcyL2Error(*darcyProblem, sol[darcyIdx]);
+
+ ////////////////////////////////////////////////////////////
+ // 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/stokesdarcy/1p_1p/convergencetest/manufactured_solution.py b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/manufactured_solution.py
new file mode 100644
index 0000000000000000000000000000000000000000..e06ece60e357bfee1da173bd577eb9da113dd39e
--- /dev/null
+++ b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/manufactured_solution.py
@@ -0,0 +1,81 @@
+#!/usr/bin/env python
+
+"""
+This script determines the source terms needed for analytical solutions for coupled Stokes/Darcy problems.
+Given an analytical solution, it evaluates the momentum and mass balance equations and outputs source terms
+such that the balance equations are fulfilled. It furthermore calculates the Darcy velocity from a given pressure solution.
+"""
+
+import sympy as sp
+import numpy as np
+from sympy import *
+
+# case = "Rybak"
+# case = "Shiue_1"
+case = "Shiue_2"
+
+# divergence of a matrix
+def div(M):
+ return np.array([sp.diff(M[0,0],x) + sp.diff(M[0,1],y) , sp.diff(M[1,0],x) + sp.diff(M[1,1],y) ])
+
+# gradient of a vector
+def grad(v):
+ return np.array([[sp.diff(v[0],x), sp.diff(v[0],y)], [sp.diff(v[1],x), sp.diff(v[1],y)]])
+
+# coordinates
+x, y = sp.symbols("x y")
+
+# analytical solutions for v and p in free-flow domain (see reference given in problems)
+def analyticalSolutionStokes(case):
+ if case == "Rybak":
+ vFF = np.array([-sp.cos(sp.pi*x)*sp.sin(sp.pi*y), sp.sin(sp.pi*x)*sp.cos(sp.pi*y)])
+ pFF = 0.5*y*sp.sin(sp.pi*x)
+ elif case == "Shiue_1":
+ vFF = np.array([-1.0/sp.pi * sp.exp(y) * sp.sin(sp.pi*x), (sp.exp(y) - sp.exp(1)) * sp.cos(sp.pi*x)])
+ pFF = 2.0*sp.exp(y) * sp.cos(sp.pi*x)
+ else: # Shiue_2
+ vFF = np.array([(y-1.0)*(y-1.0) + x*(y-1.0) + 3.0*x - 1.0, x*(x-1.0) - 0.5*(y-1.0)*(y-1.0) - 3.0*y + 1.0])
+ pFF = 2.0*x + y - 1.0
+ return [vFF, pFF]
+
+vFF, pFF = analyticalSolutionStokes(case)
+
+# individual terms of the Navier-Stokes eq.
+vvT = np.outer(vFF, vFF)
+gradV = grad(vFF)
+gradVT = grad(vFF).T
+pI = np.array([[pFF, 0], [0, pFF]])
+
+# complete momentum flux and its divergence
+#momentumFlux = vvT - (gradV + gradVT) +pI
+momentumFlux = - (gradV + gradVT) +pI # only Stokes
+divMomentumFlux = div(momentumFlux)
+
+print("Source terms for case", case)
+
+# solution for source term
+print(" \nStokes:")
+print("Source term mass:", sp.diff(vFF[0],x) + sp.diff(vFF[1], y))
+print("Source term momentum x:", divMomentumFlux[0])
+print("Source term momentum y:", divMomentumFlux[1])
+
+# analytical solutions for p in Darcy domain (see reference given in problems)
+def analyticalSolutionDarcy(case):
+ if case == "Rybak":
+ pD = 0.5*y*y*sp.sin(sp.pi*x)
+ elif case == "Shiue_1":
+ pD = (sp.exp(y) - y*sp.exp(1)) * sp.cos(sp.pi*x)
+ else: # Shiue_2
+ pD = x*(1.0-x)*(y-1.0) + pow(y-1.0, 3)/3.0 + 2.0*x + 2.0*y + 4.0
+ return pD
+
+pD = analyticalSolutionDarcy(case)
+
+gradPdK = np.array([sp.diff(pD,x), sp.diff(pD,y)])
+vD = -gradPdK
+divVd = sp.diff(vD[0],x) + sp.diff(vD[1],y)
+
+print("\nDarcy:")
+print("Source term mass:", simplify(divVd))
+print("v x:", simplify(vD[0]))
+print("v_y", simplify(vD[1]))
diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/params.input b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/params.input
new file mode 100644
index 0000000000000000000000000000000000000000..01107ccb8560e0ee992c7f0164e0a303f81a32fb
--- /dev/null
+++ b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/params.input
@@ -0,0 +1,35 @@
+[Darcy.Grid]
+UpperRight = 1 1
+Cells = 40 40
+
+[Stokes.Grid]
+LowerLeft = 0 1
+UpperRight = 1 2
+Cells = 40 40
+
+[Stokes.Problem]
+Name = stokes
+EnableInertiaTerms = false
+
+[Darcy.Problem]
+Name = darcy
+
+[Darcy.SpatialParams]
+Permeability = 1.0
+AlphaBeaversJoseph = 1.0
+
+[Vtk]
+OutputName = test_md_boundary_stokes1p_darcy1p_convergencetest
+
+[Problem]
+EnableGravity = false
+
+[Vtk]
+AddVelocity = 1
+
+[Component]
+LiquidDensity = 1.0
+LiquidKinematicViscosity = 1.0
+
+[Assembly]
+NumericDifference.BaseEpsilon = 1e-4
diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/problem_darcy.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/problem_darcy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..2b92475b442861bb05e4420428d38801c8a89655
--- /dev/null
+++ b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/problem_darcy.hh
@@ -0,0 +1,360 @@
+// -*- 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 <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief The Darcy sub-problem of coupled Stokes-Darcy convergence test
+ */
+
+#ifndef DUMUX_DARCY_SUBPROBLEM_HH
+#define DUMUX_DARCY_SUBPROBLEM_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/discretization/cctpfa.hh>
+
+#include <dumux/porousmediumflow/1p/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+
+#include "../spatialparams.hh"
+
+#include <dumux/material/components/constant.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+
+namespace Dumux {
+template <class TypeTag>
+class DarcySubProblem;
+
+namespace Properties {
+// Create new type tags
+namespace TTag {
+struct DarcyOneP { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
+} // end namespace TTag
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyOneP> { using type = Dumux::DarcySubProblem<TypeTag>; };
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyOneP>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::Constant<1, Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyOneP> { using type = Dune::YaspGrid<2>; };
+
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyOneP>
+{
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = OnePSpatialParams<GridGeometry, Scalar>;
+};
+} // end namespace Properties
+
+/*!
+ * \ingroup BoundaryTests
+ * \brief The Darcy sub-problem of coupled Stokes-Darcy convergence test
+ */
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+
+ static constexpr auto velocityXIdx = 0;
+ static constexpr auto velocityYIdx = 1;
+ static constexpr auto pressureIdx = 2;
+
+ enum class TestCase
+ {
+ ShiueExampleOne, ShiueExampleTwo, Rybak
+ };
+
+public:
+ //! export the Indices
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
+ DarcySubProblem(std::shared_ptr<const GridGeometry> gridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(gridGeometry, "Darcy"), couplingManager_(couplingManager)
+ {
+ problemName_ = getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
+
+ const auto testCaseInput = getParamFromGroup<std::string>(this->paramGroup(), "Problem.TestCase", "ShiueExampleTwo");
+ if (testCaseInput == "ShiueExampleOne")
+ testCase_ = TestCase::ShiueExampleOne;
+ else if (testCaseInput == "ShiueExampleTwo")
+ testCase_ = TestCase::ShiueExampleTwo;
+ else if (testCaseInput == "Rybak")
+ testCase_ = TestCase::Rybak;
+ else
+ DUNE_THROW(Dune::InvalidStateException, testCaseInput + " is not a valid test case");
+ }
+
+ /*!
+ * \brief The problem name.
+ */
+ const std::string& name() const
+ {
+ return problemName_;
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Returns the temperature within the domain in [K].
+ *
+ */
+ Scalar temperature() const
+ { return 273.15 + 10; } // 10°C
+ // \}
+
+ /*!
+ * \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;
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values.setAllCouplingNeumann();
+ else
+ values.setAllDirichlet();
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element for which the Dirichlet boundary condition is set
+ * \param scvf The boundary subcontrolvolumeface
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ const auto p = analyticalSolution(scvf.center())[pressureIdx];
+ return PrimaryVariables(p);
+ }
+
+ /*!
+ * \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 elemFluxVarsCache Flux variables caches for all faces in stencil
+ * \param scvf The boundary sub control volume face
+ *
+ * For this method, the \a values variable stores primary variables.
+ */
+ template<class ElementVolumeVariables, class ElementFluxVarsCache>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVarsCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ /*!
+ * \brief Evaluates the source term for all phases within a given
+ * sub control volume.
+ *
+ * \param element The element for which the source term is set
+ * \param fvGeometry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scv The sub control volume
+ */
+ NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
+ {
+ switch (testCase_)
+ {
+ case TestCase::ShiueExampleOne:
+ return rhsShiueEtAlExampleOne_(globalPos);
+ case TestCase::ShiueExampleTwo:
+ return rhsShiueEtAlExampleTwo_(globalPos);
+ case TestCase::Rybak:
+ return rhsRybak_(globalPos);
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
+ }
+ }
+
+ // \}
+
+ /*!
+ * \brief Evaluates the initial value for a control volume.
+ *
+ * \param element The element
+ *
+ * For this method, the \a priVars parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initial(const Element &element) const
+ {
+ return PrimaryVariables(0.0);
+ }
+
+ /*!
+ * \brief Returns the analytical solution of the problem at a given position.
+ *
+ * \param globalPos The global position
+ * \param time The current simulation time
+ */
+ auto analyticalSolution(const GlobalPosition& globalPos) const
+ {
+ switch (testCase_)
+ {
+ case TestCase::ShiueExampleOne:
+ return analyticalSolutionShiueEtAlExampleOne_(globalPos);
+ case TestCase::ShiueExampleTwo:
+ return analyticalSolutionShiueEtAlExampleTwo_(globalPos);
+ case TestCase::Rybak:
+ return analyticalSolutionRybak_(globalPos);
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
+ }
+ }
+
+ // \}
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+private:
+
+ // see Rybak et al., 2015: "Multirate time integration for coupled saturated/unsaturated porous medium and free flow systems"
+ Dune::FieldVector<Scalar, 3> analyticalSolutionRybak_(const GlobalPosition& globalPos) const
+ {
+ Dune::FieldVector<Scalar, 3> sol(0.0);
+ const Scalar x = globalPos[0];
+ const Scalar y = globalPos[1];
+
+ using std::exp; using std::sin; using std::cos;
+ sol[velocityXIdx] = -0.5*M_PI*y*y*cos(M_PI*x);
+ sol[velocityYIdx] = -1.0*y*sin(M_PI*x);
+ sol[pressureIdx] = 0.5*y*y*sin(M_PI*x);
+ return sol;
+ }
+
+ // see Rybak et al., 2015: "Multirate time integration for coupled saturated/unsaturated porous medium and free flow systems"
+ NumEqVector rhsRybak_(const GlobalPosition& globalPos) const
+ {
+ const Scalar x = globalPos[0];
+ const Scalar y = globalPos[1];
+ using std::sin;
+ return NumEqVector((0.5*M_PI*y*M_PI*y - 1)*sin(M_PI*x));
+ }
+
+ // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
+ Dune::FieldVector<Scalar, 3> analyticalSolutionShiueEtAlExampleOne_(const GlobalPosition& globalPos) const
+ {
+ Dune::FieldVector<Scalar, 3> sol(0.0);
+ const Scalar x = globalPos[0];
+ const Scalar y = globalPos[1];
+
+ using std::exp; using std::sin; using std::cos;
+ sol[pressureIdx] = (exp(y) - y*exp(1)) * cos(M_PI*x);
+ sol[velocityXIdx] = M_PI*(-exp(1)*y + exp(y))*sin(M_PI*x);
+ sol[velocityYIdx] = (exp(1) - exp(y))*cos(M_PI*x);
+
+ return sol;
+ }
+
+ // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
+ NumEqVector rhsShiueEtAlExampleOne_(const GlobalPosition& globalPos) const
+ {
+ const Scalar x = globalPos[0];
+ const Scalar y = globalPos[1];
+ using std::exp; using std::sin; using std::cos;
+ return NumEqVector(cos(M_PI*x) * (M_PI*M_PI*(exp(y) - y*exp(1)) - exp(y)));
+ }
+
+ // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
+ Dune::FieldVector<Scalar, 3> analyticalSolutionShiueEtAlExampleTwo_(const GlobalPosition& globalPos) const
+ {
+ Dune::FieldVector<Scalar, 3> sol(0.0);
+ const Scalar x = globalPos[0];
+ const Scalar y = globalPos[1];
+
+ sol[pressureIdx] = x*(1.0-x)*(y-1.0) + (y-1.0)*(y-1.0)*(y-1.0)/3.0 + 2.0*x + 2.0*y + 4.0;
+ sol[velocityXIdx] = x*(y - 1.0) + (x - 1.0)*(y - 1.0) - 2.0;
+ sol[velocityYIdx] = x*(x - 1.0) - 1.0*(y - 1.0)*(y - 1.0) - 2.0;
+
+ return sol;
+ }
+
+ // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
+ NumEqVector rhsShiueEtAlExampleTwo_(const GlobalPosition& globalPos) const
+ { return NumEqVector(0.0); }
+
+ static constexpr Scalar eps_ = 1e-7;
+ std::shared_ptr<CouplingManager> couplingManager_;
+ std::string problemName_;
+ TestCase testCase_;
+};
+} // end namespace Dumux
+
+#endif //DUMUX_DARCY_SUBPROBLEM_HH
diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/problem_stokes.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/problem_stokes.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3a02f723dd0ba43181d822eda83f124e36c108d4
--- /dev/null
+++ b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/problem_stokes.hh
@@ -0,0 +1,366 @@
+// -*- 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 <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoundaryTests
+ * \brief The Stokes sub-problem of coupled Stokes-Darcy convergence test
+ */
+
+#ifndef DUMUX_STOKES_SUBPROBLEM_HH
+#define DUMUX_STOKES_SUBPROBLEM_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/material/components/constant.hh>
+
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/freeflow/navierstokes/model.hh>
+
+namespace Dumux {
+template <class TypeTag>
+class StokesSubProblem;
+
+namespace Properties {
+// Create new type tags
+namespace TTag {
+struct StokesOneP { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
+} // end namespace TTag
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::StokesOneP>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::Constant<1, Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::StokesOneP> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::StokesOneP> { using type = Dumux::StokesSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableGridGeometryCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+} // end namespace Properties
+
+/*!
+ * \ingroup BoundaryTests
+ * \brief The Stokes sub-problem of coupled Stokes-Darcy convergence test
+ */
+template <class TypeTag>
+class StokesSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+
+ enum class TestCase
+ {
+ ShiueExampleOne, ShiueExampleTwo, Rybak
+ };
+
+public:
+ //! export the Indices
+ using Indices = typename ModelTraits::Indices;
+
+ StokesSubProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(gridGeometry, "Stokes"), couplingManager_(couplingManager)
+ {
+ problemName_ = getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
+
+ const auto testCaseInput = getParamFromGroup<std::string>(this->paramGroup(), "Problem.TestCase", "ShiueExampleTwo");
+ if (testCaseInput == "ShiueExampleOne")
+ testCase_ = TestCase::ShiueExampleOne;
+ else if (testCaseInput == "ShiueExampleTwo")
+ testCase_ = TestCase::ShiueExampleTwo;
+ else if (testCaseInput == "Rybak")
+ testCase_ = TestCase::Rybak;
+ else
+ DUNE_THROW(Dune::InvalidStateException, testCaseInput + " is not a valid test case");
+ }
+
+ /*!
+ * \brief The problem name.
+ */
+ const std::string& name() const
+ {
+ return problemName_;
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Returns the temperature within the domain in [K].
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature() const
+ { return 273.15 + 10; } // 10°C
+
+ /*!
+ * \brief Returns the sources within the domain.
+ *
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+ {
+ switch (testCase_)
+ {
+ case TestCase::ShiueExampleOne:
+ return rhsShiueEtAlExampleOne_(globalPos);
+ case TestCase::ShiueExampleTwo:
+ return rhsShiueEtAlExampleTwo_(globalPos);
+ case TestCase::Rybak:
+ return rhsRybak_(globalPos);
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
+ }
+ }
+
+ // \}
+
+ /*!
+ * \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 values;
+
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+
+ if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setBeaversJoseph(Indices::momentumXBalanceIdx);
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Dirichlet control volume.
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+ {
+ return analyticalSolution(globalPos);
+ }
+
+ /*!
+ * \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 ElementFaceVariables>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+ }
+ return values;
+ }
+
+ // \}
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluates the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
+ {
+ return PrimaryVariables(0.0);
+ }
+
+ /*!
+ * \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());
+ }
+
+ /*!
+ * \brief Returns the analytical solution of the problem at a given position.
+ *
+ * \param globalPos The global position
+ * \param time The current simulation time
+ */
+ PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
+ {
+ switch (testCase_)
+ {
+ case TestCase::ShiueExampleOne:
+ return analyticalSolutionShiueEtAlExampleOne_(globalPos);
+ case TestCase::ShiueExampleTwo:
+ return analyticalSolutionShiueEtAlExampleTwo_(globalPos);
+ case TestCase::Rybak:
+ return analyticalSolutionRybak_(globalPos);
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
+ }
+ }
+
+ // \}
+
+private:
+
+ // see Rybak et al., 2015: "Multirate time integration for coupled saturated/unsaturated porous medium and free flow systems"
+ PrimaryVariables analyticalSolutionRybak_(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables sol(0.0);
+ const Scalar x = globalPos[0];
+ const Scalar y = globalPos[1];
+
+ using std::sin; using std::cos;
+ sol[Indices::velocityXIdx] = -cos(M_PI*x)*sin(M_PI*y);
+ sol[Indices::velocityYIdx] = sin(M_PI*x)*cos(M_PI*y);
+ sol[Indices::pressureIdx] = 0.5*y*sin(M_PI*x);
+ return sol;
+ }
+
+ // see Rybak et al., 2015: "Multirate time integration for coupled saturated/unsaturated porous medium and free flow systems"
+ NumEqVector rhsRybak_(const GlobalPosition& globalPos) const
+ {
+ const Scalar x = globalPos[0];
+ const Scalar y = globalPos[1];
+ NumEqVector source(0.0);
+ using std::sin; using std::cos;
+ source[Indices::momentumXBalanceIdx] = 0.5*M_PI*y*cos(M_PI*x) - 2*M_PI*M_PI*cos(M_PI*x)*sin(M_PI*y);
+ source[Indices::momentumYBalanceIdx] = 0.5*sin(M_PI*x) + 2*M_PI*M_PI*sin(M_PI*x)*cos(M_PI*y);
+ return source;
+ }
+
+ // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
+ PrimaryVariables analyticalSolutionShiueEtAlExampleOne_(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables sol(0.0);
+ const Scalar x = globalPos[0];
+ const Scalar y = globalPos[1];
+
+ using std::exp; using std::sin; using std::cos;
+ sol[Indices::velocityXIdx] = -1/M_PI * exp(y) * sin(M_PI*x);
+ sol[Indices::velocityYIdx] = (exp(y) - exp(1)) * cos(M_PI*x);
+ sol[Indices::pressureIdx] = 2*exp(y) * cos(M_PI*x);
+ return sol;
+ }
+
+ // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
+ NumEqVector rhsShiueEtAlExampleOne_(const GlobalPosition& globalPos) const
+ {
+ const Scalar x = globalPos[0];
+ const Scalar y = globalPos[1];
+ using std::exp; using std::sin; using std::cos;
+ NumEqVector source(0.0);
+ source[Indices::momentumXBalanceIdx] = exp(y)*sin(M_PI*x) * (1/M_PI -3*M_PI);
+ source[Indices::momentumYBalanceIdx] = cos(M_PI*x) * (M_PI*M_PI*(exp(y)- exp(1)) + exp(y));
+ return source;
+ }
+
+ // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
+ PrimaryVariables analyticalSolutionShiueEtAlExampleTwo_(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables sol(0.0);
+ const Scalar x = globalPos[0];
+ const Scalar y = globalPos[1];
+
+ sol[Indices::velocityXIdx] = (y-1.0)*(y-1.0) + x*(y-1.0) + 3.0*x - 1.0;
+ sol[Indices::velocityYIdx] = x*(x-1.0) - 0.5*(y-1.0)*(y-1.0) - 3.0*y + 1.0;
+ sol[Indices::pressureIdx] = 2.0*x + y - 1.0;
+ return sol;
+ }
+
+ // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
+ NumEqVector rhsShiueEtAlExampleTwo_(const GlobalPosition& globalPos) const
+ { return NumEqVector(0.0); }
+
+ static constexpr Scalar eps_ = 1e-7;
+ std::string problemName_;
+ std::shared_ptr<CouplingManager> couplingManager_;
+ TestCase testCase_;
+};
+} // end namespace Dumux
+
+#endif // DUMUX_STOKES_SUBPROBLEM_HH