diff --git a/sources/user_problem.F90 b/sources/user_problem.F90 index 14540e0..a15ff2b 100644 --- a/sources/user_problem.F90 +++ b/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 !