通过符号微分加速模型收敛

作者 Walter Frei

2020年 10月 29日

你是否知道 COMSOL Multiphysics® 软件内置了一个符号微分求解器,在设置和解决非线性问题时会自动启用?此功能可确保软件在求解此类模型时具有很高的可靠性,但也要付出一些代价。今天,我们将通过几个不同的示例演示此功能,并研究其优点和弊端。

利用符号微分求解简单的非线性问题

让我们从一个相对简单的、单自由度的非线性问题(隐式方程)开始:

T=Q/(1+T^2/100)

其中, T 表示系统的集总温度,Q 代表热载荷,上式当中的分母 (1+T^2/100) 代表系统热导率。

要求解参数 Q 的不同值,我们可以通过阻尼牛顿法找到非线性问题的解,该方法可以通过在 COMSOL Multiphysics 中写残差来完成:

r(T)=T-Q/(1+T^2/100)

然后将其输入全局方程 接口,如下截图所示。

COMSOL Multiphysics中的全局方程式界面的屏幕截图,其中包含用于求解非线性方程式的残差。
用于求解非线性方程的全局方程接口。

如下图所示,我们可以针对一系列不同的 Q 值来求解此方程。

绘制具有单自由度的非线性问题的解的图。
单自由度非线性问题的解。

为了便于说明,我们在下图中画出了 T 取值范围内和选定值 Q 的残差。零交叉处的点是上图的解。

绘制不同Q值的残差的图形,它表示非线性问题中的热负荷。
不同参数 Q 值的残差图。与零交叉的 T 值就是非线性方程的解。

我们还画出了相同范围 T 的和 Q 值下的残差的对应导数(雅可比矩阵)。

绘制不同热负荷Q值时残差导数的图。
不同 Q 值的残差的导数图。

我们看到,随着 Q 值的增加,雅可比矩函数非线性增强。我们也可以(对于这种简单情况)手动写出其导数:

\frac{\partial r(T)}
{\partial T}=1+\frac{2TQ}{100(1+T^2/100)^2}

请记住,软件能够自动通过符号微分获得此导数,并在牛顿步骤中使用它。

可以看出,正是残差第二部分的分母使该表达式有些复杂(回忆下,该项可以被认为是系统的热阻)。事实证明,如果我们愿意,可以忽略该项,通过 nojac() 算子包裹分母的该项:

R(T) =(T – Q/nojac(1+T^2/100))

并修改全局方程,如下截图所示。

COMSOL Multiphysics 中的“全局方程式设置”窗口演示了如何使用nojac()运算符。
屏幕截图显示了如何使用 nojac() 算子。

结果,符号微分求解器将忽略 nojac() 算子中的所有内容。也就是说,该算子中的项对雅可比矩阵没有贡献。关于 T 的残差导数现在简化为:

\frac{\partial r(T)}{\partial T}=1

通过检查我们可以看到,这种简化的导数仅是 Q 为低值时的一个很好的近似。但是求解器不能对所有的 Q 值,都收敛到解。因此,我们通过避免获取某些项的导数(在这种情况下非常小)使问题变得更简单了,但是在某些情况下会使问题无解。

在 COMSOL Multiphysics® 中使用 nojac() 算子

我们考虑一下可以使用 nojac() 算子的其他几种方式。

如果我们将全局方程写为:

R(T) = nojac(T – Q/(1+T^2/100))

那么导数将恰好为零,并且这个问题根本无法求解,因此这种情况并无不有趣。

我们还可以尝试将 nojac() 仅应用于一个项:

R(T) = T – Q/(1+nojac(T)*T/100)

现在,我们的导数是:

\frac{\partial r(T)}

{\partial T}
=1+\frac{TQ}

{100(1+T^2/100)^2}

该导数非常接近于精确的残差,并且至少在这种情况下,即使采用雅可比矩阵公式的这种近似值,求解器将收敛于所有的 Q 值。尽管计算成本没有下降,但这里值得注意的一点是,我们可以将 nojac() 算子应用于残差项的任何子集。

那么,到目前为止我们学到了什么?

  1. 您可以使用 nojac() 算子来忽略对雅可比矩阵产生贡献的项
  2. 忽略项将影响非线性求解器的收敛
  3. 忽略项可减少计算模型的大小

但是,在这个小例子中,最后一点并不是很明显。

减少大型非线性模型的大小

现在让我们看一个更大的三维模型。我们求解了一个边长为 1cm 的正方体上的传热问题,并在此部分上应用了一个极细 网格,从而我们获得了一个较大的模型;自由度大约为 300 万。

带有参数标记的传热问题示意图。
一个用于研究极细网格的简单传热问题的示意图。

我们还将材料特性定义为温度的函数,并使用类似于我们已经看到的材料特性,通过用户自定义的表达式来定义材料的导热系数:

(1+(T[1/K])^2/100)[W/m/K]

可以直接在物理场接口中完成此操作,如下截图所示。

固体传热界面的“设置”窗口具有手动定义的导热系数设置。
手动将导热系数定义为温度的函数。

当我们使用所有默认设置解决此问题时,该模型将需要用阻尼牛顿法进行四次迭代来求解。整个过程需要大约 135s 的计算时间,并且在一台典型的台式计算机上使用超过 12GB 的内存。在求解过程中,我们还会在日志中收到一条消息,内容为:

找到非对称矩阵。

该消息是材料属性中与温度相关项的结果。也就是说,材料属性中的非线性项使雅可比矩阵不对称,因此该模型将需要更多的内存来求解。现在让我们将整个材料属性表达式包裹在 nojac() 算子中:

nojac(1+(T[1/K])^2/100)[W/m/K]

重新运行该模型,我们看到该模型仍会收敛,需要大约 10GB 内存求解,但现在需要五次默认阻尼牛顿法迭代。也就是说,内存需求下降了,但是收敛的迭代次数却增加了。最有趣的是,执行求解时间为 125s,比以前少了一点!为什么会这样?求解器需要更多的迭代来收敛,但是每次迭代花费的时间更少,因为雅可比矩阵是对称的,需要处理的数据更少。

因此,在这种情况下,尽管需要更多的迭代才能收敛到解,但使用 nojac() 算子既可以减少运行求解时间,又可以减少所需的内存!还要注意,因为使用更新的材料特性表达式,仍然可以正确评估材料特性本身,所以我们能够收敛于相同的解,

避免非局部耦合项

再看最后一个示例:如下图所示,压碎包含可压缩气体的薄壁容器。可以使用结构力学模块中可用的公式进行建模。

包含可压缩气体并承受外部载荷的薄壁储罐模型,演示如何使用符号微分来加速模型收敛
装有可压缩气体并施加有外部压碎载荷的薄壁容器。

通过利用对称性,我们可以将模型缩减为完整模型的1/8。我们施加近似于压碎的分布载荷,以及代表内部可压缩气体的压力载荷。为了表示可压缩气体,类似于在充气超弹性密封条压缩示例中所做的,我们使用了计算变形容器体积的方法。也就是说,通过使用积分耦合算子,我们在容器的整个表面上施加均匀的压力载荷,并且该压力是基于整个容器的变形来计算的。

这种非局部耦合的作用是,容器每个部分的变形都会影响压力,而压力会影响到各处的变形,因此我们的雅可比矩阵变得非常密集。使用默认的极细网格设置求解此问题时,该模型将花费大约 12GB 的内存和 400s 的时间来求解此具有约 22000 个自由度的模型。

“面载荷设置”窗口显示修改后的压力作为对破碎罐模型的输入
屏幕截图显示了压碎容器模型中修改后的压力。

但是,通过简单的将 nojac() 算子包裹到我们所施加的压力上,我们就可以忽略非局部耦合,从而降低内存需求和求解时间。这样,我们只需要 3GB 的内存和 10s 即可收敛到相同的解!

结论

默认情况下,COMSOL 软件将始终采用精确的雅可比矩阵,因为这可以提供最大的收敛可能性。通常,我们宁愿牺牲一些时间和内存来获得尽可能可靠的解。但是,我们确实有能力和灵活性来忽略雅可比矩阵中的选定项。如果要使用此功能,我们必须显式修改材料属性,修改基本的控制方程或在模型内编写自己的变量表达式,因此这是一个有点高级的用户功能。

当然,使用此功能将近似于雅可比矩阵,这有时会减慢甚至阻止收敛,因此在这种情况下请不要使用它。但是,如果模型应用 nojac() 算子时确实收敛,收敛于完整雅可比矩阵计算相同的解,并且可能会使用更少的时间和内存,那么在构建复杂的多物理场模型时,这是一个很好的特性。


评论 (4)

正在加载...
绍宽 毛
绍宽 毛
2022-09-14

您好,我在设置模型时,并没有设置nojac()算子,但是计算时报错总是显示这个算子中的分母有一个零项,而且我通过变量名搜索并未发现改项的变量名

Haoze Wang
Haoze Wang
2022-09-16 COMSOL 员工

您好,nojac算子可以通过避免填充雅可比矩阵来显著降低内存需求,部分模型中会内置使用,例如k-ε湍流模型,其中内置使用nojac算子提高性能,您可以在方程视图中检索到该算子。

磊 张
磊 张
2023-06-02

您好,我在计算的过程中显示如下错误:无法计算分离步骤 2 的初始残差和无法计算算子:nojac,导致无法计算临时符号衍生变量及表达式,但是我在无法计算表达式错误下搜索不到此算子表达式,请问我这样的情况该如何解决呢?

Qihang Lin
Qihang Lin
2023-06-05 COMSOL 员工

模型具体问题请联系技术支持:http://cn.comsol.com/support

浏览 COMSOL 博客