Dear Ali,
attached to this mail are two files (problem and model file) where
diffussive and advective fluxes over a certain area or layer are written
into the command line output at the end of every timestep. You can use
these files as an orientation.
However, this output is only useful if the layer is situated
perpendicular to the coordinate axis and cuts the axis at a certain
coordinate. If you need something else you have to adjust the function
"calculateFluxAcrossLayer(...)" in eco2model.hh to your needs.
To use this function you can derive a new model file from the
co2model.hh and include the function "calculateFluxAcrossLayer(...)"
there and then call this function from the function "postTimeStep()" in
your problem file (as done in the attached eco2problem.hh).
Alternatively you can also adjust the attached files to your problem
demands.
best regards
Alex
On 05/13/2013 06:15 AM, Ali Nowamooz wrote:
Dear Dumux,
I worked really hard on my problem but i didn't find a solution. I really need
your help.
I am working on HeterogeneousProblem (co2 injection) and I wanted to know if I
have any possibility to write advective and diffusive fluxes over faces of
sub-control-volume on an external file (for example a text file).
I think both these fluxes are calculated inside 2p2clocalresidual.hh
(computeAdvectiveFlux and computeDispersiveFlux).
I tried to call them as PrimaryVariables but i couldn't go further. My problem
is obviously due to my lack of knowledge of c++.
Can you please help me and give me some information?
best regards
Ali
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
--
*******************************************************
!!!! CMWR 2014: 9th - 13th June 2014 in Stuttgart !!!!
Please visit www.cmwr14.de
*******************************************************
Alexander Kissinger
Institut für Wasser- und Umweltsystemmodellierung
Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
Pfaffenwaldring 61
D-70569 Stuttgart
Telefon: +49 (0) 711 685-64729
E-Mail: [email protected]
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2012 by *
* Institute for Modelling Hydraulic and Environmental Systems *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
* *
* 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 Adaption of the BOX scheme to the two-phase two-component flow model without constraint solver.
*/
#ifndef DUMUX_CO2_MODEL_HH
#define DUMUX_CO2_MODEL_HH
#include <dumux/implicit/2p2c/2p2cmodel.hh>
namespace Dumux {
/*!
* \ingroup CO2Model
* \brief Adaption of the BOX scheme to the non-isothermal two-phase two-component flow model.
* See TwoPTwoCModel for reference to the equations used.
* The CO2 model is derived from the 2p2c model. In the CO2 model the phase switch criterion
* is different from the 2p2c model.
* The phase switch occurs when the equilibrium concentration
* of a component in a phase is exceeded, instead of the sum of the components in the virtual phase
* (the phase which is not present) being greater that unity as done in the 2p2c model.
* The CO2VolumeVariables do not use a constraint solver for calculating the mole fractions as is the
* case in the 2p2c model. Instead mole fractions are calculated in the FluidSystem with a given
* temperature, pressurem and salinity.
*
*/
template<class TypeTag>
class CO2Model: public TwoPTwoCModel<TypeTag> {
typedef TwoPTwoCModel<TypeTag> ParentType;typedef typename GET_PROP_TYPE(TypeTag, BaseModel) BaseType;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
enum {
numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
numComponents = GET_PROP_VALUE(TypeTag, NumComponents)
};
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum {
pressureIdx = Indices::pressureIdx,
switchIdx = Indices::switchIdx,
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
wCompIdx = Indices::wCompIdx,
nCompIdx = Indices::nCompIdx,
wPhaseOnly = Indices::wPhaseOnly,
nPhaseOnly = Indices::nPhaseOnly,
bothPhases = Indices::bothPhases,
pwSn = TwoPTwoCFormulation::pwSn,
pnSw = TwoPTwoCFormulation::pnSw,
formulation = GET_PROP_VALUE(TypeTag, Formulation),
contiWEqIdx = Indices::contiWEqIdx,
contiNEqIdx = Indices::contiNEqIdx
};
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename Element::Geometry Geometry;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
enum {
dim = GridView::dimension, dimWorld = GridView::dimensionworld
};
enum {
numEq = GET_PROP_VALUE(TypeTag, NumEq)
};
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GridView::ctype CoordScalar;
typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> LocalFiniteElementCache;
typedef typename LocalFiniteElementCache::FiniteElementType LocalFiniteElement;
typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<Scalar, dimWorld> Vector;
enum {
isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox)
};
enum {
dofCodim = isBox ? dim : 0
};
public:
/*!
* \brief Update the static data of all vertices in the grid.
*
* \param curGlobalSol The current global solution
* \param oldGlobalSol The previous global solution
*/
void updateStaticData(SolutionVector &curGlobalSol,
const SolutionVector &oldGlobalSol) {
bool wasSwitched = false;
for (unsigned i = 0; i < ParentType::staticDat_.size(); ++i)
ParentType::staticDat_[i].visited = false;
FVElementGeometry fvGeometry;
static VolumeVariables volVars;
ElementIterator elemIt = this->gridView_().template begin<0>();
const ElementIterator &elemEndIt = this->gridView_().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) {
fvGeometry.update(this->gridView_(), *elemIt);
for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx) {
int globalIdx = this->dofMapper().map(*elemIt, scvIdx,
dofCodim);
if (ParentType::staticDat_[globalIdx].visited)
continue;
ParentType::staticDat_[globalIdx].visited = true;
volVars.update(curGlobalSol[globalIdx], this->problem_(),
*elemIt, fvGeometry, scvIdx, false);
const GlobalPosition &globalPos = elemIt->geometry().corner(
scvIdx);
if (primaryVarSwitch_(curGlobalSol, volVars, globalIdx,
globalPos)) {
this->jacobianAssembler().markDofRed(globalIdx);
wasSwitched = true;
}
}
}
// make sure that if there was a variable switch in an
// other partition we will also set the switch flag
// for our partition.
if (this->gridView_().comm().size() > 1)
wasSwitched = this->gridView_().comm().max(wasSwitched);
ParentType::setSwitched_(wasSwitched);
}
/*!
* \brief Calculate the fluxes across a certain layer in the domain.
* The layer is situated perpendicular to the coordinate axis "coord" and cuts
* the axis at the value "coordVal".
*
* \param globalSol The global solution vector
* \param flux A vector to store the flux
* \param axis The dimension, perpendicular to which the layer is situated
* \param coordVal The (Scalar) coordinate on the axis, at which the layer is situated
*/
void calculateFluxAcrossLayer(PrimaryVariables &flux, int coord,
Scalar coordVal) {
Scalar massUpwindWeight = GET_PARAM_FROM_GROUP(TypeTag, Scalar,
Implicit, MassUpwindWeight);
ElementVolumeVariables elemVolVars;
FVElementGeometry fvElemGeom;
ElementIterator elemIt = this->gridView_().template begin<0>();
const ElementIterator &endit = this->gridView_().template end<0>();
GlobalPosition globalI, globalJ;
PrimaryVariables tmpFlux(0.0);
int sign;
// Loop over elements
for (; elemIt != endit; ++elemIt) {
if (elemIt->partitionType() == Dune::InteriorEntity) {
fvElemGeom.update(this->gridView_(), *elemIt);
elemVolVars.update(this->problem_(), *elemIt, fvElemGeom,
false /* oldSol? */);
if (elemIt->partitionType() != Dune::InteriorEntity)
continue;
for (int faceId = 0; faceId < fvElemGeom.numEdges; faceId++) {
int idxI = fvElemGeom.subContVolFace[faceId].i;
int idxJ = fvElemGeom.subContVolFace[faceId].j;
int flagI, flagJ;
globalI = fvElemGeom.subContVol[idxI].global;
globalJ = fvElemGeom.subContVol[idxJ].global;
// 2D case: give y or x value of the line over which flux is to be
// calculated.
// up to now only flux calculation to lines or planes (3D) parallel to
// x, y and z axis possible
// Flux across plane with z = 80 numEq
if (globalI[coord] < coordVal)
flagI = 1;
else
flagI = -1;
if (globalJ[coord] < coordVal)
flagJ = 1;
else
flagJ = -1;
if (flagI == flagJ) {
sign = 0;
} else {
if (flagI > 0)
sign = -1;
else
sign = 1;
}
// get variables
if (flagI != flagJ) {
FluxVariables vars(this->problem_(), *elemIt,
fvElemGeom, faceId, elemVolVars);
tmpFlux = 0;
//asImp_()->computeAdvectiveFlux(flux, vars);
//asImp_()->computeDiffusiveFlux(flux, vars);
////////
// advective fluxes of all components in all phases
////////
Vector tmpVec;
for (int phaseIdx = 0; phaseIdx < numPhases;
++phaseIdx) {
// calculate the flux in the normal direction of the
// current sub control volume face:
//
// v = - (K grad p) * n
//
// (the minus comes from the Darcy law which states that
// the flux is from high to low pressure potentials.)
Scalar normalFlux = -vars.kGradPNormal(phaseIdx);
// vars.intrinsicPermeability().mv(vars.potentialGrad(phaseIdx), tmpVec);
// Scalar normalFlux = - (tmpVec*vars.face().normal);
// data attached to upstream and the downstream vertices
// of the current phase
const VolumeVariables &up =
elemVolVars[vars.upstreamIdx(normalFlux)];
const VolumeVariables &dn =
elemVolVars[vars.downstreamIdx(normalFlux)];
// add advective flux of current component in current
// phase
int eqIdx =
(phaseIdx == wPhaseIdx) ?
contiWEqIdx : contiNEqIdx;
tmpFlux[eqIdx] += normalFlux
* ((massUpwindWeight) * up.density(phaseIdx)
* up.mobility(phaseIdx)
+ (1 - massUpwindWeight)
* dn.density(phaseIdx)
* dn.mobility(phaseIdx));
}
// add diffusive flux of gas component in liquid phase
tmpFlux += 0.0;
// the face normal points into the outward direction, so we
// have to multiply all fluxes with -1
tmpFlux *= -1;
// this->localResidual().computeFlux(tmpFlux, faceId);
tmpFlux *= sign;
flux += tmpFlux;
}
}
}
}
//If parallel sum of fluxes over all processors
flux = this->problem_().gridView().comm().sum(flux);
}
template<class MultiWriter>
void addOutputVtkFields(const SolutionVector &sol, MultiWriter &writer) {
unsigned numElements = this->gridView_().size(0);
ParentType::addOutputVtkFields(sol, writer);
typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
if (this->problem_().shouldWriteDataFile() == true) {
const std::string filename("PressSatData");
const std::string suffix(".dat");
std::ostringstream sstream, sstream1;
sstream
<< (this->problem_().timeManager().time()
+ this->problem_().timeManager().timeStepSize());
sstream1 << this->gridView().comm().rank();
std::string time = sstream.str();
std::string rank = sstream1.str();
std::string tmpname, tmpname2, tmpname1;
tmpname = filename;
tmpname += "_rank";
tmpname += rank;
tmpname += "_time";
tmpname += time;
tmpname += suffix;
std::ofstream data(tmpname.c_str());
data << "#x y z pw Sn \n";
int cellIndex = 0;
ElementIterator elemIt = this->gridView_().template begin<0>();
ElementIterator elemEndIt = this->gridView_().template end<0>();
FVElementGeometry fvGeometry;
VolumeVariables volVars;
for (; elemIt != elemEndIt; ++elemIt) {
if (elemIt->partitionType() == Dune::InteriorEntity) {
fvGeometry.update(this->gridView_(), *elemIt);
const Geometry& geometry = (*elemIt).geometry();
Dune::GeometryType gt = geometry.type();
const LocalFiniteElement &localFiniteElement = feCache_.get(
gt);
const Dune::FieldVector<Scalar, dim> &localPos =
ReferenceElements::general(
elemIt->geometry().type()).position(0, 0);
GlobalPosition globalPosCell = elemIt->geometry().center();
//data<<globalPos;
typedef Dune::FieldVector<Scalar, 1> ShapeValue;
std::vector < ShapeValue > shapeVal;
localFiniteElement.localBasis().evaluateFunction(localPos,
shapeVal);
int numVerts = elemIt->template count<dim>();
fvGeometry.update(this->gridView_(), *elemIt);
int cellIndex = this->elementMapper().map(*elemIt);
Scalar pW, Sn;
Scalar cellPw = 0;
Scalar cellSn = 0;
bool identifierSet = false;
std::string identifier("interior");
for (int i = 0; i < numVerts; ++i) {
GlobalPosition globalPosVertex =
elemIt->geometry().corner(i);
int globalIdx = this->vertexMapper().map(*elemIt, i,
dim);
volVars.update(sol[globalIdx], this->problem_(),
*elemIt, fvGeometry, i, false);
pW = volVars.pressure(wPhaseIdx);
Sn = volVars.saturation(nPhaseIdx);
pW *= shapeVal[0];
Sn *= shapeVal[0];
cellPw += pW;
cellSn += Sn;
if (this->problem_().shouldWriteDataFile() == true) {
if (identifierSet == false
&& this->problem_().inLeakageFeature(
globalPosVertex)) {
identifier = "leakage feature";
identifierSet = true;
}
if (identifierSet == false
&& this->problem_().atWell(
globalPosVertex)) {
identifier = "well";
identifierSet = true;
}
if (identifierSet == false
&& this->problem_().onLeftBoundary(
globalPosVertex)) {
identifier = "left boundary";
identifierSet = true;
}
if (identifierSet == false
&& this->problem_().onRightBoundary(
globalPosVertex)) {
identifier = "right boundary";
identifierSet = true;
}
if (identifierSet == false
&& this->problem_().onLowerBoundary(
globalPosVertex)) {
identifier = "lower boundary";
identifierSet = true;
}
if (identifierSet == false
&& this->problem_().onUpperBoundary(
globalPosVertex)) {
identifier = "upper boundary";
identifierSet = true;
}
if (identifierSet == false
&& this->problem_().onFrontBoundary(
globalPosVertex)) {
identifier = "front boundary";
identifierSet = true;
}
if (identifierSet == false
&& this->problem_().onBackBoundary(
globalPosVertex)) {
identifier = "back boundary";
identifierSet = true;
}
}
data << globalPosCell << " " << cellPw << " " << cellSn
<< " " << identifier << "\n";
}
}
}
data.close();
}
}
protected:
/*!
* \brief Set the old phase of all verts state to the current one.
*/
bool primaryVarSwitch_(SolutionVector &globalSol,
const VolumeVariables &volVars, int globalIdx,
const GlobalPosition &globalPos) {
typename FluidSystem::ParameterCache paramCache;
// evaluate primary variable switch
bool wouldSwitch = false;
int phasePresence = ParentType::staticDat_[globalIdx].phasePresence;
int newPhasePresence = phasePresence;
// check if a primary var switch is necessary
if (phasePresence == nPhaseOnly) {
Scalar xnw = volVars.fluidState().moleFraction(nPhaseIdx, wCompIdx);
Scalar xnwMax = FluidSystem::equilibriumMoleFraction(
volVars.fluidState(), paramCache, nPhaseIdx);
if (xnw > xnwMax)
wouldSwitch = true;
if (ParentType::staticDat_[globalIdx].wasSwitched)
xnwMax *= 1.02;
//If mole fraction is higher than the equilibrium mole fraction make a phase switch
if (xnw > xnwMax) {
// wetting phase appears
std::cout << "wetting phase appears at vertex " << globalIdx
<< ", coordinates: " << globalPos << ", xnw > xnwMax: "
<< xnw << " > " << xnwMax << std::endl;
newPhasePresence = bothPhases;
if (formulation == pnSw)
globalSol[globalIdx][switchIdx] = 0.0;
else if (formulation == pwSn)
globalSol[globalIdx][switchIdx] = 1.0;
}
} else if (phasePresence == wPhaseOnly) {
Scalar xwn = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
Scalar xwnMax = FluidSystem::equilibriumMoleFraction(
volVars.fluidState(), paramCache, wPhaseIdx);
//If mole fraction is higher than the equilibrium mole fraction make a phase switch
if (xwn > xwnMax)
wouldSwitch = true;
if (ParentType::staticDat_[globalIdx].wasSwitched)
xwnMax *= 1.02;
if (xwn > xwnMax) {
// non-wetting phase appears
std::cout << "non-wetting phase appears at vertex " << globalIdx
<< ", coordinates: " << globalPos << ", xwn > xwnMax: "
<< xwn << " > " << xwnMax << std::endl;
newPhasePresence = bothPhases;
if (formulation == pnSw)
globalSol[globalIdx][switchIdx] = 0.999;
else if (formulation == pwSn)
globalSol[globalIdx][switchIdx] = 0.001;
}
} else if (phasePresence == bothPhases) {
Scalar Smin = 0.0;
if (ParentType::staticDat_[globalIdx].wasSwitched)
Smin = -0.01;
if (volVars.saturation(nPhaseIdx) <= Smin) {
wouldSwitch = true;
// nonwetting phase disappears
std::cout << "Nonwetting phase disappears at vertex "
<< globalIdx << ", coordinates: " << globalPos
<< ", Sn: " << volVars.saturation(nPhaseIdx)
<< std::endl;
newPhasePresence = wPhaseOnly;
globalSol[globalIdx][switchIdx] =
volVars.fluidState().massFraction(wPhaseIdx, nCompIdx);
} else if (volVars.saturation(wPhaseIdx) <= Smin) {
wouldSwitch = true;
// wetting phase disappears
std::cout << "Wetting phase disappears at vertex " << globalIdx
<< ", coordinates: " << globalPos << ", Sw: "
<< volVars.saturation(wPhaseIdx) << std::endl;
newPhasePresence = nPhaseOnly;
globalSol[globalIdx][switchIdx] =
volVars.fluidState().massFraction(nPhaseIdx, wCompIdx);
}
}
ParentType::staticDat_[globalIdx].phasePresence = newPhasePresence;
ParentType::staticDat_[globalIdx].wasSwitched = wouldSwitch;
return phasePresence != newPhasePresence;
}
private:
const LocalFiniteElementCache feCache_;
}
;
}
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2012 by *
* Institute for Modelling Hydraulic and Environmental Systems *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
* *
* 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 an injection problem for the generic leakage scnearios
* of the ECO2 project
*
* Up to now five different scenarios have been developed:
* a chimney scenario
* a fault scenario
* a leaky well scenario
* a catastrophic blowout scneario
*
* This file contains the boundary conditions of the chimney scenario
* the boundary conditions for the three remaining scenarios are given in comments
*/
#ifndef DUMUX_ECO2_PROBLEM_HH
#define DUMUX_ECO2_PROBLEM_HH
#include <dune/grid/io/file/dgfparser/dgfs.hh>
#include "eco2model.hh"
#include <dumux/implicit/co2/co2volumevariables.hh>
#include <dumux/material/fluidsystems/brineco2fluidsystem.hh>
#include <dumux/implicit/common/implicitporousmediaproblem.hh>
#include <dumux/implicit/box/intersectiontovertexbc.hh>
#include "eco2spatialparams.hh"
#include "eco2co2tables.hh"
namespace Dumux
{
template <class TypeTag>
class ECO2Problem;
namespace Properties
{
NEW_TYPE_TAG(ECO2Problem, INHERITS_FROM(BoxTwoPTwoC, ECO2SpatialParams));
// Set the grid type
SET_PROP(ECO2Problem, Grid)
{
typedef Dune::ALUSimplexGrid<3,3> type;
};
// Set the problem property
SET_PROP(ECO2Problem, Problem)
{
typedef Dumux::ECO2Problem<TTAG(ECO2Problem)> type;
};
// Set fluid configuration
SET_PROP(ECO2Problem, FluidSystem)
{
typedef Dumux::BrineCO2FluidSystem<TypeTag> type;
};
// Set the CO2 table to be used; in this case not the the default table
SET_TYPE_PROP(ECO2Problem, CO2Table, Dumux::ECO2CO2Tables::CO2Tables);
// Set the salinity mass fraction of the brine in the reservoir
SET_SCALAR_PROP(ECO2Problem, ProblemSalinity, 1e-1);
//! the CO2 Model and VolumeVariables properties
SET_TYPE_PROP(ECO2Problem, Model, CO2Model<TypeTag>);
SET_TYPE_PROP(ECO2Problem, VolumeVariables, CO2VolumeVariables<TypeTag>);
// Enable gravity
SET_BOOL_PROP(ECO2Problem, ProblemEnableGravity, true);
SET_BOOL_PROP(ECO2Problem, ImplicitEnableJacobianRecycling, false);
SET_BOOL_PROP(ECO2Problem, VtkAddVelocity, true);
}
/*!
* \ingroup TwoPTwoCModel
* \ingroup ECO2Problem
* \brief Definition of an injection problem for the generic leakage scnearios of the ECO2 project
*
* To run the simulation execute the following line in shell:
* <tt>./eco2_scenarios -parameterFile ./eco2_scenarios.input</tt>
*/
template <class TypeTag = TTAG(ECO2Problem)>
class ECO2Problem : public ImplicitPorousMediaProblem<TypeTag>
{
typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
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 {
// 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,
BrineIdx = FluidSystem::BrineIdx,
CO2Idx = FluidSystem::CO2Idx,
conti0EqIdx = Indices::conti0EqIdx,
contiCO2EqIdx = conti0EqIdx + CO2Idx
};
typedef Dune::FieldVector<Scalar, dim> DimVector;
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, PTAG(CO2Table)) CO2Table;
typedef Dumux::CO2<Scalar, CO2Table> CO2;
public:
/*!
* \brief The constructor
*
* \param timeManager The time manager
* \param gridView The grid view
*/
ECO2Problem(TimeManager &timeManager,
const GridView &gridView)
: ParentType(timeManager, GridCreator::grid().leafView())
{
try
{
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, InitialConditions.DepthBOR);
depthSeafloor_ = GET_RUNTIME_PARAM(TypeTag, Scalar, InitialConditions.DepthSeafloor);
name_ = GET_RUNTIME_PARAM(TypeTag, std::string, SimulationControl.Name);
chimneyScenario_ = GET_RUNTIME_PARAM(TypeTag, bool, Problem.ChimneyScenario);
catastrophicBlowoutScenario_ = GET_RUNTIME_PARAM(TypeTag, bool, Problem.CatastrophicBlowoutScenario);
faultScenario_ = GET_RUNTIME_PARAM(TypeTag, bool, Problem.FaultScenario);
leakyWellScenario_ = GET_RUNTIME_PARAM(TypeTag, bool, Problem.LeakyWellScenario);
}
catch (Dumux::ParameterException &e) {
std::cerr << e << ". Abort!\n";
exit(1) ;
}
catch (...) {
std::cerr << "Unknown exception thrown!\n";
exit(1);
}
eps_ = 1e-6;
// initialize the tables of the fluid system
//FluidSystem::init();
FluidSystem::init(/*Tmin=*/temperatureLow_,
/*Tmax=*/temperatureHigh_,
/*nT=*/nTemperature_,
/*pmin=*/pressureLow_,
/*pmax=*/pressureHigh_,
/*np=*/nPressure_);
this->timeManager().startNextEpisode(3.1536e7);
}
/*!
* \brief Returns true if the current solution should be written to
* disk (i.e. as a VTK file)
*
* The default behaviour is to write out every the solution for
* very time step. This file is intented to be overwritten by the
* implementation.
*/
bool shouldWriteOutput() const
{
// output files are written according to episode lengths
return
this->timeManager().timeStepIndex() == 0 ||
this->timeManager().episodeWillBeOver() ||
this->timeManager().willBeFinished();
}
/*!
* \brief Returns true if a restart file should be written to
* disk.
*
* The default behaviour is to write one restart file every 5 time
* steps. This file is intented to be overwritten by the
* implementation.
*/
bool shouldWriteRestartFile() const
{
// restart files are only written for certain episodes (250s, 1, 10, 20, 30, 40, 50, 100, 150, 200 years)
return
this->timeManager().timeStepIndex() == 0 ||
(this->timeManager().episodeIndex() == 1 && this->timeManager().episodeIsOver()) ||
(this->timeManager().episodeIndex() == 10 && this->timeManager().episodeIsOver()) ||
(this->timeManager().episodeIndex() == 11 && this->timeManager().episodeIsOver()) ||
(this->timeManager().episodeIndex() == 12 && this->timeManager().episodeIsOver()) ||
(this->timeManager().episodeIndex() == 13 && this->timeManager().episodeIsOver()) ||
(this->timeManager().episodeIndex() == 14 && this->timeManager().episodeIsOver()) ||
(this->timeManager().episodeIndex() == 15 && this->timeManager().episodeIsOver()) ||
(this->timeManager().episodeIndex() == 16 && this->timeManager().episodeIsOver()) ||
this->timeManager().willBeFinished();
}
/*!
* \brief Returns true if a data file should be written to
* disk.
*
* The data files contain cell-wise pressure and saturation information
* for the geomechanical evaluation performed by TNO.
*/
bool shouldWriteDataFile() const
{
// restart files are only written for certain episodes (250s, 1, 2, 5, 10, 30, 40, 50, 100, 150, 200 years)
return
this->timeManager().timeStepIndex() == 0 ||
(this->timeManager().episodeIndex() == 1 && this->timeManager().episodeWillBeOver()) ||
(this->timeManager().episodeIndex() == 2 && this->timeManager().episodeWillBeOver()) ||
(this->timeManager().episodeIndex() == 5 && this->timeManager().episodeWillBeOver()) ||
(this->timeManager().episodeIndex() == 10 && this->timeManager().episodeWillBeOver()) ||
(this->timeManager().episodeIndex() == 11 && this->timeManager().episodeWillBeOver()) ||
(this->timeManager().episodeIndex() == 12 && this->timeManager().episodeWillBeOver()) ||
(this->timeManager().episodeIndex() == 13 && this->timeManager().episodeWillBeOver()) ||
(this->timeManager().episodeIndex() == 14 && this->timeManager().episodeWillBeOver()) ||
(this->timeManager().episodeIndex() == 15 && this->timeManager().episodeWillBeOver()) ||
(this->timeManager().episodeIndex() == 16 && this->timeManager().episodeWillBeOver()) ||
this->timeManager().willBeFinished();
}
/*!
* \brief Called directly after the time integration.
*/
void postTimeStep()
{
// Calculate storage terms
PrimaryVariables storageL, storageG;
Dune::FieldVector<Scalar, 2> flux(0.0);
this->model().globalPhaseStorage(storageL, lPhaseIdx);
this->model().globalPhaseStorage(storageG, gPhaseIdx);
//calculating flux across the layer of the coordinate 2 (Z axis)
//intercepted at coordinate value = 999m (just below the top of the cap rock)
this->model().calculateFluxAcrossLayer(flux, 2, 999);
// Write mass balance information for rank 0
if (this->gridView().comm().rank() == 0) {
std::cout<<"Storage: liquid=[" << storageL << "]"
<< " gas=[" << storageG << "]"<< "flux=[" <<flux << "]\n";
}
}
void episodeEnd()
{
// Start new episode if episode is over
// for first 10 year episode length is 1 year
if(this->timeManager().time() < 3.1536e8 ){
this->timeManager().startNextEpisode(3.1536e7);
std::cout<<"Episode index is set to: "<<this->timeManager().episodeIndex()<<std::endl;
}
// until 50 years episode length is 10 years
else if(this->timeManager().time() < 1.5768e9){
this->timeManager().startNextEpisode(3.1536e8);
std::cout<<"Episode index is set to: "<<this->timeManager().episodeIndex()<<std::endl;
}
// until end of simulation (200 years) episode length is 50 years
else{
this->timeManager().startNextEpisode(1.5768e9);
std::cout<<"Episode index is set to: "<<this->timeManager().episodeIndex()<<std::endl;
}
}
/*!
* \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_; }
/*!
* \brief Returns the temperature within the domain.
*
* This problem assumes a temperature of 10 degrees Celsius.
*/
Scalar temperatureAtPos(const GlobalPosition &pos) const
{
return 273.15 + 4 + (depthBOR_ - depthSeafloor_ - pos[2]) * 0.032;
};
void sourceAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
values = 0;
}
// \}
/*!
* \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
{
const GlobalPosition globalPos = vertex.geometry().center();
//set dirichlet on the top boundary for evaluation of flux across layer
if (onLeftBoundary(globalPos) || onRightBoundary(globalPos) || onBackBoundary(globalPos) || onUpperBoundary(globalPos))
values.setAllDirichlet();
else if(onLowerBoundary(globalPos))
values.setAllNeumann();
else
values.setAllNeumann();
}
/*!
* \brief Evaluate the boundary conditions for a dirichlet
* boundary segment.
*
* \param values The dirichlet values for the primary variables
* \param vertex The vertex for which the boundary type is set
*
* For this method, the \a values parameter stores primary variables.
*/
void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
{
const GlobalPosition globalPos = vertex.geometry().center();
initial_(values, globalPos);
}
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
* \param values The neumann values for the conservation equations
* \param element The finite element
* \param fvElemGeom The finite-volume geometry in the box scheme
* \param is The intersection between element and boundary
* \param scvIdx The local vertex index
* \param boundaryFaceIdx The index of the boundary face
*
* For this method, the \a values parameter stores the mass flux
* in normal direction of each phase. Negative values mean influx.
*/
void neumann(PrimaryVariables &values,
const Element &element,
const FVElementGeometry &fvElemGeom,
const Intersection &is,
int scvIdx,
int boundaryFaceIdx) const
{
const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
values = 0;
// injection well boundary condition for the chimney scenario
if (this->atWell(globalPos) && this->timeManager().time() <= 1.26144e9 )
{
if(chimneyScenario_)
values[contiCO2EqIdx] = -1.692912e-2; // kg/(s*m^2)
else if(catastrophicBlowoutScenario_)
values[contiCO2EqIdx] = -1.692912e-2; // kg/(s*m^2)
else if(faultScenario_)
values[contiCO2EqIdx] = -2.28557e-2; // kg/(s*m^2)
else if(leakyWellScenario_)
values[contiCO2EqIdx] = -4.14319e-2; // kg/(s*m^2)
else
std::cout<<"no scenario specified in input file!"<<std::endl;
}
}
// \}
/*!
* \name Volume terms
*/
// \{
/*!
* \brief Evaluate the initial value for a control volume.
*
* \param values The initial values for the primary variables
* \param element The finite element
* \param fvElemGeom The finite-volume geometry in the box scheme
* \param scvIdx The local vertex index
*
* For this method, the \a values parameter stores primary
* variables.
*/
void initial(PrimaryVariables &values,
const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
{
const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
initial_(values, globalPos);
}
/*!
* \brief Return the initial phase state inside a control volume.
*
* \param vert The vertex
* \param globalIdx The index of the global vertex
* \param globalPos The global position
*/
int initialPhasePresence(const Vertex &vert,
int &globalIdx,
const GlobalPosition &globalPos) const
{ return Indices::wPhaseOnly; }
bool inLeakageFeature(const GlobalPosition &globalPos) const
{
if(chimneyScenario_)
return (globalPos[0] > 3000 && globalPos[0] < 3500 && globalPos[1] < 250 && globalPos[2] > 100);
else if(catastrophicBlowoutScenario_)
return ((globalPos[0] > 2500 && globalPos[0] <= 3000 && globalPos[1] < 1.e-6 && globalPos[2] > 100 && globalPos[2] < 550) ||
(globalPos[0] > 3000 && globalPos[0] < 3500 && globalPos[1] < 250 && globalPos[2] > 100));
else if(faultScenario_)
return (globalPos[0] > 2999 && globalPos[0] < 3011) && (globalPos[2] > 100);
else if(leakyWellScenario_)
return ((globalPos[0] >2999 && globalPos[0] < 3002) && globalPos[1] < 2 && globalPos[2] > 100);
else{
std::cout<<"no scenario specified in input file!"<<std::endl;
return false;
}
}
bool atWell(const GlobalPosition &globalPos) const
{
if(chimneyScenario_)
return (globalPos[2] < 100 && globalPos[1] < eps_ && globalPos[0] > 2485 && globalPos[0] < 2501);
else if(catastrophicBlowoutScenario_)
return (globalPos[2] < 100 && globalPos[1] < eps_ && globalPos[0] > 2485 && globalPos[0] < 2501);
else if(faultScenario_)
return (globalPos[2] < 100 && globalPos[1] < eps_ && globalPos[0] > 2499 && globalPos[0] < 2508);
else if(leakyWellScenario_)
return (globalPos[2] < 100 && globalPos[1] < eps_ && globalPos[0] > 2495 && globalPos[0] < 2503);
else{
std::cout<<"no scenario specified in input file!"<<std::endl;
return false;
}
}
bool onLeftBoundary(const GlobalPosition &globalPos) const
{ return globalPos[0] < this->bboxMin()[0] + eps_; }
bool onRightBoundary(const GlobalPosition &globalPos) const
{ return globalPos[0] > this->bboxMax()[0] - eps_; }
bool onLowerBoundary(const GlobalPosition &globalPos) const
{ return globalPos[2] < this->bboxMin()[2] + eps_; }
bool onUpperBoundary(const GlobalPosition &globalPos) const
{ return globalPos[2] > this->bboxMax()[2] - eps_; }
bool onFrontBoundary(const GlobalPosition &globalPos) const
{return globalPos[1] < this->bboxMin()[1] + eps_; }
bool onBackBoundary(const GlobalPosition &globalPos) const
{return globalPos[1] > this->bboxMax()[1] - eps_; }
// \}
private:
// the internal method for the initial condition
void initial_(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
Scalar temperature = this->temperatureAtPos(globalPos);
Scalar pressure = (1.013e5 + (depthBOR_ - globalPos[2]) * 1100 * 9.81);
Scalar densityW = FluidSystem::Brine::liquidDensity(temperature, pressure);
values[Indices::pressureIdx] = (1.013e5 + (depthBOR_ - globalPos[2]) * densityW * 9.81);
values[Indices::switchIdx] = 0.0;
}
Scalar depthBOR_;
Scalar depthSeafloor_;
Scalar eps_;
bool chimneyScenario_;
bool catastrophicBlowoutScenario_;
bool faultScenario_;
bool leakyWellScenario_;
int nTemperature_;
int nPressure_;
std::string name_ ;
Scalar pressureLow_, pressureHigh_;
Scalar temperatureLow_, temperatureHigh_;
};
} //end namespace
#endif
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
