planetsim/src/sim/PlanetTectonics.cpp

#include "Planet.hpp"
#include <algorithm>
#include <cmath>
#include <cstdint>

// Deterministic, thread-independent value in [0,1). Pure function of (cell, tick,
// seed) -- so the parallel step() loops stay bit-identical for any thread count
// and this never touches Planet::rngState. `salt` selects an independent stream
// (one for the grow roll, one for the drop magnitude). splitmix64 finalizer.
static inline double peakHash(int i, uint64_t tick, uint32_t seed, uint32_t salt) {
    uint64_t x = (uint64_t)(uint32_t)i * 0x9E3779B97F4A7C15ull
               + (tick + 1) * 0xD1B54A32D192ED03ull
               + (uint64_t)seed * 0xBF58476D1CE4E5B9ull
               + (uint64_t)salt + 0x123456789ull;
    x ^= x >> 30; x *= 0xBF58476D1CE4E5B9ull;
    x ^= x >> 27; x *= 0x94D049BB133111EBull;
    x ^= x >> 31;
    return (double)(x >> 11) * (1.0 / 9007199254740992.0);  // 53-bit -> [0,1)
}

// --- Phase 1/2 tectonics: one stress->uplift->relax tick ---------------------
// Data-parallel: every pass writes only its own cell index (double-buffered
// where it reads a field it writes), so the OpenMP loops are bit-identical for
// any thread count -- determinism preserved.

double Planet::step() {
    const int n = (int)cells.size();

    // Remember elevation so we can report the largest change this tick (used by
    // the viewer to detect when the world has settled into equilibrium).
    sOldElev.resize(n);
    #pragma omp parallel for schedule(static) if(n > 20000)
    for (int i = 0; i < n; ++i) sOldElev[i] = cells[i].elevation;

    // Per-boundary-cell shaping factors (multipliers on upliftGain).
    const double collideFactor   = 1.0;  // continent interior baseline (no special boundary)
    const double oceanArcFactor  = 0.5;  // oceanic island arc (mild)
    const double trenchFactor    = 2.4;  // depth of the subduction trench
    const double riftFactor      = 0.5;  // gentle deepening of divergent rifts
    const double beltDecay       = 0.6;  // stress falloff per ring inland

    // Each pass below is data-parallel: every cell reads its own + neighbor
    // data and writes ONLY its own index (double-buffered where it reads a field
    // it also writes), so the OpenMP loops give bit-identical results for any
    // thread count -- determinism is preserved.

    // 1. Boundary stress. For each cell touching another plate, the signed
    //    convergence (positive = plates closing, negative = pulling apart),
    //    plus flags marking subduction geometry from the cell's own crust type.
    //    Scratch buffers persist across calls (sized once); cleared per step.
    sStress.assign(n, 0.0);
    sSub.assign(n, 0);                      // oceanic cell facing continental -> trench
    sOver.assign(n, 0);                     // continental cell facing oceanic -> arc
    sColl.assign(n, 0);                     // continental cell facing continental -> collision
    std::vector<double>&  stress     = sStress;
    std::vector<uint8_t>& subducting = sSub;
    std::vector<uint8_t>& overriding = sOver;
    std::vector<uint8_t>& colliding  = sColl;

    #pragma omp parallel for schedule(static) if(n > 20000)
    for (int i = 0; i < n; ++i) {
        const Cell& a = cells[i];
        Vec3 va = driftVelocity(a.plateId, a.unit);

        double convergence = 0.0;
        int boundaryNeighbors = 0;
        bool facesOceanic = false, facesContinental = false;

        for (int nb : a.neighbors) {
            const Cell& b = cells[nb];
            if (b.plateId == a.plateId) continue;
            ++boundaryNeighbors;
            Vec3 dir = (b.unit - a.unit).normalized();
            Vec3 vb  = driftVelocity(b.plateId, b.unit);
            convergence += (va - vb).dot(dir);
            if (b.oceanic) facesOceanic = true;
            else           facesContinental = true;
        }
        if (boundaryNeighbors == 0) continue;
        stress[i] = convergence / boundaryNeighbors;

        bool oceanic = a.oceanic;
        if (oceanic && facesContinental)  subducting[i] = 1;
        if (!oceanic && facesOceanic)     overriding[i] = 1;
        if (!oceanic && facesContinental) colliding[i]  = 1;  // continent-continent orogeny
    }

    // 2. Spread the convergent stress a few cell-rings into the plate interior
    //    so mountain belts have realistic width and flanks, with the peak left
    //    at the boundary (dilation with geometric decay, not blurring). Only
    //    the positive component spreads -- subduction trenches stay narrow.
    //    Boundary cells (sources) are re-anchored after each ring so they never
    //    accumulate stress from neighbouring boundaries, keeping peaks sharp.
    std::vector<double>& belt = sBelt;  belt.resize(n);
    std::vector<double>& next = sBeltNext;  next.resize(n);
    #pragma omp parallel for schedule(static) if(n > 20000)
    for (int i = 0; i < n; ++i) belt[i] = std::max(0.0, stress[i]);
    for (int it = 0; it < cfg.beltWidth; ++it) {
        #pragma omp parallel for schedule(static) if(n > 20000)
        for (int i = 0; i < n; ++i) {
            double v = belt[i];
            for (int nb : cells[i].neighbors) v = std::max(v, belt[nb] * beltDecay);
            next[i] = v;
        }
        #pragma omp parallel for schedule(static) if(n > 20000)
        for (int i = 0; i < n; ++i) if (stress[i] > 0.0) next[i] = stress[i];
        belt.swap(next);
    }

    // 3. Convert stress -> elevation change, shaped by plate geometry.
    const double gain = cfg.upliftGain;
    std::vector<double>& delta = sDelta;  delta.resize(n);
    #pragma omp parallel for schedule(static) if(n > 20000)
    for (int i = 0; i < n; ++i) {
        bool continental = !cells[i].oceanic;
        double d = 0.0;

        if (subducting[i]) {
            // Narrow, deep oceanic trench where this crust dives under a continent.
            d -= std::max(0.0, stress[i]) * gain * trenchFactor;
        } else if (belt[i] > 0.0) {
            // Orogeny boosts (collision/arc) are drift-only; forming uses the
            // original mild factors so Phase 1 settles as it always did.
            double f;
            if      (colliding[i])  f = drifting ? cfg.collisionFactor : 1.0;  // continent-continent
            else if (overriding[i]) f = drifting ? cfg.arcFactor       : 1.1;  // continental arc
            else                    f = continental ? collideFactor : oceanArcFactor;
            d += belt[i] * gain * f;
        }
        if (stress[i] < 0.0)            // divergent: rift valley / spreading
            d += stress[i] * gain * riftFactor;

        delta[i] = d;
    }

    // 4. Apply uplift, then isostatic relaxation toward the crust's base
    //    elevation. Equilibrium sits at base + delta/relax, which bounds growth
    //    well below the clamp instead of railing in a single tick.
    #pragma omp parallel for schedule(static) if(n > 20000)
    for (int i = 0; i < n; ++i) {
        double base, relaxEff;
        if (cells[i].oceanic) {
            base = oceanicBase(cells[i].geoAge);  // deepens with crustal age (always on)
            relaxEff = cfg.relax;
        } else {
            base = cfg.continentBase;
            // Drift-only: thick (high) continental crust resists isostatic
            // relaxation, so collision ranges stand and become erosion-limited
            // instead of snapping back. Forming uses full relax (original).
            if (drifting) {
                double over = std::clamp((cells[i].elevation - cfg.continentBase) / cfg.rootScale, 0.0, 1.0);
                relaxEff = cfg.relax * (1.0 - cfg.isostaticPersist * over);
            } else {
                relaxEff = cfg.relax;
            }
        }
        double d = delta[i];
        // Soft probabilistic peak cap (drift-only): above peakSoftCapStart the
        // chance that uplift takes falls linearly to 0 at peakSoftCapEnd, so peaks
        // spread across a height band instead of railing at the old hard ceiling.
        // A lost roll forfeits this tick's uplift and shaves a random 0..peakFailDrop
        // metres off (so a peak hovers near its own height). erodeIter is the tick
        // counter (step/erode run 1:1 in drift; it's saved -> exact resume).
        if (drifting && d > 0.0 && cells[i].elevation > cfg.peakSoftCapStart) {
            double span = cfg.peakSoftCapEnd - cfg.peakSoftCapStart;
            double p = span > 0.0 ? (cfg.peakSoftCapEnd - cells[i].elevation) / span : 0.0;
            p = std::clamp(p, 0.0, 1.0);
            if (peakHash(i, (uint64_t)erodeIter, cfg.seed, 0u) >= p) {   // grow roll lost
                d = 0.0;
                cells[i].elevation -= peakHash(i, (uint64_t)erodeIter, cfg.seed, 1u) * cfg.peakFailDrop;
            }
        }
        cells[i].elevation += d;
        cells[i].elevation += (base - cells[i].elevation) * relaxEff;
        // (geoAge is crust age in My, advanced by advect() in Phase 2.)
    }

    // 5. Gentle diffusion (erosion / sediment transport) so relief stays smooth.
    std::vector<double>& smoothed = sSmoothed;  smoothed.resize(n);
    #pragma omp parallel for schedule(static) if(n > 20000)
    for (int i = 0; i < n; ++i) {
        double sum = cells[i].elevation; int cnt = 1;
        for (int nb : cells[i].neighbors) { sum += cells[nb].elevation; ++cnt; }
        smoothed[i] = cells[i].elevation * 0.96 + (sum / cnt) * 0.04;
    }
    // Upper clamp tracks peakSoftCapEnd so the soft cap governs peak heights
    // (a fixed 9000 m ceiling would re-flatten everything the soft cap allows);
    // it's just a safety rail now. Lower clamp (-11000 m, trenches) is unchanged.
    const double elevCeil = cfg.peakSoftCapEnd;
    #pragma omp parallel for schedule(static) if(n > 20000)
    for (int i = 0; i < n; ++i)
        cells[i].elevation = std::clamp(smoothed[i], -11000.0, elevCeil);

    // Largest elevation change this tick -> 0 as the world reaches equilibrium.
    double maxChange = 0.0;
    #pragma omp parallel for schedule(static) reduction(max:maxChange) if(n > 20000)
    for (int i = 0; i < n; ++i)
        maxChange = std::max(maxChange, std::fabs(cells[i].elevation - sOldElev[i]));
    return maxChange;
}