diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 4bfecfedbc3cd35b637371f7b267459f76908e86..f08ce4098e1464dc53396d39d3127e6c406136c1 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -29,6 +29,10 @@ before_script:
- pip install numpy
- pip install matplotlib
+stages:
+ - test
+ - diags
+
test-advection:
stage: test
script:
@@ -141,5 +145,107 @@ test-PrescEmi:
- artifacts/PrescEmi.logs
- artifacts/PrescEmi_*.png
-
+test-2dtest-ppmw:
+ stage: test
+ script:
+ - python 2dtestcases.py --ts ppmw --windfield leveque --nx 5 --ny 5 --period 86400. --initshape smooth --chemistry ozoneco --split strang --dsize 1.e+5 --makedumps --label leveque-ozoneco-smooth-ppmw > artifacts/ppmw.log
+ artifacts:
+ paths:
+ - artifacts/CO_leveque-ozoneco-smooth-ppmw-086400.png
+ - artifacts/NO2_leveque-ozoneco-smooth-ppmw-086400.png
+ - artifacts/NO_leveque-ozoneco-smooth-ppmw-086400.png
+ - artifacts/O3_leveque-ozoneco-smooth-ppmw-086400.png
+ - artifacts/TRC_leveque-ozoneco-smooth-ppmw-086400.png
+ - artifacts/leveque-ozoneco-smooth-ppmw-*.pickle
+ - artifacts/ppmw.log
+
+test-2dtest-ppmSL:
+ stage: test
+ script:
+ - python 2dtestcases.py --ts ppmSL --windfield leveque --nx 5 --ny 5 --period 86400. --initshape smooth --chemistry ozoneco --split strang --dsize 1.e+5 --makedumps --label leveque-ozoneco-smooth-ppmSL > artifacts/ppmSL.log
+
+ artifacts:
+ paths:
+ - artifacts/CO_leveque-ozoneco-smooth-ppmSL-086400.png
+ - artifacts/NO2_leveque-ozoneco-smooth-ppmSL-086400.png
+ - artifacts/NO_leveque-ozoneco-smooth-ppmSL-086400.png
+ - artifacts/O3_leveque-ozoneco-smooth-ppmSL-086400.png
+ - artifacts/TRC_leveque-ozoneco-smooth-ppmSL-086400.png
+ - artifacts/leveque-ozoneco-smooth-ppmSL-*.pickle
+ - artifacts/ppmSL.log
+
+test-2dtest-upwind:
+ stage: test
+ script:
+ - python 2dtestcases.py --ts upwind --windfield leveque --nx 5 --ny 5 --period 86400. --initshape smooth --chemistry ozoneco --split strang --dsize 1.e+5 --makedumps --label leveque-ozoneco-smooth-upwind > artifacts/upwind.log
+ artifacts:
+ paths:
+ - artifacts/CO_leveque-ozoneco-smooth-upwind-086400.png
+ - artifacts/NO2_leveque-ozoneco-smooth-upwind-086400.png
+ - artifacts/NO_leveque-ozoneco-smooth-upwind-086400.png
+ - artifacts/O3_leveque-ozoneco-smooth-upwind-086400.png
+ - artifacts/TRC_leveque-ozoneco-smooth-upwind-086400.png
+ - artifacts/leveque-ozoneco-smooth-upwind-*.pickle
+ - artifacts/upwind.log
+
+test-2dtest-vanleerSL:
+ stage: test
+ script:
+ - python 2dtestcases.py --ts vanleerSL --windfield leveque --nx 5 --ny 5 --period 86400. --initshape smooth --chemistry ozoneco --split strang --dsize 1.e+5 --makedumps --label leveque-ozoneco-smooth-vanleerSL > artifacts/vanleerSL.log
+
+ artifacts:
+ paths:
+ - artifacts/CO_leveque-ozoneco-smooth-vanleerSL-086400.png
+ - artifacts/NO2_leveque-ozoneco-smooth-vanleerSL-086400.png
+ - artifacts/NO_leveque-ozoneco-smooth-vanleerSL-086400.png
+ - artifacts/O3_leveque-ozoneco-smooth-vanleerSL-086400.png
+ - artifacts/TRC_leveque-ozoneco-smooth-vanleerSL-086400.png
+ - artifacts/leveque-ozoneco-smooth-vanleerSL-*.pickle
+ - artifacts/vanleerSL.log
+
+test-2dtest-walcek:
+ stage: test
+ script:
+ - python 2dtestcases.py --ts walcek --windfield leveque --nx 5 --ny 5 --period 86400. --initshape smooth --chemistry ozoneco --split strang --dsize 1.e+5 --makedumps --label leveque-ozoneco-smooth-walcek > artifacts/walcek.log
+
+ artifacts:
+ paths:
+ - artifacts/CO_leveque-ozoneco-smooth-walcek-086400.png
+ - artifacts/NO2_leveque-ozoneco-smooth-walcek-086400.png
+ - artifacts/NO_leveque-ozoneco-smooth-walcek-086400.png
+ - artifacts/O3_leveque-ozoneco-smooth-walcek-086400.png
+ - artifacts/TRC_leveque-ozoneco-smooth-walcek-086400.png
+ - artifacts/leveque-ozoneco-smooth-walcek-*.pickle
+ - artifacts/walcek.log
+
+test-2dtest-null:
+ stage: test
+ script:
+ - python 2dtestcases.py --ts null --windfield leveque --nx 5 --ny 5 --period 86400. --initshape smooth --chemistry ozoneco --split strang --dsize 1.e+5 --makedumps --label leveque-ozoneco-smooth-null > artifacts/null.log
+ artifacts:
+ paths:
+ - artifacts/CO_leveque-ozoneco-smooth-null-086400.png
+ - artifacts/NO2_leveque-ozoneco-smooth-null-086400.png
+ - artifacts/NO_leveque-ozoneco-smooth-null-086400.png
+ - artifacts/O3_leveque-ozoneco-smooth-null-086400.png
+ - artifacts/TRC_leveque-ozoneco-smooth-null-086400.png
+ - artifacts/leveque-ozoneco-smooth-null-*.pickle
+ - artifacts/null.log
+
+diags-2d:
+ stage: diags
+ needs: ["test-2dtest-null", "test-2dtest-walcek", "test-2dtest-vanleerSL", "test-2dtest-upwind", "test-2dtest-ppmSL", "test-2dtest-ppmw"]
+ script: ./2d-diags.py --CIcheck > artifacts/2d-diags.log
+ artifacts:
+ paths:
+ - artifacts/norm2.tex
+ - artifacts/norm1.tex
+ - artifacts/e1.tex
+ - artifacts/2d-diags.log
+ - artifacts/pdfdist*.png
+ - artifacts/finalpdf_fancy_*.png
+ - artifacts/cminmax.png
+ - artifacts/testnminmax.png
+ - artifacts/map-*-086400.png
+
\ No newline at end of file
diff --git a/2d-diags.py b/2d-diags.py
new file mode 100755
index 0000000000000000000000000000000000000000..ed538ce931674c822ecc1b9e114e2c02536f31b1
--- /dev/null
+++ b/2d-diags.py
@@ -0,0 +1,1250 @@
+#!/bin/env python
+
+# Standard modules
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+import sys
+import os
+import pickle
+import argparse
+
+# Toyctm
+from toyctm.toyctm import ctm # Main module. Defines the ctm class
+
+# Toyctm addons
+from toyctm_addons.initialization import initialization
+from toyctm_addons.chemistry import mechanism, reaction
+from toyctm_addons.twodgrid import Latlongrid
+from toyctm_addons.threedimgrid import Threedimgrid
+from toyctm_addons.ctmtransport import TransportScheme
+from toyctm_addons.threedmassfluxes import MassFluxGenerator, MassFluxParams
+from toyctm_addons.defaults import defaults
+import toyctm.constants as constants
+
+plt.ioff()
+
+# Set the default text font size
+plt.rc("font", size=14)
+
+# Set the axes title font size
+plt.rc("axes", titlesize=14)
+
+# Set the axes labels font size
+plt.rc("axes", labelsize=14)
+
+# Set the font size for x tick labels
+plt.rc("xtick", labelsize=10)
+
+# Set the font size for y tick labels
+plt.rc("ytick", labelsize=10)
+
+# Set the legend font size
+plt.rc("legend", fontsize=8)
+
+# Set the font size of the figure title
+plt.rc("figure", titlesize=16)
+
+#############################################################
+# Initialize some stuff to avoid pickle load errors
+def makeplotchem(
+ self, sp, simlabel="toyctm", unit="ppb", spname="FAMILY"
+): # unit should be ppb or None
+ print("making plot for", sp)
+ if type(sp) is tuple:
+ sptup = sp
+ else:
+ sptup = (sp,)
+ spname = sp
+
+ if unit == "ppb":
+ factor = 1.0
+ else:
+ factor = 1.0e-9
+
+ tlabel = "{0}-{1:06d}".format(simlabel, int(self.t + 0.5))
+ field = None
+ for spec in sptup:
+ print("Extracting field", spec)
+ if field is None:
+ field = self.activeConc_ppb(spec) * factor
+ else:
+ field = field + self.activeConc_ppb(spec) * factor
+
+ ax = plt.figure(figsize=(4.5, 4.0))
+ plt.pcolormesh(
+ self.tdgrid.hgrid.activeloncorners / 1000.0,
+ self.tdgrid.hgrid.activelatcorners / 1000.0,
+ field[:, :, 0],
+ cmap=plt.cm.Reds,
+ )
+ plt.xlabel("x (km)")
+ plt.ylabel("y (km)")
+ cbar = plt.colorbar()
+ if unit == "ppb":
+ cbar.set_label(r"$\alpha{}_{" + spname + r"} \mathrm{(ppb)}$")
+ else:
+ cbar.set_label(r"$\alpha{}_{" + spname + r"}$")
+ plt.savefig(
+ "{path}map-{spname}-{tlabel}.png".format(
+ path=diagspath, spname=spname, tlabel=tlabel
+ ),
+ dpi=300,
+ bbox_inches="tight",
+ )
+ plt.savefig(
+ "{path}map-{spname}-{tlabel}.pdf".format(
+ path=diagspath, spname=spname, tlabel=tlabel
+ ),
+ bbox_inches="tight",
+ )
+ plt.close()
+
+
+ctm.makeplotchem = makeplotchem
+
+#############################################################
+# Initialize some stuff to avoid pickle load errors
+
+
+def psi_translation(inst, x, y, t, dt):
+ return inst.params.V0 * x - inst.params.U0 * y
+
+
+def psi_solid(inst, x, y, t, dt):
+ rsq = (x - inst.params.xo) ** 2.0 + (y - inst.params.yo) ** 2.0
+ return np.pi * rsq / inst.params.period
+
+
+def psi_leveque(inst, x, y, t, dt):
+ U0 = inst.params.U0
+ L0 = inst.params.L0
+ xo = inst.params.xo
+ yo = inst.params.yo
+ T = inst.params.period
+ prof = (
+ U0
+ * L0
+ * np.cos(0.5 * (x - inst.params.xo) * np.pi / L0) ** 2.0
+ * np.cos(0.5 * (y - inst.params.yo) * np.pi / L0) ** 2.0
+ )
+ if not (dt == 0.0):
+ ta = T / (np.pi * dt) * (np.sin(np.pi * (t + dt) / T) - np.sin(np.pi * t / T))
+ else:
+ ta = np.cos(np.pi * t / T)
+ return prof * ta
+
+
+def frompsi(inst, t, dt, x1, y1, zbot1, ztop1, x2, y2, zbot2, ztop2):
+ zbot = 0.5 * (zbot1 + zbot2) # Per surface du trapeze this should work
+ ztop = 0.5 * (ztop1 + ztop2)
+ return (
+ (ztop - zbot)
+ * inst.params.density
+ * (inst.psi(inst, x1, y1, t, dt) - inst.psi(inst, x2, y2, t, dt))
+ )
+
+
+def plotpsi(self, xmin, xmax, ymin, ymax, t, fname=None):
+ nxplot = 100
+ nyplot = (
+ 110 # nxplot slighty differs from nyplot to avoid any permutation in dimensions
+ )
+ xvec = np.linspace(xmin, xmax, nxplot + 1)
+ yvec = np.linspace(ymin, ymax, nyplot + 1)
+ psiar = np.zeros((nxplot + 1, nyplot + 1))
+ xar = np.zeros((nxplot + 1, nyplot + 1))
+ yar = np.zeros((nxplot + 1, nyplot + 1))
+ if fname is None:
+ fname = "psiplots/psi_" + self.params.name + ".png"
+
+ for i in range(nxplot + 1):
+ for j in range(nyplot + 1):
+ psiar[i, j] = self.psi(self, xvec[i], yvec[j], t, 0.0)
+ xar[i, j] = xvec[i]
+ yar[i, j] = yvec[j]
+
+ gradx, grady = np.gradient(psiar, xvec, yvec)
+ uspeed = -grady
+ vspeed = gradx
+ module = np.sqrt(uspeed * uspeed + vspeed * vspeed)
+ maxmod = np.max(module)
+ minmod = np.min(module)
+ meanspeed = (maxmod + minmod) / 2.0
+ ampli = np.abs((maxmod - minmod) / maxmod)
+ plt.figure(figsize=(4.5, 4.0))
+ print(xvec.shape)
+ print(yvec.shape)
+ print(psiar.shape)
+ plt.contour(xar / 1000.0, yar / 1000.0, psiar, colors="black")
+ if ampli > 0.001:
+ plt.contourf(xar / 1000.0, yar / 1000.0, module, cmap=plt.cm.viridis)
+ else:
+ plt.contourf(
+ xar / 1000.0,
+ yar / 1000.0,
+ module,
+ [meanspeed * 0.999, meanspeed * 1.001],
+ cmap=plt.cm.viridis,
+ ) # We do not want to plot numerical noise
+
+ cbar = plt.colorbar()
+ cbar.set_label("U (m/s)")
+ plt.xlabel("x (km)")
+ plt.ylabel("y (km)")
+ plt.quiver(
+ xar[::10, ::10] / 1000.0,
+ yar[::10, ::10] / 1000.0,
+ uspeed[::10, ::10],
+ vspeed[::10, ::10],
+ )
+ plt.savefig(fname + ".png", dpi=300, bbox_inches="tight")
+ plt.savefig(fname + ".pdf", bbox_inches="tight")
+ plt.close()
+
+
+# Enhance MassFluxGenerator class
+MassFluxGenerator.plotpsi = plotpsi
+
+TransportScheme.__computeFluxSL = None
+TransportScheme.__computeFluxMOL = None
+
+
+# END OF PICKLE WORKAROUNDS
+#
+
+#############################################################
+# Scatter class
+#
+
+
+class Scatter:
+ def __init__(self, x, y, t=0.0, unitx="ppb", unity="ppb", spx=None, spy=None):
+ # x, y provided in ppb
+
+ # Sort by increasing x to permit connected plots
+ xvec = np.asarray(x)
+ yvec = np.asarray(y)
+ p = xvec.argsort()
+
+ self.x = xvec[p]
+ self.y = yvec[p]
+ self.t = t
+ self.spx = spx
+ self.spy = spy
+ self.unitx = unitx
+ self.unity = unity
+ if unitx == None:
+ self.x = self.x / 1.0e9
+ if unity == None:
+ self.y = self.y / 1.0e9
+
+ def plot(self, simlabel="toyctm", color="k", refscat=None, lw=1.0):
+ fig, ax1 = plt.subplots(1, 1, figsize=(4.5, 4.0))
+ tlabel = "{0}-{1:06d}".format(simlabel, int(t + 0.5))
+ ax1.scatter(self.x, self.y, color=color, s=1)
+ ax1.set_title(nicelabel[simdict[lab].tsname])
+ if refscat is not None:
+ ax1.plot(refscat.x, refscat.y, color="k", lw=lw)
+ if self.unitx is None:
+ ax1.set_xlim(0.0, 1.0)
+ if self.unity is None:
+ ax1.set_ylim(0.0, 1.0)
+ plt.savefig(
+ "{0}{1}-{2}-{3}-{4}.png".format(diagspath, "scat", sp1, sp2, tlabel),
+ dpi=300,
+ bbox_inches="tight",
+ )
+ plt.savefig(
+ "{0}{1}-{2}-{3}-{4}.pdf".format(diagspath, "scat", sp1, sp2, tlabel),
+ bbox_inches="tight",
+ )
+ plt.close()
+
+
+#############################################################
+# Simulation class
+#
+class simulation:
+ def __init__(self, myctm):
+ self.label = myctm.label
+ self.tsname = myctm.ts.name
+ self.times = []
+ self.o3max = []
+ self.o3min = []
+ self.cmin = []
+ self.cmax = []
+ self.n2min = []
+ self.n2max = []
+ self.testnmin = []
+ self.testnmax = []
+ self.omass = []
+ self.hmass = []
+ self.nmass = []
+ self.n2mass = []
+ self.mmass = []
+ self.o3mass = []
+ self.nomass = []
+ self.no2mass = []
+ self.trcmass = []
+ self.cmass = []
+ self.trcmax = []
+ self.trchalf = []
+ self.niters = []
+ self.pdfdist = None
+ self.scatters = {}
+
+ def add_diags(self, myctm, myctmref=None):
+ self.times.append(
+ np.fix(myctm.t + 0.5)
+ ) # round up to the nearest integer second to avoid truncature errors
+
+ # Mass conservation
+ self.omass.append(
+ 2.0 * myctm.moleculesNumber("NO2")
+ + myctm.moleculesNumber("NO")
+ + myctm.moleculesNumber("O")
+ + 2.0 * myctm.moleculesNumber("O2")
+ + 3.0 * myctm.moleculesNumber("O3")
+ + myctm.moleculesNumber("O1D")
+ + myctm.moleculesNumber("H2O")
+ + 2.0 * myctm.moleculesNumber("CO2")
+ + myctm.moleculesNumber("CO")
+ + myctm.moleculesNumber("OH")
+ + 2.0 * myctm.moleculesNumber("HO2")
+ + 2.0 * myctm.moleculesNumber("H2O2")
+ + 3.0 * myctm.moleculesNumber("HNO3")
+ )
+ self.hmass.append(
+ 2.0 * myctm.moleculesNumber("H2O")
+ + myctm.moleculesNumber("OH")
+ + myctm.moleculesNumber("HO2")
+ + 2.0 * myctm.moleculesNumber("H2O2")
+ + myctm.moleculesNumber("HNO3")
+ )
+ self.nmass.append(
+ myctm.moleculesNumber("NO")
+ + myctm.moleculesNumber("NO2")
+ + myctm.moleculesNumber("HNO3")
+ )
+ self.n2mass.append(2.0 * myctm.moleculesNumber("N2"))
+ self.mmass.append(myctm.moleculesNumber("M"))
+ self.o3mass.append(myctm.moleculesNumber("O3"))
+ self.nomass.append(myctm.moleculesNumber("NO"))
+ self.no2mass.append(myctm.moleculesNumber("NO2"))
+ self.cmass.append(myctm.moleculesNumber("CO") + myctm.moleculesNumber("CO2"))
+ self.trcmass.append(myctm.moleculesNumber("TRC"))
+
+ # cmin and cmax conservation
+ carray = myctm.activeConc_ppb("CO") + myctm.activeConc_ppb("CO2")
+ cmin, cmax = np.amin(carray), np.amax(carray)
+ self.cmin.append(cmin)
+ self.cmax.append(cmax)
+
+ # testnmin and testnmax conservation
+ testnarray = (
+ myctm.activeConc_ppb("NO")
+ + myctm.activeConc_ppb("NO2")
+ + myctm.activeConc_ppb("HNO3")
+ + myctm.activeConc_ppb("TRCb")
+ )
+ testnmin, testnmax = np.amin(testnarray), np.amax(testnarray)
+ self.testnmin.append(testnmin)
+ self.testnmax.append(testnmax)
+
+ # n2min and n2max conservation
+ n2array = myctm.activeConc_ppb("N2")
+ n2min, n2max = np.amin(n2array), np.amax(n2array)
+ self.n2min.append(n2min)
+ self.n2max.append(n2max)
+
+ o3array = myctm.activeConc_ppb("O3")
+ o3minconc, o3maxconc = np.amin(o3array), np.amax(o3array)
+ self.o3max.append(o3maxconc)
+ self.o3min.append(o3minconc)
+
+ self.trcmax.append(np.amax(myctm.activeConc_ppb("TRC")))
+ self.trchalf.append(halfvolume(myctm, "TRC"))
+ if myctm.t > 0.0:
+ self.niters.append(myctm.niters)
+ else:
+ self.niters.append(np.nan)
+
+ if myctmref is not None:
+ if self.pdfdist is None:
+ self.pdfdist = {}
+ for sp in accurspecs:
+ self.pdfdist[sp] = []
+ for sp in accurspecs:
+ self.pdfdist[sp].append(pdfdist(sp, myctm, myctmref))
+
+ def addscatter(self, myctm, sp1, sp2):
+ if sp1 == "NOx":
+ self.scatters[sp1, sp2, myctm.t] = Scatter(
+ myctm.activeConc_ppb("NO").flatten()
+ + myctm.activeConc_ppb("NO2").flatten(),
+ myctm.activeConc_ppb(sp2).flatten(),
+ myctm.t,
+ unitx="ppb",
+ unity="ppb",
+ spx=sp1,
+ spy=sp2,
+ )
+ elif sp1 == "TRC1" and sp2 == "TRC2":
+ self.scatters[sp1, sp2, myctm.t] = Scatter(
+ myctm.activeConc_ppb(sp1).flatten(),
+ myctm.activeConc_ppb(sp2).flatten(),
+ myctm.t,
+ unitx=None,
+ unity=None,
+ spx=sp1,
+ spy=sp2,
+ )
+ else:
+ self.scatters[sp1, sp2, myctm.t] = Scatter(
+ myctm.activeConc_ppb(sp1),
+ myctm.activeConc_ppb(sp2).flatten(),
+ myctm.t,
+ unitx="ppb",
+ unity="ppb",
+ spx=sp1,
+ spy=sp2,
+ )
+
+
+# END SIMULATION CLASS DEFINITION
+
+
+def dict2latex(normdict, reflabel, norm="norm1"):
+ out = ""
+ spset = set()
+ labset = set()
+
+ for label in normdict:
+ lab, sp = label
+ if lab == reflabel:
+ continue
+ spset.add(sp)
+ labset.add(lab)
+ sptup = simdict
+
+ splist = sorted(list(spset), reverse=True) # we like to have TRC and O3 first
+ lablist = sorted(list(labset))
+
+ ccc = "c"
+ for sp in splist:
+ ccc = ccc + "c"
+
+ line = "\\begin{{tabular}}{{ {cols} }}".format(cols=ccc)
+ out = out + line + "\n"
+
+ line = "\\hline"
+ out = out + line + "\n"
+
+ line = " "
+ for sp in splist:
+ line = line + "&"
+ line = line + " {specname} ".format(specname=sp)
+ line = line + "\\\\"
+ out = out + line + "\n"
+
+ line = "\\hline"
+ out = out + line + "\n"
+
+ for lab in lablist:
+ line = "{label} ".format(label=lab)
+ for sp in splist:
+ line = line + "&"
+ line = line + " {:1.3} ".format(normdict[lab, sp])
+ line = line + "\\\\"
+ out = out + line + "\n"
+
+ line = "\\hline"
+ out = out + line + "\n"
+
+ line = "\\end{tabular}"
+ out = out + line + "\n"
+
+ return out
+
+
+def pdfdist(spec, myctm, myctmref):
+ conc = myctm.activeConc_ppb(spec)
+ concref = myctmref.activeConc_ppb(spec)
+ cflat = np.sort(conc, axis=None) # sort the flattened array in increasing order
+ crefflat = np.sort(concref, axis=None)
+ return np.sum(np.abs(cflat - crefflat)) / np.sum(
+ np.abs(crefflat)
+ ) # The result is adimensionnal typically between 0 and 1
+
+
+def plotpdfs(myctm, myctmref, prefix="pdfdisttribs-", fancy=False):
+
+ for sp in accurspecs:
+ conc = np.sort(
+ myctm.activeConc_ppb(sp), axis=None
+ ) # sort the flattened array in increasing order
+ concref = np.sort(myctmref.activeConc_ppb(sp), axis=None)
+ norm1err = np.sum(np.abs(conc - concref)) / np.sum(np.abs(concref))
+ pdf = ( np.cumsum(np.ones(conc.shape)) -1. ) / conc.size
+
+ # Add initial and final point
+ a = np.array([0.0])
+ z = np.array([1.0])
+ conc = np.concatenate((a, conc, np.array([conc[-1]])))
+ concref = np.concatenate((a, concref, np.array([concref[-1]])))
+ pdf = np.concatenate((a, pdf, z))
+
+ fig, ax1 = plt.subplots(1, 1, figsize=(4.5, 4.0))
+ ax1.step(
+ conc,
+ pdf,
+ label=nicelabel[myctm.ts.name],
+ color=colordict[myctm.ts.name],
+ linewidth=1.0,
+ )
+ if sp in {"TRC", "TRC1", "TRC2"}:
+ ax1.step(
+ concref,
+ pdf,
+ label=r"$\mathcal{S}^0$",
+ color=colordict[myctmref.ts.name],
+ linewidth=1.0,
+ )
+ else:
+ ax1.step(
+ concref,
+ pdf,
+ label=nicelabel[myctmref.ts.name],
+ color=colordict[myctmref.ts.name],
+ linewidth=1.0,
+ )
+
+ ax1.set_ylim(0.0, 1.0)
+ ax1.set_xlim(-0.002 * np.max(concref), np.max(concref))
+ ax1.text(
+ 0.5,
+ 0.01,
+ r"$\mathcal{{E}}_1={:1.3}$".format(norm1err),
+ verticalalignment="bottom",
+ horizontalalignment="center",
+ transform=ax1.transAxes,
+ color="green",
+ fontsize=15,
+ )
+ if fancy:
+ for i in range(conc.size - 1):
+ # An additional point for origin has beed added to conc so the
+ # actual number of cells is conc.size-1
+ ymin, ymax = pdf[i], pdf[i + 1]
+ xmin, xmax = min(conc[i], concref[i]), max(conc[i], concref[i])
+ ax1.add_patch(
+ mpl.patches.Rectangle(
+ (xmin, ymin),
+ width=xmax - xmin,
+ height=ymax - ymin,
+ edgecolor=None,
+ facecolor="orange",
+ linewidth=0.75,
+ )
+ )
+ plt.rcParams["hatch.linewidth"] = 0.30
+ ax1.fill_between(
+ concref,
+ pdf,
+ y2=1.0,
+ step="pre",
+ hatch="/////",
+ zorder=2,
+ edgecolor="black",
+ facecolor="none",
+ linewidth=0.0,
+ )
+ plt.legend()
+ plt.xlabel(r"$\alpha{}\ \mathrm{(ppb)}$")
+ plt.ylabel(r"$\mathcal{S}^t$")
+ # plt.show()
+ plt.savefig(
+ diagspath + prefix + myctm.ts.name + "-" + sp + ".png",
+ dpi=300,
+ bbox_inches="tight",
+ )
+ plt.savefig(
+ diagspath + prefix + myctm.ts.name + "-" + sp + ".pdf", bbox_inches="tight"
+ )
+ plt.close()
+
+
+def halfvolume(myctm, spec):
+ conc = myctm.activeConc_ppb(spec) * myctm.activeDensity() / 1.0e9
+ # Here we suppose that all the cell volumes are equal and we produce a volume normalized by the volume of entire box
+ cflat = np.flip(
+ np.sort(conc, axis=None)
+ ) # sort the flattened array in decreasing order
+ totconc = np.sum(cflat)
+ ncells = cflat.shape[0]
+ partconc = 0.0
+ for i in range(ncells):
+ partconc = partconc + cflat[i]
+ if partconc > 0.5 * totconc:
+ break
+ partconc = partconc - cflat[i]
+ iaccur = float(i) + (0.5 * totconc - partconc) / cflat[i]
+ return iaccur / float(ncells)
+
+
+# parameters
+
+#############################################################
+# ARGUMENTS MANAGEMENT
+parser = argparse.ArgumentParser(
+ description="""This script performs a 2d test case for horizontal advection in a square domain, for three possible advection schemes and three possible wind fields. Graphic outputs are in psiplots.""",
+ epilog="""Refer to the README for more details (maybe)""",
+)
+
+parser.add_argument(
+ "--picklespath",
+ type=str,
+ default="artifacts/",
+ help="Relative or absolute path where the pickles are to be found",
+)
+
+parser.add_argument(
+ "--diagspath",
+ type=str,
+ default="artifacts/",
+ help="Relative or absolute path to store the outputs",
+)
+
+parser.add_argument('--CIcheck',
+ action = 'store_true',
+ help = 'Perform CI check for expected convergence rates')
+
+args = parser.parse_args()
+
+
+picklespath = args.picklespath
+diagspath = args.diagspath
+
+colordict = {}
+colordict["null"] = "g"
+colordict["ppmSL"] = "b"
+colordict["vanleerSL"] = "r"
+colordict["walcek"] = "m"
+colordict["upwind"] = "c"
+colordict["ppmw"] = "y"
+
+nicelabel = {}
+nicelabel["null"] = "Base"
+nicelabel["ppmSL"] = "PPM"
+nicelabel["vanleerSL"] = "Van Leer"
+nicelabel["walcek"] = "Walcek"
+nicelabel["upwind"] = "Godunov"
+nicelabel["ppmw"] = "PPM+W"
+
+accurspecs = {"O3", "NO", "NO2", "TRC", "TRC1", "TRC2"}
+
+plt.rc("axes", titlesize=16) # Increase axis title font for readability
+
+# Read pickles, list labels and times
+files = [picklespath + f for f in os.listdir(picklespath) if not f.startswith(".") if f.endswith(".pickle")]
+files.sort() # The list of files is then sorted alphabetically.
+# Since files are named ${label}-${time} theyy are sorted by simulation and then chronological order
+
+# First loop to sort pickles
+pickdict = {}
+simdict = {}
+for f in files:
+ with open(f, "rb") as file:
+ myctm = pickle.load(file)
+ t = np.fix(
+ myctm.t + 0.5
+ ) # round up to the nearest integer second to avoid truncature errors
+ pickdict[myctm.label, t] = f
+ if myctm.label not in simdict:
+ simdict[myctm.label] = simulation(myctm)
+ if myctm.ts.name == "null":
+ reflabel = myctm.label
+
+for label, t in pickdict:
+ print(label, t)
+ f = pickdict[label, t]
+
+ with open(f, "rb") as file:
+ myctm = pickle.load(file)
+ if label == reflabel:
+ simdict[myctm.label].add_diags(myctm)
+ else:
+ fref = pickdict[reflabel, t]
+ with open(fref, "rb") as file:
+ myctmref = pickle.load(file)
+ simdict[myctm.label].add_diags(myctm, myctmref)
+ if myctm.t < 1.0: # Initial time
+ simdict[myctm.label].addscatter(myctm, "TRC1", "TRC2")
+ simdict[myctm.label].inictm = myctm
+ if myctm.t > 43199.0 and myctm.t < 43201.0: # halftime
+ simdict[myctm.label].addscatter(myctm, "TRC1", "TRC2")
+ simdict[myctm.label].addscatter(myctm, "NOx", "O3")
+ simdict[myctm.label].halfctm = myctm
+ if myctm.t > 86399.0: # last time step
+ simdict[myctm.label].addscatter(myctm, "TRC1", "TRC2")
+ simdict[myctm.label].addscatter(myctm, "NOx", "O3")
+ simdict[myctm.label].finalctm = myctm
+
+# Plot diags
+
+# plot o3min and o3max
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].o3max,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].o3min,
+ linestyle="--",
+ color=colordict[simdict[lab].tsname],
+ )
+plt.legend()
+plt.xlabel("t (s)")
+plt.ylabel(r"$\alpha{}_{\mathrm{O3}}\,\mathrm{(ppb)}$")
+# plt.show()
+plt.savefig(diagspath + "o3minmax.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "o3minmax.pdf", bbox_inches="tight")
+plt.close()
+
+# plot cmin and cmax
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].cmax,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].cmin,
+ linestyle="--",
+ color=colordict[simdict[lab].tsname],
+ )
+plt.xlabel("t (s)")
+plt.legend()
+# plt.show()
+plt.savefig(diagspath + "cminmax.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "cminmax.pdf", bbox_inches="tight")
+plt.close()
+
+# plot testnmin and testnmax
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].testnmax,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].testnmin,
+ linestyle="--",
+ color=colordict[simdict[lab].tsname],
+ )
+plt.xlabel("t (s)")
+plt.legend()
+# plt.show()
+plt.savefig(diagspath + "testnminmax.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "testnminmax.pdf", bbox_inches="tight")
+plt.close()
+
+# plot n2min and n2max
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].n2max,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].n2min,
+ linestyle="--",
+ color=colordict[simdict[lab].tsname],
+ )
+plt.xlabel("t (s)")
+plt.legend()
+# plt.show()
+plt.savefig(diagspath + "n2minmax.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "n2minmax.pdf", bbox_inches="tight")
+plt.close()
+
+# plot max TRC
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].trcmax,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+plt.legend()
+plt.xlabel("t (s)")
+plt.ylabel(r"$\alpha{}_{\mathrm{TRC}}\,\mathrm{(ppb)}$")
+# plt.show()
+plt.savefig(diagspath + "trcmax.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "trcmax.pdf", bbox_inches="tight")
+plt.close()
+
+# plot TRC half volume
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].trchalf,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+plt.legend()
+ylim = np.asarray(plt.gca().get_ylim())
+plt.gca().set_ylim((0.0, ylim[1])) # Force axis minimum to 0
+# plt.show()
+plt.xlabel("t (s)")
+plt.ylabel("Half-volume (nd)")
+plt.savefig(diagspath + "halfvolume.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "halfvolume.pdf", bbox_inches="tight")
+plt.close()
+
+# plot the mass conservation checks
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].nmass,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+plt.legend()
+# plt.show()
+plt.savefig(diagspath + "nmass.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "nmass.pdf", bbox_inches="tight")
+plt.close()
+
+
+# N2 mass
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].n2mass,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+plt.legend()
+# plt.show()
+plt.savefig(diagspath + "n2mass.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "n2mass.pdf", bbox_inches="tight")
+plt.close()
+
+
+# O3 mixing ratio
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ np.asarray(simdict[lab].o3mass) / np.asarray(simdict[lab].mmass) * 1.0e9,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+plt.xlabel("t (s)")
+plt.ylabel(r"Average $\mathrm{O}_3$ mixing ratio (ppb)")
+plt.legend()
+# plt.show()
+plt.savefig(diagspath + "o3mr.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "o3mr.pdf", bbox_inches="tight")
+plt.close()
+
+# NO mixing ratio
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ np.asarray(simdict[lab].nomass) / np.asarray(simdict[lab].mmass) * 1.0e9,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+plt.xlabel("t (s)")
+plt.ylabel(r"Average NO mixing ratio (ppb)")
+plt.legend()
+# plt.show()
+plt.savefig(diagspath + "nomr.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "nomr.pdf", bbox_inches="tight")
+plt.close()
+
+# NO2 mixing ratio
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ np.asarray(simdict[lab].no2mass) / np.asarray(simdict[lab].mmass) * 1.0e9,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+plt.xlabel("t (s)")
+plt.ylabel(r"Average NO2 mixing ratio (ppb)")
+plt.legend()
+# plt.show()
+plt.savefig(diagspath + "no2mr.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "no2mr.pdf", bbox_inches="tight")
+plt.close()
+
+
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].hmass,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+plt.legend()
+# plt.show()
+plt.savefig(diagspath + "hmass.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "hmass.pdf", bbox_inches="tight")
+plt.close()
+
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].cmass,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+plt.legend()
+# plt.show()
+plt.savefig(diagspath + "cmass.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "cmass.pdf", bbox_inches="tight")
+plt.close()
+
+plt.figure()
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].trcmass,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+plt.legend()
+# plt.show()
+plt.savefig(diagspath + "trcmass.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "trcmass.pdf", bbox_inches="tight")
+plt.close()
+
+
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].omass,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+plt.legend()
+# plt.show()
+plt.savefig(diagspath + "omass.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "omass.pdf", bbox_inches="tight")
+plt.close()
+
+
+# plot niters
+plt.figure(figsize=(4.5, 4.0))
+for lab in simdict:
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].niters,
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+
+plt.legend()
+# plt.show()
+plt.xlabel("t (s)")
+plt.savefig(diagspath + "niters.png", dpi=300, bbox_inches="tight")
+plt.savefig(diagspath + "niters.pdf", bbox_inches="tight")
+plt.close()
+
+
+# Plot pdf error with and without Godunov
+
+for plotGodunov in True, False:
+ for sp in accurspecs:
+ plt.figure(figsize=(6.4, 4.8))
+ # This is the default figure size in matplotlib, bigger than the files intended for 3column to obtain more or less uniform fonts
+ for lab in simdict:
+ if lab == reflabel:
+ continue
+ if not plotGodunov:
+ if nicelabel[simdict[lab].tsname] == "Godunov":
+ continue
+ plt.plot(
+ simdict[lab].times,
+ simdict[lab].pdfdist[sp],
+ label=nicelabel[simdict[lab].tsname],
+ color=colordict[simdict[lab].tsname],
+ )
+ plt.legend()
+ plt.xlabel("t (s)")
+ plt.ylabel(r"$\mathcal{E}_1$")
+ if plotGodunov:
+ plt.savefig(
+ diagspath + "pdfdist-" + sp + ".png", dpi=300, bbox_inches="tight"
+ )
+ plt.savefig(diagspath + "pdfdist-" + sp + ".pdf", bbox_inches="tight")
+ else:
+ plt.savefig(
+ diagspath + "pdfdist-noGod-" + sp + ".png", dpi=300, bbox_inches="tight"
+ )
+ plt.savefig(diagspath + "pdfdist-noGod-" + sp + ".pdf", bbox_inches="tight")
+
+ plt.close()
+
+
+# Plot scatter
+
+xvec = np.linspace(0.1, 0.9, 100)
+refscat_trc = Scatter(
+ xvec * 1.0e9,
+ np.asarray([0.9 - 0.8 * x ** 2 for x in xvec]) * 1.0e9,
+ unitx=None,
+ unity=None,
+ spx="TRC1",
+ spy="TRC2",
+)
+
+for lab in simdict:
+ if lab == reflabel:
+ continue
+ for key in simdict[lab].scatters:
+ sp1, sp2, t = key
+ scat = simdict[lab].scatters[key]
+ if sp1 == "TRC1" and sp2 == "TRC2":
+ scat.plot(
+ color=colordict[simdict[lab].tsname], refscat=refscat_trc, simlabel=lab
+ )
+ else:
+ scat.plot(
+ color=colordict[simdict[lab].tsname],
+ refscat=simdict[reflabel].scatters[key],
+ simlabel=lab,
+ )
+
+# Plot complete PDFs
+for lab in simdict:
+ if lab == reflabel:
+ continue
+ plotpdfs(simdict[lab].finalctm, simdict[reflabel].finalctm, prefix="finalpdf")
+ plotpdfs(simdict[lab].halfctm, simdict[reflabel].halfctm, prefix="halfpdf")
+ plotpdfs(
+ simdict[lab].finalctm,
+ simdict[reflabel].finalctm,
+ prefix="finalpdf_fancy_",
+ fancy=True,
+ )
+ plotpdfs(
+ simdict[lab].halfctm,
+ simdict[reflabel].halfctm,
+ prefix="halfpdf_fancy_",
+ fancy=True,
+ )
+
+# Make some maps
+# Map the flow
+
+for lab in simdict:
+ if lab == reflabel:
+ continue
+ simdict[lab].inictm.MassFluxFunc.plotpsi(
+ 0.0, 1.0e5, 0.0, 1.0e5, simdict[lab].inictm.t, fname=diagspath + "psi_leveque"
+ )
+ break
+
+# Initialisations for maps
+sdict = {}
+sdict["O3"] = "ppb"
+sdict["NO"] = "ppb"
+sdict["NO2"] = "ppb"
+sdict["TRC"] = "ppb"
+sdict["TRC1"] = "ppb"
+sdict["TRC2"] = "ppb"
+sdict["CO"] = "ppb"
+sdict["testN"] = "ppb"
+sdict["testC"] = "ppb"
+
+for lab in simdict:
+ ctmtup = (simdict[lab].inictm, simdict[lab].halfctm, simdict[lab].finalctm)
+
+ for chemtm in ctmtup:
+ for spec in sdict:
+ if spec == "testN":
+ chemtm.makeplotchem(
+ sp=("NO", "NO2", "HNO3", "TRCb"),
+ simlabel=lab,
+ unit=sdict[spec],
+ spname="testN",
+ )
+ elif spec == "testC":
+ chemtm.makeplotchem(
+ sp=("CO", "CO2"), simlabel=lab, unit=sdict[spec], spname="testC"
+ )
+ else:
+ chemtm.makeplotchem(sp=spec, simlabel=lab, unit=sdict[spec])
+
+# Accuracy diagnostics
+
+norm1dict = {}
+norm2dict = {}
+e1dict = {}
+
+# Norm-2
+for lab in simdict:
+ if lab == reflabel:
+ continue
+
+ for sp in accurspecs:
+ diff = simdict[lab].finalctm.activeConc_ppb(sp) - simdict[
+ reflabel
+ ].finalctm.activeConc_ppb(sp)
+ ref = simdict[reflabel].finalctm.activeConc_ppb(sp)
+ norm2error = np.sqrt(np.sum(diff * diff)) / np.sqrt(
+ np.sum(ref * ref)
+ ) # RMSE normalized by reference value
+ norm2dict[nicelabel[simdict[lab].tsname], sp] = norm2error
+
+latex2 = dict2latex(norm2dict, reflabel, norm="norm2")
+print(latex2)
+with open(diagspath + "norm2.tex", "w") as text_file:
+ text_file.write(latex2)
+
+# Norm-1
+for lab in simdict:
+ if lab == reflabel:
+ continue
+ for sp in accurspecs:
+ diff = np.abs(
+ simdict[lab].finalctm.activeConc_ppb(sp)
+ - simdict[reflabel].finalctm.activeConc_ppb(sp)
+ )
+ ref = simdict[reflabel].finalctm.activeConc_ppb(sp)
+ norm1error = np.sum(diff) / np.sum(ref) # RMSE normalized by reference value
+ norm1dict[nicelabel[simdict[lab].tsname], sp] = norm1error
+
+latex1 = dict2latex(norm1dict, reflabel, norm="norm1")
+print(latex1)
+with open(diagspath + "norm1.tex", "w") as text_file:
+ text_file.write(latex1)
+
+# e1
+for lab in simdict:
+ if lab == reflabel:
+ continue
+ for sp in accurspecs:
+ e1dict[nicelabel[simdict[lab].tsname], sp] = pdfdist(
+ sp, simdict[lab].finalctm, simdict[reflabel].finalctm
+ )
+
+latex_e1 = dict2latex(e1dict, reflabel, norm="e1")
+print(latex_e1)
+with open(diagspath + "e1.tex", "w") as text_file:
+ text_file.write(latex_e1)
+
+
+# Screen outputs for initial / final mass for C, H, N
+for lab in simdict:
+ print(
+ nicelabel[simdict[lab].tsname],
+ "C & & ",
+ simdict[lab].cmass[0],
+ "& ",
+ simdict[lab].cmass[-1],
+ "&",
+ np.abs(simdict[lab].cmass[-1] - simdict[lab].cmass[0]) / simdict[lab].cmass[0],
+ )
+ print(
+ nicelabel[simdict[lab].tsname],
+ "N & & ",
+ simdict[lab].nmass[0],
+ "& ",
+ simdict[lab].nmass[-1],
+ "&",
+ np.abs(simdict[lab].nmass[-1] - simdict[lab].nmass[0]) / simdict[lab].nmass[0],
+ )
+ print(
+ nicelabel[simdict[lab].tsname],
+ "H & & ",
+ simdict[lab].hmass[0],
+ "& ",
+ simdict[lab].hmass[-1],
+ "&",
+ np.abs(simdict[lab].hmass[-1] - simdict[lab].hmass[0]) / simdict[lab].hmass[0],
+ )
+ print(
+ nicelabel[simdict[lab].tsname],
+ "TRC & & ",
+ simdict[lab].trcmass[0],
+ "& ",
+ simdict[lab].trcmass[-1],
+ "&",
+ np.abs(simdict[lab].trcmass[-1] - simdict[lab].trcmass[0])
+ / simdict[lab].trcmass[0],
+ )
+
+# if CI test
+if(args.CIcheck):
+ CItolerance = 0.01
+ passed = True
+ targetErrorDict = {}
+ targetErrorDict['Godunov', 'TRC'] = 1.23
+ targetErrorDict['PPM', 'TRC'] = 0.986
+ targetErrorDict['PPM+W', 'TRC'] = 0.946
+ targetErrorDict['Van Leer', 'TRC'] = 1.05
+ targetErrorDict['Walcek', 'TRC'] = 0.961
+ targetErrorDict['Godunov', 'O3'] = 1.53
+ targetErrorDict['PPM', 'O3'] = 0.979
+ targetErrorDict['PPM+W', 'O3'] = 0.897
+ targetErrorDict['Van Leer', 'O3'] = 1.09
+ targetErrorDict['Walcek', 'O3'] = 0.914
+ for sp in ('TRC', 'O3'):
+ for scheme in ('Godunov', 'PPM', 'PPM+W', 'Van Leer', 'Walcek', 'Godunov'):
+ if (norm1dict[scheme, sp] < targetErrorDict[scheme, sp] + 0.01):
+ print('CI test', scheme, sp, norm1dict[scheme, sp], '<=', targetErrorDict[scheme, sp],'PASSED')
+ else:
+ passed = False
+ print('CI test', scheme, sp, norm1dict[scheme, sp], '> ', targetErrorDict[scheme, sp],'FAILED')
+
+ if(passed):
+ print('CI test sussessful. Leaving')
+ sys.exit(0)
+ else:
+ print('ERROR: CI test failed')
+ sys.exit(1)
+
diff --git a/2dtestcases.py b/2dtestcases.py
index 089813b34b12916dac3e69800894f1b1616eb9fd..334aa9a48465e0d4281f565ccdc0b0a93a6cc22d 100755
--- a/2dtestcases.py
+++ b/2dtestcases.py
@@ -106,7 +106,7 @@ parser.add_argument(
parser.add_argument(
"--picklespath",
type=str,
- default="pickles/",
+ default="artifacts/",
help="Relative or absolute path to store the pickles",
)
@@ -611,7 +611,7 @@ print("init TRC mass", myctm.mass("TRC"))
if makedumps:
file = open(
- picklespath + "{0}-{1:06d}".format(myctm.label, int(myctm.t + 0.5)), "wb"
+ picklespath + "{0}-{1:06d}.pickle".format(myctm.label, int(myctm.t + 0.5)), "wb"
)
pickle.dump(myctm, file)
file.close()
@@ -644,7 +644,7 @@ for i in range(nt):
# Dump myctm
if makedumps:
print(myctm.t, int(myctm.t + 0.5))
- file = open(picklespath + flabel, "wb")
+ file = open(picklespath + flabel + '.pickle', "wb")
pickle.dump(myctm, file)
file.close()