GEOS
ThermalCompositionalMultiphaseWellKernels.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_THERMALCOMPOSITIONALMULTIPHASEWELLKERNELS_HPP
21 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_THERMALCOMPOSITIONALMULTIPHASEWELLKERNELS_HPP
22 
25 
26 namespace geos
27 {
28 
29 namespace thermalCompositionalMultiphaseWellKernels
30 {
31 
32 using namespace constitutive;
33 
34 /******************************** TotalMassDensityKernel ****************************/
35 
42 template< integer NUM_COMP, integer NUM_PHASE >
44 {
45 public:
47  using Base::m_dCompFrac_dCompDens;
48  using Base::m_dPhaseMassDens;
49  using Base::m_dPhaseVolFrac;
50  using Base::m_dTotalMassDens;
51  using Base::m_phaseMassDens;
52  using Base::m_phaseVolFrac;
53  using Base::m_totalMassDens;
54  using Base::numComp;
55  using Base::numPhase;
56 
63  MultiFluidBase const & fluid )
64  : Base( subRegion, fluid )
65  {}
66 
73  GEOS_HOST_DEVICE inline void compute( localIndex const ei ) const
74  {
75  using Deriv = multifluid::DerivativeOffset;
76 
77  arraySlice1d< real64 const, compflow::USD_PHASE - 1 > phaseVolFrac = m_phaseVolFrac[ei];
78  arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > dPhaseVolFrac = m_dPhaseVolFrac[ei];
79  arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > phaseMassDens = m_phaseMassDens[ei][0];
80  arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > dPhaseMassDens = m_dPhaseMassDens[ei][0];
81 
82  real64 & dTotalMassDens_dT = m_dTotalMassDens[ei][Deriv::dT];
83  dTotalMassDens_dT=0.0;
84  // Call the base compute the compute the total mass density and derivatives
85  return Base::compute( ei, [&]( localIndex const ip )
86  {
87  dTotalMassDens_dT += dPhaseVolFrac[ip][Deriv::dT] * phaseMassDens[ip] + phaseVolFrac[ip] * dPhaseMassDens[ip][Deriv::dT];
88  } );
89  }
90 
91 protected:
92  // outputs
93  arrayView1d< real64 > m_dTotalMassDens_dTemp;
94 };
95 
100 {
101 public:
110  template< typename POLICY >
111  static void
112  createAndLaunch( integer const numComp,
113  integer const numPhase,
114  ObjectManagerBase & subRegion,
115  MultiFluidBase const & fluid )
116  {
117  if( numPhase == 2 )
118  {
119  isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&]( auto NC )
120  {
121  integer constexpr NUM_COMP = NC();
122  TotalMassDensityKernel< NUM_COMP, 2 > kernel( subRegion, fluid );
123  TotalMassDensityKernel< NUM_COMP, 2 >::template launch< POLICY >( subRegion.size(), kernel );
124  } );
125  }
126  else if( numPhase == 3 )
127  {
128  isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&]( auto NC )
129  {
130  integer constexpr NUM_COMP = NC();
131  TotalMassDensityKernel< NUM_COMP, 3 > kernel( subRegion, fluid );
132  TotalMassDensityKernel< NUM_COMP, 3 >::template launch< POLICY >( subRegion.size(), kernel );
133  } );
134  }
135  }
136 };
137 
138 /******************************** ResidualNormKernel ********************************/
139 
143 template< localIndex NUM_COMP >
145 {
146 public:
147 
149  static constexpr integer numComp = NUM_COMP;
150 
151 
153 
155  using Base::m_minNormalizer;
156  using Base::m_rankOffset;
157  using Base::m_localResidual;
158  using Base::m_dofNumber;
159 
160  ResidualNormKernel( globalIndex const rankOffset,
161  arrayView1d< real64 const > const & localResidual,
162  arrayView1d< globalIndex const > const & dofNumber,
163  arrayView1d< localIndex const > const & ghostRank,
164  integer const targetPhaseIndex,
165  WellElementSubRegion const & subRegion,
166  MultiFluidBase const & fluid,
167  WellControls const & wellControls,
168  real64 const time,
169  real64 const dt,
170  real64 const minNormalizer )
171  : Base( rankOffset,
172  localResidual,
173  dofNumber,
174  ghostRank,
175  minNormalizer ),
176  m_numPhases( fluid.numFluidPhases()),
177  m_targetPhaseIndex( targetPhaseIndex ),
178  m_dt( dt ),
179  m_isLocallyOwned( subRegion.isLocallyOwned() ),
180  m_iwelemControl( subRegion.getTopWellElementIndex() ),
181  m_isProducer( wellControls.isProducer() ),
182  m_currentControl( wellControls.getControl() ),
183  m_targetBHP( wellControls.getTargetBHP( time ) ),
184  m_targetTotalRate( wellControls.getTargetTotalRate( time ) ),
185  m_targetPhaseRate( wellControls.getTargetPhaseRate( time ) ),
186  m_targetMassRate( wellControls.getTargetMassRate( time ) ),
187  m_volume( subRegion.getElementVolume() ),
188  m_phaseDens_n( fluid.phaseDensity_n() ),
189  m_totalDens_n( fluid.totalDensity_n() ),
190  m_phaseVolFraction_n( subRegion.getField< fields::well::phaseVolumeFraction_n >()),
191  m_phaseInternalEnergy_n( fluid.phaseInternalEnergy_n() )
192  {}
193 
194 
196  void computeMassEnergyNormalizers( localIndex const iwelem,
197  real64 & massNormalizer,
198  real64 & energyNormalizer ) const
199  {
200  massNormalizer = LvArray::math::max( m_minNormalizer, m_totalDens_n[iwelem][0] * m_volume[iwelem] );
201 
202  for( integer ip = 0; ip < m_numPhases; ++ip )
203  {
204  energyNormalizer += m_phaseInternalEnergy_n[iwelem][0][ip] * m_phaseDens_n[iwelem][0][ip] * m_phaseVolFraction_n[iwelem][ip] * m_volume[iwelem];
205  }
206  // warning: internal energy can be negative
207  energyNormalizer = LvArray::math::max( m_minNormalizer, LvArray::math::abs( energyNormalizer ) );
208  }
209 
211  virtual void computeLinf( localIndex const iwelem,
212  LinfStackVariables & stack ) const override
213  {
214  real64 normalizer = 0.0;
215  for( integer idof = 0; idof < WJ_ROFFSET::nEqn; ++idof )
216  {
217 
218  // Step 1: compute a normalizer for the control or pressure equation
219 
220  // for the control equation, we distinguish two cases
221  if( idof == WJ_ROFFSET::CONTROL )
222  {
223 
224  // for the top well element, normalize using the current control
225  if( m_isLocallyOwned && iwelem == m_iwelemControl )
226  {
227  if( m_currentControl == WellControls::Control::BHP )
228  {
229  // the residual entry is in pressure units
230  normalizer = m_targetBHP;
231  }
232  else if( m_currentControl == WellControls::Control::TOTALVOLRATE )
233  {
234  // the residual entry is in volume / time units
235  normalizer = LvArray::math::max( LvArray::math::abs( m_targetTotalRate ), m_minNormalizer );
236  }
237  else if( m_currentControl == WellControls::Control::PHASEVOLRATE )
238  {
239  // the residual entry is in volume / time units
240  normalizer = LvArray::math::max( LvArray::math::abs( m_targetPhaseRate ), m_minNormalizer );
241  }
242  else if( m_currentControl == WellControls::Control::MASSRATE )
243  {
244  // the residual entry is in volume / time units
245  normalizer = LvArray::math::max( LvArray::math::abs( m_targetMassRate ), m_minNormalizer );
246  }
247  }
248  // for the pressure difference equation, always normalize by the BHP
249  else
250  {
251  normalizer = m_targetBHP;
252  }
253  }
254  // Step 2: compute a normalizer for the mass balance equations
255  else if( idof >= WJ_ROFFSET::MASSBAL && idof < WJ_ROFFSET::MASSBAL + numComp )
256  {
257  if( m_isProducer ) // only PHASEVOLRATE is supported for now
258  {
259  // the residual is in mass units
260  normalizer = m_dt * LvArray::math::abs( m_targetPhaseRate ) * m_phaseDens_n[iwelem][0][m_targetPhaseIndex];
261  }
262  else // Type::INJECTOR, only TOTALVOLRATE is supported for now
263  {
264  if( m_currentControl == WellControls::Control::MASSRATE )
265  {
266  normalizer = m_dt * LvArray::math::abs( m_targetMassRate );
267  }
268  else
269  {
270  // the residual is in mass units
271  normalizer = m_dt * LvArray::math::abs( m_targetTotalRate ) * m_totalDens_n[iwelem][0];
272  }
273 
274  }
275 
276  // to make sure that everything still works well if the rate is zero, we add this check
277  normalizer = LvArray::math::max( normalizer, m_volume[iwelem] * m_totalDens_n[iwelem][0] );
278  }
279  // Step 3: compute a normalizer for the volume balance equations
280  else if( idof == WJ_ROFFSET::VOLBAL )
281  {
282  if( m_isProducer ) // only PHASEVOLRATE is supported for now
283  {
284  // the residual is in volume units
285  normalizer = m_dt * LvArray::math::abs( m_targetPhaseRate );
286  }
287  else // Type::INJECTOR, only TOTALVOLRATE is supported for now
288  {
289  if( m_currentControl == WellControls::Control::MASSRATE )
290  {
291  normalizer = m_dt * LvArray::math::abs( m_targetMassRate/ m_totalDens_n[iwelem][0] );
292  }
293  else
294  {
295  normalizer = m_dt * LvArray::math::abs( m_targetTotalRate );
296  }
297 
298  }
299  // to make sure that everything still works well if the rate is zero, we add this check
300  normalizer = LvArray::math::max( normalizer, m_volume[iwelem] );
301  }
302  // step 3: energy residual
303  if( idof == WJ_ROFFSET::ENERGYBAL )
304  {
305  real64 massNormalizer = 0.0, energyNormalizer = 0.0;
306  computeMassEnergyNormalizers( iwelem, massNormalizer, energyNormalizer );
307  real64 const valEnergy = LvArray::math::abs( m_localResidual[stack.localRow + WJ_ROFFSET::ENERGYBAL] ) / energyNormalizer;
308  if( valEnergy > stack.localValue[1] )
309  {
310  stack.localValue[1] = valEnergy;
311  }
312 
313  }
314  else
315  {
316  normalizer = LvArray::math::max( m_minNormalizer, normalizer );
317  // Step 4: compute the contribution to the residual
318  real64 const val = LvArray::math::abs( m_localResidual[stack.localRow + idof] ) / normalizer;
319  if( val > stack.localValue[0] )
320  {
321  stack.localValue[0] = val;
322  }
323  }
324  }
325  }
326 
328  virtual void computeL2( localIndex const iwelem,
329  L2StackVariables & stack ) const override
330  {
331  GEOS_UNUSED_VAR( iwelem, stack );
332  GEOS_ERROR( "The L2 norm is not implemented for CompositionalMultiphaseWell" );
333  }
334 
335 
336 protected:
337 
340 
343 
345  real64 const m_dt;
346 
348  bool const m_isLocallyOwned;
349 
352 
354  bool const m_isProducer;
355 
358  real64 const m_targetBHP;
359  real64 const m_targetTotalRate;
360  real64 const m_targetPhaseRate;
361  real64 const m_targetMassRate;
362 
365 
369  arrayView2d< real64 const, compflow::USD_PHASE > const m_phaseVolFraction_n;
370  arrayView3d< real64 const, multifluid::USD_PHASE > const m_phaseInternalEnergy_n;
371 
372 };
373 
374 /*
375  *@class ResidualNormKernelFactory
376  */
378 {
379 public:
380 
397  template< typename POLICY >
398  static void
399  createAndLaunch( integer const numComp,
400  integer const targetPhaseIndex,
401  globalIndex const rankOffset,
402  string const & dofKey,
403  arrayView1d< real64 const > const & localResidual,
404  WellElementSubRegion const & subRegion,
405  MultiFluidBase const & fluid,
406  WellControls const & wellControls,
407  real64 const time,
408  real64 const dt,
409  real64 const minNormalizer,
410  real64 (& residualNorm)[2] )
411  {
412  isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch ( numComp, [&]( auto NC )
413  {
414 
415  integer constexpr NUM_COMP = NC();
416  using kernelType = ResidualNormKernel< NUM_COMP >;
417  arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey );
418  arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
419 
420  kernelType kernel( rankOffset, localResidual, dofNumber, ghostRank,
421  targetPhaseIndex, subRegion, fluid, wellControls, time, dt, minNormalizer );
422  kernelType::template launchLinf< POLICY >( subRegion.size(), kernel, residualNorm );
423  } );
424  }
425 
426 };
427 
428 /******************************** ElementBasedAssemblyKernel ********************************/
429 
436 template< localIndex NUM_COMP >
438 {
439 public:
441  using Base::m_dCompFrac_dCompDens;
442  using Base::m_dofNumber;
443  using Base::m_dPhaseCompFrac;
444  using Base::m_dPhaseDens;
445  using Base::m_dPhaseVolFrac;
446  using Base::m_dPoro_dPres;
447  using Base::m_elemGhostRank;
448  using Base::m_localMatrix;
449  using Base::m_localRhs;
450  using Base::m_numPhases;
451  using Base::m_phaseCompFrac;
452  using Base::m_phaseCompFrac_n;
453  using Base::m_phaseDens;
454  using Base::m_phaseDens_n;
455  using Base::m_phaseVolFrac;
456  using Base::m_phaseVolFrac_n;
457  using Base::m_porosity;
458  using Base::m_porosity_n;
459  using Base::m_rankOffset;
460  using Base::m_volume;
461  using Base::numComp;
462  using Base::numDof;
463  using Base::numEqn;
464 
465  using FLUID_PROP_COFFSET = multifluid::DerivativeOffsetC< NUM_COMP, 1 >;
466 
467 
479  integer const isProducer,
480  globalIndex const rankOffset,
481  string const dofKey,
482  WellElementSubRegion const & subRegion,
483  MultiFluidBase const & fluid,
484  CRSMatrixView< real64, globalIndex const > const & localMatrix,
485  arrayView1d< real64 > const & localRhs,
486  BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > const kernelFlags )
487  : Base( numPhases, isProducer, rankOffset, dofKey, subRegion, fluid, localMatrix, localRhs, kernelFlags ),
488  m_phaseInternalEnergy_n( fluid.phaseInternalEnergy_n()),
489  m_phaseInternalEnergy( fluid.phaseInternalEnergy()),
490  m_dPhaseInternalEnergy( fluid.dPhaseInternalEnergy())
491  {}
492 
493  struct StackVariables : public Base::StackVariables
494  {
495 public:
498  : Base::StackVariables()
499  {}
500  using Base::StackVariables::eqnRowIndices;
501  using Base::StackVariables::dofColIndices;
502  using Base::StackVariables::localJacobian;
503  using Base::StackVariables::localResidual;
504  using Base::StackVariables::localRow;
505  using Base::StackVariables::volume;
506 
507 
508 
509  };
516  void setup( localIndex const ei,
517  StackVariables & stack ) const
518  {
519  Base::setup( ei, stack );
520 
521 
522 
523  }
524 
533  StackVariables & stack ) const
534  {
535  using Deriv = multifluid::DerivativeOffset;
536 
537  Base::computeAccumulation( ei, stack, [&]( integer const ip
538  , real64 const & phaseAmount
539  , real64 const & phaseAmount_n
540  , real64 const (&dPhaseAmount)[FLUID_PROP_COFFSET::nDer] )
541  {
542  // We are in the loop over phases, ip provides the current phase index.
543  // We have to do two things:
544  // 1- Assemble the derivatives of the component mass balance equations with respect to temperature
545  // 2- Assemble the phase-dependent part of the accumulation term of the energy equation
546 
547  real64 dPhaseInternalEnergy_dC[numComp]{};
548 
549  // construct the slices
550  arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > dCompFrac_dCompDens = m_dCompFrac_dCompDens[ei];
551  arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > phaseInternalEnergy_n = m_phaseInternalEnergy_n[ei][0];
552  arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > phaseInternalEnergy = m_phaseInternalEnergy[ei][0];
553  arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > dPhaseInternalEnergy = m_dPhaseInternalEnergy[ei][0];
554 
555  // Step 1: assemble the phase-dependent part of the accumulation term of the energy equation
556 
557  real64 const phaseEnergy = phaseAmount * phaseInternalEnergy[ip];
558  real64 const phaseEnergy_n = phaseAmount_n * phaseInternalEnergy_n[ip];
559  real64 const dPhaseEnergy_dP = dPhaseAmount[FLUID_PROP_COFFSET::dP] * phaseInternalEnergy[ip]
560  + phaseAmount * dPhaseInternalEnergy[ip][Deriv::dP];
561  real64 const dPhaseEnergy_dT = dPhaseAmount[FLUID_PROP_COFFSET::dT] * phaseInternalEnergy[ip]
562  + phaseAmount * dPhaseInternalEnergy[ip][Deriv::dT];
563  // local accumulation
564  stack.localResidual[numEqn-1] += phaseEnergy - phaseEnergy_n;
565 
566  // derivatives w.r.t. pressure and temperature
567  stack.localJacobian[numEqn-1][0] += dPhaseEnergy_dP;
568  stack.localJacobian[numEqn-1][numDof-1] += dPhaseEnergy_dT;
569 
570  // derivatives w.r.t. component densities
571  applyChainRule( numComp, dCompFrac_dCompDens, dPhaseInternalEnergy[ip], dPhaseInternalEnergy_dC, Deriv::dC );
572  for( integer jc = 0; jc < numComp; ++jc )
573  {
574  stack.localJacobian[numEqn-1][jc + 1] += phaseInternalEnergy[ip] * dPhaseAmount[FLUID_PROP_COFFSET::dC+jc]
575  + dPhaseInternalEnergy_dC[jc] * phaseAmount;
576  }
577  } );
578 
579 
580  }
581 
590  StackVariables & stack ) const
591  {
592  Base::computeVolumeBalance( ei, stack );
593 
594  }
595 
597  void complete( localIndex const ei,
598  StackVariables & stack ) const
599  {
600  // Assemble the component mass balance equations and volume balance equations
601  // Energy balance equation updates to solver matrices included in Base class
602  Base::complete( ei, stack );
603 
604  }
605 
606 protected:
609 
614 
615 };
616 
621 {
622 public:
635  template< typename POLICY >
636  static void
637  createAndLaunch( localIndex const numComps,
638  localIndex const numPhases,
639  integer const isProducer,
640  globalIndex const rankOffset,
641  BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags,
642  string const dofKey,
643  WellElementSubRegion const & subRegion,
644  MultiFluidBase const & fluid,
645  CRSMatrixView< real64, globalIndex const > const & localMatrix,
646  arrayView1d< real64 > const & localRhs )
647  {
648  isothermalCompositionalMultiphaseBaseKernels::
649  internal::kernelLaunchSelectorCompSwitch( numComps, [&]( auto NC )
650  {
651  localIndex constexpr NUM_COMP = NC();
652 
654  kernel( numPhases, isProducer, rankOffset, dofKey, subRegion, fluid, localMatrix, localRhs, kernelFlags );
656  launch< POLICY, ElementBasedAssemblyKernel< NUM_COMP > >( subRegion.size(), kernel );
657  } );
658  }
659 };
660 
661 /******************************** FaceBasedAssemblyKernel ********************************/
662 
668 template< integer NC >
670 {
671 public:
672  static constexpr integer IS_THERMAL = 1;
674 
675  // Well jacobian column and row indicies
678 
679  using CP_Deriv = multifluid::DerivativeOffsetC< NC, IS_THERMAL >;
680 
682 
683 
684  using Base::m_isProducer;
685  using Base::m_dt;
686  using Base::m_localRhs;
687  using Base::m_localMatrix;
688  using Base::m_rankOffset;
689  using Base::maxNumElems;
690  using Base::maxStencilSize;
691  using Base::m_useTotalMassEquation;
692 
694  static constexpr integer numComp = NC;
695 
697  static constexpr integer numDof = WJ_COFFSET::nDer;
698 
700  static constexpr integer numEqn = WJ_ROFFSET::nEqn - 2;
701 
717  globalIndex const rankOffset,
718  string const wellDofKey,
719  WellControls const & wellControls,
720  WellElementSubRegion const & subRegion,
721  MultiFluidBase const & fluid,
722  CRSMatrixView< real64, globalIndex const > const & localMatrix,
723  arrayView1d< real64 > const & localRhs,
724  BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags )
725  : Base( dt
726  , rankOffset
727  , wellDofKey
728  , wellControls
729  , subRegion
730  , localMatrix
731  , localRhs
732  , kernelFlags ),
733  m_numPhases ( fluid.numFluidPhases()),
734  m_globalWellElementIndex( subRegion.getGlobalWellElementIndex() ),
735  m_phaseFraction( fluid.phaseFraction()),
736  m_dPhaseFraction( fluid.dPhaseFraction()),
737  m_phaseEnthalpy( fluid.phaseEnthalpy()),
738  m_dPhaseEnthalpy( fluid.dPhaseEnthalpy())
739  { }
740 
742  {
743 public:
744 
746  StackVariables( localIndex const size )
747  : Base::StackVariables( size )
748  {}
749 
756  };
757 
759  inline
760  void setup( localIndex const iwelem, StackVariables & stack ) const
761  {
762  Base::setup ( iwelem, stack );
763  stack.localEnergyFlux.resize( stack.numConnectedElems );
764  stack.localEnergyFluxJacobian.resize( stack.numConnectedElems, stack.stencilSize * numDof );
765  stack.localEnergyFluxJacobian_dQ.resize( stack.numConnectedElems, 1 );
766  for( integer i=0; i<stack.numConnectedElems; i++ )
767  {
768  stack.localEnergyFlux[i]=0.0;
769  stack.localEnergyFluxJacobian_dQ[i][0]=0.0;
770  for( integer j=0; j< stack.stencilSize * numDof; j++ )
771  stack.localEnergyFluxJacobian[i][j] = 0.0;
772  }
773  }
774 
776  inline
777  void complete( localIndex const iwelem, StackVariables & stack ) const
778  {
779  Base::complete ( iwelem, stack );
780 
781  using namespace compositionalMultiphaseUtilities;
782  if( stack.numConnectedElems ==1 )
783  {
784  // Setup Jacobian global row indicies for energy equation
785  globalIndex oneSidedEqnRowIndices = stack.offsetUp + WJ_ROFFSET::ENERGYBAL - m_rankOffset;
786 
787  if( oneSidedEqnRowIndices >= 0 && oneSidedEqnRowIndices < m_localMatrix.numRows() )
788  {
789 
790  if( !m_isProducer && m_globalWellElementIndex[iwelem] == 0 )
791  {
792  // For top segment energy balance eqn replaced with T(n+1) - T = 0
793  // No other energy balance derivatives
794  // Assumption is global index == 0 is top segment with fixed temp BC
795  for( integer i=0; i< CP_Deriv::nDer; i++ )
796  {
797  stack.localEnergyFluxJacobian[0][i] = 0.0;
798  }
799  stack.localEnergyFluxJacobian_dQ[0][0]=0;
800  stack.localEnergyFlux[0]=0;
801  }
802 
803 
804  // Setup Jacobian global col indicies ( Mapping from local jac order to well jac order)
805  globalIndex oneSidedDofColIndices_dRate = stack.offsetCurrent + WJ_COFFSET::dQ;
806  globalIndex oneSidedDofColIndices_dPresCompTempUp[CP_Deriv::nDer]{};
807 
808  int ioff=0;
809  oneSidedDofColIndices_dPresCompTempUp[ioff++] = stack.offsetUp + WJ_COFFSET::dP;
810  oneSidedDofColIndices_dPresCompTempUp[ioff++] = stack.offsetUp + WJ_COFFSET::dT;
811  for( integer jdof = 0; jdof < NC; ++jdof )
812  {
813  oneSidedDofColIndices_dPresCompTempUp[ioff++] = stack.offsetUp + WJ_COFFSET::dC+ jdof;
814  }
815 
816  m_localMatrix.template addToRow< parallelDeviceAtomic >( oneSidedEqnRowIndices,
817  &oneSidedDofColIndices_dRate,
818  stack.localEnergyFluxJacobian_dQ[0],
819  1 );
820  m_localMatrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( oneSidedEqnRowIndices,
821  oneSidedDofColIndices_dPresCompTempUp,
822  stack.localEnergyFluxJacobian[0],
823  CP_Deriv::nDer );
824  RAJA::atomicAdd( parallelDeviceAtomic{}, &m_localRhs[oneSidedEqnRowIndices], stack.localEnergyFlux[0] );
825  }
826  }
827  else // if ( stack.numConnectedElems == 2 )
828  {
829  globalIndex row_current = stack.offsetCurrent + WJ_ROFFSET::ENERGYBAL - m_rankOffset;
830  globalIndex row_next = stack.offsetNext + WJ_ROFFSET::ENERGYBAL - m_rankOffset;
831 
832  if( !m_isProducer )
833  {
834  if( row_next >= 0 && row_next < m_localMatrix.numRows() )
835  {
836  if( m_globalWellElementIndex[stack.iwelemNext] == 0 )
837  {
838  for( integer i=0; i<CP_Deriv::nDer; i++ )
839  stack.localEnergyFluxJacobian[TAG::NEXT][i] = 0.0;
840  stack.localEnergyFluxJacobian_dQ[TAG::NEXT][0] =0;
841  stack.localEnergyFlux[TAG::NEXT] =0;
842  }
843 
844  }
845  }
846  // Setup Jacobian global row indicies
847  // equations for COMPONENT + ENERGY balances
848  globalIndex eqnRowIndices[2]{};
849 
850  // energy balance equations
851  eqnRowIndices[TAG::CURRENT ] = row_current;
852  eqnRowIndices[TAG::NEXT ] = row_next;
853 
854 
855  // Setup Jacobian global col indicies ( Mapping from local jac order to well jac order)
856  globalIndex dofColIndices[CP_Deriv::nDer]{};
857  globalIndex dofColIndices_dRate = stack.offsetCurrent + WJ_COFFSET::dQ;
858 
859  int ioff=0;
860  // Indice storage order reflects local jac col storage order CP::Deriv order P T DENS
861  // well jacobian order is P DENS Q T
862  dofColIndices[ioff++] = stack.offsetUp + WJ_COFFSET::dP;
863 
864  if constexpr ( IS_THERMAL )
865  {
866  dofColIndices[ioff++] = stack.offsetUp + WJ_COFFSET::dT;
867  }
868  for( integer jdof = 0; jdof < NC; ++jdof )
869  {
870  dofColIndices[ioff++] = stack.offsetUp + WJ_COFFSET::dC+ jdof;
871  }
872  // Note this updates diag and offdiag
873  for( integer i = 0; i < 2; ++i )
874  {
875  if( eqnRowIndices[i] >= 0 && eqnRowIndices[i] < m_localMatrix.numRows() )
876  {
877  m_localMatrix.template addToRow< parallelDeviceAtomic >( eqnRowIndices[i],
878  &dofColIndices_dRate,
879  stack.localEnergyFluxJacobian_dQ[i],
880  1 );
881  m_localMatrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( eqnRowIndices[i],
882  dofColIndices,
883  stack.localEnergyFluxJacobian[i],
884  CP_Deriv::nDer );
885  RAJA::atomicAdd( parallelDeviceAtomic{}, &m_localRhs[eqnRowIndices[i]], stack.localEnergyFlux[i] );
886  }
887  }
888  }
889 
890  }
891 
901  inline
902  void computeFlux( localIndex const iwelem, StackVariables & stack ) const
903  {
904  Base::computeFlux ( iwelem, stack, [&] ( localIndex const & iwelemNext
905  , localIndex const & iwelemUp
906  , real64 const & currentConnRate
907  , real64 const (&dCompFrac_dCompDens)[NC][NC] )
908  {
909 
910  if( iwelemNext < 0 && !m_isProducer ) // exit connection, injector
911  {
912  real64 eflux=0;
913  real64 eflux_dq=0;
914  for( integer ip = 0; ip < m_numPhases; ++ip )
915  {
916  eflux += m_phaseEnthalpy[iwelemUp][0][ip]* m_phaseFraction[iwelemUp][0][ip];
917  eflux_dq += m_phaseEnthalpy[iwelemUp][0][ip] * m_phaseFraction[iwelemUp][0][ip];
918 
919  stack.localEnergyFluxJacobian[0] [CP_Deriv::dP] += m_phaseEnthalpy[iwelemUp][0][ip]*m_dPhaseFraction[iwelemUp][0][ip][CP_Deriv::dP]
920  + m_dPhaseEnthalpy[iwelemUp][0][ip][CP_Deriv::dP]*m_phaseFraction[iwelemUp][0][ip];
921  stack.localEnergyFluxJacobian[0] [CP_Deriv::dT] += m_phaseEnthalpy[iwelemUp][0][ip]*m_dPhaseFraction[iwelemUp][0][ip][CP_Deriv::dT]
922  + m_dPhaseEnthalpy[iwelemUp][0][ip][CP_Deriv::dT]*m_phaseFraction[iwelemUp][0][ip];
923 
924  real64 dProp1_dC[numComp]{};
925  applyChainRule( numComp, dCompFrac_dCompDens, m_dPhaseEnthalpy[iwelemUp][0][ip], dProp1_dC, CP_Deriv::dC );
926  real64 dProp2_dC[numComp]{};
927  applyChainRule( numComp, dCompFrac_dCompDens, m_dPhaseFraction[iwelemUp][0][ip], dProp2_dC, CP_Deriv::dC );
928  for( integer dof=0; dof < numComp; dof++ )
929  {
930  stack.localEnergyFluxJacobian[0] [CP_Deriv::dC+dof] += m_phaseEnthalpy[iwelemUp][0][ip]*dProp2_dC[dof]
931  + dProp1_dC[dof]*m_phaseFraction[iwelemUp][0][ip];
932  }
933  }
934  for( integer dof=0; dof < CP_Deriv::nDer; dof++ )
935  {
936  stack.localEnergyFluxJacobian[0] [dof] *= -m_dt*currentConnRate;
937  }
938  // Energy equation
939  stack.localEnergyFlux[0] = -m_dt * eflux * currentConnRate;
940  stack.localEnergyFluxJacobian_dQ[0][0] = -m_dt * eflux_dq;
941  }
942  else if( ( iwelemNext < 0 && m_isProducer ) || currentConnRate < 0 ) // exit connection, producer
943  {
944  real64 eflux=0;
945  real64 eflux_dq=0;
946  for( integer ip = 0; ip < m_numPhases; ++ip )
947  {
948  eflux += m_phaseEnthalpy[iwelemUp][0][ip]* m_phaseFraction[iwelemUp][0][ip];
949  eflux_dq += m_phaseEnthalpy[iwelemUp][0][ip] * m_phaseFraction[iwelemUp][0][ip];
950  stack.localEnergyFluxJacobian[0] [CP_Deriv::dP] += m_phaseEnthalpy[iwelemUp][0][ip]*m_dPhaseFraction[iwelemUp][0][ip][CP_Deriv::dP]
951  + m_dPhaseEnthalpy[iwelemUp][0][ip][CP_Deriv::dP]*m_phaseFraction[iwelemUp][0][ip];
952  stack.localEnergyFluxJacobian[0] [CP_Deriv::dT] += m_phaseEnthalpy[iwelemUp][0][ip]*m_dPhaseFraction[iwelemUp][0][ip][CP_Deriv::dT]
953  + m_dPhaseEnthalpy[iwelemUp][0][ip][CP_Deriv::dT]*m_phaseFraction[iwelemUp][0][ip];
954 
955  real64 dProp1_dC[numComp]{};
956  applyChainRule( numComp, dCompFrac_dCompDens, m_dPhaseEnthalpy[iwelemUp][0][ip], dProp1_dC, CP_Deriv::dC );
957  real64 dProp2_dC[numComp]{};
958  applyChainRule( numComp, dCompFrac_dCompDens, m_dPhaseFraction[iwelemUp][0][ip], dProp2_dC, CP_Deriv::dC );
959  for( integer dof=0; dof < numComp; dof++ )
960  {
961  stack.localEnergyFluxJacobian[0] [CP_Deriv::dC+dof] += m_phaseEnthalpy[iwelemUp][0][ip]*dProp2_dC[dof]
962  + dProp1_dC[dof]*m_phaseFraction[iwelemUp][0][ip];
963  }
964 
965  }
966 
967  for( integer dof=0; dof < CP_Deriv::nDer; dof++ )
968  {
969  stack.localEnergyFluxJacobian[0][dof] *= -m_dt*currentConnRate;
970  }
971  stack.localEnergyFlux[0] = -m_dt * eflux * currentConnRate;
972  stack.localEnergyFluxJacobian_dQ[0][0] = -m_dt*eflux_dq;
973  }
974  else
975  {
976  real64 eflux=0;
977  real64 eflux_dq=0;
978  for( integer ip = 0; ip < m_numPhases; ++ip )
979  {
980  eflux += m_phaseEnthalpy[iwelemUp][0][ip]* m_phaseFraction[iwelemUp][0][ip];
981  eflux_dq += m_phaseEnthalpy[iwelemUp][0][ip] * m_phaseFraction[iwelemUp][0][ip];
982 
983  real64 dprop_dp = m_phaseEnthalpy[iwelemUp][0][ip]*m_dPhaseFraction[iwelemUp][0][ip][CP_Deriv::dP]
984  + m_dPhaseEnthalpy[iwelemUp][0][ip][CP_Deriv::dP]*m_phaseFraction[iwelemUp][0][ip];
985  real64 dprop_dt = m_phaseEnthalpy[iwelemUp][0][ip]*m_dPhaseFraction[iwelemUp][0][ip][CP_Deriv::dT]
986  + m_dPhaseEnthalpy[iwelemUp][0][ip][CP_Deriv::dT]*m_phaseFraction[iwelemUp][0][ip];
987 
988  stack.localEnergyFluxJacobian[TAG::NEXT ] [CP_Deriv::dP] += dprop_dp;
989  stack.localEnergyFluxJacobian[TAG::NEXT] [CP_Deriv::dT] += dprop_dt;
990 
991  stack.localEnergyFluxJacobian[TAG::CURRENT ] [CP_Deriv::dP] += dprop_dp;
992  stack.localEnergyFluxJacobian[TAG::CURRENT] [CP_Deriv::dT] += dprop_dt;
993 
994  real64 dPE_dC[numComp]{};
995  applyChainRule( numComp, dCompFrac_dCompDens, m_dPhaseEnthalpy[iwelemUp][0][ip], dPE_dC, CP_Deriv::dC );
996  real64 dPF_dC[numComp]{};
997  applyChainRule( numComp, dCompFrac_dCompDens, m_dPhaseFraction[iwelemUp][0][ip], dPF_dC, CP_Deriv::dC );
998 
999  for( integer dof=0; dof < numComp; dof++ )
1000  {
1001  stack.localEnergyFluxJacobian[TAG::NEXT ][CP_Deriv::dC+dof] += m_phaseEnthalpy[iwelemUp][0][ip]*dPF_dC[dof]
1002  +dPE_dC[dof]*m_phaseFraction[iwelemUp][0][ip];
1003  stack.localEnergyFluxJacobian[TAG::CURRENT ][CP_Deriv::dC+dof] += m_phaseEnthalpy[iwelemUp][0][ip]*dPF_dC[dof]
1004  +dPE_dC[dof]*m_phaseFraction[iwelemUp][0][ip];
1005  }
1006  }
1007  stack.localEnergyFlux[TAG::NEXT ] = m_dt * eflux * currentConnRate;
1008  stack.localEnergyFlux[TAG::CURRENT ] = -m_dt * eflux * currentConnRate;
1009  stack.localEnergyFluxJacobian_dQ [TAG::NEXT ][0] = m_dt * eflux_dq;
1010  stack.localEnergyFluxJacobian_dQ [TAG::CURRENT][0] = -m_dt * eflux_dq;
1011  for( integer dof=0; dof < CP_Deriv::nDer; dof++ )
1012  {
1013  stack.localEnergyFluxJacobian[TAG::NEXT ][dof] *= m_dt*currentConnRate;
1014  stack.localEnergyFluxJacobian[TAG::CURRENT ][dof] *= -m_dt*currentConnRate;
1015  }
1016  }
1017 
1018  } );
1019 
1020  }
1021 
1022 
1031  template< typename POLICY, typename KERNEL_TYPE >
1032  static void
1033  launch( localIndex const numElements,
1034  KERNEL_TYPE const & kernelComponent )
1035  {
1037  forAll< POLICY >( numElements, [=] GEOS_HOST_DEVICE ( localIndex const ie )
1038  {
1039  typename KERNEL_TYPE::StackVariables stack( 1 );
1040 
1041  kernelComponent.setup( ie, stack );
1042  kernelComponent.computeFlux( ie, stack );
1043  kernelComponent.complete( ie, stack );
1044  } );
1045  }
1046 
1047 protected:
1050 
1053 
1057 
1061 
1062 
1063 };
1064 
1069 {
1070 public:
1071 
1085  template< typename POLICY >
1086  static void
1087  createAndLaunch( integer const numComps,
1088  real64 const dt,
1089  globalIndex const rankOffset,
1090  BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags,
1091  string const dofKey,
1092  WellControls const & wellControls,
1093  WellElementSubRegion const & subRegion,
1094  MultiFluidBase const & fluid,
1095  CRSMatrixView< real64, globalIndex const > const & localMatrix,
1096  arrayView1d< real64 > const & localRhs )
1097  {
1098  isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComps, [&]( auto NC )
1099  {
1100  integer constexpr NUM_COMP = NC();
1101 
1102  using kernelType = FaceBasedAssemblyKernel< NUM_COMP >;
1103  kernelType kernel( dt, rankOffset, dofKey, wellControls, subRegion, fluid, localMatrix, localRhs, kernelFlags );
1104  kernelType::template launch< POLICY >( subRegion.size(), kernel );
1105  } );
1106  }
1107 };
1108 
1109 } // end namespace thermalCompositionalMultiphaseWellKernels
1110 
1111 } // end namespace geos
1112 
1113 #endif // GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_THERMALCOMPOSITIONALMULTIPHASEWELLKERNELS_HPP
#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.
Definition: GeosxMacros.hpp:84
#define GEOS_ERROR(msg)
Raise a hard error and terminate the program.
Definition: Logger.hpp:157
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
The ObjectManagerBase is the base object of all object managers in the mesh data hierachy.
array1d< integer > const & ghostRank()
Get the ghost information of each object.
This class describes the controls used to operate a well.
This class describes a collection of local well elements and perforations.
Define the interface for the assembly kernel in charge of accumulation and volume balance.
Define the interface for the assembly kernel in charge of flux terms.
Define the interface for the property kernel in charge of computing the total mass density.
GEOS_DECLTYPE_AUTO_RETURN getReference(LOOKUP_TYPE const &lookup) const
Look up a wrapper and get reference to wrapped object.
Definition: Group.hpp:1275
localIndex size() const
Get the "size" of the group, which determines the number of elements in resizable wrappers.
Definition: Group.hpp:1317
Define the base interface for the property update kernels.
Define the base interface for the residual calculations.
static void createAndLaunch(localIndex const numComps, localIndex const numPhases, integer const isProducer, globalIndex const rankOffset, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags, string const dofKey, WellElementSubRegion const &subRegion, MultiFluidBase const &fluid, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs)
Create a new kernel and launch.
Define the interface for the assembly kernel in charge of thermal accumulation and volume balance.
arrayView2d< real64 const > const m_dPoro_dTemp
View on derivative of porosity w.r.t temperature.
GEOS_HOST_DEVICE void computeVolumeBalance(localIndex const ei, StackVariables &stack) const
Compute the local volume balance contributions to the residual and Jacobian.
arrayView3d< real64 const, multifluid::USD_PHASE > m_phaseInternalEnergy_n
Views on phase internal energy.
GEOS_HOST_DEVICE void computeAccumulation(localIndex const ei, StackVariables &stack) const
Compute the local accumulation contributions to the residual and Jacobian.
ElementBasedAssemblyKernel(localIndex const numPhases, integer const isProducer, globalIndex const rankOffset, string const dofKey, WellElementSubRegion const &subRegion, MultiFluidBase const &fluid, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > const kernelFlags)
Constructor.
GEOS_HOST_DEVICE void setup(localIndex const ei, StackVariables &stack) const
Performs the setup phase for the kernel.
static void createAndLaunch(integer const numComps, real64 const dt, globalIndex const rankOffset, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags, string const dofKey, WellControls const &wellControls, WellElementSubRegion const &subRegion, MultiFluidBase const &fluid, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs)
Create a new kernel and launch.
Define the interface for the assembly kernel in charge of flux terms.
FaceBasedAssemblyKernel(real64 const dt, globalIndex const rankOffset, string const wellDofKey, WellControls const &wellControls, WellElementSubRegion const &subRegion, MultiFluidBase const &fluid, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags)
Constructor for the kernel interface.
arrayView3d< real64 const, multifluid::USD_PHASE > m_phaseEnthalpy
Views on phase enthalpy.
GEOS_HOST_DEVICE void computeFlux(localIndex const iwelem, StackVariables &stack) const
Compute the local flux contributions to the residual and Jacobian.
arrayView1d< globalIndex const > m_globalWellElementIndex
Global index of local element.
arrayView3d< real64 const, multifluid::USD_PHASE > const m_phaseFraction
Element phase fraction.
static void launch(localIndex const numElements, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
static void createAndLaunch(integer const numComp, integer const targetPhaseIndex, globalIndex const rankOffset, string const &dofKey, arrayView1d< real64 const > const &localResidual, WellElementSubRegion const &subRegion, MultiFluidBase const &fluid, WellControls const &wellControls, real64 const time, real64 const dt, real64 const minNormalizer, real64(&residualNorm)[2])
Create a new kernel and launch.
arrayView3d< real64 const, multifluid::USD_PHASE > const m_phaseDens_n
View on phase/total density at the previous converged time step.
localIndex const m_iwelemControl
Index of the element where the control is enforced.
virtual GEOS_HOST_DEVICE void computeLinf(localIndex const iwelem, LinfStackVariables &stack) const override
Compute the local values for the Linf norm.
virtual GEOS_HOST_DEVICE void computeL2(localIndex const iwelem, L2StackVariables &stack) const override
Compute the local values and normalizer for the L2 norm.
bool const m_isProducer
Flag indicating whether the well is a producer or an injector.
bool const m_isLocallyOwned
Flag indicating whether the well is locally owned or not.
static void createAndLaunch(integer const numComp, integer const numPhase, ObjectManagerBase &subRegion, MultiFluidBase const &fluid)
Create a new kernel and launch.
Define the interface for the property kernel in charge of computing the total mass density.
GEOS_HOST_DEVICE void compute(localIndex const ei) const
Compute the total mass density in an element.
TotalMassDensityKernel(ObjectManagerBase &subRegion, MultiFluidBase const &fluid)
Constructor.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:188
StackArray< T, 2, MAXSIZE > stackArray2d
Alias for 2D stack array.
Definition: DataTypes.hpp:212
LvArray::CRSMatrixView< T, COL_INDEX, localIndex const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
Definition: DataTypes.hpp:318
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:196
ArraySlice< T, 2, USD > arraySlice2d
Alias for 2D array slice.
Definition: DataTypes.hpp:208
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
ArraySlice< T, 1, USD > arraySlice1d
Alias for 1D array slice.
Definition: DataTypes.hpp:192
std::int32_t integer
Signed integer type.
Definition: DataTypes.hpp:81
ArrayView< T, 4, USD > arrayView4d
Alias for 4D array view.
Definition: DataTypes.hpp:236
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:204
Array< T, 1 > array1d
Alias for 1D array.
Definition: DataTypes.hpp:184
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
Definition: DataTypes.hpp:220
Kernel variables (dof numbers, jacobian and residual) located on the stack.
stackArray1d< real64, maxNumElems > localEnergyFlux
Storage for the face local residual vector (energy equation)
stackArray2d< real64, maxNumElems *maxStencilSize > localEnergyFluxJacobian_dQ
Storage for the face local Jacobian matrix dQ only.
stackArray2d< real64, maxNumElems *maxStencilSize *CP_Deriv::nDer > localEnergyFluxJacobian
Storage for the face local energy Jacobian matrix dC dP dT.