Commit 45925f19 authored by Timo Koch's avatar Timo Koch
Browse files

[md][boundary] Port Darcy-Stokes convergence test to new staggered

parent 3a7688af
Pipeline #8820 passed with stages
add_subdirectory(convergence)
add_input_file_links()
add_executable(test_md_boundary_ff1p_pm1p EXCLUDE_FROM_ALL main.cc)
......
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)
// -*- 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
#!/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()
// -*- 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