steel_worldgen/noise/

aquifer.rs

//! Noise-based aquifer system for underground fluid placement.
//!
//! Matches vanilla's `Aquifer.NoiseBasedAquifer`. Divides the world into a
//! 16×12×16 grid of aquifer cells, each with a randomly-jittered center point.
//! For each non-solid block, finds the 4 nearest aquifer centers and computes
//! fluid status (water vs lava, surface level) based on noise functions.
//! Barrier pressure between neighboring aquifer cells creates solid rock
//! walls between fluid pockets.

use std::simd::i32x4;

use rustc_hash::FxHashMap;

use crate::density::{ColumnCache, DimensionNoises, NoiseSettings};
use steel_math::{clamp, map, map_clamped};
use steel_registry::{REGISTRY, vanilla_blocks};
use steel_utils::BlockStateId;
use steel_utils::random::name_hash::NameHash;
use steel_utils::random::{PositionalRandom, Random, RandomSplitter};

/// Deferred [`Aquifer`]. Used by `create_structures` so chunks where no structure
/// queries the aquifer skip its (expensive) `max_preliminary_surface_level` scan.
pub struct LazyAquifer<'a, N: DimensionNoises> {
    chunk_min_x: i32,
    chunk_min_z: i32,
    splitter: &'a RandomSplitter,
    noises: &'a N,
    inner: Option<Aquifer<N>>,
}

impl<'a, N: DimensionNoises> LazyAquifer<'a, N> {
    /// Deferred aquifer for the given chunk.
    #[must_use]
    pub const fn new(
        chunk_min_x: i32,
        chunk_min_z: i32,
        splitter: &'a RandomSplitter,
        noises: &'a N,
    ) -> Self {
        Self {
            chunk_min_x,
            chunk_min_z,
            splitter,
            noises,
            inner: None,
        }
    }

    /// Build on first call; `height_cache` is cloned into the aquifer's own cache.
    ///
    /// # Panics
    /// Never — `inner` is initialized above if it was `None`.
    pub fn ensure(&mut self, height_cache: &N::ColumnCache) -> &mut Aquifer<N> {
        if self.inner.is_none() {
            self.inner = Some(Aquifer::<N>::new(
                self.chunk_min_x,
                self.chunk_min_z,
                <N::Settings as NoiseSettings>::MIN_Y,
                <N::Settings as NoiseSettings>::HEIGHT,
                self.splitter,
                self.noises,
                height_cache.clone(),
            ));
        }
        #[expect(clippy::unwrap_used, reason = "just initialized above")]
        self.inner.as_mut().unwrap()
    }
}

// Grid spacing
const Y_SPACING: i32 = 12;

// Jitter range per cell center
const X_RANGE: i32 = 10;
const Y_RANGE: i32 = 9;
const Z_RANGE: i32 = 10;

// Anchor offsets for neighborhood lookup
const SAMPLE_OFFSET_X: i32 = -5;
const SAMPLE_OFFSET_Y: i32 = 1;
const SAMPLE_OFFSET_Z: i32 = -5;

const LAVA_LEVEL: i32 = -54;
/// Sentinel for "no fluid" — well below any real Y coordinate.
const WAY_BELOW_MIN_Y: i32 = -32512;
const FLOWING_UPDATE_SIMILARITY: f64 = 1.0 - ((12 * 12 - 10 * 10) as f64) / 25.0;

/// Chunk offsets (in chunks, ×16 for blocks) used when sampling
/// preliminary surface levels around an aquifer cell center.
const SURFACE_SAMPLING_OFFSETS: [[i32; 2]; 13] = [
    [0, 0],
    [-2, -1],
    [-1, -1],
    [0, -1],
    [1, -1],
    [-3, 0],
    [-2, 0],
    [-1, 0],
    [1, 0],
    [-2, 1],
    [-1, 1],
    [0, 1],
    [1, 1],
];

/// Fluid status at an aquifer cell center.
///
/// Matches vanilla's `Aquifer.FluidStatus` — stores the actual fluid block state
/// rather than a boolean flag, so the aquifer is agnostic to which fluids exist.
#[derive(Clone, Copy, PartialEq, Eq)]
struct FluidStatus {
    /// Y level of the fluid surface (exclusive upper bound).
    fluid_level: i32,
    /// Block state placed below `fluid_level`.
    fluid_type: BlockStateId,
}

impl FluidStatus {
    /// What block is at `block_y`? Returns the fluid type if below the surface,
    /// or `None` for air above the surface.
    const fn at(self, block_y: i32) -> Option<BlockStateId> {
        if block_y < self.fluid_level {
            Some(self.fluid_type)
        } else {
            None
        }
    }
}

/// Result of the aquifer substance check.
pub enum AquiferResult {
    /// Solid block (density > 0 or barrier makes it solid).
    Solid,
    /// Air (no block placed).
    Air,
    /// Fluid block to place.
    Fluid(BlockStateId),
}

/// Column-scan state for the 12 aquifer-neighborhood cells, stored `SoA` so the
/// per-Y distance computation can be SIMD-batched as 3× `i32x4`.
///
/// `compute_substance` is called many times with the same `(world_x, world_z)`
/// and decreasing `world_y` (innermost loop in `noise_chunk::fill`). Within a
/// stable `y_anchor` window (`Y_SPACING` blocks tall) the 12 cells of the
/// neighborhood don't change, and `dx*dx + dz*dz` only depends on x/z — only
/// `dy = loc_y - world_y` varies per call.
///
/// `y_anchor = i32::MIN` marks the cache as invalid (forces a refill).
struct AquiferColumnCache {
    world_x: i32,
    world_z: i32,
    y_anchor: i32,
    /// Per-cell unpacked Y of the aquifer-cell center (constant while cached).
    /// Padded to 16 entries so the 12 valid cells fit cleanly into 3× i32x4
    /// SIMD batches; trailing slots stay at default `0`.
    cell_loc_y: [i32; 16],
    /// Per-cell `(dx + fx)² + (dz + fz)²` (Y-independent component of distance).
    cell_xz_dist_sq: [i32; 16],
    /// Per-cell index into `location_cache` / `status_cache`.
    cell_idx: [u32; 12],
}

impl Default for AquiferColumnCache {
    fn default() -> Self {
        Self {
            world_x: 0,
            world_z: 0,
            y_anchor: i32::MIN,
            cell_loc_y: [0; 16],
            cell_xz_dist_sq: [0; 16],
            cell_idx: [0; 12],
        }
    }
}

/// Noise-based aquifer for a single chunk.
///
/// Constructed once per chunk, used throughout the fill loop.
pub struct Aquifer<N: DimensionNoises> {
    /// Packed (x, y, z) locations of aquifer cell centers.
    /// `i64::MAX` = not yet computed.
    location_cache: Vec<i64>,
    /// Lazily computed fluid statuses per grid cell.
    status_cache: Vec<Option<FluidStatus>>,
    /// Positional random for grid cell jitter.
    splitter: RandomSplitter,
    /// Column cache owned by the aquifer for density function evaluation.
    cache: N::ColumnCache,
    /// Grid bounds.
    min_grid_x: i32,
    min_grid_y: i32,
    min_grid_z: i32,
    grid_size_x: i32,
    grid_size_z: i32,
    /// Skip aquifer sampling above this Y (optimization).
    skip_sampling_above_y: i32,
    /// Sea level for this dimension.
    sea_level: i32,
    /// Precomputed `min(LAVA_LEVEL, sea_level)`. Hoisted out of `global_fluid`
    /// (which is on the per-block hot path) so we don't recompute it millions
    /// of times per chunk.
    lava_floor: i32,
    /// Block state IDs.
    water_id: BlockStateId,
    lava_id: BlockStateId,
    /// The dimension's default fluid (water for overworld, lava for nether).
    default_fluid_id: BlockStateId,
    /// Vanilla's `shouldScheduleFluidUpdate` flag from the most recent substance lookup.
    should_schedule_fluid_update: bool,
    /// 12-cell neighborhood snapshot for the current Y-column.
    /// Placed at the end so dimensions with disabled aquifers (nether/end)
    /// keep the hot fluid-id fields earlier in the struct's cache lines.
    col_cache: AquiferColumnCache,
    /// Per-quart-column cache of `preliminary_surface_level` results, matching
    /// vanilla's `NoiseBasedAquifer.preliminarySurfaceLevel` `Long2IntMap`.
    /// `compute_fluid` samples surface level 13× per aquifer cell, and each miss
    /// recomputes the entire flat `NormalNoise` router for that column via
    /// `cache.ensure`. Memoizing the `i32` result per column collapses that to
    /// one evaluation per unique column for the chunk.
    prelim_cache: FxHashMap<(i32, i32), i32>,
}

// Grid coordinate conversions
#[inline]
const fn grid_x(block: i32) -> i32 {
    block >> 4
}
#[inline]
const fn grid_z(block: i32) -> i32 {
    block >> 4
}
#[inline]
const fn grid_y(block: i32) -> i32 {
    block.div_euclid(Y_SPACING)
}
#[inline]
const fn from_grid_x(grid: i32, offset: i32) -> i32 {
    (grid << 4) + offset
}
#[inline]
const fn from_grid_y(grid: i32, offset: i32) -> i32 {
    grid * Y_SPACING + offset
}
#[inline]
const fn from_grid_z(grid: i32, offset: i32) -> i32 {
    (grid << 4) + offset
}

// BlockPos packing (matches vanilla's BlockPos.asLong / getX / getY / getZ)
const PACKED_X_MASK: i64 = 0x3FF_FFFF; // 26 bits
const PACKED_Y_MASK: i64 = 0xFFF; // 12 bits
const PACKED_Z_MASK: i64 = 0x3FF_FFFF; // 26 bits
const X_OFFSET: i32 = 38;
const Z_OFFSET: i32 = 12;

#[inline]
fn pack_pos(x: i32, y: i32, z: i32) -> i64 {
    ((i64::from(x) & PACKED_X_MASK) << X_OFFSET)
        | (i64::from(y) & PACKED_Y_MASK)
        | ((i64::from(z) & PACKED_Z_MASK) << Z_OFFSET)
}

#[inline]
const fn unpack_x(packed: i64) -> i32 {
    (packed >> X_OFFSET) as i32
}

#[inline]
const fn unpack_y(packed: i64) -> i32 {
    ((packed << 52) >> 52) as i32
}

#[inline]
const fn unpack_z(packed: i64) -> i32 {
    ((packed << 26) >> X_OFFSET) as i32
}

/// Similarity between two squared distances. Positive when the two nearest
/// aquifer cells are close together (near a boundary).
#[inline]
fn similarity(dist_sq1: i32, dist_sq2: i32) -> f64 {
    1.0 - f64::from(dist_sq2 - dist_sq1) / 25.0
}

/// Deep dark region check matching `OverworldBiomeBuilder.isDeepDarkRegion`.
fn is_deep_dark_region<N: DimensionNoises>(
    noises: &N,
    cache: &mut N::ColumnCache,
    x: i32,
    y: i32,
    z: i32,
) -> bool {
    cache.ensure(x, z, noises);
    let erosion = noises.router_erosion(cache, x, y, z);
    let depth = noises.router_depth(cache, x, y, z);
    erosion < -0.225 && depth > 0.9
}

/// Global fluid picker matching vanilla's `NoiseBasedChunkGenerator.createFluidPicker`.
///
/// Below `min(-54, sea_level)` → lava at Y=-54. Otherwise → the dimension's
/// default fluid at sea level (water for overworld, lava for nether).
const fn global_fluid(
    y: i32,
    lava_floor: i32,
    sea_level: i32,
    lava_id: BlockStateId,
    default_fluid_id: BlockStateId,
) -> FluidStatus {
    if y < lava_floor {
        FluidStatus {
            fluid_level: LAVA_LEVEL,
            fluid_type: lava_id,
        }
    } else {
        FluidStatus {
            fluid_level: sea_level,
            fluid_type: default_fluid_id,
        }
    }
}

impl<N: DimensionNoises> Aquifer<N> {
    /// Create an aquifer for a full 16×16 chunk.
    ///
    /// `chunk_min_x/z` are the block coordinates of the chunk's NW corner.
    /// `min_block_y` and `y_block_size` define the vertical range.
    /// `splitter` is the seed's positional splitter.
    /// `cache` should be a pre-initialized column cache for this chunk
    /// (avoids a redundant `init_grid` call).
    #[must_use]
    pub fn new(
        chunk_min_x: i32,
        chunk_min_z: i32,
        min_block_y: i32,
        y_block_size: i32,
        splitter: &RandomSplitter,
        noises: &N,
        cache: N::ColumnCache,
    ) -> Self {
        Self::new_sized(
            chunk_min_x,
            chunk_min_z,
            16,
            16,
            min_block_y,
            y_block_size,
            splitter,
            noises,
            cache,
        )
    }

    /// Create an aquifer with custom XZ extent (in blocks).
    /// Vanilla's iterateNoiseColumn uses width=cellWidth (4) for single-column queries.
    #[expect(
        clippy::too_many_arguments,
        reason = "mirrors vanilla's Aquifer constructor shape"
    )]
    pub fn new_sized(
        chunk_min_x: i32,
        chunk_min_z: i32,
        width_x: i32,
        width_z: i32,
        min_block_y: i32,
        y_block_size: i32,
        splitter: &RandomSplitter,
        noises: &N,
        mut cache: N::ColumnCache,
    ) -> Self {
        const AQUIFER_HASH: NameHash = NameHash::new("minecraft:aquifer");

        let sea_level = N::Settings::SEA_LEVEL;
        let water_id = REGISTRY.blocks.get_default_state_id(&vanilla_blocks::WATER);
        let lava_id = REGISTRY.blocks.get_default_state_id(&vanilla_blocks::LAVA);
        let default_fluid_id = N::Settings::default_fluid_id();

        let mut aquifer_rng = splitter.with_hash_of(&AQUIFER_HASH);
        let splitter = aquifer_rng.next_positional();

        // When aquifers are disabled (nether/end), compute_substance uses only
        // the global fluid picker — skip grid allocation and surface sampling.
        if !N::Settings::AQUIFERS_ENABLED {
            return Self {
                location_cache: Vec::new(),
                status_cache: Vec::new(),
                splitter,
                cache,
                col_cache: AquiferColumnCache::default(),
                min_grid_x: 0,
                min_grid_y: 0,
                min_grid_z: 0,
                grid_size_x: 0,
                grid_size_z: 0,
                skip_sampling_above_y: 0,
                sea_level,
                lava_floor: LAVA_LEVEL.min(sea_level),
                water_id,
                lava_id,
                default_fluid_id,
                should_schedule_fluid_update: false,
                prelim_cache: FxHashMap::default(),
            };
        }

        let chunk_max_x = chunk_min_x + width_x - 1;
        let chunk_max_z = chunk_min_z + width_z - 1;

        let min_grid_x = grid_x(chunk_min_x + SAMPLE_OFFSET_X);
        let max_grid_x = grid_x(chunk_max_x + SAMPLE_OFFSET_X) + 1;
        let grid_size_x = max_grid_x - min_grid_x + 1;

        let min_grid_y = grid_y(min_block_y + SAMPLE_OFFSET_Y) - 1;
        let max_grid_y = grid_y(min_block_y + y_block_size + SAMPLE_OFFSET_Y) + 1;
        let grid_size_y = max_grid_y - min_grid_y + 1;

        let min_grid_z = grid_z(chunk_min_z + SAMPLE_OFFSET_Z);
        let max_grid_z = grid_z(chunk_max_z + SAMPLE_OFFSET_Z) + 1;
        let grid_size_z = max_grid_z - min_grid_z + 1;

        let total = (grid_size_x * grid_size_y * grid_size_z) as usize;
        let location_cache = vec![i64::MAX; total];
        let status_cache = vec![None; total];

        // Compute skip_sampling_above_y from max preliminary surface level.
        // The scan primes `prelim_cache` for the columns `compute_fluid` reuses.
        let mut prelim_cache = FxHashMap::default();
        let max_surface = Self::max_preliminary_surface_level(
            noises,
            &mut cache,
            &mut prelim_cache,
            from_grid_x(min_grid_x, 0),
            from_grid_z(min_grid_z, 0),
            from_grid_x(max_grid_x, X_RANGE - 1),
            from_grid_z(max_grid_z, Z_RANGE - 1),
        );
        let adjusted = max_surface + 8;
        let skip_grid_y = grid_y(adjusted + 12) + 1;
        let skip_sampling_above_y = from_grid_y(skip_grid_y, Y_RANGE + 2) - 1;

        Self {
            location_cache,
            status_cache,
            splitter,
            cache,
            col_cache: AquiferColumnCache::default(),
            min_grid_x,
            min_grid_y,
            min_grid_z,
            grid_size_x,
            grid_size_z,
            skip_sampling_above_y,
            sea_level,
            lava_floor: LAVA_LEVEL.min(sea_level),
            water_id,
            lava_id,
            default_fluid_id,
            should_schedule_fluid_update: false,
            prelim_cache,
        }
    }

    fn max_preliminary_surface_level(
        noises: &N,
        cache: &mut N::ColumnCache,
        prelim_cache: &mut FxHashMap<(i32, i32), i32>,
        min_x: i32,
        min_z: i32,
        max_x: i32,
        max_z: i32,
    ) -> i32 {
        let mut max_level = i32::MIN;
        // Sample at 4-block intervals (quart-block resolution) across the chunk area
        let mut z = min_z;
        while z <= max_z {
            let mut x = min_x;
            while x <= max_x {
                let level = cached_preliminary_surface_level(noises, cache, prelim_cache, x, z);
                if level > max_level {
                    max_level = level;
                }
                x += 4;
            }
            z += 4;
        }
        max_level
    }

    #[inline]
    const fn get_index(&self, gx: i32, gy: i32, gz: i32) -> usize {
        let x = gx - self.min_grid_x;
        let y = gy - self.min_grid_y;
        let z = gz - self.min_grid_z;
        ((y * self.grid_size_z + z) * self.grid_size_x + x) as usize
    }

    /// Refill the 12-cell column cache for the current `(world_x, world_z, y_anchor)`.
    ///
    /// Iterates the same `(x1, y1, z1)` order as the inline scan so cells are
    /// stored at consistent indices, preserving tie-breaking when the per-Y
    /// `new_dist` values are compared in `compute_substance`.
    fn refill_col_cache(&mut self, world_x: i32, world_y: i32, world_z: i32) -> i32 {
        let x_anchor = grid_x(world_x + SAMPLE_OFFSET_X);
        let y_anchor = grid_y(world_y + SAMPLE_OFFSET_Y);
        let z_anchor = grid_z(world_z + SAMPLE_OFFSET_Z);

        let mut i = 0;
        for x1 in 0..=1i32 {
            for y1 in -1..=1i32 {
                for z1 in 0..=1i32 {
                    let gx = x_anchor + x1;
                    let gy = y_anchor + y1;
                    let gz = z_anchor + z1;
                    let idx = self.get_index(gx, gy, gz);

                    let loc = self.location_cache[idx];
                    let loc = if loc == i64::MAX {
                        let mut rng = self.splitter.at(gx, gy, gz);
                        let packed = pack_pos(
                            from_grid_x(gx, rng.next_i32_bounded(X_RANGE)),
                            from_grid_y(gy, rng.next_i32_bounded(Y_RANGE)),
                            from_grid_z(gz, rng.next_i32_bounded(Z_RANGE)),
                        );
                        self.location_cache[idx] = packed;
                        packed
                    } else {
                        loc
                    };

                    let dx = unpack_x(loc) - world_x;
                    let dz = unpack_z(loc) - world_z;
                    self.col_cache.cell_loc_y[i] = unpack_y(loc);
                    self.col_cache.cell_xz_dist_sq[i] = dx * dx + dz * dz;
                    self.col_cache.cell_idx[i] = idx as u32;
                    i += 1;
                }
            }
        }

        self.col_cache.world_x = world_x;
        self.col_cache.world_z = world_z;
        self.col_cache.y_anchor = y_anchor;
        y_anchor
    }

    /// Compute what block to place at this position given the interpolated density.
    #[expect(
        clippy::too_many_lines,
        reason = "splitting would hurt readability of the aquifer sampling logic"
    )]
    pub fn compute_substance(
        &mut self,
        noises: &N,
        world_x: i32,
        world_y: i32,
        world_z: i32,
        density: f64,
    ) -> AquiferResult {
        // Solid block — let the caller decide (stone or ore)
        if density > 0.0 {
            self.should_schedule_fluid_update = false;
            return AquiferResult::Solid;
        }

        // Disabled aquifers (nether/end): use global fluid picker directly,
        // matching vanilla's `Aquifer.createDisabled`.
        if !N::Settings::AQUIFERS_ENABLED {
            self.should_schedule_fluid_update = false;
            let gf = global_fluid(
                world_y,
                self.lava_floor,
                self.sea_level,
                self.lava_id,
                self.default_fluid_id,
            );
            return match gf.at(world_y) {
                Some(id) => AquiferResult::Fluid(id),
                None => AquiferResult::Air,
            };
        }

        let gf = global_fluid(
            world_y,
            self.lava_floor,
            self.sea_level,
            self.lava_id,
            self.default_fluid_id,
        );

        // Above the skip threshold: use global fluid directly
        if world_y > self.skip_sampling_above_y {
            self.should_schedule_fluid_update = false;
            return match gf.at(world_y) {
                Some(id) => AquiferResult::Fluid(id),
                None => AquiferResult::Air,
            };
        }

        // If global fluid is lava here, return lava
        if gf.fluid_type == self.lava_id && world_y < gf.fluid_level {
            self.should_schedule_fluid_update = false;
            return AquiferResult::Fluid(self.lava_id);
        }

        // Find 4 nearest aquifer cell centers from the 2×3×2 neighborhood.
        // Within a Y-column scan, the 12 cells and their `xz_dist_sq` values
        // are constant for a given `y_anchor`; only `dy` varies. The column
        // cache amortizes location lookup, splitter calls, and i64 unpacking.
        let y_anchor = grid_y(world_y + SAMPLE_OFFSET_Y);
        if self.col_cache.world_x != world_x
            || self.col_cache.world_z != world_z
            || self.col_cache.y_anchor != y_anchor
        {
            self.refill_col_cache(world_x, world_y, world_z);
        }

        // SIMD-batch the per-cell distance computation: `loc_y - world_y` then
        // `xz_dist_sq + dy²`, processing 4 cells per `i32x4` op. The 12 valid
        // cells fit in 3 batches; the 4-slot tail of `cell_loc_y` /
        // `cell_xz_dist_sq` is harmless padding (we ignore the trailing slot).
        let world_y_v = i32x4::splat(world_y);
        let mut dists = [0i32; 12];
        for batch in 0..3 {
            let base = batch * 4;
            let loc_y_v = i32x4::from_slice(&self.col_cache.cell_loc_y[base..base + 4]);
            let xz_v = i32x4::from_slice(&self.col_cache.cell_xz_dist_sq[base..base + 4]);
            let dy = loc_y_v - world_y_v;
            let dist_v = xz_v + dy * dy;
            dists[base..base + 4].copy_from_slice(&dist_v.to_array());
        }

        let mut dist_sq = [i32::MAX; 4];
        let mut closest_idx = [0usize; 4];

        for (i, &new_dist) in dists.iter().enumerate() {
            let index = self.col_cache.cell_idx[i] as usize;

            // Insert into sorted top-4
            if dist_sq[0] >= new_dist {
                dist_sq[3] = dist_sq[2];
                closest_idx[3] = closest_idx[2];
                dist_sq[2] = dist_sq[1];
                closest_idx[2] = closest_idx[1];
                dist_sq[1] = dist_sq[0];
                closest_idx[1] = closest_idx[0];
                dist_sq[0] = new_dist;
                closest_idx[0] = index;
            } else if dist_sq[1] >= new_dist {
                dist_sq[3] = dist_sq[2];
                closest_idx[3] = closest_idx[2];
                dist_sq[2] = dist_sq[1];
                closest_idx[2] = closest_idx[1];
                dist_sq[1] = new_dist;
                closest_idx[1] = index;
            } else if dist_sq[2] >= new_dist {
                dist_sq[3] = dist_sq[2];
                closest_idx[3] = closest_idx[2];
                dist_sq[2] = new_dist;
                closest_idx[2] = index;
            } else if dist_sq[3] >= new_dist {
                dist_sq[3] = new_dist;
                closest_idx[3] = index;
            }
        }

        let status1 = self.get_aquifer_status(closest_idx[0], noises);
        let fluid_at = status1.at(world_y);

        // `similarity(d1, d2) = 1 - (d2 - d1) / 25`, so `sim12 <= 0.0` is exactly
        // `d2 - d1 >= 25` in i32. Defer the f64 conversion + divide until after
        // the early-return check. Fluid-update scheduling still matches vanilla:
        // `sim12 >= FLOWING_UPDATE_SIMILARITY` is exactly `d2 - d1 <= 44`.
        let dist12_delta = dist_sq[1] - dist_sq[0];
        if dist12_delta >= 25 {
            if dist12_delta <= 12 * 12 - 10 * 10 {
                let status2 = self.get_aquifer_status(closest_idx[1], noises);
                self.should_schedule_fluid_update = status1 != status2;
            } else {
                self.should_schedule_fluid_update = false;
            }
            return match fluid_at {
                Some(id) => AquiferResult::Fluid(id),
                None => AquiferResult::Air,
            };
        }
        let sim12 = similarity(dist_sq[0], dist_sq[1]);

        // Water adjacent to global lava below → return water
        if let Some(id) = fluid_at
            && id == self.water_id
        {
            let below = global_fluid(
                world_y - 1,
                self.lava_floor,
                self.sea_level,
                self.lava_id,
                self.default_fluid_id,
            );
            if below.fluid_type == self.lava_id && (world_y - 1) < below.fluid_level {
                self.should_schedule_fluid_update = true;
                return AquiferResult::Fluid(id);
            }
        }

        // Compute barrier pressure between closest pairs
        let mut barrier_noise = f64::NAN;
        let status2 = self.get_aquifer_status(closest_idx[1], noises);
        let barrier12 = sim12
            * self.calculate_pressure(
                noises,
                world_x,
                world_y,
                world_z,
                &mut barrier_noise,
                status1,
                status2,
            );
        if density + barrier12 > 0.0 {
            self.should_schedule_fluid_update = false;
            return AquiferResult::Solid;
        }

        let status3 = self.get_aquifer_status(closest_idx[2], noises);
        let sim13 = similarity(dist_sq[0], dist_sq[2]);
        if sim13 > 0.0 {
            let barrier13 = sim12
                * sim13
                * self.calculate_pressure(
                    noises,
                    world_x,
                    world_y,
                    world_z,
                    &mut barrier_noise,
                    status1,
                    status3,
                );
            if density + barrier13 > 0.0 {
                self.should_schedule_fluid_update = false;
                return AquiferResult::Solid;
            }
        }

        let sim23 = similarity(dist_sq[1], dist_sq[2]);
        if sim23 > 0.0 {
            let barrier23 = sim12
                * sim23
                * self.calculate_pressure(
                    noises,
                    world_x,
                    world_y,
                    world_z,
                    &mut barrier_noise,
                    status2,
                    status3,
                );
            if density + barrier23 > 0.0 {
                self.should_schedule_fluid_update = false;
                return AquiferResult::Solid;
            }
        }

        let may_flow12 = status1 != status2;
        let may_flow23 = sim23 >= FLOWING_UPDATE_SIMILARITY && status2 != status3;
        let may_flow13 = sim13 >= FLOWING_UPDATE_SIMILARITY && status1 != status3;
        if may_flow12 || may_flow23 || may_flow13 {
            self.should_schedule_fluid_update = true;
        } else {
            self.should_schedule_fluid_update = sim13 >= FLOWING_UPDATE_SIMILARITY
                && similarity(dist_sq[0], dist_sq[3]) >= FLOWING_UPDATE_SIMILARITY
                && status1 != self.get_aquifer_status(closest_idx[3], noises);
        }

        // Return the closest fluid
        match fluid_at {
            Some(id) => AquiferResult::Fluid(id),
            None => AquiferResult::Air,
        }
    }

    /// Returns whether the most recent substance lookup needs postprocessing for placed fluids.
    #[must_use]
    pub const fn should_schedule_fluid_update(&self) -> bool {
        self.should_schedule_fluid_update
    }

    /// Returns the quart-quantized preliminary surface level, reusing this aquifer's
    /// density-column and result caches.
    pub fn preliminary_surface_level(&mut self, noises: &N, x: i32, z: i32) -> i32 {
        cached_preliminary_surface_level(noises, &mut self.cache, &mut self.prelim_cache, x, z)
    }

    /// Get or compute the fluid status for the aquifer cell at the given cache index.
    fn get_aquifer_status(&mut self, index: usize, noises: &N) -> FluidStatus {
        if let Some(status) = self.status_cache[index] {
            return status;
        }

        let loc = self.location_cache[index];
        let x = unpack_x(loc);
        let y = unpack_y(loc);
        let z = unpack_z(loc);
        let status = self.compute_fluid(x, y, z, noises);
        self.status_cache[index] = Some(status);
        status
    }

    /// Compute the fluid status for an aquifer cell centered at (x, y, z).
    fn compute_fluid(&mut self, x: i32, y: i32, z: i32, noises: &N) -> FluidStatus {
        let gf = global_fluid(
            y,
            self.lava_floor,
            self.sea_level,
            self.lava_id,
            self.default_fluid_id,
        );
        let mut lowest_surface = i32::MAX;
        let top_of_cell = y + Y_SPACING;
        let bottom_of_cell = y - Y_SPACING;
        let mut surface_under_global = false;

        for offset in &SURFACE_SAMPLING_OFFSETS {
            let sx = x + offset[0] * 16; // sectionToBlockCoord
            let sz = z + offset[1] * 16;

            let preliminary = cached_preliminary_surface_level(
                noises,
                &mut self.cache,
                &mut self.prelim_cache,
                sx,
                sz,
            );
            let adjusted = preliminary + 8;

            let is_center = offset[0] == 0 && offset[1] == 0;

            if is_center && bottom_of_cell > adjusted {
                return gf;
            }

            let top_pokes_above = top_of_cell > adjusted;
            if top_pokes_above || is_center {
                let gf_at_surface = global_fluid(
                    adjusted,
                    self.lava_floor,
                    self.sea_level,
                    self.lava_id,
                    self.default_fluid_id,
                );
                let has_fluid = adjusted < gf_at_surface.fluid_level;
                if has_fluid {
                    if is_center {
                        surface_under_global = true;
                    }
                    if top_pokes_above {
                        return gf_at_surface;
                    }
                }
            }

            if preliminary < lowest_surface {
                lowest_surface = preliminary;
            }
        }

        let fluid_level =
            self.compute_surface_level(x, y, z, noises, gf, lowest_surface, surface_under_global);
        let fluid_type = self.compute_fluid_type(x, y, z, noises, gf, fluid_level);
        FluidStatus {
            fluid_level,
            fluid_type,
        }
    }

    #[expect(
        clippy::too_many_arguments,
        reason = "matches vanilla NoiseBasedAquifer.computeSurface signature"
    )]
    fn compute_surface_level(
        &mut self,
        x: i32,
        y: i32,
        z: i32,
        noises: &N,
        gf: FluidStatus,
        lowest_surface: i32,
        surface_under_global: bool,
    ) -> i32 {
        let (partially_flooded, fully_flooded) =
            if is_deep_dark_region(noises, &mut self.cache, x, y, z) {
                (-1.0, -1.0)
            } else {
                let dist_below = lowest_surface + 8 - y;
                let floodedness_factor = if surface_under_global {
                    map_clamped(f64::from(dist_below), 0.0, 64.0, 1.0, 0.0)
                } else {
                    0.0
                };

                self.cache.ensure(x, z, noises);
                let floodedness_noise = clamp(
                    noises.router_fluid_level_floodedness(&mut self.cache, x, y, z),
                    -1.0,
                    1.0,
                );

                let fully_threshold = map(floodedness_factor, 1.0, 0.0, -0.3, 0.8);
                let partially_threshold = map(floodedness_factor, 1.0, 0.0, -0.8, 0.4);

                (
                    floodedness_noise - partially_threshold,
                    floodedness_noise - fully_threshold,
                )
            };

        if fully_flooded > 0.0 {
            gf.fluid_level
        } else if partially_flooded > 0.0 {
            self.compute_randomized_fluid_surface_level(x, y, z, noises, lowest_surface)
        } else {
            WAY_BELOW_MIN_Y
        }
    }

    fn compute_randomized_fluid_surface_level(
        &mut self,
        x: i32,
        y: i32,
        z: i32,
        noises: &N,
        lowest_surface: i32,
    ) -> i32 {
        let cell_x = x.div_euclid(16);
        let cell_y = y.div_euclid(40);
        let cell_z = z.div_euclid(16);
        let cell_middle_y = cell_y * 40 + 20;

        // fluid_level_spread is evaluated at grid coordinates (not block coordinates)
        self.cache.ensure(cell_x, cell_z, noises);
        let spread =
            noises.router_fluid_level_spread(&mut self.cache, cell_x, cell_y, cell_z) * 10.0;
        let spread_quantized = quantize(spread, 3);
        let target = cell_middle_y + spread_quantized;

        lowest_surface.min(target)
    }

    fn compute_fluid_type(
        &mut self,
        x: i32,
        y: i32,
        z: i32,
        noises: &N,
        gf: FluidStatus,
        fluid_level: i32,
    ) -> BlockStateId {
        if fluid_level <= -10 && fluid_level != WAY_BELOW_MIN_Y && gf.fluid_type != self.lava_id {
            let cell_x = x.div_euclid(64);
            let cell_y = y.div_euclid(40);
            let cell_z = z.div_euclid(64);
            self.cache.ensure(cell_x, cell_z, noises);
            let lava_noise = noises.router_lava(&mut self.cache, cell_x, cell_y, cell_z);
            if lava_noise.abs() > 0.3 {
                return self.lava_id;
            }
        }
        gf.fluid_type
    }

    /// Calculate barrier pressure between two aquifer cells.
    ///
    /// Matches vanilla's check: if lava meets water at this Y, return max pressure.
    #[expect(
        clippy::too_many_arguments,
        reason = "matches vanilla NoiseBasedAquifer.calculatePressure signature"
    )]
    fn calculate_pressure(
        &mut self,
        noises: &N,
        x: i32,
        y: i32,
        z: i32,
        barrier_noise: &mut f64,
        s1: FluidStatus,
        s2: FluidStatus,
    ) -> f64 {
        let f1 = s1.at(y);
        let f2 = s2.at(y);
        let f1_is_lava = f1 == Some(self.lava_id);
        let f2_is_lava = f2 == Some(self.lava_id);
        let f1_is_water = f1 == Some(self.water_id);
        let f2_is_water = f2 == Some(self.water_id);

        // Lava–water interface → max pressure
        if (f1_is_lava && f2_is_water) || (f1_is_water && f2_is_lava) {
            return 2.0;
        }

        let fluid_y_diff = (s1.fluid_level - s2.fluid_level).abs();
        if fluid_y_diff == 0 {
            return 0.0;
        }

        let avg_fluid_y = 0.5 * f64::from(s1.fluid_level + s2.fluid_level);
        let above_avg = f64::from(y) + 0.5 - avg_fluid_y;
        let base = f64::from(fluid_y_diff) / 2.0;
        let edge_dist = base - above_avg.abs();

        let gradient = if above_avg > 0.0 {
            if edge_dist > 0.0 {
                edge_dist / 1.5
            } else {
                edge_dist / 2.5
            }
        } else {
            let center = 3.0 + edge_dist;
            if center > 0.0 {
                center / 3.0
            } else {
                center / 10.0
            }
        };

        let noise_val = if !(-2.0..=2.0).contains(&gradient) {
            0.0
        } else if barrier_noise.is_nan() {
            self.cache.ensure(x, z, noises);
            let n = noises.router_barrier(&mut self.cache, x, y, z);
            *barrier_noise = n;
            n
        } else {
            *barrier_noise
        };

        2.0 * (noise_val + gradient)
    }
}

/// Quantize: snap value down to the nearest multiple of `quantum`.
#[inline]
fn quantize(value: f64, quantum: i32) -> i32 {
    let q = f64::from(quantum);
    (value / q).floor() as i32 * quantum
}

/// Evaluate preliminary surface level at quart-quantized coordinates.
///
/// Vanilla's `NoiseChunk.preliminarySurfaceLevel()` quantizes X/Z to quart
/// positions before lookup, matching `FlatCache`'s 4-block grid.
pub fn preliminary_surface_level<N: DimensionNoises>(
    noises: &N,
    cache: &mut N::ColumnCache,
    x: i32,
    z: i32,
) -> i32 {
    // Quantize to quart positions: (x >> 2) << 2
    let qx = (x >> 2) << 2;
    let qz = (z >> 2) << 2;
    cache.ensure(qx, qz, noises);
    // Vanilla uses Mth.floor(), not truncation
    noises
        .router_preliminary_surface_level(cache, qx, 0, qz)
        .floor() as i32
}

/// [`preliminary_surface_level`] with a per-quart-column result cache (vanilla's
/// `Long2IntMap`). On a hit it returns the memoized `i32` and skips the
/// expensive flat-router recompute in `cache.ensure`. Bit-identical: the cached
/// value is the same deterministic function of the quart column.
fn cached_preliminary_surface_level<N: DimensionNoises>(
    noises: &N,
    cache: &mut N::ColumnCache,
    prelim_cache: &mut FxHashMap<(i32, i32), i32>,
    x: i32,
    z: i32,
) -> i32 {
    let key = ((x >> 2) << 2, (z >> 2) << 2);
    if let Some(&level) = prelim_cache.get(&key) {
        return level;
    }
    let level = preliminary_surface_level(noises, cache, x, z);
    prelim_cache.insert(key, level);
    level
}