!!****************************************************************************** !! !! 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-2019 Grzegorz Kowal !! !! 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 . !! !!****************************************************************************** !! !! module: DOMAINS !! !! This module handles the initialization of the problem domains. !! !! !!****************************************************************************** ! module domains ! module variables are not implicit by default ! implicit none ! interfaces for procedure pointers ! abstract interface subroutine setup_domain_iface() end subroutine end interface ! pointer to the problem setup subroutine ! procedure(setup_domain_iface), pointer, save :: setup_domain => null() ! by default everything is private ! private ! declare public subroutines ! public :: initialize_domains, finalize_domains public :: setup_domain !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== !! !!*** PUBLIC SUBROUTINES ***************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine INITIALIZE_DOMAINS: ! ----------------------------- ! ! Subroutine prepares module DOMAINS. ! ! Arguments: ! ! problem - the problem name ! verbose - a logical flag turning the information printing; ! iret - an integer flag for error return value; ! !=============================================================================== ! subroutine initialize_domains(problem, verbose, iret) ! include external procedures and variables ! use parameters, only : get_parameter ! local variables are not implicit by default ! implicit none ! subroutine arguments ! character(len=64), intent(in) :: problem logical , intent(in) :: verbose integer , intent(inout) :: iret ! !------------------------------------------------------------------------------- ! ! associate the setup_domain pointer with the respective problem setup ! subroutine ! select case(trim(problem)) case default setup_domain => setup_domain_default end select !------------------------------------------------------------------------------- ! end subroutine initialize_domains ! !=============================================================================== ! ! subroutine FINALIZE_DOMAINS: ! --------------------------- ! ! Subroutine releases memory used by the module. ! ! Arguments: ! ! iret - an integer flag for error return value; ! !=============================================================================== ! subroutine finalize_domains(iret) ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(inout) :: iret ! !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- ! end subroutine finalize_domains ! !=============================================================================== !! !!*** PRIVATE SUBROUTINES **************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine SETUP_DOMAIN_DEFAULT: ! ------------------------------- ! ! Subroutine sets the default domain of N₁xN₂xN₃ blocks in the right ! configuration. ! ! !=============================================================================== ! subroutine setup_domain_default() ! include external procedures and variables ! use blocks , only : pointer_meta, block_meta, block_data use blocks , only : append_metablock, append_datablock, link_blocks use blocks , only : metablock_set_leaf, metablock_set_level use blocks , only : metablock_set_configuration use blocks , only : metablock_set_coordinates, metablock_set_bounds use blocks , only : nsides use coordinates , only : xmin, ymin, zmin, xlen, ylen, zlen use coordinates , only : ddims => domain_base_dims use coordinates , only : periodic ! local variables are not implicit by default ! implicit none ! local variables ! integer :: i, j, k, n, p, ic, jc, kc real(kind=8) :: xl, xmn, xmx, yl, ymn, ymx, zl, zmn, zmx ! local arrays ! integer, dimension(3) :: loc, del ! local pointers ! type(block_meta), pointer :: pmeta, pnext ! allocatable arrays ! integer, dimension(:,:,:), allocatable :: cfg integer, dimension(:) , allocatable :: im, jm, km integer, dimension(:) , allocatable :: ip, jp, kp ! local pointer array ! type(pointer_meta), dimension(:,:,:), allocatable :: block_array ! !------------------------------------------------------------------------------- ! ! obtain the number of blocks ! n = product(ddims(1:NDIMS)) !! PREPARE BLOCK CONFIGURATION ARRAY !! ! allocate the configuration array ! allocate(cfg(ddims(1),ddims(2),ddims(3))) ! set the block configurations ! cfg(1:ddims(1),1:ddims(2):2,1:ddims(3):2) = 12 if (ddims(2) > 1) then cfg(1:ddims(1),2:ddims(2):2,1:ddims(3):2) = 43 cfg( ddims(1),1:ddims(2) ,1:ddims(3):2) = 13 end if if (ddims(3) > 1) then cfg(1:ddims(1),1:ddims(2):2,2:ddims(3):2) = 65 if (ddims(2) > 1) then cfg(1:ddims(1),2:ddims(2):2,2:ddims(3):2) = 78 cfg( ddims(1),1:ddims(2) ,2:ddims(3):2) = 75 end if if (ddims(1) == 1 .or. mod(ddims(2),2) == 1) then cfg( ddims(1), ddims(2) ,1:ddims(3) ) = 15 else cfg( 1 , ddims(2) ,1:ddims(3) ) = 48 end if end if !! ALLOCATE AND GENERATE META BLOCK CHAIN AND SET BLOCK CONFIGURATIONS !! ! allocate the block pointer array ! allocate(block_array(ddims(1),ddims(2),ddims(3))) ! generate the gray code for a given configuration and link the block in ! the proper order ! loc(:) = (/ 0, 0, 0 /) del(:) = (/ 1, 1, 1 /) p = 1 do k = 1, ddims(3) if (del(3) == 1) loc(3) = loc(3) + del(3) do j = 1, ddims(2) if (del(2) == 1) loc(2) = loc(2) + del(2) do i = 1, ddims(1) if (del(1) == 1) loc(1) = loc(1) + del(1) ! append a new metablock ! call append_metablock(block_array(loc(1),loc(2),loc(3))%ptr) ! set the configuration type ! call metablock_set_configuration( & block_array(loc(1),loc(2),loc(3))%ptr & , cfg(loc(1),loc(2),loc(3))) ! increase the block number ! p = p + 1 if (del(1) == -1) loc(1) = loc(1) + del(1) end do if (del(2) == -1) loc(2) = loc(2) + del(2) del(1) = - del(1) end do if (del(3) == -1) loc(3) = loc(3) + del(3) del(2) = - del(2) end do ! deallocate the configuration array ! deallocate(cfg) !! FILL OUT THE REMAINING FIELDS AND ALLOCATE AND ASSOCIATE DATA BLOCKS !! ! calculate block sizes ! xl = xlen / ddims(1) yl = ylen / ddims(2) zl = zlen / ddims(3) ! fill out block structure fields ! do k = 1, ddims(3) ! claculate the block position along Z ! kc = k - 1 ! calculate the Z bounds ! zmn = zmin + kc * zl zmx = zmin + k * zl do j = 1, ddims(2) ! claculate the block position along Y ! jc = j - 1 ! calculate the Y bounds ! ymn = ymin + jc * yl ymx = ymin + j * yl do i = 1, ddims(1) ! claculate the block position along Y ! ic = i - 1 ! calculate the Z bounds ! xmn = xmin + ic * xl xmx = xmin + i * xl ! assign a pointer ! pmeta => block_array(i,j,k)%ptr ! mark it as the leaf ! call metablock_set_leaf(pmeta) ! set the level ! call metablock_set_level(pmeta, 1) ! set block coordinates ! call metablock_set_coordinates(pmeta, ic, jc, kc) ! set the bounds ! call metablock_set_bounds(pmeta, xmn, xmx, ymn, ymx, zmn, zmx) end do end do end do !! ASSIGN THE BLOCK NEIGHBORS !! ! allocate indices ! allocate(im(ddims(1)), ip(ddims(1))) allocate(jm(ddims(2)), jp(ddims(2))) #if NDIMS == 3 allocate(km(ddims(3)), kp(ddims(3))) #endif /* NDIMS == 3 */ ! generate indices ! im(:) = cshift((/(i, i = 1, ddims(1))/),-1) ip(:) = cshift((/(i, i = 1, ddims(1))/), 1) jm(:) = cshift((/(j, j = 1, ddims(2))/),-1) jp(:) = cshift((/(j, j = 1, ddims(2))/), 1) #if NDIMS == 3 km(:) = cshift((/(k, k = 1, ddims(3))/),-1) kp(:) = cshift((/(k, k = 1, ddims(3))/), 1) #endif /* NDIMS == 3 */ ! check periodicity and reset the edge indices if box is not periodic ! if (.not. periodic(1)) then im( 1) = 0 ip(ddims(1)) = 0 end if if (.not. periodic(2)) then jm( 1) = 0 jp(ddims(2)) = 0 end if #if NDIMS == 3 if (.not. periodic(3)) then km( 1) = 0 kp(ddims(3)) = 0 end if #endif /* NDIMS == 3 */ ! iterate over all initial blocks ! do k = 1, ddims(3) do j = 1, ddims(2) do i = 1, ddims(1) ! assign pmeta with the current block ! pmeta => block_array(i,j,k)%ptr #if NDIMS == 3 ! assign face neighbor pointers ! if (im(i) > 0) then pmeta%faces(1,1,1,1)%ptr => block_array(im(i),j,k)%ptr pmeta%faces(1,2,1,1)%ptr => block_array(im(i),j,k)%ptr pmeta%faces(1,1,2,1)%ptr => block_array(im(i),j,k)%ptr pmeta%faces(1,2,2,1)%ptr => block_array(im(i),j,k)%ptr end if if (ip(i) > 0) then pmeta%faces(2,1,1,1)%ptr => block_array(ip(i),j,k)%ptr pmeta%faces(2,2,1,1)%ptr => block_array(ip(i),j,k)%ptr pmeta%faces(2,1,2,1)%ptr => block_array(ip(i),j,k)%ptr pmeta%faces(2,2,2,1)%ptr => block_array(ip(i),j,k)%ptr end if if (jm(j) > 0) then pmeta%faces(1,1,1,2)%ptr => block_array(i,jm(j),k)%ptr pmeta%faces(2,1,1,2)%ptr => block_array(i,jm(j),k)%ptr pmeta%faces(1,1,2,2)%ptr => block_array(i,jm(j),k)%ptr pmeta%faces(2,1,2,2)%ptr => block_array(i,jm(j),k)%ptr end if if (jp(j) > 0) then pmeta%faces(1,2,1,2)%ptr => block_array(i,jp(j),k)%ptr pmeta%faces(2,2,1,2)%ptr => block_array(i,jp(j),k)%ptr pmeta%faces(1,2,2,2)%ptr => block_array(i,jp(j),k)%ptr pmeta%faces(2,2,2,2)%ptr => block_array(i,jp(j),k)%ptr end if if (km(k) > 0) then pmeta%faces(1,1,1,3)%ptr => block_array(i,j,km(k))%ptr pmeta%faces(2,1,1,3)%ptr => block_array(i,j,km(k))%ptr pmeta%faces(1,2,1,3)%ptr => block_array(i,j,km(k))%ptr pmeta%faces(2,2,1,3)%ptr => block_array(i,j,km(k))%ptr end if if (kp(k) > 0) then pmeta%faces(1,1,2,3)%ptr => block_array(i,j,kp(k))%ptr pmeta%faces(2,1,2,3)%ptr => block_array(i,j,kp(k))%ptr pmeta%faces(1,2,2,3)%ptr => block_array(i,j,kp(k))%ptr pmeta%faces(2,2,2,3)%ptr => block_array(i,j,kp(k))%ptr end if #endif /* NDIMS == 3 */ ! assign edge neighbor pointers ! #if NDIMS == 2 if (im(i) > 0) then pmeta%edges(1,1,2)%ptr => block_array(im(i),j,k)%ptr pmeta%edges(1,2,2)%ptr => block_array(im(i),j,k)%ptr end if if (ip(i) > 0) then pmeta%edges(2,1,2)%ptr => block_array(ip(i),j,k)%ptr pmeta%edges(2,2,2)%ptr => block_array(ip(i),j,k)%ptr end if if (jm(j) > 0) then pmeta%edges(1,1,1)%ptr => block_array(i,jm(j),k)%ptr pmeta%edges(2,1,1)%ptr => block_array(i,jm(j),k)%ptr end if if (jp(j) > 0) then pmeta%edges(1,2,1)%ptr => block_array(i,jp(j),k)%ptr pmeta%edges(2,2,1)%ptr => block_array(i,jp(j),k)%ptr end if #endif /* NDIMS == 2 */ #if NDIMS == 3 if (jm(j) > 0 .and. km(k) > 0) then pmeta%edges(1,1,1,1)%ptr => block_array(i,jm(j),km(k))%ptr pmeta%edges(2,1,1,1)%ptr => block_array(i,jm(j),km(k))%ptr end if if (jp(j) > 0 .and. km(k) > 0) then pmeta%edges(1,2,1,1)%ptr => block_array(i,jp(j),km(k))%ptr pmeta%edges(2,2,1,1)%ptr => block_array(i,jp(j),km(k))%ptr end if if (jm(j) > 0 .and. kp(k) > 0) then pmeta%edges(1,1,2,1)%ptr => block_array(i,jm(j),kp(k))%ptr pmeta%edges(2,1,2,1)%ptr => block_array(i,jm(j),kp(k))%ptr end if if (jp(j) > 0 .and. kp(k) > 0) then pmeta%edges(1,2,2,1)%ptr => block_array(i,jp(j),kp(k))%ptr pmeta%edges(2,2,2,1)%ptr => block_array(i,jp(j),kp(k))%ptr end if if (im(i) > 0 .and. km(k) > 0) then pmeta%edges(1,1,1,2)%ptr => block_array(im(i),j,km(k))%ptr pmeta%edges(1,2,1,2)%ptr => block_array(im(i),j,km(k))%ptr end if if (ip(i) > 0 .and. km(k) > 0) then pmeta%edges(2,1,1,2)%ptr => block_array(ip(i),j,km(k))%ptr pmeta%edges(2,2,1,2)%ptr => block_array(ip(i),j,km(k))%ptr end if if (im(i) > 0 .and. kp(k) > 0) then pmeta%edges(1,1,2,2)%ptr => block_array(im(i),j,kp(k))%ptr pmeta%edges(1,2,2,2)%ptr => block_array(im(i),j,kp(k))%ptr end if if (ip(i) > 0 .and. kp(k) > 0) then pmeta%edges(2,1,2,2)%ptr => block_array(ip(i),j,kp(k))%ptr pmeta%edges(2,2,2,2)%ptr => block_array(ip(i),j,kp(k))%ptr end if if (im(i) > 0 .and. jm(j) > 0) then pmeta%edges(1,1,1,3)%ptr => block_array(im(i),jm(j),k)%ptr pmeta%edges(1,1,2,3)%ptr => block_array(im(i),jm(j),k)%ptr end if if (ip(i) > 0 .and. jm(j) > 0) then pmeta%edges(2,1,1,3)%ptr => block_array(ip(i),jm(j),k)%ptr pmeta%edges(2,1,2,3)%ptr => block_array(ip(i),jm(j),k)%ptr end if if (im(i) > 0 .and. jp(j) > 0) then pmeta%edges(1,2,1,3)%ptr => block_array(im(i),jp(j),k)%ptr pmeta%edges(1,2,2,3)%ptr => block_array(im(i),jp(j),k)%ptr end if if (ip(i) > 0 .and. jp(j) > 0) then pmeta%edges(2,2,1,3)%ptr => block_array(ip(i),jp(j),k)%ptr pmeta%edges(2,2,2,3)%ptr => block_array(ip(i),jp(j),k)%ptr end if #endif /* NDIMS == 3 */ ! assign corner neighbor pointers ! #if NDIMS == 2 if (im(i) > 0 .and. jm(j) > 0) & pmeta%corners(1,1)%ptr => block_array(im(i),jm(j),k)%ptr if (ip(i) > 0 .and. jm(j) > 0) & pmeta%corners(2,1)%ptr => block_array(ip(i),jm(j),k)%ptr if (im(i) > 0 .and. jp(j) > 0) & pmeta%corners(1,2)%ptr => block_array(im(i),jp(j),k)%ptr if (ip(i) > 0 .and. jp(j) > 0) & pmeta%corners(2,2)%ptr => block_array(ip(i),jp(j),k)%ptr #endif /* NDIMS == 2 */ #if NDIMS == 3 if (im(i) > 0 .and. jm(j) > 0 .and. km(k) > 0) & pmeta%corners(1,1,1)%ptr => block_array(im(i),jm(j),km(k))%ptr if (ip(i) > 0 .and. jm(j) > 0 .and. km(k) > 0) & pmeta%corners(2,1,1)%ptr => block_array(ip(i),jm(j),km(k))%ptr if (im(i) > 0 .and. jp(j) > 0 .and. km(k) > 0) & pmeta%corners(1,2,1)%ptr => block_array(im(i),jp(j),km(k))%ptr if (ip(i) > 0 .and. jp(j) > 0 .and. km(k) > 0) & pmeta%corners(2,2,1)%ptr => block_array(ip(i),jp(j),km(k))%ptr if (im(i) > 0 .and. jm(j) > 0 .and. kp(k) > 0) & pmeta%corners(1,1,2)%ptr => block_array(im(i),jm(j),kp(k))%ptr if (ip(i) > 0 .and. jm(j) > 0 .and. kp(k) > 0) & pmeta%corners(2,1,2)%ptr => block_array(ip(i),jm(j),kp(k))%ptr if (im(i) > 0 .and. jp(j) > 0 .and. kp(k) > 0) & pmeta%corners(1,2,2)%ptr => block_array(im(i),jp(j),kp(k))%ptr if (ip(i) > 0 .and. jp(j) > 0 .and. kp(k) > 0) & pmeta%corners(2,2,2)%ptr => block_array(ip(i),jp(j),kp(k))%ptr #endif /* NDIMS == 3 */ end do ! over i end do ! over j end do ! over k ! deallocate indices ! deallocate(im, ip) deallocate(jm, jp) #if NDIMS == 3 deallocate(km, kp) #endif /* NDIMS == 3 */ ! deallocate the block pointer array ! deallocate(block_array) !------------------------------------------------------------------------------- ! end subroutine setup_domain_default !=============================================================================== ! end module domains