load_cloudy_texture.h

Go to the documentation of this file.

 
#pragma once
 
#include <cmath>  // pow, log10
#include <string>
 
#include "../cooling/texture_utilities.h"  // Bilinear_Texture
#include "../global/global.h"
#include "../io/ParameterMap.h"
#include "../utils/shared.h"
 
namespace cool_component
{
 
class CloudyHeatAndCool
{
  // our usage of SharedHandle allows for the wrapped texture objects to be shared
  // among multiple owners (i.e. multiple copies of CloudyHeatAndCool) and ensures that
  // texture objects are properly cleaned up when the number of owners go to 0
  SharedHandle<cudaTextureObject_t> coolTexObj_;
  SharedHandle<cudaTextureObject_t> heatTexObj_;
 
 public:
  __host__ explicit CloudyHeatAndCool(std::string filename);
 
  __host__ explicit CloudyHeatAndCool(ParameterMap& pmap)
      : CloudyHeatAndCool(pmap.value_or("chemistry.data_file", ""))  // delegate to other constructor
  {
  }
 
  template <bool TABLE_ONLY>
  __device__ __forceinline__ Real calc_contrib_(Real n, Real T) const
  {
    Real lambda  = 0.0;  // log cooling rate, erg s^-1 cm^3
    Real cooling = 0.0;  // cooling per unit volume, erg /s / cm^3
    Real heating = 0.0;  // heating per unit volume, erg /s / cm^3
 
    // To keep texture code simple, we use floats (which have built-in support) as opposed to doubles (which would
    // require casting)
    float log_n, log_T;
    log_n = log10(n);
    log_T = log10(T);
    // this temp.
    if ((not TABLE_ONLY) and (log10(T) > 9.0)) {
      lambda = 0.45 * log10(T) - 26.065;
    } else if (TABLE_ONLY or (log10(T) >= 1.0)) {
      // remap coordinates for texture
      // remapped = (input - TABLE_MIN_VALUE)*(1/TABLE_SPACING)
      // remapped = (input - TABLE_MIN_VALUE)*(NUM_CELLS_PER_DECADE)
      const Real remap_log_T = (log_T - 1.0) * 10;
      const Real remap_log_n = (log_n + 6.0) * 10;
 
      lambda       = Bilinear_Texture(this->coolTexObj_.get(), remap_log_T, remap_log_n);
      const Real H = Bilinear_Texture(this->heatTexObj_.get(), remap_log_T, remap_log_n);
      heating      = pow(10, H);
    } else {
      // Do nothing below 10 K
      return 0.0;
    }
 
    cooling = pow(10, lambda);
    return n * n * (cooling - heating);
  }
 
  __device__ Real operator()(Real n, Real T) const { return calc_contrib_<false>(n, T); }
};
 
}  // namespace cool_component