Possible upper/lower triangle storage inconsistency in rates() in multi-phase equilibrium solver

未关闭
  • 哪里不对？
#331 0 条评论 0 个 reaction 已指派 0 人 在 GitHub 查看

还没有人认领这个 Issue。

评估

难度
3/5
预计耗时
半天
新手友好度
50/100
Issue 类型
缺陷
描述清晰度
基本清楚
活跃度
活跃
技术栈
cpp
领域
backend

调研方向

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 小时
30 天内合并 PR
2

环境准备

  • 没有 Dockerfile 或 Docker Compose 文件
  • 没有 Pull Request 模板
  • 阅读贡献指南

从这里开始

  1. 先读完整个 Issue，再读项目的贡献指南。
  2. 在 Issue 下留言说明你要接手 —— 这能避免两个人做同样的事。
  3. Fork 仓库，在一个分支上完成修改。
  4. 提交 Pull Request，并在描述里引用这个 Issue 编号。

mutationpp/Mutationpp 的其他 Issue

查看 mutationpp/Mutationpp 的全部 Issue

相似的 Issue

更多 C++ Issue

把新 issue 发到你的邮箱

精选适合新手参与的 GitHub issue 摘要。