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

Merge branch 'master' into initialisation

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen
2014-02-10 13:24:35 +01:00
parent 4396130765 995064635f
commit b3ac9529c8
13 changed files with 158 additions and 69 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@@ -306,6 +306,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/props/rock/RockBasic.hpp opm/core/props/rock/RockBasic.hpp
opm/core/props/rock/RockCompressibility.hpp opm/core/props/rock/RockCompressibility.hpp
opm/core/props/rock/RockFromDeck.hpp opm/core/props/rock/RockFromDeck.hpp
opm/core/props/satfunc/SatFuncBase.hpp
opm/core/props/satfunc/SatFuncGwseg.hpp opm/core/props/satfunc/SatFuncGwseg.hpp
opm/core/props/satfunc/SatFuncSimple.hpp opm/core/props/satfunc/SatFuncSimple.hpp
opm/core/props/satfunc/SatFuncStone2.hpp opm/core/props/satfunc/SatFuncStone2.hpp

3
examples/compute_tof.cpp
View File

@@ -68,7 +68,6 @@ namespace
} }
void buildTracerheadsFromWells(const Wells* wells, void buildTracerheadsFromWells(const Wells* wells,
const Opm::WellState& well_state,
Opm::SparseTable<int>& tracerheads) Opm::SparseTable<int>& tracerheads)
{ {
if (wells == 0) { if (wells == 0) {
@@ -257,7 +256,7 @@ try
std::vector<double> tracer; std::vector<double> tracer;
Opm::SparseTable<int> tracerheads; Opm::SparseTable<int> tracerheads;
if (compute_tracer) { if (compute_tracer) {
buildTracerheadsFromWells(wells->c_wells(), well_state, tracerheads); buildTracerheadsFromWells(wells->c_wells(), tracerheads);
} }
if (use_dg) { if (use_dg) {
if (compute_tracer) { if (compute_tracer) {

109
opm/core/io/eclipse/EclipseWriter.cpp
View File

@@ -575,6 +575,18 @@ private:
fortio_fclose) { } fortio_fclose) { }
}; };
// in order to get RTTI for this "class" (which is just a typedef), we must
// ask the compiler to explicitly instantiate it.
template struct EclipseHandle<ecl_sum_tstep_struct>;
} // anonymous namespace
// Note: the following parts were taken out of the anonymous
// namespace, since EclipseSummary is now used as a pointer member in
// EclipseWriter and forward declared in EclipseWriter.hpp.
// forward decl. of mutually dependent type // forward decl. of mutually dependent type
struct EclipseWellReport; struct EclipseWellReport;
@@ -655,9 +667,6 @@ private:
} }
}; };
// in order to get RTTI for this "class" (which is just a typedef), we must
// ask the compiler to explicitly instantiate it.
template struct EclipseHandle<ecl_sum_tstep_struct>;
/** /**
* Summary variable that reports a characteristics of a well. * Summary variable that reports a characteristics of a well.
@@ -670,7 +679,7 @@ protected:
PhaseUsage uses, /* phases present */ PhaseUsage uses, /* phases present */
BlackoilPhases::PhaseIndex phase, /* oil, water or gas */ BlackoilPhases::PhaseIndex phase, /* oil, water or gas */
WellType type, /* prod. or inj. */ WellType type, /* prod. or inj. */
char aggregation, /* rate or total */ char aggregation, /* rate or total or BHP */
std::string unit) std::string unit)
: EclipseHandle <smspec_node_type> ( : EclipseHandle <smspec_node_type> (
ecl_sum_add_var (summary, ecl_sum_add_var (summary,
@@ -718,22 +727,26 @@ private:
char aggregation) { char aggregation) {
std::string name; std::string name;
name += 'W'; // well name += 'W'; // well
switch (phase) { if (aggregation == 'B') {
name += "BHP";
} else {
switch (phase) {
case BlackoilPhases::Aqua: name += 'W'; break; /* water */ case BlackoilPhases::Aqua: name += 'W'; break; /* water */
case BlackoilPhases::Vapour: name += 'G'; break; /* gas */ case BlackoilPhases::Vapour: name += 'G'; break; /* gas */
case BlackoilPhases::Liquid: name += 'O'; break; /* oil */ case BlackoilPhases::Liquid: name += 'O'; break; /* oil */
default: default:
OPM_THROW(std::runtime_error, OPM_THROW(std::runtime_error,
"Unknown phase used in blackoil reporting"); "Unknown phase used in blackoil reporting");
} }
switch (type) { switch (type) {
case WellType::INJECTOR: name += 'I'; break; case WellType::INJECTOR: name += 'I'; break;
case WellType::PRODUCER: name += 'P'; break; case WellType::PRODUCER: name += 'P'; break;
default: default:
OPM_THROW(std::runtime_error, OPM_THROW(std::runtime_error,
"Unknown well type used in blackoil reporting"); "Unknown well type used in blackoil reporting");
}
name += aggregation; /* rate ('R') or total ('T') */
} }
name += aggregation; /* rate ('R') or total ('T') */
return name; return name;
} }
protected: protected:
@@ -743,6 +756,13 @@ protected:
const double value = sign_ * wellState.wellRates () [index_] * convFactor; const double value = sign_ * wellState.wellRates () [index_] * convFactor;
return value; return value;
} }
double bhp (const WellState& wellstate) {
// Note that 'index_' is used here even though it is meant
// to give a (well,phase) pair.
const int num_phases = wellstate.wellRates().size() / wellstate.bhp().size();
return wellstate.bhp()[index_/num_phases];
}
}; };
/// Monitors the rate given by a well. /// Monitors the rate given by a well.
@@ -761,7 +781,7 @@ struct EclipseWellRate : public EclipseWellReport {
type, type,
'R', 'R',
"SM3/DAY" /* surf. cub. m. per day */ ) { } "SM3/DAY" /* surf. cub. m. per day */ ) { }
virtual double update (const SimulatorTimer& timer, virtual double update (const SimulatorTimer& /*timer*/,
const WellState& wellState) { const WellState& wellState) {
// TODO: Why only positive rates? // TODO: Why only positive rates?
return std::max (0., rate (wellState)); return std::max (0., rate (wellState));
@@ -790,6 +810,11 @@ struct EclipseWellTotal : public EclipseWellReport {
virtual double update (const SimulatorTimer& timer, virtual double update (const SimulatorTimer& timer,
const WellState& wellState) { const WellState& wellState) {
if (timer.currentStepNum() == 0) {
// We are at the initial state.
// No step has been taken yet.
return 0.0;
}
// TODO: Is the rate average for the timestep, or is in // TODO: Is the rate average for the timestep, or is in
// instantaneous (in which case trapezoidal or Simpson integration // instantaneous (in which case trapezoidal or Simpson integration
// would probably be better) // would probably be better)
@@ -805,6 +830,31 @@ private:
double total_; double total_;
}; };
/// Monitors the bottom hole pressure in a well.
struct EclipseWellBhp : public EclipseWellReport {
EclipseWellBhp (const EclipseSummary& summary,
const EclipseGridParser& parser,
int whichWell,
PhaseUsage uses,
BlackoilPhases::PhaseIndex phase,
WellType type)
: EclipseWellReport (summary,
parser,
whichWell,
uses,
phase,
type,
'B',
"Pascal")
{ }
virtual double update (const SimulatorTimer& /*timer*/,
const WellState& wellState)
{
return bhp(wellState);
}
};
inline void inline void
EclipseSummary::writeTimeStep (const SimulatorTimer& timer, EclipseSummary::writeTimeStep (const SimulatorTimer& timer,
const WellState& wellState) { const WellState& wellState) {
@@ -861,8 +911,30 @@ EclipseSummary::addWells (const EclipseGridParser& parser,
} }
} }
} }
// Add BHP monitors
for (int whichWell = 0; whichWell != numWells; ++whichWell) {
// In the call below: uses, phase and the well type arguments
// are not used, except to set up an index that stores the
// well indirectly. For details see the implementation of the
// EclipseWellReport constructor, and the method
// EclipseWellReport::bhp().
BlackoilPhases::PhaseIndex phase = BlackoilPhases::Liquid;
if (!uses.phase_used[BlackoilPhases::Liquid]) {
phase = BlackoilPhases::Vapour;
}
add (std::unique_ptr <EclipseWellReport> (
new EclipseWellBhp (*this,
parser,
whichWell,
uses,
phase,
WELL_TYPES[0])));
}
} }
namespace {
/// Helper method that can be used in keyword transformation (must curry /// Helper method that can be used in keyword transformation (must curry
/// the barsa argument) /// the barsa argument)
static double toBar (const double& pressure) { static double toBar (const double& pressure) {
@@ -902,11 +974,16 @@ void EclipseWriter::writeInit(const SimulatorTimer &timer,
/* Initial solution (pressure and saturation) */ /* Initial solution (pressure and saturation) */
writeSolution (timer, reservoirState, wellState); writeSolution (timer, reservoirState, wellState);
/* Create summary object (could not do it at construction time,
since it requires knowledge of the start time). */
summary_.reset(new EclipseSummary(outputDir_, baseName_, timer, *parser_));
summary_->addWells (*parser_, uses_);
} }
void EclipseWriter::writeSolution (const SimulatorTimer& timer, void EclipseWriter::writeSolution (const SimulatorTimer& timer,
const SimulatorState& reservoirState, const SimulatorState& reservoirState,
const WellState& wellState) { const WellState& /*wellState*/) {
// start writing to files // start writing to files
EclipseRestart rst (outputDir_, EclipseRestart rst (outputDir_,
baseName_, baseName_,
@@ -952,9 +1029,15 @@ void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
// (first timestep, in practice), and reused later. but how to do this // (first timestep, in practice), and reused later. but how to do this
// without keeping the complete summary in memory (which will then // without keeping the complete summary in memory (which will then
// accumulate all the timesteps)? // accumulate all the timesteps)?
EclipseSummary sum (outputDir_, baseName_, timer, *parser_); //
sum.addWells (*parser_, uses_); // Note: The answer to the question above is still not settled,
sum.writeTimeStep (timer, wellState); // but now we do keep the complete summary in memory, as a member
// variable in the EclipseWriter class, instead of creating a
// temporary EclipseSummary in this function every time it is
// called. This has been changed so that the final summary file
// will contain data from the whole simulation, instead of just
// the last step.
summary_->writeTimeStep (timer, wellState);
} }
#else #else

2
opm/core/io/eclipse/EclipseWriter.hpp
View File

@@ -28,6 +28,7 @@
#include <memory> // std::unique_ptr #include <memory> // std::unique_ptr
struct UnstructuredGrid; struct UnstructuredGrid;
struct EclipseSummary;
namespace Opm { namespace Opm {
@@ -90,6 +91,7 @@ private:
std::string outputDir_; std::string outputDir_;
std::string baseName_; std::string baseName_;
PhaseUsage uses_; // active phases in the input deck PhaseUsage uses_; // active phases in the input deck
std::shared_ptr <EclipseSummary> summary_;
/// Write solution field variables (pressure and saturation) /// Write solution field variables (pressure and saturation)
void writeSolution (const SimulatorTimer& timer, void writeSolution (const SimulatorTimer& timer,

2
opm/core/linalg/LinearSolverUmfpack.cpp
View File

@@ -56,7 +56,7 @@ namespace Opm
const_cast<double*>(sa) const_cast<double*>(sa)
}; };
call_UMFPACK(&A, rhs, solution); call_UMFPACK(&A, rhs, solution);
LinearSolverReport rep = {0}; LinearSolverReport rep = {};
rep.converged = true; rep.converged = true;
return rep; return rep;
} }

24
opm/core/props/pvt/SinglePvtLiveGas.cpp
View File

@@ -99,12 +99,12 @@ namespace Opm
} }
/// Viscosity and its derivatives as a function of p and r. /// Viscosity and its derivatives as a function of p and r.
void SinglePvtLiveGas::mu(const int n, void SinglePvtLiveGas::mu(const int /*n*/,
const double* p, const double* /*p*/,
const double* r, const double* /*r*/,
double* output_mu, double* /*output_mu*/,
double* output_dmudp, double* /*output_dmudp*/,
double* output_dmudr) const double* /*output_dmudr*/) const
{ {
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented"); OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
} }
@@ -156,12 +156,12 @@ namespace Opm
} }
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r. /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
void SinglePvtLiveGas::b(const int n, void SinglePvtLiveGas::b(const int /*n*/,
const double* p, const double* /*p*/,
const double* r, const double* /*r*/,
double* output_b, double* /*output_b*/,
double* output_dbdp, double* /*output_dbdp*/,
double* output_dbdr) const double* /*output_dbdr*/) const
{ {
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented"); OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");

4
opm/core/props/satfunc/SatFuncSimple.cpp
View File

@@ -36,7 +36,7 @@ namespace Opm
void SatFuncBase<NonuniformTableLinear<double> >::initializeTableType(NonuniformTableLinear<double> & table, void SatFuncBase<NonuniformTableLinear<double> >::initializeTableType(NonuniformTableLinear<double> & table,
const std::vector<double>& arg, const std::vector<double>& arg,
const std::vector<double>& value, const std::vector<double>& value,
const int samples) const int /*samples*/)
{ {
table = NonuniformTableLinear<double>(arg, value); table = NonuniformTableLinear<double>(arg, value);
} }
@@ -84,7 +84,7 @@ namespace Opm
} }
} }
double EPSTransforms::Transform::scaleSatDeriv(double s, double s_r, double s_cr, double s_max) const double EPSTransforms::Transform::scaleSatDeriv(double s, double /*s_r*/, double /*s_cr*/, double /*s_max*/) const
{ {
if (doNotScale) { if (doNotScale) {
return 1.0; return 1.0;

10
opm/core/props/satfunc/SatFuncSimple.hpp
View File

@@ -32,16 +32,16 @@ namespace Opm
void evalPc(const double* s, double* pc) const; void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const; void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
void evalKr(const double* s, double* kr, const EPSTransforms* epst) const void evalKr(const double* /*s*/, double* /*kr*/, const EPSTransforms* /*epst*/) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
void evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const void evalKr(const double* /*s*/, double* /*kr*/, const EPSTransforms* /*epst*/, const EPSTransforms* /*epst_hyst*/, const SatHyst* /*sat_hyst*/) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const; void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const void evalKrDeriv(const double* /*s*/, double* /*kr*/, double* /*dkrds*/, const EPSTransforms* /*epst*/, const EPSTransforms* /*epst_hyst*/, const SatHyst* /*sat_hyst*/) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
void evalPc(const double* s, double* pc, const EPSTransforms* epst) const void evalPc(const double* /*s*/, double* /*pc*/, const EPSTransforms* /*epst*/) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
void evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const void evalPcDeriv(const double* /*s*/, double* /*pc*/, double* /*dpcds*/, const EPSTransforms* /*epst*/) const
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
private: private:

12
opm/core/props/satfunc/SatFuncStone2.hpp
View File

@@ -32,17 +32,17 @@ namespace Opm
void evalPc(const double* s, double* pc) const; void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const; void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
void evalKr(const double* s, double* kr, const EPSTransforms* epst) const void evalKr(const double* /*s*/, double* /*kr*/, const EPSTransforms* /*epst*/) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
void evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const void evalKr(const double* /*s*/, double* /*kr*/, const EPSTransforms* /*epst*/, const EPSTransforms* /*epst_hyst*/, const SatHyst* /*sat_hyst*/) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const void evalKrDeriv(const double* /*s*/, double* /*kr*/, double* /*dkrds*/, const EPSTransforms* /*epst*/) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const void evalKrDeriv(const double* /*s*/, double* /*kr*/, double* /*dkrds*/, const EPSTransforms* /*epst*/, const EPSTransforms* /*epst_hyst*/, const SatHyst* /*sat_hyst*/) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
void evalPc(const double* s, double* pc, const EPSTransforms* epst) const void evalPc(const double* /*s*/, double* /*pc*/, const EPSTransforms* /*epst*/) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
void evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const void evalPcDeriv(const double* /*s*/, double* /*pc*/, double* /*dpcds*/, const EPSTransforms* /*epst*/) const
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}

5
opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
View File

@@ -100,8 +100,9 @@ namespace Opm
} }
// Saturation table scaling // Saturation table scaling
do_eps_ = false; do_hyst_ = false;
do_3pt_ = false; do_eps_ = false;
do_3pt_ = false;
if (deck.hasField("ENDSCALE")) { if (deck.hasField("ENDSCALE")) {
//if (!phase_usage_.phase_used[Aqua] || !phase_usage_.phase_used[Liquid] || phase_usage_.phase_used[Vapour]) { //if (!phase_usage_.phase_used[Aqua] || !phase_usage_.phase_used[Liquid] || phase_usage_.phase_used[Vapour]) {
// OPM_THROW(std::runtime_error, "Currently endpoint-scaling limited to oil-water systems without gas."); // OPM_THROW(std::runtime_error, "Currently endpoint-scaling limited to oil-water systems without gas.");

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

@@ -68,16 +68,18 @@ namespace Opm
{ {
static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1) static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1)
{ {
OPM_MESSAGE("Error in parameters, zero not bracketed: [a, b] = [" OPM_REPORT;
<< x0 << ", " << x1 << "] f(a) = " << f0 << " f(b) = " << f1 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; return std::fabs(f0) < std::fabs(f1) ? x0 : x1;
} }
static double handleTooManyIterations(const double x0, const double x1, const int maxiter) static double handleTooManyIterations(const double x0, const double x1, const int maxiter)
{ {
OPM_MESSAGE("Maximum number of iterations exceeded: " << maxiter OPM_REPORT;
<< ", current interval is [" << std::min(x0, x1) << ", " std::cerr << "Maximum number of iterations exceeded: " << maxiter
<< std::max(x0, x1) << "]"); << ", current interval is [" << std::min(x0, x1) << ", "
<< std::max(x0, x1) << "]";
return 0.5*(x0 + x1); return 0.5*(x0 + x1);
} }
}; };

9
opm/core/utility/Spline.hpp
View File

@@ -945,6 +945,7 @@ public:
int r = 3; int r = 3;
if (x0 < xMin()) { if (x0 < xMin()) {
static_cast<void>(extrapolate);
assert(extrapolate); assert(extrapolate);
Scalar m = evalDerivative_(xMin(), /*segmentIdx=*/0); Scalar m = evalDerivative_(xMin(), /*segmentIdx=*/0);
if (std::abs(m) < 1e-20) if (std::abs(m) < 1e-20)
@@ -1575,16 +1576,16 @@ protected:
{ return 2*3*t - 2; } { return 2*3*t - 2; }
// third derivative of the hermite basis functions // third derivative of the hermite basis functions
Scalar h00_prime3_(Scalar t) const Scalar h00_prime3_(Scalar /*t*/) const
{ return 2*3*2; } { return 2*3*2; }
Scalar h10_prime3_(Scalar t) const Scalar h10_prime3_(Scalar /*t*/) const
{ return 2*3; } { return 2*3; }
Scalar h01_prime3_(Scalar t) const Scalar h01_prime3_(Scalar /*t*/) const
{ return -2*3*2; } { return -2*3*2; }
Scalar h11_prime3_(Scalar t) const Scalar h11_prime3_(Scalar /*t*/) const
{ return 2*3; } { return 2*3; }
// returns the monotonicality of an interval of a spline segment // returns the monotonicality of an interval of a spline segment

22
tests/test_spline.cpp
View File

@@ -80,14 +80,14 @@ void testCommon(const Spline &sp,
// make sure the derivatives are consistent with the curve // make sure the derivatives are consistent with the curve
int np = 3*n; int np = 3*n;
for (int i = 0; i < np; ++i) { for (int i = 0; i < np; ++i) {
double x = sp.xMin() + (sp.xMax() - sp.xMin())*i/np; double xval = sp.xMin() + (sp.xMax() - sp.xMin())*i/np;
// first derivative // first derivative
double y1 = sp.eval(x+epsFD); double y1 = sp.eval(xval+epsFD);
double y0 = sp.eval(x); double y0 = sp.eval(xval);
double mFD = (y1 - y0)/epsFD; double mFD = (y1 - y0)/epsFD;
double m = sp.evalDerivative(x); double m = sp.evalDerivative(xval);
if (std::abs( mFD - m ) > 1000*epsFD) if (std::abs( mFD - m ) > 1000*epsFD)
OPM_THROW(std::runtime_error, OPM_THROW(std::runtime_error,
@@ -95,11 +95,11 @@ void testCommon(const Spline &sp,
" (" << mFD << " - " << m << " = " << mFD - m << ")!"); " (" << mFD << " - " << m << " = " << mFD - m << ")!");
// second derivative // second derivative
y1 = sp.evalDerivative(x+epsFD); y1 = sp.evalDerivative(xval+epsFD);
y0 = sp.evalDerivative(x); y0 = sp.evalDerivative(xval);
mFD = (y1 - y0)/epsFD; mFD = (y1 - y0)/epsFD;
m = sp.evalSecondDerivative(x); m = sp.evalSecondDerivative(xval);
if (std::abs( mFD - m ) > 1000*epsFD) if (std::abs( mFD - m ) > 1000*epsFD)
OPM_THROW(std::runtime_error, OPM_THROW(std::runtime_error,
@@ -107,11 +107,11 @@ void testCommon(const Spline &sp,
" (" << mFD << " - " << m << " = " << mFD - m << ")!"); " (" << mFD << " - " << m << " = " << mFD - m << ")!");
// Third derivative // Third derivative
y1 = sp.evalSecondDerivative(x+epsFD); y1 = sp.evalSecondDerivative(xval+epsFD);
y0 = sp.evalSecondDerivative(x); y0 = sp.evalSecondDerivative(xval);
mFD = (y1 - y0)/epsFD; mFD = (y1 - y0)/epsFD;
m = sp.evalThirdDerivative(x); m = sp.evalThirdDerivative(xval);
if (std::abs( mFD - m ) > 1000*epsFD) if (std::abs( mFD - m ) > 1000*epsFD)
OPM_THROW(std::runtime_error, OPM_THROW(std::runtime_error,
@@ -331,7 +331,7 @@ void plot()
std::cout << "\n"; std::cout << "\n";
} }
int main(int argc, char** argv) int main()
{ {
try { try {
testAll(); testAll();