From fe48378e83e389edb8f4833c23015a34739aa996 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 15 Jun 2021 20:44:40 +0200
Subject: [PATCH] [md][boundary] Port Darcy-Stokes convergence test to new
staggered
---
.../freeflowporousmedium/1p_1p/CMakeLists.txt | 2 +
.../1p_1p/convergence/CMakeLists.txt | 50 +++
.../1p_1p/convergence/analyticalsolutions.hh | 340 +++++++++++++++
.../1p_1p/convergence/convergencetest.py | 124 ++++++
.../1p_1p/convergence/main.cc | 320 ++++++++++++++
.../convergence/manufactured_solution.py | 93 +++++
.../1p_1p/convergence/params.input | 34 ++
.../1p_1p/convergence/problem_darcy.hh | 271 ++++++++++++
.../1p_1p/convergence/problem_freeflow.hh | 391 ++++++++++++++++++
.../1p_1p/convergence/properties.hh | 129 ++++++
.../1p_1p/convergence/spatialparams.hh | 129 ++++++
.../1p_1p/convergence/testcase.hh | 37 ++
12 files changed, 1920 insertions(+)
create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/CMakeLists.txt
create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/analyticalsolutions.hh
create mode 100755 test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/convergencetest.py
create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/main.cc
create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/manufactured_solution.py
create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/params.input
create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_darcy.hh
create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_freeflow.hh
create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/properties.hh
create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/spatialparams.hh
create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/testcase.hh
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/CMakeLists.txt b/test/multidomain/boundary/freeflowporousmedium/1p_1p/CMakeLists.txt
index da2b6d53d6..617ade1705 100644
--- a/test/multidomain/boundary/freeflowporousmedium/1p_1p/CMakeLists.txt
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/CMakeLists.txt
@@ -1,3 +1,5 @@
+add_subdirectory(convergence)
+
add_input_file_links()
add_executable(test_md_boundary_ff1p_pm1p EXCLUDE_FROM_ALL main.cc)
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/CMakeLists.txt b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/CMakeLists.txt
new file mode 100644
index 0000000000..8336a0b30d
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/CMakeLists.txt
@@ -0,0 +1,50 @@
+add_input_file_links()
+dune_symlink_to_source_files(FILES "convergencetest.py")
+
+add_executable(test_md_boundary_ff1p_pm1p_convtest EXCLUDE_FROM_ALL main.cc)
+
+dumux_add_test(NAME test_md_boundary_ff1p_pm1p_convtest_shiue1
+ TARGET test_md_boundary_ff1p_pm1p_convtest
+ LABELS multidomain multidomain_boundary stokesdarcy
+ TIMEOUT 1000
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ./convergencetest.py
+ CMD_ARGS test_md_boundary_ff1p_pm1p_convtest params.input
+ -Problem.TestCase ShiueExampleOne
+ -Darcy.Problem.BoundaryConditions Dirichlet
+ -FreeFlow.Problem.BoundaryConditions Dirichlet)
+
+dumux_add_test(NAME test_md_boundary_ff1p_pm1p_convtest_shiue2
+ TARGET test_md_boundary_ff1p_pm1p_convtest
+ LABELS multidomain multidomain_boundary stokesdarcy
+ TIMEOUT 1000
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ./convergencetest.py
+ CMD_ARGS test_md_boundary_ff1p_pm1p_convtest params.input
+ -Problem.TestCase ShiueExampleTwo
+ -Darcy.Problem.BoundaryConditions Mixed
+ -FreeFlow.Problem.BoundaryConditions Mixed)
+
+dumux_add_test(NAME test_md_boundary_ff1p_pm1p_convtest_rybak
+ TARGET test_md_boundary_ff1p_pm1p_convtest
+ LABELS multidomain multidomain_boundary stokesdarcy
+ TIMEOUT 1000
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ./convergencetest.py
+ CMD_ARGS test_md_boundary_ff1p_pm1p_convtest params.input
+ -Problem.TestCase Rybak
+ -Darcy.Problem.BoundaryConditions Neumann
+ -FreeFlow.Problem.BoundaryConditions Stress)
+
+dumux_add_test(NAME test_md_boundary_ff1p_pm1p_convtest_navierstokes
+ TARGET test_md_boundary_ff1p_pm1p_convtest
+ LABELS multidomain multidomain_boundary stokesdarcy
+ TIMEOUT 1000
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ./convergencetest.py
+ CMD_ARGS test_md_boundary_ff1p_pm1p_convtest params.input
+ -Problem.TestCase Schneider
+ -FreeFlow.Problem.EnableInertiaTerms true
+ -FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph true
+ -Darcy.Problem.BoundaryConditions Mixed
+ -FreeFlow.Problem.BoundaryConditions Mixed)
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/analyticalsolutions.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/analyticalsolutions.hh
new file mode 100644
index 0000000000..6df9efadd3
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/analyticalsolutions.hh
@@ -0,0 +1,340 @@
+// -*- 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 Analytical solutions and rhs for the different test cases
+ */
+
+#ifndef DUMUX_CONVERGENCE_TEST_ANALYTICAL_SOLUTIONS_HH
+#define DUMUX_CONVERGENCE_TEST_ANALYTICAL_SOLUTIONS_HH
+
+#include <cmath>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+namespace Dumux::Solution::DarcyStokes {
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// see Rybak et al., 2015: "Multirate time integration for coupled
+// saturated/unsaturated porous medium and free flow systems"
+/////////////////////////////////////////////////////////////////////////////////////////////////
+namespace Rybak {
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> darcy(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ using std::exp; using std::sin; using std::cos;
+ sol[0] = -0.5*M_PI*y*y*cos(M_PI*x);
+ sol[1] = -1.0*y*sin(M_PI*x);
+ sol[2] = 0.5*y*y*sin(M_PI*x);
+ return sol;
+}
+
+// 0: mass
+Dune::FieldVector<double, 1> darcyRHS(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ using std::sin;
+ return { (0.5*M_PI*y*M_PI*y - 1)*sin(M_PI*x) };
+}
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> stokes(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ using std::sin; using std::cos;
+ sol[0] = -cos(M_PI*x)*sin(M_PI*y);
+ sol[1] = sin(M_PI*x)*cos(M_PI*y);
+ sol[2] = 0.5*y*sin(M_PI*x);
+ return sol;
+}
+
+// 0: mom-x | 1: mom-y | 2: mass
+Dune::FieldVector<double, 3> stokesRHS(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ Dune::FieldVector<double, 3> source(0.0);
+ using std::sin; using std::cos;
+ source[0] = 0.5*M_PI*y*cos(M_PI*x) - 2*M_PI*M_PI*cos(M_PI*x)*sin(M_PI*y);
+ source[1] = 0.5*sin(M_PI*x) + 2*M_PI*M_PI*sin(M_PI*x)*cos(M_PI*y);
+ return source;
+}
+
+// sigma_ff = -sym(grad(v)) + pI
+Dune::FieldMatrix<double, 2, 2> stokesStress(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ using std::exp; using std::sin; using std::cos;
+
+ Dune::FieldMatrix<double, 2, 2> stress(0.0);
+ stress[0][0] = (0.5*y - 2*M_PI*sin(M_PI*y))*sin(M_PI*x);
+ stress[1][0] = 0.0;
+ stress[0][1] = stress[1][0]; // symmetric
+ stress[1][1] = (0.5*y - 2*M_PI*sin(M_PI*y))*sin(M_PI*x);
+ return stress;
+}
+
+} // end namespace Rybak
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
+// Example 1
+/////////////////////////////////////////////////////////////////////////////////////////////////
+namespace ShiueOne {
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> darcy(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ using std::exp; using std::sin; using std::cos;
+ sol[0] = M_PI*(-exp(1)*y + exp(y))*sin(M_PI*x);
+ sol[1] = (exp(1) - exp(y))*cos(M_PI*x);
+ sol[2] = (exp(y) - y*exp(1)) * cos(M_PI*x);
+
+ return sol;
+}
+
+// 0: mass
+Dune::FieldVector<double, 1> darcyRHS(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ using std::exp; using std::sin; using std::cos;
+ return { cos(M_PI*x) * (M_PI*M_PI*(exp(y) - y*exp(1)) - exp(y)) };
+}
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> stokes(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ using std::exp; using std::sin; using std::cos;
+ sol[0] = -1/M_PI * exp(y) * sin(M_PI*x);
+ sol[1] = (exp(y) - exp(1)) * cos(M_PI*x);
+ sol[2] = 2*exp(y) * cos(M_PI*x);
+ return sol;
+}
+
+// 0: mom-x | 1: mom-y | 2: mass
+Dune::FieldVector<double, 3> stokesRHS(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ using std::exp; using std::sin; using std::cos;
+ Dune::FieldVector<double, 3> source(0.0);
+ source[0] = exp(y)*sin(M_PI*x) * (1/M_PI -3*M_PI);
+ source[1] = cos(M_PI*x) * (M_PI*M_PI*(exp(y)- exp(1)) + exp(y));
+ return source;
+}
+
+// sigma_ff = -sym(grad(v)) + pI
+Dune::FieldMatrix<double, 2, 2> stokesStress(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ using std::exp; using std::sin; using std::cos;
+
+ Dune::FieldMatrix<double, 2, 2> stress(0.0);
+ stress[0][0] = 4.0*exp(y)*cos(M_PI*x);
+ stress[1][0] = (M_PI*M_PI*(exp(y) - exp(1)) + 1.0*exp(y))*sin(M_PI*x)/M_PI;
+ stress[0][1] = stress[1][0]; // symmetric
+ stress[1][1] = 0.0;
+ return stress;
+}
+
+} // end namespace ShiueExOne
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
+// Example 2
+/////////////////////////////////////////////////////////////////////////////////////////////////
+namespace ShiueTwo {
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> darcy(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ sol[0] = x*(y - 1.0) + (x - 1.0)*(y - 1.0) - 2.0;
+ sol[1] = x*(x - 1.0) - 1.0*(y - 1.0)*(y - 1.0) - 2.0;
+ sol[2] = 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;
+
+ return sol;
+}
+
+// 0: mass
+Dune::FieldVector<double, 1> darcyRHS(const Dune::FieldVector<double, 2>& globalPos)
+{ return { 0.0 }; }
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> stokes(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ sol[0] = (y-1.0)*(y-1.0) + x*(y-1.0) + 3.0*x - 1.0;
+ sol[1] = x*(x-1.0) - 0.5*(y-1.0)*(y-1.0) - 3.0*y + 1.0;
+ sol[2] = 2.0*x + y - 1.0;
+ return sol;
+}
+
+// 0: mom-x | 1: mom-y | 2: mass
+Dune::FieldVector<double, 3> stokesRHS(const Dune::FieldVector<double, 2>& globalPos)
+{ return { 0.0, 0.0, 0.0 }; }
+
+// sigma_ff = -sym(grad(v)) + pI
+Dune::FieldMatrix<double, 2, 2> stokesStress(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+
+ Dune::FieldMatrix<double, 2, 2> stress(0.0);
+ stress[0][0] = 2.0*x - y - 5.0;
+ stress[1][0] = -3*x - 2*y + 3.0;
+ stress[0][1] = stress[1][0]; // symmetric
+ stress[1][1] = 2.0*x + 3.0*y + 3.0;
+ return stress;
+}
+
+} // end namespace ShiueTwo
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// see Schneider et al., 2019: "Coupling staggered-grid and MPFA finite volume methods for
+// free flow/porous-medium flow problems (with c = 0)"
+/////////////////////////////////////////////////////////////////////////////////////////////////
+namespace Schneider {
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> darcy(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ static constexpr double omega = M_PI;
+ static constexpr double c = 0.0;
+ using std::exp; using std::sin; using std::cos;
+ const double sinOmegaX = sin(omega*x);
+ const double cosOmegaX = cos(omega*x);
+ static const double expTwo = exp(2);
+ const double expYPlusOne = exp(y+1);
+
+ sol[0] = c/(2*omega)*expYPlusOne*sinOmegaX*sinOmegaX
+ -omega*(expYPlusOne + 2 - expTwo)*cosOmegaX;
+ sol[1] = (0.5*c*(expYPlusOne + 2 - expTwo)*cosOmegaX
+ -(c*cosOmegaX + 1)*exp(y-1))*sinOmegaX;
+ sol[2] = (expYPlusOne + 2 - expTwo)*sinOmegaX;
+
+ return sol;
+}
+
+// 0: mass
+Dune::FieldVector<double, 1> darcyRHS(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ using std::exp; using std::sin; using std::cos;
+ static constexpr double omega = M_PI;
+ static constexpr double c = 0.0;
+ const double cosOmegaX = cos(omega*x);
+ static const double expTwo = exp(2);
+ const double expYPlusOne = exp(y+1);
+
+ return { sin(omega*x)*(-(c*cosOmegaX + 1)*exp(y - 1)
+ + 1.5*c*expYPlusOne*cosOmegaX
+ + omega*omega*(expYPlusOne - expTwo + 2)) };
+}
+
+// 0: velocity-x | 1: velocity-y | 2: pressure
+Dune::FieldVector<double, 3> stokes(const Dune::FieldVector<double, 2>& globalPos)
+{
+ Dune::FieldVector<double, 3> sol(0.0);
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ using std::sin;
+ static constexpr double omega = M_PI;
+ const double sinOmegaX = sin(omega*x);
+
+ sol[0] = y;
+ sol[1] = -y*sinOmegaX;
+ sol[2] = -y*y*sinOmegaX*sinOmegaX;
+ return sol;
+}
+
+// 0: mom-x | 1: mom-y | 2: mass
+Dune::FieldVector<double, 3> stokesRHS(const Dune::FieldVector<double, 2>& globalPos)
+{
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ using std::exp; using std::sin; using std::cos;
+ static constexpr double omega = M_PI;
+ const double sinOmegaX = sin(omega*x);
+ const double cosOmegaX = cos(omega*x);
+
+ Dune::FieldVector<double, 3> source(0.0);
+ source[0] = -2*omega*y*y*sinOmegaX*cosOmegaX
+ -2*y*sinOmegaX + omega*cosOmegaX;
+ source[1] = -omega*y*y*cosOmegaX - omega*omega*y*sinOmegaX;
+ source[2] = -sinOmegaX;
+ return source;
+}
+
+// sigma_ff = vvT + sym(grad(v)) + pI
+Dune::FieldMatrix<double, 2, 2> stokesStress(const Dune::FieldVector<double, 2>& globalPos)
+{
+ using std::sin; using std::cos;
+ const double x = globalPos[0];
+ const double y = globalPos[1];
+ static constexpr double omega = M_PI;
+ const double sinOmegaX = sin(omega*x);
+ const double cosOmegaX = cos(omega*x);
+
+ Dune::FieldMatrix<double, 2, 2> stress(0.0);
+ stress[0][0] = y*y * cosOmegaX*cosOmegaX;
+ stress[0][1] = -y*y*sinOmegaX + omega*y*cosOmegaX - 1.0;
+ stress[1][0] = stress[0][1]; // symmetric
+ stress[1][1] = 2*sinOmegaX;
+ return stress;
+}
+
+} // end namespace Schneider
+
+} // end namespace Dumux::Solution::DarcyStokes
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/convergencetest.py b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/convergencetest.py
new file mode 100755
index 0000000000..b52826c3c4
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/convergencetest.py
@@ -0,0 +1,124 @@
+#!/usr/bin/env python3
+
+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)
+
+executableName = str(sys.argv[1])
+testargs = [str(i) for i in sys.argv][2:]
+testname = testargs[testargs.index('-Problem.TestCase')+1]
+
+# remove the old log files
+subprocess.call(['rm', testname + '_freeFlow.log'])
+print("Removed old log file ({})!".format(testname + '_freeFlow.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(['./' + executableName] + testargs + ['-Grid.Refinement', str(i)])
+
+def checkRatesFreeFlow():
+ # check the rates and append them to the log file
+ with open(testname + '_freeFlow.log', "r+") as logfile:
+
+ errorP, errorVx, errorVy = [], [], []
+ for line in logfile:
+ line = line.strip("\n").strip("\[ConvergenceTest\]").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)
+
+ print("\nComputed the following convergence rates for {}:\n".format(testname))
+
+ subprocess.call(['cat', testname + '_freeFlow.log'])
+
+ return {"p" : resultsP, "v_x" : resultsVx, "v_y" : resultsVy}
+
+def checkRatesDarcy():
+ # check the rates and append them to the log file
+ with open(testname + '_darcy.log', "r+") as logfile:
+
+ 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)
+
+ print("\nComputed the following convergence rates for {}:\n".format(testname))
+
+ subprocess.call(['cat', testname + '_darcy.log'])
+
+ return {"p" : resultsP}
+
+def checkRatesFreeFlowAndDarcy():
+ resultsFreeFlow = checkRatesFreeFlow()
+ resultsDarcy = checkRatesDarcy()
+
+ def mean(numbers):
+ return float(sum(numbers)) / len(numbers)
+
+ # check the rates, we expect rates around 2
+ if mean(resultsFreeFlow["p"]) < 2.05 and mean(resultsFreeFlow["p"]) < 1.7:
+ 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(resultsFreeFlow["v_x"]) < 2.05 and mean(resultsFreeFlow["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(resultsFreeFlow["v_y"]) < 2.05 and mean(resultsFreeFlow["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)
+
+
+checkRatesFreeFlowAndDarcy()
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/main.cc b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/main.cc
new file mode 100644
index 0000000000..8a9a0044d8
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/main.cc
@@ -0,0 +1,320 @@
+// -*- 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 convergence test for the coupled FreeFlow/Darcy problem (1p).
+ */
+
+#include <config.h>
+
+#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/dumuxmessage.hh>
+#include <dumux/common/partial.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/io/format.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager_yasp.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+#include <dumux/freeflow/navierstokes/velocityoutput.hh>
+
+#include <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
+#include <test/freeflow/navierstokes/errors.hh>
+
+#include "testcase.hh"
+#include "properties.hh"
+
+/*!
+* \brief Creates analytical solution vector for output
+* Returns a tuple of the analytical solution for the pressure and velocity at cell centers
+* \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.fullAnalyticalSolution(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 MassP, class MomP, class SolutionVector>
+void printFreeFlowL2Error(std::shared_ptr<MomP> momentumProblem,
+ std::shared_ptr<MassP> massProblem,
+ const SolutionVector& sol)
+{
+ using namespace Dumux;
+
+ NavierStokesTest::Errors errors(momentumProblem, massProblem, sol);
+ const int numCellCenterDofs = massProblem->gridGeometry().numDofs();
+ const int numFaceDofs = momentumProblem->gridGeometry().numDofs();
+ const auto absL2 = errors.l2Absolute();
+ const auto relL2 = errors.l2Relative();
+
+ std::cout << Fmt::format("** L2 error (abs/rel) for {} cc dofs and {} face dofs (total: {}): ",
+ numCellCenterDofs, numFaceDofs, numCellCenterDofs + numFaceDofs)
+ << Fmt::format("L2(p) = {:.8e} / {:.8e}", absL2[0], relL2[0])
+ << Fmt::format(", L2(vx) = {:.8e} / {:.8e}", absL2[1], relL2[1])
+ << Fmt::format(", L2(vy) = {:.8e} / {:.8e}", absL2[2], relL2[2])
+ << std::endl;
+
+ // write the norm into a log file
+ std::ofstream logFile(massProblem->name() + ".log", std::ios::app);
+ logFile << "[ConvergenceTest] L2(p) = " << absL2[0]
+ << " L2(vx) = " << absL2[1]
+ << " L2(vy) = " << absL2[2]
+ << std::endl;
+}
+
+template<class Problem, class SolutionVector>
+void printDarcyL2Error(std::shared_ptr<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->fullAnalyticalSolution(scv.center())[2/*pressureIdx*/];
+ l2error += scv.volume()*(delta*delta);
+ }
+ }
+ using std::sqrt;
+ l2error = sqrt(l2error);
+
+ const auto numDofs = problem->gridGeometry().numDofs();
+ std::cout << Fmt::format("** L2 error (abs) for {} cc dofs", numDofs)
+ << Fmt::format("L2 error = {:.8e}", l2error)
+ << std::endl;
+
+ // write the norm into a log file
+ std::ofstream logFile(problem->name() + ".log", std::ios::app);
+ logFile << "[ConvergenceTest] L2(p) = " << l2error << std::endl;
+}
+
+int main(int argc, char** argv)
+{
+ 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 FreeFlowMomentumTypeTag = Properties::TTag::FreeFlowOnePMomentum;
+ using FreeFlowMassTypeTag = Properties::TTag::FreeFlowOnePMass;
+ 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 FreeFlowGridManager = Dumux::GridManager<GetPropType<FreeFlowMomentumTypeTag, Properties::Grid>>;
+ FreeFlowGridManager freeFlowGridManager;
+ freeFlowGridManager.init("FreeFlow"); // pass parameter group
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridManager.grid().leafGridView();
+ const auto& freeFlowGridView = freeFlowGridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FreeFlowMomentumGridGeometry = GetPropType<FreeFlowMomentumTypeTag, Properties::GridGeometry>;
+ auto freeFlowMomentumGridGeometry = std::make_shared<FreeFlowMomentumGridGeometry>(freeFlowGridView);
+ using FreeFlowMassGridGeometry = GetPropType<FreeFlowMassTypeTag, Properties::GridGeometry>;
+ auto freeFlowMassGridGeometry = std::make_shared<FreeFlowMassGridGeometry>(freeFlowGridView);
+ using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::GridGeometry>;
+ auto darcyGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+
+ using Traits = MultiDomainTraits<FreeFlowMomentumTypeTag, FreeFlowMassTypeTag, DarcyTypeTag>;
+ using CouplingManager = FreeFlowPorousMediumCouplingManager<Traits>;
+ auto couplingManager = std::make_shared<CouplingManager>();
+
+ // the indices
+ constexpr auto freeFlowMomentumIndex = CouplingManager::freeFlowMomentumIndex;
+ constexpr auto freeFlowMassIndex = CouplingManager::freeFlowMassIndex;
+ constexpr auto porousMediumIndex = CouplingManager::porousMediumIndex;
+
+ // the problem (initial and boundary conditions)
+ const auto testCaseName = getParam<std::string>("Problem.TestCase");
+ const auto testCase = [&]()
+ {
+ if (testCaseName == "ShiueExampleOne")
+ return DarcyStokesTestCase::ShiueExampleOne;
+ else if (testCaseName == "ShiueExampleTwo")
+ return DarcyStokesTestCase::ShiueExampleTwo;
+ else if (testCaseName == "Rybak")
+ return DarcyStokesTestCase::Rybak;
+ else if (testCaseName == "Schneider")
+ return DarcyStokesTestCase::Schneider;
+ else
+ DUNE_THROW(Dune::InvalidStateException, testCaseName + " is not a valid test case");
+ }();
+
+ // the problem (initial and boundary conditions)
+ using FreeFlowMomentumProblem = GetPropType<FreeFlowMomentumTypeTag, Properties::Problem>;
+ auto freeFlowMomentumProblem = std::make_shared<FreeFlowMomentumProblem>(freeFlowMomentumGridGeometry, couplingManager, testCase, testCaseName);
+ using FreeFlowMassProblem = GetPropType<FreeFlowMassTypeTag, Properties::Problem>;
+ auto freeFlowMassProblem = std::make_shared<FreeFlowMassProblem>(freeFlowMassGridGeometry, couplingManager, testCase, testCaseName);
+ using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
+ auto spatialParams = std::make_shared<typename DarcyProblem::SpatialParams>(darcyGridGeometry, testCase);
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyGridGeometry, couplingManager, spatialParams, testCase, testCaseName);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[freeFlowMomentumIndex].resize(freeFlowMomentumGridGeometry->numDofs());
+ sol[freeFlowMassIndex].resize(freeFlowMassGridGeometry->numDofs());
+ sol[porousMediumIndex].resize(darcyGridGeometry->numDofs());
+
+ // the grid variables
+ using FreeFlowMomentumGridVariables = GetPropType<FreeFlowMomentumTypeTag, Properties::GridVariables>;
+ auto freeFlowMomentumGridVariables = std::make_shared<FreeFlowMomentumGridVariables>(freeFlowMomentumProblem, freeFlowMomentumGridGeometry);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
+ using FreeFlowMassGridVariables = GetPropType<FreeFlowMassTypeTag, Properties::GridVariables>;
+ auto freeFlowMassGridVariables = std::make_shared<FreeFlowMassGridVariables>(freeFlowMassProblem, freeFlowMassGridGeometry);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyGridGeometry);
+
+ couplingManager->init(freeFlowMomentumProblem, freeFlowMassProblem, darcyProblem,
+ std::make_tuple(freeFlowMomentumGridVariables, freeFlowMassGridVariables, darcyGridVariables),
+ sol);
+
+ freeFlowMomentumGridVariables->init(sol[freeFlowMomentumIndex]);
+ freeFlowMassGridVariables->init(sol[freeFlowMassIndex]);
+ darcyGridVariables->init(sol[porousMediumIndex]);
+
+ // initialize the vtk output module
+ VtkOutputModule freeflowVtkWriter(*freeFlowMassGridVariables, sol[freeFlowMassIndex], freeFlowMassProblem->name());
+ GetPropType<FreeFlowMassTypeTag, Properties::IOFields>::initOutputModule(freeflowVtkWriter); // Add model specific output fields
+ freeflowVtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<FreeFlowMassGridVariables>>());
+ NavierStokesTest::AnalyticalSolutionVectors freeFlowAnalyticSol(freeFlowMomentumProblem, freeFlowMassProblem);
+ freeflowVtkWriter.addField(freeFlowAnalyticSol.analyticalPressureSolution(), "pressureExact");
+ freeflowVtkWriter.addField(freeFlowAnalyticSol.analyticalVelocitySolution(), "velocityExact");
+ freeflowVtkWriter.write(0.0);
+
+ VtkOutputModule darcyVtkWriter(*darcyGridVariables, sol[porousMediumIndex], darcyProblem->name());
+ GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
+ darcyVtkWriter.addVelocityOutput(std::make_shared<GetPropType<DarcyTypeTag, Properties::VelocityOutput>>(*darcyGridVariables));
+ const auto darcyAnalyticalSolution = createDarcyAnalyticalSolution<double>(*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(freeFlowMomentumProblem, freeFlowMassProblem, darcyProblem),
+ std::make_tuple(freeFlowMomentumGridGeometry,
+ freeFlowMassGridGeometry,
+ darcyGridGeometry),
+ std::make_tuple(freeFlowMomentumGridVariables,
+ freeFlowMassGridVariables,
+ 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);
+
+ const auto dofsVS = freeFlowMomentumGridGeometry->numDofs();
+ const auto dofsPS = freeFlowMassGridGeometry->numDofs();
+ const auto dofsPD = darcyGridGeometry->numDofs();
+ std::cout << "Stokes velocity dofs: " << dofsVS << std::endl;
+ std::cout << "Stokes pressure dofs: " << dofsPS << std::endl;
+ std::cout << "Darcy pressure dofs: " << dofsPD << std::endl;
+ std::cout << "Total number of dofs: " << dofsVS + dofsPS + dofsPD << std::endl;
+
+ // solve the non-linear system
+ nonLinearSolver.solve(sol);
+
+ // write vtk output
+ freeflowVtkWriter.write(1.0);
+ darcyVtkWriter.write(1.0);
+
+ // print L2 errors
+ printFreeFlowL2Error(freeFlowMomentumProblem, freeFlowMassProblem, partial(sol, freeFlowMomentumIndex, freeFlowMassIndex));
+ printDarcyL2Error(darcyProblem, sol[porousMediumIndex]);
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+}
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/manufactured_solution.py b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/manufactured_solution.py
new file mode 100644
index 0000000000..7bc744363b
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/manufactured_solution.py
@@ -0,0 +1,93 @@
+#!/usr/bin/env python3
+
+"""
+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 = "Schneider"
+
+# 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)
+ elif case == "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
+ else: # Schneider
+ vFF = np.array([y, -y*sp.sin(sp.pi*x)])
+ pFF = -y*y*sp.sin(sp.pi*x)*sp.sin(sp.pi*x)
+ 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
+if case == "Schneider":
+ momentumFlux = vvT - (gradV + gradVT) +pI
+else:
+ momentumFlux = - (gradV + gradVT) +pI # only Stokes
+divMomentumFlux = div(momentumFlux)
+
+print("\n\nStress tensor flux for case", case)
+print("-- [0][0]: ", sp.simplify(momentumFlux[0][0]))
+print("-- [1][1]: ", sp.simplify(momentumFlux[1][1]))
+print("-- [0][1] == [1][0]: ", sp.simplify(momentumFlux[0][1]))
+
+print("\n\nSource 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)
+ elif case == "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
+ else: # Schneider
+ pD = (sp.exp(y+1) + 2.0 - sp.exp(2))*sp.sin(sp.pi*x)
+ 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/freeflowporousmedium/1p_1p/convergence/params.input b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/params.input
new file mode 100644
index 0000000000..6d656742ed
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/params.input
@@ -0,0 +1,34 @@
+[Darcy.Grid]
+UpperRight = 1 1
+Cells = 40 40
+
+[FreeFlow.Grid]
+LowerLeft = 0 1
+UpperRight = 1 2
+Cells = 40 40
+
+[FreeFlow.Problem]
+Name = freeFlow
+EnableInertiaTerms = false
+
+[Darcy.Problem]
+Name = darcy
+
+[Darcy.SpatialParams]
+AlphaBeaversJoseph = 1.0
+
+[Problem]
+EnableGravity = false
+
+[Vtk]
+AddVelocity = 1
+
+[Component]
+LiquidDensity = 1.0
+LiquidKinematicViscosity = 1.0
+
+[Assembly]
+NumericDifference.BaseEpsilon = 1e-4
+
+[FreeFlow.Flux]
+UpwindWeight = 0.5
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_darcy.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_darcy.hh
new file mode 100644
index 0000000000..a575666dbe
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_darcy.hh
@@ -0,0 +1,271 @@
+// -*- 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 <dumux/common/boundarytypes.hh>
+#include <dumux/common/numeqvector.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
+
+#include "testcase.hh"
+#include "analyticalsolutions.hh"
+
+namespace Dumux {
+
+/*!
+ * \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 Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
+ using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using GridView = typename GridGeometry::GridView;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+
+ enum class BC {
+ dirichlet, neumann, mixed
+ };
+
+public:
+ //! export the Indices
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
+ DarcySubProblem(std::shared_ptr<const GridGeometry> gridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager,
+ std::shared_ptr<typename ParentType::SpatialParams> spatialParams,
+ const DarcyStokesTestCase testCase,
+ const std::string& name)
+ : ParentType(gridGeometry, spatialParams, "Darcy")
+ , couplingManager_(couplingManager)
+ , testCase_(testCase)
+ {
+ problemName_ = name + "_"
+ + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
+
+ auto bc = getParamFromGroup<std::string>(this->paramGroup(), "Problem.BoundaryConditions", "Dirichlet");
+ if (bc == "Dirichlet")
+ boundaryConditions_ = BC::dirichlet;
+ else if (bc == "Neumann")
+ boundaryConditions_ = BC::neumann;
+ else if (bc == "Mixed")
+ boundaryConditions_ = BC::mixed;
+ else
+ DUNE_THROW(Dune::Exception, "Wrong BC type choose: Dirichlet, Neumann or Mixed");
+
+ std::cout << "Porous medium domain: Using " << bc << " boundary conditions" << std::endl;
+ }
+
+ const std::string& name() const
+ { return problemName_; }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ 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().isCoupled(CouplingManager::porousMediumIndex, CouplingManager::freeFlowMassIndex, scvf))
+ values.setAllCouplingNeumann();
+ else
+ {
+ if (boundaryConditions_ == BC::dirichlet)
+ values.setAllDirichlet();
+ else if (boundaryConditions_ == BC::neumann)
+ values.setAllNeumann();
+ else
+ {
+ if (onLeftBoundary_(scvf.center()))
+ values.setAllNeumann();
+ 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 sub-control-volume-face
+ */
+ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ const auto p = fullAnalyticalSolution(scvf.center())[2];
+ 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
+ */
+ 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().isCoupled(CouplingManager::porousMediumIndex, CouplingManager::freeFlowMassIndex, scvf))
+ values[Indices::conti0EqIdx] = couplingManager().massCouplingCondition(
+ CouplingManager::porousMediumIndex, CouplingManager::freeFlowMassIndex,
+ fvGeometry, scvf, elemVolVars
+ );
+
+ else
+ {
+ const auto sol = fullAnalyticalSolution(scvf.center());
+ const auto n = scvf.unitOuterNormal();
+ auto v = n; v[0] = sol[0]; v[1] = sol[1];
+ values[Indices::conti0EqIdx] = v*n;
+ }
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ /*!
+ * \brief Evaluates the source term for all phases within a given
+ * sub control volume.
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
+ {
+ using namespace Solution::DarcyStokes;
+ switch (testCase_)
+ {
+ case DarcyStokesTestCase::ShiueExampleOne:
+ return ShiueOne::darcyRHS(globalPos);
+ case DarcyStokesTestCase::ShiueExampleTwo:
+ return ShiueTwo::darcyRHS(globalPos);
+ case DarcyStokesTestCase::Rybak:
+ return Rybak::darcyRHS(globalPos);
+ case DarcyStokesTestCase::Schneider:
+ return Schneider::darcyRHS(globalPos);
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
+ }
+ }
+
+ // \}
+
+ /*!
+ * \brief Evaluates the initial value for a control volume.
+ * \param element The element
+ */
+ 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
+ * Returns vector with entries: (velocity-x | velocity-y | pressure)
+ */
+ Dune::FieldVector<Scalar, 3> fullAnalyticalSolution(const GlobalPosition& globalPos) const
+ {
+ using namespace Solution::DarcyStokes;
+ switch (testCase_)
+ {
+ case DarcyStokesTestCase::ShiueExampleOne:
+ return ShiueOne::darcy(globalPos);
+ case DarcyStokesTestCase::ShiueExampleTwo:
+ return ShiueTwo::darcy(globalPos);
+ case DarcyStokesTestCase::Rybak:
+ return Rybak::darcy(globalPos);
+ case DarcyStokesTestCase::Schneider:
+ return Schneider::darcy(globalPos);
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
+ }
+ }
+
+ // \}
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+private:
+ bool onLeftBoundary_(const GlobalPosition& globalPos) const
+ { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
+
+ static constexpr Scalar eps_ = 1e-7;
+ std::shared_ptr<CouplingManager> couplingManager_;
+ std::string problemName_;
+ DarcyStokesTestCase testCase_;
+ BC boundaryConditions_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_freeflow.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_freeflow.hh
new file mode 100644
index 0000000000..9ef9fd3db9
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_freeflow.hh
@@ -0,0 +1,391 @@
+// -*- 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 free-flow sub-problem of coupled FreeFlow/Darcy convergence test
+ */
+
+#ifndef DUMUX_FREEFLOW_SUBPROBLEM_HH
+#define DUMUX_FREEFLOW_SUBPROBLEM_HH
+
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/freeflow/navierstokes/momentum/fluxhelper.hh>
+#include <dumux/freeflow/navierstokes/scalarfluxhelper.hh>
+#include <dune/common/fmatrix.hh>
+
+#include "testcase.hh"
+#include "analyticalsolutions.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup BoundaryTests
+ * \brief The Stokes sub-problem of coupled Stokes-Darcy convergence test
+ */
+template <class TypeTag>
+class FreeFlowSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using GridView = typename GridGeometry::GridView;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using PrimaryVariables = typename ParentType::PrimaryVariables;
+ using NumEqVector = typename ParentType::NumEqVector;
+ using BoundaryTypes = typename ParentType::BoundaryTypes;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+
+ enum class BC {
+ dirichlet, stress, mixed
+ };
+
+public:
+ //! export the Indices
+ using Indices = typename ModelTraits::Indices;
+
+ FreeFlowSubProblem(std::shared_ptr<const GridGeometry> gridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager,
+ const DarcyStokesTestCase testCase,
+ const std::string& name)
+ : ParentType(gridGeometry, couplingManager, "FreeFlow")
+ , couplingManager_(couplingManager)
+ , testCase_(testCase)
+ {
+ problemName_ = name + "_"
+ + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
+
+ auto bc = getParamFromGroup<std::string>(this->paramGroup(), "Problem.BoundaryConditions", "Dirichlet");
+ if (bc == "Dirichlet")
+ boundaryConditions_ = BC::dirichlet;
+ else if (bc == "Stress")
+ boundaryConditions_ = BC::stress;
+ else if (bc == "Mixed")
+ boundaryConditions_ = BC::mixed;
+ else
+ DUNE_THROW(Dune::Exception, "Wrong BC type choose: Dirichlet, Stress or Mixed");
+
+ std::cout << "Free flow domain: Using " << bc << " boundary conditions" << std::endl;
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ const std::string& name() const
+ { return problemName_; }
+
+ /*!
+ * \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
+
+ // \}
+
+ /*!
+ * \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;
+
+ if constexpr (ParentType::isMomentumProblem())
+ {
+ if (couplingManager().isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf))
+ {
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setCouplingNeumann(Indices::momentumXBalanceIdx);
+ }
+ else
+ {
+ if (boundaryConditions_ == BC::dirichlet)
+ values.setAllDirichlet();
+ else if (boundaryConditions_ == BC::stress)
+ values.setAllNeumann();
+ else
+ {
+ if (onLeftBoundary_(scvf.center()))
+ values.setAllNeumann();
+ else
+ values.setAllDirichlet();
+ }
+ }
+ }
+ else
+ {
+ if (couplingManager().isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
+ values.setAllCouplingNeumann();
+ else
+ values.setAllNeumann();
+ }
+
+ 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 ElementFluxVariablesCache>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+ const auto& globalPos = scvf.ipGlobal();
+ using FluxHelper = NavierStokesMomentumBoundaryFluxHelper;
+
+ // momentum boundary conditions
+ if constexpr (ParentType::isMomentumProblem())
+ {
+ // We need to take care here: If the integration point of an scvf lies both on the coupling interface and the
+ // domain boundary (i.e., the leftmost and rightmost points of the interface), do not evaluate the coupling or slip condition
+ // because they would require data from the boundary not available in this case (velocities for evaluating gradients,
+ // those would only be available for Dirichlet BCs). Instead, directly use the given Neumann boundary stress.
+ // TODO: Maybe couplingManager().isCoupled(...) could return false for these scvfs.
+ //
+ // | <-- left Neumann boundary
+ // |
+ // integration point --> o##### <-- scvf lying on coupling interface with integration point touching left Neumann boundary
+ //
+ if (couplingManager().isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf)
+ && scvf.ipGlobal()[0] > this->gridGeometry().bBoxMin()[0] + eps_
+ && scvf.ipGlobal()[0] < this->gridGeometry().bBoxMax()[0] - eps_)
+ {
+ values += couplingManager().momentumCouplingCondition(
+ CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex,
+ fvGeometry, scvf, elemVolVars
+ );
+
+ values += FluxHelper::slipVelocityMomentumFlux(
+ *this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache
+ );
+ }
+ else
+ {
+ const auto stress = stressTensorAtPos(globalPos);
+ const auto n = scvf.unitOuterNormal();
+ stress.mv(n, values);
+ }
+ }
+
+ // mass boundary conditions
+ else
+ {
+ if (couplingManager().isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf))
+ {
+ values = couplingManager().massCouplingCondition(
+ CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex,
+ fvGeometry, scvf, elemVolVars
+ );
+ }
+ else
+ {
+ using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>;
+ values = FluxHelper::scalarOutflowFlux(
+ *this, element, fvGeometry, scvf, elemVolVars
+ );
+ }
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Returns true if the scvf lies on a porous slip boundary
+ */
+ bool onSlipBoundary(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
+ {
+ return scvf.boundary()
+ && couplingManager().isCoupled(
+ CouplingManager::freeFlowMomentumIndex,
+ CouplingManager::porousMediumIndex,
+ scvf
+ );
+ }
+
+ // \}
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Returns the sources within the domain.
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+ {
+ const auto select = [&](const auto& rhs) -> NumEqVector {
+ if constexpr (ParentType::isMomentumProblem())
+ return { rhs[0], rhs[1] };
+ else
+ return { rhs[2] };
+ };
+
+ using namespace Solution::DarcyStokes;
+ switch (testCase_)
+ {
+ case DarcyStokesTestCase::ShiueExampleOne:
+ return select(ShiueOne::stokesRHS(globalPos));
+ case DarcyStokesTestCase::ShiueExampleTwo:
+ return select(ShiueTwo::stokesRHS(globalPos));
+ case DarcyStokesTestCase::Rybak:
+ return select(Rybak::stokesRHS(globalPos));
+ case DarcyStokesTestCase::Schneider:
+ return select(Schneider::stokesRHS(globalPos));
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
+ }
+ }
+
+ /*!
+ * \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
+ */
+ auto permeability(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
+ { return couplingManager().darcyPermeability(fvGeometry, scvf); }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the
+ Beavers-Joseph-Saffman boundary condition.
+ */
+ Scalar alphaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
+ { return couplingManager().problem(CouplingManager::porousMediumIndex).spatialParams().beaversJosephCoeffAtPos(scvf.ipGlobal()); }
+
+ /*!
+ * \brief Returns the analytical solution of the problem at a given position.
+ * \param globalPos The global position
+ * Returns vector with entries: (velocity-x | velocity-y | pressure)
+ */
+ Dune::FieldVector<Scalar, 3> fullAnalyticalSolution(const GlobalPosition& globalPos) const
+ {
+ using namespace Solution::DarcyStokes;
+ switch (testCase_)
+ {
+ case DarcyStokesTestCase::ShiueExampleOne:
+ return ShiueOne::stokes(globalPos);
+ case DarcyStokesTestCase::ShiueExampleTwo:
+ return ShiueTwo::stokes(globalPos);
+ case DarcyStokesTestCase::Rybak:
+ return Rybak::stokes(globalPos);
+ case DarcyStokesTestCase::Schneider:
+ return Schneider::stokes(globalPos);
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
+ }
+ }
+
+ /*!
+ * \brief Returns the analytical solution of the problem at a given position.
+ * \param globalPos The global position
+ * Returns analytical solution depending on the type of sub-problem
+ */
+ PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
+ {
+ using namespace Solution::DarcyStokes;
+ const auto sol = fullAnalyticalSolution(globalPos);
+ if constexpr (ParentType::isMomentumProblem())
+ return { sol[0], sol[1] };
+ else
+ return { sol[2] };
+ }
+
+ // \}
+
+private:
+ /*!
+ * \brief Returns the stress tensor within the domain.
+ * \param globalPos The global position
+ */
+ Dune::FieldMatrix<double, 2, 2> stressTensorAtPos(const GlobalPosition &globalPos) const
+ {
+ using namespace Solution::DarcyStokes;
+ switch (testCase_)
+ {
+ case DarcyStokesTestCase::ShiueExampleOne:
+ return ShiueOne::stokesStress(globalPos);
+ case DarcyStokesTestCase::ShiueExampleTwo:
+ return ShiueTwo::stokesStress(globalPos);
+ case DarcyStokesTestCase::Rybak:
+ return Rybak::stokesStress(globalPos);
+ case DarcyStokesTestCase::Schneider:
+ return Schneider::stokesStress(globalPos);
+ default:
+ DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
+ }
+ }
+
+ bool onLeftBoundary_(const GlobalPosition& globalPos) const
+ { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
+
+ static constexpr Scalar eps_ = 1e-7;
+ std::string problemName_;
+ std::shared_ptr<CouplingManager> couplingManager_;
+ DarcyStokesTestCase testCase_;
+ BC boundaryConditions_;
+};
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/properties.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/properties.hh
new file mode 100644
index 0000000000..70c62ded2f
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/properties.hh
@@ -0,0 +1,129 @@
+// -*- 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 properties for a simple Darcy test (cell-centered finite volume method)
+ */
+#ifndef DUMUX_DARCYSTOKES_PROPERTIES_HH
+#define DUMUX_DARCYSTOKES_PROPERTIES_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/material/components/constant.hh>
+
+#include <dumux/porousmediumflow/1p/model.hh>
+#include <dumux/freeflow/navierstokes/mass/1p/model.hh>
+#include <dumux/freeflow/navierstokes/momentum/model.hh>
+#include <dumux/discretization/cctpfa.hh>
+#include <dumux/discretization/fcstaggered.hh>
+
+#include <dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh>
+#include <dumux/multidomain/traits.hh>
+
+#include "spatialparams.hh"
+#include "problem_darcy.hh"
+#include "problem_freeflow.hh"
+
+namespace Dumux::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 = DarcySubProblem<TypeTag>; };
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyOneP>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, 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 = ConvergenceTestSpatialParams<GridGeometry, Scalar>;
+};
+
+// Create new type tags
+namespace TTag {
+struct FreeFlowOneP {};
+struct FreeFlowOnePMass { using InheritsFrom = std::tuple<FreeFlowOneP, NavierStokesMassOneP, CCTpfaModel>; };
+struct FreeFlowOnePMomentum { using InheritsFrom = std::tuple<FreeFlowOneP, NavierStokesMomentum, FaceCenteredStaggeredModel>; };
+} // end namespace TTag
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::FreeFlowOneP>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::FreeFlowOneP>
+{
+ using Coords = Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2>;
+ using type = Dune::YaspGrid<2, Coords>;
+};
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::FreeFlowOneP>
+{ using type = FreeFlowSubProblem<TypeTag>; };
+
+template<class TypeTag>
+struct EnableGridGeometryCache<TypeTag, TTag::FreeFlowOneP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeFlowOneP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeFlowOneP> { static constexpr bool value = true; };
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyOneP>
+{
+ using Traits = MultiDomainTraits<TTag::FreeFlowOnePMomentum, TTag::FreeFlowOnePMass, TTag::DarcyOneP>;
+ using type = FreeFlowPorousMediumCouplingManager<Traits>;
+};
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::FreeFlowOneP>
+{
+ using Traits = MultiDomainTraits<TTag::FreeFlowOnePMomentum, TTag::FreeFlowOnePMass, TTag::DarcyOneP>;
+ using type = FreeFlowPorousMediumCouplingManager<Traits>;
+};
+
+} // end namespace Dumux::Properties
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/spatialparams.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/spatialparams.hh
new file mode 100644
index 0000000000..82eb797c14
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/spatialparams.hh
@@ -0,0 +1,129 @@
+// -*- 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 spatial parameters class for the test problem using the 1p cc model.
+ */
+
+#ifndef DUMUX_1P_TEST_SPATIALPARAMS_HH
+#define DUMUX_1P_TEST_SPATIALPARAMS_HH
+
+#include <dumux/material/spatialparams/fv1p.hh>
+#include "testcase.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup BoundaryTests
+ * \brief The spatial parameters class for the test problem using the
+ * 1p cc model.
+ */
+template<class GridGeometry, class Scalar>
+class ConvergenceTestSpatialParams
+: public FVSpatialParamsOneP<GridGeometry, Scalar,
+ ConvergenceTestSpatialParams<GridGeometry, Scalar>>
+{
+ using GridView = typename GridGeometry::GridView;
+ using ParentType = FVSpatialParamsOneP<
+ GridGeometry, Scalar, ConvergenceTestSpatialParams<GridGeometry, Scalar>
+ >;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ static constexpr int dimWorld = GridView::dimensionworld;
+ using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
+
+
+public:
+ // export permeability type
+ using PermeabilityType = DimWorldMatrix;
+
+ ConvergenceTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry,
+ const DarcyStokesTestCase testCase)
+ : ParentType(gridGeometry)
+ , testCase_(testCase)
+ {
+ alphaBJ_ = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
+ }
+
+ /*!
+ * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
+ *
+ * \param element The element
+ * \param scv The sub control volume
+ * \param elemSol The element solution vector
+ * \return the intrinsic permeability
+ */
+ template<class SubControlVolume, class ElementSolution>
+ PermeabilityType permeability(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
+ {
+ PermeabilityType K(0.0);
+
+ if (testCase_ == DarcyStokesTestCase::Schneider)
+ {
+ using std::cos;
+ using std::sin;
+ using std::exp;
+
+ static constexpr Scalar c = 0.0;
+ static constexpr Scalar omega = M_PI;
+
+ const Scalar x = scv.center()[0];
+ K[0][0] = 1.0;
+ K[0][1] = -c/(2*omega) * sin(omega*x);
+ K[1][0] = K[0][1];
+ K[1][1] = exp(-2)*(1 + c*cos(omega*x));
+ }
+ else
+ {
+ const static Scalar permeability = getParam<Scalar>("Darcy.SpatialParams.Permeability", 1.0);
+ K[0][0] = permeability;
+ K[1][1] = permeability;
+ }
+
+ return K;
+ }
+
+ /*! \brief Defines the porosity in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return 0.4; }
+
+ /*! \brief Defines the Beavers-Joseph coefficient in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
+ { return alphaBJ_; }
+
+private:
+ DarcyStokesTestCase testCase_;
+ Scalar permeability_;
+ Scalar alphaBJ_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/testcase.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/testcase.hh
new file mode 100644
index 0000000000..9496ffe6ae
--- /dev/null
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/testcase.hh
@@ -0,0 +1,37 @@
+// -*- 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 different test cases
+ */
+
+#ifndef DUMUX_CONVERGENCE_TEST_TESTCASE_HH
+#define DUMUX_CONVERGENCE_TEST_TESTCASE_HH
+
+namespace Dumux {
+
+enum class DarcyStokesTestCase
+{
+ ShiueExampleOne, ShiueExampleTwo, Rybak, Schneider
+};
+
+} // end namespace Dumux
+
+#endif
--
GitLab