FluxComputeKernel.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_SINGLEPHASE_REACTIVE_FLUXCOMPUTEKERNEL_HPP
#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_FLUXCOMPUTEKERNEL_HPP
 
#include "constitutive/diffusion/DiffusionFields.hpp"
#include "constitutive/diffusion/DiffusionBase.hpp"
#include "constitutive/solid/porosity/PorosityBase.hpp"
#include "constitutive/solid/porosity/PorosityFields.hpp"
#include "constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp"
#include "constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp"
#include "physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernel.hpp"
#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/KernelLaunchSelectors.hpp"
 
 
namespace geos
{
 
namespace singlePhaseReactiveFVMKernels
{
 
template< integer NUM_SPECIES, integer NUM_EQN, integer NUM_DOF, typename STENCILWRAPPER, typename BASE_FLUID_TYPE >
class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, STENCILWRAPPER >
{
public:
 
  static constexpr integer numSpecies = NUM_SPECIES;
 
  static constexpr integer numFluxSupportPoints = 2;
 
  template< typename VIEWTYPE >
  using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
 
  using AbstractBase = singlePhaseFVMKernels::FluxComputeKernelBase;
  using DofNumberAccessor = AbstractBase::DofNumberAccessor;
  using SinglePhaseFlowAccessors = AbstractBase::SinglePhaseFlowAccessors;
  using SinglePhaseFluidAccessors = AbstractBase::SinglePhaseFluidAccessors;
  using PermeabilityAccessors = AbstractBase::PermeabilityAccessors;
 
  using AbstractBase::m_dt;
  using AbstractBase::m_rankOffset;
  using AbstractBase::m_dofNumber;
  using AbstractBase::m_gravCoef;
  using AbstractBase::m_mob;
  using AbstractBase::m_dens;
  using AbstractBase::m_dDens;
 
  using Base = singlePhaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, STENCILWRAPPER >;
  using Base::numDof;
  using Base::numEqn;
  using Base::maxNumElems;
  using Base::maxNumConns;
  using Base::maxStencilSize;
  using Base::m_stencilWrapper;
  using Base::m_seri;
  using Base::m_sesri;
  using Base::m_sei;
 
  using ReactiveSinglePhaseFlowAccessors =
    StencilAccessors< fields::flow::logPrimarySpeciesConcentration,
                      fields::flow::dMobility_dLogPrimaryConc >;
 
  using ReactiveSinglePhaseFluidAccessors =
    StencilMaterialAccessors< constitutive::reactivefluid::ReactiveSinglePhaseFluid< BASE_FLUID_TYPE >,
                              fields::reactivefluid::primarySpeciesMobileAggregateConcentration,
                              fields::reactivefluid::dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations >;
 
  using DiffusionAccessors =
    StencilMaterialAccessors< constitutive::DiffusionBase,
                              fields::diffusion::diffusivity,
                              fields::diffusion::dDiffusivity_dTemperature >;
 
  using PorosityAccessors =
    StencilMaterialAccessors< constitutive::PorosityBase,
                              fields::porosity::referencePorosity >;
 
  FluxComputeKernel( globalIndex const rankOffset,
                     STENCILWRAPPER const & stencilWrapper,
                     DofNumberAccessor const & dofNumberAccessor,
                     SinglePhaseFlowAccessors const & singlePhaseFlowAccessors,
                     ReactiveSinglePhaseFlowAccessors const & reactiveSinglePhaseFlowAccessors,
                     SinglePhaseFluidAccessors const & singlePhaseFluidAccessors,
                     ReactiveSinglePhaseFluidAccessors const & reactiveSinglePhaseFluidAccessors,
                     PermeabilityAccessors const & permeabilityAccessors,
                     DiffusionAccessors const & diffusionAccessors,
                     PorosityAccessors const & porosityAccessors,
                     integer const & hasDiffusion,
                     arrayView1d< integer const > const & mobilePrimarySpeciesFlags,
                     real64 const & solventMassFraction,
                     real64 const & dt,
                     CRSMatrixView< real64, globalIndex const > const & localMatrix,
                     arrayView1d< real64 > const & localRhs )
    : Base( rankOffset,
            stencilWrapper,
            dofNumberAccessor,
            singlePhaseFlowAccessors,
            singlePhaseFluidAccessors,
            permeabilityAccessors,
            dt,
            localMatrix,
            localRhs ),
    m_logPrimarySpeciesConc( reactiveSinglePhaseFlowAccessors.get( fields::flow::logPrimarySpeciesConcentration {} ) ),
    m_dMob_dLogPrimaryConc( reactiveSinglePhaseFlowAccessors.get( fields::flow::dMobility_dLogPrimaryConc {} ) ),
    m_primarySpeciesMobileAggregateConc( reactiveSinglePhaseFluidAccessors.get( fields::reactivefluid::primarySpeciesMobileAggregateConcentration {} ) ),
    m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc( reactiveSinglePhaseFluidAccessors.get(
                                                            fields::reactivefluid::dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations {} ) ),
    m_diffusivity( diffusionAccessors.get( fields::diffusion::diffusivity {} ) ),
    m_dDiffusivity_dTemp( diffusionAccessors.get( fields::diffusion::dDiffusivity_dTemperature {} ) ),
    m_referencePorosity( porosityAccessors.get( fields::porosity::referencePorosity {} ) ),
    m_hasDiffusion( hasDiffusion ),
    m_mobilePrimarySpeciesFlags( mobilePrimarySpeciesFlags ),
    m_solventMassFraction( solventMassFraction )
  {}
 
  struct StackVariables : public Base::StackVariables
  {
public:
 
    GEOS_HOST_DEVICE
    StackVariables( localIndex const size, localIndex numElems )
      : Base::StackVariables( size, numElems )
    {}
 
    using Base::StackVariables::stencilSize;
    using Base::StackVariables::numFluxElems;
    using Base::StackVariables::transmissibility;
    using Base::StackVariables::dTrans_dPres;
    using Base::StackVariables::dofColIndices;
    using Base::StackVariables::localFlux;
    using Base::StackVariables::localFluxJacobian;
 
    real64 diffusionTransmissibility[maxNumConns][numFluxSupportPoints]{};
    real64 dDiffusionTrans_dT[maxNumConns][numFluxSupportPoints]{};
 
  };
 
  template< typename FUNC = NoOpFunc >
  GEOS_HOST_DEVICE
  void computeFlux( localIndex const iconn,
                    StackVariables & stack,
                    FUNC && kernelOp = NoOpFunc{} ) const
  {
    // ***********************************************
    // First, we call the base computeFlux to compute:
    //  1) massFlux and its derivatives,
    //  2) speciesFlux and its derivatives
    Base::computeFlux( iconn, stack, [&] ( localIndex const (&k)[2],
                                           localIndex const (&seri)[2],
                                           localIndex const (&sesri)[2],
                                           localIndex const (&sei)[2],
                                           localIndex const connectionIndex,
                                           real64 const alpha,
                                           real64 const mobility,
                                           real64 const & potGrad,
                                           real64 const & fluxVal,
                                           real64 const (&dFlux_dP)[2] )
    {
      // Step 1: compute the derivatives of the fluid density, potential difference,
      // and the massFlux wrt log of primary species concentration (to complete)
      real64 dFlux_dLogConc[numFluxSupportPoints][numSpecies]{};
 
      GEOS_UNUSED_VAR( dFlux_dLogConc ); // Todo: to add the massFlux derivatives wrt speciesConc
 
      // Step 2: compute the speciesFlux
      real64 speciesFlux[numSpecies]{};
      real64 dSpeciesFlux_dP[numFluxSupportPoints][numSpecies]{};
      real64 dSpeciesFlux_dLogConc[numFluxSupportPoints][numSpecies][numSpecies]{};
      // real64 dSpeciesFlux_dTrans[numSpecies]{};
 
      // choose upstream cell
      localIndex const k_up = (potGrad >= 0) ? 0 : 1;
 
      localIndex const er_up  = seri[k_up];
      localIndex const esr_up = sesri[k_up];
      localIndex const ei_up  = sei[k_up];
 
      // compute species fluxes and derivatives using upstream cell concentration: the mass flux carries
      // molality * solvent mass fraction moles per kg of solution, so the density does not enter
      for( integer is = 0; is < numSpecies; ++is )
      {
        real64 const aggregateConcPerMass_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is] * m_solventMassFraction;
        speciesFlux[is] = aggregateConcPerMass_i * fluxVal;
 
        for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
        {
          dSpeciesFlux_dP[ke][is] += aggregateConcPerMass_i * dFlux_dP[ke];
        }
 
        for( integer js = 0; js < numSpecies; ++js )
        {
          real64 const dAggregateConcPerMass_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er_up][esr_up][ei_up][0][is][js] * m_solventMassFraction;
          dSpeciesFlux_dLogConc[k_up][is][js] += dAggregateConcPerMass_i_dLogConc_j * fluxVal;
        }
      }
 
      for( integer is = 0; is < numSpecies; ++is )
      {
        integer const eqIndex0 = k[0] * numEqn + numEqn - numSpecies + is;
        integer const eqIndex1 = k[1] * numEqn + numEqn - numSpecies + is;
 
        stack.localFlux[eqIndex0] +=  m_dt * speciesFlux[is] * m_mobilePrimarySpeciesFlags[is];
        stack.localFlux[eqIndex1] -=  m_dt * speciesFlux[is] * m_mobilePrimarySpeciesFlags[is];
 
        for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
        {
          localIndex const localDofIndexPres = k[ke] * numDof;
          stack.localFluxJacobian[eqIndex0][localDofIndexPres] += m_dt * dSpeciesFlux_dP[ke][is] * m_mobilePrimarySpeciesFlags[is];
          stack.localFluxJacobian[eqIndex1][localDofIndexPres] -= m_dt * dSpeciesFlux_dP[ke][is] * m_mobilePrimarySpeciesFlags[is];
 
          for( integer js = 0; js < numSpecies; ++js )
          {
            localIndex const localDofIndexSpecies = localDofIndexPres + js + numDof - numSpecies;
            stack.localFluxJacobian[eqIndex0][localDofIndexSpecies] += m_dt * dSpeciesFlux_dLogConc[ke][is][js] * m_mobilePrimarySpeciesFlags[is];
            stack.localFluxJacobian[eqIndex1][localDofIndexSpecies] -= m_dt * dSpeciesFlux_dLogConc[ke][is][js] * m_mobilePrimarySpeciesFlags[is];
          }
        }
      }
 
      // Customize the kernel with this lambda
      kernelOp( k, seri, sesri, sei, connectionIndex, alpha, mobility, potGrad, fluxVal, dFlux_dP );
    } );
  }
 
  template< typename FUNC = NoOpFunc >
  GEOS_HOST_DEVICE
  void computeDiffusion( localIndex const iconn,
                         StackVariables & stack,
                         FUNC && kernelOp = NoOpFunc{} ) const
  {
    using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >;
    if( m_hasDiffusion )
    {
      // *****************************************************
      // Computation of the diffusion term in the species flux
 
      // Step 1: compute the diffusion transmissibilities at this face
      m_stencilWrapper.computeWeights( iconn,
                                       m_diffusivity,
                                       m_dDiffusivity_dTemp,
                                       stack.diffusionTransmissibility,
                                       stack.dDiffusionTrans_dT );
 
      localIndex k[numFluxSupportPoints]{};
      localIndex connectionIndex = 0;
 
      for( k[0] = 0; k[0] < stack.numFluxElems; ++k[0] )
      {
        for( k[1] = k[0] + 1; k[1] < stack.numFluxElems; ++k[1] )
        {
          localIndex const seri[numFluxSupportPoints]  = {m_seri( iconn, k[0] ), m_seri( iconn, k[1] )};
          localIndex const sesri[numFluxSupportPoints] = {m_sesri( iconn, k[0] ), m_sesri( iconn, k[1] )};
          localIndex const sei[numFluxSupportPoints]   = {m_sei( iconn, k[0] ), m_sei( iconn, k[1] )};
 
          // clear working arrays
          real64 diffusionFlux[numSpecies]{};
          real64 speciesGrad[numSpecies]{};
          real64 dDiffusionFlux_dP[numFluxSupportPoints][numSpecies]{};
          real64 dDiffusionFlux_dLogConc[numFluxSupportPoints][numSpecies][numSpecies]{};
 
          real64 const diffusionTrans[numFluxSupportPoints] = { stack.diffusionTransmissibility[connectionIndex][0],
                                                                stack.diffusionTransmissibility[connectionIndex][1] };
 
          //***** calculation of flux *****
          // loop over primary species
          for( integer is = 0; is < numSpecies; ++is )
          {
            real64 dSpeciesGrad_i_dP[numFluxSupportPoints]{};
            real64 dSpeciesGrad_i_dLogConc[numFluxSupportPoints][numSpecies]{};
 
            // Step 2: compute the gradient of the molarity, molality * solvent mass fraction * density
            for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
            {
              localIndex const er  = seri[ke];
              localIndex const esr = sesri[ke];
              localIndex const ei  = sei[ke];
 
              real64 const aggregateConcPerMass_i = m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is] * m_solventMassFraction;
              real64 const dens = m_dens[er][esr][ei][0];
 
              speciesGrad[is] += diffusionTrans[ke] * aggregateConcPerMass_i * dens;
              dSpeciesGrad_i_dP[ke] += diffusionTrans[ke] * aggregateConcPerMass_i * m_dDens[er][esr][ei][0][DerivOffset::dP];
 
              for( integer js = 0; js < numSpecies; ++js )
              {
                real64 const dAggregateConcPerMass_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er][esr][ei][0][is][js] * m_solventMassFraction;
 
                dSpeciesGrad_i_dLogConc[ke][js] += diffusionTrans[ke] * dAggregateConcPerMass_i_dLogConc_j * dens;
              }
            }
 
            // choose upstream cell for species upwinding
            localIndex const k_up = (speciesGrad[is] >= 0) ? 0 : 1;
 
            localIndex const er_up  = seri[k_up];
            localIndex const esr_up = sesri[k_up];
            localIndex const ei_up  = sei[k_up];
 
            // computation of the upwinded species flux
            diffusionFlux[is] += m_referencePorosity[er_up][esr_up][ei_up] * speciesGrad[is];
 
            // add contributions of the derivatives of component fractions wrt pressure/component fractions
            for( integer ke = 0; ke < numFluxSupportPoints; ke++ )
            {
              dDiffusionFlux_dP[ke][is] += m_referencePorosity[er_up][esr_up][ei_up] * dSpeciesGrad_i_dP[ke];
              for( integer js = 0; js < numSpecies; ++js )
              {
                dDiffusionFlux_dLogConc[ke][is][js] += m_referencePorosity[er_up][esr_up][ei_up] * dSpeciesGrad_i_dLogConc[ke][js];
              }
            }
 
            // Add the local diffusion flux contribution to the residual and Jacobian
            integer const eqIndex0 = k[0] * numEqn + numEqn - numSpecies + is;
            integer const eqIndex1 = k[1] * numEqn + numEqn - numSpecies + is;
 
            stack.localFlux[eqIndex0] += m_dt * diffusionFlux[is] * m_mobilePrimarySpeciesFlags[is];
            stack.localFlux[eqIndex1] -= m_dt * diffusionFlux[is] * m_mobilePrimarySpeciesFlags[is];
 
            for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
            {
              localIndex const localDofIndexPres = k[ke] * numDof;
              stack.localFluxJacobian[eqIndex0][localDofIndexPres] += m_dt * dDiffusionFlux_dP[ke][is] * m_mobilePrimarySpeciesFlags[is];
              stack.localFluxJacobian[eqIndex1][localDofIndexPres] -= m_dt * dDiffusionFlux_dP[ke][is] * m_mobilePrimarySpeciesFlags[is];
 
              for( integer js = 0; js < numSpecies; ++js )
              {
                localIndex const localDofIndexSpecies = localDofIndexPres + js + numDof - numSpecies;
                stack.localFluxJacobian[eqIndex0][localDofIndexSpecies] += m_dt * dDiffusionFlux_dLogConc[ke][is][js] * m_mobilePrimarySpeciesFlags[is];
                stack.localFluxJacobian[eqIndex1][localDofIndexSpecies] -= m_dt * dDiffusionFlux_dLogConc[ke][is][js] * m_mobilePrimarySpeciesFlags[is];
              }
            }
 
            // Customize the kernel with this lambda
            kernelOp( is, k, seri, sesri, sei, connectionIndex, k_up );
          }   // loop over primary species
          connectionIndex++;
        }
      }
    }
  }
 
  template< typename FUNC = NoOpFunc >
  GEOS_HOST_DEVICE
  void complete( localIndex const iconn,
                 StackVariables & stack,
                 FUNC && kernelOp = NoOpFunc{} ) const
  {
    // Call Base::complete to assemble the total mass balance equation
    // In the lambda, add contribution to residual and jacobian into the species amount balance equation
    Base::complete( iconn, stack, [&] ( integer const i,
                                        localIndex const localRow )
    {
      // The no. of fluxes is equal to the no. of equations in m_localRhs and m_localMatrix
      for( integer is = 0; is < numSpecies; ++is )
      {
        RAJA::atomicAdd( parallelDeviceAtomic{}, &AbstractBase::m_localRhs[localRow + numEqn - numSpecies + is],
                         stack.localFlux[i * numEqn + numEqn - numSpecies + is] );
        AbstractBase::m_localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >
          ( localRow + numEqn - numSpecies + is,
          stack.dofColIndices.data(),
          stack.localFluxJacobian[i * numEqn + numEqn - numSpecies + is].dataIfContiguous(),
          stack.stencilSize * numDof );
      }
 
      // call the lambda to assemble additional terms, such as thermal terms
      kernelOp( i, localRow );
    } );
  }
 
  template< typename POLICY, typename KERNEL_TYPE >
  static void
  launch( localIndex const numConnections,
          KERNEL_TYPE const & kernelComponent )
  {
    GEOS_MARK_FUNCTION;
 
    forAll< POLICY >( numConnections, [=] GEOS_HOST_DEVICE ( localIndex const iconn )
    {
      typename KERNEL_TYPE::StackVariables stack( kernelComponent.stencilSize( iconn ),
                                                  kernelComponent.numPointsInFlux( iconn ) );
 
      kernelComponent.setup( iconn, stack );
      kernelComponent.computeFlux( iconn, stack );
      kernelComponent.computeDiffusion( iconn, stack );
      kernelComponent.complete( iconn, stack );
    } );
  }
 
protected:
 
  ElementViewConst< arrayView2d< real64 const, compflow::USD_COMP > > const m_logPrimarySpeciesConc;
 
  ElementViewConst< arrayView2d< real64 const, compflow::USD_FLUID_DC > > const m_dMob_dLogPrimaryConc;
 
  ElementViewConst< arrayView3d< real64 const, constitutive::reactivefluid::USD_SPECIES > > const m_primarySpeciesMobileAggregateConc;
 
  ElementViewConst< arrayView4d< real64 const, constitutive::reactivefluid::USD_SPECIES_DC > > const m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc;
 
  ElementViewConst< arrayView3d< real64 const > > const m_diffusivity;
  ElementViewConst< arrayView3d< real64 const > > const m_dDiffusivity_dTemp;
 
  ElementViewConst< arrayView1d< real64 const > > const m_referencePorosity;
 
  integer const m_hasDiffusion;
 
  arrayView1d< integer const > const m_mobilePrimarySpeciesFlags;
 
  real64 const m_solventMassFraction;
};
 
class FluxComputeKernelFactory
{
public:
 
  template< typename POLICY, typename STENCILWRAPPER >
  static void
  createAndLaunch( integer const numSpecies,
                   integer const hasDiffusion,
                   arrayView1d< integer const > const mobilePrimarySpeciesFlags,
                   real64 const solventMassFraction,
                   globalIndex const rankOffset,
                   string const & dofKey,
                   string const & solverName,
                   ElementRegionManager const & elemManager,
                   STENCILWRAPPER const & stencilWrapper,
                   real64 const & dt,
                   CRSMatrixView< real64, globalIndex const > const & localMatrix,
                   arrayView1d< real64 > const & localRhs )
  {
    singlePhaseReactiveBaseKernels::internal::kernelLaunchSelectorCompSwitch( numSpecies, [&]( auto NS )
    {
      integer constexpr NUM_SPECIES = NS();
      integer constexpr NUM_DOF = 1+NS();
      integer constexpr NUM_EQN = 1+NS();
 
      ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > > dofNumberAccessor =
        elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey );
      dofNumberAccessor.setName( solverName + "/accessors/" + dofKey );
 
      using KernelType = FluxComputeKernel< NUM_SPECIES, NUM_EQN, NUM_DOF, STENCILWRAPPER, constitutive::CompressibleSinglePhaseFluid >;
      typename KernelType::SinglePhaseFlowAccessors flowAccessors( elemManager, solverName );
      typename KernelType::ReactiveSinglePhaseFlowAccessors reactiveFlowAccessors( elemManager, solverName );
      typename KernelType::SinglePhaseFluidAccessors fluidAccessors( elemManager, solverName );
      typename KernelType::ReactiveSinglePhaseFluidAccessors reactiveFluidAccessors( elemManager, solverName );
      typename KernelType::PermeabilityAccessors permAccessors( elemManager, solverName );
      typename KernelType::DiffusionAccessors diffusionAccessors( elemManager, solverName );
      typename KernelType::PorosityAccessors porosityAccessors( elemManager, solverName );
 
      KernelType kernel( rankOffset, stencilWrapper, dofNumberAccessor,
                         flowAccessors, reactiveFlowAccessors, fluidAccessors, reactiveFluidAccessors,
                         permAccessors, diffusionAccessors, porosityAccessors, hasDiffusion, mobilePrimarySpeciesFlags,
                         solventMassFraction, dt, localMatrix, localRhs );
      KernelType::template launch< POLICY >( stencilWrapper.size(), kernel );
    } );
  }
};
 
} // namespace singlePhaseReactiveFVMKernels
 
} // namespace geos
 
#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_FLUXCOMPUTEKERNEL_HPP