veloren_world/sim/erosion.rs
1use super::{diffusion, downhill, uphill};
2use crate::{config::CONFIG, util::RandomField};
3use common::{
4 terrain::{
5 MapSizeLg, NEIGHBOR_DELTA, TerrainChunkSize, neighbors, uniform_idx_as_vec2,
6 vec2_as_uniform_idx,
7 },
8 vol::RectVolSize,
9};
10use tracing::{debug, error, info, warn};
11// use faster::*;
12use itertools::izip;
13use noise::NoiseFn;
14use num::{Float, Zero};
15use ordered_float::{FloatCore, NotNan};
16use rayon::prelude::*;
17use std::{
18 cmp::{Ordering, Reverse},
19 collections::BinaryHeap,
20 fmt, mem,
21 time::Instant,
22};
23use vek::*;
24
25pub type Alt = f64;
26pub type Compute = f64;
27pub type Computex8 = [Compute; 8];
28
29/* code used by sharp in future
30/// Compute the water flux at all chunks, given a list of chunk indices sorted
31/// by increasing height.
32pub fn get_drainage(
33 map_size_lg: MapSizeLg,
34 newh: &[u32],
35 downhill: &[isize],
36 _boundary_len: usize,
37) -> Box<[f32]> {
38 // FIXME: Make the below work. For now, we just use constant flux.
39 // Initially, flux is determined by rainfall. We currently treat this as the
40 // same as humidity, so we just use humidity as a proxy. The total flux
41 // across the whole map is normalize to 1.0, and we expect the average flux
42 // to be 0.5. To figure out how far from normal any given chunk is, we use
43 // its logit. NOTE: If there are no non-boundary chunks, we just set
44 // base_flux to 1.0; this should still work fine because in that case
45 // there's no erosion anyway. let base_flux = 1.0 / ((map_size_lg.chunks_len())
46 // as f32);
47 let base_flux = 1.0;
48 let mut flux = vec![base_flux; map_size_lg.chunks_len()].into_boxed_slice();
49 newh.iter().rev().for_each(|&chunk_idx| {
50 let chunk_idx = chunk_idx as usize;
51 let downhill_idx = downhill[chunk_idx];
52 if downhill_idx >= 0 {
53 flux[downhill_idx as usize] += flux[chunk_idx];
54 }
55 });
56 flux
57}
58*/
59
60/// Compute the water flux at all chunks for multiple receivers, given a list of
61/// chunk indices sorted by increasing height and weights for each receiver.
62pub fn get_multi_drainage(
63 map_size_lg: MapSizeLg,
64 mstack: &[u32],
65 mrec: &[u8],
66 mwrec: &[Computex8],
67 _boundary_len: usize,
68) -> Box<[Compute]> {
69 // FIXME: Make the below work. For now, we just use constant flux.
70 // Initially, flux is determined by rainfall. We currently treat this as the
71 // same as humidity, so we just use humidity as a proxy. The total flux
72 // across the whole map is normalize to 1.0, and we expect the average flux
73 // to be 0.5. To figure out how far from normal any given chunk is, we use
74 // its logit. NOTE: If there are no non-boundary chunks, we just set
75 // base_flux to 1.0; this should still work fine because in that case
76 // there's no erosion anyway.
77 let base_area = 1.0;
78 let mut area = vec![base_area; map_size_lg.chunks_len()].into_boxed_slice();
79 mstack.iter().for_each(|&ij| {
80 let ij = ij as usize;
81 let wrec_ij = &mwrec[ij];
82 let area_ij = area[ij];
83 mrec_downhill(map_size_lg, mrec, ij).for_each(|(k, ijr)| {
84 area[ijr] += area_ij * wrec_ij[k];
85 });
86 });
87 area
88 /*
89 a=dx*dy*precip
90 do ij=1,nn
91 ijk=mstack(ij)
92 do k =1,mnrec(ijk)
93 a(mrec(k,ijk))=a(mrec(k,ijk))+a(ijk)*mwrec(k,ijk)
94 enddo
95 enddo
96 */
97}
98
99/// Kind of water on this tile.
100#[derive(Clone, Copy, Debug, PartialEq)]
101pub enum RiverKind {
102 Ocean,
103 Lake {
104 /// In addition to a downhill node (pointing to, eventually, the bottom
105 /// of the lake), each lake also has a "pass" that identifies
106 /// the direction out of which water should flow from this lake
107 /// if it is minimally flooded. While some lakes may be too full for
108 /// this to be the actual pass used by their enclosing lake, we
109 /// still use this as a way to make sure that lake connections
110 /// to rivers flow in the correct direction.
111 neighbor_pass_pos: Vec2<i32>,
112 },
113 /// River should be maximal.
114 River {
115 /// Dimensions of the river's cross-sectional area, as m × m (rivers are
116 /// approximated as an open rectangular prism in the direction
117 /// of the velocity vector).
118 cross_section: Vec2<f32>,
119 },
120}
121
122impl RiverKind {
123 pub fn is_ocean(&self) -> bool { matches!(*self, RiverKind::Ocean) }
124
125 pub fn is_river(&self) -> bool { matches!(*self, RiverKind::River { .. }) }
126
127 pub fn is_lake(&self) -> bool { matches!(*self, RiverKind::Lake { .. }) }
128}
129
130impl PartialOrd for RiverKind {
131 fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
132 match (*self, *other) {
133 (RiverKind::Ocean, RiverKind::Ocean) => Some(Ordering::Equal),
134 (RiverKind::Ocean, _) => Some(Ordering::Less),
135 (_, RiverKind::Ocean) => Some(Ordering::Greater),
136 (RiverKind::Lake { .. }, RiverKind::Lake { .. }) => None,
137 (RiverKind::Lake { .. }, _) => Some(Ordering::Less),
138 (_, RiverKind::Lake { .. }) => Some(Ordering::Greater),
139 (RiverKind::River { .. }, RiverKind::River { .. }) => None,
140 }
141 }
142}
143
144/// From velocity and cross_section we can calculate the volumetric flow rate,
145/// as the cross-sectional area times the velocity.
146///
147/// TODO: we might choose to include a curve for the river, as long as it didn't
148/// allow it to cross more than one neighboring chunk away. For now we defer
149/// this to rendering time.
150///
151/// NOTE: This structure is 57 (or more likely 64) bytes, which is kind of big.
152#[derive(Clone, Debug, Default)]
153pub struct RiverData {
154 /// A velocity vector (in m / minute, i.e. voxels / second from a game
155 /// perspective).
156 ///
157 /// TODO: To represent this in a better-packed way, use u8s instead (as
158 /// "f8s").
159 pub(crate) velocity: Vec3<f32>,
160 /// The computed derivative for the segment of river starting at this chunk
161 /// (and flowing downhill). Should be 0 at endpoints. For rivers with
162 /// more than one incoming segment, we weight the derivatives by flux
163 /// (cross-sectional area times velocity) which is correlated
164 /// with mass / second; treating the derivative as "velocity" with respect
165 /// to length along the river, we treat the weighted sum of incoming
166 /// splines as the "momentum", and can divide it by the total incoming
167 /// mass as a way to find the velocity of the center of mass. We can
168 /// then use this derivative to find a "tangent" for the incoming river
169 /// segment at this point, and as the linear part of the interpolating
170 /// spline at this point.
171 ///
172 /// Note that we aren't going to have completely smooth curves here anyway,
173 /// so we will probably end up applying a dampening factor as well
174 /// (maybe based on the length?) to prevent extremely wild oscillations.
175 pub(crate) spline_derivative: Vec2<f32>,
176 /// If this chunk is part of a river, this should be true. We can't just
177 /// compute this from the cross section because once a river becomes
178 /// visible, we want it to stay visible until it reaches its sink.
179 pub river_kind: Option<RiverKind>,
180 /// We also have a second record for recording any rivers in nearby chunks
181 /// that manage to intersect this chunk, though this is unlikely to
182 /// happen in current gameplay. This is because river areas are allowed
183 /// to cross arbitrarily many chunk boundaries, if they are wide enough.
184 /// In such cases we may choose to render the rivers as particularly deep in
185 /// those places.
186 pub(crate) neighbor_rivers: Vec<u32>,
187}
188
189impl RiverData {
190 pub fn is_ocean(&self) -> bool {
191 self.river_kind
192 .as_ref()
193 .map(RiverKind::is_ocean)
194 .unwrap_or(false)
195 }
196
197 pub fn is_river(&self) -> bool {
198 self.river_kind
199 .as_ref()
200 .map(RiverKind::is_river)
201 .unwrap_or(false)
202 }
203
204 pub fn is_lake(&self) -> bool {
205 self.river_kind
206 .as_ref()
207 .map(RiverKind::is_lake)
208 .unwrap_or(false)
209 }
210
211 pub fn near_river(&self) -> bool { self.is_river() || !self.neighbor_rivers.is_empty() }
212
213 pub fn near_water(&self) -> bool { self.near_river() || self.is_lake() || self.is_ocean() }
214}
215
216/// Draw rivers and assign them heights, widths, and velocities. Take some
217/// liberties with the constant factors etc. in order to make it more likely
218/// that we draw rivers at all.
219pub fn get_rivers<F: fmt::Debug + Float + Into<f64>, G: Float + Into<f64>>(
220 map_size_lg: MapSizeLg,
221 continent_scale_hack: f64,
222 newh: &[u32],
223 water_alt: &[F],
224 downhill: &[isize],
225 indirection: &[i32],
226 drainage: &[G],
227) -> Box<[RiverData]> {
228 // For continuity-preserving quadratic spline interpolation, we (appear to) need
229 // to build up the derivatives from the top down. Fortunately this
230 // computation seems tractable.
231
232 let mut rivers = vec![RiverData::default(); map_size_lg.chunks_len()].into_boxed_slice();
233 let neighbor_coef = TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
234 // (Roughly) area of a chunk, times minutes per second.
235 // NOTE: Clearly, this should "actually" be 1/60 (or maybe 1/64, if you want to
236 // retain powers of 2).
237 //
238 // But since we want rivers to form more often than they do in real life, we use
239 // this as a way to control the frequency of river formation. As grid_scale
240 // increases, mins_per_sec should decrease, until it hits 1 / 60 or 1/ 64.
241 // For example, if grid_scale is multiplied by 4, mins_per_sec should be
242 // multiplied by 1 / (4.0 * 4.0).
243 let mins_per_sec = 1.0 / (continent_scale_hack * continent_scale_hack)/*1.0 / 16.0*//*1.0 / 64.0*/;
244 let chunk_area_factor = neighbor_coef.x * neighbor_coef.y * mins_per_sec;
245 // NOTE: This technically makes us discontinuous, so we should be cautious about
246 // using this.
247 let derivative_divisor = 1.0;
248 newh.iter().rev().for_each(|&chunk_idx| {
249 let chunk_idx = chunk_idx as usize;
250 let downhill_idx = downhill[chunk_idx];
251 if downhill_idx < 0 {
252 // We are in the ocean.
253 debug_assert!(downhill_idx == -2);
254 rivers[chunk_idx].river_kind = Some(RiverKind::Ocean);
255 return;
256 }
257 let downhill_idx = downhill_idx as usize;
258 let downhill_pos = uniform_idx_as_vec2(map_size_lg, downhill_idx);
259 let dxy = (downhill_pos - uniform_idx_as_vec2(map_size_lg, chunk_idx)).map(|e| e as f64);
260 let neighbor_dim = neighbor_coef * dxy;
261 // First, we calculate the river's volumetric flow rate.
262 let chunk_drainage = drainage[chunk_idx].into();
263 // Volumetric flow rate is just the total drainage area to this chunk, times
264 // rainfall height per chunk per minute, times minutes per second
265 // (needed in order to use this as a m³ volume).
266 // TODO: consider having different rainfall rates (and including this
267 // information in the computation of drainage).
268 let volumetric_flow_rate =
269 chunk_drainage * chunk_area_factor * CONFIG.rainfall_chunk_rate as f64;
270 let downhill_drainage = drainage[downhill_idx].into();
271
272 // We know the drainage to the downhill node is just chunk_drainage - 1.0 (the
273 // amount of rainfall this chunk is said to get), so we don't need to
274 // explicitly remember the incoming mass. How do we find a slope for
275 // endpoints where there is no incoming data? Currently, we just assume
276 // it's set to 0.0. TODO: Fix this when we add differing amounts of
277 // rainfall.
278 let incoming_drainage = downhill_drainage - 1.0;
279 let get_river_spline_derivative =
280 |neighbor_dim: Vec2<f64>, spline_derivative: Vec2<f32>| {
281 // "Velocity of center of mass" of splines of incoming flows.
282 let river_prev_slope = spline_derivative.map(|e| e as f64);
283 // NOTE: We need to make sure the slope doesn't get *too* crazy.
284 // ((dpx - cx) - 4 * MAX).abs() = bx
285 // NOTE: This will fail if the distance between chunks in any direction
286 // is exactly TerrainChunkSize::RECT * 4.0, but hopefully this should not be
287 // possible. NOTE: This isn't measuring actual distance, you can
288 // go farther on diagonals.
289 let max_deriv = neighbor_dim - neighbor_coef * 2.0 * 2.0.sqrt();
290 let extra_divisor = river_prev_slope
291 .map2(max_deriv, |e, f| (e / f).abs())
292 .reduce_partial_max();
293 // Set up the river's spline derivative. For each incoming river at pos with
294 // river_spline_derivative bx, we can compute our interpolated slope as:
295 // d_x = 2 * (chunk_pos - pos - bx) + bx
296 // = 2 * (chunk_pos - pos) - bx
297 //
298 // which is exactly twice what was weighted by uphill nodes to get our
299 // river_spline_derivative in the first place.
300 //
301 // NOTE: this probably implies that the distance shouldn't be normalized, since
302 // the distances aren't actually equal between x and y... we'll
303 // see what happens.
304 (if extra_divisor > 1.0 {
305 river_prev_slope / extra_divisor
306 } else {
307 river_prev_slope
308 })
309 .map(|e| e as f32)
310 };
311
312 let river = &rivers[chunk_idx];
313 let river_spline_derivative =
314 get_river_spline_derivative(neighbor_dim, river.spline_derivative);
315
316 let indirection_idx = indirection[chunk_idx];
317 // Find the lake we are flowing into.
318 let lake_idx = if indirection_idx < 0 {
319 // If we're a lake bottom, our own indirection is negative.
320 let pass_idx = (-indirection_idx) as usize;
321 // NOTE: Must exist since this lake had a downhill in the first place.
322 let neighbor_pass_idx = downhill[pass_idx] as usize/*downhill_idx*/;
323 let lake_neighbor_pass = &mut rivers[neighbor_pass_idx];
324 // We definitely shouldn't have encountered this yet!
325 debug_assert!(lake_neighbor_pass.velocity == Vec3::zero());
326 // TODO: Rethink making the lake neighbor pass always a river or lake, no matter
327 // how much incoming water there is? Sometimes it looks weird
328 // having a river emerge from a tiny pool.
329 lake_neighbor_pass.river_kind = Some(RiverKind::River {
330 cross_section: Vec2::default(),
331 });
332 chunk_idx
333 } else {
334 indirection_idx as usize
335 };
336
337 // Find the pass this lake is flowing into (i.e. water at the lake bottom gets
338 // pushed towards the point identified by pass_idx).
339 let pass_idx = if downhill[lake_idx] < 0 {
340 // Flows into nothing, so this lake is its own pass.
341 lake_idx
342 } else {
343 (-indirection[lake_idx]) as usize
344 };
345
346 // Add our spline derivative to the downhill river (weighted by the chunk's
347 // drainage). NOTE: Don't add the spline derivative to the lake side of
348 // the pass for our own lake, because we don't want to preserve weird
349 // curvature from before we hit the lake in the outflowing river (this
350 // will not apply to one-chunk lakes, which are their own pass).
351 if pass_idx != downhill_idx {
352 // TODO: consider utilizing height difference component of flux as well;
353 // currently we just discard it in figuring out the spline's slope.
354 let downhill_river = &mut rivers[downhill_idx];
355 let weighted_flow = (neighbor_dim * 2.0 - river_spline_derivative.map(|e| e as f64))
356 / derivative_divisor
357 * chunk_drainage
358 / incoming_drainage;
359 downhill_river.spline_derivative += weighted_flow.map(|e| e as f32);
360 }
361
362 let neighbor_pass_idx = downhill[pass_idx];
363 // Find our own water height.
364 let chunk_water_alt = water_alt[chunk_idx];
365 if neighbor_pass_idx >= 0 {
366 // We may be a river. But we're not sure yet, since we still could be
367 // underwater. Check the lake height and see if our own water height is within
368 // ε of it.
369 let lake_water_alt = water_alt[lake_idx];
370 if chunk_water_alt == lake_water_alt {
371 // Not a river.
372 // Check whether we we are the lake side of the pass.
373 // NOTE: Safe because this is a lake.
374 let (neighbor_pass_pos, river_spline_derivative) = if pass_idx == chunk_idx {
375 // This is a pass, so set our flow direction to point to the neighbor pass
376 // rather than downhill.
377 // NOTE: Safe because neighbor_pass_idx >= 0.
378 (
379 uniform_idx_as_vec2(map_size_lg, downhill_idx),
380 river_spline_derivative,
381 )
382 } else {
383 // Try pointing towards the lake side of the pass.
384 (
385 uniform_idx_as_vec2(map_size_lg, pass_idx),
386 river_spline_derivative,
387 )
388 };
389 let lake = &mut rivers[chunk_idx];
390 lake.spline_derivative = river_spline_derivative;
391 lake.river_kind = Some(RiverKind::Lake {
392 neighbor_pass_pos: neighbor_pass_pos
393 * TerrainChunkSize::RECT_SIZE.map(|e| e as i32),
394 });
395 return;
396 }
397 // Otherwise, we must be a river.
398 } else {
399 // We are flowing into the ocean.
400 debug_assert!(neighbor_pass_idx == -2);
401 // But we are not the ocean, so we must be a river.
402 }
403 // Now, we know we are a river *candidate*. We still don't know whether we are
404 // actually a river, though. There are two ways for that to happen:
405 // (i) We are already a river.
406 // (ii) Using the Gauckler–Manning–Strickler formula for cross-sectional
407 // average velocity of water, we establish that the river can be
408 // "big enough" to appear on the Veloren map.
409 //
410 // This is very imprecise, of course, and (ii) may (and almost certainly will)
411 // change over time.
412 //
413 // In both cases, we preemptively set our child to be a river, to make sure we
414 // have an unbroken stream. Also in both cases, we go to the effort of
415 // computing an effective water velocity vector and cross-sectional
416 // dimensions, as well as figuring out the derivative of our
417 // interpolating spline (since this percolates through the whole river
418 // network).
419 let downhill_water_alt = water_alt[downhill_idx];
420 let neighbor_distance = neighbor_dim.magnitude();
421 let dz = (downhill_water_alt - chunk_water_alt).into();
422 let slope = dz.abs() / neighbor_distance;
423 if slope == 0.0 {
424 // This is not a river--how did this even happen?
425 let pass_idx = (-indirection_idx) as usize;
426 error!(
427 "Our chunk (and downhill, lake, pass, neighbor_pass): {:?} (to {:?}, in {:?} via \
428 {:?} to {:?}), chunk water alt: {:?}, lake water alt: {:?}",
429 uniform_idx_as_vec2(map_size_lg, chunk_idx),
430 uniform_idx_as_vec2(map_size_lg, downhill_idx),
431 uniform_idx_as_vec2(map_size_lg, lake_idx),
432 uniform_idx_as_vec2(map_size_lg, pass_idx),
433 if neighbor_pass_idx >= 0 {
434 Some(uniform_idx_as_vec2(map_size_lg, neighbor_pass_idx as usize))
435 } else {
436 None
437 },
438 water_alt[chunk_idx],
439 water_alt[lake_idx]
440 );
441 panic!("Should this happen at all?");
442 }
443 let slope_sqrt = slope.sqrt();
444 // Now, we compute a quantity that is proportional to the velocity of the chunk,
445 // derived from the Manning formula, equal to
446 // volumetric_flow_rate / slope_sqrt * CONFIG.river_roughness.
447 let almost_velocity = volumetric_flow_rate / slope_sqrt * CONFIG.river_roughness as f64;
448 // From this, we can figure out the width of the chunk if we know the height.
449 // For now, we hardcode the height to 0.5, but it should almost
450 // certainly be much more complicated than this.
451 // let mut height = 0.5f32;
452 // We approximate the river as a rectangular prism. Theoretically, we need to
453 // solve the following quintic equation to determine its width from its
454 // height:
455 //
456 // h^5 * w^5 = almost_velocity^3 * (w + 2 * h)^2.
457 //
458 // This is because one of the quantities in the Manning formula (the unknown) is
459 // R_h = (area of cross-section / h)^(2/3).
460 //
461 // Unfortunately, quintic equations do not in general have algebraic solutions,
462 // and it's not clear (to me anyway) that this one does in all cases.
463 //
464 // In practice, for high ratios of width to height, we can approximate the
465 // rectangular prism's perimeter as equal to its width, so R_h as equal
466 // to height. This greatly simplifies the calculation. For simplicity,
467 // we do this even for low ratios of width to height--I found that for
468 // most real rivers, at least big ones, this approximation is
469 // "good enough." We don't need to be *that* realistic :P
470 //
471 // NOTE: Derived from a paper on estimating river width.
472 let mut width = 5.0
473 * (CONFIG.river_width_to_depth as f64
474 * (CONFIG.river_width_to_depth as f64 + 2.0).powf(2.0 / 3.0))
475 .powf(3.0 / 8.0)
476 * volumetric_flow_rate.powf(3.0 / 8.0)
477 * slope.powf(-3.0 / 16.0)
478 * (CONFIG.river_roughness as f64).powf(3.0 / 8.0);
479 width = width.max(0.0);
480
481 let mut height = if width == 0.0 {
482 CONFIG.river_min_height as f64
483 } else {
484 (almost_velocity / width).powf(3.0 / 5.0)
485 };
486
487 // We can now weight the river's drainage by its direction, which we use to help
488 // improve the slope of the downhill node.
489 let river_direction = Vec3::new(neighbor_dim.x, neighbor_dim.y, dz.signum() * dz);
490
491 // Now, we can check whether this is "really" a river.
492 // Currently, we just check that width and height are at least 0.5 and
493 // CONFIG.river_min_height.
494 let river = &rivers[chunk_idx];
495 let is_river = river.is_river() || width >= 0.5 && height >= CONFIG.river_min_height as f64;
496 let downhill_river = &mut rivers[downhill_idx];
497
498 if is_river {
499 // Provisionally make the downhill chunk a river as well.
500 downhill_river.river_kind = Some(RiverKind::River {
501 cross_section: Vec2::default(),
502 });
503
504 // Additionally, if the cross-sectional area for this river exceeds the max
505 // river width, the river is overflowing the two chunks adjacent to
506 // it, which we'd prefer to avoid since only its two immediate
507 // neighbors (orthogonal to the downhill direction) are guaranteed
508 // uphill of it. Solving this properly most likely requires
509 // modifying the erosion model to take channel width into account,
510 // which is a formidable task that likely requires rethinking the
511 // current grid-based erosion model (or at least, requires tracking some
512 // edges that aren't implied by the grid graph). For now, we will solve this
513 // problem by making the river deeper when it hits the max width,
514 // until it consumes all the available energy in this part of the
515 // river.
516 let max_width = TerrainChunkSize::RECT_SIZE.x as f64 * CONFIG.river_max_width as f64;
517 if width > max_width {
518 width = max_width;
519 height = (almost_velocity / width).powf(3.0 / 5.0);
520 }
521 }
522 // Now we can compute the river's approximate velocity magnitude as well, as
523 let velocity_magnitude =
524 1.0 / CONFIG.river_roughness as f64 * height.powf(2.0 / 3.0) * slope_sqrt;
525
526 // Set up the river's cross-sectional area.
527 let cross_section = Vec2::new(width as f32, height as f32);
528 // Set up the river's velocity vector.
529 let mut velocity = river_direction;
530 velocity.normalize();
531 velocity *= velocity_magnitude;
532
533 let river = &mut rivers[chunk_idx];
534 // NOTE: Not trying to do this more cleverly because we want to keep the river's
535 // neighbors. TODO: Actually put something in the neighbors.
536 river.velocity = velocity.map(|e| e as f32);
537 river.spline_derivative = river_spline_derivative;
538 river.river_kind = if is_river {
539 Some(RiverKind::River { cross_section })
540 } else {
541 None
542 };
543 });
544 rivers
545}
546
547/// Precompute the maximum slope at all points.
548///
549/// TODO: See if allocating in advance is worthwhile.
550fn get_max_slope(
551 map_size_lg: MapSizeLg,
552 h: &[Alt],
553 rock_strength_nz: &(impl NoiseFn<f64, 3> + Sync),
554 height_scale: impl Fn(usize) -> Alt + Sync,
555) -> Box<[f64]> {
556 let min_max_angle = (15.0 / 360.0 * 2.0 * std::f64::consts::PI).tan();
557 let max_max_angle = (60.0 / 360.0 * 2.0 * std::f64::consts::PI).tan();
558 let max_angle_range = max_max_angle - min_max_angle;
559 h.par_iter()
560 .enumerate()
561 .map(|(posi, &z)| {
562 let wposf = uniform_idx_as_vec2(map_size_lg, posi).map(|e| e as f64)
563 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
564 let height_scale = height_scale(posi);
565 let wposz = z / height_scale;
566 // Normalized to be between 6 and and 54 degrees.
567 let div_factor = (2.0 * TerrainChunkSize::RECT_SIZE.x as f64) / 8.0;
568 let rock_strength = rock_strength_nz.get([wposf.x, wposf.y, wposz * div_factor]);
569 let rock_strength = rock_strength.clamp(-1.0, 1.0) * 0.5 + 0.5;
570 // Logistic regression. Make sure x ∈ (0, 1).
571 let logit = |x: f64| x.ln() - (-x).ln_1p();
572 // 0.5 + 0.5 * tanh(ln(1 / (1 - 0.1) - 1) / (2 * (sqrt(3)/pi)))
573 let logistic_2_base = 3.0f64.sqrt() * std::f64::consts::FRAC_2_PI;
574 // Assumes μ = 0, σ = 1
575 let logistic_cdf = |x: f64| (x / logistic_2_base).tanh() * 0.5 + 0.5;
576
577 // We do log-odds against center, so that our log odds are 0 when x = 0.25,
578 // lower when x is lower, and higher when x is higher.
579 //
580 // (NOTE: below sea level, we invert it).
581 //
582 // TODO: Make all this stuff configurable... but honestly, it's so complicated
583 // that I'm not sure anyone would be able to usefully tweak them on
584 // a per-map basis? Plus it's just a hacky heuristic anyway.
585 let center = 0.4;
586 let dmin = center - 0.05;
587 let dmax = center + 0.05;
588 let log_odds = |x: f64| logit(x) - logit(center);
589 let rock_strength = logistic_cdf(
590 1.0 * logit(rock_strength.clamp(1e-7, 1.0f64 - 1e-7))
591 + 1.0
592 * log_odds(
593 (wposz / CONFIG.mountain_scale as f64)
594 .abs()
595 .clamp(dmin, dmax),
596 ),
597 );
598 // NOTE: If you want to disable varying rock strength entirely, uncomment this
599 // line. let max_slope = 3.0.sqrt() / 3.0;
600 rock_strength * max_angle_range + min_max_angle //max_slope
601 })
602 .collect::<Vec<_>>()
603 .into_boxed_slice()
604}
605
606// simd alternative
607#[cfg(not(feature = "simd"))]
608#[derive(Copy, Clone)]
609struct M32(u32);
610#[cfg(not(feature = "simd"))]
611impl M32 {
612 #[inline]
613 fn splat(x: bool) -> Self { if x { Self(u32::MAX) } else { Self(u32::MIN) } }
614
615 #[inline]
616 fn any(&self) -> bool { self.0 != 0 }
617}
618#[cfg(feature = "simd")]
619type M32 = std::simd::Mask<i32, 1>;
620
621/// Erode all chunks by amount.
622///
623/// Our equation is:
624///
625/// dh(p) / dt = uplift(p)−k * A(p)^m * slope(p)^n
626///
627/// where A(p) is the drainage area at p, m and n are constants
628/// (we choose m = 0.4 and n = 1), and k is a constant. We choose
629///
630/// k = 2.244 * uplift.max() / (desired_max_height)
631///
632/// since this tends to produce mountains of max height desired_max_height;
633/// and we set desired_max_height = 1.0 to reflect limitations of mountain
634/// scale.
635///
636/// This algorithm does this in four steps:
637///
638/// 1. Sort the nodes in h by height (so the lowest node by altitude is first in
639/// the list, and the highest node by altitude is last).
640/// 2. Iterate through the list in *reverse.* For each node, we compute its
641/// drainage area as the sum of the drainage areas of its "children" nodes
642/// (i.e. the nodes with directed edges to this node). To do this
643/// efficiently, we start with the "leaves" (the highest nodes), which have
644/// no neighbors higher than them, hence no directed edges to them. We add
645/// their area to themselves, and then to all neighbors that they flow into
646/// (their "ancestors" in the flow graph); currently, this just means the
647/// node immediately downhill of this node. As we go lower, we know that all
648/// our "children" already had their areas computed, which means that we can
649/// repeat the process in order to derive all the final areas.
650/// 3. Now, iterate through the list in *order.* Whether we used the filling
651/// method to compute a "filled" version of each depression, or used the lake
652/// connection algorithm described in [1], each node is guaranteed to have
653/// zero or one drainage edges out, representing the direction of water flow
654/// for that node. For nodes i with zero drainage edges out (boundary nodes
655/// and lake bottoms) we set the slope to 0 (so the change in altitude is
656/// uplift(i)).
657///
658/// For nodes with at least one drainage edge out, we take
659/// advantage of the fact that we are computing new heights in order and
660/// rewrite our equation as (letting j = downhill[i], A[i] be the computed
661/// area of point i, p(i) be the x-y position of point i, flux(i) = k *
662/// A[i]^m / ((p(i) - p(j)).magnitude()), and δt = 1):
663///
664/// h[i](t + dt) = h[i](t) + δt * (uplift[i] + flux(i) * h[j](t + δt)) / (1 +
665/// flux(i) * δt).
666///
667/// Since we compute heights in ascending order by height, and j is downhill
668/// from i, h[j] will always be the *new* h[j](t + δt), while h[i] will still
669/// not have been computed yet, so we only need to visit each node once.
670///
671/// Afterwards, we also apply a hillslope diffusion process using an ADI
672/// (alternating direction implicit) method:
673///
674/// https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/Diffusion.f90
675///
676/// We also borrow the implementation for sediment transport from
677///
678/// https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/StreamPowerLaw.f90
679///
680/// The approximate equation for soil production function (predicting the rate
681/// at which bedrock turns into soil, increasing the distance between the
682/// basement and altitude) is taken from equation (11) from [2]. This (among
683/// numerous other sources) also includes at least one prediction that hillslope
684/// diffusion should be nonlinear, which we sort of attempt to approximate.
685///
686/// [1] Guillaume Cordonnier, Jean Braun, Marie-Paule Cani,
687/// Bedrich Benes, Eric Galin, et al..
688/// Large Scale Terrain Generation from Tectonic Uplift and Fluvial Erosion.
689/// Computer Graphics Forum, Wiley, 2016, Proc. EUROGRAPHICS 2016, 35 (2),
690/// pp.165-175. ⟨10.1111/cgf.12820⟩. ⟨hal-01262376⟩
691///
692/// [2] William E. Dietrich, Dino G. Bellugi, Leonard S. Sklar,
693/// Jonathan D. Stock
694/// Geomorphic Transport Laws for Predicting Landscape Form and Dynamics.
695/// Prediction in Geomorphology, Geophysical Monograph 135.
696/// Copyright 2003 by the American Geophysical Union
697/// 10.1029/135GM09
698#[expect(clippy::too_many_arguments)]
699fn erode(
700 // Underlying map dimensions.
701 map_size_lg: MapSizeLg,
702 // Height above sea level of topsoil
703 h: &mut [Alt],
704 // Height above sea level of bedrock
705 b: &mut [Alt],
706 // Height above sea level of water
707 wh: &mut [Alt],
708 max_uplift: f32,
709 max_g: f32,
710 kdsed: f64,
711 _seed: &RandomField,
712 rock_strength_nz: &(impl NoiseFn<f64, 3> + Sync),
713 uplift: impl Fn(usize) -> f32 + Sync,
714 n_f: impl Fn(usize) -> f32 + Sync,
715 m_f: impl Fn(usize) -> f32 + Sync,
716 kf: impl Fn(usize) -> f64 + Sync,
717 kd: impl Fn(usize) -> f64,
718 g: impl Fn(usize) -> f32 + Sync,
719 epsilon_0: impl Fn(usize) -> f32 + Sync,
720 alpha: impl Fn(usize) -> f32 + Sync,
721 is_ocean: impl Fn(usize) -> bool + Sync,
722 // scaling factors
723 height_scale: impl Fn(f32) -> Alt + Sync,
724 k_da_scale: impl Fn(f64) -> f64,
725 threadpool: &rayon::ThreadPool,
726) {
727 let compute_stats = true;
728 debug!("Done draining...");
729 // NOTE: To experimentally allow erosion to proceed below sea level, replace 0.0
730 // with -<Alt as Float>::infinity().
731 let min_erosion_height = 0.0; // -<Alt as Float>::infinity();
732
733 // NOTE: The value being divided by here sets the effective maximum uplift rate,
734 // as everything is scaled to it!
735 let dt = max_uplift as f64 / 1e-3;
736 debug!(?dt, "");
737 // Minimum sediment thickness before we treat erosion as sediment based.
738 let sediment_thickness = |_n| /*6.25e-5*/1.0e-4 * dt;
739 let neighbor_coef = TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
740 let chunk_area = neighbor_coef.x * neighbor_coef.y;
741 let min_length = neighbor_coef.reduce_partial_min();
742 let max_stable = min_length * min_length / dt;
743
744 // Debris flow area coefficient (m^(-2q)).
745 let q = 0.2;
746 // NOTE: Set to 1.0 to make (assuming n = 1) the erosion equation linear during
747 // each stream power iteration. This will result in significant speedups,
748 // at the cost of less interesting erosion behavior (linear vs. nonlinear).
749 let q_ = 1.5;
750 let k_da = 2.5 * k_da_scale(q);
751 let nx = usize::from(map_size_lg.chunks().x);
752 let ny = usize::from(map_size_lg.chunks().y);
753 let dx = TerrainChunkSize::RECT_SIZE.x as f64;
754 let dy = TerrainChunkSize::RECT_SIZE.y as f64;
755
756 // ```ignore
757 // ε₀ and α are part of the soil production approximate
758 // equation:
759 //
760 // -∂z_b / ∂t = ε₀ * e^(-αH)
761 //
762 // where
763 // z_b is the elevation of the soil-bedrock interface (i.e. the basement),
764 // ε₀ is the production rate of exposed bedrock (H = 0),
765 // H is the soil thickness normal to the ground surface,
766 // and α is a parameter (units of 1 / length).
767 //
768 // Note that normal depth at i, for us, will be interpreted as the soil depth vector,
769 // sed_i = ((0, 0), h_i - b_i),
770 // projected onto the ground surface slope vector,
771 // ground_surface_i = ((p_i - p_j), h_i - h_j),
772 // yielding the soil depth vector
773 // H_i = sed_i - sed_i ⋅ ground_surface_i / (ground_surface_i ⋅ ground_surface_i) * ground_surface_i
774 //
775 // = ((0, 0), h_i - b_i) -
776 // (0 * ||p_i - p_j|| + (h_i - b_i) * (h_i - h_j)) / (||p_i - p_j||^2 + (h_i - h_j)^2)
777 // * (p_i - p_j, h_i - h_j)
778 // = ((0, 0), h_i - b_i) -
779 // ((h_i - b_i) * (h_i - h_j)) / (||p_i - p_j||^2 + (h_i - h_j)^2)
780 // * (p_i - p_j, h_i - h_j)
781 // = (h_i - b_i) *
782 // (((0, 0), 1) - (h_i - h_j) / (||p_i - p_j||^2 + (h_i - h_j)^2) * (p_i - p_j, h_i - h_j))
783 // H_i_fact = (h_i - h_j) / (||p_i - p_j||^2 + (h_i - h_j)^2)
784 // H_i = (h_i - b_i) * ((((0, 0), 1) - H_i_fact * (p_i - p_j, h_i - h_j)))
785 // = (h_i - b_i) * (-H_i_fact * (p_i - p_j), 1 - H_i_fact * (h_i - h_j))
786 // ||H_i|| = (h_i - b_i) * √(H_i_fact^2 * ||p_i - p_j||^2 + (1 - H_i_fact * (h_i - h_j))^2)
787 // = (h_i - b_i) * √(H_i_fact^2 * ||p_i - p_j||^2 + 1 - 2 * H_i_fact * (h_i - h_j) +
788 // H_i_fact^2 * (h_i - h_j)^2)
789 // = (h_i - b_i) * √(H_i_fact^2 * (||p_i - p_j||^2 + (h_i - h_j)^2) +
790 // 1 - 2 * H_i_fact * (h_i - h_j))
791 // = (h_i - b_i) * √((h_i - h_j)^2 / (||p_i - p_j||^2 + (h_i - h_j)^2) * (||p_i - p_j||^2 + (h_i - h_j)^2) +
792 // 1 - 2 * (h_i - h_j)^2 / (||p_i - p_j||^2 + (h_i - h_j)^2))
793 // = (h_i - b_i) * √((h_i - h_j)^2 - 2(h_i - h_j)^2 / (||p_i - p_j||^2 + (h_i - h_j)^2) + 1)
794 //
795 // where j is i's receiver and ||p_i - p_j|| is the horizontal displacement between them. The
796 // idea here is that we first compute the hypotenuse between the horizontal and vertical
797 // displacement of ground (getting the horizontal component of the triangle), and then this is
798 // taken as one of the non-hypotenuse sides of the triangle whose other non-hypotenuse side is
799 // the normal height H_i, while their square adds up to the vertical displacement (h_i - b_i).
800 // If h and b have different slopes, this may not work completely correctly, but this is
801 // probably fine as an approximation.
802 // ```
803
804 // Spatio-temporal variation in net precipitation rate ((m / year) / (m / year))
805 // (ratio of precipitation rate at chunk relative to mean precipitation rate
806 // p₀).
807 let p = 1.0;
808 // Dimensionless multiplier for stream power erosion constant when land becomes
809 // sediment.
810 let k_fs_mult_sed = 4.0;
811 // Dimensionless multiplier for G when land becomes sediment.
812 let g_fs_mult_sed = 1.0;
813 let ((dh, newh, maxh, mrec, mstack, mwrec, area), (mut max_slopes, h_t)) = threadpool.join(
814 || {
815 let mut dh = downhill(
816 map_size_lg,
817 |posi| h[posi],
818 |posi| is_ocean(posi) && h[posi] <= 0.0,
819 );
820 debug!("Computed downhill...");
821 let (boundary_len, _indirection, newh, maxh) =
822 get_lakes(map_size_lg, |posi| h[posi], &mut dh);
823 debug!("Got lakes...");
824 let (mrec, mstack, mwrec) = get_multi_rec(
825 map_size_lg,
826 |posi| h[posi],
827 &dh,
828 &newh,
829 wh,
830 nx,
831 ny,
832 dx as Compute,
833 dy as Compute,
834 maxh,
835 threadpool,
836 );
837 debug!("Got multiple receivers...");
838 // TODO: Figure out how to switch between single-receiver and multi-receiver
839 // drainage, as the former is much less computationally costly.
840 // let area = get_drainage(map_size_lg, &newh, &dh, boundary_len);
841 let area = get_multi_drainage(map_size_lg, &mstack, &mrec, &mwrec, boundary_len);
842 debug!("Got flux...");
843 (dh, newh, maxh, mrec, mstack, mwrec, area)
844 },
845 || {
846 threadpool.join(
847 || {
848 let max_slope = get_max_slope(map_size_lg, h, rock_strength_nz, |posi| {
849 height_scale(n_f(posi))
850 });
851 debug!("Got max slopes...");
852 max_slope
853 },
854 || h.to_vec().into_boxed_slice(),
855 )
856 },
857 );
858
859 assert!(h.len() == dh.len() && dh.len() == area.len());
860
861 // max angle of slope depends on rock strength, which is computed from noise
862 // function. TODO: Make more principled.
863 let mid_slope = (30.0 / 360.0 * 2.0 * std::f64::consts::PI).tan();
864
865 type SimdType = f32;
866 type MaskType = M32;
867
868 // Precompute factors for Stream Power Law.
869 let czero = <SimdType as Zero>::zero();
870 let (k_fs_fact, k_df_fact) = threadpool.join(
871 || {
872 dh.par_iter()
873 .enumerate()
874 .map(|(posi, &posj)| {
875 let mut k_tot = [czero; 8];
876 if posj < 0 {
877 // Egress with no outgoing flows, no stream power.
878 k_tot
879 } else {
880 let old_b_i = b[posi];
881 let sed = h_t[posi] - old_b_i;
882 let n = n_f(posi);
883 // Higher rock strength tends to lead to higher erodibility?
884 let kd_factor = 1.0;
885 let k_fs = kf(posi) / kd_factor;
886
887 let k = if sed > sediment_thickness(n) {
888 // Sediment
889 k_fs_mult_sed * k_fs
890 } else {
891 // Bedrock
892 k_fs
893 } * dt;
894 let n = n as f64;
895 let m = m_f(posi) as f64;
896
897 let mwrec_i = &mwrec[posi];
898 mrec_downhill(map_size_lg, &mrec, posi).for_each(|(kk, posj)| {
899 let dxy = (uniform_idx_as_vec2(map_size_lg, posi)
900 - uniform_idx_as_vec2(map_size_lg, posj))
901 .map(|e| e as f64);
902 let neighbor_distance = (neighbor_coef * dxy).magnitude();
903 let knew = (k * (p * chunk_area * (area[posi] * mwrec_i[kk])).powf(m)
904 / neighbor_distance.powf(n))
905 as SimdType;
906 k_tot[kk] = knew;
907 });
908 k_tot
909 }
910 })
911 .collect::<Vec<[SimdType; 8]>>()
912 },
913 || {
914 dh.par_iter()
915 .enumerate()
916 .map(|(posi, &posj)| {
917 let mut k_tot = [czero; 8];
918 let uplift_i = uplift(posi) as f64;
919 debug_assert!(uplift_i.is_normal() && uplift_i > 0.0 || uplift_i == 0.0);
920 if posj < 0 {
921 // Egress with no outgoing flows, no stream power.
922 k_tot
923 } else {
924 let area_i = area[posi];
925 let max_slope = max_slopes[posi];
926 let chunk_area_pow = chunk_area.powf(q);
927
928 let old_b_i = b[posi];
929 let sed = h_t[posi] - old_b_i;
930 let n = n_f(posi);
931 let g_i = if sed > sediment_thickness(n) {
932 // Sediment
933 (g_fs_mult_sed * g(posi)) as f64
934 } else {
935 // Bedrock
936 g(posi) as f64
937 };
938
939 // Higher rock strength tends to lead to higher curvature?
940 let kd_factor = (max_slope / mid_slope).powi(2);
941 let k_da = k_da * kd_factor;
942
943 let mwrec_i = &mwrec[posi];
944 mrec_downhill(map_size_lg, &mrec, posi).for_each(|(kk, posj)| {
945 let mwrec_kk = mwrec_i[kk];
946
947 #[rustfmt::skip]
948 // Working equation:
949 // U = uplift per time
950 // D = sediment deposition per time
951 // E = fluvial erosion per time
952 // 0 = U + D - E - k_df * (1 + k_da * (mrec_kk * A)^q) * (∂B/∂p)^(q_)
953 //
954 // k_df = (U + D - E) / (1 + k_da * (mrec_kk * A)^q) / (∂B/∂p)^(q_)
955 //
956 // Want: ∂B/∂p = max slope at steady state, i.e.
957 // ∂B/∂p = max_slope
958 // Then:
959 // k_df = (U + D - E) / (1 + k_da * (mrec_kk * A)^q) / max_slope^(q_)
960 // Letting
961 // k = k_df * Δt
962 // we get:
963 // k = (U + D - E)Δt / (1 + k_da * (mrec_kk * A)^q) / (ΔB)^(q_)
964 //
965 // Now ∂B/∂t under constant uplift, without debris flow (U + D - E), is
966 // U + D - E = U - E + G/(p̃A) * ∫_A((U - ∂h/∂t) * dA)
967 //
968 // Observing that at steady state ∂h/∂t should theoretically
969 // be 0, we can simplify to:
970 // U + D = U + G/(p̃A) * ∫_A(U * dA)
971 //
972 // Upper bounding this at uplift = max_uplift/∂t for the whole prior
973 // drainage area, and assuming we account for just mrec_kk of
974 // the combined uplift and deposition, we get:
975 //
976 // U + D ≤ mrec_kk * U + G/p̃ * max_uplift/∂t
977 // (U + D - E)Δt ≤ (mrec_kk * uplift_i + G/p̃ * mrec_kk * max_uplift - EΔt)
978 //
979 // therefore
980 // k * (1 + k_da * (mrec_kk * A)^q) * max_slope^(q_) ≤ (mrec_kk * (uplift_i + G/p̃ * max_uplift) - EΔt)
981 // i.e.
982 // k ≤ (mrec_kk * (uplift_i + G/p̃ * max_uplift) - EΔt) / (1 + k_da * (mrec_kk * A)^q) / max_slope^q_
983 //
984 // (eliminating EΔt maintains the sign, but it's somewhat imprecise;
985 // we can address this later, e.g. by assigning a debris flow / fluvial erosion ratio).
986 let chunk_neutral_area = 0.1e6; // 1 km^2 * (1000 m / km)^2 = 1e6 m^2
987 let k = (mwrec_kk * (uplift_i + max_uplift as f64 * g_i / p))
988 / (1.0 + k_da * (mwrec_kk * chunk_neutral_area).powf(q))
989 / max_slope.powf(q_);
990 let dxy = (uniform_idx_as_vec2(map_size_lg, posi)
991 - uniform_idx_as_vec2(map_size_lg, posj))
992 .map(|e| e as f64);
993 let neighbor_distance = (neighbor_coef * dxy).magnitude();
994
995 let knew = (k
996 * (1.0 + k_da * chunk_area_pow * (area_i * mwrec_kk).powf(q))
997 / neighbor_distance.powf(q_))
998 as SimdType;
999 k_tot[kk] = knew;
1000 });
1001 k_tot
1002 }
1003 })
1004 .collect::<Vec<[SimdType; 8]>>()
1005 },
1006 );
1007
1008 debug!("Computed stream power factors...");
1009
1010 let mut lake_water_volume: Box<[Compute]> =
1011 vec![0.0_f64; map_size_lg.chunks_len()].into_boxed_slice();
1012 let mut elev: Box<[Compute]> = vec![0_f64; map_size_lg.chunks_len()].into_boxed_slice();
1013 let mut h_p: Box<[Compute]> = vec![0_f64; map_size_lg.chunks_len()].into_boxed_slice();
1014 let mut deltah: Box<[Compute]> = vec![0_f64; map_size_lg.chunks_len()].into_boxed_slice();
1015
1016 // calculate the elevation / SPL, including sediment flux
1017 let tol1: Compute = 1.0e-4_f64 * (maxh as Compute + 1.0);
1018 let tol2: Compute = 1.0e-3_f64 * (max_uplift as Compute + 1.0);
1019 let tol = tol1.max(tol2);
1020 let mut err = 2.0 * tol;
1021
1022 // Some variables for tracking statistics, currently only for debugging
1023 // purposes.
1024 let mut minh = <Alt as Float>::infinity();
1025 let mut maxh = 0.0;
1026 let mut nland = 0usize;
1027 let mut ncorr = 0usize;
1028 let mut sums = 0.0;
1029 let mut sumh = 0.0;
1030 let mut sumsed = 0.0;
1031 let mut sumsed_land = 0.0;
1032 let mut ntherm = 0usize;
1033 // ln of product of actual slopes (only where actual is above critical).
1034 let mut prods_therm = 0.0;
1035 // ln of product of critical slopes (only where actual is above critical).
1036 let mut prodscrit_therm = 0.0;
1037 let avgz = |x, y: usize| if y == 0 { f64::INFINITY } else { x / y as f64 };
1038 let geomz = |x: f64, y: usize| {
1039 if y == 0 {
1040 f64::INFINITY
1041 } else {
1042 (x / y as f64).exp()
1043 }
1044 };
1045
1046 // Gauss-Seidel iteration
1047
1048 let mut lake_silt: Box<[Compute]> = vec![0.0_f64; map_size_lg.chunks_len()].into_boxed_slice();
1049 let mut lake_sill = vec![-1isize; map_size_lg.chunks_len()].into_boxed_slice();
1050
1051 let mut n_gs_stream_power_law = 0;
1052 // NOTE: Increasing this can theoretically sometimes be necessary for
1053 // convergence, but in practice it is fairly unlikely that you should need
1054 // to do this (especially if you stick to g ∈ [0, 1]).
1055 let max_n_gs_stream_power_law = 99;
1056 let mut mstack_inv = vec![0usize; dh.len()];
1057 let mut h_t_stack = vec![Zero::zero(); dh.len()];
1058 let mut dh_stack = vec![0isize; dh.len()];
1059 let mut h_stack = vec![Zero::zero(); dh.len()];
1060 let mut b_stack = vec![Zero::zero(); dh.len()];
1061 let mut area_stack = vec![Zero::zero(); dh.len()];
1062 assert_eq!(mstack.len(), dh.len());
1063 assert_eq!(b.len(), dh.len());
1064 assert_eq!(h_t.len(), dh.len());
1065 let mstack_inv = &mut *mstack_inv;
1066 mstack.iter().enumerate().for_each(|(stacki, &posi)| {
1067 let posi = posi as usize;
1068 mstack_inv[posi] = stacki;
1069 dh_stack[stacki] = dh[posi];
1070 h_t_stack[stacki] = h_t[posi];
1071 h_stack[stacki] = h[posi];
1072 b_stack[stacki] = b[posi];
1073 area_stack[stacki] = area[posi];
1074 });
1075
1076 while err > tol && n_gs_stream_power_law < max_n_gs_stream_power_law {
1077 debug!("Stream power iteration #{:?}", n_gs_stream_power_law);
1078
1079 // Reset statistics in each loop.
1080 maxh = 0.0;
1081 minh = <Alt as Float>::infinity();
1082 nland = 0usize;
1083 ncorr = 0usize;
1084 sums = 0.0;
1085 sumh = 0.0;
1086 sumsed = 0.0;
1087 sumsed_land = 0.0;
1088 ntherm = 0usize;
1089 prods_therm = 0.0;
1090 prodscrit_therm = 0.0;
1091
1092 let start_time = Instant::now();
1093 // Keep track of how many iterations we've gone to test for convergence.
1094 n_gs_stream_power_law += 1;
1095
1096 threadpool.join(
1097 || {
1098 // guess/update the elevation at t+Δt (k)
1099 (&mut *h_p, &*h_stack)
1100 .into_par_iter()
1101 .for_each(|(h_p, h_)| {
1102 *h_p = (*h_) as Compute;
1103 });
1104 },
1105 || {
1106 // calculate erosion/deposition of sediment at each node
1107 (&*mstack, &mut *deltah, &*h_t_stack, &*h_stack)
1108 .into_par_iter()
1109 .for_each(|(&posi, deltah, &h_t_i, &h_i)| {
1110 let posi = posi as usize;
1111 let uplift_i = uplift(posi) as Alt;
1112 let delta = (h_t_i + uplift_i - h_i) as Compute;
1113 *deltah = delta;
1114 });
1115 },
1116 );
1117 debug!(
1118 "(Done pre-computation, time={:?}ms).",
1119 start_time.elapsed().as_millis()
1120 );
1121 #[rustfmt::skip]
1122 // sum the erosion in stack order
1123 //
1124 // After:
1125 // deltah_i = Σ{j ∈ {i} ∪ upstream_i(t)}(h_j(t, FINAL) + U_j * Δt - h_i(t + Δt, k))
1126 let start_time = Instant::now();
1127 izip!(&*mstack, &*dh_stack, &h_t_stack, &*h_p)
1128 .enumerate()
1129 .for_each(|(stacki, (&posi, &posj, &h_t_i, &h_p_i))| {
1130 let posi = posi as usize;
1131 let deltah_i = deltah[stacki];
1132 if posj < 0 {
1133 lake_silt[stacki] = deltah_i;
1134 } else {
1135 let uplift_i = uplift(posi) as Alt;
1136 let uphill_deltah_i = deltah_i - ((h_t_i + uplift_i) as Compute - h_p_i);
1137 let lposi = lake_sill[stacki];
1138 if lposi == stacki as isize {
1139 if uphill_deltah_i <= 0.0 {
1140 lake_silt[stacki] = 0.0;
1141 } else {
1142 lake_silt[stacki] = uphill_deltah_i;
1143 }
1144 }
1145 let mwrec_i = &mwrec[posi];
1146 mrec_downhill(map_size_lg, &mrec, posi).for_each(|(k, posj)| {
1147 let stack_posj = mstack_inv[posj];
1148 deltah[stack_posj] += deltah_i * mwrec_i[k];
1149 });
1150 }
1151 });
1152 debug!(
1153 "(Done sediment transport computation, time={:?}ms).",
1154 start_time.elapsed().as_millis()
1155 );
1156 // ```ignore
1157 // do ij=nn,1,-1
1158 // ijk=stack(ij)
1159 // ijr=rec(ijk)
1160 // if (ijr.ne.ijk) then
1161 // dh(ijk)=dh(ijk)-(ht(ijk)-hp(ijk))
1162 // if (lake_sill(ijk).eq.ijk) then
1163 // if (dh(ijk).le.0.d0) then
1164 // lake_sediment(ijk)=0.d0
1165 // else
1166 // lake_sediment(ijk)=dh(ijk)
1167 // endif
1168 // endif
1169 // dh(ijk)=dh(ijk)+(ht(ijk)-hp(ijk))
1170 // dh(ijr)=dh(ijr)+dh(ijk)
1171 // else
1172 // lake_sediment(ijk)=dh(ijk)
1173 // endif
1174 // enddo
1175 // ```
1176
1177 let start_time = Instant::now();
1178 (
1179 &*mstack,
1180 &mut *elev,
1181 &*dh_stack,
1182 &*h_t_stack,
1183 &*area_stack,
1184 &*deltah,
1185 &*h_p,
1186 &*b_stack,
1187 )
1188 .into_par_iter()
1189 .for_each(
1190 |(&posi, elev, &dh_i, &h_t_i, &area_i, &deltah_i, &h_p_i, &b_i)| {
1191 let posi = posi as usize;
1192
1193 let uplift_i = uplift(posi) as Alt;
1194 if dh_i < 0 {
1195 *elev = (h_t_i + uplift_i) as Compute;
1196 } else {
1197 let old_h_after_uplift_i = (h_t_i + uplift_i) as Compute;
1198 let area_i = area_i as Compute;
1199 let uphill_silt_i = deltah_i - (old_h_after_uplift_i - h_p_i);
1200 let old_b_i = b_i;
1201 let sed = h_t_i - old_b_i;
1202 let n = n_f(posi);
1203 let g_i = if sed > sediment_thickness(n) {
1204 (g_fs_mult_sed * g(posi)) as Compute
1205 } else {
1206 g(posi) as Compute
1207 };
1208 // Make sure deposition coefficient doesn't result in more deposition than
1209 // there actually was material to deposit. The
1210 // current assumption is that as long as we
1211 // are storing at most as much sediment as there actually was along the
1212 // river, we are in the clear.
1213 let g_i_ratio = g_i / (p * area_i);
1214 // One side of nonlinear equation (23):
1215 //
1216 // h_i(t) + U_i * Δt + G / (p̃ * Ã_i) * Σ
1217 // {j ∈ upstream_i(t)}(h_j(t, FINAL)
1218 // + U_j * Δt - h_j(t + Δt, k))
1219 //
1220 // where
1221 //
1222 // Ã_i = A_i / (∆x∆y) = N_i,
1223 // number of cells upstream of cell i.
1224 *elev = old_h_after_uplift_i + uphill_silt_i * g_i_ratio;
1225 }
1226 },
1227 );
1228 debug!(
1229 "(Done elevation estimation, time={:?}ms).",
1230 start_time.elapsed().as_millis()
1231 );
1232
1233 let start_time = Instant::now();
1234 // TODO: Consider taking advantage of multi-receiver flow here.
1235 // Iterate in ascending height order.
1236 let mut sum_err: Compute = 0.0_f64;
1237 izip!(&*mstack, &*elev, &*b_stack, &*h_t_stack, &*dh_stack, &*h_p)
1238 .enumerate()
1239 .rev()
1240 .for_each(|(stacki, (&posi, &elev_i, &b_i, &h_t_i, &dh_i, &h_p_i))| {
1241 let iteration_error = 0.0;
1242 let posi = posi as usize;
1243 let old_elev_i = elev_i;
1244 let old_b_i = b_i;
1245 let old_ht_i = h_t_i;
1246 let sed = old_ht_i - old_b_i;
1247
1248 let posj = dh_i;
1249 if posj < 0 {
1250 if posj == -1 {
1251 panic!("Disconnected lake!");
1252 }
1253 if h_t_i > 0.0 {
1254 warn!("Ocean above zero?");
1255 }
1256 // Egress with no outgoing flows.
1257 // wh for oceans is always at least min_erosion_height.
1258 let uplift_i = uplift(posi) as Alt;
1259 wh[posi] = FloatCore::max(min_erosion_height, h_t_i + uplift_i);
1260 lake_sill[stacki] = posi as isize;
1261 lake_water_volume[stacki] = 0.0;
1262 } else {
1263 let posj = posj as usize;
1264
1265 // Has an outgoing flow edge (posi, posj).
1266 // flux(i) = k * A[i]^m / ((p(i) - p(j)).magnitude()), and δt = 1
1267 // h[i](t + dt) = (h[i](t) + δt * (uplift[i] + flux(i) * h[j](t + δt))) / (1 +
1268 // flux(i) * δt). NOTE: posj has already been computed since
1269 // it's downhill from us. Therefore, we can rely on wh being
1270 // set to the water height for that node.
1271 // let h_j = h[posj_stack] as f64;
1272 let wh_j = wh[posj];
1273 let old_h_i = h_stack[stacki];
1274 let mut new_h_i = old_h_i;
1275
1276 // Only perform erosion if we are above the water level of the previous node.
1277 // NOTE: Can replace wh_j with h_j here (and a few other places) to allow
1278 // erosion underwater, producing very different looking
1279 // maps!
1280 if old_elev_i > wh_j
1281 /* h_j */
1282 {
1283 let dtherm = 0.0f64;
1284 let h0 = old_elev_i + dtherm;
1285
1286 // hi(t + ∂t) = (hi(t) + ∂t(ui + kp^mAi^m(hj(t + ∂t)/||pi - pj||))) / (1 +
1287 // ∂t * kp^mAi^m / ||pi - pj||)
1288 let n = n_f(posi) as f64;
1289
1290 // Fluvial erosion.
1291 let k_df_fact = &k_df_fact[posi];
1292 let k_fs_fact = &k_fs_fact[posi];
1293 if (n - 1.0).abs() <= 1.0e-3 && (q_ - 1.0).abs() <= 1.0e-3 {
1294 let mut f = h0;
1295 let mut df = 1.0;
1296 mrec_downhill(map_size_lg, &mrec, posi).for_each(|(kk, posj)| {
1297 let posj_stack = mstack_inv[posj];
1298 let h_j = h_stack[posj_stack];
1299 // This can happen in cases where receiver kk is neither uphill of
1300 // nor downhill from posi's direct receiver.
1301 // NOTE: Fastscape does h_t[posi] + uplift(posi) as f64 >= h_t[posj]
1302 // + uplift(posj) as f64
1303 // NOTE: We also considered using old_elev_i > wh[posj] here.
1304 if old_elev_i > h_j {
1305 let elev_j = h_j;
1306 let fact = k_fs_fact[kk] as f64 + k_df_fact[kk] as f64;
1307 f += fact * elev_j;
1308 df += fact;
1309 }
1310 });
1311 new_h_i = f / df;
1312 } else {
1313 // Local Newton-Raphson
1314 // TODO: Work out how (if possible) to make this converge for tiny n.
1315 let omega1 = 0.875f64 * n;
1316 let omega2 = 0.875f64 / q_;
1317 let omega = omega1.max(omega2);
1318 let tolp = 1.0e-3;
1319 let mut errp = 2.0 * tolp;
1320 let mut rec_heights = [0.0; 8];
1321 let mut mask = [MaskType::splat(false); 8];
1322 mrec_downhill(map_size_lg, &mrec, posi).for_each(|(kk, posj)| {
1323 let posj_stack = mstack_inv[posj];
1324 let h_j = h_stack[posj_stack];
1325 // NOTE: Fastscape does h_t[posi] + uplift(posi) as f64 >= h_t[posj]
1326 // + uplift(posj) as f64
1327 // NOTE: We also considered using old_elev_i > wh[posj] here.
1328 if old_elev_i > h_j {
1329 mask[kk] = MaskType::splat(true);
1330 rec_heights[kk] = h_j as SimdType;
1331 }
1332 });
1333 while errp > tolp {
1334 let mut f = new_h_i - h0;
1335 let mut df = 1.0;
1336 izip!(&mask, &rec_heights, k_fs_fact, k_df_fact).for_each(
1337 |(&mask_kk, &rec_heights_kk, &k_fs_fact_kk, &k_df_fact_kk)| {
1338 if mask_kk.any() {
1339 let h_j = rec_heights_kk;
1340 let elev_j = h_j;
1341 let dh =
1342 FloatCore::max(0.0, new_h_i as SimdType - elev_j);
1343 let powf = |a: SimdType, b| a.powf(b);
1344 let dh_fs_sample = k_fs_fact_kk as SimdType
1345 * powf(dh, n as SimdType - 1.0);
1346 let dh_df_sample = k_df_fact_kk as SimdType
1347 * powf(dh, q_ as SimdType - 1.0);
1348 // Want: h_i(t+Δt) = h0 - fact * (h_i(t+Δt) -
1349 // h_j(t+Δt))^n
1350 // Goal: h_i(t+Δt) - h0 + fact * (h_i(t+Δt) -
1351 // h_j(t+Δt))^n = 0
1352 f += ((dh_fs_sample + dh_df_sample) * dh) as f64;
1353 // ∂h_i(t+Δt)/∂n = 1 + fact * n * (h_i(t+Δt) -
1354 // h_j(t+Δt))^(n - 1)
1355 df += (n as SimdType * dh_fs_sample
1356 + q_ as SimdType * dh_df_sample)
1357 as f64;
1358 }
1359 },
1360 );
1361 // hn = h_i(t+Δt, k) - (h_i(t+Δt, k) - (h0 - fact * (h_i(t+Δt, k) -
1362 // h_j(t+Δt))^n)) / ∂h_i/∂n(t+Δt, k)
1363 let hn = new_h_i - f / df;
1364 // errp = |(h_i(t+Δt, k) - (h0 - fact * (h_i(t+Δt, k) -
1365 // h_j(t+Δt))^n)) / ∂h_i/∂n(t+Δt, k)|
1366 errp = (hn - new_h_i).abs();
1367 // h_i(t+∆t, k+1) = ...
1368 new_h_i = new_h_i * (1.0 - omega) + hn * omega;
1369 }
1370 /* omega=0.875d0/n
1371 tolp=1.d-3
1372 errp=2.d0*tolp
1373 h0=elev(ijk)
1374 do while (errp.gt.tolp)
1375 f=h(ijk)-h0
1376 df=1.d0
1377 if (ht(ijk).gt.ht(ijr)) then
1378 fact = kfint(ijk)*dt*a(ijk)**m/length(ijk)**n
1379 f=f+fact*max(0.d0,h(ijk)-h(ijr))**n
1380 df=df+fact*n*max(0.d0,h(ijk)-h(ijr))**(n-1.d0)
1381 endif
1382 hn=h(ijk)-f/df
1383 errp=abs(hn-h(ijk))
1384 h(ijk)=h(ijk)*(1.d0-omega)+hn*omega
1385 enddo */
1386 }
1387
1388 lake_sill[stacki] = posi as isize;
1389 lake_water_volume[stacki] = 0.0;
1390
1391 // If we dipped below the receiver's water level, set our height to the
1392 // receiver's water level.
1393 // NOTE: If we want erosion to proceed underwater, use h_j here instead of
1394 // wh_j.
1395 if new_h_i <= wh_j
1396 /* h_j */
1397 {
1398 if compute_stats {
1399 ncorr += 1;
1400 }
1401 // NOTE: Why wh_j?
1402 // Because in the next round, if the old height is still wh_j or under,
1403 // it will be set precisely equal to the
1404 // estimated height, meaning it effectively
1405 // "vanishes" and just deposits sediment to its receiver.
1406 // (This is probably related to criteria for block Gauss-Seidel, etc.).
1407 // NOTE: If we want erosion to proceed underwater, use h_j here instead
1408 // of wh_j.
1409 new_h_i = wh_j;
1410 } else if compute_stats && new_h_i > 0.0 {
1411 let dxy = (uniform_idx_as_vec2(map_size_lg, posi)
1412 - uniform_idx_as_vec2(map_size_lg, posj))
1413 .map(|e| e as f64);
1414 let neighbor_distance = (neighbor_coef * dxy).magnitude();
1415 let dz = (new_h_i - wh_j).max(0.0);
1416 let mag_slope = dz / neighbor_distance;
1417
1418 nland += 1;
1419 sumsed_land += sed;
1420 sumh += new_h_i;
1421 sums += mag_slope;
1422 }
1423 } else {
1424 new_h_i = old_elev_i;
1425 let posj_stack = mstack_inv[posj];
1426 let lposj = lake_sill[posj_stack];
1427 lake_sill[stacki] = lposj;
1428 if lposj >= 0 {
1429 let lposj = lposj as usize;
1430 lake_water_volume[lposj] += (wh_j - old_elev_i) as Compute;
1431 }
1432 }
1433 // Set max_slope to this node's water height (max of receiver's water height and
1434 // this node's height).
1435 wh[posi] = wh_j.max(new_h_i) as Alt;
1436 h_stack[stacki] = new_h_i as Alt;
1437 }
1438 if compute_stats {
1439 sumsed += sed;
1440 let h_i = h_stack[stacki];
1441 if h_i > 0.0 {
1442 minh = h_i.min(minh);
1443 }
1444 maxh = h_i.max(maxh);
1445 }
1446
1447 // Add sum of squares of errors.
1448 sum_err +=
1449 (iteration_error + h_stack[stacki] as Compute - h_p_i as Compute).powi(2);
1450 });
1451 debug!(
1452 "(Done erosion computation, time={:?}ms)",
1453 start_time.elapsed().as_millis()
1454 );
1455
1456 err = (sum_err / mstack.len() as Compute).sqrt();
1457 debug!("(RMSE: {:?})", err);
1458 if max_g == 0.0 {
1459 err = 0.0;
1460 }
1461 if n_gs_stream_power_law == max_n_gs_stream_power_law {
1462 warn!(
1463 "Beware: Gauss-Seidel scheme not convergent: err={:?}, expected={:?}",
1464 err, tol
1465 );
1466 }
1467 }
1468
1469 (&*mstack_inv, &mut *h)
1470 .into_par_iter()
1471 .enumerate()
1472 .for_each(|(posi, (&stacki, h))| {
1473 assert_eq!(posi, mstack[stacki] as usize);
1474 *h = h_stack[stacki];
1475 });
1476
1477 // update the basement
1478 //
1479 // NOTE: Despite this not quite applying since basement order and height order
1480 // differ, we once again borrow the implicit FastScape stack order. If this
1481 // becomes a problem we can easily compute a separate stack order just for
1482 // basement. TODO: Consider taking advantage of multi-receiver flow here.
1483 b.par_iter_mut()
1484 .zip_eq(h.par_iter())
1485 .enumerate()
1486 .for_each(|(posi, (b, &h_i))| {
1487 let old_b_i = *b;
1488 let uplift_i = uplift(posi) as Alt;
1489
1490 // First, add uplift...
1491 let mut new_b_i = old_b_i + uplift_i;
1492
1493 let posj = dh[posi];
1494 // Sediment height normal to bedrock. NOTE: Currently we can actually have
1495 // sediment and bedrock slope at different heights, meaning there's
1496 // no uniform slope. There are probably more correct ways to
1497 // account for this, such as averaging, integrating, or doing things
1498 // by mass / volume instead of height, but for now we use the time-honored
1499 // technique of ignoring the problem.
1500 let vertical_sed = (h_i - new_b_i).max(0.0);
1501 let h_normal = if posj < 0 {
1502 // Egress with no outgoing flows; for now, we assume this means normal and
1503 // vertical coincide.
1504 vertical_sed
1505 } else {
1506 let posj = posj as usize;
1507 let h_j = h[posj];
1508 let dxy = (uniform_idx_as_vec2(map_size_lg, posi)
1509 - uniform_idx_as_vec2(map_size_lg, posj))
1510 .map(|e| e as f64);
1511 let neighbor_distance_squared = (neighbor_coef * dxy).magnitude_squared();
1512 let dh = h_i - h_j;
1513 // H_i_fact = (h_i - h_j) / (||p_i - p_j||^2 + (h_i - h_j)^2)
1514 let h_i_fact = dh / (neighbor_distance_squared + dh * dh);
1515 let h_i_vertical = 1.0 - h_i_fact * dh;
1516 // ||H_i|| = (h_i - b_i) * √((H_i_fact^2 * ||p_i - p_j||^2 + (1 - H_i_fact *
1517 // (h_i - h_j))^2))
1518 vertical_sed
1519 * (h_i_fact * h_i_fact * neighbor_distance_squared
1520 + h_i_vertical * h_i_vertical)
1521 .sqrt()
1522 };
1523 // Rate of sediment production: -∂z_b / ∂t = ε₀ * e^(-αH)
1524 let p_i = epsilon_0(posi) as f64 * dt * f64::exp(-alpha(posi) as f64 * h_normal);
1525
1526 new_b_i -= p_i as Alt;
1527
1528 // Clamp basement so it doesn't exceed height.
1529 new_b_i = new_b_i.min(h_i);
1530 *b = new_b_i;
1531 });
1532 debug!("Done updating basement and applying soil production...");
1533
1534 // update the height to reflect sediment flux.
1535 if max_g > 0.0 {
1536 // If max_g = 0.0, lake_silt will be too high during the first iteration since
1537 // our initial estimate for h is very poor; however, the elevation
1538 // estimate will have been unaffected by g.
1539 (&mut *h, &*mstack_inv)
1540 .into_par_iter()
1541 .enumerate()
1542 .for_each(|(posi, (h, &stacki))| {
1543 let lposi = lake_sill[stacki];
1544 if lposi >= 0 {
1545 let lposi = lposi as usize;
1546 if lake_water_volume[lposi] > 0.0 {
1547 // +max(0.d0,min(lake_sediment(lake_sill(ij)),
1548 // lake_water_volume(lake_sill(ij))))/
1549 // lake_water_volume(lake_sill(ij))*(water(ij)-h(ij))
1550 *h += (FloatCore::max(0.0, lake_silt[stacki].min(lake_water_volume[lposi]))
1551 / lake_water_volume[lposi]
1552 * (wh[posi] - *h) as Compute) as Alt;
1553 }
1554 }
1555 });
1556 }
1557 // do ij=1,nn
1558 // if (lake_sill(ij).ne.0) then
1559 // if (lake_water_volume(lake_sill(ij)).gt.0.d0) h(ij)=h(ij) &
1560 // +max(0.d0,min(lake_sediment(lake_sill(ij)),
1561 // lake_water_volume(lake_sill(ij))))/ &
1562 // lake_water_volume(lake_sill(ij))*(water(ij)-h(ij))
1563 // endif
1564 // enddo
1565
1566 debug!(
1567 "Done applying stream power (max height: {:?}) (avg height: {:?}) (min height: {:?}) (avg \
1568 slope: {:?})\n (above talus angle, geom. mean slope [actual/critical/ratio]: {:?} \
1569 / {:?} / {:?})\n (old avg sediment thickness [all/land]: {:?} / {:?})\n \
1570 (num land: {:?}) (num thermal: {:?}) (num corrected: {:?})",
1571 maxh,
1572 avgz(sumh, nland),
1573 minh,
1574 avgz(sums, nland),
1575 geomz(prods_therm, ntherm),
1576 geomz(prodscrit_therm, ntherm),
1577 geomz(prods_therm - prodscrit_therm, ntherm),
1578 avgz(sumsed, newh.len()),
1579 avgz(sumsed_land, nland),
1580 nland,
1581 ntherm,
1582 ncorr,
1583 );
1584
1585 // Apply thermal erosion.
1586 maxh = 0.0;
1587 minh = <Alt as Float>::infinity();
1588 sumh = 0.0;
1589 sums = 0.0;
1590 sumsed = 0.0;
1591 sumsed_land = 0.0;
1592 nland = 0usize;
1593 ncorr = 0usize;
1594 ntherm = 0usize;
1595 prods_therm = 0.0;
1596 prodscrit_therm = 0.0;
1597 newh.iter().for_each(|&posi| {
1598 let posi = posi as usize;
1599 let old_h_i = h[posi];
1600 let old_b_i = b[posi];
1601 let sed = old_h_i - old_b_i;
1602
1603 let max_slope = max_slopes[posi];
1604 let n = n_f(posi);
1605 max_slopes[posi] = if sed > sediment_thickness(n) && kdsed > 0.0 {
1606 // Sediment
1607 kdsed
1608 } else {
1609 // Bedrock
1610 kd(posi)
1611 };
1612
1613 let posj = dh[posi];
1614 if posj < 0 {
1615 // Egress with no outgoing flows.
1616 if posj == -1 {
1617 panic!("Disconnected lake!");
1618 }
1619 // wh for oceans is always at least min_erosion_height.
1620 wh[posi] = FloatCore::max(min_erosion_height, old_h_i as Alt);
1621 } else {
1622 let posj = posj as usize;
1623 // Find the water height for this chunk's receiver; we only apply thermal
1624 // erosion for chunks above water.
1625 let wh_j = wh[posj];
1626 let mut new_h_i = old_h_i;
1627 if wh_j < old_h_i {
1628 // NOTE: Currently assuming that talus angle is not eroded once the substance is
1629 // totally submerged in water, and that talus angle if part of the substance is
1630 // in water is 0 (or the same as the dry part, if this is set to wh_j), but
1631 // actually that's probably not true.
1632 let old_h_j = h[posj];
1633 let h_j = old_h_j;
1634 // Test the slope.
1635 // Hacky version of thermal erosion: only consider lowest neighbor, don't
1636 // redistribute uplift to other neighbors.
1637 let (posk, h_k) = (posj, h_j);
1638 let (posk, h_k) = if h_k < h_j { (posk, h_k) } else { (posj, h_j) };
1639 let dxy = (uniform_idx_as_vec2(map_size_lg, posi)
1640 - uniform_idx_as_vec2(map_size_lg, posk))
1641 .map(|e| e as f64);
1642 let neighbor_distance = (neighbor_coef * dxy).magnitude();
1643 let dz = (new_h_i - h_k).max(0.0);
1644 let mag_slope = dz / neighbor_distance;
1645 if mag_slope >= max_slope {
1646 let dtherm = 0.0;
1647 new_h_i -= dtherm;
1648 if new_h_i <= wh_j {
1649 if compute_stats {
1650 ncorr += 1;
1651 }
1652 } else if compute_stats && new_h_i > 0.0 {
1653 let dz = (new_h_i - h_j).max(0.0);
1654 let slope = dz / neighbor_distance;
1655 sums += slope;
1656 nland += 1;
1657 sumh += new_h_i;
1658 sumsed_land += sed;
1659 }
1660 if compute_stats {
1661 ntherm += 1;
1662 prodscrit_therm += max_slope.ln();
1663 prods_therm += mag_slope.ln();
1664 }
1665 } else {
1666 // Poorly emulating nonlinear hillslope transport as described by
1667 // http://eps.berkeley.edu/~bill/papers/112.pdf.
1668 // sqrt(3)/3*32*32/(128000/2)
1669 // Also Perron-2011-Journal_of_Geophysical_Research__Earth_Surface.pdf
1670 let slope_ratio = (mag_slope / max_slope).powi(2);
1671 let slope_nonlinear_factor =
1672 slope_ratio * ((3.0 - slope_ratio) / (1.0 - slope_ratio).powi(2));
1673 max_slopes[posi] += (max_slopes[posi] * slope_nonlinear_factor).min(max_stable);
1674 if compute_stats && new_h_i > 0.0 {
1675 sums += mag_slope;
1676 // Just use the computed rate.
1677 nland += 1;
1678 sumh += new_h_i;
1679 sumsed_land += sed;
1680 }
1681 }
1682 }
1683 // Set wh to this node's water height (max of receiver's water height and
1684 // this node's height).
1685 wh[posi] = wh_j.max(new_h_i) as Alt;
1686 }
1687
1688 if compute_stats {
1689 sumsed += sed;
1690 let h_i = h[posi];
1691 if h_i > 0.0 {
1692 minh = h_i.min(minh);
1693 }
1694 maxh = h_i.max(maxh);
1695 }
1696 });
1697 debug!(
1698 "Done applying thermal erosion (max height: {:?}) (avg height: {:?}) (min height: {:?}) \
1699 (avg slope: {:?})\n (above talus angle, geom. mean slope [actual/critical/ratio]: \
1700 {:?} / {:?} / {:?})\n (avg sediment thickness [all/land]: {:?} / {:?})\n \
1701 (num land: {:?}) (num thermal: {:?}) (num corrected: {:?})",
1702 maxh,
1703 avgz(sumh, nland),
1704 minh,
1705 avgz(sums, nland),
1706 geomz(prods_therm, ntherm),
1707 geomz(prodscrit_therm, ntherm),
1708 geomz(prods_therm - prodscrit_therm, ntherm),
1709 avgz(sumsed, newh.len()),
1710 avgz(sumsed_land, nland),
1711 nland,
1712 ntherm,
1713 ncorr,
1714 );
1715
1716 // Apply hillslope diffusion.
1717 diffusion(
1718 nx,
1719 ny,
1720 nx as f64 * dx,
1721 ny as f64 * dy,
1722 dt,
1723 (),
1724 h,
1725 b,
1726 |posi| max_slopes[posi],
1727 -1.0,
1728 );
1729 debug!("Done applying diffusion.");
1730 debug!("Done eroding.");
1731}
1732
1733/// The Planchon-Darboux algorithm for extracting drainage networks.
1734///
1735/// http://horizon.documentation.ird.fr/exl-doc/pleins_textes/pleins_textes_7/sous_copyright/010031925.pdf
1736///
1737/// See https://github.com/mewo2/terrain/blob/master/terrain.js
1738pub(crate) fn fill_sinks<F: Float + Send + Sync>(
1739 map_size_lg: MapSizeLg,
1740 h: impl Fn(usize) -> F + Sync,
1741 is_ocean: impl Fn(usize) -> bool + Sync,
1742) -> Box<[F]> {
1743 // NOTE: We are using the "exact" version of depression-filling, which is slower
1744 // but doesn't change altitudes.
1745 let epsilon = F::zero();
1746 let infinity = F::infinity();
1747 let range = 0..map_size_lg.chunks_len();
1748 let oldh = range
1749 .into_par_iter()
1750 .map(&h)
1751 .collect::<Vec<_>>()
1752 .into_boxed_slice();
1753 let mut newh = oldh
1754 .par_iter()
1755 .enumerate()
1756 .map(|(posi, &h)| {
1757 let is_near_edge = is_ocean(posi);
1758 if is_near_edge {
1759 debug_assert!(h <= F::zero());
1760 h
1761 } else {
1762 infinity
1763 }
1764 })
1765 .collect::<Vec<_>>()
1766 .into_boxed_slice();
1767
1768 loop {
1769 let mut changed = false;
1770 (0..newh.len()).for_each(|posi| {
1771 let nh = newh[posi];
1772 let oh = oldh[posi];
1773 if nh == oh {
1774 return;
1775 }
1776 for nposi in neighbors(map_size_lg, posi) {
1777 let onbh = newh[nposi];
1778 let nbh = onbh + epsilon;
1779 // If there is even one path downhill from this node's original height, fix
1780 // the node's new height to be equal to its original height, and break out of
1781 // the loop.
1782 if oh >= nbh {
1783 newh[posi] = oh;
1784 changed = true;
1785 break;
1786 }
1787 // Otherwise, we know this node's original height is below the new height of all
1788 // of its neighbors. Then, we try to choose the minimum new
1789 // height among all this node's neighbors that is (plus a
1790 // constant epsilon) below this node's new height.
1791 //
1792 // (If there is no such node, then the node's new height must be (minus a
1793 // constant epsilon) lower than the new height of every
1794 // neighbor, but above its original height. But this can't be
1795 // true for *all* nodes, because if this is true for any
1796 // node, it is not true of any of its neighbors. So all neighbors must either
1797 // be their original heights, or we will have another iteration
1798 // of the loop (one of its neighbors was changed to its minimum
1799 // neighbor). In the second case, in the next round, all
1800 // neighbor heights will be at most nh + epsilon).
1801 if nh > nbh && nbh > oh {
1802 newh[posi] = nbh;
1803 changed = true;
1804 }
1805 }
1806 });
1807 if !changed {
1808 return newh;
1809 }
1810 }
1811}
1812
1813/// Algorithm for finding and connecting lakes. Assumes newh and downhill have
1814/// already been computed. When a lake's value is negative, it is its own lake
1815/// root, and when it is 0, it is on the boundary of Ω.
1816///
1817/// Returns a 4-tuple containing:
1818/// - The first indirection vector (associating chunk indices with their lake's
1819/// root node).
1820/// - A list of chunks on the boundary (non-lake egress points).
1821/// - The second indirection vector (associating chunk indices with their lake's
1822/// adjacency list).
1823/// - The adjacency list (stored in a single vector), indexed by the second
1824/// indirection vector.
1825pub fn get_lakes<F: FloatCore>(
1826 map_size_lg: MapSizeLg,
1827 h: impl Fn(usize) -> F,
1828 downhill: &mut [isize],
1829) -> (usize, Box<[i32]>, Box<[u32]>, F) {
1830 // Associates each lake index with its root node (the deepest one in the lake),
1831 // and a list of adjacent lakes. The list of adjacent lakes includes the
1832 // lake index of the adjacent lake, and a node index in the adjacent lake
1833 // which has a neighbor in this lake. The particular neighbor should be the
1834 // one that generates the minimum "pass height" encountered so far, i.e. the
1835 // chosen pair should minimize the maximum of the heights of the nodes in the
1836 // pair.
1837
1838 // We start by taking steps to allocate an indirection vector to use for storing
1839 // lake indices. Initially, each entry in this vector will contain 0. We
1840 // iterate in ascending order through the sorted newh. If the node has
1841 // downhill == -2, it is a boundary node Ω and we store it in the boundary
1842 // vector. If the node has downhill == -1, it is a fresh lake, and we store 0
1843 // in it. If the node has non-negative downhill, we use the downhill index
1844 // to find the next node down; if the downhill node has a lake entry < 0,
1845 // then downhill is a lake and its entry can be negated to find an
1846 // (over)estimate of the number of entries it needs. If the downhill
1847 // node has a non-negative entry, then its entry is the lake index for this
1848 // node, so we should access that entry and increment it, then set our own
1849 // entry to it.
1850 let mut boundary = Vec::with_capacity(downhill.len());
1851 let mut indirection = vec![/*-1i32*/0i32; map_size_lg.chunks_len()].into_boxed_slice();
1852
1853 let mut newh = Vec::with_capacity(downhill.len());
1854
1855 // Now, we know that the sum of all the indirection nodes will be the same as
1856 // the number of nodes. We can allocate a *single* vector with 8 * nodes
1857 // entries, to be used as storage space, and augment our indirection vector
1858 // with the starting index, resulting in a vector of slices. As we go, we
1859 // replace each lake entry with its index in the new indirection buffer,
1860 // allowing
1861 let mut lakes = vec![(-1, 0); /*(indirection.len() - boundary.len())*/indirection.len() * 8];
1862 let mut indirection_ = vec![0u32; indirection.len()];
1863 // First, find all the lakes.
1864 let mut lake_roots = Vec::with_capacity(downhill.len()); // Test
1865 (*downhill)
1866 .iter()
1867 .enumerate()
1868 .filter(|(_, &dh_idx)| dh_idx < 0)
1869 .for_each(|(chunk_idx, &dh)| {
1870 if dh == -2 {
1871 // On the boundary, add to the boundary vector.
1872 boundary.push(chunk_idx);
1873 // Still considered a lake root, though.
1874 } else if dh == -1 {
1875 lake_roots.push(chunk_idx);
1876 } else {
1877 panic!("Impossible.");
1878 }
1879 // Find all the nodes uphill from this lake. Since there is only one outgoing
1880 // edge in the "downhill" graph, this is guaranteed never to visit a
1881 // node more than once.
1882 let start = newh.len();
1883 let indirection_idx = (start * 8) as u32;
1884 // New lake root
1885 newh.push(chunk_idx as u32);
1886 let mut cur = start;
1887 while cur < newh.len() {
1888 let node = newh[cur];
1889 uphill(map_size_lg, downhill, node as usize).for_each(|child| {
1890 // lake_idx is the index of our lake root.
1891 indirection[child] = chunk_idx as i32;
1892 indirection_[child] = indirection_idx;
1893 newh.push(child as u32);
1894 });
1895 cur += 1;
1896 }
1897 // Find the number of elements pushed.
1898 let length = (cur - start) * 8;
1899 // New lake root (lakes have indirection set to -length).
1900 indirection[chunk_idx] = -(length as i32);
1901 indirection_[chunk_idx] = indirection_idx;
1902 });
1903 assert_eq!(newh.len(), downhill.len());
1904
1905 debug!("Old lake roots: {:?}", lake_roots.len());
1906
1907 let newh = newh.into_boxed_slice();
1908 let mut maxh = -F::infinity();
1909 // Now, we know that the sum of all the indirection nodes will be the same as
1910 // the number of nodes. We can allocate a *single* vector with 8 * nodes
1911 // entries, to be used as storage space, and augment our indirection vector
1912 // with the starting index, resulting in a vector of slices. As we go, we
1913 // replace each lake entry with its index in the new indirection buffer,
1914 // allowing
1915 newh.iter().for_each(|&chunk_idx_| {
1916 let chunk_idx = chunk_idx_ as usize;
1917 let lake_idx_ = indirection_[chunk_idx];
1918 let lake_idx = lake_idx_ as usize;
1919 let height = h(chunk_idx_ as usize);
1920 // While we're here, compute the max elevation difference from zero among all
1921 // nodes.
1922 maxh = maxh.max(height.abs());
1923 // For every neighbor, check to see whether it is already set; if the neighbor
1924 // is set, its height is ≤ our height. We should search through the
1925 // edge list for the neighbor's lake to see if there's an entry; if not,
1926 // we insert, and otherwise we get its height. We do the same thing in
1927 // our own lake's entry list. If the maximum of the heights we get out
1928 // from this process is greater than the maximum of this chunk and its
1929 // neighbor chunk, we switch to this new edge.
1930 neighbors(map_size_lg, chunk_idx).for_each(|neighbor_idx| {
1931 let neighbor_height = h(neighbor_idx);
1932 let neighbor_lake_idx_ = indirection_[neighbor_idx];
1933 let neighbor_lake_idx = neighbor_lake_idx_ as usize;
1934 if neighbor_lake_idx_ < lake_idx_ {
1935 // We found an adjacent node that is not on the boundary and has already
1936 // been processed, and also has a non-matching lake. Therefore we can use
1937 // split_at_mut to get disjoint slices.
1938 let (lake, neighbor_lake) = {
1939 // println!("Okay, {:?} < {:?}", neighbor_lake_idx, lake_idx);
1940 let (neighbor_lake, lake) = lakes.split_at_mut(lake_idx);
1941 (lake, &mut neighbor_lake[neighbor_lake_idx..])
1942 };
1943
1944 // We don't actually need to know the real length here, because we've reserved
1945 // enough spaces that we should always either find a -1 (available slot) or an
1946 // entry for this chunk.
1947 'outer: for pass in lake.iter_mut() {
1948 if pass.0 == -1 {
1949 // println!("One time, in my mind, one time... (neighbor lake={:?}
1950 // lake={:?})", neighbor_lake_idx, lake_idx_);
1951 *pass = (chunk_idx_ as i32, neighbor_idx as u32);
1952 // Should never run out of -1s in the neighbor lake if we didn't find
1953 // the neighbor lake in our lake.
1954 *neighbor_lake
1955 .iter_mut()
1956 .find(|neighbor_pass| neighbor_pass.0 == -1)
1957 .unwrap() = (neighbor_idx as i32, chunk_idx_);
1958 // panic!("Should never happen; maybe didn't reserve enough space in
1959 // lakes?")
1960 break;
1961 } else if indirection_[pass.1 as usize] == neighbor_lake_idx_ {
1962 for neighbor_pass in neighbor_lake.iter_mut() {
1963 // Should never run into -1 while looping here, since (i, j)
1964 // and (j, i) should be added together.
1965 if indirection_[neighbor_pass.1 as usize] == lake_idx_ {
1966 let pass_height = h(neighbor_pass.1 as usize);
1967 let neighbor_pass_height = h(pass.1 as usize);
1968 if height.max(neighbor_height)
1969 < pass_height.max(neighbor_pass_height)
1970 {
1971 *pass = (chunk_idx_ as i32, neighbor_idx as u32);
1972 *neighbor_pass = (neighbor_idx as i32, chunk_idx_);
1973 }
1974 break 'outer;
1975 }
1976 }
1977 // Should always find a corresponding match in the neighbor lake if
1978 // we found the neighbor lake in our lake.
1979 let indirection_idx = indirection[chunk_idx];
1980 let lake_chunk_idx = if indirection_idx >= 0 {
1981 indirection_idx as usize
1982 } else {
1983 chunk_idx
1984 };
1985 let indirection_idx = indirection[neighbor_idx];
1986 let neighbor_lake_chunk_idx = if indirection_idx >= 0 {
1987 indirection_idx as usize
1988 } else {
1989 neighbor_idx
1990 };
1991 panic!(
1992 "For edge {:?} between lakes {:?}, couldn't find partner for pass \
1993 {:?}. Should never happen; maybe forgot to set both edges?",
1994 (
1995 (chunk_idx, uniform_idx_as_vec2(map_size_lg, chunk_idx)),
1996 (neighbor_idx, uniform_idx_as_vec2(map_size_lg, neighbor_idx))
1997 ),
1998 (
1999 (
2000 lake_chunk_idx,
2001 uniform_idx_as_vec2(map_size_lg, lake_chunk_idx),
2002 lake_idx_
2003 ),
2004 (
2005 neighbor_lake_chunk_idx,
2006 uniform_idx_as_vec2(map_size_lg, neighbor_lake_chunk_idx),
2007 neighbor_lake_idx_
2008 )
2009 ),
2010 (
2011 (pass.0, uniform_idx_as_vec2(map_size_lg, pass.0 as usize)),
2012 (pass.1, uniform_idx_as_vec2(map_size_lg, pass.1 as usize))
2013 ),
2014 );
2015 }
2016 }
2017 }
2018 });
2019 });
2020
2021 // Now it's time to calculate the lake connections graph T_L covering G_L.
2022 let mut candidates = BinaryHeap::with_capacity(indirection.len());
2023 // let mut pass_flows : Vec<i32> = vec![-1; indirection.len()];
2024
2025 // We start by going through each pass, deleting the ones that point out of
2026 // boundary nodes and adding ones that point into boundary nodes from
2027 // non-boundary nodes.
2028 lakes.iter_mut().for_each(|edge| {
2029 let edge: &mut (i32, u32) = edge;
2030 // Only consider valid elements.
2031 if edge.0 == -1 {
2032 return;
2033 }
2034 // Check to see if this edge points out *from* a boundary node.
2035 // Delete it if so.
2036 let from = edge.0 as usize;
2037 let indirection_idx = indirection[from];
2038 let lake_idx = if indirection_idx < 0 {
2039 from
2040 } else {
2041 indirection_idx as usize
2042 };
2043 if downhill[lake_idx] == -2 {
2044 edge.0 = -1;
2045 return;
2046 }
2047 // This edge is not pointing out from a boundary node.
2048 // Check to see if this edge points *to* a boundary node.
2049 // Add it to the candidate set if so.
2050 let to = edge.1 as usize;
2051 let indirection_idx = indirection[to];
2052 let lake_idx = if indirection_idx < 0 {
2053 to
2054 } else {
2055 indirection_idx as usize
2056 };
2057 if downhill[lake_idx] == -2 {
2058 // Find the pass height
2059 let pass = h(from).max(h(to));
2060 candidates.push(Reverse((
2061 NotNan::new(pass).unwrap(),
2062 (edge.0 as u32, edge.1),
2063 )));
2064 }
2065 });
2066
2067 let mut pass_flows_sorted: Vec<usize> = Vec::with_capacity(indirection.len());
2068
2069 // Now all passes pointing to the boundary are in candidates.
2070 // As long as there are still candidates, we continue...
2071 // NOTE: After a lake is added to the stream tree, the lake bottom's indirection
2072 // entry no longer negates its maximum number of passes, but the lake side
2073 // of the chosen pass. As such, we should make sure not to rely on using it
2074 // this way afterwards. provides information about the number of candidate
2075 // passes in a lake.
2076 while let Some(Reverse((_, (chunk_idx, neighbor_idx)))) = candidates.pop() {
2077 // We have the smallest candidate.
2078 let lake_idx = indirection_[chunk_idx as usize] as usize;
2079 let indirection_idx = indirection[chunk_idx as usize];
2080 let lake_chunk_idx = if indirection_idx >= 0 {
2081 indirection_idx as usize
2082 } else {
2083 chunk_idx as usize
2084 };
2085 if downhill[lake_chunk_idx] >= 0 {
2086 // Candidate lake has already been connected.
2087 continue;
2088 }
2089 assert_eq!(downhill[lake_chunk_idx], -1);
2090 // Candidate lake has not yet been connected, and is the lowest candidate.
2091 // Delete all other outgoing edges.
2092 let max_len = -if indirection_idx < 0 {
2093 indirection_idx
2094 } else {
2095 indirection[indirection_idx as usize]
2096 } as usize;
2097 // Add this chunk to the tree.
2098 downhill[lake_chunk_idx] = neighbor_idx as isize;
2099 // Also set the indirection of the lake bottom to the negation of the
2100 // lake side of the chosen pass (chunk_idx).
2101 // NOTE: This can't overflow i32 because map_size_lg.chunks_len() should fit
2102 // in an i32.
2103 indirection[lake_chunk_idx] = -(chunk_idx as i32);
2104 // Add this edge to the sorted list.
2105 pass_flows_sorted.push(lake_chunk_idx);
2106 // pass_flows_sorted.push((chunk_idx as u32, neighbor_idx as u32));
2107 for edge in &mut lakes[lake_idx..lake_idx + max_len] {
2108 if *edge == (chunk_idx as i32, neighbor_idx) {
2109 // Skip deleting this edge.
2110 continue;
2111 }
2112 // Delete the old edge, and remember it.
2113 let edge = mem::replace(edge, (-1, 0));
2114 if edge.0 == -1 {
2115 // Don't fall off the end of the list.
2116 break;
2117 }
2118 // Don't add incoming pointers from already-handled lakes or boundary nodes.
2119 let indirection_idx = indirection[edge.1 as usize];
2120 let neighbor_lake_idx = if indirection_idx < 0 {
2121 edge.1 as usize
2122 } else {
2123 indirection_idx as usize
2124 };
2125 if downhill[neighbor_lake_idx] != -1 {
2126 continue;
2127 }
2128 // Find the pass height
2129 let pass = h(edge.0 as usize).max(h(edge.1 as usize));
2130 // Put the reverse edge in candidates, sorted by height, then chunk idx, and
2131 // finally neighbor idx.
2132 candidates.push(Reverse((
2133 NotNan::new(pass).unwrap(),
2134 (edge.1, edge.0 as u32),
2135 )));
2136 }
2137 }
2138 debug!("Total lakes: {:?}", pass_flows_sorted.len());
2139
2140 // Perform the bfs once again.
2141 #[derive(Clone, Copy, PartialEq)]
2142 enum Tag {
2143 UnParsed,
2144 InQueue,
2145 WithRcv,
2146 }
2147 let mut tag = vec![Tag::UnParsed; map_size_lg.chunks_len()];
2148 // TODO: Combine with adding to vector.
2149 let mut filling_queue = Vec::with_capacity(downhill.len());
2150
2151 let mut newh = Vec::with_capacity(downhill.len());
2152 (*boundary)
2153 .iter()
2154 .chain(pass_flows_sorted.iter())
2155 .for_each(|&chunk_idx| {
2156 // Find all the nodes uphill from this lake. Since there is only one outgoing
2157 // edge in the "downhill" graph, this is guaranteed never to visit a
2158 // node more than once.
2159 let mut start = newh.len();
2160 // First, find the neighbor pass (assuming this is not the ocean).
2161 let neighbor_pass_idx = downhill[chunk_idx];
2162 let first_idx = if neighbor_pass_idx < 0 {
2163 // This is the ocean.
2164 newh.push(chunk_idx as u32);
2165 start += 1;
2166 chunk_idx
2167 } else {
2168 // This is a "real" lake.
2169 let neighbor_pass_idx = neighbor_pass_idx as usize;
2170 // Let's find our side of the pass.
2171 let pass_idx = -indirection[chunk_idx];
2172 // NOTE: Since only lakes are on the boundary, this should be a valid array
2173 // index.
2174 assert!(pass_idx >= 0);
2175 let pass_idx = pass_idx as usize;
2176 // Now, we should recompute flow paths so downhill nodes are contiguous.
2177
2178 /* // Carving strategy: reverse the path from the lake side of the pass to the
2179 // lake bottom, and also set the lake side of the pass's downhill to its
2180 // neighbor pass.
2181 let mut to_idx = neighbor_pass_idx;
2182 let mut from_idx = pass_idx;
2183 // NOTE: Since our side of the lake pass must be in the same basin as chunk_idx,
2184 // and chunk_idx is the basin bottom, we must reach it before we reach an ocean
2185 // node or other node with an invalid index.
2186 while from_idx != chunk_idx {
2187 // Reverse this (from, to) edge by first replacing to_idx with from_idx,
2188 // then replacing from_idx's downhill with the old to_idx, and finally
2189 // replacing from_idx with from_idx's old downhill.
2190 //
2191 // println!("Reversing (lake={:?}): to={:?}, from={:?}, dh={:?}", chunk_idx, to_idx, from_idx, downhill[from_idx]);
2192 from_idx = mem::replace(
2193 &mut downhill[from_idx],
2194 mem::replace(
2195 &mut to_idx,
2196 // NOTE: This cast should be valid since the node is either a path on the way
2197 // to a lake bottom, or a lake bottom with an actual pass outwards.
2198 from_idx
2199 ) as isize,
2200 ) as usize;
2201 }
2202 // Remember to set the actual lake's from_idx properly!
2203 downhill[from_idx] = to_idx as isize; */
2204
2205 // TODO: Enqueue onto newh simultaneously with filling (this could help prevent
2206 // needing a special tag just for checking if we are already enqueued).
2207 // Filling strategy: nodes are assigned paths based on cost.
2208 {
2209 assert!(tag[pass_idx] == Tag::UnParsed);
2210
2211 filling_queue.push(pass_idx);
2212 tag[neighbor_pass_idx] = Tag::WithRcv;
2213 tag[pass_idx] = Tag::InQueue;
2214
2215 let outflow_coords = uniform_idx_as_vec2(map_size_lg, neighbor_pass_idx);
2216 let elev = h(neighbor_pass_idx).max(h(pass_idx));
2217
2218 while let Some(node) = filling_queue.pop() {
2219 let coords = uniform_idx_as_vec2(map_size_lg, node);
2220
2221 let mut rcv = -1;
2222 let mut rcv_cost = -f64::INFINITY; /*f64::EPSILON;*/
2223 let outflow_distance = (outflow_coords - coords).map(|e| e as f64);
2224
2225 neighbors(map_size_lg, node).for_each(|ineighbor| {
2226 if indirection[ineighbor] != chunk_idx as i32
2227 && ineighbor != chunk_idx
2228 && ineighbor != neighbor_pass_idx
2229 || h(ineighbor) > elev
2230 {
2231 return;
2232 }
2233 let dxy = (uniform_idx_as_vec2(map_size_lg, ineighbor) - coords)
2234 .map(|e| e as f64);
2235 let neighbor_distance = /*neighbor_coef * */dxy;
2236 let tag = &mut tag[ineighbor];
2237 match *tag {
2238 Tag::WithRcv => {
2239 // TODO: Remove outdated comment.
2240 //
2241 // vec_to_outflow ⋅ (vec_to_neighbor / |vec_to_neighbor|) =
2242 // ||vec_to_outflow||cos Θ
2243 // where θ is the angle between vec_to_outflow and
2244 // vec_to_neighbor.
2245 //
2246 // Which is also the scalar component of vec_to_outflow in the
2247 // direction of vec_to_neighbor.
2248 let cost = outflow_distance
2249 .dot(neighbor_distance / neighbor_distance.magnitude());
2250 if cost > rcv_cost {
2251 rcv = ineighbor as isize;
2252 rcv_cost = cost;
2253 }
2254 },
2255 Tag::UnParsed => {
2256 filling_queue.push(ineighbor);
2257 *tag = Tag::InQueue;
2258 },
2259 _ => {},
2260 }
2261 });
2262 assert_ne!(rcv, -1);
2263 downhill[node] = rcv;
2264 tag[node] = Tag::WithRcv;
2265 }
2266 }
2267
2268 // Use our side of the pass as the initial node in the stack order.
2269 // TODO: Verify that this stack order will not "double reach" any lake chunks.
2270 neighbor_pass_idx
2271 };
2272 // New lake root
2273 let mut cur = start;
2274 let mut node = first_idx;
2275 loop {
2276 uphill(map_size_lg, downhill, node).for_each(|child| {
2277 // lake_idx is the index of our lake root.
2278 // Check to make sure child (flowing into us) is in the same lake.
2279 if indirection[child] == chunk_idx as i32 || child == chunk_idx {
2280 newh.push(child as u32);
2281 }
2282 });
2283
2284 if cur == newh.len() {
2285 break;
2286 }
2287 node = newh[cur] as usize;
2288 cur += 1;
2289 }
2290 });
2291 assert_eq!(newh.len(), downhill.len());
2292 (boundary.len(), indirection, newh.into_boxed_slice(), maxh)
2293}
2294
2295/// Iterate through set neighbors of multi-receiver flow.
2296pub fn mrec_downhill(
2297 map_size_lg: MapSizeLg,
2298 mrec: &[u8],
2299 posi: usize,
2300) -> impl Clone + Iterator<Item = (usize, usize)> {
2301 let pos = uniform_idx_as_vec2(map_size_lg, posi);
2302 let mrec_i = mrec[posi];
2303 NEIGHBOR_DELTA
2304 .iter()
2305 .enumerate()
2306 .filter(move |&(k, _)| (mrec_i >> k as isize) & 1 == 1)
2307 .map(move |(k, &(x, y))| {
2308 (
2309 k,
2310 vec2_as_uniform_idx(map_size_lg, Vec2::new(pos.x + x, pos.y + y)),
2311 )
2312 })
2313}
2314
2315/// Algorithm for computing multi-receiver flow.
2316///
2317/// * `map_size_lg`: Size of the underlying map.
2318/// * `h`: altitude
2319/// * `downhill`: single receiver
2320/// * `newh`: single receiver stack
2321/// * `wh`: buffer into which water height will be inserted.
2322/// * `nx`, `ny`: resolution in x and y directions.
2323/// * `dx`, `dy`: grid spacing in x- and y-directions
2324/// * `maxh`: maximum |height| among all nodes.
2325///
2326///
2327/// Updates the water height to a nearly planar surface, and returns a 3-tuple
2328/// containing:
2329/// * A bitmask representing which neighbors are downhill.
2330/// * Stack order for multiple receivers (from top to bottom).
2331/// * The weight for each receiver, for each node.
2332pub fn get_multi_rec<F: fmt::Debug + Float + Sync + Into<Compute>>(
2333 map_size_lg: MapSizeLg,
2334 h: impl Fn(usize) -> F + Sync,
2335 downhill: &[isize],
2336 newh: &[u32],
2337 wh: &mut [F],
2338 nx: usize,
2339 ny: usize,
2340 dx: Compute,
2341 dy: Compute,
2342 _maxh: F,
2343 threadpool: &rayon::ThreadPool,
2344) -> (Box<[u8]>, Box<[u32]>, Box<[Computex8]>) {
2345 let nn = nx * ny;
2346 let dxdy = Vec2::new(dx, dy);
2347
2348 /* // set bc
2349 let i1 = 0;
2350 let i2 = nx;
2351 let j1 = 0;
2352 let j2 = ny;
2353 let xcyclic = false;
2354 let ycyclic = false; */
2355 /*
2356 write (cbc,'(i4)') ibc
2357 i1=1
2358 i2=nx
2359 j1=1
2360 j2=ny
2361 if (cbc(4:4).eq.'1') i1=2
2362 if (cbc(2:2).eq.'1') i2=nx-1
2363 if (cbc(1:1).eq.'1') j1=2
2364 if (cbc(3:3).eq.'1') j2=ny-1
2365 xcyclic=.FALSE.
2366 ycyclic=.FALSE.
2367 if (cbc(4:4).ne.'1'.and.cbc(2:2).ne.'1') xcyclic=.TRUE.
2368 if (cbc(1:1).ne.'1'.and.cbc(3:3).ne.'1') ycyclic=.TRUE.
2369 */
2370 assert_eq!(nn, wh.len());
2371
2372 // fill the local minima with a nearly planar surface
2373 // See https://matthew-brett.github.io/teaching/floating_error.html;
2374 // our absolute error is bounded by β^(e-(p-1)), where e is the exponent of the
2375 // largest value we care about. In this case, since we are summing up to nn
2376 // numbers, we are bounded from above by nn * |maxh|; however, we only need
2377 // to invoke this when we actually encounter a number, so we compute it
2378 // dynamically. for nn + |maxh|
2379 // TODO: Consider that it's probably not possible to have a downhill path the
2380 // size of the whole grid... either measure explicitly (maybe in get_lakes)
2381 // or work out a more precise upper bound (since using nn * 2 * (maxh +
2382 // epsilon) makes f32 not work very well).
2383 let deltah = F::epsilon() + F::epsilon();
2384 newh.iter().for_each(|&ijk| {
2385 let ijk = ijk as usize;
2386 let h_i = h(ijk);
2387 let ijr = downhill[ijk];
2388 wh[ijk] = if ijr >= 0 {
2389 let ijr = ijr as usize;
2390 let wh_j = wh[ijr];
2391 if wh_j >= h_i {
2392 let deltah = deltah * wh_j.abs();
2393 wh_j + deltah
2394 } else {
2395 h_i
2396 }
2397 } else {
2398 h_i
2399 };
2400 });
2401
2402 let mut wrec = Vec::<Computex8>::with_capacity(nn);
2403 let mut mrec = Vec::with_capacity(nn);
2404 let mut don_vis = Vec::with_capacity(nn);
2405
2406 // loop on all nodes
2407 (0..nn)
2408 .into_par_iter()
2409 .map(|ij| {
2410 // TODO: SIMDify? Seems extremely amenable to that.
2411 let wh_ij = wh[ij];
2412 let mut mrec_ij = 0u8;
2413 let mut ndon_ij = 0u8;
2414 let neighbor_iter = |posi| {
2415 let pos = uniform_idx_as_vec2(map_size_lg, posi);
2416 NEIGHBOR_DELTA
2417 .iter()
2418 .map(move |&(x, y)| Vec2::new(pos.x + x, pos.y + y))
2419 .enumerate()
2420 .filter(move |&(_, pos)| {
2421 pos.x >= 0 && pos.y >= 0 && pos.x < nx as i32 && pos.y < ny as i32
2422 })
2423 .map(move |(k, pos)| (k, vec2_as_uniform_idx(map_size_lg, pos)))
2424 };
2425
2426 neighbor_iter(ij).for_each(|(k, ijk)| {
2427 let wh_ijk = wh[ijk];
2428 if wh_ij > wh_ijk {
2429 // Set neighboring edge lower than this one as being downhill.
2430 // NOTE: relying on at most 8 neighbors.
2431 mrec_ij |= 1 << k;
2432 } else if wh_ijk > wh_ij {
2433 // Avoiding ambiguous cases.
2434 ndon_ij += 1;
2435 }
2436 });
2437 (mrec_ij, (ndon_ij, ndon_ij))
2438 })
2439 .unzip_into_vecs(&mut mrec, &mut don_vis);
2440
2441 let czero = <Compute as Zero>::zero();
2442 let (wrec, stack) = threadpool.join(
2443 || {
2444 (0..nn)
2445 .into_par_iter()
2446 .map(|ij| {
2447 let mut sumweight = czero;
2448 let mut wrec = [czero; 8];
2449 let mut nrec = 0;
2450 mrec_downhill(map_size_lg, &mrec, ij).for_each(|(k, ijk)| {
2451 let lrec_ijk = ((uniform_idx_as_vec2(map_size_lg, ijk)
2452 - uniform_idx_as_vec2(map_size_lg, ij))
2453 .map(|e| e as Compute)
2454 * dxdy)
2455 .magnitude();
2456 let wrec_ijk = (wh[ij] - wh[ijk]).into() / lrec_ijk;
2457 // NOTE: To emulate single-direction flow, uncomment this line.
2458 // let wrec_ijk = if ijk as isize == downhill[ij] { <Compute as One>::one()
2459 // } else { <Compute as Zero>::zero() };
2460 wrec[k] = wrec_ijk;
2461 sumweight += wrec_ijk;
2462 nrec += 1;
2463 });
2464 let slope = sumweight / (nrec as Compute).max(1.0);
2465 let p_mfd_exp = 0.5 + 0.6 * slope;
2466 sumweight = czero;
2467 wrec.iter_mut().for_each(|wrec_k| {
2468 let wrec_ijk = wrec_k.powf(p_mfd_exp);
2469 sumweight += wrec_ijk;
2470 *wrec_k = wrec_ijk;
2471 });
2472 if sumweight > czero {
2473 wrec.iter_mut().for_each(|wrec_k| {
2474 *wrec_k /= sumweight;
2475 });
2476 }
2477 wrec
2478 })
2479 .collect_into_vec(&mut wrec);
2480 wrec
2481 },
2482 || {
2483 let mut stack = Vec::with_capacity(nn);
2484 let mut parse = Vec::with_capacity(nn);
2485
2486 // we go through the nodes
2487 (0..nn).for_each(|ij| {
2488 let (ndon_ij, _) = don_vis[ij];
2489 // when we find a "summit" (ie a node that has no donors)
2490 // we parse it (put it in a stack called parse)
2491 if ndon_ij == 0 {
2492 parse.push(ij);
2493 }
2494 // we go through the parsing stack
2495 while let Some(ijn) = parse.pop() {
2496 // we add the node to the stack
2497 stack.push(ijn as u32);
2498 mrec_downhill(map_size_lg, &mrec, ijn).for_each(|(_, ijr)| {
2499 let (_, ref mut vis_ijr) = don_vis[ijr];
2500 if *vis_ijr >= 1 {
2501 *vis_ijr -= 1;
2502 if *vis_ijr == 0 {
2503 parse.push(ijr);
2504 }
2505 }
2506 });
2507 }
2508 });
2509
2510 assert_eq!(stack.len(), nn);
2511 stack
2512 },
2513 );
2514
2515 (
2516 mrec.into_boxed_slice(),
2517 stack.into_boxed_slice(),
2518 wrec.into_boxed_slice(),
2519 )
2520}
2521
2522/// Perform erosion n times.
2523#[expect(clippy::too_many_arguments)]
2524pub fn do_erosion(
2525 map_size_lg: MapSizeLg,
2526 _max_uplift: f32,
2527 n_steps: usize,
2528 seed: &RandomField,
2529 rock_strength_nz: &(impl NoiseFn<f64, 3> + Sync),
2530 oldh: impl Fn(usize) -> f32 + Sync,
2531 oldb: impl Fn(usize) -> f32 + Sync,
2532 is_ocean: impl Fn(usize) -> bool + Sync,
2533 uplift: impl Fn(usize) -> f64 + Sync,
2534 n: impl Fn(usize) -> f32 + Sync,
2535 theta: impl Fn(usize) -> f32 + Sync,
2536 kf: impl Fn(usize) -> f64 + Sync,
2537 kd: impl Fn(usize) -> f64 + Sync,
2538 g: impl Fn(usize) -> f32 + Sync,
2539 epsilon_0: impl Fn(usize) -> f32 + Sync,
2540 alpha: impl Fn(usize) -> f32 + Sync,
2541 // scaling factors
2542 height_scale: impl Fn(f32) -> Alt + Sync,
2543 k_d_scale: f64,
2544 k_da_scale: impl Fn(f64) -> f64,
2545 threadpool: &rayon::ThreadPool,
2546 report_progress: &dyn Fn(f64),
2547) -> (Box<[Alt]>, Box<[Alt]> /* , Box<[Alt]> */) {
2548 debug!("Initializing erosion arrays...");
2549 let oldh_ = (0..map_size_lg.chunks_len())
2550 .into_par_iter()
2551 .map(|posi| oldh(posi) as Alt)
2552 .collect::<Vec<_>>()
2553 .into_boxed_slice();
2554 // Topographic basement (The height of bedrock, not including sediment).
2555 let mut b = (0..map_size_lg.chunks_len())
2556 .into_par_iter()
2557 .map(|posi| oldb(posi) as Alt)
2558 .collect::<Vec<_>>()
2559 .into_boxed_slice();
2560 // Stream power law slope exponent--link between channel slope and erosion rate.
2561 let n = (0..map_size_lg.chunks_len())
2562 .into_par_iter()
2563 .map(&n)
2564 .collect::<Vec<_>>()
2565 .into_boxed_slice();
2566 // Stream power law concavity index (θ = m/n), turned into an exponent on
2567 // drainage (which is a proxy for discharge according to Hack's Law).
2568 let m = (0..map_size_lg.chunks_len())
2569 .into_par_iter()
2570 .map(|posi| theta(posi) * n[posi])
2571 .collect::<Vec<_>>()
2572 .into_boxed_slice();
2573 // Stream power law erodability constant for fluvial erosion (bedrock)
2574 let kf = (0..map_size_lg.chunks_len())
2575 .into_par_iter()
2576 .map(&kf)
2577 .collect::<Vec<_>>()
2578 .into_boxed_slice();
2579 // Stream power law erodability constant for hillslope diffusion (bedrock)
2580 let kd = (0..map_size_lg.chunks_len())
2581 .into_par_iter()
2582 .map(&kd)
2583 .collect::<Vec<_>>()
2584 .into_boxed_slice();
2585 // Deposition coefficient
2586 let g = (0..map_size_lg.chunks_len())
2587 .into_par_iter()
2588 .map(&g)
2589 .collect::<Vec<_>>()
2590 .into_boxed_slice();
2591 let epsilon_0 = (0..map_size_lg.chunks_len())
2592 .into_par_iter()
2593 .map(&epsilon_0)
2594 .collect::<Vec<_>>()
2595 .into_boxed_slice();
2596 let alpha = (0..map_size_lg.chunks_len())
2597 .into_par_iter()
2598 .map(&alpha)
2599 .collect::<Vec<_>>()
2600 .into_boxed_slice();
2601 let mut wh = vec![0.0; map_size_lg.chunks_len()].into_boxed_slice();
2602 // TODO: Don't do this, maybe?
2603 // (To elaborate, maybe we should have varying uplift or compute it some other
2604 // way).
2605 let uplift = (0..oldh_.len())
2606 .into_par_iter()
2607 .map(|posi| uplift(posi) as f32)
2608 .collect::<Vec<_>>()
2609 .into_boxed_slice();
2610 let sum_uplift = uplift
2611 .into_par_iter()
2612 .cloned()
2613 .map(|e| e as f64)
2614 .sum::<f64>();
2615 debug!("Sum uplifts: {:?}", sum_uplift);
2616
2617 let max_uplift = uplift
2618 .into_par_iter()
2619 .cloned()
2620 .max_by(|a, b| a.partial_cmp(b).unwrap())
2621 .unwrap();
2622 let max_g = g
2623 .into_par_iter()
2624 .cloned()
2625 .max_by(|a, b| a.partial_cmp(b).unwrap())
2626 .unwrap();
2627 debug!("Max uplift: {:?}", max_uplift);
2628 debug!("Max g: {:?}", max_g);
2629 // Height of terrain, including sediment.
2630 let mut h = oldh_;
2631 // Bedrock transport coefficients (diffusivity) in m^2 / year, for sediment.
2632 // For now, we set them all to be equal
2633 // on land, but in theory we probably want to at least differentiate between
2634 // soil, bedrock, and sediment.
2635 let kdsed = 1.5e-2 / 4.0;
2636 let kdsed = kdsed * k_d_scale;
2637 let n = |posi: usize| n[posi];
2638 let m = |posi: usize| m[posi];
2639 let kd = |posi: usize| kd[posi];
2640 let kf = |posi: usize| kf[posi];
2641 let g = |posi: usize| g[posi];
2642 let epsilon_0 = |posi: usize| epsilon_0[posi];
2643 let alpha = |posi: usize| alpha[posi];
2644
2645 (0..n_steps).for_each(|i| {
2646 debug!("Erosion iteration #{:?}", i);
2647
2648 // Print out the percentage complete. Do this at most 20 times.
2649 if i % std::cmp::max(n_steps / 20, 1) == 0 {
2650 let pct = (i as f64 / n_steps as f64) * 100.0;
2651 report_progress(pct);
2652 info!("{:.2}% complete", pct);
2653 }
2654
2655 erode(
2656 map_size_lg,
2657 &mut h,
2658 &mut b,
2659 &mut wh,
2660 max_uplift,
2661 max_g,
2662 kdsed,
2663 seed,
2664 rock_strength_nz,
2665 |posi| uplift[posi],
2666 n,
2667 m,
2668 kf,
2669 kd,
2670 g,
2671 epsilon_0,
2672 alpha,
2673 &is_ocean,
2674 &height_scale,
2675 &k_da_scale,
2676 threadpool,
2677 );
2678 });
2679 (h, b)
2680}