diff --git a/src/blocks.F90 b/src/blocks.F90 index 0a47ec4..7222a0b 100644 --- a/src/blocks.F90 +++ b/src/blocks.F90 @@ -45,6 +45,9 @@ module blocks ! timer indices ! integer, save :: imi, ima, imu, imp, imq, imr, imd +#ifdef DEBUG + integer, save :: imc +#endif /* DEBUG */ #endif /* PROFILE */ ! MODULE PARAMETERS: @@ -54,11 +57,13 @@ module blocks ! nsides - the number of sides along each direction (2); ! nfaces - the number of faces at each side (2 for 2D, 4 for 3D); ! nchildren - the number of child blocks for each block (4 for 2D, 8 for 3D); +! mfaces - the number of faces in block (8 for 2D, 24 for 3D); ! integer(kind=4), parameter :: ndims = NDIMS integer(kind=4), parameter :: nsides = 2 integer(kind=4), parameter :: nfaces = 2**(ndims - 1) integer(kind=4), parameter :: nchildren = 2**ndims + integer(kind=4), parameter :: mfaces = nsides * nfaces * ndims ! MODULE VARIABLES: ! ================ @@ -271,11 +276,15 @@ module blocks public :: refine_block, derefine_block public :: set_last_id, get_last_id, get_mblocks, get_dblocks, get_nleafs public :: set_blocks_update + public :: change_blocks_process public :: set_neighbors_refine public :: metablock_set_id, metablock_set_process, metablock_set_level public :: metablock_set_configuration, metablock_set_refinement public :: metablock_set_position, metablock_set_coordinates public :: metablock_set_bounds, metablock_set_leaf, metablock_unset_leaf +#ifdef DEBUG + public :: check_neighbors +#endif /* DEBUG */ !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! @@ -324,6 +333,9 @@ module blocks call set_timer('blocks:: data block deallocation', imq) call set_timer('blocks:: refine' , imr) call set_timer('blocks:: derefine' , imd) +#ifdef DEBUG + call set_timer('blocks:: check neighbors' , imc) +#endif /* DEBUG */ ! start accounting time for module initialization/finalization ! @@ -2152,6 +2164,58 @@ module blocks ! !=============================================================================== ! +! subroutine CHANGE_BLOCKS_PROCESS: +! -------------------------------- +! +! Subroutine switches meta blocks which belong to old process to the new one. +! +! Arguments: +! +! npold - the old process number; +! npnew - the new process number; +! +!=============================================================================== +! + subroutine change_blocks_process(npold, npnew) + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer, intent(in) :: npold, npnew + +! local pointers +! + type(block_meta), pointer :: pmeta +! +!------------------------------------------------------------------------------- +! +! associate the pointer with the first block on the meta block list +! + pmeta => list_meta + +! iterate over all blocks in the list +! + do while(associated(pmeta)) + +! if the meta block belongs to process npold, switch it to process npnew +! + if (pmeta%process == npold) pmeta%process = npnew + +! associate the pointer with the next block on the list +! + pmeta => pmeta%next + + end do ! meta blocks + +!------------------------------------------------------------------------------- +! + end subroutine change_blocks_process +! +!=============================================================================== +! ! subroutine METABLOCK_SET_ID: ! --------------------------- ! @@ -2658,6 +2722,64 @@ module blocks !------------------------------------------------------------------------------- ! end subroutine set_neighbors_refine +#ifdef DEBUG +! +!=============================================================================== +! +! subroutine CHECK_NEIGHBORS: +! -------------------------- +! +! Subroutine iterates over all blocks and checks if the pointers to their +! neighbors are consistent. +! +!=============================================================================== +! + subroutine check_neighbors() + +! local variables are not implicit by default +! + implicit none + +! local pointers +! + type(block_meta), pointer :: pmeta +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for the neighbor consistency check +! + call start_timer(imc) +#endif /* PROFILE */ + +! associate the pointer with the first block on the meta block list +! + pmeta => list_meta + +! iterate over all blocks in the list +! + do while(associated(pmeta)) + +! check the block neighbors +! + call check_block_neighbors(pmeta) + +! associate the pointer with the next block on the list +! + pmeta => pmeta%next + + end do ! meta blocks + +#ifdef PROFILE +! stop accounting time for the neighbor consistency check +! + call stop_timer(imc) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine check_neighbors +#endif /* DEBUG */ ! !=============================================================================== !! @@ -3064,442 +3186,595 @@ module blocks type(block_meta) , pointer, intent(inout) :: pmeta procedure(reset_neighbors_update), pointer, intent(in) :: pprocedure -! local pointers +! local saved variables ! - type(block_meta), pointer :: pneigh + logical, save :: first = .true. + +! local saved arrays +! + integer, dimension(mfaces,3) , save :: fidx + integer, dimension(mfaces,2,3), save :: eidx + integer, dimension(mfaces,2,3), save :: cidx + +! local variables +! + integer :: l ! !------------------------------------------------------------------------------- ! -#if NDIMS == 3 -!! 3D CASE: WALK AROUND THE CORNERS -!! +! prepare indices +! + if (first) then + +! inicialize indices +! + fidx(: ,:) = 0 + eidx(:,:,:) = 0 + cidx(:,:,:) = 0 + +! prepare indices to get proper face, edge and corner neighbors +! +#if NDIMS == 2 ! around (0,0,0) corner ! -! (1,1,1) - X = (1,2,1) - Y = (2,1,2) +! [1,1,1]:X:(1,2,1) <- Y = [2,1,2] ! - pneigh => pmeta%neigh(1,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx( 1, :) = (/ 1, 1, 1 /) + eidx( 1,1,:) = (/ 2, 1, 2 /) - pneigh => pneigh%neigh(2,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,1,1) - Z = (3,2,1) - X = (1,1,3) +! [2,1,1]:Y:(2,2,1) <- X = [1,1,2] ! - pneigh => pmeta%neigh(3,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,1,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (2,1,1) - Y = (2,2,1) - Z = (3,1,3) - X = [1,1,4] -! - pneigh => pmeta%neigh(2,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(3,1,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,1,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if + fidx( 2, :) = (/ 2, 1, 1 /) + eidx( 2,1,:) = (/ 1, 1, 2 /) ! around (0,1,0) corner ! -! (2,2,1) + Y = (2,1,1) - X = (1,1,1) +! [1,1,2]:X:(1,2,2) -> Y = [2,2,2] ! - pneigh => pmeta%neigh(2,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx( 3, :) = (/ 1, 1, 2 /) + eidx( 3,1,:) = (/ 2, 2, 2 /) - pneigh => pneigh%neigh(1,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,1,3) - Z = (3,2,3) + Y = (2,2,3) +! [2,2,1]:Y:(2,1,2) <- X = [1,1,1] ! - pneigh => pmeta%neigh(3,1,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (1,1,2) - X = (1,2,2) - Z = (3,1,4) + Y = [2,2,4] -! - pneigh => pmeta%neigh(1,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(3,1,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,2,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if + fidx( 4, :) = (/ 2, 2, 1 /) + eidx( 4,1,:) = (/ 1, 1, 1 /) ! around (1,1,0) corner ! -! (1,2,2) + X = (1,1,2) + Y = (2,2,1) +! [1,2,2]:X:(1,1,2) -> Y = [2,2,1] ! - pneigh => pmeta%neigh(1,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx( 5, :) = (/ 1, 2, 2 /) + eidx( 5,1,:) = (/ 2, 2, 1 /) - pneigh => pneigh%neigh(2,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,1,4) - Z = (3,2,4) + X = (1,2,4) +! [2,2,2]:Y:(2,1,2) -> X = [1,2,1] ! - pneigh => pmeta%neigh(3,1,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,2,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (2,2,2) + Y = (2,1,2) - Z = (3,1,2) + X = [1,2,3] -! - pneigh => pmeta%neigh(2,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(3,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if + fidx( 6, :) = (/ 2, 2, 2 /) + eidx( 6,1,:) = (/ 1, 2, 1 /) ! around (1,0,0) corner ! -! (2,1,2) - Y = (2,2,2) + X = (1,2,2) +! [1,2,1]:X:(1,1,1) <- Y = [2,1,1] ! - pneigh => pmeta%neigh(2,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx( 7, :) = (/ 1, 2, 1 /) + eidx( 7,1,:) = (/ 2, 1, 1 /) - pneigh => pneigh%neigh(1,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,1,2) - Z = (3,2,2) - Y = (2,1,4) +! [2,1,2]:Y:(2,2,2) -> X = [1,2,2] ! - pneigh => pmeta%neigh(3,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,1,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (1,2,1) + X = (1,1,1) - Z = (3,1,1) - Y = [2,1,3] + fidx( 8, :) = (/ 2, 1, 2 /) + eidx( 8,1,:) = (/ 1, 2, 2 /) +#endif /* NDIMS == 2 */ +#if NDIMS == 3 +! around (0,0,0) corner ! - pneigh => pmeta%neigh(1,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) +! [1,1,1]:X:(1,2,1) <- Y = [2,1,2] <- Z = [3,1,4] +! [1,1,1]:X:(1,2,1) <- Z = [3,1,2] <- Y = [2,1,4] +! + fidx( 1, :) = (/ 1, 1, 1 /) + eidx( 1,1,:) = (/ 2, 1, 2 /) + eidx( 1,2,:) = (/ 3, 1, 2 /) + cidx( 1,1,:) = (/ 3, 1, 4 /) + cidx( 1,2,:) = (/ 2, 1, 4 /) - pneigh => pneigh%neigh(3,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) +! [2,1,1]:Y:(2,2,1) <- Z = [3,1,3] <- X = [1,1,4] +! [2,1,1]:Y:(2,2,1) <- X = [1,1,2] <- Z = [3,1,4] +! + fidx( 2, :) = (/ 2, 1, 1 /) + eidx( 2,1,:) = (/ 3, 1, 3 /) + eidx( 2,2,:) = (/ 1, 1, 2 /) + cidx( 2,1,:) = (/ 1, 1, 4 /) + cidx( 2,2,:) = (/ 3, 1, 4 /) - pneigh => pneigh%neigh(2,1,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if +! [3,1,1]:Z:(3,2,1) <- X = [1,1,3] <- Y = [2,1,4] +! [3,1,1]:Z:(3,2,1) <- Y = [2,1,3] <- X = [1,1,4] +! + fidx( 3, :) = (/ 3, 1, 1 /) + eidx( 3,1,:) = (/ 1, 1, 3 /) + eidx( 3,2,:) = (/ 2, 1, 3 /) + cidx( 3,1,:) = (/ 2, 1, 4 /) + cidx( 3,2,:) = (/ 1, 1, 4 /) + +! around (0,1,0) corner +! +! [1,1,2]:X:(1,2,2) -> Y = [2,2,2] <- Z = [3,1,2] +! [1,1,2]:X:(1,2,2) <- Z = [3,1,4] -> Y = [2,2,4] +! + fidx( 4, :) = (/ 1, 1, 2 /) + eidx( 4,1,:) = (/ 2, 2, 2 /) + eidx( 4,2,:) = (/ 3, 1, 4 /) + cidx( 4,1,:) = (/ 3, 1, 2 /) + cidx( 4,2,:) = (/ 2, 2, 4 /) + +! [2,2,1]:Y:(2,1,2) <- Z = [3,1,1] <- X = [1,1,3] +! [2,2,1]:Y:(2,1,2) <- X = [1,1,1] <- Z = [3,1,2] +! + fidx( 5, :) = (/ 2, 2, 1 /) + eidx( 5,1,:) = (/ 3, 1, 1 /) + eidx( 5,2,:) = (/ 1, 1, 1 /) + cidx( 5,1,:) = (/ 1, 1, 3 /) + cidx( 5,2,:) = (/ 3, 1, 2 /) + +! [3,1,3]:Z:(3,2,3) <- X = [1,1,4] -> Y = [2,2,4] +! [3,1,3]:Z:(3,2,3) -> Y = [2,2,3] <- X = [1,1,3] +! + fidx( 6, :) = (/ 3, 1, 3 /) + eidx( 6,1,:) = (/ 1, 1, 4 /) + eidx( 6,2,:) = (/ 2, 2, 3 /) + cidx( 6,1,:) = (/ 2, 2, 4 /) + cidx( 6,2,:) = (/ 1, 1, 3 /) + +! around (1,1,0) corner +! +! [1,2,2]:X:(1,1,2) -> Y = [2,2,1] <- Z = [3,1,1] +! [1,2,2]:X:(1,1,2) <- Z = [3,1,3] -> Y = [2,2,3] +! + fidx( 7, :) = (/ 1, 2, 2 /) + eidx( 7,1,:) = (/ 2, 2, 1 /) + eidx( 7,2,:) = (/ 3, 1, 3 /) + cidx( 7,1,:) = (/ 3, 1, 1 /) + cidx( 7,2,:) = (/ 2, 2, 3 /) + +! [2,2,2]:Y:(2,1,2) <- Z = [3,1,2] -> X = [1,2,3] +! [2,2,2]:Y:(2,1,2) -> X = [1,2,1] <- Z = [3,1,1] +! + fidx( 8, :) = (/ 2, 2, 2 /) + eidx( 8,1,:) = (/ 3, 1, 2 /) + eidx( 8,2,:) = (/ 1, 2, 1 /) + cidx( 8,1,:) = (/ 1, 2, 3 /) + cidx( 8,2,:) = (/ 3, 1, 1 /) + +! [3,1,4]:Z:(3,2,4) -> X = [1,2,4] -> Y = [2,2,3] +! [3,1,4]:Z:(3,2,4) -> Y = [2,2,4] -> X = [1,2,3] +! + fidx( 9, :) = (/ 3, 1, 4 /) + eidx( 9,1,:) = (/ 1, 2, 4 /) + eidx( 9,2,:) = (/ 2, 2, 4 /) + cidx( 9,1,:) = (/ 2, 2, 3 /) + cidx( 9,2,:) = (/ 1, 2, 3 /) + +! around (1,0,0) corner +! +! [1,2,1]:X:(1,1,1) <- Y = [2,1,1] <- Z = [3,1,3] +! [1,2,1]:X:(1,1,1) <- Z = [3,1,1] <- Y = [2,1,3] +! + fidx(10, :) = (/ 1, 2, 1 /) + eidx(10,1,:) = (/ 2, 1, 1 /) + eidx(10,2,:) = (/ 3, 1, 1 /) + cidx(10,1,:) = (/ 3, 1, 3 /) + cidx(10,2,:) = (/ 2, 1, 3 /) + +! [2,1,2]:Y:(2,2,2) <- Z = [3,1,4] -> X = [1,2,4] +! [2,1,2]:Y:(2,2,2) -> X = [1,2,2] <- Z = [3,1,3] +! + fidx(11, :) = (/ 2, 1, 2 /) + eidx(11,1,:) = (/ 3, 1, 4 /) + eidx(11,2,:) = (/ 1, 2, 2 /) + cidx(11,1,:) = (/ 1, 2, 4 /) + cidx(11,2,:) = (/ 3, 1, 3 /) + +! [3,1,2]:Z:(3,2,2) -> X = [1,2,3] <- Y = [2,1,3] +! [3,1,2]:Z:(3,2,2) <- Y = [2,1,4] -> X = [1,2,4] +! + fidx(12, :) = (/ 3, 1, 2 /) + eidx(12,1,:) = (/ 1, 2, 3 /) + eidx(12,2,:) = (/ 2, 1, 4 /) + cidx(12,1,:) = (/ 2, 1, 3 /) + cidx(12,2,:) = (/ 1, 2, 4 /) ! around (0,0,1) corner ! -! (1,1,3) - X = (1,2,3) + Z = (3,2,2) +! [1,1,3]:X:(1,2,3) <- Y = [2,1,4] -> Z = [3,2,4] +! [1,1,3]:X:(1,2,3) -> Z = [3,2,2] <- Y = [2,1,2] ! - pneigh => pmeta%neigh(1,1,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(13, :) = (/ 1, 1, 3 /) + eidx(13,1,:) = (/ 2, 1, 4 /) + eidx(13,2,:) = (/ 3, 2, 2 /) + cidx(13,1,:) = (/ 3, 2, 4 /) + cidx(13,2,:) = (/ 2, 1, 2 /) - pneigh => pneigh%neigh(3,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (2,1,3) - Y = (2,2,3) - X = (1,1,4) +! [2,1,3]:Y:(2,2,3) -> Z = [3,2,3] <- X = [1,1,2] +! [2,1,3]:Y:(2,2,3) <- X = [1,1,4] -> Z = [3,2,1] ! - pneigh => pmeta%neigh(2,1,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(14, :) = (/ 2, 1, 3 /) + eidx(14,1,:) = (/ 3, 2, 3 /) + eidx(14,2,:) = (/ 1, 1, 4 /) + cidx(14,1,:) = (/ 1, 1, 2 /) + cidx(14,2,:) = (/ 3, 2, 1 /) - pneigh => pneigh%neigh(1,1,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,2,1) + Z = (3,1,1) - Y = (2,1,1) - X = [1,1,2] +! [3,2,1]:Z:(3,1,1) <- X = [1,1,1] <- Y = [2,1,2] +! [3,2,1]:Z:(3,1,1) <- Y = [2,1,1] <- X = [1,1,2] ! - pneigh => pmeta%neigh(3,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if + fidx(15, :) = (/ 3, 2, 1 /) + eidx(15,1,:) = (/ 1, 1, 1 /) + eidx(15,2,:) = (/ 2, 1, 1 /) + cidx(15,1,:) = (/ 2, 1, 2 /) + cidx(15,2,:) = (/ 1, 1, 2 /) ! around (0,1,1) corner ! -! (2,2,3) + Y = (2,1,3) + Z = (3,2,1) +! [1,1,4]:X:(1,2,4) -> Y = [2,2,4] -> Z = [3,2,2] +! [1,1,4]:X:(1,2,4) -> Z = [3,2,4] -> Y = [2,2,2] ! - pneigh => pmeta%neigh(2,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(16, :) = (/ 1, 1, 4 /) + eidx(16,1,:) = (/ 2, 2, 4 /) + eidx(16,2,:) = (/ 3, 2, 4 /) + cidx(16,1,:) = (/ 3, 2, 2 /) + cidx(16,2,:) = (/ 2, 2, 2 /) - pneigh => pneigh%neigh(3,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (1,1,4) - X = (1,2,4) + Y = (2,2,4) +! [2,2,3]:Y:(2,1,3) -> Z = [3,2,1] <- X = [1,1,1] +! [2,2,3]:Y:(2,1,3) <- X = [1,1,3] -> Z = [3,2,2] ! - pneigh => pmeta%neigh(1,1,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(17, :) = (/ 2, 2, 3 /) + eidx(17,1,:) = (/ 3, 2, 1 /) + eidx(17,2,:) = (/ 1, 1, 3 /) + cidx(17,1,:) = (/ 1, 1, 1 /) + cidx(17,2,:) = (/ 3, 2, 2 /) - pneigh => pneigh%neigh(2,2,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,2,3) + Z = (3,1,3) - X = (1,1,2) + Z = [2,2,2] +! [3,2,3]:Z:(3,1,3) <- X = [1,1,2] -> Y = [2,2,2] +! [3,2,3]:Z:(3,1,3) -> Y = [2,2,1] <- X = [1,1,1] ! - pneigh => pmeta%neigh(3,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if + fidx(18, :) = (/ 3, 2, 3 /) + eidx(18,1,:) = (/ 1, 1, 2 /) + eidx(18,2,:) = (/ 2, 2, 1 /) + cidx(18,1,:) = (/ 2, 2, 2 /) + cidx(18,2,:) = (/ 1, 1, 1 /) ! around (1,1,1) corner ! -! (1,2,4) + X = (1,1,4) + Z = (3,2,3) +! [1,2,4]:X:(1,1,4) -> Y = [2,2,3] -> Z = [3,2,1] +! [1,2,4]:X:(1,1,4) -> Z = [3,2,3] -> Y = [2,2,1] ! - pneigh => pmeta%neigh(1,2,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(19, :) = (/ 1, 2, 4 /) + eidx(19,1,:) = (/ 2, 2, 3 /) + eidx(19,2,:) = (/ 3, 2, 3 /) + cidx(19,1,:) = (/ 3, 2, 1 /) + cidx(19,2,:) = (/ 2, 2, 1 /) - pneigh => pneigh%neigh(3,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (2,2,4) + Y = (2,1,4) + X = (1,2,3) +! [2,2,4]:Y:(2,1,4) -> Z = [3,2,2] -> X = [1,2,1] +! [2,2,4]:Y:(2,1,4) -> X = [1,2,3] -> Z = [3,2,1] ! - pneigh => pmeta%neigh(2,2,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(20, :) = (/ 2, 2, 4 /) + eidx(20,1,:) = (/ 3, 2, 2 /) + eidx(20,2,:) = (/ 1, 2, 3 /) + cidx(20,1,:) = (/ 1, 2, 1 /) + cidx(20,2,:) = (/ 3, 2, 1 /) - pneigh => pneigh%neigh(1,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,2,4) + Z = (3,1,4) + Y = (2,2,2) + X = [1,2,1] +! [3,2,4]:Z:(3,1,4) -> X = [1,2,2] -> Y = [2,2,1] +! [3,2,4]:Z:(3,1,4) -> Y = [2,2,2] -> X = [1,2,1] ! - pneigh => pmeta%neigh(3,2,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(1,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if + fidx(21, :) = (/ 3, 2, 4 /) + eidx(21,1,:) = (/ 1, 2, 2 /) + eidx(21,2,:) = (/ 2, 2, 2 /) + cidx(21,1,:) = (/ 2, 2, 1 /) + cidx(21,2,:) = (/ 1, 2, 1 /) ! around (1,0,1) corner ! -! (2,1,4) - Y = (2,2,4) + Z = (3,2,4) +! [1,2,3]:X:(1,1,3) <- Y = [2,1,3] -> Z = [3,2,3] +! [1,2,3]:X:(1,1,3) -> Z = [3,2,1] <- Y = [2,1,1] ! - pneigh => pmeta%neigh(2,1,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(22, :) = (/ 1, 2, 3 /) + eidx(22,1,:) = (/ 2, 1, 3 /) + eidx(22,2,:) = (/ 3, 2, 1 /) + cidx(22,1,:) = (/ 3, 2, 3 /) + cidx(22,2,:) = (/ 2, 1, 1 /) - pneigh => pneigh%neigh(3,2,4)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (3,2,2) + Z = (3,1,2) + X = (1,2,1) +! [2,1,4]:Y:(2,2,4) -> Z = [3,2,4] -> X = [1,2,2] +! [2,1,4]:Y:(2,2,4) -> X = [1,2,4] -> Z = [3,2,3] ! - pneigh => pmeta%neigh(3,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) + fidx(23, :) = (/ 2, 1, 4 /) + eidx(23,1,:) = (/ 3, 2, 4 /) + eidx(23,2,:) = (/ 1, 2, 4 /) + cidx(23,1,:) = (/ 1, 2, 2 /) + cidx(23,2,:) = (/ 3, 2, 3 /) - pneigh => pneigh%neigh(1,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (1,2,3) + X = (1,1,3) - Y = (2,1,3) + Z = [3,2,3] +! [3,2,2]:Z:(3,1,2) -> X = [1,2,1] <- Y = [2,1,3] +! [3,2,2]:Z:(3,1,2) <- Y = [2,1,2] -> X = [1,2,4] ! - pneigh => pmeta%neigh(1,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,1,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(3,2,3)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - end if - -#else /* NDIMS == 3 */ -!! 2D CASE -!! -! (0,0)-(0,½) edge and (0,0) corner -! - pneigh => pmeta%neigh(1,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (0,½)-(0,1) edge and (0,1) corner -! - pneigh => pmeta%neigh(1,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (1,0)-(1,½) edge and (1,0) corner -! - pneigh => pmeta%neigh(1,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (1,½)-(1,1) edge and (1,1) corner -! - pneigh => pmeta%neigh(1,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - - pneigh => pneigh%neigh(2,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - end if - -! (0,0)-(½,0) edge -! - pneigh => pmeta%neigh(2,1,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - -! (½,0)-(1,0) edge -! - pneigh => pmeta%neigh(2,1,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - -! (0,1)-(½,1) edge -! - pneigh => pmeta%neigh(2,2,1)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if - -! (½,1)-(1,1) edge -! - pneigh => pmeta%neigh(2,2,2)%ptr - if (associated(pneigh)) then - call pprocedure(pmeta, pneigh) - end if + fidx(24, :) = (/ 3, 2, 2 /) + eidx(24,1,:) = (/ 1, 2, 1 /) + eidx(24,2,:) = (/ 2, 1, 2 /) + cidx(24,1,:) = (/ 2, 1, 3 /) + cidx(24,2,:) = (/ 1, 2, 3 /) #endif /* NDIMS == 3 */ +! reset the first time execution flag +! + first = .false. + + end if + +! iterate over all block faces (or edges in the 2D case) +! + do l = 1, mfaces + call iterate_over_face(pmeta, pprocedure & + , fidx(l,:), eidx(l,:,:), cidx(l,:,:)) + end do + !------------------------------------------------------------------------------- ! end subroutine iterate_over_neighbors +! +!=============================================================================== +! +! subroutine ITERATE_OVER_FACE: +! ---------------------------- +! +! Subroutine iterates over all neighbors, edges and corners linked to +! the input meta block and executes a subroutine provided by the pointer. +! +! Arguments: +! +! pmeta - a pointer to the meta block which neighbors are iterated over; +! pproc - a pointer to the subroutine called with each pair (pmeta, pneigh); +! fidx - the index of face to process; +! eidx - the indices of faces connected with edges; +! cidx - the indices of faces connected with corners; +! +!=============================================================================== +! + subroutine iterate_over_face(pmeta, pprocedure, fidx, eidx, cidx) + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + type(block_meta) , pointer, intent(inout) :: pmeta + procedure(reset_neighbors_update), pointer, intent(in) :: pprocedure + integer, dimension(3) , intent(in) :: fidx + integer, dimension(2,3) , intent(in) :: eidx + integer, dimension(2,3) , intent(in) :: cidx + +! local pointers +! + type(block_meta), pointer :: pneigh, pedge, pcorner +! +!------------------------------------------------------------------------------- +! +! associate a pointer with the neighbor +! + pneigh => pmeta%neigh(fidx(1),fidx(2),fidx(3))%ptr + +! check if the neighbors is associated +! + if (associated(pneigh)) then + +! call the procedure for the face neighbor +! + call pprocedure(pmeta, pneigh) + +! associate a pointer with the first edge +! + pedge => pneigh%neigh(eidx(1,1),eidx(1,2),eidx(1,3))%ptr + +! check if the edge pointer is associated +! + if (associated(pedge)) then + +! call the procedure for the edge neighbor +! + call pprocedure(pmeta, pedge) + +#if NDIMS == 3 +! associate a pointer with the first corner +! + pcorner => pedge%neigh(cidx(1,1),cidx(1,2),cidx(1,3))%ptr + +! call the procedure for the corner neighbor if it is associated +! + if (associated(pcorner)) call pprocedure(pmeta, pcorner) +#endif /* NDIMS == 3 */ + + end if ! pedge associated + +#if NDIMS == 3 +! associate a pointer with the second edge +! + pedge => pneigh%neigh(eidx(2,1),eidx(2,2),eidx(2,3))%ptr + +! check if the edge pointer is associated +! + if (associated(pedge)) then + +! call the procedure for the edge neighbor +! + call pprocedure(pmeta, pedge) + +! associate a pointer with the second corner +! + pcorner => pedge%neigh(cidx(2,1),cidx(2,2),cidx(2,3))%ptr + +! call the procedure for the corner neighbor if it is associated +! + if (associated(pcorner)) call pprocedure(pmeta, pcorner) + + end if ! pedge associated +#endif /* NDIMS == 3 */ + + end if ! pneigh associated + +!------------------------------------------------------------------------------- +! + end subroutine iterate_over_face +#ifdef DEBUG +! +!=============================================================================== +! +! subroutine CHECK_BLOCK_NEIGHBORS: +! -------------------------------- +! +! Subroutine iterates over all neighbor blocks and checks if their pointers +! are consistent, i.e. if their corresponding neighbor pointers point to +! the current block. +! +! Arguments: +! +! pmeta - a pointer to the meta block; +! +!=============================================================================== +! + subroutine check_block_neighbors(pmeta) + +! import external procedures and variables +! + use error , only : print_warning + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + type(block_meta), pointer, intent(in) :: pmeta + +! local pointers +! + type(block_meta), pointer :: pneigh, pnneigh + +! local variables +! + integer :: i, j, k, l, m +! +!------------------------------------------------------------------------------- +! +! return if the block is not leaf +! + if (.not. pmeta%leaf) return + +! iterate over all face neighbors +! + do i = 1, ndims + do j = 1, nsides + m = 3 - j + do k = 1, nfaces + +! assign pointer with the neighbor +! + pneigh => pmeta%neigh(i,j,k)%ptr + +! check if the neighbor is associated +! + if (associated(pneigh)) then + +! check neighbors on the same levels +! + if (pmeta%level == pneigh%level) then + +! assign pointer to the neighbor of the neighbor pointing to the current meta +! block +! + pnneigh => pneigh%neigh(i,m,k)%ptr + +! check if it is associated +! + if (associated(pnneigh)) then + +! check if the pointer of the neighbor points to the current meta block +! + if (pmeta%id /= pnneigh%id) then + +! print warning, since the blocks differ +! + call print_warning("blocks::check_block_neighbors" & + , "Inconsistent same level neighbors!") + print *, 'metablock: ', pmeta%id, pnneigh%id + print *, 'neighbor : ', pneigh%id + print *, 'index : ', i, j, k + + end if + + else ! pnneigh associated + +! print warning, since the pointer should be associated +! + call print_warning("blocks::check_block_neighbors" & + , "Same level neighbor not associated!") + print *, 'metablock: ', pmeta%id, pnneigh%id + print *, 'neighbor : ', pneigh%id + print *, 'index : ', i, j, k + + end if ! pnneigh associated + + end if ! the same levels + +! check neighbors on the level higher than the meta block's level; it also +! covers the other way around, since we iterate over all neighbor faces +! + if (pmeta%level < pneigh%level) then + +! iterate over whole face of the corresponding neighbor +! + do l = 1, nfaces + +! assign pointer to the corresponding neighbor of the neighbor +! + pnneigh => pneigh%neigh(i,m,l)%ptr + +! check if it is associated +! + if (associated(pnneigh)) then + +! check if the pointer of the neighbor points to the current meta block +! + if (pmeta%id /= pnneigh%id) then + +! print warning, since the blocks differ +! + call print_warning("blocks::check_block_neighbors" & + , "Inconsistent different level neighbors!") + print *, 'metablock: ', pmeta%id, pnneigh%id + print *, 'neighbor : ', pneigh%id + print *, 'index : ', i, j, k, l + + end if + + else ! pnneigh associated + +! print warning, since the pointer should be associated +! + call print_warning("blocks::check_block_neighbors" & + , "Different level neighbor not associated!") + print *, 'metablock: ', pmeta%id, pnneigh%id + print *, 'neighbor : ', pneigh%id + print *, 'index : ', i, j, k, l + + end if ! pnneigh associated + + end do ! l = 1, nfaces + + end if ! pmeta's level < pneigh's level + + end if ! pneigh associated + + end do ! nfaces + end do ! nsides + end do ! ndims + +!------------------------------------------------------------------------------- +! + end subroutine check_block_neighbors +#endif /* DEBUG */ !=============================================================================== ! diff --git a/src/boundaries.F90 b/src/boundaries.F90 index 5b895bd..0a3f837 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -42,17 +42,18 @@ module boundaries #ifdef PROFILE ! timer indices ! - integer , save :: imi, imv, imf, ims, imc, imr, imp + integer , save :: imi, imv, imf, ims, imc, imr, imp #endif /* PROFILE */ -! module parameters for the boundary update order and boundary type +! parameters corresponding to the boundary type ! - character(len = 32), save :: xlbndry = "periodic" - character(len = 32), save :: xubndry = "periodic" - character(len = 32), save :: ylbndry = "periodic" - character(len = 32), save :: yubndry = "periodic" - character(len = 32), save :: zlbndry = "periodic" - character(len = 32), save :: zubndry = "periodic" + integer, parameter :: bnd_periodic = 0 + integer, parameter :: bnd_open = 1 + integer, parameter :: bnd_reflective = 2 + +! variable to store boundary type flags +! + integer, dimension(3,2), save :: bnd_type = bnd_periodic ! by default everything is private ! @@ -62,7 +63,7 @@ module boundaries ! public :: initialize_boundaries, finalize_boundaries public :: boundary_variables, boundary_fluxes - public :: xlbndry, ylbndry, zlbndry, xubndry, yubndry, zubndry + public :: bnd_type, bnd_periodic !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! @@ -105,6 +106,15 @@ module boundaries ! logical, intent(in) :: verbose integer, intent(inout) :: iret + +! module parameters for the boundary update order and boundary type +! + character(len = 32) :: xlbndry = "periodic" + character(len = 32) :: xubndry = "periodic" + character(len = 32) :: ylbndry = "periodic" + character(len = 32) :: yubndry = "periodic" + character(len = 32) :: zlbndry = "periodic" + character(len = 32) :: zubndry = "periodic" ! !------------------------------------------------------------------------------- ! @@ -133,6 +143,62 @@ module boundaries call get_parameter_string ("zlbndry" , zlbndry) call get_parameter_string ("zubndry" , zubndry) +! fill the boundary type flags +! + select case(xlbndry) + case("open") + bnd_type(1,1) = bnd_open + case("reflective", "reflecting", "reflect") + bnd_type(1,1) = bnd_reflective + case default + bnd_type(1,1) = bnd_periodic + end select + + select case(xubndry) + case("open") + bnd_type(1,2) = bnd_open + case("reflective", "reflecting", "reflect") + bnd_type(1,2) = bnd_reflective + case default + bnd_type(1,2) = bnd_periodic + end select + + select case(ylbndry) + case("open") + bnd_type(2,1) = bnd_open + case("reflective", "reflecting", "reflect") + bnd_type(2,1) = bnd_reflective + case default + bnd_type(2,1) = bnd_periodic + end select + + select case(yubndry) + case("open") + bnd_type(2,2) = bnd_open + case("reflective", "reflecting", "reflect") + bnd_type(2,2) = bnd_reflective + case default + bnd_type(2,2) = bnd_periodic + end select + + select case(zlbndry) + case("open") + bnd_type(3,1) = bnd_open + case("reflective", "reflecting", "reflect") + bnd_type(3,1) = bnd_reflective + case default + bnd_type(3,1) = bnd_periodic + end select + + select case(zubndry) + case("open") + bnd_type(3,2) = bnd_open + case("reflective", "reflecting", "reflect") + bnd_type(3,2) = bnd_reflective + case default + bnd_type(3,2) = bnd_periodic + end select + ! print information about the boundary conditions ! if (verbose) then @@ -148,52 +214,6 @@ module boundaries end if -#ifdef MPI -! change the internal boundaries to "exchange" type for the MPI update -! - if (pdims(1) > 1) then - if (periodic(1)) then - xlbndry = "exchange" - xubndry = "exchange" - else - if (pcoords(1) > 0 ) then - xlbndry = "exchange" - end if - if (pcoords(1) < pdims(1)-1) then - xubndry = "exchange" - end if - end if - end if - - if (pdims(2) > 1) then - if (periodic(2)) then - ylbndry = "exchange" - yubndry = "exchange" - else - if (pcoords(2) > 0 ) then - ylbndry = "exchange" - end if - if (pcoords(2) < pdims(2)-1) then - yubndry = "exchange" - end if - end if - end if - - if (pdims(3) > 1) then - if (periodic(3)) then - zlbndry = "exchange" - zubndry = "exchange" - else - if (pcoords(3) > 0 ) then - zlbndry = "exchange" - end if - if (pcoords(3) < pdims(3)-1) then - zubndry = "exchange" - end if - end if - end if -#endif /* MPI */ - #ifdef PROFILE ! stop accounting time for module initialization/finalization ! @@ -326,6 +346,10 @@ module boundaries ! call update_corners() +! convert updated primitive variables to conservative ones in all ghost cells +! + call update_ghost_cells() + #ifdef PROFILE ! stop accounting time for variable boundary update ! @@ -943,7 +967,7 @@ module boundaries ! local variables ! - integer :: n, i, j, k + integer :: i, j, k, p ! local pointers ! @@ -961,50 +985,50 @@ module boundaries ! iterate over all variables ! - do n = 1, nv + do p = 1, nv ! edges ! #if NDIMS == 3 do i = 1, im - pdata%u(n,i, 1:nh, 1:nh) = pdata%u(n,i,jbl,kbl) - pdata%u(n,i,jt:jm, 1:nh) = pdata%u(n,i,jeu,kbl) - pdata%u(n,i, 1:nh,kt:km) = pdata%u(n,i,jbl,keu) - pdata%u(n,i,jt:jm,kt:km) = pdata%u(n,i,jeu,keu) + pdata%q(p,i, 1:nh, 1:nh) = pdata%q(p,i,jbl,kbl) + pdata%q(p,i,jt:jm, 1:nh) = pdata%q(p,i,jeu,kbl) + pdata%q(p,i, 1:nh,kt:km) = pdata%q(p,i,jbl,keu) + pdata%q(p,i,jt:jm,kt:km) = pdata%q(p,i,jeu,keu) end do do j = 1, jm - pdata%u(n, 1:nh,j, 1:nh) = pdata%u(n,ibl,j,kbl) - pdata%u(n,it:im,j, 1:nh) = pdata%u(n,ieu,j,kbl) - pdata%u(n, 1:nh,j,kt:km) = pdata%u(n,ibl,j,keu) - pdata%u(n,it:im,j,kt:km) = pdata%u(n,ieu,j,keu) + pdata%q(p, 1:nh,j, 1:nh) = pdata%q(p,ibl,j,kbl) + pdata%q(p,it:im,j, 1:nh) = pdata%q(p,ieu,j,kbl) + pdata%q(p, 1:nh,j,kt:km) = pdata%q(p,ibl,j,keu) + pdata%q(p,it:im,j,kt:km) = pdata%q(p,ieu,j,keu) end do #endif /* == 3 */ do k = 1, km - pdata%u(n, 1:nh, 1:nh,k) = pdata%u(n,ibl,jbl,k) - pdata%u(n,it:im, 1:nh,k) = pdata%u(n,ieu,jbl,k) - pdata%u(n, 1:nh,jt:jm,k) = pdata%u(n,ibl,jeu,k) - pdata%u(n,it:im,jt:jm,k) = pdata%u(n,ieu,jeu,k) + pdata%q(p, 1:nh, 1:nh,k) = pdata%q(p,ibl,jbl,k) + pdata%q(p,it:im, 1:nh,k) = pdata%q(p,ieu,jbl,k) + pdata%q(p, 1:nh,jt:jm,k) = pdata%q(p,ibl,jeu,k) + pdata%q(p,it:im,jt:jm,k) = pdata%q(p,ieu,jeu,k) end do ! corners ! #if NDIMS == 3 - pdata%u(n, 1:nh, 1:nh, 1:nh) = pdata%u(n,ibl,jbl,kbl) - pdata%u(n,it:im, 1:nh, 1:nh) = pdata%u(n,ieu,jbl,kbl) - pdata%u(n, 1:nh,jt:jm, 1:nh) = pdata%u(n,ibl,jeu,kbl) - pdata%u(n,it:im,jt:jm, 1:nh) = pdata%u(n,ieu,jeu,kbl) - pdata%u(n, 1:nh, 1:nh,kt:km) = pdata%u(n,ibl,jbl,keu) - pdata%u(n,it:im, 1:nh,kt:km) = pdata%u(n,ieu,jbl,keu) - pdata%u(n, 1:nh,jt:jm,kt:km) = pdata%u(n,ibl,jeu,keu) - pdata%u(n,it:im,jt:jm,kt:km) = pdata%u(n,ieu,jeu,keu) + pdata%q(p, 1:nh, 1:nh, 1:nh) = pdata%q(p,ibl,jbl,kbl) + pdata%q(p,it:im, 1:nh, 1:nh) = pdata%q(p,ieu,jbl,kbl) + pdata%q(p, 1:nh,jt:jm, 1:nh) = pdata%q(p,ibl,jeu,kbl) + pdata%q(p,it:im,jt:jm, 1:nh) = pdata%q(p,ieu,jeu,kbl) + pdata%q(p, 1:nh, 1:nh,kt:km) = pdata%q(p,ibl,jbl,keu) + pdata%q(p,it:im, 1:nh,kt:km) = pdata%q(p,ieu,jbl,keu) + pdata%q(p, 1:nh,jt:jm,kt:km) = pdata%q(p,ibl,jeu,keu) + pdata%q(p,it:im,jt:jm,kt:km) = pdata%q(p,ieu,jeu,keu) #endif /* == 3 */ end do @@ -1021,6 +1045,110 @@ module boundaries ! !=============================================================================== ! +! subroutine UPDATE_GHOST_CELLS: +! ----------------------------- +! +! Subroutine updates conservative variables in all ghost cells from +! already updated primitive variables. +! +! +!=============================================================================== +! + subroutine update_ghost_cells() + +! include external variables +! + use blocks , only : block_data, list_data + use coordinates , only : im , jm , km , in , jn , kn + use coordinates , only : ib , jb , kb , ie , je , ke + use coordinates , only : ibl, jbl, kbl, ieu, jeu, keu + use equations , only : nv + use equations , only : prim2cons + +! local variables are not implicit by default +! + implicit none + +! local variables +! + integer :: i, j, k + +! local pointers +! + type(block_data), pointer :: pdata +! +!------------------------------------------------------------------------------- +! +! assign the pointer to the first block on the list +! + pdata => list_data + +! scan all data blocks until the last is reached +! + do while(associated(pdata)) + +! update the X and Y boundary ghost cells +! + do k = 1, km + +! update lower layers of the Y boundary +! + do j = 1, jbl + call prim2cons(im, pdata%q(1:nv,1:im,j,k), pdata%u(1:nv,1:im,j,k)) + end do ! j = 1, jbl + +! update upper layers of the Y boundary +! + do j = jeu, jm + call prim2cons(im, pdata%q(1:nv,1:im,j,k), pdata%u(1:nv,1:im,j,k)) + end do ! j = jeu, jm + +! update remaining left layers of the X boundary +! + do i = 1, ibl + call prim2cons(jn, pdata%q(1:nv,i,jb:je,k), pdata%u(1:nv,i,jb:je,k)) + end do ! i = 1, ibl + +! update remaining right layers of the X boundary +! + do i = ieu, im + call prim2cons(jn, pdata%q(1:nv,i,jb:je,k), pdata%u(1:nv,i,jb:je,k)) + end do ! i = 1, ibl + + end do ! k = 1, km + +#if NDIMS == 3 +! update the Z boundary ghost cells +! + do j = jb, je + +! update the remaining front layers of the Z boundary +! + do k = 1, kbl + call prim2cons(in, pdata%q(1:nv,ib:ie,j,k), pdata%u(1:nv,ib:ie,j,k)) + end do ! k = 1, kbl + +! update the remaining back layers of the Z boundary +! + do k = keu, km + call prim2cons(in, pdata%q(1:nv,ib:ie,j,k), pdata%u(1:nv,ib:ie,j,k)) + end do ! k = keu, km + + end do ! j = jb, je +#endif /* NDIMS == 3 */ + +! assign the pointer to the next block on the list +! + pdata => pdata%next + + end do ! data blocks + +!------------------------------------------------------------------------------- +! + end subroutine update_ghost_cells +! +!=============================================================================== +! ! subroutine SPECIFIC_BOUNDARIES: ! ------------------------------ ! @@ -1282,27 +1410,27 @@ module boundaries case(1) if (iside == 1) then call boundary_copy(pdata & - , pneigh%data%u(:,iel:ie,:,:), idir, iside) + , pneigh%data%q(:,iel:ie,:,:), idir, iside) else call boundary_copy(pdata & - , pneigh%data%u(:,ib:ibu,:,:), idir, iside) + , pneigh%data%q(:,ib:ibu,:,:), idir, iside) end if case(2) if (iside == 1) then call boundary_copy(pdata & - , pneigh%data%u(:,:,jel:je,:), idir, iside) + , pneigh%data%q(:,:,jel:je,:), idir, iside) else call boundary_copy(pdata & - , pneigh%data%u(:,:,jb:jbu,:), idir, iside) + , pneigh%data%q(:,:,jb:jbu,:), idir, iside) end if #if NDIMS == 3 case(3) if (iside == 1) then call boundary_copy(pdata & - , pneigh%data%u(:,:,:,kel:ke), idir, iside) + , pneigh%data%q(:,:,:,kel:ke), idir, iside) else call boundary_copy(pdata & - , pneigh%data%u(:,:,:,kb:kbu), idir, iside) + , pneigh%data%q(:,:,:,kb:kbu), idir, iside) end if #endif /* NDIMS == 3 */ end select @@ -1431,9 +1559,9 @@ module boundaries ! fill the buffer with data from the current block (depending on the side) ! if (pinfo%side == 1) then - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,iel:ie,:,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,iel:ie,:,:) else - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,ib:ibu,:,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,ib:ibu,:,:) end if ! associate the pointer with the next block @@ -1459,9 +1587,9 @@ module boundaries ! fill the buffer with data from the current block (depending on the side) ! if (pinfo%side == 1) then - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,jel:je,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,jel:je,:) else - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,jb:jbu,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,jb:jbu,:) end if ! associate the pointer with the next block @@ -1488,9 +1616,9 @@ module boundaries ! fill the buffer with data from the current block (depending on the side) ! if (pinfo%side == 1) then - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,:,kel:ke) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,:,kel:ke) else - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,:,kb:kbu) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,:,kb:kbu) end if ! associate the pointer with the next block @@ -1806,7 +1934,7 @@ module boundaries ! update boundaries of the current block ! call boundary_restrict(pdata & - , pneigh%data%u(:,il:iu,jl:ju,kl:ku) & + , pneigh%data%q(:,il:iu,jl:ju,kl:ku) & , idir, iside, iface) #ifdef MPI @@ -1938,7 +2066,7 @@ module boundaries ! fill the data buffer with the current block variable slices ! - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,il:iu,:,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,il:iu,:,:) ! associate the pointer with the next block ! @@ -1973,7 +2101,7 @@ module boundaries ! fill the data buffer with the current block variable slices ! - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,jl:ju,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,jl:ju,:) ! associate the pointer with the next block ! @@ -2009,7 +2137,7 @@ module boundaries ! fill the data buffer with the current block variable slices ! - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,:,kl:ku) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,:,kl:ku) ! associate the pointer with the next block ! @@ -2315,7 +2443,7 @@ module boundaries ! update boundaries of the current block from its neighbor ! call boundary_prolong(pdata & - , pneigh%data%u(:,il:iu,jl:ju,kl:ku) & + , pneigh%data%q(:,il:iu,jl:ju,kl:ku) & , idir, iside, nface) #ifdef MPI @@ -2449,7 +2577,7 @@ module boundaries ! fill the data buffer with the current block variable slices ! - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,il:iu,:,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,il:iu,:,:) ! associate the pointer with the next block ! @@ -2484,7 +2612,7 @@ module boundaries ! fill the data buffer with the current block variable slices ! - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,jl:ju,:) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,jl:ju,:) ! associate the pointer with the next block ! @@ -2520,7 +2648,7 @@ module boundaries ! fill the data buffer with the current block variable slices ! - rbuf(l,:,:,:,:) = pinfo%neigh%data%u(:,:,:,kl:ku) + rbuf(l,:,:,:,:) = pinfo%neigh%data%q(:,:,:,kl:ku) ! associate the pointer with the next block ! @@ -2664,11 +2792,12 @@ module boundaries ! import external procedures and variables ! use blocks , only : block_data - use coordinates , only : ng, im, jm, km, ib, ibl, ie, ieu, jb & - , jbl, je, jeu, kb, kbl, ke, keu + use coordinates , only : im , jm , km , ng + use coordinates , only : ib , jb , kb , ie , je , ke + use coordinates , only : ibl, jbl, kbl, ieu, jeu, keu use equations , only : nv - use equations , only : idn, imx, imy, imz, ibx, iby, ibz, ibp - use error , only : print_warning + use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp + use error , only : print_error, print_warning ! local variables are not implicit by default ! @@ -2699,25 +2828,36 @@ module boundaries ! apply selected boundary condition ! - select case(xlbndry) + select case(bnd_type(idir,iside)) - case("reflecting", "reflect") +! "open" boundary conditions +! + case(bnd_open) + + do i = 1, ng + pdata%q( :,i,:,:) = pdata%q(:,ib,:,:) + end do + +! "reflective" boundary conditions +! + case(bnd_reflective) do i = 1, ng it = ib - i is = ibl + i - pdata%u( :,it,:,:) = pdata%u( :,is,:,:) - pdata%u(imx,it,:,:) = - pdata%u(imx,is,:,:) + pdata%q( :,it,:,:) = pdata%q( :,is,:,:) + pdata%q(ivx,it,:,:) = - pdata%q(ivx,is,:,:) end do - case default ! "open" as default boundary conditions +! wrong boundary conditions +! + case default - do i = 1, ng - pdata%u( :,i,:,:) = pdata%u(:,ib,:,:) - end do + call print_error("boundaries:boundary_specific()" & + , "Wrong left X boundary type!") end select @@ -2727,23 +2867,34 @@ module boundaries ! apply selected boundary condition ! - select case(xubndry) + select case(bnd_type(idir,iside)) - case("reflecting", "reflect") +! "open" boundary conditions +! + case(bnd_open) + + do i = ieu, im + pdata%q( :,i ,:,:) = pdata%q( :,ie,:,:) + end do + +! "reflective" boundary conditions +! + case(bnd_reflective) do i = 1, ng it = ie + i is = ieu - i - pdata%u( :,it,:,:) = pdata%u( :,is,:,:) - pdata%u(imx,it,:,:) = - pdata%u(imx,is,:,:) + pdata%q( :,it,:,:) = pdata%q( :,is,:,:) + pdata%q(ivx,it,:,:) = - pdata%q(ivx,is,:,:) end do - case default ! "open" as default boundary conditions +! wrong boundary conditions +! + case default - do i = ieu, im - pdata%u( :,i,:,:) = pdata%u(:,ie,:,:) - end do + call print_error("boundaries:boundary_specific()" & + , "Wrong right X boundary type!") end select @@ -2753,23 +2904,34 @@ module boundaries ! apply selected boundary condition ! - select case(ylbndry) + select case(bnd_type(idir,iside)) - case("reflecting", "reflect") +! "open" boundary conditions +! + case(bnd_open) + + do j = 1, ng + pdata%q( :,:,j ,:) = pdata%q( :,:,jb,:) + end do + +! "reflective" boundary conditions +! + case(bnd_reflective) do j = 1, ng jt = jb - j js = jbl + j - pdata%u( :,:,jt,:) = pdata%u( :,:,js,:) - pdata%u(imy,:,jt,:) = - pdata%u(imy,:,js,:) + pdata%q( :,:,jt,:) = pdata%q( :,:,js,:) + pdata%q(ivy,:,jt,:) = - pdata%q(ivy,:,js,:) end do - case default ! "open" as default boundary conditions +! wrong boundary conditions +! + case default - do j = 1, ng - pdata%u( :,:,j,:) = pdata%u(:,:,jb,:) - end do + call print_error("boundaries:boundary_specific()" & + , "Wrong left Y boundary type!") end select @@ -2779,23 +2941,34 @@ module boundaries ! apply selected boundary condition ! - select case(yubndry) + select case(bnd_type(idir,iside)) - case("reflecting", "reflect") +! "open" boundary conditions +! + case(bnd_open) + + do j = jeu, jm + pdata%q( :,:,j ,:) = pdata%q( :,:,je,:) + end do + +! "reflective" boundary conditions +! + case(bnd_reflective) do j = 1, ng jt = je + j js = jeu - j - pdata%u( :,:,jt,:) = pdata%u( :,:,js,:) - pdata%u(imy,:,jt,:) = - pdata%u(imy,:,js,:) + pdata%q( :,:,jt,:) = pdata%q( :,:,js,:) + pdata%q(ivy,:,jt,:) = - pdata%q(ivy,:,js,:) end do - case default ! "open" as default boundary conditions +! wrong boundary conditions +! + case default - do j = jeu, jm - pdata%u( :,:,j,:) = pdata%u(:,:,je,:) - end do + call print_error("boundaries:boundary_specific()" & + , "Wrong right Y boundary type!") end select @@ -2806,23 +2979,34 @@ module boundaries ! apply selected boundary condition ! - select case(zlbndry) + select case(bnd_type(idir,iside)) - case("reflecting", "reflect") +! "open" boundary conditions +! + case(bnd_open) + + do k = 1, ng + pdata%q( :,:,:,k ) = pdata%q( :,:,:,kb) + end do + +! "reflective" boundary conditions +! + case(bnd_reflective) do k = 1, ng kt = kb - k ks = kbl + k - pdata%u( :,:,:,kt) = pdata%u( :,:,:,ks) - pdata%u(imz,:,:,kt) = - pdata%u(imz,:,:,ks) + pdata%q( :,:,:,kt) = pdata%q( :,:,:,ks) + pdata%q(ivz,:,:,kt) = - pdata%q(ivz,:,:,ks) end do - case default ! "open" as default boundary conditions +! wrong boundary conditions +! + case default - do k = 1, ng - pdata%u( :,:,:,k) = pdata%u(:,:,:,kb) - end do + call print_error("boundaries:boundary_specific()" & + , "Wrong left Z boundary type!") end select @@ -2832,23 +3016,34 @@ module boundaries ! apply selected boundary condition ! - select case(zubndry) + select case(bnd_type(idir,iside)) - case("reflecting", "reflect") +! "open" boundary conditions +! + case(bnd_open) + + do k = keu, km + pdata%q( :,:,:,k ) = pdata%q( :,:,:,ke) + end do + +! "reflective" boundary conditions +! + case(bnd_reflective) do k = 1, ng kt = ke + k ks = keu - k - pdata%u( :,:,:,kt) = pdata%u( :,:,:,ks) - pdata%u(imz,:,:,kt) = - pdata%u(imz,:,:,ks) + pdata%q( :,:,:,kt) = pdata%q( :,:,:,ks) + pdata%q(ivz,:,:,kt) = - pdata%q(ivz,:,:,ks) end do - case default ! "open" as default boundary conditions +! wrong boundary conditions +! + case default - do k = keu, km - pdata%u( :,:,:,k) = pdata%u(:,:,:,ke) - end do + call print_error("boundaries:boundary_specific()" & + , "Wrong right Z boundary type!") end select #endif /* NDIMS == 3 */ @@ -2876,13 +3071,13 @@ module boundaries ! Arguments: ! ! pdata - the pointer to modified data block; -! u - the variable array from which boundaries are updated; +! q - the variable array from which boundaries are updated; ! idir - the direction to be processed; ! iside - the side to be processed; ! !=============================================================================== ! - subroutine boundary_copy(pdata, u, idir, iside) + subroutine boundary_copy(pdata, q, idir, iside) ! import external procedures and variables ! @@ -2897,7 +3092,7 @@ module boundaries ! subroutine arguments ! type(block_data), pointer , intent(inout) :: pdata - real , dimension(:,:,:,:), intent(in) :: u + real , dimension(:,:,:,:), intent(in) :: q integer , intent(in) :: idir, iside ! !------------------------------------------------------------------------------- @@ -2909,26 +3104,26 @@ module boundaries case(1) if (iside == 1) then - pdata%u(1:nv, 1:ibl,1:jm,1:km) = u(1:nv,1:ng,1:jm,1:km) + pdata%q(1:nv, 1:ibl,1:jm,1:km) = q(1:nv,1:ng,1:jm,1:km) else - pdata%u(1:nv,ieu:im ,1:jm,1:km) = u(1:nv,1:ng,1:jm,1:km) + pdata%q(1:nv,ieu:im ,1:jm,1:km) = q(1:nv,1:ng,1:jm,1:km) end if case(2) if (iside == 1) then - pdata%u(1:nv,1:im, 1:jbl,1:km) = u(1:nv,1:im,1:ng,1:km) + pdata%q(1:nv,1:im, 1:jbl,1:km) = q(1:nv,1:im,1:ng,1:km) else - pdata%u(1:nv,1:im,jeu:jm ,1:km) = u(1:nv,1:im,1:ng,1:km) + pdata%q(1:nv,1:im,jeu:jm ,1:km) = q(1:nv,1:im,1:ng,1:km) end if #if NDIMS == 3 case(3) if (iside == 1) then - pdata%u(1:nv,1:im,1:jm, 1:kbl) = u(1:nv,1:im,1:jm,1:ng) + pdata%q(1:nv,1:im,1:jm, 1:kbl) = q(1:nv,1:im,1:jm,1:ng) else - pdata%u(1:nv,1:im,1:jm,keu:km ) = u(1:nv,1:im,1:jm,1:ng) + pdata%q(1:nv,1:im,1:jm,keu:km ) = q(1:nv,1:im,1:jm,1:ng) end if #endif /* NDIMS == 3 */ @@ -2950,12 +3145,12 @@ module boundaries ! Arguments: ! ! pdata - the input data block; -! u - the conserved array; +! q - the variable array from which boundaries are updated; ! idir, iside, iface - the positions of the neighbor block; ! !=============================================================================== ! - subroutine boundary_restrict(pdata, u, idir, iside, iface) + subroutine boundary_restrict(pdata, q, idir, iside, iface) ! import external procedures and variables ! @@ -2972,7 +3167,7 @@ module boundaries ! subroutine arguments ! type(block_data) , pointer, intent(inout) :: pdata - real(kind=8) , dimension(:,:,:,:), intent(in) :: u + real(kind=8) , dimension(:,:,:,:), intent(in) :: q integer , intent(in) :: idir, iside, iface ! local variables @@ -3112,22 +3307,22 @@ module boundaries ! update boundaries of the conserved variables ! #if NDIMS == 2 - pdata%u(:,is:it,js:jt, 1 ) = & - 2.50d-01 * ((u(1:nv,il:iu:2,jl:ju:2, 1 ) & - + u(1:nv,ip:iu:2,jp:ju:2, 1 )) & - + (u(1:nv,il:iu:2,jp:ju:2, 1 ) & - + u(1:nv,ip:iu:2,jl:ju:2, 1 ))) + pdata%q(:,is:it,js:jt, 1 ) = & + 2.50d-01 * ((q(1:nv,il:iu:2,jl:ju:2, 1 ) & + + q(1:nv,ip:iu:2,jp:ju:2, 1 )) & + + (q(1:nv,il:iu:2,jp:ju:2, 1 ) & + + q(1:nv,ip:iu:2,jl:ju:2, 1 ))) #endif /* NDIMS == 2 */ #if NDIMS == 3 - pdata%u(:,is:it,js:jt,ks:kt) = & - 1.25d-01 * (((u(1:nv,il:iu:2,jl:ju:2,kl:ku:2) & - + u(1:nv,ip:iu:2,jp:ju:2,kp:ku:2)) & - + (u(1:nv,il:iu:2,jl:ju:2,kp:ku:2) & - + u(1:nv,ip:iu:2,jp:ju:2,kl:ku:2))) & - + ((u(1:nv,il:iu:2,jp:ju:2,kp:ku:2) & - + u(1:nv,ip:iu:2,jl:ju:2,kl:ku:2)) & - + (u(1:nv,il:iu:2,jp:ju:2,kl:ku:2) & - + u(1:nv,ip:iu:2,jl:ju:2,kp:ku:2)))) + pdata%q(:,is:it,js:jt,ks:kt) = & + 1.25d-01 * (((q(1:nv,il:iu:2,jl:ju:2,kl:ku:2) & + + q(1:nv,ip:iu:2,jp:ju:2,kp:ku:2)) & + + (q(1:nv,il:iu:2,jl:ju:2,kp:ku:2) & + + q(1:nv,ip:iu:2,jp:ju:2,kl:ku:2))) & + + ((q(1:nv,il:iu:2,jp:ju:2,kp:ku:2) & + + q(1:nv,ip:iu:2,jl:ju:2,kl:ku:2)) & + + (q(1:nv,il:iu:2,jp:ju:2,kl:ku:2) & + + q(1:nv,ip:iu:2,jl:ju:2,kp:ku:2)))) #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- @@ -3146,12 +3341,12 @@ module boundaries ! Arguments: ! ! pdata - the input data block; -! u - the conserved array; +! q - the variable array from which boundaries are updated; ! idir, iside, iface - the positions of the neighbor block; ! !=============================================================================== ! - subroutine boundary_prolong(pdata, u, idir, iside, iface) + subroutine boundary_prolong(pdata, q, idir, iside, iface) ! import external procedures and variables ! @@ -3169,16 +3364,16 @@ module boundaries ! subroutine arguments ! type(block_data), pointer , intent(inout) :: pdata - real , dimension(:,:,:,:), intent(in) :: u + real , dimension(:,:,:,:), intent(in) :: q integer , intent(in) :: idir, iside, iface ! local variables ! - integer :: i, j, k, q + integer :: i, j, k, p integer :: ic, jc, kc, ip, jp, kp integer :: il, jl, kl, iu, ju, ku integer :: is, js, ks, it, jt, kt - real :: dul, dur, dux, duy, duz + real :: dql, dqr, dqx, dqy, dqz, dq1, dq2, dq3, dq4 ! !------------------------------------------------------------------------------- ! @@ -3295,37 +3490,43 @@ module boundaries ! iterate over all variables ! - do q = 1, nv + do p = 1, nv - dul = u(q,i ,j,k) - u(q,i-1,j,k) - dur = u(q,i+1,j,k) - u(q,i ,j,k) - dux = limiter(0.25d+00, dul, dur) + dql = q(p,i ,j,k) - q(p,i-1,j,k) + dqr = q(p,i+1,j,k) - q(p,i ,j,k) + dqx = limiter(0.25d+00, dql, dqr) - dul = u(q,i,j ,k) - u(q,i,j-1,k) - dur = u(q,i,j+1,k) - u(q,i,j ,k) - duy = limiter(0.25d+00, dul, dur) + dql = q(p,i,j ,k) - q(p,i,j-1,k) + dqr = q(p,i,j+1,k) - q(p,i,j ,k) + dqy = limiter(0.25d+00, dql, dqr) #if NDIMS == 3 - dul = u(q,i,j,k ) - u(q,i,j,k-1) - dur = u(q,i,j,k+1) - u(q,i,j,k ) - duz = limiter(0.25d+00, dul, dur) + dql = q(p,i,j,k ) - q(p,i,j,k-1) + dqr = q(p,i,j,k+1) - q(p,i,j,k ) + dqz = limiter(0.25d+00, dql, dqr) #endif /* NDIMS == 3 */ #if NDIMS == 2 - pdata%u(q,it,jt,kt) = u(q,i,j,k) - (dux + duy) - pdata%u(q,ip,jt,kt) = u(q,i,j,k) + (dux - duy) - pdata%u(q,it,jp,kt) = u(q,i,j,k) + (duy - dux) - pdata%u(q,ip,jp,kt) = u(q,i,j,k) + (dux + duy) + dq1 = dqx + dqy + dq2 = dqx - dqy + pdata%q(p,it,jt,kt) = q(p,i,j,k) - dq1 + pdata%q(p,ip,jt,kt) = q(p,i,j,k) + dq2 + pdata%q(p,it,jp,kt) = q(p,i,j,k) - dq2 + pdata%q(p,ip,jp,kt) = q(p,i,j,k) + dq1 #endif /* NDIMS == 2 */ #if NDIMS == 3 - pdata%u(q,it,jt,kt) = u(q,i,j,k) - (dux + duy + duz) - pdata%u(q,ip,jt,kt) = u(q,i,j,k) + (dux - duy - duz) - pdata%u(q,it,jp,kt) = u(q,i,j,k) - (dux - duy + duz) - pdata%u(q,ip,jp,kt) = u(q,i,j,k) + (dux + duy - duz) - pdata%u(q,it,jt,kp) = u(q,i,j,k) - (dux + duy - duz) - pdata%u(q,ip,jt,kp) = u(q,i,j,k) + (dux - duy + duz) - pdata%u(q,it,jp,kp) = u(q,i,j,k) - (dux - duy - duz) - pdata%u(q,ip,jp,kp) = u(q,i,j,k) + (dux + duy + duz) + dq1 = dqx + dqy + dqz + dq2 = dqx - dqy - dqz + dq3 = dqx - dqy + dqz + dq4 = dqx + dqy - dqz + pdata%q(p,it,jt,kt) = q(p,i,j,k) - dq1 + pdata%q(p,ip,jt,kt) = q(p,i,j,k) + dq2 + pdata%q(p,it,jp,kt) = q(p,i,j,k) - dq3 + pdata%q(p,ip,jp,kt) = q(p,i,j,k) + dq4 + pdata%q(p,it,jt,kp) = q(p,i,j,k) - dq4 + pdata%q(p,ip,jt,kp) = q(p,i,j,k) + dq3 + pdata%q(p,it,jp,kp) = q(p,i,j,k) - dq2 + pdata%q(p,ip,jp,kp) = q(p,i,j,k) + dq1 #endif /* NDIMS == 3 */ end do ! q - 1, nv @@ -3427,8 +3628,8 @@ module boundaries k1 = 2 * (k - kl) + kb k2 = k1 + 1 - pdata%f(idir,:,it,j,k) = 2.5d-01 * (f(:,j1,k1) + f(:,j2,k1) & - + f(:,j1,k2) + f(:,j2,k2)) + pdata%f(idir,:,it,j,k) = 2.5d-01 * ((f(:,j1,k1) + f(:,j2,k2)) & + + (f(:,j2,k1) + f(:,j1,k2))) end do #endif /* NDIMS == 3 */ end do @@ -3475,8 +3676,8 @@ module boundaries k1 = 2 * (k - kl) + kb k2 = k1 + 1 - pdata%f(idir,:,i,jt,k) = 2.5d-01 * (f(:,i1,k1) + f(:,i2,k1) & - + f(:,i1,k2) + f(:,i2,k2)) + pdata%f(idir,:,i,jt,k) = 2.5d-01 * ((f(:,i1,k1) + f(:,i2,k2)) & + + (f(:,i2,k1) + f(:,i1,k2))) end do #endif /* NDIMS == 3 */ end do @@ -3516,8 +3717,8 @@ module boundaries j1 = 2 * (j - jl) + jb j2 = j1 + 1 - pdata%f(idir,:,i,j,kt) = 2.5d-01 * (f(:,i1,j1) + f(:,i2,j1) & - + f(:,i1,j2) + f(:,i2,j2)) + pdata%f(idir,:,i,j,kt) = 2.5d-01 * ((f(:,i1,j1) + f(:,i2,j2)) & + + (f(:,i2,j1) + f(:,i1,j2))) end do end do #endif /* NDIMS == 3 */ diff --git a/src/domains.F90 b/src/domains.F90 index f32de85..8c3eda3 100644 --- a/src/domains.F90 +++ b/src/domains.F90 @@ -143,8 +143,7 @@ module domains use blocks , only : metablock_set_configuration use blocks , only : metablock_set_coordinates, metablock_set_bounds use blocks , only : nsides, nfaces - use boundaries , only : xlbndry, ylbndry, zlbndry - use boundaries , only : xubndry, yubndry, zubndry + use boundaries , only : bnd_type, bnd_periodic use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax use coordinates , only : ir, jr, kr @@ -349,7 +348,7 @@ module domains ! if periodic boundary conditions set edge block neighbors ! - if (xlbndry == 'periodic' .and. xubndry == 'periodic') then + if (bnd_type(1,1) == bnd_periodic .and. bnd_type(1,2) == bnd_periodic) then do k = 1, kr do j = 1, jr @@ -395,7 +394,7 @@ module domains ! if periodic boundary conditions set edge block neighbors ! - if (ylbndry == 'periodic' .and. yubndry == 'periodic') then + if (bnd_type(2,1) == bnd_periodic .and. bnd_type(2,2) == bnd_periodic) then do k = 1, kr do i = 1, ir @@ -442,7 +441,7 @@ module domains ! if periodic boundary conditions set edge block neighbors ! - if (zlbndry == 'periodic' .and. zubndry == 'periodic') then + if (bnd_type(3,1) == bnd_periodic .and. bnd_type(3,2) == bnd_periodic) then do j = 1, jr do i = 1, ir diff --git a/src/driver.F90 b/src/driver.F90 index a37a866..4f17446 100644 --- a/src/driver.F90 +++ b/src/driver.F90 @@ -81,7 +81,7 @@ program amun ! integer, dimension(3) :: div = 1 logical, dimension(3) :: per = .true. - integer :: nmax = 0, ndat = 1 + integer :: nmax = huge(1), ndat = 1 real :: tmax = 0.0d+00, trun = 9.999d+03, tsav = 3.0d+01 real(kind=8) :: dtnext = 0.0d+00 diff --git a/src/equations.F90 b/src/equations.F90 index bde35e8..71cd7e8 100644 --- a/src/equations.F90 +++ b/src/equations.F90 @@ -605,7 +605,8 @@ module equations ! include external procedures and variables ! - use coordinates, only : im, jm, km + use coordinates, only : im, jm, km, in, jn, kn + use coordinates, only : ib, jb, kb, ie, je, ke ! local variables are not implicit by default ! @@ -618,29 +619,21 @@ module equations ! temporary variables ! - integer :: j, k - -! temporary array to store conserved variable vector -! - real(kind=8), dimension(nv,im) :: u + integer :: j, k ! !------------------------------------------------------------------------------- ! ! update primitive variables ! - do k = 1, km - do j = 1, jm - -! copy variables to temporary array of conserved variables -! - u(1:nv,1:im) = uu(1:nv,1:im,j,k) + do k = kb, ke + do j = jb, je ! convert conserved variables to primitive ones ! - call cons2prim(im, u(1:nv,1:im), qq(1:nv,1:im,j,k)) + call cons2prim(in, uu(1:nv,ib:ie,j,k), qq(1:nv,ib:ie,j,k)) - end do ! j = 1, jm - end do ! k = 1, km + end do ! j = jb, je + end do ! k = kb, ke !------------------------------------------------------------------------------- ! diff --git a/src/evolution.F90 b/src/evolution.F90 index 9091111..eadb831 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -249,14 +249,14 @@ module evolution ! call update_mesh() -! update boundaries -! - call boundary_variables() - ! update primitive variables ! call update_variables() +! update boundaries +! + call boundary_variables() + ! set all meta blocks to be updated ! call set_blocks_update(.true.) @@ -467,14 +467,14 @@ module evolution end do -! update boundaries -! - call boundary_variables() - ! update primitive variables ! call update_variables() +! update boundaries +! + call boundary_variables() + !------------------------------------------------------------------------------- ! end subroutine evolve_euler @@ -552,14 +552,14 @@ module evolution end do -! update boundaries -! - call boundary_variables() - ! update primitive variables ! call update_variables() +! update boundaries +! + call boundary_variables() + ! update fluxes for the second step of the RK2 integration ! call update_fluxes() @@ -598,14 +598,14 @@ module evolution end do -! update boundaries -! - call boundary_variables() - ! update primitive variables ! call update_variables() +! update boundaries +! + call boundary_variables() + !------------------------------------------------------------------------------- ! end subroutine evolve_rk2 diff --git a/src/io.F90 b/src/io.F90 index cecf5a3..d9948c3 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -57,17 +57,22 @@ module io ! nrest - for job restarting, this is the number of restart snapshot; ! irest - the local counter for the restart snapshots; ! isnap - the local counter for the regular snapshots; +! ishift - the shift of the snapshot counter for restarting job with +! different snapshot interval; ! hrest - the execution time interval for restart snapshot storing -! (in hours); +! (in hours); the minimum allowed value is 3 minutes; ! hsnap - the problem time interval for regular snapshot storing; +! tsnap - the next snapshot time; ! character(len=255), save :: respath = "./" character , save :: ftype = "p" integer , save :: nrest = -1 integer(kind=4) , save :: irest = 1 - integer(kind=4) , save :: isnap = -1 - real , save :: hrest = 6.0d+00 - real , save :: hsnap = 1.0d+00 + integer(kind=4) , save :: isnap = 0 + integer(kind=4) , save :: ishift = 0 + real , save :: hrest = 6.0e+00 + real , save :: hsnap = 1.0e+00 + real , save :: tsnap = 0.0e+00 ! flags to determine the way of data writing ! @@ -283,9 +288,10 @@ module io call read_restart_snapshot_h5() #endif /* HDF5 */ -! calculate the next snapshot number +! calculate the shift of the snapshot counter, and the next snapshot time ! - isnap = int(time / hsnap) + ishift = int(time / hsnap) - isnap + 1 + tsnap = (ishift + isnap) * hsnap #ifdef PROFILE ! stop accounting time for the data reading @@ -302,14 +308,14 @@ module io ! subroutine WRITE_RESTART_SNAPSHOTS: ! ---------------------------------- ! -! Subroutine stores current snapshot in restart files. -! This is a wrapper calling specific format subroutine. +! Subroutine stores current restart snapshot files. This is a wrapper +! calling specific format subroutine. ! ! Arguments: ! -! thrs - the interval time in hours for storing the restart snapshots; -! nrun - the restart snapshot number; -! iret - the return flag to inform if subroutine succeeded or failed; +! thrs - the current execution time in hours; +! nrun - the run number; +! iret - the return flag; ! !=============================================================================== ! @@ -327,28 +333,25 @@ module io ! !------------------------------------------------------------------------------- ! +! check if conditions for storing the restart snapshot have been met +! + if (hrest < 5.0e-02 .or. thrs < irest * hrest) return + #ifdef PROFILE ! start accounting time for the data writing ! call start_timer(iow) #endif /* PROFILE */ -! store restart snapshots at constant interval or when the job is done only -! if hrest is a positive number -! - if (hrest > 0.0d+00 .and. thrs >= irest * hrest) then - #ifdef HDF5 ! store restart file ! - call write_restart_snapshot_h5(nrun, iret) + call write_restart_snapshot_h5(nrun, iret) #endif /* HDF5 */ ! increase the restart snapshot counter ! - irest = irest + 1 - - end if + irest = irest + 1 #ifdef PROFILE ! stop accounting time for the data writing @@ -366,7 +369,7 @@ module io ! ------------------------- ! ! Subroutine stores block data in snapshots. Block variables are grouped -! todether and stored in big 4D arrays separately. This is a wrapper for +! together and stored in big 4D arrays separately. This is a wrapper for ! specific format storing. ! ! @@ -384,9 +387,9 @@ module io ! !------------------------------------------------------------------------------- ! -! exit the subroutine, if the time of the next snapshot is not reached +! check if conditions for storing the regular snapshot have been met ! - if (hsnap <= 0.0d+00 .or. isnap >= (int(time / hsnap))) return + if (hsnap <= 0.0e+00 .or. time < tsnap) return #ifdef PROFILE ! start accounting time for the data writing @@ -394,16 +397,17 @@ module io call start_timer(iow) #endif /* PROFILE */ -! increase the file counter -! - isnap = isnap + 1 - #ifdef HDF5 ! store variable snapshot file ! call write_snapshot_h5() #endif /* HDF5 */ +! increase the snapshot counter and calculate the next snapshot time +! + isnap = isnap + 1 + tsnap = (ishift + isnap) * hsnap + #ifdef PROFILE ! stop accounting time for the data writing ! @@ -419,7 +423,7 @@ module io ! function NEXT_TOUT: ! ------------------ ! -! Function returns time of the next data snapshot. +! Function returns the next data snapshot time. ! ! !=============================================================================== @@ -433,7 +437,7 @@ module io !------------------------------------------------------------------------------- ! if (hsnap > 0.0d+00) then - next_tout = (isnap + 1) * hsnap + next_tout = tsnap else next_tout = huge(hsnap) end if @@ -465,6 +469,7 @@ module io ! import external procedures and variables ! + use blocks , only : change_blocks_process use error , only : print_error use hdf5 , only : hid_t use hdf5 , only : H5F_ACC_RDONLY_F @@ -540,7 +545,7 @@ module io return end if -! opent the HDF5 file +! open the HDF5 file ! call h5fopen_f(fl, H5F_ACC_RDONLY_F, fid, err) @@ -665,6 +670,10 @@ module io ! do lfile = nprocs, nfiles - 1 +! switch meta blocks from the read file to belong to the reading process +! + call change_blocks_process(lfile, nprocs - 1) + ! read the remaining files by the last process only ! if (nproc == nprocs - 1) then @@ -1026,46 +1035,48 @@ module io ! !=============================================================================== ! -! write_attributes_h5: subroutine writes attributes in the HDF5 format -! connected to the provided identificator +! subroutine WRITE_ATTRIBUTES_H5: +! ------------------------------ ! -! info: this subroutine stores only the global attributes +! Subroutine stores global attributes in the HDF5 file provided by an +! identifier. ! -! arguments: -! fid - the HDF5 file identificator +! Arguments: +! +! fid - the HDF5 file identifier; ! !=============================================================================== ! subroutine write_attributes_h5(fid) -! references to other modules +! import external procedures and variables ! - use blocks , only : get_mblocks, get_dblocks, get_nleafs - use blocks , only : get_last_id - use coordinates, only : nn, ng, in, jn, kn, minlev, maxlev, toplev, ir, jr, kr - use coordinates, only : xmin, xmax, ymin, ymax, zmin, zmax - use error , only : print_error - use evolution, only : step, time, dt, dtn - use hdf5 , only : hid_t - use hdf5 , only : h5gcreate_f, h5gclose_f - use mpitools , only : nprocs, nproc - use random , only : nseeds, get_seeds + use blocks , only : get_mblocks, get_dblocks, get_nleafs + use blocks , only : get_last_id + use coordinates , only : minlev, maxlev, toplev + use coordinates , only : nn, ng, in, jn, kn, ir, jr, kr + use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax + use error , only : print_error + use evolution , only : step, time, dt, dtn + use hdf5 , only : hid_t + use hdf5 , only : h5gcreate_f, h5gclose_f + use mpitools , only : nprocs, nproc + use random , only : nseeds, get_seeds -! declare variables +! local variables are not implicit by default ! implicit none -! input variables +! subroutine arguments ! integer(hid_t), intent(in) :: fid ! local variables ! - integer(hid_t) :: gid - integer(kind=4) :: dm(3) - integer :: err + integer(hid_t) :: gid + integer :: err -! allocatable arrays +! local allocatable arrays ! integer(kind=4), dimension(:), allocatable :: seeds ! @@ -1077,85 +1088,87 @@ module io ! check if the group has been created successfuly ! - if (err .ge. 0) then - -! store the integer attributes -! - call write_attribute_integer_h5(gid, 'ndims' , NDIMS) - call write_attribute_integer_h5(gid, 'last_id', get_last_id()) - call write_attribute_integer_h5(gid, 'mblocks', get_mblocks()) - call write_attribute_integer_h5(gid, 'dblocks', get_dblocks()) - call write_attribute_integer_h5(gid, 'nleafs' , get_nleafs()) - call write_attribute_integer_h5(gid, 'ncells' , nn) - call write_attribute_integer_h5(gid, 'nghost' , ng) - call write_attribute_integer_h5(gid, 'minlev' , minlev) - call write_attribute_integer_h5(gid, 'maxlev' , maxlev) - call write_attribute_integer_h5(gid, 'toplev' , toplev) - call write_attribute_integer_h5(gid, 'nprocs' , nprocs) - call write_attribute_integer_h5(gid, 'nproc' , nproc) - call write_attribute_integer_h5(gid, 'nseeds' , nseeds) - call write_attribute_integer_h5(gid, 'step' , step ) - -! store the real attributes -! - call write_attribute_double_h5(gid, 'xmin', xmin) - call write_attribute_double_h5(gid, 'xmax', xmax) - call write_attribute_double_h5(gid, 'ymin', ymin) - call write_attribute_double_h5(gid, 'ymax', ymax) - call write_attribute_double_h5(gid, 'zmin', zmin) - call write_attribute_double_h5(gid, 'zmax', zmax) - call write_attribute_double_h5(gid, 'time', time) - call write_attribute_double_h5(gid, 'dt' , dt ) - call write_attribute_double_h5(gid, 'dtn' , dtn ) - -! store the vector attributes -! - dm(:) = (/ in, jn, kn /) - call write_attribute_vector_integer_h5(gid, 'dims' , 3, dm(:)) - call write_attribute_vector_integer_h5(gid, 'rdims', 3, (/ ir, jr, kr /)) - -! store random number generator seed values -! - if (nseeds .gt. 0) then - -! allocate space for seeds -! - allocate(seeds(nseeds)) - -! get the seed values -! - call get_seeds(seeds) - -! store them in the current group -! - call write_attribute_vector_integer_h5(gid, 'seeds', nseeds, seeds(:)) - -! deallocate seed array -! - deallocate(seeds) - - end if ! nseeds > 0 - -! close the group -! - call h5gclose_f(gid, err) - -! check if the group has been closed successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the group -! - call print_error("io::write_attributes_h5", "Cannot close the group!") - - end if - - else + if (err < 0) then ! print error about the problem with creating the group ! call print_error("io::write_attributes_h5", "Cannot create the group!") +! return from the subroutine +! + return + + end if + +! store the integer attributes +! + call write_attribute_integer_h5(gid, 'ndims' , NDIMS ) + call write_attribute_integer_h5(gid, 'last_id', get_last_id()) + call write_attribute_integer_h5(gid, 'mblocks', get_mblocks()) + call write_attribute_integer_h5(gid, 'dblocks', get_dblocks()) + call write_attribute_integer_h5(gid, 'nleafs' , get_nleafs() ) + call write_attribute_integer_h5(gid, 'ncells' , nn ) + call write_attribute_integer_h5(gid, 'nghost' , ng ) + call write_attribute_integer_h5(gid, 'minlev' , minlev ) + call write_attribute_integer_h5(gid, 'maxlev' , maxlev ) + call write_attribute_integer_h5(gid, 'toplev' , toplev ) + call write_attribute_integer_h5(gid, 'nprocs' , nprocs ) + call write_attribute_integer_h5(gid, 'nproc' , nproc ) + call write_attribute_integer_h5(gid, 'nseeds' , nseeds ) + call write_attribute_integer_h5(gid, 'step' , step ) + call write_attribute_integer_h5(gid, 'isnap' , isnap ) + +! store the real attributes +! + call write_attribute_double_h5(gid, 'xmin', xmin) + call write_attribute_double_h5(gid, 'xmax', xmax) + call write_attribute_double_h5(gid, 'ymin', ymin) + call write_attribute_double_h5(gid, 'ymax', ymax) + call write_attribute_double_h5(gid, 'zmin', zmin) + call write_attribute_double_h5(gid, 'zmax', zmax) + call write_attribute_double_h5(gid, 'time', time) + call write_attribute_double_h5(gid, 'dt' , dt ) + call write_attribute_double_h5(gid, 'dtn' , dtn ) + +! store the vector attributes +! + call write_attribute_vector_integer_h5(gid, 'dims' , (/ in, jn, kn /)) + call write_attribute_vector_integer_h5(gid, 'rdims', (/ ir, jr, kr /)) + +! store random number generator seed values +! + if (nseeds > 0) then + +! allocate space for seeds +! + allocate(seeds(nseeds)) + +! get the seed values +! + call get_seeds(seeds) + +! store them in the current group +! + call write_attribute_vector_integer_h5(gid, 'seeds', seeds(:)) + +! deallocate seed array +! + deallocate(seeds) + + end if ! nseeds > 0 + +! close the group +! + call h5gclose_f(gid, err) + +! check if the group has been closed successfuly +! + if (err < 0) then + +! print error about the problem with closing the group +! + call print_error("io::write_attributes_h5", "Cannot close the group!") + end if !------------------------------------------------------------------------------- @@ -1164,57 +1177,55 @@ module io ! !=============================================================================== ! -! read_attributes_h5: subroutine restores attributes from an HDF5 file linked -! to the HDF5 file identificator +! subroutine READ_ATTRIBUTES_H5: +! ----------------------------- ! -! info: this subroutine restores only the global attributes +! Subroutine restores global attributes from an HDF5 file provided by its +! identifier. ! -! arguments: -! fid - the HDF5 file identificator +! Arguments: +! +! fid - the HDF5 file identifier; ! !=============================================================================== ! subroutine read_attributes_h5(fid) -! references to other modules +! import external procedures and variables ! - use blocks , only : block_meta, block_data - use blocks , only : append_metablock - use blocks , only : set_last_id, get_last_id, get_mblocks, get_dblocks & - , get_nleafs - use coordinates, only : nn, ng, in, jn, kn, maxlev, toplev, ir, jr, kr - use coordinates, only : initialize_coordinates, finalize_coordinates - use coordinates, only : xmin, xmax, ymin, ymax, zmin, zmax - use error , only : print_error - use evolution, only : step, time, dt, dtn - use hdf5 , only : hid_t, hsize_t - use hdf5 , only : h5gopen_f, h5gclose_f, h5aget_num_attrs_f & - , h5aopen_idx_f, h5aclose_f, h5aget_name_f - use mpitools , only : nprocs, nproc - use random , only : nseeds, set_seeds + use blocks , only : block_meta + use blocks , only : append_metablock + use blocks , only : set_last_id, get_last_id + use blocks , only : get_mblocks, get_dblocks, get_nleafs + use coordinates , only : nn, ng, in, jn, kn, ir, jr, kr + use coordinates , only : maxlev, toplev + use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax + use coordinates , only : initialize_coordinates, finalize_coordinates + use error , only : print_error + use evolution , only : step, time, dt, dtn + use hdf5 , only : hid_t, hsize_t + use hdf5 , only : h5gopen_f, h5gclose_f + use mpitools , only : nprocs, nproc + use random , only : nseeds, set_seeds -! declare variables +! local variables are not implicit by default ! implicit none -! input variables +! subroutine arguments ! integer(hid_t), intent(in) :: fid ! local variables ! - character(len=16) :: aname - integer(hid_t) :: gid, aid - integer(hsize_t) :: alen = 16 - integer(kind=4) :: dm(3) - integer :: err, i, l - integer :: nattrs, lndims, llast_id, lmblocks, lnleafs & - , lncells, lnghost, lnseeds, lmaxlev, lnproc + integer(hid_t) :: gid + integer :: ierr, l + integer :: lndims, lmaxlev, lmblocks, lnleafs, llast_id + integer :: lncells, lnghost, lnproc, lnseeds ! local pointers ! type(block_meta), pointer :: pmeta - type(block_data), pointer :: pdata ! allocatable arrays ! @@ -1224,574 +1235,228 @@ module io ! ! open the global attributes group ! - call h5gopen_f(fid, 'attributes', gid, err) + call h5gopen_f(fid, 'attributes', gid, ierr) ! check if the group has been opened successfuly ! - if (err .ge. 0) then + if (ierr < 0) then + call print_error("io::read_attributes_h5", "Cannot open the group!") + return + end if -! read the number of global attributes +! restore integer attributes ! - call h5aget_num_attrs_f(gid, nattrs, err) + call read_attribute_integer_h5(gid, 'ndims' , lndims ) + call read_attribute_integer_h5(gid, 'maxlev' , lmaxlev ) + call read_attribute_integer_h5(gid, 'nprocs' , nfiles ) + call read_attribute_integer_h5(gid, 'nproc' , lnproc ) + call read_attribute_integer_h5(gid, 'mblocks', lmblocks) + call read_attribute_integer_h5(gid, 'nleafs' , lnleafs ) + call read_attribute_integer_h5(gid, 'last_id', llast_id) + call read_attribute_integer_h5(gid, 'ncells' , lncells ) + call read_attribute_integer_h5(gid, 'nghost' , lnghost ) + call read_attribute_integer_h5(gid, 'nseeds' , lnseeds ) + call read_attribute_integer_h5(gid, 'step' , step ) + call read_attribute_integer_h5(gid, 'isnap' , isnap ) -! check if the number of attributes has been read successfuly +! restore double precision attributes ! - if (err .ge. 0) then + call read_attribute_double_h5(gid, 'xmin', xmin) + call read_attribute_double_h5(gid, 'xmax', xmax) + call read_attribute_double_h5(gid, 'ymin', ymin) + call read_attribute_double_h5(gid, 'ymax', ymax) + call read_attribute_double_h5(gid, 'zmin', zmin) + call read_attribute_double_h5(gid, 'zmax', zmax) + call read_attribute_double_h5(gid, 'time', time) + call read_attribute_double_h5(gid, 'dt' , dt ) + call read_attribute_double_h5(gid, 'dtn' , dtn ) -! iterate over all attributes +! check the number of dimensions ! - do i = 0, nattrs - 1 + if (lndims /= NDIMS) then + call print_error("io::read_attributes_h5" & + , "The number of dimensions does not match!") + return + end if -! open the current attribute +! check the block dimensions ! - call h5aopen_idx_f(gid, i, aid, err) + if (lncells /= nn) then + call print_error("io::read_attributes_h5" & + , "The block dimensions do not match!") + end if -! check if the attribute has been opened successfuly +! check the number of ghost layers ! - if (err .ge. 0) then + if (lnghost /= ng) then + call print_error("io::read_attributes_h5" & + , "The number of ghost layers does not match!") + end if -! obtain the attribute name +! prepare coordinates and rescaling factors if the maximum level has changed ! - call h5aget_name_f(aid, alen, aname, err) - -! depending on the attribute name use proper subroutine to read its value -! - select case(trim(aname)) - case('ndims') - call read_attribute_integer_h5(aid, aname, lndims) - -! check if the restart file and compiled program have the same number of -! dimensions -! - if (lndims .ne. NDIMS) then - call print_error("io::read_attributes_h5" & - , "File and program dimensions are incompatible!") - end if - case('maxlev') - call read_attribute_integer_h5(aid, aname, lmaxlev) - if (lmaxlev .gt. toplev) then + if (lmaxlev > toplev) then ! subtitute the new value of toplev ! - toplev = lmaxlev + toplev = lmaxlev ! regenerate coordinates ! - call finalize_coordinates(err) - call initialize_coordinates(.false., err) + call finalize_coordinates(ierr) + call initialize_coordinates(.false., ierr) ! calculate a factor to rescale the block coordinates ! - dcor = 2**(toplev - maxlev) - - else - -! calculate a factor to rescale the block coordinates -! - ucor = 2**(maxlev - lmaxlev) - end if - case('nprocs') - call read_attribute_integer_h5(aid, aname, nfiles) - case('nproc') - call read_attribute_integer_h5(aid, aname, lnproc) - case('last_id') - call read_attribute_integer_h5(aid, aname, llast_id) - case('mblocks') - call read_attribute_integer_h5(aid, aname, lmblocks) - case('nleafs') - call read_attribute_integer_h5(aid, aname, lnleafs) - case('ncells') - call read_attribute_integer_h5(aid, aname, lncells) - -! check if the block dimensions are compatible -! - if (lncells .ne. nn) then - call print_error("io::read_attributes_h5" & - , "File and program block dimensions are incompatible!") - end if - case('nghost') - call read_attribute_integer_h5(aid, aname, lnghost) - -! check if the ghost layers are compatible -! - if (lnghost .ne. ng) then - call print_error("io::read_attributes_h5" & - , "File and program block ghost layers are incompatible!") - end if - case('step') - call read_attribute_integer_h5(aid, aname, step) - case('time') - call read_attribute_double_h5(aid, aname, time) - case('dt') - call read_attribute_double_h5(aid, aname, dt) - case('dtn') - call read_attribute_double_h5(aid, aname, dtn) - case('xmin') - call read_attribute_double_h5(aid, aname, xmin) - case('xmax') - call read_attribute_double_h5(aid, aname, xmax) - case('ymin') - call read_attribute_double_h5(aid, aname, ymin) - case('ymax') - call read_attribute_double_h5(aid, aname, ymax) - case('zmin') - call read_attribute_double_h5(aid, aname, zmin) - case('zmax') - call read_attribute_double_h5(aid, aname, zmax) - case('nseeds') - call read_attribute_integer_h5(aid, aname, lnseeds) - -! ! check if the numbers of seeds are compatible -! ! -! if (lnseeds .ne. nseeds) then -! call print_error("io::read_attributes_h5" & -! , "The number of seeds from file and program are incompatible!") -! end if - case('seeds') - -! check if the numbers of seeds are compatible -! - if (lnseeds .eq. nseeds) then - -! allocate space for seeds -! - allocate(seeds(nseeds)) - -! store them in the current group -! - call read_attribute_vector_integer_h5(aid, aname, nseeds & - , seeds(:)) - -! set the seed values -! - call set_seeds(nseeds, seeds(:)) - -! deallocate seed array -! - deallocate(seeds) - - end if - case default - end select - -! close the current attribute -! - call h5aclose_f(aid, err) - - else - -! print error about the problem with obtaining the number of attributes -! - call print_error("io::read_attributes_h5", & - "Cannot open the current attribute!") - - end if - - end do - -! allocate all metablocks -! - do l = 1, lmblocks - call append_metablock(pmeta) - end do - -! check if the number of created metablocks is equal to lbmcloks -! - if (lmblocks .ne. get_mblocks()) then - call print_error("io::read_attributes_h5" & - , "Number of metablocks doesn't match!") - end if - -! allocate an array of pointers with the size llast_id -! - allocate(block_array(llast_id)) - call set_last_id(llast_id) - - else - -! print error about the problem with obtaining the number of attributes -! - call print_error("io::read_attributes_h5", & - "Cannot get the number of global attributes!") - - end if - -! close the group -! - call h5gclose_f(gid, err) - -! check if the group has been closed successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the group -! - call print_error("io::read_attributes_h5", "Cannot close the group!") - - end if + dcor = 2**(toplev - maxlev) else -! print error about the problem with creating the group +! calculate a factor to rescale the block coordinates ! - call print_error("io::read_attributes_h5", "Cannot open the group!") + ucor = 2**(maxlev - lmaxlev) end if +! allocate space for seeds +! + allocate(seeds(lnseeds)) + +! store them in the current group +! + call read_attribute_vector_integer_h5(gid, 'seeds', seeds(:)) + +! set the seed values +! + call set_seeds(lnseeds, seeds(:)) + +! deallocate seed array +! + deallocate(seeds) + +! allocate all metablocks +! + do l = 1, lmblocks + call append_metablock(pmeta) + end do + +! check if the number of created metablocks is equal to lbmcloks +! + if (lmblocks /= get_mblocks()) then + call print_error("io::read_attributes_h5" & + , "Number of metablocks does not match!") + end if + +! allocate an array of pointers with the size llast_id +! + allocate(block_array(llast_id)) + +! set the last_id +! + call set_last_id(llast_id) + +! close the group +! + call h5gclose_f(gid, ierr) + +! check if the group has been closed successfuly +! + if (ierr /= 0) then + call print_error("io::read_attributes_h5", "Cannot close the group!") + end if + !------------------------------------------------------------------------------- ! end subroutine read_attributes_h5 ! !=============================================================================== -! -! read_datablock_dims_h5: subroutine reads the data block dimensions from the -! attributes group of the file given by the file -! identificator -! -! arguments: -! fid - the HDF5 file identificator -! dm - the data block dimensions +!! +!!--- ATTRIBUTE SUBROUTINES -------------------------------------------------- +!! +!=============================================================================== ! !=============================================================================== ! - subroutine read_datablock_dims_h5(fid, dm) - -! references to other modules +! subroutine READ_ATTRIBUTE_INTEGER_H5: +! ------------------------------------ ! - use error , only : print_error - use hdf5 , only : hid_t, hsize_t - use hdf5 , only : h5gopen_f, h5gclose_f, h5aget_num_attrs_f & - , h5aopen_idx_f, h5aclose_f, h5aget_name_f - use equations, only : nv +! Subroutine reads a value of the integer attribute provided by the group +! identifier to which it is linked and its name. +! +! Arguments: +! +! gid - the group identifier to which the attribute is linked; +! aname - the attribute name; +! avalue - the attribute value; +! +!=============================================================================== +! + subroutine read_attribute_integer_h5(gid, aname, avalue) -! declare variables +! import external procedures and variables +! + use error , only : print_error + use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER + use hdf5 , only : h5aexists_by_name_f + use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f + +! local variables are not implicit by default ! implicit none -! input variables +! attribute arguments ! - integer(hid_t) , intent(in) :: fid - integer(hsize_t), dimension(5), intent(out) :: dm + integer(hid_t) , intent(in) :: gid + character(len=*), intent(in) :: aname + integer , intent(inout) :: avalue ! local variables ! - character(len=16) :: aname - integer(hid_t) :: gid, aid - integer(hsize_t) :: alen = 16 - integer :: err, i - integer :: nattrs, ldblocks, lnghost - -! local arrays -! - integer(kind=4), dimension(3) :: lm + logical :: exists = .false. + integer(hid_t) :: aid + integer(hsize_t), dimension(1) :: am = (/ 1 /) + integer :: ierr ! !------------------------------------------------------------------------------- ! -! initiate the output vector +! check if the attribute exists in the group provided by gid ! - dm(:) = 0 - dm(2) = nv - -! open the global attributes group -! - call h5gopen_f(fid, 'attributes', gid, err) - -! check if the group has been opened successfuly -! - if (err .ge. 0) then - -! read the number of global attributes -! - call h5aget_num_attrs_f(gid, nattrs, err) - -! check if the number of attributes has been read successfuly -! - if (err .ge. 0) then - -! iterate over all attributes -! - do i = 0, nattrs - 1 - -! open the current attribute -! - call h5aopen_idx_f(gid, i, aid, err) - -! check if the attribute has been opened successfuly -! - if (err .ge. 0) then - -! obtain the attribute name -! - call h5aget_name_f(aid, alen, aname, err) - -! check if the attribute name has been read successfuly -! - if (err .ge. 0) then - -! depending on the attribute name use proper subroutine to read its value -! - select case(trim(aname)) - case('dblocks') - -! obtain the number of data blocks -! - call read_attribute_integer_h5(aid, aname, ldblocks) - - case('dims') - -! obtain the block dimensions -! - call read_attribute_vector_integer_h5(aid, aname, 3, lm(:)) - - case('nghost') - -! obtain the number of data blocks -! - call read_attribute_integer_h5(aid, aname, lnghost) - - case default - end select - - else - -! print error about the problem with reading the attribute name -! - call print_error("io::read_datablock_dims_h5", & - "Cannot read the current attribute name!") - - end if - -! close the current attribute -! - call h5aclose_f(aid, err) - - else - -! print error about the problem with opening the current attribute -! - call print_error("io::read_datablock_dims_h5", & - "Cannot open the current attribute!") - - end if - - end do ! i = 0, nattrs - 1 - -! prepare the output array -! - dm(1) = ldblocks - do i = 1, 3 - if (lm(i) .gt. 1) lm(i) = lm(i) + 2 * lnghost - end do - dm(3:5) = lm(1:3) - - else - -! print error about the problem with obtaining the number of attributes -! - call print_error("io::read_datablock_dims_h5", & - "Cannot get the number of global attributes!") - - end if - -! close the group -! - call h5gclose_f(gid, err) - -! check if the group has been closed successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the group -! - call print_error("io::read_datablock_dims_h5" & - , "Cannot close the attributes group!") - - end if - - else - -! print error about the problem with creating the group -! - call print_error("io::read_datablock_dims_h5" & - , "Cannot open the attributes group!") - + call h5aexists_by_name_f(gid, '.', aname, exists, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_integer_h5" & + , "Cannot check if attribute exists :" // trim(aname)) + return + end if + if (.not. exists) then + call print_error("io::read_attribute_integer_h5" & + , "Attribute does not exist :" // trim(aname)) + return end if -!------------------------------------------------------------------------------- +! open the attribute ! - end subroutine read_datablock_dims_h5 -! -!=============================================================================== -! -! write_attribute_integer_h5: subroutine writes an attribute with the integer -! value in a group given by its indentificator -! -! arguments: -! gid - the HDF5 group identificator -! name - the string name of the attribute -! value - the attribute value -! -!=============================================================================== -! - subroutine write_attribute_integer_h5(gid, name, value) + call h5aopen_by_name_f(gid, '.', aname, aid, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_integer_h5" & + , "Cannot open attribute :" // trim(aname)) + return + end if -! references to other modules +! read attribute value ! - use error, only : print_error - use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER - use hdf5 , only : h5screate_simple_f, h5sclose_f & - , h5acreate_f, h5awrite_f, h5aclose_f - -! declare variables -! - implicit none - -! input variables -! - integer(hid_t) , intent(in) :: gid - character(len=*), intent(in) :: name - integer , intent(in) :: value - -! local variables -! - integer(hid_t) :: sid, aid - integer(hsize_t), dimension(1) :: am = (/ 1 /) - integer :: err -! -!------------------------------------------------------------------------------- -! -! create space for the attribute -! - call h5screate_simple_f(1, am, sid, err) - -! check if the space has been created successfuly -! - if (err .ge. 0) then - -! create the attribute -! - call h5acreate_f(gid, name, H5T_NATIVE_INTEGER, sid, aid, err) - -! check if the attribute has been created successfuly -! - if (err .ge. 0) then - -! write the attribute data -! - call h5awrite_f(aid, H5T_NATIVE_INTEGER, value, am, err) - -! check if the attribute data has been written successfuly -! - if (err .gt. 0) then - -! print error about the problem with storing the attribute data -! - call print_error("io::write_attribute_integer_h5" & - , "Cannot write the attribute data in :" // trim(name)) - - end if + call h5aread_f(aid, H5T_NATIVE_INTEGER, avalue, am(:), ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_integer_h5" & + , "Cannot read attribute :" // trim(aname)) + end if ! close the attribute ! - call h5aclose_f(aid, err) - -! check if the attribute has been closed successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the attribute -! - call print_error("io::write_attribute_integer_h5" & - , "Cannot close attribute :" // trim(name)) - - end if - - else - -! print error about the problem with creating the attribute -! - call print_error("io::write_attribute_integer_h5" & - , "Cannot create attribute :" // trim(name)) - - end if - -! release the space -! - call h5sclose_f(sid, err) - -! check if the space has been released successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the space -! - call print_error("io::write_attribute_integer_h5" & - , "Cannot close space for attribute :" // trim(name)) - - end if - - else - -! print error about the problem with creating the space for the attribute -! - call print_error("io::write_attribute_integer_h5" & - , "Cannot create space for attribute :" // trim(name)) - - end if - -!------------------------------------------------------------------------------- -! - end subroutine write_attribute_integer_h5 -! -!=============================================================================== -! -! read_attribute_integer_h5: subroutine read an integer value from an attribute -! given by the identificator -! -! arguments: -! aid - the HDF5 attribute identificator -! value - the attribute value -! -!=============================================================================== -! - subroutine read_attribute_integer_h5(aid, name, value) - -! references to other modules -! - use error, only : print_error - use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER - use hdf5 , only : h5aread_f - -! declare variables -! - implicit none - -! input variables -! - integer(hid_t) , intent(in) :: aid - character(len=*), intent(in) :: name - integer , intent(inout) :: value - -! local variables -! - integer(hsize_t), dimension(1) :: am = (/ 1 /) - integer :: err -! -!------------------------------------------------------------------------------- -! -! read attribute value -! - call h5aread_f(aid, H5T_NATIVE_INTEGER, value, am(:), err) - -! check if the attribute has been read successfuly -! - if (err .gt. 0) then - -! print error about the problem with reading the attribute -! + call h5aclose_f(aid, ierr) + if (ierr /= 0) then call print_error("io::read_attribute_integer_h5" & - , "Cannot read attribute :" // trim(name)) - + , "Cannot close attribute :" // trim(aname)) + return end if !------------------------------------------------------------------------------- @@ -1800,377 +1465,86 @@ module io ! !=============================================================================== ! -! write_attribute_vector_integer_h5: subroutine writes an vector intiger -! attribute in a group given by its indentificator +! subroutine READ_ATTRIBUTE_DOUBLE_H5: +! ----------------------------------- ! -! arguments: -! gid - the HDF5 group identificator -! name - the string name of the attribute -! length - the vector length -! data - the attribute value +! Subroutine reads a value of the double precision attribute provided by +! the group identifier to which it is linked and its name. +! +! Arguments: +! +! gid - the group identifier to which the attribute is linked; +! aname - the attribute name; +! avalue - the attribute value; ! !=============================================================================== ! - subroutine write_attribute_vector_integer_h5(gid, name, length, data) + subroutine read_attribute_double_h5(gid, aname, avalue) -! references to other modules +! import external procedures and variables ! - use error, only : print_error - use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER - use hdf5 , only : h5screate_simple_f, h5sclose_f & - , h5acreate_f, h5awrite_f, h5aclose_f + use error , only : print_error + use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_DOUBLE + use hdf5 , only : h5aexists_by_name_f + use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f -! declare variables +! local variables are not implicit by default ! implicit none -! input variables +! attribute arguments ! - integer(hid_t) , intent(in) :: gid - character(len=*) , intent(in) :: name - integer , intent(in) :: length - integer(kind=4) , dimension(:), intent(in) :: data + integer(hid_t) , intent(in) :: gid + character(len=*), intent(in) :: aname + real(kind=8) , intent(inout) :: avalue ! local variables ! - integer(hid_t) :: sid, aid - integer(hsize_t), dimension(1) :: am - integer :: err + logical :: exists = .false. + integer(hid_t) :: aid + integer(hsize_t), dimension(1) :: am = (/ 1 /) + integer :: ierr ! !------------------------------------------------------------------------------- ! -! prepare the space dimensions -! - am(1) = length - -! create space for the attribute -! - call h5screate_simple_f(1, am, sid, err) - -! check if the space has been created successfuly -! - if (err .ge. 0) then - -! create the attribute -! - call h5acreate_f(gid, name, H5T_NATIVE_INTEGER, sid, aid, err) - -! check if the attribute has been created successfuly -! - if (err .ge. 0) then - -! write the attribute data -! - call h5awrite_f(aid, H5T_NATIVE_INTEGER, data, am, err) - -! check if the attribute data has been written successfuly -! - if (err .gt. 0) then - -! print error about the problem with storing the attribute data -! - call print_error("io::write_attribute_vector_integer_h5" & - , "Cannot write the attribute data in :" // trim(name)) - - end if - -! close the attribute -! - call h5aclose_f(aid, err) - -! check if the attribute has been closed successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the attribute -! - call print_error("io::write_attribute_vector_integer_h5" & - , "Cannot close attribute :" // trim(name)) - - end if - - else - -! print error about the problem with creating the attribute -! - call print_error("io::write_attribute_vector_integer_h5" & - , "Cannot create attribute :" // trim(name)) - - end if - -! release the space -! - call h5sclose_f(sid, err) - -! check if the space has been released successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the space -! - call print_error("io::write_attribute_vector_integer_h5" & - , "Cannot close space for attribute :" // trim(name)) - - end if - - else - -! print error about the problem with creating the space for the attribute -! - call print_error("io::write_attribute_vector_integer_h5" & - , "Cannot create space for attribute :" // trim(name)) - - end if - -!------------------------------------------------------------------------------- -! - end subroutine write_attribute_vector_integer_h5 -! -!=============================================================================== -! -! read_attribute_vector_integer_h5: subroutine reads a vector of integer values -! from an attribute given by the identificator -! -! arguments: -! aid - the HDF5 attribute identificator -! name - the attribute name -! length - the vector length -! value - the attribute value -! -!=============================================================================== -! - subroutine read_attribute_vector_integer_h5(aid, name, length, value) - -! references to other modules -! - use error, only : print_error - use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER - use hdf5 , only : h5aread_f - -! declare variables -! - implicit none - -! input variables -! - integer(hid_t) , intent(in) :: aid - character(len=*), intent(in) :: name - integer , intent(in) :: length - integer, dimension(:), intent(inout) :: value - -! local variables -! - integer(hsize_t), dimension(1) :: am - integer :: err -! -!------------------------------------------------------------------------------- -! -! check if the length is larger than 0 -! - if (length .gt. 0) then - -! prepare dimension array -! - am(1) = length - -! read attribute value -! - call h5aread_f(aid, H5T_NATIVE_INTEGER, value(:), am(:), err) - -! check if the attribute has been read successfuly -! - if (err .gt. 0) then - -! print error about the problem with reading the attribute -! - call print_error("io::read_attribute_vector_integer_h5" & - , "Cannot read attribute :" // trim(name)) - - end if - - else ! length > 0 - -! print error about the wrong vector size -! - call print_error("io::read_attribute_vector_integer_h5" & - , "Wrong length of vector") - - end if ! length > 0 - -!------------------------------------------------------------------------------- -! - end subroutine read_attribute_vector_integer_h5 -! -!=============================================================================== -! -! write_attribute_double_h5: subroutine writes a double precision attribute in -! a group given by its indentificator -! -! arguments: -! gid - the HDF5 group identificator -! name - the string name of the attribute -! value - the attribute value -! -!=============================================================================== -! - subroutine write_attribute_double_h5(gid, name, value) - -! references to other modules -! - use error, only : print_error - use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_DOUBLE - use hdf5 , only : h5screate_simple_f, h5sclose_f, h5acreate_f, h5awrite_f & - , h5aclose_f - -! declare variables -! - implicit none - -! input variables -! - integer(hid_t) , intent(in) :: gid - character(len=*), intent(in) :: name - real(kind=8) , intent(in) :: value - -! local variables -! - integer(hid_t) :: sid, aid - integer(hsize_t), dimension(1) :: am = (/ 1 /) - integer :: err -! -!------------------------------------------------------------------------------- -! -! create space for the attribute -! - call h5screate_simple_f(1, am, sid, err) - -! check if the space has been created successfuly -! - if (err .ge. 0) then - -! create the attribute -! - call h5acreate_f(gid, name, H5T_NATIVE_DOUBLE, sid, aid, err) - -! check if the attribute has been created successfuly -! - if (err .ge. 0) then - -! write the attribute data -! - call h5awrite_f(aid, H5T_NATIVE_DOUBLE, value, am, err) - -! check if the attribute data has been written successfuly -! - if (err .gt. 0) then - -! print error about the problem with storing the attribute data -! - call print_error("io::write_attribute_double_h5" & - , "Cannot write the attribute data in :" // trim(name)) - - end if - -! close the attribute -! - call h5aclose_f(aid, err) - -! check if the attribute has been closed successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the attribute -! - call print_error("io::write_attribute_double_h5" & - , "Cannot close attribute :" // trim(name)) - - end if - - else - -! print error about the problem with creating the attribute -! - call print_error("io::write_attribute_double_h5" & - , "Cannot create attribute :" // trim(name)) - - end if - -! release the space -! - call h5sclose_f(sid, err) - -! check if the space has been released successfuly -! - if (err .gt. 0) then - -! print error about the problem with closing the space -! - call print_error("io::write_attribute_double_h5" & - , "Cannot close space for attribute :" // trim(name)) - - end if - - else - -! print error about the problem with creating the space for the attribute -! - call print_error("io::write_attribute_double_h5" & - , "Cannot create space for attribute :" // trim(name)) - - end if - -!------------------------------------------------------------------------------- -! - end subroutine write_attribute_double_h5 -! -!=============================================================================== -! -! read_attribute_double_h5: subroutine reads a double precision attribute -! -! arguments: -! aid - the HDF5 attribute identificator -! value - the attribute value -! -!=============================================================================== -! - subroutine read_attribute_double_h5(aid, name, value) - -! references to other modules -! - use error, only : print_error - use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_DOUBLE - use hdf5 , only : h5aread_f - -! declare variables -! - implicit none - -! input variables -! - integer(hid_t) , intent(in) :: aid - character(len=*), intent(in) :: name - real(kind=8) , intent(inout) :: value - -! local variables -! - integer(hsize_t), dimension(1) :: am = (/ 1 /) - integer :: err -! -!------------------------------------------------------------------------------- -! -! read attribute value -! - call h5aread_f(aid, H5T_NATIVE_DOUBLE, value, am(:), err) - -! check if the attribute has been read successfuly -! - if (err .gt. 0) then - -! print error about the problem with reading the attribute +! check if the attribute exists in the group provided by gid ! + call h5aexists_by_name_f(gid, '.', aname, exists, ierr) + if (ierr /= 0) then call print_error("io::read_attribute_double_h5" & - , "Cannot read attribute :" // trim(name)) + , "Cannot check if attribute exists :" // trim(aname)) + return + end if + if (.not. exists) then + call print_error("io::read_attribute_double_h5" & + , "Attribute does not exist :" // trim(aname)) + return + end if +! open the attribute +! + call h5aopen_by_name_f(gid, '.', aname, aid, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_double_h5" & + , "Cannot open attribute :" // trim(aname)) + return + end if + +! read attribute value +! + call h5aread_f(aid, H5T_NATIVE_DOUBLE, avalue, am(:), ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_double_h5" & + , "Cannot read attribute :" // trim(aname)) + end if + +! close the attribute +! + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_double_h5" & + , "Cannot close attribute :" // trim(aname)) + return end if !------------------------------------------------------------------------------- @@ -2179,13 +1553,405 @@ module io ! !=============================================================================== ! +! subroutine READ_ATTRIBUTE_VECTOR_INTEGER_H5: +! ------------------------------------------- +! +! Subroutine reads a vector of the integer attribute provided by the group +! identifier to which it is linked and its name. +! +! Arguments: +! +! gid - the group identifier to which the attribute is linked; +! aname - the attribute name; +! avalue - the attribute value; +! +!=============================================================================== +! + subroutine read_attribute_vector_integer_h5(gid, aname, avalue) + +! import external procedures and variables +! + use error , only : print_error + use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER + use hdf5 , only : h5aexists_by_name_f, h5aget_space_f + use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f + use hdf5 , only : h5sclose_f, h5sget_simple_extent_dims_f + +! local variables are not implicit by default +! + implicit none + +! attribute arguments +! + integer(hid_t) , intent(in) :: gid + character(len=*) , intent(in) :: aname + integer , dimension(:), intent(inout) :: avalue + +! local variables +! + logical :: exists = .false. + integer(hid_t) :: aid, sid + integer(hsize_t), dimension(1) :: am, bm + integer(hsize_t) :: alen + integer :: ierr +! +!------------------------------------------------------------------------------- +! +! check if the attribute exists in the group provided by gid +! + call h5aexists_by_name_f(gid, '.', aname, exists, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Cannot check if attribute exists :" // trim(aname)) + return + end if + if (.not. exists) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Attribute does not exist :" // trim(aname)) + return + end if + +! open the attribute +! + call h5aopen_by_name_f(gid, '.', aname, aid, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Cannot open attribute :" // trim(aname)) + return + end if + +! get the attribute space +! + call h5aget_space_f(aid, sid, ierr) + if (ierr == 0) then + call h5sget_simple_extent_dims_f(sid, am, bm, ierr) + if (ierr /= 1) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Cannot get attribute dimensions :" // trim(aname)) + end if + call h5sclose_f(sid, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Cannot close the attribute space :" // trim(aname)) + end if + else + call print_error("io::read_attribute_vector_integer_h5" & + , "Cannot get the attribute space :" // trim(aname)) + return + end if + +! check if the output array is large enough +! + if (am(1) > size(avalue)) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Attribute too large for output argument :" // trim(aname)) + return + end if + +! read attribute value +! + call h5aread_f(aid, H5T_NATIVE_INTEGER, avalue, am(:), ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Cannot read attribute :" // trim(aname)) + end if + +! close the attribute +! + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + call print_error("io::read_attribute_vector_integer_h5" & + , "Cannot close attribute :" // trim(aname)) + return + end if + +!------------------------------------------------------------------------------- +! + end subroutine read_attribute_vector_integer_h5 +! +!=============================================================================== +! +! subroutine WRITE_ATTRIBUTE_INTEGER_H5: +! ------------------------------------- +! +! Subroutine stores a value of the integer attribute in the group provided +! by an identifier and the attribute name. +! +! Arguments: +! +! gid - the group identifier to which the attribute should be linked; +! aname - the attribute name; +! avalue - the attribute value; +! +!=============================================================================== +! + subroutine write_attribute_integer_h5(gid, aname, avalue) + +! import external procedures and variables +! + use error , only : print_error + use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER + use hdf5 , only : h5screate_simple_f, h5sclose_f + use hdf5 , only : h5acreate_f, h5awrite_f, h5aclose_f + +! local variables are not implicit by default +! + implicit none + +! attribute arguments +! + integer(hid_t) , intent(in) :: gid + character(len=*), intent(in) :: aname + integer , intent(in) :: avalue + +! local variables +! + integer(hid_t) :: sid, aid + integer(hsize_t), dimension(1) :: am = (/ 1 /) + integer :: ierr +! +!------------------------------------------------------------------------------- +! +! create space for the attribute value +! + call h5screate_simple_f(1, am, sid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_integer_h5" & + , "Cannot create space for attribute :" // trim(aname)) + return + end if + +! create the attribute in the given group +! + call h5acreate_f(gid, aname, H5T_NATIVE_INTEGER, sid, aid, ierr) + if (ierr == 0) then + +! write the attribute data +! + call h5awrite_f(aid, H5T_NATIVE_INTEGER, avalue, am, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_integer_h5" & + , "Cannot write the attribute data in :" // trim(aname)) + end if + +! close the attribute +! + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_integer_h5" & + , "Cannot close attribute :" // trim(aname)) + end if + + else + call print_error("io::write_attribute_integer_h5" & + , "Cannot create attribute :" // trim(aname)) + end if + +! release the space +! + call h5sclose_f(sid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_integer_h5" & + , "Cannot close space for attribute :" // trim(aname)) + end if + +!------------------------------------------------------------------------------- +! + end subroutine write_attribute_integer_h5 +! +!=============================================================================== +! +! subroutine WRITE_ATTRIBUTE_DOUBLE_H5: +! ------------------------------------ +! +! Subroutine stores a value of the double precision attribute in the group +! provided by an identifier and the attribute name. +! +! Arguments: +! +! gid - the group identifier to which the attribute should be linked; +! aname - the attribute name; +! avalue - the attribute value; +! +!=============================================================================== +! + subroutine write_attribute_double_h5(gid, aname, avalue) + +! import external procedures and variables +! + use error , only : print_error + use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_DOUBLE + use hdf5 , only : h5screate_simple_f, h5sclose_f + use hdf5 , only : h5acreate_f, h5awrite_f, h5aclose_f + +! local variables are not implicit by default +! + implicit none + +! attribute arguments +! + integer(hid_t) , intent(in) :: gid + character(len=*), intent(in) :: aname + real(kind=8) , intent(in) :: avalue + +! local variables +! + integer(hid_t) :: sid, aid + integer(hsize_t), dimension(1) :: am = (/ 1 /) + integer :: ierr +! +!------------------------------------------------------------------------------- +! +! create space for the attribute value +! + call h5screate_simple_f(1, am, sid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_double_h5" & + , "Cannot create space for attribute :" // trim(aname)) + return + end if + +! create the attribute in the given group +! + call h5acreate_f(gid, aname, H5T_NATIVE_DOUBLE, sid, aid, ierr) + if (ierr == 0) then + +! write the attribute data +! + call h5awrite_f(aid, H5T_NATIVE_DOUBLE, avalue, am, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_double_h5" & + , "Cannot write the attribute data in :" // trim(aname)) + end if + +! close the attribute +! + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_double_h5" & + , "Cannot close attribute :" // trim(aname)) + end if + + else + call print_error("io::write_attribute_double_h5" & + , "Cannot create attribute :" // trim(aname)) + end if + +! release the space +! + call h5sclose_f(sid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_double_h5" & + , "Cannot close space for attribute :" // trim(aname)) + end if + +!------------------------------------------------------------------------------- +! + end subroutine write_attribute_double_h5 +! +!=============================================================================== +! +! subroutine WRITE_ATTRIBUTE_VECTOR_INTEGER_H5: +! -------------------------------------------- +! +! Subroutine stores a vector of the integer attribute in the group provided +! by an identifier and the attribute name. +! +! Arguments: +! +! gid - the group identifier to which the attribute should be linked; +! aname - the attribute name; +! avalue - the attribute values; +! +!=============================================================================== +! + subroutine write_attribute_vector_integer_h5(gid, aname, avalue) + +! import external procedures and variables +! + use error , only : print_error + use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_INTEGER + use hdf5 , only : h5screate_simple_f, h5sclose_f + use hdf5 , only : h5acreate_f, h5awrite_f, h5aclose_f + +! local variables are not implicit by default +! + implicit none + +! attribute arguments +! + integer(hid_t) , intent(in) :: gid + character(len=*) , intent(in) :: aname + integer , dimension(:), intent(in) :: avalue + +! local variables +! + integer(hid_t) :: sid, aid + integer(hsize_t), dimension(1) :: am = (/ 1 /) + integer :: ierr +! +!------------------------------------------------------------------------------- +! +! set the proper attribute length +! + am(1) = size(avalue) + +! create space for the attribute value +! + call h5screate_simple_f(1, am, sid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_vector_integer_h5" & + , "Cannot create space for attribute :" // trim(aname)) + return + end if + +! create the attribute in the given group +! + call h5acreate_f(gid, aname, H5T_NATIVE_INTEGER, sid, aid, ierr) + if (ierr == 0) then + +! write the attribute data +! + call h5awrite_f(aid, H5T_NATIVE_INTEGER, avalue, am, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_vector_integer_h5" & + , "Cannot write the attribute data in :" // trim(aname)) + end if + +! close the attribute +! + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_vector_integer_h5" & + , "Cannot close attribute :" // trim(aname)) + end if + + else + call print_error("io::write_attribute_vector_integer_h5" & + , "Cannot create attribute :" // trim(aname)) + end if + +! release the space +! + call h5sclose_f(sid, ierr) + if (ierr /= 0) then + call print_error("io::write_attribute_vector_integer_h5" & + , "Cannot close space for attribute :" // trim(aname)) + end if + +!------------------------------------------------------------------------------- +! + end subroutine write_attribute_vector_integer_h5 +! +!=============================================================================== +! ! write_metablocks_h5: subroutine writes metablocks in the HDF5 format connected -! to the provided identificator +! to the provided identifier ! ! info: this subroutine stores only the metablocks ! ! arguments: -! fid - the HDF5 file identificator +! fid - the HDF5 file identifier ! !=============================================================================== ! @@ -2411,7 +2177,7 @@ module io ! info: this subroutine restores metablocks only ! ! arguments: -! fid - the HDF5 file identificator +! fid - the HDF5 file identifier ! !=============================================================================== ! @@ -2555,7 +2321,7 @@ module io block_array(id(l))%ptr => pmeta call metablock_set_id (pmeta, id (l)) - call metablock_set_process (pmeta, min(lcpu, cpu(l))) + call metablock_set_process (pmeta, cpu(l)) call metablock_set_refinement (pmeta, ref(l)) call metablock_set_configuration(pmeta, cfg(l)) call metablock_set_level (pmeta, lev(l)) @@ -2655,7 +2421,7 @@ module io ! ! Arguments: ! -! fid - the HDF5 file identificator; +! fid - the HDF5 file identifier; ! !=============================================================================== ! @@ -2725,7 +2491,7 @@ module io dm(4) = jm dm(5) = km -! allocate arrays to store associated meta block identificators, conserved and +! allocate arrays to store associated meta block identifiers, conserved and ! primitive variables ! allocate(id(am(1))) @@ -2841,7 +2607,7 @@ module io ! ! Arguments: ! -! fid - the HDF5 file identificator; +! fid - the HDF5 file identifier; ! !=============================================================================== ! @@ -2852,6 +2618,7 @@ module io use blocks , only : block_meta, block_data, list_data use blocks , only : append_datablock, link_blocks use coordinates , only : im, jm, km + use equations , only : nv use error , only : print_error use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5gopen_f, h5gclose_f @@ -2872,7 +2639,7 @@ module io ! integer(hid_t) :: gid integer(kind=4) :: l - integer :: err + integer :: dblocks, ierr ! local arrays ! @@ -2885,84 +2652,90 @@ module io ! !------------------------------------------------------------------------------- ! -! get datablock array dimensions +! read the number of data blocks ! - call read_datablock_dims_h5(fid, dm(:)) - -! open the group 'datablocks' -! - call h5gopen_f(fid, 'datablocks', gid, err) - -! check if the datablock group has been opened successfuly -! - if (err >= 0) then + call h5gopen_f(fid, 'attributes', gid, ierr) + if (ierr /= 0) then + call print_error("io::read_datablocks_h5" & + , "Cannot open the attribute group!") + return + end if + call read_attribute_integer_h5(gid, 'dblocks', dblocks) + call h5gclose_f(gid, ierr) + if (ierr /= 0) then + call print_error("io::read_datablocks_h5" & + , "Cannot close the attribute group!") + return + end if ! restore data blocks only if there are any ! - if (dm(1) > 0) then + if (dblocks > 0) then + +! fill out dimensions dm(:) +! + dm(1) = dblocks + dm(2) = nv + dm(3) = im + dm(4) = jm + dm(5) = km ! allocate arrays to read data ! - allocate(id(dm(1))) - allocate(uv(dm(1),dm(2),dm(3),dm(4),dm(5))) - allocate(qv(dm(1),dm(2),dm(3),dm(4),dm(5))) + allocate(id(dm(1))) + allocate(uv(dm(1),dm(2),dm(3),dm(4),dm(5))) + allocate(qv(dm(1),dm(2),dm(3),dm(4),dm(5))) + +! open the group 'datablocks' +! + call h5gopen_f(fid, 'datablocks', gid, ierr) + if (ierr /= 0) then + call print_error("io::read_datablocks_h5" & + , "Cannot open the data block group!") + return + end if ! read array data from the HDF5 file ! - call read_vector_integer_h5(gid, 'meta', dm(1:1), id(:) ) - call read_array5_double_h5 (gid, 'uvar', dm(1:5), uv(:,:,:,:,:)) - call read_array5_double_h5 (gid, 'qvar', dm(1:5), qv(:,:,:,:,:)) + call read_vector_integer_h5(gid, 'meta', dm(1:1), id(:) ) + call read_array5_double_h5 (gid, 'uvar', dm(1:5), uv(:,:,:,:,:)) + call read_array5_double_h5 (gid, 'qvar', dm(1:5), qv(:,:,:,:,:)) + +! close the data block group +! + call h5gclose_f(gid, ierr) + if (ierr /= 0) then + call print_error("io::read_datablocks_h5" & + , "Cannot close the data block group!") + return + end if ! iterate over data blocks, allocate them and fill out their fields ! - do l = 1, dm(1) + do l = 1, dm(1) ! allocate and append to the end of the list a new datablock ! - call append_datablock(pdata) + call append_datablock(pdata) ! associate the corresponding meta block with the current data block ! - call link_blocks(block_array(id(l))%ptr, pdata) + call link_blocks(block_array(id(l))%ptr, pdata) ! fill out the array of conservative and primitive variables ! - pdata%u(:,:,:,:) = uv(l,:,:,:,:) - pdata%q(:,:,:,:) = qv(l,:,:,:,:) + pdata%u(:,:,:,:) = uv(l,:,:,:,:) + pdata%q(:,:,:,:) = qv(l,:,:,:,:) - end do ! l = 1, dm(1) + end do ! l = 1, dm(1) ! deallocate allocatable arrays ! - if (allocated(id)) deallocate(id) - if (allocated(uv)) deallocate(uv) - if (allocated(qv)) deallocate(qv) + if (allocated(id)) deallocate(id) + if (allocated(uv)) deallocate(uv) + if (allocated(qv)) deallocate(qv) - end if ! dblocks > 0 - -! close the group -! - call h5gclose_f(gid, err) - -! check if the group has been closed successfuly -! - if (err > 0) then - -! print error about the problem with closing the group -! - call print_error("io::read_datablocks_h5" & - , "Cannot close data block group!") - - end if - - else - -! print error about the problem with opening the group -! - call print_error("io::read_datablocks_h5" & - , "Cannot open data block group!") - - end if + end if ! dblocks > 0 !------------------------------------------------------------------------------- ! @@ -2977,7 +2750,7 @@ module io ! info: this subroutine stores coordinates ! ! arguments: -! fid - the HDF5 file identificator +! fid - the HDF5 file identifier; ! !=============================================================================== ! @@ -3139,11 +2912,11 @@ module io ! ! Subroutine groups each primitive variable from all data blocks and writes ! it as an array in the HDF5 dataset connected to the input HDF file -! identificator. +! identifier. ! ! Arguments: ! -! fid - the HDF5 file identificator; +! fid - the HDF5 file identifier; ! !=============================================================================== ! @@ -3311,11 +3084,11 @@ module io ! ! Subroutine groups each conservative variable from all data blocks and writes ! it as an array in the HDF5 dataset connected to the input HDF file -! identificator. +! identifier. ! ! Arguments: ! -! fid - the HDF5 file identificator; +! fid - the HDF5 file identifier; ! !=============================================================================== ! @@ -3481,7 +3254,7 @@ module io ! write_vector_integer_h5: subroutine stores a 1D integer vector in a group ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! length - the vector length ! value - the data @@ -3608,7 +3381,7 @@ module io ! read_vector_integer_h5: subroutine reads a 1D integer vector ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! length - the vector length ! value - the data @@ -3697,7 +3470,7 @@ module io ! write_array2_integer_h5: subroutine stores a 2D integer array in a group ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -3907,7 +3680,7 @@ module io ! read_array2_integer_h5: subroutine reads a 2D integer array ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -3996,7 +3769,7 @@ module io ! write_array4_integer_h5: subroutine stores a 4D integer array in a group ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -4206,7 +3979,7 @@ module io ! read_array4_integer_h5: subroutine reads a 4D integer array ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -4296,7 +4069,7 @@ module io ! a group ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! !=============================================================================== ! @@ -4420,7 +4193,7 @@ module io ! read_vector_double_h5: subroutine reads a 1D double precision vector ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the vector dimensions ! value - the data @@ -4509,7 +4282,7 @@ module io ! write_array4_float_h5: subroutine stores a 4D single precision array ! ! arguments: -! gid - the HDF5 group identificator where the dataset should be located +! gid - the HDF5 group identifier where the dataset should be located ! name - the string name representing the dataset ! dm - the dataset dimensions ! value - the dataset values @@ -4663,7 +4436,7 @@ module io ! write_array3_double_h5: subroutine stores a 3D double precision array ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -4873,7 +4646,7 @@ module io ! write_array4_double_h5: subroutine stores a 4D double precision array ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -5083,7 +4856,7 @@ module io ! write_array5_double_h5: subroutine stores a 5D double precision array ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -5293,7 +5066,7 @@ module io ! read_array5_double_h5: subroutine reads a 5D double precision array ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data @@ -5382,7 +5155,7 @@ module io ! write_array6_double_h5: subroutine stores a 6D double precision array ! ! arguments: -! gid - the HDF5 group identificator +! gid - the HDF5 group identifier ! name - the string name representing the dataset ! dm - the data dimensions ! value - the data diff --git a/src/mesh.F90 b/src/mesh.F90 index 236e14d..baa8ac2 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -462,6 +462,9 @@ module mesh use blocks , only : link_blocks, unlink_blocks, refine_block use blocks , only : get_mblocks, get_nleafs use blocks , only : set_neighbors_refine +#ifdef DEBUG + use blocks , only : check_neighbors +#endif /* DEBUG */ use coordinates , only : minlev, maxlev use domains , only : setup_domain use error , only : print_error @@ -735,6 +738,12 @@ module mesh end do ! pmeta +#ifdef DEBUG +! check if neighbors are consistent after mesh generation +! + call check_neighbors() +#endif /* DEBUG */ + #ifdef PROFILE ! stop accounting time for the initial mesh generation ! @@ -768,6 +777,9 @@ module mesh use blocks , only : refine_block, derefine_block use blocks , only : append_datablock, remove_datablock, link_blocks use blocks , only : set_neighbors_refine +#ifdef DEBUG + use blocks , only : check_neighbors +#endif /* DEBUG */ use coordinates , only : minlev, maxlev, toplev, im, jm, km use equations , only : nv use error , only : print_error @@ -819,12 +831,6 @@ module mesh call start_timer(imu) #endif /* PROFILE */ -#ifdef DEBUG -! check the mesh when debugging -! - call check_mesh('before update_mesh') -#endif /* DEBUG */ - !! DETERMINE THE REFINEMENT OF ALL DATA BLOCKS !! ! set the pointer to the first block on the data block list @@ -1264,9 +1270,9 @@ module mesh #endif /* MPI */ #ifdef DEBUG -! check mesh +! check if neighbors are consistent after mesh refinement ! - call check_mesh('after update_mesh') + call check_neighbors() #endif /* DEBUG */ #ifdef PROFILE @@ -1359,92 +1365,98 @@ module mesh ! do while (associated(pmeta)) +! consider only meta blocks which belong to active processes +! + if (pmeta%process < nprocs) then + ! check if the block belongs to another process ! - if (pmeta%process /= np) then + if (pmeta%process /= np) then ! check if the block is the leaf ! - if (pmeta%leaf) then + if (pmeta%leaf) then ! generate a tag for communication ! - itag = pmeta%process * nprocs + np + nprocs + 1 + itag = pmeta%process * nprocs + np + nprocs + 1 ! sends the block to the right process ! - if (nproc == pmeta%process) then + if (nproc == pmeta%process) then ! copy data to buffer ! - rbuf(1,:,:,:,:) = pmeta%data%u(:,:,:,:) - rbuf(2,:,:,:,:) = pmeta%data%q(:,:,:,:) + rbuf(1,:,:,:,:) = pmeta%data%u(:,:,:,:) + rbuf(2,:,:,:,:) = pmeta%data%q(:,:,:,:) ! send data ! - call send_real_array(size(rbuf), np, itag, rbuf, iret) + call send_real_array(size(rbuf), np, itag, rbuf, iret) ! remove data block from the current process ! - call remove_datablock(pmeta%data) + call remove_datablock(pmeta%data) ! send data block ! - end if ! nproc == pmeta%process + end if ! nproc == pmeta%process ! receive the block from another process ! - if (nproc == np) then + if (nproc == np) then ! allocate a new data block and link it with the current meta block ! - call append_datablock(pdata) - call link_blocks(pmeta, pdata) + call append_datablock(pdata) + call link_blocks(pmeta, pdata) ! receive the data ! - call receive_real_array(size(rbuf), pmeta%process, itag, rbuf, iret) + call receive_real_array(size(rbuf), pmeta%process, itag, rbuf, iret) ! coppy the buffer to data block ! - pmeta%data%u(:,:,:,:) = rbuf(1,:,:,:,:) - pmeta%data%q(:,:,:,:) = rbuf(2,:,:,:,:) + pmeta%data%u(:,:,:,:) = rbuf(1,:,:,:,:) + pmeta%data%q(:,:,:,:) = rbuf(2,:,:,:,:) - end if ! nproc == n + end if ! nproc == n - end if ! leaf + end if ! leaf ! set new processor number ! - pmeta%process = np + pmeta%process = np - end if ! pmeta%process /= np + end if ! pmeta%process /= np ! increase the number of blocks on the current process; if it exceeds the ! allowed number reset the counter and increase the processor number ! - if (pmeta%leaf) then + if (pmeta%leaf) then ! increase the number of leafs for the current process ! - nl = nl + 1 + nl = nl + 1 ! if the number of leafs for the current process exceeds the number of assigned ! blocks, reset the counter and increase the process number ! - if (nl >= lb(np)) then + if (nl >= lb(np)) then ! reset the leaf counter for the current process ! - nl = 0 + nl = 0 ! increase the process number ! - np = min(nprocs - 1, np + 1) + np = min(nprocs - 1, np + 1) - end if ! l >= lb(n) + end if ! l >= lb(n) - end if ! leaf + end if ! leaf + + end if ! pmeta%process < nprocs ! assign the pointer to the next meta block ! @@ -1502,7 +1514,7 @@ module mesh integer :: i, j, k, q, p integer :: il, iu, jl, ju, kl, ku integer :: ic, jc, kc, ip, jp, kp - real :: dul, dur, dux, duy, duz + real :: dul, dur, dux, duy, duz, du1, du2, du3, du4 ! local pointers ! @@ -1586,21 +1598,27 @@ module mesh #endif /* NDIMS == 3 */ #if NDIMS == 2 - u(p,ic,jc,kc) = pdata%u(p,i,j,k) - (dux + duy) - u(p,ip,jc,kc) = pdata%u(p,i,j,k) + (dux - duy) - u(p,ic,jp,kc) = pdata%u(p,i,j,k) + (duy - dux) - u(p,ip,jp,kc) = pdata%u(p,i,j,k) + (dux + duy) + du1 = dux + duy + du2 = dux - duy + u(p,ic,jc,kc) = pdata%u(p,i,j,k) - du1 + u(p,ip,jc,kc) = pdata%u(p,i,j,k) + du2 + u(p,ic,jp,kc) = pdata%u(p,i,j,k) - du2 + u(p,ip,jp,kc) = pdata%u(p,i,j,k) + du1 #endif /* NDIMS == 2 */ #if NDIMS == 3 - u(p,ic,jc,kc) = pdata%u(p,i,j,k) - dux - duy - duz - u(p,ip,jc,kc) = pdata%u(p,i,j,k) + dux - duy - duz - u(p,ic,jp,kc) = pdata%u(p,i,j,k) - dux + duy - duz - u(p,ip,jp,kc) = pdata%u(p,i,j,k) + dux + duy - duz - u(p,ic,jc,kp) = pdata%u(p,i,j,k) - dux - duy + duz - u(p,ip,jc,kp) = pdata%u(p,i,j,k) + dux - duy + duz - u(p,ic,jp,kp) = pdata%u(p,i,j,k) - dux + duy + duz - u(p,ip,jp,kp) = pdata%u(p,i,j,k) + dux + duy + duz + du1 = dux + duy + duz + du2 = dux - duy - duz + du3 = dux - duy + duz + du4 = dux + duy - duz + u(p,ic,jc,kc) = pdata%u(p,i,j,k) - du1 + u(p,ip,jc,kc) = pdata%u(p,i,j,k) + du2 + u(p,ic,jp,kc) = pdata%u(p,i,j,k) - du3 + u(p,ip,jp,kc) = pdata%u(p,i,j,k) + du4 + u(p,ic,jc,kp) = pdata%u(p,i,j,k) - du4 + u(p,ip,jc,kp) = pdata%u(p,i,j,k) + du3 + u(p,ic,jp,kp) = pdata%u(p,i,j,k) - du2 + u(p,ip,jp,kp) = pdata%u(p,i,j,k) + du1 #endif /* NDIMS == 3 */ end do end do diff --git a/src/mpitools.F90 b/src/mpitools.F90 index 0d137bc..93e74f5 100644 --- a/src/mpitools.F90 +++ b/src/mpitools.F90 @@ -141,7 +141,7 @@ module mpitools stop end if -! obtain the current process identificator +! obtain the current process identifier ! call mpi_comm_rank(mpi_comm_world, nproc , iret) diff --git a/src/problems.F90 b/src/problems.F90 index 68be81d..2cd4b59 100644 --- a/src/problems.F90 +++ b/src/problems.F90 @@ -131,6 +131,9 @@ module problems case("blast") setup_problem => setup_problem_blast + case("implosion") + setup_problem => setup_problem_implosion + case("kh", "kelvinhelmholtz", "kelvin-helmholtz") setup_problem => setup_problem_kelvin_helmholtz @@ -510,6 +513,228 @@ module problems ! !=============================================================================== ! +! subroutine SETUP_PROBLEM_IMPLOSION: +! ---------------------------------- +! +! Subroutine sets the initial conditions for the implosion problem. +! +! Arguments: +! +! pdata - pointer to the datablock structure of the currently initialized +! block; +! +!=============================================================================== +! + subroutine setup_problem_implosion(pdata) + +! include external procedures and variables +! + use blocks , only : block_data, ndims + use constants , only : d2r + use coordinates, only : im, jm, km + use coordinates, only : ax, ay, az, adx, ady, adz + use equations , only : prim2cons + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp + use parameters , only : get_parameter_real + +! local variables are not implicit by default +! + implicit none + +! input arguments +! + type(block_data), pointer, intent(inout) :: pdata + +! default parameter values +! + real(kind=8), save :: sline = 1.50d-01 + real(kind=8), save :: adens = 1.00d+00 + real(kind=8), save :: apres = 1.00d+00 + real(kind=8), save :: drat = 1.25d-01 + real(kind=8), save :: prat = 1.40d-01 + real(kind=8), save :: buni = 1.00d+00 + real(kind=8), save :: bgui = 0.00d+00 + real(kind=8), save :: angle = 0.00d+00 + +! local saved parameters +! + logical , save :: first = .true. + real(kind=8), save :: odens = 1.25d-01 + real(kind=8), save :: opres = 1.40d-01 + +! local variables +! + integer :: i, j, k + real(kind=8) :: rl, ru, dx, dy, dz, dxh, dyh, dzh, ds, dl, dr + real(kind=8) :: sn, cs + +! local arrays +! + real(kind=8), dimension(nv,im) :: q, u + real(kind=8), dimension(im) :: x, xl, xu + real(kind=8), dimension(jm) :: y, yl, yu + real(kind=8), dimension(km) :: z, zl, zu +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for the problem setup +! + call start_timer(imu) +#endif /* PROFILE */ + +! prepare problem constants during the first subroutine call +! + if (first) then + +! get problem parameters +! + call get_parameter_real("shock_line" , sline ) + call get_parameter_real("ambient_density" , adens ) + call get_parameter_real("ambient_pressure", apres ) + call get_parameter_real("density_ratio" , drat ) + call get_parameter_real("pressure_ratio" , prat ) + call get_parameter_real("buni" , buni ) + call get_parameter_real("bgui" , bgui ) + call get_parameter_real("angle" , angle ) + +! calculate parameters +! + odens = drat * adens + opres = prat * apres + +! reset the first execution flag +! + first = .false. + + end if ! first call + +! prepare block coordinates +! + x(1:im) = pdata%meta%xmin + ax(pdata%meta%level,1:im) + y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm) +#if NDIMS == 3 + z(1:km) = pdata%meta%zmin + az(pdata%meta%level,1:km) +#endif /* NDIMS == 3 */ + +! calculate mesh intervals and areas +! + dx = adx(pdata%meta%level) + dy = ady(pdata%meta%level) +#if NDIMS == 3 + dz = adz(pdata%meta%level) +#endif /* NDIMS == 3 */ + dxh = 0.5d+00 * dx + dyh = 0.5d+00 * dy +#if NDIMS == 3 + dzh = 0.5d+00 * dz +#endif /* NDIMS == 3 */ + +! calculate edge coordinates +! + xl(:) = abs(x(:)) - dxh + xu(:) = abs(x(:)) + dxh + yl(:) = abs(y(:)) - dyh + yu(:) = abs(y(:)) + dyh +#if NDIMS == 3 + zl(:) = abs(z(:)) - dzh + zu(:) = abs(z(:)) + dzh +#endif /* NDIMS == 3 */ + +! reset velocity components +! + q(ivx,:) = 0.0d+00 + q(ivy,:) = 0.0d+00 + q(ivz,:) = 0.0d+00 + +! if magnetic field is present, set it to be uniform with the desired strength +! and orientation +! + if (ibx > 0) then + +! calculate the orientation angles +! + sn = sin(d2r * angle) + cs = sqrt(1.0d+00 - sn * sn) + +! set magnetic field components +! + q(ibx,:) = buni * cs + q(iby,:) = buni * sn + q(ibz,:) = bgui + q(ibp,:) = 0.0d+00 + + end if + +! iterate over all positions +! + do k = 1, km + do j = 1, jm + do i = 1, im + +! calculate the distance from the origin +! +#if NDIMS == 3 + rl = xl(i) + yl(j) + zl(k) + ru = xu(i) + yu(j) + zu(k) +#else /* NDIMS == 3 */ + rl = xl(i) + yl(j) + ru = xu(i) + yu(j) +#endif /* NDIMS == 3 */ + +! initialize density and pressure +! + if (ru <= sline) then + q(idn,i) = odens + if (ipr > 0) q(ipr,i) = opres + else if (rl >= sline) then + q(idn,i) = adens + if (ipr > 0) q(ipr,i) = apres + else + ds = (sline - rl) / dx + if (ds <= 1.0d+00) then + dl = 5.0d-01 * ds**ndims + dr = 1.0d+00 - dl + else + ds = (ru - sline) / dx + dr = 5.0d-01 * ds**ndims + dl = 1.0d+00 - dr + end if + + q(idn,i) = adens * dl + odens * dr + if (ipr > 0) q(ipr,i) = apres * dl + opres * dr + end if + + end do ! i = 1, im + +! convert the primitive variables to conservative ones +! + call prim2cons(im, q(1:nv,1:im), u(1:nv,1:im)) + +! copy the conserved variables to the current block +! + pdata%u(1:nv,1:im,j,k) = u(1:nv,1:im) + +! copy the primitive variables to the current block +! + pdata%q(1:nv,1:im,j,k) = q(1:nv,1:im) + + end do ! j = 1, jm + end do ! k = 1, km + +#ifdef PROFILE +! stop accounting time for the problems setup +! + call stop_timer(imu) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine setup_problem_implosion +! +!=============================================================================== +! ! subroutine SETUP_PROBLEM_KELVIN_HELMHOLTZ: ! ----------------------------------------- ! diff --git a/src/random.F90 b/src/random.F90 index 125112f..cf95008 100644 --- a/src/random.F90 +++ b/src/random.F90 @@ -214,6 +214,11 @@ module random ! integer , intent(in) :: np integer(kind=4), dimension(0:np-1), intent(in) :: seed + +! local variables +! + integer :: i, l + real :: r ! !------------------------------------------------------------------------------- ! @@ -225,10 +230,34 @@ module random ! set the seeds only if the input array and seeds have the same sizes ! - if (np .eq. nseeds) then + if (np == nseeds) then seeds(0:lseed) = seed(0:lseed) + else + +! if the input array and seeds have different sizes, expand or shrink seeds +! + select case(gentype) + case('random') + l = min(lseed, np - 1) + seeds(0:l) = seed(0:l) + if (l < lseed) then + do i = l + 1, lseed + call random_number(r) + seeds(i) = 123456789 * r + end do + end if + case default + l = nseeds / 2 + do i = 0, l - 1 + seeds(i) = seed(0) + end do + do i = l, lseed + seeds(i) = seed(np-1) + end do + end select + end if #ifdef PROFILE