feat(core): add mock flow solver with conservation-based distribution

Treats each relation network as a conductance network (diameter^4 / length),
solves nodal balance per connected component, so flow is conserved at every
junction and wider pipes carry proportionally more flow.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

src/CMakeLists.txt

@@ -8,6 +8,8 @@ add_library(diagcore STATIC
     core/NetworkModel.cpp
     core/Router.h
     core/Router.cpp
+    core/FlowSolver.h
+    core/FlowSolver.cpp
 )
 target_include_directories(diagcore PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
 target_link_libraries(diagcore PUBLIC Qt6::Core Qt6::Gui)

src/core/FlowSolver.cpp

@@ -0,0 +1,177 @@
+#include "FlowSolver.h"
+
+#include "NetworkModel.h"
+
+#include <QtMath>
+#include <QVector>
+
+#include <functional>
+#include <numeric>
+
+namespace diag {
+
+namespace {
+
+// Conductance mock: wider and shorter pipes carry more flow.
+double conductance(const NetworkModel& model, QUuid edgeId)
+{
+    const double d = model.edgeProperty(edgeId, QStringLiteral("diameter")).toDouble();
+    const double len = model.edgeProperty(edgeId, QStringLiteral("length")).toDouble();
+    const double dm = d > 0 ? d : 100.0;
+    const double lm = len > 0 ? len : 10.0;
+    return qPow(dm / 100.0, 4) / lm;
+}
+
+double injection(const NetworkModel& model, QUuid nodeId)
+{
+    const double supply = model.nodeProperty(nodeId, QStringLiteral("supply")).toDouble();
+    const double demand = model.nodeProperty(nodeId, QStringLiteral("demand")).toDouble();
+    return supply - demand;
+}
+
+// Solve A x = b in-place; returns false on a singular matrix.
+bool gauss(QVector<QVector<double>>& a, QVector<double>& b, QVector<double>& x)
+{
+    const int n = b.size();
+    for (int col = 0; col < n; ++col) {
+        int pivot = col;
+        for (int row = col + 1; row < n; ++row)
+            if (qAbs(a[row][col]) > qAbs(a[pivot][col]))
+                pivot = row;
+        if (qAbs(a[pivot][col]) < 1e-12)
+            return false;
+        a.swapItemsAt(col, pivot);
+        std::swap(b[col], b[pivot]);
+        for (int row = col + 1; row < n; ++row) {
+            const double f = a[row][col] / a[col][col];
+            if (f == 0.0)
+                continue;
+            for (int k = col; k < n; ++k)
+                a[row][k] -= f * a[col][k];
+            b[row] -= f * b[col];
+        }
+    }
+    x.resize(n);
+    for (int row = n - 1; row >= 0; --row) {
+        double sum = b[row];
+        for (int k = row + 1; k < n; ++k)
+            sum -= a[row][k] * x[k];
+        x[row] = sum / a[row][row];
+    }
+    return true;
+}
+
+}  // namespace
+
+FlowResult FlowSolver::solve(const NetworkModel& model, const QString& relationId)
+{
+    FlowResult result;
+
+    // Edges of this relation and the nodes they touch.
+    QList<QUuid> edges;
+    QList<QUuid> nodes;
+    QHash<QUuid, int> nodeIndex;
+    for (QUuid id : model.edgeIds()) {
+        const Edge* e = model.edge(id);
+        if (!e || e->relation != relationId)
+            continue;
+        edges.append(id);
+        for (QUuid n : {e->fromNode, e->toNode}) {
+            if (!nodeIndex.contains(n)) {
+                nodeIndex.insert(n, nodes.size());
+                nodes.append(n);
+            }
+        }
+    }
+    if (edges.isEmpty()) {
+        result.ok = true;
+        result.message = QStringLiteral("No %1 edges").arg(relationId);
+        return result;
+    }
+
+    // Connected components (union-find).
+    QVector<int> parent(nodes.size());
+    for (int i = 0; i < parent.size(); ++i)
+        parent[i] = i;
+    std::function<int(int)> find = [&](int i) {
+        while (parent[i] != i)
+            i = parent[i] = parent[parent[i]];
+        return i;
+    };
+    for (QUuid id : edges) {
+        const Edge* e = model.edge(id);
+        parent[find(nodeIndex[e->fromNode])] = find(nodeIndex[e->toNode]);
+    }
+
+    // Per component: balance injections to zero-sum, then solve G*p = inj
+    // with the component's last node grounded (p = 0).
+    QHash<int, QList<int>> components;  // root -> node indices
+    for (int i = 0; i < nodes.size(); ++i)
+        components[find(i)].append(i);
+
+    QHash<QUuid, double> potentials;
+    for (auto it = components.begin(); it != components.end(); ++it) {
+        const QList<int>& comp = it.value();
+        const int n = comp.size();
+        QHash<int, int> local;  // global node index -> local index
+        for (int i = 0; i < n; ++i)
+            local.insert(comp.at(i), i);
+
+        QVector<double> inj(n, 0.0);
+        for (int i = 0; i < n; ++i)
+            inj[i] = injection(model, nodes.at(comp.at(i)));
+        const double imbalance = std::accumulate(inj.begin(), inj.end(), 0.0);
+        for (int i = 0; i < n; ++i)
+            inj[i] -= imbalance / n;
+
+        if (n < 2)
+            continue;
+
+        // Laplacian, reduced by grounding the last local node.
+        const int m = n - 1;
+        QVector<QVector<double>> a(m, QVector<double>(m, 0.0));
+        QVector<double> b(m, 0.0);
+        for (int i = 0; i < m; ++i)
+            b[i] = inj[i];
+        for (QUuid id : edges) {
+            const Edge* e = model.edge(id);
+            if (!local.contains(nodeIndex[e->fromNode]))
+                continue;  // edge belongs to another component
+            const int u = local[nodeIndex[e->fromNode]];
+            const int v = local[nodeIndex[e->toNode]];
+            const double g = conductance(model, id);
+            if (u < m)
+                a[u][u] += g;
+            if (v < m)
+                a[v][v] += g;
+            if (u < m && v < m) {
+                a[u][v] -= g;
+                a[v][u] -= g;
+            }
+        }
+        QVector<double> p;
+        if (!gauss(a, b, p)) {
+            result.message = QStringLiteral("Singular network (parallel zero paths)");
+            return result;
+        }
+        for (int i = 0; i < n; ++i)
+            potentials.insert(nodes.at(comp.at(i)), i < m ? p[i] : 0.0);
+    }
+
+    for (QUuid id : edges) {
+        const Edge* e = model.edge(id);
+        const double g = conductance(model, id);
+        const double flow = g * (potentials.value(e->fromNode) - potentials.value(e->toNode));
+        result.edgeFlow.insert(id, flow);
+    }
+    result.ok = true;
+    return result;
+}
+
+void FlowSolver::apply(NetworkModel& model, const FlowResult& result)
+{
+    for (auto it = result.edgeFlow.constBegin(); it != result.edgeFlow.constEnd(); ++it)
+        model.setEdgeFlow(it.key(), it.value());
+}
+
+}  // namespace diag

src/core/FlowSolver.h

@@ -0,0 +1,29 @@
+#pragma once
+
+#include "Types.h"
+
+#include <QHash>
+
+namespace diag {
+
+class NetworkModel;
+
+struct FlowResult {
+    bool ok = false;
+    QString message;
+    QHash<QUuid, double> edgeFlow;  // signed; positive flows from -> to
+};
+
+// Mock flow-distribution calculation. Treats each relation network as a
+// resistor-style network: edge conductance derives from pipe diameter and
+// length, node injections from "supply"/"demand" properties. Solves nodal
+// balance (Kirchhoff) per connected component with dense Gaussian
+// elimination, so flow is conserved at every junction.
+class FlowSolver {
+public:
+    static FlowResult solve(const NetworkModel& model, const QString& relationId);
+    // Writes flows into the model (setEdgeFlow); edges of other relations keep theirs.
+    static void apply(NetworkModel& model, const FlowResult& result);
+};
+
+}  // namespace diag

tests/CMakeLists.txt

@@ -9,3 +9,4 @@ endfunction()
 diag_add_test(tst_model diagcore)
 diag_add_test(tst_properties diagcore)
 diag_add_test(tst_router diagcore)
+diag_add_test(tst_solver diagcore)

tests/tst_solver.cpp

@@ -0,0 +1,123 @@
+#include "TestFixtures.h"
+#include "core/FlowSolver.h"
+#include "core/NetworkModel.h"
+
+#include <QSignalSpy>
+#include <QTest>
+
+using namespace diag;
+
+class TestSolver : public QObject {
+    Q_OBJECT
+private slots:
+    void init()
+    {
+        fixtures::fillRegistry(m_reg);
+        delete m_model;
+        m_model = new NetworkModel(&m_reg, this);
+    }
+
+    void chain_carriesSupplyThrough()
+    {
+        // source(10) -> pump -> consumer(10)
+        const QUuid src = m_model->addNode(QStringLiteral("source"), {0, 0});
+        const QUuid pump = m_model->addNode(QStringLiteral("pump"), {200, 0});
+        const QUuid c = m_model->addNode(QStringLiteral("consumer"), {400, 0});
+        const QUuid e1 = m_model->addEdge(src, "out", pump, "in");
+        const QUuid e2 = m_model->addEdge(pump, "out", c, "in");
+
+        const FlowResult r = FlowSolver::solve(*m_model, QStringLiteral("water"));
+        QVERIFY(r.ok);
+        QVERIFY(qAbs(r.edgeFlow.value(e1) - 10.0) < 1e-6);
+        QVERIFY(qAbs(r.edgeFlow.value(e2) - 10.0) < 1e-6);
+    }
+
+    void junction_conservesFlow()
+    {
+        // source(10) -> junction -> two consumers (arbitrary demands).
+        const QUuid src = m_model->addNode(QStringLiteral("source"), {0, 0});
+        const QUuid j = m_model->addNode(QStringLiteral("junction"), {200, 0});
+        const QUuid c1 = m_model->addNode(QStringLiteral("consumer"), {400, -100});
+        const QUuid c2 = m_model->addNode(QStringLiteral("consumer"), {400, 100});
+        m_model->setNodeProperty(c1, "demand", 4.0);
+        m_model->setNodeProperty(c2, "demand", 6.0);
+
+        const QUuid eIn = m_model->addEdge(src, "out", j, "w");
+        const QUuid eOut1 = m_model->addEdge(j, "n", c1, "in");
+        const QUuid eOut2 = m_model->addEdge(j, "s", c2, "in");
+
+        const FlowResult r = FlowSolver::solve(*m_model, QStringLiteral("water"));
+        QVERIFY(r.ok);
+        const double in = r.edgeFlow.value(eIn);
+        const double out = r.edgeFlow.value(eOut1) + r.edgeFlow.value(eOut2);
+        QVERIFY(qAbs(in - out) < 1e-6);          // Kirchhoff at the junction
+        QVERIFY(qAbs(in - 10.0) < 1e-6);         // all supply enters
+        QVERIFY(qAbs(r.edgeFlow.value(eOut1) - 4.0) < 1e-6);
+        QVERIFY(qAbs(r.edgeFlow.value(eOut2) - 6.0) < 1e-6);
+    }
+
+    void parallelPipes_widerCarriesMore()
+    {
+        // Two parallel pipes between two junctions; DN150 vs DN75.
+        const QUuid src = m_model->addNode(QStringLiteral("source"), {0, 0});
+        const QUuid j1 = m_model->addNode(QStringLiteral("junction"), {200, 0});
+        const QUuid j2 = m_model->addNode(QStringLiteral("junction"), {400, 0});
+        const QUuid c = m_model->addNode(QStringLiteral("consumer"), {600, 0});
+        m_model->addEdge(src, "out", j1, "w");
+        const QUuid wide = m_model->addEdge(j1, "n", j2, "n");
+        const QUuid narrow = m_model->addEdge(j1, "s", j2, "s");
+        m_model->addEdge(j2, "e", c, "in");
+        m_model->setEdgeProperty(wide, "diameter", 150.0);
+        m_model->setEdgeProperty(narrow, "diameter", 75.0);
+
+        const FlowResult r = FlowSolver::solve(*m_model, QStringLiteral("water"));
+        QVERIFY(r.ok);
+        const double fWide = r.edgeFlow.value(wide);
+        const double fNarrow = r.edgeFlow.value(narrow);
+        QVERIFY(fWide > 0 && fNarrow > 0);
+        QVERIFY(qAbs(fWide + fNarrow - 10.0) < 1e-6);
+        // Conductance ~ d^4: ratio (150/75)^4 = 16.
+        QVERIFY(qAbs(fWide / fNarrow - 16.0) < 1e-6);
+    }
+
+    void unbalancedInjections_areNormalized()
+    {
+        // Supply 10 vs demand 4: imbalance is spread, but flow stays conserved.
+        const QUuid src = m_model->addNode(QStringLiteral("source"), {0, 0});
+        const QUuid c = m_model->addNode(QStringLiteral("consumer"), {200, 0});
+        m_model->setNodeProperty(c, "demand", 4.0);
+        const QUuid e = m_model->addEdge(src, "out", c, "in");
+
+        const FlowResult r = FlowSolver::solve(*m_model, QStringLiteral("water"));
+        QVERIFY(r.ok);
+        QVERIFY(qAbs(r.edgeFlow.value(e) - 7.0) < 1e-6);  // (10 + 4) / 2
+    }
+
+    void otherRelation_untouched()
+    {
+        const QUuid src = m_model->addNode(QStringLiteral("source"), {0, 0});
+        const QUuid c = m_model->addNode(QStringLiteral("consumer"), {200, 0});
+        m_model->addEdge(src, "out", c, "in");
+        const FlowResult r = FlowSolver::solve(*m_model, QStringLiteral("power"));
+        QVERIFY(r.ok);
+        QVERIFY(r.edgeFlow.isEmpty());
+    }
+
+    void apply_writesFlowsAndSignals()
+    {
+        const QUuid src = m_model->addNode(QStringLiteral("source"), {0, 0});
+        const QUuid c = m_model->addNode(QStringLiteral("consumer"), {200, 0});
+        const QUuid e = m_model->addEdge(src, "out", c, "in");
+        QSignalSpy spy(m_model, &NetworkModel::flowChanged);
+        FlowSolver::apply(*m_model, FlowSolver::solve(*m_model, QStringLiteral("water")));
+        QCOMPARE(spy.count(), 1);
+        QVERIFY(m_model->edge(e)->flowRate > 0);
+    }
+
+private:
+    TypeRegistry m_reg;
+    NetworkModel* m_model = nullptr;
+};
+
+QTEST_GUILESS_MAIN(TestSolver)
+#include "tst_solver.moc"