"""OFFLINE ONLY: generate Q1 FEM heat-decay data for the static website. Run: python scripts/generate_benchmark.py Dependencies are development-only; see requirements.txt. No Python is used by the website. This script assembles an actual FEM mass/stiffness system using scikit-fem, integrates it, and measures continuous error against an exact field. """ import json from pathlib import Path import numpy as np import scipy import skfem from scipy.sparse.linalg import splu from skfem import Basis, ElementQuad1, Functional, MeshQuad, asm from skfem.models import laplace, mass LX, LY = 1.0, 0.75 # metres ALPHA = 0.01 # square metres / second AMPLITUDE = 100.0 # K above a uniform reference temperature FINAL_TIME = 1.0 # seconds DT = 0.0025 # seconds; Crank-Nicolson LEVELS = (4, 8, 16, 32, 64) OUTPUT = Path(__file__).resolve().parents[1] / "data" / "heat-benchmark.json" def exact(x, t): decay = np.exp(-ALPHA * np.pi**2 * (LX**-2 + LY**-2) * t) return AMPLITUDE * np.sin(np.pi * x[0] / LX) * np.sin(np.pi * x[1] / LY) * decay def integrate(mass_matrix, stiffness, initial, interior, dt): """M du/dt + alpha K u = 0, with eliminated zero Dirichlet nodes.""" m = mass_matrix[interior][:, interior].tocsc() k = stiffness[interior][:, interior].tocsc() left = splu(m + (ALPHA * dt / 2) * k) right = m - (ALPHA * dt / 2) * k values = initial[interior].copy() steps = round(FINAL_TIME / dt) assert abs(steps * dt - FINAL_TIME) < 1e-14 for _ in range(steps): values = left.solve(right @ values) result = np.zeros_like(initial) result[interior] = values return result @Functional def squared_error(w): return (w.uh - exact(w.x, FINAL_TIME)) ** 2 def continuous_error(mesh, values, order=8): # Integrates the Q1 field, not just differences at mesh nodes. basis = Basis(mesh, ElementQuad1(), intorder=order) error_squared = asm(squared_error, basis, uh=basis.interpolate(values)) # Exact norm squared of the separable sine benchmark on the rectangle. norm_squared = (AMPLITUDE * np.exp( -ALPHA * np.pi**2 * (LX**-2 + LY**-2) * FINAL_TIME )) ** 2 * LX * LY / 4 return float(np.sqrt(error_squared / norm_squared)) def generate(): levels = [] for n in LEVELS: mesh = MeshQuad.init_tensor(np.linspace(0, LX, n + 1), np.linspace(0, LY, n + 1)) basis = Basis(mesh, ElementQuad1(), intorder=4) m, k = asm(mass, basis), asm(laplace, basis) interior = basis.complement_dofs(basis.get_dofs()) initial = exact(basis.doflocs, 0.0) values = integrate(m, k, initial, interior, DT) half_dt_values = integrate(m, k, initial, interior, DT / 2) # Independent uniform-grid check: this sampled sine is a generalized # eigenvector of the assembled Q1 matrices. Check against its closed # form Crank-Nicolson amplification (not used to generate the data). hx, hy = LX / n, LY / n eigenvalue = sum(6 * (1 - np.cos(np.pi / n)) / (h**2 * (2 + np.cos(np.pi / n))) for h in (hx, hy)) q = ALPHA * eigenvalue * DT / 2 expected = initial * ((1 - q) / (1 + q)) ** round(FINAL_TIME / DT) assert np.max(np.abs(values - expected)) < 1e-9, "Independent eigenmode check failed" error = continuous_error(mesh, values) error_order10 = continuous_error(mesh, values, order=10) difference = half_dt_values - values exact_norm = (AMPLITUDE * np.exp( -ALPHA * np.pi**2 * (LX**-2 + LY**-2) * FINAL_TIME )) * np.sqrt(LX * LY / 4) time_sensitivity = float(np.sqrt(difference @ m @ difference) / exact_norm) # Put coefficients into rows of increasing y; x increases within a row. grid = np.zeros((n + 1, n + 1)) ix = np.rint(basis.doflocs[0] * n / LX).astype(int) iy = np.rint(basis.doflocs[1] * n / LY).astype(int) grid[iy, ix] = values rounded = np.round(grid, 6) round_difference = rounded[iy, ix] - values rounding_error = float(np.sqrt(round_difference @ m @ round_difference) / exact_norm) assert mesh.nelements == n * n assert np.max(np.abs(values[basis.get_dofs().all()])) == 0 assert np.all(np.isfinite(values)) and np.min(values) >= -1e-10 assert time_sensitivity < error * 0.001, "Temporal error is not negligible" assert abs(error - error_order10) < 1e-10, "Error quadrature has not converged" assert rounding_error < error * 0.001, "Serialization is too coarse" levels.append({ "n": n, "elements": n * n, "nodes": int(basis.N), "hx": LX / n, "hy": LY / n, "relativeL2Error": error, "dtHalvingDifference": time_sensitivity, "quadratureErrorDifference": abs(error - error_order10), "serializationRelativeL2Difference": rounding_error, "temperature": rounded.ravel().tolist(), }) print(f"{n:2} x {n:2}: relative L2 error = {error:.9e}, dt sensitivity = {time_sensitivity:.3e}") errors = [level["relativeL2Error"] for level in levels] assert all(a > b > 0 for a, b in zip(errors, errors[1:])) rates = [float(np.log2(a / b)) for a, b in zip(errors, errors[1:])] assert all(1.8 < rate < 2.2 for rate in rates), rates result = { "schemaVersion": 1, "benchmark": "Homogeneous 2D heat equation, rectangular sine-mode decay", "domain": {"lx": LX, "ly": LY, "units": "m"}, "alpha": ALPHA, "amplitude": AMPLITUDE, "time": FINAL_TIME, "temperatureUnits": "K above a uniform reference temperature", "method": "Continuous bilinear Q1 quadrilateral FEM; consistent mass matrix; nodal initial interpolation; homogeneous Dirichlet boundaries; Crank-Nicolson time integration", "dt": DT, "steps": round(FINAL_TIME / DT), "quadratureOrder": 8, "errorDefinition": "sqrt(integral_domain (T_h - T_exact)^2 dA / integral_domain T_exact^2 dA)", "normalization": "dimensionless", "layout": "Nodal values, row-major, x increasing fastest; y increases from 0 to Ly", "generator": "scripts/generate_benchmark.py", "versions": {"scikit-fem": skfem.__version__, "numpy": np.__version__, "scipy": scipy.__version__}, "observedOrders": rates, "levels": levels, } OUTPUT.parent.mkdir(exist_ok=True) OUTPUT.write_text(json.dumps(result, separators=(",", ":")) + "\n", encoding="utf-8") print(f"Saved {OUTPUT.name} ({OUTPUT.stat().st_size:,} bytes). Observed orders: {rates}") if __name__ == "__main__": generate()