Fast Gray-Scott reaction-diffusion C# component using Alea for GPU computation

[Bathsheba]

I wrote a component that computes reaction-diffusion like this, Reaction-Diffusion Tutorial, on a mesh. It builds on Vicente Soler and Laurent Delrieu’s work here reaction diffusion on triangular mesh - Grasshopper . Thanks guys!

I extended the inputs to take either a single value, or one for each mesh vertex, so f, k etc. can vary over the mesh just as the dir input for flow direction already does.

I used Alea 3.0.4, http://www.aleagpu.com/release/3_0_4/doc/, to do the computation on the GPU. It takes a while to initialize but then it goes fast.

Alea just needs the basic install, Community license is fine. I tried using Alea.Fody for automated memory management as they recommend, but it didn’t work so I just copied everything onto the GPU explicitly.

Installing Alea: download it using nuget or whatever they want. Put the DLLs in Rhino’s System directory, not in a subdirectory. Don’t forget to unblock them, don’t forget FSharp.Core, and put only the CUDA dll that matches your system, i.e. Windows 32, Windows 64 or Mac. Then you have to launch Rhino once as admin because Alea wants to create some folders in that directory. After that you can go back to launching normally. Lastly, to pull Alea into a component right-click it and do Manage Assemblies and point that at the dll’s, then also add “using” lines.

I could not properly install anything from Alea’s samples library, they’re all broken because they refer to some beta server that doesn’t exist anymore. But bits can be downloaded and there’s some source in each one, deeply nested in the tutorials folder, which is better than nothing if you want examples.

I didn’t do anything with this yet, but thought I’d post it while it’s fresh as an example of using Alea. I was surprised how easy it was to get it basically working. I don’t know much C# and don’t plan to ever write a lot of it, so the docs were hard for me to understand. It wasn’t fast until I figured out the memory management (profiling is your friend) but now it’s nice to see it run. 50K iterations? no problem.

reaction diffusion Alea.gh (2.0 MB)

[DanielPiker]

Very nice. Thanks Bathsheba!
I’ve been meaning to have a go with Alea for a while, and this should help make the jump.

[mohamad.elatab]

DLLs in Rhino’s System directory, not in a subdirectory

What you did looks really interesting !! i want to ask you if you can upload dll files i was trying to do what you said and didn’t manage to do it .
Thanks

[Bathsheba]

I think maybe not a great idea to put a copy of the Alea dll’s here, but if I could figure out how to install it maybe you can too.

On the first machine I tried it on, Windows PowerShell, “Install-Package Alea -Version 3.0.4” worked, it picked up Nuget automagically. But I don’t think that’s really supposed to work.

Let’s see…on this machine it didn’t work. I opened Visual Studio 2015 (free from Microsoft), opened a new project (random type), then did View/Other Windows/Package Manager Console. This brought up a command line window, and from there “Install-Package Alea -Version 3.0.4” worked. That is, it installed Alea and then told me where it put it. I think that’s how you’re supposed to do this.

NB I can’t provide actual support here. I don’t know beans about Windows or C# or Visual Studio, I just put on a blindfold and use the Force.

[Bathsheba]

When I said above that automated memory management with Alea.Fody [Gpu.Managed] didn’t work, I mean it didn’t improve performance in this case. It installed and compiled fine, there are probably lots of cases where it does work.

[mohamad.elatab]

Thanks Bathsheba i will try again following all the steps on a different machine!

[timothytai (tim)]

I believe you could simply extract/unzip the .nupkg file downloaded and the .dll’s would be in one of the folders?

[maxim2 (Maxim)]

Does anyone have a copy of the alea installer package saved somewhere? The project appears to be defunct and I can’t find it anywhere.

[Bathsheba]

Ah, that is sad. It looks like you can still get the package at
NuGet Gallery | Alea 3.0.4 .

[Bathsheba]

Best luck! Alea definitely still runs on Windows, it’s not quite dead code yet.

If it’s better to rebuild from scratch with a maintained code base, it can’t be that hard, and there are certainly better resources now than two years ago. Here’s an open-source implementation in Javascript Reaction-Diffusion Playground - Introduction , whatever this guy did has to be possible from Python.

[maxim2 (Maxim)]

Excellent! thanks for the links and suggestions. Very helpful indeed.
I’m pretty fresh to grasshopper and not versed in coding, so trying to figure out which existing tools I can best adapt for using reaction-diffusion algorithms.

[laurent_delrieu]

Hello
I am not sure the search for performance is the most important thing to begin with. I didn’t look closely at @Bathsheba code but if it is similar to mine it needs first a good triangular mesh with the correct size. It also depends on what you want, so many things needs to be tune. It is first lot of tests with “simple” geometry.
So the scripts on “old” links are still valid.

[Bathsheba]

GitHub - GollyGang/ready: A cross-platform implementation of various reaction-diffusion systems. is another implementation of a variety of systems on arbitrary meshes, using OpenCL for GPU support. I haven’t done more than look at it, and as Laurent says I bet the meshes have to be pretty good to make it work, but I think if I were rebuilding from scratch now, I might look hard at their code. Or if Alea doesn’t play along and you just need to quick spray some r-d on a mesh, maybe it’ll do that out of the box.

And as Laurent also says, if speed isn’t essential and you want to keep it in Grasshopper, his code without GPU support still runs just fine.

[maxim2 (Maxim)]

Laurent, you are right about speed not being crucial for the moment, just getting it to work reliably is a good start. So, far so good on that front.
However, for the purpose of fast iteration, the GPU acceleration is definitely going to be a bonus. At some point a rehash built on maintained code would be ideal of course.
I’ll let you know how I get along once I have had a chance to play with it properly.

[Bathsheba]

Last footnote here – the new grasshopper component Triremesh is a great match for all these! It produces really good-quality meshes that are perfect for this kind of thing.

[laurent_delrieu]

Or calculate something like Cotan Weight for each vertex of the mesh! I didn’t do it because some years ago Mesh Machine did also a pretty good job.

Compute Harmonic weights on a triangular mesh - Rodolphe Vaillant's homepage

Discrete Laplacian

I presume you are referring to the cotan weighting option for the Laplacian smoothing force in the old version of Kangaroo. (Kangaroo2 has instead the SoapFilm goal) MeshMachine also used cotan weighted smoothing for the relocation of the vertices. Here’s an example showing the effect of cotangent weighting on smoothing values per vertex (shown as colours). The mesh is deliberately uneven in triangle size. cotanweight_example.gh (28.0 KB) On the left is uniform weighted smoothing,…

Principal curvature approximation on meshes

Yes, using the cotan weights is a way of making the smoothing less dependent on the quality of the triangulation, compared to just taking the simple unweighted average. Remember though that here we can choose the points at which we take our initial curvature samples. The closer the distribution of these points is to an equilateral grid, the less difference there is between cotan weighted vs uniform weighted.

[DanielPiker]

Here’s a little method for getting the cotan weights at any vertex of a mesh from part of the script I posted in this thread. It takes a Plankton mesh but I think could be converted to use a normal Rhino mesh.
Though as you both point out, if the mesh quality is good, the difference between using this and just the uniform weighted Laplacian is very small, so it’s often not really needed.

  {
    int[] neighbours = P.Vertices.GetVertexNeighbours(i);
    Point3d vert = P.Vertices[i].ToPoint3d();
    int valence = P.Vertices.GetValence(i);
    Point3d[] neighbourPts = new Point3d[valence];
    Vector3d[] radial = new Vector3d[valence];
    Vector3d[] around = new Vector3d[valence];
    double[] cotWeight = new double[valence];
    double weightSum = 0;
    for (int j = 0; j < valence; j++)
    {
      neighbourPts[j] = P.Vertices[neighbours[j]].ToPoint3d();
      radial[j] = neighbourPts[j] - vert;
    }
    for (int j = 0; j < valence; j++)
    {
      around[j] = neighbourPts[(j + 1) % valence] - neighbourPts[j];
    }
    for (int j = 0; j < neighbours.Length; j++)
    {
      int previous = (j + valence - 1) % valence;
      Vector3d cross1 = Vector3d.CrossProduct(radial[previous], around[previous]);
      double dot1 = radial[previous] * around[previous];
      int next = (j + 1) % valence;
      Vector3d cross2 = Vector3d.CrossProduct(radial[next], around[j]);
      double dot2 = radial[next] * around[j];
      cotWeight[j] = Math.Abs(dot1 / cross1.Length) + Math.Abs(dot2 / cross2.Length);
      weightSum += cotWeight[j];
    }
    for (int j = 0; j < cotWeight.Length; j++) cotWeight[j] /= weightSum;
    return cotWeight;
  }

[Bathsheba]

This is a not-Rhino-related digression, but @laurent_delrieu maybe you’d like to know: I read about a neural net-based cellular automaton (https://distill.pub/2020/selforg/textures/) that synthesizes textures over a plane, and I used pretty much the same math to run it on a mesh. It’s in a colab notebook, Runtime/‘Run all’ if you wanna see it go.
@DanielPiker this called for Sobel operators as well as a fake Laplacian, and I did pretty much what you mention for that.