TL: added imex stability scripts for thomas

Author: tlunet

qmat/playgrounds/tibo/__init__.py

@@ -2,4 +2,6 @@
 - :class:`orthogonalPolynomials` : generate orthogonal polynomial values from any distribution.
 - :class:`lorenz` : application example of the generic solvers to solve the Lorenz equations.
 - :class:`imex` : starting development for the IMEX generic solvers.
+- :class:`imexStabilityRK` : investigating IMEX stability for RK methods
+- :class:`imexStabilitySDC` : investigating IMEX stability for SDC methods
 """

qmat/playgrounds/tibo/imexStabilityRK.py

@@ -0,0 +1,81 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Script investigating IMEX stability for RK methods
+"""
+import numpy as np
+import scipy.optimize as sco
+import matplotlib.pyplot as plt
+
+from qmat.qcoeff.butcher import RK_SCHEMES
+from qmat.solvers.dahlquist import DahlquistIMEX
+
+
+# Script parameters
+implScheme = "ARK443ESDIRK"
+explScheme = "ARK443ERK"
+stepUpdate = False
+
+
+# Script execution
+lamE = np.linspace(0, 6, num=501)
+ratio = np.logspace(-3, 2, num=301)
+zI = -ratio[None, :]*lamE[:, None]
+zE = 1j*lamE[:, None]
+
+problem = DahlquistIMEX(zI, zE)
+
+schemeI = RK_SCHEMES[implScheme]()
+schemeE = RK_SCHEMES[explScheme]()
+
+QI = schemeI.Q
+wI = schemeI.weights if stepUpdate else None
+QE = schemeE.Q
+wE = schemeE.weights if stepUpdate else None
+
+uNum = problem.solve(QI, wI, QE, wE)
+
+u1 = uNum[-1]
+stab = np.abs(u1)
+stab = np.clip(stab, 0, 1.2) # clip to ignore unstable area
+error = np.abs(u1 - np.exp(zI+zE))
+
+
+plt.figure("imex stability")
+plt.clf()
+
+coords = (ratio, lamE)
+plt.contourf(*coords, stab, levels = [0, 0.2, 0.4, 0.6, 0.8, 1, 1.2])
+plt.colorbar()
+plt.contour(*coords, stab, levels=[1], colors="black")
+plt.contour(*coords, error, levels=[1], colors="red", linestyles=":")
+plt.contour(*coords, error, levels=[1e-1], colors="orange", linestyles="-.")
+plt.contour(*coords, error, levels=[1e-2], colors="gray", linestyles="--")
+plt.grid(True)
+plt.xscale('log')
+plt.ylabel(r"$\lambda_E \Delta t$")
+plt.xlabel(r"advection $\quad\leftarrow\quad\frac{-\lambda_I}{\lambda_E}\quad\rightarrow\quad$ diffusion", fontsize=20)
+plt.tight_layout()
+
+
+def imStab(x):
+    return np.abs(DahlquistIMEX([0], [x*1j]).solve(QI, wI, QE, wE)[-1]) - 1
+
+try:
+    sol = sco.bisect(imStab, 1e-1, 1e2, xtol=1e-16)
+    print(f"CFL max [RK] : {sol}")
+except:
+    pass
+
+plt.figure("stability on imaginary axis")
+plt.clf()
+plt.grid(True)
+x = np.linspace(0, 6, num=501)
+plt.plot(x, [imStab(s)[0] for s in x], label="RK")
+plt.ylim(-1, 0.5)
+plt.ylabel("over-amplification")
+plt.xlabel(r"$\lambda_E \Delta t$")
+plt.legend()
+plt.tight_layout()
+
+plt.show()

qmat/playgrounds/tibo/imexStabilitySDC.py

@@ -0,0 +1,91 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Script investigating IMEX stability for SDC methods
+"""
+import numpy as np
+import scipy.optimize as sco
+import matplotlib.pyplot as plt
+
+from qmat.qcoeff.collocation import Collocation
+from qmat.qdelta import QDELTA_GENERATORS
+
+from qmat.solvers.dahlquist import DahlquistIMEX
+
+
+# Script parameters
+nNodes = 4
+nSweeps = 4
+stepUpdate = False
+implSweep = "MIN-SR-FLEX"
+explSweep = "PIC"
+
+
+# Script execution
+lamE = np.linspace(0, 6, num=501)
+ratio = np.logspace(-3, 2, num=301)
+zI = -ratio[None, :]*lamE[:, None]
+zE = 1j*lamE[:, None]
+
+problem = DahlquistIMEX(zI, zE)
+
+coll = Collocation(nNodes=nNodes, nodeType="LEGENDRE", quadType="RADAU-RIGHT")
+
+genI = QDELTA_GENERATORS[implSweep](qGen=coll)
+genE = QDELTA_GENERATORS[explSweep](qGen=coll)
+
+sweeps = [k+1 for k in range(nSweeps)]
+
+uNum = problem.solveSDC(
+    coll.Q, coll.weights if stepUpdate else None,
+    genI.genCoeffs(k=sweeps), genE.genCoeffs(k=sweeps),
+    nSweeps=nSweeps)
+
+u1 = uNum[-1]
+stab = np.abs(u1)
+stab = np.clip(stab, 0, 1.2) # clip to ignore unstable area
+error = np.abs(u1 - np.exp(zI+zE))
+
+
+plt.figure("imex stability")
+plt.clf()
+
+coords = (ratio, lamE)
+plt.contourf(*coords, stab, levels = [0, 0.2, 0.4, 0.6, 0.8, 1, 1.2])
+plt.colorbar()
+plt.contour(*coords, stab, levels=[1], colors="black")
+plt.contour(*coords, error, levels=[1], colors="red", linestyles=":")
+plt.contour(*coords, error, levels=[1e-1], colors="orange", linestyles="-.")
+plt.contour(*coords, error, levels=[1e-2], colors="gray", linestyles="--")
+plt.grid(True)
+plt.xscale('log')
+plt.ylabel(r"$\lambda_E \Delta t$")
+plt.xlabel(r"advection $\quad\leftarrow\quad\frac{-\lambda_I}{\lambda_E}\quad\rightarrow\quad$ diffusion", fontsize=20)
+plt.tight_layout()
+
+
+def imStab(x):
+    uSDC = DahlquistIMEX([0], [x*1j]).solveSDC(
+        coll.Q, coll.weights if stepUpdate else None,
+        genI.genCoeffs(k=sweeps), genE.genCoeffs(k=sweeps),
+        nSweeps=nSweeps)
+    return np.abs(uSDC[-1]) - 1
+
+try:
+    sol = sco.bisect(imStab, 1e-1, 1e2, xtol=1e-16)
+    print(f"CFL max [SDC] : {sol}")
+except:
+    pass
+
+plt.figure("stability on imaginary axis")
+plt.clf()
+plt.grid(True)
+x = np.linspace(0, 6, num=501)
+plt.plot(x, [imStab(s)[0] for s in x], label="RK")
+plt.ylim(-1, 0.5)
+plt.ylabel("over-amplification")
+plt.xlabel(r"$\lambda_E \Delta t$")
+plt.legend()
+plt.tight_layout()
+
+plt.show()