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
|
program monte_carlo
use, intrinsic :: iso_fortran_env
implicit none
!> Initialize some variables for neutrons
integer(int64), parameter :: n = 1e7
integer(int64), parameter :: bins = 800
real(real64), parameter :: eps = 1.0e-12_real64
!> And for geometry
real(real64), parameter :: L = 8.0_real64
real(real64), parameter :: dx = L / bins
real(real64), dimension(6), parameter :: boundaries = [0.0_real64, 2.0_real64, 3.0_real64, 5.0_real64, 6.0_real64, 8.0_real64]
!> These will contain the x-position and cosine of angle for each neutron
real(real64), dimension(n) :: px, mu
!> Track whether or not the neutron is alive
logical, dimension(n) :: is_alive = .true.
!> And for tallying the flux
real(real64), dimension(bins) :: flux_tally = 0.0_real64
real(real64), dimension(bins) :: thread_tally = 0.0_real64
!> Some local trackers to be updated via subroutine call
real(real64) :: sigma_t, sigma_s, q, d_flight, d_boundary, d_step, xi, next_boundary
integer(int64) :: i, region
!> Initialize RNG to be actually random
call random_seed()
call random_init(.false., .false.)
!> Provide all neutrons an initial location/angle
call initialize_neutrons(n, px, mu)
!> The main monte carlo loop
main: do i = 1,n
!> Skip if the neutron is already lost
if(.not. is_alive(i)) cycle main
!> Otherwise, simulate it until it is lost
single: do while(is_alive(i))
!> Find what region the neutron is in and get XS
region = count(px(i) .ge. boundaries)
call get_region_info(region, sigma_t, sigma_s, q)
!> Determine the distance before collision
if(sigma_t .eq. 0.0_real64) then
!> If we are in a void, make it huge
d_flight = huge(real64)
else
!> Else, calculate with d = -ln(xi)/sigma_t
call random_number(xi)
d_flight = -log(xi) / sigma_t
end if
!> Determine where the next boundary is and the distance
if(mu(i) .gt. 0.0_real64) then
d_boundary = (boundaries(region + 1) - px(i)) / mu(i)
else if(mu(i) .lt. 0.0_real64) then
d_boundary = (boundaries(region) - px(i)) / mu(i)
else
!> In the rare chance it is perfectly parallel with the walls
d_boundary = huge(real64)
is_alive(i) = .false.
end if
!> How far we will actually go
d_step = min(d_flight,d_boundary)
!> Tally the track-length in the bin
call tally_segment(px(i), mu(i), d_step, bins, dx, thread_tally)
!> Update the neutron position
px(i) = px(i) + d_step * mu(i)
!> If the neutron crosses a boundary, restart the simulation at the new location
if(d_boundary .le. d_flight) then
!> We cross a boundary
if(mu(i) .gt. 0.0_real64 .and. boundaries(region+1) .eq. L ) then
!> Bin neutron at vacuum boundary at x=8.0
is_alive(i) = .false.
cycle main
else if(mu(i) .lt. 0.0_real64 .and. boundaries(region) .eq. 0.0_real64) then
!> Reflect on specular boundary at x=0.0
mu(i) = -mu(i)
px(i) = eps
cycle single
else
!> Internal boundary
if(mu(i) .gt. 0.0_real64) then
px(i) = boundaries(region+1) + eps
else
px(i) = boundaries(region) - eps
end if
cycle single
end if
end if
!> If not, determine what type of interaction happens
call random_number(xi)
if(xi .le. sigma_s / sigma_t) then
!> Scattering interaction
!> Sample new isotropic angle
call random_number(xi)
mu(i) = 2.0_real64 * xi - 1.0_real64
else
!> Absorption interaction
is_alive(i) = .false.
end if
end do single
end do main
flux_tally = flux_tally + thread_tally
flux_tally = flux_tally / thread_tally(100)
open(unit=10, file='flux.out', status='replace')
write: do i = 1,bins
write(10,'(2ES15.7)') (i-0.5_real64)*dx, flux_tally(i)
end do write
contains
pure elemental subroutine get_region_info(region, sigma_t, sigma_s, q)
integer(int64), intent(in) :: region
real(real64), intent(out) :: sigma_t, sigma_s, q
if(region .eq. 1) then
sigma_t = 50.0_real64; sigma_s = 0.0_real64; q = 50.0_real64
else if(region .eq. 2) then
sigma_t = 5.0_real64; sigma_s = 0.0_real64; q = 5.0_real64
else if(region .eq. 3) then
sigma_t = 0.0_real64; sigma_s = 0.0_real64; q = 0.0_real64
else if(region .eq. 4) then
sigma_t = 1.0_real64; sigma_s = 0.9_real64; q = 1.0_real64
else if(region .eq. 5) then
sigma_t = 1.0_real64; sigma_s = 0.9_real64; q = 0.0_real64
end if
end subroutine get_region_info
impure subroutine initialize_neutrons(n, px, mu)
integer(int64), intent(in) :: n
real(real64), intent(out), dimension(n) :: px, mu
real(real64), dimension(n) :: r_region, r_position, r_mu
!> Make some uniform random numbers
call random_number(r_region)
call random_number(r_position)
call random_number(r_mu)
!> Consider the total integrated source-distance for regions
r_region = r_region * 106.0_real64
!> Assign a random position in each source region
where(r_region .lt. 100)
px = r_position * 2.0_real64
else where(r_region .lt. 105)
px = 2.0_real64 + r_position
else where
px = 5.0_real64 + r_position
end where
!> Assign random angles to each particle
mu = 2.0_real64 * r_mu - 1.0_real64
end subroutine initialize_neutrons
pure subroutine tally_segment(px, mu, d_step, bins, dx, thread_tally)
integer(int64), intent(in) :: bins
real(real64), intent(in) :: px, mu, d_step, dx
real(real64), intent(in out), dimension(bins) :: thread_tally
real(real64) :: x_start, x_end, x_left, x_right
integer(int64) :: i_start, i_end, i
if (d_step <= 0.0_real64) return
x_start = px
x_end = px + mu * d_step
if (x_start > x_end) call swap(x_start, x_end)
i_start = floor(x_start / dx) + 1
i_end = floor(x_end / dx) + 1
i_start = max(1, min(bins, i_start))
i_end = max(1, min(bins, i_end))
if (i_start == i_end) then
thread_tally(i_start) = thread_tally(i_start) + d_step
else
! First partial bin
x_right = i_start * dx
thread_tally(i_start) = thread_tally(i_start) + (x_right - x_start)
! Full interior bins
do i = i_start+1, i_end-1
thread_tally(i) = thread_tally(i) + dx
end do
! Last partial bin
x_left = (i_end - 1) * dx
thread_tally(i_end) = thread_tally(i_end) + (x_end - x_left)
end if
end subroutine tally_segment
pure elemental subroutine swap(x,y)
real(real64), intent(in out) :: x,y
real(real64) :: temp
temp = x
x = y
y = temp
end subroutine swap
end program monte_carlo
|
