USER_PROBLEM: Implement flux tube reconnection problem.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
cd0e6ea675
commit
4cb96ca2b3
@ -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
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user