diff --git a/toyctm/toyctm.py b/toyctm/toyctm.py
index 4b6086c7e2a353db73053ac5c687577cedd01aea..d8df99117cc03186d7d1bc6a64dde62210fa79c6 100644
--- a/toyctm/toyctm.py
+++ b/toyctm/toyctm.py
@@ -1,15 +1,14 @@
#!/bin/env python
-############################################################
-# This is the toyctm module
-# It permits to create instances of the ctm class
-# These instances can be initialized and updated
-# From the calling script
-# References:
-# Mailler et al. 2021, https://doi.org/10.5194/gmd-14-2221-2021
-# Mailler and Pennel, 2021, https://hal.archives-ouvertes.fr/hal-02933095
-# https://gitlab.in2p3.fr/sylvain.mailler/toyctm
-############################################################
+""" This is the toyctm module
+ It permits to create instances of the ctm class
+ These instances can be initialized and updated
+ From the calling script
+ References:
+ Mailler et al. 2021, https://doi.org/10.5194/gmd-14-2221-2021
+ Mailler and Pennel, 2021, https://hal.archives-ouvertes.fr/hal-02933095
+ https://gitlab.in2p3.fr/sylvain.mailler/toyctm
+"""
# import glob,os
from __future__ import print_function
@@ -21,39 +20,50 @@ from .constants import avogadro
class ctm:
"""
S. Mailler 2018-2020
- ctm class contain all the concentrations and metainformation that shall be sufficient to define entirely a ctm run, including time from the beginning of the simulation, (ctm.t), concentration of each species (0d to 3d depending on grid configuration), etc. see the init routine for the exhaustive list.
-
- The two main entry points for user are ctm.init to create a new ctm instance and ctm.update to perform the evolution of the ctm over one time step.
+ ctm class contain all the concentrations and metainformation that shall be sufficient to define
+ entirely a ctm run, including time from the beginning of the simulation, (ctm.t), concentration
+ of each species (0d to 3d depending on grid configuration), etc.
+ see the init routine for the exhaustive list.
+
+ The two main entry points for user are ctm.init to create a new ctm instance and ctm.update to
+ perform the evolution of the ctm over one time step.
The main elements of ctm are :
self.tdgrid : three-dimensional grid (instance of the Threedimgrid class).
- self.tdgrid.hgrid : 2d drid underlying tdgrid
- self.conc : concentration array. First index is species index, last index is vertical dimension, intermediate indices are horizontal and vertical space dimensions
+ self.tdgrid.hgrid : 2d drid underlying tdgrid
+ self.conc : concentration array. First index is species index,
+ last index is vertical dimension,
+ intermediate indices are horizontal and vertical space dimensions
self.mech : chemical mechanism, of class mechanism defined in chemistry.py
- self.ts : the transport scheme instance of class TransportScheme defined in ctmtransport.py
- self.t : the current time
+ self.ts : the transport scheme instance of class TransportScheme defined in
+ ctmtransport.py
+ self.t : the current time
self.dt : the time step
"""
def __init__(self, *args, **kwargs):
"""
- This routine initializes a ctm object from user specification.
- Mandatory arguments (in arbitrary order) are as follow:
+ This routine initializes a ctm object from user specification.
+ Mandatory arguments (in arbitrary order) are as follow:
- mech : object of mechanism class, decribing the reactions
- tdgrid : grid object describing the 3d ctm grid
- - inidict : initial values for all species (dictionary of instances in the initialization class)
+ - inidict : initial values for all species (dictionary of instances in the
+ initialization class)
- bgdict : background values for all species (dictionary of real values)
- - emidict : emission dictionary for all species (dictionary of instances in the emission class)
+ - emidict : emission dictionary for all species (dictionary of instances in the
+ emission class)
- ts : object of the TransportScheme class. default : TransportScheme('upwind')
- - vts : object of the TransportScheme class for vertical transport. default : TransportScheme('upwind')
+ - vts : object of the TransportScheme class for vertical transport. default :
+ TransportScheme('upwind')
- vertperiod : boolean, defaults to False.
- - prescdict : Dictionary for prescribed concentrations
+ - prescdict : Dictionary for prescribed concentrations
- dt : time step in seconds
- - sunrise, sunset : seconds after midnight for sunrise and sunset (for photochemical reactions)
+ - sunrise, sunset : seconds after midnight for sunrise and sunset (for photochemical reactions)
- MassFluxFunc : Instance of the MassFluxGenerator class (examples in threedmassfluxes.py)
- lagdict : defaults to {}
- label : string tha gives a name to the simulation
- - chemIntegrator : time scheme for chemistry. So far available, EulerBackward (strongly recommended) and EulerForward
+ - chemIntegrator : time scheme for chemistry. So far available,
+ EulerBackward (strongly recommended) and EulerForward
"""
# read Mandatory arguments
@@ -65,7 +75,7 @@ class ctm:
) # Three-dimensional grid : instance of the Threedimdgrid class
self.MassFluxFunc = kwargs.get(
"MassFluxFunc", None
- ) # Provider of mass flux and density : instance of the MassFluxGenerator class.
+ ) # Provider of mass flux and density : instance of the MassFluxGenerator class
self.ts = kwargs.get(
"ts", None
) # Horizontal trasnport scheme : instance of the TransportScheme class
@@ -83,54 +93,62 @@ class ctm:
print(
"toyctm.py ERROR: Please provide a chemical mechanism with mech= when calling ctm()"
)
- exit()
+ sys.exit(1)
if self.tdgrid is None:
print(
"toyctm.py ERROR: Please provide a grid with tdgrid= when calling ctm()"
)
- exit()
+ sys.exit(1)
if self.MassFluxFunc is None:
print(
- "toyctm.py ERROR: Please provide a Mass flux generator with MassFluxFunc= when calling ctm()"
+ """toyctm.py ERROR: Please provide a Mass flux generator with
+ MassFluxFunc= when calling ctm()"""
)
- exit()
+ sys.exit(1)
if self.ts is None:
print(
- "toyctm.py ERROR: Please provide a horizontal transport scheme with ts= when calling ctm()"
+ """toyctm.py ERROR: Please provide a horizontal transport scheme with
+ ts= when calling ctm()"""
)
- exit()
+ sys.exit(1)
if self.vts is None:
print(
- "toyctm.py ERROR: Please provide a vertical transport scheme with vts= when calling ctm()"
+ """toyctm.py ERROR: Please provide a vertical transport scheme with
+ vts= when calling ctm()"""
)
- exit()
+ sys.exit(1)
if self.PrescSet is None:
print(
- "toyctm.py ERROR: Please provide a Prescription set with prescset= when calling ctm()"
+ """toyctm.py ERROR: Please provide a Prescription set with
+ prescset= when calling ctm()"""
)
- exit()
+ sys.exit(1)
if self.dt is None:
print(
"toyctm.py ERROR: Please provide a time step with dt= when calling ctm()"
)
- exit()
+ sys.exit(1)
elif self.dt <= 0.0:
print("toyctm.py ERROR: please provide dt>0. with dt= when calling ctm()")
- exit()
+ sys.exit(1)
# read optional arguments
self.chemIntegrator = kwargs.get(
"chemIntegrator", "EulerBackward"
- ) # Chemical integrator. 'EulerBackward', 'EulerForward' or None are supported. 'EulerBackward' is strongly recommended if you have chemistry.
+ ) # Chemical integrator. 'EulerBackward', 'EulerForward' or None are supported.
+ # 'EulerBackward' is strongly recommended if you have chemistry.
inidict = kwargs.get(
"inidict", {}
- ) # Dictionnary associating an instance of the initialization class to each model species. Optional, 0. initial value is applied otherwise.
+ ) # Dictionnary associating an instance of the initialization class to each model species.
+ # Optional, 0. initial value is applied otherwise.
bgdict = kwargs.get(
"bgdict", {}
- ) # Dictionnary associating a background value to each model species. Optional, 0. otherwise.
+ ) # Dictionnary associating a background value to each model species.
+ # Optional, 0. otherwise.
self.emidict = kwargs.get(
"emidict", {}
- ) # Dictionnary associating an instance of the emission class to each model species. Optional, 0. otherwise.
+ ) # Dictionnary associating an instance of the emission class to each model species.
+ # Optional, 0. otherwise.
self.sunrise = kwargs.get(
"sunrise", 6.0 * 86400 / 24.0
) # Sunrise time (for photochemical reactions)
@@ -160,7 +178,7 @@ class ctm:
self.update = self.lieUpdate
else:
print("toyctm.py ERROR: only split=strang or split=lie supported so far")
- exit()
+ sys.exit(1)
# Check ndtchem and calculate dtchem
if self.ndtchem is None:
@@ -169,11 +187,11 @@ class ctm:
else:
self.ndtchem = 1
if self.split == "strang":
- if not (self.ndtchem % 2 == 0):
+ if not self.ndtchem % 2 == 0 :
print(
"toyctm.py ERROR: with Strang splitting ndtchem must be an even number"
)
- exit()
+ sys.exit(1)
self.dtchem = self.dt / self.ndtchem
# Some convenient shortcuts
@@ -182,8 +200,11 @@ class ctm:
self.species = list(self.mech.specdict.keys())
self.specdict = self.mech.specdict # Use by some old scripts
- self.t = 0.0 # Time in seconds since the beginning of the experiment, which also has to be time since midnight (for traffic emissions and photolytic emissions)
- self.maxcfl = 0.0 # Maximal value of the cfl number during the simulation. initialized to 0.
+ self.t = 0.0 # Time in seconds since the beginning of the experiment,
+ # which also has to be time since midnight
+ # (for traffic emissions and photolytic emissions)
+ self.maxcfl = 0.0 # Maximal value of the cfl number during the simulation.
+ # initialized to 0.
# Initialization of the main arrays
self.density = np.zeros(self.tdgrid.mass_shape())
@@ -197,6 +218,7 @@ class ctm:
# MIsc initializations
self.niters = -999
+ self.vdconc = 0.
# Real data between 0 and nhalo-1
@@ -232,7 +254,7 @@ class ctm:
s = s + "species =" + str(self.species) + "\n"
s = s + "---------------------------------------------------------------\n"
s = s + "EMISSIONS\n" # Emissions overview
- if not (self.emidict == {}):
+ if not self.emidict == {} :
for spec in self.emidict:
s = s + str(self.emidict[spec]) + "\n"
else:
@@ -267,7 +289,8 @@ class ctm:
def init_bg(self, bgdict):
"""
- Fills bgvals array with values from bgdict. bgvals is used by updateBC of the twodgrid module to set the boundary conditions
+ Fills bgvals array with values from bgdict. bgvals is used by updateBC of the
+ twodgrid module to set the boundary conditions
"""
for spec in bgdict:
@@ -276,15 +299,16 @@ class ctm:
def init_PrescSet(self, prescdict):
"""
- Creates ctm.PrescSet that gives the prescribed values for species that have a prescribed concentration value instead of chemistry
- """
+ Creates ctm.PrescSet that gives the prescribed values for species that
+ have a prescribed concentration value instead of chemistry
+ """
if prescdict is not None:
self.PrescSet = PrescriptionSet(prescdict)
def init_bound(self):
"""
- Fills the boundaries of the domain with background values from array ctm.bgvals
- """
+ Fills the boundaries of the domain with background values from array ctm.bgvals
+ """
for iz in self.tdgrid.activezIterator:
for isp in range(self.nspec):
for indices in self.grid.backgroundCellsIterator:
@@ -343,7 +367,7 @@ class ctm:
intermdensity = intermdensity + self.dt * mddensity
# 3. One complete step vertical transport
- if not (self.vts.name == "null"):
+ if not self.vts.name == "null" :
self.tdgrid.updateBC(
self.conc, self.bgvals, self.nspec
) # Indispensable in case of periodic BC
@@ -378,10 +402,10 @@ class ctm:
0.5 meridional step
0.5 zonal step
- If vertical transport is absent, simplifies to:
+ If vertical transport is absent, simplifies to:
0.5 zonal step
1. meridional step
- 0.5 zonal step
+ 0.5 zonal step
"""
intermdensity = np.zeros((1,) + self.density.shape)
@@ -556,15 +580,18 @@ class ctm:
): # emidict[sp] should be an instance of the emission class
self.emidict[sp].putemismap(
self
- ) # Takes the ctm as argument. A closed list of arguments should be provided for modularity.
+ ) # Takes the ctm as argument. A closed list of arguments should be provided for
+ # modularity.
self.prod = self.emis / self.tdgrid.celldz
# update concentration after one time step
- # To make things cleaner, emissions should be the average between t and t+dt and not an instantaneous value.
+ # To make things cleaner, emissions should be the average between t and t+dt
+ # and not an instantaneous value.
self.conc = self.conc + self.prod * self.dtchem
def chemstep(self):
+ """Perform chemistry for one timestep"""
self.prod.fill(0.0)
self.loss.fill(0.0)
@@ -603,7 +630,8 @@ class ctm:
def lagrangianupdate(self):
"""
- Routine to follow Lagrangian particles if self.MassFluxFunc.lagrangianAdvect is provided in threedmassfluxes.py for the current mass flux generator
+ Routine to follow Lagrangian particles if self.MassFluxFunc.lagrangianAdvect
+ is provided in threedmassfluxes.py for the current mass flux generator
"""
for key in self.lagdict:
coords = self.lagdict[key]
@@ -620,13 +648,13 @@ class ctm:
def lieUpdate(self):
"""
- Routine update of the ctm class performs one time step. So far this consists in
- - Integrate chemistry over dt
- - Integrate emissions over dt
- - Integrate transport over dt
- - Update Lagrangian tracers bw t and t+dt
- - t <- t + dt.
- """
+ Routine update of the ctm class performs one time step. So far this consists in
+ - Integrate chemistry over dt
+ - Integrate emissions over dt
+ - Integrate transport over dt
+ - Update Lagrangian tracers bw t and t+dt
+ - t <- t + dt.
+ """
# Advect Lagrangian particles if relevant
self.lagrangianupdate()
@@ -669,7 +697,8 @@ class ctm:
self.sunset,
)
- def chemupdate(self): # Everything except transport on time dtchem
+ def chemupdate(self):
+ """Do Everything except transport during time dtchem"""
# EMISSIONS half time step
if self.emidict: # test if disctionary is not empty
self.emi()
@@ -685,10 +714,9 @@ class ctm:
): # Typical observed machine truncature is ~ e-18
print(sp, mini, maxi, "NEGATIVE VALUES AFTER CHEMISTRY")
- def chemupdate_reverse(
- self,
- ): # Everything except transport on time dtchem in reverse order
- # dtchem for chemistry
+ def chemupdate_reverse(self):
+ """ Everything except transport on time dtchem in reverse order
+ during time dtchem"""
if self.chemIntegrator is not None:
# CHEMISTRY
self.chemstep()
@@ -705,11 +733,11 @@ class ctm:
def strangUpdate(self):
"""
Routine update of the ctm class performs one time step. So far this consists in
- - Integrate chemistry over dt
+ - Integrate chemistry over dt
- Integrate emissions over dt
- Integrate transport over dt
- Update Lagrangian tracers bw t and t+dt
- - t <- t + dt.
+ - t <- t + dt.
"""
# Advect Lagrangian over dt if relevant
@@ -757,28 +785,30 @@ class ctm:
)
def minmax(self, spec, unit=None):
- """
- returns the min and max of species concentration over the
- entire active domain in molecules per cubic meter
- This method is used within the model for some checks
- """
+ """
+ returns the min and max of species concentration over the
+ entire active domain in molecules per cubic meter
+ This method is used within the model for some checks
+ """
if unit is None:
- return self.tdgrid.minmax(self.conc[self.mech.specdict[spec], :, :, :])
+ retval = self.tdgrid.minmax(self.conc[self.mech.specdict[spec], :, :, :])
if unit == "ppb":
- return self.tdgrid.minmax(
- self.conc[self.mech.specdict[spec], :, :, :] / self.density * 1.0e9
- )
+ retval = self.tdgrid.minmax(
+ self.conc[self.mech.specdict[spec], :, :, :] / self.density * 1.0e9
+ )
+ return retval
- # Non-critical methods ("gadgets") to answer users request (give easier access to quantities inside the ctm instance)
+ # Non-critical methods ("gadgets") to answer users request (give easier access to
+ # quantities inside the ctm instance)
def moleculesNumber(self, spec, mask=1.0):
- """ Returns Number of molecules of species spec, over all the active grid area,
+ """ Returns Number of molecules of species spec, over all the active grid area,
masked by the optional mask argument"""
return self.tdgrid.integrVolumicScalar(
mask * self.conc[self.mech.specdict[spec], :, :, :]
)
def mass(self, spec, mask=1.0):
- """
+ """
Returns mass of species spec, in kg, over the entire active domain
"""
return (
@@ -788,11 +818,13 @@ class ctm:
)
def extractConc(self, spec):
- """ Returns numpy array of shape self.tdgrid.mass_shape() with the concentration field for species spec """
+ """ Returns numpy array of shape self.tdgrid.mass_shape() with the concentration field for
+ species spec """
return np.copy(self.conc[self.mech.specdict[spec], :, :, :])
def activeConc_ppb(self, spec):
- """ Returns numpy array of field shape self.tdgrid.active_field_shape with all the useful concentrations of species spec expressed in ppb"""
+ """ Returns numpy array of field shape self.tdgrid.active_field_shape with all the useful
+ concentrations of species spec expressed in ppb"""
return np.copy(
self.tdgrid.extract_active(
self.conc[self.mech.specdict[spec], :, :, :] / self.density * 1.0e9
@@ -800,11 +832,12 @@ class ctm:
)
def activeDensity(self):
- """ Returns numpy array of field shape self.tdgrid.active_field_shape with all the useful density values in molecules per unit volume"""
+ """ Returns numpy array of field shape self.tdgrid.active_field_shape with all the useful
+ density values in molecules per unit volume"""
return np.copy(self.tdgrid.extract_active(self.density))
def centerConc(self, spec):
- """ Returns the mixing ratio of species spec in ppb at the lowest model
+ """ Returns the mixing ratio of species spec in ppb at the lowest model
level and in horizontally closest to the domain center
"""
ic, jc = self.tdgrid.hgrid.centerIndices()
@@ -813,22 +846,23 @@ class ctm:
]
def centerDensity(self):
- """ Returns the mixing ratio of species spec in ppb at the lowest model
+ """ Returns the mixing ratio of species spec in ppb at the lowest model
level and in horizontally closest to the domain center
"""
ic, jc = self.tdgrid.hgrid.centerIndices()
return self.density[ic, jc, self.tdgrid.activezIterator[0]]
def centerConc_ppb(self, spec):
- """ Returns the mixing ratio of species spec in ppb at the lowest model
+ """ Returns the mixing ratio of species spec in ppb at the lowest model
level and in horizontally closest to the domain center
"""
return self.centerConc(spec) / self.centerDensity() * 1.0e9
def returnslice(self, spec, unit="number_conc"):
- """plotslice(self,spec,fname) gives a vertical slice of species spec and the corner coordinates for later use with pcolormesh
- This routine is not portable, no design effort has been made, Quick and dirty
- Actually not easy to make it clean because it mixes horizontal ad vertical properties."""
+ """plotslice(self,spec,fname) gives a vertical slice of species spec and the corner
+ coordinates for later use with pcolormesh
+ This routine is not portable, no design effort has been made, Quick and dirty
+ Actually not easy to make it clean because it mixes horizontal ad vertical properties."""
ns = self.mech.specdict[spec]
if self.grid.dim == 2:
iy = int(np.floor(self.grid.ny / 2))
@@ -862,16 +896,16 @@ class ctm:
return xcorners2d, zcorners, scalarfield
- else:
- print("Vertical slice plots not available for unstructured grids yet")
- exit("Vertical slice plots not available for unstructured grids yet")
+ print("ERROR: Vertical slice plots not available for unstructured grids yet")
+ sys.exit(1)
class ReactionConditions:
+ """this class contains the reaction conditions that may impact reaction rates"""
def __init__(self, t=None, n2=None, o2=None, airm=None, h2o=None, temperature=None):
self.t = t
self.n2 = n2
self.o2 = o2
self.airm = airm
self.h2o = h2o
- self.temperature = None
+ self.temperature = temperature