GEOS
CompositionalMultiphaseHybridFVMKernels_impl.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 #include "CompositionalMultiphaseHybridFVMKernels.hpp"
21 
22 #include "finiteVolume/mimeticInnerProducts/MimeticInnerProductBase.hpp"
23 #include "finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp"
24 #include "finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp"
25 #include "finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp"
26 
27 #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp"
28 #include "physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp"
29 #include "physicsSolvers/fluidFlow/kernels/HybridFVMHelperKernels.hpp"
30 #include "constitutive/fluid/multifluid/MultiFluidSelector.hpp"
31 #include "constitutive/relativePermeability/RelativePermeabilitySelector.hpp"
32 
33 namespace geos
34 {
35 using namespace constitutive;
36 
37 namespace compositionalMultiphaseHybridFVMKernels
38 {
39 
40 /******************************** UpwindingHelper ********************************/
41 
42 template< integer NC, integer NP >
43 GEOS_HOST_DEVICE
44 inline
45 void
46 UpwindingHelper::
47  upwindViscousCoefficient( localIndex const (&localIds)[ 3 ],
48  localIndex const (&neighborIds)[ 3 ],
49  ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseDens,
50  ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens,
51  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
52  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
53  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
54  ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_COMP > > const & phaseCompFrac,
55  ElementViewConst< arrayView5d< real64 const, multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac,
56  ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber,
57  real64 const & oneSidedVolFlux,
58  real64 ( & upwPhaseViscCoef )[ NP ][ NC ],
59  real64 ( & dUpwPhaseViscCoef_dPres )[ NP ][ NC ],
60  real64 ( & dUpwPhaseViscCoef_dCompDens )[ NP ][ NC ][ NC ],
61  globalIndex & upwViscDofNumber )
62 {
63  using Deriv = multifluid::DerivativeOffset;
64 
65  real64 dUpwMobRatio_dCompDens[ NC ]{};
66  real64 dUpwDensMobRatio_dCompDens[ NC ]{};
67  real64 dPhaseDens_dC[ NC ]{};
68  real64 dPhaseCompFrac_dC[ NC ]{};
69 
70  // 1) Upwind
71  localIndex const er = ( oneSidedVolFlux > 0 ) ? localIds[0] : neighborIds[0];
72  localIndex const esr = ( oneSidedVolFlux > 0 ) ? localIds[1] : neighborIds[1];
73  localIndex const ei = ( oneSidedVolFlux > 0 ) ? localIds[2] : neighborIds[2];
74 
75  // 2) Compute total mobility: \lambda_T = \sum_{\ell} \lambda_{\ell}
76  real64 totalMob = 0;
77  real64 dTotalMob_dPres = 0;
78  real64 dTotalMob_dCompDens[ NC ]{};
79  for( integer ip = 0; ip < NP; ++ip )
80  {
81  totalMob = totalMob + phaseMob[er][esr][ei][ip];
82  dTotalMob_dPres = dTotalMob_dPres + dPhaseMob[er][esr][ei][ip][Deriv::dP];
83  for( integer ic = 0; ic < NC; ++ic )
84  {
85  dTotalMob_dCompDens[ic] = dTotalMob_dCompDens[ic] + dPhaseMob[er][esr][ei][ip][Deriv::dC+ic];
86  }
87  }
88  real64 const totalMobInv = 1.0 / totalMob;
89 
90  for( integer ip = 0; ip < NP; ++ip )
91  {
92  // 3) Compute viscous mobility ratio: \frac{\lambda_{\ell}}{\lambda_T}
93  real64 const upwMobRatio = phaseMob[er][esr][ei][ip] * totalMobInv;
94  real64 const dUpwMobRatio_dPres = ( dPhaseMob[er][esr][ei][ip][Deriv::dP] - upwMobRatio * dTotalMob_dPres )
95  * totalMobInv;
96  for( integer ic = 0; ic < NC; ++ic )
97  {
98  dUpwMobRatio_dCompDens[ic] = ( dPhaseMob[er][esr][ei][ip][Deriv::dC+ic] - upwMobRatio * dTotalMob_dCompDens[ic] )
99  * totalMobInv;
100  }
101 
102  // 4) Multiply mobility ratio by phase density: \rho^{up}_{\ell} \frac{\lambda_{\ell}}{\lambda_T}
103  applyChainRule( NC,
104  dCompFrac_dCompDens[er][esr][ei],
105  dPhaseDens[er][esr][ei][0][ip],
106  dPhaseDens_dC,
107  Deriv::dC );
108  real64 const upwDensMobRatio = phaseDens[er][esr][ei][0][ip] * upwMobRatio;
109  real64 const dUpwDensMobRatio_dPres = dPhaseDens[er][esr][ei][0][ip][Deriv::dP] * upwMobRatio
110  + phaseDens[er][esr][ei][0][ip] * dUpwMobRatio_dPres;
111  for( integer ic = 0; ic < NC; ++ic )
112  {
113  dUpwDensMobRatio_dCompDens[ic] = dPhaseDens_dC[ic] * upwMobRatio
114  + phaseDens[er][esr][ei][0][ip] * dUpwMobRatio_dCompDens[ic];
115  }
116 
117  // 5) Multiply density mobility ratio by phase comp fraction: x_{c,\ell} \rho^{up}_{\ell} \frac{\lambda_{\ell}}{\lambda_T}
118  for( integer ic = 0; ic < NC; ++ic )
119  {
120  applyChainRule( NC,
121  dCompFrac_dCompDens[er][esr][ei],
122  dPhaseCompFrac[er][esr][ei][0][ip][ic],
123  dPhaseCompFrac_dC,
124  Deriv::dC );
125  upwPhaseViscCoef[ip][ic] = phaseCompFrac[er][esr][ei][0][ip][ic] * upwDensMobRatio;
126  dUpwPhaseViscCoef_dPres[ip][ic] = dPhaseCompFrac[er][esr][ei][0][ip][ic][Deriv::dP] * upwDensMobRatio
127  + phaseCompFrac[er][esr][ei][0][ip][ic] * dUpwDensMobRatio_dPres;
128  for( integer jc = 0; jc < NC; ++jc )
129  {
130  dUpwPhaseViscCoef_dCompDens[ip][ic][jc] = dPhaseCompFrac_dC[jc] * upwDensMobRatio
131  + phaseCompFrac[er][esr][ei][0][ip][ic] * dUpwDensMobRatio_dCompDens[jc];
132  }
133  }
134  }
135  // 6) Save the dof number of the upwind cell
136  upwViscDofNumber = elemDofNumber[er][esr][ei];
137 
138 }
139 
140 template< integer NC, integer NP >
141 GEOS_HOST_DEVICE
142 inline
143 void
144 UpwindingHelper::
145  upwindBuoyancyCoefficient( localIndex const (&localIds)[ 3 ],
146  localIndex const (&neighborIds)[ 3 ],
147  real64 const & transGravCoef,
148  ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseDens,
149  ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens,
150  ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseMassDens,
151  ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
152  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
153  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
154  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
155  ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_COMP > > const & phaseCompFrac,
156  ElementViewConst< arrayView5d< real64 const, multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac,
157  real64 ( & phaseGravTerm )[ NP ][ NP-1 ],
158  real64 ( & dPhaseGravTerm_dPres )[ NP ][ NP-1 ][ 2 ],
159  real64 ( & dPhaseGravTerm_dCompDens )[ NP ][ NP-1 ][ 2 ][ NC ],
160  real64 ( & upwPhaseGravCoef )[ NP ][ NP-1 ][ NC ],
161  real64 ( & dUpwPhaseGravCoef_dPres )[ NP ][ NP-1 ][ NC ][ 2 ],
162  real64 ( & dUpwPhaseGravCoef_dCompDens )[ NP ][ NP-1 ][ NC ][ 2 ][ NC ] )
163 {
164  using Deriv = multifluid::DerivativeOffset;
165 
166  // 1) Compute the driving force: T ( \rho^{avg}_{\ell} - \rho^{avg}_m ) g \Delta z
167  computePhaseGravTerm( localIds,
168  neighborIds,
169  transGravCoef,
170  phaseMassDens,
171  dPhaseMassDens,
172  dCompFrac_dCompDens,
173  phaseGravTerm,
174  dPhaseGravTerm_dPres,
175  dPhaseGravTerm_dCompDens );
176 
177  // 2) Compute the total mobility: \lambda_T = \sum_{\ell} \lambda_{\ell}
178  real64 totalMob = 0.0;
179  real64 dTotalMob_dPres[ 2 ]{};
180  real64 dTotalMob_dCompDens[ 2 ][ NC ]{};
181 
182  computeUpwindedTotalMobility( localIds,
183  neighborIds,
184  phaseMob,
185  dPhaseMob,
186  phaseGravTerm,
187  totalMob,
188  dTotalMob_dPres,
189  dTotalMob_dCompDens );
190  real64 const totalMobInv = 1.0 / totalMob;
191 
192  // 3) Compute the quantities \x_{up}_{c,p} \rho_p \frac{\lambda_p \lambda_m}{\lambda_T}
193  real64 dMobRatio_dPres[ 2 ]{};
194  real64 dMobRatio_dCompDens[ 2 ][ NC ]{};
195  real64 dDensMobRatio_dPres[ 2 ]{};
196  real64 dDensMobRatio_dCompDens[ 2 ][ NC ]{};
197 
198  real64 dPhaseDens_dC[ NC ]{};
199  real64 dPhaseCompFrac_dC[ NC ]{};
200 
201  for( integer ip = 0; ip < NP; ++ip )
202  {
203  localIndex k = 0;
204  for( integer jp = 0; jp < NP; ++jp )
205  {
206  if( ip == jp )
207  {
208  continue;
209  }
210 
211  // 3.a) Upwinding using the gravity term
212  localIndex eru, esru, eiu, posu; // upwind
213  localIndex erd, esrd, eid, posd; // downwind
214  setIndicesForMobilityRatioUpwinding( localIds, neighborIds,
215  phaseGravTerm[ip][k],
216  eru, esru, eiu, posu,
217  erd, esrd, eid, posd );
218 
219  // 3.b) Compute mobility ratio \frac{\lambda_l \lambda_m}{\lambda_T}
220  real64 const mobRatio = phaseMob[eru][esru][eiu][ip] * phaseMob[erd][esrd][eid][jp]
221  * totalMobInv;
222  dMobRatio_dPres[posu] = ( dPhaseMob[eru][esru][eiu][ip][Deriv::dP] * phaseMob[erd][esrd][eid][jp]
223  - mobRatio * dTotalMob_dPres[posu] ) * totalMobInv;
224  dMobRatio_dPres[posd] = ( dPhaseMob[erd][esrd][eid][jp][Deriv::dP] * phaseMob[eru][esru][eiu][ip]
225  - mobRatio * dTotalMob_dPres[posd] ) * totalMobInv;
226 
227  for( integer ic = 0; ic < NC; ++ic )
228  {
229  dMobRatio_dCompDens[posu][ic] = ( dPhaseMob[eru][esru][eiu][ip][Deriv::dC+ic] * phaseMob[erd][esrd][eid][jp]
230  - mobRatio * dTotalMob_dCompDens[posu][ic] ) * totalMobInv;
231  dMobRatio_dCompDens[posd][ic] = ( dPhaseMob[erd][esrd][eid][jp][Deriv::dC+ic] * phaseMob[eru][esru][eiu][ip]
232  - mobRatio * dTotalMob_dCompDens[posd][ic] ) * totalMobInv;
233  }
234 
235  // 3.c) Compute mobility ratio multiplied by upwinded phase density \rho_l \frac{\lambda_l \lambda_m}{\lambda_T}
236  applyChainRule( NC,
237  dCompFrac_dCompDens[eru][esru][eiu],
238  dPhaseDens[eru][esru][eiu][0][ip],
239  dPhaseDens_dC,
240  Deriv::dC );
241  real64 const densMobRatio = phaseDens[eru][esru][eiu][0][ip] * mobRatio;
242  dDensMobRatio_dPres[posu] = dPhaseDens[eru][esru][eiu][0][ip][Deriv::dP] * mobRatio
243  + phaseDens[eru][esru][eiu][0][ip] * dMobRatio_dPres[posu];
244  dDensMobRatio_dPres[posd] = phaseDens[eru][esru][eiu][0][ip] * dMobRatio_dPres[posd];
245  for( integer ic = 0; ic < NC; ++ic )
246  {
247  dDensMobRatio_dCompDens[posu][ic] = dPhaseDens_dC[ic] * mobRatio
248  + phaseDens[eru][esru][eiu][0][ip] * dMobRatio_dCompDens[posu][ic];
249  dDensMobRatio_dCompDens[posd][ic] = phaseDens[eru][esru][eiu][0][ip] * dMobRatio_dCompDens[posd][ic];
250  }
251 
252  // 3.d) Compute the final gravity coefficient \x_{up}_{c,p} \rho_p \frac{\lambda_l \lambda_m}{\lambda_T}
253  for( integer ic = 0; ic < NC; ++ic )
254  {
255  applyChainRule( NC,
256  dCompFrac_dCompDens[eru][esru][eiu],
257  dPhaseCompFrac[eru][esru][eiu][0][ip][ic],
258  dPhaseCompFrac_dC,
259  Deriv::dC );
260  upwPhaseGravCoef[ip][k][ic] = phaseCompFrac[eru][esru][eiu][0][ip][ic] * densMobRatio;
261  dUpwPhaseGravCoef_dPres[ip][k][ic][posu] = dPhaseCompFrac[eru][esru][eiu][0][ip][ic][Deriv::dP] * densMobRatio
262  + phaseCompFrac[eru][esru][eiu][0][ip][ic] * dDensMobRatio_dPres[posu];
263  dUpwPhaseGravCoef_dPres[ip][k][ic][posd] = phaseCompFrac[eru][esru][eiu][0][ip][ic] * dDensMobRatio_dPres[posd];
264 
265  for( integer jc = 0; jc < NC; ++jc )
266  {
267  dUpwPhaseGravCoef_dCompDens[ip][k][ic][posu][jc] = dPhaseCompFrac_dC[jc] * densMobRatio
268  + phaseCompFrac[eru][esru][eiu][0][ip][ic] * dDensMobRatio_dCompDens[posu][jc];
269  dUpwPhaseGravCoef_dCompDens[ip][k][ic][posd][jc] = phaseCompFrac[eru][esru][eiu][0][ip][ic] * dDensMobRatio_dCompDens[posd][jc];
270  }
271  }
272  ++k;
273  }
274  }
275 
276 }
277 
278 template< integer NC, integer NP >
279 GEOS_HOST_DEVICE
280 inline
281 void
282 UpwindingHelper::
283  computePhaseGravTerm( localIndex const (&localIds)[ 3 ],
284  localIndex const (&neighborIds)[ 3 ],
285  real64 const & transGravCoef,
286  ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseMassDens,
287  ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
288  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
289  real64 ( & phaseGravTerm )[ NP ][ NP-1 ],
290  real64 ( & dPhaseGravTerm_dPres )[ NP ][ NP-1 ][ 2 ],
291  real64 ( & dPhaseGravTerm_dCompDens )[ NP ][ NP-1 ][ 2 ][ NC ] )
292 {
293  using Deriv = multifluid::DerivativeOffset;
294 
295  localIndex const er = localIds[0];
296  localIndex const esr = localIds[1];
297  localIndex const ei = localIds[2];
298  localIndex const ern = neighborIds[0];
299  localIndex const esrn = neighborIds[1];
300  localIndex const ein = neighborIds[2];
301 
302  real64 dPhaseMassDens_dCLoc[ NC ]{};
303  real64 dPhaseMassDens_dCNeighbor[ NC ]{};
304  real64 dPhaseMassDens_dC[ NC ]{};
305 
306  for( integer ip = 0; ip < NP; ++ip )
307  {
308  applyChainRule( NC,
309  dCompFrac_dCompDens[er][esr][ei],
310  dPhaseMassDens[er][esr][ei][0][ip],
311  dPhaseMassDens_dCLoc,
312  Deriv::dC );
313  applyChainRule( NC,
314  dCompFrac_dCompDens[ern][esrn][ein],
315  dPhaseMassDens[ern][esrn][ein][0][ip],
316  dPhaseMassDens_dCNeighbor,
317  Deriv::dC );
318 
319  localIndex k = 0;
320  for( integer jp = 0; jp < NP; ++jp )
321  {
322  if( ip == jp )
323  {
324  continue;
325  }
326 
327  phaseGravTerm[ip][k] = -( phaseMassDens[er][esr][ei][0][ip] + phaseMassDens[ern][esrn][ein][0][ip] );
328  phaseGravTerm[ip][k] += ( phaseMassDens[er][esr][ei][0][jp] + phaseMassDens[ern][esrn][ein][0][jp] );
329  phaseGravTerm[ip][k] *= 0.5 * transGravCoef;
330 
331  dPhaseGravTerm_dPres[ip][k][Pos::LOCAL] = ( -dPhaseMassDens[er][esr][ei][0][ip][Deriv::dP] + dPhaseMassDens[er][esr][ei][0][jp][Deriv::dP] );
332  dPhaseGravTerm_dPres[ip][k][Pos::LOCAL] *= 0.5 * transGravCoef;
333 
334  dPhaseGravTerm_dPres[ip][k][Pos::NEIGHBOR] = ( -dPhaseMassDens[ern][esrn][ein][0][ip][Deriv::dP] + dPhaseMassDens[ern][esrn][ein][0][jp][Deriv::dP] );
335  dPhaseGravTerm_dPres[ip][k][Pos::NEIGHBOR] *= 0.5 * transGravCoef;
336 
337  for( integer ic = 0; ic < NC; ++ic )
338  {
339  dPhaseGravTerm_dCompDens[ip][k][Pos::LOCAL][ic] = -0.5 * transGravCoef * dPhaseMassDens_dCLoc[ic];
340  dPhaseGravTerm_dCompDens[ip][k][Pos::NEIGHBOR][ic] = -0.5 * transGravCoef * dPhaseMassDens_dCNeighbor[ic];
341  }
342  applyChainRule( NC,
343  dCompFrac_dCompDens[er][esr][ei],
344  dPhaseMassDens[er][esr][ei][0][jp],
345  dPhaseMassDens_dC,
346  Deriv::dC );
347  for( integer ic = 0; ic < NC; ++ic )
348  {
349  dPhaseGravTerm_dCompDens[ip][k][Pos::LOCAL][ic] += 0.5 * transGravCoef * dPhaseMassDens_dC[ic];
350  }
351  applyChainRule( NC,
352  dCompFrac_dCompDens[ern][esrn][ein],
353  dPhaseMassDens[ern][esrn][ein][0][jp],
354  dPhaseMassDens_dC,
355  Deriv::dC );
356  for( integer ic = 0; ic < NC; ++ic )
357  {
358  dPhaseGravTerm_dCompDens[ip][k][Pos::NEIGHBOR][ic] += 0.5 * transGravCoef * dPhaseMassDens_dC[ic];
359  }
360  ++k;
361  }
362  }
363 }
364 
365 template< integer NC, integer NP >
366 GEOS_HOST_DEVICE
367 inline
368 void
369 UpwindingHelper::
370  computeUpwindedTotalMobility( localIndex const (&localIds)[ 3 ],
371  localIndex const (&neighborIds)[ 3 ],
372  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
373  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
374  real64 const (&phaseGravTerm)[ NP ][ NP-1 ],
375  real64 & totalMob,
376  real64 ( & dTotalMob_dPres )[ 2 ],
377  real64 ( & dTotalMob_dCompDens )[ 2 ][ NC ] )
378 {
379  using Deriv = multifluid::DerivativeOffset;
380 
381  localIndex totalMobIds[ NP ][ 3 ]{};
382  localIndex totalMobPos[ NP ]{};
383  setIndicesForTotalMobilityUpwinding< NP >( localIds,
384  neighborIds,
385  phaseGravTerm,
386  totalMobIds,
387  totalMobPos );
388  for( integer ip = 0; ip < NP; ++ip )
389  {
390  localIndex const er = totalMobIds[ip][0];
391  localIndex const esr = totalMobIds[ip][1];
392  localIndex const ei = totalMobIds[ip][2];
393  localIndex const pos = totalMobPos[ip];
394  totalMob = totalMob + phaseMob[er][esr][ei][ip];
395  dTotalMob_dPres[pos] = dTotalMob_dPres[pos] + dPhaseMob[er][esr][ei][pos][Deriv::dP];
396  for( integer ic = 0; ic < NC; ++ic )
397  {
398  dTotalMob_dCompDens[pos][ic] = dTotalMob_dCompDens[pos][ic] + dPhaseMob[er][esr][ei][ip][Deriv::dC+ic];
399  }
400  }
401  if( totalMob < 1e-12 )
402  {
403  totalMob = 1e-12;
404  }
405 }
406 
407 GEOS_HOST_DEVICE
408 inline
409 void
410 UpwindingHelper::
411  setIndicesForMobilityRatioUpwinding( localIndex const (&localIds)[ 3 ],
412  localIndex const (&neighborIds)[ 3 ],
413  real64 const & gravTerm,
414  localIndex & eru, localIndex & esru, localIndex & eiu, localIndex & posu,
415  localIndex & erd, localIndex & esrd, localIndex & eid, localIndex & posd )
416 {
417  if( gravTerm > 0 )
418  {
419  eru = localIds[0];
420  esru = localIds[1];
421  eiu = localIds[2];
422  posu = Pos::LOCAL;
423  erd = neighborIds[0];
424  esrd = neighborIds[1];
425  eid = neighborIds[2];
426  posd = Pos::NEIGHBOR;
427  }
428  else
429  {
430  eru = neighborIds[0];
431  esru = neighborIds[1];
432  eiu = neighborIds[2];
433  posu = Pos::NEIGHBOR;
434  erd = localIds[0];
435  esrd = localIds[1];
436  eid = localIds[2];
437  posd = Pos::LOCAL;
438  }
439 }
440 
441 template< integer NP >
442 GEOS_HOST_DEVICE
443 inline
444 void
445 UpwindingHelper::
446  setIndicesForTotalMobilityUpwinding( localIndex const (&localIds)[ 3 ],
447  localIndex const (&neighborIds)[ 3 ],
448  real64 const (&gravTerm)[ NP ][ NP-1 ],
449  localIndex ( & totalMobIds )[ NP ][ 3 ],
450  localIndex ( & totalMobPos )[ NP ] )
451 {
452  if constexpr ( NP == 2 )
453  {
454  if( gravTerm[0][0] > 0 )
455  {
456  totalMobIds[0][0] = localIds[0];
457  totalMobIds[0][1] = localIds[1];
458  totalMobIds[0][2] = localIds[2];
459  totalMobPos[0] = Pos::LOCAL;
460  totalMobIds[1][0] = neighborIds[0];
461  totalMobIds[1][1] = neighborIds[1];
462  totalMobIds[1][2] = neighborIds[2];
463  totalMobPos[1] = Pos::NEIGHBOR;
464  }
465  else
466  {
467  totalMobIds[0][0] = neighborIds[0];
468  totalMobIds[0][1] = neighborIds[1];
469  totalMobIds[0][2] = neighborIds[2];
470  totalMobPos[0] = Pos::NEIGHBOR;
471  totalMobIds[1][0] = localIds[0];
472  totalMobIds[1][1] = localIds[1];
473  totalMobIds[1][2] = localIds[2];
474  totalMobPos[1] = Pos::LOCAL;
475  }
476  }
477  else if constexpr ( NP == 3 )
478  {
479  // TODO Francois: this should be improved
480  // currently this implements the algorithm proposed by SH Lee
481  for( integer ip = 0; ip < NP; ++ip )
482  {
483  if( ( gravTerm[ip][0] >= 0 && gravTerm[ip][1] >= 0 ) || // includes the no-buoyancy case
484  ( ( LvArray::math::abs( gravTerm[ip][0] ) >= LvArray::math::abs( gravTerm[ip][1] ) ) && gravTerm[ip][1] >= 0 ) ||
485  ( ( LvArray::math::abs( gravTerm[ip][1] ) >= LvArray::math::abs( gravTerm[ip][0] ) ) && gravTerm[ip][0] >= 0 ) )
486  {
487  totalMobIds[ip][0] = localIds[0];
488  totalMobIds[ip][1] = localIds[1];
489  totalMobIds[ip][2] = localIds[2];
490  totalMobPos[ip] = Pos::LOCAL;
491  }
492  else
493  {
494  totalMobIds[ip][0] = neighborIds[0];
495  totalMobIds[ip][1] = neighborIds[1];
496  totalMobIds[ip][2] = neighborIds[2];
497  totalMobPos[ip] = Pos::NEIGHBOR;
498  }
499  }
500  }
501 }
502 
503 /******************************** AssemblerKernelHelper ********************************/
504 
505 template< integer NF, integer NC, integer NP >
506 GEOS_HOST_DEVICE
507 inline
508 void
509 AssemblerKernelHelper::
510  applyGradient( arrayView1d< real64 const > const & facePres,
511  arrayView1d< real64 const > const & faceGravCoef,
512  arraySlice1d< localIndex const > const & elemToFaces,
513  real64 const & elemPres,
514  real64 const & elemGravCoef,
515  arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const & elemPhaseMassDens,
516  arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const & dElemPhaseMassDens,
517  arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const & elemPhaseMob,
518  arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > const & dElemPhaseMob,
519  arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > const & dElemCompFrac_dCompDens,
520  arraySlice2d< real64 const > const & transMatrix,
521  real64 ( & oneSidedVolFlux )[ NF ],
522  real64 ( & dOneSidedVolFlux_dPres )[ NF ],
523  real64 ( & dOneSidedVolFlux_dFacePres )[ NF ][ NF ],
524  real64 ( & dOneSidedVolFlux_dCompDens )[ NF ][ NC ] )
525 {
526  using Deriv = multifluid::DerivativeOffset;
527 
528  real64 dPhaseMassDens_dC[ NP ][ NC ]{};
529  real64 dPresDif_dCompDens[ NC ]{};
530  real64 dPhaseGravDif_dCompDens[ NC ]{};
531  real64 dPhaseMobPotDif_dCompDens[ NC ]{};
532 
533  // 0) precompute dPhaseDens_dC since it is always computed at the element center
534  for( integer ip = 0; ip < NP; ++ip )
535  {
536  applyChainRule( NC,
537  dElemCompFrac_dCompDens,
538  dElemPhaseMassDens[ip],
539  dPhaseMassDens_dC[ip],
540  Deriv::dC );
541  }
542 
543  for( integer ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc )
544  {
545  // now in the following nested loop,
546  // we compute the contribution of face jfaceLoc to the one sided total volumetric flux at face iface
547  for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc )
548  {
549 
550  // depth difference between element center and face center
551  real64 const ccGravCoef = elemGravCoef;
552  real64 const fGravCoef = faceGravCoef[elemToFaces[jfaceLoc]];
553  real64 const gravCoefDif = ccGravCoef - fGravCoef;
554 
555  for( integer ip = 0; ip < NP; ++ip )
556  {
557 
558  // 1) compute the potential diff between the cell center and the face center
559  real64 const ccPres = elemPres;
560  real64 const fPres = facePres[elemToFaces[jfaceLoc]];
561 
562  // pressure difference
563  real64 const presDif = ccPres - fPres;
564  real64 const dPresDif_dPres = 1;
565  real64 const dPresDif_dFacePres = -1;
566  for( integer ic = 0; ic < NC; ++ic )
567  {
568  dPresDif_dCompDens[ic] = 0.0; // no capillary pressure
569  }
570 
571  // gravity term
572  real64 const phaseGravDif = elemPhaseMassDens[ip] * gravCoefDif;
573  real64 const dPhaseGravDif_dPres = dElemPhaseMassDens[ip][Deriv::dP] * gravCoefDif;
574  for( localIndex ic = 0; ic < NC; ++ic )
575  {
576  dPhaseGravDif_dCompDens[ic] = dPhaseMassDens_dC[ip][ic] * gravCoefDif;
577  }
578  // no density evaluated at the face center
579 
580  // potential difference
581  real64 const phasePotDif = presDif - phaseGravDif;
582  real64 const phaseMobPotDif = elemPhaseMob[ip] * phasePotDif;
583  real64 const dPhaseMobPotDif_dPres = dElemPhaseMob[ip][Deriv::dP] * phasePotDif
584  + elemPhaseMob[ip] * (dPresDif_dPres - dPhaseGravDif_dPres);
585  real64 const dPhaseMobPotDif_dFacePres = elemPhaseMob[ip] * dPresDif_dFacePres;
586  for( integer ic = 0; ic < NC; ++ic )
587  {
588  dPhaseMobPotDif_dCompDens[ic] = dElemPhaseMob[ip][Deriv::dC+ic] * phasePotDif
589  + elemPhaseMob[ip] * (dPresDif_dCompDens[ic] - dPhaseGravDif_dCompDens[ic]);
590  }
591 
592  // this is going to store T \sum_p \lambda_p (\nabla p - \rho_p g \nabla d)
593  oneSidedVolFlux[ifaceLoc] = oneSidedVolFlux[ifaceLoc]
594  + transMatrix[ifaceLoc][jfaceLoc] * phaseMobPotDif;
595  dOneSidedVolFlux_dPres[ifaceLoc] = dOneSidedVolFlux_dPres[ifaceLoc]
596  + transMatrix[ifaceLoc][jfaceLoc] * dPhaseMobPotDif_dPres;
597  dOneSidedVolFlux_dFacePres[ifaceLoc][jfaceLoc] = dOneSidedVolFlux_dFacePres[ifaceLoc][jfaceLoc]
598  + transMatrix[ifaceLoc][jfaceLoc] * dPhaseMobPotDif_dFacePres;
599  for( localIndex ic = 0; ic < NC; ++ic )
600  {
601  dOneSidedVolFlux_dCompDens[ifaceLoc][ic] = dOneSidedVolFlux_dCompDens[ifaceLoc][ic]
602  + transMatrix[ifaceLoc][jfaceLoc] * dPhaseMobPotDif_dCompDens[ic];
603  }
604  }
605  }
606  }
607 }
608 
609 template< integer NF, integer NC, integer NP >
610 GEOS_HOST_DEVICE
611 inline
612 void
613 AssemblerKernelHelper::
614  assembleFluxDivergence( localIndex const (&localIds)[ 3 ],
615  globalIndex const rankOffset,
616  arrayView2d< localIndex const > const & elemRegionList,
617  arrayView2d< localIndex const > const & elemSubRegionList,
618  arrayView2d< localIndex const > const & elemList,
619  SortedArrayView< localIndex const > const & regionFilter,
620  arrayView1d< globalIndex const > const & faceDofNumber,
621  arrayView1d< integer const > const & isBoundaryFace,
622  arrayView1d< real64 const > const & mimFaceGravCoef,
623  arraySlice1d< localIndex const > const & elemToFaces,
624  real64 const & elemGravCoef,
625  integer const useTotalMassEquation,
626  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseDens,
627  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseDens,
628  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens,
629  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
630  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
631  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
632  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
633  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const & phaseCompFrac,
634  ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac,
635  ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber,
636  arraySlice2d< real64 const > const & transMatrixGrav,
637  real64 const (&oneSidedVolFlux)[ NF ],
638  real64 const (&dOneSidedVolFlux_dPres)[ NF ],
639  real64 const (&dOneSidedVolFlux_dFacePres)[ NF ][ NF ],
640  real64 const (&dOneSidedVolFlux_dCompDens)[ NF ][ NC ],
641  real64 const & dt,
642  CRSMatrixView< real64, globalIndex const > const & localMatrix,
643  arrayView1d< real64 > const & localRhs )
644 {
645  using namespace compositionalMultiphaseUtilities;
646  integer constexpr NDOF = NC+1;
647 
648  // dof numbers
649  globalIndex dofColIndicesElemVars[ NDOF*(NF+1) ]{};
650  globalIndex dofColIndicesFaceVars[ NF ]{};
651  integer faceIndexMap[ NF ]{}; // Maps local face index to position in compact DOF array
652  integer numNonBoundaryFaces = 0;
653  for( integer idof = 0; idof < NDOF; ++idof )
654  {
655  dofColIndicesElemVars[idof] = elemDofNumber[localIds[0]][localIds[1]][localIds[2]] + idof;
656  }
657 
658  // divergence of fluxes
659  real64 divMassFluxes[ NC ]{};
660  real64 dDivMassFluxes_dElemVars[ NC ][ NDOF*(NF+1) ]{};
661  real64 dDivMassFluxes_dFaceVars[ NC ][ NF ]{};
662 
663  // auxiliary variables for upwinding
664 
665  // upwinding phase buoyancy transport coefficients
666  real64 upwPhaseViscCoef[ NP ][ NC ]{};
667  real64 dUpwPhaseViscCoef_dPres[ NP ][ NC ]{};
668  real64 dUpwPhaseViscCoef_dCompDens[ NP ][ NC ][ NC ]{};
669  globalIndex upwViscDofNumber = 0;
670 
671  // gravity term: ( \rho_l - \rho_m ) g \Delta z
672  real64 phaseGravTerm[ NP ][ NP-1 ]{};
673  real64 dPhaseGravTerm_dPres[ NP ][ NP-1 ][ 2 ]{};
674  real64 dPhaseGravTerm_dCompDens[ NP ][ NP-1 ][ 2 ][ NC ]{};
675 
676  // upwinding phase buoyancy transport coefficients
677  real64 upwPhaseGravCoef[ NP ][ NP-1 ][ NC ]{};
678  real64 dUpwPhaseGravCoef_dPres[ NP ][ NP-1 ][ NC ][ 2 ]{};
679  real64 dUpwPhaseGravCoef_dCompDens[ NP ][ NP-1 ][ NC ][ 2 ][ NC ]{};
680 
681  // for each element, loop over the one-sided faces
682  for( integer ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc )
683  {
684 
685  // Contribute and Collect face DOF number only if this is not a boundary face
686  if( isBoundaryFace[elemToFaces[ifaceLoc]] == 0 )
687  {
688 
689  faceIndexMap[ifaceLoc] = numNonBoundaryFaces;
690  dofColIndicesFaceVars[numNonBoundaryFaces] = faceDofNumber[elemToFaces[ifaceLoc]];
691  numNonBoundaryFaces++;
692 
693 
694  // 1) Find if there is a neighbor, and if there is, grab the indices of the neighbor element
695  localIndex neighborIds[ 3 ] = { localIds[0], localIds[1], localIds[2] };
696  hybridFVMKernels::CellConnectivity::isNeighborFound( localIds,
697  ifaceLoc,
698  elemRegionList,
699  elemSubRegionList,
700  elemList,
701  regionFilter,
702  elemToFaces,
703  neighborIds );
704  localIndex const neighborDofNumber = elemDofNumber[neighborIds[0]][neighborIds[1]][neighborIds[2]];
705 
706  // 2) *************** Assemble viscous terms ******************
707 
708  // 2.a) Compute the upwinded x_{c, \ell} \rho_{\ell} \frac{\lambda_{\ell}}{\lambda_T} for each phase at this face
709  UpwindingHelper::upwindViscousCoefficient< NC, NP >( localIds,
710  neighborIds,
711  phaseDens,
712  dPhaseDens,
713  phaseMob,
714  dPhaseMob,
715  dCompFrac_dCompDens,
716  phaseCompFrac,
717  dPhaseCompFrac,
718  elemDofNumber,
719  oneSidedVolFlux[ifaceLoc],
720  upwPhaseViscCoef,
721  dUpwPhaseViscCoef_dPres,
722  dUpwPhaseViscCoef_dCompDens,
723  upwViscDofNumber );
724 
725  // 2.b) Add the \x_{c,\ell} \rho_{\ell} \frac{\lambda_{\ell}}{\lambda_T} q_T of this face to the divergence of the flux in this cell
726  assembleViscousFlux< NF, NC, NP >( ifaceLoc,
727  oneSidedVolFlux,
728  dOneSidedVolFlux_dPres,
729  dOneSidedVolFlux_dFacePres,
730  dOneSidedVolFlux_dCompDens,
731  upwPhaseViscCoef,
732  dUpwPhaseViscCoef_dPres,
733  dUpwPhaseViscCoef_dCompDens,
734  elemDofNumber[localIds[0]][localIds[1]][localIds[2]],
735  neighborDofNumber,
736  upwViscDofNumber,
737  dt,
738  divMassFluxes,
739  dDivMassFluxes_dElemVars,
740  dDivMassFluxes_dFaceVars,
741  dofColIndicesElemVars );
742 
743 
744  // 3) *************** Assemble buoyancy terms ******************
745 
746  real64 const transGravCoef = (localIds[0] != neighborIds[0] || localIds[1] != neighborIds[1] || localIds[2] != neighborIds[2])
747  * transMatrixGrav[ifaceLoc][ifaceLoc] * (elemGravCoef - mimFaceGravCoef[elemToFaces[ifaceLoc]]);
748 
749  // 3.a) Compute the upwinded x_{c, \ell} \rho_{\ell} \frac{\lambda_{\ell}\lambda_m}{\lambda_T}
750  // and (\rho_{\ell} - \rho_m) g \Delta z for each phase at this face
751  UpwindingHelper::upwindBuoyancyCoefficient< NC, NP >( localIds,
752  neighborIds,
753  transGravCoef,
754  phaseDens,
755  dPhaseDens,
756  phaseMassDens,
757  dPhaseMassDens,
758  phaseMob,
759  dPhaseMob,
760  dCompFrac_dCompDens,
761  phaseCompFrac,
762  dPhaseCompFrac,
763  phaseGravTerm,
764  dPhaseGravTerm_dPres,
765  dPhaseGravTerm_dCompDens,
766  upwPhaseGravCoef,
767  dUpwPhaseGravCoef_dPres,
768  dUpwPhaseGravCoef_dCompDens );
769 
770  // 3.b) Add the buoyancy term of this face to the divergence of the flux in this cell
771  assembleBuoyancyFlux< NF, NC, NP >( ifaceLoc,
772  phaseGravTerm,
773  dPhaseGravTerm_dPres,
774  dPhaseGravTerm_dCompDens,
775  upwPhaseGravCoef,
776  dUpwPhaseGravCoef_dPres,
777  dUpwPhaseGravCoef_dCompDens,
778  dt,
779  divMassFluxes,
780  dDivMassFluxes_dElemVars );
781 
782  }
783  else
784  {
785  faceIndexMap[ifaceLoc] = -1; // Mark boundary faces
786  }
787 
788  }
789 
790  if( useTotalMassEquation > 0 )
791  {
792  // Apply equation/variable change transformation(s)
793  real64 work[NDOF * ( NF + 1 )];
794  shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NDOF * ( NF + 1 ), dDivMassFluxes_dElemVars, work );
795  shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NF, dDivMassFluxes_dFaceVars, work );
796  shiftElementsAheadByOneAndReplaceFirstElementWithSum( NC, divMassFluxes );
797  }
798 
799  // we are ready to assemble the local flux and its derivatives
800  // no need for atomic adds - each row is assembled by a single thread
801 
802  for( integer ic = 0; ic < NC; ++ic )
803  {
804  localIndex const eqnRowLocalIndex =
805  LvArray::integerConversion< localIndex >( elemDofNumber[localIds[0]][localIds[1]][localIds[2]] + ic - rankOffset );
806 
807  GEOS_ASSERT_GE( eqnRowLocalIndex, 0 );
808  GEOS_ASSERT_GT( localMatrix.numRows(), eqnRowLocalIndex );
809 
810  // residual
811  localRhs[eqnRowLocalIndex] = localRhs[eqnRowLocalIndex] + divMassFluxes[ic];
812 
813  // jacobian -- derivative wrt elem centered vars
814  // Need to compact both DOF indices and derivatives to only include non-boundary neighbors
815  globalIndex compactElemDofs[ NDOF*(NF+1) ];
816  real64 compactElemDerivs[ NDOF*(NF+1) ];
817 
818  // Copy current element DOFs and derivatives
819  for( integer idof = 0; idof < NDOF; ++idof )
820  {
821  compactElemDofs[idof] = dofColIndicesElemVars[idof];
822  compactElemDerivs[idof] = dDivMassFluxes_dElemVars[ic][idof];
823  }
824 
825  // Copy neighbor DOFs and derivatives only for non-boundary faces
826  integer compactIdx = NDOF;
827  for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc )
828  {
829  if( faceIndexMap[jfaceLoc] >= 0 ) // non-boundary face
830  {
831  integer const neighborOffset = NDOF * (jfaceLoc + 1);
832  for( integer idof = 0; idof < NDOF; ++idof )
833  {
834  compactElemDofs[compactIdx] = dofColIndicesElemVars[neighborOffset + idof];
835  compactElemDerivs[compactIdx] = dDivMassFluxes_dElemVars[ic][neighborOffset + idof];
836  compactIdx++;
837  }
838  }
839  }
840 
841  localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( eqnRowLocalIndex,
842  compactElemDofs,
843  compactElemDerivs,
844  compactIdx );
845 
846  // jacobian -- derivatives wrt face centered vars (only non-boundary faces)
847  // Need to compact the derivatives array to only include non-boundary faces
848  if( numNonBoundaryFaces > 0 )
849  {
850  real64 compactFaceDerivs[ NF ]{};
851  for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc )
852  {
853  if( faceIndexMap[jfaceLoc] >= 0 )
854  {
855  compactFaceDerivs[faceIndexMap[jfaceLoc]] = dDivMassFluxes_dFaceVars[ic][jfaceLoc];
856  }
857  }
858 
859  localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( eqnRowLocalIndex,
860  &dofColIndicesFaceVars[0],
861  compactFaceDerivs,
862  numNonBoundaryFaces );
863  }
864  }
865 }
866 
867 template< integer NF, integer NC, integer NP >
868 GEOS_HOST_DEVICE
869 inline
870 void
871 AssemblerKernelHelper::
872  assembleViscousFlux( localIndex const ifaceLoc,
873  real64 const (&oneSidedVolFlux)[ NF ],
874  real64 const (&dOneSidedVolFlux_dPres)[ NF ],
875  real64 const (&dOneSidedVolFlux_dFacePres)[ NF ][ NF ],
876  real64 const (&dOneSidedVolFlux_dCompDens)[ NF ][ NC ],
877  real64 const (&upwPhaseViscCoef)[ NP ][ NC ],
878  real64 const (&dUpwPhaseViscCoef_dPres)[ NP ][ NC ],
879  real64 const (&dUpwPhaseViscCoef_dCompDens)[ NP ][ NC ][ NC ],
880  globalIndex const elemDofNumber,
881  globalIndex const neighborDofNumber,
882  globalIndex const upwViscDofNumber,
883  real64 const & dt,
884  real64 ( & divMassFluxes )[ NC ],
885  real64 ( & dDivMassFluxes_dElemVars )[ NC ][ (NC+1)*(NF+1) ],
886  real64 ( & dDivMassFluxes_dFaceVars )[ NC ][ NF ],
887  globalIndex ( & dofColIndicesElemVars )[ (NC+1)*(NF+1) ] )
888 {
889  integer constexpr NDOF = NC+1;
890  localIndex const elemVarsOffset = NDOF*(ifaceLoc+1);
891 
892  for( integer ip = 0; ip < NP; ++ip )
893  {
894  for( integer ic = 0; ic < NC; ++ic )
895  {
896  // compute the mass flux at the one-sided face plus its derivatives
897  // add the newly computed flux to the sum
898 
899  real64 const dt_upwPhaseViscCoef = dt * upwPhaseViscCoef[ip][ic];
900 
901  // residual
902  divMassFluxes[ic] = divMassFluxes[ic] + dt_upwPhaseViscCoef * oneSidedVolFlux[ifaceLoc];
903 
904  // local derivatives
905  dDivMassFluxes_dElemVars[ic][0] = dDivMassFluxes_dElemVars[ic][0]
906  + dt_upwPhaseViscCoef * dOneSidedVolFlux_dPres[ifaceLoc];
907  dDivMassFluxes_dElemVars[ic][0] = dDivMassFluxes_dElemVars[ic][0]
908  + ( elemDofNumber == upwViscDofNumber )
909  * dt * dUpwPhaseViscCoef_dPres[ip][ic] * oneSidedVolFlux[ifaceLoc];
910  for( integer jc = 0; jc < NC; ++jc )
911  {
912  dDivMassFluxes_dElemVars[ic][jc+1] = dDivMassFluxes_dElemVars[ic][jc+1]
913  + dt_upwPhaseViscCoef * dOneSidedVolFlux_dCompDens[ifaceLoc][jc];
914  dDivMassFluxes_dElemVars[ic][jc+1] = dDivMassFluxes_dElemVars[ic][jc+1]
915  + ( elemDofNumber == upwViscDofNumber )
916  * dt * dUpwPhaseViscCoef_dCompDens[ip][ic][jc] * oneSidedVolFlux[ifaceLoc];
917  }
918 
919  // neighbor derivatives
920  dDivMassFluxes_dElemVars[ic][elemVarsOffset] = dDivMassFluxes_dElemVars[ic][elemVarsOffset]
921  + ( elemDofNumber != upwViscDofNumber )
922  * dt * dUpwPhaseViscCoef_dPres[ip][ic] * oneSidedVolFlux[ifaceLoc];
923 
924  for( integer jc = 0; jc < NC; ++jc )
925  {
926  dDivMassFluxes_dElemVars[ic][elemVarsOffset+jc+1] = dDivMassFluxes_dElemVars[ic][elemVarsOffset+jc+1]
927  + ( elemDofNumber != upwViscDofNumber )
928  * dt * dUpwPhaseViscCoef_dCompDens[ip][ic][jc] * oneSidedVolFlux[ifaceLoc];
929  }
930 
931  for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc )
932  {
933  dDivMassFluxes_dFaceVars[ic][jfaceLoc] = dDivMassFluxes_dFaceVars[ic][jfaceLoc]
934  + dt_upwPhaseViscCoef * dOneSidedVolFlux_dFacePres[ifaceLoc][jfaceLoc];
935  }
936  }
937  }
938 
939  // collect the relevant dof numbers
940  for( integer idof = 0; idof < NDOF; ++idof )
941  {
942  dofColIndicesElemVars[elemVarsOffset+idof] = neighborDofNumber + idof;
943  }
944 }
945 
946 template< integer NF, integer NC, integer NP >
947 GEOS_HOST_DEVICE
948 inline
949 void
950 AssemblerKernelHelper::
951  assembleBuoyancyFlux( localIndex const ifaceLoc,
952  real64 const (&phaseGravTerm)[ NP ][ NP-1 ],
953  real64 const (&dPhaseGravTerm_dPres)[ NP ][ NP-1 ][ 2 ],
954  real64 const (&dPhaseGravTerm_dCompDens)[ NP ][ NP-1 ][ 2 ][ NC ],
955  real64 const (&upwPhaseGravCoef)[ NP ][ NP-1 ][ NC ],
956  real64 const (&dUpwPhaseGravCoef_dPres)[ NP ][ NP-1 ][ NC ][ 2 ],
957  real64 const (&dUpwPhaseGravCoef_dCompDens)[ NP ][ NP-1 ][ NC ][ 2 ][ NC ],
958  real64 const & dt,
959  real64 ( & divMassFluxes )[ NC ],
960  real64 ( & dDivMassFluxes_dElemVars )[ NC ][ (NC+1)*(NF+1) ] )
961 {
962  integer constexpr NDOF = NC+1;
963  localIndex const elemVarsOffset = NDOF*(ifaceLoc+1);
964 
965  for( integer ip = 0; ip < NP; ++ip )
966  {
967  for( integer jp = 0; jp < NP - 1; ++jp )
968  {
969  for( integer ic = 0; ic < NC; ++ic )
970  {
971  real64 const dt_upwPhaseGravCoef = dt * upwPhaseGravCoef[ip][jp][ic];
972 
973  // residual
974  divMassFluxes[ic] = divMassFluxes[ic] + dt_upwPhaseGravCoef * phaseGravTerm[ip][jp];
975 
976  // local derivatives
977  dDivMassFluxes_dElemVars[ic][0] = dDivMassFluxes_dElemVars[ic][0]
978  + dt * dUpwPhaseGravCoef_dPres[ip][jp][ic][Pos::LOCAL] * phaseGravTerm[ip][jp];
979  dDivMassFluxes_dElemVars[ic][0] = dDivMassFluxes_dElemVars[ic][0]
980  + dt_upwPhaseGravCoef * dPhaseGravTerm_dPres[ip][jp][Pos::LOCAL];
981 
982  for( integer jc = 0; jc < NC; ++jc )
983  {
984  dDivMassFluxes_dElemVars[ic][jc+1] = dDivMassFluxes_dElemVars[ic][jc+1]
985  + dt * dUpwPhaseGravCoef_dCompDens[ip][jp][ic][Pos::LOCAL][jc] * phaseGravTerm[ip][jp];
986  dDivMassFluxes_dElemVars[ic][jc+1] = dDivMassFluxes_dElemVars[ic][jc+1]
987  + dt_upwPhaseGravCoef * dPhaseGravTerm_dCompDens[ip][jp][Pos::LOCAL][jc];
988  }
989 
990  // neighbor derivatives
991  dDivMassFluxes_dElemVars[ic][elemVarsOffset] = dDivMassFluxes_dElemVars[ic][elemVarsOffset]
992  + dt * dUpwPhaseGravCoef_dPres[ip][jp][ic][Pos::NEIGHBOR] * phaseGravTerm[ip][jp];
993  dDivMassFluxes_dElemVars[ic][elemVarsOffset] = dDivMassFluxes_dElemVars[ic][elemVarsOffset]
994  + dt_upwPhaseGravCoef * dPhaseGravTerm_dPres[ip][jp][Pos::NEIGHBOR];
995 
996  for( integer jc = 0; jc < NC; ++jc )
997  {
998  dDivMassFluxes_dElemVars[ic][elemVarsOffset+jc+1] = dDivMassFluxes_dElemVars[ic][elemVarsOffset+jc+1]
999  + dt * dUpwPhaseGravCoef_dCompDens[ip][jp][ic][Pos::NEIGHBOR][jc] * phaseGravTerm[ip][jp];
1000  dDivMassFluxes_dElemVars[ic][elemVarsOffset+jc+1] = dDivMassFluxes_dElemVars[ic][elemVarsOffset+jc+1]
1001  + dt_upwPhaseGravCoef * dPhaseGravTerm_dCompDens[ip][jp][Pos::NEIGHBOR][jc];
1002  }
1003  }
1004  }
1005  }
1006 }
1007 
1008 template< integer NF, integer NC >
1009 GEOS_HOST_DEVICE
1010 inline
1011 void
1012 AssemblerKernelHelper::
1013  assembleFaceConstraints( arrayView1d< globalIndex const > const & faceDofNumber,
1014  arrayView1d< integer const > const & faceGhostRank,
1015  arrayView1d< integer const > const & isBoundaryFace,
1016  arraySlice1d< localIndex const > const & elemToFaces,
1017  globalIndex const elemDofNumber,
1018  globalIndex const rankOffset,
1019  real64 const (&oneSidedVolFlux)[ NF ],
1020  real64 const (&dOneSidedVolFlux_dPres)[ NF ],
1021  real64 const (&dOneSidedVolFlux_dFacePres)[ NF ][ NF ],
1022  real64 const (&dOneSidedVolFlux_dCompDens)[ NF ][ NC ],
1023  CRSMatrixView< real64, globalIndex const > const & localMatrix,
1024  arrayView1d< real64 > const & localRhs )
1025 {
1026  integer constexpr NDOF = NC+1;
1027 
1028  // fluxes
1029  real64 dFlux_dElemVars[ NDOF ]{};
1030  real64 dFlux_dFaceVars[ NF ]{};
1031 
1032  // dof numbers
1033  globalIndex dofColIndicesElemVars[ NDOF ]{};
1034  globalIndex dofColIndicesFaceVars[ NF ]{};
1035  for( integer idof = 0; idof < NDOF; ++idof )
1036  {
1037  dofColIndicesElemVars[idof] = elemDofNumber + idof;
1038  }
1039 
1040  // for each element, loop over the local (one-sided) faces
1041  for( integer ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc )
1042  {
1043  localIndex const kf = elemToFaces[ifaceLoc];
1044 
1045  // Skip ghost faces
1046  if( faceGhostRank[kf] >= 0 )
1047  {
1048  continue;
1049  }
1050 
1051  // Skip boundary faces - they have Dirichlet constraints, not flux continuity
1052  if( isBoundaryFace[kf] > 0 )
1053  {
1054  continue;
1055  }
1056 
1057  // flux at this face
1058  real64 const flux = oneSidedVolFlux[ifaceLoc];
1059  dFlux_dElemVars[0] = dOneSidedVolFlux_dPres[ifaceLoc];
1060  for( integer ic = 0; ic < NC; ++ic )
1061  {
1062  dFlux_dElemVars[ic+1] = dOneSidedVolFlux_dCompDens[ifaceLoc][ic];
1063  }
1064 
1065  for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc )
1066  {
1067  dFlux_dFaceVars[jfaceLoc] = dOneSidedVolFlux_dFacePres[ifaceLoc][jfaceLoc];
1068  dofColIndicesFaceVars[jfaceLoc] = faceDofNumber[elemToFaces[jfaceLoc]];
1069  }
1070 
1071  // dof number of this face constraint
1072  localIndex const eqnLocalRowIndex = LvArray::integerConversion< localIndex >( faceDofNumber[elemToFaces[ifaceLoc]] - rankOffset );
1073 
1074  GEOS_ASSERT_GE( eqnLocalRowIndex, 0 );
1075  GEOS_ASSERT_GT( localMatrix.numRows(), eqnLocalRowIndex );
1076 
1077  // residual
1078  RAJA::atomicAdd( parallelDeviceAtomic{}, &localRhs[eqnLocalRowIndex], flux );
1079 
1080  // jacobian -- derivatives wrt elem-centered terms
1081  localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( eqnLocalRowIndex,
1082  &dofColIndicesElemVars[0],
1083  &dFlux_dElemVars[0],
1084  NDOF );
1085 
1086  // jacobian -- derivatives wrt face pressure terms
1087  localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( eqnLocalRowIndex,
1088  &dofColIndicesFaceVars[0],
1089  &dFlux_dFaceVars[0],
1090  NF );
1091  }
1092 }
1093 
1094 /******************************** AssemblerKernel ********************************/
1095 
1096 template< integer NF, integer NC, integer NP >
1097 GEOS_HOST_DEVICE
1098 inline
1099 void
1100 AssemblerKernel::
1101  compute( localIndex const er, localIndex const esr, localIndex const ei,
1102  SortedArrayView< localIndex const > const & regionFilter,
1103  arrayView2d< localIndex const > const & elemRegionList,
1104  arrayView2d< localIndex const > const & elemSubRegionList,
1105  arrayView2d< localIndex const > const & elemList,
1106  arrayView1d< globalIndex const > const & faceDofNumber,
1107  arrayView1d< integer const > const & faceGhostRank,
1108  arrayView1d< integer const > const & isBoundaryFace,
1109  arrayView1d< real64 const > const & facePres,
1110  arrayView1d< real64 const > const & faceGravCoef,
1111  arrayView1d< real64 const > const & mimFaceGravCoef,
1112  arraySlice1d< localIndex const > const & elemToFaces,
1113  real64 const & elemPres,
1114  real64 const & elemGravCoef,
1115  integer const useTotalMassEquation,
1116  ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseDens,
1117  ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens,
1118  ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseMassDens,
1119  ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1120  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1121  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
1122  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1123  ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_COMP > > const & phaseCompFrac,
1124  ElementViewConst< arrayView5d< real64 const, multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac,
1125  ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber,
1126  integer const elemGhostRank,
1127  globalIndex const rankOffset,
1128  real64 const & dt,
1129  arraySlice2d< real64 const > const & transMatrix,
1130  arraySlice2d< real64 const > const & transMatrixGrav,
1131  CRSMatrixView< real64, globalIndex const > const & localMatrix,
1132  arrayView1d< real64 > const & localRhs )
1133 {
1134  // one sided flux
1135  real64 oneSidedVolFlux[ NF ]{};
1136  real64 dOneSidedVolFlux_dPres[ NF ]{};
1137  real64 dOneSidedVolFlux_dFacePres[ NF ][ NF ]{};
1138  real64 dOneSidedVolFlux_dCompDens[ NF ][ NC ]{};
1139 
1140  localIndex const localIds[3] = { er, esr, ei };
1141 
1142  /*
1143  * compute auxiliary quantities at the one sided faces of this element:
1144  * 1) One-sided volumetric fluxes
1145  * 2) Upwinded mobilities
1146  */
1147 
1148  // for each one-sided face of the elem,
1149  // compute the volumetric flux using transMatrix
1150  AssemblerKernelHelper::applyGradient< NF, NC, NP >( facePres,
1151  faceGravCoef,
1152  elemToFaces,
1153  elemPres,
1154  elemGravCoef,
1155  phaseMassDens[er][esr][ei][0],
1156  dPhaseMassDens[er][esr][ei][0],
1157  phaseMob[er][esr][ei],
1158  dPhaseMob[er][esr][ei],
1159  dCompFrac_dCompDens[er][esr][ei],
1160  transMatrix,
1161  oneSidedVolFlux,
1162  dOneSidedVolFlux_dPres,
1163  dOneSidedVolFlux_dFacePres,
1164  dOneSidedVolFlux_dCompDens );
1165 
1166  // at this point, we know the local flow direction in the element
1167  // so we can upwind the transport coefficients (mobilities) at the one sided faces
1168  // ** this function needs non-local information **
1169  if( elemGhostRank < 0 )
1170  {
1171  /*
1172  * perform assembly in this element in two steps:
1173  * 1) mass conservation equations
1174  * 2) face constraints
1175  */
1176 
1177  // use the computed one sided vol fluxes and the upwinded mobilities
1178  // to assemble the upwinded mass fluxes in the mass conservation eqn of the elem
1179  AssemblerKernelHelper::assembleFluxDivergence< NF, NC, NP >( localIds,
1180  rankOffset,
1181  elemRegionList,
1182  elemSubRegionList,
1183  elemList,
1184  regionFilter,
1185  faceDofNumber,
1186  isBoundaryFace,
1187  mimFaceGravCoef,
1188  elemToFaces,
1189  elemGravCoef,
1190  useTotalMassEquation,
1191  phaseDens,
1192  dPhaseDens,
1193  phaseMassDens,
1194  dPhaseMassDens,
1195  phaseMob,
1196  dPhaseMob,
1197  dCompFrac_dCompDens,
1198  phaseCompFrac,
1199  dPhaseCompFrac,
1200  elemDofNumber,
1201  transMatrixGrav,
1202  oneSidedVolFlux,
1203  dOneSidedVolFlux_dPres,
1204  dOneSidedVolFlux_dFacePres,
1205  dOneSidedVolFlux_dCompDens,
1206  dt,
1207  localMatrix,
1208  localRhs );
1209  }
1210 
1211  // use the computed one sided vol fluxes to assemble the constraints
1212  // enforcing flux continuity at this element's faces
1213  AssemblerKernelHelper::assembleFaceConstraints< NF, NC >( faceDofNumber,
1214  faceGhostRank,
1215  isBoundaryFace,
1216  elemToFaces,
1217  elemDofNumber[er][esr][ei],
1218  rankOffset,
1219  oneSidedVolFlux,
1220  dOneSidedVolFlux_dPres,
1221  dOneSidedVolFlux_dFacePres,
1222  dOneSidedVolFlux_dCompDens,
1223  localMatrix,
1224  localRhs );
1225 }
1226 
1227 /******************************** FluxKernel ********************************/
1228 
1229 template< integer NF, integer NC, integer NP, typename IP_TYPE >
1230 void
1231 FluxKernel::
1232  launch( localIndex er, localIndex esr,
1233  CellElementSubRegion const & subRegion,
1234  constitutive::PermeabilityBase const & permeabilityModel,
1235  SortedArrayView< localIndex const > const & regionFilter,
1236  arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition,
1237  arrayView2d< localIndex const > const & elemRegionList,
1238  arrayView2d< localIndex const > const & elemSubRegionList,
1239  arrayView2d< localIndex const > const & elemList,
1240  ArrayOfArraysView< localIndex const > const & faceToNodes,
1241  arrayView1d< globalIndex const > const & faceDofNumber,
1242  arrayView1d< integer const > const & faceGhostRank,
1243  arrayView1d< integer const > const & isBoundaryFace,
1244  arrayView1d< real64 const > const & facePres,
1245  arrayView1d< real64 const > const & faceGravCoef,
1246  arrayView1d< real64 const > const & mimFaceGravCoef,
1247  arrayView1d< real64 const > const & transMultiplier,
1248  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob,
1249  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
1250  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens,
1251  ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseDens,
1252  ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens,
1253  ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseMassDens,
1254  ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseMassDens,
1255  ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_COMP > > const & phaseCompFrac,
1256  ElementViewConst< arrayView5d< real64 const, multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac,
1257  ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber,
1258  globalIndex const rankOffset,
1259  real64 const lengthTolerance,
1260  real64 const dt,
1261  integer const useTotalMassEquation,
1262  CRSMatrixView< real64, globalIndex const > const & localMatrix,
1263  arrayView1d< real64 > const & localRhs )
1264 {
1265  // get the cell-centered DOF numbers and ghost rank for the assembly
1266  arrayView1d< integer const > const & elemGhostRank = subRegion.ghostRank();
1267 
1268  // get the map from elem to faces
1269  arrayView2d< localIndex const > const & elemToFaces = subRegion.faceList();
1270 
1271  // get the cell-centered pressures
1272  arrayView1d< real64 const > const & elemPres =
1273  subRegion.getReference< array1d< real64 > >( fields::flow::pressure::key() );
1274 
1275  // get the element data needed for transmissibility computation
1276  arrayView2d< real64 const > const & elemCenter =
1277  subRegion.getReference< array2d< real64 > >( CellElementSubRegion::viewKeyStruct::elementCenterString() );
1278  arrayView1d< real64 const > const & elemVolume =
1279  subRegion.getReference< array1d< real64 > >( CellElementSubRegion::viewKeyStruct::elementVolumeString() );
1280 
1281  // TODO add this dependency to the compute function
1282  //arrayView3d< real64 const > const elemdPermdPres = permeabilityModel.dPerm_dPressure();
1283 
1284  arrayView3d< real64 const > const & elemPerm = permeabilityModel.permeability();
1285 
1286  // get the cell-centered depth
1287  arrayView1d< real64 const > const & elemGravCoef =
1288  subRegion.getReference< array1d< real64 > >( fields::flow::gravityCoefficient::key() );
1289 
1290  // assemble the residual and Jacobian element by element
1291  // in this loop we assemble both equation types: mass conservation in the elements and constraints at the faces
1292  forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_DEVICE ( localIndex const ei )
1293  {
1294  // transmissibility matrix
1295  stackArray2d< real64, NF *NF > transMatrix( NF, NF );
1296  stackArray2d< real64, NF *NF > transMatrixGrav( NF, NF );
1297 
1298  real64 const perm[ 3 ] = { elemPerm[ei][0][0], elemPerm[ei][0][1], elemPerm[ei][0][2] };
1299 
1300  // recompute the local transmissibility matrix at each iteration
1301  // we can decide later to precompute transMatrix if needed
1302  IP_TYPE::template compute< NF >( nodePosition,
1303  transMultiplier,
1304  faceToNodes,
1305  elemToFaces[ei],
1306  elemCenter[ei],
1307  elemVolume[ei],
1308  perm,
1309  lengthTolerance,
1310  transMatrix );
1311 
1312  // currently the gravity term in the transport scheme is treated as in MRST, that is, always with TPFA
1313  // this is why below we have to recompute the TPFA transmissibility in addition to the transmissibility matrix above
1314  // TODO: treat the gravity term with a consistent inner product
1315  mimeticInnerProduct::TPFAInnerProduct::compute< NF >( nodePosition,
1316  transMultiplier,
1317  faceToNodes,
1318  elemToFaces[ei],
1319  elemCenter[ei],
1320  elemVolume[ei],
1321  perm,
1322  lengthTolerance,
1323  transMatrixGrav );
1324 
1325  // perform flux assembly in this element
1326  compositionalMultiphaseHybridFVMKernels::AssemblerKernel::compute< NF, NC, NP >( er, esr, ei,
1327  regionFilter,
1328  elemRegionList,
1329  elemSubRegionList,
1330  elemList,
1331  faceDofNumber,
1332  faceGhostRank,
1333  isBoundaryFace,
1334  facePres,
1335  faceGravCoef,
1336  mimFaceGravCoef,
1337  elemToFaces[ei],
1338  elemPres[ei],
1339  elemGravCoef[ei],
1340  useTotalMassEquation,
1341  phaseDens,
1342  dPhaseDens,
1343  phaseMassDens,
1344  dPhaseMassDens,
1345  phaseMob,
1346  dPhaseMob,
1347  dCompFrac_dCompDens,
1348  phaseCompFrac,
1349  dPhaseCompFrac,
1350  elemDofNumber,
1351  elemGhostRank[ei],
1352  rankOffset,
1353  dt,
1354  transMatrix,
1355  transMatrixGrav,
1356  localMatrix,
1357  localRhs );
1358  } );
1359 }
1360 
1361 /******************************** DirichletFluxKernel ********************************/
1362 
1363 template< integer NF, integer NC, integer NP, typename IP_TYPE >
1364 void
1365 DirichletFluxKernel::
1366  launch( integer const numPhases,
1367  localIndex const er,
1368  localIndex const esr,
1369  CellElementSubRegion const & subRegion,
1370  constitutive::PermeabilityBase const & permeabilityModel,
1371  SortedArrayView< localIndex const > const & boundaryFaceSet,
1372  arrayView1d< real64 const > const & facePres,
1373  arrayView1d< real64 const > const & faceTemp,
1374  arrayView2d< real64 const, compflow::USD_PHASE > const & facePhaseMob,
1375  arrayView2d< real64 const, compflow::USD_PHASE > const & facePhaseMassDens,
1376  arrayView3d< real64 const, compflow::USD_PHASE_COMP > const & facePhaseCompFrac,
1377  arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition,
1378  ArrayOfArraysView< localIndex const > const & faceToNodes,
1379  arrayView2d< localIndex const > const & elemRegionList,
1380  arrayView2d< localIndex const > const & elemSubRegionList,
1381  arrayView2d< localIndex const > const & elemList,
1382  arrayView1d< globalIndex const > const & faceDofNumber,
1383  arrayView1d< real64 const > const & faceGravCoef,
1384  arrayView1d< real64 const > const & transMultiplier,
1385  real64 const lengthTolerance,
1386  real64 const dt,
1387  globalIndex const rankOffset,
1388  integer const useTotalMassEquation,
1389  CompFlowAccessors const & compFlowAccessors,
1390  MultiFluidAccessors const & multiFluidAccessors,
1391  ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber,
1392  CRSMatrixView< real64, globalIndex const > const & localMatrix,
1393  arrayView1d< real64 > const & localRhs )
1394 {
1395  GEOS_UNUSED_VAR( faceTemp );
1396  using Deriv = multifluid::DerivativeOffset;
1397 
1398  // Get element data
1399  arrayView2d< localIndex const > const & elemToFaces = subRegion.faceList();
1400  arrayView1d< real64 const > const & elemPres = subRegion.getField< fields::flow::pressure >();
1401  arrayView1d< real64 const > const & elemGravCoef = subRegion.getField< fields::flow::gravityCoefficient >();
1402  arrayView2d< real64 const > const & elemCenter = subRegion.getReference< array2d< real64 > >( CellElementSubRegion::viewKeyStruct::elementCenterString() );
1403  arrayView1d< real64 const > const & elemVolume = subRegion.getReference< array1d< real64 > >( CellElementSubRegion::viewKeyStruct::elementVolumeString() );
1404  arrayView3d< real64 const > const & elemPerm = permeabilityModel.permeability();
1405  arrayView1d< integer const > const & elemGhostRank = subRegion.ghostRank();
1406 
1407  // Get compositional flow fields
1408  auto phaseMob = compFlowAccessors.get( fields::flow::phaseMobility{} );
1409  auto dPhaseMob = compFlowAccessors.get( fields::flow::dPhaseMobility{} );
1410  auto dCompFrac_dCompDens = compFlowAccessors.get( fields::flow::dGlobalCompFraction_dGlobalCompDensity{} );
1411 
1412  // Get multifluid fields
1413  auto phaseMassDens = multiFluidAccessors.get( fields::multifluid::phaseMassDensity{} );
1414  auto dPhaseMassDens = multiFluidAccessors.get( fields::multifluid::dPhaseMassDensity{} );
1415  auto phaseCompFrac = multiFluidAccessors.get( fields::multifluid::phaseCompFraction{} );
1416  auto dPhaseCompFrac = multiFluidAccessors.get( fields::multifluid::dPhaseCompFraction{} );
1417 
1418  // Loop over boundary faces
1419  forAll< parallelDevicePolicy<> >( boundaryFaceSet.size(), [=] GEOS_DEVICE ( localIndex const iset )
1420  {
1421  localIndex const kf = boundaryFaceSet[iset];
1422 
1423  // Find the element adjacent to this boundary face
1424  localIndex erAdj = -1, esrAdj = -1, eiAdj = -1;
1425  for( integer ke = 0; ke < elemRegionList.size( 1 ); ++ke )
1426  {
1427  if( elemRegionList[kf][ke] == er && elemSubRegionList[kf][ke] == esr )
1428  {
1429  erAdj = elemRegionList[kf][ke];
1430  esrAdj = elemSubRegionList[kf][ke];
1431  eiAdj = elemList[kf][ke];
1432  break;
1433  }
1434  }
1435 
1436  // Skip if no adjacent element in target region
1437  if( eiAdj < 0 || elemGhostRank[eiAdj] >= 0 )
1438  {
1439  return;
1440  }
1441 
1442  // Compute one-sided transmissibility
1443  stackArray2d< real64, NF * NF > transMatrix( NF, NF );
1444  real64 const perm[3] = { elemPerm[eiAdj][0][0], elemPerm[eiAdj][0][1], elemPerm[eiAdj][0][2] };
1445 
1446  IP_TYPE::template compute< NF >( nodePosition,
1447  transMultiplier,
1448  faceToNodes,
1449  elemToFaces[eiAdj],
1450  elemCenter[eiAdj],
1451  elemVolume[eiAdj],
1452  perm,
1453  lengthTolerance,
1454  transMatrix );
1455 
1456  // Find local face index
1457  integer ifaceLoc = -1;
1458  for( integer j = 0; j < NF; ++j )
1459  {
1460  if( elemToFaces[eiAdj][j] == kf )
1461  {
1462  ifaceLoc = j;
1463  break;
1464  }
1465  }
1466  if( ifaceLoc < 0 )
1467  {
1468  return;
1469  }
1470 
1471  // Get all global face indices of the cell touching the boundary
1472  stackArray1d< localIndex, NF > cellFaces( NF );
1473  for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc )
1474  {
1475  cellFaces[jfaceLoc] = elemToFaces[eiAdj][jfaceLoc];
1476  }
1477 
1478  // Component fluxes
1479  real64 compFlux[NC]{};
1480  real64 dCompFlux_dP[NC]{};
1481  real64 dCompFlux_dC[NC][NC]{};
1482  real64 dCompFlux_dFaceP[NC][NF]{}; // Derivatives w.r.t. face pressures
1483 
1484  // Loop over phases to compute component fluxes
1485  // Use face-based properties for boundary conditions (upwinding)
1486  for( integer ip = 0; ip < numPhases; ++ip )
1487  {
1488  real64 dDensMean_dC[NC]{};
1489  real64 dF_dC[NC]{};
1490  real64 dProp_dC[NC]{};
1491 
1492  // Use average of element and face phase mass density for gravity term
1493  real64 const elemDens = phaseMassDens[erAdj][esrAdj][eiAdj][0][ip];
1494  real64 const faceDens = facePhaseMassDens[kf][ip];
1495  real64 const densMean = 0.5 * (elemDens + faceDens);
1496 
1497  applyChainRule( NC,
1498  dCompFrac_dCompDens[erAdj][esrAdj][eiAdj],
1499  dPhaseMassDens[erAdj][esrAdj][eiAdj][0][ip],
1500  dProp_dC,
1501  Deriv::dC );
1502 
1503  real64 const dDensMean_dP = 0.5 * dPhaseMassDens[erAdj][esrAdj][eiAdj][0][ip][Deriv::dP];
1504  for( integer jc = 0; jc < NC; ++jc )
1505  {
1506  dDensMean_dC[jc] = 0.5 * dProp_dC[jc];
1507  }
1508 
1509  // Compute flux by looping over all connected faces with consistent transmissibility
1510  real64 f = 0.0;
1511  real64 dF_dP = 0.0;
1512  real64 dF_dFaceP[NF]{}; // Derivatives of flux w.r.t. face pressures
1513  for( integer jc = 0; jc < NC; ++jc )
1514  {
1515  dF_dC[jc] = 0.0;
1516  }
1517 
1518  // Sum contributions from all connected faces
1519  for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc )
1520  {
1521  // Each stencil face j uses its own face pressure
1522  localIndex const kfj = cellFaces[jfaceLoc];
1523  real64 const gravTimesDz_j = elemGravCoef[eiAdj] - faceGravCoef[kfj];
1524  real64 const Tij = transMatrix[ifaceLoc][jfaceLoc];
1525  real64 const potDif_j = elemPres[eiAdj] - facePres[kfj] - densMean * gravTimesDz_j;
1526 
1527  // Accumulate flux and derivatives
1528  f += Tij * potDif_j;
1529  dF_dP += Tij * ( 1.0 - dDensMean_dP * gravTimesDz_j );
1530 
1531  // Correct Jacobian: d(potDif_j)/d(facePres[kfj]) = -1, so dF/d(facePres_j) = -T_ij
1532  dF_dFaceP[jfaceLoc] = -Tij;
1533 
1534  for( integer jc = 0; jc < NC; ++jc )
1535  {
1536  dF_dC[jc] += -Tij * dDensMean_dC[jc] * gravTimesDz_j;
1537  }
1538  }
1539 
1540  // Use element mobility (simplified upwinding for Dirichlet BC)
1541  // Upwind phase component fraction based on flow direction
1542  // Use the total flux f for upwinding decision
1543  real64 const beta = (f > 0.0) ? 1.0 : 0.0;
1544  real64 const facePhaseMobility = facePhaseMob[kf][ip];
1545  real64 const phaseMobility = beta * phaseMob[erAdj][esrAdj][eiAdj][ip] + (1.0 - beta) * facePhaseMobility;
1546  real64 const phaseFlux = phaseMobility * f;
1547  real64 const dPhaseFlux_dP = phaseMobility * dF_dP + beta * dPhaseMob[erAdj][esrAdj][eiAdj][ip][Deriv::dP] * f;
1548  real64 dPhaseFlux_dC[NC];
1549  real64 dPhaseFlux_dFaceP[NF];
1550  for( integer jc = 0; jc < NC; ++jc )
1551  {
1552  dPhaseFlux_dC[jc] = phaseMobility * dF_dC[jc] + beta * dPhaseMob[erAdj][esrAdj][eiAdj][ip][Deriv::dC+jc] * f;
1553  }
1554  for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc )
1555  {
1556  dPhaseFlux_dFaceP[jfaceLoc] = phaseMobility * dF_dFaceP[jfaceLoc];
1557  }
1558 
1559  // Use face-based phase composition for boundary conditions when flow is INTO the domain (potDif < 0)
1560  // Use element phase composition when flow is OUT of the domain (potDif > 0)
1561  for( integer ic = 0; ic < NC; ++ic )
1562  {
1563  real64 const ycpElem = phaseCompFrac[erAdj][esrAdj][eiAdj][0][ip][ic];
1564  real64 const ycpFace = facePhaseCompFrac[kf][ip][ic];
1565 
1566  // Upwind phase component fraction based on flow direction
1567  real64 const ycp = beta * ycpElem + (1.0 - beta) * ycpFace;
1568 
1569  compFlux[ic] += ycp * phaseFlux;
1570 
1571  // For derivatives, use element properties (simplified for boundary conditions)
1572  dCompFlux_dP[ic] += dPhaseFlux_dP * ycp + beta * phaseFlux * dPhaseCompFrac[erAdj][esrAdj][eiAdj][0][ip][ic][Deriv::dP];
1573 
1574  // Compute dycpElem/dC for this component ic
1575  real64 dycpElem_dC[NC];
1576  applyChainRule( NC,
1577  dCompFrac_dCompDens[erAdj][esrAdj][eiAdj],
1578  dPhaseCompFrac[erAdj][esrAdj][eiAdj][0][ip][ic],
1579  dycpElem_dC,
1580  Deriv::dC );
1581  for( integer jc = 0; jc < NC; ++jc )
1582  {
1583  // d(ycp * phaseFlux)/dC = dycp/dC * phaseFlux + ycp * dPhaseFlux/dC
1584  // dycp/dC = beta * dycpElem/dC (ycpFace is BC, independent of C)
1585  dCompFlux_dC[ic][jc] += ycp * dPhaseFlux_dC[jc] + beta * phaseFlux * dycpElem_dC[jc];
1586  }
1587 
1588  // Add derivatives w.r.t. face pressures
1589  for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc )
1590  {
1591  dCompFlux_dFaceP[ic][jfaceLoc] += dPhaseFlux_dFaceP[jfaceLoc] * ycp;
1592  }
1593  }
1594  }
1595 
1596  // Assemble into residual and Jacobian
1597  real64 localFlux[NC];
1598  real64 localFluxJacobian[NC][NC+1];
1599 
1600  for( integer ic = 0; ic < NC; ++ic )
1601  {
1602  localFlux[ic] = dt * compFlux[ic];
1603  localFluxJacobian[ic][0] = dt * dCompFlux_dP[ic];
1604  for( integer jc = 0; jc < NC; ++jc )
1605  {
1606  localFluxJacobian[ic][jc+1] = dt * dCompFlux_dC[ic][jc];
1607  }
1608  }
1609 
1610  // Apply total mass equation transformation if needed
1611  if( useTotalMassEquation )
1612  {
1613  real64 work[NC+1]{};
1614  compositionalMultiphaseUtilities::shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NC+1, localFluxJacobian, work );
1615  compositionalMultiphaseUtilities::shiftElementsAheadByOneAndReplaceFirstElementWithSum( NC, localFlux );
1616  }
1617 
1618  // Add to global system
1619  globalIndex const elemDof = elemDofNumber[erAdj][esrAdj][eiAdj];
1620  localIndex const localRow = LvArray::integerConversion< localIndex >( elemDof - rankOffset );
1621 
1622  for( integer ic = 0; ic < NC; ++ic )
1623  {
1624  RAJA::atomicAdd( parallelDeviceAtomic{}, &localRhs[localRow + ic], localFlux[ic] );
1625 
1626  globalIndex dofColIndices[NC+1];
1627  dofColIndices[0] = elemDof;
1628  for( integer jc = 0; jc < NC; ++jc )
1629  {
1630  dofColIndices[jc+1] = elemDof + jc + 1;
1631  }
1632 
1633  localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( localRow + ic,
1634  dofColIndices,
1635  localFluxJacobian[ic],
1636  NC+1 );
1637 
1638  // Add contributions from face pressure derivatives
1639  for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc )
1640  {
1641  localIndex const kfj = cellFaces[jfaceLoc];
1642  globalIndex const faceDof = faceDofNumber[kfj];
1643 
1644  real64 facePressureJacobian = dt * dCompFlux_dFaceP[ic][jfaceLoc];
1645 
1646  // Apply total mass equation transformation if needed
1647  if( useTotalMassEquation && ic == 0 )
1648  {
1649  // For total mass equation, sum all component derivatives
1650  facePressureJacobian = 0.0;
1651  for( integer icc = 0; icc < NC; ++icc )
1652  {
1653  facePressureJacobian += dt * dCompFlux_dFaceP[icc][jfaceLoc];
1654  }
1655  }
1656 
1657  localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( localRow + ic,
1658  &faceDof,
1659  &facePressureJacobian,
1660  1 );
1661  }
1662 
1663  }
1664 
1665  } );
1666 }
1667 
1668 /******************************** EvaluateBCFacePropertiesKernel ********************************/
1669 
1694 template< integer NC, integer NP >
1695 void
1696 evaluateBCFaceProperties( integer const numPhases,
1697  SortedArrayView< localIndex const > const & boundaryFaceSet,
1698  arrayView1d< real64 const > const & facePres,
1699  arrayView1d< real64 const > const & faceTemp,
1700  arrayView2d< real64 const, compflow::USD_COMP > const & faceCompFrac,
1701  arrayView2d< localIndex const > const & elemRegionList,
1702  arrayView2d< localIndex const > const & elemSubRegionList,
1703  arrayView2d< localIndex const > const & elemList,
1704  localIndex const er,
1705  localIndex const esr,
1706  constitutive::MultiFluidBase & fluid,
1707  constitutive::RelativePermeabilityBase & relperm,
1708  arrayView2d< real64, compflow::USD_PHASE > const & facePhaseMob,
1709  arrayView2d< real64, compflow::USD_PHASE > const & facePhaseMassDens,
1710  arrayView3d< real64, compflow::USD_PHASE_COMP > const & facePhaseCompFrac )
1711 {
1712  // Use constitutiveUpdatePassThru to dispatch to concrete fluid and relperm types
1713  constitutiveUpdatePassThru( fluid, [&] ( auto & castedFluid )
1714  {
1715  using FluidType = TYPEOFREF( castedFluid );
1716  typename FluidType::KernelWrapper fluidWrapper = castedFluid.createKernelWrapper();
1717 
1718  constitutive::constitutiveUpdatePassThru( relperm, [&] ( auto & castedRelperm )
1719  {
1720  using RelPermType = TYPEOFREF( castedRelperm );
1721  typename RelPermType::KernelWrapper relPermWrapper = castedRelperm.createKernelWrapper();
1722  GEOS_UNUSED_VAR( relPermWrapper );
1723 
1724  // Loop over BC faces and evaluate properties at BC conditions
1725  // Note: Using serialPolicy (host) because extended device lambdas cannot be defined inside generic lambdas
1726  // The output arrays will be moved to device after this completes
1727  forAll< serialPolicy >( boundaryFaceSet.size(), [=, &facePhaseMob, &facePhaseMassDens, &facePhaseCompFrac] ( localIndex const iset )
1728  {
1729  localIndex const kf = boundaryFaceSet[iset];
1730 
1731  // Find adjacent element in target region
1732  localIndex eiAdj = -1;
1733  for( integer ke = 0; ke < elemRegionList.size( 1 ); ++ke )
1734  {
1735  if( elemRegionList[kf][ke] == er && elemSubRegionList[kf][ke] == esr )
1736  {
1737  eiAdj = elemList[kf][ke];
1738  break;
1739  }
1740  }
1741  if( eiAdj < 0 )
1742  return;
1743 
1744  // Allocate temporary storage for face constitutive properties
1745  StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseFrac( 1, 1, numPhases );
1746  StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseDens( 1, 1, numPhases );
1747  StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseMassDensLocal( 1, 1, numPhases );
1748  StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseVisc( 1, 1, numPhases );
1749  StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseEnthalpy( 1, 1, numPhases );
1750  StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseInternalEnergy( 1, 1, numPhases );
1751  StackArray< real64, 4, constitutive::MultiFluidBase::MAX_NUM_PHASES * NC,
1752  constitutive::multifluid::LAYOUT_PHASE_COMP > facePhaseCompFracLocal( 1, 1, numPhases, NC );
1753 
1754  // Initialize all temporary arrays to zero to ensure deterministic behavior on GPU
1755  for( integer ip = 0; ip < numPhases; ++ip )
1756  {
1757  facePhaseFrac[0][0][ip] = 0.0;
1758  facePhaseDens[0][0][ip] = 0.0;
1759  facePhaseMassDensLocal[0][0][ip] = 0.0;
1760  facePhaseVisc[0][0][ip] = 0.0;
1761  facePhaseEnthalpy[0][0][ip] = 0.0;
1762  facePhaseInternalEnergy[0][0][ip] = 0.0;
1763  for( integer ic = 0; ic < NC; ++ic )
1764  {
1765  facePhaseCompFracLocal[0][0][ip][ic] = 0.0;
1766  }
1767  }
1768 
1769  real64 faceTotalDens = 0.0;
1770 
1771  // Evaluate fluid properties at BC face conditions using flash calculation
1772  constitutive::MultiFluidBase::KernelWrapper::computeValues( fluidWrapper,
1773  facePres[kf],
1774  faceTemp[kf],
1775  faceCompFrac[kf],
1776  facePhaseFrac[0][0],
1777  facePhaseDens[0][0],
1778  facePhaseMassDensLocal[0][0],
1779  facePhaseVisc[0][0],
1780  facePhaseEnthalpy[0][0],
1781  facePhaseInternalEnergy[0][0],
1782  facePhaseCompFracLocal[0][0],
1783  faceTotalDens );
1784 
1785  // Evaluate relative permeability at face saturation from flash calculation
1786 // relPermWrapper.compute( facePhaseFrac[0][0], 0, 0 );
1787 
1788  // Store computed properties in output arrays
1789  for( integer ip = 0; ip < numPhases; ++ip )
1790  {
1791  // Store phase mass density from flash calculation
1792  facePhaseMassDens[kf][ip] = facePhaseMassDensLocal[0][0][ip];
1793 
1794  // Compute mobility from relative permeability evaluated at face conditions
1795  real64 const faceKr = facePhaseFrac[0][0][ip]; // phaseRelPerm[eiAdj][0][ip];
1796  real64 const mu = facePhaseVisc[0][0][ip];
1797  // Safety check: ensure faceTotalDens and mu are valid before computing mobility
1798  facePhaseMob[kf][ip] = (mu > 0 && faceTotalDens > 0) ? faceTotalDens * faceKr / mu : 0.0;
1799 
1800  // Store phase composition from flash calculation
1801  for( integer ic = 0; ic < NC; ++ic )
1802  {
1803  facePhaseCompFrac[kf][ip][ic] = facePhaseCompFracLocal[0][0][ip][ic];
1804  }
1805  }
1806  } );
1807  } ); // end nested constitutiveUpdatePassThru for relperm
1808  } ); // end constitutiveUpdatePassThru for fluid
1809 }
1810 
1811 } // namespace compositionalMultiphaseHybridFVMKernels
1812 
1813 } // namespace geos
geos::compositionalMultiphaseHybridFVMKernels::evaluateBCFaceProperties
void evaluateBCFaceProperties(integer const numPhases, SortedArrayView< localIndex const > const &boundaryFaceSet, arrayView1d< real64 const > const &facePres, arrayView1d< real64 const > const &faceTemp, arrayView2d< real64 const, compflow::USD_COMP > const &faceCompFrac, arrayView2d< localIndex const > const &elemRegionList, arrayView2d< localIndex const > const &elemSubRegionList, arrayView2d< localIndex const > const &elemList, localIndex const er, localIndex const esr, constitutive::MultiFluidBase &fluid, constitutive::RelativePermeabilityBase &relperm, arrayView2d< real64, compflow::USD_PHASE > const &facePhaseMob, arrayView2d< real64, compflow::USD_PHASE > const &facePhaseMassDens, arrayView3d< real64, compflow::USD_PHASE_COMP > const &facePhaseCompFrac)
Evaluate constitutive properties at BC face conditions.
Definition: CompositionalMultiphaseHybridFVMKernels_impl.hpp:1696
GEOS_HOST_DEVICE
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
GEOS_UNUSED_VAR
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
Definition: GeosxMacros.hpp:84
GEOS_DEVICE
#define GEOS_DEVICE
Marks a device-only function.
Definition: GeosxMacros.hpp:47
GEOS_ASSERT_GT
#define GEOS_ASSERT_GT(lhs, rhs)
Assert that one value compares greater than the other in debug builds.
Definition: Logger.hpp:853
GEOS_ASSERT_GE
#define GEOS_ASSERT_GE(lhs, rhs)
Assert that one value compares greater than or equal to the other in debug builds.
Definition: Logger.hpp:870
geos::CellElementSubRegion::faceList
FixedOneToManyRelation & faceList()
Get the element-to-face map.
Definition: CellElementSubRegion.hpp:225
geos::ObjectManagerBase::ghostRank
array1d< integer > const & ghostRank()
Get the ghost information of each object.
Definition: ObjectManagerBase.hpp:856
geos::ObjectManagerBase::getField
GEOS_DECLTYPE_AUTO_RETURN getField() const
Get a view to the field associated with a trait from this ObjectManagerBase.
Definition: ObjectManagerBase.hpp:629
geos::StencilAccessors
A struct to automatically construct and store element view accessors.
Definition: StencilAccessors.hpp:39
geos::StencilMaterialAccessors
A struct to automatically construct and store element view accessors.
Definition: StencilAccessors.hpp:98
geos::dataRepository::Group::getReference
GEOS_DECLTYPE_AUTO_RETURN getReference(LOOKUP_TYPE const &lookup) const
Look up a wrapper and get reference to wrapped object.
Definition: Group.hpp:1275
geos::arrayView1d
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:179
geos::stackArray2d
StackArray< T, 2, MAXSIZE > stackArray2d
Alias for 2D stack array.
Definition: DataTypes.hpp:203
geos::array2d
Array< T, 2, PERMUTATION > array2d
Alias for 2D array.
Definition: DataTypes.hpp:191
geos::globalIndex
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
Definition: DataTypes.hpp:87
geos::StackArray
LvArray::StackArray< T, NDIM, PERMUTATION, localIndex, MAXSIZE > StackArray
Multidimensional stack-based array type. See LvArray:StackArray for details.
Definition: DataTypes.hpp:155
geos::stackArray1d
StackArray< T, 1, MAXSIZE > stackArray1d
Alias for 1D stack array.
Definition: DataTypes.hpp:187
geos::ArrayOfArraysView
LvArray::ArrayOfArraysView< T, INDEX_TYPE const, CONST_SIZES, LvArray::ChaiBuffer > ArrayOfArraysView
View of array of variable-sized arrays. See LvArray::ArrayOfArraysView for details.
Definition: DataTypes.hpp:285
geos::arraySlice2d
ArraySlice< T, 2, USD > arraySlice2d
Alias for 2D array slice.
Definition: DataTypes.hpp:199
geos::arrayView5d
ArrayView< T, 5, USD > arrayView5d
Alias for 5D array view.
Definition: DataTypes.hpp:243
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::arraySlice1d
ArraySlice< T, 1, USD > arraySlice1d
Alias for 1D array slice.
Definition: DataTypes.hpp:183
geos::CRSMatrixView
LvArray::CRSMatrixView< T, COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
Definition: DataTypes.hpp:309
geos::SortedArrayView
LvArray::SortedArrayView< T, localIndex, LvArray::ChaiBuffer > SortedArrayView
A sorted array view of local indices.
Definition: DataTypes.hpp:270
geos::arrayView4d
ArrayView< T, 4, USD > arrayView4d
Alias for 4D array view.
Definition: DataTypes.hpp:227
geos::arrayView2d
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:195
geos::integer
int integer
Signed integer type.
Definition: DataTypes.hpp:81
geos::array1d
Array< T, 1 > array1d
Alias for 1D array.
Definition: DataTypes.hpp:175
geos::arrayView3d
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
Definition: DataTypes.hpp:211
geos::ElementSubRegionBase::viewKeyStruct::elementVolumeString
static constexpr char const * elementVolumeString()
Definition: ElementSubRegionBase.hpp:231
geos::ElementSubRegionBase::viewKeyStruct::elementCenterString
static constexpr char const * elementCenterString()
Definition: ElementSubRegionBase.hpp:229
geos::compositionalMultiphaseHybridFVMKernels::AssemblerKernelHelper::ElementViewConst
ElementRegionManager::ElementViewConst< VIEWTYPE > ElementViewConst
The type for element-based non-constitutive data parameters. Consists entirely of ArrayView's.
Definition: CompositionalMultiphaseHybridFVMKernels.hpp:268
geos::compositionalMultiphaseHybridFVMKernels::DirichletFluxKernel::ElementViewConst
ElementRegionManager::ElementViewConst< VIEWTYPE > ElementViewConst
The type for element-based non-constitutive data parameters. Consists entirely of ArrayView's.
Definition: CompositionalMultiphaseHybridFVMKernels.hpp:567
geos::compositionalMultiphaseHybridFVMKernels::UpwindingHelper::ElementViewConst
ElementRegionManager::ElementViewConst< VIEWTYPE > ElementViewConst
The type for element-based non-constitutive data parameters. Consists entirely of ArrayView's.
Definition: CompositionalMultiphaseHybridFVMKernels.hpp:66
Generated by doxygen 1.9.1