diff --git a/src/coords.F90 b/src/coords.F90 new file mode 100644 index 0000000..e9867a6 --- /dev/null +++ b/src/coords.F90 @@ -0,0 +1,234 @@ +!!****************************************************************************** +!! +!! module: COORDS - handles coordinates +!! +!! Copyright (C) 2011 Grzegorz Kowal +!! +!!****************************************************************************** +!! +!! This file is part of the AMUN code. +!! +!! This program is free software; you can redistribute it and/or +!! modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation; either version 2 +!! of the License, or (at your option) any later version. +!! +!! This program is distributed in the hope that it will be useful, +!! but WITHOUT ANY WARRANTY; without even the implied warranty of +!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +!! GNU General Public License for more details. +!! +!! You should have received a copy of the GNU General Public License +!! along with this program; if not, write to the Free Software Foundation, +!! Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +!! +!!****************************************************************************** +!! +! +module coords + + implicit none + +! the effective resolution of the full domain +! + integer,dimension(NDIMS), save :: effres = 1 + +! the block coordinates for all levels of refinement +! + real, dimension(:,:), allocatable, save :: ax , ay , az + real, dimension(: ), allocatable, save :: adx , ady , adz + real, dimension(: ), allocatable, save :: adxi, adyi, adzi + real, dimension(: ), allocatable, save :: advol + +! the block resolution in the units of effective resolution for all levels +! + integer(kind=4), dimension(:,:), allocatable, save :: res + + contains +! +!=============================================================================== +! +! init_coords: subroutine allocates and initializes coordinates +! +!=============================================================================== +! + subroutine init_coords(flag) + + use config, only : maxlev, ng, in, jn, kn, im, jm, km, rdims + use config, only : xmin, xmax, ymin, ymax, zmin, zmax + use timer , only : start_timer, stop_timer + + implicit none + +! input arguments +! + logical, intent(in) :: flag + +! local variables +! + integer :: i, j, k, l + integer :: ni, nj, nk + logical :: info + +! local arrays +! + integer(kind=4), dimension(3) :: bm, cm, rm, dm +! +!------------------------------------------------------------------------------- +! +! start the coordinate timer +! + call start_timer(7) + +! allocate space for coordinate variables and resolutions +! + allocate(ax (maxlev, im)) + allocate(ay (maxlev, jm)) + allocate(az (maxlev, km)) + allocate(adx (maxlev)) + allocate(ady (maxlev)) + allocate(adz (maxlev)) + allocate(adxi (maxlev)) + allocate(adyi (maxlev)) + allocate(adzi (maxlev)) + allocate(advol(maxlev)) + allocate(res (maxlev, 3)) + +! reset all coordinate variables to initial values +! + ax (:,:) = 0.0d0 + ay (:,:) = 0.0d0 + az (:,:) = 0.0d0 + adx (:) = 1.0d0 + ady (:) = 1.0d0 + adz (:) = 1.0d0 + adxi (:) = 1.0d0 + adyi (:) = 1.0d0 + adzi (:) = 1.0d0 + advol(:) = 1.0d0 + res (:,:) = 1 + +! generate the coordinate variables for each level +! + do l = 1, maxlev + +! calculate the block resolution at each level +! + ni = in * 2**(l - 1) + nj = jn * 2**(l - 1) + nk = kn * 2**(l - 1) + +! calculate the cell size at each level +! + adx (l) = (xmax - xmin) / (rdims(1) * ni) + ady (l) = (ymax - ymin) / (rdims(2) * nj) +#if NDIMS == 3 + adz (l) = (zmax - zmin) / (rdims(3) * nk) +#endif /* NDIMS == 3 */ + +! calculate the inverse of cell size at each level +! + adxi(l) = 1.0d0 / adx(l) + adyi(l) = 1.0d0 / ady(l) +#if NDIMS == 3 + adzi(l) = 1.0d0 / adz(l) +#endif /* NDIMS == 3 */ + +! calculate the block coordinates at each level +! + ax(l,:) = ((/(i, i = 1, im)/) - ng - 0.5d0) * adx(l) + ay(l,:) = ((/(j, j = 1, jm)/) - ng - 0.5d0) * ady(l) +#if NDIMS == 3 + az(l,:) = ((/(k, k = 1, km)/) - ng - 0.5d0) * adz(l) +#endif /* NDIMS == 3 */ + +! calculate the cell volume at each level +! + advol(l) = adx(l) * ady(l) * adz(l) + +! calculate the effective block resolution at each level +! + res(l,1) = in * 2**(maxlev - l) + res(l,2) = jn * 2**(maxlev - l) +#if NDIMS == 3 + res(l,3) = kn * 2**(maxlev - l) +#endif /* NDIMS == 3 */ + + end do + +! print general information about resolutions +! + if (flag) then + + bm( : ) = (/ in, jn, kn /) + rm( : ) = 1 + dm( : ) = 1 + + cm(1:NDIMS) = rdims(1:NDIMS) * bm(1:NDIMS) + rm(1:NDIMS) = rdims(1:NDIMS) * res(1,1:NDIMS) + dm(1:NDIMS) = rm(1:NDIMS) / bm(1:NDIMS) + + effres(1:NDIMS) = rm(1:NDIMS) + + write(*,*) + write(*,"(1x,a)" ) "Geometry:" + write(*,"(4x,a, 1x,i6)" ) "refinement to level =", maxlev + write(*,"(4x,a,3(1x,i6))") "base configuration =", rdims(1:NDIMS) + write(*,"(4x,a,3(1x,i6))") "top level blocks =", dm(1:NDIMS) + write(*,"(4x,a, 1x,i6)" ) "maxium cover blocks =", product(dm(:)) + write(*,"(4x,a,3(1x,i6))") "base resolution =", cm(1:NDIMS) + write(*,"(4x,a,3(1x,i6))") "effective resolution =", rm(1:NDIMS) + + end if ! master + +! stop the coordinate timer +! + call stop_timer(7) + +!------------------------------------------------------------------------------- +! + end subroutine init_coords +! +!=============================================================================== +! +! clears_coords: subroutine deallocates coordinate arrays +! +!=============================================================================== +! + subroutine clear_coords() + + use timer, only : start_timer, stop_timer + + implicit none + +!------------------------------------------------------------------------------- +! +! start the coordinate timer +! + call start_timer(7) + +! deallocating coordinate variables +! + if (allocated(ax) ) deallocate(ax) + if (allocated(ay) ) deallocate(ay) + if (allocated(az) ) deallocate(az) + if (allocated(adx) ) deallocate(adx) + if (allocated(ady) ) deallocate(ady) + if (allocated(adz) ) deallocate(adz) + if (allocated(adxi) ) deallocate(adxi) + if (allocated(adyi) ) deallocate(adyi) + if (allocated(adzi) ) deallocate(adzi) + if (allocated(advol)) deallocate(advol) + if (allocated(res) ) deallocate(res) + +! stop the coordinate timer +! + call stop_timer(7) + +!------------------------------------------------------------------------------- +! + end subroutine clear_coords + +!=============================================================================== +! +end module diff --git a/src/driver.F90 b/src/driver.F90 index ba3ea3f..e44be4a 100644 --- a/src/driver.F90 +++ b/src/driver.F90 @@ -34,6 +34,7 @@ program amun #ifdef FORCE use config , only : fdt #endif /* FORCE */ + use coords , only : init_coords, clear_coords use evolution, only : evolve, find_new_timestep, n, t, dt, dtn #ifdef FORCE use forcing , only : init_forcing, clear_forcing @@ -148,6 +149,10 @@ program amun ! call init_blocks() +! initialize the coordinate module +! + call init_coords(is_master()) + ! check if we initiate new problem or restart previous job ! if (nres .lt. 0) then @@ -329,6 +334,10 @@ program amun ! call clear_mesh() +! finalize the coordinate module +! + call clear_coords() + ! get total time ! tall = get_timer_total() diff --git a/src/evolution.F90 b/src/evolution.F90 index eca6176..92e7135 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -197,7 +197,7 @@ module evolution #ifdef MPI use mpitools, only : mallreducemaxr #endif /* MPI */ - use mesh , only : adx, ady, adz + use coords , only : adx, ady, adz use scheme , only : maxspeed, cmax use timer , only : start_timer, stop_timer #ifdef VISCOSITY @@ -335,7 +335,7 @@ module evolution use blocks , only : block_data use config , only : im, jm, km - use mesh , only : adxi, adyi, adzi + use coords , only : adxi, adyi, adzi #if defined MHD && defined GLM use config , only : alpha_p use scheme , only : cmax @@ -556,7 +556,7 @@ module evolution subroutine flux_euler(pblock) use blocks , only : block_data - use mesh , only : adxi, adyi, adzi + use coords , only : adxi, adyi, adzi use scheme , only : update_flux implicit none @@ -598,7 +598,7 @@ module evolution use blocks , only : block_data use config , only : im, jm, km - use mesh , only : adxi, adyi, adzi + use coords , only : adxi, adyi, adzi use scheme , only : update_flux use variables, only : nqt @@ -662,7 +662,7 @@ module evolution use blocks , only : block_data use config , only : im, jm, km - use mesh , only : adxi, adyi, adzi + use coords , only : adxi, adyi, adzi use scheme , only : update_flux use variables, only : nqt @@ -746,7 +746,7 @@ module evolution #ifdef FORCE use forcing , only : real_forcing #endif /* FORCE */ - use mesh , only : adxi, adyi, adzi + use coords , only : adxi, adyi, adzi #ifdef SHAPE use problem , only : update_shapes #endif /* SHAPE */ @@ -882,7 +882,7 @@ module evolution #ifdef FORCE use forcing , only : real_forcing #endif /* FORCE */ - use mesh , only : adxi, adyi, adzi + use coords , only : adxi, adyi, adzi #ifdef SHAPE use problem , only : update_shapes #endif /* SHAPE */ @@ -1043,7 +1043,7 @@ module evolution #ifdef FORCE use forcing , only : real_forcing #endif /* FORCE */ - use mesh , only : adxi, adyi, adzi + use coords , only : adxi, adyi, adzi #ifdef SHAPE use problem , only : update_shapes #endif /* SHAPE */ diff --git a/src/forcing.F90 b/src/forcing.F90 index 49fa455..9319bce 100644 --- a/src/forcing.F90 +++ b/src/forcing.F90 @@ -499,7 +499,7 @@ module forcing use config , only : im, jm, km, ib, ie, jb, je, kb, ke use constants, only : dpi - use mesh , only : ax, ay, az, advol + use coords , only : ax, ay, az, advol use timer , only : start_timer, stop_timer use variables, only : idn, imx, imy, imz @@ -664,7 +664,7 @@ module forcing use config , only : im, jm, km, ib, ie, jb, je, kb, ke use constants, only : dpi - use mesh , only : ax, ay, az, advol + use coords , only : ax, ay, az, advol use timer , only : start_timer, stop_timer implicit none diff --git a/src/integrals.F90 b/src/integrals.F90 index 924272e..e238b93 100644 --- a/src/integrals.F90 +++ b/src/integrals.F90 @@ -141,7 +141,7 @@ module integrals #ifdef FORCE use forcing , only : fcor, finp #endif /* FORCE */ - use mesh , only : advol + use coords , only : advol use mpitools , only : is_master, mallreducesumd use variables, only : idn, imx, imy, imz #ifdef ADI diff --git a/src/io.F90 b/src/io.F90 index f800ab2..a7aedcb 100644 --- a/src/io.F90 +++ b/src/io.F90 @@ -2240,7 +2240,7 @@ module io use error , only : print_error use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5gcreate_f, h5gclose_f - use mesh , only : adx, ady, adz, res + use coords, only : adx, ady, adz, res ! declare variables ! diff --git a/src/makefile b/src/makefile index 4a1458a..d7fe721 100644 --- a/src/makefile +++ b/src/makefile @@ -156,11 +156,11 @@ name = amun default: $(name).x -sources = blocks.F90 boundaries.F90 config.F90 constants.F90 driver.F90 \ - error.F90 evolution.F90 forcing.F90 integrals.F90 interpolation.F90 \ - io.F90 mesh.F90 mpitools.F90 problem.F90 random.F90 scheme.F90 \ - timer.F90 variables.F90 -objects = blocks.o boundaries.o config.o constants.o driver.o error.o \ +sources = blocks.F90 boundaries.F90 config.F90 constants.F90 coords.F90 \ + driver.F90 error.F90 evolution.F90 forcing.F90 integrals.F90 \ + interpolation.F90 io.F90 mesh.F90 mpitools.F90 problem.F90 \ + random.F90 scheme.F90 timer.F90 variables.F90 +objects = blocks.o boundaries.o config.o constants.o coords.o driver.o error.o \ evolution.o forcing.o integrals.o interpolation.o io.o mesh.o \ mpitools.o problem.o random.o scheme.o timer.o variables.o files = $(sources) makefile make.default config.in license.txt hosts @@ -201,20 +201,21 @@ boundaries.o : boundaries.F90 blocks.o config.o error.o interpolation.o \ mpitools.o variables.o timer.o config.o : config.F90 error.o constants.o : constants.F90 -driver.o : driver.F90 blocks.o config.o evolution.o forcing.o \ +coords.o : coords.F90 config.o timer.o +driver.o : driver.F90 blocks.o config.o coords.o evolution.o forcing.o \ integrals.o io.o mesh.o mpitools.o random.o timer.o error.o : error.F90 -evolution.o : evolution.F90 blocks.o boundaries.o config.o forcing.o \ - interpolation.o mesh.o mpitools.o problem.o scheme.o timer.o \ - variables.o -forcing.o : forcing.F90 config.o constants.o mesh.o mpitools.o random.o \ - timer.o variables.o -integrals.o : integrals.F90 blocks.o config.o evolution.o mesh.o \ +evolution.o : evolution.F90 blocks.o boundaries.o config.o coords.o \ + forcing.o interpolation.o mesh.o mpitools.o problem.o \ + scheme.o timer.o variables.o +forcing.o : forcing.F90 config.o constants.o coords.o mpitools.o \ + random.o timer.o variables.o +integrals.o : integrals.F90 blocks.o config.o coords.o evolution.o \ mpitools.o variables.o interpolation.o : interpolation.F90 blocks.o config.o variables.o -io.o : io.F90 blocks.o error.o evolution.o mesh.o mpitools.o \ +io.o : io.F90 blocks.o coords.o error.o evolution.o mpitools.o \ random.o scheme.o variables.o -mesh.o : mesh.F90 blocks.o config.o error.o interpolation.o \ +mesh.o : mesh.F90 blocks.o config.o coords.o error.o interpolation.o \ mpitools.o problem.o variables.o mpitools.o : mpitools.F90 problem.o : problem.F90 blocks.o constants.o mpitools.o random.o \ diff --git a/src/mesh.F90 b/src/mesh.F90 index e284382..6846049 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -29,21 +29,9 @@ module mesh implicit none -! spatial coordinates for all levels of refinements -! - real, dimension(:,:), allocatable, save :: ax , ay , az - real, dimension(: ), allocatable, save :: adx , ady , adz - real, dimension(: ), allocatable, save :: adxi, adyi, adzi - real, dimension(: ), allocatable, save :: advol - ! maximum number of covering blocks ! integer , save :: tblocks = 1 - integer,dimension(NDIMS), save :: effres = 1 - -! the effective resolution for all levels -! - integer(kind=4), dimension(:,:), allocatable, save :: res ! log file for the mesh statistics ! @@ -92,93 +80,10 @@ module mesh ! call datablock_set_dims(nqt, nvr, im, jm, km) -! allocating space for coordinate variables -! - allocate(ax (maxlev, im)) - allocate(ay (maxlev, jm)) - allocate(az (maxlev, km)) - allocate(adx (maxlev)) - allocate(ady (maxlev)) - allocate(adz (maxlev)) - allocate(adxi (maxlev)) - allocate(adyi (maxlev)) - allocate(adzi (maxlev)) - allocate(advol(maxlev)) - allocate(res (maxlev, 3)) - -! reset the coordinate variables -! - ax(:,:) = 0.0d0 - ay(:,:) = 0.0d0 - az(:,:) = 0.0d0 - adx(:) = 1.0d0 - ady(:) = 1.0d0 - adz(:) = 1.0d0 - adxi(:) = 1.0d0 - adyi(:) = 1.0d0 - adzi(:) = 1.0d0 - advol(:) = 1.0d0 - res(:,:) = 1 - -! generating coordinates for all levels -! - do l = 1, maxlev - - n = ncells * 2**(l - 1) - -! calculate the effective resolution at each level -! - res(l,1) = in * 2**(maxlev - l) - res(l,2) = jn * 2**(maxlev - l) -#if NDIMS == 3 - res(l,3) = kn * 2**(maxlev - l) -#endif /* NDIMS == 3 */ - - adx (l) = (xmax - xmin) / (rdims(1) * n) - ady (l) = (ymax - ymin) / (rdims(2) * n) -#if NDIMS == 3 - adz (l) = (zmax - zmin) / (rdims(3) * n) -#endif /* NDIMS == 3 */ - - ax(l,:) = ((/(i, i = 1, im)/) - ng - 0.5d0) * adx(l) - ay(l,:) = ((/(j, j = 1, jm)/) - ng - 0.5d0) * ady(l) -#if NDIMS == 3 - az(l,:) = ((/(k, k = 1, km)/) - ng - 0.5d0) * adz(l) -#endif /* NDIMS == 3 */ - - adxi(l) = 1.0d0 / adx(l) - adyi(l) = 1.0d0 / ady(l) -#if NDIMS == 3 - adzi(l) = 1.0d0 / adz(l) -#endif /* NDIMS == 3 */ - - advol(l) = adx(l) * ady(l) * adz(l) - end do - ! print general information about resolutions ! if (is_master()) then - bm( : ) = 1 - rm( : ) = 1 - dm( : ) = 1 - - bm(1:NDIMS) = rdims(1:NDIMS) * ncells - rm(1:NDIMS) = rdims(1:NDIMS) * res(1,1:NDIMS) - dm(1:NDIMS) = rm(1:NDIMS) / ncells - - effres(1:NDIMS) = rm(1:NDIMS) - tblocks = product(dm(:)) - - write(*,*) - write(*,"(1x,a)" ) "Geometry:" - write(*,"(4x,a, 1x,i6)" ) "refinement to level =", maxlev - write(*,"(4x,a,3(1x,i6))") "base configuration =", rdims(1:NDIMS) - write(*,"(4x,a,3(1x,i6))") "top level blocks =", dm(1:NDIMS) - write(*,"(4x,a, 1x,i6)" ) "maxium cover blocks =", tblocks - write(*,"(4x,a,3(1x,i6))") "base resolution =", bm(1:NDIMS) - write(*,"(4x,a,3(1x,i6))") "effective resolution =", rm(1:NDIMS) - ! prepare the file for logging mesh statistics; only the master process handles ! this part ! @@ -248,11 +153,12 @@ module mesh ! subroutine generate_mesh() - use config , only : maxlev, rdims use blocks , only : block_meta, block_data, list_meta, list_data use blocks , only : refine_block, deallocate_datablock use blocks , only : nchild, nsides, nfaces use blocks , only : get_mblocks, get_nleafs + use config , only : maxlev, rdims + use coords , only : res use error , only : print_info, print_error use mpitools, only : is_master, ncpu, ncpus use problem , only : init_domain, init_problem, check_ref @@ -526,12 +432,13 @@ module mesh ! subroutine update_mesh() - use config , only : maxlev, im, jm, km use blocks , only : block_meta, block_data, list_meta, list_data & , nchild, ndims, nsides, nfaces & , refine_block, derefine_block, append_datablock & , associate_blocks, deallocate_datablock use blocks , only : get_nleafs + use config , only : maxlev, im, jm, km + use coords , only : res use error , only : print_info, print_error #ifdef MPI use mpitools , only : ncpus, ncpu, is_master, mallreducesuml, msendf, mrecvf @@ -1022,7 +929,7 @@ module mesh ! local variables ! integer(kind=4) :: l, n -! +! ! tag for the MPI data exchange ! integer(kind=4) :: itag @@ -1452,20 +1359,6 @@ module mesh ! call clear_blocks -! deallocating coordinate variables -! - if (allocated(ax) ) deallocate(ax) - if (allocated(ay) ) deallocate(ay) - if (allocated(az) ) deallocate(az) - if (allocated(adx) ) deallocate(adx) - if (allocated(ady) ) deallocate(ady) - if (allocated(adz) ) deallocate(adz) - if (allocated(adxi) ) deallocate(adxi) - if (allocated(adyi) ) deallocate(adyi) - if (allocated(adzi) ) deallocate(adzi) - if (allocated(advol)) deallocate(advol) - if (allocated(res) ) deallocate(res) - ! close the handler of the mesh statistics file ! if (is_master()) close(funit) @@ -1489,6 +1382,7 @@ module mesh use blocks , only : block_meta, list_meta use blocks , only : get_mblocks, get_nleafs use config , only : ncells, nghost, maxlev + use coords , only : effres use mpitools, only : is_master, ncpus implicit none