https://davidar.io/post/sim-glsl

1<!DOCTYPE html><html lang="en">
2<head>
3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
4<meta name="viewport" content="width=device-width, initial-scale=1">
5<link rel="stylesheet" href="//texify.davidar.io/main.css" type="text/css">
6<title>Simulating worlds on the GPU: Four billion years in four minutes</title>
7<meta name="description" content="This post delves into the implementation of my procedural earth simulation, written entirely in GLSL fragment shaders. It simulates the complete history of an earth-like planet in a few minutes, with the simulation updating at 60 frames per second.">
8<meta name="author" content="David A Roberts">
9<link rel="icon" type="image/png" href="/favicon.png">
10<link rel="canonical" href="https://davidar.io/post/sim-glsl">
11
12<!-- Facebook Meta Tags -->
13<meta property="og:url" content="https://davidar.io/post/sim-glsl">
14<meta property="og:type" content="website">
15<meta property="og:title" content="Simulating worlds on the GPU: Four billion years in four minutes">
16<meta property="og:description" content="This post delves into the implementation of my procedural earth simulation, written entirely in GLSL fragment shaders. It simulates the complete history of an earth-like planet in a few minutes, with the simulation updating at 60 frames per second.">
17<meta property="og:image" content="https://davidar.io/img/sh18.jpg">
18
19<!-- Twitter Meta Tags -->
20<meta name="twitter:card" content="summary_large_image">
21<meta property="twitter:domain" content="davidar.io">
22<meta property="twitter:url" content="https://davidar.io/post/sim-glsl">
23<meta name="twitter:title" content="Simulating worlds on the GPU: Four billion years in four minutes">
24<meta name="twitter:description" content="This post delves into the implementation of my procedural earth simulation, written entirely in GLSL fragment shaders. It simulates the complete history of an earth-like planet in a few minutes, with the simulation updating at 60 frames per second.">
25<meta name="twitter:image" content="https://davidar.io/img/sh18.jpg">
26
27</head>
28<body>
29<header style="">
30    <h1>Simulating worlds on the GPU</h1>
31    <p>Four billion years in four minutes</p>
32    <address><a href="/">David A Roberts</a></address>
33</header>
34<main>
35<section id="abstract">
36    <h6>Abstract</h6>
37    <p>This post delves into the implementation of my <a href="https://www.shadertoy.com/view/XttcWn">procedural earth simulation</a>, written entirely in GLSL fragment shaders. It simulates the complete history of an earth-like planet in a few minutes, with the simulation updating at 60 frames per second.</p>
38</section>
39<figure>
40<div class="embed-16-9">
41    <!-- <iframe src="https://player.vimeo.com/video/283607168" width="640" height="360" frameborder="0" allow="autoplay; fullscreen; picture-in-picture" allowfullscreen></iframe> -->
42    <iframe width="1280" height="720" src="https://www.youtube.com/embed/h9mVtkzJkK4?vq=hd720p" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture; web-share" referrerpolicy="strict-origin-when-cross-origin" allowfullscreen></iframe>
43</div>
44<figcaption>A video recording of the <a href="https://www.shadertoy.com/view/XttcWn">final shader</a>.</figcaption>
45</figure>
46
47<h2 class="num">Protoplanet</h2>
48
49<blockquote>
50<p>This story begins four and a half billion years ago, with a lump of molten rock...
51</blockquote>
52
53<figure>
54<div class="embed-16-9"><iframe src="https://www.shadertoy.com/embed/XsGBDt" width="640" height="360"></iframe></div>
55</figure>
56
57<p>The early earth was a <a href="https://en.wikipedia.org/wiki/Protoplanet">protoplanet</a>, red hot and heavily cratered by asteroid impacts. As my earth simulation is <em>entirely procedurally generated</em>, with no pre-rendered textures, the first task is to generate a map of this terrain. To calculate the <code>height</code> of the terrain at a given <code>lat</code>itude and <code>lon</code>gitude, first translate to 3D cartesian coordinates:
58
59<pre><code class="glsl">
60vec3 p = 1.5 * vec3(
61    sin(lon*PI/180.) * cos(lat*PI/180.),

62    sin(lat*PI/180.),
63    cos(lon*PI/180.) * cos(lat*PI/180.));
64</code></pre>
65
66<p>Now, as asteroids come in a variety of sizes, so do the resulting craters. To accommodate this, the shader iterates over five levels of detail, layering craters of decreasing size over each other. To make the craters have a realistic rugged appearance, this is mixed with some <a href="https://iquilezles.untergrund.net/www/articles/fbm/fbm.htm">fractional Brownian motion</a> noise, and scaled so that the largest craters have the most impact on the terrain.
67
68<pre><code class="glsl">
69float height = 0.;
70for (float i = 0.; i &lt; 5.; i++) {
71    float c = craters(0.4 * pow(2.2, i) * p);
72    float noise = 0.4 * exp(-3. * c) * FBM(10. * p);
73    float w = clamp(3. * pow(0.4, i), 0., 1.);
74    height += w * (c + noise);
75}
76height = pow(height, 3.);
77</code></pre>
78
79<p>The craters themselves are generated on a 3D grid, from which a sphere is carved out for the surface terrain. To avoid visible regularity, the crater centres are given a pseudo-random offset from the grid points, using a <a href="https://www.shadertoy.com/view/4djSRW">hash function</a>. To calculate influence of a crater at a given location, take a weighted average of the craters belonging to the nearby grid points, with weights exponentially decreasing with distance from the centre. The crater rims are generated by a simple sine curve.
80
81<pre><code class="glsl">
82float craters(vec3 x) {
83    vec3 p = floor(x);
84    vec3 f = fract(x);
85    float va = 0.;
86    float wt = 0.;
87    for (int i = -2; i &lt;= 2; i++)
88     for (int j = -2; j &lt;= 2; j++)
89      for (int k = -2; k &lt;= 2; k++) {
90        vec3 g = vec3(i,j,k);
91        vec3 o = 0.8 * hash33(p + g);
92        float d = distance(f - g, o);
93        float w = exp(-4. * d);
94        va += w * sin(2.*PI * sqrt(d));
95        wt += w;
96    }
97    return abs(va / wt);
98}
99</code></pre>
100
101<p>The final procedurally generated heightmap looks like this:
102
103<figure>
104<img src="/img/protoplanet.jpg">
105</figure>
106
107<p>Although relatively simple, after filling the low-lying regions with water, this procedural terrain resembles what scientists believe the early earth actually looked like:
108
109<figure>
110<img src="/img/hadeanearth.jpg">
111<figcaption>Artistic impression of the early earth, by <a href="https://sservi.nasa.gov/articles/new-nasa-research-shows-giant-asteroids-battered-early-earth/">NASA</a>.</figcaption>
112</figure>
113
114<blockquote>
115<p>Water contained within was vaporised by the heat, which escaped and began circulating through the early atmosphere forming around the planet. As time progressed and the rock cooled, the water vapour began to condense into oceans. The flow of liquid water across the surface carved valleys in the terrain, leaving an accumulation of sediment in its wake.
116</blockquote>
117
118<h2 class="num">Tectonic plates</h2>
119
120<p>The formation of mountains, ocean trenches, and familiar continental landforms requires a model of tectonic movement. The simulation randomly generates seed locations for plates, with an initial velocity. These plates grow in size over time with a simple aggregation model, which randomly selects neighbouring points and adds them to a plate if they have not already been assigned to another plate. All of the pixels within a plate store the velocity of the plate's movement. The aggregation model is similar to that of a diffusion-limited aggregation (but without the diffusion):
121
122<figure>
123<div class="embed-16-9"><iframe src="https://www.shadertoy.com/embed/ldK3RW" width="640" height="360"></iframe></div>
124</figure>
125
126<p>Continuous movement of the plates is difficult, as it would require plate boundaries to account for movements measured in fractions of a pixel. To avoid this, the plates are instead moved at discrete time-steps, by a whole pixel either horizontally or vertically. These times are randomised for each plate such that the average velocity is maintained at the set speed and direction, and also so that it is unlikely that neighbouring plates will move simultaneously.
127
128<p>Plate collisions occur when some boundary pixels of one plate move onto a location previously occupied by pixels belonging to another plate. This causes <a href="https://en.wikipedia.org/wiki/Subduction">subduction</a>, which is modelled by simply slightly increasing the elevation of the terrain at the locations of the collision. Although this only occurs at the pixels along the boundary of a plate, the impact is gradually spread to neighbouring pixels through a simple thermal erosion model, which pu

128shes the elevation of a pixel in the direction of the average of its neighbours.
129
130<p>Altogether this provides a decent simulation of the formation of continents with mountain ranges (which will be further improved with the introduction of hydraulic erosion in the next section):
131
132<figure>
133<div class="embed-16-9"><iframe src="https://www.shadertoy.com/embed/XtffW8" width="640" height="360"></iframe></div>
134</figure>
135
136<h2 class="num">Hydraulic erosion</h2>
137
138<p>The rugged appearance of natural terrain is largely driven by the formation of river basins, which erode landscapes in a familiar branching pattern. A variety of water flow simulations are readily available for this task, but a difficulty here is that the resolution of the terrain map is quite low for an entire planet. Therefore, the model will have to be able to simulate rivers which are no more than a single pixel wide. <a href="https://arxiv.org/abs/1803.02977">Barnes (2018)</a> proposes a simple model which achieves just this.
139
140<p>Simply put, each pixel examines its eight neighbours, to determine which direction has the greatest decrease in elevation (adjusted for the fact that the diagonal neighbours are further away). This direction of greatest slope is where water flowing out of this pixel will travel. Water is initially distributed amongst cells by rainfall, which is then transported between neighbouring pixels at each time-step.
141
142<p>Erosion is driven by a <a href="https://en.wikipedia.org/wiki/Stream_power_law">stream power law</a>:
143
144<pre><code class="glsl">
145elevation -= 0.05 * pow(water, 0.8) * pow(slope, 2.);
146</code></pre>
147
148<p>Here we have the <code>elevation</code> and amount of <code>water</code> located at the current cell, along with the <code>slope</code> in the direction the water is travelling. The decrease in elevation is capped so that it doesn't become lower than the location the water is flowing to.
149
150<p>The interaction between the water flow and erosion results in the natural formation of river basins in the terrain:
151
152<figure>
153<div class="embed-16-9"><iframe src="https://www.shadertoy.com/embed/XsVBRm" width="640" height="360"></iframe></div>
154</figure>
155
156<p>By colouring connected waterways (with the colour determined by the location of the river's mouth), it's possible to produce striking visualisations reminiscent of <a href="https://imgur.com/gallery/WaEbi">real river basin maps</a>:
157
158<figure>
159<img src="/img/basin.png">
160<figcaption>Simulated river basins. <a href="https://www.shadertoy.com/view/XsVBDz">Original shader</a>.</figcaption>
161</figure>
162
163<figure>
164<img src="https://i.imgur.com/ZXLEvU3.jpg">
165<figcaption>River basins of USA, by <a href="https://www.grasshoppergeography.com/"></a>Grasshopper Geography</a>.</figcaption>
166</figure>
167
168<h2 class="num">Global climate</h2>
169
170<p>Simulating the climate system of an entire planet is a daunting task, but luckily it turns out that it can be approximated relatively easily. The driving force behind everything in my climate simulation is a procedurally generated map of the <a href="https://en.wikipedia.org/wiki/Atmospheric_pressure#Mean_sea-level_pressure">mean sea-level pressure (MSLP)</a>.
171
172<p>According to <a href="https://web.archive.org/web/20130619132254/http://jc.tech-galaxy.com/bricka/climate_cookbook.html">the Climate Cookbook</a>, the main ingredients in creating a MSLP map are where the landforms are located amidst the ocean, and the impact of latitude. In fact, if you take data from a real MSLP map of the Earth, separate out locations according to whether they are land or ocean, and plot the MSLP against latitude, you end up with two sinusoidal curves for the land and ocean with slightly different shapes.
173
174<p>By fitting the parameters appropriately, I came up with a crude model of the annual mean pressure (here the <code>lat</code>itude is measured in degrees):
175
176<pre><code class="glsl">
177if (land) {
178    mslp = 1012.5 - 6. * cos(lat*PI/45.);
179} else { // ocean
180    mslp = 1014.5 - 20. * cos(lat*PI/30.);
181}
182</code></pre>
183
184<p>Of course, this isn't quite enough to generate a realistic MSLP map, as generating values for the land and ocean separately results in sharp discontinuities at the boundaries between them. In reality, MSLP smoothly varies across the transition from ocean to land, due to the local diffusion of gas pressure. This diffusion process can be approximated quite well by simply applying a <a href="https://en.wikipedia.org/wiki/Gaussian_blur">Gaussian blur</a>

184 to the MSLP map (with a standard deviation of 10--15 degrees).
185
186<p>To allow for the climate to change along with the seasons, it's necessary to also model the difference in MSLP between January and July. Once again, terrestrial data suggests this follows a sinusoidal pattern. By fitting parameters and applying a Gaussian blur, this can be combined with the annual MSLP map to generate dynamic climate patterns which vary throughout the year.
187
188<pre><code class="glsl">
189if (land) {
190    delta = 15. * sin(lat*PI/90.);
191} else { // ocean
192    delta = 20. * sin(lat*PI/35.) * abs(lat)/90.;
193}
194</code></pre>
195
196<p>Now, with the MSLP in hand, it is possible to generate wind currents and temperatures. In reality it's the temperate which generates the pressure, but correlation is correlation. This requires a little more fiddling to generate realistic values (<code>season</code> oscillates between -1 and 1 throughout the year):
197
198<pre><code class="glsl">
199float temp = 40. * tanh(2.2 * exp(-0.5 * pow((lat + 5. * season)/30., 2.)))
200             - 15. - (mslp - 1012.) / 1.8 + 1.5 * land - 4. * elevation;
201</code></pre>
202
203<p>Wind tends to move from high-pressure to low, but at a global scale we also need to account for the <a href="https://en.wikipedia.org/wiki/Coriolis_force">Coriolis force</a>, which is responsible for causing winds to circulate <em>around</em> pressure zones (<code>grad</code> is the MSLP gradient vector):
204
205<pre><code class="glsl">
206vec2 coriolis = 15. * sin(lat*PI/180.) * vec2(-grad.y, grad.x);
207vec2 velocity = coriolis - grad;
208</code></pre>
209
210<p>Although a relatively crude simulation, this generates remarkably <a href="https://gist.github.com/davidar/229193b04bdb0dd8cba20dc31592625a">realistic</a> wind circulation patterns. If you look closely, you may notice a number of natural phenomena being replicated, including the reversal of winds over India during the monsoon season:
211
212<figure>
213<div class="embed-16-9"><iframe src="https://www.shadertoy.com/embed/MdGBWG" width="640" height="360"></iframe></div>
214</figure>
215
216<p>As a final detail, precipitation can be simulated by advecting water vapour from the ocean, through the wind vector field, and onto the land:
217
218<figure>
219<div class="embed-16-9"><iframe src="https://www.shadertoy.com/embed/MdKfWK" width="640" height="360"></iframe></div>
220</figure>
221
222<p>The advection is implemented in a similar manner to fluid simulations:
223
224<figure>
225<div class="embed-16-9"><iframe src="https://www.shadertoy.com/embed/XlsBDf" width="640" height="360"></iframe></div>
226</figure>
227
228<h2 class="num">Life</h2>
229
230<p>The climate influences the distribution of life on a planet. Rainfall patterns and temperature variation dictate rates of plant growth. As the seasons change, herbivores migrate to regions with enough vegetation to sustain them. And, as they follow the vegetation, predators follow them. All of these dynamics can be captured by a <a href="https://en.wikipedia.org/wiki/Lotka-Volterra_equations">Lotka--Volterra</a> diffusion model:
231
232<pre><code class="glsl">
233float dx = plant_growth - c.y;
234float dy = reproduction * c.x - predation * c.z - 1.;
235float dz = predation * c.y - 1.;
236float dt = 0.1;
237c.xyz += dt * c.xyz * vec3(dx, dy, dz);
238</code></pre>
239
240<p>The <code>xyz</code> elements of <code>c</code> represent the populations of vegetation, herbivores, and predators respectively. On a large scale, the dynamics of animal populations generate interesting patterns:
241
242<figure>
243<div class="embed-16-9"><iframe src="https://www.shadertoy.com/embed/Xtcyzr" width="640" height="360"></iframe></div>
244</figure>
245
246<p>In real life, these kinds of patterns are most easily seen with microbe populations in a petri dish, but the same laws govern large animal populations across the globe.
247
248<figure>
249<img src="/img/spiralmold.jpg">
250<figcaption><a href="http://www.evsc.net/projects/reaction-diffusion-2">Spiral waves in colonies of mold</a>.</figcaption>
251</figure>
252
253<h2 class="num">Humanity</h2>
254
255<blockquote>
256<p>Concluding the prelude on the early earth, the pace slows to a cycle between day and night, terrain becoming fixed as tectonic movements become imperceptible. Soon the night reveals unprecedented patterns of light, as humanity proceeds to colonise the surface of the planet.
257
258<p>This rapid expansion brings its own set of changes, as humans begin to burn large amounts of fossil fuels to power their settlements. Carbon that had lain dormant for millions of years is released into the atmosphere, and dispersed around the planet.
259
260<p>Over several hundred years, humans burn through all available fossil fuel resources, releasing five trillion tonnes of c

260arbon into the atmosphere. This strengthens the greenhouse effect, <a href="https://www.nature.com/articles/nclimate3036">raising the global average temperature by almost 10 degrees Celsius</a>. Large regions of land around the equator are rendered uninhabitable by extreme temperatures, resulting in the disappearance of humanity from a significant portion of the planet.
261</blockquote>
262
263</main>

264<script async src="//texify.davidar.io/load.js" type="71fc772daff89597ead7c8a9-text/javascript"></script>

264
265<!-- Cloudflare Pages Analytics -->

265<script defer src='https://static.cloudflareinsights.com/beacon.min.js' data-cf-beacon='{"token": "

26558aced9614cc4a259a465717f7347819

265"}' type="71fc772daff89597ead7c8a9-text/javascript"></script>

265<!-- Cloudflare Pages Analytics -->

265<script src="/cdn-cgi/scripts/7d0fa10a/cloudflare-static/rocket-loader.min.js" data-cf-settings="71fc772daff89597ead7c8a9-|49" defer></script>

265<script type="module" src="https://static.cloudflareinsights.com/beacon.min.js/v31edd6df95cf4e85bb4c19e7a9bdbcba1788362987495" integrity="sha512-iIg7k2xntmwu6/uSb5tpc/hySgZc4eoL31yB29W6tJFo2akwjPWcEqnCEdJvGexCL0KEQwVYv5BlowfhVz26hg==" data-cf-beacon='{"version":"2024.11.0","token":"10d49f17c8b24855a767c6b3d2b1b64c","spa":2}' crossorigin="anonymous"></script>

265
266</body>
267</html>