From 4cb96ca2b3e044d5c3ff777167184881055a9e6d Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 4 Oct 2022 13:10:38 -0300 Subject: [PATCH] USER_PROBLEM: Implement flux tube reconnection problem. Signed-off-by: Grzegorz Kowal --- sources/user_problem.F90 | 720 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 718 insertions(+), 2 deletions(-) diff --git a/sources/user_problem.F90 b/sources/user_problem.F90 index ed4fc60..905d6a2 100644 --- a/sources/user_problem.F90 +++ b/sources/user_problem.F90 @@ -31,6 +31,23 @@ module user_problem implicit none + real(kind=8), save :: beta = 2.00d+00 + real(kind=8), save :: zeta = 0.00d+00 + real(kind=8), save :: dens = 1.00d+00 + real(kind=8), save :: bamp = 1.00d+00 + real(kind=8), save :: bgui = 1.00d+00 + real(kind=8), save :: dlta = 3.3333d-03 + 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 :: eta = 0.00d+00 + + integer(kind=4), save :: runit = 33 + + logical, save :: resistive = .false. + private public :: initialize_user_problem, finalize_user_problem public :: user_statistics @@ -58,6 +75,12 @@ module user_problem ! subroutine initialize_user_problem(problem, rcount, verbose, status) + use boundaries , only : custom_boundary_x, custom_boundary_y + use equations , only : adiabatic_index, csnd, csnd2, ipr + use helpers , only : print_section, print_parameter, print_message + use mesh , only : setup_problem + use parameters , only : get_parameter + implicit none character(len=64), intent(in) :: problem @@ -65,8 +88,126 @@ module user_problem logical , intent(in) :: verbose integer , intent(out) :: status + character(len=64) :: append = "off", fname + logical :: flag + + character(len=*), parameter :: & + loc = 'USER_PROBLEM::initialize_user_problem()' !------------------------------------------------------------------------------- ! + call get_parameter("dens", dens) + 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 + 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("delta", dlta) + if (dlta <= 0.0d+00) then + if (verbose) & + call print_message(loc, & + "The initial thickness must be positive (delta > 0.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 + + pmag = 0.5d+00 * (bamp**2 + bgui**2) + pres = beta * pmag + ptot = pres + pmag + if (ipr > 0) then + csnd2 = adiabatic_index * pres / dens + else + csnd2 = pres / dens + end if + csnd = sqrt(csnd2) + valf = sqrt(2.0d+00 * pmag / dens) + lund = valf / max(tiny(eta), eta) + +! determine if to append or create another file +! + call get_parameter("statistics_append", append) + select case(trim(append)) + case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") + write(fname, "('energy_balance.dat')") + inquire(file = fname, exist = flag) + case default + write(fname, "('energy_balance_',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,13a25)') 'time', 'Emag total', & + '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 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 (amplitude)' , bamp) + call print_parameter(verbose, '(*) bgui (guide field)' , bgui) + call print_parameter(verbose, '(*) delta (thickness)' , dlta) + call print_parameter(verbose, '( ) pres (thermal pres.)' , pres) + call print_parameter(verbose, '( ) pmag (magnetic pres.)', pmag) + call print_parameter(verbose, '( ) csnd (sound speed)' , csnd) + call print_parameter(verbose, '( ) Valf (Alfven speed)' , valf) + if (resistive) & + call print_parameter(verbose, '( ) S (Lundquist number)', lund) + + setup_problem => setup_user_problem + custom_boundary_x => user_boundary_x + custom_boundary_y => user_boundary_y + status = 0 !------------------------------------------------------------------------------- @@ -116,14 +257,67 @@ module user_problem ! subroutine setup_user_problem(pdata) - use blocks, only : block_data + use blocks , only : block_data + use constants , only : pi, pi2 + use coordinates, only : nn => bcells + use coordinates, only : ax, ay + use equations , only : prim2cons, csnd2 + use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp implicit none type(block_data), pointer, intent(inout) :: pdata + integer :: i, j, k + real(kind=8) :: bx, by, bz + real(kind=8), dimension(nn) :: x, y + real(kind=8), dimension(nn) :: snx, csx, sny, csy, thy, ch2 + !------------------------------------------------------------------------------- ! + x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) + y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) + + snx(:) = sin(pi * x(:)) + csx(:) = cos(pi * x(:)) + sny(:) = sin(pi2 * y(:)) + csy(:) = cos(pi2 * y(:)) + thy(:) = tanh(y(:) / dlta) + ch2(:) = cosh(y(:) / dlta)**2 + + do j = 1, nn + do i = 1, nn + bx = bamp * csx(i) * (csy(j) * thy(j) + sny(j) / (pi2 * dlta * ch2(j))) + by = bamp * 5.0d-01 * snx(i) * sny(j) * thy(j) + bz = sqrt(bgui**2 + zeta**2 * max(0.0d+00, bamp**2 - bx**2 - by**2)) + + pdata%q(ibx,i,j,:) = bx + pdata%q(iby,i,j,:) = by + pdata%q(ibz,i,j,:) = bz + end do + end do + pdata%q(ibp,:,:,:) = 0.0d+00 + + if (ipr > 0) then + pdata%q(idn,:,:,:) = dens + pdata%q(ipr,:,:,:) = ptot - 5.0d-01 * sum(pdata%q(ibx:ibz,:,:,:)**2,1) + else + pdata%q(idn,:,:,:) = (ptot - 5.0d-01 * sum(pdata%q(ibx:ibz,:,:,:)**2,1)) & + / csnd2 + end if + pdata%q(ivx,:,:,:) = 0.0d+00 + pdata%q(ivy,:,:,:) = 0.0d+00 + pdata%q(ivz,:,:,:) = 0.0d+00 + +#if NDIMS == 3 + do k = 1, nn +#else /* NDIMS == 3 */ + do k = 1, 1 +#endif /* NDIMS == 3 */ + do j = 1, nn + call prim2cons(pdata%q(:,:,j,k), pdata%u(:,:,j,k), .true.) + end do + end do !------------------------------------------------------------------------------- ! @@ -212,6 +406,9 @@ module user_problem ! 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 : ivx, ibx + implicit none integer , intent(in) :: is, jl, ju, kl, ku @@ -219,8 +416,31 @@ module user_problem real(kind=8), dimension(:) , intent(in) :: x, y, z real(kind=8), dimension(:,:,:,:), intent(inout) :: qn + integer :: i, j, ir + !------------------------------------------------------------------------------- ! + if (is == 1) then + do i = nbl, 1, -1 + ir = nb + (nbl - i) + do j = jl, ju + qn( : ,i,j,:) = qn( : ,ir,j,:) + + qn(ivx,i,j,:) = - qn(ivx,ir,j,:) + qn(ibx,i,j,:) = - qn(ibx,ir,j,:) + end do + end do + else ! is == 1 + do i = neu, nn + ir = ne - (i - neu) + do j = jl, ju + qn( : ,i,j,:) = qn( : ,ir,j,:) + + qn(ivx,i,j,:) = - qn(ivx,ir,j,:) + qn(ibx,i,j,:) = - qn(ibx,ir,j,:) + end do + end do + end if !------------------------------------------------------------------------------- ! @@ -247,6 +467,9 @@ module user_problem ! 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 : ivy, iby + implicit none integer , intent(in) :: js, il, iu, kl, ku @@ -254,8 +477,31 @@ module user_problem real(kind=8), dimension(:) , intent(in) :: x, y, z real(kind=8), dimension(:,:,:,:), intent(inout) :: qn + integer :: i, j, jr + !------------------------------------------------------------------------------- ! + if (js == 1) then + do j = nbl, 1, -1 + jr = nb + (nbl - j) + do i = il, iu + qn( : ,i,j,:) = qn( : ,i,jr,:) + + qn(ivy,i,j,:) = - qn(ivy,i,jr,:) + qn(iby,i,j,:) = - qn(iby,i,jr,:) + end do + end do + else ! js == 1 + do j = neu, nn + jr = ne - (j - neu) + do i = il, iu + qn( : ,i,j,:) = qn( : ,i,jr,:) + + qn(ivy,i,j,:) = - qn(ivy,i,jr,:) + qn(iby,i,j,:) = - qn(iby,i,jr,:) + end do + end do + end if !------------------------------------------------------------------------------- ! @@ -313,12 +559,482 @@ module user_problem ! subroutine user_statistics(time) - implicit none + use blocks , only : block_meta, block_data, list_data + use coordinates, only : nc => ncells, nn => bcells, nb, ne, nbm, nep +#if NDIMS == 3 + use coordinates, only : periodic, xmin, xmax, ymin, ymax, zmin, zmax +#else /* NDIMS == 3 */ + use coordinates, only : periodic, xmin, xmax, ymin, ymax +#endif /* NDIMS == 3 */ + use coordinates, only : adx, ady, adz, advol + 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 :: nl, nu, status +#if NDIMS == 3 + real(kind=8) :: dvol, dyz, dxz, dxy +#else /* NDIMS == 3 */ + real(kind=8) :: dvol, dyz, dxz +#endif /* NDIMS == 3 */ + real(kind=8), dimension(3) :: dh + real(kind=8), dimension(16) :: rterms + + real(kind=8), dimension(:,:,:,:), pointer, save :: cur + real(kind=8), dimension(:,:,:) , pointer, save :: vv, bb, jj +!$omp threadprivate(first, nt, cur, vv, bb, jj) + + character(len=*), parameter :: loc = 'USER_PROBLEM:user_time_statistics()' !------------------------------------------------------------------------------- ! + nt = 0 +!$ nt = omp_get_thread_num() + if (first) then + if (resistive) then + nu = 3 * nn**NDIMS + 9 * nc**(NDIMS - 1) + else + nu = 3 * nn**NDIMS + 6 * nc**(NDIMS - 1) + end if + + call resize_workspace(nu, status) + if (status /= 0) then + call print_message(loc, "Could not resize the workspace!") + go to 100 + end if + + nl = 1 + nu = nl - 1 + 3 * nn**NDIMS +#if NDIMS == 3 + cur(1:3,1:nn,1:nn,1:nn) => work(nl:nu,nt) +#else /* NDIMS == 3 */ + cur(1:3,1:nn,1:nn,1: 1) => work(nl:nu,nt) +#endif /* NDIMS == 3 */ + nl = nu + 1 + nu = nl - 1 + 3 * nc**(NDIMS - 1) +#if NDIMS == 3 + vv(1:3,1:nc,1:nc) => work(nl:nu,nt) +#else /* NDIMS == 3 */ + vv(1:3,1:nc,1: 1) => work(nl:nu,nt) +#endif /* NDIMS == 3 */ + nl = nu + 1 + nu = nl - 1 + 3 * nc**(NDIMS - 1) +#if NDIMS == 3 + bb(1:3,1:nc,1:nc) => work(nl:nu,nt) +#else /* NDIMS == 3 */ + bb(1:3,1:nc,1: 1) => work(nl:nu,nt) +#endif /* NDIMS == 3 */ + if (resistive) then + nl = nu + 1 + nu = nl - 1 + 3 * nc**(NDIMS - 1) +#if NDIMS == 3 + jj(1:3,1:nc,1:nc) => work(nl:nu,nt) +#else /* NDIMS == 3 */ + jj(1:3,1:nc,1: 1) => work(nl:nu,nt) +#endif /* NDIMS == 3 */ + end if + + 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) +#if NDIMS == 3 + dxy = adx(pmeta%level) * ady(pmeta%level) +#endif /* NDIMS == 3 */ + +! calculate current density (J = ∇xB) +! + call curl(dh(:), pdata%q(ibx:ibz,:,:,:), cur(1:3,:,:,:)) + +! the integral of magnetic energy +! + rterms(1) = rterms(1) & +#if NDIMS == 3 + + sum(pdata%u(ibx:ibz,nb:ne,nb:ne,nb:ne)**2) * dvol +#else /* NDIMS == 3 */ + + sum(pdata%u(ibx:ibz,nb:ne,nb:ne, : )**2) * dvol +#endif /* NDIMS == 3 */ + + if (.not. periodic(1)) then + + if (pmeta%xmin < xmin) then + +! interpolate velocity and magnetic field components +! +#if NDIMS == 3 + vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nbm:nb,nb:ne,nb:ne), 1) + vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nbm:nb,nb:ne,nb:ne), 1) + vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nbm:nb,nb:ne,nb:ne), 1) + bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nbm:nb,nb:ne,nb:ne), 1) + bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nbm:nb,nb:ne,nb:ne), 1) + bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nbm:nb,nb:ne,nb:ne), 1) +#else /* NDIMS == 3 */ + vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nbm:nb,nb:ne, : ), 1) + vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nbm:nb,nb:ne, : ), 1) + vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nbm:nb,nb:ne, : ), 1) + bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nbm:nb,nb:ne, : ), 1) + bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nbm:nb,nb:ne, : ), 1) + bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nbm:nb,nb:ne, : ), 1) +#endif /* NDIMS == 3 */ + +! advection of magnetic energy across the lower X-boundary +! + rterms(2) = rterms(2) + sum(sum(bb * bb, 1) * vv(1,:,:)) * dyz + +! generation of magnetic energy at the lower X-boundary +! + rterms(5) = rterms(5) - sum(sum(vv * bb, 1) * bb(1,:,:)) * dyz + + if (resistive) then + +! interpolate the components of current density +! +#if NDIMS == 3 + jj(2,:,:) = 5.0d-01 * sum(cur(2,nbm:nb,nb:ne,nb:ne), 1) + jj(3,:,:) = 5.0d-01 * sum(cur(3,nbm:nb,nb:ne,nb:ne), 1) +#else /* NDIMS == 3 */ + jj(2,:,:) = 5.0d-01 * sum(cur(2,nbm:nb,nb:ne, : ), 1) + jj(3,:,:) = 5.0d-01 * sum(cur(3,nbm:nb,nb:ne, : ), 1) +#endif /* NDIMS == 3 */ + +! diffusion of magnetic energy across the lower X-boundary +! + rterms(8) = rterms(8) & + + sum(jj(2,:,:) * bb(3,:,:) - jj(3,:,:) * bb(2,:,:)) * dyz + + end if ! resistivity + end if + + if (pmeta%xmax > xmax) then + +#if NDIMS == 3 + vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,ne:nep,nb:ne,nb:ne), 1) + vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,ne:nep,nb:ne,nb:ne), 1) + vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,ne:nep,nb:ne,nb:ne), 1) + bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,ne:nep,nb:ne,nb:ne), 1) + bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,ne:nep,nb:ne,nb:ne), 1) + bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,ne:nep,nb:ne,nb:ne), 1) +#else /* NDIMS == 3 */ + vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,ne:nep,nb:ne, : ), 1) + vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,ne:nep,nb:ne, : ), 1) + vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,ne:nep,nb:ne, : ), 1) + bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,ne:nep,nb:ne, : ), 1) + bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,ne:nep,nb:ne, : ), 1) + bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,ne:nep,nb:ne, : ), 1) +#endif /* NDIMS == 3 */ + +! advection of magnetic energy across the upper X-boundary +! + rterms(2) = rterms(2) - sum(sum(bb * bb, 1) * vv(1,:,:)) * dyz + +! generation of magnetic energy at the upper X-boundary +! + rterms(5) = rterms(5) + sum(sum(vv * bb, 1) * bb(1,:,:)) * dyz + + if (resistive) then + +! interpolate the components of current density +! +#if NDIMS == 3 + jj(2,:,:) = 5.0d-01 * sum(cur(2,ne:nep,nb:ne,nb:ne), 1) + jj(3,:,:) = 5.0d-01 * sum(cur(3,ne:nep,nb:ne,nb:ne), 1) +#else /* NDIMS == 3 */ + jj(2,:,:) = 5.0d-01 * sum(cur(2,ne:nep,nb:ne, : ), 1) + jj(3,:,:) = 5.0d-01 * sum(cur(3,ne:nep,nb:ne, : ), 1) +#endif /* NDIMS == 3 */ + +! diffusion of magnetic energy across the upper X-boundary +! + rterms(8) = rterms(8) & + - sum(jj(2,:,:) * bb(3,:,:) - jj(3,:,:) * bb(2,:,:)) * dyz + + end if ! resistivity + end if + end if ! non-periodic along X + + if (.not. periodic(2)) then + + if (pmeta%ymin < ymin) then + +#if NDIMS == 3 + vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nb:ne,nbm:nb,nb:ne), 1) + vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nb:ne,nbm:nb,nb:ne), 1) + vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nb:ne,nbm:nb,nb:ne), 1) + bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nb:ne,nbm:nb,nb:ne), 1) + bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nb:ne,nbm:nb,nb:ne), 1) + bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nb:ne,nbm:nb,nb:ne), 1) +#else /* NDIMS == 3 */ + vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nb:ne,nbm:nb, : ), 1) + vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nb:ne,nbm:nb, : ), 1) + vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nb:ne,nbm:nb, : ), 1) + bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nb:ne,nbm:nb, : ), 1) + bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nb:ne,nbm:nb, : ), 1) + bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nb:ne,nbm:nb, : ), 1) +#endif /* NDIMS == 3 */ + +! advection of magnetic energy across the lower Y-boundary +! + rterms(3) = rterms(3) + sum(sum(bb * bb, 1) * vv(2,:,:)) * dxz + +! generation of magnetic energy at the lower Y-boundary +! + rterms(6) = rterms(6) - sum(sum(vv * bb, 1) * bb(2,:,:)) * dxz + + if (resistive) then + +! interpolate the components of current density +! +#if NDIMS == 3 + jj(1,:,:) = 5.0d-01 * sum(cur(1,nb:ne,nbm:nb,nb:ne), 1) + jj(3,:,:) = 5.0d-01 * sum(cur(3,nb:ne,nbm:nb,nb:ne), 1) +#else /* NDIMS == 3 */ + jj(1,:,:) = 5.0d-01 * sum(cur(1,nb:ne,nbm:nb, : ), 1) + jj(3,:,:) = 5.0d-01 * sum(cur(3,nb:ne,nbm:nb, : ), 1) +#endif /* NDIMS == 3 */ + +! diffusion of magnetic energy through the lower Y boundary +! + rterms(9) = rterms(9) & + + sum(jj(3,:,:) * bb(1,:,:) - jj(1,:,:) * bb(3,:,:)) * dxz + + end if ! resistivity + end if + + if (ymax < pmeta%ymax) then + +#if NDIMS == 3 + vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nb:ne,ne:nep,nb:ne), 1) + vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nb:ne,ne:nep,nb:ne), 1) + vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nb:ne,ne:nep,nb:ne), 1) + bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nb:ne,ne:nep,nb:ne), 1) + bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nb:ne,ne:nep,nb:ne), 1) + bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nb:ne,ne:nep,nb:ne), 1) +#else /* NDIMS == 3 */ + vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nb:ne,ne:nep, : ), 1) + vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nb:ne,ne:nep, : ), 1) + vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nb:ne,ne:nep, : ), 1) + bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nb:ne,ne:nep, : ), 1) + bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nb:ne,ne:nep, : ), 1) + bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nb:ne,ne:nep, : ), 1) +#endif /* NDIMS == 3 */ + +! advection of magnetic energy across the upper Y-boundary +! + rterms(3) = rterms(3) - sum(sum(bb * bb, 1) * vv(2,:,:)) * dxz + +! generation of magnetic energy at the upper Y-boundary +! + rterms(6) = rterms(6) + sum(sum(vv * bb, 1) * bb(2,:,:)) * dxz + + if (resistive) then + +! interpolate the components of current density +! +#if NDIMS == 3 + jj(1,:,:) = 5.0d-01 * sum(cur(1,nb:ne,ne:nep,nb:ne), 1) + jj(3,:,:) = 5.0d-01 * sum(cur(3,nb:ne,ne:nep,nb:ne), 1) +#else /* NDIMS == 3 */ + jj(1,:,:) = 5.0d-01 * sum(cur(1,nb:ne,ne:nep, : ), 1) + jj(3,:,:) = 5.0d-01 * sum(cur(3,nb:ne,ne:nep, : ), 1) +#endif /* NDIMS == 3 */ + +! diffusion of magnetic energy +! + rterms(9) = rterms(9) & + - sum(jj(3,:,:) * bb(1,:,:) - jj(1,:,:) * bb(3,:,:)) * dxz + + end if ! resistivity + end if + end if ! non-periodic along Y + +#if NDIMS == 3 + if (.not. periodic(3)) then + + if (pmeta%zmin < zmin) then + +! interpolate velocity and magnetic field components +! + vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nb:ne,nb:ne,nbm:nb), 1) + vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nb:ne,nb:ne,nbm:nb), 1) + vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nb:ne,nb:ne,nbm:nb), 1) + bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nb:ne,nb:ne,nbm:nb), 1) + bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nb:ne,nb:ne,nbm:nb), 1) + bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nb:ne,nb:ne,nbm:nb), 1) + +! advection of magnetic energy across the lower Z-boundary +! + rterms(4) = rterms(4) + sum(sum(bb * bb, 1) * vv(3,:,:)) * dxy + +! generation of magnetic energy at the lower Z-boundary +! + rterms(7) = rterms(7) - sum(sum(vv * bb, 1) * bb(3,:,:)) * dxy + + if (resistive) then + +! interpolate the components of current density +! + jj(1,:,:) = 5.0d-01 * sum(cur(1,nb:ne,nb:ne,nbm:nb), 1) + jj(2,:,:) = 5.0d-01 * sum(cur(2,nb:ne,nb:ne,nbm:nb), 1) + +! diffusion of magnetic energy across the lower Z-boundary +! + rterms(10) = rterms(10) & + + sum(jj(1,:,:) * bb(2,:,:) - jj(2,:,:) * bb(1,:,:)) * dxy + + end if ! resistivity + end if + + if (pmeta%zmax > zmax) then + + vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,nb:ne,nb:ne,ne:nep), 1) + vv(2,:,:) = 5.0d-01 * sum(pdata%q(ivy,nb:ne,nb:ne,ne:nep), 1) + vv(3,:,:) = 5.0d-01 * sum(pdata%q(ivz,nb:ne,nb:ne,ne:nep), 1) + bb(1,:,:) = 5.0d-01 * sum(pdata%q(ibx,nb:ne,nb:ne,ne:nep), 1) + bb(2,:,:) = 5.0d-01 * sum(pdata%q(iby,nb:ne,nb:ne,ne:nep), 1) + bb(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nb:ne,nb:ne,ne:nep), 1) + +! advection of magnetic energy across the upper Z-boundary +! + rterms(4) = rterms(4) - sum(sum(bb * bb, 1) * vv(3,:,:)) * dxy + +! generation of magnetic energy at the upper Z-boundary +! + rterms(7) = rterms(7) + sum(sum(vv * bb, 1) * bb(3,:,:)) * dxy + + if (resistive) then + +! interpolate the components of current density +! + jj(1,:,:) = 5.0d-01 * sum(cur(1,nb:ne,nb:ne,ne:nep), 1) + jj(2,:,:) = 5.0d-01 * sum(cur(2,nb:ne,nb:ne,ne:nep), 1) + +! diffusion of magnetic energy across the upper Z-boundary +! + rterms(10) = rterms(10) & + - sum(jj(1,:,:) * bb(2,:,:) - jj(2,:,:) * bb(1,:,:)) * dxy + + end if ! resistivity + end if + end if ! non-periodic along Z +#endif /* NDIMS == 3 */ + +! conversion between magnetic and kinetic energies +! + rterms(11) = rterms(11) & +#if NDIMS == 3 + + sum((pdata%q(ivy,nb:ne,nb:ne,nb:ne) & + * pdata%q(ibz,nb:ne,nb:ne,nb:ne) & + - pdata%q(ivz,nb:ne,nb:ne,nb:ne) & + * pdata%q(iby,nb:ne,nb:ne,nb:ne)) & + * cur(1,nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + + sum((pdata%q(ivy,nb:ne,nb:ne, : ) & + * pdata%q(ibz,nb:ne,nb:ne, : ) & + - pdata%q(ivz,nb:ne,nb:ne, : ) & + * pdata%q(iby,nb:ne,nb:ne, : )) & + * cur(1,nb:ne,nb:ne, : )) * dvol +#endif /* NDIMS == 3 */ + rterms(11) = rterms(11) & +#if NDIMS == 3 + + sum((pdata%q(ivz,nb:ne,nb:ne,nb:ne) & + * pdata%q(ibx,nb:ne,nb:ne,nb:ne) & + - pdata%q(ivx,nb:ne,nb:ne,nb:ne) & + * pdata%q(ibz,nb:ne,nb:ne,nb:ne)) & + * cur(2,nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + + sum((pdata%q(ivz,nb:ne,nb:ne, : ) & + * pdata%q(ibx,nb:ne,nb:ne, : ) & + - pdata%q(ivx,nb:ne,nb:ne, : ) & + * pdata%q(ibz,nb:ne,nb:ne, : )) & + * cur(2,nb:ne,nb:ne, : )) * dvol +#endif /* NDIMS == 3 */ + rterms(11) = rterms(11) & +#if NDIMS == 3 + + sum((pdata%q(ivx,nb:ne,nb:ne,nb:ne) & + * pdata%q(iby,nb:ne,nb:ne,nb:ne) & + - pdata%q(ivy,nb:ne,nb:ne,nb:ne) & + * pdata%q(ibx,nb:ne,nb:ne,nb:ne)) & + * cur(3,nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + + sum((pdata%q(ivx,nb:ne,nb:ne, : ) & + * pdata%q(iby,nb:ne,nb:ne, : ) & + - pdata%q(ivy,nb:ne,nb:ne, : ) & + * pdata%q(ibx,nb:ne,nb:ne, : )) & + * cur(3,nb:ne,nb:ne, : )) * dvol +#endif /* NDIMS == 3 */ + +! conversion between magnetic energy and internal energies +! + if (resistive) then + rterms(12) = rterms(12) & +#if NDIMS == 3 + - sum(cur(1:3,nb:ne,nb:ne,nb:ne)**2) * dvol +#else /* NDIMS == 3 */ + - sum(cur(1:3,nb:ne,nb:ne, : )**2) * dvol +#endif /* NDIMS == 3 */ + end if + +! calculate gradient (∇ψ) +! + call gradient(dh(:), pdata%q(ibp,:,:,:), cur(1:3,:,:,:)) + + rterms(13) = rterms(13) & +#if NDIMS == 3 + - sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nb:ne) & + * cur(1:3,nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + - sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : ) & + * cur(1:3,nb:ne,nb:ne, : )) * dvol +#endif /* NDIMS == 3 */ + + pdata => pdata%next + end do + + work_in_use(nt) = .false. + +#ifdef MPI + call reduce_sum(rterms(:)) +#endif /* MPI */ + + rterms(1) = 5.00d-01 * rterms(1) + rterms(8:10) = eta * rterms(8:10) + rterms(12) = eta * rterms(12) + + write(runit,"(14es25.15e3)") time, rterms(1:13) + call flush_and_sync(runit) + + 100 continue !------------------------------------------------------------------------------- !