Commits

Devin Silvia  committed ffd63cb

minor modifications and changes to the plotting script to produce relative abundances.

  • Participants
  • Parent commits 42c5a48

Comments (0)

Files changed (3)

File examples/write_primordial_network.py

 primordial = ChemicalNetwork()
 primordial.add_energy_term()
 primordial.write_intermediate_solutions = True
-#for ca in cooling_registry.values():
-    #if not all(sp.name in want for sp in ca.species): continue
-#    primordial.add_cooling(ca)
+for ca in cooling_registry.values():
+    if not all(sp.name in want for sp in ca.species): continue
+    primordial.add_cooling(ca)
 #primordial.add_cooling(cooling_registry["ceHI"])
 
 for i, rname in enumerate(sorted(reaction_registry)):
     s = reaction_registry[rname]
-    #if not all(sp.name in want for sp in s.species): continue
+    if not all(sp.name in want for sp in s.species): continue
     primordial.add_reaction(s)
 #s = reaction_registry["k01"]
 #primordial.add_reaction(s)
 if generate_initial_conditions:
     import numpy as na
     NCELLS = 8
-    density = 1.0e15
+    density = 1.0
     init_array = na.ones(NCELLS) 
     X = 0.5
 
     
     # set up initial temperatures values used to define ge
     temperature = na.logspace(2, 4, NCELLS)
-    temperature[:] = 1e3
+    temperature[:] = 2e4
 
     # calculate the number density
     number_density = 0.0

File plot_intermediate.py

                 else:
                     plt.loglog(self.t, self.data[s][0,:], '-', label=s, lw=1.5)
         plt.xlabel("Time")
-        plt.ylim(1e-20, 1e16)
+        plt.ylim(1e-20, 1e10)
         leg = plt.legend(loc = "best")
-        leg.set_alpha(0.1)
+        leg.legendPatch.set_alpha(0.1)
         plt.savefig("%s_time.png" % (self.network_name))
+        mpl.rcParams['axes.color_cycle'] = [list(clr) for clr in mpl.cm.spectral(np.linspace(0,1,(self.nspecies)))]
+        plt.clf()
+        for s in sorted(self.species):
+            #if s != 'ge':
+            if s == 'HII':
+                plt.loglog(self.t, self.data[s][0,:]/self.data[s][0,0], '-', label=s, lw=1.5, marker='x')
+            else:
+                plt.loglog(self.t, self.data[s][0,:]/self.data[s][0,0], '-', label=s, lw=1.5)
+        plt.xlabel("Time")
+        plt.ylim(1e-4, 1e1)
+        leg = plt.legend(loc = "best")
+        leg.legendPatch.set_alpha(0.1)
+        plt.savefig("%s_time_relative.png" % (self.network_name))
         
         # plot the initial conditions and the final solution
         mpl.rcParams['axes.color_cycle'] = [list(clr) for clr in mpl.cm.spectral(np.linspace(0,1,(self.nspecies - 2)))]
         leg = plt.legend(loc = "best")
         leg.set_alpha(0.1)
         plt.savefig("%s_ic_final.png" %(self.network_name))
+        
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser()

File templates/rates_and_rate_tables.c.template

     rhs_f f = calculate_rhs_{{solver_name}};
     jac_f jf = calculate_jacobian_{{solver_name}};
     float dtf = 3.1557e16;
-    float dt = 1.0; //dtf / 10000.0;
+    float dt = 1.e3; //dtf / 10000.0;
     int niter = 0;
     int siter = 0;
     double ttot = 0;
     double *scale = (double *) alloca(dims * N * sizeof(double));
-    for (i = 0; i < dims * N; i++) scale[i] = 1.0; //+ 10.0*i/dims/N;
+    for (i = 0; i < dims * N; i++) scale[i] = 1.0; //ics[i];
     for (i = 0; i < dims * N; i++) input[i] = prev[i] = ics[i];
     while (ttot < dtf) {
         int rv = BE_chem_solve(f, jf, input, dt, rtol, atol, dims, N, scale, (void *) data);