GEOS
PerforationFluxKernels.hpp
Go to the documentation of this file.
1 /*
2  * ------------------------------------------------------------------------------------------------------------
3  * SPDX-License-Identifier: LGPL-2.1-only
4  *
5  * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
6  * Copyright (c) 2018-2024 TotalEnergies
7  * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
8  * Copyright (c) 2023-2024 Chevron
9  * Copyright (c) 2019- GEOS/GEOSX Contributors
10  * All rights reserved
11  *
12  * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
13  * ------------------------------------------------------------------------------------------------------------
14  */
15 
20 #ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_PERFORATIONFLUXLKERNELS_HPP
21 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_PERFORATIONFLUXLKERNELS_HPP
22 
23 #include "codingUtilities/Utilities.hpp"
24 #include "common/DataTypes.hpp"
25 #include "common/GEOS_RAJA_Interface.hpp"
26 #include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
27 #include "constitutive/fluid/multifluid/MultiFluidFields.hpp"
28 #include "constitutive/relativePermeability/RelativePermeabilityBase.hpp"
29 #include "constitutive/relativePermeability/RelativePermeabilityFields.hpp"
30 #include "mesh/ElementRegionManager.hpp"
31 #include "mesh/ObjectManagerBase.hpp"
32 #include "physicsSolvers/KernelLaunchSelectors.hpp"
33 #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp"
34 #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp"
35 #include "physicsSolvers/fluidFlow/StencilAccessors.hpp"
36 #include "physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellFields.hpp"
37 #include "physicsSolvers/fluidFlow/wells/WellSolverBaseFields.hpp"
38 #include "physicsSolvers/fluidFlow/wells/WellTags.hpp"
39 #include "physicsSolvers/fluidFlow/wells/WellFields.hpp"
40 
41 
42 namespace geos
43 {
44 
45 struct NoOpStuct
46 {
47  NoOpStuct(){}
48 };
49 
50 namespace isothermalPerforationFluxKernels
51 {
52 
53 /******************************** PerforationFluxKernel ********************************/
54 
55 template< integer NC, integer NP, integer IS_THERMAL >
56 class PerforationFluxKernel
57 {
58 public:
60  static constexpr integer numComp = NC;
61 
63  static constexpr integer numPhase = NP;
64 
66  static constexpr integer isThermal = IS_THERMAL;
67 
68  using TAG = wellTags::SubRegionTag;
69 
70  using CompFlowAccessors =
71  StencilAccessors< fields::flow::pressure,
72  fields::flow::phaseVolumeFraction,
73  fields::flow::dPhaseVolumeFraction,
74  fields::flow::dGlobalCompFraction_dGlobalCompDensity >;
75 
76  using MultiFluidAccessors =
77  StencilMaterialAccessors< constitutive::MultiFluidBase,
78  fields::multifluid::phaseDensity,
79  fields::multifluid::dPhaseDensity,
80  fields::multifluid::phaseViscosity,
81  fields::multifluid::dPhaseViscosity,
82  fields::multifluid::phaseCompFraction,
83  fields::multifluid::dPhaseCompFraction >;
84 
85  using RelPermAccessors =
86  StencilMaterialAccessors< constitutive::RelativePermeabilityBase,
87  fields::relperm::phaseRelPerm,
88  fields::relperm::dPhaseRelPerm_dPhaseVolFraction >;
89 
90 
98  template< typename VIEWTYPE >
99  using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
100 
101  PerforationFluxKernel ( PerforationData * const perforationData,
102  ElementSubRegionBase const & subRegion,
103  CompFlowAccessors const & compFlowAccessors,
104  MultiFluidAccessors const & multiFluidAccessors,
105  RelPermAccessors const & relPermAccessors,
106  bool const isInjector,
107  bool const isCrossflowEnabled ):
108  m_resPres( compFlowAccessors.get( fields::flow::pressure {} )),
109  m_resPhaseVolFrac( compFlowAccessors.get( fields::flow::phaseVolumeFraction {} )),
110  m_dResPhaseVolFrac( compFlowAccessors.get( fields::flow::dPhaseVolumeFraction {} )),
111  m_dResCompFrac_dCompDens( compFlowAccessors.get( fields::flow::dGlobalCompFraction_dGlobalCompDensity {} )),
112  m_resPhaseDens( multiFluidAccessors.get( fields::multifluid::phaseDensity {} )),
113  m_dResPhaseDens( multiFluidAccessors.get( fields::multifluid::dPhaseDensity {} )),
114  m_resPhaseVisc( multiFluidAccessors.get( fields::multifluid::phaseViscosity {} )),
115  m_dResPhaseVisc( multiFluidAccessors.get( fields::multifluid::dPhaseViscosity {} )),
116  m_resPhaseCompFrac( multiFluidAccessors.get( fields::multifluid::phaseCompFraction {} )),
117  m_dResPhaseCompFrac( multiFluidAccessors.get( fields::multifluid::dPhaseCompFraction {} )),
118  m_resPhaseRelPerm( relPermAccessors.get( fields::relperm::phaseRelPerm {} )),
119  m_dResPhaseRelPerm_dPhaseVolFrac( relPermAccessors.get( fields::relperm::dPhaseRelPerm_dPhaseVolFraction {} )),
120  m_wellElemGravCoef( subRegion.getField< fields::well::gravityCoefficient >()),
121  m_wellElemPres( subRegion.getField< fields::well::pressure >()),
122  m_wellElemCompDens( subRegion.getField< fields::well::globalCompDensity >()),
123  m_wellElemTotalMassDens( subRegion.getField< fields::well::totalMassDensity >()),
124  m_dWellElemTotalMassDens( subRegion.getField< fields::well::dTotalMassDensity >()),
125  m_wellElemCompFrac( subRegion.getField< fields::well::globalCompFraction >()),
126  m_dWellElemCompFrac_dCompDens( subRegion.getField< fields::well::dGlobalCompFraction_dGlobalCompDensity >()),
127  m_perfGravCoef( perforationData->getField< fields::well::gravityCoefficient >()),
128  m_perfWellElemIndex( perforationData->getField< fields::perforation::wellElementIndex >()),
129  m_perfTrans( perforationData->getField< fields::perforation::wellTransmissibility >()),
130  m_resElementRegion( perforationData->getField< fields::perforation::reservoirElementRegion >()),
131  m_resElementSubRegion( perforationData->getField< fields::perforation::reservoirElementSubRegion >()),
132  m_resElementIndex( perforationData->getField< fields::perforation::reservoirElementIndex >()),
133  m_compPerfRate( perforationData->getField< fields::well::compPerforationRate >()),
134  m_dCompPerfRate( perforationData->getField< fields::well::dCompPerforationRate >()),
135  m_perfStatus( perforationData->getField< fields::perforation::perforationStatus >()),
136  m_isInjector( isInjector ),
137  m_isCrossflowEnabled( isCrossflowEnabled )
138  {}
139 
140  template< typename FUNC = NoOpFunc >
141  GEOS_HOST_DEVICE
142  inline
143  void
144  computeFlux( localIndex const iperf, FUNC && fluxKernelOp= NoOpFunc {} ) const
145  {
146  using CP_Deriv = constitutive::multifluid::DerivativeOffsetC< NC, IS_THERMAL >;
147 
148  if( !m_perfStatus[iperf] )
149  {
150  for( integer ic = 0; ic < NC; ++ic )
151  {
152  m_compPerfRate[iperf][ic] = 0.0;
153  for( integer ke = 0; ke < 2; ++ke )
154  {
155  for( integer jc = 0; jc < CP_Deriv::nDer; ++jc )
156  {
157  m_dCompPerfRate[iperf][ke][ic][jc] = 0.0;
158  }
159  }
160  }
161  return;
162  }
163  // get the index of the reservoir elem
164  localIndex const er = m_resElementRegion[iperf];
165  localIndex const esr = m_resElementSubRegion[iperf];
166  localIndex const ei = m_resElementIndex[iperf];
167 
168  // get the index of the well elem
169  localIndex const iwelem = m_perfWellElemIndex[iperf];
170 
171  using Deriv = constitutive::multifluid::DerivativeOffset;
172 
173 
174  // local working variables and arrays
175  real64 pres[2]{};
176  real64 multiplier[2]{};
177 
178  // local working variables - compact
179  // All derivative quantiites generated are stored in arrays using CP_Deriv offsets
180  // The input well/reservoir quantites use the Deriv offsets
181  // The arrays using the deriv offsets have extra column for dT in isothermal cases
182 
183  real64 dPres[2][CP_Deriv::nDer]{};
184  real64 dFlux[2][CP_Deriv::nDer]{};
185  real64 dMob[CP_Deriv::nDer]{};
186  real64 dPotDiff[2][CP_Deriv::nDer]{};
187  real64 dCompFrac[CP_Deriv::nDer]{};
188 
189  // Step 1: reset the perforation rates
190  for( integer ic = 0; ic < NC; ++ic )
191  {
192  m_compPerfRate[iperf][ic] = 0.0;
193  for( integer ke = 0; ke < 2; ++ke )
194  {
195  for( integer jc = 0; jc < CP_Deriv::nDer; ++jc )
196  {
197  m_dCompPerfRate[iperf][ke][ic][jc] = 0.0;
198  }
199  }
200  }
201 
202  // Step 2: copy the variables from the reservoir and well element
203 
204  // a) get reservoir variables
205 
206  pres[TAG::RES] = m_resPres[er][esr][ei];
207  dPres[TAG::RES][CP_Deriv::dP] = 1.0;
208  multiplier[TAG::RES] = 1.0;
209 
210  // Here in the absence of a buoyancy term we assume that the reservoir cell is perforated at its center
211  // TODO: add a buoyancy term for the reservoir side here
212 
213 
214  // b) get well variables
215 
216  pres[TAG::WELL] = m_wellElemPres[iwelem];
217  dPres[TAG::WELL][CP_Deriv::dP] = 1.0;
218  multiplier[TAG::WELL] = -1.0;
219 
220  real64 const gravD = ( m_perfGravCoef[iperf] - m_wellElemGravCoef[iwelem] );
221 
222  pres[TAG::WELL] += m_wellElemTotalMassDens[iwelem] * gravD;
223  // Note LHS uses CP_Deriv while RHS uses Deriv !!!
224  dPres[TAG::WELL][CP_Deriv::dP] += m_dWellElemTotalMassDens[iwelem][Deriv::dP] * gravD;
225  if constexpr ( IS_THERMAL )
226  {
227  dPres[TAG::WELL][CP_Deriv::dT] += m_dWellElemTotalMassDens[iwelem][Deriv::dT] * gravD;
228  }
229  for( integer ic = 0; ic < NC; ++ic )
230  {
231  dPres[TAG::WELL][CP_Deriv::dC+ic] += m_dWellElemTotalMassDens[iwelem][Deriv::dC+ic] * gravD;
232  }
233 
234  // Step 3: compute potential difference
235 
236  real64 potDiff = 0.0;
237  for( integer i = 0; i < 2; ++i )
238  {
239  potDiff += multiplier[i] * m_perfTrans[iperf] * pres[i];
240  // LHS & RHS both use CP_Deriv
241  for( integer ic = 0; ic < CP_Deriv::nDer; ++ic )
242  {
243  dPotDiff[i][ic] += multiplier[i] * m_perfTrans[iperf] * dPres[i][ic];
244  }
245  }
246  // Step 4: upwinding based on the flow direction
247 
248  real64 flux = 0.0;
249  if( potDiff >= 0 ) // ** reservoir cell is upstream **
250  {
251 
252  // loop over phases, compute and upwind phase flux
253  // and sum contributions to each component's perforation rate
254  for( integer ip = 0; ip < NP; ++ip )
255  {
256  // skip the rest of the calculation if the phase is absent
257  // or if crossflow is disabled for injectors
258  bool const phaseExists = (m_resPhaseVolFrac[er][esr][ei][ip] > 0);
259  if( !phaseExists || (m_isInjector && !m_isCrossflowEnabled) )
260  {
261  continue;
262  }
263 
264  // here, we have to recompute the reservoir phase mobility (not including density)
265 
266  // density
267  real64 const resDens = m_resPhaseDens[er][esr][ei][0][ip];
268  real64 dDens[CP_Deriv::nDer]{};
269 
270  dDens[CP_Deriv::dP] = m_dResPhaseDens[er][esr][ei][0][ip][Deriv::dP];
271  if constexpr ( IS_THERMAL )
272  {
273  dDens[CP_Deriv::dT] = m_dResPhaseDens[er][esr][ei][0][ip][Deriv::dT];
274  }
275  applyChainRule( NC, m_dResCompFrac_dCompDens[er][esr][ei],
276  m_dResPhaseDens[er][esr][ei][0][ip],
277  &dDens[CP_Deriv::dC],
278  Deriv::dC );
279  // viscosity
280  real64 const resVisc = m_resPhaseVisc[er][esr][ei][0][ip];
281  real64 dVisc[CP_Deriv::nDer]{};
282  dVisc[CP_Deriv::dP] = m_dResPhaseVisc[er][esr][ei][0][ip][Deriv::dP];
283  if constexpr ( IS_THERMAL )
284  {
285  dVisc[CP_Deriv::dT] = m_dResPhaseVisc[er][esr][ei][0][ip][Deriv::dT];
286  }
287 
288  applyChainRule( NC, m_dResCompFrac_dCompDens[er][esr][ei],
289  m_dResPhaseVisc[er][esr][ei][0][ip],
290  &dVisc[CP_Deriv::dC],
291  Deriv::dC );
292 
293  // relative permeability
294  real64 const resRelPerm = m_resPhaseRelPerm[er][esr][ei][0][ip];
295  real64 dRelPerm[CP_Deriv::nDer]{};
296  for( integer jc = 0; jc < CP_Deriv::nDer; ++jc )
297  {
298  dRelPerm[jc]=0;
299  }
300  for( integer jp = 0; jp < NP; ++jp )
301  {
302  real64 const dResRelPerm_dS = m_dResPhaseRelPerm_dPhaseVolFrac[er][esr][ei][0][ip][jp];
303  dRelPerm[CP_Deriv::dP] += dResRelPerm_dS * m_dResPhaseVolFrac[er][esr][ei][jp][Deriv::dP];
304  if constexpr ( IS_THERMAL )
305  {
306  dRelPerm[CP_Deriv::dT] += dResRelPerm_dS * m_dResPhaseVolFrac[er][esr][ei][jp][Deriv::dT];
307  }
308  for( integer jc = 0; jc < NC; ++jc )
309  {
310  dRelPerm[CP_Deriv::dC+jc] += dResRelPerm_dS * m_dResPhaseVolFrac[er][esr][ei][jp][Deriv::dC+jc];
311  }
312  }
313 
314  // compute the reservoir phase mobility, including phase density
315  real64 const resPhaseMob = resDens * resRelPerm / resVisc;
316 
317  // Handles all dependencies
318  for( integer jc = 0; jc < CP_Deriv::nDer; ++jc )
319  {
320  dMob[jc] = dRelPerm[jc] * resDens / resVisc
321  + resPhaseMob * (dDens[jc] / resDens - dVisc[jc] / resVisc);
322  }
323  // compute the phase flux and derivatives using upstream cell mobility
324  flux = resPhaseMob * potDiff;
325  // Handles all dependencies
326  for( integer jc = 0; jc < CP_Deriv::nDer; ++jc )
327  {
328  dFlux[TAG::RES][jc] = dMob[jc] * potDiff + resPhaseMob * dPotDiff[TAG::RES][jc];
329  dFlux[TAG::WELL][jc] = resPhaseMob * dPotDiff[TAG::WELL][jc];
330  }
331 
332  // increment component fluxes
333  for( integer ic = 0; ic < NC; ++ic )
334  {
335  // Note this needs to be uncommented out
336  m_compPerfRate[iperf][ic] += flux * m_resPhaseCompFrac[er][esr][ei][0][ip][ic];
337  dCompFrac[CP_Deriv::dP] = m_dResPhaseCompFrac[er][esr][ei][0][ip][ic][Deriv::dP];
338  if constexpr (IS_THERMAL)
339  {
340  dCompFrac[CP_Deriv::dT] = m_dResPhaseCompFrac[er][esr][ei][0][ip][ic][Deriv::dT];
341  }
342 
343  applyChainRule( NC,
344  m_dResCompFrac_dCompDens[er][esr][ei],
345  m_dResPhaseCompFrac[er][esr][ei][0][ip][ic],
346  &dCompFrac[CP_Deriv::dC],
347  Deriv::dC );
348 
349  for( integer jc = 0; jc < CP_Deriv::nDer; ++jc )
350  {
351  m_dCompPerfRate[iperf][TAG::RES][ic][jc] += dFlux[TAG::RES][jc] * m_resPhaseCompFrac[er][esr][ei][0][ip][ic];
352  m_dCompPerfRate[iperf][TAG::RES][ic][jc] += flux * dCompFrac[jc];
353  m_dCompPerfRate[iperf][TAG::WELL][ic][jc] += dFlux[TAG::WELL][jc] * m_resPhaseCompFrac[er][esr][ei][0][ip][ic];
354  }
355  }
356  if constexpr ( IS_THERMAL )
357  {
358  fluxKernelOp( iwelem, er, esr, ei, ip, potDiff, flux, dFlux );
359  }
360 
361  } // end resevoir is upstream phase loop
362 
363  }
364  else // ** well is upstream **
365  {
366 
367  real64 resTotalMob = 0.0;
368 
369  // we re-compute here the total mass (when useMass == 1) or molar (when useMass == 0) density
370  real64 wellElemTotalDens = 0;
371  for( integer ic = 0; ic < NC; ++ic )
372  {
373  wellElemTotalDens += m_wellElemCompDens[iwelem][ic];
374  }
375 
376  // first, compute the reservoir total mobility (excluding phase density)
377  for( integer ip = 0; ip < NP; ++ip )
378  {
379 
380  // skip the rest of the calculation if the phase is absent
381  // or if crossflow is disabled for non-injectors aka producers
382  bool const phaseExists = (m_resPhaseVolFrac[er][esr][ei][ip] > 0);
383  if( !phaseExists || (!m_isInjector && !m_isCrossflowEnabled) )
384  {
385  continue;
386  }
387 
388  // viscosity
389  real64 const resVisc = m_resPhaseVisc[er][esr][ei][0][ip];
390  real64 dVisc[CP_Deriv::nDer]{};
391  dVisc[CP_Deriv::dP] = m_dResPhaseVisc[er][esr][ei][0][ip][Deriv::dP];
392  if constexpr ( IS_THERMAL )
393  {
394  dVisc[CP_Deriv::dT] = m_dResPhaseVisc[er][esr][ei][0][ip][Deriv::dT];
395  }
396 
397  applyChainRule( NC, m_dResCompFrac_dCompDens[er][esr][ei],
398  m_dResPhaseVisc[er][esr][ei][0][ip],
399  &dVisc[CP_Deriv::dC],
400  Deriv::dC );
401 
402 
403  // relative permeability
404  real64 const resRelPerm = m_resPhaseRelPerm[er][esr][ei][0][ip];
405  real64 dRelPerm[CP_Deriv::nDer]{};
406  for( integer jc = 0; jc < CP_Deriv::nDer; ++jc )
407  {
408  dRelPerm[jc]=0;
409  }
410  for( integer jp = 0; jp < NP; ++jp )
411  {
412  real64 const dResRelPerm_dS = m_dResPhaseRelPerm_dPhaseVolFrac[er][esr][ei][0][ip][jp];
413  dRelPerm[CP_Deriv::dP] += dResRelPerm_dS * m_dResPhaseVolFrac[er][esr][ei][jp][Deriv::dP];
414  if constexpr ( IS_THERMAL )
415  {
416  dRelPerm[CP_Deriv::dT] += dResRelPerm_dS * m_dResPhaseVolFrac[er][esr][ei][jp][Deriv::dT];
417  }
418  for( integer jc = 0; jc < NC; ++jc )
419  {
420  dRelPerm[CP_Deriv::dC+jc] += dResRelPerm_dS * m_dResPhaseVolFrac[er][esr][ei][jp][Deriv::dC+jc];
421  }
422  }
423  // increment total mobility
424  resTotalMob += resRelPerm / resVisc;
425  // Handles all dependencies
426  for( integer jc = 0; jc < CP_Deriv::nDer; ++jc )
427  {
428  dMob[jc] += (dRelPerm[jc] *resVisc - resRelPerm * dVisc[jc] )
429  / ( resVisc * resVisc);
430  }
431  } // end well is upstream phase loop
432 
433  // compute a potdiff multiplier = wellElemTotalDens * resTotalMob
434  // wellElemTotalDens is a mass density if useMass == 1 and a molar density otherwise
435  real64 const mult = wellElemTotalDens * resTotalMob;
436 
437  real64 dMult[2][CP_Deriv::nDer]{};
438  dMult[TAG::WELL][CP_Deriv::dP] = 0.0;
439  if constexpr ( IS_THERMAL )
440  {
441  dMult[TAG::WELL][CP_Deriv::dT] = 0.0;
442  }
443  for( integer ic = 0; ic < NC; ++ic )
444  {
445  dMult[TAG::WELL][CP_Deriv::dC+ic] = resTotalMob;
446  }
447  for( integer jc = 0; jc < CP_Deriv::nDer; ++jc )
448  {
449  dMult[TAG::RES][jc] = wellElemTotalDens * dMob[jc];
450  }
451 
452 
453  // compute the volumetric flux and derivatives using upstream cell mobility
454  flux = mult * potDiff;
455 
456  for( integer ic = 0; ic < CP_Deriv::nDer; ++ic )
457  {
458  dFlux[TAG::RES][ic] = dMult[TAG::RES][ic] * potDiff + mult * dPotDiff[TAG::RES][ic];
459  dFlux[TAG::WELL][ic] = dMult[TAG::WELL][ic] * potDiff + mult * dPotDiff[TAG::WELL][ic];
460  }
461  // compute component fluxes
462  for( integer ic = 0; ic < NC; ++ic )
463  {
464  m_compPerfRate[iperf][ic] += m_wellElemCompFrac[iwelem][ic] * flux;
465  for( integer jc = 0; jc < CP_Deriv::nDer; ++jc )
466  {
467  m_dCompPerfRate[iperf][TAG::RES][ic][jc] = m_wellElemCompFrac[iwelem][ic] * dFlux[TAG::RES][jc];
468  }
469  }
470  for( integer ic = 0; ic < NC; ++ic )
471  {
472  m_dCompPerfRate[iperf][TAG::WELL][ic][CP_Deriv::dP] = m_wellElemCompFrac[iwelem][ic] * dFlux[TAG::WELL][CP_Deriv::dP];
473  if constexpr ( IS_THERMAL )
474  {
475  m_dCompPerfRate[iperf][TAG::WELL][ic][CP_Deriv::dT] = m_wellElemCompFrac[iwelem][ic] * dFlux[TAG::WELL][CP_Deriv::dT];
476  }
477  for( integer jc = 0; jc < NC; ++jc )
478  {
479  m_dCompPerfRate[iperf][TAG::WELL][ic][CP_Deriv::dC+jc] += m_wellElemCompFrac[iwelem][ic] * dFlux[TAG::WELL][CP_Deriv::dC+jc];
480  m_dCompPerfRate[iperf][TAG::WELL][ic][CP_Deriv::dC+jc] += m_dWellElemCompFrac_dCompDens[iwelem][ic][jc] * flux;
481  }
482  }
483  if constexpr ( IS_THERMAL )
484  {
485  fluxKernelOp( iwelem, er, esr, ei, -1, potDiff, flux, dFlux );
486  }
487  } // end upstream
488  }
496  template< typename POLICY, typename KERNEL_TYPE >
497  static void
498  launch( localIndex const numElements,
499  KERNEL_TYPE const & kernelComponent )
500  {
501  GEOS_MARK_FUNCTION;
502  forAll< POLICY >( numElements, [=] GEOS_HOST_DEVICE ( localIndex const iperf )
503  {
504  kernelComponent.computeFlux( iperf );
505  } );
506  }
507 
508 
509 protected:
510  ElementViewConst< arrayView1d< real64 const > > const m_resPres;
511  ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const m_resPhaseVolFrac;
512  ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const m_dResPhaseVolFrac;
513  ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const m_dResCompFrac_dCompDens;
514  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const m_resPhaseDens;
515  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const m_dResPhaseDens;
516  ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const m_resPhaseVisc;
517  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const m_dResPhaseVisc;
518  ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const m_resPhaseCompFrac;
519  ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const m_dResPhaseCompFrac;
520  ElementViewConst< arrayView3d< real64 const, constitutive::relperm::USD_RELPERM > > const m_resPhaseRelPerm;
521  ElementViewConst< arrayView4d< real64 const, constitutive::relperm::USD_RELPERM_DS > > const m_dResPhaseRelPerm_dPhaseVolFrac;
522  arrayView1d< real64 const > const m_wellElemGravCoef;
523  arrayView1d< real64 const > const m_wellElemPres;
524  arrayView2d< real64 const, compflow::USD_COMP > const m_wellElemCompDens;
525  arrayView1d< real64 const > const m_wellElemTotalMassDens;
526  arrayView2d< real64 const, compflow::USD_FLUID_DC > const m_dWellElemTotalMassDens;
527  arrayView2d< real64 const, compflow::USD_COMP > const m_wellElemCompFrac;
528  arrayView3d< real64 const, compflow::USD_COMP_DC > const m_dWellElemCompFrac_dCompDens;
529  arrayView1d< real64 const > const m_perfGravCoef;
530  arrayView1d< localIndex const > const m_perfWellElemIndex;
531  arrayView1d< real64 const > const m_perfTrans;
532  arrayView1d< localIndex const > const m_resElementRegion;
533  arrayView1d< localIndex const > const m_resElementSubRegion;
534  arrayView1d< localIndex const > const m_resElementIndex;
535  arrayView2d< real64 > const m_compPerfRate;
536  arrayView4d< real64 > const m_dCompPerfRate;
537  arrayView3d< real64 > const m_dCompPerfRate_dPres;
538  arrayView4d< real64 > const m_dCompPerfRate_dComp;
539  arrayView1d< integer > const m_perfStatus;
540 
541  bool const m_isInjector;
542  bool const m_isCrossflowEnabled;
543 
544 };
545 
549 class PerforationFluxKernelFactory
550 {
551 public:
552 
565  template< typename POLICY >
566  static void
567  createAndLaunch( integer const numComp,
568  integer const numPhases,
569  string const flowSolverName,
570  PerforationData * const perforationData,
571  ElementSubRegionBase const & subRegion,
572  ElementRegionManager & elemManager,
573  bool const isInjector,
574  bool const isCrossflowEnabled )
575  {
576  geos::internal::kernelLaunchSelectorCompPhaseSwitch( numComp, numPhases, [&]( auto NC, auto NP )
577  {
578  integer constexpr NUM_COMP = NC();
579  integer constexpr NUM_PHASE = NP();
580  integer constexpr IS_THERMAL = 0;
581 
582  using kernelType = PerforationFluxKernel< NUM_COMP, NUM_PHASE, IS_THERMAL >;
583  typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, flowSolverName );
584  typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, flowSolverName );
585  typename kernelType::RelPermAccessors relPermAccessors( elemManager, flowSolverName );
586 
587  kernelType kernel( perforationData, subRegion,
588  compFlowAccessors, multiFluidAccessors, relPermAccessors,
589  isInjector, isCrossflowEnabled );
590  kernelType::template launch< POLICY >( perforationData->size(), kernel );
591  } );
592  }
593 };
594 
595 } // end namespace isothermalPerforationFluxKernels
596 
597 namespace thermalPerforationFluxKernels
598 {
599 
600 using namespace constitutive;
601 
602 /******************************** PerforationFluxKernel ********************************/
603 
604 template< integer NC, integer NP, integer IS_THERMAL >
605 class PerforationFluxKernel : public isothermalPerforationFluxKernels::PerforationFluxKernel< NC, NP, IS_THERMAL >
606 {
607 public:
608 
609  using Base = isothermalPerforationFluxKernels::PerforationFluxKernel< NC, NP, IS_THERMAL >;
610  using Base::m_resPhaseCompFrac;
611  using Base::m_dResCompFrac_dCompDens;
612  using Base::m_dWellElemCompFrac_dCompDens;
613 
615  static constexpr integer numComp = NC;
616 
618  static constexpr integer numPhase = NP;
619 
621  static constexpr integer isThermal = IS_THERMAL;
622 
623  using TAG = typename Base::TAG;
624  using CompFlowAccessors = typename Base::CompFlowAccessors;
625  using MultiFluidAccessors = typename Base::MultiFluidAccessors;
626  using RelPermAccessors = typename Base::RelPermAccessors;
627 
628 
629  using ThermalCompFlowAccessors =
630  StencilAccessors< fields::flow::temperature >;
631 
632  using ThermalMultiFluidAccessors =
633  StencilMaterialAccessors< MultiFluidBase,
634  fields::multifluid::phaseEnthalpy,
635  fields::multifluid::dPhaseEnthalpy >;
636 
644  template< typename VIEWTYPE >
645  using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
646 
647  PerforationFluxKernel ( PerforationData * const perforationData,
648  ElementSubRegionBase const & subRegion,
649  MultiFluidBase const & wellFluid,
650  CompFlowAccessors const & compFlowAccessors,
651  MultiFluidAccessors const & multiFluidAccessors,
652  RelPermAccessors const & relPermAccessors,
653  ThermalCompFlowAccessors const & thermalCompFlowAccessors,
654  ThermalMultiFluidAccessors const & thermalMultiFluidAccessors,
655  bool const isInjector,
656  bool const isCrossflowEnabled )
657  : Base( perforationData,
658  subRegion,
659  compFlowAccessors,
660  multiFluidAccessors,
661  relPermAccessors,
662  isInjector,
663  isCrossflowEnabled ),
664  m_wellElemPhaseFrac( wellFluid.phaseFraction() ),
665  m_dWellElemPhaseFrac( wellFluid.dPhaseFraction() ),
666  m_wellElemPhaseEnthalpy( wellFluid.phaseEnthalpy()),
667  m_dWellElemPhaseEnthalpy( wellFluid.dPhaseEnthalpy()),
668  m_energyPerfFlux( perforationData->getField< fields::well::energyPerforationFlux >() ),
669  m_dEnergyPerfFlux( perforationData->getField< fields::well::dEnergyPerforationFlux >() ),
670  m_temp( thermalCompFlowAccessors.get( fields::flow::temperature {} ) ),
671  m_resPhaseEnthalpy( thermalMultiFluidAccessors.get( fields::multifluid::phaseEnthalpy {} ) ),
672  m_dResPhaseEnthalpy( thermalMultiFluidAccessors.get( fields::multifluid::dPhaseEnthalpy {} ) )
673  {}
674 
675  template< typename FUNC = NoOpFunc >
676  GEOS_HOST_DEVICE
677  inline
678  void
679  computeFlux( localIndex const iperf ) const
680  {
681  using Deriv = constitutive::multifluid::DerivativeOffset;
682  using CP_Deriv =constitutive::multifluid::DerivativeOffsetC< NC, IS_THERMAL >;
683  // initialize outputs
684  m_energyPerfFlux[iperf]=0;
685  for( integer ke = 0; ke < 2; ++ke )
686  {
687  for( integer i = 0; i < CP_Deriv::nDer; ++i )
688  {
689  m_dEnergyPerfFlux[iperf][ke][i]=0;
690  }
691  }
692 
693  Base::computeFlux ( iperf, [&]( localIndex const iwelem, localIndex const er, localIndex const esr, localIndex const ei, localIndex const ip,
694  real64 const potDiff, real64 const flux, real64 const (&dFlux)[2][CP_Deriv::nDer] )
695  {
696  if( potDiff >= 0 ) // ** reservoir cell is upstream **
697  {
698 
699  real64 const res_enthalpy = m_resPhaseEnthalpy[er][esr][ei][0][ip];
700 
701  m_energyPerfFlux[iperf] += flux * res_enthalpy;
702 
703  // energy equation derivatives WRT res P & T
704  m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dP] += dFlux[TAG::RES][CP_Deriv::dP] * res_enthalpy +
705  flux * m_dResPhaseEnthalpy[er][esr][ei][0][ip][Deriv::dP];
706  m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dT] += dFlux[TAG::RES][CP_Deriv::dT] * res_enthalpy +
707  flux * m_dResPhaseEnthalpy[er][esr][ei][0][ip][Deriv::dT];
708  // energy equation derivatives WRT well P
709  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dP] += dFlux[TAG::WELL][CP_Deriv::dP] * res_enthalpy;
710  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dT] += dFlux[TAG::WELL][CP_Deriv::dT] * res_enthalpy;
711 
712 
713  // energy equation derivatives WRT reservoir dens
714  real64 dProp_dC[numComp]{};
715  applyChainRule( NC,
716  m_dResCompFrac_dCompDens[er][esr][ei],
717  m_dResPhaseEnthalpy[er][esr][ei][0][ip],
718  dProp_dC,
719  Deriv::dC );
720 
721  for( integer jc = 0; jc < NC; ++jc )
722  {
723  m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dC+jc] += flux * dProp_dC[jc];
724  }
725  }
726  else // ** reservoir cell is downstream
727  {
728  for( integer iphase = 0; iphase < NP; ++iphase )
729  {
730  bool const phaseExists = m_wellElemPhaseFrac[iwelem][0][iphase] > 0.0;
731  if( !phaseExists )
732  continue;
733  real64 const phaseFlux = m_wellElemPhaseFrac[iwelem][0][iphase] * flux;
734  real64 const wellelem_enthalpy = m_wellElemPhaseEnthalpy[iwelem][0][iphase];
735  m_energyPerfFlux[iperf] += phaseFlux * wellelem_enthalpy;
736 
737  // energy equation derivatives WRT res P & T
738  m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dP] += dFlux[TAG::RES][CP_Deriv::dP] * wellelem_enthalpy;
739  m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dT] += dFlux[TAG::RES][CP_Deriv::dT] * wellelem_enthalpy;
740 
741  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dP] += dFlux[TAG::WELL][CP_Deriv::dP] * wellelem_enthalpy
742  + phaseFlux * m_dWellElemPhaseEnthalpy[iwelem][0][iphase][Deriv::dP]
743  + phaseFlux * wellelem_enthalpy * m_dWellElemPhaseFrac[iwelem][0][iphase][Deriv::dP];
744  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dT] += dFlux[TAG::WELL][CP_Deriv::dT] * wellelem_enthalpy
745  + phaseFlux * m_dWellElemPhaseEnthalpy[iwelem][0][iphase][Deriv::dT]
746  + phaseFlux * wellelem_enthalpy * m_dWellElemPhaseFrac[iwelem][0][iphase][Deriv::dT];
747 
748  //energy e
749  real64 dPVF_dC[numComp]{};
750  applyChainRule( NC,
751  m_dWellElemCompFrac_dCompDens[iwelem],
752  m_dWellElemPhaseFrac[iwelem][0][iphase],
753  dPVF_dC,
754  Deriv::dC );
755  for( integer ic=0; ic<NC; ic++ )
756  {
757  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dC+ic] += wellelem_enthalpy * dFlux[TAG::WELL][CP_Deriv::dC+ic] * m_wellElemPhaseFrac[iwelem][0][iphase];
758  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dC+ic] += wellelem_enthalpy * phaseFlux * dPVF_dC[ic];
759  }
760  // energy equation enthalpy derivatives WRT well dens
761  real64 dProp_dC[numComp]{};
762  applyChainRule( NC,
763  m_dWellElemCompFrac_dCompDens[iwelem],
764  m_dWellElemPhaseEnthalpy[iwelem][0][iphase],
765  dProp_dC,
766  Deriv::dC );
767 
768  for( integer jc = 0; jc < NC; ++jc )
769  {
770  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dC+jc] += phaseFlux * dProp_dC[jc];
771  }
772  }
773 
774  }
775  } );
776  }
777 
778 
779 
788  template< typename POLICY, typename KERNEL_TYPE >
789  static void
790  launch( localIndex const numElements,
791  KERNEL_TYPE const & kernelComponent )
792  {
793  GEOS_MARK_FUNCTION;
794  forAll< POLICY >( numElements, [=] GEOS_HOST_DEVICE ( localIndex const iperf )
795  {
796  kernelComponent.computeFlux( iperf );
797  } );
798  }
799 
800 protected:
801 
804  arrayView3d< real64 const, multifluid::USD_PHASE > const m_wellElemPhaseFrac;
805  arrayView4d< real64 const, multifluid::USD_PHASE_DC > const m_dWellElemPhaseFrac;
806  arrayView3d< real64 const, multifluid::USD_PHASE > const m_wellElemPhaseEnthalpy;
807  arrayView4d< real64 const, multifluid::USD_PHASE_DC > const m_dWellElemPhaseEnthalpy;
808 
810  arrayView1d< real64 > const m_energyPerfFlux;
811  arrayView3d< real64 > const m_dEnergyPerfFlux;
812 
814  ElementViewConst< arrayView1d< real64 const > > const m_temp;
815 
817  ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const m_resPhaseEnthalpy;
818  ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const m_dResPhaseEnthalpy;
819 
820 
821 };
822 
826 class PerforationFluxKernelFactory
827 {
828 public:
829 
843  template< typename POLICY >
844  static void
845  createAndLaunch( integer const numComp,
846  integer const numPhases,
847  string const flowSolverName,
848  PerforationData * const perforationData,
849  ElementSubRegionBase const & subRegion,
850  MultiFluidBase const & fluid,
851  ElementRegionManager & elemManager,
852  bool const isInjector,
853  bool const isCrossflowEnabled )
854  {
855  geos::internal::kernelLaunchSelectorCompPhaseSwitch( numComp, numPhases, [&]( auto NC, auto NP )
856  {
857  integer constexpr NUM_COMP = NC();
858  integer constexpr NUM_PHASE = NP();
859  integer constexpr IS_THERMAL = 1;
860 
861  using kernelType = PerforationFluxKernel< NUM_COMP, NUM_PHASE, IS_THERMAL >;
862  typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, flowSolverName );
863  typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, flowSolverName );
864  typename kernelType::RelPermAccessors relPermAccessors( elemManager, flowSolverName );
865  typename kernelType::ThermalCompFlowAccessors thermalCompFlowAccessors( elemManager, flowSolverName );
866  typename kernelType::ThermalMultiFluidAccessors thermalMultiFluidAccessors( elemManager, flowSolverName );
867 
868  kernelType kernel( perforationData, subRegion, fluid, compFlowAccessors, multiFluidAccessors,
869  relPermAccessors, thermalCompFlowAccessors, thermalMultiFluidAccessors,
870  isInjector, isCrossflowEnabled );
871  kernelType::template launch< POLICY >( perforationData->size(), kernel );
872  } );
873  }
874 };
875 
876 } // end namespace thermalPerforationFluxKernels
877 
878 } // end namespace geos
879 
880 #endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_PERFORATIONFLUXLKERNELS_HPP
GEOS_HOST_DEVICE
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
GEOS_MARK_FUNCTION
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
Definition: TimingMacros.hpp:134
geos::ElementRegionManager
The ElementRegionManager class provides an interface to ObjectManagerBase in order to manage ElementR...
Definition: ElementRegionManager.hpp:45
geos::ElementRegionManager::ElementViewConst
typename ElementViewAccessor< VIEWTYPE >::NestedViewTypeConst ElementViewConst
The ElementViewAccessor at the ElementRegionManager level is the type resulting from ElementViewAcces...
Definition: ElementRegionManager.hpp:83
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::size
localIndex size() const
Get the "size" of the group, which determines the number of elements in resizable wrappers.
Definition: Group.hpp:1317
geos::isothermalPerforationFluxKernels::PerforationFluxKernelFactory::createAndLaunch
static void createAndLaunch(integer const numComp, integer const numPhases, string const flowSolverName, PerforationData *const perforationData, ElementSubRegionBase const &subRegion, ElementRegionManager &elemManager, bool const isInjector, bool const isCrossflowEnabled)
Create a new kernel and launch.
Definition: PerforationFluxKernels.hpp:567
geos::isothermalPerforationFluxKernels::PerforationFluxKernel::numComp
static constexpr integer numComp
Compile time value for the number of components.
Definition: PerforationFluxKernels.hpp:60
geos::isothermalPerforationFluxKernels::PerforationFluxKernel::numPhase
static constexpr integer numPhase
Compile time value for the number of phases.
Definition: PerforationFluxKernels.hpp:63
geos::isothermalPerforationFluxKernels::PerforationFluxKernel::isThermal
static constexpr integer isThermal
Compile time value for thermal option.
Definition: PerforationFluxKernels.hpp:66
geos::isothermalPerforationFluxKernels::PerforationFluxKernel::launch
static void launch(localIndex const numElements, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
Definition: PerforationFluxKernels.hpp:498
geos::isothermalPerforationFluxKernels::PerforationFluxKernel::ElementViewConst
ElementRegionManager::ElementViewConst< VIEWTYPE > ElementViewConst
The type for element-based non-constitutive data parameters. Consists entirely of ArrayView's.
Definition: PerforationFluxKernels.hpp:99
geos::thermalPerforationFluxKernels::PerforationFluxKernelFactory::createAndLaunch
static void createAndLaunch(integer const numComp, integer const numPhases, string const flowSolverName, PerforationData *const perforationData, ElementSubRegionBase const &subRegion, MultiFluidBase const &fluid, ElementRegionManager &elemManager, bool const isInjector, bool const isCrossflowEnabled)
Create a new kernel and launch.
Definition: PerforationFluxKernels.hpp:845
geos::thermalPerforationFluxKernels::PerforationFluxKernel::m_wellElemPhaseFrac
arrayView3d< real64 const, multifluid::USD_PHASE > const m_wellElemPhaseFrac
Definition: PerforationFluxKernels.hpp:804
geos::thermalPerforationFluxKernels::PerforationFluxKernel::m_resPhaseEnthalpy
ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const m_resPhaseEnthalpy
Views on phase enthalpies.
Definition: PerforationFluxKernels.hpp:817
geos::thermalPerforationFluxKernels::PerforationFluxKernel::m_energyPerfFlux
arrayView1d< real64 > const m_energyPerfFlux
Views on energy flux.
Definition: PerforationFluxKernels.hpp:810
geos::thermalPerforationFluxKernels::PerforationFluxKernel::ElementViewConst
ElementRegionManager::ElementViewConst< VIEWTYPE > ElementViewConst
The type for element-based non-constitutive data parameters. Consists entirely of ArrayView's.
Definition: PerforationFluxKernels.hpp:645
geos::thermalPerforationFluxKernels::PerforationFluxKernel::m_temp
ElementViewConst< arrayView1d< real64 const > > const m_temp
Views on temperature.
Definition: PerforationFluxKernels.hpp:814
geos::thermalPerforationFluxKernels::PerforationFluxKernel::launch
static void launch(localIndex const numElements, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
Definition: PerforationFluxKernels.hpp:790
geos::arrayView1d
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:179
geos::real64
double real64
64-bit floating point type.
Definition: DataTypes.hpp:98
geos::localIndex
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:84
geos::arrayView4d
ArrayView< T, 4, USD > arrayView4d
Alias for 4D array view.
Definition: DataTypes.hpp:227
geos::arrayView2d
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:195
geos::integer
int integer
Signed integer type.
Definition: DataTypes.hpp:81
geos::arrayView3d
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
Definition: DataTypes.hpp:211
Generated by doxygen 1.9.1