원자로 수치해석 과제 #3

: Multi-group 1-D Eigenvalue problem

with Chebyshev Acceleration

제출 일자 : 2010/04/26

2008-11833. 반영석

Part 1. Input Processing and Composition k_inf Determination

1. input processing

cross section input file (xsec.inp) 의 형식은 다음과 같다.

Line Number	Entries
1	Number of Groups (G), Number of Compositions (C)
2	Composition Number, Description
3 ~ G+2	Group Number, Diff. Coef., Removal, nu-Fission, Chi for g
G+3 ~ 2G+2	Group Number, Row of Scattering Matrix
(2G+3) ~	Repeat Lines 2 ~ (2G+2) for other compositions

MATLAB을 이용하여 위의 형식의 xsec input file을 읽는 함수를 작성하였다.

function rxsec

global comp ng

%comp contains D,sig_a, nusig_f, chi, sig_s for each composition

fid=fopen('xsec.inp');

tline = fgetl(fid);

temp = sscanf(tline, '%d %d');

ng=temp(1); %----------# of groups

nc=temp(2); %----------# of composition

for i=1:nc

tline=fgetl(fid);

% ---------- k inf ----------

comp(i).kinf=str2double(tline(11:length(tline)));

% ---------- D,sig_r, nu sig_f, chi ----------

for j=1:ng

tline=fgetl(fid);

temp=sscanf(tline,'%f');

comp(i).D(j)=temp(2);

comp(i).sig_r(j)=temp(3);

comp(i).sig_f(j)=temp(4);

comp(i).chi(j)=temp(5);

end

% ---------- scattering xsec (GxG matrix) ----------

for j=1:ng

tline=fgetl(fid);

for k=1:ng

arr=(tline(13*k-11:13*k+1));

if strcmp(arr,'             ')

tarr(k)=0;

else

tarr(k)=str2double(arr);

end

end

comp(i).scat(j,1:ng)=tarr(1:ng);

end

end

여기서 모든 input data는 comp라는 구조체 안에 저장된다. comp는 composition 종류 수만큼 있고 각 comp값 안에는 각 composition에 대한 그룹별로 diffusion coefficient [eqn: D], removal cross section [eqn: SIGMA_r], nu* fission cross section [eqn: nu Sigma_f] , Chi [eqn: chi] 값과 scattering cross section [eqn: Sigma_g'->g]값이 저장된다. scattering cross section은 scat(g,g')으로 저장된다. 또 추가적으로 [eqn: k_inf]값을 저장하여 후에 계산할 [eqn: k_inf]의 reference값을 저장하였다.

MATLAB에서 구분자를 인식하지 못하는 문제로 STR2DOUBLE 이라는 string format을 double format으로 바꾸어주는 함수와 두 string을 비교해서 일치여부를 판별해주는 STRCMP 라는 함수를 활용하였다.

다음으로 main input file을 읽는 함수를 작성하였다. main input file의 형식은 지난 HW.2 와 같고 추가적으로 Chebyshev의 최대차수와 시작 전 power iteration 수를 포함하고 있다.

Line Number	Entries
1	Number of Region (nc), [eqn: xi], Number of Compositions
2	Composition
3	Number of nodes in one region
4	Boundary Conditions ( [eqn: alpha_L] , [eqn: alpha _R] )
5	Convergence Criterion (for eigenvalue & fission source error)
6	Maximum order of Chebyshev , # of iteration before starting Chebyshev acceleration

위를 읽어 들이는 함수는 지난 숙제에서 작성한 readinp.m 파일과 거의 동일하다.

function rmain

global carr nc xi n h maxorder startnum conv_eig conv_sol al_l al_r

fid=fopen('main.inp');

tline=fgetl(fid);

temp=sscanf(tline,'%f');

nc=temp(1);

xi=temp(2);

%---------- composition array ----------

tline=fgetl(fid);

carr=sscanf(tline,'%d');

%----------------------------------------

%---------- # of nodes in region ----------

tline=fgetl(fid);

n=sscanf(tline,'%d');

h=xi/n;

%----------------------------------------

%---------- boundary condition ----------

tline=fgetl(fid);

temp=sscanf(tline,'%f');

al_l=temp(1);

al_r=temp(2);

%----------------------------------------

%---------- convergence criteria ----------

tline=fgetl(fid);

temp=sscanf(tline,'%f');

conv_eig=temp(1);

conv_sol=temp(2);

%-----------------------------------------

tline=fgetl(fid);

temp=sscanf(tline,'%f');

maxorder=temp(1);

startnum=temp(2);

2. Determine the k_inf of each fuel composition by using cross section data and compare with given reference date.

[eqn: k_inf] 값은 무한 증배계수로서 무한한 길이의 원자로에서 좌우로 빠져나가는 leakage가 없을 때의 [eqn: k]값이다. 그 값은 주어진 diffusion equation에서 좌변의 leakage term을 zero로 둔 후 [eqn: k]값을 계산하면 된다. 즉 아래의 식을 푼다.

[eqn: SIGMA _{r,g} phi _{g} = {1} over {k _{INF }} chi _{g} psi + sum _{eqalign{g prime =1# g prime != g}} ^{G} Sigma_{s,g'->g } phi_g']

initial value로 [eqn: k =1], [eqn: psi = 1]로 둔 후 [eqn: g = 1]에서부터 순차적으로 계산을 다음과 같이 한다.

[eqn: Sigma _ r,1 phi_1 = chi_1],  [eqn: phi _1 = chi_1 over Sigma_r,1]

[eqn: Sigma_r,2 phi_2 = chi_2 + Sigma_s,1->2 phi_1], [eqn: phi_2 = {chi_2 + Sigma_s,1->2 phi_1 } over {Sigma_r,2}]

[eqn: cdots]

[eqn: g=G]까지 계산을 해준 후에 [eqn: k_inf]값은 다음과 같이 구해진다.

[eqn: k_inf = psi* = SUM_g=1 ^G nu Sigma_f,g phi_g]

이를 반복적으로 계산하는 함수를 작성하였다.

function cal_kinf

global comp cc ng

rxsec();

for c=1:cc

estkinf=0;

phi(c,1)=0;

for g=1:ng

sum=0;

for i=1:g-1

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

sum=sum+comp(c).scat(g,i)*phi(c,i);

end

phi(c,g)=(comp(c).chi(g)+sum)/comp(c).sig_r(g);

estkinf=estkinf+comp(c).sig_f(g)*phi(c,g);

end

str=sprintf('%d th comp : estimated k_inf = %f , real k_inf = %f',c,estkinf,comp(c).kinf);

disp(str);

end

출력결과는 다음과 같다.

1 th comp : estimated k_inf = 1.202909 , real k_inf = 1.202909

2 th comp : estimated k_inf = 0.921402 , real k_inf = 0.921402

3 th comp : estimated k_inf = 1.128506 , real k_inf = 1.128506

4 th comp : estimated k_inf = 1.151832 , real k_inf = 1.151832

5 th comp : estimated k_inf = 1.067965 , real k_inf = 1.067965

6 th comp : estimated k_inf = 0.000000 , real k_inf = NaN

계산을 통해 estimate 한 [eqn: k_inf]값과 input file에 포함되어있는 reference [eqn: k_inf]값을 비교해보았다. 둘에는 차이가 없었다. 6번째 composition에서는 [eqn: nu Sigma_f]값이 zero이기 때문에 [eqn: k_inf]가 zero가 될 수밖에 없고 또 그렇기 때문에 input data에서도 그 값을 주지 않았기 때문에 NaN이 출력되었다.

Part 2. Fission Source Iteration

1. Theoretical Background

이번숙제에서 풀어야 하는 것은 One dimensional, Multi-group Neutron Diffusion Equation이다. 이에 해당하는 diffusion equation은 다음과 같다.

[eqn: - {d} over {dx} [D _{g} ^{k} {d phi _{g} ^{k}} over {dx} ]+ SIGMA _{t,g} ^{k} phi _{g} ^{k} = {1} over {k _{eff}} chi _{g} ^{k} [ sum _{g prime =1} ^{G} nu SIGMA _{f,g prime } ^{k} phi _{g prime } ^{k} ]+[ sum _{g prime =1} ^{G} SIGMA _{g prime -> g} ^{k} phi _{g prime } ^{k} ]]	... (2.1)

[eqn: k]는 node number를 뜻하고 증배계수는 [eqn: k_eff]로 표기하였다. 식(2.1)에서 scattering source 항에서 self scattering 항을 지워주기 위해 removal cross section[eqn: SIGMA _{r} = Sigma_t -Sigma_g->g] 을 도입하여 식을 정리하면 식(2.2)와 같다.

[eqn: - {d} over {dx} [D _{g} ^{k} {d phi _{g} ^{k}} over {dx} ]+ SIGMA _{r,g} ^{k} phi _{g} ^{k} = {1} over {k _{eff}} chi _{g} ^{k} [ sum _{g prime =1} ^{G} nu SIGMA _{f,g prime } ^{k} phi _{g prime } ^{k} ]+[ sum _{eqalign{g prime =1# g prime != g}} ^{G} SIGMA _{g prime -> g} ^{k} phi _{g prime } ^{k} ]]	... (2.2)

지난 숙제와 같이 box scheme을 적용하여 discretization 하면 위의 식은 다음과 같이 다시 변경된다.

[eqn: -l _{g} ^{k} phi _{g} ^{k-1} +d _{g} ^{k} phi _{g} ^{k} -u _{g} ^{k} phi _{g} ^{k+1} =q _{g} ^{k}]	... (2.3)

식 (2.3)과 같이 [eqn: k]번째 node는 좌우의 node에 의해 [eqn: l]과 [eqn: u]만큼 영향을 받고 있다. 그 값은 아래와 같다.

[eqn: l_g^k = 2 over h^2 {D_g^k-1 D_g^k } over {D_g^k-1 + D_g^k}]	... (2.4-1)

[eqn: u_g^k = 2 over h^2 {D_g^k D_g^k+1 } over {D_g^k + D_g^k+1}]	... (2.4-2)

[eqn: d_g^k = l_g^k +u_g^k + Sigma_r,g^k]	... (2.4-3)

[eqn: q _{g} ^{k} ={1} over {k _{eff}} chi _{g} ^{k} [ sum _{g prime =1} ^{G} nu SIGMA _{f,g prime } ^{k} phi _{g prime } ^{k} ]+[ sum _{eqalign{g prime =1# g prime != g}} ^{G} SIGMA _{g prime -> g} ^{k} phi _{g prime } ^{k} ]]	... (2.4-4)

식 (2.3)을 사용하여 linear system [eqn: M_g PHI_g = s_g] 을 구성한다. [eqn: M_g]은 총 mesh수에 의해 size가 결정된다. nc개의 region으로 구성되어있고 각 region은 n개의 node로 구성되어있으므로 총 mesh수는 [eqn: nc times n]개 이다. 앞으로 [eqn: nn=nc times n]이라고 총 mesh수로 nn이라는 값을 사용할 것이다. 따라서 [eqn: M]은 nn by nn matrix이고 또 tri-diagonal하다. diagonal 성분에는 [eqn: d]가 바로 윗줄인 upper diagonal 성분에는 [eqn: u]가 lower diagonal 성분에는 [eqn: l]이 들어간다.

이제 문제가 되는 것은 [eqn: k=1, nn]인 boundary에서 [eqn: phi^0]와 [eqn: phi^nn+1]이 존재하지 않음에 따라 [eqn: l]과 [eqn: u]가 없다는 것이 문제가 된다. Boundary condition은 main input file에서 albedo로서 주어진다. [eqn: k=1 , nn]에서는 각각 [eqn: u]와 [eqn: l]값이 다음과 같이 다르게 주어진다.

[eqn: l_g^1 = 1 over h { 2 alpha_L over 2 D_g^1 over h } over {{alpha_L }over 2 + D_g^1 over h }]	... (2.5-1)

[eqn: u _{g} ^{nn} = {1} over {h} {2 {alpha _{R}} over {2} {D _{g} ^{nn}} over {h}} over {{alpha _{R}} over {2} + {D _{g} ^{nn}} over {h}}]	... (2.5-2)

그에 따른 식(2.3)은 다음과 같이 변경된다.

[eqn: d _{g} ^{1} phi _{g} ^{1} -u _{g} ^{1} phi _{g} ^{2} =q _{g} ^{1}]	... (2.6-1)

[eqn: -l _{g} ^{nn} phi _{g} ^{nn-1} +d _{g} ^{nn} phi _{g} ^{nn} =q _{g} ^{nn}]	... (2.6-2)

[eqn: M_g]는 하나의 group에 대한 linear system이다. 전체적으로 보면 [eqn: M]은 [eqn: M_g]를 diagonal로 갖는 Box tri-diagonal matrix이다. 그러나 문제를 풀 때는 한 group에 대해서만 nn by nn matrix를 사용하여 계산하게 된다.

이렇게 구성된 linear system은 다음과 같은 방법으로 해를 구한다.

#.0 Initialize and define variables

#.1 Determine initial fission source & save L and U from M

#.2 Determine source at each node

#.3 Solve [eqn: M phi = s] by using Power method for each group

#.4 Determine new fission source

#.5 Estimate eigenvalue

#.6 Estimate pseudo error and dominance ratio

#.7 Save and update data

#.8 Error monitoring & do over 2~8

#.9 Plotting flux and errors

새로운 eigenvalue는 update를 한 fission source로부터 [eqn: k^(l) = k^(l-1) {<psi^(l) , psi^(l) >} over {<psi^(l) , psi^(l-1) >}] 구할 수 있다. pseudo error는 참값을 모를 때 현재의 error를 이전 iteration에서의 값과 차이를 구함으로서 구하는 error값이다. 이 pseudo error의 비를 이용하여 수렴속도에 비례하는 dominance ratio를 구할 수 있다.

2. Write a program that solves this problem

** code는 part.3에서의 chebyshev acceleration을 제외한 부분과 동일하므로 여기에는 포함하지 않았습니다.

출력되는 결과는 다음과 같다.

[eqn: k_eff = 1.1347788]

[eqn: sigma =9.70024E-001]

iteration = 352

eigenvalue error = 3.567147e-011

fission source error = 9.945075e-008

그림  Flux distribution

그림  Heat distribution

그림  Error

<그림 1>을 보면 fast neutron들은 주로 composition 2~6th [ 4 4 4 3 3 ] 사이에서 분포하고 thermal neutron들은 그 밖의 구간에서 fast neutron보다는 낮은 flux를 보이고 특히 2~6th 이외의 구간에서 flux가 크다.

input data를 보면 6번 composition은 fission 가능한 fuel이 아님을 알 수 있다. [eqn: nu Sigma_f]가 존재하지 않기 때문이다. 따라서 moderator로 추정되는 6번 composition으로 구성된 양쪽 끝 구간에서는 thermal neutron flux만 급격히 존재하고 heat은 거의 zero이다. 5번 composition은 [eqn: k_inf]값이 상대적으로 작기 때문에 flux가 주변지역에 비해서 낮아졌다고 생각된다.

<그림 2>는 fission source를 plot한 것이다. 그를 통해 fission 의한 열 분포를 볼 수 있다. 조성이 바뀌는 구간에서 불연속적인 값을 가지는 것을 관찰 할 수 있었다.

Part 3. Chebyshev Acceleration

1. Theoretical Background

일반적인 Power method를 사용하면 error는 iteration 마다 dominance ratio [eqn: sigma = {lambda_2 over lambda_1]값 만큼 감소하게 된다. 이 경우에 [eqn: sigma]값이 작다면 수렴속도는 무지하게 느리다. 다른 방법을 사용하여 수렴속도를 증가시켜야하는데 지난번의 경우 Wielandt Shift를 사용하였는데 이번 문제에서는 Wielandt Shift를 사용하게 되면 matrix system이 변경되어 문제를 풀기에 유리했던 tri-diagonal하던 성질이 사라지게 된다. 따라서 Wielandt shift를 사용하는 것보다는 다른 방법을 사용하는 것이 좋다.

이번 HW.3에서는 Chebyshev acceleration을 하게 된다. 이 방법은 간략하게 요약하자면 extrapolation을 Chebyshev polynomial을 사용하여 최대 error들을 최소화 시키는 것이다. Two parameter method라는 두 변수를 사용하여 외삽하는 방법은 다음 식을 바탕으로 진행된다.

[eqn: x ^(k) = alpha^(k) x_pow^(k) + (1- alpha^(k) + beta^(k) ) x^(k-1) - beta^(k) x^(k-2)]	... (3.1)

여기서 [eqn: x_pow^(k)]는 [eqn: x^(k-1)]에 [eqn: M]을 곱하고 scaling하여 power method를 사용하여 구한 해이다. 이는 다음과 같이 표현된다.

[eqn: x_pow ^(k) = 1 over lambda^(k-1) A x^(k-1)]	... (3.2)

여기서 [eqn: x]는 [eqn: A]의 eigen vector들의 선형적인 합으로 표현될 수 있다.

[eqn: x^(0) = sum_i=1 ^ n c_i u_i]	... (3.3)

식 (3.3)에서 [eqn: A]를 곱하게 되면 다음과 같이 되고

[eqn: x^(1) = sum_i=1 ^n c_i ( alpha^(1) {lambda_i } over lambda^(0) + (1-alpha^(k) ) ) u_i =sum_i=1^n c_i p^(1) (lambda_i ) u_i]	... (3.4)

따라서 [eqn: A]를 계속 곱하게 되면 [eqn: x^(k)]는 다음과 같이 표현된다.

[eqn: x^(k) = sum_i=1^n c_i p^(k) (lambda_i ) u_i]	... (3.5)

여기서 Chebyshev polynomial을 사용하기 위해서 eigenvalue를 [eqn: gamma]라는 값으로 다음과 같이 치환을 해준다.

[eqn: gamma_i = 2 lambda_i over lambda_2 -1]	... (3.6)

이를 사용하면 식 (3.5)는 [eqn: lambda]에 대한 다항식에서 [eqn: gamma]에 대한 다항식으로 바뀌게 된다. 이를 다시 식 (3.1)에 식 (3.2)를 대입하여 바꾸면 다음과 같이 표현된다.

[eqn: x^(k) = sum_i=1 ^n c_i u_i left[ alpha^(k) sigma {1+gamma_i } over 2 +1 -alpha^(k) + beta^(k) right] eta^(k-1) (gamma_i ) -beta^(k) eta^(k-2) (gamma_i ) = sum_i=1 ^n c_i u_i eta^(k) (gamma_i )]  ... (3.7)

여기서 [eqn: eta^(k)]는 [eqn: gamma_i]에 대한 [eqn: k]차 polynomial이다. 이 polynomial이 Chebyshev polynomial이 되도록 하는 것이 목표이다. 그렇게 하면 error의 maximum 항이 Chebyshev polynomial의 특성에 의해 minimum이 되게 된다.

normalized된 Chebyshev polynomial은 다음과 같은 recurrence relation을 만족한다.

[eqn: {tilde T}_k (gamma) = 2 gamma {T_k-1 (gamma_1 )} over {T_k (gamma_1 ) }{tilde T}_k-1 (gamma) - {T_k-2 (gamma_1 ) }over {T_k (gamma_1 ) } {tilde T } _k-2 (gamma)]	... (3.8)

식 (3.7)을 식 (3.8)에 비교하여 계수를 맞춰주게 된다. 계산을 통해 extrapolating parameter [eqn: alpha], [eqn: beta]를 결정하면 다음과 같다.

[eqn: alpha^(k) = 4 over sigma {T_k-1 (gamma_1 ) }over {T_k (gamma_1 )]	... (3.9-1)

[eqn: beta^(k) = { T_k-2 (gamma_1 ) } over {T_k (gamma_1 ) }]	... (3.9-2)

여기서 우리는 [eqn: sigma]값을 알아야 하는데 Chebyshev acceleration을 하면 pseudo error로부터 dominance ratio를 구할 수 없다. 따라서 어느 정도 power method를 사용하여 [eqn: sigma]값을 추정하고 iteration 하면서 [eqn: xi]가 1보다 커져서 [eqn: sigma ']가 [eqn: sigma*]를 벗어나게 되면 다음 방법으로 [eqn: sigma]를 보정해주는 것이 필요하다.

[eqn: zeta^(l) = {tildeT }_l-1 (gamma_2 ' ) = { T_l-1 (gamma_2 ') } over { T_l-1 (gamma_1 )]	... (3.10)

여기서 [eqn: xi = zeta^(l) T_l-1 (gamma_1 ) = zeta^(l) over zeta*^(l-1)] 값을 사용하여 [eqn: gamma_2 ']을 구하고 그 값으로부터 새로운 [eqn: sigma']값을 얻는다. 식 (3.6)의 관계로부터 다음 식 (3.11)이 유도된다.

[eqn: lambda_2 ' = lambda_2 {1+gamma_2 ' } over 2]	... (3.11-1)

[eqn: sigma' = sigma { 1+gamma_2 ' } over 2]	... (3.11-2)

식 (3.11-2)는 식 (3.11-1)를 양변을 [eqn: lambda_1]로 나누어 줌으로서 구할 수 있다.

이러한 전 과정을 통하여 다음과 같은 방법으로 Chebyshev acceleration을 하게 된다.

#.0 Initialize and define variables

#.1 Determine initial fission source & save L and U from M

#.2 Determine source at each node

#.3 Solve [eqn: M phi = s] by using Power method for each group

#.4 Determine new fission source

##.0 If (condition satisfy)

condition: [eqn: xi]>1 or chebyshev order > chebyshev maximum order

##.1 reset Chebyshev

replace [eqn: sigma'] , reset order

##.0 else (condition doesn't satisfy) do Chebyshev acceleration at least 3

###.0 if chebyshev order ==1 ; do single parameter method

##.1 Determine extrapolation parameters by two parameter method

##.2 Perform fission source extrapolation

#.5 Estimate eigenvalue

#.6 Estimate pseudo error and dominance ratio

#.7 Save and update data

#.8 Error monitoring & do over #2~8

#.9 Plotting flux and errors

위의 붉은색으로 표시된 부분이 part.2에서와 다르게 Chebyshev acceleration에서 추가된 부분이다. 식 (3.5)에서와 같이 power iteration을 할 때마다 Chebyshev의 order가 한 차씩 올라가게 되는데 그렇게 되면 [eqn: gamma_1]값이 조금의 오차가 생겨도 실제 [eqn: phi]에는 큰 오차를 야기할 수 있기 때문에 maximum order를 정해주고 그 order를 넘어갈 때마다 reset해주는 방법을 사용하였다. 또 actual error [eqn: zeta^(l)]이 theoretical error [eqn: zeta*^(l-1)]보다 크게 되는 경우에는 해가 extrapolation 되어 real solution을 지나쳐  over extrapolation 된 경우라 이 때에도 reset을 해주게 된다. 이 모든 것을 종합하여 solve 해주는 함수를 작성하였다.

2. Write a program that solves this problem

%function solve_2

clear all

load reference

global comp ng

global carr nc n maxorder startnum conv_eig conv_sol

global rpsi

global L U

rxsec();rmain();

nn=n*nc;

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);

for g=1:ng

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

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

end

iter=0; order=1-startnum;

M=con_M();

% ---------- #.1 determine fission source ----------

for c=1:nc

for i=1:n

k=(c-1)*n+i;

psi_o(k)=0;

for g=1:ng

psi_o(k)=psi_o(k)+comp(carr(c)).sig_f(g)*phi_o{g}(k);

end

end

end

% L=cell(ng);% U=cell(ng);

% for i=1:ng

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

% end

load lulu

while e_eig>conv_eig || e_sol>conv_sol

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

for g=1:ng

for c=1:nc

for i=1:n

k=(c-1)*n+i;

sum=0;

for gg=1:ng

if(gg==g) continue end

% scattering source

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

sum=sum+comp(carr(c)).scat(g,gg)*phi_n{gg}(k);

end

%total source term

s(k)=lam_o*comp(carr(c)).chi(g)*psi_o(k)+sum;

end

end

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

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

phi_n{g}=solve_axb(L{g},U{g},s');

end

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

for c=1:nc

for i=1:n

k=(c-1)*n+i;

psi_n(k)=0;

for g=1:ng

psi_n(k)=psi_n(k)+comp(carr(c)).sig_f(g)*phi_n{g}(k);

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;

order=order+1;

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

e_r(iter)=e_ref;    e_rr(iter)=e_res;

r(iter)=rho;

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

arr=sprintf('iter = %d ,k-eff = %7.5f rho = %e , error_1 = %e , error_2 = %e , error_3 = %e , xi = %f , r = %f',iter,k_o,sig,e_eig,e_sol,e_ref,xi,gamma);

disp(arr);

%if(e_sol < 1.0E-3) break; end

end

draw(phi_o);

x=1:nn;

y=1:iter;

hold off

fig2=figure(2);

fig2=plot(x,psi_o);

legend('heat graph');

fig3=figure(3);

fig4=semilogy(e_e);

legend('eigenvalue error');

fig4=figure(4);

fig4=semilogy(e_s);

legend('fission source error');

fig5=figure(5);

fig5=semilogy(y,e_e,y,e_s,y,e_rr,y,e_r);

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

grid on

% save lulu L U

[eqn: k_eff = 1.1347788]

[eqn: sigma =9.70024E-001]

iteration = 66

eigenvalue error = 3.236389e-011

fission source error = 8.545804e-008

reference error = 1.166507e-007

그림   flux

그림  heat distribution

그림  error reduction

<그림 4,5> 는 앞에서 power method로 구한 solution과 동일하였다. 다만 error의 경향이 다소 차이가 있었는데 <그림 6>에서는 error를 10E-12까지 줄여서 구한 fission source로부터 계산한 reference error와 residual error, 그리고 eigenvalue error와 fission source error가 있다. residual error와 reference error는 비슷한 경향으로 줄어감을 알 수 있고 fission source error는 acceleration이 reset 될 때마다 꺾임을 볼 수 있었다. 40번째 iteration 이후에는 dominance ratio [eqn: sigma]가 정확하게 추정되어 해가 완만하게 잘 수렴하감을 볼 수 있었다. 이는 reference error의 기울기가 더욱 급해진 것으로부터 알 수 있기도 하다.

총 iteration 횟수는 66번으로 acceleration을 하지 않은 경우보다 대략 6배 빨랐다.

3. Examine various orders of Chebyshev polynomial

Acceleration을 하기 전에 [eqn: n_0]번의 power iteration을 하게 되는데 이 횟수를 3번에서 10번으로 늘려보았다.

그림  error [eqn: n_0]=10

<그림 7>로부터 10번째 iteration 까지는 완만하게 error가 감소하다가 acceleration에 의해 error가 많이 흔들리는 모양을 보인다. 그러나 총 iteration 횟수에는 변화가 없이 66번 진행되었다.

[eqn: n_0 = 1]인 경우에는 그래프의 모양은 위와 동일했지만 iteration 횟수가 2번 증가한 68번이 되었다. [eqn: n_o]값에 따라서 iteration 횟수에는 크게 변화하지 않았다.

maximum order를 변화시켜보려 했으나 maximum order를 넘어서 reset 되는 경우는 없었기 때문에 구하고자 하는 error criterion 범위에서는 관찰을 할 수가 없었다.

그림  error when initial [eqn: rho =0.9700]

다음으로는 초기 [eqn: rho]([eqn: =sigma]) 추정값을 Power method가 아닌 0.9700의 값을 주고 수렴이 어떻게 되어 가는지 관찰하였다. reset되는 수가 줄어들었고 또 iteration 수도 50회로 감소하였다. 그러나 reset 되면서 [eqn: sigma]값은 reference 값인 0.970024에서는 조금씩 멀어져갔다.

추가적으로 iteration 수에 따라서 fission source의 분포가 어떻게 수렴해 가는지를 관찰해 보았다.

그림  iteration = 1

그림  iteration = 2

그림  iteration = 10

그림  iteration = 300

교수님이 말씀하신 cosine shape를 관찰 할 수는 없었다. 두 번째 iteration에서도 불연속적인 모양이 관찰되었기 때문이다.

결론적으로는 Chebyshev acceleration을 사용하여 Power method만을 사용할 때 수렴속도가 느린 점을 보완하여 약 6배 정도의 수렴속도 증가를 관찰할 수 있었으나 초기 추정 [eqn: sigma]값은 총 iteration 횟수에 큰 영향을 주지 않음을 확인 할 수 있었다.

< 부록 - 코드 >

function M=con_M

rmain();rxsec();

global carr nc n h al_l al_r comp ng

M=cell(ng,1);

for g=1:ng

nn=n*nc;

%block matrix

B=sparse(nn);

for c=1:nc

temp=comp(carr(c)).D(g)/h^2;

for j=1:n-1

k=(c-1)*n+j;

l(k+1)=temp;

u(k)=temp;

end

end

for c=1:nc-1

k=n*c;

l(k+1)=2*(comp(carr(c)).D(g)*comp(carr(c+1)).D(g))/(h^2 * (comp(carr(c)).D(g)+comp(carr(c+1)).D(g)));

u(k)=l(k+1);

end

l(1)=2*(al_l*comp(carr(1)).D(g))/(h*al_l+2*comp(carr(1)).D(g))/h;

u(nn)=2*(al_r*comp(carr(nc)).D(g))/(h*al_r+2*comp(carr(nc)).D(g))/h;

for c=1:nc

for j=1:n

k=(c-1)*n+j;

d(k)=l(k)+u(k)+comp(carr(c)).sig_r(g);

end

end

for k=2:nn-1

B(k,k-1)=-l(k);

B(k,k)=d(k);

B(k,k+1)=-u(k);

end

B(1,2)=-u(1);    B(1,1)=l(1)+u(1)+comp(carr(1)).sig_r(g);

B(nn,nn-1)=-l(nn);    B(nn,nn)=l(nn)+u(nn)+comp(carr(nc)).sig_r(g);

M{g}=B;

end

M을 구성한 matrix입니다.