diff --git a/sources/user_problem.F90 b/sources/user_problem.F90 index c37c089..98e4b4b 100644 --- a/sources/user_problem.F90 +++ b/sources/user_problem.F90 @@ -44,8 +44,22 @@ module user_problem 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 :: prdl = 0.00d+00 + real(kind=8), save :: nu = 0.00d+00 real(kind=8), save :: eta = 0.00d+00 + integer , save :: pert = -1 + integer , save :: nper = 0 + real(kind=8), save :: vper = 0.00d+00 + real(kind=8), save :: kper = 1.00d+00 + real(kind=8), save :: kdel = 1.00d+00 + real(kind=8), save :: xcut = 1.00d+99 + real(kind=8), save :: ycut = 1.00d+99 + real(kind=8), save :: xdec = 1.00d+00 + real(kind=8), save :: ydec = 1.00d+00 + + real(kind=8), dimension(:), allocatable :: kx, ky, kz, ux, uy, uz, ph + integer(kind=4), save :: runit = 33 logical, save :: resistive = .false. @@ -78,10 +92,13 @@ module user_problem subroutine initialize_user_problem(problem, rcount, verbose, status) use boundaries , only : custom_boundary_x, custom_boundary_y + use constants , only : pi, pi2 + use coordinates, only : xlen, zlen 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 + use random , only : reset_seeds, randuni, randsym implicit none @@ -90,15 +107,19 @@ module user_problem logical , intent(in) :: verbose integer , intent(out) :: status - character(len=64) :: tmp, append = "off", fname - logical :: flag + character(len=64) :: str_param, append = "off", fname + logical :: flag, perturb_tearing + integer :: n, kd, ki, kj, kk, kv, ks + real(kind=8) :: kl2, ku2, kv2, ka + real(kind=8) :: thv, thh, tx, ty, tz character(len=*), parameter :: & loc = 'USER_PROBLEM::initialize_user_problem()' !------------------------------------------------------------------------------- ! - call get_parameter("profile", tmp) - select case(tmp) + str_param = "sincos" + call get_parameter("profile", str_param) + select case(str_param) case('bhattacharjee', 'sincos', 'trigonometric') profile = 2 case default @@ -128,6 +149,7 @@ module user_problem status = 1 return end if + call get_parameter("viscosity" , nu ) call get_parameter("resistivity", eta) if (eta < 0.0d+00) then if (verbose) & @@ -146,6 +168,8 @@ module user_problem status = 1 return end if + ycut = 4 * dlta + ydec = dlta pmag = 0.5d+00 * (bamp**2 + bgui**2) pres = beta * pmag @@ -158,6 +182,142 @@ module user_problem csnd = sqrt(csnd2) valf = bamp / sqrt(dens) lund = valf / max(tiny(eta), eta) + prdl = nu / max(tiny(eta), eta) + +! get the perturbation parameters +! + call get_parameter("vper", vper) + call get_parameter("kper", kper) + call get_parameter("kdel", kdel) + call get_parameter("xcut", xcut) + call get_parameter("ycut", ycut) + call get_parameter("xdec", xdec) + call get_parameter("ydec", ydec) + +! choose the perturbation type +! + str_param = "none" + call get_parameter("perturbation", str_param) + if (vper == 0.0d+00) str_param = "none" + select case(str_param) + case('noise', 'random') + pert = 0 + call reset_seeds('random') + case('mode', 'one mode', 'one-mode', 'one_mode') + pert = 1 + case('multi-wave', 'random waves', 'random-waves', 'random_waves') + pert = 2 + case default + pert = -1 + end select + +! enable or disable tearing modes perturbation +! + str_param = "on" + call get_parameter("perturb_tearing", str_param) + select case(str_param) + case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") + perturb_tearing = .true. + case default + perturb_tearing = .false. + end select + +! prepare the wave vectors for multi-wave perturbation +! + if (pert == 2) then + + kl2 = max(0.0d+00, kper - kdel)**2 + ku2 = (kper + kdel)**2 + + kv = int(kper) + ks = int(ceiling(1.0d+00 / (pi2 * dlta))) + if (perturb_tearing) ks = 0 + n = 0 +#if NDIMS == 3 + kd = int(xlen / zlen) + do kk = - kv + mod(kv, kd), kv - mod(kv, kd), kd +#else /* NDIMS == 3 */ + kk = 0 +#endif /* NDIMS == 3 */ + do kj = - kv, kv + do ki = ks, kv + kv2 = 1.0d+00 * (ki**2 + kj**2 + kk**2) + if (kl2 <= kv2 .and. kv2 <= ku2) then + if (.not. (ki == 0 .and. kj == 0 .and. kk < 0) .and. & + .not. (ki == 0 .and. kj < 0)) then + n = n + 1 + end if + end if + end do + end do +#if NDIMS == 3 + end do +#endif /* NDIMS == 3 */ + if (n > 0) nper = n + +! allocate phase and wave vector components +! + allocate(kx(nper), ky(nper), kz(nper), ux(nper), uy(nper), uz(nper), & + ph(nper), stat = status) + + if (status == 0) then + +! choose random wave vector directions +! + n = 0 +#if NDIMS == 3 + do kk = - kv + mod(kv, kd), kv - mod(kv, kd), kd +#else /* NDIMS == 3 */ + kk = 0 +#endif /* NDIMS == 3 */ + do kj = - kv, kv + do ki = ks, kv + kv2 = 1.0d+00 * (ki**2 + kj**2 + kk**2) + if (kl2 <= kv2 .and. kv2 <= ku2) then + if (.not. (ki == 0 .and. kj == 0 .and. kk < 0) .and. & + .not. (ki == 0 .and. kj < 0)) then + n = n + 1 + kx(n) = pi2 * ki + ky(n) = pi2 * kj + kz(n) = pi2 * kk + if (kk == 0) then + ka = sqrt(kx(n)**2 + ky(n)**2) + ux(n) = ky(n) / ka + uy(n) = - kx(n) / ka + uz(n) = 0.0d+00 + else + ka = 0.0d+00 + do while(ka < 1.0d-08) + thh = pi2 * randuni() + thv = pi * randsym() + tx = cos(thh) * cos(thv) + ty = sin(thh) * cos(thv) + tz = sin(thv) + ux(n) = ty * kz(n) - tz * ky(n) + uy(n) = tz * kx(n) - tx * kz(n) + uz(n) = tx * ky(n) - ty * kx(n) + ka = sqrt(ux(n)**2 + uy(n)**2 + uz(n)**2) + end do + ux(n) = ux(n) / ka + uy(n) = uy(n) / ka + uz(n) = uz(n) / ka + end if + ph(n) = pi2 * randuni() + end if + end if + end do + end do +#if NDIMS == 3 + end do +#endif /* NDIMS == 3 */ + + else + if (verbose) & + call print_message(loc, & + "Could not allocate space for perturbation vectors!") + return + end if + end if ! pert = 2 ! determine if to append or create another file ! @@ -215,13 +375,37 @@ module user_problem 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) + call print_parameter(verbose, '( ) Pr (Prandtl number)' , prdl) + call print_parameter(verbose, '(*) delta (thickness)' , dlta) + select case(pert) + case(0) + call print_parameter(verbose, '(*) perturbation', 'noise') + case(1) + call print_parameter(verbose, '(*) perturbation', 'one mode') + case(2) + call print_parameter(verbose, '(*) perturbation', 'waves at given k') + end select + if (pert == 2) then + if (perturb_tearing) then + call print_parameter(verbose, '(*) perturb tearing modes', 'on') + else + call print_parameter(verbose, '(*) perturb tearing modes', 'off') + end if + end if + call print_parameter(verbose, '(*) vper' , vper) + if (pert == 2) then + call print_parameter(verbose, '(*) kper' , kper) + call print_parameter(verbose, '(*) kdel' , kdel) + call print_parameter(verbose, '(*) nper' , nper) + call print_parameter(verbose, '(*) ycut' , ycut) + call print_parameter(verbose, '(*) ydec' , ydec) + end if setup_problem => setup_user_problem custom_boundary_x => user_boundary_x @@ -248,14 +432,25 @@ module user_problem ! subroutine finalize_user_problem(status) + use helpers, only : print_message + implicit none integer, intent(out) :: status + character(len=*), parameter :: loc = 'USER_PROBLEM::finalize_user_problem()' + !------------------------------------------------------------------------------- ! status = 0 + if (pert == 2) then + deallocate(kx, ky, kz, ux, uy, uz, ph, stat = status) + if (status /= 0) & + call print_message(loc, & + "Could not deallocate space for perturbation vectors!") + end if + !------------------------------------------------------------------------------- ! end subroutine finalize_user_problem @@ -279,25 +474,46 @@ module user_problem use blocks , only : block_data use constants , only : pi, pi2 use coordinates, only : nn => bcells, xmin, xmax +#if NDIMS == 3 + use coordinates, only : ax, ay, az +#else /* NDIMS == 3 */ use coordinates, only : ax, ay +#endif /* NDIMS == 3 */ use equations , only : prim2cons, csnd2 use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp + use random , only : randnorz implicit none type(block_data), pointer, intent(inout) :: pdata - integer :: i, j, k - real(kind=8) :: r1, r2, bx, by, bz + logical :: flag + integer :: i, j, k, n + real(kind=8) :: r1, r2, qx, qy, qz, fv, xp, yp, kv, va + complex(kind=8) :: c + real(kind=8), dimension(3) :: q + real(kind=8), dimension(6) :: t +#if NDIMS == 3 + real(kind=8), dimension(nn) :: x, y, z +#else /* NDIMS == 3 */ real(kind=8), dimension(nn) :: x, y +#endif /* NDIMS == 3 */ + real(kind=8), dimension(nn) :: fx, fy real(kind=8), dimension(nn) :: snx, csx, sny, csy, thy, shy - real(kind=8), dimension(nn,nn) :: az + real(kind=8), dimension(nn,nn) :: ps !------------------------------------------------------------------------------- ! + k = 1 + x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) +#if NDIMS == 3 + z(:) = pdata%meta%zmin + az(pdata%meta%level,:) +#endif /* NDIMS == 3 */ +! choose profile +! select case(profile) case(2) @@ -313,13 +529,13 @@ module user_problem do j = 1, nn do i = 1, nn - bx = bamp * csx(i) * (csy(j) * thy(j) + sny(j) * shy(j) / pi2) - by = bamp * 5.0d-01 * snx(i) * sny(j) * thy(j) - bz = sqrt(max(0.0d+00, bgui**2 + zeta * (bamp**2 - bx**2 - by**2))) + qx = bamp * csx(i) * (csy(j) * thy(j) + sny(j) * shy(j) / pi2) + qy = bamp * 5.0d-01 * snx(i) * sny(j) * thy(j) + qz = sqrt(max(0.0d+00, bgui**2 + zeta * (bamp**2 - qx**2 - qy**2))) - pdata%q(ibx,i,j,:) = bx - pdata%q(iby,i,j,:) = by - pdata%q(ibz,i,j,:) = bz + pdata%q(ibx,i,j,:) = qx + pdata%q(iby,i,j,:) = qy + pdata%q(ibz,i,j,:) = qz end do end do pdata%q(ibp,:,:,:) = 0.0d+00 @@ -329,20 +545,20 @@ module user_problem do i = 1, nn r1 = 5.0d-01 - sqrt(x(i)**2 + ((2.0d+00 * y(j) - 5.0d-01))**2) r2 = 5.0d-01 - sqrt(x(i)**2 + ((2.0d+00 * y(j) + 5.0d-01))**2) - az(i,j) = 2.5d-01 * (r1 * (tanh(r1 / dlta) + 1.0d+00) & + ps(i,j) = 2.5d-01 * (r1 * (tanh(r1 / dlta) + 1.0d+00) & + r2 * (tanh(r2 / dlta) + 1.0d+00)) end do end do do j = 2, nn - 1 do i = 2, nn - 1 - bx = (az(i,j+1) - az(i,j-1)) / (y(j+1) - y(j-1)) - by = (az(i-1,j) - az(i+1,j)) / (x(i+1) - x(i-1)) - bz = sqrt(max(0.0d+00, bgui**2 + zeta * (bamp**2 - bx**2 - by**2))) + qx = (ps(i,j+1) - ps(i,j-1)) / (y(j+1) - y(j-1)) + qy = (ps(i-1,j) - ps(i+1,j)) / (x(i+1) - x(i-1)) + qz = sqrt(max(0.0d+00, bgui**2 + zeta * (bamp**2 - qx**2 - qy**2))) - pdata%q(ibx,i,j,:) = bx - pdata%q(iby,i,j,:) = by - pdata%q(ibz,i,j,:) = bz + pdata%q(ibx,i,j,:) = qx + pdata%q(iby,i,j,:) = qy + pdata%q(ibz,i,j,:) = qz end do end do pdata%q(ibp,:,:,:) = 0.0d+00 @@ -361,6 +577,108 @@ module user_problem pdata%q(ivy,:,:,:) = 0.0d+00 pdata%q(ivz,:,:,:) = 0.0d+00 +! prepare decaying factors +! + if (pert >= 0) then + + fv = 5.0d-01 * pi + + do i = 1, nn + xp = fv * min(1.0d+00, max(0.0d+00, abs(x(i)) - xcut) / xdec) + fx(i) = cos(xp)**2 + end do ! i = 1, nn + + do j = 1, nn + yp = fv * min(1.0d+00, max(0.0d+00, abs(y(j)) - ycut) / ydec) + fy(j) = cos(yp)**2 + end do ! i = 1, nn + + end if ! pert >= 0 + +! add some noise to the velocity components +! + if (pert == 0) then + + fv = vper / sqrt(2.0d+00) + + flag = .true. +#if NDIMS == 3 + do k = 1, nn +#endif /* NDIMS == 3 */ + do j = 1, nn + if (fy(j) > 0.0d+00) then + do i = 1, nn + if (fx(i) > 0.0d+00) then + + va = fv * fx(i) * fy(j) + + if (flag) then + c = randnorz() + t(1:2) = [ real(c), aimag(c) ] + c = randnorz() + t(3:4) = [ real(c), aimag(c) ] + c = randnorz() + t(5:6) = [ real(c), aimag(c) ] + + q(1:3) = t(1:3) + else + q(1:3) = t(4:6) + end if + flag = .not. flag + + pdata%q(ivx,i,j,k) = va * q(1) + pdata%q(ivy,i,j,k) = va * q(2) + pdata%q(ivz,i,j,k) = va * q(3) + + end if ! fx > 0.0 + end do ! i = 1, nn + end if ! fy > 0.0 + end do ! j = 1, nn +#if NDIMS == 3 + end do ! k = 1, nn +#endif /* NDIMS == 3 */ + + end if ! pert == 0 + +! add the multi-wave perturbation of velocity +! + if (pert == 2) then + +#if NDIMS == 3 + do k = 1, nn +#endif /* NDIMS == 3 */ + do j = 1, nn + if (fy(j) > 0.0d+00) then + do i = 1, nn + if (fx(i) > 0.0d+00) then + + fv = vper * fx(i) * fy(j) + + do n = 1, nper +#if NDIMS == 3 + kv = kx(n) * x(i) + ky(n) * y(j) + kz(n) * z(k) + ph(n) +#else /* NDIMS == 3 */ + kv = kx(n) * x(i) + ky(n) * y(j) + ph(n) +#endif /* NDIMS == 3 */ + va = fv * sin(kv) + + pdata%q(ivx,i,j,k) = pdata%q(ivx,i,j,k) + va * ux(n) + pdata%q(ivy,i,j,k) = pdata%q(ivy,i,j,k) + va * uy(n) +#if NDIMS == 3 + pdata%q(ivz,i,j,k) = pdata%q(ivz,i,j,k) + va * uz(n) +#endif /* NDIMS == 3 */ + end do + + end if ! fx > 0.0 + end do ! i = 1, nn + end if ! fy > 0.0 + end do ! j = 1, nn +#if NDIMS == 3 + end do ! k = 1, nn +#endif /* NDIMS == 3 */ + + end if ! pert == 2 + #if NDIMS == 3 do k = 1, nn #else /* NDIMS == 3 */