前言

在分析了pennes模型下的热传导的三维形式之后，本章主要针对一个具体的问题进行分析————时空步长和实际仿真需求之间的矛盾。

p.s. 经过对加速方面进行了一些简答的调研，加上分工的一些变化，简单地使用AMP进行了运算的加速，本系列中基本就不再赘述。

FDM和MC方法中的差分形式

显式差分FDM

在上一期生物热仿真(3)：由一维到三维，由一般到特殊中有提到：显式差分不存在依赖关系，可以直接由上一时刻某些坐标的温度推导下一时刻的温度，但是由于差分的稳定性条件，细微尺度（例如mm级）下仿真的时间步长将会是空间步长的指数倍，这对于实际工程来说是不能容忍的。举个例子来说，如果我将空间区域分割为1mm大小，准备进行1min的仿真，则需要进行 $60/(0.001)^2=6 \times 10^7$ 数量级次的仿真，这样的要求使得加速带来的性能提升失去了意义。

隐式差分FDM

另一方面，使用隐式差分或者半隐式差分难以避免的问题就是结果间的相互依赖，一般的解法是通过联立方程组，将问题转换为矩阵计算，最后将结果的求解转换为大矩阵乘法。尽管存在结果上的相互依赖，但由于矩阵运算对于GPU加速的良好适应性，一种可行的FDM算法思路诞生了：使用隐式/半隐式差分+GPU加速矩阵运算获得可接受的性能。

传统Random Walk实现的显式/隐式差分MC

在上一期文章中提到的MC方法，尽管具有所谓的“适用于不同形状的物体，易于并行计算”的特点，但其中对于概率的计算方式仍然是脱胎于FDM的差分化方法，最主要的不同点则在于和边界点建立联系这一想法。这一思路的基础在于，边界点的温度往往具有较稳定的解，例如第一类边界条件下的恒温、第二第三边界条件下较易推导的温度。

然而，当内部存在热源时，单纯的通过边界计算温度则明显存在信息不足的问题。此时，我们需要将内部的热源也视为边界，虽然此种处理的直观性不足，但是合理性是不存在问题的。

p.s. 如果需要对所有点进行迭代式的非稳态仿真，传统MC相较于FDM不存在性能优势。但是MC支持单点温度求取，应用场景更多元，同时受边界限制也确实较少。

小结

不管是FDM还是MC，根本的求解思路还是对微分方程的差分化迭代求解，Random Walk虽然是利用了MC的统计思路，但理论依据还是同样的微分方程。由于显式差分显然不符合实际需求，（半）隐式差分成为了FDM和MC中的主流，而隐式差分之间存在的依赖性问题导致似乎只能通过联立方程求解。

由于懒得重新起一套方程求解的架子，我思考了一下顾及了性能之后的对MC方法的修改，主要是两种思路：

对Random Walk进行某种改进：找到了一篇Boundary-dispatch的论文做了这项工作，后文将对比传统Random Walk和他的区别
先计算隐式差分中的显式部分，然后基于统计的思路做一个修正，使得误差能达到可接受范围：这部分我做了一下简单的数学推导，可以对各方向的概率进行一个修正。

Boundary-dispatch: 更直观的Random Walk

传统的Monte Carlo流程中，每个内部点发出N个粒子进行Random Walk，最终粒子走到边界点就会被吸收，使得内部点和对应的边界点建立温度上的联系。因为边界点的温度可以通过边界条件进行确定，所以这种迭代方式并不会引起误差的积累。Boundary-dispatch和传统MC的主要区别有以下两点，下面将分小章进行说明。

区别1：发出粒子的方式不同

Boundary-dispatch方式的Monte Carlo流程中，区别于从每个内部点发出粒子，改为了从每个边界点（内部热源视为边界）发出粒子。显而易见，如果每个内部点发出 $10^2$ 级别的粒子，则改为边界点后，每个边界点应该发出的粒子数要高出几个量级。

区别2：添加了粒子的分类

Boundary-dispatch方式对于粒子进行了分类，基于分类方式的不同又分为双目（Bispecies）和多目（Multispecies）。

双目形式下，粒子被分为高能和低能粒子， $\frac{T_b}{T_{max}}$ 比例的粒子被设为高能，每个内部点 $(i,j)$ 会记录经过它的高能粒子的数量 $H_{i,j}$ ，低能粒子的数量 $L_{i,j}$ ，最后换算为温度 $T_{i, j}=T_{\max } \frac{H_{i, j}}{H_{i, j}+L_{i, j}}$ 。

多目形式下，有 $m$ 个边界点时，粒子被分类为 $m$ 类，分别对应一个边界点。记录内部点 $(i,j)$ 被第 $k$ 类的粒子经过的次数为 $V^k_{i,j}$ ，第 $k$ 个边界点的温度为 $T^k_b$ ，最终的温度为 $T_{i,j}=\frac{\sum^M_{k=1}{T^k_b V^k_{i,j}}}{\sum^M_{k=1}{V^k_{i,j}}}$ 。

小结

这个方法更加直观地反映了热从热源或边界点向外发散的过程，也存在对于内部热源点的考虑不够充分的问题，不过思路十分有趣，可以参照着修改传统MC，进而获得更好的效果。

隐式差分的工程显式化以及修正

附录：进一步分析之后发现这个思路不可行，下面的可以不用看了。。。

如前面所提，只有隐式差分可能满足小尺度下的时间步长需求，但是隐式差分的相互依赖又限制了求解的便捷性。三维空间内，坐标 $(x,y,z)$ 上第 $k+1$ 次迭代的温度值为 $T^{k+1}_{x,y,z}$ ，存在:

$T^{k+1}_{x,y,z} = T^{k}_{x,y,z} + \Delta T^{k}_{x,y,z} \tag{1}$

而温度场计算的隐式差分迭代式为：
$T^{k+1}_{x,y,z} = \sum ^6_{i=1} {F_i T^{k+1}_{N_i}} + F_0 T^k_{x,y,z} \tag{2}$
其中有 $(x,y,z)$ 的邻居集合 $N=\{(x-1,y,z), (x+1,y,z), (x,y-1,z), (x,y+1,z), (x,y,z-1), (x,y,z+1)\}$ ， $N_i$ 表示第 $i$ 个邻居， $F_i$ 表示第 $i$ 个邻居的迭代概率，具体计算方法参见上一期。

显然，通过式 $(1)$ 可对式 $(2)$ 进行简单的化简，将变化剥离开来：
$T^{k+1}_{x,y,z} = \sum ^6_{i=1} {F_i T^{k}_{N_i}} + F_0 T^k_{x,y,z} + \sum ^6_{i=1} {F_i \Delta T^{k}_{N_i}} \tag{3}$
其中的 $\sum ^6_{i=1} {F_i \Delta T^{k}_{N_i}}$ 项就是这样简单粗暴的处理产生的误差，显然这样的误差在较大时间尺度上是不可接受的，因为较大时间尺度上 $\Delta T^{k}_{N_i}$ 将会是很大的值。