Loading...
Searching...
No Matches
SimulatorUtilities.hpp
//===========================================================================
//
// File: SimulatorUtilities.hpp
//
// Created: Fri Aug 28 15:00:15 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <bard.skaflestad@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
OpenRS is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenRS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPENRS_SIMULATORUTILITIES_HEADER
#define OPENRS_SIMULATORUTILITIES_HEADER
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <dune/common/version.hh>
#include <dune/common/fvector.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/common/ErrorMacros.hpp>
#include <vector>
#include <fstream>
#include <algorithm>
#include <iterator>
namespace Opm
{
template <class GridInterface, class FlowSol>
void estimateCellVelocity(std::vector<typename GridInterface::Vector>& cell_velocity,
const GridInterface& ginterf,
const FlowSol& flow_solution)
{
// Algorithm used is same as in halfFaceFluxToCellVelocity.hpp
// in the Sintef legacy c++ code.
cell_velocity.clear();
cell_velocity.resize(ginterf.numberOfCells());
for (typename GridInterface::CellIterator c = ginterf.cellbegin(); c != ginterf.cellend(); ++c) {
int numf = 0;
typename GridInterface::Vector cell_v(0.0);
typename GridInterface::CellIterator::FaceIterator f = c->facebegin();
for (; f != c->faceend(); ++f, ++numf) {
double flux = flow_solution.outflux(f);
typename GridInterface::Vector v = f->centroid();
v -= c->centroid();
v *= flux/c->volume();
cell_v += v;
}
cell_velocity[c->index()] = cell_v;//.two_norm();
}
}
template <class GridInterface>
void estimateCellVelocitySimpleInterface(std::vector<typename GridInterface::Vector>& cell_velocity,
const GridInterface& grid,
const std::vector<double>& face_flux)
{
// Algorithm used is same as in halfFaceFluxToCellVelocity.hpp
// in the Sintef legacy c++ code.
typedef typename GridInterface::Vector Vec;
cell_velocity.clear();
cell_velocity.resize(grid.numCells(), Vec(0.0));
for (int face = 0; face < grid.numFaces(); ++face) {
int c[2] = { grid.faceCell(face, 0), grid.faceCell(face, 1) };
Vec fc = grid.faceCentroid(face);
double flux = face_flux[face];
for (int i = 0; i < 2; ++i) {
if (c[i] >= 0) {
Vec v_contrib = fc - grid.cellCentroid(c[i]);
v_contrib *= flux/grid.cellVolume(c[i]);
cell_velocity[c[i]] += (i == 0) ? v_contrib : -v_contrib;
}
}
}
}
template <class GridInterface, class FlowSol>
void estimateCellVelocity(std::vector<typename GridInterface::Vector>& cell_velocity,
const GridInterface& ginterf,
const FlowSol& flow_solution,
const std::vector<int>& partition,
const int my_partition)
{
// Algorithm used is same as in halfFaceFluxToCellVelocity.hpp
// in the Sintef legacy c++ code.
cell_velocity.clear();
cell_velocity.resize(ginterf.numberOfCells());
for (typename GridInterface::CellIterator c = ginterf.cellbegin(); c != ginterf.cellend(); ++c) {
if (partition[c->index()] != my_partition) {
cell_velocity[c->index()] = 0.0;
} else {
int numf = 0;
typename GridInterface::Vector cell_v(0.0);
typename GridInterface::CellIterator::FaceIterator f = c->facebegin();
for (; f != c->faceend(); ++f, ++numf) {
double flux = flow_solution.outflux(f);
typename GridInterface::Vector v = f->centroid();
v -= c->centroid();
v *= flux/c->volume();
cell_v += v;
}
cell_velocity[c->index()] = cell_v;//.two_norm();
}
}
}
template <class ReservoirProperty>
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)
{
assert(saturation.size() == cell_velocity.size());
int num_cells = saturation.size();
phase_velocity_water = cell_velocity;
phase_velocity_oil = cell_velocity;
for (int i = 0; i < num_cells; ++i) {
double f = res_prop.fractionalFlow(i, saturation[i]);
phase_velocity_water[i] *= f;
phase_velocity_oil[i] *= (1.0 - f);
}
}
template <class GridInterface, class FlowSol>
void getCellPressure(std::vector<double>& cell_pressure,
const GridInterface& ginterf,
const FlowSol& flow_solution)
{
cell_pressure.clear();
cell_pressure.resize(ginterf.numberOfCells());
for (typename GridInterface::CellIterator c = ginterf.cellbegin(); c != ginterf.cellend(); ++c) {
cell_pressure[c->index()] = flow_solution.pressure(c);
}
}
template <class GridInterface, class FlowSol>
void getCellPressure(std::vector<double>& cell_pressure,
const GridInterface& ginterf,
const FlowSol& flow_solution,
const std::vector<int>& partition,
const int my_partition)
{
cell_pressure.clear();
cell_pressure.resize(ginterf.numberOfCells());
for (typename GridInterface::CellIterator c = ginterf.cellbegin(); c != ginterf.cellend(); ++c) {
if (partition[c->index()] != my_partition) {
cell_pressure[c->index()] = 0.0;
} else {
cell_pressure[c->index()] = flow_solution.pressure(c);
}
}
}
template <class ReservoirProperties>
void computeCapPressure(std::vector<double>& cap_pressure,
const ReservoirProperties& rp,
const std::vector<double>& sat)
{
int num_cells = sat.size();
cap_pressure.resize(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
cap_pressure[cell] = rp.capillaryPressure(cell, sat[cell]);
}
}
template <class GridInterface, class ReservoirProperties, class FlowSol>
void writeVtkOutput(const GridInterface& ginterf,
const ReservoirProperties& rp,
const FlowSol& flowsol,
const std::vector<double>& saturation,
const std::string& filename)
{
// Extract data in proper format.
typedef typename GridInterface::Vector Vec;
std::vector<Vec> cell_velocity;
estimateCellVelocity(cell_velocity, ginterf, flowsol);
std::array<std::vector<Vec>, 2> phase_velocities;
computePhaseVelocities(phase_velocities[0], phase_velocities[1], rp, saturation, cell_velocity);
// Dune's vtk writer wants multi-component data to be flattened.
std::vector<double> cell_velocity_flat(&*cell_velocity.front().begin(),
&*cell_velocity.back().end());
std::vector<double> water_velocity_flat(&*phase_velocities[0].front().begin(),
&*phase_velocities[0].back().end());
std::vector<double> oil_velocity_flat(&*phase_velocities[1].front().begin(),
&*phase_velocities[1].back().end());
std::vector<double> cell_pressure;
getCellPressure(cell_pressure, ginterf, flowsol);
std::vector<double> cap_pressure;
computeCapPressure(cap_pressure, rp, saturation);
int num_cells = saturation.size();
// std::array<std::vector<double>, 2> phase_mobilities_;
// phase_mobilities_[0].resize(num_cells);
// phase_mobilities_[1].resize(num_cells);
std::vector<double> fractional_flow_(num_cells);
for (int i = 0; i < num_cells; ++i) {
// for (int phase = 0; phase < 2; ++phase) {
// rp.phaseMobility(phase, i, saturation[i], phase_mobilities_[phase][i]);
// }
fractional_flow_[i] = rp.fractionalFlow(i, saturation[i]);
}
// Write data.
Dune::VTKWriter<typename GridInterface::GridType::LeafGridView> vtkwriter(ginterf.grid().leafGridView());
vtkwriter.addCellData(saturation, "saturation");
vtkwriter.addCellData(cell_pressure, "pressure");
vtkwriter.addCellData(cap_pressure, "capillary pressure");
vtkwriter.addCellData(fractional_flow_, "fractional flow [water]");
// vtkwriter.addCellData(phase_mobilities_[0], "phase mobility [water]");
// vtkwriter.addCellData(phase_mobilities_[1], "phase mobility [oil]");
vtkwriter.addCellData(cell_velocity_flat, "velocity", Vec::dimension);
vtkwriter.addCellData(water_velocity_flat, "phase velocity [water]", Vec::dimension);
vtkwriter.addCellData(oil_velocity_flat, "phase velocity [oil]", Vec::dimension);
vtkwriter.write(filename, Dune::VTK::ascii);
}
inline void writeField(const std::vector<double>& field,
const std::string& filename)
{
std::ofstream os(filename.c_str());
if (!os) {
OPM_THROW(std::runtime_error, "Could not open file " + filename);
}
os << field.size() << '\n';
std::copy(field.begin(), field.end(), std::ostream_iterator<double>(os, "\n"));
}
} // namespace Opm
#endif // OPENRS_SIMULATORUTILITIES_HEADER
Opm::GridInterface
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Opm::estimateCellVelocitySimpleInterface
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:93
Opm::computeCapPressure
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:223
Opm::estimateCellVelocity
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
Opm::getCellPressure
void getCellPressure(std::vector< double > &cell_pressure, const GridInterface &ginterf, const FlowSol &flow_solution)
Definition SimulatorUtilities.hpp:182