Skip to content

Commit d49f636

Browse files
sbryngelsonsbryngelson
authored
Merge pull request #140 from lee-hyeoksu/master
2 parents 61b0ec3 + a54494a commit d49f636

File tree

1 file changed

+48
-45
lines changed

1 file changed

+48
-45
lines changed

‎src/pre_process/m_initial_condition.fpp‎

Lines changed: 48 additions & 45 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,40 @@ 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
305-
304+
Lx = x_domain%end - x_domain%beg
305+
if (p > 0) then
306+
Lz = z_domain%end - z_domain%beg
307+
end if
308+
306309
wave = 0d0
307310
wave1 = 0d0
308311
wave2 = 0d0
309312

310313
! Compute 2D waves
311-
call s_instability_wave(2*pi*4.0/Ly,0d0,tr,ti,wave_tmp,0d0)
314+
call s_instability_wave(2*pi*4.0/Lx,0d0,tr,ti,wave_tmp,0d0)
312315
wave1 = wave1 + wave_tmp
313-
call s_instability_wave(2*pi*2.0/Ly,0d0,tr,ti,wave_tmp,0d0)
316+
call s_instability_wave(2*pi*2.0/Lx,0d0,tr,ti,wave_tmp,0d0)
314317
wave1 = wave1 + wave_tmp
315-
call s_instability_wave(2*pi*1.0/Ly,0d0,tr,ti,wave_tmp,0d0)
318+
call s_instability_wave(2*pi*1.0/Lx,0d0,tr,ti,wave_tmp,0d0)
316319
wave1 = wave1 + wave_tmp
317320
wave = wave1*0.05
318-
321+
319322
if (p > 0) then
320323
! 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)
324+
call s_instability_wave(2*pi*4.0/Lx, 2*pi*4.0/Lz,tr,ti,wave_tmp,2*pi*11d0/31d0)
322325
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)
326+
call s_instability_wave(2*pi*2.0/Lx, 2*pi*2.0/Lz,tr,ti,wave_tmp,2*pi*13d0/31d0)
324327
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)
328+
call s_instability_wave(2*pi*1.0/Lx, 2*pi*1.0/Lz,tr,ti,wave_tmp,2*pi*17d0/31d0)
326329
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)
330+
call s_instability_wave(2*pi*4.0/Lx,-2*pi*4.0/Lz,tr,ti,wave_tmp,2*pi*19d0/31d0)
328331
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)
332+
call s_instability_wave(2*pi*2.0/Lx,-2*pi*2.0/Lz,tr,ti,wave_tmp,2*pi*23d0/31d0)
330333
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)
334+
call s_instability_wave(2*pi*1.0/Lx,-2*pi*1.0/Lz,tr,ti,wave_tmp,2*pi*29d0/31d0)
332335
wave2 = wave2 + wave_tmp
333336
wave = wave + 0.15*wave2
334337
end if
@@ -346,7 +349,7 @@ contains
346349
end do
347350
end do
348351

349-
end subroutine s_superposition_instability_wave ! ----------------------
352+
end subroutine s_superposition_instability_wave ! ----------------------
350353

351354
!> This subroutine computes instability waves for a given set of spatial
352355
!! wavenumbers (alpha, beta) in x and z directions.
@@ -360,7 +363,7 @@ contains
360363
real(kind(0d0)),dimension(0:n,0:n) :: d !< differential operator in y dir
361364
real(kind(0d0)),dimension(0:5*(n+1)-1,0:5*(n+1)-1) :: ar,ai,br,bi,ci !< matrices for eigenvalue problem
362365
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
366+
real(kind(0d0)),dimension(0:5*(n+1)-1) :: wr,wi !< eigenvalues
364367
real(kind(0d0)),dimension(0:5*(n+1)-1) :: fv1,fv2,fv3 !< temporary memory
365368
real(kind(0d0)) :: tr,ti !< most unstable eigenvalue
366369
real(kind(0d0)),dimension(0:5*(n+1)-1) :: vr,vi,vnr,vni !< most unstable eigenvector and normalized one
@@ -380,9 +383,9 @@ contains
380383

381384
! Assign mean profiles
382385
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)
386+
u_mean(j)=tanh(y_cc(j))
387+
t_mean(j)=1+0.5*(gam-1)*mach**2*(1-u_mean(j)**2)
388+
rho_mean(j)=1/T_mean(j)
386389
end do
387390

388391
! Compute differential operator in y-dir
@@ -393,24 +396,24 @@ contains
393396
d(1,0)=-1/(2*dy)
394397
d(1,2)= 1/(2*dy)
395398
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)
399+
d(j,j-2)= 1/(12*dy)
400+
d(j,j-1)=-8/(12*dy)
401+
d(j,j+1)= 8/(12*dy)
402+
d(j,j+2)=-1/(12*dy)
400403
end do
401404
d(n-1,n-2)=-1/(2*dy)
402405
d(n-1,n) = 1/(2*dy)
403406

404407
! Compute y-derivatives of rho, u, T
405408
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
409+
drho_mean(j)=0
410+
du_mean(j)=0
411+
dt_mean(j)=0
412+
do k=0,n
413+
drho_mean(j) = drho_mean(j)+d(j,k)*rho_mean(k)
414+
du_mean(j) = du_mean(j)+d(j,k)*u_mean(k)
415+
dt_mean(j) = dt_mean(j)+d(j,k)*t_mean(k)
416+
end do
414417
end do
415418

416419
! Compute B and C, then A = B + C
@@ -455,7 +458,7 @@ contains
455458

456459
! Compute eigenvalues and eigenvectors
457460
call cg(5*(n+1),5*(n+1),ar,ai,wr,wi,zr,zi,fv1,fv2,fv3,ierr)
458-
461+
459462
! Generate instability wave
460463
call s_generate_wave(5*(n+1),wr,wi,zr,zi,alpha,beta,wave,shift)
461464

@@ -485,9 +488,9 @@ contains
485488
end do
486489
vr = zr(:,k)
487490
vi = zi(:,k)
488-
491+
489492
! Normalize the eigenvector by its component with the largest modulus.
490-
norm = 0d0
493+
norm = 0d0
491494
do i=0,nl-1
492495
if (dsqrt(vr(i)**2+vi(i)**2) .gt. norm) then
493496
idx = i
@@ -502,7 +505,7 @@ contains
502505
vnr(i) = cr
503506
vni(i) = ci
504507
end do
505-
508+
506509
! Generate an instability wave
507510
do i=0,m
508511
do j=0,n
@@ -512,17 +515,17 @@ contains
512515
else
513516
ang = alpha*x_cc(i)+beta*z_cc(k)+shift
514517
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
518+
wave(1,i,j,k) = vnr(j)*cos(ang)-vni(j)*sin(ang) ! rho
519+
wave(2,i,j,k) = vnr((n+1)+j)*cos(ang)-vni((n+1)+j)*sin(ang) ! u
520+
wave(3,i,j,k) = vnr(2*(n+1)+j)*cos(ang)-vni(2*(n+1)+j)*sin(ang) ! v
521+
wave(4,i,j,k) = vnr(3*(n+1)+j)*cos(ang)-vni(3*(n+1)+j)*sin(ang) ! w
522+
wave(5,i,j,k) = vnr(4*(n+1)+j)*cos(ang)-vni(4*(n+1)+j)*sin(ang) ! T
520523
end do
521524
end do
522525
end do
523-
524-
end subroutine s_generate_wave
525-
526+
527+
end subroutine s_generate_wave
528+
526529
!> Deallocation procedures for the module
527530
subroutine s_finalize_initial_condition_module() ! ---------------------
528531

0 commit comments

Comments
 (0)