00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #ifndef OPM_PARALLELDEBUGOUTPUT_HEADER_INCLUDED
00020 #define OPM_PARALLELDEBUGOUTPUT_HEADER_INCLUDED
00021
00022 #include <unordered_set>
00023
00024 #include <opm/common/data/SimulationDataContainer.hpp>
00025
00026
00027 #include <opm/core/grid.h>
00028 #include <opm/core/simulator/WellState.hpp>
00029 #include <opm/core/wells/WellsManager.hpp>
00030
00031 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
00032 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
00033 #include <opm/core/wells/DynamicListEconLimited.hpp>
00034
00035 #if HAVE_OPM_GRID
00036 #include <dune/grid/common/p2pcommunicator.hh>
00037 #endif
00038
00039 namespace Opm
00040 {
00041
00042 class ParallelDebugOutputInterface
00043
00044 {
00045 protected:
00046 ParallelDebugOutputInterface () {}
00047 public:
00048 virtual ~ParallelDebugOutputInterface() {}
00049
00057 virtual bool collectToIORank( const SimulationDataContainer& localReservoirState,
00058 const WellStateFullyImplicitBlackoil& localWellState,
00059 const data::Solution& localCellData,
00060 const int wellStateStepNumber ) = 0;
00061
00062 virtual const SimulationDataContainer& globalReservoirState() const = 0 ;
00063 virtual const data::Solution& globalCellData() const = 0 ;
00064 virtual const WellStateFullyImplicitBlackoil& globalWellState() const = 0 ;
00065 virtual bool isIORank() const = 0;
00066 virtual bool isParallel() const = 0;
00067 virtual int numCells() const = 0 ;
00068 virtual const int* globalCell() const = 0;
00069 };
00070
00071 template <class GridImpl>
00072 class ParallelDebugOutput : public ParallelDebugOutputInterface
00073 {
00074 protected:
00075 const GridImpl& grid_;
00076
00077 const SimulationDataContainer* globalState_;
00078 const WellStateFullyImplicitBlackoil* wellState_;
00079 const data::Solution* globalCellData_;
00080
00081 public:
00082 ParallelDebugOutput ( const GridImpl& grid,
00083 const EclipseState& ,
00084 const int,
00085 const Opm::PhaseUsage& )
00086 : grid_( grid ) {}
00087
00088
00089 virtual bool collectToIORank( const SimulationDataContainer& localReservoirState,
00090 const WellStateFullyImplicitBlackoil& localWellState,
00091 const data::Solution& localCellData,
00092 const int )
00093 {
00094 globalState_ = &localReservoirState;
00095 wellState_ = &localWellState;
00096 globalCellData_ = &localCellData;
00097 return true ;
00098 }
00099
00100 virtual const SimulationDataContainer& globalReservoirState() const { return *globalState_; }
00101 virtual const data::Solution& globalCellData() const
00102 {
00103 return *globalCellData_;
00104 }
00105 virtual const WellStateFullyImplicitBlackoil& globalWellState() const { return *wellState_; }
00106 virtual bool isIORank () const { return true; }
00107 virtual bool isParallel () const { return false; }
00108 virtual int numCells() const { return Opm::AutoDiffGrid::numCells(grid_); }
00109 virtual const int* globalCell() const { return Opm::AutoDiffGrid::globalCell(grid_); }
00110 };
00111
00112 #if HAVE_OPM_GRID
00113 template <>
00114 class ParallelDebugOutput< Dune::CpGrid> : public ParallelDebugOutputInterface
00115 {
00116 public:
00117 typedef Dune::CpGrid Grid;
00118 typedef typename Grid :: CollectiveCommunication CollectiveCommunication;
00119
00120
00121 class GlobalCellIndex
00122 {
00123 int globalId_;
00124 int localIndex_;
00125 bool isInterior_;
00126 public:
00127 GlobalCellIndex() : globalId_(-1), localIndex_(-1), isInterior_(true) {}
00128 void setGhost() { isInterior_ = false; }
00129
00130 void setId( const int globalId ) { globalId_ = globalId; }
00131 void setIndex( const int localIndex ) { localIndex_ = localIndex; }
00132
00133 int localIndex () const { return localIndex_; }
00134 int id () const { return globalId_; }
00135 bool isInterior() const { return isInterior_; }
00136 };
00137
00138 typedef typename Dune::PersistentContainer< Grid, GlobalCellIndex > GlobalIndexContainer;
00139
00140 static const int dimension = Grid :: dimension ;
00141
00142 typedef typename Grid :: LeafGridView GridView;
00143 typedef GridView AllGridView;
00144
00145 typedef Dune :: Point2PointCommunicator< Dune :: SimpleMessageBuffer > P2PCommunicatorType;
00146 typedef typename P2PCommunicatorType :: MessageBufferType MessageBufferType;
00147
00148 typedef std::vector< GlobalCellIndex > LocalIndexMapType;
00149
00150 typedef std::vector<int> IndexMapType;
00151 typedef std::vector< IndexMapType > IndexMapStorageType;
00152
00153 class DistributeIndexMapping : public P2PCommunicatorType::DataHandleInterface
00154 {
00155 protected:
00156 const std::vector<int>& distributedGlobalIndex_;
00157 IndexMapType& localIndexMap_;
00158 IndexMapStorageType& indexMaps_;
00159 std::map< const int, const int > globalPosition_;
00160 #ifndef NDEBUG
00161 std::set< int > checkPosition_;
00162 #endif
00163
00164 public:
00165 DistributeIndexMapping( const std::vector<int>& globalIndex,
00166 const std::vector<int>& distributedGlobalIndex,
00167 IndexMapType& localIndexMap,
00168 IndexMapStorageType& indexMaps )
00169 : distributedGlobalIndex_( distributedGlobalIndex ),
00170 localIndexMap_( localIndexMap ),
00171 indexMaps_( indexMaps ),
00172 globalPosition_()
00173 {
00174 const size_t size = globalIndex.size();
00175
00176 for ( size_t index = 0; index < size; ++index )
00177 {
00178 globalPosition_.insert( std::make_pair( globalIndex[ index ], index ) );
00179 }
00180
00181
00182 if( ! indexMaps_.empty() )
00183 {
00184
00185 IndexMapType& indexMap = indexMaps_.back();
00186 const size_t localSize = localIndexMap_.size();
00187 indexMap.resize( localSize );
00188 for( size_t i=0; i<localSize; ++i )
00189 {
00190 const int id = distributedGlobalIndex_[ localIndexMap_[ i ] ];
00191 indexMap[ i ] = globalPosition_[ id ] ;
00192 #ifndef NDEBUG
00193 assert( checkPosition_.find( id ) == checkPosition_.end() );
00194 checkPosition_.insert( id );
00195 #endif
00196 }
00197 }
00198 }
00199
00200 void pack( const int link, MessageBufferType& buffer )
00201 {
00202
00203 if( link != 0 ) {
00204 OPM_THROW(std::logic_error,"link in method pack is not 0 as execpted");
00205 }
00206
00207
00208 const int size = localIndexMap_.size();
00209 buffer.write( size );
00210
00211 for( int index = 0; index < size; ++index )
00212 {
00213 const int globalIdx = distributedGlobalIndex_[ localIndexMap_[ index ] ];
00214 buffer.write( globalIdx );
00215 }
00216 }
00217
00218 void unpack( const int link, MessageBufferType& buffer )
00219 {
00220
00221 IndexMapType& indexMap = indexMaps_[ link ];
00222 assert( ! globalPosition_.empty() );
00223
00224
00225 int numCells = 0;
00226 buffer.read( numCells );
00227 indexMap.resize( numCells );
00228 for( int index = 0; index < numCells; ++index )
00229 {
00230 int globalId = -1;
00231 buffer.read( globalId );
00232 assert( globalPosition_.find( globalId ) != globalPosition_.end() );
00233 indexMap[ index ] = globalPosition_[ globalId ];
00234 #ifndef NDEBUG
00235 assert( checkPosition_.find( globalId ) == checkPosition_.end() );
00236 checkPosition_.insert( globalId );
00237 #endif
00238 }
00239 }
00240 };
00241
00242 enum { ioRank = 0 };
00243
00249 ParallelDebugOutput( const Dune::CpGrid& otherGrid,
00250 const EclipseState& eclipseState,
00251 const int numPhases,
00252 const Opm::PhaseUsage& phaseUsage)
00253 : grid_(),
00254 eclipseState_( eclipseState ),
00255 globalCellData_(new data::Solution),
00256 isIORank_(true),
00257 phaseUsage_(phaseUsage)
00258
00259 {
00260
00261 Dune::CpGrid distributed_grid = otherGrid;
00262
00263 const CollectiveCommunication& comm = otherGrid.comm();
00264 if( comm.size() > 1 )
00265 {
00266 std::set< int > send, recv;
00267 distributed_grid.switchToDistributedView();
00268 toIORankComm_ = distributed_grid.comm();
00269 isIORank_ = (distributed_grid.comm().rank() == ioRank);
00270
00271
00272 if( isIORank() )
00273 {
00274
00275 grid_.reset( new Dune::CpGrid(otherGrid ) );
00276 grid_->switchToGlobalView();
00277 Dune::CpGrid& globalGrid = *grid_;
00278
00279
00280 globalReservoirState_.reset( new SimulationDataContainer( globalGrid.numCells(), globalGrid.numFaces(), numPhases ));
00281
00282
00283 globalIndex_ = globalGrid.globalCell();
00284
00285 unsigned int count = 0;
00286 auto gridView = globalGrid.leafGridView();
00287 for( auto it = gridView.begin< 0 >(),
00288 end = gridView.end< 0 >(); it != end; ++it, ++count )
00289 {
00290 }
00291 assert( count == globalIndex_.size() );
00292
00293 for(int i=0; i<comm.size(); ++i)
00294 {
00295 if( i != ioRank )
00296 {
00297 recv.insert( i );
00298 }
00299 }
00300 }
00301 else
00302 {
00303
00304
00305 globalReservoirState_.reset( new SimulationDataContainer( 0, 0, 0));
00306 send.insert( ioRank );
00307 }
00308
00309 localIndexMap_.clear();
00310 localIndexMap_.reserve( distributed_grid.size( 0 ) );
00311
00312 unsigned int index = 0;
00313 auto localView = distributed_grid.leafGridView();
00314 for( auto it = localView.begin< 0 >(),
00315 end = localView.end< 0 >(); it != end; ++it, ++index )
00316 {
00317 const auto element = *it ;
00318
00319 if( element.partitionType() == Dune :: InteriorEntity )
00320 {
00321 localIndexMap_.push_back( index );
00322 }
00323 }
00324
00325
00326 toIORankComm_.insertRequest( send, recv );
00327
00328 if( isIORank() )
00329 {
00330
00331 indexMaps_.clear();
00332 indexMaps_.resize( comm.size() );
00333 }
00334
00335
00336 DistributeIndexMapping distIndexMapping( globalIndex_, distributed_grid.globalCell(), localIndexMap_, indexMaps_ );
00337 toIORankComm_.exchange( distIndexMapping );
00338 }
00339 else
00340 {
00341
00342 globalIndex_ = distributed_grid.globalCell();
00343 }
00344 }
00345
00346 class PackUnPackSimulationDataContainer : public P2PCommunicatorType::DataHandleInterface
00347 {
00348 const data::Solution& localCellData_;
00349 data::Solution& globalCellData_;
00350 const WellStateFullyImplicitBlackoil& localWellState_;
00351 WellStateFullyImplicitBlackoil& globalWellState_;
00352 const IndexMapType& localIndexMap_;
00353 const IndexMapStorageType& indexMaps_;
00354
00355 public:
00356 PackUnPackSimulationDataContainer( std::size_t numGlobalCells,
00357 const data::Solution& localCellData,
00358 data::Solution& globalCellData,
00359 const WellStateFullyImplicitBlackoil& localWellState,
00360 WellStateFullyImplicitBlackoil& globalWellState,
00361 const IndexMapType& localIndexMap,
00362 const IndexMapStorageType& indexMaps,
00363 const bool isIORank )
00364 : localCellData_( localCellData ),
00365 globalCellData_( globalCellData ),
00366 localWellState_( localWellState ),
00367 globalWellState_( globalWellState ),
00368 localIndexMap_( localIndexMap ),
00369 indexMaps_( indexMaps )
00370 {
00371
00372 if( isIORank )
00373 {
00374
00375 for (const auto& pair : localCellData_) {
00376 const std::string& key = pair.first;
00377 std::size_t container_size = numGlobalCells;
00378 auto ret = globalCellData_.insert(key, pair.second.dim,
00379 std::vector<double>(container_size),
00380 pair.second.target);
00381 assert(ret.second);
00382 DUNE_UNUSED_PARAMETER(ret.second);
00383 }
00384
00385 MessageBufferType buffer;
00386 pack( 0, buffer );
00387
00388 doUnpack( indexMaps.back(), buffer );
00389 }
00390 }
00391
00392
00393 void pack( const int link, MessageBufferType& buffer )
00394 {
00395
00396 if( link != 0 ) {
00397 OPM_THROW(std::logic_error,"link in method pack is not 0 as execpted");
00398 }
00399
00400
00401 for (const auto& pair : localCellData_) {
00402 const auto& data = pair.second.data;
00403
00404
00405 write( buffer, localIndexMap_, data);
00406 }
00407
00408
00409 writeWells( buffer );
00410 }
00411
00412 void doUnpack( const IndexMapType& indexMap, MessageBufferType& buffer )
00413 {
00414
00415
00416 for (auto& pair : localCellData_) {
00417 const std::string& key = pair.first;
00418 auto& data = globalCellData_.data(key);
00419
00420
00421 read( buffer, indexMap, data);
00422 }
00423
00424
00425 readWells( buffer );
00426 }
00427
00428
00429 void unpack( const int link, MessageBufferType& buffer )
00430 {
00431 doUnpack( indexMaps_[ link ], buffer );
00432 }
00433
00434 protected:
00435 template <class Vector>
00436 void write( MessageBufferType& buffer, const IndexMapType& localIndexMap,
00437 const Vector& vector,
00438 const unsigned int offset = 0, const unsigned int stride = 1 ) const
00439 {
00440 unsigned int size = localIndexMap.size();
00441 buffer.write( size );
00442 assert( vector.size() >= stride * size );
00443 for( unsigned int i=0; i<size; ++i )
00444 {
00445 const unsigned int index = localIndexMap[ i ] * stride + offset;
00446 assert( index < vector.size() );
00447 buffer.write( vector[ index ] );
00448 }
00449 }
00450
00451 template <class Vector>
00452 void read( MessageBufferType& buffer,
00453 const IndexMapType& indexMap,
00454 Vector& vector,
00455 const unsigned int offset = 0, const unsigned int stride = 1 ) const
00456 {
00457 unsigned int size = 0;
00458 buffer.read( size );
00459 assert( size == indexMap.size() );
00460 for( unsigned int i=0; i<size; ++i )
00461 {
00462 const unsigned int index = indexMap[ i ] * stride + offset;
00463 assert( index < vector.size() );
00464 buffer.read( vector[ index ] );
00465 }
00466 }
00467
00468 void writeString( MessageBufferType& buffer, const std::string& s) const
00469 {
00470 const int size = s.size();
00471 buffer.write( size );
00472 for( int i=0; i<size; ++i )
00473 {
00474 buffer.write( s[ i ] );
00475 }
00476 }
00477
00478 void readString( MessageBufferType& buffer, std::string& s) const
00479 {
00480 int size = -1;
00481 buffer.read( size );
00482 s.resize( size );
00483 for( int i=0; i<size; ++i )
00484 {
00485 buffer.read( s[ i ] );
00486 }
00487 }
00488
00489 void writeWells( MessageBufferType& buffer ) const
00490 {
00491 int nWells = localWellState_.wellMap().size();
00492 buffer.write( nWells );
00493 auto end = localWellState_.wellMap().end();
00494 for( auto it = localWellState_.wellMap().begin(); it != end; ++it )
00495 {
00496 const std::string& name = it->first;
00497 const int wellIdx = it->second[ 0 ];
00498
00499
00500 writeString( buffer, name );
00501
00502
00503 buffer.write( localWellState_.bhp()[ wellIdx ] );
00504 buffer.write( localWellState_.thp()[ wellIdx ] );
00505 const int wellRateIdx = wellIdx * localWellState_.numPhases();
00506 for( int np=0; np<localWellState_.numPhases(); ++np )
00507 buffer.write( localWellState_.wellRates()[ wellRateIdx + np ] );
00508
00509
00510 buffer.write(localWellState_.currentControls()[ wellIdx ]);
00511
00512
00513
00514
00515 const int end_con = it->second[1] + it->second[2];
00516
00517 for( int con = it->second[1]; con < end_con; ++con )
00518 {
00519 buffer.write( localWellState_.perfRates()[ con ] );
00520 }
00521
00522 for( int con = it->second[1]; con < end_con; ++con )
00523 {
00524 buffer.write( localWellState_.perfPress()[ con ] );
00525 }
00526
00527
00528 const int np = localWellState_.perfPhaseRates().size() /
00529 localWellState_.perfRates().size();
00530
00531 for( int con = it->second[1]*np; con < end_con*np; ++con )
00532 {
00533 buffer.write( localWellState_.perfPhaseRates()[ con ] );
00534 }
00535 }
00536 }
00537
00538 void readWells( MessageBufferType& buffer )
00539 {
00540 int nWells = -1;
00541 buffer.read( nWells );
00542
00543 std::string name ;
00544 for( int well = 0; well < nWells ; ++well )
00545 {
00546
00547 readString( buffer, name );
00548
00549
00550 auto it = globalWellState_.wellMap().find( name );
00551 if( it == globalWellState_.wellMap().end() )
00552 {
00553 OPM_THROW(std::logic_error,"global state does not contain well " << name );
00554 }
00555 const int wellIdx = it->second[ 0 ];
00556
00557 buffer.read( globalWellState_.bhp()[ wellIdx ] );
00558 buffer.read( globalWellState_.thp()[ wellIdx ] );
00559 const int wellRateIdx = wellIdx * globalWellState_.numPhases();
00560 for( int np=0; np<globalWellState_.numPhases(); ++np )
00561 buffer.read( globalWellState_.wellRates()[ wellRateIdx + np ] );
00562
00563
00564 buffer.read(globalWellState_.currentControls()[ wellIdx ]);
00565
00566
00567
00568
00569 const int end_con = it->second[1] + it->second[2];
00570
00571 for( int con = it->second[1]; con < end_con; ++con )
00572 {
00573 buffer.read( globalWellState_.perfRates()[ con ] );
00574 }
00575
00576 for( int con = it->second[1]; con < end_con; ++con )
00577 {
00578 buffer.read( globalWellState_.perfPress()[ con ] );
00579 }
00580
00581
00582 const int np = globalWellState_.perfPhaseRates().size() /
00583 globalWellState_.perfRates().size();
00584
00585 for( int con = it->second[1]*np; con < end_con*np; ++con )
00586 {
00587 buffer.read( globalWellState_.perfPhaseRates()[ con ] );
00588 }
00589 }
00590 }
00591 };
00592
00593
00594 bool collectToIORank( const SimulationDataContainer& ,
00595 const WellStateFullyImplicitBlackoil& localWellState,
00596 const data::Solution& localCellData,
00597 const int wellStateStepNumber )
00598 {
00599 if( isIORank() )
00600 {
00601 Dune::CpGrid& globalGrid = *grid_;
00602
00603
00604
00605 const DynamicListEconLimited dynamic_list_econ_limited;
00606
00607 WellsManager wells_manager(eclipseState_,
00608 wellStateStepNumber,
00609 Opm::UgGridHelpers::numCells( globalGrid ),
00610 Opm::UgGridHelpers::globalCell( globalGrid ),
00611 Opm::UgGridHelpers::cartDims( globalGrid ),
00612 Opm::UgGridHelpers::dimensions( globalGrid ),
00613 Opm::UgGridHelpers::cell2Faces( globalGrid ),
00614 Opm::UgGridHelpers::beginFaceCentroids( globalGrid ),
00615 dynamic_list_econ_limited,
00616 false,
00617
00618
00619
00620
00621 std::unordered_set<std::string>());
00622
00623 const Wells* wells = wells_manager.c_wells();
00624 globalWellState_.init(wells, *globalReservoirState_, globalWellState_, phaseUsage_ );
00625 globalCellData_->clear();
00626 }
00627
00628 PackUnPackSimulationDataContainer packUnpack( numCells(),
00629 localCellData, *globalCellData_,
00630 localWellState, globalWellState_,
00631 localIndexMap_, indexMaps_,
00632 isIORank() );
00633
00634
00635 toIORankComm_.exchange( packUnpack );
00636 #ifndef NDEBUG
00637
00638 toIORankComm_.barrier();
00639 #endif
00640 if( isIORank() )
00641 {
00642
00643 const std::map<std::string, std::vector<double> > no_extra_data;
00644 solutionToSim(*globalCellData_, no_extra_data, phaseUsage_, *globalReservoirState_);
00645 }
00646 return isIORank();
00647 }
00648
00649 const SimulationDataContainer& globalReservoirState() const { return *globalReservoirState_; }
00650
00651 const data::Solution& globalCellData() const
00652 {
00653 return *globalCellData_;
00654 }
00655
00656 const WellStateFullyImplicitBlackoil& globalWellState() const { return globalWellState_; }
00657
00658 bool isIORank() const
00659 {
00660 return isIORank_;
00661 }
00662
00663 bool isParallel() const
00664 {
00665 return toIORankComm_.size() > 1;
00666 }
00667
00668 int numCells () const { return globalIndex_.size(); }
00669 const int* globalCell () const
00670 {
00671 assert( ! globalIndex_.empty() );
00672 return globalIndex_.data();
00673 }
00674
00675 protected:
00676 std::unique_ptr< Dune::CpGrid > grid_;
00677 const EclipseState& eclipseState_;
00678 P2PCommunicatorType toIORankComm_;
00679 IndexMapType globalIndex_;
00680 IndexMapType localIndexMap_;
00681 IndexMapStorageType indexMaps_;
00682 std::unique_ptr<SimulationDataContainer> globalReservoirState_;
00683 std::unique_ptr<data::Solution> globalCellData_;
00684
00685 WellStateFullyImplicitBlackoil globalWellState_;
00686
00687 bool isIORank_;
00688
00689 Opm::PhaseUsage phaseUsage_;
00690 };
00691 #endif // #if HAVE_OPM_GRID
00692
00693 }
00694 #endif