ThermalAccumulationKernels.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_THERMALACCUMULATIONKERNELS_HPP
#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_THERMALACCUMULATIONKERNELS_HPP
 
#include "physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp"
 
namespace geos
{
 
namespace thermalSinglePhaseBaseKernels
{
 
/******************************** AccumulationKernel ********************************/
 
template< typename SUBREGION_TYPE, integer NUM_DOF >
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_volume;
  using Base::m_deltaVolume;
  using Base::m_porosity;
  using Base::m_dPoro_dPres;
  using Base::m_density;
  using Base::m_dDensity;
  using Base::m_localMatrix;
  using Base::m_localRhs;
 
  static constexpr integer isThermal = NUM_DOF-1;
  using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< isThermal >;
  AccumulationKernel( globalIndex const rankOffset,
                      string const dofKey,
                      SUBREGION_TYPE const & subRegion,
                      constitutive::SingleFluidBase const & fluid,
                      constitutive::CoupledSolidBase const & solid,
                      CRSMatrixView< real64, globalIndex const > const & localMatrix,
                      arrayView1d< real64 > const & localRhs )
    : Base( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ),
    m_dPoro_dTemp( solid.getDporosity_dTemperature() ),
    m_internalEnergy( fluid.internalEnergy() ),
    m_dInternalEnergy( fluid.dInternalEnergy() ),
    m_rockInternalEnergy( solid.getInternalEnergy() ),
    m_dRockInternalEnergy_dTemp( solid.getDinternalEnergy_dTemperature() ),
    m_energy_n( subRegion.template getField< fields::flow::energy_n >() )
  {}
 
  struct StackVariables : public Base::StackVariables
  {
public:
 
    GEOS_HOST_DEVICE
    StackVariables()
      : Base::StackVariables()
    {}
 
    using Base::StackVariables::poreVolume;
    using Base::StackVariables::dPoreVolume_dPres;
    using Base::StackVariables::localRow;
    using Base::StackVariables::dofIndices;
    using Base::StackVariables::localResidual;
    using Base::StackVariables::localJacobian;
 
    real64 dPoreVolume_dTemp = 0.0;
 
    // Solid energy
 
    real64 solidEnergy = 0.0;
 
    real64 dSolidEnergy_dPres = 0.0;
 
    real64 dSolidEnergy_dTemp = 0.0;
  };
 
 
  GEOS_HOST_DEVICE
  void setup( localIndex const ei,
              StackVariables & stack ) const
  {
    Base::setup( ei, stack );
 
    stack.dPoreVolume_dTemp = ( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dTemp[ei][0];
 
    // initialize the solid volume
    real64 const solidVolume = ( m_volume[ei] + m_deltaVolume[ei] ) * ( 1.0 - m_porosity[ei][0] );
    real64 const dSolidVolume_dPres = -( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dPres[ei][0];
    real64 const dSolidVolume_dTemp = -( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dTemp[ei][0];
 
    // initialize the solid internal energy
    stack.solidEnergy = solidVolume * m_rockInternalEnergy[ei][0];
    stack.dSolidEnergy_dPres = dSolidVolume_dPres * m_rockInternalEnergy[ei][0];
    stack.dSolidEnergy_dTemp = solidVolume * m_dRockInternalEnergy_dTemp[ei][0] + dSolidVolume_dTemp * m_rockInternalEnergy[ei][0];
  }
 
  GEOS_HOST_DEVICE
  void computeAccumulation( localIndex const ei,
                            StackVariables & stack ) const
  {
    stack.localResidual[numEqn-1] = -m_energy_n[ei];
 
    Base::computeAccumulation( ei, stack, [&] ()
    {
      // Step 1: assemble the derivatives of the mass balance equation w.r.t temperature
      stack.localJacobian[0][numDof-1] = stack.poreVolume * m_dDensity[ei][0][DerivOffset::dT] + stack.dPoreVolume_dTemp * m_density[ei][0];
 
      // Step 2: assemble the fluid part of the accumulation term of the energy equation
      real64 const fluidEnergy = stack.poreVolume * m_density[ei][0] * m_internalEnergy[ei][0];
 
      real64 const dFluidEnergy_dP = stack.dPoreVolume_dPres * m_density[ei][0] * m_internalEnergy[ei][0]
                                     + stack.poreVolume * m_dDensity[ei][0][DerivOffset::dP] * m_internalEnergy[ei][0]
                                     + stack.poreVolume * m_density[ei][0] * m_dInternalEnergy[ei][0][DerivOffset::dP];
 
      real64 const dFluidEnergy_dT = stack.poreVolume * m_dDensity[ei][0][DerivOffset::dT] * m_internalEnergy[ei][0]
                                     + stack.poreVolume * m_density[ei][0] * m_dInternalEnergy[ei][0][DerivOffset::dT]
                                     + stack.dPoreVolume_dTemp * m_density[ei][0] * m_internalEnergy[ei][0];
      // local accumulation
      stack.localResidual[numEqn-1] += fluidEnergy;
 
      // derivatives w.r.t. pressure and temperature
      stack.localJacobian[numEqn-1][0]        = dFluidEnergy_dP;
      stack.localJacobian[numEqn-1][numDof-1] = dFluidEnergy_dT;
    } );
 
    // Step 3: assemble the solid part of the accumulation term of the energy equation
    stack.localResidual[numEqn-1] += stack.solidEnergy;
    stack.localJacobian[numEqn-1][0] += stack.dSolidEnergy_dPres;
    stack.localJacobian[numEqn-1][numDof-1] += stack.dSolidEnergy_dTemp;
  }
 
  GEOS_HOST_DEVICE
  void complete( localIndex const ei,
                 StackVariables & stack ) const
  {
    // Step 1: assemble the mass balance equation
    Base::complete( ei, stack );
 
    // Step 2: assemble the energy equation
    m_localRhs[stack.localRow + numEqn-1] += stack.localResidual[numEqn-1];
    m_localMatrix.template addToRow< serialAtomic >( stack.localRow + numEqn-1,
                                                     stack.dofIndices,
                                                     stack.localJacobian[numEqn-1],
                                                     numDof );
  }
 
protected:
 
 
  arrayView2d< real64 const > const m_dPoro_dTemp;
 
  arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_internalEnergy;
  arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > const m_dInternalEnergy;
 
  arrayView2d< real64 const > const m_rockInternalEnergy;
  arrayView2d< real64 const > const m_dRockInternalEnergy_dTemp;
 
  arrayView1d< real64 const > const m_energy_n;
 
};
 
class SurfaceElementAccumulationKernel : public AccumulationKernel< SurfaceElementSubRegion, 2 >
{
 
public:
 
  using Base = AccumulationKernel< SurfaceElementSubRegion, 2 >;
 
  SurfaceElementAccumulationKernel( globalIndex const rankOffset,
                                    string const dofKey,
                                    SurfaceElementSubRegion const & subRegion,
                                    constitutive::SingleFluidBase const & fluid,
                                    constitutive::CoupledSolidBase const & solid,
                                    CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                    arrayView1d< real64 > const & localRhs )
    : Base( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ),
    m_creationMass( subRegion.getField< fields::flow::massCreated >() )
  {}
 
  GEOS_HOST_DEVICE
  void computeAccumulation( localIndex const ei,
                            Base::StackVariables & stack ) const
  {
    Base::computeAccumulation( ei, stack );
    if( Base::m_mass_n[ei] > 1.1 * m_creationMass[ei] )
    {
      stack.localResidual[0] += m_creationMass[ei] * 0.25;
    }
  }
 
protected:
 
  arrayView1d< real64 const > const m_creationMass;
 
};
 
class AccumulationKernelFactory
{
public:
 
  template< typename POLICY, typename SUBREGION_TYPE >
  static void
  createAndLaunch( globalIndex const rankOffset,
                   string const dofKey,
                   SUBREGION_TYPE const & subRegion,
                   constitutive::SingleFluidBase const & fluid,
                   constitutive::CoupledSolidBase const & solid,
                   CRSMatrixView< real64, globalIndex const > const & localMatrix,
                   arrayView1d< real64 > const & localRhs )
  {
    if constexpr ( std::is_base_of_v< CellElementSubRegion, SUBREGION_TYPE > )
    {
      integer constexpr NUM_DOF = 2;
      AccumulationKernel< CellElementSubRegion, NUM_DOF > kernel( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs );
      AccumulationKernel< CellElementSubRegion, NUM_DOF >::template launch< POLICY >( subRegion.size(), kernel );
    }
    else if constexpr ( std::is_base_of_v< SurfaceElementSubRegion, SUBREGION_TYPE > )
    {
      SurfaceElementAccumulationKernel kernel( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs );
      SurfaceElementAccumulationKernel::launch< POLICY >( subRegion.size(), kernel );
    }
    else
    {
      GEOS_UNUSED_VAR( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs );
      GEOS_ERROR( "Unsupported subregion type: " << typeid(SUBREGION_TYPE).name() );
    }
  }
 
};
 
} // namespace thermalSinglePhaseBaseKernels
 
} // namespace geos
 
#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_THERMALACCUMULATIONKERNELS_HPP