API tools faq
paste
Advertisement
mabruzzo

reformatted initialize_cloudy_data.c

mabruzzo
Jun 7th, 2025
84
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 10.67 KB | None | 0 0
raw download clone embed print report
  1. /***********************************************************************
  2. /
  3. / Initialize Cloudy cooling data
  4. /
  5. /
  6. / Copyright (c) 2013, Enzo/Grackle Development Team.
  7. /
  8. / Distributed under the terms of the Enzo Public Licence.
  9. /
  10. / The full license is in the file LICENSE, distributed with this
  11. / software.
  12. ************************************************************************/
  13.  
  14. #include <stdlib.h>
  15. #include <string.h>
  16. #include <math.h>
  17. #include "hdf5.h"
  18. #include "grackle.h"
  19. #include "grackle_macros.h"
  20. #include "grackle_types.h"
  21. #include "grackle_chemistry_data.h"
  22.  
  23. #define SMALL_LOG_VALUE -99.0
  24.  
  25. extern int grackle_verbose;
  26.  
  27. /**
  28.  * Initializes an empty #cloudy_data struct with zeros and NULLs.
  29.  */
  30. void initialize_empty_cloudy_data_struct(cloudy_data* my_cloudy) {
  31.   my_cloudy->grid_rank = 0LL;
  32.   for (long long i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++) {
  33.     my_cloudy->grid_dimension[i] = 0LL;
  34.     my_cloudy->grid_parameters[i] = NULL;
  35.   }
  36.   my_cloudy->heating_data = NULL;
  37.   my_cloudy->cooling_data = NULL;
  38.   my_cloudy->mmw_data = NULL;
  39.   my_cloudy->data_size = 0LL;
  40. }
  41.  
  42. // Initialize Cloudy cooling data
  43. int initialize_cloudy_data(chemistry_data* my_chemistry,
  44.                            chemistry_data_storage* my_rates,
  45.                            cloudy_data* my_cloudy, char* group_name,
  46.                            code_units* my_units, int read_data) {
  47.   long long q, w;
  48.   double* temp_data;
  49.   long long temp_int;
  50.   long long* temp_int_arr;
  51.   char parameter_name[MAX_LINE_LENGTH];
  52.  
  53.   // Initialize things (to the null-state) even if cloudy cooling is not used.
  54.   initialize_empty_cloudy_data_struct(my_cloudy);
  55.  
  56.   if (read_data == 0) {
  57.     return SUCCESS;
  58.   }
  59.  
  60.   if (grackle_verbose) {
  61.     fprintf(stdout, "Initializing Cloudy cooling: %s.\n", group_name);
  62.     fprintf(stdout, "cloudy_table_file: %s.\n",
  63.             my_chemistry->grackle_data_file);
  64.   }
  65.  
  66.   /* Get conversion units. */
  67.  
  68.   double co_length_units, co_density_units;
  69.   if (my_units->comoving_coordinates == TRUE) {
  70.     co_length_units = my_units->length_units;
  71.     co_density_units = my_units->density_units;
  72.   } else {
  73.     co_length_units =
  74.         my_units->length_units * my_units->a_value * my_units->a_units;
  75.     co_density_units =
  76.         my_units->density_units / POW(my_units->a_value * my_units->a_units, 3);
  77.   }
  78.  
  79.   double tbase1 = my_units->time_units;
  80.   double xbase1 = co_length_units / (my_units->a_value * my_units->a_units);
  81.   double dbase1 =
  82.       co_density_units * POW(my_units->a_value * my_units->a_units, 3);
  83.   double mh = 1.67e-24;
  84.   double CoolUnit = (POW(my_units->a_units, 5) * POW(xbase1, 2) * POW(mh, 2)) /
  85.                     (POW(tbase1, 3) * dbase1);
  86.  
  87.   // Read cooling data in from hdf5 file.
  88.   hid_t file_id, dset_id, attr_id;
  89.   herr_t status;
  90.   herr_t h5_error = -1;
  91.  
  92.   file_id =
  93.       H5Fopen(my_chemistry->grackle_data_file, H5F_ACC_RDONLY, H5P_DEFAULT);
  94.  
  95.   if (H5Aexists(file_id, "old_style")) {
  96.     my_rates->cloudy_data_new = 0;
  97.     if (grackle_verbose) fprintf(stdout, "Loading old-style Cloudy tables.\n");
  98.   }
  99.  
  100.   // Open cooling dataset and get grid dimensions.
  101.  
  102.   sprintf(parameter_name, "/CoolingRates/%s/Cooling", group_name);
  103.   dset_id = H5Dopen(file_id, parameter_name);
  104.   if (dset_id == h5_error) {
  105.     fprintf(stderr, "Can't open Cooling in %s.\n",
  106.             my_chemistry->grackle_data_file);
  107.     return FAIL;
  108.   }
  109.  
  110.   // Grid rank.
  111.   attr_id = H5Aopen_name(dset_id, "Rank");
  112.   if (attr_id == h5_error) {
  113.     fprintf(stderr, "Failed to open Rank attribute in Cooling dataset.\n");
  114.     return FAIL;
  115.   }
  116.   status = H5Aread(attr_id, HDF5_I8, &temp_int);
  117.   if (attr_id == h5_error) {
  118.     fprintf(stderr, "Failed to read Rank attribute in Cooling dataset.\n");
  119.     return FAIL;
  120.   }
  121.   my_cloudy->grid_rank = (long long)temp_int;
  122.   if (grackle_verbose)
  123.     fprintf(stdout, "Cloudy cooling grid rank: %lld.\n", my_cloudy->grid_rank);
  124.   status = H5Aclose(attr_id);
  125.   if (attr_id == h5_error) {
  126.     fprintf(stderr, "Failed to close Rank attribute in Cooling dataset.\n");
  127.     return FAIL;
  128.   }
  129.  
  130.   // Grid dimension.
  131.   temp_int_arr = malloc(my_cloudy->grid_rank * sizeof(long long));
  132.   attr_id = H5Aopen_name(dset_id, "Dimension");
  133.   if (attr_id == h5_error) {
  134.     fprintf(stderr, "Failed to open Dimension attribute in Cooling dataset.\n");
  135.     return FAIL;
  136.   }
  137.   status = H5Aread(attr_id, HDF5_I8, temp_int_arr);
  138.   if (attr_id == h5_error) {
  139.     fprintf(stderr, "Failed to read Dimension attribute in Cooling dataset.\n");
  140.     return FAIL;
  141.   }
  142.   if (grackle_verbose) fprintf(stdout, "Cloudy cooling grid dimensions:");
  143.   for (q = 0; q < my_cloudy->grid_rank; q++) {
  144.     my_cloudy->grid_dimension[q] = (long long)temp_int_arr[q];
  145.     if (grackle_verbose) fprintf(stdout, " %lld", my_cloudy->grid_dimension[q]);
  146.   }
  147.   if (grackle_verbose) fprintf(stdout, ".\n");
  148.   status = H5Aclose(attr_id);
  149.   if (attr_id == h5_error) {
  150.     fprintf(stderr,
  151.             "Failed to close Dimension attribute in Cooling dataset.\n");
  152.     return FAIL;
  153.   }
  154.   free(temp_int_arr);
  155.  
  156.   // Grid parameters.
  157.   for (q = 0; q < my_cloudy->grid_rank; q++) {
  158.     if (q < my_cloudy->grid_rank - 1) {
  159.       sprintf(parameter_name, "Parameter%lld", (q + 1));
  160.     } else {
  161.       sprintf(parameter_name, "Temperature");
  162.     }
  163.  
  164.     temp_data = malloc(my_cloudy->grid_dimension[q] * sizeof(double));
  165.  
  166.     attr_id = H5Aopen_name(dset_id, parameter_name);
  167.     if (attr_id == h5_error) {
  168.       fprintf(stderr, "Failed to open %s attribute in Cooling dataset.\n",
  169.               parameter_name);
  170.       return FAIL;
  171.     }
  172.     status = H5Aread(attr_id, HDF5_R8, temp_data);
  173.     if (attr_id == h5_error) {
  174.       fprintf(stderr, "Failed to read %s attribute in Cooling dataset.\n",
  175.               parameter_name);
  176.       return FAIL;
  177.     }
  178.  
  179.     my_cloudy->grid_parameters[q] =
  180.         malloc(my_cloudy->grid_dimension[q] * sizeof(double));
  181.     for (w = 0; w < my_cloudy->grid_dimension[q]; w++) {
  182.       if (q < my_cloudy->grid_rank - 1) {
  183.         my_cloudy->grid_parameters[q][w] = (double)temp_data[w];
  184.       } else {
  185.         // convert temeperature to log
  186.         my_cloudy->grid_parameters[q][w] = (double)log10(temp_data[w]);
  187.       }
  188.     }
  189.     if (grackle_verbose)
  190.       fprintf(stdout, "%s: %" GSYM " to %" GSYM " (%lld steps).\n",
  191.               parameter_name, my_cloudy->grid_parameters[q][0],
  192.               my_cloudy->grid_parameters[q][my_cloudy->grid_dimension[q] - 1],
  193.               my_cloudy->grid_dimension[q]);
  194.     status = H5Aclose(attr_id);
  195.     if (attr_id == h5_error) {
  196.       fprintf(stderr, "Failed to close %s attribute in Cooling dataset.\n",
  197.               parameter_name);
  198.       return FAIL;
  199.     }
  200.     free(temp_data);
  201.   }
  202.  
  203.   // Read Cooling data.
  204.   my_cloudy->data_size = 1;
  205.   for (q = 0; q < my_cloudy->grid_rank; q++) {
  206.     my_cloudy->data_size *= my_cloudy->grid_dimension[q];
  207.   }
  208.   temp_data = malloc(my_cloudy->data_size * sizeof(double));
  209.  
  210.   status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, temp_data);
  211.   if (grackle_verbose) fprintf(stdout, "Reading Cloudy Cooling dataset.\n");
  212.   if (status == h5_error) {
  213.     fprintf(stderr, "Failed to read Cooling dataset.\n");
  214.     return FAIL;
  215.   }
  216.  
  217.   my_cloudy->cooling_data = malloc(my_cloudy->data_size * sizeof(double));
  218.   for (q = 0; q < my_cloudy->data_size; q++) {
  219.     my_cloudy->cooling_data[q] = temp_data[q] > 0 ? (double)log10(temp_data[q])
  220.                                                   : (double)SMALL_LOG_VALUE;
  221.  
  222.     // Convert to code units.
  223.     my_cloudy->cooling_data[q] -= log10(CoolUnit);
  224.   }
  225.   free(temp_data);
  226.  
  227.   status = H5Dclose(dset_id);
  228.   if (status == h5_error) {
  229.     fprintf(stderr, "Failed to close Cooling dataset.\n");
  230.     return FAIL;
  231.   }
  232.  
  233.   // Read Heating data.
  234.   if (my_chemistry->UVbackground == 1) {
  235.     temp_data = malloc(my_cloudy->data_size * sizeof(double));
  236.  
  237.     sprintf(parameter_name, "/CoolingRates/%s/Heating", group_name);
  238.     dset_id = H5Dopen(file_id, parameter_name);
  239.     if (dset_id == h5_error) {
  240.       fprintf(stderr, "Can't open Heating in %s.\n",
  241.               my_chemistry->grackle_data_file);
  242.       return FAIL;
  243.     }
  244.  
  245.     status =
  246.         H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, temp_data);
  247.     if (grackle_verbose) fprintf(stdout, "Reading Cloudy Heating dataset.\n");
  248.     if (status == h5_error) {
  249.       fprintf(stderr, "Failed to read Heating dataset.\n");
  250.       return FAIL;
  251.     }
  252.  
  253.     my_cloudy->heating_data = malloc(my_cloudy->data_size * sizeof(double));
  254.     for (q = 0; q < my_cloudy->data_size; q++) {
  255.       my_cloudy->heating_data[q] = temp_data[q] > 0
  256.                                        ? (double)log10(temp_data[q])
  257.                                        : (double)SMALL_LOG_VALUE;
  258.  
  259.       // Convert to code units.
  260.       my_cloudy->heating_data[q] -= log10(CoolUnit);
  261.     }
  262.     free(temp_data);
  263.  
  264.     status = H5Dclose(dset_id);
  265.     if (status == h5_error) {
  266.       fprintf(stderr, "Failed to close Heating dataset.\n");
  267.       return FAIL;
  268.     }
  269.   }
  270.  
  271.   // Read MMW data.
  272.   if (my_chemistry->primordial_chemistry == 0 &&
  273.       strcmp(group_name, "Primordial") == 0) {
  274.     my_cloudy->mmw_data = malloc(my_cloudy->data_size * sizeof(double));
  275.  
  276.     sprintf(parameter_name, "/CoolingRates/%s/MMW", group_name);
  277.     dset_id = H5Dopen(file_id, parameter_name);
  278.     if (dset_id == h5_error) {
  279.       fprintf(stderr, "Can't open MMW in %s.\n",
  280.               my_chemistry->grackle_data_file);
  281.       return FAIL;
  282.     }
  283.  
  284.     status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT,
  285.                      my_cloudy->mmw_data);
  286.     if (grackle_verbose) fprintf(stdout, "Reading Cloudy MMW dataset.\n");
  287.     if (status == h5_error) {
  288.       fprintf(stderr, "Failed to read MMW dataset.\n");
  289.       return FAIL;
  290.     }
  291.  
  292.     status = H5Dclose(dset_id);
  293.     if (status == h5_error) {
  294.       fprintf(stderr, "Failed to close MMW dataset.\n");
  295.       return FAIL;
  296.     }
  297.   }
  298.  
  299.   status = H5Fclose(file_id);
  300.  
  301.   if (my_cloudy->grid_rank > GRACKLE_CLOUDY_TABLE_MAX_DIMENSION) {
  302.     fprintf(stderr,
  303.             "Error: rank of Cloudy cooling data must be less than or equal to "
  304.             "%d.\n",
  305.             GRACKLE_CLOUDY_TABLE_MAX_DIMENSION);
  306.     return FAIL;
  307.   }
  308.  
  309.   return SUCCESS;
  310. }
  311.  
  312. int _free_cloudy_data(cloudy_data* my_cloudy, chemistry_data* my_chemistry,
  313.                       int primordial) {
  314.   int i;
  315.  
  316.   for (i = 0; i < my_cloudy->grid_rank; i++) {
  317.     GRACKLE_FREE(my_cloudy->grid_parameters[i]);
  318.   }
  319.  
  320.   GRACKLE_FREE(my_cloudy->cooling_data);
  321.   if (my_chemistry->UVbackground == 1) {
  322.     GRACKLE_FREE(my_cloudy->heating_data);
  323.   }
  324.   if (my_chemistry->primordial_chemistry == 0 && primordial) {
  325.     GRACKLE_FREE(my_cloudy->mmw_data);
  326.   }
  327.   return GR_SUCCESS;
  328. }
  329.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!