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