USER_PROBLEM: Rewrite initial perturbation to inject in all modes.
Add a flag to select only the perturbation modes which are stable for tearing mode instability. Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
d4c692d09a
commit
f99a1fb0bf
@ -53,9 +53,10 @@ module user_problem
|
||||
|
||||
integer , save :: pert = 0
|
||||
integer , save :: nper = 16
|
||||
real(kind=8), save :: kper = 1.00d+00
|
||||
real(kind=8), save :: kdel = 1.00d+00
|
||||
real(kind=8), save :: bper = 0.00d+00
|
||||
real(kind=8), save :: vper = 0.00d+00
|
||||
real(kind=8), save :: kper = 1.00d+00
|
||||
real(kind=8), save :: kvec = 1.00d+00
|
||||
real(kind=8), save :: xcut = 1.00d+99
|
||||
real(kind=8), save :: ycut = 1.00d+99
|
||||
@ -106,10 +107,7 @@ module user_problem
|
||||
subroutine initialize_user_problem(problem, rcount, verbose, status)
|
||||
|
||||
use boundaries , only : custom_boundary_x, custom_boundary_y
|
||||
#if NDIMS == 3
|
||||
use constants , only : pi
|
||||
#endif /* NDIMS == 3 */
|
||||
use constants , only : pi2
|
||||
use constants , only : pi, pi2
|
||||
use coordinates, only : xlen, zlen, ymax
|
||||
use equations , only : magnetized, csnd, csnd2, qpbnd
|
||||
use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp, ipr
|
||||
@ -126,16 +124,14 @@ module user_problem
|
||||
logical , intent(in) :: verbose
|
||||
integer , intent(out) :: status
|
||||
|
||||
logical :: perturb_tearing = .true.
|
||||
character(len=32) :: stmp
|
||||
character(len=64) :: perturbation = "noise", append = "off", fname
|
||||
character(len=64) :: enable_absorption = "off"
|
||||
logical :: flag
|
||||
integer :: n, kd
|
||||
real(kind=8) :: thh, fc, ydis = 9.00d+99
|
||||
#if NDIMS == 3
|
||||
real(kind=8) :: thv, tx, ty, tz, tt
|
||||
#else /* NDIMS == 3 */
|
||||
real(kind=8) :: ka
|
||||
#endif /* NDIMS == 3 */
|
||||
integer :: n, kd, ki, kj, kk, kv, ks
|
||||
real(kind=8) :: kl2, ku2, kv2, ka, fc, ydis = 9.00d+99
|
||||
real(kind=8) :: thv, thh, tx, ty, tz
|
||||
|
||||
character(len=*), parameter :: &
|
||||
loc = 'USER_PROBLEM::initialize_user_problem()'
|
||||
@ -242,6 +238,8 @@ module user_problem
|
||||
status = 1
|
||||
return
|
||||
end if
|
||||
ycut = 2 * dlta
|
||||
ydec = dlta
|
||||
|
||||
! get the perturbation parameters
|
||||
!
|
||||
@ -250,6 +248,7 @@ module user_problem
|
||||
call get_parameter("bper" , bper)
|
||||
call get_parameter("vper" , vper)
|
||||
call get_parameter("kper" , kper)
|
||||
call get_parameter("kdel" , kdel)
|
||||
call get_parameter("xcut" , xcut)
|
||||
call get_parameter("ycut" , ycut)
|
||||
call get_parameter("xdec" , xdec)
|
||||
@ -273,6 +272,17 @@ module user_problem
|
||||
pert = 1
|
||||
end select
|
||||
|
||||
! enable or disable tearing modes perturbation
|
||||
!
|
||||
stmp = 'on'
|
||||
call get_parameter("perturb_tearing", stmp)
|
||||
select case(stmp)
|
||||
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
||||
perturb_tearing = .true.
|
||||
case default
|
||||
perturb_tearing = .false.
|
||||
end select
|
||||
|
||||
! prepare the wave vector of the perturbation
|
||||
!
|
||||
kvec = pi2 * kper
|
||||
@ -281,6 +291,35 @@ module user_problem
|
||||
!
|
||||
if (pert == 2) then
|
||||
|
||||
kl2 = max(0.0d+00, kper - kdel)**2
|
||||
ku2 = (kper + kdel)**2
|
||||
|
||||
kv = int(kper)
|
||||
ks = int(ceiling(1.0d+00 / (pi2 * dlta)))
|
||||
if (perturb_tearing) ks = 0
|
||||
n = 0
|
||||
#if NDIMS == 3
|
||||
kd = int(xlen / zlen)
|
||||
do kk = - kv + mod(kv, kd), kv - mod(kv, kd), kd
|
||||
#else /* NDIMS == 3 */
|
||||
kk = 0
|
||||
#endif /* NDIMS == 3 */
|
||||
do kj = - kv, kv
|
||||
do ki = ks, kv
|
||||
kv2 = 1.0d+00 * (ki**2 + kj**2 + kk**2)
|
||||
if (kl2 <= kv2 .and. kv2 <= ku2) then
|
||||
if (.not. (ki == 0 .and. kj == 0 .and. kk < 0) .and. &
|
||||
.not. (ki == 0 .and. kj < 0)) then
|
||||
n = n + 1
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
#if NDIMS == 3
|
||||
end do
|
||||
#endif /* NDIMS == 3 */
|
||||
if (n > 0) nper = n
|
||||
|
||||
! allocate phase and wave vector components
|
||||
!
|
||||
allocate(kx(nper), ky(nper), kz(nper), ux(nper), uy(nper), uz(nper), &
|
||||
@ -290,45 +329,52 @@ module user_problem
|
||||
|
||||
! choose random wave vector directions
|
||||
!
|
||||
kd = int(xlen / zlen)
|
||||
fc = 1.0d+00 / sqrt(1.0d+00 * nper)
|
||||
do n = 1, nper
|
||||
thh = pi2 * randuni()
|
||||
n = 0
|
||||
#if NDIMS == 3
|
||||
thv = pi * randsym()
|
||||
tx = cos(thh) * cos(thv)
|
||||
ty = sin(thh) * cos(thv)
|
||||
tz = sin(thv)
|
||||
kx(n) = pi2 * nint(kper * tx)
|
||||
ky(n) = pi2 * nint(kper * ty)
|
||||
kz(n) = pi2 * nint(kper * tz / kd) * kd
|
||||
tt = 0.0d+00
|
||||
do while(tt < 1.0d-08)
|
||||
thh = pi2 * randuni()
|
||||
thv = pi * randsym()
|
||||
tx = cos(thh) * cos(thv)
|
||||
ty = sin(thh) * cos(thv)
|
||||
tz = sin(thv)
|
||||
ux(n) = ty * kz(n) - tz * ky(n)
|
||||
uy(n) = tz * kx(n) - tx * kz(n)
|
||||
uz(n) = tx * ky(n) - ty * kx(n)
|
||||
tt = sqrt(ux(n)**2 + uy(n)**2 + uz(n)**2)
|
||||
end do
|
||||
ux(n) = fc * ux(n) / tt
|
||||
uy(n) = fc * uy(n) / tt
|
||||
uz(n) = fc * uz(n) / tt
|
||||
do kk = - kv + mod(kv, kd), kv - mod(kv, kd), kd
|
||||
#else /* NDIMS == 3 */
|
||||
kx(n) = pi2 * nint(kper * cos(thh))
|
||||
ky(n) = pi2 * nint(kper * sin(thh))
|
||||
kz(n) = 0.0d+00
|
||||
ka = sqrt(kx(n)**2 + ky(n)**2)
|
||||
ux(n) = fc * ky(n) / ka
|
||||
uy(n) = - fc * kx(n) / ka
|
||||
uz(n) = 0.0d+00
|
||||
kk = 0
|
||||
#endif /* NDIMS == 3 */
|
||||
ph(n) = pi2 * randuni()
|
||||
do kj = - kv, kv
|
||||
do ki = ks, kv
|
||||
kv2 = 1.0d+00 * (ki**2 + kj**2 + kk**2)
|
||||
if (kl2 <= kv2 .and. kv2 <= ku2) then
|
||||
if (.not. (ki == 0 .and. kj == 0 .and. kk < 0) .and. &
|
||||
.not. (ki == 0 .and. kj < 0)) then
|
||||
n = n + 1
|
||||
kx(n) = pi2 * ki
|
||||
ky(n) = pi2 * kj
|
||||
kz(n) = pi2 * kk
|
||||
if (kk == 0) then
|
||||
ka = sqrt(kx(n)**2 + ky(n)**2)
|
||||
ux(n) = ky(n) / ka
|
||||
uy(n) = - kx(n) / ka
|
||||
uz(n) = 0.0d+00
|
||||
else
|
||||
ka = 0.0d+00
|
||||
do while(ka < 1.0d-08)
|
||||
thh = pi2 * randuni()
|
||||
thv = pi * randsym()
|
||||
tx = cos(thh) * cos(thv)
|
||||
ty = sin(thh) * cos(thv)
|
||||
tz = sin(thv)
|
||||
ux(n) = ty * kz(n) - tz * ky(n)
|
||||
uy(n) = tz * kx(n) - tx * kz(n)
|
||||
uz(n) = tx * ky(n) - ty * kx(n)
|
||||
ka = sqrt(ux(n)**2 + uy(n)**2 + uz(n)**2)
|
||||
end do
|
||||
ux(n) = ux(n) / ka
|
||||
uy(n) = uy(n) / ka
|
||||
uz(n) = uz(n) / ka
|
||||
end if
|
||||
ph(n) = pi2 * randuni()
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
#if NDIMS == 3
|
||||
end do
|
||||
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
else
|
||||
if (verbose) &
|
||||
@ -446,12 +492,20 @@ module user_problem
|
||||
call print_parameter(verbose, '( ) P (Prandtl number)', prdl)
|
||||
call print_parameter(verbose, '(*) delta (thickness)', dlta)
|
||||
call print_parameter(verbose, '(*) perturbation', perturbation)
|
||||
if (pert == 2) then
|
||||
if (perturb_tearing) then
|
||||
call print_parameter(verbose, '(*) perturb tearing modes', 'on')
|
||||
else
|
||||
call print_parameter(verbose, '(*) perturb tearing modes', 'off')
|
||||
end if
|
||||
end if
|
||||
call print_parameter(verbose, '(*) bper' , bper)
|
||||
call print_parameter(verbose, '(*) vper' , vper)
|
||||
if (pert >= 1) then
|
||||
call print_parameter(verbose, '(*) kper' , kper)
|
||||
end if
|
||||
if (pert == 2) then
|
||||
call print_parameter(verbose, '(*) kdel' , kdel)
|
||||
call print_parameter(verbose, '(*) nper' , nper)
|
||||
end if
|
||||
call print_parameter(verbose, '(*) xcut' , xcut)
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user