CN101158985A - 一种超维度河流动力学自适应并行监测的方法 - Google Patents

Abstract

本发明公开一种超维度河流动力学自适应并行监测的方法，包括以下步骤：将超维度数据输入系统，并按照数据所处维度的不同进行分类；建立基于特征型高分辨率数值算法的超维度非结构网格河流动力学模型；依据超维度流体分裂格式的高效并行计算方法，进行维度内的计算和维度间的计算；将计算区域分成多个子区域，每个子区域映射到并行体系结构上的一个计算结点上，结点间的通信采用标准的消息传递接口，在自适应网格中采用计算与通信的重叠并行优化技术，变量关于空间的计算独立。本发明方法将超维度河流动力学在自适应网格中，进行高效的分裂格式的并行计算，并同时处理维度的变化，其实现河流的监测方便及时，准确度高。

Description

一种超维度河流动力学自适应并行监测的方法

技术领域

本发明主要涉及到一种水利的动力学监测方法，特别是涉及一种基于自适应网格的并行计算方法，同时也涉及到一种基于超维度河流动力学的建模预处理方法和新型的基于特征型高分辨率数值算法和超维度非结构网格河流动力学模型的河流监测方法。

背景技术

目前的数字流域领域对于河流动力学的研究主要集中在一维和二维，以产生河道、河网等的模拟，并进行规划，在三维的流体模拟方面还不成熟和深入，常用的流体三维分裂格式的高效并行计算，其流体连续方程微分形式为(ρ)(t)+·(ρμ)＝0，其中，ρ为流体密度，μ为速度。

设计算区域为Q_T∈(0，II)×(0，JJ)×(0，KK)×(0，T)，用平行线x＝x_i＝iΔx(i＝1，2，...，II)，y＝y_j＝jΔy(i＝1，2，...，JJ)，z＝z_k＝kΔz(i＝1，2，...，KK)和t＝tⁿ＝nΔt对计算区域进行分割，这里iΔx＝ii，jΔy＝jj，KΔz＝KK，NΔt＝T，I，J，K，和N为正整数，Δx，Δy，Δz和Δt分为空间网格步长和时间步长。假设等空间步长和时间步长。

在给出连续方程形式上的分裂格式，由第n时间排上的值求第n+1时间排上的值，其计算过程分为3步，以方向离散的差分格式为例，其他方向类同。

$f_{i,j.k}^n={\mathrm{\α}}_1{f}_{i+3,j,k}^n+{\mathrm{\α}}_2{f}_{i+2,j,k}^n+{\mathrm{\α}}_3{f}_{i+1,j,k}^n+{\mathrm{\α}}_4{f}_{i,j,k}^n+{\mathrm{\α}}_5{f}_{i-1,j,k}^n+{\mathrm{\α}}_6{f}_{i-2,j,k}^n+{\mathrm{\α}}_7{f}_{i-3,j,k}^n$

这里，f代表连续方程中的密度，α代表x方向上差分方程的系数。

将整个流动区域分割成N个子区域分配给N个CPU计算，把子区域的初始流场信息、几何信息(网格坐标、标识号)分别装载入各子区域对应的CPU的内存中，在每一个CPU中启动计算进程，由主进程调度各CPU的计算。在每一次全场的扫描过程中，由各CPU完成子区域的计算并在边界完成数据交换(即各CPU间的通信)，由主进程收集全场数据完成收敛准则判别，并按需要进行写盘等其它操作。

在物理模型和数值算法已定的情况下，计算速度主要取决于CPU性能、CPU个数、内存、CPU-内存访问带宽、结点互连带宽、网格质量及分区质量等。每一个特定问题、每一台特定机器对应于一个最佳分区数，大量的实践会对同一类问题总结出一个最佳网格数/CPU比值。分区数过多，CPU间通信量增大，分区数增大到一定程度会反而降低计算速度；分区数过少，没有充分利用更多的CPU参与计算，也会影响计算速度；分区质量差，各CPU负载不均匀，CPU有等待现象，也影响速度。

通常成熟的商业软件包，如FLUENT，都至少提供两种分区方法：自动分区和手动分区。自动分区的好处在于操作简便，适宜工程应用，但各CPU的负载平衡不一定能保证。手动分区恰好相反。值得注意的是，对于稳态流动，在任何一个计算时刻内存中都只保留整场的一系列数据信息，故在计算过程中只存在CPU-内存、CPU-CPU(或计算结点-计算结点)之间的通讯，无需访问磁盘。

但对于非定常问题，每一时间步均产生一系列场信息，每一时间步计算结束后必须将各计算结点的场信息会聚起来，执行内存到磁盘的写操作，此时通信量骤增，计算速度无疑会放慢。

在流动机理研究方面，区域分割并行方法仍然广为使用。如美国劳伦斯利佛莫尔实验室(Lawrence Livermore National Laboratory，USA)研究湍流及其非定常性时就采用了此法。由于湍流尺度极小，故网格分辨率必须小于湍流尺度，从而使整个计算区域的网格数达到数十亿。

随着河流动力学的发展和大流域大规模场景的要求，三维的处理已经无法满足水利工程的研究要求。需要注入超维度的信息，即通过对超维度数据进行维度划分提取，构建具有一定格式的数据表和数据库，以方便并行处理，例如：地理空间维，即河流的地理位置；时间维，反映连续变化的属性；相关信息维，包括水利信息(如水位、雨量、流量、工情资料等)和与水利相关的社经信息(如人口、经济、土地等)等等，而随着研究的深入和不同需求的特殊情况，维度的变化和增多是不可避免的。

因此，现有技术还有待于改进和发展。

发明内容

本发明的目的之一是为了克服目前常用的河流动力学监测处理的速度不够快，计算量增大下的计算有效性降低，以及对现实状况的模拟效果和能力有限的问题。

本发明的目的之二是为了克服低纬度信息的贫乏，将超维度数据经简单预处理后输入系统，并按照数据所处维度的不同进行分类，建立新型的基于特征型高分辨率数值算法的超维度非结构网格河流动力学模型，并运用维度分裂和自适应网格并行处理，以提供较为精确的监测依据。

为了实现上述目的，本发明提供了一种超维度河流动力学自适应并行监测的方法，该方法包括如下步骤：

A、将超维度数据输入系统，并按照数据所处维度的不同进行分类；

B、建立基于特征型高分辨率数值算法的超维度非结构网格河流动力学模型；

C、依据超维度流体分裂格式的高效并行计算方法，进行维度内的计算和维度间的计算；

D、将计算区域分成多个子区域，每个子区域映射到并行体系结构上的一个计算结点上，结点间的通信采用标准的消息传递接口，在自适应网格中采用计算与通信的重叠并行优化技术，变量关于空间的计算独立。

所述的方法，其中，还包括在所述步骤A之前还设置有步骤：

A0、建立多维度的基本方程模型。

所述的方法，其中，所述步骤C的并行计算过程还包括：

C1、计算每个处理机上的关于某个方向的内边界点，即定义在每个处理机的子区域的边界点，不包括求解区域的自然边界点；

C2、将该方向上的内边界点的数据进行消息传递；

C3、计算内点，即定义在每个处理机上的子区域的非边界点。

所述的方法，其中，所述步骤D中自适应网格的生成过程采用超三角形渐次插入法，包括步骤如下：

D1、形成全域超三角形；

D2、离散边界，按顺序插入边界点；

D3、搜索包含该边界点的三角形，将此三角形分解为3个三角形；重复上述D2～D3步骤，直到完成所有边界；

D4、删除所有包含一个或多个超三角形顶点的三角形；

D5、根据网格尺度控制，加密三角网，直到形成符合尺度的三角形网格；

D6、网格优化。

本发明所提供的一种超维度河流动力学自适应并行监测的方法，适用于基于超维度河流信息的有效提取处理分析，能够切实反映超维度分析下的河流变化规律；通过运用基于特征型高分辨率数值算法的超维度非结构网格河流动力学模型，进行超维度河流动力学的映射和计算，其处理速度快；将超维度河流动力学在自适应网格中，进行高效的分裂格式的并行计算，并同时处理维度的变化，其实现河流的监测方便及时，准确度高。

附图说明

图1是本发明方法的处理流程图；

图2为本发明方法的硬件系统示意图；

图3为本发明方法的逻辑结构示意图。

具体实施方式

参照附图，将对本发明所述方法做详细的说明。

本发明超维度河流动力学自适应并行监测的方法，其核心的处理构思是，将超维度数据经简单预处理后输入通用处理系统，并按照数据所处维度的不同进行分类，建立新型的基于特征型高分辨率数值算法的超维度非结构网格河流动力学模型，依据超维度流体分裂格式的高效并行计算方法，进行维度内的计算和维度间的计算，将计算区域分成多个子区域，每个子区域映射到并行体系结构上的一个计算结点上，结点间的通信采用标准的MPI，在自适应网格中采用计算与通信的重叠并行优化技术，变量关于空间的计算是独立的。

所述河流动力学模型的建立过程包括，基于卫星图片和水文地理数据的析取，将多种信息条件以边界条件的形式加入到三维流动基本方程中，从而形成该模型。

本发明方法所采用的硬件为一通用计算机或高性能计算机HPC，如图2所示，其包括一3D交互操作显示平台，高性能计算平台，大规模数据实时采集处理平台，所述显示平台是为了三维显示系统内容并可以进行实时互操作；所述高性能计算平台提供了强大的计算能力和快捷的处理速度；采集处理平台可以快速的接受检测水文地理数据等，并进行处理等。本发明方法所采用的系统逻辑结构，如图3所示，其包括并行数据和模型库，并行计算系统，以及3D操作交互系统。

如图1所示，为本发明方法的一种优选实施方式的流程图，其包括以下步骤：

步骤1，建立基本方程模型。对于河口、海岸或大型湖泊(水库)，其垂向加速度与重力加速度相比很小，满足静水压强分布的假定，可以略去垂向加速度项。静水压强假定下的三维流动基本方程可写成：

连续方程 $\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}+\frac{\∂w}{\∂z}=0$ (1)

动量方程

$\frac{\∂u}{\∂t}+\frac{\∂\mathrm{uu}}{\∂x}+\frac{\∂\mathrm{uv}}{\∂y}+\frac{\∂\mathrm{uw}}{\∂z}=\mathrm{fv}-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂x}+\frac{1}{\mathrm{\ρ}}(\frac{{\∂\mathrm{\τ}}_{\mathrm{xx}}}{\∂x}+\frac{{\∂\mathrm{\τ}}_{\mathrm{yz}}}{\∂y}+\frac{{\∂\mathrm{\τ}}_{\mathrm{zz}}}{\∂z})$ (2)

$\frac{\∂v}{\∂t}+\frac{\∂\mathrm{uv}}{\∂x}+\frac{\∂\mathrm{vv}}{\∂y}+\frac{\∂wv}{\∂z}=-\mathrm{fu}-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂y}+\frac{1}{\mathrm{\ρ}}(\frac{{\∂\mathrm{\τ}}_{\mathrm{xy}}}{\∂x}+\frac{{\∂\mathrm{\τ}}_{\mathrm{yy}}}{\∂y}+\frac{{\∂\mathrm{\τ}}_{\mathrm{zy}}}{\∂z})$ (3)

$\frac{\∂p}{\∂z}=-\mathrm{\ρg}$ (4)

式中u、v、w分别为速度矢量沿三个坐标轴x、y、z的分量；τ_ij为切应力；f为柯氏力系数，f＝2ωsinψ；ψ为当地纬度；ω为地球自转角速度。

将笛卡儿坐标转换成sigma坐标，即(x，y，z)→(x，y，σ)，可得：连续方程 $\frac{\∂\mathrm{\ξ}}{\∂t}+\frac{\∂\mathrm{hu}}{\∂x}+\frac{\∂\mathrm{hv}}{\∂y}+\frac{\∂\mathrm{\ω}}{\∂\mathrm{\σ}}=0$ (5)

动量方程

$\frac{\∂\mathrm{hu}}{\∂t}+\frac{\∂\mathrm{huu}}{\∂x}+\frac{\∂\mathrm{huv}}{\∂y}+\frac{\∂\mathrm{u\ω}}{\∂\mathrm{\σ}}=\mathrm{hfv}-\mathrm{gh}\frac{\∂\mathrm{\ξ}}{\∂x}+F_x+\frac{\∂}{\∂\mathrm{\σ}}\left(\frac{E_v}{h}\frac{\∂u}{\∂\mathrm{\σ}}\right)$ (6)

$\frac{\∂\mathrm{hv}}{\∂t}+\frac{\∂\mathrm{hvu}}{\∂x}+\frac{\∂\mathrm{hvv}}{\∂y}+\frac{\∂\mathrm{v\ω}}{\∂\mathrm{\σ}}=-\mathrm{hfu}-\mathrm{gh}\frac{\∂\mathrm{\ξ}}{\∂y}+F_y+\frac{\∂}{\∂\mathrm{\σ}}\left(\frac{E_v}{h}\frac{\∂v}{\∂\mathrm{\σ}}\right)$ (7)

式中 $\mathrm{\σ}=\frac{z-\mathrm{\ξ}}{h};$ h为水深；z为水位；ξ为自由表面高度；ω为σ坐标下的垂向流速；水平扩散项 $F_z\≡\frac{\∂}{\∂x}\left({\mathrm{h\τ}}_{\mathrm{xx}}\right)+\frac{\∂}{\∂y}\left(h{\mathrm{\τ}}_{\mathrm{xy}}\right),$ $F_y\≡\frac{\∂}{\∂x}\left({\mathrm{h\τ}}_{\mathrm{xy}}\right)+\frac{\∂}{\∂y}\left({\mathrm{h\τ}}_{\mathrm{yy}}\right);$

E_v为垂向扩散系数。

σ坐标系中的垂向速度ω与z坐标系中的垂向速度ω转换关系如下：

$w=\mathrm{\ω}+u(\mathrm{\σ}\frac{\∂h}{\∂x}+\frac{\∂\mathrm{\ξ}}{\∂x})+v(\mathrm{\σ}\frac{\∂h}{\∂y}+\frac{\∂\mathrm{\ξ}}{\∂y})+\mathrm{\σ}\frac{\∂h}{\∂t}+\frac{\∂\mathrm{\ξ}}{\∂t}$ (8)

式中ω(0)＝ω(-1)＝0。为便于后续建模，将式(5)～式(7)写成如下统一的向

量形式： $\frac{\∂U}{\∂t}+\▿F\left(U\right)=S$ (9)

式中 $U=\left(\begin{array}{c}\mathrm{\ϵ}\\ \mathrm{hu}\\ \mathrm{hv}\end{array}\right);$ $F_x=\left(\begin{array}{c}\mathrm{hu}\\ {\mathrm{hu}}^2+\frac{g}{2}({\mathrm{\ξ}}^2+2\mathrm{\ξ}{h}_s)\\ \mathrm{huv}\end{array}\right);$ $F_y=\left(\begin{array}{c}\mathrm{hv}\\ \mathrm{huv}\\ {\mathrm{hv}}^2+\frac{g}{2}({\mathrm{\ξ}}^2+2\mathrm{\ξ}{h}_s)\end{array}\right)$

$S=\left(\begin{array}{c}-\frac{\∂\mathrm{\ω}}{\∂\mathrm{\σ}}\\ -\mathrm{g\ξ}{S}_{0_z}+F_x+\frac{\∂}{\∂\mathrm{\σ}}\left(\frac{E_v}{h}\frac{\∂u}{\∂\mathrm{\σ}}\right)+\mathrm{hfv}-\frac{\∂\mathrm{u\ω}}{\∂\mathrm{\σ}}\\ -\mathrm{g\ξ}{S}_{0_y}+F_y+\frac{\∂}{\∂\mathrm{\σ}}\left(\frac{E_v}{h}\frac{\∂v}{\∂\mathrm{\σ}}\right)-\mathrm{hfv}-\frac{\∂\mathrm{v\ω}}{\∂\mathrm{\σ}}\end{array}\right).$

该方程组形式中因变量以自由表面高度ξ代替水深h，能够避免静水计算不和谐的问题。

其中，边界条件为：自由表面边界和底部边界应当满足相应的运动学和动力学条件，岸边界法向流速为零；水域开边界则给定水位或流速。

步骤2，超维度算法实现。根据模式分裂方法的思想求解此问题。三维流动的物理过程可以认为是快速移行的外重力波(外模式)和速度垂直剪变所组成。对方程组消掉全部垂直结构，获得外模式的离散方程：

$\frac{d\tilde{U}}{\mathrm{ds}}=-\frac{1}{\mathrm{\ΔV}}\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\tilde{F}}_{n,j}{A}_j+\tilde{S}$ (10)

式中上标～表示垂向平均， ${\tilde{F}}_{n,j}=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{F}_{n,j}/L,$ L为垂向网格数。

各因变量的基本方程为式(5)～式(8)，其中自由表面高度ξ按式(10)计算。采用显格式离散，算法的具体实现过程如下：

1)分层计算单元水平方向的数值净通量；

2)沿垂向平均第一步计算得到的单元水平数值净通量，根据式(10)，可计算出下一时间步上该单元所在的自由表面高程ξⁿ⁺¹；

3)将ξⁿ⁺¹代入式(6)、式(7)，根据各分层的数值净通量求出水平方向的流速uⁿ⁺¹、vⁿ⁺¹；

4)应用中心差分格式离散连续性式(5)，代入ξⁿ⁺¹、uⁿ⁺¹、vⁿ⁺¹即可求出ω，再由式(8)求出物理坐标下的垂向速度wⁿ⁺¹。

步骤3，并行计算。将计算区域分成多个子区域，如首先将大任务划分成子任务，这是任务级的划分，可以分布到计算节点；在任务可以进一步划分的情况下，还可以划分成进程级的划分，同理，映射到计算节点。另外，还可以进行数据的划分，划分之后，便可以进行并发的处理。

每个子区域映射到并行体系结构上的一个计算节点上，结点间的通信采用标准的MPI(Message Passing Interface)。采用计算与通信的重叠并行优化技术：变量关于空间的计算是独立的。在第1时间步，不妨假定计算网格的扫描顺序为：x→y→z......。对于某一方向的计算，可假定为x方向上的计算。

其计算步可分为3步：

1)计算每个处理机上的关于y方向的内边界点(内边界点：指定义在每个处理机t的子区域的边界点，不包括求解区域的自然边界点)；

2)将y方向上的内边界点的数据进行消息传递；

3)计算内点(内点：指定义在每个处理机上的子区域的非边界点)。

步骤4，自适应网格。网格自动生成技术已经在两个领域广泛地应用，一是计算力学中的有限元网格自动生成；二是数字高程模型(DEN)应用中建立地面不规则三角网(TIN)。能实现非结构网格自动生成的方法有很多种，应用较多的主要有狄洛尼方法和推进阵面法。本方法采用的网格自动生成方法总体上属于狄洛尼法，具体地说又属于超三角形渐次插入算法。

自动生成计算网格的主要步骤如下：

1)形成全域超三角形；

2)离散边界，按顺序插入边界点P；

3)搜索包含边界点P的三角形，将此三角形分解为3个三角形；

4)重复2)～3)，直到完成所有边界；

5)删除所有包含一个或多个超三角形顶点的三角形；

6)根据网格尺度控制，加密三角网，直到形成符合尺度的三角形网格；

7)网格优化。

采用这一方法首先可以实现边界自适应(步骤2)～3))，然后可以控制任意点的三角形尺度(步骤6))，为实现地形自适应的网格生成奠定了基础。

步骤5，数据和变化规律比较分析。基于非结构网格有限体积法建立了垂向sigma坐标的超维度河流动力学数值模型，模型采用四边形和三角形混合网格，不仅能够很好地拟合复杂几何边界，而且可根据地形局部加密网格。它建立在守恒型积分方程基础上，具有守恒性，能保证全局水量平衡。Roe平均的黎曼近似解的应用，使格式还具有单调逆风性和高分辨率等特性，可模拟间断型浅水流动。

本发明方法的上述步骤中，步骤1和2可以有效的将超维度的河流信息映射到数学模型和计算模型中，步骤3和4可以应用并行计算和自适应网格快速的进行数据处理，得到精细的模拟仿真结果，步骤5进行所得计算结果分析，从而使得到的结果具有现实意义：超维度河流信息在水利规划上的意义以及自适应并行计算在快速计算上的意义。

在大型水利工程计算和规划中，由于不仅需要一维二维的模拟仿真河流形状，流域，水量及河道网等，还需要三维的模拟仿真河流、河道，尤其更需要通过河流历史的数据，水质数据等进行河流的预测和监控等，这在水利规划和河流流量调配过程中，信息的维度增多可以更好的支持决策和水利规划。但是维度的增多会带来计算量的增大，尤其是在广流域情况下，于是，本发明方法采用自适应并行计算的方法来进行计算，可以实时有效快速的得到结果，这对于现实应用和数字流域的研究是非常有意义的。

本发明方法维度信息丰富，处理速度快，适于大流域处理，其模拟仿真显示效果好，动态且能精确显示；采用本发明方法后可以快速实时地将大型水利工程方案进行精细的模拟仿真，大大降低了规划的风险和成本，而且可以直观的科学的有效的进行决策支持，使得大型水利工程方案的论证进入方便快捷、低成本低风险的时代。在本发明方法支持的系统中，水利工程的模拟修改是非常快捷的，而且可以实时的显示出因此而产生的影响。本发明方法在水利工程的规划设计方面具有重要意义，基于三维的快速处理，可以减轻工程人员的设计负担，并且可以让工程设计进度加快。如未采用并行处理的情况下，设计一个大坝，可能需要1年甚至几年，但是采用本发明方法可以缩短至几个月，并且设计效果更加的精确细致。

应当理解的是，上述针对具体实施例的描述较为详细，并不能因此而理解为对本发明专利保护范围的限制，本发明的专利保护范围应以所附权利要求为准。

Claims (4)

1.一种超维度河流动力学自适应并行监测的方法，其用于一连接并行数据和模型库的并行计算系统，该并行计算系统通讯获取卫星图片和河流水文数据；所述方法包括以下步骤：

A、将超维度数据输入系统，并按照数据所处维度的不同进行分类；

B、建立基于特征型高分辨率数值算法的超维度非结构网格河流动力学模型；

C、依据超维度流体分裂格式的高效并行计算方法，进行维度内的计算和维度间的计算；

D、将计算区域分成多个子区域，每个子区域映射到并行体系结构上的一个计算结点上，结点间的通信采用标准的消息传递接口，在自适应网格中采用计算与通信的重叠并行优化技术，变量关于空间的计算独立。

2.根据权利要求1所述的方法，其特征在于，还包括在所述步骤A之前还设置有步骤：

A0、建立多维度的基本方程模型。

3.根据权利要求1所述的方法，其特征在于，所述步骤C的并行计算过程还包括：

C1、计算每个处理机上的关于某个方向的内边界点，即定义在每个处理机的子区域的边界点，不包括求解区域的自然边界点；

C2、将该方向上的内边界点的数据进行消息传递；

C3、计算内点，即定义在每个处理机上的子区域的非边界点。

4.根据权利要求1所述的方法，其特征在于，所述步骤D中自适应网格的生成过程采用超三角形渐次插入法，包括步骤如下：

D1、形成全域超三角形；

D2、离散边界，按顺序插入边界点；

D3、搜索包含该边界点的三角形，将此三角形分解为3个三角形；重复上述D2～D3步骤，直到完成所有边界；

D4、删除所有包含一个或多个超三角形顶点的三角形；

D5、根据网格尺度控制，加密三角网，直到形成符合尺度的三角形网格；

D6、网格优化。

Images