gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into reconnection

Browse Source
This commit is contained in:
Grzegorz Kowal 2021-12-15 09:47:11 -03:00
commit 0cae0eb66e
19 changed files with 1360 additions and 1001 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
CMakeLists.txt
View File

@ -64,6 +64,16 @@ if(ENABLE_MPI)
endif()
endif()
option(ENABLE_OPENMP "Enable OpenMP support" ON)
if(ENABLE_OPENMP)
find_package(OpenMP COMPONENTS Fortran)
if(OpenMP_Fortran_FOUND)
include_directories(${OpenMP_Fortran_INCLUDE_DIRS})
target_compile_options(amun.x PRIVATE "${OpenMP_Fortran_FLAGS}")
target_link_libraries(amun.x ${OpenMP_Fortran_LIBRARIES})
endif()
endif()
option(ENABLE_SIGNALS "Enables signal handler support" ON)
if(ENABLE_SIGNALS)
target_compile_definitions(amun.x PRIVATE SIGNALS)

2
README.md
View File

@ -29,7 +29,7 @@ following features are already implemented:
* support for Zstandard, LZ4, and LZMA compressions in the XML+binary format,
* support for Deflate, SZIP, Zstandard, and ZFP compressions in the HDF5 format,
* easy and consistend Python interface to read snapshots in both formats,
* MPI parallelization,
* MPI/OpenMP parallelization,
* completely written in Fortran 2008,
* simple Makefile or CMake for building the code executable,
* minimum requirements, only Fortran compiler and Python are required to

12
build/hosts/default
View File

@ -27,6 +27,9 @@ FFLAGS = -Og -g -DDEBUG -fcheck=all
else
FFLAGS = -O2 -march=native -pipe
endif
ifeq ($(OPENMP),Y)
FFLAGS += -fopenmp
endif
LDFLAGS = $(FFLAGS)
ifeq ($(STATIC),Y)
LDFLAGS += -static
@ -56,6 +59,9 @@ FFLAGS = -O0 -g -Mbounds -Mchkptr -Mcache_align -Mnovintr -Mchkstk -DDEBUG
else
FFLAGS = -fast -O3
endif
ifeq ($(OPENMP),Y)
FFLAGS += -mp
endif
LDFLAGS = $(FFLAGS)
ifeq ($(STATIC),Y)
LDFLAGS += -Bstatic
@ -85,6 +91,9 @@ FFLAGS = -O -g -DDEBUG
else
FFLAGS = -fast
endif
ifeq ($(OPENMP),Y)
FFLAGS += -mp
endif
LDFLAGS = $(FFLAGS)
ifeq ($(STATIC),Y)
LDFLAGS += -Bstatic
@ -114,6 +123,9 @@ FFLAGS = -O -g -DDEBUG
else
FFLAGS = -O2 -xHost
endif
ifeq ($(OPENMP),Y)
FFLAGS += -openmp
endif
LDFLAGS = $(FFLAGS)
ifeq ($(STATIC),Y)
LDFLAGS += -static

2
build/make.default
View File

@ -15,8 +15,10 @@ SIGNALS = N
# parallelization flags:
#
# MPI - enable or disable MPI parallelization
# OPENMP - enable or disable OpenMP parallelization
#
MPI = N
OPENMP = N
# output data format:
#

26
sources/amun.F90
View File

@ -35,15 +35,18 @@ program amun
use helpers , only : print_message
use io , only : initialize_io, finalize_io
use mpitools , only : initialize_mpitools, finalize_mpitools
use mpitools , only : verbose => master, nproc, check_status
#ifdef MPI
use mpitools , only : reduce_sum
use mpitools , only : reduce_sum, nprocs
#endif /* MPI */
use mpitools , only : verbose => master, nprocs, nproc, check_status
use parameters, only : read_parameters, finalize_parameters
use random , only : initialize_random, finalize_random
use system , only : initialize_system, finalize_system
use system , only : print_system_info, prepare_system, evolve_system
use system , only : quit, nsteps, nwork
use system , only : nsteps, nwork
#ifdef SIGNALS
use system , only : quit
#endif /* SIGNALS */
use timers , only : initialize_timers, finalize_timers
use timers , only : start_timer, stop_timer, set_timer, get_timer
use timers , only : get_timer_total, timer_enabled, timer_description
@ -62,6 +65,11 @@ program amun
integer :: iin, iev, itm
integer :: ed, eh, em, es
! OpenMP
!
integer :: nthreads = 1
!$ integer :: omp_get_num_threads
! the format string
!
character(len=80) :: sfmt
@ -128,8 +136,16 @@ program amun
call print_welcome(verbose)
#ifdef MPI
call print_section(verbose, "Parallelization")
call print_parameter(verbose, "MPI processes", nprocs)
call print_parameter(verbose, "MPI processes" , nprocs)
#else /* MPI */
!$ call print_section(verbose, "Parallelization")
#endif /* MPI */
!$omp parallel
!$omp master
!$ nthreads = omp_get_num_threads()
!$ call print_parameter(verbose, "OpenMP threads", nthreads)
!$omp end master
!$omp end parallel
! read parameters
!
@ -155,7 +171,7 @@ program amun
call print_message(loc, "Could not initialize module SYSTEM!")
if (check_status(status /= 0)) go to 3000
call initialize_workspace(nwork, status)
call initialize_workspace(nwork, nthreads, status)
if (status /= 0) &
call print_message(loc, "Could not initialize module WORKSPACE!")
if (check_status(status /= 0)) go to 2000

80
sources/blocks.F90
View File

@ -303,13 +303,18 @@ module blocks
type(block_data), pointer, save :: list_data, last_data
type(block_leaf), pointer, save :: list_leaf
! THE VECTOR OF DATA BLOCK POINTERS:
! =================================
!
type(pointer_data), dimension(:), allocatable :: data_blocks
! all variables and subroutines are private by default
!
private
! declare public pointers, structures, and variables
!
public :: pointer_meta, pointer_info
public :: pointer_meta, pointer_data, pointer_info
public :: block_meta, block_data, block_info, block_leaf
public :: list_meta, list_data, list_leaf
public :: ndims, nsides, nchildren, nregs
@ -331,7 +336,8 @@ module blocks
public :: metablock_set_configuration, metablock_set_refinement
public :: metablock_set_position, metablock_set_coordinates
public :: metablock_set_bounds, metablock_set_leaf, metablock_unset_leaf
public :: build_leaf_list
public :: build_leaf_list, build_datablock_list
public :: data_blocks
#ifdef DEBUG
public :: check_neighbors
#endif /* DEBUG */
@ -426,16 +432,29 @@ module blocks
!
subroutine finalize_blocks(status)
use iso_fortran_env, only : error_unit
implicit none
integer, intent(out) :: status
type(block_meta), pointer :: pmeta
character(len=*), parameter :: loc = 'BLOCKS::finalize_blocks()'
!-------------------------------------------------------------------------------
!
status = 0
! deallocate the vector of data blocks
!
if (allocated(data_blocks)) then
deallocate(data_blocks, stat=status)
if (status /= 0) &
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not deallocate the vector of data blocks!"
end if
! remove all metablocks
!
pmeta => last_meta
@ -3545,6 +3564,63 @@ module blocks
!
!===============================================================================
!
! subroutine BUILD_DATABLOCK_LIST:
! -------------------------------
!
! Subroutine builds the list of data blocks data_blocks.
!
! Arguments:
!
! status - the subroutine call status: 0 for success, otherwise failure;
!
!===============================================================================
!
subroutine build_datablock_list(status)
use iso_fortran_env, only : error_unit
implicit none
integer, intent(out) :: status
type(block_data), pointer :: pdata
integer :: l
character(len=*), parameter :: loc = 'BLOCKS::build_datablock_list()'
!-------------------------------------------------------------------------------
!
status = 0
if (allocated(data_blocks)) then
deallocate(data_blocks, stat=status)
if (status /= 0) &
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not deallocate the vector of data blocks!"
end if
allocate(data_blocks(dblocks), stat=status)
if (status == 0) then
l = 0
pdata => list_data
do while (associated(pdata))
l = l + 1
data_blocks(l)%ptr => pdata
pdata => pdata%next
end do
else
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not allocate the vector of data blocks!"
end if
!-------------------------------------------------------------------------------
!
end subroutine build_datablock_list
!
!===============================================================================
!
! subroutine SET_NEIGHBORS_REFINE:
! -------------------------------
!

83
sources/boundaries.F90
View File

@ -1546,7 +1546,9 @@ module boundaries
use blocks , only : block_meta, block_data, block_leaf
use blocks , only : list_leaf
use blocks , only : block_info, pointer_info
#ifdef MPI
use coordinates, only : nh => ncells_half, ng => nghosts
#endif /* MPI */
use coordinates, only : faces_gc, faces_dc
use equations , only : nv
#ifdef MPI
@ -1888,7 +1890,9 @@ module boundaries
use blocks , only : block_meta, block_data, block_leaf
use blocks , only : list_leaf
use blocks , only : block_info, pointer_info
#ifdef MPI
use coordinates, only : nh => ncells_half, ng => nghosts
#endif /* MPI */
use coordinates, only : faces_gr
use equations , only : nv
#ifdef MPI
@ -2638,7 +2642,9 @@ module boundaries
use blocks , only : block_meta, block_data, block_leaf
use blocks , only : list_leaf
use blocks , only : block_info, pointer_info
#ifdef MPI
use coordinates, only : nh => ncells_half, ng => nghosts
#endif /* MPI */
use coordinates, only : edges_gc, edges_dc
use equations , only : nv
#ifdef MPI
@ -3060,7 +3066,9 @@ module boundaries
use blocks , only : block_meta, block_data, block_leaf
use blocks , only : list_leaf
use blocks , only : block_info, pointer_info
#ifdef MPI
use coordinates, only : nh => ncells_half, ng => nghosts
#endif /* MPI */
use coordinates, only : edges_gr
use equations , only : nv
#ifdef MPI
@ -3973,7 +3981,9 @@ module boundaries
use blocks , only : block_meta, block_data, block_leaf
use blocks , only : list_leaf
use blocks , only : block_info, pointer_info
#ifdef MPI
use coordinates, only : ng => nghosts
#endif /* MPI */
use coordinates, only : corners_gc, corners_dc
use equations , only : nv
#ifdef MPI
@ -4349,7 +4359,9 @@ module boundaries
use blocks , only : block_meta, block_data, block_leaf
use blocks , only : list_leaf
use blocks , only : block_info, pointer_info
#ifdef MPI
use coordinates, only : ng => nghosts
#endif /* MPI */
use coordinates, only : corners_gr
use equations , only : nv
#ifdef MPI
@ -4706,7 +4718,9 @@ module boundaries
use blocks , only : block_meta, block_data, block_leaf
use blocks , only : list_leaf
use blocks , only : block_info, pointer_info
#ifdef MPI
use coordinates, only : ng => nghosts
#endif /* MPI */
use coordinates, only : corners_gp
use equations , only : nv
#ifdef MPI
@ -6468,93 +6482,64 @@ module boundaries
!
subroutine update_ghost_cells()
! include external variables
!
use blocks , only : block_data, list_data
use blocks , only : block_data, data_blocks, get_dblocks
use coordinates, only : nn => bcells, nb, ne, nbl, neu
use equations , only : prim2cons
! local variables are not implicit by default
!
implicit none
! local variables
!
integer :: i, j, k = 1
integer :: i, j, k, l, n
! local pointers
!
type(block_data), pointer :: pdata
!
!-------------------------------------------------------------------------------
!
! assign the pointer to the first block on the list
!
pdata => list_data
n = get_dblocks()
! scan all data blocks until the last is reached
!
do while(associated(pdata))
!$omp parallel do default(shared) private(pdata,i,j,k)
do l = 1, n
pdata => data_blocks(l)%ptr
! update the X and Y boundary ghost cells
!
#if NDIMS == 3
do k = 1, nn
#else /* NDIMS == 3 */
do k = 1, 1
#endif /* NDIMS == 3 */
! update lower layers of the Y boundary
!
do j = 1, nbl
call prim2cons(pdata%q(:,:,j,k), pdata%u(:,:,j,k), .true.)
end do ! j = 1, nbl
end do
! update upper layers of the Y boundary
!
do j = neu, nn
call prim2cons(pdata%q(:,:,j,k), pdata%u(:,:,j,k), .true.)
end do ! j = neu, nn
end do
! update remaining left layers of the X boundary
!
do i = 1, nbl
call prim2cons(pdata%q(:,i,nb:ne,k), pdata%u(:,i,nb:ne,k), .true.)
end do ! i = 1, nbl
end do
! update remaining right layers of the X boundary
!
do i = neu, nn
call prim2cons(pdata%q(:,i,nb:ne,k), pdata%u(:,i,nb:ne,k), .true.)
end do ! i = neu, nn
end do
end do
#if NDIMS == 3
end do ! k = 1, nn
#endif /* NDIMS == 3 */
#if NDIMS == 3
! update the Z boundary ghost cells
!
do j = nb, ne
! update the remaining front layers of the Z boundary
!
do k = 1, nbl
call prim2cons(pdata%q(:,nb:ne,j,k), pdata%u(:,nb:ne,j,k), .true.)
end do ! k = 1, nbl
end do
! update the remaining back layers of the Z boundary
!
do k = neu, nn
call prim2cons(pdata%q(:,nb:ne,j,k), pdata%u(:,nb:ne,j,k), .true.)
end do ! k = neu, nn
end do
end do ! j = nb, ne
end do
#endif /* NDIMS == 3 */
! assign the pointer to the next block on the list
!
pdata => pdata%next
end do ! data blocks
end do
!$omp end parallel do
!-------------------------------------------------------------------------------
!

8
sources/compression.F90
View File

@ -53,8 +53,8 @@ module compression
integer(kind=c_int) , value :: level
end function zstd_compress
integer function zstd_iserror(code) bind(C, name="ZSTD_isError")
use iso_c_binding, only: c_size_t
integer(kind=c_int) function zstd_iserror(code) bind(C, name="ZSTD_isError")
use iso_c_binding, only: c_int, c_size_t
implicit none
integer(kind=c_size_t), value :: code
end function zstd_iserror
@ -81,8 +81,8 @@ module compression
type(c_ptr) , value :: src, dst, prefsPtr
end function lz4_compress
integer function lz4_iserror(code) bind(C, name="LZ4F_isError")
use iso_c_binding, only: c_size_t
integer(kind=c_int) function lz4_iserror(code) bind(C, name="LZ4F_isError")
use iso_c_binding, only: c_int, c_size_t
implicit none
integer(kind=c_size_t), value :: code
end function lz4_iserror

1604
sources/evolution.F90
View File

File diff suppressed because it is too large Load Diff

19
sources/forcing.F90
View File

@ -1398,10 +1398,15 @@ module forcing
real(kind=8), dimension(:,:,:,:), pointer, save :: acc
real(kind=8), dimension(:,:,:) , pointer, save :: den
integer :: nt = 0
!$ integer :: omp_get_thread_num
!$omp threadprivate(first, nt, acc, den)
character(len=*), parameter :: loc = 'FORCING:inject_fmodes_block()'
!-------------------------------------------------------------------------------
!
!$ nt = omp_get_thread_num()
if (first) then
i = 3 * nn**NDIMS
j = 4 * nn**NDIMS
@ -1413,11 +1418,11 @@ module forcing
end if
#if NDIMS == 3
acc(1:3,1:nn,1:nn,1:nn) => work( 1:i)
den(1:nn,1:nn,1:nn) => work(i+1:j)
acc(1:3,1:nn,1:nn,1:nn) => work( 1:i,nt)
den(1:nn,1:nn,1:nn) => work(i+1:j,nt)
#else /* NDIMS == 3 */
acc(1:2,1:nn,1:nn,1: 1) => work(1:i)
den(1:nn,1:nn,1: 1) => work(i+1:j)
acc(1:2,1:nn,1:nn,1: 1) => work( 1:i,nt)
den(1:nn,1:nn,1: 1) => work(i+1:j,nt)
#endif /* NDIMS == 3 */
first = .false.
@ -1432,11 +1437,11 @@ module forcing
#endif /* NDIMS == 3 */
dvol = advol(pdata%meta%level)
if (work_in_use) &
if (work_in_use(nt)) &
call print_message(loc, "Workspace is being used right now! " // &
"Corruptions can occur!")
work_in_use = .true.
work_in_use(nt) = .true.
acc(1:NDIMS,:,:,:) = 0.0d+00
@ -1524,7 +1529,7 @@ module forcing
pdata%u(ien,:,:,:) = pdata%u(ien,:,:,:) + den(:,:,:)
end if
work_in_use = .false.
work_in_use(nt) = .false.
100 continue

4
sources/hash.F90
View File

@ -206,6 +206,7 @@ module hash
!-------------------------------------------------------------------------------
!
hash = 0
select case(hash_id)
case(hash_xxh64)
#ifndef XXHASH
@ -215,8 +216,7 @@ module hash
case(hash_xxh3)
hash = xxh3_lib(buffer, length)
#endif /* XXHASH */
case(hash_none)
hash = 0
case default
end select
return

84
sources/io.F90
View File

@ -107,9 +107,11 @@ module io
!
logical , save :: precise_snapshots = .false.
#ifdef HDF5
! a flag to determine if XDMF files should be generated
!
logical , save :: with_xdmf = .false.
logical , save :: with_xdmf = .false.
#endif /* HDF5 */
! the compression format and level of the XML+binary files
!
@ -186,7 +188,9 @@ module io
use compression, only : set_compression, get_compression
use hash , only : hash_info
use helpers , only : print_message
#ifdef HDF5
use mpitools , only : nproc
#endif /* HDF5 */
use parameters , only : get_parameter
implicit none
@ -343,6 +347,7 @@ module io
call print_message(loc, &
"Could not check if the filter is available!")
end if
with_xdmf = .false.
case("zstd", "zstandard")
call h5zfilter_avail_f(H5Z_ZSTANDARD, test, status)
if (status == 0) then
@ -357,6 +362,7 @@ module io
call print_message(loc, &
"Could not check if the filter is available!")
end if
with_xdmf = .false.
case("zfp")
call h5zfilter_avail_f(H5Z_ZFP, test, status)
if (status == 0) then
@ -382,6 +388,7 @@ module io
call print_message(loc, &
"Could not check if the filter is available!")
end if
with_xdmf = .false.
case default
end select
@ -1080,7 +1087,10 @@ module io
logical :: test
character(len=255) :: dname, fname, line, sname, svalue
integer :: il, iu, nl, nx, nm, nd, nv, i, j, k, l, n, p, nu, nr
integer :: il, iu, nl, nx, nm, nd, nv, i, j, l, n, p, nu, nr
#if NDIMS == 3
integer :: k
#endif /* NDIMS == 3 */
integer(kind=4) :: lndims, lnprocs, lnproc, lmblocks, lnleafs, llast_id
integer(kind=4) :: ldblocks, lncells, lnghosts, lnseeds, lnmodes
real(kind=8) :: deinj
@ -5436,29 +5446,21 @@ module io
!
subroutine write_snapshot_xdmf()
! import procedures and variables from other modules
!
use blocks , only : block_data, list_data
use blocks , only : get_dblocks
use equations , only : nv, pvars
use mpitools , only : nproc
use coordinates , only : ni => ncells, ng => nghosts
use coordinates , only : adx, ady
use blocks , only : block_data, list_data
use blocks , only : get_dblocks
use equations , only : nv, pvars
use mpitools , only : nproc
use coordinates, only : ni => ncells, ng => nghosts
use coordinates, only : adx, ady
#if NDIMS == 3
use coordinates , only : adz
use coordinates, only : adz
#endif /* NDIMS == 3 */
use evolution , only : time
use evolution , only : time
! local variables are not implicit by default
!
implicit none
! local pointers
!
type(block_data), pointer :: pdata
! local variables
!
character(len=64) :: fname, hname
character(len=128) :: stmp, ttmp, sdim, bdim, pdim
integer(kind=4) :: l, p
@ -5467,23 +5469,15 @@ module io
integer(kind=4) :: kp
#endif /* NDIMS == 3 */
! local arrays
!
integer, dimension(12) :: slab
! local parameters
!
integer, parameter :: xdmf = 101
!
!-------------------------------------------------------------------------------
!
! prepare the XDMF and HDF5 file names
write(fname, "('p',i6.6,'_',i5.5,'.xdmf')") isnap, nproc
write(hname, "('p',i6.6,'_',i5.5,'.h5' )") isnap, nproc
! open the XDMF file
!
open (unit = xdmf, file = fname, status = 'replace')
! write the header
@ -5522,20 +5516,24 @@ module io
#endif /* NDIMS == 3 */
bdim = trim(adjustl(stmp)) // ' ' // trim(adjustl(ttmp))
#if NDIMS == 3
write(stmp, "(1i8)") ni
#else /* NDIMS == 3 */
write(stmp, "(1i8)") 1
#endif /* NDIMS == 3 */
write(ttmp, "(1i8)") ni
stmp = trim(adjustl(stmp)) // ' ' // trim(adjustl(ttmp))
write(ttmp, "(1i8)") ni
stmp = trim(adjustl(stmp)) // ' ' // trim(adjustl(ttmp))
write(ttmp, "(1i8)") get_dblocks()
sdim = trim(adjustl(stmp)) // ' ' // trim(adjustl(ttmp))
sdim = trim(adjustl(ttmp)) // ' ' // trim(adjustl(stmp))
! prepare slab indices
!
#if NDIMS == 3
slab(:) = (/ ng, ng, ng, -1, 1, 1, 1, 1, ni, ni, ni, 1 /)
slab(:) = (/ -1, ng, ng, ng, 1, 1, 1, 1, 1, ni, ni, ni /)
#else /* NDIMS == 3 */
slab(:) = (/ 0, ng, ng, -1, 1, 1, 1, 1, 1, ni, ni, 1 /)
slab(:) = (/ -1, 0, ng, ng, 1, 1, 1, 1, 1, 1, ni, ni /)
#endif /* NDIMS == 3 */
! iterate over all data blocks
@ -5588,7 +5586,7 @@ module io
! convert slab dimensions to string
!
slab(4) = l
slab(1) = l
write(pdim, "(1i8)") slab(1)
do p = 2, 12
write(ttmp, "(1i8)") slab(p)
@ -5633,8 +5631,6 @@ module io
write(xdmf, "(a)") ' </Domain>'
write(xdmf, "(a)") '</Xdmf>'
! close the XDMF file
!
close(xdmf)
!-------------------------------------------------------------------------------
@ -5657,31 +5653,19 @@ module io
!
subroutine write_snapshot_xdmf_master()
! import procedures and variables from other modules
!
use mpitools , only : npmax
use mpitools, only : npmax
! local variables are not implicit by default
!
implicit none
! local variables
!
character(len=64) :: fname, pname
integer(kind=4) :: np
character(len=64) :: fname, pname
integer(kind=4) :: np
! local parameters
!
integer, parameter :: xdmf = 102
!
!-------------------------------------------------------------------------------
!
! prepare the XDMF and HDF5 file names
write(fname, "('p',i6.6,'.xdmf')") isnap
! open the XDMF file
!
open (unit = xdmf, file = fname, status = 'replace')
! write the header
@ -5707,8 +5691,6 @@ module io
write(xdmf, "(a)") ' </Domain>'
write(xdmf, "(a)") '</Xdmf>'
! close the XDMF file
!
close(xdmf)
!-------------------------------------------------------------------------------

266
sources/mesh.F90
View File

@ -202,7 +202,7 @@ module mesh
use blocks , only : link_blocks, unlink_blocks, refine_block
use blocks , only : get_mblocks, get_nleafs
use blocks , only : set_neighbors_refine
use blocks , only : build_leaf_list
use blocks , only : build_leaf_list, build_datablock_list
#ifdef DEBUG
use blocks , only : check_neighbors
#endif /* DEBUG */
@ -483,6 +483,10 @@ module mesh
!
100 continue
! update the vector of data blocks
!
call build_datablock_list(status)
!-------------------------------------------------------------------------------
!
end subroutine generate_mesh
@ -505,26 +509,21 @@ module mesh
!
subroutine update_mesh(status)
! import external procedures and variables
!
use blocks, only : build_leaf_list
use blocks, only : block_data
use blocks, only : build_leaf_list, build_datablock_list
#ifdef DEBUG
use blocks, only : check_neighbors
#endif /* DEBUG */
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer, intent(out) :: status
!-------------------------------------------------------------------------------
!
! check the refinement criterion of all data blocks at the current process
!
call check_data_block_refinement(status)
call check_block_refinement(status)
if (status /= 0) go to 100
! update neighbor refinement flags, if they need to be refined as well
@ -562,10 +561,12 @@ module mesh
call check_neighbors()
#endif /* DEBUG */
! error entry point
!
100 continue
! update the vector of data blocks
!
call build_datablock_list(status)
!-------------------------------------------------------------------------------
!
end subroutine update_mesh
@ -591,7 +592,7 @@ module mesh
use blocks , only : append_datablock, remove_datablock, link_blocks
use helpers , only : print_message
use mpitools, only : check_status
use mpitools, only : nprocs, npmax, nproc, nodes, lprocs
use mpitools, only : npmax, nproc, nodes, lprocs
use mpitools, only : send_array, receive_array
#endif /* MPI */
@ -637,58 +638,66 @@ module mesh
pmeta => list_meta
do while (associated(pmeta))
if (pmeta%process < nprocs) then
! exchange data blocks between processes
!
if (pmeta%process /= np) then
if (pmeta%leaf) then
changed = .true.
if (pmeta%process /= np) then
if (pmeta%leaf) then
changed = .true.
tag1 = pmeta%id
tag2 = pmeta%id + get_last_id()
tag1 = pmeta%id
tag2 = pmeta%id + get_last_id()
if (nproc == pmeta%process) then
if (nproc == pmeta%process) then
call send_array(np, tag1, pmeta%data%uu)
call send_array(np, tag2, pmeta%data%q)
call remove_datablock(pmeta%data, status)
if (status /= 0) then
call print_message(loc, "Could not remove data block!")
go to 1000
end if
call send_array(np, tag1, pmeta%data%uu)
call send_array(np, tag2, pmeta%data%q)
call remove_datablock(pmeta%data, status)
if (status /= 0) then
call print_message(loc, "Could not remove data block!")
go to 1000
end if
if (nproc == np) then
call append_datablock(pdata, status)
if (status /= 0) then
call print_message(loc, "Could not append new data block!")
go to 1000
end if
call link_blocks(pmeta, pdata)
call receive_array(pmeta%process, tag1, pmeta%data%uu)
call receive_array(pmeta%process, tag2, pmeta%data%q)
end if
end if
pmeta%process = np
if (nproc == np) then
call append_datablock(pdata, status)
if (status /= 0) then
call print_message(loc, "Could not append new data block!")
go to 1000
end if
call link_blocks(pmeta, pdata)
call receive_array(pmeta%process, tag1, pmeta%data%uu)
call receive_array(pmeta%process, tag2, pmeta%data%q)
end if
end if
pmeta%process = np
end if
! determine the new block distribution
!
if (pmeta%leaf) then
nl = nl + 1
if (pmeta%leaf) then
nl = nl + 1
if (nl >= lb(nr,nc)) then
if (nl >= lb(nr,nc)) then
np = min(npmax, np + 1)
nl = 0
nr = nr + 1
if (nr > lprocs) then
nr = 1
nc = nc + 1
end if
flag = nc <= nodes
if (flag) flag = lb(nr,nc) == 0
do while(flag)
np = min(npmax, np + 1)
nl = 0
nr = nr + 1
if (nr > lprocs) then
nr = 1
@ -696,20 +705,10 @@ module mesh
end if
flag = nc <= nodes
if (flag) flag = lb(nr,nc) == 0
do while(flag)
np = min(npmax, np + 1)
nr = nr + 1
if (nr > lprocs) then
nr = 1
nc = nc + 1
end if
flag = nc <= nodes
if (flag) flag = lb(nr,nc) == 0
end do
end if
end do
end if
end if
pmeta => pmeta%next
@ -1193,8 +1192,8 @@ module mesh
!
!===============================================================================
!
! subroutine CHECK_DATA_BLOCK_REFINEMENT:
! --------------------------------------
! subroutine CHECK_BLOCK_REFINEMENT:
! ---------------------------------
!
! Subroutine scans over all data blocks, gets and corrects their refinement
! flags. If the MPI is used, the refinement flags are syncronized among all
@ -1206,18 +1205,16 @@ module mesh
!
!===============================================================================
!
subroutine check_data_block_refinement(status)
subroutine check_block_refinement(status)
use blocks , only : block_meta, block_data, list_meta, list_data
use blocks , only : block_meta, list_meta
use blocks , only : block_data, data_blocks, get_dblocks
#ifdef MPI
use blocks , only : get_nleafs
#endif /* MPI */
use coordinates, only : minlev, maxlev
use helpers , only : print_message
#ifdef MPI
#ifdef DEBUG
use mpitools , only : nproc
#endif /* DEBUG */
use mpitools , only : reduce_sum
#endif /* MPI */
use refinement , only : check_refinement_criterion
@ -1229,21 +1226,22 @@ module mesh
type(block_meta), pointer :: pmeta
type(block_data), pointer :: pdata
integer :: l, n
#ifdef MPI
integer(kind=4) :: nl, l
integer(kind=4) :: nl
integer(kind=4), dimension(:), allocatable :: ibuf
#endif /* MPI */
character(len=*), parameter :: loc = 'MESH::check_data_block_refinement()'
character(len=*), parameter :: loc = 'MESH::check_block_refinement()'
#endif /* MPI */
!-------------------------------------------------------------------------------
!
status = 0
n = get_dblocks()
! 1) reset the refinement flag for all meta blocks
!
! iterate over all meta blocks
!
pmeta => list_meta
do while (associated(pmeta))
@ -1253,149 +1251,74 @@ module mesh
! 2) determine the refinement of data block from the current process
!
! iterate over all data blocks
!
pdata => list_data
do while (associated(pdata))
! assign pmeta to the meta block associated with pdata
!
!$omp parallel do default(shared) private(pdata,pmeta)
do l = 1, n
pdata => data_blocks(l)%ptr
pmeta => pdata%meta
#ifdef DEBUG
! check if pmeta is associated
!
if (associated(pmeta)) then
! check if pmeta is a leaf
!
if (pmeta%leaf) then
#endif /* DEBUG */
! check the refinement criterion for the current data block
!
pmeta%refine = check_refinement_criterion(pdata)
pmeta%refine = check_refinement_criterion(pdata)
! correct the refinement flag for the minimum and maximum levels
!
if (pmeta%level < minlev) pmeta%refine = 1
if (pmeta%level == minlev) pmeta%refine = max(0, pmeta%refine)
if (pmeta%level == maxlev) pmeta%refine = min(0, pmeta%refine)
if (pmeta%level > maxlev) pmeta%refine = -1
if (pmeta%level < minlev) pmeta%refine = 1
if (pmeta%level == minlev) pmeta%refine = max(0, pmeta%refine)
if (pmeta%level == maxlev) pmeta%refine = min(0, pmeta%refine)
if (pmeta%level > maxlev) pmeta%refine = -1
#ifdef DEBUG
else ! pmeta is a leaf
call print_message(loc, "Associated meta block is not a leaf!")
end if ! pmeta is a leaf
else ! pmeta associated
call print_message(loc, &
"No meta block associated with the current data block!")
end if ! pmeta associated
#endif /* DEBUG */
pdata => pdata%next
end do ! iterate over data blocks
end do
!$omp end parallel do
#ifdef MPI
! 3) synchronize the refinement flags between all processes
!
! get the number of leafs
!
nl = get_nleafs()
! allocate a buffer for the refinement flags
!
allocate(ibuf(nl), stat = status)
! check if the buffer was allocated successfully
!
if (status == 0) then
! reset the buffer
!
ibuf(1:nl) = 0
! reset the leaf block counter
!
l = 0
! iterate over all meta blocks
!
pmeta => list_meta
do while (associated(pmeta))
! process only leafs
!
if (pmeta%leaf) then
! increase the leaf block counter
!
l = l + 1
! store the refinement flag in the buffer
!
ibuf(l) = pmeta%refine
end if ! pmeta is a leaf
end if
pmeta => pmeta%next
end do ! iterate over meta blocks
end do
! update refinement flags across all processes
!
call reduce_sum(ibuf(1:nl))
! reset the leaf block counter
!
l = 0
! iterate over all meta blocks
!
pmeta => list_meta
do while (associated(pmeta))
! process only leafs
!
if (pmeta%leaf) then
! increase the leaf block counter
!
l = l + 1
#ifdef DEBUG
! check if the MPI update process does not change the local refinement flags
!
if (pmeta%process == nproc .and. pmeta%refine /= ibuf(l)) &
call print_message(loc, &
"Refinement flag does not match after MPI reduction!")
#endif /* DEBUG */
! restore the refinement flags
!
pmeta%refine = ibuf(l)
end if ! pmeta is a leaf
end if
pmeta => pmeta%next
end do ! iterate over meta blocks
end do
! deallocate the refinement flag buffer
!
deallocate(ibuf, stat = status)
if (status /= 0) &
call print_message(loc, &
call print_message(loc, &
"Refinement flag buffer could not be deallocated!")
else ! buffer couldn't be allocated
else
call print_message(loc, "Refinement flag buffer could not be allocated!")
end if ! buffer couldn't be allocated
end if
#endif /* MPI */
!-------------------------------------------------------------------------------
!
end subroutine check_data_block_refinement
end subroutine check_block_refinement
!
!===============================================================================
!
@ -1501,7 +1424,7 @@ module mesh
use coordinates, only : toplev
use helpers , only : print_message
#ifdef MPI
use mpitools , only : nprocs, nproc
use mpitools , only : nproc
use mpitools , only : send_array, receive_array
#endif /* MPI */
@ -1521,7 +1444,9 @@ module mesh
logical, dimension(nchildren) :: flag
#ifdef MPI
character(len=*), parameter :: loc = 'MESH::prepare_sibling_derefinement()'
#endif /* MPI */
!-------------------------------------------------------------------------------
!
@ -1619,9 +1544,9 @@ module mesh
pmeta => pmeta%next
end do
#endif /* MPI */
100 continue
#endif /* MPI */
!-------------------------------------------------------------------------------
!
@ -1935,10 +1860,15 @@ module mesh
real(kind=8), dimension(:,:,:), pointer, save :: tmp
integer :: nt = 0
!$ integer :: omp_get_thread_num
!$omp threadprivate(first, nt, tmp)
character(len=*), parameter :: loc = 'MESH::prolong_block()'
!-------------------------------------------------------------------------------
!
!$ nt = omp_get_thread_num()
status = 0
if (first) then
@ -1950,16 +1880,16 @@ module mesh
go to 100
end if
tmp(1:pm(1),1:pm(2),1:pm(3)) => work(1:n)
tmp(1:pm(1),1:pm(2),1:pm(3)) => work(1:n,nt)
first = .false.
end if
if (work_in_use) &
if (work_in_use(nt)) &
call print_message(loc, "Workspace is being used right now! " // &
"Corruptions can occur!")
work_in_use = .true.
work_in_use(nt) = .true.
#if NDIMS == 2
k = 1
@ -2063,7 +1993,7 @@ module mesh
end do ! nchildren
end do ! n = 1, nv
work_in_use = .false.
work_in_use(nt) = .false.
100 continue

10
sources/mpitools.F90
View File

@ -60,11 +60,11 @@ module mpitools
module procedure exchange_arrays_diff
module procedure exchange_arrays_same
end interface
#endif /* MPI */
! timer indices
!
integer, save :: imi, imc
#endif /* MPI */
! MPI global variables
!
@ -131,9 +131,9 @@ module mpitools
integer :: ierror
integer(kind=4), dimension(:), allocatable :: procs
#endif /* MPI */
character(len=*), parameter :: loc = 'MPITOOLS::initialize_mpitools()'
#endif /* MPI */
!-------------------------------------------------------------------------------
!
@ -275,9 +275,9 @@ module mpitools
#ifdef MPI
integer :: ierror
#endif /* MPI */
character(len=*), parameter :: loc = 'MPITOOLS::finalize_mpitools()'
#endif /* MPI */
!-------------------------------------------------------------------------------
!
@ -332,7 +332,9 @@ module mpitools
!
!-------------------------------------------------------------------------------
!
#ifdef MPI
call start_timer(imc)
#endif /* MPI */
check_status = flag
@ -342,9 +344,9 @@ module mpitools
if (ierror /= MPI_SUCCESS) &
call print_message(loc, "MPI_Allreduce of logical buffer failed!")
#endif /* MPI */
call stop_timer(imc)
#endif /* MPI */
!-------------------------------------------------------------------------------
!

39
sources/refinement.F90
View File

@ -102,7 +102,7 @@ module refinement
integer :: p
character(len=255) :: variables = "dens pres"
character(len=*), parameter :: loc = 'REFINEMENT:initialize_refinement()'
character(len=*), parameter :: loc = 'REFINEMENT::initialize_refinement()'
!-------------------------------------------------------------------------------
!
@ -217,7 +217,7 @@ module refinement
integer, intent(out) :: status
character(len=*), parameter :: loc = 'REFINEMENT:dinalize_refinement()'
character(len=*), parameter :: loc = 'REFINEMENT:finalize_refinement()'
!-------------------------------------------------------------------------------
!
@ -508,10 +508,15 @@ module refinement
real(kind=8), dimension(:,:,:,:), pointer, save :: w
character(len=*), parameter :: loc = 'REFINEMENT:vorticity_magnitude()'
integer, save :: nt = 0
!$ integer :: omp_get_thread_num
!$omp threadprivate(first,nt,dh,w)
character(len=*), parameter :: loc = 'REFINEMENT::vorticity_magnitude()'
!-------------------------------------------------------------------------------
!
!$ nt = omp_get_thread_num()
wmax = 0.0e+00
if (first) then
@ -524,9 +529,9 @@ module refinement
end if
#if NDIMS == 3
w(1:3,1:nn,1:nn,1:nn) => work(1:i)
w(1:3,1:nn,1:nn,1:nn) => work(1:i,nt)
#else /* NDIMS == 3 */
w(1:3,1:nn,1:nn,1: 1) => work(1:i)
w(1:3,1:nn,1:nn,1: 1) => work(1:i,nt)
#endif /* NDIMS == 3 */
dh(:) = 1.0d+00
@ -534,11 +539,11 @@ module refinement
first = .false.
end if
if (work_in_use) &
if (work_in_use(nt)) &
call print_message(loc, &
"Workspace is being used right now! Corruptions can occur!")
work_in_use = .true.
work_in_use(nt) = .true.
call curl(dh(:), pdata%q(ivx:ivz,:,:,:), w(:,:,:,:))
@ -558,7 +563,7 @@ module refinement
end do
#endif /* NDIMS == 3 */
work_in_use = .false.
work_in_use(nt) = .false.
wmax = sqrt(wmax)
@ -606,13 +611,19 @@ module refinement
real(kind=8), dimension(:,:,:,:), pointer, save :: w
character(len=*), parameter :: loc = 'REFINEMENT:current_density_magnitude()'
integer, save :: nt = 0
!$ integer :: omp_get_thread_num
!$omp threadprivate(first,nt,dh,w)
character(len=*), parameter :: loc = &
'REFINEMENT::current_density_magnitude()'
!-------------------------------------------------------------------------------
!
jmax = 0.0e+00
if (.not. magnetized) return
!$ nt = omp_get_thread_num()
if (first) then
i = 3 * nn**NDIMS
@ -624,9 +635,9 @@ module refinement
end if
#if NDIMS == 3
w(1:3,1:nn,1:nn,1:nn) => work(1:i)
w(1:3,1:nn,1:nn,1:nn) => work(1:i,nt)
#else /* NDIMS == 3 */
w(1:3,1:nn,1:nn,1: 1) => work(1:i)
w(1:3,1:nn,1:nn,1: 1) => work(1:i,nt)
#endif /* NDIMS == 3 */
dh(:) = 1.0d+00
@ -634,11 +645,11 @@ module refinement
first = .false.
end if
if (work_in_use) &
if (work_in_use(nt)) &
call print_message(loc, &
"Workspace is being used right now! Corruptions can occur!")
work_in_use = .true.
work_in_use(nt) = .true.
call curl(dh(:), pdata%q(ibx:ibz,:,:,:), w(:,:,:,:))
@ -658,7 +669,7 @@ module refinement
end do
#endif /* NDIMS == 3 */
work_in_use = .false.
work_in_use(nt) = .false.
jmax = sqrt(jmax)

27
sources/sources.F90
View File

@ -296,6 +296,9 @@ module sources
real(kind=8) :: dvydx, dvydy, dvydz
real(kind=8) :: dvzdx, dvzdy, dvzdz
integer :: nt = 0
!$ integer :: omp_get_thread_num
real(kind=8), dimension(3) :: ga, dh
real(kind=8), dimension(nn) :: x, y
#if NDIMS == 3
@ -305,11 +308,13 @@ module sources
#endif /* NDIMS == 3 */
real(kind=8), dimension(:,:,:) , pointer, save :: db
real(kind=8), dimension(:,:,:,:,:), pointer, save :: tmp
!$omp threadprivate(first, nt, db, tmp)
character(len=*), parameter :: loc = 'SOURCES:update_sources()'
!-------------------------------------------------------------------------------
!
!$ nt = omp_get_thread_num()
if (first) then
i = nn**NDIMS
j = 10 * i
@ -321,11 +326,11 @@ module sources
end if
#if NDIMS == 3
db(1:nn,1:nn,1:nn) => work( 1:i)
tmp(1:3,1:3,1:nn,1:nn,1:nn) => work(i+1:j)
db(1:nn,1:nn,1:nn) => work( 1:i,nt)
tmp(1:3,1:3,1:nn,1:nn,1:nn) => work(i+1:j,nt)
#else /* NDIMS == 3 */
db(1:nn,1:nn,1: 1) => work( 1:i)
tmp(1:3,1:3,1:nn,1:nn,1: 1) => work(i+1:j)
db(1:nn,1:nn,1: 1) => work( 1:i,nt)
tmp(1:3,1:3,1:nn,1:nn,1: 1) => work(i+1:j,nt)
#endif /* NDIMS == 3 */
first = .false.
@ -399,11 +404,11 @@ module sources
!
if (viscosity > 0.0d+00) then
if (work_in_use) &
if (work_in_use(nt)) &
call print_message(loc, "Workspace is being used right now! " // &
"Corruptions can occur!")
work_in_use = .true.
work_in_use(nt) = .true.
! prepare coordinate increments
!
@ -531,7 +536,7 @@ module sources
end if ! ien > 0
work_in_use = .false.
work_in_use(nt) = .false.
end if ! viscosity is not zero
@ -539,13 +544,11 @@ module sources
!
if (magnetized) then
if (work_in_use) &
if (work_in_use(nt)) &
call print_message(loc, "Workspace is being used right now! " // &
"Corruptions can occur!")
if (work_in_use) then
end if
work_in_use = .true.
work_in_use(nt) = .true.
! prepare coordinate increments
!
@ -720,7 +723,7 @@ module sources
end if ! resistivity is not zero
work_in_use = .false.
work_in_use(nt) = .false.
end if ! magnetized

38
sources/statistics.F90
View File

@ -55,7 +55,10 @@ module statistics
! the format strings
!
character(len=32), save :: mfmt1, mfmt2, mfmt3
character(len=32), save :: mfmt1, mfmt2
#ifdef MPI
character(len=32), save :: mfmt3
#endif /* MPI */
character(len=32), save :: efmt
private
@ -93,7 +96,10 @@ module statistics
use coordinates, only : toplev, ncells, bcells, nghosts, domain_base_dims
use equations , only : pvars, nf
use evolution , only : error_control
use mpitools , only : master, nprocs
use mpitools , only : master
#ifdef MPI
use mpitools , only : nprocs
#endif /* MPI */
use parameters , only : get_parameter
implicit none
@ -165,6 +171,8 @@ module statistics
#ifdef MPI
write(stmp,"(a,i0,a)") "(", max(1, 10 * toplev - 28), "x,'|',1x,a)"
write(munit,stmp) "block distribution per process"
#else /* MPI */
write(munit,"('')")
#endif /* MPI */
write(munit,"('#',75x,'|')", advance="no")
do l = 1, toplev
@ -385,11 +393,12 @@ module statistics
use forcing , only : einj, rinj, arms
use helpers , only : print_message, flush_and_sync
use mesh , only : changed
use mpitools , only : master, nprocs
use mpitools , only : master
#ifdef MPI
use mpitools , only : reduce_minimum, reduce_maximum, reduce_sum
use mpitools , only : nprocs
use mpitools , only : reduce_minimum, reduce_maximum, reduce_sum
#endif /* MPI */
use workspace , only : resize_workspace, work, work_in_use
use workspace , only : resize_workspace, work, work_in_use
implicit none
@ -414,10 +423,15 @@ module statistics
integer(kind=4), dimension(nprocs) :: cdist
#endif /* MPI */
integer :: nt = 0
!$ integer :: omp_get_thread_num
!$omp threadprivate(first, nt)
character(len=*), parameter :: loc = 'INTEGRALS:store_statistics()'
!-------------------------------------------------------------------------------
!
!$ nt = omp_get_thread_num()
! process and store the mesh statistics only on the master node
!
if (master) then
@ -483,16 +497,16 @@ module statistics
n = ni**NDIMS
l = 1
u = n
vel(1:ni,1:ni,1:nk) => work(l:u)
vel(1:ni,1:ni,1:nk) => work(l:u,nt)
l = l + n
u = u + n
mag(1:ni,1:ni,1:nk) => work(l:u)
mag(1:ni,1:ni,1:nk) => work(l:u,nt)
l = l + n
u = u + n
sqd(1:ni,1:ni,1:nk) => work(l:u)
sqd(1:ni,1:ni,1:nk) => work(l:u,nt)
l = l + n
u = u + n
tmp(1:ni,1:ni,1:nk) => work(l:u)
tmp(1:ni,1:ni,1:nk) => work(l:u,nt)
first = .false.
end if
@ -519,11 +533,11 @@ module statistics
mxarr(7) = 0.0d+00
end if
if (work_in_use) &
if (work_in_use(nt)) &
call print_message(loc, "Workspace is being used right now! " // &
"Corruptions can occur!")
work_in_use = .true.
work_in_use(nt) = .true.
! associate the pointer with the first block on the data block list
!
@ -692,7 +706,7 @@ module statistics
end do ! data blocks
work_in_use = .false.
work_in_use(nt) = .false.
#ifdef MPI
! sum the integral array from all processes

12
sources/system.F90
View File

@ -38,9 +38,11 @@ module system
!
logical, save :: resumed = .false.
#ifdef SIGNALS
! the flag indicating that the system evolution should be interrupted
!
integer, save :: quit = 0
#endif /* SIGNALS */
! the number of resumed run
!
@ -77,7 +79,10 @@ module system
public :: initialize_system, finalize_system, prepare_system, evolve_system
public :: print_system_info
public :: resumed, name, nrun, eqsys, eos, ncells, nghosts, maxlev, nsteps
public :: bdims, dbnds, rngtype, tsav, nwork, quit
public :: bdims, dbnds, rngtype, tsav, nwork
#ifdef SIGNALS
public :: quit
#endif /* SIGNALS */
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
@ -499,7 +504,10 @@ module system
if (check_status(status /= 0)) return
proceed = (nsteps <= nmax) .and. (time < tmax) .and. &
(thrs <= trun) .and. .not. check_status(quit /= 0)
(thrs <= trun)
#ifdef SIGNALS
proceed = proceed .and. .not. check_status(quit /= 0)
#endif /* SIGNALS */
if (verbose) then

35
sources/workspace.F90
View File

@ -33,13 +33,17 @@ module workspace
!
integer, save :: nwork = 0
! the last thread number
!
integer, save :: nt = 0
! the flag indicating that the workspace is in use
!
logical, save :: work_in_use = .false.
logical, dimension(:), allocatable, save :: work_in_use
! the common workspace to use for local arrays
!
real(kind=8), dimension(:), allocatable, target :: work
real(kind=8), dimension(:,:), allocatable, target :: work
private
@ -64,18 +68,19 @@ module workspace
!
! Arguments:
!
! ninit - the initial workspace size;
! status - the call status (0 for success, otherwise failure);
! ninit - the initial workspace size;
! nthreads - the number of threads;
! status - the call status (0 for success, otherwise failure);
!
!===============================================================================
!
subroutine initialize_workspace(ninit, status)
subroutine initialize_workspace(ninit, nthreads, status)
use helpers, only : print_message
implicit none
integer, intent(in) :: ninit
integer, intent(in) :: ninit, nthreads
integer, intent(out) :: status
character(len=*), parameter :: loc = 'WORKSPACE::initialize_workspace()'
@ -91,12 +96,15 @@ module workspace
end if
nwork = ninit
nt = nthreads - 1
allocate(work(nwork), stat=status)
allocate(work_in_use(0:nt), work(nwork,0:nt), stat=status)
if (status /= 0) &
call print_message(loc, "Could not allocate the common workspace!")
work_in_use = .false.
!-------------------------------------------------------------------------------
!
end subroutine initialize_workspace
@ -127,7 +135,7 @@ module workspace
status = 0
if (allocated(work)) then
deallocate(work, stat=status)
deallocate(work_in_use, work, stat=status)
if (status /= 0) &
call print_message(loc, "Could not deallocate the common workspace!")
@ -164,7 +172,13 @@ module workspace
!
status = 0
if (work_in_use) then
#ifdef _OPENMP
if (nsize > nwork) then
call print_message(loc, "Workspace is too small!")
status = 1
end if
#else /* _OPENMP */
if (any(work_in_use)) then
call print_message(loc, "Could not resize the workspace. " // &
"It is being used right now!")
@ -173,7 +187,7 @@ module workspace
if (nsize > nwork) then
deallocate(work, stat=status)
if (status == 0) then
allocate(work(nsize), stat=status)
allocate(work(nsize,0:nt), stat=status)
if (status /= 0) then
call print_message(loc, "Could not allocate a new workspace!")
status = 1
@ -188,6 +202,7 @@ module workspace
end if
end if
end if
#endif /* _OPENMP */
!-------------------------------------------------------------------------------
!