Neutron Diffusion Equation with Julia

We define the data for the problem, i.e.,

  • the cross sections for node;
  • the length of every node;
  • the boudnary conditions;
  • and the local refinement of every node;

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]:
power_iteration (generic function with 4 methods)

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


iter = 152;  lambda = 0.2877123282630023; error = 9.520800634218764e-11

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


iter = 149;  lambda = 0.2863104920180563; error = 9.807024387333124e-11

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="-")


INFO: Loading help data...
Out[8]:
PyObject <matplotlib.legend.Legend object at 0x7f4b72890810>

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


Maximum eigenvalue of L^{-1}*M = 0.28573162478748304
Minimum eigenvalue of L^{-1}*M = 0.0

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]:
PyObject <matplotlib.legend.Legend object at 0x7f4b57025d90>

In [1430]:
# L_fmfd

In [1431]:
# L_cmfd

In [1439]: