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

Cleaned up ownership (private vs public) in class definition.

Browse Source
This commit is contained in:
Xavier Raynaud 2012-08-30 09:09:23 +02:00
parent eb52ee58c1
commit 92def1d514
2 changed files with 131 additions and 124 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

232
opm/polymer/TransportModelCompressiblePolymer.cpp
View File

@ -62,6 +62,10 @@ public:
void computeGradientResS(const double* x, double* res, double* gradient) const;
void computeGradientResC(const double* x, double* res, double* gradient) const;
void computeJacobiRes(const double* x, double* dres_s_dsdc, double* dres_c_dsdc) const;
bool solveNewtonStepSC(const double*, const double*, double*) const;
bool solveNewtonStepC(const double* , const double*, double*) const;
private:
void computeResAndJacobi(const double* x, const bool if_res_s, const bool if_res_c,
@ -115,13 +119,6 @@ namespace
return std::max(std::abs(res[0]), std::abs(res[1]));
}
bool solveNewtonStepSC(const double* , const Opm::TransportModelCompressiblePolymer::ResidualEquation&,
const double*, double*);
bool solveNewtonStepC(const double* , const Opm::TransportModelCompressiblePolymer::ResidualEquation&,
const double*, double*);
// Define a piecewise linear curve along which we will look for zero of the "s" or "r" residual.
// The curve starts at "x", goes along the direction "direction" until it hits the boundary of the box of
// admissible values for "s" and "x" (which is given by "[x_min[0], x_max[0]]x[x_min[1], x_max[1]]").
@ -146,31 +143,6 @@ namespace
double x_[2];
};
// Compute the "s" residual along the curve "curve" for a given residual equation "res_eq".
// The operator() is sent to a root solver.
class ResSOnCurve
{
public:
ResSOnCurve(const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq);
double operator()(const double t) const;
CurveInSCPlane curve;
private:
const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq_;
};
// Compute the "c" residual along the curve "curve" for a given residual equation "res_eq".
// The operator() is sent to a root solver.
class ResCOnCurve
{
public:
ResCOnCurve(const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq);
double operator()(const double t) const;
CurveInSCPlane curve;
private:
const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq_;
};
}
@ -577,6 +549,60 @@ namespace Opm
}
}
// Compute the "s" residual along the curve "curve" for a given residual equation "res_eq".
// The operator() is sent to a root solver.
class TransportModelCompressiblePolymer::ResSOnCurve
{
public:
ResSOnCurve(const TransportModelCompressiblePolymer::ResidualEquation& res_eq);
double operator()(const double t) const;
CurveInSCPlane curve;
private:
const TransportModelCompressiblePolymer::ResidualEquation& res_eq_;
};
// Compute the "c" residual along the curve "curve" for a given residual equation "res_eq".
// The operator() is sent to a root solver.
class TransportModelCompressiblePolymer::ResCOnCurve
{
public:
ResCOnCurve(const TransportModelCompressiblePolymer::ResidualEquation& res_eq);
double operator()(const double t) const;
CurveInSCPlane curve;
private:
const TransportModelCompressiblePolymer::ResidualEquation& res_eq_;
};
TransportModelCompressiblePolymer::ResSOnCurve::ResSOnCurve(const TransportModelCompressiblePolymer::ResidualEquation& res_eq)
: res_eq_(res_eq)
{
}
double TransportModelCompressiblePolymer::ResSOnCurve::operator()(const double t) const
{
double x_of_t[2];
double x_c[2];
curve.computeXOfT(x_of_t, t);
res_eq_.tm.scToc(x_of_t, x_c);
return res_eq_.computeResidualS(x_c);
}
TransportModelCompressiblePolymer::ResCOnCurve::ResCOnCurve(const TransportModelCompressiblePolymer::ResidualEquation& res_eq)
: res_eq_(res_eq)
{
}
double TransportModelCompressiblePolymer::ResCOnCurve::operator()(const double t) const
{
double x_of_t[2];
double x_c[2];
curve.computeXOfT(x_of_t, t);
res_eq_.tm.scToc(x_of_t, x_c);
return res_eq_.computeResidualC(x_c);
}
void TransportModelCompressiblePolymer::solveSingleCell(const int cell)
{
switch (method_) {
@ -1472,6 +1498,62 @@ namespace Opm
}
}
bool TransportModelCompressiblePolymer::ResidualEquation::solveNewtonStepSC(const double* xx,
const double* res, double* x_new) const {
double dres_s_dsdc[2];
double dres_c_dsdc[2];
double dx=1e-8;
double x[2];
if(!(x[0]>0)){
x[0] = dx;
x[1] = 0;
} else {
x[0] = xx[0];
x[1] = xx[1];
}
computeJacobiRes(x, dres_s_dsdc, dres_c_dsdc);
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) {
return false;
} else {
x_new[0] = x[0] - (res[0]*dFy_dy - res[1]*dFx_dy)/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;
}
}
bool TransportModelCompressiblePolymer::ResidualEquation::solveNewtonStepC(const double* xx, const double* res, double* x_new) const {
double dres_s_dsdc[2];
double dres_c_dsdc[2];
double dx=1e-8;
double x[2];
if(!(x[0]>0)){
x[0] = dx;
x[1] = 0;
} else {
x[0] = xx[0];
x[1] = xx[1];
}
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];
if (std::abs(det) < 1e-8) {
return false;
} else {
x_new[0] = x[0] - (res[0]*dres_c_dsdc[1] - res[1]*dres_s_dsdc[1])/det;
x_new[1] = x[1] - (res[1]*dres_s_dsdc[0] - res[0]*dres_c_dsdc[0])/det;
return true;
}
}
} // namespace Opm
@ -1575,92 +1657,6 @@ namespace
}
}
ResSOnCurve::ResSOnCurve(const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq)
: res_eq_(res_eq)
{
}
double ResSOnCurve::operator()(const double t) const
{
double x_of_t[2];
double x_c[2];
curve.computeXOfT(x_of_t, t);
res_eq_.tm.scToc(x_of_t, x_c);
return res_eq_.computeResidualS(x_c);
}
ResCOnCurve::ResCOnCurve(const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq)
: res_eq_(res_eq)
{
}
double ResCOnCurve::operator()(const double t) const
{
double x_of_t[2];
double x_c[2];
curve.computeXOfT(x_of_t, t);
res_eq_.tm.scToc(x_of_t, x_c);
return res_eq_.computeResidualC(x_c);
}
bool solveNewtonStepSC(const double* xx, const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq,
const double* res, double* x_new) {
double dres_s_dsdc[2];
double dres_c_dsdc[2];
double dx=1e-8;
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);
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) {
return false;
} else {
x_new[0] = x[0] - (res[0]*dFy_dy - res[1]*dFx_dy)/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;
}
}
bool solveNewtonStepC(const double* xx, const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq,
const double* res, double* x_new) {
double dres_s_dsdc[2];
double dres_c_dsdc[2];
double dx=1e-8;
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);
double det = dres_s_dsdc[0]*dres_c_dsdc[1] - dres_c_dsdc[0]*dres_s_dsdc[1];
if (std::abs(det) < 1e-8) {
return false;
} else {
x_new[0] = x[0] - (res[0]*dres_c_dsdc[1] - res[1]*dres_s_dsdc[1])/det;
x_new[1] = x[1] - (res[1]*dres_s_dsdc[0] - res[0]*dres_c_dsdc[0])/det;
return true;
}
}
} // Anonymous namespace

23
opm/polymer/TransportModelCompressiblePolymer.hpp
View File

@ -29,6 +29,11 @@
class UnstructuredGrid;
namespace {
class ResSOnCurve;
class ResCOnCurve;
}
namespace Opm
{
@ -109,11 +114,8 @@ namespace Opm
std::vector<double>& cmax);
// This should be private:
class ResidualEquation;
friend class TransportModelCompressiblePolymer::ResidualEquation;
void scToc(const double* x, double* x_c) const;
//
private:
@ -162,13 +164,21 @@ namespace Opm
std::vector<int> ia_downw_;
std::vector<int> ja_downw_;
struct ResidualC;
struct ResidualS;
class ResidualCGrav;
class ResidualSGrav;
class ResidualEquation;
class ResSOnCurve;
class ResCOnCurve;
friend class TransportModelCompressiblePolymer::ResidualEquation;
friend class TransportModelCompressiblePolymer::ResSOnCurve;
friend class TransportModelCompressiblePolymer::ResCOnCurve;
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
void solveSingleCellBracketing(int cell);
@ -191,6 +201,7 @@ namespace Opm
void computeMc(double c, double& mc) const;
void computeMcWithDer(double c, double& mc, double& dmc_dc) const;
void mobility(double s, double c, int cell, double* mob) const;
void scToc(const double* x, double* x_c) const;
};
} // namespace Opm