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

SN_hom.m 4.8KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212
  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. clear all
  9. clc
  10. % INPUTS--------------------------------
  11. display('Default is ROD geometry!!');
  12. display(' ');
  13. rod_slab=input('Enter 1 if you want to change to SLAB Geometry: ');
  14. if rod_slab==1
  15. N=1;
  16. while floor(N/2)~=ceil(N/2)
  17. N=input( 'Enter the number of Discrete ordinate directions: ');
  18. end
  19. clc
  20. display('SOLVING PROBLEM IN SLAB GEOMETRY!!');
  21. display(' ')
  22. else
  23. rod_slab=0;
  24. clc
  25. display('SOLVING PROBLEM IN ROD GEOMETRY!!');
  26. display(' ');
  27. N=2;
  28. end
  29. T=input('Enter the total length of the system: ');
  30. n=1;
  31. while n<4
  32. n=input('Enter at how many points you want to calculate: ');
  33. end
  34. Et=input('Enter the total cross section Sigma_t : ');
  35. cc=10;
  36. while cc>1 || cc<0
  37. cc=input('Enter the Scattering Ratio c (between 0 and 1): ');
  38. end
  39. Es = cc*Et;
  40. yo=input('Enter the boundary value in the positive direction: ');
  41. y_=input('Enter the boundary value in the negative direction: ');
  42. Q=input('Enter the homogeneous isotropic Source: ');
  43. M=n*N; h=T/n;
  44. A=zeros(M,M); B(1:M)=Q/2;
  45. % GAUSS-LEGENDRE QUADRATURE............................
  46. beta = (1:N-1)./sqrt(4*((1:N-1)).^2-1);
  47. [w,x] = eig(diag(beta,-1)+diag(beta,1));
  48. u = diag(x);
  49. wt = 2*w(1,:)'.^2;
  50. if rod_slab~=1
  51. u(1)=-1;u(2)=1;
  52. end
  53. % Diagonal Block of matrix up to N/2.................................
  54. for t=1:N/2
  55. s=(t-1)*n;
  56. A(s+1,s+1)=(-11*u(t)/(6*h)+Et-Es*wt(t)/2);
  57. A(s+1,s+2) = 3*u(t)/h;
  58. A(s+1,s+3) = -3*u(t)/(2*h);
  59. A(s+1,s+4) = u(t)/(3*h);
  60. for i=2:n-1
  61. A(s+i,s+i-1) = -u(t)/(2*h);
  62. A(s+i,s+i)=(Et-Es*wt(t)/2);
  63. A(s+i,s+i+1)= u(t)/(2*h);
  64. end
  65. A(s+n,s+n-1) = -u(t)/(2*h);
  66. A(s+n,s+n)=(Et-Es*wt(t)/2);
  67. B(s+n)=-u(t)*y_/(2*h)+Q/2;
  68. % Remaining Blocks in same direction up to N/2..............
  69. l=t;
  70. if l==1 && N>2
  71. for p=l+1:N/2
  72. S=(p-1)*n;
  73. for i=1:n
  74. A(s+i,S+i)=-Es*wt(p)/2;
  75. end
  76. end
  77. elseif l>1 && N>2
  78. for p=1:l-1
  79. S=(p-1)*n;
  80. for i=1:n
  81. A(s+i,S+i)=-Es*wt(p)/2;
  82. end
  83. end
  84. for p=l+1:N/2
  85. S=(p-1)*n;
  86. for i=1:n
  87. A(s+i,S+i)=-Es*wt(p)/2;
  88. end
  89. end
  90. end
  91. % Blocks from N/2 to N........................................
  92. a=0;
  93. for p=1:N/2
  94. S=(N/2+p-1)*n;
  95. for i=2:n
  96. A(s+i,S+i-1)=-Es*wt(N/2+p)/2;
  97. end
  98. a=a+(Es*wt(N/2+p)*yo/2);
  99. end
  100. B(s+1)=a+Q/2;
  101. end
  102. % Diagonal Block of matrix from N/2+1 to N.........................
  103. for t=N/2+1:N
  104. s=(t-1)*n;
  105. A(s+1,s+1)=(Et-Es*wt(t)/2);
  106. A(s+1,s+1+1) = u(t)/(2*h);
  107. B(s+1)=u(t)*yo/(2*h)+Q/2;
  108. for i=2:n-1
  109. A(s+i,s+i-1) = -u(t)/(2*h);
  110. A(s+i,s+i)=(Et-Es*wt(t)/2);
  111. A(s+i,s+i+1)= u(t)/(2*h);
  112. end
  113. A(s+n,s+n-3)=-u(t)/(3*h);
  114. A(s+n,s+n-2)=3*u(t)/(2*h);
  115. A(s+n,s+n-1)=-3*u(t)/h;
  116. A(s+n,s+n) =(11*u(t)/(6*h)+Et-Es*wt(t)/2);
  117. % Remaining Blocks in same direction up to N..............
  118. l=t;
  119. if l==N/2+1 && N>2
  120. for p=l+1:N
  121. S=(p-1)*n;
  122. for i=1:n
  123. A(s+i,S+i)=-Es*wt(p)/2;
  124. end
  125. end
  126. elseif l>N/2+1 && N>2
  127. for p=N/2+1:l-1
  128. S=(p-1)*n;
  129. for i=1:n
  130. A(s+i,S+i)=-Es*wt(p)/2;
  131. end
  132. end
  133. for p=l+1:N
  134. S=(p-1)*n;
  135. for i=1:n
  136. A(s+i,S+i)=-Es*wt(p)/2;
  137. end
  138. end
  139. end
  140. % Blocks from 1 to N/2........................................
  141. a=0;
  142. for p=1:N/2
  143. S=(p-1)*n;
  144. for i=1:n-1
  145. A(s+i,S+i+1)=-Es*wt(p)/2;
  146. end
  147. a=a+(Es*wt(p)*y_/2);
  148. end
  149. B(s+n)=a+Q/2;
  150. end
  151. X=A\B';
  152. Y=zeros(M+N,1);
  153. % Adding boundary conditions to the array...............
  154. i=1;
  155. for t=1:N/2
  156. s=(t-1);
  157. for j=i:t*n+s
  158. Y(j)=X(j-s);
  159. end
  160. Y(j+1)=y_;
  161. i=j+2;
  162. end
  163. i=i-1;
  164. for t=N/2+1:N
  165. Y(i+1)=yo;
  166. i=i+1;
  167. for j=i+1:t*n+t
  168. Y(j)=X(j-t);
  169. end
  170. i=j;
  171. end
  172. % Calculating Reflection, Transmission and Scalar Flux...........
  173. RL=0;
  174. for t=1:N/2
  175. s=(t-1)*n;
  176. RL=RL+abs(wt(t)*u(t)*Y(s+t));
  177. end
  178. TR=0;
  179. for t=N/2+1:N
  180. s=(t-1)*n;
  181. S=(t-N/2-1)*n;
  182. k=N/2+1;
  183. TR=TR+wt(t)*u(t)*Y(s+n+t);
  184. end
  185. m=n+1;
  186. SCAL=zeros(m,1);
  187. for t=1:N
  188. s=(t-1)*m;
  189. for i=1:m
  190. SCAL(i)=SCAL(i)+wt(t)*Y(s+i);
  191. end
  192. end
  193. x=0:h:T;
  194. %plot(x,SCAL,'--');
  195. plot(x,SCAL,'b'); hold on