A better fluid model for Factorio.

fluid.zig

// Run with `zig run -lc -O ReleaseSafe fluid.zig`.
// On a Ryzen 3700X, I'm getting about 0.45 million edges per tick on a single core.

const do_output = false;
const live_output = false;
const prevent_overflow = true;

const c = @cImport({
    @cInclude("stdio.h");
});

const std = @import("std");
const time = std.time;

inline fn itof(x: anytype) f32 {
    return @floatFromInt(x);
}

inline fn ftoi(x: f32) i64 {
    return @intFromFloat(@trunc(x));
}

// Optimal flow speed is when pipes are at half pressure.

const Node = struct {
    // All of these must be positive, I'm just using signed types to avoid casting.
    height: i32,
    width: i32,
    contents: i64,
    potential_inflow: f32,
    potential_outflow: f32,

    pub inline fn capacity(node: Node) i64 {
        return @as(i64, @intCast(node.width)) * @as(i64, @intCast(node.height));
    }

    pub inline fn space(node: Node) i64 {
        return node.capacity() - node.contents;
    }

    pub inline fn level(node: Node) f32 {
        return itof(node.contents) / itof(node.width);
    }

    pub inline fn empty(node: *Node) void {
        node.contents = 0;
    }

    pub inline fn fill(node: *Node) void {
        node.contents = node.capacity();
    }
};

const Edge = struct {
    node_a: *Node,
    node_b: *Node,
    // Positive means flow from A to B, negative means flow from B to A.
    flow_ab: f32,
};

// Invariant: The flow through an edge is never larger than the total amount of
// fluid in its input and never larger than the total amount of space in its
// output.

// The amount of flow we want to reach eventually after stabilizing per
// difference in level between consecutive pipes.
const flow_per_slope: f32 = 10;
const flow_decay: f32 = 1 / (flow_per_slope * 2);

// Expected stable flow for a pipeline of length l is:
//      stable flow = level difference * flow_per_slope / (length + 2 * flow_per_slope)

const Model = struct {
    nodes: []Node,
    edges: []Edge,

    fn tick(model: Model) void {
        // Clear flow potentials.
        if (prevent_overflow) {
            for (model.nodes) |*node| {
                node.potential_inflow = 0;
                node.potential_outflow = 0;
            }
        }

        for (model.edges) |*edge| {
            const force_ab_from_level = (edge.node_a.level() - edge.node_b.level()) * flow_per_slope;
            edge.flow_ab = (1 - flow_decay) * edge.flow_ab + flow_decay * force_ab_from_level;
            if (prevent_overflow) {
                edge.node_a.potential_inflow += @max(0, -edge.flow_ab);
                edge.node_a.potential_outflow += @max(0, edge.flow_ab);
                edge.node_b.potential_inflow += @max(0, edge.flow_ab);
                edge.node_b.potential_outflow += @max(0, -edge.flow_ab);
            }
        }

        // Scale down flow based on available input fluid and output space.
        if (prevent_overflow) {
            for (model.edges) |*edge| {
                if (edge.flow_ab > 0) {
                    // Flow from A to B.
                    const input_limited_flow = edge.flow_ab * itof(edge.node_a.contents) / edge.node_a.potential_outflow;
                    const output_limited_flow = edge.flow_ab * itof(edge.node_b.space()) / edge.node_b.potential_inflow;
                    edge.flow_ab = @min(edge.flow_ab, @min(input_limited_flow, output_limited_flow));
                } else if (edge.flow_ab < 0) {
                    // Flow from B to A.
                    const flow_ba = -edge.flow_ab;
                    const input_limited_flow = flow_ba * itof(edge.node_b.contents) / edge.node_b.potential_outflow;
                    const output_limited_flow = flow_ba * itof(edge.node_a.space()) / edge.node_a.potential_inflow;
                    edge.flow_ab = -@min(flow_ba, @min(input_limited_flow, output_limited_flow));
                }
            }
        }

        // Apply flow. This must be a separate step from computing the flow
        // since the flow computation is based on the node contents which are
        // changed here.
        for (model.edges) |edge| {
            edge.node_a.contents -= ftoi(edge.flow_ab);
            edge.node_b.contents += ftoi(edge.flow_ab);
        }
    }
};

pub fn main() !void {
    var gpa = std.heap.GeneralPurposeAllocator(.{}){};
    const heap = gpa.allocator();
    var prng = std.rand.DefaultPrng.init(0);
    const random = prng.random();

    const pipe = Node{
        .height = 999,
        .width = 1,
        .contents = 0,
        .potential_inflow = undefined,
        .potential_outflow = undefined,
    };

    const node_count = 0x100000;
    var nodes: []Node = try heap.create([node_count]Node);
    for (nodes) |*node| {
        node.* = pipe;
    }

    const edge_count = node_count - 1;
    var edges: []Edge = try heap.create([edge_count]Edge);
    for (0..edges.len) |i| {
        edges[i] = .{
            .node_a = &nodes[i],
            .node_b = &nodes[i+1],
            .flow_ab = 0,
        };
    }

    // Shuffle edges to get more realistic performance.
    random.shuffle(Edge, edges);

    var model = Model{
        .nodes = nodes,
        .edges = edges,
    };

    // nodes[0].width = 1000;
    // nodes[nodes.len - 1].width = 1000;
    // nodes[0].fill();
    // nodes[nodes.len - 1].empty();
    nodes[8].fill();

    var timer = try time.Timer.start();
    var start: u64 = undefined;
    const tick_ns: u64 = 16_666_666;
    const iter_count = 200;
    for (0..iter_count) |t| {
        start = timer.read();
        model.tick();
        if (do_output and !live_output) {
            if (t % 1 == 0) {
                for (0..nodes.len) |i| {
                    const node = nodes[i];
                    if (i % (edges.len >> 4) == 0) {
                        _ = c.printf("%3i ", ftoi(node.level()));
                    }
                }
                _ = c.printf("\n  ");
                for (0..edges.len) |i| {
                    const edge = edges[i];
                    if (i % (edges.len >> 4) == 0) {
                        _ = c.printf("%3i ", ftoi(edge.flow_ab));
                    }
                }
                _ = c.printf("\n\n");
            }
        } else if (do_output) {
            for (0..nodes.len) |i| {
                const node = nodes[i];
                if (i % (edges.len >> 4) == 0) {
                    _ = c.printf("%3i ", ftoi(node.level()));
                }
            }
            _ = c.printf("\r");
            _ = c.fflush(c.stdout);
            const elapsed = timer.read() - start;
            if (elapsed < tick_ns) {
                time.sleep(tick_ns - elapsed);
            }
        }
    }
    _ = c.printf("\n%i edges per tick\n", ftoi(itof(iter_count) * itof(edges.len) * tick_ns / itof(timer.read())));
}

The only actual if statement can also easily be converted into branchless code due to its symmetry by the way, I just wrote it like this because it's more readable.