!!****************************************************************************** !! !! 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) 2017-2023 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: USER_PROBLEM !! !! This module provides subroutines to setup custom problem. !! !!***************************************************************************** ! module user_problem implicit none real(kind=8), save :: alpha = 5.00d-01 real(kind=8), save :: beta = 1.00d+00 real(kind=8), save :: zeta = 0.00d+00 real(kind=8), save :: eta = 0.00d+00 real(kind=8), save :: nu = 0.00d+00 real(kind=8), save :: dens = 1.00d+00 real(kind=8), save :: bamp = 1.00d+00 real(kind=8), save :: btra = 0.00d+00 real(kind=8), save :: bgui = 0.00d+00 real(kind=8), save :: pres = 5.00d-01 real(kind=8), save :: pmag = 5.00d-01 real(kind=8), save :: ptot = 1.00d+00 real(kind=8), save :: valf = 1.00d+00 real(kind=8), save :: lund = 1.00d+00 real(kind=8), save :: prdl = 0.00d+00 real(kind=8), save :: dlta = 1.00d-16 real(kind=8), save :: tdec = 1.00d+00 integer , save :: pert = 0 integer , save :: nper = 16 real(kind=8), save :: bper = 0.00d+00 real(kind=8), save :: vper = 0.00d+00 real(kind=8), save :: kper = 1.00d+00 real(kind=8), save :: kvec = 1.00d+00 real(kind=8), save :: xcut = 1.00d+99 real(kind=8), save :: ycut = 1.00d+99 real(kind=8), save :: xdec = 1.00d-01 real(kind=8), save :: ydec = 1.00d-01 real(kind=8), save :: ylo = -9.00d+99 real(kind=8), save :: yup = 9.00d+99 real(kind=8), save :: tabs = 1.00d+00 real(kind=8), save :: adif = 5.00d-01 real(kind=8), save :: acut = 1.00d+00 real(kind=8), save :: adec = 1.00d+00 real(kind=8), save :: yabs = 9.00d+99 integer(kind=4), save :: runit = 33 logical, save :: absorption = .false. logical, save :: resistive = .false. real(kind=8), dimension(:), allocatable :: kx, ky, kz, ux, uy, uz, ph private public :: initialize_user_problem, finalize_user_problem public :: user_statistics !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== ! ! subroutine INITIALIZE_USER_PROBLEM: ! ---------------------------------- ! ! Subroutine initializes user problem. It could read problem parameters which ! are used in all subroutines defining this specific problem. ! ! Arguments: ! ! problem - the problem name ! rcount - the run count for restarted jobs ! verbose - a logical flag turning the information printing; ! status - an integer flag for error return value; ! !=============================================================================== ! subroutine initialize_user_problem(problem, rcount, verbose, status) use boundaries , only : custom_boundary_x, custom_boundary_y #if NDIMS == 3 use constants , only : pi #endif /* NDIMS == 3 */ use constants , only : pi2 use coordinates, only : xlen, zlen, ymax use equations , only : magnetized, csnd, csnd2, qpbnd use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp, ipr use helpers , only : print_section, print_parameter, print_message use mesh , only : setup_problem use parameters , only : get_parameter use random , only : randuni, randsym use shapes , only : update_shapes implicit none character(len=64), intent(in) :: problem integer , intent(in) :: rcount logical , intent(in) :: verbose integer , intent(out) :: status character(len=64) :: perturbation = "noise", append = "off", fname character(len=64) :: enable_absorption = "off" logical :: flag integer :: n, kd 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()' !------------------------------------------------------------------------------- ! status = 0 if (.not. magnetized) then if (verbose) & call print_message(loc, & "The reconnection problem does not work without magnetic field.") status = 1 return end if ! get the reconnection equilibrium parameters ! 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("viscosity" , nu ) 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("btra", btra) call get_parameter("bgui", bgui) call get_parameter("boundary_decay_time", tdec) if (tdec <= 0.0d+00) then if (verbose) & call print_message(loc, "Boundary decay time must be larger than zero!") status = 1 return end if ! get the index for the current sheet thickness scaling ! call get_parameter("alpha", alpha) if (alpha <= 0.0d+00) then if (verbose) & call print_message(loc, "Current sheet thickness scaling index " // & "(alpha) must be larger than zero! " // & "Resetting it to the default value of 1/2.") alpha = 5.0d-01 end if ! 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(bamp**2 / dens) lund = valf / max(tiny(eta), eta) prdl = nu / max(tiny(eta), eta) dlta = lund**(-alpha) ! fill up the boundary values ! qpbnd(idn,2,:) = dens qpbnd(ivx,2,:) = 0.0d+00 qpbnd(ivy,2,:) = 0.0d+00 qpbnd(ivz,2,:) = 0.0d+00 qpbnd(ibx,2,:) = [ - bamp, bamp ] qpbnd(iby,2,:) = btra qpbnd(ibz,2,:) = bgui qpbnd(ibp,2,:) = 0.0d+00 if (ipr > 0) qpbnd(ipr,2,:) = pres ! 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 ! 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) ! 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('wave-spectrum', 'wave spectrum', 'wave_spectrum') pert = 3 case('external') pert = 4 case default perturbation = 'mode' pert = 1 end select ! prepare the wave vector of the perturbation ! kvec = pi2 * kper ! prepare the wave vectors for multi-wave perturbation ! 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) 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() #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) 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 else if (verbose) & call print_message(loc, & "Could not allocate space for perturbation vectors!") return end if end if ! pert = 2 ! prepare the wave vectors for wave-spectrum perturbation ! if (pert == 3) then ! allocate phase and wave vector components ! allocate(kx(nper), ph(nper), stat = status) if (status == 0) then do n = 1, nper kx(n) = pi2 * n / xlen ph(n) = pi2 * randuni() end do else if (verbose) & call print_message(loc, & "Could not allocate space for perturbation vectors!") return end if end if ! pert = 3 ! 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 ! check if the file exists; if not, create a new one, otherwise move to the end ! if (flag .and. rcount > 1) then #ifdef __INTEL_COMPILER 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' & , position = 'append') #endif /* __INTEL_COMPILER */ write(runit,"('#')") else #ifdef __INTEL_COMPILER open(newunit = runit, file = fname, form = 'formatted' & , status = 'replace', buffered = 'yes') #else /* __INTEL_COMPILER */ open(newunit = runit, file = fname, form = 'formatted' & , status = 'replace') #endif /* __INTEL_COMPILER */ 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,"('#')") 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) ! 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, '(*) btra (transverse)' , btra) 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, '( ) P (Prandtl number)', prdl) call print_parameter(verbose, '(*) delta (thickness)', dlta) 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, '(*) boundary decay time', tdec) 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 ! set procedure pointers for the problem setup subroutine, the shapes update ! for the absorption layer, and boundary conditions; ! setup_problem => setup_user_problem custom_boundary_x => user_boundary_x custom_boundary_y => user_boundary_y if (absorption) & update_shapes => update_user_shapes !------------------------------------------------------------------------------- ! end subroutine initialize_user_problem ! !=============================================================================== ! ! subroutine FINALIZE_USER_PROBLEM: ! -------------------------------- ! ! Subroutine releases memory used by the module. ! ! Arguments: ! ! status - an integer flag for error return value; ! !=============================================================================== ! 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(runit) if (pert == 2) then deallocate(kx, ky, kz, ux, uy, uz, ph, stat = status) if (status /= 0) & call print_message(loc, & "Could not deallocate space for perturbation vectors!") end if if (pert == 3) then deallocate(kx, ph, stat = status) if (status /= 0) & call print_message(loc, & "Could not deallocate space for perturbation vectors!") end if !------------------------------------------------------------------------------- ! end subroutine finalize_user_problem ! !=============================================================================== ! ! subroutine SETUP_USER_PROBLEM: ! ----------------------------- ! ! Subroutine sets the initial conditions for the user specific problem. ! ! Arguments: ! ! pdata - pointer to the datablock structure of the currently initialized ! block; ! !=============================================================================== ! subroutine setup_user_problem(pdata) use blocks , only : block_data use constants , only : pi, pi2 use coordinates, only : nn => bcells, nc => ncells, ng => nghosts use coordinates, only : ax, ay, adx, ady, ymin, ylen #if NDIMS == 3 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 helpers , only : print_message use operators , only : curl use random , only : randsym use workspace , only : resize_workspace, work, work_in_use implicit none type(block_data), pointer, intent(inout) :: pdata logical, save :: first = .true. character(len=64) :: filename integer :: i, j, k, n, io, status, jl, ju, p real(kind=8) :: xpl, xpu, ypl, ypu real(kind=8) :: xp, yp real(kind=8) :: yrat real(kind=8) :: kv, fv real(kind=8) :: vx = 0.0d+00, vy = 0.0d+00, vv, va #if NDIMS == 3 real(kind=8) :: vz = 0.0d+00 #endif /* NDIMS == 3 */ real(kind=8) :: pz = 0.0d+00, pp, ba #if NDIMS == 3 real(kind=8) :: px = 0.0d+00, py = 0.0d+00 #endif /* NDIMS == 3 */ integer, save :: nt !$ integer :: omp_get_thread_num real(kind=8), dimension(:,:,:,:), pointer, save :: qpert, pot !$omp threadprivate(first, nt, qpert, pot) real(kind=8), dimension(nv,nn) :: q, u, qprof real(kind=8), dimension(nn) :: xl, xu, xc, fx real(kind=8), dimension(nn) :: yl, yu, yc, fy #if NDIMS == 3 real(kind=8), dimension(nn) :: zc #endif /* NDIMS == 3 */ real(kind=8), dimension(3) :: dh real(kind=8), dimension(:) , allocatable :: yt real(kind=8), dimension(:,:), allocatable :: qr, qi character(len=*), parameter :: loc = 'USER_PROBLEM:setup_user_problem()' !------------------------------------------------------------------------------- ! #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ nt = 0 !$ nt = omp_get_thread_num() if (first) then n = (nv + 3) * nn**NDIMS call resize_workspace(n, status) if (status /= 0) then call print_message(loc, "Could not resize the workspace!") go to 100 end if i = nv * nn**NDIMS #if NDIMS == 3 qpert(1:nv,1:nn,1:nn,1:nn) => work( 1:i,nt) pot (1:3,1:nn,1:nn,1:nn) => work(i+1:n,nt) #else /* NDIMS == 3 */ qpert(1:nv,1:nn,1:nn,1: 1) => work( 1:i,nt) pot (1:3,1:nn,1:nn,1: 1) => work(i+1:n,nt) #endif /* NDIMS == 3 */ first = .false. end if dh(1) = adx(pdata%meta%level) dh(2) = ady(pdata%meta%level) #if NDIMS == 3 dh(3) = adz(pdata%meta%level) #endif /* NDIMS == 3 */ xc(:) = pdata%meta%xmin + ax(pdata%meta%level,:) yc(:) = pdata%meta%ymin + ay(pdata%meta%level,:) #if NDIMS == 3 zc(:) = pdata%meta%zmin + az(pdata%meta%level,:) #endif /* NDIMS == 3 */ xl(:) = pdata%meta%xmin + ax(pdata%meta%level,:) - 5.0d-01 * dh(1) xu(:) = xl(:) + dh(1) yl(:) = pdata%meta%ymin + ay(pdata%meta%level,:) - 5.0d-01 * dh(2) yu(:) = yl(:) + dh(2) ! set the equilibrium profile ! do j = 1, nn qprof(ibx,j) = bamp * max(-1.0d+00, min(1.0d+00, & dlta * (log_cosh(yu(j) / dlta) & - log_cosh(yl(j) / dlta)) / dh(2))) end do qprof(iby,:) = btra qprof(ibz,:) = zeta * sqrt(bamp**2 - qprof(ibx,:)**2) + bgui if (ipr > 0) then qprof(idn,:) = dens qprof(ipr,:) = ptot - 0.5d+00 * sum(qprof(ibx:ibz,:)**2, 1) else qprof(idn,:) = (ptot - 0.5d+00 * sum(qprof(ibx:ibz,:)**2, 1)) / csnd2 end if qprof(ivx,:) = 0.0d+00 qprof(ivy,:) = 0.0d+00 qprof(ivz,:) = 0.0d+00 qprof(ibp,:) = 0.0d+00 if (ns > 0) then qprof(isl,:) = sign(1.0d+00, yc(:)) end if ! ratio of the current sheet thickness to the cell size ! yrat = dlta / dh(2) ! prepare decaying factors ! fv = 0.5d+00 * pi do i = 1, nn xp = fv * min(1.0d+00, max(0.0d+00, abs(xc(i)) - xcut) / xdec) fx(i) = cos(xp)**2 end do ! i = 1, nn do j = 1, nn yp = fv * min(1.0d+00, max(0.0d+00, abs(yc(j)) - ycut) / ydec) fy(j) = cos(yp)**2 end do ! i = 1, nn if (work_in_use(nt)) & call print_message(loc, "Workspace is being used right now! " // & "Corruptions can occur!") work_in_use(nt) = .true. qpert(:,:,:,:) = 0.0d+00 if (pert == 0) then ! of velocity ! if (abs(vper) > 0.0d+00) then ! initiate the random velocity components ! #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn if (fy(j) > 0.0d+00) then do i = 1, nn if (fx(i) > 0.0d+00) then ! calculate the profile of perturbation amplitude ! va = vper * fx(i) * fy(j) ! get the random direction ! vv = 0.0d+00 do while(vv < 1.0d-08) vx = randsym() vy = randsym() #if NDIMS == 3 vz = randsym() #endif /* NDIMS == 3 */ #if NDIMS == 3 vv = sqrt(vx * vx + vy * vy + vz * vz) #else /* NDIMS == 3 */ vv = sqrt(vx * vx + vy * vy) #endif /* NDIMS == 3 */ end do qpert(ivx,i,j,k) = va * (vx / vv) qpert(ivy,i,j,k) = va * (vy / vv) #if NDIMS == 3 qpert(ivz,i,j,k) = va * (vz / vv) #endif /* NDIMS == 3 */ end if ! |x| < xcut end do ! i = 1, nn end if ! |y| < ycut end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ end if ! vper /= 0.0 ! of magnetic field ! if (abs(bper) > 0.0d+00) then ! reset the potential ! pot(:,:,:,:) = 0.0d+00 ! initiate the random magnetic field components ! #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn if (fy(j) > 0.0d+00) then do i = 1, nn if (fx(i) > 0.0d+00) then ! calculate the profile of perturbation amplitude ! ba = dh(1) * bper * fx(i) * fy(j) ! get the random direction ! pp = 0.0d+00 do while(pp < 1.0d-08) #if NDIMS == 3 px = randsym() py = randsym() #endif /* NDIMS == 3 */ pz = randsym() #if NDIMS == 3 pp = sqrt(px * px + py * py + pz * pz) #else /* NDIMS == 3 */ pp = abs(pz) #endif /* NDIMS == 3 */ end do #if NDIMS == 3 pot(1,i,j,k) = ba * (px / pp) pot(2,i,j,k) = ba * (py / pp) #endif /* NDIMS == 3 */ pot(3,i,j,k) = ba * (pz / pp) end if ! |x| < xcut end do ! i = 1, nn end if ! |y| < ycut end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ ! calculate magnetic field perturbation components from vector potential ! call curl(dh(1:3), pot(:,:,:,:), qpert(ibx:ibz,:,:,:)) end if ! bper /= 0.0 end if ! pert == 0 ! one-mode perturbation ! if (pert == 1) then ! of velocity ! if (abs(vper) > 0.0d+00) then ! reset the potential ! pot(:,:,:,:) = 0.0d+00 ! calculate the perturbation factor ! fv = vper * ylen / (pi2 * pi2 * pi * dh(1) * dh(2) * kper) ! calculate the perturbation profile ! #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn if (fy(j) > 0.0d+00) then ypl = pi * yl(j) / ylen ypu = pi * yu(j) / ylen do i = 1, nn if (fx(i) > 0.0d+00) then xpl = pi2 * kper * xl(i) xpu = pi2 * kper * xu(i) pot(3,i,j,k) = fv * fx(i) * fy(j) * (cos(xpl) - cos(xpu)) & * (cos(ypl) - cos(ypu)) end if ! |x| < xcut end do ! i = 1, nn end if ! |y| < ycut end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ ! calculate magnetic field components from vector potential ! call curl(dh(1:3), pot(:,:,:,:), qpert(ivx:ivz,:,:,:)) end if ! vper /= 0.0 ! of magnetic field ! if (abs(bper) > 0.0d+00) then ! reset the potential ! pot(:,:,:,:) = 0.0d+00 ! calculate the perturbation factor ! fv = bper * ylen / (pi2 * pi2 * pi * dh(1) * dh(2) * kper) ! calculate the perturbation profile ! #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn if (fy(j) > 0.0d+00) then ypl = pi * yl(j) / ylen ypu = pi * yu(j) / ylen do i = 1, nn if (fx(i) > 0.0d+00) then xpl = pi2 * kper * xl(i) xpu = pi2 * kper * xu(i) pot(3,i,j,k) = fv * fx(i) * fy(j) * (sin(xpl) - sin(xpu)) & * (sin(ypl) - sin(ypu)) end if ! |x| < xcut end do ! i = 1, nn end if ! |y| < ycut end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ ! calculate magnetic field components from vector potential ! call curl(dh(1:3), pot(:,:,:,:), qpert(ibx:ibz,:,:,:)) end if ! bper /= 0.0 end if ! pert == 1 ! prepare the random perturbation of velocity ! if (pert == 2) then if (abs(vper) > 0.0d+00) then ! iterate over the block position and initiate the velocity perturbation ! #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn if (fy(j) > 0.0d+00) then do i = 1, nn if (fx(i) > 0.0d+00) then ! calculate the velocity amplitude profile ! fv = vper * fx(i) * fy(j) ! add perturbation components ! do n = 1, nper #if NDIMS == 3 kv = kx(n) * xc(i) + ky(n) * yc(j) + kz(n) * zc(k) + ph(n) #else /* NDIMS == 3 */ kv = kx(n) * xc(i) + ky(n) * yc(j) + ph(n) #endif /* NDIMS == 3 */ va = fv * sin(kv) qpert(ivx,i,j,k) = qpert(ivx,i,j,k) + va * ux(n) qpert(ivy,i,j,k) = qpert(ivy,i,j,k) + va * uy(n) #if NDIMS == 3 qpert(ivz,i,j,k) = qpert(ivz,i,j,k) + va * uz(n) #endif /* NDIMS == 3 */ end do end if ! fx > 0.0 end do ! i = 1, nn end if ! fy > 0.0 end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ end if ! vper /= 0.0 end if ! pert == 2 ! prepare the spectrum of waves ! if (pert == 3) then ! of velocity ! if (abs(vper) > 0.0d+00) then pot(:,:,:,:) = 0.0d+00 do n = 1, nper do j = 1, nn yp = sqrt(max(0.0d+00, 1.0d+00 - tanh(yc(j) / dlta)**2)) do i = 1, nn pot(3,i,j,:) = pot(3,i,j,:) & + vper / kx(n) * cos(kx(n) * xc(i) + ph(n)) * yp end do ! i = 1, nn end do ! j = 1, nn end do ! n = 1, nper call curl(dh(1:3), pot(:,:,:,:), qpert(ivx:ivz,:,:,:)) end if ! vper /= 0.0 ! of magnetic field ! if (abs(bper) > 0.0d+00) then pot(:,:,:,:) = 0.0d+00 do n = 1, nper do j = 1, nn yp = sqrt(max(0.0d+00, 1.0d+00 - tanh(yc(j) / dlta)**2)) do i = 1, nn pot(3,i,j,:) = pot(3,i,j,:) & + bper / kx(n) * cos(kx(n) * xc(i) + ph(n)) * yp end do ! i = 1, nn end do ! j = 1, nn end do ! n = 1, nper call curl(dh(1:3), pot(:,:,:,:), qpert(ibx:ibz,:,:,:)) end if ! bper /= 0.0 end if ! pert == 3 ! load the external perturbation ! if (pert == 4) then write(filename, "('perturbations-',i2.2,'-',i2.2,'.dat')") & nc, pdata%meta%level open(newunit=io, file=filename) read(io, *, iostat=status) kvec read(io, *, iostat=status) n read(io, *, iostat=status) ypl allocate(yt(n), qr(nv,n), qi(nv,n)) qr = 0.0d+00 qi = 0.0d+00 do p = 1, n read(io, *, iostat=status) yt(p), qr(ivx, p), qi(ivx, p), & qr(ivy, p), qi(ivy, p), & qr(ibx, p), qi(ibx, p), & qr(iby, p), qi(iby, p) if (status /= 0) exit end do close(io) jl = (pdata%meta%ymin - ypl) / dh(2) - ng + 1 ju = jl + nn - 1 do p = jl, ju j = p - jl + 1 do i = 1, nn qpert(ivx,i,j,:) = qr(ivx,p) * cos(kvec * xc(i)) & + qi(ivx,p) * sin(kvec * xc(i)) qpert(ivy,i,j,:) = qr(ivy,p) * cos(kvec * xc(i)) & + qi(ivy,p) * sin(kvec * xc(i)) qpert(ibx,i,j,:) = qr(ibx,p) * cos(kvec * xc(i)) & + qi(ibx,p) * sin(kvec * xc(i)) qpert(iby,i,j,:) = qr(iby,p) * cos(kvec * xc(i)) & + qi(iby,p) * sin(kvec * xc(i)) end do end do deallocate(yt, qr, qi) end if ! pert == 4 ! iterate over all positions in the XZ plane ! #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do i = 1, nn ! set the variable profiles with perturbations ! q(:,:) = qprof(:,:) + qpert(:,i,:,k) ! convert the primitive variables to the conservative ones ! call prim2cons(q(:,:), u(:,:), .true.) ! copy the primitive variables to the current block ! pdata%q(:,i,:,k) = q(:,:) ! copy the conserved variables to the current block ! pdata%u(:,i,:,k) = u(:,:) end do ! i = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ work_in_use(nt) = .false. 100 continue !------------------------------------------------------------------------------- ! end subroutine setup_user_problem ! !=============================================================================== ! ! subroutine UPDATE_USER_SOURCES: ! ------------------------------ ! ! Subroutine adds the user defined source terms. ! ! Arguments: ! ! t, dt - the time and time increment; ! pdata - the pointer to a data block; ! !=============================================================================== ! subroutine update_user_sources(t, dt, pdata) use blocks, only : block_data implicit none real(kind=8) , intent(in) :: t, dt type(block_data), pointer, intent(inout) :: pdata !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- ! end subroutine update_user_sources ! !=============================================================================== ! ! subroutine UPDATE_USER_SHAPES: ! ----------------------------- ! ! Subroutine defines the regions updated by user. ! ! Arguments: ! ! pdata - pointer to the data block structure of the currently initialized ! block; ! time - time at the moment of update; ! dt - time step since the last update; ! !=============================================================================== ! subroutine update_user_shapes(pdata, time, dt) use blocks , only : block_data use constants , only : pi use coordinates, only : nn => bcells use coordinates, only : ay, adx, ady #if NDIMS == 3 use coordinates, only : adz #endif /* NDIMS == 3 */ 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 type(block_data), pointer, intent(inout) :: pdata real(kind=8) , intent(in) :: time, dt logical, save :: first = .true. integer :: j, k, status real(kind=8) :: fl, fr, fd, fa, fb integer, save :: nt !$ integer :: omp_get_thread_num real(kind=8), dimension(3) :: dh real(kind=8), dimension(nn) :: yc, fc real(kind=8), dimension(nv,nn) :: q, u real(kind=8), dimension(:,:,:), pointer, save :: q2 !$omp threadprivate(first, nt, q2) character(len=*), parameter :: loc = 'USER_PROBLEM:update_user_shapes()' !------------------------------------------------------------------------------- ! #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ nt = 0 !$ nt = omp_get_thread_num() if (first) then j = nn**NDIMS call resize_workspace(j, status) if (status /= 0) then call print_message(loc, "Could not resize the workspace!") go to 100 end if #if NDIMS == 3 q2(1:nn,1:nn,1:nn) => work(1:j,nt) #else /* NDIMS == 3 */ q2(1:nn,1:nn,1: 1) => work(1:j,nt) #endif /* NDIMS == 3 */ first = .false. end if if (work_in_use(nt)) & call print_message(loc, "Workspace is being used right now! " // & "Corruptions can occur!") work_in_use(nt) = .true. yc(:) = pdata%meta%ymin + ay(pdata%meta%level,:) if (yc(1) < -yabs .or. yc(nn) > yabs) then dh(1) = adx(pdata%meta%level) dh(2) = ady(pdata%meta%level) #if NDIMS == 3 dh(3) = adz(pdata%meta%level) #endif /* NDIMS == 3 */ fa = 1.0d+00 - exp(- dt / tabs) fb = 5.0d-01 * adif * minval(dh(1:NDIMS))**2 call laplace(dh, pdata%q(ivy,:,:,:), q2) fc(:) = (abs(yc(:)) - yabs) / adec do j = 1, nn if (fc(j) > 0.0d+00) then if (fc(j) < 1.0d+00) then fr = fa * sin(5.0d-01 * pi * fc(j))**2 else fr = fa end if fl = 1.0d+00 - fr fd = fr * fb #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ q(idn,:) = fl * pdata%q(idn,:,j,k) + fr * dens q(ivx,:) = fl * pdata%q(ivx,:,j,k) q(ivy,:) = pdata%q(ivy,:,j,k) + fd * q2(:,j,k) q(ivz,:) = fl * pdata%q(ivz,:,j,k) q(ibx,:) = fl * pdata%q(ibx,:,j,k) + fr * sign(bamp, yc(j)) q(iby,:) = fl * pdata%q(iby,:,j,k) + fr * btra q(ibz,:) = fl * pdata%q(ibz,:,j,k) + fr * bgui q(ibp,:) = fl * pdata%q(ibp,:,j,k) if (ipr > 0) q(ipr,:) = fl * pdata%q(ipr,:,j,k) + fr * pres call prim2cons(q(:,:), u(:,:)) pdata%q(:,:,j,k) = q(:,:) pdata%u(:,:,j,k) = u(:,:) #if NDIMS == 3 end do #endif /* NDIMS == 3 */ end if end do end if work_in_use(nt) = .false. 100 continue !------------------------------------------------------------------------------- ! end subroutine update_user_shapes ! !=============================================================================== ! ! subroutine USER_BOUNDARY_X: ! -------------------------- ! ! Subroutine updates ghost zones within the specific region along ! the X direction. ! ! Arguments: ! ! is - the block side along the X direction for the ghost zone update; ! jl, ju - the cell index limits for the Y direction; ! kl, ku - the cell index limits for the Z direction; ! t, dt - time and time increment; ! x, y, z - the block coordinates; ! qn - the array of variables to update; ! !=============================================================================== ! subroutine user_boundary_x(is, jl, ju, kl, ku, t, dt, x, y, z, qn) use coordinates , only : nn => bcells, nb, ne, nbl, neu use equations , only : magnetized, ivx, ibx, iby, ibp #if NDIMS == 3 use equations , only : ibz #endif /* NDIMS == 3 */ implicit none integer , intent(in) :: is, jl, ju, kl, ku real(kind=8) , intent(in) :: t, dt real(kind=8), dimension(:) , intent(in) :: x, y, z real(kind=8), dimension(:,:,:,:), intent(inout) :: qn integer :: im2, im1, i , ip1, ip2 integer :: jm2, jm1, j , jp1, jp2 integer :: k #if NDIMS == 3 integer :: km2, km1, kp1, kp2 #endif /* NDIMS == 3 */ real(kind=8) :: dx, dy, dxy #if NDIMS == 3 real(kind=8) :: dz, dxz #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! if (magnetized) then #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ dx = x(2) - x(1) dy = y(2) - y(1) #if NDIMS == 3 dz = z(2) - z(1) #endif /* NDIMS == 3 */ dxy = dx / dy #if NDIMS == 3 dxz = dx / dz #endif /* NDIMS == 3 */ if (is == 1) then ! iterate over left-side ghost layers ! do i = nbl, 1, -1 ! calculate neighbor cell indices ! ip1 = min(nn, i + 1) ip2 = min(nn, i + 2) ! iterate over boundary layer ! #if NDIMS == 3 do k = kl, ku km2 = max( 1, k - 2) km1 = max( 1, k - 1) kp1 = min(nn, k + 1) kp2 = min(nn, k + 2) #endif /* NDIMS == 3 */ do j = jl, ju jm2 = max( 1, j - 2) jm1 = max( 1, j - 1) jp1 = min(nn, j + 1) jp2 = min(nn, j + 2) ! make the normal derivative zero ! qn(:,i,j,k) = qn(:,nb,j,k) ! prevent the inflow ! qn(ivx,i,j,k) = min(0.0d+00, qn(ivx,nb,j,k)) ! update the normal component of magnetic field from divergence-free condition ! qn(ibx,i,j,k) = qn(ibx,ip2,j,k) & + (qn(iby,ip1,jp1,k) - qn(iby,ip1,jm1,k)) * dxy #if NDIMS == 3 qn(ibx,i,j,k) = qn(ibx,i ,j,k) & + (qn(ibz,ip1,j,kp1) - qn(ibz,ip1,j,km1)) * dxz #endif /* NDIMS == 3 */ qn(ibp,i,j,k) = 0.0d+00 end do ! j = jl, ju #if NDIMS == 3 end do ! k = kl, ku #endif /* NDIMS == 3 */ end do ! i = nbl, 1, -1 else ! is == 1 ! iterate over right-side ghost layers ! do i = neu, nn ! calculate neighbor cell indices ! im1 = max( 1, i - 1) im2 = max( 1, i - 2) ! iterate over boundary layer ! #if NDIMS == 3 do k = kl, ku km1 = max( 1, k - 1) kp1 = min(nn, k + 1) km2 = max( 1, k - 2) kp2 = min(nn, k + 2) #endif /* NDIMS == 3 */ do j = jl, ju jm1 = max( 1, j - 1) jp1 = min(nn, j + 1) jm2 = max( 1, j - 2) jp2 = min(nn, j + 2) ! make the normal derivative zero ! qn(:,i,j,k) = qn(:,ne,j,k) ! prevent the inflow ! qn(ivx,i,j,k) = max(0.0d+00, qn(ivx,ne,j,k)) ! update the normal component of magnetic field from divergence-free condition ! qn(ibx,i,j,k) = qn(ibx,im2,j,k) & + (qn(iby,im1,jm1,k) - qn(iby,im1,jp1,k)) * dxy #if NDIMS == 3 qn(ibx,i,j,k) = qn(ibx,i ,j,k) & + (qn(ibz,im1,j,km1) - qn(ibz,im1,j,kp1)) * dxz #endif /* NDIMS == 3 */ qn(ibp,i,j,k) = 0.0d+00 end do ! j = jl, ju #if NDIMS == 3 end do ! k = kl, ku #endif /* NDIMS == 3 */ end do ! i = neu, nn end if ! is == 1 else ! ibx > 0 if (is == 1) then do i = nbl, 1, -1 #if NDIMS == 3 qn( : ,i,jl:ju,kl:ku) = qn( : ,nb,jl:ju,kl:ku) qn(ivx,i,jl:ju,kl:ku) = min(qn(ivx,nb,jl:ju,kl:ku), 0.0d+00) #else /* NDIMS == 3 */ qn( : ,i,jl:ju, : ) = qn( : ,nb,jl:ju, : ) qn(ivx,i,jl:ju, : ) = min(qn(ivx,nb,jl:ju, : ), 0.0d+00) #endif /* NDIMS == 3 */ end do ! i = nbl, 1, -1 else do i = neu, nn #if NDIMS == 3 qn( : ,i,jl:ju,kl:ku) = qn( : ,ne,jl:ju,kl:ku) qn(ivx,i,jl:ju,kl:ku) = max(qn(ivx,ne,jl:ju,kl:ku), 0.0d+00) #else /* NDIMS == 3 */ qn( : ,i,jl:ju, : ) = qn( : ,ne,jl:ju, : ) qn(ivx,i,jl:ju, : ) = max(qn(ivx,ne,jl:ju, : ), 0.0d+00) #endif /* NDIMS == 3 */ end do ! i = neu, nn end if end if ! magnetized !------------------------------------------------------------------------------- ! end subroutine user_boundary_x ! !=============================================================================== ! ! subroutine USER_BOUNDARY_Y: ! -------------------------- ! ! Subroutine updates ghost zones within the specific region along ! the Y direction. ! ! Arguments: ! ! js - the block side along the Y direction for the ghost zone update; ! il, iu - the cell index limits for the X direction; ! kl, ku - the cell index limits for the Z direction; ! t, dt - time and time increment; ! x, y, z - the block coordinates; ! qn - the array of variables to update; ! !=============================================================================== ! subroutine user_boundary_y(js, il, iu, kl, ku, t, dt, x, y, z, qn) use coordinates, only : nn => bcells, nb, ne, nbl, neu use equations , only : magnetized, nv, ivy, qpbnd implicit none integer , intent(in) :: js, il, iu, kl, ku real(kind=8) , intent(in) :: t, dt real(kind=8), dimension(:) , intent(in) :: x, y, z real(kind=8), dimension(:,:,:,:), intent(inout) :: qn integer :: i, j, k, jm, jp real(kind=8) :: fl, fr !------------------------------------------------------------------------------- ! if (magnetized) then #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ fl = exp(- dt / tdec) fr = 1.0d+00 - fl if (js == 1) then do j = nbl, 1, -1 jp = j + 1 #if NDIMS == 3 do k = kl, ku #endif /* NDIMS == 3 */ do i = il, iu qn( : ,i,j,k) = fl * qn( : ,i,jp,k) + fr * qpbnd(:,2,1) qn(ivy,i,j,k) = qn(ivy,i,jp,k) end do ! i = il, iu #if NDIMS == 3 end do ! k = kl, ku #endif /* NDIMS == 3 */ end do ! j = nbl, 1, -1 else ! js = 1 do j = neu, nn jm = j - 1 #if NDIMS == 3 do k = kl, ku #endif /* NDIMS == 3 */ do i = il, iu qn( : ,i,j,k) = fl * qn( : ,i,jm,k) + fr * qpbnd(:,2,2) qn(ivy,i,j,k) = qn(ivy,i,jm,k) end do ! i = il, iu #if NDIMS == 3 end do ! k = kl, ku #endif /* NDIMS == 3 */ end do ! j = neu, nn end if ! js = 1 else ! magnetized if (js == 1) then do j = nbl, 1, -1 #if NDIMS == 3 qn(1:nv,il:iu,j,kl:ku) = qn(1:nv,il:iu,nb,kl:ku) qn(ivy ,il:iu,j,kl:ku) = min(0.0d+00, qn(ivy,il:iu,nb,kl:ku)) #else /* NDIMS == 3 */ qn(1:nv,il:iu,j, : ) = qn(1:nv,il:iu,nb, : ) qn(ivy ,il:iu,j, : ) = min(0.0d+00, qn(ivy,il:iu,nb, : )) #endif /* NDIMS == 3 */ end do ! j = nbl, 1, -1 else do j = neu, nn #if NDIMS == 3 qn(1:nv,il:iu,j,kl:ku) = qn(1:nv,il:iu,ne,kl:ku) qn(ivy ,il:iu,j,kl:ku) = max(0.0d+00, qn(ivy,il:iu,ne,kl:ku)) #else /* NDIMS == 3 */ qn(1:nv,il:iu,j, : ) = qn(1:nv,il:iu,ne, : ) qn(ivy ,il:iu,j, : ) = max(0.0d+00, qn(ivy,il:iu,ne, : )) #endif /* NDIMS == 3 */ end do ! j = neu, nn end if end if ! magnetized !------------------------------------------------------------------------------- ! end subroutine user_boundary_y ! !=============================================================================== ! ! subroutine USER_BOUNDARY_Z: ! -------------------------- ! ! Subroutine updates ghost zones within the specific region along ! the Z direction. ! ! Arguments: ! ! ks - the block side along the Z direction for the ghost zone update; ! il, iu - the cell index limits for the X direction; ! jl, ju - the cell index limits for the Y direction; ! t, dt - time and time increment; ! x, y, z - the block coordinates; ! qn - the array of variables to update; ! !=============================================================================== ! subroutine user_boundary_z(ks, il, iu, jl, ju, t, dt, x, y, z, qn) implicit none integer , intent(in) :: ks, il, iu,jl, ju real(kind=8) , intent(in) :: t, dt real(kind=8), dimension(:) , intent(in) :: x, y, z real(kind=8), dimension(:,:,:,:), intent(inout) :: qn !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- ! end subroutine user_boundary_z ! !=============================================================================== ! ! subroutine USER_STATISTICS: ! ------------------------------- ! ! Subroutine can be use to store user defined time statistics. The file to ! store these statistics should be properly created in subroutine ! initialize_user_problem() and closed in finalize_user_problem(). ! ! Arguments: ! ! time - the current simulation time; ! !=============================================================================== ! subroutine user_statistics(time) use blocks , only : block_meta, block_data, list_data use coordinates, only : nn => bcells, nb, ne, nbm, nep, periodic use coordinates, only : xmin, xmax 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 #endif /* MPI */ use operators , only : curl, gradient use workspace , only : resize_workspace, work, work_in_use implicit none real(kind=8), intent(in) :: time logical, save :: first = .true. type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata integer, save :: nt !$ integer :: omp_get_thread_num integer :: ni, nl, nu, nm, np, status real(kind=8) :: dvol, dyz, dxz real(kind=8), dimension(3) :: dh real(kind=8), dimension(nn) :: yc real(kind=8), dimension(32) :: rterms real(kind=8), dimension(:,:,:,:), pointer, save :: cur !$omp threadprivate(first, nt, cur) character(len=*), parameter :: loc = 'USER_PROBLEM:user_time_statistics()' !------------------------------------------------------------------------------- ! nt = 0 !$ nt = omp_get_thread_num() if (first) then ni = 3 * nn**NDIMS call resize_workspace(ni, status) if (status /= 0) then call print_message(loc, "Could not resize the workspace!") go to 100 end if #if NDIMS == 3 cur(1:3,1:nn,1:nn,1:nn) => work(1:ni,nt) #else /* NDIMS == 3 */ cur(1:3,1:nn,1:nn,1: 1) => work(1:ni,nt) #endif /* NDIMS == 3 */ first = .false. end if rterms(:) = 0.0d+00 if (work_in_use(nt)) & call print_message(loc, & "Workspace is being used right now! Corruptions can occur!") work_in_use(nt) = .true. pdata => list_data do while(associated(pdata)) pmeta => pdata%meta dh(1) = adx(pmeta%level) dh(2) = ady(pmeta%level) dh(3) = adz(pmeta%level) dvol = advol(pmeta%level) dyz = ady(pmeta%level) * adz(pmeta%level) dxz = adx(pmeta%level) * adz(pmeta%level) yc(:) = pmeta%ymin + ay(pmeta%level,:) ni = -1 nl = nb nu = ne do while (yc(nl) < ylo .and. nl < ne) nl = nl + 1 end do do while (yc(nu) > yup .and. nu > nb) nu = nu - 1 end do if (nl < nu) then ! calculate current density (J = ∇xB) ! call curl(dh(:), pdata%q(ibx:ibz,:,:,:), cur(1:3,:,:,:)) ! the integral of |Bx| ! rterms(1) = rterms(1) & #if NDIMS == 3 + sum(abs(pdata%q(ibx,nb:ne,nl:nu,nb:ne))) * dvol #else /* NDIMS == 3 */ + sum(abs(pdata%q(ibx,nb:ne,nl:nu, : ))) * dvol #endif /* NDIMS == 3 */ ! the integral of magnetic energy ! rterms(12) = rterms(12) & #if NDIMS == 3 + sum(pdata%u(ibx:ibz,nb:ne,nl:nu,nb:ne)**2) * dvol #else /* NDIMS == 3 */ + sum(pdata%u(ibx:ibz,nb:ne,nl:nu, : )**2) * dvol #endif /* NDIMS == 3 */ if (.not. periodic(1)) then if (pmeta%xmin <= xmin) then ! advection of magnetic energy across the lower X-boundary ! rterms(14) = rterms(14) & #if NDIMS == 3 + sum(sum(pdata%q(ibx:ibz,nbm:nb,nl:nu,nb:ne)**2,1) & * pdata%q(ivx ,nbm:nb,nl:nu,nb:ne)) * dyz #else /* NDIMS == 3 */ + sum(sum(pdata%q(ibx:ibz,nbm:nb,nl:nu, : )**2,1) & * pdata%q(ivx ,nbm:nb,nl:nu, : )) * dyz #endif /* NDIMS == 3 */ rterms(17) = rterms(17) & #if NDIMS == 3 - sum(sum(pdata%q(ivx:ivz,nbm:nb,nl:nu,nb:ne) & * pdata%q(ibx:ibz,nbm:nb,nl:nu,nb:ne),1) & * pdata%q(ibx ,nbm:nb,nl:nu,nb:ne)) * dyz #else /* NDIMS == 3 */ - sum(sum(pdata%q(ivx:ivz,nbm:nb,nl:nu, : ) & * pdata%q(ibx:ibz,nbm:nb,nl:nu, : ),1) & * pdata%q(ibx ,nbm:nb,nl:nu, : )) * dyz #endif /* NDIMS == 3 */ if (resistive) then ! diffusion of magnetic energy across the lower X-boundary ! rterms(20) = rterms(20) & #if NDIMS == 3 - sum(cur(2,nbm:nb,nl:nu,nb:ne) & * pdata%q(ibz,nbm:nb,nl:nu,nb:ne) & - cur(3,nbm:nb,nl:nu,nb:ne) & * pdata%q(iby,nbm:nb,nl:nu,nb:ne)) * dyz #else /* NDIMS == 3 */ - sum(cur(2,nbm:nb,nl:nu, : ) & * pdata%q(ibz,nbm:nb,nl:nu, : ) & - cur(3,nbm:nb,nl:nu, : ) & * pdata%q(iby,nbm:nb,nl:nu, : )) * dyz #endif /* NDIMS == 3 */ end if ! resistivity end if if (pmeta%xmax >= xmax) then ! advection of magnetic energy across the upper X-boundary ! rterms(14) = rterms(14) & #if NDIMS == 3 - sum(sum(pdata%q(ibx:ibz,ne:nep,nl:nu,nb:ne)**2,1) & * pdata%q(ivx ,ne:nep,nl:nu,nb:ne)) * dyz #else /* NDIMS == 3 */ - sum(sum(pdata%q(ibx:ibz,ne:nep,nl:nu, : )**2,1) & * pdata%q(ivx ,ne:nep,nl:nu, : )) * dyz #endif /* NDIMS == 3 */ rterms(17) = rterms(17) & #if NDIMS == 3 + sum(sum(pdata%q(ivx:ivz,ne:nep,nl:nu,nb:ne) & * pdata%q(ibx:ibz,ne:nep,nl:nu,nb:ne),1) & * pdata%q(ibx ,ne:nep,nl:nu,nb:ne)) * dyz #else /* NDIMS == 3 */ + sum(sum(pdata%q(ivx:ivz,ne:nep,nl:nu, : ) & * pdata%q(ibx:ibz,ne:nep,nl:nu, : ),1) & * pdata%q(ibx ,ne:nep,nl:nu, : )) * dyz #endif /* NDIMS == 3 */ if (resistive) then ! diffusion of magnetic energy across the upper X-boundary ! rterms(20) = rterms(20) & #if NDIMS == 3 + sum(cur(2,ne:nep,nl:nu,nb:ne) & * pdata%q(ibz,ne:nep,nl:nu,nb:ne) & - cur(3,ne:nep,nl:nu,nb:ne) & * pdata%q(iby,ne:nep,nl:nu,nb:ne)) * dyz #else /* NDIMS == 3 */ + sum(cur(2,ne:nep,nl:nu, : ) & * pdata%q(ibz,ne:nep,nl:nu, : ) & - cur(3,ne:nep,nl:nu, : ) & * pdata%q(iby,ne:nep,nl:nu, : )) * dyz #endif /* NDIMS == 3 */ end if ! resistivity end if end if if (pmeta%ymin <= ylo .and. ylo < pmeta%ymax) then ni = nl nm = ni - 1 ! get the inflow speed ! rterms(10) = rterms(10) & #if NDIMS == 3 + sum(pdata%q(ivy,nb:ne,ni,nb:ne)) * dxz #else /* NDIMS == 3 */ + sum(pdata%q(ivy,nb:ne,ni, : )) * dxz #endif /* NDIMS == 3 */ ! mean Bx at boundary ! rterms(2) = rterms(2) & #if NDIMS == 3 + sum(abs(pdata%q(ibx,nb:ne,nm:ni,nb:ne))) * dxz #else /* NDIMS == 3 */ + sum(abs(pdata%q(ibx,nb:ne,nm:ni, : ))) * dxz #endif /* NDIMS == 3 */ ! advection of Bx along Y ! rterms(3) = rterms(3) & #if NDIMS == 3 + sum(abs(pdata%q(ibx,nb:ne,nm:ni,nb:ne)) & * pdata%q(ivy,nb:ne,nm:ni,nb:ne)) * dxz #else /* NDIMS == 3 */ + sum(abs(pdata%q(ibx,nb:ne,nm:ni, : )) & * pdata%q(ivy,nb:ne,nm:ni, : )) * dxz #endif /* NDIMS == 3 */ ! shear of By along X ! rterms(5) = rterms(5) & #if NDIMS == 3 - sum(sign(pdata%q(iby,nb:ne,nm:ni,nb:ne) & * pdata%q(ivx,nb:ne,nm:ni,nb:ne), & pdata%q(ibx,nb:ne,nm:ni,nb:ne))) * dxz #else /* NDIMS == 3 */ - sum(sign(pdata%q(iby,nb:ne,nm:ni, : ) & * pdata%q(ivx,nb:ne,nm:ni, : ), & pdata%q(ibx,nb:ne,nm:ni, : ))) * dxz #endif /* NDIMS == 3 */ ! mean magnetic energy ! rterms(13) = rterms(13) & #if NDIMS == 3 + sum(pdata%q(ibx:ibz,nb:ne,nm:ni,nb:ne)**2) * dxz #else /* NDIMS == 3 */ + sum(pdata%q(ibx:ibz,nb:ne,nm:ni, : )**2) * dxz #endif /* NDIMS == 3 */ ! advection of magnetic energy across the lower Y-boundary ! rterms(15) = rterms(15) & #if NDIMS == 3 + sum(sum(pdata%q(ibx:ibz,nb:ne,nm:ni,nb:ne)**2,1) & * pdata%q(ivy ,nb:ne,nm:ni,nb:ne)) * dxz #else /* NDIMS == 3 */ + sum(sum(pdata%q(ibx:ibz,nb:ne,nm:ni, : )**2,1) & * pdata%q(ivy ,nb:ne,nm:ni, : )) * dxz #endif /* NDIMS == 3 */ ! advection of magnetic energy ! rterms(18) = rterms(18) & #if NDIMS == 3 - sum(sum(pdata%q(ivx:ivz,nb:ne,nm:ni,nb:ne) & * pdata%q(ibx:ibz,nb:ne,nm:ni,nb:ne),1) & * pdata%q(iby ,nb:ne,nm:ni,nb:ne)) * dxz #else /* NDIMS == 3 */ - sum(sum(pdata%q(ivx:ivz,nb:ne,nm:ni, : ) & * pdata%q(ibx:ibz,nb:ne,nm:ni, : ),1) & * pdata%q(iby ,nb:ne,nm:ni, : )) * dxz #endif /* NDIMS == 3 */ if (resistive) then ! diffusion of Bx through ! rterms(7) = rterms(7) & #if NDIMS == 3 + sum(sign(cur(3,nb:ne,nm:ni,nb:ne), & pdata%q(ibx,nb:ne,nm:ni,nb:ne))) * dxz #else /* NDIMS == 3 */ + sum(sign(cur(3,nb:ne,nm:ni, : ), & pdata%q(ibx,nb:ne,nm:ni, : ))) * dxz #endif /* NDIMS == 3 */ ! diffusion of magnetic energy through the lower Y boundary ! rterms(21) = rterms(21) & #if NDIMS == 3 + sum(cur(3,nb:ne,nm:ni,nb:ne) & * pdata%q(ibx,nb:ne,nm:ni,nb:ne) & - cur(1,nb:ne,nm:ni,nb:ne) & * pdata%q(ibz,nb:ne,nm:ni,nb:ne)) * dxz #else /* NDIMS == 3 */ + sum(cur(3,nb:ne,nm:ni, : ) & * pdata%q(ibx,nb:ne,nm:ni, : ) & - cur(1,nb:ne,nm:ni, : ) & * pdata%q(ibz,nb:ne,nm:ni, : )) * dxz #endif /* NDIMS == 3 */ end if ! resistivity end if if (pmeta%ymin < yup .and. yup <= pmeta%ymax) then ni = nu np = ni + 1 ! get the inflow speed ! rterms(11) = rterms(11) & #if NDIMS == 3 - sum(pdata%q(ivy,nb:ne,ni,nb:ne)) * dxz #else /* NDIMS == 3 */ - sum(pdata%q(ivy,nb:ne,ni, : )) * dxz #endif /* NDIMS == 3 */ ! mean Bx at boundary ! rterms(2) = rterms(2) & #if NDIMS == 3 + sum(abs(pdata%q(ibx,nb:ne,ni:np,nb:ne))) * dxz #else /* NDIMS == 3 */ + sum(abs(pdata%q(ibx,nb:ne,ni:np, : ))) * dxz #endif /* NDIMS == 3 */ ! advection of Bx along Y ! rterms(3) = rterms(3) & #if NDIMS == 3 - sum(abs(pdata%q(ibx,nb:ne,ni:np,nb:ne)) & * pdata%q(ivy,nb:ne,ni:np,nb:ne)) * dxz #else /* NDIMS == 3 */ - sum(abs(pdata%q(ibx,nb:ne,ni:np, : )) & * pdata%q(ivy,nb:ne,ni:np, : )) * dxz #endif /* NDIMS == 3 */ ! shear of By along X ! rterms(5) = rterms(5) & #if NDIMS == 3 + sum(sign(pdata%q(iby,nb:ne,ni:np,nb:ne) & * pdata%q(ivx,nb:ne,ni:np,nb:ne), & pdata%q(ibx,nb:ne,ni:np,nb:ne))) * dxz #else /* NDIMS == 3 */ + sum(sign(pdata%q(iby,nb:ne,ni:np, : ) & * pdata%q(ivx,nb:ne,ni:np, : ), & pdata%q(ibx,nb:ne,ni:np, : ))) * dxz #endif /* NDIMS == 3 */ ! mean magnetic energy ! rterms(13) = rterms(13) & #if NDIMS == 3 + sum(pdata%q(ibx:ibz,nb:ne,ni:np,nb:ne)**2) * dxz #else /* NDIMS == 3 */ + sum(pdata%q(ibx:ibz,nb:ne,ni:np, : )**2) * dxz #endif /* NDIMS == 3 */ ! advection of magnetic energy ! rterms(15) = rterms(15) & #if NDIMS == 3 - sum(sum(pdata%q(ibx:ibz,nb:ne,ni:np,nb:ne)**2,1) & * pdata%q(ivy ,nb:ne,ni:np,nb:ne)) * dxz #else /* NDIMS == 3 */ - sum(sum(pdata%q(ibx:ibz,nb:ne,ni:np, : )**2,1) & * pdata%q(ivy ,nb:ne,ni:np, : )) * dxz #endif /* NDIMS == 3 */ rterms(18) = rterms(18) & #if NDIMS == 3 + sum(sum(pdata%q(ivx:ivz,nb:ne,ni:np,nb:ne) & * pdata%q(ibx:ibz,nb:ne,ni:np,nb:ne),1) & * pdata%q(iby ,nb:ne,ni:np,nb:ne)) * dxz #else /* NDIMS == 3 */ + sum(sum(pdata%q(ivx:ivz,nb:ne,ni:np, : ) & * pdata%q(ibx:ibz,nb:ne,ni:np, : ),1) & * pdata%q(iby ,nb:ne,ni:np, : )) * dxz #endif /* NDIMS == 3 */ if (resistive) then ! diffusion of Bx ! rterms(7) = rterms(7) & #if NDIMS == 3 - sum(sign(cur(3,nb:ne,ni:np,nb:ne), & pdata%q(ibx,nb:ne,ni:np,nb:ne))) * dxz #else /* NDIMS == 3 */ - sum(sign(cur(3,nb:ne,ni:np, : ), & pdata%q(ibx,nb:ne,ni:np, : ))) * dxz #endif /* NDIMS == 3 */ ! diffusion of magnetic energy ! rterms(21) = rterms(21) & #if NDIMS == 3 - sum(cur(3,nb:ne,ni:np,nb:ne) & * pdata%q(ibx,nb:ne,ni:np,nb:ne) & - cur(1,nb:ne,ni:np,nb:ne) & * pdata%q(ibz,nb:ne,ni:np,nb:ne)) * dxz #else /* NDIMS == 3 */ - sum(cur(3,nb:ne,ni:np, : ) & * pdata%q(ibx,nb:ne,ni:np, : ) & - cur(1,nb:ne,ni:np, : ) & * pdata%q(ibz,nb:ne,ni:np, : )) * dxz #endif /* NDIMS == 3 */ end if ! resistivity end if ! get the conversion between the magnetic energy and kinetic and ! internal energies ! rterms(23) = rterms(23) & #if NDIMS == 3 + sum((pdata%q(ivy,nb:ne,nl:nu,nb:ne) & * pdata%q(ibz,nb:ne,nl:nu,nb:ne) & - pdata%q(ivz,nb:ne,nl:nu,nb:ne) & * pdata%q(iby,nb:ne,nl:nu,nb:ne)) & * cur(1,nb:ne,nl:nu,nb:ne)) * dvol #else /* NDIMS == 3 */ + sum((pdata%q(ivy,nb:ne,nl:nu, : ) & * pdata%q(ibz,nb:ne,nl:nu, : ) & - pdata%q(ivz,nb:ne,nl:nu, : ) & * pdata%q(iby,nb:ne,nl:nu, : )) & * cur(1,nb:ne,nl:nu, : )) * dvol #endif /* NDIMS == 3 */ rterms(23) = rterms(23) & #if NDIMS == 3 + sum((pdata%q(ivz,nb:ne,nl:nu,nb:ne) & * pdata%q(ibx,nb:ne,nl:nu,nb:ne) & - pdata%q(ivx,nb:ne,nl:nu,nb:ne) & * pdata%q(ibz,nb:ne,nl:nu,nb:ne)) & * cur(2,nb:ne,nl:nu,nb:ne)) * dvol #else /* NDIMS == 3 */ + sum((pdata%q(ivz,nb:ne,nl:nu, : ) & * pdata%q(ibx,nb:ne,nl:nu, : ) & - pdata%q(ivx,nb:ne,nl:nu, : ) & * pdata%q(ibz,nb:ne,nl:nu, : )) & * cur(2,nb:ne,nl:nu, : )) * dvol #endif /* NDIMS == 3 */ rterms(23) = rterms(23) & #if NDIMS == 3 + sum((pdata%q(ivx,nb:ne,nl:nu,nb:ne) & * pdata%q(iby,nb:ne,nl:nu,nb:ne) & - pdata%q(ivy,nb:ne,nl:nu,nb:ne) & * pdata%q(ibx,nb:ne,nl:nu,nb:ne)) & * cur(3,nb:ne,nl:nu,nb:ne)) * dvol #else /* NDIMS == 3 */ + sum((pdata%q(ivx,nb:ne,nl:nu, : ) & * pdata%q(iby,nb:ne,nl:nu, : ) & - pdata%q(ivy,nb:ne,nl:nu, : ) & * pdata%q(ibx,nb:ne,nl:nu, : )) & * cur(3,nb:ne,nl:nu, : )) * dvol #endif /* NDIMS == 3 */ if (resistive) then rterms(24) = rterms(24) & #if NDIMS == 3 - sum(cur(1:3,nb:ne,nl:nu,nb:ne)**2) * dvol #else /* NDIMS == 3 */ - sum(cur(1:3,nb:ne,nl:nu, : )**2) * dvol #endif /* NDIMS == 3 */ end if ! calculate gradient (∇ψ) ! call gradient(dh(:), pdata%q(ibp,:,:,:), cur(1:3,:,:,:)) rterms(25) = rterms(25) & #if NDIMS == 3 - sum(pdata%q(ibx:ibz,nb:ne,nl:nu,nb:ne) & * cur(1:3,nb:ne,nl:nu,nb:ne)) * dvol #else /* NDIMS == 3 */ - sum(pdata%q(ibx:ibz,nb:ne,nl:nu, : ) & * cur(1:3,nb:ne,nl:nu, : )) * dvol #endif /* NDIMS == 3 */ ! contribution to |Bx| from ∇ψ ! rterms(9) = rterms(9) & #if NDIMS == 3 - sum(sign(cur(1,nb:ne,nl:nu,nb:ne), & pdata%q(ibx,nb:ne,nl:nu,nb:ne))) * dvol #else /* NDIMS == 3 */ - sum(sign(cur(1,nb:ne,nl:nu, : ), & pdata%q(ibx,nb:ne,nl:nu, : ))) * dvol #endif /* NDIMS == 3 */ end if ! nl < nu pdata => pdata%next end do work_in_use(nt) = .false. #ifdef MPI call reduce_sum(rterms(:)) #endif /* MPI */ rterms(2) = 2.50d-01 * rterms(2) / yarea rterms(12) = 5.00d-01 * rterms(12) rterms(13) = 1.25d-01 * rterms(13) / yarea rterms(3:6) = 5.00d-01 * rterms(3:6) rterms(14:19) = 5.00d-01 * rterms(14:19) rterms(7:8) = 5.00d-01 * eta * rterms(7:8) rterms(20:22) = 5.00d-01 * eta * rterms(20:22) rterms(24) = eta * rterms(24) write(runit,"(26es25.15e3)") time, rterms(1:25) call flush_and_sync(runit) 100 continue !------------------------------------------------------------------------------- ! end subroutine user_statistics ! !=============================================================================== ! ! subroutine LOG_COSH: ! ------------------- ! ! Function calculates the logarithm of the hyperbolic cosine, which is ! the result of the integration of tanh(x). Direct calculation using ! Fortran intrinsic subroutines fails for large values of x, therefore ! the logarithm of cosh is approximated as |x| + log(1/2) for ! |x| > threshold. ! ! Arguments: ! ! x - function argument; ! !=============================================================================== ! function log_cosh(x) result(y) implicit none real(kind=8), intent(in) :: x real(kind=8) :: y real(kind=8), parameter :: th = acosh(huge(x)), lh = log(0.5d+00) !------------------------------------------------------------------------------- ! if (abs(x) < th) then y = log(cosh(x)) else y = abs(x) + lh end if !------------------------------------------------------------------------------- ! end function log_cosh !=============================================================================== ! end module user_problem