1편에서는 지구와 하부 전리층 사이의 공간을 하나의 구형 도파관으로 보고 왜 지구를 통째로 계산 영역으로 삼는지를 이야기했습니다. 이번 편에서는 그 연속적인 공간을 컴퓨터가 다룰 수 있는 유한한 셀과 모서리로 바꿉니다. 둥근 지구를 극점의 특이점 없이 표현하면서, 모서리의 방향과 셀의 크기를 계산에 어떻게 반영할까?
구 위에 점을 촘촘히 뿌리는 것만으로는 부족합니다. 유한차분 시간영역(FDTD) 알고리즘이 돌아가려면 격자에 방향이 있어야 합니다. 모서리는 어느 쪽으로 향하는지, 각 면의 둘레는 어느 방향으로 도는지, 이웃한 두 셀이 공유하는 경계에서 부호가 서로 맞는지가 정해져 있어야 합니다. 여기에 모서리 길이와 셀 면적을 반영하면 단순한 값의 차이와 둘레를 따라 더한 값으로 실제 공간 변화율을 계산할 수 있습니다.
하나의 정이십면체(icosahedron)에서 출발해 구조를 쌓아 올리면 결과물은 Yee, 1966의 고전적인 격자 방식을 구면으로 옮겨 온 것과 비슷합니다. 다만 지표면을 따라 전자기장의 차이와 순환을 계산할 때는 삼각형 격자와 그 쌍대인 오각형–육각형 격자를 함께 사용합니다(Simpson et al., 2006). 이 편에서는 공간의 이산화에 집중하고 만들어진 연산자를 맥스웰 방정식의 시간 갱신에 넣는 과정은 3편으로 넘깁니다.
정이십면체에서 구면 격자까지
정이십면체는 꼭짓점 12개, 모서리 30개, 삼각형 면 20개로 시작합니다. 한 번 세분화할 때마다 각 삼각형을 네 개로 쪼갭니다. 세 모서리의 중점을 새로 찍고 그 점들을 다시 구면 위로 끌어올린 뒤 이어 주면 됩니다. 이 과정을 반복할수록 점점 더 촘촘한 삼각형 그물이 지구를 덮습니다(Randall et al., 2002).
세분화 레벨 $L$에서 꼭짓점 수는
\[N_v = 10\cdot4^L + 2\]로 늘어납니다. 한 단계 올라갈 때마다 셀이 네 배씩 촘촘해지는 셈입니다.
그림의 4번은 정이십면체 전체에 세분화와 구면 투영을 한 차례 적용해 1번보다 레벨이 하나 높아진 삼각형 격자를 보여줍니다.
이 삼각형들이 주 격자(primal mesh)입니다. 여기서 짝이 되는 쌍대 격자(dual mesh)를 하나 더 만듭니다. 이웃한 삼각형들의 중심을 이어 주면 원래의 꼭짓점 하나를 정확히 하나의 셀이 둘러싸는 새 격자가 생깁니다. 정이십면체에서 물려받은 꼭짓점 12개는 이웃이 다섯이라 오각형이 되고 나머지 꼭짓점은 이웃이 여섯이라 육각형이 됩니다. 이는 오일러의 다면체 공식에서 따라오는 결과이며 이러한 구조를 골드버그 다면체라고 합니다.
이렇게 만든 구면 격자는 5개의 동일한 패널로 나눌 수 있고 각 패널을 논리적인 직교 배열에 담으면 컴퓨터가 규칙적인 인덱스로 훑을 수 있습니다. 위도-경도 격자와 달리 극점에 특이점이 없으니 이 배열 표현은 지구 전역에서 고르게 동작합니다. 아래는 Simpson & Taflove, 2004가 원 논문의 Fig. 2에서 제시한 패널 전개와 논리적 직교 배열입니다. 실제로 이 배열을 어떻게 다루는지는 4편에서 자세히 이야기합니다.
다만 이 그림엔 북극과 남극을 포함해 정이십면체의 꼭짓점에서 생기는 오각형 셀 12개를 육각형으로 그려 넣은 오류가 존재합니다. 구현할 때는 이 12개 셀만 이웃이 다섯인 오각형으로 따로 처리해야 합니다.
해상도가 결정하는 수치적 관측 가능성
격자를 만드는 법을 정했다고 해서 어떤 크기의 구조든 계산할 수 있는 것은 아닙니다. 정이십면체를 한 단계 세분화할 때마다 삼각형 면은 네 조각으로 나뉘고 세분화 레벨 $L$에서 표면 셀 수는
\[N_v=10\cdot4^L+2\]로 커집니다. 대략적인 규모를 표로 보면 이렇습니다.
| 레벨 | 표면 셀 수 | 대략적인 중심 간격 |
|---|---|---|
| 1 | 42 | 3,765 km |
| 2 | 162 | 1,910 km |
| 3 | 642 | 962 km |
| 6 | 40,962 | 약 120 km |
| 7 | 163,842 | 약 60 km |
| 8 | 655,362 | 약 30 km |
여기서 파장과 관측하려는 구조의 크기를 구분해야 합니다. 20 Hz 파장은 약 15,000 km지만 셀 간격이 120 km인 격자에서 80 km짜리 지질 구조는 여전히 보이지 않습니다. 물질 모델에 상세한 구조를 넣어도 격자가 그보다 성기면 수치 문제 안에서는 구조가 존재하지 않는 것과 다르지 않습니다. 그래서 명령줄 도구는 요청한 이상 구조가 격자 간격보다 작을 때 경고합니다.
낮은 기본 해상도는 설정 점검과 개발에 알맞을 뿐, 작은 지질 목표물을 해석하기 위한 설정은 아닙니다. 얕은 구조를 보려면 수평 격자를 충분히 세분화하는 동시에 근지표의 방사 간격도 1.25 km처럼 촘촘하게 나누거나 노드 위치를 직접 지정해야 합니다. 수평 해상도와 방사 해상도는 서로 다른 방향의 구조를 제한하며 둘 중 어느 하나만 세밀하다고 해결되지는 않습니다.
$\nabla \times$ 연산자의 계산
Taflove & Hagness, 2005, Chapter 3은 Yee 알고리즘을 맥스웰 방정식의 적분형으로 해석합니다. 패러데이 법칙은 한 면을 통과하는 자기선속의 변화가 그 면의 둘레를 따라 측정한 전기장의 순환으로 정해진다고 말합니다. 앙페르 법칙은 반대로 전기선속의 변화가 둘레를 도는 자기장의 순환과 그 면을 통과하는 전류로 정해진다고 말합니다. 따라서 FDTD의 한 번의 갱신은 셀의 각 변에 놓인 접선 성분에 변의 길이를 곱해 부호를 맞춰 더하고 그 합을 셀의 면적으로 나누는 계산으로 볼 수 있습니다.
직교 격자에서는 모든 셀이 같은 직사각형이라 이 기하 정보가 $\Delta x$, $\Delta y$, $\Delta z$ 안에 숨어 있습니다. 측지 격자에서는 셀마다 변의 길이와 면적이 조금씩 다르므로 이를 직접 들고 있어야 합니다. Simpson et al., 2006에서는 이 원리를 구면에 적용해 육각형·오각형 TM 셀의 중심에 $E_r$를 두고 셀 둘레의 모서리에 $H_t$를 놓았습니다. $E_r$를 갱신할 때는 셀을 반시계 방향으로 돌며 각 $H_t$에 해당 모서리 길이를 곱해 더한 뒤 셀 면적으로 나눕니다. 오각형도 항이 하나 적을 뿐 계산법은 같습니다. 짝이 되는 삼각형 TE 셀에서는 세 변의 $E_t$ 순환을 삼각형 면적으로 나눠 $H_r$를 갱신합니다.
여기서 모서리의 방향은 단순한 구현상의 표식이 아닙니다. 셀 둘레를 도는 방향과 모서리 방향이 같으면 더하고 반대면 빼야, 이웃한 두 셀이 공유하는 모서리에서 같은 장을 서로 반대 부호로 사용합니다. 격자에는 그래서 각 셀을 이루는 모서리의 순서와 부호, 실제 모서리 길이, 셀 면적이 함께 필요합니다. 전자는 누구의 값을 어떤 부호로 가져올지를 정하고 후자는 그 합을 실제 공간의 순환과 선속 밀도로 바꿉니다.
같은 구면 격자를 반경만 달리해 여러 층으로 쌓을 때는 단위 구에서 구한 중심각 $\theta$와 입체각 $\Omega$를 재사용할 수 있습니다. 반경 $r$인 층의 접선 길이는 $\ell=r\theta$, 셀 면적은 $A=r^2\Omega$가 됩니다. 다만 층 사이를 잇는 방사 방향 항은 이 관계로 대신할 수 없고 서로 엇갈려 놓인 TM-r 평면과 TE-r 평면 사이의 실제 반경 간격을 사용합니다.
격자를 만든 뒤에는 주 격자의 삼각형 면적을 모두 더한 값과 쌍대 격자의 오각형·육각형 면적을 모두 더한 값이 각각 단위 구의 표면적인 $4\pi$가 되는지 확인합니다. 어느 한쪽이라도 맞지 않으면 셀의 누락, 중복 또는 경계 순서 오류가 있다는 뜻입니다.
이제 구면은 방향 있는 주·쌍대 격자와 길이·면적으로 표현됩니다. 연속 공간의 $\nabla\times$는 접속 표에 따라 모서리 값을 모으고 부호를 맞춰 더하는 유한한 연산이 됐습니다. 아직 파동을 전파한 것은 아닙니다. 다음 3편에서는 이 공간 연산자 위에 $E_r$, $E_t$, $H_r$, $H_t$를 엇갈려 놓고 네 번의 시간 갱신으로 맥스웰 방정식을 전진시킵니다.