mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Introduced ResidualEquation class to compute value of residual and derivatives.
This commit is contained in:
Xavier Raynaud
parent
e11cea7432
commit
196c29522d
@ -24,14 +24,89 @@
|
|||||||
#include <opm/core/utility/RootFinders.hpp>
|
#include <opm/core/utility/RootFinders.hpp>
|
||||||
#include <cmath>
|
#include <cmath>
|
||||||
|
|
||||||
|
|
||||||
|
class Opm::TransportModelPolymer::ResidualEquation
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
int cell;
|
||||||
|
double s0;
|
||||||
|
double c0;
|
||||||
|
double cmax0;
|
||||||
|
double influx; // sum_j min(v_ij, 0)*f(s_j)
|
||||||
|
double influx_polymer; // sum_j min(v_ij, 0)*f(s_j)*mc(c_j)
|
||||||
|
double outflux; // sum_j max(v_ij, 0)
|
||||||
|
double porosity;
|
||||||
|
double dtpv; // dt/pv(i)
|
||||||
|
double dps;
|
||||||
|
double rhor;
|
||||||
|
double ads0;
|
||||||
|
int gradient_method;
|
||||||
|
const TransportModelPolymer& tm;
|
||||||
|
|
||||||
|
ResidualEquation(const TransportModelPolymer& tmodel, int cell_index);
|
||||||
|
void computeResidual(const double* x, double* res) const;
|
||||||
|
void computeResidual(const double* x, double* res, double& mc, double& ff) const;
|
||||||
|
double computeResidualS(const double* x) const;
|
||||||
|
double computeResidualC(const double* x) const;
|
||||||
|
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* res_s_ds_dc, double* res_c_ds_dc) const;
|
||||||
|
|
||||||
|
};
|
||||||
|
|
||||||
|
|
||||||
namespace
|
namespace
|
||||||
{
|
{
|
||||||
double norm(double* res)
|
double norm(double* res)
|
||||||
{
|
{
|
||||||
return std::max(std::abs(res[0]), std::abs(res[1]));
|
return std::max(std::abs(res[0]), std::abs(res[1]));
|
||||||
}
|
}
|
||||||
|
|
||||||
|
bool solveNewtonStep(const double* , const Opm::TransportModelPolymer::ResidualEquation& , const double*, double*);
|
||||||
|
|
||||||
|
class CurveInSCPlane{
|
||||||
|
public:
|
||||||
|
CurveInSCPlane();
|
||||||
|
void setup(const double* , const double* ,
|
||||||
|
const double* , const double* ,
|
||||||
|
const double* , double& ,
|
||||||
|
double& );
|
||||||
|
void computeXOfT(double*, const double) const;
|
||||||
|
|
||||||
|
private:
|
||||||
|
double direction_[2];
|
||||||
|
double end_point_[2];
|
||||||
|
double x_max_[2];
|
||||||
|
double x_min_[2];
|
||||||
|
double t_out_;
|
||||||
|
double t_max_; // t_max = t_out + 1
|
||||||
|
double x_out_[2];
|
||||||
|
double x_[2];
|
||||||
|
};
|
||||||
|
|
||||||
|
class ResSOnCurve
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
ResSOnCurve(const Opm::TransportModelPolymer::ResidualEquation& res_eq);
|
||||||
|
double operator()(const double t) const;
|
||||||
|
CurveInSCPlane curve;
|
||||||
|
private:
|
||||||
|
Opm::TransportModelPolymer::ResidualEquation res_eq_;
|
||||||
|
};
|
||||||
|
|
||||||
|
class ResCOnCurve
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
ResCOnCurve(const Opm::TransportModelPolymer::ResidualEquation& res_eq);
|
||||||
|
double operator()(const double t) const;
|
||||||
|
CurveInSCPlane curve;
|
||||||
|
private:
|
||||||
|
Opm::TransportModelPolymer::ResidualEquation res_eq_;
|
||||||
|
};
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
namespace Opm
|
namespace Opm
|
||||||
{
|
{
|
||||||
TransportModelPolymer::TransportModelPolymer(const UnstructuredGrid& grid,
|
TransportModelPolymer::TransportModelPolymer(const UnstructuredGrid& grid,
|
||||||
@ -241,334 +316,236 @@ namespace Opm
|
|||||||
};
|
};
|
||||||
|
|
||||||
|
|
||||||
// Residual for s and c. Includes method to compute the gradient.
|
// ResidualEquation gathers parameter to construct residual equations and to compute the value
|
||||||
// Compute residual in s (or c) for a given piecewise linear curve (with only one node) in the s-c
|
// of the residual and of its derivative.
|
||||||
// plane. The method operator() is used by a 1d root solver.
|
|
||||||
|
|
||||||
struct TransportModelPolymer::Residual
|
|
||||||
|
TransportModelPolymer::ResidualEquation::ResidualEquation(const TransportModelPolymer& tmodel, int cell_index)
|
||||||
|
: tm(tmodel)
|
||||||
{
|
{
|
||||||
int cell;
|
gradient_method = 1;
|
||||||
double s0;
|
cell = cell_index;
|
||||||
double c0;
|
s0 = tm.saturation_[cell];
|
||||||
double cmax0;
|
c0 = tm.concentration_[cell];
|
||||||
double influx; // sum_j min(v_ij, 0)*f(s_j)
|
cmax0 = tm.cmax_[cell];
|
||||||
double influx_polymer; // sum_j min(v_ij, 0)*f(s_j)*mc(c_j)
|
dps = tm.polyprops_.deadPoreVol();
|
||||||
double outflux; // sum_j max(v_ij, 0)
|
rhor = tm.polyprops_.rockDensity();
|
||||||
double porosity;
|
ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
|
||||||
double dtpv; // dt/pv(i)
|
|
||||||
double direction[2];
|
|
||||||
double end_point[2];
|
|
||||||
double x_max[2];
|
|
||||||
double x_min[2];
|
|
||||||
double t_out;
|
|
||||||
double t_max; // t_max = t_out + 1
|
|
||||||
double x_out[2];
|
|
||||||
double x[2];
|
|
||||||
bool if_res_s;
|
|
||||||
const TransportModelPolymer& tm;
|
|
||||||
|
|
||||||
Residual(const TransportModelPolymer& tmodel, int cell_index)
|
double dflux = -tm.source_[cell];
|
||||||
: tm(tmodel)
|
bool src_is_inflow = dflux < 0.0;
|
||||||
{
|
influx = src_is_inflow ? dflux : 0.0;
|
||||||
cell = cell_index;
|
influx_polymer = src_is_inflow ? dflux*tm.computeMc(tm.inflow_c_) : 0.0;
|
||||||
s0 = tm.saturation_[cell];
|
outflux = !src_is_inflow ? dflux : 0.0;
|
||||||
c0 = tm.concentration_[cell];
|
dtpv = tm.dt_/tm.porevolume_[cell];
|
||||||
cmax0 = tm.cmax_[cell];
|
porosity = tm.porosity_[cell];
|
||||||
double dflux = -tm.source_[cell];
|
for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
|
||||||
bool src_is_inflow = dflux < 0.0;
|
int f = tm.grid_.cell_faces[i];
|
||||||
influx = src_is_inflow ? dflux : 0.0;
|
double flux;
|
||||||
influx_polymer = src_is_inflow ? dflux*tm.computeMc(tm.inflow_c_) : 0.0;
|
int other;
|
||||||
outflux = !src_is_inflow ? dflux : 0.0;
|
// Compute cell flux
|
||||||
dtpv = tm.dt_/tm.porevolume_[cell];
|
if (cell == tm.grid_.face_cells[2*f]) {
|
||||||
porosity = tm.porosity_[cell];
|
flux = tm.darcyflux_[f];
|
||||||
for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
|
other = tm.grid_.face_cells[2*f+1];
|
||||||
int f = tm.grid_.cell_faces[i];
|
} else {
|
||||||
double flux;
|
flux =-tm.darcyflux_[f];
|
||||||
int other;
|
other = tm.grid_.face_cells[2*f];
|
||||||
// Compute cell flux
|
}
|
||||||
if (cell == tm.grid_.face_cells[2*f]) {
|
// Add flux to influx or outflux, if interior.
|
||||||
flux = tm.darcyflux_[f];
|
if (other != -1) {
|
||||||
other = tm.grid_.face_cells[2*f+1];
|
if (flux < 0.0) {
|
||||||
|
influx += flux*tm.fractionalflow_[other];
|
||||||
|
influx_polymer += flux*tm.fractionalflow_[other]*tm.mc_[other];
|
||||||
} else {
|
} else {
|
||||||
flux =-tm.darcyflux_[f];
|
outflux += flux;
|
||||||
other = tm.grid_.face_cells[2*f];
|
|
||||||
}
|
|
||||||
// Add flux to influx or outflux, if interior.
|
|
||||||
if (other != -1) {
|
|
||||||
if (flux < 0.0) {
|
|
||||||
influx += flux*tm.fractionalflow_[other];
|
|
||||||
influx_polymer += flux*tm.fractionalflow_[other]*tm.mc_[other];
|
|
||||||
} else {
|
|
||||||
outflux += flux;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
}
|
||||||
|
|
||||||
void computeResidual(const double* x, double* res) const
|
void TransportModelPolymer::ResidualEquation::computeResidual(const double* x, double* res) const
|
||||||
{
|
{
|
||||||
|
double s = x[0];
|
||||||
|
double c = x[1];
|
||||||
|
double ff = tm.fracFlow(s, c, cell);
|
||||||
|
double mc = tm.computeMc(c);
|
||||||
|
double dps = tm.polyprops_.deadPoreVol();
|
||||||
|
double rhor = tm.polyprops_.rockDensity();
|
||||||
|
double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
|
||||||
|
double ads = tm.polyprops_.adsorbtion(std::max(c, cmax0));
|
||||||
|
res[0] = s - s0 + dtpv*(outflux*ff + influx);
|
||||||
|
res[1] = (s - dps)*c - (s0 - dps)*c0
|
||||||
|
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
|
||||||
|
+ dtpv*(outflux*ff*mc + influx_polymer);
|
||||||
|
}
|
||||||
|
|
||||||
|
void TransportModelPolymer::ResidualEquation::computeResidual(const double* x, double* res, double& mc, double& ff) const
|
||||||
|
{
|
||||||
|
double s = x[0];
|
||||||
|
double c = x[1];
|
||||||
|
ff = tm.fracFlow(s, c, cell);
|
||||||
|
mc = tm.computeMc(c);
|
||||||
|
double dps = tm.polyprops_.deadPoreVol();
|
||||||
|
double rhor = tm.polyprops_.rockDensity();
|
||||||
|
double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
|
||||||
|
double ads = tm.polyprops_.adsorbtion(std::max(c, cmax0));
|
||||||
|
res[0] = s - s0 + dtpv*(outflux*ff + influx);
|
||||||
|
res[1] = (s - dps)*c - (s0 - dps)*c0
|
||||||
|
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
|
||||||
|
+ dtpv*(outflux*ff*mc + influx_polymer);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
double TransportModelPolymer::ResidualEquation::computeResidualS(const double* x) const
|
||||||
|
{
|
||||||
|
double s = x[0];
|
||||||
|
double c = x[1];
|
||||||
|
double ff = tm.fracFlow(s, c, cell);
|
||||||
|
return s - s0 + dtpv*(outflux*ff + influx);
|
||||||
|
}
|
||||||
|
|
||||||
|
double TransportModelPolymer::ResidualEquation::computeResidualC(const double* x) const
|
||||||
|
{
|
||||||
|
double s = x[0];
|
||||||
|
double c = x[1];
|
||||||
|
double ff = tm.fracFlow(s, c, cell);
|
||||||
|
double mc = tm.computeMc(c);
|
||||||
|
double ads = tm.polyprops_.adsorbtion(std::max(c, cmax0));
|
||||||
|
return (s - dps)*c - (s0 - dps)*c0
|
||||||
|
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
|
||||||
|
+ dtpv*(outflux*ff*mc + influx_polymer);
|
||||||
|
}
|
||||||
|
|
||||||
|
void TransportModelPolymer::ResidualEquation::computeGradientResS(const double* x, double* res, double* gradient) const
|
||||||
|
// If gradient_method == 1, use finite difference
|
||||||
|
// If gradient_method == 2, use analytic expresions
|
||||||
|
{
|
||||||
|
if (gradient_method == 1) {
|
||||||
|
double epsi = 1e-8;
|
||||||
|
double res_epsi[2];
|
||||||
|
double x_epsi[2];
|
||||||
|
computeResidual(x, res);
|
||||||
|
x_epsi[0] = x[0] + epsi;
|
||||||
|
x_epsi[1] = x[1];
|
||||||
|
computeResidual(x_epsi, res_epsi);
|
||||||
|
gradient[0] = (res_epsi[0] - res[0])/epsi;
|
||||||
|
x_epsi[0] = x[0];
|
||||||
|
x_epsi[1] = x[1] + epsi;
|
||||||
|
computeResidual(x_epsi, res_epsi);
|
||||||
|
gradient[1] = (res_epsi[0] - res[0])/epsi;
|
||||||
|
} else if (gradient_method == 2) {
|
||||||
double s = x[0];
|
double s = x[0];
|
||||||
double c = x[1];
|
double c = x[1];
|
||||||
double ff = tm.fracFlow(s, c, cell);
|
double ff_ds_dc[2];
|
||||||
double mc = tm.computeMc(c);
|
double ff = tm.fracFlowWithDer(s, c, cell, ff_ds_dc);
|
||||||
double dps = tm.polyprops_.deadPoreVol();
|
double mc_dc;
|
||||||
double rhor = tm.polyprops_.rockDensity();
|
double mc = tm.computeMcWithDer(c, &mc_dc);
|
||||||
double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
|
double ads;
|
||||||
double ads = tm.polyprops_.adsorbtion(std::max(c, cmax0));
|
double ads_dc;
|
||||||
|
if (c < cmax0) {
|
||||||
|
ads = tm.polyprops_.adsorbtion(cmax0);
|
||||||
|
ads_dc = 0;
|
||||||
|
} else {
|
||||||
|
ads = tm.polyprops_.adsorbtionWithDer(c, &ads_dc);
|
||||||
|
}
|
||||||
res[0] = s - s0 + dtpv*(outflux*ff + influx);
|
res[0] = s - s0 + dtpv*(outflux*ff + influx);
|
||||||
res[1] = (s - dps)*c - (s0 - dps)*c0
|
res[1] = (s - dps)*c - (s0 - dps)*c0
|
||||||
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
|
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
|
||||||
+ dtpv*(outflux*ff*mc + influx_polymer);
|
+ dtpv*(outflux*ff*mc + influx_polymer);
|
||||||
|
gradient[0] = 1 + dtpv*outflux*ff_ds_dc[0];
|
||||||
|
gradient[1] = dtpv*outflux*ff_ds_dc[1];
|
||||||
}
|
}
|
||||||
|
}
|
||||||
|
|
||||||
void computeResidual(const double* x, double* res, double& mc, double& ff) const
|
void TransportModelPolymer::ResidualEquation::computeGradientResC(const double* x, double* res, double* gradient) const
|
||||||
{
|
// If gradient_method == 1, use finite difference
|
||||||
|
// If gradient_method == 2, use analytic expresions
|
||||||
|
{
|
||||||
|
if (gradient_method == 1) {
|
||||||
|
double epsi = 1e-8;
|
||||||
|
double res_epsi[2];
|
||||||
|
double x_epsi[2];
|
||||||
|
computeResidual(x, res);
|
||||||
|
x_epsi[0] = x[0] + epsi;
|
||||||
|
x_epsi[1] = x[1];
|
||||||
|
computeResidual(x_epsi, res_epsi);
|
||||||
|
gradient[0] = (res_epsi[1] - res[1])/epsi;
|
||||||
|
x_epsi[0] = x[0];
|
||||||
|
x_epsi[1] = x[1] + epsi;
|
||||||
|
computeResidual(x_epsi, res_epsi);
|
||||||
|
gradient[1] = (res_epsi[1] - res[1])/epsi;
|
||||||
|
} else if (gradient_method == 2) {
|
||||||
double s = x[0];
|
double s = x[0];
|
||||||
double c = x[1];
|
double c = x[1];
|
||||||
ff = tm.fracFlow(s, c, cell);
|
double ff_ds_dc[2];
|
||||||
mc = tm.computeMc(c);
|
double ff = tm.fracFlowWithDer(s, c, cell, ff_ds_dc);
|
||||||
|
double mc_dc;
|
||||||
|
double mc = tm.computeMcWithDer(c, &mc_dc);
|
||||||
double dps = tm.polyprops_.deadPoreVol();
|
double dps = tm.polyprops_.deadPoreVol();
|
||||||
double rhor = tm.polyprops_.rockDensity();
|
double rhor = tm.polyprops_.rockDensity();
|
||||||
double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
|
double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
|
||||||
double ads = tm.polyprops_.adsorbtion(std::max(c, cmax0));
|
double ads;
|
||||||
|
double ads_dc;
|
||||||
|
if (c < cmax0) {
|
||||||
|
ads = tm.polyprops_.adsorbtion(cmax0);
|
||||||
|
ads_dc = 0;
|
||||||
|
} else {
|
||||||
|
ads = tm.polyprops_.adsorbtionWithDer(c, &ads_dc);
|
||||||
|
}
|
||||||
res[0] = s - s0 + dtpv*(outflux*ff + influx);
|
res[0] = s - s0 + dtpv*(outflux*ff + influx);
|
||||||
res[1] = (s - dps)*c - (s0 - dps)*c0
|
res[1] = (s - dps)*c - (s0 - dps)*c0
|
||||||
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
|
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
|
||||||
+ dtpv*(outflux*ff*mc + influx_polymer);
|
+ dtpv*(outflux*ff*mc + influx_polymer);
|
||||||
|
gradient[0] = c + dtpv*outflux*(ff_ds_dc[0])*mc;
|
||||||
|
gradient[1] = s - dps + rhor*((1.0 - porosity)/porosity)*ads_dc
|
||||||
|
+ dtpv*outflux*(ff_ds_dc[1]*mc + ff*mc_dc);
|
||||||
}
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void TransportModelPolymer::ResidualEquation::computeJacobiRes(const double* x, double* res_s_ds_dc, double* res_c_ds_dc) const
|
||||||
bool solveNewtonStep(const double* x, double* x_new, const int& gradient_method) {
|
{
|
||||||
|
if (gradient_method == 1) {
|
||||||
// If gradient_method == 1, use finite difference
|
double epsi = 1e-8;
|
||||||
// If gradient_method == 2, use analytic expresions
|
|
||||||
|
|
||||||
double res[2];
|
double res[2];
|
||||||
double res_s_ds_dc[2];
|
double res_epsi[2];
|
||||||
double res_c_ds_dc[2];
|
double x_epsi[2];
|
||||||
|
computeResidual(x, res);
|
||||||
if (gradient_method == 1) {
|
x_epsi[0] = x[0] + epsi;
|
||||||
double epsi = 1e-8;
|
x_epsi[1] = x[1];
|
||||||
double res_epsi[2];
|
computeResidual(x_epsi, res_epsi);
|
||||||
double x_epsi[2];
|
res_s_ds_dc[0] = (res_epsi[0] - res[0])/epsi;
|
||||||
computeResidual(x, res);
|
x_epsi[0] = x[0];
|
||||||
x_epsi[0] = x[0] + epsi;
|
x_epsi[1] = x[1] + epsi;
|
||||||
x_epsi[1] = x[1];
|
computeResidual(x_epsi, res_epsi);
|
||||||
computeResidual(x_epsi, res_epsi);
|
res_s_ds_dc[1] = (res_epsi[0] - res[0])/epsi;
|
||||||
res_s_ds_dc[0] = (res_epsi[0] - res[0])/epsi;
|
x_epsi[0] = x[0] + epsi;
|
||||||
x_epsi[0] = x[0];
|
x_epsi[1] = x[1];
|
||||||
x_epsi[1] = x[1] + epsi;
|
computeResidual(x_epsi, res_epsi);
|
||||||
computeResidual(x_epsi, res_epsi);
|
res_c_ds_dc[0] = (res_epsi[1] - res[1])/epsi;
|
||||||
res_s_ds_dc[1] = (res_epsi[0] - res[0])/epsi;
|
x_epsi[0] = x[0];
|
||||||
x_epsi[0] = x[0] + epsi;
|
x_epsi[1] = x[1] + epsi;
|
||||||
x_epsi[1] = x[1];
|
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[0] = (res_epsi[1] - res[1])/epsi;
|
} else if (gradient_method == 2) {
|
||||||
x_epsi[0] = x[0];
|
double s = x[0];
|
||||||
x_epsi[1] = x[1] + epsi;
|
double c = x[1];
|
||||||
computeResidual(x_epsi, res_epsi);
|
double ff_ds_dc[2];
|
||||||
res_c_ds_dc[1] = (res_epsi[1] - res[1])/epsi;
|
double ff = tm.fracFlowWithDer(s, c, cell, ff_ds_dc);
|
||||||
} else if (gradient_method == 2) {
|
double mc_dc;
|
||||||
double s = x[0];
|
double mc = tm.computeMcWithDer(c, &mc_dc);
|
||||||
double c = x[1];
|
double ads_dc;
|
||||||
double ff_ds_dc[2];
|
if (c < cmax0) {
|
||||||
double ff = tm.fracFlowWithDer(s, c, cell, ff_ds_dc);
|
ads_dc = 0;
|
||||||
double mc_dc;
|
|
||||||
double mc = tm.computeMcWithDer(c, &mc_dc);
|
|
||||||
double dps = tm.polyprops_.deadPoreVol();
|
|
||||||
double rhor = tm.polyprops_.rockDensity();
|
|
||||||
double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
|
|
||||||
double ads;
|
|
||||||
double ads_dc;
|
|
||||||
if (c < cmax0) {
|
|
||||||
ads = tm.polyprops_.adsorbtion(cmax0);
|
|
||||||
ads_dc = 0;
|
|
||||||
} else {
|
|
||||||
ads = tm.polyprops_.adsorbtionWithDer(c, &ads_dc);
|
|
||||||
}
|
|
||||||
res[0] = s - s0 + dtpv*(outflux*ff + influx);
|
|
||||||
res[1] = (s - dps)*c - (s0 - dps)*c0
|
|
||||||
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
|
|
||||||
+ dtpv*(outflux*ff*mc + influx_polymer);
|
|
||||||
res_s_ds_dc[0] = 1 + dtpv*outflux*ff_ds_dc[0];
|
|
||||||
res_s_ds_dc[1] = dtpv*outflux*ff_ds_dc[1];
|
|
||||||
res_c_ds_dc[0] = c + dtpv*outflux*(ff_ds_dc[0])*mc;
|
|
||||||
res_c_ds_dc[1] = s - dps + rhor*((1.0 - porosity)/porosity)*ads_dc
|
|
||||||
+ dtpv*outflux*(ff_ds_dc[1]*mc + ff*mc_dc);
|
|
||||||
}
|
|
||||||
|
|
||||||
double det = res_s_ds_dc[0]*res_c_ds_dc[1] - res_c_ds_dc[0]*res_s_ds_dc[1];
|
|
||||||
if (std::abs(det) < 1e-8) {
|
|
||||||
return false;
|
|
||||||
} else {
|
} else {
|
||||||
x_new[0] = x[0] - (res[0]*res_c_ds_dc[1] - res[1]*res_s_ds_dc[1])/det;
|
tm.polyprops_.adsorbtionWithDer(c, &ads_dc);
|
||||||
x_new[1] = x[1] - (res[1]*res_s_ds_dc[0] - res[0]*res_c_ds_dc[0])/det;
|
|
||||||
return true;
|
|
||||||
}
|
}
|
||||||
|
res_s_ds_dc[0] = 1 + dtpv*outflux*ff_ds_dc[0];
|
||||||
|
res_s_ds_dc[1] = dtpv*outflux*ff_ds_dc[1];
|
||||||
|
res_c_ds_dc[0] = c + dtpv*outflux*(ff_ds_dc[0])*mc;
|
||||||
|
res_c_ds_dc[1] = s - dps + rhor*((1.0 - porosity)/porosity)*ads_dc
|
||||||
|
+ dtpv*outflux*(ff_ds_dc[1]*mc + ff*mc_dc);
|
||||||
}
|
}
|
||||||
|
}
|
||||||
void computeGradient(const double* x, double* res, double* gradient, const int& gradient_method) const
|
|
||||||
// If gradient_method == 1, use finite difference
|
|
||||||
// If gradient_method == 2, use analytic expresions
|
|
||||||
{
|
|
||||||
if (gradient_method == 1) {
|
|
||||||
double epsi = 1e-8;
|
|
||||||
double res_epsi[2];
|
|
||||||
double x_epsi[2];
|
|
||||||
computeResidual(x, res);
|
|
||||||
if (if_res_s) {
|
|
||||||
x_epsi[0] = x[0] + epsi;
|
|
||||||
x_epsi[1] = x[1];
|
|
||||||
computeResidual(x_epsi, res_epsi);
|
|
||||||
gradient[0] = (res_epsi[0] - res[0])/epsi;
|
|
||||||
x_epsi[0] = x[0];
|
|
||||||
x_epsi[1] = x[1] + epsi;
|
|
||||||
computeResidual(x_epsi, res_epsi);
|
|
||||||
gradient[1] = (res_epsi[0] - res[0])/epsi;
|
|
||||||
} else {
|
|
||||||
x_epsi[0] = x[0] + epsi;
|
|
||||||
x_epsi[1] = x[1];
|
|
||||||
computeResidual(x_epsi, res_epsi);
|
|
||||||
gradient[0] = (res_epsi[1] - res[1])/epsi;
|
|
||||||
x_epsi[0] = x[0];
|
|
||||||
x_epsi[1] = x[1] + epsi;
|
|
||||||
computeResidual(x_epsi, res_epsi);
|
|
||||||
gradient[1] = (res_epsi[1] - res[1])/epsi;
|
|
||||||
}
|
|
||||||
} else if (gradient_method == 2) {
|
|
||||||
double s = x[0];
|
|
||||||
double c = x[1];
|
|
||||||
double ff_ds_dc[2];
|
|
||||||
double ff = tm.fracFlowWithDer(s, c, cell, ff_ds_dc);
|
|
||||||
double mc_dc;
|
|
||||||
double mc = tm.computeMcWithDer(c, &mc_dc);
|
|
||||||
double dps = tm.polyprops_.deadPoreVol();
|
|
||||||
double rhor = tm.polyprops_.rockDensity();
|
|
||||||
double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
|
|
||||||
double ads;
|
|
||||||
double ads_dc;
|
|
||||||
if (c < cmax0) {
|
|
||||||
ads = tm.polyprops_.adsorbtion(cmax0);
|
|
||||||
ads_dc = 0;
|
|
||||||
} else {
|
|
||||||
ads = tm.polyprops_.adsorbtionWithDer(c, &ads_dc);
|
|
||||||
}
|
|
||||||
res[0] = s - s0 + dtpv*(outflux*ff + influx);
|
|
||||||
res[1] = (s - dps)*c - (s0 - dps)*c0
|
|
||||||
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
|
|
||||||
+ dtpv*(outflux*ff*mc + influx_polymer);
|
|
||||||
if (if_res_s) {
|
|
||||||
gradient[0] = 1 + dtpv*outflux*ff_ds_dc[0];
|
|
||||||
gradient[1] = dtpv*outflux*ff_ds_dc[1];
|
|
||||||
} else {
|
|
||||||
gradient[0] = c + dtpv*outflux*(ff_ds_dc[0])*mc;
|
|
||||||
gradient[1] = s - dps + rhor*((1.0 - porosity)/porosity)*ads_dc
|
|
||||||
+ dtpv*outflux*(ff_ds_dc[1]*mc + ff*mc_dc);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// setup 1d function, which is called by operator()
|
|
||||||
// For a given point x=(s,c) in the s,c plane, set up a piecewise linear curve wich starts
|
|
||||||
// from "x" with slope "direction", hits the bound of the rectangle
|
|
||||||
// [s_min,s_max]x[c_min,c_max] and continue in a straight line to "end_point". The curve is
|
|
||||||
// parametrized by t in [0, t_max], t_out is equal to t when the curve hits the bounding
|
|
||||||
// rectangle, x_out=(s_out, c_out) denotes the values of s and c at that point.
|
|
||||||
|
|
||||||
void setup(const double* x_arg, const double* direction_arg, const double* end_point_arg, const double* x_min_arg, const double* x_max_arg, double& t_max_out, double& t_out_out)
|
|
||||||
{
|
|
||||||
x[0] = x_arg[0];
|
|
||||||
x[1] = x_arg[1];
|
|
||||||
x_max[0] = x_max_arg[0];
|
|
||||||
x_max[1] = x_max_arg[1];
|
|
||||||
x_min[0] = x_min_arg[0];
|
|
||||||
x_min[1] = x_min_arg[1];
|
|
||||||
direction[0] = direction_arg[0];
|
|
||||||
direction[1] = direction_arg[1];
|
|
||||||
end_point[0] = end_point_arg[0];
|
|
||||||
end_point[1] = end_point_arg[1];
|
|
||||||
if ((end_point[0]-x[0])*direction[0] + (end_point[1]-x[1])*direction[1] < 0) {
|
|
||||||
direction[0] *= -1.0;
|
|
||||||
direction[1] *= -1.0;
|
|
||||||
}
|
|
||||||
if ((std::abs(direction[0]) + std::abs(direction[0])) == 0) {
|
|
||||||
direction[0] = end_point[0]-x[0];
|
|
||||||
direction[1] = end_point[1]-x[1];
|
|
||||||
}
|
|
||||||
bool t0_exists = true;
|
|
||||||
double t0;
|
|
||||||
if (direction[0] > 0) {
|
|
||||||
t0 = (x_max[0] - x[0])/direction[0];
|
|
||||||
} else if (direction[0] < 0) {
|
|
||||||
t0 = (x_min[0] - x[0])/direction[0];
|
|
||||||
} else {
|
|
||||||
t0_exists = false;
|
|
||||||
}
|
|
||||||
bool t1_exists = true;
|
|
||||||
double t1;
|
|
||||||
if (direction[1] > 0) {
|
|
||||||
t1 = (x_max[1] - x[1])/direction[1];
|
|
||||||
} else if (direction[1] < 0) {
|
|
||||||
t1 = (x_min[1] - x[1])/direction[1];
|
|
||||||
} else {
|
|
||||||
t1_exists = false;
|
|
||||||
}
|
|
||||||
if (t0_exists) {
|
|
||||||
if (t1_exists) {
|
|
||||||
t_out = std::min(t0, t1);
|
|
||||||
} else {
|
|
||||||
t_out = t0;
|
|
||||||
}
|
|
||||||
} else if (t1_exists) {
|
|
||||||
t_out = t1;
|
|
||||||
} else {
|
|
||||||
THROW("Direction illegal: is a zero vector.");
|
|
||||||
}
|
|
||||||
x_out[0] = x[0] + t_out*direction[0];
|
|
||||||
x_out[1] = x[1] + t_out*direction[1];
|
|
||||||
t_max = t_out + 1;
|
|
||||||
t_max_out = t_max;
|
|
||||||
t_out_out = t_out;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
// Compute x=(s,c) for a given t (t is the parameter for the piecewise linear curve)
|
|
||||||
void compute_x_of_t(double* x_of_t, const double t) const {
|
|
||||||
if (t <= t_out) {
|
|
||||||
x_of_t[0] = x[0] + t*direction[0];
|
|
||||||
x_of_t[1] = x[1] + t*direction[1];
|
|
||||||
} else {
|
|
||||||
x_of_t[0] = 1/(t_max-t_out)*((t_max - t)*x_out[0] + end_point[0]*(t - t_out));
|
|
||||||
x_of_t[1] = 1/(t_max-t_out)*((t_max - t)*x_out[1] + end_point[1]*(t - t_out));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
double operator()(const double t) const
|
|
||||||
{
|
|
||||||
double x_of_t[2];
|
|
||||||
compute_x_of_t(x_of_t, t);
|
|
||||||
double s;
|
|
||||||
double c;
|
|
||||||
s = x_of_t[0];
|
|
||||||
c = x_of_t[1];
|
|
||||||
if (if_res_s) {
|
|
||||||
return s - s0 + dtpv*(outflux*tm.fracFlow(s, c, cell) + influx);
|
|
||||||
} else {
|
|
||||||
double ff = tm.fracFlow(s, c, cell);
|
|
||||||
double mc = tm.computeMc(c);
|
|
||||||
double dps = tm.polyprops_.deadPoreVol();
|
|
||||||
double rhor = tm.polyprops_.rockDensity();
|
|
||||||
double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
|
|
||||||
double ads = tm.polyprops_.adsorbtion(std::max(c, cmax0));
|
|
||||||
return (s - dps)*c - (s0 - dps)*c0
|
|
||||||
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
|
|
||||||
+ dtpv*(outflux*ff*mc + influx_polymer);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
};
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
void TransportModelPolymer::solveSingleCell(const int cell)
|
void TransportModelPolymer::solveSingleCell(const int cell)
|
||||||
{
|
{
|
||||||
@ -585,8 +562,6 @@ namespace Opm
|
|||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
void TransportModelPolymer::solveSingleCellBracketing(int cell)
|
void TransportModelPolymer::solveSingleCellBracketing(int cell)
|
||||||
{
|
{
|
||||||
ResidualC res(*this, cell);
|
ResidualC res(*this, cell);
|
||||||
@ -621,18 +596,17 @@ namespace Opm
|
|||||||
// The tolerance falsi_tol is improved by (reduc_factor_falsi_tol * "previous residual") at each step
|
// The tolerance falsi_tol is improved by (reduc_factor_falsi_tol * "previous residual") at each step
|
||||||
double falsi_tol;
|
double falsi_tol;
|
||||||
const double reduc_factor_falsi_tol = 1e-4;
|
const double reduc_factor_falsi_tol = 1e-4;
|
||||||
const double gradient_method = 2; // method to compute derivative ( 1: finite difference, 2: formulae)
|
|
||||||
int iters_used_falsi = 0;
|
int iters_used_falsi = 0;
|
||||||
const int max_iters_split = 20;
|
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.
|
||||||
Residual residual(*this, cell);
|
ResidualEquation res_eq(*this, cell);
|
||||||
double x[2] = {saturation_[cell], concentration_[cell]};
|
double x[2] = {saturation_[cell], concentration_[cell]};
|
||||||
double res[2];
|
double res[2];
|
||||||
double mc;
|
double mc;
|
||||||
double ff;
|
double ff;
|
||||||
residual.computeResidual(x, res, mc, ff);
|
res_eq.computeResidual(x, 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;
|
||||||
@ -640,6 +614,7 @@ namespace Opm
|
|||||||
return;
|
return;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
double x_min[2] = { std::max(polyprops_.deadPoreVol(), smin_[2*cell]), 0.0 };
|
double x_min[2] = { std::max(polyprops_.deadPoreVol(), smin_[2*cell]), 0.0 };
|
||||||
double x_max[2] = { 1.0, polyprops_.cMax() };
|
double x_max[2] = { 1.0, polyprops_.cMax() };
|
||||||
double t;
|
double t;
|
||||||
@ -651,20 +626,26 @@ namespace Opm
|
|||||||
bool unsuccessfull_newton_step;
|
bool unsuccessfull_newton_step;
|
||||||
double x_new[2];
|
double x_new[2];
|
||||||
double res_new[2];
|
double res_new[2];
|
||||||
|
ResSOnCurve res_s_on_curve(res_eq);
|
||||||
|
ResCOnCurve res_c_on_curve(res_eq);
|
||||||
|
bool if_res_s;
|
||||||
|
|
||||||
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
|
// We first try a Newton step
|
||||||
if (residual.solveNewtonStep(x, x_new, gradient_method)) {
|
if (solveNewtonStep(x, res_eq, res, x_new)) {
|
||||||
residual.computeResidual(x_new, res_new, mc, ff);
|
res_eq.computeResidual(x_new, res_new, mc, ff);
|
||||||
unsuccessfull_newton_step = false;
|
unsuccessfull_newton_step = false;
|
||||||
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]) {
|
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;
|
unsuccessfull_newton_step = true;
|
||||||
} else {
|
} else {
|
||||||
x[0] = x_new[0];
|
x[0] = x_new[0];
|
||||||
x[1] = x_new[1];
|
x[1] = x_new[1];
|
||||||
res[0] = res_new[0];
|
res[0] = res_new[0];
|
||||||
res[1] = res_new[1];
|
res[1] = res_new[1];
|
||||||
iters_used_split += 1;
|
iters_used_split += 1;
|
||||||
|
if (norm(res_new) > 1e-1*norm(res)) {
|
||||||
|
unsuccessfull_newton_step = true;
|
||||||
|
}
|
||||||
}
|
}
|
||||||
} else {
|
} else {
|
||||||
unsuccessfull_newton_step = true;
|
unsuccessfull_newton_step = true;
|
||||||
@ -676,50 +657,55 @@ namespace Opm
|
|||||||
if (res[0] < -falsi_tol) {
|
if (res[0] < -falsi_tol) {
|
||||||
direction[0] = x_max[0] - x[0];
|
direction[0] = x_max[0] - x[0];
|
||||||
direction[1] = x_min[1] - x[1];
|
direction[1] = x_min[1] - x[1];
|
||||||
residual.if_res_s = true;
|
if_res_s = true;
|
||||||
} else if (res[0] > falsi_tol) {
|
} else if (res[0] > falsi_tol) {
|
||||||
direction[0] = x_min[0] - x[0];
|
direction[0] = x_min[0] - x[0];
|
||||||
direction[1] = x_max[1] - x[1];
|
direction[1] = x_max[1] - x[1];
|
||||||
residual.if_res_s = true;
|
if_res_s = true;
|
||||||
} else {
|
} else {
|
||||||
residual.computeGradient(x, res, gradient, gradient_method);
|
res_eq.computeGradientResS(x, res, gradient);
|
||||||
direction[0] = -gradient[1];
|
direction[0] = -gradient[1];
|
||||||
direction[1] = gradient[0];
|
direction[1] = gradient[0];
|
||||||
residual.if_res_s = false;
|
if_res_s = false;
|
||||||
}
|
}
|
||||||
} else {
|
} else {
|
||||||
falsi_tol = std::max(reduc_factor_falsi_tol*std::abs(res[1]), tol_);
|
falsi_tol = std::max(reduc_factor_falsi_tol*std::abs(res[1]), tol_);
|
||||||
if (res[1] < -falsi_tol) {
|
if (res[1] < -falsi_tol) {
|
||||||
direction[0] = x_max[0] - x[0];
|
direction[0] = x_max[0] - x[0];
|
||||||
direction[1] = x_max[1] - x[1];
|
direction[1] = x_max[1] - x[1];
|
||||||
residual.if_res_s = false;
|
if_res_s = false;
|
||||||
} else if (res[1] > falsi_tol) {
|
} else if (res[1] > falsi_tol) {
|
||||||
direction[0] = x_min[0] - x[0];
|
direction[0] = x_min[0] - x[0];
|
||||||
direction[1] = x_min[1] - x[1];
|
direction[1] = x_min[1] - x[1];
|
||||||
residual.if_res_s = false;
|
if_res_s = false;
|
||||||
} else {
|
} else {
|
||||||
residual.computeGradient(x, res, gradient, gradient_method);
|
res_eq.computeGradientResC(x, res, gradient);
|
||||||
direction[0] = -gradient[1];
|
direction[0] = -gradient[1];
|
||||||
direction[1] = gradient[0];
|
direction[1] = gradient[0];
|
||||||
residual.if_res_s = true;
|
if_res_s = true;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
if (residual.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[0];
|
||||||
end_point[1] = x_min[1];
|
end_point[1] = x_min[1];
|
||||||
residual.setup(x, direction, end_point, x_min, x_max, t_max, t_out);
|
res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, t_max, t_out);
|
||||||
if (residual(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[0];
|
||||||
end_point[1] = x_max[1];
|
end_point[1] = x_max[1];
|
||||||
residual.setup(x, direction, end_point, x_min, x_max, t_max, t_out);
|
res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, t_max, t_out);
|
||||||
if (residual(t_out) <= 0) {
|
if (res_s_on_curve(t_out) <= 0) {
|
||||||
t_max = t_out;
|
t_max = t_out;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
t = modifiedRegulaFalsi(res_s_on_curve, 0., t_max, maxit_, falsi_tol, iters_used_falsi);
|
||||||
|
if (std::abs(res_s_on_curve(t)) > falsi_tol) {
|
||||||
|
std::cout << "modifiedRegulaFalsi did not produce result under tolerance." << std::endl;
|
||||||
|
}
|
||||||
|
res_s_on_curve.curve.computeXOfT(x, t);
|
||||||
} else {
|
} else {
|
||||||
if (res[1] < 0) {
|
if (res[1] < 0) {
|
||||||
// We update the bounding box (Here we assume that the curve res_s(s,c)=0 is
|
// We update the bounding box (Here we assume that the curve res_s(s,c)=0 is
|
||||||
@ -731,8 +717,8 @@ namespace Opm
|
|||||||
//
|
//
|
||||||
end_point[0] = x_max[0];
|
end_point[0] = x_max[0];
|
||||||
end_point[1] = x_max[1];
|
end_point[1] = x_max[1];
|
||||||
residual.setup(x, direction, end_point, x_min, x_max, t_max, t_out);
|
res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, t_max, t_out);
|
||||||
if (residual(t_out) >= 0) {
|
if (res_c_on_curve(t_out) >= 0) {
|
||||||
t_max = t_out;
|
t_max = t_out;
|
||||||
}
|
}
|
||||||
} else {
|
} else {
|
||||||
@ -744,20 +730,20 @@ namespace Opm
|
|||||||
//
|
//
|
||||||
end_point[0] = x_min[0];
|
end_point[0] = x_min[0];
|
||||||
end_point[1] = x_min[1];
|
end_point[1] = x_min[1];
|
||||||
residual.setup(x, direction, end_point, x_min, x_max, t_max, t_out);
|
res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, t_max, t_out);
|
||||||
if (residual(t_out) <= 0) {
|
if (res_c_on_curve(t_out) <= 0) {
|
||||||
t_max = t_out;
|
t_max = t_out;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
t = modifiedRegulaFalsi(res_c_on_curve, 0., t_max, maxit_, falsi_tol, iters_used_falsi);
|
||||||
|
if (std::abs(res_c_on_curve(t)) > falsi_tol) {
|
||||||
|
std::cout << "modifiedRegulaFalsi did not produce result under tolerance." << std::endl;
|
||||||
|
}
|
||||||
|
res_c_on_curve.curve.computeXOfT(x, t);
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
t = modifiedRegulaFalsi(residual, 0., t_max, maxit_, falsi_tol, iters_used_falsi);
|
res_eq.computeResidual(x, res, mc, ff);
|
||||||
if (std::abs(residual(t)) > falsi_tol) {
|
|
||||||
std::cout << "modifiedRegulaFalsi did not produce result under tolerance." << std::endl;
|
|
||||||
}
|
|
||||||
residual.compute_x_of_t(x, t);
|
|
||||||
residual.computeResidual(x, res, mc, ff);
|
|
||||||
|
|
||||||
iters_used_split += 1;
|
iters_used_split += 1;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@ -936,8 +922,133 @@ namespace Opm
|
|||||||
} // namespace Opm
|
} // namespace Opm
|
||||||
|
|
||||||
|
|
||||||
|
namespace
|
||||||
|
{
|
||||||
|
|
||||||
|
// setup 1d function, which is called by operator()
|
||||||
|
// For a given point x=(s,c) in the s,c plane, set up a piecewise linear curve wich starts
|
||||||
|
// from "x" with slope "direction", hits the bound of the rectangle
|
||||||
|
// [s_min,s_max]x[c_min,c_max] and continue in a straight line to "end_point". The curve is
|
||||||
|
// parametrized by t in [0, t_max], t_out is equal to t when the curve hits the bounding
|
||||||
|
// rectangle, x_out=(s_out, c_out) denotes the values of s and c at that point.
|
||||||
|
CurveInSCPlane::CurveInSCPlane()
|
||||||
|
{
|
||||||
|
}
|
||||||
|
|
||||||
|
void CurveInSCPlane::setup(const double* x, const double* direction,
|
||||||
|
const double* end_point, const double* x_min,
|
||||||
|
const double* x_max, double& t_max_out,
|
||||||
|
double& t_out_out)
|
||||||
|
{
|
||||||
|
x_[0] = x[0];
|
||||||
|
x_[1] = x[1];
|
||||||
|
x_max_[0] = x_max[0];
|
||||||
|
x_max_[1] = x_max[1];
|
||||||
|
x_min_[0] = x_min[0];
|
||||||
|
x_min_[1] = x_min[1];
|
||||||
|
direction_[0] = direction[0];
|
||||||
|
direction_[1] = direction[1];
|
||||||
|
end_point_[0] = end_point[0];
|
||||||
|
end_point_[1] = end_point[1];
|
||||||
|
if ((end_point_[0]-x_[0])*direction_[0] + (end_point_[1]-x_[1])*direction_[1] < 0) {
|
||||||
|
direction_[0] *= -1.0;
|
||||||
|
direction_[1] *= -1.0;
|
||||||
|
}
|
||||||
|
if ((std::abs(direction_[0]) + std::abs(direction_[0])) == 0) {
|
||||||
|
direction_[0] = end_point_[0]-x_[0];
|
||||||
|
direction_[1] = end_point_[1]-x_[1];
|
||||||
|
}
|
||||||
|
bool t0_exists = true;
|
||||||
|
double t0;
|
||||||
|
if (direction_[0] > 0) {
|
||||||
|
t0 = (x_max_[0] - x_[0])/direction_[0];
|
||||||
|
} else if (direction_[0] < 0) {
|
||||||
|
t0 = (x_min_[0] - x_[0])/direction_[0];
|
||||||
|
} else {
|
||||||
|
t0_exists = false;
|
||||||
|
}
|
||||||
|
bool t1_exists = true;
|
||||||
|
double t1;
|
||||||
|
if (direction_[1] > 0) {
|
||||||
|
t1 = (x_max_[1] - x_[1])/direction_[1];
|
||||||
|
} else if (direction[1] < 0) {
|
||||||
|
t1 = (x_min_[1] - x_[1])/direction_[1];
|
||||||
|
} else {
|
||||||
|
t1_exists = false;
|
||||||
|
}
|
||||||
|
if (t0_exists) {
|
||||||
|
if (t1_exists) {
|
||||||
|
t_out_ = std::min(t0, t1);
|
||||||
|
} else {
|
||||||
|
t_out_ = t0;
|
||||||
|
}
|
||||||
|
} else if (t1_exists) {
|
||||||
|
t_out_ = t1;
|
||||||
|
} else {
|
||||||
|
THROW("Direction illegal: is a zero vector.");
|
||||||
|
}
|
||||||
|
x_out_[0] = x_[0] + t_out_*direction_[0];
|
||||||
|
x_out_[1] = x_[1] + t_out_*direction_[1];
|
||||||
|
t_max_ = t_out_ + 1;
|
||||||
|
t_max_out = t_max_;
|
||||||
|
t_out_out = t_out_;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
// Compute x=(s,c) for a given t (t is the parameter for the piecewise linear curve)
|
||||||
|
void CurveInSCPlane::computeXOfT(double* x_of_t, const double t) const {
|
||||||
|
if (t <= t_out_) {
|
||||||
|
x_of_t[0] = x_[0] + t*direction_[0];
|
||||||
|
x_of_t[1] = x_[1] + t*direction_[1];
|
||||||
|
} else {
|
||||||
|
x_of_t[0] = 1/(t_max_-t_out_)*((t_max_ - t)*x_out_[0] + end_point_[0]*(t - t_out_));
|
||||||
|
x_of_t[1] = 1/(t_max_-t_out_)*((t_max_ - t)*x_out_[1] + end_point_[1]*(t - t_out_));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
ResSOnCurve::ResSOnCurve(const Opm::TransportModelPolymer::ResidualEquation& res_eq)
|
||||||
|
: res_eq_(res_eq)
|
||||||
|
{
|
||||||
|
}
|
||||||
|
|
||||||
|
double ResSOnCurve::operator()(const double t) const
|
||||||
|
{
|
||||||
|
double x_of_t[2];
|
||||||
|
curve.computeXOfT(x_of_t, t);
|
||||||
|
return res_eq_.computeResidualS(x_of_t);
|
||||||
|
}
|
||||||
|
|
||||||
|
ResCOnCurve::ResCOnCurve(const Opm::TransportModelPolymer::ResidualEquation& res_eq)
|
||||||
|
: res_eq_(res_eq)
|
||||||
|
{
|
||||||
|
}
|
||||||
|
|
||||||
|
double ResCOnCurve::operator()(const double t) const
|
||||||
|
{
|
||||||
|
double x_of_t[2];
|
||||||
|
curve.computeXOfT(x_of_t, t);
|
||||||
|
return res_eq_.computeResidualC(x_of_t);
|
||||||
|
}
|
||||||
|
|
||||||
|
bool solveNewtonStep(const double* x, const Opm::TransportModelPolymer::ResidualEquation& res_eq, const double* res, double* x_new) {
|
||||||
|
|
||||||
|
double res_s_ds_dc[2];
|
||||||
|
double res_c_ds_dc[2];
|
||||||
|
|
||||||
|
res_eq.computeJacobiRes(x, res_s_ds_dc, res_c_ds_dc);
|
||||||
|
|
||||||
|
double det = res_s_ds_dc[0]*res_c_ds_dc[1] - res_c_ds_dc[0]*res_s_ds_dc[1];
|
||||||
|
if (std::abs(det) < 1e-8) {
|
||||||
|
return false;
|
||||||
|
} else {
|
||||||
|
x_new[0] = x[0] - (res[0]*res_c_ds_dc[1] - res[1]*res_s_ds_dc[1])/det;
|
||||||
|
x_new[1] = x[1] - (res[1]*res_s_ds_dc[0] - res[0]*res_c_ds_dc[0])/det;
|
||||||
|
return true;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
} // Anonymous namespace
|
||||||
|
|
||||||
|
|
||||||
/* Local Variables: */
|
/* Local Variables: */
|
||||||
|
|||||||
@ -66,6 +66,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);
|
||||||
|
class ResidualEquation;
|
||||||
|
|
||||||
|
|
||||||
private:
|
private:
|
||||||
@ -94,8 +95,6 @@ namespace Opm
|
|||||||
struct ResidualC;
|
struct ResidualC;
|
||||||
struct ResidualS;
|
struct ResidualS;
|
||||||
|
|
||||||
// Residual function which is used in splitting method
|
|
||||||
struct Residual;
|
|
||||||
|
|
||||||
double fracFlow(double s, double c, int cell) const;
|
double fracFlow(double s, double c, int cell) const;
|
||||||
double fracFlowWithDer(double s, double c, int cell, double* der) const;
|
double fracFlowWithDer(double s, double c, int cell, double* der) const;
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user