diff --git a/examples/Makefile.am b/examples/Makefile.am index a6c9a677a..f17a3ce17 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -9,4 +9,6 @@ hello_SOURCES = \ SayHello.cpp polymer_reorder_SOURCES = \ -polymer_reorder.cpp +polymer_reorder.cpp \ +Utilities.hpp \ +Utilities.cpp diff --git a/examples/Utilities.cpp b/examples/Utilities.cpp new file mode 100644 index 000000000..3b5366458 --- /dev/null +++ b/examples/Utilities.cpp @@ -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 . +*/ + +#include "config.h" + +#include "Utilities.hpp" + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + + +namespace Opm +{ + + + + + + + void + compute_porevolume(const UnstructuredGrid* g, + const Opm::IncompPropertiesInterface& props, + std::vector& 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()); + } + + + void + compute_totmob(const Opm::IncompPropertiesInterface& props, + const std::vector& s, + std::vector& 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 cells(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + cells[cell] = cell; + } + std::vector 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& dims, + const std::tr1::array& cell_size, + const std::vector& pressure, + const std::vector& 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 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_ << "\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& pressure, + const std::vector& saturation, + std::ostream& os) + { + if (grid->dimensions != 3) { + THROW("Vtk output for 3d grids only"); + } + os.precision(12); + os << "\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(num_pts); + pm["NumberOfCells"] = boost::lexical_cast(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 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 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(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 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(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& sboth, std::vector& 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& sw, std::vector& 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 + diff --git a/examples/Utilities.hpp b/examples/Utilities.hpp new file mode 100644 index 000000000..5b69443be --- /dev/null +++ b/examples/Utilities.hpp @@ -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 . +*/ + +#ifndef OPM_UTILITIES_HEADER_INCLUDED +#define OPM_UTILITIES_HEADER_INCLUDED + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + + +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& zcorn = deck.getFloatingPointValue("ZCORN"); + const std::vector& coord = deck.getFloatingPointValue("COORD"); + const std::vector& actnum = deck.getIntegerValue("ACTNUM"); + std::vector 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(g); + tpfa_htrans_compute(gg, props.permeability(), &htrans_[0]); + + h_ = ifs_tpfa_construct(gg); + } + + ~PressureSolver() + { + ifs_tpfa_destroy(h_); + } + + template + void + solve(const UnstructuredGrid* g , + const ::std::vector& totmob, + const ::std::vector& src , + State& state ) + { + UnstructuredGrid* gg = const_cast(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 htrans_; + ::std::vector trans_ ; + ::std::vector gpress_; + + struct ifs_tpfa_data* h_; + }; + + + + + void + compute_porevolume(const UnstructuredGrid* g, + const Opm::IncompPropertiesInterface& props, + std::vector& porevol); + + + void + compute_totmob(const Opm::IncompPropertiesInterface& props, + const std::vector& s, + std::vector& totmob); + + + + + void writeVtkDataAllCartesian(const std::tr1::array& dims, + const std::tr1::array& cell_size, + const std::vector& pressure, + const std::vector& saturation, + std::ostream& vtk_file); + + + + void writeVtkDataGeneralGrid(const UnstructuredGrid* grid, + const std::vector& pressure, + const std::vector& saturation, + std::ostream& os); + + + + + void toWaterSat(const std::vector& sboth, std::vector& sw); + + void toBothSat(const std::vector& sw, std::vector& sboth); + + + +} // namespace Opm + + + + +#endif // OPM_UTILITIES_HEADER_INCLUDED diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp index 379904dec..685de3a96 100644 --- a/examples/polymer_reorder.cpp +++ b/examples/polymer_reorder.cpp @@ -1,21 +1,6 @@ -/*=========================================================================== -// -// File: spu_2p.cpp -// -// Created: 2011-10-05 10:29:01+0200 -// -// Authors: Ingeborg S. Ligaarden -// Jostein R. Natvig -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Bård Skaflestad -// -//==========================================================================*/ - - /* - 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 . */ + #include "config.h" -#include +#include "Utilities.hpp" #include #include @@ -49,14 +35,7 @@ #include #include -#include #include -#include -#include -#include -#include -#include -#include #include @@ -75,79 +54,11 @@ #include #include #include +#include -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& zcorn = deck.getFloatingPointValue("ZCORN"); - const std::vector& coord = deck.getFloatingPointValue("COORD"); - const std::vector& actnum = deck.getIntegerValue("ACTNUM"); - std::vector 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(g); - tpfa_htrans_compute(gg, props.permeability(), &htrans_[0]); - - h_ = ifs_tpfa_construct(gg); - } - - ~PressureSolver() - { - ifs_tpfa_destroy(h_); - } - - template - void - solve(const UnstructuredGrid* g , - const ::std::vector& totmob, - const ::std::vector& src , - State& state ) - { - UnstructuredGrid* gg = const_cast(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 htrans_; - ::std::vector trans_ ; - ::std::vector gpress_; - - struct ifs_tpfa_data* h_; -}; - - - - -static void -compute_porevolume(const UnstructuredGrid* g, - const Opm::IncompPropertiesInterface& props, - std::vector& 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()); -} - - -static void -compute_totmob(const Opm::IncompPropertiesInterface& props, - const std::vector& s, - std::vector& 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 cells(num_cells); - for (int cell = 0; cell < num_cells; ++cell) { - cells[cell] = cell; - } - std::vector 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 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 -void writeVtkDataAllCartesian(const std::tr1::array& dims, - const std::tr1::array& 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 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_ << "\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 -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 << "\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(num_pts); - pm["NumberOfCells"] = boost::lexical_cast(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 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 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(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 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(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& sboth, std::vector& 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& sw, std::vector& 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 - 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 - 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 TransportModel; - -using namespace Opm::ImplicitTransportDefault; - -typedef NewtonVectorCollection< ::std::vector > NVecColl; -typedef JacobianSystem < struct CSRMatrix, NVecColl > JacSys; - -template -class MaxNorm { -public: - static double - norm(const Vector& v) { - return AccumulationNorm ::norm(v); - } -}; - -typedef Opm::ImplicitTransport 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 totmob; @@ -766,27 +179,11 @@ main(int argc, char** argv) std::vector 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 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(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(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); }