diff --git a/sources/user_problem.F90 b/sources/user_problem.F90 index 6caf6a3..00bc3fd 100644 --- a/sources/user_problem.F90 +++ b/sources/user_problem.F90 @@ -74,6 +74,13 @@ module user_problem 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 + + integer(kind=4), save :: runit = 33 + + logical, save :: resistive = .false. + ! flag indicating if the gravitational source term is enabled ! logical, save :: user_gravity_enabled = .false. @@ -120,7 +127,7 @@ module user_problem use constants , only : pi #endif /* NDIMS == 3 */ use constants , only : pi2 - use coordinates, only : ng => nghosts, ady, xlen, zlen + use coordinates, only : ng => nghosts, ady, xlen, zlen, ymax use equations , only : magnetized, csnd, csnd2 use helpers , only : print_section, print_parameter use parameters , only : get_parameter @@ -139,9 +146,10 @@ module user_problem ! local variables ! - character(len=64) :: perturbation = "noise" + character(len=64) :: perturbation = "noise", append = "off", fname + logical :: flag integer :: n, kd - real(kind=8) :: thh, fc, ka + real(kind=8) :: thh, fc, ka, ydis = 9.00d+99 #if NDIMS == 3 real(kind=8) :: thv, tx, ty, tz, tt #endif /* NDIMS == 3 */ @@ -216,6 +224,8 @@ module user_problem write(*,*) end if status = 1 + else + resistive = .true. end if call get_parameter("dens", dens) if (dens <= 0.0d+00) then @@ -353,6 +363,64 @@ module user_problem end if end if ! status +! prepare file to store reconnection rate terms +! + if (status == 0) then + +! 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-new.dat')") + inquire(file = fname, exist = flag) + case default + write(fname, "('reconnection-new_',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 + + end if ! status + ! proceed if no errors ! if (status == 0) then @@ -370,7 +438,7 @@ module user_problem 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 (eta > 0.0d+00) then + if (resistive) then call print_parameter(verbose, '( ) S (Lundquist number)', lund) end if call print_parameter(verbose, '(*) delta (thickness)', dlta) @@ -388,6 +456,7 @@ module user_problem call print_parameter(verbose, '(*) ycut' , ycut) call print_parameter(verbose, '(*) xdec' , xdec) call print_parameter(verbose, '(*) ydec' , ydec) + call print_parameter(verbose, '(*) ydistance' , ydis) end if ! status #ifdef PROFILE @@ -435,6 +504,10 @@ module 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) @@ -1555,16 +1628,408 @@ module user_problem ! subroutine user_time_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 +#ifdef MPI + use mpitools , only : reduce_sum +#endif /* MPI */ + use operators , only : curl, gradient + implicit none real(kind=8), intent(in) :: time + type(block_meta), pointer :: pmeta + type(block_data), pointer :: pdata + + integer :: ni, nl, nu, nm, np + real(kind=8) :: dvol, dvolh, dxz + real(kind=8), dimension(3) :: dh + real(kind=8), dimension(nn) :: yc + real(kind=8), dimension(32) :: rterms + +#if NDIMS == 3 + real(kind=8), dimension(3,nn,nn,nn) :: cur +#else /* NDIMS == 3 */ + real(kind=8), dimension(3,nn,nn, 1) :: cur +#endif /* NDIMS == 3 */ + !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imt) #endif /* PROFILE */ + rterms(:) = 0.0d+00 + + 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) + dvolh = 5.0d-01 * dvol + + 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 + +! 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) * dvolh +#else /* NDIMS == 3 */ + + sum(pdata%u(ibx:ibz,nb:ne,nl:nu, : )**2) * dvolh +#endif /* NDIMS == 3 */ + + 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 +! + 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 */ + + pdata => pdata%next + end do + +#ifdef MPI + call reduce_sum(rterms(:)) +#endif /* MPI */ + + rterms(2) = 2.50d-01 * rterms(2) / yarea + 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) + #ifdef PROFILE call stop_timer(imt) #endif /* PROFILE */