Commit 1417d7c9 authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[stokesdarcy][convergencetest] Add Navier-Stokes convergence test

* based on Schneider et al. 2019:
  Coupling staggered-grid and MPFA finite volume methods for
  free flow/porous-medium flow problems
parent d688ea14
......@@ -11,3 +11,14 @@ dune_add_test(NAME test_md_boundary_darcy1p_stokes1p_convtest
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_freeflow1p_convtest params.input
-Problem.TestCase Schneider
-FreeFlow.Problem.EnableInertiaTerms true
-FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph true)
......@@ -292,16 +292,16 @@ int main(int argc, char** argv) try
const auto testCase = []()
{
const auto testCaseInput = getParam<std::string>("Problem.TestCase", "ShiueExampleTwo");
TestCase tc;
if (testCaseInput == "ShiueExampleOne")
tc = TestCase::ShiueExampleOne;
return TestCase::ShiueExampleOne;
else if (testCaseInput == "ShiueExampleTwo")
tc = TestCase::ShiueExampleTwo;
return TestCase::ShiueExampleTwo;
else if (testCaseInput == "Rybak")
tc = TestCase::Rybak;
return TestCase::Rybak;
else if (testCaseInput == "Schneider")
return TestCase::Schneider;
else
DUNE_THROW(Dune::InvalidStateException, testCaseInput + " is not a valid test case");
return tc;
}();
using FreeFlowProblem = GetPropType<FreeFlowTypeTag, Properties::Problem>;
......
......@@ -32,3 +32,6 @@ LiquidKinematicViscosity = 1.0
[Assembly]
NumericDifference.BaseEpsilon = 1e-4
[FreeFlow.Flux]
UpwindWeight = 0.5
......@@ -33,7 +33,7 @@
#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>
......@@ -70,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
......@@ -223,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");
}
......@@ -258,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");
}
......@@ -336,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_;
......
......@@ -140,6 +140,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");
}
......@@ -265,6 +267,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");
}
......@@ -343,6 +347,42 @@ 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"
PrimaryVariables analyticalSolutionSchneiderEtAl_(const GlobalPosition& globalPos) const
{
PrimaryVariables sol(0.0);
const Scalar x = globalPos[0];
const Scalar y = globalPos[1];
using std::sin;
static constexpr Scalar omega = M_PI;
const Scalar sinOmegaX = sin(omega*x);
sol[Indices::velocityXIdx] = y;
sol[Indices::velocityYIdx] = -y*sinOmegaX;
sol[Indices::pressureIdx] = -y*y*sinOmegaX*sinOmegaX;
return sol;
}
// see Schneider et al., 2019: "Coupling staggered-grid and MPFA finite volume methods for
// free flow/porous-medium flow problems"
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;
const Scalar sinOmegaX = sin(omega*x);
const Scalar cosOmegaX = cos(omega*x);
NumEqVector source(0.0);
source[Indices::conti0EqIdx] = -sinOmegaX;
source[Indices::momentumXBalanceIdx] = -2*omega*y*y*sinOmegaX*cosOmegaX
-2*y*sinOmegaX + omega*cosOmegaX;
source[Indices::momentumYBalanceIdx] = -omega*y*y*cosOmegaX - omega*omega*y*sinOmegaX;
return source;
}
static constexpr Scalar eps_ = 1e-7;
std::string problemName_;
std::shared_ptr<CouplingManager> couplingManager_;
......
// -*- 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_CONVERGENCE_TEST_SPATIALPARAMS_HH
#define DUMUX_CONVERGENCE_TEST_SPATIALPARAMS_HH
#include <cmath>
#include <dumux/common/math.hh>
#include <dumux/common/parameters.hh>
#include <dune/common/fmatrix.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 TestCase 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_ == TestCase::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");
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:
TestCase testCase_;
Scalar permeability_;
Scalar alphaBJ_;
};
} // end namespace Dumux
#endif // DUMUX_CONVERGENCE_TEST_SPATIALPARAMS_HH
......@@ -29,7 +29,7 @@ namespace Dumux {
enum class TestCase
{
ShiueExampleOne, ShiueExampleTwo, Rybak
ShiueExampleOne, ShiueExampleTwo, Rybak, Schneider
};
} // 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