///|
/// Boundary behavior for a passive scalar.
pub(all) enum ScalarBoundary {
  ZeroFlux
  Periodic
  Fixed(Double)
} derive(Debug)

///|
/// Explicit finite-difference advection-diffusion solver for one scalar.
pub(all) struct ScalarSolver {
  size : Size
  diffusivity : Double
  dt : Double
  boundary : ScalarBoundary
  values : Array[Double]
  next : Array[Double]
  source : Array[Double]
  mut step_count : Int
} derive(Debug)

///|
/// Allocate a scalar transport field.
pub fn ScalarSolver::new(
  size~ : Size,
  diffusivity~ : Double,
  dt~ : Double,
  boundary? : ScalarBoundary = ZeroFlux,
) -> ScalarSolver {
  let cells = size.width.max(0) * size.height.max(0)
  {
    size,
    diffusivity: diffusivity.max(0.0),
    dt: dt.max(0.0),
    boundary,
    values: Array::make(cells, 0.0),
    next: Array::make(cells, 0.0),
    source: Array::make(cells, 0.0),
    step_count: 0,
  }
}

///|
/// Return the field size.
pub fn ScalarSolver::domain(self : ScalarSolver) -> Size {
  self.size
}

///|
/// Set one scalar value.
pub fn ScalarSolver::set(
  self : ScalarSolver,
  x~ : Int,
  y~ : Int,
  value~ : Double,
) -> Unit {
  if x >= 0 && y >= 0 && x < self.size.width && y < self.size.height {
    self.values[y * self.size.width + x] = value
  }
}

///|
/// Read one scalar value with boundary-aware extension.
pub fn ScalarSolver::get(self : ScalarSolver, x~ : Int, y~ : Int) -> Double {
  if x >= 0 && y >= 0 && x < self.size.width && y < self.size.height {
    self.values[y * self.size.width + x]
  } else {
    match self.boundary {
      Periodic =>
        if self.size.width == 0 || self.size.height == 0 {
          0.0
        } else {
          self.values[wrap_coordinate(y, self.size.height) * self.size.width +
          wrap_coordinate(x, self.size.width)]
        }
      Fixed(value) => value
      ZeroFlux => {
        let xx = x.clamp(min=0, max=self.size.width - 1)
        let yy = y.clamp(min=0, max=self.size.height - 1)
        if self.size.width == 0 || self.size.height == 0 {
          0.0
        } else {
          self.values[yy * self.size.width + xx]
        }
      }
    }
  }
}

///|
/// Fill the scalar field.
pub fn ScalarSolver::fill(self : ScalarSolver, value~ : Double) -> Unit {
  self.values.fill(value)
}

///|
/// Set a per-cell source term.
pub fn ScalarSolver::set_source(
  self : ScalarSolver,
  x~ : Int,
  y~ : Int,
  value~ : Double,
) -> Unit {
  if x >= 0 && y >= 0 && x < self.size.width && y < self.size.height {
    self.source[y * self.size.width + x] = value
  }
}

///|
/// Clear all source terms.
pub fn ScalarSolver::clear_sources(self : ScalarSolver) -> Unit {
  self.source.fill(0.0)
}

///|
/// Set a uniform velocity and advance one advection-diffusion step.
pub fn ScalarSolver::step_with_velocity(
  self : ScalarSolver,
  velocity~ : VectorField2D,
) -> Unit {
  self.next.fill(0.0)
  for y in 0.. Unit {
  self.step_with_velocity(velocity=VectorField2D::new(size=self.size))
}

///|
/// Convert the scalar values to a reusable field object.
pub fn ScalarSolver::field(self : ScalarSolver) -> Field2D {
  { size: self.size, data: self.values.copy() }
}

///|
/// Return the mean scalar value.
pub fn ScalarSolver::mean(self : ScalarSolver) -> Double {
  mean_value(self.values)
}

///|
/// Return the minimum scalar value.
pub fn ScalarSolver::minimum(self : ScalarSolver) -> Double {
  min_max(self.values).0
}

///|
/// Return the maximum scalar value.
pub fn ScalarSolver::maximum(self : ScalarSolver) -> Double {
  min_max(self.values).1
}

///|
/// Return descriptive scalar statistics.
pub fn ScalarSolver::statistics(self : ScalarSolver) -> FieldStatistics {
  self.field().statistics()
}

///|
/// Estimate a scalar Courant number for a uniform velocity.
pub fn ScalarSolver::courant_number(
  self : ScalarSolver,
  ux~ : Double,
  uy~ : Double,
) -> Double {
  (ux.abs() + uy.abs()) * self.dt
}

///|
/// Estimate the explicit diffusion stability number.
pub fn ScalarSolver::diffusion_number(self : ScalarSolver) -> Double {
  self.diffusivity * self.dt
}

///|
/// Apply a source to every cell.
pub fn ScalarSolver::add_uniform_source(
  self : ScalarSolver,
  value~ : Double,
) -> Unit {
  for i in 0..