Tangent runs log the energy parts of run and tangent per block (kinetic, bending, linear membrane; the run's true membrane energy), since the linear membrane energy is inflated at large deflection; tangent cache keys carry TANGENT_VERSION = 2 - proof-of-concept/cymsim

proof-of-concept / cymsim

Tangent runs log the energy parts of run and tangent per block (kinetic, bending, linear membrane; the run's true membrane energy), since the linear membrane energy is inflated at large deflection; tangent cache keys carry TANGENT_VERSION = 2

Jeremy Magland <jmagland@flatironinstitute.org> committed 2 days ago commit 240ff79d0d29 parent de5d1ce Browse files

5 changed files+53−13

src/cymsim/cache.pymodified+8−1View file

@@ -44,6 +44,10 @@ from .config import CONTINUATION_FIELDS, SimulationConfig

4444 log = logging.getLogger("cymsim")

4545

4646 MODEL_VERSION = 1

47+# Tangent runs (numerics.tangent) also carry TANGENT_VERSION in their key: bump it when what a tangent

48+# run records changes, so that older tangent entries are not taken for complete runs. 2: the energy

49+# parts of run and tangent are logged (simulate.NORM_COLUMNS), not only the energy norms.

50+TANGENT_VERSION = 2

4751 MAX_ENTRY_BYTES = 100 * 1024 * 1024 # the artifact-cache limit

4852 MAX_META_BYTES = 4096

4953 MAX_SEGMENTS = 20 # segments of history kept in an entry's metadata

@@ -72,7 +76,10 @@ def key_fields(cfg: SimulationConfig) -> dict:

7276 d = _non_default(cfg)

7377 for section, name in CONTINUATION_FIELDS:

7478 d.get(section, {}).pop(name, None)

75- return {k: v for k, v in d.items() if v != {}}

79+ out = {k: v for k, v in d.items() if v != {}}

80+ if cfg.numerics.tangent:

81+ out["tangent_version"] = TANGENT_VERSION

82+ return out

7683

7784

7885 def cache_key(cfg: SimulationConfig) -> str:

src/cymsim/cymsim_file.pymodified+7−4View file

@@ -23,7 +23,7 @@ import yaml

2323 from .config import SimulationConfig

2424 from .decimation import CUTOFF, DECIMATION_SPAN, KAISER_BETA

2525 from .mesh import Mesh

26-from .simulate import (NORM_COLUMNS, AccelerationRecord, DisplacementRecord, EnergyRecord, HitSummary,

26+from .simulate import (N_NORM_COLUMNS, NORM_COLUMNS, AccelerationRecord, DisplacementRecord, EnergyRecord, HitSummary,

2727 SimulationResult, SimulatorState, TangentRecord)

2828

2929 FORMAT_VERSION = 1

@@ -141,10 +141,13 @@ class CymsimWriter:

141141 if r.tangent.displacement is not None:

142142 g.create_dataset("displacement/w", shape=(0, N), maxshape=(None, N), dtype=sig, chunks=(64, N))

143143 g.create_dataset("displacement/time", shape=(0,), maxshape=(None,), dtype=np.float64, chunks=(4096,))

144- d = g.create_dataset("norms", shape=(0, 4), maxshape=(None, 4), dtype=np.float64, chunks=(1024, 4))

145- d.attrs["columns"] = NORM_COLUMNS

144+ ncol = r.tangent.norms.shape[1] if len(r.tangent.norms) else N_NORM_COLUMNS

145+ d = g.create_dataset("norms", shape=(0, ncol), maxshape=(None, ncol), dtype=np.float64, chunks=(1024, ncol))

146+ d.attrs["columns"] = " ".join(NORM_COLUMNS.split()[:ncol])

146147 d.attrs["description"] = ("after every block: time; kinetic_energy(dv) + linear_energy(dq) of the tangent; "

147- "kinetic_energy(v) and linear_energy(q) of the run (linear_energy: strain energy "

148+ "kinetic_energy(v) and linear_energy(q) of the run; then the tangent's kinetic, "

149+ "bending and linear membrane energy and the run's bending, linear membrane and true "

150+ "membrane energy (linear_energy: strain energy "

148151 "without the von Karman terms)")

149152

150153 g = f.create_group("energy")

src/cymsim/shell.pymodified+6−0View file

@@ -188,6 +188,12 @@ class ShellModel:

188188 q[self.tri_j], self.Gx_j, self.Gy_j, self.area_j, self.z0x_j, self.z0y_j,

189189 self.Cm_j, self.Ke_j, self.nl)

190190

191+ def linear_energy_parts(self, q):

192+ """(membrane, bending) strain energy of q without the von Karman terms; both quadratic forms."""

193+ Em, Eb = jax.vmap(element_energy_parts, in_axes=(0, 0, 0, 0, 0, 0, 0, 0, None))(

194+ q[self.tri_j], self.Gx_j, self.Gy_j, self.area_j, self.z0x_j, self.z0y_j, self.Cm_j, self.Ke_j, 0.0)

195+ return Em.sum(), Eb.sum()

196+

191197 def linear_energy(self, q):

192198 """Strain energy of q with the quadratic (von Karman) strain terms dropped: a quadratic

193199 form, so that sqrt(kinetic_energy(v) + linear_energy(q)) is a norm on the state."""

src/cymsim/simulate.pymodified+22−7View file

@@ -161,7 +161,15 @@ class SimulatorState:

161161 overhang_displacement_w: np.ndarray | None = None # (Td', N)

162162

163163

164-NORM_COLUMNS = "time tangent_norm2 kinetic linear_strain"

164+# After every block of a tangent run. The first four columns: the squared energy norm of the tangent,

165+# kinetic_energy(dv) + linear_energy(dq), and the two parts of the same norm of the run. Then the parts:

166+# of the tangent its kinetic, bending and linear membrane energy; of the run its bending, linear

167+# membrane and (nonlinear, true) membrane energy. At large deflection the linear membrane energy of

168+# the run, and presumably of the tangent, is far larger than the true one: the in-plane motion that

169+# cancels the quadratic stretch shows up as linear strain.

170+NORM_COLUMNS = ("time tangent_norm2 kinetic linear_strain "

171+ "tangent_kinetic tangent_bending tangent_membrane_linear bending membrane_linear membrane")

172+N_NORM_COLUMNS = len(NORM_COLUMNS.split())

165173

166174

167175 @dataclass

@@ -170,12 +178,13 @@ class TangentRecord:

170178 respect to delta when every prescribed peak force is multiplied by (1 + delta), the limit of

171179 (twin - run) / delta for twin runs. Signals as in SimulationResult; ``norms`` holds, after

172180 every block, the time, the squared energy norm of the tangent state, kinetic_energy(dv) +

173- linear_energy(dq), and the two parts of the same norm of the run's own state."""

181+ linear_energy(dq), the two parts of the same norm of the run's own state, and the energy parts of

182+ both (NORM_COLUMNS; files written before the parts were logged have the first four columns only)."""

174183 pressure_rayleigh: np.ndarray

175184 pressure_unbaffled: np.ndarray

176185 probe_acceleration: np.ndarray | None

177186 displacement: DisplacementRecord | None

178- norms: np.ndarray # (B, 4), columns NORM_COLUMNS

187+ norms: np.ndarray # (B, N_NORM_COLUMNS), columns NORM_COLUMNS

179188 state: SimulatorState | None = None

180189

181190

@@ -559,9 +568,13 @@ class Simulator:

559568

560569 @jax.jit

561570 def norms(st, dst):

562- """Squared energy norm of the tangent state, and the two parts of the run's own."""

563- return (sh.kinetic_energy(dst[1]) + sh.linear_energy(dst[0]),

564- sh.kinetic_energy(st[1]), sh.linear_energy(st[0]))

571+ """The NORM_COLUMNS after the time: energy norms and energy parts of tangent and run."""

572+ tk = sh.kinetic_energy(dst[1])

573+ tm, tb = sh.linear_energy_parts(dst[0])

574+ k = sh.kinetic_energy(st[1])

575+ ml, b = sh.linear_energy_parts(st[0])

576+ m, _ = sh.energy_parts_per_element(st[0])

577+ return (tk + tm + tb, k, ml + b, tk, tb, tm, b, ml, m.sum())

565578

566579 self._step, self._run_block, self._energies, self._clear = step, run_block, energies, clear

567580 self._tangent_block, self._norms = tangent_block, norms

@@ -617,6 +630,8 @@ class Simulator:

617630 raise ValueError("the earlier run stopped inside a block")

618631 if self.config.numerics.tangent and (prev.tangent is None or prev.tangent.state is None):

619632 raise ValueError("the earlier run has no tangent state")

633+ if self.config.numerics.tangent and prev.tangent.norms.shape[1] != N_NORM_COLUMNS:

634+ raise ValueError("the earlier tangent run logged fewer energy columns (written by an older cymsim)")

620635

621636 def _empty_outputs(self):

622637 N, dtype = self.mesh.n_nodes, self.shell.dtype

@@ -646,7 +661,7 @@ class Simulator:

646661 if tangent is not None:

647662 dstate, tstream, norms = tangent

648663 tf, tu, tp, td, tst = tstream.split(n_done, dstate)

649- trec = TangentRecord(tf, tu, tp, td, np.array(norms, np.float64).reshape(-1, 4), tst)

664+ trec = TangentRecord(tf, tu, tp, td, np.array(norms, np.float64).reshape(len(norms), N_NORM_COLUMNS), tst)

650665 from . import __version__

651666 return SimulationResult(

652667 self.config, self.fs, self.dt, self.substeps, self.mesh, self.shell.z0_node, self.shell.h_node,

tests/test_tangent.pymodified+10−1View file

@@ -84,8 +84,17 @@ def test_tangent_run_outputs_match_central_differences(coarse_config):

8484 assert np.linalg.norm(tan - fd) / np.linalg.norm(fd) < 1e-5

8585 assert np.array_equal(tg.displacement.time, res.displacement.time)

8686 norms = tg.norms

87- assert norms.shape[1] == 4 and len(norms) == res.n_steps_done // res.config.numerics.steps_per_block

87+ from cymsim.simulate import N_NORM_COLUMNS

88+ assert norms.shape[1] == N_NORM_COLUMNS and len(norms) == res.n_steps_done // res.config.numerics.steps_per_block

8889 assert np.all(norms[:, 1:] >= 0) and norms[-1, 1] > 0

90+ # the parts add up to the norms

91+ assert np.allclose(norms[:, 1], norms[:, 4] + norms[:, 5] + norms[:, 6], rtol=1e-12)

92+ assert np.allclose(norms[:, 3], norms[:, 7] + norms[:, 8], rtol=1e-12)

93+ # the run's kinetic, bending and true membrane energy are its energy log's at the reported times

94+ for e in res.energy_log:

95+ i = int(np.argmin(np.abs(norms[:, 0] - e.time)))

96+ if abs(norms[i, 0] - e.time) < 1e-12:

97+ assert np.allclose([norms[i, 2], norms[i, 7], norms[i, 9]], [e.kinetic, e.bending, e.membrane], rtol=1e-9)

8998

9099

91100 def test_tangent_run_continues_and_cuts(coarse_config, tmp_path):