第四届全国计算爆炸力学会议论文集
三维弹塑性流体力学欧拉方法动态并行研究
李征 廖振民 梁仙红 刘军 冯其京 何长江
（北京应用物理与计算数学研究所，北京 100088）
摘 要：欧拉方法是流体力学数值模拟的主要方法，适于计算多介质大变形的实际三维问题。针对主要关心的爆轰、侵彻等问题的特点，设计了动态调整计算区域的并行方法，可以随着计算的推进，根据各处理器的计算量，动态调整计算区域以达到较好的负载平衡。本方法主要可以分成两部分，一部分是初始时刻计算区域的分割与数据的定位，另一部分是模拟过程中计算区域的扩充与数据的调整。结果表明，采用动态并行方法可以大大缩短计算时间，更适于实际的应用问题计算。并行测试中，其加速比基本呈线性增长。该动态并行方法可以自然推广到高维多介质欧拉程序的并行计算中。 *
关键词：动态并行；三维；弹塑性流体力学；欧拉方法
Dynamic Parallel Method for Eulerian Method of Elastoplastic
Hydrodynamics in Three Dimension
LI Zheng, LIAO Zhen-min, LIANG Xian-hong, LIU Jun, FENG Qi-jing, HE Chang-jiang
（Institute of Applied Physics and Computational Mathematics, Beijing 100088）
Abstract: Eulerian method is a main numerical method in hydrodynamics, which is suitable for the problems with multi-component and large deformation. As the feature of the problems to be simulated, such as detonation and penetration, the parallel method is designed to adjust the computational domain dynamically to get better load balance. Dynamic parallel method can be separated into two parts: one is division of initial computational domain and location of the data, the other is expansion of the computational domain and adjustment of the data location. Dynamic parallel code can greatly shorten computational time and be preferable in simulating actual problems. The speedup of the program is linear in parallel test. The method can be popularized to parallel program of other multi-component high dimension Eulerian methods naturally.
Key words: dynamic parallel; three dimension; elastoplastic hydrodynamics; Eulerian method
引 言
欧拉方法是流体力学数值模拟的主要方法，在高温高压物理研究、高速爆炸现象分析、燃烧机理和爆炸效应应用等方面都发挥了重要的作用。欧拉方法的网格固定，不存在因物质界面大变形而产生网格扭曲，计算结果的网格依赖性较小，适宜推广到高维进行实际问题的模拟。三维实际问题往往计算规模巨大，迫切需要研究可缩短计算周期的并行程序。 *基金项目：计算物理国家重点实验室基金（9140C690101070C69）
作者简介：李征（1979—），女，助理研究员，理学硕士，E-mail: [email protected]
三维弹塑性流体力学欧拉方法动态并行研究
武器物理所研究的包括爆轰、弹片飞散、射流的形成与侵彻等一系列复杂问题。大部分问题在物质实体周围经常有成片空网格（或者水下计算问题中的静态水），特别是侵彻问题，弹体与靶的距离比较远，弹体以外的计算空间比较大。采用欧拉方法的程序计算，如果对整个计算区域进行计算，扰动到来前对成片空网格、静态水或者靶的计算是无效的。并且，如果采用传统的固定的区域分解的并行方法，容易出现负载不平衡的情况。为了提高计算效率，尽量使并行达到负载平衡，设计了动态并行方法：初始计算区域取一个包含初始状态的最小长方体，并随着扰动的发展扩大计算区域，计算区域动态变化。选择考虑负载平衡的动态调整空间区域方法，随着计算的推进，根据计算区域动态调整各处理器的计算范围以达到较好的负载平衡。考虑到程序的可移植性和可扩展性，并行程序采用目前国际上主流的MPI （Message Passing Interface）的并行设计平台。
文章主要包括三部分的内容：第一部分，动态并行方法，具体可以分为初始计算区域的分割与数据定位，以及计算区域的扩充与数据的调整。第二部分，程序流程图。第三部分，并行测试，动态并行程序与串行程序的计算结果完全一致，且其加速比呈线性增长，并且计算时间与固定的区域分解的情况进行了比较，本方法并行程序更节省时间。
1 动态并行方法
动态并行中，计算区域随扰动的发展而动态变化，适于爆轰、射流等问题的计算。动态并行方法主要可以分为两部分：初始时刻计算区域的分割与数据定位；计算区域的扩充和数据的调整。
1.1 初始计算区域的分割与数据定位
初始时刻计算区域的范围是与所模拟的问题相关的，根据问题的不同，初始区域的选择有两个标准。第一个，以将问题的基本状态描述清楚的最小立方空间。第二个，如果速度、压力等物理量的扰动已经到达初始计算空间的边界，则需要适当扩充这个最小立方空间，扩充的大小与数值方法相关，遵循的原则是选择的计算区域要完成初始物理量的运算。
确定计算区域后需要面对两个问题：计算区域如何分配到每台处理机；以及每台处理机上的数据如何定位。下面以x 方向为例进行说明。
设模型问题在x 方向的总网格数im ，初始时刻的计算范围[ia , ib ]只是整个计算区域的一部分，网格数imx =ib -ia +1≤im 。x 方向分配了icpu 台处理器，编号为1，2，…，icpu 。如果网格数imx 是icpu 的整数倍，将网格平均分配于每台处理器；否则，每个处理器对应的网格数imn 可能有两种情况，⎢⎡imx ⎤⎡imx ⎤+1（[]表示小数舍掉的取整数方法）和⎢。为简便起见，公式可以统一记作（1）。 ⎥⎥icpu icpu ⎣⎦⎣⎦
⎧⎡imx ⎤less mesh number⎪⎢⎥⎪⎣icpu ⎦imn =⎨ −1imx ⎤+1more mesh number⎪⎡
⎪⎣icpu ⎥⎦⎩⎢ （1）
具有较多网格数的处理机为nn =mod(imx,icpu ) 个。选择靠近扩充区域的处理器具有较多网格，并在较少网格数处理机的两端设置标志ipr 1和ipr 2，扩充过程中随着各处理器网格数的变化而改变位置。图1给出了计算区域两端都是空区时的网格分布情况。这样，通过计算每个处理器的计算区域长度imn (n ) ，就确定了其负责的计算子区域范围[in (n ) +1, in (n ) +imn (n )]。其中，
if n =1⎧ia −1, in (n ) =⎨ −+−>in (n 1) imn (n 1), if n 1⎩ （2）
第四届全国计算爆炸力学会议论文集
图1 初始计算区域的分割示意图
Fig1 division of initial computational domain
（a ） （b ）
图2 数据定位示意图
Fig2 location of the data
每个处理器定义在x 方向的网格量数组[1∶imm ]，。这里，imm =(im -1)/icpu +1+mk ，其中，mk 是某一正整数，提供了动态并行扩充数据时数据移动的空间。每个处理器内部网格的范围为[ian , ibn ]，长度imn (n )=ibn -ian +1，网格量数组中空闲网格的长度Δ=imm -imm (n ) 。按照问题扩充方式不同，空闲网格的分布可以分为两种情况。 ● 求解区域一端扩充网格时，空闲网格位于靠近扩充区域的一端。图2（a ），扩充区域位于x 方向
的负方向或者正方向，空闲网格位于求解区域的负方向或者正方向。
● 求解区域两端扩充网格时，空闲网格位于区域的两端，如图2（b ），且两端空闲网格的长度随着
处理机序号改变。具体情况见表1，其中，
δ=⎢⎥ icpu +1⎣⎦⎡Δ⎤ （3） 表1 空闲网格分布表 Tab1 distribution of the empty mesh 处理器号 1 2 …… Δ-δ Δ-2δδ 2δ …… …… icpu-1 icpu Δ-(icpu-1)*δ Δ-icpu*δ d1 d2 (icpu-1)*δ icpu*δ 1.2 计算区域的扩充与数据的调整
图3 计算区域扩充示意图
Fig3 expansion of the computational domain
三维弹塑性流体力学欧拉方法动态并行研究
区域扩充可以根据问题不同扩充不同的网格，可以是空网格、静态水或者靶。扩充的标志为扰动靠近计算空间的边界，以保证计算区域内可以正确完成相应的计算。每个方向的区域扩充可以分为正负两个方向，下面以x 负方向扩充为例进行说明，如图3所示。
● 处理机1在负方向增加网格量，在正方向向相邻的处理器传送ibn 网格量。处理机2至处理机
ipr1-1在负方向接受相应的网格量，并向正向相邻处理机发送ibn 网格量，处理机ipr1在负方向接受传输的网格量。
● 计算区域扩充后的求解空间：
处理器1至处理器ipr1-1的求解空间：ian =ian -1，ibn =ibn -1
处理器ipr1的求解空间：ian =ian -1，ibn =ibn
● 更改短网格数处理器标志ipr1：ipr1=iprl +1。
扩充过程中，如果处理器1到ipr1中ian=1时，扩充无法进行，必须先调整各处理器的计算子区域的数据位置，即重新计算Δ和δ，按初始数据定位的原则重新定位，并将网格量数据移到相应的位置；然后继续扩充的过程。
2 程序流程图
图4是动态并行程序流程图，一般的弹塑性流体力学程序通过初始时刻和每一步末尾的两部分的动态并行方法实现程序的动态并行，随着计算区域的变化调整各处理器的计算范围，可以达到较好的负载平衡。
图4 动态并行程序流程图
Fig4 flow chart of dynamic parallel program
3 并行程序测试
3.1 动态并行程序与固定区域分割方式的并行程序计算比较
设计水下爆轰问题为测试模型，如图5所示，炸药内部是真空，外面包了铁壳和铝板，弹体外面是水。选择网格数为100×100×100，计算100时间步。动态并行程序初始区域包括部分计算空间，随着炸药爆轰逐步增大计算空间，其结果与串行程序完全一样。同样的问题采用固定区域分割的方式，即传统的欧拉方法（简称静态并行），在初始时刻对整个计算空间进行分割并分配到各个处理器
第四届全国计算爆炸力学会议论文集
上去，可以看到动态并行程序可以大大节省计算时间。图6给出了动态与静态并行程序在不同数目处理机上的时间比较。针对本测试模型而言，动态并行计算所需的时间仅为静态并行计算时间的1/5，计算效率提高了5倍。正如测试数据所表明的那样，动态并行程序对于很多实际物理问题的计算比静态并行程序节省很多时间，更适合三维实际问题的应用。
图5 水下爆轰问题模型
Fig5 model of detonation in the water
3.2 动态并行程序加速比测试
图7给出了不同数目处理器的并行状况。虽然动态并行方法中，扩充计算区域和调整数据位置需要增加一些计算量和数据通讯，但并行效果依然很理想，其加速比基本呈线性增长，8台处理器获得了高达7.31的加速。
图6 动态并行与静态并行计算时间比较图 7 动态并行程序加速比测试表
Fig6 time comparison of dynamic and static parallel code Fig7 speedup of the dynamic parallel code
3.3 计算结果
本水下爆轰问题计算到3.0时刻的结果如图6所示。
三维弹塑性流体力学欧拉方法动态并行研究
密度 压力
图6 水下爆轰问题结果
Fig6 results of detonation in the water
结 论
该动态并行算法适用于一步欧拉方法和两步欧拉方法，一维、二维的并行计算，也可以自然推广到高维多介质问题的计算。动态并行算法适于计算射流形成，飞片计算等从局部到整体的过程，与传统的固定区域分割的并行方法相比，可以大大缩短计算时间，提高计算效率。并行效果也比较好，加速比基本呈线性增长。
参考文献
[1] 何长江，于志鲁，冯其京等，三维整体求解弹塑性流体力学欧拉数值方法[J]，ZW-J-2002135，
2003.
He Chang-jiang, Yu Zhi-lu, Feng Qi-jing,etc., Euler numerical method of elastic-plasticity hydrodynamics in 3D[J], ZW-J-2002135，2003.
[2] 何长江，冯其京，于志鲁等，三维多介质欧拉方法混合网格的有限体积计算格式[A]，中国工程
物理研究院科技年报2003，绵阳：原子能出版社，2003：428-429.
He Chang-jiang, Feng Qi-jing, Yu Zhi-lu, etc., Finite Volume scheme of 3D Multi-material Eulerian Method on Mixed Cells[A], Annual report of China academy of engineering physics 2003, Mianyang: atomic energy publishing company, 2003: 428-429.
[3] 李征，何长江，袁仙春，三维整体求解弹塑性流体力学欧拉方法程序的优化与并行[J]，
ZW-J-2005134，2006.
Li Zheng, He Chang-jiang, Yuan Xian-chun, Parallel and optimization of elastic-plasticity hydrodynamics Eulerian method non-splitting scheme in three dimension[J], ZW-J-2005134，2006.
ZW-J-2006104，[4] 李征，廖振民，三维整体求解弹塑性流体力学欧拉方法程序的MPI 并行研究[J]，
2006.
Li Zheng, Liao Zhen-min, MPI parallel for non-splitting Eulerian method of elastoplastic hydrodynamics in three dimension[J], ZW-J-2006104, 2006.
第四届全国计算爆炸力学会议论文集
三维弹塑性流体力学欧拉方法动态并行研究
李征 廖振民 梁仙红 刘军 冯其京 何长江
（北京应用物理与计算数学研究所，北京 100088）
摘 要：欧拉方法是流体力学数值模拟的主要方法，适于计算多介质大变形的实际三维问题。针对主要关心的爆轰、侵彻等问题的特点，设计了动态调整计算区域的并行方法，可以随着计算的推进，根据各处理器的计算量，动态调整计算区域以达到较好的负载平衡。本方法主要可以分成两部分，一部分是初始时刻计算区域的分割与数据的定位，另一部分是模拟过程中计算区域的扩充与数据的调整。结果表明，采用动态并行方法可以大大缩短计算时间，更适于实际的应用问题计算。并行测试中，其加速比基本呈线性增长。该动态并行方法可以自然推广到高维多介质欧拉程序的并行计算中。 *
关键词：动态并行；三维；弹塑性流体力学；欧拉方法
Dynamic Parallel Method for Eulerian Method of Elastoplastic
Hydrodynamics in Three Dimension
LI Zheng, LIAO Zhen-min, LIANG Xian-hong, LIU Jun, FENG Qi-jing, HE Chang-jiang
（Institute of Applied Physics and Computational Mathematics, Beijing 100088）
Abstract: Eulerian method is a main numerical method in hydrodynamics, which is suitable for the problems with multi-component and large deformation. As the feature of the problems to be simulated, such as detonation and penetration, the parallel method is designed to adjust the computational domain dynamically to get better load balance. Dynamic parallel method can be separated into two parts: one is division of initial computational domain and location of the data, the other is expansion of the computational domain and adjustment of the data location. Dynamic parallel code can greatly shorten computational time and be preferable in simulating actual problems. The speedup of the program is linear in parallel test. The method can be popularized to parallel program of other multi-component high dimension Eulerian methods naturally.
Key words: dynamic parallel; three dimension; elastoplastic hydrodynamics; Eulerian method
引 言
欧拉方法是流体力学数值模拟的主要方法，在高温高压物理研究、高速爆炸现象分析、燃烧机理和爆炸效应应用等方面都发挥了重要的作用。欧拉方法的网格固定，不存在因物质界面大变形而产生网格扭曲，计算结果的网格依赖性较小，适宜推广到高维进行实际问题的模拟。三维实际问题往往计算规模巨大，迫切需要研究可缩短计算周期的并行程序。 *基金项目：计算物理国家重点实验室基金（9140C690101070C69）
作者简介：李征（1979—），女，助理研究员，理学硕士，E-mail: [email protected]
三维弹塑性流体力学欧拉方法动态并行研究
武器物理所研究的包括爆轰、弹片飞散、射流的形成与侵彻等一系列复杂问题。大部分问题在物质实体周围经常有成片空网格（或者水下计算问题中的静态水），特别是侵彻问题，弹体与靶的距离比较远，弹体以外的计算空间比较大。采用欧拉方法的程序计算，如果对整个计算区域进行计算，扰动到来前对成片空网格、静态水或者靶的计算是无效的。并且，如果采用传统的固定的区域分解的并行方法，容易出现负载不平衡的情况。为了提高计算效率，尽量使并行达到负载平衡，设计了动态并行方法：初始计算区域取一个包含初始状态的最小长方体，并随着扰动的发展扩大计算区域，计算区域动态变化。选择考虑负载平衡的动态调整空间区域方法，随着计算的推进，根据计算区域动态调整各处理器的计算范围以达到较好的负载平衡。考虑到程序的可移植性和可扩展性，并行程序采用目前国际上主流的MPI （Message Passing Interface）的并行设计平台。
文章主要包括三部分的内容：第一部分，动态并行方法，具体可以分为初始计算区域的分割与数据定位，以及计算区域的扩充与数据的调整。第二部分，程序流程图。第三部分，并行测试，动态并行程序与串行程序的计算结果完全一致，且其加速比呈线性增长，并且计算时间与固定的区域分解的情况进行了比较，本方法并行程序更节省时间。
1 动态并行方法
动态并行中，计算区域随扰动的发展而动态变化，适于爆轰、射流等问题的计算。动态并行方法主要可以分为两部分：初始时刻计算区域的分割与数据定位；计算区域的扩充和数据的调整。
1.1 初始计算区域的分割与数据定位
初始时刻计算区域的范围是与所模拟的问题相关的，根据问题的不同，初始区域的选择有两个标准。第一个，以将问题的基本状态描述清楚的最小立方空间。第二个，如果速度、压力等物理量的扰动已经到达初始计算空间的边界，则需要适当扩充这个最小立方空间，扩充的大小与数值方法相关，遵循的原则是选择的计算区域要完成初始物理量的运算。
确定计算区域后需要面对两个问题：计算区域如何分配到每台处理机；以及每台处理机上的数据如何定位。下面以x 方向为例进行说明。
设模型问题在x 方向的总网格数im ，初始时刻的计算范围[ia , ib ]只是整个计算区域的一部分，网格数imx =ib -ia +1≤im 。x 方向分配了icpu 台处理器，编号为1，2，…，icpu 。如果网格数imx 是icpu 的整数倍，将网格平均分配于每台处理器；否则，每个处理器对应的网格数imn 可能有两种情况，⎢⎡imx ⎤⎡imx ⎤+1（[]表示小数舍掉的取整数方法）和⎢。为简便起见，公式可以统一记作（1）。 ⎥⎥icpu icpu ⎣⎦⎣⎦
⎧⎡imx ⎤less mesh number⎪⎢⎥⎪⎣icpu ⎦imn =⎨ −1imx ⎤+1more mesh number⎪⎡
⎪⎣icpu ⎥⎦⎩⎢ （1）
具有较多网格数的处理机为nn =mod(imx,icpu ) 个。选择靠近扩充区域的处理器具有较多网格，并在较少网格数处理机的两端设置标志ipr 1和ipr 2，扩充过程中随着各处理器网格数的变化而改变位置。图1给出了计算区域两端都是空区时的网格分布情况。这样，通过计算每个处理器的计算区域长度imn (n ) ，就确定了其负责的计算子区域范围[in (n ) +1, in (n ) +imn (n )]。其中，
if n =1⎧ia −1, in (n ) =⎨ −+−>in (n 1) imn (n 1), if n 1⎩ （2）
第四届全国计算爆炸力学会议论文集
图1 初始计算区域的分割示意图
Fig1 division of initial computational domain
（a ） （b ）
图2 数据定位示意图
Fig2 location of the data
每个处理器定义在x 方向的网格量数组[1∶imm ]，。这里，imm =(im -1)/icpu +1+mk ，其中，mk 是某一正整数，提供了动态并行扩充数据时数据移动的空间。每个处理器内部网格的范围为[ian , ibn ]，长度imn (n )=ibn -ian +1，网格量数组中空闲网格的长度Δ=imm -imm (n ) 。按照问题扩充方式不同，空闲网格的分布可以分为两种情况。 ● 求解区域一端扩充网格时，空闲网格位于靠近扩充区域的一端。图2（a ），扩充区域位于x 方向
的负方向或者正方向，空闲网格位于求解区域的负方向或者正方向。
● 求解区域两端扩充网格时，空闲网格位于区域的两端，如图2（b ），且两端空闲网格的长度随着
处理机序号改变。具体情况见表1，其中，
δ=⎢⎥ icpu +1⎣⎦⎡Δ⎤ （3） 表1 空闲网格分布表 Tab1 distribution of the empty mesh 处理器号 1 2 …… Δ-δ Δ-2δδ 2δ …… …… icpu-1 icpu Δ-(icpu-1)*δ Δ-icpu*δ d1 d2 (icpu-1)*δ icpu*δ 1.2 计算区域的扩充与数据的调整
图3 计算区域扩充示意图
Fig3 expansion of the computational domain
三维弹塑性流体力学欧拉方法动态并行研究
区域扩充可以根据问题不同扩充不同的网格，可以是空网格、静态水或者靶。扩充的标志为扰动靠近计算空间的边界，以保证计算区域内可以正确完成相应的计算。每个方向的区域扩充可以分为正负两个方向，下面以x 负方向扩充为例进行说明，如图3所示。
● 处理机1在负方向增加网格量，在正方向向相邻的处理器传送ibn 网格量。处理机2至处理机
ipr1-1在负方向接受相应的网格量，并向正向相邻处理机发送ibn 网格量，处理机ipr1在负方向接受传输的网格量。
● 计算区域扩充后的求解空间：
处理器1至处理器ipr1-1的求解空间：ian =ian -1，ibn =ibn -1
处理器ipr1的求解空间：ian =ian -1，ibn =ibn
● 更改短网格数处理器标志ipr1：ipr1=iprl +1。
扩充过程中，如果处理器1到ipr1中ian=1时，扩充无法进行，必须先调整各处理器的计算子区域的数据位置，即重新计算Δ和δ，按初始数据定位的原则重新定位，并将网格量数据移到相应的位置；然后继续扩充的过程。
2 程序流程图
图4是动态并行程序流程图，一般的弹塑性流体力学程序通过初始时刻和每一步末尾的两部分的动态并行方法实现程序的动态并行，随着计算区域的变化调整各处理器的计算范围，可以达到较好的负载平衡。
图4 动态并行程序流程图
Fig4 flow chart of dynamic parallel program
3 并行程序测试
3.1 动态并行程序与固定区域分割方式的并行程序计算比较
设计水下爆轰问题为测试模型，如图5所示，炸药内部是真空，外面包了铁壳和铝板，弹体外面是水。选择网格数为100×100×100，计算100时间步。动态并行程序初始区域包括部分计算空间，随着炸药爆轰逐步增大计算空间，其结果与串行程序完全一样。同样的问题采用固定区域分割的方式，即传统的欧拉方法（简称静态并行），在初始时刻对整个计算空间进行分割并分配到各个处理器
第四届全国计算爆炸力学会议论文集
上去，可以看到动态并行程序可以大大节省计算时间。图6给出了动态与静态并行程序在不同数目处理机上的时间比较。针对本测试模型而言，动态并行计算所需的时间仅为静态并行计算时间的1/5，计算效率提高了5倍。正如测试数据所表明的那样，动态并行程序对于很多实际物理问题的计算比静态并行程序节省很多时间，更适合三维实际问题的应用。
图5 水下爆轰问题模型
Fig5 model of detonation in the water
3.2 动态并行程序加速比测试
图7给出了不同数目处理器的并行状况。虽然动态并行方法中，扩充计算区域和调整数据位置需要增加一些计算量和数据通讯，但并行效果依然很理想，其加速比基本呈线性增长，8台处理器获得了高达7.31的加速。
图6 动态并行与静态并行计算时间比较图 7 动态并行程序加速比测试表
Fig6 time comparison of dynamic and static parallel code Fig7 speedup of the dynamic parallel code
3.3 计算结果
本水下爆轰问题计算到3.0时刻的结果如图6所示。
三维弹塑性流体力学欧拉方法动态并行研究
密度 压力
图6 水下爆轰问题结果
Fig6 results of detonation in the water
结 论
该动态并行算法适用于一步欧拉方法和两步欧拉方法，一维、二维的并行计算，也可以自然推广到高维多介质问题的计算。动态并行算法适于计算射流形成，飞片计算等从局部到整体的过程，与传统的固定区域分割的并行方法相比，可以大大缩短计算时间，提高计算效率。并行效果也比较好，加速比基本呈线性增长。
参考文献
[1] 何长江，于志鲁，冯其京等，三维整体求解弹塑性流体力学欧拉数值方法[J]，ZW-J-2002135，
2003.
He Chang-jiang, Yu Zhi-lu, Feng Qi-jing,etc., Euler numerical method of elastic-plasticity hydrodynamics in 3D[J], ZW-J-2002135，2003.
[2] 何长江，冯其京，于志鲁等，三维多介质欧拉方法混合网格的有限体积计算格式[A]，中国工程
物理研究院科技年报2003，绵阳：原子能出版社，2003：428-429.
He Chang-jiang, Feng Qi-jing, Yu Zhi-lu, etc., Finite Volume scheme of 3D Multi-material Eulerian Method on Mixed Cells[A], Annual report of China academy of engineering physics 2003, Mianyang: atomic energy publishing company, 2003: 428-429.
[3] 李征，何长江，袁仙春，三维整体求解弹塑性流体力学欧拉方法程序的优化与并行[J]，
ZW-J-2005134，2006.
Li Zheng, He Chang-jiang, Yuan Xian-chun, Parallel and optimization of elastic-plasticity hydrodynamics Eulerian method non-splitting scheme in three dimension[J], ZW-J-2005134，2006.
ZW-J-2006104，[4] 李征，廖振民，三维整体求解弹塑性流体力学欧拉方法程序的MPI 并行研究[J]，
2006.
Li Zheng, Liao Zhen-min, MPI parallel for non-splitting Eulerian method of elastoplastic hydrodynamics in three dimension[J], ZW-J-2006104, 2006.