From b6d86db3dcabcf754faf92bff912abec95c64e29 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 17 Oct 2022 11:09:30 -0300 Subject: [PATCH] USER_PROBLEM: Implement setup for GEM magnetic reconnection challenge. See, e.g., Birn et al., Journal of Geophysical Research, 2001, 106, A3, 3715-3719, DOI:10.1029/1999JA900449. Signed-off-by: Grzegorz Kowal --- sources/user_problem.F90 | 180 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 179 insertions(+), 1 deletion(-) diff --git a/sources/user_problem.F90 b/sources/user_problem.F90 index ed4fc60..939b075 100644 --- a/sources/user_problem.F90 +++ b/sources/user_problem.F90 @@ -31,6 +31,22 @@ module user_problem implicit none + real(kind=8), save :: beta = 2.00d+00 + real(kind=8), save :: zeta = 0.00d+00 + real(kind=8), save :: dens = 1.00d+00 + real(kind=8), save :: bamp = 1.00d+00 + real(kind=8), save :: bgui = 0.00d+00 + real(kind=8), save :: psi0 = 1.00d-01 + real(kind=8), save :: dlta = 5.00d-01 + real(kind=8), save :: pres = 5.00d-01 + real(kind=8), save :: pmag = 5.00d-01 + 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 :: eta = 0.00d+00 + + logical, save :: resistive = .false. + private public :: initialize_user_problem, finalize_user_problem public :: user_statistics @@ -58,6 +74,12 @@ module user_problem ! subroutine initialize_user_problem(problem, rcount, verbose, status) + use boundaries , only : custom_boundary_y + 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 + implicit none character(len=64), intent(in) :: problem @@ -65,8 +87,83 @@ module user_problem logical , intent(in) :: verbose integer , intent(out) :: status + character(len=*), parameter :: & + loc = 'USER_PROBLEM::initialize_user_problem()' + !------------------------------------------------------------------------------- ! + call get_parameter("dens", dens) + call get_parameter("bamp", bamp) + call get_parameter("bgui", bgui) + call get_parameter("beta", beta) + if (beta <= 0.0d+00) then + if (verbose) & + call print_message(loc, "Plasma-beta must be positive (beta > 0.0)!") + status = 1 + return + end if + call get_parameter("psi0", psi0) + call get_parameter("zeta", zeta) + if (zeta < 0.0d+00 .or. zeta > 1.0d+00) then + if (verbose) & + call print_message(loc, "Parameter 'zeta' must be between 0.0 and 1.0!") + status = 1 + return + end if + if (dens <= 0.0d+00) then + if (verbose) & + call print_message(loc, "Density must be positive (dens > 0.0)!") + status = 1 + return + end if + call get_parameter("delta", dlta) + if (dlta <= 0.0d+00) then + if (verbose) & + call print_message(loc, & + "The initial thickness must be positive (delta > 0.0)!") + status = 1 + return + end if + call get_parameter("resistivity", eta) + if (eta < 0.0d+00) then + if (verbose) & + call print_message(loc, "Resistivity cannot be negative!") + status = 1 + return + else + resistive = .true. + end if + + pmag = 5.0d-01 * (bamp**2 + bgui**2) + pres = beta * pmag + ptot = pres + pmag + if (ipr > 0) then + csnd2 = adiabatic_index * pres / dens + else + csnd2 = pres / dens + end if + csnd = sqrt(csnd2) + valf = sqrt(2.0d+00 * pmag / dens) + lund = valf / max(tiny(eta), eta) + + call print_section(verbose, "Parameters (* - set, otherwise calculated)") + call print_parameter(verbose, '(*) beta (plasma-beta)' , beta) + call print_parameter(verbose, '(*) zeta' , zeta) + call print_parameter(verbose, '(*) phi0' , psi0) + 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) + + setup_problem => setup_user_problem + custom_boundary_y => user_boundary_y + status = 0 !------------------------------------------------------------------------------- @@ -116,14 +213,73 @@ module user_problem ! subroutine setup_user_problem(pdata) - use blocks, only : block_data + use blocks , only : block_data + use constants , only : pi, pi2 + use coordinates, only : nn => bcells, xlen, ylen + use coordinates, only : ax, ay + use equations , only : prim2cons, csnd2 + use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp implicit none type(block_data), pointer, intent(inout) :: pdata + integer :: i, j, k + real(kind=8), dimension(nn) :: x, y + real(kind=8), dimension(nn) :: snx, csx, sny, csy, tgy + !------------------------------------------------------------------------------- ! + x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) + y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) + +! equilibrium +! + tgy(:) = bamp * tanh(y(:) / dlta) + do j = 1, nn + do i = 1, nn + pdata%q(ibx,i,j,:) = tgy(j) + pdata%q(ibz,i,j,:) = sqrt(bgui**2 & + + zeta**2 * max(0.0d+00, bamp**2 - tgy(j)**2)) + end do + end do + pdata%q(iby,:,:,:) = 0.0d+00 + pdata%q(ibp,:,:,:) = 0.0d+00 + + if (ipr > 0) then + pdata%q(ipr,:,:,:) = ptot - 5.0d-01 * sum(pdata%q(ibx:ibz,:,:,:)**2,1) + pdata%q(idn,:,:,:) = (dens / pres) * pdata%q(ipr,:,:,:) + else + 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 + +! perturbation +! + snx(:) = sin(pi2 * x(:) / xlen) * pi2 / xlen + csx(:) = cos(pi2 * x(:) / xlen) + sny(:) = sin(pi * y(:) / ylen) * pi / ylen + csy(:) = cos(pi * y(:) / ylen) + + do j = 1, nn + do i = 1, nn + pdata%q(ibx,i,j,:) = pdata%q(ibx,i,j,:) - psi0 * csx(i) * sny(j) + pdata%q(iby,i,j,:) = pdata%q(iby,i,j,:) + psi0 * snx(i) * csy(j) + end do + end do + +#if NDIMS == 3 + do k = 1, nn +#else /* NDIMS == 3 */ + do k = 1, 1 +#endif /* NDIMS == 3 */ + do j = 1, nn + call prim2cons(pdata%q(:,:,j,k), pdata%u(:,:,j,k), .true.) + end do + end do !------------------------------------------------------------------------------- ! @@ -247,6 +403,9 @@ module user_problem ! subroutine user_boundary_y(js, il, iu, kl, ku, t, dt, x, y, z, qn) + use coordinates, only : nn => bcells, nb, ne, nbl, neu + use equations , only : ivy, iby + implicit none integer , intent(in) :: js, il, iu, kl, ku @@ -254,8 +413,27 @@ module user_problem real(kind=8), dimension(:) , intent(in) :: x, y, z real(kind=8), dimension(:,:,:,:), intent(inout) :: qn + integer :: i, j + !------------------------------------------------------------------------------- ! + if (js == 1) then + do j = nbl, 1, -1 + do i = il, iu + qn( : ,i,j,:) = qn( : ,i,nb,:) + qn(ivy,i,j,:) = 0.0d+00 + qn(iby,i,j,:) = 0.0d+00 + end do + end do + else + do j = neu, nn + do i = il, iu + qn( : ,i,j,:) = qn( : ,i,ne,:) + qn(ivy,i,j,:) = 0.0d+00 + qn(iby,i,j,:) = 0.0d+00 + end do + end do + end if !------------------------------------------------------------------------------- !