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

SN_per_bench_solver.m 5.5KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207
  1. %THIS IS THE AUXILIARY FUNCTION OF THE MAIN CODE "periodic_benchmark"
  2. %RICHARD VASQUES & NITIN KUMAR YADAV
  3. function[Z,extra,n1,B] = SN_per_bench_solver(T,m1,m2,n,N,Es1,Es2,Et1,Et2,yo,y_,Q1,Q2,u,wt,a)
  4. %Periodic Medium........
  5. interval=T/n;
  6. m12=zeros(2,1);
  7. m12(1)=m1;m12(2)=m2;
  8. mm=0;
  9. i=0;
  10. s=0;
  11. if a>1
  12. if a<m2/interval+2
  13. i=i+1;
  14. x1=(a-1)*interval;
  15. s=s+x1;
  16. x(i,1)=s;
  17. x(i,2)=2;
  18. else
  19. i=i+1;
  20. x1=(a-m2/interval-1)*interval;
  21. s=s+x1;
  22. x(i,1)=s;
  23. x(i,2)=mod(mm,2)+1;
  24. mm=mm+1;
  25. end
  26. end
  27. while s<T
  28. i=i+1;
  29. x1=m12(mod(mm,2)+1);
  30. s=s+x1;
  31. if s<=T
  32. x(i,1)=s;
  33. else
  34. x(i,1)=T;
  35. end
  36. x(i,2)=mod(mm,2)+1;
  37. mm=mm+1;
  38. end
  39. H=T/n;n1=1;i=1;j=1;
  40. extra=zeros(n,1);
  41. t1=i*H;
  42. L(1,1)=0;L(1,2)=x(j,2);
  43. if t1==x(j,1)
  44. extra(i)=1;
  45. h(n1)=x(j,1)/2;
  46. L(n1+1,1)=h(n1);
  47. L(n1+1,2)=x(j,2);
  48. h(n1+1)=h(n1);
  49. L(n1+2,1)=x(j,1);
  50. L(n1+2,2)=3;
  51. n1=n1+2;
  52. i=i+1;
  53. t1=i*H;
  54. else
  55. while t1<=x(j,1)
  56. h(n1)=H;
  57. L(n1+1,1)=i*H;
  58. if L(n1+1,1)==x(j,1) % checking is point is interface point?
  59. L(n1+1,2)=3;
  60. else
  61. L(n1+1,2)=x(j,2);
  62. end
  63. n1=n1+1;
  64. i=i+1;
  65. t1=i*H;
  66. end
  67. end
  68. j=2;
  69. if x(1,1)~=T
  70. while i<=n
  71. while t1<=x(j,1)
  72. L(n1+1,1)=i*H;
  73. if L(n1+1,1)==x(j,1) %checking is point is interface point?
  74. L(n1+1,2)=3;
  75. else
  76. L(n1+1,2)=x(j,2);
  77. end
  78. h(n1)=L(n1+1,1)-L(n1,1);
  79. n1=n1+1;
  80. i=i+1;
  81. t1=i*H;
  82. end
  83. j=j+1;
  84. if L(n1,1)==T && L(n1,2)==3 && L(n1-1,2)==3
  85. i=i-1;
  86. extra(i)=extra(i)+1;
  87. L(n1,1)=(x(j-1,1)+x(j-2,1))/2;
  88. L(n1,2)=x(j-1,2);
  89. h(n1-1)=L(n1,1)-L(n1-1,1);
  90. n1=n1+1;
  91. L(n1,1)=x(j-1,1);
  92. L(n1,2)=3;
  93. h(n1-1)=L(n1,1)-L(n1-1,1);
  94. i=i+1;
  95. end
  96. end
  97. end
  98. n1=n1-1;
  99. A=zeros(N*n1,N*n1); B=zeros(N*n1,1);
  100. Et12=zeros(2,1); Es12=zeros(2,1); Q12=zeros(2,1);
  101. Et12(1)=Et1;Et12(2)=Et2;
  102. Es12(1)=Es1;Es12(2)=Es2;
  103. Q12(1)=Q1;Q12(2)=Q2;
  104. % Diagonal Block of matrix up to N/2.................................
  105. for t=1:N/2
  106. s=(t-1)*n1;
  107. A(s+1,s+1)=-u(t)*(1/h(1))+Et12(L(1,2))-Es12(L(1,2))*wt(t)/2;
  108. A(s+1,s+2)=u(t)*(1/h(1));
  109. B(s+1)=Q12(L(1,2))/2;
  110. for i=2:n1-1
  111. if L(i,2)==3
  112. if i==n1-1
  113. 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;
  114. A(s+i,s+i+1)=u(t)*(1/h(i)+1/h(i+1));
  115. B(s+i)=u(t)*y_*(h(i)/(h(i+1)*(h(i)+h(i+1))))+Q12(L(i+1,2))/2;
  116. else
  117. 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;
  118. A(s+i,s+i+1)=u(t)*(1/h(i)+1/h(i+1));
  119. A(s+i,s+i+2)=-u(t)*(h(i)/(h(i+1)*(h(i)+h(i+1))));
  120. B(s+i)=Q12(L(i+1,2))/2;
  121. end
  122. else
  123. A(s+i,s+i-1)=-u(t)*h(i)/(h(i-1)*(h(i-1)+h(i)));
  124. 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;
  125. A(s+i,s+i+1)=u(t)*h(i-1)/(h(i)*(h(i-1)+h(i)));
  126. B(s+i)=Q12(L(i,2))/2;
  127. end
  128. end
  129. A(s+n1,s+n1-1)=-u(t)*h(n1)/(h(n1-1)*(h(n1-1)+h(n1)));
  130. 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;
  131. B(s+n1)=-u(t)*y_*h(n1-1)/(h(n1)*(h(n1-1)+h(n1)))+Q12(L(n1,2))/2;
  132. % Blocks from N/2 to N........................................
  133. a=0;
  134. for p=1:N/2
  135. S=(N/2+p-1)*n1;
  136. for i=2:n1
  137. if L(i,2)==3
  138. A(s+i,S+i-1)=-Es12(L(i+1,2))*wt(N/2+p)/2;
  139. else
  140. A(s+i,S+i-1)=-Es12(L(i,2))*wt(N/2+p)/2;
  141. end
  142. end
  143. a=a+(Es12(L(1,2))*wt(N/2+p)*yo/2);
  144. end
  145. B(s+1)=B(s+1)+a;
  146. end
  147. % Diagonal Block of matrix from N/2+1 to N.........................
  148. for t=N/2+1:N
  149. s=(t-1)*n1;
  150. 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;
  151. A(s+1,s+1+1) = u(t)*h(1)/(h(2)*(h(1)+h(2)));
  152. B(s+1)=u(t)*yo*h(2)/(h(1)*(h(1)+h(2)))+Q12(L(1,2))/2;
  153. for i=2:n1-1
  154. if L(i+1,2)==3
  155. if i==2
  156. 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;
  157. A(s+i,s+i-1)=-u(t)*(1/h(i)+1/h(i-1));
  158. B(s+i)=-u(t)*yo*(h(i)/(h(i-1)*(h(i)+h(i-1))))+Q12(L(i,2))/2;
  159. else
  160. 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;
  161. A(s+i,s+i-1)=-u(t)*(1/h(i)+1/h(i-1));
  162. A(s+i,s+i-2)=u(t)*(h(i)/(h(i-1)*(h(i)+h(i-1))));
  163. B(s+i)=Q12(L(i,2))/2;
  164. end
  165. else
  166. A(s+i,s+i-1)=-u(t)*h(i+1)/(h(i)*(h(i)+h(i+1)));
  167. 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;
  168. A(s+i,s+i+1)=u(t)*h(i)/(h(i+1)*(h(i)+h(i+1)));
  169. B(s+i)=Q12(L(i+1,2))/2;
  170. end
  171. end
  172. A(s+n1,s+n1)=u(t)*(1/h(n1))+Et12(L(n1,2))-Es12(L(n1,2))*wt(t)/2;
  173. A(s+n1,s+n1-1)=-u(t)*(1/h(n1));
  174. B(s+n1)=Q12(L(n1,2))/2;
  175. % Blocks from 1 to N/2........................................
  176. a=0;
  177. for p=1:N/2
  178. S=(p-1)*n1;
  179. for i=1:n1-1
  180. if L(i+1,2)==3
  181. A(s+i,S+i+1)=-Es12(L(i,2))*wt(p)/2;
  182. else
  183. A(s+i,S+i+1)=-Es12(L(i+1,2))*wt(p)/2;
  184. end
  185. end
  186. a=a+(Es12(L(n1,2))*wt(p)*y_/2);
  187. end
  188. B(s+n1)=B(s+n1)+a;
  189. end
  190. Z=A\B; % Solving for angular flux.
  191. return
  192. end