Skip to content

Commit 9343be6

Browse files
hyeoksu-leehyeoksu-lee
committed
Fix a bug in s_superposition_instability_wave
1 parent e76c85d commit 9343be6

File tree

1 file changed

+45
-44
lines changed

1 file changed

+45
-44
lines changed

‎src/pre_process/m_initial_condition.fpp‎

Lines changed: 45 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -27,10 +27,10 @@ module m_initial_condition
2727
use m_patches
2828

2929
use m_assign_variables
30-
31-
use m_eigen_solver ! Subroutines to solve eigenvalue problem for
30+
31+
use m_eigen_solver ! Subroutines to solve eigenvalue problem for
3232
! complex general matrix
33-
33+
3434
! ==========================================================================
3535
! ==========================================================================
3636

@@ -298,37 +298,38 @@ contains
298298
subroutine s_superposition_instability_wave() ! ------------------------
299299
real(kind(0d0)), dimension(5,0:m,0:n,0:p) :: wave,wave1,wave2,wave_tmp
300300
real(kind(0d0)) :: tr,ti
301-
real(kind(0d0)) :: Ly
301+
real(kind(0d0)) :: Lx,Lz
302302
integer :: i,j,k
303303

304-
Ly = y_domain%end - y_domain%beg
304+
Lx = x_domain%end - x_domain%beg
305+
Lz = z_domain%end - z_domain%beg
305306

306307
wave = 0d0
307308
wave1 = 0d0
308309
wave2 = 0d0
309310

310311
! Compute 2D waves
311-
call s_instability_wave(2*pi*4.0/Ly,0d0,tr,ti,wave_tmp,0d0)
312+
call s_instability_wave(2*pi*4.0/Lx,0d0,tr,ti,wave_tmp,0d0)
312313
wave1 = wave1 + wave_tmp
313-
call s_instability_wave(2*pi*2.0/Ly,0d0,tr,ti,wave_tmp,0d0)
314+
call s_instability_wave(2*pi*2.0/Lx,0d0,tr,ti,wave_tmp,0d0)
314315
wave1 = wave1 + wave_tmp
315-
call s_instability_wave(2*pi*1.0/Ly,0d0,tr,ti,wave_tmp,0d0)
316+
call s_instability_wave(2*pi*1.0/Lx,0d0,tr,ti,wave_tmp,0d0)
316317
wave1 = wave1 + wave_tmp
317318
wave = wave1*0.05
318-
319+
319320
if (p > 0) then
320321
! Compute 3D waves with phase shifts.
321-
call s_instability_wave(2*pi*4.0/Ly, 2*pi*4.0/Ly,tr,ti,wave_tmp,2*pi*11d0/31d0)
322+
call s_instability_wave(2*pi*4.0/Lx, 2*pi*4.0/Lz,tr,ti,wave_tmp,2*pi*11d0/31d0)
322323
wave2 = wave2 + wave_tmp
323-
call s_instability_wave(2*pi*2.0/Ly, 2*pi*2.0/Ly,tr,ti,wave_tmp,2*pi*13d0/31d0)
324+
call s_instability_wave(2*pi*2.0/Lx, 2*pi*2.0/Lz,tr,ti,wave_tmp,2*pi*13d0/31d0)
324325
wave2 = wave2 + wave_tmp
325-
call s_instability_wave(2*pi*1.0/Ly, 2*pi*1.0/Ly,tr,ti,wave_tmp,2*pi*17d0/31d0)
326+
call s_instability_wave(2*pi*1.0/Lx, 2*pi*1.0/Lz,tr,ti,wave_tmp,2*pi*17d0/31d0)
326327
wave2 = wave2 + wave_tmp
327-
call s_instability_wave(2*pi*4.0/Ly,-2*pi*4.0/Ly,tr,ti,wave_tmp,2*pi*19d0/31d0)
328+
call s_instability_wave(2*pi*4.0/Lx,-2*pi*4.0/Lz,tr,ti,wave_tmp,2*pi*19d0/31d0)
328329
wave2 = wave2 + wave_tmp
329-
call s_instability_wave(2*pi*2.0/Ly,-2*pi*2.0/Ly,tr,ti,wave_tmp,2*pi*23d0/31d0)
330+
call s_instability_wave(2*pi*2.0/Lx,-2*pi*2.0/Lz,tr,ti,wave_tmp,2*pi*23d0/31d0)
330331
wave2 = wave2 + wave_tmp
331-
call s_instability_wave(2*pi*1.0/Ly,-2*pi*1.0/Ly,tr,ti,wave_tmp,2*pi*29d0/31d0)
332+
call s_instability_wave(2*pi*1.0/Lx,-2*pi*1.0/Lz,tr,ti,wave_tmp,2*pi*29d0/31d0)
332333
wave2 = wave2 + wave_tmp
333334
wave = wave + 0.15*wave2
334335
end if
@@ -346,7 +347,7 @@ contains
346347
end do
347348
end do
348349

349-
end subroutine s_superposition_instability_wave ! ----------------------
350+
end subroutine s_superposition_instability_wave ! ----------------------
350351

351352
!> This subroutine computes instability waves for a given set of spatial
352353
!! wavenumbers (alpha, beta) in x and z directions.
@@ -360,7 +361,7 @@ contains
360361
real(kind(0d0)),dimension(0:n,0:n) :: d !< differential operator in y dir
361362
real(kind(0d0)),dimension(0:5*(n+1)-1,0:5*(n+1)-1) :: ar,ai,br,bi,ci !< matrices for eigenvalue problem
362363
real(kind(0d0)),dimension(0:5*(n+1)-1,0:5*(n+1)-1) :: zr,zi !< eigenvectors
363-
real(kind(0d0)),dimension(0:5*(n+1)-1) :: wr,wi !< eigenvalues
364+
real(kind(0d0)),dimension(0:5*(n+1)-1) :: wr,wi !< eigenvalues
364365
real(kind(0d0)),dimension(0:5*(n+1)-1) :: fv1,fv2,fv3 !< temporary memory
365366
real(kind(0d0)) :: tr,ti !< most unstable eigenvalue
366367
real(kind(0d0)),dimension(0:5*(n+1)-1) :: vr,vi,vnr,vni !< most unstable eigenvector and normalized one
@@ -380,9 +381,9 @@ contains
380381

381382
! Assign mean profiles
382383
do j=0,n
383-
u_mean(j)=tanh(y_cc(j))
384-
t_mean(j)=1+0.5*(gam-1)*mach**2*(1-u_mean(j)**2)
385-
rho_mean(j)=1/T_mean(j)
384+
u_mean(j)=tanh(y_cc(j))
385+
t_mean(j)=1+0.5*(gam-1)*mach**2*(1-u_mean(j)**2)
386+
rho_mean(j)=1/T_mean(j)
386387
end do
387388

388389
! Compute differential operator in y-dir
@@ -393,24 +394,24 @@ contains
393394
d(1,0)=-1/(2*dy)
394395
d(1,2)= 1/(2*dy)
395396
do j=2,n-2
396-
d(j,j-2)= 1/(12*dy)
397-
d(j,j-1)=-8/(12*dy)
398-
d(j,j+1)= 8/(12*dy)
399-
d(j,j+2)=-1/(12*dy)
397+
d(j,j-2)= 1/(12*dy)
398+
d(j,j-1)=-8/(12*dy)
399+
d(j,j+1)= 8/(12*dy)
400+
d(j,j+2)=-1/(12*dy)
400401
end do
401402
d(n-1,n-2)=-1/(2*dy)
402403
d(n-1,n) = 1/(2*dy)
403404

404405
! Compute y-derivatives of rho, u, T
405406
do j=0,n
406-
drho_mean(j)=0
407-
du_mean(j)=0
408-
dt_mean(j)=0
409-
do k=0,n
410-
drho_mean(j) = drho_mean(j)+d(j,k)*rho_mean(k)
411-
du_mean(j) = du_mean(j)+d(j,k)*u_mean(k)
412-
dt_mean(j) = dt_mean(j)+d(j,k)*t_mean(k)
413-
end do
407+
drho_mean(j)=0
408+
du_mean(j)=0
409+
dt_mean(j)=0
410+
do k=0,n
411+
drho_mean(j) = drho_mean(j)+d(j,k)*rho_mean(k)
412+
du_mean(j) = du_mean(j)+d(j,k)*u_mean(k)
413+
dt_mean(j) = dt_mean(j)+d(j,k)*t_mean(k)
414+
end do
414415
end do
415416

416417
! Compute B and C, then A = B + C
@@ -455,7 +456,7 @@ contains
455456

456457
! Compute eigenvalues and eigenvectors
457458
call cg(5*(n+1),5*(n+1),ar,ai,wr,wi,zr,zi,fv1,fv2,fv3,ierr)
458-
459+
459460
! Generate instability wave
460461
call s_generate_wave(5*(n+1),wr,wi,zr,zi,alpha,beta,wave,shift)
461462

@@ -485,9 +486,9 @@ contains
485486
end do
486487
vr = zr(:,k)
487488
vi = zi(:,k)
488-
489+
489490
! Normalize the eigenvector by its component with the largest modulus.
490-
norm = 0d0
491+
norm = 0d0
491492
do i=0,nl-1
492493
if (dsqrt(vr(i)**2+vi(i)**2) .gt. norm) then
493494
idx = i
@@ -502,7 +503,7 @@ contains
502503
vnr(i) = cr
503504
vni(i) = ci
504505
end do
505-
506+
506507
! Generate an instability wave
507508
do i=0,m
508509
do j=0,n
@@ -512,17 +513,17 @@ contains
512513
else
513514
ang = alpha*x_cc(i)+beta*z_cc(k)+shift
514515
end if
515-
wave(1,i,j,k) = vnr(j)*cos(ang)-vni(j)*sin(ang) ! rho
516-
wave(2,i,j,k) = vnr((n+1)+j)*cos(ang)-vni((n+1)+j)*sin(ang) ! u
517-
wave(3,i,j,k) = vnr(2*(n+1)+j)*cos(ang)-vni(2*(n+1)+j)*sin(ang) ! v
518-
wave(4,i,j,k) = vnr(3*(n+1)+j)*cos(ang)-vni(3*(n+1)+j)*sin(ang) ! w
519-
wave(5,i,j,k) = vnr(4*(n+1)+j)*cos(ang)-vni(4*(n+1)+j)*sin(ang) ! T
516+
wave(1,i,j,k) = vnr(j)*cos(ang)-vni(j)*sin(ang) ! rho
517+
wave(2,i,j,k) = vnr((n+1)+j)*cos(ang)-vni((n+1)+j)*sin(ang) ! u
518+
wave(3,i,j,k) = vnr(2*(n+1)+j)*cos(ang)-vni(2*(n+1)+j)*sin(ang) ! v
519+
wave(4,i,j,k) = vnr(3*(n+1)+j)*cos(ang)-vni(3*(n+1)+j)*sin(ang) ! w
520+
wave(5,i,j,k) = vnr(4*(n+1)+j)*cos(ang)-vni(4*(n+1)+j)*sin(ang) ! T
520521
end do
521522
end do
522523
end do
523-
524-
end subroutine s_generate_wave
525-
524+
525+
end subroutine s_generate_wave
526+
526527
!> Deallocation procedures for the module
527528
subroutine s_finalize_initial_condition_module() ! ---------------------
528529

0 commit comments

Comments
 (0)