Fixed typo in cooling.c which was causing SWIFT to crash

examples/UniformBox/makeIC.py

@@ -27,10 +27,10 @@ from numpy import *
 
 # Parameters
 periodic= 1           # 1 For periodic box
-boxSize = 1.0         
+boxSize = 1.
 L = int(sys.argv[1])  # Number of particles along one axis
-rho = 3.3e6           # Density in solar masses per cubic kiloparsec (0.1 hydrogen atoms per cm^3)
-P = 1.4e8             # Pressure in code units (at 10^5K)
+rho = 2.              # Density
+P = 1.                # Pressure
 gamma = 5./3.         # Gas adiabatic index
 eta = 1.2349          # 48 ngbs with cubic spline kernel
 fileName = "uniformBox.hdf5" 
@@ -38,7 +38,6 @@ fileName = "uniformBox.hdf5"
 #---------------------------------------------------
 numPart = L**3
 mass = boxSize**3 * rho / numPart
-print mass
 internalEnergy = P / ((gamma - 1.)*rho)
 
 #--------------------------------------------------

examples/UniformBox/plots/energy_plot.py

-import numpy as np
-import matplotlib.pyplot as plt
-import h5py as h5
-
-n_snaps = 101
-directory = '../'
-snap_name_base = "uniformBox_"
-plot_dir = "./"
-cooling_lambda = 1.0
-u_floor = 0.0
-delta_time = 0.1
-plot_title = r"$\Lambda \, = \, %1.3f \, \, \, u_{floor} = %1.3f$" %(cooling_lambda,u_floor)
-plot_filename = "energy_plot_creasey_lambda_1p0eminus50.png"
-e_kin_array=[]
-e_therm_array=[]
-e_total_array = []
-time = np.zeros(n_snaps-1)
-analytic_solution = np.zeros(n_snaps-1)
-for i in range(1,n_snaps):
-    snap_number = "%03d" %i
-    filename = directory + snap_name_base + snap_number + ".hdf5"
-    f = h5.File(filename)
-    print "Reading snap number %d" %i
-    header = f["Header"]
-    n_parts = header.attrs["NumPart_ThisFile"][0]
-    time[i-1] = i*delta_time
-    u = np.array(f["PartType0/InternalEnergy"])
-    u_total = np.sum(u)
-    e_therm_array = np.append(e_therm_array,u_total)
-    v = np.array(f["PartType0/Velocities"])
-    e_kin_particles = 0.5*v**2
-    e_kin_total = np.sum(e_kin_particles)
-    e_kin_array = np.append(e_kin_array,e_kin_total)
-    e_total_array = np.append(e_total_array,e_kin_total+u_total)
-    analytic_solution[i-1] = e_total_array[0] - n_parts*(time[i-1] - time[0])*cooling_lambda
-
-
-# plt.plot(e_kin_array,label = "Kinetic Energy")
-# plt.plot(e_therm_array,label = "Internal Energy")
-plt.plot(time,e_total_array,'r',label = "Total Energy")
-plt.plot(time,analytic_solution,'--k',label = "Analytic Solution")
-plt.xlabel("Time (seconds)")
-plt.ylabel("Energy (erg/g)")
-plt.ylim((0,1000))
-plt.legend(loc = "upper right")    
-plt.title(plot_title)
-full_plot_filename = "%s/%s" %(plot_dir,plot_filename)
-plt.savefig(full_plot_filename,format = "png")
-plt.close()

examples/UniformBox/plots/plot_particle_positions.py

-import h5py as h5
-import numpy as np
-import matplotlib.pyplot as plt
-
-n_snaps = 100
-directory = '../'
-snap_name_base = "uniformBox_"
-plot_directory = "./"
-plot_filename_base = "snap"
-for i in range(n_snaps):
-    snap_number = "%03d" %i
-    filename = directory + snap_name_base + snap_number + ".hdf5"
-    f = h5.File(filename)
-    #find the box size
-
-    header = f["Header"]
-    box_size = header.attrs["BoxSize"][0]
-    pos = np.array(f["PartType0/Coordinates"])
-    x = pos[:,0]
-    y = pos[:,1]
-    z = pos[:,2]
-
-    #x_projection
-    plt.plot(y,z,'bo')
-    plt.xlabel("y")
-    plt.ylabel("z")
-    plt.xlim((0.0,box_size))
-    plt.ylim((0.0,box_size))
-    plot_dir = "%s/projection_x" %plot_directory
-    plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
-    plt.savefig(plot_filename,format = "png")
-    plt.close()
-
-    #y_projection
-    plt.plot(x,z,'bo')
-    plt.xlabel("x")
-    plt.ylabel("z")
-    plt.xlim((0.0,box_size))
-    plt.ylim((0.0,box_size))
-    plot_dir = "%s/projection_y" %plot_directory
-    plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
-    plt.savefig(plot_filename,format = "png")
-    plt.close()
-
-    #x_projection
-    plt.plot(x,y,'bo')
-    plt.xlabel("x")
-    plt.ylabel("y")
-    plt.xlim((0.0,box_size))
-    plt.ylim((0.0,box_size))
-    plot_dir = "%s/projection_z" %plot_directory
-    plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
-    plt.savefig(plot_filename,format = "png")
-    plt.close()
-    
-

examples/UniformBox/run.sh

@@ -4,7 +4,7 @@
 if [ ! -e uniformBox.hdf5 ]
 then
     echo "Generating initial conditions for the uniform box example..."
-    python makeIC.py 10
+    python makeIC.py 100
 fi
 
-../swift -s -C -t 16 uniformBox.yml
+../swift -s -t 16 uniformBox.yml

examples/UniformBox/uniformBox.yml

 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     2.0e33   # Solar masses
-  UnitLength_in_cgs:   3.01e21   # Kilparsecs
-  UnitVelocity_in_cgs: 1.0e5   # Kilometres per second
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
 
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   10.    # The end time of the simulation (in internal units).
+  time_end:   1.    # The end time of the simulation (in internal units).
   dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
 
@@ -17,7 +17,7 @@ TimeIntegration:
 Snapshots:
   basename:            uniformBox # Common part of the name of output files
   time_first:          0.         # Time of the first output (in internal units)
-  delta_time:          0.1       # Time difference between consecutive outputs (in internal units)
+  delta_time:          0.01       # Time difference between consecutive outputs (in internal units)
 
 # Parameters governing the conserved quantities statistics
 Statistics:
@@ -41,12 +41,3 @@ PointMass:
   position_y:      50.
   position_z:      50.	
   mass:            1e10     # mass of external point mass in internal units
-
-  # Cooling parameters (Creasey cooling) (always in cgs)
-
-CreaseyCooling:
-  Lambda:                      1.0e-22
-  minimum_temperature:         1.0e4
-  mean_molecular_weight:       0.59
-  hydrogen_mass_abundance:     0.75
-  cooling_tstep_mult:          1.0
\ No newline at end of file

src/cooling.c

@@ -133,7 +133,7 @@ float calculate_new_thermal_energy(float u_old, float rho, float dt,
   float u_new;
   float m_p = phys_const->const_proton_mass;
   float X_H = cooling->creasey_cooling.hydrogen_mass_abundance;
-  float lambda = cooling->creasey_cooling.lambda; //this is always in cgs
+  float lambda_cgs = cooling->creasey_cooling.lambda; //this is always in cgs
   float u_floor = cooling->creasey_cooling.min_internal_energy;
 
   /*convert from internal code units to cgs*/
@@ -143,19 +143,19 @@ float calculate_new_thermal_energy(float u_old, float rho, float dt,
   float n_H_cgs = rho_cgs / m_p_cgs;
   float u_old_cgs =  u_old * units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS);
   float u_floor_cgs =  u_floor * units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS);
-  float du_dt = -lambda * n_H_cgs * n_H_cgs / rho;
+  float du_dt_cgs = -lambda_cgs * n_H_cgs * n_H_cgs / rho_cgs;
   float u_new_cgs;
 
-  if (u_old_cgs - du_dt*dt_cgs > u_floor_cgs){
-    u_new_cgs = u_old_cgs + du_dt*dt_cgs;
+  if (u_old_cgs + du_dt_cgs*dt_cgs > u_floor_cgs){
+    u_new_cgs = u_old_cgs + du_dt_cgs*dt_cgs;
   }
   else{
     u_new_cgs = u_floor_cgs;
   }
-
   /*convert back to internal code units when returning new internal energy*/
 
   u_new = u_new_cgs / units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS);  
+  
 #endif /*CREASEY_COOLING*/
   return u_new;
 }