From a39804bd8fccd64c7638b3e035f14c40cbe2c1c3 Mon Sep 17 00:00:00 2001 From: umsangdon Date: Mon, 26 Oct 2026 03:39:51 +0900 Subject: [PATCH] =?UTF-8?q?perf(E03):=202D=20LEM=20=EC=A0=88=ED=8E=B8=20?= =?UTF-8?q?=C2=B7=20=EB=B0=98=EB=B3=B5=20=EA=B3=84=EC=82=B0=20=EC=A4=84?= =?UTF-8?q?=EC=9E=84=20=E2=80=94=20=EA=B2=B0=EA=B3=BC=20=EA=B0=99=EC=9D=8C?= =?UTF-8?q?=20(532-2=202=EC=B0=A8=20=E2=91=A0)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 지층 칸: 어느 띠에도 안 쓰인 빈 구간 칸을 잘라 절편 폭 22 → 3 - Bishop · Janbu 반복: 끝난 면은 빼고 셈(남은 면이 절반 아래면 묶음 줄임) - λ = 0 길에서 λ 칸 계산 뺌 · Spencer 뉴턴 상한 40 → 25(수렴 면은 모두 27 회 안 · 넘는 면은 끝내 안 풀림) - 검증_sub4 59 측점 전후 1648 값 차 0 % · 18 프로세스 벽시계 73.4 → 44.9 s Co-Authored-By: Claude Opus 5.5 --- E03_StabSection/E03_StabSection_Engine_LEM.py | 63 ++++++++++++------- 1 file changed, 41 insertions(+), 22 deletions(-) diff --git a/E03_StabSection/E03_StabSection_Engine_LEM.py b/E03_StabSection/E03_StabSection_Engine_LEM.py index b3b05c057..ec30ffc3b 100644 --- a/E03_StabSection/E03_StabSection_Engine_LEM.py +++ b/E03_StabSection/E03_StabSection_Engine_LEM.py @@ -104,6 +104,11 @@ class Prepared: self.bot_a, self.bot_b = (np.take_along_axis(bot[k], order, 1) for k in _halves(half)) self.top_a, self.top_b = (np.take_along_axis(top[k], order, 1) for k in _halves(half)) self.mat = np.take_along_axis(mat[:half], order, 1) + # 빈 구간은 정렬 뒤 끝으로 — 어느 띠에도 안 쓰인 끝 칸을 잘라 절편 계산 폭을 줄임(값 같음) + used = ~(np.isnan(self.bot_a) & np.isnan(self.bot_b)) + k = int(np.max(np.nonzero(used.any(0))[0], initial=0)) + 1 + for name in ("bot_a", "bot_b", "top_a", "top_b", "mat"): + setattr(self, name, getattr(self, name)[:, :k]) # 지표면 꺾은선(정확) — 띠 양 끝으로 늘린 맨 윗면 · 연직 면은 같은 x 두 점 ta = np.nanmax(self.top_a, axis=1) tb = np.nanmax(self.top_b, axis=1) @@ -295,11 +300,17 @@ def _normal(sl, F, lam, grad=0): p = _pre(sl) g = 1.0 / np.asarray(F, float)[:, None] L = np.asarray(lam, float)[:, None] - P1 = p["cusa"] + L * p["cuca"] - Q1 = p["sat"] + L * p["tca"] - den = sl["ca"] - L * sl["sa"] + g * Q1 + if L.any(): + P1 = p["cusa"] + L * p["cuca"] + Q1 = p["sat"] + L * p["tca"] + den = sl["ca"] - L * sl["sa"] + g * Q1 + top = sl["W"] + L * p["lw"] + else: # λ = 0(Bishop · Janbu) — 같은 식에서 λ 칸만 뺌 + P1, Q1 = p["cusa"], p["sat"] + den = sl["ca"] + g * Q1 + top = sl["W"] safe = np.where(np.abs(den) < _TINY, _TINY, den) - N = (sl["W"] + L * p["lw"] - g * P1) / safe + N = (top - g * P1) / safe if not grad: return N, den dg = -(P1 + N * Q1) / safe @@ -360,29 +371,37 @@ def _ff(sl, F, lam, grad=0, nd=None): return (*_balance(p["tca"], sl["sa"], p["rc0"], p["fd0"], nd, F, False), nd[1]) -def _iterate(fun, F0): - """F = fun(F) 풀이(묶음 전체) — 해석 도함수 뉴턴 · 안 되면 반쯤 섞은 고정점 걸음. +def _iterate(sl, center, lam, F0, moment=True): + """F = f(F) 풀이 — f = 모멘트(F_m) 또는 힘(F_f) 평형 · 해석 도함수 뉴턴 · 안 되면 반쯤 섞은 고정점 걸음. - fun(F) → (f, df/dF, …). 단순 고정점은 급한 머리 절편에서 진동해 느림(도함수 ≈ −1) — 뉴턴은 몇 번에 끝남. + 단순 고정점은 급한 머리 절편에서 진동해 느림(도함수 ≈ −1) — 뉴턴은 몇 번에 끝남. + 끝난 면은 빼고 셈(남은 면이 절반 아래면 묶음을 줄임 · 면마다 따로 셈이라 값 같음). """ F = np.asarray(F0, float).copy() done = np.zeros(len(F), bool) iters = np.zeros(len(F), int) + idx, w, wc, wl = np.arange(len(F)), sl, center, lam for _ in range(MAX_ITER): - f1, df = fun(F)[:2] - h1 = f1 - F + Fi, di = F[idx], done[idx] + f1, df = (_fm(w, wc, Fi, wl, 1) if moment else _ff(w, Fi, wl, 1))[:2] + h1 = f1 - Fi dh = df - 1.0 - newton = F - h1 / np.where(np.abs(dh) < 1e-12, -1.0, dh) - relaxed = 0.5 * (F + f1) - use = np.isfinite(newton) & (dh < -1e-9) & (newton > 0.5 * F) & (newton < 2.0 * F) + newton = Fi - h1 / np.where(np.abs(dh) < 1e-12, -1.0, dh) + relaxed = 0.5 * (Fi + f1) + use = np.isfinite(newton) & (dh < -1e-9) & (newton > 0.5 * Fi) & (newton < 2.0 * Fi) Fn = np.where(use, newton, relaxed) Fn = np.where(np.isfinite(Fn) & (Fn > 0), Fn, np.nan) - conv = np.abs(h1) < TOL * np.maximum(1.0, np.abs(F)) - F = np.where(done | conv, F, Fn) - iters += ~done - done |= conv | ~np.isfinite(F) - if done.all(): + conv = np.abs(h1) < TOL * np.maximum(1.0, np.abs(Fi)) + F[idx] = np.where(di | conv, Fi, Fn) + iters[idx] += ~di + done[idx] = di | conv | ~np.isfinite(F[idx]) + act = ~done[idx] + if not act.any(): break + if act.sum() < 0.5 * len(idx): + idx, w = idx[act], _subset(w, act) + wc = None if wc is None else wc[act] + wl = wl[act] return F, iters, done & np.isfinite(F) @@ -406,7 +425,7 @@ def ordinary(sl: dict, center: np.ndarray) -> dict: def bishop(sl: dict, center: np.ndarray, F0=None) -> dict: F0 = np.full(len(sl["W"]), 1.5) if F0 is None else F0 lam = np.zeros(len(F0)) - F, it, conv = _iterate(lambda f: _fm(sl, center, f, lam, 1), F0) + F, it, conv = _iterate(sl, center, lam, F0) return _finish(sl, F, it, conv, lam) @@ -414,7 +433,7 @@ def janbu(sl: dict, F0=None) -> dict: """Janbu 간편(F_0) + 보정 f_0 — 지식DB 02 §2.3 · b1 = 0.31(c=0) · 0.69(φ=0) · 0.50.""" F0 = np.full(len(sl["W"]), 1.5) if F0 is None else F0 lam = np.zeros(len(F0)) - F, it, conv = _iterate(lambda f: _ff(sl, f, lam, 1), F0) + F, it, conv = _iterate(sl, None, lam, F0, moment=False) out = _finish(sl, F, it, conv, lam) out["F0"] = out["F"] out["f0"] = janbu_f0(sl) @@ -461,7 +480,7 @@ def spencer(sl: dict, center: np.ndarray, F0=None, scan: bool = True) -> dict: return out -def _spencer_newton(sl, center, F, lam, iters=40): +def _spencer_newton(sl, center, F, lam, iters=25): """(F_m − F, F_f − F) = 0 뉴턴 — 걸음 크기 제한(F 절반 · λ 0.3) · 수렴 = 두 잔차 < 1e-4 · λ 범위 안. 끝난 면은 빼고(남은 면이 절반 아래면 묶음을 줄임) · F < F_MIN 으로 빠진 면은 버림(어차피 `_finish` 가 버림). @@ -514,8 +533,8 @@ def lambda_curve(sl, center, grid=None): F0 = np.full(S, 1.5) for k, g in enumerate(grid): lam = np.full(S, g) - fm_all[k] = _iterate(lambda f: _fm(sl, center, f, lam, 1), F0)[0] - ff_all[k] = _iterate(lambda f: _ff(sl, f, lam, 1), F0)[0] + fm_all[k] = _iterate(sl, center, lam, F0)[0] + ff_all[k] = _iterate(sl, None, lam, F0, moment=False)[0] return grid, fm_all, ff_all