ThermalSolutionScalingKernel.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_COMPOSITIONAL_THERMALSOLUTIONSCALINGKERNEL_HPP
#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_THERMALSOLUTIONSCALINGKERNEL_HPP
 
#include "physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingKernel.hpp"
 
namespace geos
{
 
namespace thermalCompositionalMultiphaseBaseKernels
{
 
/******************************** SolutionScalingKernel ********************************/
 
class SolutionScalingKernel : public isothermalCompositionalMultiphaseBaseKernels::SolutionScalingKernel
{
public:
 
  using Base = isothermalCompositionalMultiphaseBaseKernels::SolutionScalingKernel;
  using Base::m_numComp;
  using Base::m_localSolution;
 
  SolutionScalingKernel( real64 const maxRelativePresChange,
                         real64 const maxAbsolutePresChange,
                         real64 const maxRelativeTempChange,
                         real64 const maxAbsoluteTempChange,
                         real64 const maxCompFracChange,
                         real64 const maxRelativeCompDensChange,
                         globalIndex const rankOffset,
                         integer const numComp,
                         string const dofKey,
                         ElementSubRegionBase const & subRegion,
                         arrayView1d< real64 const > const localSolution,
                         arrayView1d< real64 const > const pressure,
                         arrayView1d< real64 const > const temperature,
                         arrayView2d< real64 const, compflow::USD_COMP > const compDens,
                         arrayView1d< real64 > pressureScalingFactor,
                         arrayView1d< real64 > compDensScalingFactor,
                         arrayView1d< real64 > temperatureScalingFactor,
                         integer const temperatureOffset )
    : Base( maxRelativePresChange,
            maxAbsolutePresChange,
            maxCompFracChange,
            maxRelativeCompDensChange,
            rankOffset,
            numComp,
            dofKey,
            subRegion,
            localSolution,
            pressure,
            compDens,
            pressureScalingFactor,
            compDensScalingFactor ),
    m_maxRelativeTempChange( maxRelativeTempChange ),
    m_maxAbsoluteTempChange( maxAbsoluteTempChange ),
    m_temperature( temperature ),
    m_temperatureScalingFactor( temperatureScalingFactor ),
    m_temperatureOffset( temperatureOffset )
  {}
 
  GEOS_HOST_DEVICE
  void compute( localIndex const ei,
                StackVariables & stack ) const
  {
    computeScalingFactor( ei, stack );
  }
 
  GEOS_HOST_DEVICE
  void computeScalingFactor( localIndex const ei,
                             StackVariables & stack ) const
  {
    real64 constexpr eps = isothermalCompositionalMultiphaseBaseKernels::minDensForDivision;
    Base::computeScalingFactor( ei, stack, [&] ()
    {
      // compute the change in temperature
      real64 const temp = m_temperature[ei];
      real64 const absTempChange = LvArray::math::abs( m_localSolution[stack.localRow + m_temperatureOffset] );
      if( stack.localMaxDeltaTemp < absTempChange )
      {
        stack.localMaxDeltaTemp = absTempChange;
      }
 
      // compute temperature scaling factor
      real64 tempScalingFactor = 1.0;
      // the stricter of the absolute and relative change limits is applied
      if( m_maxAbsoluteTempChange > 0.0 ) // maxAbsoluteTempChange <= 0.0 means that absolute scaling is disabled
      {
        if( absTempChange > m_maxAbsoluteTempChange )
        {
          tempScalingFactor = m_maxAbsoluteTempChange / absTempChange;
        }
      }
      if( temp > eps )
      {
        real64 const relativeTempChange = absTempChange / temp;
        if( relativeTempChange > m_maxRelativeTempChange )
        {
          tempScalingFactor = LvArray::math::min( tempScalingFactor, m_maxRelativeTempChange / relativeTempChange );
        }
      }
      m_temperatureScalingFactor[ei] = tempScalingFactor;
      if( stack.localMinVal > tempScalingFactor )
      {
        stack.localMinVal = tempScalingFactor;
      }
      if( stack.localMinTempScalingFactor > tempScalingFactor )
      {
        stack.localMinTempScalingFactor = tempScalingFactor;
      }
    } );
  }
 
protected:
 
  real64 const m_maxRelativeTempChange;
  real64 const m_maxAbsoluteTempChange;
 
  arrayView1d< real64 const > const m_temperature;
 
  arrayView1d< real64 > const m_temperatureScalingFactor;
 
  integer const m_temperatureOffset;
};
 
class SolutionScalingKernelFactory
{
public:
 
  template< typename POLICY >
  static SolutionScalingKernel::StackVariables
  createAndLaunch( real64 const maxRelativePresChange,
                   real64 const maxAbsolutePresChange,
                   real64 const maxRelativeTempChange,
                   real64 const maxAbsoluteTempChange,
                   real64 const maxCompFracChange,
                   real64 const maxRelativeCompDensChange,
                   arrayView1d< real64 const > const pressure,
                   arrayView1d< real64 const > const temperature,
                   arrayView2d< real64 const, compflow::USD_COMP > const compDens,
                   arrayView1d< real64 > pressureScalingFactor,
                   arrayView1d< real64 > compDensScalingFactor,
                   arrayView1d< real64 > temperatureScalingFactor,
                   globalIndex const rankOffset,
                   integer const numComp,
                   string const dofKey,
                   ElementSubRegionBase & subRegion,
                   arrayView1d< real64 const > const localSolution,
                   integer const temperatureOffset )
  {
    SolutionScalingKernel kernel( maxRelativePresChange, maxAbsolutePresChange, maxRelativeTempChange, maxAbsoluteTempChange,
                                  maxCompFracChange, maxRelativeCompDensChange,
                                  rankOffset, numComp, dofKey, subRegion, localSolution,
                                  pressure, temperature, compDens, pressureScalingFactor,
                                  compDensScalingFactor, temperatureScalingFactor, temperatureOffset );
    return thermalCompositionalMultiphaseBaseKernels::
             SolutionScalingKernel::launch< POLICY >( subRegion.size(), kernel );
  }
 
};
 
} // namespace thermalCompositionalMultiphaseBaseKernels
 
} // namespace geos
 
 
#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_THERMALSOLUTIONSCALINGKERNEL_HPP