Manual merging of cooling_grackle (work done by @yvesrevaz and @mculpo)

configure.ac

@@ -421,6 +421,45 @@ AC_SUBST([METIS_LIBS])
 AC_SUBST([METIS_INCS])
 AM_CONDITIONAL([HAVEMETIS],[test -n "$METIS_LIBS"])
 
+# Check for grackle. 
+have_grackle="no"
+AC_ARG_WITH([grackle],
+    [AS_HELP_STRING([--with-grackle=PATH],
+       [root directory where grackle is installed @<:@yes/no@:>@]
+    )],
+    [],
+    [with_grackle="no"]
+)
+if test "x$with_grackle" != "xno"; then
+   AC_PROG_FC
+   AC_FC_LIBRARY_LDFLAGS
+   if test "x$with_grackle" != "xyes" -a "x$with_grackle" != "x"; then
+      GRACKLE_LIBS="-L$with_grackle/lib -lgrackle $FCLIBS"
+      GRACKLE_INCS="-I$with_grackle/include"
+   else
+      GRACKLE_LIBS="-lgrackle $FCLIBS"
+      GRACKLE_INCS=""
+   fi
+   have_grackle="yes"
+   AC_CHECK_LIB(
+      [grackle],
+      [initialize_grackle],
+      [AC_CHECK_HEADERS(
+         [grackle.h grackle_types.h grackle_macros.h],
+         [AC_DEFINE([HAVE_GRACKLE],1,[The GRACKLE library appears to be present.]),
+	  AC_DEFINE([CONFIG_BFLOAT_8],1,[Use doubles in grackle])
+	  ],
+         [AC_MSG_ERROR(Cannot find grackle headers!)]
+      )],
+      [AC_MSG_ERROR(Cannot find grackle library!)],
+      [$GRACKLE_LIBS]
+   )
+fi
+AC_SUBST([GRACKLE_LIBS])
+AC_SUBST([GRACKLE_INCS])
+AM_CONDITIONAL([HAVEGRACKLE],[test -n "$GRACKLE_LIBS"])
+
+
 #  Check for tcmalloc a fast malloc that is part of the gperftools.
 have_tcmalloc="no"
 AC_ARG_WITH([tcmalloc],
@@ -948,6 +987,7 @@ AC_MSG_RESULT([
    Metis enabled    : $have_metis
    FFTW3 enabled    : $have_fftw3
    libNUMA enabled  : $have_numa
+   GRACKLE enabled  : $have_grackle
    Using tcmalloc   : $have_tcmalloc
    Using jemalloc   : $have_jemalloc
    CPU profiler     : $have_profiler

examples/Makefile.am

@@ -27,8 +27,8 @@ AM_LDFLAGS = $(HDF5_LDFLAGS)
 EXTRA_LIBS = $(HDF5_LIBS) $(FFTW_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS)
 
 # MPI libraries.
-MPI_LIBS = $(METIS_LIBS) $(MPI_THREAD_LIBS)
-MPI_FLAGS = -DWITH_MPI $(METIS_INCS)
+MPI_LIBS = $(METIS_LIBS) $(GRACKLE_LIBS) $(MPI_THREAD_LIBS)
+MPI_FLAGS = -DWITH_MPI $(METIS_INCS) $(GRACKLE_INCS)
 
 # Programs.
 bin_PROGRAMS = swift

src/Makefile.am

@@ -28,8 +28,8 @@ GIT_CMD = @GIT_CMD@
 EXTRA_LIBS = $(HDF5_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS)
 
 # MPI libraries.
-MPI_LIBS = $(METIS_LIBS) $(MPI_THREAD_LIBS)
-MPI_FLAGS = -DWITH_MPI $(METIS_INCS)
+MPI_LIBS = $(METIS_LIBS) $(GRACKLE_LIB) $(MPI_THREAD_LIBS)
+MPI_FLAGS = -DWITH_MPI $(METIS_INCS) $(GRACKLE_LIB)
 
 # Build the libswiftsim library
 lib_LTLIBRARIES = libswiftsim.la

src/cooling/grackle/cooling.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_COOLING_GRACKLE_H
+#define SWIFT_COOLING_GRACKLE_H
+
+/**
+ * @file src/cooling/none/cooling.h
+ * @brief Empty infrastructure for the cases without cooling function
+ */
+
+/* Some standard headers. */
+#include <float.h>
+#include <math.h>
+
+/* Local includes. */
+#include "error.h"
+#include "hydro.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/* include the grackle wrapper */
+#include "grackle_wrapper.h"
+
+
+
+
+
+/*! This function computes the new entropy due to the cooling,
+ *  between step t0 and t1.
+ */
+
+static INLINE double DoCooling_GRACKLE(double energy, double density, double dtime, double *ne, double Z, double a_now)
+{
+
+  
+
+  /*********************************************************************
+   call to the main chemistry solver
+   *********************************************************************/
+  
+  if (wrap_do_cooling(density, &energy, dtime,Z, a_now) == 0) {
+    fprintf(stderr, "Error in do_cooling.\n");
+    return 0;
+  }
+  
+    
+  return energy;
+
+
+}
+
+
+
+
+
+
+/**
+ * @brief Apply the cooling function to a particle.
+ *
+ * We do nothing.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cooling The #cooling_function_data used in the run.
+ * @param p Pointer to the particle data.
+ * @param dt The time-step of this particle.
+ */
+__attribute__((always_inline)) INLINE static void cooling_cool_part(
+    const struct phys_const* restrict phys_const,
+    const struct UnitSystem* restrict us,
+    const struct cooling_function_data* restrict cooling,
+    struct part* restrict p, struct xpart* restrict xp, float dt) {
+    
+
+
+
+  /* Get current internal energy (dt=0) */
+  const float u_old = hydro_get_internal_energy(p, 0.f);
+  /* Get current density */
+  const float rho = hydro_get_density(p);
+  /* Actual scaling fractor */
+  const float a_now =  1. / (1. + cooling->GrackleRedshift); ;			/*  must be chaged !!! */
+
+  double ne,Z;
+
+  Z =  0.02041;					/* 0.02041 (= 1 Zsun in Grackle v2.0, but = 1.5761 Zsun in Grackle v2.1) */
+  ne=  0.0;					/* mass fraction of eletron */          	/* useless for GRACKLE_CHEMISTRY = 0 */
+ 
+ 
+  float u_new;
+  float delta_u;
+  
+  u_new = DoCooling_GRACKLE(u_old, rho, dt, &ne, Z, a_now);
+  //u_new = u_old * 0.99;
+
+
+
+
+
+  
+  //if (u_new < 0)
+  //if (p->id==50356)
+  //  printf("WARNING !!! ID=%llu  u_old=%g  u_new=%g rho=%g dt=%g ne=%g Z=%g a_now=%g\n",p->id,u_old,u_new,rho,dt,ne,Z,a_now);
+
+  
+  delta_u = u_new - u_old;
+
+  /* record energy lost */
+  xp->cooling_data.radiated_energy += -delta_u * hydro_get_mass(p);  
+  
+ 
+  
+  /* Update the internal energy */
+  hydro_set_internal_energy(p, u_new);  
+
+
+
+
+    
+}
+
+/**
+ * @brief Computes the cooling time-step.
+ *
+ * We return FLT_MAX so as to impose no limit on the time-step.
+ *
+ * @param cooling The #cooling_function_data used in the run.
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float cooling_timestep(
+    const struct cooling_function_data* restrict cooling,
+    const struct phys_const* restrict phys_const,
+    const struct UnitSystem* restrict us, const struct part* restrict p) {
+
+  return FLT_MAX;
+}
+
+/**
+ * @brief Sets the cooling properties of the (x-)particles to a valid start
+ * state.
+ *
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void cooling_init_part(
+    const struct part* restrict p, struct xpart* restrict xp) {
+
+  xp->cooling_data.radiated_energy = 0.f;    
+}
+
+/**
+ * @brief Returns the total radiated energy by this particle.
+ *
+ * @param xp The extended particle data
+ */
+__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
+    const struct xpart* restrict xp) {
+
+  return xp->cooling_data.radiated_energy;
+}
+
+/**
+ * @brief Initialises the cooling properties.
+ *
+ * @param parameter_file The parsed parameter file.
+ * @param us The current internal system of units.
+ * @param phys_const The physical constants in internal units.
+ * @param cooling The cooling properties to initialize
+ */
+static INLINE void cooling_init_backend(
+    const struct swift_params* parameter_file, const struct UnitSystem* us,
+    const struct phys_const* phys_const,
+    struct cooling_function_data* cooling) {
+    
+
+
+  char cloudytable[200];
+  double units_density, units_length, units_time;
+  int grackle_chemistry;
+  int UVbackground;
+
+
+  parser_get_param_string(parameter_file,"GrackleCooling:GrackleCloudyTable",cooling->GrackleCloudyTable);
+  cooling->UVbackground = parser_get_param_int(parameter_file, "GrackleCooling:UVbackground");
+  cooling->GrackleRedshift = parser_get_param_double(parameter_file, "GrackleCooling:GrackleRedshift");
+  cooling->GrackleHSShieldingDensityThreshold = parser_get_param_double(parameter_file, "GrackleCooling:GrackleHSShieldingDensityThreshold");
+
+  // FIXME : Why a strcpy ?
+  strcpy(cloudytable,cooling->GrackleCloudyTable);
+  
+  
+  UVbackground  =  cooling->UVbackground;
+  grackle_chemistry = 0;			/* forced to be zero : read table */ 
+  
+  units_density = us->UnitMass_in_cgs/pow(us->UnitLength_in_cgs,3);
+  units_length  = us->UnitLength_in_cgs;
+  units_time    = us->UnitTime_in_cgs;
+    
+    
+
+  printf("          ***************************************\n");
+  printf("          initializing grackle cooling function\n");
+  printf("          \n");
+  
+  printf("          CloudyTable                        = %s\n",cloudytable);
+  printf("          UVbackground                       = %d\n",UVbackground);
+  printf("          GrackleRedshift                    = %g\n",cooling->GrackleRedshift);
+  printf("          GrackleHSShieldingDensityThreshold = %g\n",cooling->GrackleHSShieldingDensityThreshold);
+  
+
+  if(wrap_init_cooling(cloudytable,UVbackground,units_density, units_length, units_time, grackle_chemistry) != 1)
+    {
+      fprintf(stderr,"Error in initialize_chemistry_data.");
+      exit(-1);
+    }
+  
+  printf("          \n");
+  printf("          ***************************************\n");
+
+
+
+    
+}
+
+/**
+ * @brief Prints the properties of the cooling model to stdout.
+ *
+ * @param cooling The properties of the cooling function.
+ */
+static INLINE void cooling_print_backend(
+    const struct cooling_function_data* cooling) {
+
+  message("Cooling function is 'Grackle'.");
+}
+
+#endif /* SWIFT_COOLING_GRACKLE_H */

src/cooling/grackle/cooling_struct.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_COOLING_STRUCT_NONE_H
+#define SWIFT_COOLING_STRUCT_NONE_H
+
+/**
+ * @file src/cooling/none/cooling_struct.h
+ * @brief Empty infrastructure for the cases without cooling function
+ */
+
+/**
+ * @brief Properties of the cooling function.
+ */
+struct cooling_function_data {
+
+  char GrackleCloudyTable[200];
+  int UVbackground;
+  double GrackleRedshift;
+  double GrackleHSShieldingDensityThreshold;
+
+};
+
+/**
+ * @brief Properties of the cooling stored in the particle data
+ */
+struct cooling_xpart_data {
+
+  /*! Energy radiated away by this particle since the start of the run */
+  float radiated_energy;
+
+
+};
+
+#endif /* SWIFT_COOLING_STRUCT_NONE_H */

src/cooling/grackle/grackle_wrapper.c

+/***********************************************************************
+/
+/ Grackle c wrapper
+/
+/
+/ Copyright (c) 2013, Enzo/Grackle Development Team.
+/
+/ Distributed under the terms of the Enzo Public Licence.
+/
+/ The full license is in the file LICENSE, distributed with this 
+/ software.
+************************************************************************/
+
+
+#include <cooling/grackle/grackle_wrapper.h>
+
+#define GRACKLE_DEBUG
+
+#ifdef GRACKLE_DEBUG
+#include <assert.h>
+#define GRACKLE_ASSERT(X) assert( (X) )
+#else
+#define GRACKLE_ASSERT(X) 
+#endif
+
+code_units my_units;
+
+// arrays passed to grackle as input and to be filled
+#define FIELD_SIZE  1
+
+gr_float HDI_density[FIELD_SIZE];
+
+// Set grid dimension and size.
+// grid_start and grid_end are used to ignore ghost zones.
+const int grid_rank = 1;
+
+int wrap_init_cooling(char* CloudyTable, int UVbackground, double udensity, double ulength, double utime, int grackle_chemistry){
+
+  double velocity_units;
+  
+  // First, set up the units system.
+  // These are conversions from code units to cgs.
+
+  my_units.a_units = 1.0; // units for the expansion factor (1/1+zi)
+
+  my_units.comoving_coordinates = 0;		/* so, according to the doc, we assume here all physical quantities to be in proper coordiname (not comobile)  */
+  my_units.density_units = udensity;
+  my_units.length_units = ulength;
+  my_units.time_units = utime;  
+  velocity_units = my_units.a_units * my_units.length_units /my_units.time_units;
+  my_units.velocity_units = velocity_units;  
+
+
+  // Second, create a chemistry object for parameters and rate data.
+  if (set_default_chemistry_parameters() == 0) {
+    fprintf(stderr, "Error in set_default_chemistry_parameters.\n");
+    return 0;
+  }
+
+
+  // Set parameter values for chemistry.
+  grackle_data.use_grackle = 1;          
+  grackle_data.with_radiative_cooling = 1;
+  grackle_data.primordial_chemistry = grackle_chemistry;   // molecular network with H, He, D
+  grackle_data.metal_cooling = 1;          // metal cooling on
+  grackle_data.UVbackground = UVbackground;            
+  grackle_data.grackle_data_file = CloudyTable;
+
+  // Finally, initialize the chemistry object.
+  // snl: a_value is not the true initial ae!! This should get set during update_UVbackground
+  double initial_redshift = 0.;
+  double a_value = 1. / (1. + initial_redshift);
+  
+  // Finally, initialize the chemistry object.
+  if (initialize_chemistry_data(&my_units, a_value) == 0) {
+    fprintf(stderr, "Error in initialize_chemistry_data.\n");
+    return 0;
+  }
+  
+  return 1;
+}
+
+//int wrap_update_UVbackground_rates(double auni) {
+//	// The UV background rates must be updated before 
+//	// calling the other functions.
+//	/* a_value = auni / my_units.a_units; */
+//	double a_value = auni / my_units.a_units;
+//	if (update_UVbackground_rates(my_units, a_value) == 0) {
+//		return 0;
+//	}
+//	return 1;
+//}
+
+int wrap_set_UVbackground_On() {
+	// The UV background rates is enabled
+        grackle_data.UVbackground = 1;
+	return 1;
+}
+
+int wrap_set_UVbackground_Off() {
+	// The UV background rates is disabled
+        grackle_data.UVbackground = 0;
+	return 1;
+}
+
+int wrap_get_cooling_time(double rho, double u,double Z, double a_now, double *coolingtime)
+{
+  gr_float den_factor = 1.0;
+  gr_float u_factor = 1.0;
+  
+  gr_float x_velocity[FIELD_SIZE] = { 0.0 };
+  gr_float y_velocity[FIELD_SIZE] = { 0.0 };
+  gr_float z_velocity[FIELD_SIZE] = { 0.0 };
+  
+  gr_float Density[FIELD_SIZE] = { rho*den_factor };
+  gr_float metal_density[FIELD_SIZE] = { Z * Density[0] };
+  gr_float energy[FIELD_SIZE] = { u * u_factor };
+  
+  gr_float cooling_time[FIELD_SIZE] = { 0.0 };
+
+  
+  int grid_dimension[3] = {1, 0, 0};
+  int grid_start[3] = {0, 0, 0};
+  int grid_end[3] = {0, 0, 0};
+  
+  if(FIELD_SIZE != 1){
+    fprintf(stderr,"field_size must currently be set to 1.\n");
+    return 0;
+  } 
+
+
+  // passed density and energy are proper
+  /*
+  if(my_units.comoving_coordinates){
+    den_factor = pow(a_now, 3);
+    u_factor   = pow(a_now, 0);
+  } else {
+    den_factor = 1.0;
+    u_factor   = 1.0;
+  }
+  */
+      
+  if (calculate_cooling_time_table(&my_units,
+    			     a_now,
+    			     grid_rank, grid_dimension,
+    			     grid_start, grid_end,
+    			     Density, energy,
+    			     x_velocity, y_velocity, z_velocity,
+    			     metal_density, 
+    			     cooling_time) == 0) {
+    fprintf(stderr,"Error in calculate_cooling_time.\n");
+    return 0;
+  }
+      
+  // return updated chemistry and energy
+  for (int i = 0; i < FIELD_SIZE; i++) {
+    *coolingtime = cooling_time[i];
+  }
+
+  return 1;
+}
+
+
+int wrap_do_cooling(double rho, double *u, double dt,double Z, double a_now)
+{
+
+  gr_float den_factor = 1.0;
+  gr_float u_factor = 1.0;
+  
+  gr_float x_velocity[FIELD_SIZE] = { 0.0 };
+  gr_float y_velocity[FIELD_SIZE] = { 0.0 };
+  gr_float z_velocity[FIELD_SIZE] = { 0.0 };
+  
+  
+  gr_float Density[FIELD_SIZE] = { rho*den_factor };
+  gr_float metal_density[FIELD_SIZE] = { Z * Density[0] };
+  gr_float energy[FIELD_SIZE] = { (*u) * u_factor };
+
+  int grid_dimension[3] = {1, 0, 0};
+  int grid_start[3] = {0, 0, 0};
+  int grid_end[3] = {0, 0, 0};
+  
+  GRACKLE_ASSERT(FIELD_SIZE == 1);
+
+  // passed density and energy are proper
+  /*
+  if(my_units.comoving_coordinates){
+    den_factor = pow(a_now, 3);
+    u_factor   = pow(a_now, 0);
+  } else {
+    den_factor = 1.0;
+    u_factor   = 1.0;
+  }
+  */
+      
+#ifdef GRACKLE_DEBUG
+  double old_value = energy[0];
+#endif
+  if (solve_chemistry_table(&my_units,
+  		      a_now, dt,
+  		      grid_rank, grid_dimension,
+  		      grid_start, grid_end,
+  		      Density, energy,
+  		      x_velocity, y_velocity, z_velocity,
+  		      metal_density) == 0) {
+    fprintf(stderr,"Error in solve_chemistry.\n");
+    return 0;
+  }  
+  GRACKLE_ASSERT(old_value != energy[0]);
+
+  // return updated chemistry and energy
+  for (int i = 0; i < FIELD_SIZE; i++) {
+    *u =  energy[i]/u_factor;
+  }
+  
+  return 1;
+}
+
+
+
+

src/cooling/grackle/grackle_wrapper.h

+/***********************************************************************
+/
+/ Grackle c wrapper
+/
+/
+/ Copyright (c) 2013, Enzo/Grackle Development Team.
+/
+/ Distributed under the terms of the Enzo Public Licence.
+/
+/ The full license is in the file LICENSE, distributed with this 
+/ software.
+************************************************************************/
+#ifndef SWIFT_COOLING_GRACKLE_WRAPPER_H
+#define SWIFT_COOLING_GRACKLE_WRAPPER_H
+
+#include <config.h>
+#include <grackle.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>