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

New_SN_bench_solver_MC.m 10KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382
  1. %THIS IS THE AUXILIARY FUNCTION OF THE MAIN CODE "New_SN_bench"
  2. %RICHARD VASQUES & NITIN KUMAR YADAV
  3. function[Z,extra,n1,randseed] = New_SN_bench_solver_MC(T,m1,m2,n,N,Es1,Es2,Et1,Et2,yo,y_,Q1,Q2,u,wt,randseed,med,a)
  4. %Building realization----------
  5. %EXP. Random Medium....
  6. if med==0
  7. rng(randseed);
  8. m12=zeros(2,1);
  9. m12(1)=m1;m12(2)=m2;
  10. xx=rand(1);
  11. if xx<=m1/(m1+m2)
  12. mm=0;
  13. else
  14. mm=1;
  15. end
  16. s=0;i=0;
  17. while s<T
  18. i=i+1;
  19. x1=exprnd(m12(mod(mm,2)+1),1,1);
  20. s=s+x1;
  21. if s<=T
  22. x(i,1)=s;
  23. else
  24. x(i,1)=T;
  25. end
  26. x(i,2)=mod(mm,2)+1;
  27. mm=mm+1;
  28. end
  29. randseed=rng;
  30. else
  31. interval=T/n;
  32. %Periodic Medium........
  33. m12=zeros(2,1);
  34. m12(1)=m1;m12(2)=m2;
  35. mm=0;
  36. i=0;
  37. s=0;
  38. if a>1
  39. if a<m1/interval+2
  40. i=i+1;
  41. x1=(a-1)*interval;
  42. s=s+x1;
  43. x(i,1)=s;
  44. x(i,2)=2;
  45. else
  46. i=i+1;
  47. x1=(a-m1/interval-1)*interval;
  48. s=s+x1;
  49. x(i,1)=s;
  50. x(i,2)=mod(mm,2)+1;
  51. mm=mm+1;
  52. end
  53. end
  54. while s<T
  55. i=i+1;
  56. x1=m12(mod(mm,2)+1);
  57. s=s+x1;
  58. if s<=T
  59. x(i,1)=s;
  60. else
  61. x(i,1)=T;
  62. end
  63. x(i,2)=mod(mm,2)+1;
  64. mm=mm+1;
  65. end
  66. end
  67. %Adding interfaces and extra points inside layers--------------
  68. %L(*,1) is the spatial point
  69. %L(*,2) indicates in which material is the spatial point (1 or 2).
  70. %If L(*,2)=3, the point is an interface.
  71. H=T/n;n1=1;i=1;j=1;
  72. extra=zeros(n,1);
  73. t1=i*H;
  74. L(1,1)=0;L(1,2)=x(j,2);
  75. if t1>x(j,1)
  76. extra(i)=2;
  77. h(n1)=x(j,1)/2;
  78. L(n1+1,1)=h(n1);
  79. L(n1+1,2)=x(j,2);
  80. h(n1+1)=h(n1);
  81. L(n1+2,1)=x(j,1);
  82. L(n1+2,2)=3;
  83. n1=n1+2;
  84. else
  85. while t1<=x(j,1)
  86. h(n1)=H;
  87. L(n1+1,1)=i*H;
  88. if L(n1+1,1)==x(j,1) % checking is point is interface point?
  89. L(n1+1,2)=3;
  90. else
  91. L(n1+1,2)=x(j,2);
  92. end
  93. n1=n1+1;
  94. i=i+1;
  95. t1=i*H;
  96. end
  97. if x(j,1)<T && L(n1)~=x(j,1)
  98. extra(i)=1;
  99. L(n1+1,1)=x(j,1);
  100. L(n1+1,2)=3;
  101. h(n1)=L(n1+1,1)-L(n1,1);
  102. n1=n1+1;
  103. end
  104. if t1==T && (x(j+1,1)-x(j,1))<H
  105. extra(i)=extra(i)+1;
  106. j=j+1;
  107. L(n1+1,1)=(x(j,1)+x(j-1,1))/2;
  108. L(n1+1,2)=x(j-1,2);
  109. h(n1)=L(n1+1,1)-L(n1,1);
  110. n1=n1+1;
  111. L(n1+1,1)=x(j,1);
  112. L(n1+1,2)=3;
  113. h(n1)=L(n1+1,1)-L(n1,1);
  114. n1=n1+1;
  115. i=i+1;
  116. end
  117. end
  118. j=2;
  119. if x(1,1)~=T
  120. while i<=n
  121. if t1>x(j,1)
  122. extra(i)=extra(i)+2;
  123. L(n1+1,1)=(x(j,1)+x(j-1,1))/2;
  124. L(n1+1,2)=x(j,2);
  125. h(n1)=L(n1+1,1)-L(n1,1);
  126. n1=n1+1;
  127. L(n1+1,1)=x(j,1);
  128. L(n1+1,2)=3;
  129. h(n1)=L(n1+1,1)-L(n1,1);
  130. n1=n1+1;
  131. else
  132. while t1<=x(j,1)
  133. L(n1+1,1)=i*H;
  134. if L(n1+1,1)==x(j,1) %checking is point is interface point?
  135. L(n1+1,2)=3;
  136. else
  137. L(n1+1,2)=x(j,2);
  138. end
  139. h(n1)=L(n1+1,1)-L(n1,1);
  140. n1=n1+1;
  141. i=i+1;
  142. t1=i*H;
  143. end
  144. if x(j,1)<T && L(n1)~=x(j,1)
  145. extra(i)=extra(i)+1;
  146. L(n1+1,1)=x(j,1);
  147. L(n1+1,2)=3;
  148. h(n1)=L(n1+1,1)-L(n1,1);
  149. n1=n1+1;
  150. end
  151. end
  152. j=j+1;
  153. if t1==T && (x(j,1)-x(j-1,1))<H
  154. extra(i)=extra(i)+1;
  155. L(n1+1,1)=(x(j,1)+x(j-1,1))/2;
  156. L(n1+1,2)=x(j-1,2);
  157. h(n1)=L(n1+1,1)-L(n1,1);
  158. n1=n1+1;
  159. L(n1+1,1)=x(j,1);
  160. L(n1+1,2)=3;
  161. h(n1)=L(n1+1,1)-L(n1,1);
  162. n1=n1+1;
  163. i=i+1;
  164. end
  165. end
  166. end
  167. n1=n1-1;
  168. A=zeros(N*n1,N*n1); B=zeros(N*n1,1);
  169. Et12=zeros(2,1); Es12=zeros(2,1); Q12=zeros(2,1);
  170. Et12(1)=Et1;Et12(2)=Et2;
  171. Es12(1)=Es1;Es12(2)=Es2;
  172. Q12(1)=Q1;Q12(2)=Q2;
  173. minleft=T/2-m1/2;
  174. maxright=T/2+m1/2;
  175. % Diagonal Block of matrix up to N/2.................................
  176. for t=1:N/2
  177. s=(t-1)*n1;
  178. A(s+1,s+1)=-u(t)*(1/h(1)+1/(h(1)+h(2)))+Et12(L(1,2))-Es12(L(1,2))*wt(t)/2;
  179. A(s+1,s+2)=u(t)*(1/h(1)+1/h(2));
  180. A(s+1,s+3)=-u(t)*(h(1)/(h(2)*(h(1)+h(2))));
  181. %B(s+1)=Q12(L(1,2))/2;
  182. B(s+1)=0.;
  183. for i=2:n1-1
  184. if L(i,2)==3
  185. if i==n1-1
  186. A(s+i,s+i)=-u(t)*(1/h(i)+1/(h(i)+h(i+1)))+Et12(L(i+1,2))-Es12(L(i+1,2))*wt(t)/2;
  187. A(s+i,s+i+1)=u(t)*(1/h(i)+1/h(i+1));
  188. if L(i+1,1)>=minleft && L(i+1,1)<=maxright
  189. B(s+i)=u(t)*y_*(h(i)/(h(i+1)*(h(i)+h(i+1))))+Q12(L(i+1,2))/2;
  190. else
  191. B(s+i)=u(t)*y_*(h(i)/(h(i+1)*(h(i)+h(i+1))));
  192. end
  193. else
  194. A(s+i,s+i)=-u(t)*(1/h(i)+1/(h(i)+h(i+1)))+Et12(L(i+1,2))-Es12(L(i+1,2))*wt(t)/2;
  195. A(s+i,s+i+1)=u(t)*(1/h(i)+1/h(i+1));
  196. A(s+i,s+i+2)=-u(t)*(h(i)/(h(i+1)*(h(i)+h(i+1))));
  197. if L(i+1,1)>=minleft && L(i+1,1)<=maxright
  198. B(s+i)=Q12(L(i+1,2))/2;
  199. else
  200. B(s+i)=0.;
  201. end
  202. end
  203. else
  204. A(s+i,s+i-1)=-u(t)*h(i)/(h(i-1)*(h(i-1)+h(i)));
  205. A(s+i,s+i)=u(t)*(h(i)-h(i-1))/(h(i)*h(i-1))+Et12(L(i,2))-Es12(L(i,2))*wt(t)/2;
  206. A(s+i,s+i+1)=u(t)*h(i-1)/(h(i)*(h(i-1)+h(i)));
  207. if L(i,1)>=minleft && L(i,1)<=maxright
  208. B(s+i)=Q12(L(i,2))/2;
  209. else
  210. B(s+i)=0.;
  211. end
  212. end
  213. end
  214. A(s+n1,s+n1-1)=-u(t)*h(n1)/(h(n1-1)*(h(n1-1)+h(n1)));
  215. A(s+n1,s+n1)=u(t)*(h(n1)-h(n1-1))/(h(n1)*h(n1-1))+Et12(L(n1,2))-Es12(L(n1,2))*wt(t)/2;
  216. %B(s+n1)=-u(t)*y_*h(n1-1)/(h(n1)*(h(n1-1)+h(n1)))+Q12(L(n1,2))/2;
  217. B(s+n1)=-u(t)*y_*h(n1-1)/(h(n1)*(h(n1-1)+h(n1)));
  218. % Remaining Blocks in same direction up to N/2..............
  219. l=t;
  220. if l==1 && N>2
  221. for p=l+1:N/2
  222. S=(p-1)*n1;
  223. for i=1:n1
  224. if L(i,2)==3
  225. A(s+i,S+i)=-Es12(L(i+1,2))*wt(p)/2;
  226. else
  227. A(s+i,S+i)=-Es12(L(i,2))*wt(p)/2;
  228. end
  229. end
  230. end
  231. elseif l>1 && N>2
  232. for p=1:l-1
  233. S=(p-1)*n1;
  234. for i=1:n1
  235. if L(i,2)==3
  236. A(s+i,S+i)=-Es12(L(i+1,2))*wt(p)/2;
  237. else
  238. A(s+i,S+i)=-Es12(L(i,2))*wt(p)/2;
  239. end
  240. end
  241. end
  242. for p=l+1:N/2
  243. S=(p-1)*n1;
  244. for i=1:n1
  245. if L(i,2)==3
  246. A(s+i,S+i)=-Es12(L(i+1,2))*wt(p)/2;
  247. else
  248. A(s+i,S+i)=-Es12(L(i,2))*wt(p)/2;
  249. end;
  250. end
  251. end
  252. end
  253. % Blocks from N/2 to N........................................
  254. a=0;
  255. for p=1:N/2
  256. S=(N/2+p-1)*n1;
  257. for i=2:n1
  258. if L(i,2)==3
  259. A(s+i,S+i-1)=-Es12(L(i+1,2))*wt(N/2+p)/2;
  260. else
  261. A(s+i,S+i-1)=-Es12(L(i,2))*wt(N/2+p)/2;
  262. end
  263. end
  264. a=a+(Es12(L(1,2))*wt(N/2+p)*yo/2);
  265. end
  266. B(s+1)=B(s+1)+a;
  267. end
  268. % Diagonal Block of matrix from N/2+1 to N.........................
  269. for t=N/2+1:N
  270. s=(t-1)*n1;
  271. A(s+1,s+1)=u(t)*(h(2)-h(1))/(h(2)*h(1))+Et12(L(1,2))-Es12(L(1,2))*wt(t)/2;
  272. A(s+1,s+1+1) = u(t)*h(1)/(h(2)*(h(1)+h(2)));
  273. %B(s+1)=u(t)*yo*h(2)/(h(1)*(h(1)+h(2)))+Q12(L(1,2))/2;
  274. B(s+1)=u(t)*yo*h(2)/(h(1)*(h(1)+h(2)));
  275. for i=2:n1-1
  276. if L(i+1,2)==3
  277. if i==2
  278. A(s+i,s+i)=u(t)*(1/h(i)+1/(h(i)+h(i-1)))+Et12(L(i,2))-Es12(L(i,2))*wt(t)/2;
  279. A(s+i,s+i-1)=-u(t)*(1/h(i)+1/h(i-1));
  280. if L(i,1)>=minleft && L(i,1)<=maxright
  281. B(s+i)=-u(t)*yo*(h(i)/(h(i-1)*(h(i)+h(i-1))))+Q12(L(i,2))/2;
  282. else
  283. B(s+i)=-u(t)*yo*(h(i)/(h(i-1)*(h(i)+h(i-1))));
  284. end
  285. else
  286. A(s+i,s+i)=u(t)*(1/h(i)+1/(h(i)+h(i-1)))+Et12(L(i,2))-Es12(L(i,2))*wt(t)/2;
  287. A(s+i,s+i-1)=-u(t)*(1/h(i)+1/h(i-1));
  288. A(s+i,s+i-2)=u(t)*(h(i)/(h(i-1)*(h(i)+h(i-1))));
  289. if L(i,1)>=minleft && L(i,1)<=maxright
  290. B(s+i)=Q12(L(i,2))/2;
  291. else
  292. B(s+i)=0.;
  293. end
  294. end
  295. else
  296. A(s+i,s+i-1)=-u(t)*h(i+1)/(h(i)*(h(i)+h(i+1)));
  297. A(s+i,s+i)=u(t)*(h(i+1)-h(i))/(h(i+1)*h(i))+Et12(L(i+1,2))-Es12(L(i+1,2))*wt(t)/2;
  298. A(s+i,s+i+1)=u(t)*h(i)/(h(i+1)*(h(i)+h(i+1)));
  299. if L(i+1,1)>=minleft && L(i+1,1)<=maxright
  300. B(s+i)=Q12(L(i+1,2))/2;
  301. else
  302. B(s+i)=0.;
  303. end
  304. end
  305. end
  306. A(s+n1,s+n1)=u(t)*(1/h(n1)+1/(h(n1)+h(n1-1)))+Et12(L(n1,2))-Es12(L(n1,2))*wt(t)/2;
  307. A(s+n1,s+n1-1)=-u(t)*(1/h(n1)+1/h(n1-1));
  308. A(s+n1,s+n1-2)=u(t)*(h(n1)/(h(n1-1)*(h(n1)+h(n1-1))));
  309. %B(s+n1)=Q12(L(n1,2))/2;
  310. B(s+n1)=0.;
  311. % Remaining Blocks in same direction up to N..............
  312. l=t;
  313. if l==N/2+1 && N>2
  314. for p=l+1:N
  315. S=(p-1)*n1;
  316. for i=1:n1
  317. if L(i+1,2)==3
  318. A(s+i,S+i)=-Es12(L(i,2))*wt(p)/2;
  319. else
  320. A(s+i,S+i)=-Es12(L(i+1,2))*wt(p)/2;
  321. end
  322. end
  323. end
  324. elseif l>N/2+1 && N>2
  325. for p=N/2+1:l-1
  326. S=(p-1)*n1;
  327. for i=1:n1
  328. if L(i+1,2)==3
  329. A(s+i,S+i)=-Es12(L(i,2))*wt(p)/2;
  330. else
  331. A(s+i,S+i)=-Es12(L(i+1,2))*wt(p)/2;
  332. end
  333. end
  334. end
  335. for p=l+1:N
  336. S=(p-1)*n1;
  337. for i=1:n1
  338. if L(i+1,2)==3
  339. A(s+i,S+i)=-Es12(L(i,2))*wt(p)/2;
  340. else
  341. A(s+i,S+i)=-Es12(L(i+1,2))*wt(p)/2;
  342. end
  343. end
  344. end
  345. end
  346. % Blocks from 1 to N/2........................................
  347. a=0;
  348. for p=1:N/2
  349. S=(p-1)*n1;
  350. for i=1:n1-1
  351. if L(i+1,2)==3
  352. A(s+i,S+i+1)=-Es12(L(i,2))*wt(p)/2;
  353. else
  354. A(s+i,S+i+1)=-Es12(L(i+1,2))*wt(p)/2;
  355. end
  356. end
  357. a=a+(Es12(L(n1,2))*wt(p)*y_/2);
  358. end
  359. B(s+n1)=B(s+n1)+a;
  360. end
  361. Z=A\B; % Solving for angular flux.
  362. return
  363. end
  364. %If you want to see Flux outgoing in positive directin see Y(n) or Z(M+1)
  365. %If you want see Flux scattering back see Y(n+1) or Z(M+2)