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