Hi Mojdeh,
it is possible, but there exists no "easy" interface. Moreover, if the
grid really is nonconforming, our methods currently might not be
suitable. Especially, the box method will fail in such situations. You
can try with a cell centered method, but most likely, it will require
quite a bit of manual fine tuning. If the grid is conforming, it is
rather straightforward.
If you want to try, you have to obtain the modules dune-cornerpoint,
opm-core and opm-porsol from
http://opm-project.org/download.php
If you work within the standard Dumux module, you have to declare the
dependencies on those modules in the file dune.module. And rebuild all
modules with dunecontrol.
Attached are a main file and a problem file that read a grid and
properties from an Eclipse input file and use them in a Dumux model.
They could get you started.
Kind regards
Bernd
On 08/27/2013 10:52 AM, mojdeh rasoulzadeh wrote:
Hi dear DuMuxes,
I am trying to import an Eclipse grid in Dumux.
Is it possible to introduce corner-point grids into dumux model?
Thank you in advance,
Best regards,
Mojdeh
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
--
_______________________________________________________________
!!!! CMWR 2014: 10th - 13th June 2014 in Stuttgart !!!!
Please visit www.cmwr14.de
_______________________________________________________________
Bernd Flemisch phone: +49 711 685 69162
IWS, Universität Stuttgart fax: +49 711 685 60430
Pfaffenwaldring 61 email: [email protected]
D-70569 Stuttgart url: www.hydrosys.uni-stuttgart.de
_______________________________________________________________
/*****************************************************************************
* Copyright (C) 2011 by Francesco Patacchini *
* Copyright (C) 2011 by Andreas Lauser *
* Institute of Hydraulic Engineering *
* 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/>. *
*****************************************************************************/
#include "config.h"
#include "spe9problem.hh"
#include <dumux/common/start.hh>
void usage(const char *progname)
{
std::cout << "usage: " << progname << " fileformat=eclipse filename=FILENAME\n"
<< "\n"
<< "Where FILENAME is the name of the file which describes the\n"
<< "SPE-9 problem in ECLiPSE format\n";
}
int main(int argc, const char** argv)
{
#ifdef NDEBUG
try {
#endif
if (argc == 1) {
usage(argv[0]);
return 0;
}
typedef TTAG(Spe9Problem) TypeTag;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
typedef Opm::ReservoirPropertyCapillary<3> ReservoirProperties;
Dumux::Properties::print<TypeTag>();
double tEnd = 90000.0*24*60*60;
//double tEnd = 90.0*24*60*60;
double dt = 1;
// create grid
// -> load the grid from file
Opm::parameter::ParameterGroup param(argc, argv);
Grid grid;
ReservoirProperties resProps;
Opm::setupGridAndProps(param, grid, resProps);
Opm::EclipseGridParser parser(param.get<std::string>("filename"));
// instantiate and run the concrete problem
TimeManager timeManager;
Problem problem(timeManager, grid.leafView(), parser);
problem.setReservoirProperties(&resProps);
timeManager.init(problem, 0, dt, tEnd, /*restart=*/false);
timeManager.run();
return 0;
#ifdef NDEBUG
}
catch (Dune::Exception &e) {
std::cerr << "Dune reported error: " << e << std::endl;
}
catch (...) {
std::cerr << "Unknown exception thrown!\n";
}
#endif
return 3;
}
/*****************************************************************************
* Copyright (C) 2011 by Francesco Patacchini *
* Copyright (C) 2011 by Andreas Lauser *
* Institute of Hydraulic Engineering *
* 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
* \ingroup TwoPTwoNBoxProblems
*
* \brief Problem where water is injected by means of a
* dirchlet condition on the lower right of the domain which have to go
* around an obstacle with \f$10^3\f$ lower permeability.
* \author Andreas Lauser
*/
#ifndef DUMUX_SPE9_PROBLEM_HH
#define DUMUX_SPE9_PROBLEM_HH
#include <dumux/implicit/mpnc/mpncmodel.hh>
#include <dumux/material/fluidstates/immisciblefluidstate.hh>
#include <dumux/material/fluidstates/compositionalfluidstate.hh>
#include <dumux/implicit/common/implicitporousmediaproblem.hh>
#include <dune/grid/uggrid.hh>
#include <opm/porsol/common/setupGridAndProps.hpp>
//#include <opm/core/common/param/ParameterGroup.hpp>
#include "spe9law.hh"
#include "spe9well.hh"
#include "spe9fluidsystem.hh"
#include "spe9compositionfromfugacities.hh"
namespace Dumux
{
template <class TypeTag>
class Spe9Problem;
namespace Properties
{
NEW_TYPE_TAG(Spe9Problem, INHERITS_FROM(BoxMPNC));
// Set the grid type
SET_TYPE_PROP(Spe9Problem, Grid, Dune::CpGrid);
//SET_PROP(Spe9Problem, Grid) { typedef Dune::ALUCubeGrid<3,3> type; };
// Set the problem property
SET_TYPE_PROP(Spe9Problem,
Problem,
Dumux::Spe9Problem<TypeTag>);
// Set fluid configuration
SET_PROP(Spe9Problem, FluidSystem)
{
private: typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
public: typedef Dumux::Spe9FluidSystem<Scalar> type;
};
// Use the SPE-9 problem specific composition from fugacities
// thermodynamic constraint solver.
SET_PROP(Spe9Problem, ConstraintSolver)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
public:
typedef Dumux::Spe9CompositionFromFugacities<Scalar, FluidSystem> type;
};
// Set the material Law
SET_PROP(Spe9Problem, MaterialLaw)
{
private:
// define the material law which is parameterized by effective
// saturations
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
public:
typedef Dumux::Spe9Law<Scalar, FluidSystem> type;
};
// Enable molecular diffusion of the components?
SET_BOOL_PROP(Spe9Problem, EnableDiffusion, false);
// Enable gravity
SET_BOOL_PROP(Spe9Problem, ProblemEnableGravity, true);
// Write Newton convergence to disk?
SET_BOOL_PROP(Spe9Problem, NewtonWriteConvergence, false);
// Enable the re-use of the jacobian matrix whenever possible?
SET_BOOL_PROP(Spe9Problem, ImplicitEnableJacobianRecycling, true);
// use forward differences to calculate the jacobian
SET_INT_PROP(Spe9Problem, ImplicitNumericDifferenceMethod, +1);
// Reassemble the jacobian matrix only where it changed?
SET_BOOL_PROP(Spe9Problem, ImplicitEnablePartialReassemble, true);
// This problem uses cheap equilibrium calculations, so we don't need
// hints
SET_BOOL_PROP(Spe9Problem, ImplicitEnableHints, false);
// Enable smooth upwinding?
SET_BOOL_PROP(Spe9Problem, ImplicitEnableSmoothUpwinding, true);
// Tolerance of the non-linear solver
SET_SCALAR_PROP(Spe9Problem, NewtonRelTolerance, 1e-6);
// Residual reduction of the linear solver
SET_SCALAR_PROP(Spe9Problem, LinearSolverResidualReduction, 1e-6);
}
/*!
* \ingroup MpNcBoxProblems
* \brief Fifth benchmark problem from the society of petrolium engineers.
*
* See: J.E. Killough, et al.: Fifth Comparative Solution Project:
* Evaluation of Miscible Flood Simulators, Ninth SPE Symposium on
* Reservoir Simulation, 1987
*/
template <class TypeTag>
class Spe9Problem
: public ImplicitPorousMediaProblem<TypeTag>
{
typedef Spe9Problem<TypeTag> ThisType;
typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLawParams)) MaterialLawParams;
enum {
// Grid and world dimension
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
};
typedef typename GridView::template Codim<0>::Entity Element;
typedef Dune::FieldVector<typename GridView::Grid::ctype, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
// copy some indices for convenience
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Indices)) Indices;
enum {
// number of phases and components
numPhases = FluidSystem::numPhases,
numComponents = FluidSystem::numComponents,
// component indices
gCompIdx = FluidSystem::gCompIdx,
wCompIdx = FluidSystem::wCompIdx,
oCompIdx = FluidSystem::oCompIdx,
// phase indices
gPhaseIdx = FluidSystem::gPhaseIdx,
oPhaseIdx = FluidSystem::oPhaseIdx,
wPhaseIdx = FluidSystem::wPhaseIdx,
// primary variable indices
fug0Idx = Indices::fug0Idx,
S0Idx = Indices::S0Idx,
p0Idx = Indices::p0Idx,
// equation indices
conti0EqIdx = Indices::conti0EqIdx,
phase0NcpIdx = Indices::phase0NcpIdx,
// number of initial phases
numPhaseStates = 3,
// epsiode types
doNothingEpisode = 0,
produceFullInjectEpisode = 1,
produceLittleInjectEpisode = 2,
};
typedef Dumux::CompositionalFluidState<Scalar, FluidSystem> CompositionalFluidState;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef Dumux::Spe9Well<TypeTag> WellModel;
typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
typedef Opm::EclipseGridParser Parser;
typedef Opm::ReservoirPropertyCapillary<3> ReservoirProperties;
typedef std::vector<std::pair<Scalar, Scalar> > SSP;
public:
Spe9Problem(TimeManager &timeManager, const GridView &gridView, const Parser& parser)
: ParentType(timeManager, gridView), parser_(parser)
{
// enable gravity
if (GET_PROP_VALUE(TypeTag, PTAG(ProblemEnableGravity)))
this->gravity_[dim-1] = 9.81;
readFromEclipseParser_(parser);
////////////////////
// initialize fluid system
////////////////////
FluidSystem::initBegin();
int regionNumber = 0;
enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
// Surface densities in [kg/m^3]. Original values in [lb/ft^3]
if (parser_.hasField("DENSITY")) {
const std::vector<double>& d = parser_.getDENSITY().densities_[regionNumber];
FluidSystem::setSurfaceDensities(d[ECL_oil], d[ECL_water], d[ECL_gas]);
std::cout << "densities: oil, water, gas\n" << d[ECL_oil]
<< ", " << d[ECL_water] << ", " << d[ECL_gas] << std::endl;
}
else {
THROW("Input is missing DENSITY\n");
}
////////////////////
// oil phase parameters
////////////////////
// set the gas formation factor in [m^3/m^3] depending on the
// pressure in [Pa]. This is equivalent to the volumetric
// amount of gas at standard conditions which can emerge from
// a given quantity of liquid oil at reservoir pressure (and
// probably temperature) if this oil is brought to the surface
// (i.e. it exhibits surface density).
std::cout << "\nreading oil PVT:\nRs\t\tp\t\tB\t\tmu\n";
typedef std::vector<std::vector<std::vector<double> > > TableT;
if (parser_.hasField("PVTO")) {
const TableT& pvto = parser_.getPVTO().pvto_;
const int sz = pvto[regionNumber].size();
SSP boTable, muTable, rsTable;
for (int i = 0; i < sz; i++)
{
typename SSP::value_type pB(pvto[regionNumber][i][1], pvto[regionNumber][i][2]);
boTable.push_back(pB);
typename SSP::value_type pMu(pvto[regionNumber][i][1], pvto[regionNumber][i][3]);
muTable.push_back(pMu);
typename SSP::value_type pRs(pvto[regionNumber][i][1], pvto[regionNumber][i][0]);
rsTable.push_back(pRs);
std::cout << std::scientific << rsTable[i].second << "\t"
<< boTable[i].first << "\t" << boTable[i].second
<< "\t" << muTable[i].second << std::endl;
}
FluidSystem::setOilFormationVolumeFactor(boTable);
FluidSystem::setOilViscosity(muTable);
FluidSystem::setGasFormationFactor(rsTable);
}
else if (parser_.hasField("PVDO")) {
const TableT& pvdo = parser_.getPVDO().pvdo_;
const int sz = pvdo[regionNumber][0].size();
SSP boTable, muTable;
for (int i = 0; i < sz; i++)
{
typename SSP::value_type pB(pvdo[regionNumber][0][i], pvdo[regionNumber][1][i]);
boTable.push_back(pB);
typename SSP::value_type pMu(pvdo[regionNumber][0][i], pvdo[regionNumber][2][i]);
muTable.push_back(pMu);
std::cout << std::scientific << boTable[i].first << "\t" << boTable[i].second
<< "\t" << muTable[i].second << std::endl;
}
FluidSystem::setOilFormationVolumeFactor(boTable);
FluidSystem::setOilViscosity(muTable);
}
else {
THROW("Input is missing PVTO or PVDO\n");
}
FluidSystem::setReferenceVolumeFactor(oPhaseIdx, 1.0);
////////////////////
// gas phase parameters
////////////////////
// set the gas formation volume factor (surface
// volume/reservoir volume). ATTENTION: This factor is NOT
// normalized to surface conditions!
// set the gas viscosity [Pa s] dependent on pressure. The
// original units are [psia] for pressure and [cP] for
// viscosity.
std::cout << "\nreading gas PVT:\np\t\tB\t\tmu\n";
typedef std::vector<std::vector<std::vector<double> > > TableT;
if (parser_.hasField("PVTG")) {
const TableT& pvtg = parser_.getPVTG().pvtg_;
const int sz = pvtg[regionNumber].size();
SSP bgTable, muTable;
for (int i = 0; i < sz; i++)
{
typename SSP::value_type pB(pvtg[regionNumber][i][0], pvtg[regionNumber][i][2]);
bgTable.push_back(pB);
typename SSP::value_type pMu(pvtg[regionNumber][i][0], pvtg[regionNumber][i][3]);
muTable.push_back(pMu);
std::cout << std::scientific << bgTable[i].first << "\t" << bgTable[i].second
<< "\t" << muTable[i].second << std::endl;
}
FluidSystem::setGasFormationVolumeFactor(bgTable);
FluidSystem::setGasViscosity(muTable);
}
else if (parser_.hasField("PVDG")) {
const TableT& pvdg = parser_.getPVDG().pvdg_;
const int sz = pvdg[regionNumber][0].size();
SSP bgTable, muTable;
for (int i = 0; i < sz; i++)
{
typename SSP::value_type pB(pvdg[regionNumber][0][i], pvdg[regionNumber][1][i]);
bgTable.push_back(pB);
typename SSP::value_type pMu(pvdg[regionNumber][0][i], pvdg[regionNumber][2][i]);
muTable.push_back(pMu);
std::cout << std::scientific << bgTable[i].first << "\t" << bgTable[i].second
<< "\t" << muTable[i].second << std::endl;
}
FluidSystem::setGasFormationVolumeFactor(bgTable);
FluidSystem::setGasViscosity(muTable);
}
else {
THROW("Input is missing PVTG or PVDG\n");
}
FluidSystem::setReferenceVolumeFactor(gPhaseIdx, 1.0);
////////////////////
// water phase parameters
////////////////////
std::cout << "\nreading water PVT:\ncompressibility\tmu\n";
typedef std::vector<std::vector<double> > TableW;
if (parser_.hasField("PVTW")) {
const TableW& pvtw = parser_.getPVTW().pvtw_;
if (pvtw.size() != 1) {
THROW("More than one PVD-region");
}
// set the water compressibility factor in [1/psi]
FluidSystem::setWaterCompressibility(pvtw[regionNumber][2]/*1e-6/6894.7573*/);
// set the dynamic water viscosity [Pa s]. Original value in [cP]
FluidSystem::setWaterViscosity(pvtw[regionNumber][3]/*0.96 * 1e-3*/);
std::cout << std::scientific << pvtw[regionNumber][2] << "\t" << pvtw[regionNumber][3] << std::endl;
}
else {
THROW("Input is missing PVTW\n");
}
std::cout.unsetf(std::ios_base::floatfield);
FluidSystem::initEnd();
std::cout << "Oil bubble pressure: " << FluidSystem::bubblePressure() << "Pa\n";
std::cout << "Oil reservoir pressure: " << pRefOil_() << "Pa\n";
////////////////////
// parameters which are not directly related to the fluid
// system
////////////////////
// temperature. this is actually unspecified by the SPE-9. We
// just assume 311 K (37.85 degrees Celsius)
temperature_ = 311;
// set the first episode to be 3 years long. in the first
// episode we don't do anything and just wait to get
// hydrostatic pressure
this->timeManager().startNextEpisode(3*365*24*60*60);
episodeType_ = doNothingEpisode;
////////////////////
// initialize fluid state and wells
////////////////////
for (int i = 0; i < numPhaseStates; ++i)
setupReservoirFluidState_(i);
// initialize the wells
initWells_();
// setup the fluid states for the water injection well
setupInjectionFluids_();
// start a water injection episode
setInjectionWell_(/*pBottomPSI=*/5e3, /*pMaxRate=*/0);
}
void episodeEnd()
{
Scalar duration, dt;
if (episodeType_ == doNothingEpisode)
{
// start a water injection episode
episodeType_ = produceFullInjectEpisode;
duration = 1e100; // 900*24*60*60;
dt = 1;
#warning hack: should be 5e3 STB/d
setInjectionWell_(/*pBottomPSI=*/5e3, /*maxRateSTB/d=*/5e9);
}
else
{
duration = 1e100;
dt = 1;
}
std::cout << "End of episode " << this->timeManager().episodeIndex()
<< ". Start next episode of " << duration/(24*3600) << " days "
<< "with time step size " << dt << std::endl;
this->model().jacobianAssembler().setMatrixReuseable(false);
this->timeManager().startNextEpisode(duration);
this->timeManager().setTimeStepSize(dt);
}
/*!
* \brief Called directly after the time integration.
*/
void postTimeStep()
{
ComponentVector injectionRate;
ComponentVector productionRate;
PhaseVector avgPressure;
gatherStatistics_(injectionRate, productionRate, avgPressure);
printRate_("injection", injectionRate);
printRate_("production", productionRate);
if (this->gridView().comm().rank() == 0)
std::cout << "Average phase pressure [Pa]: " << avgPressure << "\n";
updateProductionWell_(productionRate);
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief The problem name.
*
* This is used as a prefix for files generated by the simulation.
*/
const char *name() const
{ return "spe9"; }
/*!
* \brief Returns the temperature within the domain.
*
* This problem assumes a temperature of 30 degrees Celsius. This
* method is only called if the problem is isothermal.
*/
Scalar temperature() const
{ return temperature_; };
// \}
/*!
* \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 on the boundary for which the
* conditions needs to be specified
*/
void boundaryTypesAtPos(BoundaryTypes &values, const GlobalPosition &pos) const
{
// only no-flow boundaries. wells are modeled by the source
// term!
values.setAllNeumann();
}
void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &pos) const
{
values = 0;
};
void neumannAtPos(PrimaryVariables &values, const GlobalPosition &pos) const
{
values = 0;
};
// \}
/*!
* \name Volume terms
*/
// \{
/*!
* \brief Evaluate the source term for all phases within a given
* sub-control-volume.
*
* For this method, the \a values parameter stores the rate mass
* of a component is generated or annihilate per volume
* unit. Positive values mean that mass is created, negative ones
* mean that it vanishes.
*/
void solDependentSource(PrimaryVariables &values,
const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx,
const ElementVariables &elemVolVars) const
{
values = 0;
PrimaryVariables q;
#if 1
Valgrind::SetUndefined(q);
injectionWell_.source(q, element, fvGeometry, scvIdx, elemVolVars);
Valgrind::CheckDefined(q);
values += q;
#endif
#if 1
Valgrind::SetUndefined(q);
productionWell_.source(q, element, fvGeometry, scvIdx, elemVolVars);
Valgrind::CheckDefined(q);
values += q;
#endif
}
// \}
/*!
* \brief Write a restart file?
*/
bool shouldWriteRestartFile() const
{
return false;
//return ParentType::shouldWriteRestartFile();
}
bool shouldWriteOutput() const
{
//return false;
return true;
}
void serialize()
{ ParentType::serialize(); }
/*!
* \brief This method writes the complete state of the problem
* to the harddisk.
*
* The file will start with the prefix returned by the name()
* method, has the current time of the simulation clock in it's
* name and uses the extension <tt>.drs</tt>. (Dumux ReStart
* file.) See Dumux::Restart for details.
*
* \tparam Restarter The serializer type
*
* \param res The serializer object
*/
template <class Restarter>
void serialize(Restarter &res)
{
ParentType::serialize(res);
injectionWell_.serialize(res);
productionWell_.serialize(res);
}
/*!
* \brief This method restores the complete state of the problem
* from disk.
*
* It is the inverse of the serialize() method.
*
* \tparam Restarter The deserializer type
*
* \param res The desrializer object
*/
template <class Restarter>
void deserialize(Restarter &res)
{
ParentType::deserialize(res);
injectionWell_.deserialize(res);
productionWell_.deserialize(res);
};
/*!
* \brief Returns the intrinsic permeability tensor.
*
* \param element The current finite element
* \param fvElemGeom The current finite volume geometry of the element
* \param scvIdx The index sub-control volume where the
* intrinsic permeability is given.
*/
Scalar intrinsicPermeability(const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
{
int elemIdx = this->elementMapper().map(element);
const typename ReservoirProperties::PermTensor &permRes = resProps_->permeability(elemIdx);
Tensor K;
for (int i = 0; i < dim; i++)
for (int j = 0; j < dim; j++)
K[i][j] = permRes(i, j);
return K[0][0];
}
/*!
* \brief Define the porosity \f$[-]\f$ of the soil
*
* \param element The finite element
* \param fvElemGeom The finite volume geometry
* \param scvIdx The local index of the sub-control volume where
* the porosity needs to be defined
*/
template <class Context>
Scalar porosity(const Context &context, int spaceIdx, int timeIdx) const
{
int elemIdx = this->elementMapper().map(context.element());
return resProps_->porosity(elemIdx);
}
// return the brooks-corey context depending on the position
template <class Context>
const MaterialLawParams& materialLawParams(const Context &context, int spaceIdx, int timeIdx) const
{ return materialParams_; }
void setReservoirProperties(ReservoirProperties* resProps)
{ resProps_ = resProps; }
const ReservoirProperties& reservoirProperties() const
{ return *resProps_; }
private:
// do some statistics after a time step
void gatherStatistics_(ComponentVector &injectionRate,
ComponentVector &productionRate,
PhaseVector &avgPressure)
{
injectionRate = 0;
productionRate = 0;
avgPressure = 0;
PrimaryVariables q;
Scalar totalVolume = 0;
ElementVariables elemVars;
ElementIterator elemIt = this->gridView().template begin<0>();
const ElementIterator elemEndIt = this->gridView().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) {
FVElementGeometry fvGeometry;
fvGeometry.update(this->gridView(), *elemIt);
elemVars.update(*this, *elemIt, fvGeometry, false);
int numScv = fvGeometry.numScv;
for (int i = 0; i < numScv; ++i)
{
Scalar Vscv = fvGeometry.subContVol[i].volume;
injectionWell_.source(q, *elemIt, fvGeometry, i, elemVars);
q *= Vscv;
for (int j = 0; j < numComponents; ++j)
injectionRate[j] += q[conti0EqIdx + j];
if (productionWell_.applies(*elemIt, i)) {
if (fvGeometry.subContVol[i].global[2] >= productionWell_.wellDepth()) {
std::cout << "Pressure in the deepest production cell = "
<< elemVars[i].fluidState().pressure(oPhaseIdx)
<< " Pa\n";
}
}
productionWell_.source(q, *elemIt, fvGeometry, i, elemVars);
q *= Vscv;
for (int j = 0; j < numComponents; ++j)
productionRate[j] += q[conti0EqIdx + j];
totalVolume += Vscv;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
avgPressure[phaseIdx] +=
elemVars[i].fluidState().pressure(phaseIdx)
* Vscv;
}
};
}
injectionRate = this->gridView().comm().sum(injectionRate);
productionRate = this->gridView().comm().sum(productionRate);
avgPressure = this->gridView().comm().sum(avgPressure);
totalVolume = this->gridView().comm().sum(totalVolume);
avgPressure /= totalVolume;
};
void printRate_(const std::string &sourceType,
const ComponentVector &compRates)
{
if (this->gridView().comm().rank() != 0)
return;
std::cout << sourceType << " rates [mol/s]: " << compRates << "\n";
std::cout << sourceType << " rates [kg/s]: ";
for (int i = 0; i < numComponents; ++i)
std::cout << compRates[i]*FluidSystem::molarMass(i) << " ";
std::cout << "\n";
Scalar Vo = 0.0;
Scalar Vw = 0.0;
Scalar Vg = 0.0;
typename FluidSystem::ParameterCache paramCache;
// create a fluid state which assumes immiscibility. this
// implicitly sets the composition
ImmiscibleFluidState<Scalar, FluidSystem> fs;
fs.setTemperature((60 + 459.67)*5/9); // SPE standard temperature
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
fs.setPressure(phaseIdx, 14.7 * 6894.7573); // SPE standard pressure
paramCache.updateAll(fs);
Vg = compRates[gCompIdx] * FluidSystem::molarMass(gCompIdx) / FluidSystem::density(fs, paramCache, gPhaseIdx);
Vo = compRates[oCompIdx] * FluidSystem::molarMass(oCompIdx) / FluidSystem::density(fs, paramCache, oPhaseIdx);
Vw = compRates[wCompIdx] * FluidSystem::molarMass(wCompIdx) / FluidSystem::density(fs, paramCache, wPhaseIdx);
ComponentVector massRates(compRates);
massRates *= fs.averageMolarMass(oPhaseIdx);
std::cout << "Volumetric "<<sourceType<<" rates:\n";
if (sourceType == "injection")
std::cout << " water = " << Vw << "m^3/s (=" << Vw*(24*60*60)/0.1591132 << "stb/d)\n";
else {
std::cout << " oil = " << -Vo << "m^3/s (=" << -Vo*(24*60*60)/0.1591132 << "stb/d)\n";
std::cout << " water = " << -Vw << "m^3/s (=" << -Vw*(24*60*60)/0.1591132 << "stb/d)\n";
std::cout << " gas = " << -Vg << "m^3/s (=" << -Vg*(24*60*60) << "m^3/d)\n";
}
}
// initialize the injection and extraction wells
void initWells_()
{
GlobalPosition brange = this->bBoxMax();
brange -= this->bBoxMin();
injectionWell_.beginSpec(*this);
productionWell_.beginSpec(*this);
ComponentVector prodComp(0.0);
updateProductionWell_(prodComp);
///
// add sub control volumes to the wells
ElementIterator elemIt = this->gridView().template begin<0>();
ElementIterator elemEndIt = this->gridView().template end<0>();
ElementVariables elemVars;
Scalar tagInj = 0;
Scalar tagProd = 0;
for (; elemIt != elemEndIt; ++elemIt) {
FVElementGeometry fvGeometry;
fvGeometry.update(this->gridView(), *elemIt);
elemVars.update(*this, *elemIt, fvGeometry, false);
int n = elemIt->template count<dim>();
for (int i = 0; i < n; ++i) {
GlobalPosition pos = elemIt->geometry().corner(i);
if (pos[0] > 2103.0174*cos(10*M_PI/180) - 50 &&
pos[0] < 2103.0174*cos(10*M_PI/180) + 50 &&
pos[1] > brange[1] - 91.435538 - 50 &&
pos[1] < brange[1] - 91.435538 + 50 &&
pos[2] > 2743.0661 + 2103.0174*sin(10*M_PI/180) + 42.669918 - 1 &&
pos[2] < 2743.0661 + 2103.0174*sin(10*M_PI/180) + 42.669918 + 36.26943 + 1)
{
tagInj += 1;
if (tagInj == 1 || tagInj == 3 || tagInj == 5 || tagInj == 7 || tagInj == 9)
{
std::cout << "injection well at position " << pos << "\n";
injectionWell_.addScv(elemVars, *elemIt, fvGeometry, i);
}
}
if (pos[0] > 457.17769*cos(10*M_PI/180) - 50 &&
pos[0] < 457.17769*cos(10*M_PI/180) + 50 &&
pos[1] > brange[1] - 91.435538 - 50 &&
pos[1] < brange[1] - 91.435538 + 50 &&
pos[2] > 2743.0661 + 457.17769*sin(10*M_PI/180) + 10.667479 - 1 &&
pos[2] < 2743.0661 + 457.17769*sin(10*M_PI/180) + 10.667479 + 12.49619 + 1)
{
tagProd += 1;
if (tagProd == 1 || tagProd == 3 || tagProd == 5)
{
std::cout << "production well at position " << pos << "\n";
productionWell_.addScv(elemVars, *elemIt, fvGeometry, i);
}
}
};
injectionWell_.endSpec();
}
}
void setInjectionWell_(Scalar pBottomPsi, Scalar VmaxSTB)
{
// First, calculate the total amount of injected water
// in mol/s
Scalar Vinj =
VmaxSTB * 0.1591132/(60*60*24); // injected volume in [m^3/s]
//////////////////
// -> injection borehole
//////////////////
injectionWell_.setInjection(injectWFluidState_, Vinj);
// set the maximum bottom hole pressure as specified in SPE9
Scalar pBottom = pBottomPsi * 6894.7573; // bottom hole pressure [Pa]
injectionWell_.setBottomPressure(pBottom);
}
void updateProductionWell_(const ComponentVector &producedRates)
{
#if 1
Scalar maxVolumeRate;
// produced volume in [stdb/day] as specified by SPE9, changing with time
// scenario 1
if (episodeType_ == doNothingEpisode)
maxVolumeRate = 0.0;
else if (episodeType_ == produceFullInjectEpisode)
#warning hack: should be 1.5e3 STB/d
// first phase
maxVolumeRate = 1.5e9;
else {
assert(episodeType_ == produceLittleInjectEpisode);
// second phase
maxVolumeRate = 1.0e2;
}
// convert to [m^3/s]
maxVolumeRate *= 0.1591132/(60*60*24);
productionWell_.setProduction(maxVolumeRate);
#else
productionWell_.setProduction(0);
#endif
// set the minimum bottom hole pressure as specified in SPE9
Scalar pBottom = 1e3; // [PSI]
pBottom *= 6894.7573; // convert to [Pa]
productionWell_.setBottomPressure(pBottom);
};
void setupInjectionFluids_()
{
//////////
// set the saturations and temperatures for the water injection
//////////
auto &fs = injectWFluidState_;
// set the temperature
fs.setTemperature((60 + 459.67)*5/9); // SPE standard temperature
// 100% phase water, standard conditions
fs.setSaturation(gPhaseIdx, 0.0);
fs.setSaturation(oPhaseIdx, 0.0);
fs.setSaturation(wPhaseIdx, 1.0);
// set the remaining parameters
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
fs.setPressure(phaseIdx, 14.7 * 6894.7573); // SPE standard pressure
// 100% water
fs.setMoleFraction(phaseIdx, wCompIdx, 1.0);
fs.setMoleFraction(phaseIdx, gCompIdx, 0.0);
fs.setMoleFraction(phaseIdx, oCompIdx, 0.0);
}
typename FluidSystem::ParameterCache paramCache;
paramCache.updateAll(fs);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
fs.setDensity(phaseIdx, rho);
Scalar mu = FluidSystem::viscosity(fs, paramCache, phaseIdx);
fs.setViscosity(phaseIdx, mu);
}
};
void setupReservoirFluidState_(int regionIdx)
{
//////////
// set temperatures
//////////
Dumux::CompositionalFluidState<Scalar, FluidSystem> fs;
// constant temperature
fs.setTemperature(temperature_);
//////////
// calculate the pressures
//////////
// calculate capillary pressures
const MaterialLawParams &matParams = materialParams_;
Scalar pC[numPhases];
MaterialLaw::capillaryPressures(pC, matParams, fs);
// Calculate the pressures of the phases by adding the
// difference of their capillary pressure to the capillary
// pressure of the reference pressure
Scalar p = pRefOil_();
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
fs.setPressure(phaseIdx,
p +
(pC[phaseIdx] - pC[oPhaseIdx]));
//////////
// calculate the saturations and the composition of all phases
//////////
int referencePhaseIdx;
if (regionIdx == oPhaseIdx) {
referencePhaseIdx = oPhaseIdx;
fs.setSaturation(gPhaseIdx, 0.0);
fs.setSaturation(wPhaseIdx, 0.15109);
fs.setSaturation(oPhaseIdx, 1.0 - 0.15109);
Scalar pBubble = FluidSystem::bubblePressure();
Scalar x_oG =
0.9*
FluidSystem::fugCoefficientInGas_(gCompIdx, pBubble) /
FluidSystem::fugCoefficientInOil_(gCompIdx, pBubble);
Scalar x_oO = 1 - x_oG;
Scalar x_oW =
FluidSystem::fugCoefficientInWater_(wCompIdx, p) /
FluidSystem::fugCoefficientInOil_(wCompIdx, p);
fs.setMoleFraction(oPhaseIdx, oCompIdx, (1 - x_oW)*x_oO);
fs.setMoleFraction(oPhaseIdx, gCompIdx, (1 - x_oW)*x_oG);
fs.setMoleFraction(oPhaseIdx, wCompIdx, x_oW);
}
else if (regionIdx == wPhaseIdx) {
referencePhaseIdx = oPhaseIdx;
fs.setSaturation(gPhaseIdx, 0.0);
fs.setSaturation(wPhaseIdx, 0.881490);
fs.setSaturation(oPhaseIdx, 1.0 - fs.saturation(wPhaseIdx));
Scalar pBubble = FluidSystem::bubblePressure();
Scalar x_oG =
0.9*
FluidSystem::fugCoefficientInGas_(gCompIdx, pBubble) /
FluidSystem::fugCoefficientInOil_(gCompIdx, pBubble);
Scalar x_oO = 1 - x_oG;
Scalar x_oW =
FluidSystem::fugCoefficientInWater_(wCompIdx, p) /
FluidSystem::fugCoefficientInOil_(wCompIdx, p);
fs.setMoleFraction(oPhaseIdx, oCompIdx, (1 - x_oW)*x_oO);
fs.setMoleFraction(oPhaseIdx, gCompIdx, (1 - x_oW)*x_oG);
fs.setMoleFraction(oPhaseIdx, wCompIdx, x_oW);
}
else {
/*DUNE_THROW(Dune::InvalidStateException,
"Only gas is not a valid initial value");*/
return;
}
// compute fugacities from the oil phase composition
if (this->gridView().comm().rank() == 0)
std::cout << "Fugacities " << " ------------------------\n";
Dune::FieldVector<Scalar, numComponents> fugVec;
typename FluidSystem::ParameterCache paramCache;
paramCache.updatePhase(fs, referencePhaseIdx);
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
Scalar phi = FluidSystem::fugacityCoefficient(fs, paramCache, referencePhaseIdx, compIdx);
fs.setFugacityCoefficient(referencePhaseIdx, compIdx, phi);
fugVec[compIdx] = fs.fugacity(referencePhaseIdx, compIdx);
if (this->gridView().comm().rank() == 0)
std::cout << "f_" << compIdx << " = " << fugVec[compIdx] << "\n";
};
// compute final phase compositions assuming thermodynamic
// equilibrium
typedef Dumux::CompositionFromFugacities<Scalar, FluidSystem> CompFromFug;
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
CompFromFug::guessInitial(fs, paramCache, phaseIdx, fugVec);
CompFromFug::solve(fs, paramCache, phaseIdx, fugVec);
if (this->gridView().comm().rank() == 0) {
std::cout << FluidSystem::phaseName(phaseIdx) << "Phase ------------------------\n";
std::cout << "rho: " << fs.density(phaseIdx) << "\n";;
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
std::cout << "x_" << FluidSystem::phaseName(phaseIdx)
<< "^" << FluidSystem::componentName(compIdx)
<< " = " << fs.moleFraction(phaseIdx, compIdx) << "\n";
}
}
// Assign the result to the reservoir fluid state
reservoirFluidState_[regionIdx].assign(fs);
}
Scalar pRefOil_() const
{
#warning "HACK: reservoir pressure of oil must be larger than the bubble pressure since the oil-gas contact is above the reservoir!"
//return 4000*6894.76; // set to 4000 psi because we don't know better...
return 3600*6894.76;
}
template<class Parser>
void readFromEclipseParser_(const Parser& parser)
{
// Extract input data.
const Opm::SGOF::table_t& sgof_table = parser.getSGOF().sgof_;
const Opm::SWOF::table_t& swof_table = parser.getSWOF().swof_;
if (sgof_table.size() != 1 || swof_table.size() != 1) {
std::cerr << "We must have exactly one SWOF and one SGOF table (at the moment).\n";
throw std::logic_error("Not implemented");
}
const std::vector<double>& sw = swof_table[0][0];
const std::vector<double>& krw = swof_table[0][1];
const std::vector<double>& krow = swof_table[0][2];
const std::vector<double>& pcow = swof_table[0][3];
std::cout << "\nreading SWOF:\nSw\t\tkrw\t\tkrow\t\tpcow\n";
const int sz = sw.size();
SSP krwTable, krowTable, pcowTable;
for (int i = 0; i < sz; i++)
{
typename SSP::value_type swKrw(sw[i], krw[i]);
krwTable.push_back(swKrw);
typename SSP::value_type swKrow(sw[i], krow[i]);
krowTable.push_back(swKrow);
typename SSP::value_type swPcow(sw[i], pcow[i]);
pcowTable.push_back(swPcow);
std::cout << std::scientific << sw[i] << "\t" << krw[i] << "\t"
<< krow[i] << "\t" << pcow[i] << std::endl;
}
materialParams_.setKrw(krwTable);
materialParams_.setKrow(krowTable);
materialParams_.setPcow(pcowTable);
const std::vector<double>& sg = sgof_table[0][0];
const std::vector<double>& krg = sgof_table[0][1];
const std::vector<double>& krog = sgof_table[0][2];
const std::vector<double>& pcog = sgof_table[0][3];
std::cout << "\nreading SGOF:\nSg\t\tkrg\t\tkrog\t\tpcog\n";
const int szg = sg.size();
SSP krgTable, krogTable, pcogTable;
for (int i = 0; i < szg; i++)
{
typename SSP::value_type sgKrg(sg[i], krg[i]);
krgTable.push_back(sgKrg);
typename SSP::value_type sgKrog(sg[i], krog[i]);
krogTable.push_back(sgKrog);
typename SSP::value_type sgPcog(sg[i], pcog[i]);
pcogTable.push_back(sgPcog);
std::cout << std::scientific << sg[i] << "\t" << krg[i] << "\t"
<< krog[i] << "\t" << pcog[i] << std::endl;
}
materialParams_.setKrg(krgTable);
materialParams_.setKrl(krogTable);
materialParams_.setPcgo(pcogTable);
std::cout.unsetf(std::ios_base::floatfield);
#if 0
// plot the capillary pressure curves
GenericFluidState<Scalar, FluidSystem> fs;
int n = 100;
for (int i = -n/10; i < n + n/10; ++i) {
Scalar Sw = Scalar(i)/n;
fs.setSaturation(wPhaseIdx, Sw);
for (int j = -n/10; j < n + n/10; ++j) {
Scalar Sg = Scalar(j)/n;
fs.setSaturation(wPhaseIdx, Sw);
fs.setSaturation(oPhaseIdx, 1 - Sw - Sg);
Scalar pC[numPhases];
Scalar kr[numPhases];
MaterialLaw::pC(pC, materialParams_, fs);
MaterialLaw::kr(kr, materialParams_, fs);
std::cerr << Sw << " "
<< Sg << " "
<< kr[0] << " "
<< kr[1] << " "
<< kr[2] << "\n";
}
std::cerr << "\n";
}
exit(1);
#endif
}
/*!
* \brief Returns the initial fluid state at a given a position.
*/
const CompositionalFluidState &
initialFluidState_(const GlobalPosition &pos) const
{
if (pos[2] < 2682.11)
return reservoirFluidState_[gPhaseIdx];
else {
if (2682.11 < pos[2] && pos[2] < 3032.61)
return reservoirFluidState_[oPhaseIdx];
else
return reservoirFluidState_[wPhaseIdx];
}
};
Scalar temperature_;
CompositionalFluidState reservoirFluidState_[numPhaseStates];
// // fluid state for gas injection
// MutableParameters injectGParams_;
// fluid state for water injection
CompositionalFluidState injectWFluidState_;
WellModel injectionWell_;
WellModel productionWell_;
static constexpr Scalar eps_ = 1e-6;
const Parser& parser_;
int episodeType_;
MaterialLawParams materialParams_;
ReservoirProperties* resProps_;
};
} //end namespace
#endif
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
