diff options
| -rw-r--r-- | .gitignore | 3 | ||||
| -rw-r--r-- | Makefile | 9 | ||||
| -rw-r--r-- | main.f90 | 147 | ||||
| -rw-r--r-- | plot.gnu | 3 |
4 files changed, 127 insertions, 35 deletions
diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..9482666 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +*.x +*.out +*.mod @@ -1,10 +1,13 @@ FC = gfortran -FFLAGS = -march=native -Ofast +FFLAGS = -march=native -O3 XNAME = monte_carlo.x NPROC := $(shell nproc) -all: build run plot +all: clean build run plot + +clean: + rm -f *.x *.out *.mod build: main.f90 $(FC) $(FFLAGS) main.f90 -o $(XNAME) @@ -13,4 +16,4 @@ run: build OMP_NUM_THREADS=$(NPROC) ./$(XNAME) plot: run - echo test + gnuplot plot.gnu @@ -4,8 +4,9 @@ program monte_carlo implicit none !> Initialize some variables for neutrons - integer(int64), parameter :: n = 1e6 + 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 @@ -15,66 +16,95 @@ program monte_carlo !> 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 + 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, d_boundary, xi, next_boundary + 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 - single: do while (is_alive(i)) + !> 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) + 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(sigma_t .eq. 0.0_real64) then !> If we are in a void, make it huge - d = 1e10_real64 + d_flight = huge(real64) else !> Else, calculate with d = -ln(xi)/sigma_t call random_number(xi) - d = -log(xi) / sigma_t + 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) + 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 = 1e10_real64 + d_boundary = huge(real64) is_alive(i) = .false. end if - !> If it goes through a wall, restart the particle at the new boundary - if(d_boundary.le.d) then - px(i) = px(i) + d_boundary*mu(i) - if(px(i).le.0.0_real64) then - px(i) = 0.0_real64 + !> 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) - else if (px(i).ge.L) then - is_alive = .false. - exit single + 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 - cycle single end if !> If not, determine what type of interaction happens call random_number(xi) - if(xi.le.sigma_s/sigma_t) then + 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 + mu(i) = 2.0_real64 * xi - 1.0_real64 else !> Absorption interaction is_alive(i) = .false. @@ -82,9 +112,16 @@ program monte_carlo 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 @@ -92,15 +129,15 @@ contains integer(int64), intent(in) :: region real(real64), intent(out) :: sigma_t, sigma_s, q - if(region.eq.1) then + 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 + 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 + 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 + 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 + 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 @@ -119,9 +156,9 @@ contains r_region = r_region * 106.0_real64 !> Assign a random position in each source region - where(r_region.lt.100) + where(r_region .lt. 100) px = r_position * 2.0_real64 - else where(r_region.lt.105) + else where(r_region .lt. 105) px = 2.0_real64 + r_position else where px = 5.0_real64 + r_position @@ -132,4 +169,50 @@ contains 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 diff --git a/plot.gnu b/plot.gnu new file mode 100644 index 0000000..46079c6 --- /dev/null +++ b/plot.gnu @@ -0,0 +1,3 @@ +set term x11 persist + +plot 'flux.out' using 1:2 with lines |