forked from openmc-dev/openmc
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtally.F90
3015 lines (2499 loc) · 118 KB
/
tally.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
module tally
#ifdef MPI
use message_passing
#endif
use algorithm, only: binary_search
use constants
use error, only: fatal_error
use geometry_header
use global
use math, only: t_percentile, calc_pn, calc_rn
use mesh, only: get_mesh_bin, bin_to_mesh_indices, &
get_mesh_indices, mesh_indices_to_bin, &
mesh_intersects_1d, mesh_intersects_2d, &
mesh_intersects_3d
use mesh_header, only: RegularMesh
use output, only: header
use particle_header, only: LocalCoord, Particle
use string, only: to_str
use tally_header, only: TallyResult
use tally_filter
implicit none
integer :: position(N_FILTER_TYPES - 3) = 0 ! Tally map positioning array
!$omp threadprivate(position)
procedure(score_general_), pointer :: score_general => null()
procedure(score_analog_tally_), pointer :: score_analog_tally => null()
abstract interface
subroutine score_general_(p, t, start_index, filter_index, i_nuclide, &
atom_density, flux)
import Particle
import TallyObject
type(Particle), intent(in) :: p
type(TallyObject), intent(inout) :: t
integer, intent(in) :: start_index
integer, intent(in) :: i_nuclide
integer, intent(in) :: filter_index ! for % results
real(8), intent(in) :: flux ! flux estimate
real(8), intent(in) :: atom_density ! atom/b-cm
end subroutine score_general_
subroutine score_analog_tally_(p)
import Particle
type(Particle), intent(in) :: p
end subroutine score_analog_tally_
end interface
contains
!===============================================================================
! INIT_TALLY_ROUTINES Sets the procedure pointers needed for minimizing code
! with the CE and MG modes.
!===============================================================================
subroutine init_tally_routines()
if (run_CE) then
score_general => score_general_ce
score_analog_tally => score_analog_tally_ce
else
score_general => score_general_mg
score_analog_tally => score_analog_tally_mg
end if
end subroutine init_tally_routines
!===============================================================================
! SCORE_GENERAL* adds scores to the tally array for the given filter and
! nuclide. This function is called by all volume tallies. For analog tallies,
! the flux estimate depends on the score type so the flux argument is really
! just used for filter weights. The atom_density argument is not used for
! analog tallies.
!===============================================================================
subroutine score_general_ce(p, t, start_index, filter_index, i_nuclide, &
atom_density, flux)
type(Particle), intent(in) :: p
type(TallyObject), intent(inout) :: t
integer, intent(in) :: start_index
integer, intent(in) :: i_nuclide
integer, intent(in) :: filter_index ! for % results
real(8), intent(in) :: flux ! flux estimate
real(8), intent(in) :: atom_density ! atom/b-cm
integer :: i ! loop index for scoring bins
integer :: l ! loop index for nuclides in material
integer :: m ! loop index for reactions
integer :: q ! loop index for scoring bins
integer :: i_temp ! temperature index
integer :: i_nuc ! index in nuclides array (from material)
integer :: i_energy ! index in nuclide energy grid
integer :: score_bin ! scoring bin, e.g. SCORE_FLUX
integer :: score_index ! scoring bin index
integer :: d ! delayed neutron index
integer :: g ! delayed neutron index
integer :: k ! loop index for bank sites
integer :: d_bin ! delayed group bin index
integer :: dg_filter ! index of delayed group filter
real(8) :: yield ! delayed neutron yield
real(8) :: atom_density_ ! atom/b-cm
real(8) :: f ! interpolation factor
real(8) :: score ! analog tally score
real(8) :: E ! particle energy
i = 0
SCORE_LOOP: do q = 1, t % n_user_score_bins
i = i + 1
! determine what type of score bin
score_bin = t % score_bins(i)
! determine scoring bin index
score_index = start_index + i
!#########################################################################
! Determine appropirate scoring value.
select case(score_bin)
case (SCORE_FLUX, SCORE_FLUX_YN)
if (t % estimator == ESTIMATOR_ANALOG) then
! All events score to a flux bin. We actually use a collision
! estimator in place of an analog one since there is no way to count
! 'events' exactly for the flux
if (survival_biasing) then
! We need to account for the fact that some weight was already
! absorbed
score = p % last_wgt + p % absorb_wgt
else
score = p % last_wgt
end if
score = score / material_xs % total * flux
else
! For flux, we need no cross section
score = flux
end if
case (SCORE_TOTAL, SCORE_TOTAL_YN)
if (t % estimator == ESTIMATOR_ANALOG) then
! All events will score to the total reaction rate. We can just
! use the weight of the particle entering the collision as the
! score
if (survival_biasing) then
! We need to account for the fact that some weight was already
! absorbed
score = p % last_wgt + p % absorb_wgt * flux
else
score = p % last_wgt * flux
end if
else
if (i_nuclide > 0) then
score = micro_xs(i_nuclide) % total * atom_density * flux
else
score = material_xs % total * flux
end if
end if
case (SCORE_INVERSE_VELOCITY)
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
E = p % E
else
E = p % last_E
end if
if (t % estimator == ESTIMATOR_ANALOG) then
! All events score to an inverse velocity bin. We actually use a
! collision estimator in place of an analog one since there is no way
! to count 'events' exactly for the inverse velocity
if (survival_biasing) then
! We need to account for the fact that some weight was already
! absorbed
score = p % last_wgt + p % absorb_wgt
else
score = p % last_wgt
end if
! Score the flux weighted inverse velocity with velocity in units of
! cm/s
score = score / material_xs % total &
/ (sqrt(TWO * E / (MASS_NEUTRON_MEV)) * C_LIGHT * 100.0_8) * flux
else
! For inverse velocity, we don't need a cross section. The velocity is
! in units of cm/s.
score = flux / (sqrt(TWO * E / (MASS_NEUTRON_MEV)) * C_LIGHT * 100.0_8)
end if
case (SCORE_SCATTER, SCORE_SCATTER_N)
if (t % estimator == ESTIMATOR_ANALOG) then
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP
! Since only scattering events make it here, again we can use
! the weight entering the collision as the estimator for the
! reaction rate
score = p % last_wgt * flux
else
! Note SCORE_SCATTER_N not available for tracklength/collision.
if (i_nuclide > 0) then
score = (micro_xs(i_nuclide) % total &
- micro_xs(i_nuclide) % absorption) * atom_density * flux
else
score = (material_xs % total - material_xs % absorption) * flux
end if
end if
case (SCORE_SCATTER_PN)
! Only analog estimators are available.
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) then
i = i + t % moment_order(i)
cycle SCORE_LOOP
end if
! Since only scattering events make it here, again we can use
! the weight entering the collision as the estimator for the
! reaction rate
score = p % last_wgt * flux
case (SCORE_SCATTER_YN)
! Only analog estimators are available.
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) then
i = i + (t % moment_order(i) + 1)**2 - 1
cycle SCORE_LOOP
end if
! Since only scattering events make it here, again we can use
! the weight entering the collision as the estimator for the
! reaction rate
score = p % last_wgt * flux
case (SCORE_NU_SCATTER, SCORE_NU_SCATTER_N)
! Only analog estimators are available.
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP
! For scattering production, we need to use the pre-collision weight
! times the yield as the estimate for the number of neutrons exiting a
! reaction with neutrons in the exit channel
if (p % event_MT == ELASTIC .or. p % event_MT == N_LEVEL .or. &
(p % event_MT >= N_N1 .and. p % event_MT <= N_NC)) then
! Don't waste time on very common reactions we know have
! multiplicities of one.
score = p % last_wgt * flux
else
m = nuclides(p % event_nuclide) % reaction_index % &
get_key(p % event_MT)
! Get yield and apply to score
associate (rxn => nuclides(p % event_nuclide) % reactions(m))
score = p % last_wgt * flux &
* rxn % products(1) % yield % evaluate(p % last_E)
end associate
end if
case (SCORE_NU_SCATTER_PN)
! Only analog estimators are available.
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) then
i = i + t % moment_order(i)
cycle SCORE_LOOP
end if
! For scattering production, we need to use the pre-collision
! weight times the yield as the estimate for the number of
! neutrons exiting a reaction with neutrons in the exit channel
if (p % event_MT == ELASTIC .or. p % event_MT == N_LEVEL .or. &
(p % event_MT >= N_N1 .and. p % event_MT <= N_NC)) then
! Don't waste time on very common reactions we know have multiplicities
! of one.
score = p % last_wgt * flux
else
m = nuclides(p%event_nuclide)%reaction_index% &
get_key(p % event_MT)
! Get yield and apply to score
associate (rxn => nuclides(p % event_nuclide) % reactions(m))
score = p % last_wgt * flux &
* rxn % products(1) % yield % evaluate(p % last_E)
end associate
end if
case (SCORE_NU_SCATTER_YN)
! Only analog estimators are available.
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) then
i = i + (t % moment_order(i) + 1)**2 - 1
cycle SCORE_LOOP
end if
! For scattering production, we need to use the pre-collision
! weight times the yield as the estimate for the number of
! neutrons exiting a reaction with neutrons in the exit channel
if (p % event_MT == ELASTIC .or. p % event_MT == N_LEVEL .or. &
(p % event_MT >= N_N1 .and. p % event_MT <= N_NC)) then
! Don't waste time on very common reactions we know have multiplicities
! of one.
score = p % last_wgt * flux
else
m = nuclides(p%event_nuclide)%reaction_index% &
get_key(p % event_MT)
! Get yield and apply to score
associate (rxn => nuclides(p%event_nuclide)%reactions(m))
score = p % last_wgt * flux &
* rxn % products(1) % yield % evaluate(p % last_E)
end associate
end if
case (SCORE_ABSORPTION)
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
! No absorption events actually occur if survival biasing is on --
! just use weight absorbed in survival biasing
score = p % absorb_wgt * flux
else
! Skip any event where the particle wasn't absorbed
if (p % event == EVENT_SCATTER) cycle SCORE_LOOP
! All fission and absorption events will contribute here, so we
! can just use the particle's weight entering the collision
score = p % last_wgt * flux
end if
else
if (i_nuclide > 0) then
score = micro_xs(i_nuclide) % absorption * atom_density * flux
else
score = material_xs % absorption * flux
end if
end if
case (SCORE_FISSION)
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! fission
if (micro_xs(p % event_nuclide) % absorption > ZERO) then
score = p % absorb_wgt * micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption * flux
else
score = ZERO
end if
else
! Skip any non-absorption events
if (p % event == EVENT_SCATTER) cycle SCORE_LOOP
! All fission events will contribute, so again we can use
! particle's weight entering the collision as the estimate for the
! fission reaction rate
score = p % last_wgt * micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption * flux
end if
else
if (i_nuclide > 0) then
score = micro_xs(i_nuclide) % fission * atom_density * flux
else
score = material_xs % fission * flux
end if
end if
case (SCORE_NU_FISSION)
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing .or. p % fission) then
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
! Normally, we only need to make contributions to one scoring
! bin. However, in the case of fission, since multiple fission
! neutrons were emitted with different energies, multiple
! outgoing energy bins may have been scored to. The following
! logic treats this special case and results to multiple bins
call score_fission_eout_ce(p, t, score_index, score_bin)
cycle SCORE_LOOP
end if
end if
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! nu-fission
if (micro_xs(p % event_nuclide) % absorption > ZERO) then
score = p % absorb_wgt * micro_xs(p % event_nuclide) % &
nu_fission / micro_xs(p % event_nuclide) % absorption * flux
else
score = ZERO
end if
else
! Skip any non-fission events
if (.not. p % fission) cycle SCORE_LOOP
! If there is no outgoing energy filter, than we only need to
! score to one bin. For the score to be 'analog', we need to
! score the number of particles that were banked in the fission
! bank. Since this was weighted by 1/keff, we multiply by keff
! to get the proper score.
score = keff * p % wgt_bank * flux
end if
else
if (i_nuclide > 0) then
score = micro_xs(i_nuclide) % nu_fission * atom_density * flux
else
score = material_xs % nu_fission * flux
end if
end if
case (SCORE_PROMPT_NU_FISSION)
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing .or. p % fission) then
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
! Normally, we only need to make contributions to one scoring
! bin. However, in the case of fission, since multiple fission
! neutrons were emitted with different energies, multiple
! outgoing energy bins may have been scored to. The following
! logic treats this special case and results to multiple bins
call score_fission_eout_ce(p, t, score_index, score_bin)
cycle SCORE_LOOP
end if
end if
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! prompt-nu-fission
if (micro_xs(p % event_nuclide) % absorption > ZERO) then
score = p % absorb_wgt * micro_xs(p % event_nuclide) % fission &
* nuclides(p % event_nuclide) % nu(E, EMISSION_PROMPT) &
/ micro_xs(p % event_nuclide) % absorption
else
score = ZERO
end if
else
! Skip any non-fission events
if (.not. p % fission) cycle SCORE_LOOP
! If there is no outgoing energy filter, than we only need to
! score to one bin. For the score to be 'analog', we need to
! score the number of particles that were banked in the fission
! bank as prompt neutrons. Since this was weighted by 1/keff, we
! multiply by keff to get the proper score.
score = keff * p % wgt_bank * (ONE - sum(p % n_delayed_bank) &
/ real(p % n_bank, 8))
end if
else
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
E = p % E
else
E = p % last_E
end if
if (i_nuclide > 0) then
score = micro_xs(i_nuclide) % fission * nuclides(i_nuclide) % &
nu(E, EMISSION_PROMPT) * atom_density * flux
else
score = ZERO
! Loop over all nuclides in the current material
do l = 1, materials(p % material) % n_nuclides
! Get atom density
atom_density_ = materials(p % material) % atom_density(l)
! Get index in nuclides array
i_nuc = materials(p % material) % nuclide(l)
! Accumulate the contribution from each nuclide
score = score + micro_xs(i_nuc) % fission * nuclides(i_nuc) % &
nu(E, EMISSION_PROMPT) * atom_density_ * flux
end do
end if
end if
case (SCORE_DELAYED_NU_FISSION)
! make sure the correct energy is used
if (t % estimator == ESTIMATOR_TRACKLENGTH) then
E = p % E
else
E = p % last_E
end if
! Set the delayedgroup filter index and the number of delayed group bins
dg_filter = t % find_filter(FILTER_DELAYEDGROUP)
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing .or. p % fission) then
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
! Normally, we only need to make contributions to one scoring
! bin. However, in the case of fission, since multiple fission
! neutrons were emitted with different energies, multiple
! outgoing energy bins may have been scored to. The following
! logic treats this special case and results to multiple bins
call score_fission_eout_ce(p, t, score_index, score_bin)
cycle SCORE_LOOP
end if
end if
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! delayed-nu-fission
if (micro_xs(p % event_nuclide) % absorption > ZERO) then
! Check if the delayed group filter is present
if (dg_filter > 0) then
select type(filt => t % filters(dg_filter) % obj)
type is (DelayedGroupFilter)
! Loop over all delayed group bins and tally to them
! individually
do d_bin = 1, filt % n_bins
! Get the delayed group for this bin
d = filt % groups(d_bin)
! Compute the yield for this delayed group
yield = nuclides(p % event_nuclide) &
% nu(E, EMISSION_DELAYED, d)
! Compute the score and tally to bin
score = p % absorb_wgt * yield &
* micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
cycle SCORE_LOOP
end select
else
! If the delayed group filter is not present, compute the score
! by multiplying the absorbed weight by the fraction of the
! delayed-nu-fission xs to the absorption xs
score = p % absorb_wgt * micro_xs(p % event_nuclide) % fission &
* nuclides(p % event_nuclide) % nu(E, EMISSION_DELAYED) &
/ micro_xs(p % event_nuclide) % absorption
end if
end if
else
! Skip any non-fission events
if (.not. p % fission) cycle SCORE_LOOP
! If there is no outgoing energy filter, than we only need to
! score to one bin. For the score to be 'analog', we need to
! score the number of particles that were banked in the fission
! bank. Since this was weighted by 1/keff, we multiply by keff
! to get the proper score. Loop over the neutrons produced from
! fission and check which ones are delayed. If a delayed neutron is
! encountered, add its contribution to the fission bank to the
! score.
! Check if the delayed group filter is present
if (dg_filter > 0) then
select type(filt => t % filters(dg_filter) % obj)
type is (DelayedGroupFilter)
! Loop over all delayed group bins and tally to them
! individually
do d_bin = 1, filt % n_bins
! Get the delayed group for this bin
d = filt % groups(d_bin)
! Compute the score and tally to bin
score = keff * p % wgt_bank / p % n_bank * p % n_delayed_bank(d)
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
cycle SCORE_LOOP
end select
else
! Add the contribution from all delayed groups
score = keff * p % wgt_bank / p % n_bank * sum(p % n_delayed_bank)
end if
end if
else
! Check if tally is on a single nuclide
if (i_nuclide > 0) then
! Check if the delayed group filter is present
if (dg_filter > 0) then
select type(filt => t % filters(dg_filter) % obj)
type is (DelayedGroupFilter)
! Loop over all delayed group bins and tally to them
! individually
do d_bin = 1, filt % n_bins
! Get the delayed group for this bin
d = filt % groups(d_bin)
! Compute the yield for this delayed group
yield = nuclides(i_nuclide) % nu(E, EMISSION_DELAYED, d)
! Compute the score and tally to bin
score = micro_xs(i_nuclide) % fission * yield * &
atom_density * flux
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
cycle SCORE_LOOP
end select
else
! If the delayed group filter is not present, compute the score
! by multiplying the delayed-nu-fission macro xs by the flux
score = micro_xs(i_nuclide) % fission * nuclides(i_nuclide) % &
nu(E, EMISSION_DELAYED) * atom_density * flux
end if
! Tally is on total nuclides
else
! Check if the delayed group filter is present
if (dg_filter > 0) then
select type(filt => t % filters(dg_filter) % obj)
type is (DelayedGroupFilter)
! Loop over all nuclides in the current material
do l = 1, materials(p % material) % n_nuclides
! Get atom density
atom_density_ = materials(p % material) % atom_density(l)
! Get index in nuclides array
i_nuc = materials(p % material) % nuclide(l)
! Loop over all delayed group bins and tally to them
! individually
do d_bin = 1, filt % n_bins
! Get the delayed group for this bin
d = filt % groups(d_bin)
! Get the yield for the desired nuclide and delayed group
yield = nuclides(i_nuc) % nu(E, EMISSION_DELAYED, d)
! Compute the score and tally to bin
score = micro_xs(i_nuc) % fission * yield * atom_density_ * flux
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
end do
cycle SCORE_LOOP
end select
else
score = ZERO
! Loop over all nuclides in the current material
do l = 1, materials(p % material) % n_nuclides
! Get atom density
atom_density_ = materials(p % material) % atom_density(l)
! Get index in nuclides array
i_nuc = materials(p % material) % nuclide(l)
! Accumulate the contribution from each nuclide
score = score + micro_xs(i_nuc) % fission * nuclides(i_nuc) % &
nu(E, EMISSION_DELAYED) * atom_density_ * flux
end do
end if
end if
end if
case (SCORE_DECAY_RATE)
! Set the delayedgroup filter index
dg_filter = t % find_filter(FILTER_DELAYEDGROUP)
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! delayed-nu-fission
if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
nuclides(p % event_nuclide) % fissionable) then
! Check if the delayed group filter is present
if (dg_filter > 0) then
select type(filt => t % filters(dg_filter) % obj)
type is (DelayedGroupFilter)
! Loop over all delayed group bins and tally to them
! individually
do d_bin = 1, filt % n_bins
! Get the delayed group for this bin
d = filt % groups(d_bin)
! Compute the yield for this delayed group
yield = nuclides(p % event_nuclide) &
% nu(E, EMISSION_DELAYED, d)
associate (rxn => nuclides(p % event_nuclide) % &
reactions(nuclides(p % event_nuclide) % index_fission(1)))
! Compute the score
score = p % absorb_wgt * yield * &
micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption &
* rxn % products(1 + d) % decay_rate
end associate
! Tally to bin
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
cycle SCORE_LOOP
end select
else
! If the delayed group filter is not present, compute the score
! by accumulating the absorbed weight times the decay rate times
! the fraction of the delayed-nu-fission xs to the absorption xs
! for all delayed groups.
score = ZERO
associate (rxn => nuclides(p % event_nuclide) % &
reactions(nuclides(p % event_nuclide) % index_fission(1)))
! We need to be careful not to overshoot the number of delayed
! groups since this could cause the range of the rxn % products
! array to be exceeded. Hence, we use the size of this array
! and not the MAX_DELAYED_GROUPS constant for this loop.
do d = 1, size(rxn % products) - 2
score = score + rxn % products(1 + d) % decay_rate * &
p % absorb_wgt * micro_xs(p % event_nuclide) % fission *&
nuclides(p % event_nuclide) % nu(E, EMISSION_DELAYED, d)&
/ micro_xs(p % event_nuclide) % absorption
end do
end associate
end if
end if
else
! Skip any non-fission events
if (.not. p % fission) cycle SCORE_LOOP
! If there is no outgoing energy filter, than we only need to
! score to one bin. For the score to be 'analog', we need to
! score the number of particles that were banked in the fission
! bank. Since this was weighted by 1/keff, we multiply by keff
! to get the proper score. Loop over the neutrons produced from
! fission and check which ones are delayed. If a delayed neutron is
! encountered, add its contribution to the fission bank to the
! score.
score = ZERO
! loop over number of particles banked
do k = 1, p % n_bank
! get the delayed group
g = fission_bank(n_bank - p % n_bank + k) % delayed_group
! Case for tallying delayed emissions
if (g /= 0) then
! Accumulate the decay rate times delayed nu fission score
associate (rxn => nuclides(p % event_nuclide) % &
reactions(nuclides(p % event_nuclide) % index_fission(1)))
! determine score based on bank site weight and keff.
score = score + keff * fission_bank(n_bank - p % n_bank + k) &
% wgt * rxn % products(1 + g) % decay_rate
end associate
! if the delayed group filter is present, tally to corresponding
! delayed group bin if it exists
if (dg_filter > 0) then
! declare the delayed group filter type
select type(filt => t % filters(dg_filter) % obj)
type is (DelayedGroupFilter)
! loop over delayed group bins until the corresponding bin is
! found
do d_bin = 1, filt % n_bins
d = filt % groups(d_bin)
! check whether the delayed group of the particle is equal to
! the delayed group of this bin
if (d == g) then
call score_fission_delayed_dg(t, d_bin, score, score_index)
end if
end do
end select
! Reset the score to zero
score = ZERO
end if
end if
end do
! If the delayed group filter is present, cycle because the
! score_fission_delayed_dg(...) has already tallied the score
if (dg_filter > 0) then
cycle SCORE_LOOP
end if
end if
case (SCORE_KAPPA_FISSION)
! Determine kappa-fission cross section on the fly. The ENDF standard
! (ENDF-102) states that MT 18 stores the fission energy as the Q_value
! (fission(1))
score = ZERO
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! fission scaled by kappa-fission
associate (nuc => nuclides(p % event_nuclide))
if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
nuc % fissionable) then
score = p % absorb_wgt * &
nuc % reactions(nuc % index_fission(1)) % Q_value * &
micro_xs(p % event_nuclide) % fission / &
micro_xs(p % event_nuclide) % absorption * flux
end if
end associate
else
! Skip any non-absorption events
if (p % event == EVENT_SCATTER) cycle SCORE_LOOP
! All fission events will contribute, so again we can use
! particle's weight entering the collision as the estimate for
! the fission energy production rate
associate (nuc => nuclides(p % event_nuclide))
if (nuc % fissionable) then
score = p % last_wgt * &
nuc % reactions(nuc % index_fission(1)) % Q_value * &
micro_xs(p % event_nuclide) % fission / &
micro_xs(p % event_nuclide) % absorption * flux
end if
end associate
end if
else
if (i_nuclide > 0) then
associate (nuc => nuclides(i_nuclide))
if (nuc % fissionable) then
score = nuc % reactions(nuc % index_fission(1)) % Q_value * &
micro_xs(i_nuclide) % fission * atom_density * flux
end if
end associate
else
do l = 1, materials(p % material) % n_nuclides
! Determine atom density and index of nuclide
atom_density_ = materials(p % material) % atom_density(l)
i_nuc = materials(p % material) % nuclide(l)
! If nuclide is fissionable, accumulate kappa fission
associate(nuc => nuclides(i_nuc))
if (nuc % fissionable) then
score = score + &
nuc % reactions(nuc % index_fission(1)) % Q_value * &
micro_xs(i_nuc) % fission * atom_density_ * flux
end if
end associate
end do
end if
end if
case (SCORE_EVENTS)
! Simply count number of scoring events
score = ONE
case (ELASTIC)
if (t % estimator == ESTIMATOR_ANALOG) then
! Check if event MT matches
if (p % event_MT /= ELASTIC) cycle SCORE_LOOP
score = p % last_wgt * flux
else
if (i_nuclide > 0) then
score = micro_xs(i_nuclide) % elastic * atom_density * flux
else
score = material_xs % elastic * flux
end if
end if
case (SCORE_FISS_Q_PROMPT)
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! fission scaled by Q-value
associate (nuc => nuclides(p % event_nuclide))
if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
allocated(nuc % fission_q_prompt)) then
score = p % absorb_wgt &
* nuc % fission_q_prompt % evaluate(p % last_E) &
* micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption * flux
end if
end associate
else
! Skip any non-absorption events
if (p % event == EVENT_SCATTER) cycle SCORE_LOOP
! All fission events will contribute, so again we can use
! particle's weight entering the collision as the estimate for
! the fission energy production rate
associate (nuc => nuclides(p % event_nuclide))
if (allocated(nuc % fission_q_prompt)) then
score = p % last_wgt &
* nuc % fission_q_prompt % evaluate(p % last_E) &
* micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption * flux
end if
end associate
end if
else
if (t % estimator == ESTIMATOR_COLLISION) then
E = p % last_E
else
E = p % E
end if
if (i_nuclide > 0) then
if (allocated(nuclides(i_nuclide) % fission_q_prompt)) then
score = micro_xs(i_nuclide) % fission * atom_density * flux &
* nuclides(i_nuclide) % fission_q_prompt % evaluate(E)
else
score = ZERO
end if
else
score = ZERO
do l = 1, materials(p % material) % n_nuclides
atom_density_ = materials(p % material) % atom_density(l)
i_nuc = materials(p % material) % nuclide(l)
if (allocated(nuclides(i_nuc) % fission_q_prompt)) then
score = score + micro_xs(i_nuc) % fission * atom_density_ &
* flux &
* nuclides(i_nuc) % fission_q_prompt % evaluate(E)
end if
end do
end if
end if
case (SCORE_FISS_Q_RECOV)
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! fission scaled by Q-value
associate (nuc => nuclides(p % event_nuclide))
if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
allocated(nuc % fission_q_recov)) then
score = p % absorb_wgt &
* nuc % fission_q_recov % evaluate(p % last_E) &
* micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption * flux
end if
end associate
else
! Skip any non-absorption events
if (p % event == EVENT_SCATTER) cycle SCORE_LOOP
! All fission events will contribute, so again we can use
! particle's weight entering the collision as the estimate for
! the fission energy production rate
associate (nuc => nuclides(p % event_nuclide))
if (allocated(nuc % fission_q_recov)) then
score = p % last_wgt &
* nuc % fission_q_recov % evaluate(p % last_E) &
* micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption * flux
end if
end associate
end if
else
if (t % estimator == ESTIMATOR_COLLISION) then
E = p % last_E
else
E = p % E
end if
if (i_nuclide > 0) then
if (allocated(nuclides(i_nuclide) % fission_q_recov)) then
score = micro_xs(i_nuclide) % fission * atom_density * flux &
* nuclides(i_nuclide) % fission_q_recov % evaluate(E)
else
score = ZERO
end if
else
score = ZERO
do l = 1, materials(p % material) % n_nuclides
atom_density_ = materials(p % material) % atom_density(l)