GEOS
SinglePhasePerforationFluxKernels.hpp
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_SINGLEPHASEPERFORATIONFLUXLKERNELS_HPP
21 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_SINGLEPHASEPERFORATIONFLUXLKERNELS_HPP
22 
23 #include "codingUtilities/Utilities.hpp"
24 #include "common/DataTypes.hpp"
25 #include "common/GEOS_RAJA_Interface.hpp"
28 #include "constitutive/fluid/singlefluid/SingleFluidLayouts.hpp"
29 
30 #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp"
31 #include "constitutive/fluid/singlefluid/SingleFluidFields.hpp"
32 
36 
39 
43 
44 
45 namespace geos
46 {
47 
48 struct NoOpStuct
49 {
50  NoOpStuct(){}
51 };
52 
53 namespace isothermalSinglePhasePerforationFluxKernels
54 {
55 
56 /******************************** PerforationFluxKernel ********************************/
57 
58 template< integer IS_THERMAL >
60 {
61 public:
62 
64  static constexpr integer isThermal = IS_THERMAL;
65 
69 
70  using SingleFluidAccessors =
71  StencilMaterialAccessors< constitutive::SingleFluidBase,
72  fields::singlefluid::density,
73  fields::singlefluid::dDensity,
74  fields::singlefluid::viscosity,
75  fields::singlefluid::dViscosity >;
76 
77 
85  template< typename VIEWTYPE >
87 
88  PerforationFluxKernel ( PerforationData * const perforationData,
89  ElementSubRegionBase const & subRegion,
90  constitutive::SingleFluidBase const & fluid,
91  SinglePhaseFlowAccessors const & singlePhaseFlowAccessors,
92  SingleFluidAccessors const & singleFluidAccessors
93  ):
94  m_resPres( singlePhaseFlowAccessors.get( fields::flow::pressure {} )),
95  m_resDens( singleFluidAccessors.get( fields::singlefluid::density {} )),
96  m_dResDens( singleFluidAccessors.get( fields::singlefluid::dDensity {} )),
97  m_resVisc( singleFluidAccessors.get( fields::singlefluid::viscosity {} )),
98  m_dResVisc( singleFluidAccessors.get( fields::singlefluid::dViscosity {} )),
99  m_wellElemGravCoef( subRegion.getField< fields::well::gravityCoefficient >()),
100  m_wellElemPres( subRegion.getField< fields::well::pressure >()),
101  m_wellElemDens( fluid.density()),
102  m_dWellElemDens( fluid.dDensity()),
103  m_wellElemVisc( fluid.viscosity()),
104  m_dWellElemVisc( fluid.dViscosity()),
105  m_perfGravCoef( perforationData->getField< fields::well::gravityCoefficient >()),
106  m_perfWellElemIndex( perforationData->getField< fields::perforation::wellElementIndex >()),
107  m_perfTrans( perforationData->getField< fields::perforation::wellTransmissibility >()),
108  m_resElementRegion( perforationData->getField< fields::perforation::reservoirElementRegion >()),
109  m_resElementSubRegion( perforationData->getField< fields::perforation::reservoirElementSubRegion >()),
110  m_resElementIndex( perforationData->getField< fields::perforation::reservoirElementIndex >()),
111  m_perfRate( perforationData->getField< fields::well::perforationRate >()),
112  m_dPerfRate( perforationData->getField< fields::well::dPerforationRate >()),
113  m_perfStatus( perforationData->getField< fields::perforation::perforationStatus >())
114  {}
115 
116  template< typename FUNC = NoOpFunc >
118  inline
119  void
120  computeFlux( localIndex const iperf, FUNC && fluxKernelOp= NoOpFunc {} ) const
121  {
122 
123  using Deriv = constitutive::singlefluid::DerivativeOffset;
124  using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< IS_THERMAL >;
125 
126  m_perfRate[iperf] = 0.0;
127  for( integer i = 0; i < 2; ++i )
128  {
129  for( integer j=0; j<DerivOffset::nDer; j++ )
130  {
131  m_dPerfRate[iperf][i][j] = 0.0;
132  }
133  }
134 
135  if( !m_perfStatus[iperf] )
136  {
137  return;
138  }
139  // get the reservoir (sub)region and element indices
140  localIndex const er = m_resElementRegion[iperf];
141  localIndex const esr = m_resElementSubRegion[iperf];
142  localIndex const ei = m_resElementIndex[iperf];
143 
144  // get the local index of the well element
145  localIndex const iwelem = m_perfWellElemIndex[iperf];
146 
147  // local working variables and arrays
148  real64 pressure[2]{};
149  real64 dPressure[2][DerivOffset::nDer]{};
150 
151  real64 multiplier[2]{};
152 
153  // 1) Reservoir side
154  // get reservoir variables
155  pressure[TAG::RES] = m_resPres[er][esr][ei];
156  dPressure[TAG::RES][DerivOffset::dP] = 1.0;
157 
158  // TODO: add a buoyancy term for the reservoir side here
159 
160  // multiplier for reservoir side in the flux
161  multiplier[TAG::RES] = 1;
162 
163  // 2) Well side
164 
165  // get well variables
166  pressure[TAG::WELL] = m_wellElemPres[iwelem];
167  dPressure[TAG::WELL][DerivOffset::dP] = 1.0;
168 
169  real64 const gravD = ( m_perfGravCoef[iperf] - m_wellElemGravCoef[iwelem] );
170  pressure[TAG::WELL] += m_wellElemDens[iwelem][0] * gravD;
171  dPressure[TAG::WELL][DerivOffset::dP] += m_dWellElemDens[iwelem][0][Deriv::dP] * gravD;
172  if constexpr (IS_THERMAL)
173  {
174  dPressure[TAG::WELL][DerivOffset::dT] = m_dWellElemDens[iwelem][0][DerivOffset::dT] * gravD;
175  }
176 
177  // multiplier for well side in the flux
178  multiplier[TAG::WELL] = -1;
179 
180  // compute potential difference
181  real64 potDif = 0.0;
182  for( localIndex i = 0; i < 2; ++i )
183  {
184  potDif += multiplier[i] * m_perfTrans[iperf] * pressure[i];
185  for( integer j=0; j<DerivOffset::nDer; j++ )
186  {
187  m_dPerfRate[iperf][i][j] = multiplier[i] * m_perfTrans[iperf] * dPressure[i][j];
188  }
189  }
190 
191  // choose upstream cell based on potential difference
192  localIndex const k_up = (potDif >= 0) ? TAG::RES : TAG::WELL;
193 
194  // compute upstream density, viscosity, and mobility
195  real64 densityUp = 0.0;
196  real64 viscosityUp = 0.0;
197  real64 dDensityUp[DerivOffset::nDer]{};
198  real64 dViscosityUp[DerivOffset::nDer]{};
199  for( integer j=0; j<DerivOffset::nDer; j++ )
200  {
201  dDensityUp[j]=0.0;
202  dViscosityUp[j]=0.0;
203  }
204 
205  // upwinding the variables
206  if( k_up == TAG::RES ) // use reservoir vars
207  {
208  densityUp = m_resDens[er][esr][ei][0];
209  viscosityUp = m_resVisc[er][esr][ei][0];
210  for( integer j=0; j<DerivOffset::nDer; j++ )
211  {
212  dDensityUp[j] = m_dResDens[er][esr][ei][0][j];
213  dViscosityUp[j] = m_dResVisc[er][esr][ei][0][j];
214  }
215  }
216  else // use well vars
217  {
218  densityUp = m_wellElemDens[iwelem][0];
219  viscosityUp = m_wellElemVisc[iwelem][0];
220  for( integer j=0; j<DerivOffset::nDer; j++ )
221  {
222  dDensityUp[j] = m_dWellElemDens[iwelem][0][j];
223  dViscosityUp[j] = m_dWellElemVisc[iwelem][0][j];
224  }
225  }
226 
227  // compute mobility
228  real64 const mobilityUp = densityUp / viscosityUp;
229  real64 dMobilityUp[DerivOffset::nDer]{};
230  for( integer j=0; j<DerivOffset::nDer; j++ )
231  {
232  dMobilityUp[j] = dDensityUp[j] / viscosityUp
233  - mobilityUp / viscosityUp * dViscosityUp[j];
234  }
235 
236  m_perfRate[iperf] = mobilityUp * potDif;
237  for( localIndex ke = 0; ke < 2; ++ke )
238  {
239  for( integer j=0; j<DerivOffset::nDer; j++ )
240  {
241  m_dPerfRate[iperf][ke][j] *= mobilityUp;
242  }
243  }
244  for( integer j=0; j<DerivOffset::nDer; j++ )
245  {
246  m_dPerfRate[iperf][k_up][j] += dMobilityUp[j]* potDif;
247  }
248  if constexpr ( IS_THERMAL )
249  {
250  fluxKernelOp( iwelem, er, esr, ei, potDif, m_perfRate[iperf], m_dPerfRate[iperf] );
251  }
252  }
260  template< typename POLICY, typename KERNEL_TYPE >
261  static void
262  launch( localIndex const numElements,
263  KERNEL_TYPE const & kernelComponent )
264  {
266  forAll< POLICY >( numElements, [=] GEOS_HOST_DEVICE ( localIndex const iperf )
267  {
268  kernelComponent.computeFlux( iperf );
269  } );
270  }
271 
272 
273 protected:
274  ElementViewConst< arrayView1d< real64 const > > const m_resPres;
275  ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > > const m_resDens;
276  ElementViewConst< arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > > const m_dResDens;
277  ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > > const m_resVisc;
278  ElementViewConst< arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > > const m_dResVisc;
279  arrayView1d< real64 const > const m_wellElemGravCoef;
280  arrayView1d< real64 const > const m_wellElemPres;
285 
286  arrayView1d< real64 const > const m_perfGravCoef;
287  arrayView1d< localIndex const > const m_perfWellElemIndex;
288  arrayView1d< real64 const > const m_perfTrans;
289  arrayView1d< localIndex const > const m_resElementRegion;
290  arrayView1d< localIndex const > const m_resElementSubRegion;
291  arrayView1d< localIndex const > const m_resElementIndex;
292  arrayView1d< real64 > const m_perfRate;
293  arrayView3d< real64 > const m_dPerfRate;
294  arrayView1d< integer > const m_perfStatus;
295 
296 };
297 
302 {
303 public:
304 
315  template< typename POLICY >
316  static void
317  createAndLaunch( string const flowSolverName,
318  PerforationData * const perforationData,
319  ElementSubRegionBase const & subRegion,
320  constitutive::SingleFluidBase const & fluid,
321  ElementRegionManager & elemManager )
322  {
323  integer constexpr IS_THERMAL = 0;
324  using kernelType = PerforationFluxKernel< IS_THERMAL >;
325  typename kernelType::SinglePhaseFlowAccessors singlePhaseFlowAccessors( elemManager, flowSolverName );
326  typename kernelType::SingleFluidAccessors singleFluidAccessors( elemManager, flowSolverName );
327  kernelType kernel( perforationData, subRegion, fluid, singlePhaseFlowAccessors, singleFluidAccessors );
328  kernelType::template launch< POLICY >( perforationData->size(), kernel );
329  }
330 };
331 
332 } // end namespace isothermalPerforationFluxKernels
333 
334 namespace thermalSinglePhasePerforationFluxKernels
335 {
336 
337 using namespace constitutive;
338 
339 /******************************** PerforationFluxKernel ********************************/
340 
341 template< integer IS_THERMAL >
343 {
344 public:
345 
349 
351  static constexpr integer isThermal = IS_THERMAL;
352 
353  using TAG = typename Base::TAG;
354 
357 
359  StencilMaterialAccessors< SingleFluidBase,
360  fields::singlefluid::enthalpy,
361  fields::singlefluid::dEnthalpy >;
362 
363 
371  template< typename VIEWTYPE >
373 
374  PerforationFluxKernel ( PerforationData * const perforationData,
375  ElementSubRegionBase const & subRegion,
376  SingleFluidBase const & fluid,
377  SinglePhaseFlowAccessors const & singlePhaseFlowAccessors,
378  SingleFluidAccessors const & singleFluidAccessors,
379  ThermalSinglePhaseFlowAccessors const & thermalSinglePhaseFlowAccessors,
380  ThermalSingleFluidAccessors const & thermalSingleFluidAccessors )
381  : Base( perforationData,
382  subRegion,
383  fluid,
384  singlePhaseFlowAccessors,
385  singleFluidAccessors ),
386  m_wellElemEnthalpy( fluid.enthalpy()),
387  m_dWellElemEnthalpy( fluid.dEnthalpy()),
388  m_energyPerfFlux( perforationData->getField< fields::well::energyPerforationFlux >()),
389  m_dEnergyPerfFlux( perforationData->getField< fields::well::dEnergyPerforationFlux >()),
390  m_temp( thermalSinglePhaseFlowAccessors.get( fields::flow::temperature {} ) ),
391  m_resEnthalpy( thermalSingleFluidAccessors.get( fields::singlefluid::enthalpy {} ) ),
392  m_dResEnthalpy( thermalSingleFluidAccessors.get( fields::singlefluid::dEnthalpy {} ) )
393 
394  {}
395 
396  template< typename FUNC = NoOpFunc >
398  inline
399  void
400  computeFlux( localIndex const iperf ) const
401  {
402  using Deriv = constitutive::singlefluid::DerivativeOffset;
403  using CP_Deriv =constitutive::singlefluid::DerivativeOffsetC< IS_THERMAL >;
404  // initialize outputs
405  m_energyPerfFlux[iperf]=0;
406  for( integer ke = 0; ke < 2; ++ke )
407  {
408  for( integer i = 0; i < CP_Deriv::nDer; ++i )
409  {
410  m_dEnergyPerfFlux[iperf][ke][i]=0;
411  }
412  }
413 
414  Base::computeFlux ( iperf, [&]( localIndex const iwelem, localIndex const er, localIndex const esr, localIndex const ei,
415  real64 const potDiff, real64 const flux, arraySlice2d< real64 > const dFlux )
416  {
417  if( potDiff >= 0 ) // ** reservoir cell is upstream **
418  {
419 
420  real64 const res_enthalpy = m_resEnthalpy[er][esr][ei][0];
421 
422  m_energyPerfFlux[iperf] += flux * res_enthalpy;
423 
424  // energy equation derivatives WRT res P & T
425  m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dP] += dFlux[TAG::RES][CP_Deriv::dP] * res_enthalpy +
426  flux * m_dResEnthalpy[er][esr][ei][0][Deriv::dP];
427  m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dT] += dFlux[TAG::RES][CP_Deriv::dT] * res_enthalpy +
428  flux * m_dResEnthalpy[er][esr][ei][0][Deriv::dT];
429  // energy equation derivatives WRT well P
430  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dP] += dFlux[TAG::WELL][CP_Deriv::dP] * res_enthalpy;
431  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dT] += dFlux[TAG::WELL][CP_Deriv::dT] * res_enthalpy;
432 
433  }
434  else // ** reservoir cell is downstream
435  {
436  real64 const wellelem_enthalpy = m_wellElemEnthalpy[iwelem][0];
437  m_energyPerfFlux[iperf] += flux * wellelem_enthalpy;
438 
439  // energy equation derivatives WRT res P & T
440  m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dP] += dFlux[TAG::RES][CP_Deriv::dP] * wellelem_enthalpy;
441  m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dT] += dFlux[TAG::RES][CP_Deriv::dT] * wellelem_enthalpy;
442 
443  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dP] += dFlux[TAG::WELL][CP_Deriv::dP] * wellelem_enthalpy
444  + flux * m_dWellElemEnthalpy[iwelem][0][Deriv::dP];
445  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dT] += dFlux[TAG::WELL][CP_Deriv::dT] * wellelem_enthalpy
446  + flux * m_dWellElemEnthalpy[iwelem][0][Deriv::dT];
447  }
448  } );
449 
450  }
451 
452 
453 
462  template< typename POLICY, typename KERNEL_TYPE >
463  static void
464  launch( localIndex const numElements,
465  KERNEL_TYPE const & kernelComponent )
466  {
468  forAll< POLICY >( numElements, [=] GEOS_HOST_DEVICE ( localIndex const iperf )
469  {
470  kernelComponent.computeFlux( iperf );
471 
472  } );
473  }
474 
475 protected:
476 
480 
483  arrayView3d< real64 > const m_dEnergyPerfFlux;
484 
487 
491 
492 };
493 
498 {
499 public:
500 
511  template< typename POLICY >
512  static void
513  createAndLaunch( string const flowSolverName,
514  PerforationData * const perforationData,
515  ElementSubRegionBase const & subRegion,
516  SingleFluidBase const & fluid,
517  ElementRegionManager & elemManager )
518  {
519  integer constexpr IS_THERMAL = 1;
520  using kernelType = PerforationFluxKernel< IS_THERMAL >;
521  typename kernelType::SinglePhaseFlowAccessors singlePhaseFlowAccessors( elemManager, flowSolverName );
522  typename kernelType::SingleFluidAccessors singleFluidAccessors( elemManager, flowSolverName );
523  typename kernelType::ThermalSinglePhaseFlowAccessors thermalSinglePhaseFlowAccessors( elemManager, flowSolverName );
524  typename kernelType::ThermalSingleFluidAccessors thermalSingleFluidAccessors( elemManager, flowSolverName );
525  kernelType kernel( perforationData, subRegion, fluid, singlePhaseFlowAccessors, singleFluidAccessors, thermalSinglePhaseFlowAccessors, thermalSingleFluidAccessors );
526  kernelType::template launch< POLICY >( perforationData->size(), kernel );
527  }
528 };
529 
530 } // end namespace thermalPerforationFluxKernels
531 
532 } // end namespace geos
533 
534 #endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_SINGLEPHASEPERFORATIONFLUXLKERNELS_HPP
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
The ElementRegionManager class provides an interface to ObjectManagerBase in order to manage ElementR...
typename ElementViewAccessor< VIEWTYPE >::NestedViewTypeConst ElementViewConst
The ElementViewAccessor at the ElementRegionManager level is the type resulting from ElementViewAcces...
A struct to automatically construct and store element view accessors.
A struct to automatically construct and store element view accessors.
localIndex size() const
Get the "size" of the group, which determines the number of elements in resizable wrappers.
Definition: Group.hpp:1317
static void createAndLaunch(string const flowSolverName, PerforationData *const perforationData, ElementSubRegionBase const &subRegion, constitutive::SingleFluidBase const &fluid, ElementRegionManager &elemManager)
Create a new kernel and launch.
ElementRegionManager::ElementViewConst< VIEWTYPE > ElementViewConst
The type for element-based non-constitutive data parameters. Consists entirely of ArrayView's.
static constexpr integer isThermal
Compile time value for thermal option.
static void launch(localIndex const numElements, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
static void createAndLaunch(string const flowSolverName, PerforationData *const perforationData, ElementSubRegionBase const &subRegion, SingleFluidBase const &fluid, ElementRegionManager &elemManager)
Create a new kernel and launch.
ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > > const m_resEnthalpy
Views on phase enthalpies.
static void launch(localIndex const numElements, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
ElementRegionManager::ElementViewConst< VIEWTYPE > ElementViewConst
The type for element-based non-constitutive data parameters. Consists entirely of ArrayView's.
arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_wellElemEnthalpy
Views on well element properties.
ElementViewConst< arrayView1d< real64 const > > const m_temp
Views on temperature.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:179
double real64
64-bit floating point type.
Definition: DataTypes.hpp:98
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:84
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:195
int integer
Signed integer type.
Definition: DataTypes.hpp:81
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
Definition: DataTypes.hpp:211