solver モジュール

線形系の求解・Newton 法・疑似過渡継続を担当する。

モジュール構成

solver/
├── mod.rs               — サブモジュールの公開宣言
├── banded.rs            — 帯行列 LU 分解・求解
├── dense.rs             — 密行列 LU 分解・求解 (faer)
├── newton.rs            — 減衰 Newton 法（直線探索付き）
└── pseudo_transient.rs  — 疑似過渡継続 (PT)

banded.rs — 帯行列

BandedMatrix

LAPACK 互換の一般帯行列 (GB) 形式で係数行列を保持する。

格納レイアウト: data[col × ldab + (kl + ku + row − col)]

• ldab = 2×kl + ku + 1（部分ピボッティングの fill-in 分の kl 行を余分に確保）

fill-in の扱い: ピボット行入れ替え後、行 k の非ゼロ範囲は k+ku+kl 列まで 拡張しうる。後退代入でも同じ幅を使用することで fill-in を正確に処理する。

テスト (3 テスト、全 PASS)

dense.rs — 密行列

DenseMatrix

faer の部分ピボッティング LU で n×n 密行列を求解する。

帯行列では帯外に落ちる質量流束 M 列（解ベクトル末尾、全残差に結合）を 正確に扱うために使用する。

テスト (2 テスト、全 PASS)

newton.rs — 減衰 Newton 法

アルゴリズム概要

for iter in 0..max_iter:
    評価: F(x)
    収束判定: ‖F‖ < atol·√n + rtol·‖x‖
    Jacobian: (J_int, m_col, jfc) = numerical_jacobian(x, F)  ← 最大 max_jac_age 反復再利用
    線形系:   枠付き帯行列系を bordered_solve_m で解く
                δx_int = J_int⁻¹(-F_int + F_M / (jfc - m_col·J_int⁻¹·e_M) × m_col)
                δM    = シュア補完で決定
    ステップクリップ: 物理スケール上限を超える成分を縮小
    直線探索: α を 1 から 0.5 刻みで最大 12 回試行
    更新:     x ← x + α·step
    物理クランプ: T, Y, M を有効範囲に制限
    失速検出: best ‖F‖ が 1% 以上改善しない反復が max_stall_iter 回続いたら受理

帯行列 + 枠付きシュア補完方式: numerical_jacobian が返す (n-1)×(n-1) 帯行列 J_int、 M 列ベクトル m_col、M 行ベクトル jfc を使って bordered_solve_m が 1 回の帯 LU で δM を含む完全な Newton ステップを求める。密 Jacobian 不要で O(n·kl) のコスト。

NewtonConfig

Jacobian の強制再計算: ‖F‖ が前回評価時から 10 倍以上改善した場合、 線形化点が大きく変わったとみなして jac_age をリセットする。

ステップクリップ (step_clip_factor)

物理クランプ (clamp_physical)

テスト (2 テスト、全 PASS)

pseudo_transient.rs — 疑似過渡継続

目的

Newton 法の収束域外（点火前の急勾配・スティッフ化学反応領域）から Newton 域へ 誘導する。TWOPNT スタイルの後退オイラー法を採用。

アルゴリズム概要

for step in 0..n_steps:
    評価: F_steady(x)
    ‖F_steady‖ < 1e7 → Newton 域に達したので終了
    dt < 1e-13 → 収束断念して Newton へ

    rdt = 1/dt
    x_old = x

    内側 Newton (max_inner_iter 回):
        F_pt(x) = F_steady(x) + D/dt·(x − x_old)
        J_pt    = J_steady(x) + D/dt の対角（ρ·cp で T、ρ で Yk、1 で M）
        M 列をゼロクリア（M の疑似時間項は対角のみ）
        step = −J_pt⁻¹·F_pt
        ステップクリップ + 直線探索 + 物理クランプ

    外側受理: ‖F_steady(x_new)‖ ≤ 1.1×‖F_steady(x_old)‖
        成功 → dt *= dt_grow (最大 dt_max まで)
        失敗 → x = x_old, dt *= 0.5; 30 回連続失敗で中断

PseudoTransientConfig

設計上の注意

直線探索の初期化: 内側直線探索の best_norm を ∞ で初期化することで、 J_pt が不定値（D/dt < |λ_chem|）のとき全 α で F_pt が増加しても α = 0 に停滞しない。最小の有限 F_pt を与える α を返す。

外側受理の 10% マージン: 化学反応が非常に剛な点火フェーズでは、 良い内側 Newton ステップでも F_steady が一時的に増加することがある。 1.1 倍の許容マージンにより探索的なステップを通過させ、Newton 域への 到達を可能にする。

D/dt 対角: ρ·cp (T 方程式)、ρ (Yk 方程式)、1 (M) で各変数の 次元が整合した疑似時間スケーリングを実現する（residual.rs の PT 埋め込みと同一）。

テスト (2 テスト、全 PASS)