From 92def1d51417363b0a0bc38d4c29b5db8a705b97 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Thu, 30 Aug 2012 09:09:23 +0200 Subject: [PATCH] Cleaned up ownership (private vs public) in class definition. --- .../TransportModelCompressiblePolymer.cpp | 232 +++++++++--------- .../TransportModelCompressiblePolymer.hpp | 23 +- 2 files changed, 131 insertions(+), 124 deletions(-) diff --git a/opm/polymer/TransportModelCompressiblePolymer.cpp b/opm/polymer/TransportModelCompressiblePolymer.cpp index ec56d9e99..9b91c5fa0 100644 --- a/opm/polymer/TransportModelCompressiblePolymer.cpp +++ b/opm/polymer/TransportModelCompressiblePolymer.cpp @@ -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 diff --git a/opm/polymer/TransportModelCompressiblePolymer.hpp b/opm/polymer/TransportModelCompressiblePolymer.hpp index 9d5dfdf1d..3b39c2b00 100644 --- a/opm/polymer/TransportModelCompressiblePolymer.hpp +++ b/opm/polymer/TransportModelCompressiblePolymer.hpp @@ -29,6 +29,11 @@ class UnstructuredGrid; +namespace { + class ResSOnCurve; + class ResCOnCurve; +} + namespace Opm { @@ -109,11 +114,8 @@ namespace Opm std::vector& 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 ia_downw_; std::vector 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