SinglePointUpwindTwoPhasePolymer.hpp
Go to the documentation of this file.
1/*===========================================================================
2//
3// File: SinglePointUpwindTwoPhase.hpp
4//
5// Created: 2011-09-28 14:21:34+0200
6//
7// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
8// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
9// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
10// Atgeirr F. Rasmussen <atgeirr@sintef.no>
11// Bård Skaflestad <Bard.Skaflestad@sintef.no>
12//
13//==========================================================================*/
14
15
16/*
17 Copyright 2011 SINTEF ICT, Applied Mathematics.
18 Copyright 2011 Statoil ASA.
19
20 This file is part of the Open Porous Media Project (OPM).
21
22 OPM is free software: you can redistribute it and/or modify
23 it under the terms of the GNU General Public License as published by
24 the Free Software Foundation, either version 3 of the License, or
25 (at your option) any later version.
26
27 OPM is distributed in the hope that it will be useful,
28 but WITHOUT ANY WARRANTY; without even the implied warranty of
29 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 GNU General Public License for more details.
31
32 You should have received a copy of the GNU General Public License
33 along with OPM. If not, see <http://www.gnu.org/licenses/>.
34*/
35
36#ifndef OPM_SINGLEPOINTUPWINDTWOPHASEPOLYMER_HPP_HEADER
37#define OPM_SINGLEPOINTUPWINDTWOPHASEPOLYMER_HPP_HEADER
38
39#include <cassert>
40#include <cstddef>
41
42#include <algorithm>
43#include <vector>
44#include <iostream>
45
46namespace Opm {
47 namespace polymer_reorder {
48 class ModelParameterStorage {
49 public:
50 ModelParameterStorage(int nc, int totconn)
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),
55 trans_(0), data_()
56 {
57 size_t alloc_sz;
58
59 alloc_sz = 2 * nc; // mob_
60 alloc_sz += 2 * nc; // dmobds_
61 alloc_sz += nc; // dmobwatdc_
62 alloc_sz += nc; // mc_
63 alloc_sz += nc; // dmcdc_
64 alloc_sz += 1 * nc; // porevol_
65 alloc_sz += 1 * nc; // porosity_
66 alloc_sz += 1 * totconn; // dg_
67 alloc_sz += 1 * nc; // sw_
68 alloc_sz += 1 * nc; // c_
69 alloc_sz += 1 * nc; // cmax_
70 alloc_sz += 1 * nc; // ds_
71 alloc_sz += 1 * nc; // dsc_
72 alloc_sz += 1 * nc; // dcads_
73 alloc_sz += 1 * nc; // dcadsdc_
74 alloc_sz += 1 * nc; // pc_
75 alloc_sz += 1 * nc; // dpc_
76 alloc_sz += 1 * totconn; // trans_
77 data_.resize(alloc_sz);
78
79 mob_ = &data_[0];
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);
88 c_ = sw_ + (1 * nc );
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 );
97 }
98
99 double& drho () { return drho_ ; }
100 double drho () const { return drho_ ; }
101
102 double& rockdensity() { return rockdensity_ ; }
103 double rockdensity() const { return rockdensity_ ; }
104
105 double* mob (int cell) { return mob_ + (2*cell + 0); }
106 const double* mob (int cell) const { return mob_ + (2*cell + 0); }
107
108 double* dmobds (int cell) { return dmobds_ + (2*cell + 0); }
109 const double* dmobds (int cell) const { return dmobds_ + (2*cell + 0); }
110
111 double& dmobwatdc (int cell) { return dmobwatdc_[cell]; }
112 double dmobwatdc (int cell) const { return dmobwatdc_[cell]; }
113
114 double& mc (int cell) { return mc_[cell]; }
115 double mc (int cell) const { return mc_[cell]; }
116
117 double& dmcdc (int cell) { return dmcdc_[cell]; }
118 double dmcdc (int cell) const { return dmcdc_[cell]; }
119
120 double* porevol() { return porevol_ ; }
121 double porevol(int cell) const { return porevol_[cell] ; }
122
123 double* porosity() { return porosity_ ; }
124 double porosity(int cell) const { return porosity_[cell] ; }
125
126
127 double& dg(int i) { return dg_[i] ; }
128 double dg(int i) const { return dg_[i] ; }
129
130 double& sw(int cell) { return sw_[cell] ; }
131 double sw(int cell) const { return sw_[cell] ; }
132
133 double& c(int cell) { return c_[cell] ; }
134 double c(int cell) const { return c_[cell] ; }
135
136 double& cmax(int cell) { return cmax_[cell] ; }
137 double cmax(int cell) const { return cmax_[cell] ; }
138
139 double& ds(int cell) { return ds_[cell] ; }
140 double ds(int cell) const { return ds_[cell] ; }
141
142 double& dsc(int cell) { return dsc_[cell] ; }
143 double dsc(int cell) const { return dsc_[cell] ; }
144
145 double& dcads(int cell) { return dcads_[cell] ; }
146 double dcads(int cell) const { return dcads_[cell] ; }
147
148 double& dcadsdc(int cell) { return dcadsdc_[cell] ; }
149 double dcadsdc(int cell) const { return dcadsdc_[cell] ; }
150
151 double& pc(int cell) { return pc_[cell] ; }
152 double pc(int cell) const { return pc_[cell] ; }
153
154 double& dpc(int cell) { return dpc_[cell] ; }
155 double dpc(int cell) const { return dpc_[cell] ; }
156
157 double& trans(int f) { return trans_[f] ; }
158 double trans(int f) const { return trans_[f] ; }
159
160 private:
161 double drho_ ;
162 double rockdensity_ ;
163 double *mob_ ;
164 double *dmobds_ ;
165 double *dmobwatdc_ ;
166 double *mc_ ;
167 double *dmcdc_ ;
168 double *porevol_ ;
169 double *porosity_ ;
170 double *dg_ ;
171 double *sw_ ;
172 double *c_ ;
173 double *cmax_ ;
174 double *ds_ ;
175 double *dsc_ ;
176 double *dcads_ ; // difference of cads to compute residual
177 double *dcadsdc_ ; // derivative of cads
178 double *pc_ ;
179 double *dpc_ ;
180 double *trans_ ;
181
182 std::vector<double> data_;
183 };
184 }
185
186
187 template <class TwophaseFluidPolymer>
188 class SinglePointUpwindTwoPhasePolymer {
189 public:
190 template <class Grid>
191 SinglePointUpwindTwoPhasePolymer(const TwophaseFluidPolymer& fluid ,
192 const Grid& g ,
193 const std::vector<double>& porevol ,
194 const double* grav = 0,
195 const bool guess_previous = true)
196 : fluid_ (fluid) ,
197 gravity_ (grav) ,
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) ,
202 sat_tol_ (1e-5)
203 {
204
205 if (gravity_) {
206 store_.drho() = fluid_.density(0) - fluid_.density(1);
207 }
208
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);
213 f2hf_[2*f + p] = i;
214 }
215 }
216
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();
221 }
222
223 void makefhfQPeriodic( const std::vector<int>& p_faces,const std::vector<int>& hf_faces,
224 const std::vector<int>& nb_faces)
225 {
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];
231 }else{
232 assert(f2hf_[2*nbf+1]==-1);
233 nbhf[i] = f2hf_[2*nbf];
234 }
235 }
236 for(unsigned int i=0; i<p_faces.size(); ++i){
237
238 int f = p_faces[i];
239 int hf = hf_faces[i];
240 bool changed=false;
241
242 if(f2hf_[2*f] == hf){
243 assert(f2hf_[2*f+1]==-1);
244 }else{
245 assert(f2hf_[2*f]==-1);
246 f2hf_[2*f]=nbhf[i];
247 changed=true;
248 }
249 if(!changed){
250 if(f2hf_[2*f+1]== hf){
251 assert(f2hf_[2*f]==-1);
252 }else{
253 assert(f2hf_[2*f+1]==-1);
254 f2hf_[2*f+1]=nbhf[i];
255 changed=true;
256 }
257 }
258 assert(changed);
259 }
260 }
261
262 // -----------------------------------------------------------------
263 // System assembly innards
264 // -----------------------------------------------------------------
265
266 enum { DofPerCell = 1 };
267
268 void
269 initResidual(const int c, double* Fs, double* Fc) const {
270 (void) c; // Suppress 'unused' warning
271 *Fs = 0.0;
272 *Fc = 0.0;
273 }
274
275 template <class ReservoirState,
276 class Grid >
277 void
278 fluxConnection(const ReservoirState& state ,
279 const Grid& g ,
280 const double dt ,
281 const int cell ,
282 const int f ,
283 double* F , // F[0] = s-residual, F[1] = c-residual
284 double* dFd1 , //Jacobi matrix for residual with respect to variables in cell
285 double* dFd2 //Jacobi matrix for residual with respect to variables in OTHER cell
286 //dFd1[0]= d(F[0])/d(s1), dFd1[1]= d(F[0])/d(c1), dFd1[2]= d(F[1])/d(s1), dFd1[3]= d(F[1])/d(c1),
287 //dFd2[0]= d(F[0])/d(s2), dFd2[1]= d(F[0])/d(c2), dFd2[2]= d(F[1])/d(s2), dFd2[3]= d(F[1])/d(c2).
288
289 ) const {
290
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);
296 gflux += pcflux;
297
298 int pix[2];
299 double m[2], dmds[2], dmobwatdc;
300 double mc, dmcdc;
301 upwindMobility(dflux, gflux, n, pix, m, dmds, dmobwatdc, mc, dmcdc);
302
303 assert ((m[0] >= 0.0) && (m[1] >= 0.0));
304
305 double mt = m[0] + m[1];
306 assert (mt >= 0.0);
307
308 double sgn = 2.0*(n[0] == cell) - 1.0;
309 dflux *= sgn;
310 gflux *= sgn;
311
312
313 double f1 = m[0] / mt;
314 const double v1 = dflux + m[1]*gflux;
315
316 // Assemble residual contributions
317 F[0] += dt * f1 * v1;
318 F[1] += dt * mc * f1 * v1;
319
320 // Assemble Jacobian (J1 <-> cell, J2 <-> other)
321 double *dFsds[2];
322 double *dFsdc[2];
323 double *dFcds[2];
324 double *dFcdc[2];
325 if (n[0] == cell) {
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];
330 // sign is positive
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];
335 // We assume that the capillary pressure is independent of the polymer concentration.
336 // Hence, no more contributions.
337 } else {
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];
342 // sign is negative
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];
347 // We assume that the capillary pressure is independent of the polymer concentration.
348 // Hence, no more contributions.
349 }
350
351 // dFs/dm_1 \cdot dm_1/ds
352 *dFsds[ pix[0] ] += dt * (1 - f1) / mt * v1 * dmds[0];
353 // dFc/dm_1 \cdot dm_1/ds
354 *dFcds[ pix[0] ] += dt * (1 - f1) / mt * v1 * mc * dmds[0];
355
356 // dFs/dm_2 \cdot dm_2/ds
357 *dFsds[ pix[1] ] -= dt * f1 / mt * v1 * dmds[1];
358 *dFsds[ pix[1] ] += dt * f1 * gflux * dmds[1];
359 // dFc/dm_2 \cdot dm_2/ds
360 *dFcds[ pix[1] ] -= dt * f1 / mt * v1 * mc * dmds[1];
361 *dFcds[ pix[1] ] += dt * f1 * gflux * mc * dmds[1];
362
363 // dFs/dm_1 \cdot dm_1/dc
364 *dFsdc[ pix[0] ] += dt * (1 - f1) / mt * v1 * dmobwatdc;
365 // dFc/dm_1 \cdot dm_1/dc
366 *dFcdc[ pix[0] ] += dt * (1 - f1) / mt * v1 * mc * dmobwatdc;
367 *dFcdc[ pix[0] ] += dt * f1 * v1 * dmcdc; // Polymer is only carried by water.
368 }
369
370 template <class Grid>
371 void
372 accumulation(const Grid& g,
373 const int cell,
374 double* F, // Residual vector,
375 double* dF // Jacobian, same convention as for fluxConnection.
376 ) const {
377 (void) g;
378
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);
383
384 F[0] += pv * store_.ds(cell);
385 F[1] += pv * (1 - dps) * store_.dsc(cell) + rhor*(1 - poro)/poro*pv*store_.dcads(cell);
386 dF[0] += pv;
387 dF[1] += 0.;
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);
390 }
391
392 template <class Grid ,
393 class SourceTerms>
394 void
395 sourceTerms(const Grid& g ,
396 const SourceTerms* src,
397 const int i ,
398 const double dt ,
399 double* J ,
400 double* F ) const {
401
402 (void) g;
403
404 double dflux = -src->flux[i]; // ->flux[] is rate of *inflow*
405
406 if (dflux < 0) {
407 // src -> cell, affects residual only.
408 *F += dt * dflux * src->saturation[2*i + 0];
409 } else {
410 // cell -> src
411 const int cell = src->cell[i];
412 const double* m = store_.mob (cell);
413 const double* dm = store_.dmobds(cell);
414
415 const double mt = m[0] + m[1];
416
417 assert (! ((m[0] < 0) || (m[1] < 0)));
418 assert (mt > 0);
419
420 const double f = m[0] / mt;
421 const double df = ((1 - f)*dm[0] - f*dm[1]) / mt;
422
423 *F += dt * dflux * f;
424 *J += dt * dflux * df;
425 }
426 }
427 template <class Grid>
428 void
429 initGravityTrans(const Grid& g ,
430 const std::vector<double> & htrans) {
431
432 assert (htrans.size() ==
433 static_cast<std::vector<double>::size_type>(g.cell_facepos[ g.number_of_cells ]));
434
435 for (int f = 0; f < g.number_of_faces; ++f) {
436 store_.trans(f) = 0.0;
437 }
438
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];
442
443 assert (htrans[i] > 0.0);
444
445 store_.trans(f) += 1.0 / htrans[i];
446 }
447 }
448
449 for (int f = 0; f < g.number_of_faces; ++f) {
450 store_.trans(f) = 1.0 / store_.trans(f);
451 }
452
453 if (gravity_) {
454 this->computeStaticGravity(g);
455 }
456 }
457
458 // -----------------------------------------------------------------
459 // Newton control
460 // -----------------------------------------------------------------
461
462 template <class ReservoirState,
463 class Grid ,
464 class JacobianSystem>
465 void
466 initStep(const ReservoirState& state,
467 const Grid& g ,
468 JacobianSystem& sys ) {
469
470 (void) state; // Suppress 'unused' warning.
471
472 typename JacobianSystem::vector_type& x =
473 sys.vector().writableSolution();
474
475 assert (x.size() == (::std::size_t) (2*g.number_of_cells));
476
477 if (init_step_use_previous_sol_) {
478 std::fill(x.begin(), x.end(), 0.0);
479 } else {
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) {
484 // Impose s=0.5 at next time level as an NR initial value.
485 x[2*cell + 0] = 0.5 - s[2*cell + 0];
486 x[2*cell + 1] = 1e-5 - c[cell];
487 }
488 }
489 }
490
491 template <class ReservoirState,
492 class Grid ,
493 class JacobianSystem>
494 bool
495 initIteration(const ReservoirState& state,
496 const Grid& g ,
497 JacobianSystem& sys) {
498
499 double s[2];
500 double mob[2];
501 double dmobds[4];
502 double dmobwatdc;
503 double c, cmax;
504 double mc, dmcdc;
505 double pc, dpc;
506
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();
512
513 bool in_range = true;
514 for (int cell = 0; cell < g.number_of_cells; ++cell) {
515 // Store wat-sat, sat-change, cpoly, (sat * cpoly)-change for accumulation().
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];
520 store_.c(cell) = c;
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];
524 double dcadsdc;
525 double cads;
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);
533
534 if ( s[0] < (s_min - sat_tol_) || s[0] > (s_max + sat_tol_) ) {
535 // if (s[0] < s_min){
536 // std::cout << "Warning: s out of range, s-s_min = " << s_min-s[0] << std::endl;
537 // }
538 // if (s[0] > s_max){
539 // std::cout << "Warning: s out of range, s-s_max = " << s[0]-s_max << std::endl;
540 // }
541 in_range = false; //line search fails
542 }
543 s[0] = std::max(s_min, s[0]);
544 s[0] = std::min(s_max, s[0]);
545 s[1] = 1 - s[0];
546
547 fluid_.mobility(cell, s, c, cmax, mob, dmobds, dmobwatdc);
548 fluid_.computeMc(c, mc, dmcdc);
549 fluid_.pc(cell, s, pc, dpc);
550
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;
560 }
561 if (!in_range) {
562 std::cout << "Warning: initIteration() - s was clamped in some cells.\n";
563 }
564 return in_range;
565 }
566
567 template <class ReservoirState,
568 class Grid ,
569 class NewtonIterate >
570 void
571 finishIteration(const ReservoirState& state,
572 const Grid& g ,
573 NewtonIterate& it ) {
574 // Nothing to do at end of iteration in this model.
575 (void) state; (void) g; (void) it;
576 typedef typename NewtonIterate::vector_type vector_t;
577 }
578
579 template <class Grid ,
580 class SolutionVector,
581 class ReservoirState>
582 void
583 finishStep(const Grid& g ,
584 const SolutionVector& x ,
585 ReservoirState& state) {
586
587 double *s = &state.saturation()[0*2 + 0];
588 double *c = &state.concentration()[0*1 + 0];
589 double *cmax = &state.maxconcentration()[0*1 + 0];
590
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]);
601 s[1] = 1.0 - s[0];
602 }
603 }
604
605 private:
606 void
607 upwindMobility(const double dflux,
608 const double gflux,
609 const int* n ,
610 int* pix ,
611 double* m ,
612 double* dmds ,
613 double& dmobwatdc ,
614 double& mc,
615 double& dmcdc) const {
616 bool equal_sign = ( (! (dflux < 0)) && (! (gflux < 0)) ) ||
617 ( (! (dflux > 0)) && (! (gflux > 0)) );
618
619 if (equal_sign) {
620
621 if (! (dflux < 0) && ! (gflux < 0)) { pix[0] = 0; }
622 else { pix[0] = 1; }
623
624 m[0] = store_.mob(n[ pix[0] ]) [ 0 ];
625 mc = store_.mc(n[ pix[0] ]);
626
627 if (! (dflux - m[0]*gflux < 0)) { pix[1] = 0; }
628 else { pix[1] = 1; }
629
630 m[1] = store_.mob(n[ pix[1] ]) [ 1 ];
631
632 } else {
633
634 if (! (dflux < 0) && ! (gflux > 0)) { pix[1] = 0; }
635 else { pix[1] = 1; }
636
637 m[1] = store_.mob(n[ pix[1] ]) [ 1 ];
638
639 if (dflux + m[1]*gflux > 0) { pix[0] = 0; }
640 else { pix[0] = 1; }
641
642 m[0] = store_.mob(n[ pix[0] ]) [ 0 ];
643 mc = store_.mc(n[ pix[0] ]);
644 }
645
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] ]);
650 }
651
652 template <class Grid>
653 void
654 computeStaticGravity(const Grid& g) {
655
656 const int d = g.dimensions;
657
658 for (int c = 0, i = 0; c < g.number_of_cells; ++c) {
659 const double* cc = g.cell_centroids + (c * d);
660
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);
664
665 double dg = 0.0;
666 for (int j = 0; j < d; ++j) {
667 dg += gravity_[j] * (fc[j] - cc[j]);
668 }
669
670 store_.dg(i) = store_.trans(f) * dg;
671 }
672 }
673 }
674
675 double
676 gravityFlux(const int f) const {
677 double gflux;
678
679 if (gravity_) {
680 int i1 = f2hf_[2*f + 0];
681 int i2 = f2hf_[2*f + 1];
682
683 assert ((i1 >= 0) && (i2 >= 0));
684
685 gflux = store_.dg(i1) - store_.dg(i2);
686 gflux *= store_.drho();
687 } else {
688 gflux = 0.0;
689 }
690
691 return gflux;
692 }
693 void
694 capFlux(const int f,const int* n,double& pcflux, double* dpcflux) const {
695 //double capflux;
696 int i1 = n[0];
697 int i2 = n[1];
698 assert ((i1 >= 0) && (i2 >= 0));
699 //double sgn=-1.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);
703 }
704
705 TwophaseFluidPolymer fluid_ ;
706 const double* gravity_;
707 std::vector<int> f2hf_ ;
708 polymer_reorder::ModelParameterStorage store_ ;
709 bool init_step_use_previous_sol_;
710 double sat_tol_;
711 };
712}
713#endif /* OPM_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER */
cpoly
std::vector< double > & cpoly
Definition: GravityColumnSolverPolymer_impl.hpp:77
sat
std::vector< double > & sat
Definition: GravityColumnSolverPolymer_impl.hpp:76
cmax
std::vector< double > & cmax
Definition: GravityColumnSolverPolymer_impl.hpp:78
Opm::SinglePointUpwindTwoPhasePolymer
Definition: SinglePointUpwindTwoPhasePolymer.hpp:188
Opm::SinglePointUpwindTwoPhasePolymer::sourceTerms
void sourceTerms(const Grid &g, const SourceTerms *src, const int i, const double dt, double *J, double *F) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:395
Opm::SinglePointUpwindTwoPhasePolymer::initIteration
bool initIteration(const ReservoirState &state, const Grid &g, JacobianSystem &sys)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:495
Opm::SinglePointUpwindTwoPhasePolymer::fluxConnection
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
Opm::SinglePointUpwindTwoPhasePolymer::DofPerCell
@ DofPerCell
Definition: SinglePointUpwindTwoPhasePolymer.hpp:266
Opm::SinglePointUpwindTwoPhasePolymer::makefhfQPeriodic
void makefhfQPeriodic(const std::vector< int > &p_faces, const std::vector< int > &hf_faces, const std::vector< int > &nb_faces)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:223
Opm::SinglePointUpwindTwoPhasePolymer::finishIteration
void finishIteration(const ReservoirState &state, const Grid &g, NewtonIterate &it)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:571
Opm::SinglePointUpwindTwoPhasePolymer::finishStep
void finishStep(const Grid &g, const SolutionVector &x, ReservoirState &state)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:583
Opm::SinglePointUpwindTwoPhasePolymer::initGravityTrans
void initGravityTrans(const Grid &g, const std::vector< double > &htrans)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:429
Opm::SinglePointUpwindTwoPhasePolymer::initStep
void initStep(const ReservoirState &state, const Grid &g, JacobianSystem &sys)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:466
Opm::SinglePointUpwindTwoPhasePolymer::initResidual
void initResidual(const int c, double *Fs, double *Fc) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:269
Opm::SinglePointUpwindTwoPhasePolymer::SinglePointUpwindTwoPhasePolymer
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
Opm::SinglePointUpwindTwoPhasePolymer::accumulation
void accumulation(const Grid &g, const int cell, double *F, double *dF) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:372
Opm::polymer_reorder::ModelParameterStorage
Definition: SinglePointUpwindTwoPhasePolymer.hpp:48
Opm::polymer_reorder::ModelParameterStorage::porosity
double porosity(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:124
Opm::polymer_reorder::ModelParameterStorage::dmobds
double * dmobds(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:108
Opm::polymer_reorder::ModelParameterStorage::porevol
double porevol(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:121
Opm::polymer_reorder::ModelParameterStorage::dsc
double dsc(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:143
Opm::polymer_reorder::ModelParameterStorage::dmcdc
double & dmcdc(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:117
Opm::polymer_reorder::ModelParameterStorage::trans
double & trans(int f)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:157
Opm::polymer_reorder::ModelParameterStorage::drho
double & drho()
Definition: SinglePointUpwindTwoPhasePolymer.hpp:99
Opm::polymer_reorder::ModelParameterStorage::mc
double mc(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:115
Opm::polymer_reorder::ModelParameterStorage::dcads
double & dcads(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:145
Opm::polymer_reorder::ModelParameterStorage::dpc
double & dpc(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:154
Opm::polymer_reorder::ModelParameterStorage::sw
double & sw(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:130
Opm::polymer_reorder::ModelParameterStorage::dmobwatdc
double & dmobwatdc(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:111
Opm::polymer_reorder::ModelParameterStorage::dmobwatdc
double dmobwatdc(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:112
Opm::polymer_reorder::ModelParameterStorage::cmax
double cmax(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:137
Opm::polymer_reorder::ModelParameterStorage::c
double c(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:134
Opm::polymer_reorder::ModelParameterStorage::pc
double & pc(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:151
Opm::polymer_reorder::ModelParameterStorage::cmax
double & cmax(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:136
Opm::polymer_reorder::ModelParameterStorage::pc
double pc(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:152
Opm::polymer_reorder::ModelParameterStorage::ModelParameterStorage
ModelParameterStorage(int nc, int totconn)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:50
Opm::polymer_reorder::ModelParameterStorage::rockdensity
double & rockdensity()
Definition: SinglePointUpwindTwoPhasePolymer.hpp:102
Opm::polymer_reorder::ModelParameterStorage::ds
double & ds(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:139
Opm::polymer_reorder::ModelParameterStorage::dcads
double dcads(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:146
Opm::polymer_reorder::ModelParameterStorage::dg
double & dg(int i)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:127
Opm::polymer_reorder::ModelParameterStorage::dcadsdc
double & dcadsdc(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:148
Opm::polymer_reorder::ModelParameterStorage::drho
double drho() const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:100
Opm::polymer_reorder::ModelParameterStorage::sw
double sw(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:131
Opm::polymer_reorder::ModelParameterStorage::mob
double * mob(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:105
Opm::polymer_reorder::ModelParameterStorage::dg
double dg(int i) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:128
Opm::polymer_reorder::ModelParameterStorage::dsc
double & dsc(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:142
Opm::polymer_reorder::ModelParameterStorage::rockdensity
double rockdensity() const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:103
Opm::polymer_reorder::ModelParameterStorage::mob
const double * mob(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:106
Opm::polymer_reorder::ModelParameterStorage::trans
double trans(int f) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:158
Opm::polymer_reorder::ModelParameterStorage::mc
double & mc(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:114
Opm::polymer_reorder::ModelParameterStorage::ds
double ds(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:140
Opm::polymer_reorder::ModelParameterStorage::dmcdc
double dmcdc(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:118
Opm::polymer_reorder::ModelParameterStorage::porevol
double * porevol()
Definition: SinglePointUpwindTwoPhasePolymer.hpp:120
Opm::polymer_reorder::ModelParameterStorage::porosity
double * porosity()
Definition: SinglePointUpwindTwoPhasePolymer.hpp:123
Opm::polymer_reorder::ModelParameterStorage::dmobds
const double * dmobds(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:109
Opm::polymer_reorder::ModelParameterStorage::c
double & c(int cell)
Definition: SinglePointUpwindTwoPhasePolymer.hpp:133
Opm::polymer_reorder::ModelParameterStorage::dpc
double dpc(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:155
Opm::polymer_reorder::ModelParameterStorage::dcadsdc
double dcadsdc(int cell) const
Definition: SinglePointUpwindTwoPhasePolymer.hpp:149
Opm
Definition: CompressibleTpfaPolymer.hpp:33