TL: playground update

Author: tlunet

docs/notebooks.md

@@ -7,9 +7,9 @@ All tutorials are written in jupyter notebooks, that can be :
 - read using the [online documentation](https://qmat.readthedocs.io/en/latest/notebooks.html)
 - downloaded from the [notebook folder](https://github.com/Parallel-in-Time/qmat/tree/main/docs/notebooks) and played with
 
-> 📋 _Table of content_ :
+📋 **Table of content**
 
-1. **Basic usage** : how to generate and use basic $Q$-coefficients and $Q_\Delta$ approximations, through a step-by-step tutorial going from generic Runge-Kutta methods to SDC for simple problems.
+1. **Basic usage tutorials** : how to generate and use basic $Q$-coefficients and $Q_\Delta$ approximations, through a step-by-step tutorial going from generic Runge-Kutta methods to SDC for simple problems.
 2. **Advanced tutorials** : additional features or `qmat` to go deeper into time-integration (Node-to-Node formulation, use for non-linear problems, $\phi$-SDC, ...)
 3. **Components usage** : how to use the main utility modules, like {py:mod}`qmat.lagrange`, etc ...
 

qmat/__init__.py

@@ -9,7 +9,7 @@
 **Secondary sub-packages** 🍭
 
 - :class:`solvers` : implementations of time-integration solvers that make use of `qmat`-generated coefficients
-- :class:`playgrounds`: **non-tested but documented** codes with experiments or small applications with `qmat`
+- :class:`playgrounds`: **non-tested** codes with experiments or small applications with `qmat`
 
 **Utility modules** ⚙️
 

qmat/playgrounds/tibo/__init__.py

@@ -1,3 +1,6 @@
+"""
+First tentative for a dedicated IMEX solver class
+"""
 import numpy as np
 
 from qmat.solvers.dahlquist import DahlquistIMEX
@@ -9,7 +12,9 @@ class DiffOpIMEX(DiffOp):
     Base class for an IMEX differential operator
     """
 
-    def evalF2(self, u:np.ndarray, t:float, out:np.ndarray):
+    evalF_IM = DiffOp.evalF
+
+    def evalF_EX(self, u:np.ndarray, t:float, out:np.ndarray):
         """
         Evaluate f_EX(u,t) and store the result into out
         """
@@ -26,9 +31,11 @@ def __init__(self, diffOp, tEnd=1, nSteps=1, t0=0):
             f"DiffOpIMEX object is required for diffOp argument, not {diffOp}"
         super().__init__(diffOp, tEnd, nSteps, t0)
 
+    def evalF_IM(self, u:np.ndarray, t:float, out:np.ndarray):
+        self.diffOp.evalF_IM(u, t, out)
 
-    def evalF2(self, u:np.ndarray, t:float, out:np.ndarray):
-        self.diffOp.evalF2(u, t, out)
+    def evalF_EX(self, u:np.ndarray, t:float, out:np.ndarray):
+        self.diffOp.evalF_EX(u, t, out)
 
     def solve(self, QI, wI, QE, wE, uNum=None):
         nNodes, QI, wI, QE, wE, useWeights = DahlquistIMEX.checkCoeff(QI, wI, QE, wE)
@@ -74,8 +81,8 @@ def solve(self, QI, wI, QE, wE, uNum=None):
                     np.copyto(uNode, rhs)
 
                 # evalF on current stage
-                self.evalF(u=uNode, t=tNode, out=fEvals[m])
-                self.evalF2(u=uNode, t=tNode, out=fEvals2[m])
+                self.evalF_IM(u=uNode, t=tNode, out=fEvals[m])
+                self.evalF_EX(u=uNode, t=tNode, out=fEvals2[m])
 
             # step update (if not, uNum[i+1] is already the last stage)
             if useWeights:

qmat/playgrounds/tibo/imex.py

@@ -2,6 +2,11 @@
 # -*- coding: utf-8 -*-
 """
 Script investigating IMEX stability for advection-diffusion solved with SDC
+
+.. literalinclude:: /../qmat/playgrounds/tibo/imexStabilityAdvDiffSDC.py
+   :language: python
+   :linenos:
+   :lines: 11-
 """
 import numpy as np
 import matplotlib.pyplot as plt

qmat/playgrounds/tibo/imexStabilityAdvDiffSDC.py

@@ -2,6 +2,11 @@
 # -*- coding: utf-8 -*-
 """
 Script investigating IMEX stability for RK methods
+
+.. literalinclude:: /../qmat/playgrounds/tibo/imexStabilityRK.py
+   :language: python
+   :linenos:
+   :lines: 11-
 """
 import numpy as np
 import scipy.optimize as sco

qmat/playgrounds/tibo/imexStabilityRK.py

@@ -2,6 +2,11 @@
 # -*- coding: utf-8 -*-
 """
 Script investigating IMEX stability for SDC methods
+
+.. literalinclude:: /../qmat/playgrounds/tibo/imexStabilitySDC.py
+   :language: python
+   :linenos:
+   :lines: 11-
 """
 import numpy as np
 import scipy.optimize as sco

qmat/playgrounds/tibo/imexStabilitySDC.py

@@ -7,7 +7,7 @@
 .. literalinclude:: /../qmat/playgrounds/tibo/lorenz.py
    :language: python
    :linenos:
-   :lines: 11-
+   :lines: 12-
 """
 import numpy as np
 import matplotlib.pyplot as plt

qmat/playgrounds/tibo/lorenz.py

@@ -0,0 +1,56 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Compute accuracy and stability of a monolitic SDC scheme
+
+.. literalinclude:: /../qmat/playgrounds/tibo/stability.py
+   :language: python
+   :linenos:
+   :lines: 11-
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+from qmat.qcoeff.collocation import Collocation
+from qmat.qdelta import QDELTA_GENERATORS
+
+from qmat.solvers.dahlquist import Dahlquist
+
+# Script parameters
+nNodes = 4
+nSweeps = 1
+sweepType = "SOE"
+
+# Script execution
+re = np.linspace(-4.5, 1, num=200)
+im = np.linspace(0, 5.0, num=201)
+lam = re[None, :] + 1j*im[:, None]
+
+problem = Dahlquist(lam)
+coll = Collocation(nNodes=nNodes, nodeType="LEGENDRE", quadType="RADAU-RIGHT")
+approx = QDELTA_GENERATORS[sweepType](qGen=coll)
+
+sweeps = [k+1 for k in range(nSweeps)]
+
+uNum = problem.solveSDC(
+    coll.Q, None, approx.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(lam))
+
+plt.figure(f"{sweepType}_M{nNodes}_K{nSweeps}")
+coords = (re, im)
+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.ylabel(r"$Im(\lambda)$")
+plt.xlabel(r"$Re(\lambda)$")
+plt.gca().set_aspect('equal', 'box')
+plt.tight_layout()