USER_PROBLEM: Integrate the magnetic flux at x=0 only.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
711fc23f03
commit
8a66bc5769
@ -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
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user