From 7b36bcc0d7a5e9f7f928508b177e3e331fcb9d0f Mon Sep 17 00:00:00 2001 From: "V. Raffuzzi" Date: Fri, 13 Dec 2024 16:49:43 +0000 Subject: [PATCH] Addressing reviews and fixing github tests --- .github/workflows/build-and-test.yml | 29 ++- .gitignore | 3 +- CMakeLists.txt | 23 +- DataStructures/heapQueue_class.f90 | 9 +- .../Tests/baseMgNeutronDatabase_iTest.f90 | 9 +- .../mgNeutronData/mgNeutronMaterial_inter.f90 | 6 +- .../Tests/particleDungeon_test.f90 | 9 + ParticleObjects/particleDungeon_class.f90 | 211 +++--------------- ParticleObjects/particle_class.f90 | 52 ++++- PhysicsPackages/eigenPhysicsPackage_class.f90 | 12 +- .../fixedSourcePhysicsPackage_class.f90 | 4 +- SharedModules/errors_mod.f90 | 4 - SharedModules/mpi_func.f90 | 58 ++--- Tallies/scoreMemory_class.f90 | 8 +- docs/Running.rst | 4 +- 15 files changed, 189 insertions(+), 252 deletions(-) diff --git a/.github/workflows/build-and-test.yml b/.github/workflows/build-and-test.yml index 7bb7a6a9e..124f174ef 100644 --- a/.github/workflows/build-and-test.yml +++ b/.github/workflows/build-and-test.yml @@ -10,10 +10,11 @@ jobs: build-and-test: strategy: matrix: - compiler: [gfortran8, gfortran9, gfortran10] + compiler: [gfortran10, gfortran11, gfortran12] runs-on: ubuntu-20.04 container: image: mikolajkowalski/scone-test:${{matrix.compiler}}_pfu4 + options: --env OMPI_ALLOW_RUN_AS_ROOT=1 --env OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1 steps: - uses: actions/checkout@v3 - name: CompileAndTest @@ -22,12 +23,13 @@ jobs: cd build cmake .. make -j - make test + ctest --output-on-faliure cd - build-and-test-debug: runs-on: ubuntu-20.04 container: - image: mikolajkowalski/scone-test:gfortran10_pfu4 + image: mikolajkowalski/scone-test:gfortran12_pfu4 + options: --env OMPI_ALLOW_RUN_AS_ROOT=1 --env OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1 steps: - uses: actions/checkout@v3 - name: CompileAndTest @@ -36,12 +38,13 @@ jobs: cd build cmake -DDEBUG=ON .. make -j - make test + ctest --output-on-faliure cd - build-and-test-no-openmp: runs-on: ubuntu-20.04 container: - image: mikolajkowalski/scone-test:gfortran10_pfu4 + image: mikolajkowalski/scone-test:gfortran12_pfu4 + options: --env OMPI_ALLOW_RUN_AS_ROOT=1 --env OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1 steps: - uses: actions/checkout@v3 - name: CompileAndTest @@ -50,5 +53,19 @@ jobs: cd build cmake -DOPENMP=OFF .. make -j - make test + ctest --output-on-faliure + cd - + build-and-test-no-mpi: + runs-on: ubuntu-20.04 + container: + image: mikolajkowalski/scone-test:gfortran12_pfu4 + steps: + - uses: actions/checkout@v3 + - name: CompileAndTest + run : | + mkdir build + cd build + cmake -DMPI=OFF .. + make -j + ctest --output-on-faliure cd - diff --git a/.gitignore b/.gitignore index d51b089f5..32fd9eea9 100644 --- a/.gitignore +++ b/.gitignore @@ -7,6 +7,7 @@ cream.egg-info/ Build build -# Ignore all hidden files (except gitignore) +# Ignore all hidden files (except gitignore and the github folder) .* !/.gitignore +!/.github/* diff --git a/CMakeLists.txt b/CMakeLists.txt index 13c1f3d28..7cc3fdabd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -92,7 +92,7 @@ endif() # Add environmental variable to default search directories list(APPEND CMAKE_PREFIX_PATH $ENV{LAPACK_INSTALL}) -find_package(LAPACK REQUIRED ) +find_package(LAPACK REQUIRED) message(STATUS ${LAPACK_LIBRARIES}) # Dependencies for BUILD_TESTS @@ -144,9 +144,10 @@ get_property(SRCS GLOBAL PROPERTY SRCS_LIST) # Compile library add_library(scone STATIC ${SRCS}) -target_compile_options(scone PRIVATE ${scone_extra_flags} ) -target_link_libraries(scone PUBLIC ${LAPACK_LIBRARIES} ) +target_compile_options(scone PRIVATE ${scone_extra_flags}) +target_link_libraries(scone PUBLIC ${LAPACK_LIBRARIES}) if(MPI) + add_compile_definitions(MPI) target_link_libraries(scone PUBLIC MPI::MPI_Fortran) endif() @@ -156,8 +157,8 @@ endif() #################################################################################################### # COMPILE SOLVERS -add_executable(scone.out ./Apps/scone.f90 ) -target_link_libraries(scone.out scone ) +add_executable(scone.out ./Apps/scone.f90) +target_link_libraries(scone.out scone) #################################################################################################### # COMPILE UNIT TESTS @@ -172,16 +173,22 @@ if(BUILD_TESTS) list(APPEND UNIT_TESTS_RELATIVE ${_testPath}) endforeach() + if(MPI) + set(MAX_PES_OPTION MAX_PES 1) + else() + set(MAX_PES_OPTION "") + endif() + add_pfunit_ctest(unitTests TEST_SOURCES ${UNIT_TESTS_RELATIVE} LINK_LIBRARIES scone - MAX_PES 1 + ${MAX_PES_OPTION} ) # pFUnit may have a bug which causes a unused variable `class(Test), allocatable :: t` be # present if the suite contains only a TestCase and its methods # We need to suppress this warning for clarity - target_compile_options(unitTests PRIVATE "-Wno-unused-variable" ) + target_compile_options(unitTests PRIVATE "-Wno-unused-variable") #################################################################################################### # COMPILE INTEGRATION TESTS @@ -203,7 +210,7 @@ if(BUILD_TESTS) # pFUnit may have a bug which causes a unused variable `class(Test), allocatable :: t` be # present if the suite contains only a TestCase and its methods # We need to suppress this warning for clarity - target_compile_options(integrationTests PRIVATE "-Wno-unused-variable" ) + target_compile_options(integrationTests PRIVATE "-Wno-unused-variable") endif() diff --git a/DataStructures/heapQueue_class.f90 b/DataStructures/heapQueue_class.f90 index 0a101cf93..6daf01854 100644 --- a/DataStructures/heapQueue_class.f90 +++ b/DataStructures/heapQueue_class.f90 @@ -96,14 +96,21 @@ subroutine push(self, val) self % size = self % size + 1 self % heap(self % size) = val + ! If the heap is of size 1 there is no need to order it + ! Also, avoid test fail in debug mode, since parent would be 0 and the + ! code is looking for self % heap(parent) inside the while condition + if (self % size == 1) return + ! Shift the new value up the heap to restore the heap property child = self % size parent = child / 2 - do while (child > 1 .and. self % heap(parent) < self % heap(child)) + do while (self % heap(parent) < self % heap(child)) call swap(self % heap(parent), self % heap(child)) child = parent parent = child / 2 + ! As above: avoid error in debug mode, caused by trying to access self % heap(0) + if (parent == 0) return end do end subroutine push diff --git a/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 b/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 index 5ce32ddca..daa6784d9 100644 --- a/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 +++ b/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 @@ -7,6 +7,7 @@ module baseMgNeutronDatabase_iTest use dictionary_class, only : dictionary use dictParser_func, only : charToDict use particle_class, only : particle + use RNG_class, only : RNG ! Nuclear Data Objects & Interfaces use baseMgNeutronDatabase_class, only : baseMgNeutronDatabase, baseMgNeutronDatabase_CptrCast, & @@ -59,12 +60,12 @@ subroutine testBaseMgNeutronDatabaseWithP0() type(dictionary) :: matMenuDict type(particle) :: p type(neutronMacroXSs) :: xss + type(RNG), target :: pRNG type(baseMgNeutronMaterial),pointer :: mat class(baseMgNeutronMaterial),pointer :: matClass class(reactionHandle), pointer :: reac real(defReal),parameter :: TOL = 1.0E-6_defReal - data_ptr => database ! Load materialMenu @@ -81,6 +82,8 @@ subroutine testBaseMgNeutronDatabaseWithP0() @assertEqual(4, database % nGroups()) ! Test getting Transport XS + ! Associate pointer to pass tests in debug mode + p % pRNG => pRNG p % G = 1 @assertEqual(2.1_defReal, database % getTrackingXS(p, 1, MATERIAL_XS), TOL) @@ -185,12 +188,12 @@ subroutine testBaseMgNeutronDatabaseWithP1() type(dictionary) :: matMenuDict type(particle) :: p type(neutronMacroXSs) :: xss + type(RNG), target :: pRNG type(baseMgNeutronMaterial),pointer :: mat class(baseMgNeutronMaterial),pointer :: matClass class(reactionHandle), pointer :: reac real(defReal),parameter :: TOL = 1.0E-6_defReal - data_ptr => database ! Load materialMenu @@ -207,6 +210,8 @@ subroutine testBaseMgNeutronDatabaseWithP1() @assertEqual(4, database % nGroups()) ! Test getting Transport XS + ! Associate pointer to pass tests in debug mode + p % pRNG => pRNG p % G = 1 @assertEqual(2.1_defReal, database % getTrackingXS(p, 1, MATERIAL_XS), TOL) diff --git a/NuclearData/mgNeutronData/mgNeutronMaterial_inter.f90 b/NuclearData/mgNeutronData/mgNeutronMaterial_inter.f90 index 4b8aa0cfa..daa90c37f 100644 --- a/NuclearData/mgNeutronData/mgNeutronMaterial_inter.f90 +++ b/NuclearData/mgNeutronData/mgNeutronMaterial_inter.f90 @@ -111,11 +111,15 @@ subroutine getMacroXSs_byP(self, xss, p) class(mgNeutronMaterial), intent(in) :: self type(neutronMacroXSs), intent(out) :: xss class(particle), intent(in) :: p + integer(shortInt) :: matIdx character(100), parameter :: Here = 'getMacroXSs_byP (mgNeutronMateerial_inter.f90)' if (.not. p % isMG) call fatalError(Here, 'CE particle was given to MG data') - associate (matCache => cache_materialCache(p % matIdx())) + ! Store p % matIdx() in a dedicated variable to avoid compilation errors with gfortran >= 13.2 + matIdx = p % matIdx() + + associate (matCache => cache_materialCache(matIdx)) if (matCache % G_tail /= p % G) then ! Get cross sections diff --git a/ParticleObjects/Tests/particleDungeon_test.f90 b/ParticleObjects/Tests/particleDungeon_test.f90 index bea00ad86..7628e7746 100644 --- a/ParticleObjects/Tests/particleDungeon_test.f90 +++ b/ParticleObjects/Tests/particleDungeon_test.f90 @@ -4,7 +4,9 @@ module particleDungeon_test use RNG_class, only : RNG use particle_class, only : particle, particleState use particleDungeon_class, only : particleDungeon +#ifdef MPI use mpi_func, only : mpiInitTypes, MPI_COMM_WORLD +#endif use funit implicit none @@ -213,13 +215,18 @@ subroutine testNormPopDown() real(defReal), parameter :: TOL = 1.0E-9 character(100),parameter :: Here = 'testNormPopDown (particleDungeon_test.f90)' +#ifdef MPI call mpi_comm_size(MPI_COMM_WORLD, worldSize, ierr) if (worldSize > 1) & call fatalError(Here, 'This test cannot be run with multiple MPI processes') ! Initialise MPI types needed for this procedure + ! NOTE: This is necessary because the normalisation uses some mpi procedure + ! with data types manually defined inside mpiInitTypes. During the tests, + ! mpiInit and mpiInitTypes aren't called, so this is done manually here call mpiInitTypes() +#endif ! Initialise call dungeon % init(10) @@ -260,10 +267,12 @@ subroutine testNormPopUp() real(defReal), parameter :: TOL = 1.0E-9 character(100),parameter :: Here = 'testNormPopUp (particleDungeon_test.f90)' +#ifdef MPI call mpi_comm_size(MPI_COMM_WORLD, worldSize, ierr) if (worldSize > 1) & call fatalError(Here, 'This test cannot be run with multiple MPI processes') +#endif ! Initialise call dungeon % init(20) diff --git a/ParticleObjects/particleDungeon_class.f90 b/ParticleObjects/particleDungeon_class.f90 index 364ba1acd..b92d9fffe 100644 --- a/ParticleObjects/particleDungeon_class.f90 +++ b/ParticleObjects/particleDungeon_class.f90 @@ -3,16 +3,17 @@ module particleDungeon_class use numPrecision use errors_mod, only : fatalError use genericProcedures, only : numToChar, swap - use particle_class, only : particle, particleState + use particle_class, only : particle, particleStateData, particleState use RNG_class, only : RNG use heapQueue_class, only : heapQueue use mpi_func, only : isMPIMaster, getMPIWorldSize, getMPIRank, getOffset #ifdef MPI use mpi_func, only : mpi_gather, mpi_allgather, mpi_send, mpi_recv, & - mpi_Bcast, MPI_COMM_WORLD, MASTER_RANK, MPI_DOUBLE, & - MPI_INT, MPI_LONG_LONG, MPI_PARTICLE_STATE, & - MPI_STATUS_IGNORE, particleStateDummy + mpi_bcast, MPI_COMM_WORLD, MPI_STATUS_IGNORE, & + MASTER_RANK, MPI_PARTICLE_STATE, MPI_DEFREAL, & + MPI_SHORTINT, MPI_LONGINT + #endif implicit none @@ -102,7 +103,6 @@ module particleDungeon_class procedure :: setSize procedure :: printToFile procedure :: sortByBroodID - procedure :: samplingWoReplacement ! Private procedures procedure, private :: detain_particle @@ -429,8 +429,8 @@ subroutine normSize(self, totPop, rand) real(defReal) :: threshold, rn integer(longInt) :: seedTemp integer(shortInt) :: maxbroodID, totSites, excess, heapSize, & - n_duplicates, i, j, n_copies, count, nRanks, & - rank + n_duplicates, n_copies, count, nRanks, & + rank, i, j integer(longInt), dimension(:), allocatable :: seeds integer(shortInt), dimension(:), allocatable :: keepers, popSizes #ifdef MPI @@ -453,7 +453,7 @@ subroutine normSize(self, totPop, rand) #ifdef MPI ! Get the population sizes of all ranks into the array popSizes in master branch - call mpi_gather(self % pop, 1, MPI_INT, popSizes, 1, MPI_INT, MASTER_RANK, MPI_COMM_WORLD, error) + call mpi_gather(self % pop, 1, MPI_SHORTINT, popSizes, 1, MPI_SHORTINT, MASTER_RANK, MPI_COMM_WORLD, error) #endif ! In the master process, calculate sampling threshold for the whole population @@ -505,9 +505,9 @@ subroutine normSize(self, totPop, rand) ! Broadcast threshold, excess and random number seeds to all processes #ifdef MPI - call MPI_Bcast(threshold, 1, MPI_DOUBLE, MASTER_RANK, MPI_COMM_WORLD) - call MPI_Bcast(excess, 1, MPI_INT, MASTER_RANK, MPI_COMM_WORLD) - call MPI_Bcast(seeds, nRanks, MPI_LONG_LONG, MASTER_RANK, MPI_COMM_WORLD) + call mpi_bcast(threshold, 1, MPI_DEFREAL, MASTER_RANK, MPI_COMM_WORLD) + call mpi_bcast(excess, 1, MPI_SHORTINT, MASTER_RANK, MPI_COMM_WORLD) + call mpi_bcast(seeds, nRanks, MPI_LONGINT, MASTER_RANK, MPI_COMM_WORLD) #endif ! Get local process rank and initialise local rng with the correct seed @@ -582,7 +582,7 @@ subroutine normSize(self, totPop, rand) #ifdef MPI ! Get the updated population numbers from all processes - call mpi_allgather(self % pop, 1, MPI_INT, popSizes, 1, MPI_INT, MPI_COMM_WORLD, error) + call mpi_allgather(self % pop, 1, MPI_SHORTINT, popSizes, 1, MPI_SHORTINT, MPI_COMM_WORLD, error) #endif ! Check that normalisation worked @@ -607,16 +607,16 @@ subroutine loadBalancing(self, totPop, nRanks, rank, popSizes) integer(shortInt), dimension(:), intent(in) :: popSizes integer(shortInt), dimension(:), allocatable :: rankOffsets integer(shortInt), dimension(2) :: offset, targetOffset - integer(shortInt) :: mpiOffset, excess, error - type(particleState), dimension(:), allocatable :: stateBuffer - type(particleStateDummy), dimension(:), allocatable :: dummyBuffer + integer(shortInt) :: mpiOffset, excess, error, i + class(particleState), dimension(:), allocatable :: stateBuffer + class(particleStateData), dimension(:), allocatable :: dataBuffer ! Get expected particle population in each process via the offset mpiOffset = getOffset(totPop) ! Communicates the offset from all processes to all processes allocate(rankOffsets(nRanks)) - call mpi_allgather(mpiOffset, 1, MPI_INT, rankOffsets, 1, MPI_INT, MPI_COMM_WORLD, error) + call mpi_allgather(mpiOffset, 1, MPI_SHORTINT, rankOffsets, 1, MPI_SHORTINT, MPI_COMM_WORLD, error) ! Calculate actual and target cumulative number of sites in the processes before offset(1) = sum(popSizes(1 : rank)) @@ -634,25 +634,26 @@ subroutine loadBalancing(self, totPop, nRanks, rank, popSizes) if (excess > 0) then ! Send particles from the end of the dungeon to the rank above - stateBuffer = self % prisoners(self % pop - excess + 1 : self % pop) - call initStateDummies(stateBuffer, dummyBuffer) - call mpi_send(dummyBuffer, excess, MPI_PARTICLE_STATE, rank + 1, rank, MPI_COMM_WORLD, error) + dataBuffer = self % prisoners(self % pop - excess + 1 : self % pop) + call mpi_send(dataBuffer, excess, MPI_PARTICLE_STATE, rank + 1, rank, MPI_COMM_WORLD, error) self % pop = self % pop - excess elseif (excess < 0) then ! Receive particles from the rank above and store them at the end of the dungeon excess = -excess - allocate(dummyBuffer(excess)) - call mpi_recv(dummyBuffer, excess, MPI_PARTICLE_STATE, rank + 1, rank + 1, & + allocate(dataBuffer(excess), stateBuffer(excess)) + call mpi_recv(dataBuffer, excess, MPI_PARTICLE_STATE, rank + 1, rank + 1, & MPI_COMM_WORLD, MPI_STATUS_IGNORE, error) - call createStatesFromDummies(dummyBuffer, stateBuffer) + do i = 1, abs(excess) + stateBuffer(i) = dataBuffer(i) + end do self % prisoners(self % pop + 1 : self % pop + excess) = stateBuffer self % pop = self % pop + excess end if - if (excess /= 0) deallocate(stateBuffer, dummyBuffer) + if (allocated(dataBuffer)) deallocate(dataBuffer) ! If needed, send/receive particle states from/to the beginning of the dungeon excess = offset(1) - targetOffset(1) @@ -661,9 +662,8 @@ subroutine loadBalancing(self, totPop, nRanks, rank, popSizes) ! Send particles from the beginning of the dungeon to the rank below excess = -excess - stateBuffer = self % prisoners(1 : excess) - call initStateDummies(stateBuffer, dummyBuffer) - call mpi_send(dummyBuffer, excess, MPI_PARTICLE_STATE, rank - 1, rank, MPI_COMM_WORLD, error) + dataBuffer = self % prisoners(1 : excess) + call mpi_send(dataBuffer, excess, MPI_PARTICLE_STATE, rank - 1, rank, MPI_COMM_WORLD, error) ! Move the remaining particles to the beginning of the dungeon self % prisoners(1 : self % pop - excess) = self % prisoners(excess + 1 : self % pop) @@ -672,10 +672,12 @@ subroutine loadBalancing(self, totPop, nRanks, rank, popSizes) elseif (excess > 0) then ! Receive particles from the rank below and store them at the beginning of the dungeon - allocate(dummyBuffer(excess)) - call mpi_recv(dummyBuffer, excess, MPI_PARTICLE_STATE, rank - 1, rank - 1, & + allocate(dataBuffer(excess), stateBuffer(excess)) + call mpi_recv(dataBuffer, excess, MPI_PARTICLE_STATE, rank - 1, rank - 1, & MPI_COMM_WORLD, MPI_STATUS_IGNORE, error) - call createStatesFromDummies(dummyBuffer, stateBuffer) + do i = 1, abs(excess) + stateBuffer(i) = dataBuffer(i) + end do self % prisoners(excess + 1 : self % pop + excess) = self % prisoners(1 : self % pop) self % prisoners(1 : excess) = stateBuffer self % pop = self % pop + excess @@ -685,157 +687,6 @@ subroutine loadBalancing(self, totPop, nRanks, rank, popSizes) end subroutine loadBalancing #endif - !! - !! Helper procedure for MPI loadBalancing - !! - !! Given an input array of particleState, it allocates a vector of particleStateDummy - !! of the same length and copies all the particleState variables. It returns - !! the array of particleStateDummy - !! -#ifdef MPI - subroutine initStateDummies(states, dummies) - type(particleState), dimension(:), intent(in) :: states - type(particleStateDummy), dimension(:), allocatable, intent(out) :: dummies - integer(shortInt) :: i, N - - ! Allocate particleStateDummy array - N = size(states) - allocate(dummies(N)) - - ! Copy particleState attributes - do i = 1, N - dummies(i) % r = states(i) % r - dummies(i) % dir = states(i) % dir - end do - - dummies % wgt = states % wgt - dummies % E = states % E - dummies % G = states % G - dummies % isMG = states % isMG - dummies % type = states % type - dummies % time = states % time - dummies % matIdx = states % matIdx - dummies % uniqueID = states % uniqueID - dummies % cellIdx = states % cellIdx - dummies % collisionN = states % collisionN - dummies % broodID = states % broodID - - end subroutine initStateDummies -#endif - - !! - !! Helper procedure for MPI loadBalancing - !! - !! Given an input array of particleStateDummy, it allocates a vector of particleState - !! of the same length and copies all the particleStateDummy variables. It returns - !! the array of particleState - !! -#ifdef MPI - subroutine createStatesFromDummies(dummies, states) - type(particleStateDummy), dimension(:), intent(in) :: dummies - type(particleState), dimension(:), allocatable, intent(out) :: states - integer(shortInt) :: i, N - - ! Allocate particleState array - N = size(dummies) - allocate(states(N)) - - ! Copy particleStateDummy attributes - do i = 1, N - states(i) % r = dummies(i) % r - states(i) % dir = dummies(i) % dir - end do - - states % wgt = dummies % wgt - states % E = dummies % E - states % G = dummies % G - states % isMG = dummies % isMG - states % type = dummies % type - states % time = dummies % time - states % matIdx = dummies % matIdx - states % uniqueID = dummies % uniqueID - states % cellIdx = dummies % cellIdx - states % collisionN = dummies % collisionN - states % broodID = dummies % broodID - - end subroutine createStatesFromDummies -#endif - - - !! - !! Normalise total number of particles in the dungeon to match the provided number. - !! Randomly duplicate or remove particles to match the number, performing - !! sampling without replacement (Knuth S algorithm). - !! - !! Does not take weight of a particle into account! - !! - !! NOTE: this procedure is not currently used - !! - subroutine samplingWoReplacement(self, N, rand) - class(particleDungeon), intent(inout) :: self - integer(shortInt), intent(in) :: N - class(RNG), intent(inout) :: rand - integer(shortInt) :: excessP, n_copies, n_duplicates - integer(shortInt) :: i, idx, maxBroodID - integer(shortInt), dimension(:), allocatable :: duplicates - character(100), parameter :: Here = 'samplingWoReplacement (particleDungeon_class.f90)' - - ! Protect against invalid N - if (N > size(self % prisoners)) then - call fatalError(Here,'Requested size: '//numToChar(N) //& - 'is greater then max size: '//numToChar(size(self % prisoners))) - else if (N <= 0) then - call fatalError(Here,'Requested size: '//numToChar(N) //' is not +ve') - end if - - ! Determine the maximum brood ID and sort the dungeon - maxBroodID = maxval(self % prisoners(1:self % pop) % broodID) - call self % sortByBroodID(maxbroodID) - - ! Calculate excess particles to be removed - excessP = self % pop - N - - if (excessP > 0) then ! Reduce population with reservoir sampling - do i = N + 1, self % pop - ! Select new index. Copy data if it is in the safe zone (<= N). - idx = int(i * rand % get()) + 1 - if (idx <= N) then - self % prisoners(idx) = self % prisoners(i) - end if - end do - self % pop = N - - else if (excessP < 0) then ! Clone randomly selected particles - ! For a massive undersampling duplicate (or n-plicate) particles - excessP = -excessP - n_copies = excessP / self % pop - n_duplicates = modulo(excessP, self % pop) - - ! Copy all particles the maximum possible number of times - do i = 1, n_copies - self % prisoners(self % pop * i + 1 : self % pop * (i + 1)) = self % prisoners(1:self % pop) - end do - - ! Choose the remainder particles to duplicate without replacement - duplicates = [(i, i = 1, n_duplicates)] - do i = n_duplicates + 1, self % pop - idx = int(i * rand % get()) + 1 - if (idx <= n_duplicates) then - duplicates(idx) = i - end if - end do - self % pop = self % pop * (n_copies + 1) - - ! Copy the duplicated particles at the end - do i = 1, n_duplicates - self % prisoners(self % pop + i) = self % prisoners(duplicates(i)) - end do - self % pop = N - - end if - - end subroutine samplingWoReplacement - !! !! Reorder the dungeon so the brood ID is in the ascending order !! diff --git a/ParticleObjects/particle_class.f90 b/ParticleObjects/particle_class.f90 index f5b21c484..1338221dc 100644 --- a/ParticleObjects/particle_class.f90 +++ b/ParticleObjects/particle_class.f90 @@ -25,6 +25,9 @@ module particle_class !! !! Particle compressed for storage !! + !! particleStateData contains only the properties (useful for MPI); it is extended into + !! particleState which includes the procedures too and is used for tallying. + !! !! Public Members: !! wgt -> Weight of the particle !! r -> Global Position of the particle [cm] @@ -42,12 +45,7 @@ module particle_class !! in the particleDungeon so they can be sorted, which is necessary for reproducibility !! with OpenMP !! - !! Interface: - !! assignemnt(=) -> Build particleState from particle - !! operator(.eq.) -> Return True if particle are exactly the same - !! display -> Print debug information about the state to the console - !! - type, public :: particleState + type, public :: particleStateData real(defReal) :: wgt = ZERO ! Particle weight real(defReal),dimension(3) :: r = ZERO ! Global position real(defReal),dimension(3) :: dir = ZERO ! Global direction @@ -61,15 +59,28 @@ module particle_class integer(shortInt) :: uniqueID = -1 ! Unique id at the lowest coord level integer(shortInt) :: collisionN = 0 ! Number of collisions integer(shortInt) :: broodID = 0 ! ID of the source particle + end type particleStateData + + !! + !! Extension of particleStateData, which includes procedure + !! + !! Interface: + !! assignemnt(=) -> Build particleState from particle + !! operator(.eq.) -> Return True if particle are exactly the same + !! display -> Print debug information about the state to the console + !! + type, public, extends(particleStateData) :: particleState contains generic :: assignment(=) => fromParticle + generic :: assignment(=) => fromParticleStateData generic :: operator(.eq.) => equal_particleState procedure :: display => display_particleState - procedure :: fromParticle => particleState_fromParticle procedure :: kill => kill_particleState + procedure :: fromParticle => particleState_fromParticle + procedure :: fromParticleStateData => particleState_fromParticleStateData ! Private procedures - procedure,private :: equal_particleState + procedure, private :: equal_particleState end type particleState @@ -663,6 +674,31 @@ subroutine particleState_fromParticle(LHS,RHS) end subroutine particleState_fromParticle + !! + !! Copy particleStateData into phase coordinates + !! + subroutine particleState_fromParticleStateData(LHS,RHS) + class(particleState), intent(out) :: LHS + class(particleStateData), intent(in) :: RHS + + LHS % wgt = RHS % wgt + LHS % r = RHS % r + LHS % dir = RHS % dir + LHS % E = RHS % E + LHS % G = RHS % G + LHS % isMG = RHS % isMG + LHS % type = RHS % type + LHS % time = RHS % time + + ! Save all indexes + LHS % matIdx = RHS % matIdx + LHS % uniqueID = RHS % uniqueId + LHS % cellIdx = RHS % cellIdx + LHS % collisionN = RHS % collisionN + LHS % broodID = RHS % broodID + + end subroutine particleState_fromParticleStateData + !! !! Define equal operation on phase coordinates !! Phase coords are equal if all their components are the same diff --git a/PhysicsPackages/eigenPhysicsPackage_class.f90 b/PhysicsPackages/eigenPhysicsPackage_class.f90 index aef203e8a..670e10292 100644 --- a/PhysicsPackages/eigenPhysicsPackage_class.f90 +++ b/PhysicsPackages/eigenPhysicsPackage_class.f90 @@ -8,8 +8,8 @@ module eigenPhysicsPackage_class printSectionEnd, printSeparatorLine use mpi_func, only : isMPIMaster, getWorkshare, getOffset, getMPIRank #ifdef MPI - use mpi_func, only : MASTER_RANK, MPI_Bcast, MPI_INT, MPI_COMM_WORLD, & - MPI_DOUBLE, mpi_reduce, MPI_SUM + use mpi_func, only : mpi_reduce, mpi_bcast, MPI_COMM_WORLD, & + MASTER_RANK, MPI_SUM, MPI_DEFREAL, MPI_SHORTINT #endif use hashFunctions_func, only : FNV_1 use dictionary_class, only : dictionary @@ -290,7 +290,7 @@ subroutine cycles(self, tally, tallyAtch, N_cycles) #ifdef MPI ! Broadcast k_eff obtained in the master to all processes - call MPI_Bcast(k_new, 1, MPI_DOUBLE, MASTER_RANK, MPI_COMM_WORLD) + call mpi_bcast(k_new, 1, MPI_DEFREAL, MASTER_RANK, MPI_COMM_WORLD) #endif ! Load new k-eff estimate into next cycle dungeon @@ -310,9 +310,9 @@ subroutine cycles(self, tally, tallyAtch, N_cycles) #ifdef MPI ! Print the population numbers referred to all processes to screen - call mpi_reduce(nStart, nTemp, 1, MPI_INT, MPI_SUM, MASTER_RANK, MPI_COMM_WORLD, error) + call mpi_reduce(nStart, nTemp, 1, MPI_SHORTINT, MPI_SUM, MASTER_RANK, MPI_COMM_WORLD, error) nStart = nTemp - call mpi_reduce(nEnd, nTemp, 1, MPI_INT, MPI_SUM, MASTER_RANK, MPI_COMM_WORLD, error) + call mpi_reduce(nEnd, nTemp, 1, MPI_SHORTINT, MPI_SUM, MASTER_RANK, MPI_COMM_WORLD, error) nEnd = nTemp #endif @@ -474,7 +474,7 @@ subroutine init(self, dict) ! Broadcast seed to all processes #ifdef MPI - call MPI_Bcast(seed_temp, 1, MPI_INT, MASTER_RANK, MPI_COMM_WORLD) + call mpi_bcast(seed_temp, 1, MPI_SHORTINT, MASTER_RANK, MPI_COMM_WORLD) #endif seed = seed_temp diff --git a/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 b/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 index b3905d6d3..6698ba830 100644 --- a/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 +++ b/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 @@ -8,7 +8,7 @@ module fixedSourcePhysicsPackage_class printSeparatorLine use mpi_func, only : isMPIMaster, getWorkshare, getOffset, getMPIRank #ifdef MPI - use mpi_func, only : MASTER_RANK, MPI_Bcast, MPI_INT, MPI_COMM_WORLD + use mpi_func, only : mpi_bcast, MASTER_RANK, MPI_COMM_WORLD, MPI_SHORTINT #endif use hashFunctions_func, only : FNV_1 use dictionary_class, only : dictionary @@ -383,7 +383,7 @@ subroutine init(self, dict) ! Broadcast seed to all processes #ifdef MPI - call MPI_Bcast(seed_temp, 1, MPI_INT, MASTER_RANK, MPI_COMM_WORLD) + call mpi_bcast(seed_temp, 1, MPI_SHORTINT, MASTER_RANK, MPI_COMM_WORLD) #endif seed = seed_temp diff --git a/SharedModules/errors_mod.f90 b/SharedModules/errors_mod.f90 index f924e7ab1..09b81708b 100644 --- a/SharedModules/errors_mod.f90 +++ b/SharedModules/errors_mod.f90 @@ -9,7 +9,6 @@ module errors_mod use numPrecision use universalVariables, only : MAX_COL - !use mpi_func, only : getMPIRank implicit none @@ -37,9 +36,6 @@ subroutine fatalError(where, why) ! Upper frame write(error_unit, *) repeat('<>', MAX_COL / 2) -! #ifdef MPI -! write(error_unit, *) 'Process rank: ', getMPIRank() -! #endif write(error_unit, *) 'Fatal has occurred in:' write(error_unit, *) where, new_line('') write(error_unit, *) 'Because:' diff --git a/SharedModules/mpi_func.f90 b/SharedModules/mpi_func.f90 index da17759c7..33ab7152c 100644 --- a/SharedModules/mpi_func.f90 +++ b/SharedModules/mpi_func.f90 @@ -4,6 +4,7 @@ module mpi_func use mpi_f08 #endif use errors_mod, only : fatalError + use particle_class, only : particleStateData implicit none @@ -11,30 +12,13 @@ module mpi_func integer(shortInt), private :: rank = 0 integer(shortInt), parameter :: MASTER_RANK = 0 - !! Public type that replicates exactly particleState - !! - !! It is necessary for load balancing in the dungeon: particles have to be - !! transferred betwen processes, and MPI doesn't allow to transfer types with - !! type-bound procedures - type, public :: particleStateDummy - real(defReal) :: wgt - real(defReal),dimension(3) :: r - real(defReal),dimension(3) :: dir - real(defReal) :: E - integer(shortInt) :: G - logical(defBool) :: isMG - integer(shortInt) :: type - real(defReal) :: time - integer(shortInt) :: matIdx - integer(shortInt) :: cellIdx - integer(shortInt) :: uniqueID - integer(shortInt) :: collisionN - integer(shortInt) :: broodID - end type particleStateDummy - - !! Common MPI types + !! MPI types #ifdef MPI - type(MPI_Datatype) :: MPI_PARTICLE_STATE + type(MPI_Datatype) :: MPI_DEFREAL + type(MPI_Datatype) :: MPI_SHORTINT + type(MPI_Datatype) :: MPI_LONGINT + type(MPI_Datatype) :: MPI_DEFBOOL + type(MPI_Datatype) :: MPI_PARTICLE_STATE #endif contains @@ -46,7 +30,7 @@ module mpi_func !! subroutine mpiInit() #ifdef MPI - integer(shortInt) :: ierr + integer(shortInt) :: ierr call mpi_init(ierr) @@ -79,12 +63,28 @@ end subroutine mpiFinalise !! subroutine mpiInitTypes() #ifdef MPI - integer(shortInt) :: ierr, stateSize - type(particleStateDummy) :: state + integer(shortInt) :: ierr, stateSize + type(particleStateData) :: state integer(kind = MPI_ADDRESS_KIND), dimension(:), allocatable :: displacements integer(shortInt), dimension(:), allocatable :: blockLengths type(MPI_Datatype), dimension(:), allocatable :: types + ! Define MPI type for DEFREAL + call mpi_type_create_f90_real(precision(1.0_defReal), range(1.0_defReal), MPI_DEFREAL, ierr) + call mpi_type_commit(MPI_DEFREAL, ierr) + + ! Define MPI type for SHORTINT + call mpi_type_create_f90_integer(range(1_shortInt), MPI_SHORTINT, ierr) + call mpi_type_commit(MPI_SHORTINT, ierr) + + ! Define MPI type for LONGINT + call mpi_type_create_f90_integer(range(1_longInt), MPI_LONGINT, ierr) + call mpi_type_commit(MPI_LONGINT, ierr) + + ! Define MPI type for DEFBOOL + call MPI_type_contiguous(defBool, MPI_BYTE, MPI_DEFBOOL, ierr) + call MPI_type_commit(MPI_DEFBOOL, ierr) + ! Define MPI type for particleState ! Note that particleState has stateSize = 13 attributes; if an attribute is ! added to particleState, it had to be added here too @@ -93,8 +93,9 @@ subroutine mpiInitTypes() ! Create arrays with dimension and type of each property of particleStateDummy blockLengths = (/1, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1/) - types = (/MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_LOGICAL, MPI_INT, & - MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_INT/) + types = (/MPI_DEFREAL, MPI_DEFREAL, MPI_DEFREAL, MPI_DEFREAL, MPI_SHORTINT, & + MPI_DEFBOOL, MPI_SHORTINT, MPI_DEFREAL, MPI_SHORTINT, MPI_SHORTINT, & + MPI_SHORTINT, MPI_SHORTINT, MPI_SHORTINT/) ! Create array of memory byte displacements call mpi_get_address(state % wgt, displacements(1), ierr) @@ -115,6 +116,7 @@ subroutine mpiInitTypes() ! Define new type call mpi_type_create_struct(stateSize, blockLengths, displacements, types, MPI_PARTICLE_STATE, ierr) call mpi_type_commit(MPI_PARTICLE_STATE, ierr) + #endif end subroutine mpiInitTypes diff --git a/Tallies/scoreMemory_class.f90 b/Tallies/scoreMemory_class.f90 index bbe88a453..bf67940c5 100644 --- a/Tallies/scoreMemory_class.f90 +++ b/Tallies/scoreMemory_class.f90 @@ -2,8 +2,8 @@ module scoreMemory_class use numPrecision #ifdef MPI - use mpi_func, only : mpi_reduce, MPI_SUM, MPI_DOUBLE, MPI_COMM_WORLD, & - MASTER_RANK, isMPIMaster, MPI_INT + use mpi_func, only : mpi_reduce, isMPIMaster, MPI_SUM, & + MPI_COMM_WORLD, MASTER_RANK, MPI_DEFREAL, MPI_SHORTINT #endif use universalVariables, only : array_pad use genericProcedures, only : fatalError, numToChar @@ -449,7 +449,7 @@ subroutine collectDistributed(self) ! Reduce the batch count ! Note we need to use size 1 arrays to fit the interface of mpi_reduce, which expects ! to be given arrays - call mpi_reduce(self % batchN, buffer, 1, MPI_INT, MPI_SUM, MASTER_RANK, MPI_COMM_WORLD, error) + call mpi_reduce(self % batchN, buffer, 1, MPI_SHORTINT, MPI_SUM, MASTER_RANK, MPI_COMM_WORLD, error) if (isMPIMaster()) then self % batchN = buffer else @@ -510,7 +510,7 @@ subroutine reduceArray(data, buffer) call mpi_reduce(data(start : start + chunk - 1), & buffer(start : start + chunk - 1), & int(chunk, shortInt), & - MPI_DOUBLE, & + MPI_DEFREAL, & MPI_SUM, & MASTER_RANK, & MPI_COMM_WORLD, & diff --git a/docs/Running.rst b/docs/Running.rst index 8913299dd..6aeb21655 100644 --- a/docs/Running.rst +++ b/docs/Running.rst @@ -18,7 +18,7 @@ and to the input file:: .. admonition:: Options - OpenMPI + OpenMP Specifies the number ```` of OpenMP threads to be used:: --omp @@ -34,4 +34,6 @@ and to the input file:: /mpirun -np /scone.out + OpenMP and MPI can be run together too:: + /mpirun -np /scone.out --omp