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

periodic_benchmark.m 3.7KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143
  1. %THIS CODE GENERATES ENSEMBLE-AVERAGED BENCHMARK RESULTS FOR THE
  2. %STEADY-STATE, MONOENERGETIC TRANSPORT EQUATION IN A PERIODIC MEDIUM
  3. %COMPOSED OF TWO MATERIALS WITH ISOTROPIC SCATTERING AND ISOTROPIC SOURCE
  4. %IN ROD GEOMETRY. THE TOTAL NUMBER OF REALIZATIONS IS DEFINED BY THE
  5. %SPATIAL DISCRETIZATION.
  6. %THIS CODE CALLS THE SOLVER ROUTINE "SN_per_bench_solver", USING 2-POINT CENTRAL
  7. %DIFFERENCING (ORDER H^2).
  8. clear all
  9. clc
  10. %------------------------------------------------------------
  11. %INPUTS
  12. T= 20; % total length of the system
  13. n= 2560; % # of points (discretization)
  14. yo=0.; % left boundary condition
  15. y_=0.; % right boundary condition
  16. m1= 1; % thickness of material 1 layers
  17. m2= 1; % thickness of material 2 layers
  18. Et1= 1.0; % total cross section Sigma_t1 of material 1
  19. cc= 0.5; % Scattering Ratio c1 of material 1
  20. Q1= 1.0; % homogeneous isotropic Source of material 1
  21. Et2=0.; % total cross section Sigma_t2 of material 2
  22. cc2= 0.; % Scattering Ratio c2 of material 2
  23. Q2= 0.; % homogeneous isotropic Source of material 2
  24. %------------------------------------------------------------
  25. %------------------------------------------------------------
  26. Es1=cc*Et1;
  27. Ea1= Et1-Es1;
  28. Es2=cc2*Et2;
  29. Ea2= Et2-Es2;
  30. %weights and directions
  31. N=2;
  32. wt(1)= 1;
  33. wt(2)=1;
  34. u(1)=-1;
  35. u(2)=1;
  36. a=1;
  37. reflec=0;reflec2=0;
  38. transm=0;transm2=0;
  39. SF=zeros(n+1,1);
  40. SF2=zeros(n+1,1);
  41. cond=0;
  42. %MAIN LOOP
  43. total = (m1+m2)/(T/n);
  44. while(a<total+1)
  45. fprintf('Problem %d of %d\n', a, total);
  46. [Z,extra,n1,B] = SN_per_bench_solver(T,m1,m2,n,N,Es1,Es2,Et1,Et2,yo,y_,Q1,Q2,u,wt,a);
  47. % Adjusting points..........................
  48. X=zeros(n*N,1);
  49. for i=1:N/2
  50. X((i-1)*n+1)=Z((i-1)*n1+1);
  51. k=2;
  52. for j=2:n
  53. k=k+extra(j-1);
  54. X((i-1)*n+j)=Z((i-1)*n1+k);
  55. k=k+1;
  56. end
  57. end
  58. for i=N/2+1:N
  59. k=1;
  60. for j=1:n
  61. k=k+extra(j);
  62. X((i-1)*n+j)=Z((i-1)*n1+k);
  63. k=k+1;
  64. end
  65. end
  66. Y=zeros(N*(n+1),1);
  67. % Adding boundary conditions...............
  68. i=1;
  69. for t=1:N/2
  70. s=(t-1);
  71. for j=i:t*n+s
  72. Y(j)=X(j-s);
  73. end
  74. Y(j+1)=y_;
  75. i=j+2;
  76. end
  77. i=i-1;
  78. for t=N/2+1:N
  79. Y(i+1)=yo;
  80. i=i+1;
  81. for j=i+1:t*n+t
  82. Y(j)=X(j-t);
  83. end
  84. i=j;
  85. end
  86. % Calculating Reflection, Transmission and Scalar Flux......
  87. RL=0;
  88. for t=1:N/2
  89. s=(t-1)*n;
  90. RL=RL+abs(wt(t)*u(t)*Y(s+t));
  91. end
  92. TR=0;
  93. for t=N/2+1:N
  94. s=(t-1)*n;
  95. TR=TR+wt(t)*u(t)*Y(s+n+t);
  96. end
  97. m=n+1;
  98. SCAL=zeros(m,1);
  99. for t=1:N
  100. s=(t-1)*m;
  101. for i=1:m
  102. SCAL(i)=SCAL(i)+wt(t)*Y(s+i);
  103. end
  104. end
  105. reflec=reflec+RL;
  106. reflec2=reflec2+(RL)^2;
  107. transm=transm+TR;
  108. transm2=transm2+(TR)^2;
  109. for i=1:n+1
  110. SF(i)=SF(i)+SCAL(i);
  111. SF2(i)=SF2(i)+(SCAL(i))^2;
  112. end
  113. a=a+1;
  114. end
  115. a=a-1;
  116. SF=SF./a;
  117. SF2=SF2./a;
  118. kk=T/n;
  119. p=0:kk:T; % p is line along x-axis.
  120. plot(p,SF,'r'); hold on
  121. plot(p,SF2,'b');
  122. save per_bench.mat