Skip to content

add NL99_GC (untested) #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
BenWibking wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

add NL99_GC (untested) #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
290 changes: 290 additions & 0 deletions NL99_GC_rates.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,290 @@
# ABOUTME: Convert NL99_GC chemistry network from DESPOTIC to GPUAstroChem format
# ABOUTME: Creates SymPy-based reaction rates for C++ code generation

import sympy as sp
import numpy as np
from rate_new import ChemSpecie, SympyChemRate

# Define symbolic variables for GPUAstroChem
T = sp.symbols('Tgas', positive=True, real=True)
Te = sp.symbols('Te', positive=True, real=True)
invTe = sp.symbols('invTe', positive=True, real=True)
invT = sp.symbols('invT', positive=True, real=True)
lnTe = sp.symbols('lnTe', real=True)
T32 = sp.symbols('T32', positive=True, real=True)
invT32 = sp.symbols('invT32', positive=True, real=True)
invsqrT = sp.symbols('invsqrT', positive=True, real=True)
user_crate = sp.symbols('user_crate', positive=True, real=True)
user_Av = sp.symbols('user_Av', positive=True, real=True)
user_ionH = sp.symbols('user_ionH', positive=True, real=True)
user_ionH2 = sp.symbols('user_ionH2', positive=True, real=True)
user_dissH2 = sp.symbols('user_dissH2', positive=True, real=True)
user_ionC = sp.symbols('user_ionC', positive=True, real=True)
user_ionO = sp.symbols('user_ionO', positive=True, real=True)
user_dissCO = sp.symbols('user_dissCO', positive=True, real=True)
user_dust2gas_ratio = sp.symbols('user_dust2gas_ratio', positive=True, real=True)
Tdust = sp.symbols('Tdust', positive=True, real=True)

# Define all species used in NL99_GC network
# Note: GPUAstroChem uses different naming conventions
species_map = {
'H+': 'hp',
'H': 'h',
'H-': 'hm',
'H2': 'h2',
'H2+': 'h2p',
'H3+': 'h3p',
'He+': 'hep',
'He++': 'hepp',
'He': 'he',
'C+': 'cp',
'C': 'c',
'CHx': 'ch', # Simplified - NL99 treats CHx as single species
'OHx': 'oh', # Simplified - NL99 treats OHx as single species
'CO': 'co',
'HCO+': 'hcop',
'O': 'o',
'M+': 'mp', # Metal ion
'M': 'm', # Neutral metal
'e-': 'elec'
}

# Create rate list to store all reactions
rates = []

# Cosmic ray reactions from NL99_GC.py
# (1) cr + H -> H+ + e-
rates.append(SympyChemRate(
reactants=[ChemSpecie('h')],
products=[ChemSpecie('hp'), ChemSpecie('elec')],
rate_expr=user_crate * 1.0
))

# (2) cr + He -> He+ + e-
rates.append(SympyChemRate(
reactants=[ChemSpecie('he')],
products=[ChemSpecie('hep'), ChemSpecie('elec')],
rate_expr=user_crate * 1.1
))

# (3) cr + H2 -> H3+ + H (simplified from H2+ + H)
rates.append(SympyChemRate(
reactants=[ChemSpecie('h2')],
products=[ChemSpecie('h3p'), ChemSpecie('h')],
rate_expr=user_crate * 2.0
))

# Photoreactions from NL99_GC.py with shielding
# (1) h nu + H2 -> H + H
# Note: Shielding functions would need to be implemented separately
rates.append(SympyChemRate(
reactants=[ChemSpecie('h2')],
products=[ChemSpecie('h'), ChemSpecie('h')],
rate_expr=user_dissH2 * 3.3e-11 * 1.7 * sp.exp(-3.74 * user_Av)
))

# (2) h nu + CO -> C + O
rates.append(SympyChemRate(
reactants=[ChemSpecie('co')],
products=[ChemSpecie('c'), ChemSpecie('o')],
rate_expr=user_dissCO * 1.2e-10 * 1.7 * sp.exp(-3.53 * user_Av)
))

# (3) h nu + C -> C+ + e-
rates.append(SympyChemRate(
reactants=[ChemSpecie('c')],
products=[ChemSpecie('cp'), ChemSpecie('elec')],
rate_expr=user_ionC * 1.8e-10 * 1.7 * sp.exp(-3.0 * user_Av)
))

# (4) h nu + CHx -> C + H
rates.append(SympyChemRate(
reactants=[ChemSpecie('ch')],
products=[ChemSpecie('c'), ChemSpecie('h')],
rate_expr=1.0e-9 * sp.exp(-1.5 * user_Av)
))

# (5) h nu + OHx -> O + H
rates.append(SympyChemRate(
reactants=[ChemSpecie('oh')],
products=[ChemSpecie('o'), ChemSpecie('h')],
rate_expr=5.0e-10 * sp.exp(-1.7 * user_Av)
))

# (7) h nu + HCO+ -> CO + H+
rates.append(SympyChemRate(
reactants=[ChemSpecie('hcop')],
products=[ChemSpecie('co'), ChemSpecie('hp')],
rate_expr=1.5e-10 * sp.exp(-2.5 * user_Av)
))

# Two-body reactions from NL99_GC.py
# (1) H3+ + C -> CHx + H2
rates.append(SympyChemRate(
reactants=[ChemSpecie('h3p'), ChemSpecie('c')],
products=[ChemSpecie('ch'), ChemSpecie('h2')],
rate_expr=2.0e-9
))

# (2) H3+ + O -> OHx + H2
rates.append(SympyChemRate(
reactants=[ChemSpecie('h3p'), ChemSpecie('o')],
products=[ChemSpecie('oh'), ChemSpecie('h2')],
rate_expr=8.0e-10
))

# (3) H3+ + CO -> HCO+ + H2
rates.append(SympyChemRate(
reactants=[ChemSpecie('h3p'), ChemSpecie('co')],
products=[ChemSpecie('hcop'), ChemSpecie('h2')],
rate_expr=1.7e-9
))

# (4) He+ + H2 -> He + H + H+
rates.append(SympyChemRate(
reactants=[ChemSpecie('hep'), ChemSpecie('h2')],
products=[ChemSpecie('he'), ChemSpecie('h'), ChemSpecie('hp')],
rate_expr=7.0e-15
))

# (5) He+ + CO -> C+ + O + He
rates.append(SympyChemRate(
reactants=[ChemSpecie('hep'), ChemSpecie('co')],
products=[ChemSpecie('cp'), ChemSpecie('he'), ChemSpecie('o')],
rate_expr=1.4e-9 * (T/300.0)**(-0.5)
))

# (6) C+ + H2 -> CHx + H
rates.append(SympyChemRate(
reactants=[ChemSpecie('cp'), ChemSpecie('h2')],
products=[ChemSpecie('ch'), ChemSpecie('h')],
rate_expr=4.0e-16
))

# (7) C+ + OHx -> HCO+
rates.append(SympyChemRate(
reactants=[ChemSpecie('cp'), ChemSpecie('oh')],
products=[ChemSpecie('hcop')],
rate_expr=1.0e-9
))

# (8) O + CHx -> CO + H
rates.append(SympyChemRate(
reactants=[ChemSpecie('o'), ChemSpecie('ch')],
products=[ChemSpecie('co'), ChemSpecie('h')],
rate_expr=2.0e-10
))

# (9) C + OHx -> CO + H
rates.append(SympyChemRate(
reactants=[ChemSpecie('c'), ChemSpecie('oh')],
products=[ChemSpecie('co'), ChemSpecie('h')],
rate_expr=5.0e-12 * T**0.5
))

# (10) He+ + e- -> He
# Temperature-dependent recombination rate
logT = sp.log(T, 10)
rates.append(SympyChemRate(
reactants=[ChemSpecie('hep'), ChemSpecie('elec')],
products=[ChemSpecie('he')],
rate_expr=1.0e-11 / sp.sqrt(T) * (11.19 - 1.676*logT - 0.2852*logT**2 + 0.04433*logT**3)
))

# (11) H3+ + e- -> H2 + H
rates.append(SympyChemRate(
reactants=[ChemSpecie('h3p'), ChemSpecie('elec')],
products=[ChemSpecie('h2'), ChemSpecie('h')],
rate_expr=2.34e-8 * (T/300.0)**(-0.52)
))

# (12) H3+ + e- -> 3H
rates.append(SympyChemRate(
reactants=[ChemSpecie('h3p'), ChemSpecie('elec')],
products=[ChemSpecie('h'), ChemSpecie('h'), ChemSpecie('h')],
rate_expr=4.36e-8 * (T/300.0)**(-0.52)
))

# (13) C+ + e- -> C
rates.append(SympyChemRate(
reactants=[ChemSpecie('cp'), ChemSpecie('elec')],
products=[ChemSpecie('c')],
rate_expr=4.67e-12 * (T/300.0)**(-0.6)
))

# (14) HCO+ + e- -> CO + H
rates.append(SympyChemRate(
reactants=[ChemSpecie('hcop'), ChemSpecie('elec')],
products=[ChemSpecie('co'), ChemSpecie('h')],
rate_expr=2.76e-7 * (T/300.0)**(-0.64)
))

# (15) H+ + e- -> H
rates.append(SympyChemRate(
reactants=[ChemSpecie('hp'), ChemSpecie('elec')],
products=[ChemSpecie('h')],
rate_expr=2.753e-14 * (315614.0/T)**1.5 * (1.0+(115188.0/T)**0.407)**(-2.242)
))

# (16) H2 + H -> 3H (collisional dissociation)
# This is a complex rate with density dependence - simplified version
rates.append(SympyChemRate(
reactants=[ChemSpecie('h2'), ChemSpecie('h')],
products=[ChemSpecie('h'), ChemSpecie('h'), ChemSpecie('h')],
rate_expr=6.67e-12 * sp.sqrt(T) * sp.exp(-(1.0+63590.0/T))
))

# (17) H2 + H2 -> H2 + 2H
rates.append(SympyChemRate(
reactants=[ChemSpecie('h2'), ChemSpecie('h2')],
products=[ChemSpecie('h2'), ChemSpecie('h'), ChemSpecie('h')],
rate_expr=5.996e-30 * T**4.1881 / (1.0+6.761e-6*T)**5.6881 * sp.exp(-54657.4/T)
))

# (18) H + e- -> H+ + 2e- (collisional ionization)
rates.append(SympyChemRate(
reactants=[ChemSpecie('h'), ChemSpecie('elec')],
products=[ChemSpecie('hp'), ChemSpecie('elec'), ChemSpecie('elec')],
rate_expr=sp.exp(-32.71396786 + 13.5365560*lnTe - 5.73932875*lnTe**2 +
1.56315498*lnTe**3 - 0.2877056*lnTe**4 + 0.0348255977*lnTe**5 -
2.63197617e-3*lnTe**6 + 1.11954395e-4*lnTe**7 - 2.03914985e-6*lnTe**8)
))

# (19) He+ + H2 -> H+ + He + H
rates.append(SympyChemRate(
reactants=[ChemSpecie('hep'), ChemSpecie('h2')],
products=[ChemSpecie('hp'), ChemSpecie('he'), ChemSpecie('h')],
rate_expr=3.7e-14 * sp.exp(-35.0/T)
))

# (20) H + H -> H2 (grain catalyzed)
# Grain surface formation rate
fA = 1.0/(1.0+1.0e4*sp.exp(-600/Tdust))
rates.append(SympyChemRate(
reactants=[ChemSpecie('h'), ChemSpecie('h')],
products=[ChemSpecie('h2')],
rate_expr=3.0e-18 * T**0.5 * fA * user_dust2gas_ratio /
(1.0+0.04*(T+Tdust)**0.5 + 0.002*T+8.0e-6*T**2)
))

# (21) H+ + e- -> H (grain catalyzed)
# Simplified grain recombination rate based on NL99_GC._twobody_ratecoef[20]
# This avoids complex dependencies and provides a cleaner implementation
psi = user_Av * sp.sqrt(T) / (1e-6 + 1e-40) # Simplified psi calculation
rates.append(SympyChemRate(
reactants=[ChemSpecie('hp'), ChemSpecie('elec')],
products=[ChemSpecie('h')],
rate_expr=12.25e-14 * user_dust2gas_ratio /
(1.0 + 8.074e-6 * psi**1.378 *
(1.0 + 508.7 * T**0.01586 * psi**(-0.4723-1.102e-5*sp.log(T))))
))

# Export rates for use in GPUAstroChem
def get_NL99_GC_rates():
return rates

if __name__ == "__main__":
print(f"Created {len(rates)} reactions for NL99_GC network")
print("Species included:")
for desp_name, gpu_name in species_map.items():
print(f" {desp_name} -> {gpu_name}")
Loading