42#ifndef OPM_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER
43#define OPM_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER
92 : drho_(0.0), mob_(0), dmob_(0),
93 porevol_(0), dg_(0), ds_(0), pc_(0), dpc_(0), trans_(0),
101 alloc_sz += 1 * totconn;
105 alloc_sz += 1 * totconn;
106 data_.resize(alloc_sz);
109 dmob_ = mob_ + (2 * nc );
110 porevol_ = dmob_ + (2 * nc );
111 dg_ = porevol_ + (1 * nc );
112 ds_ = dg_ + (1 * totconn);
113 pc_ = ds_ + (1 * nc );
114 dpc_ = pc_ + (1 * nc );
115 trans_ = dpc_ + (1 * nc );
123 double&
drho () {
return drho_ ; }
128 double drho ()
const {
return drho_ ; }
136 double*
mob (
int c) {
return mob_ + (2*c + 0); }
143 const double*
mob (
int c)
const {
return mob_ + (2*c + 0); }
152 double*
dmob (
int c) {
return dmob_ + (2*c + 0); }
160 const double*
dmob (
int c)
const {
return dmob_ + (2*c + 0); }
172 double porevol(
int c)
const {
return porevol_[c] ; }
182 double&
dg(
int i) {
return dg_[i] ; }
190 double dg(
int i)
const {
return dg_[i] ; }
199 double&
ds(
int c) {
return ds_[c] ; }
207 double ds(
int c)
const {
return ds_[c] ; }
215 double&
pc(
int c) {
return pc_[c] ; }
222 double pc(
int c)
const {
return pc_[c] ; }
231 double&
dpc(
int c) {
return dpc_[c] ; }
239 double dpc(
int c)
const {
return dpc_[c] ; }
248 double&
trans(
int f) {
return trans_[f] ; }
256 double trans(
int f)
const {
return trans_[f] ; }
272 std::vector<double> data_;
277 template <
class TwophaseFlu
id>
280 template <
class Gr
id>
283 const std::vector<double>& porevol ,
284 const double* grav = 0,
285 const bool guess_previous =
true)
288 f2hf_ (2 * g.number_of_faces, -1) ,
289 store_ (g.number_of_cells,
290 g.cell_facepos[ g.number_of_cells ]),
291 init_step_use_previous_sol_(guess_previous),
296 store_.
drho() = fluid_.density(0) - fluid_.density(1);
299 for (
int c = 0, i = 0; c < g.number_of_cells; ++c) {
300 for (; i < g.cell_facepos[c + 1]; ++i) {
301 const int f = g.cell_faces[i];
302 const int p = 1 - (g.face_cells[2*f + 0] == c);
307 std::copy(porevol.begin(), porevol.end(), store_.
porevol());
312 const std::vector<int>& hf_faces,
313 const std::vector<int>& nb_faces)
315 if (p_faces.empty()) {
318 assert (p_faces.size() == hf_faces.size());
319 assert (hf_faces.size() == nb_faces.size());
321 std::vector<int> nbhf(hf_faces.size());
323 for (std::vector<int>::size_type i = 0; i < p_faces.size(); ++i) {
324 const int nbf = nb_faces[i];
326 assert (2*std::vector<int>::size_type(nbf) + 1 < f2hf_.size());
327 assert ((f2hf_[2*nbf + 0] < 0) ^ (f2hf_[2*nbf + 1] < 0));
329 const int p = (f2hf_[2*nbf + 0] < 0) ? 1 : 0;
330 nbhf[ i ] = f2hf_[2*nbf + p];
333 for (std::vector<int>::size_type i = 0; i < p_faces.size(); ++i) {
334 const int f = p_faces [i];
335 const int hf = hf_faces[i];
339 assert (2*std::vector<int>::size_type(f) + 1 < f2hf_.size());
341 assert ((f2hf_[2*f + 0] < 0 ) ^ (f2hf_[2*f + 1] < 0 ));
342 assert ((f2hf_[2*f + 0] == hf) ^ (f2hf_[2*f + 1] == hf));
344 const int p = (f2hf_[2*f + 0] == hf) ? 1 : 0;
346 f2hf_[2*f + p] = nbhf[ i ];
374 const int *n = g.face_cells + (2 * f);
376 double gflux = gravityFlux(f);
377 double pcflux,dpcflux[2];
378 capFlux(f,n, pcflux, dpcflux);
383 upwindMobility(dflux, gflux, n, pix, m, dm);
385 assert (! ((m[0] < 0) || (m[1] < 0)));
387 double mt = m[0] + m[1];
390 double sgn = 2.0*(n[0] == c) - 1.0;
395 double f1 = m[0] / mt;
396 const double v1 = dflux + m[1]*gflux;
404 J[0] = J1; J[1] = J2;
406 J1[0*2 + 0] += sgn*dt * f1 * dpcflux[0] * m[1];
407 J2[0*2 + 0] += sgn*dt * f1 * dpcflux[1] * m[1];
409 J[0] = J2; J[1] = J1;
411 J1[0*2 + 0] += sgn*dt * f1 * dpcflux[1] * m[1];
412 J2[0*2 + 0] += sgn*dt * f1 * dpcflux[0] * m[1];
416 *J[ pix[0] ] += dt * (1 - f1) / mt * v1 * dm[0];
419 *J[ pix[1] ] -= dt * f1 / mt * v1 * dm[1];
420 *J[ pix[1] ] += dt * f1 * gflux * dm[1];
423 template <
class Gr
id>
431 const double pv = store_.
porevol(c);
434 *F += pv * store_.
ds(c);
437 template <
class Grid ,
441 const SourceTerms* src,
449 double dflux = -src->flux[i];
453 *F += dt * dflux * src->saturation[2*i + 0];
456 const int c = src->cell[i];
457 const double* m = store_.
mob (c);
458 const double* dm = store_.
dmob(c);
460 const double mt = m[0] + m[1];
462 assert (! ((m[0] < 0) || (m[1] < 0)));
465 const double f = m[0] / mt;
466 const double df = ((1 - f)*dm[0] - f*dm[1]) / mt;
468 *F += dt * dflux * f;
469 *J += dt * dflux * df;
472 template <
class Gr
id>
475 const std::vector<double> & htrans) {
477 assert (htrans.size() ==
478 static_cast<std::vector<double>::size_type
>(g.cell_facepos[ g.number_of_cells ]));
480 for (
int f = 0; f < g.number_of_faces; ++f) {
481 store_.
trans(f) = 0.0;
484 for (
int c = 0, i = 0; c < g.number_of_cells; ++c) {
485 for (; i < g.cell_facepos[c + 1]; ++i) {
486 int f = g.cell_faces[i];
488 assert (htrans[i] > 0.0);
490 store_.
trans(f) += 1.0 / htrans[i];
494 for (
int f = 0; f < g.number_of_faces; ++f) {
499 this->computeStaticGravity(g);
509 class JacobianSystem>
513 JacobianSystem& sys ) {
517 typename JacobianSystem::vector_type& x =
518 sys.vector().writableSolution();
520 assert (x.size() == (::std::size_t) (g.number_of_cells));
522 if (init_step_use_previous_sol_) {
523 std::fill(x.begin(), x.end(), 0.0);
525 const std::vector<double>& s = state.
saturation();
526 for (
int c = 0, nc = g.number_of_cells; c < nc; ++c) {
528 x[c] = 0.5 - s[2*c + 0];
535 class JacobianSystem>
539 JacobianSystem& sys) {
541 double s[2], mob[2], dmob[2 * 2], pc, dpc;
543 const typename JacobianSystem::vector_type& x =
544 sys.vector().solution();
545 const ::std::vector<double>& sat = state.
saturation();
547 bool in_range =
true;
548 for (
int c = 0; c < g.number_of_cells; ++c) {
551 s[0] = sat[c*2 + 0] + x[c];
553 double s_min = fluid_.s_min(c);
554 double s_max = fluid_.s_max(c);
556 if ( s[0] < (s_min - sat_tol_) || s[0] > (s_max + sat_tol_) ) {
565 s[0] = std::max(s_min, s[0]);
569 fluid_.mobility(c, s, mob, dmob);
570 fluid_.pc(c, s, pc, dpc);
572 store_.
mob (c)[0] = mob [0];
573 store_.
mob (c)[1] = mob [1];
574 store_.
dmob(c)[0] = dmob[0*2 + 0];
575 store_.
dmob(c)[1] = -dmob[1*2 + 1];
581 std::cout <<
"Warning: initIteration() - s was clamped in some cells.\n";
589 class NewtonIterate >
593 NewtonIterate& it ) {
595 (void) state; (void) g; (void) it;
598 template <
class Grid ,
599 class SolutionVector,
603 const SolutionVector& x ,
608 for (
int c = 0; c < g.number_of_cells; ++c, s += 2) {
610 double s_min = fluid_.s_min(c);
611 double s_max = fluid_.s_max(c);
614 assert(s[0] >= s_min - sat_tol_);
615 assert(s[0] <= s_max + sat_tol_);
618 s[0] = std::max(s_min, s[0]);
626 upwindMobility(
const double dflux,
632 bool equal_sign = ( (! (dflux < 0)) && (! (gflux < 0)) ) ||
633 ( (! (dflux > 0)) && (! (gflux > 0)) );
637 if (! (dflux < 0) && ! (gflux < 0)) { pix[0] = 0; }
640 m[0] = store_.
mob(n[ pix[0] ]) [ 0 ];
642 if (! (dflux - m[0]*gflux < 0)) { pix[1] = 0; }
645 m[1] = store_.
mob(n[ pix[1] ]) [ 1 ];
649 if (! (dflux < 0) && ! (gflux > 0)) { pix[1] = 0; }
652 m[1] = store_.
mob(n[ pix[1] ]) [ 1 ];
654 if (dflux + m[1]*gflux > 0) { pix[0] = 0; }
657 m[0] = store_.
mob(n[ pix[0] ]) [ 0 ];
660 dm[0] = store_.
dmob(n[ pix[0] ]) [ 0 ];
661 dm[1] = store_.
dmob(n[ pix[1] ]) [ 1 ];
664 template <
class Gr
id>
666 computeStaticGravity(
const Grid& g) {
668 const int d = g.dimensions;
670 for (
int c = 0, i = 0; c < g.number_of_cells; ++c) {
671 const double* cc = g.cell_centroids + (c * d);
673 for (; i < g.cell_facepos[c + 1]; ++i) {
674 const int f = g.cell_faces[i];
675 const double* fc = g.face_centroids + (f * d);
678 for (
int j = 0;
j < d; ++
j) {
679 dg += gravity_[
j] * (fc[
j] - cc[
j]);
682 store_.
dg(i) = store_.
trans(f) * dg;
688 gravityFlux(
const int f)
const {
692 int i1 = f2hf_[2*f + 0];
693 int i2 = f2hf_[2*f + 1];
695 assert ((i1 >= 0) && (i2 >= 0));
697 gflux = store_.
dg(i1) - store_.
dg(i2);
698 gflux *= store_.
drho();
706 capFlux(
const int f,
const int* n,
double& pcflux,
double* dpcflux)
const {
710 assert ((i1 >= 0) && (i2 >= 0));
712 pcflux = store_.
trans(f)*(store_.
pc(i2) - store_.
pc(i1));
713 dpcflux[0] = -store_.
trans(f)*store_.
dpc(i1);
714 dpcflux[1] = store_.
trans(f)*store_.
dpc(i2);
717 TwophaseFluid fluid_ ;
718 const double* gravity_;
719 std::vector<int> f2hf_ ;
720 spu_2p::ModelParameterStorage store_ ;
721 bool init_step_use_previous_sol_;
Definition: ImplicitTransportDefs.hpp:52
::std::vector< double > & faceflux()
Definition: ImplicitTransportDefs.hpp:63
::std::vector< double > & saturation()
Definition: ImplicitTransportDefs.hpp:64
Definition: SinglePointUpwindTwoPhase.hpp:278
void initStep(const ReservoirState &state, const Grid &g, JacobianSystem &sys)
Definition: SinglePointUpwindTwoPhase.hpp:511
void accumulation(const Grid &g, const int c, double *J, double *F) const
Definition: SinglePointUpwindTwoPhase.hpp:425
void initGravityTrans(const Grid &g, const std::vector< double > &htrans)
Definition: SinglePointUpwindTwoPhase.hpp:474
void initResidual(const int c, double *F) const
Definition: SinglePointUpwindTwoPhase.hpp:357
void sourceTerms(const Grid &g, const SourceTerms *src, const int i, const double dt, double *J, double *F) const
Definition: SinglePointUpwindTwoPhase.hpp:440
SinglePointUpwindTwoPhase(const TwophaseFluid &fluid, const Grid &g, const std::vector< double > &porevol, const double *grav=0, const bool guess_previous=true)
Definition: SinglePointUpwindTwoPhase.hpp:281
void makefhfQPeriodic(const std::vector< int > &p_faces, const std::vector< int > &hf_faces, const std::vector< int > &nb_faces)
Definition: SinglePointUpwindTwoPhase.hpp:311
void fluxConnection(const ReservoirState &state, const Grid &g, const double dt, const int c, const int f, double *J1, double *J2, double *F) const
Definition: SinglePointUpwindTwoPhase.hpp:365
void finishStep(const Grid &g, const SolutionVector &x, ReservoirState &state)
Definition: SinglePointUpwindTwoPhase.hpp:602
void finishIteration(const ReservoirState &state, const Grid &g, NewtonIterate &it)
Definition: SinglePointUpwindTwoPhase.hpp:591
@ DofPerCell
Definition: SinglePointUpwindTwoPhase.hpp:354
bool initIteration(const ReservoirState &state, const Grid &g, JacobianSystem &sys)
Definition: SinglePointUpwindTwoPhase.hpp:537
Definition: SinglePointUpwindTwoPhase.hpp:82
double * porevol()
Definition: SinglePointUpwindTwoPhase.hpp:166
ModelParameterStorage(int nc, int totconn)
Definition: SinglePointUpwindTwoPhase.hpp:91
double & drho()
Definition: SinglePointUpwindTwoPhase.hpp:123
double drho() const
Definition: SinglePointUpwindTwoPhase.hpp:128
double trans(int f) const
Definition: SinglePointUpwindTwoPhase.hpp:256
double dg(int i) const
Definition: SinglePointUpwindTwoPhase.hpp:190
double * mob(int c)
Definition: SinglePointUpwindTwoPhase.hpp:136
double & pc(int c)
Definition: SinglePointUpwindTwoPhase.hpp:215
double & dpc(int c)
Definition: SinglePointUpwindTwoPhase.hpp:231
double & trans(int f)
Definition: SinglePointUpwindTwoPhase.hpp:248
double & ds(int c)
Definition: SinglePointUpwindTwoPhase.hpp:199
double pc(int c) const
Definition: SinglePointUpwindTwoPhase.hpp:222
double dpc(int c) const
Definition: SinglePointUpwindTwoPhase.hpp:239
double ds(int c) const
Definition: SinglePointUpwindTwoPhase.hpp:207
const double * mob(int c) const
Definition: SinglePointUpwindTwoPhase.hpp:143
double porevol(int c) const
Definition: SinglePointUpwindTwoPhase.hpp:172
double & dg(int i)
Definition: SinglePointUpwindTwoPhase.hpp:182
const double * dmob(int c) const
Definition: SinglePointUpwindTwoPhase.hpp:160
double * dmob(int c)
Definition: SinglePointUpwindTwoPhase.hpp:152
min[0]
Definition: elasticity_upscale_impl.hpp:146
int j
Definition: elasticity_upscale_impl.hpp:658
Definition: ImplicitAssembly.hpp:43