PyNE work done for Rachel Slaybaugh at UC Berkeley. http://pyne.io/

steady-state.py 22KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838
  1. # %THIS CODE SOLVES THE STEADY-STATE, MONOENERGETIC TRANSPORT EQUATION
  2. # %IN A HOMOGENEOUS MEDIUM WITH ISOTROPIC SCATTERING AND ISOTROPIC SOURCE
  3. # %IN ROD AND SLAB (SN) GEOMETRIES. IT USES 2-POINT CENTRAL DIFFERENCING
  4. # %(ORDER H^2) WITH 3-POINT FORWARD/BACKWARD AT THE BOUNDARIES.
  5. # %MAIN OUTPUTS: transmission "TR"; reflection "RL"
  6. # % scalar flux "SCAL"
  7. # %RICHARD VASQUES & NITIN KUMAR YADAV
  8. import numpy as np
  9. import matplotlib.pyplot as plt
  10. import os, sys
  11. import random
  12. def config_file(): #only supports homogenous right now. When refactoring code to clean out repeat functions, make generic
  13. folder = os.path.join(os.getcwd(), os.path.dirname(__file__))
  14. config = os.path.join(folder, "steady-state.config")
  15. result = []
  16. #add error checking later
  17. with open(config, 'rb') as f:
  18. for line in f:
  19. line = line.strip()
  20. if len(line) == 0:
  21. continue
  22. if line[0] == "%":
  23. continue
  24. line = line[1:len(line)-1].split(",")
  25. if len(line) == 9:
  26. if line[0]:
  27. if line[0] == "SLAB":
  28. rod_slab = 1
  29. rod_slab = 0
  30. if line[1]:
  31. if rod_slab:
  32. N = int(line[1])
  33. else:
  34. N = 2
  35. if line[2]:
  36. T = float(line[2])
  37. if line[3]:
  38. n = int(line[3])
  39. if line[4]:
  40. et = float(line[4])
  41. if line[5]:
  42. cc = float(line[5])
  43. es = cc * et
  44. if line[6]:
  45. yo = float(line[6])
  46. if line[7]:
  47. y_ = float(line[7])
  48. if line[8]:
  49. Q = float(line[8])
  50. elif len(line) == 12:
  51. if line[0]:
  52. T = float(line[0])
  53. if line[1]:
  54. n = float(line[1])
  55. if line[2]:
  56. yo = float(line[2])
  57. if line[3]:
  58. y_ = float(line[3])
  59. if line[4]:
  60. m1 = float(line[4])
  61. if line[5]:
  62. m2 = float(line[5])
  63. if line[6]:
  64. Et1 = float(line[6])
  65. if line[7]:
  66. cc = float(line[7])
  67. if line[8]:
  68. Q1 = float(line[8])
  69. if line[9]:
  70. Et2 = float(line[9])
  71. if line[10]:
  72. cc2 = float(line[10])
  73. if line[11]:
  74. Q2 = float(line[11])
  75. print [T, n, yo, y_, m1, m2, Et1, cc, Q1, Et2, cc2, Q2]
  76. else:
  77. print "Malformed input data, please check."
  78. break
  79. result.append([rod_slab, N, T, n, et, cc, es, yo, y_, Q])
  80. return result
  81. def homogenous_config_manual():
  82. print 'Default is ROD geometry'
  83. rod_slab = int(raw_input('Enter 1 if you want to change to SLAB Geometry: '))
  84. if rod_slab == 1:
  85. N = 1
  86. while N % 2 == 1: # N has to be even
  87. N = int(raw_input('Enter the number of Discrete ordinate directions: '))
  88. print 'Solving problem in SLAB Geometry'
  89. print
  90. else:
  91. rod_slab = 0
  92. print 'Solving problem in ROD Geometry'
  93. print
  94. N = 2
  95. T = float(raw_input('Enter the total length of the system: '))
  96. n = 1
  97. while n < 4:
  98. n = int(raw_input('Enter at how many points you want to calculate: '))
  99. et = float(raw_input('Enter the total cross section Sigma_t: '))
  100. cc = 10
  101. while (cc > 1) or (cc < 0):
  102. cc = float(raw_input('Enter the Scattering Ratio c (between 0 and 1): '))
  103. es = cc * et
  104. yo = float(raw_input('Enter the boundary value in the positive direction: '))
  105. y_ = float(raw_input('Enter the boundary value in the negative direction: '))
  106. Q = float(raw_input('Enter the homogeneous isotropic source: '))
  107. return [[rod_slab, N, T, n, et, cc, es, yo, y_, Q]]
  108. def frange(x, y, jump):
  109. result = []
  110. while x <= y:
  111. result.append(x)
  112. x += jump
  113. result.append(y)
  114. return result
  115. def plot(x, SCAL, save, name):
  116. try:
  117. if len(x) > SCAL.shape[0]:
  118. x.pop(-1)
  119. plt.plot(x, SCAL, 'b')
  120. if save:
  121. plt.savefig(name +".png")
  122. else:
  123. plt.show()
  124. plt.clf()
  125. except Exception as e:
  126. print e
  127. print len(x), SCAL.shape
  128. return
  129. # not compatible with virtualenv
  130. def save_csv(matrix, name):
  131. np.savetxt("csv/" + name + ".csv", matrix, delimiter=",")
  132. def homogeneous(var, save=False):
  133. rod_slab = var[0]
  134. N = var[1]
  135. T = var[2]
  136. n = var[3]
  137. et = var[4]
  138. cc = var[5]
  139. es = var[6]
  140. yo = var[7]
  141. y_ = var[8]
  142. Q = var[9]
  143. M = n*N
  144. h = T/n
  145. A = np.zeros((M,M))
  146. B = np.full((1, M), Q/2)
  147. # % GAUSS-LEGENDRE QUADRATURE
  148. beta = range(1,N)
  149. beta = np.divide(beta, np.sqrt(np.subtract(np.multiply(np.power((range(1,N)), 2), 4), 1) ))
  150. #power 2, times 4, minus 1, sqrt each term, element wise divide
  151. #should it return a vector or a value? When testing with random numbers in Matlab, it looks like value. But steps after imply that its a vector
  152. x, w = np.linalg.eig(np.add(np.diag(beta, -1), np.diag(beta, 1))) #check indexing later
  153. u = x #already in the right form
  154. wt = np.multiply(np.power(w[0,:].T, 2), 2)
  155. #first row (index 0), transpose, square, multiply
  156. #central differences
  157. #finite volume
  158. #finite differences
  159. #things to search
  160. if rod_slab != 1:
  161. u[0] = -1
  162. u[1] = 1
  163. # % Diagonal Block of matrix up to N/2
  164. for t in range(1,(N/2) + 1):
  165. s = ((t - 1) * n) - 1 # minus one to account for indexing
  166. ut = u[t-1]
  167. A[s+1, s+1] = (-11 * ut) / (6*h) + et - es * wt[t-1] / 2
  168. A[s+1, s+2] = 3 * ut / h
  169. A[s+1, s+3] = -3 * ut / (2*h)
  170. A[s+1, s+4] = ut/(3*h)
  171. for i in range(2, n):
  172. A[s+i, s+i-1] = -ut/(2*h)
  173. A[s+i, s+i] = (et-es*wt[t-1]/2)
  174. A[s+i, s+i+1] = ut/(2*h)
  175. A[s+n, s+n-1] = -ut / (2*h)
  176. A[s+n, s+n] = (et-es*wt[t-1]/2)
  177. B[0, s+n] = -ut * y_ / (2*h) + Q/2
  178. # % Remaining Blocks in same direction up to N/2
  179. l = t
  180. if (l == 1) and (N > 2):
  181. for p in range(l+1, (N/2) + 1):
  182. S = ((p-1) * n) - 1
  183. for i in range(1, n+1):
  184. A[s+i, S+i] = -es * wt[p-1] / 2
  185. elif (l > 1) and (N > 2):
  186. for p in range(1,l):
  187. S = ((p-1) * n) - 1
  188. for i in range(1, n+1):
  189. A[s+i, S+i] = -es * wt[p-1] / 2
  190. for p in range(l+1, (N/2) + 1):
  191. S = ((p-1) * n) - 1
  192. for i in range(1, n+1):
  193. A[s+i, S+i] = -es * wt[p-1] / 2
  194. # % Blocks from N/2 to N
  195. a = 0
  196. for p in range(1, (N/2) + 1):
  197. S = ((N/2 + p - 1) * n) - 1
  198. for i in range(2, n+1):
  199. A[s+i, S+i-1] = -es * wt[(N/2+p) - 1] / 2
  200. a = a + (es * wt[(N/2+p) - 1] * yo/2)
  201. B[0,s+1] = a + Q/2
  202. # % Diagonal Block of matrix from N/2+1 to N
  203. for t in range(N/2 + 1, N+1):
  204. s = ((t-1) * n) - 1
  205. A[s+1, s+1] = (et-es*wt[t-1]/2)
  206. A[s+1, s+1+1] = u[t-1] / (2 * h)
  207. B[0,s+1] = u[t-1] * yo / (2 * h) + Q/2
  208. for i in range(2, (n-1) + 1):
  209. A[s+i, s+i-1] = -u[t-1] / (2*h)
  210. A[s+i, s+i] = (et-es*wt[t-1]/2)
  211. A[s+i, s+i+1] = u[t-1]/(2*h)
  212. A[s+n, s+n-3] = -u[t-1] / (3*h)
  213. A[s+n, s+n-2] = 3 * u[t-1] / (2*h)
  214. A[s+n, s+n-1] = -3 * u[t-1] / h;
  215. A[s+n, s+n] = (11 * u[t-1] / (6*h) + et - es * wt[t-1] / 2)
  216. # % Remaining Blocks in same direction up to N
  217. l = t
  218. if (l==N/2 + 1) and (N>2):
  219. for p in range(l+1, N+1):
  220. S = ((p-1) * n) - 1
  221. for i in range(1, n+1):
  222. A[s+i, S+i] = -es * wt[p-1] / 2
  223. elif (1 > N/2 + 1) and (N>2):
  224. for p in range(N/2 + 1, (l-1) + 1):
  225. S = ((p-1) * n) - 1
  226. for i in range(1, n+1):
  227. A[s+i, S+i] = -es * wt[p-1] / 2
  228. for p in range(l+1, N+1):
  229. S = ((p-1) * n) - 1
  230. for i in range(1, n+1):
  231. A[s+i, S+i] = -es * wt[p-1] / 2
  232. # % Blocks from 1 to N/2
  233. a = 0
  234. for p in range(1, (N/2) + 1):
  235. S = ((p-1) * n) - 1
  236. for i in range(1, (n-1) + 1):
  237. A[s+i, S+i+1] = -es * wt[p-1] / 2
  238. a = a + (es * wt[p-1] * y_ / 2)
  239. B[0,s+n] = a + Q / 2
  240. save_csv(A, "A")
  241. save_csv(B, "B")
  242. #may be some accuracy solving issues algorithm wise with numpy
  243. # refer to http://stackoverflow.com/questions/25001753/numpys-linalg-solve-and-linalg-lstsq-not-giving-same-answer-as-matlabs
  244. # below are different attempts to see if any difference RESULT: no difference
  245. # AX = B.T
  246. # A1 = np.linalg.inv(A)
  247. # np.savetxt("Ainv.csv", A1, delimiter=",")
  248. # X = np.dot(A1, B.T)
  249. # np.savetxt("B.csv", B.T, delimiter=",")
  250. # np.savetxt("X_manual.csv", X, delimiter=",")
  251. # np.savetxt("X_solve.csv", X1, delimiter=",")
  252. # X2 = np.linalg.lstsq(A, B.T)[0]
  253. # np.savetxt("X_lstsq.csv", X2, delimiter=",")
  254. # print A1 == A
  255. # print np.linalg.cond(A1)
  256. X = np.linalg.solve(A, B.T) #checked values compared to the other methods of solving and it seems to be accurate enough
  257. Y = np.zeros((M+N,1))
  258. # % Adding boundary conditions to the array
  259. i = 1
  260. for t in range(1, (N/2) + 1):
  261. s = (t-1)
  262. for j in range(i, (t*n+s) + 1): #might have conflict here with the s-1
  263. Y[j-1] = X[(j-s)-1]
  264. Y[(j+1) - 1] = y_
  265. i = j + 2
  266. i -= 1
  267. for t in range((N/2) + 1, N+1):
  268. Y[(i+1) -1] = yo
  269. i += 1
  270. for j in range(i+1, (t*n+t) + 1):
  271. Y[j - 1] = X[(j-t) - 1]
  272. i = j
  273. RL = 0
  274. for t in range(1, (N/2) + 1):
  275. s = (t-1) * n
  276. RL = RL + abs(wt[t-1] * u[t-1] * Y[(s+t) - 1]) # might have conflict here
  277. TR = 0
  278. for t in range((N/2) + 1, N + 1):
  279. s = (t-1) * n
  280. S = (t - N/2 - 1) * n
  281. k = N/2 + 1
  282. TR = TR + wt[t-1] * u[t-1] * Y[(s+n+t) - 1] # might have a conflict here
  283. m = n + 1
  284. SCAL = np.zeros((m, 1))
  285. for t in range(1, N + 1):
  286. s = (t-1) * m
  287. for i in range(1, m + 1):
  288. SCAL[i-1] = SCAL[i-1] + wt[t-1] * Y[(s+i) - 1]
  289. save_csv(Y, "Y")
  290. save_csv(SCAL, "SCAL")
  291. x = frange(0, T, h)
  292. #forming name for the file save
  293. name = str(rod_slab) + "_" + str(T) + "_" + str(N) + "_" + str(n) + "_" + str(et) + "_" + str(cc) + "_" + str(es) + "_" + str(yo) + "_" + str(y_) + "_" + str(Q)
  294. plot(x, SCAL, save, name)
  295. #extra order does matter
  296. #check h ascending order
  297. def SN_per_bench_solver(T,m1,m2,n,N,Es1,Es2,Et1,Et2,yo,y_,Q1,Q2,u,wt,a):
  298. interval = T / n
  299. m12 = [0, m1, m2] #which material, so index of this matters
  300. mm = 0
  301. s = 0
  302. i = 0
  303. x = None
  304. if a > 1:
  305. if a < m2/interval + 2:
  306. i += 1
  307. x1 = (a-1) * interval
  308. s += x1
  309. tmp = [] #can shorten to declare immediately but verbose for readability
  310. tmp.append(s)
  311. tmp.append(2)
  312. x = np.asarray([tmp]) #x is offset to account for material
  313. else:
  314. i += 1
  315. x1 = (a - m2/interval - 1) * interval
  316. s += x1;
  317. tmp = []
  318. tmp.append(s)
  319. tmp.append((mm % 2) + 1)
  320. x = np.asarray([tmp])
  321. mm += 1
  322. while s < T:
  323. i += 1
  324. x1 = m12[(mm % 2) + 1] #whats the point of this? isnt this either 1 or 2?
  325. s = s + x1
  326. tmp = []
  327. if s <= T:
  328. tmp.append(s)
  329. else:
  330. tmp.append(T) #check here
  331. tmp.append((mm % 2) + 1)
  332. if x == None:
  333. x = np.asarray([tmp])
  334. else: #check logic here
  335. x = np.vstack((x, np.asarray([tmp])))
  336. mm += 1
  337. H = T / n
  338. n1 = 1
  339. i = 1
  340. j = 1
  341. extra = [0] * n
  342. t1 = i * H
  343. save_csv(x, "xsn")
  344. tmp = [0, int(x[j-1 , 2-1])]
  345. L = np.asarray([tmp])
  346. L1 = []
  347. L1.append(tmp)
  348. h = [] #is h just a list?
  349. if t1 == x[j-1, 1-1]:
  350. extra[i-1] = 1
  351. h.append(x[j-1, 1-1] / 2)
  352. #n1+1 case
  353. tmp = [h[n1-1], int(x[j-1, 2-1])]
  354. L1.append(tmp)
  355. tmp = np.asarray([tmp])
  356. L = np.vstack((L, tmp))
  357. h.append(h[n1-1])
  358. #n1+2 case
  359. tmp = np.asarray([x[j-1,1-1], 3]) #are the values of x guaranteed to be ints???
  360. L1.append([x[j-1,1-1], 3])
  361. L = np.vstack((L, tmp))
  362. n1 += 2
  363. i += 1
  364. t1 = i * H
  365. else:
  366. while t1 <= x[j-1,1-1]:
  367. h.append(H) #h(n1) = H
  368. tmp = [i * H]
  369. if tmp[0] == x[j-1,1-1]:
  370. tmp.append(3)
  371. else:
  372. tmp.append(int(x[j-1,2-1]))
  373. L1.append(tmp)
  374. L = np.vstack((L, np.asarray([tmp])))
  375. n1 += 1
  376. i += 1
  377. t1 = i * H
  378. j = 2
  379. if x[1-1,1-1] != T:
  380. while i <= n:
  381. while t1 <= x[j-1,1-1]:
  382. tmp = [i * H]
  383. if tmp[0] == x[j-1,1-1]:
  384. tmp.append(3)
  385. else:
  386. tmp.append(int(x[j-1,2-1]))
  387. L1.append(tmp)
  388. L = np.vstack((L, np.asarray([tmp])))
  389. h.append(L1[n1+1-1][1-1] - L1[n1-1][1-1])
  390. n1 += 1 #n1 isnt really necessary to keep the index
  391. i += 1
  392. t1 = i * H
  393. j += 1
  394. if L1[n1-1][1-1] == T and L1[n1-1][2-1] == 3 and L1[n1-1-1][2-1] == 3:
  395. i -= 1
  396. extra[i-1] += 1
  397. tmp = []
  398. tmp.append((x[j-1-1, 1-1] + x[j-2-1, 1-1])/2)
  399. tmp.append(x[j-1-1, 2-1])
  400. tmp[1] = int(tmp[1])
  401. L = np.vstack((L, np.asarray([tmp])))
  402. L1[n1-1] = tmp
  403. h[n1-1-1] = L1[n1-1][1-1] - L1[n1-1-1][1-1] #still within range
  404. n1 += 1
  405. tmp = [x[j-1-1, 1-1], 3]
  406. tmp[1] = int(tmp[1])
  407. L1.append(tmp)
  408. L = np.vstack((L, np.asarray([tmp])))
  409. h.append(L1[n1-1][1-1] - L1[n1-1-1][1-1]) #adds to h[] (Matlab) not indexing
  410. i += 1
  411. save_csv(L1, "L")
  412. n1 -= 1
  413. L = L1
  414. save_csv(np.asarray(h), "h")
  415. #L tells you which material, L's index doesnt matter, x's index doesnt matter, Es, Et, Q does
  416. A = np.zeros((N*n1,N*n1))
  417. B = np.zeros((N*n1,1))
  418. Et12 = [0, Et1, Et2]
  419. Es12 = [0, Es1, Es2]
  420. Q12 = [0, Q1, Q2]
  421. #% Diagonal Block of matrix up to N/2.................................
  422. for t in range(1, N/2 + 1):
  423. s = (t-1) * n1
  424. A[s+1-1,s+1-1] = -u[t-1] * (1 / h[1-1]) + Et12[L[1-1][2-1]] - Es12[L[1-1][2-1]] * wt[t-1]/2
  425. A[s+1-1,s+2-1] = u[t-1] * (1 / h[1-1])
  426. B[s+1-1, 0] = Q12[L[1-1][2-1]]/2
  427. for i in range(2, n1-1+1):
  428. if L[i-1][2-1] == 3:
  429. if i == n1 - 1:
  430. A[s+i-1,s+i-1] = -u[t-1] * (1 / h[i-1] + 1 / (h[i-1] + h[i+1-1])) + Et12[L[i+1-1][2-1]] - Es12[L[i+1-1][2-1]] * wt[t-1]/2
  431. A[s+i-1,s+i+1-1] = u[t-1] * (1/h[i-1] + 1/h[i+1-1])
  432. B[s+i-1,0] = u[t-1] * y_ * (h[i-1] / (h[i+1-1] * (h[i-1] + h[i+1-1]))) + Q12[L[i+1-1][2-1]]/2
  433. else:
  434. A[s+i-1,s+i-1] = -u[t-1] * (1/h[i-1] + 1/(h[i-1] + h[i+1-1] )) + Et12[L[i+1-1][2-1]] - Es12[L[i+1-1][2-1]] * wt[t-1]/2
  435. A[s+i-1,s+i+1-1] = u[t-1] * (1/h[i-1] + 1/h[i+1-1])
  436. A[s+i-1,s+i+2-1] = -u[t-1] * (h[i-1]/(h[i+1-1] * (h[i-1] + h[i+1-1])))
  437. B[s+i-1, 0] = Q12[L[i+1-1][2-1]]/2
  438. else:
  439. A[s+i-1,s+i-1-1] = -u[t-1] * h[i-1]/(h[i-1-1]*(h[i-1-1]+h[i-1]))
  440. A[s+i-1,s+i-1] = u[t-1] * (h[i-1] - h[i-1-1])/(h[i-1] * h[i-1-1]) + Et12[L[i-1][2-1]] - Es12[L[i-1][2-1]] * wt[t-1]/2
  441. A[s+i-1,s+i+1-1] = u[t-1] * h[i-1-1] / (h[i-1] * (h[i-1-1] + h[i-1]))
  442. B[s+i-1,0] = Q12[L[i-1][2-1]] /2
  443. A[s+n1-1,s+n1-1-1] = -u[t-1] * h[n1-1] / (h[n1-1-1] * (h[n1-1-1] + h[n1-1]))
  444. A[s+n1-1,s+n1-1] = u[t-1] * (h[n1-1] - h[n1-1-1]) / (h[n1-1] * h[n1-1-1]) + Et12[L[n1-1][2-1]] - Es12[L[n1-1][2-1]] * wt[t-1]/2
  445. B[s+n1-1,0] = -u[t-1] * y_ * h[n1-1-1] / (h[n1-1] * (h[n1-1-1] + h[n1-1] )) + Q12[L[n1-1][2-1]]/2
  446. #% Blocks from N/2 to N........................................
  447. a = 0
  448. for p in range(1, N/2 + 1):
  449. S = (N/2 + p - 1) * n1
  450. for i in range(2, n1+1):
  451. if L[i-1][2-1] == 3:
  452. A[s+i-1,S+i-1-1] = -Es12[L[i+1-1][2-1]] * wt[N/2+p-1]/2
  453. else:
  454. A[s+i-1,S+i-1-1] = -Es12[L[i-1][2-1]] * wt[N/2+p -1]/2
  455. a += (Es12[L[1-1][2-1]] * wt[N/2+p-1] * yo/2)
  456. B[s+1-1,0] = B[s+1-1,0] + a
  457. #% Diagonal Block of matrix from N/2+1 to N.........................
  458. for t in range(N/2+1, N+1):
  459. s = (t-1) * n1
  460. A[s+1-1,s+1-1] = u[t-1] * (h[2-1] - h[1-1]) / (h[2-1] * h[1-1]) + Et12[L[1-1][2-1]] - Es12[L[1-1][2-1]] * wt[t-1]/2
  461. A[s+1-1,s+1+1-1] = u[t-1] * h[1-1] / (h[2-1] * (h[1-1] + h[2-1]))
  462. B[s+1-1,0] = u[t-1] * yo * h[2-1] / (h[1-1] * (h[1-1] + h[2-1])) + Q12[L[1-1][2-1]] / 2
  463. for i in range(2, n1-1+1):
  464. if L[i+1-1][2-1] == 3:
  465. if i == 2:
  466. A[s+i-1,s+i-1] = u[t-1] * (1/h[i-1] + 1/(h[i-1] + h[i-1-1])) + Et12[L[i-1][2-1]] - Es12[L[i-1][2-1]] * wt[t-1]/2
  467. A[s+i-1,s+i-1-1] = -u[t-1] * (1/h[i-1] + 1/h[i-1-1])
  468. B[s+i-1,0] = -u[t-1] * yo * (h[i-1]/(h[i-1-1] * (h[i-1] + h[i-1-1]))) + Q12[L[i-1][2-1]]/2
  469. else:
  470. A[s+i-1,s+i-1] = u[t-1] * (1/h[i-1] + 1/(h[i-1] + h[i-1-1])) + Et12[L[i-1][2-1]] - Es12[L[i-1][2-1]] * wt[t-1]/2
  471. A[s+i-1,s+i-1-1] = -u[t-1] * (1/h[i-1] + 1/h[i-1-1])
  472. A[s+i-1,s+i-2-1] = u[t-1] * (h[i-1] / (h[i-1-1] * (h[i-1] + h[i-1-1])))
  473. B[s+i-1,0] = Q12[L[i-1][2-1]]/2
  474. else:
  475. A[s+i-1,s+i-1-1] = -u[t-1] * h[i+1-1] / (h[i-1] * (h[i-1] + h[i+1-1]))
  476. A[s+i-1,s+i-1] = u[t-1] * (h[i+1-1] - h[i-1]) / (h[i+1-1] * h[i-1]) + Et12[L[i+1-1][2-1]] - Es12[L[i+1-1][2-1]] * wt[t-1]/2
  477. A[s+i-1,s+i+1-1] = u[t-1] * h[i-1] / (h[i+1-1] * (h[i-1] + h[i+1-1]))
  478. B[s+i-1,0] = Q12[L[i+1-1][2-1]]/2
  479. A[s+n1-1,s+n1-1] = u[t-1] * (1/h[n1-1]) + Et12[L[n1-1][2-1]] - Es12[L[n1-1][2-1]] * wt[t-1]/2
  480. A[s+n1-1,s+n1-1-1] = -u[t-1] * (1/h[n1-1])
  481. B[s+n1-1,0] = Q12[L[n1-1][2-1]]/2
  482. #% Blocks from 1 to N/2........................................
  483. a = 0
  484. for p in range(1, N/2+1):
  485. S = (p-1) * n1
  486. for i in range(1,n1-1+1):
  487. if L[i+1-1][2-1] == 3:
  488. A[s+i-1,S+i+1-1] = -Es12[L[i-1][2-1]] * wt[p-1]/2
  489. else:
  490. A[s+i-1,S+i+1-1] = -Es12[L[i+1-1][2-1]] * wt[p-1]/2
  491. a += (Es12[L[n1-1][2-1]] * wt[p-1] * y_/2)
  492. B[s+n1-1, 0] += a
  493. Z = np.linalg.solve(A, B)
  494. return Z, n1, B, L, A, extra
  495. def periodic_config_manual(): #redundant, clean up to make it generic later
  496. print 'Default is ROD geometry'
  497. rod_slab = int(raw_input('Enter 1 if you want to change to SLAB Geometry: '))
  498. if rod_slab == 1:
  499. N = 1
  500. while N % 2 == 1: # N has to be even
  501. N = int(raw_input('Enter the number of Discrete ordinate directions: '))
  502. print 'Solving problem in SLAB Geometry'
  503. print
  504. else:
  505. rod_slab = 0
  506. print 'Solving problem in ROD Geometry'
  507. print
  508. N = 2
  509. T = float(raw_input('Enter the total length of the system: '))
  510. n = 1
  511. while n < 4:
  512. n = int(raw_input('Enter at how many points you want to calculate: '))
  513. et = float(raw_input('Enter the total cross section Sigma_t: '))
  514. cc = 10
  515. while (cc > 1) or (cc < 0):
  516. cc = float(raw_input('Enter the Scattering Ratio c (between 0 and 1): '))
  517. es = cc * et
  518. yo = float(raw_input('Enter the boundary value in the positive direction: '))
  519. y_ = float(raw_input('Enter the boundary value in the negative direction: '))
  520. Q = float(raw_input('Enter the homogeneous isotropic source: '))
  521. return [[rod_slab, N, T, n, et, cc, es, yo, y_, Q]]
  522. def periodic(var):
  523. #T and n have to be floats for math
  524. rod_slab = var[0]
  525. N = var[1] #has to be even add checking functions to config file function later
  526. T = var[2]
  527. n = var[3]
  528. yo = var[4]
  529. y_ = var[5]
  530. m1 = var[6]
  531. m2 = var[7]
  532. Et1 = var[8]
  533. cc = var[9]
  534. Q1 = var[10]
  535. Et2 = var[11]
  536. cc2 = var[12]
  537. Q2 = var[13]
  538. beta = range(1,N)
  539. beta = np.divide(beta, np.sqrt(np.subtract(np.multiply(np.power((range(1,N)), 2), 4), 1) ))
  540. x, w = np.linalg.eig(np.add(np.diag(beta, -1), np.diag(beta, 1))) #check indexing later
  541. u = x #already in the right form
  542. wt = np.multiply(np.power(w[0,:].T, 2), 2)
  543. if rod_slab != 1:
  544. u[0] = -1
  545. u[1] = 1
  546. #ROD case
  547. # N = 2
  548. # wt = [1, 1]
  549. # u = [-1, 1]
  550. #Add ROD SLAB code to deal with weights, weights is 1xn mat
  551. Es1 = cc * Et1
  552. Ea1 = Et1 - Es1
  553. Es2 = cc2 * Et2
  554. Ea2 = Et2 - Es2
  555. a = 1
  556. reflec = 0
  557. reflec2 = 0
  558. transm = 0
  559. transm2 = 0
  560. SF = np.zeros((n+1,1))
  561. SF2 = np.zeros((n+1,1))
  562. cond = 0
  563. total = (m1+m2) / (T/n)
  564. while (a < (total+1)):
  565. print "Problem " + str(a) + " of " + str(total)
  566. Z,n1,B,L,A, extra = SN_per_bench_solver(T,m1,m2,n,N,Es1,Es2,Et1,Et2,yo,y_,Q1,Q2,u,wt,a)
  567. save_csv(np.asarray(extra), "extra")
  568. save_csv(L, "L")
  569. save_csv(A, "A")
  570. save_csv(B, "B")
  571. save_csv(Z, "Z")
  572. #% Adjusting points..........................
  573. X = np.zeros((n*N,1))
  574. for i in range(1, (N/2) + 1):
  575. X[(i-1)*n+1 - 1] = Z[(i-1)*n1+1 - 1]
  576. k = 2 #not sure why this is needed.
  577. for j in range(2, n + 1):
  578. k += extra[j-1-1]
  579. # k is correct
  580. #there is some trickery going on with k.
  581. X[(i-1)*n+j - 1] = Z[(i-1)*n1+k - 1]
  582. k += 1
  583. for i in range(N/2 + 1, N + 1):
  584. k = 1
  585. for j in range(1, n + 1):
  586. k += extra[j-1]
  587. X[(i-1)*n+j - 1] = Z[(i-1)*n1+k - 1]
  588. k += 1
  589. Y = np.zeros((N*(n+1),1));
  590. #% Adding boundary conditions...............
  591. i = 1
  592. for t in range(1, N/2 + 1):
  593. s = t-1
  594. for j in range(i, t*n+s + 1):
  595. Y[j - 1] = X[j-s - 1]
  596. Y[j+1 - 1] = y_
  597. i = j + 2
  598. i -= 1
  599. for t in range(N/2+1, N+1):
  600. Y[i+1 - 1] = yo
  601. i += 1
  602. for j in range(i+1, t*n+t + 1):
  603. Y[j - 1] = X[j-t - 1]
  604. i = j
  605. save_csv(X, "X")
  606. save_csv(Y, "Y")
  607. #% Calculating Reflection, Transmission and Scalar Flux......
  608. RL = 0
  609. for t in range(1, N/2 + 1):
  610. s = (t-1) * n
  611. RL += abs(wt[t - 1] * u[t - 1] * Y[s+t - 1])
  612. TR = 0
  613. for t in range(N/2 + 1, N + 1):
  614. s = (t-1) * n
  615. TR += wt[t - 1] * u[t - 1] * Y[s+n+t - 1]
  616. m = n + 1
  617. SCAL = np.zeros((m,1))
  618. for t in range(1, N + 1):
  619. s = (t-1) * m
  620. for i in range(1, m + 1):
  621. SCAL[i - 1] = SCAL[i - 1] + wt[t -1] * Y[s+i - 1]
  622. reflec += RL
  623. reflec2 += (RL)**2
  624. transm += TR
  625. transm2 += (TR)**2
  626. for i in range(1, n+1 + 1):
  627. SF[i - 1] = SF[i - 1] + SCAL[i - 1]
  628. SF2[i - 1] = SF2[i - 1] + (SCAL[i - 1])**2
  629. a += 1
  630. a -= 1
  631. SF = np.divide(SF, a)
  632. SF2 = np.divide(SF2, a)
  633. kk = T/n
  634. #% p is line along x-axis.
  635. p = frange(0, T, kk)
  636. p.pop()
  637. save_csv(SF, "SF")
  638. plt.plot(p,SF,'r')
  639. plt.plot(p,SF2,'b')
  640. plt.show()
  641. def non_homogenous(var):
  642. rod_slab = var[0]
  643. N = var[1]
  644. T = var[2]
  645. m1 = var[3]
  646. m2 = var[4]
  647. n = var[5]
  648. N = var[6]
  649. Es1 = var[7]
  650. Es2 = var[8]
  651. Et1 = var[9]
  652. Et2 = var[10]
  653. yo = var[11]
  654. y_ = var[12]
  655. Q1 = var[13]
  656. Q2 = var[14]
  657. u = var[15]
  658. wt = var[16]
  659. randseed = var[17]
  660. med = var[18]
  661. a = var[19]
  662. # %Building realization----------
  663. # %EXP. Random Medium....
  664. #random python docs: https://docs.python.org/2/library/random.html
  665. if not med:
  666. random.seed(randseed)
  667. m12 = [0, m1, m2] #refers to which material, pad with 0 so index 1 == material 1 (slightly less confusing this way)
  668. xx == 0.0
  669. xx = random.random() #python random gives a float in [0.0, 1.0)
  670. #turns out both Matlab and python use deterministic randomgenerators
  671. if xx <= m1/(m1+m2):
  672. mm = 0
  673. else:
  674. mm = 1
  675. s = 0
  676. i = 0
  677. x = None #is this a guaranteed loop?
  678. while s < T:
  679. i += 1
  680. x1 = random.expovariate(1/m12[mm%2 + 1]) #exprnd returns a 1x1
  681. #why does x1 matter? its not being used in the loop unless it changes the previous value for random need to check matlab docs
  682. s += x1
  683. tmp = []
  684. if s <= T:
  685. tmp.append(s)
  686. else:
  687. tmp.append(T)
  688. tmp.append(mm%2 + 1)
  689. if not x:
  690. x = np.asarray([tmp])
  691. else:
  692. x = np.vstack((x, np.asarray([tmp])))
  693. mm += 1
  694. randseed = random.getstate() #use setstate() later
  695. else:
  696. interval = T/n
  697. m12 = [0, m1, m2]
  698. mm = 9
  699. i = 0
  700. s = 0
  701. if a > 1:
  702. if a < m1/interval+2: #check order of op here
  703. i += 1
  704. x1 = (a-1) * interval
  705. s += x1
  706. x = np.asarray([[s, 2]])
  707. else: # not sure if this is inner or outer if CHECK
  708. i += 1
  709. x1 = (a - m1/interval - 1) * interval
  710. s += x1
  711. x = np.asarray([s, mm%2 + 1])
  712. mm += 1
  713. while s < T:
  714. i += 1
  715. x1 = m12[mm%2 + 1]
  716. s += 1
  717. tmp = []
  718. if s <= T:
  719. tmp.append(s)
  720. else:
  721. tmp.append(T)
  722. tmp.append(mm%2 + 1)
  723. x = np.vstack((x, np.asarray([tmp])))
  724. mm += 1
  725. # %Adding interfaces and extra points inside layers--------------
  726. H = T/n
  727. n1 = 1
  728. i = 1
  729. j = 1
  730. extra = [] #using a list will be easier for extra follow periodic code
  731. h1 = [] # ditto comment above
  732. t1 = i * H
  733. L = np.asarray([[0, x[j-1, 2-1]]])
  734. if t1 > x[j-1, 1-1]:
  735. extra[i-1] = 2
  736. h[n1-1] = x[j-1, 1-1]/2
  737. L[n1+1-1,1-1] = h[n1-1]
  738. L[n1+1-1,2-1] = x[j-1, 2-1]
  739. h[n1+1-1]=h[n1-1]
  740. L[n1+2-1,1-1] = x[j-1, 1-1]
  741. L[n1+2-1,2-1] = 3
  742. n1 += 2
  743. else:
  744. while t1 <= x[j-1, 1-1]:
  745. h[n1-1] = H
  746. def main(sysargs):
  747. # default state is file configuration
  748. # if not len(sysargs):
  749. # save = False
  750. # elif sysargs[0].lower() == "manual":
  751. # var = homogenous_config_manual()
  752. # elif sysargs[0].lower() == "save":
  753. # save = True
  754. # var = config_file()
  755. # for v in var:
  756. # homogenous(v, save)
  757. #lots of redudant code, need to make them generic later
  758. periodic([0, 2, 20.0, 100, 0.0, 0.0, 1.0, 1.0, 1.0, 0.5, 1.0, 0.0, 0.0, 0.0])
  759. if __name__ == '__main__':
  760. main(sys.argv[1:])