Commit 1bf7d0dd authored by Martin Schneider's avatar Martin Schneider
Browse files

[md][ffpm][tests] Incorporate newIC tests into existing convergence test

parent 10e39298
......@@ -47,4 +47,36 @@ dune_add_test(NAME test_md_boundary_darcy1pbox_navierstokes1p_convtest
CMD_ARGS test_md_boundary_darcy1pbox_freeflow1p_convtest params.input
-Problem.TestCase Schneider
-FreeFlow.Problem.EnableInertiaTerms true
-FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph true)
\ No newline at end of file
-FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph true)
dune_add_test(NAME test_md_boundary_darcy1pbox_stokes1p_convtest_bj
TARGET test_md_boundary_darcy1pbox_freeflow1p_convtest
LABELS multidomain multidomain_boundary stokesdarcy
TIMEOUT 1000
CMAKE_GUARD HAVE_UMFPACK
COMMAND ./convergencetest.py
CMD_ARGS test_md_boundary_darcy1pbox_freeflow1p_convtest params.input
-Problem.TestCase BJSymmetrized
-Problem.SlipCondition BJ
-Darcy.SpatialParams.Porosity 1.0
-Darcy.SpatialParams.Permeability "1.0 1.0 2.0")
dune_add_test(NAME test_md_boundary_darcy1p_stokes1p_newic
TARGET test_md_boundary_darcy1pbox_freeflow1p_convtest
LABELS multidomain multidomain_boundary stokesdarcy
TIMEOUT 1000
CMAKE_GUARD HAVE_UMFPACK
COMMAND ./convergencetest.py
CMD_ARGS test_md_boundary_darcy1pbox_freeflow1p_convtest params.input
-Problem.TestCase NewICNonSymmetrized
-Problem.SlipCondition ER
-Problem.NewIc true
-FreeFlow.EnableUnsymmetrizedVelocityGradientForIC true
-FreeFlow.EnableUnsymmetrizedVelocityGradient true
-Darcy.InterfaceParams.EpsInterface 1.0
-Darcy.InterfaceParams.N_s_bl 0.5
-Darcy.InterfaceParams.N_1_bl -3.0
-Darcy.InterfaceParams.M_bl "1.0 1.0"
-Darcy.SpatialParams.Porosity 1.0
-Darcy.SpatialParams.Permeability "1.0 1.0 0.0"
-Component.LiquidKinematicViscosity 2.0)
\ No newline at end of file
......@@ -129,7 +129,7 @@ void printDarcyL2Error(const Problem& problem, const SolutionVector& x)
for (auto&& scv : scvs(fvGeometry))
{
const auto dofIdx = scv.dofIndex();
const Scalar delta = x[dofIdx] - problem.analyticalSolution(scv.center())[2/*pressureIdx*/];
const Scalar delta = x[dofIdx] - problem.analyticalSolution(scv.dofPosition())[2/*pressureIdx*/];
l2error += scv.volume()*(delta*delta);
}
}
......@@ -212,6 +212,10 @@ int main(int argc, char** argv)
return TestCase::Rybak;
else if (testCaseName == "Schneider")
return TestCase::Schneider;
else if (testCaseName == "BJSymmetrized")
return TestCase::BJSymmetrized;
else if (testCaseName == "NewICNonSymmetrized")
return TestCase::NewICNonSymmetrized;
else
DUNE_THROW(Dune::InvalidStateException, testCaseName + " is not a valid test case");
}();
......
......@@ -29,9 +29,9 @@
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/cctpfa.hh>
#include <dumux/discretization/box.hh>
#include <dumux/common/boundarytypes.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
......@@ -46,6 +46,39 @@ namespace Dumux {
template <class TypeTag>
class DarcySubProblem;
namespace Properties {
// Create new type tags
namespace TTag {
struct DarcyOneP { using InheritsFrom = std::tuple<OneP>; };
struct DarcyOnePBox { using InheritsFrom = std::tuple<DarcyOneP, BoxModel>; };
struct DarcyOnePCC { using InheritsFrom = std::tuple<DarcyOneP, CCTpfaModel>; };
} // end namespace TTag
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::DarcyOneP> { using type = Dumux::DarcySubProblem<TypeTag>; };
// the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::DarcyOneP>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::Constant<1, Scalar> > ;
};
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::DarcyOneP> { using type = Dune::YaspGrid<2>; };
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::DarcyOneP>
{
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = ConvergenceTestSpatialParams<GridGeometry, Scalar>;
};
} // end namespace Properties
/*!
* \ingroup BoundaryTests
* \brief The Darcy sub-problem of coupled Stokes-Darcy convergence test
......@@ -57,7 +90,7 @@ class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
......@@ -134,17 +167,48 @@ public:
return values;
}
/*!
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary control volume.
*
* \param element The element
* \param scv The boundary sub control volume
*/
BoundaryTypes boundaryTypes(const Element &element, const SubControlVolume &scv) const
{
BoundaryTypes values;
values.setAllDirichlet();
//corners will not be handled as coupled
if(onLeftBoundary_(scv.dofPosition()) || onRightBoundary_(scv.dofPosition()))
{
values.setAllDirichlet();
return values;
}
auto fvGeometry = localView(this->gridGeometry());
fvGeometry.bindElement(element);
for (auto&& scvf : scvfs(fvGeometry))
{
if (couplingManager().isCoupledEntity(CouplingManager::porousMediumIdx, element, scvf))
{
values.setAllCouplingNeumann();
}
}
return values;
}
/*!
* \brief Evaluates the boundary conditions for a Dirichlet control volume.
*
* \param element The element for which the Dirichlet boundary condition is set
* \param scvf The boundary subcontrolvolumeface
* \param globalPos The position for which the Dirichlet value is set
*
* For this method, the \a values parameter stores primary variables.
*/
PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
const auto p = analyticalSolution(scvf.center())[pressureIdx];
const auto p = analyticalSolution(globalPos)[pressureIdx];
return PrimaryVariables(p);
}
......@@ -198,6 +262,10 @@ public:
return rhsRybak_(globalPos);
case TestCase::Schneider:
return rhsSchneiderEtAl_(globalPos);
case TestCase::BJSymmetrized:
return rhsBJSymmetrized_(globalPos);
case TestCase::NewICNonSymmetrized:
return rhsNewICNonSymmetrized_(globalPos);
default:
DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
}
......@@ -235,6 +303,10 @@ public:
return analyticalSolutionRybak_(globalPos);
case TestCase::Schneider:
return analyticalSolutionSchneiderEtAl_(globalPos);
case TestCase::BJSymmetrized:
return analyticalSolutionBJSymmetrized_(globalPos);
case TestCase::NewICNonSymmetrized:
return analyticalSolutionNewICNonSymmetrized_(globalPos);
default:
DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
}
......@@ -357,6 +429,54 @@ private:
return NumEqVector(result);
}
// exact solution for BJ-IC with symmetrized stress tensor (by Elissa Eggenweiler)
Dune::FieldVector<Scalar, 3> analyticalSolutionBJSymmetrized_(const GlobalPosition& globalPos) const
{
Dune::FieldVector<Scalar, 3> sol(0.0);
const Scalar x = globalPos[0];
const Scalar y = globalPos[1];
using std::exp; using std::sin; using std::cos;
sol[velocityXIdx] = 2.0*x*x + 2.0*x*y - 6.0*x - 2.0*y*y + 3.0*y - 5.0;
sol[velocityYIdx] = 1.0*x*x + 4.0*x*y - 9.0*x - 1.0*y*y - 1.0;
sol[pressureIdx] = -x*(x-1.0)*(y-1.0) + 1.0/3.0*(y-1.0)*(y-1.0)*(y-1.0) + 4.0 + 2.0*y + x*x;
return sol;
}
// exact solution for BJ-IC with symmetrized stress tensor (by Elissa Eggenweiler)
NumEqVector rhsBJSymmetrized_(const GlobalPosition& globalPos) const
{
const Scalar x = globalPos[0];
using std::sin;
return NumEqVector(8.0*x - 6.0);
}
// exact solution for new IC with non-symmetrized stress tensor (by Elissa Eggenweiler)
Dune::FieldVector<Scalar, 3> analyticalSolutionNewICNonSymmetrized_(const GlobalPosition& globalPos) const
{
Dune::FieldVector<Scalar, 3> sol(0.0);
const Scalar x = globalPos[0];
const Scalar y = globalPos[1];
using std::exp; using std::sin; using std::cos;
sol[velocityXIdx] = x*(y-1.0) + (x-1.0)*(y-1.0) - 3.0;
sol[velocityYIdx] = x*(x-1.0) - (y-1.0)*(y-1.0) - 2.0;
sol[pressureIdx] = x*(1.0-x)*(y-1.0) + 1.0/3.0*(y-1.0)*(y-1.0)*(y-1.0) + 3.0*x + 2.0*y + 1.0;
return sol;
}
// exact solution for new IC with non-symmetrized stress tensor (by Elissa Eggenweiler)
NumEqVector rhsNewICNonSymmetrized_(const GlobalPosition& globalPos) const
{
return NumEqVector(0.0);
}
bool onLeftBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
bool onRightBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
static constexpr Scalar eps_ = 1e-7;
std::shared_ptr<CouplingManager> couplingManager_;
std::string problemName_;
......
......@@ -119,6 +119,10 @@ public:
return rhsRybak_(globalPos);
case TestCase::Schneider:
return rhsSchneiderEtAl_(globalPos);
case TestCase::BJSymmetrized:
return rhsBJSymmetrized_(globalPos);
case TestCase::NewICNonSymmetrized:
return rhsNewICNonSymmetrized_(globalPos);
default:
DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
}
......@@ -211,6 +215,15 @@ public:
return PrimaryVariables(0.0);
}
/*!
* \brief Returns the velocity in the porous medium (which is 0 by default according to Saffmann).
*/
auto porousMediumTerm(const Element& element, const SubControlVolumeFace& scvf) const
{
return couplingManager().couplingData().porousMediumVelocity(element, scvf);
}
/*!
* \brief Returns the intrinsic permeability of required as input parameter
for the Beavers-Joseph-Saffman boundary condition
......@@ -229,6 +242,42 @@ public:
return couplingManager().problem(CouplingManager::porousMediumIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
}
/*!
* \brief Returns the scale separation parameter epsilon required as input parameter for the
new coupling conditions
*/
Scalar epsInterface(const SubControlVolumeFace& scvf) const
{
return couplingManager().problem(CouplingManager::porousMediumIdx).spatialParams().epsInterfaceAtPos(scvf.center());
}
/*!
* \brief Returns the boundary layer constant N_1_bl required as input parameter for the
new coupling condition for the tangential component
*/
Scalar factorNTangential(const SubControlVolumeFace& scvf) const
{
return couplingManager().problem(CouplingManager::porousMediumIdx).spatialParams().factorNTangentialAtPos(scvf.center());
}
/*!
* \brief Returns the boundary layer constant N_s_bl required as input parameter for the
new coupling condition for the tangential component
*/
Scalar factorNMomentum(const SubControlVolumeFace& scvf) const
{
return couplingManager().problem(CouplingManager::porousMediumIdx).spatialParams().factorNMomentumAtPos(scvf.center());
}
/*!
* \brief Returns the boundary layer matrix M_bl required as input parameter for the
new coupling condition for the tangential component
*/
auto matrixNTangential(const SubControlVolumeFace& scvf) const
{
return couplingManager().problem(CouplingManager::porousMediumIdx).spatialParams().matrixNTangentialAtPos(scvf.center());
}
/*!
* \brief Returns the analytical solution of the problem at a given position.
*
......@@ -247,6 +296,10 @@ public:
return analyticalSolutionRybak_(globalPos);
case TestCase::Schneider:
return analyticalSolutionSchneiderEtAl_(globalPos);
case TestCase::BJSymmetrized:
return analyticalSolutionBJSymmetrized_(globalPos);
case TestCase::NewICNonSymmetrized:
return analyticalSolutionNewICNonSymmetrized_(globalPos);
default:
DUNE_THROW(Dune::InvalidStateException, "Invalid test case");
}
......@@ -361,6 +414,55 @@ private:
return source;
}
// exact solution for BJ-IC with symmetrized stress tensor (by Elissa Eggenweiler)
PrimaryVariables analyticalSolutionBJSymmetrized_(const GlobalPosition& globalPos) const
{
PrimaryVariables sol(0.0);
const Scalar x = globalPos[0];
const Scalar y = globalPos[1];
using std::sin; using std::cos;
sol[Indices::velocityXIdx] = (y-1.0)*(y-1.0) + x*(y-1.0) - x - 7.0 + 2.0*x*x + 2.0*y*y;
sol[Indices::velocityYIdx] = (x-1.0)*x - 0.5*(y-1.0)*(y-1.0) - 3.0*y - 3.0 -4.0*y*(x-1.0);
sol[Indices::pressureIdx] = -8.0*x + 1.0*y + 7.0 + x*x;
return sol;
}
// exact solution for BJ-IC with symmetrized stress tensor (by Elissa Eggenweiler)
NumEqVector rhsBJSymmetrized_(const GlobalPosition& globalPos) const
{
const Scalar x = globalPos[0];
NumEqVector source(0.0);
using std::sin; using std::cos;
source[Indices::momentumXBalanceIdx] = 2.0*x - 18.0;
source[Indices::momentumYBalanceIdx] = 0.0;
return source;
}
// exact solution for new IC with non-symmetrized stress tensor (by Elissa Eggenweiler)
PrimaryVariables analyticalSolutionNewICNonSymmetrized_(const GlobalPosition& globalPos) const
{
PrimaryVariables sol(0.0);
const Scalar x = globalPos[0];
const Scalar y = globalPos[1];
using std::exp; using std::sin; using std::cos;
sol[Indices::velocityXIdx] = (y-1.0)*(y-1.0) + x*(y-1.0) + 3.0*x + 1.5;
sol[Indices::velocityYIdx] = 0.5*x*(x-1.0) - 0.5*(y-1.0)*(y-1.0) - 3.0*y + 2.0;
sol[Indices::pressureIdx] = 2.0*x + y - 4.0;
return sol;
}
// exact solution for new IC with non-symmetrized stress tensor (by Elissa Eggenweiler)
NumEqVector rhsNewICNonSymmetrized_(const GlobalPosition& globalPos) const
{
using std::exp; using std::sin; using std::cos;
NumEqVector source(0.0);
source[Indices::momentumXBalanceIdx] = -2.0;
source[Indices::momentumYBalanceIdx] = 1.0;
return source;
}
static constexpr Scalar eps_ = 1e-7;
std::string problemName_;
std::shared_ptr<CouplingManager> couplingManager_;
......
......@@ -62,8 +62,24 @@ public:
ConvergenceTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry, const TestCase testCase)
: ParentType(gridGeometry)
, testCase_(testCase)
, K_(0.0)
{
alphaBJ_ = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
porosity_ = getParam<Scalar>("Darcy.SpatialParams.Porosity", 0.4);
if(testCase_ == TestCase::BJSymmetrized || testCase_ == TestCase::NewICNonSymmetrized)
{
std::vector<Scalar> permeability = getParam<std::vector<Scalar>>("Darcy.SpatialParams.Permeability");
K_[0][0] = permeability[0];
K_[1][1] = permeability[1];
K_[0][1] = K_[1][0] = permeability[2];
}
else
{
Scalar permeability = getParam<Scalar>("Darcy.SpatialParams.Permeability");
K_[0][0] = permeability;
K_[1][1] = permeability;
}
}
/*!
......@@ -79,7 +95,7 @@ public:
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
PermeabilityType K(0.0);
PermeabilityType K(K_);
if (testCase_ == TestCase::Schneider)
{
......@@ -96,12 +112,6 @@ public:
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");
K[0][0] = permeability;
K[1][1] = permeability;
}
return K;
}
......@@ -111,7 +121,7 @@ public:
* \param globalPos The global position
*/
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{ return 0.4; }
{ return porosity_; }
/*! \brief Defines the Beavers-Joseph coefficient in [-].
*
......@@ -120,10 +130,39 @@ public:
Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
{ return alphaBJ_; }
Scalar epsInterfaceAtPos(const GlobalPosition& globalPos) const
{
static const Scalar epsInterface = getParam<Scalar>("Darcy.InterfaceParams.EpsInterface");
return epsInterface;
}
Scalar factorNMomentumAtPos(const GlobalPosition& globalPos) const
{
static const Scalar N_s_bl = getParam<Scalar>("Darcy.InterfaceParams.N_s_bl");
return N_s_bl;
}
Scalar factorNTangentialAtPos(const GlobalPosition& globalPos) const
{
static const Scalar N_1_bl = getParam<Scalar>("Darcy.InterfaceParams.N_1_bl");
return N_1_bl;
}
PermeabilityType matrixNTangentialAtPos(const GlobalPosition& globalPos) const
{
static const std::vector<Scalar> M_bl = getParam<std::vector<Scalar>>("Darcy.InterfaceParams.M_bl");
PermeabilityType M(0.0);
M[0][0]= M_bl[0];
M[1][1]= M_bl[1];
return M;
}
private:
TestCase testCase_;
Scalar permeability_;
PermeabilityType K_;
Scalar alphaBJ_;
Scalar porosity_;
};
} // end namespace Dumux
......
......@@ -29,7 +29,7 @@ namespace Dumux {
enum class TestCase
{
ShiueExampleOne, ShiueExampleTwo, Rybak, Schneider
ShiueExampleOne, ShiueExampleTwo, Rybak, Schneider, BJSymmetrized, NewICNonSymmetrized
};
} // end namespace Dumux
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment