摘要:为提高传统光滑粒子动力学(SPH)方法求解高维非线性薛定谔(nonlinear Schr?dinger/Gross-Pitaevskii equation, NLS/GP)方程的数值精度和计算效率, 本文首先基于高阶时间分裂思想将非线性薛定谔方程分解成线性导数项和非线性项, 其次拓展一阶对称SPH方法对复数域上线性导数部分进行显式求解, 最后引入MPI并行技术, 结合边界施加虚粒子方法给出一种能够准确、高效地求解高维NLS/GP方程的高阶分裂修正并行SPH方法. 数值模拟中, 首先对带有周期性和Dirichlet边界条件的NLS方程进行求解, 并与解析解做对比, 准确地得到了周期边界下孤立波的奇异性, 且对提出方法的数值精度、收敛速度和计算效率进行了分析; 随后, 运用给出的高阶分裂粒子方法对复杂二维和三维NLS/GP问题进行了数值预测, 并与其他数值结果进行比较, 准确地展现了非线性孤立波传播中的奇异现象和玻色-爱因斯坦凝聚态中带外旋转项的量子涡旋变化过程.
关键词: 非线性薛定谔方程/
光滑粒子动力学/
时间分裂/
玻色-爱因斯坦凝聚态

English Abstract

全文HTML

--> --> -->

众所周知, 非线性薛定谔方程(nonlinear Schr?dinger/Gross-Pitaevskii equation, NLS/GP)常被用来描述一些非线性物理特性, 例如: 晶体中热脉冲、非线性光学、等离子和量子动力学现象等^[1-5]. NLS/GP方程的研究涉及复杂非线性项和外部旋转导数项, 导致孤立波值随时间变化过程在许多情况下难以用理论手段^[6-8]精确地得到. 已有许多数值方法被提出用于对高维NLS/GP方程进行模拟, 比如分裂格式有限差分法、有限元法、时间伪谱法、修正欧拉算法、蒙特卡罗法和无网格方法等^[3-14], 然而上述方法属于网格类或是基于背景网格, 对非均匀节点布置下或高维复杂区域问题的模拟实现都较复杂^[15-18]. 为此已有许多****对完全不依赖于网格的粒子方法(或纯无网格方法)进行了关注, 并成功将粒子方法推广应用到实数域上偏微分方程的求解, 如光滑粒子动力学(smoothed particle hydrodynamics, SPH)方法^[15-25]等. 然而上述纯无网格方法对非线性薛定谔方程的研究还处于起步阶段, 特别对高维复杂GP方程的研究在国际上尚不多见.
作为完全不依赖于网格的SPH粒子方法, 在高维非线性薛定谔方程上的模拟应用还处在初步研究阶段^[15], 其原因在于传统SPH方法的数值精度低、稳定性差^{[16-19,20-22]}. SPH方法精度和稳定性的改进与高维非线性薛定谔方程的数值研究是国际上的两个热点问题. 针对SPH方法数值精度的提高, 已有一些改进SPH方法^{[16-19,21-25]}被提出, 并被应用于许多领域^[16-25], 但它们仍存在一些自身缺陷, 因此当推广应用到非线性薛定谔方程的模拟时需要进一步改进. 粒子方法模拟高维NLS/GP方程较网格类方法的主要优势在于: 不依赖于网格可以任意布点; 便于处理复杂区域和计算空间导数; 模拟程序易于实现, 并且便于实施并行计算以提高计算效率.
基于上述分析, 本文针对含有非线性项和旋转导数项高维NLS/GP方程的准确、高效性粒子法模拟给出一种四阶分裂修正并行SPH (HSS-CPSPH)方法. 提出的HSS-CPSPH方法思想主要来源于: 首先, 采用分裂思想将NLS方程分解为线性导数项和非线性项两个微分方程; 其次, 基于Taylor展开拓展应用文献[17, 21]中修正SPH方法对含线性导数方程进行二阶精度显式求解; 再次, 运用四阶对称分裂格式对上述两个方程进行积分离散; 最后引入基于粒子搜索的MPI并行技术提高计算效率. 给出的HSS-CPSPH方法将兼具传统SPH方法和已有改进SPH以及分裂格式的优点, 也较已有分裂有限差分法^[6-8]具有更好、更灵活的模拟应用性. 在数值算例模拟中, 将HSS-CPSPH方法的数值结果与解析解或其他数值结果进行了比较, 分析了数值精度、收敛阶和计算效率. 比较分析结果表明, 给出的HSS-CPSPH方法可以准确、高效地求解周期和Dirichlet边界条件下二维或三维NLS/GP方程, 并能成功地预测二分量孤立波传播中的非线性奇异特性和带外旋转算子下量子涡旋的复杂变化过程.

量子力学中许多非线性物理现象常用Gross-Pitaevskii (GP)方程^[3,4,6,7]来描述, 其高维直角坐标系下无量纲化的控制方程为

$\begin{split}{\rm{i}}\dfrac{{\partial \phi \left( {{{x}},t} \right)}}{{\partial t}} = &- \dfrac{{\rm{1}}}{{\rm{2}}}\Delta \phi + {V_d}\left( {{{x}},t} \right)\phi + {\beta _d}{\left| \phi \right|^2}\phi - {\varOmega}{L_z}\phi ,\\ &{{x}} \in {{{R}}^d},\;t \geqslant 0,\end{split}$

其中${\rm{i}} = \sqrt { - 1} $为虚数单位; $\Delta $是拉普拉斯算子; $\phi $为待求的$d$维复数域上的函数($d = 1,2,3$); ${{x}} = x$或$\left( {x,y} \right)$或$\left( {x,y,z} \right)$; 无量纲实常数${\beta _d}$常常用来描述三次非线性项${\left| \phi \right|^2}\phi $的相互作用强度; ${V_d}$是势函数,

${V_d}\left( {x,t} \right) = \left\{ \begin{aligned}& \dfrac{1}{2}{x^2},\quad\quad\quad\quad\quad\quad\quad\quad\quad\; d = 1,\\& \left( {{x^2} + {\gamma _y}^2{y^2}} \right)/2,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;d = 2,\\& \left( {{x^2} + {\gamma _y}^2{y^2} + {\gamma _z}^2{z^2}} \right)/2,\;\;d = 3,\end{aligned} \right.$

${\gamma _y},{\gamma _z}$是实常数; ${L_z} = - {\rm{i}}\left( {x\partial y - y\partial x} \right)$是玻色-爱因斯坦凝聚态(Bose-Einstein condensates, BEC)下无量纲速度旋转项$\varOmega $的角动量算符的$z$分量.
为封闭上述方程, 考虑如下初始条件:

$\phi \left( {{{x}},0} \right) = \phi \left( {{x}} \right),\;{{x}} \in {{{R}}^d},$

边界条件

$\mathop {\lim }\limits_{\left| {{x}} \right| \leftarrow \infty } \phi \left( {{{x}},t} \right) = 0,\;t \geqslant 0,$

或周期边界条件$\phi \left( {x,t} \right) = \phi \left( {x + l,t} \right)$, 其中$l$表示方程周期, 或$\phi \left( {x,y,t} \right) = \phi \left( {x + {l_1},y,t} \right)$, $\phi \left( {x,y,t} \right) = $$\phi \left( {x,y + {l_2},t} \right)$, 其中$\left( {{l_1},{l_2}} \right)$为$\left( {x,y} \right)$方向的周期.
根据文献[6—8]可知, 上述无穷区域边界条件(4)式可以近似处理为齐次Dirichlet边界, 周期边界条件采用两边施加几层虚粒子方式进行处理(详见3.3节).

针对高维非线性GP方程的粒子方法求解, 拓展应用文献[17, 21]中给出的一阶修正SPH方法时, 随着模拟时间的延长会出现精度低和稳定性差的现象, 这是因为本文考虑的GP方程中存在非线性项和一阶导数项, 对模拟算法的精度和稳定性要求较高. 为此, 在时间高阶分裂法^[1,2]和已有改进SPH方法^[15-18]的基础上, 引入MPI并行算法, 提高高维问题模拟计算效率. 本文对二维/三维GP方程模拟给出一种HSS-CPSPH方法, 基本思想是: 首先, 采用时间分裂思想将GP方程分成线性导数和非线性项两个方程; 其次, 对线性导数方程拓展应用文献[15, 17]中的改进SPH方法进行离散; 再次, 引入文献[1, 2]中四阶分裂格式对分裂的两个方程进行交替求解; 最后, 采用基于MPI的并行算法来提高计算机模拟效率.
2

-->

3.1.时间分裂法^[1-3,6]

作为一种精确高效的算法, 时间分裂法在复数域非线性薛定谔方程求解中已被广泛应用, 根据文献[6—8], GP方程(1)可以重写为

${\rm{i}}\dfrac{{\partial \phi \left( {{{x},}t} \right)}}{{\partial t}} = \left( {A + B} \right)\phi \left( {{{x}},t} \right),\;{{x}} \in {{{R}}^d},$

其中线性微分算子$A = - \dfrac{1}{2}\Delta - \varOmega {L_z}$, 非线性算子$B = {V_d}\left( {{{x}},t} \right) + {\beta _d}\left| {\phi \left( {{{x}},t} \right)} \right|{^2}$. 由文献[26]知 (5)式在新时刻里的近似解析形式可以写成$\phi ({{x}},t + \tau) \approx$$ \exp \left[ {{\rm{i}}\tau \left( {A + B} \right)} \right]\phi \left( {{{x}},t} \right)$, 即: $\phi \left( {{{x}},t + \tau } \right) \approx \exp \left( {{\rm{i}}\tau A} \right) \cdot $$\left( {\exp \left( {{\rm{i}}\tau B} \right) \cdot \phi \left( {{{x}},t} \right)} \right)$. 根据该近似表达式和文献[26]中分解理论可知, (5)式可以分解成下面的线性问题

${\rm{i}}\dfrac{{\partial \phi \left( {{{x}},t} \right)}}{{\partial t}} = A\phi $

和非线性问题

${\rm{i}}\dfrac{{\partial \phi \left( {{{x}},t} \right)}}{{\partial t}} = B\phi .$

上述过程的成立条件和精度可见文献[26].
本文对方程(6)和(7)的求解采用如下四阶分裂法^[1,2]: 第一步求解方程(7), 第二步用第一步得到的解作为初始条件来求解方程(6), 第三步用第二步得到的解作为初始条件求解方程(7), 第四步使用第三步得到的解作为初始条件求解方程(6), 第五步用第四步得到的解作为初始条件求解方程(7), 此时得到的解作为下一时间层的数值解.
2 -->

3.2.高阶分裂修正SPH离散格式

3 -->

3.2.1.修正SPH离散格式

在传统SPH方法模拟中, 首先使用核函数进行积分插值, 其次在模拟区域$\varOmega $内, 对有限个粒子赋予密度、质量、温度等物理含义, 最后用粒子离散化方式近似得到核函数的积分形式.
在任意点${{x}} = \left( {x,y,z} \right)$处, 函数$\phi \left( {x} \right)$及其微分形式可以用核函数$W$和核函数的导数$\nabla W$表示^[19]:

$\left\langle {\phi \left( {{x}} \right)} \right\rangle = \int_\varOmega {\phi \left( {{{{x}}{'}}} \right)} W\left( {{{x}} - {{{x}}{'}},h} \right){\rm{d}}{{{x}}{'}},$

$\left\langle {\nabla \phi \left( x \right)} \right\rangle = \int_\varOmega {\phi \left( {{{{x}}{'}}} \right)\nabla W{\rm{d}}{{{x}}{'}}} - \phi \left( {{x}} \right)\int_\varOmega {\nabla W{\rm{d}}{{{x}}{'}}}, $

其中$h$为光滑长度, 决定了核函数$W$的支持域范围, $W$一般要求满足正则化、对称性、紧致性和Dirac函数性质^[19]. 本文模拟中取$h = 0.9{d_0} - 1.1{d_0}$(${d_0}$为初始粒子间距), 取分段五次样条核函数^[19], 对应支持域尺寸为${\rm{3}}h$.
引入体积$v = m/\rho $ (其中$m,\;\rho $都是实常数), 对于任意位置${{{x}}_i} = \left( {{x_i},\;{y_i},\;{z_i}} \right)$处粒子$i$, 其支持域内的相邻粒子为$j$, 则有粒子近似公式:

${\phi _i} = \sum\limits_j {{\phi _j}{W_{ij}}{v_j}} ,$

$\left( {\dfrac{{\partial \phi }}{{\partial {x_i}}}} \right) = \sum\limits_j {\left( {{\phi _j} - {\phi _i}} \right){\nabla _i}{W_{ij}}} {v_j}.$

${W_{ij}} \!= W\left( {\left| {{{{x}}_i} - {{{x}}_j}} \right|,h} \right)$, ${\nabla _i}{W_{ij}} \!= \partial W(\left| {{{{x}}_i} - {{x}}_j} \right|,h)/$$\partial {{{x}}_i}$, 且${\nabla _i}{W_{ij}}$满足关系

${\nabla _i}{W_{ij}} = - {\nabla _j}{W_{ij}}.$

为了提高传统SPH方法的精度和稳定性, 基于Taylor展开, 人们已经提出了一些改进SPH方法^[15-22]. 本文拓展应用文献[17, 21]中具有二阶精度的一阶修正对称SPH方法, 该方法首先将二阶导数项分为两个一阶导数, 然后对一阶导数基于Taylor展开进行离散, 则(11)式变为

$\left( {\dfrac{{\partial \phi }}{{\partial {{{x}}_i}}}} \right) = \sum\limits_j {\left( {{\phi _j} - {\phi _i}} \right)} \nabla _i^c {{\widehat W_{ij}}} {v_j},$

其中$\nabla _i^c {{\widehat W_{ij}}} {v_j}$通过求解下列方程组得到

${A^s}\left(\!\!\!\begin{array}{l}\partial \widehat W {_{ij}^c}/\partial {x_i}\\\partial \widehat W {_{ij}^c}/\partial {y_i}\\\partial \widehat W {_{ij}^c}/\partial {z_i}\end{array}\!\!\!\right) = \left(\!\!\begin{array}{l}{x_{ji}}{W_{ij}}\\{y_{ji}}{W_{ij}}\\{z_{ji}}{W_{ij}}\end{array}\!\!\right).$

其中${A^s} =\displaystyle \sum\limits_{j = 1}^N {\left( {{W_{ij}}\left( \begin{array}{l} {x_{ji}}{x_{ji}}\;{y_{ji}}{x_{ji}}\;{z_{ji}}{x_{ji}} \\ {x_{ji}}{y_{ji}}\;{y_{ji}}{y_{ji}}\;{z_{ji}}{y_{ji}} \\ {x_{ji}}{z_{ji}}\;{y_{ji}}{z_{ji}}\;{z_{ji}}{z_{ji}}\end{array} \right){v_j}} \right)} $, N是粒子$i$支持域内相邻粒子个数, ${x_{ji}} = {x_j} - {x_i}$, ${y_{ji}} = {y_j} - {y_i}$, ${z_{ji}} = {z_j} - {z_i}$.
3 -->

3.2.2.高阶分裂修正SPH格式

将 3.1节四阶分裂格式与 3.2.1节修正SPH格式进行结合, 对方程(1)进行离散, 可得如下高阶分裂修正SPH离散格式

$\phi _i^* = {{\rm{e}}^{ - {\rm{i}}{d_1}\left( {{V_1} + \beta {{\left| {\phi _i^n} \right|}^2}} \right){\rm{d}}t}}\phi _i^n,$

$\left( {\dfrac{{\partial \phi }}{{\partial {{x}}}}} \right)_i^{**} = {\sum\limits_j {{v_j}\phi _{ji}^{**}\left( {\dfrac{{\partial {\widehat W_{ij}^c} }}{{\partial {{x}}}}} \right)} _i},$

$\begin{split} &{\rm{i}}\dfrac{{\phi _i^{**} - \phi _i^*}}{{{\rm{d}}t}} \\=\; & {c_1}\Bigg[ { - \dfrac{1}{2}} \Bigg.\left\{ {\sum\limits_j {{v_j}\left( {\left( {\dfrac{{\partial \phi }}{{\partial x}}} \right)_j^{**}\!\!- \left( {\dfrac{{\partial \phi }}{{\partial x}}} \right)_i^{**}} \right)} }\right.\!\!{\left( {\dfrac{{\partial {\widehat W_{ij}^c} }}{{\partial x}}} \right)_i}\\ &+ \sum\limits_j {{v_j}\left( {\left( {\dfrac{{\partial \phi }}{{\partial y}}} \right)_j^{**} - \left( {\dfrac{{\partial \phi }}{{\partial y}}} \right)_i^{**}} \right)} {\left( {\dfrac{{\partial {\widehat W_{ij}^c} }}{{\partial y}}} \right)_i} \\&+ \left. {\sum\limits_j {{v_j}\left( {\left( {\dfrac{{\partial \phi }}{{\partial z}}} \right)_j^{**} - \left( {\dfrac{{\partial \phi }}{{\partial z}}} \right)_i^{**}} \right)} {{\left( {\dfrac{{\partial {\widehat W_{ij}^c} }}{{\partial z}}} \right)}_i}} \right\}\\&+ \Bigg. {{\rm{i}}\varOmega \left( {{x_i}\left( {\dfrac{{\partial \phi }}{{\partial y}}} \right)_i^{**} - {y_i}\left( {\dfrac{{\partial \phi }}{{\partial x}}} \right)_i^{**}} \right)} \Bigg],\end{split} $

$\left( {\dfrac{{\partial \phi }}{{\partial {{x}}}}} \right)_i^{**{\rm{*}}} = {\sum\limits_j {{v_j}\phi _{ji}^{**{\rm{*}}}\left( {\dfrac{{\partial {\widehat W_{ij}^c}}}{{\partial {{x}}}}} \right)} _i},$

$\begin{split} & {\rm{i}}\dfrac{{\phi _i^{**{\rm{*}}} - \phi _i^{*{\rm{*}}}}}{{{\rm{d}}t}} \\ =\; & {c_2}\Bigg[ { - \dfrac{1}{2}} \Bigg.\left\{\!{\sum\limits_j {{v_j}\left( {\left( {\dfrac{{\partial \phi }}{{\partial x}}} \right)_j^{**{\rm{*}}}\!\!- \left( {\dfrac{{\partial \phi }}{{\partial x}}} \right)_i^{**{\rm{*}}}} \right)} }\right.\!\!{\left( {\dfrac{{\partial {\widehat W_{ij}^c} }}{{\partial x}}} \right)_i} \\ & + \sum\limits_j {{v_j}\left( {\left( {\dfrac{{\partial \phi }}{{\partial y}}} \right)_j^{**{\rm{*}}} - \left( {\dfrac{{\partial \phi }}{{\partial y}}} \right)_i^{**{\rm{*}}}} \right)} {\left( {\dfrac{{\partial {\widehat W_{ij}^c} }}{{\partial y}}} \right)_i} \\& + \left. {\sum\limits_j {{v_j}\left( {\left( {\dfrac{{\partial \phi }}{{\partial z}}} \right)_j^{*{\rm{*}}*} - \left( {\dfrac{{\partial \phi }}{{\partial z}}} \right)_i^{**{\rm{*}}}} \right)} {{\left( {\dfrac{{\partial {\widehat W_{ij}^c} }}{{\partial z}}} \right)}_i}} \right\}\\ & + \Bigg. {{\rm{i}}\varOmega \left( {{x_i}\left( {\dfrac{{\partial \phi }}{{\partial y}}} \right)_i^{**{\rm{*}}} - {y_i}\left( {\dfrac{{\partial \phi }}{{\partial x}}} \right)_i^{**{\rm{*}}}} \right)} \Bigg],\end{split}$

$\phi _i^{*{\rm{***}}} = {{\rm{e}}^{ - {\rm{i}}{d_2}\left( {{V_1} + \beta {{\left| {\phi _i^{{\rm{***}}}} \right|}^2}} \right){\rm{d}}t}}\phi _i^{{\rm{***}}},$

$\left( {\dfrac{{\partial \phi }}{{\partial {{x}}}}} \right)_i^{**{\rm{***}}} = {\sum\limits_j {{v_j}\phi _{ji}^{**{\rm{***}}}\left( {\dfrac{{\partial {\widehat W_{ij}^c} }}{{\partial {{x}}}}} \right)} _i},$

$\begin{split} & {\rm{i}}\dfrac{{\phi _i^{**{\rm{***}}} - \phi _i^{*{\rm{***}}}}}{{{\rm{d}}t}}\\ =\;& {c_3}\Bigg[ { - \dfrac{1}{2}}\!\! \Bigg.\left\{\!{\sum\limits_j {{v_j}\!\!\left( {\left( {\dfrac{{\partial \phi }}{{\partial x}}} \right)_j^{**{\rm{***}}} \!\!-\!\!\left( {\dfrac{{\partial \phi }}{{\partial x}}} \right)_i^{**{\rm{***}}}} \right)}}\right.\!\!{\left( {\dfrac{{\partial {\widehat W_{ij}^c} }}{{\partial x}}} \right)_i}\\&+ \sum\limits_j {{v_j}\left( {\left( {\dfrac{{\partial \phi }}{{\partial y}}} \right)_j^{**{\rm{***}}} - \left( {\dfrac{{\partial \phi }}{{\partial y}}} \right)_i^{**{\rm{***}}}} \right)} {\left( {\dfrac{{\partial {\widehat W_{ij}^c} }}{{\partial y}}} \right)_i}\\ & + \left. {\sum\limits_j {{v_j}\left( {\left( {\dfrac{{\partial \phi }}{{\partial z}}} \right)_j^{*{\rm{*}}*{\rm{**}}} - \left( {\dfrac{{\partial \phi }}{{\partial z}}} \right)_i^{**{\rm{***}}}} \right)} {{\left( {\dfrac{{\partial {\widehat W_{ij}^c} }}{{\partial z}}} \right)}_i}} \right\}\\&+ \Bigg. {{\rm{i}}\varOmega \left( {{x_i}\left( {\dfrac{{\partial \phi }}{{\partial y}}} \right)_i^{**{\rm{***}}} - {y_i}\left( {\dfrac{{\partial \phi }}{{\partial x}}} \right)_i^{**{\rm{***}}}} \right)} \Bigg] ,\end{split}$

$\left( {\dfrac{{\partial \phi }}{{\partial {{x}}}}} \right)_i^{**{\rm{****}}} = {\sum\limits_j {{v_j}\phi _{ji}^{**{\rm{****}}}\left( {\dfrac{{\partial {\widehat W_{ij}^c}}}{{\partial {{x}}}}} \right)} _i},$

$\begin{split} & {\rm{i}}\dfrac{{\phi _i^{**{\rm{****}}} - \phi _i^{*{\rm{****}}}}}{{{\rm{d}}t}} \\ =\; & {c_4}\Bigg[ { - \dfrac{1}{2}} \!\!\Bigg.\left\{\!{\sum\limits_j {{v_j}\!\!\left( {\left( {\dfrac{{\partial \phi }}{{\partial x}}} \right)_j^{**{\rm{****}}} \!-\!\!\left( {\dfrac{{\partial \phi }}{{\partial x}}} \right)_i^{**{\rm{****}}}}\right)}} \right.\!\!{\left( {\dfrac{{\partial {\widehat W_{ij}^c} }}{{\partial x}}} \right)_i} \\ & + \sum\limits_j {{v_j}\left( {\left( {\dfrac{{\partial \phi }}{{\partial y}}} \right)_j^{**{\rm{****}}} - \left( {\dfrac{{\partial \phi }}{{\partial y}}} \right)_i^{**{\rm{****}}}} \right)} {\left( {\dfrac{{\partial {\widehat W_{ij}^c} }}{{\partial y}}} \right)_i} \\& + \left. {\sum\limits_j {{v_j}\left( {\left( {\dfrac{{\partial \phi }}{{\partial z}}} \right)_j^{*{\rm{*}}*{\rm{***}}} - \left( {\dfrac{{\partial \phi }}{{\partial z}}} \right)_i^{**{\rm{****}}}} \right)} {{\left( {\dfrac{{\partial {\widehat W_{ij}^c} }}{{\partial z}}} \right)}_i}} \right\} \\ & +\Bigg. {{\rm{i}}\varOmega \left( {{x_i}\left( {\dfrac{{\partial \phi }}{{\partial y}}} \right)_i^{**{\rm{****}}} - {y_i}\left( {\dfrac{{\partial \phi }}{{\partial x}}} \right)_i^{**{\rm{****}}}} \right)} \Bigg] ,\end{split}$

$\phi _i^{{n} + {\rm{1}}} = {{\rm{e}}^{ - {\rm{i}}{d_3}({V_1} + \beta |\phi _i^{{\rm{******}}}{|^2}){\rm{d}}t}}\phi _i^{{\rm{******}}},$

其中$i = 0,\;1,\;2,\; \cdots,M;\;n = 0,\;1,\; \cdots;$$M$是区间内所有粒子总数, $n$是时间层, ${\rm{d}}t$是时间步长,
$\phi _{ji}^{**} = \phi _j^{**} - \phi _i^{**}$, $\phi _{ji}^{**{\rm{*}}} = \phi _j^{**{\rm{*}}} - \phi _i^{**{\rm{*}}}$,
$\phi _{ji}^{**{\rm{***}}} \!=\! \phi _j^{**{\rm{***}}} \!-\! \phi _i^{**{\rm{***}}},$ $\phi _{ji}^{**{\rm{****}}} \!=\! \phi _j^{**{\rm{****}}} \!-\! \phi _i^{**{\rm{****}}},$
四阶分裂格式参数
${c_1} = {c_4} = \dfrac{1}{{2 \times \left( {2 - {2^{1/3}}} \right)}}$, ${c_2} = {c_3} = \dfrac{{1 - {2^{1/3}}}}{{2 \times \left( {2 - {2^{1/3}}} \right)}}$, ${d_1} = {d_3} = \dfrac{1}{{2 - {2^{1/3}}}}$, ${d_2} = - \dfrac{{{2^{1/3}}}}{{2 - {2^{1/3}}}}$^[2].
上述离散格式被称为高阶分裂修正SPH方法, 结合3.4节MPI并行计算技术将本文提出的方法称为“HSS-CPSPH”, 其中对线性导数方程(6)的时间离散采用二阶龙格-库塔显格式, 根据文献[16—19]中的稳定性限制条件, 时间步长取$ {\rm{d}}t \leqslant$$ 0.1d_0^2$. 本文方法与文献[15]中给出的SS-ICPSPH方法主要的不同之处在于: 对线性导数方程的离散本文采用二阶龙格-库塔显格式, 后者采用隐格式; 本文采用四阶分裂格式, 后者采用二阶分裂格式, 使得本文方法精度较文献[15]中的方法精度高(见第4节); 文献[15]运用提出的方法仅对带Dirichlet边界的GP方程进行模拟研究, 本文运用给出的方法对带周期边界和Dirichlet边界的非线性薛定谔问题都进行了模拟分析, 并且准确预测了周期边界下孤立波传播中出现的奇异现象(见5.1节).
2 -->

3.3.边界条件处理

运用上述提出的方法对高维GP方程进行模拟时, 初边值条件的准确处理对数值模拟的精度和稳定性至关重要. 对初始条件可以准确离散施加, 边界条件(4)式可以近似处理为Dirichlet边界. 以二维区域情况为例, 周期性边界条件处理如下:

$\left( {x,y} \right) \in \left[ {{X_a},{X_b}} \right] \times \left[ {{Y_a},{Y_b}} \right],\;0 \leqslant t \leqslant T,$

其中${X_a},\;{X_b},\;{Y_a},\;{Y_b}$分别为$x,\;y$方向区域的最小和最大值. 将区域沿$\left( {x,\;y} \right)$方向分别均分成$J$和$K$整数份, 令空间步长${h_1} = \left( {{X_b} - {X_a}} \right)/J$, ${h_2} = \left( {{Y_b} - {Y_a}} \right)/K$, ${x_j} = {X_a} + j{h_1}$, $0 \leqslant j \leqslant J - 1$, ${y_k} = {Y_a} + k{h_2}$, $0 \leqslant k $$\leqslant K - 1$. 为防止粒子方法处理周期边界的粒子缺失问题, 在区间外侧各扩充4个点(大于等于$3h$), 可以得到

$\begin{split}& {x_{ - 1}} = {X_a} - {h_1},\;{x_{ - 2}} = {X_a} - 2{h_1},\\&{x_{ - 3}} = {X_a} - 3{h_1},\;{x_{ - 4}} = {X_a} - 4{h_1},\\& {x_{J + 1}} = {X_a} + \left( {J + 1} \right){h_1},\;{x_{J + 2}} = {X_a} + \left( {J + 2} \right){h_1},\\&{x_{J + 3}} = {X_a} + \left( {J + 3} \right){h_1},\;{x_{J + 4}} = {X_a} + \left( {J + 4} \right){h_1},\\ & {y_{ - 1}} = {Y_a} - {h_2},\;{y_{ - 2}} = {Y_a} - 2{h_2},\\&{y_{ - 3}} = {Y_a} - 3{h_2},\;{y_{ - 4}} = {Y_a} - 4{h_2},\\& {y_{J + 1}} = {Y_a} + \left( {J + 1} \right){h_2},\;{y_{J + 2}} = {Y_a} + \left( {J + 2} \right){h_2},\\&{y_{J + 3}} = {Y_a} + \left( {J + 3} \right){h_2},\;{y_{J + 4}} = {Y_a} + \left( {J + 4} \right){h_2},\end{split}$

则具有$\left( {{l_1},\;{l_2}} \right)$周期边界条件的方程满足如下条件:
${\phi _{ - 4,k}} = {\phi _{J - 4,k}}$, ${\phi _{ - 3,k}} = {\phi _{J - 3,k}}$, ${\phi _{ - 2,k}} = {\phi _{J - 2,k}}$, ${\phi _{{\rm{ - }}1,k}} = {\phi _{J - 1,k}}$,
${\phi _{J,k}} = {\phi _{0,k}}$, ${\phi _{J + 1,k}} = {\phi _{1,k}}$, ${\phi _{J + 2,k}} = {\phi _{2,k}}$, ${\phi _{J + 3,k}} = {\phi _{3,k}}$, ${\phi _{J + 4,k}} = {\phi _{4,k}}$, $k = - 4, - 3, \cdots,K + 4$;
${\phi _{j, - 4}} = {\phi _{j,K - 4}}$, ${\phi _{j, - 3}} = {\phi _{j,K - 3}}$, ${\phi _{j, - 2}} = {\phi _{j,K - 2}}$, ${\phi _{j, - 1}} = {\phi _{j,K - 1}}$,
${\phi _{j,K}} = {\phi _{j,0}}$, ${\phi _{j,K + 1}} = {\phi _{j,1}}$, ${\phi _{j,K + 2}} = {\phi _{j,2}}$, ${\phi _{j,K + 3}} = {\phi _{j,3}}$, ${\phi _{j,K + 4}} = {\phi _{j,4}}$, $j = - 4, - 3, \cdots,J + 4$.
本文在第4和第5节通过算例模拟验证了上述边界处理的准确性和有效性.
2 -->

3.4.粒子搜索并行

SPH粒子方法模拟中, 相邻粒子搜索和高维问题下需要几十万至上百万粒子物理量更新计算, 会使得计算机模拟中占较大计算内存和较长CPU计算时间, 因此运用提出方法对高维GP进行模拟时提高计算效率是必要的. 根据模拟中粒子位置不变的特点, 引入文献[17, 21]中基于MPI粒子搜索并行技术. 即基于MPI并行算法, 首先对每个粒子支持域内的相邻粒子进行标记, 标记后对所有粒子进行物理量赋值, 主要是质量和密度, 为了提高计算效率, 在并行计算中, 不同节点间会进行密度通信. 在计算过程中节点间进行空间上的一阶导数项通信以及最后计算函数值的更新. 本文的并行计算是将所有粒子进行编号, 再根据CPU数量将粒子平均分配, 进行同时计算. 每个粒子与之相互作用的相邻粒子支持域范围通过下列核函数确定,

$\begin{split}& {W_{ij}} = W(r_{ij},h) \\= &\; {w_0} \begin{cases} (3 - q)^5 - 6 (2 - q)^5 + 15(1 - q)^5, & 0 \leqslant q < 1, \\ (3 - q)^5 - 6(2 - q)^5 , & 1 \leqslant q < 2, \\ (3 - q)^5, & 2\leqslant q < 3, \\ 0, & q \geqslant 3,\end{cases} \end{split}$

其中$q = {r_{ij}}/h$, ${r_{ij}} = \left| {{{{x}}_i} - {{{x}}_j}} \right|$, 系数${w_0}$为$120/h$ (一维), $7/(478{\text{π}}{h^2})$ (二维), $3/(359{\text{π}}{h^3})$ (三维). 本文中的支持域范围是以$3h$为半径的圆或球.

本节通过对具有解析解带不同边界条件的非线性薛定谔方程进行模拟, 对提出的HSS-CPSPH方法的精度、收敛阶以及计算效率进行分析. 分别对二维带周期边值、一维二分量带周期边值和二维带第一类边值的问题进行数值模拟, 并与解析解做比较, 分析所提方法的准确性和高效性, 与已有SS-ICPSPH方法^[10]做对比证明本文所提方法具有较高精度和较快的收敛速度.
为分析方法的数值精度和收敛速度, 定义如下最大误差范数和收敛阶:

${e_{\rm{m}}} = {\left\| {{\phi _{{\rm{exact}}}} - {\phi _{{\rm{numerical}}}}} \right\|_\infty },$

$o{r_{\rm{\alpha }}} \approx \dfrac{{\log \left( {{e_{\rm{m}}}\left( {{d_{02}}} \right)/{e_{\rm{m}}}\left( {{d_{01}}} \right)} \right)}}{{\log \left( {{d_{02}}/{d_{01}}} \right)}},$

其中${d_{01}}$和${d_{02}}$代表了两种不同的粒子初始间距.
2

-->

4.1.二维带周期边界下孤立波传播

考虑方形区域$\left[ {0,\;2{\text{π}}} \right] \times \left[ {0,\;2{\text{π}}} \right]$上带周期边界的非线性薛定谔方程^[3]

${\rm{i}}{u_t} + \Delta u + \beta {\left| u \right|^2}u = 0,$

对应解析解为

$u\left( {x,y,t} \right) = A\exp \left[ {{\rm{i}}\left( {{k_1}x + {k_2}y - wt} \right)} \right],$

其中$w = k_1^2 + k_2^2 - \beta {\left| A \right|^2}$, 系数$A$ = 1, ${k_1} = {k_2}$为常数, $\beta = - 2$. 数值模拟中初始条件可以从解析解中获得, 对应的周期长度为$\left( {2{\text{π}},\;2{\text{π}}} \right)$且

$\begin{split}& u(x,y,t) = u(x + 2{\rm{\pi }},y,t),\\& u(x,y,t) = u(x,y + 2{\rm{\pi }},t),\\& (x,y) \in R \times R,\;0 < t \leqslant T,\;T = 1. \end{split}$

该算例模拟中, ${k_1} = {k_2} = 1$时采用$129 \times 129$个均匀分布粒子, 时间步长为${\rm{d}}t = {10^{ - 4}}$(见图1); ${k_1} = {k_2} = 4$时采用$h = {\text{π}} /32$个均匀分布粒子, 时间步长为${\rm{d}}t = {10^{ - 4}}$(见图2). 图1和图2分别展示了${k_{\rm{1}}}$和${k_{\rm{2}}}$取不同系数时由HSS-CPSPH方法获得的$u\left( {x,{\text{π}}} \right)$的实部曲线图, 并将其与解析解及SS-ICPSPH结果进行对比, 可以看出, HSS-CPSPH结果与解析解相符合, 但随着时间的延长, 两种数值结果产生较小的偏差, 这与文献[3]中分裂有限差分法得到数值模拟结论类似.
图 1 ${k_{\rm{1}}} = {k_{\rm{2}}} = 1,\;h = {\text{π}}/64$时不同时刻$u\left( {x,{\text{π}}} \right)$的实部沿x轴的变化　(a) $t=1$; (b) $t = 3$
Figure1. Curve of the $\operatorname{Re} \left( {u\left( {x,{\text{π}}} \right)} \right)$ along x-axis at different time with ${k_{\rm{1}}} = {k_{\rm{2}}} = 1,\;h = {\text{π}}/64$: (a) $t=1$; (b) $t = 3$.

图 2 ${k_{\rm{1}}} = {k_{\rm{2}}} = 4,\;h = {\text{π}}/128$时不同时刻$u\left( {x,{\text{π}}} \right)$的实部沿x轴的变化　(a) t = 0.1; (b) t = 1
Figure2. Curve of the $\operatorname{Re} \left( {u\left( {x,{\text{π}}} \right)} \right)$ along x-axis at different time with ${k_{\rm{1}}} = {k_{\rm{2}}} = 4,\;h = {\text{π}}/128$: (a) t = 0.1; (b) t = 1.

为体现提出的HSS-CPSPH方法求解周期性问题的精度和收敛速度, 表1和表2分别列出了模拟较短时间内数值结果的误差和收敛阶, 由表1和表2可以看出: 1) ${e_{\rm m}}$随着粒子数的增加而减小, HSS-CPSPH方法得到的误差较SS-ICPSPH方法的小; 2) HSS-CPSPH方法较SS-ICPSPH方法具有更快的收敛速度. 为进一步体现粒子方法在粒子分布非均匀情况下数值模拟的精度, 表3列出了粒子分布均匀和两种非均匀情况下(见文献[15])两个不同时刻的最大误差${e_{\rm{m}}}$. 由表3可知, 粒子方法在粒子分布均匀和分布非均匀情况下得到的数值结果误差都比较接近, 表明了该方法易推广应用到非规则区域问题的模拟, 且保持较好的精度.

时间t	SS-ICPSPH	HSS-CPSPH
0.5	1.697 × 10–3	1.696 × 10–3
1	3.616 × 10–3	2.494 × 10–3
2	7.347 × 10–3	4.857 × 10–3

表1${k_{\rm{1}}} = {k_{\rm{2}}} = 1,h = {\text{π}}/64$时几个不同时刻里两种方法的误差${e_{\rm{m}}}$
Table1.Error ${e_{\rm{m}}}$ obtained using two different methods at different time (${k_{\rm{1}}} = {k_{\rm{2}}} = 1,h = {\text{π}}/64$).

	$h = {\text{π}}/32$			$h = {\text{π}}/64$			$h = {\text{π}}/128$
	${e_{\rm{m}}}$	$o{r_{\rm{\alpha }}}$		${e_{\rm{m}}}$	$o{r_{\rm{\alpha }}}$		${e_{\rm{m}}}$	$o{r_{\rm{\alpha }}}$
SS-ICPSPH	1.381 × 10–2	—		3.616 × 10–3	1.933		9.0412 × 10–4	2.00
HSS-CPSPH	1.381 × 10–2	—		2.494 × 10–3	2.47		4.498 × 10–4	2.47

表2${k_{\rm{1}}} = {k_{\rm{2}}} = 1$, 时间t = 1时, 两种方法在不同粒子间距下的误差和收敛阶
Table2.Error ${e_{\rm{m}}}$ and convergent order $o{r_{\rm{\alpha }}}$ obtained using two different methods at $t=1$ and different particle distance (${k_1} = {k_2} = 1$).

	均匀分布粒子			非均匀分布情形1			非均匀分布情形2
	$t = 0.1$	$t = 1$		$t = 0.1$	$t = 1$		$t = 0.1$	$t = 1$
SS-ICPSPH	2.776 × 10–4	3.616 × 10–3		2.944 × 10–4	3.818 × 10–3		3.116 × 10–4	4.082 × 10–3
HSS-CPSPH	2.774 × 10–4	2.494 × 10–3		2.886 × 10–4	2.527 × 10–3		2.967 × 10–4	2.578 × 10–3

表3${k_{\rm{1}}} = {k_{\rm{2}}} = 1,h = {\text{π}}/64$时, 粒子分布均匀或不均匀方式下, 两种方法的误差${e_{\rm{m}}}$
Table3.Error ${e_{\rm{m}}}$ obtained using different methods at different distribution (${k_1} = {k_2} = 1$,$h = {\text{π}}/64$).

2 -->

4.2.一维二分量带周期边界下孤立波传播

为进一步展示本文提出的HSS-CPSPH方法能够准确捕捉多分量多孤立波传播中的尖角现象, 本小节考虑区间$[ - 20,\;80]$上一维二分量带周期边界两种初始条件下的非线性薛定谔方程^[14]

$\left\{ \begin{aligned}& {\rm{i}}{u_t} + {\rm{i}}\alpha {u_x} + \dfrac{1}{2}{u_{xx}} + \left( {{{\left| u \right|}^2} + \beta {{\left| v \right|}^2}} \right)u = 0,\\& {\rm{i}}{v_t} - {\rm{i}}\alpha {v_x} + \dfrac{1}{2}{v_{xx}} + \left( {{{\left| v \right|}^2} + \beta {{\left| u \right|}^2}} \right)v = 0.\end{aligned} \right.$

初始条件1为

$\begin{split} u\left( {x,0} \right) =\;&\sqrt {\dfrac{{2\alpha }}{{1 + \beta }}} {\rm{sech}}\left( {\sqrt {2a} \left( {x - {x_0}} \right)} \right)\\ & \times\exp \left\{ {{\rm{i}}\left[ {\left( {c - \alpha } \right)\left( {x - {x_0}} \right)} \right]} \right\},\\ v\left( {x,0} \right) =\;&\sqrt {\dfrac{{2\alpha }}{{1 + \beta }}} {\rm{sech}}\left( {\sqrt {2a} \left( {x - {x_0}} \right)} \right)\\ &\times\exp \left\{ {{\rm{i}}\left[ {\left( {c - \alpha } \right)\left( {x - {x_0}} \right)} \right]} \right\},\end{split}$

其中参数$\alpha = 0.5,\;\beta = 2/3,\;c = 1,\;a = 1,\;{x_0} = 0$.
具有三个孤子波传播的初始条件2为

$\begin{aligned} u\left( {x,0} \right) =\; & \sum\limits_{j = 1}^3 {\sqrt {\dfrac{{2{a_j}}}{{1 + \beta }}} } {\rm{sech}}\left( {\sqrt {2{a_j}} \left( {x - {x_j}} \right)} \right)\\ & \times\exp \left\{ {{\rm{i}}\left[ {\left( {{c_j} - \alpha } \right)\left( {x - {x_j}} \right)} \right]} \right\},\\ v\left( {x,0} \right) =\; &\sum\limits_{j = 1}^3 {\sqrt {\dfrac{{2{a_j}}}{{1 + \beta }}} } {\rm{sech}}\left( {\sqrt {2{a_j}} \left( {x - {x_j}} \right)} \right)\\ & \times\exp \left\{ {{\rm{i}}\left[ {\left( {{c_j} - \alpha } \right)\left( {x - {x_j}} \right)} \right]} \right\},\end{aligned}$

其中$\alpha = 0.5$, $\beta = 2/3$, ${c_1} = 1$, ${c_2} = 0.1$, ${c_3} = - 1$, ${a_1} = 1$, ${a_2} = 0.72$, ${a_3} = 0.36$, ${x_1} = 0$, ${x_2} = 25$, ${x_3} = 50$.
周期边界条件$u\left( {x,t} \right) = u\left( {x + 100,t} \right)$, 上述两种初始条件下对应的解析解可参见文献[14].
图3和图4给出了两种不同初始条件下几个时刻里分量$u$的不同数值结果. 观察图3和图4可知, HSS-CPSPH结果与解析解相符合, 即使在较长模拟时间里; 提出的HSS-CPSPH方法能准确捕捉到孤立子传播中奇异现象; HSS-CPSPH方法可以准确得到一定时间里孤立子波的传播过程.
图 4 初始条件2下, 在4个不同时刻三孤立子波函数$\left| u \right|$的传播过程　(a) $t=0$; (b) $t = 20$; (c) $t = 30$; (d) $t = 50$
Figure4. Solitary wave propagation process of $\left| u \right|$ at different time with initial condition 2: (a) $t=0$; (b) $t = 20$; (c) $t = 30$; (d) $t = 50$.

图 3 初始条件1下, 在4个不同时刻孤立波函数$\left| u \right|$的传播过程　(a) $t = 0$; (b) $t = 20$; (c) $t = 30$; (d) $t = 50$
Figure3. Solitary wave propagation process of $\left| u \right|$ at different time with the initial condition 1: (a) $t = 0$; (b) $t = 20$; (c) $t = 30$; (d) $t = 50$.

2 -->

4.3.二维第一类边界下孤立波传播

为验证提出的HSS-CPSPH方法求解带第一类边界非线性薛定谔问题的准确性, 以及进一步体现它与文献[15]给出的粒子方法相比具有的优点. 本小节考虑区域${\left[ {0,\;2{\text{π}}} \right]^2}$内的二维非线性薛定谔方程^[8]

$\begin{split}{\rm{i}}\dfrac{{\partial u}}{{\partial t}}=& - \dfrac{1}{2}( {{u_{xx}}+{u_{yy}}} ) + V( {x,y} )u +{| u |^2}u, \\& ( {x,y} ) \!\in\!{[ {0,2{\text{π}}} ]^2},\end{split}$

初值条件为${u_0}\left( {x,y} \right) = \sin x\sin y$, 第一类边界条件可以通过方程解析解得到, 其中$V\left( {x,y} \right) = 1 -$$ {\sin ^2}x{\sin ^2}y$. 对应方程解析解为${u_{{\rm{exact}}}}\left( {x,y,t} \right) = \sin x$$\sin y\exp \left( { - {\rm{i}}2t} \right)$.
表4和表5列出了两种粒子方法模拟上述方程的误差和收敛速度. 通过表4和表5可知, 在模拟第一类边值NLS方程时, 本文提出的HSS-CPSPH方法较文献[15]中的粒子方法具有较小误差和较快收敛速度. 此外, 为了体现本文算法通过并行提高计算效率的必要性, 模拟过程中基于MPI并行算法, 采用了多个CPU进行计算机模拟, 取粒子数为66049 (${257^2}$), 时间步长为${\rm{d}}t = $$5 \times {10^{{\rm{ - }}5}}$. 通过实际计算模拟知, 用一个CPU串行计算时, 运行第1步的时间约为49.89 s, 之后100步平均每步约1.5 s; 当采用4个CPU时, 运行第1步的时间约为12.0586 s, 之后100步平均每步约0.25 s. 可以看出, 运行第1步时4个CPU的计算时间基本上是1个CPU的计算时间的1/4, 计算步数增加后, 4个CPU平均每步的计算时间略小于1个CPU平均每步计算时间的1/4. 这与理想的1/4有一定偏差, 可能原因有下面几点: 1) 同一节点上不同CPU之间通讯消耗时间非常小; 2) 此处为二维模拟区域粒子数不是很多的情况; 3) 得到的计算时间有舍入误差. 当模拟粒子数增大(如5.2节三维区域模拟中涉及几百万及以上粒子), CPU个数增加到在不同节点上时, 会使得平均每步的计算时间增加, 不同节点上通讯消耗时间延长, 此时经实际模拟得到的结论会与实际计算理论相符合, 即随着CPU数量增加, 计算效率的提高倍数要明显小于CPU增加的倍数(详见5.2节并行计算效率分析).

时间t	SS-ICPSPH	HSS-CPSPH
0.5	9.131 × 10–4	4.512 × 10–4
1	1.828 × 10–3	8.135 × 10–4
2	3.658 × 10–3	1.623 × 10–3

表4$h = {\text{π}}/64$时, 三个不同时刻两种方法的最大误差${e_{\rm{m}}}$
Table4.Error ${e_{\rm{m}}}$ obtained using two different methods at three times ($h = {\text{π}}/64$).

	$h = {\text{π}}/32$			$h = {\text{π}}/64$			$h = {\text{π}}/128$
	${e_{\rm{m}}}$	$o{r_{\rm{\alpha }}}$		${e_{\rm{m}}}$	$o{r_{\rm{\alpha }}}$		${e_{\rm{m}}}$	$o{r_{\rm{\alpha }}}$
SS-ICPSPH	7.553 × 10–3	—		1.828 × 10–3	2.046		4.316 × 10–4	2.082
HSS-CPSPH	4.534 × 10–3	—		8.135 × 10–4	2.476		1.379 × 10–4	2.560

表5t = 1时不同空间步长情况下两种粒子方法的误差和收敛阶
Table5.Error and order of convergence by different methods at t = 1 and different h.

为了更好地体现提出的基于分裂格式的粒子方法模拟NLS/GP问题的高效、准确性, 本节选取了一个二维周期性无解析解NLS问题和两个BEC下带外旋转项的GP问题, 分别预测了周期边界下孤立子传播的奇异特性和BEC下量子化涡旋变化过程, 并将模拟结果与其他数值方法^[8,15]计算结果做对比. 模拟中, 也通过三维算例的计算机模拟展示了采用基于MPI并行算法提高计算效率的必要性.
2

-->

5.1.二维周期边界无解析解NLS方程

考虑区间$\left[ {0,2{\text{π}}} \right] \times \left[ {0,2{\text{π}}} \right]$上的非线性薛定谔方程^[3]

${\rm{i}}{u_t} + \Delta u + \beta {\left| u \right|^2}u = 0,$

初值条件为$u\left( {x,y,0} \right) = \left( {1 + \sin x} \right)\left( {2 + \sin y} \right)$, 参数$\beta = 1$, 周期边界长度为$\left( {2{\text{π}},2{\text{π}}} \right)$.
上述NLS方程是带初始和周期边界下孤立子传播中出现奇异特性的典型问题, 目前该方程无解析解, 常被用来验证提出算法模拟带周期边界下孤立子传播过程的可靠性(见文献[3]). 图5给出了初始时刻和$t = 0.108$时刻三种数值方法的模拟结果. 从图5可知, 在$t = 0.108$时出现明显的孤立子波奇异特性, 三种数值结果相符合, 进一步表明本文提出的粒子方法捕捉孤立子波传播中奇异现象是可靠的.
图 5 在两个不同时刻不同数值方法得到的$\left| u \right|$等值线图 (a) $t = 0$; (b) $t = 0.108$
Figure5. Contours of $\left| u \right|$ obtained using different methods at two different times: (a) $t = 0$; (b) $t = 0.108$.

2 -->

5.2.三维Gross-Pitaevskii方程

考虑三维BEC下带旋转项的Gross-Pitaevskii方程^[15]

$\begin{aligned}{\rm{i}}\dfrac{{\partial u}}{{\partial t}} =\;& - \dfrac{1}{2}\left( {{u_{xx}} + {u_{yy}} + {u_{zz}}} \right) + V\left( {x,y,z} \right)u \\ &+ \beta {\left| u \right|^2}u - \varOmega {L_z}u,\left( {x,y,z} \right) \in {\left[ {- 10,10} \right]^3},\end{aligned}$

初始条件为
$u\left( {x,y,z,0} \right) = \dfrac{{{{\left( {{\gamma _x}{\gamma _y}{\gamma _z}} \right)}^{1/4}}}}{{2{{\text{π}}^{3/4}}}}{{\rm{e}}^{ - \left( {{\gamma _x}{x^2} + {\gamma _y}{y^2} + {\gamma _z}{z^2}} \right)/2}}$.
其中, $V\left( {x,y,z,t} \right) = \left( {{\gamma _x}{x^2} + {\gamma _y}{y^2} + {\gamma _z}{z^2}} \right)/2$, 参数$\varOmega = 0.75$, $\left( {{\gamma _x},{\gamma _y},{\gamma _z},\beta } \right) = $ $\left( {1.2,0.8,2,50} \right)$.
为表明HSS-CPSPH方法模拟三维GP问题的可靠性和展示BEC下不同截面上波函数值随时间的变化过程, 图6给出了沿$y$轴方向上两种不同数值曲线, 图7给出了由HSS-CPSPH方法得到的3个不同时刻沿$u\left( {0,y,z} \right)$和$u\left( {x,y,0} \right)$两个截面上波函数变化等值线. 图6中本文方法得到的变化曲线与HO-SSFDM方法^[15]结果一致, 体现了HSS-CPSPH方法模拟预测三维GP问题是有效的; 波函数的峰值随时间演化出现先变小后变大的趋势, 但波峰值的总体变化趋势是慢慢变小. 通过观察图7可知, 不同截面上随时间演化波函数值的分布是不同的、复杂的, 且等值线峰值随时间延长逐渐变小.
图 6 不同时刻${\left| u \right|^2}$沿y轴$\left( {x = 0,z = 0} \right)$变化曲线
Figure6. Curve of ${\left| u \right|^2}$ along y-axis $\left( {x = 0,z = 0} \right)$ at different time.

图 7 在3个不同时刻${\left| u \right|^2}$在不同截面上的等值线　(a) $\left( {0,y,z} \right)$截面; (b)$\left( {x,y,0} \right)$截面
Figure7. Contour of ${\left| u \right|^2}$ along different profile at different time: (a) $\left( {0,y,z} \right)$; (b) $\left( {x,y,0} \right)$.

为进一步体现运用本文粒子方法模拟高维GP问题时引入MPI并行算法以减少CPU消耗时间的必要性, 对固定粒子数和增加粒子数两种情况采用不同CPU个数的计算效率进行了分析(见表6和表7). 本文并行计算机采用Red Hat Enterprise Linux 5.8 × 86_64操作系统, MPI类型为Intel MPI Toolkits 4.0.3, 总共有17个IBM BladeCenter HS22计算节点, 每个节点包括12个两路六核CPU, CPU型号为Intel Xeon 6C X5650, 主频为2.67 GHz. 表6列出了固定粒子数${161^3}$ (约四百多万)增加CPU个数情况下运行不同步数时所消耗的CPU时间, 表7列出了不同粒子数下不同CPU情况下, 运行到1000步时平均每步所消耗的CPU时间. 为了更加直观地看出并行计算对计算效率的提高, 此处定义相对加速比S为单节点上并行算法的运行时间/N个节点上并行算法的运行时间 (每个节点上12个CPU). 选用表6中12, 24, 36, 72个CPU在计算到1000步时的时间作为研究对象, 分别对应着1, 2, 3, 6个节点.

CPU数量	步数			相对加速比S
	num = 1	num = 10	num = 1000
2	97805.9	107508	1174728	—
12	16716.9	18516.7	215526.7	—
24	8388.87	9404.37	120284.37	1.792
36	5603.29	6344.98	87524.98	2.462
72	2948.83	3189.24	48564.28	4.438

表6粒子数为${161^3}$时, 不同CPU个数下运行到不同步数所需时间(单位: s)
Table6.Consumed CPU time (unit: s) of different calculated time step with particle number ${161^3}$ at different CPUs.

粒子数	CPU数量
	2	12	24	36	72
${121^3}$	449.55	82.926	45.962	35.000	19.585
${161^3}$	1076.922	198.810	111.90	81.922	47.363
${181^3}$	1558.445	292.711	164.838	120.886	65.437
${201^3}$	2190.921	425.688	235.775	179.856	96.836

表7在不同粒子数下不同CPU个数下, 运行到1000步时平均每步所消耗时间(单位: s)
Table7.The average consumed CPU time (unit: s) of calculated time step 1000 with different particle number and different CPUs.

由表6和表7可以看到, 随着CPU个数的增加, 计算效率提高很快, 但计算效率提高的比值是低于CPU个数增加的比值, 这是因为CPU个数增加后不同CPU上得到结果相互通讯的时间也将增加. 从加速比值上可以明显看出, 在开始使用1, 2, 3个节点时, 计算效率的提升是十分明显的, 但略低于CPU增加的比值, 这是由于不同节点之间通讯消耗一定的时间导致的, 但是当节点数继续增加到6个节点时, 通讯消耗时间也随之增加, 使得计算效率增加减缓, 明显低于CPU增加的比值, 这与计算机相关计算理论也较为符合; 计算机模拟中运行第一步消耗CPU时间比后续运行的平均步消耗CPU时间要长, 且随粒子数增加这个消耗CPU时间的比值要变大, 这是因在计算物理量循环更新前需要标定每个粒子的相邻粒子, 此标定时间随粒子数增加而增加; 在CPU个数不变的情况下, 随着粒子数增加平均每个计算步消耗CPU时间也将变长, 且该变长时间比值与粒子数增加比值是呈非线性的. 由于本文采用的并行计算主要是通过增加CPU个数来提高计算效率, 不同CPU数下的数值近似算法不变, 使得多个CPU下的计算结果与单个CPU的计算结果是一致的, 通过不同CPU下的实际模拟结果对比发现亦是如此, 从而表明多CPU下的求解质量是可靠的.
通过上述模拟分析得知, 本文基于MPI并行计算给出的HSS-CPSPH方法模拟三维GP问题是高效、可靠的.
2 -->

5.3.二维二分量Gross-Pitaevskii方程

为进一步展示提出的HSS-CPSPH方法预测BEC下量子涡旋随时间演化过程的可靠性, 本小节考虑区间${\left[ { - 8,\;8} \right]^2}$上二维二分量BEC下带旋转项的Gross-Pitaevskii方程^[8]

$\left\{ \begin{aligned}& {\rm{i}}\dfrac{{\partial {u_1}}}{{\partial t}} + \dfrac{1}{2}\dfrac{{{\partial ^2}{u_1}}}{{\partial {x^2}}} + \dfrac{1}{2}\dfrac{{{\partial ^2}{u_1}}}{{\partial {y^2}}} + \sigma \left( {{{\left| {{u_1}} \right|}^2} + \varsigma {{\left| {{u_2}} \right|}^2}} \right){u_1} \\ &\quad\quad\, + {w_1}\left( {x,y} \right){u_1} + {\rm{i}}\varTheta \left( {y\dfrac{{\partial {u_1}}}{{\partial x}} - x\dfrac{{\partial {u_1}}}{{\partial y}}} \right) = 0,\\& {\rm{i}}\dfrac{{\partial {u_2}}}{{\partial t}} + \dfrac{1}{2}\dfrac{{{\partial ^2}{u_2}}}{{\partial {x^2}}} + \dfrac{1}{2}\dfrac{{{\partial ^2}{u_2}}}{{\partial {y^2}}} + \sigma \left( {{{\left| {{u_2}} \right|}^2} + \varsigma {{\left| {{u_1}} \right|}^2}} \right){u_2} \\ &\quad\quad\, + {w_2}\left( {x,y} \right){u_2} + {\rm{i}}\varTheta \left( {y\dfrac{{\partial {u_2}}}{{\partial x}} - x\dfrac{{\partial {u_2}}}{{\partial y}}} \right) = 0.\end{aligned} \right.$

初边值条件为

$\left\{ \begin{aligned}& {u_1}\left( {x,y,0} \right) = \dfrac{{x + {\rm{i}}y}}{{\sqrt {\text{π}} }}\exp \left[ { - \dfrac{{\left( {{x^2} + {y^2}} \right)}}{2}} \right],\;\;{\rm{on }}\;\varOmega ,\\& {u_2}\left( {x,y,0} \right) = \dfrac{{x + {\rm{i}}y}}{{\sqrt {\text{π}} }}\exp \left[ { - \dfrac{{\left( {{x^2} + {y^2}} \right)}}{2}} \right],\;\;{\rm{on }}\;\varOmega ,\\& {u_1}\left( {x,y,t} \right) = 0, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{on }}\;\partial \varOmega \times \left[ {0,\;T} \right],\\& {u_2}\left( {x,y,t} \right) = 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{on }}\;\partial \varOmega \times \left[ {0,\;T} \right],\end{aligned} \right.$

其中, ${w_i}\left( {x,y} \right) = 1.5{x^2} + 0.5{y^2},\;i = 1,\;2$, $\varTheta = 0.7, $$\sigma = - 100,\;\varsigma = 0.8 $.
图8给出了两个时刻两种不同数值方法得到的波函数沿$x$轴($y$ = 0.5)的变化, 图9给出了由HSS-CPSPH方法得到的两个时刻3个物理量(第一分量${u_1}$的实部、虚部和模)等值线. 由图8和图9知, HSS-CPSPH结果与SS-FDM结果^[7]相符, 表明本文粒子法模拟二分量GP问题是可靠的; HSS-CPSPH方法能够准确预测出具有二分量BEC下带旋转项量子涡旋随时间的变化过程, 进一步展示了所提方法高效、准确预测GP问题的能力.
图 8 两个不同时刻下$\left| {{u_1}} \right|$沿$x$轴($y$ = 0.5)的变化　(a) $t$ = 0.05; (b) $t$ = 0.25
Figure8. Curve of $\left| {{u_1}} \right|$ along x-axis ($y$ = 0.5) at two different time: (a) $t$ = 0.05; (b) $t$ = 0.25.

图 9 两个不同时刻下t = 0 (第一列)和t = 0.25 (第二列)三个物理量等值线变化　(a1), (a2) $ {\rm Re} ({u_1}) $; (b1), (b2) ${\rm Im} ({u_1})$; (c1), (c2) $\left| {{u_1}} \right|$
Figure9. Contours of three physical quantities at two different times t = 0 (the first row) and t = 0.25 (the second row): (a1), (a2) $\operatorname{Re} ({u_1})$; (b1), (b2) $\operatorname{Im} ({u_1})$; (c1), (c2) $\left| {{u_1}} \right|$.

本文为提高SPH方法模拟含孤立波非线性问题的数值精度和稳定性, 将四阶精度时间分裂格式与修正SPH方法相结合, 并引入基于粒子搜索的MPI并行计算技术, 给出一种能够高效准确地模拟高维NLS/GP方程的HSS-CPSPH方法. 数值模拟中, 考虑了粒子分布均匀和非均匀情况下带有两种不同边界条件的二维和三维NLS/GP方程, 并与解析解或其他数值结果做对比, 分析了所提方法的数值精度、收敛阶及数值预测的可靠性. 所有数值算例表明:
1)提出的HSS-CPSPH方法求解高维非线性薛定谔方程较已有粒子方法具有较高精度和较快收敛速度, 且较网格类方法具有更好、更灵活的应用性;
2)无论在粒子分布均匀还是非均匀情况下, HSS-CPSPH方法模拟NLS方程都具有较高的数值精度;
3)给出的HSS-CPSPH方法成功地预测了周期边界下二维NLS方程描述的孤立波传播中的奇异特性, 准确地展示了BEC下带外旋转项三维单分量和二维二分量GP方程中量子涡旋随时间演化过程.
值得注意的是, 目前未见文献将四阶分裂格式与修正SPH方法耦合, 并引入MPI并行技术对NLS/GP方程进行模拟研究. 本文对不同边界下高维NLS方程的模拟给出的HSS-CPSPH法较已有粒子方法具有较高精度和较快收敛速度, 较网格类方法具有更好、更灵活的拓展应用性, 为高维NLS方程的数值模拟提供了一种高效准确的粒子方法.