谱模式四维变分资料同化并行算法设计

谱模式四维变分资料同化并行算法设计
曹小林 宋君强 孙安香 （国防科技大学计算机学院，中国长沙，410073, xiaolincao@） 摘 要 四维变分资料同化并行算法的关键在于三个模式之间的数据传递。 本文充分利用切 线性模式和伴随模式与原始预报谱模式的内在联系， 对切线性模式和伴随模式采用与原始预 报模式一致的并行方案。 对于分布式存储并行计算机， 把原始预报谱模式计算出来每时间步、 每个空间的预报变量值或系数存储在本地内存或文件以实现三个模式之间的数据传递。 这包 括谱模式的 Fourier 系数 1、格点场值、Fourier 系数 2、谱系数。采用该方案可避免在切线 性模式和伴随模式中对预报量进行 Fourier 变换和勒让德变换，更重要的是避免对这些量进 行全局通信。 关键字 并行算法、切线性模式、伴随模式、原始预报谱模式、数据传递
计算梯度
−1 ∇J = B −1 (∆x k + ∆x) a
endif do i=nt,1,-1 get x i ∆x i
伴随模式反向积分 读取模式状态量
Q i = H i (x i + ∆x i ) − y i
−1 J O = J O + QT i O Qi −1 qi = qi + H* i O Qi
∆x1 a =0
∆x = 0 ∇J = 0
do k=1,ninner if(k.ne.1)then 内循环
∆x 0 = ∆x k a
do i=1,nt
2
get x i
读取模式状态量 切线性模式从 i − 1 时刻到 i 时刻积分
∆x i = M ′ i ∆x i −1
enddo
JB =
1 −1 −1 (∆x k + ∆x) T B −1 (∆x k + ∆x) a a 2
4．结论
4
四维变分同化的并行算法设计与数值预报模式并行算法设计最根本的区别就是三个模 式之间数据的传递， 而这是实现四维变分同化业务化运行的一个关键问题。 尤其是数值预报 模式和切线性模式之间、数值预报模式和伴随模式之间数据传递。这些数据量相当大，很难 用消息传递机制实际。 而只能利用切线性和伴随模式与预报模式的内在联系， 使得保存在各 个节点上的模式状态量应与该节点所要作的计算紧密相连， 尽量避免访问其他节点保存的数 据。这些需要在实际运行中反复调整，并结合高性能计算机结构特征，才能达到最佳效果。
1
信息的综合归纳的一个先验估计。B 是描述预报质量 x 0 的预报（或背景）误差协方差矩阵。

b
y i 是在 t i 时刻的观测。 O i 为 t i 时刻的观测资料误差协方差矩阵。 H i 是观测资料转换算子。
这利用到了两个假设[3]： a.预测模式可表示成各中间预测步的乘积
x i = x(t i ) = M (t 0 , t i )x 0 = M i M i −1 " M 1 x 0
q i −1 = M * i qi
enddo
伴随模式从 i 时刻到 i − 1 时刻反向积分
J = JO + JB
如果 J < eps exit 迭代终止条件 计算梯度
∇J = ∇J + q 0
∆x = α∇J
k −1 ∆x k + ∆x a = ∆x a
α 为最佳步长
更新扰动量
∆x a = ∆x k a
enddo
n +1 n xb = xb + ∆x a
内循环结束 更新背景场 外循环结束
enddo
3．并行算法设计
从上述的算法可看出， 四维变分资料同化的主要计算量在三个模式的计算。 我们原始模 [4] 式采用 T106 谱模式。由于 T106 谱模式的并行计算是比较成熟的 ，在此基础上发展起来的 切线性模式和伴随模式可以采用与 T106 谱模式一致的并行方案。根据 Fourier 空间的计算 和谱空间的计算关于波无相关性， 格点空间的计算关于纬圈和经圈无相关性。 为了减少积分
参考文献
[1]Rantakokko J., Strategies for parallel variational data assimilation, at http://www.tdb.uu.sc/ PSCI/Results/List.html, 1996. [2]Courtier, P and O Talagrand, Variational assimilation of meteorological observations with the direct and adjoint vorticity equation.Part II: Numerical results, Q. J. R. Meteorol. Soc, 113, 1329—1368, 1987. [3]Courtier P., Thepaut, J. N., and Hollingsworth, A., A strategy for operational implemtntation of 4D-Var, using an incremental approach. Quart. J. Roy. Meteor. Soc. 120, 1367—1387, 1994.. [4]孙安香，宋君强，T106 谱模式的并行计算，第 6 届全国并行计算学术论文集，234—241， 2000。 [5] Giering R. and Kaminski T. Recipes for Adjoint Code Construction. ACM Trans. On Math. Software, 24(4): 437—474, 1998.
[2]
J=
1 1 n T ∆x 0 B −1 ∆x 0 + ∑ Q i O i−1Q T i = JB + JO 2 2 i =0
（1a）
Q i = H i (M (t 0 , t i )x b + M ′(t 0 , t i )∆x) − y i min J (∆x 0 )
（1b） （1c）
其中， x 0 = x b + ∆x 0 是模式初始状态或控制变量。 x b 是在同化周期初始时刻 t 0 时的 背景模式（或预报）状态。一般由初始时刻的预报求得。它是对初始条件的优化以及对过去
[1]
∑ M' (t0 , ti )* H i*O −1Q i
n
如果采用最速下降法进行求解。那么迭代增量法具体步骤 为： 给定 x b 、 y i 、 B 、 O i
x1 b = xb
do j=1,nouter 外循环
x 0 = x bj
do i=1,nt
x i = M i x i −1
enddo
预报模式从 i − 1 时刻到 i 时刻积分
2.迭代增量方法
四维变分资料同化在给定周期内的观测资料和先验信息之间寻求一个最佳平衡 。 它通 过由动力模式和观测资料构造的泛函极小化来实现。 其核心思想是： 寻求一个最佳的模式解， 其预报演变轨迹与给定周期 T 内的观测资料拟合的最好， 并能与初始时刻的背景场信息吻合 [3] 地很好。具体数学公式 为：
b.切线性假设
（2）
x i = M (t 0 , t i )x 0 = M (t 0 , t i )(x b + ∆x) = M (t 0 , t i )x b + M ' (t 0 , t i )∆x
目标函数的梯度为
（3）
∇J B = B −1 ∆x 0
∇J o =
(4)
i =0 n * * * −1 (5) = ∑ M′ i Hi O Qi 0 "M′ i =0 * * −1 *H *O −1Q " ′ * H1*O −1Q1 + " + M ′ = M′ 0 H 0 O Q 0 + M1 n n n
3
过程中数据重分配的次数， 在同一个时间积分步里相关的 Legendre 逆变换和 Legendre 正变 换按一维波进行一维剖分将数据分配到各个处理机。Fourier 逆变换和 Fourier 正变换按纬 圈进行一维剖分将数据分配到各个处理机以实现 并行计算。对波数在各处理机上按螺旋式 分布，各纬圈则交替分布在各处理机上，使得负载平衡。这种并行计算方法在一个时间积分 需要进行两次数据重分配， 格点空间的并行计算前进行一次数据重分配； 然而谱空间按波进 行分布式并行计算， 谱空间并行计算之前必须进行第二次数据重分配。 采用上述并行计算方 3 法，一个时间步有两次数据交换，各处理机的通信复杂性为 O(M /p)，理论上可以利用的最 大处理机数为 O(M)。由于通信复杂性与处理机数目成反比，在处理机数目不大于 O(M)的情 况下，该并行计算方法具有良好的可扩展性。对于现行业务分辨率的计算规模，完全可以满 足业务的实时性要求。 Fourier 空间 1
1. 引言
四维变分资料同化为数值预报模式提供一个质量场和流场基本平衡的最佳初始场。近 几年，在国外，四维变分资料同化的业务化大大提高了数值预报准确率。国内也开展这方面 的研究， 但离业务化仍然有一段距离。 这主要是由于许多关键技术未突破或达不到业务化要 求，如伴随模式的生成、背景场协方差矩阵的形成、最优化算法、并行算法设计等。本文就 主要讨论四维变分资料同化并行算法设计方案。 四维变分资料同化计算量很大，大约为数值预报模式的 50－100 倍。如果采用迭代增 量方法，可减少到 5－25 倍[1]。因此只有采用并行算法，才能满足业务化要求。它最主要的 运算量是三个模式的运行， 包括原始数值预报模式、 数值预报模式的切线性模式及其伴随模 式。如何处理好三个模式的并行运算及其间的数据交换是四维变分资料同化并行化的关键。
5
Fourier 逆变换
Fourier 系数 1
格点空间
勒让德逆变换 Fourier 正变换