Untitled

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
pd.set_option('display.max_rows', None)
pd.set_option('display.max_columns', None)

class NeutronFluxSimulation:
    def __init__(self, E_a, D, R_1, s0, n):
        self.E_a = E_a  # cm^-1
        self.D = D  # cm
        self.R_1 = R_1  # cm
        self.s0 = s0  # Source strength
        self.n = n  # Number of nodes

        self.r = np.linspace(0, self.R_1, self.n + 1)
        self.dr = self.r[1] - self.r[0]
        self.A = np.zeros((self.n + 1, self.n + 1))
        self.b = np.zeros(self.n + 1)

    def setup_matrices(self):
        # Left boundary
        r_half = (self.r[0] + self.dr / 2)
        r_n_min = self.r[self.n] - (self.dr / 2)
        self.A[0, 0] = ((self.D * r_half) / self.dr) + (self.E_a * 0.5 * np.square(r_half))
        self.A[0, 1] = (-self.D * r_half) / self.dr
        self.b[0] = self.s0 * 0.5 * np.square(r_half)

        # Right boundary
        self.A[self.n, self.n] = (self.r[self.n] / 2) + ((self.D * r_n_min) / self.dr) + (
                    self.E_a * 0.5 * (np.square(self.r[self.n]) - np.square(r_n_min)))
        self.A[self.n, self.n - 1] = (-self.D * r_n_min) / self.dr
        self.b[self.n] = self.s0 * 0.5 * (np.square(self.r[self.n]) - np.square(r_n_min))

        # Interior nodes
        for i in range(1, self.n):  # from 1 to n-1
            r_i_min = self.r[i] - (self.dr / 2)
            r_i_plus = self.r[i] + (self.dr / 2)
            self.A[i, i] = ((self.D * (r_i_plus + r_i_min) / self.dr)) + (
                        self.E_a * 0.5 * (np.square(r_i_plus) - np.square(r_i_min)))
            self.A[i, i - 1] = -self.D * r_i_min / self.dr
            self.A[i, i + 1] = -self.D * r_i_plus / self.dr
            self.b[i] = self.s0 * 0.5 * (np.square(r_i_plus) - np.square(r_i_min))

    def gauss_seidel(self, tol=1e-6, max_iterations=100000):
        n = len(self.b)
        x = np.zeros_like(self.b)

        for iteration in range(max_iterations):
            x_old = np.copy(x)
            for i in range(n):
                sum1 = np.dot(self.A[i, :i], x[:i])
                sum2 = np.dot(self.A[i, i + 1:], x_old[i + 1:])
                x[i] = (self.b[i] - sum1 - sum2) / self.A[i, i]

            # Check for convergence
            if np.linalg.norm(abs(x - x_old), ord=np.inf) < tol:
                print(f'Gauss-Seidel: Converged after {iteration + 1} iterations')
                return x

        print('Max iterations reached without convergence')
        return x

    def jacobi(self, tol=1e-6, max_iterations=100000):
        n = len(self.b)
        x = np.zeros_like(self.b)
        x_new = np.zeros_like(self.b)

        for iteration in range(max_iterations):
            for i in range(n):
                sum1 = np.dot(self.A[i, :i], x[:i])
                sum2 = np.dot(self.A[i, i + 1:], x[i + 1:])
                x_new[i] = (self.b[i] - sum1 - sum2) / self.A[i, i]

            # Check for convergence
            if np.linalg.norm(x_new - x, ord=np.inf) < tol:
                print(f'Jacobi: Converged after {iteration + 1} iterations')
                return x_new

            x = np.copy(x_new)

        print('Max iterations reached without convergence')
        return x_new

    def calculate_rates(self):
        fi_gauss = self.gauss_seidel()
        fi_jacobi = self.jacobi()

        EaM = np.zeros(self.n + 1)
        EaM[0] = self.E_a * np.pi * (self.r[0] + self.dr / 2) ** 2 * fi_gauss[0]
        EaM[self.n] = self.E_a * np.pi * fi_gauss[self.n] * (self.r[self.n] ** 2 - (self.r[self.n] - self.dr / 2) ** 2)

        for i in range(1, self.n):
            EaM[i] = self.E_a * np.pi * ((self.r[i] + self.dr / 2) ** 2 - (self.r[i] - self.dr / 2) ** 2) * fi_gauss[i]

        absorption_rate = np.sum(EaM)


        # Calculate the leakage rate at the boundary
        leakage_rate = (fi_gauss[self.n] / 2) * 2 * np.pi * self.R_1
        source_rate =  self.s0 * np.pi * self.R_1 ** 2
        ab_leak = absorption_rate + leakage_rate
        data = {"Leakage Rate": [leakage_rate],
                "Absorption Rate": [absorption_rate],
                "Source Rate": [source_rate],
                "Absorption Rate + Leakage Rate": [ab_leak]
        }
        index_labels = ["Rate Values"]
        df = pd.DataFrame(data, index=index_labels)
        print(df)
        return fi_gauss, fi_jacobi

    def plot_results(self, fi_gauss, fi_jacobi):
        plt.figure(figsize=(10, 6))
        plt.plot(self.r[:(self.n + 1)], fi_gauss[:(self.n + 1)], label='Neutron Flux (Gauss-Seidel)', color='blue')
        plt.xlabel('x (cm)')
        plt.ylabel('Flux (fi)', labelpad=15)
        plt.title('Neutron Flux Distribution using Gauss-Seidel Method')
        plt.grid(True)
        plt.legend()
        plt.tight_layout()
        plt.show()

        plt.figure(figsize=(10, 6))
        plt.plot(self.r[:(self.n + 1)], fi_jacobi[:(self.n + 1)], label='Neutron Flux (Jacobi)', color='black')
        plt.xlabel('x (cm)')
        plt.ylabel('Flux (fi)', labelpad=15)
        plt.title('Neutron Flux Distribution using Jacobi Method')
        plt.grid(True)
        plt.legend()
        plt.tight_layout()
        plt.show()



simulation = NeutronFluxSimulation(E_a=0.1, D=1.2, R_1=40, s0=10 ** 4, n=50)
simulation.setup_matrices()
fi_gauss, fi_jacobi = simulation.calculate_rates()
simulation.plot_results(fi_gauss, fi_jacobi)