Matthew Turk avatar Matthew Turk committed 189c02b

Primordial network now only outputs 8 cells, and just the
ionization/recombination for HI. Seems to work to 1e-7. Had to make two
changes:

1) Remove de from the weight calculation
2) Fix the Jacobian's row/column ordering.

Comments (0)

Files changed (4)

dengo/chemical_network.py

     def print_mass_density(self):
         eq = 0
         for s in sorted(self.required_species):
-            if s.name != 'ge': 
+            if s.name not in ('ge', 'de'): 
                 eq += s.symbol
         return ccode(eq)
 

examples/write_oxygen_network.py

 
 if generate_initial_conditions:
     import numpy as na
-    NCELLS = 64
+    NCELLS = 8
     density = 1.0e3
     init_array = na.ones(NCELLS) 
     X = 1.0e-4

examples/write_primordial_network.py

 from dengo.known_species import *
 
 primordial = ChemicalNetwork()
-for ca in cooling_registry.values():
-    #if any(sn.name.startswith("HeI") for sn in
-           #ca.species): continue
-    primordial.add_cooling(ca)
+#for ca in cooling_registry.values():
+#    primordial.add_cooling(ca)
+primordial.add_cooling(cooling_registry["ceHI"])
 
-for i, rname in enumerate(sorted(reaction_registry)):
-    s = reaction_registry[rname]
-    if i == -14:
-        print "rhs_total[0] skipping", rname, s
-        #continue
-    if any(sn in ("H2II","HM") for sn in 
-           s.considered): continue
-    primordial.add_reaction(s)
+#for i, rname in enumerate(sorted(reaction_registry)):
+#    s = reaction_registry[rname]
+#    primordial.add_reaction(s)
+s = reaction_registry["k01"]
+primordial.add_reaction(s)
+s = reaction_registry["k02"]
+primordial.add_reaction(s)
 
 primordial.init_temperature((1e0, 1e8))
 
 
 if generate_initial_conditions:
     import numpy as na
-    NCELLS = 64
+    NCELLS = 8
     density = 1.0e3
     init_array = na.ones(NCELLS) 
-    X = 1.0e-4
+    X = 0.5
 
     init_values = dict()
     init_values['density'] = density * init_array

templates/rates_and_rate_tables.c.template

     H5LTread_dataset_double(file_id, "/{{ s.name }}", tics);
     for (j = 0; j < dims; j++) {
         ics[j * N + i] = tics[j];
-        atol[j * N + i] = tics[j] * 1e-1 * 0.0;
-        rtol[j * N + i] = 1e-1;
+        atol[j * N + i] = tics[j] * 1e-4;
+        rtol[j * N + i] = 1e-7;
         if(j==0) {
             fprintf(stderr, "{{s.name}}[0] = %0.3g, atol => % 0.16g\n",
                     tics[j], atol[j]);
 
     H5Fclose(file_id);
 
-    float dt = 1e5;
+    float dt = 1e10;
     rhs_f f = calculate_rhs_{{solver_name}};
     jac_f jf = calculate_jacobian_{{solver_name}};
     int rv = BE_chem_solve(f, jf, ics, dt, rtol, atol, dims, N, (void *) data);
         // Species: {{s1.name }}
         // 
         {%- for s2 in network.required_species | sort %}
-            // {{s1.name}} by {{s2.name}}
-            {% if s1.name == s2.name %}
-            Joutput[j] = 1.0;
-            {% else %}
-            {{ network.print_jacobian_component(s1, s2, assign_to="Joutput[j]") }}
-            {% endif %}
-	        {% if s2.name == 'ge' %}
+            // {{s2.name}} by {{s1.name}}
+            {{ network.print_jacobian_component(s2, s1, assign_to="Joutput[j]") }}
+	        {% if s1.name == 'ge' %}
             Joutput[j] *= T{{ network.energy_term.name }}[i];
             {% endif %}
             if(Joutput[j] != -100.0) {
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.