Martin burgers (#20)

* first test

* updates for playground

* burgers update

* further updates

* cleanups

* u

* updated test

* u

* added IMEX12 and IMEX21

* dahlquist arbitrary t0 integrator

* u

---------

Co-authored-by: SCHREIBER Martin <martin.schreiber@inria.fr>

Author: schreiberx

qmat/playgrounds/martin/diff_eqs/__init__.py

@@ -0,0 +1,10 @@
+
+from qmat.playgrounds.martin.diff_eqs.de_solver import DESolver
+from qmat.playgrounds.martin.diff_eqs.burgers import Burgers
+from qmat.playgrounds.martin.diff_eqs.dahlquist2 import Dahlquist2
+
+__all__ = [
+    "DESolver",
+    "Burgers",
+    "Dahlquist2",
+]

qmat/playgrounds/martin/diff_eqs/burgers.py

@@ -0,0 +1,159 @@
+import numpy as np
+from qmat.playgrounds.martin.diff_eqs.de_solver import DESolver
+
+
+class Burgers(DESolver):
+    """
+    Class to handle the 1D viscous Burgers' equation.
+    """
+
+    def __init__(self, N: int, nu: float, domain_size: float = 2.0 * np.pi):
+        # Resolution
+        self._N: int = N
+
+        # Viscosity
+        self._nu: float = nu
+
+        # Domain size
+        self._domain_size: float = domain_size
+
+        # Prepare spectral differentiation values
+        self._d_dx_ = 1j * np.fft.fftfreq(self._N, d=1.0 / self._N) * 2.0 * np.pi / self._domain_size
+
+    def _d_dx(self, u: np.ndarray) -> np.ndarray:
+        """Compute the spatial derivative of `u` using spectral methods.
+
+        Parameters
+        ----------
+        u : np.ndarray
+            Array of shape (N,) representing the solution at the current time step.
+
+        Returns
+        -------
+        du_dx : np.ndarray
+            Array of shape (N,) representing the spatial derivative of `u`.
+        """
+        u_hat = np.fft.fft(u)
+        du_dx_hat = u_hat * self._d_dx_
+        du_dx = np.fft.ifft(du_dx_hat).real
+        return du_dx
+
+    def initial_u0(self, mode: str) -> np.ndarray:
+        """Compute some initial conditions for the 1D viscous Burgers' equation."""
+
+        if mode == "sine":
+            x = np.linspace(0, self._domain_size, self._N, endpoint=False)
+            u0 = np.sin(x)
+
+        elif mode == "hat":
+            x = np.linspace(0, self._domain_size, self._N, endpoint=False)
+            u0 = np.where((x >= np.pi / 2) & (x <= 3 * np.pi / 2), 1.0, 0.0)
+
+        elif mode == "random":
+            np.random.seed(42)
+            u0 = np.random.rand(self._N)
+
+        else:
+            raise ValueError(f"Unknown initial condition mode: {mode}")
+
+        return u0
+
+    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.
+
+        Returns
+        -------
+        f : np.ndarray
+            Array of shape (N,) representing the right-hand side evaluated at `u`.
+        """
+        # Compute spatial derivatives using spectral methods
+        u_hat = np.fft.fft(u)
+        du_dx_hat = self._d_dx_ * u_hat
+        d2u_dx2_hat = (self._d_dx_**2) * u_hat
+
+        du_dx = np.fft.ifft(du_dx_hat).real
+        d2u_dx2 = np.fft.ifft(d2u_dx2_hat).real
+
+        f = -u * du_dx + self._nu * d2u_dx2
+        return f
+
+    def int_f(self, u0: np.ndarray, dt: float, t: float = 0.0) -> np.ndarray:
+        """
+        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
+
+        Parameters
+        ----------
+        u0 : np.ndarray
+            Array of shape (N,) representing the initial condition.
+        t : float
+            Time at which to evaluate the analytical solution.
+
+        Returns
+        -------
+        u_analytical : np.ndarray
+            Array of shape (N,) representing the analytical solution at time `t`.
+        """
+
+        if self._nu < 0.05:
+            print("Viscosity is very small which can lead to errors in analytical solution!")
+
+        u0_hat = np.fft.fft(u0)
+
+        # Divide by d/dx operator in spectral space
+        tmp = np.zeros_like(u0_hat, dtype=complex)
+        tmp[1:] = u0_hat[1:] / self._d_dx_[1:]
+
+        # Back to physical space
+        phi = np.fft.ifft(tmp).real
+
+        # Apply exp(...)
+        phi = np.exp(-phi / (2 * self._nu))
+
+        phi_hat = np.fft.fft(phi)
+
+        # Solve directly the heat equation in spectral space with exponential integration
+        phi_hat = phi_hat * np.exp(self._nu * self._d_dx_**2 * (t+dt))
+
+        phi = np.fft.ifft(phi_hat)
+        phi = np.log(phi)
+
+        phi_hat = np.fft.fft(phi)
+
+        u1_hat = -2.0 * self._nu * phi_hat * self._d_dx_
+        return np.fft.ifft(u1_hat).real
+
+    def test(self):
+        """
+        Run test for currently set up Burgers instance.
+        """
+        x = np.linspace(0, self._domain_size, self._N, endpoint=False)
+        w = 2.0 * np.pi / self._domain_size
+        u0 = np.sin(x * w)
+
+        u1_analytical = np.cos(x * w) * w
+        u1_num = self._d_dx(u0)
+
+        error: float = np.max(np.abs(u1_num - u1_analytical))
+
+        if error > 1e-10:
+            raise AssertionError(f"Test failed: error {error} too large for domain size {self._domain_size}.")
+
+    def run_tests(self):
+        """
+        Run basic tests to verify the correctness of the implementation.
+
+        This doesn't change the current instance, but will create new instances.
+        """
+
+        for domain_size in [2.0 * np.pi, 1.0, 9.0]:
+            burgers = Burgers(N=128, nu=0.01, domain_size=domain_size)
+            burgers.test()

qmat/playgrounds/martin/diff_eqs/dahlquist.py

@@ -0,0 +1,121 @@
+import numpy as np
+from qmat.playgrounds.martin.diff_eqs.de_solver import DESolver
+
+
+class Dahlquist(DESolver):
+    """
+    Standard Dahlquist test equation with two eigenvalues.
+    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, 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 picardF(self, u: np.ndarray, dt: float, t: float) -> np.ndarray:
+        """
+        Exactly integrate over one time step using the analytical solution
+        (=Picard).
+        """
+        lam = self.lam1 + self.lam2
+        s = self.ext_scalar
+        assert isinstance(t, float)
+        retval = np.exp(t * lam) * u + s * np.sin(t)
+
+        assert retval.shape == u.shape
+        return retval
+
+    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 evalF1(self, u: np.ndarray, t: float) -> np.ndarray:
+        lam = self.lam1
+        retval = lam * u
+
+        assert retval.shape == u.shape
+        return retval
+
+    def evalF2(self, u: np.ndarray, t: float) -> np.ndarray:
+        lam = 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 fSolve1(self, u: np.ndarray, dt: float, t: float) -> np.ndarray:
+        retval = u / (1.0 - dt * self.lam1)
+
+        assert retval.shape == u.shape
+        return retval
+
+    def fSolve2(self, u: np.ndarray, dt: float, t: float) -> np.ndarray:
+        t1 = t + dt
+        lam = 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 int_f_t0(self, u0: np.ndarray, dt: float) -> np.ndarray:
+        lam = self.lam1 + self.lam2
+        s = self.ext_scalar
+        assert isinstance(dt, float)
+        retval = np.exp(dt * lam) * u0 + s * np.sin(dt)
+
+        assert retval.shape == u0.shape
+        return retval
+
+    def int_f(self, u0: np.ndarray, dt: float, t: float = 0.0) -> np.ndarray:
+        """
+        Integrate the solution from t0 to t1.
+        """
+
+        assert isinstance(t, (float, int))
+        assert isinstance(dt, (float, int))
+
+        if t == 0:
+            return self.int_f_t0(u0, dt=dt)
+
+        # Lambda
+        lam = self.lam1 + self.lam2
+
+        s = self.ext_scalar
+
+        retval = np.exp(dt * lam) * (u0 - s*np.sin(t)) + s * np.sin(t+dt)
+
+        assert retval.shape == u0.shape
+        return retval

qmat/playgrounds/martin/diff_eqs/dahlquist2.py

@@ -0,0 +1,44 @@
+import numpy as np
+from qmat.playgrounds.martin.diff_eqs.de_solver import DESolver
+
+
+class Dahlquist2(DESolver):
+    """
+    Modified Dahlquist test equation with superposition of two
+    frequencies in the solution u(t).
+
+    u(t) = 0.5*(exp(lam1*t) + exp(lam2*t)) * u(0)
+    """
+
+    def __init__(self, lam1: complex, lam2: complex, s: float = 0.6):
+        """
+        :param lam1: First eigenvalue (complex)
+        :param lam2: Second eigenvalue (complex)
+        :param s: Weighting between the two exponentials in the solution
+        """
+        self.lam1: complex = lam1
+        self.lam2: complex = lam2
+        # Weighting between the two exponentials in the solution
+        # to avoid division by 0
+        self.s: float = s
+
+        assert 0 <= self.s <= 1, "s must be in [0,1]"
+
+    def initial_u0(self, mode: str = None) -> np.ndarray:
+        return np.array([1.0 + 0.0j], dtype=np.complex128)
+
+    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))
+            * u
+        )
+        assert retval.shape == u.shape
+        return retval
+
+    def int_f(self, u0: np.ndarray, dt: float, t: float = 0.0) -> np.ndarray:
+        assert isinstance(dt, float)
+        tfinal = t + dt
+        retval = (self.s * np.exp(tfinal * self.lam1) + (1.0 - self.s) * np.exp(tfinal * self.lam2)) * u0
+        assert retval.shape == u0.shape
+        return retval