USER_PROBLEM: Calculate reconnection rate terms in user_time_statistics().
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
62f0ca2deb
commit
e4a7890959
@ -74,6 +74,13 @@ module user_problem
|
|||||||
real(kind=8), save :: xdec = 1.00d-01
|
real(kind=8), save :: xdec = 1.00d-01
|
||||||
real(kind=8), save :: ydec = 1.00d-01
|
real(kind=8), save :: ydec = 1.00d-01
|
||||||
|
|
||||||
|
real(kind=8), save :: ylo = -9.00d+99
|
||||||
|
real(kind=8), save :: yup = 9.00d+99
|
||||||
|
|
||||||
|
integer(kind=4), save :: runit = 33
|
||||||
|
|
||||||
|
logical, save :: resistive = .false.
|
||||||
|
|
||||||
! flag indicating if the gravitational source term is enabled
|
! flag indicating if the gravitational source term is enabled
|
||||||
!
|
!
|
||||||
logical, save :: user_gravity_enabled = .false.
|
logical, save :: user_gravity_enabled = .false.
|
||||||
@ -120,7 +127,7 @@ module user_problem
|
|||||||
use constants , only : pi
|
use constants , only : pi
|
||||||
#endif /* NDIMS == 3 */
|
#endif /* NDIMS == 3 */
|
||||||
use constants , only : pi2
|
use constants , only : pi2
|
||||||
use coordinates, only : ng => nghosts, ady, xlen, zlen
|
use coordinates, only : ng => nghosts, ady, xlen, zlen, ymax
|
||||||
use equations , only : magnetized, csnd, csnd2
|
use equations , only : magnetized, csnd, csnd2
|
||||||
use helpers , only : print_section, print_parameter
|
use helpers , only : print_section, print_parameter
|
||||||
use parameters , only : get_parameter
|
use parameters , only : get_parameter
|
||||||
@ -139,9 +146,10 @@ module user_problem
|
|||||||
|
|
||||||
! local variables
|
! local variables
|
||||||
!
|
!
|
||||||
character(len=64) :: perturbation = "noise"
|
character(len=64) :: perturbation = "noise", append = "off", fname
|
||||||
|
logical :: flag
|
||||||
integer :: n, kd
|
integer :: n, kd
|
||||||
real(kind=8) :: thh, fc, ka
|
real(kind=8) :: thh, fc, ka, ydis = 9.00d+99
|
||||||
#if NDIMS == 3
|
#if NDIMS == 3
|
||||||
real(kind=8) :: thv, tx, ty, tz, tt
|
real(kind=8) :: thv, tx, ty, tz, tt
|
||||||
#endif /* NDIMS == 3 */
|
#endif /* NDIMS == 3 */
|
||||||
@ -216,6 +224,8 @@ module user_problem
|
|||||||
write(*,*)
|
write(*,*)
|
||||||
end if
|
end if
|
||||||
status = 1
|
status = 1
|
||||||
|
else
|
||||||
|
resistive = .true.
|
||||||
end if
|
end if
|
||||||
call get_parameter("dens", dens)
|
call get_parameter("dens", dens)
|
||||||
if (dens <= 0.0d+00) then
|
if (dens <= 0.0d+00) then
|
||||||
@ -353,6 +363,64 @@ module user_problem
|
|||||||
end if
|
end if
|
||||||
end if ! status
|
end if ! status
|
||||||
|
|
||||||
|
! prepare file to store reconnection rate terms
|
||||||
|
!
|
||||||
|
if (status == 0) then
|
||||||
|
|
||||||
|
! determine if to append or create another file
|
||||||
|
!
|
||||||
|
call get_parameter("reconnection_append", append)
|
||||||
|
select case(trim(append))
|
||||||
|
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
||||||
|
write(fname, "('reconnection-new.dat')")
|
||||||
|
inquire(file = fname, exist = flag)
|
||||||
|
case default
|
||||||
|
write(fname, "('reconnection-new_',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,26a25)') &
|
||||||
|
'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' &
|
||||||
|
, 'Vin lower' , 'Vin upper' &
|
||||||
|
, 'Emag', 'Emag inf' &
|
||||||
|
, '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 get_parameter("ydistance", ydis)
|
||||||
|
ydis = min(ydis, ymax)
|
||||||
|
ylo = - ydis
|
||||||
|
yup = ydis
|
||||||
|
|
||||||
|
end if ! status
|
||||||
|
|
||||||
! proceed if no errors
|
! proceed if no errors
|
||||||
!
|
!
|
||||||
if (status == 0) then
|
if (status == 0) then
|
||||||
@ -370,7 +438,7 @@ module user_problem
|
|||||||
call print_parameter(verbose, '( ) ptot (total pressure)', ptot)
|
call print_parameter(verbose, '( ) ptot (total pressure)', ptot)
|
||||||
call print_parameter(verbose, '( ) csnd (sound speed)' , csnd)
|
call print_parameter(verbose, '( ) csnd (sound speed)' , csnd)
|
||||||
call print_parameter(verbose, '( ) Valf (Alfven speed)' , valf)
|
call print_parameter(verbose, '( ) Valf (Alfven speed)' , valf)
|
||||||
if (eta > 0.0d+00) then
|
if (resistive) then
|
||||||
call print_parameter(verbose, '( ) S (Lundquist number)', lund)
|
call print_parameter(verbose, '( ) S (Lundquist number)', lund)
|
||||||
end if
|
end if
|
||||||
call print_parameter(verbose, '(*) delta (thickness)', dlta)
|
call print_parameter(verbose, '(*) delta (thickness)', dlta)
|
||||||
@ -388,6 +456,7 @@ module user_problem
|
|||||||
call print_parameter(verbose, '(*) ycut' , ycut)
|
call print_parameter(verbose, '(*) ycut' , ycut)
|
||||||
call print_parameter(verbose, '(*) xdec' , xdec)
|
call print_parameter(verbose, '(*) xdec' , xdec)
|
||||||
call print_parameter(verbose, '(*) ydec' , ydec)
|
call print_parameter(verbose, '(*) ydec' , ydec)
|
||||||
|
call print_parameter(verbose, '(*) ydistance' , ydis)
|
||||||
end if ! status
|
end if ! status
|
||||||
|
|
||||||
#ifdef PROFILE
|
#ifdef PROFILE
|
||||||
@ -435,6 +504,10 @@ module user_problem
|
|||||||
!
|
!
|
||||||
status = 0
|
status = 0
|
||||||
|
|
||||||
|
! close the reconnection file
|
||||||
|
!
|
||||||
|
close(runit)
|
||||||
|
|
||||||
! deallocate wave vector components, random directions, and random phase
|
! deallocate wave vector components, random directions, and random phase
|
||||||
!
|
!
|
||||||
if (allocated(kx)) deallocate(kx, ky, kz, stat = status)
|
if (allocated(kx)) deallocate(kx, ky, kz, stat = status)
|
||||||
@ -1555,16 +1628,408 @@ module user_problem
|
|||||||
!
|
!
|
||||||
subroutine user_time_statistics(time)
|
subroutine user_time_statistics(time)
|
||||||
|
|
||||||
|
use blocks , only : block_meta, block_data, list_data
|
||||||
|
use coordinates, only : nn => bcells, nb, ne
|
||||||
|
use coordinates, only : adx, ady, adz, advol, ay, yarea
|
||||||
|
use equations , only : ivx, ivy, ivz, ibx, iby, ibz, ibp
|
||||||
|
use helpers , only : flush_and_sync
|
||||||
|
#ifdef MPI
|
||||||
|
use mpitools , only : reduce_sum
|
||||||
|
#endif /* MPI */
|
||||||
|
use operators , only : curl, gradient
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
real(kind=8), intent(in) :: time
|
real(kind=8), intent(in) :: time
|
||||||
|
|
||||||
|
type(block_meta), pointer :: pmeta
|
||||||
|
type(block_data), pointer :: pdata
|
||||||
|
|
||||||
|
integer :: ni, nl, nu, nm, np
|
||||||
|
real(kind=8) :: dvol, dvolh, dxz
|
||||||
|
real(kind=8), dimension(3) :: dh
|
||||||
|
real(kind=8), dimension(nn) :: yc
|
||||||
|
real(kind=8), dimension(32) :: rterms
|
||||||
|
|
||||||
|
#if NDIMS == 3
|
||||||
|
real(kind=8), dimension(3,nn,nn,nn) :: cur
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
real(kind=8), dimension(3,nn,nn, 1) :: cur
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
!-------------------------------------------------------------------------------
|
!-------------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
#ifdef PROFILE
|
#ifdef PROFILE
|
||||||
call start_timer(imt)
|
call start_timer(imt)
|
||||||
#endif /* PROFILE */
|
#endif /* PROFILE */
|
||||||
|
|
||||||
|
rterms(:) = 0.0d+00
|
||||||
|
|
||||||
|
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)
|
||||||
|
dvolh = 5.0d-01 * dvol
|
||||||
|
|
||||||
|
dxz = adx(pmeta%level) * adz(pmeta%level)
|
||||||
|
|
||||||
|
yc(:) = pmeta%ymin + ay(pmeta%level,:)
|
||||||
|
ni = -1
|
||||||
|
nl = nb
|
||||||
|
nu = ne
|
||||||
|
do while (yc(nl) < ylo .and. nl < ne)
|
||||||
|
nl = nl + 1
|
||||||
|
end do
|
||||||
|
do while (yc(nu) > yup .and. nu > nb)
|
||||||
|
nu = nu - 1
|
||||||
|
end do
|
||||||
|
|
||||||
|
! calculate current density (J = ∇xB)
|
||||||
|
!
|
||||||
|
call curl(dh(:), pdata%q(ibx:ibz,:,:,:), cur(1:3,:,:,:))
|
||||||
|
|
||||||
|
! the integral of |Bx|
|
||||||
|
!
|
||||||
|
rterms(1) = rterms(1) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum(abs(pdata%q(ibx,nb:ne,nl:nu,nb:ne))) * dvol
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum(abs(pdata%q(ibx,nb:ne,nl:nu, : ))) * dvol
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! the integral of magnetic energy
|
||||||
|
!
|
||||||
|
rterms(12) = rterms(12) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum(pdata%u(ibx:ibz,nb:ne,nl:nu,nb:ne)**2) * dvolh
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum(pdata%u(ibx:ibz,nb:ne,nl:nu, : )**2) * dvolh
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
if (pmeta%ymin <= ylo .and. ylo < pmeta%ymax) then
|
||||||
|
ni = nl
|
||||||
|
nm = ni - 1
|
||||||
|
|
||||||
|
! get the inflow speed
|
||||||
|
!
|
||||||
|
rterms(10) = rterms(10) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum(pdata%q(ivy,nb:ne,ni,nb:ne)) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum(pdata%q(ivy,nb:ne,ni, : )) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! mean Bx at boundary
|
||||||
|
!
|
||||||
|
rterms(2) = rterms(2) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum(abs(pdata%q(ibx,nb:ne,nm:ni,nb:ne))) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum(abs(pdata%q(ibx,nb:ne,nm:ni, : ))) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! advection of Bx along Y
|
||||||
|
!
|
||||||
|
rterms(3) = rterms(3) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum(abs(pdata%q(ibx,nb:ne,nm:ni,nb:ne)) &
|
||||||
|
* pdata%q(ivy,nb:ne,nm:ni,nb:ne)) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum(abs(pdata%q(ibx,nb:ne,nm:ni, : )) &
|
||||||
|
* pdata%q(ivy,nb:ne,nm:ni, : )) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! shear of By along X
|
||||||
|
!
|
||||||
|
rterms(5) = rterms(5) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
- sum(sign(pdata%q(iby,nb:ne,nm:ni,nb:ne) &
|
||||||
|
* pdata%q(ivx,nb:ne,nm:ni,nb:ne), &
|
||||||
|
pdata%q(ibx,nb:ne,nm:ni,nb:ne))) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
- sum(sign(pdata%q(iby,nb:ne,nm:ni, : ) &
|
||||||
|
* pdata%q(ivx,nb:ne,nm:ni, : ), &
|
||||||
|
pdata%q(ibx,nb:ne,nm:ni, : ))) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! mean magnetic energy
|
||||||
|
!
|
||||||
|
rterms(13) = rterms(13) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum(pdata%q(ibx:ibz,nb:ne,nm:ni,nb:ne)**2) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum(pdata%q(ibx:ibz,nb:ne,nm:ni, : )**2) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! advection of magnetic energy
|
||||||
|
!
|
||||||
|
rterms(15) = rterms(15) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum(sum(pdata%q(ibx:ibz,nb:ne,nm:ni,nb:ne)**2,1) &
|
||||||
|
* pdata%q(ivy ,nb:ne,nm:ni,nb:ne)) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum(sum(pdata%q(ibx:ibz,nb:ne,nm:ni, : )**2,1) &
|
||||||
|
* pdata%q(ivy ,nb:ne,nm:ni, : )) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! advection of magnetic energy
|
||||||
|
!
|
||||||
|
rterms(18) = rterms(18) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
- sum(sum(pdata%q(ivx:ivz,nb:ne,nm:ni,nb:ne) &
|
||||||
|
* pdata%q(ibx:ibz,nb:ne,nm:ni,nb:ne),1) &
|
||||||
|
* pdata%q(iby ,nb:ne,nm:ni,nb:ne)) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
- sum(sum(pdata%q(ivx:ivz,nb:ne,nm:ni, : ) &
|
||||||
|
* pdata%q(ibx:ibz,nb:ne,nm:ni, : ),1) &
|
||||||
|
* pdata%q(iby ,nb:ne,nm:ni, : )) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
if (resistive) then
|
||||||
|
|
||||||
|
! diffusion of Bx through
|
||||||
|
!
|
||||||
|
rterms(7) = rterms(7) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum(sign(cur(3,nb:ne,nm:ni,nb:ne), &
|
||||||
|
pdata%q(ibx,nb:ne,nm:ni,nb:ne))) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum(sign(cur(3,nb:ne,nm:ni, : ), &
|
||||||
|
pdata%q(ibx,nb:ne,nm:ni, : ))) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! diffusion of magnetic energy through the lower Y boundary
|
||||||
|
!
|
||||||
|
rterms(21) = rterms(21) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum(cur(3,nb:ne,nm:ni,nb:ne) &
|
||||||
|
* pdata%q(ibx,nb:ne,nm:ni,nb:ne) &
|
||||||
|
- cur(1,nb:ne,nm:ni,nb:ne) &
|
||||||
|
* pdata%q(ibz,nb:ne,nm:ni,nb:ne)) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum(cur(3,nb:ne,nm:ni, : ) &
|
||||||
|
* pdata%q(ibx,nb:ne,nm:ni, : ) &
|
||||||
|
- cur(1,nb:ne,nm:ni, : ) &
|
||||||
|
* pdata%q(ibz,nb:ne,nm:ni, : )) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
end if ! resistivity
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
if (pmeta%ymin < yup .and. yup <= pmeta%ymax) then
|
||||||
|
ni = nu
|
||||||
|
np = ni + 1
|
||||||
|
|
||||||
|
! get the inflow speed
|
||||||
|
!
|
||||||
|
rterms(11) = rterms(11) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
- sum(pdata%q(ivy,nb:ne,ni,nb:ne)) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
- sum(pdata%q(ivy,nb:ne,ni, : )) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! mean Bx at boundary
|
||||||
|
!
|
||||||
|
rterms(2) = rterms(2) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum(abs(pdata%q(ibx,nb:ne,ni:np,nb:ne))) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum(abs(pdata%q(ibx,nb:ne,ni:np, : ))) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! advection of Bx along Y
|
||||||
|
!
|
||||||
|
rterms(3) = rterms(3) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
- sum(abs(pdata%q(ibx,nb:ne,ni:np,nb:ne)) &
|
||||||
|
* pdata%q(ivy,nb:ne,ni:np,nb:ne)) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
- sum(abs(pdata%q(ibx,nb:ne,ni:np, : )) &
|
||||||
|
* pdata%q(ivy,nb:ne,ni:np, : )) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! shear of By along X
|
||||||
|
!
|
||||||
|
rterms(5) = rterms(5) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum(sign(pdata%q(iby,nb:ne,ni:np,nb:ne) &
|
||||||
|
* pdata%q(ivx,nb:ne,ni:np,nb:ne), &
|
||||||
|
pdata%q(ibx,nb:ne,ni:np,nb:ne))) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum(sign(pdata%q(iby,nb:ne,ni:np, : ) &
|
||||||
|
* pdata%q(ivx,nb:ne,ni:np, : ), &
|
||||||
|
pdata%q(ibx,nb:ne,ni:np, : ))) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! mean magnetic energy
|
||||||
|
!
|
||||||
|
rterms(13) = rterms(13) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum(pdata%q(ibx:ibz,nb:ne,ni:np,nb:ne)**2) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum(pdata%q(ibx:ibz,nb:ne,ni:np, : )**2) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
|
||||||
|
! advection of magnetic energy
|
||||||
|
!
|
||||||
|
rterms(15) = rterms(15) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
- sum(sum(pdata%q(ibx:ibz,nb:ne,ni:np,nb:ne)**2,1) &
|
||||||
|
* pdata%q(ivy ,nb:ne,ni:np,nb:ne)) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
- sum(sum(pdata%q(ibx:ibz,nb:ne,ni:np, : )**2,1) &
|
||||||
|
* pdata%q(ivy ,nb:ne,ni:np, : )) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
rterms(18) = rterms(18) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum(sum(pdata%q(ivx:ivz,nb:ne,ni:np,nb:ne) &
|
||||||
|
* pdata%q(ibx:ibz,nb:ne,ni:np,nb:ne),1) &
|
||||||
|
* pdata%q(iby ,nb:ne,ni:np,nb:ne)) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum(sum(pdata%q(ivx:ivz,nb:ne,ni:np, : ) &
|
||||||
|
* pdata%q(ibx:ibz,nb:ne,ni:np, : ),1) &
|
||||||
|
* pdata%q(iby ,nb:ne,ni:np, : )) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
if (resistive) then
|
||||||
|
|
||||||
|
! diffusion of Bx
|
||||||
|
!
|
||||||
|
rterms(7) = rterms(7) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
- sum(sign(cur(3,nb:ne,ni:np,nb:ne), &
|
||||||
|
pdata%q(ibx,nb:ne,ni:np,nb:ne))) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
- sum(sign(cur(3,nb:ne,ni:np, : ), &
|
||||||
|
pdata%q(ibx,nb:ne,ni:np, : ))) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! diffusion of magnetic energy
|
||||||
|
!
|
||||||
|
rterms(21) = rterms(21) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
- sum(cur(3,nb:ne,ni:np,nb:ne) &
|
||||||
|
* pdata%q(ibx,nb:ne,ni:np,nb:ne) &
|
||||||
|
- cur(1,nb:ne,ni:np,nb:ne) &
|
||||||
|
* pdata%q(ibz,nb:ne,ni:np,nb:ne)) * dxz
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
- sum(cur(3,nb:ne,ni:np, : ) &
|
||||||
|
* pdata%q(ibx,nb:ne,ni:np, : ) &
|
||||||
|
- cur(1,nb:ne,ni:np, : ) &
|
||||||
|
* pdata%q(ibz,nb:ne,ni:np, : )) * dxz
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
end if ! resistivity
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! get the conversion between the magnetic energy and kinetic and
|
||||||
|
! internal energies
|
||||||
|
!
|
||||||
|
rterms(23) = rterms(23) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum((pdata%q(ivy,nb:ne,nl:nu,nb:ne) &
|
||||||
|
* pdata%q(ibz,nb:ne,nl:nu,nb:ne) &
|
||||||
|
- pdata%q(ivz,nb:ne,nl:nu,nb:ne) &
|
||||||
|
* pdata%q(iby,nb:ne,nl:nu,nb:ne)) &
|
||||||
|
* cur(1,nb:ne,nl:nu,nb:ne)) * dvol
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum((pdata%q(ivy,nb:ne,nl:nu, : ) &
|
||||||
|
* pdata%q(ibz,nb:ne,nl:nu, : ) &
|
||||||
|
- pdata%q(ivz,nb:ne,nl:nu, : ) &
|
||||||
|
* pdata%q(iby,nb:ne,nl:nu, : )) &
|
||||||
|
* cur(1,nb:ne,nl:nu, : )) * dvol
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
rterms(23) = rterms(23) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum((pdata%q(ivz,nb:ne,nl:nu,nb:ne) &
|
||||||
|
* pdata%q(ibx,nb:ne,nl:nu,nb:ne) &
|
||||||
|
- pdata%q(ivx,nb:ne,nl:nu,nb:ne) &
|
||||||
|
* pdata%q(ibz,nb:ne,nl:nu,nb:ne)) &
|
||||||
|
* cur(2,nb:ne,nl:nu,nb:ne)) * dvol
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum((pdata%q(ivz,nb:ne,nl:nu, : ) &
|
||||||
|
* pdata%q(ibx,nb:ne,nl:nu, : ) &
|
||||||
|
- pdata%q(ivx,nb:ne,nl:nu, : ) &
|
||||||
|
* pdata%q(ibz,nb:ne,nl:nu, : )) &
|
||||||
|
* cur(2,nb:ne,nl:nu, : )) * dvol
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
rterms(23) = rterms(23) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
+ sum((pdata%q(ivx,nb:ne,nl:nu,nb:ne) &
|
||||||
|
* pdata%q(iby,nb:ne,nl:nu,nb:ne) &
|
||||||
|
- pdata%q(ivy,nb:ne,nl:nu,nb:ne) &
|
||||||
|
* pdata%q(ibx,nb:ne,nl:nu,nb:ne)) &
|
||||||
|
* cur(3,nb:ne,nl:nu,nb:ne)) * dvol
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
+ sum((pdata%q(ivx,nb:ne,nl:nu, : ) &
|
||||||
|
* pdata%q(iby,nb:ne,nl:nu, : ) &
|
||||||
|
- pdata%q(ivy,nb:ne,nl:nu, : ) &
|
||||||
|
* pdata%q(ibx,nb:ne,nl:nu, : )) &
|
||||||
|
* cur(3,nb:ne,nl:nu, : )) * dvol
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
if (resistive) then
|
||||||
|
rterms(24) = rterms(24) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
- sum(cur(1:3,nb:ne,nl:nu,nb:ne)**2) * dvol
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
- sum(cur(1:3,nb:ne,nl:nu, : )**2) * dvol
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
end if
|
||||||
|
|
||||||
|
! calculate gradient (∇ψ)
|
||||||
|
!
|
||||||
|
call gradient(dh(:), pdata%q(ibp,:,:,:), cur(1:3,:,:,:))
|
||||||
|
|
||||||
|
rterms(25) = rterms(25) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
- sum(pdata%q(ibx:ibz,nb:ne,nl:nu,nb:ne) &
|
||||||
|
* cur(1:3,nb:ne,nl:nu,nb:ne)) * dvol
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
- sum(pdata%q(ibx:ibz,nb:ne,nl:nu, : ) &
|
||||||
|
* cur(1:3,nb:ne,nl:nu, : )) * dvol
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! contribution to |Bx| from ∇ψ
|
||||||
|
!
|
||||||
|
rterms(9) = rterms(9) &
|
||||||
|
#if NDIMS == 3
|
||||||
|
- sum(sign(cur(1,nb:ne,nl:nu,nb:ne), &
|
||||||
|
pdata%q(ibx,nb:ne,nl:nu,nb:ne))) * dvol
|
||||||
|
#else /* NDIMS == 3 */
|
||||||
|
- sum(sign(cur(1,nb:ne,nl:nu, : ), &
|
||||||
|
pdata%q(ibx,nb:ne,nl:nu, : ))) * dvol
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
pdata => pdata%next
|
||||||
|
end do
|
||||||
|
|
||||||
|
#ifdef MPI
|
||||||
|
call reduce_sum(rterms(:))
|
||||||
|
#endif /* MPI */
|
||||||
|
|
||||||
|
rterms(2) = 2.50d-01 * rterms(2) / yarea
|
||||||
|
rterms(13) = 1.25d-01 * rterms(13) / yarea
|
||||||
|
rterms(3:6) = 5.00d-01 * rterms(3:6)
|
||||||
|
rterms(14:19) = 5.00d-01 * rterms(14:19)
|
||||||
|
rterms(7:8) = 5.00d-01 * eta * rterms(7:8)
|
||||||
|
rterms(20:22) = 5.00d-01 * eta * rterms(20:22)
|
||||||
|
rterms(24) = eta * rterms(24)
|
||||||
|
|
||||||
|
write(runit,"(26es25.15e3)") time, rterms(1:25)
|
||||||
|
call flush_and_sync(runit)
|
||||||
|
|
||||||
#ifdef PROFILE
|
#ifdef PROFILE
|
||||||
call stop_timer(imt)
|
call stop_timer(imt)
|
||||||
#endif /* PROFILE */
|
#endif /* PROFILE */
|
||||||
|
|||||||
Loading…
x
Reference in New Issue
Block a user