gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
amun-code/sources/user_problem.F90
Grzegorz Kowal 80f1fdde06 USER_PROBLEM: Implement alternative setup for flux tubes. 
This setup assume two flux tubes formed by eliptic magnetic field lines.
The two flux tubes enter in contact in the midplane. Moreover, the
magnetic field lines change the polarization across the contact plane.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
2023-07-13 14:56:27 -03:00

935 lines
29 KiB
Fortran
Raw Blame History

!!******************************************************************************
!!
!! This file is part of the AMUN source code, a program to perform
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
!! adaptive mesh.
!!
!! Copyright (C) 2017-2023 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!! This program is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!*****************************************************************************
!!
!! module: USER_PROBLEM
!!
!! This module provides subroutines to setup custom problem.
!!
!!*****************************************************************************
!
module user_problem
implicit none
integer :: profile = 1
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 = 0.00d+00
real(kind=8), save :: dlta = 1.00d-02
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
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
contains
!
!===============================================================================
!
! subroutine INITIALIZE_USER_PROBLEM:
! ----------------------------------
!
! Subroutine initializes user problem. It could read problem parameters which
! are used in all subroutines defining this specific problem.
!
! Arguments:
!
! problem - the problem name
! rcount - the run count for restarted jobs
! verbose - a logical flag turning the information printing;
! status - an integer flag for error return value;
!
!===============================================================================
!
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
integer , intent(in) :: rcount
logical , intent(in) :: verbose
integer , intent(out) :: status
character(len=64) :: tmp, append = "off", fname
logical :: flag
character(len=*), parameter :: &
loc = 'USER_PROBLEM::initialize_user_problem()'
!-------------------------------------------------------------------------------
!
call get_parameter("profile", tmp)
select case(tmp)
case('bhattacharjee', 'sincos', 'trigonometric')
profile = 2
case default
profile = 1
end select
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("resistivity", eta)
if (eta < 0.0d+00) then
if (verbose) &
call print_message(loc, "Resistivity cannot be negative!")
status = 1
return
else
resistive = .true.
dlta = sqrt(eta)
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
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 = bamp / sqrt(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, "('magnetic-flux.dat')")
inquire(file = fname, exist = flag)
case default
write(fname, "('magnetic-flux_',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,8a25)') 'time', '|Bx| int' , &
'|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)")
select case(profile)
case(2)
call print_parameter(verbose, '(*) flux tubes profile', &
'Bhattacharjee et al. (2009)')
case default
call print_parameter(verbose, '(*) flux tubes profile', &
'Kowal et al. (2023)')
end select
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
!-------------------------------------------------------------------------------
!
end subroutine initialize_user_problem
!
!===============================================================================
!
! subroutine FINALIZE_USER_PROBLEM:
! --------------------------------
!
! Subroutine releases memory used by the module.
!
! Arguments:
!
! status - an integer flag for error return value;
!
!===============================================================================
!
subroutine finalize_user_problem(status)
implicit none
integer, intent(out) :: status
!-------------------------------------------------------------------------------
!
status = 0
!-------------------------------------------------------------------------------
!
end subroutine finalize_user_problem
!
!===============================================================================
!
! subroutine SETUP_USER_PROBLEM:
! -----------------------------
!
! Subroutine sets the initial conditions for the user specific problem.
!
! Arguments:
!
! pdata - pointer to the datablock structure of the currently initialized
! block;
!
!===============================================================================
!
subroutine setup_user_problem(pdata)
use blocks , only : block_data
use constants , only : pi, pi2
use coordinates, only : nn => bcells, xmin, xmax
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) :: r1, r2, bx, by, bz
real(kind=8), dimension(nn) :: x, y
real(kind=8), dimension(nn) :: snx, csx, sny, csy, thy, ch2
real(kind=8), dimension(nn,nn) :: az
!-------------------------------------------------------------------------------
!
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
select case(profile)
case(2)
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
if (abs(x(i)) <= 5.0d-01) then
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)
else
bx = 0.0d+00
by = 0.0d+00
end if
bz = sqrt(max(0.0d+00, bgui**2 + zeta * (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
case default
do j = 1, nn
do i = 1, nn
r1 = 5.0d-01 - sqrt(x(i)**2 + ((2.0d+00 * y(j) - 5.0d-01))**2)
r2 = 5.0d-01 - sqrt(x(i)**2 + ((2.0d+00 * y(j) + 5.0d-01))**2)
az(i,j) = 2.5d-01 * (r1 * (tanh(r1 / dlta) + 1.0d+00) &
+ r2 * (tanh(r2 / dlta) + 1.0d+00))
end do
end do
do j = 2, nn - 1
do i = 2, nn - 1
bx = (az(i,j+1) - az(i,j-1)) / (y(j+1) - y(j-1))
by = (az(i-1,j) - az(i+1,j)) / (x(i+1) - x(i-1))
bz = sqrt(max(0.0d+00, bgui**2 + zeta * (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
end select
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
!-------------------------------------------------------------------------------
!
end subroutine setup_user_problem
!
!===============================================================================
!
! subroutine UPDATE_USER_SOURCES:
! ------------------------------
!
! Subroutine adds the user defined source terms.
!
! Arguments:
!
! t, dt - the time and time increment;
! pdata - the pointer to a data block;
!
!===============================================================================
!
subroutine update_user_sources(t, dt, pdata)
use blocks, only : block_data
implicit none
real(kind=8) , intent(in) :: t, dt
type(block_data), pointer, intent(inout) :: pdata
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
!
end subroutine update_user_sources
!
!===============================================================================
!
! subroutine UPDATE_USER_SHAPES:
! -----------------------------
!
! Subroutine defines the regions updated by user.
!
! Arguments:
!
! pdata - pointer to the data block structure of the currently initialized
! block;
! time - time at the moment of update;
! dt - time step since the last update;
!
!===============================================================================
!
subroutine update_user_shapes(pdata, time, dt)
use blocks, only : block_data
implicit none
type(block_data), pointer, intent(inout) :: pdata
real(kind=8) , intent(in) :: time, dt
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
!
end subroutine update_user_shapes
!
!===============================================================================
!
! subroutine USER_BOUNDARY_X:
! --------------------------
!
! Subroutine updates ghost zones within the specific region along
! the X direction.
!
! Arguments:
!
! is - the block side along the X direction for the ghost zone update;
! jl, ju - the cell index limits for the Y direction;
! kl, ku - the cell index limits for the Z direction;
! t, dt - time and time increment;
! x, y, z - the block coordinates;
! qn - the array of variables to update;
!
!===============================================================================
!
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
real(kind=8) , intent(in) :: t, dt
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
!-------------------------------------------------------------------------------
!
end subroutine user_boundary_x
!
!===============================================================================
!
! subroutine USER_BOUNDARY_Y:
! --------------------------
!
! Subroutine updates ghost zones within the specific region along
! the Y direction.
!
! Arguments:
!
! js - the block side along the Y direction for the ghost zone update;
! il, iu - the cell index limits for the X direction;
! kl, ku - the cell index limits for the Z direction;
! t, dt - time and time increment;
! x, y, z - the block coordinates;
! qn - the array of variables to update;
!
!===============================================================================
!
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
real(kind=8) , intent(in) :: t, dt
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
!-------------------------------------------------------------------------------
!
end subroutine user_boundary_y
!
!===============================================================================
!
! subroutine USER_BOUNDARY_Z:
! --------------------------
!
! Subroutine updates ghost zones within the specific region along
! the Z direction.
!
! Arguments:
!
! ks - the block side along the Z direction for the ghost zone update;
! il, iu - the cell index limits for the X direction;
! jl, ju - the cell index limits for the Y direction;
! t, dt - time and time increment;
! x, y, z - the block coordinates;
! qn - the array of variables to update;
!
!===============================================================================
!
subroutine user_boundary_z(ks, il, iu, jl, ju, t, dt, x, y, z, qn)
implicit none
integer , intent(in) :: ks, il, iu,jl, ju
real(kind=8) , intent(in) :: t, dt
real(kind=8), dimension(:) , intent(in) :: x, y, z
real(kind=8), dimension(:,:,:,:), intent(inout) :: qn
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
!
end subroutine user_boundary_z
!
!===============================================================================
!
! subroutine USER_STATISTICS:
! -------------------------------
!
! Subroutine can be use to store user defined time statistics. The file to
! store these statistics should be properly created in subroutine
! initialize_user_problem() and closed in finalize_user_problem().
!
! Arguments:
!
! time - the current simulation time;
!
!===============================================================================
!
subroutine user_statistics(time)
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, ymin, ymax, zmin, zmax
#else /* NDIMS == 3 */
use coordinates, only : periodic, ymin, ymax
#endif /* NDIMS == 3 */
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
#endif /* MPI */
use operators , only : curl, derivative_1st
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, i
!$ integer :: omp_get_thread_num
integer :: nl, nu, status
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, i, 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 - 2)
else
nu = 3 * nn**NDIMS + 6 * nc**(NDIMS - 2)
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 - 2)
#if NDIMS == 3
vv(1:3,1:nc) => work(nl:nu,nt)
#else /* NDIMS == 3 */
vv(1:3,1: 1) => work(nl:nu,nt)
#endif /* NDIMS == 3 */
nl = nu + 1
nu = nl - 1 + 3 * nc**(NDIMS - 2)
#if NDIMS == 3
bb(1:3,1:nc) => work(nl:nu,nt)
#else /* NDIMS == 3 */
bb(1:3,1: 1) => work(nl:nu,nt)
#endif /* NDIMS == 3 */
if (resistive) then
nl = nu + 1
nu = nl - 1 + 3 * nc**(NDIMS - 2)
#if NDIMS == 3
jj(1:3,1:nc) => work(nl:nu,nt)
#else /* NDIMS == 3 */
jj(1:3,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
if (pmeta%xmin <= 0.0d+00 .and. pmeta%xmax >= 0.0d+00) then
dh(1) = adx(pmeta%level)
dh(2) = ady(pmeta%level)
dh(3) = adz(pmeta%level)
x(:) = pmeta%xmin + ax(pmeta%level,:)
i = nb
do while(abs(x(i)) > dh(1) .and. i <= ne)
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
+ 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)
!
if (resistive) &
call curl(dh(:), pdata%q(ibx:ibz,:,:,:), cur(1:3,:,:,:))
if (.not. periodic(2)) then
if (pmeta%ymin - dh(2) < ymin) then
#if NDIMS == 3
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,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 */
! advection of Bx across the Y-boundary
!
rterms(2) = rterms(2) &
+ sum(sign(bb(1,:) * vv(2,:), bb(1,:))) * dh(3)
! shear of By along the Y-boundary
!
rterms(4) = rterms(4) &
- sum(sign(bb(2,:) * vv(1,:), bb(1,:))) * dh(3)
if (resistive) then
#if NDIMS == 3
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,i,nbm:nb, : ), 1)
#endif /* NDIMS == 3 */
! diffusion of Bx at the Y-boundary
!
rterms(6) = rterms(6) &
+ sum(sign(jj(3,:), bb(1,:))) * dh(3)
end if ! resistivity
end if
if (pmeta%ymax + dh(2) > ymax) then
#if NDIMS == 3
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,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 */
! advection of Bx across the Y-boundary
!
rterms(2) = rterms(2) &
- sum(sign(bb(1,:) * vv(2,:), bb(1,:))) * dh(3)
! shear of By along the Y-boundary
!
rterms(4) = rterms(4) &
+ sum(sign(bb(2,:) * vv(1,:), bb(1,:))) * dh(3)
if (resistive) then
#if NDIMS == 3
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,i,ne:nep, : ), 1)
#endif /* NDIMS == 3 */
! diffusion of Bx at the Y-boundary
!
rterms(6) = rterms(6) &
- sum(sign(jj(3,:), bb(1,:))) * dh(3)
end if ! resistivity
end if
end if ! non-periodic along Y
#if NDIMS == 3
if (.not. periodic(3)) then
if (pmeta%zmin - dh(3) < zmin) then
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(3) = rterms(3) &
+ sum(sign(bb(1,:) * vv(3,:), bb(1,:))) * dh(2)
! shear of Bz along the Z-boundary
!
rterms(5) = rterms(5) &
- sum(sign(bb(3,:) * vv(1,:), bb(1,:))) * dh(2)
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(7) = rterms(7) &
+ sum(sign(jj(2,:), bb(1,:))) * dh(2)
end if ! resistivity
end if
if (pmeta%zmax + dh(3) > zmax) then
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(3) = rterms(3) &
- sum(sign(bb(1,:) * vv(3,:), bb(1,:))) * dh(2)
! shear of Bz along the Z-boundary
!
rterms(5) = rterms(5) &
+ sum(sign(bb(3,:) * vv(1,:), bb(1,:))) * dh(2)
if (resistive) then
jj(2,:) = 5.0d-01 * sum(cur(2,1,nb:ne,ne:nep), 1)
! diffusion of Bx at the Z-boundary
!
rterms(7) = rterms(7) &
- sum(sign(jj(2,:), bb(1,:))) * dh(2)
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,:,:,:))
rterms(8) = rterms(8) &
#if NDIMS == 3
- 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,i,nb:ne, : ), &
pdata%q(ibx,i,nb:ne, : ))) * dyz
#endif /* NDIMS == 3 */
end if
end if
pdata => pdata%next
end do
work_in_use(nt) = .false.
#ifdef MPI
call reduce_sum(rterms(:))
#endif /* MPI */
rterms(6) = eta * rterms(6)
rterms(7) = eta * rterms(7)
write(runit,"(9es25.15e3)") time, rterms(1:8)
call flush_and_sync(runit)
100 continue
!-------------------------------------------------------------------------------
!
end subroutine user_statistics
!===============================================================================
!
end module user_problem
Reference in New Issue View Git Blame Copy Permalink