BlackoilDetails.hpp
1 /*
2  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
3  Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services
4  Copyright 2014, 2015 Statoil ASA.
5  Copyright 2015 NTNU
6  Copyright 2015 IRIS AS
7 
8  This file is part of the Open Porous Media project (OPM).
9 
10  OPM is free software: you can redistribute it and/or modify
11  it under the terms of the GNU General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  OPM is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  GNU General Public License for more details.
19 
20  You should have received a copy of the GNU General Public License
21  along with OPM. If not, see <http://www.gnu.org/licenses/>.
22 */
23 
24 #ifndef OPM_BLACKOILDETAILS_HEADER_INCLUDED
25 #define OPM_BLACKOILDETAILS_HEADER_INCLUDED
26 
27 #include <opm/core/linalg/ParallelIstlInformation.hpp>
28 
29 namespace Opm {
30 namespace detail {
31 
32 
33  inline
34  std::vector<int>
35  buildAllCells(const int nc)
36  {
37  std::vector<int> all_cells(nc);
38 
39  for (int c = 0; c < nc; ++c) { all_cells[c] = c; }
40 
41  return all_cells;
42  }
43 
44 
45 
46  template <class PU>
47  std::vector<bool>
48  activePhases(const PU& pu)
49  {
50  const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
51  std::vector<bool> active(maxnp, false);
52 
53  for (int p = 0; p < pu.MaxNumPhases; ++p) {
54  active[ p ] = pu.phase_used[ p ] != 0;
55  }
56 
57  return active;
58  }
59 
60 
61 
62  template <class PU>
63  std::vector<int>
64  active2Canonical(const PU& pu)
65  {
66  const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
67  std::vector<int> act2can(maxnp, -1);
68 
69  for (int phase = 0; phase < maxnp; ++phase) {
70  if (pu.phase_used[ phase ]) {
71  act2can[ pu.phase_pos[ phase ] ] = phase;
72  }
73  }
74 
75  return act2can;
76  }
77 
78 
79 
80  inline
81  double getGravity(const double* g, const int dim) {
82  double grav = 0.0;
83  if (g) {
84  // Guard against gravity in anything but last dimension.
85  for (int dd = 0; dd < dim - 1; ++dd) {
86  assert(g[dd] == 0.0);
87  }
88  grav = g[dim - 1];
89  }
90  return grav;
91  }
92 
93 
103  template <class Iterator>
104  inline
105  double euclidianNormSquared( Iterator it, const Iterator end, int num_components, const boost::any& pinfo = boost::any() )
106  {
107  static_cast<void>(num_components); // Suppress warning in the serial case.
108  static_cast<void>(pinfo); // Suppress warning in non-MPI case.
109 #if HAVE_MPI
110  if ( pinfo.type() == typeid(ParallelISTLInformation) )
111  {
112  const ParallelISTLInformation& info =
113  boost::any_cast<const ParallelISTLInformation&>(pinfo);
114  typedef typename Iterator::value_type Scalar;
115  Scalar product = 0.0;
116  int size_per_component = (end - it);
117  size_per_component /= num_components; // two lines to supresse unused warning.
118  assert((end - it) == num_components * size_per_component);
119 
120  if( num_components == 1 )
121  {
122  auto component_container =
123  boost::make_iterator_range(it, end);
124  info.computeReduction(component_container,
125  Opm::Reduction::makeInnerProductFunctor<double>(),
126  product);
127  }
128  else
129  {
130  auto& maskContainer = info.getOwnerMask();
131  auto mask = maskContainer.begin();
132  assert(static_cast<int>(maskContainer.size()) == size_per_component);
133 
134  for(int cell = 0; cell < size_per_component; ++cell, ++mask)
135  {
136  Scalar cell_product = (*it) * (*it);
137  ++it;
138  for(int component=1; component < num_components;
139  ++component, ++it)
140  {
141  cell_product += (*it) * (*it);
142  }
143  product += cell_product * (*mask);
144  }
145  }
146  return info.communicator().sum(product);
147  }
148  else
149 #endif
150  {
151  double product = 0.0 ;
152  for( ; it != end; ++it ) {
153  product += ( *it * *it );
154  }
155  return product;
156  }
157  }
163  template<class Grid>
164  std::size_t countLocalInteriorCells(const Grid& grid)
165  {
166  if ( grid.comm().size() == 1)
167  {
168  return grid.size(0);
169  }
170  std::size_t count = 0;
171  const auto& gridView = grid.leafGridView();
172  for(auto cell = gridView.template begin<0, Dune::Interior_Partition>(),
173  endCell = gridView.template end<0, Dune::Interior_Partition>();
174  cell != endCell; ++cell)
175  {
176  ++count;
177  }
178  return count;
179  }
180 
188  template<class Grid>
189  std::size_t countGlobalCells(const Grid& grid)
190  {
191  if ( grid.comm().size() == 1)
192  {
193  return grid.size(0);
194  }
195  std::size_t count = countLocalInteriorCells(grid);
196  return grid.comm().sum(count);
197  }
198 
199 
200  } // namespace detail
201 } // namespace Opm
202 
203 #endif // OPM_BLACKOILDETAILS_HEADER_INCLUDED
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22