Untitled

# Mesh sizes to test
mesh_sizes = np.linspace(10,400,40, dtype=int)
k_eff_values = []

# Loop over mesh sizes
for I in mesh_sizes:
    # Initialize material properties for each cell
    D1 = np.zeros(I)
    D2 = np.zeros(I)
    Sig_r1 = np.zeros(I)
    Sig_r2 = np.zeros(I)
    nu_Sig_f1 = np.zeros(I)
    nu_Sig_f2 = np.zeros(I)
    Sig_s1_to_2 = np.zeros(I)

    for i, region in enumerate(regions):
        start = i * (I // len(regions))
        end = (i + 1) * (I // len(regions))
        D1[start:end] = materials[region][0]
        D2[start:end] = materials[region][1]
        Sig_r1[start:end] = materials[region][2]
        Sig_r2[start:end] = materials[region][3]
        nu_Sig_f1[start:end] = materials[region][4]
        nu_Sig_f2[start:end] = materials[region][5]
        Sig_s1_to_2[start:end] = materials[region][6]

    # Solve the two-group diffusion equations
    phi1, phi2, centers = solve_two_group_diffusion(R, I, D1, D2, Sig_r1, Sig_r2, nu_Sig_f1, nu_Sig_f2, Sig_s1_to_2, BC, geometry)

    # Calculate Delta_R for the current mesh
    Delta_R = R / I

    #  k_eff
    k_eff = calculate_k_eff(phi1, phi2, centers, nu_Sig_f1, nu_Sig_f2, Sig_r1, Sig_r2, D1, D2, Sig_s1_to_2, Delta_R, geometry)
    k_eff_values.append(k_eff)

    print(f"Mesh size: {I}, k_eff: {k_eff:.4f}")

# Plot k_eff vs mesh size
plt.plot(mesh_sizes, k_eff_values, marker='o')
plt.xlabel('Number of Cells (I)')
plt.ylabel('k_eff')
plt.title('Convergence of k_eff with Mesh Refinement')
plt.show()