AccumulationZFormulationKernel.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 Total, S.A
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_ACCUMULATIONZFORMULATIONKERNEL_HPP
#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_ACCUMULATIONZFORMULATIONKERNEL_HPP
 
#include "codingUtilities/Utilities.hpp"
#include "common/DataLayouts.hpp"
#include "common/DataTypes.hpp"
#include "common/GEOS_RAJA_Interface.hpp"
#include "constitutive/solid/CoupledSolidBase.hpp"
#include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
#include "mesh/ElementSubRegionBase.hpp"
#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp"
#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp"
#include "physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp"
#include "physicsSolvers/fluidFlow/kernels/compositional/KernelLaunchSelectors.hpp"
 
namespace geos
{
 
namespace isothermalCompositionalMultiphaseBaseKernels
{
 
static constexpr real64 minCompFracForDivision = 0;
 
/******************************** AccumulationKernel ********************************/
 
template< integer NUM_COMP, integer NUM_DOF >
class AccumulationZFormulationKernel
{
public:
 
  static constexpr integer numComp = NUM_COMP;
 
  static constexpr integer numDof = NUM_DOF;
 
  static constexpr integer numEqn = NUM_DOF;
 
  AccumulationZFormulationKernel( 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 )
    : m_numPhases( numPhases ),
    m_rankOffset( rankOffset ),
    m_dofNumber( subRegion.getReference< array1d< globalIndex > >( dofKey ) ),
    m_elemGhostRank( subRegion.ghostRank() ),
    m_volume( subRegion.getElementVolume() ),
    m_porosity( solid.getPorosity() ),
    m_dPoro_dPres( solid.getDporosity_dPressure() ),
    m_phaseVolFrac( subRegion.getField< fields::flow::phaseVolumeFraction >() ),
    m_dPhaseVolFrac( subRegion.getField< fields::flow::dPhaseVolumeFraction >() ),
    m_totalDens( fluid.totalDensity() ),
    m_dTotalDens( fluid.dTotalDensity() ),
    m_compFrac( subRegion.getField< fields::flow::globalCompFraction >() ),
    m_compAmount_n( subRegion.getField< fields::flow::compAmount_n >() ),
    m_localMatrix( localMatrix ),
    m_localRhs( localRhs ),
    m_KernelFlags( KernelFlags )
  {}
 
  struct StackVariables
  {
public:
 
    // Pore volume information (used by both accumulation and volume balance)
 
    real64 poreVolume = 0.0;
 
    real64 dPoreVolume_dPres = 0.0;
 
    // Residual information
 
    localIndex localRow = -1;
 
    globalIndex dofIndices[numDof]{};
 
    real64 localResidual[numEqn]{};
 
    real64 localJacobian[numEqn][numDof]{};
 
  };
 
  GEOS_HOST_DEVICE
  integer elemGhostRank( localIndex const ei ) const
  { return m_elemGhostRank( ei ); }
 
 
  GEOS_HOST_DEVICE
  void setup( localIndex const ei,
              StackVariables & stack ) const
  {
    // initialize the pore volume
    stack.poreVolume = m_volume[ei] * m_porosity[ei][0];
    stack.dPoreVolume_dPres = m_volume[ei] * m_dPoro_dPres[ei][0];
 
    // set row index and degrees of freedom indices for this element
    stack.localRow = m_dofNumber[ei] - m_rankOffset;
    for( integer idof = 0; idof < numDof; ++idof )
    {
      stack.dofIndices[idof] = m_dofNumber[ei] + idof;
    }
  }
 
  template< typename FUNC = NoOpFunc >
  GEOS_HOST_DEVICE
  void computeAccumulation( localIndex const ei,
                            StackVariables & stack ) const
  {
    using Deriv = constitutive::multifluid::DerivativeOffset;
 
    arraySlice1d< real64 const, compflow::USD_COMP - 1 > const compFrac = m_compFrac[ei];
    real64 const totalDensity = m_totalDens[ei][0];
    arraySlice1d< real64 const, constitutive::multifluid::USD_FLUID_DC - 2 > const dTotalDens = m_dTotalDens[ei][0];
 
    // ic - index of component whose conservation equation is assembled
    // (i.e. row number in local matrix)
    for( integer ic = 0; ic < numComp; ++ic )
    {
      real64 const compAmount = stack.poreVolume * totalDensity * compFrac[ic];
      real64 const compAmount_n = m_compAmount_n[ei][ic];
 
      stack.localResidual[ic] += compAmount - compAmount_n;
 
      // derivatives with respect to pressure (p)
      real64 const dCompAmount_dP = compFrac[ic] * (stack.dPoreVolume_dPres * totalDensity + stack.poreVolume * dTotalDens[Deriv::dP]);
      stack.localJacobian[ic][0] += dCompAmount_dP;
 
      // derivatives with respect to global component fraction (zc)
      for( integer jc = 0; jc < numComp; ++jc )
      {
        real64 dCompAmount_dC;
        if( ic == jc )
          dCompAmount_dC = stack.poreVolume * (totalDensity + dTotalDens[Deriv::dC+jc] * compFrac[ic]);
        else
          dCompAmount_dC = stack.poreVolume * (dTotalDens[Deriv::dC+jc] * compFrac[ic]);
 
        stack.localJacobian[ic][jc + 1] += dCompAmount_dC;
      }
    }
  }
 
  template< typename FUNC = NoOpFunc >
  GEOS_HOST_DEVICE
  void computeCompFracSum( localIndex const ei,
                           StackVariables & stack ) const
  {
 
    arraySlice1d< real64 const, compflow::USD_PHASE - 1 > compFrac = m_compFrac[ei];
 
    real64 oneMinusCompFracSum = 1.0;
 
    // sum contributions to component accumulation from each component
    for( integer ic = 0; ic < numComp; ++ic )
    {
      oneMinusCompFracSum -= compFrac[ic];
      // no derivatives w.r.t pressure
 
      // derivative w.r.t component fractions are unity: dzc_dzc = 1
      stack.localJacobian[numComp][ic+1] -= 1;
    }
 
    //phaseVolFractionSumKernelOp( oneMinusCompFracSum );
 
    // scale componentFraction-based volume balance by pore volume (for better scaling w.r.t. other equations)
    stack.localResidual[numComp] = stack.poreVolume * oneMinusCompFracSum;
    for( integer idof = 0; idof < numDof; ++idof )
    {
      stack.localJacobian[numComp][idof] *= stack.poreVolume;
    }
    stack.localJacobian[numComp][0] += stack.dPoreVolume_dPres * oneMinusCompFracSum;
  }
 
  GEOS_HOST_DEVICE
  void complete( localIndex const GEOS_UNUSED_PARAM( ei ),
                 StackVariables & stack ) const
  {
    using namespace compositionalMultiphaseUtilities;
 
    if( m_KernelFlags.isSet( KernelFlags::TotalMassEquation ) )
    {
      // apply equation/variable change transformation to the component mass balance equations
      real64 work[numDof]{};
      shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numDof, stack.localJacobian, work );
      shiftElementsAheadByOneAndReplaceFirstElementWithSum( numComp, stack.localResidual );
    }
 
    // add contribution to residual and jacobian into:
    // - the component mass balance equations (i = 0 to i = numComp-1)
    // - the volume balance equations (i = numComp)
    // note that numDof includes derivatives wrt temperature if this class is derived in ThermalKernels
    integer const numRows = numComp+1;
    for( integer i = 0; i < numRows; ++i )
    {
      m_localRhs[stack.localRow + i] += stack.localResidual[i];
      m_localMatrix.addToRow< serialAtomic >( stack.localRow + i,
                                              stack.dofIndices,
                                              stack.localJacobian[i],
                                              numDof );
    }
  }
 
  template< typename POLICY, typename KERNEL_TYPE >
  static void
  launch( localIndex const numElems,
          KERNEL_TYPE const & kernelComponent )
  {
    GEOS_MARK_FUNCTION;
 
    forAll< POLICY >( numElems, [=] GEOS_HOST_DEVICE ( localIndex const ei )
    {
      if( kernelComponent.elemGhostRank( ei ) >= 0 )
      {
        return;
      }
 
      typename KERNEL_TYPE::StackVariables stack;
 
      kernelComponent.setup( ei, stack );
      kernelComponent.computeAccumulation( ei, stack );
      kernelComponent.computeCompFracSum( ei, stack );
      kernelComponent.complete( ei, stack );
    } );
  }
 
protected:
 
  integer const m_numPhases;
 
  globalIndex const m_rankOffset;
 
  arrayView1d< globalIndex const > const m_dofNumber;
 
  arrayView1d< integer const > const m_elemGhostRank;
 
  arrayView1d< real64 const > const m_volume;
 
  arrayView2d< real64 const > const m_porosity;
  arrayView2d< real64 const > const m_dPoro_dPres;
 
  arrayView3d< real64 const, compflow::USD_COMP_DC > const m_dCompFrac_dCompDens;
 
  arrayView2d< real64 const, compflow::USD_PHASE > const m_phaseVolFrac;
  arrayView3d< real64 const, compflow::USD_PHASE_DC > const m_dPhaseVolFrac;
 
  arrayView2d< real64 const, constitutive::multifluid::USD_FLUID > const m_totalDens;
  arrayView3d< real64 const, constitutive::multifluid::USD_FLUID_DC > const m_dTotalDens;
 
  // View on component densities and component fractions
  arrayView2d< real64 const, compflow::USD_COMP > m_compFrac;
 
  // View on component amount (mass/moles) from previous time step
  arrayView2d< real64 const, compflow::USD_COMP > m_compAmount_n;
 
  CRSMatrixView< real64, globalIndex const > const m_localMatrix;
  arrayView1d< real64 > const m_localRhs;
 
  BitFlags< KernelFlags > const m_KernelFlags;
};
 
class AccumulationZFormulationKernelFactory
{
public:
 
  template< typename POLICY >
  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 )
  {
    internal::kernelLaunchSelectorCompSwitch( numComps, [&] ( auto NC )
    {
      integer constexpr NUM_COMP = NC();
      integer constexpr NUM_DOF = NC()+1;
 
      AccumulationZFormulationKernel< NUM_COMP, NUM_DOF > kernel( numPhases, rankOffset, dofKey, subRegion,
                                                                  fluid, solid, localMatrix, localRhs, kernelFlags );
      AccumulationZFormulationKernel< NUM_COMP, NUM_DOF >::template launch< POLICY >( subRegion.size(), kernel );
    } );
  }
 
};
 
} // namespace isothermalCompositionalMultiphaseBaseKernels
 
} // namespace geos
 
 
#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_ACCUMULATIONZFORMULATIONKERNEL_HPP