2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 1function plan = plate_fd_plan(params)
2% PLATE_FD_PLAN Grid, time step and cost of the explicit finite-difference scheme.
3%
4% plan = plate_fd_plan(params) returns
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 5% plan.nx, plan.ny, plan.dx, plan.dy grid (nodes include the edges of the bounding box)
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 6% plan.dt_max largest stable time step [s]
7% plan.dt step actually used, params.fd.cfl * dt_max
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 8% plan.nsteps, plan.ncells, plan.cell_steps problem size (ncells counts every grid node)
9% plan.nunknowns nodes the scheme actually updates
11% Stability. Von Neumann analysis of the scheme in plate_fd for a Fourier
12% mode on which the discrete Laplacian has eigenvalue -lambda gives the
13% amplification polynomial
14% (1 + sigma0 dt) z^2 - (2 - dt^2 kappa^2 lambda^2 - 2 sigma1 dt lambda) z
15% + (1 - sigma0 dt - 2 sigma1 dt lambda) = 0.
16% Both roots stay inside the unit circle for every lambda up to the largest
17% eigenvalue of the five-point Laplacian, lambda_max = 4/dx^2 + 4/dy^2, iff
18% dt <= 2 / ( lambda_max (sigma1 + sqrt(kappa^2 + sigma1^2)) ),
19% which for a square grid without sigma1 is the familiar dt <= dx^2 / (4 kappa).
20% Because dt scales with dx^2, halving the grid spacing costs eight times
21% as many cell updates for the same simulated time.
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 23% The same bound serves the free and centre-pinned plates and the ellipse:
24% the Dirichlet Laplacian of any sub-domain has its spectrum inside that of
25% the full grid, the clamped correction of plate_matrices does not raise the
26% row sums above the interior stencil's, and the largest eigenvalue of the
27% energy-method operator M \ K lies just below lambda_max^2 (about 0.97 to
28% 0.99 of it on the grids used here), so the step is safe and slightly
29% conservative there.
31P = plate_setup(params);
32plan.nx = P.fd.nx;
33plan.ny = P.fd.ny;
34if plan.nx < 5 || plan.ny < 5
35 error('plate_fd_plan: the grid needs at least 5 nodes in each direction');
36end
37plan.dx = P.plate.Lx / (plan.nx - 1);
38plan.dy = P.plate.Ly / (plan.ny - 1);
39lambda_max = 4 / plan.dx^2 + 4 / plan.dy^2;
40s1 = P.damping.sigma1;
41plan.dt_max = 2 / (lambda_max * (s1 + sqrt(P.kappa^2 + s1^2)));
42plan.dt = P.fd.cfl * plan.dt_max;
43plan.nsteps = ceil(P.fd.duration / plan.dt);
44plan.ncells = plan.nx * plan.ny;
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 45dom = plate_domain(P.plate, plan.nx, plan.ny);
46if any(strcmp(P.bc, {'free', 'pinned-center'}))
47 plan.nunknowns = dom.nunknowns_free;
48else
49 plan.nunknowns = dom.nunknowns_dirichlet;
50end
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 51plan.cell_steps = plan.ncells * plan.nsteps;
52end