perf(533-2 2차): E04 탐색 프로세스 병렬 · 솔버 수렴 묶음 빼기 · 창 배열 덜 펴기 — 1 km 띠 287 s → 32 s

- Search: 후보 덩이를 프로세스로 나눠 셈(workers · 0 = 자동 · 격자 5만 칸 이상이면 CPU − 2 · 12 까지) · 낸 차례대로 합침 — 일꾼 수와 무관하게 같은 값
- Search: 지층은 임시 파일로 일꾼에 넘김(spawn 관이면 일꾼이 차례로 떠 12 개 6 s) · 돌려받는 칸은 iy · ix · group(int32)만
- Search: 고친 노드 판정 = 반지름 시작 때 셀 최소(덩이 나눔과 무관) · 1차 SearchParams 칸 · SearchResult 꼴 그대로
- Search: 1단계 거르기 = 활동면 높이를 창 전체가 아니라 걸친 칸에서만 · 큰 창 배열 펴지 않음
- Column: Bishop · Janbu 고정점 — 기둥별 상수 미리 · 수렴한 활동체 기둥은 빼며 셈(값 같음)
- 측정(검증_sub5 57 측점 · 노드 2 m · 구 + 타원체 2): 1차 287 s → 한 프로세스 217 s → 12 프로세스 32.1 s · 셀 9.7만 F 전후 최대 차 6e-15 · 전역 최소 1.099 그대로
- 시험: 일꾼 1 · 2 같은 값

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

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
 

E04_Stab3D/E04_Stab3D_Engine_Search.py

@@ -9,6 +9,7 @@
   구는 활동체 지표 평균 낙하 방향(Scoops3D 2.1).
 - 한 구가 지표를 여러 덩이로 자르면 덩이마다 셈(ndimage.label · 8 이웃).
 - 격자 가장자리 · DEM 빈칸에 닿은 덩이는 버림(잘린 질량 막음).
+- 후보 덩이는 프로세스로 나눠 셈 · 낸 차례대로 합침 — 일꾼 수와 무관하게 같은 값.
 
 다스리는 안전율(governing):
 - 못 뚫는 층(기본 발파암 = 연암)에 잘린 바닥 면적 비 ≥ 10 % → Janbu 보정으로 전환
@@ -21,8 +22,14 @@
 from __future__ import annotations
 
 import math
+import os
+import pickle
+import tempfile
 import time
+import weakref
+from concurrent.futures import ProcessPoolExecutor
 from dataclasses import dataclass, field
+from types import SimpleNamespace
 
 import numpy as np
 from scipy import ndimage
@@ -56,6 +63,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
@@ -103,13 +111,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:
@@ -126,13 +144,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
@@ -142,9 +163,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 점
@@ -207,6 +226,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)
@@ -276,6 +331,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()
@@ -297,27 +360,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(
@@ -331,7 +387,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)))
@@ -396,6 +453,8 @@ def search(strata: Strata, params: SearchParams | None = None) -> SearchResult:
     if ell and len(crit):
         run(crit, ell)
 
+    if pool:
+        pool.shutdown()
     f_out = np.where(np.isfinite(f_min), f_min, np.nan)
     return SearchResult(
         f_min=f_out,
@@ -410,19 +469,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

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