GEOS
AccumulationKernel.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_FLUIDFLOW_COMPOSITIONAL_ACCUMULATIONKERNEL_HPP
21 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_ACCUMULATIONKERNEL_HPP
22 
23 #include "codingUtilities/Utilities.hpp"
24 #include "common/DataLayouts.hpp"
25 #include "common/DataTypes.hpp"
26 #include "common/GEOS_RAJA_Interface.hpp"
27 #include "constitutive/solid/CoupledSolidBase.hpp"
28 #include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
33 #include "physicsSolvers/fluidFlow/kernels/compositional/KernelLaunchSelectors.hpp"
34 
35 namespace geos
36 {
37 
38 namespace isothermalCompositionalMultiphaseBaseKernels
39 {
40 
41 static constexpr real64 minDensForDivision = 1e-10;
42 
43 enum class KernelFlags
44 {
45  SimpleAccumulation = 1 << 0, // 1
46  TotalMassEquation = 1 << 1, // 2
48  // Flag3 = 1 << 2, // 4
49  // Flag4 = 1 << 3, // 8
50  // Flag5 = 1 << 4, // 16
51  // Flag6 = 1 << 5, // 32
52  // Flag7 = 1 << 6, // 64
53  // Flag8 = 1 << 7 //128
54 };
55 
56 /******************************** AccumulationKernel ********************************/
57 
64 template< integer NUM_COMP, integer NUM_DOF >
66 {
67 public:
68 
70  static constexpr integer numComp = NUM_COMP;
71 
73  static constexpr integer numDof = NUM_DOF;
74 
76  static constexpr integer numEqn = NUM_DOF;
77 
89  AccumulationKernel( localIndex const numPhases,
90  globalIndex const rankOffset,
91  string const dofKey,
92  ElementSubRegionBase const & subRegion,
93  constitutive::MultiFluidBase const & fluid,
94  constitutive::CoupledSolidBase const & solid,
96  arrayView1d< real64 > const & localRhs,
97  BitFlags< KernelFlags > const kernelFlags )
98  : m_numPhases( numPhases ),
99  m_rankOffset( rankOffset ),
100  m_dofNumber( subRegion.getReference< array1d< globalIndex > >( dofKey ) ),
101  m_elemGhostRank( subRegion.ghostRank() ),
102  m_volume( subRegion.getElementVolume() ),
103  m_porosity( solid.getPorosity() ),
104  m_porosity_n( solid.getPorosity_n() ),
105  m_dPoro_dPres( solid.getDporosity_dPressure() ),
106  m_dCompFrac_dCompDens( subRegion.getField< fields::flow::dGlobalCompFraction_dGlobalCompDensity >() ),
107  m_phaseVolFrac( subRegion.getField< fields::flow::phaseVolumeFraction >() ),
108  m_dPhaseVolFrac( subRegion.getField< fields::flow::dPhaseVolumeFraction >() ),
109  m_phaseDens( fluid.phaseDensity() ),
110  m_dPhaseDens( fluid.dPhaseDensity() ),
111  m_phaseCompFrac( fluid.phaseCompFraction() ),
112  m_dPhaseCompFrac( fluid.dPhaseCompFraction() ),
113  m_compDens( subRegion.getField< fields::flow::globalCompDensity >() ),
114  m_compAmount_n( subRegion.getField< fields::flow::compAmount_n >() ),
115  m_localMatrix( localMatrix ),
116  m_localRhs( localRhs ),
117  m_kernelFlags( kernelFlags )
118  {}
119 
125  {
126 public:
127 
128  // Pore volume information (used by both accumulation and volume balance)
129 
132 
135 
136  // Residual information
137 
140 
143 
146 
149 
150  };
151 
158  integer elemGhostRank( localIndex const ei ) const
159  { return m_elemGhostRank( ei ); }
160 
161 
168  void setup( localIndex const ei,
169  StackVariables & stack ) const
170  {
171  // initialize the pore volume
172  stack.poreVolume = m_volume[ei] * m_porosity[ei][0];
173  stack.dPoreVolume_dPres = m_volume[ei] * m_dPoro_dPres[ei][0];
174 
175  // set row index and degrees of freedom indices for this element
176  stack.localRow = m_dofNumber[ei] - m_rankOffset;
177  for( integer idof = 0; idof < numDof; ++idof )
178  {
179  stack.dofIndices[idof] = m_dofNumber[ei] + idof;
180  }
181  }
182 
190  template< typename FUNC = NoOpFunc >
193  StackVariables & stack,
194  FUNC && phaseAmountKernelOp = NoOpFunc{} ) const
195  {
196  if( m_kernelFlags.isSet( KernelFlags::SimpleAccumulation ) )
197  {
198  // ic - index of component whose conservation equation is assembled
199  // (i.e. row number in local matrix)
200  for( integer ic = 0; ic < numComp; ++ic )
201  {
202  real64 const compAmount = stack.poreVolume * m_compDens[ei][ic];
203  real64 const compAmount_n = m_compAmount_n[ei][ic];
204 
205  stack.localResidual[ic] += compAmount - compAmount_n;
206 
207  // Pavel: commented below is some experiment, needs to be re-tested
208  //real64 const compDens = (ic == 0 && m_compDens[ei][ic] < 1e-6) ? 1e-3 : m_compDens[ei][ic];
209  real64 const dCompAmount_dP = stack.dPoreVolume_dPres * m_compDens[ei][ic];
210  stack.localJacobian[ic][0] += dCompAmount_dP;
211 
212  real64 const dCompAmount_dC = stack.poreVolume;
213  stack.localJacobian[ic][ic + 1] += dCompAmount_dC;
214  }
215  }
216  else
217  {
218  using Deriv = constitutive::multifluid::DerivativeOffset;
219 
220  // construct the slices for variables accessed multiple times
221  arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > dCompFrac_dCompDens = m_dCompFrac_dCompDens[ei];
222 
223  arraySlice1d< real64 const, compflow::USD_PHASE - 1 > phaseVolFrac = m_phaseVolFrac[ei];
224  arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > dPhaseVolFrac = m_dPhaseVolFrac[ei];
225 
226  arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > phaseDens = m_phaseDens[ei][0];
227  arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_DC - 2 > dPhaseDens = m_dPhaseDens[ei][0];
228 
229  arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_COMP - 2 > phaseCompFrac = m_phaseCompFrac[ei][0];
230  arraySlice3d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC - 2 > dPhaseCompFrac = m_dPhaseCompFrac[ei][0];
231 
232  // temporary work arrays
233  real64 dPhaseAmount_dC[numComp]{};
234  real64 dPhaseCompFrac_dC[numComp]{};
235 
236  // start with old time step values
237  for( integer ic = 0; ic < numComp; ++ic )
238  {
239  stack.localResidual[ic] = -m_compAmount_n[ei][ic];
240  }
241 
242  // sum contributions to component accumulation from each phase
243  for( integer ip = 0; ip < m_numPhases; ++ip )
244  {
245  real64 const phaseAmount = stack.poreVolume * phaseVolFrac[ip] * phaseDens[ip];
246 
247  real64 const dPhaseAmount_dP = stack.dPoreVolume_dPres * phaseVolFrac[ip] * phaseDens[ip]
248  + stack.poreVolume * ( dPhaseVolFrac[ip][Deriv::dP] * phaseDens[ip]
249  + phaseVolFrac[ip] * dPhaseDens[ip][Deriv::dP] );
250 
251  // assemble density dependence
252  applyChainRule( numComp, dCompFrac_dCompDens, dPhaseDens[ip], dPhaseAmount_dC, Deriv::dC );
253  for( integer jc = 0; jc < numComp; ++jc )
254  {
255  dPhaseAmount_dC[jc] = dPhaseAmount_dC[jc] * phaseVolFrac[ip]
256  + phaseDens[ip] * dPhaseVolFrac[ip][Deriv::dC + jc];
257  dPhaseAmount_dC[jc] *= stack.poreVolume;
258  }
259 
260  // ic - index of component whose conservation equation is assembled
261  // (i.e. row number in local matrix)
262  for( integer ic = 0; ic < numComp; ++ic )
263  {
264  real64 const phaseCompAmount = phaseAmount * phaseCompFrac[ip][ic];
265 
266  real64 const dPhaseCompAmount_dP = dPhaseAmount_dP * phaseCompFrac[ip][ic]
267  + phaseAmount * dPhaseCompFrac[ip][ic][Deriv::dP];
268 
269  stack.localResidual[ic] += phaseCompAmount;
270  stack.localJacobian[ic][0] += dPhaseCompAmount_dP;
271 
272  // jc - index of component w.r.t. whose compositional var the derivative is being taken
273  // (i.e. col number in local matrix)
274 
275  // assemble phase composition dependence
276  applyChainRule( numComp, dCompFrac_dCompDens, dPhaseCompFrac[ip][ic], dPhaseCompFrac_dC, Deriv::dC );
277  for( integer jc = 0; jc < numComp; ++jc )
278  {
279  real64 const dPhaseCompAmount_dC = dPhaseCompFrac_dC[jc] * phaseAmount
280  + phaseCompFrac[ip][ic] * dPhaseAmount_dC[jc];
281 
282  stack.localJacobian[ic][jc + 1] += dPhaseCompAmount_dC;
283  }
284  }
285 
286  // call the lambda in the phase loop to allow the reuse of the phase amounts and their derivatives
287  // possible use: assemble the derivatives wrt temperature, and the accumulation term of the energy equation for this phase
288  phaseAmountKernelOp( ip, phaseAmount, dPhaseAmount_dP, dPhaseAmount_dC );
289 
290  }
291  }
292  }
293 
301  template< typename FUNC = NoOpFunc >
304  StackVariables & stack,
305  FUNC && phaseVolFractionSumKernelOp = NoOpFunc{} ) const
306  {
307  using Deriv = constitutive::multifluid::DerivativeOffset;
308 
309  arraySlice1d< real64 const, compflow::USD_PHASE - 1 > phaseVolFrac = m_phaseVolFrac[ei];
310  arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > dPhaseVolFrac = m_dPhaseVolFrac[ei];
311 
312  real64 oneMinusPhaseVolFracSum = 1.0;
313 
314  // sum contributions to component accumulation from each phase
315  for( integer ip = 0; ip < m_numPhases; ++ip )
316  {
317  oneMinusPhaseVolFracSum -= phaseVolFrac[ip];
318  stack.localJacobian[numComp][0] -= dPhaseVolFrac[ip][Deriv::dP];
319 
320  for( integer jc = 0; jc < numComp; ++jc )
321  {
322  stack.localJacobian[numComp][jc+1] -= dPhaseVolFrac[ip][Deriv::dC+jc];
323  }
324  }
325 
326  // call the lambda in the phase loop to allow the reuse of the phase amounts and their derivatives
327  // possible use: assemble the derivatives wrt temperature, and use oneMinusPhaseVolFracSum if poreVolume depends on temperature
328  phaseVolFractionSumKernelOp( oneMinusPhaseVolFracSum );
329 
330  // scale saturation-based volume balance by old pore volume (for better scaling w.r.t. other equations)
331  real64 const oldPV = m_volume[ei] * m_porosity_n[ei][0];
332  stack.localResidual[numComp] = oldPV * oneMinusPhaseVolFracSum;
333  for( integer idof = 0; idof < numDof; ++idof )
334  {
335  stack.localJacobian[numComp][idof] *= oldPV;
336  }
337  }
338 
346  StackVariables & stack ) const
347  {
348  using namespace compositionalMultiphaseUtilities;
349 
350  if( m_kernelFlags.isSet( KernelFlags::TotalMassEquation ) )
351  {
352  // apply equation/variable change transformation to the component mass balance equations
353  real64 work[numDof]{};
354  shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numDof, stack.localJacobian, work );
355  shiftElementsAheadByOneAndReplaceFirstElementWithSum( numComp, stack.localResidual );
356  }
357 
358  // add contribution to residual and jacobian into:
359  // - the component mass balance equations (i = 0 to i = numComp-1)
360  // - the volume balance equations (i = numComp)
361  // note that numDof includes derivatives wrt temperature if this class is derived in ThermalKernels
362  integer const numRows = numComp+1;
363  for( integer i = 0; i < numRows; ++i )
364  {
365  m_localRhs[stack.localRow + i] += stack.localResidual[i];
366  m_localMatrix.addToRow< serialAtomic >( stack.localRow + i,
367  stack.dofIndices,
368  stack.localJacobian[i],
369  numDof );
370  }
371  }
372 
380  template< typename POLICY, typename KERNEL_TYPE >
381  static void
382  launch( localIndex const numElems,
383  KERNEL_TYPE const & kernelComponent )
384  {
386 
387  forAll< POLICY >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei )
388  {
389  if( kernelComponent.elemGhostRank( ei ) >= 0 )
390  {
391  return;
392  }
393 
394  typename KERNEL_TYPE::StackVariables stack;
395 
396  kernelComponent.setup( ei, stack );
397  kernelComponent.computeAccumulation( ei, stack );
398  kernelComponent.computeVolumeBalance( ei, stack );
399  kernelComponent.complete( ei, stack );
400  } );
401  }
402 
403 protected:
404 
407 
410 
413 
416 
419 
422  arrayView2d< real64 const > const m_porosity_n;
423  arrayView2d< real64 const > const m_dPoro_dPres;
424 
427 
431 
435 
439 
440  // View on component densities
442 
443  // View on component amount (mass/moles) from previous time step
445 
450 
451  BitFlags< KernelFlags > const m_kernelFlags;
452 };
453 
458 {
459 public:
460 
474  template< typename POLICY >
475  static void
476  createAndLaunch( integer const numComps,
477  integer const numPhases,
478  globalIndex const rankOffset,
479  BitFlags< KernelFlags > kernelFlags,
480  string const dofKey,
481  ElementSubRegionBase const & subRegion,
482  constitutive::MultiFluidBase const & fluid,
483  constitutive::CoupledSolidBase const & solid,
484  CRSMatrixView< real64, globalIndex const > const & localMatrix,
485  arrayView1d< real64 > const & localRhs )
486  {
487  internal::kernelLaunchSelectorCompSwitch( numComps, [&] ( auto NC )
488  {
489  integer constexpr NUM_COMP = NC();
490  integer constexpr NUM_DOF = NC()+1;
491 
492  AccumulationKernel< NUM_COMP, NUM_DOF > kernel( numPhases, rankOffset, dofKey, subRegion,
493  fluid, solid, localMatrix, localRhs, kernelFlags );
494  AccumulationKernel< NUM_COMP, NUM_DOF >::template launch< POLICY >( subRegion.size(), kernel );
495  } );
496  }
497 
498 };
499 
500 } // namespace isothermalCompositionalMultiphaseBaseKernels
501 
502 } // namespace geos
503 
504 
505 #endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_ACCUMULATIONKERNEL_HPP
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
#define GEOS_UNUSED_PARAM(X)
Mark an unused argument and silence compiler warnings.
Definition: GeosxMacros.hpp:72
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
localIndex size() const
Get the "size" of the group, which determines the number of elements in resizable wrappers.
Definition: Group.hpp:1317
static void createAndLaunch(integer const numComps, integer const numPhases, globalIndex const rankOffset, BitFlags< KernelFlags > kernelFlags, string const dofKey, ElementSubRegionBase const &subRegion, constitutive::MultiFluidBase const &fluid, constitutive::CoupledSolidBase const &solid, 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 accumulation and volume balance.
static constexpr integer numDof
Compute time value for the number of degrees of freedom.
arrayView2d< real64 const, compflow::USD_PHASE > const m_phaseVolFrac
Views on the phase volume fractions.
GEOS_HOST_DEVICE void complete(localIndex const GEOS_UNUSED_PARAM(ei), StackVariables &stack) const
Performs the complete phase for the kernel.
arrayView1d< integer const > const m_elemGhostRank
View on the ghost ranks.
CRSMatrixView< real64, globalIndex const > const m_localMatrix
View on the local CRS matrix.
arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > const m_phaseCompFrac
Views on the phase component fraction.
static void launch(localIndex const numElems, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
static constexpr integer numEqn
Compute time value for the number of equations.
AccumulationKernel(localIndex const numPhases, globalIndex const rankOffset, string const dofKey, ElementSubRegionBase const &subRegion, constitutive::MultiFluidBase const &fluid, constitutive::CoupledSolidBase const &solid, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs, BitFlags< KernelFlags > const kernelFlags)
Constructor.
GEOS_HOST_DEVICE void computeVolumeBalance(localIndex const ei, StackVariables &stack, FUNC &&phaseVolFractionSumKernelOp=NoOpFunc{}) const
Compute the local volume balance contributions to the residual and Jacobian.
arrayView3d< real64 const, compflow::USD_COMP_DC > const m_dCompFrac_dCompDens
Views on the derivatives of comp fractions wrt component density.
arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const m_phaseDens
Views on the phase densities.
arrayView2d< real64 const > const m_porosity
Views on the porosity.
arrayView1d< real64 > const m_localRhs
View on the local RHS.
GEOS_HOST_DEVICE void computeAccumulation(localIndex const ei, StackVariables &stack, FUNC &&phaseAmountKernelOp=NoOpFunc{}) const
Compute the local accumulation contributions to the residual and Jacobian.
GEOS_HOST_DEVICE integer elemGhostRank(localIndex const ei) const
Getter for the ghost rank of an element.
GEOS_HOST_DEVICE void setup(localIndex const ei, StackVariables &stack) const
Performs the setup phase for the kernel.
arrayView1d< real64 const > const m_volume
View on the element volumes.
static constexpr integer numComp
Compile time value for the number of components.
arrayView1d< globalIndex const > const m_dofNumber
View on the dof numbers.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:179
LvArray::CRSMatrixView< T, COL_INDEX, localIndex const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
Definition: DataTypes.hpp:309
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
Definition: DataTypes.hpp:87
ArraySlice< T, 2, USD > arraySlice2d
Alias for 2D array slice.
Definition: DataTypes.hpp:199
ArraySlice< T, 3, USD > arraySlice3d
Alias for 3D array slice.
Definition: DataTypes.hpp:215
ArrayView< T, 5, USD > arrayView5d
Alias for 5D array view.
Definition: DataTypes.hpp:243
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:183
std::int32_t integer
Signed integer type.
Definition: DataTypes.hpp:81
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
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
Kernel variables (dof numbers, jacobian and residual) located on the stack.
globalIndex dofIndices[numDof]
Indices of the matrix rows/columns corresponding to the dofs in this element.
localIndex localRow
Index of the local row corresponding to this element.
real64 localJacobian[numEqn][numDof]
C-array storage for the element local Jacobian matrix (all equations except volume balance,...
real64 localResidual[numEqn]
C-array storage for the element local residual vector (all equations except volume balance)