Diffusion in a 3D grid

[isolin]

Hi, I am building a small ecosystem simulation. As part of it, I plan to have a 3D grid structure for the static environment with a diffusion process. The environment blocks will seek equilibrium through the diffusion, but other agent types will destabilize it through their actions. Now I am wondering what is the best approach to accomplish the diffusion inside of the agents simulation?

1. Border, edge and corner grid cells have a limited number of neighbors, they follow slightly different rules. Is it fine to have a single block type in the grid and include several ifs in the kernel to handle the border cases? Not sure about CUDA, but branching in shaders has been always slow. The whole top layer actually behaves a bit differently than the rest of cells below it. Or is it better to have a separate agent type for each (basically 27 types)? Or is there a better way like using states (I don't understand states in FlameGPU yet, improving the docs or adding an example would be great)?
2. I thought I would implement the diffusion through messaging. MessageSpatial3D will certainly do the job but it requires vector distance computation and iterates . Having a 3D grid of course calls for a "3D message box" like MessageArray3D but since the diffusion sends data to all neighbors, it would violate the "1 agent writes a message to each index" rule as I understand it. (Just to check, it actually means that "Each array element should receive messages from one agent only", right?). So I am wondering if I should stick with MessageSpatial3D or maybe split the diffusion process into several agent functions, one for each direction? In fact, after looking at MesageBucket is seems that having a dedicated bucket for each grid block where it will collect all messages from its neighbors should be the best solution. Well, the index space will be large (several millions of grid cells), each with two dozens of messages, but I hope it will work. Is may way of thinking correct, or are there perhaps better solutions?

I really like FlameGPU so far. Great job & thank you!

[mondus]

Hi @isolin . Thanks for using FLAME GPU and happy to answer your questions. It sounds like an interesting model.

1. You are right in that divergence is not great for GPUs. Divergence is much worse when it requires a data movement. FLAME GPU 2 tries to hides the impact of certain divergent behaviours such as boundary conditions through the message implementations. E.g. MessageArray3D is not wrapped so message reads from edge positions will return fewer messages (see next point about use of message gathers). This avoids the need for you to implement conditions. The implementation behind the scenes will of course have some compromises regarding divergence. The implementation inside the simulator is one which balances the pro and cons of different alternatives. We may of course improve this at some point if a more performant alternative is found but this would not effect your model implementation. One of the key principles of FLAME GPU 2 is to separate the simulator implementation (and GPU optimisations) from the model.
2. I would suggest to think of your problem as your agents gathering messages rather than scattering them to specific locations. This is very common in the FLAME GPU example models where a gather is used to get the relevant information to each agent. The choice of gather operation (i..e. using the most appropriate message type to reduce the number of messages) is important. Your agent function can also apply some filter to the results. E.g. This can be a distance check or some other filter. If you design your model this way then each agent can output a MessageArray3d message and then another agent function can gather these messages from all neighbours. This would be moist similar to the Game of Life or Sugarscape examples.

[isolin]

Hi, thank you for the answer! I looked at the Game of Life example (somehow yet I only looked at the docs and ignored the examples folder) and you are right that in that case gathering works perfectly. The reason is that the binary decision of life or death requires data from the surrounding neighborhood, but in that time step it only has an effect on the processed cell. In a diffusion process a physical quantity moves from one cell to another, so you need to update both the provider and the receiver in the same time step. That makes a stable system with a constant amount of the resources (the game of life is stable only in "special" cases). Sorry for the long explanation, you probably know all this, but maybe if someone else reads this, it may be helpful :)

The Sugarscape belongs to the same problem class, as sugar rates change independently of the neighborhood and movement is just a read operation over the set of sighted cell.

Thinking further, I realized that sine diffusion is a bilateral process, the naive solution actually computes the same twice. So for now I will try two solutions paths and report on their performance in a few days:

1. using MessageArray3D in a gathering way to make advantage of the double computation and reduce the number of messages, as every cell can compute how much it gives or looses to each neighbor (a simple linear blending), so computing the sum is straight forward
2. using MessageBucket in a spread way to halve the number of computations, sa each cell considers only neighbors with higher IDs and a simple back propagation to the rest of the neighbors happens in a second step as they read their bucket

I will let you know about my findings once I get back to it, perhaps next weekend.

[Robadob]

Today we've merged a new example which demonstrates diffusion in a wrapped 2D grid. You can find it here.

These are some outputs of the example's visualisation over time. Note, the visible lines are a side-effect of the visualisation being created as a collection of stacked cubes, so the edges are sometimes visible, we don't yet have a distinct visualiser for 2D models.

We've also introduced the ability to perform Von Neumann neighbourhood iteration of MessageArray2D and MessageArray3D.
There are some details regarding usage of these in the user guide here near the bottom of the page. However, in the coming week we're hoping to make some structural changes to the user guide, so that link may become invalid.

[isolin]

Thank you very much for the example! Most important of all I learned about the 3D visualization from it. I fought a lot with my diffusion, since in opposite to your example with a sphere-like topology I have a bounding box as hard boundary. Without a proper handling of edge cells I was losing energy in the system. Moreover, my application case is precipitation entering soil, descending deeper due to gravity and later eventually evaporating and leaving the system again. So I needed to take case of vertical and horizontal diffusion separately. It's a very simplified model, but it should be fine for my next goal to grow virtual plants.

I found it useful to consider cells at z = 0 to cover everything above the surface but besides that I have just a regular 3D grid.

My final solution works in 4 steps and I found scattering easier to keep under control:

1. surface cells at z=0 receive precipitation; all cells except those in the lowest layer decide how much water will flow to the next deeper cell
2. all cells receive vertical diffusion from cells directly above them (this needs to be a separate step for the sake of sync)
3. all cells decide how much will water will flow horizontally to the direct neighbors in the same layer
4. all cells receive horizontal diffusion

Steps 3 and 4 go as follows:

FLAMEGPU_AGENT_FUNCTION(soil_hdiff, flamegpu::MessageArray3D, flamegpu::MessageArray3D) {
    //[left out] read agent variables ix, iy, iz (coords) and water (current water amount)

    const float waterSides[4] = {
       max(0.0f, (ix > 0 ? water - FLAMEGPU->message_in.at(ix - 1, iy, iz).getVariable<float>("water") : 0.0f)),
       max(0.0f, (ix < sizeX ? water - FLAMEGPU->message_in.at(ix + 1, iy, iz).getVariable<float>("water") : 0.0f)),
       max(0.0f, (iy > 0 ? water - FLAMEGPU->message_in.at(ix, iy - 1, iz).getVariable<float>("water") : 0.0f)),
       max(0.0f, (iy < sizeY ? water - FLAMEGPU->message_in.at(ix, iy + 1, iz).getVariable<float>("water") : 0.0f))
    };

    const float waterSidesSum = waterSides[0] + waterSides[1] + waterSides[2] + waterSides[3];

    if (waterSidesSum > 0.001f) //avoid micro-diffusion that would cause numerical instabilities
    {
        float sideFlow = water * sidewaysDiffusionCoef; //decide how much water can scatter to those neighbors with a lower water level
        FLAMEGPU->setVariable<float>("water", water - sideFlow);
        sideFlow /= waterSidesSum;
        for(int i = 0; i < 4; ++i) FLAMEGPU->message_out.setVariable<float, 4>("water_sides", i, waterSides[i] * sideFlow);
    }
    else //no diffusion if the difference is too low
    {
        FLAMEGPU->setVariable<float>("water", water);
        for(int i = 0; i < 4; ++i) FLAMEGPU->message_out.setVariable<float, 4>("water_sides", i, 0.0f);
    }

    FLAMEGPU->message_out.setIndex(ix, iy, iz);
    return flamegpu::ALIVE;
}

FLAMEGPU_AGENT_FUNCTION(soil_receives_diffusion, flamegpu::MessageArray3D, flamegpu::MessageArray3D) {
    //[left out] read agent variables ix, iy, iz (coords) and water (current water amount)

    //gathering the water "donated" by the neighbors
    if (ix > 0) water += FLAMEGPU->message_in.at(ix - 1, iy, iz).getVariable<float, 4>("water_sides", 1);
    if (ix < sizeX) water += FLAMEGPU->message_in.at(ix + 1, iy, iz).getVariable<float, 4>("water_sides", 0);
    if (iy > 0) water += FLAMEGPU->message_in.at(ix, iy - 1, iz).getVariable<float, 4>("water_sides", 3);
    if (iy < sizeY) water += FLAMEGPU->message_in.at(ix, iy + 1, iz).getVariable<float, 4>("water_sides", 2);

    FLAMEGPU->setVariable<float>("water", water);
    return flamegpu::ALIVE;
}