TerrainGen coverage

terraingen/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
176 / 177
99.4%
Branches:
124 / 135
91.9%
Functions:
23 / 23
100.0%
layers/craters/degradation/diffusion.cpp
Line Branch Exec Source
1 #include "terraingen/layers/craters/degradation/diffusion.h"
2
3 #include "terraingen/util/checks.h"
4
5 #include <algorithm>
6 #include <atomic>
7 #include <bit>
8 #include <cmath>
9 #include <cstdint>
10 #include <stdexcept>
11 #include <utility>
12
13 namespace terraingen {
14
15 namespace {
16
17 /// A number no other diffused profile has, so a thread's last table is never another's.
18 54 std::uint64_t NextInstance() {
19 static std::atomic<std::uint64_t> next{1};
20 54 return next.fetch_add(1);
21 }
22
23 constexpr double kFirstAge = 1e-5;
24 constexpr double kLastAge = 0.5;
25 constexpr std::size_t kAgeSteps = 160;
26 constexpr double kReferenceDiameter = 1000.0;
27
28 /// The tabulated diffusion ages: 0, then geometric from `kFirstAge` to `kLastAge`.
29 6247810 const std::vector<double>& TableAges() {
30 29 static const std::vector<double> ages = [] {
31
1/1
✓ Branch 1 taken 29 times.
 58 std::vector<double> values{0.0};
32 29 const double ratio =
33 std::pow(kLastAge / kFirstAge, 1.0 / static_cast<double>(kAgeSteps - 1));
34
2/2
✓ Branch 0 taken 4640 times.
✓ Branch 1 taken 29 times.
 4669 for (std::size_t k = 0; k < kAgeSteps; ++k) {
35 4640 values.push_back(k + 1 == kAgeSteps
36
3/3
✓ Branch 0 taken 4611 times.
✓ Branch 1 taken 29 times.
✓ Branch 3 taken 4640 times.
 9251 ? kLastAge
37 4611 : kFirstAge * std::pow(ratio, static_cast<double>(k)));
38 }
39 29 return values;
40
4/7
✓ Branch 0 taken 29 times.
✓ Branch 1 taken 6247781 times.
✓ Branch 3 taken 29 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 29 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 6247810 }();
41 6247810 return ages;
42 }
43
44 /// The row below `tau` and the weight of the row above it, interpolating in log tau (linearly in
45 /// tau before the first nonzero age).
46 6037977 std::pair<std::size_t, double> RowOf(double tau) {
47
1/1
✓ Branch 1 taken 6037977 times.
 6037977 const std::vector<double>& ages = TableAges();
48
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6037977 times.
 6037977 if (tau <= 0.0) {
49 ✗ return {0, 0.0};
50 }
51
2/2
✓ Branch 1 taken 156606 times.
✓ Branch 2 taken 5881371 times.
 6037977 if (tau >= ages.back()) {
52 156606 return {ages.size() - 1, 0.0};
53 }
54
2/2
✓ Branch 1 taken 45049 times.
✓ Branch 2 taken 5836322 times.
 5881371 if (tau < ages[1]) {
55 45049 return {0, tau / ages[1]};
56 }
57 const auto above =
58
1/1
✓ Branch 4 taken 5836322 times.
 5836322 static_cast<std::size_t>(std::upper_bound(ages.begin(), ages.end(), tau) - ages.begin());
59 5836322 const std::size_t below = above - 1;
60 5836322 return {below, std::log(tau / ages[below]) / std::log(ages[above] / ages[below])};
61 }
62
63 /// A tridiagonal system lower[i] x[i-1] + diag[i] x[i] + upper[i] x[i+1] = rhs[i] for rows i < n,
64 /// with x[n] given.
65 struct Tridiagonal {
66 std::vector<double> lower;
67 std::vector<double> diag;
68 std::vector<double> upper;
69 };
70
71 /// Solves `system` by the Thomas algorithm: `x` holds the right-hand side of rows 0 to n - 1 and
72 /// the fixed x[n], and is replaced by the solution. `c` and `d` are scratch of n values.
73 280973 void SolveTridiagonal(const Tridiagonal& system, std::vector<double>& x, std::vector<double>& c,
74 std::vector<double>& d) {
75 280973 const std::size_t n = system.diag.size();
76 280973 c[0] = system.upper[0] / system.diag[0];
77 280973 d[0] = x[0] / system.diag[0];
78
2/2
✓ Branch 0 taken 168302827 times.
✓ Branch 1 taken 280973 times.
 168583800 for (std::size_t i = 1; i < n; ++i) {
79 168302827 const double m = system.diag[i] - system.lower[i] * c[i - 1];
80 168302827 c[i] = system.upper[i] / m;
81 168302827 d[i] = (x[i] - system.lower[i] * d[i - 1]) / m;
82 }
83
2/2
✓ Branch 0 taken 168583800 times.
✓ Branch 1 taken 280973 times.
 168864773 for (std::size_t i = n; i-- > 0;) {
84 168583800 x[i] = d[i] - c[i] * x[i + 1];
85 }
86 280973 }
87
88 } // namespace
89
90 5480250 double PowerLawDiffusivity::Diffusivity(double diameter) const {
91 5480250 return coefficient * std::pow(std::fmax(diameter, minDiameter) / 1000.0, exponent);
92 }
93
94 5480245 double PowerLawDiffusivity::operator()(double diameter, double ageGyr) const {
95 5480245 return Diffusivity(diameter) * ageGyr * 1000.0;
96 }
97
98 8 void ValidatePowerLawDiffusivity(const PowerLawDiffusivity& law) {
99
1/1
✓ Branch 2 taken 5 times.
 8 RequirePositive(law.coefficient, "A diffusivity coefficient");
100
6/6
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 time.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 3 times.
 5 if (!(std::isfinite(law.exponent) && IsNonNegative(law.minDiameter))) {
101 throw std::invalid_argument(
102
1/1
✓ Branch 2 taken 2 times.
 2 "A diffusivity exponent must be finite and its minimum diameter not negative");
103 }
104 3 }
105
106 2171942 std::size_t DiffusedProfile::KeyHash::operator()(const Key& key) const {
107 2171942 std::uint64_t hash = 1469598103934665603ull;
108
2/2
✓ Branch 2 taken 6515826 times.
✓ Branch 3 taken 2171942 times.
 8687768 for (const double part : {key.depthRatio, key.rimHeightRatio, key.diameter}) {
109 6515826 hash = (hash ^ std::bit_cast<std::uint64_t>(part)) * 1099511628211ull;
110 }
111 2171942 return static_cast<std::size_t>(hash);
112 }
113
114 54 DiffusedProfile::DiffusedProfile(std::shared_ptr<const CraterProfile> fresh, double reachRadii,
115 54 std::size_t cells, std::size_t cacheTables)
116 54 : fresh_(std::move(fresh)), reachRadii_(reachRadii), cells_(cells), instance_(NextInstance()) {
117
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 53 times.
 54 if (!fresh_) {
118
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("A diffused profile needs a fresh profile");
119 }
120
6/7
✓ Branch 1 taken 53 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 53 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 51 times.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 51 times.
 53 if (!(std::isfinite(reachRadii_) && reachRadii_ >= fresh_->ReachRadii())) {
121 throw std::invalid_argument(
122
1/1
✓ Branch 2 taken 2 times.
 2 "A diffused profile must reach at least as far as the fresh one");
123 }
124
3/4
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✓ Branch 3 taken 50 times.
 51 if (cells_ < 100 || cacheTables == 0) {
125
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("A diffused profile needs at least 100 cells and a cache");
126 }
127
1/1
✓ Branch 1 taken 50 times.
 100 tables_ = std::make_shared<LruCache<Key, Table, KeyHash>>(
128 139 cacheTables, [](const Table&) { return std::size_t{1}; });
129 62 }
130
131 209744 const std::vector<double>& DiffusedProfile::Ages() const {
132 209744 return TableAges();
133 }
134
135 89 DiffusedProfile::Table DiffusedProfile::Solve(const CraterMorphometry& fresh,
136 double diameter) const {
137 89 const std::size_t n = cells_;
138 89 const double extent = 0.5 * reachRadii_;
139 89 const double dx = extent / static_cast<double>(n);
140 89 const CraterState initial{fresh, fresh, 0.0};
141
1/1
✓ Branch 1 taken 89 times.
 89 std::vector<double> h(n + 1, 0.0);
142
2/2
✓ Branch 0 taken 53400 times.
✓ Branch 1 taken 89 times.
 53489 for (std::size_t i = 0; i < n; ++i) {
143
1/1
✓ Branch 2 taken 53400 times.
 53400 h[i] = fresh_->Height(static_cast<double>(i) * dx * diameter, diameter, initial) / diameter;
144 }
145
1/1
✓ Branch 1 taken 89 times.
 89 std::vector<double> below(n + 1, 0.0);
146
1/1
✓ Branch 1 taken 89 times.
 89 std::vector<double> above(n + 1, 0.0);
147
2/2
✓ Branch 0 taken 53311 times.
✓ Branch 1 taken 89 times.
 53400 for (std::size_t i = 1; i < n; ++i) {
148 53311 const double x = static_cast<double>(i) * dx;
149 53311 below[i] = (x - 0.5 * dx) / (x * dx * dx);
150 53311 above[i] = (x + 0.5 * dx) / (x * dx * dx);
151 }
152 89 above[0] = 4.0 / (dx * dx);
153
154 89 Table table;
155
1/1
✓ Branch 1 taken 89 times.
 89 const std::vector<double>& ages = TableAges();
156
1/1
✓ Branch 2 taken 89 times.
 89 table.heights.reserve(ages.size() * (n + 1));
157
3/7
✓ Branch 1 taken 89 times.
✓ Branch 4 taken 89 times.
✓ Branch 7 taken 89 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 267 Tridiagonal system{std::vector<double>(n), std::vector<double>(n), std::vector<double>(n)};
158
1/1
✓ Branch 1 taken 89 times.
 178 std::vector<double> c(n);
159
1/1
✓ Branch 1 taken 89 times.
 89 std::vector<double> d(n);
160 280973 const auto step = [&](double dt) {
161
2/2
✓ Branch 0 taken 168583800 times.
✓ Branch 1 taken 280973 times.
 168864773 for (std::size_t i = 0; i < n; ++i) {
162 168583800 system.lower[i] = -dt * below[i];
163 168583800 system.diag[i] = 1.0 + dt * (below[i] + above[i]);
164 168583800 system.upper[i] = -dt * above[i];
165 }
166 280973 h[n] = 0.0;
167 280973 SolveTridiagonal(system, h, c, d);
168 281062 };
169 89 double tau = 0.0;
170
2/2
✓ Branch 4 taken 14329 times.
✓ Branch 5 taken 89 times.
 14418 for (const double target : ages) {
171
2/2
✓ Branch 0 taken 280973 times.
✓ Branch 1 taken 14329 times.
 295302 while (tau < target) {
172 280973 const double dt = std::fmin(target - tau, std::fmax(dx * dx, 0.002 * tau));
173 280973 step(dt);
174 280973 tau += dt;
175 }
176 14329 double rim = -INFINITY;
177 14329 double floor = INFINITY;
178
2/2
✓ Branch 0 taken 8611729 times.
✓ Branch 1 taken 14329 times.
 8626058 for (std::size_t i = 0; i <= n; ++i) {
179 8611729 const double x = static_cast<double>(i) * dx;
180
2/2
✓ Branch 0 taken 7178829 times.
✓ Branch 1 taken 1432900 times.
 8611729 if (x >= 0.25) {
181 7178829 rim = std::fmax(rim, h[i]);
182 }
183
2/2
✓ Branch 0 taken 2880129 times.
✓ Branch 1 taken 5731600 times.
 8611729 if (x <= 0.5) {
184 2880129 floor = std::fmin(floor, h[i]);
185 }
186
1/1
✓ Branch 2 taken 8611729 times.
 8611729 table.heights.push_back(static_cast<float>(h[i]));
187 }
188
1/1
✓ Branch 1 taken 14329 times.
 14329 table.morphometry.push_back({rim - floor, rim});
189 }
190 178 return table;
191 89 }
192
193 std::shared_ptr<const DiffusedProfile::Table>
194 7690144 DiffusedProfile::TableFor(const CraterMorphometry& fresh, double diameter) const {
195
3/3
✓ Branch 2 taken 7690144 times.
✓ Branch 4 taken 24900 times.
✓ Branch 5 taken 7665244 times.
 7690144 const double reference = fresh_->IsComplex(diameter) ? diameter : kReferenceDiameter;
196 7690144 const Key key{fresh.depthRatio, fresh.rimHeightRatio, reference};
197 struct Last {
198 std::uint64_t instance = 0;
199 Key key{};
200 std::shared_ptr<const Table> table;
201 };
202
2/2
✓ Branch 0 taken 31 times.
✓ Branch 1 taken 7690113 times.
 7690144 thread_local Last last;
203
6/6
✓ Branch 0 taken 7690103 times.
✓ Branch 1 taken 41 times.
✓ Branch 3 taken 6604182 times.
✓ Branch 4 taken 1085921 times.
✓ Branch 5 taken 6604182 times.
✓ Branch 6 taken 1085962 times.
 7690144 if (last.instance == instance_ && last.key == key) {
204 6604182 return last.table;
205 }
206
2/2
✓ Branch 2 taken 1085962 times.
✓ Branch 5 taken 1085962 times.
 1086051 last.table = tables_->Get(key, [&] { return Solve(fresh, reference); });
207 1085962 last.instance = instance_;
208 1085962 last.key = key;
209 1085962 return last.table;
210 }
211
212 std::shared_ptr<const std::vector<float>>
213 1652167 DiffusedProfile::HeightTable(const CraterMorphometry& fresh, double diameter) const {
214
1/1
✓ Branch 1 taken 1652167 times.
 1652167 const std::shared_ptr<const Table> table = TableFor(fresh, diameter);
215 3304334 return {table, &table->heights};
216 1652167 }
217
218 836212 double DiffusedProfile::Height(double distance, double diameter, const CraterState& state) const {
219
2/2
✓ Branch 0 taken 273412 times.
✓ Branch 1 taken 562800 times.
 836212 if (!(state.diffusion > 0.0)) {
220
1/1
✓ Branch 2 taken 273412 times.
 273412 return fresh_->Height(distance, diameter, {state.fresh, state.fresh, 0.0});
221 }
222 562800 const double x = distance / diameter;
223 562800 const double extent = 0.5 * reachRadii_;
224
2/2
✓ Branch 0 taken 26 times.
✓ Branch 1 taken 562774 times.
 562800 if (x >= extent) {
225 26 return 0.0;
226 }
227
1/1
✓ Branch 1 taken 562774 times.
 562774 const auto table = TableFor(state.fresh, diameter);
228 562774 const double cell = x / extent * static_cast<double>(cells_);
229 562774 const auto i = std::min(static_cast<std::size_t>(cell), cells_ - 1);
230 562774 const double t = cell - static_cast<double>(i);
231
1/1
✓ Branch 1 taken 562774 times.
 562774 const auto [row, weight] = RowOf(state.diffusion);
232 562774 const std::size_t width = cells_ + 1;
233 973106 const auto at = [&](std::size_t r) {
234 973106 const float* values = table->heights.data() + r * width;
235 973106 return (1.0 - t) * static_cast<double>(values[i]) + t * static_cast<double>(values[i + 1]);
236 562774 };
237 562774 double value = at(row);
238
2/2
✓ Branch 0 taken 410332 times.
✓ Branch 1 taken 152442 times.
 562774 if (weight > 0.0) {
239 410332 value += weight * (at(row + 1) - value);
240 }
241 562774 return diameter * value;
242 562774 }
243
244 5480252 CraterMorphometry DiffusedProfile::Morphometry(const CraterMorphometry& fresh, double diameter,
245 double tau) const {
246
2/2
✓ Branch 0 taken 5049 times.
✓ Branch 1 taken 5475203 times.
 5480252 if (!(tau > 0.0)) {
247 5049 return fresh;
248 }
249
1/1
✓ Branch 1 taken 5475203 times.
 5475203 const auto table = TableFor(fresh, diameter);
250
1/1
✓ Branch 1 taken 5475203 times.
 5475203 const auto [row, weight] = RowOf(tau);
251 5475203 CraterMorphometry result = table->morphometry[row];
252
2/2
✓ Branch 0 taken 5471038 times.
✓ Branch 1 taken 4165 times.
 5475203 if (weight > 0.0) {
253 5471038 const CraterMorphometry& next = table->morphometry[row + 1];
254 5471038 result.depthRatio += weight * (next.depthRatio - result.depthRatio);
255 5471038 result.rimHeightRatio += weight * (next.rimHeightRatio - result.rimHeightRatio);
256 }
257 5475203 return result;
258 5475203 }
259
260 47 DiffusionDegradation::DiffusionDegradation(FreshMorphometry fresh, DiffusionAgeModel ages,
261 47 std::shared_ptr<const DiffusedProfile> profile)
262 47 : fresh_(std::move(fresh)), ages_(std::move(ages)), profile_(std::move(profile)) {
263
7/8
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 1 time.
✓ Branch 4 taken 46 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 time.
✓ Branch 8 taken 45 times.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 45 times.
 47 if (!fresh_ || !ages_ || !profile_) {
264 throw std::invalid_argument(
265
1/1
✓ Branch 2 taken 2 times.
 2 "Diffusion degradation needs fresh morphometry, diffusion ages and a diffused profile");
266 }
267 53 }
268
269 5479640 CraterState DiffusionDegradation::State(const Crater& crater) const {
270 5479640 return WornState(crater, 1.0);
271 }
272
273 5480238 CraterState DiffusionDegradation::WornState(const Crater& crater, double rate) const {
274
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 5480237 times.
 5480238 if (!IsPositive(crater.diameter)) {
275
1/1
✓ Branch 2 taken 1 time.
 1 throw std::domain_error("Crater diameter must be positive and finite");
276 }
277
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 5480236 times.
 5480237 if (!IsNonNegative(crater.ageGyr)) {
278
1/1
✓ Branch 2 taken 1 time.
 1 throw std::domain_error("Crater age must be finite and not negative");
279 }
280
1/1
✓ Branch 1 taken 5480236 times.
 5480236 const double kappaT = rate * ages_(crater.diameter, crater.ageGyr);
281
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 5480235 times.
 5480236 if (!IsNonNegative(kappaT)) {
282
1/1
✓ Branch 2 taken 1 time.
 1 throw std::domain_error("A crater's diffusion age must be finite and not negative");
283 }
284
1/1
✓ Branch 1 taken 5480235 times.
 5480235 const CraterMorphometry fresh = fresh_(crater.diameter);
285 5480235 const double tau = kappaT / (crater.diameter * crater.diameter);
286
1/1
✓ Branch 2 taken 5480235 times.
 10960470 return {fresh, profile_->Morphometry(fresh, crater.diameter, tau), tau};
287 }
288
289 } // namespace terraingen
290