이중 푸리에 급수 분광법 역학코어의 정확도와 계산 효율성 평가

1. 서 론

집중호우 등 큰 피해를 끼치는 기상 현상들이 수백 킬로미터 이하의 중규모에서 발생하기 때문에 중규모 현상 모의 정확도를 향상시킬 수 있는 고해상도 수치예보모델링은 정확한 일기 예보를 위해 매우 중요하다. 최근 컴퓨팅 파워의 비약적인 발달로 예전에는 실행 불가능하던 해상도의 모델 적분이 가능해지자, 일기 예보뿐만 아니라 기후 예측에서도 고해상도 모델링의 중요성이 대두되고 있다.

수치 모델의 해상도를 높이면서 필연적으로 발생하는 문제점들 중의 하나는 해상도에 따라 기하급수적으로 늘어나는 연산량 증가이다. 이를 극복하기 위해 물리모수화의 최적화와 함께 다양한 역학코어(dynamical core)에 대한 연구가 진행되고 있다(Williamson, 2007; Choi and Hong, 2016; Ullrich et al., 2017; Diamantakis and Váň, 2022). 전구 수치모델의 경우, 구면 조화함수(spherical harmonics function, SPH)를 기저함수(basis function)로 사용하는 분광 역학코어가 광범위하게 사용되고 있다. 높은 안정성과 정확성을 가지기 때문이다. 그러나 SPH는 기저함수를 분해하고 또 다시 역변환하는 과정에 많은 계산 시간이 소요되는 단점을 가지고 있다. SPH는 동서방향으로는 푸리에 급수(Fourier series)를, 남북방향으로는 버금 르장드르 다항식(associated Legendre polynomial)을 활용한다. 이 중 시간이 많이 소요되는 부분은 르장드르 다항식 변환 과정이다. 변환할 방향으로의 격자점의 개수를 N이라고 했을 때, 르장드르 다항식으로의 변환은 O(N²) 만큼의 연산을 요구하며, 모델의 수평 해상도가 늘어날수록 계산량이 기하급수적으로 늘어나는 특징이 있다. 이런 단점을 보완하기 위해 르장드르 다항식 변환 속도를 신속하게 하는 방법들이 개발되기도 했다(Suda, 2005; Tygert, 2008).

SPH의 계산 속도를 개선하기 위한 한가지 방법으로 이중 푸리에 급수(double Fourier series, DFS)를 기저함수로 사용하는 분광 역학코어가 제안되었다. 구면조화함수의 동서방향 파동함수이기도 한 푸리에 급수는 변환 시 고속 푸리에 변환(fast Fourier transform, FFT) (Cooley and Tukey, 1965)이라 불리는 알고리즘을 적용할 수 있어서 소요되는 계산량이 O(N·logN) 수준이기 때문에 르장드르 다항식 변환보다 훨씬 더 빠르게 계산할 수 있다. 이를 토대로 동서방향뿐만 아니라 남북방향으로도 푸리에 급수를 적용하는 DFS가 고안되었다(Merilees, 1973; Orszag, 1974). 비록 DFS가 SPH에 비해 구면 위의 모든 곳에서 일정한 해상도를 가지지 못하며, 남북방향 미분 혹은 라플라시안(Laplacian) 연산이 더 복잡하다는 단점을 가지고 있지만, 계산 속도가 더 빠르기 때문에 이러한 문제들을 극복한다면 현업 수치모델에 유용하게 사용될 수 있다(Cheong, 2006; Koo and Hong, 2013; Park et al., 2013; Yoshimura, 2022).

구면에서 DFS를 사용하는 것에 대한 이론적 정당성은 1970년대 초기 연구에서 이미 확립되었다(Merilees, 1973; Orszag, 1974). 이들은 남북방향 기저함수를 동서방향 파수에 따라 신중하게 설정함에 따라 구면 위에서 불연속점이 없이 DFS가 작동할 수 있음을 보였다. 일례로 구면에서 포아송 방정식(Poisson equation)을 푸는 데 DFS가 적용되었다(Yee, 1981). 그러나 DFS는 시간 의존적 문제(time dependent problem)에는 거의 사용되지 않았다. 시간 의존적 방정식의 수치적 풀이에 DFS가 사용될 경우, SPH에 비해 모델 불안정성이 크게 증가하였기 때문이다(Boyd, 1978). 시간 의존적 문제는 최근에야 본격적으로 다루어지기 시작했다(Cheong 2000a; Cheong 2000b; Cheong et al., 2002; Cheong et al., 2004). 일례로 Cheong (2000a, b; 이하 C2000)는 새로운 형태의 DFS 기저함수를 소개하고, 해당 기저함수에 기반한 DFS를 이용하여 포아송 방정식뿐만 아니라 천수 방정식(shallow water equation)의 수치해를 구할 수 있음을 보였다. 특히, SPH 대비 결과에는 큰 차이가 없지만 계산 시간은 크게 줄어듦을 보고하였다. Cheong (2006)은 3차원 DFS를 개발하면서, 연직으로 빠르게 전파하는 중력파(gravity wave)에 의한 불안정 효과를 완화할 수 있는 준암시적 시간 적분(semi-implicit time integration) (Hoskins and Simmons, 1975) 기법을 제안하였다. 이 방법은 Global/Regional Integrated Model system (GRIMs) (Hong et al., 2013; Koo et al., 2023)에 최초 적용되었으며(Park et al., 2013), GRIMs-SPH와 GRIMs-DFS는 일기 예보 및 계절 모의에 있어서 유사한 성능이 확인되었다(Koo and Hong, 2013; Park et al., 2013).

최근 Yoshimura (2022, Y2022)는 새로운 방식의 DFS 남북방향 기저함수 계산법을 제시하였다. 이 방법은 기존 C2000에서 제시된 방법과 비교하였을 때 DFS에 대한 초기의 이론적 연구(Orszag, 1974; Boyd, 1978)에서 도출된 계산 안정화 조건에 더 잘 부합한다는 특징을 가지고 있지만, 추가적인 계산비용 증가가 발생한다는 단점 역시 가지고 있다. Y2022는 새로운 방법을 단일 층의 구면 천수 방정식 풀이에 적용해 봄으로써 기존 대비 비교적 개선된 안정도를 가지고 잘 작동함을 보였으나 실제 3차원 대기모델로 확장되지는 않았다.

본 연구는 고해상도 전구 수치모델링에서 DFS의 장점을 이해하기 위해 C2000과 Y2022를 GRIMs에 적용하고 기존 GRIMs-SPH (Sela, 1980; Kanamitsu, 1989)와 비교‧분석하였다. C2000은 GRIMs에 장착되어 있는 GRIMs-DFS를 사용했다. 반면, Y2022는 기존 2차원 모델을 3차원으로 확장하여 GRIMs에 이식하였다. SPH 대비 DFS의 정확도 및 계산 효율성은 두 가지 이상화 실험(idealized test)을 통해 평가하였다. 첫 번째 실험은 물질 이동 실험이며, 두 번째 실험은 경압 파동 발달 실험이다.

본 논문은 다음과 같이 구성되었다. 2장에서는 C2000과 Y2022에 소개된 DFS 변환 방법을 각각 소개하고 그 특성을 비교하였다. 3장에서는 본 연구에서 사용한 두 가지 역학코어 실험의 설계에 대해 설명하고, 4장에서는 두 실험의 결과를 제시하였다. 마지막으로 5장에서는 주요 결과를 요약하고, 이를 토대로 DFS 방법의 활용성에 대해 간단히 논의하였다.

2. DFS 방법 소개 및 비교

분광 역학코어는 구면 위의 위도, 경도 좌표계에서 정의된 함수 F(θ, ϕ)를 다음과 같이 분해하여 사용한다.

이때, θ는 여위도(90^o - 위도), ϕ는 경도, m은 동서방향 파수, n은 남북방향 파수, N은 파수 절단(truncated wavenumber truncation), F_n,m은 분광계수, 그리고 $$ P_n^{\left|m\right|}$$은 버금 르장드르 다항식을 의미한다. 즉, SPH는 남북방향의 기저함수로 버금 르장드르 다항식($$ P_n^{\left|m\right|}$$)을 사용한다. 반면, DFS는 다음과 같이 삼각함수를 사용한다(C2000).

한가지 눈여겨볼 점은 동서방향 파수 m에 따라 남북방향 기저함수 F_m을 다르게 적용한다는 것이다. 이는 SPH의 남북방향 기저함수가 가지는 대칭성을 보존하고, 극(θ = 0^o, 180^o)에서 특이성(singularity)이 생기지 않게 하기 위함이다. 예를 들어, m이 0이 아닐 때 극에서 F_m이 0이 아닌 값을 갖는다면 극 주변에서 불연속이 발생할 것이다.

DFS 방법의 실질적인 구현을 위해서는 각 기저함수에 곱해지는 분광 계수, 즉 F_n,m을 구할 수 있어야 하는데, 이러한 삼각함수 형태 급수들의 계수는 FFT를 기반으로 한 이산 푸리에 변환 알고리즘들을 이용하여 구할 수 있다. m이 2이상의 짝수 일 때는 이를 바로 적용할 수는 없으며, F_m(θ)에 sinθ를 나눈 뒤, 사인 변환을 시행하는 식으로 계수를 구한다(C2000). 한편, Y2022는 sinθ sinnθ의 계수 F_n,m을 계산하는 다른 방법을 제시하였다. 이 방법의 핵심은 원함수 F_m(θ)와 결과로 나온 sinθ sinnθ 급수 간의 차이의 제곱 적분 값, 즉 L2 오차를 최소화하는 것이다.

E를 최소화하는 계수 F_n,m을 구하기 위하여 최소 제곱법(least squares method)이 사용되며, 결과적으로는 원함수를 코사인 변환한 다음 삼중 대각행렬(tridiagonal matrix)을 역변환하는 계산 과정을 거쳐 얻을 수 있다. 이에 관한 상세한 유도 과정과 계산 과정은 모두 Y2022에 기술되어 있다. Y2022에서 제안한 방법은 삼중 대각행렬과 관련된 연산이 추가되어 계산량이 증가하지만 그 크기가 O(N) 정도이므로 전체 계산량(SPH: O(N²), DFS: O(N·logN); 제 1장 참고)과 비교하였을 때 상대적으로 소량이다. 따라서 여전히 SPH와의 비교에서는 Y2022의 DFS가 시간 상 이점이 유효하다.

본 논문에서 C2000과 같이 sinθ를 나눠서 계수를 구하는 방식을 DFS_1, Y2022와 같이 최소 제곱법을 이용하여 계수를 구하는 방식을 DFS_2로 명명하였다. Figure 1은 DFS_1과 DFS_2 변환 방법을 요약하여 나타낸 모식도이다. 참고로, Y2022가 제안한 방법은 본 연구에서 기술한 방법 이외에 남북방향 기저함수의 세분화 등 추가적인 과정이 포함되어 있다. 그러나, 가장 핵심적인 아이디어인 변환 과정에서 최소 제곱 최적화 과정은 동일하다. 기저함수의 세분화와 관련한 추가적인 정보는 Orszag (1974)에서 찾아볼 수 있다. 한편, 식(3)과 같은 형태의 기저함수가 구면 위의 변수를 나타내기에 충분한다는 근거는 C2000에 서술되어 있으며, 본 연구에서는 특히 최소 제곱 최적화의 효용성을 알아보기 위해 DFS_2에서 DFS_1과 같은 기저함수를 사용하였다.

Fig. 1.

Schematic diagram for the discretized Fourier transformation process from physical space to wave space in the (a) DFS_1 and (b) DFS_2. DFT denotes discrete Fourier transform, DST denotes discrete sine transform, and DCT denotes discrete cosine transform.

DFS_1과 DFS_2 변환 방법의 특성을 파악하기 위해 각각의 방법이 임의의 원함수를 얼마나 잘 표현하는지 비교하였다(Fig. 2). 사용한 원함수는 형태가 매끄럽고(미분 가능) 양극에서 함수 값이 0으로 수렴한 경우(Fig. 2a), 형태가 매끄럽지만 양극에서 0이 아닌 값을 가진 경우(Fig. 2b), 양극에서 함수 값이 0이지만 함수의 형태가 매끄럽지 않은 경우(Fig. 2c)이다. 각각의 변환 방법으로 sinθ sinnθ 급수로 분해한 다음, 다시 그 계수를 사용해 원함수를 재구축(reconstruction)한 결과를 빨간색과 파란색으로 나타냈다. 여기서 급수 분해를 위해 사용된 격자의 개수는 [0, π] 안에서 총 24개이다. 원함수가 급수로서 재구축된 결과를 보면 Fig. 2a에서는 DFS_1과 DFS_2 방법 모두 급수가 원함수를 잘 표현하고 있음을 볼 수 있다. 반면 Fig. 2b에서는 DFS_2는 비교적 급수가 원함수를 잘 표현했지만, DFS_1로 구한 급수는 그에 비해 훨씬 더 큰 진동을 보였다. 이는 원함수를 sinθ로 나누는 과정에서 sinθ가 0에 매우 가까워지는 극 지점 주변에서 값이 0으로부터 많이 떨어지게 되는데, 이를 사인 변환하면서 생기는 오차로 분석된다. 마지막으로 Fig. 2c에서는 양 방법이 모두 원함수의 불연속성의 영향으로 깁스 현상(Gibbs phenomenon)이 뚜렷하게 나타났다. 특히 이 경우에서, 함수 값이 1인 구간에서는 DFS_2가 적은 오차를 보이고 함수 값이 0인 구간에서는 DFS_1이 적은 오차를 보이는 등 두 방법 간의 명확한 우열을 가릴 수 없었다.

Fig. 2.

Comparison of reconstructed functions by the DFS_1 (red) and DFS_2 (blue) runs for the shapes of (a) normal cosine function, (b) trigonometric function with non-zero boundary, and (c) step function. The black dots represent the original function ranged in [0, π] with 24 gird points.

종합하면 원함수가 양극에서 0이 아닌 값을 가지는 경우에는 Y2022에서 제안된 DFS_2가 원함수를 표현하는 데에 있어서 더 고성능을 지니고 있다고 할 수 있다. 반면, 함수가 매끄럽지 못한 형태를 가지고 있는 경우에는 양 방법이 모두 구조적인 오차를 보였다. 그러나 실제 대기모델에서 남북방향으로 변환해야 할 함수가 Figs. 2b, c와 같이 극단적인 형태를 지닐 확률이 크지는 않기 때문에 위와 같은 분석 만으로 DFS_2의 성능을 속단할 수는 없다.

3. 실험 설계

이상 조건에서 3차원 모의 성능을 비교하기 위하여 Dynamical Core Model Intercomparison Project(DCMIP) (Ullrich et al., 2012)에서 제시된 테스트 베드 실험 중 두가지를 진행하였다. DCMIP은 대기 대순환 모형(general circulation model, GCM)의 역학 방정식 처리 성능을 검증하기 위해 마련된 프로젝트이다. Ullrich et al. (2012)에는 이를 위한 여러 실험들의 구성이 기술되어 있어, 각 실험의 초기 조건, 실험 옵션, 모델 적분 기간, 실험의 의도 등을 확인할 수 있다. 본 연구에서는 수평방향 움직임을 잘 알아볼 수 있는 검증사례를 선정하여 비교 실험을 진행하였다.

4. 실험 결과

4.1 실험1

Figure 3은 실험1 결과를 보여준다. 사용된 모델 해상도는 T126L30로, 동서방향 파수 절단은 126이며 연직으로 30층을 가진다. Figures 3a-c는 모델 적분 6일째의 결과이며, 바람장의 변형(deformation)이 가장 강한 시점이다. 반면 Figs. 3d-f는 모델 적분 12일째로 모든 물질이 이론상 원래의 위치로 돌아오는 시점이다. 전체적인 결과를 살펴보면 DFS_1과 DFS_2 모두 물질의 농도를 전반적으로 잘 모의하고, SPH와 비교하였을 때 큰 차이가 없다. 12일째에 대부분의 물질이 원래의 위치로 잘 돌아오며, 바람의 뒤틀림이 가장 심한 시점인 6일째에도 모든 역학코어가 올바른 형태를 보이고 있다(Lauritzen and Thuburn, 2012).

Fig. 3.

Tracer field at day (a-c) 6 and (d-f) 12 for the idealized 3D deformational flow, simulated by the (a, d) SPH, (b, e) DFS_1, and (c, f) DFS_2 runs at T126L30 resolution.

Table 1은 12일째에 물질의 농도 분포를 초기 조건과 정량적으로 비교한 수치를 보여준다. 표시된 정규화된 L₁, L₂, L_∞ 오류(normalized L₁, L₂, L_∞ error)는 다음과 같은 방식으로 계산되었다.

$\[ L_1\left(q\right)=\frac{\mathrm{\Sigma }\left|p_i-q_i\right|}{\mathrm{\Sigma }\left|q_i\right|} \]$

(5)

$\[ L_2\left(q\right)=\sqrt{\frac{\mathrm{\Sigma }{\left|p_i-q_i\right|}^2}{\mathrm{\Sigma }{\left|q_i\right|}^2}} \]$

(6)

$\[ L_{\mathrm{\infty }}\left(q\right)=\frac{\text{max}\left|p_i-q_i\right|}{\text{max}\left|q_i\right|} \]$

(7)

Table 1.

Normalized error norms for the idealized tracer field at day 12, simulated from the SPH, DFS_1, and DFS_2 runs at T126L30 resolution.

여기서 p_i는 각 점에서 초기 조건에 해당하는 물질의 농도, q_i는 각 점에서 12일째 물질의 농도를 나타낸다. Figure 3에서 각 역학코어 간의 차이가 주목할 만큼 크지 않았던 것처럼, 정규화된 오류값도 차이가 크지 않다(Table 1). SPH와 비교하여 DFS는 더 높은 L₁, L_∞ 오류를 보이는 반면 L₂ 오류는 더 낮게 나오는 등 모의 성능 차이가 뚜렷하게 나타나지 않는다. 한편, DFS_1과 DFS_2의 비교에서는 세 가지의 정규화된 오류 값이 소수점 셋째 자리까지 일치하여, 그 차이가 미미했다.

한가지 차이점은 물질의 농도가 0 이하로 내려가는 영역의 비율이다. 이는 물질이 바람에 의한 이류로 수송되는 과정을 분광법으로 모의함에 있어서 가장 민감한 문제 중 하나이다. 분광법을 사용하면 파수 절단에 의해 어쩔 수 없이 원함수에 없는 진동, 즉 깁스 현상이 생기게 되는데, 물질의 농도를 모의하는 경우에는 이것에 의해 농도가 0 이하로 떨어지는 물리적으로 타당하지 않은 결과가 나오기도 한다. 이 실험에서는 DFS_2가 DFS_1에 비해 물질 농도가 0 미만인 영역이 작게 나왔다(Table 2). SPH와 DFS의 비교에서는 6일째에서 SPH의 음수 영역 비율이 더 높고, 12일째에서 DFS_1과 DFS_2 사이의 비율을 보이는 등 명확한 우열을 가릴 수 없었다.

Table 2.

Same as Table 1 but ratio of the area with negative values to the total area at day 6 and 12.

4.2 실험2

Figure 4는 실험2의 결과를 보여준다. 사용된 모델의 해상도는 T126L30이며, Fig. 4에 표시된 결과는 모델 적분 기간 10일째의 여러 가지 변수들의 분포이다. Figures 4a-c는 각각 SPH, DFS_1, 그리고 DFS_2를 이용하여 모의한 850 hPa 기온의 동서방향 에디($$ T_{850}^{\mathrm{\text{'}}}$$)이다. 여기서 동서방향 에디는 해당 시점의 각 위도에서 그 위도의 기온 값의 평균을 제거하여 만들어졌다. Figures 4d-f는 850 hPa 남북방향 바람 속도(V₈₅₀), 그리고 Figs. 4g-i는 지면 기압(P_sfc) 모의 결과를 나타낸다. 세 실험에서 모두 (235^oE, 65^oN)을 중심으로 저기압성 소용돌이가 매우 크게 발달한 양상을 볼 수 있으며, 그로부터 서쪽으로 가면서 또 다른 저기압성 소용돌이가 발달하고 있는 등 전형적인 경압 파동의 모습이 잘 나타난다. (235^oE, 65^oN)에 위치한 가장 큰 저기압 소용돌이는 DFS에서 좀 더 강하게 모의되며, (175^oE, 55^oN)에 위치한 세 번째 저기압 소용돌이는 SPH에서 좀 더 강하게 모의되었다. 또한, 저기압 소용돌이에 연관된 남북방향 바람 속도 역시 SPH와 DFS를 비교하면 사뭇 다른 결과를 보이고 있다. 하지만 DFS_1과 DFS_2 간의 차이는 SPH와 DFS 간의 차이에 비해 매우 미미하다. DFS_1과 DFS_2 둘 간의 차를 구하면 최대 각 변수의 스케일의 10^-3배 정도의 오차가 있으며, 오차가 생기는 지점 역시 무작위로 분포해 있다. 즉, 해상도와 관계없이 DFS_1과 DFS_2 사이의 오차는 DFS가 SPH에 대해 보이는 편차에 비하면 거의 의미가 없는 수준임을 알 수 있다.

Fig. 4.

(a-c) 850-hPa temperature zonal eddy (K), (d-f) 850-hPa meridional wind speed (m s-1), and (g-h) surface pressure (hPa) at day 10 for the dry baroclinic wave development, simulated from the (a, d, g) SPH, (b, e, h) DFS_1, and (c, f, i) DFS_2 runs at T126L30 resolution.

Figure 5는 Fig. 4와 같은 검증사례이지만, T510L30의 상대적으로 높은 수평 해상도를 이용한 실험 결과이다. Figure 5에서도 역시 전형적인 경압 파동의 형태가 잘 나타나는데, 저해상도 실험인 Fig. 4와 비교하였을 때 SPH를 이용한 결과와 DFS를 이용한 결과 사이의 유사성이 상당히 증가하였다. DFS_1과 DFS_2는 Fig. 4와 마찬가지로 거의 같은 결과를 보인다. 결과적으로 세 역학코어의 결과 간의 차이는 그렇게 크지 않으며, 바람 속도 0m s^-1 선, 지면 기압 1000 hPa 선 등 매우 작은 신호에도 민감하게 반응하는 영역을 제외하고는 거의 비슷한 결과를 보인다. 저해상도 실험에서 확인할 수 있는 SPH와 DFS 간의 차이가 고해상도 실험에서는 상당 부분 제거되며, DFS_1과 DFS_2 간의 차이는 해상도에 관계없이 크지 않아 SPH와 DFS 간의 차이에 비해서 중요도가 적음을 볼 수 있다.

Fig. 5.

Same as Fig. 4 but at T510L30 resolution.

Figure 6은 실험2 검증사례를 보다 자세히 분석하기 위해 북위 60도 선에서 해상도 별로 세 역학코어 간의 실험 결과를 비교한 그림이다. Figures 4, 5와 마찬가지로 초기 상태로부터 10일 이후의 값을 표시한 것이지만, 고정된 하나의 위도에서 T126L30, T254L30, 그리고 T510L30 해상도에서의 결과도 같이 보여주고 있다. Figures 6a-c는 각각 T126L30 실험에서 나온 60 N에서의 850 hPa 남북방향 바람 속도, 850 hPa 지위 고도, 850 hPa 기온이며, Figs. 6d-f는 T254L30 실험의 동일한 변수, Figs. 6g-i는 T510L30 실험의 동일한 변수이다. 녹색 실선이 SPH의 결과이고, 적색 점선이 DFS_1, 청색 실선이 DFS_2의 결과인데, DFS_1과 DFS_2의 결과가 거의 차이가 없어 적색선과 청색선이 겹쳐서 보인다. SPH와 DFS 결과를 비교해보면, 모델 해상도에 따라 SPH와 DFS 사이에 차이가 있음을 알 수 있다. 충분한 고해상도에서는 변수에 관계없이 SPH와 DFS 간 차이가 거의 없지만, 수평해상도가 T126인 실험에서는 각 변수의 국소 피크의 강도가 과소모의 되는 등 차이를 보였다. 한편, 그에 비교하였을 때 DFS_1과 DFS_2의 차이는 거의 없었다. 이는 Fig. 2a 실험과 어느 정도 유사성을 보인다. 실험2에서 주어진 변수들도 전체 영역에서 충분히 매끄럽고 특이성이 없으며, 극 주변에서도 유의미한 신호가 있지 않기 때문에 Fig. 2a와 같이 DFS_1과 DFS_2가 거의 같은 결과를 보였다.

Fig. 6.

(a, d, g) 850-hPa meridional wind speed (m s-1), (b, e, h) 850-hPa geopotential height (m), and (c, f, i) 850-hPa temperature (K) at day 10 for the dry baroclinic wave development, simulated from the SPH (green), DFS_1 (red dotted), and DFS_2 (blue) at the resolutions of (a-c) T126, (d-f) T254, and (g-i) T510 with the 30 vertical levels.

Table 3은 실험2의 초기 상태로부터 10일 이후에 850 hPa 남북방향 바람 속도, 지위 고도, 기온이 전체 영역에서 가지는 최대·최솟값을 보여준다. 각 해상도와 변수에 대해 첫 번째 줄이 SPH , 두 번째 줄이 DFS_1, 세 번째 줄이 DFS_2를 이용한 결과를 가리키고 있다. Figure 6에서 확인한 바와 같이 DFS_1과 DFS_2는 변수의 종류에 관계 없이 거의 같은 결과를 보이며, 값의 차이는 해당 변수가 가지는 스케일의 약 10^-3배 이하이다. 한편, SPH와 DFS의 비교에서는 비교적 큰 차이가 나타나며, 특히 변수 중에서는 바람 속도에서 확연한 차이를 보인다. Figure 6과 마찬가지로 저해상도 실험에서 고해상도 실험보다 SPH와 DFS 간의 차이가 더 크게 나타나는 경향을 확인할 수 있다.

Table 3.

The global maximum (max) and minimum (min) of 850-hPa meridional wind speed (V850; m s-1), 850-hPa geopotential height (Z850; m), and 850-hPa temperature (T850; K) for the dry baroclinic wave development at day 10, simulated from the SPH (first line), DFS_1 (second line), and DFS_2 (third_line) runs at the horizontal resolution of T126, T254, and T510.

4.3 모델 실행 시간

마지막으로 Table 4는 실험2 검증사례에 대한 수행 시간 정보를 보여준다. 15일 간의 모델 적분 과정을 위해 소모된 실제 모델 실행 시간이 총 몇 초 인지를 사용한 역학코어의 종류 별로, 그리고 사용한 수평 해상도 별로 각각 ‘Total’ 항목에 기록하였다. 해당 실험들은 모두 32개의 코어가 병렬로 연결(노드 수: 2, 노드 당 CPU 수: 16)되어 있는 클러스터 컴퓨팅 환경에서 진행되었다. DFS의 총 모델 실행 시간 바로 옆의 괄호 안에는 동일 해상도의 SPH 실행 시간과 비교하여 그 비율이 몇 퍼센트인지 기록되어 있다. 총 모델 실행 시간에 있어서는 DFS_1과 DFS_2는 SPH에 비하여 명확히 빨랐으며, Cheong (2006)과 Park et al. (2013) 결과와 같이 고해상도로 갈수록 시간이 줄어드는 비율이 커져갔다. 이는 앞서 1장에서 언급한 것처럼 남북방향 변환에서 사용되는 계산량이 격자점의 개수 N에 의존하기 때문이다. 한편, DFS_1과 DFS_2의 비교에서는 DFS_2가 DFS_1에 비해 평균적으로 약 10%의 시간을 더 소모함을 확인할 수 있었다. 이는 앞서 2장에서 언급한 것처럼 DFS_2는 삼중 대각 행렬의 역행렬 연산을 추가로 더 필요로 하기 때문에 예상이 된 결과이다. 물론, 추가적으로 소모되는 시간의 양이 SPH가 DFS_1에 비해 더 소모되는 시간의 양 보다는 훨씬 더 적기 때문에 DFS_2도 SPH에 비해 여전히 계산량 감소의 우수성을 가지고 있다고 할 수 있다.

Table 4.

Comparison of wall-clock time (in second) required for 15-day forecast in the SPH, DFS_1, and DFS2 runs at the resolutions of T126, T254, and T510. The parenthesis represents the relative ratio to the SPH’s wall-clock time.

‘Total’ 외의 항목들은 역학코어 내부의 서브루틴(subroutine) 각각이 총 실행 시간 중 차지한 부분을 표시한 것이다. ‘spectral transform’ 항목은 역학코어 내부에서 매 타임스텝 마다 해당 분광법의 변환과 역변환을 진행하면서 지배 방정식을 푸는데 걸리는 시간에 해당한다. 한편, ‘semi-implicit’과 ‘diffusion’은 모델 안정화를 위해 필요한 과정이다. ‘semi-implicit’은 연직으로 매우 빠르게 전파되는 중력파와 관련된 불안정 요소를 제어하는 준암시적 시간 적분 기법과 관계된 과정을 포함하며, ‘diffusion’은 높은 파수 잡음의 누적으로 인한 불안정 요소를 막기 위하여 추가되는 수평 확산 과정을 포함한다. 그 외에 시간 필터링(Time filter) (Robert, 1966; Asselin, 1972) 과정이나 결과 파일을 작성하는 과정 등의 요소들은 따로 명시하지 않고 ‘Total’ 항목 속에만 포함시켰다. 모델 안정화를 위해 필요한 이 두 항목은 라플라시안(Laplacian) 연산을 필요로 하는데, SPH의 경우 기저함수 그 자체가 라플라시안 연산자의 고유벡터(eigenvector)로 작용하면서 관련 계산들이 매우 간단해지지만, DFS는 그렇지 않다. 따라서 실제 실행 시간에서도 ‘spectral transform’ 항목은 SPH가 DFS들에 비해 확연히 느리지만 나머지 두 항목은 SPH가 오히려 더 빠르다. DFS에 적용되는 ‘semi-implicit’에 대한 정보는 Cheong (2006)에서, ‘diffusion’에 대한 정보는 C2000, Cheong et al. (2002), Park et al. (2013) 등에서 찾을 수 있다. 한편, DFS_1과 DFS_2의 비교에서는 세 항목 모두에서 DFS_2가 DFS_1 보다 평균적으로 시간이 더 걸린 것을 확인할 수 있다. 다만, 대부분의 시간차이는 ‘spectral transform’ 항목에서 나왔으며, 나머지 두 항목에서 난 차이는 그것에 비하면 작다고 할 수 있는 수준이었다.

5. 요약 및 제언

본 연구에서는 DFS 역학코어의 정확도와 계산 효율성을 평가하기 위해 Cheong (2000a, b)의 DFS 방법(DFS_1)과 Yoshimura (2022)의 DFS 방법(DFS_2)을 모두 고려하였고 SPH 역학코어를 규준실험으로 설정하였다. DFS_1과 달리 DFS_2의 경우, 동서방향 파수 m에 따라 sinθ sinnθ와 같이 일반적이지 않은 푸리에 변환이 필요할 경우 이를 원함수와 급수 간의 차이의 제곱 적분 값을 최소화하는 방식을 사용한다. 이러한 방식을 사용하면 극 지점(θ = 0, π) 근처에서 원함수가 0이 아닌 값이 생기는 경우 그것을 표현하는 급수에서 발생하는 깁스 현상을 줄일 수 있다(Fig. 2). Yoshimura (2022)는 단일 층 구면 위에서 천수방정식을 해결하는 데에 이 방안을 사용하였으나, 본 연구에서는 3차원 대기 역학코어로 확장하여 그 성능을 검증했다. 성능 검증에 사용된 실험은 역학코어 상호 비교 프로젝트(DCMIP)의 이상화 실험 중 수평방향 이동에 깊이 연관되어 있는 DCMIP 1-1(실험1)과 DCMIP 4-1(실험2)이다.

물질 수송을 모의하는 실험1에서는 SPH와 DFS가 유의미한 성능의 차이를 보이지 않았으며, DFS_1과 DFS_2의 차이 또한 미미한 수준이었다(Table 1). 다만, DFS_1 대비 DFS_2는 깁스 현상을 약하게 완화하면서 비물리적 값인 음수 값(undershoot)이 존재하는 영역도 적절히 감소시켰다(Fig. 3). 한편, 경압 파동의 발달을 모의하는 실험2에서는 SPH와 DFS의 차이가 일부 확인되었다. 두 역학코어 간의 차이는 저해상도 실험에서 두드러졌으며, 수평 해상도가 높아질수록 그 차이가 줄어들었다. 특히 해상도에 대한 민감도는 SPH보다 DFS가 더 큰 것으로 확인하였다. DFS에서 낮은 수평 해상도를 사용하였을 때 극점의 값을 SPH보다 과소모의하는 경향이 확인되었지만, 높은 수평 해상도에서는 차이가 발생하지 않았다(Fig. 6). 한편, 실험2에서도 DFS_1과 DFS_2는 유의미한 차이를 보이지 않았다(Table 3). 다른 방식으로 변환 과정을 계산했음에도 불구하고 두 결과 간의 차이는 해당 변수가 가지는 스케일의 약 10^-3배 이하 수준 이었다.

실험2의 모델 적분 시간을 비교하였을 때, DFS는 SPH에 비해 유의하게 짧았다(Table 4). 또한 DFS_1은 DFS_2보다 계산 비용이 약 10% 적었다(Table 4). DFS_1과 DFS_2의 계산 결과가 유사함을 고려한다며, 계산 시간 증가를 감수하면서 DFS_2를 채택해야 할 필요성은 없는 것으로 보인다. DFS_1과 DFS_2 실험결과가 유사한 이유로는 역학코어의 안정화를 위해 같이 실행되는 두 가지 요소를 고려할 수 있다. 하나는 극지방 필터링으로 남북방향으로 Fig. 2b와 같은 함수가 나오는 것을 미연에 방지한다. 다른 하나는 수평 확산으로 모델 실행 중에 생성되는 높은 파수의 잡음을 제거하여 DFS_1과 DFS_2의 차이를 줄이게 된다.

결론적으로 DFS를 이용한 역학코어는 SPH와 비교하였을 때 비교적 해상도에 더 민감하게 반응하는 특성을 지니지만 SPH에 비해 향상된 계산 효율성을 보인다. 고해상도에서는 DFS와 SPH의 차이가 최소화되며, 계산 시간이 줄어드는 비율 역시 증가한다. 그러므로 고해상도 모델링에서 SPH의 대체재로서 DFS를 사용하는 것을 고려해 볼 수 있다. 한편, DFS_1과 DFS_2의 비교에서는 뚜렷한 모의 성능 차이를 확인할 수 없으며, 계산 효율성을 고려하였을 때 DFS_2보다 DFS_1을 사용하는 것이 더 적합할 것으로 보인다. 하지만 DFS 방법의 발전을 위해서는 DFS_2와 같은 시도를 지속해야 할 것이며, Y2022에 제안된 다른 방안들도 추가로 도입하여 개선효과를 살펴볼 필요가 있다.

본 연구에서는 역학코어 만을 실행하여 성능을 비교하는 실험을 진행하였다. 향후 물리 과정을 모두 포함시킨 모델 실험을 통해 SPH, DFS_1, DFS_2의 특성을 파악할 필요가 있다. 특히, 실제 GCM에 적용하는 것을 염두에 두고 새로운 DFS 방안을 개발을 이어가야 할 것이다.

Acknowledgments

본 논문의 개선을 위해 좋은 의견을 제시해 주신 두 분의 심사위원께 감사를 드립니다. 이 연구는 환경부의 재원을 지원받아 한국환경산업기술원 “신기후 체제 대응 환경기술개발사업”의 연구개발을 통해 창출되었습니다(2022003560004).

References

• Asselin, R., 1972: Frequency filter for time integrations. Mon. Wea. Rev., 100, 487-490. [https://doi.org/10.1175/1520-0493(1972)100<0487:FFFTI>2.3.CO;2]
• Boyd, J. P., 1978: The choice of spectral functions on a sphere for boundary and eigenvalue problems: A comparison of Chebyshev, Fourier and associated Legendre expansions. Mon. Wea. Rev., 106, 1184-1191. [https://doi.org/10.1175/1520-0493(1978)106<1184:TCOSFO>2.0.CO;2]
• Cheong, H.-B., 2000a: Double Fourier series on a sphere: applications to elliptic and vorticity equations. J. Comput. Phys., 157, 327-349. [https://doi.org/10.1006/jcph.1999.6385]
• Cheong, H.-B., 2000b: Application of double Fourier series to shallow water equations on a sphere. J. Comput. Phys., 165, 261-287. [https://doi.org/10.1006/jcph.2000.6615]
• Cheong, H.-B., 2006: A dynamical core with double Fourier series: Comparison with the spherical harmonics method. Mon. Wea. Rev., 134, 1299-1315. [https://doi.org/10.1175/MWR3121.1]
• Cheong, H.-B., I.-H. Kwon, and T.-Y. Goo, 2004: Further study on the high-order double-Fourier-series spectral filtering on a sphere. J. Comput. Phys., 193, 180-197. [https://doi.org/10.1016/j.jcp.2003.07.029]
• Cheong, H.-B., I.-H. Kwon, T.-Y. Goo, and M.-J. Lee, 2002: High-order harmonic spectral filter with the double Fourier series on a sphere. J. Comput. Phys., 177, 313-335. [https://doi.org/10.1006/jcph.2002.6997]
• Choi, S.-J., and S.-Y. Hong, 2016: A global non-hydrostatic dynamical core using the spectral element method on a cubed-sphere grid. Asia-Pac. J. Atmos. Sci., 52, 291-307. [https://doi.org/10.1007/s13143-016-0005-0]
• Cooley, J. W., and J. W. Tukey, 1965: An algorithm for the machine calculation of complex Fourier series. Math. Comput., 19, 297-301. [https://doi.org/10.2307/2003354]
• Diamantakis, M., and F. Váňa, 2022: A fast converging and concise algorithm for computing the departure points in semi-lagrangian weather and climate models. Quart. J. Roy. Meteor. Soc., 148, 670-684. [https://doi.org/10.1002/qj.4224]
• Hong, S.-Y., and Coauthors, 2013: The Global/Regional Integrated Model system (GRIMs). Asia-Pac. J. Atmos. Sci., 49, 219-243. [https://doi.org/10.1007/s13143-013-0023-0]
• Hortal, M., and A. J. Simmons, 1991: Use of Reduced Gaussian Grids in Spectral Models. Mon. Wea. Rev., 119, 1057-1074. [https://doi.org/10.1175/1520-0493(1991)119<1057:UORGGI>2.0.CO;2]
• Hoskins, B. J., and A. J. Simmons, 1975: A multi-layer spectral model and the semi-implicit method. Quart. J. Roy. Meteor. Soc., 101, 637-655. [https://doi.org/10.1002/qj.49710142918]
• Jablonowski, C., and D. L. Williamson, 2006: A baroclinic instability test case for atmospheric model dynamical cores. Quart. J. Roy. Meteor. Soc., 132, 2943-2975. [https://doi.org/10.1256/qj.06.12]
• Kanamitsu, M., 1989: Description of the NMC global data assimilation and forecast system. Wea. Forecasting, 4, 335-342. [https://doi.org/10.1175/1520-0434(1989)004<0335:DOTNGD>2.0.CO;2]
• Koo, M.-S., and S.-Y. Hong, 2013: Double Fourier series dynamical core with hybrid sigma-pressure vertical coordinate. Tellus A, 65, 19851. [https://doi.org/10.3402/tellusa.v65i0.19851]
• Koo, M.-S., and Coauthors, 2023: The Global/Regional Integrated Model System (GRIMs): an Update and Seasonal Evaluation. Asia-Pac. J. Atmos. Sci., 59, 113-132. [https://doi.org/10.1007/s13143-022-00297-y]
• Lauritzen, P. H., and J. Thuburn, 2012: Evaluating advection/transport schemes using interrelated tracers, scatter plots and numerical mixing diagnostics. Quart. J. Roy. Meteor. Soc., 138, 906-918. [https://doi.org/10.1002/qj.986]
• Merilees, P. E., 1973: An alternative scheme for the simulation of a series of spherical harmonics. J. Appl. Meteor. Climatol., 12, 224-227. [https://doi.org/10.1175/1520-0450(1973)012<0224:AASFTS>2.0.CO;2]
• Nair, R. D., and P. H. Lauritzen, 2010: A class of deformational flow test cases for linear transport problems on the sphere. J. Comput. Phys., 229, 8868-8887. [https://doi.org/10.1016/j.jcp.2010.08.014]
• Orszag, S. A., 1974: Fourier series on spheres. Mon. Wea. Rev., 102, 56-75. [https://doi.org/10.1175/1520-0493(1974)102<0056:FSOS>2.0.CO;2]
• Park, H., S.-Y. Hong, H.-B. Cheong, and M.-S. Koo, 2013: A Double Fourier Series (DFS) Dynamical Core in a Global Atmospheric Model with Full Physics. Mon. Wea. Rev., 141, 3052-3061. [https://doi.org/10.1175/MWR-D-12-00270.1]
• Robert, A. J., 1966: The integration of a low order spectral form of the primitive meteorological equations. J. Meteor. Soc. Japan, 44, 237-245. [https://doi.org/10.2151/jmsj1965.44.5_237]
• Sela, J. G., 1980: Spectral modeling at the National Meteorological Center. Mon. Wea. Rev., 108, 1279-1292. [https://doi.org/10.1175/1520-0493(1980)108<1279:SMATNM>2.0.CO;2]
• Suda, R., 2005: Fast spherical harmonic transform routine FLTSS applied to the shallow water test set. Mon. Wea. Rev., 133, 634-648. [https://doi.org/10.1175/MWR-2871.1]
• Tygert, M., 2008: Fast algorithms for spherical harmonic expansions, II. J. Comput. Phys., 227, 4260-4279. [https://doi.org/10.1016/j.jcp.2007.12.019]
• Ullrich, P. A., C. Jablonowski, J. Kent, P. H. Lauritzen, R. D. Nair, and M. A. Taylor, 2012: Dynamical core model intercomparison project (DCMIP) test case document, 83 pp [Available online at http://www-personal.umich.edu/~cjablono/DCMIP-2012_TestCaseDocument_v1.7.pdf, ].
• Ullrich, P. A., and Couthors, 2017: DCMIP2016: a review of non-hydrostatic dynamical core design and intercomparison of participating models. Geosci. Model Dev., 10, 4477-4509. [https://doi.org/10.5194/gmd-10-4477-2017]
• Williamson D. L., 2007: The evolution of dynamical cores for global atmospheric models. J. Meteor. Soc. Japan, 85B, 241-269. [https://doi.org/10.2151/jmsj.85B.241]
• Yee, S. Y. K., 1981: Solution of Poisson’s equation on a sphere by truncated double Fourier series. Mon. Wea. Rev., 109, 501-505. [https://doi.org/10.1175/1520-0493(1981)109<0501:SOPEOA>2.0.CO;2]
• Yoshimura, H., 2022: Improved double Fourier series on a sphere and its application to a semi-implicit semi-Lagrangian shallow water model. Geosci. Model Dev., 15, 2561-2597. [https://doi.org/10.5194/gmd-2021-168]