Commit 385b464f authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Timo Koch
Browse files

[test][multidomain][stokesdarcy] Add 1p2c 1p2c stokes darcy test

parent ff21c437
// -*- 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>("SpatialParams.Permeability");
porosity_ = getParam<Scalar>("SpatialParams.Porosity");
alphaBJ_ = getParam<Scalar>("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 porosity_; }
/*! \brief Define the Beavers-Joseph coefficient in [-].
*
* \param globalPos The global position
*/
Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
{ return alphaBJ_; }
private:
Scalar permeability_;
Scalar porosity_;
Scalar alphaBJ_;
};
} // end namespace
#endif
add_subdirectory("verticalflow")
add_subdirectory("horizontalflow")
add_input_file_links()
add_executable(test_stokes1p2cdarcy1p2chorizontal EXCLUDE_FROM_ALL test_stokes1p2cdarcy1p2chorizontal.cc)
dune_add_test(NAME test_stokes1p2cdarcy1p2chorizontal
TARGET test_stokes1p2cdarcy1p2chorizontal
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_stokes1p2cdarcy1p2chorizontal_stokes-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_stokes1p2cdarcy1p2chorizontal_stokes-00020.vtu
${CMAKE_SOURCE_DIR}/test/references/test_stokes1p2cdarcy1p2chorizontal_darcy-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_stokes1p2cdarcy1p2chorizontal_darcy-00020.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_stokes1p2cdarcy1p2chorizontal test_stokes1p2cdarcy1p2chorizontal.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/1pnc/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include "./../1pspatialparams.hh"
#include <dumux/material/fluidsystems/h2oair.hh>
#include <dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh>
namespace Dumux
{
template <class TypeTag>
class DarcySubProblem;
namespace Properties
{
NEW_TYPE_TAG(DarcyOnePTwoCTypeTag, INHERITS_FROM(CCTpfaModel, OnePNC));
// Set the problem property
SET_TYPE_PROP(DarcyOnePTwoCTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
// The fluid system
SET_TYPE_PROP(DarcyOnePTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
// Use water as phase
SET_INT_PROP(DarcyOnePTwoCTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::liquidPhaseIdx);
// Use moles
SET_BOOL_PROP(DarcyOnePTwoCTypeTag, UseMoles, true);
// Do not replace one equation with a total mass balance
SET_INT_PROP(DarcyOnePTwoCTypeTag, ReplaceCompEqIdx, 3);
//! Use a model with constant tortuosity for the effective diffusivity
SET_TYPE_PROP(DarcyOnePTwoCTypeTag, EffectiveDiffusivityModel,
DiffusivityConstantTortuosity<typename GET_PROP_TYPE(TypeTag, Scalar)>);
// Set the grid type
SET_TYPE_PROP(DarcyOnePTwoCTypeTag, Grid, Dune::YaspGrid<2>);
// Set the spatial paramaters type
SET_TYPE_PROP(DarcyOnePTwoCTypeTag, 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 FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
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);
// copy some indices for convenience
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
enum {
// grid and world dimension
dim = GridView::dimension,
dimworld = GridView::dimensionworld,
// primary variable indices
conti0EqIdx = Indices::conti0EqIdx,
pressureIdx = Indices::pressureIdx,
};
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
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)
{
pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
initialMoleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.InitialMoleFraction");
}
/*!
* \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 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 = 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
{
PrimaryVariables values(0.0);
values[pressureIdx] = pressure_;
values[conti0EqIdx + 1] = initialMoleFraction_;
return values;
}
// \}
//! 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_;
Scalar pressure_;
Scalar initialMoleFraction_;
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 Navier-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/h2oair.hh>
#include <dumux/freeflow/navierstokes/problem.hh>
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
namespace Dumux
{
template <class TypeTag>
class StokesSubProblem;
namespace Properties
{
NEW_TYPE_TAG(StokesOnePTwoCTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
// Set the fluid system
SET_TYPE_PROP(StokesOnePTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
// Set the grid type
SET_TYPE_PROP(StokesOnePTwoCTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
// Set the problem property
SET_TYPE_PROP(StokesOnePTwoCTypeTag, Problem, Dumux::StokesSubProblem<TypeTag> );
SET_BOOL_PROP(StokesOnePTwoCTypeTag, EnableFVGridGeometryCache, true);
SET_BOOL_PROP(StokesOnePTwoCTypeTag, EnableGridFluxVariablesCache, true);
SET_BOOL_PROP(StokesOnePTwoCTypeTag, EnableGridVolumeVariablesCache, true);
// Do consider intertia terms
SET_BOOL_PROP(StokesOnePTwoCTypeTag, EnableInertiaTerms, true);
// Use moles
SET_BOOL_PROP(StokesOnePTwoCTypeTag, UseMoles, true);
// Do not replace one equation with a total mass balance
SET_INT_PROP(StokesOnePTwoCTypeTag, ReplaceCompEqIdx, 3);
}
/*!
* \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 FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
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 PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
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), injectionState_(false), couplingManager_(couplingManager)
{
inletVelocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Velocity");
pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
inletMoleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.InletMoleFraction");
}
/*!
* \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))
{
values.setDirichlet(Indices::conti0EqIdx + 1);
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
}
else if(onRightBoundary_(globalPos))
{
values.setDirichlet(Indices::pressureIdx);
values.setOutflow(Indices::conti0EqIdx + 1);
}
else
{
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
values.setNeumann(Indices::conti0EqIdx);
values.setNeumann(Indices::conti0EqIdx + 1);
}
if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
{
values.setCouplingNeumann(Indices::conti0EqIdx);
values.setCouplingNeumann(Indices::conti0EqIdx + 1);
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);
// start injecting after the velocity field had enough time to initialize
if(globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_ && isInjectionPeriod())
values[Indices::conti0EqIdx + 1] = inletMoleFraction_;
return values;
}
/*!