diff options
Diffstat (limited to 'main.f90')
| -rw-r--r-- | main.f90 | 33 |
1 files changed, 23 insertions, 10 deletions
@@ -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 |