OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Reorganized and cleaned up in prep for polymer work.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2012-02-02 16:03:50 +01:00
parent 1c66958c0a
commit 11b4c91d2e
4 changed files with 671 additions and 635 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
examples/Makefile.am
View File

@ -9,4 +9,6 @@ hello_SOURCES = \
SayHello.cpp
polymer_reorder_SOURCES = \
polymer_reorder.cpp
polymer_reorder.cpp \
Utilities.hpp \
Utilities.cpp

418
examples/Utilities.cpp Normal file
View File

@ -0,0 +1,418 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include "Utilities.hpp"
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/utility/cart_grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/cpgpreprocess/cgridinterface.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/SimpleFluid2p.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
#include <opm/core/transport/reorder/twophasetransport.hpp>
#include <boost/filesystem/convenience.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/lexical_cast.hpp>
#include <cassert>
#include <cstddef>
#include <algorithm>
#include <tr1/array>
#include <functional>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <iterator>
#include <vector>
#include <numeric>
namespace Opm
{
void
compute_porevolume(const UnstructuredGrid* g,
const Opm::IncompPropertiesInterface& props,
std::vector<double>& porevol)
{
int num_cells = g->number_of_cells;
porevol.resize(num_cells);
const double* poro = props.porosity();
::std::transform(poro, poro + num_cells,
g->cell_volumes,
porevol.begin(),
::std::multiplies<double>());
}
void
compute_totmob(const Opm::IncompPropertiesInterface& props,
const std::vector<double>& s,
std::vector<double>& totmob)
{
int num_cells = props.numCells();
int num_phases = props.numPhases();
totmob.resize(num_cells);
ASSERT(int(s.size()) == num_cells*num_phases);
std::vector<int> cells(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
cells[cell] = cell;
}
std::vector<double> kr(num_cells*num_phases);
props.relperm(num_cells, &s[0], &cells[0], &kr[0], 0);
const double* mu = props.viscosity();
for (int cell = 0; cell < num_cells; ++cell) {
totmob[cell] = 0;
for (int phase = 0; phase < num_phases; ++phase) {
totmob[cell] += kr[2*cell + phase]/mu[phase];
}
}
}
void writeVtkDataAllCartesian(const std::tr1::array<int, 3>& dims,
const std::tr1::array<double, 3>& cell_size,
const std::vector<double>& pressure,
const std::vector<double>& saturation,
std::ostream& vtk_file)
{
// Dimension is hardcoded in the prototype and the next two lines,
// but the rest is flexible (allows dimension == 2 or 3).
int dimension = 3;
int num_cells = dims[0]*dims[1]*dims[2];
ASSERT(dimension == 2 || dimension == 3);
ASSERT(num_cells = dims[0]*dims[1]* (dimension == 2 ? 1 : dims[2]));
vtk_file << "# vtk DataFile Version 2.0\n";
vtk_file << "Structured Grid\n \n";
vtk_file << "ASCII \n";
vtk_file << "DATASET STRUCTURED_POINTS\n";
vtk_file << "DIMENSIONS "
<< dims[0] + 1 << " "
<< dims[1] + 1 << " ";
if (dimension == 3) {
vtk_file << dims[2] + 1;
} else {
vtk_file << 1;
}
vtk_file << "\n";
vtk_file << "ORIGIN " << 0.0 << " " << 0.0 << " " << 0.0 << "\n";
vtk_file << "SPACING " << cell_size[0] << " " << cell_size[1];
if (dimension == 3) {
vtk_file << " " << cell_size[2];
} else {
vtk_file << " " << 0.0;
}
vtk_file << "\n";
vtk_file << "CELL_DATA " << num_cells << '\n';
vtk_file << "SCALARS pressure float" << '\n';
vtk_file << "LOOKUP_TABLE pressure_table " << '\n';
for (int i = 0; i < num_cells; ++i) {
vtk_file << pressure[i] << '\n';
}
vtk_file << "SCALARS saturation float" << '\n';
vtk_file << "LOOKUP_TABLE saturation_table " << '\n';
for (int i = 0; i < num_cells; ++i) {
double s = saturation[2*i];
if (s > 1e-10) {
vtk_file << s << '\n';
} else {
vtk_file << 0.0 << '\n';
}
}
}
typedef std::map<std::string, std::string> PMap;
struct Tag
{
Tag(const std::string& tag, const PMap& props, std::ostream& os)
: name_(tag), os_(os)
{
indent(os);
os << "<" << tag;
for (PMap::const_iterator it = props.begin(); it != props.end(); ++it) {
os << " " << it->first << "=\"" << it->second << "\"";
}
os << ">\n";
++indent_;
}
Tag(const std::string& tag, std::ostream& os)
: name_(tag), os_(os)
{
indent(os);
os << "<" << tag << ">\n";
++indent_;
}
~Tag()
{
--indent_;
indent(os_);
os_ << "</" << name_ << ">\n";
}
static void indent(std::ostream& os)
{
for (int i = 0; i < indent_; ++i) {
os << " ";
}
}
private:
static int indent_;
std::string name_;
std::ostream& os_;
};
int Tag::indent_ = 0;
void writeVtkDataGeneralGrid(const UnstructuredGrid* grid,
const std::vector<double>& pressure,
const std::vector<double>& saturation,
std::ostream& os)
{
if (grid->dimensions != 3) {
THROW("Vtk output for 3d grids only");
}
os.precision(12);
os << "<?xml version=\"1.0\"?>\n";
PMap pm;
pm["type"] = "UnstructuredGrid";
Tag vtkfiletag("VTKFile", pm, os);
Tag ugtag("UnstructuredGrid", os);
int num_pts = grid->number_of_nodes;
int num_cells = grid->number_of_cells;
pm.clear();
pm["NumberOfPoints"] = boost::lexical_cast<std::string>(num_pts);
pm["NumberOfCells"] = boost::lexical_cast<std::string>(num_cells);
Tag piecetag("Piece", pm, os);
{
Tag pointstag("Points", os);
pm.clear();
pm["type"] = "Float64";
pm["Name"] = "Coordinates";
pm["NumberOfComponents"] = "3";
pm["format"] = "ascii";
Tag datag("DataArray", pm, os);
for (int i = 0; i < num_pts; ++i) {
Tag::indent(os);
os << grid->node_coordinates[3*i + 0] << ' '
<< grid->node_coordinates[3*i + 1] << ' '
<< grid->node_coordinates[3*i + 2] << '\n';
}
}
{
Tag cellstag("Cells", os);
pm.clear();
pm["type"] = "Int32";
pm["NumberOfComponents"] = "1";
pm["format"] = "ascii";
std::vector<int> cell_numpts;
cell_numpts.reserve(num_cells);
{
pm["Name"] = "connectivity";
Tag t("DataArray", pm, os);
int hf = 0;
for (int c = 0; c < num_cells; ++c) {
std::set<int> cell_pts;
for (; hf < grid->cell_facepos[c+1]; ++hf) {
int f = grid->cell_faces[hf];
const int* fnbeg = grid->face_nodes + grid->face_nodepos[f];
const int* fnend = grid->face_nodes + grid->face_nodepos[f+1];
cell_pts.insert(fnbeg, fnend);
}
cell_numpts.push_back(cell_pts.size());
Tag::indent(os);
std::copy(cell_pts.begin(), cell_pts.end(),
std::ostream_iterator<int>(os, " "));
os << '\n';
}
}
{
pm["Name"] = "offsets";
Tag t("DataArray", pm, os);
int offset = 0;
const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) {
Tag::indent(os);
}
offset += cell_numpts[c];
os << offset << ' ';
if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) {
os << '\n';
}
}
}
std::vector<int> cell_foffsets;
cell_foffsets.reserve(num_cells);
{
pm["Name"] = "faces";
Tag t("DataArray", pm, os);
const int* fp = grid->cell_facepos;
int offset = 0;
for (int c = 0; c < num_cells; ++c) {
Tag::indent(os);
os << fp[c+1] - fp[c] << '\n';
++offset;
for (int hf = fp[c]; hf < fp[c+1]; ++hf) {
int f = grid->cell_faces[hf];
const int* np = grid->face_nodepos;
int f_num_pts = np[f+1] - np[f];
Tag::indent(os);
os << f_num_pts << ' ';
++offset;
std::copy(grid->face_nodes + np[f],
grid->face_nodes + np[f+1],
std::ostream_iterator<int>(os, " "));
os << '\n';
offset += f_num_pts;
}
cell_foffsets.push_back(offset);
}
}
{
pm["Name"] = "faceoffsets";
Tag t("DataArray", pm, os);
const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) {
Tag::indent(os);
}
os << cell_foffsets[c] << ' ';
if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) {
os << '\n';
}
}
}
{
pm["type"] = "UInt8";
pm["Name"] = "types";
Tag t("DataArray", pm, os);
const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) {
Tag::indent(os);
}
os << "42 ";
if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) {
os << '\n';
}
}
}
}
{
pm.clear();
pm["Scalars"] = "saturation";
Tag celldatatag("CellData", pm, os);
pm.clear();
pm["type"] = "Int32";
pm["NumberOfComponents"] = "1";
pm["format"] = "ascii";
pm["type"] = "Float64";
{
pm["Name"] = "pressure";
Tag ptag("DataArray", pm, os);
const int num_per_line = 5;
for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) {
Tag::indent(os);
}
os << pressure[c] << ' ';
if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) {
os << '\n';
}
}
}
{
pm["Name"] = "saturation";
Tag ptag("DataArray", pm, os);
const int num_per_line = 5;
for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) {
Tag::indent(os);
}
os << saturation[2*c] << ' ';
if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) {
os << '\n';
}
}
}
}
}
void toWaterSat(const std::vector<double>& sboth, std::vector<double>& sw)
{
int num = sboth.size()/2;
sw.resize(num);
for (int i = 0; i < num; ++i) {
sw[i] = sboth[2*i];
}
}
void toBothSat(const std::vector<double>& sw, std::vector<double>& sboth)
{
int num = sw.size();
sboth.resize(2*num);
for (int i = 0; i < num; ++i) {
sboth[2*i] = sw[i];
sboth[2*i + 1] = 1.0 - sw[i];
}
}
} // namespace Opm

226
examples/Utilities.hpp Normal file
View File

@ -0,0 +1,226 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_UTILITIES_HEADER_INCLUDED
#define OPM_UTILITIES_HEADER_INCLUDED
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/utility/cart_grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/cpgpreprocess/cgridinterface.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/SimpleFluid2p.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
#include <opm/core/transport/reorder/twophasetransport.hpp>
#include <boost/filesystem/convenience.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/lexical_cast.hpp>
#include <cassert>
#include <cstddef>
#include <algorithm>
#include <tr1/array>
#include <functional>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <iterator>
#include <vector>
#include <numeric>
namespace Opm
{
/// Concrete grid class constructing a
/// corner point grid from a deck,
/// or a cartesian grid.
class Grid
{
public:
Grid(const Opm::EclipseGridParser& deck)
{
// Extract data from deck.
const std::vector<double>& zcorn = deck.getFloatingPointValue("ZCORN");
const std::vector<double>& coord = deck.getFloatingPointValue("COORD");
const std::vector<int>& actnum = deck.getIntegerValue("ACTNUM");
std::vector<int> dims;
if (deck.hasField("DIMENS")) {
dims = deck.getIntegerValue("DIMENS");
} else if (deck.hasField("SPECGRID")) {
dims = deck.getSPECGRID().dimensions;
} else {
THROW("Deck must have either DIMENS or SPECGRID.");
}
// Collect in input struct for preprocessing.
struct grdecl grdecl;
grdecl.zcorn = &zcorn[0];
grdecl.coord = &coord[0];
grdecl.actnum = &actnum[0];
grdecl.dims[0] = dims[0];
grdecl.dims[1] = dims[1];
grdecl.dims[2] = dims[2];
// Process and compute.
ug_ = preprocess(&grdecl, 0.0);
compute_geometry(ug_);
}
Grid(int nx, int ny)
{
ug_ = create_cart_grid_2d(nx, ny);
}
Grid(int nx, int ny, int nz)
{
ug_ = create_cart_grid_3d(nx, ny, nz);
}
~Grid()
{
free_grid(ug_);
}
virtual const UnstructuredGrid* c_grid() const
{
return ug_;
}
private:
// Disable copying and assignment.
Grid(const Grid& other);
Grid& operator=(const Grid& other);
struct UnstructuredGrid* ug_;
};
class PressureSolver
{
public:
PressureSolver(const UnstructuredGrid* g,
const IncompPropertiesInterface& props)
: htrans_(g->cell_facepos[ g->number_of_cells ]),
trans_ (g->number_of_faces),
gpress_(g->cell_facepos[ g->number_of_cells ])
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(g);
tpfa_htrans_compute(gg, props.permeability(), &htrans_[0]);
h_ = ifs_tpfa_construct(gg);
}
~PressureSolver()
{
ifs_tpfa_destroy(h_);
}
template <class State>
void
solve(const UnstructuredGrid* g ,
const ::std::vector<double>& totmob,
const ::std::vector<double>& src ,
State& state )
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(g);
tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
// No gravity
std::fill(gpress_.begin(), gpress_.end(), double(0.0));
ifs_tpfa_assemble(gg, &trans_[0], &src[0], &gpress_[0], h_);
using ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver;
CSRMatrixUmfpackSolver linsolve;
linsolve.solve(h_->A, h_->b, h_->x);
ifs_tpfa_press_flux(gg, &trans_[0], h_,
&state.pressure()[0],
&state.faceflux()[0]);
}
private:
::std::vector<double> htrans_;
::std::vector<double> trans_ ;
::std::vector<double> gpress_;
struct ifs_tpfa_data* h_;
};
void
compute_porevolume(const UnstructuredGrid* g,
const Opm::IncompPropertiesInterface& props,
std::vector<double>& porevol);
void
compute_totmob(const Opm::IncompPropertiesInterface& props,
const std::vector<double>& s,
std::vector<double>& totmob);
void writeVtkDataAllCartesian(const std::tr1::array<int, 3>& dims,
const std::tr1::array<double, 3>& cell_size,
const std::vector<double>& pressure,
const std::vector<double>& saturation,
std::ostream& vtk_file);
void writeVtkDataGeneralGrid(const UnstructuredGrid* grid,
const std::vector<double>& pressure,
const std::vector<double>& saturation,
std::ostream& os);
void toWaterSat(const std::vector<double>& sboth, std::vector<double>& sw);
void toBothSat(const std::vector<double>& sw, std::vector<double>& sboth);
} // namespace Opm
#endif // OPM_UTILITIES_HEADER_INCLUDED

658
examples/polymer_reorder.cpp
View File

@ -1,21 +1,6 @@
/*===========================================================================
//
// File: spu_2p.cpp
//
// Created: 2011-10-05 10:29:01+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
Copyright 2011 Statoil ASA.
Copyright 2012 SINTEF ICT, Applied Mathematics.
Copyright 2012 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
@ -32,9 +17,10 @@
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/linalg/sparse_sys.h>
#include "Utilities.hpp"
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
@ -49,14 +35,7 @@
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/transport/transport_source.h>
#include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
#include <opm/core/transport/NormSupport.hpp>
#include <opm/core/transport/ImplicitAssembly.hpp>
#include <opm/core/transport/ImplicitTransport.hpp>
#include <opm/core/transport/JacobianSystem.hpp>
#include <opm/core/transport/CSRMatrixBlockAssembler.hpp>
#include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
#include <opm/core/transport/reorder/twophasetransport.hpp>
@ -75,79 +54,11 @@
#include <fstream>
#include <iterator>
#include <vector>
#include <numeric>
namespace Opm
{
/// Concrete grid class constructing a
/// corner point grid from a deck,
/// or a cartesian grid.
class Grid
{
public:
Grid(const Opm::EclipseGridParser& deck)
{
// Extract data from deck.
const std::vector<double>& zcorn = deck.getFloatingPointValue("ZCORN");
const std::vector<double>& coord = deck.getFloatingPointValue("COORD");
const std::vector<int>& actnum = deck.getIntegerValue("ACTNUM");
std::vector<int> dims;
if (deck.hasField("DIMENS")) {
dims = deck.getIntegerValue("DIMENS");
} else if (deck.hasField("SPECGRID")) {
dims = deck.getSPECGRID().dimensions;
} else {
THROW("Deck must have either DIMENS or SPECGRID.");
}
// Collect in input struct for preprocessing.
struct grdecl grdecl;
grdecl.zcorn = &zcorn[0];
grdecl.coord = &coord[0];
grdecl.actnum = &actnum[0];
grdecl.dims[0] = dims[0];
grdecl.dims[1] = dims[1];
grdecl.dims[2] = dims[2];
// Process and compute.
ug_ = preprocess(&grdecl, 0.0);
compute_geometry(ug_);
}
Grid(int nx, int ny)
{
ug_ = create_cart_grid_2d(nx, ny);
}
Grid(int nx, int ny, int nz)
{
ug_ = create_cart_grid_3d(nx, ny, nz);
}
~Grid()
{
free_grid(ug_);
}
virtual const UnstructuredGrid* c_grid() const
{
return ug_;
}
private:
// Disable copying and assignment.
Grid(const Grid& other);
Grid& operator=(const Grid& other);
struct UnstructuredGrid* ug_;
};
} // namespace Opm
class ReservoirState {
@ -185,100 +96,6 @@ private:
class PressureSolver {
public:
PressureSolver(const UnstructuredGrid* g,
const Opm::IncompPropertiesInterface& props)
: htrans_(g->cell_facepos[ g->number_of_cells ]),
trans_ (g->number_of_faces),
gpress_(g->cell_facepos[ g->number_of_cells ])
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(g);
tpfa_htrans_compute(gg, props.permeability(), &htrans_[0]);
h_ = ifs_tpfa_construct(gg);
}
~PressureSolver()
{
ifs_tpfa_destroy(h_);
}
template <class State>
void
solve(const UnstructuredGrid* g ,
const ::std::vector<double>& totmob,
const ::std::vector<double>& src ,
State& state )
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(g);
tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
// No gravity
::std::fill(gpress_.begin(), gpress_.end(), double(0.0));
ifs_tpfa_assemble(gg, &trans_[0], &src[0], &gpress_[0], h_);
using Opm::ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver;
CSRMatrixUmfpackSolver linsolve;
linsolve.solve(h_->A, h_->b, h_->x);
ifs_tpfa_press_flux(gg, &trans_[0], h_,
&state.pressure()[0],
&state.faceflux()[0]);
}
private:
::std::vector<double> htrans_;
::std::vector<double> trans_ ;
::std::vector<double> gpress_;
struct ifs_tpfa_data* h_;
};
static void
compute_porevolume(const UnstructuredGrid* g,
const Opm::IncompPropertiesInterface& props,
std::vector<double>& porevol)
{
int num_cells = g->number_of_cells;
porevol.resize(num_cells);
const double* poro = props.porosity();
::std::transform(poro, poro + num_cells,
g->cell_volumes,
porevol.begin(),
::std::multiplies<double>());
}
static void
compute_totmob(const Opm::IncompPropertiesInterface& props,
const std::vector<double>& s,
std::vector<double>& totmob)
{
int num_cells = props.numCells();
int num_phases = props.numPhases();
totmob.resize(num_cells);
ASSERT(int(s.size()) == num_cells*num_phases);
std::vector<int> cells(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
cells[cell] = cell;
}
std::vector<double> kr(num_cells*num_phases);
props.relperm(num_cells, &s[0], &cells[0], &kr[0], 0);
const double* mu = props.viscosity();
for (int cell = 0; cell < num_cells; ++cell) {
totmob[cell] = 0;
for (int phase = 0; phase < num_phases; ++phase) {
totmob[cell] += kr[2*cell + phase]/mu[phase];
}
}
}
template <class State>
void outputState(const UnstructuredGrid* grid,
@ -292,409 +109,12 @@ void outputState(const UnstructuredGrid* grid,
if (!vtkfile) {
THROW("Failed to open " << vtkfilename.str());
}
writeVtkDataGeneralGrid(grid, state, vtkfile);
Opm::writeVtkDataGeneralGrid(grid, state.pressure(), state.saturation(), vtkfile);
}
template <class State>
void writeVtkDataAllCartesian(const std::tr1::array<int, 3>& dims,
const std::tr1::array<double, 3>& cell_size,
const State& state,
std::ostream& vtk_file)
{
// Dimension is hardcoded in the prototype and the next two lines,
// but the rest is flexible (allows dimension == 2 or 3).
int dimension = 3;
int num_cells = dims[0]*dims[1]*dims[2];
ASSERT(dimension == 2 || dimension == 3);
ASSERT(num_cells = dims[0]*dims[1]* (dimension == 2 ? 1 : dims[2]));
vtk_file << "# vtk DataFile Version 2.0\n";
vtk_file << "Structured Grid\n \n";
vtk_file << "ASCII \n";
vtk_file << "DATASET STRUCTURED_POINTS\n";
vtk_file << "DIMENSIONS "
<< dims[0] + 1 << " "
<< dims[1] + 1 << " ";
if (dimension == 3) {
vtk_file << dims[2] + 1;
} else {
vtk_file << 1;
}
vtk_file << "\n";
vtk_file << "ORIGIN " << 0.0 << " " << 0.0 << " " << 0.0 << "\n";
vtk_file << "SPACING " << cell_size[0] << " " << cell_size[1];
if (dimension == 3) {
vtk_file << " " << cell_size[2];
} else {
vtk_file << " " << 0.0;
}
vtk_file << "\n";
vtk_file << "CELL_DATA " << num_cells << '\n';
vtk_file << "SCALARS pressure float" << '\n';
vtk_file << "LOOKUP_TABLE pressure_table " << '\n';
for (int i = 0; i < num_cells; ++i) {
vtk_file << state.pressure()[i] << '\n';
}
ASSERT(state.numPhases() == 2);
vtk_file << "SCALARS saturation float" << '\n';
vtk_file << "LOOKUP_TABLE saturation_table " << '\n';
for (int i = 0; i < num_cells; ++i) {
double s = state.saturation()[2*i];
if (s > 1e-10) {
vtk_file << s << '\n';
} else {
vtk_file << 0.0 << '\n';
}
}
}
typedef std::map<std::string, std::string> PMap;
struct Tag
{
Tag(const std::string& tag, const PMap& props, std::ostream& os)
: name_(tag), os_(os)
{
indent(os);
os << "<" << tag;
for (PMap::const_iterator it = props.begin(); it != props.end(); ++it) {
os << " " << it->first << "=\"" << it->second << "\"";
}
os << ">\n";
++indent_;
}
Tag(const std::string& tag, std::ostream& os)
: name_(tag), os_(os)
{
indent(os);
os << "<" << tag << ">\n";
++indent_;
}
~Tag()
{
--indent_;
indent(os_);
os_ << "</" << name_ << ">\n";
}
static void indent(std::ostream& os)
{
for (int i = 0; i < indent_; ++i) {
os << " ";
}
}
private:
static int indent_;
std::string name_;
std::ostream& os_;
};
int Tag::indent_ = 0;
template <class State>
void writeVtkDataGeneralGrid(const UnstructuredGrid* grid,
const State& state,
std::ostream& os)
{
if (grid->dimensions != 3) {
THROW("Vtk output for 3d grids only");
}
os.precision(12);
os << "<?xml version=\"1.0\"?>\n";
PMap pm;
pm["type"] = "UnstructuredGrid";
Tag vtkfiletag("VTKFile", pm, os);
Tag ugtag("UnstructuredGrid", os);
int num_pts = grid->number_of_nodes;
int num_cells = grid->number_of_cells;
pm.clear();
pm["NumberOfPoints"] = boost::lexical_cast<std::string>(num_pts);
pm["NumberOfCells"] = boost::lexical_cast<std::string>(num_cells);
Tag piecetag("Piece", pm, os);
{
Tag pointstag("Points", os);
pm.clear();
pm["type"] = "Float64";
pm["Name"] = "Coordinates";
pm["NumberOfComponents"] = "3";
pm["format"] = "ascii";
Tag datag("DataArray", pm, os);
for (int i = 0; i < num_pts; ++i) {
Tag::indent(os);
os << grid->node_coordinates[3*i + 0] << ' '
<< grid->node_coordinates[3*i + 1] << ' '
<< grid->node_coordinates[3*i + 2] << '\n';
}
}
{
Tag cellstag("Cells", os);
pm.clear();
pm["type"] = "Int32";
pm["NumberOfComponents"] = "1";
pm["format"] = "ascii";
std::vector<int> cell_numpts;
cell_numpts.reserve(num_cells);
{
pm["Name"] = "connectivity";
Tag t("DataArray", pm, os);
int hf = 0;
for (int c = 0; c < num_cells; ++c) {
std::set<int> cell_pts;
for (; hf < grid->cell_facepos[c+1]; ++hf) {
int f = grid->cell_faces[hf];
const int* fnbeg = grid->face_nodes + grid->face_nodepos[f];
const int* fnend = grid->face_nodes + grid->face_nodepos[f+1];
cell_pts.insert(fnbeg, fnend);
}
cell_numpts.push_back(cell_pts.size());
Tag::indent(os);
std::copy(cell_pts.begin(), cell_pts.end(),
std::ostream_iterator<int>(os, " "));
os << '\n';
}
}
{
pm["Name"] = "offsets";
Tag t("DataArray", pm, os);
int offset = 0;
const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) {
Tag::indent(os);
}
offset += cell_numpts[c];
os << offset << ' ';
if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) {
os << '\n';
}
}
}
std::vector<int> cell_foffsets;
cell_foffsets.reserve(num_cells);
{
pm["Name"] = "faces";
Tag t("DataArray", pm, os);
const int* fp = grid->cell_facepos;
int offset = 0;
for (int c = 0; c < num_cells; ++c) {
Tag::indent(os);
os << fp[c+1] - fp[c] << '\n';
++offset;
for (int hf = fp[c]; hf < fp[c+1]; ++hf) {
int f = grid->cell_faces[hf];
const int* np = grid->face_nodepos;
int f_num_pts = np[f+1] - np[f];
Tag::indent(os);
os << f_num_pts << ' ';
++offset;
std::copy(grid->face_nodes + np[f],
grid->face_nodes + np[f+1],
std::ostream_iterator<int>(os, " "));
os << '\n';
offset += f_num_pts;
}
cell_foffsets.push_back(offset);
}
}
{
pm["Name"] = "faceoffsets";
Tag t("DataArray", pm, os);
const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) {
Tag::indent(os);
}
os << cell_foffsets[c] << ' ';
if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) {
os << '\n';
}
}
}
{
pm["type"] = "UInt8";
pm["Name"] = "types";
Tag t("DataArray", pm, os);
const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) {
Tag::indent(os);
}
os << "42 ";
if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) {
os << '\n';
}
}
}
}
{
pm.clear();
pm["Scalars"] = "saturation";
Tag celldatatag("CellData", pm, os);
pm.clear();
pm["type"] = "Int32";
pm["NumberOfComponents"] = "1";
pm["format"] = "ascii";
pm["type"] = "Float64";
{
pm["Name"] = "pressure";
Tag ptag("DataArray", pm, os);
const int num_per_line = 5;
for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) {
Tag::indent(os);
}
os << state.pressure()[c] << ' ';
if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) {
os << '\n';
}
}
}
{
pm["Name"] = "saturation";
Tag ptag("DataArray", pm, os);
const int num_per_line = 5;
for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) {
Tag::indent(os);
}
os << state.saturation()[2*c] << ' ';
if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) {
os << '\n';
}
}
}
}
}
static void toWaterSat(const std::vector<double>& sboth, std::vector<double>& sw)
{
int num = sboth.size()/2;
sw.resize(num);
for (int i = 0; i < num; ++i) {
sw[i] = sboth[2*i];
}
}
static void toBothSat(const std::vector<double>& sw, std::vector<double>& sboth)
{
int num = sw.size();
sboth.resize(2*num);
for (int i = 0; i < num; ++i) {
sboth[2*i] = sw[i];
sboth[2*i + 1] = 1.0 - sw[i];
}
}
// --------------- Types needed to define transport solver ---------------
class SimpleFluid2pWrappingProps
{
public:
SimpleFluid2pWrappingProps(const Opm::IncompPropertiesInterface& props)
: props_(props)
{
if (props.numPhases() != 2) {
THROW("SimpleFluid2pWrapper requires 2 phases.");
}
}
double density(int phase) const
{
return props_.density()[phase];
}
template <class Sat,
class Mob,
class DMob>
void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const
{
props_.relperm(1, &s[0], &c, &mob[0], &dmob[0]);
const double* mu = props_.viscosity();
mob[0] /= mu[0];
mob[1] /= mu[1];
// Recall that we use Fortran ordering for kr derivatives,
// therefore dmob[i*2 + j] is row j and column i of the
// matrix.
// Each row corresponds to a kr function, so which mu to
// divide by also depends on the row, j.
dmob[0*2 + 0] /= mu[0];
dmob[0*2 + 1] /= mu[1];
dmob[1*2 + 0] /= mu[0];
dmob[1*2 + 1] /= mu[1];
}
template <class Sat,
class Pcap,
class DPcap>
void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const
{
double pc[2];
double dpc[4];
props_.capPress(1, &s[0], &c, pc, dpc);
pcap = pc[0];
ASSERT(pc[1] == 0.0);
dpcap = dpc[0];
ASSERT(dpc[1] == 0.0);
ASSERT(dpc[2] == 0.0);
ASSERT(dpc[3] == 0.0);
}
/// \todo Properly implement s_min() and s_max().
/// We must think about how to do this in
/// the *Properties* classes.
double s_min(int c) const { (void) c; return 0.0; }
double s_max(int c) const { (void) c; return 1.0; }
private:
const Opm::IncompPropertiesInterface& props_;
};
typedef SimpleFluid2pWrappingProps TwophaseFluid;
typedef Opm::SinglePointUpwindTwoPhase<TwophaseFluid> TransportModel;
using namespace Opm::ImplicitTransportDefault;
typedef NewtonVectorCollection< ::std::vector<double> > NVecColl;
typedef JacobianSystem < struct CSRMatrix, NVecColl > JacSys;
template <class Vector>
class MaxNorm {
public:
static double
norm(const Vector& v) {
return AccumulationNorm <Vector, MaxAbs>::norm(v);
}
};
typedef Opm::ImplicitTransport<TransportModel,
JacSys ,
MaxNorm ,
VectorNegater ,
VectorZero ,
MatrixZero ,
VectorAssign > TransportSolver;
// ----------------- Main program -----------------
@ -709,8 +129,6 @@ main(int argc, char** argv)
const int num_psteps = param.getDefault("num_psteps", 1);
const double stepsize_days = param.getDefault("stepsize_days", 1.0);
const double stepsize = Opm::unit::convert::from(stepsize_days, Opm::unit::day);
const bool guess_old_solution = param.getDefault("guess_old_solution", false);
const bool use_reorder = param.getDefault("use_reorder", true);
const bool output = param.getDefault("output", true);
std::string output_dir;
if (output) {
@ -748,13 +166,8 @@ main(int argc, char** argv)
compute_porevolume(grid->c_grid(), *props, porevol);
double tot_porevol = std::accumulate(porevol.begin(), porevol.end(), 0.0);
// Extra fluid init for transport solver.
TwophaseFluid fluid(*props);
// Solvers init.
PressureSolver psolver(grid->c_grid(), *props);
TransportModel model (fluid, *grid->c_grid(), porevol, 0, guess_old_solution);
TransportSolver tsolver(model);
Opm::PressureSolver psolver(grid->c_grid(), *props);
// State-related and source-related variables init.
std::vector<double> totmob;
@ -766,27 +179,11 @@ main(int argc, char** argv)
std::vector<double> src (grid->c_grid()->number_of_cells, 0.0);
src[0] = flow_per_sec;
src[grid->c_grid()->number_of_cells - 1] = -flow_per_sec;
TransportSource* tsrc = create_transport_source(2, 2);
double ssrc[] = { 1.0, 0.0 };
double ssink[] = { 0.0, 1.0 };
double zdummy[] = { 0.0, 0.0 };
append_transport_source(0, 2, 0, src[0], ssrc, zdummy, tsrc);
append_transport_source(grid->c_grid()->number_of_cells - 1, 2, 0,
src.back(), ssink, zdummy, tsrc);
std::vector<double> reorder_src = src;
// Control init.
Opm::ImplicitTransportDetails::NRReport rpt;
Opm::ImplicitTransportDetails::NRControl ctrl;
double current_time = 0.0;
double total_time = stepsize*num_psteps;
ctrl.max_it = param.getDefault("max_it", 20);
ctrl.verbosity = param.getDefault("verbosity", 0);
ctrl.max_it_ls = param.getDefault("max_it_ls", 5);
// Linear solver init.
using Opm::ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver;
CSRMatrixUmfpackSolver linsolve;
// Warn if any parameters are unused.
if (param.anyUnused()) {
@ -817,28 +214,23 @@ main(int argc, char** argv)
compute_totmob(*props, state.saturation(), totmob);
psolver.solve(grid->c_grid(), totmob, src, state);
if (use_reorder) {
toWaterSat(state.saturation(), reorder_sat);
// We must treat reorder_src here,
// if we are to handle anything but simple water
// injection, since it is expected to be
// equal to total outflow (if negative)
// and water inflow (if positive).
// Also, for anything but noflow boundaries,
// boundary flows must be accumulated into
// source term following the same convention.
twophasetransport(&porevol[0],
&reorder_src[0],
stepsize,
const_cast<UnstructuredGrid*>(grid->c_grid()),
props.get(),
&state.faceflux()[0],
&reorder_sat[0]);
toBothSat(reorder_sat, state.saturation());
} else {
tsolver.solve(*grid->c_grid(), tsrc, stepsize, ctrl, state, linsolve, rpt);
std::cout << rpt;
}
Opm::toWaterSat(state.saturation(), reorder_sat);
// We must treat reorder_src here,
// if we are to handle anything but simple water
// injection, since it is expected to be
// equal to total outflow (if negative)
// and water inflow (if positive).
// Also, for anything but noflow boundaries,
// boundary flows must be accumulated into
// source term following the same convention.
twophasetransport(&porevol[0],
&reorder_src[0],
stepsize,
const_cast<UnstructuredGrid*>(grid->c_grid()),
props.get(),
&state.faceflux()[0],
&reorder_sat[0]);
Opm::toBothSat(reorder_sat, state.saturation());
current_time += stepsize;
}
@ -846,6 +238,4 @@ main(int argc, char** argv)
if (output) {
outputState(grid->c_grid(), state, num_psteps, output_dir);
}
destroy_transport_source(tsrc);
}