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