mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
move everything which is ECL specific to applications/ebos
this helps to keep the core blackoil model code lean and mean and it is also less confusing for newbies because the ECL blackoil simulator is not a "test" anymore. in case somebody wonders, "ebos" stands for "&eWoms &Black-&Oil &Simulator". I picked this name because it is short, a syllable, has not been taken by anything else (as far as I know) and "descriptive" names are rare for programs anyway: everyone who does not yet know about 'git' or 'emacs' and tells me that based on their names they must be a source-code managment system and an editor gets a crate of beer sponsored by me!
This commit is contained in:
Andreas Lauser
parent
1c35bb702f
commit
a8d5c72248
@ -1,984 +0,0 @@
|
|||||||
/*
|
|
||||||
Copyright (C) 2014 by Andreas Lauser
|
|
||||||
|
|
||||||
This file is part of the Open Porous Media project (OPM).
|
|
||||||
|
|
||||||
OPM 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.
|
|
||||||
|
|
||||||
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
/*!
|
|
||||||
* \file
|
|
||||||
*
|
|
||||||
* \copydoc Ewoms::EclProblem
|
|
||||||
*/
|
|
||||||
#ifndef EWOMS_ECL_PROBLEM_HH
|
|
||||||
#define EWOMS_ECL_PROBLEM_HH
|
|
||||||
|
|
||||||
#include <ewoms/models/blackoil/blackoilmodel.hh>
|
|
||||||
#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
|
|
||||||
#include <ewoms/io/eclgridmanager.hh>
|
|
||||||
#include <ewoms/io/eclipsesummarywriter.hh>
|
|
||||||
#include <ewoms/wells/eclwellmanager.hh>
|
|
||||||
|
|
||||||
#include <opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp>
|
|
||||||
#include <opm/material/fluidmatrixinteractions/SplineTwoPhaseMaterial.hpp>
|
|
||||||
#include <opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp>
|
|
||||||
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
|
|
||||||
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
|
|
||||||
|
|
||||||
#include <opm/core/utility/Average.hpp>
|
|
||||||
|
|
||||||
// for this simulator to make sense, dune-cornerpoint and opm-parser
|
|
||||||
// must be available
|
|
||||||
#include <dune/grid/CpGrid.hpp>
|
|
||||||
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
|
||||||
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
|
|
||||||
|
|
||||||
#include <dune/common/version.hh>
|
|
||||||
#include <dune/common/fvector.hh>
|
|
||||||
#include <dune/common/fmatrix.hh>
|
|
||||||
|
|
||||||
#include <boost/date_time.hpp>
|
|
||||||
|
|
||||||
#include <vector>
|
|
||||||
#include <string>
|
|
||||||
|
|
||||||
namespace Ewoms {
|
|
||||||
template <class TypeTag>
|
|
||||||
class EclProblem;
|
|
||||||
}
|
|
||||||
|
|
||||||
namespace Opm {
|
|
||||||
namespace Properties {
|
|
||||||
NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclGridManager));
|
|
||||||
|
|
||||||
// The temperature inside the reservoir
|
|
||||||
NEW_PROP_TAG(Temperature);
|
|
||||||
|
|
||||||
// Write all solutions for visualization, not just the ones for the
|
|
||||||
// report steps...
|
|
||||||
NEW_PROP_TAG(EnableWriteAllSolutions);
|
|
||||||
|
|
||||||
// Set the problem property
|
|
||||||
SET_TYPE_PROP(EclBaseProblem, Problem, Ewoms::EclProblem<TypeTag>);
|
|
||||||
|
|
||||||
// Select the element centered finite volume method as spatial discretization
|
|
||||||
SET_TAG_PROP(EclBaseProblem, SpatialDiscretizationSplice, EcfvDiscretization);
|
|
||||||
|
|
||||||
// Set the material Law
|
|
||||||
SET_PROP(EclBaseProblem, MaterialLaw)
|
|
||||||
{
|
|
||||||
private:
|
|
||||||
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
|
||||||
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
|
||||||
|
|
||||||
typedef Opm::TwoPhaseMaterialTraits<Scalar,
|
|
||||||
/*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
|
|
||||||
/*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx> OilWaterTraits;
|
|
||||||
|
|
||||||
typedef Opm::TwoPhaseMaterialTraits<Scalar,
|
|
||||||
/*wettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
|
|
||||||
/*nonWettingPhaseIdx=*/FluidSystem::gasPhaseIdx> GasOilTraits;
|
|
||||||
|
|
||||||
typedef Opm::ThreePhaseMaterialTraits<Scalar,
|
|
||||||
/*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
|
|
||||||
/*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
|
|
||||||
/*gasPhaseIdx=*/FluidSystem::gasPhaseIdx> Traits;
|
|
||||||
|
|
||||||
//typedef typename Opm::PiecewiseLinearTwoPhaseMaterial<OilWaterTraits> OilWaterLaw;
|
|
||||||
//typedef typename Opm::PiecewiseLinearTwoPhaseMaterial<GasOilTraits> GasOilLaw;
|
|
||||||
|
|
||||||
typedef typename Opm::SplineTwoPhaseMaterial<OilWaterTraits> OilWaterLaw;
|
|
||||||
typedef typename Opm::SplineTwoPhaseMaterial<GasOilTraits> GasOilLaw;
|
|
||||||
|
|
||||||
public:
|
|
||||||
typedef Opm::EclDefaultMaterial<Traits, GasOilLaw, OilWaterLaw> type;
|
|
||||||
};
|
|
||||||
|
|
||||||
// Enable gravity
|
|
||||||
SET_BOOL_PROP(EclBaseProblem, EnableGravity, true);
|
|
||||||
|
|
||||||
// Reuse the last linearization if possible?
|
|
||||||
SET_BOOL_PROP(EclBaseProblem, EnableLinearizationRecycling, true);
|
|
||||||
|
|
||||||
// Re-assemble the linearization only for the cells which have changed?
|
|
||||||
SET_BOOL_PROP(EclBaseProblem, EnablePartialRelinearization, true);
|
|
||||||
|
|
||||||
// only write the solutions for the report steps to disk
|
|
||||||
SET_BOOL_PROP(EclBaseProblem, EnableWriteAllSolutions, false);
|
|
||||||
|
|
||||||
// set the defaults for some problem specific properties
|
|
||||||
SET_SCALAR_PROP(EclBaseProblem, Temperature, 293.15);
|
|
||||||
|
|
||||||
// The default for the end time of the simulation [s]
|
|
||||||
//
|
|
||||||
// By default, stop it after the universe will probably have stopped
|
|
||||||
// to exist. (the ECL problem will finish the simulation explicitly
|
|
||||||
// after it simulated the last episode specified in the deck.)
|
|
||||||
SET_SCALAR_PROP(EclBaseProblem, EndTime, 1e100);
|
|
||||||
|
|
||||||
// The default for the initial time step size of the simulation [s].
|
|
||||||
//
|
|
||||||
// The chosen value means that the size of the first time step is the
|
|
||||||
// one of the initial episode (if the length of the initial episode is
|
|
||||||
// not millions of trillions of years, that is...)
|
|
||||||
SET_SCALAR_PROP(EclBaseProblem, InitialTimeStepSize, 1e100);
|
|
||||||
|
|
||||||
// Disable the VTK output by default for this problem ...
|
|
||||||
SET_BOOL_PROP(EclBaseProblem, EnableVtkOutput, false);
|
|
||||||
|
|
||||||
// ... but enable the Eclipse output by default
|
|
||||||
SET_BOOL_PROP(EclBaseProblem, EnableEclipseOutput, true);
|
|
||||||
|
|
||||||
// also enable the summary output.
|
|
||||||
SET_BOOL_PROP(EclBaseProblem, EnableEclipseSummaryOutput, true);
|
|
||||||
|
|
||||||
// The default DGF file to load
|
|
||||||
SET_STRING_PROP(EclBaseProblem, GridFile, "data/ecl.DATA");
|
|
||||||
}} // namespace Properties, Opm
|
|
||||||
|
|
||||||
namespace Ewoms {
|
|
||||||
/*!
|
|
||||||
* \ingroup TestProblems
|
|
||||||
*
|
|
||||||
* \brief This problem uses a deck in the format of the Eclipse
|
|
||||||
* simulator.
|
|
||||||
*/
|
|
||||||
template <class TypeTag>
|
|
||||||
class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
|
|
||||||
{
|
|
||||||
typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) 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;
|
|
||||||
|
|
||||||
// Grid and world dimension
|
|
||||||
enum { dim = GridView::dimension };
|
|
||||||
enum { dimWorld = GridView::dimensionworld };
|
|
||||||
|
|
||||||
// copy some indices for convenience
|
|
||||||
enum { numPhases = FluidSystem::numPhases };
|
|
||||||
enum { numComponents = FluidSystem::numComponents };
|
|
||||||
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
|
|
||||||
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
|
|
||||||
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
|
|
||||||
enum { gasCompIdx = FluidSystem::gasCompIdx };
|
|
||||||
enum { oilCompIdx = FluidSystem::oilCompIdx };
|
|
||||||
enum { waterCompIdx = FluidSystem::waterCompIdx };
|
|
||||||
|
|
||||||
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
|
|
||||||
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
|
|
||||||
typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
|
|
||||||
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
|
|
||||||
typedef typename GET_PROP_TYPE(TypeTag, BlackOilFluidState) BlackOilFluidState;
|
|
||||||
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
|
|
||||||
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
|
|
||||||
|
|
||||||
typedef Ewoms::EclipseSummaryWriter<TypeTag> EclipseSummaryWriter;
|
|
||||||
|
|
||||||
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
|
|
||||||
|
|
||||||
public:
|
|
||||||
/*!
|
|
||||||
* \copydoc FvBaseProblem::registerParameters
|
|
||||||
*/
|
|
||||||
static void registerParameters()
|
|
||||||
{
|
|
||||||
ParentType::registerParameters();
|
|
||||||
|
|
||||||
EWOMS_REGISTER_PARAM(TypeTag, Scalar, Temperature,
|
|
||||||
"The temperature [K] in the reservoir");
|
|
||||||
EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWriteAllSolutions,
|
|
||||||
"Write all solutions to disk instead of only the ones for the "
|
|
||||||
"report steps");
|
|
||||||
}
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \copydoc Doxygen::defaultProblemConstructor
|
|
||||||
*/
|
|
||||||
EclProblem(Simulator &simulator)
|
|
||||||
: ParentType(simulator)
|
|
||||||
, wellManager_(simulator)
|
|
||||||
, summaryWriter_(simulator)
|
|
||||||
{ }
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \copydoc FvBaseProblem::finishInit
|
|
||||||
*/
|
|
||||||
void finishInit()
|
|
||||||
{
|
|
||||||
ParentType::finishInit();
|
|
||||||
|
|
||||||
auto& simulator = this->simulator();
|
|
||||||
|
|
||||||
temperature_ = EWOMS_GET_PARAM(TypeTag, Scalar, Temperature);
|
|
||||||
|
|
||||||
// invert the direction of the gravity vector for ECL problems
|
|
||||||
// (z coodinates represent depth, not height.)
|
|
||||||
this->gravity_[dim - 1] *= -1;
|
|
||||||
|
|
||||||
initFluidSystem_();
|
|
||||||
readMaterialParameters_();
|
|
||||||
readInitialCondition_();
|
|
||||||
|
|
||||||
// initialize the wells. Note that this needs to be done after initializing the
|
|
||||||
// intrinsic permeabilities because the well model uses them...
|
|
||||||
wellManager_.init(simulator.gridManager().eclipseState());
|
|
||||||
|
|
||||||
// Start the first episode. For this, ask the Eclipse schedule.
|
|
||||||
Opm::TimeMapConstPtr timeMap = simulator.gridManager().schedule()->getTimeMap();
|
|
||||||
tm curTime = boost::posix_time::to_tm(timeMap->getStartTime(/*timeStepIdx=*/0));
|
|
||||||
|
|
||||||
Scalar startTime = std::mktime(&curTime);
|
|
||||||
simulator.setStartTime(startTime);
|
|
||||||
simulator.startNextEpisode(/*startTime=*/startTime,
|
|
||||||
/*length=*/timeMap->getTimeStepLength(/*timeStepIdx=*/0));
|
|
||||||
|
|
||||||
// we want the episode index to be the same as the report step
|
|
||||||
// index to make things simpler...
|
|
||||||
simulator.setEpisodeIndex(0);
|
|
||||||
|
|
||||||
// the user-specified initial time step can be shorter than
|
|
||||||
// the initial report step size given in the deck, but it
|
|
||||||
// can't be longer.
|
|
||||||
Scalar dt = simulator.timeStepSize();
|
|
||||||
if (dt > simulator.episodeLength())
|
|
||||||
simulator.setTimeStepSize(simulator.episodeLength());
|
|
||||||
}
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \brief Called by the simulator before an episode begins.
|
|
||||||
*/
|
|
||||||
void beginEpisode()
|
|
||||||
{ wellManager_.beginEpisode(this->simulator().gridManager().eclipseState()); }
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \brief Called by the simulator before each time integration.
|
|
||||||
*/
|
|
||||||
void beginTimeStep()
|
|
||||||
{ wellManager_.beginTimeStep(); }
|
|
||||||
/*!
|
|
||||||
* \brief Called by the simulator before each Newton-Raphson iteration.
|
|
||||||
*/
|
|
||||||
void beginIteration()
|
|
||||||
{ wellManager_.beginIteration(); }
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \brief Called by the simulator after each Newton-Raphson iteration.
|
|
||||||
*/
|
|
||||||
void endIteration()
|
|
||||||
{ wellManager_.endIteration(); }
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \brief Called by the simulator after each time integration.
|
|
||||||
*/
|
|
||||||
void endTimeStep()
|
|
||||||
{
|
|
||||||
wellManager_.endTimeStep();
|
|
||||||
|
|
||||||
#ifndef NDEBUG
|
|
||||||
this->model().checkConservativeness(/*tolerance=*/-1, /*verbose=*/true);
|
|
||||||
#endif // NDEBUG
|
|
||||||
}
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \brief Called by the simulator after the end of an episode.
|
|
||||||
*/
|
|
||||||
void endEpisode()
|
|
||||||
{
|
|
||||||
// first, write the summary information ...
|
|
||||||
summaryWriter_.write(wellManager_);
|
|
||||||
|
|
||||||
// ... then proceed to the next report step
|
|
||||||
Simulator &simulator = this->simulator();
|
|
||||||
Opm::EclipseStateConstPtr eclipseState = this->simulator().gridManager().eclipseState();
|
|
||||||
Opm::TimeMapConstPtr timeMap = eclipseState->getSchedule()->getTimeMap();
|
|
||||||
|
|
||||||
// TimeMap deals with points in time, so the number of time
|
|
||||||
// intervals (i.e., report steps) is one less!
|
|
||||||
int numReportSteps = timeMap->size() - 1;
|
|
||||||
|
|
||||||
// start the next episode if there are additional report
|
|
||||||
// steps, else finish the simulation
|
|
||||||
int nextEpisodeIdx = simulator.episodeIndex() + 1;
|
|
||||||
if (nextEpisodeIdx < numReportSteps) {
|
|
||||||
simulator.startNextEpisode(timeMap->getTimeStepLength(nextEpisodeIdx));
|
|
||||||
simulator.setTimeStepSize(timeMap->getTimeStepLength(nextEpisodeIdx));
|
|
||||||
}
|
|
||||||
else
|
|
||||||
simulator.setFinished(true);
|
|
||||||
}
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \brief Returns true if the current solution should be written
|
|
||||||
* to disk for visualization.
|
|
||||||
*
|
|
||||||
* For the ECL simulator we only write at the end of
|
|
||||||
* episodes/report steps...
|
|
||||||
*/
|
|
||||||
bool shouldWriteOutput()
|
|
||||||
{
|
|
||||||
if (this->simulator().timeStepIndex() < 0)
|
|
||||||
// always write the initial solution
|
|
||||||
return true;
|
|
||||||
|
|
||||||
if (EWOMS_GET_PARAM(TypeTag, bool, EnableWriteAllSolutions))
|
|
||||||
return true;
|
|
||||||
|
|
||||||
return this->simulator().episodeWillBeOver();
|
|
||||||
}
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \copydoc FvBaseMultiPhaseProblem::intrinsicPermeability
|
|
||||||
*/
|
|
||||||
template <class Context>
|
|
||||||
const DimMatrix &intrinsicPermeability(const Context &context,
|
|
||||||
int spaceIdx,
|
|
||||||
int timeIdx) const
|
|
||||||
{
|
|
||||||
int globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
|
|
||||||
return intrinsicPermeability_[globalSpaceIdx];
|
|
||||||
}
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \copydoc FvBaseMultiPhaseProblem::intersectionIntrinsicPermeability
|
|
||||||
*/
|
|
||||||
template <class Context>
|
|
||||||
void intersectionIntrinsicPermeability(DimMatrix &result,
|
|
||||||
const Context &context,
|
|
||||||
int localIntersectionIdx, int timeIdx) const
|
|
||||||
{
|
|
||||||
// calculate the intersection index
|
|
||||||
const auto &scvf = context.stencil(timeIdx).interiorFace(localIntersectionIdx);
|
|
||||||
|
|
||||||
int numElements = this->model().numGridDof();
|
|
||||||
|
|
||||||
size_t interiorElemIdx = context.globalSpaceIndex(scvf.interiorIndex(), timeIdx);
|
|
||||||
size_t exteriorElemIdx = context.globalSpaceIndex(scvf.exteriorIndex(), timeIdx);
|
|
||||||
|
|
||||||
size_t elem1Idx = std::min(interiorElemIdx, exteriorElemIdx);
|
|
||||||
size_t elem2Idx = std::max(interiorElemIdx, exteriorElemIdx);
|
|
||||||
|
|
||||||
size_t globalIntersectionIdx = elem1Idx*numElements + elem2Idx;
|
|
||||||
result = intersectionIntrinsicPermeability_.at(globalIntersectionIdx);
|
|
||||||
}
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \copydoc FvBaseMultiPhaseProblem::porosity
|
|
||||||
*/
|
|
||||||
template <class Context>
|
|
||||||
Scalar porosity(const Context &context, int spaceIdx, int timeIdx) const
|
|
||||||
{
|
|
||||||
int globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
|
|
||||||
return porosity_[globalSpaceIdx];
|
|
||||||
}
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \copydoc FvBaseMultiPhaseProblem::materialLawParams
|
|
||||||
*/
|
|
||||||
template <class Context>
|
|
||||||
const MaterialLawParams &materialLawParams(const Context &context,
|
|
||||||
int spaceIdx, int timeIdx) const
|
|
||||||
{
|
|
||||||
int tableIdx = 0;
|
|
||||||
if (materialParamTableIdx_.size() > 0) {
|
|
||||||
int globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
|
|
||||||
tableIdx = materialParamTableIdx_[globalSpaceIdx];
|
|
||||||
}
|
|
||||||
return materialParams_[tableIdx];
|
|
||||||
}
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \brief Returns the index of the relevant region for thermodynmic properties
|
|
||||||
*/
|
|
||||||
template <class Context>
|
|
||||||
int pvtRegionIndex(const Context &context, int spaceIdx, int timeIdx) const
|
|
||||||
{
|
|
||||||
Opm::DeckConstPtr deck = this->simulator().gridManager().deck();
|
|
||||||
|
|
||||||
if (!deck->hasKeyword("PVTNUM"))
|
|
||||||
return 0;
|
|
||||||
|
|
||||||
const auto &grid = this->simulator().gridManager().grid();
|
|
||||||
|
|
||||||
// this is quite specific to the ECFV discretization. But so is everything in an
|
|
||||||
// ECL deck, i.e., we don't need to care here...
|
|
||||||
int compressedDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
|
|
||||||
int cartesianDofIdx = grid.globalCell()[compressedDofIdx];
|
|
||||||
|
|
||||||
return deck->getKeyword("PVTNUM")->getIntData()[cartesianDofIdx] - 1;
|
|
||||||
}
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \name Problem parameters
|
|
||||||
*/
|
|
||||||
//! \{
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \copydoc FvBaseProblem::name
|
|
||||||
*/
|
|
||||||
std::string name() const
|
|
||||||
{ return this->simulator().gridManager().caseName(); }
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \copydoc FvBaseMultiPhaseProblem::temperature
|
|
||||||
*
|
|
||||||
* The black-oil model assumes constant temperature to define its
|
|
||||||
* parameters. Although temperature is thus not really used by the
|
|
||||||
* model, it gets written to the VTK output. Who nows, maybe we
|
|
||||||
* will need it one day?
|
|
||||||
*/
|
|
||||||
template <class Context>
|
|
||||||
Scalar temperature(const Context &context, int spaceIdx, int timeIdx) const
|
|
||||||
{ return temperature_; }
|
|
||||||
|
|
||||||
// \}
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \name Boundary conditions
|
|
||||||
*/
|
|
||||||
//! \{
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \copydoc FvBaseProblem::boundary
|
|
||||||
*
|
|
||||||
* Eclipse uses no-flow conditions for all boundaries. \todo really?
|
|
||||||
*/
|
|
||||||
template <class Context>
|
|
||||||
void boundary(BoundaryRateVector &values,
|
|
||||||
const Context &context,
|
|
||||||
int spaceIdx,
|
|
||||||
int timeIdx) const
|
|
||||||
{ values.setNoFlow(); }
|
|
||||||
|
|
||||||
//! \}
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \name Volumetric terms
|
|
||||||
*/
|
|
||||||
//! \{
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \copydoc FvBaseProblem::initial
|
|
||||||
*
|
|
||||||
* The reservoir problem uses a constant boundary condition for
|
|
||||||
* the whole domain.
|
|
||||||
*/
|
|
||||||
template <class Context>
|
|
||||||
void initial(PrimaryVariables &values, const Context &context, int spaceIdx, int timeIdx) const
|
|
||||||
{
|
|
||||||
int globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
|
|
||||||
|
|
||||||
values.assignNaive(initialFluidStates_[globalDofIdx]);
|
|
||||||
}
|
|
||||||
|
|
||||||
/*!
|
|
||||||
* \copydoc FvBaseProblem::source
|
|
||||||
*
|
|
||||||
* For this problem, the source term of all components is 0 everywhere.
|
|
||||||
*/
|
|
||||||
template <class Context>
|
|
||||||
void source(RateVector &rate,
|
|
||||||
const Context &context,
|
|
||||||
int spaceIdx,
|
|
||||||
int timeIdx) const
|
|
||||||
{
|
|
||||||
rate = 0;
|
|
||||||
wellManager_.computeTotalRatesForDof(rate, context, spaceIdx, timeIdx);
|
|
||||||
|
|
||||||
// convert the source term from the total mass rate of the
|
|
||||||
// cell to the one per unit of volume as used by the model.
|
|
||||||
int globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
|
|
||||||
rate /= this->model().dofTotalVolume(globalDofIdx);
|
|
||||||
}
|
|
||||||
|
|
||||||
//! \}
|
|
||||||
|
|
||||||
private:
|
|
||||||
void readMaterialParameters_()
|
|
||||||
{
|
|
||||||
auto deck = this->simulator().gridManager().deck();
|
|
||||||
auto eclipseState = this->simulator().gridManager().eclipseState();
|
|
||||||
const auto &grid = this->simulator().gridManager().grid();
|
|
||||||
|
|
||||||
size_t numDof = this->model().numGridDof();
|
|
||||||
|
|
||||||
intrinsicPermeability_.resize(numDof);
|
|
||||||
porosity_.resize(numDof);
|
|
||||||
materialParams_.resize(numDof);
|
|
||||||
|
|
||||||
////////////////////////////////
|
|
||||||
// permeability
|
|
||||||
|
|
||||||
// read the intrinsic permeabilities from the eclipseState. Note that all arrays
|
|
||||||
// provided by eclipseState are one-per-cell of "uncompressed" grid, whereas the
|
|
||||||
// dune-cornerpoint grid object might remove a few elements...
|
|
||||||
if (eclipseState->hasDoubleGridProperty("PERMX")) {
|
|
||||||
const std::vector<double> &permxData =
|
|
||||||
eclipseState->getDoubleGridProperty("PERMX")->getData();
|
|
||||||
std::vector<double> permyData(permxData);
|
|
||||||
if (eclipseState->hasDoubleGridProperty("PERMY"))
|
|
||||||
permyData = eclipseState->getDoubleGridProperty("PERMY")->getData();
|
|
||||||
std::vector<double> permzData(permxData);
|
|
||||||
if (eclipseState->hasDoubleGridProperty("PERMZ"))
|
|
||||||
permzData = eclipseState->getDoubleGridProperty("PERMZ")->getData();
|
|
||||||
|
|
||||||
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
|
|
||||||
int cartesianElemIdx = grid.globalCell()[dofIdx];
|
|
||||||
intrinsicPermeability_[dofIdx] = 0.0;
|
|
||||||
intrinsicPermeability_[dofIdx][0][0] = permxData[cartesianElemIdx];
|
|
||||||
intrinsicPermeability_[dofIdx][1][1] = permyData[cartesianElemIdx];
|
|
||||||
intrinsicPermeability_[dofIdx][2][2] = permzData[cartesianElemIdx];
|
|
||||||
}
|
|
||||||
|
|
||||||
// for now we don't care about non-diagonal entries
|
|
||||||
}
|
|
||||||
else
|
|
||||||
OPM_THROW(std::logic_error,
|
|
||||||
"Can't read the intrinsic permeability from the eclipse state. "
|
|
||||||
"(The PERM{X,Y,Z} keywords are missing)");
|
|
||||||
|
|
||||||
// apply the NTG keyword to the X and Y permeabilities
|
|
||||||
if (eclipseState->hasDoubleGridProperty("NTG")) {
|
|
||||||
const std::vector<double> &ntgData =
|
|
||||||
eclipseState->getDoubleGridProperty("NTG")->getData();
|
|
||||||
|
|
||||||
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
|
|
||||||
int cartesianElemIdx = grid.globalCell()[dofIdx];
|
|
||||||
intrinsicPermeability_[dofIdx][0][0] *= ntgData[cartesianElemIdx];
|
|
||||||
intrinsicPermeability_[dofIdx][1][1] *= ntgData[cartesianElemIdx];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
computeFaceIntrinsicPermeabilities_();
|
|
||||||
////////////////////////////////
|
|
||||||
|
|
||||||
|
|
||||||
////////////////////////////////
|
|
||||||
// compute the porosity
|
|
||||||
if (eclipseState->hasDoubleGridProperty("PORO")) {
|
|
||||||
const std::vector<double> &poroData =
|
|
||||||
eclipseState->getDoubleGridProperty("PORO")->getData();
|
|
||||||
|
|
||||||
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
|
|
||||||
int cartesianElemIdx = grid.globalCell()[dofIdx];
|
|
||||||
porosity_[dofIdx] = poroData[cartesianElemIdx];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
else
|
|
||||||
OPM_THROW(std::logic_error,
|
|
||||||
"Can't read the porosity from the eclipse state. "
|
|
||||||
"(The PORO keyword is missing)");
|
|
||||||
|
|
||||||
// apply the NTG keyword to the porosity
|
|
||||||
if (eclipseState->hasDoubleGridProperty("NTG")) {
|
|
||||||
const std::vector<double> &ntgData =
|
|
||||||
eclipseState->getDoubleGridProperty("NTG")->getData();
|
|
||||||
|
|
||||||
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
|
|
||||||
int cartesianElemIdx = grid.globalCell()[dofIdx];
|
|
||||||
porosity_[dofIdx] *= ntgData[cartesianElemIdx];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// apply the MULTPV keyword to the porosity
|
|
||||||
if (eclipseState->hasDoubleGridProperty("MULTPV")) {
|
|
||||||
const std::vector<double> &multpvData =
|
|
||||||
eclipseState->getDoubleGridProperty("MULTPV")->getData();
|
|
||||||
|
|
||||||
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
|
|
||||||
int cartesianElemIdx = grid.globalCell()[dofIdx];
|
|
||||||
porosity_[dofIdx] *= multpvData[cartesianElemIdx];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
////////////////////////////////
|
|
||||||
// fluid parameters
|
|
||||||
const auto& swofTables = eclipseState->getSwofTables();
|
|
||||||
const auto& sgofTables = eclipseState->getSgofTables();
|
|
||||||
|
|
||||||
// the number of tables for the SWOF and the SGOF keywords
|
|
||||||
// must be identical
|
|
||||||
assert(swofTables.size() == sgofTables.size());
|
|
||||||
|
|
||||||
size_t numSatfuncTables = swofTables.size();
|
|
||||||
materialParams_.resize(numSatfuncTables);
|
|
||||||
|
|
||||||
typedef typename MaterialLawParams::GasOilParams GasOilParams;
|
|
||||||
typedef typename MaterialLawParams::OilWaterParams OilWaterParams;
|
|
||||||
|
|
||||||
for (size_t tableIdx = 0; tableIdx < numSatfuncTables; ++ tableIdx) {
|
|
||||||
// set the parameters of the material law for a given table
|
|
||||||
OilWaterParams owParams;
|
|
||||||
GasOilParams goParams;
|
|
||||||
|
|
||||||
const auto& swofTable = swofTables[tableIdx];
|
|
||||||
const auto& sgofTable = sgofTables[tableIdx];
|
|
||||||
|
|
||||||
const auto &SwColumn = swofTable.getSwColumn();
|
|
||||||
|
|
||||||
owParams.setKrwSamples(SwColumn, swofTable.getKrwColumn());
|
|
||||||
owParams.setKrnSamples(SwColumn, swofTable.getKrowColumn());
|
|
||||||
owParams.setPcnwSamples(SwColumn, swofTable.getPcowColumn());
|
|
||||||
|
|
||||||
// convert the saturations of the SGOF keyword from gas to oil saturations
|
|
||||||
std::vector<double> SoSamples(sgofTable.numRows());
|
|
||||||
for (size_t sampleIdx = 0; sampleIdx < sgofTable.numRows(); ++ sampleIdx)
|
|
||||||
SoSamples[sampleIdx] = 1 - sgofTable.getSgColumn()[sampleIdx];
|
|
||||||
|
|
||||||
goParams.setKrwSamples(SoSamples, sgofTable.getKrogColumn());
|
|
||||||
goParams.setKrnSamples(SoSamples, sgofTable.getKrgColumn());
|
|
||||||
goParams.setPcnwSamples(SoSamples, sgofTable.getPcogColumn());
|
|
||||||
|
|
||||||
owParams.finalize();
|
|
||||||
goParams.finalize();
|
|
||||||
|
|
||||||
// compute the connate water saturation. In Eclipse decks that is defined as
|
|
||||||
// the first saturation value of the SWOF keyword.
|
|
||||||
Scalar Swco = SwColumn.front();
|
|
||||||
materialParams_[tableIdx].setConnateWaterSaturation(Swco);
|
|
||||||
|
|
||||||
materialParams_[tableIdx].setOilWaterParams(owParams);
|
|
||||||
materialParams_[tableIdx].setGasOilParams(goParams);
|
|
||||||
|
|
||||||
materialParams_[tableIdx].finalize();
|
|
||||||
}
|
|
||||||
|
|
||||||
// set the index of the table to be used
|
|
||||||
if (eclipseState->hasIntGridProperty("SATNUM")) {
|
|
||||||
const std::vector<int> &satnumData =
|
|
||||||
eclipseState->getIntGridProperty("SATNUM")->getData();
|
|
||||||
|
|
||||||
materialParamTableIdx_.resize(numDof);
|
|
||||||
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
|
|
||||||
int cartesianElemIdx = grid.globalCell()[dofIdx];
|
|
||||||
|
|
||||||
// make sure that all values are in the correct range
|
|
||||||
assert(1 <= satnumData[dofIdx]);
|
|
||||||
assert(satnumData[dofIdx] <= static_cast<int>(numSatfuncTables));
|
|
||||||
|
|
||||||
// Eclipse uses Fortran-style indices which start at
|
|
||||||
// 1, but this here is C++...
|
|
||||||
materialParamTableIdx_[dofIdx] = satnumData[cartesianElemIdx] - 1;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
else
|
|
||||||
materialParamTableIdx_.clear();
|
|
||||||
////////////////////////////////
|
|
||||||
}
|
|
||||||
|
|
||||||
void initFluidSystem_()
|
|
||||||
{
|
|
||||||
const auto deck = this->simulator().gridManager().deck();
|
|
||||||
const auto eclipseState = this->simulator().gridManager().eclipseState();
|
|
||||||
|
|
||||||
FluidSystem::initBegin();
|
|
||||||
|
|
||||||
int numRegions = deck->getKeyword("DENSITY")->size();
|
|
||||||
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
|
|
||||||
// set the reference densities
|
|
||||||
Opm::DeckRecordConstPtr densityRecord =
|
|
||||||
deck->getKeyword("DENSITY")->getRecord(regionIdx);
|
|
||||||
FluidSystem::setReferenceDensities(densityRecord->getItem("OIL")->getSIDouble(0),
|
|
||||||
densityRecord->getItem("WATER")->getSIDouble(0),
|
|
||||||
densityRecord->getItem("GAS")->getSIDouble(0),
|
|
||||||
regionIdx);
|
|
||||||
|
|
||||||
// so far, we require the presence of the PVTO, PVTW and PVDG
|
|
||||||
// keywords...
|
|
||||||
FluidSystem::setPvtoTable(eclipseState->getPvtoTables()[regionIdx], regionIdx);
|
|
||||||
FluidSystem::setPvtw(deck->getKeyword("PVTW"), regionIdx);
|
|
||||||
FluidSystem::setPvdgTable(eclipseState->getPvdgTables()[regionIdx], regionIdx);
|
|
||||||
}
|
|
||||||
|
|
||||||
FluidSystem::initEnd();
|
|
||||||
}
|
|
||||||
|
|
||||||
void readInitialCondition_()
|
|
||||||
{
|
|
||||||
const auto deck = this->simulator().gridManager().deck();
|
|
||||||
const auto &grid = this->simulator().gridManager().grid();
|
|
||||||
|
|
||||||
if (!deck->hasKeyword("SWAT") ||
|
|
||||||
!deck->hasKeyword("SGAS"))
|
|
||||||
OPM_THROW(std::runtime_error,
|
|
||||||
"So far, the Eclipse input file requires the presence of the SWAT "
|
|
||||||
"and SGAS keywords");
|
|
||||||
if (!deck->hasKeyword("PRESSURE"))
|
|
||||||
OPM_THROW(std::runtime_error,
|
|
||||||
"So far, the Eclipse input file requires the presence of the PRESSURE "
|
|
||||||
"keyword");
|
|
||||||
if (!deck->hasKeyword("DISGAS"))
|
|
||||||
OPM_THROW(std::runtime_error,
|
|
||||||
"The deck must exhibit gas dissolved in the oil phase"
|
|
||||||
" (DISGAS keyword is missing)");
|
|
||||||
if (!deck->hasKeyword("RS"))
|
|
||||||
OPM_THROW(std::runtime_error,
|
|
||||||
"The Eclipse input file requires the presence of the RS keyword");
|
|
||||||
|
|
||||||
if (deck->hasKeyword("VAPOIL"))
|
|
||||||
OPM_THROW(std::runtime_error,
|
|
||||||
"The deck must _not_ exhibit vaporized oil"
|
|
||||||
" (The VAPOIL keyword is unsupported)");
|
|
||||||
if (deck->hasKeyword("RV"))
|
|
||||||
OPM_THROW(std::runtime_error,
|
|
||||||
"The Eclipse input file requires the RV keyword to be non-present");
|
|
||||||
|
|
||||||
size_t numDof = this->model().numGridDof();
|
|
||||||
|
|
||||||
initialFluidStates_.resize(numDof);
|
|
||||||
|
|
||||||
const std::vector<double> &waterSaturationData =
|
|
||||||
deck->getKeyword("SWAT")->getSIDoubleData();
|
|
||||||
const std::vector<double> &gasSaturationData =
|
|
||||||
deck->getKeyword("SGAS")->getSIDoubleData();
|
|
||||||
const std::vector<double> &pressureData =
|
|
||||||
deck->getKeyword("PRESSURE")->getSIDoubleData();
|
|
||||||
const std::vector<double> &rsData =
|
|
||||||
deck->getKeyword("RS")->getSIDoubleData();
|
|
||||||
|
|
||||||
// make sure that the size of the data arrays is correct
|
|
||||||
#ifndef NDEBUG
|
|
||||||
const auto &cartSize = grid.logicalCartesianSize();
|
|
||||||
size_t numCartesianCells = cartSize[0] * cartSize[1] * cartSize[2];
|
|
||||||
assert(waterSaturationData.size() == numCartesianCells);
|
|
||||||
assert(gasSaturationData.size() == numCartesianCells);
|
|
||||||
assert(pressureData.size() == numCartesianCells);
|
|
||||||
assert(rsData.size() == numCartesianCells);
|
|
||||||
#endif
|
|
||||||
|
|
||||||
// calculate the initial fluid states
|
|
||||||
for (size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
|
|
||||||
auto &dofFluidState = initialFluidStates_[dofIdx];
|
|
||||||
|
|
||||||
size_t cartesianDofIdx = grid.globalCell()[dofIdx];
|
|
||||||
assert(0 <= cartesianDofIdx);
|
|
||||||
assert(cartesianDofIdx <= numCartesianCells);
|
|
||||||
|
|
||||||
//////
|
|
||||||
// set temperatures
|
|
||||||
//////
|
|
||||||
dofFluidState.setTemperature(temperature_);
|
|
||||||
|
|
||||||
//////
|
|
||||||
// set saturations
|
|
||||||
//////
|
|
||||||
dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
|
|
||||||
waterSaturationData[cartesianDofIdx]);
|
|
||||||
dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
|
|
||||||
gasSaturationData[cartesianDofIdx]);
|
|
||||||
dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
|
|
||||||
1
|
|
||||||
- waterSaturationData[cartesianDofIdx]
|
|
||||||
- gasSaturationData[cartesianDofIdx]);
|
|
||||||
|
|
||||||
//////
|
|
||||||
// set pressures
|
|
||||||
//////
|
|
||||||
Scalar oilPressure = pressureData[cartesianDofIdx];
|
|
||||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
|
||||||
dofFluidState.setPressure(phaseIdx, oilPressure);
|
|
||||||
}
|
|
||||||
|
|
||||||
//////
|
|
||||||
// set compositions
|
|
||||||
//////
|
|
||||||
|
|
||||||
// reset all mole fractions to 0
|
|
||||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
|
|
||||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
|
|
||||||
dofFluidState.setMoleFraction(phaseIdx, compIdx, 0.0);
|
|
||||||
|
|
||||||
// set compositions of the gas and water phases
|
|
||||||
dofFluidState.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
|
|
||||||
dofFluidState.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
|
|
||||||
|
|
||||||
|
|
||||||
// set the composition of the oil phase:
|
|
||||||
//
|
|
||||||
// first, retrieve the relevant black-oil parameters from
|
|
||||||
// the fluid system.
|
|
||||||
Scalar RsSat = FluidSystem::gasDissolutionFactor(oilPressure, /*regionIdx=*/0);
|
|
||||||
Scalar RsReal = rsData[cartesianDofIdx];
|
|
||||||
|
|
||||||
if (RsReal > RsSat) {
|
|
||||||
std::array<int, 3> ijk;
|
|
||||||
grid.getIJK(dofIdx, ijk);
|
|
||||||
std::cerr << "Warning: The specified amount gas (R_s = " << RsReal << ") is more"
|
|
||||||
<< " than the maximium\n"
|
|
||||||
<< " amount which can be dissolved in oil"
|
|
||||||
<< " (R_s,max=" << RsSat << ")"
|
|
||||||
<< " for cell (" << ijk[0] << ", " << ijk[1] << ", " << ijk[2] << ")."
|
|
||||||
<< " Ignoring.\n";
|
|
||||||
RsReal = RsSat;
|
|
||||||
}
|
|
||||||
|
|
||||||
// calculate composition of the real and the saturated oil phase in terms of
|
|
||||||
// mass fractions.
|
|
||||||
Scalar rhooRef = FluidSystem::referenceDensity(oilPhaseIdx, /*regionIdx=*/0);
|
|
||||||
Scalar rhogRef = FluidSystem::referenceDensity(gasPhaseIdx, /*regionIdx=*/0);
|
|
||||||
Scalar XoGReal = RsReal*rhogRef / (RsReal*rhogRef + rhooRef);
|
|
||||||
|
|
||||||
// convert mass to mole fractions
|
|
||||||
Scalar MG = FluidSystem::molarMass(gasCompIdx);
|
|
||||||
Scalar MO = FluidSystem::molarMass(oilCompIdx);
|
|
||||||
|
|
||||||
Scalar xoGReal = XoGReal * MO / ((MO - MG) * XoGReal + MG);
|
|
||||||
Scalar xoOReal = 1 - xoGReal;
|
|
||||||
|
|
||||||
// finally, set the oil-phase composition
|
|
||||||
dofFluidState.setMoleFraction(oilPhaseIdx, gasCompIdx, xoGReal);
|
|
||||||
dofFluidState.setMoleFraction(oilPhaseIdx, oilCompIdx, xoOReal);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
void computeFaceIntrinsicPermeabilities_()
|
|
||||||
{
|
|
||||||
auto eclipseState = this->simulator().gridManager().eclipseState();
|
|
||||||
const auto &grid = this->simulator().gridManager().grid();
|
|
||||||
|
|
||||||
int numElements = this->gridView().size(/*codim=*/0);
|
|
||||||
|
|
||||||
std::vector<double> multx(numElements, 1.0);
|
|
||||||
std::vector<double> multy(numElements, 1.0);
|
|
||||||
std::vector<double> multz(numElements, 1.0);
|
|
||||||
std::vector<double> multxMinus(numElements, 1.0);
|
|
||||||
std::vector<double> multyMinus(numElements, 1.0);
|
|
||||||
std::vector<double> multzMinus(numElements, 1.0);
|
|
||||||
|
|
||||||
// retrieve the transmissibility multiplier keywords. Note that we use them as
|
|
||||||
// permeability multipliers...
|
|
||||||
if (eclipseState->hasDoubleGridProperty("MULTX"))
|
|
||||||
multx = eclipseState->getDoubleGridProperty("MULTX")->getData();
|
|
||||||
if (eclipseState->hasDoubleGridProperty("MULTX-"))
|
|
||||||
multxMinus = eclipseState->getDoubleGridProperty("MULTX-")->getData();
|
|
||||||
if (eclipseState->hasDoubleGridProperty("MULTY"))
|
|
||||||
multy = eclipseState->getDoubleGridProperty("MULTY")->getData();
|
|
||||||
if (eclipseState->hasDoubleGridProperty("MULTY-"))
|
|
||||||
multyMinus = eclipseState->getDoubleGridProperty("MULTY-")->getData();
|
|
||||||
if (eclipseState->hasDoubleGridProperty("MULTZ"))
|
|
||||||
multz = eclipseState->getDoubleGridProperty("MULTZ")->getData();
|
|
||||||
if (eclipseState->hasDoubleGridProperty("MULTZ-"))
|
|
||||||
multzMinus = eclipseState->getDoubleGridProperty("MULTZ-")->getData();
|
|
||||||
|
|
||||||
// resize the hash map to a appropriate size for a conforming 3D grid
|
|
||||||
float maxLoadFactor = intersectionIntrinsicPermeability_.max_load_factor();
|
|
||||||
intersectionIntrinsicPermeability_.reserve(numElements * 6 / maxLoadFactor * 1.05 );
|
|
||||||
|
|
||||||
auto elemIt = this->gridView().template begin</*codim=*/0>();
|
|
||||||
const auto& elemEndIt = this->gridView().template end</*codim=*/0>();
|
|
||||||
for (; elemIt != elemEndIt; ++elemIt) {
|
|
||||||
auto intersectIt = elemIt->ileafbegin();
|
|
||||||
const auto &intersectEndIt = elemIt->ileafend();
|
|
||||||
for (; intersectIt != intersectEndIt; ++intersectIt) {
|
|
||||||
if (!intersectIt->neighbor())
|
|
||||||
// skip boundary intersections...
|
|
||||||
continue;
|
|
||||||
|
|
||||||
// calculate the "intersection index"
|
|
||||||
size_t interiorElemIdx = this->elementMapper().map(intersectIt->inside());
|
|
||||||
size_t exteriorElemIdx = this->elementMapper().map(intersectIt->outside());
|
|
||||||
|
|
||||||
size_t elem1Idx = std::min(interiorElemIdx, exteriorElemIdx);
|
|
||||||
size_t elem2Idx = std::max(interiorElemIdx, exteriorElemIdx);
|
|
||||||
|
|
||||||
size_t intersectIdx = elem1Idx*numElements + elem2Idx;
|
|
||||||
|
|
||||||
// do nothing if this intersection was already seen "from the other side"
|
|
||||||
if (intersectionIntrinsicPermeability_.count(intersectIdx) > 0)
|
|
||||||
continue;
|
|
||||||
|
|
||||||
auto K1 = intrinsicPermeability_[interiorElemIdx];
|
|
||||||
auto K2 = intrinsicPermeability_[exteriorElemIdx];
|
|
||||||
|
|
||||||
int interiorElemCartIdx = grid.globalCell()[interiorElemIdx];
|
|
||||||
int exteriorElemCartIdx = grid.globalCell()[exteriorElemIdx];
|
|
||||||
|
|
||||||
// local index of the face of the interior element which contains the
|
|
||||||
// intersection
|
|
||||||
int insideFaceIdx = intersectIt->indexInInside();
|
|
||||||
|
|
||||||
// take the transmissibility multipliers into account (i.e., the
|
|
||||||
// MULT[XYZ]-? keywords)
|
|
||||||
if (insideFaceIdx == 1) { // right
|
|
||||||
K1 *= multx[interiorElemCartIdx];
|
|
||||||
K2 *= multxMinus[exteriorElemCartIdx];
|
|
||||||
}
|
|
||||||
else if (insideFaceIdx == 0) { // left
|
|
||||||
K1 *= multxMinus[interiorElemCartIdx];
|
|
||||||
K2 *= multx[exteriorElemCartIdx];
|
|
||||||
}
|
|
||||||
|
|
||||||
else if (insideFaceIdx == 3) { // back
|
|
||||||
K1 *= multy[interiorElemCartIdx];
|
|
||||||
K2 *= multyMinus[exteriorElemCartIdx];
|
|
||||||
}
|
|
||||||
else if (insideFaceIdx == 2) { // front
|
|
||||||
K1 *= multyMinus[interiorElemCartIdx];
|
|
||||||
K2 *= multy[exteriorElemCartIdx];
|
|
||||||
}
|
|
||||||
|
|
||||||
else if (insideFaceIdx == 5) { // top
|
|
||||||
K1 *= multz[interiorElemCartIdx];
|
|
||||||
K2 *= multzMinus[exteriorElemCartIdx];
|
|
||||||
}
|
|
||||||
else if (insideFaceIdx == 4) { // bottom
|
|
||||||
K1 *= multzMinus[interiorElemCartIdx];
|
|
||||||
K2 *= multz[exteriorElemCartIdx];
|
|
||||||
}
|
|
||||||
|
|
||||||
// element-wise harmonic average
|
|
||||||
auto &K = intersectionIntrinsicPermeability_[intersectIdx];
|
|
||||||
K = 0.0;
|
|
||||||
for (int i = 0; i < dimWorld; ++i)
|
|
||||||
for (int j = 0; j < dimWorld; ++j)
|
|
||||||
K[i][j] = Opm::utils::harmonicAverage(K1[i][j], K2[i][j]);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
std::vector<Scalar> porosity_;
|
|
||||||
std::vector<DimMatrix> intrinsicPermeability_;
|
|
||||||
|
|
||||||
// the intrinsic permeabilities for interior faces. since grids may be
|
|
||||||
// non-conforming, and there does not seem to be a mapper for interfaces in DUNE,
|
|
||||||
// these transmissibilities are accessed via the (elementIndex1, elementIndex2) pairs
|
|
||||||
// of the interfaces where
|
|
||||||
//
|
|
||||||
// elementIndex1 = min(interiorElementIndex, exteriorElementIndex)
|
|
||||||
//
|
|
||||||
// and
|
|
||||||
//
|
|
||||||
// elementIndex2 = max(interiorElementIndex, exteriorElementIndex)
|
|
||||||
//
|
|
||||||
// To make this perform better, this is first mingled into a single index using
|
|
||||||
//
|
|
||||||
// intersectionIndex = elementIndex1*numElements + elementIndex2
|
|
||||||
//
|
|
||||||
// as the index for the hash map.
|
|
||||||
std::unordered_map<size_t, DimMatrix> intersectionIntrinsicPermeability_;
|
|
||||||
|
|
||||||
std::vector<unsigned short> materialParamTableIdx_;
|
|
||||||
std::vector<MaterialLawParams> materialParams_;
|
|
||||||
|
|
||||||
std::vector<BlackOilFluidState> initialFluidStates_;
|
|
||||||
|
|
||||||
Scalar temperature_;
|
|
||||||
|
|
||||||
EclWellManager<TypeTag> wellManager_;
|
|
||||||
EclipseSummaryWriter summaryWriter_;
|
|
||||||
};
|
|
||||||
} // namespace Ewoms
|
|
||||||
|
|
||||||
#endif
|
|
||||||
@ -1,40 +0,0 @@
|
|||||||
/*
|
|
||||||
Copyright (C) 2014 by Andreas Lauser
|
|
||||||
|
|
||||||
This file is part of the Open Porous Media project (OPM).
|
|
||||||
|
|
||||||
OPM 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.
|
|
||||||
|
|
||||||
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
/*!
|
|
||||||
* \file
|
|
||||||
*
|
|
||||||
* \brief A general-purpose simulator for Eclipse decks using the
|
|
||||||
* black-oil model.
|
|
||||||
*/
|
|
||||||
#include "config.h"
|
|
||||||
|
|
||||||
#include <ewoms/common/start.hh>
|
|
||||||
#include <ewoms/models/blackoil/blackoilmodel.hh>
|
|
||||||
#include "problems/eclproblem.hh"
|
|
||||||
|
|
||||||
namespace Opm {
|
|
||||||
namespace Properties {
|
|
||||||
NEW_TYPE_TAG(EclProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem));
|
|
||||||
}}
|
|
||||||
|
|
||||||
int main(int argc, char **argv)
|
|
||||||
{
|
|
||||||
typedef TTAG(EclProblem) ProblemTypeTag;
|
|
||||||
return Ewoms::start<ProblemTypeTag>(argc, argv);
|
|
||||||
}
|
|
||||||
Loading…
Reference in New Issue
Block a user