<commit_before>/*! \file cooling_cuda.cu
 *  \brief Functions to calculate cooling rate for a given rho, P, dt. */

#ifdef CUDA
#ifdef COOLING_GPU

#include<cuda.h>
#include<math.h>
#include"global.h"
#include"global_cuda.h"
#include"io.h"
#include"cooling_cuda.h"



/*! \fn void cooling_kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, Real dt, Real gamma)
 *  \brief When passed an array of conserved variables and a timestep, adjust the value
           of the total energy for each cell according to the specified cooling function. */
//__global__ void cooling_kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, Real dt, Real gamma, cudaTextureObject_t coolTexObj, cudaTextureObject_t heatTexObj)
__global__ void cooling_kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, Real dt, Real gamma, Real *dt_array)
{
  __shared__ Real min_dt[TPB];

  int n_cells = nx*ny*nz;
  
  Real d, E;
  Real n, T, T_init;
  Real del_T, dt_sub;
  Real mu; // mean molecular weight
  Real cool; //cooling rate per volume, erg/s/cm^3
  #ifndef DE
  Real vx, vy, vz, p;
  #endif
  #ifdef DE
  Real ge;
  #endif
  //Real T_min = 1.0e4; // minimum temperature allowed
  Real T_min = 0.0; // minimum temperature allowed

  mu = 0.6;

  // get a global thread ID
  int blockId = blockIdx.x + blockIdx.y*gridDim.x;
  int id = threadIdx.x + blockId * blockDim.x;
  int zid = id / (nx*ny);
  int yid = (id - zid*nx*ny) / nx;
  int xid = id - zid*nx*ny - yid*nx;
  // and a thread id withing hte block
  int tid = threadIdx.x;

  // set min dt to current hydro timestep
  min_dt[tid] = dt;
  __syncthreads();
  

  // all threads do the calculation
  if (xid < nx && yid < ny && zid < nz) {

    // load values of density and pressure
    d  =  dev_conserved[            id];
    E  =  dev_conserved[4*n_cells + id];
    #ifndef DE
    vx =  dev_conserved[1*n_cells + id] / d;
    vy =  dev_conserved[2*n_cells + id] / d;
    vz =  dev_conserved[3*n_cells + id] / d;
    p  = (E - 0.5*d*(vx*vx + vy*vy + vz*vz)) * (gamma - 1.0);
    p  = fmax(p, (Real) TINY_NUMBER);
    #endif
    #ifdef DE
    ge = dev_conserved[5*n_cells + id] / d;
    ge = fmax(ge, (Real) TINY_NUMBER);
    #endif
    
    // calculate the number density of the gas (in cgs)
    n = d*DENSITY_UNIT / (mu * MP);

    // calculate the temperature of the gas
    #ifndef DE
    T_init = p*PRESSURE_UNIT/ (n*KB);
    #endif
    #ifdef DE
    T_init = ge*(gamma-1.0)*SP_ENERGY_UNIT*mu*MP/KB;
    #endif

    // calculate cooling rate per volume
    T = T_init;

    // call the cooling function (could choose primoridial cool)
    //cool = Schure_cool(n, T); 
    cool = Wiersma_cool(n, T); 
    //cool = primordial_cool(n, T);
    //cool = Cloudy_cool(n, T, coolTexObj, heatTexObj);
    
    // calculate change in temperature given dt
    del_T = cool*dt*TIME_UNIT*(gamma-1.0)/(n*KB);

    // limit change in temperature to 1%
    while (del_T/T > 0.01) {
      // what dt gives del_T = 0.01*T?
      dt_sub = 0.01*T*n*KB/(cool*TIME_UNIT*(gamma-1.0));
      // apply that dt
      T -= cool*dt_sub*TIME_UNIT*(gamma-1.0)/(n*KB);
      // how much time is left from the original timestep?
      dt -= dt_sub;
      // calculate cooling again
      //cool = Schure_cool(n, T);
      cool = Wiersma_cool(n, T);
      //cool = primordial_cool(n, T);
      //cool = Cloudy_cool(n, T, coolTexObj, heatTexObj);
      // calculate new change in temperature
      del_T = cool*dt*TIME_UNIT*(gamma-1.0)/(n*KB);
    }

    // calculate final temperature
    T -= del_T;

    // set a temperature floor
    if (T < T_min) { 
      T = T_min;
    }

    // adjust value of energy based on total change in temperature
    del_T = T_init - T; // total change in T
    E -= n*KB*del_T / ((gamma-1.0)*ENERGY_UNIT);
    #ifdef DE
    ge -= KB*del_T / (mu*MP*(gamma-1.0)*SP_ENERGY_UNIT);
    #endif
    if (del_T/T_init > 0.05) {
      //printf("%3d %3d %3d Cooling over 10 percent in hydro dt. T_init: %e T: %e\n", xid, yid, zid, T_init, T);
      // what dt gives del_T = 0.1*T?
      min_dt[tid] = 0.05*T*n*KB/(cool*TIME_UNIT*(gamma-1.0));
    }
    //if (T < 1000) printf("%3d %3d %3d Low T cell. T_init: %e T: %e\n", xid, yid, zid, T_init, T);

    // and send back from kernel
    dev_conserved[4*n_cells + id] = E;<commit_msg>new timestep limit applied every step<commit_after>/*! \file cooling_cuda.cu
 *  \brief Functions to calculate cooling rate for a given rho, P, dt. */

#ifdef CUDA
#ifdef COOLING_GPU

#include<cuda.h>
#include<math.h>
#include"global.h"
#include"global_cuda.h"
#include"io.h"
#include"cooling_cuda.h"



/*! \fn void cooling_kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, Real dt, Real gamma)
 *  \brief When passed an array of conserved variables and a timestep, adjust the value
           of the total energy for each cell according to the specified cooling function. */
//__global__ void cooling_kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, Real dt, Real gamma, cudaTextureObject_t coolTexObj, cudaTextureObject_t heatTexObj)
__global__ void cooling_kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, Real dt, Real gamma, Real *dt_array)
{
  __shared__ Real min_dt[TPB];

  int n_cells = nx*ny*nz;
  
  Real d, E;
  Real n, T, T_init;
  Real del_T, dt_sub;
  Real mu; // mean molecular weight
  Real cool; //cooling rate per volume, erg/s/cm^3
  #ifndef DE
  Real vx, vy, vz, p;
  #endif
  #ifdef DE
  Real ge;
  #endif
  //Real T_min = 1.0e4; // minimum temperature allowed
  Real T_min = 0.0; // minimum temperature allowed

  mu = 0.6;

  // get a global thread ID
  int blockId = blockIdx.x + blockIdx.y*gridDim.x;
  int id = threadIdx.x + blockId * blockDim.x;
  int zid = id / (nx*ny);
  int yid = (id - zid*nx*ny) / nx;
  int xid = id - zid*nx*ny - yid*nx;
  // and a thread id withing hte block
  int tid = threadIdx.x;

  // set min dt to current hydro timestep
  min_dt[tid] = dt;
  __syncthreads();
  

  // all threads do the calculation
  if (xid < nx && yid < ny && zid < nz) {

    // load values of density and pressure
    d  =  dev_conserved[            id];
    E  =  dev_conserved[4*n_cells + id];
    #ifndef DE
    vx =  dev_conserved[1*n_cells + id] / d;
    vy =  dev_conserved[2*n_cells + id] / d;
    vz =  dev_conserved[3*n_cells + id] / d;
    p  = (E - 0.5*d*(vx*vx + vy*vy + vz*vz)) * (gamma - 1.0);
    p  = fmax(p, (Real) TINY_NUMBER);
    #endif
    #ifdef DE
    ge = dev_conserved[5*n_cells + id] / d;
    ge = fmax(ge, (Real) TINY_NUMBER);
    #endif
    
    // calculate the number density of the gas (in cgs)
    n = d*DENSITY_UNIT / (mu * MP);

    // calculate the temperature of the gas
    #ifndef DE
    T_init = p*PRESSURE_UNIT/ (n*KB);
    #endif
    #ifdef DE
    T_init = ge*(gamma-1.0)*SP_ENERGY_UNIT*mu*MP/KB;
    #endif

    // calculate cooling rate per volume
    T = T_init;

    // call the cooling function (could choose primoridial cool)
    //cool = Schure_cool(n, T); 
    cool = Wiersma_cool(n, T); 
    //cool = primordial_cool(n, T);
    //cool = Cloudy_cool(n, T, coolTexObj, heatTexObj);
    
    // calculate change in temperature given dt
    del_T = cool*dt*TIME_UNIT*(gamma-1.0)/(n*KB);

    // limit change in temperature to 1%
    while (del_T/T > 0.01) {
      // what dt gives del_T = 0.01*T?
      dt_sub = 0.01*T*n*KB/(cool*TIME_UNIT*(gamma-1.0));
      // apply that dt
      T -= cool*dt_sub*TIME_UNIT*(gamma-1.0)/(n*KB);
      // how much time is left from the original timestep?
      dt -= dt_sub;
      // calculate cooling again
      //cool = Schure_cool(n, T);
      cool = Wiersma_cool(n, T);
      //cool = primordial_cool(n, T);
      //cool = Cloudy_cool(n, T, coolTexObj, heatTexObj);
      // calculate new change in temperature
      del_T = cool*dt*TIME_UNIT*(gamma-1.0)/(n*KB);
    }

    // calculate final temperature
    T -= del_T;
    //if (T < 1000) printf("%3d %3d %3d Low T cell. T_init: %e T: %e\n", xid, yid, zid, T_init, T);

    // set a temperature floor
    if (T < T_min) { 
      T = T_min;
    }

    // adjust value of energy based on total change in temperature
    del_T = T_init - T; // total change in T
    E -= n*KB*del_T / ((gamma-1.0)*ENERGY_UNIT);
    #ifdef DE
    ge -= KB*del_T / (mu*MP*(gamma-1.0)*SP_ENERGY_UNIT);
    #endif
    if (del_T/T_init > 0.1) {
      printf("%3d %3d %3d Cooling over 10 percent in hydro dt. T_init: %e T: %e\n", xid, yid, zid, T_init, T);
    }
    // calculate cooling rate for new T
    cool = Wiersma_cool(n, T);
    // limit the timestep such that delta_T is 10% 
    min_dt[tid] = 0.1*T*n*KB/(cool*TIME_UNIT*(gamma-1.0));

    // and send back from kernel
    dev_conserved[4*n_cells + id] = E;
{'fraction_non_alphanumeric': 0.13572166724078538, 'fraction_numerical': 0.014467791939373063, 'mean_word_length': 2.588792748248867, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3471329', 'n_tokens_mistral': 3551, 'n_tokens_neox': 3326, 'n_words': 992}