원자로 수치해석 과제 #4

: 2-G, 2-D Eigenvalue problem with Iterative Solution

2008-11833. 반영석

Part 1. Linear System Setup

1. input processing

cross section input file과 main input file의 형식은 지난번과 거의 동일하다 다만 2D problem이므로 composition의 배열 정보가 2-D로서 입력된다. 이외의 경계조건 등의 다른 값들은 동일한 방법으로 받아들인다.

2. linear system setup

이번 문제는 지난숙제와 제일 크게 다른 점이 2D problem이라는 점이다. 따라서 하나의 node에서 좌우의 node의 영향만을 고려하는 것이 아니라 상하방향의 node의 효과 또한 고려되어야 한다. Diffusion equation을 먼저 차분화를 통해 linear system으로 만드는 과정은 다음과 같다.

[eqn: NABLA CDOT J _{g} + SIGMA _{rg} phi _{g} = lambda chi _{g} ( sum _{g prime =1} ^{G} nu SIGMA _{fg prime } phi _{g prime } )+( sum _{eqalign{g prime =1# g prime != g}} ^{G} SIGMA _{g prime g} phi _{g prime } )]     ... (1.1)

위와 같은 neutron balance equation에서 유도를 시작한다. RHS를 간단히 source term으로 [eqn: S_g]라고 두면 식은 다음과 같이 바뀐다.

[eqn: {partial J_x }over{partial x} +{partial J_y }over{partial y} + Sigma_r phi_g =S_g]     ... (1.2)

먼저 하나의 그룹에 대해 생각한다고 보고 아래첨자 [eqn: g]를 모두 생각해서 쓴다. 식 (1.2)를 box scheme을 적용하여 하나의 node에 대해서 양변을 부피적분을 해주면 식은 다음과 같이 바뀔 수 있다.

[eqn: (J_x^r -J_x^l )h_j^x +(J_y^d -J_y^u )h_i^y +Sigma_r phi_m V_m =SV_m]     ... (1.3)

여기서 [eqn: J,``phi,``S]는 하나의 node에서 평균된 값이다. 편의상 [eqn: bar{ `````}]를 생략하였다. 식 (1.3)에서 [eqn: V_m]은 부피로서 여기서는 [eqn: h_j^x times h_i^y]를 뜻한다. [eqn: J]의 각각 첨자들은 [eqn: x ,``y]방향의 current를 의미하고 또 [eqn: r,``l,``u,``d]은 가운데 node 즉, 자식의 [eqn: m] node에서의 방향을 말한다. 각각 right, left, up, down 에 해당한다. 여기서 linear system을 만들기 위해서 current를 flux에 대한 함수로 바꿀 필요가 있다. 따라서 여기에 Fick's law를 적용한다. 각각 방향의 current에 대하여 표면에서 정의되는 flux로부터 경계에서의 current는 다음과 같이 나타낼 수 있다.

[eqn: J_d^p = - {tilde {D}} _m^p (phi_p -phi_m )]     ... (1.4)

식 (1.4)에서 [eqn: p]는 중앙의 [eqn: m] node에서 이웃한 node의 위치를 나타내는 첨자로 이전의 [eqn: r,``l,``u,``d]에 해당하고 [eqn: d]는 [eqn: x,``y]방향을 의미한다. 또한 [eqn: {tildeD }_m^p]는 다음과 같이 정의된다.

[eqn: {tilde D}_m^p =2 { beta_m beta_p } over {beta_m times beta_p}]     ... (1.5)

[eqn: beta_i = D_i over h_i]     ... (1.6)

위의 식들을 종합하면 nodal balance equation (NBE)가 완성된다. 그는 다음과 같다.

[eqn: -{tilde D}_e h_j^y phi_e -{tilde D}_w h_j^y phi_w -{tilde D}_s h_i^x phi_s -{tilde D}_n h_i^x phi_n +({tilde D}_e h_j^y +{tilde D}_w h_j^y +{tilde D}+s h_i^x +{tilde D}_n h_i^x +Sigma_r V_m ) phi_m =SV_m]

... (1.7)

여기에 numbering scheme으로 별다른 기교를 쓰지 않고 natural ordering scheme을 적용하면 [eqn: phi_n , ``phi_s]와는 한 line에 있는 node 수만큼, 그리고 [eqn: phi_e , ``phi_w]와는 1 만큼의 node number의 차이가 있음을 알 수 있다. [eqn: phi_m]의 계수는 diagonal entry가 되고 [eqn: phi_n , ``phi_s]와 [eqn: phi_e , ``phi_w]의 계수는 off-diagonal entry가 된다. 여기서 이제 문제가 생기는 것은 problem domain에서의 boundary에서의 값이다. 예를 들어 [eqn: m] node가 West boundary에 위치해있는 node라면 [eqn: phi_w]와 [eqn: {tilde D }_w]값이 존재하지 않는다. [eqn: phi_w]의 계수는 0으로 되면 되지만 [eqn: phi_m]의 계수 내에 있는 [eqn: {tilde D }_w]값은 새로 정의를 해주어야 한다. 이는 주어진 경계조건 albedo [eqn: alpha]로부터 [eqn: {tilde D }_w]값을 정의할 수 있다.albedo의 정의로부터 시작하여 식 (1.6)의 [eqn: beta]를 다음과 같이 정의함으로서 이 문제를 해결한다.

[eqn: beta_p = alpha_p over 2]     ... (1.8)

위와 같이 경계에서의 flux가 정의되지 않는 점 때문에 migration matrix는 penta-diagonal 하지만 off-diagonal entry가 종종 끊어져있는 모양을 갖게 된다.

식 (1.1~8)을 통해 하나의 group에 대한 linear system을 구축하였다. 주어진 문제는 2-G 문제이기 때문에 위와 같은 linear system이 두 개 존재한다. 그리고 그 두 group의 flux는 scattering source로 coupling 되어있다. 이 coupling 항을 RHS로 넘겨주면 penta-diagnal matrix를 2개 갖고 있는 box-diagonal matrix가 된다.

이 문제에서는 행렬에 다섯줄 만 정보를 저장하기 때문에 sparse 형태로 저장해두는 편이 메모리 절약 면에서 유리하다. 따라서 sparse 형태로 행렬을 짜도록 하였고 또 block-diagonal한 특징을 살리기 위해 cell을 이용하였다.

반복문을 통하여 행렬 M을 구성하였는데 여기서 각 계수들은 직접 하나씩 구하는 방법을 택하였다. 먼저 각 node에서의 [eqn: beta]값을 계산하고 그로부터 각 node에서 사방에 있는 인접 node에 대한 coupling 항 계수 [eqn: tildeD]를 계산하고 그를 행렬의 각 성분에 집어넣는 방법을 사용하였다. 이렇게 코딩하여 얻는 장점은 일단 경계에서 값을 정의하기 쉬워진다. 모든 node에 대해서 같은 방법으로 계산을 하므로 경계에서 따로 추가하여 값을 지정하지 않아도 된다는 장점이 있다. 다만 homogeneous한 매질 내에서도 동일한 계산을 반복한다는 단점이 있다. 일반적인 경우에는 문제를 풀어야하는 매질이 매우 복잡한 구성을 갖고 있으므로 이 방법이 더 유리할 수 있다고 판단하여 무식하지만 일반적인 방법으로 코딩해보았다.

function M=con_M

rmain();rxsec();

global mx my hx hy carr alx_l alx_r aly_l aly_r nx ny nnx nny comp ng

M=cell(nny,nny,ng);

for k=1:ng

for i=1:nny

for j=1:nnx

M{i,j,k}=sparse(nnx,nnx);

end

end

end

betax=zeros(nny+2,nnx+2);  betay=zeros(nny+2,nnx+2);

for g=1:ng

% ------ beta 계산 ---

for ay=1:my

for j=2:ny(ay)+1

line=(ay-1)*ny(ay)+j;

for ax=1:mx

for i=2:nx(ax)+1

k=(ax-1)*nx(ax)+i;

betax(line,k)=comp(carr(ay,ax)).D(g)/hx(ax);

betay(line,k)=comp(carr(ay,ax)).D(g)/hy(ay);

end

end

end

end

for i=2:nny+1

betax(i,1)=alx_l/2;

betax(i,nnx+2)=alx_r/2;

end

for i=2:nnx+1

betay(1,i)=aly_l/2;

betay(nny+2,i)=aly_r/2;

end

% ----- 각 node의 상하좌우의 coefficient 계산 ---

T=zeros(nnx,nny,5);

for ay=1:my

for j=1:ny(ay)

y=(ay-1)*ny(ay)+j;

for ax=1:mx

for i=1:nx(ax)

x=(ax-1)*nx(ax)+i;

T(x,y,1)=2*hx(ax)*(betax(y+1,x+1)*betax(y+1,x))/(betax(y+1,x+1)+betax(y+1,x));%left                  T(x,y,2)=2*hx(ax)*(betax(y+1,x+1)*betax(y+1,x+2))/(betax(y+1,x+1)+betax(y+1,x+2));%right

T(x,y,3)=2*hy(ay)*(betay(y+1,x+1)*betay(y,x+1))/(betay(y+1,x+1)+betay(y,x+1));%up

T(x,y,4)=2*hy(ay)*(betay(y+1,x+1)*betay(y+2,x+1))/(betay(y+1,x+1)+betay(y+2,x+1));%down

T(x,y,5)=comp(carr(ay,ax)).sig_r(g)*hx(ax)*hy(ay);

for i=1:4

T(x,y,5)=T(x,y,5)+T(x,y,i);%diagonal

end

end

end

end

end

% ----- coefficient로부터 행렬 구성

for ay=1:my

for j=1:ny(ay)

line=(ay-1)*ny(ay)+j;

B=zeros(nnx);

for ax=1:mx

for i=1:nx(ax)

k=(ax-1)*nx(ax)+i;

B(k,k)=T(k,line,5);

if k+1<=nny

B(k+1,k)=-T(k+1,line,1); % lower diagonal = left node coupling

B(k,k+1)=-T(k,line,2); % upper diagonal = right node coupling

end

end

B(k,k)=T(k,line,5);

end

M{line,line,g}=B;

end

end

for ay=1:my

for j=1:ny(ay)

line=(ay-1)*ny(ay)+j;

if line>=nny

break

end

L=sparse(nnx);

U=sparse(nnx);

for ax=1:mx

for i=1:nx(ax)

k=(ax-1)*nx(ax)+i;

U(k,k)=-T(k,line,4); % upper matrix = lower node coupling

L(k,k)=-T(k,line+1,3); % lower matrix = upper node coupling

end

end

M{line,line+1,g}=U;

M{line+1,line,g}=L;

end

end

end

part 2. Point and Line SOR, and Spectral Radius

1. Point SOR

각 group의 linear system을 풀어주는 inner iteration을 수행하기위해 point SOR scheme을 사용한다. SOR scheme의 iteration matrix의 구조는 다음과 같이 유도된다. 먼저 다음과 같은 linear system을 푼다고 하자

[eqn: A phi=b]     ... (2.1)

식 (2.1)의 양변에 [eqn: w]를 곱해주고 [eqn: A]를 대각성분 [eqn: D]와 상삼각행렬 [eqn: U]와 하삼각행렬 [eqn: L]로 나눈다. 여기서 앞서 이 행렬 [eqn: A]의 대각성분을 제외하고 나머지 성분들은 음수임을 알았으므로 그를 활용하여 [eqn: A=D-L-U] 로 [eqn: L]과 [eqn: U]를 양수로 정의해준다. 그리고 적절히 나누어 [eqn: wA]를 [eqn: M-N]꼴로 바꾸면 다음과 같다.

[eqn: wA=w(D-L-U)+D-D# `````````````=[D-wL]-[wU+(1-w)D]# `````````````=M-N]     ... (2.2)

이를 식 (2.1)에 대입하여주면 다음과 같다.

[eqn: M phi - N phi = wb]     ... (2.3)

여기서 [eqn: N phi]항을 우변으로 이항해주고 좌변에는 새로운 solution을 우변에는 이전 iteration의 solution을 넣어주면 다음과 같다.

[eqn: [D-wL] phi^(l) = wb+ wU phi^(l-1) +(1-w)D phi^(l-1)]    ... (2.4)

그러므로 iteration matrix는 다음과 같이 쓸 수 있다.

[eqn: T_SOR = [D-wL]^-1 [(1-w)D +wU]]     ... (2.5)

여기서 우리가 구하는 것은 [eqn: phi^(l)]인데 이중 [eqn: L]과 곱해지는 항은 새로 구한 solution으로부터 구하는 것이 유리하므로 다음과 같이 된다.

[eqn: D phi^(l) = w[b+L phi^(l) +U phi^(l-1) ]+(1-w)D phi^(l-1)]     ... (2.6)

여기서 식 (2.5)의 양변에 [eqn: D^-1]를 곱해주고 SOR scheme 이전의 Gauss-Seidal scheme이 다음과 같다는 점을 이용하여 식 (2.8)을 얻는다.

[eqn: phi_GS^(l) = D^-1 [ b+L phi^(l) +U phi^(l-1) ]]     ... (2.7)

[eqn: phi_SOR^(l) = w phi_GS^(l) + (1-w) phi^(l-1)]     ... (2.8)

SOR scheme은 Gauss-Seidal scheme의 [eqn: w]만큼의 extrapolated된 solution을 얻는 방법이라는 것을 알 수 있다. 따라서 GS scheme을 구현한 후 그것과 이전 step에서의 solution의 외삽을 해주면 SOR scheme의 solution이 된다.

이를 코드로서 구현한다. 여기서는 아직 outer iteration이 완성되기 전이므로 이 SOR scheme을 시험하기위해 RHS에 uniform source [eqn: 1 n/sec/cm^3]를 두고 [eqn: w=1 ``,``1.25]인 경우에 residual error가 [eqn: 1E-3]이하로 수렴할 때까지 필요한 iteration 수를 결정한다. 그 결과는 다음과 같다.

g=1

w=1.00

inner iter = 197 , e_ref = 1.266873e-002 , e_sol = 2.745700e-004, e_res = 2.745700e-004

w=1.25

inner iter = 124 , e_ref = 1.046526e-002 , e_sol = 3.805773e-004, e_res = 3.044618e-004

g=2

w=1.00

inner iter = 35 , e_ref = 3.225550e-003 , e_sol = 5.648954e-004, e_res = 5.648954e-004

w=1.25

inner iter = 23 , e_ref = 1.780947e-003 , e_sol = 5.635964e-004, e_res = 4.508771e-004

2. Line SOR

Line SOR은 기본적인 배경은 Point SOR과 동일한데 다만 행렬을 식 (2.2)와같이 나눌 때 [eqn: D,``L,``U]를 element가 아닌 block matrix로 보는 것이다. 이렇게 되면 식 (2.6)에서 [eqn: D]가 tri-diagonal matrix가 되므로 [eqn: D^-1]가 바로 구해지지 않는다는 문제가 생기는데 이것은 LU-factorization을 이용하여 풀어주면 된다. 이 Line SOR의 기본 장점은 Jacobi method와 Gauss-Seidal method의 차이와 비슷하게 자신의 node를 계산할 때 자신의 line에 있는 node의 결과를 포함하여 같이 구한다는 점이다. 이렇게 함으로서 더욱 오차가 적은 계산결과를 얻을 수 있다.

2.1에서와 같은 방법으로 routine을 test해본 결과는 다음과 같다.

g=1

w=1.00

inner iter = 106 , e_ref = 9.519565e-003 , e_res = 4.164932e-004 , e_reduction = 0.958084

w=1.25

inner iter = 66 , e_ref = 8.117467e-003 , e_res = 4.819090e-004 , e_reduction = 0.931145

g=2

w=1.00

inner iter = 21 , e_ref = 1.652220e-003 , e_res = 5.827492e-004 , e_reduction = 0.739726

w=1.25

inner iter = 13 , e_ref = 1.185222e-003 , e_res = 6.541693e-004 , e_reduction = 0.596143

2.1과 2.2의 결과를 종합하여보면 다음과 같다.

Group	Point SOR		Line SOR
	w= 1.00	w= 1.25	w= 1.00	w= 1.25
1	197	124	106	66
2	35	23	21	13

group 1보다 group 2가 더 빠르게 수렴하였고 또 [eqn: w]가 1.00인 경우보다 1.25인 경우에 더 나은 수렴속도를 보여주었다. 참고로 [eqn: w=1]인 경우는 Gauss-Seidal method와 동일하므로 이것으로 SOR method가 Gauss-Seidal method보다 빠름을 볼 수 있다. 다음으로 Point SOR과 Line SOR을 비교하여보면 iteration 수가 두 배정도 차이난다. iteration수 면에서는 Line SOR이 우월함을 확인할 수 있었다. 효율성 면에서는 Line SOR 내부에서 또 다른 linear system을 푸는 routine이 들어있으므로 그것까지 고려한 FLOP을 통해 비교해보아야 한다.

3. Over relaxation factor [eqn: w]

앞에서는 [eqn: w]값을 1.00, 1.25인 경우로 고정시키고 test를 실행하였다. 이제는 더 빠른 수렴을 위해서는 정확한 [eqn: w]값이 필요하고 이 값은 [eqn: w=1]인 경우 즉, Gauss-Seidal method의 수렴속도인 spectral radius [eqn: rho]값으로부터 Young's Formula를 통하여 수렴속도가 최적이 되는 [eqn: w]값을 구할 수 있다. 이 Young's Formula는 다음과 같이 유도된다.

먼저 SOR iteration matrix와 Jacobi iteration matrix를 다음과 같이 둔다.

[eqn: T_SOR = T=(D-wL)^-1 [wU+(1-w)D]]     ... (2.9)

[eqn: T_Jacobi = T_J = D^-1 (L+U)]     ... (2.10)

이 두 행렬 [eqn: T]와 [eqn: T_J]의 eigenvalue를 각각 [eqn: gamma]와 [eqn: lambda], eigenvector를 각각 [eqn: phi]와 [eqn: psi]라고 하자. 여기서  [eqn: (```````````)^-1]항을 넘겨서 SOR iteration matrix에 대하여 다음과 같이 쓸 수 있다.

[eqn: (1-w)D phi_i +w U phi_i+1 = gamma (D phi _i -wL phi_i-1 )]     ... (2.11)

식 (2.11)에서 [eqn: D]에 대해서 항을 묶어주면 식 (2.12)가 된다. 여기서 [eqn: phi _i = alpha^i psi _i]로 두 eigenvector 사이에 [eqn: alpha]배의 관계가 있다고 가정을 하면 식 (2.13)이 얻어진다. 여기서 양변의 [eqn: alpha^i]를 약분하였다.

[eqn: wL phi_i-1 + wU phi_i+1 = (gamma+w-1)D phi_i]     ... (2.12)

[eqn: w [ gamma over alpha L psi_i-1 + U alpha psi_i+1 ] = (gamma+w -1 )D psi_i]     ... (2.13)

여기에 Jacobi iteration matrix로부터 얻어진 다음을 대입한다.

[eqn: 1 over lambda {[ L psi_i-1 + U psi _i+1 ]} = D psi_i]     ... (2.14)

식 (2.14)를 식 (2.13)우변에 대입해주고 [eqn: psi_i-1]과 [eqn: psi _i+1]의 계수를 비교하면 [eqn: w]와 [eqn: gamma], [eqn: alpha], [eqn: lambda]의 관계가 얻어진다.

[eqn: w gamma over alpha = (gamma+w -1) 1 over lambda]     ... (2.15-a)

[eqn: w alpha = (gamma +w -1 ) 1 over lambda]     ... (2.15-b)

두 식으로부터 [eqn: alpha = root gamma]임을 알 수 있다. 여기서 [eqn: w=1]인 경우 [eqn: gamma =lambda^2]이 되므로 Gauss-Seidal method가 Jacobi method보다 2배 빠름을 확인하고 넘어간다. 일단 [eqn: alpha]식을 통해 [eqn: alpha]를 소거하고 식 (2.15)로부터 [eqn: gamma]가 최소가 되는 조건을 찾으면 그 때의 [eqn: w]가 optimized 된 over-relaxation factor가 된다. [eqn: gamma(w)]식은 다음과 같다.

[eqn: gamma^2 - 2 [ (1-w) + {w^2 lambda^2 } over 2 ] gamma + (w-1)^2 =0]     ... (2.16)

여기서 근의 공식을 적용하면 두 근은 다음과 같이 표현된다.

[eqn: gamma_+ = (1-w) + {w^2 lambda^2 } over 2 +w lambda root delta]     ... (2.17-a)

[eqn: gamma_- = (1-w) + {w^2 lambda^2 } over 2 -w lambda root delta]     ... (2.17-b)

단 여기서 [eqn: delta= (1-w) + {w^2 lambda^2 } over 4]이다. 식 (2.17)로부터 이제 [eqn: |gamma|]가 최소가 되는 [eqn: w]값을 구하여야한다. 이 [eqn: delta]가 음수일 때 [eqn: gamma]는 복소수 값이 되는데 그 크기 [eqn: |gamma|]는 [eqn: w-1]이 된다. [eqn: delta]가 zero인 경우는 [eqn: w]가 2가지가 존재하는데 [eqn: w]가 증가할수록 [eqn: |gamma|]가 증가하므로 작은 해를 택하면 다음과 같다.

[eqn: w_opt = w_- = 2 over {1 + root {1-lambda^2 } } = 2 over {1+{root {1- rho_GS}]     ... (2.18)

이것이 Young's Formula 이다.

이를 통하여 다시 위에서의 test를 실행하면 iteration 수는 다음과 같다.

Group	Point SOR		Line SOR
	w	iteration	w	iteration
1	1.745421	37	1.661077	26
2	1.442765	16	1.324082	11

앞으로 inner iteration은 이 iteration 수로 고정하고 문제를 풀게 된다.

4. FLOP

코드를 어떻게 작성하느냐에 따라 계산량에 차이가 있을 수 있으므로 먼저 코드를 소개한다. 다음은 Point SOR의 코드이다.

function [new,iter]=solve_PSOR(A,B,old,w,it)

global epsin ninmax

A=sparse(A);

m=length(A);

k=1;  iter=0;  e_res=1;  temp=1;  flag=1;  e_1=1;

new=zeros(m,1);newr=zeros(m,1);

%--- for sparse matrix calculation ---

[yy,xx]=find(A');

l=length(xx);

xx=[ 1;xx ; 1];  yy=[ 1;yy ; 1];

x=xx;  y=yy;

앞부분에서는 별다른 연산을 하지 않는다. sparse matrix를 연산하기 위해서 원소의 위치를 찾아 저장하는 작업을 수행한다. 다음으로는 본격적인 point SOR 연산을 한다.

while iter<it %e_res*e_1>epsin || iter < 5 && iter < ninmax

if iter==1

e_1=1/e_r(1);

end

i=2;

while i<length(x)

sum_1=0;

sum_2=0;

while x(i)==x(i-1) || flag==0

if flag==0

flag=1;

end

j=i;

if y(j)<x(i)

sum_1=A(x(i),y(j))*new(y(j))+sum_1;

elseif y(j)>x(i)

sum_2=A(x(i),y(j))*old(y(j))+sum_2;

end

if x(i)~=x(i+1)

flag=0;

i=i+1;

break

end

i=i+1;

end

i=i-1;

newr(x(i))=(B(x(i))-sum_1-sum_2)/A(x(i),x(i));

new(x(i))=(w*newr(x(i))+(1-w)*old(x(i)));     % SOR은 이부분만 다르다.

i=i+1;

end

x=xx;

y=yy;

이 코드는 일반적인 SOR코드와 모양이 살짝 다르다. 이유는 sparse matrix의 연산을 위해 반복문의 구조를 바꾸었기 때문이다. bold체로 표시된 부분이 주된 연산을 하는 부분이다. 먼저 위의 조건문에서 off-diagonal의 수만큼 곱셈과 덧셈 연산을 하는데 그 수는 각각 10404번씩이다. (덧셈 + 10404, 곱셈 +10404)

다음으로 51*51(=2601)개의 해 벡터를 계산하기 위해서 두 번의 연산을 더 하는데 한번은 GS solution을 구하기 위해 2601*2번의 뺄셈연산과 2601번의 나눗셈 연산을 한다. ( 뺄셈 +5202, 나눗셈 +2601) 그에 이어 그 해를 통해 extrapolate하여 SOR solution을 구하기 위해 한 번의 덧셈과 두 번의 곱셈연산을 한다. (덧셈 + 2601, 곱셈 + 5202)

다음으로는 연산을 마치고 error monitoring 하는 과정을 거치는데 코드는 다음과 같다.

% ---------- #.6 estimate psuedo error ----------

iter=iter+1;

dif=new-old; % differnce btw psi_o , n

e_sol=(abs(dot(dif,dif)/dot(new,new)))^0.5;

e_s(iter)=e_sol;

dif=newr-old; % differnce btw psi_o , n

e_res=(abs(dot(dif,dif)/dot(new,new)))^0.5;    %residual error

e_r(iter)=e_res;

rho=e_res/temp;

old=new;    %다음단계를 위한 업데이트

temp=e_res;

end

여기서는 iteration check를 위해 한 번의 덧셈을 하고 (덧셈 + 1) psuedo error, residual error를 구하기 위해 뺄셈 두 번과 내적 네 번 그리고 루트 연산 두 번이 들어간다. (뺄셈 +2 덧셈 +10404, 곱셈 + 10404, 제곱근 +2) 그리고 error reduction을 보기 위해서 한 번의 나눗셈 연산을 추가한다. (나눗셈 +1 )

이를 종합하여보면 한번의 iteration에 필요한 총 연산수는 다음과 같다.

연산 종류	덧셈	뺄셈	곱셈	나눗셈	기타연산
연산수	23410	5204	26010	2602	제곱근 2
	28614

다음으로 Line SOR의 코드를 본다. Line SOR이 point SOR과 연산수에서 가장 차이가 날 부분은 linear system을 한번 풀어줘야 한다는 점이다. 따라서 LU-factorization 하는 연산이 한번 추가 되어야 한다. 이는 첫 번째 iteration에서만 수행하고 그 후에는 저장을 해두고 계속 불러와서 쓸 수 있기 때문에 iteration 수가 증가할수록 유리해진다.

function [new,iiter]=solve_LSOR(M,B,old,w,it)

global iter ninmax epsin g

iter=0;  B=s2c(B);  old=s2c(old);

n=length(M);  L=B;  U=B;

iiter=0;

new=old;  newr=old;  e_res=1;  temp=1;

if iter==0

if g==1

for i=1:n

[L{i},U{i}]=luf(M{i,i});

end

save inlu1 L U

elseif g==2

for i=1:n

[L{i},U{i}]=luf(M{i,i});

end

save inlu2 L U

end

else

if g==1

load inlu1

elseif g==2

load inlu2

end

end

newm=c2m(new);  oldm=c2m(old);  newrm=c2m(newr);

e_1=1;

다음은 위에서 사용한 LU-factorization 함수이다. 사실 LU factorization은 전체 iteration 수에 상관없이 처음에 딱 한번만 시행되고 그 이후로는 저장하고 계속 불러와서 쓰므로 FLOP에는 포함하더라도 그 수는 iterative 계산을 하는 것에 비해서 무시할만하다.

그 다음으로는 Line SOR scheme 부분이다.

while iiter<it %e_res*e_1>epsin || iiter < 5 && iiter < ninmax

if iiter==1

e_1=1/e_r(1);

end

RHS=(B{1}-M{1,2}*old{2});

newr{1}=solve_axb(L{1},U{1},RHS);

new{1}=w*newr{1}+(1-w)*old{1};

for i=2:n-1

RHS=(B{i}-M{i,i-1}*new{i-1}-M{i,i+1}*old{i+1});

newr{i}=solve_axb(L{i},U{i},RHS);

new{i}=w*newr{i}+(1-w)*old{i};

end

RHS=(B{n}-M{n,n-1}*new{n-1});

newr{n}=solve_axb(L{n},U{n},RHS);

new{n}=w*newr{n}+(1-w)*old{n};

% ---------- #.6 estimate psuedo error ----------

iiter=iiter+1;

dif=newm-oldm; % differnce btw psi_o , n

e_sol=(abs(dot(dif,dif)/dot(newm,newm)))^0.5;

e_s(iiter)=e_sol;

dif=newrm-oldm; % differnce btw psi_o , n

e_res=(abs(dot(dif,dif)/dot(newrm,newrm)))^0.5;    %residual error

e_r(iiter)=e_res;

rho=e_res/temp;

old=new;    %다음단계를 위한 업데이트

temp=e_res;

end

위에서는 RHS를 계산하는데 51개의 line에 대해서 연산을 수행해주는데 각 계산은 51개의 성분에 2번씩의 덧셈과 곱셈을 하고 (곱셈 + 102, 덧셈 + 102) Forward and backward substitution을 통해 [eqn: Ax=b]를 풀어준다. 이 계산 수는 아래의 코드로부터 51개의 성분에 대하여 forward substituiton에서는 51개의 성분에 곱셈, 덧셈, 뺄셈 한번 씩 그리고 backward substitution에서 51개의 성분에 곱셈, 나눗셈, 덧셈, 뺄셈 한번 씩 수행된다. ( 곱셈 + 102, 나눗셈 +51 , 덧셈 +102, 뺄셈 +102)

% solve AX=B by subsitution

% for tri-diagonal matrix solution

function [X]=solve_axb(L,U,B)

n=length(B);

Y=zeros(n,1);

X=zeros(n,1);

Y(1)=B(1);

for i=2:n       %Y를 forward substitution으로 구함

s=0;

for j=i-1:i-1

s=s+L(i,j)*Y(j);

end

Y(i)=B(i)-s;

end

X(n)=Y(n)/U(n,n);

for i=n-1:-1:1  %X를 backward substitution으로 구함

s=0;

for j=i+1:i+1

s=s+U(i,j)*X(j);

end

X(i)=(Y(i)-s)/U(i,i);

end

종합하면 하나의 Line에 대해서 곱셈 204번, 덧셈 204번, 나눗셈 51번, 뺄셈 102번이 실행되므로 여기에 51씩을 곱하여 준다. 그리고 error monitoring은 point SOR과 같은 방법을 사용하므로 같은 값을 사용해주어 종합하면 iteration 수는 다음과 같다.

연산 종류	덧셈	뺄셈	곱셈	나눗셈	기타연산
연산수	20655	5204	20808	2602	제곱근 2
	25859

Point SOR과 Line SOR의 FLOP 차이는 크게 나지 않았다. 그러나 수렴속도는 Line SOR이 더 빨라 group 1에 대해서 PSOR, LSOR 각각 37, 26번 수행되므로 총 수행되는 곱셈 수는  약 962,000회, 520,000회로 LSOR이 약 1.8배 빠름을 알 수 있다. 실제 계산에 소요되는 CPU time을 비교해보았을 때 각각 5.8초, 2.0초 소요되어 계산 시간에서는 약 3배 유리함을 관찰할 수 있었다.

실제로 코드에서 FLOP을 계산한 것이 아니고 손으로 계산한 결과와 CPU time을 비교하여 어느 정도 차이가 있다고 생각된다. 손으로 계산하였기에 분명 어딘가에 또 빠뜨린 항이 있어서 차이가 생겼을 것이라 생각된다.

결론적으로 같은 수렴조건을 달성하는데 필요한 계산 수는 LSOR이 2~3배 좋음을 확인할 수 있었다.

part 3. Inner and Outer Iterations

1. ChebyShev Acceleration

주어진 2D-2G 문제를 푸는 최종목적은 [eqn: k_eff]값을 구하는 것이다. 이것은 migration matrix의 eigenvalue에 해당하는 값인데 이를 구하기 위해서 앞의 숙제에서 사용한 ChebyShev acceleration을 이용한다. 이 ChebyShev acceleration을 하기위해서는 iterative method내에서 [eqn: Ax=b]꼴의 linear problem을 풀어야하는 상황이 발생한다. 이를 풀어주는 scheme으로서 part 2.에서 작성한 SOR method를 사용한다. 앞서 Line SOR method가 여러 측면에서 Point SOR method보다 우월함이 증명되었으므로 Line SOR 만을 이용한다. 따라서 종합하여보면 주어진 상황을 차분화하여 이 시스템의 eigenvalue를 계산하는데 그를 위해 ChebyShev acc.를 사용하고 그를 사용하기위해 또 그 iteration 내부에 SOR method를 포함하여 사용하는 상황이다.

ChebyShev acceleration에 대해서는 이전 숙제에서 다루었으므로 그대로 사용한다. 다만 그때는 1-dimensional 문제였기에 numbering을 하나의 line에 대해서만 했으므로 그 수를 2-dimensional problem에 맞게 확장만 시켜주면 된다. Line SOR-ChebyShev Acc code를 소개하면 다음과 같다.

%function solve_2D_2G_cheby_LSOR

arr=sprintf(' ============ solve_2D_2G_cheby_LSOR ============ ');  disp(arr);

clear all

load reference  % 10^-12 수렴시킨 해를 reference로 사용한다.

global comp ng

global carr mx my nx ny nnx nny maxorder startnum conv_eig conv_sol hx hy

global rpsi ninmax noutmax epsin g

rxsec();  rmain();

% nn=n*nc;

ntot=nnx*nny; % = (mx*nx)*(my*ny)

e_eig=1;  t_eig=1;  e_sol=1;  t_sol=1;  k_o=1;  lam_o=1;  xi=1;  flag=0;  sig=0.9; gamma=1;

rpsi=rpsi/(dot(rpsi,rpsi))^0.5;

phi_o=cell(ng,1);

phi_n=cell(ng,1);

% --- M을 구성하고 w와 iteration 수를 결정한다. ---

for g=1:ng

phi_o{g,1}=ones(ntot,1);

phi_n{g,1}=ones(ntot,1);

M{g}=M3N(g);

w(g)=Lrho(M{g});

it(g)=test_LSOR(M{g},w(g));

end

arr=sprintf(' w1 = %f , w2 = %f',w(1),w(2));   disp(arr);

arr=sprintf(' iter_1 = %d , iter_2 = %d',it(1),it(2));   disp(arr);

iter=0; order=1-startnum;

% ---------- #.1 determine fission source (for 2D) ----------

for cy=1:my

for j=1:ny(cy)

y=(cy-1)*ny(cy)+j;

for cx=1:mx

for i=1:nx(cx)

x=(cx-1)*nx(cx)+i+(y-1)*nnx;

psi_o(x)=0;

for g=1:ng

psi_o(x)=psi_o(x)+comp(carr(cy,cx)).sig_f(g)*phi_o{g,1}(x);

end

end

end

end

end

while e_eig>conv_eig || e_sol>conv_sol && iter < noutmax

% ---------- #.2 determine source at each node

for g=1:ng

for cy=1:my

for j=1:ny(cy)

y=(cy-1)*ny(cy)+j;

for cx=1:mx

for i=1:nx(cx)

x=(cx-1)*nx(cx)+i+(y-1)*nnx;

sum=0;

for gg=1:ng

if(gg==g) continue end

% scattering source

% scat(a,b); scat b->a source

sum=sum+comp(carr(cy,cx)).scat(g,gg)*phi_n{gg,1}(x);

end

%total source term

s(x)=(lam_o*comp(carr(cy,cx)).chi(g)*psi_o(x)+sum)*hx(cx)*hy(cy);

end

end

end

end

% ---------- #.3 solve for flux for group

% solve M{g,g} phi{g} = s

phi_n{g,1}=solve_LSOR(M{g},s,phi_n{g,1},w(g),it(g));  % 지난번 숙제와 다른 inner iteration

end

% ---------- #.4 determine new fission source ----------

for cy=1:my

for j=1:ny(cy)

y=(cy-1)*ny(cy)+j;

for cx=1:mx

for i=1:nx

x=(cx-1)*nx(cx)+i+(y-1)*nnx;

psi_n(x)=0;

for g=1:ng

psi_n(x)=psi_n(x)+comp(carr(cy,cx)).sig_f(g)*phi_n{g,1}(x);

end

end

end

end

end

psi_nn=psi_n;

%---------- ChebyShev Acceleration ----------

if order == maxorder || xi>1 && order >3

arr=sprintf('reset order = %d',order-1);

order=1;

disp(arr);

xi=1;

sig=sigg;

flag=1;

end

if order>0

if order ==1

if flag==0;

sig=rho;

end

a=2/(2-sig);

b=0;

psi=psi_n-psi_o; % differnce btw psi_o , n

e_sol=(abs(dot(psi,psi)/dot(psi_n,psi_n)))^0.5;

psi_n=a*psi_n+(1-a)*psi_o;

ceo=e_sol;

else

gamma=2/sig-1;

a=4/sig*(cheby(order-1,gamma)/cheby(order,gamma));

b=(1-sig/2)*a-1;

psi=psi_n-psi_o; % differnce btw psi_o , n

e_sol=(abs(dot(psi,psi)/dot(psi_n,psi_n)))^0.5;

psi_n=a*psi_n+(1-a+b)*psi_o-b*psi_oo;

zeta=e_sol/ceo;

xi=zeta*cheby(order-1,gamma);

sigg=sig*(1+cheby(1/(order-1),xi))/2;

end

end

% ---------- #.5 estimate eigenvalue ----------

k_n=k_o*(dot(psi_n,psi_n)/dot(psi_o,psi_n));

lam_n=1/k_n;

% ---------- #.6 estimate psuedo error ----------

e_eig=abs(1-k_o/k_n);

psi=psi_n-psi_o; % differnce btw psi_o , n

e_sol=(abs(dot(psi,psi)/dot(psi_n,psi_n)))^0.5;

p_n=psi_n/(dot(psi_n,psi_n))^0.5;

psi=p_n-rpsi;

e_ref=(abs(dot(psi,psi)/dot(rpsi,rpsi)))^0.5;

psi=psi_nn-psi_o; % differnce btw psi_o , n

e_res=(abs(dot(psi,psi)/dot(psi_n,psi_n)))^0.5;

% ---------- #.7 dominance ratio rho ----------

rho=e_sol/t_sol;

% ---------- #.8 save & update data ----------

t_sol=e_sol;    psi_oo=psi_o;    psi_o=psi_n;    phi_o=phi_n;

lam_o=lam_n;    k_o=k_n;    iter=iter+1;

e_e(iter)=e_eig;    e_s(iter)=e_sol;    e_f(iter)=e_ref;

e_rr(iter)=e_res;    r(iter)=rho;

order=order+1;

% ---------- #.9 monitoring error ----------

arr=sprintf('outer iter = %d ,k-eff = %8.6f, rho = %e , error_eig = %e , error_sol = %e , error_ref = %e , xi = %f , r = %f',iter,k_o,sig,e_eig,e_sol,e_ref,xi,gamma);

disp(arr);

end

draw_s(phi_o);  draw_ss(psi_o);  legend('heat graph');

y=1:iter;

fig7=figure(7);

fig7=semilogy(y,e_e,y,e_s,y,e_rr,y,e_f);

legend('eigenvalue error','fission source error','residual error','reference error');

grid on

arr=sprintf('===========end==========');

disp(arr);

2. [eqn: k_eff] and Flux Distribution

part 3.1에서의 코드를 사용하여 해당 문제의 [eqn: k_eff]값과 group별 flux 그리고 원자로 열출력에 비례하는 fission source distribution을 확인한다.

reference 값으로는 ChebyShev Acc.을 하지 않고 Line SOR을 통하여 fission source의 psuedo error가 [eqn: 10^-12]이하로 수렴할 때까지 반복하여 구한 값을 사용하였다. 그 결과로 출력된 값은 다음과 같다.

[eqn: k_eff]=9.378466e-001

다음은 Line SOR을 통한 해이다.

outer iter = 15 ,  e_eig = 5.45e-007 , e_psi = 4.78e-006 , e_ref = 2.75e-006

k-eff = 0.937847

유효숫자 여섯 자리까지 잘 일치함을 확인 할 수 있었다. 다음으로 flux와 fission source 분포를 본다.

< Group 1 Flux >

< Group 1 contour >

< Group 2 Flux >

< Group 2 contour >

< Power distribution >

< Power distribution contour>

first group 의 flux는 cos shape을 갖는 것을 확인하였다. 그러나 두 번째 group flux는 9번 assembly에서 낮은 flux를 보였다. 그리고 특히 group 1 flux가 3번 assembly(moderator)에서 급격하게 줄어들었던 것에 비해 group 2 flux는 3번 지역에서 상대적으로 많이 관찰되었다.

그리고 power distribution을 보면 출력이 가장 높은 부분은 8,9번 assembly 사이였고 평균적인 분포를 보았을 때는 내부의 8번 assembly에서의 출력이 가장 높았다. 상대적으로 leakage가 가장 심한 바깥쪽의 9번 assembly에서 출력이 가장 낮았고 그 외의 밖의 공간에서는 fission이 일어나지 않는 moderator 구간이기 때문에 아무런 출력도 관찰되지 않았다.

다음으로 iteration 수에 따라서 변화하는 flux의 수렴 모양을 관찰하였다.

위는 순서대로 [ 1 2 3 , 4 5 6 ] 번째 step에서의 group 1의 flux shape이다. 중앙에서 flux가 상대적으로 낮았다가 점차 올라가는 모양을 보였다. 그 외의 구간에서는 거의 모양이 일정하였다.

위는 group 2의 flux shape이다. 초기에 3번으로 이루어진 moderator 공간에서 주로 flux를 갖다가 후에 중앙에서 점점 높아지는 경향을 보였다.

마지막으로 power distribution을 나타낸 그림이다. 9번 영역에서 주로 열출력이 발생하는 모양은 그대로 유지되는데 이 역시 모양이 점차 중앙이 높은 모양으로 자리 잡아감을 확인 할 수 있었다.

3. Inner iteration number

inner iteration의 수렴 조건을 변경시켜가면서 최적의 iteration number를 구한다. inner iteration 수가 많으면 분명 오차가 감소하여 outer iteration 수를 줄이는 효과가 있긴 하지만 그만큼 계산이 더 많이 필요로 된다는 단점이 있다. 따라서 최적이 되는 경우를 찾는 것이 바람직하다. 앞서서는 inner iteration의 수렴조건을 residual error가 [eqn: 1E-3]이하로 수렴할 때 iteration을 종료하도록 하였고 그 때 iteration 수가 각각 26, 11회였다. 다음은 residual error convergence criteria에 대한 각 그룹별 iteration 소요 수이다. 이 test는 optimized된 over-relaxation factor에 대해서 실행되었다.

수렴 조건		[eqn: 1E-0]	[eqn: 1E-1]	[eqn: 1E-2]	[eqn: 1E-3]	[eqn: 1E-4]	[eqn: 1E-5]
inner	g1	5	6	15	26	36	44
	g2	5	5	8	11	15	19
outer		53	37	18	15	16	16
total		530	407	414	555	816	1008

outer iteration의 수는 갈수록 어느 값에 수렴함을 관찰하였다. 반대로 inner iteration의 수렴조건을 tight하게 잡을수록 그에 비례하여 inner iteration 수는 증가하였다. 이를 고려해보면 optimized 된 값은 [eqn: 1E-1]~2 라고 결론지을 수 있었다. 둘의 총 계산 수는 크게 차이 나지 않지만 [eqn: 1E-2]의 경우가 더욱 안정적이여 보였기에 이를 사용하였다.

다음으로 input file을 응용하여 두 개의 assembly 사이에 1/10 크기의 3번 assembly를 넣어보았다. 아래의 기존 결과와 비교하여 보면 감속에 의해 8, 9번 assembly 사이에서 group 2 flux가 대폭 증가함을 확인할 수 있었고 그 때문에 9번 영역에서 fission이 대폭 줄어들은 것을 관찰할 수 있었다. (아래 그림은 순서대로 group 1, 2 , fission source shape이다.)

또 추가적으로 가장 안쪽의 assembly를 바꾸어 가며 k_eff 값을 비교해 보았다.

assembly #	[eqn: k _{eff}]
1	0.833147
2	0.833154
3	0.801087
4	0.938281
5	0.935951
6	0.936227
7	0.936416
8	0.937847
9	0.948839
10	0.760753
11	0.835368

[eqn: k _{eff}]값으로 미루어 보아 1,2,3,11 번 assembly는 비슷한 특성을 갖는 물질로 추정되고 앞서 3번 assembly가 moderator 라고 추정하였으므로 위 네 개의 assembly는 moderator라고 생각된다. 4,5,6,7,8,9 번은 모두 fuel 이라고 생각되었다. 기존의 8번 assembly를 대체하여도 [eqn: k _{eff}]값에 큰 변화를 주지 못하였기 때문에 fuel 이라고 생각하였다. 마지막으로 10번 assembly는 [eqn: k _{eff}]값을 크게 감소시켰는데 그 이유로 이것은 control rod 따위의 것이라고 생각한다.

Part 4. Node Size Effect

1. Assembly-wise flux

node의 크기를 키우거나 줄이면서 그에 의한 변화를 관찰하려한다. 상대적인 변화를 관찰하기 위해서는 assembly별로 부피를 가중치로 두어 평균 낸 flux가 필요하다. 그러나 지금 경우에서는 mesh size를 동일하게 차분화 하였으므로 그냥 평균 내어도 무방하다. 앞서 Part 3.에서 계산한 결과를 합하고 node수로 나누어 평균을 내주는 routine을 새로 작성한다. 그 코드는 다음과 같다.

function [flux]=assembly(phi)

global mx my nx ny nnx hx hy

rmain();

sx=1;   sy=1;   k=0;

for i=1:my

for j=1:mx

flux(i,j)=0;

for y=sy:sy+ny(i)-1

for x=sx:sx+nx(j)-1

index=(y-1)*nnx+x;

k=k+1;

flux(i,j)=flux(i,j)+phi((y-1)*nnx+x)*hx(j)*hy(i);

end

end

flux(i,j)=flux(i,j)/(hx(j)*hy(i)*nx(j)*ny(i));

sx=x+1;

end

sy=y+1;

sx=1;

end

2. Node Size Effect

당연히 mesh size를 작게 하였을 때 error가 적을 것이다. 따라서 가장 작게 자른 경우를 기준으로 하여 그것의 절반, 또는 사분의 일로 줄였을 때의 eigenvalue와 fission source distribution을 관찰한다. 먼저 아래는 reference value로 사용한 값이다. ( residual error를 [eqn: 1E-12]이하로 수렴시킨 해 )

[eqn: k_eff]=9.378466e-001

이 값과 mesh수를 2에서 2배씩 증가시키면서 [eqn: k_eff]값의 상대오차의 거동을 살펴보았다.

mesh 수	k	relative error
2	0.935737	2.2E-03
4	0.933447	4.7E-03
8	0.935733	2.3E-03
12	0.93708	8.2E-04
16	0.937739	1.1E-04

확실히 mesh수가 증가함에 따라서 상대오차가 줄어든다는 느낌은 들었다. mesh를 2개로 나누었을 때에는 오히려 작은 오차를 보였지만 그것은 그저 우연이라고 생각된다.

k 값의 error 보다는 fission source의 error 보는 것이 더욱 정확하기 때문에 fission source error의 node size effect를 관찰해보았다. fission source error는 fission source를 normalize 한 후 그것을 assembly-wise average를 취하여 계산한 값을 가지고 계산하였다.

mesh 수	fission source error	relative error
2	0.0668378
4	0.0166191	4.02174922
8	0.0037728	4.4049893
16	0.0008763	4.3051729

node size를 반으로 줄일 때 error가 그의 제곱인 1/4씩 줄어듦을 확인 할 수 있었다.

다음으로 fission source의 분포 경향을 살펴보았다. 앞서 중심의 assembly에서 가장 출력이 높고 그다음 가장 밖의 assembly가 가장 출력이 적었다. 출력이 가장 큰 곳과 중간 곳의 비 그리고 중간 곳과 가장 낮은 곳의 비율을 살펴보았다.

mesh 수	[eqn: max over mid]	[eqn: mid over min]
2	1.389	2.033
4	1.387	2.165
8	1.317	2.241
16	1.268	2.281

mesh 수가 늘어날수록 해는 real value에 수렴을 하는데 max/mid 값은 감소하고 mid/min 값은 증가함을 관찰하였다.

mesh의 size에 따라 결과가 수렴해가는 경향성을 보고 extrapolation을 통해 true solution을 추정할 수도 있다. 그러나 위의 결과는 생각과 다르게 완벽하게 log scale로 감소하는 경향을 보이지 않아서 추정한다 하더라도 정확하지 않을 수 있다. 추정한 값은 다음과 같다.

[eqn: k_eff]= 0.937852

3. 느낀점 및 되짚어 볼 점

Line SOR이 빠르다는 것을 잘 못 느꼈으나 후에 [eqn: Ax=b] solver 코드를 수정하여 속도를 개선시키니 월등히 Line SOR이 Point SOR보다 훨씬 빠르다는 것을 알 수 있었다.

현재의 문제만을 풀 수 있는 optimized된 코드가 아닌 임의의 matrix에 대해서도 문제를 풀 수 있도록 iterative scheme을 작성하려다보니 조금 길어지기도 하고 남보다 조금 효율이 떨어지는 것 같다는 기분도 들기도 하였다.

수업시간에 단순히 SOR iterative scheme에 대해서 배우고 method 자체는 정말 간단해 보인다고 생각했었는데 실제로 2D,2G FDM 문제를 풀기위해 코드를 작성하다보니 간단하지만 실제 작성하는 데는 많은 시행착오가 필요하다는 것을 느꼈다. 특히 코드를 볼 때마다 수정해야할 부분이 나온 것이 숙제 제출이 늦어진 가장 큰 원인이었다.

숙제를 하면서 그냥 SOR scheme을 사용하는데 있어서 그 알고리즘만 작성하면 되는 것이 아니라 그것을 사용하는데 있어서 optimize 하기위해 여러 방법이 동원되는 것도 알 수 있었다. 배운 방법 이외에도 더 좋은 방법이 있지 않을까 하는 생각도 들었다.