commit dc1b9e4b50b66b2e7ee19c3a8dca0735ded1d575 Author: Halvor M. Nilsen Date: Fri Oct 7 10:54:25 2011 +0200 Move 'common' directory into 'dune'. diff --git a/dune/common/Average.hpp b/dune/common/Average.hpp new file mode 100644 index 00000000..be2321d5 --- /dev/null +++ b/dune/common/Average.hpp @@ -0,0 +1,98 @@ +//=========================================================================== +// +// File: Average.hpp<2> +// +// Created: Wed Jul 1 12:25:57 2009 +// +// Author(s): Atgeirr F Rasmussen +// Bård Skaflestad +// +// $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 . +*/ + +#ifndef OPENRS_AVERAGE_HEADER +#define OPENRS_AVERAGE_HEADER + + +#include +#include +#include + +namespace Dune { + namespace utils { + + + + /// Computes the arithmetic average: + /// a_A(t_1, t_2) = (t_1 + t_2)/2. + template + Tresult arithmeticAverage(const T& t1, const T& t2) + { + // To avoid some user errors, we disallow taking averages of + // integral values. + BOOST_STATIC_ASSERT(std::tr1::is_integral::value == false); + Tresult retval(t1); + retval += t2; + retval *= 0.5; + return retval; + } + + + + /// Computes the geometric average: + /// a_G(t_1, t_2) = \sqrt{t_1 t_2}. + /// Since we use std::sqrt(), this function will only + /// compile when T is convertible to and from double. + template + T geometricAverage(const T& t1, const T& t2) + { + // To avoid some user errors, we disallow taking averages of + // integral values. + BOOST_STATIC_ASSERT(std::tr1::is_integral::value == false); + return std::sqrt(t1*t2); + } + + + + /// Computes the harmonic average: + /// a_H(t_1, t_2) = \frac{2}{1/t_1 + 1/t_2} = \frac{2 t_1 t_2}{t_1 + t_2}. + template + T harmonicAverage(const T& t1, const T& t2) + { + // To avoid some user errors, we disallow taking averages of + // integral values. + BOOST_STATIC_ASSERT(std::tr1::is_integral::value == false); + return (2*t1*t2)/(t1 + t2); + } + + + + } // namespace utils +} // namespace Dune + + + +#endif // OPENRS_AVERAGE_HEADER diff --git a/dune/common/CornerpointChopper.hpp b/dune/common/CornerpointChopper.hpp new file mode 100644 index 00000000..cdab17d1 --- /dev/null +++ b/dune/common/CornerpointChopper.hpp @@ -0,0 +1,399 @@ +/* + Copyright 2010 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_CORNERPOINTCHOPPER_HEADER_INCLUDED +#define OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include +#include + +namespace Dune +{ + + class CornerPointChopper + { + public: + CornerPointChopper(const std::string& file) + : parser_(file, false) + { + for (int dd = 0; dd < 3; ++dd) { + dims_[dd] = parser_.getSPECGRID().dimensions[dd]; + } + + int layersz = 8*dims_[0]*dims_[1]; + const std::vector& ZCORN = parser_.getFloatingPointValue("ZCORN"); + botmax_ = *std::max_element(ZCORN.begin(), ZCORN.begin() + layersz/2); + topmin_ = *std::min_element(ZCORN.begin() + dims_[2]*layersz - layersz/2, + ZCORN.begin() + dims_[2]*layersz); + + abszmax_ = *std::max_element(ZCORN.begin(), ZCORN.end()); + abszmin_ = *std::min_element(ZCORN.begin(), ZCORN.end()); + + std::cout << "Parsed grdecl file with dimensions (" + << dims_[0] << ", " << dims_[1] << ", " << dims_[2] << ")" << std::endl; + } + + + + + const int* dimensions() const + { + return dims_; + } + + + + + const int* newDimensions() const + { + return new_dims_; + } + + + + + const std::pair zLimits() const + { + return std::make_pair(botmax_, topmin_); + } + + const std::pair abszLimits() const + { + return std::make_pair(abszmin_, abszmax_); + } + + + void verifyInscribedShoebox(int imin, int ilen, int imax, + int jmin, int jlen, int jmax, + double zmin, double zlen, double zmax) + { + if (imin < 0) { + std::cerr << "Error! imin < 0 (imin = " << imin << ")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (ilen > dims_[0]) { + std::cerr << "Error! ilen larger than grid (ilen = " << ilen <<")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (imax > dims_[0]) { + std::cerr << "Error! imax larger than input grid (imax = " << imax << ")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (jmin < 0) { + std::cerr << "Error! jmin < 0 (jmin = " << jmin << ")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (jlen > dims_[1]) { + std::cerr << "Error! jlen larger than grid (jlen = " << jlen <<")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (jmax > dims_[1]) { + std::cerr << "Error! jmax larger than input grid (jmax = " << jmax << ")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (zmin < abszmin_) { + std::cerr << "Error! zmin ("<< zmin << ") less than minimum ZCORN value ("<< abszmin_ << ")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (zmax > abszmax_) { + std::cerr << "Error! zmax ("<< zmax << ") larger than maximal ZCORN value ("<< abszmax_ << ")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (zlen > (abszmax_ - abszmin_)) { + std::cerr << "Error! zlen ("<< zlen <<") larger than maximal ZCORN (" << abszmax_ << ") minus minimal ZCORN ("<< abszmin_ <<")\n"; + throw std::runtime_error("Inconsistent user input."); + } + } + + void chop(int imin, int imax, int jmin, int jmax, double zmin, double zmax, bool resettoorigin=true) + { + new_dims_[0] = imax - imin; + new_dims_[1] = jmax - jmin; + + // Filter the coord field + const std::vector& COORD = parser_.getFloatingPointValue("COORD"); + int num_coord = COORD.size(); + if (num_coord != 6*(dims_[0] + 1)*(dims_[1] + 1)) { + std::cerr << "Error! COORD size (" << COORD.size() << ") not consistent with SPECGRID\n"; + throw std::runtime_error("Inconsistent COORD and SPECGRID."); + } + int num_new_coord = 6*(new_dims_[0] + 1)*(new_dims_[1] + 1); + double x_correction = COORD[6*((dims_[0] + 1)*jmin + imin)]; + double y_correction = COORD[6*((dims_[0] + 1)*jmin + imin) + 1]; + new_COORD_.resize(num_new_coord, 1e100); + for (int j = jmin; j < jmax + 1; ++j) { + for (int i = imin; i < imax + 1; ++i) { + int pos = (dims_[0] + 1)*j + i; + int new_pos = (new_dims_[0] + 1)*(j-jmin) + (i-imin); + // Copy all 6 coordinates for a pillar. + std::copy(COORD.begin() + 6*pos, COORD.begin() + 6*(pos + 1), new_COORD_.begin() + 6*new_pos); + if (resettoorigin) { + // Substract lowest x value from all X-coords, similarly for y, and truncate in z-direction + new_COORD_[6*new_pos] -= x_correction; + new_COORD_[6*new_pos + 1] -= y_correction; + new_COORD_[6*new_pos + 2] = 0; + new_COORD_[6*new_pos + 3] -= x_correction; + new_COORD_[6*new_pos + 4] -= y_correction; + new_COORD_[6*new_pos + 5] = zmax-zmin; + } + } + } + + // Get the z limits, check if they must be changed to make a shoe-box. + // This means that zmin must be greater than or equal to the highest + // coordinate of the bottom surface, while zmax must be less than or + // equal to the lowest coordinate of the top surface. + int layersz = 8*dims_[0]*dims_[1]; + const std::vector& ZCORN = parser_.getFloatingPointValue("ZCORN"); + int num_zcorn = ZCORN.size(); + if (num_zcorn != layersz*dims_[2]) { + std::cerr << "Error! ZCORN size (" << ZCORN.size() << ") not consistent with SPECGRID\n"; + throw std::runtime_error("Inconsistent ZCORN and SPECGRID."); + } + + zmin = std::max(zmin, botmax_); + zmax = std::min(zmax, topmin_); + if (zmin >= zmax) { + std::cerr << "Error: zmin >= zmax (zmin = " << zmin << ", zmax = " << zmax << ")\n"; + throw std::runtime_error("zmin >= zmax"); + } + std::cout << "Chopping subsample, i: (" << imin << "--" << imax << ") j: (" << jmin << "--" << jmax << ") z: (" << zmin << "--" << zmax << ")" << std::endl; + + // We must find the maximum and minimum k value for the given z limits. + // First, find the first layer with a z-coordinate strictly above zmin. + int kmin = -1; + for (int k = 0; k < dims_[2]; ++k) { + double layer_max = *std::max_element(ZCORN.begin() + k*layersz, ZCORN.begin() + (k + 1)*layersz); + if (layer_max > zmin) { + kmin = k; + break; + } + } + // Then, find the last layer with a z-coordinate strictly below zmax. + int kmax = -1; + for (int k = dims_[2]; k > 0; --k) { + double layer_min = *std::min_element(ZCORN.begin() + (k - 1)*layersz, ZCORN.begin() + k*layersz); + if (layer_min < zmax) { + kmax = k; + break; + } + } + new_dims_[2] = kmax - kmin; + + // Filter the ZCORN field, build mapping from new to old cells. + double z_origin_correction = 0.0; + if (resettoorigin) { + z_origin_correction = zmin; + } + new_ZCORN_.resize(8*new_dims_[0]*new_dims_[1]*new_dims_[2], 1e100); + new_to_old_cell_.resize(new_dims_[0]*new_dims_[1]*new_dims_[2], -1); + int cellcount = 0; + int delta[3] = { 1, 2*dims_[0], 4*dims_[0]*dims_[1] }; + int new_delta[3] = { 1, 2*new_dims_[0], 4*new_dims_[0]*new_dims_[1] }; + for (int k = kmin; k < kmax; ++k) { + for (int j = jmin; j < jmax; ++j) { + for (int i = imin; i < imax; ++i) { + new_to_old_cell_[cellcount++] = dims_[0]*dims_[1]*k + dims_[0]*j + i; + int old_ix = 2*(i*delta[0] + j*delta[1] + k*delta[2]); + int new_ix = 2*((i-imin)*new_delta[0] + (j-jmin)*new_delta[1] + (k-kmin)*new_delta[2]); + int old_indices[8] = { old_ix, old_ix + delta[0], + old_ix + delta[1], old_ix + delta[1] + delta[0], + old_ix + delta[2], old_ix + delta[2] + delta[0], + old_ix + delta[2] + delta[1], old_ix + delta[2] + delta[1] + delta[0] }; + int new_indices[8] = { new_ix, new_ix + new_delta[0], + new_ix + new_delta[1], new_ix + new_delta[1] + new_delta[0], + new_ix + new_delta[2], new_ix + new_delta[2] + new_delta[0], + new_ix + new_delta[2] + new_delta[1], new_ix + new_delta[2] + new_delta[1] + new_delta[0] }; + for (int cc = 0; cc < 8; ++cc) { + new_ZCORN_[new_indices[cc]] = std::min(zmax, std::max(zmin, ZCORN[old_indices[cc]])) - z_origin_correction; + } + } + } + } + + filterIntegerField("ACTNUM", new_ACTNUM_); + filterDoubleField("PORO", new_PORO_); + filterDoubleField("PERMX", new_PERMX_); + filterDoubleField("PERMY", new_PERMY_); + filterDoubleField("PERMZ", new_PERMZ_); + filterIntegerField("SATNUM", new_SATNUM_); + } + + + + /// Return a subparser with fields corresponding to the selected subset. + /// Note that the returned parser is NOT converted to SI, that must be done + /// by the user afterwards with the parser's convertToSI() method. + EclipseGridParser subparser() + { + if (parser_.hasField("FIELD") || parser_.hasField("LAB") || parser_.hasField("PVT-M")) { + THROW("CornerPointChopper::subparser() cannot handle any eclipse unit system other than METRIC."); + } + + EclipseGridParser sp; + boost::shared_ptr sg(new SPECGRID); + for (int dd = 0; dd < 3; ++dd) { + sg->dimensions[dd] = new_dims_[dd]; + } + sp.setSpecialField("SPECGRID", sg); + sp.setFloatingPointField("COORD", new_COORD_); + sp.setFloatingPointField("ZCORN", new_ZCORN_); + if (!new_ACTNUM_.empty()) sp.setIntegerField("ACTNUM", new_ACTNUM_); + if (!new_PORO_.empty()) sp.setFloatingPointField("PORO", new_PORO_); + if (!new_PERMX_.empty()) sp.setFloatingPointField("PERMX", new_PERMX_); + if (!new_PERMY_.empty()) sp.setFloatingPointField("PERMY", new_PERMY_); + if (!new_PERMZ_.empty()) sp.setFloatingPointField("PERMZ", new_PERMZ_); + if (!new_SATNUM_.empty()) sp.setIntegerField("SATNUM", new_SATNUM_); + sp.computeUnits(); // Always METRIC, since that is default. + return sp; + } + + + + + void writeGrdecl(const std::string& filename) + { + // Output new versions of SPECGRID, COORD, ZCORN, ACTNUM, PERMX, PORO, SATNUM. + std::ofstream out(filename.c_str()); + if (!out) { + std::cerr << "Could not open file " << filename << "\n"; + throw std::runtime_error("Could not open output file."); + } + out << "SPECGRID\n" << new_dims_[0] << ' ' << new_dims_[1] << ' ' << new_dims_[2] + << " 1 F\n/\n\n"; + out << "COORD\n"; + int num_new_coord = new_COORD_.size(); + for (int i = 0; i < num_new_coord/6; ++i) { + for (int j = 0; j < 6; ++j) { + out << " " << new_COORD_[6*i + j]; + } + out << '\n'; + } + out << "/\n\n"; + out << "ZCORN\n"; + int num_new_zcorn = new_ZCORN_.size(); + assert(num_new_zcorn%8 == 0); + for (int i = 0; i < num_new_zcorn/8; ++i) { + for (int j = 0; j < 8; ++j) { + out << " " << new_ZCORN_[8*i + j]; + } + out << '\n'; + } + out << "/\n\n"; + + outputField(out, new_ACTNUM_, "ACTNUM"); + outputField(out, new_PORO_, "PORO"); + outputField(out, new_PERMX_, "PERMX"); + outputField(out, new_PERMY_, "PERMY"); + outputField(out, new_PERMZ_, "PERMZ"); + outputField(out, new_SATNUM_, "SATNUM"); + } + + private: + EclipseGridParser parser_; + double botmax_; + double topmin_; + double abszmin_; + double abszmax_; + std::vector new_COORD_; + std::vector new_ZCORN_; + std::vector new_ACTNUM_; + std::vector new_PORO_; + std::vector new_PERMX_; + std::vector new_PERMY_; + std::vector new_PERMZ_; + std::vector new_SATNUM_; + int dims_[3]; + int new_dims_[3]; + std::vector new_to_old_cell_; + + + template + void outputField(std::ostream& os, + const std::vector& field, + const std::string& keyword) + { + if (field.empty()) return; + os << keyword << '\n'; + int sz = field.size(); + T last = std::numeric_limits::max(); + int repeats = 0; + for (int i = 0; i < sz; ++i) { + T val = field[i]; + if (val == last) { + ++repeats; + } else { + if (repeats == 1) { + os << last << '\n'; + } else if (repeats > 1) { + os << repeats << '*' << last << '\n'; + } + last = val; + repeats = 1; + } + } + if (repeats == 1) { + os << last << '\n'; + } else if (repeats > 1) { + os << repeats << '*' << last << '\n'; + } + os << "/\n\n"; + } + + + template + void filterField(const std::vector& field, + std::vector& output_field) + { + int sz = new_to_old_cell_.size(); + output_field.resize(sz); + for (int i = 0; i < sz; ++i) { + output_field[i] = field[new_to_old_cell_[i]]; + } + } + + void filterDoubleField(const std::string& keyword, std::vector& output_field) + { + if (parser_.hasField(keyword)) { + const std::vector& field = parser_.getFloatingPointValue(keyword); + filterField(field, output_field); + } + } + + void filterIntegerField(const std::string& keyword, std::vector& output_field) + { + if (parser_.hasField(keyword)) { + const std::vector& field = parser_.getIntegerValue(keyword); + filterField(field, output_field); + } + } + + }; + +} + + + + +#endif // OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED diff --git a/dune/common/EclipseGridInspector.cpp b/dune/common/EclipseGridInspector.cpp new file mode 100644 index 00000000..ea80bb85 --- /dev/null +++ b/dune/common/EclipseGridInspector.cpp @@ -0,0 +1,335 @@ +//=========================================================================== +// +// File: EclipseGridInspector.C +// +// Created: Mon Jun 2 12:17:51 2008 +// +// Author: Atgeirr F Rasmussen +// +// $Date$ +// +// $Revision$ +// +// Revision: $Id: EclipseGridInspector.C,v 1.2 2008/08/18 14:16:13 atgeirr Exp $ +// +//=========================================================================== + +/* + 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 . +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif +#include +#include +#include +#include +#include +#include "EclipseGridInspector.hpp" +#include "EclipseGridParser.hpp" +#include "SpecialEclipseFields.hpp" + +namespace Dune +{ + +EclipseGridInspector::EclipseGridInspector(const EclipseGridParser& parser) + : parser_(parser) +{ + std::vector keywords; + keywords.push_back("COORD"); + keywords.push_back("ZCORN"); + + if (!parser_.hasFields(keywords)) { + THROW("Needed field is missing in file"); + } + + if (parser_.hasField("SPECGRID")) { + const SPECGRID& sgr = parser.getSPECGRID(); + logical_gridsize_[0] = sgr.dimensions[0]; + logical_gridsize_[1] = sgr.dimensions[1]; + logical_gridsize_[2] = sgr.dimensions[2]; + } else if (parser_.hasField("DIMENS")) { + const std::vector& dim = parser.getIntegerValue("DIMENS"); + logical_gridsize_[0] = dim[0]; + logical_gridsize_[1] = dim[1]; + logical_gridsize_[2] = dim[2]; + } else { + THROW("Found neither SPECGRID nor DIMENS in file. At least one is needed."); + } + +} + +/** + Return the dip slopes for the cell relative to xy-plane in x- and y- direction. + Dip slope is average rise in positive x-direction over cell length in x-direction. + Similarly for y. + + Current implementation is for vertical pillars, but is not difficult to fix. + + @returns a std::pair with x-dip in first component and y-dip in second. +*/ +std::pair EclipseGridInspector::cellDips(int i, int j, int k) const +{ + checkLogicalCoords(i, j, k); + const std::vector& pillc = parser_.getFloatingPointValue("COORD"); + int num_pillars = (logical_gridsize_[0] + 1)*(logical_gridsize_[1] + 1); + if (6*num_pillars != int(pillc.size())) { + throw std::runtime_error("Wrong size of COORD field."); + } + const std::vector& z = parser_.getFloatingPointValue("ZCORN"); + int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2]; + if (8*num_cells != int(z.size())) { + throw std::runtime_error("Wrong size of ZCORN field"); + } + + // Pick ZCORN-value for all 8 corners of the given cell + boost::array cellz = cellZvals(i, j, k); + + // Compute rise in positive x-direction for all four edges (and then find mean) + // Current implementation is for regularly placed and vertical pillars! + int numxpill = logical_gridsize_[0] + 1; + int pix = i + j*numxpill; + double cell_xlength = pillc[6*(pix + 1)] - pillc[6*pix]; + flush(std::cout); + double xrise[4] = { (cellz[1] - cellz[0])/cell_xlength, // LLL -> HLL + (cellz[3] - cellz[2])/cell_xlength, // LHL -> HHL + (cellz[5] - cellz[4])/cell_xlength, // LLH -> HLH + (cellz[7] - cellz[6])/cell_xlength}; // LHH -> HHH + + double cell_ylength = pillc[6*(pix + numxpill) + 1] - pillc[6*pix + 1]; + double yrise[4] = { (cellz[2] - cellz[0])/cell_ylength, // LLL -> LHL + (cellz[3] - cellz[1])/cell_ylength, // HLL -> HHL + (cellz[6] - cellz[4])/cell_ylength, // LLH -> LHH + (cellz[7] - cellz[5])/cell_ylength}; // HLH -> HHH + + + // Now ignore those edges that touch the global top or bottom surface + // of the entire grdecl model. This is to avoid bias, as these edges probably + // don't follow an overall dip for the model if it exists. + int x_edges = 4; + int y_edges = 4; + boost::array gridlimits = getGridLimits(); + double zmin = gridlimits[4]; + double zmax = gridlimits[5]; + // LLL -> HLL + if ((cellz[1] == zmin) || (cellz[0] == zmin)) { + xrise[0] = 0; x_edges--; + } + // LHL -> HHL + if ((cellz[2] == zmin) || (cellz[3] == zmin)) { + xrise[1] = 0; x_edges--; + } + // LLH -> HLH + if ((cellz[4] == zmax) || (cellz[5] == zmax)) { + xrise[2] = 0; x_edges--; + } + // LHH -> HHH + if ((cellz[6] == zmax) || (cellz[7] == zmax)) { + xrise[3] = 0; x_edges--; + } + // LLL -> LHL + if ((cellz[0] == zmin) || (cellz[2] == zmin)) { + yrise[0] = 0; y_edges--; + } + // HLL -> HHL + if ((cellz[1] == zmin) || (cellz[3] == zmin)) { + yrise[1] = 0; y_edges--; + } + // LLH -> LHH + if ((cellz[6] == zmax) || (cellz[4] == zmax)) { + yrise[2] = 0; y_edges--; + } + // HLH -> HHH + if ((cellz[7] == zmax) || (cellz[5] == zmax)) { + yrise[3] = 0; y_edges--; + } + + return std::make_pair( (xrise[0] + xrise[1] + xrise[2] + xrise[3])/x_edges, + (yrise[0] + yrise[1] + yrise[2] + yrise[3])/y_edges); +} +/** + Wrapper for cellDips(i, j, k). +*/ +std::pair EclipseGridInspector::cellDips(int cell_idx) const +{ + boost::array idxs = cellIdxToLogicalCoords(cell_idx); + return cellDips(idxs[0], idxs[1], idxs[2]); +} + +boost::array EclipseGridInspector::cellIdxToLogicalCoords(int cell_idx) const +{ + + int i,j,k; // Position of cell in cell hierarchy + int horIdx = (cell_idx+1) - int(std::floor(((double)(cell_idx+1))/((double)(logical_gridsize_[0]*logical_gridsize_[1]))))*logical_gridsize_[0]*logical_gridsize_[1]; // index in the corresponding horizon + if (horIdx == 0) { + horIdx = logical_gridsize_[0]*logical_gridsize_[1]; + } + i = horIdx - int(std::floor(((double)horIdx)/((double)logical_gridsize_[0])))*logical_gridsize_[0]; + if (i == 0) { + i = logical_gridsize_[0]; + } + j = (horIdx-i)/logical_gridsize_[0]+1; + k = ((cell_idx+1)-logical_gridsize_[0]*(j-1)-1)/(logical_gridsize_[0]*logical_gridsize_[1])+1; + + boost::array a = {{i-1, j-1, k-1}}; + return a; //boost::array {{i-1, j-1, k-1}}; +} + +double EclipseGridInspector::cellVolumeVerticalPillars(int i, int j, int k) const +{ + // Checking parameters and obtaining values from parser. + checkLogicalCoords(i, j, k); + const std::vector& pillc = parser_.getFloatingPointValue("COORD"); + int num_pillars = (logical_gridsize_[0] + 1)*(logical_gridsize_[1] + 1); + if (6*num_pillars != int(pillc.size())) { + throw std::runtime_error("Wrong size of COORD field."); + } + const std::vector& z = parser_.getFloatingPointValue("ZCORN"); + int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2]; + if (8*num_cells != int(z.size())) { + throw std::runtime_error("Wrong size of ZCORN field"); + } + + // Computing the base area as half the 2d cross product of the diagonals. + int numxpill = logical_gridsize_[0] + 1; + int pix = i + j*numxpill; + double px[4] = { pillc[6*pix], + pillc[6*(pix + 1)], + pillc[6*(pix + numxpill)], + pillc[6*(pix + numxpill + 1)] }; + double py[4] = { pillc[6*pix + 1], + pillc[6*(pix + 1) + 1], + pillc[6*(pix + numxpill) + 1], + pillc[6*(pix + numxpill + 1) + 1] }; + double diag1[2] = { px[3] - px[0], py[3] - py[0] }; + double diag2[2] = { px[2] - px[1], py[2] - py[1] }; + double area = 0.5*(diag1[0]*diag2[1] - diag1[1]*diag2[0]); + + // Computing the average of the z-differences along each pillar. + int delta[3] = { 1, + 2*logical_gridsize_[0], + 4*logical_gridsize_[0]*logical_gridsize_[1] }; + int ix = 2*(i*delta[0] + j*delta[1] + k*delta[2]); + double cellz[8] = { z[ix], z[ix + delta[0]], + z[ix + delta[1]], z[ix + delta[1] + delta[0]], + z[ix + delta[2]], z[ix + delta[2] + delta[0]], + z[ix + delta[2] + delta[1]], z[ix + delta[2] + delta[1] + delta[0]] }; + double diffz[4] = { cellz[4] - cellz[0], + cellz[5] - cellz[1], + cellz[6] - cellz[2], + cellz[7] - cellz[3] }; + double averzdiff = 0.25*std::accumulate(diffz, diffz + 4, 0.0); + return averzdiff*area; +} + + +double EclipseGridInspector::cellVolumeVerticalPillars(int cell_idx) const +{ + boost::array idxs = cellIdxToLogicalCoords(cell_idx); + return cellVolumeVerticalPillars(idxs[0], idxs[1], idxs[2]); +} + +void EclipseGridInspector::checkLogicalCoords(int i, int j, int k) const +{ + if (i < 0 || i >= logical_gridsize_[0]) + throw std::runtime_error("First coordinate out of bounds"); + if (j < 0 || j >= logical_gridsize_[1]) + throw std::runtime_error("Second coordinate out of bounds"); + if (k < 0 || k >= logical_gridsize_[2]) + throw std::runtime_error("Third coordinate out of bounds"); +} + + +boost::array EclipseGridInspector::getGridLimits() const +{ + if (! (parser_.hasField("COORD") && parser_.hasField("ZCORN") && parser_.hasField("SPECGRID")) ) { + throw std::runtime_error("EclipseGridInspector: Grid does not have SPECGRID, COORD, and ZCORN, can't find dimensions."); + } + + std::vector coord = parser_.getFloatingPointValue("COORD"); + std::vector zcorn = parser_.getFloatingPointValue("ZCORN"); + + double xmin = +DBL_MAX; + double xmax = -DBL_MAX; + double ymin = +DBL_MAX; + double ymax = -DBL_MAX; + + + int pillars = (logical_gridsize_[0]+1) * (logical_gridsize_[1]+1); + + for (int pillarindex = 0; pillarindex < pillars; ++pillarindex) { + if (coord[pillarindex * 6 + 0] > xmax) + xmax = coord[pillarindex * 6 + 0]; + if (coord[pillarindex * 6 + 0] < xmin) + xmin = coord[pillarindex * 6 + 0]; + if (coord[pillarindex * 6 + 1] > ymax) + ymax = coord[pillarindex * 6 + 1]; + if (coord[pillarindex * 6 + 1] < ymin) + ymin = coord[pillarindex * 6 + 1]; + if (coord[pillarindex * 6 + 3] > xmax) + xmax = coord[pillarindex * 6 + 3]; + if (coord[pillarindex * 6 + 3] < xmin) + xmin = coord[pillarindex * 6 + 3]; + if (coord[pillarindex * 6 + 4] > ymax) + ymax = coord[pillarindex * 6 + 4]; + if (coord[pillarindex * 6 + 4] < ymin) + ymin = coord[pillarindex * 6 + 4]; + } + + boost::array gridlimits = {{ xmin, xmax, ymin, ymax, + *min_element(zcorn.begin(), zcorn.end()), + *max_element(zcorn.begin(), zcorn.end()) }}; + return gridlimits; +} + + + +boost::array EclipseGridInspector::gridSize() const +{ + boost::array retval = {{ logical_gridsize_[0], + logical_gridsize_[1], + logical_gridsize_[2] }}; + return retval; +} + + +boost::array EclipseGridInspector::cellZvals(int i, int j, int k) const +{ + // Get the zcorn field. + const std::vector& z = parser_.getFloatingPointValue("ZCORN"); + int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2]; + if (8*num_cells != int(z.size())) { + throw std::runtime_error("Wrong size of ZCORN field"); + } + + // Make the coordinate array. + int delta[3] = { 1, + 2*logical_gridsize_[0], + 4*logical_gridsize_[0]*logical_gridsize_[1] }; + int ix = 2*(i*delta[0] + j*delta[1] + k*delta[2]); + boost::array cellz = {{ z[ix], z[ix + delta[0]], + z[ix + delta[1]], z[ix + delta[1] + delta[0]], + z[ix + delta[2]], z[ix + delta[2] + delta[0]], + z[ix + delta[2] + delta[1]], z[ix + delta[2] + delta[1] + delta[0]] }}; + return cellz; +} + + +} // namespace Dune diff --git a/dune/common/EclipseGridInspector.hpp b/dune/common/EclipseGridInspector.hpp new file mode 100644 index 00000000..c314d28e --- /dev/null +++ b/dune/common/EclipseGridInspector.hpp @@ -0,0 +1,105 @@ +//=========================================================================== +// +// File: EclipseGridInspector.h +// +// Created: Mon Jun 2 09:46:08 2008 +// +// Author: Atgeirr F Rasmussen +// +// $Date$ +// +// Revision: $Id: EclipseGridInspector.h,v 1.2 2008/08/18 14:16:12 atgeirr Exp $ +// +//=========================================================================== + +/* +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 . +*/ + +#ifndef SINTEF_ECLIPSEGRIDINSPECTOR_HEADER +#define SINTEF_ECLIPSEGRIDINSPECTOR_HEADER + +#include +#include + + +namespace Dune +{ + +/** + @brief A class for inspecting the contents of an eclipse file. + + Given an EclipseGridParser that has successfully read en Eclipse + .grdecl-file, this class may be used to answer certain queries about + its contents. + + @author Atgeirr F. Rasmussen + @date 2008/06/02 09:46:08 +*/ + +class EclipseGridParser; + +class EclipseGridInspector +{ +public: + /// Constructor taking a parser as argument. + /// The parser must already have read an Eclipse file. + EclipseGridInspector(const EclipseGridParser& parser); + + /// Assuming that the pillars are vertical, compute the + /// volume of the cell given by logical coordinates (i, j, k). + double cellVolumeVerticalPillars(int i, int j, int k) const; + + /// Assuming that the pillars are vertical, compute the + /// volume of the cell given by the cell index + double cellVolumeVerticalPillars(int cell_idx) const; + + /// Compute the average dip in x- and y-direction of the + /// cell tops and bottoms relative to the xy-plane + std::pair cellDips(int i, int j, int k) const; + std::pair cellDips(int cell_idx) const; + + // Convert global cell index to logical ijk-coordinates + boost::array cellIdxToLogicalCoords(int cell_idx) const; + + /// Returns a vector with the outer limits of grid (in the grid's unit). + /// The vector contains [xmin, xmax, ymin, ymax, zmin, zmax], as + /// read from COORDS and ZCORN + boost::array getGridLimits() const; + + /// Returns the extent of the logical cartesian grid + /// as number of cells in the (i, j, k) directions. + boost::array gridSize() const; + + /// Returns the eight z-values associated with a given cell. + /// The ordering is such that i runs fastest. That is, with + /// L = low and H = high: + /// {LLL, HLL, LHL, HHL, LLH, HLH, LHH, HHH }. + boost::array cellZvals(int i, int j, int k) const; + +private: + const EclipseGridParser& parser_; + int logical_gridsize_[3]; + void checkLogicalCoords(int i, int j, int k) const; +}; + +} // namespace Dune + +#endif // SINTEF_ECLIPSEGRIDINSPECTOR_HEADER + diff --git a/dune/common/EclipseGridParser.cpp b/dune/common/EclipseGridParser.cpp new file mode 100644 index 00000000..cbe25aa0 --- /dev/null +++ b/dune/common/EclipseGridParser.cpp @@ -0,0 +1,571 @@ +//=========================================================================== +// +// File: EclipseGridParser.C +// +// Created: Thu Dec 6 08:46:05 2007 +// +// Author: Atgeirr F Rasmussen +// +// $Date$ +// +// $Revision$ +// +// Revision: $Id: EclipseGridParser.C,v 1.4 2008/08/18 14:16:14 atgeirr Exp $ +// +//=========================================================================== + +/* + 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 . +*/ +#if HAVE_CONFIG_H +#include "config.h" +#endif +#include +#include +#include +#include +#include +#include +#include "EclipseGridParser.hpp" +#include "EclipseGridParserHelpers.hpp" +#include "SpecialEclipseFields.hpp" +#include +#include +#include + +using namespace std; + +//#define VERBOSE + +namespace Dune +{ + +// ---------- List of supported keywords ---------- + +namespace EclipseKeywords +{ + string integer_fields[] = + { string("ACTNUM"), + string("SATNUM"), + string("EQLNUM"), + string("REGNUM"), + string("ROCKTYPE"), + string("DIMENS"), + string("REGDIMS"), + string("WELLDIMS"), + string("TABDIMS"), + string("FIPNUM"), + string("GRIDFILE") + }; + const int num_integer_fields = sizeof(integer_fields) / sizeof(integer_fields[0]); + + string floating_fields[] = + { string("COORD"), string("ZCORN"), string("PERMX"), + string("PERMY"), string("PERMZ"), string("PERMXX"), + string("PERMYY"), string("PERMZZ"), string("PERMXY"), + string("PERMYZ"), string("PERMZX"), string("PORO"), + string("BULKMOD"), string("YOUNGMOD"), string("LAMEMOD"), + string("SHEARMOD"), string("POISSONMOD"), string("PWAVEMOD"), + string("MULTPV"), string("PRESSURE"), string("SGAS"), + string("SWAT") + }; + const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]); + + string special_fields[] = + { string("SPECGRID"), string("FAULTS"), string("MULTFLT"), + string("TITLE"), string("START"), string("DATES"), + string("DENSITY"), string("PVDG"), string("PVDO"), + string("PVTG"), string("PVTO"), string("PVTW"), + string("SGOF"), string("SWOF"), string("ROCK"), + string("ROCKTAB"), string("WELSPECS"), string("COMPDAT"), + string("WCONINJE"), string("WCONPROD"), string("WELTARG"), + string("EQUIL"), string("PVCDO"), + // The following fields only have a dummy implementation + // that allows us to ignore them. + "SWFN", + "SOF2", + "TUNING" + }; + + const int num_special_fields = sizeof(special_fields) / sizeof(special_fields[0]); + + string ignore_with_data[] = + { string("MAPUNITS"), string("MAPAXES"), string("GRIDUNIT"), + string("NTG"), string("REGDIMS"), string("WELLDIMS"), + string("NSTACK"), string("SATNUM"), + string("RPTRST"), string("ROIP"), string("RWIP"), + string("RWSAT"), string("RPR"), string("WBHP"), + string("WOIR"), string("TSTEP"), string("BOX"), + string("COORDSYS"), string("PBVD") + }; + const int num_ignore_with_data = sizeof(ignore_with_data) / sizeof(ignore_with_data[0]); + + string ignore_no_data[] = + { string("RUNSPEC"), string("WATER"), string("OIL"), + string("METRIC"), string("FMTIN"), string("FMTOUT"), + string("GRID"), string("INIT"), string("NOECHO"), + string("ECHO"), string("EDIT"), string("PROPS"), + string("REGIONS"), string("SOLUTION"), string("SUMMARY"), + string("FPR"), string("FOIP"), string("FWIP"), + string("RUNSUM"), string("EXCEL"), string("SCHEDULE"), + string("END"), string("ENDBOX"), string("CONTINUE"), + string("NONNC"), string("GAS"), string("DISGAS"), + string("FIELD") + }; + const int num_ignore_no_data = sizeof(ignore_no_data) / sizeof(ignore_no_data[0]); + + string include_keywords[] = { string("INCLUDE") }; + const int num_include_keywords = sizeof(include_keywords) / sizeof(include_keywords[0]); + + +} // namespace EclipseKeywords + +namespace { + + enum FieldType { + Integer, + FloatingPoint, + SpecialField, + IgnoreWithData, + IgnoreNoData, + Include, + Unknown + }; + + inline FieldType classifyKeyword(const string& keyword) + { + using namespace EclipseKeywords; + if (count(integer_fields, integer_fields + num_integer_fields, keyword)) { + return Integer; + } else if (count(floating_fields, floating_fields + num_floating_fields, keyword)) { + return FloatingPoint; + } else if (count(special_fields, special_fields + num_special_fields, keyword)) { + return SpecialField; + } else if (count(ignore_with_data, ignore_with_data + num_ignore_with_data, keyword)) { + return IgnoreWithData; + } else if (count(ignore_no_data, ignore_no_data + num_ignore_no_data, keyword)) { + return IgnoreNoData; + } else if (count(include_keywords, include_keywords + num_include_keywords, keyword)) { + return Include; + } else { + return Unknown; + } + } + +} // anon namespace + + + +// ---------- Member functions ---------- + +/// Default constructor. +//--------------------------------------------------------------------------- +EclipseGridParser::EclipseGridParser() +//--------------------------------------------------------------------------- +{ +} + + +/// Constructor taking an eclipse filename. +//--------------------------------------------------------------------------- +EclipseGridParser::EclipseGridParser(const string& filename, bool convert_to_SI) +//--------------------------------------------------------------------------- +{ + // Store directory of filename + boost::filesystem::path p(filename); + directory_ = p.parent_path().string(); + ifstream is(filename.c_str()); + if (!is) { + cerr << "Unable to open file " << filename << endl; + throw exception(); + } + read(is, convert_to_SI); +} + + +/// Read the given stream, overwriting any previous data. +//--------------------------------------------------------------------------- +void EclipseGridParser::read(istream& is, bool convert_to_SI) +//--------------------------------------------------------------------------- +{ + integer_field_map_.clear(); + floating_field_map_.clear(); + special_field_map_.clear(); + + readImpl(is); + computeUnits(); + if (convert_to_SI) { + convertToSI(); + } +} + +//--------------------------------------------------------------------------- +void EclipseGridParser::readImpl(istream& is) +//--------------------------------------------------------------------------- +{ + if (!is) { + cerr << "Could not read given input stream." << endl; + throw exception(); + } + + // Make temporary maps that will at the end be swapped with the + // member maps + // NOTE: Above is no longer true, for easier implementation of + // the INCLUDE keyword. We lose the strong exception guarantee, + // though (of course retaining the basic guarantee). + map >& intmap = integer_field_map_; + map >& floatmap = floating_field_map_; + map >& specialmap = special_field_map_; + + // Actually read the data + is >> ignoreWhitespace; + while (!is.eof()) { + string keyword = readKeyword(is); +#ifdef VERBOSE + cout << "Keyword found: " << keyword << endl; +#endif + FieldType type = classifyKeyword(keyword); + switch (type) { + case Integer: + readVectorData(is, intmap[keyword]); + break; + case FloatingPoint: + readVectorData(is, floatmap[keyword]); + break; + case SpecialField: { + boost::shared_ptr sb_ptr = createSpecialField(is, keyword); + if (sb_ptr) { + specialmap[keyword] = sb_ptr; + } else { + THROW("Could not create field " << keyword); + } + break; + } + case IgnoreWithData: { + ignored_fields_.insert(keyword); + is >> ignoreSlashLine; +#ifdef VERBOSE + cout << "(ignored)" << endl; +#endif + break; + } + case IgnoreNoData: { + ignored_fields_.insert(keyword); + is >> ignoreLine; +#ifdef VERBOSE + cout << "(ignored)" << endl; +#endif + break; + } + case Include: { + string include_filename = readString(is); + if (!directory_.empty()) { + include_filename = directory_ + '/' + include_filename; + } + ifstream include_is(include_filename.c_str()); + if (!include_is) { + THROW("Unable to open INCLUDEd file " << include_filename); + } + readImpl(include_is); + is >> ignoreSlashLine; + break; + } + case Unknown: + default: + cerr << "Keyword " << keyword << " not recognized." << endl; + throw exception(); + } + is >> ignoreWhitespace; + } + +#define VERBOSE_LIST_FIELDS 0 +#if VERBOSE_LIST_FIELDS + std::cout << "\nInteger fields:\n"; + for (std::map >::iterator + i = intmap.begin(); i != intmap.end(); ++i) + std::cout << '\t' << i->first << '\n'; + + std::cout << "\nFloat fields:\n"; + for (std::map >::iterator + i = floatmap.begin(); i != floatmap.end(); ++i) + std::cout << '\t' << i->first << '\n'; + + std::cout << "\nSpecial fields:\n"; + for (std::map >::iterator + i = specialmap.begin(); i != specialmap.end(); ++i) + std::cout << '\t' << i->first << '\n'; +#endif +} + + + +//--------------------------------------------------------------------------- +void EclipseGridParser::convertToSI() +//--------------------------------------------------------------------------- +{ + // Convert all special fields. + typedef std::map >::iterator SpecialIt; + for (SpecialIt i = special_field_map_.begin(); i != special_field_map_.end(); ++i) { + i->second->convertToSI(units_); + } + + // Convert all floating point fields. + typedef std::map >::iterator FloatIt; + for (FloatIt i = floating_field_map_.begin(); i != floating_field_map_.end(); ++i) { + const std::string& key = i->first; + std::vector& field = i->second; + // Find the right unit. + double unit = 1e100; + if (key == "COORD" || key == "ZCORN") { + unit = units_.length; + } else if (key == "PERMX" || key == "PERMY" || key == "PERMZ" || + key == "PERMXX" || key == "PERMYY" || key == "PERMZZ" || + key == "PERMXY" || key == "PERMYZ" || key == "PERMZX") { + unit = units_.permeability; + } else if (key == "PORO" || key == "BULKMOD" || key == "YOUNGMOD" || + key == "LAMEMOD" || key == "SHEARMOD" || key == "POISSONMOD" || + key == "PWAVEMOD" || key == "MULTPV" || key == "PWAVEMOD" || + key == "SGAS" || key == "SWAT") { + unit = 1.0; + } else if (key == "PRESSURE") { + unit = units_.pressure; + } else { + THROW("Units for field " << key << " not specified. Cannon convert to SI."); + } + // Convert. + for (int i = 0; i < int(field.size()); ++i) { + field[i] = Dune::unit::convert::from(field[i], unit); + } + } + + // Set all units to one. + units_.setToOne(); +} + + +/// Returns true if the given keyword corresponds to a field that +/// was found in the file. +//--------------------------------------------------------------------------- +bool EclipseGridParser::hasField(const string& keyword) const +//--------------------------------------------------------------------------- +{ + string ukey = upcase(keyword); + return integer_field_map_.count(ukey) || floating_field_map_.count(ukey) || + special_field_map_.count(ukey) || ignored_fields_.count(ukey); +} + + +/// Returns true if all the given keywords correspond to fields +/// that were found in the file. +//--------------------------------------------------------------------------- +bool EclipseGridParser::hasFields(const vector& keywords) const +//--------------------------------------------------------------------------- +{ + int num_keywords = keywords.size(); + for (int i = 0; i < num_keywords; ++i) { + if (!hasField(keywords[i])) { + return false; + } + } + return true; +} + +//--------------------------------------------------------------------------- +vector EclipseGridParser::fieldNames() const +//--------------------------------------------------------------------------- +{ + vector names; + names.reserve(integer_field_map_.size() + + floating_field_map_.size() + + special_field_map_.size() + + ignored_fields_.size()); + { + map >::const_iterator it = integer_field_map_.begin(); + for (; it != integer_field_map_.end(); ++it) { + names.push_back(it->first); + } + } + { + map >::const_iterator it = floating_field_map_.begin(); + for (; it != floating_field_map_.end(); ++it) { + names.push_back(it->first); + } + } + { + map >::const_iterator it = special_field_map_.begin(); + for (; it != special_field_map_.end(); ++it) { + names.push_back(it->first); + } + } + { + set::const_iterator it = ignored_fields_.begin(); + for (; it != ignored_fields_.end(); ++it) { + names.push_back(*it); + } + } + return names; +} + +//--------------------------------------------------------------------------- +const std::vector& EclipseGridParser::getIntegerValue(const std::string& keyword) const +//--------------------------------------------------------------------------- +{ + if (keyword == "SPECGRID") { + cerr << "\nERROR. Interface has changed!\n" + << "const vector& dim = parser.getIntegerValue(""SPECGRID"") is deprecated.\n" + << "Use:\n" + << "const SPECGRID& specgrid = parser.getSPECGRID();\n" + << "const vector& dim = specgrid.dimensions;\n\n"; + throw exception(); + } + + map >::const_iterator it + = integer_field_map_.find(keyword); + if (it == integer_field_map_.end()) { + return empty_integer_field_; + } else { + return it->second; + } +} + +//--------------------------------------------------------------------------- +const std::vector& EclipseGridParser::getFloatingPointValue(const std::string& keyword) const +//--------------------------------------------------------------------------- +{ + map >::const_iterator it + = floating_field_map_.find(keyword); + if (it == floating_field_map_.end()) { + return empty_floating_field_; + } else { + return it->second; + } +} + + +//--------------------------------------------------------------------------- +const boost::shared_ptr EclipseGridParser::getSpecialValue(const std::string& keyword) const +//--------------------------------------------------------------------------- +{ + map >::const_iterator it = special_field_map_.find(keyword); + if (it == special_field_map_.end()) { + THROW("No such field: " << keyword); + } else { + return it->second; + } +} + +//--------------------------------------------------------------------------- +boost::shared_ptr +EclipseGridParser::createSpecialField(std::istream& is, + const std::string& fieldname) +//--------------------------------------------------------------------------- +{ + string ukey = upcase(fieldname); + boost::shared_ptr spec_ptr + = Factory::createObject(fieldname); + spec_ptr->read(is); + return spec_ptr; +} + +//--------------------------------------------------------------------------- +void EclipseGridParser::setIntegerField(const std::string& keyword, + const std::vector& field) +//--------------------------------------------------------------------------- +{ + integer_field_map_[keyword] = field; +} + +//--------------------------------------------------------------------------- +void EclipseGridParser::setFloatingPointField(const std::string& keyword, + const std::vector& field) +//--------------------------------------------------------------------------- +{ + floating_field_map_[keyword] = field; +} + +//--------------------------------------------------------------------------- +void EclipseGridParser::setSpecialField(const std::string& keyword, + boost::shared_ptr field) +//--------------------------------------------------------------------------- +{ + special_field_map_[keyword] = field; +} + +//--------------------------------------------------------------------------- +const EclipseUnits& EclipseGridParser::units() const +//--------------------------------------------------------------------------- +{ + return units_; +} + +//--------------------------------------------------------------------------- +void EclipseGridParser::computeUnits() +//--------------------------------------------------------------------------- +{ + // Decide unit family. + enum EclipseUnitFamily { Metric = 0, Field = 1, Lab = 2, Pvtm = 3 }; + EclipseUnitFamily unit_family = Metric; // The default. + if (hasField("FIELD")) unit_family = Field; + if (hasField("LAB")) unit_family = Lab; + if (hasField("PVT-M")) unit_family = Pvtm; + + // Set units. + using namespace prefix; + using namespace unit; + switch (unit_family) { + case Metric: + units_.length = meter; + units_.time = day; + units_.density = kilogram/cubic(meter); + units_.pressure = barsa; + units_.compressibility = 1.0/barsa; + units_.viscosity = centi*Poise; + units_.permeability = milli*darcy; + units_.liqvol_s = cubic(meter); + units_.liqvol_r = cubic(meter); + units_.gasvol_s = cubic(meter); + units_.gasvol_r = cubic(meter); + units_.transmissibility = centi*Poise * cubic(meter) / (day * barsa); + break; + case Field: + units_.length = feet; + units_.time = day; + units_.density = pound/cubic(feet); + units_.pressure = psia; + units_.compressibility = 1.0/psia; + units_.viscosity = centi*Poise; + units_.permeability = milli*darcy; + units_.liqvol_s = stb; + units_.liqvol_r = stb; + units_.gasvol_s = 1000*cubic(feet); // Prefix 'M' is 1000 + units_.gasvol_r = stb; + units_.transmissibility = centi*Poise * stb / (day * psia); + break; + case Lab: + THROW("Unhandled unit family " << unit_family); + break; + case Pvtm: + THROW("Unhandled unit family " << unit_family); + break; + default: + THROW("Unknown unit family " << unit_family); + } +} + +} // namespace Dune diff --git a/dune/common/EclipseGridParser.hpp b/dune/common/EclipseGridParser.hpp new file mode 100644 index 00000000..2b84a657 --- /dev/null +++ b/dune/common/EclipseGridParser.hpp @@ -0,0 +1,191 @@ +//=========================================================================== +// +// File: EclipseGridParser.h +// +// Created: Wed Dec 5 17:05:13 2007 +// +// Author: Atgeirr F Rasmussen +// +// $Date$ +// +// Revision: $Id: EclipseGridParser.h,v 1.3 2008/08/18 14:16:13 atgeirr Exp $ +// +//=========================================================================== + +/* +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 . +*/ + +#ifndef SINTEF_ECLIPSEGRIDPARSER_HEADER +#define SINTEF_ECLIPSEGRIDPARSER_HEADER + +#include +#include +#include +#include +#include +#include +#include "SpecialEclipseFields.hpp" +#include "EclipseUnits.hpp" +#include + +namespace Dune +{ + +/** + @brief A class for reading and parsing all fields of an eclipse file. + + This object is constructed using an Eclipse .grdecl-file. All data + fields are extracted upon construction and written to vector data + structures, which can then be read out in O(1) time afterwards via + convenience functions. + + There is also a convenience function to easily check which fields + were successfully parsed. + + @author Atgeirr F. Rasmussen + @date 2007/12/06 13:00:03 + +*/ + +class EclipseGridParser +{ +public: + /// Default constructor. + EclipseGridParser(); + /// Constructor taking an eclipse filename. Unless the second + /// argument 'convert_to_SI' is false, all fields will be + /// converted to SI units. + explicit EclipseGridParser(const std::string& filename, bool convert_to_SI = true); + + /// Read the given stream, overwriting any previous data. Unless + /// the second argument 'convert_to_SI' is false, all fields will + /// be converted to SI units. + void read(std::istream& is, bool convert_to_SI = true); + + /// Convert all data to SI units, according to unit category as + /// specified in the eclipse file (METRIC, FIELD etc.). + /// It is unnecessary to call this if constructed with + /// 'convert_to_SI' equal to true, but it is not an error. + void convertToSI(); + + /// Returns true if the given keyword corresponds to a field that + /// was found in the file. + bool hasField(const std::string& keyword) const; + /// Returns true if all the given keywords correspond to fields + /// that were found in the file. + bool hasFields(const std::vector& keywords) const; + /// The keywords/fields found in the file. + std::vector fieldNames() const; + + /// Returns a reference to a vector containing the values + /// corresponding to the given integer keyword. + const std::vector& getIntegerValue(const std::string& keyword) const; + + /// Returns a reference to a vector containing the values + /// corresponding to the given floating-point keyword. + const std::vector& getFloatingPointValue(const std::string& keyword) const; + + /// Returns a reference to a vector containing pointers to the values + /// corresponding to the given keyword when the values are not only integers + /// or floats. + const boost::shared_ptr getSpecialValue(const std::string& keyword) const; + + // This macro implements support for a special field keyword. It requires that a subclass + // of SpecialBase exists, that has the same name as the keyword. + // After using SPECIAL_FIELD(KEYWORD), the public member getKEYWORD will be available. +#define SPECIAL_FIELD(keyword) \ +private: \ + struct X##keyword { X##keyword() { Factory::addCreator(#keyword); } }; \ + X##keyword x##keyword; \ +public: \ + const keyword& get##keyword() const \ + { return dynamic_cast(*getSpecialValue(#keyword)); } + + // Support for special fields. + SPECIAL_FIELD(SPECGRID); + SPECIAL_FIELD(FAULTS); + SPECIAL_FIELD(MULTFLT); + SPECIAL_FIELD(TITLE); + SPECIAL_FIELD(START); + SPECIAL_FIELD(DATES); + SPECIAL_FIELD(DENSITY); + SPECIAL_FIELD(PVDG); + SPECIAL_FIELD(PVDO); + SPECIAL_FIELD(PVTG); + SPECIAL_FIELD(PVTO); + SPECIAL_FIELD(PVTW); + SPECIAL_FIELD(SGOF); + SPECIAL_FIELD(SWOF); + SPECIAL_FIELD(ROCK); + SPECIAL_FIELD(ROCKTAB); + SPECIAL_FIELD(WELSPECS); + SPECIAL_FIELD(COMPDAT); + SPECIAL_FIELD(WCONINJE); + SPECIAL_FIELD(WCONPROD); + SPECIAL_FIELD(WELTARG); + SPECIAL_FIELD(EQUIL); + SPECIAL_FIELD(PVCDO); + + // The following fields only have a dummy implementation + // that allows us to ignore them. + SPECIAL_FIELD(SWFN); + SPECIAL_FIELD(SOF2); + SPECIAL_FIELD(TUNING); + +#undef SPECIAL_FIELD + + + /// Sets an integer field to have a particular value. + void setIntegerField(const std::string& keyword, const std::vector& field); + + /// Sets a floating point field to have a particular value. + void setFloatingPointField(const std::string& keyword, const std::vector& field); + + /// Sets a special field to have a particular value. + void setSpecialField(const std::string& keyword, boost::shared_ptr field); + + /// Compute the units used by the deck, depending on the presence + /// of keywords such as METRIC, FIELD etc. It is an error to call + /// this after conversion to SI has taken place. + void computeUnits(); + + /// The units specified by the eclipse file read. + const EclipseUnits& units() const; + +private: + boost::shared_ptr createSpecialField(std::istream& is, const std::string& fieldname); + void readImpl(std::istream& is); + + + std::string directory_; + std::map > integer_field_map_; + std::map > floating_field_map_; + std::map > special_field_map_; + std::set ignored_fields_; + std::vector empty_integer_field_; + std::vector empty_floating_field_; + boost::shared_ptr empty_special_field_; + EclipseUnits units_; +}; + + +} // namespace Dune + +#endif // SINTEF_ECLIPSEGRIDPARSER_HEADER diff --git a/dune/common/EclipseGridParserHelpers.hpp b/dune/common/EclipseGridParserHelpers.hpp new file mode 100644 index 00000000..fdf056dc --- /dev/null +++ b/dune/common/EclipseGridParserHelpers.hpp @@ -0,0 +1,469 @@ +//=========================================================================== +// +// File: EclipseGridParserHelpers.hpp +// +// Created: Tue Dec 22 11:35:32 2009 +// +// Author(s): Atgeirr F Rasmussen +// Bård Skaflestad +// +// $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 . +*/ + +#ifndef OPENRS_ECLIPSEGRIDPARSERHELPERS_HEADER +#define OPENRS_ECLIPSEGRIDPARSERHELPERS_HEADER + +#include +#include +#include +#include +#include +#include "linInt.hpp" +#include + +namespace Dune +{ + +namespace +{ + + inline std::istream& ignoreLine(std::istream& is) + { + is.ignore(std::numeric_limits::max(), '\n'); + return is; + } + + inline std::istream& ignoreSlashLine(std::istream& is) + { + is.ignore(std::numeric_limits::max(), '/'); + is.ignore(std::numeric_limits::max(), '\n'); + return is; + } + + inline std::istream& ignoreWhitespace(std::istream& is) + { + // Getting the character type facet for is() + // We use the classic (i.e. C) locale. + const std::ctype& ct = std::use_facet< std::ctype >(std::locale::classic()); + char c; + while (is.get(c)) { + if (!ct.is(std::ctype_base::space, c)) { + is.putback(c); + break; + } + } + return is; + } + + inline std::string upcase(const std::string& s) + { + std::string us(s); + // Getting the character type facet for toupper(). + // We use the classic (i.e. C) locale. + const std::ctype& ct = std::use_facet< std::ctype >(std::locale::classic()); + for (int i = 0; i < int(s.size()); ++i) { + us[i] = ct.toupper(s[i]); + } + return us; + } + + inline std::string readKeyword(std::istream& is) + { + std::string keyword_candidate; + while (!is.eof()) { + is >> keyword_candidate; + if(keyword_candidate.find("--") == 0) { + is >> ignoreLine; // This line is a comment + } else { + return upcase(keyword_candidate); + } + } + return "CONTINUE"; // Last line in included file is a comment + } + + inline std::string readString(std::istream& is) + { + std::string string_candidate; + is >> string_candidate; + const char quote('\''); + int beg = string_candidate[0] == quote ? 1 : 0; + int len = string_candidate[0] == quote ? string_candidate.size() - 2 : string_candidate.size(); + return string_candidate.substr(beg, len); + } + + // Reads data until '/' or an error is encountered. + template + inline void readVectorData(std::istream& is, std::vector& data) + { + data.clear(); + while (is) { + T candidate; + is >> candidate; + if (is.rdstate() & std::ios::failbit) { + is.clear(is.rdstate() & ~std::ios::failbit); + is >> ignoreWhitespace; + char dummy; + is >> dummy; + if (dummy == '/') { + is >> ignoreLine; + break; + } else if (dummy == '-') { // "comment test" + is >> ignoreLine; // This line is a comment + } else { + THROW("Encountered format error while reading data values. Value = " << dummy); + } + } else { + if (is.peek() == int('*')) { + is.ignore(); // ignore the '*' + int multiplier = int(candidate); + is >> candidate; + data.insert(data.end(), multiplier, candidate); + } else { + data.push_back(candidate); + } + } + } + if (!is) { + THROW("Encountered error while reading data values."); + } + } + + + // Reads data items of type T. Not more than 'max_values' items. + // Asterisks may be used to signify 'repeat counts'. 5*3.14 will + // insert 3.14 five times. Asterisk followed by a space is used to + // signify default values. n* will default n consecutive quantities. + template + inline int readDefaultedVectorData(std::istream& is, Vec& data, int max_values) + { + ASSERT(int(data.size()) >= max_values); + const std::ctype& ct = std::use_facet< std::ctype >(std::locale::classic()); + int num_values = 0; + while (is) { + typename Vec::value_type candidate; + is >> candidate; + if (is.rdstate() & std::ios::failbit) { + is.clear(is.rdstate() & ~std::ios::failbit); + std::string dummy; + is >> dummy; + if (dummy == "/") { + is >> ignoreLine; + break; + } else if (dummy[0] == '-') { // "comment test" + is >> ignoreLine; // This line is a comment + } else { + THROW("Encountered format error while reading data values. Value = " << dummy); + } + } else { + if (is.peek() == int('*')) { + is.ignore(); // ignore the '*' + int multiplier = (int)candidate; + if (ct.is(std::ctype_base::space, is.peek())) { + num_values += multiplier; // Use default value(s) + } else { + is >> candidate; // Use candidate 'multipler' times + for (int i=0; i= max_values) { + //is >> ignoreLine; + break; + } + } + if (!is) { + THROW("Encountered error while reading data values."); + } + return num_values; + } + + + // Keywords SGOF and SWOF. Reads data until '/' or an error is encountered. + // Default values represented by 1* is replaced by -1. Use linear interpolation + // outside this function to replace -1. + template + inline void readRelPermTable(std::istream& is, std::vector& data) + { + data.clear(); + while (is) { + T candidate; + is >> candidate; + if (is.rdstate() & std::ios::failbit) { + is.clear(is.rdstate() & ~std::ios::failbit); + std::string dummy; + is >> dummy; + if (dummy == "/") { + is >> ignoreLine; + break; + } else if (dummy[0] == '-') { // "comment test" + is >> ignoreLine; // This line is a comment + } else { + THROW("Encountered format error while reading data values. Value = " << dummy); + } + } else { + if (is.peek() == int('*')) { + is.ignore(); // ignore the '*' + ASSERT(int(candidate) == 1); + data.push_back(-1); // Set new flag for interpolation. + } else { + data.push_back(candidate); + } + } + } + if (!is) { + THROW("Encountered error while reading data values."); + } + } + + + // Returns month number 1-12. Returns 0 if illegal month name. + inline int getMonthNumber(const std::string& month_name) + { + const int num_months = 12; + std::string months[num_months] = {"JAN", "FEB", "MAR", "APR", "MAY", "JUN", + "JLY", "AUG", "SEP", "OCT", "NOV", "DEC"}; + if (month_name == "JUL") { + return 7; // "JUL" is an acceptable alternative to 'JLY' + } + int m = 0; + for (int i=0; i> ignoreLine; // This line is a comment + } + int day, year; + std::string month_name; + is >> day; + month_name = readString(is); + is >> year; + ignoreSlashLine(is); + int month = getMonthNumber(month_name); + return boost::gregorian::date(boost::gregorian::greg_year(year), + boost::gregorian::greg_month(month), + boost::gregorian::greg_day(day)); + } + + // Get first character after whitespace and comments, and decide + // next action. + // NB! Will fail for negative number in column one. + inline int next_action(std::istream& is) + { + int task(0); // 0:continue 1:return 2:throw + const std::ctype& ct = + std::use_facet< std::ctype >(std::locale::classic()); + + while (!is.eof()) { + is >> ignoreWhitespace; + char c; + is.get(c); + if (is.eof()) { + return task; + } + is.putback(c); + if (ct.is(std::ctype_base::digit, c) || c== '.') { + task = 0; // Decimal digit. Read more data. + break; + } + if (ct.is(std::ctype_base::alpha, c)) { + task = 1; // Alphabetic char. Read next keyword. + break; + } + if (c == '-') { + is >> ignoreLine; // This line is a comment + } else { + task = 2; // Read error. Unexpected character. + break; + } + } + return task; + } + + // Reads keywords PVTG, PVTO + typedef std::vector > > table_t; + inline void readPvtTable(std::istream& is, table_t& pvt_table, + const std::string& field_name) + { + const std::ctype& ct = + std::use_facet< std::ctype >(std::locale::classic()); + std::vector record; + std::vector > table; + while (!is.eof()) { + record.clear(); + readVectorData(is, record); + table.push_back(record); + while (!is.eof()) { + is >> ignoreWhitespace; + char c; + is.get(c); + if (is.eof()) { + // Reached end of file, we should have pushed + // the last table, and emptied it. If not, + // we have an error. + if (!table.empty()) { + THROW("Reached EOF while still building PVT table. Missing end-of-table (slash)?"); + } + return; + } + is.putback(c); + if (ct.is(std::ctype_base::digit, c) || c== '.') { + break; // Decimal digit. Read more records. + } + if (ct.is(std::ctype_base::alpha, c)) { + return; // Alphabetic char. Read next keyword. + } + if (c == '-') { + is >> ignoreLine; // This line is a comment + continue; + } + if (c == '/') { + is >> ignoreLine; + pvt_table.push_back(table); + table.clear(); + } else { + std::ostringstream oss; + oss << "Error reading " << field_name + << ". Next character is " << (char)is.peek(); + THROW(oss.str()); + } + } + } + } + + // Reads keywords PVDG, PVDO, ROCKTAB + inline void readPvdTable(std::istream& is, table_t& pvd_table, + const std::string& field_name, int ncol) + { + std::vector record; + std::vector > table(ncol); + while (!is.eof()) { + record.clear(); + readVectorData(is, record); + const int rec_size = record.size()/ncol; + for (int k=0; k >& table, int ncol) + { + const int sz = table[0].size(); + for (int k=1; k indx; + std::vector x; + for (int i=0; i xv, yv; + for (int i=0; i record; + std::vector > table(ncol); + while (!is.eof()) { + record.clear(); + readRelPermTable(is, record); + const int rec_size = record.size()/ncol; + for (int k=0; k. +*/ + +#ifndef OPM_ECLIPSEUNITS_HEADER_INCLUDED +#define OPM_ECLIPSEUNITS_HEADER_INCLUDED + +namespace Dune +{ + struct EclipseUnits + { + double length; + double time; + double density; + double pressure; + double compressibility; + double viscosity; + double permeability; + double liqvol_s; + double liqvol_r; + double gasvol_s; + double gasvol_r; + double transmissibility; + + void setToOne() + { + length = 1.0; + time = 1.0; + density = 1.0; + pressure = 1.0; + compressibility = 1.0; + viscosity = 1.0; + permeability = 1.0; + liqvol_s = 1.0; + liqvol_r = 1.0; + gasvol_s = 1.0; + gasvol_r = 1.0; + transmissibility = 1.0; + } + }; + + +} // namespace Dune + + +#endif // OPM_ECLIPSEUNITS_HEADER_INCLUDED diff --git a/dune/common/ErrorMacros.hpp b/dune/common/ErrorMacros.hpp new file mode 100644 index 00000000..10992251 --- /dev/null +++ b/dune/common/ErrorMacros.hpp @@ -0,0 +1,107 @@ +//=========================================================================== +// +// File: ErrorMacros.hpp +// +// Created: Tue May 14 12:22:16 2002 +// +// Author(s): Atgeirr F Rasmussen +// +// $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 . +*/ + +#ifndef OPENRS_ERRORMACROS_HEADER +#define OPENRS_ERRORMACROS_HEADER + + +/// Error macros. In order to use some of them, you must also +/// include or , so they are included by +/// this file. The compile defines NDEBUG and NVERBOSE control +/// the behaviour of these macros. + +#ifndef NVERBOSE + #include +#endif +#include + +/// Usage: REPORT; +/// Usage: MESSAGE("Message string."); +#ifdef NVERBOSE // Not verbose mode +# ifndef REPORT +# define REPORT +# endif +# ifndef MESSAGE +# define MESSAGE(x) +# endif +# ifndef MESSAGE_IF +# define MESSAGE_IF(cond, m) +# endif +#else // Verbose mode +# ifndef REPORT +# define REPORT std::cerr << "\nIn file " << __FILE__ << ", line " << __LINE__ << std::endl +# endif +# ifndef MESSAGE +# define MESSAGE(x) std::cerr << "\nIn file " << __FILE__ << ", line " << __LINE__ << ": " << x << std::endl +# endif +# ifndef MESSAGE_IF +# define MESSAGE_IF(cond, m) do {if(cond) MESSAGE(m);} while(0) +# endif +#endif + + +/// Usage: THROW("Error message string."); +#ifndef THROW +# define THROW(x) do { MESSAGE(x); throw std::exception(); } while(0) +#endif + +#define ALWAYS_ERROR_IF(condition, message) do {if(condition){ THROW(message);}} while(0) + +/// Usage: ASSERT(condition) +/// Usage: ASSERT2(condition, "Error message string.") +/// Usage: DEBUG_ERROR_IF(condition, "Error message string."); +#ifdef NDEBUG // Not in debug mode +# ifndef ASSERT +# define ASSERT(x) +# endif +# ifndef ASSERT2 +# define ASSERT2(cond, x) +# endif +# ifndef DEBUG_ERROR_IF +# define DEBUG_ERROR_IF(cond, x) +# endif +#else // Debug mode +# ifndef ASSERT +# define ASSERT(cond) if (!(cond)) THROW("Assertion \'" #cond "\' failed.") +# endif +# ifndef ASSERT2 +# define ASSERT2(cond, x) do { if (!(cond)) THROW(x);} while(0) +# endif +# ifndef DEBUG_ERROR_IF +# define DEBUG_ERROR_IF(cond, x) do { if (cond) THROW(x); } while(0) +# endif +#endif + + +#endif // OPENRS_ERRORMACROS_HEADER diff --git a/dune/common/Factory.hpp b/dune/common/Factory.hpp new file mode 100644 index 00000000..949e4ea7 --- /dev/null +++ b/dune/common/Factory.hpp @@ -0,0 +1,145 @@ +//=========================================================================== +// +// File: Factory.hpp +// +// Created: Mon Nov 6 10:00:43 2000 +// +// Author(s): Atgeirr F Rasmussen +// +// $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 . +*/ + +#ifndef OPENRS_FACTORY_HEADER +#define OPENRS_FACTORY_HEADER + + +#include +#include + +namespace Dune +{ + + /** This is an object factory for creating objects of some type + * requested by the user, with a shared base class. The user + * need only interact with the factory through the static + * template member addCreator() and the static member function + * createObject(). + */ + template + class Factory + { + public: + /// The type of pointer returned by createObject(). + typedef boost::shared_ptr ProductPtr; + + /// Creates a new object of the class associated with the given type string, + /// and returns a pointer to it. + /// \param type the type string of the class that the user wants to have + /// constructed. + /// \return (smart) pointer to the created object. + static ProductPtr createObject(const std::string& type) + { + return instance().doCreateObject(type); + } + + /// Add a creator to the Factory. + /// After the call, the user may obtain new objects of the Derived type by + /// calling createObject() with the given type string as an argument. + /// \tparam Derived the class we want to add a creator for, must inherit + /// the class template parameter Base. + /// \param type the type string with which we want the Factory to associate + /// the class Derived. + template + static void addCreator(const std::string& type) + { + instance().doAddCreator(type); + } + + private: + // The method that implements the singleton pattern, + // using the Meyers singleton technique. + static Factory& instance() + { + static Factory singleton; + return singleton; + } + + // Private constructor, to keep users from creating a Factory. + Factory() + { + } + + // Abstract base class for Creators. + class Creator + { + public: + virtual ProductPtr create() = 0; + virtual ~Creator() {} + }; + + /// This is the concrete Creator subclass for generating Derived objects. + template + class ConcreteCreator : public Creator + { + public: + virtual ProductPtr create() + { + return ProductPtr(new Derived); + } + }; + + typedef boost::shared_ptr CreatorPtr; + typedef std::map CreatorMap; + // This map contains the whole factory, i.e. all the Creators. + CreatorMap string_to_creator_; + + // Actually creates the product object. + ProductPtr doCreateObject(const std::string& type) + { + typename CreatorMap::iterator it; + it = string_to_creator_.find(type); + if (it == string_to_creator_.end()) { + THROW("Creator type " << type + << " is not registered in the factory."); + } + return it->second->create(); + } + + // Actually adds the creator. + template + void doAddCreator(const std::string& type) + { + boost::shared_ptr c(new ConcreteCreator); + string_to_creator_[type] = c; + } + }; + + +} // namespace Dune + +#endif // OPENRS_FACTORY_HEADER + + diff --git a/dune/common/Makefile.am b/dune/common/Makefile.am new file mode 100644 index 00000000..aa3db8e7 --- /dev/null +++ b/dune/common/Makefile.am @@ -0,0 +1,38 @@ +# $Date$ +# $Revision$ + +SUBDIRS = . param test + +common_HEADERS = \ + Average.hpp \ + EclipseGridInspector.hpp \ + EclipseGridParser.hpp \ + EclipseGridParserHelpers.hpp \ + ErrorMacros.hpp \ + Factory.hpp \ + linInt.hpp \ + MonotCubicInterpolator.hpp \ + RootFinders.hpp \ + SparseTable.hpp \ + SparseVector.hpp \ + SpecialEclipseFields.hpp \ + StopWatch.hpp \ + Units.hpp + +commondir = $(includedir)/dune/common + +noinst_LTLIBRARIES = libcommon.la + +libcommon_la_SOURCES = \ + EclipseGridInspector.cpp \ + EclipseGridParser.cpp \ + MonotCubicInterpolator.cpp \ + StopWatch.cpp + + +libcommon_la_CPPFLAGS = $(DUNE_CPPFLAGS) $(BOOST_CPPFLAGS) +libcommon_la_LDFLAGS = $(DUNE_LDFLAGS) $(BOOST_LDFLAGS) \ + $(BOOST_DATE_TIME_LIB) + + +include $(top_srcdir)/am/global-rules diff --git a/dune/common/MonotCubicInterpolator.cpp b/dune/common/MonotCubicInterpolator.cpp new file mode 100644 index 00000000..8a20ecf8 --- /dev/null +++ b/dune/common/MonotCubicInterpolator.cpp @@ -0,0 +1,726 @@ +/* + MonotCubicInterpolator + Copyright (C) 2006 Statoil ASA + + This program 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 2 + of the License, or (at your option) any later version. + + This program 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 this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +/** + @file MonotCubicInterpolator.C + @brief Represents one dimensional function f with single valued argument x + + Class to represent a one-dimensional function with single-valued + argument. Cubic interpolation, preserving the monotonicity of the + discrete known function values + + @see MonotCubicInterpolator.h for further documentation. + +*/ + + +#include "MonotCubicInterpolator.hpp" +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +/* + + SOME DISCUSSION ON DATA STORAGE: + + Internal data structure of points and values: + + vector(s): + - Easier coding + - Faster vector operations when setting up derivatives. + - sorting slightly more complex. + - insertion of further values bad. + + vector + - easy sorting + - code complexity almost as for map. + - insertion of additional values bad + + vector over struct datapoint { x, f, d} + - nice code + - not as sortable, insertion is cumbersome. + + ** This is used currently: ** + map one for (x,f) and one for (x,d) + - Naturally sorted on x-values (done by the map-construction) + - Slower to set up, awkward loop coding (?) + - easy to add more points. + - easier to just add code to linear interpolation code + - x-data is duplicated, but that memory waste is + unlikely to represent a serious issue. + + map > + - naturally couples x-value, f-value and d-value + - even more unreadable(??) code? + - higher skills needed by programmer. + + + MONOTONE CUBIC INTERPOLATION: + + It is a local algorithm. Adding one point only incur recomputation + of values in a neighbourhood of the point (in the interval getting + divided). + + NOTE: We do not currently make use of this local fact. Upon + insertion of a new data pair, everything is recomputed. Revisit + this when needed. + +*/ + + + + +MonotCubicInterpolator:: +MonotCubicInterpolator(const vector & x, const vector & f) { + if (x.size() != f.size()) { + throw("Unable to constuct MonotCubicInterpolator from vectors.") ; + } + + // Add the contents of the input vectors to our map of values. + vector::const_iterator posx, posf; + for (posx = x.begin(), posf = f.begin(); posx != x.end(); ++posx, ++posf) { + data[*posx] = *posf ; + } + + computeInternalFunctionData(); +} + + + +bool +MonotCubicInterpolator:: +read(const std::string & datafilename, int xColumn, int fColumn) +{ + data.clear() ; + ddata.clear() ; + + ifstream datafile_fs(datafilename.c_str()); + if (!datafile_fs) { + return false ; + } + + string linestring; + while (!datafile_fs.eof()) { + getline(datafile_fs, linestring); + + // Replace commas with space: + string::size_type pos = 0; + while ( (pos = linestring.find(",", pos)) != string::npos ) { + // cout << "Found comma at position " << pos << endl; + linestring.replace(pos, 1, " "); + pos++; + } + + stringstream strs(linestring); + int columnindex = 0; + std::vector value; + if (linestring.size() > 0 && linestring.at(0) != '#') { + while (!(strs.rdstate() & std::ios::failbit)) { + double tmp; + strs >> tmp; + value.push_back(tmp); + columnindex++; + } + } + if (columnindex >= (max(xColumn, fColumn))) { + data[value[xColumn-1]] = value[fColumn-1] ; + } + } + datafile_fs.close(); + + if (data.size() == 0) { + return false ; + } + + computeInternalFunctionData(); + return true ; +} + + +void +MonotCubicInterpolator:: +addPair(double newx, double newf) throw(const char*) { + if (isnan(newx) || isinf(newx) || isnan(newf) || isinf(newf)) { + throw("MonotCubicInterpolator: addPair() received inf/nan input."); + } + data[newx] = newf ; + + // In a critical application, we should only update the + // internal function data for the offended interval, + // not for all function values over again. + computeInternalFunctionData(); +} + + +double +MonotCubicInterpolator:: +evaluate(double x) const throw(const char*){ + + if (isnan(x) || isinf(x)) { + throw("MonotCubicInterpolator: evaluate() received inf/nan input."); + } + + // xf becomes the first (xdata,fdata) pair where xdata >= x + map::const_iterator xf_iterator = data.lower_bound(x); + + // First check if we must extrapolate: + if (xf_iterator == data.begin()) { + if (data.begin()->first == x) { // Just on the interval limit + return data.begin()->second; + } + else { + // Constant extrapolation (!!) + return data.begin()->second; + } + } + if (xf_iterator == data.end()) { + // Constant extrapolation (!!) + return data.rbegin()->second; + } + + + // Ok, we have x_min < x < x_max + + pair xf2 = *xf_iterator; + pair xf1 = *(--xf_iterator); + // we now have: xf2.first > x + + // Linear interpolation if derivative data is not available: + if (ddata.size() != data.size()) { + double finterp = xf1.second + + (xf2.second - xf1.second) / (xf2.first - xf1.first) + * (x - xf1.first); + return finterp; + } + else { // Do Cubic Hermite spline + double t = (x - xf1.first)/(xf2.first - xf1.first); // t \in [0,1] + double h = xf2.first - xf1.first; + double finterp + = xf1.second * H00(t) + + ddata[xf1.first] * H10(t) * h + + xf2.second * H01(t) + + ddata[xf2.first] * H11(t) * h ; + return finterp; + } + +} + + +// double +// MonotCubicInterpolator:: +// evaluate(double x, double& errorestimate_output) { +// cout << "Error: errorestimate not implemented" << endl; +// throw("error estimate not implemented"); +// return x; +// } + +vector +MonotCubicInterpolator:: +get_xVector() const +{ + map::const_iterator xf_iterator = data.begin(); + + vector outputvector; + outputvector.reserve(data.size()); + for (xf_iterator = data.begin(); xf_iterator != data.end(); ++xf_iterator) { + outputvector.push_back(xf_iterator->first); + } + return outputvector; +} + + +vector +MonotCubicInterpolator:: +get_fVector() const +{ + + map::const_iterator xf_iterator = data.begin(); + + vector outputvector; + outputvector.reserve(data.size()); + for (xf_iterator = data.begin(); xf_iterator != data.end(); ++xf_iterator) { + outputvector.push_back(xf_iterator->second); + } + return outputvector; +} + + + +string +MonotCubicInterpolator:: +toString() const +{ + const int precision = 20; + std::string dataString; + std::stringstream dataStringStream; + for (map::const_iterator it = data.begin(); + it != data.end(); ++it) { + dataStringStream << setprecision(precision) << it->first; + dataStringStream << '\t'; + dataStringStream << setprecision(precision) << it->second; + dataStringStream << '\n'; + } + dataStringStream << "Derivative values:" << endl; + for (map::const_iterator it = ddata.begin(); + it != ddata.end(); ++it) { + dataStringStream << setprecision(precision) << it->first; + dataStringStream << '\t'; + dataStringStream << setprecision(precision) << it->second; + dataStringStream << '\n'; + } + + return dataStringStream.str(); + +} + + +pair +MonotCubicInterpolator:: +getMissingX() const throw(const char*) +{ + if( data.size() < 2) { + throw("MonotCubicInterpolator::getMissingX() only one datapoint."); + } + + // Search for biggest difference value in function-datavalues: + + map::const_iterator maxfDiffPair_iterator = data.begin();; + double maxfDiffValue = 0; + + map::const_iterator xf_iterator; + map::const_iterator xf_next_iterator; + + for (xf_iterator = data.begin(), xf_next_iterator = ++(data.begin()); + xf_next_iterator != data.end(); + ++xf_iterator, ++xf_next_iterator) { + double absfDiff = fabs((double)(*xf_next_iterator).second + - (double)(*xf_iterator).second); + if (absfDiff > maxfDiffValue) { + maxfDiffPair_iterator = xf_iterator; + maxfDiffValue = absfDiff; + } + } + + double newXvalue = ((*maxfDiffPair_iterator).first + ((*(++maxfDiffPair_iterator)).first))/2; + return make_pair(newXvalue, maxfDiffValue); + +} + + + +pair +MonotCubicInterpolator:: +getMaximumF() const throw(const char*) { + if (data.size() <= 1) { + throw ("MonotCubicInterpolator::getMaximumF() empty data.") ; + } + if (strictlyIncreasing) + return *data.rbegin(); + else if (strictlyDecreasing) + return *data.begin(); + else { + pair maxf = *data.rbegin() ; + map::const_iterator xf_iterator; + for (xf_iterator = data.begin() ; xf_iterator != data.end(); ++xf_iterator) { + if (xf_iterator->second > maxf.second) { + maxf = *xf_iterator ; + } ; + } + return maxf ; + } +} + + +pair +MonotCubicInterpolator:: +getMinimumF() const throw(const char*) { + if (data.size() <= 1) { + throw ("MonotCubicInterpolator::getMinimumF() empty data.") ; + } + if (strictlyIncreasing) + return *data.begin(); + else if (strictlyDecreasing) { + return *data.rbegin(); + } + else { + pair minf = *data.rbegin() ; + map::const_iterator xf_iterator; + for (xf_iterator = data.begin() ; xf_iterator != data.end(); ++xf_iterator) { + if (xf_iterator->second < minf.second) { + minf = *xf_iterator ; + } ; + } + return minf ; + } +} + + +void +MonotCubicInterpolator:: +computeInternalFunctionData() const { + + /* The contents of this function is meaningless if there is only one datapoint */ + if (data.size() <= 1) { + return; + } + + /* We do not check the caching flag if we are instructed + to do this computation */ + + /* We compute monotoneness and directions by assuming + monotoneness, and setting to false if the function is not for + some value */ + + map::const_iterator xf_iterator; + map::const_iterator xf_next_iterator; + + + strictlyMonotone = true; // We assume this is true, and will set to false if not + monotone = true; + strictlyDecreasing = true; + decreasing = true; + strictlyIncreasing = true; + increasing = true; + + // Increasing or decreasing?? + xf_iterator = data.begin(); + xf_next_iterator = ++(data.begin()); + /* Cater for non-strictness, search for direction for monotoneness */ + while (xf_next_iterator != data.end() && + xf_iterator->second == xf_next_iterator->second) { + /* Ok, equal values, this is not strict. */ + strictlyMonotone = false; + strictlyIncreasing = false; + strictlyDecreasing = false; + + ++xf_iterator; + ++xf_next_iterator; + } + + + if (xf_next_iterator != data.end()) { + + if ( xf_iterator->second > xf_next_iterator->second) { + // Ok, decreasing, check monotoneness: + strictlyDecreasing = true;// if strictlyMonotone == false, this one should not be trusted anyway + decreasing = true; + strictlyIncreasing = false; + increasing = false; + while(++xf_iterator, ++xf_next_iterator != data.end()) { + if ((*xf_iterator).second < (*xf_next_iterator).second) { + monotone = false; + strictlyMonotone = false; + strictlyDecreasing = false; // meaningless now + break; // out of while loop + } + if ((*xf_iterator).second <= (*xf_next_iterator).second) { + strictlyMonotone = false; + strictlyDecreasing = false; // meaningless now + } + } + } + else if (xf_iterator->second < xf_next_iterator->second) { + // Ok, assume increasing, check monotoneness: + strictlyDecreasing = false; + strictlyIncreasing = true; + decreasing = false; + increasing = true; + while(++xf_iterator, ++xf_next_iterator != data.end()) { + if ((*xf_iterator).second > (*xf_next_iterator).second) { + monotone = false; + strictlyMonotone = false; + strictlyIncreasing = false; // meaningless now + break; // out of while loop + } + if ((*xf_iterator).second >= (*xf_next_iterator).second) { + strictlyMonotone = false; + strictlyIncreasing = false; // meaningless now + } + } + } + else { + // first two values must be equal if we get + // here, but that should have been taken care + // of by the while loop above. + throw("Programming logic error.") ; + } + + } + computeSimpleDerivatives(); + + + // If our input data is monotone, we can do monotone cubic + // interpolation, so adjust the derivatives if so. + // + // If input data is not monotone, we should not touch + // the derivatives, as this code should reduce to a + // standard cubic interpolation algorithm. + if (monotone) { + adjustDerivativesForMonotoneness(); + } + + strictlyMonotoneCached = true; + monotoneCached = true; +} + +// Checks if the function curve is flat (zero derivative) at the +// endpoints, chop off endpoint data points if that is the case. +// +// The notion of "flat" is determined by the input parameter "epsilon" +// Values whose difference are less than epsilon are regarded as equal. +// +// This is implemented to be able to obtain a strictly monotone +// curve from a data set that is strictly monotone except at the +// endpoints. +// +// Example: +// The data points +// (1,3), (2,3), (3,4), (4,5), (5,5), (6,5) +// will become +// (2,3), (3,4), (4,5) +// +// Assumes at least 3 datapoints. If less than three, this function is a noop. +void +MonotCubicInterpolator:: +chopFlatEndpoints(const double epsilon) { + + if (getSize() < 3) { + return; + } + + map::iterator xf_iterator; + map::iterator xf_next_iterator; + + // Clear flags: + strictlyMonotoneCached = false; + monotoneCached = false; + + // Chop left end: + xf_iterator = data.begin(); + xf_next_iterator = ++(data.begin()); + // Erase data points that are similar to its right value from the left end. + while ((xf_next_iterator != data.end()) && + (fabs(xf_iterator->second - xf_next_iterator->second) < epsilon )) { + xf_next_iterator++; + data.erase(xf_iterator); + xf_iterator++; + } + + xf_iterator = data.end(); + xf_iterator--; // (data.end() points beyond the last element) + xf_next_iterator = xf_iterator; + xf_next_iterator--; + // Erase data points that are similar to its left value from the right end. + while ((xf_next_iterator != data.begin()) && + (fabs(xf_iterator->second - xf_next_iterator->second) < epsilon )) { + xf_next_iterator--; + data.erase(xf_iterator); + xf_iterator--; + } + + // Finished chopping, so recompute function data: + computeInternalFunctionData(); +} + + +// +// If function is monotone, but not strictly monotone, +// this function will remove datapoints from intervals +// with zero derivative so that the curves become +// strictly monotone. +// +// Example +// The data points +// (1,2), (2,3), (3,4), (4,4), (5,5), (6,6) +// will become +// (1,2), (2,3), (3,4), (5,5), (6,6) +// +// Assumes at least two datapoints, if one or zero datapoint, this is a noop. +// +// +void +MonotCubicInterpolator:: +shrinkFlatAreas(const double epsilon) { + + if (getSize() < 2) { + return; + } + + map::iterator xf_iterator; + map::iterator xf_next_iterator; + + + // Nothing to do if we already are strictly monotone + if (isStrictlyMonotone()) { + return; + } + + // Refuse to change a curve that is not monotone. + if (!isMonotone()) { + return; + } + + // Clear flags, they are not to be trusted after we modify the + // data + strictlyMonotoneCached = false; + monotoneCached = false; + + // Iterate through data values, if two data pairs + // have equal values, delete one of the data pair. + // Do not trust the source code on which data point is being + // removed (x-values of equal y-points might be averaged in the future) + xf_iterator = data.begin(); + xf_next_iterator = ++(data.begin()); + + while (xf_next_iterator != data.end()) { + //cout << xf_iterator->first << "," << xf_iterator->second << " " << xf_next_iterator->first << "," << xf_next_iterator->second << "\n"; + if (fabs(xf_iterator->second - xf_next_iterator->second) < epsilon ) { + //cout << "erasing data pair" << xf_next_iterator->first << " " << xf_next_iterator->second << "\n"; + map ::iterator xf_tobedeleted_iterator = xf_next_iterator; + xf_next_iterator++; + data.erase(xf_tobedeleted_iterator); + } + else { + xf_iterator++; + xf_next_iterator++; + } + } + +} + + +void +MonotCubicInterpolator:: +computeSimpleDerivatives() const { + + ddata.clear(); + + // Do endpoints first: + map::const_iterator xf_prev_iterator; + map::const_iterator xf_iterator; + map::const_iterator xf_next_iterator; + double diff; + + // Leftmost interval: + xf_iterator = data.begin(); + xf_next_iterator = ++(data.begin()); + diff = + (xf_next_iterator->second - xf_iterator->second) / + (xf_next_iterator->first - xf_iterator->first); + ddata[xf_iterator->first] = diff ; + + // Rightmost interval: + xf_iterator = --(--(data.end())); + xf_next_iterator = --(data.end()); + diff = + (xf_next_iterator->second - xf_iterator->second) / + (xf_next_iterator->first - xf_iterator->first); + ddata[xf_next_iterator->first] = diff ; + + // If we have more than two intervals, loop over internal points: + if (data.size() > 2) { + + map::const_iterator intpoint; + for (intpoint = ++data.begin(); intpoint != --data.end(); ++intpoint) { + /* + diff = (f2 - f1)/(x2-x1)/w + (f3-f1)/(x3-x2)/2 + + average of the forward and backward difference. + Weights are equal, should we weigh with h_i? + */ + + map::const_iterator lastpoint = intpoint; --lastpoint; + map::const_iterator nextpoint = intpoint; ++nextpoint; + + diff = (nextpoint->second - intpoint->second)/ + (2*(nextpoint->first - intpoint->first)) + + + (intpoint->second - lastpoint->second) / + (2*(intpoint->first - lastpoint->first)); + + ddata[intpoint->first] = diff ; + } + } +} + + + +void +MonotCubicInterpolator:: +adjustDerivativesForMonotoneness() const { + map::const_iterator point, dpoint; + + /* Loop over all intervals, ie. loop over all points and look + at the interval to the right of the point */ + for (point = data.begin(), dpoint = ddata.begin(); + point != --data.end(); + ++point, ++dpoint) { + map::const_iterator nextpoint, nextdpoint; + nextpoint = point; ++nextpoint; + nextdpoint = dpoint; ++nextdpoint; + + double delta = + (nextpoint->second - point->second) / + (nextpoint->first - point->first); + if (fabs(delta) < 1e-14) { + ddata[point->first] = 0.0; + ddata[nextpoint->first] = 0.0; + } else { + double alpha = ddata[point->first] / delta; + double beta = ddata[nextpoint->first] / delta; + + if (! isMonotoneCoeff(alpha, beta)) { + double tau = 3/sqrt(alpha*alpha + beta*beta); + + ddata[point->first] = tau*alpha*delta; + ddata[nextpoint->first] = tau*beta*delta; + } + } + + + } + + +} + + + + +void +MonotCubicInterpolator:: +scaleData(double factor) { + map::iterator it , itd ; + if (data.size() == ddata.size()) { + for (it = data.begin() , itd = ddata.begin() ; it != data.end() ; ++it , ++itd) { + it->second *= factor ; + itd->second *= factor ; + } ; + } else { + for (it = data.begin() ; it != data.end() ; ++it ) { + it->second *= factor ; + } + } +} diff --git a/dune/common/MonotCubicInterpolator.hpp b/dune/common/MonotCubicInterpolator.hpp new file mode 100644 index 00000000..7e1c6710 --- /dev/null +++ b/dune/common/MonotCubicInterpolator.hpp @@ -0,0 +1,578 @@ +/* -*-C++-*- */ + +#ifndef _MONOTCUBICINTERPOLATOR_H +#define _MONOTCUBICINTERPOLATOR_H + +#include +#include +#include + +/* + MonotCubicInterpolator + Copyright (C) 2006 Statoil ASA + + This program 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 2 + of the License, or (at your option) any later version. + + This program 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 this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +/** + Class to represent a one-dimensional function f with single-valued + argument x. The function is represented by a table of function + values. Interpolation between table values is cubic and monotonicity + preserving if input values are monotonous. + + Outside x_min and x_max, the class will extrapolate using the + constant f(x_min) or f(x_max). + + Extra functionality: + - Can return (x_1+x_2)/2 where x_1 and x_2 are such that + abs(f(x_1) - f(x_2)) is maximized. This is used to determine where + one should calculate a new value for increased accuracy in the + current function + + Monotonicity preserving cubic interpolation algorithm is taken + from Fritsch and Carlson, "Monotone piecewise cubic interpolation", + SIAM J. Numer. Anal. 17, 238--246, no. 2, + + $Id$ + + Algorithm also described here: + http://en.wikipedia.org/wiki/Monotone_cubic_interpolation + + + @author HÃ¥vard Berland , December 2006 + @brief Represents one dimensional function f with single valued argument x that can be interpolated using monotone cubic interpolation + +*/ + +class MonotCubicInterpolator { + public: + + /** + @param datafilename A datafile with the x values and the corresponding f(x) values + + Accepts a filename as input and parses this file for + two-column floating point data, interpreting the data as + representing function values x and f(x). + + Ignores all lines not conforming to \\\\\\ + */ + MonotCubicInterpolator(const std::string & datafilename) throw (const char*) + { + if (!read(datafilename)) { + throw("Unable to constuct MonotCubicInterpolator from file.") ; + } ; + } ; + + + /** + @param datafilename A datafile with the x values and the corresponding f(x) values + + Accepts a filename as input and parses this file for + two-column floating point data, interpreting the data as + representing function values x and f(x). + + Ignores all lines not conforming to \\\\\\ + + All commas in the file will be treated as spaces when parsing. + + */ + + MonotCubicInterpolator(const char* datafilename) throw (const char*) + { + if (!read(std::string(datafilename))) { + throw("Unable to constuct MonotCubicInterpolator from file.") ; + } ; + } ; + + + /** + @param datafilename data file + @param XColumn x values + @param fColumn f values + + Accepts a filename as input, and parses the chosen columns in + that file. + */ + MonotCubicInterpolator(const char* datafilename, int xColumn, int fColumn) throw (const char*) + { + if (!read(std::string(datafilename),xColumn,fColumn)) { + throw("Unable to constuct MonotCubicInterpolator from file.") ; + } ; + } ; + + /** + @param datafilename data file + @param XColumn x values + @param fColumn f values + + Accepts a filename as input, and parses the chosen columns in + that file. + */ + MonotCubicInterpolator(const std::string & datafilename, int xColumn, int fColumn) throw (const char*) + { + if (!read(datafilename,xColumn,fColumn)) { + throw("Unable to constuct MonotCubicInterpolator from file.") ; + } ; + } ; + + /** + @param x vector of x values + @param f vector of corresponding f values + + Accepts two equal-length vectors as input for constructing + the interpolation object. First vector is the x-values, the + second vector is the function values + */ + MonotCubicInterpolator(const std::vector & x , + const std::vector & f); + + /** + No input, an empty function object is created. + + This object must be treated with care until + populated. + */ + MonotCubicInterpolator() {;} ; + + + + /** + @param datafilename A datafile with the x values and the corresponding f(x) values + + Accepts a filename as input and parses this file for + two-column floating point data, interpreting the data as + representing function values x and f(x). + + returns true on success + + All commas in file will be treated as spaces when parsing + + Ignores all lines not conforming to \\\\\\ + */ + bool read(const std::string & datafilename) { + return read(datafilename,1,2) ; + } ; + + /** + @param datafilename data file + @param XColumn x values + @param fColumn f values + + Accepts a filename as input, and parses the chosen columns in + that file. + */ + bool read(const std::string & datafilename, int xColumn, int fColumn) ; + + + + /** + @param x x value + + Returns f(x) for given x (input). Interpolates (monotone cubic + or linearly) if necessary. + + Extrapolates using the constants f(x_min) or f(x_max) if + input x is outside (x_min, x_max) + + @return f(x) for a given x + */ + double operator () (double x) const { return evaluate(x) ; } ; + + /** + @param x x value + + Returns f(x) for given x (input). Interpolates (monotone cubic + or linearly) if necessary. + + Extrapolates using the constants f(x_min) or f(x_max) if + input x is outside (x_min, x_max) + + @return f(x) for a given x + */ + double evaluate(double x) const throw(const char*); + + /** + @param x x value + @param errorestimate_output + + Returns f(x) and an error estimate for given x (input). + + Interpolates (linearly) if necessary. + + Throws an exception if extrapolation would be necessary for + evaluation. We do not want to do extrapolation (yet). + + The error estimate for x1 < x < x2 is + (x2 - x1)^2/8 * f''(x) where f''(x) is evaluated using + the stencil (1 -2 1) using either (x0, x1, x2) or (x1, x2, x3); + + Throws an exception if the table contains only two x-values. + + NOT IMPLEMENTED YET! + */ + double evaluate(double x, double & errorestimate_output ) const ; + + /** + Minimum x-value, returns both x and f in a pair. + + @return minimum x value + @return f(minimum x value) + */ + std::pair getMinimumX() const { + // Easy since the data is sorted on x: + return *data.begin(); + } + + /** + Maximum x-value, returns both x and f in a pair. + + @return maximum x value + @return f(maximum x value) + */ + std::pair getMaximumX() const { + // Easy since the data is sorted on x: + return *data.rbegin(); + } + + /** + Maximum f-value, returns both x and f in a pair. + + @return x value corresponding to maximum f value + @return maximum f value + */ + std::pair getMaximumF() const throw(const char*) ; + + /** + Minimum f-value, returns both x and f in a pair + + @return x value corresponding to minimal f value + @return minimum f value + */ + std::pair getMinimumF() const throw(const char*) ; + + + /** + Provide a copy of the x-data as a vector + + Unspecified order, but corresponds to get_fVector. + + @return x values as a vector + */ + std::vector get_xVector() const ; + + /** + Provide a copy of tghe function data as a vector + + Unspecified order, but corresponds to get_xVector + + @return f values as a vector + + */ + std::vector get_fVector() const ; + + /** + @param factor Scaling constant + + Scale all the function value data by a constant + */ + void scaleData(double factor); + + /** + Determines if the current function-value-data is strictly + monotone. This is a utility function for outsiders if they want + to invert the data for example. + + @return True if f(x) is strictly monotone, else False + */ + bool isStrictlyMonotone() { + + /* Use cached value if it can be trusted */ + if (strictlyMonotoneCached) { + return strictlyMonotone; + } + else { + computeInternalFunctionData(); + return strictlyMonotone; + } + } + + /** + Determines if the current function-value-data is monotone. + + @return True if f(x) is monotone, else False + */ + bool isMonotone() const { + if (monotoneCached) { + return monotone; + } + else { + computeInternalFunctionData(); + return monotone; + } + } + /** + Determines if the current function-value-data is strictly + increasing. This is a utility function for outsiders if they want + to invert the data for example. + + @return True if f(x) is strictly increasing, else False + */ + bool isStrictlyIncreasing() { + + /* Use cached value if it can be trusted */ + if (strictlyMonotoneCached) { + return (strictlyMonotone && strictlyIncreasing); + } + else { + computeInternalFunctionData(); + return (strictlyMonotone && strictlyIncreasing); + } + } + + /** + Determines if the current function-value-data is monotone and increasing. + + @return True if f(x) is monotone and increasing, else False + */ + bool isMonotoneIncreasing() const { + if (monotoneCached) { + return (monotone && increasing); + } + else { + computeInternalFunctionData(); + return (monotone && increasing); + } + } + /** + Determines if the current function-value-data is strictly + decreasing. This is a utility function for outsiders if they want + to invert the data for example. + + @return True if f(x) is strictly decreasing, else False + */ + bool isStrictlyDecreasing() { + + /* Use cached value if it can be trusted */ + if (strictlyMonotoneCached) { + return (strictlyMonotone && strictlyDecreasing); + } + else { + computeInternalFunctionData(); + return (strictlyMonotone && strictlyDecreasing); + } + } + + /** + Determines if the current function-value-data is monotone and decreasing + + @return True if f(x) is monotone and decreasing, else False + */ + bool isMonotoneDecreasing() const { + if (monotoneCached) { + return (monotone && decreasing); + } + else { + computeInternalFunctionData(); + return (monotone && decreasing); + } + } + + + + /** + @param newx New x point + @param newf New f(x) point + + Adds a new datapoint to the function. + + This causes all the derivatives at all points of the functions + to be recomputed and then adjusted for monotone cubic + interpolation. If this function ever enters a critical part of + any code, the locality of the algorithm for monotone adjustment + must be exploited. + + */ + void addPair(double newx, double newf) throw(const char*); + + /** + Returns an x-value that is believed to yield the best + improvement in global accuracy for the interpolation if + computed. + + Searches for the largest jump in f-values, and returns a x + value being the average of the two x-values representing the + f-value-jump. + + @return New x value beleived to yield the best improvement in global accuracy + @return Maximal difference + */ + std::pair getMissingX() const throw(const char*) ; + + /** + Constructs a string containing the data in a table + + @return a string containing the data in a table + */ + std::string toString() const; + + /** + @return Number of datapoint pairs in this object + */ + int getSize() const { + return data.size(); + } + + /** + Checks if the function curve is flat at the endpoints, chop off + endpoint data points if that is the case. + + The notion of "flat" is determined by the input parameter "epsilon" + Values whose difference are less than epsilon are regarded as equal. + + This is implemented to be able to obtain a strictly monotone + curve from a data set that is strictly monotone except at the + endpoints. + + Example: + The data points + (1,3), (2,3), (3,4), (4,5), (5,5), (6,5) + will become + (2,3), (3,4), (4,5) + + Assumes at least 3 datapoints. If less than three, this function is a noop. + */ + void chopFlatEndpoints(const double); + + /** + Wrapper function for chopFlatEndpoints(const double) + providing a default epsilon parameter + */ + void chopFlatEndpoints() { + chopFlatEndpoints(1e-14); + } + + /** + If function is monotone, but not strictly monotone, + this function will remove datapoints from intervals + with zero derivative so that the curve become + strictly monotone. + + Example + The data points + (1,2), (2,3), (3,4), (4,4), (5,5), (6,6) + will become + (1,2), (2,3), (3,4), (5,5), (6,6) + + Assumes at least two datapoints, if one or zero datapoint, this is a noop. + */ + void shrinkFlatAreas(const double); + + /** + Wrapper function for shrinkFlatAreas(const double) + providing a default epsilon parameter + */ + void shrinkFlatAreas() { + shrinkFlatAreas(1e-14); + }; + + + +private: + + // Data structure to store x- and f-values + std::map data; + + // Data structure to store x- and d-values + mutable std::map ddata; + + + // Storage containers for precomputed interpolation data + // std::vector dvalues; // derivatives in Hermite interpolation. + + // Flag to determine whether the boolean strictlyMonotone can be + // trusted. + mutable bool strictlyMonotoneCached; + mutable bool monotoneCached; /* only monotone, not stricly montone */ + + mutable bool strictlyMonotone; + mutable bool monotone; + + // if strictlyMonotone is true (and can be trusted), the two next are meaningful + mutable bool strictlyDecreasing; + mutable bool strictlyIncreasing; + mutable bool decreasing; + mutable bool increasing; + + + /* Hermite basis functions, t \in [0,1] , + notation from: + http://en.wikipedia.org/w/index.php?title=Cubic_Hermite_spline&oldid=84495502 + */ + + double H00(double t) const { + return 2*t*t*t - 3*t*t + 1; + } + double H10(double t) const { + return t*t*t - 2*t*t + t; + } + double H01(double t) const { + return -2*t*t*t + 3*t*t; + } + double H11(double t) const { + return t*t*t - t*t; + } + + + void computeInternalFunctionData() const ; + + /** + Computes initial derivative values using centered (second order) difference + for internal datapoints, and one-sided derivative for endpoints + + The internal datastructure map ddata is populated by this method. + */ + + void computeSimpleDerivatives() const ; + + + /** + Adjusts the derivative values (ddata) so that we can guarantee that + the resulting piecewise Hermite polymial is monotone. This is + done according to the algorithm of Fritsch and Carlsson 1980, + see Section 4, especially the two last lines. + */ + void adjustDerivativesForMonotoneness() const ; + + /** + Checks if the coefficient alpha and beta is in + the region that guarantees monotoneness of the + derivative values they represent + + See Fritsch and Carlson 1980, Lemma 2, + alternatively Step 5 in Wikipedia's article + on Monotone cubic interpolation. + */ + bool isMonotoneCoeff(double alpha, double beta) const { + if ((alpha*alpha + beta*beta) <= 9) { + return true; + } else { + return false; + } + } + +}; + + +#endif diff --git a/dune/common/RootFinders.hpp b/dune/common/RootFinders.hpp new file mode 100644 index 00000000..3011be18 --- /dev/null +++ b/dune/common/RootFinders.hpp @@ -0,0 +1,178 @@ +//=========================================================================== +// +// File: RootFinders.hpp +// +// Created: Thu May 6 19:59:42 2010 +// +// Author(s): Atgeirr F Rasmussen +// Jostein R Natvig +// +// $Date$ +// +// $Revision$ +// +//=========================================================================== + +/* + Copyright 2010 SINTEF ICT, Applied Mathematics. + Copyright 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 . +*/ + +#ifndef OPENRS_ROOTFINDERS_HEADER +#define OPENRS_ROOTFINDERS_HEADER + + + +//#include +//#include +//#include +//#include +#include + +namespace Dune +{ + + inline double regulaFalsiStep(const double a, + const double b, + const double fa, + const double fb) + { + ASSERT(fa*fb < 0.0); + return (b*fa - a*fb)/(fa - fb); + } + + + /// Implements a modified regula falsi method as described in + /// "Improved algorithms of Illinois-type for the numerical + /// solution of nonlinear equations" + /// by J. A. Ford. + /// Current variant is the 'Pegasus' method. + template + inline double modifiedRegulaFalsi(const Functor& f, + const double a, + const double b, + const int max_iter, + const double tolerance, + int& iterations_used) + { + using namespace std; + const double macheps = numeric_limits::epsilon(); + const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0); + + double x0 = a; + double x1 = b; + double f0 = f(x0); + const double epsF = tolerance + macheps*max(fabs(f0), 1.0); + if (fabs(f0) < epsF) { + return x0; + } + double f1 = f(x1); + if (fabs(f1) < epsF) { + return x1; + } + if (f0*f1 > 0.0) { + THROW("Error in parameters, zero not bracketed: [a, b] = [" + << a << ", " << b << "] fa = " << f0 << " fb = " << f1); + } + iterations_used = 0; + // In every iteraton, x1 is the last point computed, + // and x0 is the last point computed that makes it a bracket. + while (fabs(x1 - x0) >= 0.95*eps) { + double xnew = regulaFalsiStep(x0, x1, f0, f1); + double fnew = f(xnew); +// cout << "xnew = " << xnew << " fnew = " << fnew << endl; + ++iterations_used; + if (iterations_used > max_iter) { + THROW("Maximum number of iterations exceeded.\n" + << "Current interval is [" << min(x0, x1) << ", " + << max(x0, x1) << "]"); + } + if (fabs(fnew) < epsF) { + return xnew; + } + // Now we must check which point we must replace. + if ((fnew > 0.0) == (f0 > 0.0)) { + // We must replace x0. + x0 = x1; + f0 = f1; + } else { + // We must replace x1, this is the case where + // the modification to regula falsi kicks in, + // by modifying f0. + // 1. The classic Illinois method +// const double gamma = 0.5; + // @afr: The next two methods do not work??!!? + // 2. The method called 'Method 3' in the paper. +// const double phi0 = f1/f0; +// const double phi1 = fnew/f1; +// const double gamma = 1.0 - phi1/(1.0 - phi0); + // 3. The method called 'Method 4' in the paper. +// const double phi0 = f1/f0; +// const double phi1 = fnew/f1; +// const double gamma = 1.0 - phi0 - phi1; +// cout << "phi0 = " << phi0 <<" phi1 = " << phi1 << +// " gamma = " << gamma << endl; + // 4. The 'Pegasus' method + const double gamma = f1/(f1 + fnew); + f0 *= gamma; + } + x1 = xnew; + f1 = fnew; + } + return 0.5*(x0 + x1); + } + + + template + inline void bracketZero(const Functor& f, + const double x0, + const double dx, + double& a, + double& b) + { + const int max_iters = 100; + double f0 = f(x0); + double cur_dx = dx; + int i = 0; + for (; i < max_iters; ++i) { + double x = x0 + cur_dx; + double f_new = f(x); + if (f0*f_new <= 0.0) { + break; + } + cur_dx = -2.0*cur_dx; + } + if (i == max_iters) { + THROW("Could not bracket zero in " << max_iters << "iterations."); + } + if (cur_dx < 0.0) { + a = x0 + cur_dx; + b = i < 2 ? x0 : x0 + 0.25*cur_dx; + } else { + a = i < 2 ? x0 : x0 + 0.25*cur_dx; + b = x0 + cur_dx; + } + } + + +} // namespace Dune + + + + +#endif // OPENRS_ROOTFINDERS_HEADER diff --git a/dune/common/SparseTable.hpp b/dune/common/SparseTable.hpp new file mode 100644 index 00000000..86414923 --- /dev/null +++ b/dune/common/SparseTable.hpp @@ -0,0 +1,239 @@ +//=========================================================================== +// +// File: SparseTable.hpp +// +// Created: Fri Apr 24 09:50:27 2009 +// +// Author(s): Atgeirr F Rasmussen +// +// $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 . +*/ + +#ifndef OPENRS_SPARSETABLE_HEADER +#define OPENRS_SPARSETABLE_HEADER + +#include +#include +#include +#include +#include "ErrorMacros.hpp" + +#include + +namespace Dune +{ + + /// A SparseTable stores a table with rows of varying size + /// as efficiently as possible. + /// It is supposed to behave similarly to a vector of vectors. + /// Its behaviour is similar to compressed row sparse matrices. + template + class SparseTable + { + public: + /// Default constructor. Yields an empty SparseTable. + SparseTable() + { + } + + /// A constructor taking all the data for the table and row sizes. + /// \param data_beg The start of the table data. + /// \param data_end One-beyond-end of the table data. + /// \param rowsize_beg The start of the row length data. + /// \param rowsize_end One beyond the end of the row length data. + template + SparseTable(DataIter data_beg, DataIter data_end, + IntegerIter rowsize_beg, IntegerIter rowsize_end) + : data_(data_beg, data_end) + { + setRowStartsFromSizes(rowsize_beg, rowsize_end); + } + + + /// Sets the table to contain the given data, organized into + /// rows as indicated by the given row sizes. + /// \param data_beg The start of the table data. + /// \param data_end One-beyond-end of the table data. + /// \param rowsize_beg The start of the row length data. + /// \param rowsize_end One beyond the end of the row length data. + template + void assign(DataIter data_beg, DataIter data_end, + IntegerIter rowsize_beg, IntegerIter rowsize_end) + { + data_.assign(data_beg, data_end); + setRowStartsFromSizes(rowsize_beg, rowsize_end); + } + + + /// Request storage for table of given size. + /// \param rowsize_beg Start of row size data. + /// \param rowsize_end One beyond end of row size data. + template + void allocate(IntegerIter rowsize_beg, IntegerIter rowsize_end) + { + typedef typename std::vector::size_type sz_t; + + sz_t ndata = std::accumulate(rowsize_beg, rowsize_end, sz_t(0)); + data_.resize(ndata); + setRowStartsFromSizes(rowsize_beg, rowsize_end); + } + + + /// Appends a row to the table. + template + void appendRow(DataIter row_beg, DataIter row_end) + { + data_.insert(data_.end(), row_beg, row_end); + if (row_start_.empty()) { + row_start_.reserve(2); + row_start_.push_back(0); + } + row_start_.push_back(data_.size()); + } + + /// True if the table contains no rows. + bool empty() const + { + return row_start_.empty(); + } + + /// Returns the number of rows in the table. + int size() const + { + return empty() ? 0 : row_start_.size() - 1; + } + + /// Allocate storage for table of expected size + void reserve(int exptd_nrows, int exptd_ndata) + { + row_start_.reserve(exptd_nrows + 1); + data_.reserve(exptd_ndata); + } + + /// Swap contents for other SparseTable + void swap(SparseTable& other) + { + row_start_.swap(other.row_start_); + data_.swap(other.data_); + } + + /// Returns the number of data elements. + int dataSize() const + { + return data_.size(); + } + + /// Returns the size of a table row. + int rowSize(int row) const + { + ASSERT(row >= 0 && row < size()); + return row_start_[row + 1] - row_start_[row]; + } + + /// Makes the table empty(). + void clear() + { + data_.clear(); + row_start_.clear(); + } + + /// Defining the row type, returned by operator[]. + typedef boost::iterator_range row_type; + typedef boost::iterator_range mutable_row_type; + + /// Returns a row of the table. + row_type operator[](int row) const + { + ASSERT(row >= 0 && row < size()); + const T* start_ptr = data_.empty() ? 0 : &data_[0]; + return row_type(start_ptr + row_start_[row], start_ptr + row_start_[row + 1]); + } + + /// Returns a mutable row of the table. + mutable_row_type operator[](int row) + { + ASSERT(row >= 0 && row < size()); + T* start_ptr = data_.empty() ? 0 : &data_[0]; + return mutable_row_type(start_ptr + row_start_[row], start_ptr + row_start_[row + 1]); + } + + /// Equality. + bool operator==(const SparseTable& other) const + { + return data_ == other.data_ && row_start_ == other.row_start_; + } + + template + void print(std::basic_ostream& os) const + { + os << "Number of rows: " << size() << '\n'; + + os << "Row starts = ["; + std::copy(row_start_.begin(), row_start_.end(), + std::ostream_iterator(os, " ")); + os << "\b]\n"; + + os << "Data values = ["; + std::copy(data_.begin(), data_.end(), + std::ostream_iterator(os, " ")); + os << "\b]\n"; + } + + private: + std::vector data_; + // Like in the compressed row sparse matrix format, + // row_start_.size() is equal to the number of rows + 1. + std::vector row_start_; + + template + void setRowStartsFromSizes(IntegerIter rowsize_beg, IntegerIter rowsize_end) + { + // Since we do not store the row sizes, but cumulative row sizes, + // we have to create the cumulative ones. + int num_rows = rowsize_end - rowsize_beg; + if (num_rows < 1) { + THROW("Must have at least one row. Got " << num_rows << " rows."); + } +#ifndef NDEBUG + if (*std::min_element(rowsize_beg, rowsize_end) < 0) { + THROW("All row sizes must be at least 0."); + } +#endif + row_start_.resize(num_rows + 1); + row_start_[0] = 0; + std::partial_sum(rowsize_beg, rowsize_end, row_start_.begin() + 1); + // Check that data_ and row_start_ match. + if (int(data_.size()) != row_start_.back()) { + THROW("End of row start indices different from data size."); + } + + } + }; + +} // namespace Dune + + +#endif // OPENRS_SPARSETABLE_HEADER diff --git a/dune/common/SparseVector.hpp b/dune/common/SparseVector.hpp new file mode 100644 index 00000000..bce407d6 --- /dev/null +++ b/dune/common/SparseVector.hpp @@ -0,0 +1,195 @@ +//=========================================================================== +// +// File: SparseVector.hpp +// +// Created: Mon Jun 29 15:28:59 2009 +// +// Author(s): Atgeirr F Rasmussen +// Bård Skaflestad +// +// $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 . +*/ + + +#ifndef OPENRS_SPARSEVECTOR_HEADER +#define OPENRS_SPARSEVECTOR_HEADER + +#include +#include +#include +#include +#include "ErrorMacros.hpp" + +namespace Dune +{ + + /// A SparseVector stores a vector with possibly many empty elements + /// as efficiently as possible. + /// It is supposed to behave similarly to a standard vector, but since + /// direct indexing is a O(log n) operation instead of O(1), we do not + /// supply it as operator[]. + template + class SparseVector + { + public: + /// Default constructor. Yields an empty SparseVector. + SparseVector() + : size_(0), default_elem_() + { + } + + /// Constructs a SparseVector with a given size, but no nonzero + /// elements. + explicit SparseVector(int size) + : size_(size), default_elem_() + { + } + + /// A constructor taking all the element data for the vector and their indices. + /// \param data_beg The start of the element data. + /// \param data_end One-beyond-end of the element data. + /// \param rowsize_beg The start of the index data. + /// \param rowsize_end One beyond the end of the index data. + template + SparseVector(int size, + DataIter data_beg, DataIter data_end, + IntegerIter index_beg, IntegerIter index_end) + : size_(size), data_(data_beg, data_end), indices_(index_beg, index_end), + default_elem_() + { +#ifndef NDEBUG + ASSERT(size >= 0); + ASSERT(indices_.size() == data_.size()); + int last_index = -1; + int num_ind = indices_.size(); + for (int i = 0; i < num_ind; ++i) { + int index = indices_[i]; + if (index <= last_index || index >= size) { + THROW("Error in SparseVector construction, index is nonincreasing or out of range."); + } + last_index = index; + } +#endif + } + + + /// Appends an element to the vector. Note that this function does not + /// increase the size() of the vector, it just adds another nonzero element. + /// Elements must be added in index order. + void addElement(const T& element, int index) + { + ASSERT(indices_.empty() || index > indices_.back()); + ASSERT(index < size_); + data_.push_back(element); + indices_.push_back(index); + } + + /// \return true if the vector has size 0. + bool empty() const + { + return size_ == 0; + } + + /// Returns the size of the vector. + /// Recall that most or all of the vector may be default/zero. + int size() const + { + return size_; + } + + /// Returns the number of nonzero data elements. + int nonzeroSize() const + { + return data_.size(); + } + + /// Makes the vector empty(). + void clear() + { + data_.clear(); + indices_.clear(); + size_ = 0; + } + + /// Equality. + bool operator==(const SparseVector& other) const + { + return size_ == other.size_ && data_ == other.data_ && indices_ == other.indices_; + } + + /// O(log n) element access. + /// \param index the proper vector index + /// \return the element with the given index, or the default element if no element in + /// the vector has the given index. + const T& element(int index) const + { + ASSERT(index >= 0); + ASSERT(index < size_); + std::vector::const_iterator lb = std::lower_bound(indices_.begin(), indices_.end(), index); + if (lb != indices_.end() && *lb == index) { + return data_[lb - indices_.begin()]; + } else { + return default_elem_; + } + } + + /// O(1) element access. + /// \param nzindex an index counting only nonzero elements. + /// \return the nzindex'th nonzero element. + const T& nonzeroElement(int nzindex) const + { + ASSERT(nzindex >= 0); + ASSERT(nzindex < nonzeroSize()); + return data_[nzindex]; + } + + /// O(1) index access. + /// \param nzindex an index counting only nonzero elements. + /// \return the index of the nzindex'th nonzero element. + int nonzeroIndex(int nzindex) const + { + ASSERT(nzindex >= 0); + ASSERT(nzindex < nonzeroSize()); + return indices_[nzindex]; + } + + private: + // The vectors data_ and indices_ are always the same size. + // The indices are supposed to be stored in increasing order, + // to be unique, and to be in [0, size_ - 1]. + // default_elem_ is returned when a default element is requested. + int size_; + std::vector data_; + std::vector indices_; + T default_elem_; + }; + +} // namespace Dune + + + + +#endif // OPENRS_SPARSEVECTOR_HEADER diff --git a/dune/common/SpecialEclipseFields.hpp b/dune/common/SpecialEclipseFields.hpp new file mode 100644 index 00000000..81373ef3 --- /dev/null +++ b/dune/common/SpecialEclipseFields.hpp @@ -0,0 +1,1473 @@ +//=========================================================================== +// +// File: SpecialEclipseFields.hpp +// +// Created: Mon Sep 21 14:09:54 2009 +// +// Author(s): Atgeirr F Rasmussen +// BÃ¥rd Skaflestad +// BjÞrn Spjelkavik +// +// $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 . +*/ + +#ifndef OPENRS_SPECIALECLIPSEFIELDS_HEADER +#define OPENRS_SPECIALECLIPSEFIELDS_HEADER + +#include +#include +#include +#include +#include +#include + +namespace Dune +{ + +// Abstract base class for special fields. +struct SpecialBase { + virtual ~SpecialBase() {} // Default destructor + //virtual std::string name() const = 0; // Keyword name + virtual void read(std::istream& is) = 0; // Reads data + //virtual void write(std::ostream& os) const = 0; // Writes data + virtual void convertToSI(const EclipseUnits&) + { + THROW("Default conversion not defined."); + } + typedef std::vector > > table_t; +}; + + + + +/// Class for keyword SPECGRID +struct SPECGRID : public SpecialBase +{ + std::vector dimensions; // Number of grid blocks in x-, y- and z-directions. + int numres; // Number of reservoirs. + char qrdial; // Coordinates. F=cartesian, T=Cylindrical(radial). + + SPECGRID() + { + dimensions.resize(3,1); + numres = 1; + qrdial = 'F'; + } + + virtual ~SPECGRID() + { + } + + virtual std::string name() const {return std::string("SPECGRID");} + + virtual void read(std::istream& is) + { + const int ndim = 3; + std::vector data(ndim+1,1); + int nread = readDefaultedVectorData(is , data, ndim+1); + int nd = std::min(nread, ndim); + copy(data.begin(), data.begin()+nd, &dimensions[0]); + numres = data[ndim]; + std::string candidate; + is >> candidate; + if (candidate == "/") { + return; + } else { + qrdial = candidate[0]; + } + + if (ignoreSlashLine(is)) { + return; + } else { + THROW("End of file reading" << name()); + } + } + + virtual void write(std::ostream& os) const + { + os << name() << std::endl; + os << dimensions[0] << " " << dimensions[1] << " " + << dimensions[2] << " " << numres << " " << qrdial << std::endl; + os << std::endl; + } + + virtual void convertToSI(const EclipseUnits&) + {} +}; + + + + +/// Class holding segment data of keyword FAULTS. +struct FaultSegment +{ + std::string fault_name; // Fault name + std::vector ijk_coord; // ijk-coordinates of segment cells + std::string face; // Fault face of cells +}; + +/// Class for keyword FAULTS. +struct FAULTS : public SpecialBase +{ + std::vector faults; + + FAULTS() + { + } + + virtual ~FAULTS() + { + } + + virtual std::string name() const {return std::string("FAULTS");} + + virtual void read(std::istream& is) + { + while(is) { + std::string name; + is >> name; + if (name[0] == '/') { + is >> ignoreLine; + break; + } + while (name.find("--") == 0) { + // This line is a comment + is >> ignoreLine >> name; + } + FaultSegment fault_segment; + fault_segment.ijk_coord.resize(6); + fault_segment.fault_name = name; + int nread = readDefaultedVectorData(is, fault_segment.ijk_coord, 6); + if (nread != 6) { + THROW("Error reading fault_segment " << name); + } + is >> fault_segment.face; + faults.push_back(fault_segment); + ignoreSlashLine(is); + } + } + + virtual void write(std::ostream& os) const + { + os << name() << std::endl; + for (int i=0; i<(int)faults.size(); ++i) { + os << faults[i].fault_name << " "; + copy(faults[i].ijk_coord.begin(), faults[i].ijk_coord.end(), + std::ostream_iterator(os, " ")); + os << faults[i].face << std::endl; + } + os << std::endl; + } + + virtual void convertToSI(const EclipseUnits&) + {} +}; + + + + +/// Class holding a data line of keyword MULTFLT +struct MultfltLine +{ + std::string fault_name; // Fault name, as in FAULTS + double transmis_multiplier; // Transmissibility multiplier + double diffusivity_multiplier; // Diffusivity multiplier; + MultfltLine() : + fault_name(""), transmis_multiplier(1.0), diffusivity_multiplier(1.0) + {} +}; + +/// Class for keyword MULFLT +struct MULTFLT : public SpecialBase +{ + std::vector multflts; + + MULTFLT() + { + } + + virtual ~MULTFLT() + {} + + virtual std::string name() const {return std::string("MULTFLT");} + + virtual void read(std::istream& is) + { + while(is) { + std::string name; + is >> name; + if (name[0] == '/') { + is >> ignoreLine; + break; + } + while (name == "--") { + // This line is a comment + is >> ignoreLine >> name; + } + MultfltLine multflt_line; + multflt_line.fault_name = name; + std::vector data(2,1.0); + if (readDefaultedVectorData(is, data, 2) == 2) { + ignoreSlashLine(is); + } + multflt_line.transmis_multiplier = data[0]; + multflt_line.diffusivity_multiplier = data[1]; + multflts.push_back(multflt_line); + } + } + + virtual void write(std::ostream& os) const + { + os << name() << std::endl; + for (int i=0; i<(int)multflts.size(); ++i) { + os << multflts[i].fault_name << " " + << multflts[i].transmis_multiplier << " " + << multflts[i].diffusivity_multiplier << std::endl; + } + os << std::endl; + } + + virtual void convertToSI(const EclipseUnits&) + {} +}; + + + + +struct TITLE : public SpecialBase +{ + std::string title; + virtual std::string name() const + { return std::string("TITLE"); } + virtual void read(std::istream& is) + { is >> ignoreLine; std::getline(is, title); } + virtual void write(std::ostream& os) const + { os << name() << '\n' << title << '\n'; } + virtual void convertToSI(const EclipseUnits&) + {} +}; + + + + +struct START : public SpecialBase +{ + boost::gregorian::date date; + virtual std::string name() const + { return std::string("START"); } + virtual void read(std::istream& is) + { date = readDate(is); } + virtual void write(std::ostream& os) const + { os << name() << '\n' << date << '\n'; } + virtual void convertToSI(const EclipseUnits&) + {} +}; + + + + +struct DATES : public SpecialBase +{ + std::vector dates; + virtual std::string name() const + { return std::string("DATES"); } + virtual void read(std::istream& is) + { + while(is) { + dates.push_back(readDate(is)); + is >> ignoreWhitespace; + if (is.peek() == int('/')) { + is >> ignoreLine; + break; + } + } + } + virtual void write(std::ostream& os) const + { + os << name() << '\n'; + copy(dates.begin(), dates.end(), + std::ostream_iterator(os, "\n")); + } + virtual void convertToSI(const EclipseUnits&) + {} +}; + + +struct DENSITY : public SpecialBase +{ + std::vector > densities_; + + virtual std::string name() const {return std::string("DENSITY");} + + virtual void read(std::istream& is) + { + while (!is.eof()) { + std::vector density(3,-1e100); + if (readDefaultedVectorData(is, density, 3) == 3) { + ignoreSlashLine(is); + } + densities_.push_back(density); + + int action = next_action(is); // 0:continue 1:return 2:throw + if (action == 1) { + return; // Alphabetic char. Read next keyword. + } else if (action == 2) { + THROW("Error reading DENSITY. Next character is " + << (char)is.peek()); + } + } + } + + virtual void write(std::ostream& os) const + { + os << name() << '\n'; + for (int i=0; i<(int)densities_.size(); ++i) { + os << densities_[i][0] << " " << densities_[i][1] << " " + << densities_[i][2] << '\n'; + } + os << '\n'; + } + + virtual void convertToSI(const EclipseUnits& units) + { + for (int i=0; i<(int)densities_.size(); ++i) { + densities_[i][0] *= units.density; + densities_[i][1] *= units.density; + densities_[i][2] *= units.density; + } + } +}; + +struct PVDG : public SpecialBase +{ + table_t pvdg_; + + virtual std::string name() const {return std::string("PVDG");} + + virtual void read(std::istream& is) + { + readPvdTable(is, pvdg_, name(), 3); + } + + virtual void write(std::ostream& os) const + { + os << name() << '\n'; + for (int rn=0; rn<(int)pvdg_.size(); ++rn) { + for (int i=0; i<(int)pvdg_[rn][0].size(); ++i) { + os << pvdg_[rn][0][i] << " " << pvdg_[rn][1][i] << " " + << pvdg_[rn][2][i] << '\n'; + } + os << '\n'; + } + os << '\n'; + } + + virtual void convertToSI(const EclipseUnits& units) + { + double volfac = units.gasvol_r/units.gasvol_s; + for (int rn=0; rn<(int)pvdg_.size(); ++rn) { + for (int i=0; i<(int)pvdg_[rn][0].size(); ++i) { + pvdg_[rn][0][i] *= units.pressure; + pvdg_[rn][1][i] *= volfac; + pvdg_[rn][2][i] *= units.viscosity; + } + } + } +}; + +struct PVDO : public SpecialBase +{ + table_t pvdo_; + + virtual std::string name() const {return std::string("PVDO");} + + virtual void read(std::istream& is) + { + readPvdTable(is, pvdo_, name(), 3); + } + + virtual void write(std::ostream& os) const + { + os << name() << '\n'; + for (int rn=0; rn<(int)pvdo_.size(); ++rn) { + for (int i=0; i<(int)pvdo_[rn][0].size(); ++i) { + os << pvdo_[rn][0][i] << " " << pvdo_[rn][1][i] << " " + << pvdo_[rn][2][i] << '\n'; + } + os << '\n'; + } + os << '\n'; + } + + virtual void convertToSI(const EclipseUnits& units) + { + double volfac = units.liqvol_r/units.liqvol_s; + for (int rn=0; rn<(int)pvdo_.size(); ++rn) { + for (int i=0; i<(int)pvdo_[rn][0].size(); ++i) { + pvdo_[rn][0][i] *= units.pressure; + pvdo_[rn][1][i] *= volfac; + pvdo_[rn][2][i] *= units.viscosity; + } + } + } +}; + +struct PVTG : public SpecialBase +{ + table_t pvtg_; + + virtual std::string name() const {return std::string("PVTG");} + + virtual void read(std::istream& is) + { + readPvtTable(is, pvtg_, name()); + } + + virtual void write(std::ostream& os) const + { + os << name() << '\n'; + for (int rn=0; rn<(int)pvtg_.size(); ++rn) { + for (int i=0; i<(int)pvtg_[rn].size(); ++i) { + int nl = (pvtg_[rn][i].size()-1) / 3; + os << pvtg_[rn][i][0] << " " << pvtg_[rn][i][1] << " " + << pvtg_[rn][i][2] << " " << pvtg_[rn][i][3] << '\n'; + for (int j=1, n=3; j > pvtw_; + + virtual std::string name() const {return std::string("PVTW");} + + virtual void read(std::istream& is) + { + while (!is.eof()) { + std::vector pvtw; + readVectorData(is, pvtw); + if (pvtw.size() == 4) { + pvtw.push_back(0.0); // Not used by frontsim + } + pvtw_.push_back(pvtw); + + int action = next_action(is); // 0:continue 1:return 2:throw + if (action == 1) { + return; // Alphabetic char. Read next keyword. + } else if (action == 2) { + THROW("Error reading PVTW. Next character is " + << (char)is.peek()); + } + } + } + + virtual void write(std::ostream& os) const + { + os << name() << '\n'; + for (int i=0; i<(int)pvtw_.size(); ++i) { + os << pvtw_[i][0] << " " << pvtw_[i][1] << " " << pvtw_[i][2] + << " " << pvtw_[i][3] << " " << pvtw_[i][4] << '\n'; + } + os << '\n'; + } + + virtual void convertToSI(const EclipseUnits& units) + { + double volfac = units.liqvol_r/units.liqvol_s; + for (int i=0; i<(int)pvtw_.size(); ++i) { + pvtw_[i][0] *= units.pressure; + pvtw_[i][1] *= volfac; + pvtw_[i][2] *= units.compressibility; + pvtw_[i][3] *= units.viscosity; + pvtw_[i][4] *= units.compressibility; + } + } +}; + + +struct ROCK : public SpecialBase +{ + std::vector > rock_compressibilities_; + + virtual std::string name() const {return std::string("ROCK");} + + virtual void read(std::istream& is) + { + while (!is.eof()) { + std::vector rock; + readVectorData(is, rock); + rock_compressibilities_.push_back(rock); + + int action = next_action(is); // 0:continue 1:return 2:throw + if (action == 1) { + return; // Alphabetic char. Read next keyword. + } else if (action == 2) { + THROW("Error reading ROCK. Next character is " + << (char)is.peek()); + } + } + } + + virtual void write(std::ostream& os) const + { + os << name() << '\n'; + for (int i=0; i<(int)rock_compressibilities_.size(); ++i) { + os << rock_compressibilities_[i][0] << " " + << rock_compressibilities_[i][1] << '\n'; + } + os << '\n'; + } + + virtual void convertToSI(const EclipseUnits& units) + { + for (int i=0; i<(int)rock_compressibilities_.size(); ++i) { + rock_compressibilities_[i][0] *= units.pressure; + rock_compressibilities_[i][1] *= units.compressibility; + } + } +}; + + +struct ROCKTAB : public SpecialBase +{ + table_t rocktab_; + + virtual std::string name() const {return std::string("ROCKTAB");} + + virtual void read(std::istream& is) + { + readPvdTable(is, rocktab_, name(), 3); + } + + virtual void write(std::ostream& os) const + { + os << name() << '\n'; + for (int rn=0; rn<(int)rocktab_.size(); ++rn) { + for (int i=0; i<(int)rocktab_[rn][0].size(); ++i) { + os << rocktab_[rn][0][i] << " " << rocktab_[rn][1][i] << " " + << rocktab_[rn][2][i] << '\n'; + } + os << '\n'; + } + os << '\n'; + } + + virtual void convertToSI(const EclipseUnits& units) + { + for (int rn=0; rn<(int)rocktab_.size(); ++rn) { + for (int i=0; i<(int)rocktab_[rn][0].size(); ++i) { + rocktab_[rn][0][i] *= units.pressure; + } + } + } +}; + + +struct SGOF : public SpecialBase +{ + table_t sgof_; + + virtual std::string name() const {return std::string("SGOF");} + + virtual void read(std::istream& is) {readSGWOF(is, sgof_, name(), 4);} + + virtual void write(std::ostream& os) const + { + os << name() << '\n'; + for (int rn=0; rn<(int)sgof_.size(); ++rn) { + for (int i=0; i<(int)sgof_[rn][0].size(); ++i) { + os << sgof_[rn][0][i] << " " << sgof_[rn][1][i] << " " + << sgof_[rn][2][i] << " " << sgof_[rn][3][i] << '\n'; + } + os << '\n'; + } + os << '\n'; + } + + virtual void convertToSI(const EclipseUnits& units) + { + for (int rn=0; rn<(int)sgof_.size(); ++rn) { + for (int i=0; i<(int)sgof_[rn][0].size(); ++i) { + sgof_[rn][3][i] *= units.pressure; + } + } + } +}; + +struct SWOF : public SpecialBase +{ + table_t swof_; + + virtual std::string name() const {return std::string("SWOF");} + + virtual void read(std::istream& is) {readSGWOF(is, swof_, name(), 4);} + + virtual void write(std::ostream& os) const + { + os << name() << '\n'; + for (int rn=0; rn<(int)swof_.size(); ++rn) { + for (int i=0; i<(int)swof_[rn][0].size(); ++i) { + os << swof_[rn][0][i] << " " << swof_[rn][1][i] << " " + << swof_[rn][2][i] << " " << swof_[rn][3][i] << '\n'; + } + os << '\n'; + } + os << '\n'; + } + + virtual void convertToSI(const EclipseUnits& units) + { + for (int rn=0; rn<(int)swof_.size(); ++rn) { + for (int i=0; i<(int)swof_[rn][0].size(); ++i) { + swof_[rn][3][i] *= units.pressure; + } + } + } +}; + +/// Class holding a data line of keyword WELSPECS +struct WelspecsLine +{ + std::string name_; // Well name + std::string group_; // Group name + int I_; // I-location of well head or heel + int J_; // J-location of well head or heel + double datum_depth_BHP_; // Datum depth for bottom hole pressure + std::string pref_phase_; // Preferred phase for the well + double drain_rad_; // Drainage radius for prod/inj index calculation + std::string spec_inflow_; // Flag for special inflow equation + std::string shut_in_; // Instructions for automatic shut-in + std::string crossflow_; // Crossflow ability flag + int pressure_table_number_; // Pressure table number for wellbore fluid properties + std::string density_calc_type_; // Type of density calculation for wellbore hydrostatic head + int fluids_in_place_reg_numb_; // Fluids in place region number + + WelspecsLine() : + datum_depth_BHP_(-1.0), drain_rad_(0.0), spec_inflow_("STD"), + shut_in_("SHUT"), crossflow_("YES"), pressure_table_number_(0), + density_calc_type_("SEG"), fluids_in_place_reg_numb_(0) + {} +}; + +/// Class for keyword WELSPECS +struct WELSPECS : public SpecialBase +{ + std::vector welspecs; + + WELSPECS() + { + } + + virtual ~WELSPECS() + {} + + virtual std::string name() const {return std::string("WELSPECS");} + + virtual void read(std::istream& is) + { + while(is) { + std::string name = readString(is); + if (name[0] == '/') { + is >> ignoreLine; + break; + } + while (name.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + name = readString(is); + } + WelspecsLine welspecs_line; + welspecs_line.name_ = name; + welspecs_line.group_ = readString(is); + std::vector int_data(2,1); + readDefaultedVectorData(is, int_data, 2); + welspecs_line.I_ = int_data[0]; + welspecs_line.J_ = int_data[1]; + std::vector double_data(1,-1.0); + readDefaultedVectorData(is, double_data, 1); + welspecs_line.datum_depth_BHP_ = double_data[0]; + welspecs_line.pref_phase_ = readString(is); + + // HACK! Ignore items 7-13. + ignoreSlashLine(is); + welspecs.push_back(welspecs_line); + + // double_data[0] = 0.0; + // readDefaultedVectorData(is, double_data, 1); + // welspecs_line.drain_rad_ = double_data[0]; + // welspecs_line.spec_inflow_ = readString(is); + // welspecs_line.shut_in_ = readString(is); + // welspecs_line.crossflow_ = readString(is); + // int_data[0] = 0; + // readDefaultedVectorData(is, int_data, 1); + // welspecs_line.pressure_table_number_ = int_data[0]; + // welspecs_line.density_calc_type_ = readString(is); + // int_data[0] = 0; + // readDefaultedVectorData(is, int_data, 1); + // welspecs_line.fluids_in_place_reg_numb_ = int_data[0]; + // welspecs.push_back(welspecs_line); + } + } + + virtual void write(std::ostream& os) const + { + os << name() << std::endl; + for (int i=0; i<(int)welspecs.size(); ++i) { + os << welspecs[i].name_ << " " + << welspecs[i].group_ << " " + << welspecs[i].I_ << " " + << welspecs[i].J_ << " " + << welspecs[i].datum_depth_BHP_ << " " + << welspecs[i].pref_phase_ << " " + << welspecs[i].drain_rad_ << " " + << welspecs[i].shut_in_ << " " + << welspecs[i].crossflow_ << " " + << welspecs[i].pressure_table_number_ << " " + << welspecs[i].density_calc_type_ << " " + << welspecs[i].fluids_in_place_reg_numb_ << " " + << std::endl; + } + os << std::endl; + } + + virtual void convertToSI(const EclipseUnits& units) + { + for (int i=0; i<(int)welspecs.size(); ++i) { + welspecs[i].datum_depth_BHP_ *= units.length; + welspecs[i].drain_rad_ *= units.length; + } + } +}; + +/// Class holding a data line of keyword COMPDAT +struct CompdatLine +{ + std::string well_; // Well name + std::vector grid_ind_; // Grid block location + std::string open_shut_flag_; // Open/shut flag of connection + int sat_table_number_; // Saturation table number + double connect_trans_fac_; // Connection transmillibilty factor + double diameter_; // Well bore internal diameter + double Kh_; // Effective Kh value of the connection + double skin_factor_; // Skin factor + double D_factor_; // D-factor, for non-Darcy flow of free gas + std::string penetration_direct_; // Penetration direction + double r0_; // Pressure equivalent radius + + // Default values + CompdatLine() : + open_shut_flag_("OPEN"), + sat_table_number_(0), + connect_trans_fac_(0.0), + diameter_(0.0), + Kh_(-1.0), + skin_factor_(0.0), + D_factor_(-1e100), + penetration_direct_("Z"), + r0_(0.0) + { + grid_ind_.resize(4); + } +}; + +/// Class for keyword COMPDAT +struct COMPDAT : public SpecialBase +{ + std::vector compdat; + + COMPDAT() + { + } + + virtual ~COMPDAT() + {} + + virtual std::string name() const {return std::string("COMPDAT");} + + virtual void read(std::istream& is) + { + while(is) { + std::string name = readString(is); + if (name[0] == '/') { + is >> ignoreLine; + break; + } + while (name.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + name = readString(is); + } + CompdatLine compdat_line; + compdat_line.well_ = name; + readDefaultedVectorData(is, compdat_line.grid_ind_, 4); + compdat_line.open_shut_flag_ = readString(is); + std::vector int_data(1,-1); + readDefaultedVectorData(is, int_data, 1); + compdat_line.sat_table_number_ = int_data[0]; + std::vector double_data(2, 0.0); + int num_to_read = 2; + int num_read = readDefaultedVectorData(is, double_data, num_to_read); + compdat_line.connect_trans_fac_ = double_data[0]; + compdat_line.diameter_ = double_data[1]; + + // HACK! Ignore items 10-14. + if (num_read == num_to_read) { + ignoreSlashLine(is); + } + compdat.push_back(compdat_line); + } + } + + virtual void write(std::ostream& os) const + { + os << name() << std::endl; + for (int i=0; i<(int)compdat.size(); ++i) { + os << compdat[i].well_ << " " + << compdat[i].grid_ind_[0] << " " + << compdat[i].grid_ind_[1] << " " + << compdat[i].grid_ind_[2] << " " + << compdat[i].grid_ind_[3] << " " + << compdat[i].open_shut_flag_ << " " + << compdat[i].sat_table_number_ << " " + << compdat[i].connect_trans_fac_ << " " + << compdat[i].diameter_ << " " + << compdat[i].Kh_ << " " + << compdat[i].skin_factor_ << " " + << compdat[i].D_factor_ << " " + << compdat[i].penetration_direct_ << " " + << compdat[i].r0_ + << std::endl; + } + os << std::endl; + } + + virtual void convertToSI(const EclipseUnits& units) + { + for (int i=0; i<(int)compdat.size(); ++i) { + compdat[i].connect_trans_fac_ *= units.transmissibility; + compdat[i].diameter_ *= units.length; + compdat[i].Kh_ *= units.permeability*units.length; + compdat[i].r0_ *= units.length; + } + } +}; + +/// Class holding a data line of keyword WCONINJE +struct WconinjeLine +{ + std::string well_; // Well name or well name root + std::string injector_type_; // Injector type + std::string open_shut_flag_; // Open/shut flag for the well + std::string control_mode_; // Control mode + double surface_flow_max_rate_; // Surface flow rate target or upper limit + double fluid_volume_max_rate_; // Reservoir fluid volume rate target or + // upper limit + double BHP_limit_; // BHP target or upper limit + double THP_limit_; // THP target or upper limit + int VFP_table_number_; // Injection well VFP table number + double concentration_; // Vaporised oil concentration in the + // injected gas, or dissolved gas + // concentration in the injected oil + + // Default values + WconinjeLine() : + open_shut_flag_("OPEN"), surface_flow_max_rate_(1.0E20), + fluid_volume_max_rate_(1.0E20), BHP_limit_(6895), THP_limit_(1.0E20), + VFP_table_number_(0), concentration_(0.0) + { + } +}; + +/// Class for keyword WCONINJE +struct WCONINJE : public SpecialBase +{ + std::vector wconinje; + + WCONINJE() + { + } + + virtual ~WCONINJE() + {} + + virtual std::string name() const {return std::string("WCONINJE");} + + virtual void read(std::istream& is) + { + while(is) { + std::string name = readString(is); + if (name[0] == '/') { + is >> ignoreLine; + break; + } + while (name.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + name = readString(is); + } + WconinjeLine wconinje_line; + wconinje_line.well_ = name; + wconinje_line.injector_type_ = readString(is); + wconinje_line.open_shut_flag_ = readString(is); + wconinje_line.control_mode_ = readString(is); + std::vector double_data(6, 1.0E20); + double_data[2] = wconinje_line.BHP_limit_; + double_data[4] = wconinje_line.VFP_table_number_; + double_data[5] = wconinje_line.concentration_; + readDefaultedVectorData(is, double_data, 6); + wconinje_line.surface_flow_max_rate_ = double_data[0]; + wconinje_line.fluid_volume_max_rate_ = double_data[1]; + wconinje_line.BHP_limit_ = double_data[2]; + wconinje_line.THP_limit_ = double_data[3]; + wconinje_line.VFP_table_number_ = (int)double_data[4]; + wconinje_line.concentration_ = double_data[5]; + wconinje.push_back(wconinje_line); + } + } + + virtual void write(std::ostream& os) const + { + os << name() << std::endl; + for (int i=0; i<(int) wconinje.size(); ++i) { + os << wconinje[i].well_ << " " + << wconinje[i].injector_type_ << " " + << wconinje[i].open_shut_flag_ << " " + << wconinje[i].control_mode_ << " " + << wconinje[i].surface_flow_max_rate_ << " " + << wconinje[i].fluid_volume_max_rate_ << " " + << wconinje[i].BHP_limit_ << " " + << wconinje[i].THP_limit_ << " " + << wconinje[i].VFP_table_number_ << " " + << wconinje[i].concentration_ + << std::endl; + } + os << std::endl; + } + + virtual void convertToSI(const EclipseUnits& units) + { + for (int i=0; i<(int) wconinje.size(); ++i) { + if (wconinje[i].injector_type_ == "GAS") { + wconinje[i].surface_flow_max_rate_ *= units.gasvol_s/units.time; + } else { + wconinje[i].surface_flow_max_rate_ *= units.liqvol_s/units.time; + } + wconinje[i].fluid_volume_max_rate_ *= units.liqvol_r/units.time; + wconinje[i].BHP_limit_ *= units.pressure; + wconinje[i].THP_limit_ *= units.pressure; + wconinje[i].concentration_ *= units.gasvol_s/units.liqvol_s; // ??? @bsp 10 + } + } +}; + +/// Class holding a data line of keyword WCONPROD +struct WconprodLine +{ + std::string well_; // Well name or well name root + std::string open_shut_flag_; // Open/shut flag for the well + std::string control_mode_; // Control mode + double oil_max_rate_; // Oil rate target or upper limit + double water_max_rate_; // Water rate target or upper limit + double gas_max_rate_; // Gas rate target or upper limit + double liquid_max_rate_; // Liquid rate target or upper limit + double fluid_volume_max_rate_; // Reservoir fluid volume rate target or + // upper limit + double BHP_limit_; // BHP target or upper limit + double THP_limit_; // THP target or upper limit + int VFP_table_number_; // Injection well VFP table number + double artif_lift_quantity_; // Artificial lift quantity in THP calculations + + // Default values + WconprodLine() : + open_shut_flag_("OPEN"), oil_max_rate_(1.0E20), water_max_rate_(1.0E20), + gas_max_rate_(1.0E20), liquid_max_rate_(1.0E20), + fluid_volume_max_rate_(1.0E20), BHP_limit_(-1.0), THP_limit_(0.0), + VFP_table_number_(0), artif_lift_quantity_(0.0) + { + } +}; + +/// Class for keyword WCONPROD +struct WCONPROD : public SpecialBase +{ + std::vector wconprod; + + WCONPROD() + { + } + + virtual ~WCONPROD() + {} + + virtual std::string name() const {return std::string("WCONPROD");} + + virtual void read(std::istream& is) + { + while(is) { + std::string name = readString(is); + if (name[0] == '/') { + is >> ignoreLine; + break; + } + while (name.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + name = readString(is); + } + WconprodLine wconprod_line; + wconprod_line.well_ = name; + wconprod_line.open_shut_flag_ = readString(is); + wconprod_line.control_mode_ = readString(is); + std::vector double_data(9, 1.0E20); + double_data[5] = wconprod_line.BHP_limit_; + double_data[6] = wconprod_line.THP_limit_; + double_data[7] = wconprod_line.VFP_table_number_; + double_data[8] = wconprod_line.artif_lift_quantity_; + readDefaultedVectorData(is, double_data, 9); + wconprod_line.oil_max_rate_ = double_data[0]; + wconprod_line.water_max_rate_ = double_data[1]; + wconprod_line.gas_max_rate_ = double_data[2]; + wconprod_line.liquid_max_rate_ = double_data[3]; + wconprod_line.fluid_volume_max_rate_ = double_data[4]; + wconprod_line.BHP_limit_ = double_data[5]; + wconprod_line.THP_limit_ = double_data[6]; + wconprod_line.VFP_table_number_ = (int)double_data[7]; + wconprod_line.artif_lift_quantity_ = double_data[8]; + wconprod.push_back(wconprod_line); + } + } + + virtual void write(std::ostream& os) const + { + os << name() << std::endl; + for (int i=0; i<(int) wconprod.size(); ++i) { + os << wconprod[i].well_ << " " + << wconprod[i].open_shut_flag_ << " " + << wconprod[i].control_mode_ << " " + << wconprod[i].oil_max_rate_ << " " + << wconprod[i].water_max_rate_ << " " + << wconprod[i].gas_max_rate_ << " " + << wconprod[i].liquid_max_rate_ << " " + << wconprod[i].fluid_volume_max_rate_ << " " + << wconprod[i].BHP_limit_ << " " + << wconprod[i].THP_limit_ << " " + << wconprod[i].VFP_table_number_ << " " + << wconprod[i].artif_lift_quantity_ + << std::endl; + } + os << std::endl; + } + + virtual void convertToSI(const EclipseUnits& units) + { + double lrat = units.liqvol_s / units.time; + double grat = units.gasvol_s / units.time; + double resv = units.liqvol_r / units.time; + for (int i=0; i<(int) wconprod.size(); ++i) { + wconprod[i].oil_max_rate_ *= lrat; + wconprod[i].water_max_rate_ *= lrat; + wconprod[i].gas_max_rate_ *= grat; + wconprod[i].liquid_max_rate_ *= lrat; + wconprod[i].fluid_volume_max_rate_ *= resv; + wconprod[i].BHP_limit_ *= units.pressure; + wconprod[i].THP_limit_ *= units.pressure; + } + } +}; + + +/// Class holding a data line of keyword WELTARG +struct WeltargLine +{ + std::string well_; // Well name or well name root + std::string control_change_; // Definition of the control or constraint + // quantity to be changed + double new_value_; // New value of this quantity + + WeltargLine() + { + } +}; + +/// Class for keyword WELTARG +struct WELTARG : public SpecialBase +{ + std::vector weltarg; + + WELTARG() + { + } + + virtual ~WELTARG() + {} + + virtual std::string name() const {return std::string("WELTARG");} + + virtual void read(std::istream& is) + { + while(is) { + std::string name = readString(is); + if (name[0] == '/') { + is >> ignoreLine; + break; + } + while (name.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + name = readString(is); + } + WeltargLine weltarg_line; + weltarg_line.well_ = name; + weltarg_line.control_change_ = readString(is); + is >> weltarg_line.new_value_; + ignoreSlashLine(is); + weltarg.push_back(weltarg_line); + } + } + + virtual void write(std::ostream& os) const + { + os << name() << std::endl; + for (int i=0; i<(int)weltarg.size(); ++i) { + os << weltarg[i].well_ << " " + << weltarg[i].control_change_ << " " + << weltarg[i].new_value_ << " " + << std::endl; + } + os << std::endl; + } + + virtual void convertToSI(const EclipseUnits& units) + { + double lrat = units.liqvol_s / units.time; + double grat = units.gasvol_s / units.time; + double resv = units.liqvol_r / units.time; + for (int i=0; i<(int) weltarg.size(); ++i) { + if (weltarg[i].control_change_[0] == 'O' || + weltarg[i].control_change_[0] == 'W' || + weltarg[i].control_change_[0] == 'L') { + weltarg[i].new_value_ *= lrat; + } else if (weltarg[i].control_change_[0] == 'G') { + weltarg[i].new_value_ *= grat; + } else if (weltarg[i].control_change_[0] == 'R') { + weltarg[i].new_value_ *= resv; + } else if (weltarg[i].control_change_[0] == 'B' || + weltarg[i].control_change_[0] == 'T') { + weltarg[i].new_value_ *= units.pressure; + } else { + THROW("WELTARG. Unknown control or constraint " + << weltarg[i].control_change_[0]); + } + } + } +}; + + +/// Class holding a data line of keyword EQUIL +struct EquilLine +{ + double datum_depth_; // Well name or well name root. + double datum_depth_pressure_; // Pressure at datum depth. + double water_oil_contact_depth_; // Depth of water oil contact. + double oil_water_cap_pressure_; // Oil-water capillary pressure at the + // water-oil contact or Gas-water capillary + // pressure at the gas-water contact contact. + double gas_oil_contact_depth_; // Depth of the gas-oil contact. + double gas_oil_cap_pressure_; // Gas-oil capillary pressure at the gas-oil + // contact. + int live_oil_table_index_; // Rs v Depth or Pb v Depth table index for + // under-saturated live oil. + int wet_gas_table_index_; // Rv v Depth or Pd v Depth table index for + // under-saturated wet gas. + int N_; // Integer defining the accuracy of the + // initial fluids in place calculation. + EquilLine() + { + } +}; + +/// Class for keyword EQUIL +struct EQUIL : public SpecialBase +{ + std::vector equil; + + EQUIL() + { + } + + virtual ~EQUIL() + {} + + virtual std::string name() const {return std::string("EQUIL");} + + virtual void read(std::istream& is) + { + // Note. This function assumes that NTEQUL = 1, and reads only one line. + int num_to_read = 9; + std::vector data(num_to_read,0); + int num_read = readDefaultedVectorData(is, data, num_to_read); + if (num_read == num_to_read) { + ignoreSlashLine(is); + } + + EquilLine equil_line; + equil_line.datum_depth_ = data[0]; + equil_line.datum_depth_pressure_ = data[1]; + equil_line.water_oil_contact_depth_ = data[2]; + equil_line.oil_water_cap_pressure_ = data[3]; + equil_line.gas_oil_contact_depth_ = data[4]; + equil_line.gas_oil_cap_pressure_ = data[5]; + equil_line.live_oil_table_index_ = int(data[6]); + equil_line.wet_gas_table_index_ = int(data[7]); + equil_line.N_ = int(data[8]); + equil.push_back(equil_line); + } + + virtual void write(std::ostream& os) const + { + os << name() << std::endl; + for (int i=0; i<(int)equil.size(); ++i) { + os << equil[i].datum_depth_ << " " + << equil[i].datum_depth_pressure_ << " " + << equil[i].water_oil_contact_depth_ << " " + << equil[i].oil_water_cap_pressure_ << " " + << equil[i].gas_oil_contact_depth_ << " " + << equil[i].gas_oil_cap_pressure_ << " " + << equil[i].live_oil_table_index_ << " " + << equil[i].wet_gas_table_index_ << " " + << equil[i].N_ << " " + + << std::endl; + } + os << std::endl; + } + + virtual void convertToSI(const EclipseUnits& units) + { + for (int i=0; i<(int)equil.size(); ++i) { + equil[i].datum_depth_ *= units.length; + equil[i].datum_depth_pressure_ *= units.pressure; + equil[i].water_oil_contact_depth_ *= units.length; + equil[i].oil_water_cap_pressure_ *= units.pressure; + equil[i].gas_oil_contact_depth_ *= units.length; + equil[i].gas_oil_cap_pressure_ *= units.pressure; + } + } +}; + +struct PVCDO : public SpecialBase +{ + std::vector > pvcdo_; + + virtual std::string name() const {return std::string("PVCDO");} + + virtual void read(std::istream& is) + { + while (!is.eof()) { + std::vector pvcdo; + readVectorData(is, pvcdo); + if (pvcdo.size() == 4) { + pvcdo.push_back(0.0); + } + pvcdo_.push_back(pvcdo); + + int action = next_action(is); // 0:continue 1:return 2:throw + if (action == 1) { + return; // Alphabetic char. Read next keyword. + } else if (action == 2) { + THROW("Error reading PVCDO. Next character is " + << (char)is.peek()); + } + } + } + + virtual void write(std::ostream& os) const + { + os << name() << '\n'; + for (int i=0; i<(int)pvcdo_.size(); ++i) { + os << pvcdo_[i][0] << " " << pvcdo_[i][1] << " " << pvcdo_[i][2] + << " " << pvcdo_[i][3] << " " << pvcdo_[i][4] << '\n'; + } + os << '\n'; + } + + virtual void convertToSI(const EclipseUnits& units) + { + double volfac = units.liqvol_r/units.liqvol_s; + for (int i=0; i<(int)pvcdo_.size(); ++i) { + pvcdo_[i][0] *= units.pressure; + pvcdo_[i][1] *= volfac; + pvcdo_[i][2] *= units.compressibility; + pvcdo_[i][3] *= units.viscosity; + pvcdo_[i][4] *= units.compressibility; + } + } +}; + +struct MultRec : public SpecialBase +{ + virtual void read(std::istream& is) + { +#ifdef VERBOSE + std::cout << "(dummy implementation)" << std::endl; +#endif + const std::ctype& ct = std::use_facet< std::ctype >(std::locale::classic()); + is >> ignoreSlashLine; + while (!is.eof()) { + is >> ignoreWhitespace; + std::streampos pos = is.tellg(); + char c; + is.get(c); + if (is.eof()) { + return; + } + if (ct.is(std::ctype_base::alpha, c)) { + std::string name; // Unquoted name or new keyword? + std::getline(is, name); + if (name.rfind('/') != std::string::npos) { + continue; // Unquoted name + } else { + is.seekg(pos); + break; // Read next keyword + } + } else if (ct.is(std::ctype_base::digit, c) || c== '.') { + is >> ignoreSlashLine; // Decimal digit. Ignore data. + continue; + } else if (c== '\'') { + is >> ignoreSlashLine; // Quote. Ignore data. + continue; + } else if(c == '-' && is.peek() == int('-')) { + is >> ignoreLine; // This line is a comment + continue; + } else if (c == '/' ) { + is >> ignoreLine; // This line is a null record. + continue; // (No data before slash) + } else { + is.putback(c); + std::string temp; + is >> temp; + std::cout << "READ ERROR! Next word is " << temp << std::endl; + } + } + } + + virtual void convertToSI(const EclipseUnits&) + {} + +}; + +// The following fields only have a dummy implementation +// that allows us to ignore them. +struct SWFN : public MultRec {}; +struct SOF2 : public MultRec {}; +struct TUNING : public MultRec {}; + + +} // End of namespace Dune + +#endif // OPENRS_SPECIALECLIPSEFIELDS_HEADER + diff --git a/dune/common/StopWatch.cpp b/dune/common/StopWatch.cpp new file mode 100644 index 00000000..1d5a2bbd --- /dev/null +++ b/dune/common/StopWatch.cpp @@ -0,0 +1,106 @@ +//=========================================================================== +// +// File: StopWatch.cpp +// +// Created: Thu Jul 2 23:04:51 2009 +// +// Author(s): Atgeirr F Rasmussen +// +// $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 . +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif +#include +#include "StopWatch.hpp" +#include "ErrorMacros.hpp" + +namespace Dune +{ + + namespace time + { + + StopWatch::StopWatch() + : state_(UnStarted) + { + } + + + void StopWatch::start() + { + start_time_ = boost::posix_time::microsec_clock::local_time(); + last_time_ = start_time_; + state_ = Running; + } + + void StopWatch::stop() + { + if (state_ != Running) { + THROW("Called stop() on a StopWatch that was not running."); + } + stop_time_ = boost::posix_time::microsec_clock::local_time(); + state_ = Stopped; + } + + double StopWatch::secsSinceLast() + { + boost::posix_time::ptime run_time; + if (state_ == Running) { + run_time = boost::posix_time::microsec_clock::local_time(); + } else if (state_ == Stopped) { + run_time = stop_time_; + } else { + ASSERT(state_ == UnStarted); + THROW("Called secsSinceLast() on a StopWatch that had not been started."); + } + boost::posix_time::time_duration dur = run_time - last_time_; + last_time_ = run_time; + return double(dur.total_microseconds())/1000000.0; + } + + double StopWatch::secsSinceStart() + { + boost::posix_time::ptime run_time; + if (state_ == Running) { + run_time = boost::posix_time::microsec_clock::local_time(); + } else if (state_ == Stopped) { + run_time = stop_time_; + } else { + ASSERT(state_ == UnStarted); + THROW("Called secsSinceStart() on a StopWatch that had not been started."); + } + boost::posix_time::time_duration dur = run_time - start_time_; + return double(dur.total_microseconds())/1000000.0; + } + + } // namespace time + +} // namespace Dune + + + diff --git a/dune/common/StopWatch.hpp b/dune/common/StopWatch.hpp new file mode 100644 index 00000000..50106413 --- /dev/null +++ b/dune/common/StopWatch.hpp @@ -0,0 +1,81 @@ +//=========================================================================== +// +// File: StopWatch.hpp +// +// Created: Thu Jul 2 23:04:17 2009 +// +// Author(s): Atgeirr F Rasmussen +// +// $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 . +*/ + +#ifndef OPENRS_STOPWATCH_HEADER +#define OPENRS_STOPWATCH_HEADER + +#include + +namespace Dune +{ + + namespace time + { + + class StopWatch + { + public: + /// Default constructor. Before the StopWatch is start()-ed, + /// it is an error to call anything other than start(). + StopWatch(); + + /// Starts the StopWatch. It is always legal to call + /// start(), even if not stop()-ped. + void start(); + /// Stops the StopWatch. The watch no longer runs, until + /// restarted by a call to start(). + void stop(); + + /// \ret the number of running seconds that have passed + /// since last call to start(), secsSinceLast() or + /// secsSinceStart() + double secsSinceLast(); + /// \ret the number of running seconds that have passed + /// since last call to start(). + double secsSinceStart(); + + private: + enum StopWatchState { UnStarted, Running, Stopped }; + + StopWatchState state_; + boost::posix_time::ptime start_time_; + boost::posix_time::ptime last_time_; + boost::posix_time::ptime stop_time_; + }; + + } // namespace time + +} // namespace Dune + +#endif // OPENRS_STOPWATCH_HEADER diff --git a/dune/common/Units.hpp b/dune/common/Units.hpp new file mode 100644 index 00000000..34e46d8c --- /dev/null +++ b/dune/common/Units.hpp @@ -0,0 +1,182 @@ +//=========================================================================== +// +// File: Units.hpp +// +// Created: Thu Jul 2 09:19:08 2009 +// +// Author(s): Halvor M Nilsen +// +// $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 . +*/ + +#ifndef OPENRS_UNITS_HEADER +#define OPENRS_UNITS_HEADER + +namespace Dune +{ + namespace prefix + { + const double micro = 1.0e-6; + const double milli = 1.0e-3; + const double centi = 1.0e-2; + const double deci = 1.0e-1; + const double kilo = 1.0e3; + const double mega = 1.0e6; + const double giga = 1.0e9; + } // namespace prefix + + namespace unit + { + // Common powers + inline double square(double v) { return v * v; } + inline double cubic (double v) { return v * v * v; } + + // -------------------------------------------------------------- + // Basic (fundamental) units and conversions + // -------------------------------------------------------------- + + // Length: + const double meter = 1; + const double inch = 2.54 * prefix::centi*meter; + const double feet = 12 * inch; + + // Time: + const double second = 1; + const double minute = 60 * second; + const double hour = 60 * minute; + const double day = 24 * hour; + const double year = 365 * day; + + // Volume + const double stb = 0.158987294928 * cubic(meter); + + + // Mass: + const double kilogram = 1; + + // http://en.wikipedia.org/wiki/Pound_(mass)#Avoirdupois_pound + const double pound = 0.45359237 * kilogram; + + // -------------------------------------------------------------- + // Standardised constants + // -------------------------------------------------------------- + + const double gravity = 9.80665 * meter/square(second); + + // -------------------------------------------------------------- + // Derived units and conversions + // -------------------------------------------------------------- + + // Force: + const double Newton = kilogram*meter / square(second); // == 1 + const double lbf = pound * gravity; // Pound-force + + // Pressure: + const double Pascal = Newton / square(meter); // == 1 + const double barsa = 100000 * Pascal; + const double atm = 101325 * Pascal; + const double psia = lbf / square(inch); + + // Viscosity: + const double Pas = Pascal * second; // == 1 + const double Poise = prefix::deci*Pas; + + // Permeability: + // + // A porous medium with a permeability of 1 darcy permits a + // flow (flux) of 1 cm³/s of a fluid with viscosity 1 cP (1 + // mPa·s) under a pressure gradient of 1 atm/cm acting across + // an area of 1 cm². + // + namespace perm_details { + const double p_grad = atm / (prefix::centi*meter); + const double area = square(prefix::centi*meter); + const double flux = cubic (prefix::centi*meter) / second; + const double velocity = flux / area; + const double visc = prefix::centi*Poise; + const double darcy = (velocity * visc) / p_grad; + // == 1e-7 [m^2] / 101325 + // == 9.869232667160130e-13 [m^2] + } + const double darcy = perm_details::darcy; + + // Unit conversion support. + // + // Note: Under the penalty of treason will you be + // + // using namespace Dune::unit::convert; + // + // I mean it! + // + namespace convert { + // Convert from external units of measurements to equivalent + // internal units of measurements. Note: The internal units + // of measurements are *ALWAYS*, and exclusively, SI. + // + // Example: Convert a double kx, containing a permeability + // value in units of milli-darcy (mD) to the equivalent + // value in SI units (m^2). + // + // using namespace Dune::unit; + // using namespace Dune::prefix; + // convert::from(kx, milli*darcy); + // + inline double from(const double q, const double unit) + { + return q * unit; + } + + // Convert from internal units of measurements to equivalent + // external units of measurements. Note: The internal units + // of measurements are *ALWAYS*, and exclusively, SI. + // + // Example: Convert a std::vector p, containing + // pressure values in the SI unit Pascal (i.e., unit::Pascal) + // to the equivalent values in Psi (unit::psia). + // + // using namespace Dune::unit; + // std::transform(p.begin(), p.end(), p.begin(), + // boost::bind(convert::to, _1, psia)); + // + inline double to(const double q, const double unit) + { + return q / unit; + } + } // namespace convert + } // namespace unit + + namespace units { + // const double MILLIDARCY = 1.0;//9.86923e-16; + // const double VISCOSITY_UNIT = 1.0;//1e-3; + // const double DAYS2SECONDS = 1.0;//86400; + const double MILLIDARCY = 9.86923e-16; + const double VISCOSITY_UNIT = 1e-3; + const double DAYS2SECONDS = 86400; + const double FEET = 0.30479999798832; + const double WELL_INDEX_UNIT = VISCOSITY_UNIT/(DAYS2SECONDS*1e5); + } // namespace units +} // namespace Dune +#endif // OPENRS_UNITS_HEADER diff --git a/dune/common/linInt.hpp b/dune/common/linInt.hpp new file mode 100644 index 00000000..f22a25bf --- /dev/null +++ b/dune/common/linInt.hpp @@ -0,0 +1,115 @@ +//=========================================================================== +// +// File: linInt.hpp +// +// Created: Tue Feb 16 14:44:10 2010 +// +// Author(s): Bjørn Spjelkavik +// Atgeirr F Rasmussen +// Bård Skaflestad +// +// $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 . +*/ + +#ifndef OPENRS_LININT_HEADER +#define OPENRS_LININT_HEADER + + +namespace Dune +{ + + inline int tableIndex(const std::vector& table, double x) + { + // Returns an index in an ordered table such that x is between + // table[j] and table[j+1]. If x is out of range, first or last + // interval is returned; Binary search. + int n = table.size() - 1; + if (n < 2) { + return 0; + } + int jl = 0; + int ju = n; + bool ascend = (table[n] > table[0]); + while (ju - jl > 1) { + int jm = (ju + jl)/2; // Compute a midpoint + if ( (x >= table[jm]) == ascend) { + jl = jm; // Replace lower limit + } else { + ju = jm; // Replace upper limit + } + } + return jl; + } + + inline double linearInterpolation(const std::vector& xv, + const std::vector& yv, double x) + { + // Returns end point if x is outside xv + std::vector::const_iterator lb = lower_bound(xv.begin(), xv.end(), x); + int ix2 = lb - xv.begin(); + if (ix2 == 0) { + return yv[0]; + } else if (ix2 == int(xv.size())) { + return yv[ix2-1]; + } + int ix1 = ix2 - 1; + return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1]; + } + + inline double linearInterpolDerivative(const std::vector& xv, + const std::vector& yv, double x) + { + int ix1 = tableIndex(xv, x); + int ix2 = ix1 + 1; + return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1]); + } + + inline double linearInterpolationExtrap(const std::vector& xv, + const std::vector& yv, double x) + { + // Extrapolates if x is outside xv + int ix1 = tableIndex(xv, x); + int ix2 = ix1 + 1; + return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1]; + } + + inline double linearInterpolationExtrap(const std::vector& xv, + const std::vector& yv, + double x, int& ix1) + { + // Extrapolates if x is outside xv + ix1 = tableIndex(xv, x); + int ix2 = ix1 + 1; + return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1]; + } + +} // namespace Dune + + + + + +#endif // OPENRS_LININT_HEADER diff --git a/dune/common/param/Makefile.am b/dune/common/param/Makefile.am new file mode 100644 index 00000000..91cbbe3c --- /dev/null +++ b/dune/common/param/Makefile.am @@ -0,0 +1,21 @@ +# $Date$ +# $Revision$ + +SUBDIRS = . test + +paramdir = $(includedir)/dune/common/param + +param_HEADERS = Parameter.hpp ParameterGroup.hpp ParameterGroup_impl.hpp \ + ParameterMapItem.hpp ParameterRequirement.hpp \ + ParameterStrings.hpp ParameterTools.hpp ParameterXML.hpp + +noinst_LTLIBRARIES = libparam.la + +libparam_la_SOURCES = Parameter.cpp ParameterGroup.cpp ParameterTools.cpp \ + ParameterXML.cpp + +libparam_la_CPPFLAGS = $(XML_CPPFLAGS) $(BOOST_CPPFLAGS) +libparam_la_LDFLAGS = $(BOOST_LDFLAGS) +libparam_la_LIBADD = $(XML_LIBS) $(BOOST_FILESYSTEM_LIB) $(BOOST_SYSTEM_LIB) + +include $(top_srcdir)/am/global-rules diff --git a/dune/common/param/Parameter.cpp b/dune/common/param/Parameter.cpp new file mode 100644 index 00000000..dfb27b9d --- /dev/null +++ b/dune/common/param/Parameter.cpp @@ -0,0 +1,74 @@ +//=========================================================================== +// +// File: Parameter.cpp +// +// Created: Tue Jun 2 19:18:25 2009 +// +// Author(s): Bård Skaflestad +// Atgeirr F Rasmussen +// +// $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 . +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif +#include +#include + +namespace Dune { + namespace parameter { + std::string + correct_parameter_tag(const ParameterMapItem& item) + { + std::string tag = item.getTag(); + if (tag != ID_xmltag__param) { + std::string error = "The XML tag was '" + + tag + "' but should be '" + + ID_xmltag__param + "'.\n"; + return error; + } else { + return ""; + } + } + + std::string + correct_type(const Parameter& parameter, + const std::string& param_type) + { + std::string type = parameter.getType(); + if ( (type != param_type) && + (type != ID_param_type__cmdline) ) { + std::string error = "The data was of type '" + type + + "' but should be of type '" + + param_type + "'.\n"; + return error; + } else { + return ""; + } + } + } // namespace parameter +} // namespace Dune diff --git a/dune/common/param/Parameter.hpp b/dune/common/param/Parameter.hpp new file mode 100644 index 00000000..078e7416 --- /dev/null +++ b/dune/common/param/Parameter.hpp @@ -0,0 +1,214 @@ +//=========================================================================== +// +// File: Parameter.hpp +// +// Created: Tue Jun 2 16:00:21 2009 +// +// Author(s): Bård Skaflestad +// Atgeirr F Rasmussen +// +// $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 . +*/ + +#ifndef OPENRS_PARAMETER_HEADER +#define OPENRS_PARAMETER_HEADER + +#include +#include + +#include +#include + +namespace Dune { + /// See ParameterGroup.hpp for how to use the parameter system + namespace parameter { + /// @brief + /// @todo Doc me! + class Parameter : public ParameterMapItem { + public: + /// @brief + /// @todo Doc me! + virtual ~Parameter() {} + /// @brief + /// @todo Doc me! + /// @return + virtual std::string getTag() const {return ID_xmltag__param;} + /// @brief + /// @todo Doc me! + /// @param + Parameter(const std::string& value, const std::string& type) + : value_(value), type_(type) {} + /// @brief + /// @todo Doc me! + /// @return + std::string getValue() const {return value_;} + /// @brief + /// @todo Doc me! + /// @return + std::string getType() const {return type_;} + private: + std::string value_; + std::string type_; + }; + + /// @brief + /// @todo Doc me! + /// @param + /// @return + std::string correct_parameter_tag(const ParameterMapItem& item); + std::string correct_type(const Parameter& parameter, + const std::string& type); + + /// @brief + /// @todo Doc me! + /// @tparam + /// @param + /// @return + template<> + struct ParameterMapItemTrait { + static int convert(const ParameterMapItem& item, + std::string& conversion_error) + { + conversion_error = correct_parameter_tag(item); + if (conversion_error != "") { + return 0; + } + const Parameter& parameter = dynamic_cast(item); + conversion_error = correct_type(parameter, ID_param_type__int); + if (conversion_error != "") { + return 0; + } + std::stringstream stream; + stream << parameter.getValue(); + int value; + stream >> value; + if (stream.fail()) { + conversion_error = "Conversion to '" + + ID_param_type__int + + "' failed. Data was '" + + parameter.getValue() + "'.\n"; + return 0; + } + return value; + } + static std::string type() {return ID_param_type__int;} + }; + + /// @brief + /// @todo Doc me! + /// @tparam + /// @param + /// @return + template<> + struct ParameterMapItemTrait { + static double convert(const ParameterMapItem& item, + std::string& conversion_error) + { + conversion_error = correct_parameter_tag(item); + if (conversion_error != "") { + return 0.0; + } + const Parameter& parameter = dynamic_cast(item); + conversion_error = correct_type(parameter, ID_param_type__float); + if (conversion_error != "") { + return 0.0; + } + std::stringstream stream; + stream << parameter.getValue(); + double value; + stream >> value; + if (stream.fail()) { + conversion_error = "Conversion to '" + + ID_param_type__float + + "' failed. Data was '" + + parameter.getValue() + "'.\n"; + return 0.0; + } + return value; + } + static std::string type() {return ID_param_type__float;} + }; + + /// @brief + /// @todo Doc me! + /// @tparam + /// @param + /// @return + template<> + struct ParameterMapItemTrait { + static bool convert(const ParameterMapItem& item, + std::string& conversion_error) + { + conversion_error = correct_parameter_tag(item); + if (conversion_error != "") { + return false; + } + const Parameter& parameter = dynamic_cast(item); + conversion_error = correct_type(parameter, ID_param_type__bool); + if (conversion_error != "") { + return false; + } + if (parameter.getValue() == ID_true) { + return true; + } else if (parameter.getValue() == ID_false) { + return false; + } else { + conversion_error = "Conversion failed. Data was '" + + parameter.getValue() + + "', but should be one of '" + + ID_true + "' or '" + ID_false + "'.\n"; + return false; + } + } + static std::string type() {return ID_param_type__bool;} + }; + + /// @brief + /// @todo Doc me! + /// @tparam + /// @param + /// @return + template<> + struct ParameterMapItemTrait { + static std::string convert(const ParameterMapItem& item, + std::string& conversion_error) + { + conversion_error = correct_parameter_tag(item); + if (conversion_error != "") { + return ""; + } + const Parameter& parameter = dynamic_cast(item); + conversion_error = correct_type(parameter, ID_param_type__string); + if (conversion_error != "") { + return ""; + } + return parameter.getValue(); + } + static std::string type() {return ID_param_type__string;} + }; + } // namespace parameter +} // namespace Dune +#endif // OPENRS_PARAMETER_HPP diff --git a/dune/common/param/ParameterGroup.cpp b/dune/common/param/ParameterGroup.cpp new file mode 100644 index 00000000..6cb729f3 --- /dev/null +++ b/dune/common/param/ParameterGroup.cpp @@ -0,0 +1,334 @@ +//=========================================================================== +// +// File: ParameterGroup.cpp +// +// Created: Tue Jun 2 19:13:17 2009 +// +// Author(s): Bård Skaflestad +// Atgeirr F Rasmussen +// +// $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 . +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif +#include + +#include +#include +#include +#include + +#include +#include + +//#include "TermColors.hpp" // from utils + +#include +#include +#include +#include + +namespace Dune { + namespace parameter { + + ParameterGroup::ParameterGroup() + : path_(ID_path_root), parent_(0), output_is_enabled_(true) + { + } + + ParameterGroup::~ParameterGroup() { + if (output_is_enabled_) { + //TermColors::Normal(); + } + } + + std::string ParameterGroup::getTag() const { + return ID_xmltag__param_grp; + } + + + bool ParameterGroup::has(const std::string& name) const + { + std::pair name_path = split(name); + map_type::const_iterator it = map_.find(name_path.first); + if (it == map_.end()) { + return false; + } else if (name_path.second == "") { + return true; + } else { + ParameterGroup& pg = dynamic_cast(*(*it).second); + return pg.has(name_path.second); + } + } + + std::string ParameterGroup::path() const { + return path_; + } + + ParameterGroup::ParameterGroup(const std::string& path, + const ParameterGroup* parent) + : path_(path), parent_(parent), output_is_enabled_(true) + { + } + + ParameterGroup + ParameterGroup::getGroup(const std::string& name) const + { + return get(name); + } + + void ParameterGroup::readXML(const std::string& xml_filename) + { + fill_xml(*this, xml_filename); + } + + namespace { + inline std::istream& + samcode_readline(std::istream& is, std::string& parameter) { + return std::getline(is, parameter); + } + } // anonymous namespace + + void ParameterGroup::readParam(const std::string& param_filename) { + std::ifstream is(param_filename.c_str()); + if (!is) { + std::cerr << "Error: Failed to open parameter file '" << param_filename << "'.\n"; + throw std::exception(); + } + std::string parameter; + int lineno = 0; + while (samcode_readline(is, parameter)) { + ++lineno; + int fpos = parameter.find(ID_delimiter_assignment); + if (fpos == int(std::string::npos)) { + std::cerr << "WARNING: No '" << ID_delimiter_assignment << "' found on line " << lineno << ".\n"; + } + int pos = fpos + ID_delimiter_assignment.size(); + int spos = parameter.find(ID_delimiter_assignment, pos); + if (spos == int(std::string::npos)) { + std::string name = parameter.substr(0, fpos); + std::string value = parameter.substr(pos, spos); + this->insertParameter(name, value); + } else { + std::cerr << "WARNING: To many '" << ID_delimiter_assignment << "' found on line " << lineno << ".\n"; + } + } + #ifdef MATLAB_MEX_FILE + fclose(is); + #endif + } + + void ParameterGroup::writeParam(const std::string& param_filename) const { + std::ofstream os(param_filename.c_str()); + if (!os) { + std::cerr << "Error: Failed to open parameter file '" + << param_filename << "'.\n"; + throw std::exception(); + } + this->writeParamToStream(os); + } + + void ParameterGroup::writeParamToStream(std::ostream& os) const + { + if (map_.empty()) { + os << this->path() << "/" << "dummy" + << ID_delimiter_assignment << "0\n"; + } + for (map_type::const_iterator it = map_.begin(); it != map_.end(); ++it) { + if ( (*it).second->getTag() == ID_xmltag__param_grp ) { + ParameterGroup& pg = dynamic_cast(*(*it).second); + pg.writeParamToStream(os); + } else if ( (*it).second->getTag() == ID_xmltag__param) { + os << this->path() << "/" << (*it).first << ID_delimiter_assignment + << dynamic_cast(*(*it).second).getValue() << "\n"; + } else { + os << this->path() << "/" << (*it).first << ID_delimiter_assignment + << "<" << (*it).second->getTag() << ">" << "\n"; + } + } + } + + void ParameterGroup::insert(const std::string& name, + const std::tr1::shared_ptr& data) + { + std::pair name_path = split(name); + map_type::const_iterator it = map_.find(name_path.first); + assert(name_path.second == ""); + if (it == map_.end()) { + map_[name] = data; + } else { + if ( (map_[name]->getTag() == data->getTag()) && + (data->getTag() == ID_xmltag__param_grp) ) { + ParameterGroup& alpha = dynamic_cast(*(*it).second); + ParameterGroup& beta = dynamic_cast(*data); + for (map_type::const_iterator + item = beta.map_.begin(); item != beta.map_.end(); ++item) { + alpha.insert((*item).first, (*item).second); + } + } else { + std::cout << "WARNING : The '" + << map_[name]->getTag() + << "' element '" + << name + << "' already exist in group '" + << this->path() + << "'. The element will be replaced by a '" + << data->getTag() + << "' element.\n"; + map_[name] = data; + } + } + } + + void ParameterGroup::insertParameter(const std::string& name, + const std::string& value) + { + std::pair name_path = split(name); + while (name_path.first == "") { + name_path = split(name_path.second); + } + map_type::const_iterator it = map_.find(name_path.first); + if (it == map_.end()) { + if (name_path.second == "") { + std::tr1::shared_ptr data(new Parameter(value, ID_param_type__cmdline)); + map_[name_path.first] = data; + } else { + std::tr1::shared_ptr data(new ParameterGroup(this->path() + ID_delimiter_path + name_path.first, + this)); + ParameterGroup& child = dynamic_cast(*data); + child.insertParameter(name_path.second, value); + map_[name_path.first] = data; + } + } else if (name_path.second == "") { + if ((*it).second->getTag() != ID_xmltag__param) { + std::cout << "WARNING : The " + << map_[name_path.first]->getTag() + << " element '" + << name + << "' already exist in group '" + << this->path() + << "'. It will be replaced by a " + << ID_xmltag__param + << " element.\n"; + } + std::tr1::shared_ptr data(new Parameter(value, ID_param_type__cmdline)); + map_[name_path.first] = data; + } else { + ParameterGroup& pg = dynamic_cast(*(*it).second); + pg.insertParameter(name_path.second, value); + } + } + +#if 0 + void ParameterGroup::display() const { + std::cout << "Begin: " << this->path() << "\n"; + for (map_type::const_iterator it = map_.begin(); it != map_.end(); ++it) { + if ( (*it).second->getTag() == ID_xmltag__param_grp ) { + ParameterGroup& pg = dynamic_cast(*(*it).second); + pg.display(); + } else if ( (*it).second->getTag() == ID_xmltag__param) { + std::cout << (*it).first + << " (" + << (*it).second->getTag() + << ") has value " + << dynamic_cast(*(*it).second).getValue() + << "\n"; + } else { + std::cout << (*it).first + << " (" + << (*it).second->getTag() + << ")\n"; + } + } + std::cout << "End: " << this->path() << "\n"; + } +#endif + + bool ParameterGroup::anyUnused() const + { + if (!this->used()) { + return true; + } + for (map_type::const_iterator it = map_.begin(); it != map_.end(); ++it) { + if (it->second->getTag() == ID_xmltag__param_grp) { + ParameterGroup& pg = dynamic_cast(*(it->second)); + if (pg.anyUnused()) { + return true; + } + } else if (it->second->getTag() == ID_xmltag__param) { + if (!it->second->used()) { + return true; + } + } + } + return false; + } + + void ParameterGroup::displayUsage(bool used_params) const + { + if (this->used() == used_params) { + std::cout << this->path() << '\n'; + } + for (map_type::const_iterator it = map_.begin(); it != map_.end(); ++it) { + if (it->second->getTag() == ID_xmltag__param_grp) { + ParameterGroup& pg = dynamic_cast(*(it->second)); + pg.displayUsage(used_params); + } else if (it->second->getTag() == ID_xmltag__param) { + if (it->second->used() == used_params) { + std::cout << path() << '/' << it->first << '\n'; + } + } + } + std::cout << std::flush; + } + + void ParameterGroup::disableOutput() { + this->recursiveSetIsOutputEnabled(false); + } + + void ParameterGroup::enableOutput() { + this->recursiveSetIsOutputEnabled(true); + } + + bool ParameterGroup::isOutputEnabled() const { + return output_is_enabled_; + } + + void ParameterGroup::recursiveSetIsOutputEnabled(bool output_is_enabled) + { + output_is_enabled_ = output_is_enabled; + for (map_type::const_iterator it = map_.begin(); it != map_.end(); ++it) { + if (it->second->getTag() == ID_xmltag__param_grp) { + ParameterGroup& pg = dynamic_cast(*(it->second)); + pg.recursiveSetIsOutputEnabled(output_is_enabled); + } + } + } + + } // namespace parameter +} // namespace Dune diff --git a/dune/common/param/ParameterGroup.hpp b/dune/common/param/ParameterGroup.hpp new file mode 100644 index 00000000..ecfe0684 --- /dev/null +++ b/dune/common/param/ParameterGroup.hpp @@ -0,0 +1,292 @@ +//=========================================================================== +// +// File: ParameterGroup.hpp +// +// Created: Tue Jun 2 19:11:11 2009 +// +// Author(s): Bård Skaflestad +// Atgeirr F Rasmussen +// +// $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 . +*/ + +#ifndef OPENRS_PARAMETERGROUP_HEADER +#define OPENRS_PARAMETERGROUP_HEADER + +#include +#include +#include +#include +#include + +#include +#include + +namespace Dune { + namespace parameter { + /// ParameterGroup is a class that is used to provide run-time paramters. + /// The standard use of the class is to call create it with the + /// (int argc, char** argv) constructor (where the arguments are those + /// given by main). This parses the command line, where each token + /// either + /// A) specifies a parameter (by a "param=value" token). + /// B) specifies a xml file to be read (by a "filename.xml" token). + /// C) specifies a param file to be read (by a "filename.param" token). + /// After the tokens are parsem they are stored in a tree structure + /// in the ParameterGroup object; it is worth mentioning that parameters + /// are inherited in this tree structure. Thus, if ``grid\_prefix'' is + /// given a value in the root node, this value will be visible in all + /// parts of the tree (unless the parameter is overwritten in a subtree). + /// Applications using this ParameterGroup objects will usually write out + /// a message for each node in the three that is used by the application; + /// this is one way to determine valid parameters. + /// + /// Parameters specified as "param=value" on the command line + /// + /// To specify a parameter on the command line, you must know where in the tree the + /// parameter resides. The syntax for specifying parameters on the command line given + /// an application called ``simulator'' is + /// simulator param1=value grp/param2=value + // for parameter ``param1'' lying at the root and ``param2'' in the group + /// ``grp''. If the same parameters are specified multiple times on the command + /// line, only the last will be used. Thus an application named ``simulator'' run with + /// the following command + /// simulator param=value1 param=value2 + /// will use value2 as the value of ``param''. + /// + /// XML parameters + /// + /// In the introduction to this section it was mentioned that the parameters for + /// the simulator are organized in a tree structure. This is mirrored in the XML + /// file by the use of groups; a group is a collection of parameters and subgroups + /// enclosed by a \fixed{}-\fixed{} pair. The only + /// attribute of relevance in a \fixed{ParameterGroup} tag is that of \fixed{name}, + /// which is used to navigate in the tree. + /// The actual parameters are specified by the \fixed{} tag. Each + /// parameter has three attributes: \fixed{name}, \fixed{type}, and \fixed{value}. + /// Both name and value should be evident. Type is one of \fixed{bool} (for things + /// that are either true or false), \fixed{int} (for integers), \fixed{double} (for + /// floating point numbers with double precision), \fixed{string} (for text + /// strings), or \fixed{file} (for files and directories relative to the location + /// of the XML file). + /// + /// param files + /// + /// A param file consists of multiple lienes, where each line consists of "param=value". + /// This syntax is identical to the one for paramters specified on the command line. + /// + /// + /// If one combines both XML files and parameters, one should note that that if a parameter + /// is specified on the command line is also found in a XML file, the parameter + /// specified on the command line is the one used. Thus, if ``parameters.xml'' + /// specifies that ``stepsize'' is 2.71828 while, the command used to run the + /// application ``simulator'' is + /// simulator stepsize=3.14159 parameters.xml + /// the simulator will run with ``stepsize'' equal to 3.14159. + /// + class ParameterGroup : public ParameterMapItem { + public: + struct NotFoundException : public std::exception {}; + struct WrongTypeException : public std::exception {}; + struct ConversionFailedException : public std::exception {}; + + template + struct RequirementFailedException : public std::exception {}; + + ParameterGroup(); + ParameterGroup(const std::string& path, const ParameterGroup* parent); + + // From ParameterMapItem + virtual ~ParameterGroup(); + virtual std::string getTag() const; + + /// \brief A constructor typically used to initialize a + /// ParameterGroup from command-line arguments. + /// + /// It is required that argv[0] is the program name, while if + /// 0 < i < argc, then argv[i] is either + /// the name of an xml file or parametername=value. + /// + /// \param argc is the number of command-line arguments, + /// including the name of the executable. + /// \param argv is an array of char*, each of which is a + /// command line argument. + template + ParameterGroup(int argc, StringArray argv); + + /// \brief This method checks if there is something with name + /// \p name in the parameter gropup. + /// + /// \param name is the name of the parameter. + /// \return true if \p name is the name of something in the parameter + /// group, false otherwise. + bool has(const std::string& name) const; + + /// \brief This method is used to read a parameter from the + /// parameter group. + /// + /// NOTE: If the reading of the parameter fails, then this method + /// throws an appropriate exception. + /// + /// \param name is the name of the parameter in question. + /// \return The value associated with then name in this parameter + /// group. + template + T get(const std::string& name) const; + + template + T get(const std::string& name, const Requirement&) const; + + /// \brief This method is used to read a parameter from the + /// parameter group. + /// + /// NOTE: If the reading of the parameter fails, then either + /// a) this method returns \p default_value if there + /// was no parameter with name \p name in + /// the parameter group + /// or + /// b) this method throws an appropriate exception. + /// + /// \param name is the name of the parameter in question. + /// \param default_value the defualt value of the parameter in + /// question. + /// \return The value associated with then name in this parameter + /// group. + template + T getDefault(const std::string& name, + const T& default_value) const; + + template + T getDefault(const std::string& name, + const T& default_value, + const Requirement& r) const; + + /// \brief This method returns the parameter group given by name, + /// i.e. it is an alias of get(). + /// + /// \param name is the name of the parameter group sought. + /// \return the parameter group sought. + ParameterGroup getGroup(const std::string& name) const; + + /// \brief Disables the output from get, getDefault and getGroup. + /// By default, such output is enabled. + void disableOutput(); + + /// \brief Enables the output from get, getDefault and getGroup. + /// By default, such output is enabled. + void enableOutput(); + + /// \brief Returs true if and only if output from get, getDefault + /// and getGroup is enabled. + /// + /// \return true if and only if output from get, getDefault and + /// getGroup is enabled. + bool isOutputEnabled() const; + + + /// \brief Reads the contents of the xml file specified by + /// xml_filename into this ParameterGroup. + /// + /// \param xml_filename is the name of a xml file. + void readXML(const std::string& xml_filename); + + /// \brief Reads the contents of the param file specified by + /// param_filename into this ParameterGroup. + /// + /// NOTE: A param file contains lines on the form 'a/b/c=d'. + // The '/' separates ParameterGroups and Parameters + /// (e.g. c is a Parameter in the ParameterGroup b, + /// and b is a ParameterGroup in the ParameterGroup a) + /// while '=' separates the name from the value (e.g. the + /// value of the parameter c above will be d). + /// NOTE: A param file does not store any type information about + /// its values. + /// + /// \param param_filename is the name of a param file. + void readParam(const std::string& param_filename); + + /// \brief Writes this ParameterGroup into a param file + /// specified by param_filename. + /// + /// \param param_filename is the name of a param file. + void writeParam(const std::string& param_filename) const; + + /// \brief Writes this ParameterGroup to a stream. + /// + /// \param stream is the stream to write to. + void writeParamToStream(std::ostream& stream) const; + + + /// vki param interface - deprecated + template + void get(const char* name, T& value, const T& default_value) const { + value = this->getDefault(name, default_value); + } + + /// vki param interface - deprecated + template + void get(const char* name, T& value) const { + value = this->get(name); + } + + /// \brief Return true if any parameters are unused. + bool anyUnused() const; + + /// \brief Shows which parameters which are used or unused. + void displayUsage(bool used_params = false) const; + + /// \brief Returns the path of the parameter group. + std::string path() const; + + /// Insert a new item into the group. + void insert(const std::string& name, + const std::tr1::shared_ptr& data); + + /// Insert a new parameter item into the group. + void insertParameter(const std::string& name, const std::string& value); + + private: + typedef std::tr1::shared_ptr data_type; + typedef std::pair pair_type; + typedef std::map map_type; + + std::string path_; + map_type map_; + const ParameterGroup* parent_; + bool output_is_enabled_; + + template + T translate(const pair_type& data, const Requirement& chk) const; + template + void parseCommandLineArguments(int argc, StringArray argv); + void recursiveSetIsOutputEnabled(bool output_is_enabled); + }; + } // namespace parameter +} // namespace Dune + +#include + +#endif // OPENRS_PARAMETERGROUP_HEADER diff --git a/dune/common/param/ParameterGroup_impl.hpp b/dune/common/param/ParameterGroup_impl.hpp new file mode 100644 index 00000000..3ef77a72 --- /dev/null +++ b/dune/common/param/ParameterGroup_impl.hpp @@ -0,0 +1,330 @@ +//=========================================================================== +// +// File: ParameterGroup_impl.hpp +// +// Created: Tue Jun 2 19:06:46 2009 +// +// Author(s): Bård Skaflestad +// Atgeirr F Rasmussen +// +// $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 . +*/ + +#ifndef OPENRS_PARAMETERGROUP_IMPL_HEADER +#define OPENRS_PARAMETERGROUP_IMPL_HEADER + +#include +#include + +//#include "TermColors.hpp" // from utils + +#include +#include +#include + +namespace Dune { + namespace parameter { + + template<> + struct ParameterMapItemTrait { + static ParameterGroup + convert(const ParameterMapItem& item, + std::string& conversion_error) + { + std::string tag = item.getTag(); + if (tag != ID_xmltag__param_grp) { + conversion_error = "The XML tag was '" + tag + + "' but should be '" + + ID_xmltag__param_grp + "'.\n"; + return ParameterGroup("", 0); + } + conversion_error = ""; + const ParameterGroup& pg = dynamic_cast(item); + return pg; + } + static std::string type() {return "ParameterGroup";} + }; + + namespace { + template + inline std::string + to_string(const T& val) + { + std::ostringstream os; + os << val; + return os.str(); + } + + inline std::string + to_string(const bool b) { + if (b) { + return ID_true; + } else { + return ID_false; + } + } + + inline std::string + to_string(const ParameterGroup&) + { + return std::string(""); + } + + std::pair + filename_split(const std::string& filename) + { + int fpos = filename.rfind('.'); + std::string name = filename.substr(0, fpos); + std::string type = filename.substr(fpos+1); + return std::make_pair(name, type); + } + + } + + + + template + ParameterGroup::ParameterGroup(int argc, StringArray argv) + : path_(ID_path_root), parent_(0), output_is_enabled_(true) + { + if (argc < 2) { + std::cerr << "Usage: " << argv[0] << " " + << "[paramfilename1.{xml,param}] " + << "[paramfilename2.{xml,param}] " + << "[overridden_arg1=value1] " + << "[overridden_arg2=value2] " + << "[...]" << std::endl; + exit(EXIT_FAILURE); + } + this->parseCommandLineArguments(argc, argv); + } + + template + void ParameterGroup::parseCommandLineArguments(int argc, StringArray argv) + { + std::vector files; + std::vector > assignments; + for (int i = 1; i < argc; ++i) { + std::string arg(argv[i]); + int fpos = arg.find(ID_delimiter_assignment); + if (fpos == int(std::string::npos)) { + std::string filename = arg.substr(0, fpos); + files.push_back(filename); + continue; + } + int pos = fpos + ID_delimiter_assignment.size(); + int spos = arg.find(ID_delimiter_assignment, pos); + if (spos == int(std::string::npos)) { + std::string name = arg.substr(0, fpos); + std::string value = arg.substr(pos, spos); + assignments.push_back(std::make_pair(name, value)); + continue; + } + std::cout << "WARNING: Too many assignements (' " + << ID_delimiter_assignment + << "') detected in argument " << i << ".\n"; + } + for (int i = 0; i < int(files.size()); ++i) { + std::pair file_type = filename_split(files[i]); + if (file_type.second == "xml") { + this->readXML(files[i]); + } else if (file_type.second == "param") { + this->readParam(files[i]); + } else { + std::cout << "WARNING: Ignoring file '" + << files[i] << "' with unknown extension.\n" + << "Valid filename extensions are 'xml' and 'param'.\n"; + } + } + for (int i = 0; i < int(assignments.size()); ++i) { + this->insertParameter(assignments[i].first, assignments[i].second); + } + } + + + template + inline T ParameterGroup::get(const std::string& name) const + { + return this->get(name, ParameterRequirementNone()); + } + + template + inline T ParameterGroup::get(const std::string& name, + const Requirement& r) const + { + setUsed(); + std::pair name_path = split(name); + map_type::const_iterator it = map_.find(name_path.first); + if (it == map_.end()) { + if (parent_ != 0) { + // If we have a parent, ask it instead. + if (output_is_enabled_) { + //TermColors::Red(); + std::cout << name; + //TermColors::Normal(); + std::cout << " not found at " + << (path() + ID_delimiter_path) + << ", asking parent." << std::endl; + } + return parent_->get(name, r); + } else { + // We are at the top, name has not been found. + std::cerr << "ERROR: The group '" + << this->path() + << "' does not contain an element named '" + << name + << "'.\n"; + throw NotFoundException(); + } + } + if (name_path.second == "") { + T val = this->translate(*it, r); + it->second->setUsed(); + if (output_is_enabled_) { + //TermColors::Green(); + std::cout << name; + //TermColors::Normal(); + std::cout << " found at " << (path() + ID_delimiter_path) + << ", value is " << to_string(val) << std::endl; + } + return val; + } else { + ParameterGroup& pg = dynamic_cast(*(*it).second); + pg.setUsed(); + return pg.get(name_path.second, r); + } + } + + template + inline T ParameterGroup::getDefault(const std::string& name, + const T& default_value) const + { + return this->getDefault(name, default_value, ParameterRequirementNone()); + } + + template + inline T ParameterGroup::getDefault(const std::string& name, + const T& default_value, + const Requirement& r) const + { + setUsed(); + std::pair name_path = split(name); + map_type::const_iterator it = map_.find(name_path.first); + if (it == map_.end()) { + if (parent_ != 0) { + // If we have a parent, ask it instead. + if (output_is_enabled_) { + //TermColors::Red(); + std::cout << name; + //TermColors::Normal(); + std::cout << " not found at " << (path() + ID_delimiter_path) + << ", asking parent." << std::endl; + } + return parent_->getDefault(name, default_value, r); + } else { + // We check the requirement for the default value + std::string requirement_result = r(default_value); + if (requirement_result != "") { + std::cerr << "ERROR: The default value for the " + << " element named '" + << name + << "' in the group '" + << this->path() + << "' failed to meet a requirenemt.\n"; + std::cerr << "The requirement enforcer returned the following message:\n" + << requirement_result + << "\n"; + throw RequirementFailedException(); + } + } + if (output_is_enabled_) { + //TermColors::Blue(); + std::cout << name; + //TermColors::Normal(); + std::cout << " not found. Using default value '" + << to_string(default_value) << "'." << std::endl; + } + return default_value; + } + if (name_path.second == "") { + T val = this->translate(*it, r); + it->second->setUsed(); + if (output_is_enabled_) { + //TermColors::Green(); + std::cout << name; + //TermColors::Normal(); + std::cout << " found at " << (path() + ID_delimiter_path) + << ", value is '" << to_string(val) << "'." << std::endl; + } + return val; + } else { + ParameterGroup& pg = dynamic_cast(*(*it).second); + pg.setUsed(); + return pg.getDefault(name_path.second, default_value, r); + } + } + + template + inline T ParameterGroup::translate(const pair_type& named_data, + const Requirement& chk) const + { + const std::string& name = named_data.first; + const data_type data = named_data.second; + std::string conversion_error; + T value = ParameterMapItemTrait::convert(*data, conversion_error); + if (conversion_error != "") { + std::cerr << "ERROR: Failed to convert the element named '" + << name + << "' in the group '" + << this->path() + << "' to the type '" + << ParameterMapItemTrait::type() + << "'.\n"; + std::cerr << "The conversion routine returned the following message:\n" + << conversion_error + << "\n"; + throw WrongTypeException(); + } + std::string requirement_result = chk(value); + if (requirement_result != "") { + std::cerr << "ERROR: The element named '" + << name + << "' in the group '" + << this->path() + << "' of type '" + << ParameterMapItemTrait::type() + << "' failed to meet a requirenemt.\n"; + std::cerr << "The requirement enforcer returned the following message:\n" + << requirement_result + << "\n"; + throw RequirementFailedException(); + } + return value; + } + } // namespace parameter +} // namespace Dune + +#endif // OPENRS_PARAMETERGROUP_IMPL_HEADER diff --git a/dune/common/param/ParameterMapItem.hpp b/dune/common/param/ParameterMapItem.hpp new file mode 100644 index 00000000..24248e61 --- /dev/null +++ b/dune/common/param/ParameterMapItem.hpp @@ -0,0 +1,73 @@ +//=========================================================================== +// +// File: ParameterMapItem.hpp +// +// Created: Tue Jun 2 19:05:54 2009 +// +// Author(s): Bård Skaflestad +// Atgeirr F Rasmussen +// +// $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 . +*/ + +#ifndef OPENRS_PARAMETERMAPITEM_HEADER +#define OPENRS_PARAMETERMAPITEM_HEADER + +#include + +namespace Dune { + namespace parameter { + /// The parameter handlig system is structured as a tree, + /// where each node inhertis from ParameterMapItem. + /// + /// The abstract virtual function getTag() is used to determine + /// which derived class the node actually is. + struct ParameterMapItem { + /// Default constructor + ParameterMapItem() : used_(false) {} + + /// Destructor + virtual ~ParameterMapItem() {} + + /// \brief This function returns a string describing + /// the ParameterMapItem. + virtual std::string getTag() const = 0; + void setUsed() const { used_ = true; } + bool used() const { return used_; } + private: + mutable bool used_; + }; + + template + struct ParameterMapItemTrait { + static T convert(const ParameterMapItem&, + std::string& conversion_error); + static std::string type(); + }; + } // namespace parameter +} // namespace Dune + +#endif // OPENRS_PARAMETERMAPITEM_HEADER diff --git a/dune/common/param/ParameterRequirement.hpp b/dune/common/param/ParameterRequirement.hpp new file mode 100644 index 00000000..70d12a38 --- /dev/null +++ b/dune/common/param/ParameterRequirement.hpp @@ -0,0 +1,243 @@ +//=========================================================================== +// +// File: ParameterRequirement.hpp +// +// Created: Tue Jun 2 19:05:02 2009 +// +// Author(s): Bård Skaflestad +// Atgeirr F Rasmussen +// +// $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 . +*/ + +#ifndef OPENRS_PARAMETERREQUIREMENT_HEADER +#define OPENRS_PARAMETERREQUIREMENT_HEADER + +#include +#include +#include +#include +#include + +namespace Dune { + namespace parameter { + /// @brief + /// @todo Doc me! + /// @tparam + /// @param + /// @return + struct ParameterRequirementNone { + template + std::string operator()(const T&) const { + return ""; + } + }; + + /// @brief + /// @todo Doc me! + /// @param + /// @return + struct ParameterRequirementProbability { + std::string operator()(double x) const { + if ( (x < 0.0) || (x > 1.0) ) { + std::ostringstream stream; + stream << "The value '" << x + << "' is not in the interval [0, 1]."; + return stream.str(); + } else { + return ""; + } + } + }; + + /// @brief + /// @todo Doc me! + /// @tparam + /// @param + /// @return + struct ParameterRequirementPositive { + template + std::string operator()(const T& x) const { + if (x > 0) { + return ""; + } else { + std::ostringstream stream; + stream << "The value '" << x << "' is not positive."; + return stream.str(); + } + } + }; + + /// @brief + /// @todo Doc me! + /// @tparam + /// @param + /// @return + struct ParameterRequirementNegative { + template + std::string operator()(const T& x) const { + if (x < 0) { + return ""; + } else { + std::ostringstream stream; + stream << "The value '" << x << "' is not negative."; + return stream.str(); + } + } + }; + + /// @brief + /// @todo Doc me! + /// @tparam + /// @param + /// @return + struct ParameterRequirementNonPositive { + template + std::string operator()(const T& x) const { + if (x > 0) { + std::ostringstream stream; + stream << "The value '" << x << "' is positive."; + return stream.str(); + } else { + return ""; + } + } + }; + + /// @brief + /// @todo Doc me! + /// @tparam + /// @param + /// @return + struct ParameterRequirementNonNegative { + template + std::string operator()(const T& x) const { + if (x < 0) { + std::ostringstream stream; + stream << "The value '" << x << "' is negative."; + return stream.str(); + } else { + return ""; + } + } + }; + + /// @brief + /// @todo Doc me! + /// @tparam + /// @param + /// @return + struct ParameterRequirementNonZero { + template + std::string operator()(const T& x) const { + if (x != 0) { + return ""; + } else { + return "The value was zero."; + } + } + }; + + /// @brief + /// @todo Doc me! + /// @param + /// @return + struct ParameterRequirementNonEmpty { + std::string operator()(const std::string& x) const { + if (x != "") { + return "The string was empty."; + } else { + return ""; + } + } + }; + + /// @brief + /// @todo Doc me! + /// @tparam + /// @param + /// @return + template + struct ParameterRequirementAnd { + ParameterRequirementAnd(const Requirement1& r1, const Requirement2& r2) : + r1_(r1), r2_(r2) { } + + template + std::string operator()(const T& t) const { + std::string e1 = r1_(t); + std::string e2 = r2_(t); + if (e1 == "") { + return e2; + } else if (e2 == "") { + return e1; + } else { + return e1 + " AND " + e2; + } + } + private: + const Requirement1 r1_; + const Requirement2 r2_; + }; + + /// @brief + /// @todo Doc me! + /// @param + struct ParameterRequirementMemberOf { + ParameterRequirementMemberOf(const std::vector& elements) + : elements_(elements) { + assert(elements_.size() > 0); + } + + /// @brief + /// @todo Doc me! + /// @param + /// @return + std::string operator()(const std::string& x) const { + if (std::find(elements_.begin(), elements_.end(), x) == elements_.end()) { + if (elements_.size() == 1) { + return "The string '" + x + "' is not '" + elements_[0] + "'."; + } + std::ostringstream stream; + stream << "The string '" << x << "' is not among '"; + for (int i = 0; i < int(elements_.size()) - 2; ++i) { + stream << elements_[i] << "', '"; + } + stream << elements_[elements_.size() - 2] + << "' and '" + << elements_[elements_.size() - 1] + << "'."; + return stream.str(); + } else { + return ""; + } + } + private: + const std::vector elements_; + }; + } // namespace parameter +} // namespace Dune + +#endif // OPENRS_PARAMETERREQUIREMENT_HEADER diff --git a/dune/common/param/ParameterStrings.hpp b/dune/common/param/ParameterStrings.hpp new file mode 100644 index 00000000..2b61fdb6 --- /dev/null +++ b/dune/common/param/ParameterStrings.hpp @@ -0,0 +1,70 @@ +//=========================================================================== +// +// File: ParameterStrings.hpp +// +// Created: Tue Jun 2 19:04:15 2009 +// +// Author(s): Bård Skaflestad +// Atgeirr F Rasmussen +// +// $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 . +*/ + +#ifndef OPENRS_PARAMETERSTRINGS_HEADER +#define OPENRS_PARAMETERSTRINGS_HEADER + +#include + +namespace Dune { + namespace parameter { + const std::string ID_true = "true"; + const std::string ID_false = "false"; + + const std::string ID_xmltag__param_grp = "ParameterGroup"; + const std::string ID_xmltag__param = "Parameter"; + const std::string ID_xmltag__file_param_grp = "ParameterGroupFromFile"; + const std::string ID_xmltag__file_params = "MergeWithFile"; + + const std::string ID_xmlatt__value = "value"; + const std::string ID_xmlatt__name = "name"; + const std::string ID_xmlatt__type = "type"; + + const std::string ID_param_type__bool = "bool"; + const std::string ID_param_type__int = "int"; + const std::string ID_param_type__float = "double"; + const std::string ID_param_type__string = "string"; + const std::string ID_param_type__file = "file"; + const std::string ID_param_type__cmdline = "cmdline"; + + // + + const std::string ID_path_root = ""; + const std::string ID_delimiter_path = "/"; + const std::string ID_delimiter_assignment = "="; + } // namespace parameter +} // namespace Dune + +#endif // OPENRS_PARAMETERSTRINGS_HEADER diff --git a/dune/common/param/ParameterTools.cpp b/dune/common/param/ParameterTools.cpp new file mode 100644 index 00000000..f034be71 --- /dev/null +++ b/dune/common/param/ParameterTools.cpp @@ -0,0 +1,55 @@ +//=========================================================================== +// +// File: ParameterTools.cpp +// +// Created: Tue Jun 2 19:03:09 2009 +// +// Author(s): Bård Skaflestad +// Atgeirr F Rasmussen +// +// $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 . +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif +#include +#include + +namespace Dune { + namespace parameter { + std::pair split(const std::string& name) + { + int pos = name.find(ID_delimiter_path); + if (pos == int(std::string::npos)) { + return std::make_pair(name, ""); + } else { + return std::make_pair(name.substr(0, pos), + name.substr(pos + ID_delimiter_path.size())); + } + } + } // namespace parameter +} // namespace Dune diff --git a/dune/common/param/ParameterTools.hpp b/dune/common/param/ParameterTools.hpp new file mode 100644 index 00000000..05701f96 --- /dev/null +++ b/dune/common/param/ParameterTools.hpp @@ -0,0 +1,48 @@ +//=========================================================================== +// +// File: ParameterTools.hpp +// +// Created: Tue Jun 2 19:02:19 2009 +// +// Author(s): Bård Skaflestad +// Atgeirr F Rasmussen +// +// $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 . +*/ + +#ifndef OPENRS_PARAMETERTOOLS_HEADER +#define OPENRS_PARAMETERTOOLS_HEADER + +#include +#include + +namespace Dune { + namespace parameter { + std::pair split(const std::string& name); + } // namespace parameter +} // namespace Dune + +#endif // OPENRS_PARAMETERTOOLS_HEADER diff --git a/dune/common/param/ParameterXML.cpp b/dune/common/param/ParameterXML.cpp new file mode 100644 index 00000000..3147cf9c --- /dev/null +++ b/dune/common/param/ParameterXML.cpp @@ -0,0 +1,140 @@ +//=========================================================================== +// +// File: ParameterXML.cpp +// +// Created: Tue Jun 2 18:57:25 2009 +// +// Author(s): Bård Skaflestad +// Atgeirr F Rasmussen +// +// $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 . +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif +#include +#include +#include + +#include +#include +#include +#include + +#include // libxml2 +#include // libxml2 + + +namespace Dune { + namespace parameter { + + namespace libxml2 { + void read_xml(ParameterGroup& pg, const std::string filename); + void fill_tree(ParameterGroup& pg, xmlNodePtr root, + const std::string& xml_file_path); + } + + void fill_xml(ParameterGroup& pg, const std::string filename) { + libxml2::read_xml(pg, filename); + } + + namespace libxml2 { + std::string getProperty(const std::string& property, + const xmlNodePtr& node_ptr) + { + const char* prop_value_ptr = (const char*) xmlGetProp(node_ptr, (const xmlChar*)property.c_str()); + std::string property_value(prop_value_ptr ? prop_value_ptr : ""); + xmlFree((xmlChar*) prop_value_ptr); + return property_value; + } + + void read_xml(ParameterGroup& pg, const std::string filename) + { + // Macro to check that the libxml version in use is compatible + // with the version the software has been compiled against + LIBXML_TEST_VERSION; // defined by libxml2 + + xmlDocPtr doc = xmlParseFile(filename.c_str()); + if (doc == NULL) { + xmlFreeDoc(doc); + xmlCleanupParser(); + std::cerr << "ERROR: Failed to open XML file '" << filename << "'\n"; + throw std::exception(); + } + xmlNodePtr root = xmlDocGetRootElement(doc); + std::string xml_file_path = boost::filesystem::path(filename).branch_path().string(); + fill_tree(pg, root, xml_file_path); + xmlFreeDoc(doc); + xmlCleanupParser(); + } + + void fill_tree(ParameterGroup& pg, xmlNodePtr root, + const std::string& xml_file_path) + { + //std::cout << "GROUP '" << value << "' BEGIN\n"; + for (xmlNodePtr cur = root->children; cur; cur = cur->next) { + if (cur->type == XML_ELEMENT_NODE) { + const char* tag_ptr = (const char*) cur->name; + std::string tag_name(tag_ptr ? tag_ptr : ""); + if (tag_name == ID_xmltag__file_params) { + std::string relative_filename = getProperty(ID_xmlatt__value, cur); + std::string filename = (boost::filesystem::path(xml_file_path) / relative_filename).string(); + fill_xml(pg, filename); + continue; + } + std::string name = getProperty(ID_xmlatt__name, cur); + std::tr1::shared_ptr data; + if (tag_name == ID_xmltag__param) { + std::string value = getProperty(ID_xmlatt__value, cur); + std::string type = getProperty(ID_xmlatt__type, cur); + if (type == ID_param_type__file) { + value = (boost::filesystem::path(xml_file_path) / value).string(); + type = ID_param_type__string; + } + data.reset(new Parameter(value, type)); + } else if (tag_name == ID_xmltag__param_grp) { + std::string child_path = pg.path() + ID_delimiter_path + name; + data.reset(new ParameterGroup(child_path, &pg)); + fill_tree(dynamic_cast(*data), cur, xml_file_path); + } else if (tag_name == ID_xmltag__file_param_grp) { + std::string child_path = pg.path() + ID_delimiter_path + name; + data.reset(new ParameterGroup(child_path, &pg)); + std::string relative_filename = getProperty(ID_xmlatt__value, cur); + std::string filename = (boost::filesystem::path(xml_file_path) / relative_filename).string(); + fill_xml(dynamic_cast(*data), filename); + } else { + std::cerr << "ERROR: '" << tag_name << "' is an unknown xml tag.\n"; + throw std::exception(); + } + pg.insert(name, data); + } + } + //std::cout << "GROUP '" << id << "' END\n"; + } + } + } // namespace parameter +} // namespace Dune diff --git a/dune/common/param/ParameterXML.hpp b/dune/common/param/ParameterXML.hpp new file mode 100644 index 00000000..0e5b58dd --- /dev/null +++ b/dune/common/param/ParameterXML.hpp @@ -0,0 +1,54 @@ +//=========================================================================== +// +// File: ParameterXML.hpp +// +// Created: Tue Jun 2 18:55:59 2009 +// +// Author(s): Bård Skaflestad +// Atgeirr F Rasmussen +// +// $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 . +*/ + +#ifndef OPENRS_PARAMETERXML_HEADER +#define OPENRS_PARAMETERXML_HEADER + +#include + +namespace Dune { + namespace parameter { + /// \brief This function fills a ParameterGroup with a XML file. + /// + /// NOTE: The function throws if there is an error during reading + /// of the XML file. + /// + /// \retval pg is the ParameterGroup to be filled. + /// \param filename is the name of an XML file. + void fill_xml(ParameterGroup& pg, const std::string filename); + } // namespace parameter +} // namespace Dune + +#endif // OPENRS_PARAMETERXML_HEADER diff --git a/dune/common/param/test/Makefile.am b/dune/common/param/test/Makefile.am new file mode 100644 index 00000000..7a561dc9 --- /dev/null +++ b/dune/common/param/test/Makefile.am @@ -0,0 +1,14 @@ +# $Date: 2009-12-08 15:21:36 +0100 (Tue, 08 Dec 2009) $ +# $Revision: 699 $ + +check_PROGRAMS = param_test + +param_test_SOURCES = param_test.cpp +param_test_CXXFLAGS = $(DUNEMPICPPFLAGS) $(BOOST_CPPFLAGS) +param_test_LDADD = $(DUNE_LDFLAGS) $(DUNEMPILDFLAGS) $(DUNEMPILIBS) $(BOOST_LDFLAGS) \ + $(DUNE_LIBS) $(BOOST_UNIT_TEST_FRAMEWORK_LIB) ../libparam.la + +TESTS = $(check_PROGRAMS) + +include $(top_srcdir)/am/global-rules + diff --git a/dune/common/param/test/param_test.cpp b/dune/common/param/test/param_test.cpp new file mode 100644 index 00000000..1c9a06d2 --- /dev/null +++ b/dune/common/param/test/param_test.cpp @@ -0,0 +1,66 @@ +//=========================================================================== +// +// File: param_test.cpp +// +// Created: Sun Dec 13 20:08:36 2009 +// +// Author(s): Atgeirr F Rasmussen +// Bård Skaflestad +// +// $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 . +*/ + + +#define BOOST_TEST_DYN_LINK +#define NVERBOSE // to suppress our messages when throwing + +#define BOOST_TEST_MODULE ParameterTest +#include + +#include "../ParameterGroup.hpp" +#include + +using namespace Dune; + +BOOST_AUTO_TEST_CASE(commandline_syntax_init) +{ + typedef const char* cp; + cp argv[] = { "program_command", + "topitem=somestring", + "/slashtopitem=anotherstring", + "/group/item=1", + "/group/anotheritem=2", + "/group/subgroup/item=3", + "/group/subgroup/anotheritem=4", + "/group/item=overridingstring" }; + const std::size_t argc = sizeof(argv)/sizeof(argv[0]); + parameter::ParameterGroup p(argc, argv); + BOOST_CHECK(p.get("topitem") == "somestring"); + + // Tests that only run in debug mode. +#ifndef NDEBUG +#endif +} diff --git a/dune/common/test/Makefile.am b/dune/common/test/Makefile.am new file mode 100644 index 00000000..1d20fee8 --- /dev/null +++ b/dune/common/test/Makefile.am @@ -0,0 +1,29 @@ +# $Date$ +# $Revision$ + +check_PROGRAMS = sparsetable_test sparsevector_test monotcubicinterpolator_test +noinst_PROGRAMS = unit_test + +sparsetable_test_SOURCES = sparsetable_test.cpp +sparsetable_test_CXXFLAGS = $(DUNEMPICPPFLAGS) $(BOOST_CPPFLAGS) +sparsetable_test_LDADD = $(DUNE_LDFLAGS) $(DUNEMPILDFLAGS) $(BOOST_LDFLAGS) \ + $(DUNE_LIBS) $(DUNEMPILIBS) $(BOOST_UNIT_TEST_FRAMEWORK_LIB) + +sparsevector_test_SOURCES = sparsevector_test.cpp +sparsevector_test_CXXFLAGS = $(DUNEMPICPPFLAGS) $(BOOST_CPPFLAGS) +sparsevector_test_LDADD = $(DUNE_LDFLAGS) $(DUNEMPILDFLAGS) $(BOOST_LDFLAGS) \ + $(DUNE_LIBS) $(DUNEMPILIBS) $(BOOST_UNIT_TEST_FRAMEWORK_LIB) + +monotcubicinterpolator_test_SOURCES = monotcubicinterpolator_test.cpp +monotcubicinterpolator_test_CXXFLAGS = $(DUNEMPICPPFLAGS) $(BOOST_CPPFLAGS) +monotcubicinterpolator_test_LDADD = $(DUNE_LDFLAGS) $(DUNEMPILDFLAGS) \ + $(DUNE_LIBS) $(DUNEMPILIBS) ../libcommon.la + +unit_test_SOURCES = unit_test.cpp +unit_test_CXXFLAGS = $(DUNEMPICPPFLAGS) $(BOOST_CPPFLAGS) +unit_test_LDADD = $(DUNE_LDFLAGS) $(DUNEMPILDFLAGS) $(DUNE_LIBS) $(DUNEMPILIBS) + +TESTS = $(check_PROGRAMS) + +include $(top_srcdir)/am/global-rules + diff --git a/dune/common/test/monotcubicinterpolator_test.cpp b/dune/common/test/monotcubicinterpolator_test.cpp new file mode 100644 index 00000000..7e7a86d7 --- /dev/null +++ b/dune/common/test/monotcubicinterpolator_test.cpp @@ -0,0 +1,53 @@ +//=========================================================================== +// +// File: monotcubicinterpolator_test.cpp +// +// Created: Tue Dec 8 12:25:30 2009 +// +// Author(s): Atgeirr F Rasmussen +// Bård Skaflestad +// +// $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 . +*/ + +#include + +int main() +{ + const int num_v = 3; + double xv[num_v] = {0.0, 1.0, 2.0}; + double fv[num_v] = {10.0, 21.0, 2.0}; + std::vector x(xv, xv + num_v); + std::vector f(fv, fv + num_v); + MonotCubicInterpolator interp(x, f); + interp.evaluate(-1.0); + interp.evaluate(0.0); + interp.evaluate(0.0001); + interp.evaluate(0.5); + interp.evaluate(1.0); + interp.evaluate(2.0); + interp.evaluate(4.0); +} diff --git a/dune/common/test/sparsetable_test.cpp b/dune/common/test/sparsetable_test.cpp new file mode 100644 index 00000000..76bc7938 --- /dev/null +++ b/dune/common/test/sparsetable_test.cpp @@ -0,0 +1,128 @@ +//=========================================================================== +// +// File: test_sparsetable.cpp +// +// Created: Thu May 28 10:01:46 2009 +// +// Author(s): Atgeirr F Rasmussen +// Bård Skaflestad +// +// $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 . +*/ + +#define BOOST_TEST_DYN_LINK +#define NVERBOSE // to suppress our messages when throwing + +#define BOOST_TEST_MODULE SparseTableTest +#include + +#include "../SparseTable.hpp" + +using namespace Dune; + +BOOST_AUTO_TEST_CASE(construction_and_queries) +{ + const SparseTable st1; + BOOST_CHECK(st1.empty()); + BOOST_CHECK_EQUAL(st1.size(), 0); + BOOST_CHECK_EQUAL(st1.dataSize(), 0); + + // This should be getting us a table like this: + // ---------------- + // 0 + // + // 1 2 + // 3 4 5 6 + // 7 8 9 + // ---------------- + const int num_elem = 10; + const int elem[num_elem] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 }; + const int num_rows = 5; + const int rowsizes[num_rows] = { 1, 0, 2, 4, 3 }; + const SparseTable st2(elem, elem + num_elem, rowsizes, rowsizes + num_rows); + BOOST_CHECK(!st2.empty()); + BOOST_CHECK_EQUAL(st2.size(), num_rows); + BOOST_CHECK_EQUAL(st2.dataSize(), num_elem); + BOOST_CHECK_EQUAL(st2[0][0], 0); + BOOST_CHECK_EQUAL(st2.rowSize(0), 1); + BOOST_CHECK(st2[1].empty()); + BOOST_CHECK_EQUAL(st2.rowSize(1), 0); + BOOST_CHECK_EQUAL(st2[3][1], 4); + BOOST_CHECK_EQUAL(st2[4][2], 9); + BOOST_CHECK(st2[4].size() == rowsizes[4]); + const SparseTable st2_again(elem, elem + num_elem, rowsizes, rowsizes + num_rows); + BOOST_CHECK(st2 == st2_again); + SparseTable st2_byassign; + st2_byassign.assign(elem, elem + num_elem, rowsizes, rowsizes + num_rows); + BOOST_CHECK(st2 == st2_byassign); + const int last_row_size = rowsizes[num_rows - 1]; + SparseTable st2_append(elem, elem + num_elem - last_row_size, rowsizes, rowsizes + num_rows - 1); + BOOST_CHECK_EQUAL(st2_append.dataSize(), num_elem - last_row_size); + st2_append.appendRow(elem + num_elem - last_row_size, elem + num_elem); + BOOST_CHECK(st2 == st2_append); + SparseTable st2_append2; + st2_append2.appendRow(elem, elem + 1); + st2_append2.appendRow(elem + 1, elem + 1); + st2_append2.appendRow(elem + 1, elem + 3); + st2_append2.appendRow(elem + 3, elem + 7); + st2_append2.appendRow(elem + 7, elem + 10); + BOOST_CHECK(st2 == st2_append2); + st2_append2.clear(); + SparseTable st_empty; + BOOST_CHECK(st2_append2 == st_empty); + + SparseTable st2_allocate; + st2_allocate.allocate(rowsizes, rowsizes + num_rows); + BOOST_CHECK_EQUAL(st2_allocate.size(), num_rows); + BOOST_CHECK_EQUAL(st2_allocate.dataSize(), num_elem); + int s = 0; + for (int i = 0; i < num_rows; ++i) { + SparseTable::mutable_row_type row = st2_allocate[i]; + for (int j = 0; j < rowsizes[i]; ++j, ++s) + row[j] = elem[s]; + } + BOOST_CHECK(st2 == st2_allocate); + + // One element too few. + BOOST_CHECK_THROW(const SparseTable st3(elem, elem + num_elem - 1, rowsizes, rowsizes + num_rows), std::exception); + + // A few elements too many. + BOOST_CHECK_THROW(const SparseTable st4(elem, elem + num_elem, rowsizes, rowsizes + num_rows - 1), std::exception); + + // Need at least one row. + BOOST_CHECK_THROW(const SparseTable st5(elem, elem + num_elem, rowsizes, rowsizes), std::exception); + + // Tests that only run in debug mode. +#ifndef NDEBUG + // Do not ask for wrong row numbers. + BOOST_CHECK_THROW(st1.rowSize(0), std::exception); + BOOST_CHECK_THROW(st2.rowSize(-1), std::exception); + BOOST_CHECK_THROW(st2.rowSize(st2.size()), std::exception); + // No negative row sizes. + const int err_rs[num_rows] = { 1, 0, -1, 7, 3 }; + BOOST_CHECK_THROW(const SparseTable st6(elem, elem + num_elem, err_rs, err_rs + num_rows), std::exception); +#endif +} diff --git a/dune/common/test/sparsevector_test.cpp b/dune/common/test/sparsevector_test.cpp new file mode 100644 index 00000000..e6ec98bf --- /dev/null +++ b/dune/common/test/sparsevector_test.cpp @@ -0,0 +1,106 @@ +//=========================================================================== +// +// File: sparsevector_test.cpp +// +// Created: Mon Jun 29 21:00:53 2009 +// +// Author(s): Atgeirr F Rasmussen +// Bård Skaflestad +// +// $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 . +*/ + + +#define BOOST_TEST_DYN_LINK +#define NVERBOSE // to suppress our messages when throwing + +#define BOOST_TEST_MODULE SparseVectorTest +#include + +#include "../SparseVector.hpp" + +using namespace Dune; + +BOOST_AUTO_TEST_CASE(construction_and_queries) +{ + const SparseVector sv1; + BOOST_CHECK(sv1.empty()); + BOOST_CHECK_EQUAL(sv1.size(), 0); + BOOST_CHECK_EQUAL(sv1.nonzeroSize(), 0); + + const int size = 100; + const int num_elem = 9; + const int elem[num_elem] = { 9, 8, 7, 6, 5, 4, 3, 2, 1 }; + const int indices[num_elem] = { 1, 2, 3, 5, 8, 13, 21, 34, 55 }; + const SparseVector sv2(size, elem, elem + num_elem, indices, indices + num_elem); + BOOST_CHECK(!sv2.empty()); + BOOST_CHECK_EQUAL(sv2.size(), size); + BOOST_CHECK_EQUAL(sv2.element(0), 0); + BOOST_CHECK_EQUAL(sv2.element(1), 9); + BOOST_CHECK_EQUAL(sv2.element(2), 8); + BOOST_CHECK_EQUAL(sv2.element(3), 7); + BOOST_CHECK_EQUAL(sv2.element(4), 0); + BOOST_CHECK_EQUAL(sv2.element(5), 6); + BOOST_CHECK_EQUAL(sv2.element(55), 1); + BOOST_CHECK_EQUAL(sv2.element(99), 0); + BOOST_CHECK_EQUAL(sv2.nonzeroSize(), num_elem); + for (int i = 0; i < num_elem; ++i) { + BOOST_CHECK_EQUAL(sv2.nonzeroElement(i), elem[i]); + } + const SparseVector sv2_again(size, elem, elem + num_elem, indices, indices + num_elem); + BOOST_CHECK(sv2 == sv2_again); + SparseVector sv2_append(size, elem, elem + num_elem - 1, indices, indices + num_elem - 1); + BOOST_CHECK_EQUAL(sv2_append.nonzeroSize(), num_elem - 1); + sv2_append.addElement(elem[num_elem - 1], indices[num_elem - 1]); + BOOST_CHECK(sv2 == sv2_append); + SparseVector sv2_append2(size); + for (int i = 0; i < num_elem; ++i) { + sv2_append2.addElement(elem[i], indices[i]); + } + BOOST_CHECK(sv2 == sv2_append2); + sv2_append2.clear(); + SparseVector sv_empty; + BOOST_CHECK(sv2_append2 == sv_empty); + + // Tests that only run in debug mode. +#ifndef NDEBUG + // One element too few. + BOOST_CHECK_THROW(const SparseVector sv3(size, elem, elem + num_elem - 1, indices, indices + num_elem), std::exception); + // One element too many. + BOOST_CHECK_THROW(const SparseVector sv4(size, elem, elem + num_elem, indices, indices + num_elem - 1), std::exception); + // Indices out of range. + BOOST_CHECK_THROW(const SparseVector sv5(4, elem, elem + num_elem, indices, indices + num_elem), std::exception); + // Indices not strictly increasing. Cheating by using the elements as indices. + BOOST_CHECK_THROW(const SparseVector sv5(size, elem, elem + num_elem, elem, elem + num_elem), std::exception); + + // Do not ask for out of range indices. + BOOST_CHECK_THROW(sv1.element(0), std::exception); + BOOST_CHECK_THROW(sv2.element(-1), std::exception); + BOOST_CHECK_THROW(sv2.element(sv2.size()), std::exception); + BOOST_CHECK_THROW(sv2.nonzeroElement(sv2.nonzeroSize()), std::exception); +#endif +} + diff --git a/dune/common/test/unit_test.cpp b/dune/common/test/unit_test.cpp new file mode 100644 index 00000000..5904236b --- /dev/null +++ b/dune/common/test/unit_test.cpp @@ -0,0 +1,78 @@ +//=========================================================================== +// +// File: unit_test.cpp +// +// Created: Fri Jul 17 11:02:33 2009 +// +// Author(s): Bård Skaflestad +// Atgeirr F Rasmussen +// +// $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 . +*/ + +#include +#include +#include +#include + +#include + +#include "../Units.hpp" + +using namespace Dune::prefix; +using namespace Dune::unit; + +int main() +{ + std::cout << "m = " << meter << '\n'; + std::cout << "kg = " << kilogram << '\n'; + std::cout << "s = " << second << '\n'; + + std::cout << "mD = " << milli*darcy << '\n'; + std::cout << "MD = " << mega*darcy << '\n'; + + std::cout << "MD = " << convert::to(mega*darcy, milli*darcy) << " mD\n"; + + std::cout << "1 bar = " << convert::to(convert::from(1.0, barsa), psia) << " psi\n"; + + std::cout << "1 atm = " << convert::to(1*atm, barsa) << " bar\n"; + + std::vector flux(10, 10000*cubic(meter)/year); + + std::cout << "flux = ["; + std::copy(flux.begin(), flux.end(), + std::ostream_iterator(std::cout, " ")); + std::cout << "\b] (m^3/s)\n"; + + std::transform(flux.begin(), flux.end(), flux.begin(), + boost::bind(convert::to, _1, cubic(meter)/year)); + std::cout << " = ["; + std::copy(flux.begin(), flux.end(), + std::ostream_iterator(std::cout, " ")); + std::cout << "\b] (m^3/year)\n"; + + return 0; +}