Martin burgers

qmat/playgrounds/martin/diff_eqs/burgers.py

@@ -1,6 +1,5 @@
 import numpy as np
 from qmat.playgrounds.martin.diff_eqs.de_solver import DESolver
-from qmat.playgrounds.martin.time_integration.rk_integration import RKIntegration
 
 
 class Burgers(DESolver):
@@ -59,8 +58,9 @@ def initial_u0(self, mode: str) -> np.ndarray:
 
         return u0
 
-    def du_dt(self, u: np.ndarray, t: float) -> np.ndarray:
-        """Evaluate the right-hand side of the 1D viscous Burgers' equation.
+    def evalF(self, u: np.ndarray, t: float) -> np.ndarray:
+        """
+        Evaluate the right-hand side of the 1D viscous Burgers' equation.
 
         Parameters
         ----------
@@ -84,7 +84,8 @@ def du_dt(self, u: np.ndarray, t: float) -> np.ndarray:
         return f
 
     def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray:
-        """Compute the analytical solution of the 1D viscous Burgers' equation at time `t`.
+        """
+        Compute the analytical solution of the 1D viscous Burgers' equation at time `t`.
 
         See
         https://gitlab.inria.fr/sweet/sweet/-/blob/6f20b19f246bf6fcc7ace1b69567326d1da78635/src/programs/_pde_burgers/time/Burgers_Cart2D_TS_ln_cole_hopf.cpp

qmat/playgrounds/martin/diff_eqs/dahlquist.py

@@ -5,26 +5,51 @@
 class Dahlquist(DESolver):
     """
     Standard Dahlquist test equation with two eigenvalues.
-    d/dt u(t) = (lam1 + lam2) * u(t)
+    Optionally, some external frequency forcing can be added which is
+    configurable through `ext_scalar`.
+
+    u(t) = exp(t*(lam1+lam2))*u(0) + s*sin(t)
+
+    d/dt u(t) = (lam1 + lam2) * (u(t) - s*sin(t)) + s*cos(t)
     """
 
-    def __init__(self, lam1: complex, lam2: complex):
+    def __init__(self, lam1: complex, lam2: complex, ext_scalar: float = 0.0):
         # Lambda 1
         self.lam1: float = lam1
 
         # Lambda 2
         self.lam2: float = lam2
 
+        # Scaling of external sin(t) frequency. Set to 0 to deactivate.
+        self.ext_scalar = ext_scalar
+
     def initial_u0(self, mode: str = None) -> np.ndarray:
         return np.array([1.0 + 0.0j], dtype=np.complex128)
 
-    def du_dt(self, u: np.ndarray, t: float) -> np.ndarray:
-        retval = (self.lam1 + self.lam2) * u
+    def evalF(self, u: np.ndarray, t: float) -> np.ndarray:
+        lam = self.lam1 + self.lam2
+        s = self.ext_scalar
+        retval = lam * (u - s * np.sin(t)) + s * np.cos(t)
+
+        assert retval.shape == u.shape
+        return retval
+
+    def fSolve(self, u: np.ndarray, dt: float, t: float) -> np.ndarray:
+        t1 = t + dt
+        lam = self.lam1 + self.lam2
+        s = self.ext_scalar
+
+        rhs = u - s * dt * (lam * np.sin(t1) - np.cos(t1))
+        retval = rhs / (1.0 - dt * lam)
+
         assert retval.shape == u.shape
         return retval
 
     def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray:
+        lam = self.lam1 + self.lam2
+        s = self.ext_scalar
         assert isinstance(t, float)
-        retval = (np.exp(t * (self.lam1 + self.lam2))) * u0
+        retval = np.exp(t * lam) * u0 + s * np.sin(t)
+
         assert retval.shape == u0.shape
         return retval

qmat/playgrounds/martin/diff_eqs/dahlquist2.py

@@ -27,7 +27,7 @@ def __init__(self, lam1: complex, lam2: complex, s: float = 0.6):
     def initial_u0(self, mode: str = None) -> np.ndarray:
         return np.array([1.0 + 0.0j], dtype=np.complex128)
 
-    def du_dt(self, u: np.ndarray, t: float) -> np.ndarray:
+    def evalF(self, u: np.ndarray, t: float) -> np.ndarray:
         retval = (
             (self.lam1 * self.s * np.exp(t * self.lam1) + self.lam2 * (1.0 - self.s) * np.exp(t * self.lam2))
             / (self.s * np.exp(t * self.lam1) + (1.0 - self.s) * np.exp(t * self.lam2))

qmat/playgrounds/martin/diff_eqs/de_solver.py

@@ -5,16 +5,81 @@
 class DESolver(ABC):
 
     @abstractmethod
-    def du_dt(self, u: np.ndarray, t: float) -> np.ndarray:
+    def evalF(self, u: np.ndarray, t: float) -> np.ndarray:
         """Evaluate the right-hand side of the 1D viscous Burgers' equation.
 
         Parameters
         ----------
         u : np.ndarray
             Array of shape (N,) representing the solution at the current time step.
-
+        t : float
+            Current timestamp.
         Returns
         -------
         f : np.ndarray
             Array of shape (N,) representing the right-hand side evaluated at `u`.
         """
+        pass
+
+    # This is optional since not every DE might have a solver for backward Euler
+    # @abstractmethod
+    def fSolve(self, rhs: np.ndarray, dt: float, t: float) -> np.ndarray:
+        """Solve the right-hand side of an equation implicitly.
+
+        # Solving this equation implicitly...
+        u_t = f(u, t)
+
+        # ... means to u_new
+        u_new - dt * F(u_new, t + dt) = rhs
+
+        Parameters
+        ----------
+        rhs : np.ndarray
+            Right hand as given above.
+        t : float
+            Current timestamp.
+            Future one will be computed as `t + dt`
+        dt : float
+            Time step size.
+
+        Returns
+        -------
+        u_new : np.ndarray
+            Array of shape (N,) representing the solution at the next time step.
+        """
+        pass
+
+    @abstractmethod
+    def initial_u0(self, mode: str) -> np.ndarray:
+        """Compute some initial conditions for the 1D viscous Burgers' equation.
+
+        Parameters
+        ----------
+        mode : str
+            The type of initial condition to generate.
+
+        Returns
+        -------
+        u0 : np.ndarray
+            Array of shape (N,) representing the initial condition.
+        """
+        pass
+
+    @abstractmethod
+    def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray:
+        """
+        Compute the (analytical) solution at time `t`.
+
+        Parameters
+        ----------
+        u0 : np.ndarray
+            Array of shape (N,) representing the initial condition.
+        t : float
+            Time at which to evaluate the solution.
+
+        Returns
+        -------
+        u_analytical : np.ndarray
+            Solution at time `t`.
+        """
+        pass

qmat/playgrounds/martin/diff_eqs/two_freq.py

@@ -22,8 +22,12 @@ def __init__(self, lam1: complex, lam2: complex, lam3: complex = None, lam4: com
         self.lam4: complex = lam4 if lam4 is not None else lam2
 
         self.L = np.array([[self.lam1, self.lam2], [self.lam3, self.lam4]], dtype=np.complex128)
+        self.L1 = np.copy(self.L)
+        self.L1[:, 1] = 0
+        self.L2 = np.copy(self.L)
+        self.L2[:, 0] = 0
 
-        # compute eigenvalues and eigenvectors of L
+        # Compute eigenvalues and eigenvectors of L
         self.eigvals, self.eigvecs = np.linalg.eig(self.L)
 
         if not np.all(np.isclose(np.real(self.eigvals), 0)):
@@ -32,16 +36,48 @@ def __init__(self, lam1: complex, lam2: complex, lam3: complex = None, lam4: com
     def initial_u0(self, mode: str = None) -> np.ndarray:
         return np.array([1, 0.1], dtype=np.complex128)
 
-    def du_dt(self, u: np.ndarray, t: float) -> np.ndarray:
+    def evalF(self, u: np.ndarray, t: float) -> np.ndarray:
         assert isinstance(t, float)
         retval = self.L @ u
         assert retval.shape == u.shape
         return retval
 
+    def fSolve(self, rhs: np.ndarray, dt: float, t: float) -> np.ndarray:
+        """
+        Solve
+            `u_t = L u`
+        implicitly, i.e., return u_new such that
+            u_new = u + dt * L u_new
+        <=> u_new - dt * L u_new = u
+        <=> (I - dt * L) u_new = u
+        """
+        assert isinstance(t, float)
+        retval = np.linalg.solve(np.eye(rhs.shape[0]) - dt*self.L, rhs)
+        assert retval.shape == rhs.shape
+        return retval
+
     def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray:
         from scipy.linalg import expm
 
         assert isinstance(t, float)
         retval = expm(self.L * t) @ u0
         assert retval.shape == u0.shape
         return retval
+
+    def evalF1(self, u: np.ndarray, t: float) -> np.ndarray:
+        """
+        Solely compute tendencies of first frequency
+        """
+        assert isinstance(t, float)
+        retval = self.L1 @ u
+        assert retval.shape == u.shape
+        return retval
+
+    def evalF2(self, u: np.ndarray, t: float) -> np.ndarray:
+        """
+        Solely compute tendencies of first frequency
+        """
+        assert isinstance(t, float)
+        retval = self.L2 @ u
+        assert retval.shape == u.shape
+        return retval

qmat/playgrounds/martin/ex_two_freq_integrate.py

@@ -8,7 +8,9 @@
 T: float = 2 * 2 * np.pi  # Time interval
 t: float = 0.0  # Starting time
 
-time_integration = "sdc"
+time_integration: str = "sdc"
+sdc_micro_time_integration: str = "irk1"
+sdc_num_sweeps: int = 4
 
 two_freq: TwoFreq = TwoFreq(lam1=1.0j, lam2=20.0j, lam3=0.5j)
 u0 = two_freq.initial_u0()
@@ -42,7 +44,9 @@
             t_ = np.append(t_, t + (n + 1) * dt)
 
     elif time_integration == "sdc":
-        sdci = SDCIntegration(num_nodes=5, node_type="LEGENDRE", quad_type="LOBATTO", num_sweeps=4)
+        sdci = SDCIntegration(
+            num_nodes=5, node_type="LEGENDRE", quad_type="LOBATTO", num_sweeps=sdc_num_sweeps, micro_time_integration=sdc_micro_time_integration
+        )
 
         for n in range(num_timesteps):
             u = sdci.integrate(u, t + n * dt, dt, two_freq)

qmat/playgrounds/martin/tests/test_dahlquist.py

@@ -13,62 +13,74 @@ def test_dahlquist():
     dahlquist: Dahlquist = Dahlquist(lam1=1.0j, lam2=0.1j)
 
     for time_integration in ["rk1", "rk2", "rk4", "sdc"]:
-        print("="*80)
-        print(f"Time integration method: {time_integration}")
-        print("="*80)
-        results = []
+        if time_integration == "sdc":
+            micro_time_integration_ = ["erk1", "irk1"]
+        else:
+            micro_time_integration_ = ["-"]
 
-        u_analytical = dahlquist.u_solution(u0, t=T)
+        for micro_time_integration in micro_time_integration_:
+            print("=" * 80)
+            print(f"Time integration method: {time_integration} ({micro_time_integration})")
+            print("=" * 80)
+            results = []
 
-        for nt in range(4):
+            u_analytical = dahlquist.u_solution(u0, t=T)
 
-            num_timesteps = 2**nt * 10
-            print(f"Running simulation with num_timesteps={num_timesteps}")
+            for nt in range(4):
 
-            dt = T / num_timesteps
+                num_timesteps = 2**nt * 10
+                print(f"Running simulation with num_timesteps={num_timesteps}")
 
-            u0 = dahlquist.initial_u0()
+                dt = T / num_timesteps
 
-            u = u0.copy()
+                u0 = dahlquist.initial_u0()
 
-            if time_integration in RKIntegration.supported_methods:
-                rki = RKIntegration(method=time_integration)
+                u = u0.copy()
 
-                u = rki.integrate_n(u, t, dt, num_timesteps, dahlquist)
+                if time_integration in RKIntegration.supported_methods:
+                    rki = RKIntegration(method=time_integration)
 
-            elif time_integration == "sdc":
-                sdci = SDCIntegration(num_nodes=3, node_type="LEGENDRE", quad_type="LOBATTO", num_sweeps=1)
-                u = sdci.integrate_n(u, t, dt, num_timesteps, dahlquist)
+                    u = rki.integrate_n(u, t, dt, num_timesteps, dahlquist)
+
+                elif time_integration == "sdc":
+                    sdci = SDCIntegration(
+                        num_nodes=3,
+                        node_type="LEGENDRE",
+                        quad_type="LOBATTO",
+                        num_sweeps=3,
+                        micro_time_integration=micro_time_integration,
+                    )
+                    u = sdci.integrate_n(u, t, dt, num_timesteps, dahlquist)
 
-            else:
-                raise Exception("TODO")
+                else:
+                    raise Exception("TODO")
 
-            error = np.max(np.abs(u - u_analytical))
-            results.append({"N": num_timesteps, "dt": dt, "error": error})
+                error = np.max(np.abs(u - u_analytical))
+                results.append({"N": num_timesteps, "dt": dt, "error": error, "mti": micro_time_integration})
 
-        prev_error = None
-        for r in results:
-            if prev_error is None:
-                conv = None
-            else:
-                conv = np.log2(prev_error / r["error"])
+            prev_error = None
+            for r in results:
+                if prev_error is None:
+                    conv = None
+                else:
+                    conv = np.log2(prev_error / r["error"])
 
-            print(f" - N={r["N"]}, dt={r["dt"]:.6e}, error={r["error"]:.6e}, conv={conv}")
-            prev_error = r["error"]
-            r["conv"] = conv
+                print(f" - N={r["N"]}, dt={r["dt"]:.6e}, error={r["error"]:.6e}, conv={conv}, mti={micro_time_integration}")
+                prev_error = r["error"]
+                r["conv"] = conv
 
-        if time_integration == "rk1":
-            assert results[-1]["error"] < 1e-2
-            assert np.abs(results[-1]["conv"]-1.0) < 1e-2
+            if time_integration == "rk1":
+                assert results[-1]["error"] < 1e-2
+                assert np.abs(results[-1]["conv"] - 1.0) < 1e-2
 
-        elif time_integration == "rk2":
-            assert results[-1]["error"] < 1e-4
-            assert np.abs(results[-1]["conv"]-2.0) < 1e-3
+            elif time_integration == "rk2":
+                assert results[-1]["error"] < 1e-4
+                assert np.abs(results[-1]["conv"] - 2.0) < 1e-3
 
-        elif time_integration == "rk4":
-            assert results[-1]["error"] < 1e-9
-            assert np.abs(results[-1]["conv"]-4.0) < 1e-4
+            elif time_integration == "rk4":
+                assert results[-1]["error"] < 1e-9
+                assert np.abs(results[-1]["conv"] - 4.0) < 1e-4
 
-        elif time_integration == "sdc":
-            assert results[-1]["error"] < 1e-4
-            assert np.abs(results[-1]["conv"]-2.0) < 1e-3
+            elif time_integration == "sdc":
+                assert results[-1]["error"] < 1e-8
+                assert np.abs(results[-1]["conv"] - 4.0) < 1e-3