USER_PROBLEM: Add new perturbation of spectrum of waves.

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

sources/user_problem.F90

@@ -52,7 +52,7 @@ module user_problem
   real(kind=8), save :: tdec  = 1.00d+00
 
   integer     , save :: pert  = 0
-  integer     , save :: nper  = 10
+  integer     , save :: nper  = 16
   real(kind=8), save :: bper  = 0.00d+00
   real(kind=8), save :: vper  = 0.00d+00
   real(kind=8), save :: kper  = 1.00d+00
@@ -264,8 +264,10 @@ module user_problem
       pert = 1
     case('multi-wave', 'random waves', 'random-waves', 'random_waves')
       pert = 2
-    case('external')
+    case('wave-spectrum', 'wave spectrum', 'wave_spectrum')
       pert = 3
+    case('external')
+      pert = 4
     case default
       perturbation = 'mode'
       pert = 1
@@ -336,6 +338,27 @@ module user_problem
       end if
     end if ! pert = 2
 
+! prepare the wave vectors for wave-spectrum perturbation
+!
+    if (pert == 3) then
+
+! allocate phase and wave vector components
+!
+      allocate(kx(nper), ph(nper), stat = status)
+
+      if (status == 0) then
+        do n = 1, nper
+          kx(n) = pi2 * n / xlen
+          ph(n) = pi2 * randuni()
+        end do
+      else
+        if (verbose) &
+          call print_message(loc, &
+                           "Could not allocate space for perturbation vectors!")
+        return
+      end if
+    end if ! pert = 3
+
 ! determine if to append or create another file
 !
     call get_parameter("reconnection_append", append)
@@ -493,6 +516,12 @@ module user_problem
         call print_message(loc, &
                          "Could not deallocate space for perturbation vectors!")
     end if
+    if (pert == 3) then
+      deallocate(kx, ph, stat = status)
+      if (status /= 0) &
+        call print_message(loc, &
+                         "Could not deallocate space for perturbation vectors!")
+    end if
 
 !-------------------------------------------------------------------------------
 !
@@ -918,10 +947,56 @@ module user_problem
 
     end if ! pert == 2
 
-! load the external perturbation
+! prepare the spectrum of waves
 !
     if (pert == 3) then
 
+! of velocity
+!
+      if (abs(vper) > 0.0d+00) then
+
+        pot(:,:,:,:) = 0.0d+00
+
+        do n = 1, nper
+          do j = 1, nn
+            yp = sqrt(max(0.0d+00, 1.0d+00 - tanh(yc(j) / dlta)**2))
+            do i = 1, nn
+              pot(3,i,j,:) = pot(3,i,j,:) &
+                           + vper / kx(n) * cos(kx(n) * xc(i) + ph(n)) * yp
+            end do ! i = 1, nn
+          end do ! j = 1, nn
+        end do ! n = 1, nper
+
+        call curl(dh(1:3), pot(:,:,:,:), qpert(ivx:ivz,:,:,:))
+
+      end if ! vper /= 0.0
+
+! of magnetic field
+!
+      if (abs(bper) > 0.0d+00) then
+
+        pot(:,:,:,:) = 0.0d+00
+
+        do n = 1, nper
+          do j = 1, nn
+            yp = sqrt(max(0.0d+00, 1.0d+00 - tanh(yc(j) / dlta)**2))
+            do i = 1, nn
+              pot(3,i,j,:) = pot(3,i,j,:) &
+                           + bper / kx(n) * cos(kx(n) * xc(i) + ph(n)) * yp
+            end do ! i = 1, nn
+          end do ! j = 1, nn
+        end do ! n = 1, nper
+
+        call curl(dh(1:3), pot(:,:,:,:), qpert(ibx:ibz,:,:,:))
+
+      end if ! bper /= 0.0
+
+    end if ! pert == 3
+
+! load the external perturbation
+!
+    if (pert == 4) then
+
       write(filename, "('perturbations-',i2.2,'-',i2.2,'.dat')") &
                                                             nc, pdata%meta%level
       open(newunit=io, file=filename)
@@ -962,7 +1037,7 @@ module user_problem
 
       deallocate(yt, qr, qi)
 
-    end if ! pert == 3
+    end if ! pert == 4
 
 ! iterate over all positions in the XZ plane
 !