Commits

Matthew Turk committed 22e97fa

A few cleanups and messups of primordial, to pass off to Dan.

Comments (0)

Files changed (3)

examples/write_primordial_network.py

 from dengo.chemistry_constants import tiny, kboltz, mh
 from dengo.known_species import *
 
+want = ("HI", "HII", "de", "ge")
+
 primordial = ChemicalNetwork()
+primordial.add_energy_term()
+primordial.write_intermediate_solutions = True
 for ca in cooling_registry.values():
-    if any(sp.name in ("H2I", "H2II", "HeI", "HeII", "HeIII") for sp in ca.species): continue
+    #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 any(sp.name in ("H2I", "H2II", "HeI", "HeII", "HeIII") 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)
     
     # set up initial temperatures values used to define ge
     temperature = na.logspace(2, 4, NCELLS)
+    temperature[:] = 2000000
 
     # calculate the number density
     number_density = 0.0

plot_intermediate.py

 
     def grab_results(self):
         # grab initial conditions
-        icfn = "output/%s_initial_conditions.h5" %(self.network_name)
+        icfn = "%s_initial_conditions.h5" %(self.network_name)
         icf = h5py.File(icfn)
         self.icspecies = icf.keys()
         self.icnspecies = len(self.icspecies)
         icf.close()
 
         # grab intermediate data
-        fns = glob.glob("output/%s_intermediate_*.h5" % (self.network_name))
+        fns = glob.glob("%s_intermediate_*.h5" % (self.network_name))
         if len(fns) == 0: raise RuntimeError
         f = h5py.File(fns[0])
         self.species = f.keys()
 
         # grab the final solution
                 # grab initial conditions
-        fsfn = "output/%s_solution.h5" %(self.network_name)
+        fsfn = "%s_solution.h5" %(self.network_name)
         fsf = h5py.File(fsfn)
         self.fsspecies = fsf.keys()
         self.fsnspecies = len(self.fsspecies)
         plt.xlabel("Time")
         plt.ylim(1e-20, 1e8)
         plt.legend(loc = "best")
-        plt.savefig("output/%s_time.png" % (self.network_name))
+        plt.savefig("%s_time.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)))]
         plt.xlabel("Temperature")
         plt.ylim(1e-3, 10)
         plt.legend(loc = "best")
-        plt.savefig("output/%s_ic_final.png" %(self.network_name))
+        plt.savefig("%s_ic_final.png" %(self.network_name))
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser()

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 = dtf;
+    float dt = dtf / 10.0;
     int niter = 0;
     int siter = 0;
     double ttot = 0;
     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, (void *) data);
+        /*
         fprintf(stderr, "Return value [%d]: %i.  %0.5g / %0.5g = %0.5g (%0.5g)\n",
                 niter, rv, ttot, dtf, ttot/dtf, dt);
         fprintf(stderr, "Value[80] = %0.5g %0.5g %0.5g\n",
                 input[80], prev[80], ics[80]);
+        */
         for (i = 0; i < dims * N; i++) {
             if (input[i] < 0) {
                 rv = 1;
         if (rv == 0) {
 	    if (siter == 9999) break;
 	    siter++;
+            fprintf(stderr, "Successful Iteration[%d]: %0.16g / %0.16g\n",
+                     siter, dt, dtf);
             ttot += dt;
             dt = min(dt * 2.0, dtf - ttot);
             for (i = 0; i < dims * N; i++) prev[i] = input[i];
     	    for (j = 0; j < dims; j++) {
             	{{ s.name }}[j] = input[j * N + i];
     	    }
-    	    fprintf(stderr, "Writing solution for /{{ s.name }}\n");
+    	    //fprintf(stderr, "Writing solution for /{{ s.name }}\n");
     	    H5LTmake_dataset_double(file_id, "/{{ s.name }}", 1, dimsarr, {{ s.name }});
     	    i++;
     	    {% endfor %}
         {%if species.name != "de" %}
           total+={{species.name}} * {{species.weight}};
         {%endif%}
-        /*fprintf(stderr, "Density[%d][{{species.name}}] = % 0.16g\n",
-                i, {{species.name}} * {{species.weight}});*/
         {%endif%}
         j++;
     {% endfor %}
-        /*fprintf(stderr, "Density[%d] = % 0.16g\n", i, total);*/
         mdensity = total * mh;
         total = 0.0;
         j = i * nchem;
         // Species: {{species.name}}
         // 
         {{ network.print_ccode(species, assign_to = "rhs[j]") }}
-        /*fprintf(stderr, "rhs[{{species.name}}][%d] = % 0.16g\n", i,
-        rhs[j]);*/
         {% if species.name != "ge" and species.name != "de" %}
             /* Already in number density, not mass density */
             total += rhs[j] * {{species.weight}};
             total_e += {{species.name}} * {{species.free_electrons}};
         {% elif species.name == "ge" %}
-            rhs[j] = 0.0;
 	    rhs[j] /= mdensity;
         {% endif %}
         {% if species.name == "de" %}
         {%endif %}
         j++;
     {% endfor %}
-        /*fprintf(stderr, "rhs_total[%d] = % 0.16g total_e = % 0.16g total_de = %0.16g\n",
-                i, total, total_e - de, total_de);
-                */
     }  
     return 0;
 }
         j = i * nchem;
 	total = mdensity = 0.0;
     {%- for species in network.required_species | sort %}
-    	{% if species.name != 'ge' %}
-        {{species.name}} = input[j];
-	{% else %}
-	{{species.name}} = input[j];
-	{% endif %}
+	    {{species.name}} = input[j];
         if ({{species.name}} < 0.0) {
           fprintf(stderr, "JNegative[%d][{{species.name}}] = % 0.16g [%d]\n",
             i, {{species.name}}, j);
 	    {% if s1.name == 'ge' %}
             Joutput[j] *= T{{ network.energy_term.name }}[i];
             {% endif %}
-            if(Joutput[j] != -100.0) {
-	    	{% if s2.name == 'geaaaa' %}
-                fprintf(stderr, "Joutput[{{s1.name}}][{{s2.name}}][%d] = % 0.16g\n",
-                    i, Joutput[j]);
-		{% else %}
-		/*fprintf(stderr, "Joutput[{{s1.name}}][{{s2.name}}][%d] = % 0.16g, T{{ network.energy_term.name }}[%d] = % 0.16g\n",
-                    i, Joutput[j], i, Tge[i]);
-                    */
-		{% endif %}
-            }
             j++;
         {%- endfor %}
     {% endfor %}
     }
-    int nzero = 0;
-    double mi = 1e30;
-    double ma = -1e30;
-    double avg = 0.0;
-    for (i = 0; i < nchem * nchem * nstrip; i++) {
-        if (Joutput[i] == 0.0) nzero++;
-        if (Joutput[i] < mi) mi = Joutput[i];
-        if (Joutput[i] > ma) ma = Joutput[i];
-        avg += Joutput[i];
-    }
-    /*fprintf(stderr, "JOUTPUTSTATS %d % 0.16g % 0.16g %d % 0.16g\n",
-        nchem*nchem*nstrip, mi, ma, nzero, avg / (nchem*nchem*nstrip));
-        */
 
     return 0;