Water Drop Tank

// Water Drop Tank: D3Q19 free-surface lattice Boltzmann water.
// Maintained directly as ordinary Hyle; validate stencil changes with the LBM tests.
// Source: Thuerey (2007), chapter 4, https://www10.cs.fau.de/publications/dissertations/Diss_2007-Thuerey.pdf
// BGK + Guo gravity, atmospheric reconstruction and conservative interface motion.
// Lattice units: dx=dt=rho_air=1, g=0.0008, tau=0.6, nu=(tau-0.5)/3.
// kind: 0 gas, 1 interface, 2 liquid. Model names record painted defaults only.
// mass is actual liquid; rho is pressure-bearing density, not volume fraction.
// mass + reserve is conserved; reserve carries signed conversion roundoff/excess.
// f0..f18 are populations. Other fields are phase scratch, also included in bakes.
// The finite +/-1e6 guards serve compiler interval checking, not stabilization.
// No surface tension, gas dynamics, calibrated SI units or smooth surface mesh.
property FreeSurfaceWater {
    mass: Float<0.0> ~1e-10;
    reserve: Float<0.0> ~1e-10;
    ux: Float<0.0> ~1e-10;
    uy: Float<0.0> ~1e-10;
    uz: Float<0.0> ~1e-10;
    share: Float<0.0> ~1e-10;
    carried: Float<0.0> ~1e-10;
    gx: Float<0.0> ~1e-10;
    gy: Float<0.0> ~1e-10;
    gz: Float<0.0> ~1e-10;
    targetMass: Float<0.0> ~1e-10;
    f0: Float<0.3333333333333333> ~1e-10;
    f1: Float<0.027777777777777776> ~1e-10;
    f2: Float<0.027777777777777776> ~1e-10;
    f3: Float<0.05555555555555555> ~1e-10;
    f4: Float<0.027777777777777776> ~1e-10;
    f5: Float<0.027777777777777776> ~1e-10;
    f6: Float<0.027777777777777776> ~1e-10;
    f7: Float<0.05555555555555555> ~1e-10;
    f8: Float<0.027777777777777776> ~1e-10;
    f9: Float<0.05555555555555555> ~1e-10;
    f10: Float<0.05555555555555555> ~1e-10;
    f11: Float<0.027777777777777776> ~1e-10;
    f12: Float<0.05555555555555555> ~1e-10;
    f13: Float<0.027777777777777776> ~1e-10;
    f14: Float<0.027777777777777776> ~1e-10;
    f15: Float<0.027777777777777776> ~1e-10;
    f16: Float<0.05555555555555555> ~1e-10;
    f17: Float<0.027777777777777776> ~1e-10;
    f18: Float<0.027777777777777776> ~1e-10;
    rho: Float<1.0> [0.000001 4.0] ~1e-10;
    gravity: Float<0.0008> [0.0 0.001] ~1e-10;
    tau: Float<0.6> [0.6 2.0] ~1e-10;
    kind: Float<0.0> [0.0 2.0];
    nextKind: Float<0.0> [0.0 2.0];
    rank: Float<0.0> [0.0 2.0];
    wet: Float<0.0> [0.0 1.0];
    fluid: Float<0.0> [0.0 1.0];
    interface: Float<0.0> [0.0 1.0];
    gas: Float<0.0> [0.0 1.0];
    fill: Float<0.0> [0.0 1.0];
    filling: Float<0.0> [0.0 1.0];
    emptying: Float<0.0> [0.0 1.0];
    recipient: Float<0.0> [0.0 1.0];
    unstable: Float<0.0> [0.0 1.0];
}
property TankGlass { thickness: Float<1.0>; }
property TankStone { hardness: Float<1.0>; }
model Glass : TankGlass;
model LbmMedium : FreeSurfaceWater;
model Stone : TankStone;
model Water : FreeSurfaceWater { mass<1.0>; kind<2.0>; fill<1.0>; wet<1.0>; fluid<1.0>; };
neighborhood D3Q19 {
    D1 = [-1, -1, 0];
    D2 = [-1, 0, -1];
    D3 = [-1, 0, 0];
    D4 = [-1, 0, 1];
    D5 = [-1, 1, 0];
    D6 = [0, -1, -1];
    D7 = [0, -1, 0];
    D8 = [0, -1, 1];
    D9 = [0, 0, -1];
    D10 = [0, 0, 1];
    D11 = [0, 1, -1];
    D12 = [0, 1, 0];
    D13 = [0, 1, 1];
    D14 = [1, -1, 0];
    D15 = [1, 0, -1];
    D16 = [1, 0, 0];
    D17 = [1, 0, 1];
    D18 = [1, 1, 0];
}
world { dimensions 3; cell Cube; neighborhood D3Q19; }

behavior RecoverMoments for FreeSurfaceWater {
    update when self.kind > 0.0 {
        let density = self.f0 + self.f1 + self.f2 + self.f3 + self.f4 + self.f5 + self.f6 + self.f7 + self.f8 + self.f9 + self.f10 + self.f11 + self.f12 + self.f13 + self.f14 + self.f15 + self.f16 + self.f17 + self.f18;
        set self.rho = clamp(density, 0.000001, 4.0);
        set self.ux = clamp((0.0 * self.f0 + -1.0 * self.f1 + -1.0 * self.f2 + -1.0 * self.f3 + -1.0 * self.f4 + -1.0 * self.f5 + 0.0 * self.f6 + 0.0 * self.f7 + 0.0 * self.f8 + 0.0 * self.f9 + 0.0 * self.f10 + 0.0 * self.f11 + 0.0 * self.f12 + 0.0 * self.f13 + 1.0 * self.f14 + 1.0 * self.f15 + 1.0 * self.f16 + 1.0 * self.f17 + 1.0 * self.f18) / max(density, 0.000001), -1000000.0, 1000000.0);
        set self.uy = clamp((0.0 * self.f0 + -1.0 * self.f1 + 0.0 * self.f2 + 0.0 * self.f3 + 0.0 * self.f4 + 1.0 * self.f5 + -1.0 * self.f6 + -1.0 * self.f7 + -1.0 * self.f8 + 0.0 * self.f9 + 0.0 * self.f10 + 1.0 * self.f11 + 1.0 * self.f12 + 1.0 * self.f13 + -1.0 * self.f14 + 0.0 * self.f15 + 0.0 * self.f16 + 0.0 * self.f17 + 1.0 * self.f18) / max(density, 0.000001), -1000000.0, 1000000.0);
        set self.uz = clamp((0.0 * self.f0 + 0.0 * self.f1 + -1.0 * self.f2 + 0.0 * self.f3 + 1.0 * self.f4 + 0.0 * self.f5 + -1.0 * self.f6 + 0.0 * self.f7 + 1.0 * self.f8 + -1.0 * self.f9 + 1.0 * self.f10 + -1.0 * self.f11 + 0.0 * self.f12 + 1.0 * self.f13 + 0.0 * self.f14 + -1.0 * self.f15 + 0.0 * self.f16 + 1.0 * self.f17 + 0.0 * self.f18) / max(density, 0.000001) - self.gravity * 0.5, -1000000.0, 1000000.0);
    }
}
behavior CollideD3Q19 for FreeSurfaceWater {
    update when self.kind > 0.0 {
        let omega = 1.0 / self.tau;
        let cu0 = 0.0;
        let equilibrium0 = 0.3333333333333333 * self.rho * (1.0 + 3.0 * (0.0) + 4.5 * (0.0) * (0.0) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force0 = 0.3333333333333333 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (0.0 - self.uz) + 9.0 * cu0 * 0.0);
        let cu1 = (-self.ux) + (-self.uy);
        let equilibrium1 = 0.027777777777777776 * self.rho * (1.0 + 3.0 * ((-self.ux) + (-self.uy)) + 4.5 * ((-self.ux) + (-self.uy)) * ((-self.ux) + (-self.uy)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force1 = 0.027777777777777776 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (0.0 - self.uz) + 9.0 * cu1 * 0.0);
        let cu2 = (-self.ux) + (-self.uz);
        let equilibrium2 = 0.027777777777777776 * self.rho * (1.0 + 3.0 * ((-self.ux) + (-self.uz)) + 4.5 * ((-self.ux) + (-self.uz)) * ((-self.ux) + (-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force2 = 0.027777777777777776 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (-1.0 - self.uz) + 9.0 * cu2 * -1.0);
        let cu3 = (-self.ux);
        let equilibrium3 = 0.05555555555555555 * self.rho * (1.0 + 3.0 * ((-self.ux)) + 4.5 * ((-self.ux)) * ((-self.ux)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force3 = 0.05555555555555555 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (0.0 - self.uz) + 9.0 * cu3 * 0.0);
        let cu4 = (-self.ux) + self.uz;
        let equilibrium4 = 0.027777777777777776 * self.rho * (1.0 + 3.0 * ((-self.ux) + self.uz) + 4.5 * ((-self.ux) + self.uz) * ((-self.ux) + self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force4 = 0.027777777777777776 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (1.0 - self.uz) + 9.0 * cu4 * 1.0);
        let cu5 = (-self.ux) + self.uy;
        let equilibrium5 = 0.027777777777777776 * self.rho * (1.0 + 3.0 * ((-self.ux) + self.uy) + 4.5 * ((-self.ux) + self.uy) * ((-self.ux) + self.uy) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force5 = 0.027777777777777776 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (0.0 - self.uz) + 9.0 * cu5 * 0.0);
        let cu6 = (-self.uy) + (-self.uz);
        let equilibrium6 = 0.027777777777777776 * self.rho * (1.0 + 3.0 * ((-self.uy) + (-self.uz)) + 4.5 * ((-self.uy) + (-self.uz)) * ((-self.uy) + (-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force6 = 0.027777777777777776 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (-1.0 - self.uz) + 9.0 * cu6 * -1.0);
        let cu7 = (-self.uy);
        let equilibrium7 = 0.05555555555555555 * self.rho * (1.0 + 3.0 * ((-self.uy)) + 4.5 * ((-self.uy)) * ((-self.uy)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force7 = 0.05555555555555555 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (0.0 - self.uz) + 9.0 * cu7 * 0.0);
        let cu8 = (-self.uy) + self.uz;
        let equilibrium8 = 0.027777777777777776 * self.rho * (1.0 + 3.0 * ((-self.uy) + self.uz) + 4.5 * ((-self.uy) + self.uz) * ((-self.uy) + self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force8 = 0.027777777777777776 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (1.0 - self.uz) + 9.0 * cu8 * 1.0);
        let cu9 = (-self.uz);
        let equilibrium9 = 0.05555555555555555 * self.rho * (1.0 + 3.0 * ((-self.uz)) + 4.5 * ((-self.uz)) * ((-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force9 = 0.05555555555555555 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (-1.0 - self.uz) + 9.0 * cu9 * -1.0);
        let cu10 = self.uz;
        let equilibrium10 = 0.05555555555555555 * self.rho * (1.0 + 3.0 * (self.uz) + 4.5 * (self.uz) * (self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force10 = 0.05555555555555555 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (1.0 - self.uz) + 9.0 * cu10 * 1.0);
        let cu11 = self.uy + (-self.uz);
        let equilibrium11 = 0.027777777777777776 * self.rho * (1.0 + 3.0 * (self.uy + (-self.uz)) + 4.5 * (self.uy + (-self.uz)) * (self.uy + (-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force11 = 0.027777777777777776 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (-1.0 - self.uz) + 9.0 * cu11 * -1.0);
        let cu12 = self.uy;
        let equilibrium12 = 0.05555555555555555 * self.rho * (1.0 + 3.0 * (self.uy) + 4.5 * (self.uy) * (self.uy) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force12 = 0.05555555555555555 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (0.0 - self.uz) + 9.0 * cu12 * 0.0);
        let cu13 = self.uy + self.uz;
        let equilibrium13 = 0.027777777777777776 * self.rho * (1.0 + 3.0 * (self.uy + self.uz) + 4.5 * (self.uy + self.uz) * (self.uy + self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force13 = 0.027777777777777776 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (1.0 - self.uz) + 9.0 * cu13 * 1.0);
        let cu14 = self.ux + (-self.uy);
        let equilibrium14 = 0.027777777777777776 * self.rho * (1.0 + 3.0 * (self.ux + (-self.uy)) + 4.5 * (self.ux + (-self.uy)) * (self.ux + (-self.uy)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force14 = 0.027777777777777776 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (0.0 - self.uz) + 9.0 * cu14 * 0.0);
        let cu15 = self.ux + (-self.uz);
        let equilibrium15 = 0.027777777777777776 * self.rho * (1.0 + 3.0 * (self.ux + (-self.uz)) + 4.5 * (self.ux + (-self.uz)) * (self.ux + (-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force15 = 0.027777777777777776 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (-1.0 - self.uz) + 9.0 * cu15 * -1.0);
        let cu16 = self.ux;
        let equilibrium16 = 0.05555555555555555 * self.rho * (1.0 + 3.0 * (self.ux) + 4.5 * (self.ux) * (self.ux) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force16 = 0.05555555555555555 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (0.0 - self.uz) + 9.0 * cu16 * 0.0);
        let cu17 = self.ux + self.uz;
        let equilibrium17 = 0.027777777777777776 * self.rho * (1.0 + 3.0 * (self.ux + self.uz) + 4.5 * (self.ux + self.uz) * (self.ux + self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force17 = 0.027777777777777776 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (1.0 - self.uz) + 9.0 * cu17 * 1.0);
        let cu18 = self.ux + self.uy;
        let equilibrium18 = 0.027777777777777776 * self.rho * (1.0 + 3.0 * (self.ux + self.uy) + 4.5 * (self.ux + self.uy) * (self.ux + self.uy) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz));
        let force18 = 0.027777777777777776 * (1.0 - omega * 0.5) * self.rho * (-self.gravity) * (3.0 * (0.0 - self.uz) + 9.0 * cu18 * 0.0);
        set self.f0 = clamp(self.f0 + omega * (equilibrium0 - self.f0) + force0, -1000000.0, 1000000.0);
        set self.f1 = clamp(self.f1 + omega * (equilibrium1 - self.f1) + force1, -1000000.0, 1000000.0);
        set self.f2 = clamp(self.f2 + omega * (equilibrium2 - self.f2) + force2, -1000000.0, 1000000.0);
        set self.f3 = clamp(self.f3 + omega * (equilibrium3 - self.f3) + force3, -1000000.0, 1000000.0);
        set self.f4 = clamp(self.f4 + omega * (equilibrium4 - self.f4) + force4, -1000000.0, 1000000.0);
        set self.f5 = clamp(self.f5 + omega * (equilibrium5 - self.f5) + force5, -1000000.0, 1000000.0);
        set self.f6 = clamp(self.f6 + omega * (equilibrium6 - self.f6) + force6, -1000000.0, 1000000.0);
        set self.f7 = clamp(self.f7 + omega * (equilibrium7 - self.f7) + force7, -1000000.0, 1000000.0);
        set self.f8 = clamp(self.f8 + omega * (equilibrium8 - self.f8) + force8, -1000000.0, 1000000.0);
        set self.f9 = clamp(self.f9 + omega * (equilibrium9 - self.f9) + force9, -1000000.0, 1000000.0);
        set self.f10 = clamp(self.f10 + omega * (equilibrium10 - self.f10) + force10, -1000000.0, 1000000.0);
        set self.f11 = clamp(self.f11 + omega * (equilibrium11 - self.f11) + force11, -1000000.0, 1000000.0);
        set self.f12 = clamp(self.f12 + omega * (equilibrium12 - self.f12) + force12, -1000000.0, 1000000.0);
        set self.f13 = clamp(self.f13 + omega * (equilibrium13 - self.f13) + force13, -1000000.0, 1000000.0);
        set self.f14 = clamp(self.f14 + omega * (equilibrium14 - self.f14) + force14, -1000000.0, 1000000.0);
        set self.f15 = clamp(self.f15 + omega * (equilibrium15 - self.f15) + force15, -1000000.0, 1000000.0);
        set self.f16 = clamp(self.f16 + omega * (equilibrium16 - self.f16) + force16, -1000000.0, 1000000.0);
        set self.f17 = clamp(self.f17 + omega * (equilibrium17 - self.f17) + force17, -1000000.0, 1000000.0);
        set self.f18 = clamp(self.f18 + omega * (equilibrium18 - self.f18) + force18, -1000000.0, 1000000.0);
    }
}
behavior StreamBulkLiquid for FreeSurfaceWater
{
    update when self.kind == 2.0 {
        let incoming1 = (neighbors[FreeSurfaceWater using D3Q19, D18].f1 else self.f18);
        let incoming2 = (neighbors[FreeSurfaceWater using D3Q19, D17].f2 else self.f17);
        let incoming3 = (neighbors[FreeSurfaceWater using D3Q19, D16].f3 else self.f16);
        let incoming4 = (neighbors[FreeSurfaceWater using D3Q19, D15].f4 else self.f15);
        let incoming5 = (neighbors[FreeSurfaceWater using D3Q19, D14].f5 else self.f14);
        let incoming6 = (neighbors[FreeSurfaceWater using D3Q19, D13].f6 else self.f13);
        let incoming7 = (neighbors[FreeSurfaceWater using D3Q19, D12].f7 else self.f12);
        let incoming8 = (neighbors[FreeSurfaceWater using D3Q19, D11].f8 else self.f11);
        let incoming9 = (neighbors[FreeSurfaceWater using D3Q19, D10].f9 else self.f10);
        let incoming10 = (neighbors[FreeSurfaceWater using D3Q19, D9].f10 else self.f9);
        let incoming11 = (neighbors[FreeSurfaceWater using D3Q19, D8].f11 else self.f8);
        let incoming12 = (neighbors[FreeSurfaceWater using D3Q19, D7].f12 else self.f7);
        let incoming13 = (neighbors[FreeSurfaceWater using D3Q19, D6].f13 else self.f6);
        let incoming14 = (neighbors[FreeSurfaceWater using D3Q19, D5].f14 else self.f5);
        let incoming15 = (neighbors[FreeSurfaceWater using D3Q19, D4].f15 else self.f4);
        let incoming16 = (neighbors[FreeSurfaceWater using D3Q19, D3].f16 else self.f3);
        let incoming17 = (neighbors[FreeSurfaceWater using D3Q19, D2].f17 else self.f2);
        let incoming18 = (neighbors[FreeSurfaceWater using D3Q19, D1].f18 else self.f1);
        set self.f1 = clamp(incoming1, -1000000.0, 1000000.0);
        set self.f2 = clamp(incoming2, -1000000.0, 1000000.0);
        set self.f3 = clamp(incoming3, -1000000.0, 1000000.0);
        set self.f4 = clamp(incoming4, -1000000.0, 1000000.0);
        set self.f5 = clamp(incoming5, -1000000.0, 1000000.0);
        set self.f6 = clamp(incoming6, -1000000.0, 1000000.0);
        set self.f7 = clamp(incoming7, -1000000.0, 1000000.0);
        set self.f8 = clamp(incoming8, -1000000.0, 1000000.0);
        set self.f9 = clamp(incoming9, -1000000.0, 1000000.0);
        set self.f10 = clamp(incoming10, -1000000.0, 1000000.0);
        set self.f11 = clamp(incoming11, -1000000.0, 1000000.0);
        set self.f12 = clamp(incoming12, -1000000.0, 1000000.0);
        set self.f13 = clamp(incoming13, -1000000.0, 1000000.0);
        set self.f14 = clamp(incoming14, -1000000.0, 1000000.0);
        set self.f15 = clamp(incoming15, -1000000.0, 1000000.0);
        set self.f16 = clamp(incoming16, -1000000.0, 1000000.0);
        set self.f17 = clamp(incoming17, -1000000.0, 1000000.0);
        set self.f18 = clamp(incoming18, -1000000.0, 1000000.0);
        set self.mass = clamp(self.mass + (incoming1 - self.f18) + (incoming2 - self.f17) + (incoming3 - self.f16) + (incoming4 - self.f15) + (incoming5 - self.f14) + (incoming6 - self.f13) + (incoming7 - self.f12) + (incoming8 - self.f11) + (incoming9 - self.f10) + (incoming10 - self.f9) + (incoming11 - self.f8) + (incoming12 - self.f7) + (incoming13 - self.f6) + (incoming14 - self.f5) + (incoming15 - self.f4) + (incoming16 - self.f3) + (incoming17 - self.f2) + (incoming18 - self.f1), -1000000.0, 1000000.0);
    }
}
behavior StreamFreeSurface for FreeSurfaceWater
{
    update when self.kind == 1.0 {
        let type1 = (neighbors[FreeSurfaceWater using D3Q19, D18].kind else 3.0);
        let wall1 = max(type1 - 2.0, 0.0);
        let gas1 = max(1.0 - type1, 0.0);
        let wet1 = min(type1, 1.0) - wall1;
        let fluid1 = max(0.0, 1.0 - abs(type1 - 2.0));
        let interface1 = max(0.0, 1.0 - abs(type1 - 1.0));
        let incoming1 = (neighbors[FreeSurfaceWater using D3Q19, D18].f1 else 0.0);
        let reconstructed1 = (0.027777777777777776 * 1.0 * (1.0 + 3.0 * ((-self.ux) + (-self.uy)) + 4.5 * ((-self.ux) + (-self.uy)) * ((-self.ux) + (-self.uy)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.027777777777777776 * 1.0 * (1.0 + 3.0 * (self.ux + self.uy) + 4.5 * (self.ux + self.uy) * (self.ux + self.uy) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f18;
        let inward1 = (-self.gx) + (-self.gy);
        let reconstruct1 = wet1 * clamp(inward1 * 1000000000000.0, 0.0, 1.0);
        let interfacePair1 = self.interface * interface1;
        let incomingMask1 = 1.0 - interfacePair1 * clamp((neighbors[FreeSurfaceWater using D3Q19, D18].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask1 = 1.0 - interfacePair1 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D18].rank else 0.0), 0.0, 1.0);
        let area1 = wet1 * max(self.fluid, fluid1) + self.interface * interface1 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D18].fill else 0.0));
        let transfer1 = area1 * (incomingMask1 * incoming1 - outgoingMask1 * self.f18);
        let type2 = (neighbors[FreeSurfaceWater using D3Q19, D17].kind else 3.0);
        let wall2 = max(type2 - 2.0, 0.0);
        let gas2 = max(1.0 - type2, 0.0);
        let wet2 = min(type2, 1.0) - wall2;
        let fluid2 = max(0.0, 1.0 - abs(type2 - 2.0));
        let interface2 = max(0.0, 1.0 - abs(type2 - 1.0));
        let incoming2 = (neighbors[FreeSurfaceWater using D3Q19, D17].f2 else 0.0);
        let reconstructed2 = (0.027777777777777776 * 1.0 * (1.0 + 3.0 * ((-self.ux) + (-self.uz)) + 4.5 * ((-self.ux) + (-self.uz)) * ((-self.ux) + (-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.027777777777777776 * 1.0 * (1.0 + 3.0 * (self.ux + self.uz) + 4.5 * (self.ux + self.uz) * (self.ux + self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f17;
        let inward2 = (-self.gx) + (-self.gz);
        let reconstruct2 = wet2 * clamp(inward2 * 1000000000000.0, 0.0, 1.0);
        let interfacePair2 = self.interface * interface2;
        let incomingMask2 = 1.0 - interfacePair2 * clamp((neighbors[FreeSurfaceWater using D3Q19, D17].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask2 = 1.0 - interfacePair2 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D17].rank else 0.0), 0.0, 1.0);
        let area2 = wet2 * max(self.fluid, fluid2) + self.interface * interface2 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D17].fill else 0.0));
        let transfer2 = area2 * (incomingMask2 * incoming2 - outgoingMask2 * self.f17);
        let type3 = (neighbors[FreeSurfaceWater using D3Q19, D16].kind else 3.0);
        let wall3 = max(type3 - 2.0, 0.0);
        let gas3 = max(1.0 - type3, 0.0);
        let wet3 = min(type3, 1.0) - wall3;
        let fluid3 = max(0.0, 1.0 - abs(type3 - 2.0));
        let interface3 = max(0.0, 1.0 - abs(type3 - 1.0));
        let incoming3 = (neighbors[FreeSurfaceWater using D3Q19, D16].f3 else 0.0);
        let reconstructed3 = (0.05555555555555555 * 1.0 * (1.0 + 3.0 * ((-self.ux)) + 4.5 * ((-self.ux)) * ((-self.ux)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.05555555555555555 * 1.0 * (1.0 + 3.0 * (self.ux) + 4.5 * (self.ux) * (self.ux) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f16;
        let inward3 = (-self.gx);
        let reconstruct3 = wet3 * clamp(inward3 * 1000000000000.0, 0.0, 1.0);
        let interfacePair3 = self.interface * interface3;
        let incomingMask3 = 1.0 - interfacePair3 * clamp((neighbors[FreeSurfaceWater using D3Q19, D16].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask3 = 1.0 - interfacePair3 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D16].rank else 0.0), 0.0, 1.0);
        let area3 = wet3 * max(self.fluid, fluid3) + self.interface * interface3 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D16].fill else 0.0));
        let transfer3 = area3 * (incomingMask3 * incoming3 - outgoingMask3 * self.f16);
        let type4 = (neighbors[FreeSurfaceWater using D3Q19, D15].kind else 3.0);
        let wall4 = max(type4 - 2.0, 0.0);
        let gas4 = max(1.0 - type4, 0.0);
        let wet4 = min(type4, 1.0) - wall4;
        let fluid4 = max(0.0, 1.0 - abs(type4 - 2.0));
        let interface4 = max(0.0, 1.0 - abs(type4 - 1.0));
        let incoming4 = (neighbors[FreeSurfaceWater using D3Q19, D15].f4 else 0.0);
        let reconstructed4 = (0.027777777777777776 * 1.0 * (1.0 + 3.0 * ((-self.ux) + self.uz) + 4.5 * ((-self.ux) + self.uz) * ((-self.ux) + self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.027777777777777776 * 1.0 * (1.0 + 3.0 * (self.ux + (-self.uz)) + 4.5 * (self.ux + (-self.uz)) * (self.ux + (-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f15;
        let inward4 = (-self.gx) + self.gz;
        let reconstruct4 = wet4 * clamp(inward4 * 1000000000000.0, 0.0, 1.0);
        let interfacePair4 = self.interface * interface4;
        let incomingMask4 = 1.0 - interfacePair4 * clamp((neighbors[FreeSurfaceWater using D3Q19, D15].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask4 = 1.0 - interfacePair4 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D15].rank else 0.0), 0.0, 1.0);
        let area4 = wet4 * max(self.fluid, fluid4) + self.interface * interface4 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D15].fill else 0.0));
        let transfer4 = area4 * (incomingMask4 * incoming4 - outgoingMask4 * self.f15);
        let type5 = (neighbors[FreeSurfaceWater using D3Q19, D14].kind else 3.0);
        let wall5 = max(type5 - 2.0, 0.0);
        let gas5 = max(1.0 - type5, 0.0);
        let wet5 = min(type5, 1.0) - wall5;
        let fluid5 = max(0.0, 1.0 - abs(type5 - 2.0));
        let interface5 = max(0.0, 1.0 - abs(type5 - 1.0));
        let incoming5 = (neighbors[FreeSurfaceWater using D3Q19, D14].f5 else 0.0);
        let reconstructed5 = (0.027777777777777776 * 1.0 * (1.0 + 3.0 * ((-self.ux) + self.uy) + 4.5 * ((-self.ux) + self.uy) * ((-self.ux) + self.uy) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.027777777777777776 * 1.0 * (1.0 + 3.0 * (self.ux + (-self.uy)) + 4.5 * (self.ux + (-self.uy)) * (self.ux + (-self.uy)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f14;
        let inward5 = (-self.gx) + self.gy;
        let reconstruct5 = wet5 * clamp(inward5 * 1000000000000.0, 0.0, 1.0);
        let interfacePair5 = self.interface * interface5;
        let incomingMask5 = 1.0 - interfacePair5 * clamp((neighbors[FreeSurfaceWater using D3Q19, D14].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask5 = 1.0 - interfacePair5 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D14].rank else 0.0), 0.0, 1.0);
        let area5 = wet5 * max(self.fluid, fluid5) + self.interface * interface5 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D14].fill else 0.0));
        let transfer5 = area5 * (incomingMask5 * incoming5 - outgoingMask5 * self.f14);
        let type6 = (neighbors[FreeSurfaceWater using D3Q19, D13].kind else 3.0);
        let wall6 = max(type6 - 2.0, 0.0);
        let gas6 = max(1.0 - type6, 0.0);
        let wet6 = min(type6, 1.0) - wall6;
        let fluid6 = max(0.0, 1.0 - abs(type6 - 2.0));
        let interface6 = max(0.0, 1.0 - abs(type6 - 1.0));
        let incoming6 = (neighbors[FreeSurfaceWater using D3Q19, D13].f6 else 0.0);
        let reconstructed6 = (0.027777777777777776 * 1.0 * (1.0 + 3.0 * ((-self.uy) + (-self.uz)) + 4.5 * ((-self.uy) + (-self.uz)) * ((-self.uy) + (-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.027777777777777776 * 1.0 * (1.0 + 3.0 * (self.uy + self.uz) + 4.5 * (self.uy + self.uz) * (self.uy + self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f13;
        let inward6 = (-self.gy) + (-self.gz);
        let reconstruct6 = wet6 * clamp(inward6 * 1000000000000.0, 0.0, 1.0);
        let interfacePair6 = self.interface * interface6;
        let incomingMask6 = 1.0 - interfacePair6 * clamp((neighbors[FreeSurfaceWater using D3Q19, D13].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask6 = 1.0 - interfacePair6 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D13].rank else 0.0), 0.0, 1.0);
        let area6 = wet6 * max(self.fluid, fluid6) + self.interface * interface6 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D13].fill else 0.0));
        let transfer6 = area6 * (incomingMask6 * incoming6 - outgoingMask6 * self.f13);
        let type7 = (neighbors[FreeSurfaceWater using D3Q19, D12].kind else 3.0);
        let wall7 = max(type7 - 2.0, 0.0);
        let gas7 = max(1.0 - type7, 0.0);
        let wet7 = min(type7, 1.0) - wall7;
        let fluid7 = max(0.0, 1.0 - abs(type7 - 2.0));
        let interface7 = max(0.0, 1.0 - abs(type7 - 1.0));
        let incoming7 = (neighbors[FreeSurfaceWater using D3Q19, D12].f7 else 0.0);
        let reconstructed7 = (0.05555555555555555 * 1.0 * (1.0 + 3.0 * ((-self.uy)) + 4.5 * ((-self.uy)) * ((-self.uy)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.05555555555555555 * 1.0 * (1.0 + 3.0 * (self.uy) + 4.5 * (self.uy) * (self.uy) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f12;
        let inward7 = (-self.gy);
        let reconstruct7 = wet7 * clamp(inward7 * 1000000000000.0, 0.0, 1.0);
        let interfacePair7 = self.interface * interface7;
        let incomingMask7 = 1.0 - interfacePair7 * clamp((neighbors[FreeSurfaceWater using D3Q19, D12].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask7 = 1.0 - interfacePair7 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D12].rank else 0.0), 0.0, 1.0);
        let area7 = wet7 * max(self.fluid, fluid7) + self.interface * interface7 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D12].fill else 0.0));
        let transfer7 = area7 * (incomingMask7 * incoming7 - outgoingMask7 * self.f12);
        let type8 = (neighbors[FreeSurfaceWater using D3Q19, D11].kind else 3.0);
        let wall8 = max(type8 - 2.0, 0.0);
        let gas8 = max(1.0 - type8, 0.0);
        let wet8 = min(type8, 1.0) - wall8;
        let fluid8 = max(0.0, 1.0 - abs(type8 - 2.0));
        let interface8 = max(0.0, 1.0 - abs(type8 - 1.0));
        let incoming8 = (neighbors[FreeSurfaceWater using D3Q19, D11].f8 else 0.0);
        let reconstructed8 = (0.027777777777777776 * 1.0 * (1.0 + 3.0 * ((-self.uy) + self.uz) + 4.5 * ((-self.uy) + self.uz) * ((-self.uy) + self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.027777777777777776 * 1.0 * (1.0 + 3.0 * (self.uy + (-self.uz)) + 4.5 * (self.uy + (-self.uz)) * (self.uy + (-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f11;
        let inward8 = (-self.gy) + self.gz;
        let reconstruct8 = wet8 * clamp(inward8 * 1000000000000.0, 0.0, 1.0);
        let interfacePair8 = self.interface * interface8;
        let incomingMask8 = 1.0 - interfacePair8 * clamp((neighbors[FreeSurfaceWater using D3Q19, D11].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask8 = 1.0 - interfacePair8 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D11].rank else 0.0), 0.0, 1.0);
        let area8 = wet8 * max(self.fluid, fluid8) + self.interface * interface8 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D11].fill else 0.0));
        let transfer8 = area8 * (incomingMask8 * incoming8 - outgoingMask8 * self.f11);
        let type9 = (neighbors[FreeSurfaceWater using D3Q19, D10].kind else 3.0);
        let wall9 = max(type9 - 2.0, 0.0);
        let gas9 = max(1.0 - type9, 0.0);
        let wet9 = min(type9, 1.0) - wall9;
        let fluid9 = max(0.0, 1.0 - abs(type9 - 2.0));
        let interface9 = max(0.0, 1.0 - abs(type9 - 1.0));
        let incoming9 = (neighbors[FreeSurfaceWater using D3Q19, D10].f9 else 0.0);
        let reconstructed9 = (0.05555555555555555 * 1.0 * (1.0 + 3.0 * ((-self.uz)) + 4.5 * ((-self.uz)) * ((-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.05555555555555555 * 1.0 * (1.0 + 3.0 * (self.uz) + 4.5 * (self.uz) * (self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f10;
        let inward9 = (-self.gz);
        let reconstruct9 = wet9 * clamp(inward9 * 1000000000000.0, 0.0, 1.0);
        let interfacePair9 = self.interface * interface9;
        let incomingMask9 = 1.0 - interfacePair9 * clamp((neighbors[FreeSurfaceWater using D3Q19, D10].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask9 = 1.0 - interfacePair9 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D10].rank else 0.0), 0.0, 1.0);
        let area9 = wet9 * max(self.fluid, fluid9) + self.interface * interface9 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D10].fill else 0.0));
        let transfer9 = area9 * (incomingMask9 * incoming9 - outgoingMask9 * self.f10);
        let type10 = (neighbors[FreeSurfaceWater using D3Q19, D9].kind else 3.0);
        let wall10 = max(type10 - 2.0, 0.0);
        let gas10 = max(1.0 - type10, 0.0);
        let wet10 = min(type10, 1.0) - wall10;
        let fluid10 = max(0.0, 1.0 - abs(type10 - 2.0));
        let interface10 = max(0.0, 1.0 - abs(type10 - 1.0));
        let incoming10 = (neighbors[FreeSurfaceWater using D3Q19, D9].f10 else 0.0);
        let reconstructed10 = (0.05555555555555555 * 1.0 * (1.0 + 3.0 * (self.uz) + 4.5 * (self.uz) * (self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.05555555555555555 * 1.0 * (1.0 + 3.0 * ((-self.uz)) + 4.5 * ((-self.uz)) * ((-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f9;
        let inward10 = self.gz;
        let reconstruct10 = wet10 * clamp(inward10 * 1000000000000.0, 0.0, 1.0);
        let interfacePair10 = self.interface * interface10;
        let incomingMask10 = 1.0 - interfacePair10 * clamp((neighbors[FreeSurfaceWater using D3Q19, D9].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask10 = 1.0 - interfacePair10 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D9].rank else 0.0), 0.0, 1.0);
        let area10 = wet10 * max(self.fluid, fluid10) + self.interface * interface10 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D9].fill else 0.0));
        let transfer10 = area10 * (incomingMask10 * incoming10 - outgoingMask10 * self.f9);
        let type11 = (neighbors[FreeSurfaceWater using D3Q19, D8].kind else 3.0);
        let wall11 = max(type11 - 2.0, 0.0);
        let gas11 = max(1.0 - type11, 0.0);
        let wet11 = min(type11, 1.0) - wall11;
        let fluid11 = max(0.0, 1.0 - abs(type11 - 2.0));
        let interface11 = max(0.0, 1.0 - abs(type11 - 1.0));
        let incoming11 = (neighbors[FreeSurfaceWater using D3Q19, D8].f11 else 0.0);
        let reconstructed11 = (0.027777777777777776 * 1.0 * (1.0 + 3.0 * (self.uy + (-self.uz)) + 4.5 * (self.uy + (-self.uz)) * (self.uy + (-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.027777777777777776 * 1.0 * (1.0 + 3.0 * ((-self.uy) + self.uz) + 4.5 * ((-self.uy) + self.uz) * ((-self.uy) + self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f8;
        let inward11 = self.gy + (-self.gz);
        let reconstruct11 = wet11 * clamp(inward11 * 1000000000000.0, 0.0, 1.0);
        let interfacePair11 = self.interface * interface11;
        let incomingMask11 = 1.0 - interfacePair11 * clamp((neighbors[FreeSurfaceWater using D3Q19, D8].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask11 = 1.0 - interfacePair11 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D8].rank else 0.0), 0.0, 1.0);
        let area11 = wet11 * max(self.fluid, fluid11) + self.interface * interface11 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D8].fill else 0.0));
        let transfer11 = area11 * (incomingMask11 * incoming11 - outgoingMask11 * self.f8);
        let type12 = (neighbors[FreeSurfaceWater using D3Q19, D7].kind else 3.0);
        let wall12 = max(type12 - 2.0, 0.0);
        let gas12 = max(1.0 - type12, 0.0);
        let wet12 = min(type12, 1.0) - wall12;
        let fluid12 = max(0.0, 1.0 - abs(type12 - 2.0));
        let interface12 = max(0.0, 1.0 - abs(type12 - 1.0));
        let incoming12 = (neighbors[FreeSurfaceWater using D3Q19, D7].f12 else 0.0);
        let reconstructed12 = (0.05555555555555555 * 1.0 * (1.0 + 3.0 * (self.uy) + 4.5 * (self.uy) * (self.uy) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.05555555555555555 * 1.0 * (1.0 + 3.0 * ((-self.uy)) + 4.5 * ((-self.uy)) * ((-self.uy)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f7;
        let inward12 = self.gy;
        let reconstruct12 = wet12 * clamp(inward12 * 1000000000000.0, 0.0, 1.0);
        let interfacePair12 = self.interface * interface12;
        let incomingMask12 = 1.0 - interfacePair12 * clamp((neighbors[FreeSurfaceWater using D3Q19, D7].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask12 = 1.0 - interfacePair12 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D7].rank else 0.0), 0.0, 1.0);
        let area12 = wet12 * max(self.fluid, fluid12) + self.interface * interface12 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D7].fill else 0.0));
        let transfer12 = area12 * (incomingMask12 * incoming12 - outgoingMask12 * self.f7);
        let type13 = (neighbors[FreeSurfaceWater using D3Q19, D6].kind else 3.0);
        let wall13 = max(type13 - 2.0, 0.0);
        let gas13 = max(1.0 - type13, 0.0);
        let wet13 = min(type13, 1.0) - wall13;
        let fluid13 = max(0.0, 1.0 - abs(type13 - 2.0));
        let interface13 = max(0.0, 1.0 - abs(type13 - 1.0));
        let incoming13 = (neighbors[FreeSurfaceWater using D3Q19, D6].f13 else 0.0);
        let reconstructed13 = (0.027777777777777776 * 1.0 * (1.0 + 3.0 * (self.uy + self.uz) + 4.5 * (self.uy + self.uz) * (self.uy + self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.027777777777777776 * 1.0 * (1.0 + 3.0 * ((-self.uy) + (-self.uz)) + 4.5 * ((-self.uy) + (-self.uz)) * ((-self.uy) + (-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f6;
        let inward13 = self.gy + self.gz;
        let reconstruct13 = wet13 * clamp(inward13 * 1000000000000.0, 0.0, 1.0);
        let interfacePair13 = self.interface * interface13;
        let incomingMask13 = 1.0 - interfacePair13 * clamp((neighbors[FreeSurfaceWater using D3Q19, D6].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask13 = 1.0 - interfacePair13 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D6].rank else 0.0), 0.0, 1.0);
        let area13 = wet13 * max(self.fluid, fluid13) + self.interface * interface13 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D6].fill else 0.0));
        let transfer13 = area13 * (incomingMask13 * incoming13 - outgoingMask13 * self.f6);
        let type14 = (neighbors[FreeSurfaceWater using D3Q19, D5].kind else 3.0);
        let wall14 = max(type14 - 2.0, 0.0);
        let gas14 = max(1.0 - type14, 0.0);
        let wet14 = min(type14, 1.0) - wall14;
        let fluid14 = max(0.0, 1.0 - abs(type14 - 2.0));
        let interface14 = max(0.0, 1.0 - abs(type14 - 1.0));
        let incoming14 = (neighbors[FreeSurfaceWater using D3Q19, D5].f14 else 0.0);
        let reconstructed14 = (0.027777777777777776 * 1.0 * (1.0 + 3.0 * (self.ux + (-self.uy)) + 4.5 * (self.ux + (-self.uy)) * (self.ux + (-self.uy)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.027777777777777776 * 1.0 * (1.0 + 3.0 * ((-self.ux) + self.uy) + 4.5 * ((-self.ux) + self.uy) * ((-self.ux) + self.uy) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f5;
        let inward14 = self.gx + (-self.gy);
        let reconstruct14 = wet14 * clamp(inward14 * 1000000000000.0, 0.0, 1.0);
        let interfacePair14 = self.interface * interface14;
        let incomingMask14 = 1.0 - interfacePair14 * clamp((neighbors[FreeSurfaceWater using D3Q19, D5].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask14 = 1.0 - interfacePair14 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D5].rank else 0.0), 0.0, 1.0);
        let area14 = wet14 * max(self.fluid, fluid14) + self.interface * interface14 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D5].fill else 0.0));
        let transfer14 = area14 * (incomingMask14 * incoming14 - outgoingMask14 * self.f5);
        let type15 = (neighbors[FreeSurfaceWater using D3Q19, D4].kind else 3.0);
        let wall15 = max(type15 - 2.0, 0.0);
        let gas15 = max(1.0 - type15, 0.0);
        let wet15 = min(type15, 1.0) - wall15;
        let fluid15 = max(0.0, 1.0 - abs(type15 - 2.0));
        let interface15 = max(0.0, 1.0 - abs(type15 - 1.0));
        let incoming15 = (neighbors[FreeSurfaceWater using D3Q19, D4].f15 else 0.0);
        let reconstructed15 = (0.027777777777777776 * 1.0 * (1.0 + 3.0 * (self.ux + (-self.uz)) + 4.5 * (self.ux + (-self.uz)) * (self.ux + (-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.027777777777777776 * 1.0 * (1.0 + 3.0 * ((-self.ux) + self.uz) + 4.5 * ((-self.ux) + self.uz) * ((-self.ux) + self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f4;
        let inward15 = self.gx + (-self.gz);
        let reconstruct15 = wet15 * clamp(inward15 * 1000000000000.0, 0.0, 1.0);
        let interfacePair15 = self.interface * interface15;
        let incomingMask15 = 1.0 - interfacePair15 * clamp((neighbors[FreeSurfaceWater using D3Q19, D4].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask15 = 1.0 - interfacePair15 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D4].rank else 0.0), 0.0, 1.0);
        let area15 = wet15 * max(self.fluid, fluid15) + self.interface * interface15 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D4].fill else 0.0));
        let transfer15 = area15 * (incomingMask15 * incoming15 - outgoingMask15 * self.f4);
        let type16 = (neighbors[FreeSurfaceWater using D3Q19, D3].kind else 3.0);
        let wall16 = max(type16 - 2.0, 0.0);
        let gas16 = max(1.0 - type16, 0.0);
        let wet16 = min(type16, 1.0) - wall16;
        let fluid16 = max(0.0, 1.0 - abs(type16 - 2.0));
        let interface16 = max(0.0, 1.0 - abs(type16 - 1.0));
        let incoming16 = (neighbors[FreeSurfaceWater using D3Q19, D3].f16 else 0.0);
        let reconstructed16 = (0.05555555555555555 * 1.0 * (1.0 + 3.0 * (self.ux) + 4.5 * (self.ux) * (self.ux) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.05555555555555555 * 1.0 * (1.0 + 3.0 * ((-self.ux)) + 4.5 * ((-self.ux)) * ((-self.ux)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f3;
        let inward16 = self.gx;
        let reconstruct16 = wet16 * clamp(inward16 * 1000000000000.0, 0.0, 1.0);
        let interfacePair16 = self.interface * interface16;
        let incomingMask16 = 1.0 - interfacePair16 * clamp((neighbors[FreeSurfaceWater using D3Q19, D3].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask16 = 1.0 - interfacePair16 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D3].rank else 0.0), 0.0, 1.0);
        let area16 = wet16 * max(self.fluid, fluid16) + self.interface * interface16 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D3].fill else 0.0));
        let transfer16 = area16 * (incomingMask16 * incoming16 - outgoingMask16 * self.f3);
        let type17 = (neighbors[FreeSurfaceWater using D3Q19, D2].kind else 3.0);
        let wall17 = max(type17 - 2.0, 0.0);
        let gas17 = max(1.0 - type17, 0.0);
        let wet17 = min(type17, 1.0) - wall17;
        let fluid17 = max(0.0, 1.0 - abs(type17 - 2.0));
        let interface17 = max(0.0, 1.0 - abs(type17 - 1.0));
        let incoming17 = (neighbors[FreeSurfaceWater using D3Q19, D2].f17 else 0.0);
        let reconstructed17 = (0.027777777777777776 * 1.0 * (1.0 + 3.0 * (self.ux + self.uz) + 4.5 * (self.ux + self.uz) * (self.ux + self.uz) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.027777777777777776 * 1.0 * (1.0 + 3.0 * ((-self.ux) + (-self.uz)) + 4.5 * ((-self.ux) + (-self.uz)) * ((-self.ux) + (-self.uz)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f2;
        let inward17 = self.gx + self.gz;
        let reconstruct17 = wet17 * clamp(inward17 * 1000000000000.0, 0.0, 1.0);
        let interfacePair17 = self.interface * interface17;
        let incomingMask17 = 1.0 - interfacePair17 * clamp((neighbors[FreeSurfaceWater using D3Q19, D2].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask17 = 1.0 - interfacePair17 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D2].rank else 0.0), 0.0, 1.0);
        let area17 = wet17 * max(self.fluid, fluid17) + self.interface * interface17 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D2].fill else 0.0));
        let transfer17 = area17 * (incomingMask17 * incoming17 - outgoingMask17 * self.f2);
        let type18 = (neighbors[FreeSurfaceWater using D3Q19, D1].kind else 3.0);
        let wall18 = max(type18 - 2.0, 0.0);
        let gas18 = max(1.0 - type18, 0.0);
        let wet18 = min(type18, 1.0) - wall18;
        let fluid18 = max(0.0, 1.0 - abs(type18 - 2.0));
        let interface18 = max(0.0, 1.0 - abs(type18 - 1.0));
        let incoming18 = (neighbors[FreeSurfaceWater using D3Q19, D1].f18 else 0.0);
        let reconstructed18 = (0.027777777777777776 * 1.0 * (1.0 + 3.0 * (self.ux + self.uy) + 4.5 * (self.ux + self.uy) * (self.ux + self.uy) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) + (0.027777777777777776 * 1.0 * (1.0 + 3.0 * ((-self.ux) + (-self.uy)) + 4.5 * ((-self.ux) + (-self.uy)) * ((-self.ux) + (-self.uy)) - 1.5 * (self.ux * self.ux + self.uy * self.uy + self.uz * self.uz))) - self.f1;
        let inward18 = self.gx + self.gy;
        let reconstruct18 = wet18 * clamp(inward18 * 1000000000000.0, 0.0, 1.0);
        let interfacePair18 = self.interface * interface18;
        let incomingMask18 = 1.0 - interfacePair18 * clamp((neighbors[FreeSurfaceWater using D3Q19, D1].rank else 0.0) - self.rank, 0.0, 1.0);
        let outgoingMask18 = 1.0 - interfacePair18 * clamp(self.rank - (neighbors[FreeSurfaceWater using D3Q19, D1].rank else 0.0), 0.0, 1.0);
        let area18 = wet18 * max(self.fluid, fluid18) + self.interface * interface18 * 0.5 * (self.fill + (neighbors[FreeSurfaceWater using D3Q19, D1].fill else 0.0));
        let transfer18 = area18 * (incomingMask18 * incoming18 - outgoingMask18 * self.f1);
        set self.f1 = clamp((wet1 - reconstruct1) * incoming1 + (gas1 + reconstruct1) * reconstructed1 + wall1 * self.f18, -1000000.0, 1000000.0);
        set self.f2 = clamp((wet2 - reconstruct2) * incoming2 + (gas2 + reconstruct2) * reconstructed2 + wall2 * self.f17, -1000000.0, 1000000.0);
        set self.f3 = clamp((wet3 - reconstruct3) * incoming3 + (gas3 + reconstruct3) * reconstructed3 + wall3 * self.f16, -1000000.0, 1000000.0);
        set self.f4 = clamp((wet4 - reconstruct4) * incoming4 + (gas4 + reconstruct4) * reconstructed4 + wall4 * self.f15, -1000000.0, 1000000.0);
        set self.f5 = clamp((wet5 - reconstruct5) * incoming5 + (gas5 + reconstruct5) * reconstructed5 + wall5 * self.f14, -1000000.0, 1000000.0);
        set self.f6 = clamp((wet6 - reconstruct6) * incoming6 + (gas6 + reconstruct6) * reconstructed6 + wall6 * self.f13, -1000000.0, 1000000.0);
        set self.f7 = clamp((wet7 - reconstruct7) * incoming7 + (gas7 + reconstruct7) * reconstructed7 + wall7 * self.f12, -1000000.0, 1000000.0);
        set self.f8 = clamp((wet8 - reconstruct8) * incoming8 + (gas8 + reconstruct8) * reconstructed8 + wall8 * self.f11, -1000000.0, 1000000.0);
        set self.f9 = clamp((wet9 - reconstruct9) * incoming9 + (gas9 + reconstruct9) * reconstructed9 + wall9 * self.f10, -1000000.0, 1000000.0);
        set self.f10 = clamp((wet10 - reconstruct10) * incoming10 + (gas10 + reconstruct10) * reconstructed10 + wall10 * self.f9, -1000000.0, 1000000.0);
        set self.f11 = clamp((wet11 - reconstruct11) * incoming11 + (gas11 + reconstruct11) * reconstructed11 + wall11 * self.f8, -1000000.0, 1000000.0);
        set self.f12 = clamp((wet12 - reconstruct12) * incoming12 + (gas12 + reconstruct12) * reconstructed12 + wall12 * self.f7, -1000000.0, 1000000.0);
        set self.f13 = clamp((wet13 - reconstruct13) * incoming13 + (gas13 + reconstruct13) * reconstructed13 + wall13 * self.f6, -1000000.0, 1000000.0);
        set self.f14 = clamp((wet14 - reconstruct14) * incoming14 + (gas14 + reconstruct14) * reconstructed14 + wall14 * self.f5, -1000000.0, 1000000.0);
        set self.f15 = clamp((wet15 - reconstruct15) * incoming15 + (gas15 + reconstruct15) * reconstructed15 + wall15 * self.f4, -1000000.0, 1000000.0);
        set self.f16 = clamp((wet16 - reconstruct16) * incoming16 + (gas16 + reconstruct16) * reconstructed16 + wall16 * self.f3, -1000000.0, 1000000.0);
        set self.f17 = clamp((wet17 - reconstruct17) * incoming17 + (gas17 + reconstruct17) * reconstructed17 + wall17 * self.f2, -1000000.0, 1000000.0);
        set self.f18 = clamp((wet18 - reconstruct18) * incoming18 + (gas18 + reconstruct18) * reconstructed18 + wall18 * self.f1, -1000000.0, 1000000.0);
        set self.mass = clamp(self.mass + transfer1 + transfer2 + transfer3 + transfer4 + transfer5 + transfer6 + transfer7 + transfer8 + transfer9 + transfer10 + transfer11 + transfer12 + transfer13 + transfer14 + transfer15 + transfer16 + transfer17 + transfer18, -1000000.0, 1000000.0);
    }
}
// Repair fluid/gas contact after painting, and bootstrap the initial surface.
behavior PrepareSurface for FreeSurfaceWater {
    update {
        set self.wet = min(self.kind, 1.0);
        set self.fluid = max(self.kind - 1.0, 0.0);
        set self.interface = 1.0 - abs(self.kind - 1.0);
        set self.gas = max(1.0 - self.kind, 0.0);
        set self.fill = clamp(self.mass / max(self.rho, 0.000001), 0.0, 1.0);
    }
}

// Gradient points into the liquid. Rank implements the symmetric interface
// cleanup fluxes in Thuerey's table 4.1: isolated, regular, enclosed.
behavior SurfaceTopology for FreeSurfaceWater
{
    update {
        let hasFluid = min(sum n in neighbors[FreeSurfaceWater using D3Q19] { n.fluid }, 1.0);
        let hasGas = min(sum n in neighbors[FreeSurfaceWater using D3Q19] { n.gas }, 1.0);
        set self.rank = clamp(hasFluid * (2.0 - hasGas), 0.0, 2.0);
        set self.gx = clamp(0.5 * ((neighbors[FreeSurfaceWater using D3Q19, +X].fill else self.fill) - (neighbors[FreeSurfaceWater using D3Q19, -X].fill else self.fill)), -1000000.0, 1000000.0);
        set self.gy = clamp(0.5 * ((neighbors[FreeSurfaceWater using D3Q19, +Y].fill else self.fill) - (neighbors[FreeSurfaceWater using D3Q19, -Y].fill else self.fill)), -1000000.0, 1000000.0);
        set self.gz = clamp(0.5 * ((neighbors[FreeSurfaceWater using D3Q19, +Z].fill else self.fill) - (neighbors[FreeSurfaceWater using D3Q19, -Z].fill else self.fill)), -1000000.0, 1000000.0);
    }
}

behavior RepairSurface for FreeSurfaceWater
{
    update when self.kind == 2.0 && (sum n in neighbors[FreeSurfaceWater using D3Q19] { n.gas }) > 0.0 {
        set self.kind = 1.0;
    }
}

// Proposals are decided from a snapshot before repairing the surrounding band.
behavior ProposeSurface for FreeSurfaceWater {
    update {
        set self.nextKind = self.kind;
        set self.filling = 0.0;
        set self.emptying = 0.0;
    }
}

behavior MarkFilled for FreeSurfaceWater
{
    update when self.kind == 1.0 && (self.mass + self.reserve > 1.001 * self.rho || ((sum n in neighbors[FreeSurfaceWater using D3Q19] { n.gas }) == 0.0 && self.mass > 0.9 * self.rho)) {
        set self.nextKind = 2.0;
        set self.filling = 1.0;
    }
}

behavior MarkEmptied for FreeSurfaceWater
{
    update when self.kind == 1.0 && (self.mass + self.reserve < -0.001 * self.rho || ((sum n in neighbors[FreeSurfaceWater using D3Q19] { n.fluid }) == 0.0 && self.mass < 0.1 * self.rho)) {
        set self.nextKind = 0.0;
        set self.emptying = 1.0;
    }
}

behavior ProtectFilledNeighbors for FreeSurfaceWater
{
    update when (self.kind == 0.0 || self.nextKind == 0.0) && (sum n in neighbors[FreeSurfaceWater using D3Q19] { n.filling }) > 0.0 {
        set self.nextKind = 1.0;
        set self.emptying = 0.0;
    }
}

behavior ExposeFluidNeighbors for FreeSurfaceWater
{
    update when self.kind == 2.0 && (sum n in neighbors[FreeSurfaceWater using D3Q19] { n.emptying }) > 0.0 {
        set self.nextKind = 1.0;
    }
}

behavior MarkRecipients for FreeSurfaceWater {
    update { set self.recipient = 1.0 - abs(self.nextKind - 1.0); }
}

// Keep an interface overshoot local until conversion; otherwise small fluxes
// never cross the transition threshold and the interface becomes pinned.
// Never discard conversion excess. A signed reserve retains it if no interface
// neighbor exists; physical mass plus reserve is the conserved quantity.
behavior PlanRedistribution for FreeSurfaceWater
{
    update {
        let target = max(self.nextKind - 1.0, 0.0) * self.rho + self.recipient * (self.mass + self.reserve);
        let excess = self.mass - target + self.reserve;
        let count = sum n in neighbors[FreeSurfaceWater using D3Q19] { n.recipient };
        let hasRecipient = min(count, 1.0);
        set self.targetMass = clamp(target, -1000000.0, 1000000.0);
        set self.share = clamp(hasRecipient * excess / max(count, 1.0), -1000000.0, 1000000.0);
        set self.carried = clamp((1.0 - hasRecipient) * excess, -1000000.0, 1000000.0);
    }
}

// Newly wet cells start at the neighbors' mean density/velocity, without adding
// water. Initializing distributions creates pressure state, not physical mass.
behavior InitializeInterface for FreeSurfaceWater
{
    update when self.kind == 0.0 && self.nextKind == 1.0 {
        let count = max(sum n in neighbors[FreeSurfaceWater using D3Q19] { n.wet }, 1.0);
        let density = clamp((sum n in neighbors[FreeSurfaceWater using D3Q19] { n.rho * n.wet }) / count, 0.000001, 4.0);
        let u = (sum n in neighbors[FreeSurfaceWater using D3Q19] { n.ux * n.wet }) / count;
        let v = (sum n in neighbors[FreeSurfaceWater using D3Q19] { n.uy * n.wet }) / count;
        let w = (sum n in neighbors[FreeSurfaceWater using D3Q19] { n.uz * n.wet }) / count;
        set self.rho = density;
        set self.ux = clamp(u, -1000000.0, 1000000.0);
        set self.uy = clamp(v, -1000000.0, 1000000.0);
        set self.uz = clamp(w, -1000000.0, 1000000.0);
        set self.f0 = clamp(0.3333333333333333 * density * (1.0 + 3.0 * (0.0) + 4.5 * (0.0) * (0.0) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f1 = clamp(0.027777777777777776 * density * (1.0 + 3.0 * ((-u) + (-v)) + 4.5 * ((-u) + (-v)) * ((-u) + (-v)) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f2 = clamp(0.027777777777777776 * density * (1.0 + 3.0 * ((-u) + (-w)) + 4.5 * ((-u) + (-w)) * ((-u) + (-w)) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f3 = clamp(0.05555555555555555 * density * (1.0 + 3.0 * ((-u)) + 4.5 * ((-u)) * ((-u)) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f4 = clamp(0.027777777777777776 * density * (1.0 + 3.0 * ((-u) + w) + 4.5 * ((-u) + w) * ((-u) + w) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f5 = clamp(0.027777777777777776 * density * (1.0 + 3.0 * ((-u) + v) + 4.5 * ((-u) + v) * ((-u) + v) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f6 = clamp(0.027777777777777776 * density * (1.0 + 3.0 * ((-v) + (-w)) + 4.5 * ((-v) + (-w)) * ((-v) + (-w)) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f7 = clamp(0.05555555555555555 * density * (1.0 + 3.0 * ((-v)) + 4.5 * ((-v)) * ((-v)) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f8 = clamp(0.027777777777777776 * density * (1.0 + 3.0 * ((-v) + w) + 4.5 * ((-v) + w) * ((-v) + w) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f9 = clamp(0.05555555555555555 * density * (1.0 + 3.0 * ((-w)) + 4.5 * ((-w)) * ((-w)) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f10 = clamp(0.05555555555555555 * density * (1.0 + 3.0 * (w) + 4.5 * (w) * (w) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f11 = clamp(0.027777777777777776 * density * (1.0 + 3.0 * (v + (-w)) + 4.5 * (v + (-w)) * (v + (-w)) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f12 = clamp(0.05555555555555555 * density * (1.0 + 3.0 * (v) + 4.5 * (v) * (v) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f13 = clamp(0.027777777777777776 * density * (1.0 + 3.0 * (v + w) + 4.5 * (v + w) * (v + w) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f14 = clamp(0.027777777777777776 * density * (1.0 + 3.0 * (u + (-v)) + 4.5 * (u + (-v)) * (u + (-v)) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f15 = clamp(0.027777777777777776 * density * (1.0 + 3.0 * (u + (-w)) + 4.5 * (u + (-w)) * (u + (-w)) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f16 = clamp(0.05555555555555555 * density * (1.0 + 3.0 * (u) + 4.5 * (u) * (u) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f17 = clamp(0.027777777777777776 * density * (1.0 + 3.0 * (u + w) + 4.5 * (u + w) * (u + w) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
        set self.f18 = clamp(0.027777777777777776 * density * (1.0 + 3.0 * (u + v) + 4.5 * (u + v) * (u + v) - 1.5 * (u * u + v * v + w * w)), -1000000.0, 1000000.0);
    }
}

behavior CommitSurface for FreeSurfaceWater
{
    update {
        let incoming = sum n in neighbors[FreeSurfaceWater using D3Q19] { n.share };
        let raw = self.targetMass + self.recipient * incoming;
        let physicalMass = clamp(raw, 0.0, self.rho);
        set self.mass = physicalMass;
        set self.reserve = clamp(self.carried + raw - physicalMass, -1000000.0, 1000000.0);
        set self.kind = self.nextKind;
        set self.fill = clamp(physicalMass / max(self.rho, 0.000001), 0.0, 1.0);
    }
}

// Latch invalid operating conditions visibly; wide interval guards must never
// turn an unstable high-Mach run into an apparently successful simulation.
behavior CheckWaterRegime for FreeSurfaceWater {
    update when self.kind > 0.0 && (self.rho < 0.9 || self.rho > 1.1 || self.ux * self.ux + self.uy * self.uy + self.uz * self.uz > 0.04) {
        set self.unstable = 1.0;
    }
}

phase prepare {
    run PrepareSurface on LbmMedium;
    run PrepareSurface on Water;
}

phase repair {
    run RepairSurface on LbmMedium;
    run RepairSurface on Water;
}

phase classify {
    run PrepareSurface on LbmMedium;
    run PrepareSurface on Water;
}

phase moments {
    run RecoverMoments on LbmMedium;
    run RecoverMoments on Water;
}

phase topology {
    run SurfaceTopology on LbmMedium;
    run SurfaceTopology on Water;
}

phase collision {
    run CollideD3Q19 on LbmMedium;
    run CollideD3Q19 on Water;
}

phase stream {
    run StreamBulkLiquid on LbmMedium;
    run StreamBulkLiquid on Water;
    run StreamFreeSurface on LbmMedium;
    run StreamFreeSurface on Water;
}

phase recover {
    run RecoverMoments on LbmMedium;
    run RecoverMoments on Water;
}

phase proposal {
    run ProposeSurface on LbmMedium;
    run ProposeSurface on Water;
}

phase fill_flags {
    run MarkFilled on LbmMedium;
    run MarkFilled on Water;
}

phase empty_flags {
    run MarkEmptied on LbmMedium;
    run MarkEmptied on Water;
}

phase protect {
    run ProtectFilledNeighbors on LbmMedium;
    run ProtectFilledNeighbors on Water;
}

phase expose {
    run ExposeFluidNeighbors on LbmMedium;
    run ExposeFluidNeighbors on Water;
}

phase recipients {
    run MarkRecipients on LbmMedium;
    run MarkRecipients on Water;
}

phase redistribute {
    run PlanRedistribution on LbmMedium;
    run PlanRedistribution on Water;
}

phase initialize {
    run InitializeInterface on LbmMedium;
    run InitializeInterface on Water;
}

phase commit {
    run CommitSurface on LbmMedium;
    run CommitSurface on Water;
}

phase regime {
    run CheckWaterRegime on LbmMedium;
    run CheckWaterRegime on Water;
}
visual LbmWater for FreeSurfaceWater {
    color {
        let speed = clamp((abs(self.ux) + abs(self.uy) + abs(self.uz)) * 8.0, 0.0, 1.0);
        // A voxel view of fill, not a reconstructed subcell water surface.
        rgba(0.03 + speed * 0.3 + self.unstable * 0.7, 0.38 + speed * 0.25, 0.78, clamp(self.fill * 2.0, 0.0, 0.88))
    }
}
visual ClearTank for TankGlass { color { rgba(0.46, 0.79, 0.90, 0.055) } }
visual PorcelainTank for TankStone { color { rgba(0.48, 0.56, 0.62, 1.0) } }
visualize Glass with ClearTank;
visualize LbmMedium with LbmWater;
visualize Stone with PorcelainTank;
visualize Water with LbmWater;