23#ifndef OPM_FLOW_GENERIC_PROBLEM_IMPL_HPP
24#define OPM_FLOW_GENERIC_PROBLEM_IMPL_HPP
26#ifndef OPM_FLOW_GENERIC_PROBLEM_HPP
31#include <dune/common/parametertree.hh>
33#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
34#include <opm/input/eclipse/EclipseState/Tables/OverburdTable.hpp>
35#include <opm/input/eclipse/EclipseState/Tables/RockwnodTable.hpp>
36#include <opm/input/eclipse/Schedule/Schedule.hpp>
37#include <opm/input/eclipse/Units/Units.hpp>
48#include <boost/date_time.hpp>
50#include <fmt/format.h>
51#include <fmt/ranges.h>
58template<
class Gr
idView,
class Flu
idSystem>
61 const Schedule& schedule,
62 const GridView& gridView)
66 , lookUpData_(gridView)
85 ? Parameters::Get<Parameters::NumPressurePointsEquil>()
86 : eclState.getTableManager().getEqldims().getNumDepthNodesP();
91template<
class Gr
idView,
class Flu
idSystem>
95 const Schedule& schedule,
96 const GridView& gridView)
111template<
class Gr
idView,
class Flu
idSystem>
122 "Usage: "+std::string(argv[0]) +
" [OPTIONS] [ECL_DECK_FILENAME]\n"
126template<
class Gr
idView,
class Flu
idSystem>
131 return briefDescription_;
134template<
class Gr
idView,
class Flu
idSystem>
137 std::function<std::array<int,3>(
const unsigned)> ijkIndex)
139 const auto& rock_config = eclState_.getSimulationConfig().rock_config();
143 const auto& comp = rock_config.comp();
145 std::transform(comp.begin(), comp.end(), std::back_inserter(rockParams_),
148 return RockParams{static_cast<Scalar>(c.pref),
149 static_cast<Scalar>(c.compressibility)};
154 if (rock_config.store()) {
155 OpmLog::warning(
"ROCKOPTS item 2 set to STORE, ROCK item 1 replaced with initial (equilibrated) pressures");
159 readRockCompactionParameters_();
161 unsigned numElem = gridView_.size(0);
162 if (eclState_.fieldProps().has_int(rock_config.rocknum_property())) {
164 std::function<void(
int,
int)> valueCheck = [&ijkIndex,&rock_config,
this](
int fieldPropValue,
int coarseElemIdx)
166 auto fmtError = [fieldPropValue, coarseElemIdx,&ijkIndex,&rock_config](
const char* type, std::size_t size)
168 return fmt::format(
"{} table index {} for elem {} read from {}"
169 " is out of bounds for number of tables {}",
170 type, fieldPropValue,
171 ijkIndex(coarseElemIdx),
172 rock_config.rocknum_property(), size);
174 if (!rockCompPoroMult_.empty() &&
175 fieldPropValue >
static_cast<int>(rockCompPoroMult_.size())) {
176 throw std::runtime_error(fmtError(
"Rock compaction",
177 rockCompPoroMult_.size()));
179 if (!rockCompPoroMultWc_.empty() &&
180 fieldPropValue >
static_cast<int>(rockCompPoroMultWc_.size())) {
181 throw std::runtime_error(fmtError(
"Rock water compaction",
182 rockCompPoroMultWc_.size()));
186 rockTableIdx_ = this->lookUpData_.template assignFieldPropsIntOnLeaf<short unsigned int>(eclState_.fieldProps(),
187 rock_config.rocknum_property(),
193 const auto& overburdTables = eclState_.getTableManager().getOverburdTables();
194 if (!overburdTables.empty() && !rock_config.store()) {
195 overburdenPressure_.resize(numElem,0.0);
196 std::size_t numRocktabTables = rock_config.num_rock_tables();
198 if (overburdTables.size() != numRocktabTables)
199 throw std::runtime_error(
std::to_string(numRocktabTables) +
" OVERBURD tables is expected, but " +
std::to_string(overburdTables.size()) +
" is provided");
201 std::vector<Tabulated1DFunction<Scalar>> overburdenTables(numRocktabTables);
202 for (std::size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
203 const OverburdTable& overburdTable = overburdTables.template getTable<OverburdTable>(regionIdx);
204 overburdenTables[regionIdx].setXYContainers(overburdTable.getDepthColumn(),overburdTable.getOverburdenPressureColumn());
207 for (std::size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) {
208 unsigned tableIdx = 0;
209 if (!rockTableIdx_.empty()) {
210 tableIdx = rockTableIdx_[elemIdx];
212 overburdenPressure_[elemIdx] =
213 overburdenTables[tableIdx].eval(cellCenterDepths[elemIdx],
true);
216 else if (!overburdTables.empty() && rock_config.store()) {
217 OpmLog::warning(
"ROCKOPTS item 2 set to STORE, OVERBURD ignored!");
221template<
class Gr
idView,
class Flu
idSystem>
225 const auto& rock_config = eclState_.getSimulationConfig().rock_config();
227 if (!rock_config.active())
230 unsigned numElem = gridView_.size(0);
231 switch (rock_config.hysteresis_mode()) {
232 case RockConfig::Hysteresis::REVERS:
234 case RockConfig::Hysteresis::IRREVERS:
237 minRefPressure_.resize(numElem, 1e99);
240 throw std::runtime_error(
"Not support ROCKOMP hysteresis option ");
243 std::size_t numRocktabTables = rock_config.num_rock_tables();
244 bool waterCompaction = rock_config.water_compaction();
246 if (!waterCompaction) {
247 const auto& rocktabTables = eclState_.getTableManager().getRocktabTables();
248 if (rocktabTables.size() != numRocktabTables)
249 throw std::runtime_error(
"ROCKCOMP is activated." +
std::to_string(numRocktabTables)
250 +
" ROCKTAB tables is expected, but " +
std::to_string(rocktabTables.size()) +
" is provided");
252 rockCompPoroMult_.resize(numRocktabTables);
253 rockCompTransMult_.resize(numRocktabTables);
254 for (std::size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
255 const auto& rocktabTable = rocktabTables.template getTable<RocktabTable>(regionIdx);
256 const auto& pressureColumn = rocktabTable.getPressureColumn();
257 const auto& poroColumn = rocktabTable.getPoreVolumeMultiplierColumn();
258 const auto& transColumn = rocktabTable.getTransmissibilityMultiplierColumn();
259 rockCompPoroMult_[regionIdx].setXYContainers(pressureColumn, poroColumn);
260 rockCompTransMult_[regionIdx].setXYContainers(pressureColumn, transColumn);
263 const auto& rock2dTables = eclState_.getTableManager().getRock2dTables();
264 const auto& rock2dtrTables = eclState_.getTableManager().getRock2dtrTables();
265 const auto& rockwnodTables = eclState_.getTableManager().getRockwnodTables();
266 maxWaterSaturation_.resize(numElem, 0.0);
268 if (rock2dTables.size() != numRocktabTables)
269 throw std::runtime_error(
"Water compation option is selected in ROCKCOMP." +
std::to_string(numRocktabTables)
270 +
" ROCK2D tables is expected, but " +
std::to_string(rock2dTables.size()) +
" is provided");
272 if (rockwnodTables.size() != numRocktabTables)
273 throw std::runtime_error(
"Water compation option is selected in ROCKCOMP." +
std::to_string(numRocktabTables)
274 +
" ROCKWNOD tables is expected, but " +
std::to_string(rockwnodTables.size()) +
" is provided");
276 rockCompPoroMultWc_.resize(numRocktabTables,
TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical));
277 for (std::size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
278 const RockwnodTable& rockwnodTable = rockwnodTables.template getTable<RockwnodTable>(regionIdx);
279 const auto& rock2dTable = rock2dTables[regionIdx];
281 if (rockwnodTable.getSaturationColumn().size() != rock2dTable.sizeMultValues())
282 throw std::runtime_error(
"Number of entries in ROCKWNOD and ROCK2D needs to match.");
284 for (std::size_t xIdx = 0; xIdx < rock2dTable.size(); ++xIdx) {
285 rockCompPoroMultWc_[regionIdx].appendXPos(rock2dTable.getPressureValue(xIdx));
286 for (std::size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx)
287 rockCompPoroMultWc_[regionIdx].appendSamplePoint(xIdx,
288 rockwnodTable.getSaturationColumn()[yIdx],
289 rock2dTable.getPvmultValue(xIdx, yIdx));
293 if (!rock2dtrTables.empty()) {
294 rockCompTransMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical));
295 for (std::size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
296 const RockwnodTable& rockwnodTable = rockwnodTables.template getTable<RockwnodTable>(regionIdx);
297 const auto& rock2dtrTable = rock2dtrTables[regionIdx];
299 if (rockwnodTable.getSaturationColumn().size() != rock2dtrTable.sizeMultValues())
300 throw std::runtime_error(
"Number of entries in ROCKWNOD and ROCK2DTR needs to match.");
302 for (std::size_t xIdx = 0; xIdx < rock2dtrTable.size(); ++xIdx) {
303 rockCompTransMultWc_[regionIdx].appendXPos(rock2dtrTable.getPressureValue(xIdx));
304 for (std::size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx)
305 rockCompTransMultWc_[regionIdx].appendSamplePoint(xIdx,
306 rockwnodTable.getSaturationColumn()[yIdx],
307 rock2dtrTable.getTransMultValue(xIdx, yIdx));
314template<
class Gr
idView,
class Flu
idSystem>
319 if (this->rockParams_.empty())
322 unsigned tableIdx = 0;
323 if (!this->rockTableIdx_.empty()) {
324 tableIdx = this->rockTableIdx_[globalSpaceIdx];
326 return this->rockParams_[tableIdx].compressibility;
329template<
class Gr
idView,
class Flu
idSystem>
332porosity(
unsigned globalSpaceIdx,
unsigned timeIdx)
const
334 return this->referencePorosity_[timeIdx][globalSpaceIdx];
337template<
class Gr
idView,
class Flu
idSystem>
340rockFraction(
unsigned elementIdx,
unsigned timeIdx)
const
346 auto porosity = this->lookUpData_.fieldPropDouble(eclState_.fieldProps(),
"PORO", elementIdx);
347 return referencePorosity(elementIdx, timeIdx) / porosity * (1 - porosity);
350template<
class Gr
idView,
class Flu
idSystem>
353updateNum(
const std::string& name, std::vector<T>& numbers, std::size_t num_regions)
355 if (!eclState_.fieldProps().has_int(name))
358 std::function<void(T,
int)> valueCheck = [num_regions,name](T fieldPropValue, [[maybe_unused]]
int fieldPropIdx) {
359 if ( fieldPropValue > (
int)num_regions) {
360 throw std::runtime_error(
"Values larger than maximum number of regions "
363 if ( fieldPropValue <= 0) {
364 throw std::runtime_error(
"zero or negative values provided for region array: " + name);
368 numbers = this->lookUpData_.template assignFieldPropsIntOnLeaf<T>(eclState_.fieldProps(), name,
372template<
class Gr
idView,
class Flu
idSystem>
376 const auto num_regions = eclState_.getTableManager().getTabdims().getNumPVTTables();
377 updateNum(
"PVTNUM", pvtnum_, num_regions);
380template<
class Gr
idView,
class Flu
idSystem>
384 const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables();
385 updateNum(
"SATNUM", satnum_, num_regions);
388template<
class Gr
idView,
class Flu
idSystem>
392 const auto num_regions = 1;
393 updateNum(
"MISCNUM", miscnum_, num_regions);
396template<
class Gr
idView,
class Flu
idSystem>
400 const auto num_regions = 1;
401 updateNum(
"PLMIXNUM", plmixnum_, num_regions);
404template<
class Gr
idView,
class Flu
idSystem>
408 const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
409 return (oilVaporizationControl.getType() == OilVaporizationProperties::OilVaporization::VAPPARS);
412template<
class Gr
idView,
class Flu
idSystem>
417 if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) {
420 std::ostringstream ss;
421 boost::posix_time::time_facet* facet =
new boost::posix_time::time_facet(
"%d-%b-%Y");
422 boost::posix_time::ptime curDateTime =
423 boost::posix_time::from_time_t(schedule_.simTime(episodeIdx));
424 ss.imbue(std::locale(std::locale::classic(), facet));
425 ss <<
"Report step " << episodeIdx + 1
426 <<
"/" << schedule_.size() - 1
427 <<
" at day " << schedule_.seconds(episodeIdx)/(24*3600)
428 <<
"/" << schedule_.seconds(schedule_.size() - 1)/(24*3600)
429 <<
", date = " << curDateTime.date()
431 OpmLog::info(ss.str());
434 const auto& events = schedule_[episodeIdx].events();
437 if (episodeIdx > 0 && enableTuning_ && events.hasEvent(ScheduleEvents::TUNING_CHANGE))
439 const auto& sched_state = schedule_[episodeIdx];
440 const auto& tuning = sched_state.tuning();
441 initialTimeStepSize_ = sched_state.max_next_tstep(enableTuning_);
442 maxTimeStepAfterWellEvent_ = tuning.TMAXWC;
449template<
class Gr
idView,
class Flu
idSystem>
459 if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) {
460 std::ostringstream ss;
461 boost::posix_time::time_facet* facet =
new boost::posix_time::time_facet(
"%d-%b-%Y");
462 boost::posix_time::ptime date = boost::posix_time::from_time_t(startTime) + boost::posix_time::milliseconds(
static_cast<long long>(time / prefix::milli));
463 ss.imbue(std::locale(std::locale::classic(), facet));
464 ss <<
"\nTime step " << timeStepIndex <<
", stepsize "
465 << unit::convert::to(timeStepSize, unit::day) <<
" days,"
466 <<
" at day " << (double)unit::convert::to(time, unit::day)
467 <<
"/" << (double)unit::convert::to(endTime, unit::day)
468 <<
", date = " << date;
469 OpmLog::info(ss.str());
473template<
class Gr
idView,
class Flu
idSystem>
477 FluidSystem::initFromState(eclState_, schedule_);
480template<
class Gr
idView,
class Flu
idSystem>
485 bool enablePolymerMolarWeight,
488 auto getArray = [](
const std::vector<double>& input)
490 if constexpr (std::is_same_v<Scalar,double>) {
493 return std::vector<Scalar>{input.begin(), input.end()};
498 if (eclState_.fieldProps().has_double(
"SSOL")) {
499 solventSaturation_ = getArray(eclState_.fieldProps().get_double(
"SSOL"));
501 solventSaturation_.resize(numDof, 0.0);
504 solventRsw_.resize(numDof, 0.0);
508 if (eclState_.fieldProps().has_double(
"SPOLY")) {
509 polymer_.concentration = getArray(eclState_.fieldProps().get_double(
"SPOLY"));
511 polymer_.concentration.resize(numDof, 0.0);
515 if (enablePolymerMolarWeight) {
516 if (eclState_.fieldProps().has_double(
"SPOLYMW")) {
517 polymer_.moleWeight = getArray(eclState_.fieldProps().get_double(
"SPOLYMW"));
519 polymer_.moleWeight.resize(numDof, 0.0);
524 if (eclState_.fieldProps().has_double(
"SMICR")) {
525 micp_.microbialConcentration = getArray(eclState_.fieldProps().get_double(
"SMICR"));
527 micp_.microbialConcentration.resize(numDof, 0.0);
529 if (eclState_.fieldProps().has_double(
"SOXYG")) {
530 micp_.oxygenConcentration = getArray(eclState_.fieldProps().get_double(
"SOXYG"));
532 micp_.oxygenConcentration.resize(numDof, 0.0);
534 if (eclState_.fieldProps().has_double(
"SUREA")) {
535 micp_.ureaConcentration = getArray(eclState_.fieldProps().get_double(
"SUREA"));
537 micp_.ureaConcentration.resize(numDof, 0.0);
539 if (eclState_.fieldProps().has_double(
"SBIOF")) {
540 micp_.biofilmConcentration = getArray(eclState_.fieldProps().get_double(
"SBIOF"));
542 micp_.biofilmConcentration.resize(numDof, 0.0);
544 if (eclState_.fieldProps().has_double(
"SCALC")) {
545 micp_.calciteConcentration = getArray(eclState_.fieldProps().get_double(
"SCALC"));
547 micp_.calciteConcentration.resize(numDof, 0.0);
552template<
class Gr
idView,
class Flu
idSystem>
557 if (maxWaterSaturation_.empty())
560 return maxWaterSaturation_[globalDofIdx];
563template<
class Gr
idView,
class Flu
idSystem>
568 if (minRefPressure_.empty())
571 return minRefPressure_[globalDofIdx];
574template<
class Gr
idView,
class Flu
idSystem>
579 if (overburdenPressure_.empty())
582 return overburdenPressure_[elementIdx];
585template<
class Gr
idView,
class Flu
idSystem>
590 if (solventSaturation_.empty())
593 return solventSaturation_[elemIdx];
596template<
class Gr
idView,
class Flu
idSystem>
601 if (solventRsw_.empty())
604 return solventRsw_[elemIdx];
609template<
class Gr
idView,
class Flu
idSystem>
614 if (polymer_.concentration.empty()) {
618 return polymer_.concentration[elemIdx];
621template<
class Gr
idView,
class Flu
idSystem>
626 if (polymer_.moleWeight.empty()) {
630 return polymer_.moleWeight[elemIdx];
633template<
class Gr
idView,
class Flu
idSystem>
638 if (micp_.microbialConcentration.empty()) {
645template<
class Gr
idView,
class Flu
idSystem>
650 if (micp_.oxygenConcentration.empty()) {
657template<
class Gr
idView,
class Flu
idSystem>
662 if (micp_.ureaConcentration.empty()) {
669template<
class Gr
idView,
class Flu
idSystem>
674 if (micp_.biofilmConcentration.empty()) {
681template<
class Gr
idView,
class Flu
idSystem>
686 if (micp_.calciteConcentration.empty()) {
693template<
class Gr
idView,
class Flu
idSystem>
700 return pvtnum_[elemIdx];
703template<
class Gr
idView,
class Flu
idSystem>
710 return satnum_[elemIdx];
713template<
class Gr
idView,
class Flu
idSystem>
717 if (miscnum_.empty())
720 return miscnum_[elemIdx];
723template<
class Gr
idView,
class Flu
idSystem>
727 if (plmixnum_.empty())
730 return plmixnum_[elemIdx];
733template<
class Gr
idView,
class Flu
idSystem>
738 if (polymer_.maxAdsorption.empty()) {
742 return polymer_.maxAdsorption[elemIdx];
745template<
class Gr
idView,
class Flu
idSystem>
755 this->micp_ == rhs.
micp_;
Defines some fundamental parameters for all models.
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition: FlowGenericProblem.hpp:61
UniformXTabulated2DFunction< Scalar > TabulatedTwoDFunction
Definition: FlowGenericProblem.hpp:64
Scalar oxygenConcentration(unsigned elemIdx) const
Returns the initial oxygen concentration for a given a cell index.
Definition: FlowGenericProblem_impl.hpp:648
Scalar microbialConcentration(unsigned elemIdx) const
Returns the initial microbial concentration for a given a cell index.
Definition: FlowGenericProblem_impl.hpp:636
PolymerSolutionContainer< Scalar > polymer_
Definition: FlowGenericProblem.hpp:334
static std::string briefDescription()
Returns a human readable description of the problem for the help message.
Definition: FlowGenericProblem_impl.hpp:129
Scalar initialTimeStepSize_
Definition: FlowGenericProblem.hpp:345
std::vector< Scalar > solventSaturation_
Definition: FlowGenericProblem.hpp:339
bool enableDriftCompensation_
Definition: FlowGenericProblem.hpp:351
static FlowGenericProblem serializationTestObject(const EclipseState &eclState, const Schedule &schedule, const GridView &gridView)
Definition: FlowGenericProblem_impl.hpp:94
FlowGenericProblem(const EclipseState &eclState, const Schedule &schedule, const GridView &gridView)
Definition: FlowGenericProblem_impl.hpp:60
bool enableTuning_
Definition: FlowGenericProblem.hpp:344
Scalar biofilmConcentration(unsigned elemIdx) const
Returns the initial biofilm concentration for a given a cell index.
Definition: FlowGenericProblem_impl.hpp:672
MICPSolutionContainer< Scalar > micp_
Definition: FlowGenericProblem.hpp:341
void readRockParameters_(const std::vector< Scalar > &cellCenterDepths, std::function< std::array< int, 3 >(const unsigned)> ijkIndex)
Definition: FlowGenericProblem_impl.hpp:136
std::vector< Scalar > maxOilSaturation_
Definition: FlowGenericProblem.hpp:335
int numPressurePointsEquil_
Definition: FlowGenericProblem.hpp:349
std::vector< Scalar > maxWaterSaturation_
Definition: FlowGenericProblem.hpp:336
void initFluidSystem_()
Definition: FlowGenericProblem_impl.hpp:475
Scalar maxTimeStepAfterWellEvent_
Definition: FlowGenericProblem.hpp:346
Scalar calciteConcentration(unsigned elemIdx) const
Returns the initial calcite concentration for a given a cell index.
Definition: FlowGenericProblem_impl.hpp:684
Scalar ureaConcentration(unsigned elemIdx) const
Returns the initial urea concentration for a given a cell index.
Definition: FlowGenericProblem_impl.hpp:660
std::vector< Scalar > minRefPressure_
Definition: FlowGenericProblem.hpp:337
std::vector< Scalar > solventRsw_
Definition: FlowGenericProblem.hpp:340
static std::string helpPreamble(int, const char **argv)
Returns the string that is printed before the list of command line parameters in the help message.
Definition: FlowGenericProblem_impl.hpp:114
std::vector< Scalar > overburdenPressure_
Definition: FlowGenericProblem.hpp:338
bool explicitRockCompaction_
Definition: FlowGenericProblem.hpp:352
typename FluidSystem::Scalar Scalar
Definition: FlowGenericProblem.hpp:63
Declare the properties used by the infrastructure code of the finite volume discretizations.
auto Get(bool errorIfNotRegistered=true)
Retrieve a runtime parameter.
Definition: parametersystem.hpp:187
Definition: blackoilboundaryratevector.hh:39
bool operator==(const aligned_allocator< T1, Alignment > &, const aligned_allocator< T2, Alignment > &) noexcept
Definition: alignedallocator.hh:200
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
This file provides the infrastructure to retrieve run-time parameters.
static MICPSolutionContainer serializationTestObject()
Definition: EclTimeSteppingParams.hpp:45
static PolymerSolutionContainer serializationTestObject()