All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
RedistributeDataHandles.hpp
1 /*
2  Copyright 2015-2016 Dr. Blatt - HPC-Simulation-Software & Services.
3  Coypright 2015 NTNU
4  Copyright 2015-2016 Statoil AS
5  Copyright 2015 IRIS AS
6 
7  This file is part of the Open Porous Media project (OPM).
8 
9  OPM is free software: you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation, either version 3 of the License, or
12  (at your option) any later version.
13 
14  OPM is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU General Public License for more details.
18 
19  You should have received a copy of the GNU General Public License
20  along with OPM. If not, see <http://www.gnu.org/licenses/>.
21 */
22 #ifndef OPM_REDISTRIBUTEDATAHANDLES_HEADER
23 #define OPM_REDISTRIBUTEDATAHANDLES_HEADER
24 
25 #include <unordered_set>
26 #include <string>
27 #include <type_traits>
28 #include <iterator>
29 
30 #include <opm/core/simulator/BlackoilState.hpp>
31 
32 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
33 #include <opm/autodiff/ExtractParallelGridInformationToISTL.hpp>
34 #include <opm/autodiff/createGlobalCellArray.hpp>
35 
36 #include<boost/any.hpp>
37 
38 namespace Opm
39 {
40 
41 template <class Grid>
42 inline std::unordered_set<std::string>
43 distributeGridAndData( Grid& ,
44  const Opm::Deck& ,
45  const EclipseState& ,
46  BlackoilState& ,
47  BlackoilPropsAdFromDeck& ,
48  DerivedGeology&,
49  std::shared_ptr<BlackoilPropsAdFromDeck::MaterialLawManager>&,
50  std::vector<double>&,
51  boost::any& ,
52  const bool )
53 {
54  return std::unordered_set<std::string>();
55 }
56 
62 template<class Iter1, class Iter2=Iter1>
64 {
65  typedef typename std::iterator_traits<Iter1>::value_type DataType2;
66 
67 public:
68  typedef typename std::iterator_traits<Iter1>::value_type DataType;
69 
73  FixedSizeIterCopyHandle(const Iter1& send_begin,
74  const Iter2& receive_begin,
75  std::size_t size = 1)
76  : send_(send_begin), receive_(receive_begin), size_(size)
77  {
78  static_assert(std::is_same<DataType,DataType2>::value,
79  "Iter1 and Iter2 have to have the same value_type!");
80  }
81 
82  template<class Buffer>
83  void gather(Buffer& buffer, std::size_t i)
84  {
85  for(auto index = i*size(i), end = (i+1)*size(i);
86  index < end; ++index)
87  {
88  buffer.write(send_[index]);
89  }
90  }
91 
92  template<class Buffer>
93  void scatter(Buffer& buffer, std::size_t i, std::size_t s OPM_OPTIM_UNUSED)
94  {
95  assert(s==size(i));
96  static_cast<void>(s);
97 
98  for(auto index = i*size(i), end = (i+1)*size(i);
99  index < end; ++index)
100  {
101  buffer.read(receive_[index]);
102  }
103  }
104 
105  bool fixedsize()
106  {
107  return true;
108  }
109 
110  std::size_t size(std::size_t)
111  {
112  return size_;
113  }
114 private:
115  Iter1 send_;
116  Iter2 receive_;
117  std::size_t size_;
118 };
119 
120 
121 #if HAVE_OPM_GRID && HAVE_MPI
122 class ThresholdPressureDataHandle
124 {
125 public:
127  typedef double DataType;
133  ThresholdPressureDataHandle(const Dune::CpGrid& sendGrid,
134  const Dune::CpGrid& recvGrid,
135  const std::vector<double>& sendPressures,
136  std::vector<double>& recvPressures)
137  : sendGrid_(sendGrid), recvGrid_(recvGrid), sendPressures_(sendPressures),
138  recvPressures_(recvPressures)
139  {}
140 
141  bool fixedsize(int /*dim*/, int /*codim*/)
142  {
143  return false;
144  }
145  template<class T>
146  std::size_t size(const T& e)
147  {
148  if ( T::codimension == 0)
149  {
150  return sendGrid_.numCellFaces(e.index());
151  }
152  else
153  {
154  OPM_THROW(std::logic_error, "Data handle can only be used for elements");
155  }
156  }
157  template<class B, class T>
158  void gather(B& buffer, const T& e)
159  {
160  assert( T::codimension == 0);
161  for ( int i=0; i< sendGrid_.numCellFaces(e.index()); ++i )
162  {
163  buffer.write(sendPressures_[sendGrid_.cellFace(e.index(), i)]);
164  }
165  }
166  template<class B, class T>
167  void scatter(B& buffer, const T& e, std::size_t /* size */)
168  {
169  assert( T::codimension == 0);
170  for ( int i=0; i< recvGrid_.numCellFaces(e.index()); ++i )
171  {
172  double val;
173  buffer.read(val);
174  recvPressures_[recvGrid_.cellFace(e.index(), i)]=val;
175  }
176  }
177  bool contains(int dim, int codim)
178  {
179  return dim==3 && codim==0;
180  }
181 
182 private:
184  const Dune::CpGrid& sendGrid_;
186  const Dune::CpGrid& recvGrid_;
188  const std::vector<double>& sendPressures_;
190  std::vector<double>& recvPressures_;
191 };
192 
194 class GeologyDataHandle
195 {
196 public:
198  typedef double DataType;
204  GeologyDataHandle(const Dune::CpGrid& sendGrid,
205  const Dune::CpGrid& recvGrid,
206  const DerivedGeology& sendGeology,
207  DerivedGeology& recvGeology)
208  : sendGrid_(sendGrid), recvGrid_(recvGrid), sendGeology_(sendGeology),
209  recvGeology_(recvGeology)
210  {}
211 
212  bool fixedsize(int /*dim*/, int /*codim*/)
213  {
214  return false;
215  }
216  template<class T>
217  std::size_t size(const T& e)
218  {
219  if ( T::codimension == 0)
220  {
221  return 1 + sendGrid_.numCellFaces(e.index());
222  }
223  else
224  {
225  OPM_THROW(std::logic_error, "Data handle can only be used for elements");
226  }
227  }
228  template<class B, class T>
229  void gather(B& buffer, const T& e)
230  {
231  assert( T::codimension == 0);
232  buffer.write(sendGeology_.poreVolume()[e.index()]);
233  for ( int i=0; i< sendGrid_.numCellFaces(e.index()); ++i )
234  {
235  buffer.write(sendGeology_.transmissibility()[sendGrid_.cellFace(e.index(), i)]);
236  }
237  }
238  template<class B, class T>
239  void scatter(B& buffer, const T& e, std::size_t /* size */)
240  {
241  assert( T::codimension == 0);
242  double val;
243  buffer.read(val);
244  recvGeology_.poreVolume()[e.index()]=val;
245  for ( int i=0; i< recvGrid_.numCellFaces(e.index()); ++i )
246  {
247  buffer.read(val);
248  recvGeology_.transmissibility()[recvGrid_.cellFace(e.index(), i)]=val;
249  }
250  }
251  bool contains(int dim, int codim)
252  {
253  return dim==3 && codim==0;
254  }
255 private:
257  const Dune::CpGrid& sendGrid_;
259  const Dune::CpGrid& recvGrid_;
261  const DerivedGeology& sendGeology_;
263  DerivedGeology& recvGeology_;
264 };
265 
267 class BlackoilStateDataHandle
268 {
269 public:
271  typedef double DataType;
277  BlackoilStateDataHandle(const Dune::CpGrid& sendGrid,
278  const Dune::CpGrid& recvGrid,
279  const BlackoilState& sendState,
280  BlackoilState& recvState)
281  : sendGrid_(sendGrid), recvGrid_(recvGrid), sendState_(sendState), recvState_(recvState)
282  {
283  // construction does not resize surfacevol and hydroCarbonState. Do it manually.
284  recvState.surfacevol().resize(recvGrid.numCells()*sendState.numPhases(),
285  std::numeric_limits<double>::max());
286  recvState.hydroCarbonState().resize(recvGrid.numCells());
287  }
288 
289  bool fixedsize(int /*dim*/, int /*codim*/)
290  {
291  return false;
292  }
293 
294  template<class T>
295  std::size_t size(const T& e)
296  {
297  if ( T::codimension == 0)
298  {
299  return 2 * sendState_.numPhases() + 5 + 2*sendGrid_.numCellFaces(e.index());
300  }
301  else
302  {
303  OPM_THROW(std::logic_error, "Data handle can only be used for elements");
304  }
305  }
306 
307  template<class B, class T>
308  void gather(B& buffer, const T& e)
309  {
310  assert( T::codimension == 0);
311 
312  for ( size_t i=0; i<sendState_.numPhases(); ++i )
313  {
314  buffer.write(sendState_.surfacevol()[e.index()*sendState_.numPhases()+i]);
315  }
316  buffer.write(sendState_.gasoilratio()[e.index()]);
317  buffer.write(sendState_.rv()[e.index()]);
318  buffer.write(sendState_.pressure()[e.index()]);
319  buffer.write(sendState_.temperature()[e.index()]);
320  //We can only send one type with this buffer. Ergo we convert the enum to a double.
321  double hydroCarbonState_ = sendState_.hydroCarbonState()[e.index()];
322  buffer.write(hydroCarbonState_);
323 
324  for ( size_t i=0; i<sendState_.numPhases(); ++i )
325  {
326  buffer.write(sendState_.saturation()[e.index()*sendState_.numPhases()+i]);
327  }
328  for ( int i=0; i<sendGrid_.numCellFaces(e.index()); ++i )
329  {
330  buffer.write(sendState_.facepressure()[sendGrid_.cellFace(e.index(), i)]);
331  }
332  for ( int i=0; i<sendGrid_.numCellFaces(e.index()); ++i )
333  {
334  buffer.write(sendState_.faceflux()[sendGrid_.cellFace(e.index(), i)]);
335  }
336  }
337  template<class B, class T>
338  void scatter(B& buffer, const T& e, std::size_t size_arg)
339  {
340  assert( T::codimension == 0);
341  assert( size_arg == 2 * recvState_.numPhases() + 5 +2*recvGrid_.numCellFaces(e.index()));
342  static_cast<void>(size_arg);
343 
344  double val;
345  for ( size_t i=0; i<recvState_.numPhases(); ++i )
346  {
347  buffer.read(val);
348  recvState_.surfacevol()[e.index()*sendState_.numPhases()+i]=val;
349  }
350  buffer.read(val);
351  recvState_.gasoilratio()[e.index()]=val;
352  buffer.read(val);
353  recvState_.rv()[e.index()]=val;
354  buffer.read(val);
355  recvState_.pressure()[e.index()]=val;
356  buffer.read(val);
357  recvState_.temperature()[e.index()]=val;
358  //We can only send one type with this buffer. Ergo we convert the enum to a double.
359  buffer.read(val);
360  recvState_.hydroCarbonState()[e.index()]=static_cast<HydroCarbonState>(val);
361 
362  for ( size_t i=0; i<recvState_.numPhases(); ++i )
363  {
364  buffer.read(val);
365  recvState_.saturation()[e.index()*sendState_.numPhases()+i]=val;
366  }
367  for ( int i=0; i<recvGrid_.numCellFaces(e.index()); ++i )
368  {
369  buffer.read(val);
370  recvState_.facepressure()[recvGrid_.cellFace(e.index(), i)]=val;
371  }
372  for ( int i=0; i<recvGrid_.numCellFaces(e.index()); ++i )
373  {
374  buffer.read(val);
375  recvState_.faceflux()[recvGrid_.cellFace(e.index(), i)]=val;
376  }
377  }
378  bool contains(int dim, int codim)
379  {
380  return dim==3 && codim==0;
381  }
382 private:
384  const Dune::CpGrid& sendGrid_;
386  const Dune::CpGrid& recvGrid_;
388  const BlackoilState& sendState_;
389  // \brief The state where we will store the received values.
390  BlackoilState& recvState_;
391 };
392 
394 class BlackoilPropsDataHandle
395 {
396 public:
398  typedef double DataType;
402  BlackoilPropsDataHandle(const BlackoilPropsAdFromDeck& sendProps,
403  BlackoilPropsAdFromDeck& recvProps)
404  : sendProps_(sendProps), recvProps_(recvProps),
405  size_(11) // full permeability tensor 9 + porosity 1 + pvt region index
406  {
407  // satOilMax might be non empty. In this case we will need to send it, too.
408  if ( sendProps.satOilMax_.size()>0 )
409  {
410 
411  // satOilMax has to have the same size as the cellPvtRegionIdx_
412  recvProps_.satOilMax_.resize(recvProps_.cellPvtRegionIdx_.size(),
413  -std::numeric_limits<double>::max());
414  ++size_;
415  }
416  }
417 
418  bool fixedsize(int /*dim*/, int /*codim*/)
419  {
420  return true;
421  }
422 
423  template<class T>
424  std::size_t size(const T&)
425  {
426  if ( T::codimension == 0)
427  {
428  return size_;
429  }
430  else
431  {
432  OPM_THROW(std::logic_error, "Data handle can only be used for elements");
433  }
434  }
435  template<class B, class T>
436  void gather(B& buffer, const T& e)
437  {
438  assert( T::codimension == 0);
439 
440  buffer.write(sendProps_.cellPvtRegionIndex()[e.index()]);
441  for( std::size_t i = 0; i < 9; ++i )
442  {
443  buffer.write(sendProps_.rock_.permeability_[e.index()*9+i]);
444  }
445  buffer.write(sendProps_.rock_.porosity_[e.index()]);
446  if ( size_ > 11 ) {
447  buffer.write(sendProps_.satOilMax_[e.index()]);
448  }
449  }
450  template<class B, class T>
451  void scatter(B& buffer, const T& e, std::size_t size_arg)
452  {
453  assert( T::codimension == 0);
454  assert( size_arg==size_ ); (void) size_arg;
455  double val;
456  buffer.read(val);
457  recvProps_.cellPvtRegionIdx_[e.index()]=val;
458  for( std::size_t i = 0; i < 9; ++i )
459  {
460  buffer.read(val);
461  recvProps_.rock_.permeability_[e.index()*9+i]
462  = val;
463  }
464  buffer.read(val);
465  recvProps_.rock_.porosity_[e.index()]=val;
466  if ( size_ > 11 ) {
467  buffer.read(val);
468  recvProps_.satOilMax_[e.index()]=val;
469  }
470  }
471  bool contains(int dim, int codim)
472  {
473  return dim==3 && codim==0;
474  }
475 private:
477  const BlackoilPropsAdFromDeck& sendProps_;
479  BlackoilPropsAdFromDeck& recvProps_;
484  std::size_t size_;
485 };
486 
487 inline
488 std::unordered_set<std::string>
489 distributeGridAndData( Dune::CpGrid& grid,
490  const Opm::Deck& deck,
491  const EclipseState& eclipseState,
492  BlackoilState& state,
493  BlackoilPropsAdFromDeck& properties,
494  DerivedGeology& geology,
495  std::shared_ptr<BlackoilPropsAdFromDeck::MaterialLawManager>& material_law_manager,
496  std::vector<double>& threshold_pressures,
497  boost::any& parallelInformation,
498  const bool useLocalPerm)
499 {
500  Dune::CpGrid global_grid ( grid );
501  global_grid.switchToGlobalView();
502 
503  // distribute the grid and switch to the distributed view
504  using std::get;
505  auto my_defunct_wells = get<1>(grid.loadBalance(&eclipseState,
506  geology.transmissibility().data()));
507  grid.switchToDistributedView();
508  std::vector<int> compressedToCartesianIdx;
509  Opm::createGlobalCellArray(grid, compressedToCartesianIdx);
510  typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager;
511  auto distributed_material_law_manager = std::make_shared<MaterialLawManager>();
512  distributed_material_law_manager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
513  // copy the values from the global to the local MaterialLawManager
514  // We should actually communicate these to be future proof. But that is
515  // really, really cumbersome for the underlying vector<shared_ptr>
516  // where the classes pointed to even have more shared_ptr stored in them.
517  typedef Dune::CpGrid::ParallelIndexSet IndexSet;
518  const IndexSet& local_indices = grid.getCellIndexSet();
519  for ( auto index : local_indices )
520  {
521  distributed_material_law_manager->materialLawParamsPointerReferenceHack(index.local()) =
522  material_law_manager->materialLawParamsPointerReferenceHack(index.global());
523 
524  distributed_material_law_manager->oilWaterScaledEpsInfoDrainagePointerReferenceHack(index.local()) =
525  material_law_manager->oilWaterScaledEpsInfoDrainagePointerReferenceHack(index.global());
526  }
527  BlackoilPropsAdFromDeck distributed_props(properties,
528  distributed_material_law_manager,
529  grid.numCells());
530  BlackoilState distributed_state(grid.numCells(), grid.numFaces(), state.numPhases());
531  BlackoilStateDataHandle state_handle(global_grid, grid,
532  state, distributed_state);
533  BlackoilPropsDataHandle props_handle(properties,
534  distributed_props);
535  grid.scatterData(state_handle);
536  grid.scatterData(props_handle);
537  // Create a distributed Geology. Some values will be updated using communication
538  // below
539  DerivedGeology distributed_geology(grid,
540  distributed_props, eclipseState,
541  useLocalPerm, geology.gravity());
542  GeologyDataHandle geo_handle(global_grid, grid,
543  geology, distributed_geology);
544  grid.scatterData(geo_handle);
545 
546  std::vector<double> distributed_pressures;
547 
548  if( !threshold_pressures.empty() ) // Might be empty if not specified
549  {
550  if( threshold_pressures.size() !=
551  static_cast<std::size_t>(UgGridHelpers::numFaces(global_grid)) )
552  {
553  OPM_THROW(std::runtime_error, "NNCs not yet supported for parallel runs. "
554  << UgGridHelpers::numFaces(grid) << " faces but " <<
555  threshold_pressures.size()<<" threshold pressure values");
556  }
557  distributed_pressures.resize(UgGridHelpers::numFaces(grid));
558  ThresholdPressureDataHandle press_handle(global_grid, grid,
559  threshold_pressures,
560  distributed_pressures);
561  grid.scatterData(press_handle);
562  }
563 
564  // copy states
565  properties = distributed_props;
566  geology = distributed_geology;
567  state = distributed_state;
568  material_law_manager = distributed_material_law_manager;
569  threshold_pressures = distributed_pressures;
570  extractParallelGridInformationToISTL(grid, parallelInformation);
571 
572  return my_defunct_wells;
573 }
574 #endif
575 
576 } // end namespace Opm
577 #endif
A handle that copies a fixed number data per index.
Definition: RedistributeDataHandles.hpp:63
void createGlobalCellArray(const Grid &grid, std::vector< int > &dest)
Create a mapping from a global cell index of a grid to the logically Cartesian index of the ECL deck...
Definition: createGlobalCellArray.hpp:31
FixedSizeIterCopyHandle(const Iter1 &send_begin, const Iter2 &receive_begin, std::size_t size=1)
Constructor.
Definition: RedistributeDataHandles.hpp:73