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-11-25 10:41:07 -03:00
commit c347f9508f
5 changed files with 277 additions and 272 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
sources/evolution.F90
View File

@ -145,7 +145,10 @@ module evolution
!
subroutine initialize_evolution(verbose, status)
use coordinates, only : toplev, adx, ady, adz
use coordinates, only : toplev, adx, ady
#if NDIMS == 3
use coordinates, only : adz
#endif /* NDIMS == 3 */
use forcing , only : initialize_forcing
use gravity , only : initialize_gravity
use helpers , only : print_message

138
sources/io.F90
View File

@ -1837,9 +1837,7 @@ module io
if (allocated(block_array)) deallocate(block_array)
#ifdef MPI
! redistribute blocks between processors
!
call redistribute_blocks()
call redistribute_blocks(status)
#endif /* MPI */
!-------------------------------------------------------------------------------
@ -3272,11 +3270,11 @@ module io
end if
end if
call redistribute_blocks()
call redistribute_blocks(status)
end do
else
call redistribute_blocks()
call redistribute_blocks(status)
end if
#endif /* MPI */
@ -3615,7 +3613,7 @@ module io
subroutine store_attributes_h5(loc_id, problem, restart, status)
use blocks , only : get_mblocks, get_dblocks, get_nleafs
use blocks , only : get_last_id, nchildren
use blocks , only : get_last_id, nregs, nchildren
use coordinates, only : minlev, maxlev
use coordinates, only : bcells, ncells, nghosts
use coordinates, only : bdims => domain_base_dims
@ -3732,6 +3730,8 @@ module io
H5T_NATIVE_INTEGER, 1, nchildren, status)
call store_attribute_h5(grp_id, 'mblocks', &
H5T_NATIVE_INTEGER, 1, get_mblocks(), status)
call store_attribute_h5(grp_id, 'nregisters', &
H5T_NATIVE_INTEGER, 1, nregs, status)
call store_attribute_h5(grp_id, 'last_id', &
H5T_NATIVE_INTEGER, 1, get_last_id(), status)
call store_attribute_h5(grp_id, 'maximum_rejections', &
@ -4255,7 +4255,8 @@ module io
subroutine restore_datablocks_h5(loc_id, status)
use blocks , only : block_meta, block_data
use blocks , only : append_datablock, link_blocks
use blocks , only : append_datablock, link_blocks, nregs
use coordinates, only : bcells, nghosts
use helpers , only : print_message
implicit none
@ -4267,10 +4268,11 @@ module io
integer(hid_t) :: grp_id, blk_id
character(len=64) :: blk_name
integer(kind=4) :: dblocks, l, id
integer(kind=4) :: dblocks, l, id, nr, nv, nm, ng, nl, nu
integer(hsize_t), dimension(4) :: pdims = 1
integer(hsize_t), dimension(5) :: cdims = 1
integer(hsize_t), dimension(5) :: dims = 1
real(kind=8), dimension(:,:,:,:,:), allocatable :: array
character(len=*), parameter :: loc = 'IO::restore_datablocks_h5()'
@ -4284,6 +4286,14 @@ module io
call restore_attribute_h5(grp_id, 'dblocks', &
H5T_NATIVE_INTEGER, 1, dblocks, status)
call restore_attribute_h5(grp_id, 'nregisters', &
H5T_NATIVE_INTEGER, 1, nr, status)
call restore_attribute_h5(grp_id, 'nvars', &
H5T_NATIVE_INTEGER, 1, nv, status)
call restore_attribute_h5(grp_id, 'bcells', &
H5T_NATIVE_INTEGER, 1, nm, status)
call restore_attribute_h5(grp_id, 'nghosts', &
H5T_NATIVE_INTEGER, 1, ng, status)
call h5gclose_f(grp_id, status)
if (status /= 0) then
@ -4299,39 +4309,89 @@ module io
return
end if
do l = 1, dblocks
write(blk_name, "('datablock_', i0)") l
#if NDIMS == 3
allocate(array(nv,nm,nm,nm,nr), stat=status)
#else /* NDIMS == 3 */
allocate(array(nv,nm,nm, 1,nr), stat=status)
#endif /* NDIMS == 3 */
if (status == 0) then
dims = shape(array)
call append_datablock(pdata, status)
if (status /= 0) then
call print_message(loc, "Could not append new datablock!")
go to 1000
end if
nr = min(nr, nregs)
call h5gopen_f(grp_id, blk_name, blk_id, status)
if (status == 0) then
call restore_attribute_h5(blk_id, 'meta', &
H5T_NATIVE_INTEGER, 1, id, status)
call link_blocks(block_array(id)%ptr, pdata)
pdims = shape(pdata%q)
cdims = shape(pdata%uu)
call read_dataset_h5(blk_id, 'primitive_variables', &
H5T_NATIVE_DOUBLE, pdims, pdata%q , status)
call read_dataset_h5(blk_id, 'conservative_variables', &
H5T_NATIVE_DOUBLE, cdims, pdata%uu, status)
call h5gclose_f(blk_id, status)
if (status /= 0) &
call print_message(loc, &
"Could not close group '" // trim(blk_name) // "'!")
if (ng >= nghosts) then
nl = 1 + (ng - nghosts)
nu = nm - (ng - nghosts)
else
call print_message(loc, &
"Could not open group '" // trim(blk_name) // "'!")
nl = 1 + (nghosts - ng)
nu = bcells - (nghosts - ng)
end if
end do
do l = 1, dblocks
write(blk_name, "('datablock_', i0)") l
call append_datablock(pdata, status)
if (status /= 0) then
call print_message(loc, "Could not append new datablock!")
go to 1000
end if
call h5gopen_f(grp_id, blk_name, blk_id, status)
if (status == 0) then
call restore_attribute_h5(blk_id, 'meta', &
H5T_NATIVE_INTEGER, 1, id, status)
call link_blocks(block_array(id)%ptr, pdata)
call read_dataset_h5(blk_id, 'primitive_variables', &
H5T_NATIVE_DOUBLE, dims(1:4), array(:,:,:,:,1), status)
if (ng >= nghosts) then
#if NDIMS == 3
pdata%q(:,:,:,:) = array(:,nl:nu,nl:nu,nl:nu,1)
#else /* NDIMS == 3 */
pdata%q(:,:,:,:) = array(:,nl:nu,nl:nu, : ,1)
#endif /* NDIMS == 3 */
else
#if NDIMS == 3
pdata%q(:,nl:nu,nl:nu,nl:nu) = array(:,:,:,:,1)
#else /* NDIMS == 3 */
pdata%q(:,nl:nu,nl:nu, : ) = array(:,:,:,:,1)
#endif /* NDIMS == 3 */
end if
call read_dataset_h5(blk_id, 'conservative_variables', &
H5T_NATIVE_DOUBLE, dims(1:5), array(:,:,:,:,:), status)
if (ng >= nghosts) then
#if NDIMS == 3
pdata%uu(:,:,:,:,1:nr) = array(:,nl:nu,nl:nu,nl:nu,1:nr)
#else /* NDIMS == 3 */
pdata%uu(:,:,:,:,1:nr) = array(:,nl:nu,nl:nu, : ,1:nr)
#endif /* NDIMS == 3 */
else
#if NDIMS == 3
pdata%uu(:,nl:nu,nl:nu,nl:nu,1:nr) = array(:,:,:,:,1:nr)
#else /* NDIMS == 3 */
pdata%uu(:,nl:nu,nl:nu, : ,1:nr) = array(:,:,:,:,1:nr)
#endif /* NDIMS == 3 */
end if
call h5gclose_f(blk_id, status)
if (status /= 0) &
call print_message(loc, &
"Could not close group '" // trim(blk_name) // "'!")
else
call print_message(loc, &
"Could not open group '" // trim(blk_name) // "'!")
end if
end do
deallocate(array, stat=status)
if (status /= 0) &
call print_message(loc, &
"Could not deallocate memory for the temporary arrays!")
else
call print_message(loc, &
"Could not allocate memory for the temporary arrays!")
end if
1000 continue

124
sources/mesh.F90
View File

@ -67,7 +67,7 @@ module mesh
public :: initialize_mesh, finalize_mesh, print_mesh
public :: generate_mesh, update_mesh
public :: redistribute_blocks
public :: setup_problem
public :: setup_domain, setup_problem
public :: changed
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
@ -104,28 +104,14 @@ module mesh
logical, intent(in) :: verbose
integer, intent(out) :: status
character(len=64) :: problem = "none"
!-------------------------------------------------------------------------------
!
status = 0
! associate the setup_domain pointer with the respective problem setup
! subroutine
!
select case(trim(problem))
setup_domain => setup_domain_default
case default
setup_domain => setup_domain_default
end select
! prepare the dimensions of the prolongation array
!
pm(1:NDIMS) = 2 * (ncells + nghosts)
! initialize refinement module
!
call initialize_refinement(verbose, status)
!-------------------------------------------------------------------------------
@ -566,9 +552,8 @@ module mesh
if (status /= 0) go to 100
#ifdef MPI
! redistribute blocks equally among all processors
!
call redistribute_blocks()
call redistribute_blocks(status)
if (status /= 0) go to 100
#endif /* MPI */
#ifdef DEBUG
@ -592,24 +577,30 @@ module mesh
!
! Subroutine redistributes data blocks between processes.
!
! Arguments:
!
! status - the call status;
!
!===============================================================================
!
subroutine redistribute_blocks()
subroutine redistribute_blocks(status)
#ifdef MPI
use blocks , only : block_meta, block_data, list_meta
use blocks , only : get_nleafs
use blocks , only : append_datablock, remove_datablock, link_blocks
use mpitools , only : nprocs, npmax, nproc, nodes, lprocs
use mpitools , only : send_array, receive_array
use blocks , only : block_meta, block_data, list_meta
use blocks , only : get_nleafs, get_last_id
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 : send_array, receive_array
#endif /* MPI */
implicit none
integer, intent(out) :: status
#ifdef MPI
logical :: flag
integer :: status
integer(kind=4) :: np, nl, nc, nr
type(block_meta), pointer :: pmeta
@ -619,12 +610,16 @@ module mesh
integer(kind=4), dimension(nodes) :: nb
integer(kind=4), dimension(lprocs,nodes) :: lb
character(len=*), parameter :: loc = 'MESH::redistribute_blocks()'
#endif /* MPI */
!-------------------------------------------------------------------------------
!
status = 0
#ifdef MPI
! calculate the new data block division between processes
! calculate the new data block distribution between processes
!
nl = mod(get_nleafs(), nodes)
nb(:) = get_nleafs() / nodes
@ -635,8 +630,6 @@ module mesh
lb(1:nl,nc) = lb(1:nl,nc) + 1
end do
! reset the processor and leaf numbers
!
nc = 1
nr = 1
np = 0
@ -644,70 +637,54 @@ module mesh
pmeta => list_meta
do while (associated(pmeta))
! consider only meta blocks which belong to active processes
!
if (pmeta%process < nprocs) then
! check if the block belongs to another process
! exchange data blocks between processes
!
if (pmeta%process /= np) then
if (pmeta%leaf) then
! indicate that the block structure will change
!
changed = .true.
! generate a tag for communication
!
tag1 = np * nprocs + pmeta%process
tag2 = tag1 + 1
tag1 = pmeta%id
tag2 = pmeta%id + get_last_id()
! send the data block arrays to the destination process, and
! remove the data block from the current process
!
if (nproc == pmeta%process) then
call send_array(np, tag1, pmeta%data%uu(:,:,:,:,:))
call send_array(np, tag2, pmeta%data%q (:,:,:,:))
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
end if ! nproc == pmeta%process
end if
! on the destination process, allocate new data block, associate it with
! the corresponding meta block, and fill it with the received arrays
!
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 (:,:,:,:))
call receive_array(pmeta%process, tag1, pmeta%data%uu)
call receive_array(pmeta%process, tag2, pmeta%data%q)
end if ! nproc == n
end if ! leaf
end if
end if
! set new processor number
!
pmeta%process = np
end if ! pmeta%process /= np
end if
! increase the number of blocks on the current process; if it exceeds the
! allowed number, reset the counter and increase the processor number
! determine the new block distribution
!
if (pmeta%leaf) then
! increase the number of leafs for the current process
!
nl = nl + 1
! if the number of leafs for the current process exceeds the number of assigned
! blocks, reset the counter and increase the process number
!
if (nl >= lb(nr,nc)) then
np = min(npmax, np + 1)
@ -730,14 +707,17 @@ module mesh
if (flag) flag = lb(nr,nc) == 0
end do
end if ! nl >= lb(nr,np)
end if
end if ! leaf
end if ! pmeta%process < nprocs
end if
end if
pmeta => pmeta%next
end do ! pmeta
end do
1000 continue
if (check_status(status /= 0)) &
call print_message(loc, "Could not redistribute blocks!")
#endif /* MPI */
!-------------------------------------------------------------------------------
@ -1898,7 +1878,11 @@ module mesh
! redistribute blocks equally among all processors
!
call redistribute_blocks()
call redistribute_blocks(status)
if (status /= 0) then
call print_message(loc, "Cannot redistribute blocks!")
go to 100
end if
#endif /* MPI */
end if ! leaf at current level selected for refinement

249
sources/refinement.F90
View File

@ -33,30 +33,35 @@ module refinement
implicit none
! refinement criterion parameters
!
abstract interface
real(kind=4) function user_criterion_iface(pdata) result(crit)
use blocks, only : block_data
type(block_data), pointer, intent(in) :: pdata
end function
end interface
procedure(user_criterion_iface), pointer, save :: user_criterion => null()
real(kind=8), save :: crefmin = 2.0d-01
real(kind=8), save :: crefmax = 8.0d-01
real(kind=8), save :: vortmin = 1.0d-03
real(kind=8), save :: vortmax = 1.0d-01
real(kind=8), save :: currmin = 1.0d-03
real(kind=8), save :: currmax = 1.0d-01
real(kind=8), save :: usermin = 1.0d+99
real(kind=8), save :: usermax = 1.0d+99
real(kind=8), save :: epsref = 1.0d-02
! flags for variable included in the refinement criterion calculation
!
logical, dimension(:), allocatable, save :: qvar_ref
logical , save :: vort_ref = .false.
logical , save :: curr_ref = .false.
logical , save :: user_ref = .false.
! by default everything is private
!
private
! declare public subroutines
!
public :: initialize_refinement, finalize_refinement, print_refinement
public :: check_refinement_criterion
public :: user_criterion
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
@ -85,6 +90,7 @@ module refinement
subroutine initialize_refinement(verbose, status)
use equations , only : magnetized, nv, pvars
use helpers , only : print_message
use parameters, only : get_parameter
implicit none
@ -96,104 +102,94 @@ module refinement
integer :: p
character(len=255) :: variables = "dens pres"
character(len=*), parameter :: loc = 'REFINEMENT:initialize_refinement()'
!-------------------------------------------------------------------------------
!
status = 0
! get the refinement parameters
!
call get_parameter("crefmin", crefmin)
call get_parameter("crefmax", crefmax)
call get_parameter("vortmin", vortmin)
call get_parameter("vortmax", vortmax)
call get_parameter("currmin", currmin)
call get_parameter("currmax", currmax)
call get_parameter("usermin", usermin)
call get_parameter("usermax", usermax)
call get_parameter("epsref" , epsref )
! get variables to include in the refinement criterion calculation
!
call get_parameter("refinement_variables", variables)
! allocate vector for indicators, which variables are taken into account in
! calculating the refinement criterion
!
allocate(qvar_ref(nv))
allocate(qvar_ref(nv), stat=status)
if (status /= 0) &
call print_message(loc, &
"Could not allocate the space for refinement criterion names!")
! check which primitive variable is used to determine the refinement criterion
!
do p = 1, nv
qvar_ref(p) = index(variables, trim(pvars(p))) > 0
end do ! p = 1, nv
end do
! turn on refinement based on vorticity if specified
!
vort_ref = index(variables, 'vort') > 0
! check if any resonable variable controls the refinement
!
test = any(qvar_ref(:)) .or. vort_ref
! turn on refinement based on current density if specified
!
if (magnetized) then
curr_ref = index(variables, 'curr') > 0 .or. index(variables, 'jabs') > 0
test = test .or. curr_ref
end if
! check if there is any variable selected to control refinement
!
user_ref = index(variables, 'user') > 0
test = test .or. user_ref
if (.not. test) then
if (verbose) then
write(*,*)
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "No or wrong variable selected to control" // &
" the refinement."
write(*,"(1x,a)") "Available control variables: any primitive" // &
" variable, or 'vort' for the magnitude of" // &
" vorticity, or 'curr' for the magnitude of" // &
" current density."
end if
if (verbose) &
call print_message(loc, "No refinement criterion has been selected!")
status = 1
return
end if
if (crefmin > crefmax .or. crefmin < 0.0d+00) then
if (verbose) then
write(*,*)
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "Wrong 'crefmin' or 'crefmax' parameters."
write(*,"(1x,a)") "Both should be positive and 'crefmin' <= 'crefmax'."
end if
status = 1
end if
if (epsref < 0.0d+00) then
if (verbose) then
write(*,*)
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "Wrong 'epsref' parameters."
write(*,"(1x,a)") "It should be positive."
end if
if (verbose) &
call print_message(loc, &
"Wrong 'crefmin' or 'crefmax' parameters. " // &
"They should be positive and 'crefmin' <= 'crefmax'.")
status = 1
return
end if
if (vortmin > vortmax .or. vortmin < 0.0d+00) then
if (verbose) then
write(*,*)
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "Wrong 'vortmin' or 'vortmax' parameters."
write(*,"(1x,a)") "Both should be positive and 'vortmin' <= 'vortmax'."
end if
if (verbose) &
call print_message(loc, &
"Wrong 'vortmin' or 'vortmax' parameters. " // &
"They should be positive and 'vortmin' <= 'vortmax'.")
status = 1
return
end if
if ((currmin > currmax .or. currmin < 0.0d+00) .and. magnetized) then
if (verbose) then
write(*,*)
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "Wrong 'currmin' or 'currmax' parameters."
write(*,"(1x,a)") "Both should be positive and 'currmin' <= 'currmax'."
end if
if (verbose) &
call print_message(loc, &
"Wrong 'currmin' or 'currmax' parameters. " // &
"They should be positive and 'currmin' <= 'currmax'.")
status = 1
return
end if
if ((usermin > usermax .or. usermin < 0.0d+00) .and. user_ref) then
if (verbose) &
call print_message(loc, &
"Wrong 'usermin' or 'usermax' parameters. " // &
"They should be positive and 'usermin' <= 'usermax'.")
status = 1
return
end if
if (epsref <= 0.0d+00) then
if (verbose) &
call print_message(loc, &
"Wrong 'epsref' parameters. It should be positive.")
status = 1
return
end if
!-------------------------------------------------------------------------------
@ -215,17 +211,24 @@ module refinement
!
subroutine finalize_refinement(status)
use helpers, only : print_message
implicit none
integer, intent(out) :: status
character(len=*), parameter :: loc = 'REFINEMENT:dinalize_refinement()'
!-------------------------------------------------------------------------------
!
status = 0
! deallocate refined variable indicators
!
if (allocated(qvar_ref)) deallocate(qvar_ref)
if (allocated(qvar_ref)) then
deallocate(qvar_ref, stat=status)
if (status /= 0) &
call print_message(loc, &
"Could not deallocate space for the refinement criterion names!")
end if
!-------------------------------------------------------------------------------
!
@ -258,30 +261,31 @@ module refinement
!-------------------------------------------------------------------------------
!
if (verbose) then
rvars = ""
do p = 1, nv
if (qvar_ref(p)) rvars = adjustl(trim(rvars) // ' ' // trim(pvars(p)))
end do
if (vort_ref) then
rvars = adjustl(trim(rvars) // ' vort')
end if
if (magnetized .and. curr_ref) then
rvars = adjustl(trim(rvars) // ' curr')
end if
call print_section(verbose, "Refinement")
call print_parameter(verbose, "refined variables", rvars)
call print_parameter(verbose, "2nd order error limits", crefmin, crefmax)
if (vort_ref) then
call print_parameter(verbose, "vorticity limits", vortmin, vortmax)
end if
if (magnetized .and. curr_ref) &
call print_parameter(verbose, "current density limits", &
currmin, currmax)
if (.not. verbose) return
rvars = ""
do p = 1, nv
if (qvar_ref(p)) rvars = adjustl(trim(rvars) // ' ' // trim(pvars(p)))
end do
if (vort_ref) then
rvars = adjustl(trim(rvars) // ' vort')
end if
if (magnetized .and. curr_ref) then
rvars = adjustl(trim(rvars) // ' curr')
end if
if (user_ref) then
rvars = adjustl(trim(rvars) // ' user')
end if
call print_section(verbose, "Refinement")
call print_parameter(verbose, "refined variables", rvars)
call print_parameter(verbose, "2nd order error limits", crefmin, crefmax)
if (vort_ref) &
call print_parameter(verbose, "vorticity limits", vortmin, vortmax)
if (magnetized .and. curr_ref) &
call print_parameter(verbose, "current density limits", currmin, currmax)
if (user_ref) &
call print_parameter(verbose, "user criterion limits", usermin, usermax)
!-------------------------------------------------------------------------------
!
@ -313,7 +317,7 @@ module refinement
type(block_data), pointer, intent(in) :: pdata
integer(kind=4) :: criterion
integer(kind=4) :: criterion
integer :: p
real(kind=8) :: cref
@ -322,8 +326,6 @@ module refinement
!
criterion = -1
! check the second derivative error from selected primitive variables
!
do p = 1, nv
if (qvar_ref(p)) then
@ -333,10 +335,8 @@ module refinement
if (cref > crefmax) criterion = max(criterion, 1)
end if
end do ! p = 1, nv
end do
! check vorticity criterion
!
if (vort_ref) then
cref = vorticity_magnitude(pdata)
@ -344,8 +344,6 @@ module refinement
if (cref > vortmax) criterion = max(criterion, 1)
end if
! check current density criterion
!
if (curr_ref) then
cref = current_density_magnitude(pdata)
@ -353,6 +351,13 @@ module refinement
if (cref > currmax) criterion = max(criterion, 1)
end if
if (user_ref .and. associated(user_criterion)) then
cref = user_criterion(pdata)
if (cref > usermin) criterion = max(criterion, 0)
if (cref > usermax) criterion = max(criterion, 1)
end if
return
!-------------------------------------------------------------------------------
@ -409,12 +414,8 @@ module refinement
!
error = 0.0e+00
! calculate local refinement criterion for variable which exists
!
if (iqt > 0) then
! find the maximum smoothness indicator over all cells
!
#if NDIMS == 3
do k = nbl, neu
km1 = k - 1
@ -427,34 +428,26 @@ module refinement
im1 = i - 1
ip1 = i + 1
! calculate the second derivative error the X direction
!
fr = pdata%q(iqt,ip1,j,k) - pdata%q(iqt,i ,j,k)
fl = pdata%q(iqt,im1,j,k) - pdata%q(iqt,i ,j,k)
fc = abs(pdata%q(iqt,ip1,j,k)) + abs(pdata%q(iqt,im1,j,k)) &
+ 2.0d+00 * abs(pdata%q(iqt,i,j,k))
fc = abs(pdata%q(iqt,ip1,j,k)) + abs(pdata%q(iqt,im1,j,k)) &
+ 2.0d+00 * abs(pdata%q(iqt,i,j,k))
fx = abs(fr + fl) / (abs(fr) + abs(fl) + epsref * fc + eps)
! calculate the second derivative error along the Y direction
!
fr = pdata%q(iqt,i,jp1,k) - pdata%q(iqt,i,j ,k)
fl = pdata%q(iqt,i,jm1,k) - pdata%q(iqt,i,j ,k)
fc = abs(pdata%q(iqt,i,jp1,k)) + abs(pdata%q(iqt,i,jm1,k)) &
+ 2.0d+00 * abs(pdata%q(iqt,i,j,k))
fc = abs(pdata%q(iqt,i,jp1,k)) + abs(pdata%q(iqt,i,jm1,k)) &
+ 2.0d+00 * abs(pdata%q(iqt,i,j,k))
fy = abs(fr + fl) / (abs(fr) + abs(fl) + epsref * fc + eps)
#if NDIMS == 3
! calculate the second derivative error along the Z direction
!
fr = pdata%q(iqt,i,j,kp1) - pdata%q(iqt,i,j,k )
fl = pdata%q(iqt,i,j,km1) - pdata%q(iqt,i,j,k )
fc = abs(pdata%q(iqt,i,j,kp1)) + abs(pdata%q(iqt,i,j,km1)) &
+ 2.0d+00 * abs(pdata%q(iqt,i,j,k))
fc = abs(pdata%q(iqt,i,j,kp1)) + abs(pdata%q(iqt,i,j,km1)) &
+ 2.0d+00 * abs(pdata%q(iqt,i,j,k))
fz = abs(fr + fl) / (abs(fr) + abs(fl) + epsref * fc + eps)
#endif /* NDIMS == 3 */
! take the maximum second derivative error
!
#if NDIMS == 2
error = max(error, fx, fy)
#endif /* NDIMS == 2 */
@ -462,16 +455,14 @@ module refinement
error = max(error, fx, fy, fz)
#endif /* NDIMS == 3 */
end do ! i = nbl, neu
end do ! j = nbl, neu
end do
end do
#if NDIMS == 3
end do ! k = nbl, neu
end do
#endif /* NDIMS == 3 */
end if ! iqt > 0
end if
! return the refinement flag
!
return
!-------------------------------------------------------------------------------
@ -544,8 +535,8 @@ module refinement
end if
if (work_in_use) &
call print_message(loc, "Workspace is being used right now! " // &
"Corruptions can occur!")
call print_message(loc, &
"Workspace is being used right now! Corruptions can occur!")
work_in_use = .true.
@ -644,8 +635,8 @@ module refinement
end if
if (work_in_use) &
call print_message(loc, "Workspace is being used right now! " // &
"Corruptions can occur!")
call print_message(loc, &
"Workspace is being used right now! Corruptions can occur!")
work_in_use = .true.
@ -661,10 +652,10 @@ module refinement
jmax = max(jmax, real(jabs, kind=4))
end do ! i = nbl, neu
end do ! j = nbl, neu
end do
end do
#if NDIMS == 3
end do ! k = nbl, neu
end do
#endif /* NDIMS == 3 */
work_in_use = .false.

33
sources/user_problem.F90
View File

@ -1053,39 +1053,6 @@ module user_problem
!
!===============================================================================
!
! subroutine USER_GRAVITATIONAL_ACCELERATION:
! ------------------------------------------
!
! Subroutine returns the user defined gravitational acceleration.
!
! Arguments:
!
! t, dt - time and the time increment;
! x, y, z - rectangular coordinates;
! acc - vector of the gravitational acceleration;
!
!===============================================================================
!
subroutine user_gravitational_acceleration(t, dt, x, y, z, acc)
use parameters, only : get_parameter
implicit none
real(kind=8) , intent(in) :: t, dt
real(kind=8) , intent(in) :: x, y, z
real(kind=8), dimension(3), intent(out) :: acc
!-------------------------------------------------------------------------------
!
acc(:) = 0.0d+00
!-------------------------------------------------------------------------------
!
end subroutine user_gravitational_acceleration
!
!===============================================================================
!
! subroutine USER_BOUNDARY_X:
! --------------------------
!