Minimal working example. No OpenMP yet. Added plotting functionality with Gnuplot

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
diff --git a/Makefile b/Makefile
index 96a64c3..bd3491f 100644
--- a/Makefile
+++ b/Makefile
@@ -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
diff --git a/main.f90 b/main.f90
index 1be36c4..ce1bf6b 100644
--- a/main.f90
+++ b/main.f90
@@ -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