GEOS
CoupledReservoirAndWellKernels.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_MULTIPHYSICS_COUPLEDRESERVOIRANDWELLS_HPP
21 #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_COUPLEDRESERVOIRANDWELLS_HPP
22 
23 
24 #include "common/DataTypes.hpp"
25 #include "common/GEOS_RAJA_Interface.hpp"
26 #include "constitutive/fluid/multifluid/Layouts.hpp"
31 namespace geos
32 {
33 
34 namespace coupledReservoirAndWellKernels
35 {
36 
37 using namespace constitutive;
38 
44 template< integer NC, integer IS_THERMAL >
46 {
47 public:
48 
50  static constexpr integer numComp = NC;
51  static constexpr integer resNumDOF = NC+1+IS_THERMAL;
52 
53  // Well jacobian column and row indicies
56 
59 
60  using CP_Deriv = multifluid::DerivativeOffsetC< NC, IS_THERMAL >;
61 
63 
64 
65 
67  static constexpr integer numDof = WJ_COFFSET::nDer;
68 
70  static constexpr integer numEqn = WJ_ROFFSET::nEqn - 2;
71 
87  globalIndex const rankOffset,
88  string const wellDofKey,
89  WellElementSubRegion const & subRegion,
91  PerforationData const * const perforationData,
92  MultiFluidBase const & fluid,
93 
94  arrayView1d< real64 > const & localRhs,
96  bool const & detectCrossflow,
97  integer & numCrossFlowPerforations,
98  BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags )
99  :
100  m_dt( dt ),
101  m_numPhases ( fluid.numFluidPhases()),
102  m_rankOffset( rankOffset ),
103  m_compPerfRate( perforationData->getField< fields::well::compPerforationRate >() ),
104  m_dCompPerfRate( perforationData->getField< fields::well::dCompPerforationRate >() ),
105  m_perfWellElemIndex( perforationData->getField< fields::perforation::wellElementIndex >() ),
106  m_perfStatus( perforationData->getField< fields::perforation::perforationStatus >() ),
107  m_wellElemDofNumber( subRegion.getReference< array1d< globalIndex > >( wellDofKey ) ),
108  m_resElemDofNumber( resDofNumber ),
109  m_resElementRegion( perforationData->getField< fields::perforation::reservoirElementRegion >() ),
110  m_resElementSubRegion( perforationData->getField< fields::perforation::reservoirElementSubRegion >() ),
111  m_resElementIndex( perforationData->getField< fields::perforation::reservoirElementIndex >() ),
112  m_localRhs( localRhs ),
113  m_localMatrix( localMatrix ),
114  m_detectCrossflow( detectCrossflow ),
115  m_numCrossFlowPerforations( numCrossFlowPerforations ),
116  m_useTotalMassEquation ( kernelFlags.isSet( isothermalCompositionalMultiphaseBaseKernels::KernelFlags::TotalMassEquation ) )
117  { }
118 
119 
128  template< typename FUNC = NoOpFunc >
130  inline
131  void computeFlux( localIndex const iperf,
132  FUNC && compFluxKernelOp = NoOpFunc{} ) const
133  {
134 
135  using namespace compositionalMultiphaseUtilities;
136  if( m_perfStatus[iperf ] )
137  {
138  // local working variables and arrays
139  stackArray1d< localIndex, 2* numComp > eqnRowIndices( 2 * numComp );
140  stackArray1d< globalIndex, 2*resNumDOF > dofColIndices( 2 * resNumDOF );
141 
142  stackArray1d< real64, 2 * numComp > localPerf( 2 * numComp );
143  stackArray2d< real64, 2 * resNumDOF * 2 * numComp > localPerfJacobian( 2 * numComp, 2 * resNumDOF );
144 
145  // get the reservoir (sub)region and element indices
146  localIndex const er = m_resElementRegion[iperf];
147  localIndex const esr = m_resElementSubRegion[iperf];
148  localIndex const ei = m_resElementIndex[iperf];
149 
150  // get the well element index for this perforation
151  localIndex const iwelem = m_perfWellElemIndex[iperf];
152  globalIndex const resOffset = m_resElemDofNumber[er][esr][ei];
153  globalIndex const wellElemOffset = m_wellElemDofNumber[iwelem];
154 
155  for( integer ic = 0; ic < numComp; ++ic )
156  {
157  eqnRowIndices[TAG::RES * numComp + ic] = LvArray::integerConversion< localIndex >( resOffset - m_rankOffset ) + ic;
158  eqnRowIndices[TAG::WELL * numComp + ic] = LvArray::integerConversion< localIndex >( wellElemOffset - m_rankOffset ) + WJ_ROFFSET::MASSBAL + ic;
159  }
160  // Note res and well have same col lineup for P and compdens
161  for( integer jdof = 0; jdof < NC+1; ++jdof )
162  {
163  dofColIndices[TAG::RES * resNumDOF + jdof] = resOffset + jdof;
164  dofColIndices[TAG::WELL * resNumDOF + jdof] = wellElemOffset + WJ_COFFSET::dP + jdof;
165  }
166  // For temp its different
167  if constexpr ( IS_THERMAL )
168  {
169  dofColIndices[TAG::RES * resNumDOF + NC+1 ] = resOffset + NC+1;
170  dofColIndices[TAG::WELL * resNumDOF + NC+1 ] = wellElemOffset + WJ_COFFSET::dT;
171  }
172  // populate local flux vector and derivatives
173  for( integer ic = 0; ic < numComp; ++ic )
174  {
175  localPerf[TAG::RES * numComp + ic] = m_dt * m_compPerfRate[iperf][ic];
176  localPerf[TAG::WELL * numComp + ic] = -m_dt * m_compPerfRate[iperf][ic];
177 
178  if( m_detectCrossflow )
179  {
180  if( m_compPerfRate[iperf][ic] > LvArray::NumericLimits< real64 >::epsilon )
181  {
182  m_numCrossFlowPerforations += 1;
183  }
184  }
185  for( integer ke = 0; ke < 2; ++ke )
186  {
187  localIndex localDofIndexPres = ke * resNumDOF;
188 
189  localPerfJacobian[TAG::RES * numComp + ic][localDofIndexPres] = m_dt * m_dCompPerfRate[iperf][ke][ic][CP_Deriv::dP];
190  localPerfJacobian[TAG::WELL * numComp + ic][localDofIndexPres] = -m_dt * m_dCompPerfRate[iperf][ke][ic][CP_Deriv::dP];
191  for( integer jc = 0; jc < numComp; ++jc )
192  {
193  localIndex const localDofIndexComp = localDofIndexPres + jc + 1;
194 
195  localPerfJacobian[TAG::RES * numComp + ic][localDofIndexComp] = m_dt * m_dCompPerfRate[iperf][ke][ic][CP_Deriv::dC+jc];
196  localPerfJacobian[TAG::WELL * numComp + ic][localDofIndexComp] = -m_dt * m_dCompPerfRate[iperf][ke][ic][CP_Deriv::dC+jc];
197  }
198  if constexpr ( IS_THERMAL )
199  {
200  localIndex localDofIndexTemp = localDofIndexPres + NC + 1;
201  localPerfJacobian[TAG::RES * numComp + ic][localDofIndexTemp] = m_dt * m_dCompPerfRate[iperf][ke][ic][CP_Deriv::dT];
202  localPerfJacobian[TAG::WELL * numComp + ic][localDofIndexTemp] = -m_dt * m_dCompPerfRate[iperf][ke][ic][CP_Deriv::dT];
203  }
204  }
205  }
206 
207  if( m_useTotalMassEquation )
208  {
209  // Apply equation/variable change transformation(s)
210  stackArray1d< real64, 2 * resNumDOF > work( 2 * resNumDOF );
211  shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numComp, resNumDOF * 2, 2, localPerfJacobian, work );
212  shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( numComp, numComp, 2, localPerf );
213  }
214 
215  for( localIndex i = 0; i < localPerf.size(); ++i )
216  {
217  if( eqnRowIndices[i] >= 0 && eqnRowIndices[i] < m_localMatrix.numRows() )
218  {
219  m_localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( eqnRowIndices[i],
220  dofColIndices.data(),
221  localPerfJacobian[i].dataIfContiguous(),
222  2 * resNumDOF );
223  RAJA::atomicAdd( parallelDeviceAtomic{}, &m_localRhs[eqnRowIndices[i]], localPerf[i] );
224  }
225  }
226  compFluxKernelOp( resOffset, wellElemOffset, dofColIndices, iwelem );
227  }
228  }
229 
230 
239  template< typename POLICY, typename KERNEL_TYPE >
240  static void
241  launch( localIndex const numElements,
242  KERNEL_TYPE const & kernelComponent )
243  {
245  forAll< POLICY >( numElements, [=] GEOS_HOST_DEVICE ( localIndex const ie )
246  {
247  kernelComponent.computeFlux( ie );
248 
249  } );
250  }
251 
252 protected:
253 
255  real64 const m_dt;
256 
259 
260  globalIndex const m_rankOffset;
261  // Perfoation variables
262  arrayView2d< real64 const > const m_compPerfRate;
263  arrayView4d< real64 const > const m_dCompPerfRate;
264  arrayView1d< localIndex const > const m_perfWellElemIndex;
265  arrayView1d< integer const > const m_perfStatus;
266  // Element region, subregion, index
267  arrayView1d< globalIndex const > const m_wellElemDofNumber;
269  arrayView1d< localIndex const > const m_resElementRegion;
270  arrayView1d< localIndex const > const m_resElementSubRegion;
271  arrayView1d< localIndex const > const m_resElementIndex;
272 
273  // RHS and Jacobian
274  arrayView1d< real64 > const m_localRhs;
276 
277  bool const m_detectCrossflow;
278  integer & m_numCrossFlowPerforations;
279  integer const m_useTotalMassEquation;
280 };
281 
286 {
287 public:
288 
302  template< typename POLICY >
303  static void
304  createAndLaunch( integer const numComps,
305  real64 const dt,
306  globalIndex const rankOffset,
307  string const wellDofKey,
308  WellElementSubRegion const & subRegion,
310  PerforationData const * const perforationData,
311  MultiFluidBase const & fluid,
312  BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags,
313  bool const & detectCrossflow,
314  integer & numCrossFlowPerforations,
315  arrayView1d< real64 > const & localRhs,
317  )
318  {
319  isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComps, [&]( auto NC )
320  {
321  integer constexpr NUM_COMP = NC();
322 
324  kernelType kernel( dt, rankOffset, wellDofKey, subRegion, resDofNumber, perforationData,
325  fluid, localRhs, localMatrix, detectCrossflow, numCrossFlowPerforations, kernelFlags );
326  kernelType::template launch< POLICY >( perforationData->size(), kernel );
327  } );
328 
329  }
330 };
331 
337 template< integer NC, integer IS_THERMAL >
339 {
340 public:
341 
343  static constexpr integer numComp = NC;
344  static constexpr integer resNumDOF = NC+1+IS_THERMAL;
345 
346  // Well jacobian column and row indicies
349 
352 
353  using CP_Deriv = multifluid::DerivativeOffsetC< NC, IS_THERMAL >;
354 
356 
357 
358 
360  static constexpr integer numDof = WJ_COFFSET::nDer;
361 
363  static constexpr integer numEqn = WJ_ROFFSET::nEqn - 2;
364 
380  globalIndex const rankOffset,
381  string const wellDofKey,
382  WellElementSubRegion const & subRegion,
383  PerforationData const * const perforationData,
384  MultiFluidBase const & fluid,
385  arrayView1d< real64 > const & localRhs,
386  CRSMatrixView< real64, globalIndex const > const & localMatrix,
387  BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags )
388  :
389  m_dt( dt ),
390  m_numPhases ( fluid.numFluidPhases()),
391  m_rankOffset( rankOffset ),
392  m_compPerfRate( perforationData->getField< fields::well::compPerforationRate >() ),
393  m_dCompPerfRate( perforationData->getField< fields::well::dCompPerforationRate >() ),
394  m_perfWellElemIndex( perforationData->getField< fields::perforation::wellElementIndex >() ),
395  m_wellElemDofNumber( subRegion.getReference< array1d< globalIndex > >( wellDofKey ) ),
396  m_localRhs( localRhs ),
397  m_localMatrix( localMatrix ),
398  m_useTotalMassEquation ( kernelFlags.isSet( isothermalCompositionalMultiphaseBaseKernels::KernelFlags::TotalMassEquation ) )
399  { }
400 
401 
410  template< typename FUNC = NoOpFunc >
412  inline
413  void computeFlux( localIndex const iperf,
414  FUNC && compFluxKernelOp = NoOpFunc{} ) const
415  {
416 
417  using namespace compositionalMultiphaseUtilities;
418  // local working variables and arrays
419  stackArray1d< localIndex, numComp > eqnRowIndices( numComp );
420  stackArray1d< globalIndex, resNumDOF > dofColIndices( resNumDOF );
421 
422  stackArray1d< real64, numComp > localPerf( numComp );
423  stackArray2d< real64, numComp *resNumDOF > localPerfJacobian( numComp, resNumDOF );
424 
425  // get the reservoir (sub)region and element indices
426  //localIndex const er = m_resElementRegion[iperf];
427  //localIndex const esr = m_resElementSubRegion[iperf];
428  //localIndex const ei = m_resElementIndex[iperf];
429 
430  // get the well element index for this perforation
431  localIndex const iwelem = m_perfWellElemIndex[iperf];
432  //globalIndex const resOffset = m_resElemDofNumber[er][esr][ei];
433  globalIndex const wellElemOffset = m_wellElemDofNumber[iwelem];
434 
435  for( integer ic = 0; ic < numComp; ++ic )
436  {
437  eqnRowIndices[ ic] = LvArray::integerConversion< localIndex >( wellElemOffset - m_rankOffset ) + WJ_ROFFSET::MASSBAL + ic;
438  }
439  for( integer jdof = 0; jdof < NC+1; ++jdof )
440  {
441  dofColIndices[ jdof] = wellElemOffset + WJ_COFFSET::dP + jdof;
442  }
443  // For temp its different
444  if constexpr ( IS_THERMAL )
445  {
446  dofColIndices[ NC+1 ] = wellElemOffset + WJ_COFFSET::dT;
447  }
448  // populate local flux vector and derivatives
449 
450  for( integer ic = 0; ic < numComp; ++ic )
451  {
452  localPerf[ic] = -m_dt * m_compPerfRate[iperf][ic];
453  }
454  for( integer ic = 0; ic < numComp; ++ic )
455  {
456  localIndex localDofIndexPres = 0;
457 
458  localPerfJacobian[ic][localDofIndexPres] = -m_dt * m_dCompPerfRate[iperf][TAG::WELL ][ic][CP_Deriv::dP];
459  for( integer jc = 0; jc < numComp; ++jc )
460  {
461  localIndex const localDofIndexComp = localDofIndexPres + jc + 1;
462 
463  localPerfJacobian[ic][localDofIndexComp] = -m_dt * m_dCompPerfRate[iperf][TAG::WELL ][ic][CP_Deriv::dC+jc];
464  }
465  if constexpr ( IS_THERMAL )
466  {
467  localIndex localDofIndexTemp = localDofIndexPres + NC + 1;
468  localPerfJacobian[ic][localDofIndexTemp] = -m_dt * m_dCompPerfRate[iperf][TAG::WELL ][ic][CP_Deriv::dT];
469  }
470  }
471 
472  if( m_useTotalMassEquation )
473  {
474  stackArray1d< real64, resNumDOF > work( resNumDOF );
475  shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numComp, resNumDOF, 1, localPerfJacobian, work );
476 
477  // Apply equation/variable change transformation(s)
478  shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( numComp, numComp, 1, localPerf );
479  }
480 
481  for( localIndex i = 0; i < localPerf.size(); ++i )
482  {
483  if( eqnRowIndices[i] >= 0 && eqnRowIndices[i] < m_localMatrix.numRows() )
484  {
485  m_localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( eqnRowIndices[i],
486  dofColIndices.data(),
487  localPerfJacobian[i].dataIfContiguous(),
488  resNumDOF );
489  RAJA::atomicAdd( parallelDeviceAtomic{}, &m_localRhs[eqnRowIndices[i]], localPerf[i] );
490  }
491  }
492 
493  compFluxKernelOp( wellElemOffset, iwelem, dofColIndices );
494 
495  }
496 
497 
506  template< typename POLICY, typename KERNEL_TYPE >
507  static void
508  launch( localIndex const numElements,
509  KERNEL_TYPE const & kernelComponent )
510  {
512  forAll< POLICY >( numElements, [=] GEOS_HOST_DEVICE ( localIndex const ie )
513  {
514  kernelComponent.computeFlux( ie );
515 
516  } );
517  }
518 
519 protected:
520 
522  real64 const m_dt;
523 
526 
527  globalIndex const m_rankOffset;
528 // Perfoation variables
529  arrayView2d< real64 const > const m_compPerfRate;
530  arrayView4d< real64 const > const m_dCompPerfRate;
531  arrayView1d< localIndex const > const m_perfWellElemIndex;
532 
533 // Element region, subregion, index
534  arrayView1d< globalIndex const > const m_wellElemDofNumber;
535 
536 // RHS and Jacobian
537  arrayView1d< real64 > const m_localRhs;
539 
540  integer const m_useTotalMassEquation;
541 };
542 
547 {
548 public:
549 
563  template< typename POLICY >
564  static void
565  createAndLaunch( integer const numComps,
566  real64 const dt,
567  globalIndex const rankOffset,
568  string const wellDofKey,
569  WellElementSubRegion const & subRegion,
570  PerforationData const * const perforationData,
571  MultiFluidBase const & fluid,
572  arrayView1d< real64 > const & localRhs,
573  CRSMatrixView< real64, globalIndex const > const & localMatrix,
574  BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags )
575  {
576  isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComps, [&]( auto NC )
577  {
578  integer constexpr NUM_COMP = NC();
579 
581  kernelType kernel( dt, rankOffset, wellDofKey, subRegion, perforationData,
582  fluid, localRhs, localMatrix, kernelFlags );
583  kernelType::template launch< POLICY >( perforationData->size(), kernel );
584  } );
585 
586  }
587 };
588 /********************************************************/
594 template< integer NC, integer IS_THERMAL >
596 {
597 public:
600  static constexpr integer numComp = NC;
601  static constexpr integer resNumDOF = NC+1+IS_THERMAL;
602 
603  // Well jacobian column and row indicies
606 
609 
610  using CP_Deriv = multifluid::DerivativeOffsetC< NC, IS_THERMAL >;
611 
613 
614  using Base::m_dt;
615  using Base::m_localRhs;
616  using Base::m_localMatrix;
617  using Base::m_rankOffset;
618 
619 
620 
622  static constexpr integer numDof = WJ_COFFSET::nDer;
623 
625  static constexpr integer numEqn = WJ_ROFFSET::nEqn - 2;
626 
642  bool const thermalEffectsEnabled,
643  integer const isProducer,
644  globalIndex const rankOffset,
645  string const wellDofKey,
646  WellElementSubRegion const & subRegion,
648  PerforationData const * const perforationData,
649  MultiFluidBase const & fluid,
650  arrayView1d< real64 > const & localRhs,
651  CRSMatrixView< real64, globalIndex const > const & localMatrix,
652  bool const & detectCrossflow,
653  integer & numCrossFlowPerforations,
654  BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags )
655  : Base( dt,
656  rankOffset,
657  wellDofKey,
658  subRegion,
659  resDofNumber,
660  perforationData,
661  fluid,
662  localRhs,
663  localMatrix,
664  detectCrossflow,
665  numCrossFlowPerforations,
666  kernelFlags ),
667  m_thermalEffectsEnabled( thermalEffectsEnabled ),
668  m_isProducer( isProducer ),
669  m_globalWellElementIndex( subRegion.getGlobalWellElementIndex() ),
670  m_energyPerfFlux( perforationData->getField< fields::well::energyPerforationFlux >()),
671  m_dEnergyPerfFlux( perforationData->getField< fields::well::dEnergyPerforationFlux >())
672 
673  { }
674 
675 
685  inline
686  void computeFlux( localIndex const iperf ) const
687  {
688  Base::computeFlux( iperf, [&] ( globalIndex const & resOffset,
689  globalIndex const & wellElemOffset,
691  localIndex const iwelem )
692  {
693  if( !m_thermalEffectsEnabled )
694  return;
695  // No energy equation if top element and Injector
696  // Top element defined by global index == 0
697  // Assumption is global index == 0 is top segment with fixed temp BC
698  if( !m_isProducer )
699  {
700  if( m_globalWellElementIndex[iwelem] == 0 )
701  return;
702  }
703  // local working variables and arrays
704  stackArray1d< localIndex, 2* numComp > eqnRowIndices( 2* numComp );
705 
707  stackArray2d< real64, 2 * resNumDOF * 2 * numComp > localPerfJacobian( 2, 2 * resNumDOF );
708 
709 
710  // equantion offsets - note res and well have different equation lineups
711  eqnRowIndices[TAG::RES ] = LvArray::integerConversion< localIndex >( resOffset - m_rankOffset ) + NC + 1;
712  eqnRowIndices[TAG::WELL ] = LvArray::integerConversion< localIndex >( wellElemOffset - m_rankOffset ) + WJ_ROFFSET::ENERGYBAL;
713 
714  // populate local flux vector and derivatives
715  localPerf[TAG::RES ] = m_dt * m_energyPerfFlux[iperf];
716  localPerf[TAG::WELL ] = -m_dt * m_energyPerfFlux[iperf];
717 
718  for( integer ke = 0; ke < 2; ++ke )
719  {
720  localIndex localDofIndexPres = ke * resNumDOF;
721  localPerfJacobian[TAG::RES ][localDofIndexPres] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP];
722  localPerfJacobian[TAG::WELL ][localDofIndexPres] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP];
723 
724  // populate local flux vector and derivatives
725  for( integer ic = 0; ic < numComp; ++ic )
726  {
727  localIndex const localDofIndexComp = localDofIndexPres + ic + 1;
728  localPerfJacobian[TAG::RES ][localDofIndexComp] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dC+ic];
729  localPerfJacobian[TAG::WELL][localDofIndexComp] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dC+ic];
730  }
731  localPerfJacobian[TAG::RES ][localDofIndexPres+NC+1] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT];
732  localPerfJacobian[TAG::WELL][localDofIndexPres+NC+1] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT];
733  }
734 
735 
736  for( localIndex i = 0; i < localPerf.size(); ++i )
737  {
738  if( eqnRowIndices[i] >= 0 && eqnRowIndices[i] < m_localMatrix.numRows() )
739  {
740  m_localMatrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( eqnRowIndices[i],
741  dofColIndices.data(),
742  localPerfJacobian[i].dataIfContiguous(),
743  2 * resNumDOF );
744  RAJA::atomicAdd( parallelDeviceAtomic{}, &m_localRhs[eqnRowIndices[i]], localPerf[i] );
745  }
746  }
747  } );
748  }
749 
750 
759  template< typename POLICY, typename KERNEL_TYPE >
760  static void
761  launch( localIndex const numElements,
762  KERNEL_TYPE const & kernelComponent )
763  {
765  forAll< POLICY >( numElements, [=] GEOS_HOST_DEVICE ( localIndex const ie )
766  {
767  kernelComponent.computeFlux( ie );
768 
769  } );
770  }
771 
772 protected:
775 
778 
781 
784  arrayView3d< real64 const > const m_dEnergyPerfFlux;
785 };
786 
791 {
792 public:
793 
807  template< typename POLICY >
808  static void
809  createAndLaunch( integer const numComps,
810  integer const thermalEffectsEnabled,
811  integer const isProducer,
812  real64 const dt,
813  globalIndex const rankOffset,
814  string const wellDofKey,
815  WellElementSubRegion const & subRegion,
817  PerforationData const * const perforationData,
818  MultiFluidBase const & fluid,
819  BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags,
820  bool const & detectCrossflow,
821  integer & numCrossFlowPerforations,
822  arrayView1d< real64 > const & localRhs,
824  )
825  {
826  isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComps, [&]( auto NC )
827  {
828  integer constexpr NUM_COMP = NC();
829 
831  kernelType kernel( dt, thermalEffectsEnabled, isProducer, rankOffset, wellDofKey, subRegion, resDofNumber, perforationData,
832  fluid, localRhs, localMatrix, detectCrossflow, numCrossFlowPerforations, kernelFlags );
833  kernelType::template launch< POLICY >( perforationData->size(), kernel );
834  } );
835 
836  }
837 };
838 
839 /********************************************************/
845 template< integer NC, integer IS_THERMAL >
847 {
848 public:
851  static constexpr integer numComp = NC;
852  static constexpr integer resNumDOF = NC+1+IS_THERMAL;
853 
854  // Well jacobian column and row indicies
857 
860 
861  using CP_Deriv = multifluid::DerivativeOffsetC< NC, IS_THERMAL >;
862 
864 
865  using Base::m_dt;
866  using Base::m_localRhs;
867  using Base::m_localMatrix;
868  using Base::m_rankOffset;
869 
870 
871 
873  static constexpr integer numDof = WJ_COFFSET::nDer;
874 
876  static constexpr integer numEqn = WJ_ROFFSET::nEqn - 2;
877 
893  bool const thermalEffectsEnabled,
894  integer const isProducer,
895  globalIndex const rankOffset,
896  string const wellDofKey,
897  WellElementSubRegion const & subRegion,
898  PerforationData const * const perforationData,
899  MultiFluidBase const & fluid,
900  arrayView1d< real64 > const & localRhs,
901  CRSMatrixView< real64, globalIndex const > const & localMatrix,
902  BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags )
903  : Base( dt,
904  rankOffset,
905  wellDofKey,
906  subRegion,
907  perforationData,
908  fluid,
909  localRhs,
910  localMatrix,
911  kernelFlags ),
912  m_thermalEffectsEnabled( thermalEffectsEnabled ),
913  m_isProducer( isProducer ),
914  m_globalWellElementIndex( subRegion.getGlobalWellElementIndex() ),
915  m_energyPerfFlux( perforationData->getField< fields::well::energyPerforationFlux >()),
916  m_dEnergyPerfFlux( perforationData->getField< fields::well::dEnergyPerforationFlux >())
917 
918  { }
919 
920 
930  inline
931  void computeFlux( localIndex const iperf ) const
932  {
933  Base::computeFlux( iperf, [&] ( globalIndex const & wellElemOffset,
934  localIndex const iwelem,
936  {
937  GEOS_UNUSED_VAR( dofColIndices );
938  if( !m_thermalEffectsEnabled ) // tjb iso
939  return;
940  // No energy equation if top element and Injector
941  // Top element defined by global index == 0
942  // Assumption is global index == 0 is top segment with fixed temp BC
943  if( !m_isProducer )
944  {
945  if( m_globalWellElementIndex[iwelem] == 0 )
946  return;
947  }
948  // local working variables and arrays
949  localIndex eqnRowIndices = LvArray::integerConversion< localIndex >( wellElemOffset - m_rankOffset ) + WJ_ROFFSET::ENERGYBAL;
950 
951  stackArray2d< real64, resNumDOF > localPerfJacobian( 1, resNumDOF );
952  // populate local flux vector and derivatives
953 
954  real64 localPerf = -m_dt * m_energyPerfFlux[iperf];
955 
956  // std::cout << "Local perf: " << iperf << " " << localPerf << std::endl;
957  localIndex localDofIndexPres = 0;
958  localPerfJacobian [0][localDofIndexPres] = -m_dt * m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dP];
959 
960  // populate local flux vector and derivatives
961  for( integer ic = 0; ic < numComp; ++ic )
962  {
963  localIndex const localDofIndexComp = localDofIndexPres + ic + 1;
964  localPerfJacobian [0][localDofIndexComp] = -m_dt * m_dEnergyPerfFlux[iperf][TAG::WELL ][CP_Deriv::dC+ic];
965  }
966  localPerfJacobian [0][localDofIndexPres+NC+1] = -m_dt * m_dEnergyPerfFlux[iperf][TAG::WELL ][CP_Deriv::dT];
967 
968  if( eqnRowIndices >= 0 && eqnRowIndices < m_localMatrix.numRows() )
969  {
970  // tjb iso
971  m_localMatrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( eqnRowIndices,
972  dofColIndices.data(),
973  localPerfJacobian[0].dataIfContiguous(),
974  resNumDOF );
975  RAJA::atomicAdd( parallelDeviceAtomic{}, &m_localRhs[eqnRowIndices], localPerf );
976  }
977  } );
978 
979 
980  }
981 
982 
991  template< typename POLICY, typename KERNEL_TYPE >
992  static void
993  launch( localIndex const numElements,
994  KERNEL_TYPE const & kernelComponent )
995  {
997  forAll< POLICY >( numElements, [=] GEOS_HOST_DEVICE ( localIndex const ie )
998  {
999  kernelComponent.computeFlux( ie );
1000 
1001  } );
1002  }
1003 
1004 protected:
1009 
1012 
1015  arrayView3d< real64 const > const m_dEnergyPerfFlux;
1016 };
1017 
1022 {
1023 public:
1024 
1038  template< typename POLICY >
1039  static void
1040  createAndLaunch( integer const numComps,
1041  WellControls const & wellControls,
1042  integer const isProducer,
1043  real64 const dt,
1044  globalIndex const rankOffset,
1045  string const wellDofKey,
1046  WellElementSubRegion const & subRegion,
1047  PerforationData const * const perforationData,
1048  MultiFluidBase const & fluid,
1049  BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags,
1050  arrayView1d< real64 > const & localRhs,
1051  CRSMatrixView< real64, globalIndex const > const & localMatrix
1052  )
1053  {
1054  isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComps, [&]( auto NC )
1055  {
1056  integer constexpr NUM_COMP = NC();
1057 
1059  kernelType kernel( dt, wellControls.thermalEffectsEnabled(), isProducer, rankOffset, wellDofKey, subRegion, perforationData,
1060  fluid, localRhs, localMatrix, kernelFlags );
1061  kernelType::template launch< POLICY >( perforationData->size(), kernel );
1062  } );
1063 
1064  }
1065 };
1066 
1067 } // end namespace coupledReservoirAndWellKernels
1068 
1069 } // end namespace geos
1070 
1071 #endif // GEOS_PHYSICSSOLVERS_MULTIPHYSICS_COUPLEDRESERVOIRANDWELLS_HPP
GEOS_HOST_DEVICE void shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum(integer const numRowsToShift, integer const numRowsInBlock, integer const numColsInBlock, integer const numBlocks, MATRIX &&mat, VEC &&work)
In each block, shift the elements from 0 to numRowsToShift-1, shifts all rows one position ahead and ...
GEOS_HOST_DEVICE void shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum(integer const numRowsToShift, integer const numRowsInBlock, integer const numBlocks, VEC &&v)
In each block, shift the elements from 0 to numRowsToShift-1 one position ahead and replaces the firs...
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
typename ElementViewAccessor< VIEWTYPE >::NestedViewTypeConst ElementViewConst
The ElementViewAccessor at the ElementRegionManager level is the type resulting from ElementViewAcces...
This class describes the controls used to operate a well.
This class describes a collection of local well elements and perforations.
static void createAndLaunch(integer const numComps, real64 const dt, globalIndex const rankOffset, string const wellDofKey, WellElementSubRegion const &subRegion, ElementRegionManager::ElementViewConst< arrayView1d< globalIndex const > > const resDofNumber, PerforationData const *const perforationData, MultiFluidBase const &fluid, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags, bool const &detectCrossflow, integer &numCrossFlowPerforations, arrayView1d< real64 > const &localRhs, CRSMatrixView< real64, globalIndex const > const &localMatrix)
Create a new kernel and launch.
GEOS_HOST_DEVICE void computeFlux(localIndex const iperf, FUNC &&compFluxKernelOp=NoOpFunc{}) const
Compute the local flux contributions to the residual and Jacobian.
IsothermalCompositionalMultiPhaseFluxKernel(real64 const dt, globalIndex const rankOffset, string const wellDofKey, WellElementSubRegion const &subRegion, ElementRegionManager::ElementViewConst< arrayView1d< globalIndex const > > const resDofNumber, PerforationData const *const perforationData, MultiFluidBase const &fluid, arrayView1d< real64 > const &localRhs, CRSMatrixView< real64, globalIndex const > const &localMatrix, bool const &detectCrossflow, integer &numCrossFlowPerforations, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags)
Constructor for the kernel interface.
static void launch(localIndex const numElements, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
static void createAndLaunch(integer const numComps, real64 const dt, globalIndex const rankOffset, string const wellDofKey, WellElementSubRegion const &subRegion, PerforationData const *const perforationData, MultiFluidBase const &fluid, arrayView1d< real64 > const &localRhs, CRSMatrixView< real64, globalIndex const > const &localMatrix, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags)
Create a new kernel and launch.
IsothermalCompositionalMultiPhaseWellFluxKernel(real64 const dt, globalIndex const rankOffset, string const wellDofKey, WellElementSubRegion const &subRegion, PerforationData const *const perforationData, MultiFluidBase const &fluid, arrayView1d< real64 > const &localRhs, CRSMatrixView< real64, globalIndex const > const &localMatrix, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags)
Constructor for the kernel interface.
GEOS_HOST_DEVICE void computeFlux(localIndex const iperf, FUNC &&compFluxKernelOp=NoOpFunc{}) const
Compute the local flux contributions to the residual and Jacobian.
static void launch(localIndex const numElements, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
static void createAndLaunch(integer const numComps, integer const thermalEffectsEnabled, integer const isProducer, real64 const dt, globalIndex const rankOffset, string const wellDofKey, WellElementSubRegion const &subRegion, ElementRegionManager::ElementViewConst< arrayView1d< globalIndex const > > const resDofNumber, PerforationData const *const perforationData, MultiFluidBase const &fluid, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags, bool const &detectCrossflow, integer &numCrossFlowPerforations, arrayView1d< real64 > const &localRhs, CRSMatrixView< real64, globalIndex const > const &localMatrix)
Create a new kernel and launch.
GEOS_HOST_DEVICE void computeFlux(localIndex const iperf) const
Compute the local flux contributions to the residual and Jacobian.
bool const m_thermalEffectsEnabled
Flag specifying whether thermal effects are enabled.
ThermalCompositionalMultiPhaseFluxKernel(real64 const dt, bool const thermalEffectsEnabled, integer const isProducer, globalIndex const rankOffset, string const wellDofKey, WellElementSubRegion const &subRegion, ElementRegionManager::ElementViewConst< arrayView1d< globalIndex const > > const resDofNumber, PerforationData const *const perforationData, MultiFluidBase const &fluid, arrayView1d< real64 > const &localRhs, CRSMatrixView< real64, globalIndex const > const &localMatrix, bool const &detectCrossflow, integer &numCrossFlowPerforations, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags)
Constructor for the kernel interface.
static void launch(localIndex const numElements, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
arrayView1d< globalIndex const > m_globalWellElementIndex
Global index of local element.
static void createAndLaunch(integer const numComps, WellControls const &wellControls, integer const isProducer, real64 const dt, globalIndex const rankOffset, string const wellDofKey, WellElementSubRegion const &subRegion, PerforationData const *const perforationData, MultiFluidBase const &fluid, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags, arrayView1d< real64 > const &localRhs, CRSMatrixView< real64, globalIndex const > const &localMatrix)
Create a new kernel and launch.
GEOS_HOST_DEVICE void computeFlux(localIndex const iperf) const
Compute the local flux contributions to the residual and Jacobian.
arrayView1d< globalIndex const > m_globalWellElementIndex
Global index of local element.
ThermalCompositionalMultiPhaseWellFluxKernel(real64 const dt, bool const thermalEffectsEnabled, integer const isProducer, globalIndex const rankOffset, string const wellDofKey, WellElementSubRegion const &subRegion, PerforationData const *const perforationData, MultiFluidBase const &fluid, arrayView1d< real64 > const &localRhs, CRSMatrixView< real64, globalIndex const > const &localMatrix, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags)
Constructor for the kernel interface.
static void launch(localIndex const numElements, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
localIndex size() const
Get the "size" of the group, which determines the number of elements in resizable wrappers.
Definition: Group.hpp:1315
bool thermalEffectsEnabled() const
Are thermal effects enabled.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:179
StackArray< T, 2, MAXSIZE > stackArray2d
Alias for 2D stack array.
Definition: DataTypes.hpp:203
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
Definition: DataTypes.hpp:87
StackArray< T, 1, MAXSIZE > stackArray1d
Alias for 1D stack array.
Definition: DataTypes.hpp:187
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
LvArray::CRSMatrixView< T, COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
Definition: DataTypes.hpp:309
ArrayView< T, 4, USD > arrayView4d
Alias for 4D array view.
Definition: DataTypes.hpp:227
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:195
int integer
Signed integer type.
Definition: DataTypes.hpp:81
Array< T, 1 > array1d
Alias for 1D array.
Definition: DataTypes.hpp:175
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
Definition: DataTypes.hpp:211