GEOS
ThermalFluxComputeKernel.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 Total, S.A
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_COMPOSITIONAL_THERMALFLUXCOMPUTEKERNEL_HPP
21 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_THERMALFLUXCOMPUTEKERNEL_HPP
22 
24 #include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
25 #include "constitutive/fluid/multifluid/MultiFluidFields.hpp"
26 #include "constitutive/thermalConductivity/MultiPhaseThermalConductivityBase.hpp"
27 #include "constitutive/thermalConductivity/ThermalConductivityFields.hpp"
28 
29 namespace geos
30 {
31 
32 namespace thermalCompositionalMultiphaseFVMKernels
33 {
34 
35 /******************************** FluxComputeKernel ********************************/
36 
44 template< integer NUM_COMP, integer NUM_DOF, typename STENCILWRAPPER >
45 class FluxComputeKernel : public isothermalCompositionalMultiphaseFVMKernels::FluxComputeKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >
46 {
47 public:
48 
55  template< typename VIEWTYPE >
57 
59  using DofNumberAccessor = AbstractBase::DofNumberAccessor;
64 
65  using AbstractBase::m_dt;
69  using AbstractBase::m_gravCoef;
71  using AbstractBase::m_dPhaseVolFrac;
73  using AbstractBase::m_dPhaseCompFrac;
75 
77  using Base::numComp;
78  using Base::numDof;
79  using Base::numEqn;
80  using Base::maxNumElems;
81  using Base::maxNumConns;
84  using Base::m_phaseMob;
85  using Base::m_dPhaseMob;
86  using Base::m_dPhaseMassDens;
87  using Base::m_dPhaseCapPressure_dPhaseVolFrac;
89  using Base::m_seri;
90  using Base::m_sesri;
91  using Base::m_sei;
92 
95 
97  StencilMaterialAccessors< constitutive::MultiFluidBase,
98  fields::multifluid::phaseEnthalpy,
99  fields::multifluid::dPhaseEnthalpy >;
100 
102  StencilMaterialAccessors< constitutive::MultiPhaseThermalConductivityBase,
103  fields::thermalconductivity::effectiveConductivity >;
104  // for now, we treat thermal conductivity explicitly
105 
124  FluxComputeKernel( integer const numPhases,
125  globalIndex const rankOffset,
126  STENCILWRAPPER const & stencilWrapper,
127  DofNumberAccessor const & dofNumberAccessor,
128  CompFlowAccessors const & compFlowAccessors,
129  ThermalCompFlowAccessors const & thermalCompFlowAccessors,
130  MultiFluidAccessors const & multiFluidAccessors,
131  ThermalMultiFluidAccessors const & thermalMultiFluidAccessors,
132  CapPressureAccessors const & capPressureAccessors,
133  PermeabilityAccessors const & permeabilityAccessors,
134  ThermalConductivityAccessors const & thermalConductivityAccessors,
135  real64 const dt,
136  CRSMatrixView< real64, globalIndex const > const & localMatrix,
137  arrayView1d< real64 > const & localRhs,
138  BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags )
139  : Base( numPhases,
140  rankOffset,
141  stencilWrapper,
142  dofNumberAccessor,
143  compFlowAccessors,
144  multiFluidAccessors,
145  capPressureAccessors,
146  permeabilityAccessors,
147  dt,
148  localMatrix,
149  localRhs,
150  kernelFlags ),
151  m_temp( thermalCompFlowAccessors.get( fields::flow::temperature {} ) ),
152  m_phaseEnthalpy( thermalMultiFluidAccessors.get( fields::multifluid::phaseEnthalpy {} ) ),
153  m_dPhaseEnthalpy( thermalMultiFluidAccessors.get( fields::multifluid::dPhaseEnthalpy {} ) ),
154  m_thermalConductivity( thermalConductivityAccessors.get( fields::thermalconductivity::effectiveConductivity {} ) )
155  {}
156 
158  {
159 public:
160 
162  StackVariables( localIndex const size, localIndex numElems )
163  : Base::StackVariables( size, numElems )
164  {}
165 
173 
174  // Thermal transmissibility (for now, no derivatives)
175 
176  real64 thermalTransmissibility[maxNumConns][2]{};
177  };
178 
185  inline
186  void computeFlux( localIndex const iconn,
187  StackVariables & stack ) const
188  {
189  using Deriv = constitutive::multifluid::DerivativeOffset;
190 
191  // ***********************************************
192  // First, we call the base computeFlux to compute:
193  // 1) compFlux and its derivatives (including derivatives wrt temperature),
194  // 2) enthalpy part of convectiveEnergyFlux and its derivatives (including derivatives wrt temperature)
195  //
196  // Computing dCompFlux_dT and the enthalpy flux requires quantities already computed in the base computeFlux,
197  // such as potGrad, phaseFlux, and the indices of the upwind cell
198  // We use the lambda below (called **inside** the phase loop of the base computeFlux) to access these variables
199  Base::computeFlux( iconn, stack, [&] ( integer const ip,
200  localIndex const (&k)[2],
201  localIndex const (&seri)[2],
202  localIndex const (&sesri)[2],
203  localIndex const (&sei)[2],
204  localIndex const connectionIndex,
205  localIndex const k_up,
206  localIndex const er_up,
207  localIndex const esr_up,
208  localIndex const ei_up,
209  real64 const potGrad,
210  real64 const phaseFlux,
211  real64 const (&dPhaseFlux_dP)[2],
212  real64 const (&dPhaseFlux_dC)[2][numComp] )
213  {
214  // We are in the loop over phases, ip provides the current phase index.
215 
216  // Step 1: compute the derivatives of the mean density at the interface wrt temperature
217 
218  real64 dDensMean_dT[numFluxSupportPoints]{};
219 
220  real64 const trans[numFluxSupportPoints] = { stack.transmissibility[connectionIndex][0],
221  stack.transmissibility[connectionIndex][1] };
222 
223  real64 convectiveEnergyFlux = 0.0;
224  real64 dConvectiveEnergyFlux_dP[numFluxSupportPoints]{};
225  real64 dConvectiveEnergyFlux_dT[numFluxSupportPoints]{};
226  real64 dConvectiveEnergyFlux_dC[numFluxSupportPoints][numComp]{};
227  real64 dCompFlux_dT[numFluxSupportPoints][numComp]{};
228 
229  integer denom = 0;
230  for( integer i = 0; i < numFluxSupportPoints; ++i )
231  {
232  localIndex const er = seri[i];
233  localIndex const esr = sesri[i];
234  localIndex const ei = sei[i];
235 
236  bool const phaseExists = (m_phaseVolFrac[er_up][esr_up][ei_up][ip] > 0);
237  if( !phaseExists )
238  {
239  continue;
240  }
241 
242  dDensMean_dT[i] = m_dPhaseMassDens[er][esr][ei][0][ip][Deriv::dT];
243  denom++;
244  }
245  if( denom > 1 )
246  {
247  for( integer i = 0; i < numFluxSupportPoints; ++i )
248  {
249  dDensMean_dT[i] /= denom;
250  }
251  }
252 
253  // Step 2: compute the derivatives of the phase potential difference wrt temperature
254  //***** calculation of flux *****
255 
256  real64 dPresGrad_dT[numFluxSupportPoints]{};
257  real64 dGravHead_dT[numFluxSupportPoints]{};
258 
259  // compute potential difference MPFA-style
260  for( integer i = 0; i < numFluxSupportPoints; ++i )
261  {
262  localIndex const er = seri[i];
263  localIndex const esr = sesri[i];
264  localIndex const ei = sei[i];
265 
266  // Step 2.1: compute derivative of capillary pressure wrt temperature
267  real64 dCapPressure_dT = 0.0;
268  if( AbstractBase::m_kernelFlags.isSet( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::CapPressure ) )
269  {
270  for( integer jp = 0; jp < m_numPhases; ++jp )
271  {
272  real64 const dCapPressure_dS = m_dPhaseCapPressure_dPhaseVolFrac[er][esr][ei][0][ip][jp];
273  dCapPressure_dT += dCapPressure_dS * m_dPhaseVolFrac[er][esr][ei][jp][Deriv::dT];
274  }
275  }
276 
277  // Step 2.2: compute derivative of phase pressure difference wrt temperature
278  dPresGrad_dT[i] -= trans[i] * dCapPressure_dT;
279  real64 const gravD = trans[i] * m_gravCoef[er][esr][ei];
280 
281  // Step 2.3: compute derivative of gravity potential difference wrt temperature
282  for( integer j = 0; j < numFluxSupportPoints; ++j )
283  {
284  dGravHead_dT[j] += dDensMean_dT[j] * gravD;
285  }
286  }
287 
288  // Step 3: compute the derivatives of the (upwinded) compFlux wrt temperature
289  // *** upwinding ***
290 
291  // note: the upwinding is done in the base class, which is in charge of
292  // computing the following quantities: potGrad, phaseFlux, k_up, er_up, esr_up, ei_up
293 
294  real64 dPhaseFlux_dT[numFluxSupportPoints]{};
295 
296  // Step 3.1: compute the derivative of phase flux wrt temperature
297  for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
298  {
299  dPhaseFlux_dT[ke] += dPresGrad_dT[ke];
300  }
301  for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
302  {
303  dPhaseFlux_dT[ke] -= dGravHead_dT[ke];
304  }
305  for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
306  {
307  dPhaseFlux_dT[ke] *= m_phaseMob[er_up][esr_up][ei_up][ip];
308  }
309  dPhaseFlux_dT[k_up] += m_dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dT] * potGrad;
310 
311  // Step 3.2: compute the derivative of component flux wrt temperature
312 
313  // slice some constitutive arrays to avoid too much indexing in component loop
314  arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE_COMP - 3 > phaseCompFracSub =
315  m_phaseCompFrac[er_up][esr_up][ei_up][0][ip];
316  arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC - 3 > dPhaseCompFracSub =
317  m_dPhaseCompFrac[er_up][esr_up][ei_up][0][ip];
318 
319  for( integer ic = 0; ic < numComp; ++ic )
320  {
321  real64 const ycp = phaseCompFracSub[ic];
322  for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
323  {
324  dCompFlux_dT[ke][ic] += dPhaseFlux_dT[ke] * ycp;
325  }
326  dCompFlux_dT[k_up][ic] += phaseFlux * dPhaseCompFracSub[ic][Deriv::dT];
327  }
328 
329  // Step 4: add dCompFlux_dTemp to localFluxJacobian
330  for( integer ic = 0; ic < numComp; ++ic )
331  {
332  integer const eqIndex0 = k[0]* numEqn + ic;
333  integer const eqIndex1 = k[1]* numEqn + ic;
334  for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
335  {
336  integer const localDofIndexTemp = k[ke] * numDof + numDof - 1;
337  stack.localFluxJacobian[eqIndex0][localDofIndexTemp] += m_dt * dCompFlux_dT[ke][ic];
338  stack.localFluxJacobian[eqIndex1][localDofIndexTemp] -= m_dt * dCompFlux_dT[ke][ic];
339  }
340  }
341 
342  // Step 5: compute the enthalpy flux
343  real64 const enthalpy = m_phaseEnthalpy[er_up][esr_up][ei_up][0][ip];
344  convectiveEnergyFlux += phaseFlux * enthalpy;
345 
346  for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
347  {
348  dConvectiveEnergyFlux_dP[ke] += dPhaseFlux_dP[ke] * enthalpy;
349  dConvectiveEnergyFlux_dT[ke] += dPhaseFlux_dT[ke] * enthalpy;
350 
351  for( integer jc = 0; jc < numComp; ++jc )
352  {
353  dConvectiveEnergyFlux_dC[ke][jc] += dPhaseFlux_dC[ke][jc] * enthalpy;
354  }
355  }
356 
357  dConvectiveEnergyFlux_dP[k_up] += phaseFlux * m_dPhaseEnthalpy[er_up][esr_up][ei_up][0][ip][Deriv::dP];
358  dConvectiveEnergyFlux_dT[k_up] += phaseFlux * m_dPhaseEnthalpy[er_up][esr_up][ei_up][0][ip][Deriv::dT];
359 
360  real64 dProp_dC[numComp]{};
361  applyChainRule( numComp,
362  m_dCompFrac_dCompDens[er_up][esr_up][ei_up],
363  m_dPhaseEnthalpy[er_up][esr_up][ei_up][0][ip],
364  dProp_dC,
365  Deriv::dC );
366  for( integer jc = 0; jc < numComp; ++jc )
367  {
368  dConvectiveEnergyFlux_dC[k_up][jc] += phaseFlux * dProp_dC[jc];
369  }
370 
371  // Step 6: add convectiveFlux and its derivatives to localFlux and localFluxJacobian
372  integer const localRowIndexEnergy0 = k[0] * numEqn + numEqn - 1;
373  integer const localRowIndexEnergy1 = k[1] * numEqn + numEqn - 1;
374  stack.localFlux[localRowIndexEnergy0] += m_dt * convectiveEnergyFlux;
375  stack.localFlux[localRowIndexEnergy1] -= m_dt * convectiveEnergyFlux;
376 
377  for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
378  {
379  integer const localDofIndexPres = k[ke] * numDof;
380  stack.localFluxJacobian[localRowIndexEnergy0][localDofIndexPres] += m_dt * dConvectiveEnergyFlux_dP[ke];
381  stack.localFluxJacobian[localRowIndexEnergy1][localDofIndexPres] -= m_dt * dConvectiveEnergyFlux_dP[ke];
382  integer const localDofIndexTemp = localDofIndexPres + numDof - 1;
383  stack.localFluxJacobian[localRowIndexEnergy0][localDofIndexTemp] += m_dt * dConvectiveEnergyFlux_dT[ke];
384  stack.localFluxJacobian[localRowIndexEnergy1][localDofIndexTemp] -= m_dt * dConvectiveEnergyFlux_dT[ke];
385 
386  for( integer jc = 0; jc < numComp; ++jc )
387  {
388  integer const localDofIndexComp = localDofIndexPres + jc + 1;
389  stack.localFluxJacobian[localRowIndexEnergy0][localDofIndexComp] += m_dt * dConvectiveEnergyFlux_dC[ke][jc];
390  stack.localFluxJacobian[localRowIndexEnergy1][localDofIndexComp] -= m_dt * dConvectiveEnergyFlux_dC[ke][jc];
391  }
392  }
393  } );
394 
395  // *****************************************************
396  // Computation of the conduction term in the energy flux
397  // Note that the phase enthalpy term in the energy was computed above
398  // Note that this term is computed using an explicit treatment of conductivity for now
399 
400  // Step 1: compute the thermal transmissibilities at this face
401  // Below, the thermal conductivity used to compute (explicitly) the thermal conducivity
402  // To avoid modifying the signature of the "computeWeights" function for now, we pass m_thermalConductivity twice
403  // TODO: modify computeWeights to accomodate explicit coefficients
404  m_stencilWrapper.computeWeights( iconn,
406  m_thermalConductivity, // we have to pass something here, so we just use thermal conductivity
407  stack.thermalTransmissibility,
408  stack.dTrans_dPres ); // again, we have to pass something here, but this is unused for now
409 
410 
411 
412  localIndex k[2]{};
413  localIndex connectionIndex = 0;
414 
415  for( k[0] = 0; k[0] < stack.numConnectedElems; ++k[0] )
416  {
417  for( k[1] = k[0] + 1; k[1] < stack.numConnectedElems; ++k[1] )
418  {
419  real64 const thermalTrans[2] = { stack.thermalTransmissibility[connectionIndex][0], stack.thermalTransmissibility[connectionIndex][1] };
420  localIndex const seri[2] = {m_seri( iconn, k[0] ), m_seri( iconn, k[1] )};
421  localIndex const sesri[2] = {m_sesri( iconn, k[0] ), m_sesri( iconn, k[1] )};
422  localIndex const sei[2] = {m_sei( iconn, k[0] ), m_sei( iconn, k[1] )};
423 
424  real64 conductiveEnergyFlux = 0.0;
425  real64 dConductiveEnergyFlux_dT[numFluxSupportPoints]{};
426 
427  // Step 2: compute temperature difference at the interface
428  for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
429  {
430  localIndex const er = seri[ke];
431  localIndex const esr = sesri[ke];
432  localIndex const ei = sei[ke];
433 
434  conductiveEnergyFlux += thermalTrans[ke] * m_temp[er][esr][ei];
435  dConductiveEnergyFlux_dT[ke] += thermalTrans[ke];
436  }
437 
438  // Step 3: add conductiveFlux and its derivatives to localFlux and localFluxJacobian
439  integer const localRowIndexEnergy0 = k[0] * numEqn + numEqn - 1;
440  integer const localRowIndexEnergy1 = k[1] * numEqn + numEqn - 1;
441  stack.localFlux[localRowIndexEnergy0] += m_dt * conductiveEnergyFlux;
442  stack.localFlux[localRowIndexEnergy1] -= m_dt * conductiveEnergyFlux;
443 
444  for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
445  {
446  integer const localDofIndexTemp = k[ke] * numDof + numDof - 1;
447  stack.localFluxJacobian[localRowIndexEnergy0][localDofIndexTemp] += m_dt * dConductiveEnergyFlux_dT[ke];
448  stack.localFluxJacobian[localRowIndexEnergy1][localDofIndexTemp] -= m_dt * dConductiveEnergyFlux_dT[ke];
449  }
450  }
451  }
452  }
453 
460  inline
461  void complete( localIndex const iconn,
462  StackVariables & stack ) const
463  {
464  // Call Case::complete to assemble the component mass balance equations (i = 0 to i = numDof-2)
465  // In the lambda, add contribution to residual and jacobian into the energy balance equation
466  Base::complete( iconn, stack, [&] ( integer const i,
467  localIndex const localRow )
468  {
469  // beware, there is volume balance eqn in m_localRhs and m_localMatrix!
470  RAJA::atomicAdd( parallelDeviceAtomic{}, &AbstractBase::m_localRhs[localRow + numEqn], stack.localFlux[i * numEqn + numEqn-1] );
471  AbstractBase::m_localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >
472  ( localRow + numEqn,
473  stack.dofColIndices.data(),
474  stack.localFluxJacobian[i * numEqn + numEqn-1].dataIfContiguous(),
475  stack.stencilSize * numDof );
476 
477  } );
478  }
479 
480 protected:
481 
484 
488 
491  // for now, we treat thermal conductivity explicitly
492 
493 };
494 
499 {
500 public:
501 
518  template< typename POLICY, typename STENCILWRAPPER >
519  static void
520  createAndLaunch( integer const numComps,
521  integer const numPhases,
522  globalIndex const rankOffset,
523  string const & dofKey,
524  integer const hasCapPressure,
525  integer const useTotalMassEquation,
526  string const & solverName,
527  ElementRegionManager const & elemManager,
528  STENCILWRAPPER const & stencilWrapper,
529  real64 const dt,
530  CRSMatrixView< real64, globalIndex const > const & localMatrix,
531  arrayView1d< real64 > const & localRhs )
532  {
533  isothermalCompositionalMultiphaseBaseKernels::
534  internal::kernelLaunchSelectorCompSwitch( numComps, [&]( auto NC )
535  {
536  integer constexpr NUM_COMP = NC();
537  integer constexpr NUM_DOF = NC() + 2;
538 
540  elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey );
541  dofNumberAccessor.setName( solverName + "/accessors/" + dofKey );
542 
543  BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags;
544  if( hasCapPressure )
545  kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::CapPressure );
546  if( useTotalMassEquation )
547  kernelFlags.set( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::TotalMassEquation );
548 
550  typename KernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName );
551  typename KernelType::ThermalCompFlowAccessors thermalCompFlowAccessors( elemManager, solverName );
552  typename KernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName );
553  typename KernelType::ThermalMultiFluidAccessors thermalMultiFluidAccessors( elemManager, solverName );
554  typename KernelType::CapPressureAccessors capPressureAccessors( elemManager, solverName );
555  typename KernelType::PermeabilityAccessors permeabilityAccessors( elemManager, solverName );
556  typename KernelType::ThermalConductivityAccessors thermalConductivityAccessors( elemManager, solverName );
557 
558  KernelType kernel( numPhases, rankOffset, stencilWrapper, dofNumberAccessor,
559  compFlowAccessors, thermalCompFlowAccessors, multiFluidAccessors, thermalMultiFluidAccessors,
560  capPressureAccessors, permeabilityAccessors, thermalConductivityAccessors,
561  dt, localMatrix, localRhs, kernelFlags );
562  KernelType::template launch< POLICY >( stencilWrapper.size(), kernel );
563  } );
564  }
565 };
566 
567 } // namespace thermalCompositionalMultiphaseFVMKernels
568 
569 } // namespace geos
570 
571 
572 #endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_THERMALFLUXCOMPUTEKERNEL_HPP
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
The ElementRegionManager class provides an interface to ObjectManagerBase in order to manage ElementR...
ElementViewAccessor< ArrayView< T const, NDIM, getUSD< PERM > > > constructArrayViewAccessor(string const &name, string const &neighborName=string()) const
This is a function to construct a ElementViewAccessor to access array data registered on the mesh.
array1d< array1d< VIEWTYPE > > ElementViewAccessor
The ElementViewAccessor at the ElementRegionManager level is an array of array of VIEWTYPE.
typename ElementViewAccessor< VIEWTYPE >::NestedViewTypeConst ElementViewConst
The ElementViewAccessor at the ElementRegionManager level is the type resulting from ElementViewAcces...
A struct to automatically construct and store element view accessors.
A struct to automatically construct and store element view accessors.
Base class for FluxComputeKernel that holds all data not dependent on template parameters (like stenc...
ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const m_phaseVolFrac
Views on phase volume fractions.
ElementViewConst< arrayView1d< globalIndex const > > const m_dofNumber
Views on dof numbers.
ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const m_phaseCompFrac
Views on phase component fractions.
CRSMatrixView< real64, globalIndex const > const m_localMatrix
View on the local CRS matrix.
ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const m_dCompFrac_dCompDens
Views on derivatives of comp fractions.
Define the interface for the assembly kernel in charge of flux terms.
static constexpr localIndex maxStencilSize
Maximum number of points in the stencil.
GEOS_HOST_DEVICE void computeFlux(localIndex const iconn, StackVariables &stack, FUNC &&compFluxKernelOp=NoOpFunc{}) const
Compute the local flux contributions to the residual and Jacobian.
GEOS_HOST_DEVICE void complete(localIndex const iconn, StackVariables &stack, FUNC &&assemblyKernelOp=NoOpFunc{}) const
Performs the complete phase for the kernel.
static constexpr integer numFluxSupportPoints
Number of flux support points (hard-coded for TFPA)
static constexpr integer numEqn
Compute time value for the number of equations (all of them, except the volume balance equation)
static constexpr integer numComp
Compile time value for the number of components.
ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const m_phaseMob
Views on phase mobilities.
STENCILWRAPPER::IndexContainerViewConstType const m_seri
Connection to element maps.
static constexpr localIndex maxNumElems
Maximum number of elements at the face.
static constexpr integer numDof
Compute time value for the number of degrees of freedom.
STENCILWRAPPER const m_stencilWrapper
Reference to the stencil wrapper.
static constexpr localIndex maxNumConns
Maximum number of connections at the face.
static void createAndLaunch(integer const numComps, integer const numPhases, globalIndex const rankOffset, string const &dofKey, integer const hasCapPressure, integer const useTotalMassEquation, string const &solverName, ElementRegionManager const &elemManager, STENCILWRAPPER const &stencilWrapper, real64 const dt, 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.
ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const m_phaseEnthalpy
Views on phase enthalpies.
ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const m_phaseVolFrac
Views on phase volume fractions.
static constexpr integer numFluxSupportPoints
Number of flux support points (hard-coded for TFPA)
ElementViewConst< arrayView3d< real64 const > > const m_thermalConductivity
View on thermal conductivity.
static constexpr integer numEqn
Compute time value for the number of equations (all of them, except the volume balance equation)
static constexpr integer numComp
Compile time value for the number of components.
ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const m_phaseMob
Views on phase mobilities.
GEOS_HOST_DEVICE void computeFlux(localIndex const iconn, StackVariables &stack) const
Compute the local flux contributions to the residual and Jacobian.
GEOS_HOST_DEVICE void complete(localIndex const iconn, StackVariables &stack) const
Performs the complete phase for the kernel.
ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const m_phaseCompFrac
Views on phase component fractions.
ElementRegionManager::ElementViewConst< VIEWTYPE > ElementViewConst
The type for element-based data. Consists entirely of ArrayView's.
STENCILWRAPPER::IndexContainerViewConstType const m_seri
Connection to element maps.
static constexpr integer numDof
Compute time value for the number of degrees of freedom.
ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const m_dCompFrac_dCompDens
Views on derivatives of comp fractions.
STENCILWRAPPER const m_stencilWrapper
Reference to the stencil wrapper.
ElementViewConst< arrayView1d< real64 const > > const m_temp
Views on temperature.
static constexpr localIndex maxNumConns
Maximum number of connections at the face.
FluxComputeKernel(integer const numPhases, globalIndex const rankOffset, STENCILWRAPPER const &stencilWrapper, DofNumberAccessor const &dofNumberAccessor, CompFlowAccessors const &compFlowAccessors, ThermalCompFlowAccessors const &thermalCompFlowAccessors, MultiFluidAccessors const &multiFluidAccessors, ThermalMultiFluidAccessors const &thermalMultiFluidAccessors, CapPressureAccessors const &capPressureAccessors, PermeabilityAccessors const &permeabilityAccessors, ThermalConductivityAccessors const &thermalConductivityAccessors, real64 const dt, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs, BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags)
Constructor for the kernel interface.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:180
LvArray::CRSMatrixView< T, COL_INDEX, localIndex const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
Definition: DataTypes.hpp:310
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
Definition: DataTypes.hpp:88
ArraySlice< T, 2, USD > arraySlice2d
Alias for 2D array slice.
Definition: DataTypes.hpp:200
double real64
64-bit floating point type.
Definition: DataTypes.hpp:99
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:85
ArraySlice< T, 1, USD > arraySlice1d
Alias for 1D array slice.
Definition: DataTypes.hpp:184
std::int32_t integer
Signed integer type.
Definition: DataTypes.hpp:82
Kernel variables (dof numbers, jacobian and residual) located on the stack.
real64 dTrans_dPres[maxNumConns][numFluxSupportPoints]
Derivatives of transmissibility with respect to pressure.
stackArray1d< globalIndex, maxNumElems *numDof > dofColIndices
Indices of the matrix rows/columns corresponding to the dofs in this face.
stackArray2d< real64, maxNumElems *numEqn *maxStencilSize *numDof > localFluxJacobian
Storage for the face local Jacobian matrix.
stackArray1d< real64, maxNumElems *numEqn > localFlux
Storage for the face local residual vector (all equations except volume balance)
localIndex const numConnectedElems
Number of elements connected at a given connection.
real64 transmissibility[maxNumConns][numFluxSupportPoints]
Transmissibility.