LongLong's Blog

08 November 2022

拓扑优化作为结构设计中优化设计的一种自动化的手段，结合增材制造技术可以实现结构的超轻量化。但由于拓扑优化的过程中需要反复进行有限元的求解，计算量极大较难实现高分辨率模型的优化。

taichi语言的作者胡渊鸣在文章Narrow-Band Topology Optimization on a Sparsely Populated Grid中给出了一种基于稀疏多重网格预条件共轭梯度法的计算方法，其可以大幅节约计算资源实现在单台PC上进行高分辨率的拓扑优化，并给出了程序的源代码，这里主要对该代码的配置和使用做简要介绍。

1. 拓扑优化

拓扑优化（topology optimization）是一种根据给定的负载情况、约束条件和性能指标，在给定的区域内对材料分布进行优化的数学方法，是结构优化的一种。其本质上就是求一个材料的分布函数，使得在给定的边界条件下结构达到某种最优的状态，一般情况下会让结构的变形能最小（一定意义上的刚度最大化），优化的过程也是基于梯度下降的方式——沿着梯度方向更新待求函数，同时为了避免出现棋盘效应，还需要对相邻的位置进行平均化。

由于优化的目标函数一般都是节点位移的函数，为此每一步优化都需要进行一次完成的有限元求解过程——计算单元刚度矩阵，集成总体刚度矩阵，求解方程组，对于规模较大的问题，进行一次求解就需要花费大量的时间，导致拓扑优化的网格分辨率无法提升到很高，从而无法对结构的微观结构进行进一步优化。

随着稀疏存储技术和多重网格预条件方法的引入，使得计算更大规模的拓扑优化问题成为可能，以此优化后的结构在细观上会呈现出一定的仿生特征。

2. 运行环境

运行环境的核心是编译遗产(legacy)版本的taichi和代码中附带的一个SPGrid求解器，其用到了Intel编译器和MKL的一些特性，为此还需要安装Intel Parallel Studio Xe。

代码说明中中推荐的系统版本为Ubuntu18.04，经测试确实是可以的，而且20.04之后的版本都是不行的，同时其建议的Intel Parallel Studio Xe为2018版，但实际中如果使用2018版会发现其不能识别Ubuntu18.04，编译会因无法识别系统而报错，虽然可以通过一些方法解决，但如果使用2019版本的Intel Parallel Studio Xe则不会存在这个问题。

目前Intel的官网已经找不到Intel Parallel Studio Xe的下载地址，取而代之的是Intel oneAPI，经测试二者并不能兼容，不过幸好我还是找到了Intel Parallel Studio Xe的下载地址，而且该地址目前依然可用，但其为商业版需要Intel提供的序列号或是证书文件。 Intel Parallel Studio Xe 2019。

3. 安装taichi遗产版

可通过提供的Python脚本之间自动化安装，需要提前安装好python3和git，但实际运行中还会缺少很多依赖库，可在报错后安装即可。

wget https://raw.githubusercontent.com/yuanming-hu/taichi/legacy/install.py
python3 install.py

4. 编译SPGrid Solver

首先将spgrid_topo_opt的代码放置到taichi源代码目录下的projects目录下，配置好Intel编译器的配置后，修改文件SPGrid_SIMD_Utilities.h，在其include部分增加

#include <immintrin.h>

同时还可以去掉Makefile中的-DHASWELL，避免编译中产生大量重复定义的报错，然后就可以开始编译SPGrid Solver，如果CPU支持AVX512（目前能够完整支持的普通桌面CPU应该只有Intel的11代Core）则应该可以使用Makefile中被注释掉的一行进行编译

source /opt/intel/bin/compilervars.sh -arch intel64
cd solver
make

5. 编译其余部分

其余部分将作为taichi的一部分，退回到spgrid_topo_opt目录下使用ti build指令进行编译，编译前需要定义三个环境变量，分别为MLK路径，GPU对应的CUDA架构版本（没有可以设为0，具体值可以在官网查表），使用double类型（这个必须设置为1否则后面编译会报错），同时如果需要使用对计算结果的可视化功能必须具有CUDA支持才可以

export TC_MKL_PATH=/opt/intel/compilers_and_libraries_2019/linux/mkl/lib/intel64_lin/
#Pascal架构为6.1，Turing架构为7.5，Ampere架构为8.6
export CUDA_ARCH=75
export TC_USE_DOUBLE=1

然后就可以运行scripts目录中的例子了，需要注意这里的例子都需要较大的内存才能运行，结果文件会被保存到taichi/outputs目录下

6. 基本使用方法

由于没有详细的说明文档，使用的方法主要可以参考scripts目录下的相关实例脚本，同时编写的脚本也需要放置于scripts目录下。

(1) 导出模型

该代码隐含要求设计域为(-0.5, -0.5, -0,5)到(0.5, 0.5, 0.5)组成体积为1的立方体区域，为此在导入设计模型是需要将设计的三维模型缩放至在导出模型的坐标系中刚好在该包络内，导出模型的格式为Wavefront的obj格式，导出方式类似stl并且需要在导出时指定坐标系。导出的模型放置于data目录下。

(2) 参数设置

求解主要的参数为网格的分辨率和体积分数，另外一般可以开启保留位移边界和载荷施加位置的单元的设置，其他参数的设置可参考其他实例脚本中的设置

from topo_opt import TopoOpt;
version = 1;
narrow_band = True;
#体积分数
volume_fraction = 0.1;
#网格分辨率
n = 3000;
#fix_cells_at_dirichlet和fix_cells_near_force分别为保留位移边界和载荷施加位置的单元
opt = TopoOpt(res=(n, n, n), version=version, volume_fraction=volume_fraction,
              grid_update_start=5 if narrow_band else 1000000,
              fix_cells_near_force=True, fix_cells_at_dirichlet=True, fixed_cell_density=0.1);

(3) 施加边界条件

通过add_dirichlet_bc和add_load及其他的衍生方法施加位移及载荷边界条件，目前支持以点为中心的区域，平面区域，长方体区域，网格表面等施加方法。

#施加位移边界条件
opt.add_dirichlet_bc((0, 0, -0.1), radius=0.05, axis='xyz', value=(0, 0, 0));
#施加载荷边界条件
opt.add_load(center=(0.0, -y + 0.01, 0.0), force=(0, -1e6, 0));

(4) 结果可视化

计算的结果将保存到taichi/outputs目录下，使用ti vd命令可以实现对计算结果的可视化，需要配置好CUDA支持才可以使用

#对fem目录下的结果文件进行可视化显示
ti vd fem

使用1,2,3键可以对结果进行对三个坐标轴的对称显示，H和L键可以在多个求解结果文件之间进行切换，J和K键可调节显示部分的最小密度，T键可进行透明化显示（类似X光片），S及Shift+S可对可视化结果进行保存（保存为png或是ply格式到/tmp目录下）。

04 February 2022

PyTorch作为一款用于深度学习的计算框架，其不但提供了常用的各种深度学习模型，同时提供了多维数组（张量）及运算，线性代数和自动微分相关的特性。使得其不但可以用于深度神经网络的建模，同时还可以直接作为优化器使用，其实深度神经网络训练的过程就是对其参数优化的过程。

1. 自动微分

微分操作可以通过手动微分、数值微分、符号微分、自动微分实现，其中自动微分通过计算过程中对变量的跟踪实现对导数的计算，在深度学习框架中使用的为反向传播自动微分。在PyTorch中需要跟踪计算导数的变量需要在声明是增加requires_grad，相当于明确该变量为函数的自变量而非参数变量

例如计算函数 \( f(x, y) = x^2 + 2xy \)在点(1, 2)处的导数（梯度） \[\nabla f\ = [\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}] = [2x + 2y, 2x]\] \[\nabla f|_{(1, 2)}=[6, 4]\] 使用PyTorch计算的程序如下

import torch;
#声明两个需要跟踪导数的变量
x = torch.tensor((1.), requires_grad=True);
y = torch.tensor((2.), requires_grad=True);
#计算函数值
z = x ** 2 + 2 * x * y;
#反向传播
z.backward();
#输出导数
print(x.grad, y.grad);

最终的输出为，与手动微分计算的结果一致

tensor(6.) tensor(2.)

2. 优化问题

所谓优化问题，本质就是求函数的极值问题，对于连续可微问题则极值点的导数为零，优化问题将转换为方程求根问题。常用的方法为梯度下降法，即沿着梯度下降的方向寻找下一个位置。PyTorch中提供了各种常用的优化算法，在torch.optim包中，其中包含了部分自适应步长的算法。

例如计算函数\( f(x, y) = x^2 + 2x + y^2 \)的极小值

import torch;

x = torch.tensor((1.), requires_grad=True);
y = torch.tensor((2.), requires_grad=True);
#定义优化器，x y为进行优化的变量
optim = torch.optim.Adam([x, y], lr=1e-1);
#迭代优化
for i in range(1, 2000):
    #计算函数值
    z = x ** 2 + 2 * x + y ** 2;
    #梯度置零
    optim.zero_grad();
    #反向传播
    z.backward();
    #跟新优化变量
    optim.step();
    print(x.item(), y.item(), z.item());

最终得到的输出为

-1.0 0.0 -1.0

对于更一般的优化问题，除了函数本身还会有一些约束条件，对于带有约束条件的问题，在数学分析中可通过增加拉格朗日乘子的方法将其转换为不含约束的优化问题。

然而如果将乘子也作为其中的一个优化变量会导致构造的带有乘子的函数梯度为零的位置不是带有乘子函数的极值，考虑以下问题

求函数的极值 \[ f(x, y) = x^2 + y^2 \] 约束条件： \[ x + y - 1 = 0 \] 若采用乘子法，则 \[L = x^2 + y^2 + \lambda(x + y - 1)\] \[\frac{\partial L}{\partial x} = 2x + \lambda = \frac{\partial L}{\partial y} = 2y + \lambda = \frac{\partial L}{\partial \lambda} = x + y - 1 = 0\] 得到\(x = y = \frac{1}{2}, \lambda = -1\)

而如果使用以下程序处理该问题，会出现无法收敛的情况，说明梯度为零的点并不是函数L的极值点

import torch;

x = torch.tensor((1.), requires_grad=True);
y = torch.tensor((2.), requires_grad=True);
a = torch.tensor((0.), requires_grad=True);
optim = torch.optim.Adam([x, y, a], lr=1e-1);
for i in range(1, 2000):
    z = x ** 2 + y ** 2 + a * (x + y - 1)
    optim.zero_grad();
    z.backward();
    optim.step();
    print(x.item(), y.item(), a.item(), z.item());

为此为了程序实现上的便利可以采用罚函数法，考虑以下函数 \[ L = x^2 + y^2 + M(x + y - 1)^2\] 其中M为一个非常大的正数，此时当约束条件不满足时，函数L就会大幅增大，迫使沿着梯度下降时会自动落入满足边界条件的区域内，如此处理后程序变为

import torch;

x = torch.tensor((0.), requires_grad=True);
y = torch.tensor((0.), requires_grad=True);
optim = torch.optim.Adam([x, y], lr=1e-1);
for i in range(1, 2000):
    #增加约束条件的罚函数
    z = x ** 2 + y ** 2 + 1e5 * (x + y - 1) ** 2;
    optim.zero_grad();
    z.backward();
    optim.step();
    print(x.item(), y.item(), z.item());

最后将得到结果，非常接近使用解析的方法得到的结果

0.49999749660491943 0.49999749660491943 0.49999749660491943