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 enum to choose gradient computation method. By default, analytic.

Browse Source
This commit is contained in:
Xavier Raynaud 2012-03-28 10:14:24 +02:00
parent 724645b3f8
commit 18e45caf4f
2 changed files with 14 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
opm/polymer/TransportModelPolymer.cpp
View File

@ -42,7 +42,7 @@ public:
double c_max_ads; double c_max_ads;
double rhor; double rhor;
double ads0; double ads0;
int gradient_method; GradientMethod gradient_method;
const TransportModelPolymer& tm; const TransportModelPolymer& tm;
ResidualEquation(const TransportModelPolymer& tmodel, int cell_index); ResidualEquation(const TransportModelPolymer& tmodel, int cell_index);
@ -338,7 +338,7 @@ namespace Opm
TransportModelPolymer::ResidualEquation::ResidualEquation(const TransportModelPolymer& tmodel, int cell_index) TransportModelPolymer::ResidualEquation::ResidualEquation(const TransportModelPolymer& tmodel, int cell_index)
: tm(tmodel) : tm(tmodel)
{ {
gradient_method = 1; gradient_method = Analytic;
cell = cell_index; cell = cell_index;
s0 = tm.saturation_[cell]; s0 = tm.saturation_[cell];
c0 = tm.concentration_[cell]; c0 = tm.concentration_[cell];
@ -433,10 +433,10 @@ namespace Opm
} }
void TransportModelPolymer::ResidualEquation::computeGradientResS(const double* x, double* res, double* gradient) const void TransportModelPolymer::ResidualEquation::computeGradientResS(const double* x, double* res, double* gradient) const
// If gradient_method == 1, use finite difference // If gradient_method == FinDif, use finite difference
// If gradient_method == 2, use analytic expresions // If gradient_method == Analytic, use analytic expresions
{ {
if (gradient_method == 1) { if (gradient_method == FinDif) {
double epsi = 1e-8; double epsi = 1e-8;
double res_epsi[2]; double res_epsi[2];
double x_epsi[2]; double x_epsi[2];
@ -449,7 +449,7 @@ namespace Opm
x_epsi[1] = x[1] + epsi; x_epsi[1] = x[1] + epsi;
computeResidual(x_epsi, res_epsi); computeResidual(x_epsi, res_epsi);
gradient[1] = (res_epsi[0] - res[0])/epsi; gradient[1] = (res_epsi[0] - res[0])/epsi;
} else if (gradient_method == 2) { } else if (gradient_method == Analytic) {
double s = x[0]; double s = x[0];
double c = x[1]; double c = x[1];
double ff_ds_dc[2]; double ff_ds_dc[2];
@ -468,10 +468,10 @@ namespace Opm
} }
void TransportModelPolymer::ResidualEquation::computeGradientResC(const double* x, double* res, double* gradient) const void TransportModelPolymer::ResidualEquation::computeGradientResC(const double* x, double* res, double* gradient) const
// If gradient_method == 1, use finite difference // If gradient_method == FinDif, use finite difference
// If gradient_method == 2, use analytic expresions // If gradient_method == Analytic, use analytic expresions
{ {
if (gradient_method == 1) { if (gradient_method == FinDif) {
double epsi = 1e-8; double epsi = 1e-8;
double res_epsi[2]; double res_epsi[2];
double x_epsi[2]; double x_epsi[2];
@ -484,7 +484,7 @@ namespace Opm
x_epsi[1] = x[1] + epsi; x_epsi[1] = x[1] + epsi;
computeResidual(x_epsi, res_epsi); computeResidual(x_epsi, res_epsi);
gradient[1] = (res_epsi[1] - res[1])/epsi; gradient[1] = (res_epsi[1] - res[1])/epsi;
} else if (gradient_method == 2) { } else if (gradient_method == Analytic) {
double s = x[0]; double s = x[0];
double c = x[1]; double c = x[1];
double ff_ds_dc[2]; double ff_ds_dc[2];
@ -509,7 +509,7 @@ namespace Opm
// Compute the Jacobian of the residual equations. // Compute the Jacobian of the residual equations.
void TransportModelPolymer::ResidualEquation::computeJacobiRes(const double* x, double* res_s_ds_dc, double* res_c_ds_dc) const void TransportModelPolymer::ResidualEquation::computeJacobiRes(const double* x, double* res_s_ds_dc, double* res_c_ds_dc) const
{ {
if (gradient_method == 1) { if (gradient_method == FinDif) {
double epsi = 1e-8; double epsi = 1e-8;
double res[2]; double res[2];
double res_epsi[2]; double res_epsi[2];
@ -531,7 +531,7 @@ namespace Opm
x_epsi[1] = x[1] + epsi; x_epsi[1] = x[1] + epsi;
computeResidual(x_epsi, res_epsi); computeResidual(x_epsi, res_epsi);
res_c_ds_dc[1] = (res_epsi[1] - res[1])/epsi; res_c_ds_dc[1] = (res_epsi[1] - res[1])/epsi;
} else if (gradient_method == 2) { } else if (gradient_method == Analytic) {
double s = x[0]; double s = x[0];
double c = x[1]; double c = x[1];
double ff_ds_dc[2]; double ff_ds_dc[2];
@ -869,7 +869,7 @@ namespace Opm
double mu_m_dc; // derivative of mu_m with respect to c double mu_m_dc; // derivative of mu_m with respect to c
double mu_m = polyprops_.viscMultWithDer(c, &mu_m_dc)*mu_w; double mu_m = polyprops_.viscMultWithDer(c, &mu_m_dc)*mu_w;
mu_m_dc *= mu_w; mu_m_dc *= mu_w;
double mu_p = polyprops_.viscMult(polyprops_.cMax())*mu_w; double mu_p = polyprops_.viscMult(c_max_limit)*mu_w;
double omega = polyprops_.mixParam(); double omega = polyprops_.mixParam();
double mu_w_e = std::pow(mu_m, omega)*std::pow(mu_w, 1 - omega); double mu_w_e = std::pow(mu_m, omega)*std::pow(mu_w, 1 - omega);
double mu_w_e_dc = omega*mu_m_dc*std::pow(mu_m, omega - 1)*std::pow(mu_w, 1 - omega); double mu_w_e_dc = omega*mu_m_dc*std::pow(mu_m, omega - 1)*std::pow(mu_w, 1 - omega);

1
opm/polymer/TransportModelPolymer.hpp
View File

@ -40,6 +40,7 @@ namespace Opm
public: public:
enum SingleCellMethod { Bracketing, Newton }; enum SingleCellMethod { Bracketing, Newton };
enum GradientMethod { Analytic, FinDif }; // Analytic is chosen (hard-coded)
/// \TODO document me, especially method. /// \TODO document me, especially method.
TransportModelPolymer(const UnstructuredGrid& grid, TransportModelPolymer(const UnstructuredGrid& grid,