amun-code/src/coordinates.F90

!!******************************************************************************
!!
!!  This file is part of the AMUN source code, a program to perform
!!  Newtonian or relativistic magnetohydrodynamical simulations on uniform or
!!  adaptive mesh.
!!
!!  Copyright (C) 2008-2014 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!!  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 3 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, see <http://www.gnu.org/licenses/>.
!!
!!******************************************************************************
!!
!! module: COORDINATES
!!
!!  This module provides variables and subroutines handling the coordinates
!!  for all refinement levels.
!!
!!******************************************************************************
!
module coordinates

! module variables are not implicit by default
!
  implicit none

! MODULE PARAMETERS:
! =================
!
! the domain block dimensions
!
  integer     , save :: nn =  8, in =  8, jn =  8, kn =  1

! the number of ghost zones
!
  integer     , save :: ng =  2, nh =  1, nd =  4

! the domain block dimensions including the ghost zones
!
  integer     , save :: im = 12, jm = 12, km =  1

! the domain block dimensions including the ghost zones
!
  integer     , save :: it = 11, jt = 11, kt =  1

! the domain division
!
  integer     , save :: ir =  1, jr =  1, kr =  1

! the limits of refinement level
!
  integer     , save :: minlev = 1, maxlev = 1, toplev = 1

! block indices
!
  integer     , save :: ih  =  6, jh  =  6, kh  =  1
  integer     , save :: ib  =  3, jb  =  3, kb  =  1
  integer     , save :: ie  = 10, je  = 10, ke  =  1
  integer     , save :: ibl =  2, jbl =  2, kbl =  1
  integer     , save :: ibu =  4, jbu =  4, kbu =  1
  integer     , save :: iel =  9, jel =  9, kel =  1
  integer     , save :: ieu = 11, jeu = 11, keu =  1

! the domain bounds
!
  real(kind=8), save :: xmin = 0.0d+00
  real(kind=8), save :: xmax = 1.0d+00
  real(kind=8), save :: xlen = 1.0d+00
  real(kind=8), save :: ymin = 0.0d+00
  real(kind=8), save :: ymax = 1.0d+00
  real(kind=8), save :: ylen = 1.0d+00
  real(kind=8), save :: zmin = 0.0d+00
  real(kind=8), save :: zmax = 1.0d+00
  real(kind=8), save :: zlen = 1.0d+00

! the block coordinates for all levels of refinement
!
  real(kind=8), dimension(:,:), allocatable, save :: ax  , ay  , az
  real(kind=8), dimension(:  ), allocatable, save :: adx , ady , adz, adr
  real(kind=8), dimension(:  ), allocatable, save :: adxi, adyi, adzi
  real(kind=8), dimension(:  ), allocatable, save :: advol

! by default everything is private
!
  public

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
  contains
!
!===============================================================================
!
! subroutine INITIALIZE_COORDINATES:
! ---------------------------------
!
!   Subroutine initializes mesh coordinates and other coordinate parameters.
!
!   Arguments:
!
!     verbose - flag determining if the subroutine should be verbose;
!     iret    - return flag of the procedure execution status;
!
!===============================================================================
!
  subroutine initialize_coordinates(verbose, iret)

! include external procedures
!
    use parameters, only : get_parameter_integer, get_parameter_real

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    logical, intent(in)    :: verbose
    integer, intent(inout) :: iret

! local variables
!
    integer :: i, j, k, l, ff
    integer :: ni, nj, nk
    logical :: info

! local arrays
!
    integer(kind=4), dimension(3) :: cm, rm, dm
!
!-------------------------------------------------------------------------------
!
! obtain the number of cells along each block dimension in the case when all
! dimensions are the same
!
    call get_parameter_integer("ncells", nn    )

! set the block dimensions
!
    in = nn
    jn = nn
#if NDIMS == 3
    kn = nn
#endif /* NDIMS == 3 */

! change individual block dimension if requested
!
    call get_parameter_integer("icells", in    )
    call get_parameter_integer("jcells", jn    )
#if NDIMS == 3
    call get_parameter_integer("kcells", kn    )
#endif /* NDIMS == 3 */

! obtain the number of ghost zones
!
    call get_parameter_integer("nghost", ng    )

! calculate half and double of the number of ghose zones
!
    nh = ng / 2
    nd = ng * 2

! calculate the block dimensions including ghost cells
!
    im = in + 2 * ng
    jm = jn + 2 * ng
#if NDIMS == 3
    km = kn + 2 * ng
#endif /* NDIMS == 3 */

! prepare indices
!
    it = im - nh + 1
    jt = jm - nh + 1
#if NDIMS == 3
    kt = km - nh + 1
#endif /* NDIMS == 3 */

! calculate block indices
!
    ih  = im / 2
    ib  = ng +  1
    ie  = ng + in
    ibl = ib - 1
    ibu = ib + ng - 1
    iel = ie - ng + 1
    ieu = ie + 1
    jh  = jm / 2
    jb  = ng +  1
    je  = ng + jn
    jbl = jb - 1
    jbu = jb + ng - 1
    jel = je - ng + 1
    jeu = je + 1
#if NDIMS == 3
    kh  = km / 2
    kb  = ng +  1
    ke  = ng + kn
    kbl = kb - 1
    kbu = kb + ng - 1
    kel = ke - ng + 1
    keu = ke + 1
#endif /* NDIMS == 3 */

! obtain the refinement level bounds
!
    call get_parameter_integer("minlev", minlev)
    call get_parameter_integer("maxlev", maxlev)

! set the top level
!
    toplev = maxlev

! obtain the domain base division
!
    call get_parameter_integer("rdimx" , ir    )
    call get_parameter_integer("rdimy" , jr    )
#if NDIMS == 3
    call get_parameter_integer("rdimz" , kr    )
#endif /* NDIMS == 3 */

! obtain the domain bounds
!
    call get_parameter_real   ("xmin"  , xmin  )
    call get_parameter_real   ("xmax"  , xmax  )
    call get_parameter_real   ("ymin"  , ymin  )
    call get_parameter_real   ("ymax"  , ymax  )
#if NDIMS == 3
    call get_parameter_real   ("zmin"  , zmin  )
    call get_parameter_real   ("zmax"  , zmax  )
#endif /* NDIMS == 3 */

! allocate space for coordinate variables
!
    allocate(ax   (toplev, im))
    allocate(ay   (toplev, jm))
    allocate(az   (toplev, km))
    allocate(adx  (toplev))
    allocate(ady  (toplev))
    allocate(adz  (toplev))
    allocate(adr  (toplev))
    allocate(adxi (toplev))
    allocate(adyi (toplev))
    allocate(adzi (toplev))
    allocate(advol(toplev))

! reset all coordinate variables to initial values
!
    ax   (:,:) = 0.0d+00
    ay   (:,:) = 0.0d+00
    az   (:,:) = 0.0d+00
    adx  (:)   = 1.0d+00
    ady  (:)   = 1.0d+00
    adz  (:)   = 1.0d+00
    adr  (:)   = 1.0d+00
    adxi (:)   = 1.0d+00
    adyi (:)   = 1.0d+00
    adzi (:)   = 1.0d+00
    advol(:)   = 1.0d+00

! generate the coordinate variables for each level
!
    do l = 1, toplev

! calculate the block resolution at each level
!
      ff      = 2**(l - 1)
      ni      = in * ff
      nj      = jn * ff
      nk      = kn * ff

! calculate the cell sizes for each level
!
      adx (l) = (xmax - xmin) / (ir * ni)
      ady (l) = (ymax - ymin) / (jr * nj)
#if NDIMS == 3
      adz (l) = (zmax - zmin) / (kr * nk)
#endif /* NDIMS == 3 */
#if NDIMS == 2
      adr (l) = sqrt(adx(l)**2 + ady(l)**2)
#endif /* NDIMS == 2 */
#if NDIMS == 3
      adr (l) = sqrt(adx(l)**2 + ady(l)**2 + adz(l)**2)
#endif /* NDIMS == 3 */

! calculate the inverse of cell size
!
      adxi(l) = 1.0d+00 / adx(l)
      adyi(l) = 1.0d+00 / ady(l)
#if NDIMS == 3
      adzi(l) = 1.0d+00 / adz(l)
#endif /* NDIMS == 3 */

! calculate the block coordinates for each level
!
      ax(l,:) = ((/(i, i = 1, im)/) - ng - 5.0d-01) * adx(l)
      ay(l,:) = ((/(j, j = 1, jm)/) - ng - 5.0d-01) * ady(l)
#if NDIMS == 3
      az(l,:) = ((/(k, k = 1, km)/) - ng - 5.0d-01) * adz(l)
#endif /* NDIMS == 3 */

! calculate the cell volume at each level
!
      advol(l) = adx(l) * ady(l) * adz(l)

    end do ! l = 1, toplev

! print general information about the level resolutions
!
    if (verbose) then

! the base resolution
!
      cm(1) = ir * in
      cm(2) = jr * jn
      cm(3) = kr * kn

! the effective resolution
!
      ff    = 2**(maxlev - 1)
      rm(1) = cm(1) * ff
      rm(2) = cm(2) * ff
      rm(3) = cm(3) * ff

! the top level block division
!
      dm(1) = rm(1) / in
      dm(2) = rm(2) / jn
      dm(3) = rm(3) / kn

! obtain the maximum number of block
!
      ff    = product(dm(1:NDIMS))

! print info
!
      write(*,"(4x,a,  1x,i6 )" ) "refinement to level    =", toplev
      write(*,"(4x,a,3(1x,i6 ))") "base configuration     =", ir, jr, kr
      write(*,"(4x,a,3(1x,i6 ))") "top level blocks       =", dm(1:NDIMS)
      write(*,"(4x,a,  3x,i18)" ) "maximum cover blocks   =", ff
      write(*,"(4x,a,3(1x,i6 ))") "base resolution        =", cm(1:NDIMS)
      write(*,"(4x,a,3(1x,i6 ))") "effective resolution   =", rm(1:NDIMS)

    end if ! verbose

!-------------------------------------------------------------------------------
!
  end subroutine initialize_coordinates
!
!===============================================================================
!
! subroutine FINALIZE_COORDINATES:
! -------------------------------
!
!   Subroutine deallocates mesh coordinates.
!
!   Arguments:
!
!     iret    - return flag of the procedure execution status;
!
!===============================================================================
!
  subroutine finalize_coordinates(iret)

! local variables are not implicit by default
!
    implicit none

! subroutine arguments
!
    integer, intent(inout) :: iret
!
!-------------------------------------------------------------------------------
!
! 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(adr)  ) deallocate(adr)
    if (allocated(adxi) ) deallocate(adxi)
    if (allocated(adyi) ) deallocate(adyi)
    if (allocated(adzi) ) deallocate(adzi)
    if (allocated(advol)) deallocate(advol)

!-------------------------------------------------------------------------------
!
  end subroutine finalize_coordinates

!===============================================================================
!
end module