diff --git a/src/problems.F90 b/src/problems.F90 index 360a22c..93a2aaf 100644 --- a/src/problems.F90 +++ b/src/problems.F90 @@ -147,9 +147,6 @@ module problems case("current_sheet") setup_problem => setup_problem_current_sheet - case("reconnection") - setup_problem => setup_problem_reconnection - case("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: ! ---------------------------- ! diff --git a/src/user_problem.F90 b/src/user_problem.F90 index 13624a6..084a551 100644 --- a/src/user_problem.F90 +++ b/src/user_problem.F90 @@ -45,8 +45,19 @@ module user_problem integer, save :: imi, imp, ims, imu, img, imb #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 ! @@ -75,7 +86,7 @@ module user_problem ! 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 ! @@ -111,6 +122,23 @@ module user_problem ! 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 ! parameters, etc. ! @@ -192,10 +220,15 @@ module user_problem ! 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 random , only : randomn ! local variables are not implicit by default ! @@ -207,7 +240,8 @@ module user_problem ! local variables ! - integer :: i, j, k + integer :: i, j, k + real(kind=8) :: yt, yp ! local arrays ! @@ -215,6 +249,8 @@ module user_problem real(kind=8), dimension(im) :: x real(kind=8), dimension(jm) :: y 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) #endif /* PROFILE */ -! set the variables +! prepare block coordinates ! - q(idn,:) = 1.0d+00 - if (ipr > 0) q(ipr,:) = 1.0d+00 - q(ivx,:) = 0.0d+00 - q(ivy,:) = 0.0d+00 - q(ivz,:) = 0.0d+00 - if (ibx > 0) then - q(ibx,:) = 0.0d+00 - q(iby,:) = 0.0d+00 - q(ibz,:) = 0.0d+00 - q(ibp,:) = 0.0d+00 - end if + 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(im, q(1:nv,1:im), u(1:nv,1:im)) - -! iterate over all positions in the YZ plane -! - do k = 1, km - do j = 1, jm + call prim2cons(jm, q(1:nv,1:jm), u(1:nv,1:jm)) ! 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 ! - 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 #ifdef PROFILE