import copy
import os
import re
from operator import xor
from os.path import exists, getsize
import numpy as np
import pandas as pd
[docs]
def create_mandchem(chemistry, mandatory_files):
"""
Create the mandatory chemistry files using the chemistry config yml file
"""
for mfile in mandatory_files:
os.system("touch {}".format(mfile))
# Chemical reactions
with open(mandatory_files[0], "w") as f:
if chemistry.nreacs > 0:
for attr in chemistry.reactions.attributes:
reac = getattr(chemistry.reactions, attr)
f.write(reac + "\n")
if hasattr(chemistry, "prescrconcs"):
with open(mandatory_files[1], "w") as f:
for attr in chemistry.prescrconcs.attributes:
f.write(attr + "\n")
if hasattr(chemistry, "prodloss3d"):
with open(mandatory_files[2], "w") as f:
for attr in chemistry.prodloss3d.attributes:
f.write(attr + "\n")
if hasattr(chemistry, "deposition"):
with open(mandatory_files[3], "w") as f:
for attr in chemistry.deposition.attributes:
f.write(attr + "\n")
[docs]
def create_optchem(chemistry, filer, fileps):
"""
Create the chemistry files ALL_SPECIES, CHEMISTRY, REACTION_RATES,
STOICHIOMETRY and ACTIVE_SPECIES
At the moment, FAMILIES is set empty
ALL_SPECIES is created
Args:
chemistry (pycif.utils.classes.chemistries.Chemistry): chemical scheme
filer (str) : path to the file with reactions
fileps (str): path to the file with prescribed species
"""
workdir = chemistry.workdir
dirchem_ref = chemistry.dirchem_ref
mecachim = chemistry.schemeid
files = "{}/STOICHIOMETRY.{}".format(dirchem_ref, mecachim)
filec = "{}/CHEMISTRY.{}".format(dirchem_ref, mecachim)
filerr = "{}/REACTION_RATES.{}".format(dirchem_ref, mecachim)
filej = "{}/PHOTO_RATES.{}".format(dirchem_ref, mecachim)
filephot = "{}/PHOTO_PARAMETERS.{}".format(dirchem_ref, mecachim)
filedpar = "{}/DEPO_PARS.{}".format(dirchem_ref, mecachim)
filewdep = "{}/WETD_SPEC.{}".format(dirchem_ref, mecachim)
filef = "{}/FAMILIES.{}".format(dirchem_ref, mecachim)
fileals = "{}/ALL_SPECIES.{}".format(dirchem_ref, mecachim)
fileas = "{}/ACTIVE_SPECIES.{}".format(dirchem_ref, mecachim)
fileout = "{}/OUTPUT_SPECIES.{}".format(dirchem_ref, mecachim)
fileanth = "{}/ANTHROPIC.{}".format(dirchem_ref, mecachim)
filebio = "{}/BIOGENIC.{}".format(dirchem_ref, mecachim)
os.system("rm -f {} {} {} {}".format(files, filec, filerr, filef))
# Read prescribed species
prescribed_species = np.array([])
if exists(fileps) and getsize(fileps) > 0:
df_prescr = pd.read_csv(
fileps,
index_col=False,
header=None,
comment="#",
engine="python",
sep="\s+",
)
prescribed_species = df_prescr[0]
nonactive = np.append(prescribed_species, ["O2", "X", None])
# Read reactions
nphoto_rates = 0
with open(filer, "r") as f:
nlines = len(f.readlines())
if nlines > 0:
df_reac = pd.read_csv(
filer,
index_col=False,
header=None,
comment="#",
engine="python",
sep="\s+",
)
# Create stoechiometry and chemistry
parts = df_reac[0].str.split("->", n=2, expand=True)
left_hand_sides = parts[0]
right_hand_sides = parts[1]
losses = left_hand_sides.str.split("+", expand=True)
nlosses = left_hand_sides.str.split("+").apply(len)
losses.insert(0, -1, nlosses)
prods = right_hand_sides.str.split("+", expand=True)
nprods = pd.Series(np.sum(~pd.isnull(prods).values, axis=1))
prods.insert(0, -1, nprods)
prods_spec = np.copy(prods.values)
stoichiometry = np.zeros((1, 6))
for (i, j), value in np.ndenumerate(prods.values[:, 1:]):
try:
prods_list = prods.values[i, j + 1].split("*")
if len(prods_list) > 1:
arr = np.array(
[
prods_list[1],
prods_list[0],
prods_list[0],
prods_list[0],
prods_list[0],
i + 1,
]
)
stoichiometry = np.append(
stoichiometry, arr[np.newaxis, ...], axis=0
)
if prods_list[-1] not in nonactive:
prods_spec[i, j + 1] = prods_list[-1]
else:
prods_spec[i, j + 1] = None
prods_spec[i, 0] -= 1
except AttributeError:
pass
file_chemistry = np.append(losses.values, prods_spec, axis=1)
file_chemistry = pd.DataFrame(file_chemistry)
stoichiometry = pd.DataFrame(stoichiometry[1:])
# Create reactions
reactions = df_reac[1].values
reactions_rates, nphoto_rates, photo_rates = read_react(
reactions, filej)
reactions_rates = pd.DataFrame(reactions_rates)
reactions_rates = reactions_rates.astype({0: 'int32'})
chemistry.stoichiometry = stoichiometry
chemistry.reactions_rates = reactions_rates
chemistry.photo_rates = photo_rates
# Create files to create
stoichiometry.to_csv(files, sep=" ", header=False, index=False)
file_chemistry.to_csv(filec, sep=" ", header=False, index=False)
reactions_rates.to_csv(filerr, sep=" ", header=False, index=False)
# If no reaction, make empty files
else:
open(files, "w").close()
open(filec, "w").close()
open(filerr, "w").close()
open(filej, "w").close()
# Make Families
if chemistry.families.attributes != []:
raise Exception("Families not implemented yet in pyCIF")
open(filef, "w").close()
# Make empty DEPO_PARS
open(filedpar, "w").close()
open(filewdep, "w").close()
# Photolysis rates
pparam = chemistry.photo_parameters
pparam.nlevphot = len(pparam.altiphot)
with open(filephot, "w") as f:
f.write(
" ".join(
map(
str,
[pparam.ntabuzen, pparam.nphot, pparam.nwavel, pparam.nlevphot]
+ pparam.altiphot,
)
)
)
# Active species
chemistry.nspecies = len(chemistry.acspecies.attributes)
# Create active_species (output_species) and all_species
active_species = np.array(chemistry.acspecies.attributes)
output_species = np.array(chemistry.outspecies.attributes)
emis_species = np.array(chemistry.emis_species.attributes)
bio_species = np.array(chemistry.bio_species.attributes)
# Check that output/emis/bio species are active species
if np.any(np.isin(emis_species, active_species, invert=True)):
raise Exception(
"Some species in emis_species are not active species. " "Check your yaml!"
)
if np.any(np.isin(bio_species, active_species, invert=True)):
raise Exception(
"Some species in bio_species are not active species. " "Check your yaml!"
)
if np.any(np.isin(output_species, active_species, invert=True)):
raise Exception(
"Some species in output_species are not active species. " "Check your yaml!"
)
# Fill active species with proper information
all_species = np.append(active_species, prescribed_species, axis=0)
type_spec = np.array([["type"]])
acinfos = np.array([0, 0, 0])
outinfos = np.array([0])
type_spec = np.broadcast_to(type_spec, (all_species.shape[0], 1))
acinfos = np.broadcast_to(acinfos, (active_species.shape[0], 3))
outinfos = np.broadcast_to(outinfos, (output_species.shape[0], 1))
all_species = np.append(all_species[:, np.newaxis], type_spec, axis=1)
active_species = np.append(active_species[:, np.newaxis], acinfos, axis=1)
output_species = np.append(output_species[:, np.newaxis], outinfos, axis=1)
# Check whether LMDZ type or CHIMERE
is_lmdz = False
for i in range(active_species.shape[0]):
acspec = getattr(chemistry.acspecies, active_species[i, 0])
if xor(hasattr(active_species, "restart_id"), hasattr(active_species, "mass")):
raise Exception(
"WARNING: the active species {} has one of the arguments: "
"restart_id / mass, but not the other, suggesting you want to "
"set the chemistry for LMDZ, but not fully consistently. "
"Please check your Yaml"
)
# For LMDZ
if hasattr(acspec, "restart_id"):
is_lmdz = True
active_species[i, 1] = str(getattr(acspec, "restart_id"))
active_species[i, 2] = str(getattr(acspec, "mass"))
# For CHIMERE
else:
if is_lmdz:
raise Exception(
"Another species was defined for LMDZ and this one is "
"for CHIMERE, please check your Yaml"
)
active_species[i, 1] = str(getattr(acspec, "htransport"))
active_species[i, 2] = str(getattr(acspec, "vtransport"))
active_species[i, 3] = str(int(getattr(acspec, "bound_dry")))
if active_species[i, 0] in output_species[:, 0]:
outspec = getattr(chemistry.outspecies, active_species[i, 0])
output_species[i, 1] = str(getattr(outspec, "output_frac"))
for i in range(all_species.shape[0]):
if i < chemistry.nspecies:
all_species[i, 1] = "ac"
else:
all_species[i, 1] = "pr"
all_species = pd.DataFrame(all_species)
output_species = pd.DataFrame(output_species)
active_species = pd.DataFrame(active_species)
emis_species = pd.DataFrame(emis_species)
bio_species = pd.DataFrame(bio_species)
# Create files to create
all_species.to_csv(fileals, sep=" ", header=False, index=False)
active_species.to_csv(
fileas,
sep=" ",
header=False,
)
output_species.to_csv(fileout, sep=" ", header=False, index=False)
emis_species.to_csv(fileanth, sep=" ", header=False)
bio_species.to_csv(filebio, sep=" ", header=False)
return all_species.shape[0], nphoto_rates
[docs]
def read_react(reactions, filej):
"""
Parse the types of reactions in the REACTIONS file and extract the data
"""
nreacts = reactions.shape[0]
reactions_rates = np.zeros((nreacts, 10))
reactions_rates.fill(None)
photo_rates = np.zeros((0, 2))
nphoto_rates = 0
for i, value in np.ndenumerate(reactions):
# Type 1 : Constant rate
if re.search(r"k=", value):
val_list = re.split(r"=", value)
reactions_rates[i, 0] = 1
reactions_rates[i, 1] = val_list[1]
# Type 2 : Simplified Arrhenius
if re.search(r"k\(T\)=Aexp\(-B\/T\)", value):
val_list = re.split(r"=|,", value)
reactions_rates[i, 0] = 2
reactions_rates[i, 1] = val_list[3]
reactions_rates[i, 2] = val_list[5]
# Type 3 : Complete Arrhenius
if re.search(r"k\(T\)=Aexp\(-B\/T\)\(300\/T\)\*\*N", value):
val_list = re.split(r"=|,", value)
reactions_rates[i, 0] = 3
reactions_rates[i, 1] = val_list[3]
reactions_rates[i, 2] = val_list[5]
reactions_rates[i, 3] = val_list[7]
# Type 4 : Relative pressure
if re.search(r"k\(P\)=A\(B\+C\*P/Pref\)", value):
val_list = re.split(r"=|,", value)
reactions_rates[i, 0] = 4
reactions_rates[i, 1] = val_list[3]
reactions_rates[i, 2] = val_list[5]
reactions_rates[i, 3] = val_list[7]
# Type 5 : Simple Photolysis
if re.search(r"J=", value):
val_list = re.split(r"=|,", value)
reactions_rates[i, 0] = 5
reactions_rates[i, 1] = val_list[1]
photo_rates = np.append(
photo_rates, [[str(i[0] + 1), "j" + val_list[1]]], axis=0
)
nphoto_rates += 1
photo_rates = pd.DataFrame(photo_rates[:])
photo_rates.to_csv(filej, sep=" ", header=False, index=False)
return reactions_rates, nphoto_rates, photo_rates
[docs]
def make_inout_react_graph(
filer
):
# Check if REACTIONS files is empty
with open(filer, "r") as f:
nlines = len(f.readlines())
if nlines == 0:
return {}
# Otherwhise, start building the graph
inout_reaction_graph = {}
with open(filer, "r") as f:
lines = f.readlines()
for ln in lines:
ln = ln.split()
nloss = int(ln[0])
losses = [s for s in ln[1:nloss + 1] if s not in ["X", "M"]]
nprods = int(ln[nloss + 1])
prods = [s for s in ln[nloss + 2:] if s not in ["X", "M"]]
# Remove smaller 'o' in front of species if any
losses = [s[1:] if s[0] == 'o' else s for s in losses]
prods = [s[1:] if s[0] == 'o' else s for s in prods]
# Now link prods and losses
for trin in losses:
if trin not in inout_reaction_graph:
inout_reaction_graph[trin] = []
inout_reaction_graph[trin].extend(prods)
inout_reaction_graph = {
tr: list(set(inout_reaction_graph[tr]))
for tr in inout_reaction_graph
}
# Propagate recursively influence from inputs to outputs
def propagate_chemistry(inout_reaction_graph):
inout_reaction_graph_out = copy.deepcopy(inout_reaction_graph)
for trin in inout_reaction_graph_out:
for trout in inout_reaction_graph_out[trin]:
inout_reaction_graph_out[trin].extend(
inout_reaction_graph_out.get(trout, []))
inout_reaction_graph_out[trin] = list(
set(inout_reaction_graph_out[trin])
)
# Compare length and carries on recursively if any changes
lengths = [
len(inout_reaction_graph_out[tr]) == len(inout_reaction_graph[tr])
for tr in inout_reaction_graph_out
]
if all(lengths):
return inout_reaction_graph_out
else:
return propagate_chemistry(inout_reaction_graph_out)
return propagate_chemistry(inout_reaction_graph)
