AccumulationKernels.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_ACCUMULATIONKERNELS_HPP
#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_ACCUMULATIONKERNELS_HPP
 
#include "common/DataLayouts.hpp"
#include "common/DataTypes.hpp"
#include "constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp"
#include "constitutive/fluid/reactivefluid/ReactiveFluidLayouts.hpp"
#include "constitutive/solid/CoupledSolidBase.hpp"
#include "physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp"
#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/KernelLaunchSelectors.hpp"
 
namespace geos
{
 
namespace singlePhaseReactiveBaseKernels
{
 
/******************************** AccumulationKernel ********************************/
 
template< typename SUBREGION_TYPE, integer NUM_DOF, integer NUM_SPECIES, typename BASE_FLUID_TYPE >
class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SUBREGION_TYPE, NUM_DOF >
{
 
public:
 
  using Base = singlePhaseBaseKernels::AccumulationKernel< SUBREGION_TYPE, NUM_DOF >;
  using Base::numDof;
  using Base::numEqn;
  using Base::m_rankOffset;
  using Base::m_dofNumber;
  using Base::m_elemGhostRank;
  using Base::m_localMatrix;
  using Base::m_localRhs;
  using Base::m_dMass;
 
  using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 0 >;
 
  static constexpr integer numSpecies = NUM_SPECIES;
 
  AccumulationKernel( globalIndex const rankOffset,
                      string const dofKey,
                      SUBREGION_TYPE const & subRegion,
                      constitutive::reactivefluid::ReactiveSinglePhaseFluid< BASE_FLUID_TYPE > const & fluid,
                      constitutive::CoupledSolidBase const & solid,
                      real64 const & dt,
                      CRSMatrixView< real64, globalIndex const > const & localMatrix,
                      arrayView1d< real64 > const & localRhs )
    : Base( rankOffset, dofKey, subRegion, localMatrix, localRhs ),
    m_dt( dt ),
    m_solventMassFraction( fluid.solventMassFraction() ),
    m_density( fluid.density() ),
    m_dDensity( fluid.dDensity() ),
    m_volume( subRegion.getElementVolume() ),
    m_deltaVolume( subRegion.template getField< fields::flow::deltaVolume >() ),
    m_porosity( solid.getPorosity() ),
    m_dPoro_dPres( solid.getDporosity_dPressure() ),
    // m_dDensity_dLogPrimaryConc( fluid.dDensity_dLogPrimaryConc() ),
    // m_dPoro_dLogPrimaryConc( solid.getDporosity_dLogPrimaryConc() ),
    m_primarySpeciesAggregateConcentration( fluid.primarySpeciesAggregateConcentration() ),
    // m_dPrimarySpeciesAggregateConcentration_dPres( fluid.dPrimarySpeciesAggregateConcentration_dPres() ),
    m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations( fluid.dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations() ),
    m_primarySpeciesAggregateKineticRate( fluid.aggregateSpeciesRates() ),
    // m_dPrimarySpeciesAggregateKineticRate_dPres( fluid.dPrimarySpeciesAggregateKineticRate_dPres() ),
    m_dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc( fluid.dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations() ),
    m_primarySpeciesAggregateMole_n( subRegion.template getField< fields::flow::primarySpeciesAggregateMole_n >() )
  {}
 
  struct StackVariables : public Base::StackVariables
  {
public:
 
    GEOS_HOST_DEVICE
    StackVariables()
      : Base::StackVariables()
    {}
 
    using Base::StackVariables::localRow;
    using Base::StackVariables::dofIndices;
    using Base::StackVariables::localResidual;
    using Base::StackVariables::localJacobian;
 
    // Pore volume information
 
    real64 poreVolume = 0.0;
 
    real64 dPoreVolume_dPres = 0.0;
 
    real64 dPoreVolume_dLogPrimaryConc[numSpecies]{};
 
  };
 
  GEOS_HOST_DEVICE
  void setup( localIndex const ei,
              StackVariables & stack ) const
  {
    // initialize the pore volume
    stack.poreVolume = ( m_volume[ei] + m_deltaVolume[ei] ) * m_porosity[ei][0];
    stack.dPoreVolume_dPres = ( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dPres[ei][0];
 
    Base::setup( ei, stack );
 
    // // is - index of the primary species
    // for( integer is = 0; is < numSpecies; ++is )
    // {
    //   stack.dPoreVolume_dLogPrimaryConc[is] = ( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dLogPrimaryConc[ei][is]
    // }
  }
 
  GEOS_HOST_DEVICE
  void computeAccumulation( localIndex const ei,
                            StackVariables & stack ) const
  {
    // Residual[is] += (primarySpeciesAggregateConcentration[is] * solventMassFraction * density * stack.poreVolume -
    // primarySpeciesAggregateMole_n[is])
    //                 - dt * m_volume * primarySpeciesKineticRate[is]
 
    Base::computeAccumulation( ei, stack );
 
    // Step 1: assemble the derivatives of the total mass equation w.r.t log of primary species concentration
    //   for( integer is = 0; is < numSpecies; ++is )
    //   {
    //     stack.localJacobian[0][is+numDof-numSpecies] = stack.poreVolume * m_dDensity_dLogPrimaryConc[ei][is] +
    // stack.dPoreVolume_dLogPrimaryConc[is] * m_density[ei][0];
    //   }
 
    arraySlice2d< real64 const,
                  constitutive::reactivefluid::USD_SPECIES_DC -
                  2 > dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations = m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[ei][0];
    arraySlice2d< real64 const, constitutive::reactivefluid::USD_SPECIES_DC - 2 > dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc = m_dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc[ei][0];
 
    for( integer is = 0; is < numSpecies; ++is )
    {
      // Step 2: assemble the accumulation term of the species mass balance equation
      // Step 2.1: residual
      // Primary species mole amount in pore volume
      stack.localResidual[is+numEqn-numSpecies] -= m_primarySpeciesAggregateMole_n[ei][is];
      // moles of species per m^3 of solution: molality * solvent mass fraction * density
      real64 const aggregateConcMolarity = m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventMassFraction * m_density[ei][0];
      stack.localResidual[is+numEqn-numSpecies] += aggregateConcMolarity * stack.poreVolume;
 
      // Reaction term
      stack.localResidual[is+numEqn-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * m_primarySpeciesAggregateKineticRate[ei][0][is];
 
      // Step 2.1: jacobian
      // Derivative of primary species amount in pore volume wrt pressure
      stack.localJacobian[is+numEqn-numSpecies][0] += stack.dPoreVolume_dPres * aggregateConcMolarity
                                                      + stack.poreVolume * m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventMassFraction * m_dDensity[ei][0][DerivOffset::dP]
                                                      /* + stack.poreVolume * m_dTotalPrimarySpeciesConcentration_dPres[ei][is] */;
      // // Derivative of reaction term wrt pressure
      // stack.localJacobian[is+numEqn-numSpecies][0] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) *
      // m_dPrimarySpeciesTotalKineticRate_dPres[is];
 
      // Derivative wrt log of primary species concentration
      for( integer js = 0; js < numSpecies; ++js )
      {
        stack.localJacobian[is+numEqn-numSpecies][js+numDof-numSpecies] = /* stack.dPoreVolume_dLogPrimaryConc[js] *
                                                                              m_primarySpeciesAggregateConcentration[ei][0][is]
                                                                           + */stack.poreVolume * dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[is][js] *
                                                                          m_solventMassFraction * m_density[ei][0];
 
        stack.localJacobian[is+numEqn-numSpecies][js+numDof-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc[is][js];
      }
    }
  }
 
  GEOS_HOST_DEVICE
  void complete( localIndex const ei,
                 StackVariables & stack ) const
  {
    // Step 1: assemble the total mass balance equation
    // - the total mass balance equations (i = 0)
    Base::complete( ei, stack );
 
    // Step 2: assemble the primary species mole amount balance equation
    // - the species mole amount balance equations (i = numEqn-numSpecies to i = numEqn-1)
    integer const beginRowSpecies = numEqn-numSpecies;
    for( integer i = 0; i < numSpecies; ++i )
    {
      m_localRhs[stack.localRow + beginRowSpecies + i] += stack.localResidual[beginRowSpecies+i];
      m_localMatrix.template addToRow< serialAtomic >( stack.localRow + beginRowSpecies + i,
                                                       stack.dofIndices,
                                                       stack.localJacobian[beginRowSpecies + i],
                                                       numDof );
    }
  }
 
protected:
 
  real64 const m_dt;
 
  real64 const m_solventMassFraction;
 
  arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_density;
  arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > const m_dDensity;
 
  arrayView1d< real64 const > const m_volume;
  arrayView1d< real64 const > const m_deltaVolume;
 
  arrayView2d< real64 const > const m_porosity;
  arrayView2d< real64 const > const m_dPoro_dPres;
 
  // // View on the derivatives of fluid density wrt log of primary species concentration
  // arrayView2d< real64 const, compflow::USD_COMP > m_dDensity_dLogPrimaryConc;
 
  // // View on the derivatives of porosity wrt log of primary species concentration
  // arrayView2d< real64 const, compflow::USD_COMP > m_dPoro_dLogPrimaryConc;
 
  // View on the total concentration of ions that contain the primary species
  arrayView3d< real64 const, constitutive::reactivefluid::USD_SPECIES > m_primarySpeciesAggregateConcentration;
 
  // // View on the derivatives of aggregate concentration for the primary species wrt pressure
  // arrayView2d< real64 const, compflow::USD_COMP > m_dPrimarySpeciesAggregateConcentration_dPres;
 
  // 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 > m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations;
 
  // View on the aggregate kinetic rate of primary species from all reactions
  arrayView3d< real64 const, constitutive::reactivefluid::USD_SPECIES > m_primarySpeciesAggregateKineticRate;
 
  // // View on the derivatives of aggregate kinetic rate of primary species wrt pressure
  // arrayView2d< real64 const, compflow::USD_COMP > m_dPrimarySpeciesAggregateKineticRate_dPres;
 
  // View on the derivatives of aggregate kinetic rate of primary species wrt log of primary species concentration
  arrayView4d< real64 const, constitutive::reactivefluid::USD_SPECIES_DC > m_dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc;
 
  // View on primary species mole amount from previous time step
  arrayView2d< real64 const, compflow::USD_COMP > m_primarySpeciesAggregateMole_n;
};
 
class AccumulationKernelFactory
{
public:
 
  template< typename POLICY, typename SUBREGION_TYPE, typename BASE_FLUID_TYPE >
  static void
  createAndLaunch( integer const numSpecies,
                   real64 const dt,
                   globalIndex const rankOffset,
                   string const dofKey,
                   SUBREGION_TYPE const & subRegion,
                   constitutive::reactivefluid::ReactiveSinglePhaseFluid< BASE_FLUID_TYPE > const & fluid,
                   constitutive::CoupledSolidBase const & solid,
                   CRSMatrixView< real64, globalIndex const > const & localMatrix,
                   arrayView1d< real64 > const & localRhs )
  {
    internal::kernelLaunchSelectorCompSwitch( numSpecies, [&] ( auto NS )
    {
      integer constexpr NUM_SPECIES = NS();
      integer constexpr NUM_DOF = 1+NS();
      AccumulationKernel< SUBREGION_TYPE, NUM_DOF, NUM_SPECIES, BASE_FLUID_TYPE > kernel( rankOffset, dofKey, subRegion, fluid, solid, dt, localMatrix, localRhs );
      AccumulationKernel< SUBREGION_TYPE, NUM_DOF, NUM_SPECIES, BASE_FLUID_TYPE >::template launch< POLICY >( subRegion.size(), kernel );
    } );
  }
};
 
} // namespace singlePhaseReactiveBaseKernels
 
} // namespace geos
 
#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_ACCUMULATIONKERNELS_HPP