2.1　基本条件

20世纪50年代晚期数字计算机的出现迎来了解决实际流体学问题的新时代，从此开始计算能力得到惊人的提高。自1960年（非常原始的IBM704时代）以来，中央处理器（CPU）的速度经历了6个数量级左右的提升，因此在现代计算机例如CRAY T90中，2～3ns的时钟频率接近传统硅片技术的物理极限（约为1ns）。这与CPU指令的并行处理相结合，使个人计算机的计算性能超过每秒十亿次运算（GF，每秒109次浮点运算）的极限而达到顶峰（理论上）。例如，个人Cray T90处理器达到了1.8GF。因此，目前多核CPU计算机已经能保持几十亿次浮点运算的速度，例如16处理器的Cray T90，一个在IBM704上进行计算要花费一个世纪时间的CFD问题在T90上计算所需时间不到1h。同时，计算成本也以相似的方式在下降，现在只需花费几百美元的计算在30年前则是极为昂贵的，需要花费上百万美元。计算机内存的成本也急剧下降，在磁芯存储器中只能存储几千个字，而现在能存储将近1Gw（约为109个字），这使得巨大的实际流体问题能够得到高效解决。计算机速度和内存的增加、计算成本的降低以及算法的发展使在高性能计算机上的大量计算成为可能。而个人工作站的发展也保持了相同的步调，十年前许多需要用高性能计算机进行大量计算的任务现在在现代工作站上就可完成，例如DEC Alpha和硅谷图形台式机，甚至是现代的个人电脑。1955年以来计算机的发展见图2.1.1。

2.1.1　体系结构

20世纪60年代的标量体系结构意味着CPU将会对每次时钟周期执行一次指令，由于CPU的任务本质上可以分为简单的操作，例如将内存位置的内容读取到CPU存储寄存器中、对寄存器中存储的单词进行算术运算、将寄存器的内容重新存储到内存中等，对一个数组的每个元素进行的算术运算将会花费几个时钟周期的时间。然而，矢量（或者管线式）CPU，例如20世纪70年代的Cray-1s，每个时钟周期能够同时执行几个指令，从而可以一次对数组的多个元素进行操作，因此，相同的时钟速度，矢量处理器的计算速度是标量处理器的几倍，当不包含条件语句时，例如使CPU每次对数组的多个元素进行运算的语句，以及主要由数组运算组成的计算的语句，那么计算速度将会变得更快。然而，CPU速度的物理限制仍然是一个瓶颈，目前的性能水平仍然在GF范围内，即使是最快的矢量CPU也如此。使用多个CPU同时对一个问题进行处理是克服该限制的唯一方式，现代矢量处理计算机例如CRAY T90就有几个CPU（该情形下为16个），执行多重任务的能力使持续运算速度达到几GF。

要突破TF运算速度障碍（百万兆次，1012次浮点运算）需要具有上百个甚至上千个速度极快的矢量CPU的大型并行计算机。这是21世纪高性能计算机的选择，所使用的模式可能为SIMD（单指令多数据）或MIMD（多指令多数据）。在SIMD计算机中，对同一数组的多个元素同步执行相同的运算，因此在这一标量运算期间，大多数CPU是闲置的。在MIMD计算机中，每个CPU对CFD问题的不同部分进行处理，因此将任务分给不同CPU（负载均衡），而计算方法的序列处理特性可能会要求一个CPU等待另一个CPU的运算结果（处理器间交流），这是重要编程任务所需要做的（光谱元素模型对并行处理环境的自适应的讨论可参考Haidvogel等，1997）。对于具有分布式（而不是共享）内存的计算机而言，这也意味着处理器间需要具有一个很大的交流带宽，从而使数据及计算结果可以被很快地按需要从一个CPU传送到另一个CPU。因此，N个CPU的并行计算速度增加的影响因素远不止是N这个因素，它很大程度还取决于CFD问题的性质、计算机结构体系以及并行指令的效率等。然而，对于真实海洋和气候模型等非常巨大的“挑战”型计算问题来说，大型并行计算机是唯一可行的解决办法。目前的CPU速度在慢慢接近1ns（1kMHz），Cray T3E-1200的时钟速度达到600 M Hz（兆赫），单处理器的最佳性能为1.2GF。对于有1000个CPU的计算机来说，将会得到1TF的最佳性能。目前最高性能的计算机是SGI ASCI（加速策略计算计划）计算机，由5040个500MHz的处理器构成，达到的最佳性能是2.52TF，这一计算机在一些极其巨大的问题的计算中具有达到持续计算水平超过1TF的潜能（Dongarra，1999；高性能计算机概述可参考Van der Steen，1999）。

联邦编程例如加速策略计算计划（ASCI）的目标是在20、21世纪之交建成3-4TF的计算机，并在2003年左右达到10TF的计算能力。另一方面，快速CPU（比如奔腾3）是个人计算机行业中的主力，它们被组合在一起从而生产具有多GF能力的大型并行计算机，这是十年前高性能计算机从未达到的水平。现代工作站例如DEC Alpha的计算能力虽然只能达到Cray 90 CPU的一小部分，但它不仅使海洋模型中复杂的前处理和后处理能够实现，而且也能解决几年前需要高性能计算机的CFD问题。虽然发展很迅速，但是要在几年的时间里对性能进行提升很困难。

2.1.2　计算误差

数值计算包含两类误差：舍入误差和截断误差。前者由计算机不能用无限精度表示一个浮点数造成，这种浮点数不可避免地要进行舍入，仅能对数学运算进行近似表示，因此，任何计算机上的整数运算都是精确的，而浮点数运算的精度取决于所用的字长以及计算机中对该数的表示方式——小数部分和指数分配多少字长。在过去的32位（字长为4字节）计算机中，只能用大约7位有效小数对一个数进行表示，而现代CFD工作是在64位计算机上进行，例如CRAY T90、CRAY T3E以及DEC Alphas，它们可以用大约15位有效小数对一个数进行表示。计算机精度是一个浮点运算在计算机上计算时能够达到的精度，定义为可以与1.0相加但是得到的结果不等于1的最小数。对于32位计算机而言，这一精度约为10-8，而64位计算机则为10-16左右。然而，具有任意基本精度的计算机上可以进行双精度运算（和任意精度的运算），利用软件指令可以使计算中的误差显著降低。所有情况下，数字计算机上进行的任何数值运算都避免不了要产生舍入误差，在对操作数进行的序列运算中，这种误差会积累起来，而很难对这种积累的舍入误差的大小进行控制，除非使用使运算速度降低的昂贵高精度机器运算。许多CFD运算要求64位的精度，例如对多年间水体块的盐度进行追踪的海洋模型。

现代CFD计算机更大的字长意味着它可以支持数字的更大动态范围。例如当ANSI/IEEE 754的标准32位浮点数的范围为±1.2×10-38～±3.4×10+38（在这个范围之外的任何数都会产生下溢和上溢负载），64位浮点数则能到±2.2×10-308～±1.8×10+308的范围。在Cray矢量处理器上，则大约能达到10-2467～102465的范围。

截断误差由数字计算机进行计算时的离散化需要产生。空间中连续的变量只能在数字计算机上用预先选定的空间中的离散点进行表达，一个变量在该离散点的寻常导数和偏导数必须用该变量在该格网点及其邻域处的值进行表达，这可以通过泰勒级数展开式得到。函数y（x）在邻域格网点j+1处（x=x0+h，h为格网间距）的值可以用一个包含y（x）中的项的无穷级数表示，它在j点（x=x0）处的导数为

这样的展开可以用来求函数本身在格网点值处的任意阶导数，此过程中要对展开进行恰当的截断，比如，在格网点j点（x=x0）处的一阶导数可以写为

用变量在两个邻近格网点的值之间的有限差分对导数进行表示，ε是截断误差，在此情形中它的阶为h，因此这一有限差分逼近被称为一阶精度形式，n阶精度形式的截断误差为o（hn），任何精度下的有限差分逼近以及任意阶导数可以利用邻近格网点处的泰勒级数展开式得到。高阶导数和高阶精度形式自然需要更多的格网点。例如，一阶导数的二阶形式要求使用3个连续格网点处的值，可写为

由于对称性，yj没有显式出现。类似地，二阶导数的二阶形式可以用下式得到

同样地，由于对称性二阶导数的二阶形式只使用3个点，而相同精度下不对称的形式则需要4个格网点。对称的有限差分被称为中央差分，例如式（2.1.3）和式（2.1.4）。式（2.1.2）是非对称的，称之为前向差分（或顺风形式，假设流的方向沿x增大的方向），相反，一阶精度形式是后向差分（或逆风形式）为

高阶精度形式可以利用更多的邻近格网点得到，例如一阶导数的中央四阶的有限差分形式为

两个或三个变量函数的偏导数的类似表达式可以便用对应的泰勒级数展开式表示。

尽管目前的计算机具有无限的精度（即舍入误差为0），截断误差也是积累的，截断误差使得包含偏微分方程和常微分方程的CFD计算只是近似的。在理论上，截断误差是建模者可以控制的，因为格网大小的减小（相当于增加模型分辨率）或者高阶形式的求取使得截断误差降低。通过令h→0，理论上任意阶形式的截断误差都可以被消除，但实际上，由于要使模型的分辨率翻倍，就会使计算时间增加一个数量级，内存需求也成倍增长，因此通过增加分辨率的方式降低截断误差通常是不切实际的。尽管可以使h变得任意小而使截断误差最小化，在精度有限的计算机上，h的减小意味着舍入误差的增大，因为这种情况下相同的计算将对应地包含更大量的运算。因此，由于任何CFD计算中的总误差是舍入误差和截断误差之和，通常情况下，当h减小时，总误差由于截断误差更为迅速地降低而开始减小，到达某个点后，h的进一步减小会引起总误差的增大，因为舍入误差的增大开始超过截断误差的减小。在大多数情况下，模型分辨率的选择取决于可用的计算资源，并且通常情况下，即使只用两种不同的分辨率对计算进行收敛试验也是不实际的。

高阶精度形式通常是不增加模型分辨率情况下减小截断误差的一种方式，但是很少看到CFD的计算超过四阶精度的空间差分形式。大多数CFD计算使用二阶精度的空间差分，因为这种情况下边界条件的指定比四阶精度形式下简单得多，大多数海洋模型也使用一阶或二阶空间差分进行表示，首选的时间差分形式多为一阶或二阶精度。如果能够得到低截断误差，那么格网分辨率的增加将会变得简单，而且能够承受额外的计算负担，而不需要使用高阶差分对已有的复杂代码用新的形式进行表达。然而，用高阶空间差分（Bender，1996；Haidvogel和Beckmann，1977）和时间差分形式（Durran，1991；Holland等，1998）对这些模型进行重新表达的研究受到越来越多的关注。

在减小截断误差并提高计算精度的过程中，是选择使用特定格网分辨率下的高阶形式还是使用高格网分辨率下的低阶形式，这需要对增加的复杂度和增加的计算需求这两方面进行折中。显然，动力学的考虑决定着所需要的最小格网分辨率。高阶形式不能够弥补丢失的动力学和补偿低分辨率带来的不足。然而，假设所选的格网能描述动力学过程，则高阶形式通常能更好表达。通过增加格网分辨率使截断误差得到减小意味着更大的计算机内存和CPU时间需求，而高阶精度形式仅会增加CPU时间需求。对于矢量处理器的体系结构来说，CUP时间需求通常不是精度形式本身精度的强函数。这些都要求采用更高精度，这种阶数只受到侧边界条件描述的复杂度限制。

海洋模型中阶数和分辨率的详细分析可以参考Sanderson（1998）。简单地说，大多数情况下阶数为n的差分形式（截断误差为ε=aΔn，其中Δ是格网大小，a是比例常数，它取决于问题和差分形式的特征以及所选的模型域）的计算成本可写为c=bnmΔ-d，其中d为问题的维数（对时间相关的三维海洋过程来说通常为4），b是与算法相关的一个比例常数，对于所有优化实现来说，系数m接近于1。为了简单起见，忽略高分辨率模型中增加的舍入误差的影响，对于期望截断误差ε，阶数为n1和n2的两个差分形式的相对计算成本为

显然，如果忽略边界条件增加的复杂度（这在一些情形中会产生重编码的负担）的话，对计算复杂度而言，高阶精度形式是比较可取的。因为截断误差和计算成本的乘积与nmΔ-d成正比，可以看成当n＜d时增大m的值比减小Δ使截断误差降低的速度快很多。Sanderson（1998）发现，在大多数情况下，如果n＜d，格网细化将不是降低计算误差的好方法。因此，阶数大于2的精度形式对海洋模拟比较有用，按照推论，实际限制也许是d+1，也就是说，对海洋问题的处理中，对大于四阶的差分进行求取不是有利的，即使这样，仅对具有高阶形式的控制方程中最重要的项进行差分会变得更简单，比如压力梯度项（McCalpin，1994）和科里奥利项，这是海洋模型中广泛使用三阶时间差分（Durran，1991）和四阶空间差分形式的主要原因（Haidvogel和Beckmann，1997）。

海洋模型应该对上述的数值计算速度和精度的限制进行考虑。