00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #ifndef OPM_REDISTRIBUTEDATAHANDLES_HEADER
00023 #define OPM_REDISTRIBUTEDATAHANDLES_HEADER
00024
00025 #include <unordered_set>
00026 #include <string>
00027 #include <type_traits>
00028 #include <iterator>
00029
00030 #include <opm/core/simulator/BlackoilState.hpp>
00031
00032 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
00033 #include <opm/autodiff/ExtractParallelGridInformationToISTL.hpp>
00034 #include <opm/autodiff/createGlobalCellArray.hpp>
00035
00036 #include<boost/any.hpp>
00037
00038 namespace Opm
00039 {
00040
00041 template <class Grid>
00042 inline std::unordered_set<std::string>
00043 distributeGridAndData( Grid& ,
00044 const Opm::Deck& ,
00045 const EclipseState& ,
00046 BlackoilState& ,
00047 BlackoilPropsAdFromDeck& ,
00048 DerivedGeology&,
00049 std::shared_ptr<BlackoilPropsAdFromDeck::MaterialLawManager>&,
00050 std::vector<double>&,
00051 boost::any& ,
00052 const bool )
00053 {
00054 return std::unordered_set<std::string>();
00055 }
00056
00062 template<class Iter1, class Iter2=Iter1>
00063 class FixedSizeIterCopyHandle
00064 {
00065 typedef typename std::iterator_traits<Iter1>::value_type DataType2;
00066
00067 public:
00068 typedef typename std::iterator_traits<Iter1>::value_type DataType;
00069
00073 FixedSizeIterCopyHandle(const Iter1& send_begin,
00074 const Iter2& receive_begin,
00075 std::size_t size = 1)
00076 : send_(send_begin), receive_(receive_begin), size_(size)
00077 {
00078 static_assert(std::is_same<DataType,DataType2>::value,
00079 "Iter1 and Iter2 have to have the same value_type!");
00080 }
00081
00082 template<class Buffer>
00083 void gather(Buffer& buffer, std::size_t i)
00084 {
00085 for(auto index = i*size(i), end = (i+1)*size(i);
00086 index < end; ++index)
00087 {
00088 buffer.write(send_[index]);
00089 }
00090 }
00091
00092 template<class Buffer>
00093 void scatter(Buffer& buffer, std::size_t i, std::size_t s OPM_OPTIM_UNUSED)
00094 {
00095 assert(s==size(i));
00096 static_cast<void>(s);
00097
00098 for(auto index = i*size(i), end = (i+1)*size(i);
00099 index < end; ++index)
00100 {
00101 buffer.read(receive_[index]);
00102 }
00103 }
00104
00105 bool fixedsize()
00106 {
00107 return true;
00108 }
00109
00110 std::size_t size(std::size_t)
00111 {
00112 return size_;
00113 }
00114 private:
00115 Iter1 send_;
00116 Iter2 receive_;
00117 std::size_t size_;
00118 };
00119
00120
00121 #if HAVE_OPM_GRID && HAVE_MPI
00123 class ThresholdPressureDataHandle
00124 {
00125 public:
00127 typedef double DataType;
00133 ThresholdPressureDataHandle(const Dune::CpGrid& sendGrid,
00134 const Dune::CpGrid& recvGrid,
00135 const std::vector<double>& sendPressures,
00136 std::vector<double>& recvPressures)
00137 : sendGrid_(sendGrid), recvGrid_(recvGrid), sendPressures_(sendPressures),
00138 recvPressures_(recvPressures)
00139 {}
00140
00141 bool fixedsize(int , int )
00142 {
00143 return false;
00144 }
00145 template<class T>
00146 std::size_t size(const T& e)
00147 {
00148 if ( T::codimension == 0)
00149 {
00150 return sendGrid_.numCellFaces(e.index());
00151 }
00152 else
00153 {
00154 OPM_THROW(std::logic_error, "Data handle can only be used for elements");
00155 }
00156 }
00157 template<class B, class T>
00158 void gather(B& buffer, const T& e)
00159 {
00160 assert( T::codimension == 0);
00161 for ( int i=0; i< sendGrid_.numCellFaces(e.index()); ++i )
00162 {
00163 buffer.write(sendPressures_[sendGrid_.cellFace(e.index(), i)]);
00164 }
00165 }
00166 template<class B, class T>
00167 void scatter(B& buffer, const T& e, std::size_t )
00168 {
00169 assert( T::codimension == 0);
00170 for ( int i=0; i< recvGrid_.numCellFaces(e.index()); ++i )
00171 {
00172 double val;
00173 buffer.read(val);
00174 recvPressures_[recvGrid_.cellFace(e.index(), i)]=val;
00175 }
00176 }
00177 bool contains(int dim, int codim)
00178 {
00179 return dim==3 && codim==0;
00180 }
00181
00182 private:
00184 const Dune::CpGrid& sendGrid_;
00186 const Dune::CpGrid& recvGrid_;
00188 const std::vector<double>& sendPressures_;
00190 std::vector<double>& recvPressures_;
00191 };
00192
00194 class GeologyDataHandle
00195 {
00196 public:
00198 typedef double DataType;
00204 GeologyDataHandle(const Dune::CpGrid& sendGrid,
00205 const Dune::CpGrid& recvGrid,
00206 const DerivedGeology& sendGeology,
00207 DerivedGeology& recvGeology)
00208 : sendGrid_(sendGrid), recvGrid_(recvGrid), sendGeology_(sendGeology),
00209 recvGeology_(recvGeology)
00210 {}
00211
00212 bool fixedsize(int , int )
00213 {
00214 return false;
00215 }
00216 template<class T>
00217 std::size_t size(const T& e)
00218 {
00219 if ( T::codimension == 0)
00220 {
00221 return 1 + sendGrid_.numCellFaces(e.index());
00222 }
00223 else
00224 {
00225 OPM_THROW(std::logic_error, "Data handle can only be used for elements");
00226 }
00227 }
00228 template<class B, class T>
00229 void gather(B& buffer, const T& e)
00230 {
00231 assert( T::codimension == 0);
00232 buffer.write(sendGeology_.poreVolume()[e.index()]);
00233 for ( int i=0; i< sendGrid_.numCellFaces(e.index()); ++i )
00234 {
00235 buffer.write(sendGeology_.transmissibility()[sendGrid_.cellFace(e.index(), i)]);
00236 }
00237 }
00238 template<class B, class T>
00239 void scatter(B& buffer, const T& e, std::size_t )
00240 {
00241 assert( T::codimension == 0);
00242 double val;
00243 buffer.read(val);
00244 recvGeology_.poreVolume()[e.index()]=val;
00245 for ( int i=0; i< recvGrid_.numCellFaces(e.index()); ++i )
00246 {
00247 buffer.read(val);
00248 recvGeology_.transmissibility()[recvGrid_.cellFace(e.index(), i)]=val;
00249 }
00250 }
00251 bool contains(int dim, int codim)
00252 {
00253 return dim==3 && codim==0;
00254 }
00255 private:
00257 const Dune::CpGrid& sendGrid_;
00259 const Dune::CpGrid& recvGrid_;
00261 const DerivedGeology& sendGeology_;
00263 DerivedGeology& recvGeology_;
00264 };
00265
00267 class BlackoilStateDataHandle
00268 {
00269 public:
00271 typedef double DataType;
00277 BlackoilStateDataHandle(const Dune::CpGrid& sendGrid,
00278 const Dune::CpGrid& recvGrid,
00279 const BlackoilState& sendState,
00280 BlackoilState& recvState)
00281 : sendGrid_(sendGrid), recvGrid_(recvGrid), sendState_(sendState), recvState_(recvState)
00282 {
00283
00284 recvState.surfacevol().resize(recvGrid.numCells()*sendState.numPhases(),
00285 std::numeric_limits<double>::max());
00286 recvState.hydroCarbonState().resize(recvGrid.numCells());
00287 }
00288
00289 bool fixedsize(int , int )
00290 {
00291 return false;
00292 }
00293
00294 template<class T>
00295 std::size_t size(const T& e)
00296 {
00297 if ( T::codimension == 0)
00298 {
00299 return 2 * sendState_.numPhases() + 5 + 2*sendGrid_.numCellFaces(e.index());
00300 }
00301 else
00302 {
00303 OPM_THROW(std::logic_error, "Data handle can only be used for elements");
00304 }
00305 }
00306
00307 template<class B, class T>
00308 void gather(B& buffer, const T& e)
00309 {
00310 assert( T::codimension == 0);
00311
00312 for ( size_t i=0; i<sendState_.numPhases(); ++i )
00313 {
00314 buffer.write(sendState_.surfacevol()[e.index()*sendState_.numPhases()+i]);
00315 }
00316 buffer.write(sendState_.gasoilratio()[e.index()]);
00317 buffer.write(sendState_.rv()[e.index()]);
00318 buffer.write(sendState_.pressure()[e.index()]);
00319 buffer.write(sendState_.temperature()[e.index()]);
00320
00321 double hydroCarbonState_ = sendState_.hydroCarbonState()[e.index()];
00322 buffer.write(hydroCarbonState_);
00323
00324 for ( size_t i=0; i<sendState_.numPhases(); ++i )
00325 {
00326 buffer.write(sendState_.saturation()[e.index()*sendState_.numPhases()+i]);
00327 }
00328 for ( int i=0; i<sendGrid_.numCellFaces(e.index()); ++i )
00329 {
00330 buffer.write(sendState_.facepressure()[sendGrid_.cellFace(e.index(), i)]);
00331 }
00332 for ( int i=0; i<sendGrid_.numCellFaces(e.index()); ++i )
00333 {
00334 buffer.write(sendState_.faceflux()[sendGrid_.cellFace(e.index(), i)]);
00335 }
00336 }
00337 template<class B, class T>
00338 void scatter(B& buffer, const T& e, std::size_t size_arg)
00339 {
00340 assert( T::codimension == 0);
00341 assert( size_arg == 2 * recvState_.numPhases() + 5 +2*recvGrid_.numCellFaces(e.index()));
00342 static_cast<void>(size_arg);
00343
00344 double val;
00345 for ( size_t i=0; i<recvState_.numPhases(); ++i )
00346 {
00347 buffer.read(val);
00348 recvState_.surfacevol()[e.index()*sendState_.numPhases()+i]=val;
00349 }
00350 buffer.read(val);
00351 recvState_.gasoilratio()[e.index()]=val;
00352 buffer.read(val);
00353 recvState_.rv()[e.index()]=val;
00354 buffer.read(val);
00355 recvState_.pressure()[e.index()]=val;
00356 buffer.read(val);
00357 recvState_.temperature()[e.index()]=val;
00358
00359 buffer.read(val);
00360 recvState_.hydroCarbonState()[e.index()]=static_cast<HydroCarbonState>(val);
00361
00362 for ( size_t i=0; i<recvState_.numPhases(); ++i )
00363 {
00364 buffer.read(val);
00365 recvState_.saturation()[e.index()*sendState_.numPhases()+i]=val;
00366 }
00367 for ( int i=0; i<recvGrid_.numCellFaces(e.index()); ++i )
00368 {
00369 buffer.read(val);
00370 recvState_.facepressure()[recvGrid_.cellFace(e.index(), i)]=val;
00371 }
00372 for ( int i=0; i<recvGrid_.numCellFaces(e.index()); ++i )
00373 {
00374 buffer.read(val);
00375 recvState_.faceflux()[recvGrid_.cellFace(e.index(), i)]=val;
00376 }
00377 }
00378 bool contains(int dim, int codim)
00379 {
00380 return dim==3 && codim==0;
00381 }
00382 private:
00384 const Dune::CpGrid& sendGrid_;
00386 const Dune::CpGrid& recvGrid_;
00388 const BlackoilState& sendState_;
00389
00390 BlackoilState& recvState_;
00391 };
00392
00394 class BlackoilPropsDataHandle
00395 {
00396 public:
00398 typedef double DataType;
00402 BlackoilPropsDataHandle(const BlackoilPropsAdFromDeck& sendProps,
00403 BlackoilPropsAdFromDeck& recvProps)
00404 : sendProps_(sendProps), recvProps_(recvProps),
00405 size_(11)
00406 {
00407
00408 if ( sendProps.satOilMax_.size()>0 )
00409 {
00410
00411
00412 recvProps_.satOilMax_.resize(recvProps_.cellPvtRegionIdx_.size(),
00413 -std::numeric_limits<double>::max());
00414 ++size_;
00415 }
00416 }
00417
00418 bool fixedsize(int , int )
00419 {
00420 return true;
00421 }
00422
00423 template<class T>
00424 std::size_t size(const T&)
00425 {
00426 if ( T::codimension == 0)
00427 {
00428 return size_;
00429 }
00430 else
00431 {
00432 OPM_THROW(std::logic_error, "Data handle can only be used for elements");
00433 }
00434 }
00435 template<class B, class T>
00436 void gather(B& buffer, const T& e)
00437 {
00438 assert( T::codimension == 0);
00439
00440 buffer.write(sendProps_.cellPvtRegionIndex()[e.index()]);
00441 for( std::size_t i = 0; i < 9; ++i )
00442 {
00443 buffer.write(sendProps_.rock_.permeability_[e.index()*9+i]);
00444 }
00445 buffer.write(sendProps_.rock_.porosity_[e.index()]);
00446 if ( size_ > 11 ) {
00447 buffer.write(sendProps_.satOilMax_[e.index()]);
00448 }
00449 }
00450 template<class B, class T>
00451 void scatter(B& buffer, const T& e, std::size_t size_arg)
00452 {
00453 assert( T::codimension == 0);
00454 assert( size_arg==size_ ); (void) size_arg;
00455 double val;
00456 buffer.read(val);
00457 recvProps_.cellPvtRegionIdx_[e.index()]=val;
00458 for( std::size_t i = 0; i < 9; ++i )
00459 {
00460 buffer.read(val);
00461 recvProps_.rock_.permeability_[e.index()*9+i]
00462 = val;
00463 }
00464 buffer.read(val);
00465 recvProps_.rock_.porosity_[e.index()]=val;
00466 if ( size_ > 11 ) {
00467 buffer.read(val);
00468 recvProps_.satOilMax_[e.index()]=val;
00469 }
00470 }
00471 bool contains(int dim, int codim)
00472 {
00473 return dim==3 && codim==0;
00474 }
00475 private:
00477 const BlackoilPropsAdFromDeck& sendProps_;
00479 BlackoilPropsAdFromDeck& recvProps_;
00484 std::size_t size_;
00485 };
00486
00487 inline
00488 std::unordered_set<std::string>
00489 distributeGridAndData( Dune::CpGrid& grid,
00490 const Opm::Deck& deck,
00491 const EclipseState& eclipseState,
00492 BlackoilState& state,
00493 BlackoilPropsAdFromDeck& properties,
00494 DerivedGeology& geology,
00495 std::shared_ptr<BlackoilPropsAdFromDeck::MaterialLawManager>& material_law_manager,
00496 std::vector<double>& threshold_pressures,
00497 boost::any& parallelInformation,
00498 const bool useLocalPerm)
00499 {
00500 Dune::CpGrid global_grid ( grid );
00501 global_grid.switchToGlobalView();
00502
00503
00504 using std::get;
00505 auto my_defunct_wells = get<1>(grid.loadBalance(&eclipseState,
00506 geology.transmissibility().data()));
00507 grid.switchToDistributedView();
00508 std::vector<int> compressedToCartesianIdx;
00509 Opm::createGlobalCellArray(grid, compressedToCartesianIdx);
00510 typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager;
00511 auto distributed_material_law_manager = std::make_shared<MaterialLawManager>();
00512 distributed_material_law_manager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
00513
00514
00515
00516
00517 typedef Dune::CpGrid::ParallelIndexSet IndexSet;
00518 const IndexSet& local_indices = grid.getCellIndexSet();
00519 for ( auto index : local_indices )
00520 {
00521 distributed_material_law_manager->materialLawParamsPointerReferenceHack(index.local()) =
00522 material_law_manager->materialLawParamsPointerReferenceHack(index.global());
00523
00524 distributed_material_law_manager->oilWaterScaledEpsInfoDrainagePointerReferenceHack(index.local()) =
00525 material_law_manager->oilWaterScaledEpsInfoDrainagePointerReferenceHack(index.global());
00526 }
00527 BlackoilPropsAdFromDeck distributed_props(properties,
00528 distributed_material_law_manager,
00529 grid.numCells());
00530 BlackoilState distributed_state(grid.numCells(), grid.numFaces(), state.numPhases());
00531 BlackoilStateDataHandle state_handle(global_grid, grid,
00532 state, distributed_state);
00533 BlackoilPropsDataHandle props_handle(properties,
00534 distributed_props);
00535 grid.scatterData(state_handle);
00536 grid.scatterData(props_handle);
00537
00538
00539 DerivedGeology distributed_geology(grid,
00540 distributed_props, eclipseState,
00541 useLocalPerm, geology.gravity());
00542 GeologyDataHandle geo_handle(global_grid, grid,
00543 geology, distributed_geology);
00544 grid.scatterData(geo_handle);
00545
00546 std::vector<double> distributed_pressures;
00547
00548 if( !threshold_pressures.empty() )
00549 {
00550 if( threshold_pressures.size() !=
00551 static_cast<std::size_t>(UgGridHelpers::numFaces(global_grid)) )
00552 {
00553 OPM_THROW(std::runtime_error, "NNCs not yet supported for parallel runs. "
00554 << UgGridHelpers::numFaces(grid) << " faces but " <<
00555 threshold_pressures.size()<<" threshold pressure values");
00556 }
00557 distributed_pressures.resize(UgGridHelpers::numFaces(grid));
00558 ThresholdPressureDataHandle press_handle(global_grid, grid,
00559 threshold_pressures,
00560 distributed_pressures);
00561 grid.scatterData(press_handle);
00562 }
00563
00564
00565 properties = distributed_props;
00566 geology = distributed_geology;
00567 state = distributed_state;
00568 material_law_manager = distributed_material_law_manager;
00569 threshold_pressures = distributed_pressures;
00570 extractParallelGridInformationToISTL(grid, parallelInformation);
00571
00572 return my_defunct_wells;
00573 }
00574 #endif
00575
00576 }
00577 #endif