///|
/// Local pressure under the isothermal D2Q9 equation of state.
pub fn Simulation::pressure_at(self : Simulation, x~ : Int, y~ : Int) -> Double {
  self.cell(x~, y~).rho / 3.0
}

///|
/// Centered divergence of a velocity field.
pub fn Simulation::divergence_at(
  self : Simulation,
  x~ : Int,
  y~ : Int,
) -> Double {
  let left = self.cell(x=(x - 1).clamp(min=0, max=self.size.width - 1), y~).ux
  let right = self.cell(x=(x + 1).clamp(min=0, max=self.size.width - 1), y~).ux
  let bottom = self.cell(x~, y=(y - 1).clamp(min=0, max=self.size.height - 1)).uy
  let top = self.cell(x~, y=(y + 1).clamp(min=0, max=self.size.height - 1)).uy
  (right - left + top - bottom) * 0.5
}

///|
/// Centered strain-rate tensor magnitude.
pub fn Simulation::strain_rate_at(
  self : Simulation,
  x~ : Int,
  y~ : Int,
) -> Double {
  let left = self.cell(x=(x - 1).clamp(min=0, max=self.size.width - 1), y~)
  let right = self.cell(x=(x + 1).clamp(min=0, max=self.size.width - 1), y~)
  let bottom = self.cell(x~, y=(y - 1).clamp(min=0, max=self.size.height - 1))
  let top = self.cell(x~, y=(y + 1).clamp(min=0, max=self.size.height - 1))
  let dux_dx = (right.ux - left.ux) * 0.5
  let duy_dy = (top.uy - bottom.uy) * 0.5
  let dux_dy = (top.ux - bottom.ux) * 0.5
  let duy_dx = (right.uy - left.uy) * 0.5
  let shear = (dux_dy + duy_dx) * 0.5
  (2.0 * dux_dx * dux_dx + 2.0 * duy_dy * duy_dy + 4.0 * shear * shear).sqrt()
}

///|
/// Return the divergence field.
pub fn Simulation::divergence_field(self : Simulation) -> Field2D {
  let result = Field2D::new(size=self.size)
  for y in 0.. Field2D {
  self.density_field().map((_, _, rho) => rho / 3.0)
}

///|
/// Return the strain-rate magnitude field.
pub fn Simulation::strain_rate_field(self : Simulation) -> Field2D {
  let result = Field2D::new(size=self.size)
  for y in 0.. Double {
  self.divergence_field().absolute_integral() /
  self.fluid_count().max(1).to_double()
}

///|
/// Return the maximum strain-rate magnitude.
pub fn Simulation::maximum_strain_rate(self : Simulation) -> Double {
  self.strain_rate_field().statistics().maximum
}