From a8d5c72248636ca16cb7b02278cc73923c71c734 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 28 Nov 2014 12:58:48 +0100 Subject: [PATCH] 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! --- examples/problems/eclproblem.hh | 984 -------------------------------- tests/ecl_blackoil.cc | 40 -- 2 files changed, 1024 deletions(-) delete mode 100644 examples/problems/eclproblem.hh delete mode 100644 tests/ecl_blackoil.cc diff --git a/examples/problems/eclproblem.hh b/examples/problems/eclproblem.hh deleted file mode 100644 index 7cc11d8a0..000000000 --- a/examples/problems/eclproblem.hh +++ /dev/null @@ -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 . -*/ -/*! - * \file - * - * \copydoc Ewoms::EclProblem - */ -#ifndef EWOMS_ECL_PROBLEM_HH -#define EWOMS_ECL_PROBLEM_HH - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -#include - -// for this simulator to make sense, dune-cornerpoint and opm-parser -// must be available -#include -#include -#include - -#include -#include -#include - -#include - -#include -#include - -namespace Ewoms { -template -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); - -// 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 OilWaterTraits; - - typedef Opm::TwoPhaseMaterialTraits GasOilTraits; - - typedef Opm::ThreePhaseMaterialTraits Traits; - - //typedef typename Opm::PiecewiseLinearTwoPhaseMaterial OilWaterLaw; - //typedef typename Opm::PiecewiseLinearTwoPhaseMaterial GasOilLaw; - - typedef typename Opm::SplineTwoPhaseMaterial OilWaterLaw; - typedef typename Opm::SplineTwoPhaseMaterial GasOilLaw; - -public: - typedef Opm::EclDefaultMaterial 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 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 EclipseSummaryWriter; - - typedef Dune::FieldMatrix 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 - const DimMatrix &intrinsicPermeability(const Context &context, - int spaceIdx, - int timeIdx) const - { - int globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx); - return intrinsicPermeability_[globalSpaceIdx]; - } - - /*! - * \copydoc FvBaseMultiPhaseProblem::intersectionIntrinsicPermeability - */ - template - 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 - Scalar porosity(const Context &context, int spaceIdx, int timeIdx) const - { - int globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx); - return porosity_[globalSpaceIdx]; - } - - /*! - * \copydoc FvBaseMultiPhaseProblem::materialLawParams - */ - template - 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 - 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 - 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 - 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 - 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 - 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 &permxData = - eclipseState->getDoubleGridProperty("PERMX")->getData(); - std::vector permyData(permxData); - if (eclipseState->hasDoubleGridProperty("PERMY")) - permyData = eclipseState->getDoubleGridProperty("PERMY")->getData(); - std::vector 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 &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 &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 &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 &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 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 &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(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 &waterSaturationData = - deck->getKeyword("SWAT")->getSIDoubleData(); - const std::vector &gasSaturationData = - deck->getKeyword("SGAS")->getSIDoubleData(); - const std::vector &pressureData = - deck->getKeyword("PRESSURE")->getSIDoubleData(); - const std::vector &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 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 multx(numElements, 1.0); - std::vector multy(numElements, 1.0); - std::vector multz(numElements, 1.0); - std::vector multxMinus(numElements, 1.0); - std::vector multyMinus(numElements, 1.0); - std::vector 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(); - const auto& elemEndIt = this->gridView().template end(); - 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 porosity_; - std::vector 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 intersectionIntrinsicPermeability_; - - std::vector materialParamTableIdx_; - std::vector materialParams_; - - std::vector initialFluidStates_; - - Scalar temperature_; - - EclWellManager wellManager_; - EclipseSummaryWriter summaryWriter_; -}; -} // namespace Ewoms - -#endif diff --git a/tests/ecl_blackoil.cc b/tests/ecl_blackoil.cc deleted file mode 100644 index 2f2ba1990..000000000 --- a/tests/ecl_blackoil.cc +++ /dev/null @@ -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 . -*/ -/*! - * \file - * - * \brief A general-purpose simulator for Eclipse decks using the - * black-oil model. - */ -#include "config.h" - -#include -#include -#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(argc, argv); -}