SourceFluxComputeKernel.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 TotalEnergies
 * 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_SINGLEPHASEREACTIVE_SOURCEFLUXCOMPUTEKERNELS_HPP
#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_SOURCEFLUXCOMPUTEKERNELS_HPP
 
#include "common/DataTypes.hpp"
#include "common/GEOS_RAJA_Interface.hpp"
#include "constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp"
#include "constitutive/fluid/reactivefluid/ReactiveFluidLayouts.hpp"
#include "constitutive/fluid/singlefluid/SingleFluidBase.hpp"
#include "constitutive/fluid/singlefluid/SingleFluidUtils.hpp"
#include "codingUtilities/Utilities.hpp"
 
#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/KernelLaunchSelectors.hpp"
 
namespace geos
{
 
namespace singlePhaseReactiveBaseKernels
{
 
/******************************** SourceFluxComputeKernel ********************************/
 
template< integer NUM_DOF, integer NUM_SPECIES, typename BASE_FLUID_TYPE >
class SourceFluxComputeKernel
{
 
public:
 
  static constexpr integer numDof = NUM_DOF;
 
  static constexpr integer numEqn = NUM_DOF;
 
  static constexpr integer numSpecies = NUM_SPECIES;
 
  using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >;
 
  SourceFluxComputeKernel( globalIndex const rankOffset,
                           arrayView1d< globalIndex const > const dofNumber,
                           arrayView1d< integer const > const elemGhostRank,
                           arrayView1d< real64 const > const rhsContributionArrayView,
                           real64 const sizeScalingFactor,
                           constitutive::reactivefluid::ReactiveSinglePhaseFluid< BASE_FLUID_TYPE > const & fluid,
                           CRSMatrixView< real64, globalIndex const > const & localMatrix,
                           arrayView1d< real64 > const & localRhs,
                           RAJA::ReduceSum< parallelDeviceReduce, real64 > massProd )
    :
    m_rankOffset( rankOffset ),
    m_dofNumber( dofNumber ),
    m_elemGhostRank( elemGhostRank ),
    m_rhsContributionArrayView( rhsContributionArrayView ),
    m_sizeScalingFactor( sizeScalingFactor ),
    m_primarySpeciesAggregateConcentration( fluid.primarySpeciesAggregateConcentration() ),
    m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations( fluid.dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations() ),
    m_density( fluid.density() ),
    m_dDensity( fluid.dDensity() ),
    m_localMatrix( localMatrix ),
    m_localRhs( localRhs ),
    m_massProd( massProd )
  {}
 
  struct StackVariables
  {
public:
 
    localIndex massRowIndex = -1;
 
    globalIndex dofIndices[numDof]{};
 
    real64 localSpeciesRhs[numSpecies]{};
 
    real64 localSpeciesJacobian[numSpecies][numDof]{};
 
    real64 totalInflowMass = 0.0;
 
  };
 
  GEOS_HOST_DEVICE
  integer elemGhostRank( localIndex const ei ) const
  { return m_elemGhostRank( ei ); }
 
  GEOS_HOST_DEVICE
  void setup( localIndex const ei,
              localIndex const a,
              StackVariables & stack ) const
  {
    // set row index and degrees of freedom indices for this element
    stack.massRowIndex = m_dofNumber[ei] - m_rankOffset;
    for( integer idof = 0; idof < numDof; ++idof )
    {
      stack.dofIndices[idof] = m_dofNumber[ei] + idof - m_rankOffset;
    }
 
    stack.totalInflowMass = m_rhsContributionArrayView[a];
  }
 
  GEOS_HOST_DEVICE
  void computeSourceFlux( localIndex const ei,
                          StackVariables & stack ) const
  {
    real64 const scaledInflowMass = stack.totalInflowMass / m_sizeScalingFactor;
 
    for( integer i = 0; i < numSpecies; ++i )
    {
      stack.localSpeciesRhs[i] += m_primarySpeciesAggregateConcentration[ei][0][i] / m_density[ei][0] * scaledInflowMass;
      stack.localSpeciesJacobian[i][0] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_dDensity[ei][0][DerivOffset::dP] / (m_density[ei][0] * m_density[ei][0]) * scaledInflowMass;
 
      for( integer j = 0; j < numSpecies; ++j )
      {
        stack.localSpeciesJacobian[i][j+numDof-numSpecies] += m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[ei][0][i][j] / m_density[ei][0] * scaledInflowMass;
      }
    }
  }
 
  GEOS_HOST_DEVICE
  void complete( localIndex const GEOS_UNUSED_PARAM( ei ),
                 StackVariables & stack ) const
  {
    m_massProd += stack.totalInflowMass / m_sizeScalingFactor;
    m_localRhs[stack.massRowIndex] += stack.totalInflowMass / m_sizeScalingFactor;
 
    if( stack.totalInflowMass > 0.0 )
    {
      for( integer i = 0; i < numSpecies; ++i )
      {
        globalIndex const speciesRowBeginIndex = stack.massRowIndex + numEqn - numSpecies;
        m_localRhs[speciesRowBeginIndex + i] += stack.localSpeciesRhs[i];
 
        // add contribution to global residual and jacobian (no need for atomics here)
        m_localMatrix.template addToRow< serialAtomic >( speciesRowBeginIndex+i,
                                                         stack.dofIndices,
                                                         stack.localSpeciesJacobian[i],
                                                         numDof );
      }
    }
  }
 
  template< typename POLICY, typename KERNEL_TYPE >
  static void
  launch( SortedArrayView< localIndex const > const targetSet,
          KERNEL_TYPE const & kernelComponent )
  {
    GEOS_MARK_FUNCTION;
 
    forAll< POLICY >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
    {
      // we need to filter out ghosts here, because targetSet may contain them
      localIndex const ei = targetSet[a];
 
      if( kernelComponent.elemGhostRank( ei ) >= 0 )
      {
        return;
      }
 
      typename KERNEL_TYPE::StackVariables stack;
 
      kernelComponent.setup( ei, a, stack );
      kernelComponent.computeSourceFlux( ei, stack );
      kernelComponent.complete( ei, stack );
    } );
  }
 
protected:
 
  globalIndex const m_rankOffset;
 
  arrayView1d< globalIndex const > const m_dofNumber;
 
  arrayView1d< integer const > const m_elemGhostRank;
 
  arrayView1d< real64 const > const m_rhsContributionArrayView;
  real64 const m_sizeScalingFactor;
 
  // View on the total concentration of ions that contain the primary species
  arrayView3d< real64 const, constitutive::reactivefluid::USD_SPECIES > const m_primarySpeciesAggregateConcentration;
  // View on the derivatives of total ion concentration for the primary species wrt log of primary species concentration
  arrayView4d< real64 const, constitutive::reactivefluid::USD_SPECIES_DC > const m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations;
 
  // View on the fluid density
  arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_density;
  // View on the derivatives of fluid density
  arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > const m_dDensity;
 
  CRSMatrixView< real64, globalIndex const > const m_localMatrix;
  arrayView1d< real64 > const m_localRhs;
 
  RAJA::ReduceSum< parallelDeviceReduce, real64 > m_massProd;
 
};
 
class SourceFluxComputeKernelFactory
{
public:
 
  template< typename POLICY, typename BASE_FLUID_TYPE >
  static void
  createAndLaunch( integer const numSpecies,
                   globalIndex const rankOffset,
                   arrayView1d< globalIndex const > const dofNumber,
                   arrayView1d< integer const > const elemGhostRank,
                   SortedArrayView< localIndex const > const targetSet,
                   arrayView1d< real64 const > const rhsContributionArrayView,
                   real64 const sizeScalingFactor,
                   constitutive::reactivefluid::ReactiveSinglePhaseFluid< BASE_FLUID_TYPE > const & fluid,
                   CRSMatrixView< real64, globalIndex const > const & localMatrix,
                   arrayView1d< real64 > const & localRhs,
                   RAJA::ReduceSum< parallelDeviceReduce, real64 > massProd )
  {
    internal::kernelLaunchSelectorCompSwitch( numSpecies, [&] ( auto NS )
    {
      integer constexpr NUM_SPECIES = NS();
      integer constexpr NUM_DOF = 1+NS();
 
      SourceFluxComputeKernel< NUM_DOF, NUM_SPECIES, BASE_FLUID_TYPE > kernel( rankOffset, dofNumber, elemGhostRank, rhsContributionArrayView, sizeScalingFactor, fluid, localMatrix, localRhs,
                                                                               massProd );
      SourceFluxComputeKernel< NUM_DOF, NUM_SPECIES, BASE_FLUID_TYPE >::template launch< POLICY >( targetSet, kernel );
    } );
  }
};
} // namespace singlePhaseReactiveBaseKernels
 
} // namespace geos
 
#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_SOURCEFLUXCOMPUTEKERNELS_HPP