full model: remove the rigid-body modes of the free plate after every batch

The von Karman strains are not invariant under a finite rotation and the model conserves neither the angular momentum about the in-plane axes nor, strictly, the in-plane momentum. After a hard strike the free plate slowly tilts and spins (0.2 rad after 2 s of a 2.6 J hit); the tilt then acts as a shell with that slope, pumps energy into the in-plane grid modes, and the sound swells again a second after the strike. The six mass-orthonormal rigid modes carry no vibration and no sound, so they are projected out of the state after every batch (remove_rigid=True, on by default).

Jeremy Magland <jmagland@flatironinstitute.org> committed 1 month ago commit d608f95fdf81 parent 20ab5aa Browse files

1 changed file+57−1

src/cymsim/full/simulation.pymodified+57−1View file

@@ -62,8 +62,10 @@ class Simulation:

6262 dtype: Literal["float32", "float64"] | torch.dtype = "float32",

6363 compile: bool | str = True,

6464 loss: tuple[float, float, float] | None = None,

65+ remove_rigid: bool = True,

6566 ):

6667 self.plate, self.strike, self.listener = plate, strike, listener

68+ self.remove_rigid = remove_rigid

6769 self.fs, self.fps, self.duration, self.movie = fs, fps, duration, movie

6870 self.device = torch.device(device) if device is not None else default_device()

6971 self.dtype = getattr(torch, dtype) if isinstance(dtype, str) else dtype

@@ -104,6 +106,12 @@ class Simulation:

104106 self.stick = {"x": 0.0, "v": -strike.velocity, "force": 0.0, "max_force": 0.0, "withdrawn": False}

105107 self.seconds = 0.0

106108

109+ # Mass-orthonormal rigid-body modes of the free plate (three transverse:

110+ # translation and two tilts; three in-plane: two translations and the

111+ # rotation about the normal), removed from the state after every batch.

112+ self._rigid_w, self._rigid_uv = self._rigid_bases()

113+ self._mass = torch.as_tensor(plate.node_mass, dtype=self.dtype, device=self.device)

114+

107115 self._step = step

108116 self._compile_pending = bool(compile)

109117 self._compile_mode = compile if isinstance(compile, str) else "default"

@@ -166,6 +174,8 @@ class Simulation:

166174 slots.append(slot)

167175 self._state, self._stick = state, stick

168176 self.total_steps += steps

177+ if self.remove_rigid:

178+ self._remove_rigid()

169179 if slots:

170180 self.pressure.append(self._audio(slots))

171181 self.samples += len(slots)

@@ -175,6 +185,48 @@ class Simulation:

175185 self.stick = {"x": float(st[0]), "v": float(st[1]), "force": float(st[2]), "max_force": float(st[3]), "withdrawn": bool(st[5] > 0)}

176186 self.seconds += time.perf_counter() - t0

177187

188+ def _rigid_bases(self):

189+ """Rigid-body modes orthonormalized under the mass inner product, as

190+ tensors ``(3, ny, nx)`` for w and ``(3, 2, ny, nx)`` for (u, v)."""

191+ g = self.plate.grid

192+ m = self.plate.node_mass

193+ cx, cy = self.plate.center

194+ X, Y, one = g.X - cx, g.Y - cy, np.ones(g.shape)

195+ z = np.zeros(g.shape)

196+

197+ def orthonormal(modes):

198+ out = []

199+ for v in modes:

200+ v = [np.where(g.mask, c, 0.0) for c in v]

201+ for b in out:

202+ c = sum((m * bi * vi).sum() for bi, vi in zip(b, v))

203+ v = [vi - c * bi for bi, vi in zip(b, v)]

204+ n = math.sqrt(sum((m * vi * vi).sum() for vi in v))

205+ out.append([vi / n for vi in v])

206+ return torch.as_tensor(np.asarray(out), dtype=self.dtype, device=self.device)

207+

208+ return orthonormal([[one], [X], [Y]])[:, 0], orthonormal([[one, z], [z, one], [-Y, X]])

209+

210+ def _remove_rigid(self) -> None:

211+ """Projects the rigid-body motion of the free plate out of the state, in place.

212+

213+ The von Karman strains are not invariant under a finite rotation, and

214+ the model conserves neither the angular momentum about the in-plane

215+ axes nor, strictly, the in-plane momentum, so a hard strike lets the

216+ free plate slowly tilt and spin. A tilt of a few degrees then acts as a

217+ shell with that slope and pumps energy into the in-plane grid modes.

218+ The rigid modes carry no vibration and, for a free plate, no sound, so

219+ they are removed after every batch."""

220+ u, v, w, wlo, du, dv, dw = self._state

221+ mass = self._mass

222+ for f in (w, dw):

223+ c = (mass * f * self._rigid_w).sum(dim=(1, 2))

224+ f.sub_((c[:, None, None] * self._rigid_w).sum(0))

225+ for a, b in ((u, v), (du, dv)):

226+ c = (mass * (a * self._rigid_uv[:, 0] + b * self._rigid_uv[:, 1])).sum(dim=(1, 2))

227+ a.sub_((c[:, None, None] * self._rigid_uv[:, 0]).sum(0))

228+ b.sub_((c[:, None, None] * self._rigid_uv[:, 1]).sum(0))

229+

178230 def _try_compile(self) -> None:

179231 """Wraps the step in ``torch.compile``, which fuses its some two hundred

180232 small tensor operations into a few kernels (about 6 times faster than

@@ -331,6 +383,7 @@ def simulate(

331383 dtype: Literal["float32", "float64"] | torch.dtype = "float32",

332384 compile: bool | str = True,

333385 progress: bool = False,

386+ remove_rigid: bool = True,

334387 ) -> Sound:

335388 """The sound of a strike by the full (geometrically nonlinear) model.

336389

@@ -343,7 +396,10 @@ def simulate(

343396 audio sample into a CUDA graph. ``fps`` is the movie frame rate and

344397 also the batch rate of the run. ``loss`` overrides the fitted

345398 ``(sigma0, sigma1, sigma_in)`` of the two-parameter loss operator.

399+ ``remove_rigid`` projects the rigid-body motion of the free plate out of

400+ the state after every batch (see ``Simulation._remove_rigid``); leave it

401+ on unless you study that motion itself.

346402 """

347- sim = Simulation(plate, strike, listener, duration=duration, fs=fs, fps=fps, movie=movie, device=device, dtype=dtype, compile=compile)

403+ sim = Simulation(plate, strike, listener, duration=duration, fs=fs, fps=fps, movie=movie, device=device, dtype=dtype, compile=compile, remove_rigid=remove_rigid)

348404 sim.run(progress=progress)

349405 return sim.sound()