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_ ; }
102 double& rockdensity() {
return rockdensity_ ; }
103 double rockdensity()
const {
return rockdensity_ ; }
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]; }
120 double* porevol() {
return porevol_ ; }
121 double porevol(
int cell)
const {
return porevol_[cell] ; }
123 double* porosity() {
return porosity_ ; }
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());
220 store_.rockdensity() = fluid.rockdensity();
223 void makefhfQPeriodic(
const std::vector<int>& p_faces,
const std::vector<int>& hf_faces,
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];
266 enum { DofPerCell = 1 };
269 initResidual(
const int c,
double* Fs,
double* Fc)
const {
275 template <
class ReservoirState,
278 fluxConnection(
const ReservoirState& state ,
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>
372 accumulation(
const Grid& g,
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 ,
395 sourceTerms(
const Grid& g ,
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>
429 initGravityTrans(
const Grid& g ,
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) {
450 store_.trans(f) = 1.0 / store_.trans(f);
454 this->computeStaticGravity(g);
462 template <
class ReservoirState,
464 class JacobianSystem>
466 initStep(
const ReservoirState& state,
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>
495 initIteration(
const ReservoirState& state,
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]);
522 store_.cmax(cell) = cmax;
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];
555 store_.dmobwatdc(cell) = dmobwatdc;
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 >
571 finishIteration(
const ReservoirState& state,
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>
583 finishStep(
const Grid& g ,
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_ ;
709 bool init_step_use_previous_sol_;
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
Definition: SinglePointUpwindTwoPhasePolymer.hpp:48
Definition: SinglePointUpwindTwoPhasePolymer.hpp:188