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

Add FIP, INIT and EGRID output

Browse Source 
Writes INIT and EGRID files initially
Adds Fip to summary output if required.
Output Fip to log (.PRT) if Opm-logger is set
This commit is contained in:
Tor Harald Sandve 2018-01-12 15:32:43 +01:00
parent 9e7857ec1c
commit f62498440d
4 changed files with 551 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
ebos/collecttoiorank.hh
View File

@ -93,16 +93,19 @@ namespace Ewoms
IndexMapType& localIndexMap_;
IndexMapStorageType& indexMaps_;
std::map< const int, const int > globalPosition_;
std::vector<int>& ranks_;
public:
DistributeIndexMapping( const std::vector<int>& globalIndex,
const std::vector<int>& distributedGlobalIndex,
IndexMapType& localIndexMap,
IndexMapStorageType& indexMaps )
IndexMapStorageType& indexMaps,
std::vector<int>& ranks)
: distributedGlobalIndex_( distributedGlobalIndex ),
localIndexMap_( localIndexMap ),
indexMaps_( indexMaps ),
globalPosition_()
globalPosition_(),
ranks_(ranks)
{
const size_t size = globalIndex.size();
// create mapping globalIndex --> localIndex
@ -114,6 +117,7 @@ namespace Ewoms
// we need to create a mapping from local to global
if( ! indexMaps_.empty() )
{
ranks_.resize( size, -1);
// for the ioRank create a localIndex to index in global state map
IndexMapType& indexMap = indexMaps_.back();
const size_t localSize = localIndexMap_.size();
@ -122,6 +126,7 @@ namespace Ewoms
{
const int id = distributedGlobalIndex_[ localIndexMap_[ i ] ];
indexMap[ i ] = globalPosition_[ id ] ;
ranks_[ indexMap[ i ] ] = ioRank;
}
}
}
@ -160,6 +165,7 @@ namespace Ewoms
buffer.read( globalId );
assert( globalPosition_.find( globalId ) != globalPosition_.end() );
indexMap[ index ] = globalPosition_[ globalId ];
ranks_[ indexMap[ index ] ] = link + 1;
}
}
};
@ -265,7 +271,7 @@ namespace Ewoms
indexMaps_.resize( comm.size() );
// distribute global id's to io rank for later association of dof's
DistributeIndexMapping distIndexMapping( globalCartesianIndex_, distributedCartesianIndex, localIndexMap_, indexMaps_ );
DistributeIndexMapping distIndexMapping( globalCartesianIndex_, distributedCartesianIndex, localIndexMap_, indexMaps_, globalRanks_);
toIORankComm_.exchange( distIndexMapping );
}
}
@ -435,14 +441,14 @@ namespace Ewoms
return toIORankComm_.size() > 1;
}
int localIdxToGlobalIdx(const unsigned localIdx) {
int localIdxToGlobalIdx(const unsigned localIdx) const
{
if ( ! isParallel() )
{
return localIdx;
}
// the last indexMap is the local one
IndexMapType& indexMap = indexMaps_.back();
const IndexMapType& indexMap = indexMaps_.back();
if( indexMap.empty() )
OPM_THROW(std::logic_error,"index map is not created on this rank");
@ -454,11 +460,17 @@ namespace Ewoms
size_t numCells () const { return globalCartesianIndex_.size(); }
const std::vector<int>& globalRanks() const
{
return globalRanks_;
}
protected:
P2PCommunicatorType toIORankComm_;
IndexMapType globalCartesianIndex_;
IndexMapType localIndexMap_;
IndexMapStorageType indexMaps_;
std::vector<int> globalRanks_;
Opm::data::Solution globalCellData_;
};

375
ebos/ecloutputblackoilmodule.hh
View File

@ -34,6 +34,7 @@
#include <opm/common/Valgrind.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/parser/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
#include <opm/output/data/Cells.hpp>
#include <opm/output/eclipse/EclipseIO.hpp>
@ -42,6 +43,8 @@
#include <type_traits>
#define ENUM_TO_STR(ENUM) std::string(#ENUM)
namespace Ewoms {
namespace Properties {
// create new type tag for the Ecl-output
@ -81,10 +84,28 @@ class EclOutputBlackOilModule
typedef std::vector<Scalar> ScalarBuffer;
struct FIPDataType {
enum FipId {
WIP = waterPhaseIdx,
OIP = oilPhaseIdx,
GIP = gasPhaseIdx,
OIPL = 3,
OIPG = 4,
GIPL = 5,
GIPG = 6,
PV = 7, //< Pore volume
PAV = 8
};
static const int fipValues = PAV + 1 ;
};
public:
EclOutputBlackOilModule(const Simulator& simulator)
: simulator_(simulator)
{ }
{
createLocalFipnum_();
}
/*!
* \brief Allocate memory for the scalar fields we would like to
@ -101,37 +122,61 @@ public:
keyValue.second = restartConfig.getKeyword(keyValue.first, reportStepNum);
}
const Opm::SummaryConfig summaryConfig = simulator_.gridManager().summaryConfig();
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx))
continue;
saturation_[phaseIdx].resize(bufferSize,0.0);
if (!substep || (phaseIdx == waterPhaseIdx && summaryConfig.require3DField("SWAT") )
|| (phaseIdx == gasPhaseIdx && summaryConfig.require3DField("SGAS") ) )
saturation_[phaseIdx].resize(bufferSize,0.0);
}
oilPressure_.resize(bufferSize,0.0);
temperature_.resize(bufferSize,0.0);
if (!substep || summaryConfig.require3DField("PRESSURE"))
oilPressure_.resize(bufferSize,0.0);
if (!substep || summaryConfig.require3DField("TEMP"))
temperature_.resize(bufferSize,0.0);
// Output the same as legacy
// TODO: Only needed if DISGAS or VAPOIL
if (true)
rs_.resize(bufferSize,0.0);
if (true)
rv_.resize(bufferSize,0.0);
if (true) {
if (!substep || summaryConfig.require3DField("RS"))
rs_.resize(bufferSize,0.0);
}
if (true) {
if (!substep || summaryConfig.require3DField("RV"))
rv_.resize(bufferSize,0.0);
}
if (GET_PROP_VALUE(TypeTag, EnableSolvent)) {
sSol_.resize(bufferSize,0.0);
if (!substep || summaryConfig.require3DField("SSOL"))
sSol_.resize(bufferSize,0.0);
}
if (GET_PROP_VALUE(TypeTag, EnablePolymer)) {
cPolymer_.resize(bufferSize,0.0);
if (!substep || summaryConfig.require3DField("POLYMER"))
cPolymer_.resize(bufferSize,0.0);
}
// Fluid in place
for (int i = 0; i<FIPDataType::fipValues; i++) {
//std::cout << stringOfEnumIndex_(i) << std::endl;
if (!substep || summaryConfig.require3DField(stringOfEnumIndex_(i))) {
//std::cout << "required " << stringOfEnumIndex_(i) << std::endl;
fip_[i].resize(bufferSize, 0.0);
}
}
if (!substep) {
// Output the same as legacy
// TODO: Only needed if Vappars or hysteresis.
soMax_.resize(bufferSize,0.0);
pcSwMdcOw_.resize(bufferSize,0.0);
krnSwMdcOw_.resize(bufferSize,0.0);
pcSwMdcGo_.resize(bufferSize,0.0);
krnSwMdcGo_.resize(bufferSize,0.0);
soMax_.resize(bufferSize,0.0);
pcSwMdcOw_.resize(bufferSize,0.0);
krnSwMdcOw_.resize(bufferSize,0.0);
pcSwMdcGo_.resize(bufferSize,0.0);
krnSwMdcGo_.resize(bufferSize,0.0);
}
// Only provide RESTART_AUXILIARY if it is asked for by the user
if (!restartConfig.getWriteRestartFile(reportStepNum) || substep)
@ -253,6 +298,7 @@ public:
for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
typedef typename std::remove_const<typename std::remove_reference<decltype(fs)>::type>::type FluidState;
unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
unsigned pvtRegionIdx = elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex();
@ -438,6 +484,62 @@ public:
intQuants.pvtRegionIndex());
}
// FIP
// calculate the pore volume of the current cell. Note that the porosity
// returned by the intensive quantities is defined as the ratio of pore
// space to total cell volume and includes all pressure dependent (->
// rock compressibility) and static modifiers (MULTPV, MULTREGP, NTG,
// PORV, MINPV and friends). Also note that because of this, the porosity
// returned by the intensive quantities can be outside of the physical
// range [0, 1] in pathetic cases.
const double pv =
elemCtx.simulator().model().dofTotalVolume(globalDofIdx)
* intQuants.porosity().value();
if (fip_[FIPDataType::PV].size() > 0)
fip_[FIPDataType::PV][globalDofIdx] = pv;
Scalar fip[FluidSystem::numPhases];
for (unsigned phase = 0; phase < FluidSystem::numPhases; ++phase) {
if (!FluidSystem::phaseIsActive(phase)) {
continue;
}
const double b = Toolbox::value(fs.invB(phase));
const double s = Toolbox::value(fs.saturation(phase));
fip[phase] = b * s * pv;
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && fip_[FIPDataType::OIP].size() > 0)
fip_[FIPDataType::OIP][globalDofIdx] = fip[FluidSystem::oilPhaseIdx];
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && fip_[FIPDataType::GIP].size() > 0)
fip_[FIPDataType::GIP][globalDofIdx] = fip[FluidSystem::gasPhaseIdx];
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && fip_[FIPDataType::WIP].size() > 0)
fip_[FIPDataType::WIP][globalDofIdx] = fip[FluidSystem::waterPhaseIdx];
// Store the pure oil and gas FIP
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && fip_[FIPDataType::OIPL].size() > 0)
fip_[FIPDataType::OIPL][globalDofIdx] = fip[FluidSystem::oilPhaseIdx];
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && fip_[FIPDataType::GIPG].size() > 0)
fip_[FIPDataType::GIPG][globalDofIdx] = fip[FluidSystem::gasPhaseIdx];
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
// Gas dissolved in oil and vaporized oil
Scalar gipl = Toolbox::value(fs.Rs()) * fip[FluidSystem::oilPhaseIdx];
Scalar oipg = Toolbox::value(fs.Rv()) * fip[FluidSystem::gasPhaseIdx];
if (fip_[FIPDataType::GIPG].size() > 0)
fip_[FIPDataType::GIPL][globalDofIdx] = gipl;
if (fip_[FIPDataType::OIPG].size() > 0)
fip_[FIPDataType::OIPG][globalDofIdx] = oipg;
// Add dissolved gas and vaporized oil to total FIP
if (fip_[FIPDataType::OIP].size() > 0)
fip_[FIPDataType::OIP][globalDofIdx] += oipg;
if (fip_[FIPDataType::GIP].size() > 0)
fip_[FIPDataType::GIP][globalDofIdx] += gipl;
}
}
}
@ -519,7 +621,7 @@ public:
return;
if ( oilPressure_.size() > 0 ) {
sol.insert( "PRESSURE", Opm::UnitSystem::measure::pressure, std::move(oilPressure_), Opm::data::TargetType::RESTART_SOLUTION);
sol.insert( "PRESSURE", Opm::UnitSystem::measure::pressure, oilPressure_, Opm::data::TargetType::RESTART_SOLUTION);
}
if ( temperature_.size() > 0 ) {
@ -527,10 +629,10 @@ public:
}
if( FluidSystem::phaseIsActive(waterPhaseIdx) && saturation_[waterPhaseIdx].size() > 0 ) {
sol.insert( "SWAT", Opm::UnitSystem::measure::identity, std::move(saturation_[waterPhaseIdx]), Opm::data::TargetType::RESTART_SOLUTION );
sol.insert( "SWAT", Opm::UnitSystem::measure::identity, saturation_[waterPhaseIdx], Opm::data::TargetType::RESTART_SOLUTION );
}
if( FluidSystem::phaseIsActive(gasPhaseIdx) && saturation_[gasPhaseIdx].size() > 0) {
sol.insert( "SGAS", Opm::UnitSystem::measure::identity, std::move(saturation_[gasPhaseIdx]), Opm::data::TargetType::RESTART_SOLUTION );
sol.insert( "SGAS", Opm::UnitSystem::measure::identity, saturation_[gasPhaseIdx], Opm::data::TargetType::RESTART_SOLUTION );
}
if ( gasDissolutionFactor_.size() > 0 ) {
@ -614,7 +716,79 @@ public:
if (bubblePointPressure_.size() > 0)
sol.insert ("PBUB", Opm::UnitSystem::measure::pressure, std::move(bubblePointPressure_), Opm::data::TargetType::RESTART_AUXILIARY);
}
// Summary FIP output
// Fluid in place
for (int i = 0; i<FIPDataType::fipValues; i++) {
if (fip_[i].size() > 0) {
sol.insert(stringOfEnumIndex_(i),
Opm::UnitSystem::measure::volume,
fip_[i] ,
Opm::data::TargetType::SUMMARY);
}
}
}
// write Fip to output log
void outputFIPLog() {
const auto& comm = simulator_.gridView().comm();
size_t ntFip = *std::max_element(fipnum_.begin(), fipnum_.end());
ntFip = comm.max(ntFip);
// sum values over each region
ScalarBuffer regionValues[FIPDataType::fipValues];
for (int i = 0; i<FIPDataType::fipValues; i++) {
regionValues[i] = FIPTotals_(fip_[i], fipnum_, ntFip);
if (isIORank_() && origRegionValues_[i].empty())
origRegionValues_[i] = regionValues[i];
}
// sum all region values to compute the field total
std::vector<int> fieldNum(ntFip, 1);
ScalarBuffer totalValues(FIPDataType::fipValues,0.0);
for (int i = 0; i<FIPDataType::fipValues; i++) {
bool comunicateSum = false; // the regionValues are already summed over all ranks.
const ScalarBuffer& tmp = FIPTotals_(regionValues[i], fieldNum, 1, comunicateSum);
totalValues[i] = tmp[0]; //
}
// compute the hydrocarbon averaged pressure over the field.
totalValues[FIPDataType::PAV] = FipPav_(-1);
// convert units and output field values
// the original Fip values are stored on the first step
// TODO: Store initial Fip in the init file and restore them
// and use them here.
if ( isIORank_() ) {
FIPUnitConvert_(simulator_.gridManager().eclState().getUnits(), totalValues);
if (origTotalValues_.empty())
origTotalValues_ = totalValues;
outputRegionFluidInPlace_(origTotalValues_, totalValues, simulator_.gridManager().eclState().getUnits(), 0);
}
// Do the same on each region.
for (size_t reg = 0; reg < ntFip; ++reg ) {
ScalarBuffer tmp(FIPDataType::fipValues,0.0);
for (int i = 0; i<FIPDataType::fipValues; i++) {
tmp[i] = regionValues[i][reg];
}
tmp[FIPDataType::PAV] = FipPav_(reg+1);
if ( isIORank_() ) {
ScalarBuffer tmpO(FIPDataType::fipValues,0.0);
for (int i = 0; i<FIPDataType::fipValues; i++) {
tmpO[i] = origRegionValues_[i][reg];
}
FIPUnitConvert_(simulator_.gridManager().eclState().getUnits(), tmp);
FIPUnitConvert_(simulator_.gridManager().eclState().getUnits(), tmpO);
outputRegionFluidInPlace_(tmpO, tmp, simulator_.gridManager().eclState().getUnits() , reg + 1);
}
}
}
void setRestart(const Opm::data::Solution& sol, unsigned elemIdx, unsigned globalDofIndex)
{
@ -767,6 +941,165 @@ private:
return comm.rank() == 0;
}
void createLocalFipnum_()
{
const std::vector<int>& fipnum_global = simulator_.gridManager().eclState().get3DProperties().getIntGridProperty("FIPNUM").getData();
// Get compressed cell fipnum.
const auto& gridView = simulator_.gridManager().gridView();
unsigned numElements = gridView.size(/*codim=*/0);
fipnum_.resize(numElements);
if (fipnum_global.empty()) {
std::fill(fipnum_.begin(), fipnum_.end(), 0);
} else {
for (size_t elemIdx = 0; elemIdx < numElements; ++elemIdx) {
fipnum_[elemIdx] = fipnum_global[simulator_.gridManager().cartesianIndex( elemIdx )];
}
}
}
// Sum Fip values over regions.
ScalarBuffer FIPTotals_(const ScalarBuffer& fip, std::vector<int>& regionId, size_t maxNumberOfRegions, bool commSum = true)
{
ScalarBuffer totals(maxNumberOfRegions, 0.0);
assert(regionId.size() == fip.size());
for (size_t j = 0; j < regionId.size(); ++j) {
const int regionIdx = regionId[j] - 1;
// the cell is not attributed to any region. ignore it!
if (regionIdx < 0)
continue;
assert(regionIdx < static_cast<int>(maxNumberOfRegions));
totals[regionIdx] += fip[j];
}
if (commSum) {
const auto& comm = simulator_.gridView().comm();
for (size_t i = 0; i < maxNumberOfRegions; ++i)
totals[i] = comm.sum(totals[i]);
}
return totals;
}
// computes the hydrocarbon volume weighted averaged pressure
// of a region (reg)
// if reg == -1, the field average is computed.
Scalar FipPav_(int reg)
{
Scalar pPvSum = 0.0;
Scalar pvSum = 0.0;
Scalar pPvHydrocarbonSum = 0.0;
Scalar pvHydrocarbonSum = 0.0;
size_t numElem = oilPressure_.size();
for (size_t elem = 0; elem < numElem; ++elem) {
if(static_cast<int>(fipnum_[elem]) == reg || -1 == reg)
{
Scalar hydrocarbon = 0.0;
if (FluidSystem::phaseIsActive(oilPhaseIdx))
hydrocarbon += saturation_[oilPhaseIdx][elem];
if (FluidSystem::phaseIsActive(gasPhaseIdx))
hydrocarbon += saturation_[gasPhaseIdx][elem];
pPvSum += oilPressure_[elem] * fip_[FIPDataType::PV][elem];
pvSum += fip_[FIPDataType::PV][elem];
pPvHydrocarbonSum += pPvSum * hydrocarbon;
pvHydrocarbonSum += pvSum * hydrocarbon;
}
}
const auto& comm = simulator_.gridView().comm();
pPvSum = comm.sum(pPvSum);
pvSum = comm.sum(pvSum);
pPvHydrocarbonSum = comm.sum(pPvHydrocarbonSum);
pvHydrocarbonSum = comm.sum(pvHydrocarbonSum);
if (pvHydrocarbonSum > 1e-10)
return pPvHydrocarbonSum / pvHydrocarbonSum;
// return the porevolume weighted pressure if no hydrocarbon
return pPvSum / pvSum;;
}
void FIPUnitConvert_(const Opm::UnitSystem& units,
ScalarBuffer& fip)
{
if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_FIELD) {
fip[FIPDataType::WIP] = Opm::unit::convert::to(fip[FIPDataType::WIP], Opm::unit::stb);
fip[FIPDataType::OIP] = Opm::unit::convert::to(fip[FIPDataType::OIP], Opm::unit::stb);
fip[FIPDataType::OIPL] = Opm::unit::convert::to(fip[FIPDataType::OIPL], Opm::unit::stb);
fip[FIPDataType::OIPG] = Opm::unit::convert::to(fip[FIPDataType::OIPG], Opm::unit::stb);
fip[FIPDataType::GIP] = Opm::unit::convert::to(fip[FIPDataType::GIP], 1000*Opm::unit::cubic(Opm::unit::feet));
fip[FIPDataType::GIPL] = Opm::unit::convert::to(fip[FIPDataType::GIPL], 1000*Opm::unit::cubic(Opm::unit::feet));
fip[FIPDataType::GIPG] = Opm::unit::convert::to(fip[FIPDataType::GIPG], 1000*Opm::unit::cubic(Opm::unit::feet));
fip[FIPDataType::PV] = Opm::unit::convert::to(fip[FIPDataType::PV], Opm::unit::stb);
fip[FIPDataType::PAV] = Opm::unit::convert::to(fip[FIPDataType::PAV], Opm::unit::psia);
}
else if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_METRIC) {
fip[FIPDataType::PAV] = Opm::unit::convert::to(fip[FIPDataType::PAV], Opm::unit::barsa);
}
else {
OPM_THROW(std::runtime_error, "Unsupported unit type for fluid in place output.");
}
}
void outputRegionFluidInPlace_(const ScalarBuffer& oip, const ScalarBuffer& cip, const Opm::UnitSystem& units, const int reg)
{
std::ostringstream ss;
if (!reg) {
ss << " ===================================================\n"
<< " : Field Totals :\n";
} else {
ss << " ===================================================\n"
<< " : FIPNUM report region "
<< std::setw(2) << reg << " :\n";
}
if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_METRIC) {
ss << " : PAV =" << std::setw(14) << cip[FIPDataType::PAV] << " BARSA :\n"
<< std::fixed << std::setprecision(0)
<< " : PORV =" << std::setw(14) << cip[FIPDataType::PV] << " RM3 :\n";
if (!reg) {
ss << " : Pressure is weighted by hydrocarbon pore volume :\n"
<< " : Porv volumes are taken at reference conditions :\n";
}
ss << " :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n";
}
if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_FIELD) {
ss << " : PAV =" << std::setw(14) << cip[FIPDataType::PAV] << " PSIA :\n"
<< std::fixed << std::setprecision(0)
<< " : PORV =" << std::setw(14) << cip[FIPDataType::PV] << " RB :\n";
if (!reg) {
ss << " : Pressure is weighted by hydrocarbon pore volume :\n"
<< " : Pore volumes are taken at reference conditions :\n";
}
ss << " :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n";
}
ss << " : Liquid Vapour Total : Total : Free Dissolved Total :" << "\n"
<< ":------------------------:------------------------------------------:----------------:------------------------------------------:" << "\n"
<< ":Currently in place :" << std::setw(14) << cip[FIPDataType::OIPL] << std::setw(14) << cip[FIPDataType::OIPG] << std::setw(14) << cip[FIPDataType::OIP] << ":"
<< std::setw(13) << cip[FIPDataType::WIP] << " :" << std::setw(14) << (cip[FIPDataType::GIPG]) << std::setw(14) << cip[FIPDataType::GIPL] << std::setw(14) << cip[FIPDataType::GIP] << ":\n"
<< ":------------------------:------------------------------------------:----------------:------------------------------------------:\n"
<< ":Originally in place :" << std::setw(14) << oip[FIPDataType::OIPL] << std::setw(14) << oip[FIPDataType::OIPG] << std::setw(14) << oip[FIPDataType::OIP] << ":"
<< std::setw(13) << oip[FIPDataType::WIP] << " :" << std::setw(14) << oip[FIPDataType::GIPG] << std::setw(14) << oip[FIPDataType::GIPL] << std::setw(14) << oip[FIPDataType::GIP] << ":\n"
<< ":========================:==========================================:================:==========================================:\n";
Opm::OpmLog::note(ss.str());
}
std::string stringOfEnumIndex_(int i) {
typedef typename FIPDataType::FipId FipId;
switch( static_cast<FipId>(i) )
{
case FIPDataType::WIP: return "WIP"; break;
case FIPDataType::OIP: return "OIP"; break;
case FIPDataType::GIP: return "GIP"; break;
case FIPDataType::OIPL:return "OIPL"; break;
case FIPDataType::OIPG:return "OIPG"; break;
case FIPDataType::GIPL:return "GIPL"; break;
case FIPDataType::GIPG:return "GIPG"; break;
case FIPDataType::PV: return "PV"; break;
case FIPDataType::PAV: return "PAV"; break;
}
return "ERROR";
}
const Simulator& simulator_;
ScalarBuffer saturation_[numPhases];
@ -794,6 +1127,10 @@ private:
ScalarBuffer dewPointPressure_;
std::vector<int> failedCellsPb_;
std::vector<int> failedCellsPd_;
std::vector<int> fipnum_;
ScalarBuffer fip_[FIPDataType::fipValues];
ScalarBuffer origTotalValues_;
ScalarBuffer origRegionValues_[FIPDataType::fipValues];
};
} // namespace Ewoms

12
ebos/eclproblem.hh
View File

@ -460,6 +460,8 @@ public:
int numElements = gridView.size(/*codim=*/0);
maxPolymerAdsorption_.resize(numElements, 0.0);
}
eclWriter_->writeInit();
}
void prefetch(const Element& elem) const
@ -702,11 +704,10 @@ public:
}
Scalar totalSolverTime = 0.0;
Scalar nextstep = this->simulator().timeStepSize();
Opm::data::Solution fip;
writeOutput(dw, t, false, totalSolverTime, nextstep, fip, verbose);
writeOutput(dw, t, false, totalSolverTime, nextstep, verbose);
}
void writeOutput(const Opm::data::Wells& dw, Scalar t, bool substep, Scalar totalSolverTime, Scalar nextstep, const Opm::data::Solution& fip, bool verbose = true)
void writeOutput(const Opm::data::Wells& dw, Scalar t, bool substep, Scalar totalSolverTime, Scalar nextstep, bool verbose = true)
{
// use the generic code to prepare the output fields and to
// write the desired VTK files.
@ -714,7 +715,7 @@ public:
// output using eclWriter if enabled
if (eclWriter_)
eclWriter_->writeOutput(dw, t, substep, totalSolverTime, nextstep, fip);
eclWriter_->writeOutput(dw, t, substep, totalSolverTime, nextstep);
}
/*!
@ -1166,9 +1167,6 @@ public:
return initialFluidStates_[globalDofIdx];
}
void setEclIO(std::unique_ptr<Opm::EclipseIO>&& eclIO)
{ eclWriter_->setEclIO(std::move(eclIO)); }
const Opm::EclipseIO& eclIO() const
{ return eclWriter_->eclIO(); }

182
ebos/eclwriter.hh
View File

@ -95,10 +95,10 @@ public:
, eclOutputModule_(simulator)
, collectToIORank_(simulator_.gridManager())
{
Grid globalGrid = simulator_.gridManager().grid();
globalGrid.switchToGlobalView();
globalGrid_ = simulator_.gridManager().grid();
globalGrid_.switchToGlobalView();
eclIO_.reset(new Opm::EclipseIO(simulator_.gridManager().eclState(),
Opm::UgGridHelpers::createEclipseGrid( globalGrid , simulator_.gridManager().eclState().getInputGrid() ),
Opm::UgGridHelpers::createEclipseGrid( globalGrid_ , simulator_.gridManager().eclState().getInputGrid() ),
simulator_.gridManager().schedule(),
simulator_.gridManager().summaryConfig()));
}
@ -106,16 +106,28 @@ public:
~EclWriter()
{ }
void setEclIO(std::unique_ptr<Opm::EclipseIO>&& eclIO)
{ eclIO_ = std::move(eclIO); }
const Opm::EclipseIO& eclIO() const
{ return *eclIO_; }
void writeInit()
{
#if !HAVE_OPM_OUTPUT
OPM_THROW(std::runtime_error,
"Opm-output must be available to write ECL output!");
#else
if (collectToIORank_.isIORank()) {
std::map<std::string, std::vector<int> > integerVectors;
if (collectToIORank_.isParallel())
integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks());
eclIO_->writeInitial(computeTrans_(), integerVectors, exportNncStructure_());
}
#endif
}
/*!
* \brief collect and pass data and pass it to eclIO writer
*/
void writeOutput(const Opm::data::Wells& dw, Scalar t, bool substep, Scalar totalSolverTime, Scalar nextstep, const Opm::data::Solution& fip)
void writeOutput(const Opm::data::Wells& dw, Scalar t, bool substep, Scalar totalSolverTime, Scalar nextstep)
{
#if !HAVE_OPM_OUTPUT
OPM_THROW(std::runtime_error,
@ -133,6 +145,9 @@ public:
const ElementIterator& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity)
continue;
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
eclOutputModule_.processElement(elemCtx);
@ -140,11 +155,14 @@ public:
eclOutputModule_.outputErrorLog();
// collect all data to I/O rank and assign to sol
Opm::data::Solution localCellData = fip;
Opm::data::Solution localCellData;
eclOutputModule_.assignToSolution(localCellData);
if (collectToIORank_.isParallel())
collectToIORank_.collect(localCellData);
if (!substep)
eclOutputModule_.outputFIPLog();
// write output on I/O rank
if (collectToIORank_.isIORank()) {
@ -170,7 +188,7 @@ public:
false);
}
#endif
#endif
}
void restartBegin()
@ -194,7 +212,7 @@ public:
unsigned episodeIdx = simulator_.episodeIndex();
const auto& gridView = simulator_.gridManager().gridView();
unsigned numElements = gridView.size(/*codim=*/0);
eclOutputModule_.allocBuffers(numElements, episodeIdx, simulator_.gridManager().eclState().getRestartConfig(), true, false);
eclOutputModule_.allocBuffers(numElements, episodeIdx, simulator_.gridManager().eclState().getRestartConfig(), false, false);
auto restart_values = eclIO_->loadRestart(solution_keys, extra_keys);
for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
@ -213,10 +231,154 @@ private:
static bool enableEclOutput_()
{ return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); }
Opm::data::Solution computeTrans_() const
{
const auto& cartMapper = simulator_.gridManager().cartesianIndexMapper();
const auto& cartDims = cartMapper.cartesianDimensions();
const int globalSize = cartDims[0]*cartDims[1]*cartDims[2];
Opm::data::CellData tranx = {Opm::UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), Opm::data::TargetType::INIT};
Opm::data::CellData trany = {Opm::UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), Opm::data::TargetType::INIT};
Opm::data::CellData tranz = {Opm::UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), Opm::data::TargetType::INIT};
for (size_t i = 0; i < tranx.data.size(); ++i) {
tranx.data[0] = 0.0;
trany.data[0] = 0.0;
tranz.data[0] = 0.0;
}
const auto& globalGridView = globalGrid_.leafGridView();
typedef typename Grid::LeafGridView GridView;
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, Dune::MCMGElementLayout> ElementMapper;
ElementMapper globalElemMapper(globalGridView);
const auto& cartesianCellIdx = globalGrid_.globalCell();
const auto* globalTrans = &(simulator_.gridManager().globalTransmissibility());
if (!collectToIORank_.isParallel()) {
// in the sequential case we must use the transmissibilites defined by
// the problem. (because in the sequential case, the grid manager does
// not compute "global" transmissibilities for performance reasons. in
// the parallel case, the problem's transmissibilities can't be used
// because this object refers to the distributed grid and we need the
// sequential version here.)
globalTrans = &simulator_.problem().eclTransmissibilities();
}
auto elemIt = globalGridView.template begin</*codim=*/0>();
const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++ elemIt) {
const auto& elem = *elemIt;
auto isIt = globalGridView.ibegin(elem);
const auto& isEndIt = globalGridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
const auto& is = *isIt;
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.
}
int gc1 = std::min(cartesianCellIdx[c1], cartesianCellIdx[c2]);
int gc2 = std::max(cartesianCellIdx[c1], cartesianCellIdx[c2]);
if (gc2 - gc1 == 1) {
tranx.data[gc1] = globalTrans->transmissibility(c1, c2);
}
if (gc2 - gc1 == cartDims[0]) {
trany.data[gc1] = globalTrans->transmissibility(c1, c2);
}
if (gc2 - gc1 == cartDims[0]*cartDims[1]) {
tranz.data[gc1] = globalTrans->transmissibility(c1, c2);
}
}
}
return {{"TRANX" , tranx},
{"TRANY" , trany} ,
{"TRANZ" , tranz}};
}
Opm::NNC exportNncStructure_() const
{
Opm::NNC nnc = eclState().getInputNNC();
int nx = eclState().getInputGrid().getNX();
int ny = eclState().getInputGrid().getNY();
const auto& globalGridView = globalGrid_.leafGridView();
typedef typename Grid::LeafGridView GridView;
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, Dune::MCMGElementLayout> ElementMapper;
ElementMapper globalElemMapper(globalGridView);
const auto* globalTrans = &(simulator_.gridManager().globalTransmissibility());
if (!collectToIORank_.isParallel()) {
// in the sequential case we must use the transmissibilites defined by
// the problem. (because in the sequential case, the grid manager does
// not compute "global" transmissibilities for performance reasons. in
// the parallel case, the problem's transmissibilities can't be used
// because this object refers to the distributed grid and we need the
// sequential version here.)
globalTrans = &simulator_.problem().eclTransmissibilities();
}
auto elemIt = globalGridView.template begin</*codim=*/0>();
const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++ elemIt) {
const auto& elem = *elemIt;
auto isIt = globalGridView.ibegin(elem);
const auto& isEndIt = globalGridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
const auto& is = *isIt;
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.
}
// TODO (?): use the cartesian index mapper to make this code work
// with grids other than Dune::CpGrid. The problem is that we need
// the a mapper for the sequential grid, not for the distributed one.
int cc1 = globalGrid_.globalCell()[c1];
int cc2 = globalGrid_.globalCell()[c2];
if (std::abs(cc1 - cc2) != 1 &&
std::abs(cc1 - cc2) != nx &&
std::abs(cc1 - cc2) != nx*ny)
{
nnc.addNNC(cc1, cc2, globalTrans->transmissibility(c1, c2));
}
}
}
return nnc;
}
const Opm::EclipseState& eclState() const
{ return simulator_.gridManager().eclState(); }
const Simulator& simulator_;
EclOutputBlackOilModule<TypeTag> eclOutputModule_;
CollectDataToIORankType collectToIORank_;
std::unique_ptr<Opm::EclipseIO> eclIO_;
Grid globalGrid_;
};
} // namespace Ewoms