Commit 0cefaafa authored by Holger Class's avatar Holger Class Committed by Dennis Gläser
Browse files

add test for el2p

parent c4e1476a
add_subdirectory("el1p")
add_subdirectory("el2p")
// -*- 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 TwoPTests
* \ingroup PoroElastic
* \brief Definition of the spatial parameters for the two-phase flow
* sub-problem in the coupled poro-mechanical elp problem.
*/
#ifndef DUMUX_2P_SUB_PROBLEM_HH
#define DUMUX_2P_SUB_PROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
#include <dumux/porousmediumflow/2p/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/fluidsystems/1pgas.hh>
#include <dumux/material/fluidsystems/2pimmiscible.hh>
#include <dumux/material/components/h2o.hh>
#include <dumux/material/components/air.hh>
#include "2pspatialparams.hh"
namespace Dumux {
// forward declaration of the problem class
template <class TypeTag>
class TwoPSubProblem;
namespace Properties {
NEW_TYPE_TAG(TwoPSubTypeTag, INHERITS_FROM(CCTpfaModel, TwoP));
// Set the fluid system for TwoPSubProblem
SET_PROP(TwoPSubTypeTag, FluidSystem)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::H2O<Scalar> >;
using NonwettingPhase = FluidSystems::OnePGas<Scalar, Components::Air<Scalar> >;
using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
};
// Set the grid type
SET_TYPE_PROP(TwoPSubTypeTag, Grid, Dune::YaspGrid<2>);
// Set the problem property
SET_TYPE_PROP(TwoPSubTypeTag, Problem, TwoPSubProblem<TypeTag> );
// Set the spatial parameters
SET_TYPE_PROP(TwoPSubTypeTag, SpatialParams, TwoPSpatialParams<TypeTag> );
} // end namespace Properties
/*!
* \ingroup MultiDomain
* \ingroup TwoPTests
* \ingroup PoroElastic
*
* \brief The two-phase sub problem in the el2p coupled problem.
*/
template <class TypeTag>
class TwoPSubProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
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,
saturationNIdx = GET_PROP_TYPE(TypeTag, ModelTraits)::Indices::saturationIdx,
waterPhaseIdx = FluidSystem::phase0Idx,
gasPhaseIdx = FluidSystem::phase1Idx };
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:
TwoPSubProblem(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
{
PrimaryVariables values;
values[pressureIdx] = 1.0e5;
values[saturationNIdx] = 0.5;
return values;
}
//! Evaluate source terms
NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
{
NumEqVector values(0.0);
static const Scalar sourceG = getParam<Scalar>("Problem.InjectionRateGas");
static const Scalar sourceW = getParam<Scalar>("Problem.InjectionRateWater");
if (globalPos[0] > 0.4 && globalPos[0] < 0.6 && globalPos[1] < 0.6 && globalPos[1] > 0.4)
{
values[gasPhaseIdx] = sourceG;
values[waterPhaseIdx] = sourceW;
}
return values;
}
/*!
* \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 TwoPTests
* \ingroup PoroElastic
* \brief The spatial parameters class for the two-phase sub problem in the el2p test problem
*/
#ifndef DUMUX_2P_TEST_SPATIALPARAMS_HH
#define DUMUX_2P_TEST_SPATIALPARAMS_HH
#include <dumux/discretization/elementsolution.hh>
#include <dumux/material/spatialparams/fv.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include <dumux/material/spatialparams/gstatrandomfield.hh>
#include <dumux/material/fluidmatrixinteractions/porositydeformation.hh>
#include <dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh>
namespace Dumux {
/*!
* \ingroup MultiDomain
* \ingroup TwoPTests
* \ingroup PoroElastic
* \brief The spatial parameters class for the two-phase sub problem in the el2p test problem
*/
template<class TypeTag>
class TwoPSpatialParams : public FVSpatialParams< typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
typename GET_PROP_TYPE(TypeTag, Scalar),
TwoPSpatialParams<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 = TwoPSpatialParams<TypeTag>;
using ParentType = FVSpatialParams<FVGridGeometry, Scalar, ThisType>;
public:
using EffectiveLaw = RegularizedVanGenuchten<Scalar>;
using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
using MaterialLawParams = typename MaterialLaw::Params;
// export permeability type
using PermeabilityType = Scalar;
TwoPSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
, initPermeability_(getParam<Scalar>("SpatialParams.Permeability"))
, initPorosity_(getParam<Scalar>("SpatialParams.InitialPorosity"))
{
myMaterialParams_.setSwr(0.0);
myMaterialParams_.setSnr(0.0);
myMaterialParams_.setVgAlpha(0.0005);
myMaterialParams_.setVgn(4.);
}
//! 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_);
}
//! Function for defining the (intrinsic) permeability \f$[m^2]\f$.
template< class ElementSolution >
PermeabilityType permeability(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
PermeabilityKozenyCarman<PermeabilityType> permLaw;
return permLaw.evaluatePermeability(initPermeability_, initPorosity_, porosity(element, scv, elemSol));
}
/*!
* \brief Returns the parameter object for the Brooks-Corey material law.
* In this test, we use element-wise distributed material parameters.
*
* \param element The current element
* \param scv The sub-control volume inside the element.
* \param elemSol The solution at the dofs connected to the element.
* \return the material parameters object
*/
template<class ElementSolution>
const MaterialLawParams& materialLawParams(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
// do not use different parameters in the test with inverted wettability
return myMaterialParams_;
}
/*!
* \brief Function for defining which phase is to be considered as the wetting phase.
*
* \return the wetting phase index
* \param globalPos The global position
*/
template<class FluidSystem>
int wettingPhaseAtPos(const GlobalPosition& globalPos) const
{
return FluidSystem::phase0Idx;
}
//! 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 initPermeability_;
Scalar initPorosity_;
MaterialLawParams myMaterialParams_;
};
} // end namespace Dumux
#endif
dune_symlink_to_source_files(FILES "el2p.input")
dune_add_test(NAME test_el2p
SOURCES test_el2p.cc
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/el2p_2p_reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/twop-00016.vtu
${CMAKE_SOURCE_DIR}/test/references/el2p_u_reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/poroelastic-00016.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_el2p el2p.input")
set(CMAKE_BUILD_TYPE Release)
#install sources
install(FILES
2pproblem.hh
2pspatialparams.hh
poroelasticproblem.hh
poroelasticspatialparams.hh
test_el2p.cc
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/multidomain/poromechanics/el2p)
[TimeLoop]
DtInitial = 1000 # [s]
TEnd = 1000000 # [s]
[Grid]
LowerLeft = 0 0
UpperRight = 1 1
Cells = 10 10
[Problem]
InjectionRateGas = 1e-3
InjectionRateWater = 0.0
EnableGravity = false
[PoroElastic.Problem]
Name = poroelastic # name passed to the output routines
[TwoP.Problem]
Name = twop # 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 el2p 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/fluidsystems/1pgas.hh>
#include <dumux/material/fluidsystems/2pimmiscible.hh>
#include <dumux/material/components/h2o.hh>
#include <dumux/material/components/air.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>);
// Set the fluid system for TwoPSubProblem
SET_PROP(PoroElasticSubTypeTag, FluidSystem)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::H2O<Scalar> >;
using NonwettingPhase = FluidSystems::OnePGas<Scalar, Components::Air<Scalar> >;
using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
};
// 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 el2p coupled problem.
*/
template<class TypeTag>
class PoroElasticSubProblem : public GeomechanicsFVProblem<TypeTag>
{
using ParentType = GeomechanicsFVProblem<TypeTag>;
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
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
const auto& context = couplingManager().poroMechanicsCouplingContext();
// here, we know that the flow problem uses cell-centered finite volumes, thus,
// we simply take the volume variables of the scv (i.e. element) to obtain fluid properties
const auto& facetVolVars = (*context.pmFlowElemVolVars)[scv.elementIndex()];
Scalar wPhaseDensity = facetVolVars.density(FluidSystem::phase0Idx);
Scalar nPhaseDensity = facetVolVars.density(FluidSystem::phase1Idx);
Scalar Sw = facetVolVars.saturation(FluidSystem::phase0Idx);
Scalar Sn = facetVolVars.saturation(FluidSystem::phase1Idx);
return (wPhaseDensity * Sw + nPhaseDensity * Sn);
}
/*!
* \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
const auto& context = couplingManager().poroMechanicsCouplingContext();
// here, we know that the flow problem uses cell-centered finite volumes, thus,
// we simply take the volume variables of the element to obtain fluid properties
const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
const auto& facetVolVars = (*context.pmFlowElemVolVars)[eIdx];
Scalar pw = facetVolVars.pressure(FluidSystem::phase0Idx);
Scalar pn = facetVolVars.pressure(FluidSystem::phase1Idx);
Scalar Sw = facetVolVars.saturation(FluidSystem::phase0Idx);
Scalar Sn = facetVolVars.saturation(FluidSystem::phase1Idx);
return (pw * Sw + pn * Sn);
}
/*!
* \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