From b5a24284a967c7cf1199c7e91e63a48618ee7176 Mon Sep 17 00:00:00 2001 From: SCHREIBER Martin Date: Fri, 17 Oct 2025 15:09:56 +0200 Subject: [PATCH 01/11] first test --- examples/burgers.py | 191 +++++++++++++++++++ examples/example_1_burgers.py | 253 +++++++++++++++++++++++++ tests/test_qdelta/test_timestepping.py | 2 +- 3 files changed, 445 insertions(+), 1 deletion(-) create mode 100644 examples/burgers.py create mode 100644 examples/example_1_burgers.py diff --git a/examples/burgers.py b/examples/burgers.py new file mode 100644 index 0000000..677e42e --- /dev/null +++ b/examples/burgers.py @@ -0,0 +1,191 @@ +import numpy as np + + +class Burgers: + """ + 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 eval_f(self, u: np.ndarray) -> 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 analytical_integration(self, u0: np.ndarray, t: float) -> 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) + + 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 step_rk1(self, u: np.ndarray, dt: float) -> np.ndarray: + # Forward Euler + return u + dt * self.eval_f(u) + + def step_rk1_n(self, u: np.ndarray, dt: float, n: int) -> np.ndarray: + for _ in range(n): + u = self.step_rk1(u, dt) + return u + + def step_rk2(self, u: np.ndarray, dt: float) -> np.ndarray: + # Heun's method + k1 = self.eval_f(u) + k2 = self.eval_f(u + 0.5 * dt * k1) + return u + dt * k2 + + def step_rk2_n(self, u: np.ndarray, dt: float, n: int) -> np.ndarray: + for _ in range(n): + u = self.step_rk2(u, dt) + return u + + def step_rk4(self, u: np.ndarray, dt: float) -> np.ndarray: + # RK4 stages + k1 = self.eval_f(u) + k2 = self.eval_f(u + 0.5 * dt * k1) + k3 = self.eval_f(u + 0.5 * dt * k2) + k4 = self.eval_f(u + dt * k3) + + # Update solution + return u + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6 + + def step_rk4_n(self, u: np.ndarray, dt: float, n: int) -> np.ndarray: + for _ in range(n): + u = self.step_rk4(u, dt) + return u + + 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() diff --git a/examples/example_1_burgers.py b/examples/example_1_burgers.py new file mode 100644 index 0000000..7e71e56 --- /dev/null +++ b/examples/example_1_burgers.py @@ -0,0 +1,253 @@ +# +# Example adopted from `test_4_sdc.py` +# + + +from itertools import product +import numpy as np +from itertools import product as _product +from burgers import Burgers +import pandas as pd + + +def solveBurgersSDC( + burgers: Burgers, u0, T, nSteps: int, nSweeps: int, Q: np.ndarray, QDelta: np.ndarray, weights=None, monitors=None +) -> float: + r""" + Solve the Burgers' problem with SDC. + + Parameters + ---------- + u0 : complex or float + The initial solution :math:`u_0`. + T : float + Final time :math:`T`. + nSteps : int + Number of time-step for the whole :math:`[0,T]` interval. + nSweeps : int + Number of SDC sweeps. + Q : np.ndarray + Quadrature matrix :math:`Q` used for SDC. + QDelta : np.ndarray + Approximate quadrature matrix :math:`Q_\Delta` used for SDC. + If three dimensional, use the first dimension for the sweep index. + weights : np.ndarray, optional + Quadrature weights to use for the prologation. + If None, prolongation is not performed. The default is None. + + Returns + ------- + uNum : np.ndarray + Array containing the `nSteps+1` solutions :math:`\{u(0), ..., u(T)\}`. + """ + nodes = Q.sum(axis=1) + nNodes = Q.shape[0] + dt = T / nSteps + times = np.linspace(0, T, nSteps + 1) + + QDelta = np.asarray(QDelta) + if QDelta.ndim == 3: + assert QDelta.shape == (nSweeps, nNodes, nNodes), "inconsistent shape for QDelta" + else: + assert QDelta.shape == (nNodes, nNodes), "inconsistent shape for QDelta" + QDelta = np.repeat(QDelta[None, ...], nSweeps, axis=0) + + # Preconditioner built for each sweeps + P = np.eye(nNodes)[None, ...] - lam * dt * QDelta + + current_u = u0 + + # Setup monitors if any + if monitors: + tmp = {} + for key in monitors: + assert key in ["residuals", "errors"], f"unknown key '{key}' for monitors" + tmp[key] = np.zeros((nSweeps + 1, nSteps, nNodes), dtype=complex) + monitors = tmp + + for i in range(nSteps): + + uNodes = np.ones(nNodes) * current_u + + for k in range(nSweeps): + b = current_u + burgers.eval_f(current_u) * dt * nodes + + b = current_u + lam * dt * (Q - QDelta[k]) @ uNodes + uNodes = np.linalg.solve(P[k], b) + + # Monitoring + if monitors: + if "residuals" in monitors: + monitors["residuals"][k + 1, i] = current_u + lam * dt * Q @ uNodes - uNodes + if "errors" in monitors: + monitors["errors"][k + 1, i] = uNodes - u0 * np.exp(lam * (times[i] + dt * nodes)) + + if weights is None: + current_u = uNodes[-1] + else: + current_u = current_u + lam * dt * weights.dot(uNodes) + + if monitors: + return current_u, monitors + else: + return current_u + + +from matplotlib import pyplot as plt + +N = 128 +nu = 0.2 + +T: float = 4.0 +domain_size: float = 2.0 * np.pi +x = np.linspace(0, domain_size, N, endpoint=False) + +burgers: Burgers = Burgers(N=N, nu=nu, domain_size=domain_size) +u0 = burgers.initial_u0("sine") + +burgers.run_tests() + +if 0: + nt: int = 1000 + dt: float = T / nt + + if 0: + for t in [_ * dt for _ in range(nt)]: + ut = burgers.analytical_integration(u0, t=t) + print(f"t={t:.3f}, ut[0]={ut[0]:.6f}") + + plt.plot(x, ut) + else: + ut = burgers.analytical_integration(u0, t=T) + plt.plot(x, u0, label="u0") + plt.plot(x, ut, label="ut") + + plt.legend() + plt.show() + + +class SDCIntegration: + def __init__( + self, num_nodes: int = 3, node_type: str = "LOBATTO", quad_type: str = "LOBATTO" + ): + from qmat.qcoeff.collocation import Collocation + import qmat.qdelta.timestepping as module + + coll = Collocation(nNodes=num_nodes, nodeType=node_type, quadType=quad_type) + + self.gen: module.FE = module.FE(coll.nodes) + + self.nodes, self.weights, self.q = coll.genCoeffs(form="N2N") + + self.q_delta: np.array = self.gen.getQDelta() + self.d_tau: np.array = self.gen.dTau + self.deltas: np.array = self.gen.deltas + + # Number of nodes + self.N = len(self.nodes) + + +def sdc_integrate(u0: np.array, dt: float, num_timesteps, burgers: Burgers) -> np.array: + sdci = SDCIntegration(num_nodes=3, node_type="LEGENDRE", quad_type="LOBATTO") + + if not np.isclose(sdci.nodes[0], 0.0): + raise Exception("SDC nodes must include the left interval boundary.") + + if not np.isclose(sdci.nodes[-1], 1.0): + raise Exception("SDC nodes must include the right interval boundary.") + + assert sdci.N == sdci.deltas.shape[0] + + n_sweeps = 2 + shape = (sdci.N,) + u0.shape + u = np.zeros(shape=shape, dtype=float) + + u[0, :] = u0 + eval_f_k0 = np.empty_like(u) + eval_f_k1 = np.empty_like(u) + + for _ in range(num_timesteps): + + # Propagate initial condition to all nodes + for m in range(0, sdci.N): + if m > 0: + u[m] = u[m-1] + dt * sdci.deltas[m] * eval_f_k0[m-1] + eval_f_k0[m] = burgers.eval_f(u[m]) + + if 1: + # Iteratively sweep over SDC nodes + for _ in range(n_sweeps): + for m in range(0, sdci.N): + + if m > 0: + qeval = sdci.q[m] @ eval_f_k0 + u[m] = u[m-1] + dt * (sdci.deltas[m] * (eval_f_k1[m-1] - eval_f_k0[m-1]) + qeval) + + eval_f_k1[m] = burgers.eval_f(u[m]) + + # Copy tendency arrays + eval_f_k0[:] = eval_f_k1[:] + + if 0: + # If we're using Radau-right, we can just use the last value + u[0] = u[-1] + else: + # Compute new starting value with quadrature on tendencies + u[0] = u[0] + dt * sdci.weights @ eval_f_k0 + + return u[0] + + +time_integration = "sdc" + +if 1: + results = [] + + u_analytical = burgers.analytical_integration(u0, t=T) + + for nt in range(4): + + num_timesteps = 2**nt * 1000 + print(f"Running simulation with num_timesteps={num_timesteps}") + + dt = T / num_timesteps + + burgers: Burgers = Burgers(N=N, nu=nu, domain_size=domain_size) + u0 = burgers.initial_u0("sine") + + u = u0.copy() + + if time_integration == "rk1": + u = burgers.step_rk1_n(u, dt, num_timesteps) + + elif time_integration == "rk2": + u = burgers.step_rk2_n(u, dt, num_timesteps) + + elif time_integration == "rk4": + u = burgers.step_rk4_n(u, dt, num_timesteps) + + elif time_integration == "sdc": + u = sdc_integrate(u, dt, num_timesteps, burgers) + + else: + raise Exception("TODO") + + plt.plot(x, u, label=f"numerical {num_timesteps}", linestyle="dashed") + + error = np.max(np.abs(u - u_analytical)) + results.append({"N": num_timesteps, "dt": dt, "error": 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"] + + plt.plot(x, u_analytical, label="analytical", linestyle="dotted") + + plt.legend() + plt.show() diff --git a/tests/test_qdelta/test_timestepping.py b/tests/test_qdelta/test_timestepping.py index 76fb11a..6a491a6 100644 --- a/tests/test_qdelta/test_timestepping.py +++ b/tests/test_qdelta/test_timestepping.py @@ -18,7 +18,7 @@ def testBE(nNodes, nodeType, quadType): coll = Collocation(nNodes, nodeType, quadType) nodes = coll.nodes - gen = module.BE(nodes) + gen = module.BE(nodes) QDelta = gen.getQDelta() assert np.allclose(np.tril(QDelta), QDelta), \ From 77577e700381156528de64b62fbc2e0e350dcbde Mon Sep 17 00:00:00 2001 From: SCHREIBER Martin Date: Mon, 20 Oct 2025 17:54:34 +0200 Subject: [PATCH 02/11] updates for playground --- examples/example_1_burgers.py | 253 ------------------ qmat/playground/diff_eqs/__init__.py | 10 + .../playground/diff_eqs}/burgers.py | 3 +- qmat/playground/diff_eqs/dahlquist.py | 12 + qmat/playground/diff_eqs/de_solver.py | 20 ++ qmat/playground/ex_burgers_integrate.py | 80 ++++++ qmat/playground/ex_burgers_plot.py | 41 +++ qmat/playground/ex_oscillation.py | 22 ++ qmat/playground/tests/__init__.py | 0 qmat/playground/tests/test_sdc_burgers.py | 76 ++++++ .../time_integration/sdc_integration.py | 83 ++++++ 11 files changed, 346 insertions(+), 254 deletions(-) delete mode 100644 examples/example_1_burgers.py create mode 100644 qmat/playground/diff_eqs/__init__.py rename {examples => qmat/playground/diff_eqs}/burgers.py (98%) create mode 100644 qmat/playground/diff_eqs/dahlquist.py create mode 100644 qmat/playground/diff_eqs/de_solver.py create mode 100644 qmat/playground/ex_burgers_integrate.py create mode 100644 qmat/playground/ex_burgers_plot.py create mode 100644 qmat/playground/ex_oscillation.py create mode 100644 qmat/playground/tests/__init__.py create mode 100644 qmat/playground/tests/test_sdc_burgers.py create mode 100644 qmat/playground/time_integration/sdc_integration.py diff --git a/examples/example_1_burgers.py b/examples/example_1_burgers.py deleted file mode 100644 index 7e71e56..0000000 --- a/examples/example_1_burgers.py +++ /dev/null @@ -1,253 +0,0 @@ -# -# Example adopted from `test_4_sdc.py` -# - - -from itertools import product -import numpy as np -from itertools import product as _product -from burgers import Burgers -import pandas as pd - - -def solveBurgersSDC( - burgers: Burgers, u0, T, nSteps: int, nSweeps: int, Q: np.ndarray, QDelta: np.ndarray, weights=None, monitors=None -) -> float: - r""" - Solve the Burgers' problem with SDC. - - Parameters - ---------- - u0 : complex or float - The initial solution :math:`u_0`. - T : float - Final time :math:`T`. - nSteps : int - Number of time-step for the whole :math:`[0,T]` interval. - nSweeps : int - Number of SDC sweeps. - Q : np.ndarray - Quadrature matrix :math:`Q` used for SDC. - QDelta : np.ndarray - Approximate quadrature matrix :math:`Q_\Delta` used for SDC. - If three dimensional, use the first dimension for the sweep index. - weights : np.ndarray, optional - Quadrature weights to use for the prologation. - If None, prolongation is not performed. The default is None. - - Returns - ------- - uNum : np.ndarray - Array containing the `nSteps+1` solutions :math:`\{u(0), ..., u(T)\}`. - """ - nodes = Q.sum(axis=1) - nNodes = Q.shape[0] - dt = T / nSteps - times = np.linspace(0, T, nSteps + 1) - - QDelta = np.asarray(QDelta) - if QDelta.ndim == 3: - assert QDelta.shape == (nSweeps, nNodes, nNodes), "inconsistent shape for QDelta" - else: - assert QDelta.shape == (nNodes, nNodes), "inconsistent shape for QDelta" - QDelta = np.repeat(QDelta[None, ...], nSweeps, axis=0) - - # Preconditioner built for each sweeps - P = np.eye(nNodes)[None, ...] - lam * dt * QDelta - - current_u = u0 - - # Setup monitors if any - if monitors: - tmp = {} - for key in monitors: - assert key in ["residuals", "errors"], f"unknown key '{key}' for monitors" - tmp[key] = np.zeros((nSweeps + 1, nSteps, nNodes), dtype=complex) - monitors = tmp - - for i in range(nSteps): - - uNodes = np.ones(nNodes) * current_u - - for k in range(nSweeps): - b = current_u + burgers.eval_f(current_u) * dt * nodes - - b = current_u + lam * dt * (Q - QDelta[k]) @ uNodes - uNodes = np.linalg.solve(P[k], b) - - # Monitoring - if monitors: - if "residuals" in monitors: - monitors["residuals"][k + 1, i] = current_u + lam * dt * Q @ uNodes - uNodes - if "errors" in monitors: - monitors["errors"][k + 1, i] = uNodes - u0 * np.exp(lam * (times[i] + dt * nodes)) - - if weights is None: - current_u = uNodes[-1] - else: - current_u = current_u + lam * dt * weights.dot(uNodes) - - if monitors: - return current_u, monitors - else: - return current_u - - -from matplotlib import pyplot as plt - -N = 128 -nu = 0.2 - -T: float = 4.0 -domain_size: float = 2.0 * np.pi -x = np.linspace(0, domain_size, N, endpoint=False) - -burgers: Burgers = Burgers(N=N, nu=nu, domain_size=domain_size) -u0 = burgers.initial_u0("sine") - -burgers.run_tests() - -if 0: - nt: int = 1000 - dt: float = T / nt - - if 0: - for t in [_ * dt for _ in range(nt)]: - ut = burgers.analytical_integration(u0, t=t) - print(f"t={t:.3f}, ut[0]={ut[0]:.6f}") - - plt.plot(x, ut) - else: - ut = burgers.analytical_integration(u0, t=T) - plt.plot(x, u0, label="u0") - plt.plot(x, ut, label="ut") - - plt.legend() - plt.show() - - -class SDCIntegration: - def __init__( - self, num_nodes: int = 3, node_type: str = "LOBATTO", quad_type: str = "LOBATTO" - ): - from qmat.qcoeff.collocation import Collocation - import qmat.qdelta.timestepping as module - - coll = Collocation(nNodes=num_nodes, nodeType=node_type, quadType=quad_type) - - self.gen: module.FE = module.FE(coll.nodes) - - self.nodes, self.weights, self.q = coll.genCoeffs(form="N2N") - - self.q_delta: np.array = self.gen.getQDelta() - self.d_tau: np.array = self.gen.dTau - self.deltas: np.array = self.gen.deltas - - # Number of nodes - self.N = len(self.nodes) - - -def sdc_integrate(u0: np.array, dt: float, num_timesteps, burgers: Burgers) -> np.array: - sdci = SDCIntegration(num_nodes=3, node_type="LEGENDRE", quad_type="LOBATTO") - - if not np.isclose(sdci.nodes[0], 0.0): - raise Exception("SDC nodes must include the left interval boundary.") - - if not np.isclose(sdci.nodes[-1], 1.0): - raise Exception("SDC nodes must include the right interval boundary.") - - assert sdci.N == sdci.deltas.shape[0] - - n_sweeps = 2 - shape = (sdci.N,) + u0.shape - u = np.zeros(shape=shape, dtype=float) - - u[0, :] = u0 - eval_f_k0 = np.empty_like(u) - eval_f_k1 = np.empty_like(u) - - for _ in range(num_timesteps): - - # Propagate initial condition to all nodes - for m in range(0, sdci.N): - if m > 0: - u[m] = u[m-1] + dt * sdci.deltas[m] * eval_f_k0[m-1] - eval_f_k0[m] = burgers.eval_f(u[m]) - - if 1: - # Iteratively sweep over SDC nodes - for _ in range(n_sweeps): - for m in range(0, sdci.N): - - if m > 0: - qeval = sdci.q[m] @ eval_f_k0 - u[m] = u[m-1] + dt * (sdci.deltas[m] * (eval_f_k1[m-1] - eval_f_k0[m-1]) + qeval) - - eval_f_k1[m] = burgers.eval_f(u[m]) - - # Copy tendency arrays - eval_f_k0[:] = eval_f_k1[:] - - if 0: - # If we're using Radau-right, we can just use the last value - u[0] = u[-1] - else: - # Compute new starting value with quadrature on tendencies - u[0] = u[0] + dt * sdci.weights @ eval_f_k0 - - return u[0] - - -time_integration = "sdc" - -if 1: - results = [] - - u_analytical = burgers.analytical_integration(u0, t=T) - - for nt in range(4): - - num_timesteps = 2**nt * 1000 - print(f"Running simulation with num_timesteps={num_timesteps}") - - dt = T / num_timesteps - - burgers: Burgers = Burgers(N=N, nu=nu, domain_size=domain_size) - u0 = burgers.initial_u0("sine") - - u = u0.copy() - - if time_integration == "rk1": - u = burgers.step_rk1_n(u, dt, num_timesteps) - - elif time_integration == "rk2": - u = burgers.step_rk2_n(u, dt, num_timesteps) - - elif time_integration == "rk4": - u = burgers.step_rk4_n(u, dt, num_timesteps) - - elif time_integration == "sdc": - u = sdc_integrate(u, dt, num_timesteps, burgers) - - else: - raise Exception("TODO") - - plt.plot(x, u, label=f"numerical {num_timesteps}", linestyle="dashed") - - error = np.max(np.abs(u - u_analytical)) - results.append({"N": num_timesteps, "dt": dt, "error": 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"] - - plt.plot(x, u_analytical, label="analytical", linestyle="dotted") - - plt.legend() - plt.show() diff --git a/qmat/playground/diff_eqs/__init__.py b/qmat/playground/diff_eqs/__init__.py new file mode 100644 index 0000000..c60dc51 --- /dev/null +++ b/qmat/playground/diff_eqs/__init__.py @@ -0,0 +1,10 @@ + +from qmat.playground.diff_eqs.de_solver import DESolver +from qmat.playground.diff_eqs.burgers import Burgers +from qmat.playground.diff_eqs.dahlquist import Dahlquist + +__all__ = [ + "DESolver", + "Burgers", + "Dahlquist", +] diff --git a/examples/burgers.py b/qmat/playground/diff_eqs/burgers.py similarity index 98% rename from examples/burgers.py rename to qmat/playground/diff_eqs/burgers.py index 677e42e..c528d96 100644 --- a/examples/burgers.py +++ b/qmat/playground/diff_eqs/burgers.py @@ -1,7 +1,8 @@ import numpy as np +from qmat.playground.diff_eqs.de_solver import DESolver -class Burgers: +class Burgers(DESolver): """ Class to handle the 1D viscous Burgers' equation. """ diff --git a/qmat/playground/diff_eqs/dahlquist.py b/qmat/playground/diff_eqs/dahlquist.py new file mode 100644 index 0000000..7647d1b --- /dev/null +++ b/qmat/playground/diff_eqs/dahlquist.py @@ -0,0 +1,12 @@ +import numpy as np +from qmat.playground.diff_eqs.de_solver import DESolver + + +class Dahlquist(DESolver): + + def __init__(self, lam1: float, lam2: float): + # Lambda 1 + self.lam1: float = lam1 + + # Lambda 2 + self.lam2: float = lam2 diff --git a/qmat/playground/diff_eqs/de_solver.py b/qmat/playground/diff_eqs/de_solver.py new file mode 100644 index 0000000..ecd1edb --- /dev/null +++ b/qmat/playground/diff_eqs/de_solver.py @@ -0,0 +1,20 @@ +from abc import ABC, abstractmethod +import numpy as np + + +class DESolver(ABC): + + @abstractmethod + def eval_f(self, u: np.ndarray) -> 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`. + """ diff --git a/qmat/playground/ex_burgers_integrate.py b/qmat/playground/ex_burgers_integrate.py new file mode 100644 index 0000000..8b581c5 --- /dev/null +++ b/qmat/playground/ex_burgers_integrate.py @@ -0,0 +1,80 @@ +# +# Example adopted from `test_4_sdc.py` +# + + +import numpy as np +from diff_eqs.burgers import Burgers +from time_integration.sdc_integration import SDCIntegration + + +from matplotlib import pyplot as plt + +N = 128 +nu = 0.2 + +T: float = 4.0 +domain_size: float = 2.0 * np.pi +t = np.linspace(0, domain_size, N, endpoint=False) + +burgers: Burgers = Burgers(N=N, nu=nu, domain_size=domain_size) +u0 = burgers.initial_u0("sine") + +burgers.run_tests() + + +time_integration = "sdc" + +if 1: + results = [] + + u_analytical = burgers.analytical_integration(u0, t=T) + + for nt in range(4): + + num_timesteps = 2**nt * 1000 + print(f"Running simulation with num_timesteps={num_timesteps}") + + dt = T / num_timesteps + + burgers: Burgers = Burgers(N=N, nu=nu, domain_size=domain_size) + u0 = burgers.initial_u0("sine") + + u = u0.copy() + + if time_integration == "rk1": + u = burgers.step_rk1_n(u, dt, num_timesteps) + + elif time_integration == "rk2": + u = burgers.step_rk2_n(u, dt, num_timesteps) + + elif time_integration == "rk4": + u = burgers.step_rk4_n(u, dt, num_timesteps) + + elif time_integration == "sdc": + sdci = SDCIntegration(num_nodes=3, node_type="LEGENDRE", quad_type="LOBATTO") + u = sdci.integrate_n(u, dt, num_timesteps, burgers) + + else: + raise Exception("TODO") + + plt.plot(t, u, label=f"numerical {num_timesteps}", linestyle="dashed") + + error = np.max(np.abs(u - u_analytical)) + results.append({"N": num_timesteps, "dt": dt, "error": 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"] + + if 1: + plt.plot(t, u_analytical, label="analytical", linestyle="dotted") + + plt.legend() + plt.show() diff --git a/qmat/playground/ex_burgers_plot.py b/qmat/playground/ex_burgers_plot.py new file mode 100644 index 0000000..9cd49fe --- /dev/null +++ b/qmat/playground/ex_burgers_plot.py @@ -0,0 +1,41 @@ +# +# Example adopted from `test_4_sdc.py` +# + + +import numpy as np +from diff_eqs.burgers import Burgers + + +from matplotlib import pyplot as plt + +N = 128 +nu = 0.2 + +T: float = 4.0 +domain_size: float = 2.0 * np.pi +t = np.linspace(0, domain_size, N, endpoint=False) + +burgers: Burgers = Burgers(N=N, nu=nu, domain_size=domain_size) +u0 = burgers.initial_u0("sine") + +burgers.run_tests() + +if 1: + nt: int = 1000 + dt: float = T / nt + + if 0: + for t in [_ * dt for _ in range(nt)]: + ut = burgers.analytical_integration(u0, t=t) + print(f"t={t:.3f}, ut[0]={ut[0]:.6f}") + + plt.plot(t, ut) + else: + ut = burgers.analytical_integration(u0, t=T) + plt.plot(t, u0, label="u0") + plt.plot(t, ut, label="ut") + + plt.legend() + plt.show() + diff --git a/qmat/playground/ex_oscillation.py b/qmat/playground/ex_oscillation.py new file mode 100644 index 0000000..c2fd05b --- /dev/null +++ b/qmat/playground/ex_oscillation.py @@ -0,0 +1,22 @@ +import numpy as np + +if 1: + N = 500 + u0 = 1.0 + t: np.array = np.linspace(0, 4*np.pi, N, endpoint=False) + lam1 = 10.j + lam2 = 1.j + + def u(u0, t) -> np.array: + return np.exp(t*lam1) * u0 + np.exp(t*lam2) * u0 + + def du_dt(u0, t) -> np.array: + return (lam1*np.exp(t*lam1) * u0 + lam2*np.exp(t*lam2)) * u0 + + u_eval: np.array = u(u0, t) + + from matplotlib import pyplot as plt + plt.plot(t, np.real(u_eval), label="Re(u)") + plt.plot(t, np.imag(u_eval), label="Im(u)") + plt.legend() + plt.show() diff --git a/qmat/playground/tests/__init__.py b/qmat/playground/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/qmat/playground/tests/test_sdc_burgers.py b/qmat/playground/tests/test_sdc_burgers.py new file mode 100644 index 0000000..97efd5c --- /dev/null +++ b/qmat/playground/tests/test_sdc_burgers.py @@ -0,0 +1,76 @@ + +import numpy as np +from qmat.playground.diff_eqs.burgers import Burgers +from qmat.playground.time_integration.sdc_integration import SDCIntegration + + +def test_sdc_burgers(): + N = 128 + nu = 0.2 + + T: float = 4.0 + domain_size: float = 2.0 * np.pi + t = np.linspace(0, domain_size, N, endpoint=False) + + burgers: Burgers = Burgers(N=N, nu=nu, domain_size=domain_size) + u0 = burgers.initial_u0("sine") + + burgers.run_tests() + + for time_integration in ["rk1", "rk2", "rk4", "sdc"]: + print("="*80) + print(f"Time integration method: {time_integration}") + print("="*80) + results = [] + + u_analytical = burgers.analytical_integration(u0, t=T) + + for nt in range(4): + + num_timesteps = 2**nt * 1000 + print(f"Running simulation with num_timesteps={num_timesteps}") + + dt = T / num_timesteps + + burgers: Burgers = Burgers(N=N, nu=nu, domain_size=domain_size) + u0 = burgers.initial_u0("sine") + + u = u0.copy() + + if time_integration == "rk1": + u = burgers.step_rk1_n(u, dt, num_timesteps) + + elif time_integration == "rk2": + u = burgers.step_rk2_n(u, dt, num_timesteps) + + elif time_integration == "rk4": + u = burgers.step_rk4_n(u, dt, num_timesteps) + + elif time_integration == "sdc": + sdci = SDCIntegration(num_nodes=3, node_type="LEGENDRE", quad_type="LOBATTO") + u = sdci.integrate_n(u, dt, num_timesteps, burgers) + + else: + raise Exception("TODO") + + error = np.max(np.abs(u - u_analytical)) + results.append({"N": num_timesteps, "dt": dt, "error": 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"] + + if time_integration == "rk1": + assert results[-1]["error"] < 1e-4 + elif time_integration == "rk2": + assert results[-1]["error"] < 1e-8 + elif time_integration == "rk4": + assert results[-1]["error"] < 1e-14 + elif time_integration == "sdc": + assert results[-1]["error"] < 1e-14 diff --git a/qmat/playground/time_integration/sdc_integration.py b/qmat/playground/time_integration/sdc_integration.py new file mode 100644 index 0000000..a50f7d1 --- /dev/null +++ b/qmat/playground/time_integration/sdc_integration.py @@ -0,0 +1,83 @@ +import numpy as np +from qmat.playground.diff_eqs.de_solver import DESolver + + +class SDCIntegration: + def __init__( + self, num_nodes: int = 3, node_type: str = "LOBATTO", quad_type: str = "LOBATTO" + ): + from qmat.qcoeff.collocation import Collocation + import qmat.qdelta.timestepping as module + + coll = Collocation(nNodes=num_nodes, nodeType=node_type, quadType=quad_type) + + self.gen: module.FE = module.FE(coll.nodes) + + self.nodes, self.weights, self.q = coll.genCoeffs(form="N2N") + + self.q_delta: np.array = self.gen.getQDelta() + self.d_tau: np.array = self.gen.dTau + self.deltas: np.array = self.gen.deltas + + # Number of nodes + self.N = len(self.nodes) + + def integrate(self, u0: np.array, dt: float, de_solver: DESolver) -> np.array: + if not np.isclose(self.nodes[0], 0.0): + raise Exception("SDC nodes must include the left interval boundary.") + + if not np.isclose(self.nodes[-1], 1.0): + raise Exception("SDC nodes must include the right interval boundary.") + + assert self.N == self.deltas.shape[0] + + n_sweeps = 2 + shape = (self.N,) + u0.shape + u = np.zeros(shape=shape, dtype=float) + + u[0, :] = u0 + eval_f_k0 = np.empty_like(u) + eval_f_k1 = np.empty_like(u) + + # Propagate initial condition to all nodes + for m in range(0, self.N): + if m > 0: + u[m] = u[m-1] + dt * self.deltas[m] * eval_f_k0[m-1] + eval_f_k0[m] = de_solver.eval_f(u[m]) + + if 1: + # Iteratively sweep over SDC nodes + for _ in range(n_sweeps): + for m in range(0, self.N): + + if m > 0: + qeval = self.q[m] @ eval_f_k0 + u[m] = u[m-1] + dt * (self.deltas[m] * (eval_f_k1[m-1] - eval_f_k0[m-1]) + qeval) + + eval_f_k1[m] = de_solver.eval_f(u[m]) + + # Copy tendency arrays + eval_f_k0[:] = eval_f_k1[:] + + if 0: + # If we're using Radau-right, we can just use the last value + u[0] = u[-1] + else: + # Compute new starting value with quadrature on tendencies + u[0] = u[0] + dt * self.weights @ eval_f_k0 + + return u[0] + + def integrate_n(self, u0: np.array, dt: float, num_timesteps, de_solver: DESolver) -> np.array: + if not np.isclose(self.nodes[0], 0.0): + raise Exception("SDC nodes must include the left interval boundary.") + + if not np.isclose(self.nodes[-1], 1.0): + raise Exception("SDC nodes must include the right interval boundary.") + + u_value = u0 + + for _ in range(num_timesteps): + u_value = self.integrate(u_value, dt, de_solver) + + return u_value From d29562db9af23f2e21902a435f30422f5bbb78cd Mon Sep 17 00:00:00 2001 From: SCHREIBER Martin Date: Tue, 21 Oct 2025 15:08:37 +0200 Subject: [PATCH 03/11] burgers update --- qmat/playground/diff_eqs/__init__.py | 4 +- qmat/playground/diff_eqs/burgers.py | 40 +--------- qmat/playground/diff_eqs/dahlquist.py | 20 ++++- qmat/playground/diff_eqs/dahlquist2.py | 36 +++++++++ qmat/playground/diff_eqs/de_solver.py | 2 +- qmat/playground/ex_burgers_integrate.py | 31 ++++---- qmat/playground/ex_burgers_plot.py | 3 +- qmat/playground/ex_dahlquist2_integrate.py | 78 +++++++++++++++++++ qmat/playground/ex_dahlquist2_plot.py | 19 +++++ qmat/playground/ex_oscillation.py | 22 ------ .../{test_sdc_burgers.py => test_burgers.py} | 30 +++---- qmat/playground/tests/test_dahlquist.py | 75 ++++++++++++++++++ qmat/playground/tests/test_dahlquist2.py | 76 ++++++++++++++++++ .../time_integration/rk_integration.py | 49 ++++++++++++ .../time_integration/sdc_integration.py | 51 ++++++------ 15 files changed, 418 insertions(+), 118 deletions(-) create mode 100644 qmat/playground/diff_eqs/dahlquist2.py create mode 100644 qmat/playground/ex_dahlquist2_integrate.py create mode 100644 qmat/playground/ex_dahlquist2_plot.py delete mode 100644 qmat/playground/ex_oscillation.py rename qmat/playground/tests/{test_sdc_burgers.py => test_burgers.py} (72%) create mode 100644 qmat/playground/tests/test_dahlquist.py create mode 100644 qmat/playground/tests/test_dahlquist2.py create mode 100644 qmat/playground/time_integration/rk_integration.py diff --git a/qmat/playground/diff_eqs/__init__.py b/qmat/playground/diff_eqs/__init__.py index c60dc51..77a6843 100644 --- a/qmat/playground/diff_eqs/__init__.py +++ b/qmat/playground/diff_eqs/__init__.py @@ -1,10 +1,10 @@ from qmat.playground.diff_eqs.de_solver import DESolver from qmat.playground.diff_eqs.burgers import Burgers -from qmat.playground.diff_eqs.dahlquist import Dahlquist +from qmat.playground.diff_eqs.dahlquist2 import Dahlquist2 __all__ = [ "DESolver", "Burgers", - "Dahlquist", + "Dahlquist2", ] diff --git a/qmat/playground/diff_eqs/burgers.py b/qmat/playground/diff_eqs/burgers.py index c528d96..118c602 100644 --- a/qmat/playground/diff_eqs/burgers.py +++ b/qmat/playground/diff_eqs/burgers.py @@ -1,5 +1,6 @@ import numpy as np from qmat.playground.diff_eqs.de_solver import DESolver +from qmat.playground.time_integration.rk_integration import RKIntegration class Burgers(DESolver): @@ -58,7 +59,7 @@ def initial_u0(self, mode: str) -> np.ndarray: return u0 - def eval_f(self, u: np.ndarray) -> np.ndarray: + def du_dt(self, u: np.ndarray, t: float) -> np.ndarray: """Evaluate the right-hand side of the 1D viscous Burgers' equation. Parameters @@ -82,7 +83,7 @@ def eval_f(self, u: np.ndarray) -> np.ndarray: f = -u * du_dx + self._nu * d2u_dx2 return f - def analytical_integration(self, u0: np.ndarray, t: float) -> np.ndarray: + def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray: """Compute the analytical solution of the 1D viscous Burgers' equation at time `t`. See @@ -129,41 +130,6 @@ def analytical_integration(self, u0: np.ndarray, t: float) -> np.ndarray: u1_hat = -2.0 * self._nu * phi_hat * self._d_dx_ return np.fft.ifft(u1_hat).real - def step_rk1(self, u: np.ndarray, dt: float) -> np.ndarray: - # Forward Euler - return u + dt * self.eval_f(u) - - def step_rk1_n(self, u: np.ndarray, dt: float, n: int) -> np.ndarray: - for _ in range(n): - u = self.step_rk1(u, dt) - return u - - def step_rk2(self, u: np.ndarray, dt: float) -> np.ndarray: - # Heun's method - k1 = self.eval_f(u) - k2 = self.eval_f(u + 0.5 * dt * k1) - return u + dt * k2 - - def step_rk2_n(self, u: np.ndarray, dt: float, n: int) -> np.ndarray: - for _ in range(n): - u = self.step_rk2(u, dt) - return u - - def step_rk4(self, u: np.ndarray, dt: float) -> np.ndarray: - # RK4 stages - k1 = self.eval_f(u) - k2 = self.eval_f(u + 0.5 * dt * k1) - k3 = self.eval_f(u + 0.5 * dt * k2) - k4 = self.eval_f(u + dt * k3) - - # Update solution - return u + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6 - - def step_rk4_n(self, u: np.ndarray, dt: float, n: int) -> np.ndarray: - for _ in range(n): - u = self.step_rk4(u, dt) - return u - def test(self): """ Run test for currently set up Burgers instance. diff --git a/qmat/playground/diff_eqs/dahlquist.py b/qmat/playground/diff_eqs/dahlquist.py index 7647d1b..321726d 100644 --- a/qmat/playground/diff_eqs/dahlquist.py +++ b/qmat/playground/diff_eqs/dahlquist.py @@ -3,10 +3,28 @@ class Dahlquist(DESolver): + """ + Standard Dahlquist test equation with two eigenvalues. + d/dt u(t) = (lam1 + lam2) * u(t) + """ - def __init__(self, lam1: float, lam2: float): + def __init__(self, lam1: complex, lam2: complex): # Lambda 1 self.lam1: float = lam1 # Lambda 2 self.lam2: float = lam2 + + def initial_u0(self, mode: str = None) -> np.ndarray: + return np.array([1.0 + 0.0j]) + + def du_dt(self, u: np.ndarray, t: float) -> np.ndarray: + retval = (self.lam1 + self.lam2) * u + assert retval.shape == u.shape + return retval + + def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray: + assert isinstance(t, float) + retval = (np.exp(t * (self.lam1 + self.lam2))) * u0 + assert retval.shape == u0.shape + return retval diff --git a/qmat/playground/diff_eqs/dahlquist2.py b/qmat/playground/diff_eqs/dahlquist2.py new file mode 100644 index 0000000..c6a10d1 --- /dev/null +++ b/qmat/playground/diff_eqs/dahlquist2.py @@ -0,0 +1,36 @@ +import numpy as np +from qmat.playground.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): + # Lambda 1 + self.lam1: float = lam1 + + # Lambda 2 + self.lam2: float = lam2 + + def initial_u0(self, mode: str = None) -> np.ndarray: + return np.array([1.0 + 0.0j]) + + def du_dt(self, u: np.ndarray, t: float) -> np.ndarray: + retval = ( + (self.lam1 * np.exp(t * self.lam1) + self.lam2 * np.exp(t * self.lam2)) + / (np.exp(t * self.lam1) + np.exp(t * self.lam2)) + * u + ) + assert retval.shape == u.shape + return retval + + def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray: + assert isinstance(t, float) + retval = 0.5*(np.exp(t * self.lam1) + np.exp(t * self.lam2)) * u0 + assert retval.shape == u0.shape + return retval diff --git a/qmat/playground/diff_eqs/de_solver.py b/qmat/playground/diff_eqs/de_solver.py index ecd1edb..b057554 100644 --- a/qmat/playground/diff_eqs/de_solver.py +++ b/qmat/playground/diff_eqs/de_solver.py @@ -5,7 +5,7 @@ class DESolver(ABC): @abstractmethod - def eval_f(self, u: np.ndarray) -> np.ndarray: + def du_dt(self, u: np.ndarray, t: float) -> np.ndarray: """Evaluate the right-hand side of the 1D viscous Burgers' equation. Parameters diff --git a/qmat/playground/ex_burgers_integrate.py b/qmat/playground/ex_burgers_integrate.py index 8b581c5..28caae4 100644 --- a/qmat/playground/ex_burgers_integrate.py +++ b/qmat/playground/ex_burgers_integrate.py @@ -4,18 +4,20 @@ import numpy as np -from diff_eqs.burgers import Burgers -from time_integration.sdc_integration import SDCIntegration +from qmat.playground.diff_eqs.burgers import Burgers +from qmat.playground.time_integration.sdc_integration import SDCIntegration +from qmat.playground.time_integration.rk_integration import RKIntegration from matplotlib import pyplot as plt N = 128 nu = 0.2 - T: float = 4.0 +t: float = 0.0 # Starting point in time + domain_size: float = 2.0 * np.pi -t = np.linspace(0, domain_size, N, endpoint=False) +t_ = np.linspace(0, domain_size, N, endpoint=False) burgers: Burgers = Burgers(N=N, nu=nu, domain_size=domain_size) u0 = burgers.initial_u0("sine") @@ -23,12 +25,12 @@ burgers.run_tests() -time_integration = "sdc" +time_integration = "rk2" if 1: results = [] - u_analytical = burgers.analytical_integration(u0, t=T) + u_analytical = burgers.u_solution(u0, t=T) for nt in range(4): @@ -42,23 +44,18 @@ u = u0.copy() - if time_integration == "rk1": - u = burgers.step_rk1_n(u, dt, num_timesteps) - - elif time_integration == "rk2": - u = burgers.step_rk2_n(u, dt, num_timesteps) - - elif time_integration == "rk4": - u = burgers.step_rk4_n(u, dt, num_timesteps) + if time_integration in RKIntegration.supported_methods: + rki = RKIntegration(method=time_integration) + u = rki.integrate_n(u, t, dt, num_timesteps, burgers) elif time_integration == "sdc": sdci = SDCIntegration(num_nodes=3, node_type="LEGENDRE", quad_type="LOBATTO") - u = sdci.integrate_n(u, dt, num_timesteps, burgers) + u = sdci.integrate_n(u, t, dt, num_timesteps, burgers) else: raise Exception("TODO") - plt.plot(t, u, label=f"numerical {num_timesteps}", linestyle="dashed") + plt.plot(t_, u, label=f"numerical {num_timesteps}", linestyle="dashed") error = np.max(np.abs(u - u_analytical)) results.append({"N": num_timesteps, "dt": dt, "error": error}) @@ -74,7 +71,7 @@ prev_error = r["error"] if 1: - plt.plot(t, u_analytical, label="analytical", linestyle="dotted") + plt.plot(t_, u_analytical, label="analytical", linestyle="dotted") plt.legend() plt.show() diff --git a/qmat/playground/ex_burgers_plot.py b/qmat/playground/ex_burgers_plot.py index 9cd49fe..2752209 100644 --- a/qmat/playground/ex_burgers_plot.py +++ b/qmat/playground/ex_burgers_plot.py @@ -32,10 +32,9 @@ plt.plot(t, ut) else: - ut = burgers.analytical_integration(u0, t=T) + ut = burgers.u_solution(u0, t=T) plt.plot(t, u0, label="u0") plt.plot(t, ut, label="ut") plt.legend() plt.show() - diff --git a/qmat/playground/ex_dahlquist2_integrate.py b/qmat/playground/ex_dahlquist2_integrate.py new file mode 100644 index 0000000..8a8683b --- /dev/null +++ b/qmat/playground/ex_dahlquist2_integrate.py @@ -0,0 +1,78 @@ +import numpy as np +from qmat.playground.diff_eqs.dahlquist2 import Dahlquist2 +from qmat.playground.diff_eqs.dahlquist import Dahlquist +from time_integration.sdc_integration import SDCIntegration +from qmat.playground.time_integration.rk_integration import RKIntegration +from matplotlib import pyplot as plt + + +u0 = np.array([1.0]) # Initial condition +T: float = 4 * np.pi # Time interval +T: float = 0.5 # Time interval +t: float = 0.0 # Starting time + +time_integration = "rk4" + + +results = [] + +for nt in range(4): + + num_timesteps = 2**nt * 10 + print(f"Running simulation with num_timesteps={num_timesteps}") + + dt = T / num_timesteps + + dahlquist2: Dahlquist2 = Dahlquist2(lam1=1.0j, lam2=1.0j) + u_analytical_fin = dahlquist2.u_solution(u0, t=T) + u0 = dahlquist2.initial_u0() + + u = u0.copy() + + u_: np.ndarray = np.array([u]) + t_: np.ndarray = np.array([t]) + + if time_integration in RKIntegration.supported_methods: + rki = RKIntegration(method=time_integration) + + for n in range(num_timesteps): + u = rki.integrate(u, t + n * dt, dt, dahlquist2) + + u_ = np.append(u_, u) + t_ = np.append(t_, t + (n + 1) * dt) + + elif time_integration == "sdc": + sdci = SDCIntegration(num_nodes=3, node_type="LEGENDRE", quad_type="LOBATTO") + + for n in range(num_timesteps): + u = sdci.integrate(u, t + n * dt, dt, dahlquist2) + + u_ = np.append(u_, u) + t_ = np.append(t_, t + (n + 1) * dt) + + else: + raise Exception("TODO") + + plt.plot(t_, np.real(u_), label=f"ndt={num_timesteps}, real", linestyle="dashed") + plt.plot(t_, np.imag(u_), label=f"ndt={num_timesteps}, imag", linestyle="solid") + + error = np.max(np.abs(u - u_analytical_fin)) + results.append({"N": num_timesteps, "dt": dt, "error": 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"] + +if 1: + u_analytical_fin = np.array([dahlquist2.u_solution(u0, t) for t in t_]) + plt.plot(t_, np.real(u_analytical_fin), label="analytical, real", linestyle="dotted") + plt.plot(t_, np.imag(u_analytical_fin), label="analytical, imag", linestyle="dotted") + + plt.legend() + plt.show() diff --git a/qmat/playground/ex_dahlquist2_plot.py b/qmat/playground/ex_dahlquist2_plot.py new file mode 100644 index 0000000..9376ee8 --- /dev/null +++ b/qmat/playground/ex_dahlquist2_plot.py @@ -0,0 +1,19 @@ +import numpy as np +from qmat.playground.diff_eqs.dahlquist2 import Dahlquist2 + + +N = 500 +u0 = np.array([1.0]) +t: np.array = np.linspace(0, 4*np.pi, N, endpoint=False) + +dahlquist = Dahlquist2(lam1=10.j, lam2=1.j) + +u_eval = np.array([dahlquist.u_solution(u0, _) for _ in t]) + + +if 1: + from matplotlib import pyplot as plt + plt.plot(t, np.real(u_eval), label="Re(u)") + plt.plot(t, np.imag(u_eval), label="Im(u)") + plt.legend() + plt.show() diff --git a/qmat/playground/ex_oscillation.py b/qmat/playground/ex_oscillation.py deleted file mode 100644 index c2fd05b..0000000 --- a/qmat/playground/ex_oscillation.py +++ /dev/null @@ -1,22 +0,0 @@ -import numpy as np - -if 1: - N = 500 - u0 = 1.0 - t: np.array = np.linspace(0, 4*np.pi, N, endpoint=False) - lam1 = 10.j - lam2 = 1.j - - def u(u0, t) -> np.array: - return np.exp(t*lam1) * u0 + np.exp(t*lam2) * u0 - - def du_dt(u0, t) -> np.array: - return (lam1*np.exp(t*lam1) * u0 + lam2*np.exp(t*lam2)) * u0 - - u_eval: np.array = u(u0, t) - - from matplotlib import pyplot as plt - plt.plot(t, np.real(u_eval), label="Re(u)") - plt.plot(t, np.imag(u_eval), label="Im(u)") - plt.legend() - plt.show() diff --git a/qmat/playground/tests/test_sdc_burgers.py b/qmat/playground/tests/test_burgers.py similarity index 72% rename from qmat/playground/tests/test_sdc_burgers.py rename to qmat/playground/tests/test_burgers.py index 97efd5c..4f4ec1a 100644 --- a/qmat/playground/tests/test_sdc_burgers.py +++ b/qmat/playground/tests/test_burgers.py @@ -2,15 +2,16 @@ import numpy as np from qmat.playground.diff_eqs.burgers import Burgers from qmat.playground.time_integration.sdc_integration import SDCIntegration +from qmat.playground.time_integration.rk_integration import RKIntegration -def test_sdc_burgers(): +def test_burgers(): N = 128 nu = 0.2 T: float = 4.0 + t: float = 0.0 # Starting point in time domain_size: float = 2.0 * np.pi - t = np.linspace(0, domain_size, N, endpoint=False) burgers: Burgers = Burgers(N=N, nu=nu, domain_size=domain_size) u0 = burgers.initial_u0("sine") @@ -23,11 +24,11 @@ def test_sdc_burgers(): print("="*80) results = [] - u_analytical = burgers.analytical_integration(u0, t=T) + u_analytical = burgers.u_solution(u0, t=T) for nt in range(4): - num_timesteps = 2**nt * 1000 + num_timesteps = 2**nt * 500 print(f"Running simulation with num_timesteps={num_timesteps}") dt = T / num_timesteps @@ -37,18 +38,13 @@ def test_sdc_burgers(): u = u0.copy() - if time_integration == "rk1": - u = burgers.step_rk1_n(u, dt, num_timesteps) - - elif time_integration == "rk2": - u = burgers.step_rk2_n(u, dt, num_timesteps) - - elif time_integration == "rk4": - u = burgers.step_rk4_n(u, dt, num_timesteps) + if time_integration in RKIntegration.supported_methods: + rki = RKIntegration(method=time_integration) + u = rki.integrate_n(u, t, dt, num_timesteps, burgers) elif time_integration == "sdc": sdci = SDCIntegration(num_nodes=3, node_type="LEGENDRE", quad_type="LOBATTO") - u = sdci.integrate_n(u, dt, num_timesteps, burgers) + u = sdci.integrate_n(u, t, dt, num_timesteps, burgers) else: raise Exception("TODO") @@ -65,12 +61,18 @@ def test_sdc_burgers(): print(f" - N={r["N"]}, dt={r["dt"]:.6e}, error={r["error"]:.6e}, conv={conv}") prev_error = r["error"] + r["conv"] = conv if time_integration == "rk1": assert results[-1]["error"] < 1e-4 + assert np.abs(results[-1]["conv"]-1.0) < 1e-3 + elif time_integration == "rk2": - assert results[-1]["error"] < 1e-8 + assert results[-1]["error"] < 1e-7 + assert np.abs(results[-1]["conv"]-2.0) < 1e-3 + elif time_integration == "rk4": assert results[-1]["error"] < 1e-14 + elif time_integration == "sdc": assert results[-1]["error"] < 1e-14 diff --git a/qmat/playground/tests/test_dahlquist.py b/qmat/playground/tests/test_dahlquist.py new file mode 100644 index 0000000..0c44fec --- /dev/null +++ b/qmat/playground/tests/test_dahlquist.py @@ -0,0 +1,75 @@ +import numpy as np +from qmat.playground.diff_eqs.dahlquist import Dahlquist +from time_integration.sdc_integration import SDCIntegration +from qmat.playground.time_integration.rk_integration import RKIntegration +from matplotlib import pyplot as plt + + +def test_dahlquist2(): + u0 = np.array([1.0]) # Initial condition + T: float = 4 * np.pi # Time interval + T: float = 0.5 # Time interval + t: float = 0.0 # Starting time + + dahlquist: Dahlquist = Dahlquist(lam1=1.0j, lam2=1.0j) + + for time_integration in ["rk1", "rk2", "rk4", "sdc"]: + print("="*80) + print(f"Time integration method: {time_integration}") + print("="*80) + results = [] + + u_analytical = dahlquist.u_solution(u0, t=T) + + for nt in range(4): + + num_timesteps = 2**nt * 10 + print(f"Running simulation with num_timesteps={num_timesteps}") + + dt = T / num_timesteps + + u0 = dahlquist.initial_u0() + + u = u0.copy() + + if time_integration in RKIntegration.supported_methods: + rki = RKIntegration(method=time_integration) + + 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=1) + u = sdci.integrate_n(u, t, dt, num_timesteps, dahlquist) + + else: + raise Exception("TODO") + + error = np.max(np.abs(u - u_analytical)) + results.append({"N": num_timesteps, "dt": dt, "error": 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 + + 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 == "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 diff --git a/qmat/playground/tests/test_dahlquist2.py b/qmat/playground/tests/test_dahlquist2.py new file mode 100644 index 0000000..3504cc4 --- /dev/null +++ b/qmat/playground/tests/test_dahlquist2.py @@ -0,0 +1,76 @@ +import numpy as np +from qmat.playground.diff_eqs.dahlquist2 import Dahlquist2 +from qmat.playground.diff_eqs.dahlquist import Dahlquist +from time_integration.sdc_integration import SDCIntegration +from qmat.playground.time_integration.rk_integration import RKIntegration +from matplotlib import pyplot as plt + + +def test_dahlquist2(): + u0 = np.array([1.0]) # Initial condition + T: float = 4 * np.pi # Time interval + T: float = 0.5 # Time interval + t: float = 0.0 # Starting time + + dahlquist2: Dahlquist2 = Dahlquist2(lam1=1.0j, lam2=1.0j) + + for time_integration in ["rk1", "rk2", "rk4", "sdc"]: + print("="*80) + print(f"Time integration method: {time_integration}") + print("="*80) + results = [] + + u_analytical = dahlquist2.u_solution(u0, t=T) + + for nt in range(4): + + num_timesteps = 2**nt * 10 + print(f"Running simulation with num_timesteps={num_timesteps}") + + dt = T / num_timesteps + + u0 = dahlquist2.initial_u0() + + u = u0.copy() + + if time_integration in RKIntegration.supported_methods: + rki = RKIntegration(method=time_integration) + + u = rki.integrate_n(u, t, dt, num_timesteps, dahlquist2) + + 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, dahlquist2) + + else: + raise Exception("TODO") + + error = np.max(np.abs(u - u_analytical)) + results.append({"N": num_timesteps, "dt": dt, "error": 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 + + 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-5 + assert np.abs(results[-1]["conv"]-2.0) < 1e-3 + + elif time_integration == "rk4": + assert results[-1]["error"] < 1e-11 + assert np.abs(results[-1]["conv"]-4.0) < 1e-4 + + elif time_integration == "sdc": + assert results[-1]["error"] < 1e-5 + assert np.abs(results[-1]["conv"]-2.0) < 1e-3 diff --git a/qmat/playground/time_integration/rk_integration.py b/qmat/playground/time_integration/rk_integration.py new file mode 100644 index 0000000..ac77b35 --- /dev/null +++ b/qmat/playground/time_integration/rk_integration.py @@ -0,0 +1,49 @@ +from typing import List +import numpy as np +from qmat.playground.diff_eqs.de_solver import DESolver + + +class RKIntegration: + + # List of supported time integration methods + supported_methods: List[str] = ["rk1", "rk2", "rk4"] + + def __init__(self, method: str): + assert method in self.supported_methods, "Unsupported RK method" + self.method = method + + def integrate(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> np.array: + u = u0 + + if self.method == "rk1": + # RK1 (Forward Euler) + return u + dt * de_solver.du_dt(u, t) + + elif self.method == "rk2": + # RK2: Heun's method + k1 = de_solver.du_dt(u, t) + k2 = de_solver.du_dt(u + 0.5 * dt * k1, t + 0.5 * dt) + return u + dt * k2 + + elif self.method == "rk4": + # Classical RK4 method + k1 = de_solver.du_dt(u, t) + k2 = de_solver.du_dt(u + 0.5 * dt * k1, t + 0.5 * dt) + k3 = de_solver.du_dt(u + 0.5 * dt * k2, t + 0.5 * dt) + k4 = de_solver.du_dt(u + dt * k3, t + dt) + + # Update solution + return u + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6 + + else: + raise Exception("TODO") + + return u + + def integrate_n(self, u0: np.array, t: float, dt: float, num_timesteps, de_solver: DESolver) -> np.array: + u_value = u0 + + for n in range(num_timesteps): + u_value = self.integrate(u_value, t + n * dt, dt, de_solver) + + return u_value diff --git a/qmat/playground/time_integration/sdc_integration.py b/qmat/playground/time_integration/sdc_integration.py index a50f7d1..4de1b52 100644 --- a/qmat/playground/time_integration/sdc_integration.py +++ b/qmat/playground/time_integration/sdc_integration.py @@ -4,7 +4,7 @@ class SDCIntegration: def __init__( - self, num_nodes: int = 3, node_type: str = "LOBATTO", quad_type: str = "LOBATTO" + self, num_nodes: int = 3, node_type: str = "LOBATTO", quad_type: str = "LOBATTO", num_sweeps: int = None ): from qmat.qcoeff.collocation import Collocation import qmat.qdelta.timestepping as module @@ -22,7 +22,14 @@ def __init__( # Number of nodes self.N = len(self.nodes) - def integrate(self, u0: np.array, dt: float, de_solver: DESolver) -> np.array: + if num_sweeps is None: + self.num_sweeps = len(self.nodes) + else: + self.num_sweeps = num_sweeps + + assert self.num_sweeps >= 1 + + def integrate(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> np.array: if not np.isclose(self.nodes[0], 0.0): raise Exception("SDC nodes must include the left interval boundary.") @@ -31,44 +38,44 @@ def integrate(self, u0: np.array, dt: float, de_solver: DESolver) -> np.array: assert self.N == self.deltas.shape[0] - n_sweeps = 2 shape = (self.N,) + u0.shape - u = np.zeros(shape=shape, dtype=float) + u = np.zeros_like(u0, shape=shape) u[0, :] = u0 - eval_f_k0 = np.empty_like(u) - eval_f_k1 = np.empty_like(u) + du_dt_k0 = np.empty_like(u) + du_dt_k1 = np.empty_like(u) # Propagate initial condition to all nodes for m in range(0, self.N): if m > 0: - u[m] = u[m-1] + dt * self.deltas[m] * eval_f_k0[m-1] - eval_f_k0[m] = de_solver.eval_f(u[m]) + u[m] = u[m-1] + dt * self.deltas[m] * du_dt_k0[m-1] + du_dt_k0[m] = de_solver.du_dt(u[m], t + dt*self.nodes[m]) - if 1: - # Iteratively sweep over SDC nodes - for _ in range(n_sweeps): - for m in range(0, self.N): + # Iteratively sweep over SDC nodes + for _ in range(1, self.num_sweeps): + for m in range(0, self.N): - if m > 0: - qeval = self.q[m] @ eval_f_k0 - u[m] = u[m-1] + dt * (self.deltas[m] * (eval_f_k1[m-1] - eval_f_k0[m-1]) + qeval) + if m > 0: + qeval = self.q[m] @ du_dt_k0 + u[m] = u[m-1] + dt * (self.deltas[m] * (du_dt_k1[m-1] - du_dt_k0[m-1]) + qeval) - eval_f_k1[m] = de_solver.eval_f(u[m]) + du_dt_k1[m] = de_solver.du_dt(u[m], t + dt*self.nodes[m]) - # Copy tendency arrays - eval_f_k0[:] = eval_f_k1[:] + # Copy tendency arrays + du_dt_k0[:] = du_dt_k1[:] if 0: # If we're using Radau-right, we can just use the last value u[0] = u[-1] + else: # Compute new starting value with quadrature on tendencies - u[0] = u[0] + dt * self.weights @ eval_f_k0 + u[0] = u[0] + dt * self.weights @ du_dt_k0 + assert u0.shape == u[0].shape return u[0] - def integrate_n(self, u0: np.array, dt: float, num_timesteps, de_solver: DESolver) -> np.array: + def integrate_n(self, u0: np.array, t: float, dt: float, num_timesteps, de_solver: DESolver) -> np.array: if not np.isclose(self.nodes[0], 0.0): raise Exception("SDC nodes must include the left interval boundary.") @@ -77,7 +84,7 @@ def integrate_n(self, u0: np.array, dt: float, num_timesteps, de_solver: DESolve u_value = u0 - for _ in range(num_timesteps): - u_value = self.integrate(u_value, dt, de_solver) + for n in range(num_timesteps): + u_value = self.integrate(u_value, t + n * dt, dt, de_solver) return u_value From 1d334b78058a80d800eb4ecd1779d1b5b38338bd Mon Sep 17 00:00:00 2001 From: SCHREIBER Martin Date: Tue, 21 Oct 2025 17:18:46 +0200 Subject: [PATCH 04/11] further updates --- qmat/playground/diff_eqs/dahlquist.py | 2 +- qmat/playground/diff_eqs/dahlquist2.py | 27 ++++--- qmat/playground/diff_eqs/two_freq.py | 47 ++++++++++++ qmat/playground/ex_burgers_integrate.py | 2 +- qmat/playground/ex_dahlquist2_integrate.py | 27 ++++--- qmat/playground/ex_dahlquist2_plot.py | 15 ++-- qmat/playground/ex_dahlquist_plot.py | 18 +++++ qmat/playground/ex_two_freq_integrate.py | 87 ++++++++++++++++++++++ qmat/playground/ex_two_freq_plot.py | 19 +++++ qmat/playground/tests/test_burgers.py | 2 +- qmat/playground/tests/test_dahlquist.py | 5 +- qmat/playground/tests/test_dahlquist2.py | 6 +- qmat/playground/tests/test_two_freq.py | 74 ++++++++++++++++++ 13 files changed, 292 insertions(+), 39 deletions(-) create mode 100644 qmat/playground/diff_eqs/two_freq.py create mode 100644 qmat/playground/ex_dahlquist_plot.py create mode 100644 qmat/playground/ex_two_freq_integrate.py create mode 100644 qmat/playground/ex_two_freq_plot.py create mode 100644 qmat/playground/tests/test_two_freq.py diff --git a/qmat/playground/diff_eqs/dahlquist.py b/qmat/playground/diff_eqs/dahlquist.py index 321726d..73ae92f 100644 --- a/qmat/playground/diff_eqs/dahlquist.py +++ b/qmat/playground/diff_eqs/dahlquist.py @@ -16,7 +16,7 @@ def __init__(self, lam1: complex, lam2: complex): self.lam2: float = lam2 def initial_u0(self, mode: str = None) -> np.ndarray: - return np.array([1.0 + 0.0j]) + 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 diff --git a/qmat/playground/diff_eqs/dahlquist2.py b/qmat/playground/diff_eqs/dahlquist2.py index c6a10d1..9ed0ecc 100644 --- a/qmat/playground/diff_eqs/dahlquist2.py +++ b/qmat/playground/diff_eqs/dahlquist2.py @@ -10,20 +10,27 @@ class Dahlquist2(DESolver): u(t) = 0.5*(exp(lam1*t) + exp(lam2*t)) * u(0) """ - def __init__(self, lam1: complex, lam2: complex): - # Lambda 1 - self.lam1: float = lam1 - - # Lambda 2 - self.lam2: float = lam2 + 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]) + return np.array([1.0 + 0.0j], dtype=np.complex128) def du_dt(self, u: np.ndarray, t: float) -> np.ndarray: retval = ( - (self.lam1 * np.exp(t * self.lam1) + self.lam2 * np.exp(t * self.lam2)) - / (np.exp(t * self.lam1) + np.exp(t * self.lam2)) + (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 @@ -31,6 +38,6 @@ def du_dt(self, u: np.ndarray, t: float) -> np.ndarray: def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray: assert isinstance(t, float) - retval = 0.5*(np.exp(t * self.lam1) + np.exp(t * self.lam2)) * u0 + retval = (self.s * np.exp(t * self.lam1) + (1.0 - self.s) * np.exp(t * self.lam2)) * u0 assert retval.shape == u0.shape return retval diff --git a/qmat/playground/diff_eqs/two_freq.py b/qmat/playground/diff_eqs/two_freq.py new file mode 100644 index 0000000..658fc39 --- /dev/null +++ b/qmat/playground/diff_eqs/two_freq.py @@ -0,0 +1,47 @@ +import numpy as np +from qmat.playground.diff_eqs.de_solver import DESolver + + +class TwoFreq(DESolver): + """ + Test equation using a matrix to generate a + superposition of two frequencies in the solution. + """ + + def __init__(self, lam1: complex, lam2: complex, lam3: complex = None, lam4: complex = None): + """ + :param lam1: L[0,0] + :param lam2: L[0,1] + :param lam3: L[1,0] + :param lam4: L[1,1] + """ + + self.lam1: complex = lam1 + self.lam2: complex = lam2 + self.lam3: complex = lam3 if lam3 is not None else 0 + 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) + + # 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)): + raise Exception("Dahlquist3 matrix L must have purely imaginary eigenvalues") + + 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: + assert isinstance(t, float) + retval = self.L @ u + assert retval.shape == u.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 diff --git a/qmat/playground/ex_burgers_integrate.py b/qmat/playground/ex_burgers_integrate.py index 28caae4..886b5dd 100644 --- a/qmat/playground/ex_burgers_integrate.py +++ b/qmat/playground/ex_burgers_integrate.py @@ -53,7 +53,7 @@ u = sdci.integrate_n(u, t, dt, num_timesteps, burgers) else: - raise Exception("TODO") + raise Exception(f"Unsupported time integration method '{time_integration}'") plt.plot(t_, u, label=f"numerical {num_timesteps}", linestyle="dashed") diff --git a/qmat/playground/ex_dahlquist2_integrate.py b/qmat/playground/ex_dahlquist2_integrate.py index 8a8683b..c63d938 100644 --- a/qmat/playground/ex_dahlquist2_integrate.py +++ b/qmat/playground/ex_dahlquist2_integrate.py @@ -1,29 +1,30 @@ import numpy as np from qmat.playground.diff_eqs.dahlquist2 import Dahlquist2 -from qmat.playground.diff_eqs.dahlquist import Dahlquist +#from qmat.playground.diff_eqs.dahlquist import Dahlquist from time_integration.sdc_integration import SDCIntegration from qmat.playground.time_integration.rk_integration import RKIntegration from matplotlib import pyplot as plt u0 = np.array([1.0]) # Initial condition -T: float = 4 * np.pi # Time interval -T: float = 0.5 # Time interval +T: float = 2 * np.pi # Time interval t: float = 0.0 # Starting time time_integration = "rk4" +dahlquist2: Dahlquist2 = Dahlquist2(lam1=20.0j, lam2=1.0j, s=0.1) results = [] -for nt in range(4): +plt.close() - num_timesteps = 2**nt * 10 +for nt in range(1): + + num_timesteps = 2**nt * 300 print(f"Running simulation with num_timesteps={num_timesteps}") dt = T / num_timesteps - dahlquist2: Dahlquist2 = Dahlquist2(lam1=1.0j, lam2=1.0j) u_analytical_fin = dahlquist2.u_solution(u0, t=T) u0 = dahlquist2.initial_u0() @@ -38,7 +39,7 @@ for n in range(num_timesteps): u = rki.integrate(u, t + n * dt, dt, dahlquist2) - u_ = np.append(u_, u) + u_ = np.concatenate((u_, np.expand_dims(u, axis=0))) t_ = np.append(t_, t + (n + 1) * dt) elif time_integration == "sdc": @@ -47,11 +48,11 @@ for n in range(num_timesteps): u = sdci.integrate(u, t + n * dt, dt, dahlquist2) - u_ = np.append(u_, u) + u_ = np.concatenate((u_, np.expand_dims(u, axis=0))) t_ = np.append(t_, t + (n + 1) * dt) else: - raise Exception("TODO") + raise Exception(f"Unsupported time integration method '{time_integration}'") plt.plot(t_, np.real(u_), label=f"ndt={num_timesteps}, real", linestyle="dashed") plt.plot(t_, np.imag(u_), label=f"ndt={num_timesteps}, imag", linestyle="solid") @@ -70,9 +71,13 @@ prev_error = r["error"] if 1: + t_ = np.linspace(0, T, 1000) u_analytical_fin = np.array([dahlquist2.u_solution(u0, t) for t in t_]) - plt.plot(t_, np.real(u_analytical_fin), label="analytical, real", linestyle="dotted") - plt.plot(t_, np.imag(u_analytical_fin), label="analytical, imag", linestyle="dotted") + plt.plot(t_, np.real(u_analytical_fin), linestyle="dotted", color="gray") + plt.plot(t_, np.imag(u_analytical_fin), linestyle="dotted", color="gray") + +if 1: + plt.ylim(-1.5, 1.5) plt.legend() plt.show() diff --git a/qmat/playground/ex_dahlquist2_plot.py b/qmat/playground/ex_dahlquist2_plot.py index 9376ee8..a430c6a 100644 --- a/qmat/playground/ex_dahlquist2_plot.py +++ b/qmat/playground/ex_dahlquist2_plot.py @@ -1,4 +1,5 @@ import numpy as np +from matplotlib import pyplot as plt from qmat.playground.diff_eqs.dahlquist2 import Dahlquist2 @@ -6,14 +7,12 @@ u0 = np.array([1.0]) t: np.array = np.linspace(0, 4*np.pi, N, endpoint=False) -dahlquist = Dahlquist2(lam1=10.j, lam2=1.j) +dahlquist2: Dahlquist2 = Dahlquist2(lam1=20.0j, lam2=1.0j, s=0.1) -u_eval = np.array([dahlquist.u_solution(u0, _) for _ in t]) +u_eval = np.array([dahlquist2.u_solution(u0, _) for _ in t]) -if 1: - from matplotlib import pyplot as plt - plt.plot(t, np.real(u_eval), label="Re(u)") - plt.plot(t, np.imag(u_eval), label="Im(u)") - plt.legend() - plt.show() +plt.plot(t, np.real(u_eval), label="Re(u)") +plt.plot(t, np.imag(u_eval), label="Im(u)") +plt.legend() +plt.show() diff --git a/qmat/playground/ex_dahlquist_plot.py b/qmat/playground/ex_dahlquist_plot.py new file mode 100644 index 0000000..436b5de --- /dev/null +++ b/qmat/playground/ex_dahlquist_plot.py @@ -0,0 +1,18 @@ +import numpy as np +from matplotlib import pyplot as plt +from qmat.playground.diff_eqs.dahlquist import Dahlquist + + +N = 500 +u0 = np.array([1.0]) +t: np.array = np.linspace(0, 4*np.pi, N, endpoint=False) + +dahlquist: Dahlquist = Dahlquist(lam1=20.0j, lam2=1.0j) + +u_eval = np.array([dahlquist.u_solution(u0, _) for _ in t]) + + +plt.plot(t, np.real(u_eval), label="Re(u)") +plt.plot(t, np.imag(u_eval), label="Im(u)") +plt.legend() +plt.show() diff --git a/qmat/playground/ex_two_freq_integrate.py b/qmat/playground/ex_two_freq_integrate.py new file mode 100644 index 0000000..0fccc64 --- /dev/null +++ b/qmat/playground/ex_two_freq_integrate.py @@ -0,0 +1,87 @@ +import numpy as np +from qmat.playground.diff_eqs.two_freq import TwoFreq +from time_integration.sdc_integration import SDCIntegration +from qmat.playground.time_integration.rk_integration import RKIntegration +from matplotlib import pyplot as plt + + +T: float = 2 * 2 * np.pi # Time interval +t: float = 0.0 # Starting time + +time_integration = "sdc" + +two_freq: TwoFreq = TwoFreq(lam1=1.0j, lam2=20.0j, lam3=0.5j) +u0 = two_freq.initial_u0() + +results = [] + +plt.close() + +for nt in range(1): + + num_timesteps = 2**nt * 200 + print(f"Running simulation with num_timesteps={num_timesteps}") + + dt = T / num_timesteps + + u_analytical_fin = two_freq.u_solution(u0, t=T) + u0 = two_freq.initial_u0() + + u = u0.copy() + + u_: np.ndarray = np.array([u]) + t_: np.ndarray = np.array([t]) + + if time_integration in RKIntegration.supported_methods: + rki = RKIntegration(method=time_integration) + + for n in range(num_timesteps): + u = rki.integrate(u, t + n * dt, dt, two_freq) + + u_ = np.concatenate((u_, np.expand_dims(u, axis=0))) + 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) + + for n in range(num_timesteps): + u = sdci.integrate(u, t + n * dt, dt, two_freq) + + u_ = np.concatenate((u_, np.expand_dims(u, axis=0))) + t_ = np.append(t_, t + (n + 1) * dt) + + else: + raise Exception(f"Unsupported time integration method '{time_integration}'") + + for i in range(u_.shape[1]): + plt.plot(t_, np.real(u_[:, i]), label=f"u[{i}].real", linestyle="dashed") + plt.plot(t_, np.imag(u_[:, i]), label=f"u[{i}].imag", linestyle="solid") + + error = np.max(np.abs(u - u_analytical_fin)) + results.append({"N": num_timesteps, "dt": dt, "error": error}) + + +if 1: + # Plot analytical solution + t_ = np.linspace(0, T, 1000) + u_analytical_fin = np.array([two_freq.u_solution(u0, t) for t in t_]) + for i in range(u_analytical_fin.shape[1]): + plt.plot(t_, np.real(u_analytical_fin[:, i]), linestyle="dotted", color="black") + plt.plot(t_, np.imag(u_analytical_fin[:, i]), linestyle="dotted", color="black") + + +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"] + +if 1: + # plt.ylim(-1.5, 1.5) + + plt.legend() + plt.show() diff --git a/qmat/playground/ex_two_freq_plot.py b/qmat/playground/ex_two_freq_plot.py new file mode 100644 index 0000000..6429541 --- /dev/null +++ b/qmat/playground/ex_two_freq_plot.py @@ -0,0 +1,19 @@ +import numpy as np +from matplotlib import pyplot as plt +from qmat.playground.diff_eqs.two_freq import TwoFreq + + +N = 500 +t: np.array = np.linspace(0, 4*np.pi, N, endpoint=False) + +two_freq: TwoFreq = TwoFreq(lam1=1.0j, lam2=20.0j, lam3=0.5j) + +u0 = two_freq.initial_u0() +u_eval = np.array([two_freq.u_solution(u0, _) for _ in t]) + +for i in range(2): + plt.plot(t, np.real(u_eval[:, i]), label=f"Re(u[{i}])") + plt.plot(t, np.imag(u_eval[:, i]), label=f"Im(u[{i}])") + +plt.legend() +plt.show() diff --git a/qmat/playground/tests/test_burgers.py b/qmat/playground/tests/test_burgers.py index 4f4ec1a..09e2195 100644 --- a/qmat/playground/tests/test_burgers.py +++ b/qmat/playground/tests/test_burgers.py @@ -26,7 +26,7 @@ def test_burgers(): u_analytical = burgers.u_solution(u0, t=T) - for nt in range(4): + for nt in range(2, 4): num_timesteps = 2**nt * 500 print(f"Running simulation with num_timesteps={num_timesteps}") diff --git a/qmat/playground/tests/test_dahlquist.py b/qmat/playground/tests/test_dahlquist.py index 0c44fec..70bd1c1 100644 --- a/qmat/playground/tests/test_dahlquist.py +++ b/qmat/playground/tests/test_dahlquist.py @@ -2,16 +2,15 @@ from qmat.playground.diff_eqs.dahlquist import Dahlquist from time_integration.sdc_integration import SDCIntegration from qmat.playground.time_integration.rk_integration import RKIntegration -from matplotlib import pyplot as plt -def test_dahlquist2(): +def test_dahlquist(): u0 = np.array([1.0]) # Initial condition T: float = 4 * np.pi # Time interval T: float = 0.5 # Time interval t: float = 0.0 # Starting time - dahlquist: Dahlquist = Dahlquist(lam1=1.0j, lam2=1.0j) + dahlquist: Dahlquist = Dahlquist(lam1=1.0j, lam2=0.1j) for time_integration in ["rk1", "rk2", "rk4", "sdc"]: print("="*80) diff --git a/qmat/playground/tests/test_dahlquist2.py b/qmat/playground/tests/test_dahlquist2.py index 3504cc4..d2246b5 100644 --- a/qmat/playground/tests/test_dahlquist2.py +++ b/qmat/playground/tests/test_dahlquist2.py @@ -1,9 +1,7 @@ import numpy as np from qmat.playground.diff_eqs.dahlquist2 import Dahlquist2 -from qmat.playground.diff_eqs.dahlquist import Dahlquist from time_integration.sdc_integration import SDCIntegration from qmat.playground.time_integration.rk_integration import RKIntegration -from matplotlib import pyplot as plt def test_dahlquist2(): @@ -12,7 +10,7 @@ def test_dahlquist2(): T: float = 0.5 # Time interval t: float = 0.0 # Starting time - dahlquist2: Dahlquist2 = Dahlquist2(lam1=1.0j, lam2=1.0j) + dahlquist2: Dahlquist2 = Dahlquist2(lam1=1.0j, lam2=0.1j) for time_integration in ["rk1", "rk2", "rk4", "sdc"]: print("="*80) @@ -69,7 +67,7 @@ def test_dahlquist2(): elif time_integration == "rk4": assert results[-1]["error"] < 1e-11 - assert np.abs(results[-1]["conv"]-4.0) < 1e-4 + assert np.abs(results[-1]["conv"]-4.0) < 1e-2 elif time_integration == "sdc": assert results[-1]["error"] < 1e-5 diff --git a/qmat/playground/tests/test_two_freq.py b/qmat/playground/tests/test_two_freq.py new file mode 100644 index 0000000..12df2c5 --- /dev/null +++ b/qmat/playground/tests/test_two_freq.py @@ -0,0 +1,74 @@ +import numpy as np +from qmat.playground.diff_eqs.two_freq import TwoFreq +from time_integration.sdc_integration import SDCIntegration +from qmat.playground.time_integration.rk_integration import RKIntegration + + +def test_dahlquist2(): + T: float = 4 * np.pi # Time interval + T: float = 0.5 # Time interval + t: float = 0.0 # Starting time + + two_freq: TwoFreq = TwoFreq(lam1=1.0j, lam2=0.1j) + u0 = two_freq.initial_u0() + + for time_integration in ["rk1", "rk2", "rk4", "sdc"]: + print("="*80) + print(f"Time integration method: {time_integration}") + print("="*80) + results = [] + + u_analytical = two_freq.u_solution(u0, t=T) + + for nt in range(4): + + num_timesteps = 2**nt * 10 + print(f"Running simulation with num_timesteps={num_timesteps}") + + dt = T / num_timesteps + + u0 = two_freq.initial_u0() + + u = u0.copy() + + if time_integration in RKIntegration.supported_methods: + rki = RKIntegration(method=time_integration) + + u = rki.integrate_n(u, t, dt, num_timesteps, two_freq) + + elif time_integration == "sdc": + sdci = SDCIntegration(num_nodes=3, node_type="LEGENDRE", quad_type="LOBATTO", num_sweeps=2) + u = sdci.integrate_n(u, t, dt, num_timesteps, two_freq) + + else: + raise Exception("TODO") + + error = np.max(np.abs(u - u_analytical)) + results.append({"N": num_timesteps, "dt": dt, "error": 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 + + 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-5 + assert np.abs(results[-1]["conv"]-2.0) < 1e-3 + + elif time_integration == "rk4": + assert results[-1]["error"] < 1e-11 + assert np.abs(results[-1]["conv"]-4.0) < 1e-2 + + elif time_integration == "sdc": + assert results[-1]["error"] < 1e-5 + assert np.abs(results[-1]["conv"]-3.0) < 1e-3 From ba9a00813a7aadd607607acf669861e2f272e5fe Mon Sep 17 00:00:00 2001 From: SCHREIBER Martin Date: Tue, 21 Oct 2025 17:22:25 +0200 Subject: [PATCH 05/11] cleanups --- qmat/playground/diff_eqs/__init__.py | 10 ---------- qmat/playgrounds/martin/diff_eqs/__init__.py | 10 ++++++++++ .../martin}/diff_eqs/burgers.py | 4 ++-- .../martin}/diff_eqs/dahlquist.py | 2 +- .../martin}/diff_eqs/dahlquist2.py | 2 +- .../martin}/diff_eqs/de_solver.py | 0 .../martin}/diff_eqs/two_freq.py | 2 +- .../martin}/ex_burgers_integrate.py | 6 +++--- .../martin}/ex_burgers_plot.py | 0 .../martin}/ex_dahlquist2_integrate.py | 6 +++--- .../martin}/ex_dahlquist2_plot.py | 2 +- .../martin}/ex_dahlquist_plot.py | 2 +- .../martin}/ex_two_freq_integrate.py | 4 ++-- .../martin}/ex_two_freq_plot.py | 2 +- .../martin}/tests/__init__.py | 0 .../martin}/tests/test_burgers.py | 6 +++--- .../martin}/tests/test_dahlquist.py | 4 ++-- .../martin}/tests/test_dahlquist2.py | 4 ++-- .../martin}/tests/test_two_freq.py | 4 ++-- .../martin}/time_integration/rk_integration.py | 2 +- .../martin}/time_integration/sdc_integration.py | 2 +- 21 files changed, 37 insertions(+), 37 deletions(-) delete mode 100644 qmat/playground/diff_eqs/__init__.py create mode 100644 qmat/playgrounds/martin/diff_eqs/__init__.py rename qmat/{playground => playgrounds/martin}/diff_eqs/burgers.py (97%) rename qmat/{playground => playgrounds/martin}/diff_eqs/dahlquist.py (92%) rename qmat/{playground => playgrounds/martin}/diff_eqs/dahlquist2.py (95%) rename qmat/{playground => playgrounds/martin}/diff_eqs/de_solver.py (100%) rename qmat/{playground => playgrounds/martin}/diff_eqs/two_freq.py (95%) rename qmat/{playground => playgrounds/martin}/ex_burgers_integrate.py (89%) rename qmat/{playground => playgrounds/martin}/ex_burgers_plot.py (100%) rename qmat/{playground => playgrounds/martin}/ex_dahlquist2_integrate.py (91%) rename qmat/{playground => playgrounds/martin}/ex_dahlquist2_plot.py (85%) rename qmat/{playground => playgrounds/martin}/ex_dahlquist_plot.py (85%) rename qmat/{playground => playgrounds/martin}/ex_two_freq_integrate.py (94%) rename qmat/{playground => playgrounds/martin}/ex_two_freq_plot.py (87%) rename qmat/{playground => playgrounds/martin}/tests/__init__.py (100%) rename qmat/{playground => playgrounds/martin}/tests/test_burgers.py (90%) rename qmat/{playground => playgrounds/martin}/tests/test_dahlquist.py (94%) rename qmat/{playground => playgrounds/martin}/tests/test_dahlquist2.py (94%) rename qmat/{playground => playgrounds/martin}/tests/test_two_freq.py (94%) rename qmat/{playground => playgrounds/martin}/time_integration/rk_integration.py (95%) rename qmat/{playground => playgrounds/martin}/time_integration/sdc_integration.py (97%) diff --git a/qmat/playground/diff_eqs/__init__.py b/qmat/playground/diff_eqs/__init__.py deleted file mode 100644 index 77a6843..0000000 --- a/qmat/playground/diff_eqs/__init__.py +++ /dev/null @@ -1,10 +0,0 @@ - -from qmat.playground.diff_eqs.de_solver import DESolver -from qmat.playground.diff_eqs.burgers import Burgers -from qmat.playground.diff_eqs.dahlquist2 import Dahlquist2 - -__all__ = [ - "DESolver", - "Burgers", - "Dahlquist2", -] diff --git a/qmat/playgrounds/martin/diff_eqs/__init__.py b/qmat/playgrounds/martin/diff_eqs/__init__.py new file mode 100644 index 0000000..d33d182 --- /dev/null +++ b/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", +] diff --git a/qmat/playground/diff_eqs/burgers.py b/qmat/playgrounds/martin/diff_eqs/burgers.py similarity index 97% rename from qmat/playground/diff_eqs/burgers.py rename to qmat/playgrounds/martin/diff_eqs/burgers.py index 118c602..7e8f42b 100644 --- a/qmat/playground/diff_eqs/burgers.py +++ b/qmat/playgrounds/martin/diff_eqs/burgers.py @@ -1,6 +1,6 @@ import numpy as np -from qmat.playground.diff_eqs.de_solver import DESolver -from qmat.playground.time_integration.rk_integration import RKIntegration +from qmat.playgrounds.martin.diff_eqs.de_solver import DESolver +from qmat.playgrounds.martin.time_integration.rk_integration import RKIntegration class Burgers(DESolver): diff --git a/qmat/playground/diff_eqs/dahlquist.py b/qmat/playgrounds/martin/diff_eqs/dahlquist.py similarity index 92% rename from qmat/playground/diff_eqs/dahlquist.py rename to qmat/playgrounds/martin/diff_eqs/dahlquist.py index 73ae92f..74a77dd 100644 --- a/qmat/playground/diff_eqs/dahlquist.py +++ b/qmat/playgrounds/martin/diff_eqs/dahlquist.py @@ -1,5 +1,5 @@ import numpy as np -from qmat.playground.diff_eqs.de_solver import DESolver +from qmat.playgrounds.martin.diff_eqs.de_solver import DESolver class Dahlquist(DESolver): diff --git a/qmat/playground/diff_eqs/dahlquist2.py b/qmat/playgrounds/martin/diff_eqs/dahlquist2.py similarity index 95% rename from qmat/playground/diff_eqs/dahlquist2.py rename to qmat/playgrounds/martin/diff_eqs/dahlquist2.py index 9ed0ecc..de435c5 100644 --- a/qmat/playground/diff_eqs/dahlquist2.py +++ b/qmat/playgrounds/martin/diff_eqs/dahlquist2.py @@ -1,5 +1,5 @@ import numpy as np -from qmat.playground.diff_eqs.de_solver import DESolver +from qmat.playgrounds.martin.diff_eqs.de_solver import DESolver class Dahlquist2(DESolver): diff --git a/qmat/playground/diff_eqs/de_solver.py b/qmat/playgrounds/martin/diff_eqs/de_solver.py similarity index 100% rename from qmat/playground/diff_eqs/de_solver.py rename to qmat/playgrounds/martin/diff_eqs/de_solver.py diff --git a/qmat/playground/diff_eqs/two_freq.py b/qmat/playgrounds/martin/diff_eqs/two_freq.py similarity index 95% rename from qmat/playground/diff_eqs/two_freq.py rename to qmat/playgrounds/martin/diff_eqs/two_freq.py index 658fc39..a7a4a37 100644 --- a/qmat/playground/diff_eqs/two_freq.py +++ b/qmat/playgrounds/martin/diff_eqs/two_freq.py @@ -1,5 +1,5 @@ import numpy as np -from qmat.playground.diff_eqs.de_solver import DESolver +from qmat.playgrounds.martin.diff_eqs.de_solver import DESolver class TwoFreq(DESolver): diff --git a/qmat/playground/ex_burgers_integrate.py b/qmat/playgrounds/martin/ex_burgers_integrate.py similarity index 89% rename from qmat/playground/ex_burgers_integrate.py rename to qmat/playgrounds/martin/ex_burgers_integrate.py index 886b5dd..d92d141 100644 --- a/qmat/playground/ex_burgers_integrate.py +++ b/qmat/playgrounds/martin/ex_burgers_integrate.py @@ -4,9 +4,9 @@ import numpy as np -from qmat.playground.diff_eqs.burgers import Burgers -from qmat.playground.time_integration.sdc_integration import SDCIntegration -from qmat.playground.time_integration.rk_integration import RKIntegration +from qmat.playgrounds.martin.diff_eqs.burgers import Burgers +from qmat.playgrounds.martin.time_integration.sdc_integration import SDCIntegration +from qmat.playgrounds.martin.time_integration.rk_integration import RKIntegration from matplotlib import pyplot as plt diff --git a/qmat/playground/ex_burgers_plot.py b/qmat/playgrounds/martin/ex_burgers_plot.py similarity index 100% rename from qmat/playground/ex_burgers_plot.py rename to qmat/playgrounds/martin/ex_burgers_plot.py diff --git a/qmat/playground/ex_dahlquist2_integrate.py b/qmat/playgrounds/martin/ex_dahlquist2_integrate.py similarity index 91% rename from qmat/playground/ex_dahlquist2_integrate.py rename to qmat/playgrounds/martin/ex_dahlquist2_integrate.py index c63d938..f453b0a 100644 --- a/qmat/playground/ex_dahlquist2_integrate.py +++ b/qmat/playgrounds/martin/ex_dahlquist2_integrate.py @@ -1,8 +1,8 @@ import numpy as np -from qmat.playground.diff_eqs.dahlquist2 import Dahlquist2 -#from qmat.playground.diff_eqs.dahlquist import Dahlquist +from qmat.playgrounds.martin.diff_eqs.dahlquist2 import Dahlquist2 +#from qmat.playgrounds.martin.diff_eqs.dahlquist import Dahlquist from time_integration.sdc_integration import SDCIntegration -from qmat.playground.time_integration.rk_integration import RKIntegration +from qmat.playgrounds.martin.time_integration.rk_integration import RKIntegration from matplotlib import pyplot as plt diff --git a/qmat/playground/ex_dahlquist2_plot.py b/qmat/playgrounds/martin/ex_dahlquist2_plot.py similarity index 85% rename from qmat/playground/ex_dahlquist2_plot.py rename to qmat/playgrounds/martin/ex_dahlquist2_plot.py index a430c6a..168c56a 100644 --- a/qmat/playground/ex_dahlquist2_plot.py +++ b/qmat/playgrounds/martin/ex_dahlquist2_plot.py @@ -1,6 +1,6 @@ import numpy as np from matplotlib import pyplot as plt -from qmat.playground.diff_eqs.dahlquist2 import Dahlquist2 +from qmat.playgrounds.martin.diff_eqs.dahlquist2 import Dahlquist2 N = 500 diff --git a/qmat/playground/ex_dahlquist_plot.py b/qmat/playgrounds/martin/ex_dahlquist_plot.py similarity index 85% rename from qmat/playground/ex_dahlquist_plot.py rename to qmat/playgrounds/martin/ex_dahlquist_plot.py index 436b5de..232a22b 100644 --- a/qmat/playground/ex_dahlquist_plot.py +++ b/qmat/playgrounds/martin/ex_dahlquist_plot.py @@ -1,6 +1,6 @@ import numpy as np from matplotlib import pyplot as plt -from qmat.playground.diff_eqs.dahlquist import Dahlquist +from qmat.playgrounds.martin.diff_eqs.dahlquist import Dahlquist N = 500 diff --git a/qmat/playground/ex_two_freq_integrate.py b/qmat/playgrounds/martin/ex_two_freq_integrate.py similarity index 94% rename from qmat/playground/ex_two_freq_integrate.py rename to qmat/playgrounds/martin/ex_two_freq_integrate.py index 0fccc64..92b1b45 100644 --- a/qmat/playground/ex_two_freq_integrate.py +++ b/qmat/playgrounds/martin/ex_two_freq_integrate.py @@ -1,7 +1,7 @@ import numpy as np -from qmat.playground.diff_eqs.two_freq import TwoFreq +from qmat.playgrounds.martin.diff_eqs.two_freq import TwoFreq from time_integration.sdc_integration import SDCIntegration -from qmat.playground.time_integration.rk_integration import RKIntegration +from qmat.playgrounds.martin.time_integration.rk_integration import RKIntegration from matplotlib import pyplot as plt diff --git a/qmat/playground/ex_two_freq_plot.py b/qmat/playgrounds/martin/ex_two_freq_plot.py similarity index 87% rename from qmat/playground/ex_two_freq_plot.py rename to qmat/playgrounds/martin/ex_two_freq_plot.py index 6429541..54df8ce 100644 --- a/qmat/playground/ex_two_freq_plot.py +++ b/qmat/playgrounds/martin/ex_two_freq_plot.py @@ -1,6 +1,6 @@ import numpy as np from matplotlib import pyplot as plt -from qmat.playground.diff_eqs.two_freq import TwoFreq +from qmat.playgrounds.martin.diff_eqs.two_freq import TwoFreq N = 500 diff --git a/qmat/playground/tests/__init__.py b/qmat/playgrounds/martin/tests/__init__.py similarity index 100% rename from qmat/playground/tests/__init__.py rename to qmat/playgrounds/martin/tests/__init__.py diff --git a/qmat/playground/tests/test_burgers.py b/qmat/playgrounds/martin/tests/test_burgers.py similarity index 90% rename from qmat/playground/tests/test_burgers.py rename to qmat/playgrounds/martin/tests/test_burgers.py index 09e2195..de22a34 100644 --- a/qmat/playground/tests/test_burgers.py +++ b/qmat/playgrounds/martin/tests/test_burgers.py @@ -1,8 +1,8 @@ import numpy as np -from qmat.playground.diff_eqs.burgers import Burgers -from qmat.playground.time_integration.sdc_integration import SDCIntegration -from qmat.playground.time_integration.rk_integration import RKIntegration +from qmat.playgrounds.martin.diff_eqs.burgers import Burgers +from qmat.playgrounds.martin.time_integration.sdc_integration import SDCIntegration +from qmat.playgrounds.martin.time_integration.rk_integration import RKIntegration def test_burgers(): diff --git a/qmat/playground/tests/test_dahlquist.py b/qmat/playgrounds/martin/tests/test_dahlquist.py similarity index 94% rename from qmat/playground/tests/test_dahlquist.py rename to qmat/playgrounds/martin/tests/test_dahlquist.py index 70bd1c1..cd57623 100644 --- a/qmat/playground/tests/test_dahlquist.py +++ b/qmat/playgrounds/martin/tests/test_dahlquist.py @@ -1,7 +1,7 @@ import numpy as np -from qmat.playground.diff_eqs.dahlquist import Dahlquist +from qmat.playgrounds.martin.diff_eqs.dahlquist import Dahlquist from time_integration.sdc_integration import SDCIntegration -from qmat.playground.time_integration.rk_integration import RKIntegration +from qmat.playgrounds.martin.time_integration.rk_integration import RKIntegration def test_dahlquist(): diff --git a/qmat/playground/tests/test_dahlquist2.py b/qmat/playgrounds/martin/tests/test_dahlquist2.py similarity index 94% rename from qmat/playground/tests/test_dahlquist2.py rename to qmat/playgrounds/martin/tests/test_dahlquist2.py index d2246b5..f022ba7 100644 --- a/qmat/playground/tests/test_dahlquist2.py +++ b/qmat/playgrounds/martin/tests/test_dahlquist2.py @@ -1,7 +1,7 @@ import numpy as np -from qmat.playground.diff_eqs.dahlquist2 import Dahlquist2 +from qmat.playgrounds.martin.diff_eqs.dahlquist2 import Dahlquist2 from time_integration.sdc_integration import SDCIntegration -from qmat.playground.time_integration.rk_integration import RKIntegration +from qmat.playgrounds.martin.time_integration.rk_integration import RKIntegration def test_dahlquist2(): diff --git a/qmat/playground/tests/test_two_freq.py b/qmat/playgrounds/martin/tests/test_two_freq.py similarity index 94% rename from qmat/playground/tests/test_two_freq.py rename to qmat/playgrounds/martin/tests/test_two_freq.py index 12df2c5..597f845 100644 --- a/qmat/playground/tests/test_two_freq.py +++ b/qmat/playgrounds/martin/tests/test_two_freq.py @@ -1,7 +1,7 @@ import numpy as np -from qmat.playground.diff_eqs.two_freq import TwoFreq +from qmat.playgrounds.martin.diff_eqs.two_freq import TwoFreq from time_integration.sdc_integration import SDCIntegration -from qmat.playground.time_integration.rk_integration import RKIntegration +from qmat.playgrounds.martin.time_integration.rk_integration import RKIntegration def test_dahlquist2(): diff --git a/qmat/playground/time_integration/rk_integration.py b/qmat/playgrounds/martin/time_integration/rk_integration.py similarity index 95% rename from qmat/playground/time_integration/rk_integration.py rename to qmat/playgrounds/martin/time_integration/rk_integration.py index ac77b35..017a592 100644 --- a/qmat/playground/time_integration/rk_integration.py +++ b/qmat/playgrounds/martin/time_integration/rk_integration.py @@ -1,6 +1,6 @@ from typing import List import numpy as np -from qmat.playground.diff_eqs.de_solver import DESolver +from qmat.playgrounds.martin.diff_eqs.de_solver import DESolver class RKIntegration: diff --git a/qmat/playground/time_integration/sdc_integration.py b/qmat/playgrounds/martin/time_integration/sdc_integration.py similarity index 97% rename from qmat/playground/time_integration/sdc_integration.py rename to qmat/playgrounds/martin/time_integration/sdc_integration.py index 4de1b52..890ba01 100644 --- a/qmat/playground/time_integration/sdc_integration.py +++ b/qmat/playgrounds/martin/time_integration/sdc_integration.py @@ -1,5 +1,5 @@ import numpy as np -from qmat.playground.diff_eqs.de_solver import DESolver +from qmat.playgrounds.martin.diff_eqs.de_solver import DESolver class SDCIntegration: From 21e0e5b76365470816210cb769b99a329889782a Mon Sep 17 00:00:00 2001 From: SCHREIBER Martin Date: Wed, 22 Oct 2025 13:13:02 +0200 Subject: [PATCH 06/11] u --- qmat/playgrounds/martin/diff_eqs/burgers.py | 9 +- qmat/playgrounds/martin/diff_eqs/dahlquist.py | 35 ++++- .../playgrounds/martin/diff_eqs/dahlquist2.py | 2 +- qmat/playgrounds/martin/diff_eqs/de_solver.py | 69 ++++++++- qmat/playgrounds/martin/diff_eqs/two_freq.py | 40 ++++- .../martin/ex_two_freq_integrate.py | 8 +- .../martin/tests/test_dahlquist.py | 96 ++++++------ .../martin/tests/test_dahlquist_ext_freq.py | 94 ++++++++++++ .../playgrounds/martin/tests/test_two_freq.py | 96 ++++++------ .../martin/time_integration/rk_integration.py | 28 +++- .../time_integration/sdc_integration.py | 140 ++++++++++++++++-- 11 files changed, 495 insertions(+), 122 deletions(-) create mode 100644 qmat/playgrounds/martin/tests/test_dahlquist_ext_freq.py diff --git a/qmat/playgrounds/martin/diff_eqs/burgers.py b/qmat/playgrounds/martin/diff_eqs/burgers.py index 7e8f42b..786eda2 100644 --- a/qmat/playgrounds/martin/diff_eqs/burgers.py +++ b/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 diff --git a/qmat/playgrounds/martin/diff_eqs/dahlquist.py b/qmat/playgrounds/martin/diff_eqs/dahlquist.py index 74a77dd..377bbe3 100644 --- a/qmat/playgrounds/martin/diff_eqs/dahlquist.py +++ b/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 diff --git a/qmat/playgrounds/martin/diff_eqs/dahlquist2.py b/qmat/playgrounds/martin/diff_eqs/dahlquist2.py index de435c5..643d900 100644 --- a/qmat/playgrounds/martin/diff_eqs/dahlquist2.py +++ b/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)) diff --git a/qmat/playgrounds/martin/diff_eqs/de_solver.py b/qmat/playgrounds/martin/diff_eqs/de_solver.py index b057554..efca7de 100644 --- a/qmat/playgrounds/martin/diff_eqs/de_solver.py +++ b/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 diff --git a/qmat/playgrounds/martin/diff_eqs/two_freq.py b/qmat/playgrounds/martin/diff_eqs/two_freq.py index a7a4a37..776059d 100644 --- a/qmat/playgrounds/martin/diff_eqs/two_freq.py +++ b/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,12 +36,26 @@ 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 @@ -45,3 +63,21 @@ def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray: 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 diff --git a/qmat/playgrounds/martin/ex_two_freq_integrate.py b/qmat/playgrounds/martin/ex_two_freq_integrate.py index 92b1b45..d11bf39 100644 --- a/qmat/playgrounds/martin/ex_two_freq_integrate.py +++ b/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) diff --git a/qmat/playgrounds/martin/tests/test_dahlquist.py b/qmat/playgrounds/martin/tests/test_dahlquist.py index cd57623..7abaebb 100644 --- a/qmat/playgrounds/martin/tests/test_dahlquist.py +++ b/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 diff --git a/qmat/playgrounds/martin/tests/test_dahlquist_ext_freq.py b/qmat/playgrounds/martin/tests/test_dahlquist_ext_freq.py new file mode 100644 index 0000000..93ee8a7 --- /dev/null +++ b/qmat/playgrounds/martin/tests/test_dahlquist_ext_freq.py @@ -0,0 +1,94 @@ +import numpy as np +from qmat.playgrounds.martin.diff_eqs.dahlquist import Dahlquist +from time_integration.sdc_integration import SDCIntegration +from qmat.playgrounds.martin.time_integration.rk_integration import RKIntegration + + +def test_dahlquist(): + u0 = np.array([1.0]) # Initial condition + T: float = 4 * np.pi # Time interval + T: float = 0.5 # Time interval + t: float = 0.0 # Starting time + + dahlquist: Dahlquist = Dahlquist(lam1=1.0j, lam2=0.1j, ext_scalar=1.0) + + for time_integration in ["rk1", "rk2", "rk4", "irk1", "irk2", "sdc"]: + if time_integration == "sdc": + micro_time_integration_ = ["erk1", "irk1"] + else: + micro_time_integration_ = ["-"] + + for micro_time_integration in micro_time_integration_: + print("=" * 80) + print(f"Time integration method: {time_integration} ({micro_time_integration})") + print("=" * 80) + results = [] + + u_analytical = dahlquist.u_solution(u0, t=T) + + for nt in range(4): + + num_timesteps = 2**nt * 10 + print(f"Running simulation with num_timesteps={num_timesteps}") + + dt = T / num_timesteps + + u0 = dahlquist.initial_u0() + + u = u0.copy() + + if time_integration in RKIntegration.supported_methods: + rki = RKIntegration(method=time_integration) + + 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") + + 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"]) + + 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 + + elif time_integration == "rk2": + assert results[-1]["error"] < 1e-4 + assert np.abs(results[-1]["conv"] - 2.0) < 1e-2 + + elif time_integration == "rk4": + assert results[-1]["error"] < 1e-9 + assert np.abs(results[-1]["conv"] - 4.0) < 1e-2 + + elif time_integration == "irk1": + assert results[-1]["error"] < 1e-2 + assert np.abs(results[-1]["conv"] - 1.0) < 1e-2 + + elif time_integration == "irk2": + assert results[-1]["error"] < 1e-2 + assert np.abs(results[-1]["conv"] - 2.0) < 1e-2 + + elif time_integration == "sdc": + assert results[-1]["error"] < 1e-8 + assert np.abs(results[-1]["conv"] - 4.0) < 1e-1 diff --git a/qmat/playgrounds/martin/tests/test_two_freq.py b/qmat/playgrounds/martin/tests/test_two_freq.py index 597f845..4842dfd 100644 --- a/qmat/playgrounds/martin/tests/test_two_freq.py +++ b/qmat/playgrounds/martin/tests/test_two_freq.py @@ -13,62 +13,74 @@ def test_dahlquist2(): u0 = two_freq.initial_u0() 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 = two_freq.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 = two_freq.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 = two_freq.initial_u0() + dt = T / num_timesteps - u = u0.copy() + u0 = two_freq.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, two_freq) + 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=2) - u = sdci.integrate_n(u, t, dt, num_timesteps, two_freq) + u = rki.integrate_n(u, t, dt, num_timesteps, two_freq) + + 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, two_freq) - 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}) - 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-5 - assert np.abs(results[-1]["conv"]-2.0) < 1e-3 + elif time_integration == "rk2": + assert results[-1]["error"] < 1e-5 + assert np.abs(results[-1]["conv"] - 2.0) < 1e-3 - elif time_integration == "rk4": - assert results[-1]["error"] < 1e-11 - assert np.abs(results[-1]["conv"]-4.0) < 1e-2 + elif time_integration == "rk4": + assert results[-1]["error"] < 1e-11 + assert np.abs(results[-1]["conv"] - 4.0) < 1e-2 - elif time_integration == "sdc": - assert results[-1]["error"] < 1e-5 - assert np.abs(results[-1]["conv"]-3.0) < 1e-3 + elif time_integration == "sdc": + assert results[-1]["error"] < 1e-8 + assert np.abs(results[-1]["conv"] - 4.0) < 1e-3 diff --git a/qmat/playgrounds/martin/time_integration/rk_integration.py b/qmat/playgrounds/martin/time_integration/rk_integration.py index 017a592..5ffde40 100644 --- a/qmat/playgrounds/martin/time_integration/rk_integration.py +++ b/qmat/playgrounds/martin/time_integration/rk_integration.py @@ -6,7 +6,7 @@ class RKIntegration: # List of supported time integration methods - supported_methods: List[str] = ["rk1", "rk2", "rk4"] + supported_methods: List[str] = ["rk1", "rk2", "rk4", "irk1", "irk2"] def __init__(self, method: str): assert method in self.supported_methods, "Unsupported RK method" @@ -17,24 +17,36 @@ def integrate(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> n if self.method == "rk1": # RK1 (Forward Euler) - return u + dt * de_solver.du_dt(u, t) + return u + dt * de_solver.evalF(u, t) elif self.method == "rk2": # RK2: Heun's method - k1 = de_solver.du_dt(u, t) - k2 = de_solver.du_dt(u + 0.5 * dt * k1, t + 0.5 * dt) + k1 = de_solver.evalF(u, t) + k2 = de_solver.evalF(u + 0.5 * dt * k1, t + 0.5 * dt) return u + dt * k2 elif self.method == "rk4": # Classical RK4 method - k1 = de_solver.du_dt(u, t) - k2 = de_solver.du_dt(u + 0.5 * dt * k1, t + 0.5 * dt) - k3 = de_solver.du_dt(u + 0.5 * dt * k2, t + 0.5 * dt) - k4 = de_solver.du_dt(u + dt * k3, t + dt) + k1 = de_solver.evalF(u, t) + k2 = de_solver.evalF(u + 0.5 * dt * k1, t + 0.5 * dt) + k3 = de_solver.evalF(u + 0.5 * dt * k2, t + 0.5 * dt) + k4 = de_solver.evalF(u + dt * k3, t + dt) # Update solution return u + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6 + elif self.method == "irk1": + # IRK1 (Implicit/backward Euler) + return de_solver.fSolve(u, dt, t) + + elif self.method == "irk2": + # IRK2 (Crank-Nicolson method) + # Implicit step + u = de_solver.fSolve(u, 0.5*dt, t) + # Forward step + u += 0.5 * dt * de_solver.evalF(u, t + 0.5*dt) + return u + else: raise Exception("TODO") diff --git a/qmat/playgrounds/martin/time_integration/sdc_integration.py b/qmat/playgrounds/martin/time_integration/sdc_integration.py index 890ba01..210b233 100644 --- a/qmat/playgrounds/martin/time_integration/sdc_integration.py +++ b/qmat/playgrounds/martin/time_integration/sdc_integration.py @@ -4,7 +4,12 @@ class SDCIntegration: def __init__( - self, num_nodes: int = 3, node_type: str = "LOBATTO", quad_type: str = "LOBATTO", num_sweeps: int = None + self, + num_nodes: int = 3, + node_type: str = "LOBATTO", # Basic node types to generate SDC nodes + quad_type: str = "LOBATTO", # 'LOBATTO': Always include {0, 1} in quadrature points. Add them if they don't exist. + num_sweeps: int = None, + micro_time_integration: str = None, # 'erk1' = explicit Euler, 'irk1' = implicit Euler, 'imex' = implicit-explicit ): from qmat.qcoeff.collocation import Collocation import qmat.qdelta.timestepping as module @@ -15,21 +20,35 @@ def __init__( self.nodes, self.weights, self.q = coll.genCoeffs(form="N2N") - self.q_delta: np.array = self.gen.getQDelta() - self.d_tau: np.array = self.gen.dTau - self.deltas: np.array = self.gen.deltas + self.q_delta: np.ndarray = self.gen.getQDelta() + # Deltas are the \tau + self.deltas: np.ndarray = self.gen.deltas # Number of nodes self.N = len(self.nodes) + print(f"self.nodes: {self.nodes}") + print(f"self.deltas: {self.deltas}") + if num_sweeps is None: self.num_sweeps = len(self.nodes) else: self.num_sweeps = num_sweeps + # Time integration to be used within SDC sweeps + # 'erk1' = explicit Euler + # 'irk1' = implicit Euler + # 'imex' = implicit-explicit => 'imex12' + # 'imex12' = implicit-explicit: 1st term treated implicitly, 2nd term explicitly + # 'imex21' = implicit-explicit: 2nd term treated implicitly, 1st term explicitly + self.time_integration_method = micro_time_integration if micro_time_integration is not None else "erk1" + + if self.time_integration_method == "imex": + self.time_integration_method = "imex12" + assert self.num_sweeps >= 1 - def integrate(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> np.array: + def integrate_erk1(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> np.array: if not np.isclose(self.nodes[0], 0.0): raise Exception("SDC nodes must include the left interval boundary.") @@ -42,27 +61,31 @@ def integrate(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> n u = np.zeros_like(u0, shape=shape) u[0, :] = u0 - du_dt_k0 = np.empty_like(u) - du_dt_k1 = np.empty_like(u) + evalF_k0 = np.empty_like(u) + evalF_k1 = np.empty_like(u) + # # Propagate initial condition to all nodes + # for m in range(0, self.N): if m > 0: - u[m] = u[m-1] + dt * self.deltas[m] * du_dt_k0[m-1] - du_dt_k0[m] = de_solver.du_dt(u[m], t + dt*self.nodes[m]) + u[m] = u[m - 1] + dt * self.deltas[m] * evalF_k0[m - 1] + evalF_k0[m] = de_solver.evalF(u[m], t + dt * self.nodes[m]) + # # Iteratively sweep over SDC nodes + # for _ in range(1, self.num_sweeps): for m in range(0, self.N): if m > 0: - qeval = self.q[m] @ du_dt_k0 - u[m] = u[m-1] + dt * (self.deltas[m] * (du_dt_k1[m-1] - du_dt_k0[m-1]) + qeval) + qeval = self.q[m] @ evalF_k0 + u[m] = u[m - 1] + dt * (self.deltas[m] * (evalF_k1[m - 1] - evalF_k0[m - 1]) + qeval) - du_dt_k1[m] = de_solver.du_dt(u[m], t + dt*self.nodes[m]) + evalF_k1[m] = de_solver.evalF(u[m], t + dt * self.nodes[m]) # Copy tendency arrays - du_dt_k0[:] = du_dt_k1[:] + evalF_k0[:] = evalF_k1[:] if 0: # If we're using Radau-right, we can just use the last value @@ -70,11 +93,100 @@ def integrate(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> n else: # Compute new starting value with quadrature on tendencies - u[0] = u[0] + dt * self.weights @ du_dt_k0 + u[0] = u[0] + dt * self.weights @ evalF_k0 assert u0.shape == u[0].shape return u[0] + def integrate_irk1(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> np.array: + if not np.isclose(self.nodes[0], 0.0): + raise Exception("SDC nodes must include the left interval boundary.") + + if not np.isclose(self.nodes[-1], 1.0): + raise Exception("SDC nodes must include the right interval boundary.") + + assert self.N == self.deltas.shape[0] + + shape = (self.N,) + u0.shape + u = np.zeros_like(u0, shape=shape) + + u[0, :] = u0 + evalF_k0 = np.empty_like(u) + evalF_k1 = np.empty_like(u) + + # Backup integrator contribution I[...] of previous iteration + ISolves = np.empty_like(u) + + # + # Propagate initial condition to all nodes + # + for m in range(0, self.N): + if m > 0: + # + # Solve implicit step: + # + # u^n+1 = u^n + dt * delta * F(u^n+1) + # <=> u^n+1 - dt * delta * F(u^n+1) = u^n + # <=> (I - dt * delta * F) * u^n+1 = u^n + # + rhs = u[m - 1] + u[m] = de_solver.fSolve(rhs, dt * self.deltas[m], t + dt * self.nodes[m]) + # Compute I[...] term + # u^n+1 = u^n + dt * delta * F(u^n+1) + # dt * delta * F(u^n+1) = u^n+1 - u^n + ISolves[m] = u[m] - u[m - 1] + + evalF_k0[m] = de_solver.evalF(u[m], t + dt * self.nodes[m]) + + # + # Iteratively sweep over SDC nodes + # + for _ in range(1, self.num_sweeps): + for m in range(0, self.N): + if m > 0: + # + # Solve implicit step: + # + # u^n+1 = u^n + dt * delta * (F(u^n+1)) - I(u^n) + dt * Q * F(u^n) + # <=> u^n+1 - dt * delta * F(u^n+1) = u^n - I(u^n) + dt * Q * F(u^n) + # <=> (I - dt * delta * F) * u^n+1 = u^n - I(u^n) + dt * Q * F(u^n) + # + # rhs = u^n - I(u^n) + dt * Q * F(u^n) + # + qeval = self.q[m] @ evalF_k0 + rhs = u[m - 1] - ISolves[m] + dt * qeval + + u[m] = de_solver.fSolve(rhs, dt * self.deltas[m], t + dt * self.nodes[m]) + + # Update I[...] term for next correction + # <=> - dt * delta * F(u^n+1) = u^n - I(u^n) + dt * Q * F(u^n) - u^n+1 + # <=> dt * delta * F(u^n+1) = u^n+1 - rhs + ISolves[m] = u[m] - rhs + + evalF_k1[m] = de_solver.evalF(u[m], t + dt * self.nodes[m]) + + # Copy tendency arrays + evalF_k0[:] = evalF_k1[:] + + if 0: + # If we're using Radau-right, we can just use the last value + u[0] = u[-1] + + else: + # Compute new starting value with quadrature on tendencies + u[0] = u[0] + dt * self.weights @ evalF_k0 + + assert u0.shape == u[0].shape + return u[0] + + def integrate(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> np.array: + if self.time_integration_method == "erk1": + return self.integrate_erk1(u0, t, dt, de_solver) + elif self.time_integration_method == "irk1": + return self.integrate_irk1(u0, t, dt, de_solver) + else: + raise Exception(f"Unsupported time integration within SDC: '{self.time_integration_method}'") + def integrate_n(self, u0: np.array, t: float, dt: float, num_timesteps, de_solver: DESolver) -> np.array: if not np.isclose(self.nodes[0], 0.0): raise Exception("SDC nodes must include the left interval boundary.") From d0da73d43b5e7115c2d36f4cc5823f5a49a65088 Mon Sep 17 00:00:00 2001 From: SCHREIBER Martin Date: Wed, 22 Oct 2025 13:18:10 +0200 Subject: [PATCH 07/11] updated test --- qmat/playgrounds/martin/tests/test_dahlquist.py | 3 --- qmat/playgrounds/martin/tests/test_two_freq.py | 10 +++++++++- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/qmat/playgrounds/martin/tests/test_dahlquist.py b/qmat/playgrounds/martin/tests/test_dahlquist.py index 7abaebb..e4f7a3c 100644 --- a/qmat/playgrounds/martin/tests/test_dahlquist.py +++ b/qmat/playgrounds/martin/tests/test_dahlquist.py @@ -32,14 +32,11 @@ def test_dahlquist(): print(f"Running simulation with num_timesteps={num_timesteps}") dt = T / num_timesteps - u0 = dahlquist.initial_u0() - u = u0.copy() if time_integration in RKIntegration.supported_methods: rki = RKIntegration(method=time_integration) - u = rki.integrate_n(u, t, dt, num_timesteps, dahlquist) elif time_integration == "sdc": diff --git a/qmat/playgrounds/martin/tests/test_two_freq.py b/qmat/playgrounds/martin/tests/test_two_freq.py index 4842dfd..c390ef5 100644 --- a/qmat/playgrounds/martin/tests/test_two_freq.py +++ b/qmat/playgrounds/martin/tests/test_two_freq.py @@ -12,7 +12,7 @@ def test_dahlquist2(): two_freq: TwoFreq = TwoFreq(lam1=1.0j, lam2=0.1j) u0 = two_freq.initial_u0() - for time_integration in ["rk1", "rk2", "rk4", "sdc"]: + for time_integration in ["rk1", "rk2", "rk4", "irk1", "irk2", "sdc"]: if time_integration == "sdc": micro_time_integration_ = ["erk1", "irk1"] else: @@ -81,6 +81,14 @@ def test_dahlquist2(): assert results[-1]["error"] < 1e-11 assert np.abs(results[-1]["conv"] - 4.0) < 1e-2 + elif time_integration == "irk1": + assert results[-1]["error"] < 1e-2 + assert np.abs(results[-1]["conv"] - 1.0) < 1e-2 + + elif time_integration == "irk2": + assert results[-1]["error"] < 1e-2 + assert np.abs(results[-1]["conv"] - 2.0) < 1e-2 + elif time_integration == "sdc": assert results[-1]["error"] < 1e-8 assert np.abs(results[-1]["conv"] - 4.0) < 1e-3 From f4305c7cf0ea3d3c67dfca01ab8730e586fb81a2 Mon Sep 17 00:00:00 2001 From: SCHREIBER Martin Date: Wed, 22 Oct 2025 14:45:07 +0200 Subject: [PATCH 08/11] u --- qmat/playgrounds/martin/diff_eqs/de_solver.py | 3 +-- qmat/playgrounds/martin/ex_two_freq_integrate.py | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/qmat/playgrounds/martin/diff_eqs/de_solver.py b/qmat/playgrounds/martin/diff_eqs/de_solver.py index efca7de..a329edb 100644 --- a/qmat/playgrounds/martin/diff_eqs/de_solver.py +++ b/qmat/playgrounds/martin/diff_eqs/de_solver.py @@ -22,7 +22,6 @@ def evalF(self, u: np.ndarray, t: float) -> np.ndarray: 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. @@ -47,7 +46,7 @@ def fSolve(self, rhs: np.ndarray, dt: float, t: float) -> np.ndarray: u_new : np.ndarray Array of shape (N,) representing the solution at the next time step. """ - pass + raise Exception("TODO: Implicit solver not implemented for this DE solver.") @abstractmethod def initial_u0(self, mode: str) -> np.ndarray: diff --git a/qmat/playgrounds/martin/ex_two_freq_integrate.py b/qmat/playgrounds/martin/ex_two_freq_integrate.py index d11bf39..0d9fd14 100644 --- a/qmat/playgrounds/martin/ex_two_freq_integrate.py +++ b/qmat/playgrounds/martin/ex_two_freq_integrate.py @@ -10,7 +10,7 @@ time_integration: str = "sdc" sdc_micro_time_integration: str = "irk1" -sdc_num_sweeps: int = 4 +sdc_num_sweeps: int = 2 two_freq: TwoFreq = TwoFreq(lam1=1.0j, lam2=20.0j, lam3=0.5j) u0 = two_freq.initial_u0() From 8512d59332e5fc59996a004412d8ad86e4450ff3 Mon Sep 17 00:00:00 2001 From: SCHREIBER Martin Date: Thu, 23 Oct 2025 11:26:34 +0200 Subject: [PATCH 09/11] added IMEX12 and IMEX21 --- qmat/playgrounds/martin/diff_eqs/dahlquist.py | 32 +++ qmat/playgrounds/martin/diff_eqs/de_solver.py | 89 +++++++-- qmat/playgrounds/martin/diff_eqs/two_freq.py | 46 +++-- .../martin/ex_two_freq_integrate.py | 9 +- qmat/playgrounds/martin/tests/test_burgers.py | 21 +- .../martin/tests/test_dahlquist.py | 25 ++- .../martin/tests/test_dahlquist2.py | 21 +- .../martin/tests/test_dahlquist_ext_freq.py | 27 +-- .../playgrounds/martin/tests/test_two_freq.py | 29 +-- .../time_integration/sdc_integration.py | 188 ++++++++++++++++++ 10 files changed, 395 insertions(+), 92 deletions(-) diff --git a/qmat/playgrounds/martin/diff_eqs/dahlquist.py b/qmat/playgrounds/martin/diff_eqs/dahlquist.py index 377bbe3..e18bffc 100644 --- a/qmat/playgrounds/martin/diff_eqs/dahlquist.py +++ b/qmat/playgrounds/martin/diff_eqs/dahlquist.py @@ -34,6 +34,21 @@ def evalF(self, u: np.ndarray, t: float) -> np.ndarray: 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 @@ -45,6 +60,23 @@ def fSolve(self, u: np.ndarray, dt: float, t: float) -> np.ndarray: 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 u_solution(self, u0: np.ndarray, t: float) -> np.ndarray: lam = self.lam1 + self.lam2 s = self.ext_scalar diff --git a/qmat/playgrounds/martin/diff_eqs/de_solver.py b/qmat/playgrounds/martin/diff_eqs/de_solver.py index a329edb..ceed5f0 100644 --- a/qmat/playgrounds/martin/diff_eqs/de_solver.py +++ b/qmat/playgrounds/martin/diff_eqs/de_solver.py @@ -6,7 +6,7 @@ class DESolver(ABC): @abstractmethod def evalF(self, u: np.ndarray, t: float) -> np.ndarray: - """Evaluate the right-hand side of the 1D viscous Burgers' equation. + """Evaluate the right-hand side of the equation. Parameters ---------- @@ -14,13 +14,33 @@ def evalF(self, u: np.ndarray, t: float) -> 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 + def evalF1(self, u: np.ndarray, t: float) -> np.ndarray: + """Evaluate the 1st term of the equation + + Parameters + ---------- + u : np.ndarray + Array of shape (N,) representing the solution at the current time step. + t : float + Current timestamp. + """ + raise Exception("TODO: `evalF1` not implemented for this DE solver.") + + def evalF2(self, u: np.ndarray, t: float) -> np.ndarray: + """Evaluate the 2nd term of the equation + + Parameters + ---------- + u : np.ndarray + Array of shape (N,) representing the solution at the current time step. + t : float + Current timestamp. + """ + raise Exception("TODO: `evalF2` not implemented for this DE solver.") + # This is optional since not every DE might have a solver for backward Euler def fSolve(self, rhs: np.ndarray, dt: float, t: float) -> np.ndarray: """Solve the right-hand side of an equation implicitly. @@ -40,27 +60,59 @@ def fSolve(self, rhs: np.ndarray, dt: float, t: float) -> np.ndarray: 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. """ raise Exception("TODO: Implicit solver not implemented for this DE solver.") + def fSolve1(self, rhs: np.ndarray, dt: float, t: float) -> np.ndarray: + """Solve the right-hand side of the 1st 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. + """ + raise Exception("TODO: `fSolve1` not implemented for this DE solver.") + + def fSolve2(self, rhs: np.ndarray, dt: float, t: float) -> np.ndarray: + """Solve the right-hand side of the 1st 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. + """ + raise Exception("TODO: `fSolve2` not implemented for this DE solver.") + @abstractmethod def initial_u0(self, mode: str) -> np.ndarray: - """Compute some initial conditions for the 1D viscous Burgers' equation. + """Compute some initial conditions for the equation. Parameters ---------- mode : str The type of initial condition to generate. - - Returns - ------- - u0 : np.ndarray - Array of shape (N,) representing the initial condition. """ pass @@ -75,10 +127,5 @@ def u_solution(self, u0: np.ndarray, t: float) -> 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 diff --git a/qmat/playgrounds/martin/diff_eqs/two_freq.py b/qmat/playgrounds/martin/diff_eqs/two_freq.py index 776059d..88c0ac5 100644 --- a/qmat/playgrounds/martin/diff_eqs/two_freq.py +++ b/qmat/playgrounds/martin/diff_eqs/two_freq.py @@ -42,6 +42,18 @@ def evalF(self, u: np.ndarray, t: float) -> np.ndarray: assert retval.shape == u.shape return retval + def evalF1(self, u: np.ndarray, t: float) -> np.ndarray: + 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: + assert isinstance(t, float) + retval = self.L2 @ u + assert retval.shape == u.shape + return retval + def fSolve(self, rhs: np.ndarray, dt: float, t: float) -> np.ndarray: """ Solve @@ -56,28 +68,32 @@ def fSolve(self, rhs: np.ndarray, dt: float, t: float) -> np.ndarray: assert retval.shape == rhs.shape return retval - def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray: - from scipy.linalg import expm - + def fSolve1(self, rhs: np.ndarray, dt: float, t: float) -> np.ndarray: + """ + Solve + `u_t = L1 u` + implicitly + """ assert isinstance(t, float) - retval = expm(self.L * t) @ u0 - assert retval.shape == u0.shape + retval = np.linalg.solve(np.eye(rhs.shape[0]) - dt*self.L1, rhs) + assert retval.shape == rhs.shape return retval - def evalF1(self, u: np.ndarray, t: float) -> np.ndarray: + def fSolve2(self, rhs: np.ndarray, dt: float, t: float) -> np.ndarray: """ - Solely compute tendencies of first frequency + Solve + `u_t = L2 u` + implicitly """ assert isinstance(t, float) - retval = self.L1 @ u - assert retval.shape == u.shape + retval = np.linalg.solve(np.eye(rhs.shape[0]) - dt*self.L2, rhs) + assert retval.shape == rhs.shape return retval - def evalF2(self, u: np.ndarray, t: float) -> np.ndarray: - """ - Solely compute tendencies of first frequency - """ + def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray: + from scipy.linalg import expm + assert isinstance(t, float) - retval = self.L2 @ u - assert retval.shape == u.shape + retval = expm(self.L * t) @ u0 + assert retval.shape == u0.shape return retval diff --git a/qmat/playgrounds/martin/ex_two_freq_integrate.py b/qmat/playgrounds/martin/ex_two_freq_integrate.py index 0d9fd14..a519a79 100644 --- a/qmat/playgrounds/martin/ex_two_freq_integrate.py +++ b/qmat/playgrounds/martin/ex_two_freq_integrate.py @@ -10,7 +10,8 @@ time_integration: str = "sdc" sdc_micro_time_integration: str = "irk1" -sdc_num_sweeps: int = 2 +sdc_micro_time_integration: str = "imex12" +sdc_num_sweeps: int = 4 two_freq: TwoFreq = TwoFreq(lam1=1.0j, lam2=20.0j, lam3=0.5j) u0 = two_freq.initial_u0() @@ -45,7 +46,11 @@ elif time_integration == "sdc": sdci = SDCIntegration( - num_nodes=5, node_type="LEGENDRE", quad_type="LOBATTO", num_sweeps=sdc_num_sweeps, micro_time_integration=sdc_micro_time_integration + 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): diff --git a/qmat/playgrounds/martin/tests/test_burgers.py b/qmat/playgrounds/martin/tests/test_burgers.py index de22a34..2ed80be 100644 --- a/qmat/playgrounds/martin/tests/test_burgers.py +++ b/qmat/playgrounds/martin/tests/test_burgers.py @@ -18,9 +18,9 @@ def test_burgers(): burgers.run_tests() - for time_integration in ["rk1", "rk2", "rk4", "sdc"]: + for time_integration_method in ["rk1", "rk2", "rk4", "sdc"]: print("="*80) - print(f"Time integration method: {time_integration}") + print(f"Time integration method: {time_integration_method}") print("="*80) results = [] @@ -38,11 +38,11 @@ def test_burgers(): u = u0.copy() - if time_integration in RKIntegration.supported_methods: - rki = RKIntegration(method=time_integration) + if time_integration_method in RKIntegration.supported_methods: + rki = RKIntegration(method=time_integration_method) u = rki.integrate_n(u, t, dt, num_timesteps, burgers) - elif time_integration == "sdc": + elif time_integration_method == "sdc": sdci = SDCIntegration(num_nodes=3, node_type="LEGENDRE", quad_type="LOBATTO") u = sdci.integrate_n(u, t, dt, num_timesteps, burgers) @@ -63,16 +63,19 @@ def test_burgers(): prev_error = r["error"] r["conv"] = conv - if time_integration == "rk1": + if time_integration_method == "rk1": assert results[-1]["error"] < 1e-4 assert np.abs(results[-1]["conv"]-1.0) < 1e-3 - elif time_integration == "rk2": + elif time_integration_method == "rk2": assert results[-1]["error"] < 1e-7 assert np.abs(results[-1]["conv"]-2.0) < 1e-3 - elif time_integration == "rk4": + elif time_integration_method == "rk4": assert results[-1]["error"] < 1e-14 - elif time_integration == "sdc": + elif time_integration_method == "sdc": assert results[-1]["error"] < 1e-14 + + else: + raise Exception(f"TODO for {time_integration_method}") diff --git a/qmat/playgrounds/martin/tests/test_dahlquist.py b/qmat/playgrounds/martin/tests/test_dahlquist.py index e4f7a3c..9f2d897 100644 --- a/qmat/playgrounds/martin/tests/test_dahlquist.py +++ b/qmat/playgrounds/martin/tests/test_dahlquist.py @@ -12,15 +12,15 @@ def test_dahlquist(): dahlquist: Dahlquist = Dahlquist(lam1=1.0j, lam2=0.1j) - for time_integration in ["rk1", "rk2", "rk4", "sdc"]: - if time_integration == "sdc": - micro_time_integration_ = ["erk1", "irk1"] + for time_integration_method in ["rk1", "rk2", "rk4", "sdc"]: + if time_integration_method == "sdc": + micro_time_integration_ = ["erk1", "irk1", "imex12", "imex21"] else: micro_time_integration_ = ["-"] for micro_time_integration in micro_time_integration_: print("=" * 80) - print(f"Time integration method: {time_integration} ({micro_time_integration})") + print(f"Time integration method: {time_integration_method} ({micro_time_integration})") print("=" * 80) results = [] @@ -35,11 +35,11 @@ def test_dahlquist(): u0 = dahlquist.initial_u0() u = u0.copy() - if time_integration in RKIntegration.supported_methods: - rki = RKIntegration(method=time_integration) + if time_integration_method in RKIntegration.supported_methods: + rki = RKIntegration(method=time_integration_method) u = rki.integrate_n(u, t, dt, num_timesteps, dahlquist) - elif time_integration == "sdc": + elif time_integration_method == "sdc": sdci = SDCIntegration( num_nodes=3, node_type="LEGENDRE", @@ -66,18 +66,21 @@ def test_dahlquist(): prev_error = r["error"] r["conv"] = conv - if time_integration == "rk1": + if time_integration_method == "rk1": assert results[-1]["error"] < 1e-2 assert np.abs(results[-1]["conv"] - 1.0) < 1e-2 - elif time_integration == "rk2": + elif time_integration_method == "rk2": assert results[-1]["error"] < 1e-4 assert np.abs(results[-1]["conv"] - 2.0) < 1e-3 - elif time_integration == "rk4": + elif time_integration_method == "rk4": assert results[-1]["error"] < 1e-9 assert np.abs(results[-1]["conv"] - 4.0) < 1e-4 - elif time_integration == "sdc": + elif time_integration_method == "sdc": assert results[-1]["error"] < 1e-8 assert np.abs(results[-1]["conv"] - 4.0) < 1e-3 + + else: + raise Exception(f"TODO for {time_integration_method}") diff --git a/qmat/playgrounds/martin/tests/test_dahlquist2.py b/qmat/playgrounds/martin/tests/test_dahlquist2.py index f022ba7..b534131 100644 --- a/qmat/playgrounds/martin/tests/test_dahlquist2.py +++ b/qmat/playgrounds/martin/tests/test_dahlquist2.py @@ -12,9 +12,9 @@ def test_dahlquist2(): dahlquist2: Dahlquist2 = Dahlquist2(lam1=1.0j, lam2=0.1j) - for time_integration in ["rk1", "rk2", "rk4", "sdc"]: + for time_integration_method in ["rk1", "rk2", "rk4", "sdc"]: print("="*80) - print(f"Time integration method: {time_integration}") + print(f"Time integration method: {time_integration_method}") print("="*80) results = [] @@ -31,12 +31,12 @@ def test_dahlquist2(): u = u0.copy() - if time_integration in RKIntegration.supported_methods: - rki = RKIntegration(method=time_integration) + if time_integration_method in RKIntegration.supported_methods: + rki = RKIntegration(method=time_integration_method) u = rki.integrate_n(u, t, dt, num_timesteps, dahlquist2) - elif time_integration == "sdc": + elif time_integration_method == "sdc": sdci = SDCIntegration(num_nodes=3, node_type="LEGENDRE", quad_type="LOBATTO", num_sweeps=1) u = sdci.integrate_n(u, t, dt, num_timesteps, dahlquist2) @@ -57,18 +57,21 @@ def test_dahlquist2(): prev_error = r["error"] r["conv"] = conv - if time_integration == "rk1": + if time_integration_method == "rk1": assert results[-1]["error"] < 1e-2 assert np.abs(results[-1]["conv"]-1.0) < 1e-2 - elif time_integration == "rk2": + elif time_integration_method == "rk2": assert results[-1]["error"] < 1e-5 assert np.abs(results[-1]["conv"]-2.0) < 1e-3 - elif time_integration == "rk4": + elif time_integration_method == "rk4": assert results[-1]["error"] < 1e-11 assert np.abs(results[-1]["conv"]-4.0) < 1e-2 - elif time_integration == "sdc": + elif time_integration_method == "sdc": assert results[-1]["error"] < 1e-5 assert np.abs(results[-1]["conv"]-2.0) < 1e-3 + + else: + raise Exception(f"TODO for {time_integration_method}") diff --git a/qmat/playgrounds/martin/tests/test_dahlquist_ext_freq.py b/qmat/playgrounds/martin/tests/test_dahlquist_ext_freq.py index 93ee8a7..81174a2 100644 --- a/qmat/playgrounds/martin/tests/test_dahlquist_ext_freq.py +++ b/qmat/playgrounds/martin/tests/test_dahlquist_ext_freq.py @@ -12,15 +12,15 @@ def test_dahlquist(): dahlquist: Dahlquist = Dahlquist(lam1=1.0j, lam2=0.1j, ext_scalar=1.0) - for time_integration in ["rk1", "rk2", "rk4", "irk1", "irk2", "sdc"]: - if time_integration == "sdc": + for time_integration_method in ["rk1", "rk2", "rk4", "irk1", "irk2", "sdc"]: + if time_integration_method == "sdc": micro_time_integration_ = ["erk1", "irk1"] else: micro_time_integration_ = ["-"] for micro_time_integration in micro_time_integration_: print("=" * 80) - print(f"Time integration method: {time_integration} ({micro_time_integration})") + print(f"Time integration method: {time_integration_method} ({micro_time_integration})") print("=" * 80) results = [] @@ -37,12 +37,12 @@ def test_dahlquist(): u = u0.copy() - if time_integration in RKIntegration.supported_methods: - rki = RKIntegration(method=time_integration) + if time_integration_method in RKIntegration.supported_methods: + rki = RKIntegration(method=time_integration_method) u = rki.integrate_n(u, t, dt, num_timesteps, dahlquist) - elif time_integration == "sdc": + elif time_integration_method == "sdc": sdci = SDCIntegration( num_nodes=3, node_type="LEGENDRE", @@ -69,26 +69,29 @@ def test_dahlquist(): prev_error = r["error"] r["conv"] = conv - if time_integration == "rk1": + if time_integration_method == "rk1": assert results[-1]["error"] < 1e-2 assert np.abs(results[-1]["conv"] - 1.0) < 1e-2 - elif time_integration == "rk2": + elif time_integration_method == "rk2": assert results[-1]["error"] < 1e-4 assert np.abs(results[-1]["conv"] - 2.0) < 1e-2 - elif time_integration == "rk4": + elif time_integration_method == "rk4": assert results[-1]["error"] < 1e-9 assert np.abs(results[-1]["conv"] - 4.0) < 1e-2 - elif time_integration == "irk1": + elif time_integration_method == "irk1": assert results[-1]["error"] < 1e-2 assert np.abs(results[-1]["conv"] - 1.0) < 1e-2 - elif time_integration == "irk2": + elif time_integration_method == "irk2": assert results[-1]["error"] < 1e-2 assert np.abs(results[-1]["conv"] - 2.0) < 1e-2 - elif time_integration == "sdc": + elif time_integration_method == "sdc": assert results[-1]["error"] < 1e-8 assert np.abs(results[-1]["conv"] - 4.0) < 1e-1 + + else: + raise Exception(f"TODO for {time_integration_method}") diff --git a/qmat/playgrounds/martin/tests/test_two_freq.py b/qmat/playgrounds/martin/tests/test_two_freq.py index c390ef5..d080479 100644 --- a/qmat/playgrounds/martin/tests/test_two_freq.py +++ b/qmat/playgrounds/martin/tests/test_two_freq.py @@ -12,15 +12,15 @@ def test_dahlquist2(): two_freq: TwoFreq = TwoFreq(lam1=1.0j, lam2=0.1j) u0 = two_freq.initial_u0() - for time_integration in ["rk1", "rk2", "rk4", "irk1", "irk2", "sdc"]: - if time_integration == "sdc": - micro_time_integration_ = ["erk1", "irk1"] + for time_integration_method in ["rk1", "rk2", "rk4", "irk1", "irk2", "sdc"]: + if time_integration_method == "sdc": + micro_time_integration_ = ["erk1", "irk1", "imex12", "imex21"] else: micro_time_integration_ = ["-"] for micro_time_integration in micro_time_integration_: print("=" * 80) - print(f"Time integration method: {time_integration} ({micro_time_integration})") + print(f"Time integration method: {time_integration_method} ({micro_time_integration})") print("=" * 80) results = [] @@ -37,12 +37,12 @@ def test_dahlquist2(): u = u0.copy() - if time_integration in RKIntegration.supported_methods: - rki = RKIntegration(method=time_integration) + if time_integration_method in RKIntegration.supported_methods: + rki = RKIntegration(method=time_integration_method) u = rki.integrate_n(u, t, dt, num_timesteps, two_freq) - elif time_integration == "sdc": + elif time_integration_method == "sdc": sdci = SDCIntegration( num_nodes=3, node_type="LEGENDRE", @@ -69,26 +69,29 @@ def test_dahlquist2(): prev_error = r["error"] r["conv"] = conv - if time_integration == "rk1": + if time_integration_method == "rk1": assert results[-1]["error"] < 1e-2 assert np.abs(results[-1]["conv"] - 1.0) < 1e-2 - elif time_integration == "rk2": + elif time_integration_method == "rk2": assert results[-1]["error"] < 1e-5 assert np.abs(results[-1]["conv"] - 2.0) < 1e-3 - elif time_integration == "rk4": + elif time_integration_method == "rk4": assert results[-1]["error"] < 1e-11 assert np.abs(results[-1]["conv"] - 4.0) < 1e-2 - elif time_integration == "irk1": + elif time_integration_method == "irk1": assert results[-1]["error"] < 1e-2 assert np.abs(results[-1]["conv"] - 1.0) < 1e-2 - elif time_integration == "irk2": + elif time_integration_method == "irk2": assert results[-1]["error"] < 1e-2 assert np.abs(results[-1]["conv"] - 2.0) < 1e-2 - elif time_integration == "sdc": + elif time_integration_method == "sdc": assert results[-1]["error"] < 1e-8 assert np.abs(results[-1]["conv"] - 4.0) < 1e-3 + + else: + raise Exception(f"TODO for {time_integration_method}") diff --git a/qmat/playgrounds/martin/time_integration/sdc_integration.py b/qmat/playgrounds/martin/time_integration/sdc_integration.py index 210b233..feb03df 100644 --- a/qmat/playgrounds/martin/time_integration/sdc_integration.py +++ b/qmat/playgrounds/martin/time_integration/sdc_integration.py @@ -179,11 +179,199 @@ def integrate_irk1(self, u0: np.array, t: float, dt: float, de_solver: DESolver) assert u0.shape == u[0].shape return u[0] + def integrate_imex21(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> np.array: + """ + IMEX SDC where the first term is treated implicitly and the second term explicitly. + """ + if not np.isclose(self.nodes[0], 0.0): + raise Exception("SDC nodes must include the left interval boundary.") + + if not np.isclose(self.nodes[-1], 1.0): + raise Exception("SDC nodes must include the right interval boundary.") + + assert self.N == self.deltas.shape[0] + + shape = (self.N,) + u0.shape + u = np.zeros_like(u0, shape=shape) + + u[0, :] = u0 + evalF_k0 = np.empty_like(u) + evalF_k1 = np.empty_like(u) + + # Backup integrator contribution I[...] of previous iteration + ISolves = np.empty_like(u) + + # + # Propagate initial condition to all nodes + # + for m in range(0, self.N): + if m > 0: + # + # Solve explicit step first + # u* = u^n + dt * delta * F1(u^n) + # + rhs = u[m - 1] + dt * self.deltas[m] * de_solver.evalF1(u[m - 1], t + dt * self.nodes[m]) + + # + # Solve implicit step next (see integrate_irk1) + # u^{n+1} = u* + dt * delta * F2(u^{n+1}) + # + u[m] = de_solver.fSolve2(rhs, dt * self.deltas[m], t + dt * self.nodes[m]) + + # + # Compute I[...] term for implicit and explicit parts + # + # u^n+1 = u^n + dt * delta * F1(u^n) + dt * delta * F2(u^n+1) + # dt * delta * F1(u^n) + dt * delta * F2(u^n+1) = u^n+1 - u^n + # + ISolves[m] = u[m] - u[m - 1] + + evalF_k0[m] = de_solver.evalF(u[m], t + dt * self.nodes[m]) + + # + # Iteratively sweep over SDC nodes + # + for _ in range(1, self.num_sweeps): + for m in range(0, self.N): + if m > 0: + # + # Solve IMEX step: + # + # u^n+1 = u^n + dt * delta * (F1(u^n) + F2(u^n+1)) - I(u^n) + dt * Q * F(u^n) + # <=> u^n+1 - dt * delta * F2(u^n+1) = u^n + dt * delta * F1(u^n) - I(u^n) + dt * Q * F(u^n) + # <=> (I - dt * delta * F) * u^n+1 = u^n + dt * delta * F1(u^n) - I(u^n) + dt * Q * F(u^n) + # + # rhs = u^n + dt * delta * F1(u^n) - I(u^n) + dt * Q * F(u^n) + # + qeval = self.q[m] @ evalF_k0 + euler = dt * self.deltas[m] * de_solver.evalF1(u[m - 1], t + dt * self.nodes[m]) + rhs = u[m - 1] + euler - ISolves[m] + dt * qeval + + u[m] = de_solver.fSolve2(rhs, dt * self.deltas[m], t + dt * self.nodes[m]) + + # + # Update I[...] term for next correction + # <=> dt * delta * (F1(u^n) + F2(u^n+1)) = u^n+1 - u^n + I(u^n) - dt * Q * F(u^n) + # + ISolves[m] = u[m] - u[m - 1] + ISolves[m] - dt * qeval + + evalF_k1[m] = de_solver.evalF(u[m], t + dt * self.nodes[m]) + + # Copy tendency arrays + evalF_k0[:] = evalF_k1[:] + + if 0: + # If we're using Radau-right, we can just use the last value + u[0] = u[-1] + + else: + # Compute new starting value with quadrature on tendencies + u[0] = u[0] + dt * self.weights @ evalF_k0 + + assert u0.shape == u[0].shape + return u[0] + + def integrate_imex12(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> np.array: + """ + IMEX SDC where the first term is treated implicitly and the second term explicitly. + """ + if not np.isclose(self.nodes[0], 0.0): + raise Exception("SDC nodes must include the left interval boundary.") + + if not np.isclose(self.nodes[-1], 1.0): + raise Exception("SDC nodes must include the right interval boundary.") + + assert self.N == self.deltas.shape[0] + + shape = (self.N,) + u0.shape + u = np.zeros_like(u0, shape=shape) + + u[0, :] = u0 + evalF_k0 = np.empty_like(u) + evalF_k1 = np.empty_like(u) + + # Backup integrator contribution I[...] of previous iteration + ISolves = np.empty_like(u) + + # + # Propagate initial condition to all nodes + # + for m in range(0, self.N): + if m > 0: + # + # Solve explicit step first + # u* = u^n + dt * delta * F1(u^n) + # + rhs = u[m - 1] + dt * self.deltas[m] * de_solver.evalF2(u[m - 1], t + dt * self.nodes[m]) + + # + # Solve implicit step next (see integrate_irk1) + # u^{n+1} = u* + dt * delta * F2(u^{n+1}) + # + u[m] = de_solver.fSolve1(rhs, dt * self.deltas[m], t + dt * self.nodes[m]) + + # + # Compute I[...] term for implicit and explicit parts + # + # u^n+1 = u^n + dt * delta * F1(u^n) + dt * delta * F2(u^n+1) + # dt * delta * F1(u^n) + dt * delta * F2(u^n+1) = u^n+1 - u^n + # + ISolves[m] = u[m] - u[m - 1] + + evalF_k0[m] = de_solver.evalF(u[m], t + dt * self.nodes[m]) + + # + # Iteratively sweep over SDC nodes + # + for _ in range(1, self.num_sweeps): + for m in range(0, self.N): + if m > 0: + # + # Solve IMEX step: + # + # u^n+1 = u^n + dt * delta * (F1(u^n) + F2(u^n+1)) - I(u^n) + dt * Q * F(u^n) + # <=> u^n+1 - dt * delta * F2(u^n+1) = u^n + dt * delta * F1(u^n) - I(u^n) + dt * Q * F(u^n) + # <=> (I - dt * delta * F) * u^n+1 = u^n + dt * delta * F1(u^n) - I(u^n) + dt * Q * F(u^n) + # + # rhs = u^n + dt * delta * F1(u^n) - I(u^n) + dt * Q * F(u^n) + # + qeval = self.q[m] @ evalF_k0 + euler = dt * self.deltas[m] * de_solver.evalF2(u[m - 1], t + dt * self.nodes[m]) + rhs = u[m - 1] + euler - ISolves[m] + dt * qeval + + u[m] = de_solver.fSolve1(rhs, dt * self.deltas[m], t + dt * self.nodes[m]) + + # + # Update I[...] term for next correction + # <=> dt * delta * (F1(u^n) + F2(u^n+1)) = u^n+1 - u^n + I(u^n) - dt * Q * F(u^n) + # + ISolves[m] = u[m] - u[m - 1] + ISolves[m] - dt * qeval + + evalF_k1[m] = de_solver.evalF(u[m], t + dt * self.nodes[m]) + + # Copy tendency arrays + evalF_k0[:] = evalF_k1[:] + + if 0: + # If we're using Radau-right, we can just use the last value + u[0] = u[-1] + + else: + # Compute new starting value with quadrature on tendencies + u[0] = u[0] + dt * self.weights @ evalF_k0 + + assert u0.shape == u[0].shape + return u[0] + def integrate(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> np.array: if self.time_integration_method == "erk1": return self.integrate_erk1(u0, t, dt, de_solver) elif self.time_integration_method == "irk1": return self.integrate_irk1(u0, t, dt, de_solver) + elif self.time_integration_method == "imex12": + return self.integrate_imex12(u0, t, dt, de_solver) + elif self.time_integration_method == "imex21": + return self.integrate_imex21(u0, t, dt, de_solver) else: raise Exception(f"Unsupported time integration within SDC: '{self.time_integration_method}'") From 15abdac9e1ae5cf837a371f64cba8c6db0898127 Mon Sep 17 00:00:00 2001 From: SCHREIBER Martin Date: Sat, 25 Oct 2025 17:51:48 +0200 Subject: [PATCH 10/11] dahlquist arbitrary t0 integrator --- qmat/playgrounds/martin/diff_eqs/burgers.py | 2 +- qmat/playgrounds/martin/diff_eqs/dahlquist.py | 40 +++++++++++++++++-- .../playgrounds/martin/diff_eqs/dahlquist2.py | 2 +- qmat/playgrounds/martin/diff_eqs/de_solver.py | 2 +- qmat/playgrounds/martin/diff_eqs/two_freq.py | 2 +- .../martin/ex_burgers_integrate.py | 2 +- qmat/playgrounds/martin/ex_burgers_plot.py | 2 +- .../martin/ex_dahlquist2_integrate.py | 8 ++-- qmat/playgrounds/martin/ex_dahlquist2_plot.py | 2 +- qmat/playgrounds/martin/ex_dahlquist_plot.py | 2 +- .../martin/ex_two_freq_integrate.py | 8 ++-- qmat/playgrounds/martin/ex_two_freq_plot.py | 2 +- qmat/playgrounds/martin/tests/test_burgers.py | 2 +- .../martin/tests/test_dahlquist.py | 20 +++++++++- .../martin/tests/test_dahlquist2.py | 2 +- .../martin/tests/test_dahlquist_ext_freq.py | 2 +- .../playgrounds/martin/tests/test_two_freq.py | 2 +- .../martin/time_integration/rk_integration.py | 4 +- .../time_integration/sdc_integration.py | 4 +- 19 files changed, 80 insertions(+), 30 deletions(-) diff --git a/qmat/playgrounds/martin/diff_eqs/burgers.py b/qmat/playgrounds/martin/diff_eqs/burgers.py index 786eda2..3109b5e 100644 --- a/qmat/playgrounds/martin/diff_eqs/burgers.py +++ b/qmat/playgrounds/martin/diff_eqs/burgers.py @@ -83,7 +83,7 @@ def evalF(self, u: np.ndarray, t: float) -> np.ndarray: f = -u * du_dx + self._nu * d2u_dx2 return f - def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray: + def int_f(self, u0: np.ndarray, t: float) -> np.ndarray: """ Compute the analytical solution of the 1D viscous Burgers' equation at time `t`. diff --git a/qmat/playgrounds/martin/diff_eqs/dahlquist.py b/qmat/playgrounds/martin/diff_eqs/dahlquist.py index e18bffc..2ef8283 100644 --- a/qmat/playgrounds/martin/diff_eqs/dahlquist.py +++ b/qmat/playgrounds/martin/diff_eqs/dahlquist.py @@ -26,6 +26,19 @@ def __init__(self, lam1: complex, lam2: complex, ext_scalar: float = 0.0): 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 @@ -77,11 +90,32 @@ def fSolve2(self, u: np.ndarray, dt: float, t: float) -> np.ndarray: assert retval.shape == u.shape return retval - def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray: + def int_f_t0(self, u0: np.ndarray, dt: float) -> np.ndarray: lam = self.lam1 + self.lam2 s = self.ext_scalar - assert isinstance(t, float) - retval = np.exp(t * lam) * u0 + s * np.sin(t) + 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, t0: float = 0.0) -> np.ndarray: + """ + Integrate the solution from t0 to t1. + """ + + assert isinstance(t0, (float, int)) + assert isinstance(dt, (float, int)) + + if t0 == 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(t0)) + s * np.sin(t0+dt) assert retval.shape == u0.shape return retval diff --git a/qmat/playgrounds/martin/diff_eqs/dahlquist2.py b/qmat/playgrounds/martin/diff_eqs/dahlquist2.py index 643d900..30bdbb1 100644 --- a/qmat/playgrounds/martin/diff_eqs/dahlquist2.py +++ b/qmat/playgrounds/martin/diff_eqs/dahlquist2.py @@ -36,7 +36,7 @@ def evalF(self, u: np.ndarray, t: float) -> np.ndarray: assert retval.shape == u.shape return retval - def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray: + def int_f(self, u0: np.ndarray, t: float) -> np.ndarray: assert isinstance(t, float) retval = (self.s * np.exp(t * self.lam1) + (1.0 - self.s) * np.exp(t * self.lam2)) * u0 assert retval.shape == u0.shape diff --git a/qmat/playgrounds/martin/diff_eqs/de_solver.py b/qmat/playgrounds/martin/diff_eqs/de_solver.py index ceed5f0..f01ca31 100644 --- a/qmat/playgrounds/martin/diff_eqs/de_solver.py +++ b/qmat/playgrounds/martin/diff_eqs/de_solver.py @@ -117,7 +117,7 @@ def initial_u0(self, mode: str) -> np.ndarray: pass @abstractmethod - def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray: + def int_f(self, u0: np.ndarray, t: float) -> np.ndarray: """ Compute the (analytical) solution at time `t`. diff --git a/qmat/playgrounds/martin/diff_eqs/two_freq.py b/qmat/playgrounds/martin/diff_eqs/two_freq.py index 88c0ac5..0a9ea6a 100644 --- a/qmat/playgrounds/martin/diff_eqs/two_freq.py +++ b/qmat/playgrounds/martin/diff_eqs/two_freq.py @@ -90,7 +90,7 @@ def fSolve2(self, rhs: np.ndarray, dt: float, t: float) -> np.ndarray: assert retval.shape == rhs.shape return retval - def u_solution(self, u0: np.ndarray, t: float) -> np.ndarray: + def int_f(self, u0: np.ndarray, t: float) -> np.ndarray: from scipy.linalg import expm assert isinstance(t, float) diff --git a/qmat/playgrounds/martin/ex_burgers_integrate.py b/qmat/playgrounds/martin/ex_burgers_integrate.py index d92d141..6df3ee9 100644 --- a/qmat/playgrounds/martin/ex_burgers_integrate.py +++ b/qmat/playgrounds/martin/ex_burgers_integrate.py @@ -30,7 +30,7 @@ if 1: results = [] - u_analytical = burgers.u_solution(u0, t=T) + u_analytical = burgers.int_f(u0, t=T) for nt in range(4): diff --git a/qmat/playgrounds/martin/ex_burgers_plot.py b/qmat/playgrounds/martin/ex_burgers_plot.py index 2752209..e15a6e9 100644 --- a/qmat/playgrounds/martin/ex_burgers_plot.py +++ b/qmat/playgrounds/martin/ex_burgers_plot.py @@ -32,7 +32,7 @@ plt.plot(t, ut) else: - ut = burgers.u_solution(u0, t=T) + ut = burgers.int_f(u0, t=T) plt.plot(t, u0, label="u0") plt.plot(t, ut, label="ut") diff --git a/qmat/playgrounds/martin/ex_dahlquist2_integrate.py b/qmat/playgrounds/martin/ex_dahlquist2_integrate.py index f453b0a..db0449f 100644 --- a/qmat/playgrounds/martin/ex_dahlquist2_integrate.py +++ b/qmat/playgrounds/martin/ex_dahlquist2_integrate.py @@ -25,7 +25,7 @@ dt = T / num_timesteps - u_analytical_fin = dahlquist2.u_solution(u0, t=T) + u_analytical_fin = dahlquist2.int_f(u0, t=T) u0 = dahlquist2.initial_u0() u = u0.copy() @@ -37,7 +37,7 @@ rki = RKIntegration(method=time_integration) for n in range(num_timesteps): - u = rki.integrate(u, t + n * dt, dt, dahlquist2) + u = rki.int_f(u, t + n * dt, dt, dahlquist2) u_ = np.concatenate((u_, np.expand_dims(u, axis=0))) t_ = np.append(t_, t + (n + 1) * dt) @@ -46,7 +46,7 @@ sdci = SDCIntegration(num_nodes=3, node_type="LEGENDRE", quad_type="LOBATTO") for n in range(num_timesteps): - u = sdci.integrate(u, t + n * dt, dt, dahlquist2) + u = sdci.int_f(u, t + n * dt, dt, dahlquist2) u_ = np.concatenate((u_, np.expand_dims(u, axis=0))) t_ = np.append(t_, t + (n + 1) * dt) @@ -72,7 +72,7 @@ if 1: t_ = np.linspace(0, T, 1000) - u_analytical_fin = np.array([dahlquist2.u_solution(u0, t) for t in t_]) + u_analytical_fin = np.array([dahlquist2.int_f(u0, t) for t in t_]) plt.plot(t_, np.real(u_analytical_fin), linestyle="dotted", color="gray") plt.plot(t_, np.imag(u_analytical_fin), linestyle="dotted", color="gray") diff --git a/qmat/playgrounds/martin/ex_dahlquist2_plot.py b/qmat/playgrounds/martin/ex_dahlquist2_plot.py index 168c56a..6da0e78 100644 --- a/qmat/playgrounds/martin/ex_dahlquist2_plot.py +++ b/qmat/playgrounds/martin/ex_dahlquist2_plot.py @@ -9,7 +9,7 @@ dahlquist2: Dahlquist2 = Dahlquist2(lam1=20.0j, lam2=1.0j, s=0.1) -u_eval = np.array([dahlquist2.u_solution(u0, _) for _ in t]) +u_eval = np.array([dahlquist2.int_f(u0, _) for _ in t]) plt.plot(t, np.real(u_eval), label="Re(u)") diff --git a/qmat/playgrounds/martin/ex_dahlquist_plot.py b/qmat/playgrounds/martin/ex_dahlquist_plot.py index 232a22b..bd6dacc 100644 --- a/qmat/playgrounds/martin/ex_dahlquist_plot.py +++ b/qmat/playgrounds/martin/ex_dahlquist_plot.py @@ -9,7 +9,7 @@ dahlquist: Dahlquist = Dahlquist(lam1=20.0j, lam2=1.0j) -u_eval = np.array([dahlquist.u_solution(u0, _) for _ in t]) +u_eval = np.array([dahlquist.int_f(u0, _) for _ in t]) plt.plot(t, np.real(u_eval), label="Re(u)") diff --git a/qmat/playgrounds/martin/ex_two_freq_integrate.py b/qmat/playgrounds/martin/ex_two_freq_integrate.py index a519a79..00fb13a 100644 --- a/qmat/playgrounds/martin/ex_two_freq_integrate.py +++ b/qmat/playgrounds/martin/ex_two_freq_integrate.py @@ -27,7 +27,7 @@ dt = T / num_timesteps - u_analytical_fin = two_freq.u_solution(u0, t=T) + u_analytical_fin = two_freq.int_f(u0, t=T) u0 = two_freq.initial_u0() u = u0.copy() @@ -39,7 +39,7 @@ rki = RKIntegration(method=time_integration) for n in range(num_timesteps): - u = rki.integrate(u, t + n * dt, dt, two_freq) + u = rki.int_f(u, t + n * dt, dt, two_freq) u_ = np.concatenate((u_, np.expand_dims(u, axis=0))) t_ = np.append(t_, t + (n + 1) * dt) @@ -54,7 +54,7 @@ ) for n in range(num_timesteps): - u = sdci.integrate(u, t + n * dt, dt, two_freq) + u = sdci.int_f(u, t + n * dt, dt, two_freq) u_ = np.concatenate((u_, np.expand_dims(u, axis=0))) t_ = np.append(t_, t + (n + 1) * dt) @@ -73,7 +73,7 @@ if 1: # Plot analytical solution t_ = np.linspace(0, T, 1000) - u_analytical_fin = np.array([two_freq.u_solution(u0, t) for t in t_]) + u_analytical_fin = np.array([two_freq.int_f(u0, t) for t in t_]) for i in range(u_analytical_fin.shape[1]): plt.plot(t_, np.real(u_analytical_fin[:, i]), linestyle="dotted", color="black") plt.plot(t_, np.imag(u_analytical_fin[:, i]), linestyle="dotted", color="black") diff --git a/qmat/playgrounds/martin/ex_two_freq_plot.py b/qmat/playgrounds/martin/ex_two_freq_plot.py index 54df8ce..117b29a 100644 --- a/qmat/playgrounds/martin/ex_two_freq_plot.py +++ b/qmat/playgrounds/martin/ex_two_freq_plot.py @@ -9,7 +9,7 @@ two_freq: TwoFreq = TwoFreq(lam1=1.0j, lam2=20.0j, lam3=0.5j) u0 = two_freq.initial_u0() -u_eval = np.array([two_freq.u_solution(u0, _) for _ in t]) +u_eval = np.array([two_freq.int_f(u0, _) for _ in t]) for i in range(2): plt.plot(t, np.real(u_eval[:, i]), label=f"Re(u[{i}])") diff --git a/qmat/playgrounds/martin/tests/test_burgers.py b/qmat/playgrounds/martin/tests/test_burgers.py index 2ed80be..931a004 100644 --- a/qmat/playgrounds/martin/tests/test_burgers.py +++ b/qmat/playgrounds/martin/tests/test_burgers.py @@ -24,7 +24,7 @@ def test_burgers(): print("="*80) results = [] - u_analytical = burgers.u_solution(u0, t=T) + u_analytical = burgers.int_f(u0, t=T) for nt in range(2, 4): diff --git a/qmat/playgrounds/martin/tests/test_dahlquist.py b/qmat/playgrounds/martin/tests/test_dahlquist.py index 9f2d897..69612b1 100644 --- a/qmat/playgrounds/martin/tests/test_dahlquist.py +++ b/qmat/playgrounds/martin/tests/test_dahlquist.py @@ -4,7 +4,7 @@ from qmat.playgrounds.martin.time_integration.rk_integration import RKIntegration -def test_dahlquist(): +def test_dahlquist_integration(): u0 = np.array([1.0]) # Initial condition T: float = 4 * np.pi # Time interval T: float = 0.5 # Time interval @@ -24,7 +24,7 @@ def test_dahlquist(): print("=" * 80) results = [] - u_analytical = dahlquist.u_solution(u0, t=T) + u_analytical = dahlquist.int_f(u0, dt=T) for nt in range(4): @@ -84,3 +84,19 @@ def test_dahlquist(): else: raise Exception(f"TODO for {time_integration_method}") + + +def test_dahlquist_integration_of_solution(): + u0 = np.array([1.0]) # Initial condition + T: float = 4 * np.pi # Time interval + + dahlquist: Dahlquist = Dahlquist(lam1=1.0j, lam2=0.1j) + + u0 = dahlquist.initial_u0() + + u_analytical = dahlquist.int_f_t0(u0, dt=T) + + u1 = dahlquist.int_f(u0, dt=T*0.5, t0=0.0) + u2 = dahlquist.int_f(u1, dt=T*0.5, t0=T*0.5) + + assert np.all(np.isclose(u_analytical, u2)) diff --git a/qmat/playgrounds/martin/tests/test_dahlquist2.py b/qmat/playgrounds/martin/tests/test_dahlquist2.py index b534131..efcb179 100644 --- a/qmat/playgrounds/martin/tests/test_dahlquist2.py +++ b/qmat/playgrounds/martin/tests/test_dahlquist2.py @@ -18,7 +18,7 @@ def test_dahlquist2(): print("="*80) results = [] - u_analytical = dahlquist2.u_solution(u0, t=T) + u_analytical = dahlquist2.int_f(u0, t=T) for nt in range(4): diff --git a/qmat/playgrounds/martin/tests/test_dahlquist_ext_freq.py b/qmat/playgrounds/martin/tests/test_dahlquist_ext_freq.py index 81174a2..403e19a 100644 --- a/qmat/playgrounds/martin/tests/test_dahlquist_ext_freq.py +++ b/qmat/playgrounds/martin/tests/test_dahlquist_ext_freq.py @@ -24,7 +24,7 @@ def test_dahlquist(): print("=" * 80) results = [] - u_analytical = dahlquist.u_solution(u0, t=T) + u_analytical = dahlquist.int_f(u0, dt=T) for nt in range(4): diff --git a/qmat/playgrounds/martin/tests/test_two_freq.py b/qmat/playgrounds/martin/tests/test_two_freq.py index d080479..2aa3a30 100644 --- a/qmat/playgrounds/martin/tests/test_two_freq.py +++ b/qmat/playgrounds/martin/tests/test_two_freq.py @@ -24,7 +24,7 @@ def test_dahlquist2(): print("=" * 80) results = [] - u_analytical = two_freq.u_solution(u0, t=T) + u_analytical = two_freq.int_f(u0, t=T) for nt in range(4): diff --git a/qmat/playgrounds/martin/time_integration/rk_integration.py b/qmat/playgrounds/martin/time_integration/rk_integration.py index 5ffde40..1cb2fd2 100644 --- a/qmat/playgrounds/martin/time_integration/rk_integration.py +++ b/qmat/playgrounds/martin/time_integration/rk_integration.py @@ -12,7 +12,7 @@ def __init__(self, method: str): assert method in self.supported_methods, "Unsupported RK method" self.method = method - def integrate(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> np.array: + def int_f(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> np.array: u = u0 if self.method == "rk1": @@ -56,6 +56,6 @@ def integrate_n(self, u0: np.array, t: float, dt: float, num_timesteps, de_solve u_value = u0 for n in range(num_timesteps): - u_value = self.integrate(u_value, t + n * dt, dt, de_solver) + u_value = self.int_f(u_value, t + n * dt, dt, de_solver) return u_value diff --git a/qmat/playgrounds/martin/time_integration/sdc_integration.py b/qmat/playgrounds/martin/time_integration/sdc_integration.py index feb03df..7a6bb09 100644 --- a/qmat/playgrounds/martin/time_integration/sdc_integration.py +++ b/qmat/playgrounds/martin/time_integration/sdc_integration.py @@ -363,7 +363,7 @@ def integrate_imex12(self, u0: np.array, t: float, dt: float, de_solver: DESolve assert u0.shape == u[0].shape return u[0] - def integrate(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> np.array: + def int_f(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> np.array: if self.time_integration_method == "erk1": return self.integrate_erk1(u0, t, dt, de_solver) elif self.time_integration_method == "irk1": @@ -385,6 +385,6 @@ def integrate_n(self, u0: np.array, t: float, dt: float, num_timesteps, de_solve u_value = u0 for n in range(num_timesteps): - u_value = self.integrate(u_value, t + n * dt, dt, de_solver) + u_value = self.int_f(u_value, t + n * dt, dt, de_solver) return u_value From fa00525d005663399eade2f94bd058ca75a443a5 Mon Sep 17 00:00:00 2001 From: SCHREIBER Martin Date: Wed, 5 Nov 2025 10:40:41 +0100 Subject: [PATCH 11/11] u --- qmat/playgrounds/martin/diff_eqs/burgers.py | 4 +- qmat/playgrounds/martin/diff_eqs/dahlquist.py | 8 +- .../playgrounds/martin/diff_eqs/dahlquist2.py | 7 +- qmat/playgrounds/martin/diff_eqs/de_solver.py | 8 +- qmat/playgrounds/martin/diff_eqs/two_freq.py | 6 +- .../martin/ex_burgers_integrate.py | 2 +- qmat/playgrounds/martin/ex_burgers_plot.py | 2 +- .../martin/ex_dahlquist2_integrate.py | 2 +- qmat/playgrounds/martin/ex_dahlquist_plot.py | 2 +- .../martin/ex_two_freq_integrate.py | 4 +- qmat/playgrounds/martin/ex_two_freq_plot.py | 2 +- .../martin/nyquist_test/ex_slow_fast.py | 104 ++++++++++++++++++ qmat/playgrounds/martin/tests/test_burgers.py | 2 +- .../martin/tests/test_dahlquist.py | 4 +- .../martin/tests/test_dahlquist2.py | 2 +- .../playgrounds/martin/tests/test_two_freq.py | 2 +- .../time_integration/sdc_integration.py | 48 ++++---- 17 files changed, 162 insertions(+), 47 deletions(-) create mode 100644 qmat/playgrounds/martin/nyquist_test/ex_slow_fast.py diff --git a/qmat/playgrounds/martin/diff_eqs/burgers.py b/qmat/playgrounds/martin/diff_eqs/burgers.py index 3109b5e..f1771c4 100644 --- a/qmat/playgrounds/martin/diff_eqs/burgers.py +++ b/qmat/playgrounds/martin/diff_eqs/burgers.py @@ -83,7 +83,7 @@ def evalF(self, u: np.ndarray, t: float) -> np.ndarray: f = -u * du_dx + self._nu * d2u_dx2 return f - def int_f(self, u0: np.ndarray, t: float) -> np.ndarray: + 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`. @@ -121,7 +121,7 @@ def int_f(self, u0: np.ndarray, t: float) -> np.ndarray: 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) + phi_hat = phi_hat * np.exp(self._nu * self._d_dx_**2 * (t+dt)) phi = np.fft.ifft(phi_hat) phi = np.log(phi) diff --git a/qmat/playgrounds/martin/diff_eqs/dahlquist.py b/qmat/playgrounds/martin/diff_eqs/dahlquist.py index 2ef8283..8fc7447 100644 --- a/qmat/playgrounds/martin/diff_eqs/dahlquist.py +++ b/qmat/playgrounds/martin/diff_eqs/dahlquist.py @@ -99,15 +99,15 @@ def int_f_t0(self, u0: np.ndarray, dt: float) -> np.ndarray: assert retval.shape == u0.shape return retval - def int_f(self, u0: np.ndarray, dt: float, t0: float = 0.0) -> np.ndarray: + def int_f(self, u0: np.ndarray, dt: float, t: float = 0.0) -> np.ndarray: """ Integrate the solution from t0 to t1. """ - assert isinstance(t0, (float, int)) + assert isinstance(t, (float, int)) assert isinstance(dt, (float, int)) - if t0 == 0: + if t == 0: return self.int_f_t0(u0, dt=dt) # Lambda @@ -115,7 +115,7 @@ def int_f(self, u0: np.ndarray, dt: float, t0: float = 0.0) -> np.ndarray: s = self.ext_scalar - retval = np.exp(dt * lam) * (u0 - s*np.sin(t0)) + s * np.sin(t0+dt) + retval = np.exp(dt * lam) * (u0 - s*np.sin(t)) + s * np.sin(t+dt) assert retval.shape == u0.shape return retval diff --git a/qmat/playgrounds/martin/diff_eqs/dahlquist2.py b/qmat/playgrounds/martin/diff_eqs/dahlquist2.py index 30bdbb1..1fa42be 100644 --- a/qmat/playgrounds/martin/diff_eqs/dahlquist2.py +++ b/qmat/playgrounds/martin/diff_eqs/dahlquist2.py @@ -36,8 +36,9 @@ def evalF(self, u: np.ndarray, t: float) -> np.ndarray: assert retval.shape == u.shape return retval - def int_f(self, u0: np.ndarray, t: float) -> np.ndarray: - assert isinstance(t, float) - retval = (self.s * np.exp(t * self.lam1) + (1.0 - self.s) * np.exp(t * self.lam2)) * u0 + 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 diff --git a/qmat/playgrounds/martin/diff_eqs/de_solver.py b/qmat/playgrounds/martin/diff_eqs/de_solver.py index f01ca31..e31b890 100644 --- a/qmat/playgrounds/martin/diff_eqs/de_solver.py +++ b/qmat/playgrounds/martin/diff_eqs/de_solver.py @@ -117,15 +117,17 @@ def initial_u0(self, mode: str) -> np.ndarray: pass @abstractmethod - def int_f(self, u0: np.ndarray, t: float) -> np.ndarray: + def int_f(self, u0: np.ndarray, dt: float, t: float = 0.0) -> np.ndarray: """ - Compute the (analytical) solution at time `t`. + Compute the (analytical) solution at time `dt` + t. Parameters ---------- u0 : np.ndarray Array of shape (N,) representing the initial condition. t : float - Time at which to evaluate the solution. + Time stamp of u0 + dt : float + Time step size """ pass diff --git a/qmat/playgrounds/martin/diff_eqs/two_freq.py b/qmat/playgrounds/martin/diff_eqs/two_freq.py index 0a9ea6a..0f045ee 100644 --- a/qmat/playgrounds/martin/diff_eqs/two_freq.py +++ b/qmat/playgrounds/martin/diff_eqs/two_freq.py @@ -90,10 +90,10 @@ def fSolve2(self, rhs: np.ndarray, dt: float, t: float) -> np.ndarray: assert retval.shape == rhs.shape return retval - def int_f(self, u0: np.ndarray, t: float) -> np.ndarray: + def int_f(self, u0: np.ndarray, dt: float, t: float = 0.0) -> np.ndarray: from scipy.linalg import expm - assert isinstance(t, float) - retval = expm(self.L * t) @ u0 + assert isinstance(dt, float) + retval = expm(self.L * (t + dt)) @ u0 assert retval.shape == u0.shape return retval diff --git a/qmat/playgrounds/martin/ex_burgers_integrate.py b/qmat/playgrounds/martin/ex_burgers_integrate.py index 6df3ee9..92d71ca 100644 --- a/qmat/playgrounds/martin/ex_burgers_integrate.py +++ b/qmat/playgrounds/martin/ex_burgers_integrate.py @@ -30,7 +30,7 @@ if 1: results = [] - u_analytical = burgers.int_f(u0, t=T) + u_analytical = burgers.int_f(u0, dt=T) for nt in range(4): diff --git a/qmat/playgrounds/martin/ex_burgers_plot.py b/qmat/playgrounds/martin/ex_burgers_plot.py index e15a6e9..1f854b4 100644 --- a/qmat/playgrounds/martin/ex_burgers_plot.py +++ b/qmat/playgrounds/martin/ex_burgers_plot.py @@ -32,7 +32,7 @@ plt.plot(t, ut) else: - ut = burgers.int_f(u0, t=T) + ut = burgers.int_f(u0, dt=T) plt.plot(t, u0, label="u0") plt.plot(t, ut, label="ut") diff --git a/qmat/playgrounds/martin/ex_dahlquist2_integrate.py b/qmat/playgrounds/martin/ex_dahlquist2_integrate.py index db0449f..b4291b3 100644 --- a/qmat/playgrounds/martin/ex_dahlquist2_integrate.py +++ b/qmat/playgrounds/martin/ex_dahlquist2_integrate.py @@ -25,7 +25,7 @@ dt = T / num_timesteps - u_analytical_fin = dahlquist2.int_f(u0, t=T) + u_analytical_fin = dahlquist2.int_f(u0, dt=T) u0 = dahlquist2.initial_u0() u = u0.copy() diff --git a/qmat/playgrounds/martin/ex_dahlquist_plot.py b/qmat/playgrounds/martin/ex_dahlquist_plot.py index bd6dacc..f9daa64 100644 --- a/qmat/playgrounds/martin/ex_dahlquist_plot.py +++ b/qmat/playgrounds/martin/ex_dahlquist_plot.py @@ -9,7 +9,7 @@ dahlquist: Dahlquist = Dahlquist(lam1=20.0j, lam2=1.0j) -u_eval = np.array([dahlquist.int_f(u0, _) for _ in t]) +u_eval = np.array([dahlquist.int_f(u0, dt=_) for _ in t]) plt.plot(t, np.real(u_eval), label="Re(u)") diff --git a/qmat/playgrounds/martin/ex_two_freq_integrate.py b/qmat/playgrounds/martin/ex_two_freq_integrate.py index 00fb13a..669dfa9 100644 --- a/qmat/playgrounds/martin/ex_two_freq_integrate.py +++ b/qmat/playgrounds/martin/ex_two_freq_integrate.py @@ -1,6 +1,6 @@ import numpy as np from qmat.playgrounds.martin.diff_eqs.two_freq import TwoFreq -from time_integration.sdc_integration import SDCIntegration +from qmat.playgrounds.martin.time_integration.sdc_integration import SDCIntegration from qmat.playgrounds.martin.time_integration.rk_integration import RKIntegration from matplotlib import pyplot as plt @@ -27,7 +27,7 @@ dt = T / num_timesteps - u_analytical_fin = two_freq.int_f(u0, t=T) + u_analytical_fin = two_freq.int_f(u0, dt=T) u0 = two_freq.initial_u0() u = u0.copy() diff --git a/qmat/playgrounds/martin/ex_two_freq_plot.py b/qmat/playgrounds/martin/ex_two_freq_plot.py index 117b29a..c695e8c 100644 --- a/qmat/playgrounds/martin/ex_two_freq_plot.py +++ b/qmat/playgrounds/martin/ex_two_freq_plot.py @@ -9,7 +9,7 @@ two_freq: TwoFreq = TwoFreq(lam1=1.0j, lam2=20.0j, lam3=0.5j) u0 = two_freq.initial_u0() -u_eval = np.array([two_freq.int_f(u0, _) for _ in t]) +u_eval = np.array([two_freq.int_f(u0, dt=_) for _ in t]) for i in range(2): plt.plot(t, np.real(u_eval[:, i]), label=f"Re(u[{i}])") diff --git a/qmat/playgrounds/martin/nyquist_test/ex_slow_fast.py b/qmat/playgrounds/martin/nyquist_test/ex_slow_fast.py new file mode 100644 index 0000000..d65aee1 --- /dev/null +++ b/qmat/playgrounds/martin/nyquist_test/ex_slow_fast.py @@ -0,0 +1,104 @@ +import numpy as np +from qmat.playgrounds.martin.diff_eqs.two_freq import TwoFreq +from qmat.playgrounds.martin.diff_eqs.dahlquist import Dahlquist +from qmat.playgrounds.martin.time_integration.sdc_integration import SDCIntegration +from qmat.playgrounds.martin.time_integration.rk_integration import RKIntegration +from matplotlib import pyplot as plt + + +s: int = 1 +T: float = 1.0 * np.pi * s # Time interval +t: float = 0.0 # Starting time + +time_integration: str = "sdc" +sdc_micro_time_integration: str = "irk1" +sdc_micro_time_integration: str = "erk1" +sdc_micro_time_integration: str = "imex21" +sdc_num_sweeps: int = 6 + +de_solver: Dahlquist = Dahlquist(lam1=1.0j, lam2=10.0j, ext_scalar=0.5) +u0 = de_solver.initial_u0() + +results = [] + +plt.close() + +if 1: + num_timesteps = 1 + print(f"Running simulation with num_timesteps={num_timesteps}") + + dt = T / num_timesteps + + u_analytical_fin = de_solver.int_f(u0, dt=T, t=0) + u0 = de_solver.initial_u0() + + u = u0.copy() + + u_: np.ndarray = np.array([u]) + t_: np.ndarray = np.array([t]) + + if time_integration in RKIntegration.supported_methods: + rki = RKIntegration(method=time_integration) + + for n in range(num_timesteps): + u = rki.int_f(u, t + n * dt, dt, de_solver) + + u_ = np.concatenate((u_, np.expand_dims(u, axis=0))) + t_ = np.append(t_, t + (n + 1) * dt) + + elif time_integration == "sdc": + sdci = SDCIntegration( + num_nodes=6, + # node_type="LEGENDRE", + node_type="EQUID", + quad_type="LOBATTO", + num_sweeps=sdc_num_sweeps, + micro_time_integration=sdc_micro_time_integration, + use_quadrature=False, + ) + + for n in range(num_timesteps): + u = sdci.int_f(u, t + n * dt, dt, de_solver) + + u_ = np.concatenate((u_, np.expand_dims(u, axis=0))) + t_ = np.append(t_, t + (n + 1) * dt) + + else: + raise Exception(f"Unsupported time integration method '{time_integration}'") + + # for i in range(u_.shape[1]): + for i in [0]: + # plt.plot(t_, np.real(u_[:, i]), label=f"u[{i}].real", linestyle="dashed") + plt.plot(t_, np.imag(u_[:, i]), label=f"u[{i}].imag", linestyle="solid") + + for t in t_: + plt.vlines(t, ymin=-1, ymax=1, colors="gray", linewidth=1) + + error = np.max(np.abs(u - u_analytical_fin)) + results.append({"N": num_timesteps, "dt": dt, "error": error}) + + +if 1: + # Plot analytical solution + t_ = np.linspace(0, T, 1000) + u_analytical_fin = np.array([de_solver.int_f(u0, dt=t, t=0) for t in t_]) + # for i in range(u_analytical_fin.shape[1]): + for i in [0]: + # plt.plot(t_, np.real(u_analytical_fin[:, i]), linestyle="dotted", color="black") + plt.plot(t_, np.imag(u_analytical_fin[:, i]), linestyle="dotted", color="black") + + +if 0: + 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"] + +if 1: + plt.legend() + plt.show() diff --git a/qmat/playgrounds/martin/tests/test_burgers.py b/qmat/playgrounds/martin/tests/test_burgers.py index 931a004..4517504 100644 --- a/qmat/playgrounds/martin/tests/test_burgers.py +++ b/qmat/playgrounds/martin/tests/test_burgers.py @@ -24,7 +24,7 @@ def test_burgers(): print("="*80) results = [] - u_analytical = burgers.int_f(u0, t=T) + u_analytical = burgers.int_f(u0, dt=T) for nt in range(2, 4): diff --git a/qmat/playgrounds/martin/tests/test_dahlquist.py b/qmat/playgrounds/martin/tests/test_dahlquist.py index 69612b1..e8ec065 100644 --- a/qmat/playgrounds/martin/tests/test_dahlquist.py +++ b/qmat/playgrounds/martin/tests/test_dahlquist.py @@ -96,7 +96,7 @@ def test_dahlquist_integration_of_solution(): u_analytical = dahlquist.int_f_t0(u0, dt=T) - u1 = dahlquist.int_f(u0, dt=T*0.5, t0=0.0) - u2 = dahlquist.int_f(u1, dt=T*0.5, t0=T*0.5) + u1 = dahlquist.int_f(u0, dt=T*0.5, t=0.0) + u2 = dahlquist.int_f(u1, dt=T*0.5, t=T*0.5) assert np.all(np.isclose(u_analytical, u2)) diff --git a/qmat/playgrounds/martin/tests/test_dahlquist2.py b/qmat/playgrounds/martin/tests/test_dahlquist2.py index efcb179..8936ee4 100644 --- a/qmat/playgrounds/martin/tests/test_dahlquist2.py +++ b/qmat/playgrounds/martin/tests/test_dahlquist2.py @@ -18,7 +18,7 @@ def test_dahlquist2(): print("="*80) results = [] - u_analytical = dahlquist2.int_f(u0, t=T) + u_analytical = dahlquist2.int_f(u0, dt=T) for nt in range(4): diff --git a/qmat/playgrounds/martin/tests/test_two_freq.py b/qmat/playgrounds/martin/tests/test_two_freq.py index 2aa3a30..58a056e 100644 --- a/qmat/playgrounds/martin/tests/test_two_freq.py +++ b/qmat/playgrounds/martin/tests/test_two_freq.py @@ -24,7 +24,7 @@ def test_dahlquist2(): print("=" * 80) results = [] - u_analytical = two_freq.int_f(u0, t=T) + u_analytical = two_freq.int_f(u0, dt=T) for nt in range(4): diff --git a/qmat/playgrounds/martin/time_integration/sdc_integration.py b/qmat/playgrounds/martin/time_integration/sdc_integration.py index 7a6bb09..a40f45d 100644 --- a/qmat/playgrounds/martin/time_integration/sdc_integration.py +++ b/qmat/playgrounds/martin/time_integration/sdc_integration.py @@ -10,6 +10,7 @@ def __init__( quad_type: str = "LOBATTO", # 'LOBATTO': Always include {0, 1} in quadrature points. Add them if they don't exist. num_sweeps: int = None, micro_time_integration: str = None, # 'erk1' = explicit Euler, 'irk1' = implicit Euler, 'imex' = implicit-explicit + use_quadrature: bool = True, # Use final quadrature to compute final solution ): from qmat.qcoeff.collocation import Collocation import qmat.qdelta.timestepping as module @@ -46,6 +47,8 @@ def __init__( if self.time_integration_method == "imex": self.time_integration_method = "imex12" + self.use_quadrature: bool = use_quadrature + assert self.num_sweeps >= 1 def integrate_erk1(self, u0: np.array, t: float, dt: float, de_solver: DESolver) -> np.array: @@ -87,13 +90,12 @@ def integrate_erk1(self, u0: np.array, t: float, dt: float, de_solver: DESolver) # Copy tendency arrays evalF_k0[:] = evalF_k1[:] - if 0: - # If we're using Radau-right, we can just use the last value - u[0] = u[-1] - - else: + if self.use_quadrature: # Compute new starting value with quadrature on tendencies u[0] = u[0] + dt * self.weights @ evalF_k0 + else: + # If we're using Radau-right, we can just use the last value + u[0] = u[-1] assert u0.shape == u[0].shape return u[0] @@ -168,14 +170,14 @@ def integrate_irk1(self, u0: np.array, t: float, dt: float, de_solver: DESolver) # Copy tendency arrays evalF_k0[:] = evalF_k1[:] - if 0: - # If we're using Radau-right, we can just use the last value - u[0] = u[-1] - - else: + if self.use_quadrature: # Compute new starting value with quadrature on tendencies u[0] = u[0] + dt * self.weights @ evalF_k0 + else: + # If we're using Radau-right, we can just use the last value + u[0] = u[-1] + assert u0.shape == u[0].shape return u[0] @@ -201,6 +203,8 @@ def integrate_imex21(self, u0: np.array, t: float, dt: float, de_solver: DESolve # Backup integrator contribution I[...] of previous iteration ISolves = np.empty_like(u) + print("="*80) + print("INITIAL") # # Propagate initial condition to all nodes # @@ -227,11 +231,15 @@ def integrate_imex21(self, u0: np.array, t: float, dt: float, de_solver: DESolve ISolves[m] = u[m] - u[m - 1] evalF_k0[m] = de_solver.evalF(u[m], t + dt * self.nodes[m]) + print(m, "eval_f", evalF_k0[m]) + print(" ", "u", u[m]) # # Iteratively sweep over SDC nodes # for _ in range(1, self.num_sweeps): + print("="*80) + print(f"SWEET {_}") for m in range(0, self.N): if m > 0: # @@ -256,17 +264,18 @@ def integrate_imex21(self, u0: np.array, t: float, dt: float, de_solver: DESolve ISolves[m] = u[m] - u[m - 1] + ISolves[m] - dt * qeval evalF_k1[m] = de_solver.evalF(u[m], t + dt * self.nodes[m]) + print(m, "eval_f", evalF_k1[m]) + print(" ", "u", u[m]) # Copy tendency arrays evalF_k0[:] = evalF_k1[:] - if 0: - # If we're using Radau-right, we can just use the last value - u[0] = u[-1] - - else: + if self.use_quadrature: # Compute new starting value with quadrature on tendencies u[0] = u[0] + dt * self.weights @ evalF_k0 + else: + # If we're using Radau-right, we can just use the last value + u[0] = u[-1] assert u0.shape == u[0].shape return u[0] @@ -352,13 +361,12 @@ def integrate_imex12(self, u0: np.array, t: float, dt: float, de_solver: DESolve # Copy tendency arrays evalF_k0[:] = evalF_k1[:] - if 0: - # If we're using Radau-right, we can just use the last value - u[0] = u[-1] - - else: + if self.use_quadrature: # Compute new starting value with quadrature on tendencies u[0] = u[0] + dt * self.weights @ evalF_k0 + else: + # If we're using Radau-right, we can just use the last value + u[0] = u[-1] assert u0.shape == u[0].shape return u[0]