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