GEOS
SinglePhaseWellKernels.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_SINGLEPHASEWELLKERNELS_HPP
21 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_SINGLEPHASEWELLKERNELS_HPP
22 
23 #include "common/DataTypes.hpp"
24 #include "common/GEOS_RAJA_Interface.hpp"
25 #include "physicsSolvers/PhysicsSolverBaseKernels.hpp"
26 
27 #include "constitutive/fluid/singlefluid/SingleFluidFields.hpp"
28 #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp"
29 #include "constitutive/fluid/singlefluid/SingleFluidLayouts.hpp"
30 #include "constitutive/fluid/singlefluid/SingleFluidLayouts.hpp"
31 
32 #include "mesh/ElementRegionManager.hpp"
33 #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp"
34 #include "physicsSolvers/fluidFlow/StencilAccessors.hpp"
35 #include "physicsSolvers/fluidFlow/wells/WellControls.hpp"
36 #include "physicsSolvers/fluidFlow/wells/SinglePhaseWellFields.hpp"
37 
38 #include "physicsSolvers/KernelLaunchSelectors.hpp"
39 
40 namespace geos
41 {
42 
43 namespace singlePhaseWellKernels
44 {
45 
46 // tag to access well and reservoir elements in perforation rates computation
47 struct SubRegionTag
48 {
49  static constexpr integer RES = 0;
50  static constexpr integer WELL = 1;
51 };
52 
53 // tag to access the next and current well elements of a connection
54 struct ElemTag
55 {
56  static constexpr integer CURRENT = 0;
57  static constexpr integer NEXT = 1;
58 };
59 
60 // define the column offset of the derivatives
61 struct ColOffset
62 {
63  static constexpr integer DPRES = 0;
64  static constexpr integer DRATE = 1;
65 };
66 
67 template< integer IS_THERMAL >
68 struct ColOffset_WellJac;
69 
70 template<>
71 struct ColOffset_WellJac< 0 >
72 {
73  static constexpr integer dP = 0;
74  static constexpr integer dQ = dP + 1;
75  static integer constexpr nDer = dQ + 1;
76 
77 };
78 
79 template<>
80 struct ColOffset_WellJac< 1 >
81 {
82  static constexpr integer dP = 0;
83  static constexpr integer dQ = dP + 1;
84  static constexpr integer dT = dQ+1;
86  static integer constexpr nDer = dT + 1;
87 };
88 
89 // define the row offset of the residual equations
90 struct RowOffset
91 {
92  static constexpr integer CONTROL = 0;
93  static constexpr integer MASSBAL = 1;
94 };
95 
96 template< integer IS_THERMAL >
97 struct RowOffset_WellJac;
98 
99 template<>
100 struct RowOffset_WellJac< 0 >
101 {
102  static constexpr integer CONTROL = 0;
103  static constexpr integer MASSBAL = 1;
104  static constexpr integer nEqn = MASSBAL+1;
105 };
106 template<>
107 struct RowOffset_WellJac< 1 >
108 {
109  static constexpr integer CONTROL = 0;
110  static constexpr integer MASSBAL = 1;
111  static constexpr integer ENERGYBAL = MASSBAL+1;
112  static constexpr integer nEqn = ENERGYBAL+1;
113 
114 };
115 /******************************** ControlEquationHelper ********************************/
116 
117 struct ControlEquationHelper
118 {
119 
120  using ROFFSET = singlePhaseWellKernels::RowOffset;
121  using COFFSET = singlePhaseWellKernels::ColOffset;
122 
123  // add an epsilon to the checks to avoid control changes due to tiny pressure/rate updates
124  static constexpr real64 EPS = 1e-15;
125 
126  GEOS_HOST_DEVICE
127  inline
128  static
129  void
130  switchControl( bool const isProducer,
131  WellControls::Control const & currentControl,
132  real64 const & targetBHP,
133  real64 const & targetRate,
134  real64 const & currentBHP,
135  real64 const & currentVolRate,
136  WellControls::Control & newControl );
137 
138  template< integer IS_THERMAL >
139  GEOS_HOST_DEVICE
140  inline
141  static
142  void
143  compute( globalIndex const rankOffset,
144  WellControls::Control const currentControl,
145  real64 const & targetBHP,
146  real64 const & targetRate,
147  real64 const & currentBHP,
148  arrayView1d< real64 const > const & dCurrentBHP,
149  real64 const & currentVolRate,
150  arrayView1d< real64 const > const & dCurrentVolRate,
151  globalIndex const dofNumber,
152  CRSMatrixView< real64, globalIndex const > const & localMatrix,
153  arrayView1d< real64 > const & localRhs );
154 
155 };
156 
157 
158 /******************************** FluxKernel ********************************/
159 
160 struct FluxKernel
161 {
162 
163  using ROFFSET = singlePhaseWellKernels::RowOffset;
164  using COFFSET = singlePhaseWellKernels::ColOffset;
165  using TAG = singlePhaseWellKernels::ElemTag;
166 
167  template< integer IS_THERMAL >
168  static void
169  launch( localIndex const size,
170  globalIndex const rankOffset,
171  arrayView1d< globalIndex const > const & wellElemDofNumber,
172  arrayView1d< localIndex const > const & nextWellElemIndex,
173  arrayView1d< real64 const > const & connRate,
174  real64 const & dt,
175  CRSMatrixView< real64, globalIndex const > const & localMatrix,
176  arrayView1d< real64 > const & localRhs );
177 
178 };
179 
180 
181 /******************************** PressureRelationKernel ********************************/
182 
183 struct PressureRelationKernel
184 {
185 
186  using ROFFSET = singlePhaseWellKernels::RowOffset;
187  using COFFSET = singlePhaseWellKernels::ColOffset;
188  using TAG = singlePhaseWellKernels::ElemTag;
189 
190  template< integer IS_THERMAL >
191  static void
192  launch( localIndex const size,
193  globalIndex const rankOffset,
194  arrayView1d< globalIndex const > const & wellElemDofNumber,
195  arrayView1d< real64 const > const & wellElemGravCoef,
196  arrayView1d< localIndex const > const & nextWellElemIndex,
197  arrayView1d< real64 const > const & wellElemPressure,
198  arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const & wellElemDensity,
199  arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > const & dWellElemDensity,
200  CRSMatrixView< real64, globalIndex const > const & localMatrix,
201  arrayView1d< real64 > const & localRhs );
202 
203 };
204 
205 
206 /******************************** PerforationKernel ********************************/
207 
208 struct PerforationKernel
209 {
210 
211  using TAG = singlePhaseWellKernels::SubRegionTag;
212 
213  using SinglePhaseFlowAccessors =
214  StencilAccessors< fields::flow::pressure >;
215 
216  using SingleFluidAccessors =
217  StencilMaterialAccessors< constitutive::SingleFluidBase,
218  fields::singlefluid::density,
219  fields::singlefluid::dDensity,
220  fields::singlefluid::viscosity,
221  fields::singlefluid::dViscosity >;
222 
230  template< typename VIEWTYPE >
231  using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
232 
233  template< integer IS_THERMAL >
234  GEOS_HOST_DEVICE
235  inline
236  static
237  void
238  compute( real64 const & resPressure,
239  real64 const & resDensity,
240  arraySlice1d< real64 const > const & dResDensity,
241  real64 const & resViscosity,
242  arraySlice1d< real64 const > const & dResViscosity,
243  real64 const & wellElemGravCoef,
244  real64 const & wellElemPressure,
245  real64 const & wellElemDensity,
246  arraySlice1d< real64 const > const & dWellElemDensity,
247  real64 const & wellElemViscosity,
248  arraySlice1d< real64 const > const & dWellElemViscosity,
249  real64 const & perfGravCoef,
250  real64 const & trans,
251  real64 & perfRate,
252  arraySlice2d< real64 > const & dPerfRate,
253  arraySlice1d< real64 > const & dPerfRate_dPres );
254 
255  template< integer IS_THERMAL >
256  static void
257  launch( localIndex const size,
258  ElementViewConst< arrayView1d< real64 const > > const & resPressure,
259  ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > > const & resDensity,
260  ElementViewConst< arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > > const & dResDensity,
261  ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > > const & resViscosity,
262  ElementViewConst< arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > > const & dResViscosity,
263  arrayView1d< real64 const > const & wellElemGravCoef,
264  arrayView1d< real64 const > const & wellElemPressure,
265  arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const & wellElemDensity,
266  arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > const & dWellElemDensity,
267  arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const & wellElemViscosity,
268  arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > const & dWellElemViscosity,
269  arrayView1d< real64 const > const & perfGravCoef,
270  arrayView1d< localIndex const > const & perfWellElemIndex,
271  arrayView1d< real64 const > const & perfTransmissibility,
272  arrayView1d< localIndex const > const & resElementRegion,
273  arrayView1d< localIndex const > const & resElementSubRegion,
274  arrayView1d< localIndex const > const & resElementIndex,
275  arrayView1d< real64 > const & perfRate,
276  arrayView3d< real64 > const & dPerfRate );
277 
278 };
279 
280 /******************************** AccumulationKernel ********************************/
281 
282 struct AccumulationKernel
283 {
284 
285  using ROFFSET = singlePhaseWellKernels::RowOffset;
286  using COFFSET = singlePhaseWellKernels::ColOffset;
287 
288  static void
289  launch( localIndex const size,
290  globalIndex const rankOffset,
291  arrayView1d< globalIndex const > const & wellElemDofNumber,
292  arrayView1d< integer const > const & wellElemGhostRank,
293  arrayView1d< real64 const > const & wellElemVolume,
294  arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const & wellElemDensity,
295  arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > const & dWellElemDensity,
296  arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const & wellElemDensity_n,
297  CRSMatrixView< real64, globalIndex const > const & localMatrix,
298  arrayView1d< real64 > const & localRhs );
299 
300 };
301 
302 /******************************** ElementBasedAssemblyKernel ********************************/
303 
309 template< integer IS_THERMAL >
310 class ElementBasedAssemblyKernel
311 {
312 public:
313 
314  // Well jacobian column and row indicies
315  using FLUID_PROP_COFFSET = constitutive::singlefluid::DerivativeOffsetC< IS_THERMAL >;
316  using WJ_COFFSET = singlePhaseWellKernels::ColOffset_WellJac< IS_THERMAL >;
317  using WJ_ROFFSET = singlePhaseWellKernels::RowOffset_WellJac< IS_THERMAL >;
318 
320  static constexpr integer numDof = 1 + IS_THERMAL; // tjb review
321 
323  static constexpr integer numEqn = 2 + IS_THERMAL;
324 
325 
336  ElementBasedAssemblyKernel( globalIndex const rankOffset,
337  string const dofKey,
338  WellElementSubRegion const & subRegion,
339  constitutive::SingleFluidBase const & fluid,
340  CRSMatrixView< real64, globalIndex const > const & localMatrix,
341  arrayView1d< real64 > const & localRhs )
342  : m_rankOffset( rankOffset ),
343  m_iwelemControl( subRegion.getTopWellElementIndex() ),
344  m_dofNumber( subRegion.getReference< array1d< globalIndex > >( dofKey ) ),
345  m_elemGhostRank( subRegion.ghostRank() ),
346  m_wellElemVolume( subRegion.getElementVolume() ),
347  m_wellElemDensity( fluid.density() ),
348  m_dWellElemDensity( fluid.dDensity() ),
349  m_wellElemDensity_n( fluid.density_n() ),
350  m_localMatrix( localMatrix ),
351  m_localRhs( localRhs )
352  {}
353 
354 
360  GEOS_HOST_DEVICE
361  integer elemGhostRank( localIndex const ei ) const
362  { return m_elemGhostRank( ei ); }
363 
364 
365 
373  template< typename FUNC = NoOpFunc >
374  GEOS_HOST_DEVICE
375  void computeAccumulation( localIndex const iwelem,
376  FUNC && kernelOp = NoOpFunc{} ) const
377  {
378 
379  localIndex const eqnRowIndex = m_dofNumber[iwelem] + WJ_ROFFSET::MASSBAL - m_rankOffset;
380  globalIndex const presDofColIndex = m_dofNumber[iwelem] + WJ_COFFSET::dP;
381 
382  real64 const localAccum = m_wellElemVolume[iwelem] * ( m_wellElemDensity[iwelem][0] - m_wellElemDensity_n[iwelem][0] );
383  real64 const localAccumDP = m_wellElemVolume[iwelem] * m_dWellElemDensity[iwelem][0][FLUID_PROP_COFFSET::dP];
384 
385  // add contribution to global residual and jacobian (no need for atomics here)
386  m_localMatrix.addToRow< serialAtomic >( eqnRowIndex, &presDofColIndex, &localAccumDP, 1 );
387  m_localRhs[eqnRowIndex] += localAccum;
388 
389  if constexpr ( IS_THERMAL )
390  {
391  real64 const localAccumDT = m_wellElemVolume[iwelem] * m_dWellElemDensity[iwelem][0][FLUID_PROP_COFFSET::dT];
392  globalIndex const tempDofColIndex = m_dofNumber[iwelem] + WJ_COFFSET::dT;
393  m_localMatrix.addToRow< serialAtomic >( eqnRowIndex, &tempDofColIndex, &localAccumDT, 1 );
394  kernelOp( presDofColIndex, tempDofColIndex );
395  }
396 
397  }
398 
399 
400 
408  template< typename POLICY, typename KERNEL_TYPE >
409  static void
410  launch( localIndex const numElems,
411  KERNEL_TYPE const & kernelComponent )
412  {
413  GEOS_MARK_FUNCTION;
414 
415  forAll< POLICY >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const iwelem )
416  {
417  if( kernelComponent.elemGhostRank( iwelem ) >= 0 )
418  {
419  return;
420  }
421  kernelComponent.computeAccumulation( iwelem );
422  } );
423  }
424 
425 protected:
426 
428  globalIndex const m_rankOffset;
429 
431  localIndex const m_iwelemControl;
432 
434  arrayView1d< globalIndex const > const m_dofNumber;
435 
437  arrayView1d< integer const > const m_elemGhostRank;
438 
440  arrayView1d< real64 const > const m_wellElemVolume;
441 
443  arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_wellElemDensity;
444  arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > const m_dWellElemDensity;
445  arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_wellElemDensity_n;
446 
448  CRSMatrixView< real64, globalIndex const > const m_localMatrix;
450  arrayView1d< real64 > const m_localRhs;
451 
452 };
453 
457 class ElementBasedAssemblyKernelFactory
458 {
459 public:
470  template< typename POLICY >
471  static void
472  createAndLaunch( globalIndex const rankOffset,
473  string const dofKey,
474  WellElementSubRegion const & subRegion,
475  constitutive::SingleFluidBase const & fluid,
476  CRSMatrixView< real64, globalIndex const > const & localMatrix,
477  arrayView1d< real64 > const & localRhs )
478  {
479  integer constexpr isThermal=0;
480 
481  ElementBasedAssemblyKernel< isThermal >
482  kernel( rankOffset, dofKey, subRegion, fluid, localMatrix, localRhs );
483  ElementBasedAssemblyKernel< isThermal >::template
484  launch< POLICY, ElementBasedAssemblyKernel< isThermal > >( subRegion.size(), kernel );
485  }
486 };
487 
488 /******************************** PressureTemperatyrInitializationKernel ********************************/
489 
490 // tjb make this templated on thermal
491 struct PresTempInitializationKernel
492 {
493 
494  using SinglePhaseFlowAccessors =
495  StencilAccessors< fields::flow::pressure,
496  fields::flow::temperature >;
497 
498  using SingleFluidAccessors =
499  StencilMaterialAccessors< constitutive::SingleFluidBase,
500  fields::singlefluid::density >;
501 
509  template< typename VIEWTYPE >
510  using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
511 
512  static void
513  launch( integer const isThermal,
514  localIndex const perforationSize,
515  localIndex const subRegionSize,
516  localIndex const numPerforations,
517  WellControls const & wellControls,
518  real64 const & currentTime,
519  ElementViewConst< arrayView1d< real64 const > > const & resPressure,
520  ElementViewConst< arrayView1d< real64 const > > const & resTemperature,
521  ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > > const & resDensity,
522  arrayView1d< localIndex const > const & resElementRegion,
523  arrayView1d< localIndex const > const & resElementSubRegion,
524  arrayView1d< localIndex const > const & resElementIndex,
525  arrayView1d< real64 const > const & perfGravCoef,
526  arrayView1d< real64 const > const & wellElemGravCoef,
527  arrayView1d< real64 > const & wellElemPressure,
528  arrayView1d< real64 > const & wellElemTemperature );
529 
530 };
531 
532 /******************************** FaceBasedAssemblyKernel ********************************/
540 template< integer IS_THERMAL >
541 class FaceBasedAssemblyKernel
542 {
543 public:
544 
545  using COFFSET = singlePhaseWellKernels::ColOffset;
546  using ROFFSET = singlePhaseWellKernels::RowOffset;
547  using TAG = singlePhaseWellKernels::ElemTag;
548 
549  using FLUID_PROP_COFFSET = constitutive::singlefluid::DerivativeOffsetC< IS_THERMAL >;
550  using WJ_COFFSET = singlePhaseWellKernels::ColOffset_WellJac< IS_THERMAL >;
551  using WJ_ROFFSET = singlePhaseWellKernels::RowOffset_WellJac< IS_THERMAL >;
552 
553  using CP_Deriv = constitutive::singlefluid::DerivativeOffsetC< IS_THERMAL >;
554 
556  static constexpr integer numDof = WJ_COFFSET::nDer; // tjb revisit
557 
559  static constexpr integer numEqn = WJ_COFFSET::nDer;// tjb revisit
560 
561  static constexpr integer maxNumElems = 2;
562  static constexpr integer maxStencilSize = 2;
573  FaceBasedAssemblyKernel( real64 const dt,
574  globalIndex const rankOffset,
575  string const wellDofKey,
576  WellControls const & wellControls,
577  WellElementSubRegion const & subRegion,
578  CRSMatrixView< real64, globalIndex const > const & localMatrix,
579  arrayView1d< real64 > const & localRhs )
580  :
581  m_dt( dt ),
582  m_rankOffset( rankOffset ),
583  m_wellElemDofNumber ( subRegion.getReference< array1d< globalIndex > >( wellDofKey ) ),
584  m_nextWellElemIndex ( subRegion.getReference< array1d< localIndex > >( WellElementSubRegion::viewKeyStruct::nextWellElementIndexString()) ),
585  m_connRate ( subRegion.getField< fields::well::connectionRate >() ),
586  m_localMatrix( localMatrix ),
587  m_localRhs ( localRhs ),
588  m_isProducer ( wellControls.isProducer() )
589  {}
590 
591 
599  template< typename FUNC = NoOpFunc >
600  GEOS_HOST_DEVICE
601  inline
602  void computeFlux( localIndex const iwelem,
603  FUNC && compFluxKernelOp = NoOpFunc{} ) const
604  {
605  // 1) Compute the flux and its derivatives
606 
607  /* currentConnRate < 0 flow from iwelem to iwelemNext
608  * currentConnRate > 0 flow from iwelemNext to iwelem
609  * With this convention, currentConnRate < 0 at the last connection for a producer
610  * currentConnRate > 0 at the last connection for a injector
611  */
612 
613  // get next well element index
614  localIndex const iwelemNext = m_nextWellElemIndex[iwelem];
615 
616  // there is nothing to upwind for single-phase flow
617  real64 const currentConnRate = m_connRate[iwelem];
618  real64 const flux = m_dt * currentConnRate;
619  real64 const dFlux_dRate = m_dt;
620 
621  // 2) Assemble the flux into residual and Jacobian
622  if( iwelemNext < 0 )
623  {
624  // flux terms
625  real64 const oneSidedLocalFlux = -flux;
626  real64 const oneSidedLocalFluxJacobian_dRate = -dFlux_dRate;
627 
628  // jacobian indices
629  globalIndex const oneSidedEqnRowIndex = m_wellElemDofNumber[iwelem] + ROFFSET::MASSBAL - m_rankOffset;
630  globalIndex const oneSidedDofColIndex_dRate = m_wellElemDofNumber[iwelem] + COFFSET::DRATE;
631 
632  if( oneSidedEqnRowIndex >= 0 && oneSidedEqnRowIndex < m_localMatrix.numRows() )
633  {
634  m_localMatrix.addToRow< parallelDeviceAtomic >( oneSidedEqnRowIndex,
635  &oneSidedDofColIndex_dRate,
636  &oneSidedLocalFluxJacobian_dRate,
637  1 );
638  RAJA::atomicAdd( parallelDeviceAtomic{}, &m_localRhs[oneSidedEqnRowIndex], oneSidedLocalFlux );
639  }
640  }
641  else
642  {
643  // local working variables and arrays
644  globalIndex eqnRowIndices[2]{};
645 
646  real64 localFlux[2]{};
647  real64 localFluxJacobian_dRate[2]{};
648 
649  // flux terms
650  localFlux[TAG::NEXT] = flux;
651  localFlux[TAG::CURRENT] = -flux;
652 
653  localFluxJacobian_dRate[TAG::NEXT] = dFlux_dRate;
654  localFluxJacobian_dRate[TAG::CURRENT] = -dFlux_dRate;
655 
656  // indices
657  eqnRowIndices[TAG::CURRENT] = m_wellElemDofNumber[iwelem] + ROFFSET::MASSBAL - m_rankOffset;
658  eqnRowIndices[TAG::NEXT] = m_wellElemDofNumber[iwelemNext] + ROFFSET::MASSBAL - m_rankOffset;
659  globalIndex const dofColIndex_dRate = m_wellElemDofNumber[iwelem] + COFFSET::DRATE;
660 
661  for( localIndex i = 0; i < 2; ++i )
662  {
663  if( eqnRowIndices[i] >= 0 && eqnRowIndices[i] < m_localMatrix.numRows() )
664  {
665  m_localMatrix.addToRow< parallelDeviceAtomic >( eqnRowIndices[i],
666  &dofColIndex_dRate,
667  &localFluxJacobian_dRate[i],
668  1 );
669  RAJA::atomicAdd( parallelDeviceAtomic{}, &m_localRhs[eqnRowIndices[i]], localFlux[i] );
670  }
671  }
672  }
673  compFluxKernelOp( iwelemNext, currentConnRate );
674 
675  }
676 
677 
685  template< typename POLICY, typename KERNEL_TYPE >
686  static void
687  launch( localIndex const numElements,
688  KERNEL_TYPE const & kernelComponent )
689  {
690  GEOS_MARK_FUNCTION;
691  forAll< POLICY >( numElements, [=] GEOS_HOST_DEVICE ( localIndex const ie )
692  {
693  kernelComponent.computeFlux( ie );
694  } );
695  }
696 
697 protected:
699  real64 const m_dt;
701  integer const m_rankOffset;
702 
704  arrayView1d< globalIndex const > const m_wellElemDofNumber;
706  arrayView1d< localIndex const > const m_nextWellElemIndex;
707 
709  arrayView1d< real64 const > const m_connRate;
710 
712  arrayView2d< real64 const, compflow::USD_COMP > const m_wellElemCompFrac;
713 
715  CRSMatrixView< real64, globalIndex const > const m_localMatrix;
717  arrayView1d< real64 > const m_localRhs;
718 
720  bool const m_isProducer;
721 
722 };
723 
727 class FaceBasedAssemblyKernelFactory
728 {
729 public:
730 
742  template< typename POLICY >
743  static void
744  createAndLaunch( real64 const dt,
745  globalIndex const rankOffset,
746  string const dofKey,
747  WellControls const & wellControls,
748  WellElementSubRegion const & subRegion,
749  CRSMatrixView< real64, globalIndex const > const & localMatrix,
750  arrayView1d< real64 > const & localRhs )
751  {
752  integer isThermal=0;
753  geos::internal::kernelLaunchSelectorThermalSwitch( isThermal, [&]( auto IS_THERMAL )
754  {
755 
756  integer constexpr istherm = IS_THERMAL();
757 
758  using kernelType = FaceBasedAssemblyKernel< istherm >;
759  kernelType kernel( dt, rankOffset, dofKey, wellControls, subRegion, localMatrix, localRhs );
760  kernelType::template launch< POLICY >( subRegion.size(), kernel );
761  } );
762  }
763 };
764 /******************************** RateInitializationKernel ********************************/
765 
766 struct RateInitializationKernel
767 {
768 
769  static void
770  launch( localIndex const subRegionSize,
771  WellControls const & wellControls,
772  real64 const & currentTime,
773  arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const & wellElemDens,
774  arrayView1d< real64 > const & connRate );
775 
776 };
777 
778 
779 /******************************** ResidualNormKernel ********************************/
780 
784 class ResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBase< 1 >
785 {
786 public:
787 
788  using Base = physicsSolverBaseKernels::ResidualNormKernelBase< 1 >;
789  using Base::m_minNormalizer;
790  using Base::m_rankOffset;
791  using Base::m_localResidual;
792  using Base::m_dofNumber;
793 
794  ResidualNormKernel( globalIndex const rankOffset,
795  arrayView1d< real64 const > const & localResidual,
796  arrayView1d< globalIndex const > const & dofNumber,
797  arrayView1d< localIndex const > const & ghostRank,
798  WellElementSubRegion const & subRegion,
799  constitutive::SingleFluidBase const & fluid,
800  WellControls const & wellControls,
801  real64 const time,
802  real64 const dt,
803  real64 const minNormalizer )
804  : Base( rankOffset,
805  localResidual,
806  dofNumber,
807  ghostRank,
808  minNormalizer ),
809  m_dt( dt ),
810  m_isLocallyOwned( subRegion.isLocallyOwned() ),
811  m_iwelemControl( subRegion.getTopWellElementIndex() ),
812  m_currentControl( wellControls.getControl() ),
813  m_constraintValue ( wellControls.getCurrentConstraint()->getConstraintValue( time )),
814  m_volume( subRegion.getElementVolume() ),
815  m_density_n( fluid.density_n() )
816  {
817  if( wellControls.isProducer() )
818  {
819  m_targetBHP = wellControls.getMinBHPConstraint()->getConstraintValue( time );
820  }
821  else
822  {
823  m_targetBHP = wellControls.getMaxBHPConstraint()->getConstraintValue( time );
824  }
825  }
826 
827  GEOS_HOST_DEVICE
828  virtual void computeLinf( localIndex const iwelem,
829  LinfStackVariables & stack ) const override
830  {
831  for( localIndex idof = 0; idof < 2; ++idof )
832  {
833  real64 normalizer = 0.0;
834  if( idof == singlePhaseWellKernels::RowOffset::CONTROL )
835  {
836  // for the top well element, normalize using the current control
837  if( m_isLocallyOwned && iwelem == m_iwelemControl )
838  {
839  if( m_currentControl == WellControls::Control::BHP )
840  {
841  // this residual entry is in pressure units
842  normalizer = m_targetBHP;
843  }
844  else if( m_currentControl == WellControls::Control::TOTALVOLRATE )
845  {
846  // this residual entry is in volume / time units
847  normalizer = LvArray::math::max( LvArray::math::abs( m_constraintValue ), m_minNormalizer );
848  }
849  else if( m_currentControl == WellControls::Control::MASSRATE )
850  {
851  // the residual entry is in volume / time units
852  normalizer = LvArray::math::max( LvArray::math::abs( m_constraintValue ), m_minNormalizer );
853  }
854  }
855  // for the pressure difference equation, always normalize by the BHP
856  else
857  {
858  // this residual is in pressure units
859  normalizer = m_targetBHP;
860  }
861  }
862  else // SinglePhaseWell::RowOffset::MASSBAL
863  {
864  // this residual entry is in mass units
865  normalizer = m_dt * LvArray::math::abs( m_constraintValue ) * m_density_n[iwelem][0];
866 
867  // to make sure that everything still works well if the rate is zero, we add this check
868  normalizer = LvArray::math::max( normalizer, m_volume[iwelem] * m_density_n[iwelem][0] );
869  }
870 
871  // we have the normalizer now, we can compute a dimensionless Linfty norm contribution
872  real64 const val = LvArray::math::abs( m_localResidual[stack.localRow + idof] ) / normalizer;
873  if( val > stack.localValue[0] )
874  {
875  stack.localValue[0] = val;
876  }
877 
878  }
879  }
880 
881  GEOS_HOST_DEVICE
882  virtual void computeL2( localIndex const iwelem,
883  L2StackVariables & stack ) const override
884  {
885  GEOS_UNUSED_VAR( iwelem, stack );
886  GEOS_ERROR( "The L2 norm is not implemented for SinglePhaseWell" );
887  }
888 
889 
890 protected:
891 
893  real64 const m_dt;
894 
896  bool const m_isLocallyOwned;
897 
899  localIndex const m_iwelemControl;
900 
902  WellControls::Control const m_currentControl;
903  real64 const m_constraintValue;
904  real64 m_targetBHP;
905 
906 
908  arrayView1d< real64 const > const m_volume;
909 
911  arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_density_n;
912 
913 };
914 
918 class ResidualNormKernelFactory
919 {
920 public:
921 
935  template< typename POLICY >
936  static void
937  createAndLaunch( globalIndex const rankOffset,
938  string const & dofKey,
939  arrayView1d< real64 const > const & localResidual,
940  WellElementSubRegion const & subRegion,
941  constitutive::SingleFluidBase const & fluid,
942  WellControls const & wellControls,
943  real64 const time,
944  real64 const dt,
945  real64 const minNormalizer,
946  real64 (& residualNorm)[1] )
947  {
948  arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey );
949  arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
950 
951  ResidualNormKernel kernel( rankOffset, localResidual, dofNumber, ghostRank, subRegion, fluid, wellControls, time, dt, minNormalizer );
952  ResidualNormKernel::launchLinf< POLICY >( subRegion.size(), kernel, residualNorm );
953  }
954 
955 };
956 
957 } // end namespace singlePhaseWellKernels
958 
959 } // end namespace geos
960 
961 #endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_SINGLEPHASEWELLKERNELS_HPP
GEOS_HOST_DEVICE
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
GEOS_UNUSED_VAR
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
Definition: GeosxMacros.hpp:84
GEOS_ERROR
#define GEOS_ERROR(msg)
Raise a hard error and terminate the program.
Definition: Logger.hpp:157
GEOS_MARK_FUNCTION
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
Definition: TimingMacros.hpp:143
geos::ElementRegionManager::ElementViewConst
typename ElementViewAccessor< VIEWTYPE >::NestedViewTypeConst ElementViewConst
The ElementViewAccessor at the ElementRegionManager level is the type resulting from ElementViewAcces...
Definition: ElementRegionManager.hpp:83
geos::ElementSubRegionBase::getElementVolume
arrayView1d< real64 const > getElementVolume() const
Get the volume of each element in this subregion.
Definition: ElementSubRegionBase.hpp:160
geos::ObjectManagerBase::ghostRank
array1d< integer > const & ghostRank()
Get the ghost information of each object.
Definition: ObjectManagerBase.hpp:856
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::WellConstraintBase::getConstraintValue
real64 getConstraintValue(real64 const &currentTime) const
Get the target bottom hole pressure value.
Definition: WellConstraintsBase.hpp:156
geos::WellControls
This class describes the controls used to operate a well.
Definition: WellControls.hpp:54
geos::WellControls::isProducer
bool isProducer() const
Is the well a producer?
Definition: WellControls.hpp:356
geos::WellControls::getControl
Control getControl() const
Get the control type for the well.
Definition: WellControls.hpp:262
geos::WellControls::getMinBHPConstraint
MinimumBHPConstraint * getMinBHPConstraint()
Getters for constraints.
Definition: WellControls.hpp:521
geos::WellElementSubRegion
This class describes a collection of local well elements and perforations.
Definition: WellElementSubRegion.hpp:32
geos::WellElementSubRegion::isLocallyOwned
bool isLocallyOwned() const
Check if well is owned by current rank.
geos::WellElementSubRegion::getTopWellElementIndex
localIndex getTopWellElementIndex() const
Get for the top element index.
Definition: WellElementSubRegion.hpp:134
geos::dataRepository::Group::getReference
GEOS_DECLTYPE_AUTO_RETURN getReference(LOOKUP_TYPE const &lookup) const
Look up a wrapper and get reference to wrapped object.
Definition: Group.hpp:1275
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::physicsSolverBaseKernels::ResidualNormKernelBase
Define the base interface for the residual calculations.
Definition: PhysicsSolverBaseKernels.hpp:42
geos::physicsSolverBaseKernels::ResidualNormKernelBase< 1 >::m_minNormalizer
real64 const m_minNormalizer
Value used to make sure that normalizers are never zero.
Definition: PhysicsSolverBaseKernels.hpp:236
geos::physicsSolverBaseKernels::ResidualNormKernelBase< 1 >::m_dofNumber
arrayView1d< globalIndex const > const m_dofNumber
View on the dof numbers.
Definition: PhysicsSolverBaseKernels.hpp:230
geos::physicsSolverBaseKernels::ResidualNormKernelBase< 1 >::ghostRank
GEOS_HOST_DEVICE integer ghostRank(localIndex const i) const
Getter for the ghost rank.
Definition: PhysicsSolverBaseKernels.hpp:91
geos::physicsSolverBaseKernels::ResidualNormKernelBase< 1 >::m_localResidual
arrayView1d< real64 const > const m_localResidual
View on the local residual.
Definition: PhysicsSolverBaseKernels.hpp:227
geos::singlePhaseWellKernels::ElementBasedAssemblyKernelFactory::createAndLaunch
static void createAndLaunch(globalIndex const rankOffset, string const dofKey, WellElementSubRegion const &subRegion, constitutive::SingleFluidBase const &fluid, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs)
Create a new kernel and launch.
Definition: SinglePhaseWellKernels.hpp:472
geos::singlePhaseWellKernels::ElementBasedAssemblyKernel
Define the interface for the assembly kernel in charge of accumulation.
Definition: SinglePhaseWellKernels.hpp:311
geos::singlePhaseWellKernels::ElementBasedAssemblyKernel::m_wellElemVolume
arrayView1d< real64 const > const m_wellElemVolume
View on the element volumes.
Definition: SinglePhaseWellKernels.hpp:440
geos::singlePhaseWellKernels::ElementBasedAssemblyKernel::elemGhostRank
GEOS_HOST_DEVICE integer elemGhostRank(localIndex const ei) const
Getter for the ghost rank of an element.
Definition: SinglePhaseWellKernels.hpp:361
geos::singlePhaseWellKernels::ElementBasedAssemblyKernel::m_localRhs
arrayView1d< real64 > const m_localRhs
View on the local RHS.
Definition: SinglePhaseWellKernels.hpp:450
geos::singlePhaseWellKernels::ElementBasedAssemblyKernel::ElementBasedAssemblyKernel
ElementBasedAssemblyKernel(globalIndex const rankOffset, string const dofKey, WellElementSubRegion const &subRegion, constitutive::SingleFluidBase const &fluid, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs)
Constructor.
Definition: SinglePhaseWellKernels.hpp:336
geos::singlePhaseWellKernels::ElementBasedAssemblyKernel::m_wellElemDensity
arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_wellElemDensity
Views on the density.
Definition: SinglePhaseWellKernels.hpp:443
geos::singlePhaseWellKernels::ElementBasedAssemblyKernel::numEqn
static constexpr integer numEqn
Compute time value for the number of equations mass bal + momentum + energy bal.
Definition: SinglePhaseWellKernels.hpp:323
geos::singlePhaseWellKernels::ElementBasedAssemblyKernel::computeAccumulation
GEOS_HOST_DEVICE void computeAccumulation(localIndex const iwelem, FUNC &&kernelOp=NoOpFunc{}) const
Compute the local accumulation contributions to the residual and Jacobian.
Definition: SinglePhaseWellKernels.hpp:375
geos::singlePhaseWellKernels::ElementBasedAssemblyKernel::numDof
static constexpr integer numDof
Number of Dof's set in this kernal - no dQ in accum.
Definition: SinglePhaseWellKernels.hpp:320
geos::singlePhaseWellKernels::ElementBasedAssemblyKernel::m_localMatrix
CRSMatrixView< real64, globalIndex const > const m_localMatrix
View on the local CRS matrix.
Definition: SinglePhaseWellKernels.hpp:448
geos::singlePhaseWellKernels::ElementBasedAssemblyKernel::m_dofNumber
arrayView1d< globalIndex const > const m_dofNumber
View on the dof numbers.
Definition: SinglePhaseWellKernels.hpp:434
geos::singlePhaseWellKernels::ElementBasedAssemblyKernel::launch
static void launch(localIndex const numElems, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
Definition: SinglePhaseWellKernels.hpp:410
geos::singlePhaseWellKernels::ElementBasedAssemblyKernel::m_elemGhostRank
arrayView1d< integer const > const m_elemGhostRank
View on the ghost ranks.
Definition: SinglePhaseWellKernels.hpp:437
geos::singlePhaseWellKernels::ElementBasedAssemblyKernel::m_iwelemControl
localIndex const m_iwelemControl
Index of the element where the control is enforced.
Definition: SinglePhaseWellKernels.hpp:431
geos::singlePhaseWellKernels::FaceBasedAssemblyKernelFactory::createAndLaunch
static void createAndLaunch(real64 const dt, globalIndex const rankOffset, string const dofKey, WellControls const &wellControls, WellElementSubRegion const &subRegion, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs)
Create a new kernel and launch.
Definition: SinglePhaseWellKernels.hpp:744
geos::singlePhaseWellKernels::FaceBasedAssemblyKernel
Define the interface for the assembly kernel in charge of flux terms.
Definition: SinglePhaseWellKernels.hpp:542
geos::singlePhaseWellKernels::FaceBasedAssemblyKernel::computeFlux
GEOS_HOST_DEVICE void computeFlux(localIndex const iwelem, FUNC &&compFluxKernelOp=NoOpFunc{}) const
Compute the local flux contributions to the residual and Jacobian.
Definition: SinglePhaseWellKernels.hpp:602
geos::singlePhaseWellKernels::FaceBasedAssemblyKernel::m_connRate
arrayView1d< real64 const > const m_connRate
Connection rate.
Definition: SinglePhaseWellKernels.hpp:709
geos::singlePhaseWellKernels::FaceBasedAssemblyKernel::FaceBasedAssemblyKernel
FaceBasedAssemblyKernel(real64 const dt, globalIndex const rankOffset, string const wellDofKey, WellControls const &wellControls, WellElementSubRegion const &subRegion, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs)
Constructor for the kernel interface.
Definition: SinglePhaseWellKernels.hpp:573
geos::singlePhaseWellKernels::FaceBasedAssemblyKernel::m_wellElemCompFrac
arrayView2d< real64 const, compflow::USD_COMP > const m_wellElemCompFrac
Element component fraction.
Definition: SinglePhaseWellKernels.hpp:712
geos::singlePhaseWellKernels::FaceBasedAssemblyKernel::m_localMatrix
CRSMatrixView< real64, globalIndex const > const m_localMatrix
View on the local CRS matrix.
Definition: SinglePhaseWellKernels.hpp:715
geos::singlePhaseWellKernels::FaceBasedAssemblyKernel::m_localRhs
arrayView1d< real64 > const m_localRhs
View on the local RHS.
Definition: SinglePhaseWellKernels.hpp:717
geos::singlePhaseWellKernels::FaceBasedAssemblyKernel::m_rankOffset
integer const m_rankOffset
Rank offset for calculating row/col Jacobian indices.
Definition: SinglePhaseWellKernels.hpp:701
geos::singlePhaseWellKernels::FaceBasedAssemblyKernel::m_nextWellElemIndex
arrayView1d< localIndex const > const m_nextWellElemIndex
Next element index, needed since iterating over element nodes, not edges.
Definition: SinglePhaseWellKernels.hpp:706
geos::singlePhaseWellKernels::FaceBasedAssemblyKernel::numDof
static constexpr integer numDof
Number of Dof's set in this kernal.
Definition: SinglePhaseWellKernels.hpp:556
geos::singlePhaseWellKernels::FaceBasedAssemblyKernel::launch
static void launch(localIndex const numElements, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
Definition: SinglePhaseWellKernels.hpp:687
geos::singlePhaseWellKernels::FaceBasedAssemblyKernel::numEqn
static constexpr integer numEqn
Compile time value for the number of equations except rate, momentum, energy.
Definition: SinglePhaseWellKernels.hpp:559
geos::singlePhaseWellKernels::FaceBasedAssemblyKernel::m_wellElemDofNumber
arrayView1d< globalIndex const > const m_wellElemDofNumber
Reference to the degree-of-freedom numbers.
Definition: SinglePhaseWellKernels.hpp:704
geos::singlePhaseWellKernels::ResidualNormKernelFactory::createAndLaunch
static void createAndLaunch(globalIndex const rankOffset, string const &dofKey, arrayView1d< real64 const > const &localResidual, WellElementSubRegion const &subRegion, constitutive::SingleFluidBase const &fluid, WellControls const &wellControls, real64 const time, real64 const dt, real64 const minNormalizer, real64(&residualNorm)[1])
Create a new kernel and launch.
Definition: SinglePhaseWellKernels.hpp:937
geos::singlePhaseWellKernels::ResidualNormKernel::m_isLocallyOwned
bool const m_isLocallyOwned
Flag indicating whether the well is locally owned or not.
Definition: SinglePhaseWellKernels.hpp:896
geos::singlePhaseWellKernels::ResidualNormKernel::m_volume
arrayView1d< real64 const > const m_volume
View on the volume.
Definition: SinglePhaseWellKernels.hpp:908
geos::singlePhaseWellKernels::ResidualNormKernel::m_density_n
arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_density_n
View on total density at the previous converged time step.
Definition: SinglePhaseWellKernels.hpp:911
geos::singlePhaseWellKernels::ResidualNormKernel::computeL2
virtual GEOS_HOST_DEVICE void computeL2(localIndex const iwelem, L2StackVariables &stack) const override
Compute the local values and normalizer for the L2 norm.
Definition: SinglePhaseWellKernels.hpp:882
geos::singlePhaseWellKernels::ResidualNormKernel::computeLinf
virtual GEOS_HOST_DEVICE void computeLinf(localIndex const iwelem, LinfStackVariables &stack) const override
Compute the local values for the Linf norm.
Definition: SinglePhaseWellKernels.hpp:828
geos::singlePhaseWellKernels::ResidualNormKernel::m_iwelemControl
localIndex const m_iwelemControl
Index of the element where the control is enforced.
Definition: SinglePhaseWellKernels.hpp:899
geos::singlePhaseWellKernels::ResidualNormKernel::m_currentControl
WellControls::Control const m_currentControl
Controls.
Definition: SinglePhaseWellKernels.hpp:902
geos::arrayView1d
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:179
geos::globalIndex
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
Definition: DataTypes.hpp:87
geos::arraySlice2d
ArraySlice< T, 2, USD > arraySlice2d
Alias for 2D array slice.
Definition: DataTypes.hpp:199
geos::real64
double real64
64-bit floating point type.
Definition: DataTypes.hpp:98
geos::localIndex
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:84
geos::arraySlice1d
ArraySlice< T, 1, USD > arraySlice1d
Alias for 1D array slice.
Definition: DataTypes.hpp:183
geos::CRSMatrixView
LvArray::CRSMatrixView< T, COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
Definition: DataTypes.hpp:309
geos::arrayView2d
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:195
geos::integer
int integer
Signed integer type.
Definition: DataTypes.hpp:81
geos::array1d
Array< T, 1 > array1d
Alias for 1D array.
Definition: DataTypes.hpp:175
geos::arrayView3d
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
Definition: DataTypes.hpp:211
geos::singlePhaseWellKernels::ColOffset_WellJac
Definition: SinglePhaseWellKernels.hpp:68
geos::singlePhaseWellKernels::PerforationKernel::ElementViewConst
ElementRegionManager::ElementViewConst< VIEWTYPE > ElementViewConst
The type for element-based non-constitutive data parameters. Consists entirely of ArrayView's.
Definition: SinglePhaseWellKernels.hpp:231
geos::singlePhaseWellKernels::PresTempInitializationKernel::ElementViewConst
ElementRegionManager::ElementViewConst< VIEWTYPE > ElementViewConst
The type for element-based non-constitutive data parameters. Consists entirely of ArrayView's.
Definition: SinglePhaseWellKernels.hpp:510
geos::singlePhaseWellKernels::RowOffset_WellJac
Definition: SinglePhaseWellKernels.hpp:97
Generated by doxygen 1.9.1