OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

changed: remove some generic utility classes

Browse Source 
moved to opm-common
This commit is contained in:
Arne Morten Kvarving 2018-01-11 16:58:43 +01:00
parent 6c9365ac17
commit 9603fbd57f
8 changed files with 0 additions and 2267 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
CMakeLists_files.cmake
View File

@ -74,7 +74,6 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/transport/reorder/TransportSolverTwophaseReorder.cpp opm/core/transport/reorder/TransportSolverTwophaseReorder.cpp
opm/core/transport/reorder/reordersequence.cpp opm/core/transport/reorder/reordersequence.cpp
opm/core/transport/reorder/tarjan.c opm/core/transport/reorder/tarjan.c
opm/core/utility/MonotCubicInterpolator.cpp
opm/core/utility/VelocityInterpolation.cpp opm/core/utility/VelocityInterpolation.cpp
opm/core/utility/WachspressCoord.cpp opm/core/utility/WachspressCoord.cpp
opm/core/utility/compressedToCartesian.cpp opm/core/utility/compressedToCartesian.cpp
@ -97,7 +96,6 @@ list (APPEND TEST_SOURCE_FILES
tests/test_dgbasis.cpp tests/test_dgbasis.cpp
tests/test_cubic.cpp tests/test_cubic.cpp
tests/test_flowdiagnostics.cpp tests/test_flowdiagnostics.cpp
tests/test_nonuniformtablelinear.cpp
tests/test_parallelistlinformation.cpp tests/test_parallelistlinformation.cpp
tests/test_sparsevector.cpp tests/test_sparsevector.cpp
tests/test_velocityinterpolation.cpp tests/test_velocityinterpolation.cpp
@ -255,18 +253,13 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/transport/reorder/tarjan.h opm/core/transport/reorder/tarjan.h
opm/core/utility/CompressedPropertyAccess.hpp opm/core/utility/CompressedPropertyAccess.hpp
opm/core/utility/initHydroCarbonState.hpp opm/core/utility/initHydroCarbonState.hpp
opm/core/utility/MonotCubicInterpolator.hpp
opm/core/utility/NonuniformTableLinear.hpp
opm/core/utility/RegionMapping.hpp opm/core/utility/RegionMapping.hpp
opm/core/utility/RootFinders.hpp
opm/core/utility/SparseVector.hpp
opm/core/utility/UniformTableLinear.hpp opm/core/utility/UniformTableLinear.hpp
opm/core/utility/VelocityInterpolation.hpp opm/core/utility/VelocityInterpolation.hpp
opm/core/utility/WachspressCoord.hpp opm/core/utility/WachspressCoord.hpp
opm/core/utility/buildUniformMonotoneTable.hpp opm/core/utility/buildUniformMonotoneTable.hpp
opm/core/utility/compressedToCartesian.hpp opm/core/utility/compressedToCartesian.hpp
opm/core/utility/extractPvtTableIndex.hpp opm/core/utility/extractPvtTableIndex.hpp
opm/core/utility/linearInterpolation.hpp
opm/core/utility/miscUtilities.hpp opm/core/utility/miscUtilities.hpp
opm/core/utility/miscUtilitiesBlackoil.hpp opm/core/utility/miscUtilitiesBlackoil.hpp
opm/core/utility/miscUtilities_impl.hpp opm/core/utility/miscUtilities_impl.hpp

733
opm/core/utility/MonotCubicInterpolator.cpp
View File

@ -1,733 +0,0 @@
/*
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 "config.h"
#include <opm/core/utility/MonotCubicInterpolator.hpp>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <map>
#include <cmath>
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<double,double>
- 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<double, double> 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<double, <double, double> >
- 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.
*/
namespace Opm
{
MonotCubicInterpolator::
MonotCubicInterpolator(const vector<double> & x, const vector<double> & 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<double>::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<double> 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) {
if (std::isnan(newx) || std::isinf(newx) || std::isnan(newf) || std::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 {
if (std::isnan(x) || std::isinf(x)) {
throw("MonotCubicInterpolator: evaluate() received inf/nan input.");
}
// xf becomes the first (xdata,fdata) pair where xdata >= x
map<double,double>::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<double,double> xf2 = *xf_iterator;
pair<double,double> 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<double>
MonotCubicInterpolator::
get_xVector() const
{
map<double,double>::const_iterator xf_iterator = data.begin();
vector<double> 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<double>
MonotCubicInterpolator::
get_fVector() const
{
map<double,double>::const_iterator xf_iterator = data.begin();
vector<double> 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<double,double>::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<double,double>::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<double,double>
MonotCubicInterpolator::
getMissingX() const
{
if( data.size() < 2) {
throw("MonotCubicInterpolator::getMissingX() only one datapoint.");
}
// Search for biggest difference value in function-datavalues:
map<double,double>::const_iterator maxfDiffPair_iterator = data.begin();;
double maxfDiffValue = 0;
map<double,double>::const_iterator xf_iterator;
map<double,double>::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<double,double>
MonotCubicInterpolator::
getMaximumF() const {
if (data.size() <= 1) {
throw ("MonotCubicInterpolator::getMaximumF() empty data.") ;
}
if (strictlyIncreasing)
return *data.rbegin();
else if (strictlyDecreasing)
return *data.begin();
else {
pair<double,double> maxf = *data.rbegin() ;
map<double,double>::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<double,double>
MonotCubicInterpolator::
getMinimumF() const {
if (data.size() <= 1) {
throw ("MonotCubicInterpolator::getMinimumF() empty data.") ;
}
if (strictlyIncreasing)
return *data.begin();
else if (strictlyDecreasing) {
return *data.rbegin();
}
else {
pair<double,double> minf = *data.rbegin() ;
map<double,double>::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<double,double>::const_iterator xf_iterator;
map<double,double>::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<double,double>::iterator xf_iterator;
map<double,double>::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<double,double>::iterator xf_iterator;
map<double,double>::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 <double,double>::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<double,double>::const_iterator xf_prev_iterator;
map<double,double>::const_iterator xf_iterator;
map<double,double>::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<double,double>::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<double,double>::const_iterator lastpoint = intpoint; --lastpoint;
map<double,double>::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<double,double>::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<double,double>::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<double,double>::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 ;
}
}
}
} // namespace Opm

583
opm/core/utility/MonotCubicInterpolator.hpp
View File

@ -1,583 +0,0 @@
/* -*-C++-*- */
#ifndef _MONOTCUBICINTERPOLATOR_H
#define _MONOTCUBICINTERPOLATOR_H
#include <vector>
#include <map>
#include <string>
/*
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.
*/
namespace Opm
{
/**
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 <havb (at) statoil.com>, 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 \<whitespace\>\<float\>\<whitespace\>\<float\>\<whatever\>\<newline\>
*/
MonotCubicInterpolator(const std::string & datafilename)
{
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 \<whitespace\>\<float\>\<whitespace\>\<float\>\<whatever\>\<newline\>
All commas in the file will be treated as spaces when parsing.
*/
MonotCubicInterpolator(const char* datafilename)
{
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)
{
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)
{
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<double> & x ,
const std::vector<double> & 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 \<whitespace\>\<float\>\<whitespace\>\<float\>\<whatever\>\<newline\>
*/
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;
/**
@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<double,double> 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<double,double> 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<double,double> getMaximumF() const ;
/**
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<double,double> getMinimumF() const ;
/**
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<double> 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<double> 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);
/**
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<double,double> getMissingX() const;
/**
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<double, double> data;
// Data structure to store x- and d-values
mutable std::map<double, double> ddata;
// Storage containers for precomputed interpolation data
// std::vector<double> 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<double,double> 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;
}
}
};
} // namespace Opm
#endif

208
opm/core/utility/NonuniformTableLinear.hpp
View File

@ -1,208 +0,0 @@
/*
Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010, 2011, 2012 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED
#define OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED
#include <cmath>
#include <exception>
#include <vector>
#include <utility>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
namespace Opm
{
/// @brief This class uses linear interpolation to compute the value
/// (and its derivative) of a function f sampled at possibly
/// nonuniform points. If values outside the domain are sought,
/// values will be extrapolated linearly.
/// @tparam T the range type of the function (should be an algebraic ring type)
template<typename T>
class NonuniformTableLinear
{
public:
/// @brief Default constructor.
NonuniformTableLinear();
/// @brief Construct from vectors of x and y values.
/// @param x_values vector of domain values
/// @param y_values vector of corresponding range values.
template<class XContainer, class YContainer>
NonuniformTableLinear(const XContainer& x_values,
const YContainer& y_values);
/// @brief Get the domain.
/// @return the domain as a pair of doubles.
std::pair<double, double> domain();
/// @brief Rescale the domain.
/// @param new_domain the new domain as a pair of doubles.
void rescaleDomain(std::pair<double, double> new_domain);
/// @brief Evaluate the value at x.
/// @param x a domain value
/// @return f(x)
double operator()(const double x) const;
/// @brief Evaluate the derivative at x.
/// @param x a domain value
/// @return f'(x)
double derivative(const double x) const;
/// @brief Evaluate the inverse at y. Requires T to be a double.
/// @param y a range value
/// @return f^{-1}(y)
double inverse(const double y) const;
/// @brief Equality operator.
/// @param other another NonuniformTableLinear.
/// @return true if they are represented exactly alike.
bool operator==(const NonuniformTableLinear& other) const;
protected:
std::vector<double> x_values_;
std::vector<T> y_values_;
mutable std::vector<T> x_values_reversed_;
mutable std::vector<T> y_values_reversed_;
};
// A utility function
/// @brief Detect if a sequence is nondecreasing.
/// @tparam FI a forward iterator whose value type has operator< defined.
/// @param beg start of sequence
/// @param end one-beyond-end of sequence
/// @return false if there exists two consecutive values (v1, v2) in the sequence
/// for which v2 < v1, else returns true.
template <typename FI>
bool isNondecreasing(const FI beg, const FI end)
{
if (beg == end) return true;
FI it = beg;
++it;
FI prev = beg;
for (; it != end; ++it, ++prev) {
if (*it < *prev) {
return false;
}
}
return true;
}
// Member implementations.
template<typename T>
inline
NonuniformTableLinear<T>
::NonuniformTableLinear()
{
}
template<typename T>
template<class XContainer, class YContainer>
inline
NonuniformTableLinear<T>
::NonuniformTableLinear(const XContainer& x_column,
const YContainer& y_column)
: x_values_( x_column.begin() , x_column.end()),
y_values_( y_column.begin() , y_column.end())
{
assert(isNondecreasing(x_values_.begin(), x_values_.end()));
}
template<typename T>
inline std::pair<double, double>
NonuniformTableLinear<T>
::domain()
{
return std::make_pair(x_values_[0], x_values_.back());
}
template<typename T>
inline void
NonuniformTableLinear<T>
::rescaleDomain(std::pair<double, double> new_domain)
{
const double a = x_values_[0];
const double b = x_values_.back();
const double c = new_domain.first;
const double d = new_domain.second;
// x in [a, b] -> x in [c, d]
for (int i = 0; i < int(x_values_.size()); ++i) {
x_values_[i] = (x_values_[i] - a)*(d - c)/(b - a) + c;
}
}
template<typename T>
inline double
NonuniformTableLinear<T>
::operator()(const double x) const
{
return Opm::linearInterpolation(x_values_, y_values_, x);
}
template<typename T>
inline double
NonuniformTableLinear<T>
::derivative(const double x) const
{
return Opm::linearInterpolationDerivative(x_values_, y_values_, x);
}
template<typename T>
inline double
NonuniformTableLinear<T>
::inverse(const double y) const
{
if (y_values_.front() < y_values_.back()) {
return Opm::linearInterpolation(y_values_, x_values_, y);
} else {
if (y_values_reversed_.empty()) {
y_values_reversed_ = y_values_;
std::reverse(y_values_reversed_.begin(), y_values_reversed_.end());
assert(isNondecreasing(y_values_reversed_.begin(), y_values_reversed_.end()));
x_values_reversed_ = x_values_;
std::reverse(x_values_reversed_.begin(), x_values_reversed_.end());
}
return Opm::linearInterpolation(y_values_reversed_, x_values_reversed_, y);
}
}
template<typename T>
inline bool
NonuniformTableLinear<T>
::operator==(const NonuniformTableLinear<T>& other) const
{
return x_values_ == other.x_values_
&& y_values_ == other.y_values_;
}
} // namespace Opm
#endif // OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED

340
opm/core/utility/RootFinders.hpp
View File

@ -1,340 +0,0 @@
//===========================================================================
//
// File: RootFinders.hpp
//
// Created: Thu May 6 19:59:42 2010
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Jostein R Natvig <jostein.r.natvig@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
Copyright 2010 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ROOTFINDERS_HEADER
#define OPM_ROOTFINDERS_HEADER
#include <opm/common/ErrorMacros.hpp>
#include <algorithm>
#include <limits>
#include <cmath>
#include <iostream>
namespace Opm
{
struct ThrowOnError
{
static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1)
{
OPM_THROW(std::runtime_error, "Error in parameters, zero not bracketed: [a, b] = ["
<< x0 << ", " << x1 << "] f(a) = " << f0 << " f(b) = " << f1);
return -1e100; // Never reached.
}
static double handleTooManyIterations(const double x0, const double x1, const int maxiter)
{
OPM_THROW(std::runtime_error, "Maximum number of iterations exceeded: " << maxiter << "\n"
<< "Current interval is [" << std::min(x0, x1) << ", "
<< std::max(x0, x1) << "]");
return -1e100; // Never reached.
}
};
struct WarnAndContinueOnError
{
static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1)
{
OPM_REPORT;
std::cerr << "Error in parameters, zero not bracketed: [a, b] = ["
<< x0 << ", " << x1 << "] f(a) = " << f0 << " f(b) = " << f1
<< "";
return std::fabs(f0) < std::fabs(f1) ? x0 : x1;
}
static double handleTooManyIterations(const double x0, const double x1, const int maxiter)
{
OPM_REPORT;
std::cerr << "Maximum number of iterations exceeded: " << maxiter
<< ", current interval is [" << std::min(x0, x1) << ", "
<< std::max(x0, x1) << "]";
return 0.5*(x0 + x1);
}
};
struct ContinueOnError
{
static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1)
{
return std::fabs(f0) < std::fabs(f1) ? x0 : x1;
}
static double handleTooManyIterations(const double x0, const double x1, const int /*maxiter*/)
{
return 0.5*(x0 + x1);
}
};
template <class ErrorPolicy = ThrowOnError>
class RegulaFalsi
{
public:
/// 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 <class Functor>
inline static double solve(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<double>::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) {
return ErrorPolicy::handleBracketingFailure(a, b, f0, 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) >= 1e-9*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) {
return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
}
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);
}
/// 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.
/// This version takes an extra parameter for the initial guess.
template <class Functor>
inline static double solve(const Functor& f,
const double initial_guess,
const double a,
const double b,
const int max_iter,
const double tolerance,
int& iterations_used)
{
using namespace std;
const double macheps = numeric_limits<double>::epsilon();
const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
double f_initial = f(initial_guess);
const double epsF = tolerance + macheps*max(fabs(f_initial), 1.0);
if (fabs(f_initial) < epsF) {
return initial_guess;
}
double x0 = a;
double x1 = b;
double f0 = f_initial;
double f1 = f_initial;
if (x0 != initial_guess) {
f0 = f(x0);
if (fabs(f0) < epsF) {
return x0;
}
}
if (x1 != initial_guess) {
f1 = f(x1);
if (fabs(f1) < epsF) {
return x1;
}
}
if (f0*f_initial < 0.0) {
x1 = initial_guess;
f1 = f_initial;
} else {
x0 = initial_guess;
f0 = f_initial;
}
if (f0*f1 > 0.0) {
return ErrorPolicy::handleBracketingFailure(a, b, f0, 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) >= 1e-9*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) {
return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
}
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);
}
private:
inline static 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);
}
};
/// Attempts to find an interval bracketing a zero by successive
/// enlargement of search interval.
template <class Functor>
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) {
OPM_THROW(std::runtime_error, "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 Opm
#endif // OPM_ROOTFINDERS_HEADER

198
opm/core/utility/SparseVector.hpp
View File

@ -1,198 +0,0 @@
//===========================================================================
//
// File: SparseVector.hpp
//
// Created: Mon Jun 29 15:28:59 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <bard.skaflestad@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SPARSEVECTOR_HEADER
#define OPM_SPARSEVECTOR_HEADER
#include <vector>
#include <numeric>
#include <algorithm>
#include <opm/common/ErrorMacros.hpp>
namespace Opm
{
/// 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 <typename T>
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 sz)
: size_(sz), 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 <typename DataIter, typename IntegerIter>
SparseVector(int sz,
DataIter data_beg, DataIter data_end,
IntegerIter index_beg, IntegerIter index_end)
: size_(sz), data_(data_beg, data_end), indices_(index_beg, index_end),
default_elem_()
{
#ifndef NDEBUG
OPM_ERROR_IF(sz < 0, "The size of a SparseVector must be non-negative");
OPM_ERROR_IF(indices_.size() != data_.size(), "The number of indices of a SparseVector must equal to the number of entries");
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 >= sz) {
OPM_THROW(std::logic_error, "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& elem, int index)
{
assert(indices_.empty() || index > indices_.back());
assert(index < size_);
data_.push_back(elem);
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
{
#ifndef NDEBUG
OPM_ERROR_IF(index < 0, "The index of a SparseVector must be non-negative (is " << index << ")");
OPM_ERROR_IF(index >= size_, "The index of a SparseVector must be smaller than the maximum value (is " << index << ", max value: " << size_ <<")");
#endif
std::vector<int>::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
{
#ifndef NDEBUG
OPM_ERROR_IF(nzindex < 0, "The index of a SparseVector must be non-negative (is " << nzindex << ")");
OPM_ERROR_IF(nzindex >= nonzeroSize(), "The index of a SparseVector must be smaller than the maximum value (is " << nzindex << ", max value: " << nonzeroSize() <<")");
#endif
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<T> data_;
std::vector<int> indices_;
T default_elem_;
};
} // namespace Opm
#endif // OPM_SPARSEVECTOR_HEADER

107
opm/core/utility/linearInterpolation.hpp
View File

@ -1,107 +0,0 @@
/*
Copyright 2009, 2010, 2013 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_LINEARINTERPOLATION_HEADER_INCLUDED
#define OPM_LINEARINTERPOLATION_HEADER_INCLUDED
#include <vector>
#include <algorithm>
namespace Opm
{
inline int tableIndex(const std::vector<double>& 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 linearInterpolationDerivative(const std::vector<double>& xv,
const std::vector<double>& 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]);
}
inline double linearInterpolation(const std::vector<double>& xv,
const std::vector<double>& 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 linearInterpolationNoExtrapolation(const std::vector<double>& xv,
const std::vector<double>& yv, double x)
{
// Return end values if x is outside xv
if (x < xv.front()) {
return yv.front();
}
if (x > xv.back()) {
return yv.back();
}
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 linearInterpolation(const std::vector<double>& xv,
const std::vector<double>& 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 Opm
#endif // OPM_LINEARINTERPOLATION_HEADER_INCLUDED

91
tests/test_nonuniformtablelinear.cpp
View File

@ -1,91 +0,0 @@
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
OpenRS is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenRS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#if defined(HAVE_DYNAMIC_BOOST_TEST)
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // to suppress our messages when throwing
#define BOOST_TEST_MODULE NonuniformTableLinearTests
#include <boost/test/unit_test.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
BOOST_AUTO_TEST_CASE(utility_functions)
{
// Test isNondecreasing().
using Opm::isNondecreasing;
double xva1[] = { -1.0, 2.0, 2.2, 3.0, 5.0 };
const int numvals1 = sizeof(xva1)/sizeof(xva1[0]);
BOOST_CHECK(isNondecreasing(xva1, xva1 + numvals1));
double xva2[] = { -1.0, 2.0, 2.0, 2.0, 5.0 };
const int numvals2 = sizeof(xva2)/sizeof(xva2[0]);
BOOST_CHECK(isNondecreasing(xva2, xva2 + numvals2));
double xva3[] = { -1.0, 2.0, 1.9, 3.0, 5.0 };
const int numvals3 = sizeof(xva3)/sizeof(xva3[0]);
BOOST_CHECK(!isNondecreasing(xva3, xva3 + numvals3));
}
BOOST_AUTO_TEST_CASE(table_operations)
{
// Make a simple table.
double xva[] = { -1.0, 2.0, 2.2, 3.0, 5.0 };
const int numvals = sizeof(xva)/sizeof(xva[0]);
std::vector<double> xv(xva, xva + numvals);
double yva[numvals] = { 1.0, 2.0, 3.0, 4.0, 2.0 };
std::vector<double> yv(yva, yva + numvals);
Opm::NonuniformTableLinear<double> t1(xv, yv);
Opm::NonuniformTableLinear<double> t1_copy1(xv, yv);
Opm::NonuniformTableLinear<double> t1_copy2(t1);
// Check equality.
BOOST_CHECK(t1 == t1_copy1);
BOOST_CHECK(t1 == t1_copy2);
// Check some evaluations.
for (int i = 0; i < numvals; ++i) {
BOOST_CHECK_EQUAL(t1(xv[i]), yv[i]);
}
BOOST_CHECK_EQUAL(t1(2.6), 3.5);
BOOST_CHECK_EQUAL(t1(4.0), 3.0);
BOOST_CHECK_EQUAL(t1.derivative(4.0), -1.0);
// Derivatives at endpoints.
BOOST_CHECK_CLOSE(t1.derivative(-1.0), 1.0/3.0, 1e-13);
BOOST_CHECK_CLOSE(t1.derivative(5.0), -1.0, 1e-13);
// Extrapolation of values.
BOOST_CHECK_CLOSE(t1(xv[0] - 1.0), 2.0/3.0, 1e-13);
BOOST_CHECK_CLOSE(t1(xv.back() + 1.0), 1.0, 1e-13);
// Domains.
BOOST_CHECK_EQUAL(t1.domain().first, xv[0]);
BOOST_CHECK_EQUAL(t1.domain().second, xv.back());
std::pair<double, double> new_domain(-100.0, 20.0);
t1.rescaleDomain(new_domain);
BOOST_CHECK_EQUAL(t1.domain().first, new_domain.first);
BOOST_CHECK_EQUAL(t1.domain().second, new_domain.second);
for (int i = 0; i < numvals; ++i) {
BOOST_CHECK_EQUAL(t1((xv[i] + 1.0)*20.0 - 100.0), yv[i]);
}
BOOST_CHECK_EQUAL(t1(0.0), 3.0);
BOOST_CHECK(std::fabs(t1.derivative(0.0) + 1.0/20.0) < 1e-11);
}