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