【读MFiX源代码】4 run_dem(des_time_march)计算颗粒相（主要步骤DEM_TIME_STEP）(源文件在des_time_march.f)

位置

被谁调用（入口）

run_mfix(dt_loop)(mfix.f)

mfix.f文件中第146行，注意观察输出

输出仅为这三行

run_dem

进入run_dem(仍然在mfix.f)

可以发现主要步骤有三步：
des_time_init
des_time_step
des_time_end

运行完des_time_init会输出第一行
运行完des_time_step会输出第二行
运行完des_time_end会输出第三行

源代码位置

这三个子程序全部在DES_TIME_MARCH这个module里
源文件：des_time_march

des_time_march.f文件结构

结构很简单
只有三个subroutine，而且是contains的
即
des_time_init
des_time_step
des_time_end

源代码分析(des_time_march.f)

module 全局变量

NP是颗粒的总数，从0开始，并且定义成了静态变量（SAVE），即下一次再调用这个module的时候，NP这个变量不会重新定义为0，而是从上一次保存的值开始

FACTOR，这个很关键，用来控制dem_time_step外部循环（run_dem中的）的循环上限。

DTSOLID_TMP临时的固相时间
TMP_WALL临时的壁钟时间（wall time，即现实中的时间）

EXIT_LOOP 逻辑变量，表示是否退出循环

des_time_init

des_time_step（136-246行）

子程序名：SUBROUTINE DES_TIME_STEP(NN)

因为太长人为地分割成三部分
第一部分截止到调用UDF,
第二部分截止到调用邻域搜索
第三部分到结束

第1部分 函数头到call USR1_DES( 136-187)

函数头和变量定义

源代码是04年写的，证明04年之前都不能用DEM(所以96年的 theory guide必定只能用TFM)

SUBROUTINE DES_TIME_STEP(NN)
只引入一个整数NN，且只读

mod_assertion是一个逻辑变量，为true代表NN=1或者NN是NEIGHBOR_SEARCH_N的整数倍

目前mode_assertion=true，因为NN=1,证明进行邻域搜索。

user guide 20.1.0 P207

dem-doc-2012 P16

大致意思就是每隔NEIGHBOR_SEARCH_N个des迭代之后，就会进行一次邻域搜索。

注：dem-doc-2012中说默认值是25，但是目前（19.3.1之后）版本来看，默认值是20

当然还有其他条件来指定是否发生邻域搜索（比如颗粒移动距离）
GUI中的位置：solids->DEM->max steps between neighbor search

153-165行颗粒、流体时间同步的if endif

153-165行的if endif
控制dem循环的时间用的
S_TIME代表离散相的时刻
TIME表示流体相的时刻
DTSOLID表示固相的时间步长
DT表示流体相的时间不长

153行的if表示如果固体相的时刻本身都超过流体相的时刻了，那么直接退出。

163行的if表示如果固体相时刻加上下一个固体相的时间步长，就会超出流体相的下一个时刻，那么直接修剪DTSOLID。以此同步两个时刻

168行if optflag1 .eq. 1整理网格和颗粒数组

DES_SORT_PARTICLE_ARRAYS
DES_GETVALID_FLUID_CELLS
这两个看名字应该分别是整理颗粒数组和得到有效网格，暂时按下不表，后面再研究

optflag1来源于run

173-178行计算碰撞力和曳力等（核心步骤）

CALC_FORCE_DEM是计算颗粒间碰撞力的，重点

CALC_DRAG_DES是计算颗粒与流体之间的曳力的，重点

CALC_DEM_THERMO_WITH_WALL_STL计算颗粒与墙壁之间的热传导

180-187行

CFUPDATEDOLD 更新颗粒速度和位置值(目前未执行)
含义应该是 cf update old。至于cf是什么意思目前不知道。
注意这个DO_OLD目前是false，也就意味着它没有执行。所以真正更新颗粒速度位置的应该是更下面的CFNEWVALUES.

CALC_THERMO_DES计算颗粒上的热源，与反应相关，重点

USR1_DES是在颗粒相计算中的第一个用户自定义程序(还有第0个)
可以看出 CALL_USR这个逻辑变量应该是可以在deck中写的

GUI(model选项卡)中勾选即将其置为true

第2部分 CFNEWVALUES到 CALL NEIGHBOUR（189-211）

更新颗粒位置速度温度 ，核心

CFNEWVALUES 真正的更新颗粒速度和位置，核心步骤

DES_THERMO_NEWVALUES 更新颗粒的温度，核心步骤

195行-199行的DO_NSEARCH感觉是多余的代码，没什么用
另外.neqv.是逻辑不等于（也不知道和普通的.ne.有什么区别，可能逻辑变量只能用这个吧）

附Fortran中的6个关系运算符和5个逻辑运算符
http://blog.sciencenet.cn/blog-783377-903615.html

通过BC注入、移除颗粒，重点

DEM_BCMI和DEM_BCMO为正值的时候，代表用颗粒注入、飞出
MASS_INFLOW_DEM和MASS_OUTFLOW_DEM分别处理因BC造成的颗粒的添加和删除

颗粒交换（？待研究）DESGRID_PIC和DES_PAR_EXCHANGE

与幽灵颗粒有关（待研究）
numPEs是代表并行运算的处理器核心数（逻辑核心）

看条件：进行邻域搜索或者并行或者有周期壁面边界的时候，才运行这两个子程序

邻域搜索NEIGHBOUR，核心

待细致研究，是核心步骤之一。

第3部分 从封装（bin）颗粒到结束

封装颗粒（216-224）

尽管第1个IF在本次运行中跳过了，但是很多情况下我们是不打开DES_EXPLICITYLY_COUPLED的（？用户指南说只有纯动力流动才能用显性的耦合，但是实际运行中发现并非如此）

注释说显示耦合不需要每次都把颗粒封装到网格里（从而大大节约了运算成本），但是这样做无疑让网格信息和插值权重变旧了（stale:旧）

下面三个语句
分别是封装颗粒到网格里
计算插值权重
计算平均场

有待研究

更新固相时刻（226-237）

227行
S_TIME是固相的时刻
这意味着固相的计算基本结束了

229-237行，可以跳过，只对纯颗粒流。也是和时间有关

固相中的第2个UDF: USR2_DES

没什么好说的，就是计算一次固相时间步的最后调用的UDF

结束

IF(ADJUST_PARTITION)这句也可以跳过

des_time_end

待续