Initial OpenMP support

2 files changed, 25 insertions, 12 deletions

diff --git a/Makefile b/Makefile
index bd3491f..83f9d39 100644
--- a/Makefile
+++ b/Makefile
@@ -1,5 +1,5 @@
 FC = gfortran
-FFLAGS = -march=native -O3
+FFLAGS = -march=native -O3 -fopenmp
 XNAME = monte_carlo.x
 NPROC := $(shell nproc)
 
@@ -13,7 +13,7 @@ build: main.f90
 	$(FC) $(FFLAGS) main.f90 -o $(XNAME)
 
 run: build
-	OMP_NUM_THREADS=$(NPROC) ./$(XNAME)
+	ulimit -s unlimited && OMP_NUM_THREADS=$(NPROC) && ./$(XNAME)
 
 plot: run
 	gnuplot plot.gnu
diff --git a/main.f90 b/main.f90
index ce1bf6b..18eb47a 100644
--- a/main.f90
+++ b/main.f90
@@ -1,10 +1,11 @@
 program monte_carlo
 
+    use omp_lib
     use, intrinsic :: iso_fortran_env
     implicit none
 
     !> Initialize some variables for neutrons
-    integer(int64), parameter :: n = 1e7
+    integer(int64), parameter :: n = 1e8
     integer(int64), parameter :: bins = 800
     real(real64), parameter :: eps = 1.0e-12_real64
 
@@ -27,12 +28,16 @@ program monte_carlo
 
     !> Initialize RNG to be actually random
     call random_seed()
-    call random_init(.false., .false.)
+    call random_init(.false., .true.)
 
     !> Provide all neutrons an initial location/angle
     call initialize_neutrons(n, px, mu)
 
+    !$OMP PARALLEL PRIVATE(thread_tally, sigma_t, sigma_s, q, d_flight, d_boundary, d_step, xi, region, i)
+    thread_tally=0.0_real64
+
     !> The main monte carlo loop 
+    !$OMP DO
     main: do i = 1,n
 
         !> Skip if the neutron is already lost
@@ -112,14 +117,20 @@ program monte_carlo
 
         end do single
 
+
     end do main
+    !$OMP END DO
 
+    !$OMP CRITICAL
     flux_tally = flux_tally + thread_tally
-    flux_tally = flux_tally / thread_tally(100)
+    !$OMP END CRITICAL
+    !$OMP END PARALLEL
+    flux_tally = flux_tally / flux_tally(bins/8)
+    !flux_tally = flux_tally / (real(n, real64) * dx)
 
     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)
+        write(10,'(2E15.7)') (i-0.5_real64)*dx, flux_tally(i)
     end do write
 
 
@@ -156,12 +167,12 @@ contains
         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
+        where(r_region .lt. 100.0_real64)
+            px = eps + r_position * (2.0_real64 - 2.0_real64*eps)
+        else where(r_region .lt. 105.0_real64)
+            px = 2.0_real64 + eps + r_position * (1.0_real64 - 2.0_real64*eps)
+        else where(r_region .lt. 106.0_real64)
+            px = 5.0_real64 + eps + r_position * (1.0_real64 - 2.0_real64*eps)
         end where
 
         !> Assign random angles to each particle
@@ -176,8 +187,10 @@ contains
         real(real64) :: x_start, x_end, x_left, x_right
         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