{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": { "language": "markdown" }, "source": [ "# Maxwell's equations: Light Wave and Coaxial Waveguide\n", "\n", "This combined tutorial covers two complementary verification scenarios for the `Maxwell` model in Struphy:\n", "\n", "1. **1D Light wave in vacuum**: verifies free-space propagation and the linear dispersion relation $\\omega = c k$.\n", "2. **2D Coaxial waveguide mode**: verifies structured electromagnetic eigenmodes in hollow-cylinder geometry.\n", "\n", "Before running either test, we print the PDE system in Struphy units to connect the numerical setup to the governing physics." ] }, { "cell_type": "code", "execution_count": null, "id": "1", "metadata": { "language": "python" }, "outputs": [], "source": [ "from struphy.models import Maxwell\n", "\n", "# Display Maxwell equations in Struphy normalization.\n", "Maxwell.pde()" ] }, { "cell_type": "code", "execution_count": null, "id": "2", "metadata": {}, "outputs": [], "source": [ "# Display normalization\n", "Maxwell.normalization()" ] }, { "cell_type": "markdown", "id": "3", "metadata": { "language": "markdown" }, "source": [ "## How The Two Tests Relate To The PDEs\n", "\n", "- In the light-wave test, the PDEs reduce to wave propagation in a homogeneous medium, so the dominant verification target is the dispersion relation.\n", "- In the coaxial test, the same PDEs are solved with conducting-boundary geometry constraints, so the dominant verification target is mode structure and boundary-consistent field patterns." ] }, { "cell_type": "markdown", "id": "4", "metadata": {}, "source": [ "## Verification of Electromagnetic Wave Propagation\n", "\n", "### 1D Light Wave in Vacuum\n", "\n", "This tutorial demonstrates verification of the dispersion relation for electromagnetic waves in vacuum using the `Maxwell` model. We initialize a broadband spectrum of electric field perturbations and extract the dominant wave frequencies using FFT power spectrum analysis.\n", "\n", "### Physical Setup\n", "\n", "In vacuum (or non-magnetic plasma), Maxwell's equations yield the free-space wave equation:\n", "\n", "$$\\nabla^2 \\mathbf{E} - \\frac{1}{c^2} \\frac{\\partial^2 \\mathbf{E}}{\\partial t^2} = 0,$$\n", "\n", "where $c$ is the speed of light (normalized to unity in Struphy units). Plane wave solutions have the dispersion relation\n", "\n", "$$\\omega = c \\, k,$$\n", "\n", "where $\\omega$ is the angular frequency and $k$ is the wavenumber.\n", "\n", "We verify this linear dispersion by:\n", "1. Initializing small random perturbations to the electric field\n", "2. Running a transient simulation\n", "3. Extracting the dominant frequency via FFT\n", "4. Fitting the dispersion relation to extract phase velocity\n", "5. Comparing against the theoretical speed $c = 1$" ] }, { "cell_type": "code", "execution_count": null, "id": "5", "metadata": {}, "outputs": [], "source": [ "import logging\n", "import os\n", "import shutil\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import cunumpy as xp\n", "\n", "from struphy import (\n", " DerhamOptions,\n", " EnvironmentOptions,\n", " Simulation,\n", " Time,\n", " domains,\n", " grids,\n", " perturbations,\n", ")\n", "from struphy.diagnostics.diagn_tools import power_spectrum_2d\n", "from struphy.models import Maxwell\n", "\n", "logger = logging.getLogger(\"struphy\")" ] }, { "cell_type": "markdown", "id": "6", "metadata": {}, "source": [ "### Model and Time-Stepping Configuration\n", "\n", "Create a Maxwell model instance and configure the time-stepping propagator. We use implicit time-stepping for stability over the long simulation window required to collect many wave cycles." ] }, { "cell_type": "code", "execution_count": null, "id": "7", "metadata": {}, "outputs": [], "source": [ "# Model instance\n", "model = Maxwell()\n", "\n", "# Choose time-stepping algorithm\n", "algo = \"implicit\" # Alternative: \"explicit\"\n", "\n", "# Propagator options\n", "model.propagators.maxwell.options = model.propagators.maxwell.Options(algo=algo)\n", "\n", "print(f\"Time-stepping algorithm: {algo}\")" ] }, { "cell_type": "markdown", "id": "8", "metadata": {}, "source": [ "### Domain and Numerical Discretization\n", "\n", "Set up a 1D periodic domain in the $z$-direction with a high-resolution grid. Electromagnetic waves require fine spatial resolution to represent the wave profile accurately." ] }, { "cell_type": "code", "execution_count": null, "id": "9", "metadata": {}, "outputs": [], "source": [ "# 1D domain (extended in z for long-wavelength waves)\n", "domain = domains.Cuboid(r3=20.0) # z ∈ [0, 20)\n", "\n", "# High-resolution grid in z (128 elements for wave smoothness)\n", "grid = grids.TensorProductGrid(num_elements=(1, 1, 128))\n", "\n", "# Derham options: cubic edges in z for accurate wave representation\n", "derham_opts = DerhamOptions(degree=(1, 1, 3))\n", "\n", "print(f\"Domain: z ∈ [0, {domain.params['r3']})\")\n", "print(f\"Grid elements: {grid.num_elements}\")\n", "print(f\"Derham degree: {derham_opts.degree}\")" ] }, { "cell_type": "markdown", "id": "10", "metadata": {}, "source": [ "### Initial Conditions\n", "\n", "Add small random noise perturbations to the $x$ and $y$ components of the electric field. These broadband perturbations excite a spectrum of electromagnetic waves that we will analyze via FFT." ] }, { "cell_type": "code", "execution_count": null, "id": "11", "metadata": {}, "outputs": [], "source": [ "# Electric field perturbations (small noise to excite waves)\n", "model.em_fields.e_field.add_perturbation(perturbations.Noise(amp=0.1, comp=0, seed=123))\n", "model.em_fields.e_field.add_perturbation(perturbations.Noise(amp=0.1, comp=1, seed=123))" ] }, { "cell_type": "markdown", "id": "12", "metadata": {}, "source": [ "### Simulation Setup and Execution\n", "\n", "Configure the simulation environment and run the transient dynamics. We evolve for a long time to collect sufficient wave cycles (Tend=50) for accurate FFT-based dispersion analysis." ] }, { "cell_type": "code", "execution_count": null, "id": "13", "metadata": {}, "outputs": [], "source": [ "# Environment and file management\n", "test_folder = os.path.join(os.getcwd(), \"struphy_verification_tests\")\n", "out_folders = os.path.join(test_folder, \"Maxwell\")\n", "env = EnvironmentOptions(out_folders=out_folders, sim_folder=\"light_wave_1d\")\n", "\n", "# Time-stepping: small dt for stability and accuracy\n", "time_opts = Time(dt=0.05, Tend=50.0)\n", "\n", "# Instantiate and run simulation\n", "sim = Simulation(\n", " model=model,\n", " env=env,\n", " time_opts=time_opts,\n", " domain=domain,\n", " grid=grid,\n", " derham_opts=derham_opts,\n", ")\n", "\n", "print(f\"Running simulation: dt={time_opts.dt}, Tend={time_opts.Tend}\")\n", "sim.run()\n", "print(\"Simulation complete.\")\n", "\n", "# Post-processing\n", "sim.pproc()\n", "print(\"Post-processing complete.\")" ] }, { "cell_type": "markdown", "id": "14", "metadata": {}, "source": [ "### Diagnostics: Dispersion Relation Verification\n", "\n", "Extract the electric field time evolution, compute the power spectrum via FFT, and fit the dominant frequency peaks. From the fitted frequency and spatial wavenumber, we extract the phase velocity and compare against the theoretical speed $c = 1$.\n", "\n", "**Expected result**: The fitted wave speed should be approximately 1.0 (the speed of light in normalized units)." ] }, { "cell_type": "code", "execution_count": null, "id": "15", "metadata": {}, "outputs": [], "source": [ "# Load plotting data\n", "sim.load_plotting_data()\n", "\n", "# Extract electric field time-series\n", "E_of_t = sim.spline_values.em_fields.e_field_log.data\n", "\n", "# Compute power spectrum and fit dispersion relation\n", "print(\"\\n=== Light Wave Dispersion Analysis ===\")\n", "_1, _2, _3, coeffs = power_spectrum_2d(\n", " E_of_t,\n", " \"e_field_log\",\n", " grids=sim.grids_log,\n", " grids_mapped=sim.grids_phy,\n", " component=0,\n", " slice_at=[0, 0, None],\n", " do_plot=True,\n", " disp_name=\"Maxwell1D\",\n", " fit_branches=1,\n", " noise_level=0.5,\n", " extr_order=10,\n", " fit_degree=(1,),\n", ")\n", "\n", "# Extract fitted wave speed\n", "c_light_speed = 1.0 # Theoretical\n", "c_fit = float(coeffs[0][0])\n", "\n", "print(f\"\\nTheoretical wave speed (c): {c_light_speed:.6f}\")\n", "print(f\"Fitted wave speed: {c_fit:.6f}\")\n", "print(f\"Absolute error: {abs(c_fit - c_light_speed):.6f}\")\n", "print(f\"Relative error: {abs(c_fit - c_light_speed) / c_light_speed * 100:.3f}%\")\n", "\n", "# Verify against tolerance\n", "error = xp.abs(coeffs[0][0] - c_light_speed)\n", "tolerance = 0.02\n", "assert error < tolerance, f\"Wave speed error {error:.4f} exceeds tolerance {tolerance}\"\n", "print(f\"\\n✓ Light wave verification passed (error < {tolerance}).\")" ] }, { "cell_type": "markdown", "id": "16", "metadata": {}, "source": [ "### Conclusion\n", "\n", "This tutorial successfully verified the electromagnetic wave dispersion relation in vacuum:\n", "\n", "$$\\omega = c \\, k \\quad \\Rightarrow \\quad v_\\text{phase} = c,$$\n", "\n", "where the wave propagation speed is exactly the speed of light (normalized to 1.0). The FFT-based spectral analysis method provides a robust and quantitative test of the model's ability to correctly implement Maxwell's equations and represent wave dynamics at the grid scale.\n", "\n", "This verification validates:\n", "- Correct implementation of electromagnetic propagators\n", "- Numerical accuracy of the finite-element spatial discretization\n", "- Stability and fidelity of the time-integration algorithm\n", "\n", "Small deviations from perfect $c = 1$ arise from numerical dispersion in the finite-element representation, which decreases with increasing spatial resolution." ] }, { "cell_type": "code", "execution_count": null, "id": "17", "metadata": {}, "outputs": [], "source": [ "# Cleanup temporary simulation folder\n", "if False: # Set to True to enable cleanup\n", " try:\n", " shutil.rmtree(test_folder)\n", " print(f\"Cleaned up {test_folder}\")\n", " except Exception as e:\n", " print(f\"Could not remove {test_folder}: {e}\")" ] }, { "cell_type": "markdown", "id": "18", "metadata": {}, "source": [ "## Verification of Electromagnetic Modes in a Coaxial Waveguide\n", "\n", "### 2D Hollow Cylinder Geometry\n", "\n", "This tutorial demonstrates verification of electromagnetic wave modes in a coaxial waveguide (hollow cylinder) by comparing numerical solutions against analytical expressions involving Bessel functions. We use the `Maxwell` model with a 2D hollow-cylinder domain.\n", "\n", "### Physical Setup\n", "\n", "A coaxial waveguide consists of two concentric conducting cylinders with inner radius $a_1$ and outer radius $a_2$. Electromagnetic modes propagate in the axial ($z$) direction with spatial structure in the radial ($r$) and azimuthal ($\\theta$) directions.\n", "\n", "For a TE-like mode with azimuthal wavenumber $m$, the magnetic field in the $z$-direction has the form\n", "\n", "$$B_z(r, \\theta, t) = [J_m(r) - c_1\\, Y_m(r)]\\cos(m\\theta - t),$$\n", "\n", "while the electric field components transform via cylindrical coordinate relations:\n", "\n", "$$E_r(r, \\theta, t) = -\\frac{m}{r}[J_m(r) - c_1\\, Y_m(r)]\\cos(m\\theta - t),$$\n", "$$E_\\theta(r, \\theta, t) = [J_m'(r) - c_1\\, Y_m'(r)]\\sin(m\\theta - t),$$\n", "\n", "where $J_m$ and $Y_m$ are Bessel functions of the first and second kinds, respectively.\n", "\n", "We verify these modes by comparing numerical simulation against the analytical expressions." ] }, { "cell_type": "code", "execution_count": null, "id": "19", "metadata": {}, "outputs": [], "source": [ "import logging\n", "import os\n", "import shutil\n", "\n", "import cunumpy as xp\n", "from scipy.special import jv, yn # Bessel functions\n", "\n", "from struphy import (\n", " DerhamOptions,\n", " EnvironmentOptions,\n", " Simulation,\n", " Time,\n", " domains,\n", " equils,\n", " grids,\n", " perturbations,\n", ")\n", "from struphy.models import Maxwell\n", "\n", "logger = logging.getLogger(\"struphy\")" ] }, { "cell_type": "markdown", "id": "20", "metadata": {}, "source": [ "### Coaxial Waveguide Geometry\n", "\n", "Define the inner and outer cylinder radii and the waveguide length. We choose radii that match the zero-crossing of Bessel function combinations used in the mode matching." ] }, { "cell_type": "code", "execution_count": null, "id": "21", "metadata": {}, "outputs": [], "source": [ "# Coaxial geometry\n", "a1 = 2.326744 # Inner cylinder radius\n", "a2 = 3.686839 # Outer cylinder radius\n", "Lz = 2.0 # Waveguide length (axial direction)\n", "\n", "# Azimuthal mode number\n", "m = 3\n", "\n", "# Bessel function coefficient (matching inner boundary condition)\n", "# This coefficient is precomputed based on the mode matching\n", "c1_coeff = 0.28\n", "\n", "print(\"Coaxial waveguide geometry:\")\n", "print(f\" Inner radius (a1): {a1:.6f}\")\n", "print(f\" Outer radius (a2): {a2:.6f}\")\n", "print(f\" Axial length (Lz): {Lz:.6f}\")\n", "print(f\" Azimuthal mode (m): {m}\")\n", "print(f\" Bessel coeff (c1): {c1_coeff}\")" ] }, { "cell_type": "markdown", "id": "22", "metadata": {}, "source": [ "### Model and Propagator Configuration\n", "\n", "Create a Maxwell model and configure the implicit time-stepping propagator. Implicit stepping is preferred for stability in 2D domains." ] }, { "cell_type": "code", "execution_count": null, "id": "23", "metadata": {}, "outputs": [], "source": [ "# Model instance\n", "model = Maxwell()\n", "\n", "# Propagator options (implicit for stability)\n", "model.propagators.maxwell.options = model.propagators.maxwell.Options(algo=\"implicit\")\n", "\n", "# Equilibrium (homogeneous slab, used for initialization)\n", "equil = equils.HomogenSlab()" ] }, { "cell_type": "markdown", "id": "24", "metadata": {}, "source": [ "### Domain and Discretization\n", "\n", "Set up a 2D hollow-cylinder domain using curvilinear coordinates. The grid is fine in both radial ($r$) and azimuthal ($\\theta$) directions to resolve the mode structure." ] }, { "cell_type": "code", "execution_count": null, "id": "25", "metadata": {}, "outputs": [], "source": [ "# Domain: hollow cylinder with inner radius a1, outer radius a2\n", "domain = domains.HollowCylinder(a1=a1, a2=a2, Lz=Lz)\n", "\n", "# Grid: 2D in (r, θ) with finer azimuthal resolution\n", "# 32 radial × 64 azimuthal elements for mode representation\n", "grid = grids.TensorProductGrid(num_elements=(32, 64, 1))\n", "\n", "# Derham options for 2D:\n", "# - degree=(3, 3, 1): cubic in radial and azimuthal\n", "# - bcs: Dirichlet on both radial boundaries (conducting cylinders)\n", "derham_opts = DerhamOptions(\n", " degree=(3, 3, 1),\n", " bcs=((\"dirichlet\", \"dirichlet\"), None, None),\n", ")\n", "\n", "print(f\"Domain: HollowCylinder(a1={a1}, a2={a2}, Lz={Lz})\")\n", "print(f\"Grid elements: {grid.num_elements}\")\n", "print(f\"Derham degree: {derham_opts.degree}\")\n", "print(\"Boundary conditions: Dirichlet at r=a1 and r=a2\")" ] }, { "cell_type": "markdown", "id": "26", "metadata": {}, "source": [ "### Initial Conditions: Coaxial Mode Initialization\n", "\n", "Initialize the electric and magnetic fields with the spatial structure of a TE-like coaxial mode. The `CoaxialWaveguide*` perturbation classes provide the analytic mode shapes." ] }, { "cell_type": "code", "execution_count": null, "id": "27", "metadata": {}, "outputs": [], "source": [ "# Electric field components of the coaxial mode\n", "model.em_fields.e_field.add_perturbation(\n", " perturbations.CoaxialWaveguideElectric_r(m=m, a1=a1, a2=a2)\n", ")\n", "model.em_fields.e_field.add_perturbation(\n", " perturbations.CoaxialWaveguideElectric_theta(m=m, a1=a1, a2=a2)\n", ")\n", "\n", "# Magnetic field component of the coaxial mode\n", "model.em_fields.b_field.add_perturbation(\n", " perturbations.CoaxialWaveguideMagnetic(m=m, a1=a1, a2=a2)\n", ")\n", "\n", "print(f\"Coaxial mode (m={m}) initialized with analytic perturbations.\")" ] }, { "cell_type": "markdown", "id": "28", "metadata": {}, "source": [ "### Simulation Setup and Execution\n", "\n", "Configure the simulation environment and run the transient dynamics. We run for a finite time (Tend=10) to allow the mode dynamics to develop." ] }, { "cell_type": "code", "execution_count": null, "id": "29", "metadata": {}, "outputs": [], "source": [ "# Environment and file management\n", "test_folder = os.path.join(os.getcwd(), \"struphy_verification_tests\")\n", "out_folders = os.path.join(test_folder, \"Maxwell\")\n", "env = EnvironmentOptions(out_folders=out_folders, sim_folder=\"coaxial\")\n", "\n", "# Time-stepping\n", "time_opts = Time(dt=0.05, Tend=10.0)\n", "\n", "# Instantiate and run simulation\n", "sim = Simulation(\n", " model=model,\n", " env=env,\n", " time_opts=time_opts,\n", " domain=domain,\n", " equil=equil,\n", " grid=grid,\n", " derham_opts=derham_opts,\n", ")\n", "\n", "print(f\"Running coaxial mode simulation: dt={time_opts.dt}, Tend={time_opts.Tend}\")\n", "sim.run()\n", "print(\"Simulation complete.\")\n", "\n", "# Post-processing (with physical=True to extract physical fields)\n", "sim.pproc(physical=True)\n", "print(\"Post-processing complete.\")" ] }, { "cell_type": "markdown", "id": "30", "metadata": {}, "source": [ "### Diagnostics: Mode Field Verification\n", "\n", "Extract the numerical electric and magnetic fields at the final time step and compare against the analytical Bessel function expressions. We compute relative errors for $E_r$, $E_\\theta$, and $B_z$ components." ] }, { "cell_type": "code", "execution_count": null, "id": "31", "metadata": {}, "outputs": [], "source": [ "# Load plotting data\n", "sim.load_plotting_data()\n", "\n", "# Extract time and field data\n", "t_grid = sim.t_grid\n", "grids_phy = sim.grids_phy\n", "e_field_phy = sim.spline_values.em_fields.e_field_phy.data\n", "b_field_phy = sim.spline_values.em_fields.b_field_phy.data\n", "\n", "# Extract coordinate arrays in the first (r-θ) plane\n", "X = grids_phy[0][:, :, 0] # Radial coordinate (Cartesian x for plotting)\n", "Y = grids_phy[1][:, :, 0] # Azimuthal coordinate (Cartesian y for plotting)\n", "\n", "print(f\"Grid shape: {X.shape}\")\n", "print(f\"Time grid length: {len(t_grid)}\")\n", "print(f\"Final time: t = {t_grid[-1]:.3f}\")" ] }, { "cell_type": "markdown", "id": "32", "metadata": {}, "source": [ "### Analytical Bessel Function Solutions\n", "\n", "Define the analytical field components in terms of Bessel functions of the first and second kinds. These are the reference solutions against which we compare the numerical results." ] }, { "cell_type": "code", "execution_count": null, "id": "33", "metadata": {}, "outputs": [], "source": [ "# Analytical solution using Bessel functions\n", "\n", "def B_z_analytic(X, Y, Z, m, t):\n", " \"\"\"Magnetic field B_z in coaxial cable.\"\"\"\n", " r = (X**2 + Y**2) ** 0.5\n", " theta = xp.arctan2(Y, X)\n", " return (jv(m, r) - c1_coeff * yn(m, r)) * xp.cos(m * theta - t)\n", "\n", "def E_r_analytic(X, Y, Z, m, t):\n", " \"\"\"Electric field E_r in coaxial cable (radial).\"\"\"\n", " r = (X**2 + Y**2) ** 0.5\n", " theta = xp.arctan2(Y, X)\n", " return -m / r * (jv(m, r) - c1_coeff * yn(m, r)) * xp.cos(m * theta - t)\n", "\n", "def E_theta_analytic(X, Y, Z, m, t):\n", " \"\"\"Electric field E_θ in coaxial cable (azimuthal).\"\"\"\n", " r = (X**2 + Y**2) ** 0.5\n", " theta = xp.arctan2(Y, X)\n", " dJ = (m / r * jv(m, r) - jv(m + 1, r))\n", " dY = (m / r * yn(m, r) - yn(m + 1, r))\n", " return (dJ - c1_coeff * dY) * xp.sin(m * theta - t)\n", "\n", "def to_E_r_cartesian(X, Y, E_x, E_y):\n", " \"\"\"Convert Cartesian E_x, E_y to cylindrical E_r.\"\"\"\n", " r = (X**2 + Y**2) ** 0.5\n", " theta = xp.arctan2(Y, X)\n", " return xp.cos(theta) * E_x + xp.sin(theta) * E_y\n", "\n", "def to_E_theta_cartesian(X, Y, E_x, E_y):\n", " \"\"\"Convert Cartesian E_x, E_y to cylindrical E_θ.\"\"\"\n", " r = (X**2 + Y**2) ** 0.5\n", " theta = xp.arctan2(Y, X)\n", " return -xp.sin(theta) * E_x + xp.cos(theta) * E_y\n", "\n", "print(\"Analytical solution functions defined.\")" ] }, { "cell_type": "markdown", "id": "34", "metadata": {}, "source": [ "### Field Component Comparison\n", "\n", "Extract fields at the final time, convert Cartesian components to cylindrical coordinates, and compute relative errors for each field component." ] }, { "cell_type": "code", "execution_count": null, "id": "35", "metadata": {}, "outputs": [], "source": [ "# Get fields at the final time step\n", "t_end = t_grid[-1]\n", "\n", "# Numerical fields (Cartesian components)\n", "Ex_num = e_field_phy[t_end][0][:, :, 0]\n", "Ey_num = e_field_phy[t_end][1][:, :, 0]\n", "Bz_num = b_field_phy[t_end][2][:, :, 0]\n", "\n", "# Analytical fields\n", "Er_analytic = E_r_analytic(X, Y, grids_phy[0], m, t_end)\n", "Etheta_analytic = E_theta_analytic(X, Y, grids_phy[0], m, t_end)\n", "Bz_analytic = B_z_analytic(X, Y, grids_phy[0], m, t_end)\n", "\n", "# Convert numerical Cartesian E-field to cylindrical\n", "Er_num = to_E_r_cartesian(X, Y, Ex_num, Ey_num)\n", "Etheta_num = to_E_theta_cartesian(X, Y, Ex_num, Ey_num)\n", "\n", "# Compute errors\n", "error_Er = xp.max(xp.abs(Er_num - Er_analytic))\n", "error_Etheta = xp.max(xp.abs(Etheta_num - Etheta_analytic))\n", "error_Bz = xp.max(xp.abs(Bz_num - Bz_analytic))\n", "\n", "# Relative errors\n", "rel_err_Er = error_Er / xp.max(xp.abs(Er_analytic))\n", "rel_err_Etheta = error_Etheta / xp.max(xp.abs(Etheta_analytic))\n", "rel_err_Bz = error_Bz / xp.max(xp.abs(Bz_analytic))\n", "\n", "print(\"\\n=== Coaxial Mode Error Analysis ===\")\n", "print(\"\\nRadial Electric Field (E_r):\")\n", "print(f\" Absolute error: {error_Er:.6f}\")\n", "print(f\" Relative error: {rel_err_Er:.6f} ({rel_err_Er*100:.4f}%)\")\n", "\n", "print(\"\\nAzimuthal Electric Field (E_θ):\")\n", "print(f\" Absolute error: {error_Etheta:.6f}\")\n", "print(f\" Relative error: {rel_err_Etheta:.6f} ({rel_err_Etheta*100:.4f}%)\")\n", "\n", "print(\"\\nAxial Magnetic Field (B_z):\")\n", "print(f\" Absolute error: {error_Bz:.6f}\")\n", "print(f\" Relative error: {rel_err_Bz:.6f} ({rel_err_Bz*100:.4f}%)\")" ] }, { "cell_type": "markdown", "id": "36", "metadata": {}, "source": [ "### 2D Slice Visualization of the Final Fields\n", "\n", "Plot 2D cross-sections of the numerical field components at the final time and compare them against the analytical Bessel-function reference. The third column shows the pointwise difference, which makes localized structure and residual discrepancies easy to inspect." ] }, { "cell_type": "code", "execution_count": null, "id": "37", "metadata": {}, "outputs": [], "source": [ "# Plot final-time 2D slices of the coaxial mode fields\n", "field_slices = [\n", " (Er_num, Er_analytic, r\"$E_r$\"),\n", " (Etheta_num, Etheta_analytic, r\"$E_\\theta$\"),\n", " (Bz_num, Bz_analytic, r\"$B_z$\"),\n", "]\n", "\n", "fig, axes = plt.subplots(3, 3, figsize=(15, 12), constrained_layout=True)\n", "\n", "for row, (numerical_field, analytic_field, field_label) in enumerate(field_slices):\n", " numerical_np = np.asarray(numerical_field)\n", " analytic_np = np.asarray(analytic_field)\n", " difference_np = numerical_np - analytic_np\n", "\n", " amplitude_scale = np.max(np.abs(np.stack([numerical_np, analytic_np])))\n", " if amplitude_scale == 0:\n", " amplitude_scale = 1.0\n", "\n", " difference_scale = np.max(np.abs(difference_np))\n", " if difference_scale == 0:\n", " difference_scale = 1.0\n", "\n", " numerical_plot = axes[row, 0].pcolormesh(\n", " X, Y, numerical_np, shading=\"auto\", cmap=\"RdBu_r\",\n", " vmin=-amplitude_scale, vmax=amplitude_scale,\n", " )\n", " analytic_plot = axes[row, 1].pcolormesh(\n", " X, Y, analytic_np, shading=\"auto\", cmap=\"RdBu_r\",\n", " vmin=-amplitude_scale, vmax=amplitude_scale,\n", " )\n", " difference_plot = axes[row, 2].pcolormesh(\n", " X, Y, difference_np, shading=\"auto\", cmap=\"RdBu_r\",\n", " vmin=-difference_scale, vmax=difference_scale,\n", " )\n", "\n", " axes[row, 0].set_title(f\"Numerical {field_label}\")\n", " axes[row, 1].set_title(f\"Analytical {field_label}\")\n", " axes[row, 2].set_title(f\"Difference in {field_label}\")\n", "\n", " fig.colorbar(numerical_plot, ax=axes[row, 0:2], shrink=0.9, label=f\"{field_label} amplitude\")\n", " fig.colorbar(difference_plot, ax=axes[row, 2], shrink=0.9, label=f\"{field_label} error\")\n", "\n", "for ax in axes.flat:\n", " ax.set_aspect(\"equal\")\n", " ax.set_xlabel(\"x\")\n", " ax.set_ylabel(\"y\")\n", "\n", "fig.suptitle(f\"Coaxial mode field slices at t = {t_end:.3f}\", fontsize=16)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "38", "metadata": {}, "source": [ "### Verification: Error Tolerance Check\n", "\n", "Verify that the relative errors for all field components fall below the acceptable tolerance of 0.0021 (0.21%)." ] }, { "cell_type": "code", "execution_count": null, "id": "39", "metadata": {}, "outputs": [], "source": [ "# Tolerance for verification\n", "tolerance = 0.0021\n", "\n", "print(f\"\\n=== Verification Against Tolerance (< {tolerance*100:.2f}%) ===\")\n", "\n", "# Check each component\n", "assertions_passed = []\n", "\n", "try:\n", " assert rel_err_Bz < tolerance, f\"B_z error {rel_err_Bz:.6f} exceeds tolerance\"\n", " print(f\"✓ B_z verification passed ({rel_err_Bz*100:.4f}% < {tolerance*100:.2f}%)\")\n", " assertions_passed.append(True)\n", "except AssertionError as e:\n", " print(f\"✗ {e}\")\n", " assertions_passed.append(False)\n", "\n", "try:\n", " assert rel_err_Etheta < tolerance, f\"E_θ error {rel_err_Etheta:.6f} exceeds tolerance\"\n", " print(f\"✓ E_θ verification passed ({rel_err_Etheta*100:.4f}% < {tolerance*100:.2f}%)\")\n", " assertions_passed.append(True)\n", "except AssertionError as e:\n", " print(f\"✗ {e}\")\n", " assertions_passed.append(False)\n", "\n", "try:\n", " assert rel_err_Er < tolerance, f\"E_r error {rel_err_Er:.6f} exceeds tolerance\"\n", " print(f\"✓ E_r verification passed ({rel_err_Er*100:.4f}% < {tolerance*100:.2f}%)\")\n", " assertions_passed.append(True)\n", "except AssertionError as e:\n", " print(f\"✗ {e}\")\n", " assertions_passed.append(False)\n", "\n", "if all(assertions_passed):\n", " print(\"\\n✓ All coaxial mode components verified successfully.\")\n", "else:\n", " print(\"\\n✗ Some components failed verification.\")" ] }, { "cell_type": "markdown", "id": "40", "metadata": {}, "source": [ "### Conclusion\n", "\n", "This tutorial successfully verified electromagnetic wave modes in a coaxial waveguide by comparing numerical solutions against analytical Bessel function expressions. The verification demonstrates:\n", "\n", "1. **Correct implementation** of the 2D hollow-cylinder domain and curvilinear coordinate transformations\n", "2. **Accurate field representation** of azimuthally-dependent modes in finite-element space\n", "3. **Proper boundary conditions** at the inner and outer conducting cylinders\n", "4. **Stable time-integration** of Maxwell's equations in 2D\n", "\n", "The TE-like coaxial mode ($m=3$) exhibits excellent agreement between numerical and analytical field components, validating the model's ability to represent complex 3D electromagnetic phenomena in reduced geometries." ] }, { "cell_type": "code", "execution_count": null, "id": "41", "metadata": {}, "outputs": [], "source": [ "# Cleanup temporary simulation folder\n", "if False: # Set to True to enable cleanup\n", " try:\n", " shutil.rmtree(test_folder)\n", " print(f\"Cleaned up {test_folder}\")\n", " except Exception as e:\n", " print(f\"Could not remove {test_folder}: {e}\")" ] } ], "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 }