Fixed timing code so it works with openmp. Added checkpoint ability. …

…Currently its limited to saving jmean and number of photons run...

res/scat_test.toml

@@ -25,4 +25,7 @@ overwrite = true
 
 [simulation]
 iseed = 123456789
-absorb=true
+absorb=true
+load_checkpoint=false
+checkpoint_file="check.ckpt"
+checkpoint_every_n=10000

src/kernelsMod.f90

@@ -181,6 +181,7 @@ subroutine pathlength_scatter(input_file)
         use utils,         only : pbar
         use vec4_class,    only : vec4
         use vector_class,  only : vector
+        use writer_mod,    only : checkpoint
 
         !external deps
         use tev_mod, only : tevipc
@@ -204,37 +205,68 @@ subroutine pathlength_scatter(input_file)
         type(tevipc)                  :: tev
         type(seq)                     :: seqs(2)
         type(spectrum_t)              :: spectrum
-
-        call setup(input_file, tev, dects, array, packet, spectrum, dict, distances, image, nscatt, start)
+        real :: tic, toc
+
+        integer :: nphotons_run,pos
+        character(len=128) :: line
+        character(len=:), allocatable :: checkpt_input_file
+
+        if(state%loadckpt)then
+            call setup(input_file, tev, dects, array, packet, spectrum, dict, distances, image, nscatt, start, .false.)
+            open(newunit=j,file=state%ckptfile, access="stream", form="formatted")
+            read(j,"(a)")line
+            pos = scan(line, "=")
+            checkpt_input_file = trim(line(pos+1:))
+
+            read(j,"(a)")line
+            pos = scan(line, "=")
+            read(line(pos+1:),*) nphotons_run
+
+            inquire(j,pos=pos)
+            close(j)
+
+            open(newunit=j,file=state%ckptfile, access="stream", form="unformatted")
+            read(j,pos=pos)jmean
+            close(j)
+
+            call setup(checkpt_input_file, tev, dects, array, packet, spectrum, dict, distances, image, nscatt, start, .true.)
+            state%iseed=state%iseed*101
+            state%nphotons = state%nphotons - nphotons_run
+        else
+            call setup(input_file, tev, dects, array, packet, spectrum, dict, distances, image, nscatt, start, .true.)
+        end if
 
 #ifdef _OPENMP
-        !is state%seed private, i dont think so...
-        !$omp parallel default(none) shared(dict, array, numproc, start, state, bar, jmean, phasor, tev, dects, spectrum)&
-        !$omp& private(ran, id, distances, image, dir, hpoint, history, seqs) reduction(+:nscatt) firstprivate(packet)
+        tic=omp_get_wtime()
+!$omp parallel default(none) shared(dict, array, numproc, start, bar, jmean, input_file, phasor, tev, dects, spectrum)&
+        !$omp& private(ran, id, distances, image, dir, hpoint, history, seqs) reduction(+:nscatt) firstprivate(state, packet)
         numproc = omp_get_num_threads()
         id = omp_get_thread_num()
         if(numproc > state%nphotons .and. id == 0)print*,"Warning, simulation may be underministic due to low photon count!"
         if(state%trackHistory)history = history_stack_t(state%historyFilename, id)
 #elif MPI
     !nothing
 #else
+        call cpu_time(tic)
         numproc = 1
         id = 0
         if(state%trackHistory)history = history_stack_t(state%historyFilename, id)
 #endif
         if(id == 0)print("(a,I3.1,a)"),'Photons now running on', numproc,' cores.'
-
+        state%iseed = state%iseed + id
         ! set seed for rnd generator. id to change seed for each process
-        call init_rng(spread(state%iseed+id, 1, 8), fwd=.true.)
+        call init_rng(spread(state%iseed, 1, 8), fwd=.true.)
         seqs = [seq((id+1)*(state%nphotons/numproc), 2),&
                 seq((id+1)*(state%nphotons/numproc), 3)]
 
         bar = pbar(state%nphotons/ 10)
+
         !$OMP BARRIER
         !$OMP do
-        !loop over photons 
+        !loop over photons
         do j = 1, state%nphotons
             if(mod(j, 10) == 0)call bar%progress()
+            if(mod(j, state%ckptfreq) == 0 .and. id==0)call checkpoint(input_file, state%ckptfile, j, .true.)
 
             ! Release photon from point source
             call packet%emit(spectrum, dict, seqs)
@@ -291,7 +323,12 @@ subroutine pathlength_scatter(input_file)
 #ifdef _OPENMP
 !$OMP end  do
 !$OMP end parallel
+toc=omp_get_wtime()
+#else
+    call cpu_time(toc)
 #endif
+    print*,"Photons/s: ",(state%nphotons / (toc - tic))
+
     call finalise(dict, dects, nscatt, start, history)
     end subroutine pathlength_scatter
 
@@ -413,7 +450,7 @@ end subroutine test_kernel
 
 !####################################################################################################
 !                           Setup and break down helper routines
-    subroutine setup(input_file, tev, dects, array, packet, spectrum, dict, distances, image, nscatt, start)
+    subroutine setup(input_file, tev, dects, array, packet, spectrum, dict, distances, image, nscatt, start, display)
         !! setup simulation by reading in setting file, and setup variables to be used.
 
         !shared data
@@ -451,11 +488,20 @@ subroutine setup(input_file, tev, dects, array, packet, spectrum, dict, distance
         real(kind=wp),    allocatable, intent(out) :: distances(:), image(:,:,:)
         real(kind=wp),                 intent(out) :: nscatt, start
         type(spectrum_t),              intent(out) :: spectrum
+        !> flag to display simulation init settings
+        logical, optional,             intent(in) :: display
 
         ! mpi/mp variables
         integer       :: id
         real(kind=wp) :: chance, threshold
         type(toml_error), allocatable :: error
+        logical :: disp
+
+        if(present(display))then
+            disp = display
+        else
+            disp = .true.
+        end if
 
         chance = 1._wp/10._wp
         threshold = 1e-6_wp
@@ -470,7 +516,7 @@ subroutine setup(input_file, tev, dects, array, packet, spectrum, dict, distance
         end if
         allocate(image(state%grid%nxg,state%grid%nzg,1))
 
-        call display_settings(state, input_file, packet, "Pathlength")
+        if(disp)call display_settings(state, input_file, packet, "Pathlength")
 
         if(state%tev)then
             !init TEV link

src/parse/parse.f90

@@ -202,6 +202,11 @@ subroutine parse_simulation(table, error)
             call get_value(child, "iseed", state%iseed, 123456789)
             call get_value(child, "tev", state%tev, .false.)
             call get_value(child, "absorb", state%absorb, .false.)
+
+            call get_value(child, "load_checkpoint", state%loadckpt, .false.)
+            call get_value(child, "checkpoint_file", state%ckptfile, "check.ckpt")
+            call get_value(child, "checkpoint_every_n", state%ckptfreq, 1000000)
+
         else
             call make_error(error, "Need simulation table in input param file", -1)
             return

src/photon.f90

@@ -157,17 +157,22 @@ type(photon) function init_source(choice)
         end function init_source
 
         subroutine slm(this, spectrum, dict, seqs)
-
+            !! image source
+            ! TODO fix hardcoded size of 200
+            ! investigate +2 error
             use piecewiseMod
             use tomlf,         only : toml_table, get_value
             use random,        only : ran2, seq
             use sim_state_mod, only : state
             use constants,     only : TWOPI
 
             class(photon) :: this
-            type(spectrum_t), intent(in) :: spectrum
+            !> Input image to sample position from
+            type(spectrum_t),           intent(in)    :: spectrum
+            !> Metadata dictionary
             type(toml_table), optional, intent(inout) :: dict
-            type(seq), optional, intent(inout) :: seqs(2)
+            !> random numbers from a sequence. Quasi-Monte Carlo
+            type(seq),        optional, intent(inout) :: seqs(2)
 
             integer :: cell(3)
             real(kind=wp) :: x, y

src/setup.f90

@@ -82,7 +82,12 @@ subroutine directory()
             inquire(directory=trim(fileplace)//"/detectors", exist=detectorsExists)
             inquire(directory=trim(fileplace)//"/phasor", exist=phasorExists)
 #else 
-    error stop "Compiler not supported!"
+    dataExists=.true.
+    jmeanExists=.true.
+    depositExists=.true.
+    detectorsExists=.true.
+    phasorExists=.true.
+    ! error stop "Compiler not supported!"
 #endif
             if(.not. dataExists)then
                 call create_directory("", dataExists, "", .false.)
@@ -151,6 +156,8 @@ subroutine alloc_array(nxg, nyg, nzg)
             !> grid size
             integer, intent(IN) :: nxg, nyg, nzg
 
+            if(allocated(phasor))call dealloc_array()
+
             allocate(phasor(nxg, nyg, nzg), phasorGLOBAL(nxg, nyg, nzg))
             allocate(jmean(nxg, nyg, nzg), jmeanGLOBAL(nxg, nyg, nzg))
             allocate(absorb(nxg, nyg, nzg), absorbGLOBAL(nxg, nyg, nzg))

src/sim_state.f90

@@ -36,6 +36,9 @@ module sim_state_mod
         logical :: trackHistory
         !> Boolean to indicate whether to store absoption data.
         logical :: absorb
+        integer :: ckptfreq
+        logical :: loadckpt
+        character(len=:), allocatable :: ckptfile
     end type settings_t
 
     !> global var that stores simulation state

src/writer.f90

@@ -19,7 +19,7 @@ module writer_mod
     end interface raw_write
 
     private
-    public :: normalise_fluence, write_data, write_detected_photons
+    public :: normalise_fluence, write_data, write_detected_photons, checkpoint
 
     contains
         subroutine normalise_fluence(grid, array, nphotons)
@@ -332,4 +332,37 @@ subroutine write_3d_r4_nrrd(array, filename, overwrite, dict)
             close(u)
 
         end subroutine write_3d_r4_nrrd
+
+        subroutine checkpoint(toml_filename, filename, nphotons_run, overwrite)
+
+            use iarray, only : jmean
+
+            !> filename of toml file used in simulation
+            character(*), intent(IN)    :: toml_filename
+            !> name of checkpoint file to be saved
+            character(*), intent(IN)    :: filename
+            !> flag which determines if file is to be overwritten or adjusted
+            logical,      intent(IN) :: overwrite
+            !> number of photons run up to checkpoint
+            integer,      intent(IN) :: nphotons_run
+
+            character(len=:), allocatable :: file
+            integer :: u
+
+            if(check_file(filename) .and. .not. overwrite)then
+                file = get_new_file_name(filename)
+            else
+                file = filename
+            end if
+
+            open(newunit=u,file=file)
+            write(u,"(a,a)")"tomlfile=",toml_filename
+            write(u,"(a,i0)")"photons_run=",nphotons_run
+            close(u)
+
+            open(newunit=u,file=file,access="stream",form="unformatted",position="append")
+            write(u)jmean
+            close(u)
+
+        end subroutine checkpoint
 end module writer_mod