diff --git a/E04_Stab3D/E04_Stab3D_Engine_Column.py b/E04_Stab3D/E04_Stab3D_Engine_Column.py index 2b901a17d..4402df635 100644 --- a/E04_Stab3D/E04_Stab3D_Engine_Column.py +++ b/E04_Stab3D/E04_Stab3D_Engine_Column.py @@ -244,6 +244,27 @@ def _finish(res: np.ndarray, drive: np.ndarray) -> np.ndarray: return out +def _fixed_point(st: dict, f: np.ndarray, coef: dict, step) -> tuple[np.ndarray, np.ndarray]: + """F = step(기둥 묶음, F) 고정점 — 수렴한 활동체의 기둥은 빼며 셈(값은 빼지 않을 때와 같음). + + coef = 기둥별 상수 배열(group 포함) · 반환 (F, 수렴 여부). + """ + n = st["n"] + sub = coef + conv = np.zeros(n, bool) + for _ in range(MAX_ITER): + f_new = step(sub, f) + done = np.abs(f_new - f) < TOL + f = np.where(conv, f, f_new) + conv |= done + if conv.all(): + break + live = ~conv[sub["group"]] + if live.mean() < 0.6: # 남은 기둥만 — 매번 줄이면 줄이는 값이 더 듦 + sub = {k: v[live] for k, v in sub.items()} + return f, conv + + def _bishop(st: dict, f0: np.ndarray) -> tuple[np.ndarray, np.ndarray]: grp, n = st["group"], st["n"] wq = st["w"] + st["q"] @@ -252,19 +273,29 @@ def _bishop(st: dict, f0: np.ndarray) -> tuple[np.ndarray, np.ndarray]: lift = np.zeros(n) np.maximum.at(lift, grp, neg) f = np.where(np.isfinite(f0), np.maximum(f0, 0.05), 1.0) + lift - conv = np.zeros(n, bool) - for _ in range(MAX_ITER): - s_res, _ = _s_res(st, f) - s = s_res / np.maximum(f[grp], 1e-3) - nn = (wq - s * st["sin_a"]) / st["cos_psi"] # 연직 평형에서 N - res = _gsum(grp, st["arm_s"] * s_res, n) - drive = _gsum(grp, -st["p"] * wq + st["e"] * st["kh"] * st["w"] + nn * st["m_n"], n) - f_new = _finish(res, drive) - done = np.abs(f_new - f) < TOL - f = np.where(conv, f, f_new) - conv |= done - if conv.all(): - break + # S_res = num / mα · N = (W+Q − S sinα)/cosψ → 구동 = d1 − Σ S_res·k_n / F (기둥별 상수 미리) + num = st["c"] * st["a_h"] + np.maximum(wq - st["u_ah"], 0.0) * st["tan_phi"] + d1 = _gsum( + grp, -st["p"] * wq + st["e"] * st["kh"] * st["w"] + wq * st["m_n"] / st["cos_psi"], n + ) + coef = { + "group": grp, + "num": num, + "cos_psi": st["cos_psi"], + "k": st["sin_a"] * st["tan_phi"], + "arm_s": st["arm_s"], + "k_n": st["sin_a"] * st["m_n"] / st["cos_psi"], + } + + def step(c: dict, f: np.ndarray) -> np.ndarray: + g = c["group"] + fg = np.maximum(f, 1e-3) + s_res = c["num"] / (c["cos_psi"] + c["k"] / fg[g]) + res = _gsum(g, c["arm_s"] * s_res, n) + drive = d1 - _gsum(g, s_res * c["k_n"], n) / fg + return _finish(res, drive) + + f, conv = _fixed_point(st, f, coef, step) _, m_a = _s_res(st, f) ma_min = _gmin(grp, m_a, n) # 못 쓰는 면 — 수렴 안 함 · mα ≤ 0(바닥 수직력이 음) · F ≤ 0 → nan(후보에서 빠짐) @@ -280,17 +311,23 @@ def _gmin(grp: np.ndarray, values: np.ndarray, n: int) -> np.ndarray: def _janbu(st: dict, f0: np.ndarray) -> np.ndarray: grp, n = st["group"], st["n"] - drive = _gsum(grp, (st["w"] + st["q"]) * st["tan_a"] + st["kh"] * st["w"], n) + wq = st["w"] + st["q"] + drive = _gsum(grp, wq * st["tan_a"] + st["kh"] * st["w"], n) f = np.where(np.isfinite(f0), np.maximum(f0, 0.05), 1.0) - conv = np.zeros(n, bool) - for _ in range(MAX_ITER): - s_res, _ = _s_res(st, f) - f_new = _finish(_gsum(grp, s_res / st["cos_a"], n), drive) - done = np.abs(f_new - f) < TOL - f = np.where(conv, f, f_new) - conv |= done - if conv.all(): - break + num = st["c"] * st["a_h"] + np.maximum(wq - st["u_ah"], 0.0) * st["tan_phi"] + coef = { + "group": grp, + "num": num / st["cos_a"], + "cos_psi": st["cos_psi"], + "k": st["sin_a"] * st["tan_phi"], + } + + def step(c: dict, f: np.ndarray) -> np.ndarray: + g = c["group"] + s_res = c["num"] / (c["cos_psi"] + c["k"] / np.maximum(f, 1e-3)[g]) + return _finish(_gsum(g, s_res, n), drive) + + f, conv = _fixed_point(st, f, coef, step) f[~conv | (f <= 0)] = np.nan return f diff --git a/E04_Stab3D/E04_Stab3D_Engine_Search.py b/E04_Stab3D/E04_Stab3D_Engine_Search.py index 1716245f1..2dda11d70 100644 --- a/E04_Stab3D/E04_Stab3D_Engine_Search.py +++ b/E04_Stab3D/E04_Stab3D_Engine_Search.py @@ -9,6 +9,7 @@ 구는 활동체 지표 평균 낙하 방향(Scoops3D 2.1). - 한 구가 지표를 여러 덩이로 자르면 덩이마다 셈(ndimage.label · 8 이웃). - 격자 가장자리 · DEM 빈칸에 닿은 덩이는 버림(잘린 질량 막음). +- 후보 덩이는 프로세스로 나눠 셈 · 낸 차례대로 합침 — 일꾼 수와 무관하게 같은 값. 다스리는 안전율(governing): - 못 뚫는 층(기본 발파암 = 연암)에 잘린 바닥 면적 비 ≥ 10 % → Janbu 보정으로 전환 @@ -21,9 +22,15 @@ from __future__ import annotations import math +import os +import pickle +import tempfile import time +import weakref from collections.abc import Callable +from concurrent.futures import ProcessPoolExecutor from dataclasses import dataclass, field +from types import SimpleNamespace import numpy as np from scipy import ndimage @@ -57,6 +64,7 @@ class SearchParams: a_max: float = math.inf ratios: tuple[float, ...] = (1.0, 2.0) # 구(늘 봄) + 타원체 축 방향 배수(구 임계 노드에서) sub: int = 2 # 셀 안 표본 sub×sub(부분 기둥) + workers: int = 0 # 프로세스 수 · 0 = 자동(격자 5만 칸 이상이면 CPU − 2 · 12 까지) refine_rounds: int = 4 # 가장 고운 간격에서 더 내려가 보는 횟수 상한 # 고운 탐색 · 타원체는 F < max(이 값, 전역 최소 × 1.2) 로 고친 노드만(높은 곳 지도는 거친 값) refine_below: float = 2.0 @@ -107,13 +115,23 @@ def trial_columns( ixg = ix0[:, None, None] + off[None, None, :] iyg = iy0[:, None, None] + off[None, :, None] inb = (ixg >= 0) & (ixg < nx) & (iyg >= 0) & (iyg < ny) - ixc, iyc, inb = np.broadcast_arrays(np.clip(ixg, 0, nx - 1), np.clip(iyg, 0, ny - 1), inb) - shape = ixc.shape + nw = 2 * w + 1 + shape = (n_t, nw, nw) d = None if ratio == 1.0 else np.asarray(directions, float) k2 = 1.0 / (ratio * ratio) + def remainder(t, xs, ys): + """r² − (중심에서 수평 거리 꼴)² — 0 보다 크면 활동면이 그 점 위를 지남.""" + dxp = xs - centers[t, 0] + dyp = ys - centers[t, 1] + if d is None: + return r * r - dxp * dxp - dyp * dyp + pp = dxp * d[t, 0] + dyp * d[t, 1] + qq = -dxp * d[t, 1] + dyp * d[t, 0] + return r * r - pp * pp - k2 * qq * qq + def surface(t, xs, ys): - """후보 t 의 활동면 — (rem, z_s, ∂z/∂x, ∂z/∂y). 원소별(납작 · 창 모두).""" + """후보 t 의 활동면 — (rem, z_s, ∂z/∂x, ∂z/∂y). 원소별(납작).""" dxp = xs - centers[t, 0] dyp = ys - centers[t, 1] if d is None: @@ -130,13 +148,16 @@ def trial_columns( return rem, centers[t, 2] - root, gx, gy # 1) 셀 중심으로 거르기 — 걸친 셀과 그 이웃만 셀 안 표본을 봄(창 대부분은 허공 · 땅속 밖) - t_grid = np.broadcast_to(np.arange(n_t)[:, None, None], shape) - xs0 = strata.x[0] + ixg * dx - ys0 = strata.y[0] + iyg * dx - rem0, zs0, _, _ = surface(t_grid, xs0, ys0) + # 창 칸 번호 → (후보, 셀) — 큰 창 배열을 펴지 않고 셈 + def cell_of(flat): + t = flat // (nw * nw) + return t, iy0[t] + off[(flat // nw) % nw], ix0[t] + off[flat % nw] + + rem0 = remainder(np.arange(n_t)[:, None, None], strata.x[0] + ixg * dx, strata.y[0] + iyg * dx) in0 = np.flatnonzero(inb & (rem0 > 0)) - yy0, xx0 = iyc.ravel()[in0], ixc.ravel()[in0] - zb0 = np.maximum(zs0.ravel()[in0], strata.hard_top[yy0, xx0]) + t0, yy0, xx0 = cell_of(in0) + zs0 = centers[t0, 2] - np.sqrt(rem0.ravel()[in0]) + zb0 = np.maximum(zs0, strata.hard_top[yy0, xx0]) act0 = np.zeros(shape, bool) act0.ravel()[in0[zb0 < strata.z_top[yy0, xx0]]] = True cand = act0 if sub == 1 else ndimage.binary_dilation(act0, structure=_STRUCT) & inb @@ -146,9 +167,7 @@ def trial_columns( empty = (Columns(e, e, e, np.zeros(0), np.zeros(0), np.zeros(0), np.zeros((0, 3))), e, big) if not len(ci): return empty - tt = t_grid.ravel()[ci] - cy = iyc.ravel()[ci] - cx = ixc.ravel()[ci] + tt, cy, cx = cell_of(ci) hard = strata.hard_top[cy, cx] # 2) 부분 기둥 — sub×sub 점 @@ -211,6 +230,42 @@ def trial_columns( return cols, trial, big +_PART_OPTS = ("v_min", "v_max", "a_min", "a_max", "sub", "groundwater", "seismic_kh") +_PART: dict = {} # 일꾼 프로세스의 지층 · 조건 + + +def _part_init(strata: Strata | str, opts: dict) -> None: + if isinstance(strata, str): + with open(strata, "rb") as fh: + strata = pickle.load(fh) + _PART.update(strata=strata, **opts) + + +def _open_pool(strata: Strata, opts: dict, n_workers: int) -> ProcessPoolExecutor: + """일꾼 풀 — 지층은 임시 파일로 넘김(spawn 관으로 넘기면 일꾼이 차례로 뜸 · 12 개 6 s).""" + fd, path = tempfile.mkstemp(suffix=".e04.pkl") + with os.fdopen(fd, "wb") as fh: + pickle.dump(strata, fh, protocol=pickle.HIGHEST_PROTOCOL) + pool = ProcessPoolExecutor(n_workers, initializer=_part_init, initargs=(path, opts)) + weakref.finalize(pool, os.remove, path) # 멈춤 · 오류로 빠져나가도 지움 + return pool + + +def _part(job: tuple, k: dict | None = None) -> tuple: + """후보 덩이 하나 — 기둥 + 솔버. 셀 합치기에 쓰는 칸(iy · ix · group)만 돌려줌(옮길 양 줄임).""" + xyz, dirs, r, ratio = job + k = k or _PART + st = k["strata"] + cols, trial, big = trial_columns( + st, xyz, r, ratio, dirs, k["v_min"], k["v_max"], k["a_min"], k["a_max"], k["sub"] + ) + res = solve(st, cols, k["groundwater"], k["seismic_kh"]) if len(trial) else None + light = SimpleNamespace( + iy=cols.iy.astype(np.int32), ix=cols.ix.astype(np.int32), group=cols.group.astype(np.int32) + ) + return light, trial, big, res + + def governing(res: dict) -> tuple[np.ndarray, np.ndarray]: """다스리는 안전율 · 방법(0 Bishop · 1 Janbu 보정) — 머리 주석 규칙.""" jc = np.where(np.isfinite(res["janbu_corrected"]), res["janbu_corrected"], np.inf) @@ -287,6 +342,14 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult: z = strata.z_top[cy, cx] + p.h_min + ih * p.h_step return np.column_stack([x, y, z]), cy, cx, p.h_min + ih * p.h_step + # 후보 덩이 = 프로세스로 나눠 셈(_part) · 결과는 낸 차례대로 합침 → 일꾼 수와 무관한 같은 값 + n_workers = p.workers or ( + min(12, max(1, (os.cpu_count() or 2) - 2)) if ny * nx >= 50_000 else 1 + ) + opts = {k: getattr(p, k) for k in _PART_OPTS} + pool = _open_pool(strata, opts, n_workers) if n_workers > 1 else None + local = {"strata": strata, **opts} + def run(ids: np.ndarray, ratios: tuple[float, ...] = (1.0,)) -> set[int]: """노드들의 모든 반지름 · 비 — 셀 최소를 고친 노드 번호를 돌려줌.""" improved: set[int] = set() @@ -310,27 +373,20 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult: if not len(live): continue w_now = (2 * math.ceil(r * max(ratio, 1) / strata.dx) + 3) ** 2 - chunk = max(1, BATCH_CELLS // max(w_now, 1)) - for s in range(0, len(live), chunk): - part = live[s : s + chunk] - d_part = dirs[part] if ratio != 1.0 else None - cols, trial, big = trial_columns( - strata, - xyz[part], - r, - ratio, - d_part, - p.v_min, - p.v_max, - p.a_min, - p.a_max, - p.sub, - ) + # 메모리 = 일꾼 수 합쳐 BATCH_CELLS · 일꾼마다 한 덩이 이상(작은 고운 탐색도 나눔) + chunk = max(1, BATCH_CELLS // (w_now * n_workers)) + chunk = min(chunk, -(-len(live) // n_workers)) + parts = [live[s : s + chunk] for s in range(0, len(live), chunk)] + jobs = [(xyz[q], dirs[q] if ratio != 1.0 else None, r, ratio) for q in parts] + outs = pool.map(_part, jobs) if pool else (_part(j, local) for j in jobs) + # 고친 노드 판정 기준 = 이 반지름 시작 값(덩이 나눔과 무관) + ref = f_min.copy() + thr = max(p.refine_below, 1.2 * best["governing"]) + for part, (cols, trial, big, res) in zip(parts, outs): done[part[big > p.v_max]] = True counters["trials"] += len(part) if not len(trial): continue - res = solve(strata, cols, p.groundwater, p.seismic_kh) counters["evaluated"] += len(trial) f, meth = governing(res) improved |= _update( @@ -344,7 +400,8 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult: volume, method, node_map, - max(p.refine_below, 1.2 * best["governing"]), + thr, + ref, ) k = ( int(np.nanargmin(np.where(np.isfinite(f), f, np.inf))) @@ -368,50 +425,57 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult: ) return improved - step = max(1, int(p.coarse)) - grid = np.stack( - np.meshgrid( - np.arange(0, nhs, step), np.arange(0, nys, step), np.arange(0, nxs, step), indexing="ij" - ), - -1, - ).reshape(-1, 3) - ids = (grid[:, 0] * nys + grid[:, 1]) * nxs + grid[:, 2] - seen.update(ids.tolist()) - seeds = run(ids) - rounds = passes = 0 - while seeds and (step > 1 or rounds < p.refine_rounds): - if step == 1: - rounds += 1 - step = max(1, step // 2) - passes += 1 - stage[:] = [min(40 + 8 * (passes - 1), 77), 8, f"고운 탐색 {passes}"] - cand = set() - for nid in seeds: - ih, rem = divmod(nid, nys * nxs) - iy, ix = divmod(rem, nxs) - for dh in (-step, 0, step): - for dy in (-step, 0, step): - for dxs in (-step, 0, step): - h2, y2, x2 = ih + dh, iy + dy, ix + dxs - if 0 <= h2 < nhs and 0 <= y2 < nys and 0 <= x2 < nxs: - cid = (h2 * nys + y2) * nxs + x2 - if cid not in seen: - cand.add(cid) - if not cand: + try: + step = max(1, int(p.coarse)) + grid = np.stack( + np.meshgrid( + np.arange(0, nhs, step), + np.arange(0, nys, step), + np.arange(0, nxs, step), + indexing="ij", + ), + -1, + ).reshape(-1, 3) + ids = (grid[:, 0] * nys + grid[:, 1]) * nxs + grid[:, 2] + seen.update(ids.tolist()) + seeds = run(ids) + rounds = passes = 0 + while seeds and (step > 1 or rounds < p.refine_rounds): if step == 1: - break - continue - seen.update(cand) - seeds = run(np.fromiter(cand, np.int64)) + rounds += 1 + step = max(1, step // 2) + passes += 1 + stage[:] = [min(40 + 8 * (passes - 1), 77), 8, f"고운 탐색 {passes}"] + cand = set() + for nid in seeds: + ih, rem = divmod(nid, nys * nxs) + iy, ix = divmod(rem, nxs) + for dh in (-step, 0, step): + for dy in (-step, 0, step): + for dxs in (-step, 0, step): + h2, y2, x2 = ih + dh, iy + dy, ix + dxs + if 0 <= h2 < nhs and 0 <= y2 < nys and 0 <= x2 < nxs: + cid = (h2 * nys + y2) * nxs + x2 + if cid not in seen: + cand.add(cid) + if not cand: + if step == 1: + break + continue + seen.update(cand) + seeds = run(np.fromiter(cand, np.int64)) - # 타원체 — 구 탐색이 찾은 임계 노드(F < refine_below)에서만(시간 1/3 · 노선 1 km 실측) - ell = tuple(r for r in p.ratios if r != 1.0) - limit = max(p.refine_below, 1.2 * best["governing"]) - crit = np.unique(node_map[(node_map >= 0) & (f_min < limit)]) - if ell and len(crit): - stage[:] = [85, 14, "타원체"] - run(crit, ell) + # 타원체 — 구 탐색이 찾은 임계 노드(F < refine_below)에서만(시간 1/3 · 노선 1 km 실측) + ell = tuple(r for r in p.ratios if r != 1.0) + limit = max(p.refine_below, 1.2 * best["governing"]) + crit = np.unique(node_map[(node_map >= 0) & (f_min < limit)]) + if ell and len(crit): + stage[:] = [85, 14, "타원체"] + run(crit, ell) + finally: + if pool: # 끊기(InterruptedError) · 오류에도 일꾼을 닫음 + pool.shutdown(cancel_futures=True) f_out = np.where(np.isfinite(f_min), f_min, np.nan) return SearchResult( f_min=f_out, @@ -426,19 +490,35 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult: def _update( - f, meth, res, cols, trial, node_ids, f_min, volume, method, node_map, refine_below=np.inf + f, + meth, + res, + cols, + trial, + node_ids, + f_min, + volume, + method, + node_map, + refine_below=np.inf, + ref=None, ) -> set[int]: - """셀마다 더 낮은 F 면 바꿈 — 같은 셀 여럿이면 가장 낮은 것(내림차순으로 써 마지막이 이김).""" + """셀마다 더 낮은 F 면 바꿈 — 같은 셀 여럿이면 가장 낮은 것(내림차순으로 써 마지막이 이김). + + 돌려주는 노드 = ref(없으면 f_min) 보다 낮춘 후보의 노드 중 F < refine_below. + """ fc = f[cols.group] ok = np.isfinite(fc) iy, ix, grp, fc = cols.iy[ok], cols.ix[ok], cols.group[ok], fc[ok] order = np.argsort(-fc, kind="stable") iy, ix, grp, fc = iy[order], ix[order], grp[order], fc[order] + ref = f_min if ref is None else ref + hit = (fc < ref[iy, ix] * (1 - 1e-6)) & (fc < refine_below) + seeds = set(np.unique(node_ids[trial[grp[hit]]]).tolist()) better = fc < f_min[iy, ix] * (1 - 1e-6) iy, ix, grp, fc = iy[better], ix[better], grp[better], fc[better] f_min[iy, ix] = fc volume[iy, ix] = res["volume"][grp] method[iy, ix] = meth[grp] - nid = node_ids[trial[grp]] - node_map[iy, ix] = nid - return set(np.unique(nid[fc < refine_below]).tolist()) + node_map[iy, ix] = node_ids[trial[grp]] + return seeds diff --git a/resources/tester/test_533_2_stab3d_search.py b/resources/tester/test_533_2_stab3d_search.py index 1e3c408bd..aed8858d9 100644 --- a/resources/tester/test_533_2_stab3d_search.py +++ b/resources/tester/test_533_2_stab3d_search.py @@ -85,3 +85,15 @@ def test_wet_cone_pore_pressure_ratio(): groundwater={"mode": "ru", "ru": 0.4}, ) assert search(st, p).best["governing"] == pytest.approx(1.626, rel=0.03) + + +def test_workers_same_result(): + """프로세스 나눔 — 일꾼 수와 무관하게 셀별 최소 F · 전역 최소가 같음.""" + st = _embankment(10 / np.tan(np.radians(30)), 19, 20, 8.64) + kw = dict(box=(-2, 10, -4, 4), h_min=4, h_max=20, spacing=2.0, coarse=2, max_depth=16) + kw.update(v_min=50, v_max=2e4, ratios=(1.0, 2.0)) + one = search(st, SearchParams(workers=1, **kw)) + two = search(st, SearchParams(workers=2, **kw)) + np.testing.assert_array_equal(one.f_min, two.f_min) + assert one.best["governing"] == two.best["governing"] + assert one.evaluated == two.evaluated