By Xiaobing Shen, Miao Chen

Updated:

CitCom(for California Institute of Technology Convection in the Mantle)是利用有限元方法解决二维笛卡尔坐标系下地幔热对流模拟问题，编写环境为C语言。CitCom最早版本是由Louis Moresi在1990年编写的，随后在钟时杰、Clint Conrad、Eh Tan等人的努力下不断优化。

在学习CitCom代码之初，大部分初学者都被要求完成CitCom二维地幔对流benchmark测试以作为对代码及如何运行代码的基本检验，这也是本文的主要内容。

准备

CitCom代码编译

对于本文中源码压缩包名为2Dcompressible.tar，如果源码压缩包名称不同，请自行修改命令对象的名称。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 当前目录
$ pwd
/home/shen/benchmark_test

$ ls
2Dcompressible.tar

# 解压
$ tar -xvf 2Dcompressible.tar
$ mkdir Test

# 将input文件移动到Test文件夹内
$ mv input_example Test/

# 进行CitCom源码进行编译，生成可执行文件citcom.x
$ make citcom.x

正常情况下运行 make citcom.x 命令后会出现一大段C语言链接、编译的输出语句，可能会出现几个warning警告，这只是因为代码中存在小的问题，但并不会影响程序的运行，成功生成citcom.x可执行文件后，对于源代码的编译完成。

input文件参数设置

对于输入文件的参数设置，在这里只简单介绍其中在测试benchmark中会被修改的几个参数含义及其设置，其他参数再以后的博文中再做详细的介绍。

1
2
3
$ cd Test
$ vim input_example
# 查看input文件内容

对input文件中较为重要或需要修改的参数及其含义如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
rayleigh=1.0e7         # 地幔模拟计算中最主要的控制参数，非常重要

datafile="CASE_test"   # 运算结果文件将存储于该路径下的文件夹中
maxstep=100000         # 设置最大运行步数
storage_spacing=50     # 设置输出步数间隔

restart=0              # 重启设置当需接着之前某一已经运行过的case时，将其值有"0"改为"1"
restart_timesteps=0    # 设置重新运行时由之前的哪一步开始
oldfile="CASE_TV11H"   # 因重启运算读取该文件夹路径下存储的原计算结果数据

dimenx=2.0             # 该值与dimenz组合，可设置对流box的length：high比值
dimenz=1.0

z_grid_layers=4        # Z方向网格分布，主要用来网格加密，在benchmark中是均匀网格，可改为2
zz=0.,0.08,0.92,1.0    # 若为均匀网格，则改为： 0,1.0
nz=1,17,113,129        # 若为均匀网格，则改为： 1,129

x_grid_layers=2        # 同上，但为X方向的网格分布
xx=0.,2.0
nx=1,257
mgunitx=4
mgunitz=2
levels=7

Viscosity=system
VISC_UPDATE=on
rheol=4                # visvosity计算公式的选择，在Viscosity_structures.c中可查看
TDEPV=on
CHEMDEPV=off

num_mat=4              # Z方向地层物质分层数，一般设为4，即自地表至核幔边界认为存
                       # 在4层主要分层：岩石圈；岩石圈-410km；410km-660km；
                       # 660km-CMB。

viscE=11.51293,11.51293,11.51293,11.51293 # 活化能对viscosity影响系数
viscZ=0.0,0.0,0.0,0.0                     # 深度对viscosity影响系数
viscT=0.0,0.0,0.0,0.0                     # 温度对viscosity影响系数
visc0=1.0e0,1.e0,1.e0,1.e0                
VMIN=on visc_min=1.0e-4
VMAX=on visc_max=1.0e04

# BOUNDARY CONDITIONS

	topvbc=0			# 顶部速度边界条件开关
		topvbxval=0.0	# 速度边界条件X方向值
		topvbyval=0.0	# 速度边界条件Y方向值
	botvbc=0 	        # 同上，不过为底部边界条件
		botvbxval=0.0
		botvbyval=0.0	
#
	toptbc=1 bottbc=1		# 温度边界条件，顶部为0，底部为1，均为无量纲化温度
	toptbcval=0.0 bottbcval=1.0	#

# DIMENSIONAL INFORMATION (BENCHMARK)

	Ts=273                # 表面量纲化温度，单位开尔文
	ReferenceT=3000.0     # CMB superadiabatic temperature in K 
	refvisc=5.0e21        # 参考viscosity
	density=3400.0        # 密度
	thermdiff=1.0e-4      # 热扩散系数
	gravacc=9.8           # 重力梯度
	thermexp=3.0e-5       # 热膨胀系数
	cp=1200               # 比热容
	wdensity=0.0
	layerd=2870.0 # km    # 量纲化的box深度

Rayleigh数设置

rayleigh数是CitCom程序中重要的控制参数，其物理公式为：

其中\(\rho\),g,\(\alpha\),\({\Delta}T\),D,\(\kappa\),\(\eta\)分别是密度，重力梯度，热膨胀系数，上下表面温差，box深度，热扩散系数，参考viscosity。
以case 1a为例，修改Ra=1e4即可。

分辨率设置

分辨率设置主要影响计算可信度和计算效率，以Z方向为例分辨率设置参数nz，mgunitz，levels要满足公式：

Viscosity设置

在viscosity_structures.c文件中有多个viscosity计算公式可供选择，可通过改变’rheol’的值改变选择的viscosity公式，以最简单的受温、压影响的viscosity公式：

可通过修改visc0，viscT，viscE参数改变viscosity剖面。对于case 2需要修改viscosity参数，根据文章给定的参数进行修改。

量纲化参数设置

根据文章中Table1中修改ReferenceT，refvisc，density，thermdiff，thermexp，layerd等参数。

运行CitCom

1
2
3
4
5
6
7
8
9
10
11
12
# 当前目录
$ pwd
/home/shen/benchmark_test/Test

# citcom.x在"/home/shen/benchmark_test/"中
# input_example在"/home/shen/benchmark_test/Test"中

# 新建结果文件夹
$ mkdir CASE_test

# 运行程序
$ nohup ../citcom.x ./input_example &

运行上面的命令后程序就在运行中，可在生成的”nohup.out”文件查看运行日志。

按照上述方法，修改相应参数，重复出Case1a，Case1b，Case1C，Case2a，Case2b结果。

输出文件格式

CitCom输出文件主要包括：log,coord,ave,temp,heating等文件。
文件命名规则：”filename.timestep”,例如：“ave.1000”代表第1000步计算输出的ave文件。
log文件主要是文件计算运行过程中输出的信息；
coord文件主要是坐标参数等信息；
temp文件主要是温度、粘度、不同方向速度信息；
ave文件主要是温度、速度，粘度、密度等各个参数水平均值信息，绘制各参数随Z方向变化曲线图数据即从该文件中提取，故仅详细介绍ave文件内容格式：

ave文件：（以X,Z分辨率nx，nz为例）注：参数值均为无量纲化值
第1行：nz timestep time 表面热流(Nu值) 底面热流 放射性生热量 底部温度边界条件开关
第2至nz+1行(前5列)：Z方向坐标 水平平均温度 水平平均速度 水平平均粘度 水平平均密度
第nz+2至nz+nx+1行：X方向坐标及box表面各个参数值
第nz+nx+2行：略

运行结果

从运行结果文件中提取Nu、velo等参数，利用GMT绘制图件，并与文章中给出结果进行对比，误差在3%以内则基本运行正确。
Case1:

Ra=1e4

Ra=1e5

Ra=1e6




Case2a:

Error~Grid:

Ra~Nu, Ra~Vrm关系的定量化计算结果

Reference
Blankenbach, B., et al. “A benchmark comparison for mantle convection codes.” Geophysical Journal International 98.1 (1989): 23-38.