Commits

Devin Silvia  committed 9df37a4

added more metals and a write_combined_network.py file.

the combined network creates an atomic only network with
all of the species currently included in Dengo.

  • Participants
  • Parent commits 2ed2fcf

Comments (0)

Files changed (16)

File dengo/magnesium_cooling.py

+"""
+Author: Devin Silvia <devin.silvia@gmail.com>
+Affiliation: UC Boulder
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2012 Matthew Turk.  All Rights Reserved.
+
+  This file is part of the dengo package.
+
+  This file is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from reaction_classes import Species, ion_cooling_rate, species_registry
+import docutils.utils.roman as roman
+
+# Note: the atomic/species properties have to be hard-coded for this
+# we may want to come up with a better solution here...
+atomicSymbol = 'Mg'
+atomicNumber = 12
+atomicWeight = 24
+nIons = atomicNumber + 1
+
+for i in range(nIons):
+    ion_state = i + 1
+    speciesName = "%s%s" %(atomicSymbol, roman.toRoman(ion_state))
+    # Check if the species already exists
+    # in the species registry, if it does
+    # we don't want to create it again
+    if (speciesName in species_registry) == False:
+        s = Species(speciesName, atomicNumber, atomicWeight, i)
+    else:
+        s = species_registry[speciesName]
+    ion_cooling_rate(s)

File dengo/magnesium_rates.py

+"""
+Author: Devin Silvia <devin.silvia@gmail.com>
+Affiliation: UC Boulder
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2012 Matthew Turk.  All Rights Reserved.
+
+  This file is part of the dengo package.
+
+  This file is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from reaction_classes import Species, chianti_rate, species_registry
+import docutils.utils.roman as roman
+
+# Note: the atomic/species properties have to be hard-coded for this
+# we may want to come up with a better solution here...
+atomicSymbol = 'Mg'
+atomicNumber = 12
+atomicWeight = 24
+nIons = atomicNumber + 1
+
+for i in range(nIons):
+    ion_state = i + 1
+    speciesName = "%s%s" %(atomicSymbol, roman.toRoman(ion_state))
+    # Check if the species already exists
+    # in the species registry, if it does
+    # we don't want to create it again
+    if (speciesName in species_registry) == False:
+        s = Species(speciesName, atomicNumber, atomicWeight, i)
+    else:
+        s = species_registry[speciesName]
+
+    if ion_state != nIons:
+        # we need to do this to make sure the 'ion_state + 1' species
+        # exists when chianti_rate is called
+        speciesNamePlusOne = "%s%s" % (atomicSymbol, roman.toRoman(ion_state+1))
+        if (speciesNamePlusOne in species_registry) == False:
+            splusone = Species(speciesNamePlusOne, atomicNumber, atomicWeight, i+1)
+
+    chianti_rate(s)

File dengo/mg_ion_by_ion_cooling.h5

Binary file added.

File dengo/n_ion_by_ion_cooling.h5

Binary file added.

File dengo/ne_ion_by_ion_cooling.h5

Binary file added.

File dengo/neon_cooling.py

+"""
+Author: Devin Silvia <devin.silvia@gmail.com>
+Affiliation: UC Boulder
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2012 Matthew Turk.  All Rights Reserved.
+
+  This file is part of the dengo package.
+
+  This file is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from reaction_classes import Species, ion_cooling_rate, species_registry
+import docutils.utils.roman as roman
+
+# Note: the atomic/species properties have to be hard-coded for this
+# we may want to come up with a better solution here...
+atomicSymbol = 'Ne'
+atomicNumber = 10
+atomicWeight = 20
+nIons = atomicNumber + 1
+
+for i in range(nIons):
+    ion_state = i + 1
+    speciesName = "%s%s" %(atomicSymbol, roman.toRoman(ion_state))
+    # Check if the species already exists
+    # in the species registry, if it does
+    # we don't want to create it again
+    if (speciesName in species_registry) == False:
+        s = Species(speciesName, atomicNumber, atomicWeight, i)
+    else:
+        s = species_registry[speciesName]
+    ion_cooling_rate(s)

File dengo/neon_rates.py

+"""
+Author: Devin Silvia <devin.silvia@gmail.com>
+Affiliation: UC Boulder
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2012 Matthew Turk.  All Rights Reserved.
+
+  This file is part of the dengo package.
+
+  This file is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from reaction_classes import Species, chianti_rate, species_registry
+import docutils.utils.roman as roman
+
+# Note: the atomic/species properties have to be hard-coded for this
+# we may want to come up with a better solution here...
+atomicSymbol = 'Ne'
+atomicNumber = 10
+atomicWeight = 20
+nIons = atomicNumber + 1
+
+for i in range(nIons):
+    ion_state = i + 1
+    speciesName = "%s%s" %(atomicSymbol, roman.toRoman(ion_state))
+    # Check if the species already exists
+    # in the species registry, if it does
+    # we don't want to create it again
+    if (speciesName in species_registry) == False:
+        s = Species(speciesName, atomicNumber, atomicWeight, i)
+    else:
+        s = species_registry[speciesName]
+
+    if ion_state != nIons:
+        # we need to do this to make sure the 'ion_state + 1' species
+        # exists when chianti_rate is called
+        speciesNamePlusOne = "%s%s" % (atomicSymbol, roman.toRoman(ion_state+1))
+        if (speciesNamePlusOne in species_registry) == False:
+            splusone = Species(speciesNamePlusOne, atomicNumber, atomicWeight, i+1)
+
+    chianti_rate(s)

File dengo/nitrogen_cooling.py

+"""
+Author: Devin Silvia <devin.silvia@gmail.com>
+Affiliation: UC Boulder
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2012 Matthew Turk.  All Rights Reserved.
+
+  This file is part of the dengo package.
+
+  This file is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from reaction_classes import Species, ion_cooling_rate, species_registry
+import docutils.utils.roman as roman
+
+# Note: the atomic/species properties have to be hard-coded for this
+# we may want to come up with a better solution here...
+atomicSymbol = 'N'
+atomicNumber = 7
+atomicWeight = 14
+nIons = atomicNumber + 1
+
+for i in range(nIons):
+    ion_state = i + 1
+    speciesName = "%s%s" %(atomicSymbol, roman.toRoman(ion_state))
+    # Check if the species already exists
+    # in the species registry, if it does
+    # we don't want to create it again
+    if (speciesName in species_registry) == False:
+        s = Species(speciesName, atomicNumber, atomicWeight, i)
+    else:
+        s = species_registry[speciesName]
+    ion_cooling_rate(s)

File dengo/nitrogen_rates.py

+"""
+Author: Devin Silvia <devin.silvia@gmail.com>
+Affiliation: UC Boulder
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2012 Matthew Turk.  All Rights Reserved.
+
+  This file is part of the dengo package.
+
+  This file is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from reaction_classes import Species, chianti_rate, species_registry
+import docutils.utils.roman as roman
+
+# Note: the atomic/species properties have to be hard-coded for this
+# we may want to come up with a better solution here...
+atomicSymbol = 'N'
+atomicNumber = 7
+atomicWeight = 14
+nIons = atomicNumber + 1
+
+for i in range(nIons):
+    ion_state = i + 1
+    speciesName = "%s%s" %(atomicSymbol, roman.toRoman(ion_state))
+    # Check if the species already exists
+    # in the species registry, if it does
+    # we don't want to create it again
+    if (speciesName in species_registry) == False:
+        s = Species(speciesName, atomicNumber, atomicWeight, i)
+    else:
+        s = species_registry[speciesName]
+
+    if ion_state != nIons:
+        # we need to do this to make sure the 'ion_state + 1' species
+        # exists when chianti_rate is called
+        speciesNamePlusOne = "%s%s" % (atomicSymbol, roman.toRoman(ion_state+1))
+        if (speciesNamePlusOne in species_registry) == False:
+            splusone = Species(speciesNamePlusOne, atomicNumber, atomicWeight, i+1)
+
+    chianti_rate(s)

File dengo/s_ion_by_ion_cooling.h5

Binary file added.

File dengo/si_ion_by_ion_cooling.h5

Binary file added.

File dengo/silicon_cooling.py

+"""
+Author: Devin Silvia <devin.silvia@gmail.com>
+Affiliation: UC Boulder
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2012 Matthew Turk.  All Rights Reserved.
+
+  This file is part of the dengo package.
+
+  This file is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from reaction_classes import Species, ion_cooling_rate, species_registry
+import docutils.utils.roman as roman
+
+# Note: the atomic/species properties have to be hard-coded for this
+# we may want to come up with a better solution here...
+atomicSymbol = 'Si'
+atomicNumber = 14
+atomicWeight = 28
+nIons = atomicNumber + 1
+
+for i in range(nIons):
+    ion_state = i + 1
+    speciesName = "%s%s" %(atomicSymbol, roman.toRoman(ion_state))
+    # Check if the species already exists
+    # in the species registry, if it does
+    # we don't want to create it again
+    if (speciesName in species_registry) == False:
+        s = Species(speciesName, atomicNumber, atomicWeight, i)
+    else:
+        s = species_registry[speciesName]
+    ion_cooling_rate(s)

File dengo/silicon_rates.py

+"""
+Author: Devin Silvia <devin.silvia@gmail.com>
+Affiliation: UC Boulder
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2012 Matthew Turk.  All Rights Reserved.
+
+  This file is part of the dengo package.
+
+  This file is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from reaction_classes import Species, chianti_rate, species_registry
+import docutils.utils.roman as roman
+
+# Note: the atomic/species properties have to be hard-coded for this
+# we may want to come up with a better solution here...
+atomicSymbol = 'Si'
+atomicNumber = 14
+atomicWeight = 28
+nIons = atomicNumber + 1
+
+for i in range(nIons):
+    ion_state = i + 1
+    speciesName = "%s%s" %(atomicSymbol, roman.toRoman(ion_state))
+    # Check if the species already exists
+    # in the species registry, if it does
+    # we don't want to create it again
+    if (speciesName in species_registry) == False:
+        s = Species(speciesName, atomicNumber, atomicWeight, i)
+    else:
+        s = species_registry[speciesName]
+
+    if ion_state != nIons:
+        # we need to do this to make sure the 'ion_state + 1' species
+        # exists when chianti_rate is called
+        speciesNamePlusOne = "%s%s" % (atomicSymbol, roman.toRoman(ion_state+1))
+        if (speciesNamePlusOne in species_registry) == False:
+            splusone = Species(speciesNamePlusOne, atomicNumber, atomicWeight, i+1)
+
+    chianti_rate(s)

File dengo/sulfur_cooling.py

+"""
+Author: Devin Silvia <devin.silvia@gmail.com>
+Affiliation: UC Boulder
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2012 Matthew Turk.  All Rights Reserved.
+
+  This file is part of the dengo package.
+
+  This file is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from reaction_classes import Species, ion_cooling_rate, species_registry
+import docutils.utils.roman as roman
+
+# Note: the atomic/species properties have to be hard-coded for this
+# we may want to come up with a better solution here...
+atomicSymbol = 'S'
+atomicNumber = 16
+atomicWeight = 32
+nIons = atomicNumber + 1
+
+for i in range(nIons):
+    ion_state = i + 1
+    speciesName = "%s%s" %(atomicSymbol, roman.toRoman(ion_state))
+    # Check if the species already exists
+    # in the species registry, if it does
+    # we don't want to create it again
+    if (speciesName in species_registry) == False:
+        s = Species(speciesName, atomicNumber, atomicWeight, i)
+    else:
+        s = species_registry[speciesName]
+    ion_cooling_rate(s)

File dengo/sulfur_rates.py

+"""
+Author: Devin Silvia <devin.silvia@gmail.com>
+Affiliation: UC Boulder
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2012 Matthew Turk.  All Rights Reserved.
+
+  This file is part of the dengo package.
+
+  This file is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from reaction_classes import Species, chianti_rate, species_registry
+import docutils.utils.roman as roman
+
+# Note: the atomic/species properties have to be hard-coded for this
+# we may want to come up with a better solution here...
+atomicSymbol = 'S'
+atomicNumber = 16
+atomicWeight = 32
+nIons = atomicNumber + 1
+
+for i in range(nIons):
+    ion_state = i + 1
+    speciesName = "%s%s" %(atomicSymbol, roman.toRoman(ion_state))
+    # Check if the species already exists
+    # in the species registry, if it does
+    # we don't want to create it again
+    if (speciesName in species_registry) == False:
+        s = Species(speciesName, atomicNumber, atomicWeight, i)
+    else:
+        s = species_registry[speciesName]
+
+    if ion_state != nIons:
+        # we need to do this to make sure the 'ion_state + 1' species
+        # exists when chianti_rate is called
+        speciesNamePlusOne = "%s%s" % (atomicSymbol, roman.toRoman(ion_state+1))
+        if (speciesNamePlusOne in species_registry) == False:
+            splusone = Species(speciesNamePlusOne, atomicNumber, atomicWeight, i+1)
+
+    chianti_rate(s)

File examples/write_combined_network.py

+from dengo.chemical_network import \
+    ChemicalNetwork, \
+    reaction_registry, \
+    cooling_registry
+import dengo.primordial_rates, dengo.primordial_cooling
+import dengo.carbon_rates, dengo.carbon_cooling
+import dengo.nitrogen_rates, dengo.nitrogen_cooling
+import dengo.oxygen_rates, dengo.oxygen_cooling
+import dengo.neon_rates, dengo.neon_cooling
+import dengo.magnesium_rates, dengo.magnesium_cooling
+import dengo.silicon_rates, dengo.silicon_cooling
+import dengo.sulfur_rates, dengo.sulfur_cooling
+from dengo.chemistry_constants import tiny, kboltz, mh
+import numpy as np
+
+NCELLS = 1
+density = 1.0
+temperature = np.logspace(4, 6.7, NCELLS)
+temperature[:] = 5e6
+X = 1e-3
+
+combined = ChemicalNetwork()
+combined.add_energy_term()
+
+for ca in cooling_registry.values():
+    if ca.name.startswith("C") \
+    or ca.name.startswith("N") \
+    or ca.name.startswith("O") \
+    or ca.name.startswith("Ne") \
+    or ca.name.startswith("Mg") \
+    or ca.name.startswith("Si") \
+    or ca.name.startswith("S"):
+       combined.add_cooling(ca)
+
+combined.add_cooling("brem")
+combined.add_cooling("reHII")
+combined.add_cooling("reHeIII")
+combined.add_cooling("ceHI")
+combined.add_cooling("reHeII2")
+combined.add_cooling("reHeII1")
+combined.add_cooling("ciHeIS")
+combined.add_cooling("ceHeII")
+combined.add_cooling("ciHI")
+combined.add_cooling("ceHeI")
+combined.add_cooling("ciHeI")
+combined.add_cooling("ciHeII")
+
+for r in reaction_registry.values():
+    if r.name.startswith("C") \
+    or r.name.startswith("N") \
+    or r.name.startswith("O") \
+    or r.name.startswith("Ne") \
+    or r.name.startswith("Mg") \
+    or r.name.startswith("Si") \
+    or r.name.startswith("S"): 
+        combined.add_reaction(r)
+
+combined.add_reaction("k01")
+combined.add_reaction("k02")
+combined.add_reaction("k03")
+combined.add_reaction("k04")
+combined.add_reaction("k05")
+combined.add_reaction("k06")
+
+# This defines the temperature range for the rate tables
+combined.init_temperature((1e0, 1e8))
+
+# Want intermediate output?
+combined.write_intermediate_solutions = True
+
+tiny = 1e-10
+
+init_array = np.ones(NCELLS) * density
+init_values = dict()
+
+# set up initial temperatures values used to define ge
+init_values['T'] = temperature
+
+start_neutral = False
+
+if start_neutral:
+    init_values['OII']     = X * init_array
+    init_values['OIII']    = init_array * X
+    init_values['OIV']     = init_array * X
+    init_values['OV']      = init_array * X
+    init_values['OVI']     = init_array * X
+    init_values['OVII']    = init_array * X
+    init_values['OVIII']   = init_array * X
+    init_values['OIX']    = init_array * X
+    init_values['de']      = init_array * 0.0
+
+    total_density = combined.calculate_total_density(init_values, ("OI",))
+    init_values["OI"] = init_array.copy() - total_density
+    init_values = combined.convert_to_mass_density(init_values)
+else:
+    # start CIE
+    import chianti.core as ch
+    import chianti.util as chu
+
+    for s in sorted(combined.required_species):
+            if s.name != 'ge':
+                if s.name == 'de':
+                    continue
+                else:
+                    print s.name, s.number, s.free_electrons + 1
+                    ion_name = chu.zion2name(np.int(s.number),
+                                             np.int(s.free_electrons + 1))
+                    ion = ch.ion(ion_name, temperature=init_values['T'])
+                    ion.ioneqOne()
+                    ion_frac = ion.IoneqOne
+                    init_values[s.name] = ion_frac * init_array * ion.Abundance
+                
+                # in case something is negative or super small:
+                init_values[s.name][init_values[s.name] < tiny] = tiny
+
+    init_values['de']      = init_array * 0.0
+#    total_density = combined.calculate_total_density(init_values, ("OI",))
+#    init_values["OI"] = init_array.copy() - total_density
+    init_values = combined.convert_to_mass_density(init_values)
+
+init_values['de'] = combined.calculate_free_electrons(init_values)
+init_values['density'] = combined.calculate_total_density(init_values)
+number_density = combined.calculate_number_density(init_values)
+
+# calculate ge (very crudely)
+gamma = 5.0/3.0
+init_values['ge'] = ((temperature * number_density * kboltz)
+                     / (init_values['density'] * mh * (gamma - 1)))
+
+
+    # Write the initial conditions file
+combined.write_solver("combined", output_dir = ".",
+                    init_values=init_values)