mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
eclwriter: split in typetag dependent and typetag-independent parts
This commit is contained in:
Arne Morten Kvarving
parent
80cdcbb4aa
commit
70ece6d25a
@ -29,6 +29,7 @@ list (APPEND MAIN_SOURCE_FILES
|
||||
ebos/eclgenericthresholdpressure.cc
|
||||
ebos/eclgenerictracermodel.cc
|
||||
ebos/eclgenericvanguard.cc
|
||||
ebos/eclgenericwriter.cc
|
||||
ebos/ecltransmissibility.cc
|
||||
opm/core/props/phaseUsageFromDeck.cpp
|
||||
opm/core/props/satfunc/RelpermDiagnostics.cpp
|
||||
|
||||
538
ebos/eclgenericwriter.cc
Normal file
538
ebos/eclgenericwriter.cc
Normal file
@ -0,0 +1,538 @@
|
||||
// -*- 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.
|
||||
*/
|
||||
|
||||
#include <config.h>
|
||||
#include <ebos/eclgenericwriter.hh>
|
||||
|
||||
#include <opm/grid/CpGrid.hpp>
|
||||
#include <opm/grid/cpgrid/GridHelpers.hpp>
|
||||
#include <opm/grid/polyhedralgrid.hh>
|
||||
#include <opm/grid/utility/cartesianToCompressed.hpp>
|
||||
|
||||
#include <opm/output/eclipse/EclipseIO.hpp>
|
||||
#include <opm/output/eclipse/RestartValue.hpp>
|
||||
#include <opm/output/eclipse/Summary.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Schedule/Action/State.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Schedule/SummaryState.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQState.hpp>
|
||||
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
|
||||
|
||||
#include <dune/grid/common/mcmgmapper.hh>
|
||||
|
||||
#if HAVE_DUNE_FEM
|
||||
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
|
||||
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
|
||||
#include <ebos/femcpgridcompat.hh>
|
||||
#endif
|
||||
|
||||
#if HAVE_MPI
|
||||
#include <mpi.h>
|
||||
#endif
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
struct EclWriteTasklet : public Opm::TaskletInterface
|
||||
{
|
||||
Opm::Action::State actionState_;
|
||||
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::SummaryState& summaryState,
|
||||
const Opm::UDQState& udqState,
|
||||
Opm::EclipseIO& eclIO,
|
||||
int reportStepNum,
|
||||
bool isSubStep,
|
||||
double secondsElapsed,
|
||||
Opm::RestartValue restartValue,
|
||||
bool writeDoublePrecision)
|
||||
: actionState_(actionState)
|
||||
, summaryState_(summaryState)
|
||||
, udqState_(udqState)
|
||||
, eclIO_(eclIO)
|
||||
, reportStepNum_(reportStepNum)
|
||||
, isSubStep_(isSubStep)
|
||||
, secondsElapsed_(secondsElapsed)
|
||||
, restartValue_(restartValue)
|
||||
, writeDoublePrecision_(writeDoublePrecision)
|
||||
{ }
|
||||
|
||||
// callback to eclIO serial writeTimeStep method
|
||||
void run()
|
||||
{
|
||||
eclIO_.writeTimeStep(actionState_,
|
||||
summaryState_,
|
||||
udqState_,
|
||||
reportStepNum_,
|
||||
isSubStep_,
|
||||
secondsElapsed_,
|
||||
restartValue_,
|
||||
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,
|
||||
const TransmissibilityType& globalTrans,
|
||||
bool enableAsyncOutput)
|
||||
: collectToIORank_(grid,
|
||||
equilGrid,
|
||||
gridView,
|
||||
cartMapper,
|
||||
equilCartMapper)
|
||||
, grid_(grid)
|
||||
, gridView_(gridView)
|
||||
, schedule_(schedule)
|
||||
, eclState_(eclState)
|
||||
, summaryConfig_(summaryConfig)
|
||||
, globalTrans_(globalTrans)
|
||||
, cartMapper_(cartMapper)
|
||||
, equilCartMapper_(equilCartMapper)
|
||||
, equilGrid_(equilGrid)
|
||||
{
|
||||
if (collectToIORank_.isIORank()) {
|
||||
eclIO_.reset(new EclipseIO(eclState_,
|
||||
UgGridHelpers::createEclipseGrid(*equilGrid, eclState_.getInputGrid()),
|
||||
schedule_,
|
||||
summaryConfig_));
|
||||
|
||||
const auto& wbp_calculators = eclIO_->summary().wbp_calculators( schedule.size() - 1 );
|
||||
wbp_index_list_ = wbp_calculators.index_list();
|
||||
}
|
||||
if (collectToIORank_.isParallel()) {
|
||||
const auto& comm = grid_.comm();
|
||||
unsigned long size = wbp_index_list_.size();
|
||||
comm.broadcast(&size, 1, collectToIORank_.ioRank);
|
||||
if (!collectToIORank_.isIORank())
|
||||
wbp_index_list_.resize( size );
|
||||
comm.broadcast(wbp_index_list_.data(), size, collectToIORank_.ioRank);
|
||||
}
|
||||
// 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;
|
||||
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()
|
||||
{
|
||||
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), integerVectors, exportNncStructure_(cartMap));
|
||||
}
|
||||
}
|
||||
|
||||
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
|
||||
{
|
||||
const auto& cartMapper = *equilCartMapper_;
|
||||
const auto& cartDims = cartMapper.cartesianDimensions();
|
||||
const int globalSize = cartDims[0]*cartDims[1]*cartDims[2];
|
||||
|
||||
data::CellData tranx = {UnitSystem::measure::transmissibility, std::vector<double>(globalSize), data::TargetType::INIT};
|
||||
data::CellData trany = {UnitSystem::measure::transmissibility, std::vector<double>(globalSize), data::TargetType::INIT};
|
||||
data::CellData tranz = {UnitSystem::measure::transmissibility, std::vector<double>(globalSize), 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;
|
||||
}
|
||||
|
||||
using GlobalGridView = typename EquilGrid::LeafGridView;
|
||||
const GlobalGridView& globalGridView = equilGrid_->leafGridView();
|
||||
using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView>;
|
||||
GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
|
||||
|
||||
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.
|
||||
|
||||
// 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);
|
||||
|
||||
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::size_t nx = eclState_.getInputGrid().getNX();
|
||||
std::size_t ny = eclState_.getInputGrid().getNY();
|
||||
auto nncData = eclState_.getInputNNC().input();
|
||||
const auto& unitSystem = eclState_.getDeckUnitSystem();
|
||||
std::vector<NNCdata> outputNnc;
|
||||
std::size_t index = 0;
|
||||
|
||||
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);
|
||||
}
|
||||
++index;
|
||||
}
|
||||
|
||||
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();
|
||||
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.
|
||||
|
||||
std::size_t cc1 = equilCartMapper.cartesianIndex( c1 ); //globalIOGrid_.globalCell()[c1];
|
||||
std::size_t cc2 = equilCartMapper.cartesianIndex( c2 ); //globalIOGrid_.globalCell()[c2];
|
||||
|
||||
if ( cc2 < cc1 )
|
||||
std::swap(cc1, cc2);
|
||||
|
||||
auto cellDiff = cc2 - cc1;
|
||||
|
||||
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,
|
||||
const Action::State& actionState,
|
||||
const UDQState& udqState,
|
||||
const SummaryState& summaryState,
|
||||
const std::vector<Scalar>& thresholdPressure,
|
||||
Scalar curTime,
|
||||
Scalar nextStepSize,
|
||||
bool doublePrecision)
|
||||
{
|
||||
const auto isParallel = this->collectToIORank_.isParallel();
|
||||
|
||||
RestartValue restartValue {
|
||||
isParallel ? this->collectToIORank_.globalCellData()
|
||||
: std::move(localCellData),
|
||||
|
||||
isParallel ? this->collectToIORank_.globalWellData()
|
||||
: std::move(localWellData),
|
||||
|
||||
isParallel ? this->collectToIORank_.globalGroupAndNetworkData()
|
||||
: std::move(localGroupAndNetworkData)
|
||||
};
|
||||
|
||||
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));
|
||||
}
|
||||
|
||||
// first, create a tasklet to write the data for the current time
|
||||
// step to disk
|
||||
auto eclWriteTasklet = std::make_shared<EclWriteTasklet>(
|
||||
actionState, 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(int reportStepNum,
|
||||
Scalar curTime,
|
||||
const std::map<std::size_t, double>& wbpData,
|
||||
const data::Wells& localWellData,
|
||||
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,
|
||||
SummaryState& summaryState,
|
||||
UDQState& udqState,
|
||||
const Inplace& inplace,
|
||||
const Inplace& initialInPlace)
|
||||
{
|
||||
std::vector<char> buffer;
|
||||
if (collectToIORank_.isIORank()) {
|
||||
const auto& summary = eclIO_->summary();
|
||||
auto wbp_calculators = summary.wbp_calculators(reportStepNum);
|
||||
|
||||
for (const auto& [global_index, pressure] : wbpData)
|
||||
wbp_calculators.add_pressure( global_index, pressure );
|
||||
|
||||
const auto& wellData = this->collectToIORank_.isParallel()
|
||||
? this->collectToIORank_.globalWellData()
|
||||
: localWellData;
|
||||
|
||||
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,
|
||||
groupAndNetworkData,
|
||||
miscSummaryData,
|
||||
initialInPlace,
|
||||
inplace,
|
||||
wbp_calculators,
|
||||
regionData,
|
||||
blockData,
|
||||
aquiferData);
|
||||
|
||||
/*
|
||||
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.
|
||||
*/
|
||||
auto udq_step = reportStepNum - 1;
|
||||
const auto& udq_config = schedule_.getUDQConfig(udq_step);
|
||||
udq_config.eval( udq_step, schedule_.wellMatcher(udq_step), summaryState, udqState);
|
||||
|
||||
buffer = summaryState.serialize();
|
||||
}
|
||||
if (collectToIORank_.isParallel()) {
|
||||
#ifdef HAVE_MPI
|
||||
unsigned long buffer_size = buffer.size();
|
||||
MPI_Bcast(&buffer_size, 1, MPI_UNSIGNED_LONG, collectToIORank_.ioRank, MPI_COMM_WORLD);
|
||||
if (!collectToIORank_.isIORank())
|
||||
buffer.resize( buffer_size );
|
||||
|
||||
MPI_Bcast(buffer.data(), buffer_size, MPI_CHAR, collectToIORank_.ioRank, MPI_COMM_WORLD);
|
||||
if (!collectToIORank_.isIORank()) {
|
||||
SummaryState& st = summaryState;
|
||||
st.deserialize(buffer);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
#if HAVE_DUNE_FEM
|
||||
template class EclGenericWriter<Dune::CpGrid,
|
||||
Dune::CpGrid,
|
||||
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
|
||||
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>>,
|
||||
double>;
|
||||
#else
|
||||
template class EclGenericWriter<Dune::CpGrid,
|
||||
Dune::CpGrid,
|
||||
Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
|
||||
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>, Dune::Impl::MCMGFailLayout>,
|
||||
double>;
|
||||
#endif
|
||||
|
||||
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)>>, Dune::Impl::MCMGFailLayout>,
|
||||
double>;
|
||||
|
||||
} // namespace Opm
|
||||
129
ebos/eclgenericwriter.hh
Normal file
129
ebos/eclgenericwriter.hh
Normal file
@ -0,0 +1,129 @@
|
||||
// -*- 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.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::EclWriter
|
||||
*/
|
||||
#ifndef EWOMS_ECL_GENERIC_WRITER_HH
|
||||
#define EWOMS_ECL_GENERIC_WRITER_HH
|
||||
|
||||
#include "collecttoiorank.hh"
|
||||
#include <ebos/ecltransmissibility.hh>
|
||||
|
||||
#include <opm/models/parallel/tasklets.hh>
|
||||
|
||||
#include <cstddef>
|
||||
#include <map>
|
||||
#include <string>
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
namespace Action { class State; }
|
||||
class EclipseIO;
|
||||
class EclipseState;
|
||||
class Inplace;
|
||||
struct NNCdata;
|
||||
class Schedule;
|
||||
class SummaryConfig;
|
||||
class SummaryState;
|
||||
class UDQState;
|
||||
|
||||
template <class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
|
||||
class EclGenericWriter
|
||||
{
|
||||
using CollectDataToIORankType = CollectDataToIORank<Grid,EquilGrid,GridView>;
|
||||
using TransmissibilityType = EclTransmissibility<Grid,GridView,ElementMapper,Scalar>;
|
||||
|
||||
public:
|
||||
// The Simulator object should preferably have been const - the
|
||||
// only reason that is not the case is due to the SummaryState
|
||||
// object owned deep down by the vanguard.
|
||||
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,
|
||||
const TransmissibilityType& globalTrans,
|
||||
bool enableAsyncOutput);
|
||||
|
||||
const EclipseIO& eclIO() const;
|
||||
|
||||
void writeInit();
|
||||
|
||||
protected:
|
||||
void doWriteOutput(const int reportStepNum,
|
||||
const bool isSubStep,
|
||||
data::Solution&& localCellData,
|
||||
data::Wells&& localWellData,
|
||||
data::GroupAndNetworkValues&& localGroupAndNetworkData,
|
||||
const Action::State& actionState,
|
||||
const UDQState& udqState,
|
||||
const SummaryState& summaryState,
|
||||
const std::vector<Scalar>& thresholdPressure,
|
||||
Scalar curTime,
|
||||
Scalar nextStepSize,
|
||||
bool doublePrecision);
|
||||
|
||||
void evalSummary(int reportStepNum,
|
||||
Scalar curTime,
|
||||
const std::map<std::size_t, double>& wbpData,
|
||||
const data::Wells& localWellData,
|
||||
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,
|
||||
SummaryState& summaryState,
|
||||
UDQState& udqState,
|
||||
const Inplace& inplace,
|
||||
const Inplace& initialInPlace);
|
||||
|
||||
CollectDataToIORankType collectToIORank_;
|
||||
const Grid& grid_;
|
||||
const GridView& gridView_;
|
||||
const Schedule& schedule_;
|
||||
const EclipseState& eclState_;
|
||||
const SummaryConfig& summaryConfig_;
|
||||
std::unique_ptr<EclipseIO> eclIO_;
|
||||
std::unique_ptr<TaskletRunner> taskletRunner_;
|
||||
Scalar restartTimeStepSize_;
|
||||
const TransmissibilityType& globalTrans_;
|
||||
const Dune::CartesianIndexMapper<Grid>& cartMapper_;
|
||||
const Dune::CartesianIndexMapper<EquilGrid>* equilCartMapper_;
|
||||
const EquilGrid* equilGrid_;
|
||||
std::vector<std::size_t> wbp_index_list_;
|
||||
|
||||
private:
|
||||
data::Solution computeTrans_(const std::unordered_map<int,int>& cartesianToActive) const;
|
||||
std::vector<NNCdata> exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive) const;
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
@ -95,6 +95,7 @@
|
||||
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Schedule/Action/ActionContext.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Schedule/Action/ActionX.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Schedule/Action/State.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Tables/RockwnodTable.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Tables/OverburdTable.hpp>
|
||||
@ -838,7 +839,7 @@ public:
|
||||
ExtboModule::initFromState(vanguard.eclState());
|
||||
|
||||
// create the ECL writer
|
||||
eclWriter_.reset(new EclWriterType(simulator));
|
||||
eclWriter_.reset(new EclWriterType(simulator, *this));
|
||||
|
||||
enableDriftCompensation_ = EWOMS_GET_PARAM(TypeTag, bool, EclEnableDriftCompensation);
|
||||
|
||||
|
||||
@ -31,39 +31,13 @@
|
||||
#include "collecttoiorank.hh"
|
||||
#include "ecloutputblackoilmodule.hh"
|
||||
|
||||
#include <opm/models/blackoil/blackoilmodel.hh>
|
||||
|
||||
#include <opm/simulators/wells/BlackoilWellModel.hpp>
|
||||
|
||||
#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
|
||||
#include <opm/models/io/baseoutputwriter.hh>
|
||||
#include <opm/models/parallel/tasklets.hh>
|
||||
|
||||
#include <opm/output/eclipse/EclipseIO.hpp>
|
||||
|
||||
#include <opm/output/eclipse/RestartValue.hpp>
|
||||
#include <opm/output/eclipse/Summary.hpp>
|
||||
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Schedule/Action/State.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQConfig.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQState.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Schedule/Well/PAvgCalculatorCollection.hpp>
|
||||
|
||||
#include <opm/simulators/utils/ParallelRestart.hpp>
|
||||
#include <opm/grid/GridHelpers.hpp>
|
||||
#include <opm/grid/utility/cartesianToCompressed.hpp>
|
||||
|
||||
#include <opm/material/common/Valgrind.hpp>
|
||||
#include <opm/common/OpmLog/OpmLog.hpp>
|
||||
#include <ebos/eclgenericwriter.hh>
|
||||
|
||||
#include <list>
|
||||
#include <utility>
|
||||
#include <string>
|
||||
#include <chrono>
|
||||
|
||||
#ifdef HAVE_MPI
|
||||
#include <mpi.h>
|
||||
#endif
|
||||
|
||||
namespace Opm::Properties {
|
||||
|
||||
@ -84,60 +58,9 @@ struct EclOutputDoublePrecision {
|
||||
|
||||
namespace Opm {
|
||||
|
||||
template <class TypeTag>
|
||||
class EclWriter;
|
||||
|
||||
template <class TypeTag>
|
||||
class EclOutputBlackOilModule;
|
||||
|
||||
/*!
|
||||
* \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.
|
||||
*/
|
||||
inline
|
||||
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;
|
||||
}
|
||||
namespace Action { class State; }
|
||||
class EclipseIO;
|
||||
class UDQState;
|
||||
|
||||
/*!
|
||||
* \ingroup EclBlackOilSimulator
|
||||
@ -155,7 +78,11 @@ bool directVerticalNeighbors(const std::array<int, 3>& cartDims,
|
||||
* centered finite volume discretization.
|
||||
*/
|
||||
template <class TypeTag>
|
||||
class EclWriter
|
||||
class EclWriter : public EclGenericWriter<GetPropType<TypeTag, Properties::Grid>,
|
||||
GetPropType<TypeTag, Properties::EquilGrid>,
|
||||
GetPropType<TypeTag, Properties::GridView>,
|
||||
GetPropType<TypeTag, Properties::ElementMapper>,
|
||||
GetPropType<TypeTag, Properties::Scalar>>
|
||||
{
|
||||
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
|
||||
using Vanguard = GetPropType<TypeTag, Properties::Vanguard>;
|
||||
@ -166,16 +93,13 @@ class EclWriter
|
||||
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
|
||||
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
|
||||
using Element = typename GridView::template Codim<0>::Entity;
|
||||
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
|
||||
using ElementIterator = typename GridView::template Codim<0>::Iterator;
|
||||
|
||||
using CollectDataToIORankType = CollectDataToIORank<Grid,EquilGrid,GridView>;
|
||||
|
||||
typedef std::vector<Scalar> ScalarBuffer;
|
||||
using BaseType = EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>;
|
||||
|
||||
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
|
||||
enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
|
||||
|
||||
|
||||
public:
|
||||
static void registerParameters()
|
||||
{
|
||||
@ -188,73 +112,32 @@ public:
|
||||
// The Simulator object should preferably have been const - the
|
||||
// only reason that is not the case is due to the SummaryState
|
||||
// object owned deep down by the vanguard.
|
||||
EclWriter(Simulator& simulator)
|
||||
: simulator_(simulator)
|
||||
, collectToIORank_(simulator_.vanguard().grid(),
|
||||
simulator_.vanguard().grid().comm().rank() == 0 ?
|
||||
&simulator_.vanguard().equilGrid() : nullptr,
|
||||
simulator_.vanguard().gridView(),
|
||||
simulator_.vanguard().cartesianIndexMapper(),
|
||||
simulator_.vanguard().grid().comm().rank() == 0 ?
|
||||
&simulator_.vanguard().equilCartesianIndexMapper() : nullptr)
|
||||
|
||||
template<class Problem>
|
||||
EclWriter(Simulator& simulator, const Problem& problem)
|
||||
: BaseType(simulator.vanguard().schedule(),
|
||||
simulator.vanguard().eclState(),
|
||||
simulator.vanguard().summaryConfig(),
|
||||
simulator.vanguard().grid(),
|
||||
simulator.vanguard().grid().comm().rank() == 0 ? &simulator.vanguard().equilGrid() : nullptr,
|
||||
simulator.vanguard().gridView(),
|
||||
simulator.vanguard().cartesianIndexMapper(),
|
||||
simulator.vanguard().grid().comm().rank() == 0 ? &simulator.vanguard().equilCartesianIndexMapper() : nullptr,
|
||||
simulator.vanguard().grid().comm().size() > 1 ? simulator.vanguard().globalTransmissibility() : problem.eclTransmissibilities(),
|
||||
EWOMS_GET_PARAM(TypeTag, bool, EnableAsyncEclOutput))
|
||||
, simulator_(simulator)
|
||||
{
|
||||
std::vector<std::size_t> wbp_index_list;
|
||||
if (collectToIORank_.isIORank()) {
|
||||
const auto& schedule = simulator_.vanguard().schedule();
|
||||
eclIO_.reset(new EclipseIO(simulator_.vanguard().eclState(),
|
||||
UgGridHelpers::createEclipseGrid(globalGrid(), simulator_.vanguard().eclState().getInputGrid()),
|
||||
schedule,
|
||||
simulator_.vanguard().summaryConfig()));
|
||||
|
||||
const auto& wbp_calculators = eclIO_->summary().wbp_calculators( schedule.size() - 1 );
|
||||
wbp_index_list = wbp_calculators.index_list();
|
||||
}
|
||||
if (collectToIORank_.isParallel()) {
|
||||
const auto& comm = simulator_.vanguard().grid().comm();
|
||||
unsigned long size = wbp_index_list.size();
|
||||
comm.broadcast(&size, 1, collectToIORank_.ioRank);
|
||||
if (!collectToIORank_.isIORank())
|
||||
wbp_index_list.resize( size );
|
||||
comm.broadcast(wbp_index_list.data(), size, collectToIORank_.ioRank);
|
||||
}
|
||||
// create output thread if enabled and rank is I/O rank
|
||||
// async output is enabled by default if pthread are enabled
|
||||
bool enableAsyncOutput = EWOMS_GET_PARAM(TypeTag, bool, EnableAsyncEclOutput);
|
||||
int numWorkerThreads = 0;
|
||||
if (enableAsyncOutput && collectToIORank_.isIORank())
|
||||
numWorkerThreads = 1;
|
||||
taskletRunner_.reset(new TaskletRunner(numWorkerThreads));
|
||||
|
||||
this->eclOutputModule_ = std::make_unique<EclOutputBlackOilModule<TypeTag>>(simulator, wbp_index_list, this->collectToIORank_);
|
||||
this->eclOutputModule_ = std::make_unique<EclOutputBlackOilModule<TypeTag>>(simulator, this->wbp_index_list_, this->collectToIORank_);
|
||||
this->wbp_index_list_.clear();
|
||||
}
|
||||
|
||||
~EclWriter()
|
||||
{ }
|
||||
|
||||
const EclipseIO& eclIO() const
|
||||
{
|
||||
assert(eclIO_);
|
||||
return *eclIO_;
|
||||
}
|
||||
|
||||
const EquilGrid& globalGrid() const
|
||||
{
|
||||
return simulator_.vanguard().equilGrid();
|
||||
}
|
||||
|
||||
void writeInit()
|
||||
{
|
||||
if (collectToIORank_.isIORank()) {
|
||||
std::map<std::string, std::vector<int> > integerVectors;
|
||||
if (collectToIORank_.isParallel())
|
||||
integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks());
|
||||
auto cartMap = cartesianToCompressed(globalGrid().size(0),
|
||||
UgGridHelpers::globalCell(globalGrid()));
|
||||
eclIO_->writeInitial(computeTrans_(cartMap), integerVectors, exportNncStructure_(cartMap));
|
||||
}
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief collect and pass data and pass it to eclIO writer
|
||||
*/
|
||||
@ -297,13 +180,13 @@ public:
|
||||
|
||||
this->prepareLocalCellData(isSubStep, reportStepNum);
|
||||
|
||||
if (collectToIORank_.isParallel())
|
||||
collectToIORank_.collect({},
|
||||
eclOutputModule_->getBlockData(),
|
||||
eclOutputModule_->getWBPData(),
|
||||
localWellData,
|
||||
localGroupAndNetworkData,
|
||||
localAquiferData);
|
||||
if (this->collectToIORank_.isParallel())
|
||||
this->collectToIORank_.collect({},
|
||||
eclOutputModule_->getBlockData(),
|
||||
eclOutputModule_->getWBPData(),
|
||||
localWellData,
|
||||
localGroupAndNetworkData,
|
||||
localAquiferData);
|
||||
|
||||
|
||||
std::map<std::string, double> miscSummaryData;
|
||||
@ -317,81 +200,25 @@ public:
|
||||
eclOutputModule_->outputInjLog(reportStepNum, isSubStep, forceDisableInjOutput);
|
||||
eclOutputModule_->outputCumLog(reportStepNum, isSubStep, forceDisableCumOutput);
|
||||
|
||||
|
||||
std::vector<char> buffer;
|
||||
if (this->collectToIORank_.isIORank()) {
|
||||
const auto& summary = eclIO_->summary();
|
||||
auto wbp_calculators = summary.wbp_calculators(reportStepNum);
|
||||
const auto& wbpData
|
||||
= this->collectToIORank_.isParallel()
|
||||
? this->collectToIORank_.globalWBPData()
|
||||
: this->eclOutputModule_->getWBPData();
|
||||
|
||||
for (const auto& [global_index, pressure] : wbpData)
|
||||
wbp_calculators.add_pressure( global_index, pressure );
|
||||
|
||||
// Add TCPU
|
||||
if (totalCpuTime != 0.0) {
|
||||
miscSummaryData["TCPU"] = totalCpuTime;
|
||||
}
|
||||
|
||||
const auto& wellData = this->collectToIORank_.isParallel()
|
||||
? this->collectToIORank_.globalWellData()
|
||||
: localWellData;
|
||||
|
||||
const auto& groupAndNetworkData = this->collectToIORank_.isParallel()
|
||||
? this->collectToIORank_.globalGroupAndNetworkData()
|
||||
: localGroupAndNetworkData;
|
||||
|
||||
const auto& aquiferData = this->collectToIORank_.isParallel()
|
||||
? this->collectToIORank_.globalAquiferData()
|
||||
: localAquiferData;
|
||||
|
||||
const auto& blockData
|
||||
= this->collectToIORank_.isParallel()
|
||||
? this->collectToIORank_.globalBlockData()
|
||||
: this->eclOutputModule_->getBlockData();
|
||||
|
||||
summary.eval(summaryState(),
|
||||
reportStepNum,
|
||||
curTime,
|
||||
wellData,
|
||||
groupAndNetworkData,
|
||||
miscSummaryData,
|
||||
eclOutputModule_->initialInplace(),
|
||||
inplace,
|
||||
wbp_calculators,
|
||||
regionData,
|
||||
blockData,
|
||||
aquiferData);
|
||||
|
||||
/*
|
||||
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.
|
||||
*/
|
||||
auto udq_step = reportStepNum - 1;
|
||||
const auto& udq_config = schedule().getUDQConfig(udq_step);
|
||||
udq_config.eval( udq_step, schedule().wellMatcher(udq_step), summaryState(), udqState() );
|
||||
|
||||
buffer = summaryState().serialize();
|
||||
// Add TCPU
|
||||
if (totalCpuTime != 0.0) {
|
||||
miscSummaryData["TCPU"] = totalCpuTime;
|
||||
}
|
||||
|
||||
if (collectToIORank_.isParallel()) {
|
||||
#ifdef HAVE_MPI
|
||||
unsigned long buffer_size = buffer.size();
|
||||
MPI_Bcast(&buffer_size, 1, MPI_UNSIGNED_LONG, collectToIORank_.ioRank, MPI_COMM_WORLD);
|
||||
if (!collectToIORank_.isIORank())
|
||||
buffer.resize( buffer_size );
|
||||
|
||||
MPI_Bcast(buffer.data(), buffer_size, MPI_CHAR, collectToIORank_.ioRank, MPI_COMM_WORLD);
|
||||
if (!collectToIORank_.isIORank()) {
|
||||
SummaryState& st = summaryState();
|
||||
st.deserialize(buffer);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
this->evalSummary(reportStepNum, curTime,
|
||||
this->collectToIORank_.isParallel() ?
|
||||
this->collectToIORank_.globalWBPData() :
|
||||
this->eclOutputModule_->getWBPData(),
|
||||
localWellData,
|
||||
localGroupAndNetworkData,
|
||||
localAquiferData,
|
||||
this->collectToIORank_.isParallel() ?
|
||||
this->collectToIORank_.globalBlockData() :
|
||||
this->eclOutputModule_->getBlockData(),
|
||||
miscSummaryData, regionData,
|
||||
summaryState(), udqState(),
|
||||
inplace,
|
||||
eclOutputModule_->initialInplace());
|
||||
}
|
||||
|
||||
|
||||
@ -416,19 +243,27 @@ public:
|
||||
}
|
||||
|
||||
if (this->collectToIORank_.isParallel()) {
|
||||
collectToIORank_.collect(localCellData,
|
||||
eclOutputModule_->getBlockData(),
|
||||
eclOutputModule_->getWBPData(),
|
||||
localWellData,
|
||||
localGroupAndNetworkData,
|
||||
{});
|
||||
this->collectToIORank_.collect(localCellData,
|
||||
eclOutputModule_->getBlockData(),
|
||||
eclOutputModule_->getWBPData(),
|
||||
localWellData,
|
||||
localGroupAndNetworkData,
|
||||
{});
|
||||
}
|
||||
|
||||
if (this->collectToIORank_.isIORank()) {
|
||||
this->writeOutput(reportStepNum, isSubStep,
|
||||
std::move(localCellData),
|
||||
std::move(localWellData),
|
||||
std::move(localGroupAndNetworkData));
|
||||
const Scalar curTime = simulator_.time() + simulator_.timeStepSize();
|
||||
const Scalar nextStepSize = simulator_.problem().nextTimeStepSize();
|
||||
this->doWriteOutput(reportStepNum, isSubStep,
|
||||
std::move(localCellData),
|
||||
std::move(localWellData),
|
||||
std::move(localGroupAndNetworkData),
|
||||
this->actionState(),
|
||||
this->udqState(),
|
||||
this->summaryState(),
|
||||
simulator_.problem().thresholdPressure().data(),
|
||||
curTime, nextStepSize,
|
||||
EWOMS_GET_PARAM(TypeTag, bool, EclOutputDoublePrecision));
|
||||
}
|
||||
}
|
||||
|
||||
@ -470,10 +305,10 @@ public:
|
||||
{
|
||||
SummaryState& summaryState = simulator_.vanguard().summaryState();
|
||||
Action::State& actionState = simulator_.vanguard().actionState();
|
||||
auto restartValues = loadParallelRestart(eclIO_.get(), actionState, summaryState, solutionKeys, extraKeys,
|
||||
auto restartValues = loadParallelRestart(this->eclIO_.get(), actionState, summaryState, solutionKeys, extraKeys,
|
||||
gridView.grid().comm());
|
||||
for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
|
||||
unsigned globalIdx = collectToIORank_.localIdxToGlobalIdx(elemIdx);
|
||||
unsigned globalIdx = this->collectToIORank_.localIdxToGlobalIdx(elemIdx);
|
||||
eclOutputModule_->setRestart(restartValues.solution, elemIdx, globalIdx);
|
||||
}
|
||||
|
||||
@ -502,247 +337,10 @@ public:
|
||||
Scalar restartTimeStepSize() const
|
||||
{ return restartTimeStepSize_; }
|
||||
|
||||
|
||||
private:
|
||||
static bool enableEclOutput_()
|
||||
{ return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); }
|
||||
|
||||
data::Solution computeTrans_(const std::unordered_map<int,int>& cartesianToActive) const
|
||||
{
|
||||
const auto& cartMapper = simulator_.vanguard().equilCartesianIndexMapper();
|
||||
const auto& cartDims = cartMapper.cartesianDimensions();
|
||||
const int globalSize = cartDims[0]*cartDims[1]*cartDims[2];
|
||||
|
||||
data::CellData tranx = {UnitSystem::measure::transmissibility, std::vector<double>(globalSize), data::TargetType::INIT};
|
||||
data::CellData trany = {UnitSystem::measure::transmissibility, std::vector<double>(globalSize), data::TargetType::INIT};
|
||||
data::CellData tranz = {UnitSystem::measure::transmissibility, std::vector<double>(globalSize), 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;
|
||||
}
|
||||
|
||||
typedef typename EquilGrid :: LeafGridView GlobalGridView;
|
||||
const GlobalGridView& globalGridView = globalGrid().leafGridView();
|
||||
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView> ElementMapper;
|
||||
ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
|
||||
|
||||
using TransmissibilityType = typename Vanguard::TransmissibilityType;
|
||||
|
||||
const TransmissibilityType* globalTrans;
|
||||
|
||||
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();
|
||||
}
|
||||
else
|
||||
{
|
||||
globalTrans = &(simulator_.vanguard().globalTransmissibility());
|
||||
}
|
||||
|
||||
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.
|
||||
|
||||
// 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);
|
||||
|
||||
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}};
|
||||
}
|
||||
|
||||
std::vector<NNCdata> exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive) const
|
||||
{
|
||||
std::size_t nx = eclState().getInputGrid().getNX();
|
||||
std::size_t ny = eclState().getInputGrid().getNY();
|
||||
auto nncData = eclState().getInputNNC().input();
|
||||
const auto& unitSystem = simulator_.vanguard().eclState().getDeckUnitSystem();
|
||||
std::vector<NNCdata> outputNnc;
|
||||
std::size_t index = 0;
|
||||
|
||||
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);
|
||||
}
|
||||
++index;
|
||||
}
|
||||
|
||||
typedef typename EquilGrid :: LeafGridView GlobalGridView;
|
||||
const GlobalGridView& globalGridView = globalGrid().leafGridView();
|
||||
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView> ElementMapper;
|
||||
ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
|
||||
|
||||
using TransmissibilityType = typename Vanguard::TransmissibilityType;
|
||||
const TransmissibilityType* globalTrans;
|
||||
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();
|
||||
}
|
||||
else
|
||||
{
|
||||
globalTrans = &(simulator_.vanguard().globalTransmissibility());
|
||||
}
|
||||
|
||||
// Cartesian index mapper for the serial I/O grid
|
||||
const auto& equilCartMapper = simulator_.vanguard().equilCartesianIndexMapper();
|
||||
const auto& cartDims = simulator_.vanguard().cartesianIndexMapper().cartesianDimensions();
|
||||
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.
|
||||
|
||||
std::size_t cc1 = equilCartMapper.cartesianIndex( c1 ); //globalIOGrid_.globalCell()[c1];
|
||||
std::size_t cc2 = equilCartMapper.cartesianIndex( c2 ); //globalIOGrid_.globalCell()[c2];
|
||||
|
||||
if ( cc2 < cc1 )
|
||||
std::swap(cc1, cc2);
|
||||
|
||||
auto cellDiff = cc2 - cc1;
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
struct EclWriteTasklet
|
||||
: public TaskletInterface
|
||||
{
|
||||
Action::State actionState_;
|
||||
SummaryState summaryState_;
|
||||
UDQState udqState_;
|
||||
EclipseIO& eclIO_;
|
||||
int reportStepNum_;
|
||||
bool isSubStep_;
|
||||
double secondsElapsed_;
|
||||
RestartValue restartValue_;
|
||||
bool writeDoublePrecision_;
|
||||
|
||||
explicit EclWriteTasklet(const Action::State& actionState,
|
||||
const SummaryState& summaryState,
|
||||
const UDQState& udqState,
|
||||
EclipseIO& eclIO,
|
||||
int reportStepNum,
|
||||
bool isSubStep,
|
||||
double secondsElapsed,
|
||||
RestartValue restartValue,
|
||||
bool writeDoublePrecision)
|
||||
: actionState_(actionState)
|
||||
, summaryState_(summaryState)
|
||||
, udqState_(udqState)
|
||||
, eclIO_(eclIO)
|
||||
, reportStepNum_(reportStepNum)
|
||||
, isSubStep_(isSubStep)
|
||||
, secondsElapsed_(secondsElapsed)
|
||||
, restartValue_(restartValue)
|
||||
, writeDoublePrecision_(writeDoublePrecision)
|
||||
{ }
|
||||
|
||||
// callback to eclIO serial writeTimeStep method
|
||||
void run()
|
||||
{
|
||||
eclIO_.writeTimeStep(actionState_,
|
||||
summaryState_,
|
||||
udqState_,
|
||||
reportStepNum_,
|
||||
isSubStep_,
|
||||
secondsElapsed_,
|
||||
restartValue_,
|
||||
writeDoublePrecision_);
|
||||
}
|
||||
};
|
||||
|
||||
const EclipseState& eclState() const
|
||||
{ return simulator_.vanguard().eclState(); }
|
||||
|
||||
@ -763,7 +361,7 @@ private:
|
||||
{
|
||||
const auto& gridView = simulator_.vanguard().gridView();
|
||||
const int numElements = gridView.size(/*codim=*/0);
|
||||
const bool log = collectToIORank_.isIORank();
|
||||
const bool log = this->collectToIORank_.isIORank();
|
||||
|
||||
eclOutputModule_->allocBuffers(numElements, reportStepNum,
|
||||
isSubStep, log, /*isRestart*/ false);
|
||||
@ -782,59 +380,8 @@ private:
|
||||
}
|
||||
}
|
||||
|
||||
void writeOutput(const int reportStepNum,
|
||||
const bool isSubStep,
|
||||
::Opm::data::Solution&& localCellData,
|
||||
::Opm::data::Wells&& localWellData,
|
||||
::Opm::data::GroupAndNetworkValues&& localGroupAndNetworkData)
|
||||
{
|
||||
const Scalar curTime = simulator_.time() + simulator_.timeStepSize();
|
||||
const Scalar nextStepSize = simulator_.problem().nextTimeStepSize();
|
||||
const auto isParallel = this->collectToIORank_.isParallel();
|
||||
|
||||
RestartValue restartValue {
|
||||
isParallel ? this->collectToIORank_.globalCellData()
|
||||
: std::move(localCellData),
|
||||
|
||||
isParallel ? this->collectToIORank_.globalWellData()
|
||||
: std::move(localWellData),
|
||||
|
||||
isParallel ? this->collectToIORank_.globalGroupAndNetworkData()
|
||||
: std::move(localGroupAndNetworkData)
|
||||
};
|
||||
|
||||
if (simulator_.vanguard().eclState().getSimulationConfig().useThresholdPressure()) {
|
||||
restartValue.addExtra("THRESHPR", UnitSystem::measure::pressure,
|
||||
simulator_.problem().thresholdPressure().data());
|
||||
}
|
||||
|
||||
// Add suggested next timestep to extra data.
|
||||
if (! isSubStep) {
|
||||
restartValue.addExtra("OPMEXTRA", std::vector<double>(1, nextStepSize));
|
||||
}
|
||||
|
||||
// first, create a tasklet to write the data for the current time
|
||||
// step to disk
|
||||
auto eclWriteTasklet = std::make_shared<EclWriteTasklet>(
|
||||
this->actionState(), this->summaryState(), this->udqState(), *this->eclIO_,
|
||||
reportStepNum, isSubStep, curTime, std::move(restartValue),
|
||||
EWOMS_GET_PARAM(TypeTag, bool, EclOutputDoublePrecision)
|
||||
);
|
||||
|
||||
// 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));
|
||||
}
|
||||
|
||||
Simulator& simulator_;
|
||||
CollectDataToIORankType collectToIORank_;
|
||||
std::unique_ptr<EclOutputBlackOilModule<TypeTag>> eclOutputModule_;
|
||||
std::unique_ptr<EclipseIO> eclIO_;
|
||||
std::unique_ptr<TaskletRunner> taskletRunner_;
|
||||
Scalar restartTimeStepSize_;
|
||||
};
|
||||
} // namespace Opm
|
||||
|
||||
@ -159,7 +159,7 @@ BOOST_AUTO_TEST_CASE(Summary)
|
||||
|
||||
typedef Opm::EclWriter<TypeTag> EclWriterType;
|
||||
// create the actual ECL writer
|
||||
std::unique_ptr<EclWriterType> eclWriter = std::unique_ptr<EclWriterType>(new EclWriterType(*simulator));
|
||||
std::unique_ptr<EclWriterType> eclWriter = std::unique_ptr<EclWriterType>(new EclWriterType(*simulator, simulator->problem()));
|
||||
|
||||
simulator->model().applyInitialSolution();
|
||||
Opm::data::Wells dw;
|
||||
|
||||
Loading…
Reference in New Issue
Block a user