I like toying around with World Machine and other procedural terrain generation tools. I’m also trying to learn more about general-purpose GPU programming, so I thought I would see if I could recreate one of my favorite parts of World Machine on the GPU: hydraulic erosion.
I intend to replicate this basic functionality, but with an increase in speed. This is not meant to be boasting or impressive, just a consequence of using the GPU on a parallelizable problem. However I would like to improve upon it with a new feature: meanders and oxbow lakes. http://en.wikipedia.org/wiki/Meander Meanders are when rivers curve significantly, due to increased water flow and erosive force on the exterior of a bend in the river. Excessive curves may cause the river to erode an extra connection to itself, cutting off part of the river from a water source. This results in horse-shoe shaped lakes next to rivers (which may or may not dry up).
I would like to see if I can get this to happen as an emergent behavior from a properly constructed cellular automata. In order to see this sort of behavior I would need a model that allowed (1) water flowing next to a cell (which is higher than the water) to cause erosion in that cell, (2) water to flow faster on the outside of a bend, (3) faster moving water to be more erosive, (4) slower water to deposit more soil. I’m not entirely sure how to achieve that. Locally computing the gradient of water altitudes to determine flow rate each time step is the simplest method, but it also fails to achieve #2, which is an important property.
Does anyone have any ideas? I feel like I am too close to the problem and am overlooking something simple. How does World Machine do its erosion? (how much can you tell me, versus how much is the “secret sauce”?)
PS: I understand that what I am talking about is, strictly speaking, a finite element method instead of a proper CA. However, I have repeatedly seen this referred to as a CA elsewhere, so I do so here as well.
Also, to be clear I am asking for insight into erosion first, and cellular automata as a means to an end. I am not asking for help with the GPGPU model of programming, I think I have a good grasp on that.
It definitely sounds like an interesting attempt and I am very curious about what you can come up with!
In terms of how World Machine does its erosion, most of the “secret sauce” is actually in choice of parameters and deciding what to include and exclude – as a completely physically accurate model is going to be not only very difficult implement but also extremely slow to simulate. Of the four you mentioned World Machine implements all but #2 – which is also why you don’t see meanders forming in WM water courses. The 3 implemented properties are all that is really needed to implement “alpine” erosion (meanders being most important in the slow/lazy low gradient river areas). I would love to see a model that includes the second behavior as it would indeed most likely form natural meanders.
Thanks! I’ve been a big fan of World Machine. I never developed my creative skills and that is something I have always regretted/been self conscious of. However, World Machine lets me make amazing stuff regardless! I know it sounds really silly, but the sentiment is heartfelt.
I’ve found it difficult to simplify my assumptions of what the CA update rule should do, but I found that paper as a sign that I was heading in the correct direction. While their model is definitely involved, some aspects related to flow and erosion were simpler than mine.
I have the serial form of my algorithm partially done. It shows promise, and I’ve gotten some surprising results (that I will post later). Unfortunately I wont be able to work on it until after midterms.
I have lots of questions I would like to ask you while I have your attention, but I won’t be able to properly right now. So in “brief”: could we get a divergence free noise generator (or possibly a fault-algorithm generator, or a… I could go on), is GPU support through OpenCL support planned, would you add a node that lets us run arbitrary (image processing) convolution kernels, and are we ever going to be able to code our own nodes into WM? (macros don’t lend themselves well to many approaches, especially without an easy way to do iteration)
And if OpenCL support is not planned, would me submitting parallelized versions of extant WM functions help change your mind?
You can create your own nodes (devices) in World Machine by using the publicly-available PDK (http://www.world-machine.com/support.php?page=pdk), as long as you have MSVC 2010 to compile with. The free Express Edition of that compiler should work as well, although there are some annoying hoops to jump through to get 64bit support with it. You should be able to implement most anything you wish as the PDK exposes the same interface that WM devices themselves use.
1b) With regards to the lack of iteration; I agree! An “iteration macro” or “graph cycle” device that repeats the contents of a macro some N times (piping the outputs back to the inputs) has existed as a test device forever but being able to control it intuitively is a challenge; it is kept out of the main device set. But it would be extremely useful indeed, and I plan to include it when it can be exploited better…
A user-modifiable OpenCL device would be a very useful thing although obviously restricted to a more niche developer-type user. OpenCL support is not yet on the immediate schedule for WM, but it is something I am keenly interested in – this is an area I’ve been learning more about in preparations for being able to put it to real use.
There are lots more fundamental generators planned. Can’t give more details yet
I’d encourage you to check out the PDK and see if you can drop your algorithm into a WM device; it’s pretty straightforward and you will then be able to have it play with the entire host of abilities available in WM!
Sorry I did not get back to you in a timely manner. I have a bad habit of talking instead of doing, so I made a point to work on my project before I talked too much about it.
That said, I changed my project to investigate more general computational fluid dynamics techniques. By doing so I inadvertently gave myself the tools I need to complete my original project! I will post more later, but suffice to say that I am now rather confident that I know how to add meanders to WM.
csp256, this looks superb :shock: , especially the ability to generate such deltas and meanders. If this could be implemted somehow into WM, it would really raise the realism to new heights. Another approach to this whole river problem you might be interested in reading: http://hpcg.purdue.edu/?page=publication&id=170
Wow, that is amazing work! I knew I would regret not going to SIGGRAPH this year. Hopefully I will be able to go next year.
I have to finish writing a paper today, and either late tonight or sometime tomorrow I will expand upon what I have been doing.
There may be a problem with scaling the method up to the order of kilometers. The method is usually designed to operate on the order of millimeters or centimeters. We will see. We might be able to ‘fake it’. Thankfully we aren’t really going for accuracy, just what I call ‘pseudorealism’: aesthetically acceptable to a reasonably keen eye.
PS: Also note that I am trying to turn my work into an academic paper. If you stumble upon this, please play nice. If you have to use the material I have/will post, get in touch with me so you can cite me properly. It will be my first paper, so a lot is riding on it. Thank you.
I originally followed Wenda111287’s tutorial that involved World Machine. In doing so I accidentally created a large flat angled subsection of the terrain. When I tried to erode this, a straight channel was created. So I researched why, and started to think about ways I could fix this. Real erosion in those circumstances is chaotic: very minor perturbations of the state space cause significant long term differences to emerge.
I found this paper, and started to formulate ideas about how and why meanders form. It was crucial for the development of what I am about to detail. The only other resource on CA being used to describe erosion was an academic paper from an Italian research team, that had a very complicated model and a greater focus on physical accuracy.
It is my understanding now that WM has everything needed for meanders to form, except the ability for water to flow faster on the exterior of a bend [but no clear way of getting this effect; keep reading]. So I thought I would add momentum to the cells, and that this would fix it. I developed a system where the flow rate between cells was integrated over time, subject to some small exponential decay.
Each cell maintained a continuous, non negative quantity of water. Each link between von Neumann adjacent cells (roughly equal to twice the number of cells) contained a continuous flow rate. This flow rate was increment in proportion to the difference of the total altitudes (terrain height + water depth) of the cells the link joined, and then decreased by a small percentage to assure stability [this works like a low pass filter].
In the case that there was insufficient fluid for outflow, the cell is completely drained and each other cell receives an equivalent percentage of the fluid due to it. If there is no fluid on either cell of a link, then the link’s flow rate is set to 0. A border of cells 1 wide is created around the simulation, with their altitude and water depth are set independently by a first order linear interpolation of the exterior of the simulation proper. (the corner case is not needed due to the von Neumann neighborhood being use) This prevented pooling on the exterior (naïve implementations of many erosion algorithms encourage lakes to form on the border) and allowed unbiased outflow (that is to say, water that flowed to those cells were removed from the simulation).
I implemented this in Unity3d and it worked fairly okay. It showed basic wave mechanics (interference and diffraction) but it was fundamentally unstable and still expressed a tendency to prefer the interior bend. Flow on the exterior of a bend WAS slightly higher (don’t have a video showing this on me, sorry; it does have support for non level terrain), but it was typically much less than 1%. My test case for this was a 3, 5, and 7 wide right angle bend, with even terrain on the bend, and flat, sloped terrain (that is to say, planar) terrain leading into and out of the bend, with obstacles on the exterior.
I should see if I can get some videos of that, the ones with terrain are a bit more visually impressive.
Oscillations would build up if the damping coefficient was removed. I initially tried using a number (2% damping per time step) that I had empirically found performant when filtering the I term of a PID controller (for line following and other motion planning related problems) in a children’s robotics class I was teaching a couple years ago. Further experimentation showed the oscillatory behavior of the fluid simulation system was very similar to what I observed with the robot on a variety of damping coefficients, and 2% damping showed the best response while remaining stable in both systems. It is worth noting that the maximal (dimensionless) error in both cases was about +/– 100 [and when I altered the gain on the flow rate update, the optimal dampening coefficient changed in proportion]. I might investigate this later, but for now I will remember that there exists some rule of thumb for such systems.
I decided that the reason that the bend interior saw so much flow was because while there was momentum on the flow between cells (not that the fluid itself had inertia or momentum!) there was no mechanism to cause any test unit of volume to continue to propagate in a certain direction. Once it had flowed into a cell, it was not biased towards any of the three links out of a cell. I decided that I would need to create a way that each cell would contain a vector of four continuous quantities describing the amount of water to flow in each direction. These values could then be updated by the links, or some other rule could be implemented to cause the water to spread out (strict following of this rule would cause a unit of water to only travel on an axis!).
So I submitted my project proposal having no real idea how I was going to do what I wanted to, but a vague notion that I needed to keep track of portions of fluid moving in each direction, and some way to cause them to spill over into other bins. As I was taking an almost full load in this summer, I was very relieved to finish my project proposal in time. I treated myself with some undirected personal research (this being how I know about WM, why I was teaching children robotics, how I learned about CA and even how to program). I had come across a term, the Lattice Boltzmann Method (LBM), that I knew created very pretty pictures of turbulent fluids and used a discrete grid, but had no idea how it worked. Two hours later, I had discerned that I needed to change my project and abandoned direct work on the erosion simulator.
The LBM is very similar to my method, but is grounded in statistical mechanics (the branch of physics which uses the limiting, statistical properties of a very large number of particles undergoing interactions in a purely classical way to describe phenomenon which are impossible to approach directly; like fluid dynamics and thermodynamics). Boltzmann made huge contributions to this field and was a pretty cool guy (with an incredible beard). Historically the LBM evolved out of a simulation of discrete particles on a discrete grid with deterministic (or partially deterministic) update rules. There was a need to use statistics to average the discrete particles to discern the properties of the fluid at that position. The key insight in the LBM was to store the densities of particles which would move between cells and update those continuous values directly. This is directly analogous to my idea that I would need to store the continuous amounts of fluid which would transition in each direction, to conserve momentum globally.
I later stumbled upon that master’s thesis and it is what really made the LBM clear to me. Unfortunately I only found it after I had I started to implement the two dimensional (D2Q9) LBM in CUDA. I was directly porting the MATLAB code on this page:
I got about 343 million lattice updates per second, so a grid of 1000 by 1000 can be updated 343 times per second. The MATLAB code claims about 134 thousand lattice updates per second (about 1 million when compiled), while the optimized C++ implementation got around 3.15 million (on an old computer). There are stats elsewhere on their page that destroy my implementation in terms of performance, but use ~64 CPU cores to do so (while I am using a single Fermi-architecture M2070).
I tried to improve the performance using a tiled, multi-step per kernel implementation that allows tuning the arithmetic intensity by amortizing memory access times. Unfortunately I was dumb about the way I did it and mucked it up. I saw about 25% lower performance. I expect to see performance improve once I fix my mistakes. That is what I have planned for today.
That video in particular shows the increased flow around a bend. The plane waves you see coming off of the obstacles are because I initialized the fluid to moving right as if it was in an open channel. This is incorrect, but it ultimately doesn’t matter too much. The Reynolds number basically controls the “viscosity” of the simulation [or rather what the lay person would think of as viscosity]. The LBM is only stable for a range of Reynolds numbers. It also loses some physical accuracy for flow through very tiny openings (less than 5 lattice points).
There are a few ways you can tune the LBM. Grid resolution (changing the scale effectively changes the behavior of the entire simulation!), Reynolds number, speed of sound, and “relaxation time” (which basically shows how much time is passing per time step). The LBM is designed to be stable for mesoscopic scales (millimeters and centimeters). It is not usually used for simulating meters or kilometers, but thankfully we don’t need real physical accuracy. We just have to be able to fake it.
You can also model the effects of gravity on the fluids. I am going to investigate altering this technique so that each two dimensional cell gets the effects of gravity in direction and proportional to the gradient of the terrain at this location. It might be necessary to upsample the terrain for this purpose, and simulate multiple fluid update cycles per erosion update. We will see. I still don’t know exactly how this will work. Implementing gravity [probably from Stokup’s book on LBM] is what I will do immediately after fixing the improved (work-redundant) method so that it is actually an improvement.
I am also looking at sheer, stress, and strain measurements, discrete element modeling, and immersed boundary methods. I have not had a class on statistical mechanics or CFD yet, so I am really just poking around in the dark. However, I will continue to learn more about CFD, LBM or not, and have hopes of applying it to this field.
