1. 서 론

페리다이나믹(peridynamics)은 비국부 연속체 기반 해석 기술로서, 운동방정식의 응력 미분을 대신하여 유한영역 내 절점 간 상호작용력의 적분을 사용한다(Silling, 2000). 따라서 별도의 이론적 혹은 수치적 보완이 없이도 불연속 현상을 동일한 지배방정식으로 해석이 가능하도록 개발되었다. 다양한 연구를 통해 페리다이나믹이 단일 물체에 대한 재료 모델링 및 동적 파괴 해석을 통해 재료 균열 및 손상 등의 불연속 현상을 분석하는데 유용함이 확인되었다(Ha and Bobaru, 2010; Ha and Cho, 2011; Hu et al., 2012; Bobaru et al., 2012; Ha et al., 2015). Bobaru 등(2012)은 페리다이나믹을 통해 다중적층 유리구조물의 손상 패턴 분석을 시도하였다. 적층구조물은 적층판들 사이에 상대적으로 매우 얇은 층간구조가 반복되는 복합구조물로서, 적층구조물 전체를 정밀하게 수치모델링하려면 막대한 계산용량과 해석 시간이 요구된다. Bobaru 등(2012)에서는 목적하는 다중 적층구조물의 층간구조를 모두 생략하고 개별 적층구조가 직접 접촉하도록 간략하게 모델링하여 실험과 다른 동적 파괴 패턴도 함께 관측되었다. 최근에는 적층구조물에 대한 다물체 충돌 파괴 모델링을 위한 연구가 새롭게 수행되었다(Ha and Ahn, 2018; ,Ha, 2019; Ha, 2020). 해당 연구들에서는 가상 절점을 이용한 가상 층간 구조 모델링을 통해 층간 재료를 생략하지 않고 다중 적층구조물을 효율적으로 모델링하였다. 또한 충돌체와 적층구조물 등을 별도의 물체로 모델링하고 비국부 접촉 모델을 통해 충돌 파괴 해석을 수행하였다.

페리다이나믹 해석법은 복잡하고 다양한 동적 파괴 현상을 분석하기에 유용한 절점 기반 비국부 해석법이지만, 유한요소법과 같은 국부 해석법에 비해 많은 계산 시간이 소요된다. 또한 고속 충돌 시 발생하는 고주파수 영역의 응답을 해석하기 위해서는 매우 세밀한 절점 간격과 짧은 시간 간격이 요구되어 계산 시간 부담이 더욱 가중된다. 본 연구에서는 Ha(2019)와 Ha(2020) 등에서 제안된 다중 적층구조물의 고속 충돌 파괴 해석 기술을 수치적으로 구현하고 최적화 및 병렬 처리를 통해 현실적인 모델링 및 분석이 가능한 다물체 페리다이나믹 해석 코드 개발을 시도하였다. 해석 코드는 Fortran 90/95를 이용해서 개발되었으며 병렬 처리 및 실행을 위해 한국과학기술정보연구원(KISTI)의 슈퍼컴퓨팅센터 슈퍼컴 5호기 누리온(Nurion)의 Intel Fortran MPI compiler와 OpenMP를 사용하였다. 다물체 모델링 및 가상 층간 결합 모델링 등이 고려된 다물체 해석 코드를 최적화 하기 위한 핵심 요소들을 분석하고, 모델 의존성 발생 서브루틴 분석 및 프로세스 통신 데이터 분별을 통해 최적성능을 가지는 MPI-OpenMP hybrid 병렬 처리 구조를 적용하였다. 다물체 충돌 파괴 현상을 시뮬레이션하는 테스트 모델을 통해 병렬 처리 코드의 성능을 확인하였다.

2. 다물체 페리다이나믹 해석

본 장에서는 병렬화 프로그래밍을 적용한 다물체 페리다이나믹 해석 코드에 적용된 페리다이나믹 이론 및 가상 층간 결합 모델과 비국부 접촉 모델에 간략히 설명한다. 해석코드는 절점 기반 외연적 시간 적분(explicit time integration)을 통해 구현되었다.

2.1 페리다이나믹 이론(Peridynamics)

시간 $t$에서 물체 $\Omega$에 속하는 절점 $\mathit{x}$의 위치를 $\mathit{y}(\mathit{x},t)$로 정의할 때, 이에 대한 페리다이나믹 운동방정식(감쇠 무시)은 다음과 같이 정의된다(Silling, 2000).

여기서, $\rho$는 재료 밀도이고 $\mathit{b}$는 체적력 밀도를 의미한다. $H_x$는 절점 $\mathit{x}$와 일정 거리($\delta$) 내에서 상호작용하는 $\mathit{x}\text{'}$의 집합으로 식 (2)와 같이 정의되는 비국부(nonlocal) 영역이다.

페리다이나믹 이론은 $H_x$ 내에서 유효한 상호작용력을 적분하여 운동방정식 (1)에서 응력의 미분을 대신한다. 결합 기반 페리다이나믹 모델에서 미소탄성 균질 등방성 재료에 대한 상호작용력은 다음과 같이 정의될 수 있다.

여기서, 미소재료상수는 3차원 모델에서 $c=\frac{6E}{\pi {\delta }^4(1-2\nu )}$로 유도 되며 결합 신장량 $S=\frac{\left|\mathit{y}\mathbf{\text{'}}-\mathit{y}\right|-\left|\mathit{x}\mathbf{\text{'}}-\mathit{x}\right|}{\left|\mathit{x}\mathbf{\text{'}}-\mathit{x}\right|}$이다(Silling and Askari, 2005). PMB(Prototype micro-brittle) 모델(Silling and Askari, 2005)에서 각 결합은 최대 신장량 $s_0=\sqrt{\frac{5G_0}{6E\delta }}$ 이상에서 손상된다 (식 (4)). 여기서 $G_0$는 최대 에너지 해방률(critical energy release rate)이다.

따라서 절점 $\mathit{x}$에서 국부 손상은 다음과 같다.

결합 기반 페리다이나믹 상세 이론은 Silling(2000), 상태 기반 이론은 Silling 등(2007)을 참조할 수 있다.

2.2 가상 층간 결합 모델 및 비국부 접촉 모델

Fig. 1.

Ghost interlayer model(Ha, 2019)

두 물체 A와 B에 대한 충돌/접촉 현상을 해석하기 위해 다음과 같은 비국부 접촉(nonlocal contact) 모델을 도입한다.

여기서 $d_c$는 두 물체의 최소 절점 거리 중 큰 값으로 두 물체의 접촉 거리를 판별하는 기준이 된다. ${\gamma }_c$는 무차원수이며 $C_{AB}$는 두 물체의 미소재료상수의 평균이다. $H\left(a\right)$는 $a>1$에서 1이 고 다른 경우는 0인 계단함수이다. 접촉력의 방향은 두 물체의 사이의 방향 벡터 {\boldsymbol M}_{\boldsymbol A\boldsymbol B}=\frac{{\mathbf y}_{\mathbf B}\boldsymbol-{\mathbf y}_{\mathbf A}}{\boldsymbol\vert{\mathbf y}_{\mathbf B}\boldsymbol-{\mathbf y}_{\mathbf A}\boldsymbol\vert}로 정해진다.

비국부 접촉 모델은 인접한 두 물체가 서로 독립적으로 움직일 때만 적용하며, 층간 구조로 결합된 영역에는 비국부 접촉 모델을 계산하지 않는다. 전산코드 구현 측면에서는 층간 재료와 상호 작용하는 영역에 있는 절점이 해당한다. 이를 제외한 내부 절점들은 재료 손상 및 관통 등으로 인해 상대 물체에 노출이 되는 경우에 비국부 접촉 모델이 적용된다. 다만 이들에 대해서는 접촉 거리를 절점 간격이 아닌 층간 재료 두께로 설정하였다. 가상 층간 결합 모델 및 비국부 접촉 모델에 대한 보다 상세한 내용 및 수치적 검증은 Ha(2020)에 수록되어 있다.

2.3 다물체 페리다이나믹 해석 코드 구성

앞서 소개된 핵심 모듈들의 상호 구성을 다음의 의사코드(pseudocode)로 표현하였다(Algorithm 1). 비국부 접촉 상호작용(nonlocal contact interactions)을 계산하기 위해서는 먼저 변형된 위치 $\mathit{y}$에서 한 모델의 절점의 접촉 영역에 속하는 상대 모델의 이웃(neighbor) 절점들을 구성한 후, 식 (7)의 접촉력을 계산한다. 가상 층간 결합 상호관계(ghost interlayer interactions)는 각 모델 별로 가상 층간 절점들과 상호작용을 계산한다. 이때 가상 층간 절점들은 모든 모델에서 공유하게 된다. 마지막으로 각 모델에서 개별 절점들에 대해 비국부 영역 $H_x$ 내에서 속한 이웃 절점들을 모두 찾고, 이들의 상호작용을 계산한다. 이 때 이웃 절점 리스트는 절점들의 초기 위치 $\mathit{x}$에 근거한다. 일반적으로는 저장 공간 절약을 위해 반복 계산 순환문 내에서 개별 절점 별로 이웃 절점들을 찾아서 계산에 사용한다. 계산 부담을 줄이기 위해 계산 도메인을 적절히 분할해 두고, 개별 절점 인근의 분할 영역들에 대해서만 근접 검색(proximity search)을 수행한다.

Algorithm 1

Pseudocode for the core peridynamic computations

3. MPI-OpenMP 혼합 병렬 처리

본 장에서는 다물체 페리다이나믹 해석 코드 개발을 위한 핵심 최적화 요소에 대해 분석하고 MPI-OpenMP 혼합 병렬 처리에 대해 설명한다. 해석 코드는 Fortran 90/95를 이용해서 개발되었으며, 한국과학기술정보연구원(KISTI)의 슈퍼컴퓨팅센터 Nurion의 Intel Fortran MPI compiler와 OpenMP를 사용하였다. 개발 코드의 성능 평가를 위해 Intel Vtune Amplifier를 사용하였다.

3.1 OpenMP 병렬 처리 및 최적화

다물체 페리다이나믹 해석 코드의 핵심 최적화 요소를 분석하고 OpenMP 병렬 처리를 수행하였다. MPI-OpenMP 혼합 병렬 처리에 대해서는 다음 장에서 다룬다. 절점 계산 순환문에 대해 전반적으로 OpenMP 병렬 처리를 적용하였다. Algorithm 1에서는 5~6번째 줄, 9~11번째 줄, 14~17번째 줄, 21~22번째 줄 등이 병렬 처리되었다(original OpenMP code). OpenMP는 기존의 Serial 코드에 간단한 명령어 추가만으로 병렬화가 가능하다. 병렬화를 적용하고자 하는 구문의 전후에 !$OMP PARALLEL와 !$OMP END PARALLEL을 선언해 주고, 순환문의 시작과 끝에 !$OMP DO와 !$OMP END DO를 선언해준다. PUBLIC와 PRIVATE 변수선언 및 기타 상세한 OpenMP 병렬화에 대한 내용은 본 논문에 수록하지 않았다. Fig. 2는 OpenMP 병렬 처리 코드의 성능을 intel Vtune Amplifier를 통해 분석한 결과이다. 여기에 사용된 테스트 모델 1(Model 1)은 3개의 물체가 총 187,011개의 절점으로 구성된 3차원 다물체 충돌 파괴 해석 모델(Fig. 3a)로서, 각각 88,641개 절점을 가진 두께와 폭이 0.4m이고 길이 4m인 두 보 구조체가 0.04m의 층간 구조(가상 절점 4,221개)로 결합되어 있으며, 양 끝단이 고정되어 있다. 총탄(9,729개 절점)의 재료물성치는 $E=10Pa,$$\rho =4kg/m^3,G_0=1J/m^2$이며 초기속도 $0.1m/s$로 보에 수직 충돌한다. 보와 층간 구조의 영률($E$)은 각각 $1Pa$와 $0.01Pa$이며 재료밀도와 에너지방출률은 $\rho =1kg/m^3,G_0=0.00001J/m^2$로 동일하다. 사용된 재료는 물리적 의미를 가지지는 않으며 수치해석적 검토를 위해 임의로 배정하였다. 고정된 시간 간격 0.001초를 사용하였으며, 기타 상세한 모델 설명은 본 논문의 목적에 무관하므로 생략한다. 모델 1의 다물체 충돌 파괴 해석 결과(4초)를 Fig. 3b에 표시하였다. ‘DAM’은 손상지수이며 식 (5)를 통해서 0(no damage)과 1(fully damaged) 사이의 값으로 표현된다. 8개의 OpenMP 쓰레드(thread)를 사용하여 4,000회의 반복 계산을 수행하는데 약 3,453초의 시간이 소요되었다.

Fig. 2.

Application performance analysis for the original OpenMP code

Fig. 3.

Multi-body peridynamic simulations for Model 1

다물체 동적 해석에서는 모델 수 및 절점 수의 증가로 이웃 절점 검색에서 소요되는 시간 부담이 가중되어 전체 해석의 효율성을 저하시킨다. Algorithm 1의 15번째 줄에 기술된 이웃 절점 검색은 반복 계산 순환문에 독립적인 절점 초기 위치 정보($\mathit{x}$)에만 의존하므로, 반복 계산 순환문 이전에 모든 모델의 모든 절점들에 대해 이웃 절점들을 미리 찾아서 저장하도록 수정한다. 반복 계산 순환문 내에서는 이웃 절점 목록 갱신 없이 저장된 목록을 사용하도록 코드를 최적화할 수 있다. 이를 위해 각 모델 별로 (총 절점 수)×(이웃 절점 수)의 정수 배열(integer array)을 위한 추가 저장 공간이 요구된다. Algorithm 1의 4번째 줄의 이웃 절점 목록 계산도 최적화가 가능하다. 다만 이는 변형된 위치($\mathit{y}$)에 의존하여 매 반복 계산마다 갱신이 되어야 하기 때문에 저장 공간 할당에 비해 성능 향상 정도가 크지 않아서 최적화하지 않는다. 최적화 코드(optimized OpenMP code 1)는 동일한 조건에서 약 1,663초의 계산 시간(48%)이 소요되어 상당한 성능 개선이 이루어졌음에도 불구하고 OpenMP 불균형(imbalance)이 여전히 높게 나타난다(Fig. 4). 이는 병렬 처리된 절점 계산 순환문의 내부 계산량이 각 절점 별 이웃 절점들의 수에 따라 달라질 수 있음에도 불구하고 순환문 작업을 균등 분배(OpenMP 기본설정)하여 부하 불균형이 발생한 것으로 파악된다. “schedule(runtime)” 설정을 해당 순환문에 추가하고, 여러 스케쥴 정책들을 테스트한 결과 “guided” 혹은 “auto” 정책에서 좋은 성능을 보였다. Fig. 5에서 추가 최적화된 코드(optimized OpenMP code 2)가 1,386초로 가장 좋은 성능을 나타내고 부하 불균형도 대폭 개선된 것을 확인할 수 있다.

Fig. 4.

Application performance analysis for the optimized OpenMP code 1

Fig. 5.

Application performance analysis for the optimized OpenMP code 2

3.2 MPI-OpenMP 혼합 병렬 처리

단일 물체에 대한 페리다이나믹 해석 코드는 비국부 영역 $H_x$를 통해 상호작용하는 절점 의존성만 존재하기 때문에 반복 계산 1회 당 절점 계산 순환문을 병렬 처리해서 성능 개선을 의도할 수 있다. 그러나 다물체 페리다이나믹 해석 코드에서는 절점과 모델 의존성이 동시에 존재하기 때문에 절점 계산 순환문을 MPI 병렬 처리해서 얻은 성능 이득과 더불어 모델 간 의존성 때문에 필요한 데이터 통신으로 인한 성능 저하가 필연적이다. 예를 들어 비국부 접촉 상호작용을 계산하기 위해서는 두 물체가 서로 상대방 절점의 속도 및 변위 정보를 알아야 한다. 가상 층간 결합 모델의 가상 절점 정보도 구속된 두 물체에 의존적이다. 따라서 모델 간 데이터 통신을 최대한 제한하기 위해 반복 계산 순환문에서 각 모델 계산은 개별 프로세스에 할당(MPI 병렬화)하고 모델 내부 절점 계산은 OpenMP 쓰레드를 이용하여 병렬 계산을 한다. 또한 반복 계산 순환문 이전의 계산은 모든 프로세스에서 동시 실행하여 통신 부하를 최대한 방지하였다. 각 프로세스에서 계산된 절점 정보(변위, 속도, 손상지수 등)는 계산을 담당하는 프로세스로 수합(reduce)되어 다음 시간 정보로 업데이트가 된 후 각 프로세스로 할당된다. 이 때 변위, 속도, 손상지수를 동시에 전달(broadcasting)하지 않고 ‘mpi_waitall’을 이용하여 전달 게이트를 열어두고 계산이 먼저 완료되는 절점 정보(본 해석 코드에서는 손상지수, 속도, 변위 순)부터 다른 계산을 기다리지 않고 논블로킹(non-blocking) 통신을 이용하여 각 통신을 계산과 동시에 처리하여 통신 효율을 더욱 높이도록 하였다.

반복 계산을 통해 얻어진 수많은 동적해석 결과는 사용자가 지정한 수준의 간격 별로 기록되어야 한다. 예를 들어, 4,000회의 반복 계산 결과를 100회마다 기록하면 총 40개 해석 결과가 시간 순서에 따라 저장된다. MPI 병렬 코드에서는 일반적으로 쓰기를 담당한 프로세스가 각 프로세스가 계산한 데이터를 통신해서 모은 후 파일 쓰기를 한다. 따라서 여러 모델에 대한 많은 양의 데이터가 발생하는 다물체 모델의 경우 데이터를 모으기 위해 요구되는 통신 부하가 많이 발생한다. 본 병렬 코드에서는 결과를 기록하는 파일을 각 프로세스에서 추가(‘append’) 방식으로 개별적으로 열고 기록하게 설정하였다. 이는 개별 프로세스가 개별 모델의 계산을 담당하여 해석 결과도 개별 프로세스 별로 구분이 되기 때문에 가능하다. 이를 통해 먼저 계산을 완료한 프로세스부터 다른 프로세스와 통신 없이 파일에 순차적으로 쓰게 만들어서 효율적으로 결과 파일을 작성할 수 있었다.

3개의 물체로 구성된 테스트 모델 1에 대해 MPI-OpenMP hybrid 병렬 처리로 개선된 해석 코드의 성능을 비교하였다(Table 1). 하나의 프로세스 당 OpenMP 쓰레드는 64개로 사용하였으며, MPI 프로세스는 모델 수인 3개가 사용되었다. OpenMP 병렬 최적화를 통해 약 1.9배의 성능향상을 얻을 수 있었고, MPI-OpenMP hybrid 병렬 처리를 통해 약 2.7배의 성능향상을 얻을 수 있었다. 테스트 모델 2는 각 1,365,912개 절점으로 구성되는 (가로, 세로, 높이)=(0.267m, 0.0534m, 0.0073m)인 소다라임 유리판 4개가 0.0006m 두께의 층간 재료(폴리우레탄)를 통해 적층되었으며, 최하단 유리판은 동일 평면 구조와 0.0058m 두께의 폴리카보네이트 탄성판(1,062,376개 절점)에 0.0006m 두께 층간재료로 적층되어 있다(Fig. 6a). 최상단 유리판은 강철 충돌체(16,400개 절점)와 고속(880m/s)으로 충돌한다. 따라서 모델 2는 6개의 물체로 구성되며 모델 1에 비해 약 35배 많은 총 6,542,424개 절점으로 구성된다. 균일한 시간 간격 2.0e-8초를 사용하였으며, 상세한 모델 설명 및 고속 충돌 파괴 현상에 대한 분석은 본 논문과 연관성이 떨어지므로 생략한다. 300회의 반복 계산을 수행한 성능을 비교하였으며, 해석 결과는 Fig. 6b와 같다. 개별 프로세스 당 OpenMP 쓰레드는 68개를 사용하고, MPI 프로세스는 6개가 사용되었다. 참고로 본 연구에서 사용한 Nurion 슈퍼컴퓨터는 2020년 현재 각 프로세스 당 최대 68개의 쓰레드를 사용할 수 있다. 최적화 전 OpenMP 코드(original)와 비교하여 MPI-OpenMP 코드를 통해 약 5.6배 성능 향상이 확인되었다(Table 2). 또한 결과 파일을 1회 기록하는데 최적화 전 OpenMP 코드는 약 130초가 소요된 반면, MPI-OpenMP 코드는 약 29초만 필요했다. 따라서 본 연구에서 적용한 MPI-OpenMP 병렬 처리는 다물체 해석에 사용된 모델의 수가 많을수록 높은 성능 향상이 기대된다.

Fig. 6.

Multi-body peridynamic simulations for Model 2

4. 결 론

본 연구에서는 다물체 충돌 파괴 시뮬레이션을 위한 페리다이나믹 해석 코드를 최적화하고 병렬 처리하였다. 핵심 계산 모듈에 대해 소개하고 의사코드를 제공하였다. 매우 많은 수의 반복 계산이 필요한 고속 충돌 해석 코드를 위해 반복 계산 순환문 내 계산을 가능한 최소화하고, OpenMP 스케쥴 정책을 조정하여 절점 간 비국부 영역 계산량이 달라서 발생하는 부하 불균형을 해결하였다. 다물체 해석 코드에서는 모델 간 의존성이 강화되어 모델 간 데이터 통신 부담이 급증하는 절점 기반 MPI 병렬화 대신 각 모델 계산을 개별 프로세스에 할당하는 모델 기반 MPI 병렬화를 수행하였다. 또한 개별 모델 내 절점 계산은 OpenMP 쓰레드를 최대한 활용하는 MPI-OpenMP hybrid 병렬 처리 구조를 선택하였다. 다물체 해석을 위해 불가피한 프로세스 간 통신 효율을 높이기 위해 논블로킹 통신을 이용하였다. 테스트 모델들을 통해 개발된 병렬 처리 해석 코드의 성능을 검토하고, 다물체 해석에 사용된 모델의 수가 많을수록 더욱 성능이 향상되는 것을 확인하였다.