A wave front meeting a harbor breakwater, diffracting through its gap into the sheltered water behind.
uv run python examples/cli.py run wave
The functions that pose and solve the problem. The figures are below the fold.
"""A wave front meeting a harbor breakwater, diffracting through its gap. `run` meshes the basin, sets a front travelling toward the wall, and steps the wave equation by Newmark, returning a `HarborStudy` of plain results. Nothing here draws: `figures.py` does that from the `HarborStudy`, and this file is what the gallery shows. """ from dataclasses import dataclass import numpy as np from fem.algebra.integrators import NewmarkMethod from fem.conditions import Conditions, Initial from fem.field import NodalField from fem.mesh.mesh import Mesh from fem.mesh.outline import Outline from fem.physics.equations import Wave from fem.post.solution import TransientSolution def harbor_outline(length: float = 6.0, width: float = 4.0, wall_x: float = 2.5, wall_thickness: float = 0.15, gap: float = 1.2) -> Outline: """A rectangular basin crossed by a breakwater with one gap, as a single loop. Open water lies left of the wall at `wall_x`, the sheltered harbor to its right. The two wall arms grow inward from the top and bottom edges, leaving `gap` open at mid-width. """ x0, x1 = wall_x, wall_x + wall_thickness y0, y1 = (width - gap) / 2, (width + gap) / 2 outline = np.array([ [0.0, 0.0], [x0, 0.0], [x0, y0], [x1, y0], [x1, 0.0], [length, 0.0], [length, width], [x1, width], [x1, y1], [x0, y1], [x0, width], [0.0, width], ]) return Outline.from_polygons([outline]) WALL_X, WALL_THICKNESS = 2.5, 0.15 @dataclass class HarborStudy: """Everything `run` computed, for the figures to read.""" mesh: Mesh u_initial: NodalField dudt_initial: NodalField solution: TransientSolution @property def u_values(self) -> np.ndarray: return self.solution.dofs @property def t_values(self) -> np.ndarray: return self.solution.t @property def harbor(self) -> np.ndarray: """Mask of the vertices on the sheltered side of the breakwater.""" return self.mesh.vertices[:, 0] > WALL_X + WALL_THICKNESS def run(c=1.0, front_x=1.0, front_width=0.25, dt=0.02, steps=400, min_angle=28, max_area=0.04, uniform_rounds=2) -> HarborStudy: """Mesh the basin, launch a front at the breakwater, and step it by Newmark.""" pslg = harbor_outline(wall_x=WALL_X, wall_thickness=WALL_THICKNESS) # Ruppert's meshes the outline coarsely; uniform red refinement then supplies the # resolution the front needs, keeping the angle bound at a fraction of the cost. mesh = pslg.mesh(min_angle=min_angle, max_area=max_area) for _ in range(uniform_rounds): mesh = mesh.refined() # A straight front on the open side, travelling toward the wall. Given d'Alembert's # pairing u = g(x - ct), du/dt = -c g'(x), so it moves one way instead of splitting. def profile(p): return np.exp(-((p[:, 0] - front_x) / front_width) ** 2) # No boundary conditions, so every edge is a wall: the natural du/dn = 0 reflects # a wave the same way up. bc = Conditions(Initial(profile, v0=lambda p: 2 * c * (p[:, 0] - front_x) / front_width**2 * profile(p))) wave = Wave(stiffness=c**2).problem(mesh, bc) solution = NewmarkMethod(dt=dt, steps=steps).solve(wave) return HarborStudy(mesh, wave.u0, wave.v0, solution)
"""The figures of the harbor breakwater demo, drawn from a `HarborStudy`.""" import numpy as np from demo_registry import Demo, DemoResult, Figure from demos.wave import physics from demos.wave.physics import HarborStudy, run from fem.plot.plotter import Plotter def _snapshot_steps(s: HarborStudy, n_shown) -> list[int]: """The steps the snapshot panels show, spread over the run once the front is under way.""" n = len(s.u_values) return [int(i) for i in np.linspace(n // 8, n - 1, n_shown)] def _colour_limits(s: HarborStudy, shown) -> tuple[float, float]: """One colour scale, set by the harbor side, so the diffracted wave reads even though it is far lower than the front that made it (which doubles again when it reflects off the far wall).""" span = float(max(abs(s.u_values[i][s.harbor]).max() for i in shown)) return (-span, span) def _animation_figure(s: HarborStudy, clim) -> Figure: animation = Plotter(1, 1, figsize=(7.4, 4.8)) animation.plot_animation(s.mesh, s.u_values, mode='colored', clim=clim, label='height', cmap='RdBu_r', titles=[f'Harbor breakwater t={t:.2f}' for t in s.t_values], idx=(0, 0)) return Figure(animation, 'Newmark time integration of the front.', 'animation') def _snapshots_figure(s: HarborStudy, shown, clim) -> Figure: snapshots = Plotter(2, 4, figsize=(18.0, 6.4), title='Diffraction through the gap') for panel, i in enumerate(shown): snapshots.plot(s.mesh, s.u_values[i], mode='colored', idx=divmod(panel, 4), title=f't={s.t_values[i]:.2f}', clim=clim, colorbar=panel == 7, cmap='RdBu_r', label='height') return Figure( snapshots, 'The front reaches the breakwater, reflects off the wall, and passes the ' 'gap, where it spreads into the harbor as a circular wave centred on the ' 'opening, lower than the front that made it. The later frames show that wave ' 'reflecting around the harbor while the front, reflected off the wall and ' 'then the far edge, comes back through the gap.', 'snapshots') def _setup_figure(s: HarborStudy) -> Figure: setup = Plotter(1, 3, figsize=(15.0, 3.8)) setup.plot(s.mesh, mode='mesh', idx=(0, 0), title='Basin and breakwater') setup.plot(s.mesh, s.u_initial, mode='colored', idx=(0, 1), label='height', title='Initial height u(x, 0)') setup.plot(s.mesh, s.dudt_initial, mode='colored', idx=(0, 2), label='velocity', title='Initial velocity, a front moving right') return Figure( setup, 'A basin with a breakwater across it, open on the left and sheltered on the ' 'right. The initial height and velocity together make a front travelling ' 'right; every edge is a wall, reflecting the wave the same way up.', 'conditions', setup=True) def demo(**kwargs) -> DemoResult: """A wave front meeting a harbor breakwater, diffracting through its gap into the sheltered water behind.""" s = run(**kwargs) shown = _snapshot_steps(s, 8) clim = _colour_limits(s, shown) return DemoResult([ _animation_figure(s, clim), _snapshots_figure(s, shown, clim), _setup_figure(s), ]) # Builds its own harbor basin, so it takes no domain. DEMO = Demo('wave', demo, section='Meshing & solving PDEs', show_source=physics, smoke_kwargs={'steps': 6, 'max_area': 0.5, 'uniform_rounds': 0})