Commit 687e7ea9 authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

[poromechanics] Add elastic one-phase test using multidomain

parent 2cf55624
add_subdirectory(poromechanics)
// -*- 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 MultiDomain
* \ingroup OnePTests
* \ingroup PoroElastic
* \brief Definition of the spatial parameters for the single-phase flow
* sub-problem in the coupled poro-mechanical el1p problem.
*/
#ifndef DUMUX_1P_SUB_PROBLEM_HH
#define DUMUX_1P_SUB_PROBLEM_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 <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/components/constant.hh>
#include "1pspatialparams.hh"
namespace Dumux {
// forward declaration of the problem class
template <class TypeTag>
class OnePSubProblem;
namespace Properties {
NEW_TYPE_TAG(OnePSubTypeTag, INHERITS_FROM(CCTpfaModel, OneP));
// The fluid phase consists of one constant component
SET_TYPE_PROP(OnePSubTypeTag,
FluidSystem,
Dumux::FluidSystems::OnePLiquid< typename GET_PROP_TYPE(TypeTag, Scalar),
Dumux::Components::Constant<0, typename GET_PROP_TYPE(TypeTag, Scalar)> >);
// Set the grid type
SET_TYPE_PROP(OnePSubTypeTag, Grid, Dune::YaspGrid<2>);
// Set the problem property
SET_TYPE_PROP(OnePSubTypeTag, Problem, OnePSubProblem<TypeTag> );
// Set the spatial parameters
SET_TYPE_PROP(OnePSubTypeTag, SpatialParams, OnePSpatialParams<TypeTag> );
} // end namespace Properties
/*!
* \ingroup MultiDomain
* \ingroup OnePTests
* \ingroup PoroElastic
*
* \brief The single-phase sub problem in the el1p coupled problem.
*/
template <class TypeTag>
class OnePSubProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
// copy pressure index for convenience
enum { pressureIdx = GET_PROP_TYPE(TypeTag, ModelTraits)::Indices::pressureIdx };
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
public:
OnePSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, const std::string& paramGroup = "")
: ParentType(fvGridGeometry, paramGroup)
{}
//! Return the temperature within the domain in [K].
Scalar temperature() const
{ return 273.15 + 10; } // 10C
//! Evaluate the boundary conditions for a Dirichlet boundary segment.
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{ return initialAtPos(globalPos); }
//! Evaluate the initial value for a control volume.
PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{ return PrimaryVariables(1.0e5); }
//! Evaluate source terms
NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
{
static const Scalar source = getParam<Scalar>("Problem.InjectionRate");
if (globalPos[0] > 0.4 && globalPos[0] < 0.6 && globalPos[1] < 0.6 && globalPos[1] > 0.4)
return NumEqVector(source);
return NumEqVector(0.0);
}
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
* \param globalPos The global position
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
values.setAllDirichlet();
return values;
}
private:
static constexpr Scalar eps_ = 1.0e-6;
};
} //end namespace Dumux
#endif
// -*- 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 MultiDomain
* \ingroup OnePTests
* \ingroup PoroElastic
* \brief The spatial parameters class for the test problem using the
* 1p box model
*/
#ifndef DUMUX_1P_TEST_SPATIALPARAMS_HH
#define DUMUX_1P_TEST_SPATIALPARAMS_HH
#include <dumux/discretization/elementsolution.hh>
#include <dumux/material/spatialparams/fv1p.hh>
#include <dumux/material/spatialparams/gstatrandomfield.hh>
#include <dumux/material/fluidmatrixinteractions/porositydeformation.hh>
namespace Dumux {
/*!
* \ingroup MultiDomain
* \ingroup OnePTests
* \ingroup PoroElastic
* \brief The spatial parameters class for the test problem using the
* 1p box model
*/
template<class TypeTag>
class OnePSpatialParams : public FVSpatialParamsOneP< typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
typename GET_PROP_TYPE(TypeTag, Scalar),
OnePSpatialParams<TypeTag> >
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using SubControlVolume = typename FVGridGeometry::SubControlVolume;
using GridView = typename FVGridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using ThisType = OnePSpatialParams<TypeTag>;
using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, ThisType>;
public:
// export permeability type
using PermeabilityType = Scalar;
OnePSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
, permeability_(getParam<Scalar>("SpatialParams.Permeability"))
, initPorosity_(getParam<Scalar>("SpatialParams.InitialPorosity"))
{}
//! Return the permeability at a given position
PermeabilityType permeabilityAtPos(const GlobalPosition& globalPoss) const
{ return permeability_; }
//! Return the porosity for a sub-control volume
template<class ElementSolution>
Scalar porosity(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
static constexpr auto poroMechId = CouplingManager::poroMechId;
const auto& poroMechGridGeom = couplingManagerPtr_->template problem<poroMechId>().fvGridGeometry();
const auto poroMechElemSol = elementSolution(element, couplingManagerPtr_->curSol()[poroMechId], poroMechGridGeom);
// evaluate the deformation-dependent porosity at the scv center
return PorosityDeformation<Scalar>::evaluatePorosity(poroMechGridGeom, element, scv.center(), poroMechElemSol, initPorosity_);
}
//! sets the pointer to the coupling manager.
void setCouplingManager(std::shared_ptr<const CouplingManager> cm)
{ couplingManagerPtr_ = cm; }
//! returns reference to the coupling manager.
const CouplingManager& couplingManager() const
{ return *couplingManagerPtr_; }
private:
std::shared_ptr<const CouplingManager> couplingManagerPtr_;
Scalar permeability_;
Scalar initPorosity_;
};
} // end namespace Dumux
#endif
dune_symlink_to_source_files(FILES "el1p.input")
dune_add_test(NAME test_el1p
SOURCES test_el1p.cc
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/el1p_p_reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/onep.vtu
${CMAKE_SOURCE_DIR}/test/references/el1p_u_reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/poromech.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_el1p el1p.input")
set(CMAKE_BUILD_TYPE Release)
#install sources
install(FILES
1pproblem.hh
1pspatialparams.hh
poroelasticproblem.hh
poroelasticspatialparams.hh
test_el1p.cc
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/multidomain/poromechanics/el1p)
[Grid]
LowerLeft = 0 0
UpperRight = 1 1
Cells = 10 10
[Problem]
InjectionRate = 1e-3
EnableGravity = false
[PoroElastic.Problem]
Name = poroelastic # name passed to the output routines
[OneP.Problem]
Name = onep # name passed to the output routines
[SpatialParams]
Permeability = 1e-12 # [m^2]
InitialPorosity = 0.2 # [-]
[LinearSolver]
ResidualReduction = 1e-20
MaxIterations = 2000
[Component]
SolidDensity = 2700
// -*- 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 MultiDomain
* \ingroup Geomechanics
* \ingroup PoroElastic
* \brief The poro-elastic sub-problem in the el1p coupled problem.
*/
#ifndef DUMUX_POROELASTIC_SUBPROBLEM_HH
#define DUMUX_POROELASTIC_SUBPROBLEM_HH
#include <dune/common/fmatrix.hh>
#include <dumux/discretization/box/properties.hh>
#include <dumux/geomechanics/poroelastic/model.hh>
#include <dumux/geomechanics/fvproblem.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/components/constant.hh>
#include "poroelasticspatialparams.hh"
namespace Dumux {
// forward declaration of the problem class
template <class TypeTag>
class PoroElasticSubProblem;
namespace Properties {
NEW_TYPE_TAG(PoroElasticSubTypeTag, INHERITS_FROM(BoxModel, PoroElastic));
// Set the grid type
SET_TYPE_PROP(PoroElasticSubTypeTag, Grid, Dune::YaspGrid<2>);
// Set the problem property
SET_TYPE_PROP(PoroElasticSubTypeTag, Problem, Dumux::PoroElasticSubProblem<TypeTag>);
// The fluid phase consists of one constant component
SET_TYPE_PROP(PoroElasticSubTypeTag,
FluidSystem,
Dumux::FluidSystems::OnePLiquid< typename GET_PROP_TYPE(TypeTag, Scalar),
Dumux::Components::Constant<0, typename GET_PROP_TYPE(TypeTag, Scalar)> >);
// The spatial parameters property
SET_TYPE_PROP(PoroElasticSubTypeTag, SpatialParams, PoroElasticSpatialParams< typename GET_PROP_TYPE(TypeTag, Scalar),
typename GET_PROP_TYPE(TypeTag, FVGridGeometry) >);
} // end namespace Properties
/*!
* \ingroup MultiDomain
* \ingroup Geomechanics
* \ingroup PoroElastic
*
* \brief The poro-elastic sub-problem in the el1p coupled problem.
*/
template<class TypeTag>
class PoroElasticSubProblem : public GeomechanicsFVProblem<TypeTag>
{
using ParentType = GeomechanicsFVProblem<TypeTag>;
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 CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVGridGeometry::SubControlVolume;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
using GradU = Dune::FieldMatrix<Scalar, dim, dimWorld>;
public:
//! The constructor
PoroElasticSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, const std::string& paramGroup = "")
: ParentType(fvGridGeometry, paramGroup) {}
//! The temperature in the domain
static constexpr Scalar temperature()
{ return 273.15; }
//! Evaluate the initial value for a control volume.
PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{ return PrimaryVariables(0.0); }
//! Evaluate the boundary conditions for a Dirichlet boundary segment.
PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{ return PrimaryVariables(0.0); }
//! Evaluate the boundary conditions for a Neumannboundary segment.
PrimaryVariables neumannAtPos(const GlobalPosition& globalPos) const
{ return PrimaryVariables(0.0); }
/*!
* \brief Returns the effective fluid density
*/
Scalar effectiveFluidDensity(const Element& element,
const SubControlVolume& scv) const
{
// get context from coupling manager
// here, we know that the flow problem uses cell-centered finite volumes,
// thus, we simply take the volume variables of the element and return the density
const auto& context = couplingManager().poroMechanicsCouplingContext();
return (*context.pmFlowElemVolVars)[scv.elementIndex()].density();
}
/*!
* \brief Returns the effective pore pressure
*/
template< class FluxVarsCache >
Scalar effectivePorePressure(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const FluxVarsCache& fluxVarsCache) const
{
// get context from coupling manager
// here, we know that the flow problem uses cell-centered finite volumes,
// thus, we simply take the volume variables of the element and return the pressure
const auto& context = couplingManager().poroMechanicsCouplingContext();
const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
return (*context.pmFlowElemVolVars)[eIdx].pressure();
}
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
* \param globalPos The global position
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
{
BoundaryTypes values;
values.setAllDirichlet();
return values;
}
/*!
* \brief Evaluate the source term for all phases within a given
* sub-control-volume.
*/
PrimaryVariables source(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume& scv) const
{ return PrimaryVariables(0.0); }
//! sets the pointer to the coupling manager.
void setCouplingManager(std::shared_ptr<const CouplingManager> cm)
{ couplingManagerPtr_ = cm; }
//! returns reference to the coupling manager.
const CouplingManager& couplingManager() const
{ return *couplingManagerPtr_; }
private:
std::shared_ptr<const CouplingManager> couplingManagerPtr_;
static constexpr Scalar eps_ = 3e-6;
};
} //end namespace
#endif
// -*- 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 MultiDomain
* \ingroup Geomechanics
* \ingroup PoroElastic
* \brief Definition of the spatial parameters for the poro-elastic problem
*/
#ifndef DUMUX_POROELASTIC_SPATIAL_PARAMS_HH
#define DUMUX_POROELASTIC_SPATIAL_PARAMS_HH
#include <dumux/geomechanics/lameparams.hh>
#include <dumux/material/spatialparams/fvporoelastic.hh>
#include <dumux/material/fluidmatrixinteractions/porositydeformation.hh>
namespace Dumux {
/*!
* \ingroup MultiDomain
* \ingroup Geomechanics
* \ingroup PoroElastic
* \brief Definition of the spatial parameters for the poro-elastic
* sub-problem in the coupled poro-mechanical el1p problem.
*/
template<class Scalar, class FVGridGeometry>
class PoroElasticSpatialParams : public FVSpatialParamsPoroElastic< Scalar,
FVGridGeometry,
PoroElasticSpatialParams<Scalar, FVGridGeometry> >
{
using ThisType = PoroElasticSpatialParams<Scalar, FVGridGeometry>;
using ParentType = FVSpatialParamsPoroElastic<Scalar, FVGridGeometry, ThisType>;
using SubControlVolume = typename FVGridGeometry::SubControlVolume;
using GridView = typename FVGridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
//! export the type of the lame parameters
using LameParams = Dumux::LameParams<Scalar>;
//! The constructor
PoroElasticSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
, initPorosity_(getParam<Scalar>("SpatialParams.InitialPorosity"))
{
// Young's modulus [Pa]
Scalar E = 6.e9;
// Poisson's ratio [-]
Scalar nu = 0.2;
// Lame parameters [Pa]
lameParams_.setLambda( (E * nu) / ((1 + nu)*(1 - 2 * nu)) );
lameParams_.setMu( E / (2 * (1 + nu)) );
}
//! Define the Lame parameters
const LameParams& lameParamsAtPos(const GlobalPosition& globalPos) const
{ return lameParams_; }
//! Return the porosity of the porous medium
template<class ElementSolution>
Scalar porosity(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
return PorosityDeformation<Scalar>::evaluatePorosity(this->fvGridGeometry(), element, scv, elemSol, initPorosity_);
}
//! Return the biot coefficient of the porous medium
Scalar biotCoefficientAtPos(const GlobalPosition& globalPos) const
{ return 1.0; }
private:
Scalar initPorosity_;
LameParams lameParams_;
};
} // end namespace Dumux
#endif
// -*- 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 *