diff options
| -rw-r--r-- | .gitignore | 2 | ||||
| -rw-r--r-- | Makefile | 3 | ||||
| -rw-r--r-- | main.f90 | 55 | ||||
| -rw-r--r-- | plot.gnu | 13 |
4 files changed, 46 insertions, 27 deletions
@@ -1,3 +1,5 @@ *.x *.out *.mod +*.tex +*.pdf @@ -1,7 +1,6 @@ FC = gfortran FFLAGS = -march=native -O3 XNAME = monte_carlo.x -NPROC := $(shell nproc) all: clean build run plot @@ -13,7 +12,7 @@ build: main.f90 $(FC) $(FFLAGS) main.f90 -o $(XNAME) run: build - OMP_NUM_THREADS=$(NPROC) ./$(XNAME) + ./$(XNAME) plot: run gnuplot plot.gnu @@ -1,11 +1,10 @@ program monte_carlo - use omp_lib use, intrinsic :: iso_fortran_env implicit none !> Initialize some variables for neutrons - integer(int64), parameter :: n = 1e8 + integer(int64), parameter :: n = 1e7 integer(int64), parameter :: bins = 800 real(real64), parameter :: eps = 1.0e-6_real64 real(real64), parameter :: pi = 4.0_real64 * atan(1.0_real64) @@ -20,16 +19,18 @@ program monte_carlo !> Track whether or not the neutron is alive logical, allocatable, dimension(:) :: is_alive !> And for tallying the flux - real(real64), allocatable, dimension(:) :: flux_tally + real(real64), allocatable, dimension(:) :: neutron_bin_length, flux_tally, flux_sq_tally, flux_std_dev !> Some local trackers to be updated via subroutine call 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), u(n), is_alive(n), flux_tally(bins)) + allocate(x(n), u(n), is_alive(n), neutron_bin_length(bins), flux_tally(bins), flux_sq_tally(bins), flux_std_dev(bins)) is_alive = .true. flux_tally = 0.0_real64 + flux_sq_tally = 0.0_real64 + flux_std_dev = 0.0_real64 !> Initialize RNG to be actually random call random_seed() @@ -40,6 +41,8 @@ program monte_carlo !> The main monte carlo loop main: do i = 1,n + !> Initialize empty neutron length + neutron_bin_length = 0.0_real64 !> Simulate the neutron until it is lost single: do while(is_alive(i)) @@ -49,7 +52,7 @@ program monte_carlo 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. + !> Put neutrons back inside if they get lost u(i) = eps end if @@ -79,7 +82,7 @@ program monte_carlo d_step = min(d_flight,d_boundary) !> Tally the track-length contributions for the bin - call tally_segment(x(i), u(i), d_step, bins, bin_dx, flux_tally) + call tally_segment(x(i), u(i), d_step, bins, bin_dx, neutron_bin_length) !> Update the neutron position x(i) = x(i) + d_step * u(i) @@ -90,7 +93,6 @@ program monte_carlo 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(u(i) .lt. 0.0_real64 .and. x(i) .lt. eps) then !> Reflect on specular boundary at x=0.0 u(i) = -u(i) @@ -113,23 +115,30 @@ program monte_carlo else !> The particle is absorbed and lost is_alive(i) = .false. - cycle main end if end if end do single + !> Update tallies + flux_tally = flux_tally + neutron_bin_length + flux_sq_tally = flux_sq_tally + neutron_bin_length**2 + end do main + !> sqrt(x²/n - (x/n)²)/(dx²) + flux_std_dev = sqrt((flux_sq_tally/real(n, real64) - (flux_tally/real(n, real64))**2)) / bin_dx**2 + !> 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 + flux_std_dev = flux_std_dev / avg_r1_flux open(unit=10, file='flux.out', status='replace') write: do i = 1,bins - write(10,'(2E15.7)') (i-0.5_real64)*bin_dx, flux_tally(i) + write(10,'(3E15.7)') (i-0.5_real64)*bin_dx, flux_tally(i), flux_std_dev(i) end do write !> Deallocate everyting @@ -181,12 +190,12 @@ contains end subroutine initialize_neutrons - impure subroutine tally_segment(x, u, d_step, bins, bin_dx, flux_tally) + pure subroutine tally_segment(x, u, d_step, bins, bin_dx, flux_tally) integer(int64), intent(in) :: bins 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 + integer(int64) :: start_bin, end_bin, bin x_start = x x_end = x + d_step*u @@ -194,28 +203,28 @@ contains call swap(x_start, x_end) end if - i_start = ceiling(x_start/bin_dx) - i_end = ceiling(x_end/bin_dx) + start_bin = ceiling(x_start/bin_dx) + end_bin = ceiling(x_end/bin_dx) !> Clamp to valid bin range - i_start = max(1, min(bins, i_start)) - i_end = max(1, min(bins, i_end)) + start_bin = max(1, min(bins, start_bin)) + end_bin = max(1, min(bins, end_bin)) - if (i_start .eq. i_end) then - flux_tally(i_start) = flux_tally(i_start) + d_step + if (start_bin .eq. end_bin) then + flux_tally(start_bin) = flux_tally(start_bin) + d_step else !> First partial bin - dx = (i_start*bin_dx - x_start)/abs(u) - flux_tally(i_start) = flux_tally(i_start) + dx + dx = (start_bin*bin_dx - x_start)/abs(u) + flux_tally(start_bin) = flux_tally(start_bin) + dx !> Full bins in between - do i = i_start+1, i_end-1 - flux_tally(i) = flux_tally(i) + bin_dx/abs(u) + do bin = start_bin+1, end_bin-1 + flux_tally(bin) = flux_tally(bin) + bin_dx/abs(u) end do !> Last partial bin - dx = (x_end - (i_end-1)*bin_dx)/abs(u) - flux_tally(i_end) = flux_tally(i_end) + dx + dx = (x_end - (end_bin-1)*bin_dx)/abs(u) + flux_tally(end_bin) = flux_tally(end_bin) + dx end if end subroutine tally_segment @@ -1,3 +1,12 @@ -set term x11 persist +#set term x11 persist +set term cairolatex pdf size 5in,3in +set output "flux.tex" -plot 'flux.out' using 1:2 with lines +set key reverse +set grid +set xlabel "Distance [cm]" +set ylabel "$\\phi$" + + +plot 'flux.out' using 1:($2 + $3):($2 - $3) with filledcurves lc rgb "#CCCCFF" title "$\\pm 1 \\sigma$", \ + 'flux.out' using 1:2 with lines lc rgb "blue" title "$\\overbar{\\phi}$" |