When I run my program it does nothing, am trying for all print statements to print, please help!!

import math import cmath import numpy as np import matplotlib.pyplot as plt import time import pandas as pd import sys from numpy.linalg import inv import scipy.linalg from scipy.linalg import lu_factor, lu_solve

def inversePowerBlock(M11, M21, M22, P11, P12,epsilon=1.0e-6,LOUD=False): """Solve the generalized eigenvalue problem (M11 0 ) (phi_1) = l (P11 P12) using inverse power iteration (M21 M22) (phi_2) (0 0 ) Inputs Mij: An LHS matrix (must be invertible) P1j: A fission matrix epsilon: tolerance on eigenvalue Outputs: l: the smallest eigenvalue of the problem x1: the associated eigenvector for the first block x2: the associated eigenvector for the second block """ N, M = M11. shape assert (N==M) # generate initial guess x1 = mp. random. random ( (N)) x2 = mp. random. random ( (N)) l_old = np.linalog. norm (np. concatenate ( (x1, x2))) X1 = X1/l_Old x2 = x2/l_Old Converged = 0 # compute LU factorization of M11 row_order_11 = LU_factor (M11, LOUD=False) # compute LU factorization of M22 row_order_22 = LU_factor (M22, LOUD=False) iteration = 1; while not(converged): #Solve for bl b1 = LU_solve (M11, np.dot (P11, X1) + np.dot (P12, x2), row_order_11) #Solve for b2 b2 = LU_solve (M22, np.dot (-M21, b1), row_order_22) #eigenvalue estimate is norm of combined vectors l = np.linalog. norm (np. concatenate ( (b1, b2))) X1 = b1/l X2 = b2/l converged = (np.fabs (l - l_old) < epsilon) l_Old = l if (LOUD): print("Iteration: ", iteration, "\t|Magnitude of l =", 1.0/l) iteration += 1 return 1.0/l, X1, X2 print (x1, x2)

#define A M11 = np.identity (2) M11 [0,0] = 10.0 M11 [1, 1] = 0.5 M22 = np. identity (2) M22 [1, 1] = 0.1 M21 = - np. identity (2) #define P P11 = np. identity (2) P12 = np. identity (2) l, x1, x2 = inversePowerBlock (M11, M21, M22, P11, P12, epsilon=1.0e-8, LOUD=True)

def TwoGroupEigenvalue(R,I,D1,D2,Sig_r1,Sig_r2, nu_Sigf1, nu_Sigf2,Sig_s12, BC1,BC2, geometry,epsilon = 1.0e-8): """Solve a neutron diffusion eigenvalue problem in a 1-D geometry using cell-averaged unknowns Args: R: size of domain I: number of cells Dg: name of function that returns diffusion coefficient for a given r Sig_rg: name of function that returns Sigma_rg for a given r nuSig_fg: name of function that returns nu Sigma_fg for a given r Sig_s12: name of function that returns Sigma_s12 for a given r BC1: Boundary Value of fast phi at r=R in form [A,B] BC2: Boundary Value of thermal phi at r=R in form [A,B] geometry: shape of problem 0 for slab 1 for cylindrical 2 for spherical Returns: k: the multiplication factor of the system phi_fast: the fast flux fundamental mode with norm 1 phi_thermal: the thermal flux fundamental mode with norm 1 centers: position at cell centers """ #create the grid Delta_r, centers, edges = create_grid(R,I) M11 = np.zeros((I+1,I+1)) M21 = np.zeros((I+1,I+1)) M22 = np.zeros((I+1,I+1)) P11 = np.zeros((I+1,I+1)) P12 = np. zeros((I-1, I-1))

#define surface areas and Volumes assert ((geometry==0) or (geometry == 1) or (geometry == 2)) if (geometry == 0): # in Slab itS 1 everywhere except at the left edge S = 0.0 *edges + 1 S[0] = 0.0 # to enforce Ref 7 BC # in Slab its dr V = 0.0 * centers + Delta_r elif (geometry == 1): # in cylinder it is 2 pi r S = 2.0 * np.pi*edges # in cylinder its pi (r^2 - r^2) V = np.pi* ( edges [1: (I+1)] **2 - edges [0: I ] **2)

elif (geometry == 2): # in sphere it is 4 pi r^2 S = 4.0 * np.pi*edges ** k2 # in sphere its 4/3 pi (r^3 - r^3) V = 4.0/3.0 * np.pi *( edges [1: (I+1)]**3 - edges [0: I] **3 ) #Set up BC at R M11 [I, I] = (BC1 [0] *0.5 + BC1 [1]/Delta_r) M11 [I, I -1] = (BC1 [0]*0.5 - BC1 [1]/Delta_r) M22 [I, I] = (BC2[0]*0.5 + BC2 [1]/Delta_r) M22 [ I, I -1] = (BC2[0]*0.5 - BC2 [1]/Delta_r) #fi l l in rest of matrix for i in range (I): r = centers [i] M11 [i, i] = (0.5/(Delta_r * V[i])* ( (D1 (r)+D1 (r+Delta_r))*S[i+1]) +Sig_r1(r)) M22 [i, i] = (0.5/(Delta_r * V[i])* ( (D2(r)+D2 (r+Delta_r))*S[i+1]) + Sig_r2(r)) M21 [i, i] = - Sig_S12(r) P11 [i, i] = nu_Sigf1(r) P12[i, i] = nu_Sigf2(r) if (i> 0): M11 [i, i-1] = -0.5* (D1 (r)+D1 (r-Delta_r))/(Delta_r * V[i]) * S[i] M11 [i, i] += 0.5/(Delta_r * V[i]) * ( (D1 (r)+D1 (r-Delta_r))*S[i]) M22 [i, i-1] = -0.5* (D2(r)+D2 (r-Delta_r))/(Delta_r * V[i]) * S[i] M22 [i, i] += 0.5/(Delta_r * V[i]) * ( (D2(r)+D2 (r-Delta_r))* S[i]) M11 [i, i+1] = -0.5* (D1 (r)+D1 (r-Delta_r))/(Delta_r < V[i]) * S[i+1] M22 [i, i+1] = -0.5* (D2(r)+D2 (r-Delta_r))/(Delta_r < V[i]) * S[i+1] #find eigenvalue l, phi1, phi2 = inversePowerBlock (M11, M21, M22, P11, P12, epsilon) k = 1.0/l # remove last element of phi because it is outside the domain phil = phil [0: I] phi2 = phi2 [O: I] return k, phil, phi2, centers

print (k, phi1, phi2, centers) nuSigmaf_func = lambda r: 0.1570 D_func = lambda r: 3.850204978408833 Sigmaa_func = lambda r: 0.1532 Sigmas_func = lambda r: 0.0 R = 100 I = 20 #solution in spherical geometry with 100 cells k, phi_f,phi_t,centers = TwoGroupEigenvalue(R,I,D_func,D_func, Sigmaa_func,Sigmaa_func, nuSigmaf_func,nuSigmaf_func, Sigmas_func,[1,0,0],[1,0,0], 2, epsilon=1.0e-10)

print (k, phi1, phi2, centers)

nuSigmaf1_func = lambda r : 0.0085 nuSigmaf2_func =lambda r : 0.185 D_func = lambda r : 0.1 Sigmar1_func = lambda r : 0.0362 Sigmar2_func = lambda r : 0.121 Sigmar12_func = lambda r : 0.0241 R = 5 I = 50 k, phi_f, phi_t, centers = TwoGroupEigenvalue (R, I, D_func, D_func, Sigmar1_func, Sigmar2_func, nuSigmaf1_func, nuSigmaf2_func, SigmaS12_func, [0, 1,0], [0, 1, 0], 2, epsilon=1.0e-8)

print (k, phi1, phi2, centers)