1import type { MethodDefinition } from "./types";
2import {
3 CURVED_EDGE_EQUATIONS,
4 FREE_EDGE_EQUATIONS,
5 GOVERNING_EQUATIONS,
6 PINNED_CENTER_EQUATIONS,
7 RADIATION_EQUATIONS,
8 STRESS_EQUATIONS,
9} from "./common";
11export const FD_METHOD: MethodDefinition = {
12 id: "fd",
13 name: "Explicit finite differences",
14 shortName: "Explicit FD",
15 summary: "Explicit second-order time stepping of the 13-point biharmonic stencil.",
16 solverFunction: "plate_fd",
17 planFunction: "plate_fd_plan",
18 supports: {
19 shape: ["rectangle", "ellipse"],
20 bc: ["simply-supported", "clamped", "free", "pinned-center"],
21 excitation: ["impulse", "harmonic"],
22 audio: "clip",
23 },
24 readingOrder: [
25 "plate_fd.m",
26 "plate_fd_setup.m",
27 "plate_fd_loop.m",
28 "plate_fd_output.m",
29 "plate_fd_plan.m",
30 "plate_setup.m",
31 "plate_domain.m",
32 "plate_pad.m",
33 "plate_laplacian.m",
34 "plate_curvature.m",
35 "plate_matrices.m",
36 "plate_rigid_modes.m",
37 "plate_curvature_free.m",
38 "plate_curvature_domain.m",
39 "plate_force_time.m",
40 "plate_force_footprint.m",
41 "plate_stress.m",
42 "plate_listener.m",
43 "plate_spectrum.m",
44 "plate_result_schema.m",
45 ],
46 equations: [
47 ...GOVERNING_EQUATIONS,
48 {
49 tex: String.raw`(1+\sigma_0\Delta t)\,w^{n+1} = 2w^n - (1-\sigma_0\Delta t)\,w^{n-1} - \kappa^2\Delta t^2\, \mathsf L\mathsf L\, w^n + 2\sigma_1\Delta t\,\mathsf L\!\left(w^n - w^{n-1}\right) + \frac{\Delta t^2}{\mu} F^n s`,
50 caption:
51 "Explicit scheme, second order in space and time; L is the five-point Laplacian and LL the resulting 13-point biharmonic stencil.",
52 },
53 {
54 tex: String.raw`\mathsf L w_{ij} = \frac{w_{i,j+1} - 2w_{ij} + w_{i,j-1}}{\Delta x^2} + \frac{w_{i+1,j} - 2w_{ij} + w_{i-1,j}}{\Delta y^2}`,
55 caption:
56 "On a rectangle, supported and clamped edges enter through one ghost layer: odd reflection for simply supported (w = w_nn = 0), even reflection for clamped (w = w_n = 0). Free edges, and every condition on an ellipse, take the sparse operators below.",
57 },
58 {
59 tex: String.raw`\Delta t \le \frac{2}{\lambda_{\max}\left(\sigma_1 + \sqrt{\kappa^2 + \sigma_1^2}\right)}, \qquad \lambda_{\max} = \frac{4}{\Delta x^2} + \frac{4}{\Delta y^2}`,
60 caption:
61 "Stability bound from von Neumann analysis; for a square grid without sigma1 it reads dt <= dx^2 / (4 kappa). Because dt scales with dx^2, halving the spacing costs eight times more work. The sparse operators of the free plate and of the ellipse stay just under the same bound.",
62 },
63 ...FREE_EDGE_EQUATIONS,
64 ...PINNED_CENTER_EQUATIONS,
65 ...CURVED_EDGE_EQUATIONS,
66 ...STRESS_EQUATIONS,
67 ...RADIATION_EQUATIONS,
68 ],
69 validity: [
70 "Second-order truncation error in dx, dy and dt; the top modes on the grid are dispersive and their frequencies are too low.",
71 "The probe is snapped to the nearest grid node with a full stencil inside the plate.",
72 "Stress fields are evaluated only at the saved frames, so the envelope samples the true peak.",
73 "The time-step limit makes long simulated times expensive: cost grows as (1/dx)^4 for a fixed duration.",
74 "Audio is a short clip covering only the simulated window, box-averaged onto audio samples.",
75 "With the GPU option the time loop runs in single precision (a WebGPU port of plate_fd_loop.m); the setup, the stresses, the audio summation and the spectra stay in double precision in numbl. Against the double-precision loop the displacement differs by a relative L2 error of about 1e-6 over 50 ms on a 61 x 41 grid, the stress and the sound by about 1e-4.",
76 "Free edges: the moment and shear conditions are natural conditions of a discrete bending energy, met to first order at the boundary and second order inside; the three rigid-body motions are removed after every step, so the displacement is that of a plate on soft supports, not one drifting in space.",
77 "Pinned centre: the support is one grid node, so the stress right at it is resolution-dependent, and the two rigid rotations about it are removed rather than restrained by a modelled washer.",
78 "Ellipse: the edge is a staircase of grid nodes, so the geometry, the edge conditions and the stress along the edge are first order in the grid spacing. The stress on the outermost ring of nodes is copied from the nearest node with a full stencil, which puts the peak stress of a clamped edge about one spacing inside. For simply supported edges the curvature term that the staircase loses is added back (the plate paradox), so the disk's frequencies converge to the exact ones.",
79 ],
80 durationOf: (spec) => spec.fd.duration,
81 gridOf: (spec) => [spec.fd.nx, spec.fd.ny],
82 settingsKey: "fd",
83 usPerCellStep: 0.35,
84 gpu: {
85 loopFile: "plate_fd_loop.m",
86 setupFunction: "plate_fd_setup",
87 outputFunction: "plate_fd_output",
88 usPerCellStep: 0.02,
89 },
90};