257 lines
13 KiB
Python
257 lines
13 KiB
Python
#!/usr/bin/env python3
|
|
"""Screen alternative burst pipes near existing pressure meters using a frozen INP."""
|
|
from __future__ import annotations
|
|
|
|
import argparse
|
|
import asyncio
|
|
import heapq
|
|
import json
|
|
import math
|
|
import sys
|
|
import os
|
|
from collections import defaultdict, deque
|
|
from pathlib import Path
|
|
from uuid import uuid4
|
|
|
|
|
|
def nearest_pressure_distances(
|
|
graph: dict[str, list[tuple[str, float]]], devices: list[dict], radius_m: float,
|
|
) -> dict[str, tuple[float, str]]:
|
|
"""Measure to the nearest pressure meter over the network, not just a selected meter."""
|
|
distances = {}
|
|
pending = []
|
|
for device in devices:
|
|
node, device_id = device["node_id"], device["device_id"]
|
|
value = (0.0, device_id)
|
|
if value < distances.get(node, (math.inf, "")):
|
|
distances[node] = value
|
|
heapq.heappush(pending, (0.0, device_id, node))
|
|
while pending:
|
|
distance, device_id, node = heapq.heappop(pending)
|
|
if distances[node] != (distance, device_id):
|
|
continue
|
|
for neighbor, length in graph.get(node, []):
|
|
cost = distance + length
|
|
value = (cost, device_id)
|
|
if cost <= radius_m and value < distances.get(neighbor, (math.inf, "")):
|
|
distances[neighbor] = value
|
|
heapq.heappush(pending, (cost, device_id, neighbor))
|
|
return distances
|
|
|
|
|
|
def main() -> int:
|
|
parser = argparse.ArgumentParser(description=__doc__)
|
|
parser.add_argument("--calibration-dir", required=True, type=Path)
|
|
parser.add_argument("--limit", type=int, default=24)
|
|
parser.add_argument("--offset", type=int, default=0)
|
|
parser.add_argument("--radius-m", type=float, default=400)
|
|
parser.add_argument("--min-network-distance-m", type=float, default=0)
|
|
parser.add_argument("--min-diameter", type=float, default=150)
|
|
parser.add_argument("--require-isolation", action="store_true")
|
|
parser.add_argument("--exclude-existing-km", type=float, default=0)
|
|
parser.add_argument("--skip-screening", type=Path)
|
|
args = parser.parse_args()
|
|
os.umask(0o077)
|
|
sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
|
|
from scripts.benchmark.calibration import configure_inp, digest, run_file, save_json
|
|
from scripts.benchmark.simulation_design import BurstSimulationDesign
|
|
from scripts.check_burst_simulation_design import readonly, resolve_project
|
|
|
|
if (not 1 <= args.limit <= 100 or args.offset < 0 or not 100 <= args.radius_m <= 5000
|
|
or not 0 <= args.min_network_distance_m < args.radius_m
|
|
or not 50 <= args.min_diameter <= 800 or not 0 <= args.exclude_existing_km <= 10):
|
|
parser.error("invalid screen limit, offset, radius, diameter or exclusion distance")
|
|
folder = args.calibration_dir.resolve()
|
|
calibration = json.loads((folder / "calibration.json").read_text())
|
|
design = BurstSimulationDesign.model_validate(calibration["design"])
|
|
snapshot = folder / "source.inp"
|
|
if digest(snapshot) != calibration["model_sha256"]:
|
|
raise ValueError("model snapshot fingerprint changed")
|
|
biz, _ = asyncio.run(resolve_project(design))
|
|
with readonly(biz) as connection:
|
|
devices = connection.execute(
|
|
"SELECT device_id,device_type,node_id,link_id FROM asset.scada_devices "
|
|
"WHERE device_type IN ('pressure','pipe_flow') ORDER BY device_id"
|
|
).fetchall()
|
|
pressure = [device for device in devices if device["device_type"] == "pressure"]
|
|
flow = [device for device in devices if device["device_type"] == "pipe_flow" and device["device_id"] != "547063"]
|
|
source = snapshot.read_text()
|
|
graph = defaultdict(list)
|
|
pipes = {}
|
|
valve_links = {}
|
|
source_nodes = set()
|
|
closed_links = set()
|
|
junctions = set()
|
|
existing_emitters = set()
|
|
section = ""
|
|
for raw in source.splitlines():
|
|
line = raw.split(";", 1)[0].strip()
|
|
if line.startswith("["):
|
|
section = line.strip("[]").upper()
|
|
continue
|
|
tokens = line.split()
|
|
if not tokens:
|
|
continue
|
|
if section == "JUNCTIONS":
|
|
junctions.add(tokens[0])
|
|
elif section in {"RESERVOIRS", "TANKS"}:
|
|
source_nodes.add(tokens[0])
|
|
elif section == "STATUS" and tokens[-1].upper() == "CLOSED":
|
|
closed_links.add(tokens[0])
|
|
elif section == "EMITTERS":
|
|
existing_emitters.add(tokens[0])
|
|
elif section in {"PIPES", "VALVES"}:
|
|
link, left, right = tokens[:3]
|
|
length = float(tokens[3]) if section == "PIPES" else 0.01
|
|
if section == "PIPES":
|
|
if tokens[-1].upper() == "CLOSED":
|
|
continue
|
|
pipes[link] = {"pipe_id": link, "start_node_id": left, "emitter_node_id": right,
|
|
"length_m": length, "diameter_mm": float(tokens[4])}
|
|
else:
|
|
valve_links[link] = (left, right)
|
|
graph[left].append((right, max(length, 0.01)))
|
|
graph[right].append((left, max(length, 0.01)))
|
|
# Initial link overrides also apply to pipe and valve connectivity.
|
|
graph.clear()
|
|
pipe_graph = defaultdict(list)
|
|
for link, pipe in pipes.items():
|
|
if link in closed_links:
|
|
continue
|
|
left, right = pipe["start_node_id"], pipe["emitter_node_id"]
|
|
length = max(pipe["length_m"], 0.01)
|
|
graph[left].append((right, length))
|
|
graph[right].append((left, length))
|
|
pipe_graph[left].append(right)
|
|
pipe_graph[right].append(left)
|
|
for link, (left, right) in valve_links.items():
|
|
if link not in closed_links:
|
|
graph[left].append((right, 0.01))
|
|
graph[right].append((left, 0.01))
|
|
component_by_node = {}
|
|
component_data = {}
|
|
|
|
def isolation(node):
|
|
if node not in component_by_node:
|
|
seen = {node}
|
|
pending = deque([node])
|
|
while pending:
|
|
for neighbor in pipe_graph[pending.popleft()]:
|
|
if neighbor not in seen:
|
|
seen.add(neighbor)
|
|
pending.append(neighbor)
|
|
component_id = min(seen)
|
|
for member in seen:
|
|
component_by_node[member] = component_id
|
|
boundaries = sorted(link for link, (left, right) in valve_links.items()
|
|
if link not in closed_links and (left in seen) != (right in seen))
|
|
component_data[component_id] = {"region_id": component_id, "node_count": len(seen),
|
|
"boundary_valve_ids": boundaries,
|
|
"source_inside": bool(seen & source_nodes)}
|
|
return component_data[component_by_node[node]]
|
|
# Only real pressure-meter model nodes present in this frozen snapshot are used.
|
|
pressure = [device for device in pressure if device["node_id"] in junctions]
|
|
flow = [device for device in flow if device["link_id"] in pipes]
|
|
distances = nearest_pressure_distances(graph, pressure, args.radius_m)
|
|
candidates = []
|
|
for pipe in pipes.values():
|
|
node = pipe["emitter_node_id"]
|
|
if (args.min_diameter <= pipe["diameter_mm"] <= 800 and pipe["length_m"] >= 50
|
|
and pipe["pipe_id"] not in closed_links
|
|
and node in junctions and node not in existing_emitters and node in distances):
|
|
distance, device = distances[node]
|
|
if distance < args.min_network_distance_m:
|
|
continue
|
|
region = isolation(node)
|
|
if args.require_isolation and (region["source_inside"] or not 1 <= len(region["boundary_valve_ids"]) <= 12):
|
|
continue
|
|
candidates.append({**pipe, "nearest_pressure_device": device, "network_distance_m": distance,
|
|
"screen_priority": pipe["length_m"] / pipe["diameter_mm"] ** 2,
|
|
"isolation": region})
|
|
skipped = set()
|
|
if args.skip_screening:
|
|
skipped = {entry["pipe_id"] for entry in json.loads(args.skip_screening.read_text())["results"]}
|
|
with readonly(biz) as connection:
|
|
coords = {row["node_id"]: row for row in connection.execute(
|
|
"SELECT node_id,ST_X(ST_Transform(geom,4326)) longitude,"
|
|
"ST_Y(ST_Transform(geom,4326)) latitude FROM gis.node_geometries WHERE node_id=ANY(%s)",
|
|
(list({pipe["emitter_node_id"] for pipe in candidates}),),
|
|
)}
|
|
existing_positions = connection.execute(
|
|
"SELECT ST_X(ST_Transform(geom,4326)) longitude,ST_Y(ST_Transform(geom,4326)) latitude "
|
|
"FROM gis.node_geometries WHERE node_id=ANY(%s)",
|
|
([pipe.emitter_node_id for pipe in design.pipes],),
|
|
).fetchall()
|
|
filtered = []
|
|
for candidate in candidates:
|
|
position = coords.get(candidate["emitter_node_id"])
|
|
if position is None or candidate["pipe_id"] in skipped:
|
|
continue
|
|
distances_km = [111.32 * math.hypot((position["longitude"] - old["longitude"]) *
|
|
math.cos(math.radians(position["latitude"])),
|
|
position["latitude"] - old["latitude"])
|
|
for old in existing_positions]
|
|
if distances_km and min(distances_km) < args.exclude_existing_km:
|
|
continue
|
|
candidate.update(longitude=position["longitude"], latitude=position["latitude"],
|
|
nearest_existing_km=min(distances_km) if distances_km else None)
|
|
filtered.append(candidate)
|
|
candidates = filtered
|
|
# Cover different meter neighborhoods before taking additional pipes at the same meter.
|
|
candidates.sort(key=lambda pipe: pipe["screen_priority"], reverse=True)
|
|
used = set()
|
|
unique, remaining = [], []
|
|
for candidate in candidates:
|
|
device = candidate["nearest_pressure_device"]
|
|
(remaining if device in used else unique).append(candidate)
|
|
used.add(device)
|
|
selected = (unique + remaining)[args.offset:args.offset + args.limit]
|
|
node_ids = sorted({device["node_id"] for device in pressure} | {pipe["emitter_node_id"] for pipe in selected})
|
|
link_ids = sorted({device["link_id"] for device in flow})
|
|
at = max(design.event_times)
|
|
output = folder / f"screening-{uuid4().hex[:8]}"
|
|
output.mkdir(mode=0o700)
|
|
report = {"version": "burst-pipe-screening@1", "model_sha256": digest(snapshot),
|
|
"at": at.isoformat(), "candidate_count": len(candidates), "results": [],
|
|
"screen_options": {key: str(value) if isinstance(value, Path) else value for key, value in vars(args).items()},
|
|
"device_bindings": devices,
|
|
"scope": "daytime single-frame screening only; requires full paired calibration and isolation review"}
|
|
|
|
def simulate(label, emitter=None):
|
|
path = output / f"{label}.inp"
|
|
path.write_text(configure_inp(source, at, 0, emitter=emitter, pda_options=calibration["pda_options"]))
|
|
data = run_file(path, node_ids, link_ids, timeout=240)
|
|
save_json(output / f"{label}.json", data)
|
|
return data
|
|
|
|
normal = simulate("normal")
|
|
if normal["status"] != "completed" or normal["warnings"]:
|
|
raise ValueError("normal screening baseline failed")
|
|
for index, pipe in enumerate(selected):
|
|
area = math.pi * (pipe["diameter_mm"] / 10) ** 2 / 4 / 8
|
|
burst = simulate(f"candidate-{index:02d}", (pipe["emitter_node_id"], area))
|
|
result = {**pipe, "orifice_area_cm2": area, "status": burst["status"], "output_file": f"candidate-{index:02d}.json"}
|
|
if burst["status"] == "completed":
|
|
drops = {device["device_id"]: 100 * (normal["nodes"][device["node_id"]][0]["pressure"] -
|
|
burst["nodes"][device["node_id"]][0]["pressure"]) /
|
|
max(abs(normal["nodes"][device["node_id"]][0]["pressure"]), 1e-6) for device in pressure}
|
|
changes = {device["device_id"]: 100 * abs(normal["links"][device["link_id"]][0]["flow"] -
|
|
burst["links"][device["link_id"]][0]["flow"]) /
|
|
abs(normal["links"][device["link_id"]][0]["flow"]) for device in flow
|
|
if abs(normal["links"][device["link_id"]][0]["flow"]) > 1e-6}
|
|
ranked = sorted(drops.items(), key=lambda pair: pair[1], reverse=True)
|
|
negative = any(row[0]["pressure"] < 0 for row in burst["nodes"].values())
|
|
result.update(pressure_drop_percent=dict(ranked), flow_change_percent=changes,
|
|
possible_candidate=not burst["warnings"] and not negative and ranked[0][1] >= 30
|
|
and (ranked[1][1] >= 10 or max(changes.values(), default=0) >= 20))
|
|
print(f"{index+1}/{len(selected)} {pipe['pipe_id']}: max pressure drop={ranked[0][1]:.2f}%, "
|
|
f"possible candidate={result['possible_candidate']}", flush=True)
|
|
report["results"].append(result)
|
|
save_json(output / "screening.json", report)
|
|
print(f"screening complete; report: {output / 'screening.json'}", flush=True)
|
|
return 0
|
|
|
|
|
|
if __name__ == "__main__":
|
|
raise SystemExit(main())
|