Merge branch 'master' into gem-reconnection-challenge

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>

problems/gem.in

@@ -0,0 +1,58 @@
+# problem name and parameters
+#
+problem              = "gem"
+
+# diffusion coefficients
+#
+viscosity            = 1.00d-03
+resistivity          = 1.00d-03
+glm_source_terms     = "kepes"
+
+# physics
+#
+equation_system      = "mhd"
+equation_of_state    = "iso"
+
+# methods
+#
+time_advance         = "ssprk3(2)4"
+riemann_solver       = "kepes"
+reconstruction       = "ocmp5"
+
+# mesh parameters
+#
+xblocks              =  2
+yblocks              =  1
+xmin                 = -1.28d+01
+xmax                 =  1.28d+01
+ymin                 = -6.40d+00
+ymax                 =  6.40d+00
+
+# refinement control
+#
+ncells               = 16
+maxlev               = 8
+refinement_variables = "vort jabs"
+currmin              = 5.0d-02
+currmax              = 2.0d-01
+vortmin              = 5.0d-02
+vortmax              = 2.0d-01
+
+# boundary conditions
+#
+xlbndry              = "periodic"
+xubndry              = "periodic"
+ylbndry              = "user"
+yubndry              = "user"
+
+# runtime control parameters
+#
+tmax                 = 2.0d+02
+cfl                  = 5.0d-01
+
+# data output control
+#
+snapshot_interval    = 1.0d+00
+compression_format   = "zstd"
+compression_level    = 19
+restart_number       = -1

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,31 @@ module user_problem
     real(kind=8), dimension(:)      , intent(in)    :: x, y, z
     real(kind=8), dimension(:,:,:,:), intent(inout) :: qn
 
+    integer :: i, j, jr
+
 !-------------------------------------------------------------------------------
 !
+    if (js == 1) then
+      do j = nbl, 1, -1
+        jr = nb + (nbl - j)
+        do i = il, iu
+          qn( : ,i,j,:) =   qn( : ,i,jr,:)
+
+          qn(ivy,i,j,:) = - qn(ivy,i,jr,:)
+          qn(iby,i,j,:) = - qn(iby,i,jr,:)
+        end do
+      end do
+    else
+      do j = neu, nn
+        jr = ne - (j - neu)
+        do i = il, iu
+          qn( : ,i,j,:) =   qn( : ,i,jr,:)
+
+          qn(ivy,i,j,:) = - qn(ivy,i,jr,:)
+          qn(iby,i,j,:) = - qn(iby,i,jr,:)
+        end do
+      end do
+    end if
 
 !-------------------------------------------------------------------------------
 !