From 7ae74da6d6e4fef27f775d5338c289d253efea1c Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 24 Nov 2021 20:20:10 -0300 Subject: [PATCH] USER_PROBLEM: Rewrite slightly, remote useless comments. Signed-off-by: Grzegorz Kowal --- sources/user_problem.F90 | 704 +++++++++++++++------------------------ 1 file changed, 266 insertions(+), 438 deletions(-) diff --git a/sources/user_problem.F90 b/sources/user_problem.F90 index 900ffb6..b562cc0 100644 --- a/sources/user_problem.F90 +++ b/sources/user_problem.F90 @@ -31,14 +31,6 @@ module user_problem 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 :: zeta = 0.00d+00 real(kind=8), save :: eta = 0.00d+00 @@ -80,16 +72,8 @@ module user_problem logical, save :: absorption = .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 -! export subroutines -! private public :: initialize_user_problem, finalize_user_problem 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) -! include external procedures and variables -! #if NDIMS == 3 use constants , only : pi #endif /* NDIMS == 3 */ use constants , only : pi2 use coordinates, only : ng => nghosts, ady, xlen, zlen, ymax 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 random , only : randuni, randsym @@ -138,343 +120,289 @@ module user_problem logical , intent(in) :: verbose integer , intent(out) :: status -! local variables -! character(len=64) :: perturbation = "noise", append = "off", fname character(len=64) :: enable_absorption = "off" logical :: flag 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 real(kind=8) :: thv, tx, ty, tz, tt +#else /* NDIMS == 3 */ + real(kind=8) :: ka #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 -! this problem does not work with not magnetized set of equations -! if (.not. magnetized) then - if (verbose) then - write(*,*) - write(*,"(1x,a)") "ERROR!" - write(*,"(1x,a)") "The problem " // trim(problem) // & - " requires magnetized set of equations:" // & - " 'MHD', or 'SR-MHD'." - write(*,*) - end if + if (verbose) & + call print_message(loc, & + "The reconnection problem does not work without magnetic field.") status = 1 + return end if -! proceed if no errors -! - if (status == 0) then - ! get the reconnection equilibrium parameters ! - call get_parameter("beta", beta) - if (beta <= 0.0d+00) then - if (verbose) then - write(*,*) - write(*,"(1x,a)") "ERROR!" - write(*,"(1x,a)") "Parameter 'beta' must be positive (beta > 0.0)!" - write(*,*) - end if - status = 1 - end if - call get_parameter("zeta", zeta) - if (zeta < 0.0d+00 .or. zeta > 1.0d+00) then - if (verbose) then - write(*,*) - write(*,"(1x,a)") "ERROR!" - write(*,"(1x,a)") "Parameter 'zeta' must be between 0.0 and 1.0!" - write(*,*) - end if - status = 1 - end if - call get_parameter("resistivity", eta) - if (eta < 0.0d+00) then - if (verbose) then - write(*,*) - write(*,"(1x,a)") "ERROR!" - write(*,"(1x,a)") "Resistivity cannot be negative!" - write(*,*) - end if - status = 1 - else - resistive = .true. - end if - call get_parameter("dens", dens) - if (dens <= 0.0d+00) then - if (verbose) then - write(*,*) - write(*,"(1x,a)") "ERROR!" - write(*,"(1x,a)") "Parameter 'dens' must be positive (dens > 0.0)!" - write(*,*) - end if - status = 1 - end if - call get_parameter("bamp", bamp) - call get_parameter("bgui", bgui) + call get_parameter("beta", beta) + if (beta <= 0.0d+00) then + if (verbose) & + call print_message(loc, "Plasma-beta must be positive (beta > 0.0)!") + status = 1 + return + end if + call get_parameter("zeta", zeta) + if (zeta < 0.0d+00 .or. zeta > 1.0d+00) then + if (verbose) & + call print_message(loc, "Parameter 'zeta' must be between 0.0 and 1.0!") + status = 1 + return + end if + call get_parameter("resistivity", eta) + if (eta < 0.0d+00) then + if (verbose) & + call print_message(loc, "Resistivity cannot be negative!") + status = 1 + return + else + resistive = .true. + end if + call get_parameter("dens", dens) + if (dens <= 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("bamp", bamp) + 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-β ! parameters, and the sound speed in the case of isothermal equations of state ! - pmag = 0.5d+00 * (bamp**2 + bgui**2) - pres = beta * pmag - ptot = pres + pmag - csnd2 = pres / dens - csnd = sqrt(csnd2) - valf = sqrt(2.0d+00 * pmag / dens) - lund = valf / max(tiny(eta), eta) - 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)) - - end if ! status - -! proceed if no errors -! - if (status == 0) then + pmag = 0.5d+00 * (bamp**2 + bgui**2) + pres = beta * pmag + ptot = pres + pmag + csnd2 = pres / dens + csnd = sqrt(csnd2) + valf = sqrt(2.0d+00 * pmag / dens) + lund = valf / max(tiny(eta), eta) + dlta = lund**(- 1.0d+00 / 3.0d+00) + blim = max(blim, ng * ady(1)) ! get the perturbation parameters ! - call get_parameter("perturbation", perturbation) - call get_parameter("nper" , nper) - call get_parameter("bper" , bper) - call get_parameter("vper" , vper) - call get_parameter("kper" , kper) - call get_parameter("xcut" , xcut) - call get_parameter("ycut" , ycut) - call get_parameter("xdec" , xdec) - call get_parameter("ydec" , ydec) + call get_parameter("perturbation", perturbation) + call get_parameter("nper" , nper) + call get_parameter("bper" , bper) + call get_parameter("vper" , vper) + call get_parameter("kper" , kper) + call get_parameter("xcut" , xcut) + call get_parameter("ycut" , ycut) + call get_parameter("xdec" , xdec) + call get_parameter("ydec" , ydec) ! choose the perturbation type ! - select case(perturbation) - case('noise', 'random') - pert = 0 - case('mode', 'one mode', 'one-mode', 'one_mode') - pert = 1 - case('multi-wave', 'random waves', 'random-waves', 'random_waves') - pert = 2 - case default - perturbation = 'mode' - pert = 1 - end select + select case(perturbation) + case('noise', 'random') + pert = 0 + case('mode', 'one mode', 'one-mode', 'one_mode') + pert = 1 + case('multi-wave', 'random waves', 'random-waves', 'random_waves') + pert = 2 + case default + perturbation = 'mode' + pert = 1 + end select ! prepare the wave vector of the perturbation ! - kvec = pi2 * kper + kvec = pi2 * kper ! prepare the wave vectors for multi-wave perturbation ! - if (pert == 2) then + if (pert == 2) then ! allocate phase and wave vector components ! - allocate(kx(nper), ky(nper), kz(nper), ux(nper), uy(nper), uz(nper), & - ph(nper), stat = status) + allocate(kx(nper), ky(nper), kz(nper), ux(nper), uy(nper), uz(nper), & + ph(nper), stat = status) - if (status == 0) then + if (status == 0) then ! choose random wave vector directions ! - kd = int(xlen / zlen) - fc = 1.0d+00 / sqrt(1.0d+00 * nper) - do n = 1, nper - thh = pi2 * randuni() + kd = int(xlen / zlen) + fc = 1.0d+00 / sqrt(1.0d+00 * nper) + do n = 1, nper + thh = pi2 * randuni() #if NDIMS == 3 + thv = pi * randsym() + tx = cos(thh) * cos(thv) + ty = sin(thh) * cos(thv) + tz = sin(thv) + kx(n) = pi2 * nint(kper * tx) + ky(n) = pi2 * nint(kper * ty) + kz(n) = pi2 * nint(kper * tz / kd) * kd + tt = 0.0d+00 + do while(tt < 1.0d-08) + thh = pi2 * randuni() thv = pi * randsym() tx = cos(thh) * cos(thv) ty = sin(thh) * cos(thv) tz = sin(thv) - kx(n) = pi2 * nint(kper * tx) - ky(n) = pi2 * nint(kper * ty) - kz(n) = pi2 * nint(kper * tz / kd) * kd - tt = 0.0d+00 - do while(tt < 1.0d-08) - thh = pi2 * randuni() - thv = pi * randsym() - tx = cos(thh) * cos(thv) - ty = sin(thh) * cos(thv) - tz = sin(thv) - ux(n) = ty * kz(n) - tz * ky(n) - uy(n) = tz * kx(n) - tx * kz(n) - uz(n) = tx * ky(n) - ty * kx(n) - tt = sqrt(ux(n)**2 + uy(n)**2 + uz(n)**2) - end do - ux(n) = fc * ux(n) / tt - uy(n) = fc * uy(n) / tt - uz(n) = fc * uz(n) / tt -#else /* NDIMS == 3 */ - kx(n) = pi2 * nint(kper * cos(thh)) - ky(n) = pi2 * nint(kper * sin(thh)) - kz(n) = 0.0d+00 - ka = sqrt(kx(n)**2 + ky(n)**2) - ux(n) = fc * ky(n) / ka - uy(n) = - fc * kx(n) / ka - uz(n) = 0.0d+00 -#endif /* NDIMS == 3 */ - ph(n) = pi2 * randuni() + ux(n) = ty * kz(n) - tz * ky(n) + uy(n) = tz * kx(n) - tx * kz(n) + uz(n) = tx * ky(n) - ty * kx(n) + tt = sqrt(ux(n)**2 + uy(n)**2 + uz(n)**2) end do + ux(n) = fc * ux(n) / tt + uy(n) = fc * uy(n) / tt + uz(n) = fc * uz(n) / tt +#else /* NDIMS == 3 */ + kx(n) = pi2 * nint(kper * cos(thh)) + ky(n) = pi2 * nint(kper * sin(thh)) + kz(n) = 0.0d+00 + ka = sqrt(kx(n)**2 + ky(n)**2) + ux(n) = fc * ky(n) / ka + uy(n) = - fc * kx(n) / ka + uz(n) = 0.0d+00 +#endif /* NDIMS == 3 */ + ph(n) = pi2 * randuni() + end do - end if ! status + + else + if (verbose) & + call print_message(loc, & + "Could not allocate space for perturbation vectors!") + return end if - end if ! status - -! prepare file to store reconnection rate terms -! - if (status == 0) then + end if ! pert = 2 ! determine if to append or create another file ! - call get_parameter("reconnection_append", append) - select case(trim(append)) - case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") - write(fname, "('reconnection.dat')") - inquire(file = fname, exist = flag) - case default - write(fname, "('reconnection_',i2.2,'.dat')") rcount - flag = .false. - end select + call get_parameter("reconnection_append", append) + select case(trim(append)) + case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") + write(fname, "('reconnection.dat')") + inquire(file = fname, exist = flag) + case default + write(fname, "('reconnection_',i2.2,'.dat')") rcount + flag = .false. + end select ! check if the file exists; if not, create a new one, otherwise move to the end ! - if (flag .and. rcount > 1) then + if (flag .and. rcount > 1) then #ifdef __INTEL_COMPILER - open(newunit = runit, file = fname, form = 'formatted', status = 'old' & + open(newunit = runit, file = fname, form = 'formatted', status = 'old' & , position = 'append', buffered = 'yes') #else /* __INTEL_COMPILER */ - open(newunit = runit, file = fname, form = 'formatted', status = 'old' & + open(newunit = runit, file = fname, form = 'formatted', status = 'old' & , position = 'append') #endif /* __INTEL_COMPILER */ - write(runit,"('#')") - else + write(runit,"('#')") + else #ifdef __INTEL_COMPILER - open(newunit = runit, file = fname, form = 'formatted' & + open(newunit = runit, file = fname, form = 'formatted' & , status = 'replace', buffered = 'yes') #else /* __INTEL_COMPILER */ - open(newunit = runit, file = fname, form = 'formatted' & + open(newunit = runit, file = fname, form = 'formatted' & , status = 'replace') #endif /* __INTEL_COMPILER */ - end if + end if ! write the integral file header ! - write(runit,'("#",a24,26a25)') & - 'time', '|Bx| int', '|Bx| inf' & - , '|Bx| y-adv', '|Bx| z-adv', '|Bx| y-shr', '|Bx| z-shr' & - , '|Bx| y-dif', '|Bx| z-dif', '|Bx| Psi' & - , 'Vin lower' , 'Vin upper' & - , 'Emag', 'Emag inf' & - , 'Emag x-adv', 'Emag y-adv', 'Emag z-adv' & - , 'Emag x-v.b', 'Emag y-v.b', 'Emag z-v.b' & - , 'Emag x-dif', 'Emag y-dif', 'Emag z-dif' & - , 'Emag-Ekin', 'Emag-Eint', 'Emag-Psi' - write(runit,"('#')") + write(runit,'("#",a24,26a25)') & + 'time', '|Bx| int', '|Bx| inf' & + , '|Bx| y-adv', '|Bx| z-adv', '|Bx| y-shr', '|Bx| z-shr' & + , '|Bx| y-dif', '|Bx| z-dif', '|Bx| Psi' & + , 'Vin lower' , 'Vin upper' & + , 'Emag', 'Emag inf' & + , 'Emag x-adv', 'Emag y-adv', 'Emag z-adv' & + , 'Emag x-v.b', 'Emag y-v.b', 'Emag z-v.b' & + , 'Emag x-dif', 'Emag y-dif', 'Emag z-dif' & + , 'Emag-Ekin', 'Emag-Eint', 'Emag-Psi' + write(runit,"('#')") - call get_parameter("ydistance", ydis) - ydis = min(ydis, ymax) - ylo = - ydis - yup = ydis + call get_parameter("ydistance", ydis) + ydis = min(ydis, ymax) + ylo = - ydis + yup = ydis ! absorption parameters ! - call get_parameter("enable_shapes", enable_absorption) - select case(trim(enable_absorption)) - case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") - absorption = .true. - case default - absorption = .false. - end select - call get_parameter("tabs", tabs) - call get_parameter("adif", adif) - call get_parameter("acut", acut) - call get_parameter("adec", adec) - yabs = ymax - abs(acut) - - end if ! status - -! proceed if no errors -! - if (status == 0) then + call get_parameter("enable_shapes", enable_absorption) + select case(trim(enable_absorption)) + case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") + absorption = .true. + case default + absorption = .false. + end select + call get_parameter("tabs", tabs) + call get_parameter("adif", adif) + call get_parameter("acut", acut) + call get_parameter("adec", adec) + yabs = ymax - abs(acut) ! print information about the user problem setup ! - call print_section(verbose, "Parameters (* - set, otherwise calculated)") - call print_parameter(verbose, '(*) beta (plasma-beta)' , beta) - call print_parameter(verbose, '(*) zeta' , zeta) - call print_parameter(verbose, '(*) dens' , dens) - call print_parameter(verbose, '(*) bamp' , bamp) - call print_parameter(verbose, '(*) bgui (guide field)' , bgui) - call print_parameter(verbose, '( ) pres (thermal pres.)' , pres) - call print_parameter(verbose, '( ) pmag (magnetic pres.)', pmag) - call print_parameter(verbose, '( ) ptot (total pressure)', ptot) - call print_parameter(verbose, '( ) csnd (sound speed)' , csnd) - call print_parameter(verbose, '( ) Valf (Alfven speed)' , valf) - if (resistive) then - call print_parameter(verbose, '( ) S (Lundquist number)', lund) - end if - call print_parameter(verbose, '(*) delta (thickness)', dlta) - call print_parameter(verbose, '(*) blim' , blim) - call print_parameter(verbose, '(*) perturbation', perturbation) - call print_parameter(verbose, '(*) bper' , bper) - call print_parameter(verbose, '(*) vper' , vper) - if (pert >= 1) then - call print_parameter(verbose, '(*) kper' , kper) - end if - if (pert == 2) then - call print_parameter(verbose, '(*) nper' , nper) - end if - call print_parameter(verbose, '(*) xcut' , xcut) - call print_parameter(verbose, '(*) ycut' , ycut) - call print_parameter(verbose, '(*) xdec' , xdec) - call print_parameter(verbose, '(*) ydec' , ydec) - call print_parameter(verbose, '(*) ydistance (Vrec)', ydis) - if (absorption) then - call print_parameter(verbose, '(*) tabs (absorption)', tabs) - call print_parameter(verbose, '(*) adif (diffusion)' , adif) - call print_parameter(verbose, '(*) acut' , acut) - call print_parameter(verbose, '(*) adec' , adec) - call print_parameter(verbose, '( ) yabs' , yabs) - end if - end if ! status + call print_section(verbose, "Parameters (* - set, otherwise calculated)") + call print_parameter(verbose, '(*) beta (plasma-beta)' , beta) + call print_parameter(verbose, '(*) zeta' , zeta) + call print_parameter(verbose, '(*) dens' , dens) + call print_parameter(verbose, '(*) bamp' , bamp) + call print_parameter(verbose, '(*) bgui (guide field)' , bgui) + call print_parameter(verbose, '( ) pres (thermal pres.)' , pres) + call print_parameter(verbose, '( ) pmag (magnetic pres.)', pmag) + call print_parameter(verbose, '( ) ptot (total pressure)', ptot) + call print_parameter(verbose, '( ) csnd (sound speed)' , csnd) + call print_parameter(verbose, '( ) Valf (Alfven speed)' , valf) + if (resistive) then + call print_parameter(verbose, '( ) S (Lundquist number)', lund) + end if + call print_parameter(verbose, '(*) delta (thickness)', dlta) + call print_parameter(verbose, '(*) blim' , blim) + call print_parameter(verbose, '(*) perturbation', perturbation) + call print_parameter(verbose, '(*) bper' , bper) + call print_parameter(verbose, '(*) vper' , vper) + if (pert >= 1) then + call print_parameter(verbose, '(*) kper' , kper) + end if + if (pert == 2) then + call print_parameter(verbose, '(*) nper' , nper) + end if + call print_parameter(verbose, '(*) xcut' , xcut) + call print_parameter(verbose, '(*) ycut' , ycut) + call print_parameter(verbose, '(*) xdec' , xdec) + call print_parameter(verbose, '(*) ydec' , ydec) + call print_parameter(verbose, '(*) ydistance (Vrec)', ydis) + if (absorption) then + call print_parameter(verbose, '(*) tabs (absorption)', tabs) + call print_parameter(verbose, '(*) adif (diffusion)' , adif) + call print_parameter(verbose, '(*) acut' , acut) + call print_parameter(verbose, '(*) adec' , adec) + call print_parameter(verbose, '( ) yabs' , yabs) + end if !------------------------------------------------------------------------------- ! @@ -495,29 +423,23 @@ module user_problem ! subroutine finalize_user_problem(status) + use helpers, only : print_message + implicit none integer, intent(out) :: status + character(len=*), parameter :: loc = 'USER_PROBLEM::finalize_user_problem()' + !------------------------------------------------------------------------------- ! status = 0 -! close the reconnection file -! close(runit) -! deallocate wave vector components, random directions, and random phase -! - if (allocated(kx)) deallocate(kx, ky, kz, stat = status) - 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 */ + deallocate(kx, ky, kz, ux, uy, uz, ph, stat = status) + call print_message(loc, & + "Could not deallocate space for perturbation vectors!") !------------------------------------------------------------------------------- ! @@ -539,21 +461,21 @@ module user_problem ! subroutine setup_user_problem(pdata) - use blocks , only : block_data - use constants , only : pi, pi2 - use coordinates , only : nn => bcells - use coordinates , only : ax, ay, adx, ady, ylen + use blocks , only : block_data + use constants , only : pi, pi2 + use coordinates, only : nn => bcells + use coordinates, only : ax, ay, adx, ady, ylen #if NDIMS == 3 - use coordinates , only : az, adz + use coordinates, only : az, adz #endif /* NDIMS == 3 */ - use equations , only : prim2cons - use equations , only : nv, ns - use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl - use equations , only : csnd2 - use iso_fortran_env, only : error_unit - use operators , only : curl - use random , only : randsym - use workspace , only : resize_workspace, work, work_in_use + use equations , only : prim2cons + use equations , only : nv, ns + use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl + use equations , only : csnd2 + use helpers , only : print_message + use operators , only : curl + use random , only : randsym + use workspace , only : resize_workspace, work, work_in_use implicit none @@ -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 n = (nv + 3) * nn**NDIMS call resize_workspace(n, status) if (status /= 0) then - write(error_unit,"('[',a,']: ',a)") trim(loc), & - "Could not resize the workspace!" + call print_message(loc, "Could not resize the workspace!") go to 100 end if @@ -617,8 +532,6 @@ module user_problem first = .false. end if -! prepare cell sizes and block coordinates -! dh(1) = adx(pdata%meta%level) dh(2) = ady(pdata%meta%level) #if NDIMS == 3 @@ -675,20 +588,14 @@ module user_problem fy(j) = cos(yp)**2 end do ! i = 1, nn - if (work_in_use) then - write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc), & - "Workspace is being used right now! " // & - "Corruptions can occur!" - end if + if (work_in_use) & + call print_message(loc, "Workspace is being used right now! " // & + "Corruptions can occur!") work_in_use = .true. -! reset the perturbation matrix -! qpert(:,:,:,:) = 0.0d+00 -! the random perturbation -! if (pert == 0) then ! of velocity @@ -976,10 +883,6 @@ module user_problem 100 continue -#ifdef PROFILE - call stop_timer(imp) -#endif /* PROFILE */ - !------------------------------------------------------------------------------- ! end subroutine setup_user_problem @@ -1034,16 +937,16 @@ module user_problem ! subroutine update_user_shapes(pdata, time, dt) - use blocks , only : block_data - use constants , only : pi - use coordinates , only : nn => bcells - use coordinates , only : ymax, ay, adx, ady, adz - use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp - use equations , only : prim2cons - use iso_fortran_env, only : error_unit - use operators , only : laplace - use workspace , only : resize_workspace, work, work_in_use + use blocks , only : block_data + use constants , only : pi + use coordinates, only : nn => bcells + use coordinates, only : ay, adx, ady, adz + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp + use equations , only : prim2cons + use helpers , only : print_message + use operators , only : laplace + use workspace , only : resize_workspace, work, work_in_use implicit none @@ -1065,17 +968,12 @@ module user_problem !------------------------------------------------------------------------------- ! -#ifdef PROFILE - call start_timer(ims) -#endif /* PROFILE */ - if (first) then j = nn**NDIMS call resize_workspace(j, status) if (status /= 0) then - write(error_unit,"('[',a,']: ',a)") trim(loc), & - "Could not resize the workspace!" + call print_message(loc, "Could not resize the workspace!") go to 100 end if @@ -1088,11 +986,9 @@ module user_problem first = .false. end if - if (work_in_use) then - write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc), & - "Workspace is being used right now! " // & - "Corruptions can occur!" - end if + if (work_in_use) & + call print_message(loc, "Workspace is being used right now! " // & + "Corruptions can occur!") work_in_use = .true. @@ -1151,10 +1047,6 @@ module user_problem 100 continue -#ifdef PROFILE - call stop_timer(ims) -#endif /* PROFILE */ - !------------------------------------------------------------------------------- ! 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) -! import external procedures and variables -! 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 use equations , only : ibz #endif /* NDIMS == 3 */ -! local variables are not implicit by default -! implicit none integer , intent(in) :: ic @@ -1234,8 +1122,6 @@ module user_problem real(kind=8), dimension(:) , intent(in) :: z real(kind=8), dimension(:,:,:,:), intent(inout) :: qn -! local variables -! integer :: im2, im1, i , ip1, ip2 integer :: jm2, jm1, j , jp1, jp2 integer :: k = 1 @@ -1249,18 +1135,8 @@ module user_problem !------------------------------------------------------------------------------- ! -#ifdef PROFILE -! start accounting time for the boundary update -! - call start_timer(imb) -#endif /* PROFILE */ + if (magnetized) then -! 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) dy = y(2) - y(1) #if NDIMS == 3 @@ -1271,8 +1147,6 @@ module user_problem dxz = dx / dz #endif /* NDIMS == 3 */ -! process left and right side boundary separatelly -! if (ic == 1) then ! iterate over left-side ghost layers @@ -1393,13 +1267,7 @@ module user_problem #endif /* NDIMS == 3 */ end do ! i = neu, nn end if - end if ! ibx > 0 - -#ifdef PROFILE -! stop accounting time for the boundary update -! - call stop_timer(imb) -#endif /* PROFILE */ + end if ! magnetized !------------------------------------------------------------------------------- ! @@ -1426,14 +1294,10 @@ module user_problem ! 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 equations , only : nv - use equations , only : idn, ivy, ipr, ibx, iby, ibz, ibp + use coordinates, only : nn => bcells, nb, ne, nbl, neu + use equations , only : magnetized, nv + use equations , only : idn, ivy, ipr, ibx, iby, ibz, ibp -! local variables are not implicit by default -! implicit none integer , intent(in) :: jc @@ -1445,8 +1309,6 @@ module user_problem real(kind=8), dimension(:) , intent(in) :: z real(kind=8), dimension(:,:,:,:), intent(inout) :: qn -! local variables -! integer :: i, im1, ip1 integer :: j, jm1, jp1, jm2, jp2 integer :: k = 1 @@ -1461,18 +1323,8 @@ module user_problem !------------------------------------------------------------------------------- ! -#ifdef PROFILE -! start accounting time for the boundary update -! - call start_timer(imb) -#endif /* PROFILE */ + if (magnetized) then -! 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) dy = y(2) - y(1) #if NDIMS == 3 @@ -1620,12 +1472,6 @@ module user_problem end if end if ! ibx > 0 -#ifdef PROFILE -! stop accounting time for the boundary update -! - call stop_timer(imb) -#endif /* PROFILE */ - !------------------------------------------------------------------------------- ! end subroutine user_boundary_y @@ -1686,17 +1532,16 @@ module user_problem ! subroutine user_statistics(time) - use blocks , only : block_meta, block_data, list_data - use coordinates , only : nn => bcells, nb, ne - use coordinates , only : adx, ady, adz, advol, ay, yarea - use equations , only : ivx, ivy, ivz, ibx, iby, ibz, ibp - use helpers , only : flush_and_sync - use iso_fortran_env, only : error_unit + use blocks , only : block_meta, block_data, list_data + use coordinates, only : nn => bcells, nb, ne + use coordinates, only : adx, ady, adz, advol, ay, yarea + use equations , only : ivx, ivy, ivz, ibx, iby, ibz, ibp + use helpers , only : print_message, flush_and_sync #ifdef MPI - use mpitools , only : reduce_sum + use mpitools , only : reduce_sum #endif /* MPI */ - use operators , only : curl, gradient - use workspace , only : resize_workspace, work, work_in_use + use operators , only : curl, gradient + use workspace , only : resize_workspace, work, work_in_use implicit none @@ -1719,17 +1564,12 @@ module user_problem !------------------------------------------------------------------------------- ! -#ifdef PROFILE - call start_timer(imt) -#endif /* PROFILE */ - if (first) then ni = 3 * nn**NDIMS call resize_workspace(ni, status) if (status /= 0) then - write(error_unit,"('[',a,']: ',a)") trim(loc), & - "Could not resize the workspace!" + call print_message(loc, "Could not resize the workspace!") go to 100 end if @@ -1744,11 +1584,9 @@ module user_problem rterms(:) = 0.0d+00 - if (work_in_use) then - write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc), & - "Workspace is being used right now! " // & - "Corruptions can occur!" - end if + if (work_in_use) & + call print_message(loc, "Workspace is being used right now! " // & + "Corruptions can occur!") work_in_use = .true. @@ -2121,10 +1959,6 @@ module user_problem 100 continue -#ifdef PROFILE - call stop_timer(imt) -#endif /* PROFILE */ - !------------------------------------------------------------------------------- ! end subroutine user_statistics @@ -2148,17 +1982,11 @@ module user_problem ! function log_cosh(x) result(y) -! local variables are not implicit by default -! implicit none -! function arguments -! real(kind=8), intent(in) :: x real(kind=8) :: y -! local parameters -! real(kind=8), parameter :: th = acosh(huge(x)), lh = log(0.5d+00) !-------------------------------------------------------------------------------