diff --git a/sources/user_problem.F90 b/sources/user_problem.F90 index d28f870..07cf65b 100644 --- a/sources/user_problem.F90 +++ b/sources/user_problem.F90 @@ -153,10 +153,10 @@ module user_problem 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')") + write(fname, "('reconnection.dat')") inquire(file = fname, exist = flag) case default - write(fname, "('energy_balance_',i2.2,'.dat')") rcount + write(fname, "('reconnection_',i2.2,'.dat')") rcount flag = .false. end select @@ -183,11 +183,11 @@ module user_problem ! 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,'("#",a24,9a25)') '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' write(runit,"('#')") call print_section(verbose, "Parameters (* - set, otherwise calculated)") @@ -572,7 +572,7 @@ module user_problem #ifdef MPI use mpitools , only : reduce_sum #endif /* MPI */ - use operators , only : curl, gradient + use operators , only : curl, derivative_1st use workspace , only : resize_workspace, work, work_in_use implicit none @@ -681,108 +681,11 @@ module user_problem ! rterms(1) = rterms(1) & #if NDIMS == 3 - + sum(pdata%u(ibx:ibz,nb:ne,nb:ne,nb:ne)**2) * dvol + + sum(abs(pdata%u(ibx,nb:ne,nb:ne,nb:ne))) * dvol #else /* NDIMS == 3 */ - + sum(pdata%u(ibx:ibz,nb:ne,nb:ne, : )**2) * dvol + + sum(abs(pdata%u(ibx,nb:ne,nb:ne, : ))) * 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 @@ -790,43 +693,43 @@ module user_problem #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 +! mean Bx at the Y-boundary ! - rterms(3) = rterms(3) + sum(sum(bb * bb, 1) * vv(2,:,:)) * dxz + rterms(2) = rterms(2) + sum(abs(bb(1,:,:))) * dxz -! generation of magnetic energy at the lower Y-boundary +! advection of Bx across the Y-boundary ! - rterms(6) = rterms(6) - sum(sum(vv * bb, 1) * bb(2,:,:)) * dxz + rterms(3) = rterms(3) & + + sum(sign(bb(1,:,:) * vv(2,:,:), bb(1,:,:))) * dxz + +! shear of By along the Y-boundary +! + rterms(5) = rterms(5) & + - sum(sign(bb(2,:,:) * vv(1,:,:), bb(1,:,:))) * dxz if (resistive) then -! interpolate the components of current density +! interpolate Jz ! #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 +! diffusion of Bx at the Y-boundary ! - rterms(9) = rterms(9) & - + sum(jj(3,:,:) * bb(1,:,:) - jj(1,:,:) * bb(3,:,:)) * dxz + rterms(7) = rterms(7) & + + sum(sign(jj(3,:,:), bb(1,:,:))) * dxz end if ! resistivity end if @@ -836,43 +739,43 @@ module user_problem #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 +! mean Bx at the boundary ! - rterms(3) = rterms(3) - sum(sum(bb * bb, 1) * vv(2,:,:)) * dxz + rterms(2) = rterms(2) + sum(abs(bb(1,:,:))) * dxz -! generation of magnetic energy at the upper Y-boundary +! advection of Bx across the Y-boundary ! - rterms(6) = rterms(6) + sum(sum(vv * bb, 1) * bb(2,:,:)) * dxz + rterms(3) = rterms(3) & + - sum(sign(bb(1,:,:) * vv(2,:,:), bb(1,:,:))) * dxz + +! shear of By along the Y-boundary +! + rterms(5) = rterms(5) & + + sum(sign(bb(2,:,:) * vv(1,:,:), bb(1,:,:))) * dxz if (resistive) then -! interpolate the components of current density +! interpolate Jz ! #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 +! diffusion of Bx at the Y-boundary ! - rterms(9) = rterms(9) & - - sum(jj(3,:,:) * bb(1,:,:) - jj(1,:,:) * bb(3,:,:)) * dxz + rterms(7) = rterms(7) & + - sum(sign(jj(3,:,:), bb(1,:,:))) * dxz end if ! resistivity end if @@ -886,31 +789,30 @@ module user_problem ! 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 +! advection of Bx across the Z-boundary ! - rterms(4) = rterms(4) + sum(sum(bb * bb, 1) * vv(3,:,:)) * dxy + rterms(4) = rterms(4) & + + sum(sign(bb(1,:,:) * vv(3,:,:), bb(1,:,:))) * dxy -! generation of magnetic energy at the lower Z-boundary +! shear of Bz along the Z-boundary ! - rterms(7) = rterms(7) - sum(sum(vv * bb, 1) * bb(3,:,:)) * dxy + rterms(6) = rterms(6) & + - sum(sign(bb(3,:,:) * vv(1,:,:), bb(1,:,:))) * dxy if (resistive) then -! interpolate the components of current density +! interpolate Jy ! - 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 +! diffusion of Bx at the Z-boundary ! - rterms(10) = rterms(10) & - + sum(jj(1,:,:) * bb(2,:,:) - jj(2,:,:) * bb(1,:,:)) * dxy + rterms(8) = rterms(8) & + + sum(sign(jj(2,:,:), bb(1,:,:))) * dxy end if ! resistivity end if @@ -918,104 +820,47 @@ module user_problem 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 +! advection of Bx across the Z-boundary ! - rterms(4) = rterms(4) - sum(sum(bb * bb, 1) * vv(3,:,:)) * dxy + rterms(4) = rterms(4) & + - sum(sign(bb(1,:,:) * vv(3,:,:), bb(1,:,:))) * dxy -! generation of magnetic energy at the upper Z-boundary +! shear of Bz along the Z-boundary ! - rterms(7) = rterms(7) + sum(sum(vv * bb, 1) * bb(3,:,:)) * dxy + rterms(6) = rterms(6) & + + sum(sign(bb(3,:,:) * vv(1,:,:), bb(1,:,:))) * 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 +! diffusion of Bx at the Z-boundary ! - rterms(10) = rterms(10) & - - sum(jj(1,:,:) * bb(2,:,:) - jj(2,:,:) * bb(1,:,:)) * dxy + rterms(8) = rterms(8) & + - sum(sign(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 +! calculate X-derivative of ψ ! - 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 */ + call derivative_1st(1, dh(1), pdata%q(ibp,:,:,:), cur(1,:,:,:)) -! conversion between magnetic energy and internal energies -! - if (resistive) then - rterms(12) = rterms(12) & + rterms(9) = rterms(9) & #if NDIMS == 3 - - sum(cur(1:3,nb:ne,nb:ne,nb:ne)**2) * dvol + - sum(sign(cur(1,nb:ne,nb:ne,nb:ne), & + pdata%q(ibx,nb:ne,nb:ne,nb:ne))) * 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 + - sum(sign(cur(1,nb:ne,nb:ne, : ), & + pdata%q(ibx,nb:ne,nb:ne, : ))) * dvol #endif /* NDIMS == 3 */ pdata => pdata%next @@ -1027,11 +872,11 @@ module user_problem 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) + rterms(2) = 5.00d-01 * rterms(2) + rterms(7) = eta * rterms(7) + rterms(8) = eta * rterms(8) - write(runit,"(14es25.15e3)") time, rterms(1:13) + write(runit,"(10es25.15e3)") time, rterms(1:9) call flush_and_sync(runit) 100 continue