planetsim/test_volcano.cpp

// Headless test for Live World volcanoes (placement by tectonic context + eruption / island
// building). No display needed.
//
//   g++ -std=c++17 -O2 -Isrc/sim test_volcano.cpp src/sim/IcoSphere.cpp src/sim/Planet.cpp \
//       src/sim/PlanetTectonics.cpp src/sim/PlanetDrift.cpp src/sim/PlanetErosion.cpp \
//       src/sim/PlanetHydrology.cpp src/sim/PlanetBiomes.cpp src/sim/PlanetClimate.cpp \
//       src/sim/PlanetLive.cpp src/sim/PlanetOcean.cpp src/sim/PlanetWeather.cpp \
//       src/sim/PlanetVolcano.cpp src/sim/PlanetBiota.cpp src/sim/PlanetFloraGen.cpp \
//       src/sim/PlanetFaunaGen.cpp src/sim/PlanetFungiGen.cpp src/sim/PlanetIO.cpp -o /tmp/tv && /tmp/tv
//
// Verifies: placement is deterministic + isolated from the tectonic RNG; the context classification
// (ridge / border / interior) drives where vents land and respects the probabilities; a submarine
// vent's built height is a monotonic PURE FUNCTION of liveTime that breaches sea level into an island
// and recedes when the clock steps back; and an eruption injects ash cloud at the vent.

#include "Planet.hpp"
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <vector>
#include <sstream>

static int failures = 0;
static void check(bool cond, const char* what) {
    std::printf("  [%s] %s\n", cond ? "PASS" : "FAIL", what);
    if (!cond) ++failures;
}

static void settle(Planet& p, int maxSteps = 800) {
    int run = 0;
    for (int s = 0; s < maxSteps; ++s) {
        double mc = p.step();
        if (mc < 2.0) { if (++run >= 3) break; } else run = 0;
    }
    p.computeClimate();
    p.classifyBiomes();
}
static void drift(Planet& p, int iters) {
    p.drifting = true;
    for (int k = 0; k < iters; ++k) { double dt = p.cflDtMy(); p.advect(dt); p.step(); p.erode(dt); }
}

// Replicate placeVolcanoes()'s context classification: 0 = ridge (baby plate / neighbour baby),
// 1 = plate border (a differing-plate neighbour), 2 = interior.
static int classify(const Planet& p, int i) {
    auto isBaby = [&](int pid) { return pid >= 0 && pid < (int)p.plates.size() && p.plates[pid].baby; };
    int pid = p.cells[i].plateId;
    bool ridge = isBaby(pid), border = false;
    for (int j : p.cells[i].neighbors) { int pj = p.cells[j].plateId; if (pj != pid) border = true; if (isBaby(pj)) ridge = true; }
    return ridge ? 0 : (border ? 1 : 2);
}
static bool sameVolcanoes(const std::vector<Volcano>& a, const std::vector<Volcano>& b) {
    if (a.size() != b.size()) return false;
    for (size_t i = 0; i < a.size(); ++i)
        if (a[i].id != b[i].id || a[i].cell != b[i].cell || a[i].kind != b[i].kind
            || a[i].submarine != b[i].submarine || a[i].activity != b[i].activity
            || a[i].baseElev != b[i].baseElev || a[i].tStart != b[i].tStart) return false;
    return true;
}

int main() {
    PlanetConfig cfg; cfg.subdivisions = 5; cfg.seed = 9090;
    Planet p; p.generate(cfg); settle(p); drift(p, 120);
    const int n = (int)p.cells.size();

    std::printf("Volcanoes: determinism\n");
    p.placeVolcanoes(0.0); std::vector<Volcano> first = p.volcanoes;
    p.placeVolcanoes(0.0);
    check(!first.empty(), "placeVolcanoes places a non-empty set");
    check(sameVolcanoes(first, p.volcanoes), "placeVolcanoes is deterministic (re-run identical)");

    std::printf("Volcanoes: RNG isolation from tectonics\n");
    Planet a; a.generate(cfg); settle(a);
    Planet b; b.generate(cfg); settle(b);
    for (int k = 0; k < 40; ++k) {
        double dta = a.cflDtMy(); a.advect(dta); a.step(); a.erode(dta);
        double dtb = b.cflDtMy(); b.advect(dtb); b.step(); b.erode(dtb);
        if (k == 20) b.placeVolcanoes(0.0);   // must not touch the tectonic RNG stream
    }
    bool terrainSame = true;
    for (int i = 0; i < n; ++i) if (std::fabs(a.cells[i].elevation - b.cells[i].elevation) > 1e-9) terrainSame = false;
    check(terrainSame, "placeVolcanoes never perturbs tectonic evolution");

    std::printf("Volcanoes: context classification\n");
    int eligRidge = 0, eligBorder = 0, eligInterior = 0;
    for (int i = 0; i < n; ++i) { int k = classify(p, i); if (k == 0) ++eligRidge; else if (k == 1) ++eligBorder; else ++eligInterior; }
    std::printf("    eligible cells: ridge %d, border %d, interior %d\n", eligRidge, eligBorder, eligInterior);
    // probs ridge=border=1, interior=0, no cap -> exactly the ridge+border cells, none interior.
    p.cfg.volcanoProbRidge = 1.0; p.cfg.volcanoProbBorder = 1.0; p.cfg.volcanoProbInterior = 0.0;
    p.cfg.volcanoMaxCount = 1000000;
    p.placeVolcanoes(0.0);
    bool noInterior = true; for (const Volcano& v : p.volcanoes) if (v.kind == 2) noInterior = false;
    check(noInterior, "interior prob 0 places no interior vents");
    check((int)p.volcanoes.size() == eligRidge + eligBorder, "prob 1 fills exactly the ridge+border cells");

    std::printf("Volcanoes: probability ordering (border > interior)\n");
    p.cfg.volcanoProbRidge = 1.0; p.cfg.volcanoProbBorder = 0.30; p.cfg.volcanoProbInterior = 0.05;
    p.placeVolcanoes(0.0);
    int gotRidge = 0, gotBorder = 0, gotInterior = 0;
    for (const Volcano& v : p.volcanoes) { if (v.kind == 0) ++gotRidge; else if (v.kind == 1) ++gotBorder; else ++gotInterior; }
    double rB = eligBorder ? (double)gotBorder / eligBorder : 0.0;
    double rI = eligInterior ? (double)gotInterior / eligInterior : 0.0;
    std::printf("    placement rate: border %.3f, interior %.3f\n", rB, rI);
    check(rB > rI, "border cells are far likelier to host a volcano than interior cells");
    if (eligRidge > 0) {
        double rR = (double)gotRidge / eligRidge;
        std::printf("    placement rate: ridge %.3f\n", rR);
        check(rR >= rB, "young-ridge cells are the likeliest of all");
    } else std::printf("    (no young-ridge cells this seed -- ridge rate not asserted)\n");

    std::printf("Volcanoes: build is a pure function of liveTime; submarine vent breaches into an island\n");
    Planet q; q.generate(cfg); settle(q); drift(q, 120);
    q.cfg.volcanoMaxHeight = 9000.0; q.cfg.volcanoBuildStep = 130.0; q.cfg.volcanoEruptFreq = 0.05;
    q.placeVolcanoes(0.0);
    int vi = -1; double best = -1e18;
    for (size_t k = 0; k < q.volcanoes.size(); ++k)
        if (q.volcanoes[k].submarine && q.volcanoes[k].baseElev > best) { best = q.volcanoes[k].baseElev; vi = (int)k; }
    check(vi >= 0, "at least one submarine volcano was placed");
    if (vi >= 0) {
        const Volcano v = q.volcanoes[vi];
        double b0 = q.volcanoBuilt(v, 0.0), b1 = q.volcanoBuilt(v, 5000.0),
               b2 = q.volcanoBuilt(v, 50000.0), b3 = q.volcanoBuilt(v, 500000.0);
        check(b0 <= b1 && b1 <= b2 && b2 <= b3, "built height is monotonic in liveTime");
        check(b3 > b0, "a submarine vent builds up over time");
        check(q.volcanoBuilt(v, 5000.0) == b1, "volcanoBuilt is deterministic (pure function of t)");
        q.stepVolcanoes(1.0, 500000.0);
        check(q.cells[v.cell].elevation > q.cfg.seaLevel, "submarine volcano breaches sea level into an island");
        check(!q.cells[v.cell].oceanic, "the breached island is land crust");
        // Step the clock back to the start: the island must recede (pure function of liveTime).
        q.stepVolcanoes(0.0, 0.0);
        check(q.cells[v.cell].elevation <= q.cfg.seaLevel + 1e-6, "stepping the clock back recedes the island");
        check(std::fabs(q.cells[v.cell].elevation - (v.baseElev + q.volcanoBuilt(v, 0.0))) < 1e-6,
              "vent elevation = baseElev + built(liveTime)");
    }

    std::printf("Volcanoes: an eruption injects ash cloud\n");
    Planet w; w.generate(cfg); settle(w);
    w.initWeather();
    w.computeInsolation(0.25, 0.3);
    w.placeVolcanoes(0.0);   // tStart = 0 -> at liveTime 0 every vent is at peak eruption intensity
    check(!w.volcanoes.empty(), "volcanoes placed for the ash test");
    if (!w.volcanoes.empty()) {
        std::vector<double> before = w.cloud();
        w.stepVolcanoes(1.0, 0.0);   // dtHours > 0 -> inject ash
        bool rose = false;
        for (const Volcano& vv : w.volcanoes)
            if (w.cloud()[vv.cell] > before[vv.cell] + 1e-9) rose = true;
        check(rose, "an erupting vent thickens the cloud at its cell");
    }

    std::printf("Volcanoes: save v14 round-trip\n");
    {
        std::stringstream ss(std::ios::in | std::ios::out | std::ios::binary);
        q.writeState(ss);
        Planet r;
        bool ok = r.readState(ss, true, true, true, true, true, true);
        check(ok, "readState accepts a v14 stream");
        check(sameVolcanoes(q.volcanoes, r.volcanoes), "volcano set round-trips through save");
    }

    std::printf(failures ? "\nFAILURES: %d\n" : "\nALL VOLCANO CHECKS PASSED\n", failures);
    return failures ? 1 : 0;
}