Put brackets for nonlinear solvers into NonlinearSolverCtrl struct.

This commit is contained in:
Atgeirr Flø Rasmussen 2012-02-08 13:37:52 +01:00
parent 08fb3df82c
commit c0ded27f8d
4 changed files with 13 additions and 10 deletions

View File

@ -90,11 +90,11 @@ ridders (double (*G)(double, void*), void *data, struct NonlinearSolverCtrl *ctr
G2 = G(s2, data); G2 = G(s2, data);
if (fabs(G2) < ctrl->nltolerance) { return s2; } if (fabs(G2) < ctrl->nltolerance) { return s2; }
s0 = 0.0; s0 = ctrl->min_bracket;
G0 = G(s0, data); G0 = G(s0, data);
if (fabs(G0) < ctrl->nltolerance) { return s0; } if (fabs(G0) < ctrl->nltolerance) { return s0; }
s1 = 1.0; s1 = ctrl->max_bracket;
G1 = G(s1, data); G1 = G(s1, data);
if (fabs(G1) < ctrl->nltolerance) { return s1; } if (fabs(G1) < ctrl->nltolerance) { return s1; }
@ -190,12 +190,12 @@ regulafalsi (double (*G)(double, void*), void *data, struct NonlinearSolverCtrl
Gn = G(sn, data); Gn = G(sn, data);
if (fabs(Gn) < ctrl->nltolerance) { return sn; } if (fabs(Gn) < ctrl->nltolerance) { return sn; }
/* Initial guess is interval [0,1] of course */ /* Initial guess is interval [s0,s1] */
s0 = 0.0; s0 = ctrl->min_bracket;
G0 = G(s0, data); G0 = G(s0, data);
if (fabs(G0) < ctrl->nltolerance) { return s0; } if (fabs(G0) < ctrl->nltolerance) { return s0; }
s1 = 1.0; s1 = ctrl->max_bracket;
G1 = G(s1, data); G1 = G(s1, data);
if (fabs(G1) < ctrl->nltolerance) { return s1; } 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); Gn = G(sn, data);
if (fabs(Gn) < ctrl->nltolerance) { return sn; } if (fabs(Gn) < ctrl->nltolerance) { return sn; }
/* Initial guess is interval [0,1] of course */ /* Initial guess is interval [s0,s1] */
s0 = 0.0; s0 = ctrl->max_bracket;
G0 = G(s0, data); G0 = G(s0, data);
if (fabs(G0) < ctrl->nltolerance) { return s0; } if (fabs(G0) < ctrl->nltolerance) { return s0; }
s1 = 1.0; s1 = ctrl->max_bracket;
G1 = G(s1, data); G1 = G(s1, data);
if (fabs(G1) < ctrl->nltolerance) { return s1; } if (fabs(G1) < ctrl->nltolerance) { return s1; }

View File

@ -29,9 +29,11 @@ extern "C" {
struct NonlinearSolverCtrl struct NonlinearSolverCtrl
{ {
enum Method {RIDDERS, REGULAFALSI, BISECTION} method; enum Method {RIDDERS, REGULAFALSI, BISECTION} method;
double nltolerance; double nltolerance;
int maxiterations; int maxiterations;
double min_bracket;
double max_bracket;
double initialguess; double initialguess;
int iterations; /* set by solver */ int iterations; /* set by solver */
double residual; /* set by solver */ double residual; /* set by solver */

View File

@ -78,7 +78,6 @@ void solvecell(void *data, struct NonlinearSolverCtrl *ctrl, int cell)
{ {
struct SolverData *d = (struct SolverData*) data; struct SolverData *d = (struct SolverData*) data;
struct Parameters prm = get_parameters(d, cell); struct Parameters prm = get_parameters(d, cell);
d->saturation[cell] = find_zero(residual, &prm, ctrl); d->saturation[cell] = find_zero(residual, &prm, ctrl);
// double ff1 = fluxfun_props(d->saturation[cell], cell, d->props); // double ff1 = fluxfun_props(d->saturation[cell], cell, d->props);
// double ff2 = fluxfun(d->saturation[cell], -999); // double ff2 = fluxfun(d->saturation[cell], -999);

View File

@ -48,6 +48,8 @@ void twophasetransport(
ctrl.nltolerance = 1e-9; ctrl.nltolerance = 1e-9;
ctrl.method = NonlinearSolverCtrl::REGULAFALSI; ctrl.method = NonlinearSolverCtrl::REGULAFALSI;
ctrl.initialguess = 0.5; ctrl.initialguess = 0.5;
ctrl.min_bracket = 0.0;
ctrl.max_bracket = 1.0;
/* Assume all strong components are single-cell domains. */ /* Assume all strong components are single-cell domains. */
for(i=0; i<grid->number_of_cells; ++i) for(i=0; i<grid->number_of_cells; ++i)