GEOS
HU2PhaseFlux.hpp
1 /*
2  * ------------------------------------------------------------------------------------------------------------
3  * SPDX-License-Identifier: LGPL-2.1-only
4  *
5  * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
6  * Copyright (c) 2018-2024 Total, S.A
7  * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
8  * Copyright (c) 2023-2024 Chevron
9  * Copyright (c) 2019- GEOS/GEOSX Contributors
10  * All rights reserved
11  *
12  * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
13  * ------------------------------------------------------------------------------------------------------------
14  */
15 
20 #ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_IHU2PHASEFLUX_HPP
21 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_IHU2PHASEFLUX_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 /*** HU 2 phase simplified version ***/
42 
43 struct HU2PhaseFlux
44 {
45 
46  static constexpr double minTotMob = 1e-12;
47 
75  template< integer numComp, integer numFluxSupportPoints >
76  GEOS_HOST_DEVICE
77  static void
78  compute( integer const numPhase,
79  integer const ip,
80  integer const hasCapPressure,
81  integer const checkPhasePresenceInGravity,
82  localIndex const ( &seri )[numFluxSupportPoints],
83  localIndex const ( &sesri )[numFluxSupportPoints],
84  localIndex const ( &sei )[numFluxSupportPoints],
85  real64 const ( &trans )[2],
86  real64 const ( &dTrans_dPres )[2],
87  ElementViewConst< arrayView1d< real64 const > > const & pres,
88  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
89  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
90  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
91  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
92  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
93  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
94  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
95  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
96  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
97  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
98  real64 & GEOS_UNUSED_PARAM( potGrad ),
99  real64 ( &phaseFlux ),
100  real64 ( & dPhaseFlux_dP )[numFluxSupportPoints],
101  real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp] )
102  {
103  // viscous part
104  computeViscousFlux< numComp, numFluxSupportPoints >( ip, numPhase,
105  hasCapPressure,
106  checkPhasePresenceInGravity,
107  seri, sesri, sei,
108  trans, dTrans_dPres,
109  pres, gravCoef,
110  dCompFrac_dCompDens,
111  phaseMassDens, dPhaseMassDens,
112  phaseMob, dPhaseMob,
113  phaseVolFrac, dPhaseVolFrac,
114  phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac,
115  phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
116  // gravity part
117  computeGravityFlux< numComp, numFluxSupportPoints >( ip, numPhase,
118  checkPhasePresenceInGravity,
119  seri, sesri, sei,
120  trans, dTrans_dPres, gravCoef,
121  phaseVolFrac,
122  phaseMob, dPhaseMob,
123  dCompFrac_dCompDens,
124  phaseMassDens, dPhaseMassDens,
125  phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
126 
127  // capillary part
128  if( hasCapPressure )
129  {
130  computeCapillaryFlux< numComp, numFluxSupportPoints >( ip, numPhase,
131  seri, sesri, sei,
132  trans, dTrans_dPres,
133  phaseMob, dPhaseMob,
134  dPhaseVolFrac,
135  phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac,
136  phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
137  }
138  }
139 
140 protected:
141 
142  template< localIndex numComp, localIndex numFluxSupportPoints >
143  GEOS_HOST_DEVICE
144  static void
145  computeViscousFlux( integer const & ip,
146  integer const & numPhase,
147  integer const & hasCapPressure,
148  integer const & checkPhasePresenceInGravity,
149  localIndex const (&seri)[numFluxSupportPoints],
150  localIndex const (&sesri)[numFluxSupportPoints],
151  localIndex const (&sei)[numFluxSupportPoints],
152  real64 const (&trans)[2], real64 const (&dTrans_dPres)[2],
153  ElementViewConst< arrayView1d< real64 const > > const & pres,
154  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
155  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
156  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
157  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
158  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
159  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
160  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
161  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
162  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
163  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
164  real64 ( &phaseFlux ),
165  real64 ( & dPhaseFlux_dP )[numFluxSupportPoints],
166  real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp] )
167  {
168  // form total velocity and derivatives (TODO move it OUT!)
169  real64 totFlux{};
170  real64 dTotFlux_dP[numFluxSupportPoints]{};
171  real64 dTotFlux_dC[numFluxSupportPoints][numComp]{};
172 
173  computeTotalFlux( numPhase,
174  hasCapPressure,
175  checkPhasePresenceInGravity,
176  seri, sesri, sei,
177  trans, dTrans_dPres,
178  pres, gravCoef,
179  phaseMob, dPhaseMob,
180  phaseVolFrac, dPhaseVolFrac,
181  dCompFrac_dCompDens,
182  phaseMassDens, dPhaseMassDens,
183  phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac,
184  totFlux, dTotFlux_dP, dTotFlux_dC );
185 
186  localIndex k_up = -1;
187  real64 mob{};
188  real64 dMob_dP{};
189  real64 dMob_dC[numComp]{};
190 
191  real64 totMob{};
192  real64 dTotMob_dP{};
193  real64 dTotMob_dC[numComp]{};
194 
195  for( localIndex jp = 0; jp < numPhase; ++jp )
196  {
197  if( jp == ip ) // ip will be computed later
198  continue;
199 
200  // upwind based on totFlux sign
201  upwindMobility< numComp, numFluxSupportPoints >( jp,
202  seri,
203  sesri,
204  sei,
205  totFlux,
206  phaseMob,
207  dPhaseMob,
208  k_up,
209  mob,
210  dMob_dP,
211  dMob_dC );
212  // accumulate total mobility
213  UpwindHelpers::addToValueAndDerivatives( mob, dMob_dP, dMob_dC, totMob, dTotMob_dP, dTotMob_dC );
214  }
215 
216  // upwind based on totFlux sign
217  upwindMobility< numComp, numFluxSupportPoints >( ip,
218  seri,
219  sesri,
220  sei,
221  totFlux,
222  phaseMob,
223  dPhaseMob,
224  k_up,
225  mob,
226  dMob_dP,
227  dMob_dC );
228  // accumulate total mobility
229  UpwindHelpers::addToValueAndDerivatives( mob, dMob_dP, dMob_dC, totMob, dTotMob_dP, dTotMob_dC );
230 
231  // safeguard
232  totMob = LvArray::math::max( totMob, minTotMob );
233  real64 const invTotMob = 1 / totMob;
234 
235  // fractional flow for viscous part as \lambda_i^{up}/\sum_{NP}(\lambda_j^{up})
236 
237  real64 const fractionalFlow = mob * invTotMob;
238  real64 dFractionalFlow_dP{};
239  real64 dFractionalFlow_dC[numComp]{};
240  UpwindHelpers::addToDerivativesScaled( dMob_dP, dMob_dC, invTotMob, dFractionalFlow_dP, dFractionalFlow_dC );
241  UpwindHelpers::addToDerivativesScaled( dTotMob_dP, dTotMob_dC, -fractionalFlow * invTotMob, dFractionalFlow_dP, dFractionalFlow_dC );
242 
244 
245  real64 const viscousPhaseFlux = fractionalFlow * totFlux;
246  real64 dViscousPhaseFlux_dP[numFluxSupportPoints]{};
247  real64 dViscousPhaseFlux_dC[numFluxSupportPoints][numComp]{};
248 
249  // fractionalFlow derivatives
250  UpwindHelpers::addToDerivativesScaled( dFractionalFlow_dP, dFractionalFlow_dC, totFlux, dViscousPhaseFlux_dP[k_up], dViscousPhaseFlux_dC[k_up] );
251 
252  // Ut derivatives
253  UpwindHelpers::addToDerivativesScaled( dTotFlux_dP, dTotFlux_dC, fractionalFlow, dViscousPhaseFlux_dP, dViscousPhaseFlux_dC );
254 
255  // accumulate in the flux and its derivatives (need to be very careful doing that)
256  UpwindHelpers::addToValueAndDerivatives( viscousPhaseFlux, dViscousPhaseFlux_dP, dViscousPhaseFlux_dC,
257  phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
258  }
259 
260  template< localIndex numComp, localIndex numFluxSupportPoints >
261  GEOS_HOST_DEVICE
262  static void
263  computeGravityFlux( integer const & ip, integer const & numPhase,
264  integer const checkPhasePresenceInGravity,
265  localIndex const (&seri)[numFluxSupportPoints],
266  localIndex const (&sesri)[numFluxSupportPoints],
267  localIndex const (&sei)[numFluxSupportPoints],
268  real64 const (&trans)[2], real64 const (&dTrans_dPres)[2],
269  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
270  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
271  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
272  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
273  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
274  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
275  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
276  real64 & phaseFlux,
277  real64 (& dPhaseFlux_dP)[numFluxSupportPoints],
278  real64 (& dPhaseFlux_dC)[numFluxSupportPoints][numComp] )
279  {
281  real64 gravPhaseFlux{};
282  real64 dGravPhaseFlux_dP[numFluxSupportPoints]{};
283  real64 dGravPhaseFlux_dC[numFluxSupportPoints][numComp]{};
284 
285  real64 pot_i{};
286  real64 dPot_i_dP[numFluxSupportPoints]{};
287  real64 dPot_i_dC[numFluxSupportPoints][numComp]{};
288  computeGravityPotential< numComp, numFluxSupportPoints >( ip,
289  seri,
290  sesri,
291  sei,
292  checkPhasePresenceInGravity,
293  trans,
294  dTrans_dPres,
295  gravCoef,
296  dCompFrac_dCompDens,
297  phaseVolFrac,
298  phaseMassDens,
299  dPhaseMassDens,
300  pot_i,
301  dPot_i_dP,
302  dPot_i_dC );
303 
304  for( localIndex jp = 0; jp < numPhase; ++jp )
305  {
306  if( ip != jp )
307  {
308  real64 pot_j{};
309  real64 dPot_j_dP[numFluxSupportPoints]{};
310  real64 dPot_j_dC[numFluxSupportPoints][numComp]{};
311  computeGravityPotential< numComp, numFluxSupportPoints >( jp,
312  seri,
313  sesri,
314  sei,
315  checkPhasePresenceInGravity,
316  trans,
317  dTrans_dPres,
318  gravCoef,
319  dCompFrac_dCompDens,
320  phaseVolFrac,
321  phaseMassDens,
322  dPhaseMassDens,
323  pot_j,
324  dPot_j_dP,
325  dPot_j_dC );
326 
327  computePotDiffFlux( ip, jp,
328  seri, sesri, sei,
329  pot_i, dPot_i_dP, dPot_i_dC,
330  pot_j, dPot_j_dP, dPot_j_dC,
331  phaseMob, dPhaseMob,
332  gravPhaseFlux, dGravPhaseFlux_dP, dGravPhaseFlux_dC );
333  }
334  }
335 
336  // update phaseFlux from gravitational
337  UpwindHelpers::addToValueAndDerivatives( gravPhaseFlux, dGravPhaseFlux_dP, dGravPhaseFlux_dC,
338  phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
339  }
340 
341  template< localIndex numComp, localIndex numFluxSupportPoints >
342  GEOS_HOST_DEVICE
343  static void
344  computeCapillaryFlux( integer const & ip, integer const & numPhase,
345  localIndex const (&seri)[numFluxSupportPoints],
346  localIndex const (&sesri)[numFluxSupportPoints],
347  localIndex const (&sei)[numFluxSupportPoints],
348  real64 const (&trans)[2], real64 const (&dTrans_dPres)[2],
349  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
350  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
351  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
352  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
353  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
354  real64 & phaseFlux,
355  real64 (& dPhaseFlux_dP)[numFluxSupportPoints],
356  real64 (& dPhaseFlux_dC)[numFluxSupportPoints][numComp] )
357  {
359  real64 capPhaseFlux{};
360  real64 dCapPhaseFlux_dP[numFluxSupportPoints]{};
361  real64 dCapPhaseFlux_dC[numFluxSupportPoints][numComp]{};
362 
363  real64 pot_i{};
364  real64 dPot_i_dP[numFluxSupportPoints]{};
365  real64 dPot_i_dC[numFluxSupportPoints][numComp]{};
366 
367  computeCapillaryPotential< numComp, numFluxSupportPoints >( ip,
368  numPhase,
369  seri,
370  sesri,
371  sei,
372  trans,
373  dTrans_dPres,
374  dPhaseVolFrac,
375  phaseCapPressure,
376  dPhaseCapPressure_dPhaseVolFrac,
377  pot_i,
378  dPot_i_dP,
379  dPot_i_dC );
380 
381  for( localIndex jp = 0; jp < numPhase; ++jp )
382  {
383  if( ip != jp )
384  {
385  real64 pot_j{};
386  real64 dPot_j_dP[numFluxSupportPoints]{};
387  real64 dPot_j_dC[numFluxSupportPoints][numComp]{};
388  computeCapillaryPotential< numComp, numFluxSupportPoints >( jp,
389  numPhase,
390  seri,
391  sesri,
392  sei,
393  trans,
394  dTrans_dPres,
395  dPhaseVolFrac,
396  phaseCapPressure,
397  dPhaseCapPressure_dPhaseVolFrac,
398  pot_j,
399  dPot_j_dP,
400  dPot_j_dC );
401 
402  computePotDiffFlux( ip, jp,
403  seri, sesri, sei,
404  pot_i, dPot_i_dP, dPot_i_dC,
405  pot_j, dPot_j_dP, dPot_j_dC,
406  phaseMob, dPhaseMob,
407  capPhaseFlux, dCapPhaseFlux_dP, dCapPhaseFlux_dC );
408  }
409  }
410 
411  // update phaseFlux from capillary flux
412  UpwindHelpers::addToValueAndDerivatives( capPhaseFlux, dCapPhaseFlux_dP, dCapPhaseFlux_dC,
413  phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC );
414  }
415 
416  template< localIndex numComp, localIndex numFluxSupportPoints >
417  GEOS_HOST_DEVICE
418  static void
419  computeTotalFlux( integer const & numPhase,
420  const integer & hasCapPressure,
421  const integer & checkPhasePresenceInGravity,
422  localIndex const (&seri)[numFluxSupportPoints],
423  localIndex const (&sesri)[numFluxSupportPoints],
424  localIndex const (&sei)[numFluxSupportPoints],
425  real64 const (&trans)[2], real64 const (&dTrans_dPres)[2],
426  ElementViewConst< arrayView1d< real64 const > > const & pres,
427  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
428  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
429  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
430  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
431  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
432  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
433  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
434  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
435  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
436  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
437  real64 & totFlux, real64 (& dTotFlux_dP)[numFluxSupportPoints], real64 (& dTotFlux_dC)[numFluxSupportPoints][numComp] )
438  {
439  for( integer jp = 0; jp < numPhase; ++jp )
440  {
441  // working arrays for phase flux
442  real64 potGrad{};
443  real64 phaseFlux{};
444  real64 dPhaseFlux_dP[numFluxSupportPoints]{};
445  real64 dPhaseFlux_dC[numFluxSupportPoints][numComp]{};
446  real64 dPhaseFlux_dTrans{};
447  PPUPhaseFlux::compute( numPhase, jp,
448  hasCapPressure,
449  checkPhasePresenceInGravity,
450  seri, sesri, sei,
451  trans, dTrans_dPres,
452  pres, gravCoef,
453  phaseMob, dPhaseMob,
454  phaseVolFrac, dPhaseVolFrac,
455  dCompFrac_dCompDens,
456  phaseMassDens, dPhaseMassDens,
457  phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac,
458  potGrad, phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, dPhaseFlux_dTrans );
459 
460  UpwindHelpers::addToValueAndDerivatives( phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC,
461  totFlux, dTotFlux_dP, dTotFlux_dC );
462  }
463  }
464 
465  template< localIndex numComp, localIndex numFluxSupportPoints >
466  GEOS_HOST_DEVICE
467  static void
468  upwindMobility( localIndex const ip,
469  localIndex const (&seri)[numFluxSupportPoints],
470  localIndex const (&sesri)[numFluxSupportPoints],
471  localIndex const (&sei)[numFluxSupportPoints],
472  real64 const pot,
473  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
474  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
475  localIndex & upwindDir,
476  real64 & mobility,
477  real64 & dMobility_dP,
478  real64 ( & dMobility_dC)[numComp] )
479  {
480  upwindDir = (pot > 0) ? 0 : 1;
481  UpwindHelpers::assignToZero( mobility, dMobility_dP, dMobility_dC );
482  UpwindHelpers::assignMobilityAndDerivatives( ip, upwindDir, seri, sesri, sei, phaseMob, dPhaseMob, mobility, dMobility_dP, dMobility_dC );
483  }
484 
485  template< localIndex numComp, localIndex numFluxSupportPoints >
486  GEOS_HOST_DEVICE
487  static void computeGravityPotential( localIndex const ip,
488  localIndex const (&seri)[numFluxSupportPoints],
489  localIndex const (&sesri)[numFluxSupportPoints],
490  localIndex const (&sei)[numFluxSupportPoints],
491  integer const checkPhasePresenceInGravity,
492  real64 const (&trans)[2],
493  real64 const (&dTrans_dP)[2],
494  ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
495  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
496  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
497  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
498  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
499  real64 & gravPot,
500  real64 ( & dGravPot_dP )[numFluxSupportPoints],
501  real64 ( & dGravPot_dC )[numFluxSupportPoints][numComp] )
502  {
503  // init
504  UpwindHelpers::assignToZero( gravPot, dGravPot_dP, dGravPot_dC );
505 
506  // get average density
507  real64 densMean{};
508  real64 dDensMean_dP[numFluxSupportPoints]{};
509  real64 dDensMean_dC[numFluxSupportPoints][numComp]{};
510  isothermalCompositionalMultiphaseFVMKernels::helpers::
511  calculateMeanDensity( ip, seri, sesri, sei,
512  checkPhasePresenceInGravity,
513  phaseVolFrac, dCompFrac_dCompDens,
514  phaseMassDens, dPhaseMassDens,
515  densMean, dDensMean_dP, dDensMean_dC );
516 
517  // compute potential difference MPFA-style
518  for( localIndex i = 0; i < numFluxSupportPoints; ++i )
519  {
520  localIndex const er = seri[i];
521  localIndex const esr = sesri[i];
522  localIndex const ei = sei[i];
523 
524  real64 const gravD = trans[i] * gravCoef[er][esr][ei];
525  real64 const dGravD_dP = dTrans_dP[i] * gravCoef[er][esr][ei];
526  gravPot += densMean * gravD;
527  dGravPot_dP[i] += densMean * dGravD_dP;
528 
529  // need to add contributions from both cells the mean density depends on
530  UpwindHelpers::addToDerivativesScaled( dDensMean_dP, dDensMean_dC, gravD, dGravPot_dP, dGravPot_dC );
531  }
532 
533  }
534 
535  template< localIndex numComp, localIndex numFluxSupportPoints >
536  GEOS_HOST_DEVICE
537  static void computeCapillaryPotential( localIndex const ip,
538  localIndex const numPhase,
539  localIndex const (&seri)[numFluxSupportPoints],
540  localIndex const (&sesri)[numFluxSupportPoints],
541  localIndex const (&sei)[numFluxSupportPoints],
542  real64 const (&trans)[2],
543  real64 const (&dTrans_dPres)[2],
544  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
545  ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
546  ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
547  real64 & capPot,
548  real64 ( & dCapPot_dP )[numFluxSupportPoints],
549  real64 ( & dCapPot_dC )[numFluxSupportPoints][numComp] )
550  {
551  // init
552  UpwindHelpers::assignToZero( capPot, dCapPot_dP, dCapPot_dC );
553 
554  // need to add contributions from both cells
555  for( localIndex i = 0; i < numFluxSupportPoints; ++i )
556  {
557  localIndex const er = seri[i];
558  localIndex const esr = sesri[i];
559  localIndex const ei = sei[i];
560 
561  real64 const capPressure = phaseCapPressure[er][esr][ei][0][ip];
562  real64 dCapPressure_dP{};
563  real64 dCapPressure_dC[numComp]{};
564  for( integer jp = 0; jp < numPhase; ++jp )
565  {
566  real64 const dCapPressure_dS = dPhaseCapPressure_dPhaseVolFrac[er][esr][ei][0][ip][jp];
567  dCapPressure_dP += dCapPressure_dS * dPhaseVolFrac[er][esr][ei][jp][Deriv::dP];
568 
569  for( integer jc = 0; jc < numComp; ++jc )
570  {
571  dCapPressure_dC[jc] += dCapPressure_dS * dPhaseVolFrac[er][esr][ei][jp][Deriv::dC+jc];
572  }
573  }
574 
575  capPot += trans[i] * capPressure;
576  dCapPot_dP[i] += trans[i] * dCapPressure_dP + dTrans_dPres[i] * capPressure;
577  for( integer jc = 0; jc < numComp; ++jc )
578  {
579  dCapPot_dC[i][jc] += trans[i] * dCapPressure_dC[jc];
580  }
581  }
582  }
583 
584  template< localIndex numComp, localIndex numFluxSupportPoints >
585  GEOS_HOST_DEVICE
586  static void
587  computePotDiffFlux( integer const & ip, integer const & jp,
588  localIndex const (&seri)[numFluxSupportPoints],
589  localIndex const (&sesri)[numFluxSupportPoints],
590  localIndex const (&sei)[numFluxSupportPoints],
591  real64 const & pot_i, real64 const ( &dPot_i_dP )[numFluxSupportPoints], real64 const (&dPot_i_dC )[numFluxSupportPoints][numComp],
592  real64 const & pot_j, real64 const ( &dPot_j_dP )[numFluxSupportPoints], real64 const ( &dPot_j_dC)[numFluxSupportPoints][numComp],
593  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
594  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
595  real64 & phaseFlux, real64 ( & dPhaseFlux_dP )[numFluxSupportPoints], real64 ( & dPhaseFlux_dC)[numFluxSupportPoints][numComp] )
596  {
597  // upwind based on pot diff sign
598  real64 const potDiff = pot_j - pot_i;
599  real64 dPotDiff_dP[numFluxSupportPoints]{};
600  real64 dPotDiff_dC[numFluxSupportPoints][numComp]{};
601  UpwindHelpers::addToDerivativesScaled( dPot_j_dP, dPot_j_dC, 1.0, dPotDiff_dP, dPotDiff_dC );
602  UpwindHelpers::addToDerivativesScaled( dPot_i_dP, dPot_i_dC, -1.0, dPotDiff_dP, dPotDiff_dC );
603 
604  localIndex k_up_i = -1;
605  real64 mob_i{};
606  real64 dMob_i_dP{};
607  real64 dMob_i_dC[numComp]{};
608  upwindMobility< numComp, numFluxSupportPoints >( ip,
609  seri,
610  sesri,
611  sei,
612  potDiff,
613  phaseMob,
614  dPhaseMob,
615  k_up_i,
616  mob_i,
617  dMob_i_dP,
618  dMob_i_dC );
619  localIndex k_up_j = -1;
620  real64 mob_j{};
621  real64 dMob_j_dP{};
622  real64 dMob_j_dC[numComp]{};
623  upwindMobility< numComp, numFluxSupportPoints >( jp,
624  seri,
625  sesri,
626  sei,
627  -potDiff,
628  phaseMob,
629  dPhaseMob,
630  k_up_j,
631  mob_j,
632  dMob_j_dP,
633  dMob_j_dC );
634 
635  // safeguard
636  real64 const mobTot = LvArray::math::max( mob_i + mob_j, minTotMob );
637  real64 const mobTotInv = 1 / mobTot;
638  real64 dMobTot_dP[numFluxSupportPoints]{};
639  real64 dMobTot_dC[numFluxSupportPoints][numComp]{};
640  UpwindHelpers::addToDerivatives( dMob_i_dP, dMob_i_dC, dMobTot_dP[k_up_i], dMobTot_dC[k_up_i] );
641  UpwindHelpers::addToDerivatives( dMob_j_dP, dMob_j_dC, dMobTot_dP[k_up_j], dMobTot_dC[k_up_j] );
642 
643  // Assembling flux phase-wise as \phi_{i,g} = \sum_{k\nei} \lambda_k^{up,g} f_k^{up,g} (Pot_i - Pot_k)
644  phaseFlux += mob_i * mob_j * mobTotInv * potDiff;
645 
646  // mob_i derivatives
647  UpwindHelpers::addToDerivativesScaled( dMob_i_dP, dMob_i_dC, mob_j * mobTotInv * potDiff, dPhaseFlux_dP[k_up_i], dPhaseFlux_dC[k_up_i] );
648 
649  // mob_j derivatives
650  UpwindHelpers::addToDerivativesScaled( dMob_j_dP, dMob_j_dC, mob_i * mobTotInv * potDiff, dPhaseFlux_dP[k_up_j], dPhaseFlux_dC[k_up_j] );
651 
652  // mobTot derivatives
653  real64 const mobTotInv2 = mobTotInv * mobTotInv;
654  UpwindHelpers::addToDerivativesScaled( dMobTot_dP, dMobTot_dC, -mob_i * mob_j * mobTotInv2 * potDiff, dPhaseFlux_dP, dPhaseFlux_dC );
655 
656  // potDiff derivatives
657  UpwindHelpers::addToDerivativesScaled( dPotDiff_dP, dPotDiff_dC, mob_i * mob_j * mobTotInv, dPhaseFlux_dP, dPhaseFlux_dC );
658  }
659 
660 };
661 
662 } // namespace isothermalCompositionalMultiPhaseFVMKernelUtilities
663 
664 } // namespace geos
665 
666 #endif // GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_IHU2PHASEFLUX_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::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::HU2PhaseFlux::computeGravityFlux
static GEOS_HOST_DEVICE void computeGravityFlux(integer const &ip, integer const &numPhase, 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 &gravCoef, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseVolFrac, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseMob, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &dPhaseMob, 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, real64 &phaseFlux, real64(&dPhaseFlux_dP)[numFluxSupportPoints], real64(&dPhaseFlux_dC)[numFluxSupportPoints][numComp])
Definition: HU2PhaseFlux.hpp:263
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::HU2PhaseFlux::computeCapillaryFlux
static GEOS_HOST_DEVICE void computeCapillaryFlux(integer const &ip, integer const &numPhase, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], real64 const (&trans)[2], real64 const (&dTrans_dPres)[2], ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const &phaseMob, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const &dPhaseMob, 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 &phaseFlux, real64(&dPhaseFlux_dP)[numFluxSupportPoints], real64(&dPhaseFlux_dC)[numFluxSupportPoints][numComp])
Definition: HU2PhaseFlux.hpp:344
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::HU2PhaseFlux::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 &GEOS_UNUSED_PARAM(potGrad), real64(&phaseFlux), real64(&dPhaseFlux_dP)[numFluxSupportPoints], real64(&dPhaseFlux_dC)[numFluxSupportPoints][numComp])
Simplified 2-phase version of hybrid upwinding.
Definition: HU2PhaseFlux.hpp:78
geos::isothermalCompositionalMultiphaseFVMKernelUtilities::HU2PhaseFlux::computeViscousFlux
static GEOS_HOST_DEVICE void computeViscousFlux(integer const &ip, integer const &numPhase, 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< 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 &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, constitutive::cappres::USD_CAPPRES > > const &phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const &dPhaseCapPressure_dPhaseVolFrac, real64(&phaseFlux), real64(&dPhaseFlux_dP)[numFluxSupportPoints], real64(&dPhaseFlux_dC)[numFluxSupportPoints][numComp])
Definition: HU2PhaseFlux.hpp:145
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
Generated by doxygen 1.9.1