{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4d1c19d9",
   "metadata": {},
   "source": [
    "# Tutorial: 3D Poisson with `volumential`\n",
    "\n",
    "This notebook mirrors the structure of the 2D Poisson tutorial, but runs a true 3D manufactured Poisson solve on a box using `volumential` volume FMM.\n",
    "\n",
    "Compared to the 2D tutorial, this one emphasizes:\n",
    "\n",
    "- larger 3D workloads,\n",
    "- uniform co-refinement in 3D,\n",
    "- richer visualization (orthogonal slices, 3D point-cloud error, interactive isosurfaces).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3dde01f",
   "metadata": {},
   "source": [
    "## Roadmap\n",
    "\n",
    "1. **Problem definition**: manufactured exact solution and forcing `f = -\\Delta u`.\n",
    "2. **Infrastructure setup**: OpenCL context + near-field table cache.\n",
    "3. **Solver module**: reusable `run_poisson3d_case(...)`.\n",
    "4. **Single-resolution walkthrough**: inspect diagnostics and plots.\n",
    "5. **Uniform co-refinement study**: increase grid levels and quantify error/cost trends.\n",
    "6. **Interactive 3D visualization**: isosurfaces for potential and error hot spots.\n",
    "7. **Takeaways**: what to tune next for your own runs.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6611bed3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import annotations\n",
    "\n",
    "from functools import partial\n",
    "from pathlib import Path\n",
    "import os\n",
    "import time\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pyopencl as cl\n",
    "import pyopencl.array  # noqa: F401\n",
    "\n",
    "from pytools.obj_array import new_1d as obj_array_1d\n",
    "\n",
    "import volumential.meshgen as mg\n",
    "from volumential.meshgen import build_geometry_info\n",
    "from volumential.table_manager import NearFieldInteractionTableManager\n",
    "from volumential.nearfield_potential_table import DuffyBuildConfig\n",
    "from volumential.expansion_wrangler_fpnd import (\n",
    "    FPNDExpansionWrangler,\n",
    "    FPNDTreeIndependentDataForWrangler,\n",
    ")\n",
    "from volumential.volume_fmm import drive_volume_fmm, interpolate_volume_potential\n",
    "\n",
    "from sumpy.expansion import DefaultExpansionFactory\n",
    "from sumpy.kernel import LaplaceKernel\n",
    "\n",
    "plt.style.use('seaborn-v0_8-whitegrid')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8adca79",
   "metadata": {},
   "source": [
    "## 1) Problem Definition\n",
    "\n",
    "We use a smooth manufactured exact solution made of two shifted Gaussian bumps:\n",
    "\n",
    "- `u_exact(x, y, z) = g1 + c * g2`\n",
    "- `f(x, y, z) = -\\Delta u_exact`\n",
    "\n",
    "`u_exact` solves the Poisson equation in the whole space, while the volume potential integrates `f` over the box only, so the two differ by the potential of the source outside the box. Both Gaussians use `alpha = 240`: each center is at least 0.41 from the nearest face, which puts both below `exp(-40)` on the boundary. That difference is then negligible, and the errors reported below measure discretization and FMM rather than truncation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2fde747",
   "metadata": {},
   "outputs": [],
   "source": [
    "smoke_mode = os.environ.get('VOLUMENTIAL_EXAMPLE_SMOKE', '').lower() in {'1', 'true', 'yes'}\n",
    "\n",
    "dim = 3\n",
    "bbox_a, bbox_b = -0.5, 0.5\n",
    "table_root_extent = 2.0\n",
    "\n",
    "if smoke_mode:\n",
    "    q_order = 3\n",
    "    default_nlevels = 2\n",
    "    fmm_order = 10\n",
    "    regular_quad_order = 10\n",
    "    radial_quad_order = 35\n",
    "    slice_side = 72\n",
    "    volume_side = 20\n",
    "else:\n",
    "    q_order = 5\n",
    "    default_nlevels = 4\n",
    "    fmm_order = 18\n",
    "    regular_quad_order = 16\n",
    "    radial_quad_order = 80\n",
    "    slice_side = 170\n",
    "    volume_side = 44\n",
    "\n",
    "table_filename = (\n",
    "    'nft_poisson3d_notebook_smoke.sqlite' if smoke_mode\n",
    "    else 'nft_poisson3d_notebook.sqlite'\n",
    ")\n",
    "\n",
    "artifact_dir = Path('poisson3d_notebook_output')\n",
    "artifact_dir.mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "# Both Gaussians are below exp(-40) on the box boundary (see section 1).\n",
    "alpha_1 = 240.0\n",
    "alpha_2 = 240.0\n",
    "coeff_2 = -0.65\n",
    "\n",
    "def _gaussian(alpha, dx, dy, dz):\n",
    "    r2 = dx * dx + dy * dy + dz * dz\n",
    "    return np.exp(-alpha * r2), r2\n",
    "\n",
    "def _minus_laplacian_gaussian(alpha, r2):\n",
    "    return (2 * dim * alpha - 4 * alpha * alpha * r2)\n",
    "\n",
    "def u_exact(x, y, z):\n",
    "    g1, _ = _gaussian(alpha_1, x + 0.08, y - 0.06, z + 0.05)\n",
    "    g2, _ = _gaussian(alpha_2, x - 0.09, y + 0.07, z - 0.08)\n",
    "    return g1 + coeff_2 * g2\n",
    "\n",
    "def rhs_f(x, y, z):\n",
    "    g1, r1 = _gaussian(alpha_1, x + 0.08, y - 0.06, z + 0.05)\n",
    "    g2, r2 = _gaussian(alpha_2, x - 0.09, y + 0.07, z - 0.08)\n",
    "    return _minus_laplacian_gaussian(alpha_1, r1) * g1 + coeff_2 * _minus_laplacian_gaussian(alpha_2, r2) * g2\n",
    "\n",
    "print('smoke_mode:', smoke_mode)\n",
    "print('q_order:', q_order, '| default_nlevels:', default_nlevels, '| fmm_order:', fmm_order)\n",
    "print('table cache:', table_filename)\n",
    "print('artifact dir:', artifact_dir)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9505c1b0",
   "metadata": {},
   "source": [
    "## 2) Infrastructure Setup\n",
    "\n",
    "Build OpenCL context and near-field table once so multiple runs can reuse cache.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a418ea3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ctx = cl.create_some_context()\n",
    "queue = cl.CommandQueue(ctx)\n",
    "\n",
    "tm = NearFieldInteractionTableManager(\n",
    "    table_filename,\n",
    "    root_extent=table_root_extent,\n",
    "    queue=queue,\n",
    ")\n",
    "build_config = DuffyBuildConfig(\n",
    "    radial_rule='tanh-sinh-fast',\n",
    "    regular_quad_order=regular_quad_order,\n",
    "    radial_quad_order=radial_quad_order,\n",
    ")\n",
    "\n",
    "nftable, _ = tm.get_table(\n",
    "    dim,\n",
    "    'Laplace',\n",
    "    q_order,\n",
    "    force_recompute=False,\n",
    "    queue=queue,\n",
    "    build_config=build_config,\n",
    ")\n",
    "\n",
    "print('near-field table ready:', table_filename)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "939f64cd",
   "metadata": {},
   "source": [
    "## 3) Solver Module\n",
    "\n",
    "`run_poisson3d_case(...)` executes one full solve and returns diagnostics.\n",
    "\n",
    "For each run it:\n",
    "\n",
    "1. Builds a 3D tensor-product quadrature mesh.\n",
    "2. Evaluates source values from `f = -\\Delta u_exact`.\n",
    "3. Builds a mesh-aligned source-only tree (`targets=None`) from mesh boxtree geometry so source nodes stay consistent with box extents.\n",
    "4. Builds boxtree traversal and FMM wrangler.\n",
    "5. Executes volume FMM and compares against the exact field at source quadrature nodes.\n",
    "\n",
    "Using mesh-aligned coincident trees avoids unnecessary source-to-target interpolation and removes root-extent drift that otherwise degrades 3D convergence.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51195c0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _to_obj_array(arrays):\n",
    "    return obj_array_1d([cl.array.to_device(queue, np.ascontiguousarray(a)) for a in arrays])\n",
    "\n",
    "def run_poisson3d_case(*, nlevels, fmm_order_case=None, return_fields=False):\n",
    "    if fmm_order_case is None:\n",
    "        fmm_order_case = fmm_order\n",
    "\n",
    "    mesh = mg.MeshGen3D(q_order, nlevels, bbox_a, bbox_b, queue=queue)\n",
    "    q_points, source_weights, tree, trav = build_geometry_info(\n",
    "        ctx,\n",
    "        queue,\n",
    "        dim,\n",
    "        q_order,\n",
    "        mesh,\n",
    "        bbox=np.array([[bbox_a, bbox_b]] * dim, dtype=np.float64),\n",
    "    )\n",
    "\n",
    "    assert tree.sources_are_targets\n",
    "\n",
    "    source_coords_host = np.array([coords.get() for coords in q_points])\n",
    "    source_points_host = np.ascontiguousarray(source_coords_host.T)\n",
    "    source_vals_host = np.ascontiguousarray(\n",
    "        rhs_f(source_coords_host[0], source_coords_host[1], source_coords_host[2])\n",
    "    )\n",
    "    source_vals = cl.array.to_device(queue, source_vals_host)\n",
    "\n",
    "    knl = LaplaceKernel(dim)\n",
    "    expn_factory = DefaultExpansionFactory()\n",
    "    local_expn_class = expn_factory.get_local_expansion_class(knl)\n",
    "    mpole_expn_class = expn_factory.get_multipole_expansion_class(knl)\n",
    "\n",
    "    tree_indep = FPNDTreeIndependentDataForWrangler(\n",
    "        ctx,\n",
    "        partial(mpole_expn_class, knl),\n",
    "        partial(local_expn_class, knl),\n",
    "        [knl],\n",
    "        exclude_self=True,\n",
    "    )\n",
    "\n",
    "    target_to_source = np.arange(tree.ntargets, dtype=np.int32)\n",
    "    wrangler = FPNDExpansionWrangler(\n",
    "        tree_indep=tree_indep,\n",
    "        traversal=trav,\n",
    "        near_field_table=nftable,\n",
    "        dtype=np.float64,\n",
    "        fmm_level_to_order=lambda kernel, kernel_args, tree, lev: fmm_order_case,\n",
    "        quad_order=q_order,\n",
    "        queue=queue,\n",
    "        self_extra_kwargs={'target_to_source': target_to_source},\n",
    "    )\n",
    "\n",
    "    queue.finish()\n",
    "    t0 = time.time()\n",
    "    (pot,) = drive_volume_fmm(\n",
    "        trav,\n",
    "        wrangler,\n",
    "        source_vals * source_weights,\n",
    "        source_vals,\n",
    "        list1_only=False,\n",
    "    )\n",
    "    queue.finish()\n",
    "    elapsed = time.time() - t0\n",
    "\n",
    "    approx_host = pot.get()\n",
    "    exact_host = u_exact(source_points_host[:, 0], source_points_host[:, 1], source_points_host[:, 2])\n",
    "    abs_err = np.abs(approx_host - exact_host)\n",
    "    rel_l2_err = np.linalg.norm(abs_err) / max(np.linalg.norm(exact_host), 1.0e-15)\n",
    "\n",
    "    result = {\n",
    "        'nlevels': int(nlevels),\n",
    "        'h': float((bbox_b - bbox_a) / (2 ** (nlevels - 1))),\n",
    "        'n_sources': int(source_points_host.shape[0]),\n",
    "        'tree_nlevels': int(tree.nlevels),\n",
    "        'fmm_order': int(fmm_order_case),\n",
    "        'elapsed_sec': float(elapsed),\n",
    "        'points_per_sec': float(source_points_host.shape[0] / max(elapsed, 1.0e-15)),\n",
    "        'rel_l2_err': float(rel_l2_err),\n",
    "        'max_abs_err': float(abs_err.max()),\n",
    "        'p99_abs_err': float(np.percentile(abs_err, 99.0)),\n",
    "    }\n",
    "\n",
    "    if return_fields:\n",
    "        result.update({\n",
    "            'source_points_host': source_points_host,\n",
    "            'exact_host': exact_host,\n",
    "            'approx_host': approx_host,\n",
    "            'abs_err': abs_err,\n",
    "            'trav': trav,\n",
    "            'wrangler': wrangler,\n",
    "            'pot': pot,\n",
    "        })\n",
    "\n",
    "    return result\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb1466a1",
   "metadata": {},
   "source": [
    "## 4) Single-Resolution Walkthrough\n",
    "\n",
    "Run one representative case, print summary metrics, and render 2D/3D diagnostics.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aecd3b83",
   "metadata": {},
   "outputs": [],
   "source": [
    "single = run_poisson3d_case(nlevels=default_nlevels, return_fields=True)\n",
    "\n",
    "for key in ['nlevels', 'n_sources', 'fmm_order', 'elapsed_sec', 'points_per_sec', 'rel_l2_err', 'max_abs_err', 'p99_abs_err']:\n",
    "    print(f\"{key:>14}: {single[key]}\")\n",
    "\n",
    "def _plane_coords(plane_name, line, level=0.0):\n",
    "    uu, vv = np.meshgrid(line, line, indexing='xy')\n",
    "    if plane_name == 'xy':\n",
    "        xx, yy, zz = uu, vv, np.full_like(uu, level)\n",
    "        plot_u, plot_v = xx, yy\n",
    "        labels = ('x', 'y')\n",
    "    elif plane_name == 'xz':\n",
    "        xx, yy, zz = uu, np.full_like(uu, level), vv\n",
    "        plot_u, plot_v = xx, zz\n",
    "        labels = ('x', 'z')\n",
    "    elif plane_name == 'yz':\n",
    "        xx, yy, zz = np.full_like(uu, level), uu, vv\n",
    "        plot_u, plot_v = yy, zz\n",
    "        labels = ('y', 'z')\n",
    "    else:\n",
    "        raise ValueError(plane_name)\n",
    "    return labels, plot_u, plot_v, xx, yy, zz\n",
    "\n",
    "line = np.linspace(bbox_a, bbox_b, slice_side)\n",
    "slices = {}\n",
    "for plane_name in ('xy', 'xz', 'yz'):\n",
    "    labels, plot_u, plot_v, xx, yy, zz = _plane_coords(plane_name, line)\n",
    "    targets = _to_obj_array([xx.ravel(), yy.ravel(), zz.ravel()])\n",
    "    approx = interpolate_volume_potential(targets, single['trav'], single['wrangler'], single['pot']).get().reshape(xx.shape)\n",
    "    exact = u_exact(xx, yy, zz)\n",
    "    slices[plane_name] = {\n",
    "        'labels': labels,\n",
    "        'u': plot_u,\n",
    "        'v': plot_v,\n",
    "        'exact': exact,\n",
    "        'approx': approx,\n",
    "        'abs_err': np.abs(approx - exact),\n",
    "    }\n",
    "\n",
    "fig, axes = plt.subplots(3, 3, figsize=(14, 12), constrained_layout=True)\n",
    "for irow, plane_name in enumerate(('xy', 'xz', 'yz')):\n",
    "    data = slices[plane_name]\n",
    "    uu = data['u']\n",
    "    vv = data['v']\n",
    "    exact = data['exact']\n",
    "    approx = data['approx']\n",
    "    logerr = np.log10(data['abs_err'] + 1.0e-18)\n",
    "\n",
    "    extent = (uu.min(), uu.max(), vv.min(), vv.max())\n",
    "    vmin = min(exact.min(), approx.min())\n",
    "    vmax = max(exact.max(), approx.max())\n",
    "\n",
    "    im0 = axes[irow, 0].imshow(exact, origin='lower', extent=extent, cmap='viridis', vmin=vmin, vmax=vmax)\n",
    "    axes[irow, 1].imshow(approx, origin='lower', extent=extent, cmap='viridis', vmin=vmin, vmax=vmax)\n",
    "    im2 = axes[irow, 2].imshow(logerr, origin='lower', extent=extent, cmap='magma')\n",
    "\n",
    "    axes[irow, 0].set_title(f'Exact ({plane_name})')\n",
    "    axes[irow, 1].set_title(f'FMM ({plane_name})')\n",
    "    axes[irow, 2].set_title(f'log10 abs err ({plane_name})')\n",
    "\n",
    "    labels = data['labels']\n",
    "    for icol in range(3):\n",
    "        axes[irow, icol].set_xlabel(labels[0])\n",
    "    axes[irow, 0].set_ylabel(labels[1])\n",
    "\n",
    "    fig.colorbar(im0, ax=[axes[irow, 0], axes[irow, 1]], shrink=0.8)\n",
    "    fig.colorbar(im2, ax=axes[irow, 2], shrink=0.8)\n",
    "\n",
    "fig.suptitle('Poisson 3D slices: exact vs FMM', fontsize=13)\n",
    "slice_path = artifact_dir / 'single_slices.png'\n",
    "fig.savefig(slice_path, dpi=220)\n",
    "print('wrote', slice_path)\n",
    "plt.show()\n",
    "\n",
    "max_plot_points = 50000\n",
    "npts = single['n_sources']\n",
    "stride = max(1, int(np.ceil(npts / max_plot_points)))\n",
    "pts = single['source_points_host'][::stride]\n",
    "err = single['abs_err'][::stride]\n",
    "\n",
    "fig = plt.figure(figsize=(10, 7))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "sc = ax.scatter(\n",
    "    pts[:, 0], pts[:, 1], pts[:, 2],\n",
    "    c=np.log10(err + 1.0e-18),\n",
    "    cmap='inferno',\n",
    "    s=2,\n",
    "    alpha=0.5,\n",
    "    linewidths=0,\n",
    ")\n",
    "ax.set_title('3D quadrature-node error cloud (log10 abs err)')\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('y')\n",
    "ax.set_zlabel('z')\n",
    "fig.colorbar(sc, ax=ax, shrink=0.72)\n",
    "cloud_path = artifact_dir / 'single_error_cloud.png'\n",
    "fig.savefig(cloud_path, dpi=220)\n",
    "print('wrote', cloud_path)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3437080",
   "metadata": {},
   "source": [
    "## 5) Uniform Co-Refinement Study\n",
    "\n",
    "Refine volume resolution (`nlevels`) while keeping the physics fixed.\n",
    "\n",
    "This reports both accuracy and runtime scaling.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a58c7378",
   "metadata": {},
   "outputs": [],
   "source": [
    "study_levels = [2, 3] if smoke_mode else [3, 4, 5]\n",
    "study = []\n",
    "for lev in study_levels:\n",
    "    row = run_poisson3d_case(nlevels=lev)\n",
    "    study.append(row)\n",
    "    print(\n",
    "        f\"level={row['nlevels']} | sources={row['n_sources']} | elapsed={row['elapsed_sec']:.2f}s | rel_l2={row['rel_l2_err']:.3e} | max_abs={row['max_abs_err']:.3e}\"\n",
    "    )\n",
    "\n",
    "hs = np.array([row['h'] for row in study])\n",
    "rel_l2 = np.array([row['rel_l2_err'] for row in study])\n",
    "max_abs = np.array([row['max_abs_err'] for row in study])\n",
    "n_sources = np.array([row['n_sources'] for row in study])\n",
    "elapsed = np.array([row['elapsed_sec'] for row in study])\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(16, 4.5), constrained_layout=True)\n",
    "axes[0].loglog(hs, rel_l2, 'o-', label='relative L2')\n",
    "axes[0].loglog(hs, max_abs, 's-', label='max abs')\n",
    "axes[0].invert_xaxis()\n",
    "axes[0].set_xlabel('h')\n",
    "axes[0].set_ylabel('error')\n",
    "axes[0].set_title('Error vs mesh spacing')\n",
    "axes[0].legend()\n",
    "\n",
    "axes[1].plot(n_sources, elapsed, 'o-')\n",
    "axes[1].set_xlabel('number of sources')\n",
    "axes[1].set_ylabel('runtime (s)')\n",
    "axes[1].set_title('Runtime scaling')\n",
    "\n",
    "axes[2].plot(n_sources, np.array([row['points_per_sec'] for row in study]), 'o-')\n",
    "axes[2].set_xlabel('number of sources')\n",
    "axes[2].set_ylabel('points / second')\n",
    "axes[2].set_title('Throughput')\n",
    "\n",
    "refine_path = artifact_dir / 'refinement_study.png'\n",
    "fig.savefig(refine_path, dpi=220)\n",
    "print('wrote', refine_path)\n",
    "plt.show()\n",
    "\n",
    "if len(study) >= 2:\n",
    "    slope_rel_l2 = np.polyfit(np.log(hs), np.log(rel_l2), 1)[0]\n",
    "    slope_max_abs = np.polyfit(np.log(hs), np.log(max_abs), 1)[0]\n",
    "    print(f'estimated order from rel_l2: {slope_rel_l2:.3f}')\n",
    "    print(f'estimated order from max_abs: {slope_max_abs:.3f}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a5a5a49",
   "metadata": {},
   "source": [
    "## 6) Interactive 3D Isosurfaces\n",
    "\n",
    "Sample a regular 3D grid and render:\n",
    "\n",
    "- potential isosurfaces (Viridis),\n",
    "- error hot-spot isosurfaces (Reds).\n",
    "\n",
    "If `plotly` is not available, this cell prints a skip message.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b3b33ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "vol_line = np.linspace(bbox_a, bbox_b, volume_side)\n",
    "xx, yy, zz = np.meshgrid(vol_line, vol_line, vol_line, indexing='ij')\n",
    "targets = _to_obj_array([xx.ravel(), yy.ravel(), zz.ravel()])\n",
    "\n",
    "approx_vol = interpolate_volume_potential(\n",
    "    targets, single['trav'], single['wrangler'], single['pot']\n",
    ").get().reshape(xx.shape)\n",
    "exact_vol = u_exact(xx, yy, zz)\n",
    "abs_err = np.abs(approx_vol - exact_vol)\n",
    "\n",
    "print('volume grid shape:', approx_vol.shape, '| max abs err:', float(abs_err.max()))\n",
    "\n",
    "volume_npz = artifact_dir / 'volume_samples.npz'\n",
    "np.savez_compressed(\n",
    "    volume_npz,\n",
    "    line=vol_line,\n",
    "    approx=approx_vol,\n",
    "    exact=exact_vol,\n",
    "    abs_err=abs_err,\n",
    ")\n",
    "print('wrote', volume_npz)\n",
    "\n",
    "try:\n",
    "    import plotly.graph_objects as go\n",
    "\n",
    "    pot_min = float(np.percentile(approx_vol, 35))\n",
    "    pot_max = float(np.percentile(approx_vol, 98))\n",
    "    err_min = float(np.percentile(abs_err, 92))\n",
    "    err_max = float(np.percentile(abs_err, 99.8))\n",
    "\n",
    "    fig = go.Figure()\n",
    "    fig.add_trace(\n",
    "        go.Isosurface(\n",
    "            x=xx.ravel(), y=yy.ravel(), z=zz.ravel(),\n",
    "            value=approx_vol.ravel(),\n",
    "            isomin=pot_min, isomax=pot_max,\n",
    "            surface_count=8,\n",
    "            colorscale='Viridis',\n",
    "            opacity=0.62,\n",
    "            caps={'x_show': False, 'y_show': False, 'z_show': False},\n",
    "            name='Potential',\n",
    "        )\n",
    "    )\n",
    "    fig.add_trace(\n",
    "        go.Isosurface(\n",
    "            x=xx.ravel(), y=yy.ravel(), z=zz.ravel(),\n",
    "            value=abs_err.ravel(),\n",
    "            isomin=err_min, isomax=err_max,\n",
    "            surface_count=3,\n",
    "            colorscale='Reds',\n",
    "            opacity=0.18,\n",
    "            showscale=False,\n",
    "            caps={'x_show': False, 'y_show': False, 'z_show': False},\n",
    "            name='Error hot spots',\n",
    "        )\n",
    "    )\n",
    "\n",
    "    fig.update_layout(\n",
    "        title='Poisson 3D: potential + error hot spots',\n",
    "        template='plotly_white',\n",
    "        scene={\n",
    "            'xaxis_title': 'x',\n",
    "            'yaxis_title': 'y',\n",
    "            'zaxis_title': 'z',\n",
    "            'aspectmode': 'cube',\n",
    "        },\n",
    "    )\n",
    "    iso_html = artifact_dir / 'single_isosurface.html'\n",
    "    fig.write_html(str(iso_html), include_plotlyjs='cdn')\n",
    "    print('wrote', iso_html)\n",
    "    fig.show()\n",
    "except ImportError:\n",
    "    print('plotly is not installed; skipping interactive isosurface')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8a32143",
   "metadata": {},
   "source": [
    "## 7) Takeaways\n",
    "\n",
    "- 3D volume FMM gives a practical path to Poisson solves with manufactured verification.\n",
    "- Orthogonal slices are still the fastest sanity-check for 3D fields.\n",
    "- Point-cloud and isosurface views are useful to locate localized error hot spots.\n",
    "- For larger production runs, tune `q_order`, `nlevels`, and FMM order jointly based on your error/runtime target.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 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.13.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}