36#ifndef OPENRS_SIMULATORUTILITIES_HEADER
37#define OPENRS_SIMULATORUTILITIES_HEADER
40#include <opm/common/utility/platform_dependent/disable_warnings.h>
42#include <dune/common/version.hh>
43#include <dune/common/fvector.hh>
44#include <dune/grid/io/file/vtk/vtkwriter.hh>
46#include <opm/common/utility/platform_dependent/reenable_warnings.h>
48#include <opm/common/ErrorMacros.hpp>
64 template <
class Gr
idInterface,
class FlowSol>
66 const GridInterface& ginterf,
67 const FlowSol& flow_solution)
71 cell_velocity.clear();
72 cell_velocity.resize(ginterf.numberOfCells());
73 for (
typename GridInterface::CellIterator c = ginterf.cellbegin(); c != ginterf.cellend(); ++c) {
76 typename GridInterface::CellIterator::FaceIterator f = c->facebegin();
77 for (; f != c->faceend(); ++f, ++numf) {
78 double flux = flow_solution.outflux(f);
81 v *= flux/c->volume();
84 cell_velocity[c->index()] = cell_v;
94 template <
class Gr
idInterface>
96 const GridInterface& grid,
97 const std::vector<double>& face_flux)
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) {
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;
127 template <
class Gr
idInterface,
class FlowSol>
129 const GridInterface& ginterf,
130 const FlowSol& flow_solution,
131 const std::vector<int>& partition,
132 const int my_partition)
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;
144 typename GridInterface::CellIterator::FaceIterator f = c->facebegin();
145 for (; f != c->faceend(); ++f, ++numf) {
146 double flux = flow_solution.outflux(f);
149 v *= flux/c->volume();
152 cell_velocity[c->index()] = cell_v;
158 template <
class ReservoirProperty>
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)
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);
183 template <
class Gr
idInterface,
class FlowSol>
185 const GridInterface& ginterf,
186 const FlowSol& flow_solution)
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);
199 template <
class Gr
idInterface,
class FlowSol>
201 const GridInterface& ginterf,
202 const FlowSol& flow_solution,
203 const std::vector<int>& partition,
204 const int my_partition)
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;
212 cell_pressure[c->index()] = flow_solution.pressure(c);
224 template <
class ReservoirProperties>
226 const ReservoirProperties& rp,
227 const std::vector<double>& sat)
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]);
238 template <
class Gr
idInterface,
class ReservoirProperties,
class FlowSol>
240 const ReservoirProperties& rp,
241 const FlowSol& flowsol,
242 const std::vector<double>& saturation,
243 const std::string& filename)
247 std::vector<Vec> cell_velocity;
249 std::array<std::vector<Vec>, 2> phase_velocities;
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;
260 std::vector<double> cap_pressure;
262 int num_cells = saturation.size();
266 std::vector<double> fractional_flow_(num_cells);
267 for (
int i = 0; i < num_cells; ++i) {
271 fractional_flow_[i] = rp.fractionalFlow(i, saturation[i]);
275 Dune::VTKWriter<typename GridInterface::GridType::LeafGridView> vtkwriter(ginterf.grid().leafGridView());
276 vtkwriter.addCellData(saturation,
"saturation");
277 vtkwriter.addCellData(cell_pressure,
"pressure");
278 vtkwriter.addCellData(cap_pressure,
"capillary pressure");
279 vtkwriter.addCellData(fractional_flow_,
"fractional flow [water]");
282 vtkwriter.addCellData(cell_velocity_flat,
"velocity", Vec::dimension);
283 vtkwriter.addCellData(water_velocity_flat,
"phase velocity [water]", Vec::dimension);
284 vtkwriter.addCellData(oil_velocity_flat,
"phase velocity [oil]", Vec::dimension);
285 vtkwriter.write(filename, Dune::VTK::ascii);
290 const std::string& filename)
292 std::ofstream os(filename.c_str());
294 OPM_THROW(std::runtime_error,
"Could not open file " + filename);
296 os << field.size() <<
'\n';
297 std::copy(field.begin(), field.end(), std::ostream_iterator<double>(os,
"\n"));
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition: matrixops.hpp:33
Definition: ImplicitAssembly.hpp:43
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:159
void writeField(const std::vector< double > &field, const std::string &filename)
Definition: SimulatorUtilities.hpp:289
void writeVtkOutput(const GridInterface &ginterf, const ReservoirProperties &rp, const FlowSol &flowsol, const std::vector< double > &saturation, const std::string &filename)
Definition: SimulatorUtilities.hpp:239
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
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
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 getCellPressure(std::vector< double > &cell_pressure, const GridInterface &ginterf, const FlowSol &flow_solution)
Definition: SimulatorUtilities.hpp:184