SimulatorUtilities.hpp
1 //===========================================================================
2 //
3 // File: SimulatorUtilities.hpp
4 //
5 // Created: Fri Aug 28 15:00:15 2009
6 //
7 // Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8 // Bård Skaflestad <bard.skaflestad@sintef.no>
9 //
10 // $Date$
11 //
12 // $Revision$
13 //
14 //===========================================================================
15 
16 /*
17  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
18  Copyright 2009, 2010 Statoil ASA.
19 
20  This file is part of The Open Reservoir Simulator Project (OpenRS).
21 
22  OpenRS is free software: you can redistribute it and/or modify
23  it under the terms of the GNU General Public License as published by
24  the Free Software Foundation, either version 3 of the License, or
25  (at your option) any later version.
26 
27  OpenRS is distributed in the hope that it will be useful,
28  but WITHOUT ANY WARRANTY; without even the implied warranty of
29  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30  GNU General Public License for more details.
31 
32  You should have received a copy of the GNU General Public License
33  along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
34 */
35 
36 #ifndef OPENRS_SIMULATORUTILITIES_HEADER
37 #define OPENRS_SIMULATORUTILITIES_HEADER
38 
39 
40 #include <opm/common/utility/platform_dependent/disable_warnings.h>
41 
42 #include <dune/common/version.hh>
43 #include <dune/common/fvector.hh>
44 #include <dune/grid/io/file/vtk/vtkwriter.hh>
45 
46 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
47 
48 #include <opm/common/ErrorMacros.hpp>
49 #include <vector>
50 #include <fstream>
51 #include <algorithm>
52 #include <iterator>
53 
54 namespace Opm
55 {
56 
57 
64  template <class GridInterface, class FlowSol>
65  void estimateCellVelocity(std::vector<typename GridInterface::Vector>& cell_velocity,
66  const GridInterface& ginterf,
67  const FlowSol& flow_solution)
68  {
69  // Algorithm used is same as in halfFaceFluxToCellVelocity.hpp
70  // in the Sintef legacy c++ code.
71  cell_velocity.clear();
72  cell_velocity.resize(ginterf.numberOfCells());
73  for (typename GridInterface::CellIterator c = ginterf.cellbegin(); c != ginterf.cellend(); ++c) {
74  int numf = 0;
75  typename GridInterface::Vector cell_v(0.0);
76  typename GridInterface::CellIterator::FaceIterator f = c->facebegin();
77  for (; f != c->faceend(); ++f, ++numf) {
78  double flux = flow_solution.outflux(f);
79  typename GridInterface::Vector v = f->centroid();
80  v -= c->centroid();
81  v *= flux/c->volume();
82  cell_v += v;
83  }
84  cell_velocity[c->index()] = cell_v;//.two_norm();
85  }
86  }
87 
94  template <class GridInterface>
95  void estimateCellVelocitySimpleInterface(std::vector<typename GridInterface::Vector>& cell_velocity,
96  const GridInterface& grid,
97  const std::vector<double>& face_flux)
98  {
99  // Algorithm used is same as in halfFaceFluxToCellVelocity.hpp
100  // in the Sintef legacy c++ code.
101  typedef typename GridInterface::Vector Vec;
102  cell_velocity.clear();
103  cell_velocity.resize(grid.numCells(), Vec(0.0));
104  for (int face = 0; face < grid.numFaces(); ++face) {
105  int c[2] = { grid.faceCell(face, 0), grid.faceCell(face, 1) };
106  Vec fc = grid.faceCentroid(face);
107  double flux = face_flux[face];
108  for (int i = 0; i < 2; ++i) {
109  if (c[i] >= 0) {
110  Vec v_contrib = fc - grid.cellCentroid(c[i]);
111  v_contrib *= flux/grid.cellVolume(c[i]);
112  cell_velocity[c[i]] += (i == 0) ? v_contrib : -v_contrib;
113  }
114  }
115  }
116  }
117 
118 
127  template <class GridInterface, class FlowSol>
128  void estimateCellVelocity(std::vector<typename GridInterface::Vector>& cell_velocity,
129  const GridInterface& ginterf,
130  const FlowSol& flow_solution,
131  const std::vector<int>& partition,
132  const int my_partition)
133  {
134  // Algorithm used is same as in halfFaceFluxToCellVelocity.hpp
135  // in the Sintef legacy c++ code.
136  cell_velocity.clear();
137  cell_velocity.resize(ginterf.numberOfCells());
138  for (typename GridInterface::CellIterator c = ginterf.cellbegin(); c != ginterf.cellend(); ++c) {
139  if (partition[c->index()] != my_partition) {
140  cell_velocity[c->index()] = 0.0;
141  } else {
142  int numf = 0;
143  typename GridInterface::Vector cell_v(0.0);
144  typename GridInterface::CellIterator::FaceIterator f = c->facebegin();
145  for (; f != c->faceend(); ++f, ++numf) {
146  double flux = flow_solution.outflux(f);
147  typename GridInterface::Vector v = f->centroid();
148  v -= c->centroid();
149  v *= flux/c->volume();
150  cell_v += v;
151  }
152  cell_velocity[c->index()] = cell_v;//.two_norm();
153  }
154  }
155  }
156 
157 
158  template <class ReservoirProperty>
159  void computePhaseVelocities(std::vector<Dune::FieldVector<double, 3> >& phase_velocity_water,
160  std::vector<Dune::FieldVector<double, 3> >& phase_velocity_oil,
161  const ReservoirProperty& res_prop,
162  const std::vector<double>& saturation,
163  const std::vector<Dune::FieldVector<double, 3> >& cell_velocity)
164  {
165  assert(saturation.size() == cell_velocity.size());
166  int num_cells = saturation.size();
167  phase_velocity_water = cell_velocity;
168  phase_velocity_oil = cell_velocity;
169  for (int i = 0; i < num_cells; ++i) {
170  double f = res_prop.fractionalFlow(i, saturation[i]);
171  phase_velocity_water[i] *= f;
172  phase_velocity_oil[i] *= (1.0 - f);
173  }
174  }
175 
176 
177 
178 
183  template <class GridInterface, class FlowSol>
184  void getCellPressure(std::vector<double>& cell_pressure,
185  const GridInterface& ginterf,
186  const FlowSol& flow_solution)
187  {
188  cell_pressure.clear();
189  cell_pressure.resize(ginterf.numberOfCells());
190  for (typename GridInterface::CellIterator c = ginterf.cellbegin(); c != ginterf.cellend(); ++c) {
191  cell_pressure[c->index()] = flow_solution.pressure(c);
192  }
193  }
194 
199  template <class GridInterface, class FlowSol>
200  void getCellPressure(std::vector<double>& cell_pressure,
201  const GridInterface& ginterf,
202  const FlowSol& flow_solution,
203  const std::vector<int>& partition,
204  const int my_partition)
205  {
206  cell_pressure.clear();
207  cell_pressure.resize(ginterf.numberOfCells());
208  for (typename GridInterface::CellIterator c = ginterf.cellbegin(); c != ginterf.cellend(); ++c) {
209  if (partition[c->index()] != my_partition) {
210  cell_pressure[c->index()] = 0.0;
211  } else {
212  cell_pressure[c->index()] = flow_solution.pressure(c);
213  }
214  }
215  }
216 
217 
218 
224  template <class ReservoirProperties>
225  void computeCapPressure(std::vector<double>& cap_pressure,
226  const ReservoirProperties& rp,
227  const std::vector<double>& sat)
228  {
229  int num_cells = sat.size();
230  cap_pressure.resize(num_cells);
231  for (int cell = 0; cell < num_cells; ++cell) {
232  cap_pressure[cell] = rp.capillaryPressure(cell, sat[cell]);
233  }
234  }
235 
236 
238  template <class GridInterface, class ReservoirProperties, class FlowSol>
239  void writeVtkOutput(const GridInterface& ginterf,
240  const ReservoirProperties& rp,
241  const FlowSol& flowsol,
242  const std::vector<double>& saturation,
243  const std::string& filename)
244  {
245  // Extract data in proper format.
246  typedef typename GridInterface::Vector Vec;
247  std::vector<Vec> cell_velocity;
248  estimateCellVelocity(cell_velocity, ginterf, flowsol);
249  std::array<std::vector<Vec>, 2> phase_velocities;
250  computePhaseVelocities(phase_velocities[0], phase_velocities[1], rp, saturation, cell_velocity);
251  // Dune's vtk writer wants multi-component data to be flattened.
252  std::vector<double> cell_velocity_flat(&*cell_velocity.front().begin(),
253  &*cell_velocity.back().end());
254  std::vector<double> water_velocity_flat(&*phase_velocities[0].front().begin(),
255  &*phase_velocities[0].back().end());
256  std::vector<double> oil_velocity_flat(&*phase_velocities[1].front().begin(),
257  &*phase_velocities[1].back().end());
258  std::vector<double> cell_pressure;
259  getCellPressure(cell_pressure, ginterf, flowsol);
260  std::vector<double> cap_pressure;
261  computeCapPressure(cap_pressure, rp, saturation);
262  int num_cells = saturation.size();
263 // std::array<std::vector<double>, 2> phase_mobilities_;
264 // phase_mobilities_[0].resize(num_cells);
265 // phase_mobilities_[1].resize(num_cells);
266  std::vector<double> fractional_flow_(num_cells);
267  for (int i = 0; i < num_cells; ++i) {
268 // for (int phase = 0; phase < 2; ++phase) {
269 // rp.phaseMobility(phase, i, saturation[i], phase_mobilities_[phase][i]);
270 // }
271  fractional_flow_[i] = rp.fractionalFlow(i, saturation[i]);
272  }
273 
274  // Write data.
275 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
276  Dune::VTKWriter<typename GridInterface::GridType::LeafGridView> vtkwriter(ginterf.grid().leafGridView());
277 #else
278  Dune::VTKWriter<typename GridInterface::GridType::LeafGridView> vtkwriter(ginterf.grid().leafView());
279 #endif
280  vtkwriter.addCellData(saturation, "saturation");
281  vtkwriter.addCellData(cell_pressure, "pressure");
282  vtkwriter.addCellData(cap_pressure, "capillary pressure");
283  vtkwriter.addCellData(fractional_flow_, "fractional flow [water]");
284 // vtkwriter.addCellData(phase_mobilities_[0], "phase mobility [water]");
285 // vtkwriter.addCellData(phase_mobilities_[1], "phase mobility [oil]");
286  vtkwriter.addCellData(cell_velocity_flat, "velocity", Vec::dimension);
287  vtkwriter.addCellData(water_velocity_flat, "phase velocity [water]", Vec::dimension);
288  vtkwriter.addCellData(oil_velocity_flat, "phase velocity [oil]", Vec::dimension);
289  vtkwriter.write(filename, Dune::VTK::ascii);
290  }
291 
292 
293  inline void writeField(const std::vector<double>& field,
294  const std::string& filename)
295  {
296  std::ofstream os(filename.c_str());
297  if (!os) {
298  OPM_THROW(std::runtime_error, "Could not open file " << filename);
299  }
300  os << field.size() << '\n';
301  std::copy(field.begin(), field.end(), std::ostream_iterator<double>(os, "\n"));
302  }
303 
304 } // namespace Opm
305 
306 
307 #endif // OPENRS_SIMULATORUTILITIES_HEADER
void estimateCellVelocity(std::vector< typename GridInterface::Vector > &cell_velocity, const GridInterface &ginterf, const FlowSol &flow_solution)
Estimates a scalar cell velocity from outgoing fluxes.
Definition: SimulatorUtilities.hpp:65
void computeCapPressure(std::vector< double > &cap_pressure, const ReservoirProperties &rp, const std::vector< double > &sat)
Computes the capillary pressure in each cell from the cell saturations.
Definition: SimulatorUtilities.hpp:225
Class for immiscible dead oil and dry gas.
Definition: applier.hpp:18
void getCellPressure(std::vector< double > &cell_pressure, const GridInterface &ginterf, const FlowSol &flow_solution)
Definition: SimulatorUtilities.hpp:184
void estimateCellVelocitySimpleInterface(std::vector< typename GridInterface::Vector > &cell_velocity, const GridInterface &grid, const std::vector< double > &face_flux)
Estimates a scalar cell velocity from face fluxes.
Definition: SimulatorUtilities.hpp:95