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