课程摘要上一节检查了材料场离散、KL 模态和抽样次数本节补上另一个关键问题有限元方程组是否真正求解到位我们用二维三角形单元建立一个非均匀材料的反平面剪切模型分别采用稀疏直接法、共轭梯度法CG和 Jacobi 预条件 CG 求解。通过位移场、真实残差和解的差异学习判断求解器是否收敛以及为什么“残差很小”仍不能代替网格与物理模型验证。一、有限元组装完成后还没有“算出结果”无论前面的模型是拉杆、热传导还是二维弹性体离散后通常都会得到\[ K\boldsymbol u\boldsymbol f \]其中 \(K\) 是整体刚度矩阵\(\boldsymbol u\) 是未知节点位移\(\boldsymbol f\) 是载荷。真实网格中大多数节点只与邻近节点相连所以 \(K\) 的绝大部分元素为零——它是稀疏矩阵。第四十三节提醒我们检查空间离散和随机抽样本节再增加一层检查线性方程组是否求解充分。如果迭代法过早停止即使材料场和网格设置正确输出仍会包含求解器误差。不过也不能走向另一个极端把线性残差压到极小不会自动消除不合适的边界条件或粗网格误差。二、本节的力学算例考虑单位正方形截面上的反平面剪切位移 \(w(x,y)\)。材料的剪切系数 \(\mu(x,y)\) 在中心附近明显增大方程为\[ -\nabla\cdot\bigl(\mu(x,y)\nabla w\bigr)q \]本节取 \(q1\)四周边界 \(w0\)并用以下正值函数构造非均匀材料\[ \mu(x,y) 199\exp\left[ -\frac{(x-0.5)^2(y-0.5)^2}{0.015} \right] \]这是一个无量纲教学模型因此后面报告的 \(w\) 不标作毫米。乘以测试函数 \(v\) 并积分弱形式为\[ \int_{\Omega} \mu\nabla w\cdot\nabla v\,\mathrm d\Omega \int_{\Omega}qv\,\mathrm d\Omega \]采用一阶三角形单元时每个单元的梯度为常数。代码在三角形中心计算一次 \(\mu\)形成单元矩阵并组装稀疏整体矩阵。施加零位移边界条件后待求矩阵是对称正定的因此可以使用 CGCG 对矩阵及预条件器的对称正定要求也见 PETSc 官方说明。三、三种求解方式分别做什么方法本节用途要观察的量稀疏直接法为同一离散方程组提供数值对照解与残差CG不显式做完整分解逐步逼近解迭代次数、残差Jacobi 预条件 CG用刚度矩阵对角项改善迭代系统是否减少迭代本节预条件器做的事很简单\[ M^{-1}\boldsymbol r \operatorname{diag}(K)^{-1}\boldsymbol r \]它不是修改物理模型也不改变最终要解的 \(K\boldsymbol u\boldsymbol f\)它只是帮助迭代法更有效地找到解。预条件效果依赖具体矩阵不能保证每个问题都出现相同收益。SciPy CG 文档四、完整可运行代码第四十四节二维有限元稀疏方程组的直接法与预条件 CG。 from pathlib import Path import matplotlib matplotlib.use(Agg) import matplotlib.pyplot as plt import matplotlib.tri as mtri import numpy as np from scipy.sparse import coo_matrix from scipy.sparse.linalg import LinearOperator, cg, spsolve # 单位正方形上的反平面剪切 # -div(mu grad w)1边界 w0。 N 40 grid np.linspace(0.0, 1.0, N 1) xx, yy np.meshgrid(grid, grid) nodes np.column_stack((xx.ravel(), yy.ravel())) triangles [] for j in range(N): for i in range(N): a j * (N 1) i b, c, d a 1, a N 1, a N 2 triangles.extend(((a, b, d), (a, d, c))) triangles np.asarray(triangles, dtypeint) rows, cols, data [], [], [] load np.zeros(len(nodes)) for tri in triangles: xy nodes[tri] x, y xy[:, 0], xy[:, 1] area np.linalg.det( np.column_stack((np.ones(3), x, y)) ) / 2 assert area 0 grad np.column_stack(( np.roll(y, -1) - np.roll(y, 1), np.roll(x, 1) - np.roll(x, -1) )) / (2 * area) center xy.mean(axis0) mu 1 99 * np.exp( -np.sum((center - 0.5)**2) / 0.015 ) ke mu * area * grad grad.T rows.extend(np.repeat(tri, 3)) cols.extend(np.tile(tri, 3)) data.extend(ke.ravel()) load[tri] area / 3 global_k coo_matrix( (data, (rows, cols)), shape(len(nodes), len(nodes)) ).tocsr() boundary ( np.isclose(nodes[:, 0], 0) | np.isclose(nodes[:, 0], 1) | np.isclose(nodes[:, 1], 0) | np.isclose(nodes[:, 1], 1) ) free np.flatnonzero(~boundary) k global_k[free][:, free].tocsr() f load[free] f_norm np.linalg.norm(f) # 对同一个离散系统求解先取得直接法数值对照。 direct spsolve(k, f) plain_history, jacobi_history [], [] def record(history): return lambda x: history.append( np.linalg.norm(f - k x) / f_norm ) plain, plain_info cg( k, f, rtol1e-8, atol0.0, maxiter5000, callbackrecord(plain_history) ) inverse_diagonal 1 / k.diagonal() jacobi LinearOperator( k.shape, matveclambda x: inverse_diagonal * x ) preconditioned, jacobi_info cg( k, f, Mjacobi, rtol1e-8, atol0.0, maxiter5000, callbackrecord(jacobi_history) ) assert plain_info jacobi_info 0 assert np.all(k.diagonal() 0) assert ( np.max(np.abs((k - k.T).data), initial0.0) 1e-10 ) np.testing.assert_allclose( plain, direct, rtol1e-5, atol1e-8 ) np.testing.assert_allclose( preconditioned, direct, rtol1e-5, atol1e-8 ) def residual(solution): return ( np.linalg.norm(f - k solution) / f_norm ) full_solution np.zeros(len(nodes)) full_solution[free] direct print( f节点/三角形/自由度: f{len(nodes)}/{len(triangles)}/{len(free)} ) print(f稀疏矩阵非零元数: {k.nnz}) print( f最大反平面位移: f{full_solution.max():.6f} ) print( f直接法相对残差: f{residual(direct):.3e} ) print( fCG: {len(plain_history)} 次迭代 f相对残差 {residual(plain):.3e} ) print( fJacobi-PCG: {len(jacobi_history)} 次迭代 f相对残差 {residual(preconditioned):.3e} ) print( fCG 与直接法最大绝对差: f{np.max(np.abs(plain - direct)):.3e} ) print( fJacobi-PCG 与直接法最大绝对差: f{np.max(np.abs(preconditioned - direct)):.3e} ) print( 检查通过矩阵对称、迭代收敛 且三种解一致。 ) out Path(__file__).parent fig, ax plt.subplots(figsize(7, 6), dpi160) mesh mtri.Triangulation( nodes[:, 0], nodes[:, 1], triangles ) field ax.tripcolor( mesh, full_solution, shadinggouraud, cmapviridis ) fig.colorbar( field, axax, labelOut-of-plane displacement w ) ax.set( xlabelx, ylabely, titleFinite-element displacement field ) ax.set_aspect(equal) fig.tight_layout() fig.savefig(out / lesson44_displacement_field.png) plt.close(fig) fig, ax plt.subplots( figsize(8.5, 4.8), dpi160 ) ax.semilogy( range(1, len(plain_history) 1), plain_history, labelCG, linewidth2 ) ax.semilogy( range(1, len(jacobi_history) 1), jacobi_history, labelJacobi-preconditioned CG, linewidth2 ) ax.axhline( 1e-8, color#c2410c, linestyle--, labelRelative tolerance ) ax.set( xlabelIteration, ylabelTrue relative residual, titleIterative solver convergence ) ax.grid(alpha0.2, whichboth) ax.legend(frameonFalse) fig.tight_layout() fig.savefig(out / lesson44_solver_residuals.png) plt.close(fig) print( 图片已保存lesson44_displacement_field.png、 lesson44_solver_residuals.png )本课代码使用cg(..., rtol..., atol...)的较新接口若运行环境较旧请先核对所装 SciPy 的函数签名。稀疏直接求解与 CG 接口见 SciPy 官方文档。实际运行输出节点/三角形/自由度: 1681/3200/1521 稀疏矩阵非零元数: 10337 最大反平面位移: 0.056282 直接法相对残差: 9.504e-13 CG: 478 次迭代相对残差 7.722e-09 Jacobi-PCG: 94 次迭代相对残差 7.175e-09 CG 与直接法最大绝对差: 1.203e-11 Jacobi-PCG 与直接法最大绝对差: 7.255e-12 检查通过矩阵对称、迭代收敛且三种解一致。 图片已保存lesson44_displacement_field.png、lesson44_solver_residuals.png五、看图结果长什么样求解过程又怎样位移场在四周边界为零内部出现最大位移。中心附近较高的 \(\mu\) 会改变位移场的空间分布颜色图展示的是位移不是材料系数或应力。下图记录每次迭代后重新计算的“真实相对残差”\[ \rho_k \frac{\|\boldsymbol f-K\boldsymbol u_k\|_2} {\|\boldsymbol f\|_2} \]这里 CG 曲线局部上升、起伏并不自动表示程序出错应检查求解器返回状态和最终残差。图中 Jacobi 预条件使本例从478次迭代降至94次但这仅是迭代次数比较不能未经计时就称运行速度提高了相同倍数。SciPy 关于停止准则与预条件器的说明六、怎样正确理解“残差足够小”本例直接法、CG 和预条件 CG 的解高度一致说明它们对同一个离散方程组给出了相符结果。但要严格区分三件事**方程残差小**当前数值解较好地满足已组装的 \(K\boldsymbol u\boldsymbol f\)。**离散误差小**当前网格和单元近似能够充分表示连续问题这需要另做网格收敛检查。**物理模型可信**材料参数、载荷与边界条件确实对应目标工程问题这需要物理证据和验证。小残差只直接支持第一项。尤其当矩阵条件较差时残差与解误差也不能简单画等号。因此本课除了检查残差还比较三种求解结果在真实项目中则应结合关注的位移、应力或反力检查。七、放回 FEniCSx 工作流在 FEniCSx 中我们通常先用 UFL 写出弱形式后由 PETSc 处理组装后的线性系统。课程第三十六至三十九节使用的直接求解设置属于这一层本课帮助理解其中“直接法”“迭代法”和“预条件”的含义。当前 DOLFINx 的LinearProblem文档列出了通过 PETSc 选项配置线性求解器的方式。DOLFINx LinearProblem 文档选择求解器前先确认矩阵性质。像本课这样施加了适当位移约束的对称正定问题适合 CG不能见到“有限元矩阵”四个字就对任意非对称或不定问题直接套用 CG。八、课后练习将 CG 的rtol从 \(10^{-8}\) 改为 \(10^{-4}\)。观察迭代次数、最终真实残差和“与直接法的最大绝对差”如何变化。将网格参数N从40改为20和60。记录自由度、非零元数与两种 CG 迭代次数不要只凭一次运行时间判断哪个算法更优。暂时把中心高刚度系数99改为0使材料均匀。比较两条残差曲线并思考预条件效果为什么随问题改变。这一节完成了从“有限元离散”到“可靠求解线性系统”的连接。下一步进入 Abaqus 脚本化时同样要关注作业是否真正完成、求解是否收敛以及结果是否通过物理检查。
【零基础学智能仿真-44】有限元线性方程组——直接法、共轭梯度法与预条件
课程摘要上一节检查了材料场离散、KL 模态和抽样次数本节补上另一个关键问题有限元方程组是否真正求解到位我们用二维三角形单元建立一个非均匀材料的反平面剪切模型分别采用稀疏直接法、共轭梯度法CG和 Jacobi 预条件 CG 求解。通过位移场、真实残差和解的差异学习判断求解器是否收敛以及为什么“残差很小”仍不能代替网格与物理模型验证。一、有限元组装完成后还没有“算出结果”无论前面的模型是拉杆、热传导还是二维弹性体离散后通常都会得到\[ K\boldsymbol u\boldsymbol f \]其中 \(K\) 是整体刚度矩阵\(\boldsymbol u\) 是未知节点位移\(\boldsymbol f\) 是载荷。真实网格中大多数节点只与邻近节点相连所以 \(K\) 的绝大部分元素为零——它是稀疏矩阵。第四十三节提醒我们检查空间离散和随机抽样本节再增加一层检查线性方程组是否求解充分。如果迭代法过早停止即使材料场和网格设置正确输出仍会包含求解器误差。不过也不能走向另一个极端把线性残差压到极小不会自动消除不合适的边界条件或粗网格误差。二、本节的力学算例考虑单位正方形截面上的反平面剪切位移 \(w(x,y)\)。材料的剪切系数 \(\mu(x,y)\) 在中心附近明显增大方程为\[ -\nabla\cdot\bigl(\mu(x,y)\nabla w\bigr)q \]本节取 \(q1\)四周边界 \(w0\)并用以下正值函数构造非均匀材料\[ \mu(x,y) 199\exp\left[ -\frac{(x-0.5)^2(y-0.5)^2}{0.015} \right] \]这是一个无量纲教学模型因此后面报告的 \(w\) 不标作毫米。乘以测试函数 \(v\) 并积分弱形式为\[ \int_{\Omega} \mu\nabla w\cdot\nabla v\,\mathrm d\Omega \int_{\Omega}qv\,\mathrm d\Omega \]采用一阶三角形单元时每个单元的梯度为常数。代码在三角形中心计算一次 \(\mu\)形成单元矩阵并组装稀疏整体矩阵。施加零位移边界条件后待求矩阵是对称正定的因此可以使用 CGCG 对矩阵及预条件器的对称正定要求也见 PETSc 官方说明。三、三种求解方式分别做什么方法本节用途要观察的量稀疏直接法为同一离散方程组提供数值对照解与残差CG不显式做完整分解逐步逼近解迭代次数、残差Jacobi 预条件 CG用刚度矩阵对角项改善迭代系统是否减少迭代本节预条件器做的事很简单\[ M^{-1}\boldsymbol r \operatorname{diag}(K)^{-1}\boldsymbol r \]它不是修改物理模型也不改变最终要解的 \(K\boldsymbol u\boldsymbol f\)它只是帮助迭代法更有效地找到解。预条件效果依赖具体矩阵不能保证每个问题都出现相同收益。SciPy CG 文档四、完整可运行代码第四十四节二维有限元稀疏方程组的直接法与预条件 CG。 from pathlib import Path import matplotlib matplotlib.use(Agg) import matplotlib.pyplot as plt import matplotlib.tri as mtri import numpy as np from scipy.sparse import coo_matrix from scipy.sparse.linalg import LinearOperator, cg, spsolve # 单位正方形上的反平面剪切 # -div(mu grad w)1边界 w0。 N 40 grid np.linspace(0.0, 1.0, N 1) xx, yy np.meshgrid(grid, grid) nodes np.column_stack((xx.ravel(), yy.ravel())) triangles [] for j in range(N): for i in range(N): a j * (N 1) i b, c, d a 1, a N 1, a N 2 triangles.extend(((a, b, d), (a, d, c))) triangles np.asarray(triangles, dtypeint) rows, cols, data [], [], [] load np.zeros(len(nodes)) for tri in triangles: xy nodes[tri] x, y xy[:, 0], xy[:, 1] area np.linalg.det( np.column_stack((np.ones(3), x, y)) ) / 2 assert area 0 grad np.column_stack(( np.roll(y, -1) - np.roll(y, 1), np.roll(x, 1) - np.roll(x, -1) )) / (2 * area) center xy.mean(axis0) mu 1 99 * np.exp( -np.sum((center - 0.5)**2) / 0.015 ) ke mu * area * grad grad.T rows.extend(np.repeat(tri, 3)) cols.extend(np.tile(tri, 3)) data.extend(ke.ravel()) load[tri] area / 3 global_k coo_matrix( (data, (rows, cols)), shape(len(nodes), len(nodes)) ).tocsr() boundary ( np.isclose(nodes[:, 0], 0) | np.isclose(nodes[:, 0], 1) | np.isclose(nodes[:, 1], 0) | np.isclose(nodes[:, 1], 1) ) free np.flatnonzero(~boundary) k global_k[free][:, free].tocsr() f load[free] f_norm np.linalg.norm(f) # 对同一个离散系统求解先取得直接法数值对照。 direct spsolve(k, f) plain_history, jacobi_history [], [] def record(history): return lambda x: history.append( np.linalg.norm(f - k x) / f_norm ) plain, plain_info cg( k, f, rtol1e-8, atol0.0, maxiter5000, callbackrecord(plain_history) ) inverse_diagonal 1 / k.diagonal() jacobi LinearOperator( k.shape, matveclambda x: inverse_diagonal * x ) preconditioned, jacobi_info cg( k, f, Mjacobi, rtol1e-8, atol0.0, maxiter5000, callbackrecord(jacobi_history) ) assert plain_info jacobi_info 0 assert np.all(k.diagonal() 0) assert ( np.max(np.abs((k - k.T).data), initial0.0) 1e-10 ) np.testing.assert_allclose( plain, direct, rtol1e-5, atol1e-8 ) np.testing.assert_allclose( preconditioned, direct, rtol1e-5, atol1e-8 ) def residual(solution): return ( np.linalg.norm(f - k solution) / f_norm ) full_solution np.zeros(len(nodes)) full_solution[free] direct print( f节点/三角形/自由度: f{len(nodes)}/{len(triangles)}/{len(free)} ) print(f稀疏矩阵非零元数: {k.nnz}) print( f最大反平面位移: f{full_solution.max():.6f} ) print( f直接法相对残差: f{residual(direct):.3e} ) print( fCG: {len(plain_history)} 次迭代 f相对残差 {residual(plain):.3e} ) print( fJacobi-PCG: {len(jacobi_history)} 次迭代 f相对残差 {residual(preconditioned):.3e} ) print( fCG 与直接法最大绝对差: f{np.max(np.abs(plain - direct)):.3e} ) print( fJacobi-PCG 与直接法最大绝对差: f{np.max(np.abs(preconditioned - direct)):.3e} ) print( 检查通过矩阵对称、迭代收敛 且三种解一致。 ) out Path(__file__).parent fig, ax plt.subplots(figsize(7, 6), dpi160) mesh mtri.Triangulation( nodes[:, 0], nodes[:, 1], triangles ) field ax.tripcolor( mesh, full_solution, shadinggouraud, cmapviridis ) fig.colorbar( field, axax, labelOut-of-plane displacement w ) ax.set( xlabelx, ylabely, titleFinite-element displacement field ) ax.set_aspect(equal) fig.tight_layout() fig.savefig(out / lesson44_displacement_field.png) plt.close(fig) fig, ax plt.subplots( figsize(8.5, 4.8), dpi160 ) ax.semilogy( range(1, len(plain_history) 1), plain_history, labelCG, linewidth2 ) ax.semilogy( range(1, len(jacobi_history) 1), jacobi_history, labelJacobi-preconditioned CG, linewidth2 ) ax.axhline( 1e-8, color#c2410c, linestyle--, labelRelative tolerance ) ax.set( xlabelIteration, ylabelTrue relative residual, titleIterative solver convergence ) ax.grid(alpha0.2, whichboth) ax.legend(frameonFalse) fig.tight_layout() fig.savefig(out / lesson44_solver_residuals.png) plt.close(fig) print( 图片已保存lesson44_displacement_field.png、 lesson44_solver_residuals.png )本课代码使用cg(..., rtol..., atol...)的较新接口若运行环境较旧请先核对所装 SciPy 的函数签名。稀疏直接求解与 CG 接口见 SciPy 官方文档。实际运行输出节点/三角形/自由度: 1681/3200/1521 稀疏矩阵非零元数: 10337 最大反平面位移: 0.056282 直接法相对残差: 9.504e-13 CG: 478 次迭代相对残差 7.722e-09 Jacobi-PCG: 94 次迭代相对残差 7.175e-09 CG 与直接法最大绝对差: 1.203e-11 Jacobi-PCG 与直接法最大绝对差: 7.255e-12 检查通过矩阵对称、迭代收敛且三种解一致。 图片已保存lesson44_displacement_field.png、lesson44_solver_residuals.png五、看图结果长什么样求解过程又怎样位移场在四周边界为零内部出现最大位移。中心附近较高的 \(\mu\) 会改变位移场的空间分布颜色图展示的是位移不是材料系数或应力。下图记录每次迭代后重新计算的“真实相对残差”\[ \rho_k \frac{\|\boldsymbol f-K\boldsymbol u_k\|_2} {\|\boldsymbol f\|_2} \]这里 CG 曲线局部上升、起伏并不自动表示程序出错应检查求解器返回状态和最终残差。图中 Jacobi 预条件使本例从478次迭代降至94次但这仅是迭代次数比较不能未经计时就称运行速度提高了相同倍数。SciPy 关于停止准则与预条件器的说明六、怎样正确理解“残差足够小”本例直接法、CG 和预条件 CG 的解高度一致说明它们对同一个离散方程组给出了相符结果。但要严格区分三件事**方程残差小**当前数值解较好地满足已组装的 \(K\boldsymbol u\boldsymbol f\)。**离散误差小**当前网格和单元近似能够充分表示连续问题这需要另做网格收敛检查。**物理模型可信**材料参数、载荷与边界条件确实对应目标工程问题这需要物理证据和验证。小残差只直接支持第一项。尤其当矩阵条件较差时残差与解误差也不能简单画等号。因此本课除了检查残差还比较三种求解结果在真实项目中则应结合关注的位移、应力或反力检查。七、放回 FEniCSx 工作流在 FEniCSx 中我们通常先用 UFL 写出弱形式后由 PETSc 处理组装后的线性系统。课程第三十六至三十九节使用的直接求解设置属于这一层本课帮助理解其中“直接法”“迭代法”和“预条件”的含义。当前 DOLFINx 的LinearProblem文档列出了通过 PETSc 选项配置线性求解器的方式。DOLFINx LinearProblem 文档选择求解器前先确认矩阵性质。像本课这样施加了适当位移约束的对称正定问题适合 CG不能见到“有限元矩阵”四个字就对任意非对称或不定问题直接套用 CG。八、课后练习将 CG 的rtol从 \(10^{-8}\) 改为 \(10^{-4}\)。观察迭代次数、最终真实残差和“与直接法的最大绝对差”如何变化。将网格参数N从40改为20和60。记录自由度、非零元数与两种 CG 迭代次数不要只凭一次运行时间判断哪个算法更优。暂时把中心高刚度系数99改为0使材料均匀。比较两条残差曲线并思考预条件效果为什么随问题改变。这一节完成了从“有限元离散”到“可靠求解线性系统”的连接。下一步进入 Abaqus 脚本化时同样要关注作业是否真正完成、求解是否收敛以及结果是否通过物理检查。