ThermalSinglePhaseWellKernels.hpp

/*
 * ------------------------------------------------------------------------------------------------------------
 * 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_WELLS_THERMALSINGLEPHASEWELLKERNELS_HPP
#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_THERMALSINGLEPHASEWELLKERNELS_HPP
 
#include "constitutive/fluid/singlefluid/SingleFluidFields.hpp"
#include "constitutive/fluid/singlefluid/SingleFluidBase.hpp"
#include "common/DataTypes.hpp"
#include "common/GEOS_RAJA_Interface.hpp"
#include "mesh/ElementRegionManager.hpp"
#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp"
#include "physicsSolvers/fluidFlow/StencilAccessors.hpp"
#include "physicsSolvers/fluidFlow/wells/WellControls.hpp"
#include "physicsSolvers/PhysicsSolverBaseKernels.hpp"
 
namespace geos
{
 
namespace thermalSinglePhaseWellKernels
{
 
 
 
/******************************** ElementBasedAssemblyKernel ********************************/
 
template< integer NUM_DOF >
class ElementBasedAssemblyKernel : public singlePhaseWellKernels::ElementBasedAssemblyKernel< NUM_DOF >
{
public:
  using Base = singlePhaseWellKernels::ElementBasedAssemblyKernel< NUM_DOF >;
  using Base::m_rankOffset;
  using Base::m_wellElemDofNumber;
  using Base::m_elemGhostRank;
  using Base::m_wellElemVolume;
  using Base::m_wellElemDensity;
  using Base::m_wellElemDensity_n;
  using Base::m_dWellElemDensity_dPressure;
  using Base::m_localMatrix;
  using Base::m_localRhs;
  using ROFFSET = singlePhaseWellKernels::RowOffset;
  using COFFSET = singlePhaseWellKernels::ColOffset;
 
  static constexpr integer numDof = NUM_DOF;
 
  static constexpr integer numEqn = NUM_DOF;
 
  ElementBasedAssemblyKernel( globalIndex const rankOffset,
                              string const dofKey,
                              ElementSubRegionBase const & subRegion,
                              constitutive::SingleFluidBase const & fluid,
                              CRSMatrixView< real64, globalIndex const > const & localMatrix,
                              arrayView1d< real64 > const & localRhs )
    :    Base( rankOffset, dofKey, subRegion, fluid, localMatrix, localRhs ),
    m_dWellElemDensity_dTemperature( fluid.dDensity_dTemperature()  ),
    m_internalEnergy( fluid.internalEnergy() ),
    m_internalEnergy_n( fluid.internalEnergy_n() ),
    m_dInternalEnergy_dPres( fluid.dInternalEnergy_dPressure() ),
    m_dInternalEnergy_dTemp( fluid.dInternalEnergy_dTemperature() )
  {}
 
  struct StackVariables : public Base::StackVariables
  {
public:
    GEOS_HOST_DEVICE
    StackVariables()
      : Base::StackVariables()
    {}
    using Base::StackVariables::eqnRowIndices;
    using Base::StackVariables::dofColIndices;
    using Base::StackVariables::localJacobian;
    using Base::StackVariables::localResidual;
    using Base::StackVariables::localRow;
    using Base::StackVariables::volume;
    using Base::StackVariables::density;
    using Base::StackVariables::density_n;
    using Base::StackVariables::dDensity_dPres;
 
  };
  GEOS_HOST_DEVICE
  integer elemGhostRank( localIndex const ei ) const
  { return m_elemGhostRank( ei ); }
 
 
 
  template< typename FUNC = NoOpFunc >
  GEOS_HOST_DEVICE
  void computeAccumulation( localIndex const iwelem,
                            StackVariables & stack ) const
  {
    Base::computeAccumulation( iwelem, stack, [&]( )
    {
 
      // Step 1: assemble the derivatives of the mass balance equation w.r.t temperature
      stack.localJacobian[0][numDof-1] =  stack.volume * m_dWellElemDensity_dTemperature[iwelem][0];
 
      // Step 2: assemble the fluid part of the accumulation term of the energy equation
      real64 const fluidEnergy = stack.volume   * stack.density  * m_internalEnergy[iwelem][0];
      real64 const fluidEnergy_n = stack.volume   * stack.density_n  * m_internalEnergy_n[iwelem][0];
 
      real64 const dFluidEnergy_dP =  stack.volume   * stack.dDensity_dPres  * m_internalEnergy[iwelem][0]
                                     + stack.volume   * stack.density  * m_dInternalEnergy_dPres[iwelem][0];
 
 
      real64 const dFluidEnergy_dT = stack.volume   * m_dWellElemDensity_dTemperature[iwelem][0] * m_internalEnergy[iwelem][0]
                                     + stack.volume  * stack.density  * m_dInternalEnergy_dTemp[iwelem][0];
 
      // local accumulation
      stack.localResidual[numEqn-1] = fluidEnergy - fluidEnergy_n;
 
      // derivatives w.r.t. pressure and temperature
      stack.localJacobian[numEqn-1][0]        = dFluidEnergy_dP;
      stack.localJacobian[numEqn-1][numDof-1] = dFluidEnergy_dT;
    } );
  }
 
 
 
  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 iwelem )
    {
      if( kernelComponent.elemGhostRank( iwelem ) >= 0 )
      {
        return;
      }
      typename KERNEL_TYPE::StackVariables stack;
      kernelComponent.setup( iwelem, stack );
      kernelComponent.computeAccumulation( iwelem, stack );
      kernelComponent.complete( iwelem, stack );
 
    } );
  }
 
protected:
 
  arrayView2d< real64 const > const m_dWellElemDensity_dTemperature;
 
  arrayView2d< real64 const > const m_internalEnergy;
  arrayView2d< real64 const > const m_internalEnergy_n;
  arrayView2d< real64 const > const m_dInternalEnergy_dPres;
  arrayView2d< real64 const > const m_dInternalEnergy_dTemp;
 
};
 
 
class ElementBasedAssemblyKernelFactory
{
public:
  template< typename POLICY >
  static void
  createAndLaunch( globalIndex const rankOffset,
                   string const dofKey,
                   ElementSubRegionBase const & subRegion,
                   constitutive::SingleFluidBase const & fluid,
                   CRSMatrixView< real64, globalIndex const > const & localMatrix,
                   arrayView1d< real64 > const & localRhs )
  {
    integer constexpr NUM_DOF = 2;
    ElementBasedAssemblyKernel< NUM_DOF >
    kernel( rankOffset, dofKey, subRegion, fluid, localMatrix, localRhs );
    ElementBasedAssemblyKernel< NUM_DOF >::template
    launch< POLICY, ElementBasedAssemblyKernel< NUM_DOF > >( subRegion.size(), kernel );
 
  }
};
} // end namespace singlePhaseWellKernels
 
} // end namespace geos
 
#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_SINGLEPHASEWELLKERNELS_HPP