From c2501ff04e411099fc06de472bbb009168ea0536 Mon Sep 17 00:00:00 2001 From: Ilya Ashikhmin Date: Thu, 2 Jul 2026 23:14:11 +0200 Subject: [PATCH] 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 --- src/CMakeLists.txt | 2 + src/core/FlowSolver.cpp | 177 ++++++++++++++++++++++++++++++++++++++++ src/core/FlowSolver.h | 29 +++++++ tests/CMakeLists.txt | 1 + tests/tst_solver.cpp | 123 ++++++++++++++++++++++++++++ 5 files changed, 332 insertions(+) create mode 100644 src/core/FlowSolver.cpp create mode 100644 src/core/FlowSolver.h create mode 100644 tests/tst_solver.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 7415783..5933ef3 100644 --- a/src/CMakeLists.txt +++ b/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) diff --git a/src/core/FlowSolver.cpp b/src/core/FlowSolver.cpp new file mode 100644 index 0000000..4e4cd9c --- /dev/null +++ b/src/core/FlowSolver.cpp @@ -0,0 +1,177 @@ +#include "FlowSolver.h" + +#include "NetworkModel.h" + +#include +#include + +#include +#include + +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>& a, QVector& b, QVector& 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 edges; + QList nodes; + QHash 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 parent(nodes.size()); + for (int i = 0; i < parent.size(); ++i) + parent[i] = i; + std::function 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> components; // root -> node indices + for (int i = 0; i < nodes.size(); ++i) + components[find(i)].append(i); + + QHash potentials; + for (auto it = components.begin(); it != components.end(); ++it) { + const QList& comp = it.value(); + const int n = comp.size(); + QHash local; // global node index -> local index + for (int i = 0; i < n; ++i) + local.insert(comp.at(i), i); + + QVector 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> a(m, QVector(m, 0.0)); + QVector 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 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 diff --git a/src/core/FlowSolver.h b/src/core/FlowSolver.h new file mode 100644 index 0000000..ec3ae8b --- /dev/null +++ b/src/core/FlowSolver.h @@ -0,0 +1,29 @@ +#pragma once + +#include "Types.h" + +#include + +namespace diag { + +class NetworkModel; + +struct FlowResult { + bool ok = false; + QString message; + QHash 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 diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 07a3064..08b1540 100644 --- a/tests/CMakeLists.txt +++ b/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) diff --git a/tests/tst_solver.cpp b/tests/tst_solver.cpp new file mode 100644 index 0000000..5a4ce96 --- /dev/null +++ b/tests/tst_solver.cpp @@ -0,0 +1,123 @@ +#include "TestFixtures.h" +#include "core/FlowSolver.h" +#include "core/NetworkModel.h" + +#include +#include + +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"