diff options
| -rw-r--r-- | Makefile | 2 | ||||
| -rw-r--r-- | main.f90 | 235 | ||||
| -rw-r--r-- | plot.gnu | 3 |
3 files changed, 80 insertions, 160 deletions
@@ -1,5 +1,5 @@ FC = gfortran -FFLAGS = -march=native -O3 -fopenmp -fcheck=all +FFLAGS = -march=native -O3 XNAME = monte_carlo.x NPROC := $(shell nproc) @@ -5,61 +5,53 @@ program monte_carlo implicit none !> Initialize some variables for neutrons - integer(int64), parameter :: n = 1e5 + integer(int64), parameter :: n = 1e8 integer(int64), parameter :: bins = 800 - real(real64), parameter :: eps = 1.0e-12_real64 + real(real64), parameter :: eps = 1.0e-6_real64 + real(real64), parameter :: pi = 4.0_real64 * atan(1.0_real64) !> And for geometry real(real64), parameter :: L = 8.0_real64 - real(real64), parameter :: dx = L / bins + real(real64), parameter :: bin_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 position and cosines of direction vector - real(real64), allocatable, dimension(:) :: x, y, z, u, v, w + !> These will contain the position and cosine of the neutron direction vector + real(real64), allocatable, dimension(:) :: x, u !> Track whether or not the neutron is alive logical, allocatable, dimension(:) :: is_alive !> And for tallying the flux - real(real64), allocatable, dimension(:) :: flux_tally, collision_tally - real(real64), allocatable, dimension(:) :: thread_tally, col_thread_tally + real(real64), allocatable, dimension(:) :: flux_tally !> 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 + real(real64) :: sigma_t, sigma_s, d_flight, d_boundary, d_step, xi, avg_r1_flux integer(int64) :: i, region !> Allocate arrays - allocate(x(n), y(n), z(n), u(n), v(n), w(n), is_alive(n), flux_tally(n), collision_tally(n), thread_tally(n), col_thread_tally(n)) + allocate(x(n), u(n), is_alive(n), flux_tally(bins)) is_alive = .true. flux_tally = 0.0_real64 - collision_tally = 0.0_real64 - thread_tally = 0.0_real64 - col_thread_tally = 0.0_real64 !> Initialize RNG to be actually random call random_seed() call random_init(.false., .true.) !> Provide all neutrons an initial location/angle - call initialize_neutrons(n, x, y, z, u, v, w) - - !$OMP PARALLEL PRIVATE(thread_tally, col_thread_tally, sigma_t, sigma_s, q, d_flight, d_boundary, d_step, xi, region, i) - thread_tally=0.0_real64 + call initialize_neutrons(n, x, u) !> The main monte carlo loop - !$OMP DO main: do i = 1,n - !write(*,'("n = ",i0)') i - !write(*,'("Before: x=",E9.3," y=",E9.3," z=",E9.3)') x(i), y(i), z(i) - !write(*,'("Before: u=",E9.3," v=",E9.3," w=",E9.3)') u(i), v(i), w(i) - - !> Skip if the neutron is already lost - if(.not. is_alive(i)) cycle main - !> Otherwise, simulate it until it is lost + !> Simulate the neutron until it is lost single: do while(is_alive(i)) !> Find what region the neutron is in and get XS region = count(x(i) .ge. boundaries) - call get_region_info(region, sigma_t, sigma_s, q) + call get_region_info(region, sigma_t, sigma_s) + + if(region .ge. size(boundaries) .or. region .le. 0) then + !> We are outside of the problem. Not good. Skip. + u(i) = eps + end if !> Determine the distance before collision if(sigma_t .eq. 0.0_real64) then @@ -73,8 +65,10 @@ program monte_carlo !> Determine where the next boundary is and the distance if(u(i) .gt. 0.0_real64) then + !> Moving in the +x direction d_boundary = (boundaries(region + 1) - x(i)) / u(i) else if(u(i) .lt. 0.0_real64) then + !> Moving in the -x direction d_boundary = (boundaries(region) - x(i)) / u(i) else !> In the rare chance it is perfectly parallel with the walls @@ -83,183 +77,145 @@ program monte_carlo !> How far we will actually go d_step = min(d_flight,d_boundary) - !write(*,'("Step distance = ",E9.3)') d_step !> Tally the track-length contributions for the bin - call tally_segment(x(i), u(i), d_step, bins, dx, thread_tally) + call tally_segment(x(i), u(i), d_step, bins, bin_dx, flux_tally) !> Update the neutron position x(i) = x(i) + d_step * u(i) - y(i) = y(i) + d_step * v(i) - z(i) = z(i) + d_step * w(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(abs(x(i) - L) .lt. eps) then - !write(*,'(A)') "Escaped!" + if(u(i) .gt. 0.0_real64 .and. abs(x(i)-L) .lt. eps) then !> Bin neutron at vacuum boundary at x=8.0 is_alive(i) = .false. cycle main - else if(x(i) .lt. eps) then - !write(*,'(A)') "Reflected!" + else if(u(i) .lt. 0.0_real64 .and. x(i) .lt. eps) then !> Reflect on specular boundary at x=0.0 u(i) = -u(i) - v(i) = -v(i) - w(i) = -w(i) - x(i) = eps*0.1_real64 - cycle single + x(i) = -x(i) else - !write(*,'(A)') "Internal boundary!" !> Internal boundary if(u(i) .gt. 0.0_real64) then - !write(*,'(A)') "Forwards reflection!" - x(i) = boundaries(region+1) + 2*eps + x(i) = boundaries(region+1) + 2.0_real64*eps else - !write(*,'(A)') "Backwards reflection!" - x(i) = boundaries(region) - 2*eps + x(i) = boundaries(region) - 2.0_real64*eps end if - cycle single end if - end if - - !> If not, a collision happens. Tally it: - !> The current 'bin' is ceiling(x/dx) - col_thread_tally(max(ceiling(x/dx),0)) = col_thread_tally(ceiling(x/dx)) + 1 - - !> And determine the type: - call random_number(xi) - if(xi .le. sigma_s / sigma_t) then - !> Scattering interaction - !> Sample new isotropic angle - call scatter(u(i),v(i),w(i)) else - !> Absorption interaction - is_alive(i) = .false. + !> A collision happens + call random_number(xi) + if(xi .le. sigma_s/sigma_t) then + !> This is a scattering interaction + call random_number(xi) + u(i) = 2.0_real64 * xi - 1.0_real64 + else + !> The particle is absorbed and lost + is_alive(i) = .false. + cycle main + end if end if - !write(*,'("x=",E9.3," y=",E9.3," z=",E9.3)') x(i), y(i), z(i) - !write(*,'("u=",E9.3," v=",E9.3," w=",E9.3)') u(i), v(i), w(i) end do single end do main - !$OMP END DO - !$OMP CRITICAL - flux_tally = flux_tally + thread_tally - collision_tally = collision_tally + thread_tally - !$OMP END CRITICAL - !$OMP END PARALLEL - !flux_tally = flux_tally / flux_tally(bins/8) - flux_tally = flux_tally / (real(n, real64) * dx) - collision_tally = collision_tally / (real(n, real64)*4000) + !> Average region 1 flux + avg_r1_flux = sum(flux_tally(1:floor(boundaries(2)/bin_dx)))/real(floor(boundaries(2)/bin_dx), real64) + + !> Normalize to region 1 + flux_tally = flux_tally / avg_r1_flux open(unit=10, file='flux.out', status='replace') - open(unit=11, file='collision.out', status='replace') write: do i = 1,bins - write(10,'(2E15.7)') (i-0.5_real64)*dx, flux_tally(i) - write(11,'(2E15.7)') (i-0.5_real64)*dx, collision_tally(i) + write(10,'(2E15.7)') (i-0.5_real64)*bin_dx, flux_tally(i) end do write !> Deallocate everyting - deallocate(x, y, z, u, v, w, is_alive, flux_tally, collision_tally, thread_tally, col_thread_tally) + deallocate(x, u, is_alive, flux_tally) contains - pure elemental subroutine get_region_info(region, sigma_t, sigma_s, q) + pure elemental subroutine get_region_info(region, sigma_t, sigma_s) integer(int64), intent(in) :: region - real(real64), intent(out) :: sigma_t, sigma_s, q + real(real64), intent(out) :: sigma_t, sigma_s if(region .eq. 1) then - sigma_t = 50.0_real64; sigma_s = 0.0_real64; q = 50.0_real64 + sigma_t = 50.0_real64; sigma_s = 0.0_real64 else if(region .eq. 2) then - sigma_t = 5.0_real64; sigma_s = 0.0_real64; q = 5.0_real64 + sigma_t = 5.0_real64; sigma_s = 0.0_real64 else if(region .eq. 3) then - sigma_t = 0.0_real64; sigma_s = 0.0_real64; q = 0.0_real64 + sigma_t = 0.0_real64; sigma_s = 0.0_real64 else if(region .eq. 4) then - sigma_t = 1.0_real64; sigma_s = 0.9_real64; q = 1.0_real64 + sigma_t = 1.0_real64; sigma_s = 0.9_real64 else if(region .eq. 5) then - sigma_t = 1.0_real64; sigma_s = 0.9_real64; q = 0.0_real64 + sigma_t = 1.0_real64; sigma_s = 0.9_real64 end if end subroutine get_region_info - impure subroutine initialize_neutrons(n, x, y, z, u, v, w) + impure subroutine initialize_neutrons(n, x, u) integer(int64), intent(in) :: n - real(real64), intent(out), dimension(n) :: x, y, z, u, v, w - real(real64), dimension(n) :: r_region, r_position, r_mu, r_phi + real(real64), intent(out), dimension(n) :: x, u + real(real64), dimension(n) :: r_region, r_position, r_mu real(real64), parameter :: pi = 4.0_real64 * atan(1.0_real64) !> Make some uniform random numbers call random_number(r_region) call random_number(r_position) call random_number(r_mu) - call random_number(r_phi) !> Consider the total integrated source-distance for regions - r_region = r_region * 106.0_real64 + r_region = r_region * 101.0_real64 !> Assign a random position in each source region where(r_region .lt. 100.0_real64) x = eps + r_position * (2.0_real64 - 2.0_real64*eps) - else where(r_region .lt. 105.0_real64) - x = 2.0_real64 + eps + r_position * (1.0_real64 - 2.0_real64*eps) - else where(r_region .lt. 106.0_real64) + else where(r_region .lt. 101.0_real64) x = 5.0_real64 + eps + r_position * (1.0_real64 - 2.0_real64*eps) end where - !> Start on the x-axis - y = 0.0_real64 - z = 0.0_real64 - !> Assign random angles to each particle - u = sqrt(max(1.0_real64 - r_mu**2.0_real64, 0.0_real64))*cos(2.0_real64*pi*r_phi) - v = sqrt(max(1.0_real64 - r_mu**2.0_real64, 0.0_real64))*sin(2.0_real64*pi*r_phi) - w = r_mu + u = 2.0_real64 * r_mu - 1.0_real64 end subroutine initialize_neutrons - pure subroutine tally_segment(px, mu, d_step, bins, dx, thread_tally) + impure subroutine tally_segment(x, u, d_step, bins, bin_dx, flux_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 + real(real64), intent(in) :: x, u, d_step, bin_dx + real(real64), intent(in out), dimension(bins) :: flux_tally + real(real64) :: x_start, x_end, dx integer(int64) :: i_start, i_end, i - !> If no distance was travelled - if (d_step <= 0.0_real64) return - - !> Start and end position - x_start = px - x_end = px + mu * d_step - - !> Swap them if we're going backwards - if (x_start > x_end) call swap(x_start, x_end) + x_start = x + x_end = x + d_step*u + if (x_start .gt. x_end) then + call swap(x_start, x_end) + end if - !> Start and end indices for bins - i_start = ceiling(x_start / dx) - i_end = ceiling(x_end / dx) + i_start = ceiling(x_start/bin_dx) + i_end = ceiling(x_end/bin_dx) - !> Clamp them too + !> Clamp to valid bin range i_start = max(1, min(bins, i_start)) i_end = max(1, min(bins, i_end)) - if (i_start == i_end) then - !> No bins were crossed. Add the length - thread_tally(i_start) = thread_tally(i_start) + d_step + if (i_start .eq. i_end) then + flux_tally(i_start) = flux_tally(i_start) + d_step else - !> First partial bin exited - x_right = i_start * dx - thread_tally(i_start) = thread_tally(i_start) + (x_right - x_start)/abs(mu) + !> First partial bin + dx = (i_start*bin_dx - x_start)/abs(u) + flux_tally(i_start) = flux_tally(i_start) + dx - !> Going across the interior bins + !> Full bins in between do i = i_start+1, i_end-1 - thread_tally(i) = thread_tally(i) + dx/abs(mu) + flux_tally(i) = flux_tally(i) + bin_dx/abs(u) end do - !> Last partial bin entered - x_left = (i_end - 1) * dx - thread_tally(i_end) = thread_tally(i_end) + (x_end - x_left)/abs(mu) + !> Last partial bin + dx = (x_end - (i_end-1)*bin_dx)/abs(u) + flux_tally(i_end) = flux_tally(i_end) + dx end if end subroutine tally_segment @@ -272,39 +228,4 @@ contains y = temp end subroutine swap - impure elemental subroutine scatter(u, v, w) - real(real64), intent(in out) :: u, v, w - real(real64) :: cos_theta, phi, sin_theta, cos_phi, sin_phi - real(real64) :: up, vp, wp - real(real64) :: sqrt_1_minus_w2 - real(real64) :: uprime, vprime, wprime - real(real64) :: xi - real(real64), parameter :: pi = 4.0_real64 * atan(1.0_real64) - - !> Generate a random mu - call random_number(xi) - cos_theta = 2.0_real64 * xi - 1.0_real64 - - !> And a random phi - call random_number(xi) - phi = 2.0_real64 * pi * xi - - !> Introduce some quantities - sin_theta = sqrt(1.0_real64 - cos_theta**2) - cos_phi = cos(phi) - sin_phi = sin(phi) - sqrt_1_minus_w2 = sqrt(1.0_real64 - w**2) - - !> Lewis and Miller, Eq. 7-163 - up = sin_theta / sqrt_1_minus_w2 * (v*sin_phi - v*u*cos_phi) + u*cos_theta - vp = sin_theta / sqrt_1_minus_w2 * (-u*sin_phi - w*v*cos_phi) + v*cos_theta - wp = sin_theta * sqrt_1_minus_w2 * cos_phi + w*cos_theta - - !> Update new direction cosines - u = up - v = vp - w = wp - - end subroutine scatter - end program monte_carlo @@ -1,4 +1,3 @@ set term x11 persist -plot 'flux.out' using 1:2 with lines, \ - 'collision.out' using 1:2 with lines +plot 'flux.out' using 1:2 with lines |