diff --git a/CHANGELOG.md b/CHANGELOG.md
index 46f345d..89f601e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,4 +1,29 @@
-# Unreleased
+# v0.8.0
+- Add the `GpuBoundaryCondition::non_reflecting` (absorbing) boundary condition, based on
+  Lysmer-Kuhlemeyer viscous dashpots. It lets outgoing elastic waves leave the domain instead of
+  being reflected back into it, emulating an unbounded medium. See the new `non_reflecting2` 2D
+  demo for a side-by-side comparison with a reflecting boundary.
+- Add `ParticleModel::absorbing_pml`, a perfectly-matched-layer absorbing material after
+  Kurima, Chandra & Soga (arXiv:2407.02790). Particles carrying it form a layer around the region
+  of interest whose coordinates are stretched, so outgoing waves slow and spread instead of
+  returning; pair it with `ParticleDynamics::damping` over the same particles to dissipate them.
+  `models::pml_stretch` computes the per-particle stretch from the layer geometry. It absorbs
+  better than the dashpot boundary above (~0.2% vs ~1% of a reflecting wall's residual motion on
+  a 2D impulse test) at the cost of the extra particles the layer needs.
+- The PML and the two items below sit behind the new `pml` cargo feature, off by default, mirrored
+  in the shaders as `SLOSH_PML`. It shifts `DefaultParticleModelType::LEN` (5 on, 4 off).
+- With `pml` on, grid nodes carry a per-direction mass (`Node.directional_mass`), so a material can
+  rescale its own inertia per axis via `ModelUpdateResult::mass_scale`. The momentum update divides
+  by it while gravity keeps acting on the real mass. Ordinary materials report a scale of one. A
+  P2G hook replacing the built-in transfer must write this field too.
+- Add `ParticleDynamics::stiffness_damping`, the stiffness-proportional half of Rayleigh damping.
+  It adds a viscous stress `a_K * C : sym(grad v)` from each material's own elastic tensor, so
+  unlike the mass-proportional `damping` it is blind to rigid-body motion. Supported by every
+  built-in model. It tightens the explicit stability bound, which `WgTimestepBounds` accounts for.
+- Grid nodes now get a collision reported up to `COLLISION_REPORT_CELLS` (6.5) cells from a
+  collider instead of 1.5, so the absorbing boundary above can grade its damping over a band
+  several cells deep. The contact boundary conditions gate themselves on the distance to the
+  surface and are unaffected.
 - Update to Rapier 0.32. This migrates most public APIs and internals to use `glam` instead of `nalgebra`.
 - Fix a GPU validation error / panic on simulations with more than ~4.19M particles, caused by
   compute kernels dispatching more than 65535 workgroups along a single dimension. The affected 
diff --git a/Cargo.toml b/Cargo.toml
index 0946f4d..c8a22ad 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -26,7 +26,7 @@ stensor = "0.4.2"
 
 [workspace.lints]
 rust.unexpected_cfgs = { level = "warn", check-cfg = [
-    'cfg(feature, values("dim2", "dim3", "cpic", "node_particle_lists"))',
+    'cfg(feature, values("dim2", "dim3", "cpic", "node_particle_lists", "pml"))',
 ] }
 
 [patch.crates-io]
diff --git a/crates/slosh2d/Cargo.toml b/crates/slosh2d/Cargo.toml
index f0d8aa3..9047457 100644
--- a/crates/slosh2d/Cargo.toml
+++ b/crates/slosh2d/Cargo.toml
@@ -3,7 +3,7 @@ name = "slosh2d"
 authors = ["Sébastien Crozet <sebcrozet@dimforge.com>"]
 description = "Cross-platform GPU 2D Material Point Method implementation."
 repository = "https://github.com/dimforge/slosh"
-version = "0.7.0"
+version = "0.8.0"
 edition = "2024"
 license = "Apache-2.0"
 
@@ -15,11 +15,18 @@ name = "slosh2d"
 path = "../../src/lib.rs"
 required-features = ["dim2"]
 
+# The other examples are auto-discovered; this one is declared only to carry `required-features`.
+[[example]]
+name = "non_reflecting2"
+path = "examples/non_reflecting2.rs"
+required-features = ["pml"]
+
 [features]
 default = ["dim2", "cpic", "node_particle_lists"]
 dim2 = []
 cpic = []
 node_particle_lists = []
+pml = []
 
 comptime = ["slosh_testbed2d/comptime", "stensor/comptime"]
 runtime = ["slosh_testbed2d/runtime", "stensor/runtime"]
diff --git a/crates/slosh2d/examples/non_reflecting2.rs b/crates/slosh2d/examples/non_reflecting2.rs
new file mode 100644
index 0000000..347f8a7
--- /dev/null
+++ b/crates/slosh2d/examples/non_reflecting2.rs
@@ -0,0 +1,252 @@
+//! Three ways of terminating a domain, side by side. Three identical elastic blocks get the same
+//! radial pulse at their centre; what differs is what happens when the wave reaches the edge:
+//!
+//! - **left**: [`GpuBoundaryCondition::slip`] walls, a perfect reflector, the reference.
+//! - **middle**: [`GpuBoundaryCondition::non_reflecting_for_material`] walls, Lysmer-Kuhlemeyer
+//!   viscous dashpots graded over a band a few cells deep. Costs nothing but the boundary itself.
+//! - **right**: a perfectly-matched-layer skirt of [`ParticleModel::absorbing_pml`] particles
+//!   enclosed by fixed walls, after Kurima, Chandra & Soga (arXiv:2407.02790). Pays for it in
+//!   particles: the skirt more than doubles the block's particle count.
+//!
+//! Residual motion in the interior settles at roughly 1% of the reflecting block for the dashpot
+//! and 0.3% for the PML.
+//!
+//! Switch the testbed's render mode to **volume** to see the waves; they are well under 1% strain
+//! and invisible in the default mode. Compression reads orange, dilation cyan, shear white.
+
+use slosh_testbed2d::{RapierData, slosh};
+
+use glam::{Vec2, Vec4, vec2};
+use rapier2d::prelude::{
+    ColliderBuilder, ColliderHandle, ColliderSet, RigidBodyBuilder, RigidBodySet,
+};
+use slang_hal::backend::WebGpu;
+use slosh::{
+    pipeline::MpmData,
+    solver::{Particle, SimulationParams},
+};
+use slosh_testbed2d::{AppState, PhysicsContext};
+use slosh2d::models::{DEFAULT_PML_MAX_STRETCH, pml_stretch};
+use slosh2d::solver::{GpuBoundaryCondition, ParticleModel};
+
+#[allow(dead_code)]
+fn main() {
+    panic!("Run the `testbed2` example instead.");
+}
+
+/// Size of each block's region of interest.
+const BLOCK_WIDTH: f32 = 14.0;
+const BLOCK_HEIGHT: f32 = 10.0;
+/// Horizontal gap between blocks.
+const BLOCK_GAP: f32 = 5.0;
+
+/// Depth of the absorbing skirt around the third block: 15 cells, enough for this small domain
+/// (the paper uses layers tens of cells deep).
+const PML_THICKNESS: f32 = 3.0;
+
+/// Both halves of the Rayleigh damping applied over the skirt, ramped with the stretch so they
+/// stay continuous at the interface. The two are complementary in frequency and beat either alone
+/// here; raising either from these values makes absorption worse, and 0.010 stiffness damping
+/// would already drag the timestep bound below this demo's fixed timestep.
+const PML_DAMPING: f32 = 40.0;
+const PML_STIFFNESS_DAMPING: f32 = 0.008;
+
+const CELL_WIDTH: f32 = 0.2;
+const PARTICLES_PER_CELL_DIM: usize = 2;
+
+const YOUNG_MODULUS: f32 = 8.0e5;
+const POISSON_RATIO: f32 = 0.2;
+const DENSITY: f32 = 1000.0;
+
+/// Radius of the initial velocity pulse, and its peak velocity. The induced strain is roughly
+/// `PULSE_VELOCITY / c_p`, well within the linear elastic regime.
+const PULSE_RADIUS: f32 = 1.5;
+const PULSE_VELOCITY: f32 = 0.25;
+
+#[derive(Copy, Clone, PartialEq)]
+enum Termination {
+    Reflecting,
+    Dashpot,
+    Pml,
+}
+
+impl Termination {
+    /// Depth of absorbing material outside the region of interest.
+    fn skirt(self) -> f32 {
+        match self {
+            Termination::Pml => PML_THICKNESS,
+            _ => 0.0,
+        }
+    }
+
+    fn wall_condition(self) -> GpuBoundaryCondition {
+        match self {
+            Termination::Reflecting => GpuBoundaryCondition::slip(),
+            Termination::Dashpot => GpuBoundaryCondition::non_reflecting_for_material(
+                YOUNG_MODULUS,
+                POISSON_RATIO,
+                DENSITY,
+                0.0,
+            ),
+            // The paper encloses the absorbing layer in fixed displacement boundaries.
+            Termination::Pml => GpuBoundaryCondition::stick(),
+        }
+    }
+}
+
+pub fn non_reflecting_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
+    let mut rapier_data = RapierData::default();
+
+    let diameter = CELL_WIDTH / PARTICLES_PER_CELL_DIM as f32;
+    let mut particles = vec![];
+    let mut materials = vec![];
+
+    let mut x = 0.0;
+    for termination in [
+        Termination::Reflecting,
+        Termination::Dashpot,
+        Termination::Pml,
+    ] {
+        let skirt = termination.skirt();
+        x += skirt;
+
+        let mins = vec2(x, skirt);
+        let maxs = mins + vec2(BLOCK_WIDTH, BLOCK_HEIGHT);
+        let center = (mins + maxs) / 2.0;
+
+        let ni = ((BLOCK_WIDTH + 2.0 * skirt) / diameter).ceil() as usize;
+        let nj = ((BLOCK_HEIGHT + 2.0 * skirt) / diameter).ceil() as usize;
+
+        for i in 0..ni {
+            for j in 0..nj {
+                let position =
+                    mins - Vec2::splat(skirt) + vec2(i as f32 + 0.5, j as f32 + 0.5) * diameter;
+
+                // Outside the region of interest, particles become absorbing: stretched more the
+                // deeper into the skirt, and carrying the damping that dissipates the wave.
+                let stretch =
+                    pml_stretch(position, mins, maxs, PML_THICKNESS, DEFAULT_PML_MAX_STRETCH);
+                let absorbing = termination == Termination::Pml && stretch.length_squared() > 0.0;
+
+                let model = if absorbing {
+                    ParticleModel::absorbing_pml(YOUNG_MODULUS, POISSON_RATIO, stretch)
+                } else {
+                    ParticleModel::elastic(YOUNG_MODULUS, POISSON_RATIO)
+                };
+
+                let mut particle = Particle::new(position, diameter / 2.0, DENSITY, model);
+                if absorbing {
+                    let ramp = stretch.max_element() / DEFAULT_PML_MAX_STRETCH;
+                    particle.dynamics.damping = PML_DAMPING * ramp;
+                    particle.dynamics.stiffness_damping = PML_STIFFNESS_DAMPING * ramp;
+                }
+
+                // Radial pulse with a smooth (Hann) taper so it doesn't inject high frequencies
+                // the grid can't resolve.
+                let dpos = position - center;
+                let dist = dpos.length();
+
+                if dist < PULSE_RADIUS && dist > 1.0e-6 {
+                    let taper = 0.5 * (1.0 + (std::f32::consts::PI * dist / PULSE_RADIUS).cos());
+                    particle.dynamics.velocity = dpos / dist * PULSE_VELOCITY * taper;
+                }
+
+                particles.push(particle);
+            }
+        }
+
+        let walls = build_walls(
+            &mut rapier_data.bodies,
+            &mut rapier_data.colliders,
+            mins - Vec2::splat(skirt),
+            maxs + Vec2::splat(skirt),
+        );
+        let condition = termination.wall_condition();
+        materials.extend(walls.into_iter().map(|handle| (handle, condition)));
+
+        x += BLOCK_WIDTH + skirt + BLOCK_GAP;
+    }
+
+    let total_width = x - BLOCK_GAP;
+
+    if !app_state.restarting {
+        // Fixed substepping: the CFL bound for this material is `0.5 * h / c_p` ≈ 3.3ms, so
+        // 8 substeps per 1/60s frame leaves a comfortable margin.
+        app_state.min_num_substeps = 8;
+        app_state.max_num_substeps = 8;
+        // A pure wave propagation test: neither absorbing termination has the static stiffness to
+        // hold the blocks up against gravity.
+        app_state.gravity_factor = 0.0;
+        app_state.initial_camera2d_at = Some([total_width / 2.0, BLOCK_HEIGHT / 2.0]);
+        app_state.initial_camera2d_zoom = Some(950.0 / total_width);
+    }
+
+    // Uniform base color so the blocks look alike in the default render mode (the `volume` mode
+    // derives its colors from the strain and ignores this).
+    app_state.particle_colors = Some(vec![Vec4::new(0.2, 0.25, 0.3, 1.0); particles.len()]);
+
+    let params = SimulationParams {
+        gravity: vec2(0.0, -9.81) * app_state.gravity_factor,
+        dt: 1.0 / 60.0,
+        padding: 0.0,
+    };
+
+    let data = MpmData::new(
+        backend,
+        params,
+        &particles,
+        &rapier_data.bodies,
+        &rapier_data.colliders,
+        &materials,
+        CELL_WIDTH,
+        30_000,
+    )
+    .unwrap();
+
+    PhysicsContext {
+        data,
+        rapier_data,
+        callbacks: vec![],
+        hooks_state: None,
+    }
+}
+
+/// Encloses the box `[mins, maxs]` in four fixed walls whose inner faces are flush with it.
+fn build_walls(
+    bodies: &mut RigidBodySet,
+    colliders: &mut ColliderSet,
+    mins: Vec2,
+    maxs: Vec2,
+) -> [ColliderHandle; 4] {
+    const THICKNESS: f32 = 1.0;
+
+    let center = (mins + maxs) / 2.0;
+    let half = (maxs - mins) / 2.0;
+
+    let walls = [
+        // Bottom, top, left, right.
+        (
+            vec2(center.x, mins.y - THICKNESS),
+            vec2(half.x + THICKNESS, THICKNESS),
+        ),
+        (
+            vec2(center.x, maxs.y + THICKNESS),
+            vec2(half.x + THICKNESS, THICKNESS),
+        ),
+        (
+            vec2(mins.x - THICKNESS, center.y),
+            vec2(THICKNESS, half.y + THICKNESS),
+        ),
+        (
+            vec2(maxs.x + THICKNESS, center.y),
+            vec2(THICKNESS, half.y + THICKNESS),
+        ),
+    ];
+
+    walls.map(|(translation, half_extents)| {
+        let rb = RigidBodyBuilder::fixed().translation(translation);
+        let rb_handle = bodies.insert(rb);
+        let co = ColliderBuilder::cuboid(half_extents.x, half_extents.y);
+        colliders.insert_with_parent(co, rb_handle, bodies)
+    })
+}
diff --git a/crates/slosh2d/examples/testbed2.rs b/crates/slosh2d/examples/testbed2.rs
index d218299..2c81edb 100644
--- a/crates/slosh2d/examples/testbed2.rs
+++ b/crates/slosh2d/examples/testbed2.rs
@@ -1,15 +1,26 @@
 mod centilever_beam2;
 mod elastic_cut2;
 mod elasticity2;
+#[cfg(feature = "pml")]
+mod non_reflecting2;
 mod sand2;
 
 #[kiss3d::main]
 pub async fn main() {
-    slosh_testbed2d::run(vec![
+    #[allow(unused_mut)]
+    let mut scenes: slosh_testbed2d::SceneBuilders<_> = vec![
         ("centilever beam".to_string(), centilever_beam2::beam_demo),
         ("sand".to_string(), sand2::sand_demo),
         ("elasticity".to_string(), elasticity2::elasticity_demo),
         ("elastic_cut".to_string(), elastic_cut2::elastic_cut_demo),
-    ])
-    .await;
+    ];
+
+    // The absorbing skirt needs `ParticleModel::absorbing_pml`.
+    #[cfg(feature = "pml")]
+    scenes.push((
+        "non-reflecting boundary".to_string(),
+        non_reflecting2::non_reflecting_demo,
+    ));
+
+    slosh_testbed2d::run(scenes).await;
 }
diff --git a/crates/slosh3d/Cargo.toml b/crates/slosh3d/Cargo.toml
index 1cb8ad4..f6af925 100644
--- a/crates/slosh3d/Cargo.toml
+++ b/crates/slosh3d/Cargo.toml
@@ -3,7 +3,7 @@ name = "slosh3d"
 authors = ["Sébastien Crozet <sebcrozet@dimforge.com>"]
 description = "Cross-platform GPU 3D Material Point Method implementation."
 repository = "https://github.com/dimforge/slosh"
-version = "0.7.0"
+version = "0.8.0"
 edition = "2024"
 license = "Apache-2.0"
 
@@ -20,6 +20,7 @@ default = ["dim3", "cpic", "node_particle_lists"]
 dim3 = []
 cpic = []
 node_particle_lists = []
+pml = []
 
 comptime = ["stensor/comptime", "slosh_testbed3d/comptime"]
 runtime = ["stensor/runtime", "slosh_testbed3d/runtime"]
diff --git a/crates/slosh3d/examples/heightfield3.rs b/crates/slosh3d/examples/heightfield3.rs
index 5ad0346..8976bc1 100644
--- a/crates/slosh3d/examples/heightfield3.rs
+++ b/crates/slosh3d/examples/heightfield3.rs
@@ -40,8 +40,8 @@ pub fn heightfield_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCo
     }
 
     if !app_state.restarting {
-        app_state.min_num_substeps = 10;
-        app_state.max_num_substeps = 10;
+        app_state.min_num_substeps = 20;
+        app_state.max_num_substeps = 20;
         app_state.gravity_factor = 1.0;
     };
 
diff --git a/crates/slosh_testbed2d/Cargo.toml b/crates/slosh_testbed2d/Cargo.toml
index dcd0a87..1ac9c99 100644
--- a/crates/slosh_testbed2d/Cargo.toml
+++ b/crates/slosh_testbed2d/Cargo.toml
@@ -2,7 +2,7 @@
 name = "slosh_testbed2d"
 authors = ["Sébastien Crozet <sebcrozet@dimforge.com>"]
 description = "Testbed for the slosh GPU-based MPM physics simulation library (2D version)."
-version = "0.7.0"
+version = "0.8.0"
 license = "Apache-2.0"
 edition = "2021"
 
@@ -38,7 +38,7 @@ futures-test = "0.3"
 serial_test = "3"
 approx = "0.5"
 async-std = { version = "1", features = ["attributes"] }
-slosh2d = { version = "0.7", path = "../slosh2d" }
+slosh2d = { version = "0.8", path = "../slosh2d" }
 regex = "1"
 web-time = "1"
 
diff --git a/crates/slosh_testbed3d/Cargo.toml b/crates/slosh_testbed3d/Cargo.toml
index 2d7f963..b166def 100644
--- a/crates/slosh_testbed3d/Cargo.toml
+++ b/crates/slosh_testbed3d/Cargo.toml
@@ -2,7 +2,7 @@
 name = "slosh_testbed3d"
 authors = ["Sébastien Crozet <sebcrozet@dimforge.com>"]
 description = "Testbed for the slosh GPU-based MPM physics simulation library (3D version)."
-version = "0.7.0"
+version = "0.8.0"
 license = "Apache-2.0"
 edition = "2021"
 
@@ -38,7 +38,7 @@ futures-test = "0.3"
 serial_test = "3"
 approx = "0.5"
 async-std = { version = "1", features = ["attributes"] }
-slosh3d = { version = "0.7", path = "../slosh3d" }
+slosh3d = { version = "0.8", path = "../slosh3d" }
 regex = "1"
 web-time = "1"
 
diff --git a/shaders/slosh/collision/collide.slang b/shaders/slosh/collision/collide.slang
index 6d6d32e..b2d0df1 100644
--- a/shaders/slosh/collision/collide.slang
+++ b/shaders/slosh/collision/collide.slang
@@ -84,6 +84,11 @@ public struct GridCollision {
 }
 
 
+// How far (in cells) from a collider’s surface a grid node still gets a collision reported.
+// Must be at least 1.5. Larger values (e.g. 6.5) increase the effectiveness of the Dashpots
+// non-reflective boundaries but are slightly more expensive and have more effect on movent.
+public static const float COLLISION_REPORT_CELLS = 1.5;
+
 public func collide(
     StructuredBuffer<Shape> collision_shapes,
     StructuredBuffer<Pose> collision_shape_poses,
@@ -97,7 +102,7 @@ public func collide(
     static const float MAX_FLT = 1.0e10; // Is the f32::MAX constant defined somewhere in WGSL?
     var collision = GridCollision(Vect<float>(0.0), Vect<float>(0.0), 0, false);
 
-    let dist_cap = Vect<float>(cell_width * 1.5);
+    let dist_cap = Vect<float>(cell_width * COLLISION_REPORT_CELLS);
     var best_dist = MAX_FLT;
 
     // TODO: don’t  rely on the array length, e.g., if the user wants to
diff --git a/shaders/slosh/features.h b/shaders/slosh/features.h
new file mode 100644
index 0000000..0549ae5
--- /dev/null
+++ b/shaders/slosh/features.h
@@ -0,0 +1,37 @@
+// Defaults for the feature macros that gate parts of the GPU data layout.
+//
+// `register_shaders` sets all of these from the cargo features of the same (lowercase) name, so
+// the shader structs match their Rust counterparts. The defaults below are what a standalone
+// `slangc` invocation, or a build that never calls `register_shaders`, gets; they mirror slosh's
+// own defaults. Slang warns on an undefined identifier in `#if`, so every macro needs one.
+//
+// Include this from any module that reads one of them:
+//
+//     #include "slosh/features.h"
+//
+// This file holds preprocessor directives only and carries no `module` declaration, so it is a
+// plain textual include rather than a slang module (the build helpers only compile `*.slang`).
+
+#ifndef SLOSH_FEATURES_H
+#define SLOSH_FEATURES_H
+
+// Gates the CPIC parts of Node: the incompatible-momentum lane and the cdf collision field.
+// Keep in lockstep with `GpuGridNode` on the Rust side.
+#ifndef SLOSH_CPIC
+#define SLOSH_CPIC 1
+#endif
+
+// Gates the per-node particle linked lists (built in finalize_particles_sort, cleared in reset).
+// The bound buffers and the Rust `GridArgs` fields share it so everything stays in sync.
+#ifndef SLOSH_NODE_PARTICLE_LISTS
+#define SLOSH_NODE_PARTICLE_LISTS 1
+#endif
+
+// Gates everything the absorbing perfectly-matched layer adds: the per-axis mass on Node, the
+// per-particle mass scale and stiffness damping, and the model itself. Off by default, like the
+// cargo feature.
+#ifndef SLOSH_PML
+#define SLOSH_PML 0
+#endif
+
+#endif // SLOSH_FEATURES_H
diff --git a/shaders/slosh/grid/grid.slang b/shaders/slosh/grid/grid.slang
index 6b771a4..12e0069 100644
--- a/shaders/slosh/grid/grid.slang
+++ b/shaders/slosh/grid/grid.slang
@@ -2,19 +2,7 @@ module grid;
 
 import slosh.aliases;
 
-// SLOSH_CPIC gates the CPIC parts of Node: the incompatible-momentum lane and the cdf collision
-// field. register_shaders sets it from the cpic cargo feature so this layout matches the Rust
-// GpuGridNode. Default on so a compile that never sets the macro still builds the full layout.
-#ifndef SLOSH_CPIC
-#define SLOSH_CPIC 1
-#endif
-
-// SLOSH_NODE_PARTICLE_LISTS gates the per-node particle linked lists (built in
-// finalize_particles_sort, cleared in reset). Set from the node_particle_lists cargo feature; the
-// bound buffers and the Rust GridArgs fields share it so everything stays in sync.
-#ifndef SLOSH_NODE_PARTICLE_LISTS
-#define SLOSH_NODE_PARTICLE_LISTS 1
-#endif
+#include "slosh/features.h"
 
 // TODO: a lot if what exposed from this module should be methods of Grid (or other structs)
 //       rather than free function. (For now they have just been converted from WGSL to Slang.)
@@ -323,6 +311,11 @@ public struct Node {
     /// (depending on the context). The fourth component contains the cell’s mass.
     // TODO: maybe we don’t really need to pack ourself.
     public vector<float, DIM + 1> momentum_velocity_mass;
+#if SLOSH_PML
+    /// Per-axis mass, `sum_p mass_scale_j * m_p * N_Ip`, which the momentum update divides by.
+    /// Equal to the scalar mass above unless a material rescales its own inertia (the PML does).
+    public vector<float, DIM> directional_mass;
+#endif
 #if SLOSH_CPIC
     // If a particle is incompatible with this node (as pepr CPIC’s concept of compatibility
     // based on affinities), it will contribute to this field instead of `momentum_velocity_mass`.
@@ -470,6 +463,9 @@ func reset(
 //   for (var i = invocation_id.x; i < num_nodes; i += num_threads) {
     if (i < num_nodes) {
        nodes[i].momentum_velocity_mass = vector<float, DIM + 1>(0.0);
+#if SLOSH_PML
+       nodes[i].directional_mass = vector<float, DIM>(0.0);
+#endif
 #if SLOSH_CPIC
        nodes[i].momentum_velocity_mass_incompatible = vector<float, DIM + 1>(0.0);
        nodes[i].cdf = NodeCdf(0.0, 0, NONE);
diff --git a/shaders/slosh/grid/sort.slang b/shaders/slosh/grid/sort.slang
index c25389d..815df76 100644
--- a/shaders/slosh/grid/sort.slang
+++ b/shaders/slosh/grid/sort.slang
@@ -3,11 +3,7 @@ module sort;
 import slosh.grid.grid;
 import slosh.solver.particle;
 
-// SLOSH_NODE_PARTICLE_LISTS gates the per-node particle linked-list build (see grid.slang).
-// Default on so a standalone compile keeps the full behavior.
-#ifndef SLOSH_NODE_PARTICLE_LISTS
-#define SLOSH_NODE_PARTICLE_LISTS 1
-#endif
+#include "slosh/features.h"
 
 // Returns the within-block sort bucket for a particle counted/inserted into its primary
 // block: one bucket per associated-cell slab along the slowest-varying node axis (y in 2D,
diff --git a/shaders/slosh/models/default.slang b/shaders/slosh/models/default.slang
index 839c0e4..d7fd122 100644
--- a/shaders/slosh/models/default.slang
+++ b/shaders/slosh/models/default.slang
@@ -1,10 +1,16 @@
 module default_model;
 
+#include "slosh/features.h"
+
 import slosh.models.interfaces;
 import slosh.models.linear_elasticity;
 import slosh.models.neo_hookean_elasticity;
 import slosh.models.drucker_prager;
+#if SLOSH_PML
+import slosh.models.pml;
+#endif
 import slosh.solver.particle;
+import slosh.models.utils;
 import slosh.aliases;
 
 
@@ -13,8 +19,16 @@ public enum class DefaultParticleModelType {
     ElasticNeoHookean = 1,
     SandLinear = 2,
     SandNeoHookean = 3,
+#if SLOSH_PML
+    // Absorbing (perfectly-matched-layer) particle; see models/pml.slang.
+    AbsorbingPml = 4,
+#endif
     // Provided to help downstream users extend with their own enum.
+#if SLOSH_PML
+    LEN = 5,
+#else
     LEN = 4,
+#endif
 }
 
 public struct SandModel<ElasticModel> {
@@ -60,12 +74,18 @@ public struct DefaultParticleModel<let BYTES_STRIDE: uint = DEFAULT_MODEL_BYTES_
         switch (DefaultParticleModelType(model_tag)) {
             case DefaultParticleModelType::ElasticLinear: {
                 let model = models.Load<LinearElasticModel>(model_offset);
-                let stress = model.kirchoff_stress(def_grad);
+                var stress = model.kirchoff_stress(def_grad);
+            #if SLOSH_PML
+                stress += viscous_stress(model.lambda, model.mu, data.velocity_gradient, data.stiffness_damping);
+            #endif
                 return ModelUpdateResult(stress);
             }
             case DefaultParticleModelType::ElasticNeoHookean: {
                 let model = models.Load<NeoHookeanModel>(model_offset);
-                let stress = model.kirchoff_stress(def_grad);
+                var stress = model.kirchoff_stress(def_grad);
+            #if SLOSH_PML
+                stress += viscous_stress(model.lambda, model.mu, data.velocity_gradient, data.stiffness_damping);
+            #endif
                 return ModelUpdateResult(stress);
             }
             case DefaultParticleModelType::SandLinear: {
@@ -73,7 +93,10 @@ public struct DefaultParticleModel<let BYTES_STRIDE: uint = DEFAULT_MODEL_BYTES_
                 let projection = model.plastic.project(model.plastic_state, def_grad);
                 let state_offset = model_offset + model.state_offset();
                 models.Store<DruckerPragerPlasticState>(state_offset, projection.state);
-                let stress = model.elastic.kirchoff_stress(projection.deformation_gradient);
+                var stress = model.elastic.kirchoff_stress(projection.deformation_gradient);
+            #if SLOSH_PML
+                stress += viscous_stress(model.elastic.lambda, model.elastic.mu, data.velocity_gradient, data.stiffness_damping);
+            #endif
                 def_grad = projection.deformation_gradient;
                 return ModelUpdateResult(stress);
             }
@@ -82,10 +105,21 @@ public struct DefaultParticleModel<let BYTES_STRIDE: uint = DEFAULT_MODEL_BYTES_
                 let projection = model.plastic.project(model.plastic_state, def_grad);
                 let state_offset = model_offset + model.state_offset();
                 models.Store<DruckerPragerPlasticState>(state_offset, projection.state);
-                let stress = model.elastic.kirchoff_stress(projection.deformation_gradient);
+                var stress = model.elastic.kirchoff_stress(projection.deformation_gradient);
+            #if SLOSH_PML
+                stress += viscous_stress(model.elastic.lambda, model.elastic.mu, data.velocity_gradient, data.stiffness_damping);
+            #endif
                 def_grad = projection.deformation_gradient;
                 return ModelUpdateResult(stress);
             }
+        #if SLOSH_PML
+            case DefaultParticleModelType::AbsorbingPml: {
+                let model = models.Load<PmlModel>(model_offset);
+                let stress = model.kirchoff_stress(def_grad)
+                    + model.viscous_stress(data.velocity_gradient, data.stiffness_damping);
+                return ModelUpdateResult(stress, -1.0, model.mass_scale());
+            }
+        #endif
             default:
                 return ModelUpdateResult(Mat<float>(0.0));
         }
@@ -121,6 +155,12 @@ public struct DefaultParticleModel<let BYTES_STRIDE: uint = DEFAULT_MODEL_BYTES_
                 let model = models.Load<SandModel<NeoHookeanModel>>(model_offset);
                 return model.elastic.timestep_bound(particle_density0, particle_velocity, def_grad_det, 1.0, cell_width);
             }
+        #if SLOSH_PML
+            case DefaultParticleModelType::AbsorbingPml: {
+                let model = models.Load<PmlModel>(model_offset);
+                return model.timestep_bound(particle_density0, particle_velocity, def_grad_det, cell_width);
+            }
+        #endif
             default:
                 return 0.0;
         }
diff --git a/shaders/slosh/models/interfaces.slang b/shaders/slosh/models/interfaces.slang
index 266c535..39432c0 100644
--- a/shaders/slosh/models/interfaces.slang
+++ b/shaders/slosh/models/interfaces.slang
@@ -1,5 +1,7 @@
 module interfaces;
 
+#include "slosh/features.h"
+
 import slosh.aliases;
 import slosh.solver.particle;
 
@@ -7,16 +9,35 @@ public struct ModelUpdateResult {
     public Mat<float> kirchoff_stress;
     // Negative means unchanged. If >= 0, the model modified init_volume (e.g., swell).
     public float init_volume;
+#if SLOSH_PML
+    // Per-axis multiplier on this material's inertia; one for anything ordinary. P2G scatters it
+    // into the node's directional mass, so a heavier axis also falls slower along it.
+    public Vect<float> mass_scale;
+#endif
 
     public __init(Mat stress) {
         this.kirchoff_stress = stress;
         this.init_volume = -1.0;
+#if SLOSH_PML
+        this.mass_scale = Vect<float>(1.0);
+#endif
     }
 
     public __init(Mat stress, float init_volume) {
         this.kirchoff_stress = stress;
         this.init_volume = init_volume;
+#if SLOSH_PML
+        this.mass_scale = Vect<float>(1.0);
+#endif
     }
+
+#if SLOSH_PML
+    public __init(Mat stress, float init_volume, Vect<float> mass_scale) {
+        this.kirchoff_stress = stress;
+        this.init_volume = init_volume;
+        this.mass_scale = mass_scale;
+    }
+#endif
 }
 
 public enum class ModelFlags {
@@ -31,6 +52,10 @@ public struct ParticleUpdateData {
     public matrix<float, DIM, DIM> velocity_gradient;
     public float init_volume;
     public float phase;
+#if SLOSH_PML
+    // Rayleigh stiffness-proportional damping coefficient, in seconds. See `viscous_stress`.
+    public float stiffness_damping;
+#endif
 
     public __init(float dt, float cell_width, uint particle_id, matrix<float, DIM, DIM> velocity_gradient, float init_volume, float phase) {
         this.dt = dt;
@@ -39,7 +64,22 @@ public struct ParticleUpdateData {
         this.velocity_gradient = velocity_gradient;
         this.init_volume = init_volume;
         this.phase = phase;
+#if SLOSH_PML
+        this.stiffness_damping = 0.0;
+#endif
+    }
+
+#if SLOSH_PML
+    public __init(float dt, float cell_width, uint particle_id, matrix<float, DIM, DIM> velocity_gradient, float init_volume, float phase, float stiffness_damping) {
+        this.dt = dt;
+        this.cell_width = cell_width;
+        this.particle_id = particle_id;
+        this.velocity_gradient = velocity_gradient;
+        this.init_volume = init_volume;
+        this.phase = phase;
+        this.stiffness_damping = stiffness_damping;
     }
+#endif
 }
 
 // TODO: keep separate traits for constitutive, plastic, fracture models,
diff --git a/shaders/slosh/models/pml.slang b/shaders/slosh/models/pml.slang
new file mode 100644
index 0000000..22cc105
--- /dev/null
+++ b/shaders/slosh/models/pml.slang
@@ -0,0 +1,102 @@
+//! Perfectly-matched-layer absorbing particle, after Kurima, Chandra & Soga, "Absorbing boundary
+//! conditions in material point method adopting perfectly matched layer theory" (arXiv:2407.02790).
+//!
+//! The real coordinate stretch only slows the layer-normal wave; the energy is removed by damping
+//! applied within the layer (`ParticleDynamics::damping`).
+
+module pml;
+
+#include "slosh/features.h"
+
+#if SLOSH_PML
+
+import slosh.aliases;
+import slosh.models.interfaces;
+import slosh.models.linear_elasticity;
+import stensor.utils.mat;
+
+public struct PmlModel {
+    public LinearElasticModel elastic;
+    // Per-axis coordinate stretching C'_j, zero outside the absorbing layer. The third entry is
+    // unused in 2D, kept so the layout matches the host side and avoids vector alignment padding.
+    public float stretch_x;
+    public float stretch_y;
+    public float stretch_z;
+
+    public func stretch() -> Vect<float> {
+    #if DIM == 2
+        return float2(this.stretch_x, this.stretch_y);
+    #else
+        return float3(this.stretch_x, this.stretch_y, this.stretch_z);
+    #endif
+    }
+
+    // Kirchhoff stress of the PML medium, transposed (the result is not symmetric, and p2g
+    // applies `M^T`). Implements the paper's Eq. 11; the inertia scaling is left to `mass_scale`.
+    public func kirchoff_stress(def_grad: Mat<float>) -> Mat<float> {
+        // `grad_u[k][j]` is du_j/dx_k. Small strain, like the paper: the layer stays near-rigid.
+        return this.stress_from_gradient(def_grad - diag(Vect<float>(1.0)));
+    }
+
+    // Rayleigh stiffness-proportional damping stress: the elastic operator applied to the strain
+    // rate. Blind to rigid-body motion, like `viscous_stress` in models/utils.slang.
+    public func viscous_stress(velocity_gradient: Mat<float>, stiffness_damping: float) -> Mat<float> {
+        if (stiffness_damping <= 0.0) {
+            return Mat<float>(0.0);
+        }
+
+        return this.stress_from_gradient(velocity_gradient * stiffness_damping);
+    }
+
+    func stress_from_gradient(grad_u: Mat<float>) -> Mat<float> {
+        let s = Vect<float>(1.0) + this.stretch();
+        let lambda = this.elastic.lambda;
+        let mu = this.elastic.mu;
+
+        var trace = 0.0;
+        [[ForceUnroll]]
+        for (var k = 0; k < DIM; k++) {
+            trace += grad_u[k][k];
+        }
+
+        var result = Mat<float>(0.0);
+
+        [[ForceUnroll]]
+        for (var a = 0; a < DIM; a++) {
+            // result[a][a] = sig_aa.
+            result[a][a] = (lambda + 2.0 * mu) * grad_u[a][a]
+                + s[a] * lambda * (trace - grad_u[a][a]);
+
+            [[ForceUnroll]]
+            for (var b = 0; b < DIM; b++) {
+                if (a != b) {
+                    // result[a][b] = sig_ba, with g_ab = grad_u[b][a] and g_ba = grad_u[a][b].
+                    result[a][b] = mu * (s[b] * grad_u[b][a] + s[b] * s[b] * grad_u[a][b]);
+                }
+            }
+        }
+
+        return result;
+    }
+
+    // The layer's direction-dependent inertia, `rho s_j^2` (the paper's Eq. 16). Gravity keeps
+    // acting on the real mass, so the stretched layer falls slower instead of getting heavier.
+    public func mass_scale() -> Vect<float> {
+        let s = Vect<float>(1.0) + this.stretch();
+        return s * s;
+    }
+
+    public func timestep_bound(
+        particle_density0: float,
+        particle_velocity: Vect<float>,
+        def_grad_det: float,
+        cell_width: float,
+    ) -> float {
+        // The stretch scales inertia at least as fast as stiffness, so every PML wave speed is
+        // bounded by the unstretched one and the elastic bound is conservative here.
+        return this.elastic.timestep_bound(
+            particle_density0, particle_velocity, def_grad_det, 1.0, cell_width);
+    }
+}
+
+#endif
diff --git a/shaders/slosh/models/utils.slang b/shaders/slosh/models/utils.slang
index a8a399a..2d11a5c 100644
--- a/shaders/slosh/models/utils.slang
+++ b/shaders/slosh/models/utils.slang
@@ -1,5 +1,7 @@
 module utils;
 
+#include "slosh/features.h"
+
 import stensor.utils.mat;
 import slosh.aliases;
 
@@ -133,4 +135,31 @@ public struct ElasticitySoundSpeedTimestepBound {
         let max_speed = max(length(velocity), sound_speed);
         return this.alpha * cell_width / max_speed;
     }
-}
\ No newline at end of file
+}
+
+#if SLOSH_PML
+// Rayleigh stiffness-proportional damping stress: `a_K * C : sym(grad v)`, with `C` the isotropic
+// elastic tensor. Unlike the mass-proportional damping in particle_update, it responds only to the
+// strain rate, so rigid-body motion is left alone.
+public func viscous_stress(
+    float lambda,
+    float mu,
+    Mat<float> velocity_gradient,
+    float stiffness_damping,
+) -> Mat<float> {
+    if (stiffness_damping <= 0.0) {
+        return Mat<float>(0.0);
+    }
+
+    // Rate of deformation, the symmetric part of the velocity gradient.
+    let rate = (velocity_gradient + transpose(velocity_gradient)) * 0.5;
+
+    var trace = 0.0;
+    [[ForceUnroll]]
+    for (var k = 0; k < DIM; k++) {
+        trace += rate[k][k];
+    }
+
+    return (rate * (2.0 * mu) + diag(Vect<float>(lambda * trace))) * stiffness_damping;
+}
+#endif
diff --git a/shaders/slosh/solver/boundary_condition.slang b/shaders/slosh/solver/boundary_condition.slang
index ea3db87..60d692b 100644
--- a/shaders/slosh/solver/boundary_condition.slang
+++ b/shaders/slosh/solver/boundary_condition.slang
@@ -1,6 +1,14 @@
 module boundary_condition;
 
 import slosh.aliases;
+import slosh.collision.collide;
+
+// Depth, in cells, of the `NonReflecting` absorbing band, graded linearly from full strength at
+// the collider surface. A deep, gentle band absorbs much better than a shallow, aggressive one,
+// and grading makes absorption nearly insensitive to how the surface aligns with the node lattice.
+//
+// Must not exceed `COLLISION_REPORT_CELLS`, past which no collision is reported at all.
+static const float ABSORBING_LAYERS = 6.0;
 
 public enum BoundaryConditionType {
     // Stick condition.
@@ -22,13 +30,26 @@ public enum BoundaryConditionType {
 public struct BoundaryCondition {
     BoundaryConditionType type;
     float friction;
+    // Pressure and shear wave speeds of the material in contact with the boundary.
+    // Only read by the NonReflecting boundary condition.
+    float wave_speed_p;
+    float wave_speed_s;
 
     public __init(BoundaryConditionType type, float friction) {
         this.type = type;
         this.friction = friction;
+        this.wave_speed_p = 0.0;
+        this.wave_speed_s = 0.0;
+    }
+
+    public __init(BoundaryConditionType type, float friction, float wave_speed_p, float wave_speed_s) {
+        this.type = type;
+        this.friction = friction;
+        this.wave_speed_p = wave_speed_p;
+        this.wave_speed_s = wave_speed_s;
     }
 
-    public func project_velocity(vel: Vect, n: Vect, float dist, float dt) -> Vect {
+    public func project_velocity(vel: Vect, n: Vect, float dist, float dt, float cell_width) -> Vect {
         switch (this.type) {
             case BoundaryConditionType::Stick: {
                 if (dist <= 0.0) {
@@ -49,37 +70,19 @@ public struct BoundaryCondition {
             }
             case BoundaryConditionType::FrictionZUp:
             case BoundaryConditionType::Separate: {
-                let normal_vel = dot(vel, n);
-
-                if (normal_vel < 0.0) {
-                    if (dist <= 0.0) {
-                        #if DIM == 3
-                        let apply_friction = this.type == BoundaryConditionType::Separate ||
-                            (this.type == BoundaryConditionType::FrictionZUp && n.z >= 0.0);
-                        #else
-                        let apply_friction = this.type == BoundaryConditionType::Separate ||
-                            (this.type == BoundaryConditionType::FrictionZUp && n.y >= 0.0);
-                        #endif
-                        var friction = this.friction;
-                        if (!apply_friction) {
-                            friction = 0.0;
-                        }
-                        let tangent_vel = vel - n * normal_vel;
-                        let tangent_vel_len = length(tangent_vel);
-                        let tangent_vel_dir = select(tangent_vel_len > 1.0e-8, tangent_vel / tangent_vel_len, Vect<float>(0.0));
-                        return tangent_vel_dir * max(0.0, tangent_vel_len + this.friction * normal_vel);
-                    } else if (-normal_vel * dt > dist) {
-                        // Predictive
-                       return vel - (dist / dt + normal_vel) * n;
-                   } else {
-                       return vel;
-                   }
-                } else {
-                    return vel;
-                }
+                #if DIM == 3
+                let apply_friction = this.type == BoundaryConditionType::Separate || n.z >= 0.0;
+                #else
+                let apply_friction = this.type == BoundaryConditionType::Separate || n.y >= 0.0;
+                #endif
+                return this.separate_velocity(vel, n, dist, dt, cell_width, apply_friction, false);
             }
             case BoundaryConditionType::NonReflecting: {
-                return vel;
+                // Absorb first, then enforce contact, so the dashpot traction cannot push material
+                // back through a surface the contact just cleared. `friction` is read as in
+                // `Separate`, so zero leaves a pure absorbing wall.
+                let absorbed = this.absorb_velocity(vel, n, dist, dt, cell_width);
+                return this.separate_velocity(absorbed, n, dist, dt, cell_width, true, true);
             }
             case BoundaryConditionType::Disabled:
             default: {
@@ -87,7 +90,60 @@ public struct BoundaryCondition {
             }
         }
     }
-}
 
+    // One-sided non-penetration with Coulomb friction: the `Separate` response, factored out so
+    // `NonReflecting` can layer it on top of the dashpot.
+    func separate_velocity(vel: Vect, n: Vect, float dist, float dt, float cell_width, bool apply_friction, bool keep_normal) -> Vect {
+        let normal_vel = dot(vel, n);
 
+        if (normal_vel >= 0.0) {
+            return vel;
+        }
+
+        if (dist <= 0.0) {
+            var friction = this.friction;
+            var normal_coeff = 0.0;
+            if (keep_normal) {
+                normal_coeff = 1.0;
+            }
+            if (!apply_friction) {
+                friction = 0.0;
+            }
+            let tangent_vel = vel - n * normal_vel;
+            let tangent_vel_len = length(tangent_vel);
+            let tangent_vel_dir = select(tangent_vel_len > 1.0e-8, tangent_vel / tangent_vel_len, Vect<float>(0.0));
+            return n * (normal_vel * normal_coeff) + tangent_vel_dir * max(0.0, tangent_vel_len + friction * normal_vel);
+        } else if (dist <= cell_width * 0.5 && -normal_vel * dt > dist) {
+            // Predictive. The `dist` bound keeps this to half a cell, so the deeper reporting
+            // band `NonReflecting` needs (see `COLLISION_REPORT_CELLS`) does not extend it.
+            return vel - (dist / dt + normal_vel) * n;
+        } else {
+            return vel;
+        }
+    }
 
+    // Lysmer-Kuhlemeyer viscous dashpots: the boundary applies the traction a semi-infinite
+    // continuation of the material would apply, `t_n = -rho c_p v_n` and `t_t = -rho c_s v_t`.
+    //
+    // The node mass cancels out of the resulting velocity change, leaving the wave's CFL number
+    // `c * dt / h` as the damping factor; the clamp only guards against inconsistent wave speeds.
+    //
+    // A zero `wave_speed_s` absorbs the pressure wave only: the shear dashpot is indiscriminate
+    // and would drag any material resting inside the band to a halt.
+    func absorb_velocity(vel: Vect, n: Vect, float dist, float dt, float cell_width) -> Vect {
+        // `dist` is the distance to the surface minus one cell. Full strength on and inside the
+        // surface, fading out over `ABSORBING_LAYERS` cells.
+        let depth = dist + cell_width;
+        let weight = clamp(1.0 - depth / (ABSORBING_LAYERS * cell_width), 0.0, 1.0);
+
+        if (weight <= 0.0) {
+            return vel;
+        }
+
+        let normal_vel = dot(vel, n);
+        let tangent_vel = vel - n * normal_vel;
+        let damping_p = clamp(weight * this.wave_speed_p * dt / cell_width, 0.0, 1.0);
+        let damping_s = clamp(weight * this.wave_speed_s * dt / cell_width, 0.0, 1.0);
+        return n * (normal_vel * (1.0 - damping_p)) + tangent_vel * (1.0 - damping_s);
+    }
+}
diff --git a/shaders/slosh/solver/grid_update.slang b/shaders/slosh/solver/grid_update.slang
index d818cb5..86533ac 100644
--- a/shaders/slosh/solver/grid_update.slang
+++ b/shaders/slosh/solver/grid_update.slang
@@ -7,12 +7,19 @@ import slosh.solver.params;
 import slosh.collision.collide;
 import slosh.solver.boundary_condition;
 
-// SLOSH_CPIC gates the incompatible-momentum lane on Node (see grid.slang). Default on so a
-// standalone compile matches the crate's default layout.
-#ifndef SLOSH_CPIC
-#define SLOSH_CPIC 1
+#include "slosh/features.h"
+
+#if SLOSH_PML
+typealias NodeMass = Vect<float>;
+#else
+typealias NodeMass = float;
 #endif
 
+// Reciprocal node mass, zero where the node has none.
+func inv_mass_of(mass: NodeMass) -> NodeMass {
+    return select(mass > NodeMass(0.0), NodeMass(1.0) / mass, NodeMass(0.0));
+}
+
 #if DIM == 2
 static const uint WORKGROUP_SIZE_X = 8;
 static const uint WORKGROUP_SIZE_Y = 8;
@@ -79,10 +86,18 @@ func grid_update(
      * Step 1: momentum -> velocity update (gravity + clamping).
      */
     let momentum_velocity_mass = nodes[global_id].momentum_velocity_mass;
-    var velocity_mass = update_single_cell(sim_params, cell_width, cell_pos, momentum_velocity_mass);
+#if SLOSH_PML
+    let inv_mass = inv_mass_of(nodes[global_id].directional_mass);
+#else
+    let inv_mass = inv_mass_of(momentum_velocity_mass[DIM]);
+#endif
+    var velocity_mass = update_single_cell(sim_params, cell_width, cell_pos, momentum_velocity_mass, inv_mass);
 #if SLOSH_CPIC
+    // The incompatible lane carries no directional mass of its own: CPIC is not meant to overlap
+    // with the materials that need one.
     let momentum_velocity_mass_incompatible = nodes[global_id].momentum_velocity_mass_incompatible;
-    let velocity_mass_incompatible = update_single_cell(sim_params, cell_width, cell_pos, momentum_velocity_mass_incompatible);
+    let inv_mass_incompatible = inv_mass_of(NodeMass(momentum_velocity_mass_incompatible[DIM]));
+    let velocity_mass_incompatible = update_single_cell(sim_params, cell_width, cell_pos, momentum_velocity_mass_incompatible, inv_mass_incompatible);
     nodes[global_id].momentum_velocity_mass_incompatible = velocity_mass_incompatible;
 #endif
 
@@ -127,12 +142,12 @@ func grid_update(
         #if DIM == 2
             let node_velocity = velocity_mass.xy;
             let node_mass = velocity_mass.z;
-            let projected = material.project_velocity(node_velocity, collision.normal, dist_with_margin, sim_params.dt);
+            let projected = material.project_velocity(node_velocity, collision.normal, dist_with_margin, sim_params.dt, cell_width);
             velocity_mass = float3(projected, node_mass);
         #else
             let node_velocity = velocity_mass.xyz;
             let node_mass = velocity_mass.w;
-            let projected = material.project_velocity(node_velocity, collision.normal, dist_with_margin, sim_params.dt);
+            let projected = material.project_velocity(node_velocity, collision.normal, dist_with_margin, sim_params.dt, cell_width);
             velocity_mass = float4(projected, node_mass);
         #endif
     }
@@ -141,9 +156,10 @@ func grid_update(
 }
 
 #if DIM == 2
-func update_single_cell(sim_params: SimulationParams, cell_width: float, cell_pos: float2, momentum_velocity_mass: float3) -> float3 {
+func update_single_cell(sim_params: SimulationParams, cell_width: float, cell_pos: float2, momentum_velocity_mass: float3, inv_mass: NodeMass) -> float3 {
+    // Gravity stays proportional to the real mass, so a material that rescales its own inertia
+    // (the PML) falls at `g / scale` instead of carrying that much weight.
     let mass = momentum_velocity_mass.z;
-    let inv_mass = select(mass > 0.0, 1.0 / mass, 0.0);
     var velocity = (momentum_velocity_mass.xy + mass * sim_params.gravity * sim_params.dt) * inv_mass;
     // Clamp the velocity so it doesn't exceed 1 grid cell in one step.
     let vel_limit = float2(cell_width / sim_params.dt);
@@ -151,9 +167,9 @@ func update_single_cell(sim_params: SimulationParams, cell_width: float, cell_po
     return float3(velocity, mass);
 }
 #else
-func update_single_cell(sim_params: SimulationParams, cell_width: float, cell_pos: float3, momentum_velocity_mass: float4) -> float4 {
+func update_single_cell(sim_params: SimulationParams, cell_width: float, cell_pos: float3, momentum_velocity_mass: float4, inv_mass: NodeMass) -> float4 {
+    // See the 2D overload above.
     let mass = momentum_velocity_mass.w;
-    let inv_mass = select(mass > 0.0, 1.0 / mass, 0.0);
     var velocity = (momentum_velocity_mass.xyz + mass * sim_params.gravity * sim_params.dt) * inv_mass;
 
     // Clamp the velocity so it doesn't exceed 1 grid cell in one step.
diff --git a/shaders/slosh/solver/p2g.slang b/shaders/slosh/solver/p2g.slang
index cdab0a6..88f4f64 100644
--- a/shaders/slosh/solver/p2g.slang
+++ b/shaders/slosh/solver/p2g.slang
@@ -1,5 +1,7 @@
 module p2g;
 
+#include "slosh/features.h"
+
 import slosh.solver.params;
 import slosh.solver.particle;
 import slosh.solver.boundary_condition;
@@ -17,6 +19,9 @@ static const uint NUM_SHARED_CELLS = 10 * 10; // block-size plus 2 from adjacent
 groupshared float3 shared_vel_mass[NUM_SHARED_CELLS];
 groupshared float2x2 shared_affine[NUM_SHARED_CELLS];
 groupshared float2 shared_force_dt[NUM_SHARED_CELLS];
+    #if SLOSH_PML
+groupshared float2 shared_mass_scale[NUM_SHARED_CELLS];
+    #endif
 #else
 static const uint WORKGROUP_SIZE_X = 4;
 static const uint WORKGROUP_SIZE_Y = 4;
@@ -25,6 +30,9 @@ static const uint NUM_SHARED_CELLS = 6 * 6 * 6; // block-size plus 2 from adjace
 groupshared float4 shared_vel_mass[NUM_SHARED_CELLS];
 groupshared float3x3 shared_affine[NUM_SHARED_CELLS];
 groupshared float3 shared_force_dt[NUM_SHARED_CELLS];
+    #if SLOSH_PML
+groupshared float3 shared_mass_scale[NUM_SHARED_CELLS];
+    #endif
 #endif
 groupshared SharedNode shared_nodes[NUM_SHARED_CELLS];
 groupshared Position shared_pos[NUM_SHARED_CELLS];
@@ -44,8 +52,14 @@ struct SharedNode {
 struct P2GStepResult {
 #if DIM == 2
     float3 new_momentum_velocity_mass;
+    #if SLOSH_PML
+    float2 new_directional_mass;
+    #endif
 #else
     float4 new_momentum_velocity_mass;
+    #if SLOSH_PML
+    float3 new_directional_mass;
+    #endif
 #endif
 }
 
@@ -111,6 +125,9 @@ func p2g(
 
         let partial_result = p2g_step(packed_cell_index_in_block, grid[0].cell_width); // , node_affinities, collider_id);
         total_result.new_momentum_velocity_mass += partial_result.new_momentum_velocity_mass;
+#if SLOSH_PML
+        total_result.new_directional_mass += partial_result.new_directional_mass;
+#endif
     }
 
     // Grid update.
@@ -128,6 +145,9 @@ func p2g(
 
     // Write the node state to global memory.
     nodes[global_id].momentum_velocity_mass = total_result.new_momentum_velocity_mass;
+#if SLOSH_PML
+    nodes[global_id].directional_mass = total_result.new_directional_mass;
+#endif
 }
 
 func p2g_step(
@@ -145,9 +165,15 @@ func p2g_step(
 #if DIM == 2
     let bottommost_contributing_node = flatten_shared_shift(2u, 2u);
     var new_momentum_velocity_mass = float3(0.0);
+    #if SLOSH_PML
+    var new_directional_mass = float2(0.0);
+    #endif
 #else
     let bottommost_contributing_node = flatten_shared_shift(2u, 2u, 2u);
     var new_momentum_velocity_mass = float4(0.0);
+    #if SLOSH_PML
+    var new_directional_mass = float3(0.0);
+    #endif
 #endif
     for (var i = 0u; i < NBH_LEN; i += 1u) {
         let packed_shift = NBH_SHIFTS_SHARED[i];
@@ -164,14 +190,22 @@ func p2g_step(
         let particle_vel = particle_vel_mass.xy;
         let particle_mass = particle_vel_mass.z;
         let shift = uint2(2u, 2) - NBH_SHIFTS[i];
+    #if SLOSH_PML
+        let momentum = shared_mass_scale[nbh_shared_index] * particle_vel * particle_mass + particle_force_dt;
+    #else
         let momentum = particle_vel * particle_mass + particle_force_dt;
+    #endif
         let dpt = ref_elt_pos_minus_particle_pos + float2(shift) * cell_width; // cell_pos - particle_pos
         let weight = w[0][shift.x] * w[1][shift.y];
 #else
         let particle_vel = particle_vel_mass.xyz;
         let particle_mass = particle_vel_mass.w;
         let shift = uint3(2u, 2, 2) - NBH_SHIFTS[i];
+    #if SLOSH_PML
+        let momentum = shared_mass_scale[nbh_shared_index] * particle_vel * particle_mass + particle_force_dt;
+    #else
         let momentum = particle_vel * particle_mass + particle_force_dt;
+    #endif
         let dpt = ref_elt_pos_minus_particle_pos + float3(shift) * cell_width; // cell_pos - particle_pos
         let weight = w[0][shift.x] * w[1][shift.y] * w[2][shift.z];
 #endif
@@ -179,9 +213,16 @@ func p2g_step(
 //        let particle_affinity = shared_affinities[nbh_shared_index];
         let contribution = vector<float, DIM + 1>(mul(dpt, particle_affine) + momentum, particle_mass) * weight;
         new_momentum_velocity_mass += contribution;
+#if SLOSH_PML
+        new_directional_mass += shared_mass_scale[nbh_shared_index] * particle_mass * weight;
+#endif
     }
 
+#if SLOSH_PML
+    return P2GStepResult(new_momentum_velocity_mass, new_directional_mass);
+#else
     return P2GStepResult(new_momentum_velocity_mass);
+#endif
 }
 
 #if DIM == 2
@@ -340,6 +381,9 @@ func fetch_next_particle(
                     shared_pos[shared_flat_index] = particles_pos[curr_particle_id];
                     shared_affine[shared_flat_index] = particles_kin[curr_particle_id].affine;
                     shared_force_dt[shared_flat_index] = particles_kin[curr_particle_id].force_dt;
+#if SLOSH_PML
+                    shared_mass_scale[shared_flat_index] = particles_kin[curr_particle_id].mass_scale;
+#endif
 
 #if DIM == 2
                     shared_vel_mass[shared_flat_index] = float3(particles_kin[curr_particle_id].velocity, particles_kin[curr_particle_id].mass);
@@ -357,11 +401,17 @@ func fetch_next_particle(
                     shared_affine[shared_flat_index] = float2x2(float2(0.0), float2(0.0));
                     shared_vel_mass[shared_flat_index] = float3(0.0);
                     shared_force_dt[shared_flat_index] = float2(0.0);
+    #if SLOSH_PML
+                    shared_mass_scale[shared_flat_index] = float2(0.0);
+    #endif
 #else
                     shared_pos[shared_flat_index].pt = float3(0.0);
                     shared_affine[shared_flat_index] = float3x3(float3(0.0), float3(0.0), float3(0.0));
                     shared_vel_mass[shared_flat_index] = float4(0.0);
                     shared_force_dt[shared_flat_index] = float3(0.0);
+    #if SLOSH_PML
+                    shared_mass_scale[shared_flat_index] = float3(0.0);
+    #endif
 #endif
                 }
 
diff --git a/shaders/slosh/solver/p2g_scatter_style.slang b/shaders/slosh/solver/p2g_scatter_style.slang
index d11f4db..a765969 100644
--- a/shaders/slosh/solver/p2g_scatter_style.slang
+++ b/shaders/slosh/solver/p2g_scatter_style.slang
@@ -1,5 +1,7 @@
 module p2g;
 
+#include "slosh/features.h"
+
 import slosh.solver.params;
 import slosh.solver.particle;
 import slosh.solver.boundary_condition;
@@ -26,10 +28,16 @@ static const uint WORKGROUP_SIZE = WORKGROUP_SIZE_X * WORKGROUP_SIZE_Y * WORKGRO
 #if DIM == 2
 groupshared float2 shared_pos[WORKGROUP_SIZE];
 groupshared float2 shared_momentum[WORKGROUP_SIZE];
+    #if SLOSH_PML
+groupshared float2 shared_dir_mass[WORKGROUP_SIZE];
+    #endif
 groupshared float2x2 shared_affine[WORKGROUP_SIZE];
 #else
 groupshared float3 shared_pos[WORKGROUP_SIZE];
 groupshared float3 shared_momentum[WORKGROUP_SIZE];
+    #if SLOSH_PML
+groupshared float3 shared_dir_mass[WORKGROUP_SIZE];
+    #endif
 groupshared float3x3 shared_affine[WORKGROUP_SIZE];
 #endif
 groupshared float shared_mass[WORKGROUP_SIZE];
@@ -71,12 +79,18 @@ func p2g_scatter_style(
     let local_cell = int2(int(tid % 8u), int(tid / 8u));
     let cell_pos = float2(block_vid * 8 + local_cell) * cell_width;
     var acc = float3(0.0);
+    #if SLOSH_PML
+    var dir_mass_acc = float2(0.0);
+    #endif
     // This thread's node slab along the sort axis (the slowest-varying node axis).
     let node_slab = local_cell.y;
 #else
     let local_cell = int3(int(tid % 4u), int((tid / 4u) % 4u), int(tid / 16u));
     let cell_pos = float3(block_vid * 4 + local_cell) * cell_width;
     var acc = float4(0.0);
+    #if SLOSH_PML
+    var dir_mass_acc = float3(0.0);
+    #endif
     let node_slab = local_cell.z;
 #endif
 
@@ -93,7 +107,14 @@ func p2g_scatter_style(
             shared_pos[tid] = pos;
             shared_mass[tid] = kin.mass;
             shared_affine[tid] = kin.affine;
+        #if SLOSH_PML
+            // Only the inertial momentum carries the per-axis mass scale; `force_dt` is an
+            // external force and, like gravity, is left alone.
+            shared_momentum[tid] = kin.mass_scale * kin.velocity * kin.mass + kin.force_dt;
+            shared_dir_mass[tid] = kin.mass_scale * kin.mass;
+        #else
             shared_momentum[tid] = kin.velocity * kin.mass + kin.force_dt;
+        #endif
 
             // Slab key along the sort axis, relative to this block; must match the bucket key
             // used by the sort (same associated-cell rounding, same clamp at -2) so the shared
@@ -144,6 +165,9 @@ func p2g_scatter_style(
             if (weight != 0.0) {
                 let momentum = mul(dpt, shared_affine[p]) + shared_momentum[p];
                 acc += vector<float, DIM + 1>(momentum, shared_mass[p]) * weight;
+            #if SLOSH_PML
+                dir_mass_acc += shared_dir_mass[p] * weight;
+            #endif
             }
         }
     }
@@ -152,4 +176,7 @@ func p2g_scatter_style(
     // (no atomics, no inter-block races), zeroing cells that received no contribution.
     let global_chunk_id = block_header_id_to_physical_id(BlockHeaderId(bid));
     nodes[global_chunk_id.id + tid].momentum_velocity_mass = acc;
+#if SLOSH_PML
+    nodes[global_chunk_id.id + tid].directional_mass = dir_mass_acc;
+#endif
 }
diff --git a/shaders/slosh/solver/particle.slang b/shaders/slosh/solver/particle.slang
index 2b1eb35..1f014fc 100644
--- a/shaders/slosh/solver/particle.slang
+++ b/shaders/slosh/solver/particle.slang
@@ -1,5 +1,7 @@
 module particle;
 
+#include "slosh/features.h"
+
 public struct Position {
     public vector<float, DIM> pt;
 }
@@ -21,6 +23,11 @@ public struct Kinematics {
     //       the force so we don't have to have access to dt in p2g
     //       (adding that access would exceed 10 bindings…)
     public vector<float, DIM> force_dt;
+#if SLOSH_PML
+    // Per-axis multiplier on this particle's inertia, one for ordinary materials (the PML uses
+    // `s_j^2`). Written by particle_update, so it stays one over the very first substep.
+    public vector<float, DIM> mass_scale;
+#endif
     public float vel_grad_det;
     public float mass;
     public uint enabled;
@@ -35,6 +42,9 @@ public struct ParticleProperties {
     public float init_volume;
     public float init_radius;
     public float damping;
+#if SLOSH_PML
+    public float stiffness_damping;
+#endif
     public float phase;
     public uint fixed;
 }
diff --git a/shaders/slosh/solver/particle_update.slang b/shaders/slosh/solver/particle_update.slang
index bbbc0f3..e692ed1 100644
--- a/shaders/slosh/solver/particle_update.slang
+++ b/shaders/slosh/solver/particle_update.slang
@@ -1,5 +1,7 @@
 module particle_update;
 
+#include "slosh/features.h"
+
 import slosh.solver.params;
 import slosh.solver.particle;
 import slosh.solver.boundary_condition;
@@ -123,7 +125,11 @@ func particle_update(
         /*
          * Constitutive model.
          */
+    #if SLOSH_PML
+        let update_data = ParticleUpdateData(dt, cell_width, particle_id, kin.affine, props.init_volume, props.phase, props.stiffness_damping);
+    #else
         let update_data = ParticleUpdateData(dt, cell_width, particle_id, kin.affine, props.init_volume, props.phase);
+    #endif
         let update_result = model.update(particles_model, update_data, def_grad);
 
         // Write back init_volume if the model modified it (e.g., swell).
@@ -137,7 +143,15 @@ func particle_update(
          */
         let inv_d = QuadraticKernel::inv_d(cell_width);
         // NOTE: the velocity gradient was stored in the affine buffer.
-        let affine = kin.affine * kin.mass - update_result.kirchoff_stress * (effective_init_volume * inv_d * dt);
+    #if SLOSH_PML
+        // P2G contracts this as `sum_i dpt_i * A_ij`, so scaling column `j` applies the per-axis
+        // inertia to the inertial half only; the stress half is an internal force and is left alone.
+        kin.mass_scale = update_result.mass_scale;
+        let inertia = mul(kin.affine * kin.mass, diag(update_result.mass_scale));
+    #else
+        let inertia = kin.affine * kin.mass;
+    #endif
+        let affine = inertia - update_result.kirchoff_stress * (effective_init_volume * inv_d * dt);
 
         /*
          * Write back the new particle properties.
diff --git a/shaders/slosh/solver/timestep_bound.slang b/shaders/slosh/solver/timestep_bound.slang
index 18809ce..c533528 100644
--- a/shaders/slosh/solver/timestep_bound.slang
+++ b/shaders/slosh/solver/timestep_bound.slang
@@ -1,5 +1,7 @@
 module timestep_bound;
 
+#include "slosh/features.h"
+
 import slosh.aliases;
 import slosh.grid.grid;
 import slosh.models.utils;
@@ -68,6 +70,20 @@ func estimate_timestep_bound(
         let model = ParticleModel();
         var dt = model.timestep_bound(particles_model, particle_id, density0, def_grad, velocity, cell_width);
 
+#if SLOSH_PML
+        /*
+         * Rayleigh stiffness-proportional damping tightens the explicit bound below the CFL one.
+         * For central differencing the critical step is `(2 / w) * (sqrt(1 + x^2) - x)` with
+         * `x = a_K w / 2`, using the model's own bound as the highest resolved frequency
+         * `w ~ 1 / dt`. Overshooting it silently disables diverging particles instead of blowing up.
+         */
+        let stiffness_damping = particles_props[particle_id].stiffness_damping;
+        if (stiffness_damping > 0.0 && dt > 0.0) {
+            let x = stiffness_damping / (2.0 * dt);
+            dt = min(dt, 2.0 * dt * (sqrt(1.0 + x * x) - x));
+        }
+#endif
+
         /*
          * Velocity-based restrictions (section 4.2)
          */
diff --git a/shaders_testbed/slosh_testbed/prep_readback2.slang b/shaders_testbed/slosh_testbed/prep_readback2.slang
index c1a622e..2e84782 100644
--- a/shaders_testbed/slosh_testbed/prep_readback2.slang
+++ b/shaders_testbed/slosh_testbed/prep_readback2.slang
@@ -27,6 +27,13 @@ static const uint CDF_NORMALS = 4;
 static const uint CDF_DISTANCES = 5;
 static const uint CDF_SIGNS = 6;
 
+// Principal strain that saturates the VOLUME render mode.
+static const float VOLUME_STRAIN_SCALE = 0.0015;
+// Colors the VOLUME render mode saturates to. Pure shear lights up both ramps, so it reads
+// as white.
+static const float3 VOLUME_COMPRESSION_COLOR = float3(1.0, 0.45, 0.1);
+static const float3 VOLUME_DILATION_COLOR = float3(0.1, 0.65, 1.0);
+
 [shader("compute")]
 [numthreads(64, 1, 1)]
 func prep_readback(
@@ -62,9 +69,17 @@ func prep_readback(
             let vel = particles_kin[particle_id].velocity;
             instances[particle_id].color = float4(abs(vel) * dt * 100.0 + float2(0.2), color.z, color.w);
         } else if (mode == VOLUME) {
+            // Principal strains, positive under compression and negative under dilation.
             let svd = Svd2(def_grad);
-            let color_xy = (float2(1.0) - svd.S) / 0.005 + float2(0.2);
-            instances[particle_id].color = float4(color_xy, color.z, color.w);
+            let strain = (float2(1.0) - svd.S) / VOLUME_STRAIN_SCALE;
+            // One ramp out of black per sign, rather than the two ends of a single diverging
+            // ramp: that keeps the full brightness range, and dilation visible.
+            let compression = saturate(max(strain.x, strain.y));
+            let dilation = saturate(-min(strain.x, strain.y));
+            // Square root: lifts the small strains without moving the saturation point.
+            let rgb = sqrt(compression) * VOLUME_COMPRESSION_COLOR
+                + sqrt(dilation) * VOLUME_DILATION_COLOR;
+            instances[particle_id].color = float4(rgb, color.w);
         } else if (mode == PHASE) {
             let phase = particles_props[particle_id].phase;
             instances[particle_id].color = float4(0.0, 0.4 * phase, 0.4 * (1.0 - phase), color.w);
diff --git a/shaders_testbed/slosh_testbed/prep_readback3.slang b/shaders_testbed/slosh_testbed/prep_readback3.slang
index c3527bc..0426702 100644
--- a/shaders_testbed/slosh_testbed/prep_readback3.slang
+++ b/shaders_testbed/slosh_testbed/prep_readback3.slang
@@ -27,6 +27,13 @@ static const uint CDF_NORMALS = 4;
 static const uint CDF_DISTANCES = 5;
 static const uint CDF_SIGNS = 6;
 
+// Principal strain that saturates the VOLUME render mode.
+static const float VOLUME_STRAIN_SCALE = 0.0015;
+// Colors the VOLUME render mode saturates to. Pure shear lights up both ramps, so it reads
+// as white.
+static const float3 VOLUME_COMPRESSION_COLOR = float3(1.0, 0.45, 0.1);
+static const float3 VOLUME_DILATION_COLOR = float3(0.1, 0.65, 1.0);
+
 [shader("compute")]
 [numthreads(64, 1, 1)]
 func prep_readback(
@@ -63,9 +70,17 @@ func prep_readback(
             let vel = particles_kin[particle_id].velocity;
             instances[particle_id].color = float4(abs(vel) * dt * 100.0 + float3(0.2), color.w);
         } else if (mode == VOLUME) {
+            // Principal strains, positive under compression and negative under dilation.
             let svd = Svd3(def_grad);
-            let color_xyz = (float3(1.0) - svd.S) / 0.005 + float3(0.2);
-            instances[particle_id].color = float4(color_xyz, color.w);
+            let strain = (float3(1.0) - svd.S) / VOLUME_STRAIN_SCALE;
+            // One ramp out of black per sign, rather than the two ends of a single diverging
+            // ramp: that keeps the full brightness range, and dilation visible.
+            let compression = saturate(max(max(strain.x, strain.y), strain.z));
+            let dilation = saturate(-min(min(strain.x, strain.y), strain.z));
+            // Square root: lifts the small strains without moving the saturation point.
+            let rgb = sqrt(compression) * VOLUME_COMPRESSION_COLOR
+                + sqrt(dilation) * VOLUME_DILATION_COLOR;
+            instances[particle_id].color = float4(rgb, color.w);
         } else if (mode == PHASE) {
             let phase = particles_props[particle_id].phase;
             instances[particle_id].color = float4(0.0, 0.4 * phase, 0.4 * (1.0 - phase), color.w);
diff --git a/src/grid/grid.rs b/src/grid/grid.rs
index 58bdcce..5b0d98f 100644
--- a/src/grid/grid.rs
+++ b/src/grid/grid.rs
@@ -258,6 +258,11 @@ pub struct GpuGridMetadata {
 #[repr(C)]
 pub struct GpuGridNode {
     momentum_velocity_mass: glam::Vec4,
+    // Per-axis mass; see `Node.directional_mass` in grid.slang.
+    #[cfg(all(feature = "pml", feature = "dim2"))]
+    directional_mass: glam::Vec2,
+    #[cfg(all(feature = "pml", feature = "dim3"))]
+    directional_mass: glam::Vec4,
     // Gated by the cpic feature (default on). Turning cpic off drops the node from 48 to 16
     // bytes, cutting grid traffic ~3x, and gives up collision-detection-field support. Keep in
     // lockstep with Node in grid.slang, which SLOSH_CPIC gates the same way.
diff --git a/src/lib.rs b/src/lib.rs
index e978134..92d4191 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -89,9 +89,10 @@ pub const SLANG_SRC_DIR: include_dir::Dir<'_> =
 pub fn register_shaders(compiler: &mut SlangCompiler) {
     stensor::register_shaders(compiler);
     compiler.add_dir(SLANG_SRC_DIR.clone());
-    // Mirror the cpic and node_particle_lists cargo features into shader macros so the shader
-    // Node layout and sort kernels match GpuGridNode and GridArgs on the Rust side. slosh sets
-    // these itself (unlike DIM, which the consumer picks) since they're crate-layout invariants.
+    // Mirror the cpic, node_particle_lists and pml cargo features into shader macros so the
+    // shader Node/particle layouts and sort kernels match GpuGridNode, Kinematics,
+    // ParticleProperties and GridArgs on the Rust side. slosh sets these itself (unlike DIM,
+    // which the consumer picks) since they're crate-layout invariants.
     compiler.set_global_macro("SLOSH_CPIC", if cfg!(feature = "cpic") { 1 } else { 0 });
     compiler.set_global_macro(
         "SLOSH_NODE_PARTICLE_LISTS",
@@ -101,6 +102,7 @@ pub fn register_shaders(compiler: &mut SlangCompiler) {
             0
         },
     );
+    compiler.set_global_macro("SLOSH_PML", if cfg!(feature = "pml") { 1 } else { 0 });
 }
 
 /// Mathematical types and utilities for physics simulation.
diff --git a/src/models/mod.rs b/src/models/mod.rs
index d5b77be..f306576 100644
--- a/src/models/mod.rs
+++ b/src/models/mod.rs
@@ -3,9 +3,13 @@
 //! This module provides material models that define how particles respond to deformation:
 //! - [`ElasticCoefficients`]: Linear elasticity using Lamé parameters
 //! - [`DruckerPrager`]: Drucker-Prager plasticity model for granular materials (sand, soil)
+//! - `PmlModel`: Absorbing (perfectly-matched-layer) material for far-field boundaries,
+//!   behind the `pml` feature
 //!
 //! Material models are used by particles to compute stress from deformation gradients.
 
+#[cfg(feature = "pml")]
+use crate::math::Vector;
 use bytemuck::{Pod, Zeroable};
 pub use drucker_prager::{DruckerPrager, DruckerPragerPlasticState};
 
@@ -42,6 +46,45 @@ pub struct ElasticCoefficients {
     pub cfl_coeff: f32,
 }
 
+/// Absorbing (perfectly-matched-layer) material, after Kurima, Chandra & Soga
+/// ([arXiv:2407.02790](https://arxiv.org/abs/2407.02790)).
+///
+/// Linear elasticity with per-axis stretched coordinates, which slows outgoing waves instead of
+/// reflecting them. The stretch does not dissipate: pair it with damping over the same region.
+#[cfg(feature = "pml")]
+#[derive(Copy, Clone, PartialEq, Debug, Pod, Zeroable)]
+#[repr(C)]
+pub struct PmlModel {
+    /// Elastic coefficients, which should match the material the layer is absorbing for.
+    pub elastic: ElasticCoefficients,
+    /// Per-axis coordinate stretching `C'_j`, zero outside the absorbing layer. Always three
+    /// entries so the GPU layout is identical in 2D and 3D; the third is unused in 2D.
+    pub stretch: [f32; 3],
+}
+
+/// Maximum stretch `α` at the outer edge of an absorbing layer. The paper's parameter study
+/// settles on 4, with little further gain beyond ~3.2.
+#[cfg(feature = "pml")]
+pub const DEFAULT_PML_MAX_STRETCH: f32 = 4.0;
+
+/// Computes the PML coordinate stretching for a point inside an absorbing layer.
+///
+/// The layer wraps the box `[interior_mins, interior_maxs]`, `thickness` deep on every side, and
+/// the stretch ramps linearly from zero at its inner boundary to `max_stretch` at the outer one.
+#[cfg(feature = "pml")]
+pub fn pml_stretch(
+    position: Vector,
+    interior_mins: Vector,
+    interior_maxs: Vector,
+    thickness: f32,
+    max_stretch: f32,
+) -> Vector {
+    let below = (interior_mins - position).max(Vector::ZERO);
+    let above = (position - interior_maxs).max(Vector::ZERO);
+    let depth = ((below + above) / thickness.max(1.0e-6)).min(Vector::ONE);
+    depth * max_stretch
+}
+
 impl ElasticCoefficients {
     /// Creates elastic coefficients from engineering parameters.
     ///
diff --git a/src/solver/boundary_condition.rs b/src/solver/boundary_condition.rs
index 2d18d2a..fa93982 100644
--- a/src/solver/boundary_condition.rs
+++ b/src/solver/boundary_condition.rs
@@ -1,13 +1,26 @@
+use crate::models::ElasticCoefficients;
 use bytemuck::{Pod, Zeroable};
 use slang_hal::BufferUsages;
 use slang_hal::backend::Backend;
 use stensor::tensor::GpuVector;
 
+/// Boundary condition applied to the grid nodes in contact with a collider.
+///
+/// The memory layout must match the shader-side `BoundaryCondition` struct in
+/// `shaders/slosh/solver/boundary_condition.slang`.
 #[derive(Copy, Clone, Debug, PartialEq, Pod, Zeroable)]
 #[repr(C)]
 pub struct GpuBoundaryCondition {
     pub ty: u32,
     pub friction: f32,
+    /// Pressure (dilatational) wave speed of the material in contact with the boundary (m/s).
+    ///
+    /// Only read by the [`Self::NON_REFLECTING`] boundary condition.
+    pub wave_speed_p: f32,
+    /// Shear wave speed of the material in contact with the boundary (m/s).
+    ///
+    /// Only read by the [`Self::NON_REFLECTING`] boundary condition.
+    pub wave_speed_s: f32,
 }
 
 impl GpuBoundaryCondition {
@@ -19,42 +32,96 @@ impl GpuBoundaryCondition {
     pub const DISABLED: u32 = 5u32;
 
     pub fn stick() -> GpuBoundaryCondition {
-        Self {
-            ty: Self::STICK,
-            friction: 0.0,
-        }
+        Self::new(Self::STICK, 0.0)
     }
 
     pub fn slip() -> GpuBoundaryCondition {
-        Self {
-            ty: Self::SLIP,
-            friction: 0.0,
-        }
+        Self::new(Self::SLIP, 0.0)
     }
 
     pub fn separate(friction: f32) -> GpuBoundaryCondition {
-        Self {
-            ty: Self::SEPARATE,
-            friction,
-        }
+        Self::new(Self::SEPARATE, friction)
     }
 
     pub fn friction_z_up(friction: f32) -> GpuBoundaryCondition {
-        Self {
-            ty: Self::FRICTION_Z_UP,
-            friction,
-        }
+        Self::new(Self::FRICTION_Z_UP, friction)
     }
 
     pub fn disabled() -> GpuBoundaryCondition {
+        Self::new(Self::DISABLED, 0.0)
+    }
+
+    /// An absorbing (non-reflecting) boundary based on Lysmer-Kuhlemeyer viscous dashpots.
+    ///
+    /// The boundary applies the traction a semi-infinite continuation of the material would,
+    /// `-ρ·c_p·v_n` along the normal and `-ρ·c_s·v_t` along the tangent, so a wave at normal
+    /// incidence is absorbed rather than reflected. Absorption degrades away from that incidence.
+    ///
+    /// The traction is graded over a band several cells deep (see `ABSORBING_LAYERS` in
+    /// `shaders/slosh/solver/boundary_condition.slang`), which the domain must have room for.
+    ///
+    /// It is layered on top of the [`Self::separate`] contact response, so `friction` behaves as
+    /// it does there. Note the shear dashpot damps tangential velocity all through the band, so
+    /// material resting inside it is dragged to a halt; pass a zero `wave_speed_s` to leave
+    /// sliding alone.
+    ///
+    /// # Arguments
+    ///
+    /// * `wave_speed_p` - Pressure wave speed `sqrt((λ + 2μ) / ρ)` of the material in contact (m/s)
+    /// * `wave_speed_s` - Shear wave speed `sqrt(μ / ρ)` of the material in contact (m/s), or
+    ///   zero to absorb the pressure wave only
+    /// * `friction` - Coulomb friction coefficient of the contact response, as in [`Self::separate`]
+    ///
+    /// See [`Self::non_reflecting_for_material`] for computing the wave speeds from engineering
+    /// parameters.
+    pub fn non_reflecting(
+        wave_speed_p: f32,
+        wave_speed_s: f32,
+        friction: f32,
+    ) -> GpuBoundaryCondition {
         Self {
-            ty: Self::DISABLED,
-            friction: 0.0,
+            ty: Self::NON_REFLECTING,
+            friction,
+            wave_speed_p,
+            wave_speed_s,
         }
     }
 
-    pub fn non_reflecting() -> GpuBoundaryCondition {
-        todo!();
+    /// Same as [`Self::non_reflecting`], but derives the wave speeds from the elastic
+    /// properties of the material in contact with the boundary.
+    ///
+    /// # Arguments
+    ///
+    /// * `young_modulus` - Young’s modulus E (Pa) of the material in contact
+    /// * `poisson_ratio` - Poisson’s ratio ν of the material in contact
+    /// * `density` - Density ρ of the material in contact (kg/m³, or kg/m² in 2D)
+    /// * `friction` - Coulomb friction coefficient of the contact response
+    pub fn non_reflecting_for_material(
+        young_modulus: f32,
+        poisson_ratio: f32,
+        density: f32,
+        friction: f32,
+    ) -> GpuBoundaryCondition {
+        let (p, s) = Self::wave_speeds(young_modulus, poisson_ratio, density);
+        Self::non_reflecting(p, s, friction)
+    }
+
+    /// Pressure and shear wave speeds `(c_p, c_s)` of an isotropic linear elastic material.
+    pub fn wave_speeds(young_modulus: f32, poisson_ratio: f32, density: f32) -> (f32, f32) {
+        let coeffs = ElasticCoefficients::from_young_modulus(young_modulus, poisson_ratio);
+        (
+            ((coeffs.lambda + 2.0 * coeffs.mu) / density).sqrt(),
+            (coeffs.mu / density).sqrt(),
+        )
+    }
+
+    fn new(ty: u32, friction: f32) -> GpuBoundaryCondition {
+        Self {
+            ty,
+            friction,
+            wave_speed_p: 0.0,
+            wave_speed_s: 0.0,
+        }
     }
 }
 
diff --git a/src/solver/particle.rs b/src/solver/particle.rs
index 952be49..e088323 100644
--- a/src/solver/particle.rs
+++ b/src/solver/particle.rs
@@ -45,7 +45,19 @@ pub struct ParticleDynamics {
     ///
     /// Applies a damping force proportional to velocity: F_damp = -damping * m * v.
     /// Typical values: 0.0 (no damping) to 10.0 (heavy damping).
+    ///
+    /// Opposes absolute velocity, so it damps rigid-body motion too: a body in free flight is
+    /// slowed like everything else. Use `stiffness_damping` where that matters.
     pub damping: f32,
+    /// Rayleigh stiffness-proportional damping coefficient (s).
+    ///
+    /// Adds a viscous stress `a_K * C : sym(grad v)`. Unlike [`Self::damping`] it is blind to
+    /// rigid-body motion, so it attenuates waves without dragging on bulk movement.
+    ///
+    /// It tightens the explicit stability bound, which [`crate::solver::WgTimestepBounds`]
+    /// accounts for; keeping it below the timestep leaves that bound untouched.
+    #[cfg(feature = "pml")]
+    pub stiffness_damping: f32,
     /// The particle phase (used by materials that can break).
     pub phase: f32,
     /// Whether this particle is active (1) or disabled (0).
@@ -76,6 +88,8 @@ impl ParticleDynamics {
             init_radius: radius,
             mass: init_volume * density,
             damping: 0.0,
+            #[cfg(feature = "pml")]
+            stiffness_damping: 0.0,
             cdf: Cdf::default(),
             phase: 1.0,
             enabled: 1,
@@ -88,11 +102,17 @@ impl ParticleDynamics {
         self.fixed = fixed as u32;
     }
 
-    /// Sets the damping coefficient for this particle.
+    /// Sets the mass-proportional damping coefficient for this particle.
     pub fn set_damping(&mut self, damping: f32) {
         self.damping = damping;
     }
 
+    /// Sets the stiffness-proportional damping coefficient for this particle.
+    #[cfg(feature = "pml")]
+    pub fn set_stiffness_damping(&mut self, stiffness_damping: f32) {
+        self.stiffness_damping = stiffness_damping;
+    }
+
     /// Updates the particle mass based on a new density.
     ///
     /// Keeps the initial volume constant and recomputes mass = volume × density.
@@ -106,6 +126,8 @@ impl ParticleDynamics {
             affine: self.affine,
             velocity: self.velocity,
             force_dt: self.force_dt,
+            #[cfg(feature = "pml")]
+            mass_scale: Vector::ONE,
             vel_grad_det: self.vel_grad_det,
             mass: self.mass,
             enabled: self.enabled,
@@ -118,6 +140,8 @@ impl ParticleDynamics {
             init_volume: self.init_volume,
             init_radius: self.init_radius,
             damping: self.damping,
+            #[cfg(feature = "pml")]
+            stiffness_damping: self.stiffness_damping,
             phase: self.phase,
             fixed: self.fixed,
         }
@@ -138,6 +162,10 @@ pub struct Kinematics {
     pub velocity: Vector,
     /// Additional force * dt applied to the particle.
     pub force_dt: Vector,
+    /// Per-axis multiplier on this particle's inertia, one for ordinary materials (the PML uses
+    /// `s_j²`). Recomputed from the material model at every particle update.
+    #[cfg(feature = "pml")]
+    pub mass_scale: Vector,
     /// Determinant of velocity gradient (for volume change tracking).
     pub vel_grad_det: f32,
     /// Particle mass (kg).
@@ -160,6 +188,9 @@ pub struct ParticleProperties {
     pub init_radius: f32,
     /// Rayleigh mass-proportional damping coefficient (1/s).
     pub damping: f32,
+    /// Rayleigh stiffness-proportional damping coefficient (s).
+    #[cfg(feature = "pml")]
+    pub stiffness_damping: f32,
     /// The particle phase (used by materials that can break).
     pub phase: f32,
     /// Whether this particle is fixed (1) or dynamic (0).
diff --git a/src/solver/particle_model.rs b/src/solver/particle_model.rs
index 859e3d5..379e356 100644
--- a/src/solver/particle_model.rs
+++ b/src/solver/particle_model.rs
@@ -1,3 +1,7 @@
+#[cfg(feature = "pml")]
+use crate::math::Vector;
+#[cfg(feature = "pml")]
+use crate::models::PmlModel;
 use crate::models::{DruckerPrager, DruckerPragerPlasticState, ElasticCoefficients};
 use bytemuck::{NoUninit, Pod, Zeroable};
 
@@ -15,6 +19,9 @@ pub enum ParticleModel {
     SandLinear(SandModel),
     /// Sand with Neo-Hookean elasticity and Drucker-Prager plasticity.
     SandNeoHookean(SandModel),
+    /// Absorbing (perfectly-matched-layer) material for far-field boundaries.
+    #[cfg(feature = "pml")]
+    AbsorbingPml(PmlModel),
 }
 
 impl Default for ParticleModel {
@@ -77,6 +84,46 @@ impl ParticleModel {
             elastic: ElasticCoefficients::from_young_modulus(young_modulus, poisson_ratio),
         })
     }
+
+    /// Creates an absorbing (perfectly-matched-layer) material for a far-field boundary.
+    ///
+    /// Particles carrying this model form a layer around the region of interest that lets outgoing
+    /// waves leave instead of reflecting them. The elastic parameters should match the surrounded
+    /// material; [`crate::models::pml_stretch`] computes `stretch` from the layer's geometry.
+    ///
+    /// The stretch redirects energy but does not dissipate it, so pair it with
+    /// [`crate::solver::ParticleDynamics::damping`] over the same particles (the paper uses
+    /// `α_M = 1`; what matters is `α_M·L/c`, the attenuation per transit of the layer). Prefer
+    /// [`crate::solver::ParticleDynamics::stiffness_damping`] when the domain itself moves, since
+    /// mass-proportional damping would resist that motion.
+    ///
+    /// # Gravity
+    ///
+    /// The layer's inertia is `s_j²` times the real mass while gravity acts on the real one, so it
+    /// falls `s_j²` times slower along a stretched axis instead of carrying that much weight.
+    ///
+    /// The stretched operator's static stiffness still differs from the elastic one, so a layer
+    /// under sustained load settles more than the material would (about twice, for a 4 m column on
+    /// a 2 m layer). The paper removes this with a geo-static pre-step, not implemented here;
+    /// keeping absorbing layers out of the load path avoids it.
+    ///
+    /// # Arguments
+    ///
+    /// * `young_modulus` - Young's modulus E (Pa) of the surrounded material
+    /// * `poisson_ratio` - Poisson's ratio ν of the surrounded material
+    /// * `stretch` - Per-axis coordinate stretching `C'_j` (zero means plain linear elasticity)
+    #[cfg(feature = "pml")]
+    pub fn absorbing_pml(young_modulus: f32, poisson_ratio: f32, stretch: Vector) -> Self {
+        #[cfg(feature = "dim2")]
+        let stretch = [stretch.x, stretch.y, 0.0];
+        #[cfg(feature = "dim3")]
+        let stretch = [stretch.x, stretch.y, stretch.z];
+
+        ParticleModel::AbsorbingPml(PmlModel {
+            elastic: ElasticCoefficients::from_young_modulus(young_modulus, poisson_ratio),
+            stretch,
+        })
+    }
 }
 
 /// GPU-compatible version of [`ParticleModel`] with explicit padding.
@@ -95,6 +142,9 @@ pub enum GpuParticleModel {
     SandLinear(SandModel) = 2,
     /// Sand with Neo-Hookean elasticity and Drucker-Prager plasticity.
     SandNeoHookean(SandModel) = 3,
+    /// Absorbing (perfectly-matched-layer) material with padding for GPU alignment.
+    #[cfg(feature = "pml")]
+    AbsorbingPml(PmlModel, [u32; 6]) = 4,
 }
 
 // IMPORTANT: this assertions is here to reduce risks of `GpuParticleModel` from mismatching
@@ -114,6 +164,8 @@ impl From<ParticleModel> for GpuParticleModel {
             ParticleModel::SandNeoHookean(sand_neo_hookean) => {
                 GpuParticleModel::SandNeoHookean(sand_neo_hookean)
             }
+            #[cfg(feature = "pml")]
+            ParticleModel::AbsorbingPml(pml) => GpuParticleModel::AbsorbingPml(pml, [0; _]),
         }
     }
 }
@@ -131,6 +183,8 @@ impl From<GpuParticleModel> for ParticleModel {
             GpuParticleModel::SandNeoHookean(sand_neo_hookean) => {
                 ParticleModel::SandNeoHookean(sand_neo_hookean)
             }
+            #[cfg(feature = "pml")]
+            GpuParticleModel::AbsorbingPml(pml, _) => ParticleModel::AbsorbingPml(pml),
         }
     }
 }
diff --git a/src_testbed/data.rs b/src_testbed/data.rs
index b0332fe..f54e66c 100644
--- a/src_testbed/data.rs
+++ b/src_testbed/data.rs
@@ -20,6 +20,9 @@ pub struct AppState<GpuModel: GpuParticleModelData = GpuParticleModel> {
     pub num_substeps: u32,
     pub gravity_factor: f32,
     pub restarting: bool,
+    /// Set from a scene's extra UI to rebuild the current scene at the end of the frame, for
+    /// settings that change what the scene is made of. Same path as the "Restart" button.
+    pub restart_requested: bool,
     // pub hot_reload: HotReloadState,
     pub show_rigid_particles: bool,
     pub cell_width: f32,
diff --git a/src_testbed/lib.rs b/src_testbed/lib.rs
index 619404d..caf603a 100644
--- a/src_testbed/lib.rs
+++ b/src_testbed/lib.rs
@@ -43,8 +43,10 @@ use std::rc::Rc;
 use stensor::tensor::GpuTensor;
 use wgpu::Limits;
 
-type SceneBuilders<GpuModel> = Vec<(String, SceneBuildFn<GpuModel>)>;
-type SceneBuildFn<GpuModel> = fn(&WebGpu, &mut AppState<GpuModel>) -> PhysicsContext<GpuModel>;
+/// The scene list a testbed is run with: a display name and a builder per scene.
+pub type SceneBuilders<GpuModel> = Vec<(String, SceneBuildFn<GpuModel>)>;
+/// Builds one scene, called on startup and on every restart.
+pub type SceneBuildFn<GpuModel> = fn(&WebGpu, &mut AppState<GpuModel>) -> PhysicsContext<GpuModel>;
 
 /// GPU-construction options for the testbed.
 pub struct TestbedConfig {
@@ -138,6 +140,7 @@ impl<GpuModel: GpuParticleModelData> Stage<GpuModel> {
             num_substeps: 1,
             gravity_factor: 1.0,
             restarting: false,
+            restart_requested: false,
             show_rigid_particles: false,
             cell_width: 0.01,
             particle_colors: None,
@@ -398,7 +401,7 @@ pub async fn run_with_hooks_and_config<GpuModel: GpuParticleModelData>(
         hooks,
         config,
         scene_builders,
-        |_, _, _, _| None,
+        |_, _, _, _, _| None,
         #[cfg(feature = "dim3")]
         up_axis,
     )
@@ -412,6 +415,7 @@ pub async fn run_with_hooks_and_ui<GpuModel: GpuParticleModelData>(
     scene_builders: SceneBuilders<GpuModel>,
     mut extra_ui: impl FnMut(
         &egui::Context,
+        &mut AppState<GpuModel>,
         &PhysicsContext<GpuModel>,
         &SimulationStepResult,
         bool,
@@ -623,13 +627,23 @@ pub async fn run_with_hooks_and_ui<GpuModel: GpuParticleModelData>(
                 }
             });
 
-            ui_run_state = extra_ui(ctx, &stage.physics, &stage.step_result, stepped);
+            ui_run_state = extra_ui(
+                ctx,
+                &mut stage.app_state,
+                &stage.physics,
+                &stage.step_result,
+                stepped,
+            );
         });
 
         if let Some(run_state) = ui_run_state {
             stage.app_state.run_state = run_state;
         }
 
+        if std::mem::take(&mut stage.app_state.restart_requested) {
+            new_selected_demo = Some(stage.selected_demo);
+        }
+
         if let Some(demo) = new_selected_demo {
             stage.set_demo(demo);
             render_colliders(&mut window, &stage.physics, &mut colliders_gfx);