Further Improvement of Direct Solution-based FETI Algorithm

Seung-Hoon Kang; DuHyun Gong; SangJoon Shin

doi:10.7734/COSEIK.2022.35.5.249

Preview

Research Paper

Journal of the Computational Structural Engineering Institute of Korea. 31 October 2022. 249-257
https://doi.org/10.7734/COSEIK.2022.35.5.249

Further Improvement of Direct Solution-based FETI Algorithm

직접해법 기반의 FETI 알고리즘의 개선

Seung-Hoon Kang¹

DuHyun Gong¹

SangJoon Shin²^*

강 승훈¹

공 두현¹

신 상준²^*

¹Graduate Student, Department of Aerospace Engineering, Seoul National University, Seoul, 08826, Korea

²Professor, Institute of Advanced Aerospace Technology, Seoul National University, Seoul, 08826, Korea

¹서울대학교 항공우주공학과 박사과정

²서울대학교 항공우주신기술연구소 교수

^{*Corresponding Author}

ABSTRACT

This paper presents an improved computational framework for the direct-solution-based finite element tearing and interconnecting (FETI) algorithm. The FETI-local algorithm is further improved herein, and localized Lagrange multipliers are used to define the interface among its subdomains. Selective inverse entry computation, using a property of the Boolean matrix, is employed for the computation of the subdomain interface stiffness and load, in which the original FETI-local algorithm requires a full matrix inverse computation of a high computational cost. In the global interface computation step, the original serial computation is replaced by a parallel multi-frontal method. The performance of the improved FETI-local algorithm was evaluated using a numerical example with 64 million degrees of freedom (DOFs). The computational time was reduced by up to 97.8% compared to that of the original algorithm. In addition, further stable and improved scalability was obtained in terms of a speed-up indicator. Furthermore, a performance comparison was conducted to evaluate the differences between the proposed algorithm and commercial software ANSYS using a large-scale computation with 432 million DOFs. Although ANSYS is superior in terms of computational time, the proposed algorithm has an advantage in terms of the speed-up increase per processor increase.

Keywords

finite element tearing and interconnecting (FETI)

localized lagrange multiplier

direct solution

multi-frontal method

selective inverse entry computation

본 논문은 직접해법 기반 FETI 알고리즘의 개선 방안을 제시하였다. 개선 대상은 FETI-local로, 해당 알고리즘은 국부 Lagrange 승수를 통해 부영역 간 경계 문제를 정의한다. 부영역 경계 강성 및 하중 계산 단계의 경우, 전체 역행렬 계산 등 과도한 비용을 요구했던 기존 알고리즘을 Boolean 행렬 특성을 활용한 선택적 역행렬 성분 계산으로 개선하였다. 전역 경계 행렬식 계산 단계의 경우, 기존 단일 프로세서 연산을 다중 프론탈 기법 기반 병렬 연산으로 대체하였다. 제시된 FETI-local 알고리즘의 성능 개선은 64만 자유도 수치 예제를 통해 검증되었으며, 기존 대비 최대 97.8%의 계산 시간 감소가 달성되었다. 또한, 기존 대비 안정적이고 개선된 확장성이 가속 지표를 통해 확인되었다. 추가로, 432만 자유도의 대용량 계산 성능 비교가 제시된 알고리즘과 상용 프로그램인 ANSYS 간에 수행되었다. 그 결과, 계산 시간 측면에선 ANSYS가 우수하였으나, 프로세서 수에 따른 가속 성능 증가율 측면에선 제시된 알고리즘이 우수한 것이 확인되었다.

키워드

유한요소 분할 및 합성법(FETI)

국부 Lagrange 승수

직접해법

다중 프론탈 기법

선택적 역행렬 성분 계산

MAIN

1. 서 론
2. 본 론
2.1 FETI-local 알고리즘
2.2 선택적 역행렬 성분 계산
2.3 다중 프론탈 기법 적용
2.4 개선된 FETI-local 알고리즘
2.5 수치예제 결과 및 고찰
3. 결 론

1. 서 론

병렬 계산 장비의 발달에 힘입어 대자유도 유한요소 구조해석에 대한 수요는 점진적으로 증가해 왔다. 그러나 유한요소 구조해석에서 통상적으로 활용되는 행렬 연산 알고리즘 중 상당수는 대자유도 해석 병렬화에 있어 한계를 보인다. 예를 들어, 직접해법(direct solution) 기반 희박행렬 연산 알고리즘의 계산 비용은 2차원 판(plate), 3차원 솔리드(solid) 요소일 경우 각각 자유도 수의 1.5, 2제곱에 비례한다(Davis, 2006). 이러한 기하급수적인 계산 비용의 증가는 병렬 대자유도 해석의 확장성(scalability)을 저해한다. 따라서 유한요소 방정식의 별도 정식화를 통한 병렬 연산 알고리즘 및 적용 방법이 제시되었다.

영역분할(domain decomposition) 기법은 전체 에너지 방정식을 각 부영역(subdomain) 에너지 및 부영역 간 경계(interface) 에너지의 합으로 정식화하며, 부영역의 독립성으로 인해 병렬 계산에 적극적으로 활용되었다. Farhat과 Roux(1991)은 비중첩(non-overlapping) 영역분할 기법의 일종인 유한요소 분할 및 합성법(FETI: finite element tearing and interconnecting)을 제시하였으며, 전역(global) Lagrange 승수를 통해 부영역 간 경계문제를 정의하였다. 이후 Park 등(1997)은 국부(local) Lagrange 승수 기반의 A(algebraically partitioned)-FETI 알고리즘을 제시하였다. 이어서 성긴 절점(coarse node) 및 전역 Lagrange 승수에 기반한 FETI-DP(dual-primal) 알고리즘이 개발되었으며(Farhat et al., 2001), 이는 현재에도 가장 각광 받는 영역분할 기법 중의 하나다.

한편, 국부 Lagrange 승수 기반 FETI 알고리즘의 개발도 지속되었다. Bauchau(2010)는 국부 Lagrange 승수 및 성긴 절점 기반의 부영역 간 경계문제에 벌칙 계수(penalty coefficient)를 추가하여 부영역 랭크(rank) 부족 해소 및 경계문제의 예조건화(preconditioning)를 도모하였다. 이는 이후 FETI-local로 명명되었으며, 직접해법 기반 알고리즘을 바탕으로 판, 쉘(shell) 유한요소 병렬 해석에 활용되었다(Gong et al., 2017; Kwak et al., 2012; 2014). 해당 문헌은 부영역별 행렬 연산을 병렬 프로세서(processor)가 분담하여 독립적으로 수행하고 경계 자유도 연산을 호스트(host) 프로세서에서 수행하는 방식을 취하였다. 그러나 이러한 접근은 3차원 솔리드 유한요소 해석으로의 확장 시 대폭 증가하는 경계 자유도로 인해 상당한 계산 비용을 초래하였다(Joo et al., 2020). 최근 Gong 등(2022)은 직접해법을 Jacobi 예조건화 켤레 기울기법(conjugate gradient method) 기반 반복해법(iterative solution)으로 변경하여 계산 시간 및 메모리의 절감을 도모했으나, 예조건 행렬의 고도화 등 필요 추가 연구가 산적해 있다.

최근 영역분할 연구 동향은 반복해법 기반 알고리즘에 집중되어 있으나, 희박행렬 인수분해(factorization) 등 직접해법의 내용을 상당 부분 차용한다. 또한, 직접해법 기반 알고리즘은 반복해법 기반 알고리즘보다 행렬 조건 수(condition number)의 제약을 적게 받는 등의 장점으로 인해 특정 상황에서 여전히 선호된다(Tak et al., 2013). FETI-local 알고리즘 고도화의 일환으로, 본 논문은 직접해법 기반 FETI-local 병렬 해석의 개선 방안을 소개한다. 선택적 역행렬 성분 계산 등 최신 직접해법 기반 희박행렬 연산 알고리즘의 기능을 적극적으로 도입하였다. 또한, 병렬 다중 프론탈(multi-frontal) 기법을 사용하여 기존 단일(serial) 프로세서 기반 경계 자유도 해석 단계를 병렬 프로세서 과정으로 변경하였다. 이후 개선 알고리즘의 성능을 확인하였다.

2. 본 론

이 장에서는 정적 구조물을 대상으로 한 FETI-local 알고리즘의 지배방정식 및 흐름도를 간단히 소개한다. 이후 선택적 역행렬 성분 계산, 병렬 다중 프론탈 기법 도입 등 직접해법 기반 FETI-local 해석 과정 개선 방안을 제시한다.

2.1 FETI-local 알고리즘

2.1.1 정적 구조 지배방정식

FETI-local 알고리즘에 의한 부영역 간 경계문제 정의는 Fig. 1과 같다. $i$ 번째 부영역 $Ω_{i}$ 의 경계 변위 $S_{i} u_{i}$ 와 이와 연관된 성긴 절점 변위 $C_{i} c$ 간에 연속조건이 정의된다. 위 두 성분은 국부 Lagrange 승수에 의한 하중 $λ_{i}$ 을 주고받는다. 이때 $u_{i}$ 는 부영역 $Ω_{i}$ 의 변위이며, $S_{i}$ 는 $u_{i}$ 내의 경계 변위를 추출하는 Boolean 행렬이다. $c$ 는 전체 성긴 절점 변위이며, $C_{i}$ 는 $c$ 중 부영역 $Ω_{i}$ 와 연관된 성분을 추출하는 Boolean 행렬이다. 에너지 방정식 서술 시, FETI-local은 연속 조건 및 국부 Lagrange 승수에 의한 에너지 성분에 벌칙 계수 $p_{i}^{c}, p_{i}^{λ}$ 를 각각 부여한다. 이를 기반으로 한 부영역 전체 퍼텐셜 에너지 $Π_{i}$ 의 합 $Π$ 은 아래와 같다(Bauchau, 2010; Gong et al., 2022).

(1)

Π = \sum_{i = 1}^{n^{s}} Π_{i}

(2)

Π_{i} = \frac{1}{2} u_{i}^{T} K_{i} u_{i} + \frac{1}{2} p_{i}^{c} (S_{i} u_{i} - C_{i} c)^{T} (S_{i} u_{i} - C_{i} c) + p_{i}^{λ} λ_{i}^{T} (S_{i} u_{i} - C_{i} c) - u_{i}^{T} f_{i}

$n^{s}$ 는 전체 부영역 수이며 $f_{i}$ 는 부영역 $Ω_{i}$ 에 작용하는 외력 하중이다. 식 (1)에 최소 퍼텐셜 에너지 원리를 적용하여 아래 식 (3)과 같은 정적 구조 평형 방정식을 도출한다. 본 논문에선 벌칙 계수에 대해 $p_{i}^{c} = p_{i}^{λ} = p_{i}$ 임을 가정하였다.

(3)

[\begin{matrix} {\hat{K}}_{1} & - {\hat{S}}_{1}^{T} C_{1} \\ ⋱ & ⋮ \\ {\hat{K}}_{n^{s}} & - {\hat{S}}_{n^{s}}^{T} C_{n^{s}} \\ - C_{1}^{T} {\hat{S}}_{1} & \dots & - C_{n^{s}}^{T} {\hat{S}}_{n^{s}} & \sum_{i = 1}^{n^{s}} p_{i} C_{i}^{T} C_{i} \end{matrix}] \{\begin{matrix} {\hat{u}}_{1} \\ ⋮ \\ {\hat{u}}_{n^{s}} \\ c \end{matrix}\} = \{\begin{matrix} {\hat{f}}_{1} \\ ⋮ \\ {\hat{f}}_{n^{s}} \\ 0 \end{matrix}\}

(4)

{\hat{K}}_{i} = [\begin{matrix} K_{i} + p_{i} S_{i}^{T} S_{i} p_{i} S_{i}^{T} \\ p_{i} S_{i} 0 \end{matrix}]

이때, ${\hat{S}}_{i} = p_{i} [S_{i} I_{n_{i}^{λ}}], {\hat{u}}_{i}^{T} = \{u_{i}^{T} λ_{i}^{T}\}, {\hat{f}}_{i}^{T} = \{f_{i}^{T} 0\}$ 이다. $I_{n_{i}^{λ}}$ 는 크기 $n_{i}^{λ}$ 의 단위 정방행렬이며, $n_{i}^{λ}$ 는 부영역 $Ω_{i}$ 와 연결된 성긴 절점 자유도 수이다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-05/N0040350501/images/Figure_jcoseik_35_05_01_F1.jpg

Fig. 1.

Schematic description of the interface problem for the present FETI-local

2.1.2 경계 자유도 연산(Step I, II)

블록(block) 행렬 연산을 통해 식 (3)을 성긴 절점 변위 $c$ 에 대해 정리하면 아래와 같다.

(5)

\sum_{i = 1}^{n^{s}} (C_{i}^{T} K_{i}^{c} C_{i}) c = \sum_{i = 1}^{n^{s}} (C_{i}^{T} f_{i}^{c})

식 (5)는 FETI-local 알고리즘의 경계문제 정의이다. 이때, $K_{i}^{c} = p_{i} I_{n_{i}^{λ}} - {\hat{S}}_{i} {\hat{K}}_{i}^{- 1} {\hat{S}}_{i}^{T}$ 이며 $f_{i}^{c} = {\hat{S}}_{i} {\hat{K}}_{i}^{- 1} {\hat{f}}_{i}$ 이다. 이는 요소 강성행렬 $K_{i}^{c}$ 및 요소 하중 벡터 $f_{i}^{c}$ 을 갖는 $n^{s}$ 개 초요소(super- element)의 정적 구조 평형 방정식과 같다. 이때, Boolean 행렬 $C_{i}$ 는 연결 행렬(connectivity matrix)의 역할을 한다.

식 (5)의 $c$ 를 직접해법에 기반하여 계산할 경우, 부영역별 경계 강성행렬 $K_{i}^{c}$ 을 도출하는 선행 과정이 필요하다. 이는 희박행렬 ${\hat{K}}_{i}$ 의 인수분해 및 역행렬 성분 계산을 요구하며, 상당한 비용을 요구한다. 부영역별 경계 하중 벡터 $f_{i}^{c}$ 의 경우 ${\hat{K}}_{i}$ 의 인수분해 결과를 활용할 경우 비교적 수월하게 계산할 수 있다. 이러한 $K_{i}^{c}$ 및 $f_{i}^{c}$ 의 계산 과정을 본 논문에선 Step I이라 명명하였다.

이후 부영역별 경계 강성행렬 및 하중 벡터를 조립(assemblage)하여 성긴 절점 변위 $c$ 를 계산하며, 이 과정을 Step II이라 명명하였다. 앞서 언급하였듯, 판, 쉘 요소 기반의 영역분할 해석을 동일 전역 자유도 수를 지니는 3차원 솔리드 요소 기반 해석으로 확장할 경우 경계 자유도 수는 대폭 증가한다(Joo et al., 2020).

2.1.3 부영역별 자유도 연산(Step III)

전역 성긴 절점 변위 $c$ 를 계산한 이후엔, 부영역별 변위 $u_{i}$ 를 부영역별 성긴 절점 변위 $C_{i} c$ 를 통해 독립적으로 계산한다. 식 (3)에 따르면, 부영역 $Ω_{i}$ 내 변위 계산 과정은 아래와 같으며, 본 논문에선 이 과정을 Step III이라 명명하였다.

(6)

{\hat{K}}_{i} {\hat{u}}_{i} = {\hat{f}}_{i} + {\hat{S}}_{i}^{T} C_{i} c

앞서 Step I에서 행렬 ${\hat{K}}_{i}$ 의 인수분해를 수행했으므로, 이를 활용할 경우 Step III의 계산 비용은 다른 과정에 비해 미미하다.

2.2 선택적 역행렬 성분 계산

식 (5)에서 확인할 수 있듯, 직접해법 기반 FETI-local 알고리즘은 부영역별로 ${\hat{S}}_{i} {\hat{K}}_{i}^{- 1} {\hat{S}}_{i}^{T}$ 을 계산하는 것을 요구한다. 그러나, 해당 행렬의 도출을 위해 강성행렬의 역행렬 ${\hat{K}}_{i}^{- 1}$ 내 모든 성분을 계산하는 것은 상당히 비효율적이다. 행렬 ${\hat{S}}_{i}$ 가 Boolean 행렬 $S_{i}$ 및 단위행렬 $I_{n_{i}^{λ}}$ 로 구성된 점을 활용하면 역행렬 ${\hat{K}}_{i}^{- 1}$ 중 필요 성분을 파악할 수 있다.

${\hat{K}}_{i}^{- 1}$ 가 대칭행렬일 경우, 이를 내부 및 경계 자유도 성분으로 나누면 아래와 같다.

(7)

{\hat{K}}_{i}^{- 1} = [\begin{matrix} A_{i} & B_{i}^{T} \\ B_{i} & H_{i} \end{matrix}], A_{i} ∊ ℝ^{n_{i}^{u} \times n_{i}^{u}}, B_{i} ∊ ℝ^{n_{i}^{λ} \times n_{i}^{u}}, H_{i} ∊ ℝ^{n_{i}^{λ} \times n_{i}^{λ}}

$n_{i}^{u}$ 는 부영역 $Ω_{i}$ 의 내부 변위 자유도 수이다. 식 (6)을 통해 행렬 ${\hat{S}}_{i} {\hat{K}}_{i}^{- 1} {\hat{S}}_{i}^{T}$ 을 표현하면 아래와 같다.

(8)

{\hat{S}}_{i} {\hat{K}}_{i}^{- 1} {\hat{S}}_{i}^{T} = p_{i}^{2} [S_{i} I_{n_{i}^{λ}}] [\begin{matrix} A_{i} & B_{i}^{T} \\ B_{i} & H_{i} \end{matrix}] [\begin{matrix} S_{i}^{T} \\ I_{n_{i}^{λ}} \end{matrix}] = p_{i}^{2} (S_{i} A_{i} S_{i}^{T} + B_{i}^{T} S_{i}^{T} + S_{i} B_{i} + H_{i})

Boolean 행렬 $S_{i}$ 는 대상 행렬 앞에서 작용할 경우( $S_{i} X$ ), 행렬( $X$ )의 행(row) 중 성긴 절점 자유도와 연관된 행을 추출한다. 뒤에서 작용할 경우( $Y S_{i}^{T}$ ), 행렬( $Y$ )의 열(column)을 마찬가지로 추출한다. Fig. 2는 $n_{i}^{λ} = 10, n_{i}^{λ} = 3$ 일 때 Boolean 행렬 연산 예시이다. 각 곱셈 연산 시 추출되는 행렬 성분 중 상삼각(upper-triangle) 성분이 표시되어 있다. 예시의 경우 ${\hat{K}}_{i}^{- 1}$ 의 169개 성분 중 21개 성분만으로 행렬 ${\hat{S}}_{i} {\hat{K}}_{i}^{- 1} {\hat{S}}_{i}^{T}$ 를 도출한다. 일반적인 경우, 총 $(n_{i}^{u} + n_{i}^{λ})^{2}$ 개 성분 중 $2 (n_{i}^{λ})^{2} + n_{i}^{λ}$ 개 성분만을 통해 행렬 ${\hat{S}}_{i} {\hat{K}}_{i}^{- 1} {\hat{S}}_{i}^{T}$ 의 계산이 가능하다. 따라서, 본 논문에선 ${\hat{K}}_{i}^{- 1}$ 내 $2 (n_{i}^{λ})^{2} + n_{i}^{λ}$ 개 성분을 선택적으로 계산하여 ${\hat{S}}_{i} {\hat{K}}_{i}^{- 1} {\hat{S}}_{i}^{T}$ 을 도출하는 방법으로 Step I의 개선을 도모하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-05/N0040350501/images/Figure_jcoseik_35_05_01_F2.jpg

Fig. 2.

Examples of Boolean matrix multiplication

최근 PARDISO(Bollhöfer et al., 2019), MUMPS(Amestoy et al., 2019) 등 많은 수의 직접해법 기반 희박행렬 라이브러리가 선택적 역행렬 성분 계산 기능을 지원하며, 이에 관한 연구가 지속되고 있다. 본 논문은 MUMPS 5.5.0을 도입하여 해당 기능을 구현하였다.

Step I에서 MUMPS는 단일 프로세서별로 실행되었으며, 우선 행렬 ${\hat{K}}_{i}$ 의 인수분해를 수행하였다. 이때, 희박행렬 ${\hat{K}}_{i}$ 는 성분별로 행렬 색인을 지니는 좌표 형식(coordinate format)으로 저장되었다. 이후 인수분해 결과를 바탕으로 필요 역행렬 성분을 순차적으로 계산하였고 부영역 경계 강성행렬 $K_{i}^{c}$ 을 구성하였다. 인수분해 결과는 Step I의 부영역 경계 하중 벡터 $f_{i}^{c}$ 의 계산 및 Step III의 최종 변위 계산에서 또한 활용되었다. 자세한 과정은 2.4절에 설명되어 있다.

2.3 다중 프론탈 기법 적용

직접해법 기반 FETI-local 알고리즘의 Step II는 부영역 경계 강성행렬 및 하중 벡터를 바탕으로 전역 경계 강성행렬 및 하중 벡터를 조립하고 성긴 절점 변위 $c$ 를 계산하는 과정이다. 기존 Step II는 각 프로세서에서 계산한 부영역 정보를 호스트 프로세서로 전송한 뒤, 해당 프로세서에서 조립 및 성긴 절점 변위 계산을 수행하는 방식이었다(Gong et al., 2017; Joo et al., 2020; Kwak et al., 2012; 2014). 그러나 이러한 방식은 호스트 프로세서의 계산량 및 메모리 과다를 유발한다. 본 논문에선 단일 프로세서의 계산량 집중 완화 및 Step II의 병렬화를 위해 MUMPS의 다중 프론탈 행렬 연산 기능을 도입하였다.

다중 프론탈 기법 적용 시 분산 저장 형식을 사용하였다. 우선 Step I에서 계산한 부영역 경계 강성행렬 및 하중 벡터의 행렬 및 벡터 색인을 생성하였다. 이때 색인은 Step I과 같이 좌표 형식이지만 전역 경계 행렬식 (5)을 기준으로 하였다. 분산 저장 형식은 동일한 행렬 혹은 벡터 색인이 발생하면 이를 합 연산으로 취급한다. 따라서 호스트 프로세서에 부영역 정보를 전송하지 않고 병렬 프로세서상에서 전역 경계 강성행렬 및 하중 벡터 조립을 구현하는 것이 가능하다. 이러한 분산 저장 형식은 호스트 프로세서의 메모리 요구량을 절감하며, MUMPS뿐만 아니라 PARDISO 등 많은 수의 희박행렬 라이브러리에서 지원한다. Fig. 3는 2개의 프로세서가 각각 2개의 부영역 경계 강성행렬 및 하중 벡터를 지닐 때 전역 경계 강성행렬 및 하중 벡터 분산 저장 예시를 설명한다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-05/N0040350501/images/Figure_jcoseik_35_05_01_F3.jpg

Fig. 3.

Examples of the distributed storage format

전역 경계 강성 및 하중의 분산 저장 이후 MUMPS를 통해 식 (5)를 병렬 프로세서상에서 계산하였다. 이때, 그래프 분할 라이브러리인 METIS 5.1.0(Karypis and Kumar, 1998)를 별도로 연동하여 MUMPS의 행렬 성분 분석 단계에서 사용하였다. 계산된 성긴 절점 변위 $c$ 는 각 부영역에 $C_{i} c$ 꼴로 할당되어 Step III의 부영역별 변위 계산에 사용되었다.

2.4 개선된 FETI-local 알고리즘

앞선 내용을 바탕으로 개선된 FETI-local 알고리즘의 흐름도는 Fig. 4와 같다. 우선, 유한요소로 이산화된 해석 대상의 영역을 METIS를 이용해 분할하였다. Fig. 4의 경우 부영역 개수는 16개이다. 이후 부영역 정보를 프로세서에 균일하게 할당하였다. Fig. 4는 프로세서가 8개이며, 각 프로세서당 2개의 부영역이 할당된 경우이다.

이후 Step I에 진입하였다. Step I은 프로세서 간 통신 없이 독립적으로 이루어졌다. 우선 부영역 행렬 ${\hat{K}}_{i}$ 을 인수분해하였다. 이후 2.2절의 내용과 같이 선택적 역행렬 성분 계산을 수행하여 ${\hat{S}}_{i} {\hat{K}}_{i}^{- 1} {\hat{S}}_{i}^{T}$ 을 도출하였다. 부영역 경계 강성행렬 $K_{i}^{c}$ 는 $p_{i} I_{n_{i}^{λ}} - {\hat{S}}_{i} {\hat{K}}_{i}^{- 1} {\hat{S}}_{i}^{T}$ 을 연산함으로써 계산하였다. 부영역 경계 하중 벡터 $f_{i}^{c}$ 는 앞선 ${\hat{K}}_{i}$ 의 인수분해 결과를 바탕으로 방정식 ${\hat{K}}_{i} x = {\hat{f}}_{i}$ 을 계산한 후, $f_{i}^{c} = {\hat{S}}_{i} x$ 와 같이 도출하였다. Step I의 인수분해, 선택적 역행렬 성분 계산, $f_{i}^{c}$ 계산에 필요한 행렬 연산은 각 프로세서에서 단일 MUMPS를 통해 수행되었다. 이때 단일 MUMPS의 실행은 개별 MPI(message passing interface) 커뮤니케이터(MPI_COMM_SELF)를 입력하여 구현되었다.

Step II에선 Fig. 3과 같이 $K_{i}^{c}$ 및 $f_{i}^{c}$ 의 분산 저장 양식 색인을 생성함으로써 전역 경계 방정식 (5)에 필요한 항을 구축하였다. 이후 전역 경계 강성행렬 $\sum_{i = 1}^{n^{s}} (C_{i}^{T} K_{i}^{c} C_{i})$ 의 인수분해를 수행하고, 성긴 절점 변위 $c$ 를 계산하였다. Step II의 인수분해 및 성긴 절점 변위 $c$ 의 계산은 병렬 MUMPS를 통해 수행되었다. 이때, 병렬 MUMPS의 실행은 전역 MPI 커뮤니케이터(MPI_ COMM_WORLD)를 입력하여 구현되었다. 즉, Step II에서 프로세서 간 통신이 발생하였다.

Step III에서 프로세서별로 할당된 부영역 성긴 절점 변위 $C_{i} c$ 를 이용해 부영역 정적 구조 방정식 (6)을 계산하였다. 이때, 앞서 Step I에서 수행한 행렬 ${\hat{K}}_{i}$ 의 인수분해 결과를 활용하였다.

2.5 수치예제 결과 및 고찰

이 장에서는 기존 직접해법 기반 FETI-local 알고리즘 대비 본 논문에서 제시한 알고리즘의 계산 성능 개선을 검증한다. 또한, 대자유도 구조 해석 예제를 통해 알고리즘의 적용성을 확인한다.

성능 검증을 위한 프로그램은 Window 10 운영체제상에서 FORTRAN90 언어 기반으로 구현되었다. 실수 및 정수 자료형의 경우 각각 배수 정밀도(64bit) 실수 및 32bit 정수가 사용되었다. 계산 사양은 데스크톱 컴퓨터 사양으로, Intel® Core™ i9-10900KF 3.70GHz 프로세서 및 4개의 삼성 DDR4 PC4- 25600 32GB 메모리가 장착되었다. MPI 통신은 Intel MPI를 통해 구현되었다.

2.5.1 계산 성능 검증

기존 및 개선 알고리즘의 동등한 비교를 위해, 기존 알고리즘의 프로그램은 선택적 역행렬 성분 계산 및 Step II 병렬화 외에 개선 알고리즘과 동일하게 구성하였다. 즉, Fig. 4의 선택적 역행렬 성분 계산은 기존의 전체 역행렬 성분 계산으로 대체되었다. Step II는 호스트 프로세서에 전역 경계 행렬식 (5)의 정보를 전송하고 해당 프로세서 내에서 연산을 수행하는 방식으로 대체되었다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-05/N0040350501/images/Figure_jcoseik_35_05_01_F4.jpg

Fig. 4.

Framework of the improved FETI-local algorithm

Fig. 5와 같이 균열이 있는 탄소 섬유/에폭시 시편 예제를 대상으로 FETI-local 알고리즘의 계산 성능을 확인하였다. 시편의 한쪽 끝은 고정단이며, 반대쪽엔 인장 변위 경계조건이 가해졌다. 또한, 옆면에는 1/4 대칭 해석을 위한 변위 경계조건이 부여되었다. 유한요소 이산화는 약 64만 개의 변위 자유도 및 약 118만 개의 사면체 솔리드 요소로 구성되었다. 예제에 대한 자세한 설명은 참고문헌(Gong et al. 2022)에 기술되어 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-05/N0040350501/images/Figure_jcoseik_35_05_01_F5.jpg

Fig. 5.

Numerical example: carbon fiber/epoxy specimen (longitudinal stress distribution, All units in Pascal)

Fig. 6과 같이 부영역 수 64, 128개인 경우에 대해 계산이 수행되었으며, 각 경우의 성긴 절점 자유도 개수는 각각 약 8만, 11만 개이다. 프로세서 수를 1개부터 8개까지 변경하여 계산 시간을 확인하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-05/N0040350501/images/Figure_jcoseik_35_05_01_F6.jpg

Fig. 6.

Carbon fiber/epoxy specimen: finite element mesh partitioning by 64 (a) and 128 (b) subdomains

Table 1은 Step I부터 III까지의 총 계산 시간을 부영역 수, 프로세서 수별로 정리한 것이다. 본 논문에서 개선한 FETI-local 알고리즘은 기존 대비 최대 97.8%의 계산 시간 감소를 달성함을 확인하였다. 전반적으론 97% 내외의 계산 시간 감소가 이루어졌다.

Table 1.

Carbon fiber/epoxy specimen: total computational time and reduction

Number of subdomains	Number of processors	Total computational time
Number of subdomains	Number of processors	Present(s)	Original(s)	Reduction(%)
64	1	275.5	10,322.5	97.3
64	2	162.0	7,283.4	97.8
64	4	98.0	3,950.5	97.5
64	8	61.3	2,288.3	97.3
128	1	267.1	7,066.0	96.2
128	2	160.7	5,263.1	96.9
128	4	92.0	2,820.7	96.7
128	8	60.7	1,624.9	96.3

Fig. 7은 총 계산 시간 중 각 단계에 소요되는 시간을 분류한 것이다. 알고리즘 개선에 따른 계산 속도 향상이 가장 많이 이루어진 과정은 Step I으로, 이는 선택적 역행렬 성분 계산의 성능을 방증한다. 기존 알고리즘의 Step II는 단일 프로세서 연산이므로, 프로세서 수와 관계없이 일정한 계산 시간을 보였다. 반면에. 본 논문에서 제시한 FETI-local 알고리즘의 Step II는 다중 프론탈 기법 기반 병렬 프로세서 과정이므로, 프로세서 수에 따라 계산 시간이 감소하였다. 그러나 Step I 및 III와는 달리 Step II는 프로세서 간 통신이 발생하는 과정이므로, 프로세서 수가 증가함에 따라 계산 속도 향상 정도가 감소함을 확인하였다. Step III의 경우 기존 및 본 논문에서 제시한 FETI-local 알고리즘상에서 동일하므로, 유사한 계산 소요시간을 보였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-05/N0040350501/images/Figure_jcoseik_35_05_01_F7.jpg

Fig. 7.

Carbon fiber/epoxy specimen: computational time for each step (The values in the parenthesis indicates the number of the subdomains)

부영역 수에 따른 계산 시간의 경우, 부영역 수 증가에 따른 성긴 절점 변위 및 전역 경계 행렬식 (5)의 자유도 증가로 인해 Step II 계산 시간이 증가하였다. 반면 Step I은 부영역 수 증가에 따라 계산 시간이 감소하는 양상을 보였다. 이는 프로세서별 부영역 수 증가에 따른 계산 비용 증가보다 부영역별 역행렬 계산 비용 감소의 영향이 크다는 것을 의미한다.

계산 시간 감소에 이어 확장성의 개선 여부를 확인하였다. 확장성은 프로세서 수 증가에 따른 계산 시간의 감소 능력을 의미하며, 가속 지표(speed-up)을 통해 확인할 수 있다. 이때 가속 지표는 단일 프로세서 계산 시간과 병렬 프로세서 계산 시간 간의 비율로 정의된다. 사용자는 프로그램의 확장성을 통해 프로세서 추가 도입 시 계산 시간의 감소 정도를 예상할 수 있다.

Fig. 8은 가속 지표를 부영역 수 및 프로세서 수별로 정리한 것이다. 본 논문에서 개선한 FETI-local 알고리즘은 기존 대비 최대 23.8%의 가속 지표 상승을 달성함을 확인하였다. 부영역 수에 따른 유의미한 가속 지표의 변화는 확인되지 않았다. Table 2는 프로세서 수 증가에 따른 가속 지표 증가 양상을 정리한 것이다. 기존 알고리즘의 경우 프로세서 수 증가에 따른 가속 지표 증가 양상이 불규칙적인데 반하여, 본 논문에서 개선한 알고리즘은 비교적 안정적인 가속 지표 증가 양상을 보임을 확인하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-05/N0040350501/images/Figure_jcoseik_35_05_01_F8.jpg

Fig. 8.

Carbon fiber/epoxy specimen: comparison of the scalability (The values in the parenthesis indicates the number of the subdomains)

Table 2.

Carbon fiber/epoxy specimen: average speed-up and its increment in terms of the number of the processors

Number of processors	Average speed-up (increment)
Number of processors	Present	Original
1	1.00	1.00
2	1.68(68.1%)	1.38(38.0%)
4	2.86(70.0%)	2.56(85.4%)
8	4.45(55.6%)	4.42(73.1%)

본 논문 예제의 계산에 사용된 프로세서 수 내에선 Step I의 계산 시간이 Step II에 비교하여 많은 비중을 차지하였다. Fig. 7에서 확인할 수 있듯, 프로세서 수의 증가에 따라 Step II보다 Step I에서 많은 계산 시간 감소가 발생할 것을 기대할 수 있다. 따라서, 일정 수준 이상의 프로세서 수에선 Step II의 소요 시간이 Step I보다 많을 것이다. 또한, 부영역 수가 증가함에 따라 Step II의 계산 비용이 증가할 것이다. 이러한 점은 대량의 프로세서 및 부영역 해석에서 확장성 양상의 변화를 야기할 것이며, 추후 수십 개 이상의 프로세서를 지원하는 클러스터 컴퓨팅 환경에서 확인될 예정이다.

2.5.2 대자유도 예제 적용

본 논문에서 제시한 알고리즘을 대자유도 해석 예제에 적용하였다. 대상으로 Fig. 9와 같은 형상 및 Table 3의 물성치를 지닌 평직 탄소 섬유 다발/에폭시 인장 시편 예제가 제작되었다. 이때, $E, G, ν$ 는 각각 탄성/전단 계수 및 Poisson 비이며, 아래 첨자는 이방성 재료 방향을 의미한다. $L, w, t$ 는 각각 시편의 길이, 너비, 두께이다. 에폭시의 물성치 및 평직 섬유 다발의 기하 정보는 Abot 등(2011)이 기술한 내용에 기반하였다. 동일 참고문헌의 탄소 섬유 및 에폭시 정보를 기반으로 Gong 등(2022)이 대표 체적요소(representative volume element) 해석을 수행했던 바 있으며, 본 예제의 섬유 다발 물성치는 이를 참고하여 구성되었다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-05/N0040350501/images/Figure_jcoseik_35_05_01_F9.jpg

Fig. 9.

Computer-aided drawing of the plain woven fabric specimen (a) and its fiber bundle (b)

Table 3.

Material properties for the tensile specimen of the plain woven fabric

Material	Property	Value
Fiber bundle (Gong et al., 2022)	$E_{11}$	200GPa
	$E_{22}, E_{33}$	46.5GPa
	$G_{12}, G_{13}$	22.5GPa
	$G_{23}$	20.5GPa
	$ν_{12}, ν_{13}$	0.294
	$ν_{23}$	0.355
Epoxy matrix (Abot et al., 2011)	$E$	2.9GPa
Epoxy matrix (Abot et al., 2011)	$ν$	0.37

인장 시편의 한쪽 끝은 고정단이며, 반대쪽엔 인장 변위 경계조건이 가해졌다. 유한요소 이산화는 약 432만 개의 변위 자유도 및 약 784만 개의 사면체 솔리드 요소로 구성되었다. 유한요소망은 Fig. 10과 같이 512개 부영역으로 분할되었으며. 부영역 간 경계는 약 54만 자유도의 성긴 절점으로 구성되었다. 앞선 예제와 마찬가지로 사용 프로세서 수가 1개부터 8개까지인 경우에 대해 계산 성능을 확인하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-05/N0040350501/images/Figure_jcoseik_35_05_01_F10.jpg

Fig. 10.

Plain woven fabric specimen: finite element mesh partitioning by 512 subdomains

또한, 동일한 예제의 해석이 상용 프로그램인 ANSYS Mechanical 2021 R1(Anonymous, 2021)을 통해 진행되었다. 이 경우 ANSYS 내 분산 메모리 병렬 연산 및 희박행렬 직접해법 기능이 사용되었으며, FETI-local 알고리즘 해석과 동일한 장비 상에서 계산이 수행되었다.

Fig. 11은 본 논문의 알고리즘 및 ANSYS를 통한 해석의 계산 성능을 비교한 것이다. 본 논문에서 제시된 내용을 통해 FETI-local 알고리즘의 계산 속도 개선이 상당히 이루어졌으나, ANSYS의 계산 속도엔 미치지 못하는 것으로 확인되었다. 사용된 프로세서 수를 기준으로, ANSYS가 약 1.9에서 2.5배 신속한 계산을 수행함을 확인하였다. 이는 상용 프로그램인 ANSYS가 사용하는 희박행렬 라이브러리와 무료 라이브러리인 MUMPS의 성능 차이에 기인할 수 있다. 그러나 프로세서 수가 4개에서 8개로 증가할 경우, ANSYS의 가속 지표 증가율은 16.8%로 현저히 낮은 것으로 확인되었다. 반면 본 논문의 FETI-local 알고리즘의 경우, 동일 프로세서 기준 48.1%의 가속 성능 증가율이 확인되었다. 이는 수십 개 이상의 다중 프로세서 연산에서 본 논문의 FETI-local 알고리즘이 ANSYS에 비해 뛰어난 계산 속도를 발현할 것임을 기대할 수 있는 근거이며, 추후 클러스터 컴퓨팅 환경에서 확인될 예정이다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-05/N0040350501/images/Figure_jcoseik_35_05_01_F11.jpg

Fig. 11.

Plain woven fabric specimen: comparison of the computational time (a) and scalability (b)

또한, 본 논문의 알고리즘과 ANSYS는 동일한 구조 거동을 도출함이 확인되었다. 그중, 시편의 길이 방향 탄성 계수( $E_{11}$ )가 Abot 등(2011)의 인장 실험 결과와 비교되었다. 이때 유한요소 해석의 탄성 계수는 고정단의 길이 방향 반력(reaction force) 계산 결과 $R_{1}$ , 인장 경계조건의 변위 $Δ d$ 를 토대로 식 (9)와 같이 계산되었다.

(9)

E_{11} = \frac{R_{1} L}{w t Δ d}

Table 4에서 확인할 수 있듯, 유한요소 해석의 탄성 계수 계산 결과는 실험 대비 10.3%의 차이를 보였다. 이는 제시된 FETI-local 알고리즘이 복합재료 실험을 보완 및 대체하여 실험 비용을 절감할 수 있음을 보인다. 이처럼 FETI-local 알고리즘은 대자유도를 취급하는 복합재료 상세 다중 스케일(multi- scale) 해석 분야에서 활용 가능할 것으로 전망된다.

Table 4.

Plain woven fabric specimen: uniaxial elastic modulus

Finite element analysis		Experiment (Abot et al., 2011)	Discrepancy
Present	ANSYS
70.7GPa	70.7GPa	64.1GPa	10.3%

3. 결 론

본 논문은 직접해법 기반 FETI-local 알고리즘의 개선 방안을 제시하였다. 직접해법 기반 FETI-local 알고리즘은 부영역별 경계 강성행렬 및 하중벡터 계산 단계(Step I), 전역 경계 행렬식 계산 단계(Step II), 부영역별 변위 계산 단계(Step III)로 구성되었다. 이 중 Step I에서 요구하는 역행렬 계산 과정을 선택적 역행렬 성분 계산을 도입하여 개선하였다. 또한, 기존 단일 프로세서 연산으로 수행되었던 Step II를 다중 프론탈 기법 기반 병렬 프로세서 연산으로 대체하였다. 해당 기법들은 직접해법 기반 희박행렬 연산 라이브러리인 MUMPS를 통해 구현되었다.

제시된 FETI-local 알고리즘의 계산 성능 개선 여부는 3차원 이방성 탄소 섬유/에폭시 시편 수치 예제를 통해 검증되었다. 제시된 알고리즘을 도입할 경우 기존 대비 최대 97.8%의 계산 시간 절감이 달성됨을 확인하였다. 계산 단계별 소요 시간을 확인한 결과, Step I의 계산 시간이 확연하게 감소하였다. 기존 알고리즘에서 프로세서 수와 관계없이 일정한 계산 시간을 보인 Step II는 개선된 알고리즘의 병렬화에 따라 계산 시간 절감이 달성되었다. 확장성의 경우, 제시된 알고리즘은 기존 대비 최대 23.8%가 상승한 가속 지표를 보였다. 또한, 프로세서 수에 따라 안정적인 가속 지표 증가 양상을 확인하였다.

또한, 평직 탄소 섬유 다발/에폭시 인장 시편 수치 예제를 통해 대자유도 계산 성능을 ANSYS와 비교하였다. 사용 프로세서 수가 1개부터 8개까지일 때 계산 속도는 ANSYS가 본 논문의 FETI-local 알고리즘보다 우수하였다. 그러나 본 논문의 알고리즘은 프로세서 증가에 따른 가속 성능 증가율 측면에서 우수하여, 향후 클러스터 컴퓨팅 환경에서의 계산 시 유리할 것이 예상되었다. 시편 길이 방향 탄성 계수 해석 결과는 본 논문의 알고리즘과 ANSYS의 결과가 동일하였으며, 실험 결과와는 10% 이내의 차이가 발생하였다.

본 논문의 계산 알고리즘 검증은 데스크톱 사양 컴퓨터상에서 수행되었으며, 향후 수십 개 이상의 프로세서를 지원하는 클러스터 컴퓨팅 환경에서의 검증이 진행될 예정이다. 또한, 제시된 직접해법 기반 FETI-local 알고리즘의 개선기법은 향후 반복해법 기반 FETI-local 알고리즘의 고도화에 사용될 예정이다.

Acknowledgements

본 논문은 2022년도 정부(과학기술통신부)의 재원으로 한국연구재단의 지원을 받아 수행된 연구입니다(No. 2021R1A5 A1031868, 2021R1A2C1007352).

References

Abot, J.L., Gabbai, R.D., Harsley, K. (2011) Effect of Woven Fabric Architecture on Interlaminar Mechanical Response of Composite Materials: An Experimental Study, J. Reinf. Plast. & Compos., 30(24), pp.2003~2014. 10.1177/0731684411431118

Amestoy, P.R., Buttari, A., L'Excellent, J.Y., Mary, T. (2019) Performance and Scalability of the Block Low-Rank Multifrontal Factorization on Multicore Architectures, ACM Trans. Math. Softw., 45(1), pp.1~26. 10.1145/3242094

Anonymous (2021) ANSYS Mechanical User's Guide; Release 2021 R1, ANSYS Inc., Canonsburg.

Bauchau, O.A. (2010) Parallel Computation Approaches for Flexible Multibody Dynamics Simulations, J. Frankl. Inst., 347(1), pp.53~68. 10.1016/j.jfranklin.2009.10.001

Bollhöfer, M., Eftekhari, A., Scheidegger, S., Schenk, O. (2019) Large-scale Sparse Inverse Covariance Matrix Estimation, SIAM J. Sci. Comput., 41(1), pp.A380~401. 10.1137/17M1147615

Davis, T.A. (2006) Direct Methods for Sparse Linear Systems, Society for Industrial and Applied Mathematics, Philadelphia, p.130. 10.1137/1.9780898718881

Farhat, C., Lesoinne, M., LeTallec, P., Pierson, K., Rixen, D. (2001) FETI-DP: A Dual-primal Unified FETI Method-Part I: A Faster Alternative to the Two-level FETI Method, Int. J. Numer. Methods Eng., 50(7), pp.1523~1544. 10.1002/nme.76

Farhat, C., Roux, F.-X. (1991) A Method of Finite Element Tearing and Interconnecting and its Parallel Solution Algorithm, Int. J. Numer. Methods Eng., 32(6), pp.1205~1227. 10.1002/nme.1620320604

Gong, D.H., Joo, H.S., Shin, S.J. (2017) Three-dimensional Large-size Structural Analysis using FETI-local Method, COSEIK Academic Symposium, Computational Structural Engineering Institute of Korea, p.3.

Gong, D.H., Kang, S., Joo, H.S., Cho, H., Chun, T., Shin, S.J. (2022) Detailed Structural Analysis for Fiber-reinforced Polymer with Singularities via FETI Domain Decomposition, J. Compos. Mater., 56(11), pp.1753~1764. 10.1177/00219983221074511

Joo, H.S., Gong, D.H., Kang, S., Chun, T., Shin, S.J. (2020) Parallel Computational Approach for Three-dimensional Solid Element using Extra Shape Function based on Domain Decomposition Approach, J. Korean Soc. Indus. & Appl. Mech., 24(2), pp.199~214.

Karypis, G., Kumar, V. (1998) A Fast and High Quality Multilevel Scheme for Partitioning Irregular Graphs, SIAM J. Sci. Comput., 20(1), pp.359~392. 10.1137/S1064827595287997

Kwak, J.Y., Cho, H., Shin, S.J., Bauchau, O.A. (2012) Development of Finite Element Domain Decomposition Method using Local and Mixed Lagrange Multipliers, J. Comput. Struct. Eng. Inst. Korea, 25(6), pp.469~476. 10.7734/COSEIK.2012.25.6.469

Kwak, J.Y., Chun, T., Shin, S.J., Bauchau, O.A. (2014) Domain Decomposition Approach to Flexible Multibody Dynamics, Comput. Mech., 53, pp.147~158. 10.1007/s00466-013-0898-8

Park, K.C., Justino, Jr., M.R., Felippa, C.A. (1997) An Algebraically Partitioned FETI Method for Parallel Structural Analysis: Algorithm Description, Int. J. Numer. Methods Eng., 40(15), pp.2717~2737. 10.1002/(SICI)1097-0207(19970815)40:15<2717::AID-NME185>3.0.CO;2-B

Tak, M., Song, S., Jeon, H.K., Park, T. (2013) Accuracy Analysis of Parallel Method based on Non-overlapping Domain Decomposition Method, J. Comput. Struct. Eng. Inst. Korea, 26(4), pp.301~308. 10.7734/COSEIK.2013.26.4.301

Journal of the Computational Structural Engineering Institute of Korea ISSN:1229-3059(Print) 2287-2302(Online) 한국전산구조공학회 논문집

Preview

Further Improvement of Direct Solution-based FETI Algorithm

ABSTRACT

MAIN

(1)

(2)

(3)

(4)

Fig. 1.

Schematic description of the interface problem for the present FETI-local

(5)

(6)

(7)

(8)

Fig. 2.

Examples of Boolean matrix multiplication

Fig. 3.

Examples of the distributed storage format

Fig. 4.

Framework of the improved FETI-local algorithm

Fig. 5.

Numerical example: carbon fiber/epoxy specimen (longitudinal stress distribution, All units in Pascal)

Fig. 6.

Carbon fiber/epoxy specimen: finite element mesh partitioning by 64 (a) and 128 (b) subdomains

Table 1.

Carbon fiber/epoxy specimen: total computational time and reduction

Fig. 7.

Carbon fiber/epoxy specimen: computational time for each step (The values in the parenthesis indicates the number of the subdomains)

Fig. 8.

Carbon fiber/epoxy specimen: comparison of the scalability (The values in the parenthesis indicates the number of the subdomains)

Table 2.

Carbon fiber/epoxy specimen: average speed-up and its increment in terms of the number of the processors

Fig. 9.

Computer-aided drawing of the plain woven fabric specimen (a) and its fiber bundle (b)

Table 3.

Material properties for the tensile specimen of the plain woven fabric

Fig. 10.

Plain woven fabric specimen: finite element mesh partitioning by 512 subdomains

Fig. 11.

Plain woven fabric specimen: comparison of the computational time (a) and scalability (b)

(9)

Table 4.

Plain woven fabric specimen: uniaxial elastic modulus

Acknowledgements

References