TB: added two embedded methods

* Added two embedded methods

* Got rid of separate class for embedded RK methods

* "Embedded Q" but with bugs

* Removed embedded stuff from Q and added default secondary order

* Implemented Thibaut's suggestions

qmat/qcoeff/__init__.py

@@ -29,13 +29,24 @@ def Q(self):
     def weights(self):
         raise NotImplementedError("mouahahah")
 
+    @property
+    def weightsEmbedded(self):
+        """
+        These weights can be used to construct a secondary lower order method from the same stages.
+        """
+        raise NotImplementedError("No embedded weights implemented for {type(self).__name__}")
+
     @property
     def nNodes(self):
         return self.nodes.size
 
+    @property
+    def rightIsNode(self):
+        return self.nodes[-1] == 1.
+
     @property
     def T(self):
-        """Transfert matrix from zero-to-nodes to node-to-node"""
+        """Transfer matrix from zero-to-nodes to node-to-node"""
         M = self.Q.shape[0]
         T = np.eye(M)
         T[1:,:-1][np.diag_indices(M-1)] = -1
@@ -51,7 +62,7 @@ def S(self):
 
     @property
     def Tinv(self):
-        """Transfert matrix from node-to-node to zero-to-node"""
+        """Transfer matrix from node-to-node to zero-to-node"""
         M = self.Q.shape[0]
         return np.tri(M)
 
@@ -60,22 +71,31 @@ def hCoeffs(self):
         approx = LagrangeApproximation(self.nodes)
         return approx.getInterpolationMatrix([1]).ravel()
 
-    def genCoeffs(self, withS=False, hCoeffs=False):
+    def genCoeffs(self, withS=False, hCoeffs=False, embedded=False):
         out = [self.nodes, self.weights, self.Q]
+
+        if embedded:
+            out[1] = np.vstack([out[1], self.weightsEmbedded])
         if withS:
             out.append(self.S)
         if hCoeffs:
             out.append(self.hCoeffs)
         return out
 
-
     @property
     def order(self):
         raise NotImplementedError("mouahahah")
 
-    def solveDahlquist(self, lam, u0, T, nSteps):
+    @property
+    def orderEmbedded(self):
+        return self.order - 1
+
+    def solveDahlquist(self, lam, u0, T, nSteps, useEmbeddedWeights=False):
         nodes, weights, Q = self.nodes, self.weights, self.Q
 
+        if useEmbeddedWeights:
+            weights = self.weightsEmbedded
+
         uNum = np.zeros(nSteps+1, dtype=complex)
         uNum[0] = u0
 
@@ -88,9 +108,9 @@ def solveDahlquist(self, lam, u0, T, nSteps):
 
         return uNum
 
-    def errorDahlquist(self, lam, u0, T, nSteps, uNum=None):
+    def errorDahlquist(self, lam, u0, T, nSteps, uNum=None, useEmbeddedWeights=False):
         if uNum is None:
-            uNum = self.solveDahlquist(lam, u0, T, nSteps)
+            uNum = self.solveDahlquist(lam, u0, T, nSteps, useEmbeddedWeights=useEmbeddedWeights)
         times = np.linspace(0, T, nSteps+1)
         uExact = u0 * np.exp(lam*times)
         return np.linalg.norm(uNum-uExact, ord=np.inf)
@@ -116,18 +136,18 @@ def register(cls:QGenerator)->QGenerator:
             cls(**params)
         except:
             raise TypeError(
-                f"{cls.__name__} could not be instanciated with DEFAULT_PARAMS")
+                f"{cls.__name__} could not be instantiated with DEFAULT_PARAMS")
     # Store class (and aliases)
     storeClass(cls, Q_GENERATORS)
     return cls
 
-def genQCoeffs(qType, withS=False, hCoeffs=False, **params):
+def genQCoeffs(qType, withS=False, hCoeffs=False, embedded=False, **params):
     try:
         Generator = Q_GENERATORS[qType]
     except KeyError:
-        raise ValueError(f"qType={qType} is not available")
+        raise ValueError(f"{qType=!r} is not available")
     gen = Generator(**params)
-    return gen.genCoeffs(withS, hCoeffs)
+    return gen.genCoeffs(withS, hCoeffs, embedded)
 
 
 # Import all local submodules

qmat/qcoeff/butcher.py

@@ -17,20 +17,28 @@
 import numpy as np
 
 from qmat.qcoeff import QGenerator, register
-from qmat.utils import storeClass
+from qmat.utils import storeClass, checkOverriding
 
 
 class RK(QGenerator):
     A = None
     b = None
     c = None
+    b2 = None  # for embedded methods
 
     @property
     def nodes(self): return self.c
 
     @property
     def weights(self): return self.b
 
+    @property
+    def weightsEmbedded(self):
+        if self.b2 is None:
+            raise NotImplementedError(f'Kindly direct your request for an embedded version of {type(self).__name__!r} to the Mermathematicians on Europa.')
+        else:
+            return self.b2
+
     @property
     def Q(self): return self.A
 
@@ -41,8 +49,12 @@ def checkAndStore(cls:RK)->RK:
     cls.A = np.array(cls.A, dtype=float)
     cls.b = np.array(cls.b, dtype=float)
     cls.c = np.array(cls.c, dtype=float)
-    assert cls.b.size == cls.c.size, \
-        f"b (size {cls.b.size}) and c (size {cls.b.size})" + \
+
+    if cls.b2 is not None:
+        cls.b2 = np.array(cls.b2, dtype=float)
+
+    assert cls.b.shape[-1] == cls.c.size, \
+        f"b (size {cls.b.shape[-1]}) and c (size {cls.c.size})" + \
         f" have not the same size in {cls.__name__}"
     assert cls.A.shape == (cls.c.size, cls.c.size), \
         f"A (shape {cls.A.shape}) and c (size {cls.b.size})" + \
@@ -132,7 +144,7 @@ def order(self): return 1
 
 @registerRK
 class RK2(RK):
-    """Classical Runge Kutta method of order 2 (cf Wikipedia)"""
+    """Classical Runge-Kutta method of order 2 (cf Wikipedia)"""
     aliases = ["ERK2"]
     A = [[0, 0],
          [1/2, 0]]
@@ -300,3 +312,96 @@ class SDIRK54(RK):
 
     @property
     def order(self): return 4
+
+@registerRK
+class HeunEuler(RK):
+    """
+    Second order explicit embedded Runge-Kutta method.
+    """
+    A = [[0, 0],
+         [1, 0]]
+    b = [0.5, 0.5]
+    c = [0, 1]
+    b2 = [1, 0]
+
+    @property
+    def order(self): return 2
+
+@registerRK
+class CashKarp(RK):
+    """
+    Fifth order explicit embedded Runge-Kutta. See [here](https://doi.org/10.1145/79505.79507).
+    """
+    c = [0, 0.2, 0.3, 0.6, 1.0, 7.0 / 8.0]
+    b = [37.0 / 378.0, 0.0, 250.0 / 621.0, 125.0 / 594.0, 0.0, 512.0 / 1771.0]
+    b2 = [2825.0 / 27648.0, 0.0, 18575.0 / 48384.0, 13525.0 / 55296.0, 277.0 / 14336.0, 1.0 / 4.0]
+    A = np.zeros((6, 6))
+    A[1, 0] = 1.0 / 5.0
+    A[2, :2] = [3.0 / 40.0, 9.0 / 40.0]
+    A[3, :3] = [0.3, -0.9, 1.2]
+    A[4, :4] = [-11.0 / 54.0, 5.0 / 2.0, -70.0 / 27.0, 35.0 / 27.0]
+    A[5, :5] = [1631.0 / 55296.0, 175.0 / 512.0, 575.0 / 13824.0, 44275.0 / 110592.0, 253.0 / 4096.0]
+
+    @property
+    def order(self): return 5
+
+    CONV_TEST_NSTEPS = [32, 64, 128]
+
+@registerRK
+class ESDIRK53(RK):
+    """
+    A-stable embedded RK pair of orders 5 and 3, ESDIRK5(3)6L[2]SA.
+    Taken from [here](https://ntrs.nasa.gov/citations/20160005923).
+    """
+    c = [0, 4024571134387.0 / 7237035672548.0, 14228244952610.0 / 13832614967709.0, 1.0 / 10.0, 3.0 / 50.0, 1.0]
+    A = np.zeros((6, 6))
+    A[1, :2] = [3282482714977.0 / 11805205429139.0, 3282482714977.0 / 11805205429139.0]
+    A[2, :3] = [
+        606638434273.0 / 1934588254988,
+        2719561380667.0 / 6223645057524,
+        3282482714977.0 / 11805205429139.0,
+    ]
+    A[3, :4] = [
+        -651839358321.0 / 6893317340882,
+        -1510159624805.0 / 11312503783159,
+        235043282255.0 / 4700683032009.0,
+        3282482714977.0 / 11805205429139.0,
+    ]
+    A[4, :5] = [
+        -5266892529762.0 / 23715740857879,
+        -1007523679375.0 / 10375683364751,
+        521543607658.0 / 16698046240053.0,
+        514935039541.0 / 7366641897523.0,
+        3282482714977.0 / 11805205429139.0,
+    ]
+    A[5, :] = [
+        -6225479754948.0 / 6925873918471,
+        6894665360202.0 / 11185215031699,
+        -2508324082331.0 / 20512393166649,
+        -7289596211309.0 / 4653106810017.0,
+        39811658682819.0 / 14781729060964.0,
+        3282482714977.0 / 11805205429139,
+    ]
+
+    b = [
+            -6225479754948.0 / 6925873918471,
+            6894665360202.0 / 11185215031699.0,
+            -2508324082331.0 / 20512393166649,
+            -7289596211309.0 / 4653106810017,
+            39811658682819.0 / 14781729060964.0,
+            3282482714977.0 / 11805205429139,
+        ]
+    b2 = [
+           -2512930284403.0 / 5616797563683,
+           5849584892053.0 / 8244045029872,
+           -718651703996.0 / 6000050726475.0,
+           -18982822128277.0 / 13735826808854.0,
+           23127941173280.0 / 11608435116569.0,
+           2847520232427.0 / 11515777524847.0,
+        ]
+
+    @property
+    def orderEmbedded(self): return 3
+
+    @property
+    def order(self): return 5

tests/test_qcoeff/test_base.py

@@ -24,7 +24,7 @@ def testGeneration(name):
     assert n1.size == w1.size, \
         f"nodes and weights for {name} don't have the same size : {n1} / {w1}"
     assert Q1.shape == (n1.size, n1.size), \
-        f"Q for {name} has unconsistent shape : {Q1.shape}"
+        f"Q for {name} has inconsistent shape : {Q1.shape}"
     assert gen.nNodes == n1.size, \
         f"nNodes property from {name} is not equal to node size !"
 
@@ -51,7 +51,7 @@ def testAdditionalCoeffs(name):
         f"hCoeffs for {name} are not np.ndarray but {type(h1)}"
 
     assert S1.shape == (nodes.size, nodes.size), \
-        f"S for {name} has unconsistent shape : {S1.shape}"
+        f"S for {name} has inconsistent shape : {S1.shape}"
     assert h1.ndim == 1, f"hCoeffs for {name} are not 1D : {h1}"
     assert h1.size == nodes.size, \
         f"hCoeffs for {name} has inconsistent size : {h1.size}"
@@ -67,6 +67,18 @@ def testAdditionalCoeffs(name):
         f"OOP hCoeffs {h1} and PP hCoeffs {h2} are not equals for {name}"
 
 
+    try:
+        try:
+            _, b, _  = genQCoeffs(name, embedded=True)
+        except TypeError:
+            _, b, _  = genQCoeffs(name, embedded=True, **GENERATORS[name].DEFAULT_PARAMS)
+
+        assert type(b) == np.ndarray
+        assert b.ndim == 2
+    except NotImplementedError:
+        pass
+
+
 @pytest.mark.parametrize("name", GENERATORS.keys())
 def testS(name):
     gen = GENERATORS[name].getInstance()

tests/test_qcoeff/test_convergence.py

@@ -7,21 +7,23 @@
 
 SCHEMES = getClasses(Q_GENERATORS)
 
-def nStepsForTest(scheme):
+def nStepsForTest(scheme, useEmbeddedWeights=False):
     try:
         nSteps = scheme.CONV_TEST_NSTEPS
     except AttributeError:
-        if scheme.order == 1:
+        expectedOrder = scheme.orderEmbedded if useEmbeddedWeights else scheme.order
+
+        if expectedOrder == 1:
             nSteps = [64, 128, 256]
-        elif scheme.order == 2:
+        elif expectedOrder == 2:
             nSteps = [32, 64, 128]
-        elif scheme.order == 3:
+        elif expectedOrder == 3:
             nSteps = [16, 32, 64]
-        elif scheme.order in [4, 5]:
+        elif expectedOrder in [4, 5]:
             nSteps = [8, 16, 32]
-        elif scheme.order in [6, 7]:
+        elif expectedOrder in [6, 7]:
             nSteps = [4, 8, 16]
-        elif scheme.order in [8, 9]:
+        elif expectedOrder in [8, 9]:
             nSteps = [2, 4, 8]
         else:
             nSteps = [1, 2, 4]  # default value (very high order methods)
@@ -33,28 +35,47 @@ def nStepsForTest(scheme):
 
 
 @pytest.mark.parametrize("scheme", SCHEMES.keys())
-def testDahlquist(scheme):
+@pytest.mark.parametrize("useEmbeddedWeights", [True, False])
+def testDahlquist(scheme, useEmbeddedWeights):
     gen = SCHEMES[scheme].getInstance()
-    nSteps = nStepsForTest(gen)
-    err = [gen.errorDahlquist(lam, u0, T, nS) for nS in nSteps]
+
+    if useEmbeddedWeights:
+        try:
+            gen.weightsEmbedded
+        except NotImplementedError:
+            return None
+
+    expectedOrder = gen.orderEmbedded if useEmbeddedWeights else gen.order
+    nSteps = nStepsForTest(gen, useEmbeddedWeights)
+    err = [gen.errorDahlquist(lam, u0, T, nS, useEmbeddedWeights=useEmbeddedWeights) for nS in nSteps]
     order, rmse = numericalOrder(nSteps, err)
     assert rmse < 0.02, f"rmse to high ({rmse}) for {scheme}"
-    assert abs(order-gen.order) < 0.1, f"wrong numerical order ({order}) for {scheme}"
+    assert abs(order-expectedOrder) < 0.1, f"Expected order {expectedOrder:.2f}, but got {order:.2f} for {scheme}"
 
 
 @pytest.mark.parametrize("quadType", QUAD_TYPES)
 @pytest.mark.parametrize("nodesType", NODE_TYPES)
 @pytest.mark.parametrize("nNodes", [2, 3, 4])
-def testDahlquistCollocation(nNodes, nodesType, quadType):
+def testDahlquistCollocation(nNodes, nodesType, quadType, useEmbeddedWeights=False):
     gen = SCHEMES["Collocation"](nNodes, nodesType, quadType)
+    if useEmbeddedWeights:
+        try:
+            gen.weightsEmbedded
+        except NotImplementedError:
+            return None
     scheme = f"Collocation({nNodes}, {nodesType}, {quadType})"
-    nSteps = nStepsForTest(gen)
+    nSteps = nStepsForTest(gen, useEmbeddedWeights)
     tEnd = T
-    err = [gen.errorDahlquist(lam, u0, tEnd, nS) for nS in nSteps]
+    err = [gen.errorDahlquist(lam, u0, tEnd, nS, useEmbeddedWeights=useEmbeddedWeights) for nS in nSteps]
     order, rmse = numericalOrder(nSteps, err)
+    expectedOrder = gen.orderEmbedded if useEmbeddedWeights else gen.order
+
     assert rmse < 0.02, f"rmse to high ({rmse}) for {scheme} : {err}"
     if nNodes < 4:
         eps = 0.1
     else:
         eps = 0.25  # less constraining conditions for higher order
-    assert abs(order-gen.order) < eps, f"wrong numerical order ({order}) for {scheme} : {err}"
+    assert abs(order-expectedOrder) < eps, f"Expected order {expectedOrder:.2f}, but got {order:.2f} for {scheme} {err}"
+
+    if hasattr(gen, 'quadType'):
+        assert gen.rightIsNode == (gen.quadType in ['LOBATTO', 'RADAU-RIGHT'])