PROBLEMS: Move reconnection problem to USER_PROBLEM module.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2017-03-08 13:11:48 -03:00
parent 2ec54cb0f1
commit 6a9ef98fd6
2 changed files with 152 additions and 266 deletions

View File

@ -147,9 +147,6 @@ module problems
case("current_sheet") case("current_sheet")
setup_problem => setup_problem_current_sheet setup_problem => setup_problem_current_sheet
case("reconnection")
setup_problem => setup_problem_reconnection
case("jet") case("jet")
setup_problem => setup_problem_jet setup_problem => setup_problem_jet
@ -1556,245 +1553,6 @@ module problems
! !
!=============================================================================== !===============================================================================
! !
! subroutine SETUP_PROBLEM_RECONNECTION:
! -------------------------------------
!
! Subroutine sets the initial conditions for the reconnection problem.
!
! Arguments:
!
! pdata - pointer to the datablock structure of the currently initialized
! block;
!
!===============================================================================
!
subroutine setup_problem_reconnection(pdata)
! include external procedures and variables
!
use blocks , only : block_data
use constants , only : pi2
use coordinates, only : im, jm, km
use coordinates, only : ax, ay, adx, ady, adz
use equations , only : prim2cons
use equations , only : nv
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
use equations , only : csnd2
use operators , only : curl
use parameters , only : get_parameter_real
use random , only : randomn
! 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 :: dens = 1.00d+00
real(kind=8), save :: pres = 1.00d+00
real(kind=8), save :: bamp = 1.00d+00
real(kind=8), save :: bper = 0.00d+00
real(kind=8), save :: bgui = 0.00d+00
real(kind=8), save :: vper = 1.00d-02
real(kind=8), save :: xcut = 4.00d-01
real(kind=8), save :: ycut = 1.00d-02
real(kind=8), save :: yth = 1.00d-16
real(kind=8), save :: pth = 1.00d-02
real(kind=8), save :: pmag = 5.00d-01
! local saved parameters
!
logical , save :: first = .true.
! local variables
!
integer :: i, j, k
real(kind=8) :: yt, yp
! local arrays
!
real(kind=8), dimension(nv,jm) :: q, u
real(kind=8), dimension(im) :: x
real(kind=8), dimension(jm) :: y, pm
real(kind=8), dimension(3) :: dh
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for the problem setup
!
call start_timer(imu)
#endif /* PROFILE */
! prepare problem constants during the first subroutine call
!
if (first) then
! get problem parameters
!
call get_parameter_real("dens", dens)
call get_parameter_real("pres", pres)
call get_parameter_real("bamp", bamp)
call get_parameter_real("bper", bper)
call get_parameter_real("bgui", bgui)
call get_parameter_real("vper", vper)
call get_parameter_real("xcut", xcut)
call get_parameter_real("ycut", ycut)
call get_parameter_real("yth" , yth )
call get_parameter_real("pth" , pth )
! calculate the maximum magnetic pressure
!
pmag = 0.5d+00 * (bamp**2 + bgui**2)
! reset the first execution flag
!
first = .false.
end if ! first call
! prepare block coordinates
!
x(1:im) = pdata%meta%xmin + ax(pdata%meta%level,1:im)
y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm)
! prepare cell sizes
!
dh(1) = adx(pdata%meta%level)
dh(2) = ady(pdata%meta%level)
dh(3) = adz(pdata%meta%level)
! calculate the perturbation of magnetic field
!
if (bper /= 0.0d+00) then
! initiate the vector potential (we use velocity components to store vector
! potential temporarily, and we store magnetic field perturbation in U)
!
do k = 1, km
do j = 1, jm
do i = 1, im
yt = abs(y(j) / yth)
yt = log(exp(yt) + exp(- yt))
yp = y(j) / pth
pdata%q(ivx,i,j,k) = 0.0d+00
pdata%q(ivy,i,j,k) = 0.0d+00
pdata%q(ivz,i,j,k) = bper * cos(pi2 * x(i)) * exp(- yp * yp) / pi2
end do ! i = 1, im
end do ! j = 1, jm
end do ! k = 1, km
! calculate magnetic field components from vector potential
!
call curl(dh(1:3), pdata%q(ivx:ivz,1:im,1:jm,1:km) &
, pdata%q(ibx:ibz,1:im,1:jm,1:km))
else
! reset magnetic field components
!
pdata%q(ibx:ibz,:,:,:) = 0.0d+00
end if ! bper /= 0.0
! iterate over all positions in the XZ plane
!
do k = 1, km
do i = 1, im
! if magnetic field is present, set its initial configuration
!
if (ibx > 0) then
! set antiparallel magnetic field component
!
do j = 1, jm
q(ibx,j) = bamp * tanh(y(j) / yth)
end do ! j = 1, jm
! set tangential magnetic field components
!
q(iby,1:jm) = 0.0d+00
q(ibz,1:jm) = bgui
q(ibp,1:jm) = 0.0d+00
! calculate local magnetic pressure
!
pm(1:jm) = 0.5d+00 * sum(q(ibx:ibz,:) * q(ibx:ibz,:), 1)
! add magnetic field perturbation
!
if (bper /= 0.0d+00) then
q(ibx,1:jm) = q(ibx,1:jm) + pdata%q(ibx,i,1:jm,k)
q(iby,1:jm) = q(iby,1:jm) + pdata%q(iby,i,1:jm,k)
q(ibz,1:jm) = q(ibz,1:jm) + pdata%q(ibz,i,1:jm,k)
end if ! bper /= 0.0
end if ! ibx > 0
! set the uniform density and pressure
!
if (ipr > 0) then
q(idn,1:jm) = dens
q(ipr,1:jm) = pres + (pmag - pm(1:jm))
else
q(idn,1:jm) = dens + (pmag - pm(1:jm)) / csnd2
end if
! reset velocity components
!
q(ivx,1:jm) = 0.0d+00
q(ivy,1:jm) = 0.0d+00
q(ivz,1:jm) = 0.0d+00
! set the random velocity field in a layer near current sheet
!
if (abs(x(i)) <= xcut) then
do j = 1, jm
if (abs(y(j)) <= ycut) then
q(ivx,j) = vper * randomn()
q(ivy,j) = vper * randomn()
#if NDIMS == 3
q(ivz,j) = vper * randomn()
#endif /* NDIMS == 3 */
end if ! |y| < ycut
end do ! j = 1, jm
end if ! |x| < xcut
! convert the primitive variables to conservative ones
!
call prim2cons(jm, q(1:nv,1:jm), u(1:nv,1:jm))
! copy the conserved variables to the current block
!
pdata%u(1:nv,i,1:jm,k) = u(1:nv,1:jm)
! copy the primitive variables to the current block
!
pdata%q(1:nv,i,1:jm,k) = q(1:nv,1:jm)
end do ! i = 1, im
end do ! k = 1, km
#ifdef PROFILE
! stop accounting time for the problems setup
!
call stop_timer(imu)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine setup_problem_reconnection
!
!===============================================================================
!
! subroutine SETUP_PROBLEM_JET: ! subroutine SETUP_PROBLEM_JET:
! ---------------------------- ! ----------------------------
! !

View File

@ -45,8 +45,19 @@ module user_problem
integer, save :: imi, imp, ims, imu, img, imb integer, save :: imi, imp, ims, imu, img, imb
#endif /* PROFILE */ #endif /* PROFILE */
! default problem parameter values are defined here ! default problem parameter values
! !
real(kind=8), save :: dens = 1.00d+00
real(kind=8), save :: pres = 1.00d+00
real(kind=8), save :: bamp = 1.00d+00
real(kind=8), save :: bper = 0.00d+00
real(kind=8), save :: bgui = 0.00d+00
real(kind=8), save :: vper = 1.00d-02
real(kind=8), save :: xcut = 4.00d-01
real(kind=8), save :: ycut = 1.00d-02
real(kind=8), save :: yth = 1.00d-16
real(kind=8), save :: pth = 1.00d-02
real(kind=8), save :: pmag = 5.00d-01
! flag indicating if the gravitational source term is enabled ! flag indicating if the gravitational source term is enabled
! !
@ -75,7 +86,7 @@ module user_problem
! include external procedures and variables ! include external procedures and variables
! !
use parameters, only : get_parameter_string use parameters, only : get_parameter_string, get_parameter_real
! local variables are not implicit by default ! local variables are not implicit by default
! !
@ -111,6 +122,23 @@ module user_problem
! !
call get_parameter_string("problem", problem_name) call get_parameter_string("problem", problem_name)
! get the reconnection problem parameters
!
call get_parameter_real("dens", dens)
call get_parameter_real("pres", pres)
call get_parameter_real("bamp", bamp)
call get_parameter_real("bper", bper)
call get_parameter_real("bgui", bgui)
call get_parameter_real("vper", vper)
call get_parameter_real("xcut", xcut)
call get_parameter_real("ycut", ycut)
call get_parameter_real("yth" , yth )
call get_parameter_real("pth" , pth )
! calculate the maximum magnetic pressure
!
pmag = 0.5d+00 * (bamp**2 + bgui**2)
! print information about the user problem such as problem name, its ! print information about the user problem such as problem name, its
! parameters, etc. ! parameters, etc.
! !
@ -192,10 +220,15 @@ module user_problem
! include external procedures and variables ! include external procedures and variables
! !
use blocks , only : block_data use blocks , only : block_data
use constants , only : pi2
use coordinates, only : im, jm, km use coordinates, only : im, jm, km
use coordinates, only : ax, ay, adx, ady, adz
use equations , only : prim2cons use equations , only : prim2cons
use equations , only : nv use equations , only : nv
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
use equations , only : csnd2
use operators , only : curl
use random , only : randomn
! local variables are not implicit by default ! local variables are not implicit by default
! !
@ -207,7 +240,8 @@ module user_problem
! local variables ! local variables
! !
integer :: i, j, k integer :: i, j, k
real(kind=8) :: yt, yp
! local arrays ! local arrays
! !
@ -215,6 +249,8 @@ module user_problem
real(kind=8), dimension(im) :: x real(kind=8), dimension(im) :: x
real(kind=8), dimension(jm) :: y real(kind=8), dimension(jm) :: y
real(kind=8), dimension(km) :: z real(kind=8), dimension(km) :: z
real(kind=8), dimension(jm) :: pm
real(kind=8), dimension(3) :: dh
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -224,38 +260,130 @@ module user_problem
call start_timer(imp) call start_timer(imp)
#endif /* PROFILE */ #endif /* PROFILE */
! set the variables ! prepare block coordinates
! !
q(idn,:) = 1.0d+00 x(1:im) = pdata%meta%xmin + ax(pdata%meta%level,1:im)
if (ipr > 0) q(ipr,:) = 1.0d+00 y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm)
q(ivx,:) = 0.0d+00
q(ivy,:) = 0.0d+00 ! prepare cell sizes
q(ivz,:) = 0.0d+00 !
if (ibx > 0) then dh(1) = adx(pdata%meta%level)
q(ibx,:) = 0.0d+00 dh(2) = ady(pdata%meta%level)
q(iby,:) = 0.0d+00 dh(3) = adz(pdata%meta%level)
q(ibz,:) = 0.0d+00
q(ibp,:) = 0.0d+00 ! calculate the perturbation of magnetic field
end if !
if (bper /= 0.0d+00) then
! initiate the vector potential (we use velocity components to store vector
! potential temporarily, and we store magnetic field perturbation in U)
!
do k = 1, km
do j = 1, jm
do i = 1, im
yt = abs(y(j) / yth)
yt = log(exp(yt) + exp(- yt))
yp = y(j) / pth
pdata%q(ivx,i,j,k) = 0.0d+00
pdata%q(ivy,i,j,k) = 0.0d+00
pdata%q(ivz,i,j,k) = bper * cos(pi2 * x(i)) * exp(- yp * yp) / pi2
end do ! i = 1, im
end do ! j = 1, jm
end do ! k = 1, km
! calculate magnetic field components from vector potential
!
call curl(dh(1:3), pdata%q(ivx:ivz,1:im,1:jm,1:km) &
, pdata%q(ibx:ibz,1:im,1:jm,1:km))
else
! reset magnetic field components
!
pdata%q(ibx:ibz,:,:,:) = 0.0d+00
end if ! bper /= 0.0
! iterate over all positions in the XZ plane
!
do k = 1, km
do i = 1, im
! if magnetic field is present, set its initial configuration
!
if (ibx > 0) then
! set antiparallel magnetic field component
!
do j = 1, jm
q(ibx,j) = bamp * tanh(y(j) / yth)
end do ! j = 1, jm
! set tangential magnetic field components
!
q(iby,1:jm) = 0.0d+00
q(ibz,1:jm) = bgui
q(ibp,1:jm) = 0.0d+00
! calculate local magnetic pressure
!
pm(1:jm) = 0.5d+00 * sum(q(ibx:ibz,:) * q(ibx:ibz,:), 1)
! add magnetic field perturbation
!
if (bper /= 0.0d+00) then
q(ibx,1:jm) = q(ibx,1:jm) + pdata%q(ibx,i,1:jm,k)
q(iby,1:jm) = q(iby,1:jm) + pdata%q(iby,i,1:jm,k)
q(ibz,1:jm) = q(ibz,1:jm) + pdata%q(ibz,i,1:jm,k)
end if ! bper /= 0.0
end if ! ibx > 0
! set the uniform density and pressure
!
if (ipr > 0) then
q(idn,1:jm) = dens
q(ipr,1:jm) = pres + (pmag - pm(1:jm))
else
q(idn,1:jm) = dens + (pmag - pm(1:jm)) / csnd2
end if
! reset velocity components
!
q(ivx,1:jm) = 0.0d+00
q(ivy,1:jm) = 0.0d+00
q(ivz,1:jm) = 0.0d+00
! set the random velocity field in a layer near current sheet
!
if (abs(x(i)) <= xcut) then
do j = 1, jm
if (abs(y(j)) <= ycut) then
q(ivx,j) = vper * randomn()
q(ivy,j) = vper * randomn()
#if NDIMS == 3
q(ivz,j) = vper * randomn()
#endif /* NDIMS == 3 */
end if ! |y| < ycut
end do ! j = 1, jm
end if ! |x| < xcut
! convert the primitive variables to conservative ones ! convert the primitive variables to conservative ones
! !
call prim2cons(im, q(1:nv,1:im), u(1:nv,1:im)) call prim2cons(jm, q(1:nv,1:jm), u(1:nv,1:jm))
! iterate over all positions in the YZ plane
!
do k = 1, km
do j = 1, jm
! copy the primitive variables to the current block ! copy the primitive variables to the current block
! !
pdata%q(1:nv,1:im,j,k) = q(1:nv,1:im) pdata%q(1:nv,i,1:jm,k) = q(1:nv,1:jm)
! copy the conserved variables to the current block ! copy the conserved variables to the current block
! !
pdata%u(1:nv,1:im,j,k) = u(1:nv,1:im) pdata%u(1:nv,i,1:jm,k) = u(1:nv,1:jm)
end do ! j = 1, jm end do ! i = 1, im
end do ! k = 1, km end do ! k = 1, km
#ifdef PROFILE #ifdef PROFILE