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 2014-11-28 12:58:48 +01:00
parent 08c97f478f
commit 47eafa47f4
5 changed files with 373 additions and 419 deletions

View File

@ -0,0 +1,264 @@
/*
Copyright (C) 2011-2013 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::EclOutputBlackOilModule
*/
#ifndef EWOMS_ECL_OUTPUT_BLACK_OIL_MODULE_HH
#define EWOMS_ECL_OUTPUT_BLACK_OIL_MODULE_HH
#include "eclwriter.hh"
#include <ewoms/io/baseoutputmodule.hh>
#include <opm/core/utility/PropertySystem.hpp>
#include <ewoms/common/parametersystem.hh>
#include <dune/common/fvector.hh>
#include <type_traits>
namespace Opm {
namespace Properties {
// create new type tag for the VTK multi-phase output
NEW_TYPE_TAG(EclOutputBlackOil);
// create the property tags needed for the multi phase module
NEW_PROP_TAG(EclOutputWriteSaturations);
NEW_PROP_TAG(EclOutputWritePressures);
NEW_PROP_TAG(EclOutputWriteGasDissolutionFactor);
NEW_PROP_TAG(EclOutputWriteGasFormationVolumeFactor);
NEW_PROP_TAG(EclOutputWriteOilFormationVolumeFactor);
NEW_PROP_TAG(EclOutputWriteOilSaturationPressure);
// set default values for what quantities to output
SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteSaturations, true);
SET_BOOL_PROP(EclOutputBlackOil, EclOutputWritePressures, true);
SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteGasDissolutionFactor, true);
SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteGasFormationVolumeFactor, true);
SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteOilFormationVolumeFactor, true);
SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteOilSaturationPressure, true);
}} // namespace Opm, Properties
namespace Ewoms {
// forward declaration
template <class TypeTag>
class EcfvDiscretization;
/*!
* \ingroup EclOutput
*
* \brief Output module for the results black oil model writing in
* ECL binary format.
*/
template <class TypeTag>
class EclOutputBlackOilModule : public BaseOutputModule<TypeTag>
{
typedef BaseOutputModule<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef Ewoms::EclWriter<TypeTag> EclWriter;
enum { numPhases = FluidSystem::numPhases };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
enum { gasCompIdx = FluidSystem::gasCompIdx };
typedef typename ParentType::ScalarBuffer ScalarBuffer;
public:
EclOutputBlackOilModule(const Simulator &simulator)
: ParentType(simulator)
{ }
/*!
* \brief Register all run-time parameters for the multi-phase VTK output
* module.
*/
static void registerParameters()
{
EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteSaturations,
"Include the saturations of all fluid phases in the "
"ECL output files");
EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWritePressures,
"Include the absolute pressures of all fluid phases in the "
"ECL output files");
EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteGasDissolutionFactor,
"Include the gas dissolution factor in the "
"ECL output files");
EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteGasFormationVolumeFactor,
"Include the gas formation volume factor in the "
"ECL output files");
EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteOilFormationVolumeFactor,
"Include the oil formation volume factor of saturated oil "
"in the ECL output files");
EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteOilSaturationPressure,
"Include the saturation pressure of oil in the "
"ECL output files");
}
/*!
* \brief Allocate memory for the scalar fields we would like to
* write to the VTK file.
*/
void allocBuffers()
{
if (!std::is_same<Discretization, Ewoms::EcfvDiscretization<TypeTag> >::value)
return;
auto bufferType = ParentType::ElementBuffer;
if (saturationsOutput_()) {
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx)
this->resizeScalarBuffer_(saturation_[phaseIdx], bufferType);
}
if (pressuresOutput_()) {
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx)
this->resizeScalarBuffer_(pressure_[phaseIdx], bufferType);
}
if (gasDissolutionFactorOutput_())
this->resizeScalarBuffer_(gasDissolutionFactor_, bufferType);
if (gasFormationVolumeFactorOutput_())
this->resizeScalarBuffer_(gasFormationVolumeFactor_, bufferType);
if (saturatedOilFormationVolumeFactorOutput_())
this->resizeScalarBuffer_(saturatedOilFormationVolumeFactor_, bufferType);
if (oilSaturationPressureOutput_())
this->resizeScalarBuffer_(oilSaturationPressure_, bufferType);
}
/*!
* \brief Modify the internal buffers according to the intensive quanties relevant
* for an element
*/
void processElement(const ElementContext &elemCtx)
{
if (!std::is_same<Discretization, Ewoms::EcfvDiscretization<TypeTag> >::value)
return;
for (int dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
const auto &fs = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState();
int globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
int regionIdx = elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex();
Scalar po = fs.pressure(oilPhaseIdx);
Scalar XoG = fs.massFraction(oilPhaseIdx, gasCompIdx);
if (saturationsOutput_()) {
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
saturation_[phaseIdx][globalDofIdx] = fs.saturation(phaseIdx);
Valgrind::CheckDefined(saturation_[phaseIdx][globalDofIdx]);
}
}
if (pressuresOutput_()) {
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
pressure_[phaseIdx][globalDofIdx] = fs.pressure(phaseIdx) / 1e5;
Valgrind::CheckDefined(pressure_[phaseIdx][globalDofIdx]);
}
}
if (gasDissolutionFactorOutput_()) {
gasDissolutionFactor_[globalDofIdx] =
FluidSystem::gasDissolutionFactor(po, regionIdx);
Valgrind::CheckDefined(gasDissolutionFactor_[globalDofIdx]);
}
if (gasFormationVolumeFactorOutput_()) {
gasFormationVolumeFactor_[globalDofIdx] =
FluidSystem::gasFormationVolumeFactor(po, regionIdx);
Valgrind::CheckDefined(gasFormationVolumeFactor_[globalDofIdx]);
}
if (saturatedOilFormationVolumeFactorOutput_()) {
saturatedOilFormationVolumeFactor_[globalDofIdx] =
FluidSystem::saturatedOilFormationVolumeFactor(po, regionIdx);
Valgrind::CheckDefined(saturatedOilFormationVolumeFactor_[globalDofIdx]);
}
if (oilSaturationPressureOutput_()) {
oilSaturationPressure_[globalDofIdx] =
FluidSystem::oilSaturationPressure(XoG, regionIdx);
Valgrind::CheckDefined(oilSaturationPressure_[globalDofIdx]);
}
}
}
/*!
* \brief Add all buffers to the VTK output writer.
*/
void commitBuffers(BaseOutputWriter &writer)
{
if (!std::is_same<Discretization, Ewoms::EcfvDiscretization<TypeTag> >::value)
return;
if (!dynamic_cast<EclWriter*>(&writer))
return; // this module only consideres ecl writers...
typename ParentType::BufferType bufferType = ParentType::ElementBuffer;
if (pressuresOutput_()) {
this->commitScalarBuffer_(writer, "PRESSURE", pressure_[oilPhaseIdx], bufferType);
this->commitScalarBuffer_(writer, "PGAS", pressure_[gasPhaseIdx], bufferType);
this->commitScalarBuffer_(writer, "PWAT", pressure_[waterPhaseIdx], bufferType);
}
if (saturationsOutput_()) {
this->commitScalarBuffer_(writer, "SWAT", saturation_[waterPhaseIdx], bufferType);
this->commitScalarBuffer_(writer, "SGAS", saturation_[gasPhaseIdx], bufferType);
// the oil saturation is _NOT_ written to disk. Instead, it is calculated by
// the visualization tool. Wondering why is probably a waste of time...
}
if (gasDissolutionFactorOutput_())
this->commitScalarBuffer_(writer, "RS", gasDissolutionFactor_, bufferType);
if (gasFormationVolumeFactorOutput_())
this->commitScalarBuffer_(writer, "BG", gasFormationVolumeFactor_, bufferType);
if (saturatedOilFormationVolumeFactorOutput_())
this->commitScalarBuffer_(writer, "BOSAT", saturatedOilFormationVolumeFactor_, bufferType);
if (oilSaturationPressureOutput_())
this->commitScalarBuffer_(writer, "PSAT", oilSaturationPressure_);
}
private:
static bool saturationsOutput_()
{ return EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteSaturations); }
static bool pressuresOutput_()
{ return EWOMS_GET_PARAM(TypeTag, bool, EclOutputWritePressures); }
static bool gasDissolutionFactorOutput_()
{ return EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteGasDissolutionFactor); }
static bool gasFormationVolumeFactorOutput_()
{ return EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteGasFormationVolumeFactor); }
static bool saturatedOilFormationVolumeFactorOutput_()
{ return EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteOilFormationVolumeFactor); }
static bool oilSaturationPressureOutput_()
{ return EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteOilSaturationPressure); }
ScalarBuffer saturation_[numPhases];
ScalarBuffer pressure_[numPhases];
ScalarBuffer gasDissolutionFactor_;
ScalarBuffer gasFormationVolumeFactor_;
ScalarBuffer saturatedOilFormationVolumeFactor_;
ScalarBuffer oilSaturationPressure_;
};
} // namespace Ewoms
#endif

View File

@ -24,11 +24,14 @@
#ifndef EWOMS_ECL_PROBLEM_HH
#define EWOMS_ECL_PROBLEM_HH
#include "eclgridmanager.hh"
#include "eclwellmanager.hh"
#include "eclwriter.hh"
#include "eclsummarywriter.hh"
#include "ecloutputblackoilmodule.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>
@ -60,7 +63,7 @@ class EclProblem;
namespace Opm {
namespace Properties {
NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclGridManager));
NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclGridManager, EclOutputBlackOil));
// The temperature inside the reservoir
NEW_PROP_TAG(Temperature);
@ -137,11 +140,11 @@ 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);
// ... but enable the ECL output by default
SET_BOOL_PROP(EclBaseProblem, EnableEclOutput, true);
// also enable the summary output.
SET_BOOL_PROP(EclBaseProblem, EnableEclipseSummaryOutput, true);
SET_BOOL_PROP(EclBaseProblem, EnableEclSummaryOutput, true);
// The default DGF file to load
SET_STRING_PROP(EclBaseProblem, GridFile, "data/ecl.DATA");
@ -151,8 +154,8 @@ namespace Ewoms {
/*!
* \ingroup TestProblems
*
* \brief This problem uses a deck in the format of the Eclipse
* simulator.
* \brief This problem simulates an input file given in the data format used by the
* commercial ECLiPSE simulator.
*/
template <class TypeTag>
class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
@ -185,7 +188,7 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
typedef Ewoms::EclipseSummaryWriter<TypeTag> EclipseSummaryWriter;
typedef Ewoms::EclSummaryWriter<TypeTag> EclSummaryWriter;
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
@ -197,11 +200,16 @@ public:
{
ParentType::registerParameters();
Ewoms::EclOutputBlackOilModule<TypeTag>::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");
EWOMS_REGISTER_PARAM(TypeTag, bool, EnableEclOutput,
"Write binary output which is compatible with the commercial "
"Eclipse simulator");
}
/*!
@ -210,8 +218,12 @@ public:
EclProblem(Simulator &simulator)
: ParentType(simulator)
, wellManager_(simulator)
, eclWriter_(simulator)
, summaryWriter_(simulator)
{ }
{
// add the output module for the Ecl binary output
simulator.model().addOutputModule(new Ewoms::EclOutputBlackOilModule<TypeTag>(simulator));
}
/*!
* \copydoc FvBaseProblem::finishInit
@ -234,9 +246,9 @@ public:
// 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());
wellManager_.init(simulator.gridManager().eclState());
// Start the first episode. For this, ask the Eclipse schedule.
// Start the first episode. For this, ask the ECL schedule.
Opm::TimeMapConstPtr timeMap = simulator.gridManager().schedule()->getTimeMap();
tm curTime = boost::posix_time::to_tm(timeMap->getStartTime(/*timeStepIdx=*/0));
@ -261,7 +273,7 @@ public:
* \brief Called by the simulator before an episode begins.
*/
void beginEpisode()
{ wellManager_.beginEpisode(this->simulator().gridManager().eclipseState()); }
{ wellManager_.beginEpisode(this->simulator().gridManager().eclState()); }
/*!
* \brief Called by the simulator before each time integration.
@ -302,8 +314,8 @@ public:
// ... then proceed to the next report step
Simulator &simulator = this->simulator();
Opm::EclipseStateConstPtr eclipseState = this->simulator().gridManager().eclipseState();
Opm::TimeMapConstPtr timeMap = eclipseState->getSchedule()->getTimeMap();
Opm::EclipseStateConstPtr eclState = this->simulator().gridManager().eclState();
Opm::TimeMapConstPtr timeMap = eclState->getSchedule()->getTimeMap();
// TimeMap deals with points in time, so the number of time
// intervals (i.e., report steps) is one less!
@ -339,6 +351,29 @@ public:
return this->simulator().episodeWillBeOver();
}
/*!
* \brief Write the requested quantities of the current solution into the output
* files.
*/
void writeOutput(bool verbose = true)
{
// calculate the time _after_ the time was updated
Scalar t = this->simulator().time() + this->simulator().timeStepSize();
// prepare the ECL and the VTK writers
if (enableEclOutput_())
eclWriter_.beginWrite(t);
// use the generic code to prepare the output fields and to
// write the desired VTK files.
ParentType::writeOutput(verbose);
if (enableEclOutput_()) {
this->model().appendOutputFields(eclWriter_);
eclWriter_.endWrite();
}
}
/*!
* \copydoc FvBaseMultiPhaseProblem::intrinsicPermeability
*/
@ -453,7 +488,7 @@ public:
/*!
* \copydoc FvBaseProblem::boundary
*
* Eclipse uses no-flow conditions for all boundaries. \todo really?
* ECLiPSE uses no-flow conditions for all boundaries. \todo really?
*/
template <class Context>
void boundary(BoundaryRateVector &values,
@ -506,10 +541,13 @@ public:
//! \}
private:
static bool enableEclOutput_()
{ return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); }
void readMaterialParameters_()
{
auto deck = this->simulator().gridManager().deck();
auto eclipseState = this->simulator().gridManager().eclipseState();
auto eclState = this->simulator().gridManager().eclState();
const auto &grid = this->simulator().gridManager().grid();
size_t numDof = this->model().numGridDof();
@ -521,18 +559,18 @@ private:
////////////////////////////////
// permeability
// read the intrinsic permeabilities from the eclipseState. Note that all arrays
// provided by eclipseState are one-per-cell of "uncompressed" grid, whereas the
// read the intrinsic permeabilities from the eclState. Note that all arrays
// provided by eclState are one-per-cell of "uncompressed" grid, whereas the
// dune-cornerpoint grid object might remove a few elements...
if (eclipseState->hasDoubleGridProperty("PERMX")) {
if (eclState->hasDoubleGridProperty("PERMX")) {
const std::vector<double> &permxData =
eclipseState->getDoubleGridProperty("PERMX")->getData();
eclState->getDoubleGridProperty("PERMX")->getData();
std::vector<double> permyData(permxData);
if (eclipseState->hasDoubleGridProperty("PERMY"))
permyData = eclipseState->getDoubleGridProperty("PERMY")->getData();
if (eclState->hasDoubleGridProperty("PERMY"))
permyData = eclState->getDoubleGridProperty("PERMY")->getData();
std::vector<double> permzData(permxData);
if (eclipseState->hasDoubleGridProperty("PERMZ"))
permzData = eclipseState->getDoubleGridProperty("PERMZ")->getData();
if (eclState->hasDoubleGridProperty("PERMZ"))
permzData = eclState->getDoubleGridProperty("PERMZ")->getData();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
int cartesianElemIdx = grid.globalCell()[dofIdx];
@ -546,13 +584,13 @@ private:
}
else
OPM_THROW(std::logic_error,
"Can't read the intrinsic permeability from the eclipse state. "
"Can't read the intrinsic permeability from the ecl state. "
"(The PERM{X,Y,Z} keywords are missing)");
// apply the NTG keyword to the X and Y permeabilities
if (eclipseState->hasDoubleGridProperty("NTG")) {
if (eclState->hasDoubleGridProperty("NTG")) {
const std::vector<double> &ntgData =
eclipseState->getDoubleGridProperty("NTG")->getData();
eclState->getDoubleGridProperty("NTG")->getData();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
int cartesianElemIdx = grid.globalCell()[dofIdx];
@ -567,9 +605,9 @@ private:
////////////////////////////////
// compute the porosity
if (eclipseState->hasDoubleGridProperty("PORO")) {
if (eclState->hasDoubleGridProperty("PORO")) {
const std::vector<double> &poroData =
eclipseState->getDoubleGridProperty("PORO")->getData();
eclState->getDoubleGridProperty("PORO")->getData();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
int cartesianElemIdx = grid.globalCell()[dofIdx];
@ -577,14 +615,14 @@ private:
}
}
else
OPM_THROW(std::logic_error,
"Can't read the porosity from the eclipse state. "
OPM_THROW(std::runtime_error,
"Can't read the porosity from the ECL state object. "
"(The PORO keyword is missing)");
// apply the NTG keyword to the porosity
if (eclipseState->hasDoubleGridProperty("NTG")) {
if (eclState->hasDoubleGridProperty("NTG")) {
const std::vector<double> &ntgData =
eclipseState->getDoubleGridProperty("NTG")->getData();
eclState->getDoubleGridProperty("NTG")->getData();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
int cartesianElemIdx = grid.globalCell()[dofIdx];
@ -593,9 +631,9 @@ private:
}
// apply the MULTPV keyword to the porosity
if (eclipseState->hasDoubleGridProperty("MULTPV")) {
if (eclState->hasDoubleGridProperty("MULTPV")) {
const std::vector<double> &multpvData =
eclipseState->getDoubleGridProperty("MULTPV")->getData();
eclState->getDoubleGridProperty("MULTPV")->getData();
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
int cartesianElemIdx = grid.globalCell()[dofIdx];
@ -605,8 +643,8 @@ private:
////////////////////////////////
// fluid parameters
const auto& swofTables = eclipseState->getSwofTables();
const auto& sgofTables = eclipseState->getSgofTables();
const auto& swofTables = eclState->getSwofTables();
const auto& sgofTables = eclState->getSgofTables();
// the number of tables for the SWOF and the SGOF keywords
// must be identical
@ -644,7 +682,7 @@ private:
owParams.finalize();
goParams.finalize();
// compute the connate water saturation. In Eclipse decks that is defined as
// compute the connate water saturation. In ECL decks that is defined as
// the first saturation value of the SWOF keyword.
Scalar Swco = SwColumn.front();
materialParams_[tableIdx].setConnateWaterSaturation(Swco);
@ -656,9 +694,9 @@ private:
}
// set the index of the table to be used
if (eclipseState->hasIntGridProperty("SATNUM")) {
if (eclState->hasIntGridProperty("SATNUM")) {
const std::vector<int> &satnumData =
eclipseState->getIntGridProperty("SATNUM")->getData();
eclState->getIntGridProperty("SATNUM")->getData();
materialParamTableIdx_.resize(numDof);
for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
@ -668,7 +706,7 @@ private:
assert(1 <= satnumData[dofIdx]);
assert(satnumData[dofIdx] <= static_cast<int>(numSatfuncTables));
// Eclipse uses Fortran-style indices which start at
// ECL uses Fortran-style indices which start at
// 1, but this here is C++...
materialParamTableIdx_[dofIdx] = satnumData[cartesianElemIdx] - 1;
}
@ -681,7 +719,7 @@ private:
void initFluidSystem_()
{
const auto deck = this->simulator().gridManager().deck();
const auto eclipseState = this->simulator().gridManager().eclipseState();
const auto eclState = this->simulator().gridManager().eclState();
FluidSystem::initBegin();
@ -697,9 +735,9 @@ private:
// so far, we require the presence of the PVTO, PVTW and PVDG
// keywords...
FluidSystem::setPvtoTable(eclipseState->getPvtoTables()[regionIdx], regionIdx);
FluidSystem::setPvtoTable(eclState->getPvtoTables()[regionIdx], regionIdx);
FluidSystem::setPvtw(deck->getKeyword("PVTW"), regionIdx);
FluidSystem::setPvdgTable(eclipseState->getPvdgTables()[regionIdx], regionIdx);
FluidSystem::setPvdgTable(eclState->getPvdgTables()[regionIdx], regionIdx);
}
FluidSystem::initEnd();
@ -713,11 +751,11 @@ private:
if (!deck->hasKeyword("SWAT") ||
!deck->hasKeyword("SGAS"))
OPM_THROW(std::runtime_error,
"So far, the Eclipse input file requires the presence of the SWAT "
"So far, the ECL 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 "
"So far, the ECL input file requires the presence of the PRESSURE "
"keyword");
if (!deck->hasKeyword("DISGAS"))
OPM_THROW(std::runtime_error,
@ -725,7 +763,7 @@ private:
" (DISGAS keyword is missing)");
if (!deck->hasKeyword("RS"))
OPM_THROW(std::runtime_error,
"The Eclipse input file requires the presence of the RS keyword");
"The ECL input file requires the presence of the RS keyword");
if (deck->hasKeyword("VAPOIL"))
OPM_THROW(std::runtime_error,
@ -733,7 +771,7 @@ private:
" (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");
"The ECL input file requires the RV keyword to be non-present");
size_t numDof = this->model().numGridDof();
@ -845,7 +883,7 @@ private:
void computeFaceIntrinsicPermeabilities_()
{
auto eclipseState = this->simulator().gridManager().eclipseState();
auto eclState = this->simulator().gridManager().eclState();
const auto &grid = this->simulator().gridManager().grid();
int numElements = this->gridView().size(/*codim=*/0);
@ -859,18 +897,18 @@ private:
// 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();
if (eclState->hasDoubleGridProperty("MULTX"))
multx = eclState->getDoubleGridProperty("MULTX")->getData();
if (eclState->hasDoubleGridProperty("MULTX-"))
multxMinus = eclState->getDoubleGridProperty("MULTX-")->getData();
if (eclState->hasDoubleGridProperty("MULTY"))
multy = eclState->getDoubleGridProperty("MULTY")->getData();
if (eclState->hasDoubleGridProperty("MULTY-"))
multyMinus = eclState->getDoubleGridProperty("MULTY-")->getData();
if (eclState->hasDoubleGridProperty("MULTZ"))
multz = eclState->getDoubleGridProperty("MULTZ")->getData();
if (eclState->hasDoubleGridProperty("MULTZ-"))
multzMinus = eclState->getDoubleGridProperty("MULTZ-")->getData();
// resize the hash map to a appropriate size for a conforming 3D grid
float maxLoadFactor = intersectionIntrinsicPermeability_.max_load_factor();
@ -977,7 +1015,9 @@ private:
Scalar temperature_;
EclWellManager<TypeTag> wellManager_;
EclipseSummaryWriter summaryWriter_;
EclWriter<TypeTag> eclWriter_;
EclSummaryWriter summaryWriter_;
};
} // namespace Ewoms

View File

@ -24,7 +24,8 @@
#ifndef EWOMS_ECL_WELL_MANAGER_HH
#define EWOMS_ECL_WELL_MANAGER_HH
#include <ewoms/wells/eclpeacemanwell.hh>
#include "eclpeacemanwell.hh"
#include <ewoms/disc/common/fvbaseproperties.hh>
#include <opm/parser/eclipse/Deck/Deck.hpp>
@ -469,9 +470,9 @@ protected:
auto& model = simulator_.model();
model.clearAuxiliaryModules();
auto eclState = simulator_.gridManager().eclipseState();
const auto &deckSchedule = eclState->getSchedule();
const Grid &grid = simulator_.gridManager().grid();
auto eclStatePtr = simulator_.gridManager().eclState();
const auto& deckSchedule = eclStatePtr->getSchedule();
const Grid& grid = simulator_.gridManager().grid();
const GridView gridView = simulator_.gridManager().gridView();
const std::vector<Opm::WellConstPtr>& deckWells = deckSchedule->getWells(reportStepIdx);
for (size_t deckWellIdx = 0; deckWellIdx < deckWells.size(); ++deckWellIdx) {

View File

@ -1,351 +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::EclipseWriter
*/
#ifndef EWOMS_ECLIPSE_WRITER_HH
#define EWOMS_ECLIPSE_WRITER_HH
#include "baseoutputwriter.hh"
#include "ertwrappers.hh"
#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
#include <opm/material/Valgrind.hpp>
#include <boost/algorithm/string.hpp>
#include <list>
#include <utility>
#include <string>
#include <limits>
#include <sstream>
#include <fstream>
#include <type_traits>
namespace Opm {
namespace Properties {
NEW_PROP_TAG(EnableEclipseOutput);
}}
namespace Ewoms {
template <class TypeTag>
class EclipseWriter;
template <class TypeTag>
class EclGridManager;
// we need to use specialization here because we need to call
// EclGridManager-specific methods.
// this is the stub class for non-cornerpoint grids
template <class TypeTag, class GridManagerType>
class EclipseWriterHelper
{
friend class EclipseWriter<TypeTag>;
static void writeHeaders_(EclipseWriter<TypeTag> &writer)
{
OPM_THROW(std::logic_error,
"Eclipse binary output requires to use Dune::EclGridManager (is: "
<< Dune::className<GridManagerType>());
}
};
// this is the "real" code for cornerpoint grids
template <class TypeTag>
class EclipseWriterHelper<TypeTag, Ewoms::EclGridManager<TypeTag> >
{
friend class EclipseWriter<TypeTag>;
static void writeHeaders_(EclipseWriter<TypeTag> &writer)
{
typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
if (!std::is_same<Discretization, Ewoms::EcfvDiscretization<TypeTag> >::value)
OPM_THROW(std::logic_error,
"Eclipse binary output only works for the element centered "
"finite volume discretization.");
#if ! HAVE_ERT
OPM_THROW(std::logic_error,
"Eclipse binary output requires the ERT libraries");
#else
// set the index of the first time step written to 0...
writer.reportStepIdx_ = 0;
char* egridRawFileName = ecl_util_alloc_filename(/*outputDir=*/"./",
writer.caseName().c_str(),
ECL_EGRID_FILE,
/*formatted=*/false, // -> write binary output
writer.reportStepIdx_);
std::string egridFileName(egridRawFileName);
std::free(egridRawFileName);
ErtGrid ertGrid(writer.simulator_.gridManager().eclipseGrid());
ertGrid.write(egridFileName, writer.reportStepIdx_);
#endif
}
};
/*!
* \brief Implements writing Eclipse binary output files.
*
* Caveats:
* - For this class to do do anything meaningful, you need to have the
* ERT libraries with development headers installed and the ERT
* build system test must pass sucessfully.
* - The only DUNE grid which is currently supported is Dune::CpGrid
* from the OPM module "dune-cornerpoint". Using another grid won't
* fail at compile time but you will provoke a fatal exception as
* soon as you try to write an Eclipse output file.
* - This class requires to use the black oil model with the element
* centered finite volume discretization.
* - MPI-parallel computations are not (yet?) supported.
*/
template <class TypeTag>
class EclipseWriter : public BaseOutputWriter
{
typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, GridManager) GridManager;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef BaseOutputWriter::ScalarBuffer ScalarBuffer;
typedef BaseOutputWriter::VectorBuffer VectorBuffer;
typedef BaseOutputWriter::TensorBuffer TensorBuffer;
friend class EclipseWriterHelper<TypeTag, GridManager>;
public:
EclipseWriter(const Simulator &simulator)
: simulator_(simulator)
, gridView_(simulator_.gridView())
, elementMapper_(gridView_)
, vertexMapper_(gridView_)
{
reportStepIdx_ = 0;
}
~EclipseWriter()
{ }
/*!
* \brief Returns the name of the simulation.
*
* This is the prefix of the files written to disk.
*/
std::string caseName() const
{ return boost::to_upper_copy(simulator_.problem().name()); }
/*!
* \brief Updates the internal data structures after mesh
* refinement.
*
* If the grid changes between two calls of beginWrite(), this
* method _must_ be called before the second beginWrite()!
*/
void gridChanged()
{
elementMapper_.update();
vertexMapper_.update();
}
/*!
* \brief Called whenever a new time step must be written.
*/
void beginWrite(double t)
{
if (enableEclipseOutput_() && reportStepIdx_ == 0)
EclipseWriterHelper<TypeTag, GridManager>::writeHeaders_(*this);
}
/*
* \brief Add a vertex-centered scalar field to the output.
*
* For the EclipseWriter, this method is a no-op which throws a
* std::logic_error exception
*/
void attachScalarVertexData(ScalarBuffer &buf, std::string name)
{
OPM_THROW(std::logic_error,
"The EclipseWriter can only write element based quantities!");
}
/*
* \brief Add a vertex-centered vector field to the output.
*
* For the EclipseWriter, this method is a no-op which throws a
* std::logic_error exception
*/
void attachVectorVertexData(VectorBuffer &buf, std::string name)
{
OPM_THROW(std::logic_error,
"The EclipseWriter can only write element based quantities!");
}
/*
* \brief Add a vertex-centered tensor field to the output.
*/
void attachTensorVertexData(TensorBuffer &buf, std::string name)
{
OPM_THROW(std::logic_error,
"The EclipseWriter can only write element based quantities!");
}
/*!
* \brief Add a scalar quantity to the output.
*
* The buffer must exist at least until the call to endWrite()
* finishes. Modifying the buffer between the call to this method
* and endWrite() results in _undefined behavior_.
*/
void attachScalarElementData(ScalarBuffer &buf, std::string name)
{
attachedBuffers_.push_back(std::pair<std::string, ScalarBuffer*>(name, &buf));
}
/*
* \brief Add a element-centered vector field to the output.
*
* For the EclipseWriter, this method is a no-op which throws a
* std::logic_error exception
*/
void attachVectorElementData(VectorBuffer &buf, std::string name)
{
OPM_THROW(std::logic_error,
"Currently, the EclipseWriter can only write scalar quantities!");
}
/*
* \brief Add a element-centered tensor field to the output.
*/
void attachTensorElementData(TensorBuffer &buf, std::string name)
{
OPM_THROW(std::logic_error,
"Currently, the EclipseWriter can only write scalar quantities!");
}
/*!
* \brief Finalizes the current writer.
*
* This means that everything will be written to disk, except if
* the onlyDiscard argument is true. In this case, no output is
* written but the 'janitorial' jobs at the end of a time step are
* still done.
*/
void endWrite(bool onlyDiscard = false)
{
if (onlyDiscard || !enableEclipseOutput_()) {
// detach all buffers
attachedBuffers_.resize(0);
return;
}
#if !HAVE_ERT
OPM_THROW(std::runtime_error,
"The ERT libraries must be available to write Eclipse output!");
#else
ErtRestartFile restartFile(simulator_, reportStepIdx_);
restartFile.writeHeader(simulator_, reportStepIdx_);
ErtSolution solution(restartFile);
auto bufIt = attachedBuffers_.begin();
const auto &bufEndIt = attachedBuffers_.end();
for (; bufIt != bufEndIt; ++ bufIt) {
const std::string &name = bufIt->first;
const ScalarBuffer &buffer = *bufIt->second;
std::shared_ptr<const ErtKeyword<float>>
bufKeyword(new ErtKeyword<float>(name, buffer));
solution.add(bufKeyword);
}
// detach all buffers
attachedBuffers_.resize(0);
// next time we take the next report step
++ reportStepIdx_;
#endif
}
/*!
* \brief Write the multi-writer's state to a restart file.
*/
template <class Restarter>
void serialize(Restarter &res)
{
res.serializeSectionBegin("EclipseWriter");
res.serializeStream() << reportStepIdx_ << "\n";
res.serializeSectionEnd();
}
/*!
* \brief Read the multi-writer's state from a restart file.
*/
template <class Restarter>
void deserialize(Restarter &res)
{
res.deserializeSectionBegin("EclipseWriter");
res.deserializeStream() >> reportStepIdx_;
res.deserializeSectionEnd();
}
private:
static bool enableEclipseOutput_()
{ return EWOMS_GET_PARAM(TypeTag, bool, EnableEclipseOutput); }
// make sure the field is well defined if running under valgrind
// and make sure that all values can be displayed by paraview
void sanitizeBuffer_(std::vector<float> &b)
{
static bool warningPrinted = false;
for (size_t i = 0; i < b.size(); ++i) {
Valgrind::CheckDefined(b[i]);
if (!warningPrinted && !std::isfinite(b[i])) {
std::cerr << "WARNING: data field written to disk contains non-finite entries!\n";
warningPrinted = true;
}
// set values which are too small to 0 to avoid possible
// problems
if (std::abs(b[i]) < std::numeric_limits<float>::min()) {
b[i] = 0.0;
}
}
}
const Simulator& simulator_;
const GridView gridView_;
ElementMapper elementMapper_;
VertexMapper vertexMapper_;
double curTime_;
int reportStepIdx_;
std::list<std::pair<std::string, ScalarBuffer*> > attachedBuffers_;
};
} // namespace Ewoms
#endif