1. 서 론

우리나라에서 발생하는 대부분의 자연재해는 여름철 집중호우와 태풍 등으로 인한 수재해이다(Kim et al., 2013a). 최근 들어, 집중호우로 인한 범람과 침수 등의 현상이 증가하고 있는 추세이다. 이와 같은 수재해를 예측하고 방지하기 위한 대책 중 하나로써 강우유출모형에 대한 연구가 지속적으로 활발하게 수행되고 있다(Esteves et al., 2000; Tayfur and Kavvas, 1994; Rousseau et al., 2015). 특히 컴퓨터의 연산성능 향상으로 인해 강우와 지형 및 토양특성의 공간분포를 고려하여 강우유출을 계산할 수 있는 물리기반의 분포형 강우유출모형에 대한 활용이 증가하고 있다(Park et al., 2007). 그러나 분포형모형은 집중형모형에 비해 입력자료의 구축에 많은 시간과 노력이 필요하며, 특히 유출 계산 시 컴퓨터의 연산 수행시간이 상대적으로 오래 걸리는 단점이 있다. 이와 같은 이유로 인하여 분포형모형을 대유역에 적용하는 경우 격자의 해상도를 충분히 높이지 않고 사용되어 왔다(Chung et al., 2010).

지속적인 컴퓨터의 연산성능 향상으로 인하여, 전산유체역학 분야에서는 유체흐름을 보다 사실적으로 계산하고 표현하려는 연구가 활발하게 수행되고 있다. 수치해석을 이용하여 이와 같은 목적을 달성하기 위해서는 고해상도로 계산영역을 분할해야 하므로, 연산수행 시간을 단축시키기 위한 연구가 매우 활발하게 진행되어왔다. 기존에는 여러 개의 중앙처리장치(Central Processing Unit, CPU)를 이용하는 MPI 기법이나 다수의 Thread를 이용하는 OpenMP 기법과 같이, CPU 기반의 병렬처리 기법이 활발하게 연구 및 활용되어왔다. 최근에는 그래픽처리장치(Graphic Processing Unit, GPU)를 이용한 병렬처리기법에 관한 연구가 활발히 진행 중이다.

CPU 내부에는 다양한 작업을 수행하기 위해 물리적 구조가 복잡한 산술논리장치(Arithmetic Logic Unit, ALU)가 있다. 반면 GPU 내부에는 비교적 단순한 ALU가 있지만 CPU보다 다수의 ALU로 구성되어 있어서, CPU보다 한 번에 더 많은 연산수행이 가능하다. GPU는 CPU와 하드웨어적 차이로 인해 CPU와는 다른 방식의 프로그래밍 기법이 필요하다. GPU를 이용한 프로그래밍은 NVIDIA사의 CUDA (Compute Unified Device Architecture)와 Khronos Group의 OpenCL (Open Computing Language)이 대표적이다(Vanderbauwhede and Takemi, 2013). 이들 중 GPU 제조업체인 NVIDIA는 컴파일러(compiler) 개발사인 The Portland Group과 공동으로 GPU를 연산에 이용할 수 있는 컴파일러를 개발했으며, 그중 하나가 CUDA fortran 컴파일러이다. 지난 수년에 걸쳐 전산유체역학 분야에서는 GPU를 이용한 병렬연산에 대한 연구가 진행되어왔다. Gomez Gesteira et al. (2012a; 2012b)은 자유표면 유체해석이 가능한 GPU 기반 smoothed particle hydrodynamics 모형인 SPHysics (www.sphysics.org)를 개발하였다. Vanderbauwhede and Takemi (2013)는 OpenCL 기반 Weather Research and Forecasting (WRF) 모형을 개발했으며, GPU를 이용한 모형의 연산 수행시간은 CPU를 이용한 모형의 연산 수행시간에 비해 O (10²)배 단축되었다. Kim et al. (2013b)은 CUDA fortran 기반 WRF 모형을 개발하였으며, GPU를 이용한 연산 수행시간이 CPU를 이용한 연산 수행시간에 비해 약 5배 단축되었다. Chang et al. (2014)은 고차정확도 수치기법기반의 압축성 유동 해석모형 개발에 CUDA를 이용하였다. GPU를 이용한 연산수행시간은 CPU의 12개 코어를 사용한 OpenMP 병렬연산 수행시간보다 최대 약 2배 단축되었다.

본 연구에서는 GPU 가속기법을 이용하여 강우유출을 물리적으로 해석하는데 활발하게 연구되고 있는 운동파모형을 개발하였다. 또한 개발된 GPU 가속 운동파모형의 정확성과 연산수행시간 성능에 대한 검토 결과를 제시하였다.

2. 지배방정식

2.1 운동파모형

지표유출해석을 위한 운동파모형의 연속방정식과 운동방정식은 Eqs. (1) and (2)와 같다(Liu et al., 2004).

여기서, (x,y)는 평면공간좌표의 방향을 나타내고, h는 수심, t는 시간이다. q_x와 q_y는 x와 y축 방향 단위폭 당 유량이다. r은 강우강도, f는 침투율이다. S_0x와 S_0y는 x와 y축 방향 바닥경사이며, S_fx와 S_fy는 x와 y축 방향 마찰경사이다. Eqs. (2a) and (2b)의 마찰경사는 조도계수를 이용하여 Eqs. (3a) and (3b)로 나타낼 수 있다(Govindaraju et al., 1992).

여기서, n은 조도계수, u와 v는 x와 y축 방향 유속이다. Eqs. (2) and (3)을 정리하여 단위폭 당 유량으로 나타내면 Eqs. (4a) and (4b)와 같이 정리할 수 있다.

Eq. (1)을 유한체적법으로 이산화하면 Eq. (5)와 같다(Liu et al., 2004).

여기서, (i,j)는 이산화된 공간격자의 순서를 나타낸다. △t는 시간간격, k는 시간단계, △x와 △y는 x축과 y축 방향 격자크기이다. △t 동안의 x와 y축 방향 유입량은 Eqs. (6) and (7)과 같다.

△t동안의 x축과 y축 방향 유출량은 Eqs. (8) and (9)와 같으며, 수치해석기법에 대한 보다 자세한 내용은 Liu et al. (2004)에 기술되어있다.

2.2 침투 모형

균질한 흙을 통한 연직방향 침투량을 산정하기 위한 다수의 모형들이 제안되었다(Green and Ampt, 1911; Horton, 1939; Mein and Larson, 1973). 실제 침투는 흙의 함수특성, 균질성, 포화도 등에 영향을 받는다. 따라서 본 연구에서는 다양한 침투모형 중에서 물리기반의 침투모형인 Green-Ampt모형을 이용하였다. Green-Ampt모형은 1차원 연직흐름, 침윤전선에서의 일정한 모관흡수력, 그리고 침윤전선 위 부분에서의 함수량 결손과 투수계수가 일정하다는 가정을 하고 있다. Mein and Larson (1973)은 Green-Ampt모형을 일부 수정하여 담수시간에 따른 강우강도에 대한 침투율을 Eq. (10)과 같이 제시하였다.

여기서, 담수시간은 t_p = F_p / r이고, F_p = h_f$\mathit{∆}\theta$/(r/K_s－1), K_s는 투수계수, h_f는 습윤전선의 모관흡인수두이다. Eq. (10)의 t ＞ t_p에 대한 Green-Ampt모형은 초기 포화조건에서 적용되므로 누가침투량은 다음과 같이 유도된다.

여기서, t_s는 t = 0부터 침투율과 강우강도가 같아질 때까지의 시간이고 $\mathit{∆}\theta$는 토양수분부족량이다. Eq. (11)의 비선형항의 해를 구하기 위해 Newton-Raphson기법을 이용하였다.

3. GPU 가속 기법

3.1 GPU의 하드웨어 구조와 특징

GPU는 처리할 수 있는 인스트럭션(instruction)이 단순한 산술연산에 국한되어 있어 CPU에 비해 단순한 제어부와 훨씬 작은 크기의 캐시 메모리를 가지고 있다. 하지만, GPU는 수천 개의 스트림 프로세싱(stream processing)을 통해 동일한 형태의 산술연산을 동시다발적으로 처리하는데 매우 우수한 장점이 있다(Ruetsch and Fatica, 2013). Fig. 1은 CPU와 GPU의 하드웨어 구조를 간략하게 나타내고 있다. 제어장치(Control unit)는 프로세서의 조작을 지시하는 장치로 입출력 장치 간 통신 및 조율을 제어한다. ALU는 산술연산과 논리연산을 계산하는 디지털 회로이다. Fig. 1과 같이 CPU는 소수의 ALU로 구성되어있는 반면, GPU는 다수의 ALU를 가지고 있다. 본 연구에 사용된 장비를 기준으로 CPU (Intel Xeon E5-2620 v2 @ 2.10GHz)에는 ALU를 포함하고 있는 코어가 6개이며, GPU (NVIDIA GeForce GTX TITAN Z)에는 5760개의 코어가 있다.

Fig. 1.

Hardware architecture of CPU and GPU (Kim, 2019)

Fig. 2는 CUDA의 데이터 처리 흐름을 나타낸다. CPU와 GPU는 자신이 제어할 수 있는 독립적인 메모리를 가지고 있으며, 상대방의 메모리는 제어할 수 없다. 따라서 먼저 CPU 메모리에서 GPU 메모리로 데이터를 전송하고, CPU가 GPU에 데이터 처리를 명령한다. 데이터 처리 명령을 받은 GPU는 각 코어에서 데이터를 병렬로 처리하고, 결과를 GPU 메모리를 거쳐 CPU 메모리로 전송시킨다. GPU를 사용한 데이터 처리 흐름은 CPU와 GPU의 하드웨어 구조의 차이로 인해 CPU만을 이용한 데이터 처리 흐름과 비교 시 CPU와 GPU간의 데이터 전송이라는 추가적인 수행 과정이 발생하게 된다. 따라서 이와 같은 CPU와 GPU간의 데이터 전송으로 인한 추가 소요시간이 발생한다.

Fig. 2.

CUDA data processing flow (Kim, 2019)

3.2 GPU 가속 프로그래밍

Fig. 3은 CPU 기반의 순차연산기법과 GPU 기반의 병렬연산기법을 활용하여 2개 배열의 합을 구하는 코드의 한 가지 예이다(Kim, 2019). 본 연구에서 GPU 가속모형의 개발을 위해 사용한 CUDA는 계산영역내의 격자 개수에 대한 그리드(Grid)와 블록(Block)을 생성한다. 그리드는 블록의 모임이고, 블록은 쓰레드(thread)의 모임이며, 쓰레드는 GPU 내에서 작동되는 하나의 코어를 의미한다. CUDA의 병렬연산은 코어 하나에 한 개의 명령어가 할당되고, 한 블록 내의 쓰레드들이 동시에 병렬로 연산되는 방식이다. GPU 프로그램에서는 Fig. 3(b)와 같이 i와 j를 사용한 반복문으로 표현하지 않으며, 필요시 반복문을 표현 할 수도 있다. Fig. 3과 같이 쓰레드의 배열은 blockDim 변수에 저장되고, threadidx 변수는 index가 주어진다. 또한 블록의 배열은 gridDim 변수에 저장되고, blockidx 변수는 index를 주어진다(NVIDIA, 2011).

Fig. 3.

Example of GPU programming

연산 수행시간을 단축시키기 위해서 GPU 가속 수치모형 구현 시 CPU와 GPU간 대용량 데이터 전송이 빈번하게 발생하지 않도록 프로그래밍을 해야 한다. 따라서 대용량 연산시에는 CPU와 GPU간의 자료전송 빈도를 줄이고 GPU에서 연속적으로 연산이 수행되도록 프로그래밍을 해야 한다. Figs. 4 and 5는 본 연구에서 사용된 GPU 가속 운동파모형의 프로그램체계를 간략히 나타낸 것으로, 초기조건의 입력과 매우 단순한 연산을 제외한 주요부분이 GPU내에서 연산이 연속적으로 수행됨을 나타낸다. 예를 들어, 해당 코드의 일부를 보여주는 Fig. 5의 3행은 시간(t)을 GPU로 전송하는 과정을 나타내는데, 이는 하나의 변수로써 처리시간이 매우 짧아 수치모의의 전체시간에 영향을 거의 주지 않는다. Fig. 5의 4~8행은 2장에 제시되어 있는 운동모형의 이산화식과 침투모형을 풀기 위한 부분으로 대용량의 데이터와 연산과정이 필요하여 GPU 내에서 연속적으로 수행되도록 하였다.

Fig. 4.

Structure of gpu-accelerated kinematic wave model. Solid line: data transfer. Dashed line: flow sequence

Fig. 5.

GPU-accelerated code for kinematic wave model

4. 수치모의결과

4.1 GPU 가속모형 모의결과 검토

4.1.1 Delestre et al. (2008)의 수리실험

Delestre et al. (2008)의 수리실험 결과와 GPU 가속 운동파모형을 이용한 수치모의 결과를 비교하여 제시된 모형의 정확성을 검토하였다. Fig. 6은 수리실험의 측면도로써, 바닥의 길이는 L_x= 4 m이고, 바닥경사는 S_0x= 0.0496이다. 강우강도는 하류단으로부터 0.05 m까지를 제외한 전 영역에 모의 시작 후 5초부터 125초까지 r = 50.76 mm / hr의 강우조건을 부여하였다. 수치모의에 사용된 공간격자의 크기는 $\mathit{∆}x=\mathit{∆}y=$ 0.01 m, 시간간격은 $\mathit{∆}t$ = 0.01 s이고 조도계수는 n = 0.013 s/m^1/3을 이용하였다.

Fig. 6.

Schematic diagram of rainfall runoff experiment (Delestre et al., 2008)

Fig. 7은 수로 하류단에서의 유출수문곡선으로써, GPU 가속 운동파모형의 계산결과는 Delestre et al. (2008)의 실험결과와 잘 일치하고 있다. 또한 GPU 가속 운동파모형과 동일한 수치기법을 이용하여 개발된 CPU 기반 운동파모형의 수치모의 결과와도 매우 유사하다.

Fig. 7.

Measured and computed runoff hydrographs at outlet

4.1.2 Fernández Pato et al. (2016)의 수치모의

침투가 발생하는 일차원 유역에서의 검증을 수행하기 위하여, Fernández Pato et al. (2016)이 제시한 천수방정식(Shallow water equation)모형의 수치모의 결과와 본 연구에서 제시한 GPU 가속 운동파모형의 수치모의 결과를 비교하였다. 유역의 길이는 2000 m이고 바닥경사는 전구간에 걸쳐 S₀ = 0.005이다. 유역의 전구간에 걸쳐 r = 450 m/hr의 강우가 125분 동안 지속되며, 침투관련 매개변수는 K_s = 3.272×10^-5m/s, h_f = 0.0495 m, $\mathit{∆}\theta$ = 0.38이고, 조도계수는 n = 0.03 s/m^1/3이다. 격자크기는 $\mathit{∆}x=\mathit{∆}y=$1 m, 시간간격은 $\mathit{∆}t$ = 0.1 s를 이용하였다. Fig. 8은 유역출구에서의 유출수문곡선으로 GPU 가속 운동파모형의 모의결과가 CPU 기반 운동파모형의 모의결과 및 Fernández Pato et al. (2016)이 제시한 천수방정식을 이용한 수치모의 결과와 잘 일치하고 있음을 보여준다.

Fig. 8.

Computed results of kinematic wave model and shallow water equation model

4.1.3 V-형 유역에 대한 수치모의

2차원 유역에서의 검증을 위하여, 본 연구에서 개발된 GPU 가속 운동파모형을 V-형태의 지형 (Fig. 9)에 적용하고 Di Giammarco et al. (1996)이 제시한 동역학파모형의 수치모의 결과와 비교하였다. Fig. 10은 V-형 유역으로, 중앙수로는 20 m × 1000 m의 평면이고 좌우측 유역은 800 m × 1000 m인 평면으로 이루어졌다. 바닥경사는 중앙수로와 수직방향으로는 S_0x = 0.05이고, 수로의 종방향으로는 S_0y = 0.02이다. 강우강도는 전 영역에서 r = 10.8 mm/hr이다. 수치모의 시 격자크기는 $\mathit{∆}x=\mathit{∆}y=$10 m, 시간간격은 $\mathit{∆}t$ = 1.0 s이다. 조도계수는 중앙수로에서 n = 0.15 s/m^1/3이고, 좌우 유역에서 n = 0.015 s/m^1/3이다.

Fig. 9.

Geometry of V-shaped catchment

Fig. 10.

Computed results of GPU-accelerated kinematic wave model. Vector: velocity. Colormap: discharge

Fig. 10은 시간에 따른 GPU 가속 운동파모형 수치모의 결과로 유량과 유속을 나타낸다. Fig. 11은 유역 하류단에서의 유출수문곡선으로 GPU 가속 운동파모형, CPU 기반 운동파모형 및 동역학파모형(Di Giammarco et al., 1996)의 수치모의 결과이다. Fig. 11과 Figs. 10(a)∼10(c)와 같이 강우가 시작된 후 t = 20 min∼t = 50 min 동안 중앙 수로로 물이 모이고 하류단으로 이동하여, 하류단에서의 유량이 증가하고 있음을 나타내고 있다. 또한 t = 50 min∼t = 90 min에는 수로의 하류부에서 유량이 일정하게 유지되고 있으며, 이는 Fig. 10(d)가 유량이 일정하게 유지되고 있을 때의 유량과 유속을 나타낸다. Fig. 11에서 관찰 할 수 있듯이 GPU 가속 운동파모형의 계산결과는 동역학파모형의 결과 및 CPU 기반 운동파모형의 결과와 잘 일치하고 있다.

Fig. 11.

Computed hydrographs of kinematic wave and dynamic wave models

4.2 GPU 가속 운동파모형의 연산속도

GPU 가속 운동파모형의 연산 수행시간에 대한 성능을 조사하기 위해 GPU 가속 운동파모형의 전체 프로그램 수행시간과 CPU 기반 운동파모형의 전체 프로그램 수행시간을 측정하여 비교하였다. GPU 가속 운동파모형의 연산 수행시간에 대한 성능은 속도향상의 척도인 Speedup을 이용하여 계산하였다.

여기서, GPU time은 GPU 가속 수치모형의 전체 프로그램 수행시간이며, GPU time은 CPU 기반 수치모형의 전체 프로그램 수행시간이다. GPU 가속 운동파모형의 수치모의 시 동일한 예제에서 격자 개수를 다양하게 구성하여 격자 개수와 Speedup의 관계를 확인하였다. Fig. 12는 Fig. 9의 V-shaped 유역에 대한 GPU 가속 운동파모형의 Speedup 결과로써, 격자 개수를 1.6×10⁴부터 6.5×10⁶까지 증가시킬 때, Speedup도 증가하고 있음을 확인하였다. GPU 가속 운동파모형은 CPU 기반 운동파모형보다 최대 450배 정도 Speedup 되었으며, 이는 대규모의 연산 수행에 있어 GPU를 이용하는 것이 1개의 CPU만을 이용하는 것보다 연산 수행시간을 단축시킬 수 있음을 나타낸다. 본 연구에서는 운동파모형 수행시간에 대한 Speedup 결과 중 하나의 결과만을 제시하였지만, 다른 유역을 대상으로 한 모의시에도 대략적으로 비슷한 Speedup 결과가 도출되었다.

Fig. 12.

Speedup of GPU-accelerated kinematic wave model

5. 결 론

본 연구에서는 강우유출모의를 위한 운동파모형의 연산수행시간을 단축시키기 위해 CUDA fortran 컴파일러를 사용한 GPU 가속 운동파모형을 개발하고 이상적인 유역을 대상으로 수치모의 한 유출수문곡선의 정확성을 검증하며 연산 수행시간에 대한 성능을 조사하였다. 본 연구에서 제시한 GPU 가속 운동파모형을 이용한 수치모의 결과는 선행된 실험 및 수치모의 연구 결과와 잘 일치하고 있다. GPU 가속 운동파모형은 격자의 개수가 증가 할수록 Speedup이 최대 약 450배까지 증가하였다. 그러므로 대용량 연산수행에 있어 GPU 가속 운동파모형은 기존 CPU 기반 운동파모형에 비해 연산 수행시간에 대한 성능이 높으며, 지속적인 개발을 수행할 경우에는 실제 대규모 유역에 적용 가능할 것으로 판단된다.

본 연구에서 제시된 결과는 GPU 가속 운동파모형의 정확성과 연산 수행시간에 대한 성능을 검토하기 위한 초기단계의 성과이다. 따라서 GPU 가속 운동파모형의 정확성과 적용성에 대한 평가를 위해서는 실제유역을 포함한 보다 다양한 지형에 대한 적용결과의 평가가 필요하다. 또한 TVD 기법과 같은 보다 수치적으로 안정적인 기법을 이용하여 GPU 가속 운동파모형을 개발을 시도한다면, 보다 안정적이고 적용성이 높은 분포형 강우유출모형의 개발이 가능할 것으로 판단된다.