Merge remote-tracking branch 'upstream/master' into opm-parser-integrate

This commit is contained in:
Joakim Hove
2014-02-05 23:06:54 +01:00
10 changed files with 70 additions and 61 deletions

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) {

View File

@@ -575,6 +575,12 @@ 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 } // anonymous namespace
// Note: the following parts were taken out of the anonymous // Note: the following parts were taken out of the anonymous
@@ -662,10 +668,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.
*/ */
@@ -779,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));
@@ -808,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)
@@ -841,7 +848,7 @@ struct EclipseWellBhp : public EclipseWellReport {
"Pascal") "Pascal")
{ } { }
virtual double update (const SimulatorTimer& timer, virtual double update (const SimulatorTimer& /*timer*/,
const WellState& wellState) const WellState& wellState)
{ {
return bhp(wellState); return bhp(wellState);
@@ -976,7 +983,7 @@ void EclipseWriter::writeInit(const SimulatorTimer &timer,
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_,

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;
} }

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");

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;

View File

@@ -31,23 +31,23 @@ namespace Opm
void evalKrDeriv(const double* s, double* kr, double* dkrds) const; void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
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:
}; };
typedef SatFuncSimple<UniformTableLinear<double> > SatFuncSimpleUniform; typedef SatFuncSimple<UniformTableLinear<double> > SatFuncSimpleUniform;
typedef SatFuncSimple<NonuniformTableLinear<double> > SatFuncSimpleNonuniform; typedef SatFuncSimple<NonuniformTableLinear<double> > SatFuncSimpleNonuniform;
@@ -185,7 +185,7 @@ namespace Opm
double dkrww = _dsdsw*epst->wat.scaleKrDeriv(s[wpos], this->krw_.derivative(_sw)); double dkrww = _dsdsw*epst->wat.scaleKrDeriv(s[wpos], this->krw_.derivative(_sw));
double krg = epst->gas.scaleKr(s[gpos], this->krg_(_sg), this->krgr_); double krg = epst->gas.scaleKr(s[gpos], this->krg_(_sg), this->krgr_);
double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(s[gpos], this->krg_.derivative(_sg)); double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(s[gpos], this->krg_.derivative(_sg));
// TODO Check the arguments to the krow- and krog-tables below... // TODO Check the arguments to the krow- and krog-tables below...
double krow = epst->watoil.scaleKr(1.0-s[wpos]-s[gpos], this->krow_(1.0-_sow-this->smin_[gpos]), this->krorw_); // ???? double krow = epst->watoil.scaleKr(1.0-s[wpos]-s[gpos], this->krow_(1.0-_sow-this->smin_[gpos]), this->krorw_); // ????
double dkrow = _dsdsow*epst->watoil.scaleKrDeriv(1.0-s[wpos]-s[gpos], this->krow_.derivative(1.0-_sow-this->smin_[gpos])); // ???? double dkrow = _dsdsow*epst->watoil.scaleKrDeriv(1.0-s[wpos]-s[gpos], this->krow_.derivative(1.0-_sow-this->smin_[gpos])); // ????
//double krog = epst->gasoil.scaleKr(this->krog_(1.0-_sog-this->smin_[wpos]), 1.0-s[wpos]-s[gpos], this->krorg_); // ???? //double krog = epst->gasoil.scaleKr(this->krog_(1.0-_sog-this->smin_[wpos]), 1.0-s[wpos]-s[gpos], this->krorg_); // ????

View File

@@ -31,25 +31,25 @@ namespace Opm
void evalKrDeriv(const double* s, double* kr, double* dkrds) const; void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
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 ...");}
private: private:
}; };
typedef SatFuncStone2<UniformTableLinear<double> > SatFuncStone2Uniform; typedef SatFuncStone2<UniformTableLinear<double> > SatFuncStone2Uniform;
typedef SatFuncStone2<NonuniformTableLinear<double> > SatFuncStone2Nonuniform; typedef SatFuncStone2<NonuniformTableLinear<double> > SatFuncStone2Nonuniform;

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);
} }
}; };

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

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();