36 #ifndef OPM_SINGLEPOINTUPWINDTWOPHASEPOLYMER_HPP_HEADER
37 #define OPM_SINGLEPOINTUPWINDTWOPHASEPOLYMER_HPP_HEADER
47 namespace polymer_reorder {
51 : drho_(0.0), rockdensity_(0.0), mob_(0),
52 dmobds_(0), dmobwatdc_(0), mc_(0),
53 dmcdc_(0), porevol_(0), porosity_(0), dg_(0), sw_(0), c_(0), cmax_(0),
54 ds_(0), dsc_(0), dcads_(0), dcadsdc_(0), pc_(0), dpc_(0),
66 alloc_sz += 1 * totconn;
76 alloc_sz += 1 * totconn;
77 data_.resize(alloc_sz);
80 dmobds_ = mob_ + (2 * nc );
81 dmobwatdc_ = dmobds_ + (2 * nc );
82 mc_ = dmobwatdc_ + (1 * nc );
83 dmcdc_ = mc_ + (1 * nc );
84 porevol_ = dmcdc_ + (1 * nc );
85 porosity_ = porevol_ + (1 * nc );
86 dg_ = porosity_ + (1 * nc );
87 sw_ = dg_ + (1 * totconn);
89 cmax_ = c_ + (1 * nc );
90 ds_ = cmax_ + (1 * nc );
91 dsc_ = ds_ + (1 * nc );
92 dcads_ = dsc_ + (1 * nc );
93 dcadsdc_ = dcads_ + (1 * nc );
94 pc_ = dcadsdc_ + (1 * nc );
95 dpc_ = pc_ + (1 * nc );
96 trans_ = dpc_ + (1 * nc );
99 double&
drho () {
return drho_ ; }
100 double drho ()
const {
return drho_ ; }
105 double*
mob (
int cell) {
return mob_ + (2*cell + 0); }
106 const double*
mob (
int cell)
const {
return mob_ + (2*cell + 0); }
108 double*
dmobds (
int cell) {
return dmobds_ + (2*cell + 0); }
109 const double*
dmobds (
int cell)
const {
return dmobds_ + (2*cell + 0); }
111 double&
dmobwatdc (
int cell) {
return dmobwatdc_[cell]; }
112 double dmobwatdc (
int cell)
const {
return dmobwatdc_[cell]; }
114 double&
mc (
int cell) {
return mc_[cell]; }
115 double mc (
int cell)
const {
return mc_[cell]; }
117 double&
dmcdc (
int cell) {
return dmcdc_[cell]; }
118 double dmcdc (
int cell)
const {
return dmcdc_[cell]; }
121 double porevol(
int cell)
const {
return porevol_[cell] ; }
124 double porosity(
int cell)
const {
return porosity_[cell] ; }
127 double&
dg(
int i) {
return dg_[i] ; }
128 double dg(
int i)
const {
return dg_[i] ; }
130 double&
sw(
int cell) {
return sw_[cell] ; }
131 double sw(
int cell)
const {
return sw_[cell] ; }
133 double&
c(
int cell) {
return c_[cell] ; }
134 double c(
int cell)
const {
return c_[cell] ; }
136 double&
cmax(
int cell) {
return cmax_[cell] ; }
137 double cmax(
int cell)
const {
return cmax_[cell] ; }
139 double&
ds(
int cell) {
return ds_[cell] ; }
140 double ds(
int cell)
const {
return ds_[cell] ; }
142 double&
dsc(
int cell) {
return dsc_[cell] ; }
143 double dsc(
int cell)
const {
return dsc_[cell] ; }
145 double&
dcads(
int cell) {
return dcads_[cell] ; }
146 double dcads(
int cell)
const {
return dcads_[cell] ; }
148 double&
dcadsdc(
int cell) {
return dcadsdc_[cell] ; }
149 double dcadsdc(
int cell)
const {
return dcadsdc_[cell] ; }
151 double&
pc(
int cell) {
return pc_[cell] ; }
152 double pc(
int cell)
const {
return pc_[cell] ; }
154 double&
dpc(
int cell) {
return dpc_[cell] ; }
155 double dpc(
int cell)
const {
return dpc_[cell] ; }
157 double&
trans(
int f) {
return trans_[f] ; }
158 double trans(
int f)
const {
return trans_[f] ; }
162 double rockdensity_ ;
182 std::vector<double> data_;
187 template <
class TwophaseFlu
idPolymer>
190 template <
class Gr
id>
193 const std::vector<double>& porevol ,
194 const double* grav = 0,
195 const bool guess_previous =
true)
198 f2hf_ (2 * g.number_of_faces, -1) ,
199 store_ (g.number_of_cells,
200 g.cell_facepos[ g.number_of_cells ]),
201 init_step_use_previous_sol_(guess_previous) ,
206 store_.
drho() = fluid_.density(0) - fluid_.density(1);
209 for (
int c = 0, i = 0; c < g.number_of_cells; ++c) {
210 for (; i < g.cell_facepos[c + 1]; ++i) {
211 const int f = g.cell_faces[i];
212 const int p = 1 - (g.face_cells[2*f + 0] == c);
217 std::copy(porevol.begin(), porevol.end(), store_.
porevol());
218 const double* poro = fluid.porosity();
219 std::copy(poro, poro + g.number_of_cells, store_.
porosity());
224 const std::vector<int>& nb_faces)
226 std::vector<int> nbhf(hf_faces.size());
227 for(
unsigned int i=0; i<p_faces.size(); ++i){
228 int nbf = nb_faces[i];
229 if(f2hf_[2*nbf] == -1){
230 nbhf[i] = f2hf_[2*nbf+1];
232 assert(f2hf_[2*nbf+1]==-1);
233 nbhf[i] = f2hf_[2*nbf];
236 for(
unsigned int i=0; i<p_faces.size(); ++i){
239 int hf = hf_faces[i];
242 if(f2hf_[2*f] == hf){
243 assert(f2hf_[2*f+1]==-1);
245 assert(f2hf_[2*f]==-1);
250 if(f2hf_[2*f+1]== hf){
251 assert(f2hf_[2*f]==-1);
253 assert(f2hf_[2*f+1]==-1);
254 f2hf_[2*f+1]=nbhf[i];
275 template <
class ReservoirState,
291 const int *n = g.face_cells + (2 * f);
292 double dflux = state.faceflux()[f];
293 double gflux = gravityFlux(f);
294 double pcflux, dpcflux[2];
295 capFlux(f, n, pcflux, dpcflux);
299 double m[2], dmds[2], dmobwatdc;
301 upwindMobility(dflux, gflux, n, pix, m, dmds, dmobwatdc, mc, dmcdc);
303 assert ((m[0] >= 0.0) && (m[1] >= 0.0));
305 double mt = m[0] + m[1];
308 double sgn = 2.0*(n[0] == cell) - 1.0;
313 double f1 = m[0] / mt;
314 const double v1 = dflux + m[1]*gflux;
317 F[0] += dt * f1 * v1;
318 F[1] += dt * mc * f1 * v1;
326 dFsds[0] = &dFd1[0]; dFsds[1] = &dFd2[0];
327 dFsdc[0] = &dFd1[1]; dFsdc[1] = &dFd2[1];
328 dFcds[0] = &dFd1[2]; dFcds[1] = &dFd2[2];
329 dFcdc[0] = &dFd1[3]; dFcdc[1] = &dFd2[3];
331 dFd1[0] += sgn*dt * f1 * dpcflux[0] * m[1];
332 dFd2[0] += sgn*dt * f1 * dpcflux[1] * m[1];
333 dFd1[2] += sgn*dt * f1 * mc * dpcflux[0] * m[1];
334 dFd2[2] += sgn*dt * f1 * mc * dpcflux[1] * m[1];
338 dFsds[0] = &dFd2[0]; dFsds[1] = &dFd1[0];
339 dFsdc[0] = &dFd2[1]; dFsdc[1] = &dFd1[1];
340 dFcds[0] = &dFd2[2]; dFcds[1] = &dFd1[2];
341 dFcdc[0] = &dFd2[3]; dFcdc[1] = &dFd1[3];
343 dFd1[0] += sgn*dt * f1 * dpcflux[1] * m[1];
344 dFd2[0] += sgn*dt * f1 * dpcflux[0] * m[1];
345 dFd1[2] += sgn*dt * f1 * mc * dpcflux[1] * m[1];
346 dFd2[2] += sgn*dt * f1 * mc * dpcflux[0] * m[1];
352 *dFsds[ pix[0] ] += dt * (1 - f1) / mt * v1 * dmds[0];
354 *dFcds[ pix[0] ] += dt * (1 - f1) / mt * v1 * mc * dmds[0];
357 *dFsds[ pix[1] ] -= dt * f1 / mt * v1 * dmds[1];
358 *dFsds[ pix[1] ] += dt * f1 * gflux * dmds[1];
360 *dFcds[ pix[1] ] -= dt * f1 / mt * v1 * mc * dmds[1];
361 *dFcds[ pix[1] ] += dt * f1 * gflux * mc * dmds[1];
364 *dFsdc[ pix[0] ] += dt * (1 - f1) / mt * v1 * dmobwatdc;
366 *dFcdc[ pix[0] ] += dt * (1 - f1) / mt * v1 * mc * dmobwatdc;
367 *dFcdc[ pix[0] ] += dt * f1 * v1 * dmcdc;
370 template <
class Gr
id>
379 const double pv = store_.
porevol(cell);
380 const double dps = fluid_.deadporespace();
381 const double rhor = fluid_.rockdensity();
382 const double poro = store_.
porosity(cell);
384 F[0] += pv * store_.
ds(cell);
385 F[1] += pv * (1 - dps) * store_.
dsc(cell) + rhor*(1 - poro)/poro*pv*store_.
dcads(cell);
388 dF[2] += pv * (1 - dps) * store_.
c(cell);
389 dF[3] += pv * (1 - dps) * store_.
sw(cell) + rhor*(1 - poro)/poro*pv*store_.
dcadsdc(cell);
392 template <
class Grid ,
396 const SourceTerms* src,
404 double dflux = -src->flux[i];
408 *F += dt * dflux * src->saturation[2*i + 0];
411 const int cell = src->cell[i];
412 const double* m = store_.
mob (cell);
413 const double* dm = store_.
dmobds(cell);
415 const double mt = m[0] + m[1];
417 assert (! ((m[0] < 0) || (m[1] < 0)));
420 const double f = m[0] / mt;
421 const double df = ((1 - f)*dm[0] - f*dm[1]) / mt;
423 *F += dt * dflux * f;
424 *J += dt * dflux * df;
427 template <
class Gr
id>
430 const std::vector<double> & htrans) {
432 assert (htrans.size() ==
433 static_cast<std::vector<double>::size_type
>(g.cell_facepos[ g.number_of_cells ]));
435 for (
int f = 0; f < g.number_of_faces; ++f) {
436 store_.
trans(f) = 0.0;
439 for (
int c = 0, i = 0; c < g.number_of_cells; ++c) {
440 for (; i < g.cell_facepos[c + 1]; ++i) {
441 int f = g.cell_faces[i];
443 assert (htrans[i] > 0.0);
445 store_.
trans(f) += 1.0 / htrans[i];
449 for (
int f = 0; f < g.number_of_faces; ++f) {
454 this->computeStaticGravity(g);
462 template <
class ReservoirState,
464 class JacobianSystem>
468 JacobianSystem& sys ) {
472 typename JacobianSystem::vector_type& x =
473 sys.vector().writableSolution();
475 assert (x.size() == (::std::size_t) (2*g.number_of_cells));
477 if (init_step_use_previous_sol_) {
478 std::fill(x.begin(), x.end(), 0.0);
480 std::fill(x.begin(), x.end(), 0.0);
481 const std::vector<double>& s = state.saturation();
482 const std::vector<double>& c = state.concentration();
483 for (
int cell = 0, ncell = g.number_of_cells; cell < ncell; ++cell) {
485 x[2*cell + 0] = 0.5 - s[2*cell + 0];
486 x[2*cell + 1] = 1e-5 - c[cell];
491 template <
class ReservoirState,
493 class JacobianSystem>
497 JacobianSystem& sys) {
507 const typename JacobianSystem::vector_type& x =
508 sys.vector().solution();
509 const ::std::vector<double>&
sat = state.saturation();
510 const ::std::vector<double>&
cpoly = state.concentration();
511 const ::std::vector<double>& cmaxpoly = state.maxconcentration();
513 bool in_range =
true;
514 for (
int cell = 0; cell < g.number_of_cells; ++cell) {
516 store_.
ds(cell) = x[2*cell + 0];
517 s[0] = sat[cell*2 + 0] + x[2*cell + 0];
518 c = cpoly[cell] + x[2*cell + 1];
519 store_.
sw(cell) = s[0];
521 cmax = std::max(c, cmaxpoly[cell]);
523 store_.
dsc(cell) = s[0]*c - sat[cell*2 + 0]*cpoly[cell];
526 fluid_.adsorption(cpoly[cell], cmaxpoly[cell], cads, dcadsdc);
527 store_.
dcads(cell) = -cads;
528 fluid_.adsorption(c, cmax, cads, dcadsdc);
529 store_.
dcads(cell) += cads;
530 store_.
dcadsdc(cell) = dcadsdc;
531 double s_min = fluid_.s_min(cell);
532 double s_max = fluid_.s_max(cell);
534 if ( s[0] < (s_min - sat_tol_) || s[0] > (s_max + sat_tol_) ) {
543 s[0] = std::max(s_min, s[0]);
544 s[0] = std::min(s_max, s[0]);
547 fluid_.mobility(cell, s, c, cmax, mob, dmobds, dmobwatdc);
548 fluid_.computeMc(c, mc, dmcdc);
549 fluid_.pc(cell, s, pc, dpc);
551 store_.
mob (cell)[0] = mob [0];
552 store_.
mob (cell)[1] = mob [1];
553 store_.
dmobds(cell)[0] = dmobds[0*2 + 0];
554 store_.
dmobds(cell)[1] = -dmobds[1*2 + 1];
556 store_.
mc(cell) = mc;
557 store_.
dmcdc(cell) = dmcdc;
558 store_.
pc(cell) = pc;
559 store_.
dpc(cell) = dpc;
562 std::cout <<
"Warning: initIteration() - s was clamped in some cells.\n";
567 template <
class ReservoirState,
569 class NewtonIterate >
573 NewtonIterate& it ) {
575 (void) state; (void) g; (void) it;
576 typedef typename NewtonIterate::vector_type vector_t;
579 template <
class Grid ,
580 class SolutionVector,
581 class ReservoirState>
584 const SolutionVector& x ,
585 ReservoirState& state) {
587 double *s = &state.saturation()[0*2 + 0];
588 double *c = &state.concentration()[0*1 + 0];
589 double *
cmax = &state.maxconcentration()[0*1 + 0];
591 for (
int cell = 0; cell < g.number_of_cells; ++cell, s += 2, c += 1, cmax +=1) {
592 s[0] += x[2*cell + 0];
593 c[0] += x[2*cell + 1];
594 cmax[0] = std::max(c[0], cmax[0]);
595 double s_min = fluid_.s_min(cell);
596 double s_max = fluid_.s_max(cell);
597 assert(s[0] >= s_min - sat_tol_);
598 assert(s[0] <= s_max + sat_tol_);
599 s[0] = std::max(s_min, s[0]);
600 s[0] = std::min(s_max, s[0]);
607 upwindMobility(
const double dflux,
615 double& dmcdc)
const {
616 bool equal_sign = ( (! (dflux < 0)) && (! (gflux < 0)) ) ||
617 ( (! (dflux > 0)) && (! (gflux > 0)) );
621 if (! (dflux < 0) && ! (gflux < 0)) { pix[0] = 0; }
624 m[0] = store_.
mob(n[ pix[0] ]) [ 0 ];
625 mc = store_.
mc(n[ pix[0] ]);
627 if (! (dflux - m[0]*gflux < 0)) { pix[1] = 0; }
630 m[1] = store_.
mob(n[ pix[1] ]) [ 1 ];
634 if (! (dflux < 0) && ! (gflux > 0)) { pix[1] = 0; }
637 m[1] = store_.
mob(n[ pix[1] ]) [ 1 ];
639 if (dflux + m[1]*gflux > 0) { pix[0] = 0; }
642 m[0] = store_.
mob(n[ pix[0] ]) [ 0 ];
643 mc = store_.
mc(n[ pix[0] ]);
646 dmds[0] = store_.
dmobds(n[ pix[0] ]) [ 0 ];
647 dmds[1] = store_.
dmobds(n[ pix[1] ]) [ 1 ];
648 dmobwatdc = store_.
dmobwatdc(n[ pix[0] ]);
649 dmcdc = store_.
dmcdc(n[ pix[0] ]);
652 template <
class Gr
id>
654 computeStaticGravity(
const Grid& g) {
656 const int d = g.dimensions;
658 for (
int c = 0, i = 0; c < g.number_of_cells; ++c) {
659 const double* cc = g.cell_centroids + (c * d);
661 for (; i < g.cell_facepos[c + 1]; ++i) {
662 const int f = g.cell_faces[i];
663 const double* fc = g.face_centroids + (f * d);
666 for (
int j = 0; j < d; ++j) {
667 dg += gravity_[j] * (fc[j] - cc[j]);
670 store_.
dg(i) = store_.
trans(f) * dg;
676 gravityFlux(
const int f)
const {
680 int i1 = f2hf_[2*f + 0];
681 int i2 = f2hf_[2*f + 1];
683 assert ((i1 >= 0) && (i2 >= 0));
685 gflux = store_.
dg(i1) - store_.
dg(i2);
686 gflux *= store_.
drho();
694 capFlux(
const int f,
const int* n,
double& pcflux,
double* dpcflux)
const {
698 assert ((i1 >= 0) && (i2 >= 0));
700 pcflux = store_.
trans(f)*(store_.
pc(i2) - store_.
pc(i1));
701 dpcflux[0] = -store_.
trans(f)*store_.
dpc(i1);
702 dpcflux[1] = store_.
trans(f)*store_.
dpc(i2);
705 TwophaseFluidPolymer fluid_ ;
706 const double* gravity_;
707 std::vector<int> f2hf_ ;
708 polymer_reorder::ModelParameterStorage store_ ;
709 bool init_step_use_previous_sol_;
void finishIteration(const ReservoirState &state, const Grid &g, NewtonIterate &it)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:571
double porosity(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:124
double porevol(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:121
double drho() const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:100
double dsc(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:143
double & ds(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:139
Definition: CompressibleTpfaPolymer.hpp:32
std::vector< double > & sat
Definition: GravityColumnSolverPolymer_impl.hpp:76
double & drho()
Definition: SinglePointUpwindTwoPhasePolymer.hpp:99
double * porevol()
Definition: SinglePointUpwindTwoPhasePolymer.hpp:120
void makefhfQPeriodic(const std::vector< int > &p_faces, const std::vector< int > &hf_faces, const std::vector< int > &nb_faces)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:223
Definition: SinglePointUpwindTwoPhasePolymer.hpp:266
std::vector< double > & cmax
Definition: GravityColumnSolverPolymer_impl.hpp:78
double c(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:134
double * dmobds(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:108
double trans(int f) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:158
void finishStep(const Grid &g, const SolutionVector &x, ReservoirState &state)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:583
double dmobwatdc(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:112
double ds(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:140
double & mc(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:114
double & pc(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:151
double dg(int i) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:128
const double * dmobds(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:109
double * porosity()
Definition: SinglePointUpwindTwoPhasePolymer.hpp:123
double & dcads(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:145
void initStep(const ReservoirState &state, const Grid &g, JacobianSystem &sys)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:466
double dcadsdc(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:149
std::vector< double > & cpoly
Definition: GravityColumnSolverPolymer_impl.hpp:77
void accumulation(const Grid &g, const int cell, double *F, double *dF) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:372
void sourceTerms(const Grid &g, const SourceTerms *src, const int i, const double dt, double *J, double *F) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:395
double dcads(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:146
bool initIteration(const ReservoirState &state, const Grid &g, JacobianSystem &sys)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:495
double & cmax(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:136
double dmcdc(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:118
Definition: SinglePointUpwindTwoPhasePolymer.hpp:188
Definition: SinglePointUpwindTwoPhasePolymer.hpp:48
double & c(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:133
double dpc(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:155
double pc(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:152
void initGravityTrans(const Grid &g, const std::vector< double > &htrans)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:429
double rockdensity() const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:103
const double * mob(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:106
double & dmcdc(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:117
double & dsc(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:142
ModelParameterStorage(int nc, int totconn)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:50
double * mob(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:105
void initResidual(const int c, double *Fs, double *Fc) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:269
void fluxConnection(const ReservoirState &state, const Grid &g, const double dt, const int cell, const int f, double *F, double *dFd1, double *dFd2) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:278
double & dg(int i)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:127
double & trans(int f)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:157
double & dcadsdc(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:148
double & sw(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:130
double mc(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:115
double & dpc(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:154
SinglePointUpwindTwoPhasePolymer(const TwophaseFluidPolymer &fluid, const Grid &g, const std::vector< double > &porevol, const double *grav=0, const bool guess_previous=true)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:191
double sw(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:131
double cmax(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:137
double & dmobwatdc(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:111
double & rockdensity()
Definition: SinglePointUpwindTwoPhasePolymer.hpp:102