해양모델링 2장5~18 2012.7.27 12:26 AM 페이지5 02 모델의 시작 요약 이 장에서는 감쇠 문제를 이용하여 여러분을 수치 모델링 세계로 인도한다. 유한 차분법 의 양해법과 음해법 그리고 일관성, 정확도, 안정도, 효율성 등을 설명한다. 첫 번째 수치 모델의 작성과 결과를 그림으로 보기 위해 FORTRAN 프로그램과 SciLab 스크립트가 사용된다. 2.1 감쇠 문제 2.1.1 문제 소위 말하는 감쇠 문제(decay problem)를 예로 들어 수치 모델링을 소개해보겠다. 이 문제는 수학적으로 다음과 같이 쓸 수 있다. (2.1) 여기서 C는 물질의 농도, t는 시간, k(그리스 문자 카파 )는 양의 상수이다. 기호 d는 어떤 변수의 다른 매개변수(위의 방정식에서 예를 든다면 시간)에 대한 변화를 뜻한다.
해양모델링 2장5~18 2012.7.27 12:26 AM 페이지6 6 오픈소스 소프트웨어를 이용한 해양 모델링 2.1.2 물리적 해석 식 (2.1)의 좌변은 어떤 물질의 단위 시간당 변화율을 나타내며, 우변은 그 양을 나타낸 다. k 5 0이면 C는 처음 값 그대로 농도의 변화는 없다. k fi 0이면 우변이 음수가 되 며 농도는 항상 양수이다. 따라서 C는 시간에 따라 점차적으로 자신의 농도에 대해 일 정한 비율로 감소한다. 수학적으로 식 (2.1)은 1차 상미분 방정식(first-order ordinary differential equation) 이다. 이의 완전한 해를 구하려면 C의 초기값을 알아야 하기 때문에 초기값 문제(initialvalue problem)라고도 한다. 2.1.3 예제 할머니께서는 창고에 10리터짜리 술통에 포도주를 가득 채워두시곤 했는데 나는 그 포 도주의 유혹을 이기지 못해 매일 밤 몰래 창고로 가 술통에서 포도주를 1리터씩 퍼내 고는 할머니께서 눈치 채지 못하시게 빈자리를 물로 채워두었다. 그렇지만 내가 발각 되지 않았다고 해서 할머니께서 포도주의 맛이나 색깔이 변한 사실을 알지 못하셨을 까? 어떻든 여러분은 나의 옛날 술버릇 때문에 포도주의 시간에 따른 농도 변화를 계산 해낼 수 있을 것이다. 농도의 변화를 계산해보고자 할 때는 매일 포도주가 10%씩 없어 진다는 기본적인 사실만 알면 된다. 다음 표는 이 계산의 결과를 나타낸 것이다. 날짜 포도주 양(리터) 포도주 농도(%) 0 10 100 1 9.0 90.0 2 8.1 81.0 3 7.29 72.9 4 6.56 65.6 5 5.9 59.0 6 5.31 53.1 7 4.78 47.8 8 4.3 43.0 9 3.87 38.7 10 3.49 34.9 11 3.14 31.4 날짜 포도주 양(리터) 포도주 농도(%) 12 2.82 28.2 13 2.54 25.4 14 2.29 22.9 15 2.06 20.6 16 1.85 18.5 17 1.67 16.7 18 1.5 15.0 19 1.35 13.5 20 1.22 12.2 21 1.09 10.9 22 0.98 9.8
해양모델링 2장5~18 2012.7.27 12:26 AM 페이지7 제 2 장 모델의 시작 7 포도주 양(리터) 시간(일) 그림 2.1 포도주 농도 변화 그림 2.1은 이 첫 번째 수치 예보 모델의 결과를 단지 펜과 종이만을 이용하여 그림 으로 나타낸 것이다. 예상한 대로 날이 지나면서 포도주 양은 줄어들고 있다. 이 감소 경향은 대략 지수함수적이라 할 수 있다. 2.1.4 SciLab으로 간단한 그래프 그리기 이 실습은 그림 2.1을 SciLab으로 그리는 방법을 보여준다. 1단계:적절한 문서 편집기를 열어 2.1.3절에 있는 표의 3줄(앞머리 항목 이름은 제외하 고)을 입력하고, 예를 들면 winethief.txt 라는 파일에 저장한다. 2단계:SciLab을 실행시켜 방금 만든 파일이 있는 폴더로 이동한다. 3단계:다음 명령어를 입력한다(명령행 끝에는 <return>이 있어야 한다). x=read( winethief.txt,21,3); plot(x(:,1),x(:,2)); xtitle(, Time (days), Wine content (liters), ) read 명령어는 자료를 파일에서 읽어들인다. 매개변수 21 은 자료의 길이를 모른 다는 뜻이다. 3 은 자료의 처음 3칸을 읽으라는 뜻이며, 21 은 칸을 지칭하는 것 은 아니다. plot 명령어는 첫 번째 줄의 값에 대한 두 번째 줄의 값을 그림으로 나
해양모델링 2장5~18 2012.7.27 12:26 AM 페이지8 8 오픈소스 소프트웨어를 이용한 해양 모델링 타낸다. 마지막 줄은 이 그림의 축에 이름을 붙이는 일을 한다. ; 은 SciLab 명령 창에서 출력이 나오지 않게 한다. 4단계:그림창에 있는 File 메뉴의 Edit Figure Properties 나 그림창에서 직접 GED 를 클릭하여 선 굵기, 글자체 등을 바꾼다. 여기서 여러 가지 사용 가능한 옵 션을 이용해보는 것도 바람직하다. 5단계:작성된 그림을 적절한 포맷으로 저장한다. PostScript 포맷으로 출력한 다음 ImageMagick을 이용하여 PNG(Portable Network Graphics) 포맷으로 변경하면 편 리할 것이다. 2.2 유한 차분의 첫 단계 2.2.1 유한 시간 간격과 시간 단계 식 (2.1)을 시간 간격(time step)에 대해 다시 쓰면 다음과 같이 쓸 수 있다. (2.2) 여기서 n은 시간 단계(time level)를 의미하는데, 거듭제곱과 혼동하지 않도록 주의하자. n 5 0이면 계산을 시작하기 전의 초기 상태를 뜻하며, n51은 하나의 시간 간격(Dt) 후의 상태를, n 5 2는 2개의 시간 간격 후의 상태를 뜻한다. 2.2.2 양해법 시간 계산 식 (2.2)에서 미지수는 좌변으로, 기지수는 우변으로 정리하면 (2.3) 여기서 C n50 은 초기 농도로, k 및 Dt와 함께 주어져야 한다. 이 반복 계산 방법은 다음 n 1 1 시간 단계에서의 C 값을 계산하기 위해 이미 알고 있는 n 시간 단계에서의 값을 이용하게 된다. 이를 양해법 전방 시간 계산(explicit time-forward iteration)이라 한다.
해양모델링 2장5~18 2012.7.27 12:26 AM 페이지9 제 2 장 모델의 시작 9 2.2.3 양해법의 수치 안정도 조건 식 (2.3)에서 보는 바와 같이 농도는 시간이 지남에 따라 조금씩 감소하는데, 이 감소하 는 정도는 k { Dt로 주어진다. 이 값이 1을 넘으면 계산되는 농도가 음의 값을 갖게 되므 로 이 값은 절대 1을 넘을 수는 없다.k { Dt > 2가 되면 농도는 음으로 증가하기까지 한 다. 따라서 필요한 조건 (2.4) 을 수치 안정도 조건(condition of numerical stability)이라 한다. 식 (2.3)은 식 (2.4)의 조건하에서만 안정적으로 계산된다. 따라서 시간 간격의 최대값은 k의 값에 따라 결정 되어야 한다. 2.2.4 음해법 시간 계산 식 (2.1)은 다음과 같이 다른 방법으로 차분화할 수 있다. (2.5) 여기서 우변의 농도는 n 1 1 시간 단계의 값으로 주어졌다. 이러한 방법은 독자들에게 약간 이상하게 보이겠지만 미지수와 기지수를 나누어 다시 써보면 다음과 같이 된다. (2.6) 양해법에 비해 이 음해법 전방 시간 계산(implicit time-forward iteration)의 장점은 어 떤 Dt의 값에 대해서도 수치적으로 안정하다는 점이다. 우변의 분모 값은 항상 1보다 크므로 농도는 시간에 따라 감소하고 절대 음의 값이 될 수는 없다. 2.2.5 복합 방법 양해법과 음해법을 함께 사용하여 다음과 같은 차분식을 구성할 수 있다. (2.7)
해양모델링 2장5~18 2012.7.27 12:26 AM 페이지10 10 오픈소스 소프트웨어를 이용한 해양 모델링 여기서 a(그리스 문자 알파 )는 가중치로 0과 1 사이의 값을 갖는다. a 5 1이면 완전 한 음해법이며, a 5 0이면 완전한 양해법이 된다. a 5 0.5인 경우가 소위 말하는 반음 해법(semi-implicit scheme)이다. 2.2.6 그 밖의 방법 여기서 언급되지는 않지만 룽게-쿠타 방법(Runge-Kutta scheme) 이나 애덤스-배쉬 포스 방법(Adams-Bashforth scheme) 처럼 현재 시간 단계와 다음 시간 단계 사이의 중간 시간 간격을 고려하는 방법도 있는데, 이는 대체로 정확도나 효율성이 매우 높다. 2.2.7 일관성 조건 초기 농도 C o 에 대한 식 (2.1)의 정확한 이론해는 다음과 같이 주어진다. (2.8) 여기서 exp 는 지수함수를 뜻한다. 수치 시간 간격 값을 아주 작게 할수록 유한 차분식 의 수치해가 이론해에 더 가까워져 가면 이 수치 모델은 일관성(consistency)이 있다고 한다. 즉, 시간 간격 Dt를 작게 할수록 예보된 농도가 참값에 가까워져 간다는 뜻이다. 2.2.8 정확도 조건 미분 방정식을 차분 방정식으로 전환하는 과정에서 잘려나가는 항에 의한 오류(truncation error)가 발생한다. 또한 컴퓨터 계산에서는 제한된 수치 계산밖에 못하는 데서 발생하 는 반올림 오류(round-off error)도 발생한다. 이 두 오류는 당연히 계산 과정 내내 충분 히 작아야 한다. 2.2.9 효율성 조건 큰 프로그램은 계산 과정에서 자료의 출력이나 보관을 위해 막대한 양의 컴퓨터 공간(메 모리)이 필요하며 계산 시간도 매우 길어진다. 따라서 프로그램의 작성은 컴퓨터가 막대 한 자료에 치여 버벅거리지 않고 효율적으로 작동하여 적절한 시간 내에 작업을 완료할 수 있도록 작성해야 한다.
해양모델링 2장5~18 2012.7.27 12:26 AM 페이지11 제 2 장 모델의 시작 11 2.2.10 프로그램 코드의 작동 과정 컴파일러는 FORTRAN 95 프로그램을 위에서부터 아래로 한 줄씩 컴퓨터 언어로 번 역하는데, 매개변수는 반드시 사용되기 전에 값이 주어져야 한다. 매개변수의 특성에는 정수, 실수, 행렬, 문자, 논리변수 등이 있다. 프로그램의 각 줄에는 x 5 b 1 c 와 같이 좌변에는 하나의 미지수만 올 수 있고, 또 b와 c는 사전에 정의되어 있어야 하며, 최종적으로 x가 결정된다. 2.2.11 첫 번째 FORTRAN 프로그램 모델러 사이의 오래된 전통인, 컴퓨터 모니터에 Hellow World 라고 표시하는 첫 번 째 프로그램을 작성해보자. PROGRAM first write(6,*) Hellow World END PROGRAM first FORTRAN 프로그램은 PROGRAM 이름 문장으로 시작하여 END PROGRAM 이름 문장으로 끝난다. 이 프로그램의 이름은 first 이지만, 여러분은 이름을 마음대로 붙일 수 있다. 이 프로그램을 first.f95 라는 이름의 파일로 저장하자. 2.2.12 FORTRAN 프로그램의 컴파일과 실행 과정 먼저 윈도우의 명령창(시작 => 모든 프로그램 => 보조 프로그램 => 명령 프롬프트)을 열 어 FORTRAN 프로그램 파일이 있는 폴더로 이동한다. 1단계 : 다음의 명령으로 컴파일한다. g95 -o first.exe first.f95 여기서 -o 는 실행 파일의 이름을 정해주는 역할을 한다. 2단계 : 컴파일러가 더 이상 오류 메시지를 출력하지 않을 때까지 오류를 수정한다. 3단계 : 컴파일이 성공적으로 끝났다면 새로 생성된 first.exe 파일을 파일창에서 더블 클릭하거나 명령창에서 first<enter>를 치는 것으로 실행시킬 수 있다.
해양모델링 2장5~18 2012.7.27 12:26 AM 페이지12 12 오픈소스 소프트웨어를 이용한 해양 모델링 결과로 명령창에 Hellow World 가 표시될 것이다. 축하! 2.2.13 FORTRAN 요약 모든 상수, 매개변수, 변수, 행렬은 사용되기 전에 값이 주어져야 하며 정수(23, 0, 1, 3 등)는 실수(1.2, 4.2 25.23 등)와 구별되어야 한다. 논리 연산자(참 또는 거짓)와 문자도 정의되어야 한다. 상수 매개변수는 프로그램의 첫머리에서 다음의 예처럼 정의되어야 한다. INTEGER, PARAMETER :: nx = 11! horizontal dimension INTEGER, PARAMETER :: nz = 5! vertical dimension REAL, PARAMETER :: G = 9.81! acceleration due to gravity REAL, PARAMETER :: RHOREF = 1028.0! reference density REAL, PARAMETER :: PI = 3.14159265359! pi 느낌표 이후의 문자는 단순한 주석으로 컴파일 과정에서 무시된다. 주석은 꼭 필요 한 것은 아니지만 프로그램의 구조를 확인하거나 뒤에 이를 되새기는 데 매우 유용하 게 쓰인다. 매개변수는 프로그램 실행 과정에서 그 값을 변경할 수 있으며 다음과 같이 선언한다. REAL :: wspeed! wind speed INTEGER :: k! grid index CHARACTER(3) :: txt 위의 예에서 txt 는 세 글자로 구성된 문자 자료이다. 1차원 또는 2차원 행렬은 다음 과 같이 정의할 수 있다. REAL :: eta(0:nx+1)! sea-level elevation REAL :: w(0:nz+1,0:nx+1)! vertical velocity nx 5 11 그리고 nz 5 5라면 eta 는 11 1 2 5 13, 즉 아직 값이 부여되지 않은 원 소 eta(0), eta(1),..., eta(11), eta(12)를 가지며 w 는 7행 13열의 2차원 행렬로 선언된 다. 선언된 행렬에는 다음과 같이 DO 루프(loop)를 이용하여 값을 부여할 수 있다.
해양모델링 2장5~18 2012.7.27 12:26 AM 페이지13 제 2 장 모델의 시작 13 DO k = 0,nx+1 IF(k > 50) THEN eta(k) = 1.0 ELSE eta(k) = 0.0 END IF END DO 이 DO 루프는 지수 k가 0에서부터 1씩 증가되어 nx+1이 될 때까지 반복 수행된다. 다른 방법으로 다음처럼 할 수도 있다. DO k = nx+1,0,21 IF(k > 50) THEN eta(k) = 1.0 ELSE eta(k) = 0.0 END IF END DO IF 문의 옵션으로는 > (보다 큼), < (보다 작음), == (같음), >= (보다 크거나 같음), <= (보다 작거나 같음), /= (같지 않음) 등이 있다. IF 문이 한 줄로 끝나면 아래의 예 처럼 위의 예에 있는 ELSE 나 END IF 등은 필요 없어진다. DO k = nx+1,0,21 eta(k) = 0.0 IF(k > 50) eta(k) = 1.0 END DO 출력을 위한 파일은 다음과 같이 하여 생성할 수 있다. OPEN(10, file = Ex1a.txt, form = formatted, status = unknown ) 첫 칸의 10 은 WRITE 나 READ 문에 연결되는 장치 번호이며, file 엔트리는
해양모델링 2장5~18 2012.7.27 12:26 AM 페이지14 14 오픈소스 소프트웨어를 이용한 해양 모델링 파일 이름을, form 엔트리는 출력 파일이 ASCII 형식인지 바이너리(binary) 형식인 지를 구분한다. 여기서는 ASCII 형식의 출력을 택했다. status 엔트리에서 unknown 은 만약 그 이름의 파일이 없으면 새로 만들어서 출력하고, 만약 같은 이름 의 파일이 있으면 그 위에 덮어 쓰라는 뜻이다. 다음에 필요한 파일에 덮어 쓰지 않도 록 주의해야 한다. status 엔트리의 기타 옵션으로는 new 와 old 가 있다. 자료의 출력은 다음과 같이 WRITE 문으로 이루어진다. WRITE(10,*) G 여기서 장치 번호(여기서는 10)는 앞서의 OPEN 문에서의 장치 번호와 연결되며, * 는 표준형 출력 양식을 뜻한다. 장치 번호 6은 앞의 첫 FORTRAN 프로그램에서 본 바와 같이 모니터 출력으로 예약되어 있다. 마찬가지로 READ(5,*) 는 자판에서 읽어들임 을 뜻한다. 다음과 같이 한 번에 여러 개를 출력할 수 있다. WRITE(10,*) eta(10), eta(20), eta(30) 이 문장을 계속 실행시키면 3열의 자료가 출력될 것이다. 다음과 같이 하여 행렬의 한 열을 전부 다 출력할 수도 있다. WRITE(10,*) (eta(k), k = 1, nx) 파일이 더 이상 필요 없어지면 다음과 같이 하여 그 파일을 닫아야 한다. CLOSE(10) 2.3 예제 1 : 감쇠 문제 2.3.1 목표 이 예제의 목표는 FORTRAN 프로그램을 이용하여 양해법이나 음해법으로 식 (2.1)에 따른 물질의 붕괴를 예측해보는 것이다.
해양모델링 2장5~18 2012.7.27 12:26 AM 페이지15 제 2 장 모델의 시작 15 2.3.2 작업 개요 물질의 초기 농도를 100%, 붕괴 상수를 초당 0.0001, 즉 k 5 0.0001 s 21 으로 두고 시 작하자. 시간 간격의 길이를 변경해가며 시간 간격의 길이에 따라 해가 얼마나 정확해 지는지를 알아보기로 하며, 양해법과 음해법 모두 검토해보자. 2.3.3 주의사항 문서 편집기로 FORTRAN 프로그램을 작성하여 파일 이름을 Exercise1.f95 로 하여 저장한다. 파일 이름으로 빈 공간이나 특수문자 등은 사용할 수 없다. 파일 이름을 다 르게 하는 것은 괜찮으나 너무 길면 좋지 않고 이 예제와 관련 있는 이름이 좋을 것이 다. 확장자 f95 는 이 파일이 FORTRAN 95의 원시 프로그램임을 나타낸다. 2.3.4 FORTRAN 프로그램 이 예제의 FORTRAN 프로그램인 winethief.f95 는 부록 CD의 Exercise 1 폴더에 있다. 2.3.5 결과 모델 계산 결과는 출력 파일 output1.txt 와 output2.txt 로 나올 것이다. 프로그램에 서 변수 MODE에 따라 양해법과 음해법의 계산이 이루어진다. 컴파일을 두 번 하는 것을 방지하기 위해 mode 값은 READ(5,*) mode 를 통해 자판에서 읽어들인다. 그림 2.2는 Dt 5 3600 s로 두고 식 (2.3)의 양해법과 식 (2.6)의 음해법 해를 나타낸 것이다. 그림에서 보는 바와 같이 양해법의 해는 참값보다 약간 작게 나오고 음해법의 해는 약간 크게 나왔다. 반음해법을 쓰면 더 좋은 결과를 얻을 수 있는데 이는 독자들 의 몫으로 남겨놓겠다. 시간 간격을 3600초로 하면 PC에서 수 초 이내에 계산이 완료된다. 시간 간격을 Dt 5 1 s와 같이 더 짧게 한다면 계산 결과는 더욱 향상될 것이다. 각자 확인해보자. 2.3.6 추가 예제 이 예제를 복합 방법[식 (2.7)]에서 a를 0.25, 0.5, 0.75 등으로 변경해가면서 계산해보자.
해양모델링 2장5~18 2012.7.27 12:26 AM 페이지16 16 오픈소스 소프트웨어를 이용한 해양 모델링 그림 2.2 시간에 따른 농도(%) 변화. 위는 음해법의 결과이고, 아래는 양해법의 결과이며, 가운데 곡선은 식 (2.8)에 따른 이론해이다. 2.4 오류의 검토와 제거 2.4.1 오류 메시지 FORTRAN 프로그램에 오류가 있으면 컴파일러는 하나 또는 여러 개의 오류 메시지를 내놓는데, 이 메시지를 검토하여 오류를 제거할 때 주의해야 할 사항이 있다. 2.4.2 차례차례 오류 제거하기 처음 오류 메시지를 보면 한 번에 하나씩 오류를 수정하는 게 좋다. 여러 오류 메시지 가 처음의 오류에 연관된 것일 경우가 더러 있기 때문이다. 오류를 수정한 후 확인하지 도 않고 또 다른 오류를 찾는 것은 별로 좋은 방법이 아니다. 오류를 하나 고쳤으면 바 로 컴파일해보도록 하자. 때로는 정리되지 못한 결과로 오류가 발생할 수도 있다. 2.4.3 오류 메시지 무시하기 때로는 오류 메시지가 애매하거나 잘못된 내용일 수도 있다. 이런 경우는 오류 메시지 내용보다 이 메시지에 관련된 프로그램의 문장 줄 번호를 더 주의 깊게 볼 필요가 있다.
해양모델링 2장5~18 2012.7.27 12:26 AM 페이지17 제 2 장 모델의 시작 17 2.4.4 흔히 범하는 오류 흔히 오류는 숫자 0 을 문자 o 로 잘못 치거나 그 반대의 경우, 또 숫자 l (일)을 문자 l (엘)과 혼동하는 경우처럼 단순한 타이핑 실수로 일어나기도 한다. 가끔은 END DO 나 END IF 등을 빠뜨려 오류가 발생하기도 한다. 2.4.5 컴파일러를 믿어라 어떤 경우이든 컴파일러를 신뢰하는 것이 좋다. 오류 메시지가 나오면 여하튼 프로그 램에 오류가 있다는 뜻이다. 오류 메시지 내용은 바르지 못할 수도 있지만, 오류 메시 지가 나왔다면 그곳에는 최소한 1개 이상의 오류가 있음이 틀림없다. 고쳐도 고쳐도 끝 없이 오류가 나오면 휴식을 좀 취하거나 산보를 나가거나 아니면 한숨 자고 나서 다시 보는 것이 바람직하다. 휴식은 언제나 중요하다! 2.4.6 경고 메시지 컴파일 명령어에 -Wall (모든 경고 보여주기)을 추가하면 모든 경고 메시지를 볼 수 있 다. 예제 1에 다음과 같이 해볼 수 있다. g95 -Wall -o winethief.exe winethief.f95 경고 메시지는 프로그램 안에 다른 오류가 있는지를 확인하기 위해 살펴볼 필요가 있다.
해양모델링 2장5~18 2012.7.27 12:26 AM 페이지18