diff --git a/CItest-chem.py b/CItest-chem.py
index 14875aca85529d052ecdfe837a265298cf51c7b1..340f138edee2befc2891d92d49167b6deb215bad 100755
--- a/CItest-chem.py
+++ b/CItest-chem.py
@@ -47,7 +47,6 @@ l = chemtest_reacs(
1.0, 1.0, 1.0e-3 / (1.0e-9 * defaultairdensity)
) # Lifetime about 1s for A and B, about 1e3sec for C
mech = mechanism(l)
-print(mech)
inidict = {}
inidict["A"] = initialization("A", "unifppb", (1.5,))
@@ -93,6 +92,7 @@ for exp in experimentlist:
ts_c = ((myctm.t, myctm.activeConc_ppb("C")[0, 0, 0]),)
print("running for", exp[2], "....")
+ print(myctm.mech)
while myctm.t < myctm.finaltime - 0.00001:
myctm.update()
ts_a = ts_a + ((myctm.t, myctm.activeConc_ppb("A")[0, 0, 0]),)
diff --git a/toyctm/toyctm.py b/toyctm/toyctm.py
index d8df99117cc03186d7d1bc6a64dde62210fa79c6..451c26c6277ae40deb31573c73c09dfe1c713f45 100644
--- a/toyctm/toyctm.py
+++ b/toyctm/toyctm.py
@@ -216,7 +216,12 @@ class ctm:
self.zdconc = np.zeros(self.conc.shape)
self.mdconc = np.zeros(self.conc.shape)
- # MIsc initializations
+ # Initialize the chemistry solver
+ if self.chemIntegrator is not None :
+ self.mech.define_integrator(self.chemIntegrator,
+ errortarget_relative = self.ebitolerance,
+ nitersmax = self.ebinitersmax)
+ # Misc initializations
self.niters = -999
self.vdconc = 0.
@@ -595,8 +600,7 @@ class ctm:
self.prod.fill(0.0)
self.loss.fill(0.0)
- if self.chemIntegrator == "EulerForward":
- self.mech.eulerForwardTimestep(
+ self.mech.integrator(
self.t,
self.dtchem,
self.conc,
@@ -609,24 +613,9 @@ class ctm:
self.sunrise,
self.sunset,
)
- elif self.chemIntegrator == "EulerBackward":
- self.mech.eulerBackwardTimestep(
- self.t,
- self.dtchem,
- self.ebinitersmax,
- self.ebitolerance,
- self.conc,
- self.tdgrid.activezIterator,
- self.tdgrid.hgrid.activeCellsIterator,
- ReactionConditions(self.t + self.dtchem),
- self.prod,
- self.loss,
- self.PrescSet,
- self.sunrise,
- self.sunset,
- self.niters,
- np.max(self.density) / 1.0e23,
- )
+
+ if 'niters' in self.mech.integrator_feedback :
+ self.niters = self.mech.integrator_feedback['niters']
def lagrangianupdate(self):
"""
diff --git a/toyctm_addons/chemistry.py b/toyctm_addons/chemistry.py
index 223098b2eae6c498ad18bf2861495174d216e6b4..51459ae1b24ecfc7446ce051df8dcf443184f387 100644
--- a/toyctm_addons/chemistry.py
+++ b/toyctm_addons/chemistry.py
@@ -26,11 +26,13 @@ class ReactionParameters:
def __init__(self):
pass
-
class StandardReactionParameters(ReactionParameters):
def __init__(self, k):
self.k = k
+ def __str__(self):
+ return f'reaction rate constant {self.k}'
+
class PhotoReactionParameters(ReactionParameters):
def __init__(self, peakval=None, sunrise=86400.0 / 4, sunset=3.0 * 86400.0 / 4.0):
@@ -39,6 +41,8 @@ class PhotoReactionParameters(ReactionParameters):
self.sunset = sunset
self.noon = (self.sunrise + self.sunset) / 2.0
+ def __str__(self):
+ return f'peakval: {self.peakval}, sunrise: {self.sunrise}, sunset:{self.sunset}'
def createReaction(reactant_tup, products_tup, reaction_type, parameters):
# reactant_tup : a tuple of strings
@@ -129,7 +133,7 @@ class Reaction:
return st
-class StandardReaction:
+class StandardReaction(Reaction):
def __init__(self, reactant_tup, products_tup, reaction_type, parameters):
Reaction.__init__(self, reactant_tup, products_tup, reaction_type, parameters)
@@ -137,7 +141,7 @@ class StandardReaction:
return self.parameters.k
-class PhotoReaction:
+class PhotoReaction(Reaction):
def __init__(self, reactant_tup, products_tup, reaction_type, parameters):
Reaction.__init__(self, reactant_tup, products_tup, reaction_type, parameters)
@@ -262,13 +266,53 @@ class mechanism:
self.species = {}
for sp in self.specdict:
self.species[sp] = speciesInstance(sp)
+ self.integrator = None
+ self.integrator_feedback = None
def __str__(self):
- s = ""
+ s = '-----------------------------------\n'
+ s = "List of reactions in the mechanism:\n"
for r in self.reactions:
- s = s + str(r) + "\n"
+ s = s + r.__str__() + "\n"
+
+ s = s + "List of species in the mechanism:\n"
+ for i in sorted(self.revspecdict):
+ s = s + "{i} : {spname}\n".format( i = i, spname = self.revspecdict[i] )
+
+ s = s + "Number of species: {nspec}\n".format(nspec = self.nspec)
+ if self.integrator is not None:
+ s = s + "A time integrator has been defined in this mechanism:\n"
+ s = s + self.integrator_name + '\n'
+ if self.integrator_name == 'EulerBackward' :
+ s = s + "Max number of iterations {nitersmax}\n".format(nitersmax = self.nitersmax)
+ s = s + "Target for relative error {et}\n".format(et = self.errortarget_relative)
+ s = s + "Absolute error target (concentration) {ae}\n".format(ae = self.errortarget_absolute)
+ s = s + '-----------------------------------\n'
return s
+ def define_integrator(self,integrator_name,
+ nitersmax = 1000,
+ errortarget_relative = 1.0e-2,
+ minimalconcentration = 1.e-23
+ ):
+ """Initialise an integrator. Still work in progress.
+ the purpose is to make chemistry into a standalone package
+ able to integrate chemistry."""
+ self.integrator_name = integrator_name
+ self.integrator_feedback = {}
+ if integrator_name == 'EulerBackward' :
+ self.integrator = self.eulerBackwardTimestep
+ self.integrator_feedback['niters'] = -999
+ self.nitersmax = nitersmax
+ self.errortarget_relative = errortarget_relative
+ self.minimalconcentration = minimalconcentration
+ self.errortarget_absolute = self.minimalconcentration
+ elif integrator_name == 'EulerForward' :
+ self.integrator = self.eulerForwardTimestep
+ else :
+ print('ERROR : integrator', integrator_name, 'not available')
+ sys.exit(1)
+
def count_species(self):
"""
Counts the active species in the mechanism (without double-count of species that appear more than once).
@@ -418,8 +462,6 @@ class mechanism:
self,
t,
dtchem,
- nitersmax,
- errortarget_relative,
conc,
activezIterator,
activeCellsIterator,
@@ -429,8 +471,6 @@ class mechanism:
prescriptionSet,
sunrise,
sunset,
- niters,
- minimalconcentration,
):
"""
Integrate over one time step dtchem with Euler Backward strategy
@@ -443,8 +483,6 @@ class mechanism:
# set convergence criteria
- errortarget_absolute = minimalconcentration
-
concini = np.copy(
conc
) # concini serves all along the fixed point search and should not be destroyed
@@ -473,7 +511,7 @@ class mechanism:
np.seterr(
invalid="ignore", divide="ignore"
) # A division by zero is possible where conc=0 but the resulting values are ignored due to 'where'. We do not want these warnings on screen.
- for i in range(nitersmax):
+ for i in range(self.nitersmax):
prod.fill(0.0)
loss.fill(0.0)
self.chem(
@@ -500,25 +538,28 @@ class mechanism:
error = np.max(
np.abs(
np.where(
- conc > minimalconcentration,
- (concnew - conc) / conc / errortarget_relative,
- np.abs(concnew - conc) / errortarget_absolute,
+ conc > self.minimalconcentration,
+ (concnew - conc) / conc / self.errortarget_relative,
+ np.abs(concnew - conc) / self.errortarget_absolute,
)
)
)
conc[:, :, :, :] = concnew[:, :, :, :]
if error < 1 and ntheo >= 2:
+ self.integrator_feedback['niters'] = i
break
- if i > nitersmax - 10:
+ if i > self.nitersmax - 10:
print(i, error)
np.seterr(invalid="warn", divide="ignore")
- if i == nitersmax - 1:
- exit(
- "ERROR: toyctm.py: Backward Euler integration failed to converge. Try reducing time step"
- )
+ if i == self.nitersmax - 1:
+ print("""ERROR: toyctm.py: Backward Euler integration
+ failed to converge. Try reducing time step""")
+ sys.exit(1)
+
+
prescriptionSet.apply(
activeCellsIterator,
activezIterator,