From f045e063cac3c722b63c47c84a1e6962465036c6 Mon Sep 17 00:00:00 2001 From: lewis mcmillan Date: Fri, 15 Nov 2024 14:44:46 +0000 Subject: [PATCH] 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 | 5 +++- src/kernelsMod.f90 | 66 ++++++++++++++++++++++++++++++++++++++------- src/parse/parse.f90 | 5 ++++ src/photon.f90 | 11 +++++--- src/setup.f90 | 9 ++++++- src/sim_state.f90 | 3 +++ src/writer.f90 | 35 +++++++++++++++++++++++- 7 files changed, 118 insertions(+), 16 deletions(-) diff --git a/res/scat_test.toml b/res/scat_test.toml index 254afd90..8ce68fca 100644 --- a/res/scat_test.toml +++ b/res/scat_test.toml @@ -25,4 +25,7 @@ overwrite = true [simulation] iseed = 123456789 -absorb=true \ No newline at end of file +absorb=true +load_checkpoint=false +checkpoint_file="check.ckpt" +checkpoint_every_n=10000 \ No newline at end of file diff --git a/src/kernelsMod.f90 b/src/kernelsMod.f90 index 3ee363a9..21cc52e8 100644 --- a/src/kernelsMod.f90 +++ b/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,13 +205,41 @@ 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!" @@ -218,23 +247,26 @@ subroutine pathlength_scatter(input_file) #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 diff --git a/src/parse/parse.f90 b/src/parse/parse.f90 index 56d8efb8..1c28b00b 100644 --- a/src/parse/parse.f90 +++ b/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 diff --git a/src/photon.f90 b/src/photon.f90 index 1c7d9123..ffbaf501 100644 --- a/src/photon.f90 +++ b/src/photon.f90 @@ -157,7 +157,9 @@ 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 @@ -165,9 +167,12 @@ subroutine slm(this, spectrum, dict, seqs) 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 diff --git a/src/setup.f90 b/src/setup.f90 index 5019dc82..3c8d44db 100644 --- a/src/setup.f90 +++ b/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)) diff --git a/src/sim_state.f90 b/src/sim_state.f90 index 9022b632..89d17e4f 100644 --- a/src/sim_state.f90 +++ b/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 diff --git a/src/writer.f90 b/src/writer.f90 index 5888dcd9..f88a4b6c 100644 --- a/src/writer.f90 +++ b/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