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）。

海洋模型應該對上述的數值計算速度和精度的限制進行考慮。