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>
This commit is contained in:
parent
f7daa3f60d
commit
80f1fdde06
@ -31,6 +31,8 @@ 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
|
||||
@ -88,13 +90,21 @@ module user_problem
|
||||
logical , intent(in) :: verbose
|
||||
integer , intent(out) :: status
|
||||
|
||||
character(len=64) :: append = "off", fname
|
||||
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)
|
||||
@ -192,6 +202,14 @@ module user_problem
|
||||
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)
|
||||
@ -260,7 +278,7 @@ module user_problem
|
||||
|
||||
use blocks , only : block_data
|
||||
use constants , only : pi, pi2
|
||||
use coordinates, only : nn => bcells
|
||||
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
|
||||
@ -270,34 +288,68 @@ module user_problem
|
||||
type(block_data), pointer, intent(inout) :: pdata
|
||||
|
||||
integer :: i, j, k
|
||||
real(kind=8) :: bx, by, bz
|
||||
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,:)
|
||||
|
||||
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
|
||||
select case(profile)
|
||||
|
||||
do j = 1, nn
|
||||
do i = 1, nn
|
||||
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)
|
||||
bz = sqrt(max(0.0d+00, bgui**2 + zeta * (bamp**2 - bx**2 - by**2)))
|
||||
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
|
||||
|
||||
pdata%q(ibx,i,j,:) = bx
|
||||
pdata%q(iby,i,j,:) = by
|
||||
pdata%q(ibz,i,j,:) = bz
|
||||
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
|
||||
end do
|
||||
pdata%q(ibp,:,:,:) = 0.0d+00
|
||||
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
|
||||
@ -306,6 +358,7 @@ module user_problem
|
||||
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
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user