GEOS
IHUPhaseFlux.hpp
Go to the documentation of this file.
1 /*
2  * ------------------------------------------------------------------------------------------------------------
3  * SPDX-License-Identifier: LGPL-2.1-only
4  *
5  * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
6  * Copyright (c) 2018-2024 Total, S.A
7  * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
8  * Copyright (c) 2023-2024 Chevron
9  * Copyright (c) 2019- GEOS/GEOSX Contributors
10  * All rights reserved
11  *
12  * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
13  * ------------------------------------------------------------------------------------------------------------
14  */
15 
20 #ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_IHUPHASEFLUX_HPP
21 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_IHUPHASEFLUX_HPP
22 
23 #include "common/DataLayouts.hpp"
24 #include "common/DataTypes.hpp"
25 #include "constitutive/fluid/multifluid/Layouts.hpp"
26 #include "constitutive/capillaryPressure/layouts.hpp"
27 #include "mesh/ElementRegionManager.hpp"
28 
29 
30 namespace geos
31 {
32 
33 namespace isothermalCompositionalMultiphaseFVMKernelUtilities
34 {
35 
36 template< typename VIEWTYPE >
37 using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
38 
39 using Deriv = constitutive::multifluid::DerivativeOffset;
40 
41 /************************* HELPERS ******************/
42 namespace UpwindHelpers
43 {
44 
45 static constexpr double minTotMob = 1e-12;
46 
47 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
48 GEOS_HOST_DEVICE
49 static void
50 upwindMobilityViscous( localIndex const numPhase,
51  localIndex const ip,
52  localIndex const (&seri)[numFluxSupportPoints],
53  localIndex const (&sesri)[numFluxSupportPoints],
54  localIndex const (&sei)[numFluxSupportPoints],
55  real64 const (&transmissibility)[2],
56  real64 const (&dTrans_dPres)[2],
57  real64 const totFlux, //in fine should be a ElemnetViewConst once seq form are in place
58  ElementViewConst< arrayView1d< real64 const > > const & pres,
59  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
60  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
61  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
62  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
63  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
64  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
65  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
66  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
67  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
68  integer const capPressureFlag,
69  localIndex & upwindDir,
70  real64 & mobility,
71  real64( &dMobility_dP),
72  real64 ( & dMobility_dC)[numComp]
73  )
74 {
75 
76  //reinit
77  mobility = 0.0;
78  dMobility_dP = 0.0;
79  for( localIndex ic = 0; ic < numComp; ++ic )
80  {
81  dMobility_dC[ic] = 0.0;
82  }
83 
84  UPWIND scheme;
85  scheme.template getUpwindDirectionViscous< numComp, numFluxSupportPoints, UPWIND >( numPhase,
86  ip,
87  seri,
88  sesri,
89  sei,
90  transmissibility,
91  dTrans_dPres,
92  totFlux,
93  pres,
94  gravCoef,
95  phaseMob,
96  dCompFrac_dCompDens,
97  phaseMassDens,
98  dPhaseMassDens,
99  dPhaseVolFrac,
100  phaseCapPressure,
101  dPhaseCapPressure_dPhaseVolFrac,
102  capPressureFlag,
103  upwindDir );
104 
105  localIndex const er_up = seri[upwindDir];
106  localIndex const esr_up = sesri[upwindDir];
107  localIndex const ei_up = sei[upwindDir];
108 
109  if( std::fabs( phaseMob[er_up][esr_up][ei_up][ip] ) > 1e-20 )
110  {
111  mobility = phaseMob[er_up][esr_up][ei_up][ip];
112  dMobility_dP = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dP];
113  for( localIndex ic = 0; ic < numComp; ++ic )
114  {
115  dMobility_dC[ic] = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dC + ic];
116  }
117  }
118 }
119 
120 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
121 GEOS_HOST_DEVICE
122 static void
123 upwindMobilityGravity( localIndex const numPhase,
124  localIndex const ip,
125  localIndex const (&seri)[numFluxSupportPoints],
126  localIndex const (&sesri)[numFluxSupportPoints],
127  localIndex const (&sei)[numFluxSupportPoints],
128  real64 const (&transmissibility)[2],
129  real64 const (&dTrans_dPres)[2],
130  real64 const totFlux, //in fine should be a ElemnetViewConst once seq form are in place
131  ElementViewConst< arrayView1d< real64 const > > const & pres,
132  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
133  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
134  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
135  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
136  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
137  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
138  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
139  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
140  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
141  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
142  integer const capPressureFlag,
143  localIndex & upwindDir,
144  real64 & mobility,
145  real64( &dMobility_dP),
146  real64 ( & dMobility_dC)[numComp]
147  )
148 {
149 
150  //reinit
151  mobility = 0.0;
152  dMobility_dP = 0.0;
153  for( localIndex ic = 0; ic < numComp; ++ic )
154  {
155  dMobility_dC[ic] = 0.0;
156  }
157 
158  UPWIND scheme;
159  scheme.template getUpwindDirectionGravity< numComp, numFluxSupportPoints, UPWIND >( numPhase,
160  ip,
161  seri,
162  sesri,
163  sei,
164  transmissibility,
165  dTrans_dPres,
166  totFlux,
167  pres,
168  gravCoef,
169  phaseMob,
170  dCompFrac_dCompDens,
171  phaseMassDens,
172  dPhaseMassDens,
173  phaseVolFrac,
174  dPhaseVolFrac,
175  phaseCapPressure,
176  dPhaseCapPressure_dPhaseVolFrac,
177  capPressureFlag,
178  upwindDir );
179 
180  localIndex const er_up = seri[upwindDir];
181  localIndex const esr_up = sesri[upwindDir];
182  localIndex const ei_up = sei[upwindDir];
183 
184  if( std::fabs( phaseMob[er_up][esr_up][ei_up][ip] ) > 1e-20 )
185  {
186  mobility = phaseMob[er_up][esr_up][ei_up][ip];
187  dMobility_dP = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dP];
188  for( localIndex ic = 0; ic < numComp; ++ic )
189  {
190  dMobility_dC[ic] = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dC + ic];
191  }
192  }
193 }
194 
195 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
196 GEOS_HOST_DEVICE
197 static void
198 upwindMobilityCapillary( localIndex const numPhase,
199  localIndex const ip,
200  localIndex const (&seri)[numFluxSupportPoints],
201  localIndex const (&sesri)[numFluxSupportPoints],
202  localIndex const (&sei)[numFluxSupportPoints],
203  real64 const (&transmissibility)[2],
204  real64 const (&dTrans_dPres)[2],
205  real64 const totFlux, //in fine should be a ElemnetViewConst once seq form are in place
206  ElementViewConst< arrayView1d< real64 const > > const & pres,
207  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
208  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
209  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
210  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
211  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
212  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
213  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
214  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
215  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
216  integer const capPressureFlag,
217  localIndex & upwindDir,
218  real64 & mobility,
219  real64( &dMobility_dP),
220  real64 ( & dMobility_dC)[numComp]
221  )
222 {
223 
224  //reinit
225  mobility = 0.0;
226  dMobility_dP = 0.0;
227  for( localIndex ic = 0; ic < numComp; ++ic )
228  {
229  dMobility_dC[ic] = 0.0;
230  }
231 
232  UPWIND scheme;
233  scheme.template getUpwindDirectionCapillary< numComp, numFluxSupportPoints, UPWIND >( numPhase,
234  ip,
235  seri,
236  sesri,
237  sei,
238  transmissibility,
239  dTrans_dPres,
240  totFlux,
241  pres,
242  gravCoef,
243  phaseMob,
244  dCompFrac_dCompDens,
245  phaseMassDens,
246  dPhaseMassDens,
247  dPhaseVolFrac,
248  phaseCapPressure,
249  dPhaseCapPressure_dPhaseVolFrac,
250  capPressureFlag,
251  upwindDir );
252 
253  localIndex const er_up = seri[upwindDir];
254  localIndex const esr_up = sesri[upwindDir];
255  localIndex const ei_up = sei[upwindDir];
256 
257  if( std::fabs( phaseMob[er_up][esr_up][ei_up][ip] ) > 1e-20 )
258  {
259  mobility = phaseMob[er_up][esr_up][ei_up][ip];
260  dMobility_dP = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dP];
261  for( localIndex ic = 0; ic < numComp; ++ic )
262  {
263  dMobility_dC[ic] = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dC + ic];
264  }
265  }
266 }
267 
268 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
269 GEOS_HOST_DEVICE
270 static void
271 computeFractionalFlowViscous( localIndex const numPhase,
272  localIndex const ip,
273  localIndex const (&seri)[numFluxSupportPoints],
274  localIndex const (&sesri)[numFluxSupportPoints],
275  localIndex const (&sei)[numFluxSupportPoints],
276  real64 const (&transmissibility)[2],
277  real64 const (&dTrans_dPres)[2],
278  localIndex const & k_up_ppu,
279  real64 const totFlux,
280  real64 const totMob,
281  real64 const (&dTotMob_dP)[numFluxSupportPoints],
282  real64 const (&dTotMob_dC)[numFluxSupportPoints][numComp],
283  ElementViewConst< arrayView1d< real64 const > > const & pres,
284  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
285  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
286  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
287  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
288  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
289  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
290  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
291  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
292  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
293  integer const capPressureFlag,
294  localIndex & k_up_main,
295  real64 & fractionalFlow,
296  real64 ( & dFractionalFlow_dP)[numFluxSupportPoints],
297  real64 ( & dFractionalFlow_dC)[numFluxSupportPoints][numComp] )
298 {
299  // get var to memorized the numerator mobility properly upwinded
300  real64 mainMob{};
301  real64 dMMob_dP{};
302  real64 dMMob_dC[numComp]{};
303 
304 // real64 totMob{};
305 // real64 dTotMob_dP[numFluxSupportPoints]{};
306 // real64 dTotMob_dC[numFluxSupportPoints][numComp]{};
307 
308  //reinit
309  //fractional flow too low to let the upstream phase flow
310  k_up_main = -1; //to throw error if unmodified
311  fractionalFlow = 0;
312  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
313  {
314  dFractionalFlow_dP[ke] = 0;
315  for( localIndex jc = 0; jc < numComp; ++jc )
316  {
317  dFractionalFlow_dC[ke][jc] = 0;
318  }
319  }
320 
321  localIndex k_up;
322  real64 mob{};
323  real64 dMob_dP{};
324  real64 dMob_dC[numComp]{};
325 
326  upwindMobilityViscous< numComp, numFluxSupportPoints, UPWIND >( numPhase,
327  ip,
328  seri,
329  sesri,
330  sei,
331  transmissibility,
332  dTrans_dPres,
333  totFlux,
334  pres,
335  gravCoef,
336  dCompFrac_dCompDens,
337  phaseMassDens,
338  dPhaseMassDens,
339  phaseMob,
340  dPhaseMob,
341  dPhaseVolFrac,
342  phaseCapPressure,
343  dPhaseCapPressure_dPhaseVolFrac,
344  capPressureFlag,
345  k_up,
346  mob,
347  dMob_dP,
348  dMob_dC );
349 
350  k_up_main = k_up;
351  mainMob = mob;
352  dMMob_dP = dMob_dP;
353  for( localIndex ic = 0; ic < numComp; ++ic )
354  {
355  dMMob_dC[ic] = dMob_dC[ic];
356  }
357 
358  //guard against no flow region
359  if( std::fabs( mainMob ) > 1e-20 )
360  {
361  fractionalFlow = mainMob / LvArray::math::max( totMob, minTotMob );
362  dFractionalFlow_dP[k_up_main] = dMMob_dP / LvArray::math::max( totMob, minTotMob );
363  for( localIndex jc = 0; jc < numComp; ++jc )
364  {
365  dFractionalFlow_dC[k_up_main][jc] = dMMob_dC[jc] / totMob;
366 
367  }
368 
369  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
370  {
371  dFractionalFlow_dP[ke] -= fractionalFlow * dTotMob_dP[k_up_ppu] / LvArray::math::max( totMob, minTotMob );
372 
373  for( localIndex jc = 0; jc < numComp; ++jc )
374  {
375  dFractionalFlow_dC[ke][jc] -= fractionalFlow * dTotMob_dC[k_up_ppu][jc] / LvArray::math::max( totMob, minTotMob );
376  }
377  }
378  }
379 }
380 
381 
382 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
383 GEOS_HOST_DEVICE
384 static void
385 computeFractionalFlowGravity( localIndex const numPhase,
386  localIndex const ip,
387  localIndex const (&seri)[numFluxSupportPoints],
388  localIndex const (&sesri)[numFluxSupportPoints],
389  localIndex const (&sei)[numFluxSupportPoints],
390  real64 const (&transmissibility)[2],
391  real64 const (&dTrans_dPres)[2],
392  localIndex const & k_up_ppu,
393  real64 const totFlux,
394  real64 const totMob,
395  real64 const (&dTotMob_dP)[numFluxSupportPoints],
396  real64 const (&dTotMob_dC)[numFluxSupportPoints][numComp],
397  ElementViewConst< arrayView1d< real64 const > > const & pres,
398  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
399  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
400  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
401  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
402  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
403  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
404  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
405  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
406  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
407  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
408  integer const capPressureFlag,
409  localIndex & k_up_main,
410  real64 & fractionalFlow,
411  real64 ( & dFractionalFlow_dP)[numFluxSupportPoints],
412  real64 ( & dFractionalFlow_dC)[numFluxSupportPoints][numComp]
413  )
414 {
415  // get var to memorized the numerator mobility properly upwinded
416  real64 mainMob{};
417  real64 dMMob_dP{};
418  real64 dMMob_dC[numComp]{};
419 
420  //reinit
421  //fractional flow too low to let the upstream phase flow
422  k_up_main = -1; //to throw error if unmodified
423  fractionalFlow = 0;
424  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
425  {
426  dFractionalFlow_dP[ke] = 0;
427  for( localIndex jc = 0; jc < numComp; ++jc )
428  {
429  dFractionalFlow_dC[ke][jc] = 0;
430  }
431  }
432 
433 
434  localIndex k_up;
435  real64 mob{};
436  real64 dMob_dP{};
437  real64 dMob_dC[numComp]{};
438 
439  upwindMobilityGravity< numComp, numFluxSupportPoints, UPWIND >( numPhase,
440  ip,
441  seri,
442  sesri,
443  sei,
444  transmissibility,
445  dTrans_dPres,
446  totFlux,
447  pres,
448  gravCoef,
449  dCompFrac_dCompDens,
450  phaseMassDens,
451  dPhaseMassDens,
452  phaseMob,
453  dPhaseMob,
454  phaseVolFrac,
455  dPhaseVolFrac,
456  phaseCapPressure,
457  dPhaseCapPressure_dPhaseVolFrac,
458  capPressureFlag,
459  k_up,
460  mob,
461  dMob_dP,
462  dMob_dC );
463 
464  k_up_main = k_up;
465  mainMob = mob;
466  dMMob_dP = dMob_dP;
467  for( localIndex ic = 0; ic < numComp; ++ic )
468  {
469  dMMob_dC[ic] = dMob_dC[ic];
470  }
471 
472  //guard against no flow region
473  if( std::fabs( mainMob ) > 1e-20 )
474  {
475  fractionalFlow = mainMob / LvArray::math::max( totMob, minTotMob );
476  dFractionalFlow_dP[k_up_main] = dMMob_dP / LvArray::math::max( totMob, minTotMob );
477  for( localIndex jc = 0; jc < numComp; ++jc )
478  {
479  dFractionalFlow_dC[k_up_main][jc] = dMMob_dC[jc] / totMob;
480 
481  }
482 
483  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
484  {
485  dFractionalFlow_dP[ke] -= fractionalFlow * dTotMob_dP[k_up_ppu] / LvArray::math::max( totMob, minTotMob );
486 
487  for( localIndex jc = 0; jc < numComp; ++jc )
488  {
489  dFractionalFlow_dC[ke][jc] -= fractionalFlow * dTotMob_dC[k_up_ppu][jc] / LvArray::math::max( totMob, minTotMob );
490  }
491  }
492  }
493 }
494 
495 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
496 GEOS_HOST_DEVICE
497 static void
498 computeFractionalFlowCapillary( localIndex const numPhase,
499  localIndex const ip,
500  localIndex const (&seri)[numFluxSupportPoints],
501  localIndex const (&sesri)[numFluxSupportPoints],
502  localIndex const (&sei)[numFluxSupportPoints],
503  real64 const (&transmissibility)[2],
504  real64 const (&dTrans_dPres)[2],
505  localIndex const & k_up_ppu,
506  real64 const totFlux,
507  real64 const totMob,
508  real64 const (&dTotMob_dP)[numFluxSupportPoints],
509  real64 const (&dTotMob_dC)[numFluxSupportPoints][numComp],
510  ElementViewConst< arrayView1d< real64 const > > const & pres,
511  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
512  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
513  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
514  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
515  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
516  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
517  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
518  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
519  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
520  integer const capPressureFlag,
521  localIndex & k_up_main,
522  real64 & fractionalFlow,
523  real64 ( & dFractionalFlow_dP)[numFluxSupportPoints],
524  real64 ( & dFractionalFlow_dC)[numFluxSupportPoints][numComp]
525  )
526 {
527  // get var to memorized the numerator mobility properly upwinded
528  real64 mainMob{};
529  real64 dMMob_dP{};
530  real64 dMMob_dC[numComp]{};
531 
532  //reinit
533  //fractional flow too low to let the upstream phase flow
534  k_up_main = -1; //to throw error if unmodified
535  fractionalFlow = 0;
536  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
537  {
538  dFractionalFlow_dP[ke] = 0;
539  for( localIndex jc = 0; jc < numComp; ++jc )
540  {
541  dFractionalFlow_dC[ke][jc] = 0;
542  }
543  }
544 
545  localIndex k_up;
546  real64 mob{};
547  real64 dMob_dP{};
548  real64 dMob_dC[numComp]{};
549 
550  upwindMobilityCapillary< numComp, numFluxSupportPoints, UPWIND >( numPhase,
551  ip,
552  seri,
553  sesri,
554  sei,
555  transmissibility,
556  dTrans_dPres,
557  totFlux,
558  pres,
559  gravCoef,
560  dCompFrac_dCompDens,
561  phaseMassDens,
562  dPhaseMassDens,
563  phaseMob,
564  dPhaseMob,
565  dPhaseVolFrac,
566  phaseCapPressure,
567  dPhaseCapPressure_dPhaseVolFrac,
568  capPressureFlag,
569  k_up,
570  mob,
571  dMob_dP,
572  dMob_dC );
573 
574  k_up_main = k_up;
575  mainMob = mob;
576  dMMob_dP = dMob_dP;
577  for( localIndex ic = 0; ic < numComp; ++ic )
578  {
579  dMMob_dC[ic] = dMob_dC[ic];
580  }
581 
582  //guard against no flow region
583  if( std::fabs( mainMob ) > 1e-20 )
584  {
585  fractionalFlow = mainMob / LvArray::math::max( totMob, minTotMob );
586  dFractionalFlow_dP[k_up_main] = dMMob_dP / LvArray::math::max( totMob, minTotMob );
587  for( localIndex jc = 0; jc < numComp; ++jc )
588  {
589  dFractionalFlow_dC[k_up_main][jc] = dMMob_dC[jc] / LvArray::math::max( totMob, minTotMob );
590 
591  }
592 
593  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
594  {
595  dFractionalFlow_dP[ke] -= fractionalFlow * dTotMob_dP[k_up_ppu] / LvArray::math::max( totMob, minTotMob );
596 
597  for( localIndex jc = 0; jc < numComp; ++jc )
598  {
599  dFractionalFlow_dC[ke][jc] -= fractionalFlow * dTotMob_dC[k_up_ppu][jc] / LvArray::math::max( totMob, minTotMob );
600  }
601  }
602  }
603 }
604 
612 struct computePotentialViscous
613 {
614 
615  template< localIndex numComp, localIndex numFluxSupportPoints >
616  GEOS_HOST_DEVICE
617  static void compute( localIndex const GEOS_UNUSED_PARAM( numPhase ),
618  localIndex const GEOS_UNUSED_PARAM( ip ),
619  localIndex const (&GEOS_UNUSED_PARAM( seri ))[numFluxSupportPoints],
620  localIndex const (&GEOS_UNUSED_PARAM( sesri ))[numFluxSupportPoints],
621  localIndex const (&GEOS_UNUSED_PARAM( sei ))[numFluxSupportPoints],
622  real64 const (&GEOS_UNUSED_PARAM( transmissibility ))[2],
623  real64 const (&GEOS_UNUSED_PARAM( dTrans_dPres ))[2],
624  real64 const totFlux,
625  ElementViewConst< arrayView1d< real64 const > > const & GEOS_UNUSED_PARAM( gravCoef ),
626  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const &
627  GEOS_UNUSED_PARAM( dCompFrac_dCompDens ),
628  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const &
629  GEOS_UNUSED_PARAM( phaseMassDens ),
630  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const &
631  GEOS_UNUSED_PARAM( dPhaseMassDens ),
632  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &
633  GEOS_UNUSED_PARAM( dPhaseVolFrac ),
634  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const &
635  GEOS_UNUSED_PARAM( phaseCapPressure ),
636  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const &
637  GEOS_UNUSED_PARAM( dPhaseCapPressure_dPhaseVolFrac ),
638  real64 & pot,
639  real64( &GEOS_UNUSED_PARAM( dPot_dPres ))[numFluxSupportPoints],
640  real64( &GEOS_UNUSED_PARAM( dPot_dComp ))[numFluxSupportPoints][numComp],
641  real64( &GEOS_UNUSED_PARAM( dProp_dComp ))[numComp] )
642  {
643  pot = totFlux;
644  //could be relevant for symmetry to include derivative
645 
646  }
647 };
648 
651 struct computePotentialGravity
652 {
657  template< localIndex numComp, localIndex numFluxSupportPoints >
658  GEOS_HOST_DEVICE
659  static void compute( localIndex const GEOS_UNUSED_PARAM( numPhase ),
660  localIndex const ip,
661  localIndex const (&seri)[numFluxSupportPoints],
662  localIndex const (&sesri)[numFluxSupportPoints],
663  localIndex const (&sei)[numFluxSupportPoints],
664  real64 const (&transmissibility)[2],
665  real64 const (&dTrans_dPres)[2],
666  real64 const GEOS_UNUSED_PARAM( totFlux ),
667  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
668  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
669  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
670  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
671  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
672  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &
673  GEOS_UNUSED_PARAM( dPhaseVolFrac ),
674  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const &
675  GEOS_UNUSED_PARAM( phaseCapPressure ),
676  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const &
677  GEOS_UNUSED_PARAM( dPhaseCapPressure_dPhaseVolFrac ),
678  real64 & pot,
679  real64 ( & dPot_dPres )[numFluxSupportPoints],
680  real64 (& dPot_dComp )[numFluxSupportPoints][numComp],
681  real64 ( & dProp_dComp )[numComp] )
682  {
683  //working arrays
684  real64 densMean{};
685  real64 dDensMean_dPres[numFluxSupportPoints]{};
686  real64 dDensMean_dComp[numFluxSupportPoints][numComp]{};
687 
688  //init
689  pot = 0.0;
690  for( localIndex i = 0; i < numFluxSupportPoints; ++i )
691  {
692  dPot_dPres[i] = 0.0;
693  for( localIndex jc = 0; jc < numComp; ++jc )
694  {
695  dPot_dComp[i][jc] = 0.0;
696  dProp_dComp[jc] = 0.0;
697  }
698  }
699 
700  //inner loop to get average density
701  integer denom = 0;
702  for( localIndex i = 0; i < numFluxSupportPoints; ++i )
703  {
704  localIndex const er = seri[i];
705  localIndex const esr = sesri[i];
706  localIndex const ei = sei[i];
707 
708  bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0);
709  if( !phaseExists )
710  {
711  continue;
712  }
713 
714  // density
715  real64 const density = phaseMassDens[er][esr][ei][0][ip];
716  real64 const dDens_dPres = dPhaseMassDens[er][esr][ei][0][ip][Deriv::dP];
717 
718  applyChainRule( numComp,
719  dCompFrac_dCompDens[er][esr][ei],
720  dPhaseMassDens[er][esr][ei][0][ip],
721  dProp_dComp,
722  Deriv::dC );
723 
724  // average density and derivatives
725  densMean += density;
726  dDensMean_dPres[i] = dDens_dPres;
727  for( localIndex jc = 0; jc < numComp; ++jc )
728  {
729  dDensMean_dComp[i][jc] = dProp_dComp[jc];
730  }
731  denom++;
732  }
733  if( denom > 1 )
734  {
735  densMean /= denom;
736  for( localIndex i = 0; i < numFluxSupportPoints; ++i )
737  {
738  dDensMean_dPres[i] /= denom;
739  for( integer jc = 0; jc < numComp; ++jc )
740  {
741  dDensMean_dComp[i][jc] /= denom;
742  }
743  }
744  }
745 
746  // compute potential difference MPFA-style
747  for( localIndex i = 0; i < numFluxSupportPoints; ++i )
748  {
749  localIndex const er = seri[i];
750  localIndex const esr = sesri[i];
751  localIndex const ei = sei[i];
752 
753  real64 const gravD = transmissibility[i] * gravCoef[er][esr][ei];
754  real64 const dGravD_dP = dTrans_dPres[i] * gravCoef[er][esr][ei];
755  pot += densMean * gravD;
756 
757  // need to add contributions from both cells the mean density depends on
758  for( localIndex j = 0; j < numFluxSupportPoints; ++j )
759  {
760  dPot_dPres[j] += dDensMean_dPres[j] * gravD + densMean * dGravD_dP;
761  for( localIndex jc = 0; jc < numComp; ++jc )
762  {
763  dPot_dComp[j][jc] += dDensMean_dComp[j][jc] * gravD;
764  }
765  }
766  }
767 
768  }
769 };
770 
773 struct computePotentialCapillary
774 {
779  template< localIndex numComp, localIndex numFluxSupportPoints >
780  GEOS_HOST_DEVICE
781  static void compute( localIndex const numPhase,
782  localIndex const ip,
783  localIndex const (&seri)[numFluxSupportPoints],
784  localIndex const (&sesri)[numFluxSupportPoints],
785  localIndex const (&sei)[numFluxSupportPoints],
786  real64 const (&transmissibility)[2],
787  real64 const (&dTrans_dPres)[2],
788  real64 const GEOS_UNUSED_PARAM( totFlux ),
789  ElementViewConst< arrayView1d< real64 const > > const & GEOS_UNUSED_PARAM( gravCoef ),
790  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const &
791  GEOS_UNUSED_PARAM( dCompFrac_dCompDens ),
792  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const &
793  GEOS_UNUSED_PARAM( phaseMassDens ),
794  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const &
795  GEOS_UNUSED_PARAM( dPhaseMassDens ),
796  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
797  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
798  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
799  real64 & pot,
800  real64 ( & dPot_dPres)[numFluxSupportPoints],
801  real64 (& dPot_dComp)[numFluxSupportPoints][numComp],
802  real64( &GEOS_UNUSED_PARAM( dProp_dComp ))[numComp] )
803  {
804 
805 
806  for( localIndex i = 0; i < numFluxSupportPoints; ++i )
807  {
808  localIndex const er = seri[i];
809  localIndex const esr = sesri[i];
810  localIndex const ei = sei[i];
811 
812  pot += transmissibility[i] * phaseCapPressure[er][esr][ei][0][ip];
813  // need to add contributions from both cells
814  for( localIndex jp = 0; jp < numPhase; ++jp )
815  {
816 
817  real64 const dCapPressure_dS = dPhaseCapPressure_dPhaseVolFrac[er][esr][ei][0][ip][jp];
818  dPot_dPres[i] +=
819  transmissibility[i] * dCapPressure_dS * dPhaseVolFrac[er][esr][ei][jp][Deriv::dP]
820  + dTrans_dPres[i] * phaseCapPressure[er][esr][ei][0][jp];
821 
822  for( localIndex jc = 0; jc < numComp; ++jc )
823  {
824  dPot_dComp[i][jc] += transmissibility[i] * dCapPressure_dS *
825  dPhaseVolFrac[er][esr][ei][jp][Deriv::dC + jc];
826  }
827  }
828  }
829  }
830 };
831 
833 
834 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
835 GEOS_HOST_DEVICE
836 static void computePotentialFluxesGravity( localIndex const numPhase,
837  localIndex const ip,
838  localIndex const (&seri)[numFluxSupportPoints],
839  localIndex const (&sesri)[numFluxSupportPoints],
840  localIndex const (&sei)[numFluxSupportPoints],
841  real64 const (&transmissibility)[2],
842  real64 const (&dTrans_dPres)[2],
843  localIndex const & k_up_ppu,
844  real64 const totFlux,
845  real64 const totMob,
846  real64 const (&dTotMob_dP)[numFluxSupportPoints],
847  real64 const (&dTotMob_dC)[numFluxSupportPoints][numComp],
848  ElementViewConst< arrayView1d< real64 const > > const & pres,
849  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
850  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
851  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
852  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
853  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
854  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
855  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
856  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
857  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
858  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
859  localIndex const capPressureFlag,
860  localIndex( &k_up),
861  localIndex (&k_up_o),
862  real64 & phaseFlux,
863  real64 (& dPhaseFlux_dP)[numFluxSupportPoints],
864  real64 ( & dPhaseFlux_dC)[numFluxSupportPoints][numComp] )
865 {
866 
867  real64 fflow{};
868  real64 dFflow_dP[numFluxSupportPoints]{};
869  real64 dFflow_dC[numFluxSupportPoints][numComp]{};
870 
871  real64 pot{};
872  real64 dPot_dP[numFluxSupportPoints]{};
873  real64 dPot_dC[numFluxSupportPoints][numComp]{};
874  real64 dProp_dC[numComp]{};
875 
876  //
877  UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >( numPhase,
878  ip,
879  seri,
880  sesri,
881  sei,
882  transmissibility,
883  dTrans_dPres,
884  totFlux,
885  gravCoef,
886  dCompFrac_dCompDens,
887  phaseMassDens,
888  dPhaseMassDens,
889  phaseVolFrac,
890  dPhaseVolFrac,
891  phaseCapPressure,
892  dPhaseCapPressure_dPhaseVolFrac,
893  pot,
894  dPot_dP,
895  dPot_dC,
896  dProp_dC );
897 
898  // and the fractional flow for gravitational part as \lambda_i^{up}/\sum_{numPhase}(\lambda_k^{up}) with up decided upon
899  // the Upwind strategy
900  UpwindHelpers::computeFractionalFlowGravity< numComp, numFluxSupportPoints, UPWIND >( numPhase,
901  ip,
902  seri,
903  sesri,
904  sei,
905  transmissibility,
906  dTrans_dPres,
907  k_up_ppu,
908  totFlux,
909  totMob,
910  dTotMob_dP,
911  dTotMob_dC,
912  pres,
913  gravCoef,
914  dCompFrac_dCompDens,
915  phaseMassDens,
916  dPhaseMassDens,
917  phaseMob,
918  dPhaseMob,
919  phaseVolFrac,
920  dPhaseVolFrac,
921  phaseCapPressure,
922  dPhaseCapPressure_dPhaseVolFrac,
923  capPressureFlag,
924  k_up,
925  fflow,
926  dFflow_dP,
927  dFflow_dC );
928 
929 
930  for( localIndex jp = 0; jp < numPhase; ++jp )
931  {
932  if( ip != jp )
933  {
934 
935  real64 potOther{};
936  real64 dPotOther_dP[numFluxSupportPoints]{};
937  real64 dPotOther_dC[numFluxSupportPoints][numComp]{};
938  real64 dPropOther_dC[numComp]{};
939 
940  //Fetch pot for phase j!=i defined as \rho_j g dz/dx
941  UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >( numPhase,
942  jp,
943  seri,
944  sesri,
945  sei,
946  transmissibility,
947  dTrans_dPres,
948  totFlux,
949  gravCoef,
950  dCompFrac_dCompDens,
951  phaseMassDens,
952  dPhaseMassDens,
953  phaseVolFrac,
954  dPhaseVolFrac,
955  phaseCapPressure,
956  dPhaseCapPressure_dPhaseVolFrac,
957  potOther,
958  dPotOther_dP,
959  dPotOther_dC,
960  dPropOther_dC );
961 
962  //Eventually get the mobility of the second phase
963  real64 mobOther{};
964  real64 dMobOther_dP{};
965  real64 dMobOther_dC[numComp]{};
966 
967  // and the other mobility for gravitational part as \lambda_j^{up} with up decided upon
968  // the Upwind strategy - Note that it should be the same as the gravitational fractional flow
969 
970  UpwindHelpers::upwindMobilityGravity< numComp, numFluxSupportPoints, UPWIND >( numPhase,
971  jp,
972  seri,
973  sesri,
974  sei,
975  transmissibility,
976  dTrans_dPres,
977  totFlux,
978  pres,
979  gravCoef,
980  dCompFrac_dCompDens,
981  phaseMassDens,
982  dPhaseMassDens,
983  phaseMob,
984  dPhaseMob,
985  phaseVolFrac,
986  dPhaseVolFrac,
987  phaseCapPressure,
988  dPhaseCapPressure_dPhaseVolFrac,
989  capPressureFlag,
990  k_up_o,
991  mobOther,
992  dMobOther_dP,
993  dMobOther_dC );
994 
995 
996  // Assembling gravitational flux phase-wise as \phi_{i,g} = \sum_{k\nei} \lambda_k^{up,g} f_k^{up,g} (G_i - G_k)
997  phaseFlux -= fflow * mobOther * (pot - potOther);
998  dPhaseFlux_dP[k_up_o] -= fflow * dMobOther_dP * (pot - potOther);
999  for( localIndex jc = 0; jc < numComp; ++jc )
1000  {
1001  dPhaseFlux_dC[k_up_o][jc] -= fflow * dMobOther_dC[jc] * (pot - potOther);
1002  }
1003 
1004  //mob related part of dFflow_dP is only upstream defined but totMob related is defined everywhere
1005  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
1006  {
1007  dPhaseFlux_dP[ke] -= dFflow_dP[ke] * mobOther * (pot - potOther);
1008 
1009  for( localIndex jc = 0; jc < numComp; ++jc )
1010  {
1011  dPhaseFlux_dC[ke][jc] -= dFflow_dC[ke][jc] * mobOther * (pot - potOther);
1012  }
1013  }
1014 
1015  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
1016  {
1017  dPhaseFlux_dP[ke] -= fflow * mobOther * (dPot_dP[ke] - dPotOther_dP[ke]);
1018  for( localIndex jc = 0; jc < numComp; ++jc )
1019  {
1020  dPhaseFlux_dC[ke][jc] -= fflow * mobOther * (dPot_dC[ke][jc] - dPotOther_dC[ke][jc]);
1021  }
1022  }
1023  }
1024  }
1025 
1026 }
1027 
1028 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
1029 GEOS_HOST_DEVICE
1030 static void computePotentialFluxesCapillary( localIndex const numPhase,
1031  localIndex const ip,
1032  localIndex const (&seri)[numFluxSupportPoints],
1033  localIndex const (&sesri)[numFluxSupportPoints],
1034  localIndex const (&sei)[numFluxSupportPoints],
1035  real64 const (&transmissibility)[2],
1036  real64 const (&dTrans_dPres)[2],
1037  localIndex const & k_up_ppu,
1038  real64 const totFlux,
1039  real64 const totMob,
1040  real64 const (&dTotMob_dP)[numFluxSupportPoints],
1041  real64 const (&dTotMob_dC)[numFluxSupportPoints][numComp],
1042  ElementViewConst< arrayView1d< real64 const > > const & pres,
1043  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1044  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1045  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
1046  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1047  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1048  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1049  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1050  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1051  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1052  localIndex const capPressureFlag,
1053  localIndex( &k_up),
1054  localIndex (&k_up_o),
1055  real64 & phaseFlux,
1056  real64 (& dPhaseFlux_dP)[numFluxSupportPoints],
1057  real64 ( & dPhaseFlux_dC)[numFluxSupportPoints][numComp] )
1058 {
1059 
1060  real64 fflow{};
1061  real64 dFflow_dP[numFluxSupportPoints]{};
1062  real64 dFflow_dC[numFluxSupportPoints][numComp]{};
1063 
1064  real64 pot{};
1065  real64 dPot_dP[numFluxSupportPoints]{};
1066  real64 dPot_dC[numFluxSupportPoints][numComp]{};
1067  real64 dProp_dC[numComp]{};
1068 
1069  UpwindHelpers::computePotentialCapillary::compute< numComp, numFluxSupportPoints >( numPhase,
1070  ip,
1071  seri,
1072  sesri,
1073  sei,
1074  transmissibility,
1075  dTrans_dPres,
1076  totFlux,
1077  gravCoef,
1078  dCompFrac_dCompDens,
1079  phaseMassDens,
1080  dPhaseMassDens,
1081  dPhaseVolFrac,
1082  phaseCapPressure,
1083  dPhaseCapPressure_dPhaseVolFrac,
1084  pot,
1085  dPot_dP,
1086  dPot_dC,
1087  dProp_dC );
1088 
1089  // and the fractional flow for gravitational part as \lambda_i^{up}/\sum_{numPhase}(\lambda_k^{up}) with up decided upon
1090  // the Upwind strategy
1091  UpwindHelpers::computeFractionalFlowCapillary< numComp, numFluxSupportPoints, UPWIND >( numPhase,
1092  ip,
1093  seri,
1094  sesri,
1095  sei,
1096  transmissibility,
1097  dTrans_dPres,
1098  k_up_ppu,
1099  totFlux,
1100  totMob,
1101  dTotMob_dP,
1102  dTotMob_dC,
1103  pres,
1104  gravCoef,
1105  dCompFrac_dCompDens,
1106  phaseMassDens,
1107  dPhaseMassDens,
1108  phaseMob,
1109  dPhaseMob,
1110  dPhaseVolFrac,
1111  phaseCapPressure,
1112  dPhaseCapPressure_dPhaseVolFrac,
1113  capPressureFlag,
1114  k_up,
1115  fflow,
1116  dFflow_dP,
1117  dFflow_dC );
1118 
1119 
1120  for( localIndex jp = 0; jp < numPhase; ++jp )
1121  {
1122  if( ip != jp )
1123  {
1124 
1125  real64 potOther{};
1126  real64 dPotOther_dP[numFluxSupportPoints]{};
1127  real64 dPotOther_dC[numFluxSupportPoints][numComp]{};
1128  real64 dPropOther_dC[numComp]{};
1129 
1130  //Fetch pot for phase j!=i defined as \rho_j g dz/dx
1131  UpwindHelpers::computePotentialCapillary::compute< numComp, numFluxSupportPoints >( numPhase,
1132  jp,
1133  seri,
1134  sesri,
1135  sei,
1136  transmissibility,
1137  dTrans_dPres,
1138  totFlux,
1139  gravCoef,
1140  dCompFrac_dCompDens,
1141  phaseMassDens,
1142  dPhaseMassDens,
1143  dPhaseVolFrac,
1144  phaseCapPressure,
1145  dPhaseCapPressure_dPhaseVolFrac,
1146  potOther,
1147  dPotOther_dP,
1148  dPotOther_dC,
1149  dPropOther_dC );
1150 
1151  //Eventually get the mobility of the second phase
1152  real64 mobOther{};
1153  real64 dMobOther_dP{};
1154  real64 dMobOther_dC[numComp]{};
1155 
1156  // and the other mobility for gravitational part as \lambda_j^{up} with up decided upon
1157  // the Upwind strategy - Note that it should be the same as the gravitational fractional flow
1158 
1159  UpwindHelpers::upwindMobilityCapillary< numComp, numFluxSupportPoints, UPWIND >( numPhase,
1160  jp,
1161  seri,
1162  sesri,
1163  sei,
1164  transmissibility,
1165  dTrans_dPres,
1166  totFlux,
1167  pres,
1168  gravCoef,
1169  dCompFrac_dCompDens,
1170  phaseMassDens,
1171  dPhaseMassDens,
1172  phaseMob,
1173  dPhaseMob,
1174  dPhaseVolFrac,
1175  phaseCapPressure,
1176  dPhaseCapPressure_dPhaseVolFrac,
1177  capPressureFlag,
1178  k_up_o,
1179  mobOther,
1180  dMobOther_dP,
1181  dMobOther_dC );
1182 
1183 
1184  // Assembling gravitational flux phase-wise as \phi_{i,g} = \sum_{k\nei} \lambda_k^{up,g} f_k^{up,g} (G_i - G_k)
1185  phaseFlux -= fflow * mobOther * (pot - potOther);
1186  dPhaseFlux_dP[k_up_o] -= fflow * dMobOther_dP * (pot - potOther);
1187  for( localIndex jc = 0; jc < numComp; ++jc )
1188  {
1189  dPhaseFlux_dC[k_up_o][jc] -= fflow * dMobOther_dC[jc] * (pot - potOther);
1190  }
1191 
1192  //mob related part of dFflow_dP is only upstream defined but totMob related is defined everywhere
1193  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
1194  {
1195  dPhaseFlux_dP[ke] -= dFflow_dP[ke] * mobOther * (pot - potOther);
1196 
1197  for( localIndex jc = 0; jc < numComp; ++jc )
1198  {
1199  dPhaseFlux_dC[ke][jc] -= dFflow_dC[ke][jc] * mobOther * (pot - potOther);
1200  }
1201  }
1202 
1203  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
1204  {
1205  dPhaseFlux_dP[ke] -= fflow * mobOther * (dPot_dP[ke] - dPotOther_dP[ke]);
1206  for( localIndex jc = 0; jc < numComp; ++jc )
1207  {
1208  dPhaseFlux_dC[ke][jc] -= fflow * mobOther * (dPot_dC[ke][jc] - dPotOther_dC[ke][jc]);
1209  }
1210  }
1211  }
1212  }
1213 }
1214 
1215 }//end of struct UpwindHelpers
1216 
1217 /************************* UPWIND ******************/
1218 
1223 class UpwindScheme
1224 {
1225 
1226 public:
1227 
1228  //default ctor
1229  UpwindScheme() = default;
1230 
1231  //usual copy ctor
1232  UpwindScheme( UpwindScheme const & scheme ) = default;
1233 
1234  //default move ctor
1235  UpwindScheme( UpwindScheme && ) = default;
1236 
1237  //deleted copy and move assignement
1238  UpwindScheme & operator=( UpwindScheme const & ) = delete;
1239 
1240  UpwindScheme & operator=( UpwindScheme && ) = delete;
1241 
1242  virtual ~UpwindScheme() = default;
1243 
1244  template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
1245  GEOS_HOST_DEVICE
1246  inline
1247  void getUpwindDirectionViscous( localIndex const numPhase,
1248  localIndex const ip,
1249  localIndex const (&seri)[numFluxSupportPoints],
1250  localIndex const (&sesri)[numFluxSupportPoints],
1251  localIndex const (&sei)[numFluxSupportPoints],
1252  real64 const (&transmissibility)[2],
1253  real64 const (&dTrans_dPres)[2],
1254  real64 const totFlux, //in fine should be a ElemnetViewConst once seq form are in place
1255  ElementViewConst< arrayView1d< real64 const > > const & pres,
1256  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1257  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1258  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1259  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1260  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1261  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1262  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1263  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1264  integer const capPressureFlag,
1265  localIndex & upwindDir
1266  )
1267  {
1268  real64 pot{};
1269 
1272  UPWIND::template computePotentialViscous< numComp, numFluxSupportPoints >( numPhase,
1273  ip,
1274  seri,
1275  sesri,
1276  sei,
1277  transmissibility,
1278  dTrans_dPres,
1279  totFlux,
1280  pres,
1281  gravCoef,
1282  phaseMob,
1283  dCompFrac_dCompDens,
1284  phaseMassDens,
1285  dPhaseMassDens,
1286  dPhaseVolFrac,
1287  phaseCapPressure,
1288  dPhaseCapPressure_dPhaseVolFrac,
1289  capPressureFlag,
1290  pot );
1291 
1292  //all definition has been changed to fit pot>0 => first cell is upstream
1293  upwindDir = (pot > 0) ? 0 : 1;
1294  }
1295 
1296 
1297  template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
1298  GEOS_HOST_DEVICE
1299  void getUpwindDirectionGravity( localIndex const numPhase,
1300  localIndex const ip,
1301  localIndex const (&seri)[numFluxSupportPoints],
1302  localIndex const (&sesri)[numFluxSupportPoints],
1303  localIndex const (&sei)[numFluxSupportPoints],
1304  real64 const (&transmissibility)[2],
1305  real64 const (&dTrans_dPres)[2],
1306  real64 const totFlux,
1307  ElementViewConst< arrayView1d< real64 const > > const & pres,
1308  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1309  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1310  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1311  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1312  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1313  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
1314  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1315  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1316  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1317  integer const capPressureFlag,
1318  localIndex & upwindDir
1319  )
1320  {
1321  real64 pot{};
1322 
1325  UPWIND::template computePotentialGravity< numComp, numFluxSupportPoints >( numPhase,
1326  ip,
1327  seri,
1328  sesri,
1329  sei,
1330  transmissibility,
1331  dTrans_dPres,
1332  totFlux,
1333  pres,
1334  gravCoef,
1335  phaseMob,
1336  dCompFrac_dCompDens,
1337  phaseMassDens,
1338  dPhaseMassDens,
1339  phaseVolFrac,
1340  dPhaseVolFrac,
1341  phaseCapPressure,
1342  dPhaseCapPressure_dPhaseVolFrac,
1343  capPressureFlag,
1344  pot );
1345 
1346  //all definition has been changed to fit pot>0 => first cell is upstream
1347  upwindDir = (pot >= 0) ? 0 : 1;
1348  }
1349 
1350 
1351  template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
1352  GEOS_HOST_DEVICE
1353  void getUpwindDirectionCapillary( localIndex const numPhase,
1354  localIndex const ip,
1355  localIndex const (&seri)[numFluxSupportPoints],
1356  localIndex const (&sesri)[numFluxSupportPoints],
1357  localIndex const (&sei)[numFluxSupportPoints],
1358  real64 const (&transmissibility)[2],
1359  real64 const (&dTrans_dPres)[2],
1360  real64 const totFlux, //in fine should be a ElemnetViewConst once seq form are in
1361  // place
1362  ElementViewConst< arrayView1d< real64 const > > const & pres,
1363  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1364  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1365  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1366  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1367  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1368  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1369  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1370  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1371  integer const capPressureFlag,
1372  localIndex & upwindDir
1373  )
1374  {
1375  real64 pot{};
1376 
1377  // each derived concrete class has to define a computePotential method that is calling UpwindScheme::potential method with a specific
1378  // lamda defining how to get these potentials
1379  UPWIND::template computePotentialCapillary< numComp, numFluxSupportPoints >( numPhase,
1380  ip,
1381  seri,
1382  sesri,
1383  sei,
1384  transmissibility,
1385  dTrans_dPres,
1386  totFlux,
1387  pres,
1388  gravCoef,
1389  phaseMob,
1390  dCompFrac_dCompDens,
1391  phaseMassDens,
1392  dPhaseMassDens,
1393  dPhaseVolFrac,
1394  phaseCapPressure,
1395  dPhaseCapPressure_dPhaseVolFrac,
1396  capPressureFlag,
1397  pot );
1398 
1399  //all definition has been changed to fit pot>0 => first cell is upstream
1400  upwindDir = (pot >= 0) ? 0 : 1;
1401  }
1402 
1403 
1404 
1405  // templated way of evaluating the potential (to the exception of viscous one) which relies on
1406  // up-or-downwinded mobility terms pre-multiplying potential differences
1407  template< localIndex numComp, localIndex numFluxSupportPoints, typename LAMBDA >
1408  GEOS_HOST_DEVICE
1409  static void potential( localIndex numPhase,
1410  localIndex ip,
1411  localIndex const (&seri)[numFluxSupportPoints],
1412  localIndex const (&sesri)[numFluxSupportPoints],
1413  localIndex const (&sei)[numFluxSupportPoints],
1414  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1415  real64 & weightedPotential,
1416  LAMBDA && fn )
1417  {
1418  //getPhase Pot
1419  real64 pot{};
1420  real64 pot_dP[numFluxSupportPoints]{};
1421  real64 pot_dC[numFluxSupportPoints][numComp]{};
1422  real64 dProp_dC[numComp]{};
1423 
1424  fn( ip, pot, pot_dP, pot_dC, dProp_dC );
1425 
1426  localIndex const k_up = 0;
1427  localIndex const k_dw = 1;
1428 
1429  //loop other other phases to form
1430  for( localIndex jp = 0; jp < numPhase; ++jp )
1431  {
1432  if( jp != ip )
1433  {
1434  localIndex const er_up = seri[k_up];
1435  localIndex const esr_up = sesri[k_up];
1436  localIndex const ei_up = sei[k_up];
1437 
1438  localIndex const er_dw = seri[k_dw];
1439  localIndex const esr_dw = sesri[k_dw];
1440  localIndex const ei_dw = sei[k_dw];
1441 
1442  real64 potOther{};
1443  real64 potOther_dP[numFluxSupportPoints]{};
1444  real64 potOther_dC[numFluxSupportPoints][numComp]{};
1445  real64 dPropOther_dC[numComp]{};
1446 
1447  fn( jp, potOther, potOther_dP, potOther_dC, dPropOther_dC );
1448 
1449  real64 const mob_up = phaseMob[er_up][esr_up][ei_up][jp];
1450  real64 const mob_dw = phaseMob[er_dw][esr_dw][ei_dw][jp];
1451 
1452  weightedPotential += (pot - potOther >= 0) ? mob_dw * (potOther - pot) : mob_up * (potOther - pot);
1453 
1454  }
1455  }
1456  }
1457 
1458 };
1459 
1464 class HybridUpwind : public UpwindScheme
1465 {
1466 
1467 public:
1468  template< localIndex numComp, localIndex numFluxSupportPoints >
1469  GEOS_HOST_DEVICE
1470  static
1471  void computePotentialViscous( localIndex const numPhase,
1472  localIndex const ip,
1473  localIndex const (&seri)[numFluxSupportPoints],
1474  localIndex const (&sesri)[numFluxSupportPoints],
1475  localIndex const (&sei)[numFluxSupportPoints],
1476  real64 const (&transmissibility)[2],
1477  real64 const (&dTrans_dPres)[2],
1478  real64 const totalFlux,
1479  ElementViewConst< arrayView1d< real64 const > > const & GEOS_UNUSED_PARAM( pres ),
1480  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1481  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &
1482  GEOS_UNUSED_PARAM( phaseMob ),
1483  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1484  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1485  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1486  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1487  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1488  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1489  integer const GEOS_UNUSED_PARAM( capPressureFlag ),
1490  real64 & potential
1491  )
1492  {
1493  real64 dPot_dP[numFluxSupportPoints]{};
1494  real64 dPot_dC[numFluxSupportPoints][numComp]{};
1495  real64 dProp_dC[numComp]{};
1496 
1497 
1498  UpwindHelpers::computePotentialViscous::compute< numComp, numFluxSupportPoints >(
1499  numPhase,
1500  ip,
1501  seri,
1502  sesri,
1503  sei,
1504  transmissibility,
1505  dTrans_dPres,
1506  totalFlux,
1507  gravCoef,
1508  dCompFrac_dCompDens,
1509  phaseMassDens,
1510  dPhaseMassDens,
1511  dPhaseVolFrac,
1512  phaseCapPressure,
1513  dPhaseCapPressure_dPhaseVolFrac,
1514  potential,
1515  dPot_dP,
1516  dPot_dC,
1517  dProp_dC );
1518  }
1519 
1520  template< localIndex numComp, localIndex numFluxSupportPoints >
1521  GEOS_HOST_DEVICE
1522  static
1523  void computePotentialGravity( localIndex const numPhase,
1524  localIndex const ip,
1525  localIndex const (&seri)[numFluxSupportPoints],
1526  localIndex const (&sesri)[numFluxSupportPoints],
1527  localIndex const (&sei)[numFluxSupportPoints],
1528  real64 const (&transmissibility)[2],
1529  real64 const (&dTrans_dPres)[2],
1530  real64 const totalFlux,
1531  ElementViewConst< arrayView1d< real64 const > > const & GEOS_UNUSED_PARAM( pres ),
1532  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1533  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1534  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1535  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1536  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1537  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
1538  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1539  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1540  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1541  integer const GEOS_UNUSED_PARAM( capPressureFlag ),
1542  real64 & potential
1543  )
1544  {
1545 
1546  //Form total velocity
1547  potential = 0;
1548 
1549  //the arg lambda allows us to access some genericity
1550  UpwindScheme::template potential< numComp, numFluxSupportPoints >( numPhase, ip, seri, sesri, sei,
1551  phaseMob, potential,
1552  [&]( localIndex ipp,
1553  real64 & potential_,
1554  real64 (& dPotential_dP_)[numFluxSupportPoints],
1555  real64 (& dPotential_dC_)[numFluxSupportPoints][numComp],
1556  real64 (& dProp_dC)[numComp] ) {
1557 
1558  UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >(
1559  numPhase,
1560  ipp,
1561  seri,
1562  sesri,
1563  sei,
1564  transmissibility,
1565  dTrans_dPres,
1566  totalFlux,
1567  gravCoef,
1568  dCompFrac_dCompDens,
1569  phaseMassDens,
1570  dPhaseMassDens,
1571  phaseVolFrac,
1572  dPhaseVolFrac,
1573  phaseCapPressure,
1574  dPhaseCapPressure_dPhaseVolFrac,
1575  potential_,
1576  dPotential_dP_,
1577  dPotential_dC_,
1578  dProp_dC );
1579 
1580  } );
1581  }
1582 
1583 
1584  template< localIndex numComp, localIndex numFluxSupportPoints >
1585  GEOS_HOST_DEVICE
1586  static
1587  void computePotentialCapillary( localIndex const numPhase,
1588  localIndex const ip,
1589  localIndex const (&seri)[numFluxSupportPoints],
1590  localIndex const (&sesri)[numFluxSupportPoints],
1591  localIndex const (&sei)[numFluxSupportPoints],
1592  real64 const (&transmissibility)[2],
1593  real64 const (&dTrans_dPres)[2],
1594  real64 const totalFlux,
1595  ElementViewConst< arrayView1d< real64 const > > const & GEOS_UNUSED_PARAM( pres ),
1596  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1597  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &
1598  phaseMob,
1599  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1600  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1601  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1602  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1603  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1604  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1605  integer const GEOS_UNUSED_PARAM( capPressureFlag ),
1606  real64 & potential )
1607 
1608  {
1609 
1610  //Form total velocity
1611  potential = 0;
1612 
1613  //the arg lambda allows us to access some genericity
1614  UpwindScheme::template potential< numComp, numFluxSupportPoints >( numPhase, ip, seri, sesri, sei,
1615  phaseMob, potential,
1616  [&]( localIndex ipp,
1617  real64 & potential_,
1618  real64 (& dPotential_dP_)[numFluxSupportPoints],
1619  real64 (& dPotential_dC_)[numFluxSupportPoints][numComp],
1620  real64 (& dProp_dC)[numComp] ) {
1621 
1622  UpwindHelpers::computePotentialCapillary::compute< numComp, numFluxSupportPoints >(
1623  numPhase,
1624  ipp,
1625  seri,
1626  sesri,
1627  sei,
1628  transmissibility,
1629  dTrans_dPres,
1630  totalFlux,
1631  gravCoef,
1632  dCompFrac_dCompDens,
1633  phaseMassDens,
1634  dPhaseMassDens,
1635  dPhaseVolFrac,
1636  phaseCapPressure,
1637  dPhaseCapPressure_dPhaseVolFrac,
1638  potential_,
1639  dPotential_dP_,
1640  dPotential_dC_,
1641  dProp_dC );
1642  } );
1643  }
1644 
1645 };
1646 
1647 /*** IHU ***/
1648 
1649 struct IHUPhaseFlux
1650 {
1651 
1652  using UPWIND_SCHEME = HybridUpwind;
1653 
1682  template< integer numComp, integer numFluxSupportPoints >
1683  GEOS_HOST_DEVICE
1684  static void
1685  compute( integer const numPhase,
1686  integer const ip,
1687  integer const hasCapPressure,
1688  localIndex const ( &seri )[numFluxSupportPoints],
1689  localIndex const ( &sesri )[numFluxSupportPoints],
1690  localIndex const ( &sei )[numFluxSupportPoints],
1691  real64 const ( &trans )[2],
1692  real64 const ( &dTrans_dPres )[2],
1693  ElementViewConst< arrayView1d< real64 const > > const & pres,
1694  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1695  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1696  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
1697  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
1698  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1699  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const & phaseCompFrac,
1700  ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac,
1701  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1702  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1703  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1704  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1705  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1706  localIndex & k_up,
1707  real64 & potGrad,
1708  real64 ( &phaseFlux ),
1709  real64 ( & dPhaseFlux_dP )[numFluxSupportPoints],
1710  real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp],
1711  real64 ( & compFlux )[numComp],
1712  real64 ( & dCompFlux_dP )[numFluxSupportPoints][numComp],
1713  real64 ( & dCompFlux_dC )[numFluxSupportPoints][numComp][numComp] )
1714  {
1715 
1716  //loop over all phases to form total velocity
1717  real64 totFlux{};
1718  real64 dTotFlux_dP[numFluxSupportPoints]{};
1719  real64 dTotFlux_dC[numFluxSupportPoints][numComp]{};
1720 
1721  //store totMob upwinded by PPU for later schemes
1722  real64 totMob{};
1723  real64 dTotMob_dP[numFluxSupportPoints]{};
1724  real64 dTotMob_dC[numFluxSupportPoints][numComp]{};
1725  localIndex k_up_ppu = -1;
1726 
1727  //unelegant but need dummy when forming PPU total velocity
1728  real64 dummy[numComp];
1729  real64 dDummy_dP[numFluxSupportPoints][numComp];
1730  real64 dDummy_dC[numFluxSupportPoints][numComp][numComp];
1731 
1732 
1733  for( integer jp = 0; jp < numPhase; ++jp )
1734  {
1735  PPUPhaseFlux::compute( numPhase, jp, hasCapPressure,
1736  seri, sesri, sei,
1737  trans, dTrans_dPres,
1738  pres, gravCoef,
1739  phaseMob, dPhaseMob,
1740  phaseVolFrac, dPhaseVolFrac,
1741  phaseCompFrac, dPhaseCompFrac,
1742  dCompFrac_dCompDens,
1743  phaseMassDens, dPhaseMassDens,
1744  phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac,
1745  k_up_ppu, potGrad,
1746  phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC,
1747  dummy, dDummy_dP, dDummy_dC );
1748 
1749  totFlux += phaseFlux;
1750 
1751  phaseFlux = 0.;
1752  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
1753  {
1754  dTotFlux_dP[ke] += dPhaseFlux_dP[ke];
1755  totMob += phaseMob[seri[ke]][sesri[ke]][sei[ke]][jp];
1756  dTotMob_dP[ke] += dPhaseMob[seri[ke]][sesri[ke]][sei[ke]][jp][Deriv::dP];
1757  dPhaseFlux_dP[ke] = 0.;
1758 
1759  for( localIndex jc = 0; jc < numComp; ++jc )
1760  {
1761  dTotFlux_dC[ke][jc] += dPhaseFlux_dC[ke][jc];
1762  dTotMob_dC[ke][jc] += dPhaseMob[seri[ke]][sesri[ke]][sei[ke]][jp][Deriv::dC + jc];
1763  dPhaseFlux_dC[ke][jc] = 0.;
1764  }
1765  }
1766  }
1767 
1768  //fractional flow loop with IHU
1769  //maybe needed to have density out for upwinding
1770 
1771  // choose upstream cell
1772  // create local work arrays
1773  real64 viscousPhaseFlux{};
1774  real64 dViscousPhaseFlux_dP[numFluxSupportPoints]{};
1775  real64 dViscousPhaseFlux_dC[numFluxSupportPoints][numComp]{};
1776 
1777  real64 fractionalFlow{};
1778  real64 dFractionalFlow_dP[numFluxSupportPoints]{};
1779  real64 dFractionalFlow_dC[numFluxSupportPoints][numComp]{};
1780 
1781  // and the fractional flow for viscous part as \lambda_i^{up}/\sum_{NP}(\lambda_j^{up}) with up decided upon
1782  // the Upwind strategy
1783  UpwindHelpers::computeFractionalFlowViscous< numComp, numFluxSupportPoints,
1784  UPWIND_SCHEME >( numPhase,
1785  ip,
1786  seri,
1787  sesri,
1788  sei,
1789  trans,
1790  dTrans_dPres,
1791  k_up_ppu,
1792  totFlux,
1793  totMob,
1794  dTotMob_dP,
1795  dTotMob_dC,
1796  pres,
1797  gravCoef,
1798  dCompFrac_dCompDens,
1799  phaseMassDens,
1800  dPhaseMassDens,
1801  phaseMob,
1802  dPhaseMob,
1803  dPhaseVolFrac,
1804  phaseCapPressure,
1805  dPhaseCapPressure_dPhaseVolFrac,
1806  hasCapPressure,
1807  k_up,
1808  fractionalFlow,
1809  dFractionalFlow_dP,
1810  dFractionalFlow_dC );
1811 
1812 
1814  viscousPhaseFlux = fractionalFlow * totFlux;
1815  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
1816  {
1817  dViscousPhaseFlux_dP[ke] += dFractionalFlow_dP[ke] * totFlux;
1818 
1819 
1820  for( localIndex jc = 0; jc < numComp; ++jc )
1821  {
1822  dViscousPhaseFlux_dC[ke][jc] += dFractionalFlow_dC[ke][jc] * totFlux;
1823  }
1824  }
1825 
1826  //NON-FIXED UT -- to be canceled out if considered fixed
1827  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
1828  {
1829  dViscousPhaseFlux_dP[ke] += fractionalFlow * dTotFlux_dP[ke];
1830 
1831 
1832  for( localIndex jc = 0; jc < numComp; ++jc )
1833  {
1834  dViscousPhaseFlux_dC[ke][jc] += fractionalFlow * dTotFlux_dC[ke][jc];
1835  }
1836  }
1837  //distribute on phaseComponentFlux here
1838  PhaseComponentFlux::compute( ip, k_up,
1839  seri, sesri, sei,
1840  phaseCompFrac, dPhaseCompFrac, dCompFrac_dCompDens,
1841  viscousPhaseFlux, dViscousPhaseFlux_dP, dViscousPhaseFlux_dC,
1842  compFlux, dCompFlux_dP, dCompFlux_dC );
1843 
1844  // accumulate in the flux and its derivatives
1845  phaseFlux += viscousPhaseFlux;
1846  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
1847  {
1848  dPhaseFlux_dP[ke] += dViscousPhaseFlux_dP[ke];
1849 
1850 
1851  for( localIndex ic = 0; ic < numComp; ++ic )
1852  dPhaseFlux_dC[ke][ic] += dViscousPhaseFlux_dC[ke][ic];
1853  }
1854 
1856  localIndex k_up_g = -1;
1857  localIndex k_up_og = -1;
1858 
1859  real64 gravitationalPhaseFlux{};
1860  real64 gravitationalPhaseFlux_dP[numFluxSupportPoints]{};
1861  real64 gravitationalPhaseFlux_dC[numFluxSupportPoints][numComp]{};
1862 
1863  UpwindHelpers::computePotentialFluxesGravity< numComp,
1864  numFluxSupportPoints, UPWIND_SCHEME >(
1865  numPhase,
1866  ip,
1867  seri,
1868  sesri,
1869  sei,
1870  trans,
1871  dTrans_dPres,
1872  k_up_ppu,
1873  totFlux,
1874  totMob,
1875  dTotMob_dP,
1876  dTotMob_dC,
1877  pres,
1878  gravCoef,
1879  phaseMob,
1880  dPhaseMob,
1881  phaseVolFrac,
1882  dPhaseVolFrac,
1883  dCompFrac_dCompDens,
1884  phaseMassDens,
1885  dPhaseMassDens,
1886  phaseCapPressure,
1887  dPhaseCapPressure_dPhaseVolFrac,
1888  hasCapPressure,
1889  k_up_g,
1890  k_up_og,
1891  gravitationalPhaseFlux,
1892  gravitationalPhaseFlux_dP,
1893  gravitationalPhaseFlux_dC );
1894 
1895 
1896 
1897  //distribute on phaseComponentFlux here
1898  PhaseComponentFlux::compute( ip, k_up_g,
1899  seri, sesri, sei,
1900  phaseCompFrac, dPhaseCompFrac, dCompFrac_dCompDens,
1901  gravitationalPhaseFlux, gravitationalPhaseFlux_dP, gravitationalPhaseFlux_dC,
1902  compFlux, dCompFlux_dP, dCompFlux_dC );
1903 
1904 
1905  //update phaseFlux from gravitational
1906  phaseFlux += gravitationalPhaseFlux;
1907  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
1908  {
1909  dPhaseFlux_dP[ke] += gravitationalPhaseFlux_dP[ke];
1910  for( localIndex ic = 0; ic < numComp; ++ic )
1911  dPhaseFlux_dC[ke][ic] += gravitationalPhaseFlux_dC[ke][ic];
1912  }
1913 
1914 
1915  if( hasCapPressure )
1916  {
1918  localIndex k_up_pc = -1;
1919  localIndex k_up_opc = -1;
1920 
1921  real64 capillaryPhaseFlux{};
1922  real64 capillaryPhaseFlux_dP[numFluxSupportPoints]{};
1923  real64 capillaryPhaseFlux_dC[numFluxSupportPoints][numComp]{};
1924 
1925  UpwindHelpers::computePotentialFluxesCapillary< numComp,
1926  numFluxSupportPoints, UPWIND_SCHEME >(
1927  numPhase,
1928  ip,
1929  seri,
1930  sesri,
1931  sei,
1932  trans,
1933  dTrans_dPres,
1934  k_up_ppu,
1935  totFlux,
1936  totMob,
1937  dTotMob_dP,
1938  dTotMob_dC,
1939  pres,
1940  gravCoef,
1941  phaseMob,
1942  dPhaseMob,
1943  dPhaseVolFrac,
1944  dCompFrac_dCompDens,
1945  phaseMassDens,
1946  dPhaseMassDens,
1947  phaseCapPressure,
1948  dPhaseCapPressure_dPhaseVolFrac,
1949  hasCapPressure,
1950  k_up_pc,
1951  k_up_opc,
1952  capillaryPhaseFlux,
1953  capillaryPhaseFlux_dP,
1954  capillaryPhaseFlux_dC );
1955 
1956  //distribute on phaseComponentFlux here
1957  PhaseComponentFlux::compute( ip, k_up_pc,
1958  seri, sesri, sei,
1959  phaseCompFrac, dPhaseCompFrac, dCompFrac_dCompDens,
1960  capillaryPhaseFlux, capillaryPhaseFlux_dP, capillaryPhaseFlux_dC,
1961  compFlux, dCompFlux_dP, dCompFlux_dC );
1962 
1963 
1964  //update phaseFlux from capillary
1965  phaseFlux += capillaryPhaseFlux;
1966  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
1967  {
1968  dPhaseFlux_dP[ke] += capillaryPhaseFlux_dP[ke];
1969  for( localIndex ic = 0; ic < numComp; ++ic )
1970  dPhaseFlux_dC[ke][ic] += capillaryPhaseFlux_dC[ke][ic];
1971 
1972  }
1973 
1974  }//end if cappres
1975 
1976  }
1977 
1978 };
1979 
1980 } // namespace isothermalCompositionalMultiPhaseFVMKernelUtilities
1981 
1982 } // namespace geos
1983 
1984 
1985 #endif // GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_IHUPHASEFLUX_HPP
GEOS_HOST_DEVICE
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
GEOS_UNUSED_PARAM
#define GEOS_UNUSED_PARAM(X)
Mark an unused argument and silence compiler warnings.
Definition: GeosxMacros.hpp:72
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::UpwindHelpers::computePotentialFluxesGravity
static GEOS_HOST_DEVICE void computePotentialFluxesGravity(localIndex const numPhase, localIndex const ip, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], real64 const (&transmissibility)[2], real64 const (&dTrans_dPres)[2], localIndex const &k_up_ppu, real64 const totFlux, real64 const totMob, real64 const (&dTotMob_dP)[numFluxSupportPoints], real64 const (&dTotMob_dC)[numFluxSupportPoints][numComp], ElementViewConst< arrayView1d< real64 const > > const &pres, ElementViewConst< arrayView1d< real64 const > > const &gravCoef, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseMob, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &dPhaseMob, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseVolFrac, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const &dCompFrac_dCompDens, ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const &phaseMassDens, ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const &dPhaseMassDens, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const &phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const &dPhaseCapPressure_dPhaseVolFrac, localIndex const capPressureFlag, localIndex(&k_up), localIndex(&k_up_o), real64 &phaseFlux, real64(&dPhaseFlux_dP)[numFluxSupportPoints], real64(&dPhaseFlux_dC)[numFluxSupportPoints][numComp])
Form potential-related parts of fluxes.
Definition: IHUPhaseFlux.hpp:836
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::HybridUpwind
Class describing the Hybrid Upwind scheme as defined in "Consistent upwinding for sequential fully im...
Definition: IHUPhaseFlux.hpp:1465
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::UpwindScheme::getUpwindDirectionViscous
GEOS_HOST_DEVICE void getUpwindDirectionViscous(localIndex const numPhase, localIndex const ip, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], real64 const (&transmissibility)[2], real64 const (&dTrans_dPres)[2], real64 const totFlux, ElementViewConst< arrayView1d< real64 const > > const &pres, ElementViewConst< arrayView1d< real64 const > > const &gravCoef, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseMob, ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const &dCompFrac_dCompDens, ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const &phaseMassDens, ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const &dPhaseMassDens, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const &phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const &dPhaseCapPressure_dPhaseVolFrac, integer const capPressureFlag, localIndex &upwindDir)
Definition: IHUPhaseFlux.hpp:1247
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::UpwindScheme::getUpwindDirectionGravity
GEOS_HOST_DEVICE void getUpwindDirectionGravity(localIndex const numPhase, localIndex const ip, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], real64 const (&transmissibility)[2], real64 const (&dTrans_dPres)[2], real64 const totFlux, ElementViewConst< arrayView1d< real64 const > > const &pres, ElementViewConst< arrayView1d< real64 const > > const &gravCoef, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseMob, ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const &dCompFrac_dCompDens, ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const &phaseMassDens, ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const &dPhaseMassDens, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseVolFrac, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const &phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const &dPhaseCapPressure_dPhaseVolFrac, integer const capPressureFlag, localIndex &upwindDir)
Definition: IHUPhaseFlux.hpp:1299
geos::arrayView1d
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:180
geos::arrayView5d
ArrayView< T, 5, USD > arrayView5d
Alias for 5D array view.
Definition: DataTypes.hpp:244
geos::real64
double real64
64-bit floating point type.
Definition: DataTypes.hpp:99
geos::localIndex
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:85
geos::integer
std::int32_t integer
Signed integer type.
Definition: DataTypes.hpp:82
geos::arrayView4d
ArrayView< T, 4, USD > arrayView4d
Alias for 4D array view.
Definition: DataTypes.hpp:228
geos::arrayView2d
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:196
geos::arrayView3d
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
Definition: DataTypes.hpp:212
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::IHUPhaseFlux::compute
static GEOS_HOST_DEVICE void compute(integer const numPhase, integer const ip, integer const hasCapPressure, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], real64 const (&trans)[2], real64 const (&dTrans_dPres)[2], ElementViewConst< arrayView1d< real64 const > > const &pres, ElementViewConst< arrayView1d< real64 const > > const &gravCoef, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseMob, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &dPhaseMob, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseVolFrac, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &dPhaseVolFrac, ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const &phaseCompFrac, ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const &dPhaseCompFrac, ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const &dCompFrac_dCompDens, ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const &phaseMassDens, ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const &dPhaseMassDens, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const &phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const &dPhaseCapPressure_dPhaseVolFrac, localIndex &k_up, real64 &potGrad, real64(&phaseFlux), real64(&dPhaseFlux_dP)[numFluxSupportPoints], real64(&dPhaseFlux_dC)[numFluxSupportPoints][numComp], real64(&compFlux)[numComp], real64(&dCompFlux_dP)[numFluxSupportPoints][numComp], real64(&dCompFlux_dC)[numFluxSupportPoints][numComp][numComp])
Form the Implicit Hybrid Upwind from pressure gradient and gravitational head.
Definition: IHUPhaseFlux.hpp:1685
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::PPUPhaseFlux::compute
static GEOS_HOST_DEVICE void compute(integer const numPhase, integer const ip, integer const hasCapPressure, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], real64 const (&trans)[2], real64 const (&dTrans_dPres)[2], ElementViewConst< arrayView1d< real64 const > > const &pres, ElementViewConst< arrayView1d< real64 const > > const &gravCoef, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseMob, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &dPhaseMob, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseVolFrac, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &dPhaseVolFrac, ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const &phaseCompFrac, ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const &dPhaseCompFrac, ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const &dCompFrac_dCompDens, ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const &phaseMassDens, ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const &dPhaseMassDens, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const &phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const &dPhaseCapPressure_dPhaseVolFrac, localIndex &k_up, real64 &potGrad, real64(&phaseFlux), real64(&dPhaseFlux_dP)[numFluxSupportPoints], real64(&dPhaseFlux_dC)[numFluxSupportPoints][numComp], real64(&compFlux)[numComp], real64(&dCompFlux_dP)[numFluxSupportPoints][numComp], real64(&dCompFlux_dC)[numFluxSupportPoints][numComp][numComp])
Form the PhasePotentialUpwind from pressure gradient and gravitational head.
Definition: PPUPhaseFlux.hpp:75
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::PhaseComponentFlux::compute
static GEOS_HOST_DEVICE void compute(localIndex const ip, localIndex const k_up, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const &phaseCompFrac, ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const &dPhaseCompFrac, ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const &dCompFrac_dCompDens, real64 const &phaseFlux, real64 const (&dPhaseFlux_dP)[numFluxSupportPoints], real64 const (&dPhaseFlux_dC)[numFluxSupportPoints][numComp], real64(&compFlux)[numComp], real64(&dCompFlux_dP)[numFluxSupportPoints][numComp], real64(&dCompFlux_dC)[numFluxSupportPoints][numComp][numComp])
Compute the component flux for a given phase.
Definition: PhaseComponentFlux.hpp:64
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::UpwindHelpers::computePotentialCapillary::compute
static GEOS_HOST_DEVICE void compute(localIndex const numPhase, localIndex const ip, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], real64 const (&transmissibility)[2], real64 const (&dTrans_dPres)[2], real64 const GEOS_UNUSED_PARAM(totFlux), ElementViewConst< arrayView1d< real64 const > > const &GEOS_UNUSED_PARAM(gravCoef), ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const &GEOS_UNUSED_PARAM(dCompFrac_dCompDens), ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const &GEOS_UNUSED_PARAM(phaseMassDens), ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const &GEOS_UNUSED_PARAM(dPhaseMassDens), ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const &phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const &dPhaseCapPressure_dPhaseVolFrac, real64 &pot, real64(&dPot_dPres)[numFluxSupportPoints], real64(&dPot_dComp)[numFluxSupportPoints][numComp], real64(&GEOS_UNUSED_PARAM(dProp_dComp))[numComp])
Definition: IHUPhaseFlux.hpp:781
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::UpwindHelpers::computePotentialGravity::compute
static GEOS_HOST_DEVICE void compute(localIndex const GEOS_UNUSED_PARAM(numPhase), localIndex const ip, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], real64 const (&transmissibility)[2], real64 const (&dTrans_dPres)[2], real64 const GEOS_UNUSED_PARAM(totFlux), ElementViewConst< arrayView1d< real64 const > > const &gravCoef, ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const &dCompFrac_dCompDens, ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const &phaseMassDens, ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const &dPhaseMassDens, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseVolFrac, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &GEOS_UNUSED_PARAM(dPhaseVolFrac), ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const &GEOS_UNUSED_PARAM(phaseCapPressure), ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const &GEOS_UNUSED_PARAM(dPhaseCapPressure_dPhaseVolFrac), real64 &pot, real64(&dPot_dPres)[numFluxSupportPoints], real64(&dPot_dComp)[numFluxSupportPoints][numComp], real64(&dProp_dComp)[numComp])
Definition: IHUPhaseFlux.hpp:659
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::UpwindHelpers::computePotentialViscous
Struct defining formation of potential from different Physics (flagged by enum type T) to be used in ...
Definition: IHUPhaseFlux.hpp:613
Generated by doxygen 1.9.1