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

USER_PROBLEM: Implement a few types of the initial perturbation.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2024-08-17 17:55:22 -03:00
parent b24ab38389
commit 4429843c95
1 changed files with 339 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

360
sources/user_problem.F90
View File

@ -44,8 +44,22 @@ module user_problem
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 :: prdl = 0.00d+00
real(kind=8), save :: nu = 0.00d+00
real(kind=8), save :: eta = 0.00d+00
integer , save :: pert = -1
integer , save :: nper = 0
real(kind=8), save :: vper = 0.00d+00
real(kind=8), save :: kper = 1.00d+00
real(kind=8), save :: kdel = 1.00d+00
real(kind=8), save :: xcut = 1.00d+99
real(kind=8), save :: ycut = 1.00d+99
real(kind=8), save :: xdec = 1.00d+00
real(kind=8), save :: ydec = 1.00d+00
real(kind=8), dimension(:), allocatable :: kx, ky, kz, ux, uy, uz, ph
integer(kind=4), save :: runit = 33
logical, save :: resistive = .false.
@ -78,10 +92,13 @@ module user_problem
subroutine initialize_user_problem(problem, rcount, verbose, status)
use boundaries , only : custom_boundary_x, custom_boundary_y
use constants , only : pi, pi2
use coordinates, only : xlen, zlen
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
use random , only : reset_seeds, randuni, randsym
implicit none
@ -90,15 +107,19 @@ module user_problem
logical , intent(in) :: verbose
integer , intent(out) :: status
character(len=64) :: tmp, append = "off", fname
logical :: flag
character(len=64) :: str_param, append = "off", fname
logical :: flag, perturb_tearing
integer :: n, kd, ki, kj, kk, kv, ks
real(kind=8) :: kl2, ku2, kv2, ka
real(kind=8) :: thv, thh, tx, ty, tz
character(len=*), parameter :: &
loc = 'USER_PROBLEM::initialize_user_problem()'
!-------------------------------------------------------------------------------
!
call get_parameter("profile", tmp)
select case(tmp)
str_param = "sincos"
call get_parameter("profile", str_param)
select case(str_param)
case('bhattacharjee', 'sincos', 'trigonometric')
profile = 2
case default
@ -128,6 +149,7 @@ module user_problem
status = 1
return
end if
call get_parameter("viscosity" , nu )
call get_parameter("resistivity", eta)
if (eta < 0.0d+00) then
if (verbose) &
@ -146,6 +168,8 @@ module user_problem
status = 1
return
end if
ycut = 4 * dlta
ydec = dlta
pmag = 0.5d+00 * (bamp**2 + bgui**2)
pres = beta * pmag
@ -158,6 +182,142 @@ module user_problem
csnd = sqrt(csnd2)
valf = bamp / sqrt(dens)
lund = valf / max(tiny(eta), eta)
prdl = nu / max(tiny(eta), eta)
! get the perturbation parameters
!
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)
call get_parameter("ydec", ydec)
! choose the perturbation type
!
str_param = "none"
call get_parameter("perturbation", str_param)
if (vper == 0.0d+00) str_param = "none"
select case(str_param)
case('noise', 'random')
pert = 0
call reset_seeds('random')
case('mode', 'one mode', 'one-mode', 'one_mode')
pert = 1
case('multi-wave', 'random waves', 'random-waves', 'random_waves')
pert = 2
case default
pert = -1
end select
! enable or disable tearing modes perturbation
!
str_param = "on"
call get_parameter("perturb_tearing", str_param)
select case(str_param)
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
perturb_tearing = .true.
case default
perturb_tearing = .false.
end select
! prepare the wave vectors for multi-wave perturbation
!
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), &
ph(nper), stat = status)
if (status == 0) then
! choose random wave vector directions
!
n = 0
#if NDIMS == 3
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
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) &
call print_message(loc, &
"Could not allocate space for perturbation vectors!")
return
end if
end if ! pert = 2
! determine if to append or create another file
!
@ -215,13 +375,37 @@ module user_problem
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)
call print_parameter(verbose, '( ) Pr (Prandtl number)' , prdl)
call print_parameter(verbose, '(*) delta (thickness)' , dlta)
select case(pert)
case(0)
call print_parameter(verbose, '(*) perturbation', 'noise')
case(1)
call print_parameter(verbose, '(*) perturbation', 'one mode')
case(2)
call print_parameter(verbose, '(*) perturbation', 'waves at given k')
end select
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, '(*) vper' , vper)
if (pert == 2) then
call print_parameter(verbose, '(*) kper' , kper)
call print_parameter(verbose, '(*) kdel' , kdel)
call print_parameter(verbose, '(*) nper' , nper)
call print_parameter(verbose, '(*) ycut' , ycut)
call print_parameter(verbose, '(*) ydec' , ydec)
end if
setup_problem => setup_user_problem
custom_boundary_x => user_boundary_x
@ -248,14 +432,25 @@ module user_problem
!
subroutine finalize_user_problem(status)
use helpers, only : print_message
implicit none
integer, intent(out) :: status
character(len=*), parameter :: loc = 'USER_PROBLEM::finalize_user_problem()'
!-------------------------------------------------------------------------------
!
status = 0
if (pert == 2) then
deallocate(kx, ky, kz, ux, uy, uz, ph, stat = status)
if (status /= 0) &
call print_message(loc, &
"Could not deallocate space for perturbation vectors!")
end if
!-------------------------------------------------------------------------------
!
end subroutine finalize_user_problem
@ -279,25 +474,46 @@ module user_problem
use blocks , only : block_data
use constants , only : pi, pi2
use coordinates, only : nn => bcells, xmin, xmax
#if NDIMS == 3
use coordinates, only : ax, ay, az
#else /* NDIMS == 3 */
use coordinates, only : ax, ay
#endif /* NDIMS == 3 */
use equations , only : prim2cons, csnd2
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
use random , only : randnorz
implicit none
type(block_data), pointer, intent(inout) :: pdata
integer :: i, j, k
real(kind=8) :: r1, r2, bx, by, bz
logical :: flag
integer :: i, j, k, n
real(kind=8) :: r1, r2, qx, qy, qz, fv, xp, yp, kv, va
complex(kind=8) :: c
real(kind=8), dimension(3) :: q
real(kind=8), dimension(6) :: t
#if NDIMS == 3
real(kind=8), dimension(nn) :: x, y, z
#else /* NDIMS == 3 */
real(kind=8), dimension(nn) :: x, y
#endif /* NDIMS == 3 */
real(kind=8), dimension(nn) :: fx, fy
real(kind=8), dimension(nn) :: snx, csx, sny, csy, thy, shy
real(kind=8), dimension(nn,nn) :: az
real(kind=8), dimension(nn,nn) :: ps
!-------------------------------------------------------------------------------
!
k = 1
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
#if NDIMS == 3
z(:) = pdata%meta%zmin + az(pdata%meta%level,:)
#endif /* NDIMS == 3 */
! choose profile
!
select case(profile)
case(2)
@ -313,13 +529,13 @@ module user_problem
do j = 1, nn
do i = 1, nn
bx = bamp * csx(i) * (csy(j) * thy(j) + sny(j) * shy(j) / pi2)
by = bamp * 5.0d-01 * snx(i) * sny(j) * thy(j)
bz = sqrt(max(0.0d+00, bgui**2 + zeta * (bamp**2 - bx**2 - by**2)))
qx = bamp * csx(i) * (csy(j) * thy(j) + sny(j) * shy(j) / pi2)
qy = bamp * 5.0d-01 * snx(i) * sny(j) * thy(j)
qz = sqrt(max(0.0d+00, bgui**2 + zeta * (bamp**2 - qx**2 - qy**2)))
pdata%q(ibx,i,j,:) = bx
pdata%q(iby,i,j,:) = by
pdata%q(ibz,i,j,:) = bz
pdata%q(ibx,i,j,:) = qx
pdata%q(iby,i,j,:) = qy
pdata%q(ibz,i,j,:) = qz
end do
end do
pdata%q(ibp,:,:,:) = 0.0d+00
@ -329,20 +545,20 @@ module user_problem
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) &
ps(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)))
qx = (ps(i,j+1) - ps(i,j-1)) / (y(j+1) - y(j-1))
qy = (ps(i-1,j) - ps(i+1,j)) / (x(i+1) - x(i-1))
qz = sqrt(max(0.0d+00, bgui**2 + zeta * (bamp**2 - qx**2 - qy**2)))
pdata%q(ibx,i,j,:) = bx
pdata%q(iby,i,j,:) = by
pdata%q(ibz,i,j,:) = bz
pdata%q(ibx,i,j,:) = qx
pdata%q(iby,i,j,:) = qy
pdata%q(ibz,i,j,:) = qz
end do
end do
pdata%q(ibp,:,:,:) = 0.0d+00
@ -361,6 +577,108 @@ module user_problem
pdata%q(ivy,:,:,:) = 0.0d+00
pdata%q(ivz,:,:,:) = 0.0d+00
! prepare decaying factors
!
if (pert >= 0) then
fv = 5.0d-01 * pi
do i = 1, nn
xp = fv * min(1.0d+00, max(0.0d+00, abs(x(i)) - xcut) / xdec)
fx(i) = cos(xp)**2
end do ! i = 1, nn
do j = 1, nn
yp = fv * min(1.0d+00, max(0.0d+00, abs(y(j)) - ycut) / ydec)
fy(j) = cos(yp)**2
end do ! i = 1, nn
end if ! pert >= 0
! add some noise to the velocity components
!
if (pert == 0) then
fv = vper / sqrt(2.0d+00)
flag = .true.
#if NDIMS == 3
do k = 1, nn
#endif /* NDIMS == 3 */
do j = 1, nn
if (fy(j) > 0.0d+00) then
do i = 1, nn
if (fx(i) > 0.0d+00) then
va = fv * fx(i) * fy(j)
if (flag) then
c = randnorz()
t(1:2) = [ real(c), aimag(c) ]
c = randnorz()
t(3:4) = [ real(c), aimag(c) ]
c = randnorz()
t(5:6) = [ real(c), aimag(c) ]
q(1:3) = t(1:3)
else
q(1:3) = t(4:6)
end if
flag = .not. flag
pdata%q(ivx,i,j,k) = va * q(1)
pdata%q(ivy,i,j,k) = va * q(2)
pdata%q(ivz,i,j,k) = va * q(3)
end if ! fx > 0.0
end do ! i = 1, nn
end if ! fy > 0.0
end do ! j = 1, nn
#if NDIMS == 3
end do ! k = 1, nn
#endif /* NDIMS == 3 */
end if ! pert == 0
! add the multi-wave perturbation of velocity
!
if (pert == 2) then
#if NDIMS == 3
do k = 1, nn
#endif /* NDIMS == 3 */
do j = 1, nn
if (fy(j) > 0.0d+00) then
do i = 1, nn
if (fx(i) > 0.0d+00) then
fv = vper * fx(i) * fy(j)
do n = 1, nper
#if NDIMS == 3
kv = kx(n) * x(i) + ky(n) * y(j) + kz(n) * z(k) + ph(n)
#else /* NDIMS == 3 */
kv = kx(n) * x(i) + ky(n) * y(j) + ph(n)
#endif /* NDIMS == 3 */
va = fv * sin(kv)
pdata%q(ivx,i,j,k) = pdata%q(ivx,i,j,k) + va * ux(n)
pdata%q(ivy,i,j,k) = pdata%q(ivy,i,j,k) + va * uy(n)
#if NDIMS == 3
pdata%q(ivz,i,j,k) = pdata%q(ivz,i,j,k) + va * uz(n)
#endif /* NDIMS == 3 */
end do
end if ! fx > 0.0
end do ! i = 1, nn
end if ! fy > 0.0
end do ! j = 1, nn
#if NDIMS == 3
end do ! k = 1, nn
#endif /* NDIMS == 3 */
end if ! pert == 2
#if NDIMS == 3
do k = 1, nn
#else /* NDIMS == 3 */