Move coordinate variables from MESH to new COORDS module.

- a new module COORDS handles the mesh variables which needed to be
   separated from the MESH module since they are used in PROBLEM module,
   which is required by MESH module; this created a circular dependency;
   by introducing a new COORDS module we removed that problem;
This commit is contained in:
Grzegorz Kowal 2011-05-11 15:32:01 -03:00
parent 81ba4935d2
commit 0fc7717100
8 changed files with 276 additions and 138 deletions

234
src/coords.F90 Normal file
View File

@ -0,0 +1,234 @@
!!******************************************************************************
!!
!! module: COORDS - handles coordinates
!!
!! Copyright (C) 2011 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!!******************************************************************************
!!
!! 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

View File

@ -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()

View File

@ -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 */

View File

@ -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

View File

@ -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

View File

@ -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
!

View File

@ -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 \

View File

@ -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
@ -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