diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs index 05599ed..cc7017e 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/exchange.rs @@ -140,4 +140,74 @@ impl Mask { } Some(force) } + + /// The closure lag of a moving body's pressure correction: the + /// corrector applies `p'` on the STEP apertures while the operator + /// route reads the summed pressure on the END apertures, so the exact + /// discrete force carries `Σ_c p'_c (W_step,c − W_end,c)` (a force on + /// the body) that the route lacks. Zero for a body at rest. + #[must_use] + pub fn closure_lag(&self, p_prime: &[f64]) -> [f64; 3] { + let g = self.grid; + let area = [g.dy * g.dz, g.dx * g.dz, g.dx * g.dy]; + let mut lag = [0.0; 3]; + if self.step_apertures.is_none() { + return lag; + } + for k in 0..g.nz { + for j in 0..g.ny { + for i in 0..g.nx { + let idx = g.cell(k, j, i); + if !self.cell_active(idx) { + continue; + } + let pp = p_prime[idx]; + if pp == 0.0 { + continue; + } + // W_c = −Σ A_f n_f: x-part −(α_e − α_w) A_x, etc. + let (ue, uw) = (g.uface(k, j, i + 1), g.uface(k, j, i)); + let (vn, vs) = (g.vface(k, j + 1, i), g.vface(k, j, i)); + let (wt, wb) = (g.wface(k + 1, j, i), g.wface(k, j, i)); + let d = [ + (self.au_step(ue) - self.a_u(ue)) - (self.au_step(uw) - self.a_u(uw)), + (self.av_step(vn) - self.a_v(vn)) - (self.av_step(vs) - self.a_v(vs)), + (self.aw_step(wt) - self.a_w(wt)) - (self.aw_step(wb) - self.a_w(wb)), + ]; + for c in 0..3 { + lag[c] += pp * (-d[c]) * area[c]; + } + } + } + } + lag + } + + /// The force the discrete momentum equation actually applied over the + /// last step, read post-step: the pressure part on the end field, the + /// implicit wall shear on the PREDICTED velocities `u*` (the corrector + /// moves `u` without re-applying the shear) and the explicit wall + /// exchange on the OLD velocities. On a body at rest at a steady state + /// this equals `cut_wall_force`; on a moving body it is the number the + /// box route should reproduce. + pub fn cut_wall_force_applied( + &self, + body: &Body, + f: &Field, + mu: f64, + rho: f64, + t: f64, + ) -> Option<[f64; 3]> { + let mut star = f.clone(); + star.u.copy_from_slice(&f.u_star); + star.v.copy_from_slice(&f.v_star); + star.w.copy_from_slice(&f.w_star); + let (p, s) = self.cut_wall_force_parts(body, &star, mu, t)?; + let mut old = f.clone(); + old.u.copy_from_slice(&f.u_old); + old.v.copy_from_slice(&f.v_old); + old.w.copy_from_slice(&f.w_old); + let x = self.cut_wall_exchange_force(body, &old, mu, rho, t, None)?; + Some([p[0] + s[0] + x[0], p[1] + s[1] + x[1], p[2] + s[2] + x[2]]) + } } diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs index 72b186f..11bdac4 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/loads.rs @@ -288,6 +288,25 @@ impl Mask { /// without it (slip or periodic sides) the z faces carry nothing. #[allow(clippy::too_many_arguments)] pub fn control_volume_force_with_walls( + &self, + f: &Field, + dt: f64, + rho: f64, + mu: f64, + source: Option<&dyn Fn(f64, f64, f64) -> (f64, f64, f64)>, + bx: (usize, usize, usize, usize, usize, usize), + no_slip_z: bool, + ) -> [f64; 3] { + self.control_volume_force_moving(f, dt, rho, mu, source, bx, no_slip_z, None) + } + + /// The box route on a moving body: the unsteady term is the change of + /// the momentum `ρ V u` over the step with the previous step's cell + /// volumes `vol_old` (`Solver::previous_volumes`), so the fluid the + /// wall sweeps counts (`Σ ρ V (u − uⁿ)/dt` misses `ρ (V − Vⁿ) uⁿ/dt`). + /// `None` falls back to the fixed-volume form. + #[allow(clippy::too_many_arguments)] + pub fn control_volume_force_moving( &self, f: &Field, dt: f64, @@ -296,6 +315,7 @@ impl Mask { source: Option<&dyn Fn(f64, f64, f64) -> (f64, f64, f64)>, (i0, i1, j0, j1, k0, k1): (usize, usize, usize, usize, usize, usize), no_slip_z: bool, + vol_old: Option<&[f64]>, ) -> [f64; 3] { let g = self.grid(); let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz); @@ -427,19 +447,28 @@ impl Mask { for j in j0..j1 { for i in i0..i1 { let idx = g.cell(k, j, i); - if !self.is_fluid_cell(idx) { + let v_new = if self.is_fluid_cell(idx) { + self.vol(idx) + } else { + 0.0 + }; + let v_old = vol_old.map_or(v_new, |vo| vo[idx]); + if v_new == 0.0 && v_old == 0.0 { continue; } - let dv = dv * self.vol(idx); let (fu0, fu1) = (g.uface(k, j, i), g.uface(k, j, i + 1)); let (fv0, fv1) = (g.vface(k, j, i), g.vface(k, j + 1, i)); let (fw0, fw1) = (g.wface(k, j, i), g.wface(k + 1, j, i)); - let du = 0.5 * ((u[fu0] - f.u_old[fu0]) + (u[fu1] - f.u_old[fu1])); - let dvv = 0.5 * ((v[fv0] - f.v_old[fv0]) + (v[fv1] - f.v_old[fv1])); - let dw = 0.5 * ((w[fw0] - f.w_old[fw0]) + (w[fw1] - f.w_old[fw1])); - force[0] -= rho * du / dt * dv; - force[1] -= rho * dvv / dt * dv; - force[2] -= rho * dw / dt * dv; + let cen = |a: &[f64], b: &[f64], f0: usize, f1: usize| { + (0.5 * (a[f0] + a[f1]), 0.5 * (b[f0] + b[f1])) + }; + let (un, uo) = cen(u, &f.u_old, fu0, fu1); + let (vn, vo) = cen(v, &f.v_old, fv0, fv1); + let (wn, wo) = cen(w, &f.w_old, fw0, fw1); + force[0] -= rho * (v_new * un - v_old * uo) / dt * dv; + force[1] -= rho * (v_new * vn - v_old * vo) / dt * dv; + force[2] -= rho * (v_new * wn - v_old * wo) / dt * dv; + let dv = dv * v_new; if let Some(s) = source { let (sx, sy, sz) = s( (i as f64 + 0.5) * dx, diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs index 7798d10..f984203 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/cut_predictor.rs @@ -25,6 +25,7 @@ impl Solver { let periodic = self.params.boundaries.periodic_z(); let lat = mask.lattice(); let w_range = if periodic { 0..nz } else { 1..nz }; + let mut floor = [0.0; 3]; for c in 0..3 { let (ir, jr, kr) = match c { 0 => (1..nx, 0..ny, 0..nz), @@ -45,7 +46,9 @@ impl Solver { } let p = [i as i64, j as i64, k as i64]; let idx = lat.face(c, p).expect("face"); - updates.push((idx, self.cut_face_update(field, c, p, dt, t_old))); + let (val, src) = self.cut_face_update(field, c, p, dt, t_old); + floor[c] += src; + updates.push((idx, val)); } } } @@ -58,6 +61,7 @@ impl Solver { out[idx] = val; } } + self.floor_source.set(floor); if periodic { for j in 0..ny { for i in 0..nx { @@ -68,9 +72,17 @@ impl Solver { } /// The predicted value of the unknown face of component `c` at lattice - /// `p` from the old field. + /// `p` from the old field, with the momentum the inertia floor added + /// on this face, `ρ (V_eff − V_α)(u* − uⁿ)/dt`. #[allow(clippy::too_many_lines)] - fn cut_face_update(&self, field: &Field, c: usize, p: [i64; 3], dt: f64, t_old: f64) -> f64 { + fn cut_face_update( + &self, + field: &Field, + c: usize, + p: [i64; 3], + dt: f64, + t_old: f64, + ) -> (f64, f64) { let mask = self.mask.as_ref().expect("cut mask"); let body = self.body.as_ref().expect("body"); let g = field.grid; @@ -223,6 +235,9 @@ impl Solver { }; let v_eff = fraction.max(INERTIA_FLOOR) * h[c] * area[c]; let inertia = rho * v_eff / dt; - (inertia * u0 - conv + diff + pressure + source + shear * ub) / (inertia + shear) + let u_star = + (inertia * u0 - conv + diff + pressure + source + shear * ub) / (inertia + shear); + let v_alpha = fraction * h[c] * area[c]; + (u_star, rho * (v_eff - v_alpha) * (u_star - u0) / dt) } } diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs index 3aea502..e8299ce 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/mod.rs @@ -7,6 +7,7 @@ mod cut_predictor; #[cfg(feature = "cuda")] pub mod device; +mod moving; mod predictor; mod projection; @@ -128,6 +129,17 @@ pub struct Solver { pub fluid: Fluid, pub params: Parameters, pub(super) momentum_source: Option, + /// The momentum the small-cell inertia floor added over the last + /// step, `Σ ρ (V_eff − V_α)(u* − uⁿ)/dt` over the unknown faces (host + /// cut predictor only): the box route counts it, the operator route + /// does not. + pub(super) floor_source: std::cell::Cell<[f64; 3]>, + /// The closure lag of the pressure corrections over the last step + /// (`Mask::closure_lag` summed over the correctors; host path). + pub(super) pressure_lag: std::cell::Cell<[f64; 3]>, + /// The previous mask's cell volumes (moving bodies; the box route's + /// unsteady term). + pub(super) vol_old: Vec, boundary_velocity: Option, body: Option, /// The body moves: the mask is rebuilt at every step's new time. @@ -158,6 +170,9 @@ impl Solver { fluid, params, momentum_source: None, + floor_source: std::cell::Cell::new([0.0; 3]), + pressure_lag: std::cell::Cell::new([0.0; 3]), + vol_old: Vec::new(), boundary_velocity: None, body: None, moving: false, @@ -221,6 +236,27 @@ impl Solver { .expect("embedded mask") } + /// The small-cell floor's momentum source over the last step (a force + /// on the fluid; host cut predictor only, zero otherwise). + #[must_use] + pub fn floor_source(&self) -> [f64; 3] { + self.floor_source.get() + } + + /// The pressure corrections' closure lag over the last step (a force + /// on the body the operator route lacks on a moving body; host path). + #[must_use] + pub fn pressure_lag(&self) -> [f64; 3] { + self.pressure_lag.get() + } + + /// The previous step's cell volumes (fraction of the cell; 0 for a + /// solid cell), once a moving body's mask has been rebuilt. + #[must_use] + pub fn previous_volumes(&self) -> Option<&[f64]> { + (!self.vol_old.is_empty()).then_some(self.vol_old.as_slice()) + } + #[must_use] pub fn body(&self) -> Option<&Body> { self.body.as_ref() @@ -506,69 +542,6 @@ impl Solver { &self.wall_fluxes } - /// The moving body's mask at the end-of-step geometry `t_new`: the - /// pressure of the cells that just became fluid refilled from their - /// neighbours (fluid in both masks), the new mask's prescribed and - /// ghost values imposed from the previous corrected field, the - /// step-averaged apertures and the GCL wall-flux table (cut wall). - /// Returns the fresh-cell count. `field` holds the predicted field. - pub fn rebuild_moving_mask(&mut self, field: &mut Field, dt: f64, t_new: f64) -> usize { - let Some(body) = &self.body else { - return 0; - }; - let mut fresh_cells = 0; - let mut new_mask = self.build_mask(body, field.grid, t_new, dt); - if let Some(old_mask) = &self.mask { - fresh_cells = refill_fresh_cells(old_mask, &new_mask, field); - let n_in = self.params.aperture_substeps; - if n_in == 0 { - new_mask.set_step_apertures(old_mask); - } else { - // The intermediate geometries, each within the band of the - // previous one. - let g = field.grid; - let h = g.dx.min(g.dy).min(g.dz); - let t_old = t_new - dt; - let mut inner: Vec = Vec::with_capacity(n_in); - for m in 1..=n_in { - let tm = t_old + dt * m as f64 / (n_in + 1) as f64; - let prev: Option<(&CutGeometry, f64, f64)> = match ( - self.params.max_surface_speed, - inner.last().or(old_mask.cut()), - ) { - (Some(speed), Some(c)) => { - Some((c, 3.0 * h, speed * dt / (n_in + 1) as f64)) - } - _ => None, - }; - inner.push(CutGeometry::build_from(body, g, tm, prev)); - } - let refs: Vec<&CutGeometry> = inner.iter().collect(); - new_mask.set_step_apertures_with(old_mask, &refs); - } - } - new_mask.impose_from( - body, - &field.u_old, - &field.v_old, - &field.w_old, - &mut field.u, - &mut field.v, - &mut field.w, - t_new, - ); - if new_mask.cut().is_some() { - let (table, correction) = match &self.mask { - Some(old_mask) => new_mask.gcl_flux_table(old_mask, dt), - None => new_mask.wall_flux_table(body, t_new), - }; - self.wall_fluxes = table; - self.last_ghost_correction = correction; - } - self.mask = Some(new_mask); - fresh_cells - } - /// One step of `dt`: predictor, correctors, clock. pub fn advance(&mut self, field: &mut Field, dt: f64) -> StepResult { assert!( @@ -581,6 +554,7 @@ impl Solver { let t_old = self.time; let t_new = t_old + dt; field.update_old_values(); + self.pressure_lag.set([0.0; 3]); self.momentum_predictor(field, dt, t_old); self.apply_boundary_normals(field, t_new); // A moving body: the mask at the end-of-step geometry. @@ -607,6 +581,14 @@ impl Solver { let mut poisson_iterations = 0; for corrector in 0..self.params.corrector_steps.max(1) { let sol = self.solve_correction(field, dt, corrector == 0); + if let Some(m) = self.mask.as_ref() { + let l = m.closure_lag(&field.p_prime); + let mut acc = self.pressure_lag.get(); + for c in 0..3 { + acc[c] += l[c]; + } + self.pressure_lag.set(acc); + } poisson_iterations += sol.iterations; let mass_residual = self.apply_correction(field, dt); if std::env::var_os("RTX_E3_DEBUG").is_some() { @@ -637,61 +619,3 @@ impl Solver { } } } - -/// Refill the pressure of the cells fluid in `new` and not in `old` from -/// their face neighbours fluid in both; returns their count. -fn refill_fresh_cells(old: &Mask, new: &Mask, field: &mut Field) -> usize { - let g = field.grid; - let (nx, ny, nz) = (g.nx, g.ny, g.nz); - let periodic = new.periodic_z(); - let mut fresh = 0; - let mut refills = Vec::new(); - for k in 0..nz { - for j in 0..ny { - for i in 0..nx { - let idx = g.cell(k, j, i); - if !(new.is_fluid_cell(idx) && !old.is_fluid_cell(idx)) { - continue; - } - fresh += 1; - let mut sum = 0.0; - let mut count = 0usize; - let mut visit = |nb: usize| { - if new.is_fluid_cell(nb) && old.is_fluid_cell(nb) { - sum += field.p[nb]; - count += 1; - } - }; - if i + 1 < nx { - visit(g.cell(k, j, i + 1)); - } - if i > 0 { - visit(g.cell(k, j, i - 1)); - } - if j + 1 < ny { - visit(g.cell(k, j + 1, i)); - } - if j > 0 { - visit(g.cell(k, j - 1, i)); - } - if k + 1 < nz { - visit(g.cell(k + 1, j, i)); - } else if periodic && nz > 1 { - visit(g.cell(0, j, i)); - } - if k > 0 { - visit(g.cell(k - 1, j, i)); - } else if periodic && nz > 1 { - visit(g.cell(nz - 1, j, i)); - } - if count > 0 { - refills.push((idx, sum / count as f64)); - } - } - } - } - for (idx, p) in refills { - field.p[idx] = p; - } - fresh -} diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/moving.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/moving.rs new file mode 100644 index 0000000..5fe38ad --- /dev/null +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/moving.rs @@ -0,0 +1,141 @@ +//! The moving body's per-step mask rebuild: fresh-cell pressure refill, +//! imposition from history, step-averaged apertures, the GCL wall-flux +//! table and the previous volumes for the box route. +use super::super::cut::CutGeometry; +use super::super::field::Field; +use super::super::wall::{Mask, WallScheme}; +use super::Solver; + +impl Solver { + /// The moving body's mask at the end-of-step geometry `t_new`: the + /// pressure of the cells that just became fluid refilled from their + /// neighbours (fluid in both masks), the new mask's prescribed and + /// ghost values imposed from the previous corrected field, the + /// step-averaged apertures and the GCL wall-flux table (cut wall). + /// Returns the fresh-cell count. `field` holds the predicted field. + pub fn rebuild_moving_mask(&mut self, field: &mut Field, dt: f64, t_new: f64) -> usize { + let Some(body) = &self.body else { + return 0; + }; + let mut fresh_cells = 0; + let mut new_mask = self.build_mask(body, field.grid, t_new, dt); + if let Some(old_mask) = &self.mask { + fresh_cells = refill_fresh_cells(old_mask, &new_mask, field); + let n_in = self.params.aperture_substeps; + if n_in == 0 { + new_mask.set_step_apertures(old_mask); + } else { + // The intermediate geometries, each within the band of the + // previous one. + let g = field.grid; + let h = g.dx.min(g.dy).min(g.dz); + let t_old = t_new - dt; + let mut inner: Vec = Vec::with_capacity(n_in); + for m in 1..=n_in { + let tm = t_old + dt * m as f64 / (n_in + 1) as f64; + let prev: Option<(&CutGeometry, f64, f64)> = match ( + self.params.max_surface_speed, + inner.last().or(old_mask.cut()), + ) { + (Some(speed), Some(c)) => { + Some((c, 3.0 * h, speed * dt / (n_in + 1) as f64)) + } + _ => None, + }; + inner.push(CutGeometry::build_from(body, g, tm, prev)); + } + let refs: Vec<&CutGeometry> = inner.iter().collect(); + new_mask.set_step_apertures_with(old_mask, &refs); + } + } + new_mask.impose_from( + body, + &field.u_old, + &field.v_old, + &field.w_old, + &mut field.u, + &mut field.v, + &mut field.w, + t_new, + ); + if new_mask.cut().is_some() { + let (table, correction) = match &self.mask { + Some(old_mask) => new_mask.gcl_flux_table(old_mask, dt), + None => new_mask.wall_flux_table(body, t_new), + }; + self.wall_fluxes = table; + self.last_ghost_correction = correction; + } + if let Some(old) = &self.mask { + self.vol_old = (0..field.grid.cells()) + .map(|i| { + if old.is_fluid_cell(i) { + old.vol(i) + } else { + 0.0 + } + }) + .collect(); + } + self.mask = Some(new_mask); + fresh_cells + } +} + +/// Refill the pressure of the cells fluid in `new` and not in `old` from +/// their face neighbours fluid in both; returns their count. +fn refill_fresh_cells(old: &Mask, new: &Mask, field: &mut Field) -> usize { + let g = field.grid; + let (nx, ny, nz) = (g.nx, g.ny, g.nz); + let periodic = new.periodic_z(); + let mut fresh = 0; + let mut refills = Vec::new(); + for k in 0..nz { + for j in 0..ny { + for i in 0..nx { + let idx = g.cell(k, j, i); + if !(new.is_fluid_cell(idx) && !old.is_fluid_cell(idx)) { + continue; + } + fresh += 1; + let mut sum = 0.0; + let mut count = 0usize; + let mut visit = |nb: usize| { + if new.is_fluid_cell(nb) && old.is_fluid_cell(nb) { + sum += field.p[nb]; + count += 1; + } + }; + if i + 1 < nx { + visit(g.cell(k, j, i + 1)); + } + if i > 0 { + visit(g.cell(k, j, i - 1)); + } + if j + 1 < ny { + visit(g.cell(k, j + 1, i)); + } + if j > 0 { + visit(g.cell(k, j - 1, i)); + } + if k + 1 < nz { + visit(g.cell(k + 1, j, i)); + } else if periodic && nz > 1 { + visit(g.cell(0, j, i)); + } + if k > 0 { + visit(g.cell(k - 1, j, i)); + } else if periodic && nz > 1 { + visit(g.cell(nz - 1, j, i)); + } + if count > 0 { + refills.push((idx, sum / count as f64)); + } + } + } + } + for (idx, p) in refills { + field.p[idx] = p; + } + fresh +} diff --git a/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs b/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs index 8b50a72..5404855 100644 --- a/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs +++ b/crates/specialized/rtx-cfd/tests/embedded3_falsifier.rs @@ -204,6 +204,12 @@ struct Record { /// Load per unit span by the control-volume route (a box of whole /// cells around the body, reading no near-wall value). fy_cv: f64, + /// The inertia floor's momentum source per unit span (cut wall). + fy_floor: f64, + /// The pressure corrections' closure lag per unit span (cut wall). + fy_lag: f64, + /// The force the momentum equation applied (shear on u*, exchange on uⁿ). + fy_applied: f64, fresh: usize, skipped: usize, p_far: f64, @@ -294,14 +300,25 @@ fn run(scheme: WallScheme, moving: bool, dt: f64, t_end: f64) -> Run { } } fy /= lz; + let fy_floor = solver.floor_source()[1] / lz; + let fy_lag = solver.pressure_lag()[1] / lz; + let fy_applied = if scheme == WallScheme::CutCell { + mask.cut_wall_force_applied(body, &field, MU, RHO, t) + .expect("applied")[1] + / lz + } else { + 0.0 + }; let margin = 8; - let fy_cv = mask.control_volume_force( + let fy_cv = mask.control_volume_force_moving( &field, dt, RHO, MU, None, (margin, N - margin, margin, N - margin, 0, nz), + false, + solver.previous_volumes(), )[1] / lz; let p_far = field.p[g.cell(kp, jp, ip)]; let mut ke = 0.0; @@ -337,12 +354,36 @@ fn run(scheme: WallScheme, moving: bool, dt: f64, t_end: f64) -> Run { t, fy, fy_cv, + fy_floor, + fy_lag, + fy_applied, fresh: result.fresh_cells, skipped, p_far, ke, }); } + if scheme == WallScheme::CutCell { + let n = records.len().max(1) as f64; + let rms = |f: &dyn Fn(&Record) -> f64| { + (records.iter().map(|r| f(r) * f(r)).sum::() / n).sqrt() + }; + let mean = |f: &dyn Fn(&Record) -> f64| records.iter().map(f).sum::() / n; + println!( + " route residual (wall − box) RMS {:.3e} → with the floor source {:.3e}, with the closure lag {:.3e}, both {:.3e}, APPLIED form {:.3e} (mean {:+.3e}); means wall {:+.3e} box {:+.3e} floor {:+.3e} lag {:+.3e}; box RMS {:.3e}", + rms(&|r| r.fy - r.fy_cv), + rms(&|r| r.fy + r.fy_floor - r.fy_cv), + rms(&|r| r.fy + r.fy_lag - r.fy_cv), + rms(&|r| r.fy + r.fy_floor + r.fy_lag - r.fy_cv), + rms(&|r| r.fy_applied - r.fy_cv), + mean(&|r| r.fy_applied), + mean(&|r| r.fy), + mean(&|r| r.fy_cv), + mean(&|r| r.fy_floor), + mean(&|r| r.fy_lag), + rms(&|r| r.fy_cv) + ); + } Run { records, energy_per_flip,