谱模式四维变分资料同化并行算法设计
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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
α 为最佳步长
更新扰动量
∆来自百度文库 a = ∆x k a
enddo
n +1 n xb = xb + ∆x a
内循环结束 更新背景场 外循环结束
enddo
3．并行算法设计
从上述的算法可看出， 四维变分资料同化的主要计算量在三个模式的计算。 我们原始模 [4] 式采用 T106 谱模式。由于 T106 谱模式的并行计算是比较成熟的 ，在此基础上发展起来的 切线性模式和伴随模式可以采用与 T106 谱模式一致的并行方案。根据 Fourier 空间的计算 和谱空间的计算关于波无相关性， 格点空间的计算关于纬圈和经圈无相关性。 为了减少积分
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
5
[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 时刻积分
参考文献
[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.
Fourier 逆变换
Fourier 系数 1
格点空间
勒让德逆变换 Fourier 正变换
格点值
Fourier 空间 2
勒让德正变换
Fourier 系数 2
谱空间 图 1 数值预报模式中间变量存储
谱系数
以上是三个模式单独运行时的主要并行策略。但实际上，在四维变分资料同化中，三个 模式的运行是紧密联系在一起的。 伴随模式的运行是由观测资料残差驱动的， 这就需要原始 预报模式计算出的各个时刻模式状态量 x i 和切线性模式计算出的各个时刻增量 ∆x i 。另外， 从伴随理论和上面的叠代增量算法可知， 伴随代码运行过程中， 需要原始状态量和与原始状 [5] 态量有数据依赖关系的中间变量。这些变量值的获取有两种 ：重新计算和从内存或文件中 读取。对于数值预报模式这样大的计算量而言，单独使用任何一种方法均不可取，而应该是 在综合考虑计算量和数据存取速度之上的折中。于是，我们采取的方案为：把原始预报模式 中 4 个空间（见图 1）计算出的预报变量值存取在本地内存或文件中，在伴随模式运行到该 空间时， 再把它提取出来。 而该空间中其它中间变量值则重新计算。 对切线性模式也是如此。 这样作的优越性主要有： （1）由于在各个空间之间需要进行 Fourier 变换或勒让德变换，计算量很大。采用该 方案后，可以避免对中间预报变量进行变换计算。 （2）原始模式的并行方案尽管只需要两次数据重分配，但全局通信需要花费时间多， 每个处理机都必须与其它处理机进行数据交换， 这会降低并行效率。 而现在把每个空间 的原始模式计算出来的中间变量均存放在本地内存或文件中， 就避免了与其它处理机进 行数据交换。 （3〕三个模式中，各个空间的计算采用同样的数据分配方案，并行效率高，所以在各 空间中的其它中间变量可以通过重新计算获取。
计算梯度
−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
4．结论
4
四维变分同化的并行算法设计与数值预报模式并行算法设计最根本的区别就是三个模 式之间数据的传递， 而这是实现四维变分同化业务化运行的一个关键问题。 尤其是数值预报 模式和切线性模式之间、数值预报模式和伴随模式之间数据传递。这些数据量相当大，很难 用消息传递机制实际。 而只能利用切线性和伴随模式与预报模式的内在联系， 使得保存在各 个节点上的模式状态量应与该节点所要作的计算紧密相连， 尽量避免访问其他节点保存的数 据。这些需要在实际运行中反复调整，并结合高性能计算机结构特征，才能达到最佳效果。
谱模式四维变分资料同化并行算法设计
曹小林 宋君强 孙安香 （国防科技大学计算机学院，中国长沙，410073, xiaolincao@163.com） 摘 要 四维变分资料同化并行算法的关键在于三个模式之间的数据传递。 本文充分利用切 线性模式和伴随模式与原始预报谱模式的内在联系， 对切线性模式和伴随模式采用与原始预 报模式一致的并行方案。 对于分布式存储并行计算机， 把原始预报谱模式计算出来每时间步、 每个空间的预报变量值或系数存储在本地内存或文件以实现三个模式之间的数据传递。 这包 括谱模式的 Fourier 系数 1、格点场值、Fourier 系数 2、谱系数。采用该方案可避免在切线 性模式和伴随模式中对预报量进行 Fourier 变换和勒让德变换，更重要的是避免对这些量进 行全局通信。 关键字 并行算法、切线性模式、伴随模式、原始预报谱模式、数据传递
1. 引言
四维变分资料同化为数值预报模式提供一个质量场和流场基本平衡的最佳初始场。近 几年，在国外，四维变分资料同化的业务化大大提高了数值预报准确率。国内也开展这方面 的研究， 但离业务化仍然有一段距离。 这主要是由于许多关键技术未突破或达不到业务化要 求，如伴随模式的生成、背景场协方差矩阵的形成、最优化算法、并行算法设计等。本文就 主要讨论四维变分资料同化并行算法设计方案。 四维变分资料同化计算量很大，大约为数值预报模式的 50－100 倍。如果采用迭代增 量方法，可减少到 5－25 倍[1]。因此只有采用并行算法，才能满足业务化要求。它最主要的 运算量是三个模式的运行， 包括原始数值预报模式、 数值预报模式的切线性模式及其伴随模 式。如何处理好三个模式的并行运算及其间的数据交换是四维变分资料同化并行化的关键。
2.迭代增量方法
四维变分资料同化在给定周期内的观测资料和先验信息之间寻求一个最佳平衡 。 它通 过由动力模式和观测资料构造的泛函极小化来实现。 其核心思想是： 寻求一个最佳的模式解， 其预报演变轨迹与给定周期 T 内的观测资料拟合的最好， 并能与初始时刻的背景场信息吻合 [3] 地很好。具体数学公式 为：
∆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
[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 时的 背景模式（或预报）状态。一般由初始时刻的预报求得。它是对初始条件的优化以及对过去
3
过程中数据重分配的次数， 在同一个时间积分步里相关的 Legendre 逆变换和 Legendre 正变 换按一维波进行一维剖分将数据分配到各个处理机。Fourier 逆变换和 Fourier 正变换按纬 圈进行一维剖分将数据分配到各个处理机以实现 并行计算。对波数在各处理机上按螺旋式 分布，各纬圈则交替分布在各处理机上，使得负载平衡。这种并行计算方法在一个时间积分 需要进行两次数据重分配， 格点空间的并行计算前进行一次数据重分配； 然而谱空间按波进 行分布式并行计算， 谱空间并行计算之前必须进行第二次数据重分配。 采用上述并行计算方 3 法，一个时间步有两次数据交换，各处理机的通信复杂性为 O(M /p)，理论上可以利用的最 大处理机数为 O(M)。由于通信复杂性与处理机数目成反比，在处理机数目不大于 O(M)的情 况下，该并行计算方法具有良好的可扩展性。对于现行业务分辨率的计算规模，完全可以满 足业务的实时性要求。 Fourier 空间 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
α 为最佳步长
更新扰动量
∆来自百度文库 a = ∆x k a
enddo
n +1 n xb = xb + ∆x a
内循环结束 更新背景场 外循环结束
enddo
3．并行算法设计
从上述的算法可看出， 四维变分资料同化的主要计算量在三个模式的计算。 我们原始模 [4] 式采用 T106 谱模式。由于 T106 谱模式的并行计算是比较成熟的 ，在此基础上发展起来的 切线性模式和伴随模式可以采用与 T106 谱模式一致的并行方案。根据 Fourier 空间的计算 和谱空间的计算关于波无相关性， 格点空间的计算关于纬圈和经圈无相关性。 为了减少积分
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
5
[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 时刻积分
参考文献
[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.
Fourier 逆变换
Fourier 系数 1
格点空间
勒让德逆变换 Fourier 正变换
格点值
Fourier 空间 2
勒让德正变换
Fourier 系数 2
谱空间 图 1 数值预报模式中间变量存储
谱系数
以上是三个模式单独运行时的主要并行策略。但实际上，在四维变分资料同化中，三个 模式的运行是紧密联系在一起的。 伴随模式的运行是由观测资料残差驱动的， 这就需要原始 预报模式计算出的各个时刻模式状态量 x i 和切线性模式计算出的各个时刻增量 ∆x i 。另外， 从伴随理论和上面的叠代增量算法可知， 伴随代码运行过程中， 需要原始状态量和与原始状 [5] 态量有数据依赖关系的中间变量。这些变量值的获取有两种 ：重新计算和从内存或文件中 读取。对于数值预报模式这样大的计算量而言，单独使用任何一种方法均不可取，而应该是 在综合考虑计算量和数据存取速度之上的折中。于是，我们采取的方案为：把原始预报模式 中 4 个空间（见图 1）计算出的预报变量值存取在本地内存或文件中，在伴随模式运行到该 空间时， 再把它提取出来。 而该空间中其它中间变量值则重新计算。 对切线性模式也是如此。 这样作的优越性主要有： （1）由于在各个空间之间需要进行 Fourier 变换或勒让德变换，计算量很大。采用该 方案后，可以避免对中间预报变量进行变换计算。 （2）原始模式的并行方案尽管只需要两次数据重分配，但全局通信需要花费时间多， 每个处理机都必须与其它处理机进行数据交换， 这会降低并行效率。 而现在把每个空间 的原始模式计算出来的中间变量均存放在本地内存或文件中， 就避免了与其它处理机进 行数据交换。 （3〕三个模式中，各个空间的计算采用同样的数据分配方案，并行效率高，所以在各 空间中的其它中间变量可以通过重新计算获取。
计算梯度
−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
4．结论
4
四维变分同化的并行算法设计与数值预报模式并行算法设计最根本的区别就是三个模 式之间数据的传递， 而这是实现四维变分同化业务化运行的一个关键问题。 尤其是数值预报 模式和切线性模式之间、数值预报模式和伴随模式之间数据传递。这些数据量相当大，很难 用消息传递机制实际。 而只能利用切线性和伴随模式与预报模式的内在联系， 使得保存在各 个节点上的模式状态量应与该节点所要作的计算紧密相连， 尽量避免访问其他节点保存的数 据。这些需要在实际运行中反复调整，并结合高性能计算机结构特征，才能达到最佳效果。
谱模式四维变分资料同化并行算法设计
曹小林 宋君强 孙安香 （国防科技大学计算机学院，中国长沙，410073, xiaolincao@163.com） 摘 要 四维变分资料同化并行算法的关键在于三个模式之间的数据传递。 本文充分利用切 线性模式和伴随模式与原始预报谱模式的内在联系， 对切线性模式和伴随模式采用与原始预 报模式一致的并行方案。 对于分布式存储并行计算机， 把原始预报谱模式计算出来每时间步、 每个空间的预报变量值或系数存储在本地内存或文件以实现三个模式之间的数据传递。 这包 括谱模式的 Fourier 系数 1、格点场值、Fourier 系数 2、谱系数。采用该方案可避免在切线 性模式和伴随模式中对预报量进行 Fourier 变换和勒让德变换，更重要的是避免对这些量进 行全局通信。 关键字 并行算法、切线性模式、伴随模式、原始预报谱模式、数据传递
1. 引言
四维变分资料同化为数值预报模式提供一个质量场和流场基本平衡的最佳初始场。近 几年，在国外，四维变分资料同化的业务化大大提高了数值预报准确率。国内也开展这方面 的研究， 但离业务化仍然有一段距离。 这主要是由于许多关键技术未突破或达不到业务化要 求，如伴随模式的生成、背景场协方差矩阵的形成、最优化算法、并行算法设计等。本文就 主要讨论四维变分资料同化并行算法设计方案。 四维变分资料同化计算量很大，大约为数值预报模式的 50－100 倍。如果采用迭代增 量方法，可减少到 5－25 倍[1]。因此只有采用并行算法，才能满足业务化要求。它最主要的 运算量是三个模式的运行， 包括原始数值预报模式、 数值预报模式的切线性模式及其伴随模 式。如何处理好三个模式的并行运算及其间的数据交换是四维变分资料同化并行化的关键。
2.迭代增量方法
四维变分资料同化在给定周期内的观测资料和先验信息之间寻求一个最佳平衡 。 它通 过由动力模式和观测资料构造的泛函极小化来实现。 其核心思想是： 寻求一个最佳的模式解， 其预报演变轨迹与给定周期 T 内的观测资料拟合的最好， 并能与初始时刻的背景场信息吻合 [3] 地很好。具体数学公式 为：
∆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
[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 时的 背景模式（或预报）状态。一般由初始时刻的预报求得。它是对初始条件的优化以及对过去
3
过程中数据重分配的次数， 在同一个时间积分步里相关的 Legendre 逆变换和 Legendre 正变 换按一维波进行一维剖分将数据分配到各个处理机。Fourier 逆变换和 Fourier 正变换按纬 圈进行一维剖分将数据分配到各个处理机以实现 并行计算。对波数在各处理机上按螺旋式 分布，各纬圈则交替分布在各处理机上，使得负载平衡。这种并行计算方法在一个时间积分 需要进行两次数据重分配， 格点空间的并行计算前进行一次数据重分配； 然而谱空间按波进 行分布式并行计算， 谱空间并行计算之前必须进行第二次数据重分配。 采用上述并行计算方 3 法，一个时间步有两次数据交换，各处理机的通信复杂性为 O(M /p)，理论上可以利用的最 大处理机数为 O(M)。由于通信复杂性与处理机数目成反比，在处理机数目不大于 O(M)的情 况下，该并行计算方法具有良好的可扩展性。对于现行业务分辨率的计算规模，完全可以满 足业务的实时性要求。 Fourier 空间 1