We define the data for the problem, i.e.,
In [1]:
# --------------------------------------------------
# lenght of the reactor
L = 18
# Cross sections
# D Sigma_a nu Sigma_f Delta_x Ref
xs = [1/(3*0.416667) 0.334 1.0 0.178 2.4 5 ;
1/(3*0.370370) 0.334 1.0 0.000 2.7 5 ]
# Boundary conditions
# bc = 0 implies zero flux boundary condition
# bc = 1 implies reflective boundary condition (zero current)
bc_ini = 0
bc_end = 0
# Defining the material of each cell
mat_coarse = [2 1 2 1 2 1 2]
n_coarse = length(mat_coarse)
# Defining the number of refinements for each cell (does not need to be the same for all cells).
ref_coarse = fill(20, n_coarse)
;
Then we perform the refinement of the nodes, and resample the cross sections to the new fine mesh structure.
In [2]:
h = []
diff = []
siga = []
nu = []
sigf = []
x = []
# --------------------------------------------------
# we set up the coarse cross sections
diff_coarse = [xs[i,1] for i in mat_coarse]
siga_coarse = [xs[i,2] for i in mat_coarse]
nu_coarse = [xs[i,3] for i in mat_coarse]
sigf_coarse = [xs[i,4] for i in mat_coarse]
h_coarse = [xs[i,5] for i in mat_coarse]
for i = 1:n_coarse
h = [h,fill(h_coarse[i]/ref_coarse[i],ref_coarse[i])]
diff = [diff,fill(diff_coarse[i],ref_coarse[i])]
siga = [siga,fill(siga_coarse[i],ref_coarse[i])]
nu = [nu,fill(nu_coarse[i],ref_coarse[i])]
sigf = [sigf,fill(sigf_coarse[i],ref_coarse[i])]
if i == 1
x_ini = 1/ref_coarse[i]*h_coarse[i]/2
else
x_ini = sum(h_coarse[1:i-1])+1/ref_coarse[i]*h_coarse[i]/2
end
x = [x,x_ini+[0:ref_coarse[i]-1]/ref_coarse[i]*h_coarse[i]]
end
n = sum(ref_coarse)
;
Functions to build the matrices for the fine mesh finite difference (fmfd) method
In [3]:
#L = -d_{i-1}+2d_i-d_{i-1}/(h^2) + siga_i
function build_L_fmfd(diff,siga,h,n)
L = zeros(n,n);
L[1,1 ] = (2-bc_ini)*diff[1]/(h[1]^2) + siga[1]
L[1,2 ] = -diff[2]/(h[1]^2)
for i = 2:n-1
L[i,i-1] = -diff[i-1]/(h[i]^2)
L[i,i ] = 2diff[i ]/(h[i]^2)+siga[i]
L[i,i+1] = -diff[i+1]/(h[i]^2)
end
L[n,n-1] = -diff[n-1]/(h[n]^2)
L[n,n ] = (2-bc_end)*diff[n]/(h[n]^2) + siga[n]
return L
end
#M = nu_i*sigf_i
function build_M_fmfd(nu,sigf)
return diagm(nu.*sigf)
end
;
Functions to build the matrices for the coarse mesh finite difference (cmfd) method
In [4]:
function build_diff_hat(diff,h,n)
# initialize the vector
diff_hat = zeros(n+1);
# first node with boundary conditions
if (bc_ini == 0)
diff_hat[1] = 2diff[1]/(h[1]+diff[1]*h[1])
elseif (bc_ini == 1)
diff_hat[1] = 0
else
diff_hat[1] = 0 # insert albedo here
end
for i = 2:n
diff_hat[i] = 2diff[i-1]*diff[i]/(diff[i-1]*h[i-1]+diff[i]*h[i])
end
# last node with boundary conditions
if (bc_end == 0)
diff_hat[n+1] = 2diff[n]/(h[n]+diff[n]*h[n])
elseif (bc_end == 1)
diff_hat[n+1] = 0
else
diff_hat[n+1] = 0 # insert albedo here
end
# return
return diff_hat
end
function build_L_cmfd(diff_hat,siga,h,n)
# initialize matrix
L = zeros(n,n);
# first node with boundary conditions
L[1,1] = (diff_hat[1]+diff_hat[2 ])/h[1] + siga[1]
L[1,2] = -diff_hat[2]/h[1]
# internal nodes
for i = 2:n-1
L[i,i-1] = -diff_hat[i ]/h[i]
L[i,i ] = (diff_hat[i ]+diff_hat[i+1])/h[i]+siga[i]
L[i,i+1] = -diff_hat[i+1]/h[i]
end
# last node with boundary conditions
L[n,n-1] = -diff_hat[n]/h[n]
L[n,n ] = (diff_hat[n]+diff_hat[n+1])/h[n] + siga[n]
return L
end
function build_M_cmfd(nu,sigf)
return diagm(nu.*sigf)
end
;
Function to apply the power iteration method for a given set of matrices L, M
In [5]:
function power_iteration(L,M,x0,tol=1.e-10,max_it=1000, verbose = false)
LU = lufact(L)
error = 1
it = 0
phi_new = x0
lambda_new = norm(M*phi_new)
while error > tolerance && it < max_it
# setting the variable to the old name
it = it+1
phi_old = phi_new
lambda_old = lambda_new
tmp = M*phi_old
phi_new = LU[:U]\(LU[:L]\(tmp))
#lambda_new = (phi_new⋅tmp)/(phi_old⋅tmp)
lambda_new = norm(phi_new)
phi_new = phi_new/lambda_new
error = max(norm(phi_new-phi_old),abs(lambda_new - lambda_old)/lambda_new)
if (mod(it,5) == 0 && verbose)
println("iter = ", it, "; lambda = ",lambda_new, "; error = ", error)
end
end
println("iter = ", it, "; lambda = ",lambda_new, "; error = ", error)
return (lambda_new,phi_new)
end
#function julia_power_iteration(L,M,x0,tol_=1.e-10,max_it=1000)
# LU = lufact(L)
# function A(x)
# tmp = M*x
# y = LU[:U]\(LU[:L]\tmp)
# return y
# end
# (lambda,phi,:) = eigs(A; nev=1,which="LM",tol=tol_, maxiter=max_it, v0=x0)
# return (lambda,phi)
#end
Out[5]:
We apply the power iteration to the fmfd matrices
In [6]:
L_fmfd = build_L_fmfd (diff,siga,h,n)
M_fmfd = build_M_fmfd(nu,sigf)
x0 = 1+rand(n) #fill(1.0,n)
tolerance = 1.e-10
max_it = 1000
(lambda_fmfd,phi_fmfd) = power_iteration(L_fmfd,M_fmfd,x0,tolerance,max_it)
;
#using PyPlot
#plot(x, phi_fmfd, color="red", linewidth=2.0, linestyle="--", "o")
We apply the power iteration to the cmfd matrices
In [7]:
diff_hat = build_diff_hat(diff,h,n)
L_cmfd = build_L_cmfd (diff_hat,siga,h,n)
M_cmfd = build_M_cmfd(nu,sigf)
x0 = 1+rand(n) #fill(1.0,n)
tolerance = 1.e-10
max_it = 1000
(lambda_cmfd,phi_cmfd) = power_iteration(L_cmfd,M_cmfd,x0,tolerance,max_it)
;
#using PyPlot
#plot(x, phi_cmfd, color="blue", linewidth=2.0, linestyle="--", "o")
In [8]:
using PyPlot
plot(x, phi_fmfd, color="red", linewidth=3.0, linestyle="--", "o")
plot(x, phi_cmfd, color="blue", linewidth=3.0, linestyle="--", "x")
legend(["phi_1 FMFD", "phi_1 CMFD"], loc="upper left")
# phi_ref = sin(x*pi/L)/norm(sin(x*pi/L))
# plot(x, phi_ref, color="black", linewidth=1.0, linestyle="-")
Out[8]:
In [9]:
# different built-in function to solve the same
# considering a standard eigenvalue problem
lambda_max = eigmax(inv(L_cmfd)*M_cmfd)
lambda_min = eigmin(inv(L_cmfd)*M_cmfd)
println("Maximum eigenvalue of L^{-1}*M = ", lambda_max)
println("Minimum eigenvalue of L^{-1}*M = ", lambda_min)
# considering a generalized eigenvalue problem
#(lambda,phi,:) = eigs(M_fmfd,L_fmfd; nev=1,tol=1.e-6)
verbose = false
if (verbose)
LM = inv(L_cmfd)*M_cmfd;
(lambda_eigs,phi,:) = eigs(LM; nev=1,tol=1.e-6)
phi = phi/sum(phi)
phi = phi/norm(phi)
println("Generalized eigenvalue of L,M = ", lambda_eigs[1])
using PyPlot
figure, plot(x, phi, color="blue", linewidth=2.0, linestyle="--", "o")
end
In [10]:
(all_fmfd, phi_fmfd) = eig(M_fmfd,L_fmfd)
(all_cmfd, phi_cmfd) = eig(M_cmfd,L_cmfd)
perm_fmfd = sortperm(all_fmfd, by=abs, rev=true)
phi_fmfd = phi_fmfd[:,perm_fmfd[1]]
phi_fmfd = phi_fmfd/sum(phi_fmfd)
phi_fmfd = phi_fmfd/norm(phi_fmfd)
perm_cmfd = sortperm(all_cmfd, by=abs, rev=true)
phi_cmfd = phi_cmfd[:,perm_fmfd[1]]
phi_cmfd = phi_cmfd/sum(phi_cmfd)
phi_cmfd = phi_cmfd/norm(phi_cmfd)
subplot2grid((4, 1), (0, 0), rowspan=2)
plot(x, phi_fmfd, color="red", linewidth=2.0, linestyle="--", "o")
plot(x, phi_cmfd, color="blue", linewidth=2.0, linestyle="--", "x")
legend(["phi_1 FMFD", "phi_1 CMFD"], loc="upper left")
#phi_ref = sin(x*pi/L)/norm(sin(x*pi/L))
#plot(x, phi_ref, color="black", linewidth=1.0, linestyle="-")
# ploting the spectrum
all_fmfd = eig(M_fmfd,L_fmfd)[:1]
all_cmfd = eig(M_cmfd,L_cmfd)[:1]
subplot2grid((4, 1), (2, 0), rowspan=1)
semilogy(all_fmfd, color="red", "o")
semilogy(all_cmfd, color="blue", "x")
legend(["lambda_1 FMFD", "lambda_1 CMFD"], loc="upper left")
subplot2grid((4, 1), (3, 0), rowspan=1)
semilogy(-all_fmfd, color="red", "o")
semilogy(-all_cmfd, color="blue", "x")
legend(["-lambda_1 FMFD", "-lambda_1 CMFD"], loc="lower left")
#println((all_fmfd',all_cmfd'))
Out[10]:
In [1430]:
# L_fmfd
In [1431]:
# L_cmfd
In [1439]: