OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4765 from hnil/change_poly_alugrid

Browse Source 
removed use of hidden private defines for poly and alugrid
This commit is contained in:
Arne Morten Kvarving
2023-08-09 11:15:02 +02:00
committed by GitHub
parent 60b71dac25 e55008b336
commit c6f1aa0110
24 changed files with 4163 additions and 4307 deletions
Download Patch File Download Diff File Expand all files Collapse all files

17
CMakeLists.txt
View File

@@ -559,12 +559,6 @@ foreach(OBJ ${FLOW_VARIANT_MODELS})
set_property(TARGET flow_lib${OBJ} PROPERTY EXCLUDE_FROM_ALL ${FLOW_VARIANTS_DEFAULT_ENABLE_IF})
endforeach()
add_library(flow_libblackoil_poly OBJECT flow/flow_ebos_blackoil_poly.cpp)
if(TARGET fmt::fmt)
target_link_libraries(flow_libblackoil_poly fmt::fmt)
endif()
set_property(TARGET flow_libblackoil_poly PROPERTY EXCLUDE_FROM_ALL ${FLOW_POLY_ONLY_DEFAULT_ENABLE_IF})
opm_add_test(flow
ONLY_COMPILE
ALWAYS_ENABLE
@@ -577,17 +571,15 @@ opm_add_test(flow
$<TARGET_OBJECTS:moduleVersion>
)
opm_add_test(flow_poly
opm_add_test(flow_blackoil_polyhedralgrid
ONLY_COMPILE
ALWAYS_ENABLE
DEFAULT_ENABLE_IF ${FLOW_POLY_ONLY_DEFAULT_ENABLE_IF}
DEPENDS opmsimulators
LIBRARIES opmsimulators
SOURCES
flow/flow_blackoil_poly.cpp
$<TARGET_OBJECTS:flow_libblackoil_poly>
flow/flow_blackoil_polyhedralgrid.cpp
$<TARGET_OBJECTS:moduleVersion>)
target_compile_definitions(flow_poly PRIVATE USE_POLYHEDRALGRID)
opm_add_test(flow_distribute_z
ONLY_COMPILE
@@ -608,7 +600,7 @@ if(dune-alugrid_FOUND)
set(FLOW_ALUGRID_ONLY_DEFAULT_ENABLE_IF "TRUE")
endif()
opm_add_test(flow_alugrid
opm_add_test(flow_blackoil_alugrid
ONLY_COMPILE
ALWAYS_ENABLE
DEFAULT_ENABLE_IF ${FLOW_ALUGRID_ONLY_DEFAULT_ENABLE_IF}
@@ -617,7 +609,6 @@ if(dune-alugrid_FOUND)
SOURCES
flow/flow_blackoil_alugrid.cpp
$<TARGET_OBJECTS:moduleVersion>)
target_compile_definitions(flow_alugrid PRIVATE USE_ALUGRID)
endif()
if (BUILD_FLOW)
@@ -685,5 +676,3 @@ endif()
install(DIRECTORY doc/man1 DESTINATION ${CMAKE_INSTALL_MANDIR}
FILES_MATCHING PATTERN "*.1")

1144
ebos/collecttoiorank.cc
View File

File diff suppressed because it is too large Load Diff

1128
ebos/collecttoiorank_impl.hh Normal file
View File

File diff suppressed because it is too large Load Diff

847
ebos/eclgenericproblem.cc
View File

@@ -20,30 +20,17 @@
module for the precise wording of the license and the list of
copyright holders.
*/
#include <config.h>
#include <ebos/eclgenericproblem.hh>
#include <ebos/eclgenericproblem_impl.hh>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Tables/OverburdTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/RockwnodTable.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <dune/grid/common/defaultgridview.hh>
#include <dune/grid/common/gridview.hh>
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/polyhedralgrid.hh>
#include <dune/common/parametertree.hh>
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/3d/gridview.hh>
#include "alucartesianindexmapper.hh"
#endif // HAVE_DUNE_ALUGRID
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#if HAVE_DUNE_FEM
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
@@ -51,788 +38,8 @@
#include <ebos/femcpgridcompat.hh>
#endif // HAVE_DUNE_FEM
#include <boost/date_time.hpp>
#include <algorithm>
#include <limits>
#include <stdexcept>
#include <iostream>
namespace Opm {
int eclPositionalParameter(Dune::ParameterTree& tree,
std::set<std::string>& seenParams,
std::string& errorMsg,
const char** argv,
int paramIdx)
{
std::string param = argv[paramIdx];
size_t i = param.find('=');
if (i != std::string::npos) {
std::string oldParamName = param.substr(0, i);
std::string oldParamValue = param.substr(i+1);
std::string newParamName = "--" + oldParamName;
std::replace(newParamName.begin(),
newParamName.end(), '_' , '-');
errorMsg =
"The old syntax to specify parameters on the command line is no longer supported: "
"Try replacing '" + oldParamName + "=" + oldParamValue + "' with "+
"'" + newParamName + "=" + oldParamValue + "'!";
return 0;
}
if (seenParams.count("EclDeckFileName") > 0) {
errorMsg =
"Parameter 'EclDeckFileName' specified multiple times"
" as a command line parameter";
return 0;
}
tree["EclDeckFileName"] = argv[paramIdx];
seenParams.insert("EclDeckFileName");
return 1;
}
template<class GridView, class FluidSystem, class Scalar>
EclGenericProblem<GridView,FluidSystem,Scalar>::
EclGenericProblem(const EclipseState& eclState,
const Schedule& schedule,
const GridView& gridView)
: eclState_(eclState)
, schedule_(schedule)
, gridView_(gridView)
{
}
template<class GridView, class FluidSystem, class Scalar>
EclGenericProblem<GridView,FluidSystem,Scalar>
EclGenericProblem<GridView,FluidSystem,Scalar>::
serializationTestObject(const EclipseState& eclState,
const Schedule& schedule,
const GridView& gridView)
{
EclGenericProblem result(eclState, schedule, gridView);
result.maxOilSaturation_ = {1.0, 2.0};
result.maxWaterSaturation_ = {6.0};
result.minOilPressure_ = {7.0, 8.0, 9.0, 10.0};
result.overburdenPressure_ = {11.0};
result.solventSaturation_ = {15.0};
result.polymer_ = PolymerSolutionContainer<Scalar>::serializationTestObject();
result.micp_ = MICPSolutionContainer<Scalar>::serializationTestObject();
result.lastRv_ = {21.0};
result.maxDRv_ = {22.0, 23.0};
result.convectiveDrs_ = {24.0, 25.0, 26.0};
result.lastRs_ = {27.0};
result.maxDRs_ = {28.0};
result.dRsDtOnlyFreeGas_ = {false, true};
return result;
}
template<class GridView, class FluidSystem, class Scalar>
std::string
EclGenericProblem<GridView,FluidSystem,Scalar>::
helpPreamble(int,
const char **argv)
{
std::string desc = EclGenericProblem::briefDescription();
if (!desc.empty())
desc = desc + "\n";
return
"Usage: "+std::string(argv[0]) + " [OPTIONS] [ECL_DECK_FILENAME]\n"
+ desc;
}
template<class GridView, class FluidSystem, class Scalar>
std::string
EclGenericProblem<GridView,FluidSystem,Scalar>::
briefDescription()
{
if (briefDescription_.empty())
return
"The Ecl-deck Black-Oil reservoir Simulator (ebos); a hydrocarbon "
"reservoir simulation program that processes ECL-formatted input "
"files that is part of the Open Porous Media project "
"(https://opm-project.org).\n"
"\n"
"THE GOAL OF THE `ebos` SIMULATOR IS TO CATER FOR THE NEEDS OF "
"DEVELOPMENT AND RESEARCH. No guarantees are made for production use!";
else
return briefDescription_;
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
readRockParameters_(const std::vector<Scalar>& cellCenterDepths)
{
const auto& rock_config = eclState_.getSimulationConfig().rock_config();
// read the rock compressibility parameters
{
const auto& comp = rock_config.comp();
rockParams_.clear();
for (const auto& c : comp)
rockParams_.push_back( { c.pref, c.compressibility } );
}
// read the parameters for water-induced rock compaction
readRockCompactionParameters_();
unsigned numElem = gridView_.size(0);
if (eclState_.fieldProps().has_int(rock_config.rocknum_property())) {
rockTableIdx_.resize(numElem);
const auto& num = eclState_.fieldProps().get_int(rock_config.rocknum_property());
for (size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) {
rockTableIdx_[elemIdx] = num[elemIdx] - 1;
}
}
// Store overburden pressure pr element
const auto& overburdTables = eclState_.getTableManager().getOverburdTables();
if (!overburdTables.empty()) {
overburdenPressure_.resize(numElem,0.0);
size_t numRocktabTables = rock_config.num_rock_tables();
if (overburdTables.size() != numRocktabTables)
throw std::runtime_error(std::to_string(numRocktabTables) +" OVERBURD tables is expected, but " + std::to_string(overburdTables.size()) +" is provided");
std::vector<Tabulated1DFunction<Scalar>> overburdenTables(numRocktabTables);
for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
const OverburdTable& overburdTable = overburdTables.template getTable<OverburdTable>(regionIdx);
overburdenTables[regionIdx].setXYContainers(overburdTable.getDepthColumn(),overburdTable.getOverburdenPressureColumn());
}
for (size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) {
unsigned tableIdx = 0;
if (!rockTableIdx_.empty()) {
tableIdx = rockTableIdx_[elemIdx];
}
overburdenPressure_[elemIdx] = overburdenTables[tableIdx].eval(cellCenterDepths[elemIdx], /*extrapolation=*/true);
}
}
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
readRockCompactionParameters_()
{
const auto& rock_config = eclState_.getSimulationConfig().rock_config();
if (!rock_config.active())
return; // deck does not enable rock compaction
unsigned numElem = gridView_.size(0);
switch (rock_config.hysteresis_mode()) {
case RockConfig::Hysteresis::REVERS:
break;
case RockConfig::Hysteresis::IRREVERS:
// interpolate the porv volume multiplier using the minimum pressure in the cell
// i.e. don't allow re-inflation.
minOilPressure_.resize(numElem, 1e99);
break;
default:
throw std::runtime_error("Not support ROCKOMP hysteresis option ");
}
size_t numRocktabTables = rock_config.num_rock_tables();
bool waterCompaction = rock_config.water_compaction();
if (!waterCompaction) {
const auto& rocktabTables = eclState_.getTableManager().getRocktabTables();
if (rocktabTables.size() != numRocktabTables)
throw std::runtime_error("ROCKCOMP is activated." + std::to_string(numRocktabTables)
+" ROCKTAB tables is expected, but " + std::to_string(rocktabTables.size()) +" is provided");
rockCompPoroMult_.resize(numRocktabTables);
rockCompTransMult_.resize(numRocktabTables);
for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
const auto& rocktabTable = rocktabTables.template getTable<RocktabTable>(regionIdx);
const auto& pressureColumn = rocktabTable.getPressureColumn();
const auto& poroColumn = rocktabTable.getPoreVolumeMultiplierColumn();
const auto& transColumn = rocktabTable.getTransmissibilityMultiplierColumn();
rockCompPoroMult_[regionIdx].setXYContainers(pressureColumn, poroColumn);
rockCompTransMult_[regionIdx].setXYContainers(pressureColumn, transColumn);
}
} else {
const auto& rock2dTables = eclState_.getTableManager().getRock2dTables();
const auto& rock2dtrTables = eclState_.getTableManager().getRock2dtrTables();
const auto& rockwnodTables = eclState_.getTableManager().getRockwnodTables();
maxWaterSaturation_.resize(numElem, 0.0);
if (rock2dTables.size() != numRocktabTables)
throw std::runtime_error("Water compation option is selected in ROCKCOMP." + std::to_string(numRocktabTables)
+" ROCK2D tables is expected, but " + std::to_string(rock2dTables.size()) +" is provided");
if (rockwnodTables.size() != numRocktabTables)
throw std::runtime_error("Water compation option is selected in ROCKCOMP." + std::to_string(numRocktabTables)
+" ROCKWNOD tables is expected, but " + std::to_string(rockwnodTables.size()) +" is provided");
//TODO check size match
rockCompPoroMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical));
for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
const RockwnodTable& rockwnodTable = rockwnodTables.template getTable<RockwnodTable>(regionIdx);
const auto& rock2dTable = rock2dTables[regionIdx];
if (rockwnodTable.getSaturationColumn().size() != rock2dTable.sizeMultValues())
throw std::runtime_error("Number of entries in ROCKWNOD and ROCK2D needs to match.");
for (size_t xIdx = 0; xIdx < rock2dTable.size(); ++xIdx) {
rockCompPoroMultWc_[regionIdx].appendXPos(rock2dTable.getPressureValue(xIdx));
for (size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx)
rockCompPoroMultWc_[regionIdx].appendSamplePoint(xIdx,
rockwnodTable.getSaturationColumn()[yIdx],
rock2dTable.getPvmultValue(xIdx, yIdx));
}
}
if (!rock2dtrTables.empty()) {
rockCompTransMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical));
for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
const RockwnodTable& rockwnodTable = rockwnodTables.template getTable<RockwnodTable>(regionIdx);
const auto& rock2dtrTable = rock2dtrTables[regionIdx];
if (rockwnodTable.getSaturationColumn().size() != rock2dtrTable.sizeMultValues())
throw std::runtime_error("Number of entries in ROCKWNOD and ROCK2DTR needs to match.");
for (size_t xIdx = 0; xIdx < rock2dtrTable.size(); ++xIdx) {
rockCompTransMultWc_[regionIdx].appendXPos(rock2dtrTable.getPressureValue(xIdx));
for (size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx)
rockCompTransMultWc_[regionIdx].appendSamplePoint(xIdx,
rockwnodTable.getSaturationColumn()[yIdx],
rock2dtrTable.getTransMultValue(xIdx, yIdx));
}
}
}
}
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
rockCompressibility(unsigned globalSpaceIdx) const
{
if (this->rockParams_.empty())
return 0.0;
unsigned tableIdx = 0;
if (!this->rockTableIdx_.empty()) {
tableIdx = this->rockTableIdx_[globalSpaceIdx];
}
return this->rockParams_[tableIdx].compressibility;
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
rockReferencePressure(unsigned globalSpaceIdx) const
{
if (this->rockParams_.empty())
return 1e5;
unsigned tableIdx = 0;
if (!this->rockTableIdx_.empty()) {
tableIdx = this->rockTableIdx_[globalSpaceIdx];
}
return this->rockParams_[tableIdx].referencePressure;
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
porosity(unsigned globalSpaceIdx, unsigned timeIdx) const
{
return this->referencePorosity_[timeIdx][globalSpaceIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
rockFraction(unsigned elementIdx, unsigned timeIdx) const
{
const auto& fp = eclState_.fieldProps();
const std::vector<double>& poro = fp.get_double("PORO");
// the reference porosity is defined as the accumulated pore volume divided by the
// geometric volume of the element. Note that it can
// be larger than 1.0 if porevolume multipliers are used
// to for instance implement larger boundary cells
Scalar porosity = poro[elementIdx];
return referencePorosity(elementIdx, timeIdx) / porosity * (1 - porosity);
}
template<class GridView, class FluidSystem, class Scalar>
template<class T>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
updateNum(const std::string& name, std::vector<T>& numbers, size_t num_regions)
{
if (!eclState_.fieldProps().has_int(name))
return;
const auto& numData = eclState_.fieldProps().get_int(name);
unsigned numElems = gridView_.size(/*codim=*/0);
numbers.resize(numElems);
for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) {
if (numData[elemIdx] > (int)num_regions) {
throw std::runtime_error("Values larger than maximum number of regions " + std::to_string(num_regions) + " provided in " + name);
} else if (numData[elemIdx] > 0) {
numbers[elemIdx] = static_cast<T>(numData[elemIdx]) - 1;
} else {
throw std::runtime_error("zero or negative values provided for region array: " + name);
}
}
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
updatePvtnum_()
{
const auto num_regions = eclState_.getTableManager().getTabdims().getNumPVTTables();
updateNum("PVTNUM", pvtnum_, num_regions);
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
updateSatnum_()
{
const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables();
updateNum("SATNUM", satnum_, num_regions);
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
updateMiscnum_()
{
const auto num_regions = 1; // we only support single region
updateNum("MISCNUM", miscnum_, num_regions);
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
updatePlmixnum_()
{
const auto num_regions = 1; // we only support single region
updateNum("PLMIXNUM", plmixnum_, num_regions);
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
updateKrnum_()
{
const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables();
updateNum("KRNUMX", krnumx_, num_regions);
updateNum("KRNUMY", krnumy_, num_regions);
updateNum("KRNUMZ", krnumz_, num_regions);
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
updateImbnum_()
{
const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables();
updateNum("IMBNUMX", imbnumx_, num_regions);
updateNum("IMBNUMY", imbnumy_, num_regions);
updateNum("IMBNUMZ", imbnumz_, num_regions);
}
template<class GridView, class FluidSystem, class Scalar>
bool EclGenericProblem<GridView,FluidSystem,Scalar>::
vapparsActive(int episodeIdx) const
{
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
return (oilVaporizationControl.getType() == OilVaporizationProperties::OilVaporization::VAPPARS);
}
template<class GridView, class FluidSystem, class Scalar>
bool EclGenericProblem<GridView,FluidSystem,Scalar>::
drsdtActive_(int episodeIdx) const
{
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
return (oilVaporizationControl.drsdtActive() && bothOilGasActive);
}
template<class GridView, class FluidSystem, class Scalar>
bool EclGenericProblem<GridView,FluidSystem,Scalar>::
drvdtActive_(int episodeIdx) const
{
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
return (oilVaporizationControl.drvdtActive() && bothOilGasActive);
}
template<class GridView, class FluidSystem, class Scalar>
bool EclGenericProblem<GridView,FluidSystem,Scalar>::
drsdtConvective_(int episodeIdx) const
{
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
return (oilVaporizationControl.drsdtConvective() && bothOilGasActive);
}
template<class GridView, class FluidSystem, class Scalar>
bool EclGenericProblem<GridView,FluidSystem,Scalar>::
beginEpisode_(bool enableExperiments,
int episodeIdx)
{
if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) {
// print some useful information in experimental mode. (the production
// simulator does this externally.)
std::ostringstream ss;
boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y");
boost::posix_time::ptime curDateTime =
boost::posix_time::from_time_t(schedule_.simTime(episodeIdx));
ss.imbue(std::locale(std::locale::classic(), facet));
ss << "Report step " << episodeIdx + 1
<< "/" << schedule_.size() - 1
<< " at day " << schedule_.seconds(episodeIdx)/(24*3600)
<< "/" << schedule_.seconds(schedule_.size() - 1)/(24*3600)
<< ", date = " << curDateTime.date()
<< "\n ";
OpmLog::info(ss.str());
}
const auto& events = schedule_[episodeIdx].events();
// react to TUNING changes
if (episodeIdx > 0 && enableTuning_ && events.hasEvent(ScheduleEvents::TUNING_CHANGE))
{
const auto& sched_state = schedule_[episodeIdx];
const auto& tuning = sched_state.tuning();
initialTimeStepSize_ = sched_state.max_next_tstep();
maxTimeStepAfterWellEvent_ = tuning.TMAXWC;
maxTimeStepSize_ = tuning.TSMAXZ;
restartShrinkFactor_ = 1./tuning.TSFCNV;
minTimeStepSize_ = tuning.TSMINZ;
return true;
}
return false;
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
beginTimeStep_(bool enableExperiments,
int episodeIdx,
int timeStepIndex,
Scalar startTime,
Scalar time,
Scalar timeStepSize,
Scalar endTime)
{
if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) {
std::ostringstream ss;
boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y");
boost::posix_time::ptime date = boost::posix_time::from_time_t(startTime) + boost::posix_time::milliseconds(static_cast<long long>(time / prefix::milli));
ss.imbue(std::locale(std::locale::classic(), facet));
ss <<"\nTime step " << timeStepIndex << ", stepsize "
<< unit::convert::to(timeStepSize, unit::day) << " days,"
<< " at day " << (double)unit::convert::to(time, unit::day)
<< "/" << (double)unit::convert::to(endTime, unit::day)
<< ", date = " << date;
OpmLog::info(ss.str());
}
// update explicit quantities between timesteps.
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
if (drsdtActive_(episodeIdx))
// DRSDT is enabled
for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx)
maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*timeStepSize;
if (drvdtActive_(episodeIdx))
// DRVDT is enabled
for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx)
maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*timeStepSize;
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
initFluidSystem_()
{
FluidSystem::initFromState(eclState_, schedule_);
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
readBlackoilExtentionsInitialConditions_(size_t numDof,
bool enableSolvent,
bool enablePolymer,
bool enablePolymerMolarWeight,
bool enableMICP)
{
if (enableSolvent) {
if (eclState_.fieldProps().has_double("SSOL"))
solventSaturation_ = eclState_.fieldProps().get_double("SSOL");
else
solventSaturation_.resize(numDof, 0.0);
}
if (enablePolymer) {
if (eclState_.fieldProps().has_double("SPOLY"))
polymer_.concentration = eclState_.fieldProps().get_double("SPOLY");
else
polymer_.concentration.resize(numDof, 0.0);
}
if (enablePolymerMolarWeight) {
if (eclState_.fieldProps().has_double("SPOLYMW"))
polymer_.moleWeight = eclState_.fieldProps().get_double("SPOLYMW");
else
polymer_.moleWeight.resize(numDof, 0.0);
}
if (enableMICP) {
if (eclState_.fieldProps().has_double("SMICR"))
micp_.microbialConcentration = eclState_.fieldProps().get_double("SMICR");
else
micp_.microbialConcentration.resize(numDof, 0.0);
if (eclState_.fieldProps().has_double("SOXYG"))
micp_.oxygenConcentration = eclState_.fieldProps().get_double("SOXYG");
else
micp_.oxygenConcentration.resize(numDof, 0.0);
if (eclState_.fieldProps().has_double("SUREA"))
micp_.ureaConcentration = eclState_.fieldProps().get_double("SUREA");
else
micp_.ureaConcentration.resize(numDof, 0.0);
if (eclState_.fieldProps().has_double("SBIOF"))
micp_.biofilmConcentration = eclState_.fieldProps().get_double("SBIOF");
else
micp_.biofilmConcentration.resize(numDof, 0.0);
if (eclState_.fieldProps().has_double("SCALC"))
micp_.calciteConcentration = eclState_.fieldProps().get_double("SCALC");
else
micp_.calciteConcentration.resize(numDof, 0.0);
}
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
maxWaterSaturation(unsigned globalDofIdx) const
{
if (maxWaterSaturation_.empty())
return 0.0;
return maxWaterSaturation_[globalDofIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
minOilPressure(unsigned globalDofIdx) const
{
if (minOilPressure_.empty())
return 0.0;
return minOilPressure_[globalDofIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
overburdenPressure(unsigned elementIdx) const
{
if (overburdenPressure_.empty())
return 0.0;
return overburdenPressure_[elementIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
solventSaturation(unsigned elemIdx) const
{
if (solventSaturation_.empty())
return 0;
return solventSaturation_[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
drsdtcon(unsigned elemIdx, int episodeIdx) const
{
if (convectiveDrs_.empty())
return 0;
// The episode index is set to -1 in the initialization phase.
// Output drsdt value for index 0
episodeIdx = std::max(episodeIdx, 0);
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
int pvtRegionIdx = pvtRegionIndex(elemIdx);
return oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*convectiveDrs_[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
polymerConcentration(unsigned elemIdx) const
{
if (polymer_.concentration.empty()) {
return 0;
}
return polymer_.concentration[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
polymerMolecularWeight(const unsigned elemIdx) const
{
if (polymer_.moleWeight.empty()) {
return 0.0;
}
return polymer_.moleWeight[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
microbialConcentration(unsigned elemIdx) const
{
if (micp_.microbialConcentration.empty()) {
return 0;
}
return micp_.microbialConcentration[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
oxygenConcentration(unsigned elemIdx) const
{
if (micp_.oxygenConcentration.empty()) {
return 0;
}
return micp_.oxygenConcentration[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
ureaConcentration(unsigned elemIdx) const
{
if (micp_.ureaConcentration.empty()) {
return 0;
}
return micp_.ureaConcentration[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
biofilmConcentration(unsigned elemIdx) const
{
if (micp_.biofilmConcentration.empty()) {
return 0;
}
return micp_.biofilmConcentration[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
calciteConcentration(unsigned elemIdx) const
{
if (micp_.calciteConcentration.empty()) {
return 0;
}
return micp_.calciteConcentration[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
unsigned EclGenericProblem<GridView,FluidSystem,Scalar>::
pvtRegionIndex(unsigned elemIdx) const
{
if (pvtnum_.empty())
return 0;
return pvtnum_[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
unsigned EclGenericProblem<GridView,FluidSystem,Scalar>::
satnumRegionIndex(unsigned elemIdx) const
{
if (satnum_.empty())
return 0;
return satnum_[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
unsigned EclGenericProblem<GridView,FluidSystem,Scalar>::
miscnumRegionIndex(unsigned elemIdx) const
{
if (miscnum_.empty())
return 0;
return miscnum_[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
unsigned EclGenericProblem<GridView,FluidSystem,Scalar>::
plmixnumRegionIndex(unsigned elemIdx) const
{
if (plmixnum_.empty())
return 0;
return plmixnum_[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
maxPolymerAdsorption(unsigned elemIdx) const
{
if (polymer_.maxAdsorption.empty()) {
return 0;
}
return polymer_.maxAdsorption[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
initDRSDT_(size_t numDof,
int episodeIdx)
{
// deal with DRSDT
unsigned ntpvt = eclState_.runspec().tabdims().getNumPVTTables();
//TODO We may want to only allocate these properties only if active.
//But since they may be activated at later time we need some more
//intrastructure to handle it
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
maxDRv_.resize(ntpvt, 1e30);
lastRv_.resize(numDof, 0.0);
maxDRs_.resize(ntpvt, 1e30);
dRsDtOnlyFreeGas_.resize(ntpvt, false);
lastRs_.resize(numDof, 0.0);
maxDRv_.resize(ntpvt, 1e30);
lastRv_.resize(numDof, 0.0);
maxOilSaturation_.resize(numDof, 0.0);
if (drsdtConvective_(episodeIdx)) {
convectiveDrs_.resize(numDof, 1.0);
}
}
}
template<class GridView, class FluidSystem, class Scalar>
bool EclGenericProblem<GridView,FluidSystem,Scalar>::
operator==(const EclGenericProblem& rhs) const
{
return this->maxOilSaturation_ == rhs.maxOilSaturation_ &&
this->maxWaterSaturation_ == rhs.maxWaterSaturation_ &&
this->minOilPressure_ == rhs.minOilPressure_ &&
this->overburdenPressure_ == rhs.overburdenPressure_ &&
this->polymer_ == rhs.polymer_ &&
this->solventSaturation_ == rhs.solventSaturation_ &&
this->micp_ == rhs.micp_ &&
this->lastRv_ == rhs.lastRv_ &&
this->maxDRv_ == rhs.maxDRv_ &&
this->convectiveDrs_ == rhs.convectiveDrs_ &&
this->lastRs_ == rhs.lastRs_ &&
this->maxDRs_ == rhs.maxDRs_ &&
this->dRsDtOnlyFreeGas_ == rhs.dRsDtOnlyFreeGas_;
}
#if HAVE_DUNE_FEM
template class EclGenericProblem<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
@@ -842,52 +49,12 @@ template class EclGenericProblem<Dune::Fem::GridPart2GridViewImpl<
Dune::CpGrid,
Dune::PartitionIteratorType(4),
false> >,
Opm::BlackOilFluidSystem<
double,
Opm::BlackOilDefaultIndexTraits>,
double>;
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericProblem<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>,
BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
double>;
template class EclGenericProblem<Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<
ALUGrid3CN,
Dune::PartitionIteratorType(4),
false>>,
Opm::BlackOilFluidSystem<
double,
Opm::BlackOilDefaultIndexTraits>,
double>;
#endif //HAVE_DUNE_ALUGRID
#else
template class EclGenericProblem<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
double>;
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = const Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = const Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
#endif
template class EclGenericProblem<Dune::GridView<Dune::ALU3dLeafGridViewTraits<ALUGrid3CN, Dune::PartitionIteratorType(4)>>,
BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
double>;
#endif //HAVE_DUNE_ALUGRID
#endif //HAVE_DUNE_FEM
template class EclGenericProblem<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>,
BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
double>;
} // namespace Opm
} // end namespace Opm

829
ebos/eclgenericproblem_impl.hh Normal file
View File

@@ -0,0 +1,829 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#ifndef EWOMS_GENERIC_ECL_PROBLEM_IMPL_HH
#define EWOMS_GENERIC_ECL_PROBLEM_IPML_HH
#include <ebos/eclgenericproblem.hh>
#include <dune/common/parametertree.hh>
#include <ebos/eclsolutioncontainers.hh>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Tables/OverburdTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/RockwnodTable.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <boost/date_time.hpp>
#include <algorithm>
#include <limits>
#include <stdexcept>
#include <iostream>
namespace Opm {
int eclPositionalParameter(Dune::ParameterTree& tree,
std::set<std::string>& seenParams,
std::string& errorMsg,
const char** argv,
int paramIdx)
{
std::string param = argv[paramIdx];
size_t i = param.find('=');
if (i != std::string::npos) {
std::string oldParamName = param.substr(0, i);
std::string oldParamValue = param.substr(i+1);
std::string newParamName = "--" + oldParamName;
std::replace(newParamName.begin(),
newParamName.end(), '_' , '-');
errorMsg =
"The old syntax to specify parameters on the command line is no longer supported: "
"Try replacing '" + oldParamName + "=" + oldParamValue + "' with "+
"'" + newParamName + "=" + oldParamValue + "'!";
return 0;
}
if (seenParams.count("EclDeckFileName") > 0) {
errorMsg =
"Parameter 'EclDeckFileName' specified multiple times"
" as a command line parameter";
return 0;
}
tree["EclDeckFileName"] = argv[paramIdx];
seenParams.insert("EclDeckFileName");
return 1;
}
template<class GridView, class FluidSystem, class Scalar>
EclGenericProblem<GridView,FluidSystem,Scalar>::
EclGenericProblem(const EclipseState& eclState,
const Schedule& schedule,
const GridView& gridView)
: eclState_(eclState)
, schedule_(schedule)
, gridView_(gridView)
{
}
template<class GridView, class FluidSystem, class Scalar>
EclGenericProblem<GridView,FluidSystem,Scalar>
EclGenericProblem<GridView,FluidSystem,Scalar>::
serializationTestObject(const EclipseState& eclState,
const Schedule& schedule,
const GridView& gridView)
{
EclGenericProblem result(eclState, schedule, gridView);
result.maxOilSaturation_ = {1.0, 2.0};
result.maxWaterSaturation_ = {6.0};
result.minOilPressure_ = {7.0, 8.0, 9.0, 10.0};
result.overburdenPressure_ = {11.0};
result.solventSaturation_ = {15.0};
result.polymer_ = PolymerSolutionContainer<Scalar>::serializationTestObject();
result.micp_ = MICPSolutionContainer<Scalar>::serializationTestObject();
result.lastRv_ = {21.0};
result.maxDRv_ = {22.0, 23.0};
result.convectiveDrs_ = {24.0, 25.0, 26.0};
result.lastRs_ = {27.0};
result.maxDRs_ = {28.0};
result.dRsDtOnlyFreeGas_ = {false, true};
return result;
}
template<class GridView, class FluidSystem, class Scalar>
std::string
EclGenericProblem<GridView,FluidSystem,Scalar>::
helpPreamble(int,
const char **argv)
{
std::string desc = EclGenericProblem::briefDescription();
if (!desc.empty())
desc = desc + "\n";
return
"Usage: "+std::string(argv[0]) + " [OPTIONS] [ECL_DECK_FILENAME]\n"
+ desc;
}
template<class GridView, class FluidSystem, class Scalar>
std::string
EclGenericProblem<GridView,FluidSystem,Scalar>::
briefDescription()
{
if (briefDescription_.empty())
return
"The Ecl-deck Black-Oil reservoir Simulator (ebos); a hydrocarbon "
"reservoir simulation program that processes ECL-formatted input "
"files that is part of the Open Porous Media project "
"(https://opm-project.org).\n"
"\n"
"THE GOAL OF THE `ebos` SIMULATOR IS TO CATER FOR THE NEEDS OF "
"DEVELOPMENT AND RESEARCH. No guarantees are made for production use!";
else
return briefDescription_;
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
readRockParameters_(const std::vector<Scalar>& cellCenterDepths)
{
const auto& rock_config = eclState_.getSimulationConfig().rock_config();
// read the rock compressibility parameters
{
const auto& comp = rock_config.comp();
rockParams_.clear();
for (const auto& c : comp)
rockParams_.push_back( { c.pref, c.compressibility } );
}
// read the parameters for water-induced rock compaction
readRockCompactionParameters_();
unsigned numElem = gridView_.size(0);
if (eclState_.fieldProps().has_int(rock_config.rocknum_property())) {
rockTableIdx_.resize(numElem);
const auto& num = eclState_.fieldProps().get_int(rock_config.rocknum_property());
for (size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) {
rockTableIdx_[elemIdx] = num[elemIdx] - 1;
}
}
// Store overburden pressure pr element
const auto& overburdTables = eclState_.getTableManager().getOverburdTables();
if (!overburdTables.empty()) {
overburdenPressure_.resize(numElem,0.0);
size_t numRocktabTables = rock_config.num_rock_tables();
if (overburdTables.size() != numRocktabTables)
throw std::runtime_error(std::to_string(numRocktabTables) +" OVERBURD tables is expected, but " + std::to_string(overburdTables.size()) +" is provided");
std::vector<Tabulated1DFunction<Scalar>> overburdenTables(numRocktabTables);
for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
const OverburdTable& overburdTable = overburdTables.template getTable<OverburdTable>(regionIdx);
overburdenTables[regionIdx].setXYContainers(overburdTable.getDepthColumn(),overburdTable.getOverburdenPressureColumn());
}
for (size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) {
unsigned tableIdx = 0;
if (!rockTableIdx_.empty()) {
tableIdx = rockTableIdx_[elemIdx];
}
overburdenPressure_[elemIdx] = overburdenTables[tableIdx].eval(cellCenterDepths[elemIdx], /*extrapolation=*/true);
}
}
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
readRockCompactionParameters_()
{
const auto& rock_config = eclState_.getSimulationConfig().rock_config();
if (!rock_config.active())
return; // deck does not enable rock compaction
unsigned numElem = gridView_.size(0);
switch (rock_config.hysteresis_mode()) {
case RockConfig::Hysteresis::REVERS:
break;
case RockConfig::Hysteresis::IRREVERS:
// interpolate the porv volume multiplier using the minimum pressure in the cell
// i.e. don't allow re-inflation.
minOilPressure_.resize(numElem, 1e99);
break;
default:
throw std::runtime_error("Not support ROCKOMP hysteresis option ");
}
size_t numRocktabTables = rock_config.num_rock_tables();
bool waterCompaction = rock_config.water_compaction();
if (!waterCompaction) {
const auto& rocktabTables = eclState_.getTableManager().getRocktabTables();
if (rocktabTables.size() != numRocktabTables)
throw std::runtime_error("ROCKCOMP is activated." + std::to_string(numRocktabTables)
+" ROCKTAB tables is expected, but " + std::to_string(rocktabTables.size()) +" is provided");
rockCompPoroMult_.resize(numRocktabTables);
rockCompTransMult_.resize(numRocktabTables);
for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
const auto& rocktabTable = rocktabTables.template getTable<RocktabTable>(regionIdx);
const auto& pressureColumn = rocktabTable.getPressureColumn();
const auto& poroColumn = rocktabTable.getPoreVolumeMultiplierColumn();
const auto& transColumn = rocktabTable.getTransmissibilityMultiplierColumn();
rockCompPoroMult_[regionIdx].setXYContainers(pressureColumn, poroColumn);
rockCompTransMult_[regionIdx].setXYContainers(pressureColumn, transColumn);
}
} else {
const auto& rock2dTables = eclState_.getTableManager().getRock2dTables();
const auto& rock2dtrTables = eclState_.getTableManager().getRock2dtrTables();
const auto& rockwnodTables = eclState_.getTableManager().getRockwnodTables();
maxWaterSaturation_.resize(numElem, 0.0);
if (rock2dTables.size() != numRocktabTables)
throw std::runtime_error("Water compation option is selected in ROCKCOMP." + std::to_string(numRocktabTables)
+" ROCK2D tables is expected, but " + std::to_string(rock2dTables.size()) +" is provided");
if (rockwnodTables.size() != numRocktabTables)
throw std::runtime_error("Water compation option is selected in ROCKCOMP." + std::to_string(numRocktabTables)
+" ROCKWNOD tables is expected, but " + std::to_string(rockwnodTables.size()) +" is provided");
//TODO check size match
rockCompPoroMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical));
for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
const RockwnodTable& rockwnodTable = rockwnodTables.template getTable<RockwnodTable>(regionIdx);
const auto& rock2dTable = rock2dTables[regionIdx];
if (rockwnodTable.getSaturationColumn().size() != rock2dTable.sizeMultValues())
throw std::runtime_error("Number of entries in ROCKWNOD and ROCK2D needs to match.");
for (size_t xIdx = 0; xIdx < rock2dTable.size(); ++xIdx) {
rockCompPoroMultWc_[regionIdx].appendXPos(rock2dTable.getPressureValue(xIdx));
for (size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx)
rockCompPoroMultWc_[regionIdx].appendSamplePoint(xIdx,
rockwnodTable.getSaturationColumn()[yIdx],
rock2dTable.getPvmultValue(xIdx, yIdx));
}
}
if (!rock2dtrTables.empty()) {
rockCompTransMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical));
for (size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
const RockwnodTable& rockwnodTable = rockwnodTables.template getTable<RockwnodTable>(regionIdx);
const auto& rock2dtrTable = rock2dtrTables[regionIdx];
if (rockwnodTable.getSaturationColumn().size() != rock2dtrTable.sizeMultValues())
throw std::runtime_error("Number of entries in ROCKWNOD and ROCK2DTR needs to match.");
for (size_t xIdx = 0; xIdx < rock2dtrTable.size(); ++xIdx) {
rockCompTransMultWc_[regionIdx].appendXPos(rock2dtrTable.getPressureValue(xIdx));
for (size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx)
rockCompTransMultWc_[regionIdx].appendSamplePoint(xIdx,
rockwnodTable.getSaturationColumn()[yIdx],
rock2dtrTable.getTransMultValue(xIdx, yIdx));
}
}
}
}
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
rockCompressibility(unsigned globalSpaceIdx) const
{
if (this->rockParams_.empty())
return 0.0;
unsigned tableIdx = 0;
if (!this->rockTableIdx_.empty()) {
tableIdx = this->rockTableIdx_[globalSpaceIdx];
}
return this->rockParams_[tableIdx].compressibility;
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
rockReferencePressure(unsigned globalSpaceIdx) const
{
if (this->rockParams_.empty())
return 1e5;
unsigned tableIdx = 0;
if (!this->rockTableIdx_.empty()) {
tableIdx = this->rockTableIdx_[globalSpaceIdx];
}
return this->rockParams_[tableIdx].referencePressure;
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
porosity(unsigned globalSpaceIdx, unsigned timeIdx) const
{
return this->referencePorosity_[timeIdx][globalSpaceIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
rockFraction(unsigned elementIdx, unsigned timeIdx) const
{
const auto& fp = eclState_.fieldProps();
const std::vector<double>& poro = fp.get_double("PORO");
// the reference porosity is defined as the accumulated pore volume divided by the
// geometric volume of the element. Note that it can
// be larger than 1.0 if porevolume multipliers are used
// to for instance implement larger boundary cells
Scalar porosity = poro[elementIdx];
return referencePorosity(elementIdx, timeIdx) / porosity * (1 - porosity);
}
template<class GridView, class FluidSystem, class Scalar>
template<class T>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
updateNum(const std::string& name, std::vector<T>& numbers, size_t num_regions)
{
if (!eclState_.fieldProps().has_int(name))
return;
const auto& numData = eclState_.fieldProps().get_int(name);
unsigned numElems = gridView_.size(/*codim=*/0);
numbers.resize(numElems);
for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) {
if (numData[elemIdx] > (int)num_regions) {
throw std::runtime_error("Values larger than maximum number of regions " + std::to_string(num_regions) + " provided in " + name);
} else if (numData[elemIdx] > 0) {
numbers[elemIdx] = static_cast<T>(numData[elemIdx]) - 1;
} else {
throw std::runtime_error("zero or negative values provided for region array: " + name);
}
}
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
updatePvtnum_()
{
const auto num_regions = eclState_.getTableManager().getTabdims().getNumPVTTables();
updateNum("PVTNUM", pvtnum_, num_regions);
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
updateSatnum_()
{
const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables();
updateNum("SATNUM", satnum_, num_regions);
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
updateMiscnum_()
{
const auto num_regions = 1; // we only support single region
updateNum("MISCNUM", miscnum_, num_regions);
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
updatePlmixnum_()
{
const auto num_regions = 1; // we only support single region
updateNum("PLMIXNUM", plmixnum_, num_regions);
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
updateKrnum_()
{
const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables();
updateNum("KRNUMX", krnumx_, num_regions);
updateNum("KRNUMY", krnumy_, num_regions);
updateNum("KRNUMZ", krnumz_, num_regions);
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
updateImbnum_()
{
const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables();
updateNum("IMBNUMX", imbnumx_, num_regions);
updateNum("IMBNUMY", imbnumy_, num_regions);
updateNum("IMBNUMZ", imbnumz_, num_regions);
}
template<class GridView, class FluidSystem, class Scalar>
bool EclGenericProblem<GridView,FluidSystem,Scalar>::
vapparsActive(int episodeIdx) const
{
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
return (oilVaporizationControl.getType() == OilVaporizationProperties::OilVaporization::VAPPARS);
}
template<class GridView, class FluidSystem, class Scalar>
bool EclGenericProblem<GridView,FluidSystem,Scalar>::
drsdtActive_(int episodeIdx) const
{
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
return (oilVaporizationControl.drsdtActive() && bothOilGasActive);
}
template<class GridView, class FluidSystem, class Scalar>
bool EclGenericProblem<GridView,FluidSystem,Scalar>::
drvdtActive_(int episodeIdx) const
{
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
return (oilVaporizationControl.drvdtActive() && bothOilGasActive);
}
template<class GridView, class FluidSystem, class Scalar>
bool EclGenericProblem<GridView,FluidSystem,Scalar>::
drsdtConvective_(int episodeIdx) const
{
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
return (oilVaporizationControl.drsdtConvective() && bothOilGasActive);
}
template<class GridView, class FluidSystem, class Scalar>
bool EclGenericProblem<GridView,FluidSystem,Scalar>::
beginEpisode_(bool enableExperiments,
int episodeIdx)
{
if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) {
// print some useful information in experimental mode. (the production
// simulator does this externally.)
std::ostringstream ss;
boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y");
boost::posix_time::ptime curDateTime =
boost::posix_time::from_time_t(schedule_.simTime(episodeIdx));
ss.imbue(std::locale(std::locale::classic(), facet));
ss << "Report step " << episodeIdx + 1
<< "/" << schedule_.size() - 1
<< " at day " << schedule_.seconds(episodeIdx)/(24*3600)
<< "/" << schedule_.seconds(schedule_.size() - 1)/(24*3600)
<< ", date = " << curDateTime.date()
<< "\n ";
OpmLog::info(ss.str());
}
const auto& events = schedule_[episodeIdx].events();
// react to TUNING changes
if (episodeIdx > 0 && enableTuning_ && events.hasEvent(ScheduleEvents::TUNING_CHANGE))
{
const auto& sched_state = schedule_[episodeIdx];
const auto& tuning = sched_state.tuning();
initialTimeStepSize_ = sched_state.max_next_tstep();
maxTimeStepAfterWellEvent_ = tuning.TMAXWC;
maxTimeStepSize_ = tuning.TSMAXZ;
restartShrinkFactor_ = 1./tuning.TSFCNV;
minTimeStepSize_ = tuning.TSMINZ;
return true;
}
return false;
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
beginTimeStep_(bool enableExperiments,
int episodeIdx,
int timeStepIndex,
Scalar startTime,
Scalar time,
Scalar timeStepSize,
Scalar endTime)
{
if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) {
std::ostringstream ss;
boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y");
boost::posix_time::ptime date = boost::posix_time::from_time_t(startTime) + boost::posix_time::milliseconds(static_cast<long long>(time / prefix::milli));
ss.imbue(std::locale(std::locale::classic(), facet));
ss <<"\nTime step " << timeStepIndex << ", stepsize "
<< unit::convert::to(timeStepSize, unit::day) << " days,"
<< " at day " << (double)unit::convert::to(time, unit::day)
<< "/" << (double)unit::convert::to(endTime, unit::day)
<< ", date = " << date;
OpmLog::info(ss.str());
}
// update explicit quantities between timesteps.
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
if (drsdtActive_(episodeIdx))
// DRSDT is enabled
for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx)
maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*timeStepSize;
if (drvdtActive_(episodeIdx))
// DRVDT is enabled
for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx)
maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*timeStepSize;
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
initFluidSystem_()
{
FluidSystem::initFromState(eclState_, schedule_);
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
readBlackoilExtentionsInitialConditions_(size_t numDof,
bool enableSolvent,
bool enablePolymer,
bool enablePolymerMolarWeight,
bool enableMICP)
{
if (enableSolvent) {
if (eclState_.fieldProps().has_double("SSOL"))
solventSaturation_ = eclState_.fieldProps().get_double("SSOL");
else
solventSaturation_.resize(numDof, 0.0);
}
if (enablePolymer) {
if (eclState_.fieldProps().has_double("SPOLY")) {
polymer_.concentration = eclState_.fieldProps().get_double("SPOLY");
} else {
polymer_.concentration.resize(numDof, 0.0);
}
}
if (enablePolymerMolarWeight) {
if (eclState_.fieldProps().has_double("SPOLYMW")) {
polymer_.moleWeight = eclState_.fieldProps().get_double("SPOLYMW");
} else {
polymer_.moleWeight.resize(numDof, 0.0);
}
}
if (enableMICP) {
if (eclState_.fieldProps().has_double("SMICR")) {
micp_.microbialConcentration = eclState_.fieldProps().get_double("SMICR");
} else {
micp_.microbialConcentration.resize(numDof, 0.0);
}
if (eclState_.fieldProps().has_double("SOXYG")) {
micp_.oxygenConcentration = eclState_.fieldProps().get_double("SOXYG");
} else {
micp_.oxygenConcentration.resize(numDof, 0.0);
}
if (eclState_.fieldProps().has_double("SUREA")) {
micp_.ureaConcentration = eclState_.fieldProps().get_double("SUREA");
} else {
micp_.ureaConcentration.resize(numDof, 0.0);
}
if (eclState_.fieldProps().has_double("SBIOF")) {
micp_.biofilmConcentration = eclState_.fieldProps().get_double("SBIOF");
} else {
micp_.biofilmConcentration.resize(numDof, 0.0);
}
if (eclState_.fieldProps().has_double("SCALC")) {
micp_.calciteConcentration = eclState_.fieldProps().get_double("SCALC");
} else {
micp_.calciteConcentration.resize(numDof, 0.0);
}
}
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
maxWaterSaturation(unsigned globalDofIdx) const
{
if (maxWaterSaturation_.empty())
return 0.0;
return maxWaterSaturation_[globalDofIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
minOilPressure(unsigned globalDofIdx) const
{
if (minOilPressure_.empty())
return 0.0;
return minOilPressure_[globalDofIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
overburdenPressure(unsigned elementIdx) const
{
if (overburdenPressure_.empty())
return 0.0;
return overburdenPressure_[elementIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
solventSaturation(unsigned elemIdx) const
{
if (solventSaturation_.empty())
return 0;
return solventSaturation_[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
drsdtcon(unsigned elemIdx, int episodeIdx) const
{
if (convectiveDrs_.empty())
return 0;
// The episode index is set to -1 in the initialization phase.
// Output drsdt value for index 0
episodeIdx = std::max(episodeIdx, 0);
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
int pvtRegionIdx = pvtRegionIndex(elemIdx);
return oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*convectiveDrs_[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
polymerConcentration(unsigned elemIdx) const
{
if (polymer_.concentration.empty()) {
return 0;
}
return polymer_.concentration[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
polymerMolecularWeight(const unsigned elemIdx) const
{
if (polymer_.moleWeight.empty()) {
return 0.0;
}
return polymer_.moleWeight[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
microbialConcentration(unsigned elemIdx) const
{
if (micp_.microbialConcentration.empty()) {
return 0;
}
return micp_.microbialConcentration[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
oxygenConcentration(unsigned elemIdx) const
{
if (micp_.oxygenConcentration.empty()) {
return 0;
}
return micp_.oxygenConcentration[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
ureaConcentration(unsigned elemIdx) const
{
if (micp_.ureaConcentration.empty()) {
return 0;
}
return micp_.ureaConcentration[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
biofilmConcentration(unsigned elemIdx) const
{
if (micp_.biofilmConcentration.empty()) {
return 0;
}
return micp_.biofilmConcentration[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
calciteConcentration(unsigned elemIdx) const
{
if (micp_.calciteConcentration.empty()) {
return 0;
}
return micp_.calciteConcentration[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
unsigned EclGenericProblem<GridView,FluidSystem,Scalar>::
pvtRegionIndex(unsigned elemIdx) const
{
if (pvtnum_.empty())
return 0;
return pvtnum_[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
unsigned EclGenericProblem<GridView,FluidSystem,Scalar>::
satnumRegionIndex(unsigned elemIdx) const
{
if (satnum_.empty())
return 0;
return satnum_[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
unsigned EclGenericProblem<GridView,FluidSystem,Scalar>::
miscnumRegionIndex(unsigned elemIdx) const
{
if (miscnum_.empty())
return 0;
return miscnum_[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
unsigned EclGenericProblem<GridView,FluidSystem,Scalar>::
plmixnumRegionIndex(unsigned elemIdx) const
{
if (plmixnum_.empty())
return 0;
return plmixnum_[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
maxPolymerAdsorption(unsigned elemIdx) const
{
if (polymer_.maxAdsorption.empty()) {
return 0;
}
return polymer_.maxAdsorption[elemIdx];
}
template<class GridView, class FluidSystem, class Scalar>
void EclGenericProblem<GridView,FluidSystem,Scalar>::
initDRSDT_(size_t numDof,
int episodeIdx)
{
// deal with DRSDT
unsigned ntpvt = eclState_.runspec().tabdims().getNumPVTTables();
//TODO We may want to only allocate these properties only if active.
//But since they may be activated at later time we need some more
//intrastructure to handle it
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
maxDRv_.resize(ntpvt, 1e30);
lastRv_.resize(numDof, 0.0);
maxDRs_.resize(ntpvt, 1e30);
dRsDtOnlyFreeGas_.resize(ntpvt, false);
lastRs_.resize(numDof, 0.0);
maxDRv_.resize(ntpvt, 1e30);
lastRv_.resize(numDof, 0.0);
maxOilSaturation_.resize(numDof, 0.0);
if (drsdtConvective_(episodeIdx)) {
convectiveDrs_.resize(numDof, 1.0);
}
}
}
template<class GridView, class FluidSystem, class Scalar>
bool EclGenericProblem<GridView,FluidSystem,Scalar>::
operator==(const EclGenericProblem& rhs) const
{
return this->maxOilSaturation_ == rhs.maxOilSaturation_ &&
this->maxWaterSaturation_ == rhs.maxWaterSaturation_ &&
this->minOilPressure_ == rhs.minOilPressure_ &&
this->overburdenPressure_ == rhs.overburdenPressure_ &&
this->solventSaturation_ == rhs.solventSaturation_ &&
this->polymer_ == rhs.polymer_ &&
this->micp_ == rhs.micp_ &&
this->lastRv_ == rhs.lastRv_ &&
this->maxDRv_ == rhs.maxDRv_ &&
this->convectiveDrs_ == rhs.convectiveDrs_ &&
this->lastRs_ == rhs.lastRs_ &&
this->maxDRs_ == rhs.maxDRs_ &&
this->dRsDtOnlyFreeGas_ == rhs.dRsDtOnlyFreeGas_;
}
} // namespace Opm
#endif // EWOMS_GENERIC_ECL_PROBLEM_IMPL_HH

324
ebos/eclgenericthresholdpressure.cc
View File

@@ -24,22 +24,9 @@
#include <config.h>
#include <limits>
#include <ebos/eclgenericthresholdpressure.hh>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
#include <opm/input/eclipse/EclipseState/Tables/Eqldims.hpp>
#include <opm/input/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
#include <opm/input/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
#include <dune/grid/common/mcmgmapper.hh>
#include <ebos/eclgenericthresholdpressure_impl.hh>
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/polyhedralgrid.hh>
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/3d/gridview.hh>
#include "alucartesianindexmapper.hh"
#endif // HAVE_DUNE_ALUGRID
#if HAVE_DUNE_FEM
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
@@ -47,272 +34,8 @@
#include <ebos/femcpgridcompat.hh>
#endif // HAVE_DUNE_FEM
#include <fmt/format.h>
#include <algorithm>
#include <cassert>
#include <stdexcept>
namespace Opm {
template<class Grid, class GridView, class ElementMapper, class Scalar>
EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
EclGenericThresholdPressure(const CartesianIndexMapper& cartMapper,
const GridView& gridView,
const ElementMapper& elementMapper,
const EclipseState& eclState)
: cartMapper_(cartMapper)
, gridView_(gridView)
, elementMapper_(elementMapper)
, eclState_(eclState)
{
}
template<class Grid, class GridView, class ElementMapper,class Scalar>
Scalar EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
thresholdPressure(int elem1Idx, int elem2Idx) const
{
if (!enableThresholdPressure_)
return 0.0;
// threshold pressure accross faults
if (!thpresftValues_.empty()) {
int cartElem1Idx = cartMapper_.cartesianIndex(elem1Idx);
int cartElem2Idx = cartMapper_.cartesianIndex(elem2Idx);
assert(0 <= cartElem1Idx && static_cast<int>(cartElemFaultIdx_.size()) > cartElem1Idx);
assert(0 <= cartElem2Idx && static_cast<int>(cartElemFaultIdx_.size()) > cartElem2Idx);
int fault1Idx = cartElemFaultIdx_[cartElem1Idx];
int fault2Idx = cartElemFaultIdx_[cartElem2Idx];
if (fault1Idx != -1 && fault1Idx == fault2Idx)
// inside a fault there's no threshold pressure, even accross EQUIL
// regions.
return 0.0;
if (fault1Idx != fault2Idx) {
// TODO: which value if a cell is part of multiple faults? we take
// the maximum here.
Scalar val1 = (fault1Idx >= 0) ? thpresftValues_[fault1Idx] : 0.0;
Scalar val2 = (fault2Idx >= 0) ? thpresftValues_[fault2Idx] : 0.0;
return std::max(val1, val2);
}
}
// threshold pressure accross EQUIL regions
auto equilRegion1Idx = elemEquilRegion_[elem1Idx];
auto equilRegion2Idx = elemEquilRegion_[elem2Idx];
if (equilRegion1Idx == equilRegion2Idx)
return 0.0;
return thpres_[equilRegion1Idx*numEquilRegions_ + equilRegion2Idx];
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
finishInit()
{
unsigned numElements = gridView_.size(/*codim=*/0);
const auto& simConfig = eclState_.getSimulationConfig();
enableThresholdPressure_ = simConfig.useThresholdPressure();
if (!enableThresholdPressure_)
return;
numEquilRegions_ = eclState_.getTableManager().getEqldims().getNumEquilRegions();
const decltype(numEquilRegions_) maxRegions =
std::numeric_limits<std::decay_t<decltype(elemEquilRegion_[0])>>::max();
if (numEquilRegions_ > maxRegions) {
// make sure that the index of an equilibration region can be stored
// in the vector
OPM_THROW(std::invalid_argument,
(fmt::format("The maximum number of supported "
"equilibration regions by OPM flow is {}, but "
"{} are used!",
maxRegions, numEquilRegions_)));
}
if (numEquilRegions_ > 2048) {
// warn about performance
OpmLog::warning(fmt::format("Number of equilibration regions is {}, which is "
"rather large. Note, that this might "
"have a negative impact on performance "
"and memory consumption.", numEquilRegions_));
}
// internalize the data specified using the EQLNUM keyword
const auto& fp = eclState_.fieldProps();
const auto& equilRegionData = fp.get_int("EQLNUM");
elemEquilRegion_.resize(numElements, 0);
for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
elemEquilRegion_[elemIdx] = equilRegionData[elemIdx] - 1;
}
/*
If this is a restart run the ThresholdPressure object will be active,
but it will *not* be properly initialized with numerical values. The
values must instead come from the THPRES vector in the restart file.
*/
if (simConfig.getThresholdPressure().restart())
return;
// allocate the array which specifies the threshold pressures
thpres_.resize(numEquilRegions_*numEquilRegions_, 0.0);
thpresDefault_.resize(numEquilRegions_*numEquilRegions_, 0.0);
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
applyExplicitThresholdPressures_()
{
const SimulationConfig& simConfig = eclState_.getSimulationConfig();
const auto& thpres = simConfig.getThresholdPressure();
// set the threshold pressures for all EQUIL region boundaries which have a
// intersection in the grid
for (const auto& elem : elements(gridView_, Dune::Partitions::interior)) {
for (const auto& intersection : intersections(gridView_, elem)) {
if (intersection.boundary())
continue; // ignore boundary intersections for now (TODO?)
else if (!intersection.neighbor()) //processor boundary but not domain boundary
continue;
const auto& inside = intersection.inside();
const auto& outside = intersection.outside();
unsigned insideElemIdx = elementMapper_.index(inside);
unsigned outsideElemIdx = elementMapper_.index(outside);
auto equilRegionInside = elemEquilRegion_[insideElemIdx];
auto equilRegionOutside = elemEquilRegion_[outsideElemIdx];
if (thpres.hasRegionBarrier(equilRegionInside + 1, equilRegionOutside + 1)) {
Scalar pth = 0.0;
if (thpres.hasThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1)) {
// threshold pressure explicitly specified
pth = thpres.getThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1);
}
else {
// take the threshold pressure from the initial condition
unsigned offset = equilRegionInside*numEquilRegions_ + equilRegionOutside;
pth = thpresDefault_[offset];
}
unsigned offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside;
unsigned offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside;
thpres_[offset1] = pth;
thpres_[offset2] = pth;
}
}
}
// apply threshold pressures across faults
if (thpres.ftSize() > 0)
configureThpresft_();
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
configureThpresft_()
{
// retrieve the faults collection.
const FaultCollection& faults = eclState_.getFaults();
const SimulationConfig& simConfig = eclState_.getSimulationConfig();
const auto& thpres = simConfig.getThresholdPressure();
// extract the multipliers
int numFaults = faults.size();
int numCartesianElem = eclState_.getInputGrid().getCartesianSize();
thpresftValues_.resize(numFaults, -1.0);
cartElemFaultIdx_.resize(numCartesianElem, -1);
for (size_t faultIdx = 0; faultIdx < faults.size(); faultIdx++) {
auto& fault = faults.getFault(faultIdx);
thpresftValues_[faultIdx] = thpres.getThresholdPressureFault(faultIdx);
for (const FaultFace& face : fault)
// "face" is a misnomer because the object describes a set of cell
// indices, but we go with the conventions of the parser here...
for (size_t cartElemIdx : face)
cartElemFaultIdx_[cartElemIdx] = faultIdx;
}
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
std::vector<Scalar>
EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
getRestartVector() const
{
if (!enableThresholdPressure_)
return {};
std::vector<Scalar> result(numEquilRegions_ * numEquilRegions_, 0.0);
const auto& simConfig = eclState_.getSimulationConfig();
const auto& thpres = simConfig.getThresholdPressure();
std::size_t idx = 0;
for (unsigned j = 1; j <= numEquilRegions_; ++j) {
for (unsigned i = 1; i <= numEquilRegions_; ++i, ++idx) {
if (thpres.hasRegionBarrier(i, j)) {
if (thpres.hasThresholdPressure(i, j)) {
result[idx] = thpres.getThresholdPressure(i, j);
} else {
result[idx] = this->thpresDefault_[idx];
}
}
}
}
return result;
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void
EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
logPressures()
{
if (!enableThresholdPressure_)
return;
auto lineFormat = [this](unsigned i, unsigned j, double val)
{
const auto& units = eclState_.getUnits();
return fmt::format("{:4}{:>6}{:23}{:>6}{:24}{:>11.07}{:7}{}\n",
" ", i,
" ", j,
" ", units.from_si(UnitSystem::measure::pressure, val),
" ", units.name(UnitSystem::measure::pressure));
};
auto lineFormatS = [](auto s1, auto s2, auto s3)
{
return fmt::format("{:4}{:^16}{:13}{:^9}{:21}{:^18}\n",
" ", s1, " ", s2, " ", s3);
};
std::string str = "\nLIST OF ALL NON-ZERO THRESHOLD PRESSURES\n"
"----------------------------------------\n"
"\n";
str += lineFormatS("FLOW FROM REGION", "TO REGION", "THRESHOLD PRESSURE");
str += lineFormatS(std::string(16, '-'), std::string(9, '-'), std::string(18, '-'));
const auto& simConfig = eclState_.getSimulationConfig();
const auto& thpres = simConfig.getThresholdPressure();
for (unsigned i = 1; i <= numEquilRegions_; ++i) {
for (unsigned j = (thpres.irreversible() ? 1 : i); j <= numEquilRegions_; ++j) {
if (thpres.hasRegionBarrier(i, j)) {
if (thpres.hasThresholdPressure(i, j)) {
str += lineFormat(i, j, thpres.getThresholdPressure(j, i));
} else {
std::size_t idx = (j - 1) * numEquilRegions_ + (i - 1);
str += lineFormat(i, j, this->thpresDefault_[idx]);
}
}
}
}
str += lineFormatS(std::string(16, '-'), std::string(9, '-'), std::string(18, '-'));
OpmLog::note(str);
}
#if HAVE_DUNE_FEM
template class EclGenericThresholdPressure<Dune::CpGrid,
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
@@ -331,57 +54,12 @@ template class EclGenericThresholdPressure<Dune::CpGrid,
Dune::PartitionIteratorType(4),
false>>>,
double>;
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericThresholdPressure<ALUGrid3CN,
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>>,
double>;
template class EclGenericThresholdPressure<ALUGrid3CN,
Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<
ALUGrid3CN,
Dune::PartitionIteratorType(4),
false> >,
Dune::MultipleCodimMultipleGeomTypeMapper<
Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<
ALUGrid3CN,
Dune::PartitionIteratorType(4),
false>>>,
double>;
#endif //HAVE_DUNE_ALUGRID
#else
template class EclGenericThresholdPressure<Dune::CpGrid,
Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>,
double>;
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericThresholdPressure<ALUGrid3CN,
Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,Dune::PartitionIteratorType(4)>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,Dune::PartitionIteratorType(4)>>>,
double>;
#endif //HAVE_DUNE_ALUGRID
#endif
template class EclGenericThresholdPressure<Dune::PolyhedralGrid<3,3,double>,
Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>>,
double>;
} // namespace Opm

307
ebos/eclgenericthresholdpressure_impl.hh Normal file
View File

@@ -0,0 +1,307 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#ifndef EWOMS_ECL_GENERIC_THRESHOLD_PRESSURE_IMPL_HH
#define EWOMS_ECL_GENERIC_THRESHOLD_PRESSURE_IMPL_HH
#include <ebos/eclgenericthresholdpressure.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/rangegenerators.hh>
#include <opm/common/ErrorMacros.hpp>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
#include <opm/input/eclipse/EclipseState/Tables/Eqldims.hpp>
#include <opm/input/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
#include <opm/input/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
#include <fmt/format.h>
#include <algorithm>
#include <cassert>
#include <limits>
#include <stdexcept>
namespace Opm {
template<class Grid, class GridView, class ElementMapper, class Scalar>
EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
EclGenericThresholdPressure(const CartesianIndexMapper& cartMapper,
const GridView& gridView,
const ElementMapper& elementMapper,
const EclipseState& eclState)
: cartMapper_(cartMapper)
, gridView_(gridView)
, elementMapper_(elementMapper)
, eclState_(eclState)
{
}
template<class Grid, class GridView, class ElementMapper,class Scalar>
Scalar EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
thresholdPressure(int elem1Idx, int elem2Idx) const
{
if (!enableThresholdPressure_)
return 0.0;
// threshold pressure accross faults
if (!thpresftValues_.empty()) {
int cartElem1Idx = cartMapper_.cartesianIndex(elem1Idx);
int cartElem2Idx = cartMapper_.cartesianIndex(elem2Idx);
assert(0 <= cartElem1Idx && static_cast<int>(cartElemFaultIdx_.size()) > cartElem1Idx);
assert(0 <= cartElem2Idx && static_cast<int>(cartElemFaultIdx_.size()) > cartElem2Idx);
int fault1Idx = cartElemFaultIdx_[cartElem1Idx];
int fault2Idx = cartElemFaultIdx_[cartElem2Idx];
if (fault1Idx != -1 && fault1Idx == fault2Idx)
// inside a fault there's no threshold pressure, even accross EQUIL
// regions.
return 0.0;
if (fault1Idx != fault2Idx) {
// TODO: which value if a cell is part of multiple faults? we take
// the maximum here.
Scalar val1 = (fault1Idx >= 0) ? thpresftValues_[fault1Idx] : 0.0;
Scalar val2 = (fault2Idx >= 0) ? thpresftValues_[fault2Idx] : 0.0;
return std::max(val1, val2);
}
}
// threshold pressure accross EQUIL regions
auto equilRegion1Idx = elemEquilRegion_[elem1Idx];
auto equilRegion2Idx = elemEquilRegion_[elem2Idx];
if (equilRegion1Idx == equilRegion2Idx)
return 0.0;
return thpres_[equilRegion1Idx*numEquilRegions_ + equilRegion2Idx];
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
finishInit()
{
unsigned numElements = gridView_.size(/*codim=*/0);
const auto& simConfig = eclState_.getSimulationConfig();
enableThresholdPressure_ = simConfig.useThresholdPressure();
if (!enableThresholdPressure_)
return;
numEquilRegions_ = eclState_.getTableManager().getEqldims().getNumEquilRegions();
const decltype(numEquilRegions_) maxRegions =
std::numeric_limits<std::decay_t<decltype(elemEquilRegion_[0])>>::max();
if (numEquilRegions_ > maxRegions) {
// make sure that the index of an equilibration region can be stored
// in the vector
OPM_THROW(std::invalid_argument,
(fmt::format("The maximum number of supported "
"equilibration regions by OPM flow is {}, but "
"{} are used!",
maxRegions, numEquilRegions_)));
}
if (numEquilRegions_ > 2048) {
// warn about performance
OpmLog::warning(fmt::format("Number of equilibration regions is {}, which is "
"rather large. Note, that this might "
"have a negative impact on performance "
"and memory consumption.", numEquilRegions_));
}
// internalize the data specified using the EQLNUM keyword
const auto& fp = eclState_.fieldProps();
const auto& equilRegionData = fp.get_int("EQLNUM");
elemEquilRegion_.resize(numElements, 0);
for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
elemEquilRegion_[elemIdx] = equilRegionData[elemIdx] - 1;
}
/*
If this is a restart run the ThresholdPressure object will be active,
but it will *not* be properly initialized with numerical values. The
values must instead come from the THPRES vector in the restart file.
*/
if (simConfig.getThresholdPressure().restart())
return;
// allocate the array which specifies the threshold pressures
thpres_.resize(numEquilRegions_*numEquilRegions_, 0.0);
thpresDefault_.resize(numEquilRegions_*numEquilRegions_, 0.0);
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
applyExplicitThresholdPressures_()
{
const SimulationConfig& simConfig = eclState_.getSimulationConfig();
const auto& thpres = simConfig.getThresholdPressure();
// set the threshold pressures for all EQUIL region boundaries which have a
// intersection in the grid
for (const auto& elem : elements(gridView_, Dune::Partitions::interior)) {
for (const auto& intersection : intersections(gridView_, elem)) {
if (intersection.boundary())
continue; // ignore boundary intersections for now (TODO?)
else if (!intersection.neighbor()) //processor boundary but not domain boundary
continue;
const auto& inside = intersection.inside();
const auto& outside = intersection.outside();
unsigned insideElemIdx = elementMapper_.index(inside);
unsigned outsideElemIdx = elementMapper_.index(outside);
auto equilRegionInside = elemEquilRegion_[insideElemIdx];
auto equilRegionOutside = elemEquilRegion_[outsideElemIdx];
if (thpres.hasRegionBarrier(equilRegionInside + 1, equilRegionOutside + 1)) {
Scalar pth = 0.0;
if (thpres.hasThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1)) {
// threshold pressure explicitly specified
pth = thpres.getThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1);
}
else {
// take the threshold pressure from the initial condition
unsigned offset = equilRegionInside*numEquilRegions_ + equilRegionOutside;
pth = thpresDefault_[offset];
}
unsigned offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside;
unsigned offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside;
thpres_[offset1] = pth;
thpres_[offset2] = pth;
}
}
}
// apply threshold pressures across faults
if (thpres.ftSize() > 0)
configureThpresft_();
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
configureThpresft_()
{
// retrieve the faults collection.
const FaultCollection& faults = eclState_.getFaults();
const SimulationConfig& simConfig = eclState_.getSimulationConfig();
const auto& thpres = simConfig.getThresholdPressure();
// extract the multipliers
int numFaults = faults.size();
int numCartesianElem = eclState_.getInputGrid().getCartesianSize();
thpresftValues_.resize(numFaults, -1.0);
cartElemFaultIdx_.resize(numCartesianElem, -1);
for (size_t faultIdx = 0; faultIdx < faults.size(); faultIdx++) {
auto& fault = faults.getFault(faultIdx);
thpresftValues_[faultIdx] = thpres.getThresholdPressureFault(faultIdx);
for (const FaultFace& face : fault)
// "face" is a misnomer because the object describes a set of cell
// indices, but we go with the conventions of the parser here...
for (size_t cartElemIdx : face)
cartElemFaultIdx_[cartElemIdx] = faultIdx;
}
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
std::vector<Scalar>
EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
getRestartVector() const
{
if (!enableThresholdPressure_)
return {};
std::vector<Scalar> result(numEquilRegions_ * numEquilRegions_, 0.0);
const auto& simConfig = eclState_.getSimulationConfig();
const auto& thpres = simConfig.getThresholdPressure();
std::size_t idx = 0;
for (unsigned j = 1; j <= numEquilRegions_; ++j) {
for (unsigned i = 1; i <= numEquilRegions_; ++i, ++idx) {
if (thpres.hasRegionBarrier(i, j)) {
if (thpres.hasThresholdPressure(i, j)) {
result[idx] = thpres.getThresholdPressure(i, j);
} else {
result[idx] = this->thpresDefault_[idx];
}
}
}
}
return result;
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void
EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
logPressures()
{
if (!enableThresholdPressure_)
return;
auto lineFormat = [this](unsigned i, unsigned j, double val)
{
const auto& units = eclState_.getUnits();
return fmt::format("{:4}{:>6}{:23}{:>6}{:24}{:>11.07}{:7}{}\n",
" ", i,
" ", j,
" ", units.from_si(UnitSystem::measure::pressure, val),
" ", units.name(UnitSystem::measure::pressure));
};
auto lineFormatS = [](auto s1, auto s2, auto s3)
{
return fmt::format("{:4}{:^16}{:13}{:^9}{:21}{:^18}\n",
" ", s1, " ", s2, " ", s3);
};
std::string str = "\nLIST OF ALL NON-ZERO THRESHOLD PRESSURES\n"
"----------------------------------------\n"
"\n";
str += lineFormatS("FLOW FROM REGION", "TO REGION", "THRESHOLD PRESSURE");
str += lineFormatS(std::string(16, '-'), std::string(9, '-'), std::string(18, '-'));
const auto& simConfig = eclState_.getSimulationConfig();
const auto& thpres = simConfig.getThresholdPressure();
for (unsigned i = 1; i <= numEquilRegions_; ++i) {
for (unsigned j = (thpres.irreversible() ? 1 : i); j <= numEquilRegions_; ++j) {
if (thpres.hasRegionBarrier(i, j)) {
if (thpres.hasThresholdPressure(i, j)) {
str += lineFormat(i, j, thpres.getThresholdPressure(j, i));
} else {
std::size_t idx = (j - 1) * numEquilRegions_ + (i - 1);
str += lineFormat(i, j, this->thpresDefault_[idx]);
}
}
}
}
str += lineFormatS(std::string(16, '-'), std::string(9, '-'), std::string(18, '-'));
OpmLog::note(str);
}
} // namespace Opm
#endif

49
ebos/eclgenerictracermodel.cc
View File

@@ -23,7 +23,14 @@
#include <config.h>
#include "eclgenerictracermodel_impl.hh"
#if HAVE_DUNE_FEM
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
#include <ebos/femcpgridcompat.hh>
#endif // HAVE_DUNE_FEM
namespace Opm {
#if HAVE_DUNE_FEM
template class EclGenericTracerModel<Dune::CpGrid,
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
@@ -39,54 +46,12 @@ template class EclGenericTracerModel<Dune::CpGrid,
Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false> >,
false, false>,
double>;
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericTracerModel<ALUGrid3CN, Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>>, Opm::EcfvStencil<double,Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>,false,false>,
double>;
template class EclGenericTracerModel<ALUGrid3CN,
Dune::Fem::GridPart2GridViewImpl<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false> >,
Dune::MultipleCodimMultipleGeomTypeMapper<
Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false> > >,
Opm::EcfvStencil<double, Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false> >,
false, false>,
double>;
#endif //HAVE_DUNE_ALUGRID
#else
template class EclGenericTracerModel<Dune::CpGrid,
Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>,
Opm::EcfvStencil<double,Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,false,false>,
double>;
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericTracerModel<ALUGrid3CN,
Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN, Dune::PartitionIteratorType(4)>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,
Dune::PartitionIteratorType(4)>>>,
Opm::EcfvStencil<double,Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,
Dune::PartitionIteratorType(4)>>,false,false>,
double>;
#endif //HAVE_DUNE_ALUGRID
#endif //HAVE_DUNE_FEM
template class EclGenericTracerModel<Dune::PolyhedralGrid<3,3,double>,
Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>>,
Opm::EcfvStencil<double, Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>,false,false>,
double>;
} // namespace Opm

13
ebos/eclgenerictracermodel_impl.hh
View File

@@ -39,7 +39,6 @@
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/polyhedralgrid.hh>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Phase.hpp>
@@ -53,18 +52,6 @@
#include <opm/simulators/linalg/PropertyTree.hpp>
#include <opm/simulators/linalg/FlexibleSolver.hpp>
#if HAVE_DUNE_FEM
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
#include <ebos/femcpgridcompat.hh>
#endif // HAVE_DUNE_FEM
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/3d/gridview.hh>
#include "alucartesianindexmapper.hh"
#endif // HAVE_DUNE_ALUGRID
#include <fmt/format.h>
#include <array>

618
ebos/eclgenericwriter.cc
View File

@@ -23,40 +23,12 @@
#include <config.h>
#include <ebos/eclgenericwriter.hh>
// need to include these before impl because of specializations.
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/cpgrid/GridHelpers.hpp>
#include <opm/grid/polyhedralgrid.hh>
#include <opm/grid/utility/cartesianToCompressed.hpp>
#if HAVE_DUNE_ALUGRID
#include "eclalugridvanguard.hh"
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/3d/gridview.hh>
#endif // HAVE_DUNE_ALUGRID
#include <opm/output/eclipse/EclipseIO.hpp>
#include <opm/output/eclipse/RestartValue.hpp>
#include <opm/output/eclipse/Summary.hpp>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
#include <opm/input/eclipse/Schedule/Action/State.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/Schedule/SummaryState.hpp>
#include <opm/input/eclipse/Schedule/UDQ/UDQConfig.hpp>
#include <opm/input/eclipse/Schedule/UDQ/UDQState.hpp>
#include <opm/input/eclipse/Schedule/Well/WellMatcher.hpp>
#include <opm/input/eclipse/Units/UnitSystem.hpp>
#include <dune/grid/common/mcmgmapper.hh>
#if HAVE_MPI
#include <ebos/eclmpiserializer.hh>
#endif
#include <ebos/eclgenericwriter.hh>
#include <ebos/eclgenericwriter_impl.hh>
#if HAVE_DUNE_FEM
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
@@ -64,543 +36,8 @@
#include <ebos/femcpgridcompat.hh>
#endif // HAVE_DUNE_FEM
#if HAVE_MPI
#include <mpi.h>
#endif
#include <algorithm>
#include <array>
#include <cassert>
#include <functional>
#include <map>
#include <memory>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
namespace {
/*!
* \brief Detect whether two cells are direct vertical neighbours.
*
* I.e. have the same i and j index and all cartesian cells between them
* along the vertical column are inactive.
*
* \tparam CM The type of the cartesian index mapper.
* \param cartMapper The mapper onto cartesian indices.
* \param cartesianToActive The mapping of cartesian indices to active indices.
* \param smallGlobalIndex The cartesian cell index of the cell with smaller index
* \param largeGlobalIndex The cartesian cell index of the cell with larger index
* \return True if the cells have the same i and j indices and all cartesian cells
* between them are inactive.
*/
bool directVerticalNeighbors(const std::array<int, 3>& cartDims,
const std::unordered_map<int,int>& cartesianToActive,
int smallGlobalIndex, int largeGlobalIndex)
{
assert(smallGlobalIndex <= largeGlobalIndex);
std::array<int, 3> ijk1, ijk2;
auto globalToIjk = [cartDims](int gc) {
std::array<int, 3> ijk;
ijk[0] = gc % cartDims[0];
gc /= cartDims[0];
ijk[1] = gc % cartDims[1];
ijk[2] = gc / cartDims[1];
return ijk;
};
ijk1 = globalToIjk(smallGlobalIndex);
ijk2 = globalToIjk(largeGlobalIndex);
assert(ijk2[2]>=ijk1[2]);
if ( ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1)
{
assert((largeGlobalIndex-smallGlobalIndex)%(cartDims[0]*cartDims[1])==0);
for ( int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; gi < largeGlobalIndex;
gi += cartDims[0] * cartDims[1] )
{
if ( cartesianToActive.find( gi ) != cartesianToActive.end() )
{
return false;
}
}
return true;
} else
return false;
}
std::unordered_map<std::string, Opm::data::InterRegFlowMap>
getInterRegFlowsAsMap(const Opm::EclInterRegFlowMap& map)
{
auto maps = std::unordered_map<std::string, Opm::data::InterRegFlowMap>{};
const auto& regionNames = map.names();
auto flows = map.getInterRegFlows();
const auto nmap = regionNames.size();
maps.reserve(nmap);
for (auto mapID = 0*nmap; mapID < nmap; ++mapID) {
maps.emplace(regionNames[mapID], std::move(flows[mapID]));
}
return maps;
}
struct EclWriteTasklet : public Opm::TaskletInterface
{
Opm::Action::State actionState_;
Opm::WellTestState wtestState_;
Opm::SummaryState summaryState_;
Opm::UDQState udqState_;
Opm::EclipseIO& eclIO_;
int reportStepNum_;
bool isSubStep_;
double secondsElapsed_;
Opm::RestartValue restartValue_;
bool writeDoublePrecision_;
explicit EclWriteTasklet(const Opm::Action::State& actionState,
const Opm::WellTestState& wtestState,
const Opm::SummaryState& summaryState,
const Opm::UDQState& udqState,
Opm::EclipseIO& eclIO,
int reportStepNum,
bool isSubStep,
double secondsElapsed,
Opm::RestartValue restartValue,
bool writeDoublePrecision)
: actionState_(actionState)
, wtestState_(wtestState)
, summaryState_(summaryState)
, udqState_(udqState)
, eclIO_(eclIO)
, reportStepNum_(reportStepNum)
, isSubStep_(isSubStep)
, secondsElapsed_(secondsElapsed)
, restartValue_(std::move(restartValue))
, writeDoublePrecision_(writeDoublePrecision)
{}
// callback to eclIO serial writeTimeStep method
void run()
{
this->eclIO_.writeTimeStep(this->actionState_,
this->wtestState_,
this->summaryState_,
this->udqState_,
this->reportStepNum_,
this->isSubStep_,
this->secondsElapsed_,
std::move(this->restartValue_),
this->writeDoublePrecision_);
}
};
}
namespace Opm {
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
EclGenericWriter(const Schedule& schedule,
const EclipseState& eclState,
const SummaryConfig& summaryConfig,
const Grid& grid,
const EquilGrid* equilGrid,
const GridView& gridView,
const Dune::CartesianIndexMapper<Grid>& cartMapper,
const Dune::CartesianIndexMapper<EquilGrid>* equilCartMapper,
bool enableAsyncOutput,
bool enableEsmry )
: collectToIORank_(grid,
equilGrid,
gridView,
cartMapper,
equilCartMapper,
summaryConfig.fip_regions_interreg_flow())
, grid_ (grid)
, gridView_ (gridView)
, schedule_ (schedule)
, eclState_ (eclState)
, summaryConfig_ (summaryConfig)
, cartMapper_ (cartMapper)
, equilCartMapper_(equilCartMapper)
, equilGrid_ (equilGrid)
{
if (this->collectToIORank_.isIORank()) {
this->eclIO_ = std::make_unique<EclipseIO>
(this->eclState_,
UgGridHelpers::createEclipseGrid(*equilGrid, eclState_.getInputGrid()),
this->schedule_, this->summaryConfig_, "", enableEsmry);
}
// create output thread if enabled and rank is I/O rank
// async output is enabled by default if pthread are enabled
int numWorkerThreads = 0;
if (enableAsyncOutput && collectToIORank_.isIORank()) {
numWorkerThreads = 1;
}
this->taskletRunner_.reset(new TaskletRunner(numWorkerThreads));
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
const EclipseIO& EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
eclIO() const
{
assert(eclIO_);
return *eclIO_;
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
writeInit(const std::function<unsigned int(unsigned int)>& map)
{
if (collectToIORank_.isIORank()) {
std::map<std::string, std::vector<int>> integerVectors;
if (collectToIORank_.isParallel()) {
integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks());
}
auto cartMap = cartesianToCompressed(equilGrid_->size(0), UgGridHelpers::globalCell(*equilGrid_));
eclIO_->writeInitial(computeTrans_(cartMap, map),
integerVectors,
exportNncStructure_(cartMap, map));
}
#if HAVE_MPI
if (collectToIORank_.isParallel()) {
const auto& comm = grid_.comm();
Opm::EclMpiSerializer ser(comm);
ser.broadcast(outputNnc_);
}
#endif
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
data::Solution EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
computeTrans_(const std::unordered_map<int,int>& cartesianToActive,
const std::function<unsigned int(unsigned int)>& map) const
{
const auto& cartMapper = *equilCartMapper_;
const auto& cartDims = cartMapper.cartesianDimensions();
const int globalSize = cartDims[0] * cartDims[1] * cartDims[2];
auto tranx = data::CellData {
UnitSystem::measure::transmissibility,
std::vector<double>(globalSize, 0.0),
data::TargetType::INIT
};
auto trany = tranx;
auto tranz = tranx;
using GlobalGridView = typename EquilGrid::LeafGridView;
const GlobalGridView& globalGridView = equilGrid_->leafGridView();
using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView>;
GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
for (const auto& elem : elements(globalGridView)) {
for (const auto& is : intersections(globalGridView, elem)) {
if (!is.neighbor())
continue; // intersection is on the domain boundary
unsigned c1 = globalElemMapper.index(is.inside());
unsigned c2 = globalElemMapper.index(is.outside());
if (c1 > c2)
continue; // we only need to handle each connection once, thank you.
// Ordering of compressed and uncompressed index should be the same
const int cartIdx1 = cartMapper.cartesianIndex( c1 );
const int cartIdx2 = cartMapper.cartesianIndex( c2 );
// Ordering of compressed and uncompressed index should be the same
assert(cartIdx1 <= cartIdx2);
int gc1 = std::min(cartIdx1, cartIdx2);
int gc2 = std::max(cartIdx1, cartIdx2);
// Re-ordering in case of non-empty mapping between equilGrid to grid
if (map) {
c1 = map(c1); // equilGridToGrid map
c2 = map(c2);
}
if (gc2 - gc1 == 1 && cartDims[0] > 1 ) {
tranx.data[gc1] = globalTrans().transmissibility(c1, c2);
continue; // skip other if clauses as they are false, last one needs some computation
}
if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) {
trany.data[gc1] = globalTrans().transmissibility(c1, c2);
continue; // skipt next if clause as it needs some computation
}
if ( gc2 - gc1 == cartDims[0]*cartDims[1] ||
directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2))
tranz.data[gc1] = globalTrans().transmissibility(c1, c2);
}
}
return {
{"TRANX", tranx},
{"TRANY", trany},
{"TRANZ", tranz},
};
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
std::vector<NNCdata> EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive, const std::function<unsigned int(unsigned int)>& map) const
{
std::size_t nx = eclState_.getInputGrid().getNX();
std::size_t ny = eclState_.getInputGrid().getNY();
auto nncData = eclState_.getInputNNC().input();
const auto& unitSystem = eclState_.getDeckUnitSystem();
for( const auto& entry : nncData ) {
// test whether NNC is not a neighboring connection
// cell2>=cell1 holds due to sortNncAndApplyEditnnc
assert( entry.cell2 >= entry.cell1 );
auto cellDiff = entry.cell2 - entry.cell1;
if (cellDiff != 1 && cellDiff != nx && cellDiff != nx*ny) {
auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, entry.trans);
// Eclipse ignores NNCs (with EDITNNC applied) that are small. Seems like the threshold is 1.0e-6
if ( tt >= 1.0e-6 )
outputNnc_.emplace_back(entry.cell1, entry.cell2, entry.trans);
}
}
using GlobalGridView = typename EquilGrid::LeafGridView;
const GlobalGridView& globalGridView = equilGrid_->leafGridView();
using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView>;
GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
// Cartesian index mapper for the serial I/O grid
const auto& equilCartMapper = *equilCartMapper_;
const auto& cartDims = cartMapper_.cartesianDimensions();
for (const auto& elem : elements(globalGridView)) {
for (const auto& is : intersections(globalGridView, elem)) {
if (!is.neighbor())
continue; // intersection is on the domain boundary
unsigned c1 = globalElemMapper.index(is.inside());
unsigned c2 = globalElemMapper.index(is.outside());
if (c1 > c2)
continue; // we only need to handle each connection once, thank you.
std::size_t cc1 = equilCartMapper.cartesianIndex( c1 );
std::size_t cc2 = equilCartMapper.cartesianIndex( c2 );
if ( cc2 < cc1 )
std::swap(cc1, cc2);
auto cellDiff = cc2 - cc1;
// Re-ordering in case of non-empty mapping between equilGrid to grid
if (map) {
c1 = map(c1); // equilGridToGrid map
c2 = map(c2);
}
if (cellDiff != 1 &&
cellDiff != nx &&
cellDiff != nx*ny &&
!directVerticalNeighbors(cartDims, cartesianToActive, cc1, cc2)) {
// We need to check whether an NNC for this face was also specified
// via the NNC keyword in the deck (i.e. in the first origNncSize entries.
auto t = globalTrans().transmissibility(c1, c2);
auto candidate = std::lower_bound(nncData.begin(), nncData.end(), NNCdata(cc1, cc2, 0.0));
while ( candidate != nncData.end() && candidate->cell1 == cc1
&& candidate->cell2 == cc2) {
t -= candidate->trans;
++candidate;
}
// eclipse ignores NNCs with zero transmissibility (different threshold than for NNC
// with corresponding EDITNNC above). In addition we do set small transmissibilties
// to zero when setting up the simulator. These will be ignored here, too.
auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, std::abs(t));
if ( tt > 1e-12 )
outputNnc_.push_back({cc1, cc2, t});
}
}
}
return outputNnc_;
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
doWriteOutput(const int reportStepNum,
const bool isSubStep,
data::Solution&& localCellData,
data::Wells&& localWellData,
data::GroupAndNetworkValues&& localGroupAndNetworkData,
data::Aquifers&& localAquiferData,
WellTestState&& localWTestState,
const Action::State& actionState,
const UDQState& udqState,
const SummaryState& summaryState,
const std::vector<Scalar>& thresholdPressure,
Scalar curTime,
Scalar nextStepSize,
bool doublePrecision,
bool isFlowsn,
std::array<std::pair<std::string, std::pair<std::vector<int>, std::vector<double>>>, 3>&& flowsn,
bool isFloresn,
std::array<std::pair<std::string, std::pair<std::vector<int>, std::vector<double>>>, 3>&& floresn)
{
const auto isParallel = this->collectToIORank_.isParallel();
const bool needsReordering = this->collectToIORank_.doesNeedReordering();
RestartValue restartValue {
(isParallel || needsReordering)
? this->collectToIORank_.globalCellData()
: std::move(localCellData),
isParallel ? this->collectToIORank_.globalWellData()
: std::move(localWellData),
isParallel ? this->collectToIORank_.globalGroupAndNetworkData()
: std::move(localGroupAndNetworkData),
isParallel ? this->collectToIORank_.globalAquiferData()
: std::move(localAquiferData)
};
if (eclState_.getSimulationConfig().useThresholdPressure()) {
restartValue.addExtra("THRESHPR", UnitSystem::measure::pressure,
thresholdPressure);
}
// Add suggested next timestep to extra data.
if (! isSubStep) {
restartValue.addExtra("OPMEXTRA", std::vector<double>(1, nextStepSize));
}
// Add nnc flows and flores.
if (isFlowsn) {
const auto flowsn_global = isParallel ? this->collectToIORank_.globalFlowsn() : std::move(flowsn);
for (const auto& flows : flowsn_global) {
if (flows.first.empty())
continue;
if (flows.first == "FLOGASN+") {
restartValue.addExtra(flows.first, UnitSystem::measure::gas_surface_rate, flows.second.second);
}
else {
restartValue.addExtra(flows.first, UnitSystem::measure::liquid_surface_rate, flows.second.second);
}
}
}
if (isFloresn) {
const auto floresn_global = isParallel ? this->collectToIORank_.globalFloresn() : std::move(floresn);
for (const auto& flores : floresn_global) {
if (flores.first.empty())
continue;
restartValue.addExtra(flores.first, UnitSystem::measure::rate, flores.second.second);
}
}
// first, create a tasklet to write the data for the current time
// step to disk
auto eclWriteTasklet = std::make_shared<EclWriteTasklet>(
actionState,
isParallel ? this->collectToIORank_.globalWellTestState() : std::move(localWTestState),
summaryState, udqState, *this->eclIO_,
reportStepNum, isSubStep, curTime, std::move(restartValue), doublePrecision);
// then, make sure that the previous I/O request has been completed
// and the number of incomplete tasklets does not increase between
// time steps
this->taskletRunner_->barrier();
// finally, start a new output writing job
this->taskletRunner_->dispatch(std::move(eclWriteTasklet));
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
evalSummary(const int reportStepNum,
const Scalar curTime,
const data::Wells& localWellData,
const data::WellBlockAveragePressures& localWBPData,
const data::GroupAndNetworkValues& localGroupAndNetworkData,
const std::map<int,data::AquiferData>& localAquiferData,
const std::map<std::pair<std::string, int>, double>& blockData,
const std::map<std::string, double>& miscSummaryData,
const std::map<std::string, std::vector<double>>& regionData,
const Inplace& inplace,
const Inplace& initialInPlace,
const EclInterRegFlowMap& interRegFlows,
SummaryState& summaryState,
UDQState& udqState)
{
if (collectToIORank_.isIORank()) {
const auto& summary = eclIO_->summary();
const auto& wellData = this->collectToIORank_.isParallel()
? this->collectToIORank_.globalWellData()
: localWellData;
const auto& wbpData = this->collectToIORank_.isParallel()
? this->collectToIORank_.globalWBPData()
: localWBPData;
const auto& groupAndNetworkData = this->collectToIORank_.isParallel()
? this->collectToIORank_.globalGroupAndNetworkData()
: localGroupAndNetworkData;
const auto& aquiferData = this->collectToIORank_.isParallel()
? this->collectToIORank_.globalAquiferData()
: localAquiferData;
summary.eval(summaryState,
reportStepNum,
curTime,
wellData,
wbpData,
groupAndNetworkData,
miscSummaryData,
initialInPlace,
inplace,
regionData,
blockData,
aquiferData,
getInterRegFlowsAsMap(interRegFlows));
// Off-by-one-fun: The reportStepNum argument corresponds to the
// report step these results will be written to, whereas the
// argument to UDQ function evaluation corresponds to the report
// step we are currently on.
const auto udq_step = reportStepNum - 1;
this->schedule_.getUDQConfig(udq_step)
.eval(udq_step,
this->schedule_,
this->schedule_.wellMatcher(udq_step),
summaryState,
udqState);
}
#if HAVE_MPI
if (collectToIORank_.isParallel()) {
EclMpiSerializer ser(grid_.comm());
ser.append(summaryState);
}
#endif
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
const typename EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::TransmissibilityType&
EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
globalTrans() const
{
assert (globalTrans_);
return *globalTrans_;
}
#if HAVE_DUNE_FEM
template class EclGenericWriter<Dune::CpGrid,
Dune::CpGrid,
@@ -620,59 +57,12 @@ template class EclGenericWriter<Dune::CpGrid,
Dune::PartitionIteratorType(4),
false>>>,
double>;
#ifdef HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericWriter<ALUGrid3CN,
Dune::CpGrid,
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>>,
double>;
template class EclGenericWriter<ALUGrid3CN,
Dune::CpGrid,
Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<
ALUGrid3CN,
Dune::PartitionIteratorType(4),
false>>,
Dune::MultipleCodimMultipleGeomTypeMapper<
Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<
ALUGrid3CN,
Dune::PartitionIteratorType(4),
false>>>,
double>;
#endif // HAVE_DUNE_ALUGRID
#else // !HAVE_DUNE_FEM
template class EclGenericWriter<Dune::CpGrid,
Dune::CpGrid,
Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>,
double>;
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericWriter<ALUGrid3CN,
Dune::CpGrid,
Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,Dune::PartitionIteratorType(4)>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN, Dune::PartitionIteratorType(4)>>>,
double>;
#endif // HAVE_DUNE_ALUGRID
#endif // !HAVE_DUNE_FEM
template class EclGenericWriter<Dune::PolyhedralGrid<3,3,double>,
Dune::PolyhedralGrid<3,3,double>,
Dune::GridView<Dune::PolyhedralGridViewTraits<3, 3, double, Dune::PartitionIteratorType(4)>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>>,
double>;
#endif // HAVE_DUNE_FEM
} // namespace Opm

590
ebos/eclgenericwriter_impl.hh Normal file
View File

@@ -0,0 +1,590 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#ifndef EWOMS_ECL_GENERIC_WRITER_IMPL_HH
#define EWOMS_ECL_GENERIC_WRITER_IMPL_HH
#include <ebos/eclgenericwriter.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <opm/grid/GridHelpers.hpp>
#include <opm/grid/utility/cartesianToCompressed.hpp>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
#include <opm/input/eclipse/Schedule/Action/State.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/input/eclipse/Schedule/SummaryState.hpp>
#include <opm/input/eclipse/Schedule/UDQ/UDQConfig.hpp>
#include <opm/input/eclipse/Schedule/UDQ/UDQState.hpp>
#include <opm/input/eclipse/Schedule/Well/WellMatcher.hpp>
#include <opm/input/eclipse/Units/UnitSystem.hpp>
#include <opm/output/eclipse/EclipseIO.hpp>
#include <opm/output/eclipse/RestartValue.hpp>
#include <opm/output/eclipse/Summary.hpp>
#if HAVE_MPI
#include <ebos/eclmpiserializer.hh>
#endif
#if HAVE_MPI
#include <mpi.h>
#endif
#include <algorithm>
#include <array>
#include <cassert>
#include <functional>
#include <map>
#include <memory>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
namespace {
/*!
* \brief Detect whether two cells are direct vertical neighbours.
*
* I.e. have the same i and j index and all cartesian cells between them
* along the vertical column are inactive.
*
* \tparam CM The type of the cartesian index mapper.
* \param cartMapper The mapper onto cartesian indices.
* \param cartesianToActive The mapping of cartesian indices to active indices.
* \param smallGlobalIndex The cartesian cell index of the cell with smaller index
* \param largeGlobalIndex The cartesian cell index of the cell with larger index
* \return True if the cells have the same i and j indices and all cartesian cells
* between them are inactive.
*/
bool directVerticalNeighbors(const std::array<int, 3>& cartDims,
const std::unordered_map<int,int>& cartesianToActive,
int smallGlobalIndex, int largeGlobalIndex)
{
assert(smallGlobalIndex <= largeGlobalIndex);
std::array<int, 3> ijk1, ijk2;
auto globalToIjk = [cartDims](int gc) {
std::array<int, 3> ijk;
ijk[0] = gc % cartDims[0];
gc /= cartDims[0];
ijk[1] = gc % cartDims[1];
ijk[2] = gc / cartDims[1];
return ijk;
};
ijk1 = globalToIjk(smallGlobalIndex);
ijk2 = globalToIjk(largeGlobalIndex);
assert(ijk2[2]>=ijk1[2]);
if ( ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1)
{
assert((largeGlobalIndex-smallGlobalIndex)%(cartDims[0]*cartDims[1])==0);
for ( int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; gi < largeGlobalIndex;
gi += cartDims[0] * cartDims[1] )
{
if ( cartesianToActive.find( gi ) != cartesianToActive.end() )
{
return false;
}
}
return true;
} else
return false;
}
std::unordered_map<std::string, Opm::data::InterRegFlowMap>
getInterRegFlowsAsMap(const Opm::EclInterRegFlowMap& map)
{
auto maps = std::unordered_map<std::string, Opm::data::InterRegFlowMap>{};
const auto& regionNames = map.names();
auto flows = map.getInterRegFlows();
const auto nmap = regionNames.size();
maps.reserve(nmap);
for (auto mapID = 0*nmap; mapID < nmap; ++mapID) {
maps.emplace(regionNames[mapID], std::move(flows[mapID]));
}
return maps;
}
struct EclWriteTasklet : public Opm::TaskletInterface
{
Opm::Action::State actionState_;
Opm::WellTestState wtestState_;
Opm::SummaryState summaryState_;
Opm::UDQState udqState_;
Opm::EclipseIO& eclIO_;
int reportStepNum_;
bool isSubStep_;
double secondsElapsed_;
Opm::RestartValue restartValue_;
bool writeDoublePrecision_;
explicit EclWriteTasklet(const Opm::Action::State& actionState,
const Opm::WellTestState& wtestState,
const Opm::SummaryState& summaryState,
const Opm::UDQState& udqState,
Opm::EclipseIO& eclIO,
int reportStepNum,
bool isSubStep,
double secondsElapsed,
Opm::RestartValue restartValue,
bool writeDoublePrecision)
: actionState_(actionState)
, wtestState_(wtestState)
, summaryState_(summaryState)
, udqState_(udqState)
, eclIO_(eclIO)
, reportStepNum_(reportStepNum)
, isSubStep_(isSubStep)
, secondsElapsed_(secondsElapsed)
, restartValue_(std::move(restartValue))
, writeDoublePrecision_(writeDoublePrecision)
{}
// callback to eclIO serial writeTimeStep method
void run()
{
this->eclIO_.writeTimeStep(this->actionState_,
this->wtestState_,
this->summaryState_,
this->udqState_,
this->reportStepNum_,
this->isSubStep_,
this->secondsElapsed_,
std::move(this->restartValue_),
this->writeDoublePrecision_);
}
};
}
namespace Opm {
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
EclGenericWriter(const Schedule& schedule,
const EclipseState& eclState,
const SummaryConfig& summaryConfig,
const Grid& grid,
const EquilGrid* equilGrid,
const GridView& gridView,
const Dune::CartesianIndexMapper<Grid>& cartMapper,
const Dune::CartesianIndexMapper<EquilGrid>* equilCartMapper,
bool enableAsyncOutput,
bool enableEsmry )
: collectToIORank_(grid,
equilGrid,
gridView,
cartMapper,
equilCartMapper,
summaryConfig.fip_regions_interreg_flow())
, grid_ (grid)
, gridView_ (gridView)
, schedule_ (schedule)
, eclState_ (eclState)
, summaryConfig_ (summaryConfig)
, cartMapper_ (cartMapper)
, equilCartMapper_(equilCartMapper)
, equilGrid_ (equilGrid)
{
if (this->collectToIORank_.isIORank()) {
this->eclIO_ = std::make_unique<EclipseIO>
(this->eclState_,
UgGridHelpers::createEclipseGrid(*equilGrid, eclState_.getInputGrid()),
this->schedule_, this->summaryConfig_, "", enableEsmry);
}
// create output thread if enabled and rank is I/O rank
// async output is enabled by default if pthread are enabled
int numWorkerThreads = 0;
if (enableAsyncOutput && collectToIORank_.isIORank()) {
numWorkerThreads = 1;
}
this->taskletRunner_.reset(new TaskletRunner(numWorkerThreads));
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
const EclipseIO& EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
eclIO() const
{
assert(eclIO_);
return *eclIO_;
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
writeInit(const std::function<unsigned int(unsigned int)>& map)
{
if (collectToIORank_.isIORank()) {
std::map<std::string, std::vector<int>> integerVectors;
if (collectToIORank_.isParallel()) {
integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks());
}
auto cartMap = cartesianToCompressed(equilGrid_->size(0), UgGridHelpers::globalCell(*equilGrid_));
eclIO_->writeInitial(computeTrans_(cartMap, map),
integerVectors,
exportNncStructure_(cartMap, map));
}
#if HAVE_MPI
if (collectToIORank_.isParallel()) {
const auto& comm = grid_.comm();
Opm::EclMpiSerializer ser(comm);
ser.broadcast(outputNnc_);
}
#endif
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
data::Solution EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
computeTrans_(const std::unordered_map<int,int>& cartesianToActive,
const std::function<unsigned int(unsigned int)>& map) const
{
const auto& cartMapper = *equilCartMapper_;
const auto& cartDims = cartMapper.cartesianDimensions();
const int globalSize = cartDims[0] * cartDims[1] * cartDims[2];
auto tranx = data::CellData {
UnitSystem::measure::transmissibility,
std::vector<double>(globalSize, 0.0),
data::TargetType::INIT
};
auto trany = tranx;
auto tranz = tranx;
using GlobalGridView = typename EquilGrid::LeafGridView;
const GlobalGridView& globalGridView = equilGrid_->leafGridView();
using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView>;
GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
for (const auto& elem : elements(globalGridView)) {
for (const auto& is : intersections(globalGridView, elem)) {
if (!is.neighbor())
continue; // intersection is on the domain boundary
unsigned c1 = globalElemMapper.index(is.inside());
unsigned c2 = globalElemMapper.index(is.outside());
if (c1 > c2)
continue; // we only need to handle each connection once, thank you.
// Ordering of compressed and uncompressed index should be the same
const int cartIdx1 = cartMapper.cartesianIndex( c1 );
const int cartIdx2 = cartMapper.cartesianIndex( c2 );
// Ordering of compressed and uncompressed index should be the same
assert(cartIdx1 <= cartIdx2);
int gc1 = std::min(cartIdx1, cartIdx2);
int gc2 = std::max(cartIdx1, cartIdx2);
// Re-ordering in case of non-empty mapping between equilGrid to grid
if (map) {
c1 = map(c1); // equilGridToGrid map
c2 = map(c2);
}
if (gc2 - gc1 == 1 && cartDims[0] > 1 ) {
tranx.data[gc1] = globalTrans().transmissibility(c1, c2);
continue; // skip other if clauses as they are false, last one needs some computation
}
if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) {
trany.data[gc1] = globalTrans().transmissibility(c1, c2);
continue; // skipt next if clause as it needs some computation
}
if ( gc2 - gc1 == cartDims[0]*cartDims[1] ||
directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2))
tranz.data[gc1] = globalTrans().transmissibility(c1, c2);
}
}
return {
{"TRANX", tranx},
{"TRANY", trany},
{"TRANZ", tranz},
};
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
std::vector<NNCdata> EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive, const std::function<unsigned int(unsigned int)>& map) const
{
std::size_t nx = eclState_.getInputGrid().getNX();
std::size_t ny = eclState_.getInputGrid().getNY();
auto nncData = eclState_.getInputNNC().input();
const auto& unitSystem = eclState_.getDeckUnitSystem();
for( const auto& entry : nncData ) {
// test whether NNC is not a neighboring connection
// cell2>=cell1 holds due to sortNncAndApplyEditnnc
assert( entry.cell2 >= entry.cell1 );
auto cellDiff = entry.cell2 - entry.cell1;
if (cellDiff != 1 && cellDiff != nx && cellDiff != nx*ny) {
auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, entry.trans);
// Eclipse ignores NNCs (with EDITNNC applied) that are small. Seems like the threshold is 1.0e-6
if ( tt >= 1.0e-6 )
outputNnc_.emplace_back(entry.cell1, entry.cell2, entry.trans);
}
}
using GlobalGridView = typename EquilGrid::LeafGridView;
const GlobalGridView& globalGridView = equilGrid_->leafGridView();
using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView>;
GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
// Cartesian index mapper for the serial I/O grid
const auto& equilCartMapper = *equilCartMapper_;
const auto& cartDims = cartMapper_.cartesianDimensions();
for (const auto& elem : elements(globalGridView)) {
for (const auto& is : intersections(globalGridView, elem)) {
if (!is.neighbor())
continue; // intersection is on the domain boundary
unsigned c1 = globalElemMapper.index(is.inside());
unsigned c2 = globalElemMapper.index(is.outside());
if (c1 > c2)
continue; // we only need to handle each connection once, thank you.
std::size_t cc1 = equilCartMapper.cartesianIndex( c1 );
std::size_t cc2 = equilCartMapper.cartesianIndex( c2 );
if ( cc2 < cc1 )
std::swap(cc1, cc2);
auto cellDiff = cc2 - cc1;
// Re-ordering in case of non-empty mapping between equilGrid to grid
if (map) {
c1 = map(c1); // equilGridToGrid map
c2 = map(c2);
}
if (cellDiff != 1 &&
cellDiff != nx &&
cellDiff != nx*ny &&
!directVerticalNeighbors(cartDims, cartesianToActive, cc1, cc2)) {
// We need to check whether an NNC for this face was also specified
// via the NNC keyword in the deck (i.e. in the first origNncSize entries.
auto t = globalTrans().transmissibility(c1, c2);
auto candidate = std::lower_bound(nncData.begin(), nncData.end(), NNCdata(cc1, cc2, 0.0));
while ( candidate != nncData.end() && candidate->cell1 == cc1
&& candidate->cell2 == cc2) {
t -= candidate->trans;
++candidate;
}
// eclipse ignores NNCs with zero transmissibility (different threshold than for NNC
// with corresponding EDITNNC above). In addition we do set small transmissibilties
// to zero when setting up the simulator. These will be ignored here, too.
auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, std::abs(t));
if ( tt > 1e-12 )
outputNnc_.push_back({cc1, cc2, t});
}
}
}
return outputNnc_;
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
doWriteOutput(const int reportStepNum,
const bool isSubStep,
data::Solution&& localCellData,
data::Wells&& localWellData,
data::GroupAndNetworkValues&& localGroupAndNetworkData,
data::Aquifers&& localAquiferData,
WellTestState&& localWTestState,
const Action::State& actionState,
const UDQState& udqState,
const SummaryState& summaryState,
const std::vector<Scalar>& thresholdPressure,
Scalar curTime,
Scalar nextStepSize,
bool doublePrecision,
bool isFlowsn,
std::array<std::pair<std::string, std::pair<std::vector<int>, std::vector<double>>>, 3>&& flowsn,
bool isFloresn,
std::array<std::pair<std::string, std::pair<std::vector<int>, std::vector<double>>>, 3>&& floresn)
{
const auto isParallel = this->collectToIORank_.isParallel();
const bool needsReordering = this->collectToIORank_.doesNeedReordering();
RestartValue restartValue {
(isParallel || needsReordering)
? this->collectToIORank_.globalCellData()
: std::move(localCellData),
isParallel ? this->collectToIORank_.globalWellData()
: std::move(localWellData),
isParallel ? this->collectToIORank_.globalGroupAndNetworkData()
: std::move(localGroupAndNetworkData),
isParallel ? this->collectToIORank_.globalAquiferData()
: std::move(localAquiferData)
};
if (eclState_.getSimulationConfig().useThresholdPressure()) {
restartValue.addExtra("THRESHPR", UnitSystem::measure::pressure,
thresholdPressure);
}
// Add suggested next timestep to extra data.
if (! isSubStep) {
restartValue.addExtra("OPMEXTRA", std::vector<double>(1, nextStepSize));
}
// Add nnc flows and flores.
if (isFlowsn) {
const auto flowsn_global = isParallel ? this->collectToIORank_.globalFlowsn() : std::move(flowsn);
for (const auto& flows : flowsn_global) {
if (flows.first.empty())
continue;
if (flows.first == "FLOGASN+") {
restartValue.addExtra(flows.first, UnitSystem::measure::gas_surface_rate, flows.second.second);
}
else {
restartValue.addExtra(flows.first, UnitSystem::measure::liquid_surface_rate, flows.second.second);
}
}
}
if (isFloresn) {
const auto floresn_global = isParallel ? this->collectToIORank_.globalFloresn() : std::move(floresn);
for (const auto& flores : floresn_global) {
if (flores.first.empty())
continue;
restartValue.addExtra(flores.first, UnitSystem::measure::rate, flores.second.second);
}
}
// first, create a tasklet to write the data for the current time
// step to disk
auto eclWriteTasklet = std::make_shared<EclWriteTasklet>(
actionState,
isParallel ? this->collectToIORank_.globalWellTestState() : std::move(localWTestState),
summaryState, udqState, *this->eclIO_,
reportStepNum, isSubStep, curTime, std::move(restartValue), doublePrecision);
// then, make sure that the previous I/O request has been completed
// and the number of incomplete tasklets does not increase between
// time steps
this->taskletRunner_->barrier();
// finally, start a new output writing job
this->taskletRunner_->dispatch(std::move(eclWriteTasklet));
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
evalSummary(const int reportStepNum,
const Scalar curTime,
const data::Wells& localWellData,
const data::WellBlockAveragePressures& localWBPData,
const data::GroupAndNetworkValues& localGroupAndNetworkData,
const std::map<int,data::AquiferData>& localAquiferData,
const std::map<std::pair<std::string, int>, double>& blockData,
const std::map<std::string, double>& miscSummaryData,
const std::map<std::string, std::vector<double>>& regionData,
const Inplace& inplace,
const Inplace& initialInPlace,
const EclInterRegFlowMap& interRegFlows,
SummaryState& summaryState,
UDQState& udqState)
{
if (collectToIORank_.isIORank()) {
const auto& summary = eclIO_->summary();
const auto& wellData = this->collectToIORank_.isParallel()
? this->collectToIORank_.globalWellData()
: localWellData;
const auto& wbpData = this->collectToIORank_.isParallel()
? this->collectToIORank_.globalWBPData()
: localWBPData;
const auto& groupAndNetworkData = this->collectToIORank_.isParallel()
? this->collectToIORank_.globalGroupAndNetworkData()
: localGroupAndNetworkData;
const auto& aquiferData = this->collectToIORank_.isParallel()
? this->collectToIORank_.globalAquiferData()
: localAquiferData;
summary.eval(summaryState,
reportStepNum,
curTime,
wellData,
wbpData,
groupAndNetworkData,
miscSummaryData,
initialInPlace,
inplace,
regionData,
blockData,
aquiferData,
getInterRegFlowsAsMap(interRegFlows));
// Off-by-one-fun: The reportStepNum argument corresponds to the
// report step these results will be written to, whereas the
// argument to UDQ function evaluation corresponds to the report
// step we are currently on.
const auto udq_step = reportStepNum - 1;
this->schedule_.getUDQConfig(udq_step)
.eval(udq_step,
this->schedule_,
this->schedule_.wellMatcher(udq_step),
summaryState,
udqState);
}
#if HAVE_MPI
if (collectToIORank_.isParallel()) {
EclMpiSerializer ser(grid_.comm());
ser.append(summaryState);
}
#endif
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
const typename EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::TransmissibilityType&
EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
globalTrans() const
{
assert (globalTrans_);
return *globalTrans_;
}
} // namespace Opm
#endif

43
ebos/eclpolyhedralgridvanguard.hh
View File

@@ -92,12 +92,14 @@ public:
using Grid = GetPropType<TypeTag, Properties::Grid>;
using EquilGrid = GetPropType<TypeTag, Properties::EquilGrid>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
using EquilCartesianIndexMapper = Dune::CartesianIndexMapper<EquilGrid>;
static constexpr int dimension = Grid::dimension;
static constexpr int dimensionworld = Grid::dimensionworld;
private:
using GridPointer = Grid*;
using EquilGridPointer = EquilGrid*;
using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
using CartesianIndexMapperPointer = std::unique_ptr<CartesianIndexMapper>;
public:
using TransmissibilityType = EclTransmissibility<Grid, GridView, ElementMapper,
@@ -108,11 +110,17 @@ public:
simulator_( simulator )
{
this->callImplementationInit();
// add a copy in standard vector format to fullfill new interface
const int* globalcellorg = this->grid().globalCell();
int num_cells = this->gridView().size(0);
globalcell_.resize(num_cells);
for(int i=0; i < num_cells; ++i){
globalcell_[i] = globalcellorg[i];
}
}
~EclPolyhedralGridVanguard()
{
delete grid_;
}
/*!
@@ -155,7 +163,8 @@ public:
* (For parallel simulation runs.)
*/
void loadBalance()
{ /* do nothing: PolyhedralGrid is not parallel! */ }
{ /* do nothing: PolyhedralGrid is not parallel! */
}
/*!
* \brief Returns the object which maps a global element index of the simulation grid
@@ -173,6 +182,18 @@ public:
const CartesianIndexMapper& equilCartesianIndexMapper() const
{ return *cartesianIndexMapper_; }
const std::vector<int>& globalCell()
{
return globalcell_;
}
unsigned int gridEquilIdxToGridIdx(unsigned int elemIndex) const {
return elemIndex;
}
unsigned int gridIdxToEquilGridIdx(unsigned int elemIndex) const {
return elemIndex;
}
/*!
* \brief Free the memory occupied by the global transmissibility object.
@@ -204,11 +225,17 @@ public:
return this->cellCentroids_(this->cartesianIndexMapper());
}
std::vector<int> cellPartition() const
{
// not required for this type of grid yet (only from bdaBridge??)
return {};
}
protected:
void createGrids_()
{
grid_ = new Grid(this->eclState().getInputGrid(), this->eclState().fieldProps().porv(true));
grid_ = std::make_unique<Grid>(this->eclState().getInputGrid(), this->eclState().fieldProps().porv(true));
cartesianIndexMapper_ = std::make_unique<CartesianIndexMapper>(*grid_);
this->updateGridView_();
this->updateCartesianToCompressedMapping_();
this->updateCellDepths_();
}
@@ -220,10 +247,12 @@ protected:
Simulator& simulator_;
GridPointer grid_;
CartesianIndexMapperPointer cartesianIndexMapper_;
std::unique_ptr<Grid> grid_;
std::unique_ptr<CartesianIndexMapper> cartesianIndexMapper_;
//CartesianIndexMapperPointer cartesianIndexMapper_;
std::unordered_set<std::string> defunctWellNames_;
std::vector<int> globalcell_;
};
} // namespace Opm

45
ebos/eclproblem.hh
View File

@@ -34,6 +34,7 @@
#include <ebos/eclactionhandler.hh>
#include <ebos/eclbaseaquifermodel.hh>
#include <ebos/eclcpgridvanguard.hh>
#include <ebos/ecldummygradientcalculator.hh>
#include <ebos/eclequilinitializer.hh>
#include <ebos/eclfluxmodule.hh>
@@ -84,15 +85,6 @@
#include <opm/common/OpmLog/OpmLog.hpp>
#if USE_ALUGRID
#define DISABLE_ALUGRID_SFC_ORDERING 1
#include <ebos/eclalugridvanguard.hh>
#elif USE_POLYHEDRALGRID
#include <ebos/eclpolyhedralgridvanguard.hh>
#else
#include <ebos/eclcpgridvanguard.hh>
#endif
#include <algorithm>
#include <functional>
#include <set>
@@ -108,19 +100,9 @@ namespace Opm::Properties {
namespace TTag {
#if USE_ALUGRID
struct EclBaseProblem {
using InheritsFrom = std::tuple<VtkEclTracer, EclOutputBlackOil, EclAluGridVanguard>;
};
#elif USE_POLYHEDRALGRID
struct EclBaseProblem {
using InheritsFrom = std::tuple<VtkEclTracer, EclOutputBlackOil, EclPolyhedralGridVanguard>;
};
#else
struct EclBaseProblem {
using InheritsFrom = std::tuple<VtkEclTracer, EclOutputBlackOil, EclCpGridVanguard>;
};
#endif
}
// The class which deals with ECL wells
@@ -207,6 +189,7 @@ struct NumPressurePointsEquil {
};
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::EclBaseProblem> {
@@ -531,6 +514,10 @@ struct EnableTemperature<TypeTag, TTag::EclBaseProblem> {
static constexpr bool value = true;
};
template<class TypeTag>
struct EnableMech<TypeTag, TTag::EclBaseProblem> {
static constexpr bool value = false;
};
// disable all extensions supported by black oil model. this should not really be
// necessary but it makes things a bit more explicit
template<class TypeTag>
@@ -891,12 +878,9 @@ public:
readThermalParameters_();
// Re-ordering in case of ALUGrid
std::function<unsigned int(unsigned int)> gridToEquilGrid;
#if USE_ALUGRID
gridToEquilGrid = [&simulator](unsigned int i) {
std::function<unsigned int(unsigned int)> gridToEquilGrid = [&simulator](unsigned int i) {
return simulator.vanguard().gridIdxToEquilGridIdx(i);
};
#endif // USE_ALUGRID
transmissibilities_.finishInit(gridToEquilGrid);
const auto& initconfig = eclState.getInitConfig();
@@ -936,12 +920,9 @@ public:
eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());
// Re-ordering in case of ALUGrid
std::function<unsigned int(unsigned int)> equilGridToGrid;
#if USE_ALUGRID
equilGridToGrid = [&simulator](unsigned int i) {
std::function<unsigned int(unsigned int)> equilGridToGrid = [&simulator](unsigned int i) {
return simulator.vanguard().gridEquilIdxToGridIdx(i);
};
#endif // USE_ALUGRID
eclWriter_->writeInit(equilGridToGrid);
}
@@ -1030,12 +1011,9 @@ public:
eclBroadcast(cc, eclState.getTransMult() );
// Re-ordering in case of ALUGrid
std::function<unsigned int(unsigned int)> equilGridToGrid;
#if USE_ALUGRID
equilGridToGrid = [&simulator](unsigned int i) {
std::function<unsigned int(unsigned int)> equilGridToGrid = [&simulator](unsigned int i) {
return simulator.vanguard().gridEquilIdxToGridIdx(i);
};
#endif // USE_ALUGRID
// re-compute all quantities which may possibly be affected.
transmissibilities_.update(true, equilGridToGrid);
@@ -1177,12 +1155,9 @@ public:
int episodeIdx = this->episodeIndex();
// Re-ordering in case of Alugrid
std::function<unsigned int(unsigned int)> gridToEquilGrid;
#if USE_ALUGRID
gridToEquilGrid = [&simulator](unsigned int i) {
std::function<unsigned int(unsigned int)> gridToEquilGrid = [&simulator](unsigned int i) {
return simulator.vanguard().gridIdxToEquilGridIdx(i);
};
#endif // USE_ALUGRID
std::function<void(bool)> transUp =
[this,gridToEquilGrid](bool global) {

1122
ebos/ecltransmissibility.cc
View File

File diff suppressed because it is too large Load Diff

1097
ebos/ecltransmissibility_impl.hh Normal file
View File

File diff suppressed because it is too large Load Diff

46
ebos/equil/initstateequil.cc
View File

@@ -21,13 +21,22 @@
copyright holders.
*/
#include "config.h"
#include "initstateequil.hh"
#include "initstateequil_impl.hh"
#include <config.h>
#include <ebos/equil/initstateequil.hh>
#include <ebos/equil/initstateequil_impl.hh>
#include <opm/grid/CpGrid.hpp>
#if HAVE_DUNE_FEM
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
#include <ebos/femcpgridcompat.hh>
#endif
namespace Opm {
namespace EQUIL {
namespace DeckDependent {
#if HAVE_DUNE_FEM
using GridView = Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<
@@ -59,37 +68,6 @@ template InitialStateComputer<BlackOilFluidSystem<double>,
const double,
const int,
const bool);
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGridComm = Dune::ALUGridMPIComm;
#else
using ALUGridComm = Dune::ALUGridNoComm;
#endif //HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, ALUGridComm>;
using ALUGridView = Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN, Dune::PartitionIteratorType(4)>>;
using ALUGridMapper = Dune::MultipleCodimMultipleGeomTypeMapper<ALUGridView>;
template class InitialStateComputer<BlackOilFluidSystem<double>,
ALUGrid3CN,
ALUGridView,
ALUGridMapper,
Dune::CartesianIndexMapper<ALUGrid3CN>>;
template InitialStateComputer<BlackOilFluidSystem<double>,
ALUGrid3CN,
ALUGridView,
ALUGridMapper,
Dune::CartesianIndexMapper<ALUGrid3CN>>::
InitialStateComputer(MatLaw&,
const EclipseState&,
const ALUGrid3CN&,
const ALUGridView&,
const Dune::CartesianIndexMapper<ALUGrid3CN>&,
const double,
const int,
const bool);
#endif //HAVE_DUNE_ALUGRID
} // namespace DeckDependent
namespace Details {

12
ebos/equil/initstateequil_impl.hh
View File

@@ -29,7 +29,6 @@
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/utility/RegionMapping.hpp>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
@@ -54,17 +53,6 @@
#include <limits>
#include <stdexcept>
#if HAVE_DUNE_FEM
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
#include <ebos/femcpgridcompat.hh>
#endif
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/3d/gridview.hh>
#endif // HAVE_DUNE_ALUGRID
namespace Opm {
namespace EQUIL {

38
flow/flow_blackoil_alugrid.cpp
View File

@@ -16,9 +16,26 @@
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <config.h>
#include <dune/alugrid/grid.hh>
#include <ebos/eclalugridvanguard.hh>
#include <opm/simulators/flow/Main.hpp>
// for equilgrid in writer
// need to include this before eclgenericwriter_impl.hh due to specializations.
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/cpgrid/GridHelpers.hpp>
// these are not explicitly instanced in library
#include <ebos/collecttoiorank_impl.hh>
#include <ebos/eclgenericproblem_impl.hh>
#include <ebos/eclgenericthresholdpressure_impl.hh>
#include <ebos/eclgenerictracermodel_impl.hh>
#include <ebos/eclgenericwriter_impl.hh>
#include <ebos/ecltransmissibility_impl.hh>
#include <ebos/equil/initstateequil_impl.hh>
namespace Opm {
namespace Properties {
namespace TTag {
@@ -26,6 +43,25 @@ struct EclFlowProblemAlugrid {
using InheritsFrom = std::tuple<EclFlowProblem>;
};
}
template<class TypeTag>
struct Grid<TypeTag, TTag::EclFlowProblemAlugrid> {
static const int dim = 3;
#if HAVE_MPI
using type = Dune::ALUGrid<dim, dim, Dune::cube, Dune::nonconforming,Dune::ALUGridMPIComm>;
#else
using type = Dune::ALUGrid<dim, dim, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif
};
// alugrid need cp grid as equilgrid
template<class TypeTag>
struct EquilGrid<TypeTag, TTag::EclFlowProblemAlugrid> {
using type = Dune::CpGrid;
};
template<class TypeTag>
struct Vanguard<TypeTag, TTag::EclFlowProblemAlugrid> {
using type = Opm::EclAluGridVanguard<TypeTag>;
};
template<class TypeTag>
struct EclEnableAquifers<TypeTag, TTag::EclFlowProblemAlugrid> {
static constexpr bool value = false;

25
flow/flow_blackoil_poly.cpp
View File

@@ -1,25 +0,0 @@
/*
Copyright 2020, NORCE AS
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 3 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/>.
*/
#include "config.h"
#include <flow/flow_ebos_blackoil_poly.hpp>
int main(int argc, char** argv)
{
return Opm::flowEbosBlackoilPolyMain(argc, argv);
}

78
flow/flow_blackoil_polyhedralgrid.cpp Normal file
View File

@@ -0,0 +1,78 @@
/*
Copyright 2020, NORCE AS
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 3 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/>.
*/
#include <config.h>
#include <ebos/eclpolyhedralgridvanguard.hh>
#include <opm/grid/polyhedralgrid.hh>
#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
#include <opm/models/discretization/common/tpfalinearizer.hh>
#include <opm/simulators/flow/Main.hpp>
// these are not explicitly instanced in library
#include <ebos/collecttoiorank_impl.hh>
#include <ebos/eclgenericproblem_impl.hh>
#include <ebos/eclgenericthresholdpressure_impl.hh>
#include <ebos/eclgenerictracermodel_impl.hh>
#include <ebos/ecltransmissibility_impl.hh>
#include <ebos/eclgenericwriter_impl.hh>
#include <ebos/equil/initstateequil_impl.hh>
namespace Opm {
namespace Properties {
namespace TTag {
struct EclFlowProblemPoly {
using InheritsFrom = std::tuple<EclFlowProblem>;
};
}
template<class TypeTag>
struct Linearizer<TypeTag, TTag::EclFlowProblemPoly> { using type = TpfaLinearizer<TypeTag>; };
template<class TypeTag>
struct LocalResidual<TypeTag, TTag::EclFlowProblemPoly> { using type = BlackOilLocalResidualTPFA<TypeTag>; };
template<class TypeTag>
struct EnableDiffusion<TypeTag, TTag::EclFlowProblemPoly> { static constexpr bool value = false; };
template<class TypeTag>
struct Grid<TypeTag, TTag::EclFlowProblemPoly> {
using type = Dune::PolyhedralGrid<3, 3>;
};
template<class TypeTag>
struct EquilGrid<TypeTag, TTag::EclFlowProblemPoly> {
//using type = Dune::CpGrid;
using type = GetPropType<TypeTag, Properties::Grid>;
};
template<class TypeTag>
struct Vanguard<TypeTag, TTag::EclFlowProblemPoly> {
using type = Opm::EclPolyhedralGridVanguard<TypeTag>;
};
}
}
int main(int argc, char** argv)
{
using TypeTag = Opm::Properties::TTag::EclFlowProblemPoly;
auto mainObject = std::make_unique<Opm::Main>(argc, argv);
auto ret = mainObject->runStatic<TypeTag>();
// Destruct mainObject as the destructor calls MPI_Finalize!
mainObject.reset();
return ret;
}

41
flow/flow_ebos_blackoil_poly.cpp
View File

@@ -1,41 +0,0 @@
/*
Copyright 2020, NORCE AS
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 3 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/>.
*/
#include "config.h"
#include <opm/simulators/flow/Main.hpp>
namespace Opm {
namespace Properties {
namespace TTag {
struct EclFlowProblemPoly {
using InheritsFrom = std::tuple<EclFlowProblem>;
};
}
}
int flowEbosBlackoilPolyMain(int argc, char** argv)
{
using TypeTag = Opm::Properties::TTag::EclFlowProblemPoly;
auto mainObject = std::make_unique<Opm::Main>(argc, argv);
auto ret = mainObject->runStatic<TypeTag>();
// Destruct mainObject as the destructor calls MPI_Finalize!
mainObject.reset();
return ret;
}
}

29
flow/flow_ebos_blackoil_poly.hpp
View File

@@ -1,29 +0,0 @@
/*
Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics.
Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015, 2017 IRIS AS
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 3 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/>.
*/
#ifndef FLOW_BLACKOIL_POLY_HPP
#define FLOW_BLACKOIL_POLY_HPP
namespace Opm {
int flowEbosBlackoilPolyMain(int argc, char** argv);
}
#endif

8
opm/simulators/aquifers/BlackoilAquiferModel.hpp
View File

@@ -37,10 +37,8 @@
#include <opm/simulators/aquifers/AquiferNumerical.hpp>
#include <opm/grid/CpGrid.hpp>
#ifdef USE_POLYHEDRALGRID
#include <opm/grid/polyhedralgrid.hh>
#endif
#if USE_ALUGRID
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#endif
@@ -65,14 +63,12 @@ class SupportsFaceTag<Dune::CpGrid>
{};
#ifdef USE_POLYHEDRALGRID
template<>
class SupportsFaceTag<Dune::PolyhedralGrid<3, 3>>
: public std::bool_constant<true>
{};
#endif
#if USE_ALUGRID
#if HAVE_DUNE_ALUGRID
template<>
class SupportsFaceTag<Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming>>
: public std::bool_constant<true>

6
opm/simulators/linalg/ISTLSolverEbos.cpp
View File

@@ -43,6 +43,8 @@
#include <ebos/alucartesianindexmapper.hh>
#endif // HAVE_DUNE_ALUGRID
#include <opm/grid/polyhedralgrid.hh>
namespace Opm {
namespace detail {
@@ -387,7 +389,7 @@ using CommunicationType = Dune::CollectiveCommunication<int>;
const std::vector<Well>&, \
const std::vector<int>&, \
const size_t, const bool);
using PolyHedralGrid3D = Dune::PolyhedralGrid<3, 3>;
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
@@ -398,11 +400,13 @@ using CommunicationType = Dune::CollectiveCommunication<int>;
template struct BdaSolverInfo<BM<Dim>,BV<Dim>>; \
INSTANCE_GRID(Dim,Dune::CpGrid) \
INSTANCE_GRID(Dim,ALUGrid3CN) \
INSTANCE_GRID(Dim,PolyHedralGrid3D) \
INSTANCE_FLEX(Dim)
#else
#define INSTANCE(Dim) \
template struct BdaSolverInfo<BM<Dim>,BV<Dim>>; \
INSTANCE_GRID(Dim,Dune::CpGrid) \
INSTANCE_GRID(Dim,PolyHedralGrid3D) \
INSTANCE_FLEX(Dim)
#endif
#else