Commit 921757a2 authored by Martin Schneider's avatar Martin Schneider
Browse files

Merge branch 'feature/navier-stokes-darcy-convtest' into 'master'

Feature/navier stokes darcy convtest

See merge request !2128
parents 141fe13d d8932fc1
......@@ -25,6 +25,7 @@
#define DUMUX_NAVIERSTOKES_PROBLEM_HH
#include <dune/common/exceptions.hh>
#include <dune/common/typetraits.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/staggeredfvproblem.hh>
#include <dumux/discretization/method.hh>
......@@ -204,8 +205,19 @@ public:
}
/*!
* \brief Returns the beta value, or the alpha value divided by the square root of the intrinsic permeability.
* \brief Returns the beta value which is the alpha value divided by the square root of the (scalar-valued) interface permeability.
*/
Scalar betaBJ(const Element& element, const SubControlVolumeFace& scvf, const GlobalPosition& tangentialVector) const
{
const Scalar interfacePermeability = interfacePermeability_(element, scvf, tangentialVector);
using std::sqrt;
return asImp_().alphaBJ(scvf) / sqrt(interfacePermeability);
}
/*!
* \brief Returns the beta value which is the alpha value divided by the square root of the (scalar-valued) interface permeability.
*/
[[deprecated("Use betaBJ with tangential vector instead. Will be removed after 3.3")]]
Scalar betaBJ(const Element& element, const SubControlVolumeFace& scvf) const
{
using std::sqrt;
......@@ -220,7 +232,10 @@ public:
return VelocityVector(0.0);
}
//! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph condition is used
/*!
* \brief Returns the slip velocity at a porous boundary based on the Beavers-Joseph(-Saffman) condition.
*/
const Scalar beaversJosephVelocity(const Element& element,
const SubControlVolume& scv,
const SubControlVolumeFace& ownScvf,
......@@ -228,21 +243,33 @@ public:
const Scalar velocitySelf,
const Scalar tangentialVelocityGradient) const
{
// du/dy + dv/dx = alpha/sqrt(K) * (u_boundary-uPM)
// beta = alpha/sqrt(K)
const Scalar betaBJ = asImp_().betaBJ(element, faceOnPorousBoundary);
const Scalar distanceNormalToBoundary = (faceOnPorousBoundary.center() - scv.center()).two_norm();
// create a unit normal vector oriented in positive coordinate direction
GlobalPosition orientation = ownScvf.unitOuterNormal();
orientation[ownScvf.directionIndex()] = 1.0;
// du/dy + dv/dx = alpha/sqrt(K) * (u_boundary-uPM)
// beta = alpha/sqrt(K)
const Scalar betaBJ = asImp_().betaBJ(element, faceOnPorousBoundary, orientation);
const Scalar distanceNormalToBoundary = (faceOnPorousBoundary.center() - scv.center()).two_norm();
return (tangentialVelocityGradient*distanceNormalToBoundary
+ asImp_().porousMediumVelocity(element, faceOnPorousBoundary) * orientation * betaBJ * distanceNormalToBoundary
+ velocitySelf) / (betaBJ*distanceNormalToBoundary + 1.0);
}
private:
//! Returns a scalar permeability value at the coupling interface
Scalar interfacePermeability_(const Element& element, const SubControlVolumeFace& scvf, const GlobalPosition& tangentialVector) const
{
const auto& K = asImp_().permeability(element, scvf);
// use t*K*t for permeability tensors
if constexpr (Dune::IsNumber<std::decay_t<decltype(K)>>::value)
return K;
else
return vtmv(tangentialVector, K, tangentialVector);
}
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
{ return *static_cast<Implementation *>(this); }
......
......@@ -283,7 +283,7 @@ public:
/*!
* \brief Returns the intrinsic permeability of the coupled Darcy element.
*/
Scalar darcyPermeability(const Element<stokesIdx>& element, const SubControlVolumeFace<stokesIdx>& scvf) const
auto darcyPermeability(const Element<stokesIdx>& element, const SubControlVolumeFace<stokesIdx>& scvf) const
{
const auto& stokesContext = couplingManager().stokesCouplingContext(element, scvf);
return stokesContext.volVars.permeability();
......
add_input_file_links()
dune_symlink_to_source_files(FILES "convergencetest.py")
add_executable(test_md_boundary_darcy1p_freeflow1p_convtest EXCLUDE_FROM_ALL main.cc)
dune_add_test(NAME test_md_boundary_darcy1p_stokes1p_convtest
SOURCES main.cc
TARGET test_md_boundary_darcy1p_freeflow1p_convtest
LABELS multidomain multidomain_boundary stokesdarcy
TIMEOUT 1000
CMAKE_GUARD HAVE_UMFPACK
COMMAND ./convergencetest.py
CMD_ARGS test_md_boundary_darcy1p_freeflow1p_convtest params.input
-Darcy.SpatialParams.Permeability 1.0)
dune_add_test(NAME test_md_boundary_darcy1p_navierstokes1p_convtest
TARGET test_md_boundary_darcy1p_freeflow1p_convtest
LABELS multidomain multidomain_boundary stokesdarcy
TIMEOUT 1000
CMAKE_GUARD HAVE_UMFPACK
COMMAND ./convergencetest.py
CMD_ARGS test_md_boundary_darcy1p_stokes1p_convtest params.input)
CMD_ARGS test_md_boundary_darcy1p_freeflow1p_convtest params.input
-Problem.TestCase Schneider
-FreeFlow.Problem.EnableInertiaTerms true
-FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph true)
......@@ -12,8 +12,8 @@ testname = str(sys.argv[1])
testargs = [str(i) for i in sys.argv][2:]
# remove the old log files
subprocess.call(['rm', testname + '_stokes.log'])
print("Removed old log file ({})!".format(testname + '_stokes.log'))
subprocess.call(['rm', testname + '_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'))
......@@ -22,9 +22,9 @@ for i in [0, 1, 2]:
subprocess.call(['./' + testname] + testargs + ['-Grid.Refinement', str(i),
'-Vtk.OutputName', testname])
def checkRatesStokes():
def checkRatesFreeFlow():
# check the rates and append them to the log file
logfile = open(testname + '_stokes.log', "r+")
logfile = open(testname + '_freeFlow.log', "r+")
errorP = []
errorVx = []
......@@ -64,7 +64,7 @@ def checkRatesStokes():
logfile.close()
print("\nComputed the following convergence rates for {}:\n".format(testname))
subprocess.call(['cat', testname + '_stokes.log'])
subprocess.call(['cat', testname + '_freeFlow.log'])
return {"p" : resultsP, "v_x" : resultsVx, "v_y" : resultsVy}
......@@ -104,23 +104,23 @@ def checkRatesDarcy():
return {"p" : resultsP}
def checkRatesStokesAndDarcy():
resultsStokes = checkRatesStokes()
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(resultsStokes["p"]) < 2.05 and mean(resultsStokes["p"]) < 1.84:
if mean(resultsFreeFlow["p"]) < 2.05 and mean(resultsFreeFlow["p"]) < 1.84:
sys.stderr.write("*"*70 + "\n" + "The convergence rates for pressure were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
sys.exit(1)
if mean(resultsStokes["v_x"]) < 2.05 and mean(resultsStokes["v_x"]) < 1.95:
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(resultsStokes["v_y"]) < 2.05 and mean(resultsStokes["v_y"]) < 1.95:
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)
......@@ -129,4 +129,4 @@ def checkRatesStokesAndDarcy():
sys.exit(1)
checkRatesStokesAndDarcy()
checkRatesFreeFlowAndDarcy()
......@@ -19,9 +19,7 @@
/*!
* \file
* \ingroup BoundaryTests
* \brief A test problem for the coupled Stokes/Darcy problem (1p).
* The analytical solution is given in Shiue et al., 2018:
* "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem"
* \brief A test problem for the coupled FreeFlow/Darcy problem (1p).
*/
#include <config.h>
......@@ -51,6 +49,7 @@
#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
#include "testcase.hh"
#include "problem_darcy.hh"
#include "problem_stokes.hh"
......@@ -58,7 +57,7 @@ namespace Dumux {
namespace Properties {
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::StokesOneP>
struct CouplingManager<TypeTag, TTag::FreeFlowOneP>
{
using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOneP>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
......@@ -67,7 +66,7 @@ struct CouplingManager<TypeTag, TTag::StokesOneP>
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::DarcyOneP>
{
using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesOneP, Properties::TTag::StokesOneP, TypeTag>;
using Traits = StaggeredMultiDomainTraits<Properties::TTag::FreeFlowOneP, Properties::TTag::FreeFlowOneP, TypeTag>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
......@@ -80,7 +79,7 @@ struct CouplingManager<TypeTag, TTag::DarcyOneP>
* \param problem the problem for which to evaluate the analytical solution
*/
template<class Scalar, class Problem>
auto createStokesAnalyticalSolution(const Problem& problem)
auto createFreeFlowAnalyticalSolution(const Problem& problem)
{
const auto& gridGeometry = problem.gridGeometry();
using GridView = typename std::decay_t<decltype(gridGeometry)>::GridView;
......@@ -171,11 +170,11 @@ auto createDarcyAnalyticalSolution(const Problem& problem)
}
template<class Problem, class SolutionVector>
void printStokesL2Error(const Problem& problem, const SolutionVector& x)
void printFreeFlowL2Error(const Problem& problem, const SolutionVector& x)
{
using namespace Dumux;
using Scalar = double;
using TypeTag = Properties::TTag::StokesOneP;
using TypeTag = Properties::TTag::FreeFlowOneP;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
......@@ -253,7 +252,7 @@ int main(int argc, char** argv) try
Parameters::init(argc, argv);
// Define the sub problem type tags
using StokesTypeTag = Properties::TTag::StokesOneP;
using FreeFlowTypeTag = Properties::TTag::FreeFlowOneP;
using DarcyTypeTag = Properties::TTag::DarcyOneP;
// try to create a grid (from the given grid file or the input file)
......@@ -262,67 +261,83 @@ int main(int argc, char** argv) try
DarcyGridManager darcyGridManager;
darcyGridManager.init("Darcy"); // pass parameter group
using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
StokesGridManager stokesGridManager;
stokesGridManager.init("Stokes"); // pass parameter group
using FreeFlowGridManager = Dumux::GridManager<GetPropType<FreeFlowTypeTag, 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& stokesGridView = stokesGridManager.grid().leafGridView();
const auto& freeFlowGridView = freeFlowGridManager.grid().leafGridView();
// create the finite volume grid geometry
using StokesGridGeometry = GetPropType<StokesTypeTag, Properties::GridGeometry>;
auto stokesGridGeometry = std::make_shared<StokesGridGeometry>(stokesGridView);
stokesGridGeometry->update();
using FreeFlowGridGeometry = GetPropType<FreeFlowTypeTag, Properties::GridGeometry>;
auto freeFlowGridGeometry = std::make_shared<FreeFlowGridGeometry>(freeFlowGridView);
freeFlowGridGeometry->update();
using DarcyGridGeometry = GetPropType<DarcyTypeTag, Properties::GridGeometry>;
auto darcyGridGeometry = std::make_shared<DarcyGridGeometry>(darcyGridView);
darcyGridGeometry->update();
using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
using Traits = StaggeredMultiDomainTraits<FreeFlowTypeTag, FreeFlowTypeTag, DarcyTypeTag>;
// the coupling manager
using CouplingManager = StokesDarcyCouplingManager<Traits>;
auto couplingManager = std::make_shared<CouplingManager>(stokesGridGeometry, darcyGridGeometry);
auto couplingManager = std::make_shared<CouplingManager>(freeFlowGridGeometry, darcyGridGeometry);
// the indices
constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
constexpr auto freeFlowCellCenterIdx = CouplingManager::stokesCellCenterIdx;
constexpr auto freeFlowFaceIdx = CouplingManager::stokesFaceIdx;
constexpr auto darcyIdx = CouplingManager::darcyIdx;
// the problem (initial and boundary conditions)
using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
auto stokesProblem = std::make_shared<StokesProblem>(stokesGridGeometry, couplingManager);
const auto testCase = []()
{
const auto testCaseInput = getParam<std::string>("Problem.TestCase", "ShiueExampleTwo");
if (testCaseInput == "ShiueExampleOne")
return TestCase::ShiueExampleOne;
else if (testCaseInput == "ShiueExampleTwo")
return TestCase::ShiueExampleTwo;
else if (testCaseInput == "Rybak")
return TestCase::Rybak;
else if (testCaseInput == "Schneider")
return TestCase::Schneider;
else
DUNE_THROW(Dune::InvalidStateException, testCaseInput + " is not a valid test case");
}();
using FreeFlowProblem = GetPropType<FreeFlowTypeTag, Properties::Problem>;
auto freeFlowProblem = std::make_shared<FreeFlowProblem>(freeFlowGridGeometry, couplingManager, testCase);
using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
auto darcyProblem = std::make_shared<DarcyProblem>(darcyGridGeometry, couplingManager);
auto spatialParams = std::make_shared<typename DarcyProblem::SpatialParams>(darcyGridGeometry, testCase);
auto darcyProblem = std::make_shared<DarcyProblem>(darcyGridGeometry, couplingManager, spatialParams, testCase);
// the solution vector
Traits::SolutionVector sol;
sol[stokesCellCenterIdx].resize(stokesGridGeometry->numCellCenterDofs());
sol[stokesFaceIdx].resize(stokesGridGeometry->numFaceDofs());
sol[freeFlowCellCenterIdx].resize(freeFlowGridGeometry->numCellCenterDofs());
sol[freeFlowFaceIdx].resize(freeFlowGridGeometry->numFaceDofs());
sol[darcyIdx].resize(darcyGridGeometry->numDofs());
// get a solution vector storing references to the two Stokes solution vectors
auto stokesSol = partial(sol, stokesFaceIdx, stokesCellCenterIdx);
// get a solution vector storing references to the two FreeFlow solution vectors
auto freeFlowSol = partial(sol, freeFlowFaceIdx, freeFlowCellCenterIdx);
couplingManager->init(stokesProblem, darcyProblem, sol);
couplingManager->init(freeFlowProblem, darcyProblem, sol);
// the grid variables
using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesGridGeometry);
stokesGridVariables->init(stokesSol);
using FreeFlowGridVariables = GetPropType<FreeFlowTypeTag, Properties::GridVariables>;
auto freeFlowGridVariables = std::make_shared<FreeFlowGridVariables>(freeFlowProblem, freeFlowGridGeometry);
freeFlowGridVariables->init(freeFlowSol);
using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyGridGeometry);
darcyGridVariables->init(sol[darcyIdx]);
// intialize the vtk output module
using Scalar = typename Traits::Scalar;
StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name());
GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
const auto stokesAnalyticalSolution = createStokesAnalyticalSolution<Scalar>(*stokesProblem);
stokesVtkWriter.addField(std::get<0>(stokesAnalyticalSolution), "pressureExact");
stokesVtkWriter.addField(std::get<1>(stokesAnalyticalSolution), "velocityExact");
stokesVtkWriter.addFaceField(std::get<2>(stokesAnalyticalSolution), "faceVelocityExact");
stokesVtkWriter.write(0.0);
StaggeredVtkOutputModule<FreeFlowGridVariables, decltype(freeFlowSol)> freeFlowVtkWriter(*freeFlowGridVariables, freeFlowSol, freeFlowProblem->name());
GetPropType<FreeFlowTypeTag, Properties::IOFields>::initOutputModule(freeFlowVtkWriter);
const auto freeFlowAnalyticalSolution = createFreeFlowAnalyticalSolution<Scalar>(*freeFlowProblem);
freeFlowVtkWriter.addField(std::get<0>(freeFlowAnalyticalSolution), "pressureExact");
freeFlowVtkWriter.addField(std::get<1>(freeFlowAnalyticalSolution), "velocityExact");
freeFlowVtkWriter.addFaceField(std::get<2>(freeFlowAnalyticalSolution), "faceVelocityExact");
freeFlowVtkWriter.write(0.0);
VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyProblem->name());
using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
......@@ -335,12 +350,12 @@ int main(int argc, char** argv) try
// the assembler for a stationary problem
using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
std::make_tuple(stokesGridGeometry->faceFVGridGeometryPtr(),
stokesGridGeometry->cellCenterFVGridGeometryPtr(),
auto assembler = std::make_shared<Assembler>(std::make_tuple(freeFlowProblem, freeFlowProblem, darcyProblem),
std::make_tuple(freeFlowGridGeometry->faceFVGridGeometryPtr(),
freeFlowGridGeometry->cellCenterFVGridGeometryPtr(),
darcyGridGeometry),
std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
stokesGridVariables->cellCenterGridVariablesPtr(),
std::make_tuple(freeFlowGridVariables->faceGridVariablesPtr(),
freeFlowGridVariables->cellCenterGridVariablesPtr(),
darcyGridVariables),
couplingManager);
......@@ -356,10 +371,10 @@ int main(int argc, char** argv) try
nonLinearSolver.solve(sol);
// write vtk output
stokesVtkWriter.write(1.0);
freeFlowVtkWriter.write(1.0);
darcyVtkWriter.write(1.0);
printStokesL2Error(*stokesProblem, stokesSol);
printFreeFlowL2Error(*freeFlowProblem, freeFlowSol);
printDarcyL2Error(*darcyProblem, sol[darcyIdx]);
////////////////////////////////////////////////////////////
......
......@@ -2,24 +2,23 @@
UpperRight = 1 1
Cells = 40 40
[Stokes.Grid]
[FreeFlow.Grid]
LowerLeft = 0 1
UpperRight = 1 2
Cells = 40 40
[Stokes.Problem]
Name = stokes
[FreeFlow.Problem]
Name = freeFlow
EnableInertiaTerms = false
[Darcy.Problem]
Name = darcy
[Darcy.SpatialParams]
Permeability = 1.0
AlphaBeaversJoseph = 1.0
[Vtk]
OutputName = test_md_boundary_stokes1p_darcy1p_convergencetest
OutputName = test_md_boundary_freeflow1p_darcy1p_convergencetest
[Problem]
EnableGravity = false
......@@ -33,3 +32,6 @@ LiquidKinematicViscosity = 1.0
[Assembly]
NumericDifference.BaseEpsilon = 1e-4
[FreeFlow.Flux]
UpwindWeight = 0.5
......@@ -33,7 +33,8 @@
#include <dumux/porousmediumflow/1p/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include "../spatialparams.hh"
#include "spatialparams.hh"
#include "testcase.hh"
#include <dumux/material/components/constant.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
......@@ -69,7 +70,7 @@ struct SpatialParams<TypeTag, TTag::DarcyOneP>
{
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = OnePSpatialParams<GridGeometry, Scalar>;
using type = ConvergenceTestSpatialParams<GridGeometry, Scalar>;
};
} // end namespace Properties
......@@ -100,30 +101,19 @@ class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
static constexpr auto velocityYIdx = 1;
static constexpr auto pressureIdx = 2;
enum class TestCase
{
ShiueExampleOne, ShiueExampleTwo, Rybak
};
public:
//! export the Indices
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
DarcySubProblem(std::shared_ptr<const GridGeometry> gridGeometry,
std::shared_ptr<CouplingManager> couplingManager)
: ParentType(gridGeometry, "Darcy"), couplingManager_(couplingManager)
std::shared_ptr<CouplingManager> couplingManager,
std::shared_ptr<typename ParentType::SpatialParams> spatialParams,
const TestCase testCase)
: ParentType(gridGeometry, spatialParams, "Darcy")
, couplingManager_(couplingManager)
, testCase_(testCase)
{
problemName_ = getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
const auto testCaseInput = getParamFromGroup<std::string>(this->paramGroup(), "Problem.TestCase", "ShiueExampleTwo");
if (testCaseInput == "ShiueExampleOne")
testCase_ = TestCase::ShiueExampleOne;
else if (testCaseInput == "ShiueExampleTwo")
testCase_ = TestCase::ShiueExampleTwo;
else if (testCaseInput == "Rybak")
testCase_ = TestCase::Rybak;
else
DUNE_THROW(Dune::InvalidStateException, testCaseInput + " is not a valid test case");
}
/*!
......@@ -233,6 +223,8 @@ public:
return rhsShiueEtAlExampleTwo_(globalPos);
case TestCase::Rybak:
return rhsRybak_(globalPos);
case TestCase::Schneider:
return rhsSchneiderEtAl_(globalPos);
default:
DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
}
......@@ -268,6 +260,8 @@ public:
return analyticalSolutionShiueEtAlExampleTwo_(globalPos);
case TestCase::Rybak:
return analyticalSolutionRybak_(globalPos);
case TestCase::Schneider:
return analyticalSolutionSchneiderEtAl_(globalPos);
default:
DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
}
......@@ -346,6 +340,50 @@ private:
NumEqVector rhsShiueEtAlExampleTwo_(const GlobalPosition& globalPos) const
{ return NumEqVector(0.0); }
// see Schneider et al., 2019: "Coupling staggered-grid and MPFA finite volume methods for
// free flow/porous-medium flow problems"
Dune::FieldVector<Scalar, 3> analyticalSolutionSchneiderEtAl_(const GlobalPosition& globalPos) const
{
Dune::FieldVector<Scalar, 3> sol(0.0);
const Scalar x = globalPos[0];
const Scalar y = globalPos[1];
static constexpr Scalar omega = M_PI;
static constexpr Scalar c = 0.0;
using std::exp; using std::sin; using std::cos;
const Scalar sinOmegaX = sin(omega*x);
const Scalar cosOmegaX = cos(omega*x);
static const Scalar expTwo = exp(2);
const Scalar expYPlusOne = exp(y+1);
sol[pressureIdx] = (expYPlusOne + 2 - expTwo)*sinOmegaX;
sol[velocityXIdx] = c/(2*omega)*expYPlusOne*sinOmegaX*sinOmegaX
-omega*(expYPlusOne + 2 - expTwo)*cosOmegaX;
sol[velocityYIdx] = (0.5*c*(expYPlusOne + 2 - expTwo)*cosOmegaX
-(c*cosOmegaX + 1)*exp(y-1))*sinOmegaX;
return sol;
}
// see Schneider et al., 2019: "Coupling staggered-grid and MPFA finite volume methods for
// free flow/porous-medium flow problems (with c = 0)"
NumEqVector rhsSchneiderEtAl_(const GlobalPosition& globalPos) const
{
const Scalar x = globalPos[0];
const Scalar y = globalPos[1];
using std::exp; using std::sin; using std::cos;
static constexpr Scalar omega = M_PI;
static constexpr Scalar c = 0.0;
const Scalar cosOmegaX = cos(omega*x);
static const Scalar expTwo = exp(2);
const Scalar expYPlusOne = exp(y+1);
const Scalar result = (-(c*cosOmegaX + 1)*exp(y - 1)
+ 1.5*c*expYPlusOne*cosOmegaX
+ omega*omega*(expYPlusOne - expTwo + 2))
*sin(omega*x);
return NumEqVector(result);
}
static constexpr Scalar eps_ = 1e-7;
std::shared_ptr<CouplingManager> couplingManager_;
std::string problemName_;
......
......@@ -19,11 +19,11 @@
/*!
* \file
* \ingroup BoundaryTests
* \brief The Stokes sub-problem of coupled Stokes-Darcy convergence test
* \brief The free-flow sub-problem of coupled FreeFlow/Darcy convergence test
*/
#ifndef DUMUX_STOKES_SUBPROBLEM_HH
#define DUMUX_STOKES_SUBPROBLEM_HH
#ifndef DUMUX_FREEFLOW_SUBPROBLEM_HH
#define DUMUX_FREEFLOW_SUBPROBLEM_HH
#include <dune/common/fvector.hh>
#include <dune/grid/yaspgrid.hh>
......@@ -34,20 +34,21 @@
#include <dumux/freeflow/navierstokes/problem.hh>
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/freeflow/navierstokes/model.hh>
#include "testcase.hh"