Possible upper/lower triangle storage inconsistency in rates() in multi-phase equilibrium solver
まだ誰も着手していません。
評価
調査の方向性
The suspect code is in rates() in src/thermo/MultiPhaseEquilSolver.cpp: the loop assembling A(m,p) only fills the lower triangle, and the solve near line 1136 uses selfadjointView<Upper>(). Read the assembly loop and that solve to confirm whether the upper triangle is filled elsewhere before the solve. Done means either switching to selfadjointView<Lower>() or populating the upper triangle, with a maintainer confirming which one is intended and the solver results checked against a known equilibrium case.
索引モデルが issue の本文から書いたものです。
説明
I believe there may be an inconsistency in the rates() routine between the way the matrix A is assembled and the triangular part subsequently used by Eigen's self-adjoint solver.
A is populated with:
for (int p = 0; p < npr; ++p) {
for (int m = p; m < npr; ++m) {
for (int j = p_sizes[m]; j < p_sizes[m+1]; ++j) {
sum = 0.0;
for (int i = 0; i < ncr; ++i)
sum += H(j,i)*dlamy(i,p);
A(m,p) += y[j]*sum;
}
}
}
Since m >= p, this explicitly computes the lower triangular part of A, i.e. A(m,p) for m >= p. This appears intentional, since A is expected to be symmetric.
However, the resulting system is then solved with
https://github.com/mutationpp/Mutationpp/blob/7607ef4d67f4de115705f78f467b1ea48bf02c26/src/thermo/MultiPhaseEquilSolver.cpp#L1136
As I understand Eigen's SelfAdjointView, selfadjointView<Upper>() treats the upper triangular portion of the matrix as the stored representation of the self-adjoint matrix and does not use the corresponding entries in the lower triangle. I believe this is therefore causing the off-diagonal elements to be discarded in the calculation of dx.
In other words, since the assembly produces
[ a 0 0
b c 0
d e f ]
the intended symmetric matrix is presumably
[ a b d
b c e
d e f ]
However, since the upper-triangular off-diagonal elements remain zero, selfadjointView<Upper>() represents the matrix as
[ a 0 0
0 c 0
0 0 f ]
instead.
Should selfadjointView<Lower>() be used here instead? Alternatively, is there some reason the upper triangle is expected to be populated elsewhere before this solve?
- 主要言語
- C++
- スター
- 145
- フォーク
- 79
- 平均マージ
- 1日 15時間
- マージ済み PR（30日）
- 2
環境構築
- Dockerfile・Docker Compose ファイルなし
- プルリクエストのテンプレートなし
- コントリビューションガイドを読む
はじめの一歩
- issue を最後まで読み、次にプロジェクトのコントリビューションガイドを読みます。
- 着手することを issue にコメントします — 二人が同じ作業をするのを防げます。
- リポジトリをフォークし、ブランチを切って変更します。
- issue 番号を参照したプルリクエストを送ります。
mutationpp/Mutationpp のほかの issue
-
難易度 2/5 1〜3時間 初心者へのやさしさ 68/100
mutationpp/Mutationpp#337 ·
-
難易度 1/5 1時間未満 初心者へのやさしさ 74/100
mutationpp/Mutationpp#302 ·
-
難易度 4/5 3〜5日 初心者へのやさしさ 51/100
mutationpp/Mutationpp#335 ·
-
難易度 3/5 1〜2日 初心者へのやさしさ 42/100
mutationpp/Mutationpp#332 ·
-
Prevent Newton solver hangs in the multi-phase equilibrium solver対応中かも @capriatim が 3 日前に担当しました。 オープン
難易度 4/5 3〜5日 初心者へのやさしさ 40/100
mutationpp/Mutationpp#330 · 担当者 1 名 ·
mutationpp/Mutationpp の issue をすべて見る
似ている issue
-
難易度 2/5 1〜3時間 初心者へのやさしさ 72/100
KhronosGroup/Vulkan-Tutorial#534 ·
-
難易度 2/5 1〜3時間 初心者へのやさしさ 65/100
メンテナーはふだん 2 日以内に返信
-
bug IO
難易度 2/5 1〜3時間 初心者へのやさしさ 65/100
MrNeRF/LichtFeld-Studio#3824 ·
メンテナーはふだん 1 日以内に返信
-
pnnx: manual tanh GELU fusion loses the approximation mode対応中かも @ktz03 が今日担当しました。 オープン
難易度 2/5 1〜3時間 初心者へのやさしさ 78/100
Tencent/ncnn#7175 · リアクション 1 件 ·
メンテナーはふだん 1 日以内に返信
-
難易度 2/5 1〜3時間 初心者へのやさしさ 77/100
throneproj/Throne#2043 ·
メンテナーはふだん 2 日以内に返信