gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

USER_PROBLEM: Rewrite slightly, remote useless comments.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2021-11-24 20:20:10 -03:00
parent 8482edee56
commit 7ae74da6d6
1 changed files with 266 additions and 438 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

312
sources/user_problem.F90
View File

@ -31,14 +31,6 @@ module user_problem
implicit none implicit none
#ifdef PROFILE
! timer indices
!
integer, save :: imi, imp, ims, imu, img, imb
#endif /* PROFILE */
! default problem parameter values
!
real(kind=8), save :: beta = 1.00d+00 real(kind=8), save :: beta = 1.00d+00
real(kind=8), save :: zeta = 0.00d+00 real(kind=8), save :: zeta = 0.00d+00
real(kind=8), save :: eta = 0.00d+00 real(kind=8), save :: eta = 0.00d+00
@ -80,16 +72,8 @@ module user_problem
logical, save :: absorption = .false. logical, save :: absorption = .false.
logical, save :: resistive = .false. logical, save :: resistive = .false.
! flag indicating if the gravitational source term is enabled
!
logical, save :: user_gravity_enabled = .false.
! allocatable arrays for velocity perturbation
!
real(kind=8), dimension(:), allocatable :: kx, ky, kz, ux, uy, uz, ph real(kind=8), dimension(:), allocatable :: kx, ky, kz, ux, uy, uz, ph
! export subroutines
!
private private
public :: initialize_user_problem, finalize_user_problem public :: initialize_user_problem, finalize_user_problem
public :: setup_user_problem, update_user_sources, update_user_shapes public :: setup_user_problem, update_user_sources, update_user_shapes
@ -119,15 +103,13 @@ module user_problem
! !
subroutine initialize_user_problem(problem, rcount, verbose, status) subroutine initialize_user_problem(problem, rcount, verbose, status)
! include external procedures and variables
!
#if NDIMS == 3 #if NDIMS == 3
use constants , only : pi use constants , only : pi
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
use constants , only : pi2 use constants , only : pi2
use coordinates, only : ng => nghosts, ady, xlen, zlen, ymax use coordinates, only : ng => nghosts, ady, xlen, zlen, ymax
use equations , only : magnetized, csnd, csnd2 use equations , only : magnetized, csnd, csnd2
use helpers , only : print_section, print_parameter use helpers , only : print_section, print_parameter, print_message
use parameters , only : get_parameter use parameters , only : get_parameter
use random , only : randuni, randsym use random , only : randuni, randsym
@ -138,103 +120,79 @@ module user_problem
logical , intent(in) :: verbose logical , intent(in) :: verbose
integer , intent(out) :: status integer , intent(out) :: status
! local variables
!
character(len=64) :: perturbation = "noise", append = "off", fname character(len=64) :: perturbation = "noise", append = "off", fname
character(len=64) :: enable_absorption = "off" character(len=64) :: enable_absorption = "off"
logical :: flag logical :: flag
integer :: n, kd integer :: n, kd
real(kind=8) :: thh, fc, ka, ydis = 9.00d+99 real(kind=8) :: thh, fc, ydis = 9.00d+99
#if NDIMS == 3 #if NDIMS == 3
real(kind=8) :: thv, tx, ty, tz, tt real(kind=8) :: thv, tx, ty, tz, tt
#else /* NDIMS == 3 */
real(kind=8) :: ka
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
character(len=*), parameter :: &
loc = 'USER_PROBLEM::initialize_user_problem()'
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
!
#ifdef PROFILE
! set timer descriptions
!
call set_timer('user_problem:: initialize' , imi)
call set_timer('user_problem:: problem setup', imp)
call set_timer('user_problem:: shape' , ims)
call set_timer('user_problem:: sources' , imu)
call set_timer('user_problem:: gravity' , img)
call set_timer('user_problem:: boundaries' , imb)
! start accounting time for module initialization/finalization
!
call start_timer(imi)
#endif /* PROFILE */
! reset the status flag
! !
status = 0 status = 0
! this problem does not work with not magnetized set of equations
!
if (.not. magnetized) then if (.not. magnetized) then
if (verbose) then if (verbose) &
write(*,*) call print_message(loc, &
write(*,"(1x,a)") "ERROR!" "The reconnection problem does not work without magnetic field.")
write(*,"(1x,a)") "The problem " // trim(problem) // &
" requires magnetized set of equations:" // &
" 'MHD', or 'SR-MHD'."
write(*,*)
end if
status = 1 status = 1
return
end if end if
! proceed if no errors
!
if (status == 0) then
! get the reconnection equilibrium parameters ! get the reconnection equilibrium parameters
! !
call get_parameter("beta", beta) call get_parameter("beta", beta)
if (beta <= 0.0d+00) then if (beta <= 0.0d+00) then
if (verbose) then if (verbose) &
write(*,*) call print_message(loc, "Plasma-beta must be positive (beta > 0.0)!")
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "Parameter 'beta' must be positive (beta > 0.0)!"
write(*,*)
end if
status = 1 status = 1
return
end if end if
call get_parameter("zeta", zeta) call get_parameter("zeta", zeta)
if (zeta < 0.0d+00 .or. zeta > 1.0d+00) then if (zeta < 0.0d+00 .or. zeta > 1.0d+00) then
if (verbose) then if (verbose) &
write(*,*) call print_message(loc, "Parameter 'zeta' must be between 0.0 and 1.0!")
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "Parameter 'zeta' must be between 0.0 and 1.0!"
write(*,*)
end if
status = 1 status = 1
return
end if end if
call get_parameter("resistivity", eta) call get_parameter("resistivity", eta)
if (eta < 0.0d+00) then if (eta < 0.0d+00) then
if (verbose) then if (verbose) &
write(*,*) call print_message(loc, "Resistivity cannot be negative!")
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "Resistivity cannot be negative!"
write(*,*)
end if
status = 1 status = 1
return
else else
resistive = .true. resistive = .true.
end if end if
call get_parameter("dens", dens) call get_parameter("dens", dens)
if (dens <= 0.0d+00) then if (dens <= 0.0d+00) then
if (verbose) then if (verbose) &
write(*,*) call print_message(loc, "Density must be positive (dens > 0.0)!")
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "Parameter 'dens' must be positive (dens > 0.0)!"
write(*,*)
end if
status = 1 status = 1
return
end if end if
call get_parameter("bamp", bamp) call get_parameter("bamp", bamp)
call get_parameter("bgui", bgui) call get_parameter("bgui", bgui)
! get the geometry parameters
!
call get_parameter("delta", dlta)
if (dlta < 0.0d+00) then
if (verbose) &
call print_message(loc, "Density must be positive (dens > 0.0)!")
status = 1
return
end if
call get_parameter("blimit", blim)
! calculate the maximum magnetic pressure, thermal pressure from the plasma-β ! calculate the maximum magnetic pressure, thermal pressure from the plasma-β
! parameters, and the sound speed in the case of isothermal equations of state ! parameters, and the sound speed in the case of isothermal equations of state
! !
@ -246,32 +204,8 @@ module user_problem
valf = sqrt(2.0d+00 * pmag / dens) valf = sqrt(2.0d+00 * pmag / dens)
lund = valf / max(tiny(eta), eta) lund = valf / max(tiny(eta), eta)
dlta = lund**(- 1.0d+00 / 3.0d+00) dlta = lund**(- 1.0d+00 / 3.0d+00)
! get the geometry parameters
!
call get_parameter("delta", dlta)
if (dlta < 0.0d+00) then
if (verbose) then
write(*,*)
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "Parameter 'delta' must be equal or bigger than zero!"
write(*,*)
end if
status = 1
end if
call get_parameter("blimit", blim)
! lower limit for blim
!
blim = max(blim, ng * ady(1)) blim = max(blim, ng * ady(1))
end if ! status
! proceed if no errors
!
if (status == 0) then
! get the perturbation parameters ! get the perturbation parameters
! !
call get_parameter("perturbation", perturbation) call get_parameter("perturbation", perturbation)
@ -354,13 +288,14 @@ module user_problem
ph(n) = pi2 * randuni() ph(n) = pi2 * randuni()
end do end do
end if ! status
end if
end if ! status
! prepare file to store reconnection rate terms else
! if (verbose) &
if (status == 0) then call print_message(loc, &
"Could not allocate space for perturbation vectors!")
return
end if
end if ! pert = 2
! determine if to append or create another file ! determine if to append or create another file
! !
@ -429,12 +364,6 @@ module user_problem
call get_parameter("adec", adec) call get_parameter("adec", adec)
yabs = ymax - abs(acut) yabs = ymax - abs(acut)
end if ! status
! proceed if no errors
!
if (status == 0) then
! print information about the user problem setup ! print information about the user problem setup
! !
call print_section(verbose, "Parameters (* - set, otherwise calculated)") call print_section(verbose, "Parameters (* - set, otherwise calculated)")
@ -474,7 +403,6 @@ module user_problem
call print_parameter(verbose, '(*) adec' , adec) call print_parameter(verbose, '(*) adec' , adec)
call print_parameter(verbose, '( ) yabs' , yabs) call print_parameter(verbose, '( ) yabs' , yabs)
end if end if
end if ! status
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -495,29 +423,23 @@ module user_problem
! !
subroutine finalize_user_problem(status) subroutine finalize_user_problem(status)
use helpers, only : print_message
implicit none implicit none
integer, intent(out) :: status integer, intent(out) :: status
character(len=*), parameter :: loc = 'USER_PROBLEM::finalize_user_problem()'
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
status = 0 status = 0
! close the reconnection file
!
close(runit) close(runit)
! deallocate wave vector components, random directions, and random phase deallocate(kx, ky, kz, ux, uy, uz, ph, stat = status)
! call print_message(loc, &
if (allocated(kx)) deallocate(kx, ky, kz, stat = status) "Could not deallocate space for perturbation vectors!")
if (allocated(ux)) deallocate(ux, uy, uz, stat = status)
if (allocated(ph)) deallocate(ph, stat = status)
#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
call stop_timer(imi)
#endif /* PROFILE */
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -550,7 +472,7 @@ module user_problem
use equations , only : nv, ns use equations , only : nv, ns
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl
use equations , only : csnd2 use equations , only : csnd2
use iso_fortran_env, only : error_unit use helpers , only : print_message
use operators , only : curl use operators , only : curl
use random , only : randsym use random , only : randsym
use workspace , only : resize_workspace, work, work_in_use use workspace , only : resize_workspace, work, work_in_use
@ -589,19 +511,12 @@ module user_problem
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
#ifdef PROFILE
! start accounting time for the problem setup
!
call start_timer(imp)
#endif /* PROFILE */
if (first) then if (first) then
n = (nv + 3) * nn**NDIMS n = (nv + 3) * nn**NDIMS
call resize_workspace(n, status) call resize_workspace(n, status)
if (status /= 0) then if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), & call print_message(loc, "Could not resize the workspace!")
"Could not resize the workspace!"
go to 100 go to 100
end if end if
@ -617,8 +532,6 @@ module user_problem
first = .false. first = .false.
end if end if
! prepare cell sizes and block coordinates
!
dh(1) = adx(pdata%meta%level) dh(1) = adx(pdata%meta%level)
dh(2) = ady(pdata%meta%level) dh(2) = ady(pdata%meta%level)
#if NDIMS == 3 #if NDIMS == 3
@ -675,20 +588,14 @@ module user_problem
fy(j) = cos(yp)**2 fy(j) = cos(yp)**2
end do ! i = 1, nn end do ! i = 1, nn
if (work_in_use) then if (work_in_use) &
write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc), & call print_message(loc, "Workspace is being used right now! " // &
"Workspace is being used right now! " // & "Corruptions can occur!")
"Corruptions can occur!"
end if
work_in_use = .true. work_in_use = .true.
! reset the perturbation matrix
!
qpert(:,:,:,:) = 0.0d+00 qpert(:,:,:,:) = 0.0d+00
! the random perturbation
!
if (pert == 0) then if (pert == 0) then
! of velocity ! of velocity
@ -976,10 +883,6 @@ module user_problem
100 continue 100 continue
#ifdef PROFILE
call stop_timer(imp)
#endif /* PROFILE */
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end subroutine setup_user_problem end subroutine setup_user_problem
@ -1037,11 +940,11 @@ module user_problem
use blocks , only : block_data use blocks , only : block_data
use constants , only : pi use constants , only : pi
use coordinates, only : nn => bcells use coordinates, only : nn => bcells
use coordinates , only : ymax, ay, adx, ady, adz use coordinates, only : ay, adx, ady, adz
use equations , only : nv use equations , only : nv
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
use equations , only : prim2cons use equations , only : prim2cons
use iso_fortran_env, only : error_unit use helpers , only : print_message
use operators , only : laplace use operators , only : laplace
use workspace , only : resize_workspace, work, work_in_use use workspace , only : resize_workspace, work, work_in_use
@ -1065,17 +968,12 @@ module user_problem
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
#ifdef PROFILE
call start_timer(ims)
#endif /* PROFILE */
if (first) then if (first) then
j = nn**NDIMS j = nn**NDIMS
call resize_workspace(j, status) call resize_workspace(j, status)
if (status /= 0) then if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), & call print_message(loc, "Could not resize the workspace!")
"Could not resize the workspace!"
go to 100 go to 100
end if end if
@ -1088,11 +986,9 @@ module user_problem
first = .false. first = .false.
end if end if
if (work_in_use) then if (work_in_use) &
write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc), & call print_message(loc, "Workspace is being used right now! " // &
"Workspace is being used right now! " // & "Corruptions can occur!")
"Corruptions can occur!"
end if
work_in_use = .true. work_in_use = .true.
@ -1151,10 +1047,6 @@ module user_problem
100 continue 100 continue
#ifdef PROFILE
call stop_timer(ims)
#endif /* PROFILE */
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end subroutine update_user_shapes end subroutine update_user_shapes
@ -1213,16 +1105,12 @@ module user_problem
! !
subroutine user_boundary_x(ic, jl, ju, kl, ku, t, dt, x, y, z, qn) subroutine user_boundary_x(ic, jl, ju, kl, ku, t, dt, x, y, z, qn)
! import external procedures and variables
!
use coordinates , only : nn => bcells, nb, ne, nbl, neu use coordinates , only : nn => bcells, nb, ne, nbl, neu
use equations , only : ivx, ibx, iby, ibp use equations , only : magnetized, ivx, ibx, iby, ibp
#if NDIMS == 3 #if NDIMS == 3
use equations , only : ibz use equations , only : ibz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
! local variables are not implicit by default
!
implicit none implicit none
integer , intent(in) :: ic integer , intent(in) :: ic
@ -1234,8 +1122,6 @@ module user_problem
real(kind=8), dimension(:) , intent(in) :: z real(kind=8), dimension(:) , intent(in) :: z
real(kind=8), dimension(:,:,:,:), intent(inout) :: qn real(kind=8), dimension(:,:,:,:), intent(inout) :: qn
! local variables
!
integer :: im2, im1, i , ip1, ip2 integer :: im2, im1, i , ip1, ip2
integer :: jm2, jm1, j , jp1, jp2 integer :: jm2, jm1, j , jp1, jp2
integer :: k = 1 integer :: k = 1
@ -1249,18 +1135,8 @@ module user_problem
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
#ifdef PROFILE if (magnetized) then
! start accounting time for the boundary update
!
call start_timer(imb)
#endif /* PROFILE */
! process case with magnetic field, otherwise revert to standard outflow
!
if (ibx > 0) then
! get the cell sizes and their ratios
!
dx = x(2) - x(1) dx = x(2) - x(1)
dy = y(2) - y(1) dy = y(2) - y(1)
#if NDIMS == 3 #if NDIMS == 3
@ -1271,8 +1147,6 @@ module user_problem
dxz = dx / dz dxz = dx / dz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
! process left and right side boundary separatelly
!
if (ic == 1) then if (ic == 1) then
! iterate over left-side ghost layers ! iterate over left-side ghost layers
@ -1393,13 +1267,7 @@ module user_problem
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
end do ! i = neu, nn end do ! i = neu, nn
end if end if
end if ! ibx > 0 end if ! magnetized
#ifdef PROFILE
! stop accounting time for the boundary update
!
call stop_timer(imb)
#endif /* PROFILE */
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -1426,14 +1294,10 @@ module user_problem
! !
subroutine user_boundary_y(jc, il, iu, kl, ku, t, dt, x, y, z, qn) subroutine user_boundary_y(jc, il, iu, kl, ku, t, dt, x, y, z, qn)
! import external procedures and variables
!
use coordinates, only : nn => bcells, nb, ne, nbl, neu use coordinates, only : nn => bcells, nb, ne, nbl, neu
use equations , only : nv use equations , only : magnetized, nv
use equations , only : idn, ivy, ipr, ibx, iby, ibz, ibp use equations , only : idn, ivy, ipr, ibx, iby, ibz, ibp
! local variables are not implicit by default
!
implicit none implicit none
integer , intent(in) :: jc integer , intent(in) :: jc
@ -1445,8 +1309,6 @@ module user_problem
real(kind=8), dimension(:) , intent(in) :: z real(kind=8), dimension(:) , intent(in) :: z
real(kind=8), dimension(:,:,:,:), intent(inout) :: qn real(kind=8), dimension(:,:,:,:), intent(inout) :: qn
! local variables
!
integer :: i, im1, ip1 integer :: i, im1, ip1
integer :: j, jm1, jp1, jm2, jp2 integer :: j, jm1, jp1, jm2, jp2
integer :: k = 1 integer :: k = 1
@ -1461,18 +1323,8 @@ module user_problem
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
#ifdef PROFILE if (magnetized) then
! start accounting time for the boundary update
!
call start_timer(imb)
#endif /* PROFILE */
! process case with magnetic field, otherwise revert to standard outflow
!
if (ibx > 0) then
! get the cell sizes and their ratios
!
dx = x(2) - x(1) dx = x(2) - x(1)
dy = y(2) - y(1) dy = y(2) - y(1)
#if NDIMS == 3 #if NDIMS == 3
@ -1620,12 +1472,6 @@ module user_problem
end if end if
end if ! ibx > 0 end if ! ibx > 0
#ifdef PROFILE
! stop accounting time for the boundary update
!
call stop_timer(imb)
#endif /* PROFILE */
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end subroutine user_boundary_y end subroutine user_boundary_y
@ -1690,8 +1536,7 @@ module user_problem
use coordinates, only : nn => bcells, nb, ne use coordinates, only : nn => bcells, nb, ne
use coordinates, only : adx, ady, adz, advol, ay, yarea use coordinates, only : adx, ady, adz, advol, ay, yarea
use equations , only : ivx, ivy, ivz, ibx, iby, ibz, ibp use equations , only : ivx, ivy, ivz, ibx, iby, ibz, ibp
use helpers , only : flush_and_sync use helpers , only : print_message, flush_and_sync
use iso_fortran_env, only : error_unit
#ifdef MPI #ifdef MPI
use mpitools , only : reduce_sum use mpitools , only : reduce_sum
#endif /* MPI */ #endif /* MPI */
@ -1719,17 +1564,12 @@ module user_problem
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
#ifdef PROFILE
call start_timer(imt)
#endif /* PROFILE */
if (first) then if (first) then
ni = 3 * nn**NDIMS ni = 3 * nn**NDIMS
call resize_workspace(ni, status) call resize_workspace(ni, status)
if (status /= 0) then if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), & call print_message(loc, "Could not resize the workspace!")
"Could not resize the workspace!"
go to 100 go to 100
end if end if
@ -1744,11 +1584,9 @@ module user_problem
rterms(:) = 0.0d+00 rterms(:) = 0.0d+00
if (work_in_use) then if (work_in_use) &
write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc), & call print_message(loc, "Workspace is being used right now! " // &
"Workspace is being used right now! " // & "Corruptions can occur!")
"Corruptions can occur!"
end if
work_in_use = .true. work_in_use = .true.
@ -2121,10 +1959,6 @@ module user_problem
100 continue 100 continue
#ifdef PROFILE
call stop_timer(imt)
#endif /* PROFILE */
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end subroutine user_statistics end subroutine user_statistics
@ -2148,17 +1982,11 @@ module user_problem
! !
function log_cosh(x) result(y) function log_cosh(x) result(y)
! local variables are not implicit by default
!
implicit none implicit none
! function arguments
!
real(kind=8), intent(in) :: x real(kind=8), intent(in) :: x
real(kind=8) :: y real(kind=8) :: y
! local parameters
!
real(kind=8), parameter :: th = acosh(huge(x)), lh = log(0.5d+00) real(kind=8), parameter :: th = acosh(huge(x)), lh = log(0.5d+00)
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------