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  {}
114 
115  template< typename FUNC = NoOpFunc >
117  inline
118  void
119  computeFlux( localIndex const iperf, FUNC && fluxKernelOp= NoOpFunc {} ) const
120  {
121 
122  // get the reservoir (sub)region and element indices
123  localIndex const er = m_resElementRegion[iperf];
124  localIndex const esr = m_resElementSubRegion[iperf];
125  localIndex const ei = m_resElementIndex[iperf];
126 
127  // get the local index of the well element
128  localIndex const iwelem = m_perfWellElemIndex[iperf];
129 
130  using Deriv = constitutive::singlefluid::DerivativeOffset;
131  using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< IS_THERMAL >;
132 
133  // local working variables and arrays
134  real64 pressure[2]{};
135  real64 dPressure[2][DerivOffset::nDer]{};
136 
137  real64 multiplier[2]{};
138 
139  m_perfRate[iperf] = 0.0;
140  for( integer i = 0; i < 2; ++i )
141  {
142  for( integer j=0; j<DerivOffset::nDer; j++ )
143  {
144  m_dPerfRate[iperf][i][j] = 0.0;
145  }
146  }
147 
148  // 1) Reservoir side
149  // get reservoir variables
150  pressure[TAG::RES] = m_resPres[er][esr][ei];
151  dPressure[TAG::RES][DerivOffset::dP] = 1.0;
152 
153  // TODO: add a buoyancy term for the reservoir side here
154 
155  // multiplier for reservoir side in the flux
156  multiplier[TAG::RES] = 1;
157 
158  // 2) Well side
159 
160  // get well variables
161  pressure[TAG::WELL] = m_wellElemPres[iwelem];
162  dPressure[TAG::WELL][DerivOffset::dP] = 1.0;
163 
164  real64 const gravD = ( m_perfGravCoef[iperf] - m_wellElemGravCoef[iwelem] );
165  pressure[TAG::WELL] += m_wellElemDens[iwelem][0] * gravD;
166  dPressure[TAG::WELL][DerivOffset::dP] += m_dWellElemDens[iwelem][0][Deriv::dP] * gravD;
167  if constexpr (IS_THERMAL)
168  {
169  dPressure[TAG::WELL][DerivOffset::dT] = m_dWellElemDens[iwelem][0][DerivOffset::dT] * gravD;
170  }
171 
172  // multiplier for well side in the flux
173  multiplier[TAG::WELL] = -1;
174 
175  // compute potential difference
176  real64 potDif = 0.0;
177  for( localIndex i = 0; i < 2; ++i )
178  {
179  potDif += multiplier[i] * m_perfTrans[iperf] * pressure[i];
180  for( integer j=0; j<DerivOffset::nDer; j++ )
181  {
182  m_dPerfRate[iperf][i][j] = multiplier[i] * m_perfTrans[iperf] * dPressure[i][j];
183  }
184  }
185 
186  // choose upstream cell based on potential difference
187  localIndex const k_up = (potDif >= 0) ? TAG::RES : TAG::WELL;
188 
189  // compute upstream density, viscosity, and mobility
190  real64 densityUp = 0.0;
191  real64 viscosityUp = 0.0;
192  real64 dDensityUp[DerivOffset::nDer]{};
193  real64 dViscosityUp[DerivOffset::nDer]{};
194  for( integer j=0; j<DerivOffset::nDer; j++ )
195  {
196  dDensityUp[j]=0.0;
197  dViscosityUp[j]=0.0;
198  }
199 
200  // upwinding the variables
201  if( k_up == TAG::RES ) // use reservoir vars
202  {
203  densityUp = m_resDens[er][esr][ei][0];
204  viscosityUp = m_resVisc[er][esr][ei][0];
205  for( integer j=0; j<DerivOffset::nDer; j++ )
206  {
207  dDensityUp[j] = m_dResDens[er][esr][ei][0][j];
208  dViscosityUp[j] = m_dResVisc[er][esr][ei][0][j];
209  }
210  }
211  else // use well vars
212  {
213  densityUp = m_wellElemDens[iwelem][0];
214  viscosityUp = m_wellElemVisc[iwelem][0];
215  for( integer j=0; j<DerivOffset::nDer; j++ )
216  {
217  dDensityUp[j] = m_dWellElemDens[iwelem][0][j];
218  dViscosityUp[j] = m_dWellElemVisc[iwelem][0][j];
219  }
220  }
221 
222  // compute mobility
223  real64 const mobilityUp = densityUp / viscosityUp;
224  real64 dMobilityUp[DerivOffset::nDer]{};
225  for( integer j=0; j<DerivOffset::nDer; j++ )
226  {
227  dMobilityUp[j] = dDensityUp[j] / viscosityUp
228  - mobilityUp / viscosityUp * dViscosityUp[j];
229  }
230 
231  m_perfRate[iperf] = mobilityUp * potDif;
232  for( localIndex ke = 0; ke < 2; ++ke )
233  {
234  for( integer j=0; j<DerivOffset::nDer; j++ )
235  {
236  m_dPerfRate[iperf][ke][j] *= mobilityUp;
237  }
238  }
239  for( integer j=0; j<DerivOffset::nDer; j++ )
240  {
241  m_dPerfRate[iperf][k_up][j] += dMobilityUp[j]* potDif;
242  }
243  if constexpr ( IS_THERMAL )
244  {
245  fluxKernelOp( iwelem, er, esr, ei, potDif, m_perfRate[iperf], m_dPerfRate[iperf] );
246  }
247  }
255  template< typename POLICY, typename KERNEL_TYPE >
256  static void
257  launch( localIndex const numElements,
258  KERNEL_TYPE const & kernelComponent )
259  {
261  forAll< POLICY >( numElements, [=] GEOS_HOST_DEVICE ( localIndex const iperf )
262  {
263  kernelComponent.computeFlux( iperf );
264  } );
265  }
266 
267 
268 protected:
269  ElementViewConst< arrayView1d< real64 const > > const m_resPres;
270  ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > > const m_resDens;
271  ElementViewConst< arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > > const m_dResDens;
272  ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > > const m_resVisc;
273  ElementViewConst< arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > > const m_dResVisc;
274  arrayView1d< real64 const > const m_wellElemGravCoef;
275  arrayView1d< real64 const > const m_wellElemPres;
280 
281  arrayView1d< real64 const > const m_perfGravCoef;
282  arrayView1d< localIndex const > const m_perfWellElemIndex;
283  arrayView1d< real64 const > const m_perfTrans;
284  arrayView1d< localIndex const > const m_resElementRegion;
285  arrayView1d< localIndex const > const m_resElementSubRegion;
286  arrayView1d< localIndex const > const m_resElementIndex;
287  arrayView1d< real64 > const m_perfRate;
288  arrayView3d< real64 > const m_dPerfRate;
289 
290 };
291 
296 {
297 public:
298 
309  template< typename POLICY >
310  static void
311  createAndLaunch( string const flowSolverName,
312  PerforationData * const perforationData,
313  ElementSubRegionBase const & subRegion,
314  constitutive::SingleFluidBase const & fluid,
315  ElementRegionManager & elemManager )
316  {
317  integer constexpr IS_THERMAL = 0;
318  using kernelType = PerforationFluxKernel< IS_THERMAL >;
319  typename kernelType::SinglePhaseFlowAccessors singlePhaseFlowAccessors( elemManager, flowSolverName );
320  typename kernelType::SingleFluidAccessors singleFluidAccessors( elemManager, flowSolverName );
321  kernelType kernel( perforationData, subRegion, fluid, singlePhaseFlowAccessors, singleFluidAccessors );
322  kernelType::template launch< POLICY >( perforationData->size(), kernel );
323  }
324 };
325 
326 } // end namespace isothermalPerforationFluxKernels
327 
328 namespace thermalSinglePhasePerforationFluxKernels
329 {
330 
331 using namespace constitutive;
332 
333 /******************************** PerforationFluxKernel ********************************/
334 
335 template< integer IS_THERMAL >
337 {
338 public:
339 
343 
345  static constexpr integer isThermal = IS_THERMAL;
346 
347  using TAG = typename Base::TAG;
348 
351 
353  StencilMaterialAccessors< SingleFluidBase,
354  fields::singlefluid::enthalpy,
355  fields::singlefluid::dEnthalpy >;
356 
357 
365  template< typename VIEWTYPE >
367 
368  PerforationFluxKernel ( PerforationData * const perforationData,
369  ElementSubRegionBase const & subRegion,
370  SingleFluidBase const & fluid,
371  SinglePhaseFlowAccessors const & singlePhaseFlowAccessors,
372  SingleFluidAccessors const & singleFluidAccessors,
373  ThermalSinglePhaseFlowAccessors const & thermalSinglePhaseFlowAccessors,
374  ThermalSingleFluidAccessors const & thermalSingleFluidAccessors )
375  : Base( perforationData,
376  subRegion,
377  fluid,
378  singlePhaseFlowAccessors,
379  singleFluidAccessors ),
380  m_wellElemEnthalpy( fluid.enthalpy()),
381  m_dWellElemEnthalpy( fluid.dEnthalpy()),
382  m_energyPerfFlux( perforationData->getField< fields::well::energyPerforationFlux >()),
383  m_dEnergyPerfFlux( perforationData->getField< fields::well::dEnergyPerforationFlux >()),
384  m_temp( thermalSinglePhaseFlowAccessors.get( fields::flow::temperature {} ) ),
385  m_resEnthalpy( thermalSingleFluidAccessors.get( fields::singlefluid::enthalpy {} ) ),
386  m_dResEnthalpy( thermalSingleFluidAccessors.get( fields::singlefluid::dEnthalpy {} ) )
387 
388  {}
389 
390  template< typename FUNC = NoOpFunc >
392  inline
393  void
394  computeFlux( localIndex const iperf ) const
395  {
396  using Deriv = constitutive::singlefluid::DerivativeOffset;
397  using CP_Deriv =constitutive::singlefluid::DerivativeOffsetC< IS_THERMAL >;
398  // initialize outputs
399  m_energyPerfFlux[iperf]=0;
400  for( integer ke = 0; ke < 2; ++ke )
401  {
402  for( integer i = 0; i < CP_Deriv::nDer; ++i )
403  {
404  m_dEnergyPerfFlux[iperf][ke][i]=0;
405  }
406  }
407 
408  Base::computeFlux ( iperf, [&]( localIndex const iwelem, localIndex const er, localIndex const esr, localIndex const ei,
409  real64 const potDiff, real64 const flux, arraySlice2d< real64 > const dFlux )
410  {
411  if( potDiff >= 0 ) // ** reservoir cell is upstream **
412  {
413 
414  real64 const res_enthalpy = m_resEnthalpy[er][esr][ei][0];
415 
416  m_energyPerfFlux[iperf] += flux * res_enthalpy;
417 
418  // energy equation derivatives WRT res P & T
419  m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dP] += dFlux[TAG::RES][CP_Deriv::dP] * res_enthalpy +
420  flux * m_dResEnthalpy[er][esr][ei][0][Deriv::dP];
421  m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dT] += dFlux[TAG::RES][CP_Deriv::dT] * res_enthalpy +
422  flux * m_dResEnthalpy[er][esr][ei][0][Deriv::dT];
423  // energy equation derivatives WRT well P
424  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dP] += dFlux[TAG::WELL][CP_Deriv::dP] * res_enthalpy;
425  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dT] += dFlux[TAG::WELL][CP_Deriv::dT] * res_enthalpy;
426 
427  }
428  else // ** reservoir cell is downstream
429  {
430  real64 const wellelem_enthalpy = m_wellElemEnthalpy[iwelem][0];
431  m_energyPerfFlux[iperf] += flux * wellelem_enthalpy;
432 
433  // energy equation derivatives WRT res P & T
434  m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dP] += dFlux[TAG::RES][CP_Deriv::dP] * wellelem_enthalpy;
435  m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dT] += dFlux[TAG::RES][CP_Deriv::dT] * wellelem_enthalpy;
436 
437  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dP] += dFlux[TAG::WELL][CP_Deriv::dP] * wellelem_enthalpy
438  + flux * m_dWellElemEnthalpy[iwelem][0][Deriv::dP];
439  m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dT] += dFlux[TAG::WELL][CP_Deriv::dT] * wellelem_enthalpy
440  + flux * m_dWellElemEnthalpy[iwelem][0][Deriv::dT];
441  }
442  } );
443 
444  }
445 
446 
447 
456  template< typename POLICY, typename KERNEL_TYPE >
457  static void
458  launch( localIndex const numElements,
459  KERNEL_TYPE const & kernelComponent )
460  {
462  forAll< POLICY >( numElements, [=] GEOS_HOST_DEVICE ( localIndex const iperf )
463  {
464  kernelComponent.computeFlux( iperf );
465 
466  } );
467  }
468 
469 protected:
470 
474 
477  arrayView3d< real64 > const m_dEnergyPerfFlux;
478 
481 
485 
486 };
487 
492 {
493 public:
494 
505  template< typename POLICY >
506  static void
507  createAndLaunch( string const flowSolverName,
508  PerforationData * const perforationData,
509  ElementSubRegionBase const & subRegion,
510  SingleFluidBase const & fluid,
511  ElementRegionManager & elemManager )
512  {
513  integer constexpr IS_THERMAL = 1;
514  using kernelType = PerforationFluxKernel< IS_THERMAL >;
515  typename kernelType::SinglePhaseFlowAccessors singlePhaseFlowAccessors( elemManager, flowSolverName );
516  typename kernelType::SingleFluidAccessors singleFluidAccessors( elemManager, flowSolverName );
517  typename kernelType::ThermalSinglePhaseFlowAccessors thermalSinglePhaseFlowAccessors( elemManager, flowSolverName );
518  typename kernelType::ThermalSingleFluidAccessors thermalSingleFluidAccessors( elemManager, flowSolverName );
519  kernelType kernel( perforationData, subRegion, fluid, singlePhaseFlowAccessors, singleFluidAccessors, thermalSinglePhaseFlowAccessors, thermalSingleFluidAccessors );
520  kernelType::template launch< POLICY >( perforationData->size(), kernel );
521  }
522 };
523 
524 } // end namespace thermalPerforationFluxKernels
525 
526 } // end namespace geos
527 
528 #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:188
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
std::int32_t integer
Signed integer type.
Definition: DataTypes.hpp:81
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:204
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
Definition: DataTypes.hpp:220