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 template< localIndex numComp >
46 GEOS_HOST_DEVICE
47 static void assignToZero( real64 & deriv_dP, real64 ( & deriv_dC )[numComp] )
48 {
49  deriv_dP = 0.0;
50  for( localIndex ic = 0; ic < numComp; ++ic )
51  {
52  deriv_dC[ic] = 0.0;
53  }
54 }
55 
56 template< localIndex numComp >
57 GEOS_HOST_DEVICE
58 static void assignToZero( real64 & value, real64 & deriv_dP, real64 ( & deriv_dC )[numComp] )
59 {
60  value = 0.0;
61  assignToZero( deriv_dP, deriv_dC );
62 }
63 
64 template< localIndex numComp, localIndex numFluxSupportPoints >
65 GEOS_HOST_DEVICE
66 static void assignToZero( real64 & value, real64 ( & deriv_dP )[numFluxSupportPoints], real64 ( & deriv_dC )[numFluxSupportPoints][numComp] )
67 {
68  value = 0;
69  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
70  {
71  assignToZero( deriv_dP[ke], deriv_dC[ke] );
72  }
73 }
74 
75 template< localIndex numComp >
76 GEOS_HOST_DEVICE
77 static void addToDerivativesScaled( real64 const ( &dDeriv_dP ), real64 const ( &dDeriv_dC )[numComp],
78  real64 const & factor,
79  real64 ( &deriv_dP ), real64 ( & deriv_dC )[numComp] )
80 {
81  deriv_dP += dDeriv_dP * factor;
82  for( localIndex ic = 0; ic < numComp; ++ic )
83  deriv_dC[ic] += dDeriv_dC[ic] * factor;
84 }
85 
86 template< localIndex numComp, localIndex numFluxSupportPoints >
87 GEOS_HOST_DEVICE
88 static void addToDerivativesScaled( real64 const ( &dDeriv_dP )[numFluxSupportPoints], real64 const ( &dDeriv_dC )[numFluxSupportPoints][numComp],
89  real64 const & factor,
90  real64 ( & deriv_dP )[numFluxSupportPoints], real64 ( & deriv_dC )[numFluxSupportPoints][numComp] )
91 {
92  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
93  {
94  addToDerivativesScaled( dDeriv_dP[ke], dDeriv_dC[ke], factor, deriv_dP[ke], deriv_dC[ke] );
95  }
96 }
97 
98 template< localIndex numComp >
99 GEOS_HOST_DEVICE
100 static void addToDerivatives( real64 const ( &dDeriv_dP ), real64 const ( &dDeriv_dC )[numComp],
101  real64 ( &deriv_dP ), real64 ( & deriv_dC )[numComp] )
102 {
103  addToDerivativesScaled( dDeriv_dP, dDeriv_dC, 1.0, deriv_dP, deriv_dC );
104 }
105 
106 template< localIndex numComp, localIndex numFluxSupportPoints >
107 GEOS_HOST_DEVICE
108 static void addToDerivatives( real64 const ( &dDeriv_dP )[numFluxSupportPoints], real64 const ( &dDeriv_dC )[numFluxSupportPoints][numComp],
109  real64 ( & deriv_dP )[numFluxSupportPoints], real64 ( & deriv_dC )[numFluxSupportPoints][numComp] )
110 {
111  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
112  {
113  addToDerivativesScaled( dDeriv_dP, dDeriv_dC, 1.0, deriv_dP, deriv_dC );
114  }
115 }
116 
117 template< localIndex numComp >
118 GEOS_HOST_DEVICE
119 static void addToValueAndDerivatives( real64 const & dValue, real64 const ( &dDeriv_dP ), real64 const ( &dDeriv_dC )[numComp],
120  real64 & value, real64 ( &deriv_dP ), real64 ( & deriv_dC )[numComp] )
121 {
122  value += dValue;
123  addToDerivatives( dDeriv_dP, dDeriv_dC, deriv_dP, deriv_dC );
124 }
125 
126 template< localIndex numComp, localIndex numFluxSupportPoints >
127 GEOS_HOST_DEVICE
128 static void addToValueAndDerivatives( real64 const & dValue, real64 const ( &dDeriv_dP )[numFluxSupportPoints], real64 const ( &dDeriv_dC )[numFluxSupportPoints][numComp],
129  real64 & value, real64 ( & deriv_dP )[numFluxSupportPoints], real64 ( & deriv_dC )[numFluxSupportPoints][numComp] )
130 {
131  value += dValue;
132  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
133  {
134  addToDerivatives( dDeriv_dP[ke], dDeriv_dC[ke], deriv_dP[ke], deriv_dC[ke] );
135  }
136 }
137 
138 template< localIndex numComp, localIndex numFluxSupportPoints >
139 GEOS_HOST_DEVICE
140 static void assignMobilityAndDerivatives( localIndex const & ip, localIndex const & upwindDir,
141  localIndex const (&seri)[numFluxSupportPoints],
142  localIndex const (&sesri)[numFluxSupportPoints],
143  localIndex const (&sei)[numFluxSupportPoints],
144  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
145  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
146  real64 & mobility, real64 & dMobility_dP, real64 ( & dMobility_dC )[numComp] )
147 {
148  localIndex const er_up = seri[upwindDir];
149  localIndex const esr_up = sesri[upwindDir];
150  localIndex const ei_up = sei[upwindDir];
151 
152  if( std::fabs( phaseMob[er_up][esr_up][ei_up][ip] ) > 1e-20 )
153  {
154  mobility += phaseMob[er_up][esr_up][ei_up][ip];
155  dMobility_dP += dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dP];
156  for( localIndex ic = 0; ic < numComp; ++ic )
157  {
158  dMobility_dC[ic] += dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dC + ic];
159  }
160  }
161 }
162 
163 template< localIndex numComp, localIndex numFluxSupportPoints >
164 GEOS_HOST_DEVICE
165 static void computeFractionalFlowAndDerivatives( localIndex const & k_up, real64 const & mob, real64 const & dMob_dP, real64 const ( &dMob_dC )[numComp],
166  real64 const & totMob, real64 const ( &dTotMob_dP )[numFluxSupportPoints], real64 const ( &dTotMob_dC )[numFluxSupportPoints][numComp],
167  real64 & fractionalFlow, real64 ( & dFractionalFlow_dP )[numFluxSupportPoints], real64 ( & dFractionalFlow_dC )[numFluxSupportPoints][numComp] )
168 {
169  // guard against no flow region
170  // fractional flow too low to let the upstream phase flow
171  if( std::fabs( mob ) > 1e-20 )
172  {
173  real64 const invTotMob = 1 / totMob;
174 
175  fractionalFlow = mob * invTotMob;
176 
177  addToDerivativesScaled( dMob_dP, dMob_dC, invTotMob, dFractionalFlow_dP[k_up], dFractionalFlow_dC[k_up] );
178 
179  addToDerivativesScaled( dTotMob_dP, dTotMob_dC, -fractionalFlow * invTotMob, dFractionalFlow_dP, dFractionalFlow_dC );
180  }
181 }
182 
183 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
184 GEOS_HOST_DEVICE
185 static void
186 upwindMobilityViscous( localIndex const numPhase,
187  localIndex const ip,
188  localIndex const (&seri)[numFluxSupportPoints],
189  localIndex const (&sesri)[numFluxSupportPoints],
190  localIndex const (&sei)[numFluxSupportPoints],
191  real64 const (&transmissibility)[2],
192  real64 const (&dTrans_dPres)[2],
193  real64 const totFlux, // in fine should be a ElemnetViewConst once seq form are in place
194  ElementViewConst< arrayView1d< real64 const > > const & pres,
195  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
196  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
197  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
198  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
199  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
200  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
201  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
202  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
203  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
204  integer const hasCapPressure,
205  localIndex & upwindDir,
206  real64 & mobility,
207  real64 (&dMobility_dP),
208  real64 (& dMobility_dC)[numComp] )
209 {
210  UPWIND scheme;
211  scheme.template getUpwindDirectionViscous< numComp, numFluxSupportPoints, UPWIND >( numPhase,
212  ip,
213  seri,
214  sesri,
215  sei,
216  transmissibility,
217  dTrans_dPres,
218  totFlux,
219  pres,
220  gravCoef,
221  phaseMob,
222  dCompFrac_dCompDens,
223  phaseMassDens,
224  dPhaseMassDens,
225  dPhaseVolFrac,
226  phaseCapPressure,
227  dPhaseCapPressure_dPhaseVolFrac,
228  hasCapPressure,
229  upwindDir );
230 
231  //reinit
232  assignToZero( mobility, dMobility_dP, dMobility_dC );
233  assignMobilityAndDerivatives( ip, upwindDir, seri, sesri, sei, phaseMob, dPhaseMob, mobility, dMobility_dP, dMobility_dC );
234 }
235 
236 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
237 GEOS_HOST_DEVICE
238 static void
239 upwindMobilityGravity( localIndex const numPhase,
240  localIndex const ip,
241  localIndex const (&seri)[numFluxSupportPoints],
242  localIndex const (&sesri)[numFluxSupportPoints],
243  localIndex const (&sei)[numFluxSupportPoints],
244  real64 const (&transmissibility)[2],
245  real64 const (&dTrans_dPres)[2],
246  real64 const totFlux, //in fine should be a ElemnetViewConst once seq form are in place
247  ElementViewConst< arrayView1d< real64 const > > const & pres,
248  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
249  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
250  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
251  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
252  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
253  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
254  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
255  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
256  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
257  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
258  integer const hasCapPressure,
259  integer const checkPhasePresenceInGravity,
260  localIndex & upwindDir,
261  real64 & mobility,
262  real64 ( &dMobility_dP ),
263  real64 ( & dMobility_dC )[numComp] )
264 {
265  UPWIND scheme;
266  scheme.template getUpwindDirectionGravity< numComp, numFluxSupportPoints, UPWIND >( numPhase,
267  ip,
268  seri,
269  sesri,
270  sei,
271  transmissibility,
272  dTrans_dPres,
273  totFlux,
274  pres,
275  gravCoef,
276  phaseMob,
277  dCompFrac_dCompDens,
278  phaseMassDens,
279  dPhaseMassDens,
280  phaseVolFrac,
281  dPhaseVolFrac,
282  phaseCapPressure,
283  dPhaseCapPressure_dPhaseVolFrac,
284  hasCapPressure,
285  checkPhasePresenceInGravity,
286  upwindDir );
287 
288  //reinit
289  assignToZero( mobility, dMobility_dP, dMobility_dC );
290  assignMobilityAndDerivatives( ip, upwindDir, seri, sesri, sei, phaseMob, dPhaseMob, mobility, dMobility_dP, dMobility_dC );
291 }
292 
293 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
294 GEOS_HOST_DEVICE
295 static void
296 upwindMobilityCapillary( localIndex const numPhase,
297  localIndex const ip,
298  localIndex const (&seri)[numFluxSupportPoints],
299  localIndex const (&sesri)[numFluxSupportPoints],
300  localIndex const (&sei)[numFluxSupportPoints],
301  real64 const (&transmissibility)[2],
302  real64 const (&dTrans_dPres)[2],
303  real64 const totFlux, //in fine should be a ElemnetViewConst once seq form are in place
304  ElementViewConst< arrayView1d< real64 const > > const & pres,
305  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
306  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
307  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
308  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
309  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
310  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
311  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
312  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
313  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
314  integer const hasCapPressure,
315  localIndex & upwindDir,
316  real64 & mobility,
317  real64 ( &dMobility_dP ),
318  real64 ( & dMobility_dC )[numComp] )
319 {
320  UPWIND scheme;
321  scheme.template getUpwindDirectionCapillary< numComp, numFluxSupportPoints, UPWIND >( numPhase,
322  ip,
323  seri,
324  sesri,
325  sei,
326  transmissibility,
327  dTrans_dPres,
328  totFlux,
329  pres,
330  gravCoef,
331  phaseMob,
332  dCompFrac_dCompDens,
333  phaseMassDens,
334  dPhaseMassDens,
335  dPhaseVolFrac,
336  phaseCapPressure,
337  dPhaseCapPressure_dPhaseVolFrac,
338  hasCapPressure,
339  upwindDir );
340 
341  //reinit
342  assignToZero( mobility, dMobility_dP, dMobility_dC );
343  assignMobilityAndDerivatives( ip, upwindDir, seri, sesri, sei, phaseMob, dPhaseMob, mobility, dMobility_dP, dMobility_dC );
344 }
345 
346 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
347 GEOS_HOST_DEVICE
348 static void
349 computeFractionalFlowViscous( localIndex const numPhase,
350  localIndex const ip,
351  localIndex const (&seri)[numFluxSupportPoints],
352  localIndex const (&sesri)[numFluxSupportPoints],
353  localIndex const (&sei)[numFluxSupportPoints],
354  real64 const (&transmissibility)[2],
355  real64 const (&dTrans_dPres)[2],
356  real64 const totFlux,
357  real64 const totMob,
358  real64 const (&dTotMob_dP)[numFluxSupportPoints],
359  real64 const (&dTotMob_dC)[numFluxSupportPoints][numComp],
360  ElementViewConst< arrayView1d< real64 const > > const & pres,
361  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
362  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
363  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
364  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
365  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
366  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
367  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
368  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
369  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
370  integer const hasCapPressure,
371  real64 & fractionalFlow,
372  real64 (& dFractionalFlow_dP)[numFluxSupportPoints],
373  real64 (& dFractionalFlow_dC)[numFluxSupportPoints][numComp] )
374 {
375  localIndex k_up;
376  real64 mob{};
377  real64 dMob_dP{};
378  real64 dMob_dC[numComp]{};
379 
380  upwindMobilityViscous< numComp, numFluxSupportPoints, UPWIND >( numPhase,
381  ip,
382  seri,
383  sesri,
384  sei,
385  transmissibility,
386  dTrans_dPres,
387  totFlux,
388  pres,
389  gravCoef,
390  dCompFrac_dCompDens,
391  phaseMassDens,
392  dPhaseMassDens,
393  phaseMob,
394  dPhaseMob,
395  dPhaseVolFrac,
396  phaseCapPressure,
397  dPhaseCapPressure_dPhaseVolFrac,
398  hasCapPressure,
399  k_up,
400  mob,
401  dMob_dP,
402  dMob_dC );
403 
404  // reinit
405  assignToZero( fractionalFlow, dFractionalFlow_dP, dFractionalFlow_dC );
406  computeFractionalFlowAndDerivatives( k_up, mob, dMob_dP, dMob_dC,
407  totMob, dTotMob_dP, dTotMob_dC,
408  fractionalFlow, dFractionalFlow_dP, dFractionalFlow_dC );
409 }
410 
411 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
412 GEOS_HOST_DEVICE
413 static void
414 computeFractionalFlowGravity( localIndex const numPhase,
415  localIndex const ip,
416  localIndex const (&seri)[numFluxSupportPoints],
417  localIndex const (&sesri)[numFluxSupportPoints],
418  localIndex const (&sei)[numFluxSupportPoints],
419  real64 const (&transmissibility)[2],
420  real64 const (&dTrans_dPres)[2],
421  real64 const totFlux,
422  real64 const totMob,
423  real64 const (&dTotMob_dP)[numFluxSupportPoints],
424  real64 const (&dTotMob_dC)[numFluxSupportPoints][numComp],
425  ElementViewConst< arrayView1d< real64 const > > const & pres,
426  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
427  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
428  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
429  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
430  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
431  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
432  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
433  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
434  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
435  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
436  integer const hasCapPressure,
437  integer const checkPhasePresenceInGravity,
438  real64 & fractionalFlow,
439  real64 (& dFractionalFlow_dP)[numFluxSupportPoints],
440  real64 (& dFractionalFlow_dC)[numFluxSupportPoints][numComp] )
441 {
442  localIndex k_up;
443  real64 mob{};
444  real64 dMob_dP{};
445  real64 dMob_dC[numComp]{};
446 
447  upwindMobilityGravity< numComp, numFluxSupportPoints, UPWIND >( numPhase,
448  ip,
449  seri,
450  sesri,
451  sei,
452  transmissibility,
453  dTrans_dPres,
454  totFlux,
455  pres,
456  gravCoef,
457  dCompFrac_dCompDens,
458  phaseMassDens,
459  dPhaseMassDens,
460  phaseMob,
461  dPhaseMob,
462  phaseVolFrac,
463  dPhaseVolFrac,
464  phaseCapPressure,
465  dPhaseCapPressure_dPhaseVolFrac,
466  hasCapPressure,
467  checkPhasePresenceInGravity,
468  k_up,
469  mob,
470  dMob_dP,
471  dMob_dC );
472 
473  // reinit
474  assignToZero( fractionalFlow, dFractionalFlow_dP, dFractionalFlow_dC );
475  computeFractionalFlowAndDerivatives( k_up, mob, dMob_dP, dMob_dC,
476  totMob, dTotMob_dP, dTotMob_dC,
477  fractionalFlow, dFractionalFlow_dP, dFractionalFlow_dC );
478 }
479 
480 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
481 GEOS_HOST_DEVICE
482 static void
483 computeFractionalFlowCapillary( localIndex const numPhase,
484  localIndex const ip,
485  localIndex const (&seri)[numFluxSupportPoints],
486  localIndex const (&sesri)[numFluxSupportPoints],
487  localIndex const (&sei)[numFluxSupportPoints],
488  real64 const (&transmissibility)[2],
489  real64 const (&dTrans_dPres)[2],
490  real64 const totFlux,
491  real64 const totMob,
492  real64 const (&dTotMob_dP)[numFluxSupportPoints],
493  real64 const (&dTotMob_dC)[numFluxSupportPoints][numComp],
494  ElementViewConst< arrayView1d< real64 const > > const & pres,
495  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
496  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
497  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
498  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
499  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
500  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
501  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
502  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
503  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
504  integer const hasCapPressure,
505  real64 & fractionalFlow,
506  real64 ( & dFractionalFlow_dP)[numFluxSupportPoints],
507  real64 ( & dFractionalFlow_dC)[numFluxSupportPoints][numComp]
508  )
509 {
510  localIndex k_up;
511  real64 mob{};
512  real64 dMob_dP{};
513  real64 dMob_dC[numComp]{};
514 
515  upwindMobilityCapillary< numComp, numFluxSupportPoints, UPWIND >( numPhase,
516  ip,
517  seri,
518  sesri,
519  sei,
520  transmissibility,
521  dTrans_dPres,
522  totFlux,
523  pres,
524  gravCoef,
525  dCompFrac_dCompDens,
526  phaseMassDens,
527  dPhaseMassDens,
528  phaseMob,
529  dPhaseMob,
530  dPhaseVolFrac,
531  phaseCapPressure,
532  dPhaseCapPressure_dPhaseVolFrac,
533  hasCapPressure,
534  k_up,
535  mob,
536  dMob_dP,
537  dMob_dC );
538 
539  // reinit
540  assignToZero( fractionalFlow, dFractionalFlow_dP, dFractionalFlow_dC );
541  computeFractionalFlowAndDerivatives( k_up, mob, dMob_dP, dMob_dC,
542  totMob, dTotMob_dP, dTotMob_dC,
543  fractionalFlow, dFractionalFlow_dP, dFractionalFlow_dC );
544 }
545 
553 struct computePotentialViscous
554 {
555  template< localIndex numComp, localIndex numFluxSupportPoints >
556  GEOS_HOST_DEVICE
557  static void compute( localIndex const GEOS_UNUSED_PARAM( numPhase ),
558  localIndex const GEOS_UNUSED_PARAM( ip ),
559  localIndex const (&GEOS_UNUSED_PARAM( seri ))[numFluxSupportPoints],
560  localIndex const (&GEOS_UNUSED_PARAM( sesri ))[numFluxSupportPoints],
561  localIndex const (&GEOS_UNUSED_PARAM( sei ))[numFluxSupportPoints],
562  real64 const (&GEOS_UNUSED_PARAM( transmissibility ))[2],
563  real64 const (&GEOS_UNUSED_PARAM( dTrans_dPres ))[2],
564  real64 const totFlux,
565  ElementViewConst< arrayView1d< real64 const > > const & GEOS_UNUSED_PARAM( gravCoef ),
566  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const &
567  GEOS_UNUSED_PARAM( dCompFrac_dCompDens ),
568  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const &
569  GEOS_UNUSED_PARAM( phaseMassDens ),
570  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const &
571  GEOS_UNUSED_PARAM( dPhaseMassDens ),
572  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &
573  GEOS_UNUSED_PARAM( dPhaseVolFrac ),
574  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const &
575  GEOS_UNUSED_PARAM( phaseCapPressure ),
576  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const &
577  GEOS_UNUSED_PARAM( dPhaseCapPressure_dPhaseVolFrac ),
578  real64 & pot,
579  real64( &GEOS_UNUSED_PARAM( dPot_dPres ))[numFluxSupportPoints],
580  real64( &GEOS_UNUSED_PARAM( dPot_dComp ))[numFluxSupportPoints][numComp] )
581  {
582  pot = totFlux;
583  //could be relevant for symmetry to include derivative
584  }
585 };
586 
589 struct computePotentialGravity
590 {
595  template< localIndex numComp, localIndex numFluxSupportPoints >
596  GEOS_HOST_DEVICE
597  static void compute( localIndex const GEOS_UNUSED_PARAM( numPhase ),
598  localIndex const ip,
599  integer const checkPhasePresenceInGravity,
600  localIndex const (&seri)[numFluxSupportPoints],
601  localIndex const (&sesri)[numFluxSupportPoints],
602  localIndex const (&sei)[numFluxSupportPoints],
603  real64 const (&transmissibility)[2],
604  real64 const (&dTrans_dPres)[2],
605  real64 const GEOS_UNUSED_PARAM( totFlux ),
606  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
607  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
608  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
609  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
610  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
611  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &
612  GEOS_UNUSED_PARAM( dPhaseVolFrac ),
613  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const &
614  GEOS_UNUSED_PARAM( phaseCapPressure ),
615  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const &
616  GEOS_UNUSED_PARAM( dPhaseCapPressure_dPhaseVolFrac ),
617  real64 & pot,
618  real64 ( & dPot_dPres )[numFluxSupportPoints],
619  real64 ( & dPot_dComp )[numFluxSupportPoints][numComp] )
620  {
621  //init
622  assignToZero( pot, dPot_dPres, dPot_dComp );
623 
624  //working arrays
625  real64 densMean{};
626  real64 dDensMean_dPres[numFluxSupportPoints]{};
627  real64 dDensMean_dComp[numFluxSupportPoints][numComp]{};
628  isothermalCompositionalMultiphaseFVMKernels::helpers::
629  calculateMeanDensity( ip, seri, sesri, sei,
630  checkPhasePresenceInGravity,
631  phaseVolFrac, dCompFrac_dCompDens,
632  phaseMassDens, dPhaseMassDens,
633  densMean, dDensMean_dPres, dDensMean_dComp );
634 
635  // compute potential difference MPFA-style
636  for( localIndex i = 0; i < numFluxSupportPoints; ++i )
637  {
638  localIndex const er = seri[i];
639  localIndex const esr = sesri[i];
640  localIndex const ei = sei[i];
641 
642  real64 const gravD = transmissibility[i] * gravCoef[er][esr][ei];
643  real64 const dGravD_dP = dTrans_dPres[i] * gravCoef[er][esr][ei];
644  pot += densMean * gravD;
645  dPot_dPres[i] += densMean * dGravD_dP;
646 
647  // need to add contributions from both cells the mean density depends on
648  addToDerivativesScaled( dDensMean_dPres, dDensMean_dComp, gravD, dPot_dPres, dPot_dComp );
649  }
650  }
651 
652 };
653 
656 struct computePotentialCapillary
657 {
662  template< localIndex numComp, localIndex numFluxSupportPoints >
663  GEOS_HOST_DEVICE
664  static void compute( localIndex const numPhase,
665  localIndex const ip,
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 & GEOS_UNUSED_PARAM( gravCoef ),
673  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const &
674  GEOS_UNUSED_PARAM( dCompFrac_dCompDens ),
675  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const &
676  GEOS_UNUSED_PARAM( phaseMassDens ),
677  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const &
678  GEOS_UNUSED_PARAM( dPhaseMassDens ),
679  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
680  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
681  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
682  real64 & pot,
683  real64 ( & dPot_dPres)[numFluxSupportPoints],
684  real64 (& dPot_dComp)[numFluxSupportPoints][numComp] )
685  {
686  //init
687  assignToZero( pot, dPot_dPres, dPot_dComp );
688 
689  for( localIndex i = 0; i < numFluxSupportPoints; ++i )
690  {
691  localIndex const er = seri[i];
692  localIndex const esr = sesri[i];
693  localIndex const ei = sei[i];
694 
695  real64 const capPres = phaseCapPressure[er][esr][ei][0][ip];
696 
697  pot += transmissibility[i] * capPres;
698 
699  dPot_dPres[i] += dTrans_dPres[i] * capPres;
700 
701  // need to add contributions from both cells
702  for( localIndex jp = 0; jp < numPhase; ++jp )
703  {
704  real64 const dCapPressure_dS = dPhaseCapPressure_dPhaseVolFrac[er][esr][ei][0][ip][jp];
705 
706  dPot_dPres[i] += transmissibility[i] * dCapPressure_dS * dPhaseVolFrac[er][esr][ei][jp][Deriv::dP];
707 
708  for( localIndex jc = 0; jc < numComp; ++jc )
709  {
710  dPot_dComp[i][jc] += transmissibility[i] * dCapPressure_dS * dPhaseVolFrac[er][esr][ei][jp][Deriv::dC + jc];
711  }
712  }
713  }
714  }
715 };
716 
718 
719 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
720 GEOS_HOST_DEVICE
721 static void computePotentialFluxesGravity( localIndex const numPhase,
722  localIndex const ip,
723  localIndex const (&seri)[numFluxSupportPoints],
724  localIndex const (&sesri)[numFluxSupportPoints],
725  localIndex const (&sei)[numFluxSupportPoints],
726  real64 const (&transmissibility)[2],
727  real64 const (&dTrans_dPres)[2],
728  real64 const totFlux,
729  real64 const totMob,
730  real64 const (&dTotMob_dP)[numFluxSupportPoints],
731  real64 const (&dTotMob_dC)[numFluxSupportPoints][numComp],
732  ElementViewConst< arrayView1d< real64 const > > const & pres,
733  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
734  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
735  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
736  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
737  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
738  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
739  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
740  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
741  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
742  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
743  localIndex const hasCapPressure,
744  integer const checkPhasePresenceInGravity,
745  real64 & phaseFlux,
746  real64 (& dPhaseFlux_dP)[numFluxSupportPoints],
747  real64 ( & dPhaseFlux_dC)[numFluxSupportPoints][numComp] )
748 {
749 
750  real64 fflow{};
751  real64 dFflow_dP[numFluxSupportPoints]{};
752  real64 dFflow_dC[numFluxSupportPoints][numComp]{};
753 
754  real64 pot{};
755  real64 dPot_dP[numFluxSupportPoints]{};
756  real64 dPot_dC[numFluxSupportPoints][numComp]{};
757 
758  //
759  UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >( numPhase,
760  ip,
761  checkPhasePresenceInGravity,
762  seri,
763  sesri,
764  sei,
765  transmissibility,
766  dTrans_dPres,
767  totFlux,
768  gravCoef,
769  dCompFrac_dCompDens,
770  phaseMassDens,
771  dPhaseMassDens,
772  phaseVolFrac,
773  dPhaseVolFrac,
774  phaseCapPressure,
775  dPhaseCapPressure_dPhaseVolFrac,
776  pot,
777  dPot_dP,
778  dPot_dC );
779 
780  // and the fractional flow for gravitational part as \lambda_i^{up}/\sum_{numPhase}(\lambda_k^{up}) with up decided upon
781  // the Upwind strategy
782  UpwindHelpers::computeFractionalFlowGravity< numComp, numFluxSupportPoints, UPWIND >( numPhase,
783  ip,
784  seri,
785  sesri,
786  sei,
787  transmissibility,
788  dTrans_dPres,
789  totFlux,
790  totMob,
791  dTotMob_dP,
792  dTotMob_dC,
793  pres,
794  gravCoef,
795  dCompFrac_dCompDens,
796  phaseMassDens,
797  dPhaseMassDens,
798  phaseMob,
799  dPhaseMob,
800  phaseVolFrac,
801  dPhaseVolFrac,
802  phaseCapPressure,
803  dPhaseCapPressure_dPhaseVolFrac,
804  hasCapPressure,
805  checkPhasePresenceInGravity,
806  fflow,
807  dFflow_dP,
808  dFflow_dC );
809 
810 
811  for( localIndex jp = 0; jp < numPhase; ++jp )
812  {
813  if( ip != jp )
814  {
815 
816  real64 potOther{};
817  real64 dPotOther_dP[numFluxSupportPoints]{};
818  real64 dPotOther_dC[numFluxSupportPoints][numComp]{};
819 
820  //Fetch pot for phase j!=i defined as \rho_j g dz/dx
821  UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >( numPhase,
822  jp,
823  checkPhasePresenceInGravity,
824  seri,
825  sesri,
826  sei,
827  transmissibility,
828  dTrans_dPres,
829  totFlux,
830  gravCoef,
831  dCompFrac_dCompDens,
832  phaseMassDens,
833  dPhaseMassDens,
834  phaseVolFrac,
835  dPhaseVolFrac,
836  phaseCapPressure,
837  dPhaseCapPressure_dPhaseVolFrac,
838  potOther,
839  dPotOther_dP,
840  dPotOther_dC );
841 
842  //Eventually get the mobility of the second phase
843  localIndex k_up_o = -1;
844  real64 mobOther{};
845  real64 dMobOther_dP{};
846  real64 dMobOther_dC[numComp]{};
847 
848  // and the other mobility for gravitational part as \lambda_j^{up} with up decided upon
849  // the Upwind strategy - Note that it should be the same as the gravitational fractional flow
850 
851  UpwindHelpers::upwindMobilityGravity< numComp, numFluxSupportPoints, UPWIND >( numPhase,
852  jp,
853  seri,
854  sesri,
855  sei,
856  transmissibility,
857  dTrans_dPres,
858  totFlux,
859  pres,
860  gravCoef,
861  dCompFrac_dCompDens,
862  phaseMassDens,
863  dPhaseMassDens,
864  phaseMob,
865  dPhaseMob,
866  phaseVolFrac,
867  dPhaseVolFrac,
868  phaseCapPressure,
869  dPhaseCapPressure_dPhaseVolFrac,
870  hasCapPressure,
871  checkPhasePresenceInGravity,
872  k_up_o,
873  mobOther,
874  dMobOther_dP,
875  dMobOther_dC );
876 
877 
878  // Assembling gravitational flux phase-wise as \phi_{i,g} = \sum_{k\nei} \lambda_k^{up,g} f_k^{up,g} (G_i - G_k)
879  phaseFlux -= fflow * mobOther * (pot - potOther);
880 
881  // mobOther part
882  UpwindHelpers::addToDerivativesScaled( dMobOther_dP, dMobOther_dC, -fflow * (pot - potOther), dPhaseFlux_dP[k_up_o], dPhaseFlux_dC[k_up_o] );
883 
884  // mob related part of dFflow_dP is only upstream defined but totMob related is defined everywhere
885  UpwindHelpers::addToDerivativesScaled( dFflow_dP, dFflow_dC, -mobOther * (pot - potOther), dPhaseFlux_dP, dPhaseFlux_dC );
886 
887  // potential difference part
888  UpwindHelpers::addToDerivativesScaled( dPot_dP, dPot_dC, -fflow * mobOther, dPhaseFlux_dP, dPhaseFlux_dC );
889  UpwindHelpers::addToDerivativesScaled( dPotOther_dP, dPotOther_dC, fflow * mobOther, dPhaseFlux_dP, dPhaseFlux_dC );
890  }
891  }
892 
893 }
894 
895 template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
896 GEOS_HOST_DEVICE
897 static void computePotentialFluxesCapillary( localIndex const numPhase,
898  localIndex const ip,
899  localIndex const (&seri)[numFluxSupportPoints],
900  localIndex const (&sesri)[numFluxSupportPoints],
901  localIndex const (&sei)[numFluxSupportPoints],
902  real64 const (&transmissibility)[2],
903  real64 const (&dTrans_dPres)[2],
904  real64 const totFlux,
905  real64 const totMob,
906  real64 const (&dTotMob_dP)[numFluxSupportPoints],
907  real64 const (&dTotMob_dC)[numFluxSupportPoints][numComp],
908  ElementViewConst< arrayView1d< real64 const > > const & pres,
909  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
910  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
911  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
912  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
913  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
914  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
915  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
916  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
917  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
918  localIndex const hasCapPressure,
919  real64 & phaseFlux,
920  real64 (& dPhaseFlux_dP)[numFluxSupportPoints],
921  real64 ( & dPhaseFlux_dC)[numFluxSupportPoints][numComp] )
922 {
923  real64 fflow{};
924  real64 dFflow_dP[numFluxSupportPoints]{};
925  real64 dFflow_dC[numFluxSupportPoints][numComp]{};
926 
927  real64 pot{};
928  real64 dPot_dP[numFluxSupportPoints]{};
929  real64 dPot_dC[numFluxSupportPoints][numComp]{};
930 
931  computePotentialCapillary::compute< numComp, numFluxSupportPoints >( numPhase,
932  ip,
933  seri,
934  sesri,
935  sei,
936  transmissibility,
937  dTrans_dPres,
938  totFlux,
939  gravCoef,
940  dCompFrac_dCompDens,
941  phaseMassDens,
942  dPhaseMassDens,
943  dPhaseVolFrac,
944  phaseCapPressure,
945  dPhaseCapPressure_dPhaseVolFrac,
946  pot,
947  dPot_dP,
948  dPot_dC );
949 
950  // and the fractional flow for gravitational part as \lambda_i^{up}/\sum_{numPhase}(\lambda_k^{up}) with up decided upon
951  // the Upwind strategy
952  UpwindHelpers::computeFractionalFlowCapillary< numComp, numFluxSupportPoints, UPWIND >( numPhase,
953  ip,
954  seri,
955  sesri,
956  sei,
957  transmissibility,
958  dTrans_dPres,
959  totFlux,
960  totMob,
961  dTotMob_dP,
962  dTotMob_dC,
963  pres,
964  gravCoef,
965  dCompFrac_dCompDens,
966  phaseMassDens,
967  dPhaseMassDens,
968  phaseMob,
969  dPhaseMob,
970  dPhaseVolFrac,
971  phaseCapPressure,
972  dPhaseCapPressure_dPhaseVolFrac,
973  hasCapPressure,
974  fflow,
975  dFflow_dP,
976  dFflow_dC );
977 
978 
979  for( localIndex jp = 0; jp < numPhase; ++jp )
980  {
981  if( ip != jp )
982  {
983  real64 potOther{};
984  real64 dPotOther_dP[numFluxSupportPoints]{};
985  real64 dPotOther_dC[numFluxSupportPoints][numComp]{};
986 
987  //Fetch pot for phase j!=i defined as \rho_j g dz/dx
988  computePotentialCapillary::compute< numComp, numFluxSupportPoints >( numPhase,
989  jp,
990  seri,
991  sesri,
992  sei,
993  transmissibility,
994  dTrans_dPres,
995  totFlux,
996  gravCoef,
997  dCompFrac_dCompDens,
998  phaseMassDens,
999  dPhaseMassDens,
1000  dPhaseVolFrac,
1001  phaseCapPressure,
1002  dPhaseCapPressure_dPhaseVolFrac,
1003  potOther,
1004  dPotOther_dP,
1005  dPotOther_dC );
1006 
1007  //Eventually get the mobility of the second phase
1008  localIndex k_up_o = -1;
1009  real64 mobOther{};
1010  real64 dMobOther_dP{};
1011  real64 dMobOther_dC[numComp]{};
1012 
1013  // and the other mobility for gravitational part as \lambda_j^{up} with up decided upon
1014  // the Upwind strategy - Note that it should be the same as the gravitational fractional flow
1015 
1016  UpwindHelpers::upwindMobilityCapillary< numComp, numFluxSupportPoints, UPWIND >( numPhase,
1017  jp,
1018  seri,
1019  sesri,
1020  sei,
1021  transmissibility,
1022  dTrans_dPres,
1023  totFlux,
1024  pres,
1025  gravCoef,
1026  dCompFrac_dCompDens,
1027  phaseMassDens,
1028  dPhaseMassDens,
1029  phaseMob,
1030  dPhaseMob,
1031  dPhaseVolFrac,
1032  phaseCapPressure,
1033  dPhaseCapPressure_dPhaseVolFrac,
1034  hasCapPressure,
1035  k_up_o,
1036  mobOther,
1037  dMobOther_dP,
1038  dMobOther_dC );
1039 
1040 
1041  // Assembling gravitational flux phase-wise as \phi_{i,g} = \sum_{k\nei} \lambda_k^{up,g} f_k^{up,g} (G_i - G_k)
1042  phaseFlux -= fflow * mobOther * (pot - potOther);
1043 
1044  // mobOther part
1045  UpwindHelpers::addToDerivativesScaled( dMobOther_dP, dMobOther_dC, -fflow * (pot - potOther), dPhaseFlux_dP[k_up_o], dPhaseFlux_dC[k_up_o] );
1046 
1047  // mob related part of dFflow_dP is only upstream defined but totMob related is defined everywhere
1048  UpwindHelpers::addToDerivativesScaled( dFflow_dP, dFflow_dC, -mobOther * (pot - potOther), dPhaseFlux_dP, dPhaseFlux_dC );
1049 
1050  // potential difference part
1051  UpwindHelpers::addToDerivativesScaled( dPot_dP, dPot_dC, -fflow * mobOther, dPhaseFlux_dP, dPhaseFlux_dC );
1052  UpwindHelpers::addToDerivativesScaled( dPotOther_dP, dPotOther_dC, fflow * mobOther, dPhaseFlux_dP, dPhaseFlux_dC );
1053  }
1054  }
1055 }
1056 
1057 }//end of struct UpwindHelpers
1058 
1059 /************************* UPWIND ******************/
1060 
1065 class UpwindScheme
1066 {
1067 
1068 public:
1069 
1070  //default ctor
1071  UpwindScheme() = default;
1072 
1073  //usual copy ctor
1074  UpwindScheme( UpwindScheme const & scheme ) = default;
1075 
1076  //default move ctor
1077  UpwindScheme( UpwindScheme && ) = default;
1078 
1079  //deleted copy and move assignement
1080  UpwindScheme & operator=( UpwindScheme const & ) = delete;
1081 
1082  UpwindScheme & operator=( UpwindScheme && ) = delete;
1083 
1084  GEOS_HOST_DEVICE
1085  virtual ~UpwindScheme() = default;
1086 
1087  template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
1088  GEOS_HOST_DEVICE
1089  inline
1090  void getUpwindDirectionViscous( localIndex const numPhase,
1091  localIndex const ip,
1092  localIndex const (&seri)[numFluxSupportPoints],
1093  localIndex const (&sesri)[numFluxSupportPoints],
1094  localIndex const (&sei)[numFluxSupportPoints],
1095  real64 const (&transmissibility)[2],
1096  real64 const (&dTrans_dPres)[2],
1097  real64 const totFlux, //in fine should be a ElemnetViewConst once seq form are in place
1098  ElementViewConst< arrayView1d< real64 const > > const & pres,
1099  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1100  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1101  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1102  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1103  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1104  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1105  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1106  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1107  integer const hasCapPressure,
1108  localIndex & upwindDir
1109  )
1110  {
1111  real64 pot{};
1114  UPWIND::template computePotentialViscous< numComp, numFluxSupportPoints >( numPhase,
1115  ip,
1116  seri,
1117  sesri,
1118  sei,
1119  transmissibility,
1120  dTrans_dPres,
1121  totFlux,
1122  pres,
1123  gravCoef,
1124  phaseMob,
1125  dCompFrac_dCompDens,
1126  phaseMassDens,
1127  dPhaseMassDens,
1128  dPhaseVolFrac,
1129  phaseCapPressure,
1130  dPhaseCapPressure_dPhaseVolFrac,
1131  hasCapPressure,
1132  pot );
1133  //all definition has been changed to fit pot>0 => first cell is upstream
1134  upwindDir = (pot > 0) ? 0 : 1;
1135  }
1136 
1137 
1138  template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
1139  GEOS_HOST_DEVICE
1140  void getUpwindDirectionGravity( localIndex const numPhase,
1141  localIndex const ip,
1142  localIndex const (&seri)[numFluxSupportPoints],
1143  localIndex const (&sesri)[numFluxSupportPoints],
1144  localIndex const (&sei)[numFluxSupportPoints],
1145  real64 const (&transmissibility)[2],
1146  real64 const (&dTrans_dPres)[2],
1147  real64 const totFlux,
1148  ElementViewConst< arrayView1d< real64 const > > const & pres,
1149  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1150  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1151  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1152  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1153  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1154  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
1155  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1156  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1157  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1158  integer const hasCapPressure,
1159  integer const checkPhasePresenceInGravity,
1160  localIndex & upwindDir
1161  )
1162  {
1163  real64 pot{};
1164 
1167  UPWIND::template computePotentialGravity< numComp, numFluxSupportPoints >( numPhase,
1168  ip,
1169  seri,
1170  sesri,
1171  sei,
1172  transmissibility,
1173  dTrans_dPres,
1174  totFlux,
1175  pres,
1176  gravCoef,
1177  phaseMob,
1178  dCompFrac_dCompDens,
1179  phaseMassDens,
1180  dPhaseMassDens,
1181  phaseVolFrac,
1182  dPhaseVolFrac,
1183  phaseCapPressure,
1184  dPhaseCapPressure_dPhaseVolFrac,
1185  hasCapPressure,
1186  checkPhasePresenceInGravity,
1187  pot );
1188 
1189  //all definition has been changed to fit pot>0 => first cell is upstream
1190  upwindDir = (pot >= 0) ? 0 : 1;
1191  }
1192 
1193 
1194  template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
1195  GEOS_HOST_DEVICE
1196  void getUpwindDirectionCapillary( localIndex const numPhase,
1197  localIndex const ip,
1198  localIndex const (&seri)[numFluxSupportPoints],
1199  localIndex const (&sesri)[numFluxSupportPoints],
1200  localIndex const (&sei)[numFluxSupportPoints],
1201  real64 const (&transmissibility)[2],
1202  real64 const (&dTrans_dPres)[2],
1203  real64 const totFlux, //in fine should be a ElemnetViewConst once seq form are in place
1204  ElementViewConst< arrayView1d< real64 const > > const & pres,
1205  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1206  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1207  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1208  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1209  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1210  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1211  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1212  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1213  integer const hasCapPressure,
1214  localIndex & upwindDir
1215  )
1216  {
1217  real64 pot{};
1218 
1219  // each derived concrete class has to define a computePotential method that is calling UpwindScheme::potential method with a specific
1220  // lamda defining how to get these potentials
1221  UPWIND::template computePotentialCapillary< numComp, numFluxSupportPoints >( numPhase,
1222  ip,
1223  seri,
1224  sesri,
1225  sei,
1226  transmissibility,
1227  dTrans_dPres,
1228  totFlux,
1229  pres,
1230  gravCoef,
1231  phaseMob,
1232  dCompFrac_dCompDens,
1233  phaseMassDens,
1234  dPhaseMassDens,
1235  dPhaseVolFrac,
1236  phaseCapPressure,
1237  dPhaseCapPressure_dPhaseVolFrac,
1238  hasCapPressure,
1239  pot );
1240 
1241  //all definition has been changed to fit pot>0 => first cell is upstream
1242  upwindDir = (pot >= 0) ? 0 : 1;
1243  }
1244 
1245 
1246 
1247  // templated way of evaluating the potential (to the exception of viscous one) which relies on
1248  // up-or-downwinded mobility terms pre-multiplying potential differences
1249  template< localIndex numComp, localIndex numFluxSupportPoints, typename LAMBDA >
1250  GEOS_HOST_DEVICE
1251  static void potential( localIndex numPhase,
1252  localIndex ip,
1253  localIndex const (&seri)[numFluxSupportPoints],
1254  localIndex const (&sesri)[numFluxSupportPoints],
1255  localIndex const (&sei)[numFluxSupportPoints],
1256  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1257  real64 & weightedPotential,
1258  LAMBDA && fn )
1259  {
1260  //getPhase Pot
1261  real64 pot{};
1262  real64 pot_dP[numFluxSupportPoints]{};
1263  real64 pot_dC[numFluxSupportPoints][numComp]{};
1264 
1265  fn( ip, pot, pot_dP, pot_dC );
1266 
1267  localIndex const k_up = 0;
1268  localIndex const k_dw = 1;
1269 
1270  //loop other other phases to form
1271  for( localIndex jp = 0; jp < numPhase; ++jp )
1272  {
1273  if( jp != ip )
1274  {
1275  localIndex const er_up = seri[k_up];
1276  localIndex const esr_up = sesri[k_up];
1277  localIndex const ei_up = sei[k_up];
1278 
1279  localIndex const er_dw = seri[k_dw];
1280  localIndex const esr_dw = sesri[k_dw];
1281  localIndex const ei_dw = sei[k_dw];
1282 
1283  real64 potOther{};
1284  real64 potOther_dP[numFluxSupportPoints]{};
1285  real64 potOther_dC[numFluxSupportPoints][numComp]{};
1286 
1287  fn( jp, potOther, potOther_dP, potOther_dC );
1288 
1289  real64 const mob_up = phaseMob[er_up][esr_up][ei_up][jp];
1290  real64 const mob_dw = phaseMob[er_dw][esr_dw][ei_dw][jp];
1291 
1292  weightedPotential += (pot - potOther >= 0) ? mob_dw * (potOther - pot) : mob_up * (potOther - pot);
1293 
1294  }
1295  }
1296  }
1297 
1298 };
1299 
1304 class HybridUpwind : public UpwindScheme
1305 {
1306 
1307 public:
1308  template< localIndex numComp, localIndex numFluxSupportPoints >
1309  GEOS_HOST_DEVICE
1310  static
1311  void computePotentialViscous( localIndex const numPhase,
1312  localIndex const ip,
1313  localIndex const (&seri)[numFluxSupportPoints],
1314  localIndex const (&sesri)[numFluxSupportPoints],
1315  localIndex const (&sei)[numFluxSupportPoints],
1316  real64 const (&transmissibility)[2],
1317  real64 const (&dTrans_dPres)[2],
1318  real64 const totalFlux,
1319  ElementViewConst< arrayView1d< real64 const > > const & GEOS_UNUSED_PARAM( pres ),
1320  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1321  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &
1322  GEOS_UNUSED_PARAM( phaseMob ),
1323  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1324  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1325  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1326  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1327  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1328  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1329  integer const GEOS_UNUSED_PARAM( hasCapPressure ),
1330  real64 & potential )
1331  {
1332  real64 dPot_dP[numFluxSupportPoints]{};
1333  real64 dPot_dC[numFluxSupportPoints][numComp]{};
1334  UpwindHelpers::computePotentialViscous::compute< numComp, numFluxSupportPoints >(
1335  numPhase,
1336  ip,
1337  seri,
1338  sesri,
1339  sei,
1340  transmissibility,
1341  dTrans_dPres,
1342  totalFlux,
1343  gravCoef,
1344  dCompFrac_dCompDens,
1345  phaseMassDens,
1346  dPhaseMassDens,
1347  dPhaseVolFrac,
1348  phaseCapPressure,
1349  dPhaseCapPressure_dPhaseVolFrac,
1350  potential,
1351  dPot_dP,
1352  dPot_dC );
1353  }
1354 
1355  template< localIndex numComp, localIndex numFluxSupportPoints >
1356  GEOS_HOST_DEVICE
1357  static
1358  void computePotentialGravity( localIndex const numPhase,
1359  localIndex const ip,
1360  localIndex const (&seri)[numFluxSupportPoints],
1361  localIndex const (&sesri)[numFluxSupportPoints],
1362  localIndex const (&sei)[numFluxSupportPoints],
1363  real64 const (&transmissibility)[2],
1364  real64 const (&dTrans_dPres)[2],
1365  real64 const totalFlux,
1366  ElementViewConst< arrayView1d< real64 const > > const & GEOS_UNUSED_PARAM( pres ),
1367  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1368  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1369  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1370  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1371  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1372  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
1373  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1374  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1375  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1376  integer const GEOS_UNUSED_PARAM( hasCapPressure ),
1377  integer const checkPhasePresenceInGravity,
1378  real64 & potential )
1379  {
1380 
1381  //Form total velocity
1382  potential = 0;
1383 
1384  //the arg lambda allows us to access some genericity
1385  UpwindScheme::template potential< numComp, numFluxSupportPoints >( numPhase, ip, seri, sesri, sei,
1386  phaseMob, potential,
1387  [&]( localIndex ipp,
1388  real64 & potential_,
1389  real64 (& dPotential_dP_)[numFluxSupportPoints],
1390  real64 (& dPotential_dC_)[numFluxSupportPoints][numComp] ) {
1391 
1392  UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >(
1393  numPhase,
1394  ipp,
1395  checkPhasePresenceInGravity,
1396  seri,
1397  sesri,
1398  sei,
1399  transmissibility,
1400  dTrans_dPres,
1401  totalFlux,
1402  gravCoef,
1403  dCompFrac_dCompDens,
1404  phaseMassDens,
1405  dPhaseMassDens,
1406  phaseVolFrac,
1407  dPhaseVolFrac,
1408  phaseCapPressure,
1409  dPhaseCapPressure_dPhaseVolFrac,
1410  potential_,
1411  dPotential_dP_,
1412  dPotential_dC_ );
1413 
1414  } );
1415  }
1416 
1417 
1418  template< localIndex numComp, localIndex numFluxSupportPoints >
1419  GEOS_HOST_DEVICE
1420  static
1421  void computePotentialCapillary( localIndex const numPhase,
1422  localIndex const ip,
1423  localIndex const (&seri)[numFluxSupportPoints],
1424  localIndex const (&sesri)[numFluxSupportPoints],
1425  localIndex const (&sei)[numFluxSupportPoints],
1426  real64 const (&transmissibility)[2],
1427  real64 const (&dTrans_dPres)[2],
1428  real64 const totalFlux,
1429  ElementViewConst< arrayView1d< real64 const > > const & GEOS_UNUSED_PARAM( pres ),
1430  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1431  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &
1432  phaseMob,
1433  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1434  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1435  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1436  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1437  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1438  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1439  integer const GEOS_UNUSED_PARAM( hasCapPressure ),
1440  real64 & potential )
1441 
1442  {
1443 
1444  //Form total velocity
1445  potential = 0;
1446 
1447  //the arg lambda allows us to access some genericity
1448  UpwindScheme::template potential< numComp, numFluxSupportPoints >( numPhase, ip, seri, sesri, sei,
1449  phaseMob, potential,
1450  [&]( localIndex ipp,
1451  real64 & potential_,
1452  real64 (& dPotential_dP_)[numFluxSupportPoints],
1453  real64 (& dPotential_dC_)[numFluxSupportPoints][numComp] )
1454  {
1455 
1456  UpwindHelpers::computePotentialCapillary::compute< numComp, numFluxSupportPoints >(
1457  numPhase,
1458  ipp,
1459  seri,
1460  sesri,
1461  sei,
1462  transmissibility,
1463  dTrans_dPres,
1464  totalFlux,
1465  gravCoef,
1466  dCompFrac_dCompDens,
1467  phaseMassDens,
1468  dPhaseMassDens,
1469  dPhaseVolFrac,
1470  phaseCapPressure,
1471  dPhaseCapPressure_dPhaseVolFrac,
1472  potential_,
1473  dPotential_dP_,
1474  dPotential_dC_ );
1475  } );
1476  }
1477 
1478 };
1479 
1480 /*** IHU ***/
1481 
1482 struct IHUPhaseFlux
1483 {
1484 
1485  using UPWIND_SCHEME = HybridUpwind;
1486 
1487  static constexpr double minTotMob = 1e-12;
1488 
1516  template< integer numComp, integer numFluxSupportPoints >
1517  GEOS_HOST_DEVICE
1518  static void
1519  compute( integer const numPhase,
1520  integer const ip,
1521  integer const hasCapPressure,
1522  integer const checkPhasePresenceInGravity,
1523  localIndex const ( &seri )[numFluxSupportPoints],
1524  localIndex const ( &sesri )[numFluxSupportPoints],
1525  localIndex const ( &sei )[numFluxSupportPoints],
1526  real64 const ( &trans )[2],
1527  real64 const ( &dTrans_dPres )[2],
1528  ElementViewConst< arrayView1d< real64 const > > const & pres,
1529  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1530  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1531  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
1532  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
1533  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
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< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1538  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1539  real64 & potGrad,
1540  real64 ( &phaseFlux ),
1541  real64 ( & dPhaseFlux_dP )[numFluxSupportPoints],
1542  real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp] )
1543  {
1544 
1545  //loop over all phases to form total velocity
1546  real64 totFlux{};
1547  real64 dTotFlux_dP[numFluxSupportPoints]{};
1548  real64 dTotFlux_dC[numFluxSupportPoints][numComp]{};
1549 
1550  //store totMob upwinded by PPU for later schemes
1551  real64 totMob{};
1552  real64 dTotMob_dP[numFluxSupportPoints]{};
1553  real64 dTotMob_dC[numFluxSupportPoints][numComp]{};
1554 
1555  //unelegant but need dummy when forming PPU total velocity
1556  real64 dPhaseFlux_dTrans;
1557 
1558  for( integer jp = 0; jp < numPhase; ++jp )
1559  {
1560  PPUPhaseFlux::compute( numPhase, jp,
1561  hasCapPressure, checkPhasePresenceInGravity,
1562  seri, sesri, sei,
1563  trans, dTrans_dPres,
1564  pres, gravCoef,
1565  phaseMob, dPhaseMob,
1566  phaseVolFrac, dPhaseVolFrac,
1567  dCompFrac_dCompDens,
1568  phaseMassDens, dPhaseMassDens,
1569  phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac,
1570  potGrad, phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, dPhaseFlux_dTrans );
1571 
1572  // accumulate into total flux
1573  UpwindHelpers::addToValueAndDerivatives( phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC,
1574  totFlux, dTotFlux_dP, dTotFlux_dC );
1575 /*
1576  // old wrong way:
1577  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
1578  {
1579  totMob += phaseMob[seri[ke]][sesri[ke]][sei[ke]][jp];
1580  dTotMob_dP[ke] += dPhaseMob[seri[ke]][sesri[ke]][sei[ke]][jp][Deriv::dP];
1581  for( localIndex jc = 0; jc < numComp; ++jc )
1582  {
1583  dTotMob_dC[ke][jc] += dPhaseMob[seri[ke]][sesri[ke]][sei[ke]][jp][Deriv::dC + jc];
1584  }
1585  }
1586  */
1587 
1588  // new correct way:
1589  // choose upstream cell
1590  localIndex const k_up = (phaseFlux >= 0) ? 0 : 1;
1591  // accumulate into total mobility
1592  UpwindHelpers::assignMobilityAndDerivatives( jp, k_up, seri, sesri, sei,
1593  phaseMob, dPhaseMob,
1594  totMob, dTotMob_dP[k_up], dTotMob_dC[k_up] );
1595  }
1596 
1597  // safeguard
1598  totMob = LvArray::math::max( totMob, minTotMob );
1599 
1600  // assign to zero
1601  UpwindHelpers::assignToZero( phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
1602 
1603  // fractional flow loop with IHU
1604  // maybe needed to have density out for upwinding
1605 
1606  // choose upstream cell
1607  // create local work arrays
1608  real64 viscousPhaseFlux{};
1609  real64 dViscousPhaseFlux_dP[numFluxSupportPoints]{};
1610  real64 dViscousPhaseFlux_dC[numFluxSupportPoints][numComp]{};
1611 
1612  real64 fractionalFlow{};
1613  real64 dFractionalFlow_dP[numFluxSupportPoints]{};
1614  real64 dFractionalFlow_dC[numFluxSupportPoints][numComp]{};
1615 
1616  // and the fractional flow for viscous part as \lambda_i^{up}/\sum_{NP}(\lambda_j^{up}) with up decided upon
1617  // the Upwind strategy
1618  UpwindHelpers::computeFractionalFlowViscous< numComp, numFluxSupportPoints,
1619  UPWIND_SCHEME >( numPhase,
1620  ip,
1621  seri,
1622  sesri,
1623  sei,
1624  trans,
1625  dTrans_dPres,
1626  totFlux,
1627  totMob,
1628  dTotMob_dP,
1629  dTotMob_dC,
1630  pres,
1631  gravCoef,
1632  dCompFrac_dCompDens,
1633  phaseMassDens,
1634  dPhaseMassDens,
1635  phaseMob,
1636  dPhaseMob,
1637  dPhaseVolFrac,
1638  phaseCapPressure,
1639  dPhaseCapPressure_dPhaseVolFrac,
1640  hasCapPressure,
1641  fractionalFlow,
1642  dFractionalFlow_dP,
1643  dFractionalFlow_dC );
1644 
1646  viscousPhaseFlux = fractionalFlow * totFlux;
1647 
1648  UpwindHelpers::addToDerivativesScaled( dFractionalFlow_dP, dFractionalFlow_dC,
1649  totFlux,
1650  dViscousPhaseFlux_dP, dViscousPhaseFlux_dC );
1651 
1652  //NON-FIXED UT -- to be canceled out if considered fixed
1653  UpwindHelpers::addToDerivativesScaled( dTotFlux_dP, dTotFlux_dC,
1654  fractionalFlow,
1655  dViscousPhaseFlux_dP, dViscousPhaseFlux_dC );
1656 
1657  // accumulate in the flux and its derivatives
1658  UpwindHelpers::addToValueAndDerivatives( viscousPhaseFlux, dViscousPhaseFlux_dP, dViscousPhaseFlux_dC,
1659  phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
1660 
1662 
1663  real64 gravitationalPhaseFlux{};
1664  real64 gravitationalPhaseFlux_dP[numFluxSupportPoints]{};
1665  real64 gravitationalPhaseFlux_dC[numFluxSupportPoints][numComp]{};
1666 
1667  UpwindHelpers::computePotentialFluxesGravity< numComp,
1668  numFluxSupportPoints, UPWIND_SCHEME >(
1669  numPhase,
1670  ip,
1671  seri,
1672  sesri,
1673  sei,
1674  trans,
1675  dTrans_dPres,
1676  totFlux,
1677  totMob,
1678  dTotMob_dP,
1679  dTotMob_dC,
1680  pres,
1681  gravCoef,
1682  phaseMob,
1683  dPhaseMob,
1684  phaseVolFrac,
1685  dPhaseVolFrac,
1686  dCompFrac_dCompDens,
1687  phaseMassDens,
1688  dPhaseMassDens,
1689  phaseCapPressure,
1690  dPhaseCapPressure_dPhaseVolFrac,
1691  hasCapPressure,
1692  checkPhasePresenceInGravity,
1693  gravitationalPhaseFlux,
1694  gravitationalPhaseFlux_dP,
1695  gravitationalPhaseFlux_dC );
1696 
1697  //update phaseFlux from gravitational
1698  UpwindHelpers::addToValueAndDerivatives( gravitationalPhaseFlux, gravitationalPhaseFlux_dP, gravitationalPhaseFlux_dC,
1699  phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
1700 
1701  if( hasCapPressure )
1702  {
1704 
1705  real64 capillaryPhaseFlux{};
1706  real64 capillaryPhaseFlux_dP[numFluxSupportPoints]{};
1707  real64 capillaryPhaseFlux_dC[numFluxSupportPoints][numComp]{};
1708 
1709  UpwindHelpers::computePotentialFluxesCapillary< numComp,
1710  numFluxSupportPoints, UPWIND_SCHEME >(
1711  numPhase,
1712  ip,
1713  seri,
1714  sesri,
1715  sei,
1716  trans,
1717  dTrans_dPres,
1718  totFlux,
1719  totMob,
1720  dTotMob_dP,
1721  dTotMob_dC,
1722  pres,
1723  gravCoef,
1724  phaseMob,
1725  dPhaseMob,
1726  dPhaseVolFrac,
1727  dCompFrac_dCompDens,
1728  phaseMassDens,
1729  dPhaseMassDens,
1730  phaseCapPressure,
1731  dPhaseCapPressure_dPhaseVolFrac,
1732  hasCapPressure,
1733  capillaryPhaseFlux,
1734  capillaryPhaseFlux_dP,
1735  capillaryPhaseFlux_dC );
1736 
1737  //update phaseFlux from capillary
1738  UpwindHelpers::addToValueAndDerivatives( capillaryPhaseFlux, capillaryPhaseFlux_dP, capillaryPhaseFlux_dC,
1739  phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
1740 
1741  }//end if cappres
1742  }
1743 
1744 };
1745 
1746 } // namespace isothermalCompositionalMultiPhaseFVMKernelUtilities
1747 
1748 } // namespace geos
1749 
1750 
1751 #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:97
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], 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, real64 &phaseFlux, real64(&dPhaseFlux_dP)[numFluxSupportPoints], real64(&dPhaseFlux_dC)[numFluxSupportPoints][numComp])
Form potential-related parts of fluxes.
Definition: IHUPhaseFlux.hpp:721
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::HybridUpwind
Class describing the Hybrid Upwind scheme as defined in "Consistent upwinding for sequential fully im...
Definition: IHUPhaseFlux.hpp:1305
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:1090
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:1140
geos::arrayView1d
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:179
geos::real64
double real64
64-bit floating point type.
Definition: DataTypes.hpp:98
geos::localIndex
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:84
geos::arrayView4d
ArrayView< T, 4, USD > arrayView4d
Alias for 4D array view.
Definition: DataTypes.hpp:227
geos::arrayView2d
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:195
geos::integer
int integer
Signed integer type.
Definition: DataTypes.hpp:81
geos::arrayView3d
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
Definition: DataTypes.hpp:211
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< 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, real64 &potGrad, real64(&phaseFlux), real64(&dPhaseFlux_dP)[numFluxSupportPoints], real64(&dPhaseFlux_dC)[numFluxSupportPoints][numComp])
Form the Implicit Hybrid Upwind from pressure gradient and gravitational head.
Definition: IHUPhaseFlux.hpp:1519
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< 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, real64 &potGrad, real64 &phaseFlux, real64(&dPhaseFlux_dP)[numFluxSupportPoints], real64(&dPhaseFlux_dC)[numFluxSupportPoints][numComp], real64 &dPhaseFlux_dTrans)
Form the PhasePotentialUpwind from pressure gradient and gravitational head.
Definition: PPUPhaseFlux.hpp:75
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])
Definition: IHUPhaseFlux.hpp:664
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])
Definition: IHUPhaseFlux.hpp:597
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::UpwindHelpers::computePotentialViscous
Struct defining formation of potential from different Physics (flagged by enum type T) to be used in ...
Definition: IHUPhaseFlux.hpp:554
Generated by doxygen 1.9.1