注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

AlexYoung

做好大家都必须要做的事,发展自己感兴趣的事,结束圣人所不齿的事...

 
 
 

日志

 
 

二维圆柱(R,Theta)几何的单群SN(N=4,6,8,10,12,14,16)计算程序  

2011-06-20 00:33:46|  分类: 课程练习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

!***********************************************************************
!*              二维圆柱(R,Theta)几何的单群SN计算程序                *
!*                         ——《核反应堆数值分析与模拟》期末大作业报告*
!*                                                                     *
!*                                       /杨万奎  2011.6.20            *
!***********************************************************************
!有参考解的计算网格为R方向5个网格,Theta方向4个网格,故在调试阶段做相同设
!置;采用S6计算,飞行方向坐标系的每个象限6个离散方向;采用B.G.Carlson等权
!EQn求积集
,由其1971年提供的从N=4,6,8,10,12,14,16共七组求积集可供使用。单
!群调试成功后试图!加上扩散综合加速(DSA),这需要在增加网格点数的基础上才
!能体现出来加速效!果。相关参数符号定义如下:
!
!                 R1=116cm---------裂变区半径;
!                 R2=145cm---------反射层外半径;
!                 NI、NJ-----------分别为r和Theta方向的离散网格数;
!                 NM---------------飞行方向坐标系一个坐标上的离散点数;
!                 NF---------------为反射层的网格数;
!                 NT---------------每个象限的离散点数;
program SN2DCyl_Uniform
implicit none
integer,parameter::NF=1, NI=5*NF,  NJ=4,  NM=16, NT=NM*(NM+2)/8
real,parameter::pi=3.14159265, R1=116,  R2=145,  delta_theta=2*pi/NJ
real,dimension(NI,NJ)::sigma_s,  sigma_t,  vsigma_f, phi,phiold=10,keff=1
real:: mu(4*NT),eta(4*NT), Psi_c(NI,NJ,4*NT),Psi_br(NI+1,NJ,4*NT), &
  Psi_btheta(NI,NJ+1,4*NT),Psi_bangle(NI,NJ,4*NT+1), &
  S(NI,NJ,4*NT),E,F,G,r_1,r_2,tempkeff, totalkeff,guiyi, &
  alpha(4*NT+1), V(NI,NJ), A(NI+1), B(NI),totalkeffold=1,totalkeff2=1, &
  fissionsold,fissions,tol=1,time_begin, time_end, P_weight
real:: p(6),q(8),t(11),x(14),y(17)
integer::i,j,k,kk,ds,ds1, m,countd=1,ct
integer,dimension(NM/2)::m31,m32,m41,m42,m21,m22,m11,m12
!=======================================================================
call cpu_time(time_begin)
select case(NM)
 case(4)
 mu(1:12)=[-0.3500212, -0.8688903, -0.3500212, &
    0.3500212, 0.3500212, 0.8688903, &
    -0.3500212, -0.8688903, -0.3500212, &
    0.3500212, 0.3500212, 0.8688903]
 eta(1:12)=[-0.3500212,-0.3500212, -0.8688903, &
    -0.3500212,-0.3500212, -0.8688903, &
    0.3500212,0.3500212, 0.8688903, &
    0.3500212,0.3500212, 0.8688903]
 P_weight=1.0/3
 case(6)
mu(1:24)=[-0.2561429,-0.6815646,-0.2663443,-0.9320846,-0.6815646,-0.2561429, &
   0.2561429,0.2663443,0.6815646,0.2561429,0.6815646,0.9320846, &
  -0.2561429,-0.6815646,-0.2663443,-0.9320846,-0.6815646,-0.2561429, &
   0.2561429,0.2663443,0.6815646,0.2561429,0.6815646,0.9320846]
eta(1:24)=[-0.2561429,-0.2663443,-0.6815646,-0.2561429,-0.6815646,-0.9320846, &
   -0.2561429,-0.6815646,-0.2663443,-0.9320846,-0.6815646,-0.2561429, &
   0.2561429,0.2663443,0.6815646,0.2561429,0.6815646,0.9320846, &
   0.2561429,0.6815646,0.2663443,0.9320846,0.6815646,0.2561429]
 P_weight=1.0/6
 case(8)
 p(1:6)=[0.1971380,0.2133981,0.5512958,0.5773503,0.8065570,0.9603506]
 mu(1:40)=[-p(1),-p(3),-p(2),-p(5),-p(3),-p(2),-p(6),-p(5),-p(3),-p(1), &
    p(1),p(2),p(3),p(2),p(3),p(5),p(1),p(3),p(5),p(6), &
    -p(1),-p(3),-p(2),-p(5),-p(3),-p(2),-p(6),-p(5),-p(3),-p(1), &
    p(1),p(2),p(3),p(2),p(3),p(5),p(1),p(3),p(5),p(6)]
 eta(1:40)=[-p(1),-p(2),-p(3),-p(2),-p(3),-p(5),-p(1),-p(3),-p(5),-p(6), &
    -p(1),-p(3),-p(2),-p(5),-p(3),-p(2),-p(6),-p(5),-p(3),-p(1), &
    p(1),p(2),p(3),p(2),p(3),p(5),p(1),p(3),p(5),p(6), &
    p(1),p(3),p(2),p(5),p(3),p(2),p(6),p(5),p(3),p(1)]
 P_weight=1.0/10
 case(10)
 q(1:8)=[0.1631408,0.1755273,0.4567576,0.4897749,0.6961286,0.7212773,0.8721024,0.9730212]
 mu(1:60)=[-q(1),-q(3),-q(2),-q(5),-q(4),-q(2),-q(7),-q(6),-q(4),-q(2),-q(8),-q(7),-q(5),-q(3),-q(1), &
    q(1),q(2),q(3),q(2),q(4),q(5),q(2),q(4),q(6),q(7),q(1),q(3),q(5),q(7),q(8), &
    -q(1),-q(3),-q(2),-q(5),-q(4),-q(2),-q(7),-q(6),-q(4),-q(2),-q(8),-q(7),-q(5),-q(3),-q(1), &
    q(1),q(2),q(3),q(2),q(4),q(5),q(2),q(4),q(6),q(7),q(1),q(3),q(5),q(7),q(8)]
 eta(1:60)=[-q(1),-q(2),-q(3),-q(2),-q(4),-q(5),-q(2),-q(4),-q(6),-q(7),-q(1),-q(3),-q(5),-q(7),-q(8), &
    -q(1),-q(3),-q(2),-q(5),-q(4),-q(2),-q(7),-q(6),-q(4),-q(2),-q(8),-q(7),-q(5),-q(3),-q(1), &
    q(1),q(2),q(3),q(2),q(4),q(5),q(2),q(4),q(6),q(7),q(1),q(3),q(5),q(7),q(8), &
    q(1),q(3),q(2),q(5),q(4),q(2),q(7),q(6),q(4),q(2),q(8),q(7),q(5),q(3),q(1)]
 P_weight=1.0/15
 case(12)
 t(1:11)=[0.1370611,0.1497456,0.3911744,0.4213515,0.4249785,0.6040252,0.6400755,0.7827706,0.8030727,0.9080522,0.9810344]
 mu(1:84)=[-t(1),-t(3),-t(2),-t(6),-t(4),-t(2),-t(8),-t(7),-t(5),-t(2),-t(10),-t(9),-t(7),-t(4),-t(2),-t(11),-t(10),-t(8),-t(6),-t(3),-t(1), &
    t(1),t(2),t(3),t(2),t(4),t(6),t(2),t(5),t(7),t(8),t(2),t(4),t(7),t(9),t(10),t(1),t(3),t(6),t(8),t(10),t(11), &
    -t(1),-t(3),-t(2),-t(6),-t(4),-t(2),-t(8),-t(7),-t(5),-t(2),-t(10),-t(9),-t(7),-t(4),-t(2),-t(11),-t(10),-t(8),-t(6),-t(3),-t(1), &
    t(1),t(2),t(3),t(2),t(4),t(6),t(2),t(5),t(7),t(8),t(2),t(4),t(7),t(9),t(10),t(1),t(3),t(6),t(8),t(10),t(11)]
 eta(1:84)=[-t(1),-t(2),-t(3),-t(2),-t(4),-t(6),-t(2),-t(5),-t(7),-t(8),-t(2),-t(4),-t(7),-t(9),-t(10),-t(1),-t(3),-t(6),-t(8),-t(10),-t(11),&
    -t(1),-t(3),-t(2),-t(6),-t(4),-t(2),-t(8),-t(7),-t(5),-t(2),-t(10),-t(9),-t(7),-t(4),-t(2),-t(11),-t(10),-t(8),-t(6),-t(3),-t(1), &
    t(1),t(2),t(3),t(2),t(4),t(6),t(2),t(5),t(7),t(8),t(2),t(4),t(7),t(9),t(10),t(1),t(3),t(6),t(8),t(10),t(11),&
    t(1),t(3),t(2),t(6),t(4),t(2),t(8),t(7),t(5),t(2), t(10),t(9),t(7),t(4),t(2),t(11),t(10),t(8),t(6),t(3),t(1)]
 P_weight=1.0/21
 case(14)
 x(1:14)=[0.1196230,0.1301510,0.3399238,0.3700559,0.3736108,0.5326134,0.5691823,0.5773503,0.7010923,0.7324250,0.8362916,0.8521252,0.9314035,0.9855865]
 mu(1:112)=[-x(1),-x(3),-x(2),-x(6),-x(4),-x(2),-x(9),-x(7),-x(5),-x(2),-x(11),-x(10),-x(8),-x(5),-x(2),-x(13),-x(12),-x(10),-x(7),-x(4),-x(2),-x(14),-x(13),-x(11),-x(9),-x(6),-x(3),-x(1), &
    x(1),x(2),x(3),x(2),x(4),x(6),x(2),x(5),x(7),x(9),x(2),x(5),x(8),x(10),x(11),x(2),x(4),x(7),x(10),x(12),x(13),  x(1),x(3),x(6),x(9),x(11),x(13),x(14), &
    -x(1),-x(3),-x(2),-x(6),-x(4),-x(2),-x(9),-x(7),-x(5),-x(2),-x(11),-x(10),-x(8),-x(5),-x(2),-x(13),-x(12),-x(10),-x(7),-x(4),-x(2),-x(14),-x(13),-x(11),-x(9),-x(6),-x(3),-x(1), &
    x(1),x(2),x(3),x(2),x(4),x(6),x(2),x(5),x(7),x(9),x(2),x(5),x(8),x(10),x(11),x(2),x(4),x(7),x(10),x(12),x(13),x(1),x(3),x(6),x(9),x(11),x(13),x(14)]
 eta(1:112)=[-x(1),-x(2),-x(3),-x(2),-x(4),-x(6),-x(2),-x(5),-x(7),-x(9),-x(2),-x(5),-x(8),-x(10),-x(11),-x(2),-x(4),-x(7),-x(10),-x(12),-x(13),-x(1),-x(3),-x(6),-x(9),-x(11),-x(13),-x(14), &
    -x(1),-x(3),-x(2),-x(6),-x(4),-x(2),-x(9),-x(7),-x(5),-x(2),-x(11),-x(10),-x(8),-x(5),-x(2),-x(13),-x(12),-x(10),-x(7),-x(4),-x(2),-x(14),-x(13),-x(11),-x(9),-x(6),-x(3),-x(1), &
    x(1),x(2),x(3),x(2),x(4),x(6),x(2),x(5),x(7),x(9),x(2),x(5),x(8),x(10),x(11),x(2),x(4),x(7),x(10),x(12),x(13),x(1),x(3),x(6),x(9),x(11),x(13),x(14), &
    x(1),x(3),x(2),x(6),x(4),x(2),x(9),x(7),x(5),x(2),x(11),x(10),x(8),x(5),x(2),x(13),x(12),x(10),x(7),x(4),x(2),x(14),x(13),x(11),x(9),x(6),x(3),x(1)]
 P_weight=1.0/28
 case(16)
 y(1:17)=[0.1050159,0.1152880,0.3016701,0.3284315,0.3332906,0.4743525,0.5107319,0.5215431,0.6327389,0.6666774,0.6752671,0.7657351,0.7925089,0.8727534,0.8855877,0.9464163,0.9889102]
 mu(1:144)=[-y(1),-y(3),-y(2),-y(6),-y(4),-y(2),-y(9),-y(7),-y(5),-y(2),-y(12),-y(10),-y(8),-y(5),-y(2),-y(14),-y(13),-y(11),-y(8),-y(5),-y(2),-y(16),-y(15),-y(13),-y(10),-y(7),-y(4),-y(2),-y(17),-y(16),-y(14),-y(12),-y(9),-y(6),-y(3),-y(1), &
    y(1),y(2),y(3),y(2),y(4),y(6),y(2),y(5),y(7),y(9),y(2),y(5),y(8),y(10),y(12),y(2),y(5),y(8),y(11),y(13),y(14),y(2),y(4),y(7),y(10),y(13),y(15),y(16),y(1),y(3),y(6),y(9),y(12),y(14),y(16),y(17), &
    -y(1),-y(3),-y(2),-y(6),-y(4),-y(2),-y(9),-y(7),-y(5),-y(2),-y(12),-y(10),-y(8),-y(5),-y(2),-y(14),-y(13),-y(11),-y(8),-y(5),-y(2),-y(16),-y(15),-y(13),-y(10),-y(7),-y(4),-y(2),-y(17),-y(16),-y(14),-y(12),-y(9),-y(6),-y(3),-y(1), &
    y(1),y(2),y(3),y(2),y(4),y(6),y(2),y(5),y(7),y(9),y(2),y(5),y(8),y(10),y(12),y(2),y(5),y(8),y(11),y(13),y(14),y(2),y(4),y(7),y(10),y(13),y(15),y(16),y(1),y(3),y(6),y(9),y(12),y(14),y(16),y(17)]
 eta(1:144)=[y(1),y(2),y(3),y(2),y(4),y(6),y(2),y(5),y(7),y(9),y(2),y(5),y(8),y(10),y(12),y(2),y(5),y(8),y(11),y(13),y(14),y(2),y(4),y(7),y(10),y(13),y(15),y(16),y(1),y(3),y(6),y(9),y(12),y(14),y(16),y(17), &
    -y(1),-y(3),-y(2),-y(6),-y(4),-y(2),-y(9),-y(7),-y(5),-y(2),-y(12),-y(10),-y(8),-y(5),-y(2),-y(14),-y(13),-y(11),-y(8),-y(5),-y(2),-y(16),-y(15),-y(13),-y(10),-y(7),-y(4),-y(2),-y(17),-y(16),-y(14),-y(12),-y(9),-y(6),-y(3),-y(1), &
    y(1),y(2),y(3),y(2),y(4),y(6),y(2),y(5),y(7),y(9),y(2),y(5),y(8),y(10),y(12),y(2),y(5),y(8),y(11),y(13),y(14),y(2),y(4),y(7),y(10),y(13),y(15),y(16),y(1),y(3),y(6),y(9),y(12),y(14),y(16),y(17), &
    y(1),y(3),y(2),y(6),y(4),y(2),y(9),y(7),y(5),y(2),y(12),y(10),y(8),y(5),y(2),y(14),y(13),y(11),y(8),y(5),y(2),y(16),y(15),y(13),y(10),y(7),y(4),y(2),y(17),y(16),y(14),y(12),y(9),y(6),y(3),y(1)]
 P_weight=1.0/36
 case default
 print*,'The value of NM must be 4、6、8、10、12、14 or 16 !'
 stop
end select

!Parameters Initinization!
m31(1)=1;m32(1)=1;
do k=2,NM/2
 m31(k)=m31(k-1)+k-1
 m32(k)=m32(k-1)+k
end do
m41=m31+NT;m42=m32+NT;
m21=m41+NT;m22=m42+NT;
m11=m21+NT;m12=m22+NT;
!Region I CrossSection
sigma_s(1:NI-NF,:)=0.8804
sigma_t(1:NI-NF,:)=1.0386
vsigma_f(1:NI-NF,:)=0.2902
!Region II CrossSection
sigma_s(NI-NF+1:NI,:)=1.6245
sigma_t(NI-NF+1:NI,:)=1.6445
vsigma_f(NI-NF+1:NI,:)=0.0
!alpha
alpha(1)=0.0
alpha(4*NT+1)=0.0
do m=2,4*NT
 alpha(m)=alpha(m-1)-mu(m-1)*P_weight
end do
!A,B,V
do j=1,NJ
 do i=1,NI
 r_1=(i-1)*R2/NI
 r_2=i*R2/NI
 A(i)=r_1*delta_theta
 B(i)=r_2-r_1
 !注意教材公式(2.98)是错误的,应用theta角度间隔!
 V(i,j)=delta_theta*(r_2**2-r_1**2)/2
 end do
end do
 A(NI+1)=R2*delta_theta
!Vaccum boundary condition (exact and correct value)
Psi_br(NI+1,:,1:NT)=0.0
Psi_br(NI+1,:,2*NT+1 : 3*NT)=0.0
!first iteration setup (source & psi_btheta())
!S=1
psi_c=1.0
do j=1,NJ
do i=1,NI
 phi(i,j)=2.0*4.0*pi*sum(psi_c(i,j,:))*P_weight/8.0
 do m=1,4*NT
 S(i,j,m)=(sigma_s(i,j)+vsigma_f(i,j)/totalkeff2)*phi(i,j)/(4.0*pi)
 end do
end do
end do
fissionsold=sum(vsigma_f(:,:)*phi(:,:)*V(:,:))
!=======================================================================
do while(countd>0)
!do while(tol>=1e-6)
!do ds1=1,3000
 !----------------------第III象限--------------------
 do ds=1,10  !让其迭代赋值10次
 do j=NJ,1,-1
 do i=NI,1,-1
 do m=1,NT
  E=abs(mu(m))*(A(i)+A(i+1))
  F=2.0*abs(eta(m))*B(i)
  G=(A(i+1)-A(i))*(alpha(m)+alpha(m+1))*6.0
   do k=1,NM/2
   if(m==m31(k)) then
   G=0
  Psi_c(i,j,m)=(E*Psi_br(i+1,j,m)+F*Psi_btheta(i,j+1,m)+G* &
   Psi_bangle(i,j,m)+V(i,j)*S(i,j,m))/(E+F+G+V(i,j)*sigma_t(i,j))
  Psi_bangle(i,j,m+1)=max(Psi_c(i,j,m), 0.0)
   goto 101
   end if
   enddo
  Psi_c(i,j,m)=(E*Psi_br(i+1,j,m)+F*Psi_btheta(i,j+1,m)+G* &
   Psi_bangle(i,j,m)+V(i,j)*S(i,j,m))/(E+F+G+V(i,j)*sigma_t(i,j))
  Psi_bangle(i,j,m+1)=max(2.0*Psi_c(i,j,m)-Psi_bangle(i,j,m), 0.0)
101  Psi_br(i,j,m)=max(2.0*Psi_c(i,j,m)-Psi_br(i+1,j,m), 0.0)
  Psi_btheta(i,j,m)=max(2.0*Psi_c(i,j,m)-Psi_btheta(i,j+1,m), 0.0)
  end do
 end do
 end do
  forall(i=1:NI,m=1:NT) Psi_btheta(i,NJ+1,m)=Psi_btheta(i,1,m) !周期边界条件
 end do

 !传给第IV象限 Psi_bangle
 do j=NJ,1,-1
 do i=NI,1,-1
  do k=NM/2,1,-1 !注意这个指标必须从大到小
  Psi_bangle(i,j,m41(k))=Psi_bangle(i,j,m32(k)+1)
  enddo
 end do
 end do

 !第IV象限计算初值 Psi_br,由圆柱中心点处的反射边界条件给出!(必须要的)
 !看教材编号图上各点的对称关系!
 do j=1,NJ
  do k=1,NM/2
  do kk=0,k-1
  Psi_br(1,j,m42(k)-kk)=Psi_br(1,j,m31(k)+kk)
  enddo
  enddo
 end do

 !--------------------------第IV象限-------------------------------
 do ds=1,10
 do j=NJ,1,-1
 do i=1,NI
 do m=NT+1,2*NT
  E=abs(mu(m))*(A(i)+A(i+1))
  F=2.0*abs(eta(m))*B(i)
  G=(A(i+1)-A(i))*(alpha(m)+alpha(m+1))*6.0
  Psi_c(i,j,m)=(E*Psi_br(i,j,m)+F*Psi_btheta(i,j+1,m)+ &
 G*Psi_bangle(i,j,m)+V(i,j)*S(i,j,m))/(E+F+G+V(i,j)*sigma_t(i,j))
   do k=1,NM/2
   if(m==m42(k)) then
    goto 110
   end if
   enddo
  Psi_bangle(i,j,m+1)=max(2*Psi_c(i,j,m)-Psi_bangle(i,j,m), 0.0)
110  Psi_br(i+1,j,m)=max(2.0*Psi_c(i,j,m)-Psi_br(i,j,m), 0.0)
  Psi_btheta(i,j,m)=max(2.0*Psi_c(i,j,m)-Psi_btheta(i,j+1,m), 0.0)
 end do
 end do
 end do
  forall(i=1:NI,m=NT+1:2*NT) Psi_btheta(i,NJ+1,m)= Psi_btheta(i,1,m)
 end do


 !-----------------------第II象限--------------------
 do ds=1,10
 do j=1,NJ
 do i=NI,1,-1
 do m=2*NT+1,3*NT
  E=abs(mu(m))*(A(i)+A(i+1))
  F=2.0*abs(eta(m))*B(i)
  G=(A(i+1)-A(i))*(alpha(m)+alpha(m+1))*6.0
   do k=1,NM/2
   if(m==m21(k)) then
   G=0
  Psi_c(i,j,m)=(E*Psi_br(i+1,j,m)+F*Psi_btheta(i,j,m)+ &
 G*Psi_bangle(i,j,m)+V(i,j)*S(i,j,m))/(E+F+G+V(i,j)*sigma_t(i,j))
  Psi_bangle(i,j,m+1)=max(Psi_c(i,j,m), 0.0)
   goto 102
   end if
   enddo
  Psi_c(i,j,m)=(E*Psi_br(i+1,j,m)+F*Psi_btheta(i,j,m)+ &
 G*Psi_bangle(i,j,m)+V(i,j)*S(i,j,m))/(E+F+G+V(i,j)*sigma_t(i,j))
  Psi_bangle(i,j,m+1)=max(2.0*Psi_c(i,j,m)-Psi_bangle(i,j,m), 0.0)
102  Psi_br(i,j,m)=max(2.0*Psi_c(i,j,m)-Psi_br(i+1,j,m), 0.0)
  Psi_btheta(i,j+1,m)=max(2.0*Psi_c(i,j,m)-Psi_btheta(i,j,m), 0.0)
  end do
 end do
 end do
  forall(i=1:NI,m=2*NT+1:3*NT) Psi_btheta(i,1,m)= Psi_btheta(i,NJ+1,m)
   !注这与III、IV的不同
 end do


 !传给第I卦限 Psi_bangle
 do j=NJ,1,-1
 do i=NI,1,-1
  do k=NM/2,1,-1
  Psi_bangle(i,j,m11(k))=Psi_bangle(i,j,m22(k)+1)
  enddo
 end do
 end do


 !传给第I卦限 Psi_br
 do j=1,NJ
  do k=1,NM/2
  do kk=0,k
  Psi_br(1,j,m12(k)-kk)=Psi_br(1,j,m21(k)+kk)
  enddo
  enddo
 end do

 !-----------------------第I象限---------------------
 do ds=1,10
 do j=1,NJ
 do i=1,NI                                                                                                                                                                                                                                                                 
 do m=3*NT+1,4*NT
  E=abs(mu(m))*(A(i)+A(i+1))
  F=2.0*abs(eta(m))*B(i)
  G=(A(i+1)-A(i))*(alpha(m)+alpha(m+1))*6.0
  Psi_c(i,j,m)=(E*Psi_br(i,j,m)+F*Psi_btheta(i,j,m)+ &
 G*Psi_bangle(i,j,m)+V(i,j)*S(i,j,m))/(E+F+G+V(i,j)*sigma_t(i,j))
   do k=1,NM/2
   if(m==m12(k)) then 
    goto 111 
   endif
   enddo
  Psi_bangle(i,j,m+1)=max(2*Psi_c(i,j,m)-Psi_bangle(i,j,m),0.0)
111  Psi_br(i+1,j,m)=max(2.0*Psi_c(i,j,m)-Psi_br(i,j,m), 0.0)
  Psi_btheta(i,j+1,m)=max(2.0*Psi_c(i,j,m)-Psi_btheta(i,j,m),0.0)
 end do
 end do
 end do
  forall(i=1:NI,m=3*NT+1:4*NT) Psi_btheta(i,1,m)= Psi_btheta(i,NJ+1,m)
 end do

!求标通量及keff(两类keff:临界keff和真实keff)
countd=0
do j=1,NJ
do i=1,NI
 phi(i,j)=2.0*4.0*pi*sum(Psi_c(i,j,:))*P_weight/8.0
 tempkeff=phi(i,j)/phiold(i,j)
 !这里的收敛标准非常重要,为什么要与1进行比较?因为对裂变源项
 !除以keff之后,即认为系统是临界的,这就是我们可以用临界方程解非
 !临界稳态问题的理论基础......
 if (abs(tempkeff-1)/1>=1e-4) then
  countd=countd+1
 end if
 keff(i,j)=tempkeff
end do
end do
totalkeff=sum(phi)/sum(phiold)
phiold=phi
!totalkeff-----为临界方程中加入keff本征值之后的临界系统的有效增值系数,
!              可用于表征通量收敛。
!totalkeff2----为真实的有效增值系数,可用于表征keff收敛。
fissions=sum(vsigma_f(:,:)*phi(:,:)*V(:,:))
totalkeff2=fissions/fissionsold*totalkeffold
tol=abs(totalkeff2-totalkeffold)/totalkeff2
totalkeffold=totalkeff2
fissionsold=fissions

ct=ct+1
write(*,*) "iteration times=",ct,"   totalkeff=", totalkeff2,totalkeff
!----------------------------------------------------------------------
!Calculate Source
do i=1,NI
do j=1,NJ
do m=1,4*NT
S(i,j,m)=(sigma_s(i,j)+vsigma_f(i,j)/totalkeff2)*phi(i,j)/(4.0*pi)
end do
end do
end do

end do

!================================输出结果===============================
write(*,*) phi(:,1)
open(unit=11,file='Results.txt')
open(unit=12,file='scalar_flux.txt')
write(12,'(5f12.6)') ((phi(i,j),i=1,NI),j=1,NJ)

write(11,*)"scalar flux"
write(11,'(5f12.6)')((phi(i,j),i=1,NI),j=1,NJ)

write(11,*)"vector flux(Psi_c)"
write(11,'(5f12.6)') (((Psi_c(i,j,k),i=1,NI),j=1,NJ),k=1,4*NT)

write(11,*)"keff = ", totalkeff2

call cpu_time(time_end)
print *,'clapsed time is ',time_end - time_begin,'    secondes'
print *,'Press return carriage key to exit...'
read(*,*) !Wait for a return carriage key to exit current window.
stop
end

  评论这张
 
阅读(552)| 评论(11)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017