#!/usr/bin/env python3 """ Scientific ML Studio -- reference code (written by hand, not generated by the Studio) Licence : Creative Commons Attribution-NonCommercial 4.0 International (CC BY-NC 4.0) https://creativecommons.org/licenses/by-nc/4.0/ (c) 2026 Scientific ML Studio. All rights are reserved to the Scientific ML Studio platform except the permissions granted by the licence above. You may share and adapt this code for NON-COMMERCIAL purposes, provided you credit "Scientific ML Studio" with a link to https://scimlstudio.com/, and indicate any changes you make. Commercial use, paid client work and production deployment need prior written permission, which you can ask for through the site's Contact page. ------------------------------------------------------------------------------- beam_compare.py - check a Scientific ML Studio PINN for a simply supported Euler-Bernoulli beam against the exact solution and a finite element solution. The problem (non-dimensional, x in [0, 1]): EI w'''' = q uniform load q, length L, bending stiffness EI w(0) = w(1) = 0 pinned / roller supports (no deflection) w''(0) = w''(1) = 0 no bending moment at the supports Because the beam is statically determinate, statics gives the bending moment M(x) = q x (L - x) / 2, and the moment-curvature law EI w'' = -M(x) turns the fourth-order problem into a second-order one. That second-order problem is what the Studio notebook solves (a Poisson block). The deflection is scaled by its exact maximum, w_max = 5 q L^4 / (384 EI), so the answer peaks at 1: u = w / w_max, u'' = -38.4 x (1 - x), u(0) = u(1) = 0 Studio source: -38.4*x*(1-x) This script then checks the trained network three ways: 1. against the exact solution u = 3.2 x (1 - 2x^2 + x^3), 2. against a finite element solution of the full fourth-order equation (cubic Hermite beam elements, written here in NumPy/SciPy), and 3. on quantities the network was never trained on directly: the slope, and the shear force, which needs the THIRD derivative of the network. Usage ----- # FEM against exact only (no PyTorch needed): python beam_compare.py --verify # one or more trained runs of the Studio script: python beam_compare.py --script simply-supported-beam.py --checkpoint run1.pt run2.pt run3.pt # optionally also the table the Studio script wrote: python beam_compare.py --script simply-supported-beam.py --checkpoint run1.pt run2.pt run3.pt --results results.csv (results.csv is compared with the LAST checkpoint listed, since every run rewrites it) Outputs go to ./beam_results/ : figures (PNG), summary.md, summary.json and comparison.csv. Needs numpy, scipy and matplotlib; torch only for the PINN part. Licence: CC BY-NC 4.0, Scientific ML Studio. """ from __future__ import annotations import argparse import csv import importlib.util import json import math import os import sys import time from pathlib import Path import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as spla # ---------------------------------------------------------------------------- # Problem constants # ---------------------------------------------------------------------------- W_MAX_BAR = 5.0 / 384.0 # max of w*EI/(q L^4) for this beam SCALE = 1.0 / W_MAX_BAR # 76.8: u = w / w_max = SCALE * w_bar N_GRID = 2001 # uniform evaluation grid on [0, 1] N_FRESH = 5000 # fresh random points for the residual check FRESH_SEED = 12345 # different from any training seed # Colours (site coral darkened for contrast, site blue, ink for the exact line). C_PINN = "#E5533D" C_FEM = "#187EE4" C_EXACT = "#1F2328" C_MUTED = "#6B7280" C_GRID = "#E5E7EB" BLUE_RAMP = ["#9CC7F4", "#5DA3EC", "#187EE4", "#0B4F96"] # w, slope, M, V LINESTYLES = ["-", "--", "-.", ":", (0, (5, 1)), (0, (3, 1, 1, 1))] # ---------------------------------------------------------------------------- # Exact solution (scaled: u = w / w_max) # ---------------------------------------------------------------------------- def exact(x: np.ndarray, scale: float = SCALE) -> dict: """u and its first three derivatives, plus slope, moment and shear. In the scaled units the moment is M = -u'' and the shear is V = dM/dx = -u'''. Multiply by w_max * EI etc. to get dimensional values (see summary.md). """ wb = (x - 2 * x**3 + x**4) / 24.0 wb1 = (1 - 6 * x**2 + 4 * x**3) / 24.0 wb2 = -x * (1 - x) / 2.0 wb3 = -(1 - 2 * x) / 2.0 u, u1, u2, u3 = (scale * a for a in (wb, wb1, wb2, wb3)) return {"u": u, "slope": u1, "u2": u2, "u3": u3, "M": -u2, "V": -u3} def source(x: np.ndarray, scale: float = SCALE) -> np.ndarray: """The f in u'' = f that the Studio's Poisson block should hold.""" return -scale * x * (1 - x) / 2.0 # ---------------------------------------------------------------------------- # Finite elements: cubic Hermite beam elements for EI w'''' = q # ---------------------------------------------------------------------------- class BeamFEM: """Euler-Bernoulli beam, N equal elements on [0, 1], EI = q = 1. Unknowns per node: deflection w and slope w'. Supports fix w at both ends; the moment-free condition w'' = 0 is natural and needs no constraint. """ def __init__(self, n_el: int): self.n_el = n_el self.h = 1.0 / n_el self.nodes = np.linspace(0.0, 1.0, n_el + 1) self.d = self._solve() def _solve(self) -> np.ndarray: n, h = self.n_el, self.h ndof = 2 * (n + 1) ke = np.array([[12, 6 * h, -12, 6 * h], [6 * h, 4 * h * h, -6 * h, 2 * h * h], [-12, -6 * h, 12, -6 * h], [6 * h, 2 * h * h, -6 * h, 4 * h * h]]) / h**3 fe = h * np.array([0.5, h / 12.0, 0.5, -h / 12.0]) # consistent load, q = 1 rows, cols, vals = [], [], [] f = np.zeros(ndof) for e in range(n): dofs = [2 * e, 2 * e + 1, 2 * e + 2, 2 * e + 3] for a in range(4): f[dofs[a]] += fe[a] for b in range(4): rows.append(dofs[a]); cols.append(dofs[b]); vals.append(ke[a, b]) K = sp.csr_matrix((vals, (rows, cols)), shape=(ndof, ndof)) fixed = [0, 2 * n] # w(0) = w(1) = 0 free = np.setdiff1d(np.arange(ndof), fixed) d = np.zeros(ndof) d[free] = spla.spsolve(K[free][:, free].tocsc(), f[free]) return d def evaluate(self, x: np.ndarray, scale: float = SCALE) -> dict: """w and derivatives at any x, in the same scaled units as exact().""" h = self.h e = np.clip((x / h).astype(int), 0, self.n_el - 1) s = (x - e * h) / h w1, t1 = self.d[2 * e], self.d[2 * e + 1] w2, t2 = self.d[2 * e + 2], self.d[2 * e + 3] N = [1 - 3 * s**2 + 2 * s**3, h * (s - 2 * s**2 + s**3), 3 * s**2 - 2 * s**3, h * (-s**2 + s**3)] dN = [(-6 * s + 6 * s**2) / h, 1 - 4 * s + 3 * s**2, (6 * s - 6 * s**2) / h, -2 * s + 3 * s**2] ddN = [(-6 + 12 * s) / h**2, (-4 + 6 * s) / h, (6 - 12 * s) / h**2, (-2 + 6 * s) / h] dddN = [np.full_like(s, 12 / h**3), np.full_like(s, 6 / h**2), np.full_like(s, -12 / h**3), np.full_like(s, 6 / h**2)] dof = [w1, t1, w2, t2] u, u1, u2, u3 = (scale * sum(Ni * di for Ni, di in zip(NN, dof)) for NN in (N, dN, ddN, dddN)) return {"u": u, "slope": u1, "u2": u2, "u3": u3, "M": -u2, "V": -u3} def nodal_error(self, scale: float = SCALE) -> float: """Largest error of the nodal deflections against the exact solution.""" return float(np.max(np.abs(scale * self.d[0::2] - exact(self.nodes, scale)["u"]))) # ---------------------------------------------------------------------------- # Error measures # ---------------------------------------------------------------------------- def rel_l2(a: np.ndarray, ref: np.ndarray) -> float: return float(np.linalg.norm(a - ref) / np.linalg.norm(ref)) def errors_vs_exact(x: np.ndarray, sol: dict, ex: dict) -> dict: out = {} for k in ("u", "slope", "M", "V"): err = np.abs(sol[k] - ex[k]) i = int(np.argmax(err)) out[k] = {"rel_l2": rel_l2(sol[k], ex[k]), "max_abs": float(err[i]), "x_at_max": float(x[i]), "rel_max": float(err[i] / np.max(np.abs(ex[k])))} return out def fem_convergence(scale: float = SCALE, ns=(2, 4, 8, 16, 32, 64)) -> list: x = np.linspace(0, 1, 8001) ex = exact(x, scale) rows = [] for n in ns: fem = BeamFEM(n) e = errors_vs_exact(x, fem.evaluate(x, scale), ex) rows.append({"elements": n, "h": 1.0 / n, "nodal_max_err_u": fem.nodal_error(scale), **{f"rel_l2_{k}": e[k]["rel_l2"] for k in e}}) for k in ("u", "slope", "M", "V"): # observed order between successive meshes for a, b in zip(rows[:-1], rows[1:]): ea, eb = a[f"rel_l2_{k}"], b[f"rel_l2_{k}"] b[f"order_{k}"] = (math.log(ea / eb) / math.log(a["h"] / b["h"]) if ea > 1e-13 and eb > 1e-13 else float("nan")) return rows # ---------------------------------------------------------------------------- # The Studio's network # ---------------------------------------------------------------------------- def load_studio_module(script: Path): """Import the generated script. Importing it does not train (it trains only under `if __name__ == "__main__"`), it only builds the network `net`.""" name = "studio_" + "".join(ch if ch.isalnum() else "_" for ch in script.stem) spec = importlib.util.spec_from_file_location(name, script) mod = importlib.util.module_from_spec(spec) sys.modules[name] = mod spec.loader.exec_module(mod) return mod def find_net(mod): import torch net = getattr(mod, "net", None) if isinstance(net, torch.nn.Module): return net for v in vars(mod).values(): if isinstance(v, torch.nn.Module): return v raise RuntimeError("No network found in the script (expected a variable called `net`).") def load_checkpoint(net, path: Path): import torch dev = next(net.parameters()).device try: ckpt = torch.load(path, map_location=dev, weights_only=True) except Exception: ckpt = torch.load(path, map_location=dev, weights_only=False) # your own file state = ckpt["model"] if isinstance(ckpt, dict) and "model" in ckpt else ckpt net.load_state_dict(state) net.eval() return ckpt if isinstance(ckpt, dict) else {} def network_output(mod, net, pts): """u at the points. Uses the script's own `field(pts, 0)` if it has one, so any output transform the generator adds is respected; otherwise net(pts)[:, 0].""" field = getattr(mod, "field", None) if callable(field): try: return field(pts, 0) except Exception: pass return net(pts)[:, 0:1] def pinn_derivatives(mod, net, x: np.ndarray) -> dict: """u, u', u'', u''' of the trained network by automatic differentiation.""" import torch p = next(net.parameters()) pts = torch.tensor(x.reshape(-1, 1), dtype=p.dtype, device=p.device, requires_grad=True) u = network_output(mod, net, pts) ones = torch.ones_like(u) d1 = torch.autograd.grad(u, pts, ones, create_graph=True)[0] d2 = torch.autograd.grad(d1, pts, ones, create_graph=True)[0] d3 = torch.autograd.grad(d2, pts, ones, create_graph=False)[0] u, d1, d2, d3 = (t.detach().cpu().numpy().ravel().astype(float) for t in (u, d1, d2, d3)) return {"u": u, "slope": d1, "u2": d2, "u3": d3, "M": -d2, "V": -d3} def studio_residual(mod, net, x: np.ndarray): """The script's own interior residual (u_xx - f), if it exposes residual_pde1.""" import torch fn = getattr(mod, "residual_pde1", None) if not callable(fn): return None p = next(net.parameters()) pts = torch.tensor(x.reshape(-1, 1), dtype=p.dtype, device=p.device, requires_grad=True) try: r = fn(pts) except Exception: return None return r.detach().cpu().numpy().ravel().astype(float) def read_results_csv(path: Path): """The Studio's results.csv: header row starting with x, then x,u rows, then a blank line and a summary block of lines starting with '#'.""" xs, us = [], [] with open(path, newline="") as f: for row in csv.reader(f): if not row or row[0].strip().startswith("#") or row[0].strip() == "x": continue try: xs.append(float(row[0])); us.append(float(row[1])) except (ValueError, IndexError): continue return np.array(xs), np.array(us) # ---------------------------------------------------------------------------- # Figures # ---------------------------------------------------------------------------- def setup_mpl(): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({ "figure.dpi": 150, "savefig.dpi": 150, "savefig.bbox": "tight", "font.size": 10, "axes.titlesize": 11, "axes.labelsize": 10, "axes.edgecolor": C_MUTED, "axes.labelcolor": C_EXACT, "axes.linewidth": 0.8, "xtick.color": C_MUTED, "ytick.color": C_MUTED, "axes.grid": True, "grid.color": C_GRID, "grid.linewidth": 0.6, "axes.spines.top": False, "axes.spines.right": False, "legend.frameon": False, "lines.linewidth": 1.8, }) return plt def fig_sketch(out: Path): plt = setup_mpl() fig, ax = plt.subplots(figsize=(7, 2.4)) ax.axis("off") ax.plot([0, 1], [0, 0], color=C_EXACT, lw=5, solid_capstyle="butt") # pin (left) and roller (right) ax.fill([0, -0.04, 0.04], [0, -0.12, -0.12], color="white", ec=C_EXACT, lw=1.4) ax.plot([-0.07, 0.07], [-0.13, -0.13], color=C_EXACT, lw=1.4) ax.fill([1, 0.96, 1.04], [0, -0.10, -0.10], color="white", ec=C_EXACT, lw=1.4) for cx in (0.975, 1.025): ax.add_patch(plt.Circle((cx, -0.122), 0.018, fill=False, ec=C_EXACT, lw=1.2)) ax.plot([0.93, 1.07], [-0.15, -0.15], color=C_EXACT, lw=1.4) # uniform load for xa in np.linspace(0, 1, 11): ax.annotate("", xy=(xa, 0.035), xytext=(xa, 0.28), arrowprops=dict(arrowstyle="-|>", color=C_PINN, lw=1.3)) ax.plot([0, 1], [0.28, 0.28], color=C_PINN, lw=1.3) ax.text(0.5, 0.33, "uniform load q", ha="center", va="bottom", color=C_EXACT) ax.text(0.5, -0.07, "EI", ha="center", va="top", color=C_EXACT) # dimension line ax.annotate("", xy=(0, -0.27), xytext=(1, -0.27), arrowprops=dict(arrowstyle="<->", color=C_MUTED, lw=1)) ax.text(0.5, -0.29, "L", ha="center", va="top", color=C_MUTED) ax.annotate("", xy=(0.17, -0.05), xytext=(0.0, -0.05), arrowprops=dict(arrowstyle="-|>", color=C_MUTED, lw=1)) ax.text(0.18, -0.05, "x", va="center", color=C_MUTED) ax.set_xlim(-0.12, 1.12); ax.set_ylim(-0.38, 0.42); ax.set_aspect("auto") fig.savefig(out / "beam_sketch.png"); plt.close(fig) def fig_deflection(out, x, ex, fem, fem_n, runs): plt = setup_mpl() fig, (a1, a2) = plt.subplots(2, 1, figsize=(7, 6.2), sharex=True, gridspec_kw={"height_ratios": [1.3, 1]}) a1.plot(x, ex["u"], color=C_EXACT, ls="--", lw=1.4, label="exact", zorder=3) for i, r in enumerate(runs): a1.plot(x, r["sol"]["u"], color=C_PINN, ls=LINESTYLES[i % 6], lw=1.8, label=f"PINN, {r['label']}", zorder=2) a1.plot(fem.nodes, SCALE * fem.d[0::2], "o", color=C_FEM, ms=6, mec="white", mew=1.2, label=f"FEM nodes ({fem_n} elements)", zorder=4) a1.invert_yaxis() a1.set_ylabel("deflection w / w_max (down)") a1.set_title("Deflection of the simply supported beam") a1.legend(loc="upper center", ncol=2) for i, r in enumerate(runs): a2.semilogy(x, np.maximum(np.abs(r["sol"]["u"] - ex["u"]), 1e-12), color=C_PINN, ls=LINESTYLES[i % 6], label=f"PINN, {r['label']}") a2.semilogy(x, np.maximum(np.abs(fem.evaluate(x, SCALE)["u"] - ex["u"]), 1e-12), color=C_FEM, label=f"FEM ({fem_n} elements)") a2.set_ylim(1e-12, None) a2.set_ylabel("|error| in w / w_max") a2.set_xlabel("x / L") a2.legend(loc="lower center", ncol=2) fig.savefig(out / "deflection.png"); plt.close(fig) def fig_derived(out, x, ex, fem_sol, fem_n, runs): plt = setup_mpl() rows = [("slope", "slope dw/dx (scaled)"), ("M", "bending moment M (scaled)"), ("V", "shear force V (scaled)")] fig, axes = plt.subplots(3, 1, figsize=(7, 8), sharex=True) for ax, (k, lab) in zip(axes, rows): ax.plot(x, fem_sol[k], color=C_FEM, lw=1.2, label=f"FEM ({fem_n} elements)") for i, r in enumerate(runs): ax.plot(x, r["sol"][k], color=C_PINN, ls=LINESTYLES[i % 6], label=f"PINN, {r['label']}") ax.plot(x, ex[k], color=C_EXACT, ls="--", lw=1.4, label="exact") ax.set_ylabel(lab) lo, hi = float(np.min(ex[k])), float(np.max(ex[k])) pad = 0.25 * (hi - lo) ax.set_ylim(lo - pad, hi + pad) if any(np.any((r["sol"][k] < lo - pad) | (r["sol"][k] > hi + pad)) for r in runs): ax.text(0.01, 0.03, "PINN curve leaves the plotted range near an end (clipped)", transform=ax.transAxes, fontsize=8, color=C_MUTED) axes[0].set_title("Slope, moment and shear: first, second and third derivative") axes[0].legend(loc="lower left", ncol=2, fontsize=8) axes[-1].set_xlabel("x / L") fig.savefig(out / "slope_moment_shear.png"); plt.close(fig) def fig_derived_errors(out, x, ex, fem_sol, fem_n, runs): plt = setup_mpl() names = [("u", "deflection w (0th derivative)"), ("slope", "slope (1st)"), ("M", "moment (2nd)"), ("V", "shear (3rd)")] fig, axes = plt.subplots(2, 2, figsize=(9, 6.2), sharex=True) for ax, (k, title) in zip(axes.ravel(), names): scale_k = np.max(np.abs(ex[k])) for i, r in enumerate(runs): ax.semilogy(x, np.maximum(np.abs(r["sol"][k] - ex[k]) / scale_k, 1e-12), color=C_PINN, ls=LINESTYLES[i % 6], label=f"PINN, {r['label']}") ax.semilogy(x, np.maximum(np.abs(fem_sol[k] - ex[k]) / scale_k, 1e-12), color=C_FEM, label=f"FEM ({fem_n} el.)") ax.set_title(title) ax.set_ylim(1e-12, None) for ax in axes[:, 0]: ax.set_ylabel("|error| / max |exact|") for ax in axes[-1]: ax.set_xlabel("x / L") axes[0, 0].legend(loc="lower center", fontsize=8) fig.suptitle("Pointwise error against the exact solution", y=1.0) fig.tight_layout() fig.savefig(out / "derivative_errors.png"); plt.close(fig) def fig_convergence(out, conv): plt = setup_mpl() fig, ax = plt.subplots(figsize=(6.4, 4.6)) h = np.array([r["h"] for r in conv]) labels = {"u": "deflection", "slope": "slope", "M": "moment", "V": "shear"} markers = {"u": "o", "slope": "s", "M": "^", "V": "D"} for c, k in zip(BLUE_RAMP, ("u", "slope", "M", "V")): e = np.array([r[f"rel_l2_{k}"] for r in conv]) order = conv[-1].get(f"order_{k}", float("nan")) ax.loglog(h, e, color=c, marker=markers[k], ms=6, mec="white", mew=1, label=f"{labels[k]} (order {order:.1f})") ax.set_xticks(h) ax.set_xticklabels([f"1/{r['elements']}" for r in conv]) ax.minorticks_off() ax.set_xlabel("element size h / L") ax.set_ylabel("relative L2 error") ax.set_title("FEM (cubic Hermite elements): mesh convergence") ax.legend(loc="lower right") fig.savefig(out / "fem_convergence.png"); plt.close(fig) def fig_residual(out, runs): plt = setup_mpl() fig, ax = plt.subplots(figsize=(7, 3.4)) for i, r in enumerate(runs): o = np.argsort(r["fresh_x"]) ax.plot(r["fresh_x"][o], r["fresh_res"][o], color=C_PINN, ls=LINESTYLES[i % 6], lw=1.2, label=f"{r['label']} (RMS {r['metrics']['residual_rms_fresh']:.1e})") ax.axhline(0, color=C_MUTED, lw=0.8) ax.set_xlabel("x / L"); ax.set_ylabel("u'' - f") ax.set_title(f"PINN residual at {N_FRESH} new random points") ax.legend(loc="upper right", fontsize=8) fig.savefig(out / "residual.png"); plt.close(fig) # ---------------------------------------------------------------------------- # Reports # ---------------------------------------------------------------------------- def fmt(v, p=2): if v is None or (isinstance(v, float) and math.isnan(v)): return "-" return f"{v:.{p}e}" def write_reports(out, conv, fem_n, fem_err, runs, notes, scale): keys = [("u", "deflection w"), ("slope", "slope"), ("M", "moment M"), ("V", "shear V")] L = [] L.append("# Simply supported beam: PINN vs FEM vs exact\n") L.append(f"Generated {time.strftime('%Y-%m-%d %H:%M')} by beam_compare.py. " f"Scaled units: u = w / w_max with w_max = 5qL^4/(384EI); scale factor {scale:g}.\n") L.append("Errors are measured against the exact solution on a uniform grid of " f"{N_GRID} points in [0, 1]. Relative L2 = ||a - exact|| / ||exact||; " "relative max = max|a - exact| / max|exact|.\n") if notes: L.append("## Checks and warnings\n") L += [f"- {n}" for n in notes] L.append("") if runs: L.append("## PINN runs\n") hdr = "| quantity | metric | " + " | ".join(r["label"] for r in runs) hdr += f" | FEM ({fem_n} el.) |" L.append(hdr) L.append("|" + "---|" * (3 + len(runs))) for k, name in keys: for m, mname in (("rel_l2", "rel. L2"), ("rel_max", "rel. max")): row = f"| {name} | {mname} | " + " | ".join( fmt(r["metrics"]["err"][k][m]) for r in runs) row += f" | {fmt(fem_err[k][m])} |" L.append(row) L.append("") L.append("| check | exact | " + " | ".join(r["label"] for r in runs) + " |") L.append("|" + "---|" * (2 + len(runs))) ex_vals = {"u(0.5)": 1.0, "slope(0)": scale / 24, "slope(1)": -scale / 24, "M(0.5)": scale / 8, "V(0)": scale / 2, "u(0)": 0.0, "u(1)": 0.0} for name, ev in ex_vals.items(): L.append(f"| {name} | {ev:.6g} | " + " | ".join( f"{r['metrics']['points'][name]:.6g}" for r in runs) + " |") L.append("| residual RMS, new points | 0 | " + " | ".join( fmt(r["metrics"]["residual_rms_fresh"]) for r in runs) + " |") L.append("| residual max, new points | 0 | " + " | ".join( fmt(r["metrics"]["residual_max_fresh"]) for r in runs) + " |") if any("results_csv_rel_l2" in r["metrics"] for r in runs): L.append("| rel. L2 at results.csv points | 0 | " + " | ".join( fmt(r["metrics"].get("results_csv_rel_l2")) for r in runs) + " |") L.append("") if len(runs) > 1: L.append("### Spread over runs (relative L2)\n") L.append("| quantity | min | median | max |") L.append("|---|---|---|---|") for k, name in keys: v = np.array([r["metrics"]["err"][k]["rel_l2"] for r in runs]) L.append(f"| {name} | {fmt(v.min())} | {fmt(float(np.median(v)))} | {fmt(v.max())} |") L.append("") L.append("Largest PINN error location (deflection): " + ", ".join( f"{r['label']} x = {r['metrics']['err']['u']['x_at_max']:.3f}" for r in runs) + "\n") L.append("## FEM mesh convergence (cubic Hermite elements, fourth-order equation)\n") L.append("| elements | max nodal error in u | rel. L2 u | order | rel. L2 slope | order | " "rel. L2 M | order | rel. L2 V | order |") L.append("|---|---|---|---|---|---|---|---|---|---|") for r in conv: L.append(f"| {r['elements']} | {fmt(r['nodal_max_err_u'])} | " + " | ".join(f"{fmt(r[f'rel_l2_{k}'])} | " + (f"{r.get(f'order_{k}', float('nan')):.2f}" if not math.isnan(r.get(f'order_{k}', float('nan'))) else "-") for k in ("u", "slope", "M", "V")) + " |") L.append("") L.append("Dimensional values: w = u * 5qL^4/(384EI); slope = (scaled slope) * 5qL^3/(384EI); " "M = (scaled M) * 5qL^2/384; V = (scaled V) * 5qL/384.\n") (out / "summary.md").write_text("\n".join(L)) payload = {"scale": scale, "grid_points": N_GRID, "fem_elements": fem_n, "fem_errors": fem_err, "fem_convergence": conv, "notes": notes, "runs": [{"label": r["label"], "checkpoint": r.get("checkpoint"), "metrics": r["metrics"]} for r in runs]} (out / "summary.json").write_text(json.dumps(payload, indent=2)) def write_csv(out, x, ex, fem_sol, runs): cols = {"x": x} for k in ("u", "slope", "M", "V"): cols[f"exact_{k}"] = ex[k] cols[f"fem_{k}"] = fem_sol[k] for r in runs: cols[f"{r['label'].replace(' ', '_')}_{k}"] = r["sol"][k] with open(out / "comparison.csv", "w", newline="") as f: w = csv.writer(f) w.writerow(cols.keys()) for row in zip(*cols.values()): w.writerow([f"{v:.10e}" for v in row]) # ---------------------------------------------------------------------------- # Main # ---------------------------------------------------------------------------- def run_metrics(x, sol, ex, fresh_x, fresh_res, scale): def at(k, xv): return float(np.interp(xv, x, sol[k])) return { "err": errors_vs_exact(x, sol, ex), "points": {"u(0.5)": at("u", 0.5), "slope(0)": at("slope", 0.0), "slope(1)": at("slope", 1.0), "M(0.5)": at("M", 0.5), "V(0)": at("V", 0.0), "u(0)": at("u", 0.0), "u(1)": at("u", 1.0)}, "residual_rms_fresh": float(np.sqrt(np.mean(fresh_res**2))), "residual_max_fresh": float(np.max(np.abs(fresh_res))), } def synthetic_runs(x, fresh_x, scale): """--selftest only: fake 'networks' = exact + a smooth wiggle. NOT a PINN.""" runs = [] for i, eps in enumerate((2e-4, 6e-4)): k = 2 * math.pi * (3 + i) def make(xx): ex = exact(xx, scale) g = eps * np.sin(math.pi * xx) * np.sin(k * xx) # derivatives of g, by finite differences on a fine grid then interpolated xf = np.linspace(0, 1, 40001); gf = eps * np.sin(math.pi * xf) * np.sin(k * xf) g1 = np.gradient(gf, xf); g2 = np.gradient(g1, xf); g3 = np.gradient(g2, xf) ip = lambda arr: np.interp(xx, xf, arr) return {"u": ex["u"] + g, "slope": ex["slope"] + ip(g1), "u2": ex["u2"] + ip(g2), "u3": ex["u3"] + ip(g3), "M": ex["M"] - ip(g2), "V": ex["V"] - ip(g3)} sol = make(x) fs = make(fresh_x) res = fs["u2"] - source(fresh_x, scale) runs.append({"label": f"synthetic {i + 1}", "sol": sol, "fresh_x": fresh_x, "fresh_res": res, "checkpoint": None}) return runs def main(): ap = argparse.ArgumentParser(description=__doc__.split("\n\n")[0]) ap.add_argument("--script", type=Path, help="the .py file generated by the Studio") ap.add_argument("--checkpoint", type=Path, nargs="*", default=[], help="one or more trained weight files (default: model.pt next to the script)") ap.add_argument("--results", type=Path, help="optional results.csv from the Studio script") ap.add_argument("--scale", type=float, default=SCALE, help="factor used in the Studio source (default 76.8 -> -38.4*x*(1-x))") ap.add_argument("--fem-elements", type=int, default=16, help="FEM mesh for the comparison") ap.add_argument("--out", type=Path, default=Path("beam_results")) ap.add_argument("--verify", action="store_true", help="FEM vs exact only, no PINN") ap.add_argument("--selftest", action="store_true", help="exercise the figures with synthetic data (not a PINN)") args = ap.parse_args() scale = args.scale out = args.out if not args.selftest else Path(str(args.out) + "_selftest") out.mkdir(parents=True, exist_ok=True) notes = [] x = np.linspace(0.0, 1.0, N_GRID) ex = exact(x, scale) rng = np.random.default_rng(FRESH_SEED) fresh_x = np.sort(rng.uniform(0.0, 1.0, N_FRESH)) # --- FEM --------------------------------------------------------------- conv = fem_convergence(scale) fem = BeamFEM(args.fem_elements) fem_sol = fem.evaluate(x, scale) fem_err = errors_vs_exact(x, fem_sol, ex) print(f"FEM, {args.fem_elements} elements: nodal deflection error " f"{fem.nodal_error(scale):.1e}; rel. L2 u {fem_err['u']['rel_l2']:.2e}, " f"M {fem_err['M']['rel_l2']:.2e}, V {fem_err['V']['rel_l2']:.2e}") # --- PINN -------------------------------------------------------------- runs = [] if args.selftest: runs = synthetic_runs(x, fresh_x, scale) notes.append("SELF-TEST: the 'runs' are synthetic curves (exact + a wiggle), not networks.") elif not args.verify: if args.script is None: ap.error("give --script (or use --verify for the FEM check only)") import torch script = args.script.resolve() cwd = os.getcwd() os.chdir(script.parent) # the script may look for files next to itself try: mod = load_studio_module(script) finally: os.chdir(cwd) net = find_net(mod) ckpts = args.checkpoint or [script.parent / "model.pt"] src_expected = source(fresh_x, scale) for i, ck in enumerate(ckpts): ck = ck if ck.is_absolute() else Path(cwd) / ck load_checkpoint(net, ck) sol = pinn_derivatives(mod, net, x) fs = pinn_derivatives(mod, net, fresh_x) res = fs["u2"] - src_expected label = f"run {i + 1}" # Did the Studio solve the equation we think it did? r_script = studio_residual(mod, net, fresh_x) if r_script is not None: gap = float(np.max(np.abs(r_script - res)) / np.max(np.abs(src_expected))) if gap > 1e-4: notes.append(f"{label}: the script's residual differs from u'' - " f"(-{scale / 2:g}*x*(1-x)) by {gap:.1e} (relative). Check the " "Source typed in the PDE block, and --scale.") elif i == 0: notes.append("Source check passed: the script's residual_pde1 equals " f"u'' + {scale / 2:g}*x*(1-x) (difference {gap:.1e}).") elif i == 0: notes.append("The script has no residual_pde1(); the source could not be cross-checked.") r = {"label": label, "checkpoint": str(ck), "sol": sol, "fresh_x": fresh_x, "fresh_res": res} r["metrics"] = run_metrics(x, sol, ex, fresh_x, res, scale) if args.results and i == len(ckpts) - 1: xs, us = read_results_csv(args.results) if len(xs): r["metrics"]["results_csv_rel_l2"] = rel_l2(us, exact(xs, scale)["u"]) r["metrics"]["results_csv_points"] = int(len(xs)) notes.append(f"results.csv ({len(xs)} points) is compared with {label}, " "the last checkpoint given: the script rewrites the table on every run.") runs.append(r) e = r["metrics"]["err"] print(f"{label} ({ck.name}): rel. L2 u {e['u']['rel_l2']:.2e}, slope " f"{e['slope']['rel_l2']:.2e}, M {e['M']['rel_l2']:.2e}, V {e['V']['rel_l2']:.2e}; " f"residual RMS {r['metrics']['residual_rms_fresh']:.2e}") if abs(r["metrics"]["points"]["u(0.5)"] + 1) < 0.2: notes.append(f"{label}: the deflection is upside down (u(0.5) close to -1). " "The source sign is wrong: it should be -38.4*x*(1-x).") for r in runs: if "metrics" not in r: r["metrics"] = run_metrics(x, r["sol"], ex, r["fresh_x"], r["fresh_res"], scale) # --- figures and reports ------------------------------------------------- fig_sketch(out) fig_convergence(out, conv) if runs: fig_deflection(out, x, ex, fem, args.fem_elements, runs) fig_derived(out, x, ex, fem_sol, args.fem_elements, runs) fig_derived_errors(out, x, ex, fem_sol, args.fem_elements, runs) fig_residual(out, runs) write_reports(out, conv, args.fem_elements, fem_err, runs, notes, scale) write_csv(out, x, ex, fem_sol, runs) print(f"Wrote {out}/ (summary.md, summary.json, comparison.csv and figures).") for n in notes: print("note:", n) if __name__ == "__main__": main()