Commits

Devin Silvia committed 579f5b9

instead of a linear extrapolation beyond cooling rate
tables values, the code now does a gaussian fall off
to "tiny".

Also, cooling is not turned back on for the oxygen
network as it no longer stalls out when it hits the
temperature "floor"

Comments (0)

Files changed (2)

dengo/reaction_classes.py

 
     def cooling_rate(network):
         # Read in cooling rates from Gnat & Ferland 2012
-        # and do linear interpolation with extrapolation on the ends
+        # and do linear interpolation and then recompute
+        # the ends with either an extrapolation or falloff
         f = h5py.File('dengo/%s_ion_by_ion_cooling.h5' %(element_name))
         data = f['Table']
         
-        # interpolation
+        ### Intepolate values within table values ###
         vals = np.interp(network.T, data['T'], data['%s' %(ion_name)])
         
-        # extrapolation in logspace
-        vals = np.log10(vals)
-        logT = np.log10(network.T)
-        logdataT = np.log10(data['T'])
-        logdataS = np.log10(data['%s' %(ion_name)])
-        extrapdown = logdataS[0] + \
-            (logT - logdataT[0]) * (logdataS[0] - logdataS[1]) \
-            / (logdataT[0] - logdataT[1])
-        vals[logT < logdataT[0]] = extrapdown[logT < logdataT[0]]
-        extrapup = logdataS[-1] + \
-            (logT - logdataT[-1]) * (logdataS[-1] - logdataS[-2]) \
-            / (logdataT[-1] - logdataT[-2])
-        vals[logT > logdataT[-1]] = extrapdown[logT > logdataT[-1]]
+        end_method = 1 # 0 = extrapolation, 1 = gaussian falloff
 
-        # vals[logT < logdataT[0]] = np.log10(tiny)
-        # vals[logT > logdataT[-1]] = np.log10(tiny)
+        if end_method == 0:
+            ### Extrapolation in logspace ###
+            # convert to log space
+            vals = np.log10(vals)
+            logT = np.log10(network.T)
+            logdataT = np.log10(data['T'])
+            logdataS = np.log10(data['%s' %(ion_name)])
 
-        vals = 10.0**vals
-        
+            # extrapolate
+            extrapdown = logdataS[0] + \
+                (logT - logdataT[0]) * (logdataS[0] - logdataS[1]) \
+                / (logdataT[0] - logdataT[1])
+            vals[logT < logdataT[0]] = extrapdown[logT < logdataT[0]]
+            extrapup = logdataS[-1] + \
+                (logT - logdataT[-1]) * (logdataS[-1] - logdataS[-2]) \
+                / (logdataT[-1] - logdataT[-2])
+            vals[logT > logdataT[-1]] = extrapdown[logT > logdataT[-1]]
+
+            # convert back to linear
+            vals = 10.0**vals
+                
+        if end_method == 1:
+            ### Gaussian falloff when values extend beyond table values ###
+            # rename some variables to symplify code
+            T = network.T
+            dataT = data['T']
+            dataS = data['%s' %(ion_name)]
+
+            # compute gaussian tails
+            gaussdown = dataS[0] * (tiny/dataS[0])**(((T - dataT[0])/(T[0] - dataT[0])))**2
+            vals[T < dataT[0]] = gaussdown[T < dataT[0]]
+            gaussup = dataS[-1] * (tiny/dataS[-1])**(((T - dataT[-1])/(T[-1] - dataT[-1])))**2
+            vals[T > dataT[-1]] = gaussup[T > dataT[-1]]
+
         f.close()
         return vals
 

examples/run_oxygen_network.py

 oxygen = ChemicalNetwork()
 oxygen.add_energy_term()
 
-# for ca in cooling_registry.values():
-#     # The following line can be used to specify a subset of species
-#     #if all(sp.name in want for sp in ca.species):
-#     if ca.name.startswith("O"):
-#        oxygen.add_cooling(ca)
+for ca in cooling_registry.values():
+    # The following line can be used to specify a subset of species
+    #if all(sp.name in want for sp in ca.species):
+    if ca.name.startswith("O"):
+       oxygen.add_cooling(ca)
 
 for s in reaction_registry.values():
     # The following line can be used to specify a subset of species
         oxygen.add_reaction(s)
 
 # This defines the temperature range for the rate tables
-oxygen.init_temperature((1e4, 1e8))
+oxygen.init_temperature((1e0, 1e8))
 
 # Set to false if you don't want intermediate solution output
 oxygen.write_intermediate_solutions = True