A better fluid model for Factorio.

fluid.zig

// 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: i32,
    potential_inflow: i32,
    potential_outflow: i32,
    
    pub inline fn capacity(node: Node) i32 {
        return node.width * node.height;
    }

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

    pub inline fn level(node: Node) i32 {
        return @divTrunc(node.contents, node.width);
    }
};

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: i32,
    node_a_min_level: i32,
    node_b_min_level: i32,
    // TODO: Allow top-up valves by introducing a max_level or min_space.
};

inline fn muldiv(x: i32, multiplier: i32, divisor: i32) i32 {
    return @intCast(@divTrunc(@as(i64, x) * @as(i64, multiplier), @as(i64, divisor)));
}

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

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

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

        // Compute ideal next flow level, based on level difference.
        for (model.edges) |*edge| {
            const effective_node_a_level = @max(0, edge.node_a.level() - edge.node_a_min_level);
            const effective_node_b_level = @max(0, edge.node_b.level() - edge.node_b_min_level);
            const force_ab_from_level = effective_node_a_level - effective_node_b_level;
            // Increase the flow by a portion of the level difference to model inertia.
            // Don't go higher than 1/2.
            edge.flow_ab += muldiv(force_ab_from_level, 1, 10);
            // If desired, at this point you can add resistance to edges, which
            // can reduce the maximum amount of flow based on the number edges
            // being traversed.
            // For example, to lower the speed exponentially with distance:
            // edge.flow_ab = muldiv(edge.flow_ab, 99, 100);
            // Or to lower it linearly with distance:
            // edge.flow_ab = @max(0, edge.flow_ab - 10);
            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.
        for (model.edges) |*edge| {
            if (edge.flow_ab > 0) {
                // Flow from A to B.
                const input_limited_flow = muldiv(edge.flow_ab, edge.node_a.contents, edge.node_a.potential_outflow);
                const output_limited_flow = muldiv(edge.flow_ab, 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 = muldiv(flow_ba, edge.node_b.contents, edge.node_b.potential_outflow);
                const output_limited_flow = muldiv(flow_ba, 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 -= edge.flow_ab;
            edge.node_b.contents += edge.flow_ab;
        }
    }
};

pub fn main() !void {
    var model = Model{
        .nodes = &.{},
        .edges = &.{},
    };
    model.tick();
}

Have you tried any kind of benchmarks of the code? In the FF thread you asserted it to be fast but how fast/slow would it be compared to a system that is coming in 2.0? My personal educated guess is this implementation here would be suffering heavily from having a pointer-based graph instead of having stuff aligned in continuous chunks in memory. Also having a lot of conditionals in the inner loop with all those min/max things is going to hurt pretty bad as well. That probably could be fixed but I highly doubt it'd be anywhere near as fast as the proposed implementation.

Unless I'm mistaken this one doesn't simulate individual pieces of pipe as the 1.1 algorithm does so it isn't comparable to that either but is a wholly different way of doing things that isn't really similar to anything in factorio. In terms of realism, it makes zero sense to have peak throughput at 50% capacity.

On the "doesn't depend on the build order of entities", how would the graph of the nodes be built at a pipe junction? In what order would branches be processed?

You're asking some great questions!

Have you tried any kind of benchmarks of the code? In the FF thread you
asserted it to be fast but how fast/slow would it be compared to a system that
is coming in 2.0?

I haven't, no, I've only thought about the complexity and what a maximum
throughput scenario would look like. Building realistic scenarios without a
graphical editor is kind of a pain.

My personal educated guess is this implementation here would be suffering
heavily from having a pointer-based graph instead of having stuff aligned in
continuous chunks in memory.

The fluid system itself is a graph, so this is unavoidable. My implementation
still loops over a list of all the edges consecutively, but I don't think it's
possible to avoid dereferencing completely.

Also having a lot of conditionals in the inner loop with all those min/max
things is going to hurt pretty bad as well.

I expect these mins and maxes to compile to conditional move instructions, which
do not suffer from the same performance cost due to mispredictions that
conditional jumps suffer from.

That probably could be fixed but I highly doubt it'd be anywhere near as fast
as the proposed implementation.

Unless I'm mistaken this one doesn't simulate individual pieces of pipe as the
1.1 algorithm does so it isn't comparable to that either but is a wholly
different way of doing things that isn't really similar to anything in
factorio.

The 2.0 implementation gains its speed from not simulating individual pieces,
but whole connected sets of them at the same time instead. Actually, in this
algorithm, the definition of nodes and edges is intentionally left a bit vague:

If you want a simulation of individual pieces like in 1.1, you can assign a Node
to each pipe or building and an Edge to each interface between them.

If you want a bit more efficiency, you can combine non-branching lines of pipes
into a single wider Node. You'd have much fewer Nodes and Edges to work with in
long pipelines, at the cost of a bit of realism. You can still apply flow
resistance relative to the length of the pipeline, but there would no longer be
a delay from water being available at the start of the pipeline to it being
available at the end.

You can also do what the 2.0 system does and assign a single Node to a much
larger network of pipes, including branches. This has the disadvantage that you
can no longer tell in which direction the fluid in a pipe is flowing and that it
might seem to teleport from start to end.

In terms of realism, it makes zero sense to have peak throughput at 50% capacity.

I guess you're right, it just makes the algorithm a lot easier. Although if you
consider it 100% full but at 50% of the maximum pressure rather than 50% full,
it makes a bit more sense. You don't want to push a fluid through a pipe at the
maximum supported pressure because if the exit is suddenly closed, the pipe will
go over pressure and explode. Half pressure pipes having the optimum flow allows
me to avoid simulating such over pressure scenarios.

Of course, in realistic fluids, pressure doesn't scale linearly with density,
but that level of realism might be too much to expect from a game.

Regarding realism, my goals were primarily to simulate momentum, a configurable
amount of resistance, and to make sure pipes can show in which direction the
fluid flows, which is a feature that seems to have been removed in the 2.0
model.

On the "doesn't depend on the build order of entities", how would the graph of
the nodes be built at a pipe junction? In what order would branches be
processed?

A junction can be represented as a single node with more than two edges to other
nodes. Each of the four loops in the algorithm avoids reading from any field
that gets set in the same loop, so no matter in which order the list of nodes or
edges is, the results still only depend on the results of the previous loop.
This also means you can arbitrarily parallelize each loop, as long as there are
sync points between the loops.

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.