PROBLEMS: Add resistive tearing instability test problem.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
e6e73dd4f1
commit
66192401aa
65
problems/tearing.in
Normal file
65
problems/tearing.in
Normal file
@ -0,0 +1,65 @@
|
||||
# problem name and parameters
|
||||
#
|
||||
problem = "tearing"
|
||||
resistivity = 1.00d-05
|
||||
viscosity = 1.00d-05
|
||||
kper = 1.00d+00
|
||||
alpha = 1.00d+01
|
||||
beta = 1.60d+00
|
||||
vper = 0.00d+00
|
||||
bper = 1.00d-02
|
||||
zeta = 1.00d+00
|
||||
|
||||
# physics
|
||||
#
|
||||
equation_system = "mhd"
|
||||
equation_of_state = "adi"
|
||||
|
||||
# methods
|
||||
#
|
||||
time_advance = "ssprk(4,3)"
|
||||
riemann_solver = "hll"
|
||||
reconstruction = "crmp7"
|
||||
limiter = "mc"
|
||||
fix_positivity = "off"
|
||||
|
||||
# mesh parameters
|
||||
#
|
||||
xblocks = 1
|
||||
yblocks = 1
|
||||
zblocks = 1
|
||||
xmin = -5.0d-01
|
||||
xmax = 5.0d-01
|
||||
ymin = -5.0d-01
|
||||
ymax = 5.0d-01
|
||||
zmin = -5.0d-01
|
||||
zmax = 5.0d-01
|
||||
|
||||
# refinement control
|
||||
#
|
||||
ncells = 16
|
||||
minlev = 1
|
||||
maxlev = 10
|
||||
refinement_variables = "dens pres curr"
|
||||
|
||||
# boundary conditions
|
||||
#
|
||||
xlbndry = "open"
|
||||
xubndry = "open"
|
||||
ylbndry = "periodic"
|
||||
yubndry = "periodic"
|
||||
zlbndry = "periodic"
|
||||
zubndry = "periodic"
|
||||
|
||||
# runtime control parameters
|
||||
#
|
||||
tmax = 1.0d+01
|
||||
cfl = 3.0d-01
|
||||
|
||||
# data output control
|
||||
#
|
||||
precise_snapshots = "on"
|
||||
snapshot_type = "p"
|
||||
snapshot_interval = 1.0d-01
|
||||
restart_number = -1
|
||||
integrals_interval = 10
|
||||
@ -165,6 +165,9 @@ module problems
|
||||
case("current-sheet", "current_sheet")
|
||||
setup_problem => setup_problem_current_sheet
|
||||
|
||||
case("tearing-mode", "tearing_mode", "tearing")
|
||||
setup_problem => setup_problem_tearing
|
||||
|
||||
case default
|
||||
setup_problem => setup_problem_user
|
||||
|
||||
@ -2061,6 +2064,213 @@ module problems
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine setup_problem_current_sheet
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine SETUP_PROBLEM_TEARING:
|
||||
! --------------------------------
|
||||
!
|
||||
! Subroutine sets the initial conditions for the resistive tearing instability
|
||||
! test problem.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
! pdata - pointer to the datablock structure of the currently initialized
|
||||
! block;
|
||||
!
|
||||
! References:
|
||||
!
|
||||
! [1] Lapenta, G.,
|
||||
! "Self-Feeding Turbulent Magnetic Reconnection on Macroscopic Scales",
|
||||
! Physical Review Letters, 2008, vol. 100, id. 235001,
|
||||
! http://dx.doi.org/10.1103/PhysRevLett.100.235001
|
||||
!
|
||||
! [2] Landi, S., Del Zanna, L., Papini, E., Pucci, F., and Velli, M.
|
||||
! "Resistive Magnetohydrodynamics Simulations of the Ideal Tearing Mode",
|
||||
! The Astrophysical Journal, 2015, vol. 806, id. 131,
|
||||
! http://dx.doi.org/10.1088/0004-637X/806/1/131
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine setup_problem_tearing(pdata)
|
||||
|
||||
! include external procedures and variables
|
||||
!
|
||||
use blocks , only : block_data
|
||||
use constants , only : pi2
|
||||
use coordinates, only : nn => bcells
|
||||
use coordinates, only : ax, ay, adx, ady
|
||||
use equations , only : prim2cons
|
||||
use equations , only : nv, magnetized, csnd2
|
||||
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
|
||||
use parameters , only : get_parameter
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! input arguments
|
||||
!
|
||||
type(block_data), pointer, intent(inout) :: pdata
|
||||
|
||||
! default parameter values
|
||||
!
|
||||
real(kind=8), save :: beta = 1.00d-01
|
||||
real(kind=8), save :: delta = 1.00d-08
|
||||
real(kind=8), save :: eta = 1.00d-08
|
||||
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 :: bper = 1.00d-02
|
||||
real(kind=8), save :: vper = 1.00d-03
|
||||
real(kind=8), save :: kper = 1.00d+00
|
||||
real(kind=8), save :: ss = 1.00d+00
|
||||
|
||||
! local saved parameters
|
||||
!
|
||||
logical, save :: first = .true.
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: i, j, k = 1
|
||||
real(kind=8) :: dx, dy, t
|
||||
|
||||
! local arrays
|
||||
!
|
||||
real(kind=8), dimension(nv,nn) :: q, u
|
||||
real(kind=8), dimension(nn) :: xc, xl, xu, yl, yu, xi
|
||||
real(kind=8), dimension(nn) :: vx, vy, bx, by, b0
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for the problem setup
|
||||
!
|
||||
call start_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! retrieve the problem parameters during the first subroutine call
|
||||
!
|
||||
if (first) then
|
||||
call get_parameter("resistivity", eta)
|
||||
call get_parameter("beta" , beta)
|
||||
call get_parameter("zeta" , zeta)
|
||||
call get_parameter("density" , dens)
|
||||
call get_parameter("bamp" , bamp)
|
||||
call get_parameter("bper" , bper)
|
||||
call get_parameter("vper" , vper)
|
||||
call get_parameter("kper" , kper)
|
||||
|
||||
ss = abs(bamp) / (sqrt(dens) * max(1.0d-16, eta))
|
||||
delta = ss**(-1.0d+00/3.0d+00)
|
||||
if (bper /= 0.0d+00) vper = 0.0d+00
|
||||
|
||||
first = .false.
|
||||
end if
|
||||
|
||||
! prepare the block coordinates
|
||||
!
|
||||
dx = adx(pdata%meta%level)
|
||||
dy = ady(pdata%meta%level)
|
||||
xc(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
|
||||
xl(:) = pdata%meta%xmin + ax(pdata%meta%level,:) - 0.5d+00 * dx
|
||||
xu(:) = xl(:) + dx
|
||||
yl(:) = pdata%meta%ymin + ay(pdata%meta%level,:) - 0.5d+00 * dy
|
||||
yu(:) = yl(:) + dy
|
||||
|
||||
! prepare the vector of velocity perturbation (varying along Y)
|
||||
!
|
||||
if (vper /= 0.0d+00) then
|
||||
xi(:) = xc(:) * sqrt(ss)
|
||||
vx(:) = vper * tanh(xi(:)) * exp(- xi(:)**2)
|
||||
vy(:) = vper * (2.0d+00 * xi(:) * tanh(xi(:)) - 1.0d+00 / cosh(xi(:))**2)&
|
||||
* exp(- xi(:)**2) * sqrt(ss) / (pi2 * kper)
|
||||
end if
|
||||
if (bper /= 0.0d+00) then
|
||||
bx(:) = bper * (sin(pi2 * xl(:)) - sin(pi2 * xu(:))) / (pi2 * dx)
|
||||
by(:) = bper * (cos(pi2 * xl(:)) - cos(pi2 * xu(:))) / (pi2 * dx * kper)
|
||||
end if
|
||||
|
||||
! set the fields, along X, which do not vary along the YZ planes
|
||||
!
|
||||
q(idn,:) = dens
|
||||
q(ivx,:) = 0.0d+00
|
||||
q(ivy,:) = 0.0d+00
|
||||
q(ivz,:) = 0.0d+00
|
||||
if (ipr > 0) q(ipr,:) = 0.5d+00 * beta
|
||||
|
||||
! if magnetic field is present, set its antiparallel configuration
|
||||
!
|
||||
if (magnetized) then
|
||||
|
||||
! set magnetic field components
|
||||
!
|
||||
q(ibx,:) = 0.0d+00
|
||||
do i = 1, nn
|
||||
t = delta * (log(cosh(xu(i) / delta)) &
|
||||
- log(cosh(xl(i) / delta))) / dx
|
||||
q(iby,i) = bamp * max(-1.0d+00, min(1.0d+00, t))
|
||||
q(ibz,i) = zeta * sqrt(bamp**2 - q(iby,i)**2)
|
||||
end do
|
||||
q(ibp,:) = 0.0d+00
|
||||
b0(:) = q(iby,:)
|
||||
|
||||
! correct the gas pressure due to the variance of magnetic pressure
|
||||
!
|
||||
if (ipr > 0) then
|
||||
q(ipr,:) = q(ipr,:) + 0.5d+00 * (bamp**2 - sum(q(ibx:ibz,:)**2,1))
|
||||
else
|
||||
q(idn,:) = 0.5d+00 * (beta + bamp**2 - sum(q(ibx:ibz,:)**2,1)) / csnd2
|
||||
end if
|
||||
|
||||
end if
|
||||
|
||||
! iterate over all positions in the YZ planes
|
||||
!
|
||||
#if NDIMS == 3
|
||||
do k = 1, nn
|
||||
#endif /* NDIMS == 3 */
|
||||
do j = 1, nn
|
||||
|
||||
! set the perturbation
|
||||
!
|
||||
if (vper /= 0.0d+00) then
|
||||
q(ivx,:) = vx(:) * (sin(pi2 * kper * yl(j)) &
|
||||
- sin(pi2 * kper * yu(j))) / (pi2 * kper * dy)
|
||||
q(ivy,:) = vy(:) * (cos(pi2 * kper * yu(j)) &
|
||||
- cos(pi2 * kper * yl(j))) / (pi2 * kper * dy)
|
||||
end if
|
||||
if (bper /= 0.0d+00) then
|
||||
q(ibx,:) = bx(:) * (cos(pi2 * kper * yu(j)) &
|
||||
- cos(pi2 * kper * yl(j))) / (pi2 * kper * dy)
|
||||
q(iby,:) = by(:) * (sin(pi2 * kper * yl(j)) &
|
||||
- sin(pi2 * kper * yu(j))) / (pi2 * kper * dy) &
|
||||
+ b0(:)
|
||||
end if
|
||||
|
||||
! convert the primitive variables to conservative ones
|
||||
!
|
||||
call prim2cons(q(:,:), u(:,:))
|
||||
|
||||
! copy the conserved and primitive variables to the current block
|
||||
!
|
||||
pdata%u(:,:,j,k) = u(:,:)
|
||||
pdata%q(:,:,j,k) = q(:,:)
|
||||
|
||||
end do ! j = 1, nn
|
||||
#if NDIMS == 3
|
||||
end do ! k = 1, nn
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for the problems setup
|
||||
!
|
||||
call stop_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine setup_problem_tearing
|
||||
|
||||
!===============================================================================
|
||||
!
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user