Untitled

import numpy as np

def LU_factor(A):
    '''
    Perform LU factorization on matrix A.
    L has implied 1s along its diagonal, and U is stored in the upper triangular part of A.
    The lower triangular part of A stores the multipliers for L.
    '''
    N = A.shape[0]

    for col in range(N):
        for row in range(col + 1, N):
            mod_row = np.copy(A[row])
            factor = mod_row[col] / A[col, col]

            mod_row -= factor * A[col, :]

            # Store the factor in the modified row
            mod_row[col] = factor
            # Store only the part of the modified row that makes U
            mod_row = mod_row[col:N]

            A[row, col:N] = mod_row

    return A

def LU_solve(A, b):
    '''
    Solve the system Ax = b using the LU factorization stored in A.
    '''
    N = A.shape[0]
    x = np.zeros(N)
    y = np.zeros(N)  # dummy variable for L^-1 b

    # Forward solve for y in Ly = b
    for row in range(N):
        RHS = b[row]
        for col in range(row):
            RHS -= y[col] * A[row, col]
        y[row] = RHS

    # Backward solve for x in Ux = y
    for row in range(N - 1, -1, -1):
        RHS = y[row]
        for col in range(row + 1, N):
            RHS -= x[col] * A[row, col]
        x[row] = RHS / A[row, row]

    return x

def inverse_power(A, B, epsilon=1e-6):
    '''
    Solve the generalized eigenvalue problem Ax = lB using the inverse power iteration algorithm.

    Args:
        A: LHS matrix (cannot be singular)
        B: RHS matrix
        epsilon: tolerance on eigenvalue

    Output:
        l: smallest eigenvalue of the problem
        x: associated eigenvector of smallest eigenvalue
    '''
    Nrows, Ncols = A.shape
    # Generate initial guess
    x = np.random.random(Nrows)
    x = x / np.linalg.norm(x)  # normalize the initial guess

    l_old = 0
    converged = False
    # Compute LU factorization of A
    A_factored = LU_factor(np.copy(A))
    iteration = 0
    b_0s = []

    while not converged:
        iteration += 1
        b = LU_solve(A_factored, np.dot(B, x))
        b_0s.append(b[0])
        l = np.linalg.norm(b)
        x = b / l

        converged = (np.fabs(l - l_old) < epsilon)
        #print(f"Iteration: {iteration}, magnitude of l: {1/l:.4f}, epsilon: {np.fabs(l - l_old):.3e}")
        l_old = l

    sign = b_0s[iteration - 1] / b_0s[iteration - 2]
    return sign / l, x

def create_grid(R, I):
    """
    Create cell edges and centers for a domain of
    size R and for I cells

    Args:
        R: size of domain
        I: number of cells

    Returns:
        Delta_r: width of each cell
        centers: cell centers of the grid
        edges: cell edges of the grid
    """
    Delta_r = float(R) / I  # divide size by # of cells
    # Calculate the centers by getting each Delta_r
    # and adding 0.5 * Delta_r
    centers = np.arange(I) * Delta_r + 0.5 * Delta_r
    # Get edges by going beyond one cell
    edges = np.arange(I + 1) * Delta_r

    return Delta_r, centers, edges
def diffusion_setup(R, I, D, Sig_a, nuSig_f, BC, geometry):
    '''
    Set up the diffusion equation in a 1-D geometry using cell-averaged quantities.

    Args:
        R: size of domain
        I: number of cells
        D: function, D(r), returns diffusion coefficient at r
        Sig_a: function, Sig_a(r), returns macroscopic absorption cross-section at r
        nuSig_f: function, nuSig_f(r), returns nu * macroscopic fission cross-section at r
        BC: boundary conditions at r=R in form [A, B, C]
        geometry:
            0 = slab
            1 = cylindrical
            2 = spherical

    Returns:
        centers: cell centers of grid
        A: LHS matrix
        B: RHS matrix
    '''
    Delta_R, centers, edges = create_grid(R, I)
    A = np.zeros((I + 1, I + 1))
    B = np.zeros((I + 1, I + 1))

    # Define surface areas S at cell edges and volumes V at dr
    if geometry == 0:  # slab
        S = np.zeros_like(edges) + 1  # 1-D slab surface area
        S[0] = 0.  # reflective BC at center
        V = np.zeros_like(edges) + Delta_R

    elif geometry == 1:  # cylinder
        S = 2. * np.pi * edges  # cylinder surface area
        V = np.pi * (edges[1:I + 1]**2 - edges[0:I]**2)  # cylinder volume

    elif geometry == 2:  # sphere
        S = 4 * np.pi * edges**2  # sphere surface area
        V = 4/3 * np.pi * (edges[1:I + 1]**3 - edges[0:I]**3)  # sphere volume

    # Setup the BC at R
    A[I, I] = (BC[0] / 2 + BC[1] / Delta_R)
    A[I, I - 1] = (BC[0] / 2 - BC[1] / Delta_R)

    # Fill A matrix
    Dplus = 0
    for i in range(I):
        r = centers[i]
        Dminus = Dplus
        Dplus = 2 * (D(r) * D(r + Delta_R)) / (D(r) + D(r + Delta_R))
        A[i, i] = 1 / (Delta_R * V[i]) * Dplus * S[i + 1] + Sig_a(r)
        B[i, i] = nuSig_f(r)

        if i > 0:
            A[i, i - 1] = -1 * Dminus / (Delta_R * V[i]) * S[i]
            A[i, i] += 1. / (Delta_R * V[i]) * Dminus * S[i]

        A[i, i + 1] = -Dplus / (Delta_R * V[i]) * S[i + 1]

    return centers, A, B

def DiffusionEigenvalue(R, I, D, Sig_a, nuSig_f, BC, geometry, epsilon=1e-8):
    '''
    Solve the diffusion eigenvalue problem.

    Args:
        R: size of domain
        I: number of cells
        D: function, D(r), returns diffusion coefficient at r
        Sig_a: function, Sig_a(r), returns macroscopic absorption cross-section at r
        nuSig_f: function, nuSig_f(r), returns nu * macroscopic fission cross-section at r
        BC: boundary conditions at r=R in form [A, B, C]
        geometry:
            0 = slab
            1 = cylindrical
            2 = spherical
        epsilon: tolerance on eigenvalue

    Returns:
        k: effective multiplication factor
        phi: eigenvector (flux)
        centers: cell centers of grid
    '''
    centers, A, B = diffusion_setup(R, I, D, Sig_a, nuSig_f, BC, geometry)
    l, phi = inverse_power(A, B, epsilon)
    k = 1 / l
    phi = phi[0:I]
    return k, phi, centers

# Define the nuclear data for Pu-239
sigma_f = 1.85  # barns
sigma_a = 2.11  # barns
sigma_tr = 6.8  # barns
nu = 2.98
density = 19.74  # g/cm^3
N_A = 6.02214076e23  # Avogadro's number
M = 239  # molar mass of Pu-239

# Calculate macroscopic cross-sections
N = density * N_A / M  # number density in atoms/cm^3
Sigma_f = sigma_f * 1e-24 * N  # macroscopic fission cross-section
Sigma_a = sigma_a * 1e-24 * N  # macroscopic absorption cross-section
Sigma_tr = sigma_tr * 1e-24 * N  # macroscopic transport cross-section
D = 1 / (3 * Sigma_tr)  # diffusion coefficient

# Define the functions for D, Sigma_a, and nuSigma_f
def D_func(r):
    if r < 2:  # inner hollow region
        return 100  # cm
    else:
        return D

def Sigma_a_func(r):
    if r < 2:  # inner hollow region
        return 0
    else:
        return Sigma_a

def nuSigma_f_func(r):
    if r < 2:  # inner hollow region
        return 0
    else:
        return nu * Sigma_f

# Boundary conditions: reflective at r=0, vacuum at r=R
BC = [1, 0, 0]  # reflective at r=0, vacuum at r=R

# Geometry: spherical
geometry = 2

# Radius of the reactor
R = 7  # cm

# Number of cells
I = 100

# Compute k_eff with the Pu slug inserted (solid sphere)
k_eff_solid, phi_solid, centers_solid = DiffusionEigenvalue(R, I, D_func, Sigma_a_func, nuSigma_f_func, BC, geometry)

# Compute k_eff without the Pu slug inserted (hollow sphere)
k_eff_hollow, phi_hollow, centers_hollow = DiffusionEigenvalue(R, I, D_func, Sigma_a_func, nuSigma_f_func, BC, geometry)

print(f"k_eff with Pu slug inserted: {k_eff_solid:.4f}")
print(f"k_eff without Pu slug inserted: {k_eff_hollow:.4f}")