mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-11-26 03:00:17 -06:00
Cleaned up code for Splitting method for residual zero finding.
This commit is contained in:
Xavier Raynaud
parent
a26656c6e5
commit
5c78dc8960
@ -220,6 +220,9 @@ namespace Opm
|
||||
|
||||
|
||||
// Residual for s and c. Includes method to compute the gradient.
|
||||
// Compute residual in s (or c) for a given piecewise linear curve (with only one node) in the s-c
|
||||
// plane. The method operator() is used by a 1d root solver.
|
||||
|
||||
struct TransportModelPolymer::Residual
|
||||
{
|
||||
int cell;
|
||||
@ -231,6 +234,15 @@ namespace Opm
|
||||
double outflux; // sum_j max(v_ij, 0)
|
||||
double porosity;
|
||||
double dtpv; // dt/pv(i)
|
||||
double direction[2];
|
||||
double end_point[2];
|
||||
double x_max[2];
|
||||
double x_min[2];
|
||||
double t_out;
|
||||
double t_max; // t_max = t_out + 1
|
||||
double x_out[2];
|
||||
double x[2];
|
||||
bool if_res_s;
|
||||
const TransportModelPolymer& tm;
|
||||
|
||||
Residual(const TransportModelPolymer& tmodel, int cell_index)
|
||||
@ -288,8 +300,9 @@ namespace Opm
|
||||
+ dtpv*(outflux*ff*mc + influx_polymer);
|
||||
}
|
||||
|
||||
|
||||
// Compute gradient using finite difference.
|
||||
void computeGradient(const double* x, double* res, double* gradient, bool if_res_s, const int& method) const
|
||||
void computeGradient(const double* x, double* res, double* gradient, const int& method) const
|
||||
// If if_res_s == true, compute the gradient of s-residual, otherwise, compute gradient of c-residual.
|
||||
// If method == 1, use finite diference
|
||||
// If method == 2, use analytic expresions
|
||||
@ -352,81 +365,17 @@ namespace Opm
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
// Compute residual in s for a given piecewise linear curve (with only one node) in the s-c
|
||||
// plane. The method operator() is used by a 1d root solver.
|
||||
|
||||
struct TransportModelPolymer::ResidualSDir
|
||||
{
|
||||
int cell;
|
||||
double s0;
|
||||
double c0;
|
||||
double cmax0;
|
||||
double influx; // sum_j min(v_ij, 0)*f(s_j)
|
||||
double influx_polymer; // sum_j min(v_ij, 0)*f(s_j)*mc(c_j)
|
||||
double outflux; // sum_j max(v_ij, 0)
|
||||
double porosity;
|
||||
double dtpv; // dt/pv(i)
|
||||
double direction[2];
|
||||
double end_point[2];
|
||||
double x_max[2];
|
||||
double x_min[2];
|
||||
double t_out;
|
||||
double t_max; // t_max = t_out + 1
|
||||
double x_out[2];
|
||||
double x[2];
|
||||
const TransportModelPolymer& tm;
|
||||
|
||||
ResidualSDir(const TransportModelPolymer& tmodel, int cell_index)
|
||||
: tm(tmodel)
|
||||
{
|
||||
cell = cell_index;
|
||||
s0 = tm.saturation_[cell];
|
||||
c0 = tm.concentration_[cell];
|
||||
cmax0 = tm.cmax_[cell];
|
||||
double dflux = -tm.source_[cell];
|
||||
bool src_is_inflow = dflux < 0.0;
|
||||
influx = src_is_inflow ? dflux : 0.0;
|
||||
influx_polymer = src_is_inflow ? dflux*tm.computeMc(tm.inflow_c_) : 0.0;
|
||||
outflux = !src_is_inflow ? dflux : 0.0;
|
||||
dtpv = tm.dt_/tm.porevolume_[cell];
|
||||
porosity = tm.porosity_[cell];
|
||||
|
||||
for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
|
||||
int f = tm.grid_.cell_faces[i];
|
||||
double flux;
|
||||
int other;
|
||||
// Compute cell flux
|
||||
if (cell == tm.grid_.face_cells[2*f]) {
|
||||
flux = tm.darcyflux_[f];
|
||||
other = tm.grid_.face_cells[2*f+1];
|
||||
} else {
|
||||
flux =-tm.darcyflux_[f];
|
||||
other = tm.grid_.face_cells[2*f];
|
||||
}
|
||||
// Add flux to influx or outflux, if interior.
|
||||
if (other != -1) {
|
||||
if (flux < 0.0) {
|
||||
influx += flux*tm.fractionalflow_[other];
|
||||
influx_polymer += flux*tm.fractionalflow_[other]*tm.mc_[other];
|
||||
} else {
|
||||
outflux += flux;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// setup 1d function, which is called by operator()
|
||||
// For a given point x=(s,c) in the s,c plane, set up a piecewise linear curve wich starts
|
||||
// from "x" with slope "direction", hits the bound of the rectangle
|
||||
// [s_min,s_max]x[c_min,c_max] and continue in a straight line to "end_point". The curve is
|
||||
// parametrized by t in [0, t_max], t_out is equal to t when the curve hits the bounding
|
||||
// rectangle, x_out=(s_out, c_out) denotes the values of s and c at that point.
|
||||
|
||||
void setup(const double* x_arg, const double* direction_arg, const double* end_point_arg, const double* x_min_arg, const double* x_max_arg, double& t_max_out, double& t_out_out)
|
||||
void setup(const double* x_arg, const double* direction_arg, const double* end_point_arg, const double* x_min_arg, const double* x_max_arg, const bool& if_res_s_arg, double& t_max_out, double& t_out_out)
|
||||
{
|
||||
double t0, t1;
|
||||
if_res_s = if_res_s_arg;
|
||||
x[0] = x_arg[0];
|
||||
x[1] = x_arg[1];
|
||||
x_max[0] = x_max_arg[0];
|
||||
@ -441,31 +390,31 @@ namespace Opm
|
||||
direction[0] *= -1;
|
||||
direction[1] *= -1;
|
||||
}
|
||||
bool test_dir0 = true;
|
||||
bool test_dir1 = true;
|
||||
bool if_t0 = true;
|
||||
bool if_t1 = true;
|
||||
if (direction[0] > 0) {
|
||||
t0 = (x_max[0] - x[0])/direction[0];
|
||||
} else if (direction[0] < 0) {
|
||||
t0 = (x_min[0] - x[0])/direction[0];
|
||||
} else {
|
||||
test_dir0 = false;
|
||||
if_t0 = false;
|
||||
}
|
||||
if (direction[1] > 0) {
|
||||
t1 = (x_max[1] - x[1])/direction[1];
|
||||
} else if (direction[1] < 0) {
|
||||
t1 = (x_min[1] - x[1])/direction[1];
|
||||
} else {
|
||||
test_dir1 = false;
|
||||
if_t1 = false;
|
||||
}
|
||||
if (test_dir0) {
|
||||
if (test_dir1) {
|
||||
if (if_t0) {
|
||||
if (if_t1) {
|
||||
t_out = std::min(t0, t1);
|
||||
}
|
||||
else {
|
||||
t_out = t0;
|
||||
}
|
||||
} else {
|
||||
if (test_dir1) {
|
||||
if (if_t1) {
|
||||
t_out = t1;
|
||||
}
|
||||
}
|
||||
@ -488,8 +437,6 @@ namespace Opm
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
double operator()(double t) const
|
||||
{
|
||||
double s;
|
||||
@ -501,310 +448,9 @@ namespace Opm
|
||||
s = 1/(t_max-t_out)*((t_max - t)*x_out[0] + end_point[0]*(t - t_out));
|
||||
c = 1/(t_max-t_out)*((t_max - t)*x_out[1] + end_point[1]*(t - t_out));
|
||||
}
|
||||
return s - s0 + dtpv*(outflux*tm.fracFlow(s, c, cell) + influx);
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
|
||||
// Same as ResidualSDir but for the residual in c
|
||||
struct TransportModelPolymer::ResidualCDir
|
||||
{
|
||||
int cell;
|
||||
double s0;
|
||||
double c0;
|
||||
double cmax0;
|
||||
double influx; // sum_j min(v_ij, 0)*f(s_j)
|
||||
double influx_polymer; // sum_j min(v_ij, 0)*f(s_j)*mc(c_j)
|
||||
double outflux; // sum_j max(v_ij, 0)
|
||||
double porosity;
|
||||
double dtpv; // dt/pv(i)
|
||||
double direction[2];
|
||||
double end_point[2];
|
||||
double t_out;
|
||||
double t_max; // t_max = t_out + 1
|
||||
double x_out[2];
|
||||
double x_min[2];
|
||||
double x_max[2];
|
||||
double x[2];
|
||||
const TransportModelPolymer& tm;
|
||||
|
||||
ResidualCDir(const TransportModelPolymer& tmodel, int cell_index)
|
||||
: tm(tmodel)
|
||||
{
|
||||
cell = cell_index;
|
||||
s0 = tm.saturation_[cell];
|
||||
c0 = tm.concentration_[cell];
|
||||
cmax0 = tm.cmax_[cell];
|
||||
double dflux = -tm.source_[cell];
|
||||
bool src_is_inflow = dflux < 0.0;
|
||||
influx = src_is_inflow ? dflux : 0.0;
|
||||
influx_polymer = src_is_inflow ? dflux*tm.computeMc(tm.inflow_c_) : 0.0;
|
||||
outflux = !src_is_inflow ? dflux : 0.0;
|
||||
dtpv = tm.dt_/tm.porevolume_[cell];
|
||||
porosity = tm.porosity_[cell];
|
||||
|
||||
for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
|
||||
int f = tm.grid_.cell_faces[i];
|
||||
double flux;
|
||||
int other;
|
||||
// Compute cell flux
|
||||
if (cell == tm.grid_.face_cells[2*f]) {
|
||||
flux = tm.darcyflux_[f];
|
||||
other = tm.grid_.face_cells[2*f+1];
|
||||
} else {
|
||||
flux =-tm.darcyflux_[f];
|
||||
other = tm.grid_.face_cells[2*f];
|
||||
}
|
||||
// Add flux to influx or outflux, if interior.
|
||||
if (other != -1) {
|
||||
if (flux < 0.0) {
|
||||
influx += flux*tm.fractionalflow_[other];
|
||||
influx_polymer += flux*tm.fractionalflow_[other]*tm.mc_[other];
|
||||
} else {
|
||||
outflux += flux;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void compute_new_x(double* x_new, const double t) {
|
||||
if (t <= t_out) {
|
||||
x_new[0] = x[0] + t*direction[0];
|
||||
x_new[1] = x[1] + t*direction[1];
|
||||
} else {
|
||||
x_new[0] = 1/(t_max-t_out)*((t_max - t)*x_out[0] + end_point[0]*(t - t_out));
|
||||
x_new[1] = 1/(t_max-t_out)*((t_max - t)*x_out[1] + end_point[1]*(t - t_out));
|
||||
}
|
||||
}
|
||||
|
||||
void setup(const double* x_arg, const double* direction_arg, const double* end_point_arg, const double* x_min_arg, const double* x_max_arg, double& t_max_out, double& t_out_out)
|
||||
{
|
||||
bool if_t0 = true;
|
||||
bool if_t1 = true;
|
||||
double t0, t1;
|
||||
x[0] = x_arg[0];
|
||||
x[1] = x_arg[1];
|
||||
x_max[0] = x_max_arg[0];
|
||||
x_max[1] = x_max_arg[1];
|
||||
x_min[0] = x_min_arg[0];
|
||||
x_min[1] = x_min_arg[1];
|
||||
direction[0] = direction_arg[0];
|
||||
direction[1] = direction_arg[1];
|
||||
end_point[0] = end_point_arg[0];
|
||||
end_point[1] = end_point_arg[1];
|
||||
if ((end_point[0]-x[0])*direction[0] + (end_point[1]-x[1])*direction[1] < 0) {
|
||||
direction[0] *= -1;
|
||||
direction[1] *= -1;
|
||||
}
|
||||
if (direction[0] == 0) {
|
||||
if_t0 = false;
|
||||
} else {
|
||||
if (direction[0] > 0) {
|
||||
t0 = (x_max[0] - x[0])/direction[0];
|
||||
} else {
|
||||
t0 = (x_min[0] - x[0])/direction[0];
|
||||
}
|
||||
}
|
||||
if (direction[1] == 0) {
|
||||
if_t1 = false;
|
||||
} else {
|
||||
if (direction[1] > 0) {
|
||||
t1 = (x_max[1] - x[1])/direction[1];
|
||||
} else {
|
||||
t1 = (x_min[1] - x[1])/direction[1];
|
||||
}
|
||||
}
|
||||
if (if_t0 && if_t1) {
|
||||
t_out = std::min(t0, t1);
|
||||
} else {
|
||||
if (if_t0) {
|
||||
t_out = t0;
|
||||
} else {
|
||||
t_out = t1;
|
||||
}
|
||||
}
|
||||
x_out[0] = x[0] + t_out*direction[0];
|
||||
x_out[1] = x[1] + t_out*direction[1];
|
||||
t_max = t_out + 1;
|
||||
t_max_out = t_max;
|
||||
t_out_out = t_out;
|
||||
}
|
||||
|
||||
double operator()(double t) const
|
||||
{
|
||||
double s;
|
||||
double c;
|
||||
if (t <= t_out) {
|
||||
s = x[0] + t*direction[0];
|
||||
c = x[1] + t*direction[1];
|
||||
} else {
|
||||
s = 1/(t_max-t_out)*((t_max - t)*x_out[0] + end_point[0]*(t - t_out));
|
||||
c = 1/(t_max-t_out)*((t_max - t)*x_out[1] + end_point[1]*(t - t_out));
|
||||
}
|
||||
double ff = tm.fracFlow(s, c, cell);
|
||||
double mc = tm.computeMc(c);
|
||||
double dps = tm.polyprops_.dps;
|
||||
double rhor = tm.polyprops_.rhor;
|
||||
double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
|
||||
double ads = tm.polyprops_.adsorbtion(std::max(c, cmax0));
|
||||
return (s - dps)*c - (s0 - dps)*c0
|
||||
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
|
||||
+ dtpv*(outflux*ff*mc + influx_polymer);
|
||||
}
|
||||
};
|
||||
|
||||
struct TransportModelPolymer::ResidualDir
|
||||
{
|
||||
int cell;
|
||||
int s_or_c; // s_or_c = 0 if s direction, s_or_c = 1 if c direction,
|
||||
double s0;
|
||||
double c0;
|
||||
double cmax0;
|
||||
double influx; // sum_j min(v_ij, 0)*f(s_j)
|
||||
double influx_polymer; // sum_j min(v_ij, 0)*f(s_j)*mc(c_j)
|
||||
double outflux; // sum_j max(v_ij, 0)
|
||||
double porosity;
|
||||
double dtpv; // dt/pv(i)
|
||||
double direction[2];
|
||||
double end_point[2];
|
||||
double x_max[2];
|
||||
double x_min[2];
|
||||
double t_out;
|
||||
double t_max; // t_max = t_out + 1
|
||||
double x_out[2];
|
||||
double x[2];
|
||||
const TransportModelPolymer& tm;
|
||||
|
||||
ResidualDir(const TransportModelPolymer& tmodel, int cell_index)
|
||||
: tm(tmodel)
|
||||
{
|
||||
cell = cell_index;
|
||||
s0 = tm.saturation_[cell];
|
||||
c0 = tm.concentration_[cell];
|
||||
cmax0 = tm.cmax_[cell];
|
||||
double dflux = -tm.source_[cell];
|
||||
bool src_is_inflow = dflux < 0.0;
|
||||
influx = src_is_inflow ? dflux : 0.0;
|
||||
influx_polymer = src_is_inflow ? dflux*tm.computeMc(tm.inflow_c_) : 0.0;
|
||||
outflux = !src_is_inflow ? dflux : 0.0;
|
||||
dtpv = tm.dt_/tm.porevolume_[cell];
|
||||
porosity = tm.porosity_[cell];
|
||||
|
||||
for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
|
||||
int f = tm.grid_.cell_faces[i];
|
||||
double flux;
|
||||
int other;
|
||||
// Compute cell flux
|
||||
if (cell == tm.grid_.face_cells[2*f]) {
|
||||
flux = tm.darcyflux_[f];
|
||||
other = tm.grid_.face_cells[2*f+1];
|
||||
} else {
|
||||
flux =-tm.darcyflux_[f];
|
||||
other = tm.grid_.face_cells[2*f];
|
||||
}
|
||||
// Add flux to influx or outflux, if interior.
|
||||
if (other != -1) {
|
||||
if (flux < 0.0) {
|
||||
influx += flux*tm.fractionalflow_[other];
|
||||
influx_polymer += flux*tm.fractionalflow_[other]*tm.mc_[other];
|
||||
} else {
|
||||
outflux += flux;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// For a given point x=(s,c) in the s,c plane, set up a piecewise linear curve wich starts
|
||||
// from "x" with slope "direction", hits the bound of the rectangle
|
||||
// [s_min,s_max]x[c_min,c_max] and continue in a straight line to "end_point". The curve is
|
||||
// parametrized by t in [0, t_max], t_out is equal to t when the curve hits the bounding
|
||||
// rectangle, x_out=(s_out, c_out) denotes the values of s and c at that point.
|
||||
|
||||
void setup(const double* x_arg, const double* direction_arg, const double* end_point_arg, const double* x_min_arg, const double* x_max_arg, const int& s_or_c_arg, double& t_max_out, double& t_out_out)
|
||||
{
|
||||
bool if_t0 = true;
|
||||
bool if_t1 = true;
|
||||
double t0, t1;
|
||||
|
||||
s_or_c = s_or_c_arg;
|
||||
x[0] = x_arg[0];
|
||||
x[1] = x_arg[1];
|
||||
x_max[0] = x_max_arg[0];
|
||||
x_max[1] = x_max_arg[1];
|
||||
x_min[0] = x_min_arg[0];
|
||||
x_min[1] = x_min_arg[1];
|
||||
direction[0] = direction_arg[0];
|
||||
direction[1] = direction_arg[1];
|
||||
end_point[0] = end_point_arg[0];
|
||||
end_point[1] = end_point_arg[1];
|
||||
if ((end_point[0]-x[0])*direction[0] + (end_point[1]-x[1])*direction[1] < 0) {
|
||||
direction[0] *= -1;
|
||||
direction[1] *= -1;
|
||||
}
|
||||
if (direction[0] == 0) {
|
||||
if_t0 = false;
|
||||
} else {
|
||||
if (direction[0] > 0) {
|
||||
t0 = (x_max[0] - x[0])/direction[0];
|
||||
} else {
|
||||
t0 = (x_min[0] - x[0])/direction[0];
|
||||
}
|
||||
}
|
||||
if (direction[1] == 0) {
|
||||
if_t1 = false;
|
||||
} else {
|
||||
if (direction[1] > 0) {
|
||||
t1 = (x_max[1] - x[1])/direction[1];
|
||||
} else {
|
||||
t1 = (x_min[1] - x[1])/direction[1];
|
||||
}
|
||||
}
|
||||
if (if_t0 && if_t1) {
|
||||
t_out = std::min(t0, t1);
|
||||
} else {
|
||||
if (if_t0) {
|
||||
t_out = t0;
|
||||
} else {
|
||||
t_out = t1;
|
||||
}
|
||||
}
|
||||
x_out[0] = x[0] + t_out*direction[0];
|
||||
x_out[1] = x[1] + t_out*direction[1];
|
||||
t_max = t_out + 1;
|
||||
t_max_out = t_max;
|
||||
t_out_out = t_out;
|
||||
}
|
||||
|
||||
|
||||
// Compute x=(s,c) for a given t (t is the parameter for the piecewise linear curve)
|
||||
void compute_new_x(double* x_new, const double t) {
|
||||
if (t <= t_out) {
|
||||
x_new[0] = x[0] + t*direction[0];
|
||||
x_new[1] = x[1] + t*direction[1];
|
||||
} else {
|
||||
x_new[0] = 1/(t_max-t_out)*((t_max - t)*x_out[0] + end_point[0]*(t - t_out));
|
||||
x_new[1] = 1/(t_max-t_out)*((t_max - t)*x_out[1] + end_point[1]*(t - t_out));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
double operator()(double t) const
|
||||
{
|
||||
double s;
|
||||
double c;
|
||||
if (t <= t_out) {
|
||||
s = x[0] + t*direction[0];
|
||||
c = x[1] + t*direction[1];
|
||||
} else {
|
||||
s = 1/(t_max-t_out)*((t_max - t)*x_out[0] + end_point[0]*(t - t_out));
|
||||
c = 1/(t_max-t_out)*((t_max - t)*x_out[1] + end_point[1]*(t - t_out));
|
||||
}
|
||||
if (s_or_c == 0) {
|
||||
if (if_res_s) {
|
||||
return s - s0 + dtpv*(outflux*tm.fracFlow(s, c, cell) + influx);
|
||||
} else if (s_or_c == 1) {
|
||||
} else {
|
||||
double ff = tm.fracFlow(s, c, cell);
|
||||
double mc = tm.computeMc(c);
|
||||
double dps = tm.polyprops_.dps;
|
||||
@ -814,12 +460,12 @@ namespace Opm
|
||||
return (s - dps)*c - (s0 - dps)*c0
|
||||
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
|
||||
+ dtpv*(outflux*ff*mc + influx_polymer);
|
||||
} else {
|
||||
std::cout << "problem!" << std::endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
void TransportModelPolymer::solveSingleCell(const int cell)
|
||||
{
|
||||
if (method_ == 1) {
|
||||
@ -857,8 +503,8 @@ namespace Opm
|
||||
int iters_used_split = 0;
|
||||
|
||||
Residual residual(*this, cell);
|
||||
ResidualSDir residual_s_dir(*this, cell);
|
||||
ResidualCDir residual_c_dir(*this, cell);
|
||||
// ResidualSDir residual_s_dir(*this, cell);
|
||||
// ResidualCDir residual_c_dir(*this, cell);
|
||||
// const int sdir = 0;
|
||||
// const int cdir = 1;
|
||||
// ResidualDir residual_dir(*this, cell);
|
||||
@ -916,8 +562,8 @@ namespace Opm
|
||||
//
|
||||
end_point[0] = x_max[0];
|
||||
end_point[1] = x_max[1];
|
||||
residual_c_dir.setup(x, direction, end_point, x_min, x_max, t_max, t_out);
|
||||
if (residual_c_dir(t_out) >= 0) {
|
||||
residual.setup(x, direction, end_point, x_min, x_max, false, t_max, t_out);
|
||||
if (residual(t_out) >= 0) {
|
||||
t_max = t_out;
|
||||
}
|
||||
} else {
|
||||
@ -929,17 +575,17 @@ namespace Opm
|
||||
//
|
||||
end_point[0] = x_min[0];
|
||||
end_point[1] = x_min[1];
|
||||
residual_c_dir.setup(x, direction, end_point, x_min, x_max, t_max, t_out);
|
||||
if (residual_c_dir(t_out) <= 0) {
|
||||
residual.setup(x, direction, end_point, x_min, x_max, false, t_max, t_out);
|
||||
if (residual(t_out) <= 0) {
|
||||
t_max = t_out;
|
||||
}
|
||||
}
|
||||
t = modifiedRegulaFalsi(residual_c_dir, 0., t_max, max_iters_falsi, tol, iters_used_falsi);
|
||||
if (std::abs(residual_c_dir(t)) > tol) {
|
||||
t = modifiedRegulaFalsi(residual, 0., t_max, max_iters_falsi, tol, iters_used_falsi);
|
||||
if (std::abs(residual(t)) > tol) {
|
||||
std::cout << "modifiedRegulaFalsi did not produce result under tolerance." << std::endl;
|
||||
}
|
||||
residual_c_dir.compute_new_x(x, t);
|
||||
residual.computeGradient(x, res, gradient, false, 1);
|
||||
residual.compute_new_x(x, t);
|
||||
residual.computeGradient(x, res, gradient, 1);
|
||||
direction[0] = gradient[1];
|
||||
direction[1] = -gradient[0];
|
||||
res_s_done = false;
|
||||
@ -947,24 +593,24 @@ namespace Opm
|
||||
if (res[0] < 0) {
|
||||
end_point[0] = x_max[0];
|
||||
end_point[1] = x_min[1];
|
||||
residual_s_dir.setup(x, direction, end_point, x_min, x_max, t_max, t_out);
|
||||
if (residual_s_dir(t_out) >= 0) {
|
||||
residual.setup(x, direction, end_point, x_min, x_max, true, t_max, t_out);
|
||||
if (residual(t_out) >= 0) {
|
||||
t_max = t_out;
|
||||
}
|
||||
} else {
|
||||
end_point[0] = x_min[0];
|
||||
end_point[1] = x_max[1];
|
||||
residual_s_dir.setup(x, direction, end_point, x_min, x_max, t_max, t_out);
|
||||
if (residual_s_dir(t_out) <= 0) {
|
||||
residual.setup(x, direction, end_point, x_min, x_max, true, t_max, t_out);
|
||||
if (residual(t_out) <= 0) {
|
||||
t_max = t_out;
|
||||
}
|
||||
}
|
||||
t = modifiedRegulaFalsi(residual_s_dir, 0., t_max, max_iters_falsi, tol, iters_used_falsi);
|
||||
if (std::abs(residual_s_dir(t)) > tol) {
|
||||
t = modifiedRegulaFalsi(residual, 0., t_max, max_iters_falsi, tol, iters_used_falsi);
|
||||
if (std::abs(residual(t)) > tol) {
|
||||
std::cout << "modifiedRegulaFalsi did not produce result under tolerance." << std::endl;
|
||||
}
|
||||
residual_s_dir.compute_new_x(x, t);
|
||||
residual.computeGradient(x, res, gradient, true, 1);
|
||||
residual.compute_new_x(x, t);
|
||||
residual.computeGradient(x, res, gradient, 1);
|
||||
res_s_done = true;
|
||||
direction[0] = -gradient[1];
|
||||
direction[1] = gradient[0];
|
||||
|
||||
Loading…
Reference in New Issue
Block a user