Commits

Britton Smith  committed f1ae283

Cleaning up example scripts and setting them both up to produce plots.

  • Participants
  • Parent commits f09d12e

Comments (0)

Files changed (2)

File src/python/examples/cooling_rate.py

 
 import copy
 
-def check_convergence(fc1, fc2, fields=None, tol=0.1):
+def check_convergence(fc1, fc2, fields=None, tol=0.10):
     if fields is None:
         fields = ["HI", "HII", "HM", "HeI", "HeII", "HeIII",
                   "H2I", "H2II", "DI", "DII", "HDI", "de"]
     for field in fields:
-        if (np.abs(fc1[field] - fc2[field]) / fc1[field] > tol).any():
+        convergence = np.abs(fc1[field] - fc2[field]) / fc1[field]
+        if (convergence > tol).any():
+            print "Max change %s: %e." % (field, convergence.max())
             return False
     return True
 
 hplanck     = 6.6260693e-27
 ev2erg      = 1.60217653e-12
 c_light     = 2.99792458e10
-GravConst   = 6.67428e-8
+GravConst   = 6.6726e-8
 sigma_sb    = 5.670373e-5
 SolarMass   = 1.9891e33
 Mpc         = 3.0857e24
 kpc         = 3.0857e21
 pc          = 3.0857e18
+yr_to_s     = 3.15569e7
 
 from pygrackle.grackle_wrapper import *
 from pygrackle.fluid_container import FluidContainer
 my_chemistry.comoving_coordinates = 0
 my_chemistry.density_units = 1.67e-24
 my_chemistry.length_units = 1.0e16
-my_chemistry.time_units = 3.15569e7 * 1e6
+my_chemistry.time_units = yr_to_s * 1e6
 my_chemistry.a_units = 1.0
 
 energy_units = (my_chemistry.length_units /
 dt = 0.1
 my_time = 0.0
 while True:
-    print "t = %e." % my_time
+    print "t = %.2f Myr." % (my_time * my_chemistry.time_units /
+                           (yr_to_s * 1e6))
     for field in ["HI", "HII", "HM", "HeI", "HeII", "HeIII",
                   "H2I", "H2II", "DI", "DII", "HDI", "de"]:
         fc_last[field] = np.copy(fc[field])
 
 calculate_temperature(fc)
 calculate_cooling_time(fc, a_value)
-cooling_rate = fc["energy"] / fc["cooling_time"]
+
+cooling_rate = fc["density"] * my_chemistry.density_units * \
+  fc["energy"] * energy_units / \
+  (fc["cooling_time"] * my_chemistry.time_units)
 
 from matplotlib import pyplot
 pyplot.loglog(fc['temperature'], cooling_rate)
-pyplot.savefig('cooling_rate.png')
+pyplot.xlabel('T [K]')
+pyplot.ylabel('$\\Lambda$ [erg s$^{-1}$ cm$^{-3}$]')
+output_file = 'cooling_rate.png'
+print "Writing %s." % output_file
+pyplot.savefig(output_file)

File src/python/examples/freefall.py

 hplanck     = 6.6260693e-27
 ev2erg      = 1.60217653e-12
 c_light     = 2.99792458e10
-GravConst   = 6.67428e-8
+GravConst   = 6.6726e-8
 sigma_sb    = 5.670373e-5
 SolarMass   = 1.9891e33
 Mpc         = 3.0857e24
 kpc         = 3.0857e21
 pc          = 3.0857e18
+yr_to_s     = 3.15569e7
 
 from pygrackle.grackle_wrapper import *
 from pygrackle.fluid_container import FluidContainer
 energy_units = (my_chemistry.length_units /
                 my_chemistry.time_units)**2.0
 
-gravitational_constant = (4.0 * 3.1415926 * 6.6726e-8 * 
+gravitational_constant = (4.0 * pi_val * GravConst * 
   my_chemistry.density_units * my_chemistry.time_units**2)
 
 a_value = 1.0
 
 my_chemistry.initialize(a_value)
 
-
 my_chemistry.UVbackground = 1;
 my_chemistry.UVbackground_file = "UVB_rates_HM2012.hdf5";
 my_chemistry.initialize_UVbackground()
 my_chemistry.update_UVbackground(a_value)
 
-
 fc = FluidContainer(my_chemistry, 1)
 fc["density"][:] = 1.0
 fc["HI"][:] = 0.76 * fc["density"]
 fc["y-velocity"][:] = 0.0
 fc["z-velocity"][:] = 0.0
 
+density_values = []
+temperature_values = []
+
 timestep_fraction = 0.1
 current_time = 0.0
 while fc["density"][0] < 1.e10:
     dt = timestep_fraction * \
       np.power(((3. * np.pi) / 
-                (32. * gravitational_constant * fc["density"][0])), 0.5)
+                (32. * gravitational_constant * 
+                 fc["density"][0])), 0.5)
 
     calculate_temperature(fc)
-    print fc["density"][0], fc["temperature"][0]
+
+    density_values.append(fc["density"][0] * my_chemistry.density_units)
+    temperature_values.append(fc["temperature"][0])
+
+    print "t: %e yr, rho: %e g/cm^3, T: %e [K]" % \
+      ((current_time * my_chemistry.time_units / yr_to_s),
+       (fc["density"][0] * my_chemistry.density_units),
+       fc["temperature"][0])
     
     density_ratio = np.power((freefall_constant - 
-                              (0.5 * freefall_time_constant * current_time)), -2.) / \
+                              (0.5 * freefall_time_constant * 
+                               current_time)), -2.) / \
                               fc["density"][0]
                               
     for field in fc.density_fields:
       freefall_time_constant * np.power(fc["density"][0], 0.5) * dt
 
     solve_chemistry(fc, a_value, dt)
+    
+    current_time += dt
 
-    current_time += dt
+from matplotlib import pyplot
+pyplot.loglog(density_values, temperature_values)
+pyplot.xlabel('$\\rho$ [g cm$^{-3}$]')
+pyplot.ylabel('T [K]')
+output_file = 'freefall.png'
+print "Writing %s." % output_file
+pyplot.savefig(output_file)