Commit fa239dc2 authored by Sina Ackermann's avatar Sina Ackermann Committed by Timo Koch
Browse files

[test][multidomain][stokesdarcy] Add 1p/1p stokesdarcy test problems

* vertical flow
* horizontal flow
parent 41bc6c93
add_subdirectory(darcydarcy)
add_subdirectory(stokesdarcy)
// -*- 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 2 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 OnePTests
* \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>
namespace Dumux
{
/*!
* \ingroup OnePModel
* \ingroup ImplicitTestProblems
*
* \brief The spatial parameters class for the test problem using the
* 1p cc model
*/
template<class TypeTag>
class OnePSpatialParams
: public FVSpatialParamsOneP<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
typename GET_PROP_TYPE(TypeTag, Scalar),
OnePSpatialParams<TypeTag>>
{
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePSpatialParams<TypeTag>>;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
// export permeability type
using PermeabilityType = Scalar;
OnePSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
permeability_ = getParam<Scalar>("Darcy.SpatialParams.Permeability");
alphaBJ_ = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
}
/*!
* \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
*
* \param globalPos The global position
* \return the intrinsic permeability
*/
PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
{ return permeability_; }
/*! \brief Define the porosity in [-].
*
* \param globalPos The global position
*/
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{ return 0.4; }
/*! \brief Define the Beavers-Joseph coefficient in [-].
*
* \param globalPos The global position
*/
Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
{ return alphaBJ_; }
private:
Scalar permeability_;
Scalar alphaBJ_;
};
} // end namespace
#endif
add_subdirectory("verticalflow")
add_subdirectory("horizontalflow")
add_input_file_links()
dune_add_test(NAME test_stokes1pdarcy1phorizontal
SOURCES test_stokes1pdarcy1phorizontal.cc
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_stokes1pdarcy1phorizontal_stokes-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_stokes1pdarcy1phorizontal_stokes-00001.vtu
${CMAKE_SOURCE_DIR}/test/references/test_stokes1pdarcy1phorizontal_darcy-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_stokes1pdarcy1phorizontal_darcy-00001.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_stokes1pdarcy1phorizontal test_stokes1pdarcy1phorizontal.input")
// -*- 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 2 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
*
* \brief A simple Darcy test problem (cell-centered finite volume method).
*/
#ifndef DUMUX_DARCY_SUBPROBLEM_HH
#define DUMUX_DARCY_SUBPROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
#include <dumux/porousmediumflow/1p/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include "./../1pspatialparams.hh"
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
namespace Dumux
{
template <class TypeTag>
class DarcySubProblem;
namespace Properties
{
NEW_TYPE_TAG(DarcyOnePTypeTag, INHERITS_FROM(CCTpfaModel, OneP));
// Set the problem property
SET_TYPE_PROP(DarcyOnePTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
// the fluid system
SET_PROP(DarcyOnePTypeTag, FluidSystem)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
};
// Set the grid type
SET_TYPE_PROP(DarcyOnePTypeTag, Grid, Dune::YaspGrid<2>);
SET_TYPE_PROP(DarcyOnePTypeTag, SpatialParams, OnePSpatialParams<TypeTag>);
}
template <class TypeTag>
class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
public:
DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
std::shared_ptr<CouplingManager> couplingManager)
: ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
{}
/*!
* \name Simulation steering
*/
// \{
/*!
* \brief Returns true if a restart file should be written to
* disk.
*/
bool shouldWriteRestartFile() const
{ return false; }
/*!
* \name Problem parameters
*/
// \{
bool shouldWriteOutput() const // define output
{ return true; }
/*!
* \brief Return the temperature within the domain in [K].
*
*/
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;
values.setAllNeumann();
if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
values.setAllCouplingNeumann();
return values;
}
/*!
* \brief Evaluate 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
*
* For this method, the \a values parameter stores primary variables.
*/
PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
{
PrimaryVariables values(0.0);
values = initial(element);
return values;
}
/*!
* \brief Evaluate the boundary conditions for a Neumann control volume.
*
* \param element The element for which the Neumann boundary condition is set
* \param fvGeomentry The fvGeometry
* \param elemVolVars The element volume variables
* \param scvf The boundary sub control volume face
*
* For this method, the \a values variable stores primary variables.
*/
template<class ElementVolumeVariables>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
return values;
}
// \}
/*!
* \name Volume terms
*/
// \{
/*!
* \brief Evaluate the source term for all phases within a given
* sub-control-volume.
*
* \param element The element for which the source term is set
* \param fvGeomentry The fvGeometry
* \param elemVolVars The element volume variables
* \param scv The subcontrolvolume
*/
template<class ElementVolumeVariables>
NumEqVector source(const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
{ return NumEqVector(0.0); }
// \}
/*!
* \brief Evaluate the initial value for a control volume.
*
* \param element The element
*
* For this method, the \a priVars parameter stores primary
* variables.
*/
PrimaryVariables initial(const Element &element) const
{
return PrimaryVariables(0.0);
}
// \}
//! Set the coupling manager
void setCouplingManager(std::shared_ptr<CouplingManager> cm)
{ couplingManager_ = cm; }
//! Get the coupling manager
const CouplingManager& couplingManager() const
{ return *couplingManager_; }
private:
bool onLeftBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
bool onRightBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
bool onLowerBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
bool onUpperBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
Scalar eps_;
std::shared_ptr<CouplingManager> couplingManager_;
};
} //end namespace
#endif //DUMUX_DARCY_SUBPROBLEM_HH
// -*- 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 2 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 NavierStokesTests
* \brief A simple Stokes test problem for the staggered grid (Navier-)Stokes model.
*/
#ifndef DUMUX_STOKES_SUBPROBLEM_HH
#define DUMUX_STOKES_SUBPROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/freeflow/navierstokes/problem.hh>
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/freeflow/navierstokes/model.hh>
namespace Dumux
{
template <class TypeTag>
class StokesSubProblem;
namespace Properties
{
NEW_TYPE_TAG(StokesOnePTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokes));
// the fluid system
SET_PROP(StokesOnePTypeTag, FluidSystem)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
};
// Set the grid type
SET_TYPE_PROP(StokesOnePTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
// Set the problem property
SET_TYPE_PROP(StokesOnePTypeTag, Problem, Dumux::StokesSubProblem<TypeTag> );
SET_BOOL_PROP(StokesOnePTypeTag, EnableFVGridGeometryCache, true);
SET_BOOL_PROP(StokesOnePTypeTag, EnableGridFluxVariablesCache, true);
SET_BOOL_PROP(StokesOnePTypeTag, EnableGridVolumeVariablesCache, true);
// #if ENABLE_NAVIERSTOKES
// SET_BOOL_PROP(StokesOnePTypeTag, EnableInertiaTerms, true);
// #else
SET_BOOL_PROP(StokesOnePTypeTag, EnableInertiaTerms, false);
// #endif
}
/*!
* \ingroup NavierStokesTests
* \brief Test problem for the one-phase (Navier-) Stokes problem.
*
* Horizontal flow from left to right with a parabolic velocity profile.
*/
template <class TypeTag>
class StokesSubProblem : public NavierStokesProblem<TypeTag>
{
using ParentType = NavierStokesProblem<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
public:
StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
: ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
{
deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference");
}
/*!
* \name Problem parameters
*/
// \{
bool shouldWriteRestartFile() const
{ return false; }
/*!
* \brief Return 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
/*!
* \brief Return the sources within the domain.
*
* \param globalPos The global position
*/
NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
{ return NumEqVector(0.0); }
// \}
/*!
* \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;
const auto& globalPos = scvf.dofPosition();
if(onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
values.setDirichlet(Indices::pressureIdx);
else
{
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
}
if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
{
values.setCouplingNeumann(Indices::conti0EqIdx);
values.setCouplingNeumann(Indices::momentumYBalanceIdx);
values.setBJS(Indices::momentumXBalanceIdx);
}
return values;
}
/*!
* \brief Evaluate the boundary conditions for a Dirichlet control volume.
*
* \param element The element
* \param scvf The sub control volume face
*/
PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{
PrimaryVariables values(0.0);
values = initialAtPos(globalPos);
return values;
}
/*!
* \brief Evaluate the boundary conditions for a Neumann control volume.
*
* \param element The element for which the Neumann boundary condition is set
* \param fvGeomentry 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 ElementFaceVariables>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
{
values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
}
return values;
}
// \}
//! Set the coupling manager
void setCouplingManager(std::shared_ptr<CouplingManager> cm)
{ couplingManager_ = cm; }
//! Get the coupling manager