加入逆向运动粒子后的仿真

文中使用的是分层截断T样条，但是使用此方法定义的基函数不在单元上，而是跨越由几个单元组成的参数空间，因而编程复杂且无法嵌入现有的有限元法计算框架及相应算法 。文中使用了Bézier提取，将 T样条函数分解成伯恩斯坦多项式的线性组合 ，从而实现把 NURBS 单元分解为 C 0 C^0 C0 连续 Bézier 单元 ，这些单元与 Lagrange 单元相似 ，以便程序实现。

经过推导可以知道，以得到T样条基函数和Bézier 基函数（伯恩斯坦多项式）的关系为 N ( ξ ) = C B ( ξ ) N(\xi)=CB(\xi) N(ξ)=CB(ξ)，其中 C C C是Bézier 提取矩阵。

而程序的大致流程是，先使用伯恩斯坦求导，再利用伯恩斯坦基函数与T样条基函数的转换关系得到T样条基函数，以及其导数、曲线的表达和曲线的导数。之后进行的操作都是基于T样条曲线以及基函数进行的。其中T样条使用的控制顶点个数为48个。

程序中下面的代码就是由Bézier 基函数基函数得到T样条基函数的过程。其中cmat是转换矩阵。

在求解扩散方程的transport_src代码中，求解的扩散方程为
{ ∂ n 0 ∂ t − D ∇ 2 n 0 = − ( k + + k − ) n 0 + k + ′ n + + k − ′ n − in  Ω ∂ n + ∂ t + v + ⋅ ∇ n − = k + n 0 − k + ′ n + in  Ω ∂ n − ∂ t + v − ⋅ ∇ n − = k − n 0 − k − ′ n − in  Ω (1) \left\{\begin{array}{ll} \frac{\partial n_{0}}{\partial t}-D \nabla^{2} n_{0}=-\left(k_{+}+k_{-}\right) n_{0}+k_{+}^{\prime} n_{+}+k_{-}^{\prime} n_{-} & \text {in } \Omega \\ \frac{\partial n_{+}}{\partial t}+v_{+} \cdot \nabla n_{-}=k_{+} n_{0}-k_{+}^{\prime} n_{+} & \text {in } \Omega \\ \frac{\partial n_{-}}{\partial t}+v_{-} \cdot \nabla n_{-}=k_{-} n_{0}-k_{-}^{\prime} n_{-} & \text {in } \Omega \end{array}\right. \tag{1} ⎩⎨⎧​∂t∂n0​​−D∇2n0​=−(k+​+k−​)n0​+k+′​n+​+k−′​n−​∂t∂n+​​+v+​⋅∇n−​=k+​n0​−k+′​n+​∂t∂n−​​+v−​⋅∇n−​=k−​n0​−k−′​n−​​in Ωin Ωin Ω​(1)

使用SUPG方法之后，将对流项部分的权函数换为SUPG方法中求得的稳定参数 τ \tau τ然后积分，得到式（2）

∫ Ω ω h ∂ n 0 h ∂ t d Ω + ∫ Ω ∇ ω h ⋅ D ∇ n 0 h d Ω + ∫ Ω ω h ( k + + k − ) n 0 h d Ω − ∫ Ω ω h ( k + ′ n + h + k − ′ n − h ) d Ω = 0 ∫ Ω ω h ( ∂ n + h ∂ t + v + h ⋅ ∇ n + h − k + n 0 h + k + ′ n + h ) d Ω + ∑ i = 1 n e l e ∫ Ω i τ S U P G v + h ⋅ ∇ ω h ( ∂ n + h ∂ t + v + ⋅ ∇ n + h − k + n 0 h + k + ′ n + h ) d Ω = 0 ∫ Ω ω h ( ∂ n − h ∂ t + v − h ⋅ ∇ n − h − k − n 0 h + k − ′ n − h ) d Ω + ∑ i = 1 n e l e ∫ Ω i τ S U P G v − h ⋅ ∇ ω h ( ∂ n − h ∂ t + v − ⋅ ∇ n − h − k − n 0 h + k − ′ n − h ) d Ω = 0 (2) \int_{\Omega} \omega^{h} \frac{\partial n_{0}^{h}}{\partial t} d \Omega+\int_{\Omega} \nabla \omega^{h} \cdot D \nabla n_{0}^{h} d \Omega+\int_{\Omega} \omega^{h}\left(k_{+}+k_{-}\right) n_{0}^{h} d \Omega-\int_{\Omega} \omega^{h}\left(k_{+}^{\prime} n_{+}^{h}+k_{-}^{\prime} n_{-}^{h}\right) d \Omega=0 \\ \int_{\Omega} \omega^{h}\left(\frac{\partial n_{+}^{h}}{\partial t}+v_{+}^{h} \cdot \nabla n_{+}^{h}-k_{+} n_{0}^{h}+k_{+}^{\prime} n_{+}^{h}\right) d \Omega+\sum_{i=1}^{n_{e l e}} \int_{\Omega_{i}} \tau_{S U P G} v_{+}^{h} \cdot \nabla \omega^{h}\left(\frac{\partial n_{+}^{h}}{\partial t}+v_{+} \cdot \nabla n_{+}^{h}-k_{+} n_{0}^{h}+k_{+}^{\prime} n_{+}^{h}\right) d \Omega=0 \\ \int_{\Omega} \omega^{h}\left(\frac{\partial n_{-}^{h}}{\partial t}+v_{-}^{h} \cdot \nabla n_{-}^{h}-k_{-} n_{0}^{h}+k_{-}^{\prime} n_{-}^{h}\right) d \Omega+\sum_{i=1}^{n_{e l e}} \int_{\Omega_{i}} \tau_{S U P G} v_{-}^{h} \cdot \nabla \omega^{h}\left(\frac{\partial n_{-}^{h}}{\partial t}+v_{-} \cdot \nabla n_{-}^{h}-k_{-} n_{0}^{h}+k_{-}^{\prime} n_{-}^{h}\right) d \Omega=0 \tag{2} ∫Ω​ωh∂t∂n0h​​dΩ+∫Ω​∇ωh⋅D∇n0h​dΩ+∫Ω​ωh(k+​+k−​)n0h​dΩ−∫Ω​ωh(k+′​n+h​+k−′​n−h​)dΩ=0∫Ω​ωh(∂t∂n+h​​+v+h​⋅∇n+h​−k+​n0h​+k+′​n+h​)dΩ+i=1∑nele​​∫Ωi​​τSUPG​v+h​⋅∇ωh(∂t∂n+h​​+v+​⋅∇n+h​−k+​n0h​+k+′​n+h​)dΩ=0∫Ω​ωh(∂t∂n−h​​+v−h​⋅∇n−h​−k−​n0h​+k−′​n−h​)dΩ+i=1∑nele​​∫Ωi​​τSUPG​v−h​⋅∇ωh(∂t∂n−h​​+v−​⋅∇n−h​−k−​n0h​+k−′​n−h​)dΩ=0(2)

在SUPG方法中权函数选择的是 N + τ v ∇ N N+\tau v \nabla N N+τv∇N，即式（2）中的后两项中 w h = N + τ v ∇ N w^h=N+\tau v \nabla N wh=N+τv∇N，其中 v v v是速度，可取 v + , v − v_+,v_- v+​,v−​。而第一项中的权函数依旧取基函数，即 w h = N w^h=N wh=N

根据导数定义，当 Δ t \Delta t Δt足够小时 ∂ n h ∂ t \frac{\partial n^{h}}{\partial t} ∂t∂nh​可写成 n h ( t − 1 ) − n h ( t ) Δ t \frac{n^{h}(t-1)-n^{h}(t)}{\Delta t} Δtnh(t−1)−nh(t)​，整理之后的单元上的线性方程组为
[ ∑ i = 0 n e n ∑ j = 0 n e n ( 1 + Δ t ( k + + k − ) N i N j + Δ t D ∇ N i ∇ N j ) ∑ i = 0 n e n ∑ j = 0 n e n − k + ′ Δ t N i N j ∑ j = 0 n e n − k + ′ Δ t w + h Δ N j ∑ j = 0 n e n ( 1 + Δ t k + ′ ) w + h N j + Δ t w + h v + ∇ N j ] [ n 0 h ( t ) n + h ( t ) ] = [ n 0 ( t − 1 ) N i n + ( t − 1 ) w + h ] (3) \begin{bmatrix} \sum\limits_{i=0}^{nen}\sum\limits_{j=0}^{nen}(1+\Delta t(k_++k_-)N_iN_j+\Delta tD\nabla N_i \nabla N_j) & \sum\limits_{i=0}^{nen}\sum\limits_{j=0}^{nen}-k^{\prime}_+\Delta t N_iN_j \\ \sum\limits_{j=0}^{nen}-k^{\prime}_+\Delta tw^h_+\Delta N_j & \sum\limits_{j=0}^{nen}(1+\Delta tk^{\prime}_+)w^h_+N_j+\Delta tw^h_+v_+\nabla N_j \end{bmatrix} \begin{bmatrix} n_{0}^{h}(t) \\ n_{+}^{h}(t) \end{bmatrix} =\begin{bmatrix} n_0(t-1)N_i\\ n_+(t-1)w^h_+ \end{bmatrix} \tag{3} ⎣⎢⎢⎡​i=0∑nen​j=0∑nen​(1+Δt(k+​+k−​)Ni​Nj​+ΔtD∇Ni​∇Nj​)j=0∑nen​−k+′​Δtw+h​ΔNj​​i=0∑nen​j=0∑nen​−k+′​ΔtNi​Nj​j=0∑nen​(1+Δtk+′​)w+h​Nj​+Δtw+h​v+​∇Nj​​⎦⎥⎥⎤​[n0h​(t)n+h​(t)​]=[n0​(t−1)Ni​n+​(t−1)w+h​​](3)

式中 N N N是T样条基函数； Δ t \Delta t Δt是时间变化率这里取0.1， t − 1 t-1 t−1是时间 t t t的上一个状态。

代码实现如下

右端项的代码实现如下

每次使用到的 n 0 , n + n_0,n_+ n0​,n+​都是上一次迭代求解得到的，当然，第一次迭代时候使用的是模拟参数中的初始值。

模型文件2M，花费时间4个小时，耗费内存40GB

求解出的速度场结果

得到速度场后进一步求解运输模型。在对时间采样次数为100次时，得到的结果如下

在阅读代码后发现，虽然文中提出的模型有涉及对逆向运动粒子的计算，但代码中并没有对其进行体现。如果考虑逆向运动粒子的作用，则相应的线性系统将表示为

[ ∑ i = 0 n e n ∑ j = 0 n e n ( 1 + Δ t ( k + + k − ) N i N j + Δ t D ∇ N i ∇ N j ) ∑ i = 0 n e n ∑ j = 0 n e n − k + ′ Δ t N i N j ∑ i = 0 n e n ∑ j = 0 n e n − k − ′ Δ t N i N j ∑ j = 0 n e n − k + ′ Δ t w + h Δ N j ∑ j = 0 n e n ( 1 + Δ t k + ′ ) w + h N j + Δ t w + h v + ∇ N j 0 ∑ j = 0 n e n − k − ′ Δ t w − h Δ N j 0 ∑ j = 0 n e n ( 1 + Δ t k − ′ ) w − h N j + Δ t w − h v − ∇ N j ] [ n 0 h ( t ) n + h ( t ) n − h ( t ) ] = [ n 0 ( t − 1 ) N i n + ( t − 1 ) w + h n − ( t − 1 ) w − h ] (4) \begin{bmatrix} \sum\limits_{i=0}^{nen}\sum\limits_{j=0}^{nen}(1+\Delta t(k_++k_-)N_iN_j+\Delta tD\nabla N_i \nabla N_j) & \sum\limits_{i=0}^{nen}\sum\limits_{j=0}^{nen}-k^{\prime}_+\Delta t N_iN_j &\sum\limits_{i=0}^{nen}\sum\limits_{j=0}^{nen}-k^{\prime}_-\Delta t N_iN_j \\ \sum\limits_{j=0}^{nen}-k^{\prime}_+\Delta tw^h_+\Delta N_j & \sum\limits_{j=0}^{nen}(1+\Delta tk^{\prime}_+)w^h_+N_j+\Delta tw^h_+v_+\nabla N_j & 0 \\ \sum\limits_{j=0}^{nen}-k^{\prime}_-\Delta tw^h_-\Delta N_j & 0 & \sum\limits_{j=0}^{nen}(1+\Delta tk^{\prime}_-)w^h_-N_j+\Delta tw^h_-v_-\nabla N_j \end{bmatrix} \begin{bmatrix} n_{0}^{h}(t) \\ n_{+}^{h}(t) \\ n_{-}^{h}(t) \end{bmatrix} =\begin{bmatrix} n_0(t-1)N_i \\ n_+(t-1)w^h_+ \\ n_-(t-1)w^h_- \end{bmatrix} \tag{4} ⎣⎢⎢⎢⎢⎢⎢⎡​i=0∑nen​j=0∑nen​(1+Δt(k+​+k−​)Ni​Nj​+ΔtD∇Ni​∇Nj​)j=0∑nen​−k+′​Δtw+h​ΔNj​j=0∑nen​−k−′​Δtw−h​ΔNj​​i=0∑nen​j=0∑nen​−k+′​ΔtNi​Nj​j=0∑nen​(1+Δtk+′​)w+h​Nj​+Δtw+h​v+​∇Nj​0​i=0∑nen​j=0∑nen​−k−′​ΔtNi​Nj​0j=0∑nen​(1+Δtk−′​)w−h​Nj​+Δtw−h​v−​∇Nj​​⎦⎥⎥⎥⎥⎥⎥⎤​⎣⎡​n0h​(t)n+h​(t)n−h​(t)​⎦⎤​=⎣⎡​n0​(t−1)Ni​n+​(t−1)w+h​n−​(t−1)w−h​​⎦⎤​(4)

下面我模仿着对正向运动粒子的运动，在代码中加入了对逆向运动粒子的计算。其中，刚度矩阵的计算部分代码如下

修改之后的模拟中，逆向运动粒子的初始浓度 n − n_- n−​设置为0.2，其他模拟参数与之前完全相同。对时间进行100次采样，得到的结果如下所示。

从结果上看，仿真的运输结果与考虑逆向运动粒子之前的结果差别不大，可能是加入的初始逆向运动粒子浓度不高（0.2），所以结果相差不明显。

在仿真程序运行过程中，需要消耗大量资源与时间。以该模型为例，仿真所需的时间加起来需要10个小时，消耗内存达50GB！！而模型文件只有2M。