-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathA3_coreF.f90
1483 lines (1366 loc) · 58.3 KB
/
A3_coreF.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!
! Fortran 95 solver of eigenvalue probem
!
subroutine solve_eigenvalue_problem(meth, geom, nz, nx, ny, nt, ng, nmix, &
& flux, flux_a, imap, &
& sigt, sigtra, sigp, &
& nsigsn, fsigsn, tsigsn, sigsn, &
& nsign2n, fsign2n, tsign2n, sign2n, &
& chi, pitch, dz, keff, keff_a)
use omp_lib
implicit none
! method flag ('MC', 'DI')
character*2 meth
! geometry flag ('squar', 'hex01', 'hex06', 'hex24')
character*5 geom
! number of nodes in z, y, x dimensions and number of triangles per hexagon, number of groups and number of mixes
integer nz, nx, ny, nt, ng, nmix
! neutron flux array: flux(nz,nx,ny,nt,ng)
real*8 flux(:,:,:,:,:)
! adjoint flux array: flux_a(nz,nx,ny,nt,ng)
real*8 flux_a(:,:,:,:,:)
! map of material indexes: imap(nz,nx,ny)
integer imap(:,:,:)
! total cross section: sigt(nmix,ng)
real*8 sigt(:,:)
! transport cross section: sigtra(nmix,ng)
real*8 sigtra(:,:)
! production cross section: sigp(nmix,ng)
real*8 sigp(:,:)
! number of entries in full scattering cross section matrix: nsigsn(nmix)
integer nsigsn(:)
! index of energy group from which scattering occurs: fsigsn(nlgndr,nmix,max(nsigsn))
integer fsigsn(:,:,:)
! index of energy group to which scattering occurs: tsigsn(nlgndr,nmix,max(nsigsn))
integer tsigsn(:,:,:)
! full scattering cross section matrix: sigsn(nlgndr,nmix,max(nsigsn)))
real*8 sigsn(:,:,:)
! number of entries in n2n cross section matrix: sign2n(nmix)
integer nsign2n(:)
! index of energy group from which n2n occurs: fsign2n(nmix,max(nsign2n))
integer fsign2n(:,:)
! index of energy group to which n2n occurs: tsign2n(nmix,max(nsign2n))
integer tsign2n(:,:)
! n2n cross section matrix: sign2n(nmix,max(nsign2n)))
real*8 sign2n(:,:)
! fission spectrum: chi(nmix,ng)
real*8 chi(:,:)
! subassembly pitch
real*8 pitch
! axial nodalization: dz(nz-2)
real*8 dz(:)
! convergence flags
logical converge_k, converge_flux
! from and to indices for scattering matrix
integer f, t
! for do loop
integer iy, ix, iz, it, ig, indx, i, j
! index of mix in the node
integer imix
! counter of inner (flux) iterations
integer niteri
! counter of outer (fission source) iterations
integer nitero
! number of threads in case OMP is used
integer :: nthreads = 0
! absolute tolerance
real*8 atol
! node axial area-to-volume ratio
real*8 az_over_v
! node side area-to-volume ratio
real*8 aside_over_v
! sum of contributions from diffusion terms from neighbouring nodes
real*8 dif
! flux for evaluating the difference from the previous iteration
real*8 fluxnew
! multiplication factor: keff
real*8 keff(:)
! adjoint multiplication factor: keff_a
real*8 keff_a(:)
! multiplication factor for evaluating the difference from the previous iteration
real*8 knew
! multiplier at flux
real*8 mlt
! fission source
real*8 qf(nz,nx,ny,nt)
! fission source for adjoint problem
real*8 qff(nz,nx,ny,nt)
! total fission source
real*8 tfs
! fission source
real*8 qfis
! n2n source
real*8 qn2n
! scattering source
real*8 qs
! relative tolerance
real*8 rtol
! removal cross section
real*8 sigr
! sum of sigp for keff_a
real*8 sum_sigp
! Monte Carlo specific variables
! number of neutrons at the end of cycle and neutrons born at the beginning of cycle, number of inactive source cycles to skip before starting k-eff accumulation and number of active cycles for k-eff accumulation (to be done input parameters)
integer num_neutrons, num_neutrons_born, num_cycles_inactive, num_cycles_active
parameter(num_neutrons_born = 50000, num_cycles_inactive = 100, num_cycles_active = 1000)
! cycle index
integer icycle
! group number of neutron
integer igroup(num_neutrons_born*2)
! number of Legendre coefficients
integer :: nlgndr = 8
! number of scattering cosine bins
integer num_mubin
parameter(num_mubin = 200)
! neutron free path
real*8 free_path
! k-effective in a cycle
real*8 :: keff_active_cycle(num_cycles_active) = 1.0
! averaged k-effective
real*8 :: keff_expected(num_cycles_active) = 1.0
! scattering cosine
real*8 mu
! mu dependent pdf of scattering cross section from -> to
real*8 pdf_mu(nmix,num_mubin,ng,ng)
! Legendre coefficients
real*8 pl(8)
! pi number
real*8 :: pi = dacos(-1.0d0)
! k-effective standard deviation
real*8 :: sigma_keff(num_cycles_active) = 0.0d0
! total n2n cross section
real*8 sign2n_f(nmix,ng)
! n2n cross section for group ig
real*8 sign2n_ft(nmix,ng,ng)
! total scattering cross section for group f
real*8 sigs_f(nmix,ng)
! scattering cross section from -> to
real*8 sigs_ft(nmix,ng,ng),sigs1_ft(nmix,ng,ng)
! majorant
real*8 sigtmax(ng)
! virtual cross section
real*8 sigv
! neutron terminate probability
real*8 terminate_p
! angles and directions]
real*8 teta, phi, dir_x, dir_y, dir_z
! weight of neutrons
real*8 :: weight(num_neutrons_born*2) = 1.0d0
real*8 :: weight0(num_neutrons_born*2) = 1.0d0
! coordinates of hex centres
real*8 ynode(ny,nx), xnode(nx), zb(nz-1)
! coordinates of neutrons
real*8 y(num_neutrons_born*2), x(num_neutrons_born*2), z(num_neutrons_born*2)
logical absorbed
logical isotropic
logical virtual_collision
! temporal variables
integer itmp(num_neutrons_born*2), num_new, N, iactive, iyp
real*8 r, tmp, xtmp(num_neutrons_born*2), ytmp(num_neutrons_born*2), &
ztmp(num_neutrons_born*2), wtmp(num_neutrons_born*2), keff_cycle, siga, &
phi_prime, teta_prime, dmu, dcos_phi, dsin_phi, r1, &
dir_x_tmp, dir_y_tmp, dir_z_tmp
!$ call OMP_set_dynamic(.true.)
! change from python style to fortran style
imap = imap + 1
! scattering xs from group f to group t for mix imix
pdf_mu = 0.0d0
sigs_ft = 0.0d0
sigs1_ft = 0.0d0
sigs_f = 0.0d0
sign2n_ft = 0.0d0
sign2n_f = 0.0d0
do imix = 1,nmix
! scattering matrix
do indx = 1,nsigsn(imix)
f = fsigsn(2,imix,indx)+1
t = tsigsn(2,imix,indx)+1
sigs1_ft(imix,f,t) = sigs1_ft(imix,f,t) + sigsn(2,imix,indx)
f = fsigsn(1,imix,indx)+1
t = tsigsn(1,imix,indx)+1
sigs_ft(imix,f,t) = sigs_ft(imix,f,t) + sigsn(1,imix,indx)
! angular distribution of scattering pdf
mu = -1.0d0
dmu = 2.0d0/(num_mubin-1)
tmp = 0.0d0
do i = 1,num_mubin
do j = 1,nlgndr
if(j == 1)then
pl(1) = 1.0d0
else if(j == 2)then
pl(2) = mu
else
pl(j) = ( (2.0d0*dble(j)-1.0d0)*mu*pl(j-1) - (dble(j)-1.0d0)*pl(j-2) )/ dble(j)
end if
pdf_mu(imix,i,f,t) = pdf_mu(imix,i,f,t) + (dble(j) - 0.5d0)*sigsn(j,imix,indx)*pl(j)
end do
if(pdf_mu(imix,i,f,t) < 0.0d0)pdf_mu(imix,i,f,t) = 0.0d0
tmp = tmp + pdf_mu(imix,i,f,t)
mu = mu + dmu
end do
do i = 1,num_mubin
pdf_mu(imix,i,f,t) = pdf_mu(imix,i,f,t)/tmp
end do
end do
do ig = 1,ng
do j = 1,ng
sigs_f(imix,ig) = sigs_f(imix,ig) + sigs_ft(imix,ig,j)
end do
end do
! n2n matrix
do indx = 1,nsign2n(imix)
f = fsign2n(imix,indx)+1
t = tsign2n(imix,indx)+1
sign2n_ft(imix,f,t) = sign2n_ft(imix,f,t) + sign2n(imix,indx)
end do
do ig = 1,ng
do j = 1,ng
sign2n_f(imix,ig) = sign2n_f(imix,ig) + sign2n_ft(imix,ig,j)
end do
end do
end do
!open(1967,file='angular.txt')
!do imix = 1,nmix
! do indx = 1,nsigsn(imix)
! f = fsigsn(1,imix,indx)+1
! t = tsigsn(1,imix,indx)+1
! do i=1,num_mubin
! if(pdf_mu(imix,i,f,t) .ne. 0.0d0)write(1967,*)imix,f,t,pdf_mu(imix,i,f,t)
! end do
! write(1967,*)
! end do
!end do
!close(1967)
if(meth == 'MC')then
! MONTE CARLO SOLVER
! prepare flux array
do iz = 1,nz
do ix = 1,nx
do iy = 1,ny
do ig = 1,ng
flux(iz,ix,iy,1,ig) = 0.0d0
end do
end do
end do
end do
! find coordinate of the hexagon centres
do ix = 1,nx
xnode(ix) = dble(ix-1)*dsqrt(3.0d0)*pitch/2.0d0
do iy = 1,ny
if(mod(ix,2) == 0)then ! even row
ynode(iy,ix) = dble(iy-1)*pitch + pitch/2.0d0
else ! odd row
ynode(iy,ix) = dble(iy-1)*pitch
end if
end do
end do
zb(1) = 0.0d0
do iz = 2,nz-1
zb(iz) = zb(iz-1) + dz(iz-1)
end do
! define the majorant: the maximum total cross section vector
do ig = 1,ng
sigtmax(ig) = 0.0d0
do imix = 1,nmix
sigtmax(ig) = dmax1(sigtmax(ig), sigt(imix,ig))
end do
end do
! neutrons are assumed born randomly distributed in the core with weight 1 with sampled fission energy spectrum
num_neutrons = num_neutrons_born
do i = 1,num_neutrons_born
do
! generate random x coordinate
call random_number(r)
x(i) = (xnode(nx)+pitch)*r - 0.5d0*pitch
ix = minloc(dabs(xnode - x(i)),1)
! generate random y coordinate
call random_number(r)
y(i) = (ynode(ny,ix)+pitch)*r - 0.5d0*pitch
iy = minloc(dabs(ynode(:,ix) - y(i)),1)
! generate random z coordinate
call random_number(r)
z(i) = zb(nz-1)*r
iz = minloc(dabs(zb - z(i)),1)
if(z(i) > zb(iz)) iz = iz + 1
imix = imap(iz,ix,iy)
if(imix > 0 .and. sum(chi(imix,:)) > 0.0d0)then
weight(i) = 1.0d0
! sample group by comparing cumulative sum of fission spectrum with random number
call random_number(r)
tmp = 0.0d0
ig = ng + 1
do while(tmp <= r)
ig = ig - 1
tmp = tmp + chi(imix,ig)
end do
igroup(i) = ig
exit
end if
end do
end do
! main batch cycle
do icycle = 1,(num_cycles_inactive + num_cycles_active)
! normalize the weights of the neutrons to make the total weight equal to num_neutrons_born (equivalent to division by keff_cycle)
tmp = 0.0d0
do i = 1,num_neutrons
tmp = tmp + weight(i)
end do
do i = 1,num_neutrons
weight(i) = weight(i) / tmp * num_neutrons_born
weight0(i) = weight(i)
end do
!$omp parallel do default(shared) &
!$omp private(absorbed,isotropic,virtual_collision,r,free_path,teta,phi,mu, &
!$omp dir_x,dir_y,dir_z,dir_x_tmp,dir_y_tmp,dir_z_tmp, &
!$omp phi_prime,teta_prime,dcos_phi,dsin_phi, &
!$omp ix,iy,iyp,iz,imix,sigv,tmp,ig,siga,j,r1)
! loop over neutrons
do i = 1,num_neutrons
absorbed = .false.
isotropic = .true.
virtual_collision = .false.
! neutron random walk cycle: from emission to absorption
do while(.not. absorbed)
! sample free path length according to the Woodcock method
call random_number(r)
free_path = -dlog(r)/sigtmax(igroup(i))
if(.not. virtual_collision)then
if(isotropic)then
! sample the direction of neutron flight assuming isotropic fission
!call random_number(r)
!phi = 2.0d0*pi*r
!call random_number(r)
!teta = pi*r
!dir_x = dsin(teta)*dcos(phi)
!dir_y = dsin(teta)*dsin(phi)
!dir_z = dcos(teta)
! sample the direction of neutron flight assuming anisotropic scattering
r = 1.0d0
r1 = 1.0d0
do while(r**2+r1**2 > 1.0d0)
call random_number(r)
r = 2.d0*r - 1.0d0
call random_number(r1)
r1 = 2.d0*r1 - 1.0d0
end do
dir_x = 2.d0*r**2 + 2.d0*r1**2 - 1.0d0
dir_y = r*dsqrt((1.0d0-dir_x**2)/(r**2+r1**2))
dir_z = r1*dsqrt((1.0d0-dir_x**2)/(r**2+r1**2))
!dcos_phi = dcos(phi)
!dsin_phi = dsin(phi)
else
!! sample the direction of neutron flight assuming anisotropic scattering
!call random_number(r)
!phi_prime = 2.0d0*pi*r
!call random_number(r)
!teta_prime = pi*r
!!teta_prime = dacos(mu)
!
!tmp = dcos(teta_prime)*(1.0d0 - dcos(phi_prime))
!! Rodrigues formula (this is important for convergence not to have explicitly dcos(phi) and dsin(phi) in calculations)
!dir_x = dsin(teta+teta_prime)*dcos_phi*dcos(phi_prime) - &
! & dsin_phi*dsin(teta_prime)*dsin(phi_prime) + &
! & dir_x*tmp
!dir_y = dsin(teta+teta_prime)*dsin_phi*dcos(phi_prime) + &
! & dcos_phi*dsin(teta_prime)*dsin(phi_prime) + &
! & dir_y*tmp
!dir_z = dcos(teta+teta_prime)*dcos(phi_prime) + &
! & dir_z*tmp
!if(i == 100)write(*,*)dir_x**2 + dir_y**2 + dir_z**2
!! for next calculation of dir_x, dir_y and dir_z by formulas above
!teta = dacos(dir_z)
!dcos_phi = dir_x/dsin(teta)
!dsin_phi = dir_y/dsin(teta)
! sample the direction of neutron flight assuming anisotropic scattering
r = 1.0d0
r1 = 1.0d0
do while(r**2+r1**2 > 1.0d0)
call random_number(r)
r = 2.d0*r - 1.0d0
call random_number(r1)
r1 = 2.d0*r1 - 1.0d0
end do
dir_x_tmp = dir_x*mu + dsqrt(1.0d0 - mu**2)*(r*dir_x*dir_z - r1*dir_y)/dsqrt(r**2+r1**2)/dsqrt(1.d0-dir_z**2)
dir_y_tmp = dir_y*mu + dsqrt(1.0d0 - mu**2)*(r*dir_y*dir_z + r1*dir_x)/dsqrt(r**2+r1**2)/dsqrt(1.d0-dir_z**2)
dir_z_tmp = dir_z*mu - dsqrt(1.0d0 - mu**2)*r*dsqrt(1.d0-dir_z**2)/dsqrt(r**2+r1**2)
dir_x = dir_x_tmp
dir_y = dir_y_tmp
dir_z = dir_z_tmp
!if(i == 100)write(*,*)dir_x**2 + dir_y**2 + dir_z**2
end if
end if
! fly
x(i) = x(i) + free_path * dir_x
y(i) = y(i) + free_path * dir_y
z(i) = z(i) + free_path * dir_z
! find neutron position in structured mesh (ix, iy, iz)
ix = minloc(dabs(xnode - x(i)),1)
iy = minloc(dabs(ynode(:,ix) - y(i)),1)
if(xnode(ix) - x(i) > 0.5d0*pitch/dsqrt(3.0d0))then
iyp = minloc(dabs(ynode(:,ix-1) - y(i)),1)
if(dabs(ynode(iyp,ix-1) - y(i)) < dabs(ynode(iy,ix) - y(i)))then
ix = ix - 1
iy = iyp
end if
else if(x(i) - xnode(ix) > 0.5d0*pitch/dsqrt(3.0d0))then
iyp = minloc(dabs(ynode(:,ix+1) - y(i)),1)
if(dabs(ynode(iyp,ix+1) - y(i)) < dabs(ynode(iy,ix) - y(i)))then
ix = ix + 1
iy = iyp
end if
end if
iz = minloc(dabs(zb - z(i)),1)
if(z(i) > zb(iz)) iz = iz + 1
! mixture index
imix = imap(iz,ix,iy)
if(imix == 0)then
! neutron is out of core
weight(i) = 0.0d0
absorbed = .true.
else
! virtual cross section
sigv = sigtmax(igroup(i)) - sigt(imix,igroup(i))
! sample the type of the collision: virtual (do nothing) or real
call random_number(r)
if(sigv/sigtmax(igroup(i)) >= r)then
! virtual collision
virtual_collision = .true.
else
! real collision
virtual_collision = .false.
! score reactions with account for weight divided by the total cross section to estimate flux
if(icycle > num_cycles_inactive)then
flux(iz,ix,iy,1,igroup(i)) = flux(iz,ix,iy,1,igroup(i)) + weight(i)
end if
! sample type of the collision: scattering or absorption]
call random_number(r)
if(sigs_f(imix,igroup(i))/sigt(imix,igroup(i)) >= r)then
! anisotropic scattering
isotropic = .false.
! sample group for the secondary neutron by comparing cumulative sum of partial scattering xs with random number
call random_number(r)
tmp = 0.0d0
ig = ng + 1
do while(tmp <= r)
ig = ig - 1
tmp = tmp + sigs_ft(imix,igroup(i),ig)/sigs_f(imix,igroup(i))
end do
! sample scattering angle cosine mu
!call random_number(r)
!tmp = 0.0d0
!j = 0
!do while(tmp <= r)
! j = j + 1
! tmp = tmp + pdf_mu(imix,j,igroup(i),ig)
!end do
!mu = dble(j - 1)*2.0d0/dble(num_mubin - 1) - 1.0d0
! average scattering cosine
mu = sigs1_ft(imix,igroup(i),ig)/sigs_ft(imix,igroup(i),ig)
! sample scattering cosine mu isotropically
!call random_number(r)
!mu = 2.0d0*r - 1.0d0
igroup(i) = ig
else ! absorption
absorbed = .true.
siga = sigt(imix,igroup(i)) - sigs_f(imix,igroup(i))
! sample type of the collision: n2n or non-n2n absorption
call random_number(r)
if(sign2n_f(imix,igroup(i))/siga >= r)then ! isotropic n2n
! neutron is converted to the new neutron with the weight increased by 2
weight(i) = weight(i) * 2.0d0
! sample group for the secondary neutrons by comparing cumulative sum of partial n2n xs with random number
call random_number(r)
tmp = 0.0d0
ig = ng + 1
do while(tmp <= r)
ig = ig - 1
tmp = tmp + sign2n_ft(imix,igroup(i),ig)/sign2n_f(imix,igroup(i))
end do
igroup(i) = ig
else
! neutron is converted to the new fission neutron with the weight increased by eta
weight(i) = weight(i) * sigp(imix,igroup(i))/(siga - sign2n_f(imix,igroup(i)))
if(weight(i) > 0.0)then
! sample group for the new-born neutron by comparing cumulative sum of fission spectrum with random number
call random_number(r)
tmp = 0.0d0
ig = ng + 1
do while(tmp <= r)
ig = ig - 1
tmp = tmp + chi(imix,ig)
end do
igroup(i) = ig
end if
end if ! n2n or non-n2n absorption
end if ! scattering or absorption
end if ! virtual or real
end if ! imix <= 0 or imix > 0
end do ! random walk cycle
end do ! neutron cycle
!$omp end parallel do
! Russian roulette
do i = 1,num_neutrons
if(weight(i) > 0.0d0 .and. weight(i) < 0.5d0 .and. weight(i) < weight0(i))then
terminate_p = 1.0d0 - weight(i)/weight0(i)
call random_number(r)
if(terminate_p >= r)then
! kill
weight(i) = 0.0d0
else
! restore the weight
weight(i) = weight0(i)
end if
end if
end do
! clean up absorbed or killed neutrons
j = 0
do i = 1,num_neutrons
if(weight(i) > 0.0d0)then
j = j + 1
xtmp(j) = x(i)
ytmp(j) = y(i)
ztmp(j) = z(i)
itmp(j) = igroup(i)
wtmp(j) = weight(i)
end if
end do
num_neutrons = j
x = 0.0
y = 0.0
z = 0.0
igroup = 0
weight = 0.0
do i = 1,num_neutrons
x(i) = xtmp(i)
y(i) = ytmp(i)
z(i) = ztmp(i)
igroup(i) = itmp(i)
weight(i) = wtmp(i)
end do
! split too "heavy" neutrons
num_new = 0
do i = 1,num_neutrons
if(weight(i) > 1.0d0)then
! truncated integer value of the neutron weight
N = floor(weight(i))
! sample the number of split neutrons
call random_number(r)
if(weight(i)-dble(N) >= r) N = N + 1
if(N > 1)then
! change the weight of the split neutron
weight(i) = weight(i)/dble(N)
! introduce new neutrons
do j = 1,N-1
num_new = num_new + 1
x(num_neutrons + num_new) = x(i)
y(num_neutrons + num_new) = y(i)
z(num_neutrons + num_new) = z(i)
weight(num_neutrons + num_new) = weight(i)
igroup(num_neutrons + num_new) = igroup(i)
end do
end if
end if
end do
! increase the number of neutrons
num_neutrons = num_neutrons + num_new
! k-eff in a cycle equals the total weight of the new generation over the total weight of the old generation (the old generation weight = num_neutrons_born)
keff_cycle = 0.0d0
do i = 1,num_neutrons
keff_cycle = keff_cycle + weight(i)
end do
keff_cycle = keff_cycle / num_neutrons_born
iactive = icycle - num_cycles_inactive
if(iactive <= 0)then
write(*,"('Inactive cycle = ',i4,'/',i4,' | k-eff cycle = ',f8.5,' | num_neutrons = ',i6)")&
icycle, num_cycles_inactive, keff_cycle, num_neutrons
else
! k-effective of the cycle
keff_active_cycle(iactive) = keff_cycle
! k-effective of the problem
keff_expected(iactive) = 0.0d0
do i = 1, iactive
keff_expected(iactive) = keff_expected(iactive) + keff_active_cycle(i)
end do
keff_expected(iactive) = keff_expected(iactive)/iactive
! standard deviation of k-effective
sigma_keff(iactive) = 0.0d0
do i = 1, iactive
sigma_keff(iactive) = sigma_keff(iactive) + ( keff_active_cycle(i) - keff_expected(iactive) )**2
end do
sigma_keff(iactive) = dsqrt(sigma_keff(iactive) / dble(max(iactive-1,1)) / dble(iactive) )
write(*,"('Active cycle = ',i6,'/',i6,' | k-eff cycle = ',f8.5,' | num_neutrons = ',i6,&
& ' | k-eff expected = ',f9.5,' | sigma = ',f9.5)")&
& icycle-num_cycles_inactive, num_cycles_active, keff_cycle, num_neutrons,&
& keff_expected(iactive), sigma_keff(iactive)
end if
end do ! batch cycle
do iz = 1,nz
do ix = 1,nx
do iy = 1,ny
imix = imap(iz,ix,iy)
if(imix > 0)then
do ig = 1,ng
flux(iz,ix,iy,1,ig) = flux(iz,ix,iy,1,ig) / sigt(imix,ig) / dz(iz-1)
end do
end if
end do
end do
end do
else
! DIFFUSION SOLVER
! verification test homogeneous cube
ng = 2
sigtra(1,1) = 0.2468
sigtra(1,2) = 0.3084
nsigsn(1) = 1
fsigsn(1,1,1) = 0
tsigsn(1,1,1) = 1
sigsn(1,1,1) = 2.3e-3
nsign2n(1) = 0
sigp(1,1) = 2.41*2.42e-4
sigp(1,2) = 2.41*4.08e-3
sigt(1,1) = 1.382e-3 + sigsn(1,1,1)
sigt(1,2) = 5.4869e-3
chi(1,1) = 1
chi(1,2) = 0
! relative tolerance
rtol = 1.0e-8
! absolute tolerance
atol = 1.0e-8
! initialize fission source
qf = 1.
! side area to volume ratio of control volume
if(geom == 'squar')then
aside_over_v = 1./pitch
else if(geom == 'hex01')then
aside_over_v = 2./(3.*pitch)
else if(geom == 'hex06')then
aside_over_v = 4./pitch
else ! geom == 'hex24'
aside_over_v = 8./pitch
end if
keff(1) = 1.
! convergence flag
converge_k = .false.
! outer iteration counter
nitero = 0
do while(.not. converge_k .and. nitero < 1000)
nitero = nitero + 1
! initialize flux convergence flag
converge_flux = .false.
! inner iteration counter
niteri = 0
do while(.not. converge_flux)
niteri = niteri + 1
converge_flux = .true.
!$omp parallel do default(shared) &
!$omp private(imix,az_over_v,mlt,dif,sigr,qs,f,t,qn2n,qfis,fluxnew)
do iz = 1,nz
do ix = 1,nx
do iy = 1,ny
!$ if(omp_get_thread_num() == 0 .and. iz == 1 .and. ix == 1 .and. iy == 1) nthreads = OMP_get_num_threads()
! if (iz, ix, iy) is not a boundary condition node, i.e. not 0 (vac) and not -1 (ref)
imix = imap(iz,ix,iy)
if(imix > 0)then
! node axial area-to-volume ratio
az_over_v = 1./dz(iz-1)
do it = 1,nt
do ig = 1,ng
mlt = 0.
dif = 0.
! diffusion terms (mlt and dif) in different nodalizations and different directions
call mltdif(mlt, dif, iz, ix, iy, it, ig, imix, imap, nz, nx, ny, nt, ng, nmix, &
pitch, dz, sigtra, flux, aside_over_v, az_over_v, geom)
! removal xs
sigr = sigt(imix,ig)
do indx = 1,nsigsn(imix)
f = fsigsn(1,imix,indx)+1
t = tsigsn(1,imix,indx)+1
if(f == ig .and. t == ig)then
sigr = sigr - sigsn(1,imix,indx)
end if
end do
mlt = mlt + sigr
! scattering source
qs = 0.
do indx = 1,nsigsn(imix)
f = fsigsn(1,imix,indx)+1
t = tsigsn(1,imix,indx)+1
if(f .ne. ig .and. t == ig)then
qs = qs + sigsn(1,imix,indx) * flux(iz,ix,iy,it,f)
end if
end do
! n2n source
qn2n = 0.
do indx = 1,nsign2n(imix)
f = fsign2n(imix,indx)+1
t = tsign2n(imix,indx)+1
if(f .ne. ig .and. t == ig)then
qn2n = qn2n + 2.*sign2n(imix,indx)*flux(iz,ix,iy,it,f)
end if
end do
! fission source
qfis = chi(imix,ig)*qf(iz,ix,iy,it)/keff(1)
! neutron flux
fluxnew = (dif + qs + qn2n + qfis)/mlt
if(converge_flux)then
converge_flux = abs(fluxnew - flux(iz,ix,iy,it,ig)) < rtol*abs(fluxnew) + atol
end if
flux(iz,ix,iy,it,ig) = fluxnew
end do
end do
end if
end do
end do
end do
!$omp end parallel do
end do
! calculate node-wise fission source qf
!$omp parallel do default(shared) private(imix)
do iz = 1,nz
do ix = 1,nx
do iy = 1,ny
! if (iz, ix, iy) is not a boundary condition node, i.e. not 0 (vac) and not -1 (ref)
imix = imap(iz,ix,iy)
if(imix > 0)then
do it = 1,nt
qf(iz,ix,iy,it) = 0.
do ig = 1,ng
qf(iz,ix,iy,it) = qf(iz,ix,iy,it) + sigp(imix,ig)*flux(iz,ix,iy,it,ig)
end do
end do
end if
end do
end do
end do
!$omp end parallel do
! calculate total fission source tfs
tfs = 0.
do iz = 1,nz
do ix = 1,nx
do iy = 1,ny
! if (iz, ix, iy) is not a boundary condition node, i.e. not 0 (vac) and not -1 (ref)
imix = imap(iz,ix,iy)
if(imix > 0)then
do it = 1,nt
tfs = tfs + qf(iz,ix,iy,it)
end do
end if
end do
end do
end do
! new k-effective is the ratio of total fission sources at the current (tfs) and previous (1.0) iterations
knew = tfs
converge_k = abs(knew - keff(1)) < rtol*abs(knew) + atol
keff(1) = knew
if(nthreads == 0)then
write(*,'("keff: ",f13.6, " | niteri: ", i3, " | nitero: ", i3, " | ")') &
keff(1),niteri,nitero
else
write(*,'("keff: ",f13.6, " | niteri: ", i3, " | nitero: ", i3, " | OMPthreads: ", i3, " | ")') &
keff(1),niteri,nitero,nthreads
end if
if(isnan(keff(1))) stop
end do
! ADJOINT DIFFUSION SOLVER
! eigenvalue keff_a equal to ratio of total fission source at two iterations.
! flux is normalise to total fission source = 1 at previous iteration
keff_a(1) = 1.
! convergence flag
converge_k = .false.
! outer iteration counter
nitero = 0
do while(.not. converge_k .and. nitero < 1000)
nitero = nitero + 1
! initialize flux convergence flag
converge_flux = .false.
! inner iteration counter
niteri = 0
do while(.not. converge_flux)
niteri = niteri + 1
converge_flux = .true.
!$omp parallel do default(shared) &
!$omp private(imix,az_over_v,mlt,dif,sigr,qs,f,t,qn2n,qfis,fluxnew)
do iz = 1,nz
do ix = 1,nx
do iy = 1,ny
!$ if(omp_get_thread_num() == 0 .and. iz == 1 .and. ix == 1 .and. iy == 1) nthreads = OMP_get_num_threads()
! if (iz, ix, iy) is not a boundary condition node, i.e. not 0 (vac) and not -1 (ref)
imix = imap(iz,ix,iy)
if(imix > 0)then
! node axial area-to-volume ratio
az_over_v = 1./dz(iz-1)
do it = 1,nt
do ig = 1,ng
mlt = 0.
dif = 0.
! diffusion terms (mlt and dif) in different nodalizations and different directions
call mltdif(mlt, dif, iz, ix, iy, it, ig, imix, imap, nz, nx, ny, nt, ng, nmix, &
pitch, dz, sigtra, flux_a, aside_over_v, az_over_v, geom)
! removal xs
sigr = sigt(imix,ig)
do indx = 1,nsigsn(imix)
f = fsigsn(1,imix,indx)+1
t = tsigsn(1,imix,indx)+1
if(f == ig .and. t == ig)then
sigr = sigr - sigsn(1,imix,indx)
end if
end do
mlt = mlt + sigr
! scattering source
qs = 0.
do indx = 1,nsigsn(imix)
f = fsigsn(1,imix,indx)+1
t = tsigsn(1,imix,indx)+1
if(t .ne. ig .and. f == ig)then
qs = qs + sigsn(1,imix,indx) * flux_a(iz,ix,iy,it,t)
end if
end do
! n2n source
qn2n = 0.
do indx = 1,nsign2n(imix)
f = fsign2n(imix,indx)+1
t = tsign2n(imix,indx)+1
if(t .ne. ig .and. f == ig)then
qn2n = qn2n + 2.*sign2n(imix,indx)*flux_a(iz,ix,iy,it,t)
end if
end do
! fission source
qfis = sigp(imix,ig)*qf(iz,ix,iy,it)/keff_a(1)
! adjoint flux
fluxnew = (dif + qs + qn2n + qfis)/mlt
if(converge_flux)then
converge_flux = abs(fluxnew - flux_a(iz,ix,iy,it,ig)) < rtol*abs(fluxnew) + atol
end if
flux_a(iz,ix,iy,it,ig) = fluxnew
end do
end do
end if
end do
end do
end do
!$omp end parallel do
end do
! calculate node-wise fission source qf and total fission source tfs
!$omp parallel do default(shared) private(imix)
do iz = 1,nz
do ix = 1,nx
do iy = 1,ny
! if (iz, ix, iy) is not a boundary condition node, i.e. not 0 (vac) and not -1 (ref)
imix = imap(iz,ix,iy)
if(imix > 0)then
do it = 1,nt
qf(iz,ix,iy,it) = 0.
sum_sigp = 0.
do ig = 1,ng
qf(iz,ix,iy,it) = qf(iz,ix,iy,it) + chi(imix,ig)*flux_a(iz,ix,iy,it,ig)
sum_sigp = sum_sigp + sigp(imix,ig)
end do
qff(iz,ix,iy,it) = qf(iz,ix,iy,it)*sum_sigp
end do
end if
end do
end do
end do
!$omp end parallel do
! calculate total fission source tfs
tfs = 0.
do iz = 1,nz
do ix = 1,nx
do iy = 1,ny
! if (iz, ix, iy) is not a boundary condition node, i.e. not 0 (vac) and not -1 (ref)
imix = imap(iz,ix,iy)
if(imix > 0)then
do it = 1,nt
tfs = tfs + qff(iz,ix,iy,it)
end do
end if
end do
end do
end do
! new k-effective is the ratio of total fission sources at the current (tfs) and previous (1.0) iterations
knew = tfs
converge_k = abs(knew - keff_a(1)) < rtol*abs(knew) + atol
keff_a(1) = knew
if(nthreads == 0)then
write(*,'("keff_a: ",f13.6, " | niteri: ", i3, " | nitero: ", i3, " | ")') &
keff_a(1),niteri,nitero
else
write(*,'("keff_a: ",f13.6, " | niteri: ", i3, " | nitero: ", i3, " | OMPthreads: ", i3, " | ")') &
keff_a(1),niteri,nitero,nthreads
end if
if(isnan(keff_a(1))) stop
end do
end if
! change from fortran style to python style
imap = imap - 1
end subroutine
!--------------------------------------------------------------------------------------------------
! Fortran 95 solver of reactor kinetics problem
subroutine solve_kinetic_problem(keff, geom, nz, nx, ny, nt, ng, nmix, flux, dfidt, imap, &
& sigt, sigtra, sigp, &
& nsigsn, fsigsn, tsigsn, sigsn, &
& nsign2n, fsign2n, tsign2n, sign2n, &
& chi, pitch, dz, cdnp, dcdnpdt, betaeff, dnplmb, ndnp)
use omp_lib
implicit none
! geometry flag ('squar', 'hex01', 'hex06', 'hex24')
character*5 geom
! number of nodes in z, y, x dimensions and number of triangles per hexagon, number of groups and number of mixes
integer nz, nx, ny, nt, ng, nmix
! map of material indexes: imap(nz,nx,ny)
integer imap(:,:,:)
! subassembly pitch
real*8 pitch
! axial nodalization: dz(nz-2)
real*8 dz(:)
! for do loop
integer iy, ix, iz, it, ig, indx, i, j, im
! index of mix in the node