{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# Solving the Poisson equation\n", "\n", "## Overview\n", "\n", "The **Poisson equation** is a fundamental elliptic partial differential equation that appears throughout physics and engineering:\n", "- Electrostatics (Coulomb's law)\n", "- Magnetostatics \n", "- Gravitational potential\n", "- Fluid dynamics (pressure Poisson equation)\n", "- And many coupled multiphysics systems\n", "\n", "This tutorial demonstrates how to solve various manufactured Poisson problems using Struphy's `Poisson` model, with verification against analytical solutions.\n", "\n", "## Solving 1D Poisson with periodic boundary conditions\n", "\n", "We begin with the simplest case: a **1D periodic Poisson problem**. The model solves the stabilized equation:\n", "\n", "$$\n", "-\\frac{d^2\\phi}{dx^2} + \\varepsilon\\phi = \\rho(x),\n", "$$\n", "\n", "where:\n", "- $\\phi(x)$ is the potential\n", "- $\\rho(x)$ is the source (charge/forcing)\n", "- $\\varepsilon$ is a small stabilization parameter (needed for uniqueness in periodic domains)\n", "\n", "For a manufactured test, we choose the analytical solution:\n", "\n", "$$\n", "\\phi_\\mathrm{exact}(x) = \\cos(kx),\n", "$$\n", "\n", "where $kx$ is the wavenumber. This determines the required source as:\n", "\n", "$$\n", "\\rho(x) = (k^2 + \\varepsilon)\\cos(kx).\n", "$$\n", "\n", "We solve on a periodic domain $x \\in [0, L)$, run a Struphy `Simulation`, and compare numerical and exact solutions." ] }, { "cell_type": "code", "execution_count": null, "id": "1", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "from struphy import Simulation, DerhamOptions, domains, grids, perturbations\n", "from struphy.initial.base import GenericPerturbation\n", "from struphy.models import Poisson\n", "\n", "Poisson.pde()" ] }, { "cell_type": "code", "execution_count": null, "id": "2", "metadata": {}, "outputs": [], "source": [ "# Manufactured periodic test case in 1D (aligned with verification test pattern)\n", "Lx = 2.0 * np.pi\n", "mode = 2\n", "k = mode * 2.0 * np.pi / Lx\n", "\n", "# Tiny stabilization makes the periodic problem uniquely solvable\n", "stab_eps = 1e-8\n", "\n", "# Build the source through the model variable as in test_verif_Poisson.py\n", "model = Poisson()\n", "model.propagators.poisson.options = model.propagators.poisson.Options(\n", " stab_eps=stab_eps,\n", " )\n", "\n", "source_amp = k**2 + stab_eps\n", "model.em_fields.source.add_perturbation(\n", " perturbations.ModesCos(ls=(mode,), amps=(source_amp,))\n", " )\n", "\n", "phi_exact = lambda e1, e2, e3: np.cos(k * e1)" ] }, { "cell_type": "markdown", "id": "3", "metadata": {}, "source": [ "### Setting up the 1D periodic case\n", "\n", "To run a Poisson solve in Struphy, we need to:\n", "\n", "1. **Create a model** with the appropriate equation parameters\n", "2. **Define a manufactured source** that will be used as the right-hand side\n", "3. **Specify solver options** (stabilization, numerical parameters)\n", "4. **Create a simulation** with domain, grid, and solver settings\n", "5. **Execute the solve** and post-process results\n", "\n", "Below, `Poisson.Options` controls the elliptic solve. We specify:\n", "\n", "- `stab_eps`: a tiny stabilization to remove the constant null-space in the periodic case\n", "\n", "The source term (`rho`) is wired up automatically when the `Poisson` model is constructed (it passes `em_fields.source` to the propagator internally).\n", "All other option fields use defaults (solver choice, preconditioner, tolerances, etc.)." ] }, { "cell_type": "code", "execution_count": null, "id": "4", "metadata": {}, "outputs": [], "source": [ "# 1D periodic setup: periodic bcs are the default (None in each direction)\n", "domain = domains.Cuboid(l1=0.0, r1=Lx)\n", "grid = grids.TensorProductGrid(num_elements=(64, 1, 1))\n", "derham_opts = DerhamOptions()\n", "\n", "sim = Simulation(\n", " model=model,\n", " domain=domain,\n", " grid=grid,\n", " derham_opts=derham_opts,\n", " )\n", "\n", "# For a stationary Poisson solve, one step is enough\n", "sim.run(one_time_step=True)\n", "sim.pproc()\n", "sim.load_plotting_data()" ] }, { "cell_type": "code", "execution_count": null, "id": "5", "metadata": {}, "outputs": [], "source": [ "# Extract 1D line data and compare to analytic solution\n", "x = sim.grids_phy[0][:, 0, 0]\n", "\n", "t_last = max(sim.spline_values.em_fields.phi_log.data.keys())\n", "phi_num = sim.spline_values.em_fields.phi_log.data[t_last][0][:, 0, 0]\n", "phi_ref = phi_exact(x, 0.0, 0.0)\n", "\n", "err = phi_num - phi_ref\n", "err_max = np.max(np.abs(err))\n", "\n", "plt.figure(figsize=(9, 4))\n", "plt.subplot(1, 2, 1)\n", "plt.plot(x, phi_ref, \"k--\", label=\"analytic\")\n", "plt.plot(x, phi_num, \"o\", ms=3, label=\"Struphy\")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(r\"$\\phi(x)$\")\n", "plt.title(\"Poisson solution\")\n", "plt.legend()\n", "plt.grid(alpha=0.3)\n", "\n", "plt.subplot(1, 2, 2)\n", "plt.plot(x, err, \"r\", lw=1.5)\n", "plt.xlabel(\"x\")\n", "plt.ylabel(r\"$\\phi_h - \\phi_\\mathrm{exact}$\")\n", "plt.title(\"Pointwise error\")\n", "plt.grid(alpha=0.3)\n", "\n", "plt.tight_layout()\n", "print(f\"max-norm error ||phi_h - phi_exact||_inf = {err_max:.3e}\")\n", "\n", "# Electric field comparison: E = -grad(phi)\n", "E_num = -np.gradient(phi_num, x, edge_order=2)\n", "E_ref = k * np.sin(k * x)\n", "E_err = E_num - E_ref\n", "E_err_max = np.max(np.abs(E_err))\n", "\n", "plt.figure(figsize=(9, 4))\n", "plt.subplot(1, 2, 1)\n", "plt.plot(x, E_ref, \"k--\", label=r\"analytic $-\\nabla\\phi$\")\n", "plt.plot(x, E_num, \"o\", ms=3, label=r\"Struphy $-\\nabla\\phi_h$\")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(r\"$E_x(x)$\")\n", "plt.title(\"Electric field comparison\")\n", "plt.legend()\n", "plt.grid(alpha=0.3)\n", "\n", "plt.subplot(1, 2, 2)\n", "plt.plot(x, E_err, \"m\", lw=1.5)\n", "plt.xlabel(\"x\")\n", "plt.ylabel(r\"$E_{x,h} - E_{x,\\mathrm{exact}}$\")\n", "plt.title(\"Electric field pointwise error\")\n", "plt.grid(alpha=0.3)\n", "\n", "plt.tight_layout()\n", "print(f\"max-norm error ||E_h - E_exact||_inf = {E_err_max:.3e}\")" ] }, { "cell_type": "markdown", "id": "6", "metadata": {}, "source": [ "## 2D Manufactured Test on a Rectangle\n", "\n", "We now solve a **2D Poisson problem on a rectangular domain** with different side lengths. This tests:\n", "- Multi-dimensional elliptic solves\n", "- Non-square domains\n", "- Homogeneous Dirichlet (zero) boundary conditions\n", "\n", "The rectangular domain is:\n", "\n", "$$\n", "(x,y) \\in [0,L_x] \\times [0,L_y], \\quad L_x \\neq L_y.\n", "$$\n", "\n", "We choose a separable manufactured solution that vanishes on all four boundaries:\n", "\n", "$$\n", "\\phi_\\mathrm{exact}(x,y) = \\sin\\left(\\frac{2\\pi x}{L_x}\\right)\\sin\\left(\\frac{2\\pi y}{L_y}\\right).\n", "$$\n", "\n", "For the Poisson equation $-\\Delta\\phi + \\varepsilon\\phi = \\rho$, the required manufactured source is:\n", "\n", "$$\n", "\\rho(x,y) = \\left[\\left(\\frac{2\\pi}{L_x}\\right)^2 + \\left(\\frac{2\\pi}{L_y}\\right)^2 + \\varepsilon\\right] \\phi_\\mathrm{exact}(x,y).\n", "$$\n", "\n", "Unlike the periodic case, we now enforce **Dirichlet boundary conditions** ($\\phi = 0$ on all boundaries).\n", "The solution is unique without any stabilization ($\\varepsilon$ can be very small).\n", "\n", "We construct this 2D test through model perturbations (`ModesSinSin`) and compare solution fields at select points using contour plots." ] }, { "cell_type": "code", "execution_count": null, "id": "7", "metadata": {}, "outputs": [], "source": [ "# 2D manufactured Poisson setup and solve\n", "Lx2 = 2.0\n", "Ly2 = 1.0\n", "eps2 = 1e-8\n", "\n", "kx2 = 2.0 * np.pi / Lx2\n", "ky2 = 2.0 * np.pi / Ly2\n", "\n", "phi2_exact = lambda e1, e2, e3: np.sin(kx2 * e1) * np.sin(ky2 * e2)\n", "\n", "model2 = Poisson()\n", "model2.propagators.poisson.options = model2.propagators.poisson.Options(\n", " stab_eps=eps2,\n", " )\n", "\n", "source_amp2 = kx2**2 + ky2**2 + eps2\n", "model2.em_fields.source.add_perturbation(\n", " perturbations.ModesSinSin(ls=(1,), ms=(1,), amps=(source_amp2,))\n", " )\n", "\n", "domain2 = domains.Cuboid(l1=0.0, r1=Lx2, l2=0.0, r2=Ly2)\n", "grid2 = grids.TensorProductGrid(num_elements=(48, 32, 1))\n", "derham_opts2 = DerhamOptions(\n", " degree=(3, 3, 1),\n", " bcs=((\"dirichlet\", \"dirichlet\"), (\"dirichlet\", \"dirichlet\"), None),\n", " )\n", "\n", "sim2 = Simulation(\n", " model=model2,\n", " domain=domain2,\n", " grid=grid2,\n", " derham_opts=derham_opts2,\n", " )\n", "\n", "sim2.run(one_time_step=True)\n", "sim2.pproc()\n", "sim2.load_plotting_data()" ] }, { "cell_type": "code", "execution_count": null, "id": "8", "metadata": {}, "outputs": [], "source": [ "# 2D diagnostics and plots\n", "t2_last = max(sim2.spline_values.em_fields.phi_log.data.keys())\n", "X = sim2.grids_phy[0][:, :, 0]\n", "Y = sim2.grids_phy[1][:, :, 0]\n", "\n", "phi2_num = sim2.spline_values.em_fields.phi_log.data[t2_last][0][:, :, 0]\n", "phi2_ref = phi2_exact(X, Y, 0.0)\n", "err2 = phi2_num - phi2_ref\n", "err2_max = np.max(np.abs(err2))\n", "\n", "fig, axs = plt.subplots(1, 3, figsize=(15, 4.5), constrained_layout=True)\n", "\n", "im0 = axs[0].contourf(X, Y, phi2_num, levels=40, cmap=\"viridis\")\n", "axs[0].set_title(\"Numerical $\\\\phi_h(x,y)$\")\n", "axs[0].set_xlabel(\"x\")\n", "axs[0].set_ylabel(\"y\")\n", "plt.colorbar(im0, ax=axs[0])\n", "\n", "im1 = axs[1].contourf(X, Y, phi2_ref, levels=40, cmap=\"viridis\")\n", "axs[1].set_title(\"Analytic $\\\\phi_{\\\\mathrm{exact}}(x,y)$\")\n", "axs[1].set_xlabel(\"x\")\n", "axs[1].set_ylabel(\"y\")\n", "plt.colorbar(im1, ax=axs[1])\n", "\n", "im2 = axs[2].contourf(X, Y, err2, levels=40, cmap=\"coolwarm\")\n", "axs[2].set_title(\"Error $\\\\phi_h-\\\\phi_{\\\\mathrm{exact}}$\")\n", "axs[2].set_xlabel(\"x\")\n", "axs[2].set_ylabel(\"y\")\n", "plt.colorbar(im2, ax=axs[2])\n", "\n", "print(f\"2D max-norm error ||phi_h - phi_exact||_inf = {err2_max:.3e}\")" ] }, { "cell_type": "markdown", "id": "9", "metadata": {}, "source": [ "## 2D Manufactured Test on a Thin Annulus (Polar Coordinates)\n", "\n", "As a second 2D geometry test, we solve Poisson on a **thin annulus** (a disc-shaped domain with a small hole).\n", "This demonstrates:\n", "- Non-Cartesian geometries (polar coordinates via conformal mapping)\n", "- Dirichlet boundary conditions at both inner and outer radii\n", "- Manufactured solutions expressed in polar coordinates\n", "\n", "The geometry is a thin annulus with radii $a_1 < r < a_2$. We use the manufactured solution in polar coordinates $(r, \\theta)$:\n", "\n", "$$\n", "\\phi_\\mathrm{exact}(r,\\theta) = \\sin\\left(\\pi\\frac{r-a_1}{a_2-a_1}\\right)\\sin(m\\theta),\n", "$$\n", "\n", "where $m$ is an azimuthal mode number. This solution automatically satisfies:\n", "- **Inner boundary** ($r=a_1$): $\\phi = 0$ (since $\\sin(0) = 0$)\n", "- **Outer boundary** ($r=a_2$): $\\phi = 0$ (since $\\sin(\\pi) = 0$)\n", "\n", "The source term $\\rho$ is computed analytically in polar coordinates, accounting for the Laplacian in polar form:\n", "\n", "$$\n", "\\Delta\\phi = \\frac{d^2\\phi}{dr^2} + \\frac{1}{r}\\frac{d\\phi}{dr} + \\frac{1}{r^2}\\frac{d^2\\phi}{d\\theta^2}.\n", "$$\n", "\n", "This manufactured solution is injected via a `GenericPerturbation` directly in physical space.\n", "Plots are shown in physical $(x,y)$ coordinates for visualization." ] }, { "cell_type": "code", "execution_count": null, "id": "10", "metadata": {}, "outputs": [], "source": [ "# HollowCylinder manufactured annulus test (physical-space definition)\n", "a1 = 0.8\n", "a2 = 1.0\n", "Lz3 = 1.0\n", "m3 = 3\n", "eps3 = 1e-8\n", "w3 = a2 - a1\n", "alpha3 = np.pi / w3\n", "\n", "def phi3_exact(x, y, z):\n", " r = np.sqrt(x**2 + y**2)\n", " theta = np.arctan2(y, x)\n", " s = alpha3 * (r - a1)\n", " return np.sin(s) * np.sin(m3 * theta)\n", "\n", "def rho3_exact(x, y, z):\n", " r = np.sqrt(x**2 + y**2)\n", " theta = np.arctan2(y, x)\n", " s = alpha3 * (r - a1)\n", "\n", " sin_s = np.sin(s)\n", " cos_s = np.cos(s)\n", " sin_mt = np.sin(m3 * theta)\n", "\n", " # -Delta(phi) + eps*phi for phi(r,theta)=sin(s)sin(m theta)\n", " term = (alpha3**2) * sin_s - (alpha3 / r) * cos_s + (m3**2 / r**2) * sin_s + eps3 * sin_s\n", " return term * sin_mt\n", "\n", "model3 = Poisson()\n", "model3.propagators.poisson.options = model3.propagators.poisson.Options(\n", " stab_eps=eps3,\n", " )\n", "\n", "model3.em_fields.source.add_perturbation(\n", " GenericPerturbation(rho3_exact, given_in_basis=\"physical\")\n", " )\n", "\n", "domain3 = domains.HollowCylinder(a1=a1, a2=a2, Lz=Lz3)\n", "grid3 = grids.TensorProductGrid(num_elements=(36, 96, 1))\n", "derham_opts3 = DerhamOptions(\n", " degree=(3, 3, 1),\n", " bcs=((\"dirichlet\", \"dirichlet\"), None, None),\n", " )\n", "\n", "sim3 = Simulation(\n", " model=model3,\n", " domain=domain3,\n", " grid=grid3,\n", " derham_opts=derham_opts3,\n", " )\n", "\n", "sim3.run(one_time_step=True)\n", "sim3.pproc()\n", "sim3.load_plotting_data()" ] }, { "cell_type": "code", "execution_count": null, "id": "11", "metadata": {}, "outputs": [], "source": [ "# Annulus diagnostics and plots in physical coordinates only\n", "t3_last = max(sim3.spline_values.em_fields.phi_log.data.keys())\n", "X3 = sim3.grids_phy[0][:, :, 0]\n", "Y3 = sim3.grids_phy[1][:, :, 0]\n", "\n", "phi3_num = sim3.spline_values.em_fields.phi_log.data[t3_last][0][:, :, 0]\n", "phi3_ref = phi3_exact(X3, Y3, 0.0)\n", "err3 = phi3_num - phi3_ref\n", "err3_max = np.max(np.abs(err3))\n", "\n", "fig3, axs3 = plt.subplots(1, 3, figsize=(15, 4.8), constrained_layout=True)\n", "\n", "im30 = axs3[0].contourf(X3, Y3, phi3_num, levels=40, cmap=\"viridis\")\n", "axs3[0].set_title(\"Numerical $\\\\phi_h(x,y)$ on annulus\")\n", "axs3[0].set_xlabel(\"x\")\n", "axs3[0].set_ylabel(\"y\")\n", "axs3[0].set_aspect(\"equal\")\n", "plt.colorbar(im30, ax=axs3[0])\n", "\n", "im31 = axs3[1].contourf(X3, Y3, phi3_ref, levels=40, cmap=\"viridis\")\n", "axs3[1].set_title(\"Analytic $\\\\phi_{\\\\mathrm{exact}}(x,y)$ on annulus\")\n", "axs3[1].set_xlabel(\"x\")\n", "axs3[1].set_ylabel(\"y\")\n", "axs3[1].set_aspect(\"equal\")\n", "plt.colorbar(im31, ax=axs3[1])\n", "\n", "im32 = axs3[2].contourf(X3, Y3, err3, levels=40, cmap=\"coolwarm\")\n", "axs3[2].set_title(\"Error $\\\\phi_h-\\\\phi_{\\\\mathrm{exact}}$\")\n", "axs3[2].set_xlabel(\"x\")\n", "axs3[2].set_ylabel(\"y\")\n", "axs3[2].set_aspect(\"equal\")\n", "plt.colorbar(im32, ax=axs3[2])\n", "\n", "print(f\"Annulus max-norm error ||phi_h - phi_exact||_inf = {err3_max:.3e}\")" ] }, { "cell_type": "markdown", "id": "12", "metadata": {}, "source": [ "## 1D Time-Dependent Source Test Case\n", "\n", "So far we have solved stationary Poisson problems where the source term is independent of time.\n", "Now we consider a **time-dependent source** where the right-hand side of the Poisson equation evolves in time.\n", "\n", "This scenario arises in many applications, such as:\n", "- Plasma physics simulations where the source term varies with time due to wave activity\n", "- Electromagnetic problems with oscillating charge or current distributions\n", "- Coupled multiphysics systems where sources depend on time-evolving quantities\n", "\n", "### Mathematical Setup\n", "\n", "We solve the time-dependent Poisson equation on a 1D periodic domain with a **sinusoidal source**:\n", "\n", "$$\n", "-\\phi'' + \\varepsilon\\phi = \\rho(x,t),\n", "$$\n", "\n", "where the manufactured source is:\n", "\n", "$$\n", "\\rho(x,t) = A\\cos\\left(kx\\right)\\cos(\\omega t),\n", "$$\n", "\n", "with wavenumber $k$, spatial amplitude $A$, and angular frequency $\\omega$.\n", "\n", "The exact solution is:\n", "\n", "$$\n", "\\phi_\\mathrm{exact}(x,t) = \\frac{A}{k^2 + \\varepsilon}\\cos(kx)\\cos(\\omega t).\n", "$$\n", "\n", "The source varies periodically in time, and the solution follows this oscillation. We use Struphy's `TimeDependentSource` propagator to evolve both the source and the solution over multiple time steps, comparing numerical and exact solutions at each time." ] }, { "cell_type": "code", "execution_count": null, "id": "13", "metadata": {}, "outputs": [], "source": [ "# Time-dependent source setup\n", "# Parameters: domain extent, wavenumber, frequency, and stabilization\n", "Lx4 = 10.0 # Domain: x in [-Lx4/2, Lx4/2]\n", "l4 = 2 # Wavenumber mode: k = 2*pi*l/Lx\n", "k4 = l4 * 2.0 * np.pi / Lx4\n", "omega4 = 2.0 * np.pi # Angular frequency (one period in [0,1] of normalized time)\n", "amp4 = 0.1 # Spatial amplitude of source\n", "eps4 = 1e-8 # Tiny stabilization for periodic domain\n", "\n", "# Create a Poisson model with time-dependent source enabled\n", "model4 = Poisson(with_t_dep_source=True)\n", "\n", "# Configure the time-dependent source propagator with the frequency\n", "model4.propagators.source.options = model4.propagators.source.Options(omega=omega4)\n", "\n", "# Configure the Poisson solver: set the source field as the RHS\n", "model4.propagators.poisson.options = model4.propagators.poisson.Options(\n", " stab_eps=eps4,\n", ")\n", "\n", "# Initialize the source with a cosine perturbation\n", "# The source will be: rho(x,t) = amp * cos(k*x) * cos(omega*t)\n", "model4.em_fields.source.add_perturbation(\n", " perturbations.ModesCos(ls=(l4,), amps=(amp4,))\n", ")\n", "\n", "# Define exact solutions for comparison\n", "def rho4_exact(x, t):\n", " \"\"\"Exact time-dependent source: amplitude * cos(kx) * cos(omega*t)\"\"\"\n", " return amp4 * np.cos(k4 * x) * np.cos(omega4 * t)\n", "\n", "def phi4_exact(x, t):\n", " \"\"\"Exact solution: (amplitude / k^2) * cos(kx) * cos(omega*t)\"\"\"\n", " # The Poisson equation -phi'' + eps*phi = rho has solution\n", " # phi(x,t) = rho(x,t) / (k^2 + eps) for this cosine mode\n", " return (amp4 / (k4**2 + eps4)) * np.cos(k4 * x) * np.cos(omega4 * t)" ] }, { "cell_type": "code", "execution_count": null, "id": "14", "metadata": {}, "outputs": [], "source": [ "# Set up domain, grid, and simulation for time-dependent case\n", "l1_4 = -Lx4 / 2.0\n", "r1_4 = Lx4 / 2.0\n", "domain4 = domains.Cuboid(l1=l1_4, r1=r1_4)\n", "\n", "# Use a moderate grid resolution for reasonable compute time\n", "grid4 = grids.TensorProductGrid(num_elements=(48, 1, 1))\n", "\n", "# Time stepping: evolve from t=0 to t=2.0 with dt=0.1 (20 time steps)\n", "# This covers 2 full periods of oscillation (since omega = 2*pi)\n", "from struphy import Time\n", "time_opts4 = Time(dt=0.1, Tend=2.0)\n", "\n", "derham_opts4 = DerhamOptions()\n", "\n", "# Create and run the simulation\n", "sim4 = Simulation(\n", " model=model4,\n", " domain=domain4,\n", " grid=grid4,\n", " derham_opts=derham_opts4,\n", " time_opts=time_opts4,\n", ")\n", "\n", "# Run the full time-dependent simulation\n", "sim4.run()\n", "sim4.pproc()\n", "sim4.load_plotting_data()" ] }, { "cell_type": "code", "execution_count": null, "id": "15", "metadata": {}, "outputs": [], "source": [ "# Extract and visualize time-dependent results\n", "x4 = sim4.grids_phy[0][:, 0, 0]\n", "phi4_log = sim4.spline_values.em_fields.phi_log.data\n", "source4_log = sim4.spline_values.em_fields.source_log.data\n", "t_times = sorted(phi4_log.keys())\n", "\n", "print(f\"Solution saved at {len(t_times)} time points\")\n", "\n", "# Compute maximum error over all time steps\n", "err_max_global = 0.0\n", "for t in t_times:\n", " phi_h = phi4_log[t][0][:, 0, 0] # Numerical solution\n", " phi_e = phi4_exact(x4, t) # Exact solution\n", " err_local = np.max(np.abs(phi_h - phi_e))\n", " if err_local > err_max_global:\n", " err_max_global = err_local\n", "\n", "# Normalize by the solution amplitude for a relative error measure\n", "phi_amplitude = amp4 / (k4**2 + eps4)\n", "rel_err = err_max_global / phi_amplitude\n", "\n", "print(f\"Global max-norm error: {err_max_global:.3e}\")\n", "print(f\"Relative error: {rel_err:.3e} (normalized by solution amplitude {phi_amplitude:.3e})\")\n", "\n", "# Create time snapshots at select time steps\n", "fig_snaps, axs_snaps = plt.subplots(2, 3, figsize=(15, 8), constrained_layout=True)\n", "axs_snaps = axs_snaps.flatten()\n", "\n", "select_times = [t_times[i] for i in [0, len(t_times)//5, 2*len(t_times)//5, \n", " 3*len(t_times)//5, 4*len(t_times)//5, -1]]\n", "\n", "for idx, t in enumerate(select_times):\n", " ax = axs_snaps[idx]\n", " \n", " phi_num = phi4_log[t][0][:, 0, 0]\n", " phi_ref = phi4_exact(x4, t)\n", " \n", " ax.plot(x4, phi_ref, \"k--\", linewidth=2, label=\"Exact\")\n", " ax.plot(x4, phi_num, \"bo\", markersize=4, alpha=0.6, label=\"Struphy\")\n", " ax.set_xlabel(\"x\")\n", " ax.set_ylabel(r\"$\\phi(x,t)$\")\n", " ax.set_title(f\"Time $t = {t:.2f}$\")\n", " ax.legend()\n", " ax.grid(alpha=0.3)\n", " ax.set_ylim(-phi_amplitude * 1.2, phi_amplitude * 1.2)" ] }, { "cell_type": "code", "execution_count": null, "id": "16", "metadata": {}, "outputs": [], "source": [ "# Plot error evolution over time\n", "fig_err, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 4.5), constrained_layout=True)\n", "\n", "errors = []\n", "times_array = np.array(t_times)\n", "for t in t_times:\n", " phi_h = phi4_log[t][0][:, 0, 0]\n", " phi_e = phi4_exact(x4, t)\n", " err = np.max(np.abs(phi_h - phi_e))\n", " errors.append(err)\n", "\n", "ax1.plot(times_array, errors, \"ro-\", linewidth=2, markersize=6)\n", "ax1.set_xlabel(\"Time $t$\")\n", "ax1.set_ylabel(\"Max-norm error $||\\\\phi_h - \\\\phi_{\\\\mathrm{exact}}||_\\\\infty$\")\n", "ax1.set_title(\"Pointwise error evolution in time\")\n", "ax1.grid(alpha=0.3)\n", "\n", "# Plot source and solution at a single point in space (center)\n", "center_idx = len(x4) // 2\n", "phi_center = [phi4_log[t][0][center_idx, 0, 0] for t in t_times]\n", "phi_exact_center = [phi4_exact(x4[center_idx], t) for t in t_times]\n", "source_center = [source4_log[t][0][center_idx, 0, 0] for t in t_times]\n", "source_exact_center = [rho4_exact(x4[center_idx], t) for t in t_times]\n", "\n", "ax2.plot(times_array, phi_exact_center, \"k--\", linewidth=2.5, label=\"$\\\\phi_{\\\\mathrm{exact}}$ at center\")\n", "ax2.plot(times_array, phi_center, \"bo-\", markersize=4, alpha=0.7, label=\"$\\\\phi_h$ at center (Struphy)\")\n", "ax2.set_xlabel(\"Time $t$\")\n", "ax2.set_ylabel(r\"$\\phi$ at $x = 0$\")\n", "ax2.set_title(\"Solution time evolution at domain center\")\n", "ax2.legend()\n", "ax2.grid(alpha=0.3)" ] } ], "metadata": { "kernelspec": { "display_name": "env (3.12.3)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }