Giter VIP home page Giter VIP logo

grackle's Introduction

Grackle

Users' Mailing List CircleCI Documentation Status

Grackle is a chemistry and radiative cooling library for astrophysical simulations and models. Grackle has interfaces for C, C++, Fortran, and Python codes and provides:

  • two options for primordial chemistry and cooling:

    1. non-equilibrium primordial chemistry network for atomic H, D, and He as well as H2 and HD, including H2 formation on dust grains.

    2. tabulated H and He cooling rates calculated with the photo-ionization code, Cloudy.

  • tabulated metal cooling rates calculated with Cloudy.

  • photo-heating and photo-ionization from two UV backgrounds with optional self-shielding corrections:

    1. Faucher-Giguere et al. (2009).

    2. Haardt & Madau (2012).

  • support for user-provided arrays of volumetric and specific heating rates.

The Grackle provides functions to update chemistry species; solve radiative cooling and update internal energy; and calculate cooling time, temperature, pressure, and ratio of specific heats (gamma).

For more information on features, installation, and integration with simulation codes and models, see our online documentation.

Software that uses Grackle

A (non-exhaustive) list of software that provides out-of-the-box support for using Grackle includes:

We welcome PRs to add your simulation code to this list. We also welcome the inclusion of python modules that depend on Pygrackle.

Resources

grackle's People

Contributors

aemerick avatar andreaincatasciato avatar avifriedlander avatar bkoh1986 avatar brittonsmith avatar bvillasen avatar choi-junhwan avatar chummels avatar clairekope avatar ewanbjones98 avatar galtay avatar gregbryan avatar h3jia avatar hyschive avatar jwise77 avatar leonardromano avatar loikki avatar lq3552 avatar mabruzzo avatar matthewturk avatar mladenivkovic avatar peeples avatar qobilidop avatar scog1234 avatar sleitner avatar snigdaa avatar tiapac avatar weiguangcui avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

grackle's Issues

Additional parallelism in solve_rate_cool

Currently, the scale_fields, ceiling_fields, and make_consistent functions are not parallelized the way the rest of solve_rate_cool is. These are all fully local and should be straightforward to do with OpenMP in a similar strategy to the rest of the code. This might increase the parallel efficiency of solve_chemistry.

Line length issues

Commit d062c19 from pull request #94 produces lines longer than 72 characters at a couple of places in the fortran code. One is line 170 of cool1d_multi_g.F:

   anydust = (idust .gt. 0) .or. (idustall .gt. 0) .or. (idustrec .gt. 0)

The other is line 12 of solve_rate_cool_g.F:

 &                is, js, ks, ie, je, ke, ih2co, ipiht, idustrec, igammah, 

This breaks compilation if you're using the default line length for compiling fixed format files (i.e. if you're not passing the compiler the flag '-ffixed-line-length-132'). I'm raising this as an issue rather than just patching it because I'm not sure whether it's better to fix the line lengths or to simply mandate the use of the '-ffixed-line-length-132' flag.

Doc clarification

In order to avoid confusion, I suggest to clearly write in the doc that the variables density, HI_density, HII_density, etc. are mass density. This is important to distinguish them from e_density which, if I'm correct represents the number density of electron. I can do that if you agree.
I don't find the note:

Note, the electron mass density should be scaled by the ratio of the proton mass to the electron mass such that the electron density in the code is the electron number density times the proton mass.

clear enough...

Grackle with Float

Originally reported by loic hausammann (Bitbucket: loichausammann, GitHub: Unknown)


Hi,

I am trying to use grackle, but I am getting some large errors when trying to reproduce the graphs in S. Britton 2017 (see attachments).

I have been able to reproduce them (with some differences that I suppose come from the code improvements) with double, but with floats, the computation fails at large temperature.

My graphs where done with the HM2012 UV background and primordial_chemistry = 0.

Thank you for your help and time.


Incorrect scaling in low density H2 cooling rate?

Originally reported by Simon Glover (Bitbucket: scog, GitHub: scog)


In cool1d_multi_g.F, we use the following code to compute the total H2 cooling rate in the low density limit when we're using the Glover & Abel (2008) treatment:

        galdl(i) = gaHI(i) * HI(i,j,k)  + gaH2(i) * H2I(i,j,k)
 &               + gaHe(i) * HeI(i,j,k) + gaHp(i) * HII(i,j,k)
 &               + gael(i) * de(i,j,k)

Since HI, H2I, HeI etc. are mass fractions, don't we need to divide the H2 term by 2 and the He term by 4? Or am I misunderstanding something?


Can only link against lgrackle-2.0.1 on makefile

Originally reported by Nathan Goldbaum (Bitbucket: ngoldbaum, GitHub: ngoldbaum)


For some reason, I can no longer use the idiom "-lgrackle", and must instead use "-lgrackle-2.0.1" to link against grackle from enzo. There is a valid symlink in my grackle install directory from libgrackle.dylib to libgrackle-2.0.1.dylib.

This is on OSX, compiling enzo with clang. I can share my grackle and enzo makefiles if that would help.


Rename master branch to main

This has been done for most other projects. I was hoping to wait until all open PRs have been merged, in particular for the 3.2 release. Just creating this issue so I don't forget.

Issue with metal correction to mmw for tabulated cooling

I think I may have stumbled onto a small error in the equation used to compute the correction to the mean molecular weight to account for the contribution of metals, when primordial_chemistry =0. Specifically, I'm talking about the logic described in section 4.1.2 of the Grackle Paper (I confirmed that the exact same logic is used in the code).

At risk of repeating information you readily remember, I'll provide some context. When using Grackle in Tabulated mode, it iteratively solves for mean molecular weight and temperature in the following equation (equation 12 of the paper):
e = kB * T / ((gamma - 1) * μi(nh, T) * mh)
where μi(nh,T) is the tabulated mean molecular weight of the primordial elements.

Because μi only accounts for the primordial elements, a correction is required to account for the contribution of the metals' mean molecular weight, μZ=16. The correction used by Grackle is
μT = ρT / (ρT / μi + ρZ / μZ)
where μT is the total mean molecular mass of the gas, ρT is the total mass density of the gas, and ρZ is the mass density of the metals. This equation is provided in equation 14 of the paper.

This is slightly different from the correction I would naively derive. I would start from:
μT = ρT / (mh * (ni + nZ) )
where mh is the mass of hydrogen, ni is the number density of the primordial elements, and nZ is the number density of the metals. Using ni=(ρT - ρZ)/(mh*μi) and nZ=ρZ/(mh*μZ). I find:
μT = ρT / ((ρT - ρZ) / μi + ρZ / μZ).

At first, I thought that the version of the correction used by Grackle is just an approximation of my version of the correction when ρZ << ρT. However, it turns out that Grackle's version leads to μT<=μi while my version leads to μT>=μi.

This becomes more apparent if we manipulate the form of the Grackle correction. This is shown below and uses x=ρZ/ρT:
μT = ρT / (ρT / μi + ρZ / μZ)
μT = 1 / (1 / μi + x / μZ)
μT = 1 / ( (μZ + μi * x) / (μi * μZ) )
μT = μi * μZ / (μZ + μi * x)
Finally, we can take the ratio of μT to μi:
μT / μi = μZ / (μZ + μi * x)
As one would expect, the ratio is 1 when x=0. However, the ratio is less than one for all positive x.

The comparable form of my version of the correction is:
μT / μi = μZ / ( (1 - x) * μZ + μi * x)
Again, the ratio is 1 when x=0. Because μi < μZ, the ratio is always larger than 1 for 0 < x < 1.

Am I missing something obvious here?

Updated Cloudy tables for Self-shielding Approximation

Originally reported by Andrew Emerick (Bitbucket: aemerick, GitHub: aemerick)


Currently, using self-shielding approximation with the current cloudy cooling tables is highly inconsistent, and can lead to substantial overestimate of the total cooling rate in density regimes where self-shielding is important. This is because the Cloudy tables are computed as optically thin, leading to an overestimation of the metal line cooling rates. This can be a substantial difference.

I have already computed new Cloudy tables which account for self-shielding which fix this issue. The to-do is to move these tables into a more friendly format. This means:

  1. Take the current optically thin tables and remove all associated self-shielding data from them (the grey averaged cross sections)

  2. Put together both the HM2012 self-shielding table as separate from the above (done)

  3. Put together an equivalent FG2011 table (lower priority)

  4. Add into the code catches for the user if they are attempting to use the self-shielding approximations with the old tables.


Tabulated UVB/cooling needs refactoring

There is so low hanging refactoring fruit in the tabulated UV background and cooling code that yield some minor optimizations, cut down on duplicate code, and standardize a few things. The list I can think of so far includes:

  • refactor initialize_cloudy_data.c to use the read_dataset function implemented in initialize_UVbackground_data.c
  • UVB rates should be converted to code units and logspace upon read in initialize_UVbackground_data.c like the cooling tables. This would prevent having to do this every time in update_UVbackground_rates.c
  • update_UVbackground_rates.c could use a more efficient of calculating the interpolation index than the brute force for loop. For tables evenly spaced in log(1 + z), we can do this with algebra. Otherwise, bisection would also be an improvement.

Recombination cooling

Originally reported by Simon Glover (Bitbucket: scog, GitHub: scog)


Grackle's recombination cooling rates currently assume case A, regardless of whether or not case B recombination has been selected by the user.

I plan to fix this myself once I have time, but I'm opening an issue to make sure I don't forget.


Compile error : wrong auto_get_version.c

In grackle 3.2dev2, a new method to calculate the version information is adopted.
However, this will cause a compile error if the user directly downloads grackle from the website.
This error doesn't appear if the user downloads it through git pull.
Need another way to generate auto_get_version.c or temp.show-version.

Narrative docs for new rate functions

PR #87 added individual functions for each of Grackle's non-equilibrium chemistry and cooling rates. Some documentation should be added demonstrating how these can be used.

Temperature increasing from Radiative Cooling

I've found a weird bug in which the Temperature of a ~7.5e5 K gas at a mass density of 1.3557e-27 g/cm^3 increases with time (the UV background has been turned off). I first encountered this problem while using Enzo-E and have been able to reproduce the error with the python bindings of the latest version. Specifically, the bug can be reproduced with this script.

The script evolves the gas over a single time-step (which is ~10 orders of magnitude shorter than the cooling time). Interestingly, the metal density appears to increase over the timestep by 3 orders of magnitudes.

As an aside: if I switch my densities to be in units of hydrogen masses, the Temperature does decrease over the timestep, but the pressure is 14 orders of magnitudes higher (However, this might be because I am misunderstanding how the units work).

Cooling at high redshift

Hi,

It seems that we discovered a bug in the cooling at high redshift (before the reionization). When running a cosmological simulation, the cooling converges to a heating at 1000 K if the initial temperature is high (see following image).

high_temp

We can see it in the following image of the cooling rate as function of the temperature. Either we converge to the real "minimum" or we converge to a local one at 1000K.

image

H2 self-shielding only works for 3D grid codes

Originally reported by Britton Smith (Bitbucket: brittonsmith, GitHub: brittonsmith)


The calculation of the length-scale for H2 self-shielding needs to be updated to work with non-cartesian grid codes.

Below is @gbryan's comment from PR #51:

However, there is an important point -- the pseudo-Sobelov method here only works for cartesian grid codes, and Grackle is in wide use in codes that violate both of those assumptions. Ideally, we would instead pass in a "length" field (the pseudo-Sobelov method is basically using the local gradient to define a length scale) which would be computed by the code and passed in. Obviously, that is quite a bit more work, so I'm happy to accept this PR but I think (i) the documentation should indicate the H2 self-shielding only works for Cartesian grid codes, and (ii) the code should at least error-check that the we are using 3D; that won't find all calling errors but at least will flag 0D and particle codes (typically 1D).

(An intermediate possibility would be to the H2 self-shielding have multiple non-zero values, and have a variety of ways to compute the length scale - for example, one common way is to use the Jeans length -- we tested that in Wolcott-Green et al. and found that it was good at high and low densities but overestimated the shielding at intermediate densities. However, again, this is more work and could be reasonably kicked down the road...).


gr_float not defined in simulation codes

Originally reported by Hsi-Yu Schive (Bitbucket: hyschive, GitHub: hyschive)


If I understand it correctly, simulation codes in C/C++ need to include "grackle.h", which then includes "grackle_types.h". However, simulation codes, in general, may not define either CONFIG_BFLOAT_4 nor CONFIG_BFLOAT_8, and thus the "gr_float" type in "grackle_types.h" becomes undefined.

This is not an issue for either Enzo or Grackle test codes since they all define CONFIG_BFLOAT_4(8) properly.


calculate_cooling_time.c:31:13: error: ‘FORTRAN_NAME’ declared as function returning a function

Originally reported by Anonymous


I'm building this for the first time and running into this error from "make":

libtool: compile: gcc -c -DH5_USE_16_API -DCONFIG_BFLOAT_8 -I. calculate_cooling_time.c -DDLL_EXPORT -DPIC -o .libs/calculate_cooling_time.o
calculate_cooling_time.c:31:13: error: ‘FORTRAN_NAME’ declared as function returning a function
extern void FORTRAN_NAME(cool_multi_time_g)(
^~~~~~~~~~~~
calculate_cooling_time.c:70:9: warning: parameter names (without types) in function declaration
int *iVheat, int *iMheat, gr_float *Vheat, gr_float *Mheat);
^~~
calculate_cooling_time.c: In function ‘_calculate_cooling_time’:
calculate_cooling_time.c:127:18: error: ‘cool_multi_time_g’ undeclared (first use in this function); did you mean ‘cooling_time’?
FORTRAN_NAME(cool_multi_time_g)(
^~~~~~~~~~~~~~~~~
cooling_time
calculate_cooling_time.c:127:18: note: each undeclared identifier is reported only once for each function it appears in
make: *** [Makefile:175: calculate_cooling_time.lo] Error 1

I'm running this using Cygwin (which I just installed) on Windows 10 (up-to-date). I'm not very familiar with this environment, unfortunately, as I usually use MacOS but my GPU card is in this Windows box.

Is something supposed to substitute for FORTRAN_NAME, which failed earlier? "make" started compiling the C routines and has not worked on the Fortran files yet. I have gfortran and g95 installed.

Thanks for your help.


itmax should be a runtime parameter

Currently, the maximum iterations in solve_rate_cool is set with a constant. There are situations where increasing this value from the default of 10000 is warranted. This would be easier if itmax was a settable parameter.

Inconsistency in UVB rates with self-shielding

Originally reported by Andrew Emerick (Bitbucket: aemerick, GitHub: aemerick)


I've uncovered a fairly significant issue when including H_2 formation with the self-shielding prescription turned on that leads to runaway H_2 formation at number densities above about n = 10 cm^-3 (the exact point is somewhat metallicity dependent).

The issue boils down to the fact that the rate coefficients from the UVB are split into independent bands (k24 - k31). Although we account for shielding of HI, HeI, and HeII radiation in (k24, k25, k26) with the Rahmati et. al. prescription, and account for H2 self-shielding against LW radiation (k31), the rest of the bands are not modified in any fashion.

This becomes a problem in self-shielding regions that start forming molecular hydrogen. Some of this H2I becomes ionized to H2II through rate k29, generating free electrons that go into making H-, which in turn makes more H2. The H2II quickly collides with the abundant neutral HI, making H2I and HII. The net result of this is a continual source of free electrons (H2I->H2II) that drives further H2 formation via the gas phase. This spikes the electron fraction to ~percent levels, and the H2I fraction to ~50-70% (by mass).

The reality, however, is that direct H2I ionization (15.2 eV) should be shielded significantly by the neutral HI. In addition (though not related to this problem) is the H2II dissociation radiation (30 eV) should also be shielded, mostly by HeI shielding.

I propose implementing shielding of these bands (rates k29, k30, and k28) using the same scaling provided by the Rahmati et. al. presectption, attenuating the 15.2 eV band using the same reduction factor as HI, and attenuating the 30 eV band using the same reduction factor as HeI shielding.

I'm including two plots showing the evolution of all H mass fractions over time in a one-zone cooling cell model, with an initial number density of 10, low metal mass fraction (2% solar), including H2 self-shielding from LW, and with ``self_shielding_method = 3''. The top is what happens with our current settings, the bottom shows what happens when including the shielding as I proposed above.

fiducial_H2_overproduction.png

fiducial_additional_shielding.png


Docs test failing due to sphinx issue

It appears that there is an issue with sphinx (we are using sphinx<3 - not sure why?) and a new version of jinja2. It looks like switching to jinja2<3.1 might fix this?

Introducing an API for dynamically configuring ``chemistry_data``?

I've been thinking about this for a while, and I wanted to gauge how people would feel about introducing a small handful of functions to the public API that allow external programs to dynamically configure the chemistry_data struct, without directly accessing the struct fields.

New Public Functions

To be a little more concrete, I'm suggesting the introduction of something like the following 3 functions:

  // we could also define versions of the following functions that operate on
  // the global grackle_data variable

  int local_chemistry_data_assign_double(chemistry_data* my_chemistry,
				         const char* target_name, double val);
  int local_chemistry_data_assign_int(chemistry_data* my_chemistry,
				      const char* target_name, int val);
  // this last one may not be worth adding:
  int local_chemistry_data_assign_string(chemistry_data* my_chemistry,
				         const char* target_name, char* val);

The above functions would all generally do the same thing (Each has behavior specialized for a different datatype). In short, the above functions would each:

  • Check if the chemistry_data struct has a field called target_name. It would also check if that field matches the type that the function is specialized to handle.

  • If one of the above conditions isn't satisfied, the function returns FAIL.

  • Otherwise, the function will perform my_chemistry->target_name = val and
    returns SUCCESS.

Benefits

This API generally provides 2 benefits to external simulation codes that depend on Grackle:

  1. This makes it possible for simulation codes to directly forward user specified options from a parameter file directly to Grackle.

    • I'm not sure about other codes, but currently in Enzo-E, users configure Grackle by specifying the name of fields in the chemistry_data struct parameter and the associated values. When Enzo-E initializes the chemistry_data struct, there is code that manually goes through (just about) every field in the chemistry_data struct, checks for a parameter named after the field, and then modifies the value of the chemistry_data struct, accordingly.

    • This change would let us replace the majority of this manual code, with much simpler code that simply forward the name, value pairs directly to Grackle. An additional benefit of this direct forwarding is that we would no longer need to introduce new code for configuring newly added chemistry_data fields.

  2. This lets simulation codes retain backwards compatibility with older versions of Grackle. Currently, when a code adds support for configuring newly added fields of chemistry_data, it generally loses the ability to be compiled with older versions of Grackle (since you can't compile have code that tries to access a member of a struct that doesn't exist). There are definitely occasions where I would have personally found this to be useful (when I was using a personal branch of Grackle).

Implementation Challenges

I suspect that the main argument against introducing these functions to the API is that they would introduce another place in the codebase that needs to be updated any time a new field is introduced to chemistry_data.

With that in mind, I wanted to point out that there is a fairly concise way to do this using a well-established pattern called X-Macros (see here or here for more details). This approach involves maintaining a central list of each struct member that gets reused for implementing multiple functions that requires knowledge of every struct field (A lot of examples show how it can be used for serialization of a struct or printing all the fields of a struct).

I have actually prototyped a solution that leverages this approach. It requires that we maintain a central list of each struct member (and some of its properties) that can be reused in various contexts. You can find an example of this list here for a subset of the fields of chemistry_data. Each entry in the example holds:

  • the name of the field
  • a description of that field's type (i.e. whether it's an integer, floating point value or a string)
  • the default value of that struct-field (this isn't strictly necessary)

I have provided some examples below that make use of this list. The examples assume that the list is defined in a file called grackle_chemistry_data_fields.def (with some minor tweaks, we could also include the list directly within the same file where the functions that use it are defined). The example API functions are shown in the following spoiler tag.

Example Implementation

A faster implementation that does less work at runtime (and doesn't involve casts to/from void*) is definitely possible (it would just involve more code). Since configuration just happens once, this probably doesn't need to be super fast...

#define GET_FIELD_PTR(FIELD,FIELD_TYPE,target,target_type,chem_ptr,out_ptr)   \
  if ((strcmp(#FIELD, target_name) == 0) &                                    \
      (strcmp(#FIELD_TYPE, target_type) == 0)){                               \
    out_ptr = (void*)&(chem_ptr->FIELD);                                      \
  }

// retrieves a pointer to the field of my_chemistry that's named target_name.
//
// This returns a NULL pointer if: my_chemistry is NULL, the field doesn't
// exist, or the field doesn't have the expected type (INT, DOUBLE, STRING)
void* _chem_data_field_ptr(chemistry_data* my_chemistry,
			   const char* target_name, const char* target_type)
{
  if (my_chemistry == NULL) { return NULL; }

  void* field_ptr = NULL;
  // now try to retrieve field_ptr from chemistry_data
  #define ENTRY(FIELD, TYPE, DEFAULT_VAL) \
    GET_FIELD_PTR (FIELD,TYPE,target_name,target_type,my_chemistry,field_ptr)
  #include "grackle_chemistry_data_fields.def"
  #undef ENTRY

  return field_ptr;
}

int local_chemistry_data_assign_double(chemistry_data* my_chemistry,
				       const char* target_name, double val)
{
  void* field_ptr = _chem_data_field_ptr(my_chemistry, target_name, "DOUBLE");
  if (field_ptr == NULL) { return FAIL; }
  *((double*)field_ptr) = val;
  return SUCCESS;
}

int local_chemistry_data_assign_int(chemistry_data* my_chemistry,
				    const char* target_name, int val)
{
  void* field_ptr = _chem_data_field_ptr(my_chemistry, target_name, "INT");
  if (field_ptr == NULL) { return FAIL; }
  *((int*)field_ptr) = val;
  return SUCCESS;
}

While you still need to keep the X-Macro List up to date with all the fields in chemistry_data, you can use this list to drastically simplify the implementation of some other existing functions The 2 main examples are provided below:

Simplified _set_default_chemistry_parameters function

The current implementation can be found here

chemistry_data _set_default_chemistry_parameters(void)
{
  chemistry_data my_chemistry;

  #define ENTRY(FIELD, TYPE, DEFAULT_VAL) my_chemistry.FIELD = DEFAULT_VAL;
  #include "grackle_chemistry_data_fields.def"
  #undef ENTRY

  return my_chemistry;
}
Simplified show_parameters function

The current implementation can be found here

void _show_field_INT(FILE *fp, const char* field, int val)
{ fprintf(fp, "%-33s = %d\n", field, val); }
void _show_field_DOUBLE(FILE *fp, const char* field, double val)
{ fprintf(fp, "%-33s = %g\n", field, val); }
void _show_field_STRING(FILE *fp, const char* field, const char* val)
{ fprintf(fp, "%-33s = %s\n", field, val); }

void show_parameters(FILE *fp, chemistry_data *my_chemistry){
  #define ENTRY(FIELD, TYPE, DEFAULT_VAL) \
    _show_field_ ## TYPE (fp, #FIELD, my_chemistry->FIELD);
  #include "grackle_chemistry_data_fields.def"
  #undef ENTRY
}

In case you find it helpful, here is a link that put's the above code snippets together into a toy executable (Just untar the folder and call gcc examples.c to compile the executable).

I'm happy to personally create the PR introducing this change, but I wanted to gauge interest in this before I do so (it wouldn't involve much more work). In particular, I realize this X-Macro approach may be a little controversial and that there could be other objections to expanding the API.

Photoheating from radiative transfer

Originally reported by Yves Revaz (Bitbucket: revaz, GitHub: revaz)


Hi all,

Currently, when using the heating rates from radiative transfer solutions, only the heating due to HI is included through the variable photogamma:

in cool1d_multi_g.F
! Photoheating from radiative transfer
edot(i) = edot(i) + real(ipiht, DKIND) * photogamma(i,j,k) / coolunit * HI(i,j,k) /dom

However, when computing the radiative heating due to the UB background through the rates provided by the cloudy table, grackle also include the He heating:

in cool1d_multi_g.F
! Photoionization heating
edot(i) = edot(i) + real(ipiht, DKIND)*(
& piHI *HI (i,j,k) ! pi of HI
& + piHeI *HeI (i,j,k)0.25_DKIND ! pi of HeI
& + piHeII
HeII(i,j,k)*0.25_DKIND ! pi of HeII
& )/dom

For the sake of consistency, in the radiative transfert mode, I suggest to let the possibility to provide piHI, piHeI and piHeII, replacing photogamma.

Doing so will still be compatible with the current version, if setting piHeI=0 and piHeII=0.

I think there is very few lines to change and if this is fine, I can certainly to the modifications.

yves


Wrong numbering of the rates

Hi,

Here you are saying that Grackle is using the same convention than in Abel et al. 1997 https://github.com/grackle-project/grackle/blob/master/src/clib/calc_rates_g.F#L78

If I compare the equations written in the comment below, they do not correspond with the paper mentioned (both the arxiv and published paper). For example for the photoionization and dissociation processes, only the equation 28 corresponds.

Is it just an error in the comments or should we expect a large bug?

If it is only the comment how do you want to fix it? I see three possibilities. We can remove the mention to the paper, add a comment on each wrong line or fix the coefficients everywhere.

If you are fine with one of the two first propositions, I can do the merge request.

Best,

Incorrect pygrackle fluid container temperature

Found an issue within pygrackle:
The convenience.py function setup_fluid_container sets the energy within a fluid_container object using a supplied temperature and the calculated mean molecular weight.
However, if you try to calculate the mean molecular weight before the energy has been set, the mean molecular weight defaults to 1.
This results in temperatures within the fluid_container that are lower than what the user supplied to setup_fluid_container, e.g.:
input temp [ 10000. 1000000. ]
calc'd temp [ 6324.14017413 632414.0174129 ]

Rather than setting the mean molecular weight to 1 if the energy has not been set, I think there should be a check for species densities that mu can be calculated from instead.
I'll be addressing this, but I wanted to make an issue to solicit any suggestions or comments.

Numerical Error from Unit Conversion

Changing the a_value in the chemistry_data struct creates changes in the results of the one-zone collapse of freefall.py that isn't due to any change in physics. For example this image shows the temperature during a collapse using a slightly modified version of freefall.py. The modifications are:

  1. velocity_units are changed to remove the dependence on a_value
  2. An approximate UV background is added so that H2 should be suppressed. The UV background is not redshift dependent and just changes the two photorates listed in the title of the plot

Example_Issue

The only physics difference between the z=0 and z=15 runs should be from Compton cooling however turning off Compton cooling does not change either of the results shown.

I believe the redshift dependence is introduced through the conversion of proper coordinates to comoving coordinates and back again. When chemistry_data is setup in terms of proper coordinates, the units are converted to comoving coordinates in solve_chemistry.c (lines 147-150) and then converted back to proper coordinates in solve_rate_g.F (lines 323-328) and cool_1D_multi.F (lines 184-189). When I tested this by overwriting those conversions so that the units always stay as proper coordinates the redshift dependence goes away and the z=15 line in the attached plot matches the z=0 line.

My best guess of what is happening is that the conversion between comoving and proper coordinates introduces a small floating point precision error that propagats over many iterations builds up and becomes a larger numerical issue. If this is the actual issue (I am not 100% positive) then the solution may be to add a comoving_coordinates flag so that gets passed through the code along with the units so that the units do not have to be converted back and forth.

Photoionizing sources in rate equations for 6-species network

Hi,

I am having difficulty understanding precisely what sources of photons are included in the photoionization terms of the rate equations for atomic H, He chemistry (i.e. the 6-species network).

Smith et al. (2017) says that the extragalactic UV background is included, but are there any other sources of photoelectrons in the chemistry? Are there, e.g., soft X-ray background photons/stellar EUV photons included, since those dominate H+ formation in the solar neighbourhood warm neutral medium, according to Wolfire et al. (2003)?

Global structure and shared memory

Originally reported by loic hausammann (Bitbucket: loichausammann, GitHub: Unknown)


Hi,

I am trying to run grackle in parallel (shared memory). But unfortunately, the code is not able to deal with my parallelization due to the global structure. Would it be possible to expose grackle_data and grackle_rates in this function https://bitbucket.org/grackle/grackle/src/7f7306c55f7bd7530c10a915cbeb809d7542cc2f/src/clib/solve_chemistry.c#solve_chemistry.c-237?

I have done a test and just exposing them is enough to solve the bug.

Thank you for your time and help.


UVbackground & Self-shielding

Originally reported by Valentin Perret (Bitbucket: vperret, GitHub: vperret)


Hi,

I recently implemented Grackle 3.0 in the ramses code. I noticed a couple of bugs during this implementation:

  • if grackle is initialized with UVbackground and self-shielding at a redshift higher than the one covered by the Cloudy tables, the k2* rates stay to 0 and the code produces NaNs.

  • the user parameters UVbackground_redshift_on and UVbackground_redshift_fullon are overridden in initialize_UVbackgroud.c

Thank you for your hard work!

Cheers,

Valentin Perret


Modify some dates in Haardt & Madau (2012) date

Originally reported by Edward Leong (Bitbucket: Edward0402, GitHub: Unknown)


Excuse me , I have a small problem, probably a trivial problem.
I try to do a small modification in the CloudyData_UVB=HM2012.h5. I hope to increase the HeII heating rate which comes from Haardt & Madau (2012). Since the file is binary , I don't know how to change it. Is there a easy way to achieve the goal?
Sincerely


k16 old/new rate

Originally reported by John Regan (Bitbucket: john_regan, GitHub: Unknown)


k16 rate in calc_rates_g.F in the old rates #define section

This is probably just me being stupid but this is an "old" rate based on Croft et al. (1999) - which is the same as quoted in Glover & Savin (2009). However the new rate in coll_rates_g.F is from Abel et al. 1997 (based on Delgarno & Lepp (1987) - the rate isn't exactly as quoted in the paper but its pretty close. Is this correct?


Free allocated memory

Dear Britton,

I am currently writing a piece of code that needs to initialize multiple times Grackle and quickly run into a lack of RAM.

I tried to look at the code, but it seems that the C API does not have a function to free the memory allocated by Grackle. Is it possible to add a function to cleanup the memory?

I am able to "fix" this problem with a dirty trick, therefore there is no rush for this issue.

Thank you very much for your time.

Best regards,

FATAL error (2) in MULTI_COOL

Originally reported by Romeel Davé (Bitbucket: romeeld, GitHub: romeeld)


When I try to run cosmological sims at high resolution or dynamic range, I inevitably see lots of errors like the ones below from Grackle. Unfortunately it doesn't happen in small test runs so it's hard to debug quickly.

This obviously comes from solve_rate_cool(), but I am not sure how to interpret the values, or tell why grackle is having a hard time with these particular particles. It may be that cooling in those particles could be ignored (eg if they are dense enough to be artificially pressurized), but I can't tell which ones they are.

Can anyone offer thoughts on how one might debug this, or why this might be happening? Thanks.

inside if statement solve rate cool: 0 0
FATAL error (2) in MULTI_COOL
dt = 1.391E-05 ttmin = 7.688E-06
6.6E-11
7.7E-06
7.9E+09
T
inside if statement solve rate cool: 0 0
FATAL error (2) in MULTI_COOL
dt = 1.391E-05 ttmin = 7.431E-06
4.3E-10
7.4E-06
2.1E+09
T


Cloudy version for the table

Hi,

I am currently trying to reproduce the table given in the repository (Cloudy_UVB=HM2012.h5) and I am not able to find the version of Cloudy used. Do you remember it?

Thank you for your help and time.
Best regards,

Replacing For-Loops with a Macro in C functions

This issue documents a way to potentially enhance the codebase that briefly discussed in PR #106.

Overview

The idea is to replace nested for-loop in (C/C++) code with a macro of the format: FULL_LOOP(index, index_helper, omp_pragma), where

  • index is the name of the variable that is being iterated over
  • index_helper is a previously initialized instance of grackle_index_helper
  • omp_pragma is a string-literal holding the OpenMP-related pragma arguments. This is conditionally applied (when _OPENMP is defined) to the outer-loop using the C99 _Pragma keyword.

A major motivation for doing something like this is to facilitate the testing of different loop structures for GPU parallelism in the future.

For concreteness, I'll show how one of the loops in local_calculate_pressure would change.

Effective current implementation (with changes from PR #106)
  double tiny_number = 1.e-20;
  const grackle_index_helper ind_helper = _build_index_helper(my_fields);

# ifdef _OPENMP
# pragma omp parallel for schedule( runtime )
# endif
  for (int outer_ind = 0; outer_ind < ind_helper.outer_ind_size; outer_ind++){

    const grackle_index_range range = _inner_range(outer_ind, &ind_helper);

    for (int index = range.start; index <= range.end; index++) {

      pressure[index] = ((my_chemistry->Gamma - 1.0) *
			 my_fields->density[index] *
			 my_fields->internal_energy[index]);
 
      if (pressure[index] < tiny_number)
        pressure[index] = tiny_number;
    } // end: loop over i
  } // end: loop over outer_ind
Equivalent implementation with FULL_LOOP macro
  double tiny_number = 1.e-20;
  const grackle_index_helper ind_helper = _build_index_helper(my_fields);

  FULL_LOOP(index, ind_helper, "omp parallel for schedule( runtime )") {

    pressure[index] = ((my_chemistry->Gamma - 1.0) *
                       my_fields->density[index] *
                       my_fields->internal_energy[index]);
 
    if (pressure[index] < tiny_number)
      pressure[index] = tiny_number;
  }

Sample Implementation

An example implementation of this macro is provided below (I have confirmed that implementation works - with & without OpenMP):

#ifdef _OPENMP
#define APPLY_OMP_PRAGMA(statement) _Pragma(statement)
#else
#define APPLY_OMP_PRAGMA(statement) /* ... */
#endif

#define GET_J_(outer_index, index_helper)                             \
  ((outer_index % index_helper.num_j_inds) + index_helper.j_start)
#define GET_K_(outer_index, index_helper)                             \
  ((outer_index / index_helper.num_j_inds) + index_helper.k_start)

#define FULL_LOOP(index, index_helper, omp_pragma)                        \
  APPLY_OMP_PRAGMA(omp_pragma)                                            \
  for (int oind_ = 0; oind_ < ind_helper.outer_ind_size; oind_++)         \
    for (int iind_ = index_helper.i_start,                                \
         index = (index_helper.i_start + index_helper.i_dim *             \
                      (GET_J_(oind_, index_helper) + index_helper.j_dim * \
                       GET_K_(oind_, index_helper)));                     \
         iind_ <= index_helper.i_end; iind_++, index++)

I'm not a huge fan of declaring/incrementing 2 variables inside the inner-loop, but doing it this way makes the application of the loop a lot cleaner. I have confirmed that the above macro definitely works (both with and without openmp), by converting both loops in local_calculate_pressure.

At this point, it would not take much effort to finish implementing this change, so let me know if you want me to do that (in PR #106 or a separate PR). However, it may make more sense to wait until we really need this for GPU-experimentation

H number density in self-shielding ignores H2

Originally reported by Andrew Emerick (Bitbucket: aemerick, GitHub: aemerick)


In computing the total hydrogen number density for the Rahamti et. al. self-shielding method, H2 is ignored. This would only be problematic in regions that are almost entirely molecular.

This should be:

#!fortran

nratio = (HI(i,j,k) + HII(i,j,k) + 2.0 * (H2I(i,j,k) + H2II(i,j,k))) * dom / nSSh

but is currently

#!fortran

nratio = (HI(i,j,k) + HII(i,j,k)) * dom / nSSh

Happy to do the PR. For full consistency we should probably also include HM and HD, but they should always be a small component and we can likely ignore them.


`local_calculate_pressure` ignores `grid_start` & `grid_end`

Currently, local_calculate_pressure ignores the values stored in the grid_start and grid_end members of the grackle_field_data . This is a departure from the way that most other functions work and does not appear to be documented.

Is this an intentional choice? Would there be any opposition to changing this? I have a use-case where this change would actually be useful.

If there isn't any opposition to changing this I'm happy to make a PR to address this

does not compile/link on macOS Big Sur?

I cannot successfully build Grackle on macOS Big Sur using the included Make.mach.darwin machine file.

I do ./configure in the root directory, then cd src/clib, then make machine-darwin:

clib git:(3d55b4a) ✗ make show-config

   MACHINE: Darwin (OSX) 
   MACHINE-NAME: darwin

   CONFIG_PRECISION  [precision-{32,64}]                     : 64
   CONFIG_OPT  [opt-{warn,debug,high,aggressive}]            : high
   CONFIG_OMP  [omp-{on,off}]                                : off

Then I type make and it fails when linking:

clib git:(3d55b4a) ✗ cat out.compile
glibtool: compile:  gcc -c -DLINUX -DH5_USE_16_API -DCONFIG_BFLOAT_8 -O2 -I/usr/local/hdf5/include -I. update_UVbackground_rates.c  -fno-common -DPIC -o .libs/update_UVbackground_rates.o
glibtool: compile:  gcc -c -DLINUX -DH5_USE_16_API -DCONFIG_BFLOAT_8 -O2 -I/usr/local/hdf5/include -I. update_UVbackground_rates.c -o update_UVbackground_rates.o >/dev/null 2>&1
glibtool: link: clang -dynamiclib -Wl,-undefined -Wl,dynamic_lookup -o .libs/libgrackle-3.2.dev2.dylib  .libs/auto_show_config.o .libs/auto_show_flags.o .libs/auto_show_version.o .libs/calculate_cooling_time.o .libs/calculate_dust_temperature.o .libs/calculate_gamma.o .libs/calculate_pressure.o .libs/calculate_temperature.o .libs/calc_temp1d_cloudy_g.o .libs/calc_temp_cloudy_g.o .libs/calc_rates_g.o .libs/calc_tdust_1d_g.o .libs/calc_tdust_3d_g.o .libs/cie_thin_cooling_rate_g.o .libs/colh2diss_g.o .libs/coll_rates_g.o .libs/cool1d_cloudy_g.o .libs/cool1d_cloudy_old_tables_g.o .libs/cool1d_multi_g.o .libs/cool_multi_time_g.o .libs/initialize_chemistry_data.o .libs/initialize_cloudy_data.o .libs/initialize_UVbackground_data.o .libs/interpolators_g.o .libs/set_default_chemistry_parameters.o .libs/solve_chemistry.o .libs/solve_rate_cool_g.o .libs/update_UVbackground_rates.o   -L/usr/local/hdf5/lib -lhdf5 -L/usr/local/lib -lgfortran    -install_name  /Users/benwibking/local/lib/libgrackle-3.2.dev2.dylib  -Wl,-single_module
ld: warning: directory not found for option '-L/usr/local/hdf5/lib'
ld: library not found for -lgfortran
clang: error: linker command failed with exit code 1 (use -v to see invocation)

Add a CONTRIBUTING file and code of conduct

Originally reported by Nathan Goldbaum (Bitbucket: ngoldbaum, GitHub: ngoldbaum)


I just noticed that there isn't any developer documentation.

We should consider adding a CONTRIBUTING file. There's some discussion in the Github documentation about what this file is:

https://help.github.com/articles/setting-guidelines-for-repository-contributors/

Here's a blog post from Sarah Sharp (former linux kernel maintainer) discussing ways projects can make themselves more welcoming:

http://sarah.thesharps.us/2015/10/06/what-makes-a-good-community/

Adding a code of conduct and basic contribution guidelines are level 0 and level 1 on her cultural levels towards making a project more diverse and welcoming.


pygrackle needs proper versioning

The pygrackle package has had a version of 0.1 since its inception. This should be updated for the next release and incremented properly when changed.

Disagreement between documentation and code

Originally reported by Simon Glover (Bitbucket: scog, GitHub: scog)


In Parameters.rst, the default values for the h2_optical_depth_approximation and cie_cooling variables are stated to be 0. However, the default values actually set in set_default_chemistry_parameters.c are both 1. Should we adjust the code to match the docs here, or the docs to match the code? (i.e. what do we actually want the defaults to be?)


calc_temp1d_cloudy_g can return negative molecular weights

This is with grackle-3.1.1 (commit 00c1724)
This using the CloudyData_UVB=HM2012_shielded.h5 data file for the recent AGORA code comparison runs. I stuck in if (munew .le. 0.0) call abort at line 161 in calc_temp1d_cloudy_g.F, the abort got tripped with the following parameters:

(gdb) p log_n_h(1)
$69 = 4.1023389670398664
(gdb) p zr
$70 = 18.36997029960439
(gdb) p log10tem(1)
$71 = 0.34712524596478833

This gives munew = -0.10427047396950628.
This looks like the interpolation is not correctly handling values beyond the cloudy grid: the largest grid log_n_h is 4, the largest grid z is 18, and the smallest grid log10tem is 1.

Apologies for not submitting a full script to reproduce this. But I can give more information if needed.

Branching error picked up by valgrind

I was running valgrind to track down unrelated issues and it picked up a potential bug in calc_tdust_1d_g.F. I've pasted the error message here to keep this brief: https://pastebin.com/Eg5p3XPJ

In short, I believe it is saying that the if statements in calc_kappa_gr_g are branching based on an uninitialized value of tdust, stemming from the call in line 217/218. Here, bi_t_mid is passed as tdust and bi_itmask is used as the masking flag field. I believe the issue is that certain indeces of tdustnow (which is used to calculate bi_t_mid in line 212) are not set under certain conditions in the if-statements in lines 99-112.

Frankly this is all a bit abstruse to me, as I don't really have too much of an idea as to what is going on in this routine. Is this an issue? Or is valgrind being oversensitive?

Edit: forgot to add link... added now...

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.