gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

USER_PROBLEM: Implement setup for GEM magnetic reconnection challenge.

Browse Source 
See, e.g., Birn et al., Journal of Geophysical Research, 2001, 106, A3,
3715-3719, DOI:10.1029/1999JA900449.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2022-10-17 11:09:30 -03:00
parent a09bcfeccb
commit b6d86db3dc
1 changed files with 179 additions and 1 deletions

180
sources/user_problem.F90

@ -31,6 +31,22 @@ 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 = 0.00d+00
real(kind=8), save :: psi0 = 1.00d-01
real(kind=8), save :: dlta = 5.00d-01
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
logical, save :: resistive = .false.
private
public :: initialize_user_problem, finalize_user_problem
public :: user_statistics
@ -58,6 +74,12 @@ module user_problem
!
subroutine initialize_user_problem(problem, rcount, verbose, status)
use boundaries , only : 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 +87,83 @@ module user_problem
logical , intent(in) :: verbose
integer , intent(out) :: status
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("psi0", psi0)
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 = 5.0d-01 * (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)
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, '(*) phi0' , psi0)
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_y => user_boundary_y
status = 0
!-------------------------------------------------------------------------------
@ -116,14 +213,73 @@ 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, xlen, ylen
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), dimension(nn) :: x, y
real(kind=8), dimension(nn) :: snx, csx, sny, csy, tgy
!-------------------------------------------------------------------------------
!
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
! equilibrium
!
tgy(:) = bamp * tanh(y(:) / dlta)
do j = 1, nn
do i = 1, nn
pdata%q(ibx,i,j,:) = tgy(j)
pdata%q(ibz,i,j,:) = sqrt(bgui**2 &
+ zeta**2 * max(0.0d+00, bamp**2 - tgy(j)**2))
end do
end do
pdata%q(iby,:,:,:) = 0.0d+00
pdata%q(ibp,:,:,:) = 0.0d+00
if (ipr > 0) then
pdata%q(ipr,:,:,:) = ptot - 5.0d-01 * sum(pdata%q(ibx:ibz,:,:,:)**2,1)
pdata%q(idn,:,:,:) = (dens / pres) * pdata%q(ipr,:,:,:)
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
! perturbation
!
snx(:) = sin(pi2 * x(:) / xlen) * pi2 / xlen
csx(:) = cos(pi2 * x(:) / xlen)
sny(:) = sin(pi * y(:) / ylen) * pi / ylen
csy(:) = cos(pi * y(:) / ylen)
do j = 1, nn
do i = 1, nn
pdata%q(ibx,i,j,:) = pdata%q(ibx,i,j,:) - psi0 * csx(i) * sny(j)
pdata%q(iby,i,j,:) = pdata%q(iby,i,j,:) + psi0 * snx(i) * csy(j)
end do
end do
#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
!-------------------------------------------------------------------------------
!
@ -247,6 +403,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 +413,27 @@ module user_problem
real(kind=8), dimension(:) , intent(in) :: x, y, z
real(kind=8), dimension(:,:,:,:), intent(inout) :: qn
integer :: i, j
!-------------------------------------------------------------------------------
!
if (js == 1) then
do j = nbl, 1, -1
do i = il, iu
qn( : ,i,j,:) = qn( : ,i,nb,:)
qn(ivy,i,j,:) = 0.0d+00
qn(iby,i,j,:) = 0.0d+00
end do
end do
else
do j = neu, nn
do i = il, iu
qn( : ,i,j,:) = qn( : ,i,ne,:)
qn(ivy,i,j,:) = 0.0d+00
qn(iby,i,j,:) = 0.0d+00
end do
end do
end if
!-------------------------------------------------------------------------------
!