diff --git a/opm/core/transport/reorder/nlsolvers.c b/opm/core/transport/reorder/nlsolvers.c index 2f958702..02672117 100644 --- a/opm/core/transport/reorder/nlsolvers.c +++ b/opm/core/transport/reorder/nlsolvers.c @@ -90,11 +90,11 @@ ridders (double (*G)(double, void*), void *data, struct NonlinearSolverCtrl *ctr G2 = G(s2, data); if (fabs(G2) < ctrl->nltolerance) { return s2; } - s0 = 0.0; + s0 = ctrl->min_bracket; G0 = G(s0, data); if (fabs(G0) < ctrl->nltolerance) { return s0; } - s1 = 1.0; + s1 = ctrl->max_bracket; G1 = G(s1, data); if (fabs(G1) < ctrl->nltolerance) { return s1; } @@ -190,12 +190,12 @@ regulafalsi (double (*G)(double, void*), void *data, struct NonlinearSolverCtrl Gn = G(sn, data); if (fabs(Gn) < ctrl->nltolerance) { return sn; } - /* Initial guess is interval [0,1] of course */ - s0 = 0.0; + /* Initial guess is interval [s0,s1] */ + s0 = ctrl->min_bracket; G0 = G(s0, data); if (fabs(G0) < ctrl->nltolerance) { return s0; } - s1 = 1.0; + s1 = ctrl->max_bracket; G1 = G(s1, data); if (fabs(G1) < ctrl->nltolerance) { return s1; } @@ -275,12 +275,12 @@ bisection (double (*G)(double, void*), void *data, struct NonlinearSolverCtrl *c Gn = G(sn, data); if (fabs(Gn) < ctrl->nltolerance) { return sn; } - /* Initial guess is interval [0,1] of course */ - s0 = 0.0; + /* Initial guess is interval [s0,s1] */ + s0 = ctrl->max_bracket; G0 = G(s0, data); if (fabs(G0) < ctrl->nltolerance) { return s0; } - s1 = 1.0; + s1 = ctrl->max_bracket; G1 = G(s1, data); if (fabs(G1) < ctrl->nltolerance) { return s1; } diff --git a/opm/core/transport/reorder/nlsolvers.h b/opm/core/transport/reorder/nlsolvers.h index a21b1333..2f839d81 100644 --- a/opm/core/transport/reorder/nlsolvers.h +++ b/opm/core/transport/reorder/nlsolvers.h @@ -29,9 +29,11 @@ extern "C" { struct NonlinearSolverCtrl { - enum Method {RIDDERS, REGULAFALSI, BISECTION} method; + enum Method {RIDDERS, REGULAFALSI, BISECTION} method; double nltolerance; int maxiterations; + double min_bracket; + double max_bracket; double initialguess; int iterations; /* set by solver */ double residual; /* set by solver */ diff --git a/opm/core/transport/reorder/twophase.cpp b/opm/core/transport/reorder/twophase.cpp index a42c2c2f..a4d93bc3 100644 --- a/opm/core/transport/reorder/twophase.cpp +++ b/opm/core/transport/reorder/twophase.cpp @@ -78,7 +78,6 @@ void solvecell(void *data, struct NonlinearSolverCtrl *ctrl, int cell) { struct SolverData *d = (struct SolverData*) data; struct Parameters prm = get_parameters(d, cell); - d->saturation[cell] = find_zero(residual, &prm, ctrl); // double ff1 = fluxfun_props(d->saturation[cell], cell, d->props); // double ff2 = fluxfun(d->saturation[cell], -999); diff --git a/opm/core/transport/reorder/twophasetransport.cpp b/opm/core/transport/reorder/twophasetransport.cpp index 3b8aafb6..3c4f364b 100644 --- a/opm/core/transport/reorder/twophasetransport.cpp +++ b/opm/core/transport/reorder/twophasetransport.cpp @@ -48,6 +48,8 @@ void twophasetransport( ctrl.nltolerance = 1e-9; ctrl.method = NonlinearSolverCtrl::REGULAFALSI; ctrl.initialguess = 0.5; + ctrl.min_bracket = 0.0; + ctrl.max_bracket = 1.0; /* Assume all strong components are single-cell domains. */ for(i=0; inumber_of_cells; ++i)