From 8a66bc5769c1e513ada7ba85d492e6d74b5ddd73 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 4 Nov 2022 14:07:35 -0300 Subject: [PATCH] USER_PROBLEM: Integrate the magnetic flux at x=0 only. Signed-off-by: Grzegorz Kowal --- sources/user_problem.F90 | 280 +++++++++++++++++++-------------------- 1 file changed, 135 insertions(+), 145 deletions(-) diff --git a/sources/user_problem.F90 b/sources/user_problem.F90 index 07cf65b..51bd607 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, "('reconnection.dat')") + write(fname, "('magnetic_flux.dat')") inquire(file = fname, exist = flag) case default - write(fname, "('reconnection_',i2.2,'.dat')") rcount + write(fname, "('magnetic_flux_',i2.2,'.dat')") rcount flag = .false. end select @@ -183,7 +183,7 @@ module user_problem ! write the integral file header ! - write(runit,'("#",a24,9a25)') 'time', '|Bx| int' , '|Bx| inf' , & + write(runit,'("#",a24,8a25)') 'time', '|Bx| int' , & '|Bx| y-adv', '|Bx| z-adv', & '|Bx| y-shr', '|Bx| z-shr', & '|Bx| y-dif', '|Bx| z-dif', & @@ -562,12 +562,16 @@ module user_problem 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 + use coordinates, only : periodic, ymin, ymax, zmin, zmax #else /* NDIMS == 3 */ - use coordinates, only : periodic, xmin, xmax, ymin, ymax + use coordinates, only : periodic, ymin, ymax #endif /* NDIMS == 3 */ - use coordinates, only : adx, ady, adz, advol + use coordinates, only : ax, adx, ady, adz +#if NDIMS == 3 use equations , only : ivx, ivy, ivz, ibx, iby, ibz, ibp +#else /* NDIMS == 3 */ + use equations , only : ivx, ivy, ibx, iby, ibz, ibp +#endif /* NDIMS == 3 */ use helpers , only : print_message, flush_and_sync #ifdef MPI use mpitools , only : reduce_sum @@ -582,20 +586,17 @@ module user_problem type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata - integer, save :: nt + integer, save :: nt, i !$ 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) :: dyz real(kind=8), dimension(3) :: dh real(kind=8), dimension(16) :: rterms + real(kind=8), dimension(nn) :: x real(kind=8), dimension(:,:,:,:), pointer, save :: cur - real(kind=8), dimension(:,:,:) , pointer, save :: vv, bb, jj -!$omp threadprivate(first, nt, cur, vv, bb, jj) + real(kind=8), dimension(:,:) , pointer, save :: vv, bb, jj +!$omp threadprivate(first, nt, i, cur, vv, bb, jj) character(len=*), parameter :: loc = 'USER_PROBLEM:user_time_statistics()' @@ -605,9 +606,9 @@ module user_problem !$ nt = omp_get_thread_num() if (first) then if (resistive) then - nu = 3 * nn**NDIMS + 9 * nc**(NDIMS - 1) + nu = 3 * nn**NDIMS + 9 * nc**(NDIMS - 2) else - nu = 3 * nn**NDIMS + 6 * nc**(NDIMS - 1) + nu = 3 * nn**NDIMS + 6 * nc**(NDIMS - 2) end if call resize_workspace(nu, status) @@ -624,26 +625,26 @@ module user_problem 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) + nu = nl - 1 + 3 * nc**(NDIMS - 2) #if NDIMS == 3 - vv(1:3,1:nc,1:nc) => work(nl:nu,nt) + vv(1:3,1:nc) => work(nl:nu,nt) #else /* NDIMS == 3 */ - vv(1:3,1:nc,1: 1) => work(nl:nu,nt) + vv(1:3,1: 1) => work(nl:nu,nt) #endif /* NDIMS == 3 */ nl = nu + 1 - nu = nl - 1 + 3 * nc**(NDIMS - 1) + nu = nl - 1 + 3 * nc**(NDIMS - 2) #if NDIMS == 3 - bb(1:3,1:nc,1:nc) => work(nl:nu,nt) + bb(1:3,1:nc) => work(nl:nu,nt) #else /* NDIMS == 3 */ - bb(1:3,1:nc,1: 1) => work(nl:nu,nt) + bb(1:3,1: 1) => work(nl:nu,nt) #endif /* NDIMS == 3 */ if (resistive) then nl = nu + 1 - nu = nl - 1 + 3 * nc**(NDIMS - 1) + nu = nl - 1 + 3 * nc**(NDIMS - 2) #if NDIMS == 3 - jj(1:3,1:nc,1:nc) => work(nl:nu,nt) + jj(1:3,1:nc) => work(nl:nu,nt) #else /* NDIMS == 3 */ - jj(1:3,1:nc,1: 1) => work(nl:nu,nt) + jj(1:3,1: 1) => work(nl:nu,nt) #endif /* NDIMS == 3 */ end if @@ -661,208 +662,198 @@ module user_problem do while(associated(pdata)) pmeta => pdata%meta - dh(1) = adx(pmeta%level) - dh(2) = ady(pmeta%level) - dh(3) = adz(pmeta%level) + if (pmeta%xmin <= 0.0d+00 .and. pmeta%xmax >= 0.0d+00) then - dvol = advol(pmeta%level) + dh(1) = adx(pmeta%level) + dh(2) = ady(pmeta%level) + dh(3) = adz(pmeta%level) - dyz = ady(pmeta%level) * adz(pmeta%level) - dxz = adx(pmeta%level) * adz(pmeta%level) + x(:) = pmeta%xmin + ax(pmeta%level,:) + + i = nb + do while(abs(x(i)) > dh(1)) + i = i + 1 + end do + + if (nb <= i .and. i <= ne) then + + dyz = ady(pmeta%level) * adz(pmeta%level) + +! the integral of |Bx| +! + rterms(1) = rterms(1) & #if NDIMS == 3 - dxy = adx(pmeta%level) * ady(pmeta%level) + + sum(abs(pdata%u(ibx,i,nb:ne,nb:ne))) * dyz +#else /* NDIMS == 3 */ + + sum(abs(pdata%u(ibx,i,nb:ne, : ))) * dyz #endif /* NDIMS == 3 */ ! calculate current density (J = ∇xB) ! - call curl(dh(:), pdata%q(ibx:ibz,:,:,:), cur(1:3,:,:,:)) + if (resistive) & + call curl(dh(:), pdata%q(ibx:ibz,:,:,:), cur(1:3,:,:,:)) -! the integral of magnetic energy -! - rterms(1) = rterms(1) & -#if NDIMS == 3 - + sum(abs(pdata%u(ibx,nb:ne,nb:ne,nb:ne))) * dvol -#else /* NDIMS == 3 */ - + sum(abs(pdata%u(ibx,nb:ne,nb:ne, : ))) * dvol -#endif /* NDIMS == 3 */ + if (.not. periodic(2)) then - if (.not. periodic(2)) then - - if (pmeta%ymin < ymin) 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) - 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) + vv(1,:) = 5.0d-01 * sum(pdata%q(ivx,i,nbm:nb,nb:ne), 1) + vv(2,:) = 5.0d-01 * sum(pdata%q(ivy,i,nbm:nb,nb:ne), 1) + bb(1,:) = 5.0d-01 * sum(pdata%q(ibx,i,nbm:nb,nb:ne), 1) + bb(2,:) = 5.0d-01 * sum(pdata%q(iby,i,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) - 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) + vv(1,:) = 5.0d-01 * sum(pdata%q(ivx,i,nbm:nb, : ), 1) + vv(2,:) = 5.0d-01 * sum(pdata%q(ivy,i,nbm:nb, : ), 1) + bb(1,:) = 5.0d-01 * sum(pdata%q(ibx,i,nbm:nb, : ), 1) + bb(2,:) = 5.0d-01 * sum(pdata%q(iby,i,nbm:nb, : ), 1) #endif /* NDIMS == 3 */ -! mean Bx at the Y-boundary -! - rterms(2) = rterms(2) + sum(abs(bb(1,:,:))) * dxz - ! advection of Bx across the Y-boundary ! - rterms(3) = rterms(3) & - + sum(sign(bb(1,:,:) * vv(2,:,:), bb(1,:,:))) * dxz + rterms(2) = rterms(2) & + + sum(sign(bb(1,:) * vv(2,:), bb(1,:))) * dh(3) ! shear of By along the Y-boundary ! - rterms(5) = rterms(5) & - - sum(sign(bb(2,:,:) * vv(1,:,:), bb(1,:,:))) * dxz + rterms(4) = rterms(4) & + - sum(sign(bb(2,:) * vv(1,:), bb(1,:))) * dh(3) - if (resistive) then + if (resistive) then -! interpolate Jz -! #if NDIMS == 3 - jj(3,:,:) = 5.0d-01 * sum(cur(3,nb:ne,nbm:nb,nb:ne), 1) + jj(3,:) = 5.0d-01 * sum(cur(3,i,nbm:nb,nb:ne), 1) #else /* NDIMS == 3 */ - jj(3,:,:) = 5.0d-01 * sum(cur(3,nb:ne,nbm:nb, : ), 1) + jj(3,:) = 5.0d-01 * sum(cur(3,i,nbm:nb, : ), 1) #endif /* NDIMS == 3 */ ! diffusion of Bx at the Y-boundary ! - rterms(7) = rterms(7) & - + sum(sign(jj(3,:,:), bb(1,:,:))) * dxz + rterms(6) = rterms(6) & + + sum(sign(jj(3,:), bb(1,:))) * dh(3) - end if ! resistivity - end if + end if ! resistivity + end if - if (ymax < pmeta%ymax) then + 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) - 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) + vv(1,:) = 5.0d-01 * sum(pdata%q(ivx,i,ne:nep,nb:ne), 1) + vv(2,:) = 5.0d-01 * sum(pdata%q(ivy,i,ne:nep,nb:ne), 1) + bb(1,:) = 5.0d-01 * sum(pdata%q(ibx,i,ne:nep,nb:ne), 1) + bb(2,:) = 5.0d-01 * sum(pdata%q(iby,i,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) - 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) + vv(1,:) = 5.0d-01 * sum(pdata%q(ivx,i,ne:nep, : ), 1) + vv(2,:) = 5.0d-01 * sum(pdata%q(ivy,i,ne:nep, : ), 1) + bb(1,:) = 5.0d-01 * sum(pdata%q(ibx,i,ne:nep, : ), 1) + bb(2,:) = 5.0d-01 * sum(pdata%q(iby,i,ne:nep, : ), 1) #endif /* NDIMS == 3 */ -! mean Bx at the boundary -! - rterms(2) = rterms(2) + sum(abs(bb(1,:,:))) * dxz - ! advection of Bx across the Y-boundary ! - rterms(3) = rterms(3) & - - sum(sign(bb(1,:,:) * vv(2,:,:), bb(1,:,:))) * dxz + rterms(2) = rterms(2) & + - sum(sign(bb(1,:) * vv(2,:), bb(1,:))) * dh(3) ! shear of By along the Y-boundary ! - rterms(5) = rterms(5) & - + sum(sign(bb(2,:,:) * vv(1,:,:), bb(1,:,:))) * dxz + rterms(4) = rterms(4) & + + sum(sign(bb(2,:) * vv(1,:), bb(1,:))) * dh(3) - if (resistive) then + if (resistive) then -! interpolate Jz -! #if NDIMS == 3 - jj(3,:,:) = 5.0d-01 * sum(cur(3,nb:ne,ne:nep,nb:ne), 1) + jj(3,:) = 5.0d-01 * sum(cur(3,i,ne:nep,nb:ne), 1) #else /* NDIMS == 3 */ - jj(3,:,:) = 5.0d-01 * sum(cur(3,nb:ne,ne:nep, : ), 1) + jj(3,:) = 5.0d-01 * sum(cur(3,i,ne:nep, : ), 1) #endif /* NDIMS == 3 */ ! diffusion of Bx at the Y-boundary ! - rterms(7) = rterms(7) & - - sum(sign(jj(3,:,:), bb(1,:,:))) * dxz + rterms(6) = rterms(6) & + - sum(sign(jj(3,:), bb(1,:))) * dh(3) - end if ! resistivity - end if - end if ! non-periodic along Y + end if ! resistivity + end if + end if ! non-periodic along Y #if NDIMS == 3 - if (.not. periodic(3)) then + if (.not. periodic(3)) then - if (pmeta%zmin < zmin) 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(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(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nb:ne,nb:ne,nbm:nb), 1) + vv(1,:) = 5.0d-01 * sum(pdata%q(ivx,i,nb:ne,nbm:nb), 1) + vv(3,:) = 5.0d-01 * sum(pdata%q(ivz,i,nb:ne,nbm:nb), 1) + bb(1,:) = 5.0d-01 * sum(pdata%q(ibx,i,nb:ne,nbm:nb), 1) + bb(3,:) = 5.0d-01 * sum(pdata%q(ibz,i,nb:ne,nbm:nb), 1) ! advection of Bx across the Z-boundary ! - rterms(4) = rterms(4) & - + sum(sign(bb(1,:,:) * vv(3,:,:), bb(1,:,:))) * dxy + rterms(3) = rterms(3) & + + sum(sign(bb(1,:) * vv(3,:), bb(1,:))) * dh(2) ! shear of Bz along the Z-boundary ! - rterms(6) = rterms(6) & - - sum(sign(bb(3,:,:) * vv(1,:,:), bb(1,:,:))) * dxy + rterms(5) = rterms(5) & + - sum(sign(bb(3,:) * vv(1,:), bb(1,:))) * dh(2) - if (resistive) then - -! interpolate Jy -! - jj(2,:,:) = 5.0d-01 * sum(cur(2,nb:ne,nb:ne,nbm:nb), 1) + if (resistive) then + jj(2,:) = 5.0d-01 * sum(cur(2,i,nb:ne,nbm:nb), 1) ! diffusion of Bx at the Z-boundary ! - rterms(8) = rterms(8) & - + sum(sign(jj(2,:,:), bb(1,:,:))) * dxy + rterms(7) = rterms(7) & + + sum(sign(jj(2,:), bb(1,:))) * dh(2) - end if ! resistivity - end if + end if ! resistivity + end if - if (pmeta%zmax > zmax) then + if (pmeta%zmax > zmax) then - vv(1,:,:) = 5.0d-01 * sum(pdata%q(ivx,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(3,:,:) = 5.0d-01 * sum(pdata%q(ibz,nb:ne,nb:ne,ne:nep), 1) + vv(1,:) = 5.0d-01 * sum(pdata%q(ivx,i,nb:ne,ne:nep), 1) + vv(3,:) = 5.0d-01 * sum(pdata%q(ivz,i,nb:ne,ne:nep), 1) + bb(1,:) = 5.0d-01 * sum(pdata%q(ibx,i,nb:ne,ne:nep), 1) + bb(3,:) = 5.0d-01 * sum(pdata%q(ibz,i,nb:ne,ne:nep), 1) ! advection of Bx across the Z-boundary ! - rterms(4) = rterms(4) & - - sum(sign(bb(1,:,:) * vv(3,:,:), bb(1,:,:))) * dxy + rterms(3) = rterms(3) & + - sum(sign(bb(1,:) * vv(3,:), bb(1,:))) * dh(2) ! shear of Bz along the Z-boundary ! - rterms(6) = rterms(6) & - + sum(sign(bb(3,:,:) * vv(1,:,:), bb(1,:,:))) * dxy + rterms(5) = rterms(5) & + + sum(sign(bb(3,:) * vv(1,:), bb(1,:))) * dh(2) - if (resistive) then + if (resistive) then -! interpolate the components of current density -! - jj(2,:,:) = 5.0d-01 * sum(cur(2,nb:ne,nb:ne,ne:nep), 1) + jj(2,:) = 5.0d-01 * sum(cur(2,1,nb:ne,ne:nep), 1) ! diffusion of Bx at the Z-boundary ! - rterms(8) = rterms(8) & - - sum(sign(jj(2,:,:), bb(1,:,:))) * dxy + rterms(7) = rterms(7) & + - sum(sign(jj(2,:), bb(1,:))) * dh(2) - end if ! resistivity - end if - end if ! non-periodic along Z + end if ! resistivity + end if + end if ! non-periodic along Z #endif /* NDIMS == 3 */ ! calculate X-derivative of ψ ! - call derivative_1st(1, dh(1), pdata%q(ibp,:,:,:), cur(1,:,:,:)) + call derivative_1st(1, dh(1), pdata%q(ibp,:,:,:), cur(1,:,:,:)) - rterms(9) = rterms(9) & + rterms(8) = rterms(8) & #if NDIMS == 3 - - sum(sign(cur(1,nb:ne,nb:ne,nb:ne), & - pdata%q(ibx,nb:ne,nb:ne,nb:ne))) * dvol + - sum(sign(cur(1,i,nb:ne,nb:ne), & + pdata%q(ibx,i,nb:ne,nb:ne))) * dyz #else /* NDIMS == 3 */ - - sum(sign(cur(1,nb:ne,nb:ne, : ), & - pdata%q(ibx,nb:ne,nb:ne, : ))) * dvol + - sum(sign(cur(1,i,nb:ne, : ), & + pdata%q(ibx,i,nb:ne, : ))) * dyz #endif /* NDIMS == 3 */ + end if + end if + pdata => pdata%next end do @@ -872,11 +863,10 @@ module user_problem call reduce_sum(rterms(:)) #endif /* MPI */ - rterms(2) = 5.00d-01 * rterms(2) - rterms(7) = eta * rterms(7) - rterms(8) = eta * rterms(8) + rterms(7) = eta * rterms(6) + rterms(8) = eta * rterms(7) - write(runit,"(10es25.15e3)") time, rterms(1:9) + write(runit,"(9es25.15e3)") time, 5.0d-01 * rterms(1:8) call flush_and_sync(runit) 100 continue