mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
cleaned up TransportModelPolymer (default branch)
This commit is contained in:
Xavier Raynaud
@@ -127,7 +127,7 @@ namespace
|
|||||||
CurveInSCPlane();
|
CurveInSCPlane();
|
||||||
void setup(const double* x, const double* direction,
|
void setup(const double* x, const double* direction,
|
||||||
const double* end_point, const double* x_min,
|
const double* end_point, const double* x_min,
|
||||||
const double* x_max, const bool adjust_dir, const double tol,
|
const double* x_max, const double tol,
|
||||||
double& t_max_out, double& t_out_out);
|
double& t_max_out, double& t_out_out);
|
||||||
void computeXOfT(double*, const double) const;
|
void computeXOfT(double*, const double) const;
|
||||||
|
|
||||||
@@ -627,6 +627,9 @@ namespace Opm
|
|||||||
case Newton:
|
case Newton:
|
||||||
solveSingleCellNewton(cell);
|
solveSingleCellNewton(cell);
|
||||||
break;
|
break;
|
||||||
|
case Gradient:
|
||||||
|
solveSingleCellGradient(cell);
|
||||||
|
break;
|
||||||
default:
|
default:
|
||||||
THROW("Unknown method " << method_);
|
THROW("Unknown method " << method_);
|
||||||
}
|
}
|
||||||
@@ -663,7 +666,7 @@ namespace Opm
|
|||||||
// Newton method, where we first try a Newton step. Then, if it does not work well, we look for
|
// Newton method, where we first try a Newton step. Then, if it does not work well, we look for
|
||||||
// the zero of either the residual in s or the residual in c along a specified piecewise linear
|
// the zero of either the residual in s or the residual in c along a specified piecewise linear
|
||||||
// curve. In these cases, we can use a robust 1d solver.
|
// curve. In these cases, we can use a robust 1d solver.
|
||||||
void TransportModelPolymer::solveSingleCellNewtonGradient(int cell)
|
void TransportModelPolymer::solveSingleCellGradient(int cell)
|
||||||
{
|
{
|
||||||
int iters_used_falsi = 0;
|
int iters_used_falsi = 0;
|
||||||
const int max_iters_split = maxit_;
|
const int max_iters_split = maxit_;
|
||||||
@@ -672,7 +675,6 @@ namespace Opm
|
|||||||
// Check if current state is an acceptable solution.
|
// Check if current state is an acceptable solution.
|
||||||
ResidualEquation res_eq(*this, cell);
|
ResidualEquation res_eq(*this, cell);
|
||||||
double x[2] = {saturation_[cell], concentration_[cell]};
|
double x[2] = {saturation_[cell], concentration_[cell]};
|
||||||
// double x[2] = {0.5, polyprops_.cMax()/2.0};
|
|
||||||
double res[2];
|
double res[2];
|
||||||
double mc;
|
double mc;
|
||||||
double ff;
|
double ff;
|
||||||
@@ -682,25 +684,6 @@ namespace Opm
|
|||||||
fractionalflow_[cell] = ff;
|
fractionalflow_[cell] = ff;
|
||||||
mc_[cell] = mc;
|
mc_[cell] = mc;
|
||||||
return;
|
return;
|
||||||
} else {
|
|
||||||
// Can be advantageous to reset saturation and concentration when we
|
|
||||||
// are close to the end points.
|
|
||||||
/*
|
|
||||||
x[0] = saturation_[]-res[0];
|
|
||||||
if((x[0]>1) || (x[0]<0)){
|
|
||||||
x[0] = 0.5;
|
|
||||||
}
|
|
||||||
if(x[0]>0){
|
|
||||||
x[1] = concentration_[cell]*saturation_[cell]-res[1];
|
|
||||||
x[1] = x[1]/x[0];
|
|
||||||
}else{
|
|
||||||
x[1]=0;
|
|
||||||
}
|
|
||||||
*/
|
|
||||||
/*
|
|
||||||
x[0]=0.5;x[1]=polyprops_.cMax()/2.0;
|
|
||||||
res_eq.computeResidual(x, res, mc, ff);
|
|
||||||
*/
|
|
||||||
}
|
}
|
||||||
|
|
||||||
double x_min[2] = { 0.0, 0.0 };
|
double x_min[2] = { 0.0, 0.0 };
|
||||||
@@ -711,163 +694,102 @@ namespace Opm
|
|||||||
double direction[2];
|
double direction[2];
|
||||||
double end_point[2];
|
double end_point[2];
|
||||||
double gradient[2];
|
double gradient[2];
|
||||||
bool unsuccessfull_newton_step = true;
|
|
||||||
double x_new[2];
|
|
||||||
double res_new[2];
|
|
||||||
ResSOnCurve res_s_on_curve(res_eq);
|
ResSOnCurve res_s_on_curve(res_eq);
|
||||||
ResCOnCurve res_c_on_curve(res_eq);
|
ResCOnCurve res_c_on_curve(res_eq);
|
||||||
bool if_res_s;
|
bool if_res_s;
|
||||||
int counter_drop_newton = 0;
|
|
||||||
bool not_so_successfull_newton_step = false;
|
|
||||||
|
|
||||||
|
|
||||||
while ((norm(res) > tol_) && (iters_used_split < max_iters_split)) {
|
while ((norm(res) > tol_) && (iters_used_split < max_iters_split)) {
|
||||||
// We first try a Newton step
|
if (std::abs(res[0]) < std::abs(res[1])) {
|
||||||
if (counter_drop_newton == 0 && solveNewtonStep(x, res_eq, res, x_new)) {
|
if (res[0] < -tol_) {
|
||||||
check_interval(x_new, x_min, x_max);
|
direction[0] = x_max[0] - x[0];
|
||||||
res_eq.computeResidual(x_new, res_new, mc, ff);
|
direction[1] = x_min[1] - x[1];
|
||||||
|
if_res_s = true;
|
||||||
|
} else if (res[0] > tol_) {
|
||||||
|
direction[0] = x_min[0] - x[0];
|
||||||
|
direction[1] = x_max[1] - x[1];
|
||||||
|
if_res_s = true;
|
||||||
|
} else {
|
||||||
|
res_eq.computeGradientResS(x, res, gradient);
|
||||||
|
direction[0] = -gradient[1];
|
||||||
|
direction[1] = gradient[0];
|
||||||
|
if_res_s = false;
|
||||||
|
}
|
||||||
|
} else {
|
||||||
|
if (res[1] < -tol_) {
|
||||||
|
direction[0] = x_max[0] - x[0];
|
||||||
|
direction[1] = x_max[1] - x[1];
|
||||||
|
if_res_s = false;
|
||||||
|
} else if (res[1] > tol_) {
|
||||||
|
direction[0] = x_min[0] - x[0];
|
||||||
|
direction[1] = x_min[1] - x[1];
|
||||||
|
if_res_s = false;
|
||||||
|
} else {
|
||||||
|
res_eq.computeGradientResC(x, res, gradient);
|
||||||
|
direction[0] = -gradient[1];
|
||||||
|
direction[1] = gradient[0];
|
||||||
|
if_res_s = true;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (if_res_s) {
|
||||||
|
if (res[0] < 0) {
|
||||||
|
end_point[0] = x_max[0];
|
||||||
|
end_point[1] = x_min[1];
|
||||||
|
res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out);
|
||||||
|
if (res_s_on_curve(t_out) >= 0) {
|
||||||
|
t_max = t_out;
|
||||||
|
}
|
||||||
|
} else {
|
||||||
|
end_point[0] = x_min[0];
|
||||||
|
end_point[1] = x_max[1];
|
||||||
|
res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out);
|
||||||
|
if (res_s_on_curve(t_out) <= 0) {
|
||||||
|
t_max = t_out;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
// Note: In some experiments modifiedRegularFalsi does not yield a result under the given tolerance.
|
||||||
|
t = RootFinder::solve(res_s_on_curve, 0., t_max, maxit_, tol_, iters_used_falsi);
|
||||||
|
res_s_on_curve.curve.computeXOfT(x, t);
|
||||||
|
} else {
|
||||||
|
if (res[1] < 0) {
|
||||||
|
end_point[0] = x_max[0];
|
||||||
|
end_point[1] = x_max[1];
|
||||||
|
res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out);
|
||||||
|
if (res_c_on_curve(t_out) >= 0) {
|
||||||
|
t_max = t_out;
|
||||||
|
}
|
||||||
|
} else {
|
||||||
|
end_point[0] = x_min[0];
|
||||||
|
end_point[1] = x_min[1];
|
||||||
|
res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out);
|
||||||
|
if (res_c_on_curve(t_out) <= 0) {
|
||||||
|
t_max = t_out;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
t = RootFinder::solve(res_c_on_curve, 0., t_max, maxit_, tol_, iters_used_falsi);
|
||||||
|
res_c_on_curve.curve.computeXOfT(x, t);
|
||||||
|
|
||||||
unsuccessfull_newton_step = false;
|
}
|
||||||
not_so_successfull_newton_step = false;
|
res_eq.computeResidual(x, res, mc, ff);
|
||||||
|
iters_used_split += 1;
|
||||||
if (norm(res_new) > norm(res) ||
|
}
|
||||||
x_new[0] < x_min[0] ||
|
|
||||||
x_new[1] < x_min[1] ||
|
|
||||||
x_new[0] > x_max[0] ||
|
|
||||||
x_new[1] > x_max[1]) {
|
|
||||||
unsuccessfull_newton_step = true;
|
|
||||||
} else {
|
|
||||||
x[0] = x_new[0];
|
|
||||||
x[1] = x_new[1];
|
|
||||||
if (norm(res_new) > 1e-1*norm(res)) {
|
|
||||||
// We are close to the solution and Newton does not perform well.
|
|
||||||
// Then, we drop Newton for a given number of iterations.
|
|
||||||
not_so_successfull_newton_step = true;
|
|
||||||
counter_drop_newton = 4;
|
|
||||||
}
|
|
||||||
res[0] = res_new[0];
|
|
||||||
res[1] = res_new[1];
|
|
||||||
iters_used_split += 1;
|
|
||||||
}
|
|
||||||
} else {
|
|
||||||
unsuccessfull_newton_step = true;
|
|
||||||
}
|
|
||||||
|
|
||||||
if (not_so_successfull_newton_step || unsuccessfull_newton_step) {
|
|
||||||
// Newton was not satisfactory. We start 1d solvers.
|
|
||||||
if (not_so_successfull_newton_step) {
|
|
||||||
counter_drop_newton -= 1;
|
|
||||||
}
|
|
||||||
|
|
||||||
// General comment on the zero curves of the s and c residuals:
|
|
||||||
// Typically res_s(s,c)=0 defines an increasing curve in the s-c plane while
|
|
||||||
// res_c(s,c)=0 defines a decreasing curve. However, we do not assume that in the algorithm.
|
|
||||||
// We know that res_s(x_top_left)<0, res_s(x_bottom_right)>0
|
|
||||||
// and res_c(x_bottom_left)<0, res_c(x_top_right)>0
|
|
||||||
// Here, "top", "bottom", ... refer to the corner of the admissible box of (s,c) values.
|
|
||||||
// We use these results to construct a 1d curve for which we are sure that res_s or res_c change sign
|
|
||||||
// and which can therefore be used by a 1d solver.
|
|
||||||
|
|
||||||
// We start with the zero curve of the s and r residual we are closest to.
|
|
||||||
bool adjust_dir = true;
|
|
||||||
if (std::abs(res[0]) < std::abs(res[1])) {
|
|
||||||
if (res[0] < -tol_) {
|
|
||||||
direction[0] = x_max[0] - x[0];
|
|
||||||
direction[1] = x_min[1] - x[1];
|
|
||||||
adjust_dir = true;
|
|
||||||
if_res_s = true;
|
|
||||||
} else if (res[0] > tol_) {
|
|
||||||
direction[0] = x_min[0] - x[0];
|
|
||||||
direction[1] = x_max[1] - x[1];
|
|
||||||
adjust_dir = true;
|
|
||||||
if_res_s = true;
|
|
||||||
} else {
|
|
||||||
res_eq.computeGradientResS(x, res, gradient);
|
|
||||||
direction[0] = -gradient[1];
|
|
||||||
direction[1] = gradient[0];
|
|
||||||
adjust_dir = true;
|
|
||||||
if_res_s = false;
|
|
||||||
}
|
|
||||||
} else {
|
|
||||||
if (res[1] < -tol_) {
|
|
||||||
direction[0] = x_max[0] - x[0];
|
|
||||||
direction[1] = x_max[1] - x[1];
|
|
||||||
adjust_dir = true;
|
|
||||||
if_res_s = false;
|
|
||||||
} else if (res[1] > tol_) {
|
|
||||||
direction[0] = x_min[0] - x[0];
|
|
||||||
direction[1] = x_min[1] - x[1];
|
|
||||||
adjust_dir = true;
|
|
||||||
if_res_s = false;
|
|
||||||
} else {
|
|
||||||
res_eq.computeGradientResC(x, res, gradient);
|
|
||||||
direction[0] = -gradient[1];
|
|
||||||
direction[1] = gradient[0];
|
|
||||||
adjust_dir = true;
|
|
||||||
if_res_s = true;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if (if_res_s) {
|
|
||||||
if (res[0] < 0) {
|
|
||||||
end_point[0] = x_max[0];
|
|
||||||
end_point[1] = x_min[1];
|
|
||||||
res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, adjust_dir, tol_, t_max, t_out);
|
|
||||||
if (res_s_on_curve(t_out) >= 0) {
|
|
||||||
t_max = t_out;
|
|
||||||
}
|
|
||||||
} else {
|
|
||||||
end_point[0] = x_min[0];
|
|
||||||
end_point[1] = x_max[1];
|
|
||||||
res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, adjust_dir, tol_, t_max, t_out);
|
|
||||||
if (res_s_on_curve(t_out) <= 0) {
|
|
||||||
t_max = t_out;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
// Note: In some experiments modifiedRegularFalsi does not yield a result under the given tolerance.
|
|
||||||
t = RootFinder::solve(res_s_on_curve, 0., t_max, maxit_, tol_, iters_used_falsi);
|
|
||||||
res_s_on_curve.curve.computeXOfT(x, t);
|
|
||||||
} else {
|
|
||||||
if (res[1] < 0) {
|
|
||||||
end_point[0] = x_max[0];
|
|
||||||
end_point[1] = x_max[1];
|
|
||||||
res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, adjust_dir, tol_, t_max, t_out);
|
|
||||||
if (res_c_on_curve(t_out) >= 0) {
|
|
||||||
t_max = t_out;
|
|
||||||
}
|
|
||||||
} else {
|
|
||||||
end_point[0] = x_min[0];
|
|
||||||
end_point[1] = x_min[1];
|
|
||||||
res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, adjust_dir, tol_, t_max, t_out);
|
|
||||||
if (res_c_on_curve(t_out) <= 0) {
|
|
||||||
t_max = t_out;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
t = RootFinder::solve(res_c_on_curve, 0., t_max, maxit_, tol_, iters_used_falsi);
|
|
||||||
res_c_on_curve.curve.computeXOfT(x, t);
|
|
||||||
|
|
||||||
}
|
|
||||||
res_eq.computeResidual(x, res, mc, ff);
|
|
||||||
iters_used_split += 1;
|
|
||||||
}
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
if ((iters_used_split >= max_iters_split) && (norm(res) > tol_)) {
|
|
||||||
MESSAGE("Newton for single cell did not work in cell number " << cell);
|
if ((iters_used_split >= max_iters_split) && (norm(res) > tol_)) {
|
||||||
solveSingleCellBracketing(cell);
|
MESSAGE("Newton for single cell did not work in cell number " << cell);
|
||||||
} else {
|
solveSingleCellBracketing(cell);
|
||||||
concentration_[cell] = x[1];
|
} else {
|
||||||
cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
|
concentration_[cell] = x[1];
|
||||||
saturation_[cell] = x[0];
|
cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
|
||||||
fractionalflow_[cell] = ff;
|
saturation_[cell] = x[0];
|
||||||
mc_[cell] = mc;
|
fractionalflow_[cell] = ff;
|
||||||
}
|
mc_[cell] = mc;
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
void TransportModelPolymer::solveSingleCellNewton(int cell)
|
void TransportModelPolymer::solveSingleCellNewton(int cell)
|
||||||
{
|
{
|
||||||
const int max_iters_split = maxit_;
|
const int max_iters_split = maxit_;
|
||||||
int iters_used_split = 0;
|
int iters_used_split = 0;
|
||||||
|
|
||||||
// Check if current state is an acceptable solution.
|
// Check if current state is an acceptable solution.
|
||||||
@@ -883,19 +805,6 @@ namespace Opm
|
|||||||
mc_[cell] = mc;
|
mc_[cell] = mc;
|
||||||
return;
|
return;
|
||||||
} else {
|
} else {
|
||||||
/*
|
|
||||||
x[0] = saturation_[]-res[0];
|
|
||||||
if((x[0]>1) || (x[0]<0)){
|
|
||||||
x[0] = 0.5;
|
|
||||||
}
|
|
||||||
if(x[0]>0){
|
|
||||||
x[1] = concentration_[cell]*saturation_[cell]-res[1];
|
|
||||||
x[1] = x[1]/x[0];
|
|
||||||
}else{
|
|
||||||
x[1]=0;
|
|
||||||
}
|
|
||||||
*/
|
|
||||||
|
|
||||||
x[0]=0.5; x[1]=polyprops_.cMax()/2.0;
|
x[0]=0.5; x[1]=polyprops_.cMax()/2.0;
|
||||||
res_eq.computeResidual(x, res, mc, ff);
|
res_eq.computeResidual(x, res, mc, ff);
|
||||||
}
|
}
|
||||||
@@ -914,7 +823,7 @@ namespace Opm
|
|||||||
double dres_s_dsdc[2];
|
double dres_s_dsdc[2];
|
||||||
double dres_c_dsdc[2];
|
double dres_c_dsdc[2];
|
||||||
res_eq.computeJacobiRes(x, dres_s_dsdc, dres_c_dsdc);
|
res_eq.computeJacobiRes(x, dres_s_dsdc, dres_c_dsdc);
|
||||||
// double det = dres_s_dsdc[0]*dres_c_dsdc[1] - dres_c_dsdc[0]*dres_s_dsdc[1];
|
// double det = dres_s_dsdc[0]*dres_c_dsdc[1] - dres_c_dsdc[0]*dres_s_dsdc[1];
|
||||||
double dFx_dx=(dres_s_dsdc[0]-x[1]*dres_s_dsdc[1]/x[0]);
|
double dFx_dx=(dres_s_dsdc[0]-x[1]*dres_s_dsdc[1]/x[0]);
|
||||||
double dFx_dy= (dres_s_dsdc[1]/x[0]);
|
double dFx_dy= (dres_s_dsdc[1]/x[0]);
|
||||||
double dFy_dx=(dres_c_dsdc[0]-x[1]*dres_c_dsdc[1]/x[0]);
|
double dFy_dx=(dres_c_dsdc[0]-x[1]*dres_c_dsdc[1]/x[0]);
|
||||||
@@ -1390,7 +1299,7 @@ namespace
|
|||||||
// rectangle. x_out=(s_out, c_out) denotes the values of s and c at that point.
|
// rectangle. x_out=(s_out, c_out) denotes the values of s and c at that point.
|
||||||
void CurveInSCPlane::setup(const double* x, const double* direction,
|
void CurveInSCPlane::setup(const double* x, const double* direction,
|
||||||
const double* end_point, const double* x_min,
|
const double* end_point, const double* x_min,
|
||||||
const double* x_max, const bool adjust_dir, const double tol,
|
const double* x_max, const double tol,
|
||||||
double& t_max_out, double& t_out_out)
|
double& t_max_out, double& t_out_out)
|
||||||
{
|
{
|
||||||
x_[0] = x[0];
|
x_[0] = x[0];
|
||||||
@@ -1407,12 +1316,9 @@ namespace
|
|||||||
if (size_direction < tol) {
|
if (size_direction < tol) {
|
||||||
direction_[0] = end_point_[0]-x_[0];
|
direction_[0] = end_point_[0]-x_[0];
|
||||||
direction_[1] = end_point_[1]-x_[1];
|
direction_[1] = end_point_[1]-x_[1];
|
||||||
}
|
} else if ((end_point_[0]-x_[0])*direction_[0] + (end_point_[1]-x_[1])*direction_[1] < 0) {
|
||||||
if (adjust_dir || size_direction < 1e-5) {
|
direction_[0] *= -1.0;
|
||||||
if ((end_point_[0]-x_[0])*direction_[0] + (end_point_[1]-x_[1])*direction_[1] < 0) {
|
direction_[1] *= -1.0;
|
||||||
direction_[0] *= -1.0;
|
|
||||||
direction_[1] *= -1.0;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
bool t0_exists = true;
|
bool t0_exists = true;
|
||||||
double t0 = 0; // dummy default value (so that compiler does not complain).
|
double t0 = 0; // dummy default value (so that compiler does not complain).
|
||||||
@@ -1427,7 +1333,7 @@ namespace
|
|||||||
double t1 = 0; // dummy default value.
|
double t1 = 0; // dummy default value.
|
||||||
if (direction_[1] > 0) {
|
if (direction_[1] > 0) {
|
||||||
t1 = (x_max_[1] - x_[1])/direction_[1];
|
t1 = (x_max_[1] - x_[1])/direction_[1];
|
||||||
} else if (direction[1] < 0) {
|
} else if (direction_[1] < 0) {
|
||||||
t1 = (x_min_[1] - x_[1])/direction_[1];
|
t1 = (x_min_[1] - x_[1])/direction_[1];
|
||||||
} else {
|
} else {
|
||||||
t1_exists = false;
|
t1_exists = false;
|
||||||
@@ -1487,20 +1393,34 @@ namespace
|
|||||||
return res_eq_.computeResidualC(x_of_t);
|
return res_eq_.computeResidualC(x_of_t);
|
||||||
}
|
}
|
||||||
|
|
||||||
bool solveNewtonStep(const double* x, const Opm::TransportModelPolymer::ResidualEquation& res_eq,
|
bool solveNewtonStep(const double* xx, const Opm::TransportModelPolymer::ResidualEquation& res_eq,
|
||||||
const double* res, double* x_new) {
|
const double* res, double* x_new) {
|
||||||
|
|
||||||
double dres_s_dsdc[2];
|
double dres_s_dsdc[2];
|
||||||
double dres_c_dsdc[2];
|
double dres_c_dsdc[2];
|
||||||
|
double dx=1e-11;
|
||||||
|
double x[2];
|
||||||
|
if(!(x[0]>0)){
|
||||||
|
x[0] = dx;
|
||||||
|
x[1] = 0;
|
||||||
|
} else {
|
||||||
|
x[0] = xx[0];
|
||||||
|
x[1] = xx[1];
|
||||||
|
}
|
||||||
res_eq.computeJacobiRes(x, dres_s_dsdc, dres_c_dsdc);
|
res_eq.computeJacobiRes(x, dres_s_dsdc, dres_c_dsdc);
|
||||||
|
|
||||||
double det = dres_s_dsdc[0]*dres_c_dsdc[1] - dres_c_dsdc[0]*dres_s_dsdc[1];
|
double dFx_dx=(dres_s_dsdc[0]-x[1]*dres_s_dsdc[1]);
|
||||||
|
double dFx_dy= (dres_s_dsdc[1]/x[0]);
|
||||||
|
double dFy_dx=(dres_c_dsdc[0]-x[1]*dres_c_dsdc[1]);
|
||||||
|
double dFy_dy= (dres_c_dsdc[1]/x[0]);
|
||||||
|
|
||||||
|
double det = dFx_dx*dFy_dy - dFy_dx*dFx_dy;
|
||||||
if (std::abs(det) < 1e-8) {
|
if (std::abs(det) < 1e-8) {
|
||||||
return false;
|
return false;
|
||||||
} else {
|
} else {
|
||||||
x_new[0] = x[0] - (res[0]*dres_c_dsdc[1] - res[1]*dres_s_dsdc[1])/det;
|
x_new[0] = x[0] - (res[0]*dFy_dy - res[1]*dFx_dy)/det;
|
||||||
x_new[1] = x[1] - (res[1]*dres_s_dsdc[0] - res[0]*dres_c_dsdc[0])/det;
|
x_new[1] = x[1]*x[0] - (res[1]*dFx_dx - res[0]*dFy_dx)/det;
|
||||||
|
x_new[1] = (x_new[0]>0) ? x_new[1]/x_new[0] : 0.0;
|
||||||
return true;
|
return true;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|||||||
@@ -39,7 +39,7 @@ namespace Opm
|
|||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
|
|
||||||
enum SingleCellMethod { Bracketing, Newton };
|
enum SingleCellMethod { Bracketing, Newton, Gradient };
|
||||||
enum GradientMethod { Analytic, FinDif }; // Analytic is chosen (hard-coded)
|
enum GradientMethod { Analytic, FinDif }; // Analytic is chosen (hard-coded)
|
||||||
|
|
||||||
/// Construct solver.
|
/// Construct solver.
|
||||||
@@ -104,7 +104,7 @@ namespace Opm
|
|||||||
virtual void solveMultiCell(const int num_cells, const int* cells);
|
virtual void solveMultiCell(const int num_cells, const int* cells);
|
||||||
void solveSingleCellBracketing(int cell);
|
void solveSingleCellBracketing(int cell);
|
||||||
void solveSingleCellNewton(int cell);
|
void solveSingleCellNewton(int cell);
|
||||||
void solveSingleCellNewtonGradient(int cell);
|
void solveSingleCellGradient(int cell);
|
||||||
class ResidualEquation;
|
class ResidualEquation;
|
||||||
|
|
||||||
void initGravity(const double* grav);
|
void initGravity(const double* grav);
|
||||||
|
|||||||
Reference in New Issue
Block a user