模型降阶方法之张量方法应用举例

序章

POD 方法参考：
https://blog.csdn.net/lusongno1/article/details/125944587

DEIM 方法参考：
https://blog.csdn.net/lusongno1/article/details/125955245

张量方法参考：
https://blog.csdn.net/lusongno1/article/details/126013084

很多同学看完了张量方法表示，算法流程是 clear 的，但是并不知道怎么应用。这里举两个简单的例子，热传导方程和对流扩散方程。

我们可以看到，张量方法是应用于 ODEs 的，那么如果从 PDE 得到 ODEs 呢？简单地说， PDE 在空间上做离散，比如说有限元离散或者有限差分离散，那么就把问题转为了向量值 ODEs，就可以用 ROMs 方法去处理。 如果用的是有限元方法处理的，那么向量值解的每个分量其实就表示每个基函数前面的系数（自由度）。如果是用有限差分方法处理的，那么向量值解的每个分量其实就表示每个节点上的值。

下面我们考虑这两个例子。有一些符号的含义是显而易见的，所以就不多加解释了。

热传导方程

方法

考虑热传导方程，
$$w_t=\Delta w$$
它的边界条件是，
$${(\mathbf{n}\cdot \nabla w+q(\mathbf{x},\mathbf{\alpha })w)|}_{\partial \Omega }=g(\mathbf{x},\mathbf{\alpha })$$
其中，$\mathbf{\alpha }\in {\mathbb{R}}^4$ 是个四维的参数。这里的求解区域为长方形挖三个洞，

用标准的有限元方法，在空间上做离散，就得到这样一个时间连续的方程，
$$\mathrm{M}{\mathbf{u}}_t+(\mathrm{K}+\mathrm{Q}(\mathbf{\alpha }))\mathbf{u}=\mathbf{g}(\mathbf{\alpha })$$
其中，$M$ 是质量矩阵，$K$ 是刚度矩阵，
$$(\mathrm{Q}{)}_{ij}(\mathbf{\alpha })={\int }_{\partial \Omega }q(\mathbf{x},\mathbf{\alpha }){\theta }_j(\mathbf{x}){\theta }_i(\mathbf{x})d{s}_{\mathbf{x}},(\mathbf{g}{)}_j(\mathbf{\alpha })={\int }_{\partial \Omega }g(\mathbf{x},\mathbf{\alpha }){\theta }_j(\mathbf{x})d{s}_{\mathbf{x}},$$
那么，我们就可以将其写的 ODEs 的形式，
$${\mathbf{u}}_t=F(t,\mathbf{u},\mathbf{\alpha })$$
其中，
$$F(t,\mathbf{u},\mathbf{\alpha })=-{\mathrm{M}}^{-1}(\mathrm{K}+\mathrm{Q}(\mathbf{\alpha }))\mathbf{u}+{\mathrm{M}}^{-1}\mathbf{g}(\mathbf{\alpha })$$
有了这个，我们就可以用所谓的张量方法来寻求 reduced basis 了。首先，我们用 CN 格式得到一些快照，
$${\mathbf{\varphi }}_k=\mathbf{u}\left(t_k,\mathbf{\alpha }\right)$$
函数解事实上可由解向量表示出来，即
$$w\left(t_k,\mathbf{x},\mathbf{\alpha }\right)=\Theta (\mathbf{x})\mathbf{u}\left(t_k,\mathbf{\alpha }\right)$$
用之前提到的张量方法（不管是 CP、HOSVD 还是 TT），一波操作之后，我们便得到了一组 reduced 基，
$$\mathrm{Z}=\left[{\mathbf{z}}_1,\dots ,{\mathbf{z}}_n\right]\in {\mathbb{R}}^{M\times n}$$

不清楚怎么得到的，回看张量方法。事实上，$\mathrm{Z}={\mathrm{UU}}_c$，${\mathrm{U}}_c$ 取的是左奇异向量的前几列。

接着，类似于 POD 方法的投影方式，我们可以把这个方程投影到低阶的子空间中去做，得到
$$\tilde{\mathrm{M}}{\tilde{\mathbf{u}}}_t+(\tilde{\mathrm{K}}+\tilde{\mathrm{Q}}(\mathbf{\alpha }))\tilde{\mathbf{u}}=\tilde{\mathbf{g}}(\mathbf{\alpha })$$
其中，
$$\begin{array}{rl}& \tilde{\mathrm{M}}={\mathrm{Z}}^T\mathrm{MZ}\in {\mathbb{R}}^{n\times n},\tilde{\mathrm{K}}={\mathrm{Z}}^T\mathrm{KZ}\in {\mathbb{R}}^{n\times n}\text{, }\\ & \tilde{\mathrm{Q}}(\mathbf{\alpha })={\mathrm{Z}}^T\mathrm{Q}(\mathbf{\alpha })\mathrm{Z}\in {\mathbb{R}}^{n\times n},\tilde{\mathbf{g}}(\mathbf{\alpha })={\mathrm{Z}}^T\mathbf{g}(\mathbf{\alpha })\in {\mathbb{R}}^n\text{, }\end{array}$$
这个问题就被简化了，得到了低维空间的解，要恢复原来的函数，可以这样做：
$$\tilde{w}\left(t_k,\mathbf{x},\mathbf{\alpha }\right)=\Theta (\mathbf{x})\mathrm{Z}\tilde{\mathbf{u}}\left(t_k,\mathbf{\alpha }\right)\approx w\left(t_k,\mathbf{x},\mathbf{\alpha }\right)$$

数值结果

下面几乎不加解释变量定义（下同）地 post 一些数值实验的结果，变量感兴趣的同学可以参考论文：

Mamonov A V, Olshanskii M A. Interpolatory tensorial reduced order models for parametric dynamical systems[J]. Computer Methods in Applied Mechanics and Engineering, 2022, 397: 115122.

这两个例子都是从这篇论文中抄出来的。

上述的误差定义为：
$$E_{L^2(\hat{\mathcal{A}})}={\left(\frac{1}{MNK}\sum _{j=1}^K{\Vert \left(\mathrm{I}-{\mathrm{ZZ}}^T\right){\Phi }_e\left({\hat{\mathbf{\alpha }}}_j\right)\Vert }_F^2\right)}^{1/2}$$
$$E_{L^{\infty }(\mathcal{A})}=\underset{\mathbf{\alpha }\in \mathcal{A}}{\sup }{\left(\frac{1}{MN}{\Vert \left(\mathrm{I}-{\mathrm{ZZ}}^T\right){\Phi }_e(\mathbf{\alpha })\Vert }_F^2\right)}^{1/2}$$

压缩因子 CF 的定义为，
$$\mathrm{CF}=\frac{\#(\mathbf{\Phi })}{\#(\text{ online }(\tilde{\mathbf{\Phi }}))}$$
$\#$ 指代的是存储所需要的浮点数目。别的就，hmm，看文章，文章里面写得很清楚。

简单地说，TABLE 5.1~5.2 笛卡尔网格方法和通常参数采样方法的投影误差和离线在线阶段需要传输的内存衡量的。
FIG 5.2 说的是各种方法下（不同分解、笛卡尔网格和通常参数）压缩因子 CF 和 canonical 秩随着分解 tolerance 的变化而变化的图。
FIG 5.3 说的是样本外预测的投影误差随着 reduced basis 的个数、参数的采样间隔、张量分解精度的变化而变化的情况。这里的 $p$ 表示的是选择该维度上最靠近采样点的 $p$ 个点。
TABLE 5.3-5.5 不同的参数离散、不同的 reduced basis 维数、不同的投影方式下，不同方法的相对增益，相对增益指的是相对于 POD 而言的相对误差，
$$G_{\mathrm{X}}^{(r)}=\frac{R_{\mathrm{POD}}\left({\mathbf{\alpha }}^{(r)}\right)}{R_{\mathrm{X}}\left({\mathbf{\alpha }}^{(r)}\right)},r=1,2,\dots ,N_r$$
而相对误差不再是投影的误差，而是回复到原来的函数解所定义的误差，
$$\begin{array}{rl}{R}_X(\mathbf{\alpha })& =\frac{{\max }_{k=1,\dots ,N}{\Vert \tilde{w}\left(t_k,\mathbf{x},\mathbf{\alpha }\right)-w\left(t_k,\mathbf{x},\mathbf{\alpha }\right)\Vert }_{L^2(\Omega )}}{{\max }_{k=1,\dots ,N}{\Vert w\left(t_k,\mathbf{x},\mathbf{\alpha }\right)\Vert }_{L^2(\Omega )}}\\ & \approx \frac{{\sup }_{t\in [0,T]}\Vert \tilde{w}(t,\mathbf{x},\mathbf{\alpha })-w(t,\mathbf{x},\mathbf{\alpha }){\Vert }_{L^2(\Omega )}}{{\sup }_{t\in [0,T]}\Vert w(t,\mathbf{x},\mathbf{\alpha }){\Vert }_{L^2(\Omega )}}\end{array}$$

一般我们考虑投影误差，或者恢复出来的函数误差就可以了。更多细节处的解释，建议读论文。

对流扩散方程

方法

类似于热传导方程，我们下面考虑对流扩散方程，
$$w_t=\nu \Delta w-\mathbf{\eta }(\mathbf{x},\mathbf{\alpha })\cdot \nabla w+f(\mathbf{x})$$
其中，
$$f(\mathbf{x})=\frac{1}{2\pi {\sigma }_s^2}\exp \left(-\frac{{\left(x_1-x_1^s\right)}^2+{\left(x_2-x_2^s\right)}^2}{2{\sigma }_s^2}\right)$$
求解区域，
$$\Omega =[0,1]\times [0,1]\subset {\mathbb{R}}^2$$
初边值条件，
$${(\mathbf{n}\cdot \nabla w)|}_{\partial \Omega }=0,w(0,\mathbf{x},\mathbf{\alpha })=0$$
$$\mathbf{\eta }(\mathbf{x},\mathbf{\alpha })=\left(\begin{array}{c}{\eta }_1(\mathbf{x},\mathbf{\alpha })\\ {\eta }_2(\mathbf{x},\mathbf{\alpha })\end{array}\right)=\left(\begin{array}{c}\cos {\alpha }_9\\ \sin {\alpha }_9\end{array}\right)+\frac{1}{\pi }\left(\begin{array}{c}{\partial }_{x_2}h(\mathbf{x},\mathbf{\alpha })\\ -{\partial }_{x_1}h(\mathbf{x},\mathbf{\alpha })\end{array}\right)$$
$$\begin{array}{rl}h(\mathbf{x},\mathbf{\alpha })=& {\alpha }_1\cos \left(\pi x_1\right)+{\alpha }_2\cos \left(\pi x_2\right)+{\alpha }_3\cos \left(\pi x_1\right)\cos \left(\pi x_2\right)\\ & +{\alpha }_4\cos \left(2\pi x_1\right)+{\alpha }_5\cos \left(2\pi x_2\right)+{\alpha }_6\cos \left(2\pi x_1\right)\cos \left(\pi x_2\right)\\ & +{\alpha }_7\cos \left(\pi x_1\right)\cos \left(2\pi x_2\right)+{\alpha }_8\cos \left(2\pi x_1\right)\cos \left(2\pi x_2\right).\end{array}$$
可以看得出来，这里实际上有 9 个参数，用有限元方法对空间导数做离散，我们得到类似于热传导的空间离散方程，其中 $H$ 是对流矩阵，
$$M{\mathbf{u}}_t+(K+H(\mathbf{\alpha }))\mathbf{u}=\mathbf{f}$$
那么，ODEs 的右端项为，
$$F(t,\mathbf{u},\mathbf{\alpha })=-{\mathrm{M}}^{-1}(\mathrm{K}+\mathrm{H}(\mathbf{\alpha }))\mathbf{u}+{\mathrm{M}}^{-1}\mathbf{f}$$
同样地，用某种格式，比如 CN 格式得到 snapshots，之后用张量方法得到 reduced basis，再对原方程进行降阶，得到了，
$$\tilde{\mathrm{M}}{\tilde{\mathbf{u}}}_t+(\tilde{\mathrm{K}}+\tilde{\mathrm{H}}(\mathbf{\alpha }))\tilde{\mathbf{u}}=\tilde{\mathbf{f}}$$

数值结果

同样地，我们也不加解释变量含义地看一些结果，

各变量含义同上个例子解释。