From 6599e4e54ff8227f8496b59b442704d3c23020d3 Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Tue, 28 Mar 2017 23:23:50 -0300
Subject: [PATCH] RECONNECTION: Implement one more type of perturbation.

This perturbation is similar to the magnetic field one, but
applied to the velocity field.

Also allow the user to choose the perturbation type through the
parameter 'perturbation'.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
 src/user_problem.F90 | 73 +++++++++++++++++++++++++++++++-------------
 1 file changed, 52 insertions(+), 21 deletions(-)

diff --git a/src/user_problem.F90 b/src/user_problem.F90
index 218c9a4..e7b292b 100644
--- a/src/user_problem.F90
+++ b/src/user_problem.F90
@@ -53,6 +53,7 @@ module user_problem
   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 :: kper = 1.00d+00
   real(kind=8), save :: xcut = 4.00d-01
   real(kind=8), save :: ycut = 1.00d-02
   real(kind=8), save :: yth  = 1.00d-16
@@ -60,6 +61,7 @@ module user_problem
   real(kind=8), save :: pmag = 5.00d-01
   real(kind=8), save :: blim = 1.00d+00
   real(kind=8), save :: zeta = 0.00d+00
+  integer     , save :: pert = 1
 
 ! flag indicating if the gravitational source term is enabled
 !
@@ -104,6 +106,7 @@ module user_problem
 ! local variables
 !
     character(len=64) :: problem_name   = "none"
+    character(len=64) :: perturbation   = "noise"
 !
 !-------------------------------------------------------------------------------
 !
@@ -124,7 +127,19 @@ module user_problem
 
 ! get the problem name
 !
-    call get_parameter_string("problem", problem_name)
+    call get_parameter_string("problem"     , problem_name)
+    call get_parameter_string("perturbation", perturbation)
+
+! choose the perturbation type
+!
+    select case(perturbation)
+    case('noise', 'random')
+      pert = 1
+    case('vel', 'velocity_mode')
+      pert = 2
+    case('mag', 'magnetic_mode')
+      pert = 3
+    end select
 
 ! get the reconnection problem parameters
 !
@@ -134,6 +149,7 @@ module user_problem
     call get_parameter_real("bper"  , bper)
     call get_parameter_real("bgui"  , bgui)
     call get_parameter_real("vper"  , vper)
+    call get_parameter_real("kper"  , kper)
     call get_parameter_real("xcut"  , xcut)
     call get_parameter_real("ycut"  , ycut)
     call get_parameter_real("yth"   , yth )
@@ -238,7 +254,7 @@ module user_problem
     use equations  , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
     use equations  , only : csnd2
     use operators  , only : curl
-    use random     , only : randomn
+    use random     , only : randomn, randomu
 
 ! local variables are not implicit by default
 !
@@ -251,7 +267,7 @@ module user_problem
 ! local variables
 !
     integer      :: i, j, k
-    real(kind=8) :: yt, yp, yl, yu
+    real(kind=8) :: xp, yp, yt, yl, yu
     real(kind=8) :: yrat, itanh
     real(kind=8) :: vx, vy, vz, vv
 
@@ -289,7 +305,7 @@ module user_problem
 
 ! calculate the perturbation of magnetic field
 !
-    if (bper /= 0.0d+00) then
+    if (bper /= 0.0d+00 .and. pert == 3) then
 
 ! initiate the vector potential (we use velocity components to store vector
 ! potential temporarily, and we store magnetic field perturbation in U)
@@ -350,7 +366,7 @@ module user_problem
 
 ! add the magnetic field perturbation
 !
-          if (bper /= 0.0d+00) then
+          if (bper /= 0.0d+00 .and. pert == 3) 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)
@@ -384,34 +400,49 @@ module user_problem
 
 ! set the random velocity field near the current sheet
 !
-        if (abs(x(i)) <= xcut) then
-          do j = 1, jm
-            if (abs(y(j)) <= ycut) then
+        if (pert == 1) then
+          if (abs(x(i)) <= xcut) then
+            do j = 1, jm
+              if (abs(y(j)) <= ycut) then
 
-              vv = 0.0d+00
-              do while(vv == 0.0d+00)
-                vx = randomn()
-                vy = randomn()
+                vv = 0.0d+00
+                do while(vv == 0.0d+00)
+                  vx = randomn()
+                  vy = randomn()
 #if NDIMS == 3
-                vz = randomn()
+                  vz = randomn()
 #endif /* NDIMS == 3 */
 
 #if NDIMS == 3
-                vv = sqrt(vx * vx + vy * vy + vz * vz)
+                  vv = sqrt(vx * vx + vy * vy + vz * vz)
 #else /* NDIMS == 3 */
-                vv = sqrt(vx * vx + vy * vy)
+                  vv = sqrt(vx * vx + vy * vy)
 #endif /* NDIMS == 3 */
-              end do
+                end do
 
-              q(ivx,j) = vper * (vx / vv)
-              q(ivy,j) = vper * (vy / vv)
+                q(ivx,j) = vper * (vx / vv)
+                q(ivy,j) = vper * (vy / vv)
 #if NDIMS == 3
-              q(ivz,j) = vper * (vz / vv)
+                q(ivz,j) = vper * (vz / vv)
 #endif /* NDIMS == 3 */
 
-            end if ! |y| < ycut
+              end if ! |y| < ycut
+            end do ! j = 1, jm
+          end if ! |x| < xcut
+        end if ! pert == 1
+
+        if (pert == 2) then
+          xp = pi2 * kper * (x(i) + 0.5d+00)
+          do j = 1, jm
+            yp = y(j) / pth
+            vx = (2.0d+00 * yp * tanh(yp) - 1.0d+00 / cosh(yp)**2)             &
+               * exp(- yp**2) * sin(xp) / (pi2 * kper * pth)
+            vy = tanh(yp) * exp(- yp**2) * cos(xp)
+
+            q(ivx,j) = vper * vx
+            q(ivy,j) = vper * vy
           end do ! j = 1, jm
-        end if ! |x| < xcut
+        end if ! pert == 2
 
 ! convert the primitive variables to the conservative ones
 !