Hello,
I'm trying to implement a random distribution of K in my domain. I followed
the example of 1p provided in DuMux-test. But I am having the following
error :
*/home/ggl/dumux/dumux/test/porousmediumflow/co2/implicit/test_boxco2.cc:79:58:
required from here/usr/include/c++/7/ext/new_allocator.h:136:4: error: no
matching function for call to
‘Dumux::HeterogeneousSpatialParams<Dumux::Properties::TTag::HeterogeneousBoxProblem>::HeterogeneousSpatialParams(const
Dune::GridView<Dune::ALU3dLeafGridViewTraits<const Dune::ALUGrid<2, 2,
(Dune::ALUGridElementType)1, (Dune::ALUGridRefinementType)1>,
(Dune::PartitionIteratorType)4> >&)’ { ::new((void *)__p)
_Up(std::forward<_Args>(__args)...); }*
Do you have any suggestions on how can I fixe it.
PS : Please find attached the problem and spatialparameters files.
Regards,
Kenza
// -*- 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 Definition of a problem, where CO2 is injected in a reservoir.
*/
#ifndef DUMUX_HETEROGENEOUS_PROBLEM_HH
#define DUMUX_HETEROGENEOUS_PROBLEM_HH
#include <fstream>
#include <dumux/porousmediumflow/co2/implicit/model.hh>
#include <dumux/porousmediumflow/co2/implicit/volumevariables.hh>
#include <dumux/material/fluidsystems/brineco2.hh>
#include <dumux/porousmediumflow/implicit/problem.hh>
#include <dumux/implicit/box/intersectiontovertexbc.hh>
#include <dumux/linear/amgbackend.hh>
#include "heterogeneousspatialparameters.hh"
#include "heterogeneousco2tables.hh"
namespace Dumux
{
template <class TypeTag>
class HeterogeneousProblem;
namespace Properties
{
NEW_TYPE_TAG(HeterogeneousProblem, INHERITS_FROM(TwoPTwoC, HeterogeneousSpatialParams));
NEW_TYPE_TAG(HeterogeneousBoxProblem, INHERITS_FROM(BoxModel, HeterogeneousProblem));
NEW_TYPE_TAG(HeterogeneousCCProblem, INHERITS_FROM(CCModel, HeterogeneousProblem));
// Set the grid type
#if HAVE_DUNE_ALUGRID
SET_TYPE_PROP(HeterogeneousProblem, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>);
#else
SET_TYPE_PROP(HeterogeneousProblem, Grid, Dune::YaspGrid<2>);
#endif
// Set the problem property
SET_TYPE_PROP(HeterogeneousProblem, Problem, HeterogeneousProblem<TypeTag>);
// Set fluid configuration
SET_TYPE_PROP(HeterogeneousProblem, FluidSystem, BrineCO2FluidSystem<TypeTag>);
// Set the CO2 table to be used; in this case not the the default table
SET_TYPE_PROP(HeterogeneousProblem, CO2Table, HeterogeneousCO2Tables::CO2Tables);
// Set the salinity mass fraction of the brine in the reservoir
SET_SCALAR_PROP(HeterogeneousProblem, ProblemSalinity, 1e-3);
//! the CO2 Model and VolumeVariables properties
SET_TYPE_PROP(HeterogeneousProblem, Model, CO2Model<TypeTag>);
SET_TYPE_PROP(HeterogeneousProblem, VolumeVariables, CO2VolumeVariables<TypeTag>);
// Solver settings for the tests using AMG
SET_TYPE_PROP(HeterogeneousBoxProblem, LinearSolver, Dumux::AMGBackend<TypeTag> );
SET_TYPE_PROP(HeterogeneousCCProblem, LinearSolver, Dumux::AMGBackend<TypeTag> );
// Use Moles
SET_BOOL_PROP(HeterogeneousProblem, UseMoles, false);
}
/*!
* \ingroup CO2Model
* \ingroup ImplicitTestProblems
* \brief Definition of a problem, where CO2 is injected in a reservoir.
*
* The domain is sized 200m times 100m and consists of four layers, a
* permeable reservoir layer at the bottom, a barrier rock layer with reduced permeability, another reservoir layer
* and at the top a barrier rock layer with a very low permeablility.
*
* CO2 is injected at the permeable bottom layer
* from the left side. The domain is initially filled with brine.
*
* The grid is unstructered and permeability and porosity for the elements are read in from the grid file. The grid file
* also contains so-called boundary ids which can be used assigned during the grid creation in order to differentiate
* between different parts of the boundary.
* These boundary ids can be imported into the problem where the boundary conditions can then be assigned accordingly.
*
* The model is able to use either mole or mass fractions. The property useMoles can be set to either true or false in the
* problem file. Make sure that the according units are used in the problem setup. The default setting for useMoles is false.
*
* To run the simulation execute the following line in shell (works with the box and cell centered spatial discretization method):
* <tt>./test_ccco2 </tt> or <tt>./test_boxco2 </tt>
*/
template <class TypeTag >
class HeterogeneousProblem : public ImplicitPorousMediaProblem<TypeTag>
{
typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
enum {
// Grid and world dimension
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
// copy some indices for convenience
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum {
lPhaseIdx = Indices::wPhaseIdx,
gPhaseIdx = Indices::nPhaseIdx
};
enum {
wCompIdx = FluidSystem::wCompIdx,
nCompIdx = FluidSystem::nCompIdx
};
enum {
BrineIdx = FluidSystem::BrineIdx,
CO2Idx = FluidSystem::CO2Idx
};
enum {
conti0EqIdx = Indices::conti0EqIdx,
contiCO2EqIdx = conti0EqIdx + CO2Idx
};
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::Intersection Intersection;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef typename GET_PROP_TYPE(TypeTag, CO2Table) CO2Table;
typedef Dumux::CO2<Scalar, CO2Table> CO2;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
//! property that defines whether mole or mass fractions are used
static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
public:
/*,
: ParentType(timeManager, GridCreator::grid().leafGridView())/*,
//Boundary Id Setup:
injectionTop_(1),
injectionBottom_(2),
dirichletBoundary_(3),
noFlowBoundary_(4),
intersectionToVertexBC_(*this)
*/
/*!
* \brief The constructor
*
* \param timeManager The time manager
* \param gridView The grid view
*/
HeterogeneousProblem(TimeManager &timeManager,
const GridView &gridView)
: ParentType(timeManager, GridCreator::grid().leafGridView()) /*,
//Boundary Id Setup:
injectionTop_(1),
injectionBottom_(2),
dirichletBoundary_(3),
noFlowBoundary_(4),
intersectionToVertexBC_(*this) */
{
nTemperature_ = GET_RUNTIME_PARAM(TypeTag, int, FluidSystem.NTemperature);
nPressure_ = GET_RUNTIME_PARAM(TypeTag, int, FluidSystem.NPressure);
pressureLow_ = GET_RUNTIME_PARAM(TypeTag, Scalar, FluidSystem.PressureLow);
pressureHigh_ = GET_RUNTIME_PARAM(TypeTag, Scalar, FluidSystem.PressureHigh);
temperatureLow_ = GET_RUNTIME_PARAM(TypeTag, Scalar, FluidSystem.TemperatureLow);
temperatureHigh_ = GET_RUNTIME_PARAM(TypeTag, Scalar, FluidSystem.TemperatureHigh);
depthBOR_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.DepthBOR);
name_ = GET_RUNTIME_PARAM(TypeTag, std::string, Problem.Name);
injectionRate_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.InjectionRate);
timepulsedebut1_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.TimePulseDebut1);
timepulsefin1_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.TimePulseFin1);
timepulsedebut2_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.TimePulseDebut2);
timepulsefin2_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.TimePulseFin2);
resultsTime1_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.ResultsTime1);
resultsTime2_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.ResultsTime2);
resultsTime3_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.ResultsTime3);
injecttop_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.InjectTop);
injectbot_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.InjectBot);
injectright_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.InjectRight);
injectleft_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.InjectLeft);
celltemp_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.CellTemp);
hzdp_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.HZdp);
/* Alternative syntax:
* typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
* const Dune::ParameterTree &tree = ParameterTree::tree();
* nTemperature_ = tree.template get<int>("FluidSystem.nTemperature");
*
* + We see what we do
* - Reporting whether it was used does not work
* - Overwriting on command line not possible
*/
// set the spatial parameters by reading the DGF grid file
// this->spatialParams().setParams();
// initialize the tables of the fluid system
FluidSystem::init(/*Tmin=*/temperatureLow_,
/*Tmax=*/temperatureHigh_,
/*nT=*/nTemperature_,
/*pmin=*/pressureLow_,
/*pmax=*/pressureHigh_,
/*np=*/nPressure_);
//stating in the console whether mole or mass fractions are used
if(useMoles)
{
std::cout<<"problem uses mole fractions"<<std::endl;
}
else
{
std::cout<<"problem uses mass fractions"<<std::endl;
}
this->timeManager().startNextEpisode(resultsTime1_);
}
/*!
* \brief User defined output after the time integration
*
* Will be called diretly after the time integration.
*/
void postTimeStep()
{
// Calculate storage terms
PrimaryVariables storageL, storageG;
this->model().globalPhaseStorage(storageL, lPhaseIdx);
this->model().globalPhaseStorage(storageG, gPhaseIdx);
// Write mass balance information for rank 0
if (this->gridView().comm().rank() == 0) {
std::cout<<"Storage: liquid=[" << storageL << "]"
<< " gas=[" << storageG << "]\n";
}
}
/*!
* \brief Append all quantities of interest which can be derived
* from the solution of the current time step to the VTK
* writer.
*/
void addOutputVtkFields()
{
typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
// get the number of degrees of freedom
unsigned numDofs = this->model().numDofs();
unsigned numElements = this->gridView().size(0);
//create required scalar fields
ScalarField *Kxx = this->resultWriter().allocateManagedBuffer(numElements);
ScalarField *cellPorosity = this->resultWriter().allocateManagedBuffer(numElements);
ScalarField *boxVolume = this->resultWriter().allocateManagedBuffer(numDofs);
(*boxVolume) = 0;
//Fill the scalar fields with values
ScalarField *rank = this->resultWriter().allocateManagedBuffer(numElements);
FVElementGeometry fvGeometry;
VolumeVariables volVars;
for (const auto& element : elements(this->gridView()))
{
int eIdx = this->elementMapper().index(element);
(*rank)[eIdx] = this->gridView().comm().rank();
fvGeometry.update(this->gridView(), element);
for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
{
int dofIdxGlobal = this->model().dofMapper().subIndex(element, scvIdx, dofCodim);
volVars.update(this->model().curSol()[dofIdxGlobal],
*this,
element,
fvGeometry,
scvIdx,
false);
(*boxVolume)[dofIdxGlobal] += fvGeometry.subContVol[scvIdx].volume;
}
(*Kxx)[eIdx] = this->spatialParams().intrinsicPermeability(element, fvGeometry, /*element data*/ 0);
(*cellPorosity)[eIdx] = this->spatialParams().porosity(element, fvGeometry, /*element data*/ 0);
}
//pass the scalar fields to the vtkwriter
this->resultWriter().attachDofData(*Kxx, "Kxx", false); //element data
this->resultWriter().attachDofData(*cellPorosity, "cellwisePorosity", false); //element data
this->resultWriter().attachDofData(*boxVolume, "boxVolume", isBox);
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief The problem name.
*
* This is used as a prefix for files generated by the simulation.
*/
const std::string name() const
{ return name_; }
void episodeEnd()
{
Scalar sim_time_ = this->timeManager().time();
// Start new episode if episode is over
if(sim_time_ < resultsTime2_){
this->timeManager().startNextEpisode(resultsTime1_);
std::cout<<"Episode index is set to: "<<this->timeManager().episodeIndex()<<std::endl;
std::cout<<"Next episode length is: "<<resultsTime1_<<std::endl;
std::cout<<"Next epidode time is : "<<sim_time_+resultsTime1_<<std::endl;
}else{
this->timeManager().startNextEpisode(resultsTime3_);
std::cout<<"Episode index is set to: "<<this->timeManager().episodeIndex()<<std::endl;
std::cout<<"Next episode length is: "<<resultsTime3_<<std::endl;
std::cout<<"Next epidode time is : "<<sim_time_+resultsTime3_<<std::endl;
}
}
bool shouldWriteOutput() const
{
return
this->timeManager().timeStepIndex() == 0 || this->timeManager().episodeWillBeFinished() ||
this->timeManager().willBeFinished();
}
bool shouldWriteRestartFile() const
{
return
this->timeManager().timeStepIndex() == 0 || this->timeManager().episodeIsFinished() ||
this->timeManager().willBeFinished();
}
/*!
* \brief Returns the temperature within the domain.
*
* \param globalPos The position
*
* This problem assumes a geothermal gradient with
* a surface temperature of 10 degrees Celsius.
*/
Scalar temperatureAtPos(const GlobalPosition &globalPos) const
{
return temperature_(globalPos);
}
/*!
* \brief Returns the source term
*
* \param values Stores the source values for the conservation equations in
* \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
* \param globalPos The global position
*/
void sourceAtPos (PrimaryVariables & values,
const GlobalPosition &globalPos) const
{
values = 0;
Scalar sim_time_ = this->timeManager().time();
if (((globalPos[1] < injecttop_ + eps_ && globalPos[1] > injectbot_ - eps_) && (globalPos[0] < injectright_ + eps_ && globalPos[0] > injectleft_ - eps_))
&& (sim_time_ > timepulsedebut1_ && sim_time_ < timepulsefin1_))
{
values[contiCO2EqIdx] = injectionRate_;
}
}
/*
void Source (PrimaryVariables &values,
const Element &element,
const FVElementGeometry &fvGeometry,
const Intersection &intersection,
int scvIdx,
int boundaryFaceIdx) const
{
GlobalPosition globalPos;
if (isBox)
globalPos = element.geometry().corner(scvIdx);
else
globalPos = intersection.geometry().center();
values = 0;
Scalar sim_time_ = this->timeManager().time();
if (((globalPos[1] < injecttop_ + eps_ && globalPos[1] > injectbot_ - eps_) && (globalPos[0] < injectright_ + eps_ && globalPos[0] > injectleft_ - eps_))
&& (sim_time_ > timepulsedebut1_ && sim_time_ < timepulsefin1_))
{
values[contiCO2EqIdx] = injectionRate_;
}
}
*/
// \}
/*!
* \name Boundary conditions
*/
// \{
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
* \param values The boundary types for the conservation equations
* \param vertex The vertex for which the boundary type is set
*/
/* void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
{
intersectionToVertexBC_.boundaryTypes(values, vertex);
}
*/
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
* \param values The boundary types for the conservation equations
* \param intersection specifies the intersection at which boundary
* condition is to set
*/
/*
void boundaryTypes(BoundaryTypes &values, const Intersection &intersection) const
{
#if DUNE_VERSION_NEWER(DUNE_ALUGRID, 2, 6)
int boundaryId = intersection.impl().boundaryId();
#else
int boundaryId = intersection.boundaryId();
#endif
if (boundaryId < 1 || boundaryId > 4)
{
std::cout<<"invalid boundaryId: "<<boundaryId<<std::endl;
DUNE_THROW(Dune::InvalidStateException, "Invalid " << boundaryId);
}
if (boundaryId == dirichletBoundary_)
values.setAllDirichlet();
else
values.setAllNeumann();
}
*/
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment
*
* \param values Stores the value of the boundary type
* \param globalPos The global position
*/
void boundaryTypesAtPos(BoundaryTypes &values,
const GlobalPosition &globalPos) const
{
/* if (globalPos[1] < injecttop_ + eps_ && globalPos[1] > injectbot_ - eps_ && globalPos[0] < eps_)
values.setAllNeumann();
*/
if (globalPos[0] < eps_ || globalPos[0] > 1.5 - eps_)
values.setAllDirichlet();
else
values.setAllNeumann();
}
/*!
* \brief Evaluates the boundary conditions for a Dirichlet
* boundary segment
*
* \param values Stores the Dirichlet values for the conservation equations in
* \f$ [ \textnormal{unit of primary variable} ] \f$
* \param globalPos The global position
*/
void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
initial_(values, globalPos);
}
/*!
* \brief Evaluate the boundary conditions for a Neumann
* boundary segment.
*
* \param values Stores the Neumann values for the conservation equations in
* \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
* \param element The finite element
* \param fvGeometry The finite volume geometry of the element
* \param intersection The intersection between element and boundary
* \param scvIdx The local index of the sub-control volume
* \param boundaryFaceIdx The index of the boundary face
*
* The \a values store the mass flux of each phase normal to the boundary.
* Negative values indicate an inflow.
*
* Depending on whether useMoles is set on true or false, the flux has to be given either in
* kg/(m^2*s) or mole/(m^2*s) in the input file!! Convertion with molar mass obtained from fluid system FluidSystem::molarMass(nCompIdx)
*/
void neumann(PrimaryVariables &values,
const Element &element,
const FVElementGeometry &fvGeometry,
const Intersection &intersection,
int scvIdx,
int boundaryFaceIdx) const
{
values = 0;
/* GlobalPosition globalPos;
if (isBox)
globalPos = element.geometry().corner(scvIdx);
else
globalPos = intersection.geometry().center();
Scalar sim_time_ = this->timeManager().time();
if (globalPos[1] < injecttop_ + eps_ && globalPos[1] > injectbot_ - eps_ && globalPos[0] < eps_ &&
((sim_time_ > timepulsedebut1_ && sim_time_ < timepulsefin1_) || (sim_time_ > timepulsedebut2_ && sim_time_ < timepulsefin2_)))
{
values[contiCO2EqIdx] = -1*injectionRate_; // kg/(s*m)
// std::cout<<"Injection rate will be : "<<-1*injectionRate_<<std::endl;
} */
}
/*
void neumannAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
values = 0;
Scalar sim_time_ = this->timeManager().time();
if (globalPos[1] < injecttop_ + eps_ && globalPos[1] > injectbot_ - eps_ && globalPos[0] < eps_ &&
((sim_time_ > timepulsedebut1_ && sim_time_ < timepulsefin1_) || (sim_time_ > timepulsedebut2_ && sim_time_ < timepulsefin2_)))
{
values[contiCO2EqIdx] = -1*injectionRate_; // kg/(s*m)
std::cout<<"Injection rate will be : "<<-1*injectionRate_<<std::endl;
}
}
*/
// \}
/*!
* \name Volume terms
*/
// \{
/*!
* \brief Evaluates the initial values for a control volume
*
* \param values Stores the initial values for the conservation equations in
* \f$ [ \textnormal{unit of primary variables} ] \f$
* \param globalPos The global position
*/
void initialAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
initial_(values, globalPos);
}
/*!
* \brief Return the initial phase state inside a control volume.
*
* \param globalPos The global position
*/
int initialPhasePresenceAtPos(const GlobalPosition &globalPos) const
{ return Indices::wPhaseOnly; }
// \}
private:
/*!
* \brief Evaluates the initial values for a control volume
*
* The internal method for the initial condition
*
* \param values Stores the initial values for the conservation equations in
* \f$ [ \textnormal{unit of primary variables} ] \f$
* \param globalPos The global position
*/
void initial_(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
Scalar temp = temperature_(globalPos);
Scalar densityW = FluidSystem::Brine::liquidDensity(temp, 1e5);
Scalar pl = 1e5 + densityW*this->gravity()[dimWorld-1]*(globalPos[dimWorld-1] - depthBOR_);
if (globalPos[0] < eps_ )
pl = pl + hzdp_;
Scalar moleFracLiquidCO2 = 0.00;
Scalar moleFracLiquidBrine = 1.0 - moleFracLiquidCO2;
Scalar meanM =
FluidSystem::molarMass(BrineIdx)*moleFracLiquidBrine +
FluidSystem::molarMass(CO2Idx)*moleFracLiquidCO2;
if(useMoles) // mole-fraction formulation
{
values[Indices::switchIdx] = moleFracLiquidCO2;
}
else // mass-fraction formulation
{
Scalar massFracLiquidCO2 = moleFracLiquidCO2*FluidSystem::molarMass(CO2Idx)/meanM;
values[Indices::switchIdx] = massFracLiquidCO2;
}
values[Indices::pressureIdx] = pl;
}
Scalar temperature_(const GlobalPosition globalPos) const
{
Scalar T = celltemp_ ;//+ (depthBOR_ - globalPos[dimWorld-1])*0.03;
return T;
}
Scalar depthBOR_;
Scalar injectionRate_;
Scalar timepulsefin1_;
Scalar timepulsedebut1_;
Scalar timepulsefin2_;
Scalar timepulsedebut2_;
Scalar injecttop_;
Scalar injectbot_;
Scalar injectright_;
Scalar injectleft_;
Scalar celltemp_;
Scalar hzdp_;
Scalar resultsTime1_;
Scalar resultsTime2_;
Scalar resultsTime3_;
int nTemperature_;
int nPressure_;
std::string name_ ;
Scalar pressureLow_, pressureHigh_;
Scalar temperatureLow_, temperatureHigh_;
/*
int injectionTop_;
int injectionBottom_;
int dirichletBoundary_;
int noFlowBoundary_;
*/
static constexpr Scalar eps_ = 1.5e-7;
// const IntersectionToVertexBC<TypeTag> intersectionToVertexBC_;
};
} //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
*
* \brief Definition of the spatial parameters for the heterogeneous
* problem which uses the non-isothermal or isothermal CO2
* fully implicit model.
*/
#ifndef DUMUX_HETEROGENEOUS_SPATIAL_PARAMS_HH
#define DUMUX_HETEROGENEOUS_SPATIAL_PARAMS_HH
#include <dumux/material/spatialparams/implicit.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include <dumux/porousmediumflow/co2/implicit/model.hh>
#include <random>
namespace Dumux {
//forward declaration
template<class TypeTag>
class HeterogeneousSpatialParams;
namespace Properties
{
// The spatial parameters TypeTag
NEW_TYPE_TAG(HeterogeneousSpatialParams);
// Set the spatial parameters
SET_TYPE_PROP(HeterogeneousSpatialParams, SpatialParams, HeterogeneousSpatialParams<TypeTag>);
// Set the material Law
SET_PROP(HeterogeneousSpatialParams, MaterialLaw)
{
private:
// define the material law which is parameterized by effective
// saturations
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef RegularizedBrooksCorey<Scalar> EffMaterialLaw;
public:
// define the material law parameterized by absolute saturations
typedef EffToAbsLaw<EffMaterialLaw> type;
};
}
/*!
* \ingroup CO2Model
* \ingroup ImplicitTestProblems
* \brief Definition of the spatial parameters for the heterogeneous
* problem which uses the non-isothermal or isothermal CO2
* fully implicit model.
*/
template<class TypeTag>
class HeterogeneousSpatialParams : public ImplicitSpatialParams<TypeTag>
{
typedef ImplicitSpatialParams<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
// enum { dim=GridView::dimension };
enum {
dim=GridView::dimension,
dimWorld=GridView::dimensionworld
};
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GridView::template Codim<0>::Entity Element;
public:
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename MaterialLaw::Params MaterialLawParams;
// using PermeabilityType = Scalar;
HeterogeneousSpatialParams(std::shared_ptr<const FVElementGeometry> fvGeometry)
: ParentType(fvGeometry), K_(fvGeometry->gridView().size(0), 0.0)
/*!
* \brief The constructor
*
* \param gridView The grid view
*/
// HeterogeneousSpatialParams(const GridView &gridView)
// : ParentType(gridView), gridView_(gridView)
{
// /*
// * Layer Index Setup:
// */
// barrierTop_ = 1;
// barrierMiddle_ = 2;
// reservoir_ = 3;
// heat conductivity of granite
lambdaSolid_ = 2.8;
//Set the permeability for the layers
// barrierTopK_ = 1e-17; //sqm
coucheK_ = 2.5e-10; //sqm
reservoirK_ = 5e-10; //sqm
std::mt19937 rand(0);
std::lognormal_distribution<Scalar> K(std::log(reservoirK_), std::log(reservoirK_)*0.1);
std::lognormal_distribution<Scalar> KLens(std::log(coucheK_), std::log(coucheK_)*0.1);
for (const auto& element : elements(fvGeometry->gridView()))
{
const auto eIdx = fvGeometry->elementMapper().index(element);
const auto globalPos = element.geometry().center();
K_[eIdx] = barrierMiddle_(globalPos) || barrierTop_(globalPos) ? KLens(rand) : K(rand);
}
//Set the effective porosity of the layers
// barrierTopPorosity_ = 0.001;
couchePorosity_ = 0.36;
reservoirPorosity_ = 0.36;
// Same material parameters for every layer
// materialParams_.setSwr(0.2);
// materialParams_.setSnr(0.05);
// materialParams_.setLambda(2.0);
// materialParams_.setPe(1e4);
// parameters for the Brooks-Corey law
coucheParams_.setPe(849.55);
reservoirParams_.setPe(531.7);
coucheParams_.setLambda(5.57);
reservoirParams_.setLambda(3.94);
coucheParams_.setSwr(0.0);
coucheParams_.setSnr(0.0);
reservoirParams_.setSwr(0.0);
reservoirParams_.setSnr(0.0);
}
~HeterogeneousSpatialParams()
{}
/*!
* \brief Reads layer information from the grid
*
*/
/* void setParams()
{
int numElements = gridView_.size(0);
paramIdx_.resize(numElements);
for (const auto& element : elements(gridView_))
{
int eIdx = gridView_.indexSet().index(element);
int param = GridCreator::parameters(element)[0];
paramIdx_[eIdx] = param;
}
}
*/
/*!
* \brief Returns the scalar intrinsic permeability \f$[m^2]\f$
*
* \param element The finite element
* \param fvGeometry The finite volume geometry of the element
* \param scvIdx The local index of the sub-control volume
*/
const Scalar intrinsicPermeability(const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
{
const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
if (barrierMiddle_(globalPos) || barrierTop_(globalPos)) //modif par kenza barrierMiddle_ et barrierTop_(globalPos)
return coucheK_;
return reservoirK_;
}
Scalar porosity(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
if (barrierMiddle_(globalPos) || barrierTop_(globalPos))
return couchePorosity_;
return reservoirPorosity_;
}
/* {
//Get the global index of the element
int eIdx = gridView_.indexSet().index(element);
// if (paramIdx_[eIdx] == barrierTop_)
// return barrierTopK_;
else if (paramIdx_[eIdx] == barrierMiddle_)
return barrierMiddleK_;
else
return reservoirK_;
}
*/
/*!
* \brief Returns the porosity \f$[-]\f$
*
* \param element The finite element
* \param fvGeometry The finite volume geometry of the element
* \param scvIdx The local index of the sub-control volume
*/
/* Scalar porosity(const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
{
//Get the global index of the element
int eIdx = gridView_.indexSet().index(element);
// if (paramIdx_[eIdx] == barrierTop_)
// return barrierTopPorosity_;
else if (paramIdx_[eIdx] == barrierMiddle_)
return barrierMiddlePorosity_;
else
return reservoirPorosity_;
}
*/
/*!
* \brief Returns the parameter object for the Brooks-Corey material law
*
* \param element The finite element
* \param fvGeometry The finite volume geometry of the element
* \param scvIdx The local index of the sub-control volume
*/
const MaterialLawParams& materialLawParams(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
if (barrierMiddle_(globalPos) || barrierTop_(globalPos))
return coucheParams_;
return reservoirParams_;
}
/*
const MaterialLawParams& materialLawParams(const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
{
return materialParams_;
}
*/
/*!
* \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
*
* This is only required for non-isothermal models.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume
*/
Scalar solidHeatCapacity(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return 790; // specific heat capacity of granite [J / (kg K)]
}
/*!
* \brief Returns the mass density \f$[kg / m^3]\f$ of the rock matrix.
*
* This is only required for non-isothermal models.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume
*/
Scalar solidDensity(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return 2700; // density of granite [kg/m^3]
}
/*!
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid
*
* This is only required for non-isothermal models.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry of the element
* \param scvIdx The local index of the sub-control volume
*/
Scalar solidThermalConductivity(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return lambdaSolid_;
}
const std::vector<Scalar>& getKField() const
{ return K_; }
private:
//bool isFineMaterial_(const GlobalPosition &globalPos) const
//{ return globalPos[dimWorld-1] > layerBottom_; }
bool barrierTop_(const GlobalPosition &globalPos) const
{
// if (globalPos[1] < 9 + eps_ || globalPos[1] > 18 - eps_ || globalPos[0] < 40 + eps_)
if (globalPos[1] < globalPos[0] * 0.02 + 1.35) // couche fine avec une epaisseure non uniforme
return false;
return true;
}
bool barrierMiddle_(const GlobalPosition &globalPos) const
{
if (globalPos[1] < 0.7 + eps_ || globalPos[1] > 0.78 - eps_ || globalPos[0] > 0.8 - eps_)
// if (globalPos[1] < globalPos[0] * 0.05 + 25) // couche fine avec une epaisseure non uniforme
return false;
return true;
}
private:
// int barrierTop_;
int couche_;
int reservoir_;
// Scalar barrierTopPorosity_;
Scalar couchePorosity_;
Scalar reservoirPorosity_;
// Scalar barrierTopK_;
Scalar coucheK_;
Scalar reservoirK_;
Scalar lambdaSolid_;
MaterialLawParams coucheParams_;
MaterialLawParams reservoirParams_;
std::vector<Scalar> K_;
static constexpr Scalar eps_ = 1.5e-7;
const GridView gridView_;
std::vector<int> paramIdx_;
};
}
#endif
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
