Hi all,
I was able to include nonisothermal effects in the fractures setup but I
am now stuck trying to use the TwoPNCNI or TwoPNC models instead of the
TwoPmodel. The code compile but the Newton solver doesn't converge and
the run stop.
I tried to modify the boundary condition in implementing some injection
as in the basic exercise but I get the same problem.
What else do I need to modify to use several components in this setup ?
Thank you for your help,
Suzon
On 20/03/2024 13:33, Martin Schneider wrote:
Dear Suzon,
did you change the the model in properties.hh?
If you want to include nonisothermal effects you have to use the
TwoPNI model instead of the TwoP model.
Best,
Martin
On 19.03.24 17:36, Suzon Jammes wrote:
Dear users,
To learn how to use dumux I am trying to modify a bit the exercise
about fractures
(https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/-/blob/master/exercises/exercise-fractures).
Today I wanted to make it nonisothermal as explained in
exercice-basic but I am getting the following error :
error: ‘temperatureIdx’ is not a member of
‘Dumux::MatrixSubProblem<Dumux::Properties::TTag::MatrixProblem>::Indices’
{aka ‘Dumux::TwoPIndices’}
Is there anybody who knows how to solve this error ?
Thank you for your help,
Suzon
_______________________________________________
DuMux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
// -*- 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 MultiDomain
* \ingroup MultiDomainFacet
* \ingroup TwoPTests
* \brief The sub-problem for the fracture domain in the exercise on two-phase flow in fractured porous media.
*/
#ifndef DUMUX_SUZON_FRACTURESTEST1_FRACTURE_PROBLEM_HH
#define DUMUX_SUZON_FRACTURESTEST1_FRACTURE_PROBLEM_HH
// include the base problem and properties we inherit from
//#include <dumux/porousmediumflow/nonisothermal/indices.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/boundarytypes.hh>
#include <dumux/common/numeqvector.hh>
namespace Dumux {
/*!
* \ingroup MultiDomain
* \ingroup MultiDomainFacet
* \ingroup TwoPTests
* \brief The sub-problem for the fracture domain in the exercise on two-phase flow in fractured porous media.
*/
template<class TypeTag>
class FractureSubProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using PrimaryVariables = typename GridVariables::PrimaryVariables;
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
using Scalar = typename GridVariables::Scalar;
using GridGeometry = typename GridVariables::GridGeometry;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using GridView = typename GridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
// some indices for convenience
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
static constexpr int pressureIdx = Indices::pressureIdx;
static constexpr int saturationIdx = Indices::saturationIdx;
static constexpr int temperatureIdx = Indices::temperatureIdx;
public:
//! The constructor
FractureSubProblem(std::shared_ptr<const GridGeometry> gridGeometry,
std::shared_ptr<typename ParentType::SpatialParams> spatialParams,
const std::string& paramGroup)
: ParentType(gridGeometry, spatialParams, paramGroup)
, aperture_(getParamFromGroup<Scalar>(paramGroup, "SpatialParams.Aperture"))
, fracturetip_(getParamFromGroup<Scalar>(paramGroup, "SpatialParams.Fracturetip"))
, isExercisePartA_(getParamFromGroup<bool>(paramGroup, "Problem.IsExercisePartA"))
, isExercisePartB_(getParamFromGroup<bool>(paramGroup, "Problem.IsExercisePartB"))
, isExercisePartC_(getParamFromGroup<bool>(paramGroup, "Problem.IsExercisePartC"))
{
// initialize the fluid system, i.e. the tabulation
// of water properties. Use the default p/T ranges.
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
FluidSystem::init();
}
//! Specifies the type of boundary condition at a given position
BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
{
BoundaryTypes values;
// We only use no-flow boundary conditions for all immersed fractures
// in the domain (fracture tips that do not touch the domain boundary)
// Otherwise, we would lose mass leaving across the fracture tips.
values.setAllNeumann();
// However, there is one fracture reaching the top boundary. For that
// fracture tip we set Dirichlet Bcs - only in the unmodified state of
// the exercise though. In parts a, b & c we consider Neumann here.
//if (!isExercisePartA_ && !isExercisePartB_ && !isExercisePartC_)
if (globalPos[1] > (this->gridGeometry().bBoxMax()[1]+fracturetip_) - 1e-6)
values.setAllDirichlet();
return values;
}
//! Evaluate the source term in a sub-control volume of an element
NumEqVector source(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume& scv) const
{
// evaluate sources from bulk domain using the function in the coupling manager
auto source = couplingManagerPtr_->evalSourcesFromBulk(element, fvGeometry, elemVolVars, scv);
// these sources are in kg/s, divide by volume and extrusion to have it in kg/s/m³
source /= scv.volume()*elemVolVars[scv].extrusionFactor();
return source;
}
//! evaluates the Dirichlet boundary condition for a given position
PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{ return initialAtPos(globalPos); }
//! evaluate the initial conditions
PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{
// For the grid used here, the height of the domain is equal
// to the maximum y-coordinate
const auto domainHeight = (this->gridGeometry().bBoxMax()[1]+fracturetip_);
// we assume a constant water density of 1000 for initial conditions!
const auto& g = this->spatialParams().gravity(globalPos);
PrimaryVariables values;
values.setState(1);
Scalar densityW = 1000.0;
values[pressureIdx] = 1e5 - (domainHeight - globalPos[1])*densityW*g[1];
//std::cout << Fmt::format("pressure {:.2g}.\n", values[pressureIdx]);
//std::cout << Fmt::format("position {:.2g}.\n", globalPos[1]);
//std::cout << Fmt::format("gravity {:.2g}.\n", g[1]);
values[saturationIdx] = 0.0;
//values[temperatureIdx] = 283.0 + (domainHeight - globalPos[1])*0.026;
return values;
}
//! sets the pointer to the coupling manager.
void setCouplingManager(std::shared_ptr<CouplingManager> cm)
{ couplingManagerPtr_ = cm; }
//! returns reference to the coupling manager.
const CouplingManager& couplingManager() const
{ return *couplingManagerPtr_; }
private:
std::shared_ptr<CouplingManager> couplingManagerPtr_;
Scalar aperture_;
Scalar fracturetip_;
bool isExercisePartA_;
bool isExercisePartB_;
bool isExercisePartC_;
};
} // 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 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 MultiDomain
* \ingroup MultiDomainFacet
* \ingroup TwoPTests
* \brief The spatial parameters for the fracture sub-domain in the exercise
* on two-phase flow in fractured porous media.
*/
#ifndef DUMUX_SUZON_FRACTURESTEST1_FRACTURE_SPATIALPARAMS_HH
#define DUMUX_SUZON_FRACTURESTEST1_FRACTURE_SPATIALPARAMS_HH
#include <dumux/io/grid/griddata.hh>
#include <dumux/porousmediumflow/fvspatialparamsmp.hh>
#include <dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh>
namespace Dumux {
/*!
* \ingroup MultiDomain
* \ingroup MultiDomainFacet
* \ingroup TwoPTests
* \brief The spatial params the two-phase facet coupling test
*/
template< class GridGeometry, class Scalar >
class FractureSpatialParams
: public FVPorousMediumFlowSpatialParamsMP< GridGeometry, Scalar, FractureSpatialParams<GridGeometry, Scalar> >
{
using ThisType = FractureSpatialParams< GridGeometry, Scalar >;
using ParentType = FVPorousMediumFlowSpatialParamsMP< GridGeometry, Scalar, ThisType >;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using GridView = typename GridGeometry::GridView;
using Grid = typename GridView::Grid;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
// use a van-genuchten fluid matrix interaction relationship
using PcKrSwCurve = FluidMatrix::VanGenuchtenDefault<Scalar>;
// we identify those fractures as barriers, that have a domain marker
// of 2. This is what is set in the grid file (see grids/complex.geo)
static constexpr int barriersDomainMarker = 2;
public:
//! export the type used for permeabilities
using PermeabilityType = Scalar;
//! the constructor
FractureSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry,
std::shared_ptr<const Dumux::GridData<Grid>> gridData,
const std::string& paramGroup)
: ParentType(gridGeometry)
, gridDataPtr_(gridData)
, isExercisePartA_(getParamFromGroup<bool>(paramGroup, "Problem.IsExercisePartA"))
, isExercisePartB_(getParamFromGroup<bool>(paramGroup, "Problem.IsExercisePartB"))
, isExercisePartC_(getParamFromGroup<bool>(paramGroup, "Problem.IsExercisePartC"))
, pcKrSwCurve_("Fracture.SpatialParams")
, PcKrSwCurveBarrier_("Fracture.SpatialParams.Barrier")
{
porosity_ = getParamFromGroup<Scalar>(paramGroup, "SpatialParams.Porosity");
porosityBarrier_ = getParamFromGroup<Scalar>(paramGroup, "SpatialParams.PorosityBarrier");
permeability_ = getParamFromGroup<Scalar>(paramGroup, "SpatialParams.Permeability");
permeabilityBarrier_ = getParamFromGroup<Scalar>(paramGroup, "SpatialParams.PermeabilityBarrier");
aperture_ = getParamFromGroup<Scalar>(paramGroup, "SpatialParams.Aperture");
}
//! 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
{
// in the unmodified state and exercise part a return the non-barrier parameters
if (!isExercisePartB_ && !isExercisePartC_)
return permeability_;
// in exercise part b always return the barrier parameters
else if (isExercisePartB_)
return permeabilityBarrier_;
// in exercise part 3 return parameters depending on domain marker
else
{
if (getElementDomainMarker(element) == barriersDomainMarker)
return permeabilityBarrier_;
else
return permeability_;
}
}
//! Return the porosity
template< class ElementSolution >
Scalar porosity(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
// in the unmodified state and exercise part a return the non-barrier parameters
if (!isExercisePartB_ && !isExercisePartC_)
return porosity_;
// in exercise part b always return the barrier parameters
else if (isExercisePartB_)
return porosityBarrier_;
// in exercise part 3 return parameters depending on domain marker
else
{
if (getElementDomainMarker(element) == barriersDomainMarker)
return porosityBarrier_;
else
return porosity_;
}
}
/*!
* \brief Returns the fluid-matrix interaction law for the sub-control volume
*
* \param element The current finite element
* \param scv The sub-control volume
* \param elemSol The current element solution
*/
template<class ElementSolution>
auto fluidMatrixInteraction(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
// in the unmodified state and exercise part a return the non-barrier parameters
if (!isExercisePartB_ && !isExercisePartC_)
return makeFluidMatrixInteraction(pcKrSwCurve_);
// in exercise part b always return the barrier parameters
else if (isExercisePartB_)
return makeFluidMatrixInteraction(PcKrSwCurveBarrier_);
// in exercise part 3 return parameters depending on domain marker
else
{
if (getElementDomainMarker(element) == barriersDomainMarker)
return makeFluidMatrixInteraction(PcKrSwCurveBarrier_);
else
return makeFluidMatrixInteraction(pcKrSwCurve_);
}
}
//! Water is the wetting phase
template< class FluidSystem >
int wettingPhaseAtPos(const GlobalPosition& globalPos) const
{
// we set water as the wetting phase here
// which is phase0Idx in the H2oN2 fluid system
return FluidSystem::phase0Idx;
}
/*!
* \brief Returns the temperature at the domain at the given position
* \param globalPos The position in global coordinates where the temperature should be specified
*/
Scalar temperatureAtPos(const GlobalPosition& globalPos) const
{
return 283.15;
}
//! Set the aperture as extrusion factor.
Scalar extrusionFactorAtPos(const GlobalPosition& globalPos) const
{
// We treat the fractures as lower-dimensional in the grid,
// but we have to give it the aperture as extrusion factor
// such that the dimensions are correct in the end.
return aperture_;
}
//! returns the domain marker for an element
int getElementDomainMarker(const Element& element) const
{ return gridDataPtr_->getElementDomainMarker(element); }
private:
//! pointer to the grid data (contains domain markers)
std::shared_ptr<const Dumux::GridData<Grid>> gridDataPtr_;
bool isExercisePartA_;
bool isExercisePartB_;
bool isExercisePartC_;
const PcKrSwCurve pcKrSwCurve_;
const PcKrSwCurve PcKrSwCurveBarrier_;
Scalar porosity_;
Scalar porosityBarrier_;
Scalar aperture_;
PermeabilityType permeability_;
PermeabilityType permeabilityBarrier_;
};
} // 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 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
* \brief Test for the exercise on two-phase flow in fractured porous media.
*/
#include <config.h>
#include <iostream>
// include the properties header
#include "properties.hh"
#include <dumux/common/initialize.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/timeloop.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/linear/istlsolvers.hh>
#include <dumux/linear/linearalgebratraits.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <dumux/multidomain/newtonsolver.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/traits.hh>
#include <dumux/multidomain/facet/gridmanager.hh>
#include <dumux/multidomain/facet/couplingmapper.hh>
#include <dumux/multidomain/facet/couplingmanager.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/restart.hh>
// Define some types for this test so that we can set them as properties below and
// reuse them again in the main function with only one single definition of them here
using MatrixTypeTag = Dumux::Properties::TTag::MatrixProblem;
using FractureTypeTag = Dumux::Properties::TTag::FractureProblem;
using MatrixGridGeometry = Dumux::GetPropType<MatrixTypeTag, Dumux::Properties::GridGeometry>;
using FractureGridGeometry = Dumux::GetPropType<FractureTypeTag, Dumux::Properties::GridGeometry>;
using TheMultiDomainTraits = Dumux::MultiDomainTraits<MatrixTypeTag, FractureTypeTag>;
using TheCouplingMapper = Dumux::FacetCouplingMapper<MatrixGridGeometry, FractureGridGeometry>;
using TheCouplingManager = Dumux::FacetCouplingManager<TheMultiDomainTraits, TheCouplingMapper>;
// set the coupling manager property in the sub-problems
namespace Dumux::Properties {
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::MatrixProblem> { using type = TheCouplingManager; };
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::FractureProblem> { using type = TheCouplingManager; };
} // end namespace Dumux::Properties
// main program
int main(int argc, char** argv)
{
using namespace Dumux;
// initialize MPI+x, finalize is done automatically on exit
Dumux::initialize(argc, argv);
// initialize parameter tree
Parameters::init(argc, argv);
// We use the grid manager from the facet coupling framework (see alias above main)
// This requires the grids used to be passed as template arguments, where
// they are assumed to be ordered in descending grid dimension. Thus,
// we pass the matrix grid as first and the fracture grid as second argument.
using MatrixGrid = GetPropType<MatrixTypeTag, Properties::Grid>;
using FractureGrid = GetPropType<FractureTypeTag, Properties::Grid>;
using GridManager = Dumux::FacetCouplingGridManager<MatrixGrid, FractureGrid>;
GridManager gridManager;
// try to create a grid (from the grid file in the input file)
// Init() creates the grid from the grid file specified in the input file.
// This works with a single grid file in which in addition to the matrix
// (d-dimensional) elements also the (d-1)-dimensional elements are specified
// which are interpreted as the fracture elements. See the .geo file in the
// grids folder on how to create such grids using gmsh. Note that currently
// only the Gmsh mesh format (.msh) is supported!
gridManager.init();
gridManager.loadBalance();
// we compute on the leaf grid views (get them from grid manager)
// the grid ids correspond to the order of the grids passed to the manager (see above)
static constexpr std::size_t matrixGridId = 0;
static constexpr std::size_t fractureGridId = 1;
const auto& matrixGridView = gridManager.template grid<matrixGridId>().leafGridView();
const auto& fractureGridView = gridManager.template grid<fractureGridId>().leafGridView();
// create the finite volume grid geometries
auto matrixFvGridGeometry = std::make_shared<MatrixGridGeometry>(matrixGridView);
auto fractureFvGridGeometry = std::make_shared<FractureGridGeometry>(fractureGridView);
// the problems (boundary/initial conditions etc)
using MatrixProblem = GetPropType<MatrixTypeTag, Properties::Problem>;
using FractureProblem = GetPropType<FractureTypeTag, Properties::Problem>;
// pass the model parameter group to the spatial params so that they obtain the right
// values from the input file since we use groups for matrix and fracture
// We also want to use domain markers in this exercise. For this reason, we also pass
// the grid data object from the grid manager to them, so that they have access to the
// domain markers that were specified in the grid file.
auto matrixGridData = gridManager.getGridData()->template getSubDomainGridData<matrixGridId>();
auto matrixSpatialParams = std::make_shared<typename MatrixProblem::SpatialParams>(matrixFvGridGeometry, matrixGridData, "Matrix");
auto matrixProblem = std::make_shared<MatrixProblem>(matrixFvGridGeometry, matrixSpatialParams, "Matrix");
// extract domain height from the matrix problem and pass to fracture problem (needed for right initial pressure distribution)
auto fractureGridData = gridManager.getGridData()->template getSubDomainGridData<fractureGridId>();
auto fractureSpatialParams = std::make_shared<typename FractureProblem::SpatialParams>(fractureFvGridGeometry, fractureGridData, "Fracture");
auto fractureProblem = std::make_shared<FractureProblem>(fractureFvGridGeometry, fractureSpatialParams, "Fracture");
// the solution vector
using SolutionVector = typename TheMultiDomainTraits::SolutionVector;
SolutionVector x, xOld;
// The domain ids within the multi-domain framework.
// They do not necessarily have to be the same as the grid ids
// in case you have more subdomains involved. We domain ids
// correspond to the order of the type tags passed to the multidomain
// traits (see definition of the traits class at the beginning of this file)
static const auto matrixDomainId = typename TheMultiDomainTraits::template SubDomain<0>::Index();
static const auto fractureDomainId = typename TheMultiDomainTraits::template SubDomain<1>::Index();
// resize the solution vector and write initial solution to it
x[matrixDomainId].resize(matrixFvGridGeometry->numDofs());
x[fractureDomainId].resize(fractureFvGridGeometry->numDofs());
matrixProblem->applyInitialSolution(x[matrixDomainId]);
fractureProblem->applyInitialSolution(x[fractureDomainId]);
// instantiate the class holding the coupling maps between the domains
// this needs the information on embeddings (connectivity between matrix
// and fracture domain). This information is extracted directly from the
// grid during file read and can therefore be obtained from the grid manager.
const auto embeddings = gridManager.getEmbeddings();
auto couplingMapper = std::make_shared<TheCouplingMapper>();
couplingMapper->update(*matrixFvGridGeometry, *fractureFvGridGeometry, embeddings);
// the coupling manager (needs the coupling mapper)
auto couplingManager = std::make_shared<TheCouplingManager>();
couplingManager->init(matrixProblem, fractureProblem, couplingMapper, x);
// we have to set coupling manager pointer in sub-problems
// they also have to be made accessible in them (see e.g. matrixproblem.hh)
matrixProblem->setCouplingManager(couplingManager);
fractureProblem->setCouplingManager(couplingManager);
// the grid variables
using MatrixGridVariables = GetPropType<MatrixTypeTag, Properties::GridVariables>;
using FractureGridVariables = GetPropType<FractureTypeTag, Properties::GridVariables>;
auto matrixGridVariables = std::make_shared<MatrixGridVariables>(matrixProblem, matrixFvGridGeometry);
auto fractureGridVariables = std::make_shared<FractureGridVariables>(fractureProblem, fractureFvGridGeometry);
xOld = x;
matrixGridVariables->init(x[matrixDomainId]);
fractureGridVariables->init(x[fractureDomainId]);
// initialize the vtk output modules
using MatrixVtkOutputModule = VtkOutputModule<MatrixGridVariables, GetPropType<MatrixTypeTag, Properties::SolutionVector>>;
using FractureVtkOutputModule = VtkOutputModule<FractureGridVariables, GetPropType<FractureTypeTag, Properties::SolutionVector>>;
MatrixVtkOutputModule matrixVtkWriter(*matrixGridVariables, x[matrixDomainId], matrixProblem->name());
FractureVtkOutputModule fractureVtkWriter(*fractureGridVariables, x[fractureDomainId], fractureProblem->name());
// Add model specific output fields
using MatrixIOFields = GetPropType<MatrixTypeTag, Properties::IOFields>;
using FractureIOFields = GetPropType<FractureTypeTag, Properties::IOFields>;
MatrixIOFields::initOutputModule(matrixVtkWriter);
FractureIOFields::initOutputModule(fractureVtkWriter);
// add domain markers to output
std::vector<int> matrixDomainMarkers(matrixFvGridGeometry->gridView().size(0));
for (const auto& element : elements(matrixFvGridGeometry->gridView()))
matrixDomainMarkers[matrixFvGridGeometry->elementMapper().index(element)] = matrixProblem->spatialParams().getElementDomainMarker(element);
matrixVtkWriter.addField(matrixDomainMarkers, "domainMarker");
std::vector<int> fractureDomainMarkers(fractureFvGridGeometry->gridView().size(0));
for (const auto& element : elements(fractureFvGridGeometry->gridView()))
fractureDomainMarkers[fractureFvGridGeometry->elementMapper().index(element)] = fractureProblem->spatialParams().getElementDomainMarker(element);
fractureVtkWriter.addField(fractureDomainMarkers, "domainMarker");
// write out initial solution
matrixVtkWriter.write(0.0);
fractureVtkWriter.write(0.0);
// get some time loop parameters
const auto tEnd = getParam<double>("TimeLoop.TEnd");
const auto maxDt = getParam<double>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<double>("TimeLoop.DtInitial");
// instantiate time loop
auto timeLoop = std::make_shared< TimeLoop<double> >(/*startTime*/0.0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
matrixProblem->setTime(timeLoop->time()+timeLoop->timeStepSize() );
matrixProblem->setTimeStepSize(timeLoop->timeStepSize());
// the assembler for the coupled problem
using Assembler = MultiDomainFVAssembler<TheMultiDomainTraits, TheCouplingManager, DiffMethod::numeric, /*implicit?*/true>;
auto assembler = std::make_shared<Assembler>( std::make_tuple(matrixProblem, fractureProblem),
std::make_tuple(matrixFvGridGeometry, fractureFvGridGeometry),
std::make_tuple(matrixGridVariables, fractureGridVariables),
couplingManager,
timeLoop, xOld);
// the linear solver
using LinearSolver = ILUBiCGSTABIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonSolver = Dumux::MultiDomainNewtonSolver<Assembler, LinearSolver, TheCouplingManager>;
auto newtonSolver = std::make_shared<NewtonSolver>(assembler, linearSolver, couplingManager);
// time loop
timeLoop->start(); do
{
matrixProblem->setTime(timeLoop->time()+timeLoop->timeStepSize() );
matrixProblem->setTimeStepSize(timeLoop->timeStepSize());
// solve the non-linear system with time step control
newtonSolver->solve(x, *timeLoop);
// make the new solution the old solution
xOld = x;
matrixGridVariables->advanceTimeStep();
fractureGridVariables->advanceTimeStep();
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
// write vtk output
matrixVtkWriter.write(timeLoop->time());
fractureVtkWriter.write(timeLoop->time());
// report statistics of this time step
timeLoop->reportTimeStep();
// set new dt as suggested by the Newton solver
timeLoop->setTimeStepSize(newtonSolver->suggestTimeStepSize(timeLoop->timeStepSize()));
} while (!timeLoop->finished());
// output some Newton statistics
newtonSolver->report();
// report time loop statistics
timeLoop->finalize();
Parameters::print();
return 0;
}// end main
// -*- 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 MultiDomain
* \ingroup MultiDomainFacet
* \ingroup TwoPTests
* \brief The sub-problem for the matrix domain in the exercise on two-phase flow in fractured porous media.
*/
#ifndef DUMUX_SUZON_FRACTURESTEST1_MATRIX_PROBLEM_HH
#define DUMUX_SUZON_FRACTURESTEST1_MATRIX_PROBLEM_HH
// we need this in this test in order to define the domain
// id of the fracture problem (see function interiorBoundaryTypes())
#include <dune/common/indices.hh>
// include the base problem and properties we inherit from
#include <dumux/porousmediumflow/problem.hh>
//#include <dumux/porousmediumflow/nonisothermal/indices.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/boundarytypes.hh>
#include <dumux/common/numeqvector.hh>
namespace Dumux {
/*!
* \ingroup MultiDomain
* \ingroup MultiDomainFacet
* \ingroup TwoPTests
* \brief The sub-problem for the matrix domain in the exercise on two-phase flow in fractured porous media.
*/
template<class TypeTag>
class MatrixSubProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using PrimaryVariables = typename GridVariables::PrimaryVariables;
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
using Scalar = typename GridVariables::Scalar;
using GridGeometry = typename GridVariables::GridGeometry;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using GridView = typename GridGeometry::GridView;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
// some indices for convenience
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
static constexpr int pressureIdx = Indices::pressureIdx;
static constexpr int saturationIdx = Indices::saturationIdx;
static constexpr int temperatureIdx = Indices::temperatureIdx;
public:
//! The constructor
MatrixSubProblem(std::shared_ptr<const GridGeometry> gridGeometry,
std::shared_ptr<typename ParentType::SpatialParams> spatialParams,
const std::string& paramGroup = "")
: ParentType(gridGeometry, spatialParams, paramGroup)
, boundaryOverPressure_(getParamFromGroup<Scalar>(paramGroup, "Problem.BoundaryOverPressure"))
, boundarySaturation_(getParamFromGroup<Scalar>(paramGroup, "Problem.BoundarySaturation"))
, boundaryTimeChange_(getParamFromGroup<Scalar>(paramGroup, "TimeLoop.BoundaryTimeChange"))
{
// initialize the fluid system, i.e. the tabulation
// of water properties. Use the default p/T ranges.
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
FluidSystem::init();
}
//! Specifies the type of boundary condition at a given position
BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
{
BoundaryTypes values;
// in the unmodified state, we consider buoancy-driven upwards migration
// of hydrogen and set Dirichlet BCs on the top and bottom boundary
values.setAllNeumann();
if (globalPos[1] > this->gridGeometry().bBoxMax()[1] - 1e-6)
values.setAllDirichlet();
return values;
}
//! Specifies the type of interior boundary condition to be used on a sub-control volume face
BoundaryTypes interiorBoundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const
{
BoundaryTypes values;
// Here we set the type of condition to be used on faces that coincide
// with a fracture. If Neumann is specified, a flux continuity condition
// on the basis of the normal fracture permeability is evaluated. If this
// permeability is lower than that of the matrix, this approach is able to
// represent the resulting pressure jump across the fracture. If Dirichlet is set,
// the pressure jump across the fracture is neglected and the pressure inside
// the fracture is directly applied at the interface between fracture and matrix.
// This assumption is justified for highly-permeable fractures, but lead to erroneous
// results for low-permeable fractures.
// Here, we consider "open" fractures for which we cannot define a normal permeability
// and for which the pressure jump across the fracture is neglectable. Thus, we set
// the interior boundary conditions to Dirichlet.
// IMPORTANT: Note that you will never be asked to set any values at the interior boundaries!
// This simply chooses a different interface condition!
//values.setAllDirichlet();
values.setAllNeumann();
return values;
}
//! evaluates the Dirichlet boundary condition for a given position
PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{
return initialAtPos(globalPos);
}
/*!
* \brief Evaluate the boundary conditions for a Neumann
* boundary segment.
*
* \param globalPos The position of the integration point of the boundary segment.
*/
NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
{
// initialize values to zero, i.e. no-flow Neumann boundary conditions
NumEqVector values(0.0);
// if we are inside the injection zone set inflow Neumann boundary conditions
if (onInjectionBoundary(globalPos))
{
// inject hydrogen. negative values mean injection
// convert from units kg/(s*m^2) to mole/(s*m^2)
values[Indices::conti0EqIdx + FluidSystem::H2Idx] = -1e-4/FluidSystem::molarMass(FluidSystem::H2Idx);
values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0;
}
return values;
}
//! evaluate the initial conditions
PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{
// For the grid used here, the height of the domain is equal
// to the maximum y-coordinate
const auto domainHeight = this->gridGeometry().bBoxMax()[1];
// we assume a constant water density of 1000 for initial conditions!
const auto& g = this->spatialParams().gravity(globalPos);
PrimaryVariables values;
values.setState(Indices::firstPhaseOnly);
//values.setState(1);
Scalar densityW = 1000.0;
values[pressureIdx] = 1e5 - (domainHeight - globalPos[1])*densityW*g[1];
//values[saturationIdx] = 0.0;
//values[temperatureIdx] = 283.0 + (domainHeight - globalPos[1])*0.026;
// initially we have some nitrogen dissolved
// saturation mole fraction would be
//const Scalar moleFracLiquidH2 = densityW*0.95/BinaryCoeff::H2O_H2::henry(this->spatialParams().temperatureAtPos(globalPos));
// note that because we start with a single phase system the primary variables
// are pl and x^w_H2. This will switch as soon after we start injecting to a two
// phase system so the primary variables will be pl and Sn (nonwetting saturation).
//values[Indices::switchIdx] = moleFracLiquidH2;
return values;
}
//! sets the pointer to the coupling manager.
void setCouplingManager(std::shared_ptr<CouplingManager> cm)
{ couplingManagerPtr_ = cm; }
//! returns reference to the coupling manager.
const CouplingManager& couplingManager() const
{ return *couplingManagerPtr_; }
void setTime(const Scalar t)
{
time_ = t;
}
// We need the time step size to regularize the reactive source terms.
void setTimeStepSize(const Scalar dt)
{ timeStepSize_ = dt; }
//! Return true if the given position is in the injection boundary region
bool onInjectionBoundary(const GlobalPosition& globalPos) const
{
return globalPos[0] < 750.0
&& globalPos[0] > 250.0
&& globalPos[1] < eps_;
}
private:
std::shared_ptr<CouplingManager> couplingManagerPtr_;
static constexpr Scalar eps_ = 1e-6;
Scalar boundaryOverPressure_;
Scalar boundarySaturation_;
Scalar time_ = 0.0;
Scalar timeStepSize_ = 0.0;
Scalar boundaryTimeChange_;
};
} // 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 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 MultiDomain
* \ingroup MultiDomainFacet
* \ingroup TwoPTests
* \brief The spatial parameters for the fracture sub-domain in the exercise
* on two-phase flow in fractured porous media.
*/
#ifndef DUMUX_SUZON_FRACTURESTEST1_MATRIX_SPATIALPARAMS_HH
#define DUMUX_SUZON_FRACTURESTEST1_MATRIX_SPATIALPARAMS_HH
#include <dumux/io/grid/griddata.hh>
#include <dumux/porousmediumflow/fvspatialparamsmp.hh>
#include <dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh>
namespace Dumux {
/*!
* \ingroup MultiDomain
* \ingroup MultiDomainFacet
* \ingroup TwoPTests
* \brief The spatial params the two-phase facet coupling test
*/
template< class GridGeometry, class Scalar >
class MatrixSpatialParams
: public FVPorousMediumFlowSpatialParamsMP< GridGeometry, Scalar, MatrixSpatialParams<GridGeometry, Scalar> >
{
using ThisType = MatrixSpatialParams< GridGeometry, Scalar >;
using ParentType = FVPorousMediumFlowSpatialParamsMP< GridGeometry, Scalar, ThisType >;
using GridView = typename GridGeometry::GridView;
using Grid = typename GridView::Grid;
// get the dimensions of the simulation domain from GridView
static const int dimWorld = GridView::dimensionworld;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
// use a van-genuchten fluid matrix interaction relationship
using PcKrSwCurve = FluidMatrix::VanGenuchtenDefault<Scalar>;
public:
//! export the type used for permeabilities
using PermeabilityType = Scalar;
//! the constructor
MatrixSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry,
std::shared_ptr<const Dumux::GridData<Grid>> gridData,
const std::string& paramGroup)
: ParentType(gridGeometry)
, gridDataPtr_(gridData)
, aquitardPcKrSwCurve_("Matrix.SpatialParams.Aquitard")
, aquiferPcKrSwCurve_("Matrix.SpatialParams.Aquifer")
{
aquiferHeightFromBottom_ = 500.0;
// intrinsic permeabilities
aquitardK_ = getParamFromGroup<Scalar>(paramGroup, "SpatialParams.PermeabilityAquitard");
aquiferK_ = getParamFromGroup<Scalar>(paramGroup, "SpatialParams.PermeabilityAquifer");
// porosities
aquitardPorosity_ = getParamFromGroup<Scalar>(paramGroup, "SpatialParams.PorosityAquitard");
aquiferPorosity_ = getParamFromGroup<Scalar>(paramGroup, "SpatialParams.PorosityAquifer");
//porosity_ = getParamFromGroup<Scalar>(paramGroup, "SpatialParams.Porosity");
//permeability_ = getParamFromGroup<Scalar>(paramGroup, "SpatialParams.Permeability");
}
//! Function for defining the (intrinsic) permeability \f$[m^2]\f$.
PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
{
// here, either aquitard or aquifer permeability are returned, depending on the global position
if (isInAquitard_(globalPos))
return aquitardK_;
return aquiferK_;
}
//! Return the porosity
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{
// here, either aquitard or aquifer porosity are returned, depending on the global position
if (isInAquitard_(globalPos))
return aquitardPorosity_;
return aquiferPorosity_;
}
/*!
* \brief Returns the fluid-matrix interaction law at a given location
* \param globalPos The global coordinates for the given location
*/
auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const
{
if (isInAquitard_(globalPos))
return makeFluidMatrixInteraction(aquitardPcKrSwCurve_);
return makeFluidMatrixInteraction(aquiferPcKrSwCurve_);
}
//! Water is the wetting phase
template< class FluidSystem >
int wettingPhaseAtPos(const GlobalPosition& globalPos) const
{
// we set water as the wetting phase here
// which is phase0Idx in the H2oH2 fluid system
return FluidSystem::phase0Idx;
}
/*!
* \brief Returns the temperature at the domain at the given position
* \param globalPos The position in global coordinates where the temperature should be specified
*/
Scalar temperatureAtPos(const GlobalPosition& globalPos) const
{
return 283.0 ;
}
//! returns the domain marker for an element
int getElementDomainMarker(const Element& element) const
{ return gridDataPtr_->getElementDomainMarker(element); }
private:
static constexpr Scalar eps_ = 1e-6;
// provides a convenient way distinguishing whether a given location is inside the aquitard
bool isInAquitard_(const GlobalPosition &globalPos) const
{
// globalPos[dimWorld-1] is the y direction for 2D grids or the z direction for 3D grids
return globalPos[dimWorld-1] > aquiferHeightFromBottom_ + eps_;
}
private:
//! pointer to the grid data (contains domain markers)
std::shared_ptr<const Dumux::GridData<Grid>> gridDataPtr_;
Scalar aquitardK_;
Scalar aquiferK_;
Scalar aquiferHeightFromBottom_;
Scalar aquitardPorosity_;
Scalar aquiferPorosity_;
const PcKrSwCurve aquitardPcKrSwCurve_;
const PcKrSwCurve aquiferPcKrSwCurve_;
};
} // end namespace Dumux
#endif
[TimeLoop]
TEnd = 15000000 # [s]
DtInitial = 10 # [s]
MaxTimeStepSize = 2500 # [s]
BoundaryTimeChange = 7.5e6
[Problem]
EnableGravity = true
IsExercisePartA = false
IsExercisePartB = false
IsExercisePartC = false
[Grid]
File = ./grids/complex_TS1.msh
DomainMarkers = true # enable domain markers
Grid.Overlap = true
[Matrix]
Problem.Name = AmatrixTS1_CH4_test
Problem.BoundaryOverPressure = 5e5
Problem.BoundarySaturation = 0.2
SpatialParams.PermeabilityAquitard = 1e-15
SpatialParams.PermeabilityAquifer = 1e-12
SpatialParams.PorosityAquitard = 0.2
SpatialParams.PorosityAquifer = 0.4
SpatialParams.VanGenuchtenAlpha = 1e-3
SpatialParams.VanGenuchtenN = 3
SpatialParams.Snr = 0.0
SpatialParams.Swr = 0.0
[Fracture]
Problem.Name = AfracturesTS1_CH4_test
SpatialParams.Aperture = 1e-1
SpatialParams.Fracturetip = 0
SpatialParams.Permeability = 1e-7
SpatialParams.PermeabilityBarrier = 1e-16
SpatialParams.Porosity = 0.85
SpatialParams.PorosityBarrier = 0.05
SpatialParams.VanGenuchtenAlpha = 1e-1
SpatialParams.VanGenuchtenN = 2
SpatialParams.Swr = 0.0
SpatialParams.Snr = 0.0
SpatialParams.Barrier.VanGenuchtenAlpha = 1e-4
SpatialParams.Barrier.VanGenuchtenN = 2.5
SpatialParams.Barrier.Snr = 0.0
SpatialParams.Barrier.Swr = 0.0
[Component]
SolidDensity = 2700 # solid density of granite
SolidThermalConductivity = 2.8 # solid thermal conducitivity of granite
SolidHeatCapacity = 790 # solid heat capacity of granite
// -*- 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 MultiDomain
* \ingroup MultiDomainFacet
* \ingroup TwoPTests
* \brief The properties file for exercise on two-phase flow in fractured porous media.
*/
#ifndef DUMUX_SUZON_FRACTURESTESTA_PROPERTIES_HH
#define DUMUX_SUZON_FRACTURESTESTA_PROPERTIES_HH
// Both sub-problems
// include the model we inherit from
#include <dumux/porousmediumflow/2pnc/model.hh>
// we want to simulate hydrogen gas transport in a water-saturated medium
#include <dumux/material/fluidsystems/h2oh2.hh>
//#include <dumux/material/fluidsystems/h2och4.hh>
//#include <dumux/material/fluidsystems/h2on2o2.hh>
// Fracture sub-problem
// we use foam grid for the discretization of the fracture domain
// as this grid manager is able to represent network/surface grids
#include <dune/foamgrid/foamgrid.hh>
// we use a cell-centered finite volume scheme with tpfa here
#include <dumux/discretization/cctpfa.hh>
// the spatial parameters (permeabilities, material parameters etc.)
#include "fracturespatialparams.hh"
// the fracture sub-problem problem file
#include "fractureproblem.hh"
// Matrix sub-problem
// the spatial parameters (permeabilities, material parameters etc.)
#include "matrixspatialparams.hh"
// the matrix sub-problem problem file
#include "matrixproblem.hh"
// we use alu grid for the discretization of the matrix domain
#include <dune/alugrid/grid.hh>
// We are using the framework for models that consider coupling
// across the element facets of the bulk domain. This has some
// properties defined, which we have to inherit here. In this
// exercise we want to use a cell-centered finite volume scheme
// with tpfa.
#include <dumux/multidomain/facet/cellcentered/tpfa/properties.hh>
namespace Dumux::Properties {
// create the type tag node for the matrix and fracture sub-problems
namespace TTag {
struct MatrixProblem { using InheritsFrom = std::tuple<CCTpfaFacetCouplingModel, TwoPNC>; };
struct FractureProblem { using InheritsFrom = std::tuple<TwoPNC, CCTpfaModel>; };
} // end namespace TTag
// Set the grid type for the matrix and fracture sub-domains
template<class TypeTag>
struct Grid<TypeTag, TTag::MatrixProblem> { using type = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>; };
template<class TypeTag>
struct Grid<TypeTag, TTag::FractureProblem> { using type = Dune::FoamGrid<1, 2>; };
// Set the problem type for the matrix and fracture sub-domains
template<class TypeTag>
struct Problem<TypeTag, TTag::MatrixProblem> { using type = MatrixSubProblem<TypeTag>; };
template<class TypeTag>
struct Problem<TypeTag, TTag::FractureProblem> { using type = FractureSubProblem<TypeTag>; };
// set the spatial params for the matrix and fracture sub-domains
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::MatrixProblem>
{
using type = MatrixSpatialParams< GetPropType<TypeTag, Properties::GridGeometry>,
GetPropType<TypeTag, Properties::Scalar> >;
};
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::FractureProblem>
{
using type = FractureSpatialParams< GetPropType<TypeTag, Properties::GridGeometry>,
GetPropType<TypeTag, Properties::Scalar> >;
};
// the fluid system for the matrix and fracture sub-domains
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::MatrixProblem>
{
using type = Dumux::FluidSystems::H2OH2< GetPropType<TypeTag, Properties::Scalar>,
FluidSystems::H2OH2DefaultPolicy</*fastButSimplifiedRelations=*/true> >;
};
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::FractureProblem>
{
using type = Dumux::FluidSystems::H2OH2< GetPropType<TypeTag, Properties::Scalar>,
FluidSystems::H2OH2DefaultPolicy</*fastButSimplifiedRelations=*/true> >;
};
} // end namespace Dumux::Properties
#endif
_______________________________________________
DuMux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
