36#ifndef OPENRS_SIMULATORUTILITIES_HEADER
37#define OPENRS_SIMULATORUTILITIES_HEADER
40#include <opm/common/ErrorMacros.hpp>
41#include <dune/common/version.hh>
42#include <dune/common/fvector.hh>
43#include <dune/grid/io/file/vtk/vtkwriter.hh>
59 template <
class Gr
idInterface,
class FlowSol>
61 const GridInterface& ginterf,
62 const FlowSol& flow_solution)
66 cell_velocity.clear();
67 cell_velocity.resize(ginterf.numberOfCells());
68 for (
typename GridInterface::CellIterator c = ginterf.cellbegin(); c != ginterf.cellend(); ++c) {
70 typename GridInterface::Vector cell_v(0.0);
71 typename GridInterface::CellIterator::FaceIterator f = c->facebegin();
72 for (; f != c->faceend(); ++f, ++numf) {
73 double flux = flow_solution.outflux(f);
74 typename GridInterface::Vector v = f->centroid();
76 v *= flux/c->volume();
79 cell_velocity[c->index()] = cell_v;
89 template <
class Gr
idInterface>
91 const GridInterface& grid,
92 const std::vector<double>& face_flux)
96 typedef typename GridInterface::Vector Vec;
97 cell_velocity.clear();
98 cell_velocity.resize(grid.numCells(), Vec(0.0));
99 for (
int face = 0; face < grid.numFaces(); ++face) {
100 int c[2] = { grid.faceCell(face, 0), grid.faceCell(face, 1) };
101 Vec fc = grid.faceCentroid(face);
102 double flux = face_flux[face];
103 for (
int i = 0; i < 2; ++i) {
105 Vec v_contrib = fc - grid.cellCentroid(c[i]);
106 v_contrib *= flux/grid.cellVolume(c[i]);
107 cell_velocity[c[i]] += (i == 0) ? v_contrib : -v_contrib;
122 template <
class Gr
idInterface,
class FlowSol>
124 const GridInterface& ginterf,
125 const FlowSol& flow_solution,
126 const std::vector<int>& partition,
127 const int my_partition)
131 cell_velocity.clear();
132 cell_velocity.resize(ginterf.numberOfCells());
133 for (
typename GridInterface::CellIterator c = ginterf.cellbegin(); c != ginterf.cellend(); ++c) {
134 if (partition[c->index()] != my_partition) {
135 cell_velocity[c->index()] = 0.0;
138 typename GridInterface::Vector cell_v(0.0);
139 typename GridInterface::CellIterator::FaceIterator f = c->facebegin();
140 for (; f != c->faceend(); ++f, ++numf) {
141 double flux = flow_solution.outflux(f);
142 typename GridInterface::Vector v = f->centroid();
144 v *= flux/c->volume();
147 cell_velocity[c->index()] = cell_v;
153 template <
class ReservoirProperty>
155 std::vector<Dune::FieldVector<double, 3> >& phase_velocity_oil,
156 const ReservoirProperty& res_prop,
157 const std::vector<double>& saturation,
158 const std::vector<Dune::FieldVector<double, 3> >& cell_velocity)
160 assert(saturation.size() == cell_velocity.size());
161 int num_cells = saturation.size();
162 phase_velocity_water = cell_velocity;
163 phase_velocity_oil = cell_velocity;
164 for (
int i = 0; i < num_cells; ++i) {
165 double f = res_prop.fractionalFlow(i, saturation[i]);
166 phase_velocity_water[i] *= f;
167 phase_velocity_oil[i] *= (1.0 - f);
178 template <
class Gr
idInterface,
class FlowSol>
180 const GridInterface& ginterf,
181 const FlowSol& flow_solution)
183 cell_pressure.clear();
184 cell_pressure.resize(ginterf.numberOfCells());
185 for (
typename GridInterface::CellIterator c = ginterf.cellbegin(); c != ginterf.cellend(); ++c) {
186 cell_pressure[c->index()] = flow_solution.pressure(c);
194 template <
class Gr
idInterface,
class FlowSol>
196 const GridInterface& ginterf,
197 const FlowSol& flow_solution,
198 const std::vector<int>& partition,
199 const int my_partition)
201 cell_pressure.clear();
202 cell_pressure.resize(ginterf.numberOfCells());
203 for (
typename GridInterface::CellIterator c = ginterf.cellbegin(); c != ginterf.cellend(); ++c) {
204 if (partition[c->index()] != my_partition) {
205 cell_pressure[c->index()] = 0.0;
207 cell_pressure[c->index()] = flow_solution.pressure(c);
219 template <
class ReservoirProperties>
221 const ReservoirProperties& rp,
222 const std::vector<double>& sat)
224 int num_cells = sat.size();
225 cap_pressure.resize(num_cells);
226 for (
int cell = 0; cell < num_cells; ++cell) {
227 cap_pressure[cell] = rp.capillaryPressure(cell, sat[cell]);
233 template <
class Gr
idInterface,
class ReservoirProperties,
class FlowSol>
235 const ReservoirProperties& rp,
236 const FlowSol& flowsol,
237 const std::vector<double>& saturation,
238 const std::string& filename)
241 typedef typename GridInterface::Vector Vec;
242 std::vector<Vec> cell_velocity;
244 std::array<std::vector<Vec>, 2> phase_velocities;
247 std::vector<double> cell_velocity_flat(&*cell_velocity.front().begin(),
248 &*cell_velocity.back().end());
249 std::vector<double> water_velocity_flat(&*phase_velocities[0].front().begin(),
250 &*phase_velocities[0].back().end());
251 std::vector<double> oil_velocity_flat(&*phase_velocities[1].front().begin(),
252 &*phase_velocities[1].back().end());
253 std::vector<double> cell_pressure;
255 std::vector<double> cap_pressure;
257 int num_cells = saturation.size();
261 std::vector<double> fractional_flow_(num_cells);
262 for (
int i = 0; i < num_cells; ++i) {
266 fractional_flow_[i] = rp.fractionalFlow(i, saturation[i]);
270#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
271 Dune::VTKWriter<typename GridInterface::GridType::LeafGridView> vtkwriter(ginterf.grid().leafGridView());
273 Dune::VTKWriter<typename GridInterface::GridType::LeafGridView> vtkwriter(ginterf.grid().leafView());
275 vtkwriter.addCellData(saturation,
"saturation");
276 vtkwriter.addCellData(cell_pressure,
"pressure");
277 vtkwriter.addCellData(cap_pressure,
"capillary pressure");
278 vtkwriter.addCellData(fractional_flow_,
"fractional flow [water]");
281 vtkwriter.addCellData(cell_velocity_flat,
"velocity", Vec::dimension);
282 vtkwriter.addCellData(water_velocity_flat,
"phase velocity [water]", Vec::dimension);
283 vtkwriter.addCellData(oil_velocity_flat,
"phase velocity [oil]", Vec::dimension);
284 vtkwriter.write(filename, Dune::VTK::ascii);
289 const std::string& filename)
291 std::ofstream os(filename.c_str());
293 OPM_THROW(std::runtime_error,
"Could not open file " << filename);
295 os << field.size() <<
'\n';
296 std::copy(field.begin(), field.end(), std::ostream_iterator<double>(os,
"\n"));
Definition: BlackoilFluid.hpp:32
void computePhaseVelocities(std::vector< Dune::FieldVector< double, 3 > > &phase_velocity_water, std::vector< Dune::FieldVector< double, 3 > > &phase_velocity_oil, const ReservoirProperty &res_prop, const std::vector< double > &saturation, const std::vector< Dune::FieldVector< double, 3 > > &cell_velocity)
Definition: SimulatorUtilities.hpp:154
void writeField(const std::vector< double > &field, const std::string &filename)
Definition: SimulatorUtilities.hpp:288
void writeVtkOutput(const GridInterface &ginterf, const ReservoirProperties &rp, const FlowSol &flowsol, const std::vector< double > &saturation, const std::string &filename)
Definition: SimulatorUtilities.hpp:234
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:90
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:220
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:60
void getCellPressure(std::vector< double > &cell_pressure, const GridInterface &ginterf, const FlowSol &flow_solution)
Definition: SimulatorUtilities.hpp:179