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  virtual ~UpwindScheme() = default;
1085 
1086  template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
1087  GEOS_HOST_DEVICE
1088  inline
1089  void getUpwindDirectionViscous( localIndex const numPhase,
1090  localIndex const ip,
1091  localIndex const (&seri)[numFluxSupportPoints],
1092  localIndex const (&sesri)[numFluxSupportPoints],
1093  localIndex const (&sei)[numFluxSupportPoints],
1094  real64 const (&transmissibility)[2],
1095  real64 const (&dTrans_dPres)[2],
1096  real64 const totFlux, //in fine should be a ElemnetViewConst once seq form are in place
1097  ElementViewConst< arrayView1d< real64 const > > const & pres,
1098  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1099  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1100  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1101  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1102  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1103  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1104  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1105  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1106  integer const hasCapPressure,
1107  localIndex & upwindDir
1108  )
1109  {
1110  real64 pot{};
1113  UPWIND::template computePotentialViscous< numComp, numFluxSupportPoints >( numPhase,
1114  ip,
1115  seri,
1116  sesri,
1117  sei,
1118  transmissibility,
1119  dTrans_dPres,
1120  totFlux,
1121  pres,
1122  gravCoef,
1123  phaseMob,
1124  dCompFrac_dCompDens,
1125  phaseMassDens,
1126  dPhaseMassDens,
1127  dPhaseVolFrac,
1128  phaseCapPressure,
1129  dPhaseCapPressure_dPhaseVolFrac,
1130  hasCapPressure,
1131  pot );
1132  //all definition has been changed to fit pot>0 => first cell is upstream
1133  upwindDir = (pot > 0) ? 0 : 1;
1134  }
1135 
1136 
1137  template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
1138  GEOS_HOST_DEVICE
1139  void getUpwindDirectionGravity( localIndex const numPhase,
1140  localIndex const ip,
1141  localIndex const (&seri)[numFluxSupportPoints],
1142  localIndex const (&sesri)[numFluxSupportPoints],
1143  localIndex const (&sei)[numFluxSupportPoints],
1144  real64 const (&transmissibility)[2],
1145  real64 const (&dTrans_dPres)[2],
1146  real64 const totFlux,
1147  ElementViewConst< arrayView1d< real64 const > > const & pres,
1148  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1149  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1150  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1151  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1152  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1153  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
1154  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1155  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1156  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1157  integer const hasCapPressure,
1158  integer const checkPhasePresenceInGravity,
1159  localIndex & upwindDir
1160  )
1161  {
1162  real64 pot{};
1163 
1166  UPWIND::template computePotentialGravity< numComp, numFluxSupportPoints >( numPhase,
1167  ip,
1168  seri,
1169  sesri,
1170  sei,
1171  transmissibility,
1172  dTrans_dPres,
1173  totFlux,
1174  pres,
1175  gravCoef,
1176  phaseMob,
1177  dCompFrac_dCompDens,
1178  phaseMassDens,
1179  dPhaseMassDens,
1180  phaseVolFrac,
1181  dPhaseVolFrac,
1182  phaseCapPressure,
1183  dPhaseCapPressure_dPhaseVolFrac,
1184  hasCapPressure,
1185  checkPhasePresenceInGravity,
1186  pot );
1187 
1188  //all definition has been changed to fit pot>0 => first cell is upstream
1189  upwindDir = (pot >= 0) ? 0 : 1;
1190  }
1191 
1192 
1193  template< localIndex numComp, localIndex numFluxSupportPoints, class UPWIND >
1194  GEOS_HOST_DEVICE
1195  void getUpwindDirectionCapillary( localIndex const numPhase,
1196  localIndex const ip,
1197  localIndex const (&seri)[numFluxSupportPoints],
1198  localIndex const (&sesri)[numFluxSupportPoints],
1199  localIndex const (&sei)[numFluxSupportPoints],
1200  real64 const (&transmissibility)[2],
1201  real64 const (&dTrans_dPres)[2],
1202  real64 const totFlux, //in fine should be a ElemnetViewConst once seq form are in place
1203  ElementViewConst< arrayView1d< real64 const > > const & pres,
1204  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1205  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1206  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1207  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1208  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1209  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1210  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1211  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1212  integer const hasCapPressure,
1213  localIndex & upwindDir
1214  )
1215  {
1216  real64 pot{};
1217 
1218  // each derived concrete class has to define a computePotential method that is calling UpwindScheme::potential method with a specific
1219  // lamda defining how to get these potentials
1220  UPWIND::template computePotentialCapillary< numComp, numFluxSupportPoints >( numPhase,
1221  ip,
1222  seri,
1223  sesri,
1224  sei,
1225  transmissibility,
1226  dTrans_dPres,
1227  totFlux,
1228  pres,
1229  gravCoef,
1230  phaseMob,
1231  dCompFrac_dCompDens,
1232  phaseMassDens,
1233  dPhaseMassDens,
1234  dPhaseVolFrac,
1235  phaseCapPressure,
1236  dPhaseCapPressure_dPhaseVolFrac,
1237  hasCapPressure,
1238  pot );
1239 
1240  //all definition has been changed to fit pot>0 => first cell is upstream
1241  upwindDir = (pot >= 0) ? 0 : 1;
1242  }
1243 
1244 
1245 
1246  // templated way of evaluating the potential (to the exception of viscous one) which relies on
1247  // up-or-downwinded mobility terms pre-multiplying potential differences
1248  template< localIndex numComp, localIndex numFluxSupportPoints, typename LAMBDA >
1249  GEOS_HOST_DEVICE
1250  static void potential( localIndex numPhase,
1251  localIndex ip,
1252  localIndex const (&seri)[numFluxSupportPoints],
1253  localIndex const (&sesri)[numFluxSupportPoints],
1254  localIndex const (&sei)[numFluxSupportPoints],
1255  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1256  real64 & weightedPotential,
1257  LAMBDA && fn )
1258  {
1259  //getPhase Pot
1260  real64 pot{};
1261  real64 pot_dP[numFluxSupportPoints]{};
1262  real64 pot_dC[numFluxSupportPoints][numComp]{};
1263 
1264  fn( ip, pot, pot_dP, pot_dC );
1265 
1266  localIndex const k_up = 0;
1267  localIndex const k_dw = 1;
1268 
1269  //loop other other phases to form
1270  for( localIndex jp = 0; jp < numPhase; ++jp )
1271  {
1272  if( jp != ip )
1273  {
1274  localIndex const er_up = seri[k_up];
1275  localIndex const esr_up = sesri[k_up];
1276  localIndex const ei_up = sei[k_up];
1277 
1278  localIndex const er_dw = seri[k_dw];
1279  localIndex const esr_dw = sesri[k_dw];
1280  localIndex const ei_dw = sei[k_dw];
1281 
1282  real64 potOther{};
1283  real64 potOther_dP[numFluxSupportPoints]{};
1284  real64 potOther_dC[numFluxSupportPoints][numComp]{};
1285 
1286  fn( jp, potOther, potOther_dP, potOther_dC );
1287 
1288  real64 const mob_up = phaseMob[er_up][esr_up][ei_up][jp];
1289  real64 const mob_dw = phaseMob[er_dw][esr_dw][ei_dw][jp];
1290 
1291  weightedPotential += (pot - potOther >= 0) ? mob_dw * (potOther - pot) : mob_up * (potOther - pot);
1292 
1293  }
1294  }
1295  }
1296 
1297 };
1298 
1303 class HybridUpwind : public UpwindScheme
1304 {
1305 
1306 public:
1307  template< localIndex numComp, localIndex numFluxSupportPoints >
1308  GEOS_HOST_DEVICE
1309  static
1310  void computePotentialViscous( localIndex const numPhase,
1311  localIndex const ip,
1312  localIndex const (&seri)[numFluxSupportPoints],
1313  localIndex const (&sesri)[numFluxSupportPoints],
1314  localIndex const (&sei)[numFluxSupportPoints],
1315  real64 const (&transmissibility)[2],
1316  real64 const (&dTrans_dPres)[2],
1317  real64 const totalFlux,
1318  ElementViewConst< arrayView1d< real64 const > > const & GEOS_UNUSED_PARAM( pres ),
1319  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1320  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &
1321  GEOS_UNUSED_PARAM( phaseMob ),
1322  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1323  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1324  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1325  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1326  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1327  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1328  integer const GEOS_UNUSED_PARAM( hasCapPressure ),
1329  real64 & potential )
1330  {
1331  real64 dPot_dP[numFluxSupportPoints]{};
1332  real64 dPot_dC[numFluxSupportPoints][numComp]{};
1333  UpwindHelpers::computePotentialViscous::compute< numComp, numFluxSupportPoints >(
1334  numPhase,
1335  ip,
1336  seri,
1337  sesri,
1338  sei,
1339  transmissibility,
1340  dTrans_dPres,
1341  totalFlux,
1342  gravCoef,
1343  dCompFrac_dCompDens,
1344  phaseMassDens,
1345  dPhaseMassDens,
1346  dPhaseVolFrac,
1347  phaseCapPressure,
1348  dPhaseCapPressure_dPhaseVolFrac,
1349  potential,
1350  dPot_dP,
1351  dPot_dC );
1352  }
1353 
1354  template< localIndex numComp, localIndex numFluxSupportPoints >
1355  GEOS_HOST_DEVICE
1356  static
1357  void computePotentialGravity( localIndex const numPhase,
1358  localIndex const ip,
1359  localIndex const (&seri)[numFluxSupportPoints],
1360  localIndex const (&sesri)[numFluxSupportPoints],
1361  localIndex const (&sei)[numFluxSupportPoints],
1362  real64 const (&transmissibility)[2],
1363  real64 const (&dTrans_dPres)[2],
1364  real64 const totalFlux,
1365  ElementViewConst< arrayView1d< real64 const > > const & GEOS_UNUSED_PARAM( pres ),
1366  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1367  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1368  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1369  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1370  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1371  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
1372  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1373  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1374  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1375  integer const GEOS_UNUSED_PARAM( hasCapPressure ),
1376  integer const checkPhasePresenceInGravity,
1377  real64 & potential )
1378  {
1379 
1380  //Form total velocity
1381  potential = 0;
1382 
1383  //the arg lambda allows us to access some genericity
1384  UpwindScheme::template potential< numComp, numFluxSupportPoints >( numPhase, ip, seri, sesri, sei,
1385  phaseMob, potential,
1386  [&]( localIndex ipp,
1387  real64 & potential_,
1388  real64 (& dPotential_dP_)[numFluxSupportPoints],
1389  real64 (& dPotential_dC_)[numFluxSupportPoints][numComp] ) {
1390 
1391  UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >(
1392  numPhase,
1393  ipp,
1394  checkPhasePresenceInGravity,
1395  seri,
1396  sesri,
1397  sei,
1398  transmissibility,
1399  dTrans_dPres,
1400  totalFlux,
1401  gravCoef,
1402  dCompFrac_dCompDens,
1403  phaseMassDens,
1404  dPhaseMassDens,
1405  phaseVolFrac,
1406  dPhaseVolFrac,
1407  phaseCapPressure,
1408  dPhaseCapPressure_dPhaseVolFrac,
1409  potential_,
1410  dPotential_dP_,
1411  dPotential_dC_ );
1412 
1413  } );
1414  }
1415 
1416 
1417  template< localIndex numComp, localIndex numFluxSupportPoints >
1418  GEOS_HOST_DEVICE
1419  static
1420  void computePotentialCapillary( localIndex const numPhase,
1421  localIndex const ip,
1422  localIndex const (&seri)[numFluxSupportPoints],
1423  localIndex const (&sesri)[numFluxSupportPoints],
1424  localIndex const (&sei)[numFluxSupportPoints],
1425  real64 const (&transmissibility)[2],
1426  real64 const (&dTrans_dPres)[2],
1427  real64 const totalFlux,
1428  ElementViewConst< arrayView1d< real64 const > > const & GEOS_UNUSED_PARAM( pres ),
1429  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1430  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &
1431  phaseMob,
1432  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1433  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1434  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1435  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1436  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1437  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1438  integer const GEOS_UNUSED_PARAM( hasCapPressure ),
1439  real64 & potential )
1440 
1441  {
1442 
1443  //Form total velocity
1444  potential = 0;
1445 
1446  //the arg lambda allows us to access some genericity
1447  UpwindScheme::template potential< numComp, numFluxSupportPoints >( numPhase, ip, seri, sesri, sei,
1448  phaseMob, potential,
1449  [&]( localIndex ipp,
1450  real64 & potential_,
1451  real64 (& dPotential_dP_)[numFluxSupportPoints],
1452  real64 (& dPotential_dC_)[numFluxSupportPoints][numComp] )
1453  {
1454 
1455  UpwindHelpers::computePotentialCapillary::compute< numComp, numFluxSupportPoints >(
1456  numPhase,
1457  ipp,
1458  seri,
1459  sesri,
1460  sei,
1461  transmissibility,
1462  dTrans_dPres,
1463  totalFlux,
1464  gravCoef,
1465  dCompFrac_dCompDens,
1466  phaseMassDens,
1467  dPhaseMassDens,
1468  dPhaseVolFrac,
1469  phaseCapPressure,
1470  dPhaseCapPressure_dPhaseVolFrac,
1471  potential_,
1472  dPotential_dP_,
1473  dPotential_dC_ );
1474  } );
1475  }
1476 
1477 };
1478 
1479 /*** IHU ***/
1480 
1481 struct IHUPhaseFlux
1482 {
1483 
1484  using UPWIND_SCHEME = HybridUpwind;
1485 
1486  static constexpr double minTotMob = 1e-12;
1487 
1515  template< integer numComp, integer numFluxSupportPoints >
1516  GEOS_HOST_DEVICE
1517  static void
1518  compute( integer const numPhase,
1519  integer const ip,
1520  integer const hasCapPressure,
1521  integer const checkPhasePresenceInGravity,
1522  localIndex const ( &seri )[numFluxSupportPoints],
1523  localIndex const ( &sesri )[numFluxSupportPoints],
1524  localIndex const ( &sei )[numFluxSupportPoints],
1525  real64 const ( &trans )[2],
1526  real64 const ( &dTrans_dPres )[2],
1527  ElementViewConst< arrayView1d< real64 const > > const & pres,
1528  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
1529  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1530  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
1531  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
1532  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
1533  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1534  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
1535  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1536  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
1537  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
1538  real64 & potGrad,
1539  real64 ( &phaseFlux ),
1540  real64 ( & dPhaseFlux_dP )[numFluxSupportPoints],
1541  real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp] )
1542  {
1543 
1544  //loop over all phases to form total velocity
1545  real64 totFlux{};
1546  real64 dTotFlux_dP[numFluxSupportPoints]{};
1547  real64 dTotFlux_dC[numFluxSupportPoints][numComp]{};
1548 
1549  //store totMob upwinded by PPU for later schemes
1550  real64 totMob{};
1551  real64 dTotMob_dP[numFluxSupportPoints]{};
1552  real64 dTotMob_dC[numFluxSupportPoints][numComp]{};
1553 
1554  //unelegant but need dummy when forming PPU total velocity
1555  real64 dPhaseFlux_dTrans;
1556 
1557  for( integer jp = 0; jp < numPhase; ++jp )
1558  {
1559  PPUPhaseFlux::compute( numPhase, jp,
1560  hasCapPressure, checkPhasePresenceInGravity,
1561  seri, sesri, sei,
1562  trans, dTrans_dPres,
1563  pres, gravCoef,
1564  phaseMob, dPhaseMob,
1565  phaseVolFrac, dPhaseVolFrac,
1566  dCompFrac_dCompDens,
1567  phaseMassDens, dPhaseMassDens,
1568  phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac,
1569  potGrad, phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, dPhaseFlux_dTrans );
1570 
1571  // accumulate into total flux
1572  UpwindHelpers::addToValueAndDerivatives( phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC,
1573  totFlux, dTotFlux_dP, dTotFlux_dC );
1574 /*
1575  // old wrong way:
1576  for( localIndex ke = 0; ke < numFluxSupportPoints; ++ke )
1577  {
1578  totMob += phaseMob[seri[ke]][sesri[ke]][sei[ke]][jp];
1579  dTotMob_dP[ke] += dPhaseMob[seri[ke]][sesri[ke]][sei[ke]][jp][Deriv::dP];
1580  for( localIndex jc = 0; jc < numComp; ++jc )
1581  {
1582  dTotMob_dC[ke][jc] += dPhaseMob[seri[ke]][sesri[ke]][sei[ke]][jp][Deriv::dC + jc];
1583  }
1584  }
1585  */
1586 
1587  // new correct way:
1588  // choose upstream cell
1589  localIndex const k_up = (phaseFlux >= 0) ? 0 : 1;
1590  // accumulate into total mobility
1591  UpwindHelpers::assignMobilityAndDerivatives( jp, k_up, seri, sesri, sei,
1592  phaseMob, dPhaseMob,
1593  totMob, dTotMob_dP[k_up], dTotMob_dC[k_up] );
1594  }
1595 
1596  // safeguard
1597  totMob = LvArray::math::max( totMob, minTotMob );
1598 
1599  // assign to zero
1600  UpwindHelpers::assignToZero( phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
1601 
1602  // fractional flow loop with IHU
1603  // maybe needed to have density out for upwinding
1604 
1605  // choose upstream cell
1606  // create local work arrays
1607  real64 viscousPhaseFlux{};
1608  real64 dViscousPhaseFlux_dP[numFluxSupportPoints]{};
1609  real64 dViscousPhaseFlux_dC[numFluxSupportPoints][numComp]{};
1610 
1611  real64 fractionalFlow{};
1612  real64 dFractionalFlow_dP[numFluxSupportPoints]{};
1613  real64 dFractionalFlow_dC[numFluxSupportPoints][numComp]{};
1614 
1615  // and the fractional flow for viscous part as \lambda_i^{up}/\sum_{NP}(\lambda_j^{up}) with up decided upon
1616  // the Upwind strategy
1617  UpwindHelpers::computeFractionalFlowViscous< numComp, numFluxSupportPoints,
1618  UPWIND_SCHEME >( numPhase,
1619  ip,
1620  seri,
1621  sesri,
1622  sei,
1623  trans,
1624  dTrans_dPres,
1625  totFlux,
1626  totMob,
1627  dTotMob_dP,
1628  dTotMob_dC,
1629  pres,
1630  gravCoef,
1631  dCompFrac_dCompDens,
1632  phaseMassDens,
1633  dPhaseMassDens,
1634  phaseMob,
1635  dPhaseMob,
1636  dPhaseVolFrac,
1637  phaseCapPressure,
1638  dPhaseCapPressure_dPhaseVolFrac,
1639  hasCapPressure,
1640  fractionalFlow,
1641  dFractionalFlow_dP,
1642  dFractionalFlow_dC );
1643 
1645  viscousPhaseFlux = fractionalFlow * totFlux;
1646 
1647  UpwindHelpers::addToDerivativesScaled( dFractionalFlow_dP, dFractionalFlow_dC,
1648  totFlux,
1649  dViscousPhaseFlux_dP, dViscousPhaseFlux_dC );
1650 
1651  //NON-FIXED UT -- to be canceled out if considered fixed
1652  UpwindHelpers::addToDerivativesScaled( dTotFlux_dP, dTotFlux_dC,
1653  fractionalFlow,
1654  dViscousPhaseFlux_dP, dViscousPhaseFlux_dC );
1655 
1656  // accumulate in the flux and its derivatives
1657  UpwindHelpers::addToValueAndDerivatives( viscousPhaseFlux, dViscousPhaseFlux_dP, dViscousPhaseFlux_dC,
1658  phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
1659 
1661 
1662  real64 gravitationalPhaseFlux{};
1663  real64 gravitationalPhaseFlux_dP[numFluxSupportPoints]{};
1664  real64 gravitationalPhaseFlux_dC[numFluxSupportPoints][numComp]{};
1665 
1666  UpwindHelpers::computePotentialFluxesGravity< numComp,
1667  numFluxSupportPoints, UPWIND_SCHEME >(
1668  numPhase,
1669  ip,
1670  seri,
1671  sesri,
1672  sei,
1673  trans,
1674  dTrans_dPres,
1675  totFlux,
1676  totMob,
1677  dTotMob_dP,
1678  dTotMob_dC,
1679  pres,
1680  gravCoef,
1681  phaseMob,
1682  dPhaseMob,
1683  phaseVolFrac,
1684  dPhaseVolFrac,
1685  dCompFrac_dCompDens,
1686  phaseMassDens,
1687  dPhaseMassDens,
1688  phaseCapPressure,
1689  dPhaseCapPressure_dPhaseVolFrac,
1690  hasCapPressure,
1691  checkPhasePresenceInGravity,
1692  gravitationalPhaseFlux,
1693  gravitationalPhaseFlux_dP,
1694  gravitationalPhaseFlux_dC );
1695 
1696  //update phaseFlux from gravitational
1697  UpwindHelpers::addToValueAndDerivatives( gravitationalPhaseFlux, gravitationalPhaseFlux_dP, gravitationalPhaseFlux_dC,
1698  phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
1699 
1700  if( hasCapPressure )
1701  {
1703 
1704  real64 capillaryPhaseFlux{};
1705  real64 capillaryPhaseFlux_dP[numFluxSupportPoints]{};
1706  real64 capillaryPhaseFlux_dC[numFluxSupportPoints][numComp]{};
1707 
1708  UpwindHelpers::computePotentialFluxesCapillary< numComp,
1709  numFluxSupportPoints, UPWIND_SCHEME >(
1710  numPhase,
1711  ip,
1712  seri,
1713  sesri,
1714  sei,
1715  trans,
1716  dTrans_dPres,
1717  totFlux,
1718  totMob,
1719  dTotMob_dP,
1720  dTotMob_dC,
1721  pres,
1722  gravCoef,
1723  phaseMob,
1724  dPhaseMob,
1725  dPhaseVolFrac,
1726  dCompFrac_dCompDens,
1727  phaseMassDens,
1728  dPhaseMassDens,
1729  phaseCapPressure,
1730  dPhaseCapPressure_dPhaseVolFrac,
1731  hasCapPressure,
1732  capillaryPhaseFlux,
1733  capillaryPhaseFlux_dP,
1734  capillaryPhaseFlux_dC );
1735 
1736  //update phaseFlux from capillary
1737  UpwindHelpers::addToValueAndDerivatives( capillaryPhaseFlux, capillaryPhaseFlux_dP, capillaryPhaseFlux_dC,
1738  phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
1739 
1740  }//end if cappres
1741  }
1742 
1743 };
1744 
1745 } // namespace isothermalCompositionalMultiPhaseFVMKernelUtilities
1746 
1747 } // namespace geos
1748 
1749 
1750 #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], 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:1304
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:1089
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:1139
geos::arrayView1d
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:188
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::integer
std::int32_t integer
Signed integer type.
Definition: DataTypes.hpp:81
geos::arrayView4d
ArrayView< T, 4, USD > arrayView4d
Alias for 4D array view.
Definition: DataTypes.hpp:236
geos::arrayView2d
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:204
geos::arrayView3d
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
Definition: DataTypes.hpp:220
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:1518
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