OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Added sc formulation for gradient method.

Browse Source
This commit is contained in:
Xavier Raynaud 2012-06-14 13:37:08 +02:00
parent 9df87db47f
commit f3eb14e528
2 changed files with 80 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

107
opm/polymer/TransportModelPolymer.cpp
View File

@ -685,11 +685,13 @@ 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], saturation_[cell]*concentration_[cell]};
double res[2]; double res[2];
double mc; double mc;
double ff; double ff;
res_eq.computeResidual(x, res, mc, ff); double x_c[2];
scToc(x, x_c);
res_eq.computeResidual(x_c, res, mc, ff);
if (norm(res) <= tol_) { if (norm(res) <= tol_) {
cmax_[cell] = std::max(cmax_[cell], concentration_[cell]); cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
fractionalflow_[cell] = ff; fractionalflow_[cell] = ff;
@ -699,12 +701,16 @@ namespace Opm
double x_min[2] = { 0.0, 0.0 }; double x_min[2] = { 0.0, 0.0 };
double x_max[2] = { 1.0, polyprops_.cMax()*1.1 }; double x_max[2] = { 1.0, polyprops_.cMax()*1.1 };
double x_min_res_s[2] = { x_min[0], x_min[1] };
double x_max_res_s[2] = { x_max[0], x_min[0] };
double x_min_res_sc[2] = { x_min[0], x_min[1] };
double x_max_res_sc[2] = { x_max[0], x_max[1] };
double t; double t;
double t_max; double t_max;
double t_out; double t_out;
double direction[2]; double direction[2];
double end_point[2]; double end_point[2];
double gradient[2]; double gradient[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;
@ -712,46 +718,75 @@ namespace Opm
while ((norm(res) > tol_) && (iters_used_split < max_iters_split)) { while ((norm(res) > tol_) && (iters_used_split < max_iters_split)) {
if (std::abs(res[0]) < std::abs(res[1])) { if (std::abs(res[0]) < std::abs(res[1])) {
if (res[0] < -tol_) { if (res[0] < -tol_) {
direction[0] = x_max[0] - x[0]; direction[0] = x_max_res_s[0] - x[0];
direction[1] = x_min[1] - x[1]; direction[1] = x_max_res_s[1] - x[1];
if_res_s = true; if_res_s = true;
} else if (res[0] > tol_) { } else if (res[0] > tol_) {
direction[0] = x_min[0] - x[0]; direction[0] = x_min_res_s[0] - x[0];
direction[1] = x_max[1] - x[1]; direction[1] = x_min_res_s[1] - x[1];
if_res_s = true; if_res_s = true;
} else { } else {
res_eq.computeGradientResS(x, res, gradient); scToc(x, x_c);
direction[0] = -gradient[1]; res_eq.computeGradientResS(x_c, res, gradient);
direction[1] = gradient[0]; // dResS/d(s_) = dResS/ds - c/s*dResS/ds
// dResS/d(sc_) = -1/s*dResS/dc
if (x[0] < 1e-2*tol_) {
// We take 1.0/s*gradient[1]: wrong for linear permeability,
// acceptable for nonlinear relative permeability.
direction[0] = 0.0;
direction[1] = gradient[0];
} else {
// With s,c variables, we should have
// direction[0] = -gradient[1];
// direction[1] = gradient[0];
// With s, sc variables, we get:
scToc(x, x_c);
direction[0] = 1.0/x[0]*gradient[1];
direction[1] = gradient[0] - x_c[1]/x[0]*gradient[1];
}
if_res_s = false; if_res_s = false;
} }
} else { } else {
if (res[1] < -tol_) { if (res[1] < -tol_) {
direction[0] = x_max[0] - x[0]; direction[0] = x_max_res_sc[0] - x[0];
direction[1] = x_max[1] - x[1]; direction[1] = x_max_res_sc[1] - x[1];
if_res_s = false; if_res_s = false;
} else if (res[1] > tol_) { } else if (res[1] > tol_) {
direction[0] = x_min[0] - x[0]; direction[0] = x_min_res_sc[0] - x[0];
direction[1] = x_min[1] - x[1]; direction[1] = x_min_res_sc[1] - x[1];
if_res_s = false; if_res_s = false;
} else { } else {
res_eq.computeGradientResC(x, res, gradient); res_eq.computeGradientResC(x, res, gradient);
direction[0] = -gradient[1]; // dResC/d(s_) = dResC/ds - c/s*dResC/ds
direction[1] = gradient[0]; // dResC/d(sc_) = -1/s*dResC/dc
if (x[0] < 1e-2*tol_) {
// We take 1.0/s*gradient[1]: wrong for linear permeability,
// acceptable for nonlinear relative permeability.
direction[0] = 0.0;
direction[1] = gradient[0];
} else {
// With s,c variables, we should have
// direction[0] = -gradient[1];
// direction[1] = gradient[0];
// With s, sc variables, we get:
scToc(x, x_c);
direction[0] = 1.0/x[0]*gradient[1];
direction[1] = gradient[0] - x_c[1]/x[0]*gradient[1];
}
if_res_s = true; if_res_s = true;
} }
} }
if (if_res_s) { if (if_res_s) {
if (res[0] < 0) { if (res[0] < 0) {
end_point[0] = x_max[0]; end_point[0] = x_max_res_s[0];
end_point[1] = x_min[1]; end_point[1] = x_max_res_s[1];
res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out); 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) { if (res_s_on_curve(t_out) >= 0) {
t_max = t_out; t_max = t_out;
} }
} else { } else {
end_point[0] = x_min[0]; end_point[0] = x_min_res_s[0];
end_point[1] = x_max[1]; end_point[1] = x_min_res_s[1];
res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out); 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) { if (res_s_on_curve(t_out) <= 0) {
t_max = t_out; t_max = t_out;
@ -762,15 +797,15 @@ namespace Opm
res_s_on_curve.curve.computeXOfT(x, t); res_s_on_curve.curve.computeXOfT(x, t);
} else { } else {
if (res[1] < 0) { if (res[1] < 0) {
end_point[0] = x_max[0]; end_point[0] = x_max_res_sc[0];
end_point[1] = x_max[1]; end_point[1] = x_max_res_sc[1];
res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out); 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) { if (res_c_on_curve(t_out) >= 0) {
t_max = t_out; t_max = t_out;
} }
} else { } else {
end_point[0] = x_min[0]; end_point[0] = x_min_res_sc[0];
end_point[1] = x_min[1]; end_point[1] = x_min_res_sc[1];
res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out); 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) { if (res_c_on_curve(t_out) <= 0) {
t_max = t_out; t_max = t_out;
@ -780,7 +815,8 @@ namespace Opm
res_c_on_curve.curve.computeXOfT(x, t); res_c_on_curve.curve.computeXOfT(x, t);
} }
res_eq.computeResidual(x, res, mc, ff); scToc(x, x_c);
res_eq.computeResidual(x_c, res, mc, ff);
iters_used_split += 1; iters_used_split += 1;
} }
@ -790,7 +826,8 @@ namespace Opm
MESSAGE("Newton for single cell did not work in cell number " << cell); MESSAGE("Newton for single cell did not work in cell number " << cell);
solveSingleCellBracketing(cell); solveSingleCellBracketing(cell);
} else { } else {
concentration_[cell] = x[1]; scToc(x, x_c);
concentration_[cell] = x_c[1];
cmax_[cell] = std::max(cmax_[cell], concentration_[cell]); cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
saturation_[cell] = x[0]; saturation_[cell] = x[0];
fractionalflow_[cell] = ff; fractionalflow_[cell] = ff;
@ -1276,6 +1313,15 @@ namespace Opm
toBothSat(saturation_, saturation); toBothSat(saturation_, saturation);
} }
void TransportModelPolymer::scToc(const double* x, double* x_c) const {
x_c[0] = x[0];
if (x[0] < 1e-2*tol_) {
x_c[1] = polyprops_.cMax();
} else {
x_c[1] = x[1]/x[0];
}
}
} // namespace Opm } // namespace Opm
@ -1388,8 +1434,10 @@ namespace
double ResSOnCurve::operator()(const double t) const double ResSOnCurve::operator()(const double t) const
{ {
double x_of_t[2]; double x_of_t[2];
double x_c[2];
curve.computeXOfT(x_of_t, t); curve.computeXOfT(x_of_t, t);
return res_eq_.computeResidualS(x_of_t); res_eq_.tm.scToc(x_of_t, x_c);
return res_eq_.computeResidualS(x_c);
} }
ResCOnCurve::ResCOnCurve(const Opm::TransportModelPolymer::ResidualEquation& res_eq) ResCOnCurve::ResCOnCurve(const Opm::TransportModelPolymer::ResidualEquation& res_eq)
@ -1400,8 +1448,10 @@ namespace
double ResCOnCurve::operator()(const double t) const double ResCOnCurve::operator()(const double t) const
{ {
double x_of_t[2]; double x_of_t[2];
double x_c[2];
curve.computeXOfT(x_of_t, t); curve.computeXOfT(x_of_t, t);
return res_eq_.computeResidualC(x_of_t); res_eq_.tm.scToc(x_of_t, x_c);
return res_eq_.computeResidualC(x_c);
} }
bool solveNewtonStep(const double* xx, const Opm::TransportModelPolymer::ResidualEquation& res_eq, bool solveNewtonStep(const double* xx, const Opm::TransportModelPolymer::ResidualEquation& res_eq,
@ -1436,6 +1486,7 @@ namespace
} }
} }
} // Anonymous namespace } // Anonymous namespace

1
opm/polymer/TransportModelPolymer.hpp
View File

@ -132,6 +132,7 @@ namespace Opm
std::list<Newton_Iter> res_counts; std::list<Newton_Iter> res_counts;
void scToc(const double* x, double* x_c) const;
private: private:
const UnstructuredGrid& grid_; const UnstructuredGrid& grid_;