From db228151ecad701506b9eb3efa85382d419e4971 Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Sat, 25 Mar 2023 09:35:15 -0700 Subject: [PATCH 01/17] mixing layer with instability waves --- examples/3D_turb_mixing/case.py | 136 +++ src/pre_process/m_assign_patches.fpp | 6 + src/pre_process/m_create_patches.fpp | 1216 +++++++++++++++++++++++ src/pre_process/m_global_parameters.fpp | 8 + src/pre_process/m_initial_condition.fpp | 4 +- src/pre_process/m_mpi_proxy.fpp | 4 +- src/pre_process/m_start_up.fpp | 2 +- toolchain/mfc/run/case_dicts.py | 3 +- 8 files changed, 1374 insertions(+), 5 deletions(-) create mode 100644 examples/3D_turb_mixing/case.py diff --git a/examples/3D_turb_mixing/case.py b/examples/3D_turb_mixing/case.py new file mode 100644 index 0000000000..e1862ccbe8 --- /dev/null +++ b/examples/3D_turb_mixing/case.py @@ -0,0 +1,136 @@ +#!/usr/bin/env python3 + +import math +import json + +# SURROUNDING FLOW ============================================================= +# Nondimensional parameters +Re0 = 50. # Reynolds number +M0 = 0.2 # Mach number + +# Fluid properties +gamma = 1.4 +pi_inf = 0. + +# Free stream velocity & pressure +u0 = 1. +pres0 = 1./(gamma*M0**2) + +# Domain size +Lx = 59.0 +Ly = 59.0 +Lz = 59.0 + +# Number of grid cells +Nx = 191 +Ny = 191 +Nz = 191 + +# Grid spacing +dx = Lx/float(Nx) +dy = Ly/float(Ny) +dz = Lz/float(Nz) + +# Time advancement +cfl = 0.1 +T = 100. +dt = cfl*dx/u0 +Ntfinal = int(T/dt) +Ntstart = 0 +Nfiles = 100 +t_save = int(math.ceil((Ntfinal-0)/float(Nfiles))) +Nt = t_save*Nfiles +t_step_start = Ntstart +t_step_stop = int(Nt) + +# ============================================================================== + + +# Configuring case dictionary +print(json.dumps({ + # Logistics ================================================================ + 'case_dir' : '\'.\'', + 'run_time_info' : 'T', + # ========================================================================== + + # Computational Domain Parameters ========================================== + 'x_domain%beg' : 0., + 'x_domain%end' : Lx, + 'y_domain%beg' : -Ly/2., + 'y_domain%end' : Ly/2., + 'z_domain%beg' : 0., + 'z_domain%end' : Lz, + 'm' : Nx, + 'n' : Ny, + 'p' : Nz, + 'dt' : dt, + 't_step_start' : t_step_start, + 't_step_stop' : t_step_stop, + 't_step_save' : t_save, + # ========================================================================== + + # Simulation Algorithm Parameters ========================================== + 'num_patches' : 1, + 'model_eqns' : 2, + 'num_fluids' : 1, + 'adv_alphan' : 'T', + 'time_stepper' : 3, + 'weno_vars' : 2, + 'weno_order' : 5, + 'weno_eps' : 1.E-16, + 'weno_Re_flux' : 'T', + 'mapped_weno' : 'T', + 'riemann_solver' : 2, + 'wave_speeds' : 1, + 'avg_state' : 2, + 'bc_x%beg' : -1, + 'bc_x%end' : -1, + 'bc_y%beg' : -5, + 'bc_y%end' : -5, + 'bc_z%beg' : -1, + 'bc_z%end' : -1, + # ========================================================================== + + # Formatted Database Files Structure Parameters ============================ + 'format' : 1, + 'precision' : 2, + 'cons_vars_wrt' :'T', + 'prim_vars_wrt' :'T', + 'parallel_io' :'T', + 'fd_order' : 1, + 'omega_wrt(1)' :'T', + 'omega_wrt(2)' :'T', + 'omega_wrt(3)' :'T', + 'qm_wrt' :'T', + # ========================================================================== + + # Patch 1 ================================================================== + 'patch_icpp(1)%geometry' : 9, + 'patch_icpp(1)%x_centroid' : Lx/2., + 'patch_icpp(1)%y_centroid' : 0., + 'patch_icpp(1)%z_centroid' : Lz/2., + 'patch_icpp(1)%length_x' : Lx, + 'patch_icpp(1)%length_y' : Ly, + 'patch_icpp(1)%length_z' : Lz, + 'patch_icpp(1)%alpha_rho(1)' : 1., + 'patch_icpp(1)%alpha(1)' : 1., + 'patch_icpp(1)%vel(1)' : u0, + 'patch_icpp(1)%vel(2)' : 0., + 'patch_icpp(1)%vel(3)' : 0., + 'patch_icpp(1)%pres' : pres0, + # ========================================================================== + + # Mixing layer ============================================================= + 'vel_profile' : 'T', + 'instability_wave' : 'T', + # ========================================================================== + + # Fluids Physical Parameters =============================================== + # Surrounding liquid + 'fluid_pp(1)%gamma' : 1./(gamma-1.), + 'fluid_pp(1)%pi_inf' : 0., + 'fluid_pp(1)%Re(1)' : Re0, + # ========================================================================= +})) + +# ============================================================================== diff --git a/src/pre_process/m_assign_patches.fpp b/src/pre_process/m_assign_patches.fpp index d88bf6a09b..26b717b663 100644 --- a/src/pre_process/m_assign_patches.fpp +++ b/src/pre_process/m_assign_patches.fpp @@ -641,6 +641,12 @@ contains + (1d0 - eta)*orig_prim_vf(i + cont_idx%end)) end do + if (vel_profile) then + q_prim_vf(1 + cont_idx%end)%sf(j, k, l) = & + (eta*patch_icpp(patch_id)%vel(1)*tanh(y_cc(k)) & + + (1d0 - eta)*orig_prim_vf(1 + cont_idx%end)) + end if + ! Pressure q_prim_vf(E_idx)%sf(j, k, l) = & (eta*patch_icpp(patch_id)%pres & diff --git a/src/pre_process/m_create_patches.fpp b/src/pre_process/m_create_patches.fpp index 8699fc44ce..cd32081eff 100644 --- a/src/pre_process/m_create_patches.fpp +++ b/src/pre_process/m_create_patches.fpp @@ -124,6 +124,1222 @@ contains end subroutine s_perturb_surrounding_flow ! ---------------------------- +!==================================================================== + subroutine s_superposition_instability_wave() ! ------------------------------ + real(kind(0d0)), dimension(5,0:m,0:n,0:p) :: wave,wave1,wave2,wave_tmp + real(kind(0d0)) :: tr,ti + integer :: i,j,k + + write(*,*) "generate instability waves ..." + + wave = 0d0 + wave1 = 0d0 + wave2 = 0d0 + if (p .eq. 0) then + call s_instability_wave(2*pi*4.0/59.0,0d0,tr,ti,wave_tmp,0d0) + wave = wave + wave_tmp + call s_instability_wave(2*pi*2.0/59.0,0d0,tr,ti,wave_tmp,0d0) + wave = wave + wave_tmp + call s_instability_wave(2*pi*1.0/59.0,0d0,tr,ti,wave_tmp,0d0) + wave = wave + wave_tmp + wave = wave*0.05 + else + call s_instability_wave(2*pi*4.0/59.0,0d0,tr,ti,wave_tmp,0d0) + wave1 = wave1 + wave_tmp + call s_instability_wave(2*pi*2.0/59.0,0d0,tr,ti,wave_tmp,0d0) + wave1 = wave1 + wave_tmp + call s_instability_wave(2*pi*1.0/59.0,0d0,tr,ti,wave_tmp,0d0) + wave1 = wave1 + wave_tmp + + call s_instability_wave(2*pi*4.0/59.0, 2*pi*4.0/59.0,tr,ti,wave_tmp,11d0/31d0*2*pi) + wave2 = wave2 + wave_tmp + call s_instability_wave(2*pi*2.0/59.0, 2*pi*2.0/59.0,tr,ti,wave_tmp,13d0/31d0*2*pi) + wave2 = wave2 + wave_tmp + call s_instability_wave(2*pi*1.0/59.0, 2*pi*1.0/59.0,tr,ti,wave_tmp,17d0/31d0*2*pi) + wave2 = wave2 + wave_tmp + call s_instability_wave(2*pi*4.0/59.0,-2*pi*4.0/59.0,tr,ti,wave_tmp,19d0/31d0*2*pi) + wave2 = wave2 + wave_tmp + call s_instability_wave(2*pi*2.0/59.0,-2*pi*2.0/59.0,tr,ti,wave_tmp,23d0/31d0*2*pi) + wave2 = wave2 + wave_tmp + call s_instability_wave(2*pi*1.0/59.0,-2*pi*1.0/59.0,tr,ti,wave_tmp,29d0/31d0*2*pi) + wave2 = wave2 + wave_tmp + ! call s_instability_wave(2*pi*4.0/59.0, 2*pi*4.0/59.0,tr,ti,wave_tmp,0.8147d0*2*pi) + ! wave2 = wave2 + wave_tmp + ! call s_instability_wave(2*pi*2.0/59.0, 2*pi*2.0/59.0,tr,ti,wave_tmp,0.9058d0*2*pi) + ! wave2 = wave2 + wave_tmp + ! call s_instability_wave(2*pi*1.0/59.0, 2*pi*1.0/59.0,tr,ti,wave_tmp,0.0270d0*2*pi) + ! wave2 = wave2 + wave_tmp + ! call s_instability_wave(2*pi*4.0/59.0,-2*pi*4.0/59.0,tr,ti,wave_tmp,0.5324d0*2*pi) + ! wave2 = wave2 + wave_tmp + ! call s_instability_wave(2*pi*2.0/59.0,-2*pi*2.0/59.0,tr,ti,wave_tmp,0.2975d0*2*pi) + ! wave2 = wave2 + wave_tmp + ! call s_instability_wave(2*pi*1.0/59.0,-2*pi*1.0/59.0,tr,ti,wave_tmp,0.7634d0*2*pi) + ! wave2 = wave2 + wave_tmp + + wave = 0.05*wave1+0.15*wave2 + end if + + do k = 0, p + do j = 0, n + do i = 0, m + q_prim_vf(mom_idx%beg )%sf(i,j,k) = q_prim_vf(mom_idx%beg )%sf(i,j,k)+wave(2,i,j,k) + q_prim_vf(mom_idx%beg+1)%sf(i,j,k) = q_prim_vf(mom_idx%beg+1)%sf(i,j,k)+wave(3,i,j,k) + if (p .gt. 0) then + q_prim_vf(mom_idx%beg+2)%sf(i,j,k) = q_prim_vf(mom_idx%beg+2)%sf(i,j,k)+wave(4,i,j,k) + end if + end do + end do + end do + + end subroutine s_superposition_instability_wave + +!==================================================================== + subroutine s_instability_wave(alpha,beta,tr,ti,wave,shift) + real(kind(0d0)),intent(in) :: alpha, beta !< spatial wavenumbers + real(kind(0d0)),dimension(0:n) :: rho_mean, u_mean, t_mean !< mean profiles + real(kind(0d0)),dimension(0:n) :: drho_mean, du_mean, dt_mean !< mean profiles + real(kind(0d0)),dimension(0:n,0:n) :: d + real(kind(0d0)),dimension(0:5*(n+1)-1,0:5*(n+1)-1) :: ar,ai,br,bi,ci + real(kind(0d0)),dimension(0:5*(n+1)-1,0:5*(n+1)-1) :: zr,zi + real(kind(0d0)),dimension(0:5*(n+1)-1) :: wr,wi,fv1,fv2,fv3 + real(kind(0d0)) :: tr,ti + real(kind(0d0)),dimension(0:5*(n+1)-1) :: vr,vi,vnr,vni + real(kind(0d0)),dimension(5,0:m,0:n,0:p) :: wave + real(kind(0d0)) :: shift + + integer :: ierr + integer :: j, k, l !< generic loop iterators + integer :: ii, jj !< block matrix indicies + + real(kind(0d0)) :: gam,pi_inf,rho1,mach,c1 + + gam = 1.+1./fluid_pp(1)%gamma + pi_inf = fluid_pp(1)%pi_inf*(gam-1.)/gam + rho1 = patch_icpp(1)%alpha_rho(1)/patch_icpp(1)%alpha(1) + c1 = sqrt((gam*(patch_icpp(1)%pres+pi_inf))/rho1) + mach = 1./c1 + + ! Assign mean profiles + do j=0,n + u_mean(j)=tanh(y_cc(j)) + t_mean(j)=1+0.5*(gam-1)*mach**2*(1-u_mean(j)**2) + rho_mean(j)=1/T_mean(j) + end do + + ! Compute differential operator + dy = y_cc(1)-y_cc(0) + d=0d0 + d(1,0)=-1/(2*dy) + d(1,2)= 1/(2*dy) + do j=2,n-2 + d(j,j-2)= 1/(12*dy) + d(j,j-1)=-8/(12*dy) + d(j,j+1)= 8/(12*dy) + d(j,j+2)=-1/(12*dy) + end do + d(n-1,n-2)=-1/(2*dy) + d(n-1,n) = 1/(2*dy) + + ! Compute y-derivatives of rho, u, T + do j=0,n + drho_mean(j)=0 + du_mean(j)=0 + dt_mean(j)=0 + do k=0,n + drho_mean(j) = drho_mean(j)+d(j,k)*rho_mean(k) + du_mean(j) = du_mean(j)+d(j,k)*u_mean(k) + dt_mean(j) = dt_mean(j)+d(j,k)*t_mean(k) + end do + end do + + ! Compute B and C -> A=B+C + br=0d0 + bi=0d0 + ci=0d0 + do j=0,n + ii = 1; jj = 1; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); + ii = 1; jj = 2; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*rho_mean(j); + ii = 1; jj = 3; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -drho_mean(j); + ii = 1; jj = 4; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = beta*rho_mean(j); + + ii = 2; jj = 1; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*t_mean(j)/(rho_mean(j)*gam*mach**2); + ii = 2; jj = 2; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); + ii = 2; jj = 3; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -du_mean(j); + ii = 2; jj = 5; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha/(gam*mach**2); + + ii = 3; jj = 1; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -dt_mean(j)/(rho_mean(j)*gam*mach**2); + ii = 3; jj = 3; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); + ii = 3; jj = 5; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -drho_mean(j)/(rho_mean(j)*gam*mach**2); + + ii = 4; jj = 1; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = beta*t_mean(j)/(rho_mean(j)*gam*mach**2); + ii = 4; jj = 4; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); + ii = 4; jj = 5; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = beta/(gam*mach**2); + + ii = 5; jj = 2; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = (gam-1)*alpha/rho_mean(j); + ii = 5; jj = 3; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -dt_mean(j); + ii = 5; jj = 4; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = (gam-1)*beta/rho_mean(j); + ii = 5; jj = 5; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); + + do k=0,n + ii = 1; jj = 3; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -rho_mean(j)*d(j,k); + ii = 3; jj = 1; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -t_mean(j)*d(j,k)/(rho_mean(j)*gam*mach**2); + ii = 3; jj = 5; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -d(j,k)/(gam*mach**2); + ii = 5; jj = 3; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -(gam-1)*d(j,k)/rho_mean(j); + end do + end do + ar = br + ai = bi+ci + + ! Compute Eigenvalues & Eigenfunctions + call cg(5*(n+1),5*(n+1),ar,ai,wr,wi,zr,zi,fv1,fv2,fv3,ierr) + + ! Find the most unstable wave and normalize it + call find_unstable_mode(5*(n+1),wr,wi,zr,zi,tr,ti,vr,vi) + + ! Normalize + call normalize_eigvec(5*(n+1),vr,vi,vnr,vni) + + ! Generate instability waves + call generate_wave(5*(n+1),vnr,vni,alpha,beta,wave,shift) + + end subroutine s_instability_wave ! ---------------------------- + +!==================================================================== + subroutine generate_wave(nl,vnr,vni,alpha,beta,wave,shift) + integer nl + real(kind(0d0)) :: alpha,beta,ang + real(kind(0d0)), dimension(0:nl-1) :: vnr,vni + real(kind(0d0)), dimension(5,0:m,0:n,0:p) :: wave + real(kind(0d0)) :: shift + integer i,j,k + + do i=0,m + do j=0,n + do k=0,p + if (beta .eq. 0) then + ang = alpha*x_cc(i) + else + ang = alpha*x_cc(i)+beta*z_cc(k)+shift + end if + + wave(1,i,j,k) = vnr(j)*cos(ang) & + -vni(j)*sin(ang) ! rho + wave(2,i,j,k) = vnr((n+1)+j)*cos(ang) & + -vni((n+1)+j)*sin(ang) ! u + wave(3,i,j,k) = vnr(2*(n+1)+j)*cos(ang) & + -vni(2*(n+1)+j)*sin(ang) ! v + wave(4,i,j,k) = vnr(3*(n+1)+j)*cos(ang) & + -vni(3*(n+1)+j)*sin(ang) ! w + wave(5,i,j,k) = vnr(4*(n+1)+j)*cos(ang) & + -vni(4*(n+1)+j)*sin(ang) ! T + end do + end do + end do + + end subroutine generate_wave + +!==================================================================== + subroutine normalize_eigvec(nl,vr,vi,vnr,vni) + integer nl + real(kind(0d0)), dimension(nl) :: vr,vi,vnr,vni + real(kind(0d0)) :: norm + real(kind(0d0)) :: tr,ti,cr,ci + integer i,idx + + norm = 0d0 + do i=1,nl + if (dsqrt(vr(i)**2+vi(i)**2) .gt. dabs(norm)) then + idx = i + if (vr(i) .gt. 0) then + norm = sqrt(vr(i)**2+vi(i)**2) + else + norm =-sqrt(vr(i)**2+vi(i)**2) + end if + end if + end do + + vnr = vr/norm + vni = vi/norm + + tr = vnr(idx) + ti = vni(idx) + do i=1,nl + call cdiv(vnr(i),vni(i),tr,ti,cr,ci) + vnr(i) = cr + vni(i) = ci + end do + + end subroutine normalize_eigvec + +!==================================================================== + subroutine find_unstable_mode(nl,wr,wi,zr,zi,tr,ti,vr,vi) + integer nl + real(kind(0d0)), dimension(nl) :: wr,wi + real(kind(0d0)), dimension(nl,nl) :: zr,zi + real(kind(0d0)) :: tr,ti + real(kind(0d0)), dimension(nl) :: vr,vi + integer i,k + + k=1 + do i=2,nl + if (wi(i) .gt. wi(k)) then + k = i + end if + end do + + tr = wr(k) + ti = wi(k) + vr = zr(:,k) + vi = zi(:,k) + + end subroutine find_unstable_mode + +!==================================================================== +! Subroutines & functions for computing eigenvalues and eigenvectors +! which are modified from EISPACK (https://netlib.org/eispack/) +! 1) cg +! 2) cbal +! 3) corth +! 4) comqr2 +! 5) csroot +! 6) cdiv +! 7) pythag +!==================================================================== + subroutine cg(nm,nl,ar,ai,wr,wi,zr,zi,fv1,fv2,fv3,ierr) +! this subroutine calls the recommended sequence of +! subroutines from the eigensystem subroutine package (eispack) +! to find the eigenvalues and eigenvectors (if desired) +! of a complex general matrix. +! +! on input +! +! nm must be set to the row dimension of the two-dimensional +! array parameters as declared in the calling program +! dimension statement. +! +! nl is the order of the matrix a=(ar,ai). +! +! ar and ai contain the real and imaginary parts, +! respectively, of the complex general matrix. +! +! on output +! +! wr and wi contain the real and imaginary parts, +! respectively, of the eigenvalues. +! +! zr and zi contain the real and imaginary parts, +! respectively, of the eigenvectors. +! +! ierr is an integer output variable set equal to an error +! completion code described in the documentation for comqr +! and comqr2. the normal completion code is zero. +! +! fv1, fv2, and fv3 are temporary storage arrays. +! +! questions and comments should be directed to burton s. garbow, +! mathematics and computer science div, argonne national laboratory +! +! this version dated august 1983. +! +! ------------------------------------------------------------------ + integer nm,nl,is1,is2,ierr + real(kind(0d0)), dimension(nm,nl) :: ar,ai,zr,zi + real(kind(0d0)), dimension(nl) :: wr,wi,fv1,fv2,fv3 + + if (nl .le. nm) go to 10 + ierr = 10*nl + go to 50 + +10 call cbal(nm,nl,ar,ai,is1,is2,fv1) + call corth(nm,nl,is1,is2,ar,ai,fv2,fv3) + call comqr2(nm,nl,is1,is2,fv2,fv3,ar,ai,wr,wi,zr,zi,ierr) + if (ierr .ne. 0) go to 50 + call cbabk2(nm,nl,is1,is2,fv1,nl,zr,zi) + +50 return + end subroutine cg + +!==================================================================== + subroutine cbal(nm,nl,ar,ai,low,igh,scale) +! this subroutine is a translation of the algol procedure +! cbalance, which is a complex version of balance, +! num. math. 13, 293-304(1969) by parlett and reinsch. +! handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). +! +! this subroutine balances a complex matrix and isolates +! eigenvalues whenever possible. +! +! on input +! +! nm must be set to the row dimension of two-dimensional +! array parameters as declared in the calling program +! dimension statement. +! +! nl is the order of the matrix. +! +! ar and ai contain the real and imaginary parts, +! respectively, of the complex matrix to be balanced. +! +! on output +! +! ar and ai contain the real and imaginary parts, +! respectively, of the balanced matrix. +! +! low and igh are two integers such that ar(i,j) and ai(i,j) +! are equal to zero if +! (1) i is greater than j and +! (2) j=1,...,low-1 or i=igh+1,...,nl. +! +! scale contains information determining the +! permutations and scaling factors used. +! +! suppose that the principal submatrix in rows low through igh +! has been balanced, that p(j) denotes the index interchanged +! with j during the permutation step, and that the elements +! of the diagonal matrix used are denoted by d(i,j). then +! scale(j) = p(j), for j = 1,...,low-1 +! = d(j,j) j = low,...,igh +! = p(j) j = igh+1,...,nl. +! the order in which the interchanges are made is nl to igh+1, +! then 1 to low-1. +! +! note that 1 is returned for igh if igh is zero formally. +! +! the algol procedure exc contained in cbalance appears in +! cbal in line. (note that the algol roles of identifiers +! k,l have been reversed.) +! +! arithmetic is real throughout. +! +! questions and comments should be directed to burton s. garbow, +! mathematics and computer science div, argonne national laboratory +! +! this version dated august 1983. +! +! ------------------------------------------------------------------ + integer i,j,k,l,ml,nl,jj,nm,igh,low,iexc + real(kind(0d0)), dimension(nm,nl) :: ar,ai + real(kind(0d0)), dimension(nl) :: scale + real(kind(0d0)) :: c,f,g,r,s,b2,radix + logical noconv + + radix = 16.0d0 + + b2 = radix * radix + k = 1 + l = nl + go to 100 +! .......... in-line procedure for row and +! column exchange .......... +20 scale(ml) = j + if (j .eq. ml) go to 50 + + do 30 i = 1, l + f = ar(i,j) + ar(i,j) = ar(i,ml) + ar(i,ml) = f + f = ai(i,j) + ai(i,j) = ai(i,ml) + ai(i,ml) = f +30 continue + + do 40 i = k, nl + f = ar(j,i) + ar(j,i) = ar(ml,i) + ar(ml,i) = f + f = ai(j,i) + ai(j,i) = ai(ml,i) + ai(ml,i) = f +40 continue + +50 go to (80,130), iexc +! .......... search for rows isolating an eigenvalue +! and push them down .......... +80 if (l .eq. 1) go to 280 + l = l - 1 +! .......... for j=l step -1 until 1 do -- .......... +100 do 120 jj = 1, l + j = l + 1 - jj + + do 110 i = 1, l + if (i .eq. j) go to 110 + if (ar(j,i) .ne. 0.0d0 .or. ai(j,i) .ne. 0.0d0) go to 120 +110 continue + + ml = l + iexc = 1 + go to 20 +120 continue + + go to 140 +! .......... search for columns isolating an eigenvalue +! and push them left .......... +130 k = k + 1 + +140 do 170 j = k, l + + do 150 i = k, l + if (i .eq. j) go to 150 + if (ar(i,j) .ne. 0.0d0 .or. ai(i,j) .ne. 0.0d0) go to 170 +150 continue + + ml = k + iexc = 2 + go to 20 +170 continue +! .......... now balance the submatrix in rows k to l .......... + do 180 i = k, l + scale(i) = 1.0d0 +180 continue +! .......... iterative loop for norm reduction .......... +190 noconv = .false. + + do 270 i = k, l + c = 0.0d0 + r = 0.0d0 + + do 200 j = k, l + if (j .eq. i) go to 200 + c = c + dabs(ar(j,i)) + dabs(ai(j,i)) + r = r + dabs(ar(i,j)) + dabs(ai(i,j)) +200 continue +! .......... guard against zero c or r due to underflow .......... + if (c .eq. 0.0d0 .or. r .eq. 0.0d0) go to 270 + g = r / radix + f = 1.0d0 + s = c + r +210 if (c .ge. g) go to 220 + f = f * radix + c = c * b2 + go to 210 +220 g = r * radix +230 if (c .lt. g) go to 240 + f = f / radix + c = c / b2 + go to 230 +! .......... now balance .......... +240 if ((c + r) / f .ge. 0.95d0 * s) go to 270 + g = 1.0d0 / f + scale(i) = scale(i) * f + noconv = .true. + + do 250 j = k, nl + ar(i,j) = ar(i,j) * g + ai(i,j) = ai(i,j) * g +250 continue + + do 260 j = 1, l + ar(j,i) = ar(j,i) * f + ai(j,i) = ai(j,i) * f +260 continue + +270 continue + + if (noconv) go to 190 + +280 low = k + igh = l + return + end subroutine cbal + +!==================================================================== + subroutine corth(nm,nl,low,igh,ar,ai,ortr,orti) +! this subroutine is a translation of a complex analogue of +! the algol procedure orthes, num. math. 12, 349-368(1968) +! by martin and wilkinson. +! handbook for auto. comp., vol.ii-linear algebra, 339-358(1971). +! +! given a complex general matrix, this subroutine +! reduces a submatrix situated in rows and columns +! low through igh to upper hessenberg form by +! unitary similarity transformations. +! +! on input +! +! nm must be set to the row dimension of two-dimensional +! array parameters as declared in the calling program +! dimension statement. +! +! nl is the order of the matrix. +! +! low and igh are integers determined by the balancing +! subroutine cbal. if cbal has not been used, +! set low=1, igh=nl. +! +! ar and ai contain the real and imaginary parts, +! respectively, of the complex input matrix. +! +! on output +! +! ar and ai contain the real and imaginary parts, +! respectively, of the hessenberg matrix. information +! about the unitary transformations used in the reduction +! is stored in the remaining triangles under the +! hessenberg matrix. +! +! ortr and orti contain further information about the +! transformations. only elements low through igh are used. +! +! calls pythag for dsqrt(a*a + b*b) . +! +! questions and comments should be directed to burton s. garbow, +! mathematics and computer science div, argonne national laboratory +! +! this version dated august 1983. +! +! ------------------------------------------------------------------ + integer i,j,ml,nl,ii,jj,la,mp,nm,igh,kp1,low + real(kind(0d0)),dimension(nm,nl) :: ar,ai + real(kind(0d0)),dimension(igh) :: ortr,orti + real(kind(0d0)) :: f,g,h,fi,fr,scale,c + + integer mll + mll = 6 + + la = igh - 1 + kp1 = low + 1 + if (la .lt. kp1) go to 200 + + do 180 ml = kp1, la + h = 0.0d0 + ortr(ml) = 0.0d0 + orti(ml) = 0.0d0 + scale = 0.0d0 +! .......... scale column (algol tol then not needed) .......... + do 90 i = ml, igh + scale = scale + dabs(ar(i,ml-1)) + dabs(ai(i,ml-1)) +90 continue + if (scale .eq. 0d0) go to 180 + mp = ml + igh +! .......... for i=igh step -1 until ml do -- .......... + do 100 ii = ml, igh + i = mp - ii + ortr(i) = ar(i,ml-1) / scale + orti(i) = ai(i,ml-1) / scale + h = h + ortr(i) * ortr(i) + orti(i) * orti(i) +100 continue +! + g = dsqrt(h) + call pythag(ortr(ml),orti(ml),f) + if (f .eq. 0d0) go to 103 + h = h + f * g + g = g / f + ortr(ml) = (1.0d0 + g) * ortr(ml) + orti(ml) = (1.0d0 + g) * orti(ml) + go to 105 + +103 ortr(ml) = g + ar(ml,ml-1) = scale +! .......... form (i-(u*ut)/h) * a .......... +105 do 130 j = ml, nl + fr = 0.0d0 + fi = 0.0d0 +! .......... for i=igh step -1 until ml do -- .......... + do 110 ii = ml, igh + i = mp - ii + fr = fr + ortr(i) * ar(i,j) + orti(i) * ai(i,j) + fi = fi + ortr(i) * ai(i,j) - orti(i) * ar(i,j) +110 continue +! + fr = fr / h + fi = fi / h +! + do 120 i = ml, igh + ar(i,j) = ar(i,j) - fr * ortr(i) + fi * orti(i) + ai(i,j) = ai(i,j) - fr * orti(i) - fi * ortr(i) +120 continue +! +130 continue +! .......... form (i-(u*ut)/h)*a*(i-(u*ut)/h) .......... + do 160 i = 1, igh + fr = 0.0d0 + fi = 0.0d0 +! .......... for j=igh step -1 until ml do -- .......... + do 140 jj = ml, igh + j = mp - jj + fr = fr + ortr(j) * ar(i,j) - orti(j) * ai(i,j) + fi = fi + ortr(j) * ai(i,j) + orti(j) * ar(i,j) +140 continue +! + fr = fr / h + fi = fi / h +! + do 150 j = ml, igh + ar(i,j) = ar(i,j) - fr * ortr(j) - fi * orti(j) + ai(i,j) = ai(i,j) + fr * orti(j) - fi * ortr(j) +150 continue +! +160 continue +! + ortr(ml) = scale * ortr(ml) + orti(ml) = scale * orti(ml) + ar(ml,ml-1) = -g * ar(ml,ml-1) + ai(ml,ml-1) = -g * ai(ml,ml-1) +180 continue +! +200 return + end subroutine corth + +!==================================================================== + subroutine comqr2(nm,nl,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr) +! MESHED overflow control WITH vectors of isolated roots (10/19/89 BSG) +! MESHED overflow control WITH triangular multiply (10/30/89 BSG) +! +! this subroutine is a translation of a unitary analogue of the +! algol procedure comlr2, num. math. 16, 181-204(1970) by peters +! and wilkinson. +! handbook for auto. comp., vol.ii-linear algebra, 372-395(1971). +! the unitary analogue substitutes the qr algorithm of francis +! (comp. jour. 4, 332-345(1962)) for the lr algorithm. +! +! this subroutine finds the eigenvalues and eigenvectors +! of a complex upper hessenberg matrix by the qr +! method. the eigenvectors of a complex general matrix +! can also be found if corth has been used to reduce +! this general matrix to hessenberg form. +! +! on input +! +! nm must be set to the row dimension of two-dimensional +! array parameters as declared in the calling program +! dimension statement. +! +! nl is the order of the matrix. +! +! low and igh are integers determined by the balancing +! subroutine cbal. if cbal has not been used, +! set low=1, igh=nl. +! +! ortr and orti contain information about the unitary trans- +! formations used in the reduction by corth, if performed. +! only elements low through igh are used. if the eigenvectors +! of the hessenberg matrix are desired, set ortr(j) and +! orti(j) to 0.0d0 for these elements. +! +! hr and hi contain the real and imaginary parts, +! respectively, of the complex upper hessenberg matrix. +! their lower triangles below the subdiagonal contain further +! information about the transformations which were used in the +! reduction by corth, if performed. if the eigenvectors of +! the hessenberg matrix are desired, these elements may be +! arbitrary. +! +! on output +! +! ortr, orti, and the upper hessenberg portions of hr and hi +! have been destroyed. +! +! wr and wi contain the real and imaginary parts, +! respectively, of the eigenvalues. if an error +! exit is made, the eigenvalues should be correct +! for indices ierr+1,...,nl. +! +! zr and zi contain the real and imaginary parts, +! respectively, of the eigenvectors. the eigenvectors +! are unnormalized. if an error exit is made, none of +! the eigenvectors has been found. +! +! ierr is set to +! zero for normal return, +! j if the limit of 30*nl iterations is exhausted +! while the j-th eigenvalue is being sought. +! +! calls cdiv for complex division. +! calls csroot for complex square root. +! calls pythag for dsqrt(a*a + b*b) . +! +! questions and comments should be directed to burton s. garbow, +! mathematics and computer science div, argonne national laboratory +! +! this version dated october 1989. +! +! ------------------------------------------------------------------ + integer i,j,k,l,ml,nl,en,ii,jj,ll,nm,nn,igh,ip1,& + itn,its,low,lp1,enm1,iend,ierr + real(kind(0d0)),dimension(nm,nl) :: hr,hi,zr,zi + real(kind(0d0)),dimension(nl) :: wr,wi + real(kind(0d0)),dimension(igh) :: ortr,orti + real(kind(0d0)) :: si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr,& + norm,tst1,tst2,c,d +! + ierr = 0 +! .......... initialize eigenvector matrix .......... + do 101 j = 1, nl +! + do 100 i = 1, nl + zr(i,j) = 0.0d0 + zi(i,j) = 0.0d0 +100 continue + zr(j,j) = 1.0d0 +101 continue +! .......... form the matrix of accumulated transformations +! from the information left by corth .......... + iend = igh - low - 1 + if (iend .lt. 0) go to 180 + if (iend .eq. 0) go to 150 + if (iend .gt. 0) go to 105 +! .......... for i=igh-1 step -1 until low+1 do -- .......... +105 do 140 ii = 1, iend + i = igh - ii + if (dabs(ortr(i)) .eq. 0d0 .and. dabs(orti(i)) .eq. 0d0) go to 140 + if (dabs(hr(i,i-1)) .eq. 0d0 .and. dabs(hi(i,i-1)) .eq. 0d0) go to 140 +! .......... norm below is negative of h formed in corth .......... + norm = hr(i,i-1) * ortr(i) + hi(i,i-1) * orti(i) + ip1 = i + 1 + + do 110 k = ip1, igh + ortr(k) = hr(k,i-1) + orti(k) = hi(k,i-1) +110 continue +! + do 130 j = i, igh + sr = 0.0d0 + si = 0.0d0 +! + do 115 k = i, igh + sr = sr + ortr(k) * zr(k,j) + orti(k) * zi(k,j) + si = si + ortr(k) * zi(k,j) - orti(k) * zr(k,j) +115 continue +! + sr = sr / norm + si = si / norm +! + do 120 k = i, igh + zr(k,j) = zr(k,j) + sr * ortr(k) - si * orti(k) + zi(k,j) = zi(k,j) + sr * orti(k) + si * ortr(k) +120 continue +! +130 continue +! +140 continue +! .......... create real subdiagonal elements .......... +150 l = low + 1 +! + do 170 i = l, igh + ll = min0(i+1,igh) + if (dabs(hi(i,i-1)) .eq. 0d0) go to 170 + call pythag(hr(i,i-1),hi(i,i-1),norm) + yr = hr(i,i-1) / norm + yi = hi(i,i-1) / norm + hr(i,i-1) = norm + hi(i,i-1) = 0.0d0 +! + do 155 j = i, nl + si = yr * hi(i,j) - yi * hr(i,j) + hr(i,j) = yr * hr(i,j) + yi * hi(i,j) + hi(i,j) = si +155 continue +! + do 160 j = 1, ll + si = yr * hi(j,i) + yi * hr(j,i) + hr(j,i) = yr * hr(j,i) - yi * hi(j,i) + hi(j,i) = si +160 continue +! + do 165 j = low, igh + si = yr * zi(j,i) + yi * zr(j,i) + zr(j,i) = yr * zr(j,i) - yi * zi(j,i) + zi(j,i) = si +165 continue + +170 continue +! .......... store roots isolated by cbal .......... +180 do 200 i = 1, nl + if (i .ge. low .and. i .le. igh) go to 200 + wr(i) = hr(i,i) + wi(i) = hi(i,i) +200 continue +! + en = igh + tr = 0.0d0 + ti = 0.0d0 + itn = 30*nl +! .......... search for next eigenvalue .......... +220 if (en .lt. low) go to 680 + its = 0 + enm1 = en - 1 +! .......... look for single small sub-diagonal element +! for l=en step -1 until low do -- .......... +240 do 260 ll = low, en + l = en + low - ll + if (l .eq. low) go to 300 + tst1 = dabs(hr(l-1,l-1)) + dabs(hi(l-1,l-1)) & + + dabs(hr(l,l)) + dabs(hi(l,l)) + tst2 = tst1 + dabs(hr(l,l-1)) + if (tst2 .eq. tst1) go to 300 +260 continue +! .......... form shift .......... +300 if (l .eq. en) go to 660 + if (itn .eq. 0) go to 1000 + if (its .eq. 10 .or. its .eq. 20) go to 320 + sr = hr(en,en) + si = hi(en,en) + xr = hr(enm1,en) * hr(en,enm1) + xi = hi(enm1,en) * hr(en,enm1) + if (xr .eq. 0.0d0 .and. xi .eq. 0.0d0) go to 340 + yr = (hr(enm1,enm1) - sr) / 2.0d0 + yi = (hi(enm1,enm1) - si) / 2.0d0 + call csroot(yr**2-yi**2+xr,2.0d0*yr*yi+xi,zzr,zzi) + if (yr * zzr + yi * zzi .ge. 0.0d0) go to 310 + zzr = -zzr + zzi = -zzi +310 call cdiv(xr,xi,yr+zzr,yi+zzi,xr,xi) + sr = sr - xr + si = si - xi + go to 340 +! .......... form exceptional shift .......... +320 sr = dabs(hr(en,enm1)) + dabs(hr(enm1,en-2)) + si = 0.0d0 +! +340 do 360 i = low, en + hr(i,i) = hr(i,i) - sr + hi(i,i) = hi(i,i) - si +360 continue +! + tr = tr + sr + ti = ti + si + its = its + 1 + itn = itn - 1 +! .......... reduce to triangle (rows) .......... + lp1 = l + 1 +! + do 500 i = lp1, en + sr = hr(i,i-1) + hr(i,i-1) = 0.0d0 + call pythag(hr(i-1,i-1),hi(i-1,i-1),c) + call pythag(c,sr,norm) + xr = hr(i-1,i-1) / norm + wr(i-1) = xr + xi = hi(i-1,i-1) / norm + wi(i-1) = xi + hr(i-1,i-1) = norm + hi(i-1,i-1) = 0.0d0 + hi(i,i-1) = sr / norm +! + do 490 j = i, nl + yr = hr(i-1,j) + yi = hi(i-1,j) + zzr = hr(i,j) + zzi = hi(i,j) + hr(i-1,j) = xr * yr + xi * yi + hi(i,i-1) * zzr + hi(i-1,j) = xr * yi - xi * yr + hi(i,i-1) * zzi + hr(i,j) = xr * zzr - xi * zzi - hi(i,i-1) * yr + hi(i,j) = xr * zzi + xi * zzr - hi(i,i-1) * yi +490 continue +! +500 continue +! + si = hi(en,en) + if (dabs(si) .eq. 0d0) go to 540 + call pythag(hr(en,en),si,norm) + sr = hr(en,en) / norm + si = si / norm + hr(en,en) = norm + hi(en,en) = 0.0d0 + if (en .eq. nl) go to 540 + ip1 = en + 1 +! + do 520 j = ip1, nl + yr = hr(en,j) + yi = hi(en,j) + hr(en,j) = sr * yr + si * yi + hi(en,j) = sr * yi - si * yr +520 continue +! .......... inverse operation (columns) .......... +540 do 600 j = lp1, en + xr = wr(j-1) + xi = wi(j-1) +! + do 580 i = 1, j + yr = hr(i,j-1) + yi = 0.0d0 + zzr = hr(i,j) + zzi = hi(i,j) + if (i .eq. j) go to 560 + yi = hi(i,j-1) + hi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi +560 hr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr + hr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr + hi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi +580 continue +! + do 590 i = low, igh + yr = zr(i,j-1) + yi = zi(i,j-1) + zzr = zr(i,j) + zzi = zi(i,j) + zr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr + zi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi + zr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr + zi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi +590 continue + +600 continue +! + if (dabs(si) .eq. 0d0) go to 240 +! + do 630 i = 1, en + yr = hr(i,en) + yi = hi(i,en) + hr(i,en) = sr * yr - si * yi + hi(i,en) = sr * yi + si * yr +630 continue +! + do 640 i = low, igh + yr = zr(i,en) + yi = zi(i,en) + zr(i,en) = sr * yr - si * yi + zi(i,en) = sr * yi + si * yr +640 continue +! + go to 240 +! .......... a root found .......... +660 hr(en,en) = hr(en,en) + tr + wr(en) = hr(en,en) + hi(en,en) = hi(en,en) + ti + wi(en) = hi(en,en) + en = enm1 + go to 220 +! .......... all roots found. backsubstitute to find +! vectors of upper triangular form .......... +680 norm = 0.0d0 +! + do i = 1, nl + do j = i, nl + tr = dabs(hr(i,j)) + dabs(hi(i,j)) + if (tr .gt. norm) norm = tr + end do + end do +! + if (nl .eq. 1 .or. norm .eq. 0d0) go to 1001 +! .......... for en=nl step -1 until 2 do -- .......... + do 800 nn = 2, nl + en = nl + 2 - nn + xr = wr(en) + xi = wi(en) + hr(en,en) = 1.0d0 + hi(en,en) = 0.0d0 + enm1 = en - 1 +! .......... for i=en-1 step -1 until 1 do -- .......... + do 780 ii = 1, enm1 + i = en - ii + zzr = 0.0d0 + zzi = 0.0d0 + ip1 = i + 1 + + do 740 j = ip1, en + zzr = zzr + hr(i,j) * hr(j,en) - hi(i,j) * hi(j,en) + zzi = zzi + hr(i,j) * hi(j,en) + hi(i,j) * hr(j,en) +740 continue +! + yr = xr - wr(i) + yi = xi - wi(i) + if (yr .ne. 0.0d0 .or. yi .ne. 0.0d0) go to 765 + tst1 = norm + yr = tst1 +760 yr = 0.01d0 * yr + tst2 = norm + yr + if (tst2 .gt. tst1) go to 760 +765 continue + call cdiv(zzr,zzi,yr,yi,hr(i,en),hi(i,en)) +! .......... overflow control .......... + tr = dabs(hr(i,en)) + dabs(hi(i,en)) + if (tr .eq. 0.0d0) go to 780 + tst1 = tr + tst2 = tst1 + 1.0d0/tst1 + if (tst2 .gt. tst1) go to 780 + do 770 j = i, en + hr(j,en) = hr(j,en)/tr + hi(j,en) = hi(j,en)/tr +770 continue +! +780 continue +! +800 continue +! .......... end backsubstitution .......... +! .......... vectors of isolated roots .......... + do 840 i = 1, nl + if (i .ge. low .and. i .le. igh) go to 840 +! + do 820 j = I, nl + zr(i,j) = hr(i,j) + zi(i,j) = hi(i,j) +820 continue +! +840 continue +! .......... multiply by transformation matrix to give +! vectors of original full matrix. +! for j=nl step -1 until low do -- .......... + do jj = low, nl + j = nl + low - jj + ml = min0(j,igh) +! + do i = low, igh + zzr = 0.0d0 + zzi = 0.0d0 +! + do 860 k = low, ml + zzr = zzr + zr(i,k) * hr(k,j) - zi(i,k) * hi(k,j) + zzi = zzi + zr(i,k) * hi(k,j) + zi(i,k) * hr(k,j) +860 continue +! + zr(i,j) = zzr + zi(i,j) = zzi + end do + end do +! + go to 1001 +! .......... set error -- all eigenvalues have not +! converged after 30*nl iterations .......... +1000 ierr = en +1001 return + end subroutine comqr2 + +!==================================================================== + subroutine cbabk2(nm,nl,low,igh,scale,ml,zr,zi) +! + integer i,j,k,ml,nl,ii,nm,igh,low + double precision scale(nl),zr(nm,ml),zi(nm,ml) + double precision s +! +! this subroutine is a translation of the algol procedure +! cbabk2, which is a complex version of balbak, +! num. math. 13, 293-304(1969) by parlett and reinsch. +! handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). +! +! this subroutine forms the eigenvectors of a complex general +! matrix by back transforming those of the corresponding +! balanced matrix determined by cbal. +! +! on input +! +! nm must be set to the row dimension of two-dimensional +! array parameters as declared in the calling program +! dimension statement. +! +! nl is the order of the matrix. +! +! low and igh are integers determined by cbal. +! +! scale contains information determining the permutations +! and scaling factors used by cbal. +! +! ml is the number of eigenvectors to be back transformed. +! +! zr and zi contain the real and imaginary parts, +! respectively, of the eigenvectors to be +! back transformed in their first ml columns. +! +! on output +! +! zr and zi contain the real and imaginary parts, +! respectively, of the transformed eigenvectors +! in their first ml columns. +! +! questions and comments should be directed to burton s. garbow, +! mathematics and computer science div, argonne national laboratory +! +! this version dated august 1983. +! +! ------------------------------------------------------------------ +! + if (ml .eq. 0) go to 200 + if (igh .eq. low) go to 120 +! + do 110 i = low, igh + s = scale(i) +! .......... left hand eigenvectors are back transformed +! if the foregoing statement is replaced by +! s=1.0d0/scale(i). .......... + do 100 j = 1, ml + zr(i,j) = zr(i,j) * s + zi(i,j) = zi(i,j) * s +100 continue +! +110 continue +! .......... for i=low-1 step -1 until 1, +! igh+1 step 1 until nl do -- .......... +120 do 140 ii = 1, nl + i = ii + if (i .ge. low .and. i .le. igh) go to 140 + if (i .lt. low) i = low - ii + k = scale(i) + if (k .eq. i) go to 140 +! + do 130 j = 1, ml + s = zr(i,j) + zr(i,j) = zr(k,j) + zr(k,j) = s + s = zi(i,j) + zi(i,j) = zi(k,j) + zi(k,j) = s +130 continue +! +140 continue +! +200 return + end subroutine cbabk2 + +!==================================================================== + subroutine csroot(xr,xi,yr,yi) + real(kind(0d0)) :: xr,xi,yr,yi +! +! (yr,yi) = complex dsqrt(xr,xi) +! branch chosen so that yr .ge. 0.0 and sign(yi) .eq. sign(xi) +! + real(kind(0d0)) :: s,tr,ti,c + tr = xr + ti = xi + call pythag(tr,ti,c) + s = dsqrt(0.5d0*(c + dabs(tr))) + if (tr .ge. 0.0d0) yr = s + if (ti .lt. 0.0d0) s = -s + if (tr .le. 0.0d0) yi = s + if (tr .lt. 0.0d0) yr = 0.5d0*(ti/yi) + if (tr .gt. 0.0d0) yi = 0.5d0*(ti/yr) + return + end subroutine csroot + +!==================================================================== + subroutine cdiv(ar,ai,br,bi,cr,ci) + real(kind(0d0)) :: ar,ai,br,bi,cr,ci +! +! complex division, (cr,ci) = (ar,ai)/(br,bi) +! + real(kind(0d0)) :: s,ars,ais,brs,bis + s = dabs(br) + dabs(bi) + ars = ar/s + ais = ai/s + brs = br/s + bis = bi/s + s = brs**2 + bis**2 + cr = (ars*brs + ais*bis)/s + ci = (ais*brs - ars*bis)/s + return + end subroutine cdiv + +!==================================================================== + subroutine pythag(a,b,c) + real(kind(0d0)) :: a,b,c +! +! finds dsqrt(a**2+b**2) without overflow or destructive underflow +! + real(kind(0d0)) :: p,r,s,t,u + p = dmax1(dabs(a),dabs(b)) + if (p .eq. 0.0d0) go to 20 + r = (dmin1(dabs(a),dabs(b))/p)**2 +10 continue + t = 4.0d0 + r + if (t .eq. 4.0d0) go to 20 + s = r/t + u = 1.0d0 + 2.0d0*s + p = u*p + r = (s/u)**2 * r + go to 10 +20 c = p + return + end subroutine pythag + !> The line segment patch is a 1D geometry that may be used, !! for example, in creating a Riemann problem. The geometry !! of the patch is well-defined when its centroid and length diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 3bbf491b59..9375770c97 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -103,6 +103,12 @@ module m_global_parameters logical :: parallel_io !< Format of the data files integer :: precision !< Precision of output files + ! Hypertangent velocity profile + logical :: vel_profile + + ! Add instability waves to surrounding fluid + logical :: instability_wave + ! Perturb density of surrounding air so as to break symmetry of grid logical :: perturb_flow integer :: perturb_flow_fluid !< Fluid to be perturbed with perturb_flow flag @@ -252,6 +258,8 @@ contains parallel_io = .false. precision = 2 + vel_profile = .false. + instability_wave = .false. perturb_flow = .false. perturb_flow_fluid = dflt_int perturb_sph = .false. diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp index f0e36f2ef2..7a7b3caf1c 100644 --- a/src/pre_process/m_initial_condition.fpp +++ b/src/pre_process/m_initial_condition.fpp @@ -42,7 +42,7 @@ contains subroutine s_generate_initial_condition() ! ---------------------------- integer :: i !< Generic loop operator - + ! Converting the conservative variables to the primitive ones given ! preexisting initial condition data files were read in on start-up if (old_ic) then @@ -162,6 +162,8 @@ contains if (perturb_flow) call s_perturb_surrounding_flow() if (perturb_sph) call s_perturb_sphere() + if (instability_wave) call s_superposition_instability_wave() + ! Converting the primitive variables to the conservative ones call s_convert_primitive_to_conservative_variables(q_prim_vf, & q_cons_vf) diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index f47d4a37d3..a1660189a8 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/src/pre_process/m_mpi_proxy.fpp @@ -54,8 +54,8 @@ contains #:for VAR in [ 'old_grid','old_ic','stretch_x','stretch_y','stretch_z',& & 'cyl_coord','adv_alphan','mpp_lim','hypoelasticity', & - & 'parallel_io', 'perturb_flow','perturb_sph', 'bubbles', & - & 'polytropic', 'polydisperse', 'qbmm' ] + & 'parallel_io', 'vel_profile', 'instability_wave', 'perturb_flow',& + & 'perturb_sph','bubbles', 'polytropic', 'polydisperse', 'qbmm' ] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor call MPI_BCAST(fluid_rho(1), num_fluids_max, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 16d5c45f7e..7fad3a4586 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -93,7 +93,7 @@ contains adv_alphan, mpp_lim, & weno_order, bc_x, bc_y, bc_z, num_patches, & hypoelasticity, patch_icpp, fluid_pp, & - precision, parallel_io, & + precision, parallel_io, vel_profile, instability_wave, & perturb_flow, perturb_flow_fluid, & perturb_sph, perturb_sph_fluid, fluid_rho, & cyl_coord, loops_x, loops_y, loops_z, & diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index 3be1e0cd42..f9617fbd62 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -12,7 +12,8 @@ PRE_PROCESS = COMMON + [ - 'old_grid', 'old_ic', 't_step_old', 'perturb_flow', 'perturb_flow_fluid', + 'old_grid', 'old_ic', 't_step_old', 'vel_profile', 'instability_wave', + 'perturb_flow', 'perturb_flow_fluid', 'perturb_sph', 'perturb_sph_fluid', 'fluid_rho', 'num_patches', 'qbmm', 'dist_type', 'R0_type', 'sigR', 'sigV', 'rhoRV' ] From fc6718120142a602f261bb054474fa5cc41a0e57 Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Sat, 25 Mar 2023 12:00:00 -0700 Subject: [PATCH 02/17] modify case.py for 3d_turb_mixing --- examples/3D_turb_mixing/case.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/3D_turb_mixing/case.py b/examples/3D_turb_mixing/case.py index f21044c16c..4ca8c18778 100644 --- a/examples/3D_turb_mixing/case.py +++ b/examples/3D_turb_mixing/case.py @@ -22,9 +22,9 @@ Lz = 59.0 # Number of grid cells -Nx = 127 -Ny = 127 -Nz = 127 +Nx = 255 +Ny = 255 +Nz = 255 # Grid spacing dx = Lx/float(Nx) From 4b1e401d8dbf76a829bb76601541ad00f7a8e7e5 Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Sat, 25 Mar 2023 19:39:46 -0700 Subject: [PATCH 03/17] Use global parameters for domain length instead of a hard-coded value --- src/pre_process/m_initial_condition.fpp | 29 ++++++++++++++----------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp index 559c5a840b..e2557cc419 100644 --- a/src/pre_process/m_initial_condition.fpp +++ b/src/pre_process/m_initial_condition.fpp @@ -289,40 +289,43 @@ contains subroutine s_superposition_instability_wave() ! ------------------------------ real(kind(0d0)), dimension(5,0:m,0:n,0:p) :: wave,wave1,wave2,wave_tmp real(kind(0d0)) :: tr,ti + real(kind(0d0)) :: Ly integer :: i,j,k + Ly = y_domain%end - y_domain%beg + write(*,*) "generate instability waves ..." wave = 0d0 wave1 = 0d0 wave2 = 0d0 if (p .eq. 0) then - call s_instability_wave(2*pi*4.0/59.0,0d0,tr,ti,wave_tmp,0d0) + call s_instability_wave(2*pi*4.0/Ly,0d0,tr,ti,wave_tmp,0d0) wave = wave + wave_tmp - call s_instability_wave(2*pi*2.0/59.0,0d0,tr,ti,wave_tmp,0d0) + call s_instability_wave(2*pi*2.0/Ly,0d0,tr,ti,wave_tmp,0d0) wave = wave + wave_tmp - call s_instability_wave(2*pi*1.0/59.0,0d0,tr,ti,wave_tmp,0d0) + call s_instability_wave(2*pi*1.0/Ly,0d0,tr,ti,wave_tmp,0d0) wave = wave + wave_tmp wave = wave*0.05 else - call s_instability_wave(2*pi*4.0/59.0,0d0,tr,ti,wave_tmp,0d0) + call s_instability_wave(2*pi*4.0/Ly,0d0,tr,ti,wave_tmp,0d0) wave1 = wave1 + wave_tmp - call s_instability_wave(2*pi*2.0/59.0,0d0,tr,ti,wave_tmp,0d0) + call s_instability_wave(2*pi*2.0/Ly,0d0,tr,ti,wave_tmp,0d0) wave1 = wave1 + wave_tmp - call s_instability_wave(2*pi*1.0/59.0,0d0,tr,ti,wave_tmp,0d0) + call s_instability_wave(2*pi*1.0/Ly,0d0,tr,ti,wave_tmp,0d0) wave1 = wave1 + wave_tmp - call s_instability_wave(2*pi*4.0/59.0, 2*pi*4.0/59.0,tr,ti,wave_tmp,11d0/31d0*2*pi) + call s_instability_wave(2*pi*4.0/Ly, 2*pi*4.0/Ly,tr,ti,wave_tmp,2*pi*0.8147d0) wave2 = wave2 + wave_tmp - call s_instability_wave(2*pi*2.0/59.0, 2*pi*2.0/59.0,tr,ti,wave_tmp,13d0/31d0*2*pi) + call s_instability_wave(2*pi*2.0/Ly, 2*pi*2.0/Ly,tr,ti,wave_tmp,2*pi*0.9058d0) wave2 = wave2 + wave_tmp - call s_instability_wave(2*pi*1.0/59.0, 2*pi*1.0/59.0,tr,ti,wave_tmp,17d0/31d0*2*pi) + call s_instability_wave(2*pi*1.0/Ly, 2*pi*1.0/Ly,tr,ti,wave_tmp,2*pi*0.0270d0) wave2 = wave2 + wave_tmp - call s_instability_wave(2*pi*4.0/59.0,-2*pi*4.0/59.0,tr,ti,wave_tmp,19d0/31d0*2*pi) + call s_instability_wave(2*pi*4.0/Ly,-2*pi*4.0/Ly,tr,ti,wave_tmp,2*pi*0.5324d0) wave2 = wave2 + wave_tmp - call s_instability_wave(2*pi*2.0/59.0,-2*pi*2.0/59.0,tr,ti,wave_tmp,23d0/31d0*2*pi) + call s_instability_wave(2*pi*2.0/Ly,-2*pi*2.0/Ly,tr,ti,wave_tmp,2*pi*0.2975d0) wave2 = wave2 + wave_tmp - call s_instability_wave(2*pi*1.0/59.0,-2*pi*1.0/59.0,tr,ti,wave_tmp,29d0/31d0*2*pi) + call s_instability_wave(2*pi*1.0/Ly,-2*pi*1.0/Ly,tr,ti,wave_tmp,2*pi*0.7634d0) wave2 = wave2 + wave_tmp wave = 0.05*wave1+0.15*wave2 @@ -441,7 +444,7 @@ contains ! Compute Eigenvalues & Eigenfunctions call cg(5*(n+1),5*(n+1),ar,ai,wr,wi,zr,zi,fv1,fv2,fv3,ierr) - ! Find the most unstable wave and normalize it + ! Find the most unstable wave call find_unstable_mode(5*(n+1),wr,wi,zr,zi,tr,ti,vr,vi) ! Normalize From a95e0138fc4238c73c6ac71fe871db1b6b1ae5f9 Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Sat, 25 Mar 2023 19:41:33 -0700 Subject: [PATCH 04/17] Change the number of grid points --- examples/3D_turb_mixing/case.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/3D_turb_mixing/case.py b/examples/3D_turb_mixing/case.py index 4ca8c18778..acfa7daab9 100644 --- a/examples/3D_turb_mixing/case.py +++ b/examples/3D_turb_mixing/case.py @@ -22,9 +22,9 @@ Lz = 59.0 # Number of grid cells -Nx = 255 -Ny = 255 -Nz = 255 +Nx = 319 +Ny = 319 +Nz = 319 # Grid spacing dx = Lx/float(Nx) From b13d7ca69598786be41f27fe2000414042cd4dfd Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Tue, 28 Mar 2023 22:00:38 -0700 Subject: [PATCH 05/17] Separate subroutines for solving eigenvalue problem to a module --- src/common/m_eigen_solver.f90 | 944 ++++++++++++++++++++++++++++++++++ 1 file changed, 944 insertions(+) create mode 100644 src/common/m_eigen_solver.f90 diff --git a/src/common/m_eigen_solver.f90 b/src/common/m_eigen_solver.f90 new file mode 100644 index 0000000000..84181f8c89 --- /dev/null +++ b/src/common/m_eigen_solver.f90 @@ -0,0 +1,944 @@ +!> +!! @file m_eigen_solver.f90 +!! @brief Contains module m_eigen_solver + +!> @brief The purpose of the module is to solve an eigenvalue problem +!! for a complex general matrix. Subroutines are imported +!! from EISPACK (https://netlib.org/eispack/) with minor +!! modifications for compatibility. +module m_eigen_solver + + implicit none + + private; public :: cg,cbal,corth,comqr2,csroot,cdiv,pythag + +contains + + subroutine cg(nm,nl,ar,ai,wr,wi,zr,zi,fv1,fv2,fv3,ierr) +! this subroutine calls the recommended sequence of +! subroutines from the eigensystem subroutine package (eispack) +! to find the eigenvalues and eigenvectors (if desired) +! of a complex general matrix. +! +! on input +! +! nm must be set to the row dimension of the two-dimensional +! array parameters as declared in the calling program +! dimension statement. +! +! nl is the order of the matrix a=(ar,ai). +! +! ar and ai contain the real and imaginary parts, +! respectively, of the complex general matrix. +! +! on output +! +! wr and wi contain the real and imaginary parts, +! respectively, of the eigenvalues. +! +! zr and zi contain the real and imaginary parts, +! respectively, of the eigenvectors. +! +! ierr is an integer output variable set equal to an error +! completion code described in the documentation for comqr +! and comqr2. the normal completion code is zero. +! +! fv1, fv2, and fv3 are temporary storage arrays. +! +! questions and comments should be directed to burton s. garbow, +! mathematics and computer science div, argonne national laboratory +! +! this version dated august 1983. +! +! ------------------------------------------------------------------ + integer nm,nl,is1,is2,ierr + real(kind(0d0)), dimension(nm,nl) :: ar,ai,zr,zi + real(kind(0d0)), dimension(nl) :: wr,wi,fv1,fv2,fv3 + + if (nl .le. nm) go to 10 + ierr = 10*nl + go to 50 + +10 call cbal(nm,nl,ar,ai,is1,is2,fv1) + call corth(nm,nl,is1,is2,ar,ai,fv2,fv3) + call comqr2(nm,nl,is1,is2,fv2,fv3,ar,ai,wr,wi,zr,zi,ierr) + if (ierr .ne. 0) go to 50 + call cbabk2(nm,nl,is1,is2,fv1,nl,zr,zi) + +50 return + end subroutine cg + + subroutine cbal(nm,nl,ar,ai,low,igh,scale) +! this subroutine is a translation of the algol procedure +! cbalance, which is a complex version of balance, +! num. math. 13, 293-304(1969) by parlett and reinsch. +! handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). +! +! this subroutine balances a complex matrix and isolates +! eigenvalues whenever possible. +! +! on input +! +! nm must be set to the row dimension of two-dimensional +! array parameters as declared in the calling program +! dimension statement. +! +! nl is the order of the matrix. +! +! ar and ai contain the real and imaginary parts, +! respectively, of the complex matrix to be balanced. +! +! on output +! +! ar and ai contain the real and imaginary parts, +! respectively, of the balanced matrix. +! +! low and igh are two integers such that ar(i,j) and ai(i,j) +! are equal to zero if +! (1) i is greater than j and +! (2) j=1,...,low-1 or i=igh+1,...,nl. +! +! scale contains information determining the +! permutations and scaling factors used. +! +! suppose that the principal submatrix in rows low through igh +! has been balanced, that p(j) denotes the index interchanged +! with j during the permutation step, and that the elements +! of the diagonal matrix used are denoted by d(i,j). then +! scale(j) = p(j), for j = 1,...,low-1 +! = d(j,j) j = low,...,igh +! = p(j) j = igh+1,...,nl. +! the order in which the interchanges are made is nl to igh+1, +! then 1 to low-1. +! +! note that 1 is returned for igh if igh is zero formally. +! +! the algol procedure exc contained in cbalance appears in +! cbal in line. (note that the algol roles of identifiers +! k,l have been reversed.) +! +! arithmetic is real throughout. +! +! questions and comments should be directed to burton s. garbow, +! mathematics and computer science div, argonne national laboratory +! +! this version dated august 1983. +! +! ------------------------------------------------------------------ + integer i,j,k,l,ml,nl,jj,nm,igh,low,iexc + real(kind(0d0)), dimension(nm,nl) :: ar,ai + real(kind(0d0)), dimension(nl) :: scale + real(kind(0d0)) :: c,f,g,r,s,b2,radix + logical noconv + + radix = 16.0d0 + + b2 = radix * radix + k = 1 + l = nl + go to 100 +! .......... in-line procedure for row and +! column exchange .......... +20 scale(ml) = j + if (j .eq. ml) go to 50 + + do 30 i = 1, l + f = ar(i,j) + ar(i,j) = ar(i,ml) + ar(i,ml) = f + f = ai(i,j) + ai(i,j) = ai(i,ml) + ai(i,ml) = f +30 continue + + do 40 i = k, nl + f = ar(j,i) + ar(j,i) = ar(ml,i) + ar(ml,i) = f + f = ai(j,i) + ai(j,i) = ai(ml,i) + ai(ml,i) = f +40 continue + +50 go to (80,130), iexc +! .......... search for rows isolating an eigenvalue +! and push them down .......... +80 if (l .eq. 1) go to 280 + l = l - 1 +! .......... for j=l step -1 until 1 do -- .......... +100 do 120 jj = 1, l + j = l + 1 - jj + + do 110 i = 1, l + if (i .eq. j) go to 110 + if (ar(j,i) .ne. 0.0d0 .or. ai(j,i) .ne. 0.0d0) go to 120 +110 continue + + ml = l + iexc = 1 + go to 20 +120 continue + + go to 140 +! .......... search for columns isolating an eigenvalue +! and push them left .......... +130 k = k + 1 + +140 do 170 j = k, l + + do 150 i = k, l + if (i .eq. j) go to 150 + if (ar(i,j) .ne. 0.0d0 .or. ai(i,j) .ne. 0.0d0) go to 170 +150 continue + + ml = k + iexc = 2 + go to 20 +170 continue +! .......... now balance the submatrix in rows k to l .......... + do 180 i = k, l + scale(i) = 1.0d0 +180 continue +! .......... iterative loop for norm reduction .......... +190 noconv = .false. + + do 270 i = k, l + c = 0.0d0 + r = 0.0d0 + + do 200 j = k, l + if (j .eq. i) go to 200 + c = c + dabs(ar(j,i)) + dabs(ai(j,i)) + r = r + dabs(ar(i,j)) + dabs(ai(i,j)) +200 continue +! .......... guard against zero c or r due to underflow .......... + if (c .eq. 0.0d0 .or. r .eq. 0.0d0) go to 270 + g = r / radix + f = 1.0d0 + s = c + r +210 if (c .ge. g) go to 220 + f = f * radix + c = c * b2 + go to 210 +220 g = r * radix +230 if (c .lt. g) go to 240 + f = f / radix + c = c / b2 + go to 230 +! .......... now balance .......... +240 if ((c + r) / f .ge. 0.95d0 * s) go to 270 + g = 1.0d0 / f + scale(i) = scale(i) * f + noconv = .true. + + do 250 j = k, nl + ar(i,j) = ar(i,j) * g + ai(i,j) = ai(i,j) * g +250 continue + + do 260 j = 1, l + ar(j,i) = ar(j,i) * f + ai(j,i) = ai(j,i) * f +260 continue + +270 continue + + if (noconv) go to 190 + +280 low = k + igh = l + return + end subroutine cbal + + subroutine corth(nm,nl,low,igh,ar,ai,ortr,orti) +! this subroutine is a translation of a complex analogue of +! the algol procedure orthes, num. math. 12, 349-368(1968) +! by martin and wilkinson. +! handbook for auto. comp., vol.ii-linear algebra, 339-358(1971). +! +! given a complex general matrix, this subroutine +! reduces a submatrix situated in rows and columns +! low through igh to upper hessenberg form by +! unitary similarity transformations. +! +! on input +! +! nm must be set to the row dimension of two-dimensional +! array parameters as declared in the calling program +! dimension statement. +! +! nl is the order of the matrix. +! +! low and igh are integers determined by the balancing +! subroutine cbal. if cbal has not been used, +! set low=1, igh=nl. +! +! ar and ai contain the real and imaginary parts, +! respectively, of the complex input matrix. +! +! on output +! +! ar and ai contain the real and imaginary parts, +! respectively, of the hessenberg matrix. information +! about the unitary transformations used in the reduction +! is stored in the remaining triangles under the +! hessenberg matrix. +! +! ortr and orti contain further information about the +! transformations. only elements low through igh are used. +! +! calls pythag for dsqrt(a*a + b*b) . +! +! questions and comments should be directed to burton s. garbow, +! mathematics and computer science div, argonne national laboratory +! +! this version dated august 1983. +! +! ------------------------------------------------------------------ + integer i,j,ml,nl,ii,jj,la,mp,nm,igh,kp1,low + real(kind(0d0)),dimension(nm,nl) :: ar,ai + real(kind(0d0)),dimension(igh) :: ortr,orti + real(kind(0d0)) :: f,g,h,fi,fr,scale,c + + integer mll + mll = 6 + + la = igh - 1 + kp1 = low + 1 + if (la .lt. kp1) go to 200 + + do 180 ml = kp1, la + h = 0.0d0 + ortr(ml) = 0.0d0 + orti(ml) = 0.0d0 + scale = 0.0d0 +! .......... scale column (algol tol then not needed) .......... + do 90 i = ml, igh + scale = scale + dabs(ar(i,ml-1)) + dabs(ai(i,ml-1)) +90 continue + if (scale .eq. 0d0) go to 180 + mp = ml + igh +! .......... for i=igh step -1 until ml do -- .......... + do 100 ii = ml, igh + i = mp - ii + ortr(i) = ar(i,ml-1) / scale + orti(i) = ai(i,ml-1) / scale + h = h + ortr(i) * ortr(i) + orti(i) * orti(i) +100 continue +! + g = dsqrt(h) + call pythag(ortr(ml),orti(ml),f) + if (f .eq. 0d0) go to 103 + h = h + f * g + g = g / f + ortr(ml) = (1.0d0 + g) * ortr(ml) + orti(ml) = (1.0d0 + g) * orti(ml) + go to 105 + +103 ortr(ml) = g + ar(ml,ml-1) = scale +! .......... form (i-(u*ut)/h) * a .......... +105 do 130 j = ml, nl + fr = 0.0d0 + fi = 0.0d0 +! .......... for i=igh step -1 until ml do -- .......... + do 110 ii = ml, igh + i = mp - ii + fr = fr + ortr(i) * ar(i,j) + orti(i) * ai(i,j) + fi = fi + ortr(i) * ai(i,j) - orti(i) * ar(i,j) +110 continue +! + fr = fr / h + fi = fi / h +! + do 120 i = ml, igh + ar(i,j) = ar(i,j) - fr * ortr(i) + fi * orti(i) + ai(i,j) = ai(i,j) - fr * orti(i) - fi * ortr(i) +120 continue +! +130 continue +! .......... form (i-(u*ut)/h)*a*(i-(u*ut)/h) .......... + do 160 i = 1, igh + fr = 0.0d0 + fi = 0.0d0 +! .......... for j=igh step -1 until ml do -- .......... + do 140 jj = ml, igh + j = mp - jj + fr = fr + ortr(j) * ar(i,j) - orti(j) * ai(i,j) + fi = fi + ortr(j) * ai(i,j) + orti(j) * ar(i,j) +140 continue +! + fr = fr / h + fi = fi / h +! + do 150 j = ml, igh + ar(i,j) = ar(i,j) - fr * ortr(j) - fi * orti(j) + ai(i,j) = ai(i,j) + fr * orti(j) - fi * ortr(j) +150 continue +! +160 continue +! + ortr(ml) = scale * ortr(ml) + orti(ml) = scale * orti(ml) + ar(ml,ml-1) = -g * ar(ml,ml-1) + ai(ml,ml-1) = -g * ai(ml,ml-1) +180 continue +! +200 return + end subroutine corth + + subroutine comqr2(nm,nl,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr) +! MESHED overflow control WITH vectors of isolated roots (10/19/89 BSG) +! MESHED overflow control WITH triangular multiply (10/30/89 BSG) +! +! this subroutine is a translation of a unitary analogue of the +! algol procedure comlr2, num. math. 16, 181-204(1970) by peters +! and wilkinson. +! handbook for auto. comp., vol.ii-linear algebra, 372-395(1971). +! the unitary analogue substitutes the qr algorithm of francis +! (comp. jour. 4, 332-345(1962)) for the lr algorithm. +! +! this subroutine finds the eigenvalues and eigenvectors +! of a complex upper hessenberg matrix by the qr +! method. the eigenvectors of a complex general matrix +! can also be found if corth has been used to reduce +! this general matrix to hessenberg form. +! +! on input +! +! nm must be set to the row dimension of two-dimensional +! array parameters as declared in the calling program +! dimension statement. +! +! nl is the order of the matrix. +! +! low and igh are integers determined by the balancing +! subroutine cbal. if cbal has not been used, +! set low=1, igh=nl. +! +! ortr and orti contain information about the unitary trans- +! formations used in the reduction by corth, if performed. +! only elements low through igh are used. if the eigenvectors +! of the hessenberg matrix are desired, set ortr(j) and +! orti(j) to 0.0d0 for these elements. +! +! hr and hi contain the real and imaginary parts, +! respectively, of the complex upper hessenberg matrix. +! their lower triangles below the subdiagonal contain further +! information about the transformations which were used in the +! reduction by corth, if performed. if the eigenvectors of +! the hessenberg matrix are desired, these elements may be +! arbitrary. +! +! on output +! +! ortr, orti, and the upper hessenberg portions of hr and hi +! have been destroyed. +! +! wr and wi contain the real and imaginary parts, +! respectively, of the eigenvalues. if an error +! exit is made, the eigenvalues should be correct +! for indices ierr+1,...,nl. +! +! zr and zi contain the real and imaginary parts, +! respectively, of the eigenvectors. the eigenvectors +! are unnormalized. if an error exit is made, none of +! the eigenvectors has been found. +! +! ierr is set to +! zero for normal return, +! j if the limit of 30*nl iterations is exhausted +! while the j-th eigenvalue is being sought. +! +! calls cdiv for complex division. +! calls csroot for complex square root. +! calls pythag for dsqrt(a*a + b*b) . +! +! questions and comments should be directed to burton s. garbow, +! mathematics and computer science div, argonne national laboratory +! +! this version dated october 1989. +! +! ------------------------------------------------------------------ + integer i,j,k,l,ml,nl,en,ii,jj,ll,nm,nn,igh,ip1,& + itn,its,low,lp1,enm1,iend,ierr + real(kind(0d0)),dimension(nm,nl) :: hr,hi,zr,zi + real(kind(0d0)),dimension(nl) :: wr,wi + real(kind(0d0)),dimension(igh) :: ortr,orti + real(kind(0d0)) :: si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr,& + norm,tst1,tst2,c,d +! + ierr = 0 +! .......... initialize eigenvector matrix .......... + do 101 j = 1, nl +! + do 100 i = 1, nl + zr(i,j) = 0.0d0 + zi(i,j) = 0.0d0 +100 continue + zr(j,j) = 1.0d0 +101 continue +! .......... form the matrix of accumulated transformations +! from the information left by corth .......... + iend = igh - low - 1 + if (iend .lt. 0) go to 180 + if (iend .eq. 0) go to 150 + if (iend .gt. 0) go to 105 +! .......... for i=igh-1 step -1 until low+1 do -- .......... +105 do 140 ii = 1, iend + i = igh - ii + if (dabs(ortr(i)) .eq. 0d0 .and. dabs(orti(i)) .eq. 0d0) go to 140 + if (dabs(hr(i,i-1)) .eq. 0d0 .and. dabs(hi(i,i-1)) .eq. 0d0) go to 140 +! .......... norm below is negative of h formed in corth .......... + norm = hr(i,i-1) * ortr(i) + hi(i,i-1) * orti(i) + ip1 = i + 1 + + do 110 k = ip1, igh + ortr(k) = hr(k,i-1) + orti(k) = hi(k,i-1) +110 continue +! + do 130 j = i, igh + sr = 0.0d0 + si = 0.0d0 +! + do 115 k = i, igh + sr = sr + ortr(k) * zr(k,j) + orti(k) * zi(k,j) + si = si + ortr(k) * zi(k,j) - orti(k) * zr(k,j) +115 continue +! + sr = sr / norm + si = si / norm +! + do 120 k = i, igh + zr(k,j) = zr(k,j) + sr * ortr(k) - si * orti(k) + zi(k,j) = zi(k,j) + sr * orti(k) + si * ortr(k) +120 continue +! +130 continue +! +140 continue +! .......... create real subdiagonal elements .......... +150 l = low + 1 +! + do 170 i = l, igh + ll = min0(i+1,igh) + if (dabs(hi(i,i-1)) .eq. 0d0) go to 170 + call pythag(hr(i,i-1),hi(i,i-1),norm) + yr = hr(i,i-1) / norm + yi = hi(i,i-1) / norm + hr(i,i-1) = norm + hi(i,i-1) = 0.0d0 +! + do 155 j = i, nl + si = yr * hi(i,j) - yi * hr(i,j) + hr(i,j) = yr * hr(i,j) + yi * hi(i,j) + hi(i,j) = si +155 continue +! + do 160 j = 1, ll + si = yr * hi(j,i) + yi * hr(j,i) + hr(j,i) = yr * hr(j,i) - yi * hi(j,i) + hi(j,i) = si +160 continue +! + do 165 j = low, igh + si = yr * zi(j,i) + yi * zr(j,i) + zr(j,i) = yr * zr(j,i) - yi * zi(j,i) + zi(j,i) = si +165 continue + +170 continue +! .......... store roots isolated by cbal .......... +180 do 200 i = 1, nl + if (i .ge. low .and. i .le. igh) go to 200 + wr(i) = hr(i,i) + wi(i) = hi(i,i) +200 continue +! + en = igh + tr = 0.0d0 + ti = 0.0d0 + itn = 30*nl +! .......... search for next eigenvalue .......... +220 if (en .lt. low) go to 680 + its = 0 + enm1 = en - 1 +! .......... look for single small sub-diagonal element +! for l=en step -1 until low do -- .......... +240 do 260 ll = low, en + l = en + low - ll + if (l .eq. low) go to 300 + tst1 = dabs(hr(l-1,l-1)) + dabs(hi(l-1,l-1)) & + + dabs(hr(l,l)) + dabs(hi(l,l)) + tst2 = tst1 + dabs(hr(l,l-1)) + if (tst2 .eq. tst1) go to 300 +260 continue +! .......... form shift .......... +300 if (l .eq. en) go to 660 + if (itn .eq. 0) go to 1000 + if (its .eq. 10 .or. its .eq. 20) go to 320 + sr = hr(en,en) + si = hi(en,en) + xr = hr(enm1,en) * hr(en,enm1) + xi = hi(enm1,en) * hr(en,enm1) + if (xr .eq. 0.0d0 .and. xi .eq. 0.0d0) go to 340 + yr = (hr(enm1,enm1) - sr) / 2.0d0 + yi = (hi(enm1,enm1) - si) / 2.0d0 + call csroot(yr**2-yi**2+xr,2.0d0*yr*yi+xi,zzr,zzi) + if (yr * zzr + yi * zzi .ge. 0.0d0) go to 310 + zzr = -zzr + zzi = -zzi +310 call cdiv(xr,xi,yr+zzr,yi+zzi,xr,xi) + sr = sr - xr + si = si - xi + go to 340 +! .......... form exceptional shift .......... +320 sr = dabs(hr(en,enm1)) + dabs(hr(enm1,en-2)) + si = 0.0d0 +! +340 do 360 i = low, en + hr(i,i) = hr(i,i) - sr + hi(i,i) = hi(i,i) - si +360 continue +! + tr = tr + sr + ti = ti + si + its = its + 1 + itn = itn - 1 +! .......... reduce to triangle (rows) .......... + lp1 = l + 1 +! + do 500 i = lp1, en + sr = hr(i,i-1) + hr(i,i-1) = 0.0d0 + call pythag(hr(i-1,i-1),hi(i-1,i-1),c) + call pythag(c,sr,norm) + xr = hr(i-1,i-1) / norm + wr(i-1) = xr + xi = hi(i-1,i-1) / norm + wi(i-1) = xi + hr(i-1,i-1) = norm + hi(i-1,i-1) = 0.0d0 + hi(i,i-1) = sr / norm +! + do 490 j = i, nl + yr = hr(i-1,j) + yi = hi(i-1,j) + zzr = hr(i,j) + zzi = hi(i,j) + hr(i-1,j) = xr * yr + xi * yi + hi(i,i-1) * zzr + hi(i-1,j) = xr * yi - xi * yr + hi(i,i-1) * zzi + hr(i,j) = xr * zzr - xi * zzi - hi(i,i-1) * yr + hi(i,j) = xr * zzi + xi * zzr - hi(i,i-1) * yi +490 continue +! +500 continue +! + si = hi(en,en) + if (dabs(si) .eq. 0d0) go to 540 + call pythag(hr(en,en),si,norm) + sr = hr(en,en) / norm + si = si / norm + hr(en,en) = norm + hi(en,en) = 0.0d0 + if (en .eq. nl) go to 540 + ip1 = en + 1 +! + do 520 j = ip1, nl + yr = hr(en,j) + yi = hi(en,j) + hr(en,j) = sr * yr + si * yi + hi(en,j) = sr * yi - si * yr +520 continue +! .......... inverse operation (columns) .......... +540 do 600 j = lp1, en + xr = wr(j-1) + xi = wi(j-1) +! + do 580 i = 1, j + yr = hr(i,j-1) + yi = 0.0d0 + zzr = hr(i,j) + zzi = hi(i,j) + if (i .eq. j) go to 560 + yi = hi(i,j-1) + hi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi +560 hr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr + hr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr + hi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi +580 continue +! + do 590 i = low, igh + yr = zr(i,j-1) + yi = zi(i,j-1) + zzr = zr(i,j) + zzi = zi(i,j) + zr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr + zi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi + zr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr + zi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi +590 continue + +600 continue +! + if (dabs(si) .eq. 0d0) go to 240 +! + do 630 i = 1, en + yr = hr(i,en) + yi = hi(i,en) + hr(i,en) = sr * yr - si * yi + hi(i,en) = sr * yi + si * yr +630 continue +! + do 640 i = low, igh + yr = zr(i,en) + yi = zi(i,en) + zr(i,en) = sr * yr - si * yi + zi(i,en) = sr * yi + si * yr +640 continue +! + go to 240 +! .......... a root found .......... +660 hr(en,en) = hr(en,en) + tr + wr(en) = hr(en,en) + hi(en,en) = hi(en,en) + ti + wi(en) = hi(en,en) + en = enm1 + go to 220 +! .......... all roots found. backsubstitute to find +! vectors of upper triangular form .......... +680 norm = 0.0d0 +! + do i = 1, nl + do j = i, nl + tr = dabs(hr(i,j)) + dabs(hi(i,j)) + if (tr .gt. norm) norm = tr + end do + end do +! + if (nl .eq. 1 .or. norm .eq. 0d0) go to 1001 +! .......... for en=nl step -1 until 2 do -- .......... + do 800 nn = 2, nl + en = nl + 2 - nn + xr = wr(en) + xi = wi(en) + hr(en,en) = 1.0d0 + hi(en,en) = 0.0d0 + enm1 = en - 1 +! .......... for i=en-1 step -1 until 1 do -- .......... + do 780 ii = 1, enm1 + i = en - ii + zzr = 0.0d0 + zzi = 0.0d0 + ip1 = i + 1 + + do 740 j = ip1, en + zzr = zzr + hr(i,j) * hr(j,en) - hi(i,j) * hi(j,en) + zzi = zzi + hr(i,j) * hi(j,en) + hi(i,j) * hr(j,en) +740 continue +! + yr = xr - wr(i) + yi = xi - wi(i) + if (yr .ne. 0.0d0 .or. yi .ne. 0.0d0) go to 765 + tst1 = norm + yr = tst1 +760 yr = 0.01d0 * yr + tst2 = norm + yr + if (tst2 .gt. tst1) go to 760 +765 continue + call cdiv(zzr,zzi,yr,yi,hr(i,en),hi(i,en)) +! .......... overflow control .......... + tr = dabs(hr(i,en)) + dabs(hi(i,en)) + if (tr .eq. 0.0d0) go to 780 + tst1 = tr + tst2 = tst1 + 1.0d0/tst1 + if (tst2 .gt. tst1) go to 780 + do 770 j = i, en + hr(j,en) = hr(j,en)/tr + hi(j,en) = hi(j,en)/tr +770 continue +! +780 continue +! +800 continue +! .......... end backsubstitution .......... +! .......... vectors of isolated roots .......... + do 840 i = 1, nl + if (i .ge. low .and. i .le. igh) go to 840 +! + do 820 j = I, nl + zr(i,j) = hr(i,j) + zi(i,j) = hi(i,j) +820 continue +! +840 continue +! .......... multiply by transformation matrix to give +! vectors of original full matrix. +! for j=nl step -1 until low do -- .......... + do jj = low, nl + j = nl + low - jj + ml = min0(j,igh) +! + do i = low, igh + zzr = 0.0d0 + zzi = 0.0d0 +! + do 860 k = low, ml + zzr = zzr + zr(i,k) * hr(k,j) - zi(i,k) * hi(k,j) + zzi = zzi + zr(i,k) * hi(k,j) + zi(i,k) * hr(k,j) +860 continue +! + zr(i,j) = zzr + zi(i,j) = zzi + end do + end do +! + go to 1001 +! .......... set error -- all eigenvalues have not +! converged after 30*nl iterations .......... +1000 ierr = en +1001 return + end subroutine comqr2 + + subroutine cbabk2(nm,nl,low,igh,scale,ml,zr,zi) +! this subroutine is a translation of the algol procedure +! cbabk2, which is a complex version of balbak, +! num. math. 13, 293-304(1969) by parlett and reinsch. +! handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). +! +! this subroutine forms the eigenvectors of a complex general +! matrix by back transforming those of the corresponding +! balanced matrix determined by cbal. +! +! on input +! +! nm must be set to the row dimension of two-dimensional +! array parameters as declared in the calling program +! dimension statement. +! +! nl is the order of the matrix. +! +! low and igh are integers determined by cbal. +! +! scale contains information determining the permutations +! and scaling factors used by cbal. +! +! ml is the number of eigenvectors to be back transformed. +! +! zr and zi contain the real and imaginary parts, +! respectively, of the eigenvectors to be +! back transformed in their first ml columns. +! +! on output +! +! zr and zi contain the real and imaginary parts, +! respectively, of the transformed eigenvectors +! in their first ml columns. +! +! questions and comments should be directed to burton s. garbow, +! mathematics and computer science div, argonne national laboratory +! +! this version dated august 1983. +! +! ------------------------------------------------------------------ +! + integer i,j,k,ml,nl,ii,nm,igh,low + double precision scale(nl),zr(nm,ml),zi(nm,ml) + double precision s + + if (ml .eq. 0) go to 200 + if (igh .eq. low) go to 120 +! + do 110 i = low, igh + s = scale(i) +! .......... left hand eigenvectors are back transformed +! if the foregoing statement is replaced by +! s=1.0d0/scale(i). .......... + do 100 j = 1, ml + zr(i,j) = zr(i,j) * s + zi(i,j) = zi(i,j) * s +100 continue +! +110 continue +! .......... for i=low-1 step -1 until 1, +! igh+1 step 1 until nl do -- .......... +120 do 140 ii = 1, nl + i = ii + if (i .ge. low .and. i .le. igh) go to 140 + if (i .lt. low) i = low - ii + k = scale(i) + if (k .eq. i) go to 140 +! + do 130 j = 1, ml + s = zr(i,j) + zr(i,j) = zr(k,j) + zr(k,j) = s + s = zi(i,j) + zi(i,j) = zi(k,j) + zi(k,j) = s +130 continue +! +140 continue +! +200 return + end subroutine cbabk2 + + subroutine csroot(xr,xi,yr,yi) + real(kind(0d0)) :: xr,xi,yr,yi +! +! (yr,yi) = complex dsqrt(xr,xi) +! branch chosen so that yr .ge. 0.0 and sign(yi) .eq. sign(xi) +! + real(kind(0d0)) :: s,tr,ti,c + tr = xr + ti = xi + call pythag(tr,ti,c) + s = dsqrt(0.5d0*(c + dabs(tr))) + if (tr .ge. 0.0d0) yr = s + if (ti .lt. 0.0d0) s = -s + if (tr .le. 0.0d0) yi = s + if (tr .lt. 0.0d0) yr = 0.5d0*(ti/yi) + if (tr .gt. 0.0d0) yi = 0.5d0*(ti/yr) + return + end subroutine csroot + + subroutine cdiv(ar,ai,br,bi,cr,ci) + real(kind(0d0)) :: ar,ai,br,bi,cr,ci +! +! complex division, (cr,ci) = (ar,ai)/(br,bi) +! + real(kind(0d0)) :: s,ars,ais,brs,bis + s = dabs(br) + dabs(bi) + ars = ar/s + ais = ai/s + brs = br/s + bis = bi/s + s = brs**2 + bis**2 + cr = (ars*brs + ais*bis)/s + ci = (ais*brs - ars*bis)/s + return + end subroutine cdiv + + subroutine pythag(a,b,c) + real(kind(0d0)) :: a,b,c +! +! finds dsqrt(a**2+b**2) without overflow or destructive underflow +! + real(kind(0d0)) :: p,r,s,t,u + p = dmax1(dabs(a),dabs(b)) + if (p .eq. 0.0d0) go to 20 + r = (dmin1(dabs(a),dabs(b))/p)**2 +10 continue + t = 4.0d0 + r + if (t .eq. 4.0d0) go to 20 + s = r/t + u = 1.0d0 + 2.0d0*s + p = u*p + r = (s/u)**2 * r + go to 10 +20 c = p + return + end subroutine pythag + +end module m_eigen_solver From f73d53b953cc8f63e9f33b2a42b099dd5af0fa38 Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Tue, 28 Mar 2023 22:07:17 -0700 Subject: [PATCH 06/17] Add documentation and comments, make indent consistent, add checker to prevent the use of vel_profile and instability_wave in 1D case --- src/pre_process/m_assign_variables.f90 | 16 +- src/pre_process/m_checker.f90 | 12 + src/pre_process/m_global_parameters.fpp | 11 +- src/pre_process/m_initial_condition.fpp | 1411 ++++------------------- 4 files changed, 254 insertions(+), 1196 deletions(-) diff --git a/src/pre_process/m_assign_variables.f90 b/src/pre_process/m_assign_variables.f90 index 9c26feef94..0b36a13f07 100644 --- a/src/pre_process/m_assign_variables.f90 +++ b/src/pre_process/m_assign_variables.f90 @@ -485,6 +485,13 @@ subroutine s_assign_patch_species_primitive_variables_bubbles(patch_id, j, k, l, + (1d0 - eta)*orig_prim_vf(i + cont_idx%end)) end do + ! Set streamwise velocity to hypertangent function of y + if (vel_profile) then + q_prim_vf(1 + cont_idx%end)%sf(j, k, l) = & + (eta*patch_icpp(patch_id)%vel(1)*tanh(y_cc(k)) & + + (1d0 - eta)*orig_prim_vf(1 + cont_idx%end)) + end if + ! Smoothed bubble variables if (bubbles) then do i = 1, nb @@ -679,11 +686,12 @@ subroutine s_assign_patch_species_primitive_variables(patch_id, j, k, l, & + (1d0 - eta)*orig_prim_vf(i + cont_idx%end)) end do - if (vel_profile) then - q_prim_vf(1 + cont_idx%end)%sf(j, k, l) = & + ! Set streamwise velocity to hypertangent function of y + if (vel_profile) then + q_prim_vf(1 + cont_idx%end)%sf(j, k, l) = & (eta*patch_icpp(patch_id)%vel(1)*tanh(y_cc(k)) & - + (1d0 - eta)*orig_prim_vf(1 + cont_idx%end)) - end if + + (1d0 - eta)*orig_prim_vf(1 + cont_idx%end)) + end if ! Pressure q_prim_vf(E_idx)%sf(j, k, l) = & diff --git a/src/pre_process/m_checker.f90 b/src/pre_process/m_checker.f90 index fb67f2dd56..9102781686 100644 --- a/src/pre_process/m_checker.f90 +++ b/src/pre_process/m_checker.f90 @@ -565,6 +565,18 @@ subroutine s_check_inputs() end do end if + ! Constraints on the hypertangent velocity profile + if ((vel_profile .eqv. .true.) .and. (n == 0)) then + call s_mpi_abort('Unsupported choices of the combination of values for '//& + 'vel_profile and n. Exiting ...') + end if + + ! Constraints on the instability wave + if ((instability_wave .eqv. .true.) .and. (n == 0)) then + call s_mpi_abort('Unsupported choices of the combination of values for '//& + 'instability_wave and n. Exiting ...') + end if + ! Constraints on the stiffened equation of state fluids parameters do i = 1, num_fluids call s_int_to_str(i, iStr) diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 903a1d9877..31c33deb15 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -97,11 +97,9 @@ module m_global_parameters logical :: parallel_io !< Format of the data files integer :: precision !< Precision of output files - ! Hypertangent velocity profile - logical :: vel_profile + logical :: vel_profile !< Set hypertangent streamwise velocity profile - ! Add instability waves to surrounding fluid - logical :: instability_wave + logical :: instability_wave !< Superimpose instability waves to surrounding fluid flow ! Perturb density of surrounding air so as to break symmetry of grid logical :: perturb_flow @@ -250,8 +248,8 @@ contains parallel_io = .false. precision = 2 - vel_profile = .false. - instability_wave = .false. + vel_profile = .false. + instability_wave = .false. perturb_flow = .false. perturb_flow_fluid = dflt_int perturb_sph = .false. @@ -293,6 +291,7 @@ contains patch_icpp(i)%p0 = dflt_real patch_icpp(i)%m0 = dflt_real + end do ! Tait EOS diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp index e2557cc419..b4426073ee 100644 --- a/src/pre_process/m_initial_condition.fpp +++ b/src/pre_process/m_initial_condition.fpp @@ -27,6 +27,10 @@ module m_initial_condition use m_patches use m_assign_variables + + use m_eigen_solver ! Subroutines to solve eigenvalue problem for + ! complex general matrix + ! ========================================================================== ! ========================================================================== @@ -211,7 +215,6 @@ contains if (perturb_flow) call s_perturb_surrounding_flow() if (perturb_sph) call s_perturb_sphere() - if (instability_wave) call s_superposition_instability_wave() ! Converting the primitive variables to the conservative ones @@ -285,1205 +288,241 @@ contains end do end subroutine s_perturb_surrounding_flow ! ---------------------------- - - subroutine s_superposition_instability_wave() ! ------------------------------ - real(kind(0d0)), dimension(5,0:m,0:n,0:p) :: wave,wave1,wave2,wave_tmp - real(kind(0d0)) :: tr,ti - real(kind(0d0)) :: Ly - integer :: i,j,k - Ly = y_domain%end - y_domain%beg + !> This subroutine computes velocity perturbations for a temporal mixing + !! layer with hypertangent mean streamwise velocity profile + !! obtained from linear stability analysis. For a 2D case, + !! instability waves with spatial wavenumbers, (4,0), (2,0), + !! and (1,0) are superposed. For a 3D waves, (4,4), (4,-4), + !! (2,2), (2,-2), (1,1), (1,-1) areadded on top of 2D waves. + subroutine s_superposition_instability_wave() ! ------------------------ + real(kind(0d0)), dimension(5,0:m,0:n,0:p) :: wave,wave1,wave2,wave_tmp + real(kind(0d0)) :: tr,ti + real(kind(0d0)) :: Ly + integer :: i,j,k + + Ly = y_domain%end - y_domain%beg + + wave = 0d0 + wave1 = 0d0 + wave2 = 0d0 + + ! Compute 2D waves + call s_instability_wave(2*pi*4.0/Ly,0d0,tr,ti,wave_tmp,0d0) + wave1 = wave1 + wave_tmp + call s_instability_wave(2*pi*2.0/Ly,0d0,tr,ti,wave_tmp,0d0) + wave1 = wave1 + wave_tmp + call s_instability_wave(2*pi*1.0/Ly,0d0,tr,ti,wave_tmp,0d0) + wave1 = wave1 + wave_tmp + wave = wave1*0.05 - write(*,*) "generate instability waves ..." - - wave = 0d0 - wave1 = 0d0 - wave2 = 0d0 - if (p .eq. 0) then - call s_instability_wave(2*pi*4.0/Ly,0d0,tr,ti,wave_tmp,0d0) - wave = wave + wave_tmp - call s_instability_wave(2*pi*2.0/Ly,0d0,tr,ti,wave_tmp,0d0) - wave = wave + wave_tmp - call s_instability_wave(2*pi*1.0/Ly,0d0,tr,ti,wave_tmp,0d0) - wave = wave + wave_tmp - wave = wave*0.05 - else - call s_instability_wave(2*pi*4.0/Ly,0d0,tr,ti,wave_tmp,0d0) - wave1 = wave1 + wave_tmp - call s_instability_wave(2*pi*2.0/Ly,0d0,tr,ti,wave_tmp,0d0) - wave1 = wave1 + wave_tmp - call s_instability_wave(2*pi*1.0/Ly,0d0,tr,ti,wave_tmp,0d0) - wave1 = wave1 + wave_tmp - - call s_instability_wave(2*pi*4.0/Ly, 2*pi*4.0/Ly,tr,ti,wave_tmp,2*pi*0.8147d0) - wave2 = wave2 + wave_tmp - call s_instability_wave(2*pi*2.0/Ly, 2*pi*2.0/Ly,tr,ti,wave_tmp,2*pi*0.9058d0) - wave2 = wave2 + wave_tmp - call s_instability_wave(2*pi*1.0/Ly, 2*pi*1.0/Ly,tr,ti,wave_tmp,2*pi*0.0270d0) - wave2 = wave2 + wave_tmp - call s_instability_wave(2*pi*4.0/Ly,-2*pi*4.0/Ly,tr,ti,wave_tmp,2*pi*0.5324d0) - wave2 = wave2 + wave_tmp - call s_instability_wave(2*pi*2.0/Ly,-2*pi*2.0/Ly,tr,ti,wave_tmp,2*pi*0.2975d0) - wave2 = wave2 + wave_tmp - call s_instability_wave(2*pi*1.0/Ly,-2*pi*1.0/Ly,tr,ti,wave_tmp,2*pi*0.7634d0) - wave2 = wave2 + wave_tmp - - wave = 0.05*wave1+0.15*wave2 - end if - + if (p > 0) then + ! Compute 3D waves with phase shifts. + call s_instability_wave(2*pi*4.0/Ly, 2*pi*4.0/Ly,tr,ti,wave_tmp,2*pi*11d0/31d0) + wave2 = wave2 + wave_tmp + call s_instability_wave(2*pi*2.0/Ly, 2*pi*2.0/Ly,tr,ti,wave_tmp,2*pi*13d0/31d0) + wave2 = wave2 + wave_tmp + call s_instability_wave(2*pi*1.0/Ly, 2*pi*1.0/Ly,tr,ti,wave_tmp,2*pi*17d0/31d0) + wave2 = wave2 + wave_tmp + call s_instability_wave(2*pi*4.0/Ly,-2*pi*4.0/Ly,tr,ti,wave_tmp,2*pi*19d0/31d0) + wave2 = wave2 + wave_tmp + call s_instability_wave(2*pi*2.0/Ly,-2*pi*2.0/Ly,tr,ti,wave_tmp,2*pi*23d0/31d0) + wave2 = wave2 + wave_tmp + call s_instability_wave(2*pi*1.0/Ly,-2*pi*1.0/Ly,tr,ti,wave_tmp,2*pi*29d0/31d0) + wave2 = wave2 + wave_tmp + wave = wave + 0.15*wave2 + end if + + ! Superpose velocity perturbuations (instability waves) to the velocity field do k = 0, p do j = 0, n do i = 0, m - q_prim_vf(mom_idx%beg )%sf(i,j,k) = q_prim_vf(mom_idx%beg )%sf(i,j,k)+wave(2,i,j,k) - q_prim_vf(mom_idx%beg+1)%sf(i,j,k) = q_prim_vf(mom_idx%beg+1)%sf(i,j,k)+wave(3,i,j,k) - if (p .gt. 0) then - q_prim_vf(mom_idx%beg+2)%sf(i,j,k) = q_prim_vf(mom_idx%beg+2)%sf(i,j,k)+wave(4,i,j,k) - end if - end do - end do - end do - - end subroutine s_superposition_instability_wave - + q_prim_vf(mom_idx%beg )%sf(i,j,k) = q_prim_vf(mom_idx%beg )%sf(i,j,k)+wave(2,i,j,k) + q_prim_vf(mom_idx%beg+1)%sf(i,j,k) = q_prim_vf(mom_idx%beg+1)%sf(i,j,k)+wave(3,i,j,k) + if (p .gt. 0) then + q_prim_vf(mom_idx%beg+2)%sf(i,j,k) = q_prim_vf(mom_idx%beg+2)%sf(i,j,k)+wave(4,i,j,k) + end if + end do + end do + end do + + end subroutine s_superposition_instability_wave ! ---------------------- + + !> This subroutine computes instability waves corresponding to the spatial + !! wavenumbers (alpha, beta) for x and z directions, + !! respectively. The eigenvalue problem is derived from the + !! linearized Euler equations with parallel mean flow + !! assumption (See Sandham 1989 PhD thesis for details). subroutine s_instability_wave(alpha,beta,tr,ti,wave,shift) - real(kind(0d0)),intent(in) :: alpha, beta !< spatial wavenumbers - real(kind(0d0)),dimension(0:n) :: rho_mean, u_mean, t_mean !< mean profiles - real(kind(0d0)),dimension(0:n) :: drho_mean, du_mean, dt_mean !< mean profiles - real(kind(0d0)),dimension(0:n,0:n) :: d - real(kind(0d0)),dimension(0:5*(n+1)-1,0:5*(n+1)-1) :: ar,ai,br,bi,ci - real(kind(0d0)),dimension(0:5*(n+1)-1,0:5*(n+1)-1) :: zr,zi - real(kind(0d0)),dimension(0:5*(n+1)-1) :: wr,wi,fv1,fv2,fv3 - real(kind(0d0)) :: tr,ti - real(kind(0d0)),dimension(0:5*(n+1)-1) :: vr,vi,vnr,vni - real(kind(0d0)),dimension(5,0:m,0:n,0:p) :: wave - real(kind(0d0)) :: shift + real(kind(0d0)),intent(in) :: alpha, beta !< spatial wavenumbers + real(kind(0d0)),dimension(0:n) :: rho_mean, u_mean, t_mean !< mean profiles + real(kind(0d0)),dimension(0:n) :: drho_mean, du_mean, dt_mean !< y-derivatives of mean profiles + real(kind(0d0)),dimension(0:n,0:n) :: d !< differential operator in y dir + real(kind(0d0)),dimension(0:5*(n+1)-1,0:5*(n+1)-1) :: ar,ai,br,bi,ci !< matrices for eigenvalue problem + real(kind(0d0)),dimension(0:5*(n+1)-1,0:5*(n+1)-1) :: zr,zi !< eigenvectors + real(kind(0d0)),dimension(0:5*(n+1)-1) :: wr,wi !< eigenvalues + real(kind(0d0)),dimension(0:5*(n+1)-1) :: fv1,fv2,fv3 !< temporary memory + real(kind(0d0)) :: tr,ti !< most unstable eigenvalue + real(kind(0d0)),dimension(0:5*(n+1)-1) :: vr,vi,vnr,vni !< most unstable eigenvector and normalized one + real(kind(0d0)),dimension(5,0:m,0:n,0:p) :: wave !< instability wave + real(kind(0d0)) :: shift !< phase shift + real(kind(0d0)) :: gam,pi_inf,rho1,mach,c1 + integer :: ierr + integer :: j, k, l !< generic loop iterators + integer :: ii, jj !< block matrix indicies + + ! Set fluid flow properties + gam = 1.+1./fluid_pp(1)%gamma + pi_inf = fluid_pp(1)%pi_inf*(gam-1.)/gam + rho1 = patch_icpp(1)%alpha_rho(1)/patch_icpp(1)%alpha(1) + c1 = sqrt((gam*(patch_icpp(1)%pres+pi_inf))/rho1) + mach = 1./c1 + + ! Assign mean profiles + do j=0,n + u_mean(j)=tanh(y_cc(j)) + t_mean(j)=1+0.5*(gam-1)*mach**2*(1-u_mean(j)**2) + rho_mean(j)=1/T_mean(j) + end do + + ! Compute differential operator in y-dir + ! based on 4th order central difference (inner) + ! and 2nd order central difference (near boundaries) + dy = y_cc(1)-y_cc(0) + d=0d0 + d(1,0)=-1/(2*dy) + d(1,2)= 1/(2*dy) + do j=2,n-2 + d(j,j-2)= 1/(12*dy) + d(j,j-1)=-8/(12*dy) + d(j,j+1)= 8/(12*dy) + d(j,j+2)=-1/(12*dy) + end do + d(n-1,n-2)=-1/(2*dy) + d(n-1,n) = 1/(2*dy) + + ! Compute y-derivatives of rho, u, T + do j=0,n + drho_mean(j)=0 + du_mean(j)=0 + dt_mean(j)=0 + do k=0,n + drho_mean(j) = drho_mean(j)+d(j,k)*rho_mean(k) + du_mean(j) = du_mean(j)+d(j,k)*u_mean(k) + dt_mean(j) = dt_mean(j)+d(j,k)*t_mean(k) + end do + end do + + ! Compute B and C, then A = B + C + ! B includes terms without differential operator, and + ! C includes terms with differential operator + br=0d0 + bi=0d0 + ci=0d0 + do j=0,n + ii = 1; jj = 1; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); + ii = 1; jj = 2; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*rho_mean(j); + ii = 1; jj = 3; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -drho_mean(j); + ii = 1; jj = 4; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = beta*rho_mean(j); + + ii = 2; jj = 1; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*t_mean(j)/(rho_mean(j)*gam*mach**2); + ii = 2; jj = 2; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); + ii = 2; jj = 3; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -du_mean(j); + ii = 2; jj = 5; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha/(gam*mach**2); + + ii = 3; jj = 1; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -dt_mean(j)/(rho_mean(j)*gam*mach**2); + ii = 3; jj = 3; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); + ii = 3; jj = 5; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -drho_mean(j)/(rho_mean(j)*gam*mach**2); + + ii = 4; jj = 1; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = beta*t_mean(j)/(rho_mean(j)*gam*mach**2); + ii = 4; jj = 4; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); + ii = 4; jj = 5; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = beta/(gam*mach**2); + + ii = 5; jj = 2; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = (gam-1)*alpha/rho_mean(j); + ii = 5; jj = 3; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -dt_mean(j); + ii = 5; jj = 4; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = (gam-1)*beta/rho_mean(j); + ii = 5; jj = 5; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); + + do k=0,n + ii = 1; jj = 3; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -rho_mean(j)*d(j,k); + ii = 3; jj = 1; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -t_mean(j)*d(j,k)/(rho_mean(j)*gam*mach**2); + ii = 3; jj = 5; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -d(j,k)/(gam*mach**2); + ii = 5; jj = 3; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -(gam-1)*d(j,k)/rho_mean(j); + end do + end do + ar = br + ai = bi+ci + + ! Compute eigenvalues and eigenvectors + call cg(5*(n+1),5*(n+1),ar,ai,wr,wi,zr,zi,fv1,fv2,fv3,ierr) - integer :: ierr - integer :: j, k, l !< generic loop iterators - integer :: ii, jj !< block matrix indicies - - real(kind(0d0)) :: gam,pi_inf,rho1,mach,c1 - - gam = 1.+1./fluid_pp(1)%gamma - pi_inf = fluid_pp(1)%pi_inf*(gam-1.)/gam - rho1 = patch_icpp(1)%alpha_rho(1)/patch_icpp(1)%alpha(1) - c1 = sqrt((gam*(patch_icpp(1)%pres+pi_inf))/rho1) - mach = 1./c1 - - ! Assign mean profiles - do j=0,n - u_mean(j)=tanh(y_cc(j)) - t_mean(j)=1+0.5*(gam-1)*mach**2*(1-u_mean(j)**2) - rho_mean(j)=1/T_mean(j) - end do - - ! Compute differential operator - dy = y_cc(1)-y_cc(0) - d=0d0 - d(1,0)=-1/(2*dy) - d(1,2)= 1/(2*dy) - do j=2,n-2 - d(j,j-2)= 1/(12*dy) - d(j,j-1)=-8/(12*dy) - d(j,j+1)= 8/(12*dy) - d(j,j+2)=-1/(12*dy) - end do - d(n-1,n-2)=-1/(2*dy) - d(n-1,n) = 1/(2*dy) - - ! Compute y-derivatives of rho, u, T - do j=0,n - drho_mean(j)=0 - du_mean(j)=0 - dt_mean(j)=0 - do k=0,n - drho_mean(j) = drho_mean(j)+d(j,k)*rho_mean(k) - du_mean(j) = du_mean(j)+d(j,k)*u_mean(k) - dt_mean(j) = dt_mean(j)+d(j,k)*t_mean(k) - end do - end do - - ! Compute B and C -> A=B+C - br=0d0 - bi=0d0 - ci=0d0 - do j=0,n - ii = 1; jj = 1; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); - ii = 1; jj = 2; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*rho_mean(j); - ii = 1; jj = 3; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -drho_mean(j); - ii = 1; jj = 4; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = beta*rho_mean(j); - - ii = 2; jj = 1; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*t_mean(j)/(rho_mean(j)*gam*mach**2); - ii = 2; jj = 2; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); - ii = 2; jj = 3; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -du_mean(j); - ii = 2; jj = 5; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha/(gam*mach**2); - - ii = 3; jj = 1; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -dt_mean(j)/(rho_mean(j)*gam*mach**2); - ii = 3; jj = 3; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); - ii = 3; jj = 5; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -drho_mean(j)/(rho_mean(j)*gam*mach**2); - - ii = 4; jj = 1; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = beta*t_mean(j)/(rho_mean(j)*gam*mach**2); - ii = 4; jj = 4; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); - ii = 4; jj = 5; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = beta/(gam*mach**2); - - ii = 5; jj = 2; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = (gam-1)*alpha/rho_mean(j); - ii = 5; jj = 3; bi((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = -dt_mean(j); - ii = 5; jj = 4; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = (gam-1)*beta/rho_mean(j); - ii = 5; jj = 5; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); - - do k=0,n - ii = 1; jj = 3; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -rho_mean(j)*d(j,k); - ii = 3; jj = 1; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -t_mean(j)*d(j,k)/(rho_mean(j)*gam*mach**2); - ii = 3; jj = 5; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -d(j,k)/(gam*mach**2); - ii = 5; jj = 3; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -(gam-1)*d(j,k)/rho_mean(j); - end do - end do - ar = br - ai = bi+ci - - ! Compute Eigenvalues & Eigenfunctions - call cg(5*(n+1),5*(n+1),ar,ai,wr,wi,zr,zi,fv1,fv2,fv3,ierr) - - ! Find the most unstable wave - call find_unstable_mode(5*(n+1),wr,wi,zr,zi,tr,ti,vr,vi) - - ! Normalize - call normalize_eigvec(5*(n+1),vr,vi,vnr,vni) - - ! Generate instability waves - call generate_wave(5*(n+1),vnr,vni,alpha,beta,wave,shift) + ! Generate instability wave + call s_generate_wave(5*(n+1),wr,wi,zr,zi,alpha,beta,wave,shift) end subroutine s_instability_wave - subroutine generate_wave(nl,vnr,vni,alpha,beta,wave,shift) - integer nl - real(kind(0d0)) :: alpha,beta,ang - real(kind(0d0)), dimension(0:nl-1) :: vnr,vni - real(kind(0d0)), dimension(5,0:m,0:n,0:p) :: wave - real(kind(0d0)) :: shift - integer i,j,k + !> This subroutine generates an instability wave using the most unstable + !! eigenvalue and corresponding eigenvector among the + !! given set of eigenvalues and eigenvectors. + subroutine s_generate_wave(nl,wr,wi,zr,zi,alpha,beta,wave,shift) + integer nl + real(kind(0d0)), dimension(0:nl-1) :: wr,wi !< eigenvalues + real(kind(0d0)), dimension(0:nl-1,0:nl-1) :: zr,zi !< eigenvectors + real(kind(0d0)), dimension(0:nl-1) :: vr,vi,vnr,vni !< most unstable eigenvector + real(kind(0d0)), dimension(5,0:m,0:n,0:p) :: wave + real(kind(0d0)) :: alpha,beta,ang,shift + real(kind(0d0)) :: norm + real(kind(0d0)) :: tr,ti,cr,ci !< temporary memory + integer idx + integer i,j,k + + ! Find the most unstable eigenvalue and corresponding eigenvector + k=0 + do i=1,nl-1 + if (wi(i) .gt. wi(k)) then + k = i + end if + end do + vr = zr(:,k) + vi = zi(:,k) - do i=0,m - do j=0,n - do k=0,p - if (beta .eq. 0) then - ang = alpha*x_cc(i) - else - ang = alpha*x_cc(i)+beta*z_cc(k)+shift - end if - - wave(1,i,j,k) = vnr(j)*cos(ang) & - -vni(j)*sin(ang) ! rho - wave(2,i,j,k) = vnr((n+1)+j)*cos(ang) & - -vni((n+1)+j)*sin(ang) ! u - wave(3,i,j,k) = vnr(2*(n+1)+j)*cos(ang) & - -vni(2*(n+1)+j)*sin(ang) ! v - wave(4,i,j,k) = vnr(3*(n+1)+j)*cos(ang) & - -vni(3*(n+1)+j)*sin(ang) ! w - wave(5,i,j,k) = vnr(4*(n+1)+j)*cos(ang) & - -vni(4*(n+1)+j)*sin(ang) ! T - end do - end do - end do - - end subroutine generate_wave - - subroutine normalize_eigvec(nl,vr,vi,vnr,vni) - integer nl - real(kind(0d0)), dimension(nl) :: vr,vi,vnr,vni - real(kind(0d0)) :: norm - real(kind(0d0)) :: tr,ti,cr,ci - integer i,idx - - norm = 0d0 - do i=1,nl - if (dsqrt(vr(i)**2+vi(i)**2) .gt. dabs(norm)) then - idx = i - if (vr(i) .gt. 0) then - norm = sqrt(vr(i)**2+vi(i)**2) - else - norm =-sqrt(vr(i)**2+vi(i)**2) - end if - end if - end do - - vnr = vr/norm - vni = vi/norm - - tr = vnr(idx) - ti = vni(idx) - do i=1,nl - call cdiv(vnr(i),vni(i),tr,ti,cr,ci) - vnr(i) = cr - vni(i) = ci - end do - - end subroutine normalize_eigvec - - subroutine find_unstable_mode(nl,wr,wi,zr,zi,tr,ti,vr,vi) - integer nl - real(kind(0d0)), dimension(nl) :: wr,wi - real(kind(0d0)), dimension(nl,nl) :: zr,zi - real(kind(0d0)) :: tr,ti - real(kind(0d0)), dimension(nl) :: vr,vi - integer i,k + ! Normalize the eigenvector by its component with the largest modulus. + norm = 0d0 + do i=0,nl-1 + if (dsqrt(vr(i)**2+vi(i)**2) .gt. norm) then + idx = i + norm = dsqrt(vr(i)**2+vi(i)**2) + end if + end do + + tr = vr(idx) + ti = vi(idx) + do i=0,nl-1 + call cdiv(vr(i),vi(i),tr,ti,cr,ci) + vnr(i) = cr + vni(i) = ci + end do - k=1 - do i=2,nl - if (wi(i) .gt. wi(k)) then - k = i - end if - end do - - tr = wr(k) - ti = wi(k) - vr = zr(:,k) - vi = zi(:,k) - - end subroutine find_unstable_mode - -!==================================================================== -! Subroutines & functions for computing eigenvalues and eigenvectors -! which are modified from EISPACK (https://netlib.org/eispack/) -! 1) cg -! 2) cbal -! 3) corth -! 4) comqr2 -! 5) csroot -! 6) cdiv -! 7) pythag -!==================================================================== - subroutine cg(nm,nl,ar,ai,wr,wi,zr,zi,fv1,fv2,fv3,ierr) -! this subroutine calls the recommended sequence of -! subroutines from the eigensystem subroutine package (eispack) -! to find the eigenvalues and eigenvectors (if desired) -! of a complex general matrix. -! -! on input -! -! nm must be set to the row dimension of the two-dimensional -! array parameters as declared in the calling program -! dimension statement. -! -! nl is the order of the matrix a=(ar,ai). -! -! ar and ai contain the real and imaginary parts, -! respectively, of the complex general matrix. -! -! on output -! -! wr and wi contain the real and imaginary parts, -! respectively, of the eigenvalues. -! -! zr and zi contain the real and imaginary parts, -! respectively, of the eigenvectors. -! -! ierr is an integer output variable set equal to an error -! completion code described in the documentation for comqr -! and comqr2. the normal completion code is zero. -! -! fv1, fv2, and fv3 are temporary storage arrays. -! -! questions and comments should be directed to burton s. garbow, -! mathematics and computer science div, argonne national laboratory -! -! this version dated august 1983. -! -! ------------------------------------------------------------------ - integer nm,nl,is1,is2,ierr - real(kind(0d0)), dimension(nm,nl) :: ar,ai,zr,zi - real(kind(0d0)), dimension(nl) :: wr,wi,fv1,fv2,fv3 - - if (nl .le. nm) go to 10 - ierr = 10*nl - go to 50 - -10 call cbal(nm,nl,ar,ai,is1,is2,fv1) - call corth(nm,nl,is1,is2,ar,ai,fv2,fv3) - call comqr2(nm,nl,is1,is2,fv2,fv3,ar,ai,wr,wi,zr,zi,ierr) - if (ierr .ne. 0) go to 50 - call cbabk2(nm,nl,is1,is2,fv1,nl,zr,zi) - -50 return - end subroutine cg - - subroutine cbal(nm,nl,ar,ai,low,igh,scale) -! this subroutine is a translation of the algol procedure -! cbalance, which is a complex version of balance, -! num. math. 13, 293-304(1969) by parlett and reinsch. -! handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). -! -! this subroutine balances a complex matrix and isolates -! eigenvalues whenever possible. -! -! on input -! -! nm must be set to the row dimension of two-dimensional -! array parameters as declared in the calling program -! dimension statement. -! -! nl is the order of the matrix. -! -! ar and ai contain the real and imaginary parts, -! respectively, of the complex matrix to be balanced. -! -! on output -! -! ar and ai contain the real and imaginary parts, -! respectively, of the balanced matrix. -! -! low and igh are two integers such that ar(i,j) and ai(i,j) -! are equal to zero if -! (1) i is greater than j and -! (2) j=1,...,low-1 or i=igh+1,...,nl. -! -! scale contains information determining the -! permutations and scaling factors used. -! -! suppose that the principal submatrix in rows low through igh -! has been balanced, that p(j) denotes the index interchanged -! with j during the permutation step, and that the elements -! of the diagonal matrix used are denoted by d(i,j). then -! scale(j) = p(j), for j = 1,...,low-1 -! = d(j,j) j = low,...,igh -! = p(j) j = igh+1,...,nl. -! the order in which the interchanges are made is nl to igh+1, -! then 1 to low-1. -! -! note that 1 is returned for igh if igh is zero formally. -! -! the algol procedure exc contained in cbalance appears in -! cbal in line. (note that the algol roles of identifiers -! k,l have been reversed.) -! -! arithmetic is real throughout. -! -! questions and comments should be directed to burton s. garbow, -! mathematics and computer science div, argonne national laboratory -! -! this version dated august 1983. -! -! ------------------------------------------------------------------ - integer i,j,k,l,ml,nl,jj,nm,igh,low,iexc - real(kind(0d0)), dimension(nm,nl) :: ar,ai - real(kind(0d0)), dimension(nl) :: scale - real(kind(0d0)) :: c,f,g,r,s,b2,radix - logical noconv - - radix = 16.0d0 - - b2 = radix * radix - k = 1 - l = nl - go to 100 -! .......... in-line procedure for row and -! column exchange .......... -20 scale(ml) = j - if (j .eq. ml) go to 50 - - do 30 i = 1, l - f = ar(i,j) - ar(i,j) = ar(i,ml) - ar(i,ml) = f - f = ai(i,j) - ai(i,j) = ai(i,ml) - ai(i,ml) = f -30 continue - - do 40 i = k, nl - f = ar(j,i) - ar(j,i) = ar(ml,i) - ar(ml,i) = f - f = ai(j,i) - ai(j,i) = ai(ml,i) - ai(ml,i) = f -40 continue - -50 go to (80,130), iexc -! .......... search for rows isolating an eigenvalue -! and push them down .......... -80 if (l .eq. 1) go to 280 - l = l - 1 -! .......... for j=l step -1 until 1 do -- .......... -100 do 120 jj = 1, l - j = l + 1 - jj - - do 110 i = 1, l - if (i .eq. j) go to 110 - if (ar(j,i) .ne. 0.0d0 .or. ai(j,i) .ne. 0.0d0) go to 120 -110 continue - - ml = l - iexc = 1 - go to 20 -120 continue - - go to 140 -! .......... search for columns isolating an eigenvalue -! and push them left .......... -130 k = k + 1 - -140 do 170 j = k, l - - do 150 i = k, l - if (i .eq. j) go to 150 - if (ar(i,j) .ne. 0.0d0 .or. ai(i,j) .ne. 0.0d0) go to 170 -150 continue - - ml = k - iexc = 2 - go to 20 -170 continue -! .......... now balance the submatrix in rows k to l .......... - do 180 i = k, l - scale(i) = 1.0d0 -180 continue -! .......... iterative loop for norm reduction .......... -190 noconv = .false. - - do 270 i = k, l - c = 0.0d0 - r = 0.0d0 - - do 200 j = k, l - if (j .eq. i) go to 200 - c = c + dabs(ar(j,i)) + dabs(ai(j,i)) - r = r + dabs(ar(i,j)) + dabs(ai(i,j)) -200 continue -! .......... guard against zero c or r due to underflow .......... - if (c .eq. 0.0d0 .or. r .eq. 0.0d0) go to 270 - g = r / radix - f = 1.0d0 - s = c + r -210 if (c .ge. g) go to 220 - f = f * radix - c = c * b2 - go to 210 -220 g = r * radix -230 if (c .lt. g) go to 240 - f = f / radix - c = c / b2 - go to 230 -! .......... now balance .......... -240 if ((c + r) / f .ge. 0.95d0 * s) go to 270 - g = 1.0d0 / f - scale(i) = scale(i) * f - noconv = .true. - - do 250 j = k, nl - ar(i,j) = ar(i,j) * g - ai(i,j) = ai(i,j) * g -250 continue - - do 260 j = 1, l - ar(j,i) = ar(j,i) * f - ai(j,i) = ai(j,i) * f -260 continue - -270 continue - - if (noconv) go to 190 - -280 low = k - igh = l - return - end subroutine cbal + ! Generate an instability wave + do i=0,m + do j=0,n + do k=0,p + if (beta .eq. 0) then + ang = alpha*x_cc(i) + else + ang = alpha*x_cc(i)+beta*z_cc(k)+shift + end if + wave(1,i,j,k) = vnr(j)*cos(ang)-vni(j)*sin(ang) ! rho + wave(2,i,j,k) = vnr((n+1)+j)*cos(ang)-vni((n+1)+j)*sin(ang) ! u + wave(3,i,j,k) = vnr(2*(n+1)+j)*cos(ang)-vni(2*(n+1)+j)*sin(ang) ! v + wave(4,i,j,k) = vnr(3*(n+1)+j)*cos(ang)-vni(3*(n+1)+j)*sin(ang) ! w + wave(5,i,j,k) = vnr(4*(n+1)+j)*cos(ang)-vni(4*(n+1)+j)*sin(ang) ! T + end do + end do + end do - subroutine corth(nm,nl,low,igh,ar,ai,ortr,orti) -! this subroutine is a translation of a complex analogue of -! the algol procedure orthes, num. math. 12, 349-368(1968) -! by martin and wilkinson. -! handbook for auto. comp., vol.ii-linear algebra, 339-358(1971). -! -! given a complex general matrix, this subroutine -! reduces a submatrix situated in rows and columns -! low through igh to upper hessenberg form by -! unitary similarity transformations. -! -! on input -! -! nm must be set to the row dimension of two-dimensional -! array parameters as declared in the calling program -! dimension statement. -! -! nl is the order of the matrix. -! -! low and igh are integers determined by the balancing -! subroutine cbal. if cbal has not been used, -! set low=1, igh=nl. -! -! ar and ai contain the real and imaginary parts, -! respectively, of the complex input matrix. -! -! on output -! -! ar and ai contain the real and imaginary parts, -! respectively, of the hessenberg matrix. information -! about the unitary transformations used in the reduction -! is stored in the remaining triangles under the -! hessenberg matrix. -! -! ortr and orti contain further information about the -! transformations. only elements low through igh are used. -! -! calls pythag for dsqrt(a*a + b*b) . -! -! questions and comments should be directed to burton s. garbow, -! mathematics and computer science div, argonne national laboratory -! -! this version dated august 1983. -! -! ------------------------------------------------------------------ - integer i,j,ml,nl,ii,jj,la,mp,nm,igh,kp1,low - real(kind(0d0)),dimension(nm,nl) :: ar,ai - real(kind(0d0)),dimension(igh) :: ortr,orti - real(kind(0d0)) :: f,g,h,fi,fr,scale,c - - integer mll - mll = 6 - - la = igh - 1 - kp1 = low + 1 - if (la .lt. kp1) go to 200 - - do 180 ml = kp1, la - h = 0.0d0 - ortr(ml) = 0.0d0 - orti(ml) = 0.0d0 - scale = 0.0d0 -! .......... scale column (algol tol then not needed) .......... - do 90 i = ml, igh - scale = scale + dabs(ar(i,ml-1)) + dabs(ai(i,ml-1)) -90 continue - if (scale .eq. 0d0) go to 180 - mp = ml + igh -! .......... for i=igh step -1 until ml do -- .......... - do 100 ii = ml, igh - i = mp - ii - ortr(i) = ar(i,ml-1) / scale - orti(i) = ai(i,ml-1) / scale - h = h + ortr(i) * ortr(i) + orti(i) * orti(i) -100 continue -! - g = dsqrt(h) - call pythag(ortr(ml),orti(ml),f) - if (f .eq. 0d0) go to 103 - h = h + f * g - g = g / f - ortr(ml) = (1.0d0 + g) * ortr(ml) - orti(ml) = (1.0d0 + g) * orti(ml) - go to 105 - -103 ortr(ml) = g - ar(ml,ml-1) = scale -! .......... form (i-(u*ut)/h) * a .......... -105 do 130 j = ml, nl - fr = 0.0d0 - fi = 0.0d0 -! .......... for i=igh step -1 until ml do -- .......... - do 110 ii = ml, igh - i = mp - ii - fr = fr + ortr(i) * ar(i,j) + orti(i) * ai(i,j) - fi = fi + ortr(i) * ai(i,j) - orti(i) * ar(i,j) -110 continue -! - fr = fr / h - fi = fi / h -! - do 120 i = ml, igh - ar(i,j) = ar(i,j) - fr * ortr(i) + fi * orti(i) - ai(i,j) = ai(i,j) - fr * orti(i) - fi * ortr(i) -120 continue -! -130 continue -! .......... form (i-(u*ut)/h)*a*(i-(u*ut)/h) .......... - do 160 i = 1, igh - fr = 0.0d0 - fi = 0.0d0 -! .......... for j=igh step -1 until ml do -- .......... - do 140 jj = ml, igh - j = mp - jj - fr = fr + ortr(j) * ar(i,j) - orti(j) * ai(i,j) - fi = fi + ortr(j) * ai(i,j) + orti(j) * ar(i,j) -140 continue -! - fr = fr / h - fi = fi / h -! - do 150 j = ml, igh - ar(i,j) = ar(i,j) - fr * ortr(j) - fi * orti(j) - ai(i,j) = ai(i,j) + fr * orti(j) - fi * ortr(j) -150 continue -! -160 continue -! - ortr(ml) = scale * ortr(ml) - orti(ml) = scale * orti(ml) - ar(ml,ml-1) = -g * ar(ml,ml-1) - ai(ml,ml-1) = -g * ai(ml,ml-1) -180 continue -! -200 return - end subroutine corth + end subroutine s_generate_wave - subroutine comqr2(nm,nl,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr) -! MESHED overflow control WITH vectors of isolated roots (10/19/89 BSG) -! MESHED overflow control WITH triangular multiply (10/30/89 BSG) -! -! this subroutine is a translation of a unitary analogue of the -! algol procedure comlr2, num. math. 16, 181-204(1970) by peters -! and wilkinson. -! handbook for auto. comp., vol.ii-linear algebra, 372-395(1971). -! the unitary analogue substitutes the qr algorithm of francis -! (comp. jour. 4, 332-345(1962)) for the lr algorithm. -! -! this subroutine finds the eigenvalues and eigenvectors -! of a complex upper hessenberg matrix by the qr -! method. the eigenvectors of a complex general matrix -! can also be found if corth has been used to reduce -! this general matrix to hessenberg form. -! -! on input -! -! nm must be set to the row dimension of two-dimensional -! array parameters as declared in the calling program -! dimension statement. -! -! nl is the order of the matrix. -! -! low and igh are integers determined by the balancing -! subroutine cbal. if cbal has not been used, -! set low=1, igh=nl. -! -! ortr and orti contain information about the unitary trans- -! formations used in the reduction by corth, if performed. -! only elements low through igh are used. if the eigenvectors -! of the hessenberg matrix are desired, set ortr(j) and -! orti(j) to 0.0d0 for these elements. -! -! hr and hi contain the real and imaginary parts, -! respectively, of the complex upper hessenberg matrix. -! their lower triangles below the subdiagonal contain further -! information about the transformations which were used in the -! reduction by corth, if performed. if the eigenvectors of -! the hessenberg matrix are desired, these elements may be -! arbitrary. -! -! on output -! -! ortr, orti, and the upper hessenberg portions of hr and hi -! have been destroyed. -! -! wr and wi contain the real and imaginary parts, -! respectively, of the eigenvalues. if an error -! exit is made, the eigenvalues should be correct -! for indices ierr+1,...,nl. -! -! zr and zi contain the real and imaginary parts, -! respectively, of the eigenvectors. the eigenvectors -! are unnormalized. if an error exit is made, none of -! the eigenvectors has been found. -! -! ierr is set to -! zero for normal return, -! j if the limit of 30*nl iterations is exhausted -! while the j-th eigenvalue is being sought. -! -! calls cdiv for complex division. -! calls csroot for complex square root. -! calls pythag for dsqrt(a*a + b*b) . -! -! questions and comments should be directed to burton s. garbow, -! mathematics and computer science div, argonne national laboratory -! -! this version dated october 1989. -! -! ------------------------------------------------------------------ - integer i,j,k,l,ml,nl,en,ii,jj,ll,nm,nn,igh,ip1,& - itn,its,low,lp1,enm1,iend,ierr - real(kind(0d0)),dimension(nm,nl) :: hr,hi,zr,zi - real(kind(0d0)),dimension(nl) :: wr,wi - real(kind(0d0)),dimension(igh) :: ortr,orti - real(kind(0d0)) :: si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr,& - norm,tst1,tst2,c,d -! - ierr = 0 -! .......... initialize eigenvector matrix .......... - do 101 j = 1, nl -! - do 100 i = 1, nl - zr(i,j) = 0.0d0 - zi(i,j) = 0.0d0 -100 continue - zr(j,j) = 1.0d0 -101 continue -! .......... form the matrix of accumulated transformations -! from the information left by corth .......... - iend = igh - low - 1 - if (iend .lt. 0) go to 180 - if (iend .eq. 0) go to 150 - if (iend .gt. 0) go to 105 -! .......... for i=igh-1 step -1 until low+1 do -- .......... -105 do 140 ii = 1, iend - i = igh - ii - if (dabs(ortr(i)) .eq. 0d0 .and. dabs(orti(i)) .eq. 0d0) go to 140 - if (dabs(hr(i,i-1)) .eq. 0d0 .and. dabs(hi(i,i-1)) .eq. 0d0) go to 140 -! .......... norm below is negative of h formed in corth .......... - norm = hr(i,i-1) * ortr(i) + hi(i,i-1) * orti(i) - ip1 = i + 1 - - do 110 k = ip1, igh - ortr(k) = hr(k,i-1) - orti(k) = hi(k,i-1) -110 continue -! - do 130 j = i, igh - sr = 0.0d0 - si = 0.0d0 -! - do 115 k = i, igh - sr = sr + ortr(k) * zr(k,j) + orti(k) * zi(k,j) - si = si + ortr(k) * zi(k,j) - orti(k) * zr(k,j) -115 continue -! - sr = sr / norm - si = si / norm -! - do 120 k = i, igh - zr(k,j) = zr(k,j) + sr * ortr(k) - si * orti(k) - zi(k,j) = zi(k,j) + sr * orti(k) + si * ortr(k) -120 continue -! -130 continue -! -140 continue -! .......... create real subdiagonal elements .......... -150 l = low + 1 -! - do 170 i = l, igh - ll = min0(i+1,igh) - if (dabs(hi(i,i-1)) .eq. 0d0) go to 170 - call pythag(hr(i,i-1),hi(i,i-1),norm) - yr = hr(i,i-1) / norm - yi = hi(i,i-1) / norm - hr(i,i-1) = norm - hi(i,i-1) = 0.0d0 -! - do 155 j = i, nl - si = yr * hi(i,j) - yi * hr(i,j) - hr(i,j) = yr * hr(i,j) + yi * hi(i,j) - hi(i,j) = si -155 continue -! - do 160 j = 1, ll - si = yr * hi(j,i) + yi * hr(j,i) - hr(j,i) = yr * hr(j,i) - yi * hi(j,i) - hi(j,i) = si -160 continue -! - do 165 j = low, igh - si = yr * zi(j,i) + yi * zr(j,i) - zr(j,i) = yr * zr(j,i) - yi * zi(j,i) - zi(j,i) = si -165 continue - -170 continue -! .......... store roots isolated by cbal .......... -180 do 200 i = 1, nl - if (i .ge. low .and. i .le. igh) go to 200 - wr(i) = hr(i,i) - wi(i) = hi(i,i) -200 continue -! - en = igh - tr = 0.0d0 - ti = 0.0d0 - itn = 30*nl -! .......... search for next eigenvalue .......... -220 if (en .lt. low) go to 680 - its = 0 - enm1 = en - 1 -! .......... look for single small sub-diagonal element -! for l=en step -1 until low do -- .......... -240 do 260 ll = low, en - l = en + low - ll - if (l .eq. low) go to 300 - tst1 = dabs(hr(l-1,l-1)) + dabs(hi(l-1,l-1)) & - + dabs(hr(l,l)) + dabs(hi(l,l)) - tst2 = tst1 + dabs(hr(l,l-1)) - if (tst2 .eq. tst1) go to 300 -260 continue -! .......... form shift .......... -300 if (l .eq. en) go to 660 - if (itn .eq. 0) go to 1000 - if (its .eq. 10 .or. its .eq. 20) go to 320 - sr = hr(en,en) - si = hi(en,en) - xr = hr(enm1,en) * hr(en,enm1) - xi = hi(enm1,en) * hr(en,enm1) - if (xr .eq. 0.0d0 .and. xi .eq. 0.0d0) go to 340 - yr = (hr(enm1,enm1) - sr) / 2.0d0 - yi = (hi(enm1,enm1) - si) / 2.0d0 - call csroot(yr**2-yi**2+xr,2.0d0*yr*yi+xi,zzr,zzi) - if (yr * zzr + yi * zzi .ge. 0.0d0) go to 310 - zzr = -zzr - zzi = -zzi -310 call cdiv(xr,xi,yr+zzr,yi+zzi,xr,xi) - sr = sr - xr - si = si - xi - go to 340 -! .......... form exceptional shift .......... -320 sr = dabs(hr(en,enm1)) + dabs(hr(enm1,en-2)) - si = 0.0d0 -! -340 do 360 i = low, en - hr(i,i) = hr(i,i) - sr - hi(i,i) = hi(i,i) - si -360 continue -! - tr = tr + sr - ti = ti + si - its = its + 1 - itn = itn - 1 -! .......... reduce to triangle (rows) .......... - lp1 = l + 1 -! - do 500 i = lp1, en - sr = hr(i,i-1) - hr(i,i-1) = 0.0d0 - call pythag(hr(i-1,i-1),hi(i-1,i-1),c) - call pythag(c,sr,norm) - xr = hr(i-1,i-1) / norm - wr(i-1) = xr - xi = hi(i-1,i-1) / norm - wi(i-1) = xi - hr(i-1,i-1) = norm - hi(i-1,i-1) = 0.0d0 - hi(i,i-1) = sr / norm -! - do 490 j = i, nl - yr = hr(i-1,j) - yi = hi(i-1,j) - zzr = hr(i,j) - zzi = hi(i,j) - hr(i-1,j) = xr * yr + xi * yi + hi(i,i-1) * zzr - hi(i-1,j) = xr * yi - xi * yr + hi(i,i-1) * zzi - hr(i,j) = xr * zzr - xi * zzi - hi(i,i-1) * yr - hi(i,j) = xr * zzi + xi * zzr - hi(i,i-1) * yi -490 continue -! -500 continue -! - si = hi(en,en) - if (dabs(si) .eq. 0d0) go to 540 - call pythag(hr(en,en),si,norm) - sr = hr(en,en) / norm - si = si / norm - hr(en,en) = norm - hi(en,en) = 0.0d0 - if (en .eq. nl) go to 540 - ip1 = en + 1 -! - do 520 j = ip1, nl - yr = hr(en,j) - yi = hi(en,j) - hr(en,j) = sr * yr + si * yi - hi(en,j) = sr * yi - si * yr -520 continue -! .......... inverse operation (columns) .......... -540 do 600 j = lp1, en - xr = wr(j-1) - xi = wi(j-1) -! - do 580 i = 1, j - yr = hr(i,j-1) - yi = 0.0d0 - zzr = hr(i,j) - zzi = hi(i,j) - if (i .eq. j) go to 560 - yi = hi(i,j-1) - hi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi -560 hr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr - hr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr - hi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi -580 continue -! - do 590 i = low, igh - yr = zr(i,j-1) - yi = zi(i,j-1) - zzr = zr(i,j) - zzi = zi(i,j) - zr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr - zi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi - zr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr - zi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi -590 continue - -600 continue -! - if (dabs(si) .eq. 0d0) go to 240 -! - do 630 i = 1, en - yr = hr(i,en) - yi = hi(i,en) - hr(i,en) = sr * yr - si * yi - hi(i,en) = sr * yi + si * yr -630 continue -! - do 640 i = low, igh - yr = zr(i,en) - yi = zi(i,en) - zr(i,en) = sr * yr - si * yi - zi(i,en) = sr * yi + si * yr -640 continue -! - go to 240 -! .......... a root found .......... -660 hr(en,en) = hr(en,en) + tr - wr(en) = hr(en,en) - hi(en,en) = hi(en,en) + ti - wi(en) = hi(en,en) - en = enm1 - go to 220 -! .......... all roots found. backsubstitute to find -! vectors of upper triangular form .......... -680 norm = 0.0d0 -! - do i = 1, nl - do j = i, nl - tr = dabs(hr(i,j)) + dabs(hi(i,j)) - if (tr .gt. norm) norm = tr - end do - end do -! - if (nl .eq. 1 .or. norm .eq. 0d0) go to 1001 -! .......... for en=nl step -1 until 2 do -- .......... - do 800 nn = 2, nl - en = nl + 2 - nn - xr = wr(en) - xi = wi(en) - hr(en,en) = 1.0d0 - hi(en,en) = 0.0d0 - enm1 = en - 1 -! .......... for i=en-1 step -1 until 1 do -- .......... - do 780 ii = 1, enm1 - i = en - ii - zzr = 0.0d0 - zzi = 0.0d0 - ip1 = i + 1 - - do 740 j = ip1, en - zzr = zzr + hr(i,j) * hr(j,en) - hi(i,j) * hi(j,en) - zzi = zzi + hr(i,j) * hi(j,en) + hi(i,j) * hr(j,en) -740 continue -! - yr = xr - wr(i) - yi = xi - wi(i) - if (yr .ne. 0.0d0 .or. yi .ne. 0.0d0) go to 765 - tst1 = norm - yr = tst1 -760 yr = 0.01d0 * yr - tst2 = norm + yr - if (tst2 .gt. tst1) go to 760 -765 continue - call cdiv(zzr,zzi,yr,yi,hr(i,en),hi(i,en)) -! .......... overflow control .......... - tr = dabs(hr(i,en)) + dabs(hi(i,en)) - if (tr .eq. 0.0d0) go to 780 - tst1 = tr - tst2 = tst1 + 1.0d0/tst1 - if (tst2 .gt. tst1) go to 780 - do 770 j = i, en - hr(j,en) = hr(j,en)/tr - hi(j,en) = hi(j,en)/tr -770 continue -! -780 continue -! -800 continue -! .......... end backsubstitution .......... -! .......... vectors of isolated roots .......... - do 840 i = 1, nl - if (i .ge. low .and. i .le. igh) go to 840 -! - do 820 j = I, nl - zr(i,j) = hr(i,j) - zi(i,j) = hi(i,j) -820 continue -! -840 continue -! .......... multiply by transformation matrix to give -! vectors of original full matrix. -! for j=nl step -1 until low do -- .......... - do jj = low, nl - j = nl + low - jj - ml = min0(j,igh) -! - do i = low, igh - zzr = 0.0d0 - zzi = 0.0d0 -! - do 860 k = low, ml - zzr = zzr + zr(i,k) * hr(k,j) - zi(i,k) * hi(k,j) - zzi = zzi + zr(i,k) * hi(k,j) + zi(i,k) * hr(k,j) -860 continue -! - zr(i,j) = zzr - zi(i,j) = zzi - end do - end do -! - go to 1001 -! .......... set error -- all eigenvalues have not -! converged after 30*nl iterations .......... -1000 ierr = en -1001 return - end subroutine comqr2 - - subroutine cbabk2(nm,nl,low,igh,scale,ml,zr,zi) -! - integer i,j,k,ml,nl,ii,nm,igh,low - double precision scale(nl),zr(nm,ml),zi(nm,ml) - double precision s -! -! this subroutine is a translation of the algol procedure -! cbabk2, which is a complex version of balbak, -! num. math. 13, 293-304(1969) by parlett and reinsch. -! handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). -! -! this subroutine forms the eigenvectors of a complex general -! matrix by back transforming those of the corresponding -! balanced matrix determined by cbal. -! -! on input -! -! nm must be set to the row dimension of two-dimensional -! array parameters as declared in the calling program -! dimension statement. -! -! nl is the order of the matrix. -! -! low and igh are integers determined by cbal. -! -! scale contains information determining the permutations -! and scaling factors used by cbal. -! -! ml is the number of eigenvectors to be back transformed. -! -! zr and zi contain the real and imaginary parts, -! respectively, of the eigenvectors to be -! back transformed in their first ml columns. -! -! on output -! -! zr and zi contain the real and imaginary parts, -! respectively, of the transformed eigenvectors -! in their first ml columns. -! -! questions and comments should be directed to burton s. garbow, -! mathematics and computer science div, argonne national laboratory -! -! this version dated august 1983. -! -! ------------------------------------------------------------------ -! - if (ml .eq. 0) go to 200 - if (igh .eq. low) go to 120 -! - do 110 i = low, igh - s = scale(i) -! .......... left hand eigenvectors are back transformed -! if the foregoing statement is replaced by -! s=1.0d0/scale(i). .......... - do 100 j = 1, ml - zr(i,j) = zr(i,j) * s - zi(i,j) = zi(i,j) * s -100 continue -! -110 continue -! .......... for i=low-1 step -1 until 1, -! igh+1 step 1 until nl do -- .......... -120 do 140 ii = 1, nl - i = ii - if (i .ge. low .and. i .le. igh) go to 140 - if (i .lt. low) i = low - ii - k = scale(i) - if (k .eq. i) go to 140 -! - do 130 j = 1, ml - s = zr(i,j) - zr(i,j) = zr(k,j) - zr(k,j) = s - s = zi(i,j) - zi(i,j) = zi(k,j) - zi(k,j) = s -130 continue -! -140 continue -! -200 return - end subroutine cbabk2 - - subroutine csroot(xr,xi,yr,yi) - real(kind(0d0)) :: xr,xi,yr,yi -! -! (yr,yi) = complex dsqrt(xr,xi) -! branch chosen so that yr .ge. 0.0 and sign(yi) .eq. sign(xi) -! - real(kind(0d0)) :: s,tr,ti,c - tr = xr - ti = xi - call pythag(tr,ti,c) - s = dsqrt(0.5d0*(c + dabs(tr))) - if (tr .ge. 0.0d0) yr = s - if (ti .lt. 0.0d0) s = -s - if (tr .le. 0.0d0) yi = s - if (tr .lt. 0.0d0) yr = 0.5d0*(ti/yi) - if (tr .gt. 0.0d0) yi = 0.5d0*(ti/yr) - return - end subroutine csroot - - subroutine cdiv(ar,ai,br,bi,cr,ci) - real(kind(0d0)) :: ar,ai,br,bi,cr,ci -! -! complex division, (cr,ci) = (ar,ai)/(br,bi) -! - real(kind(0d0)) :: s,ars,ais,brs,bis - s = dabs(br) + dabs(bi) - ars = ar/s - ais = ai/s - brs = br/s - bis = bi/s - s = brs**2 + bis**2 - cr = (ars*brs + ais*bis)/s - ci = (ais*brs - ars*bis)/s - return - end subroutine cdiv - - subroutine pythag(a,b,c) - real(kind(0d0)) :: a,b,c -! -! finds dsqrt(a**2+b**2) without overflow or destructive underflow -! - real(kind(0d0)) :: p,r,s,t,u - p = dmax1(dabs(a),dabs(b)) - if (p .eq. 0.0d0) go to 20 - r = (dmin1(dabs(a),dabs(b))/p)**2 -10 continue - t = 4.0d0 + r - if (t .eq. 4.0d0) go to 20 - s = r/t - u = 1.0d0 + 2.0d0*s - p = u*p - r = (s/u)**2 * r - go to 10 -20 c = p - return - end subroutine pythag - - - - !> Deallocation procedures for the module subroutine s_finalize_initial_condition_module() ! --------------------- From 9959b37d7459a8e40128de3e323f4c3fb52ceec9 Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Tue, 28 Mar 2023 22:16:19 -0700 Subject: [PATCH 07/17] make indent correct --- src/pre_process/m_checker.f90 | 4 +- src/pre_process/m_global_parameters.fpp | 2 - src/pre_process/m_initial_condition.fpp | 54 ++++++++++++------------- 3 files changed, 29 insertions(+), 31 deletions(-) diff --git a/src/pre_process/m_checker.f90 b/src/pre_process/m_checker.f90 index 9102781686..a3a0cfa1a7 100644 --- a/src/pre_process/m_checker.f90 +++ b/src/pre_process/m_checker.f90 @@ -569,13 +569,13 @@ subroutine s_check_inputs() if ((vel_profile .eqv. .true.) .and. (n == 0)) then call s_mpi_abort('Unsupported choices of the combination of values for '//& 'vel_profile and n. Exiting ...') - end if + end if ! Constraints on the instability wave if ((instability_wave .eqv. .true.) .and. (n == 0)) then call s_mpi_abort('Unsupported choices of the combination of values for '//& 'instability_wave and n. Exiting ...') - end if + end if ! Constraints on the stiffened equation of state fluids parameters do i = 1, num_fluids diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 31c33deb15..e8df85723a 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -98,7 +98,6 @@ module m_global_parameters integer :: precision !< Precision of output files logical :: vel_profile !< Set hypertangent streamwise velocity profile - logical :: instability_wave !< Superimpose instability waves to surrounding fluid flow ! Perturb density of surrounding air so as to break symmetry of grid @@ -291,7 +290,6 @@ contains patch_icpp(i)%p0 = dflt_real patch_icpp(i)%m0 = dflt_real - end do ! Tait EOS diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp index b4426073ee..32e5e3cbb0 100644 --- a/src/pre_process/m_initial_condition.fpp +++ b/src/pre_process/m_initial_condition.fpp @@ -96,7 +96,7 @@ contains subroutine s_generate_initial_condition() ! ---------------------------- integer :: i !< Generic loop operator - + ! Converting the conservative variables to the primitive ones given ! preexisting initial condition data files were read in on start-up if (old_ic) then @@ -215,7 +215,7 @@ contains if (perturb_flow) call s_perturb_surrounding_flow() if (perturb_sph) call s_perturb_sphere() - if (instability_wave) call s_superposition_instability_wave() + if (instability_wave) call s_superposition_instability_wave() ! Converting the primitive variables to the conservative ones call s_convert_primitive_to_conservative_variables(q_prim_vf, & @@ -292,9 +292,9 @@ contains !> This subroutine computes velocity perturbations for a temporal mixing !! layer with hypertangent mean streamwise velocity profile !! obtained from linear stability analysis. For a 2D case, - !! instability waves with spatial wavenumbers, (4,0), (2,0), - !! and (1,0) are superposed. For a 3D waves, (4,4), (4,-4), - !! (2,2), (2,-2), (1,1), (1,-1) areadded on top of 2D waves. + !! instability waves with spatial wavenumbers, (4,0), (2,0), + !! and (1,0) are superposed. For a 3D waves, (4,4), (4,-4), + !! (2,2), (2,-2), (1,1), (1,-1) areadded on top of 2D waves. subroutine s_superposition_instability_wave() ! ------------------------ real(kind(0d0)), dimension(5,0:m,0:n,0:p) :: wave,wave1,wave2,wave_tmp real(kind(0d0)) :: tr,ti @@ -307,7 +307,7 @@ contains wave1 = 0d0 wave2 = 0d0 - ! Compute 2D waves + ! Compute 2D waves call s_instability_wave(2*pi*4.0/Ly,0d0,tr,ti,wave_tmp,0d0) wave1 = wave1 + wave_tmp call s_instability_wave(2*pi*2.0/Ly,0d0,tr,ti,wave_tmp,0d0) @@ -317,7 +317,7 @@ contains wave = wave1*0.05 if (p > 0) then - ! Compute 3D waves with phase shifts. + ! Compute 3D waves with phase shifts. call s_instability_wave(2*pi*4.0/Ly, 2*pi*4.0/Ly,tr,ti,wave_tmp,2*pi*11d0/31d0) wave2 = wave2 + wave_tmp call s_instability_wave(2*pi*2.0/Ly, 2*pi*2.0/Ly,tr,ti,wave_tmp,2*pi*13d0/31d0) @@ -333,7 +333,7 @@ contains wave = wave + 0.15*wave2 end if - ! Superpose velocity perturbuations (instability waves) to the velocity field + ! Superpose velocity perturbuations (instability waves) to the velocity field do k = 0, p do j = 0, n do i = 0, m @@ -349,10 +349,10 @@ contains end subroutine s_superposition_instability_wave ! ---------------------- !> This subroutine computes instability waves corresponding to the spatial - !! wavenumbers (alpha, beta) for x and z directions, - !! respectively. The eigenvalue problem is derived from the - !! linearized Euler equations with parallel mean flow - !! assumption (See Sandham 1989 PhD thesis for details). + !! wavenumbers (alpha, beta) for x and z directions, + !! respectively. The eigenvalue problem is derived from the + !! linearized Euler equations with parallel mean flow + !! assumption (See Sandham 1989 PhD thesis for details). subroutine s_instability_wave(alpha,beta,tr,ti,wave,shift) real(kind(0d0)),intent(in) :: alpha, beta !< spatial wavenumbers real(kind(0d0)),dimension(0:n) :: rho_mean, u_mean, t_mean !< mean profiles @@ -387,7 +387,7 @@ contains ! Compute differential operator in y-dir ! based on 4th order central difference (inner) - ! and 2nd order central difference (near boundaries) + ! and 2nd order central difference (near boundaries) dy = y_cc(1)-y_cc(0) d=0d0 d(1,0)=-1/(2*dy) @@ -414,8 +414,8 @@ contains end do ! Compute B and C, then A = B + C - ! B includes terms without differential operator, and - ! C includes terms with differential operator + ! B includes terms without differential operator, and + ! C includes terms with differential operator br=0d0 bi=0d0 ci=0d0 @@ -443,7 +443,7 @@ contains ii = 5; jj = 4; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = (gam-1)*beta/rho_mean(j); ii = 5; jj = 5; br((ii-1)*(n+1)+j,(jj-1)*(n+1)+j) = alpha*u_mean(j); - do k=0,n + do k=0,n ii = 1; jj = 3; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -rho_mean(j)*d(j,k); ii = 3; jj = 1; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -t_mean(j)*d(j,k)/(rho_mean(j)*gam*mach**2); ii = 3; jj = 5; ci((ii-1)*(n+1)+j,(jj-1)*(n+1)+k) = -d(j,k)/(gam*mach**2); @@ -462,21 +462,21 @@ contains end subroutine s_instability_wave !> This subroutine generates an instability wave using the most unstable - !! eigenvalue and corresponding eigenvector among the - !! given set of eigenvalues and eigenvectors. + !! eigenvalue and corresponding eigenvector among the + !! given set of eigenvalues and eigenvectors. subroutine s_generate_wave(nl,wr,wi,zr,zi,alpha,beta,wave,shift) integer nl - real(kind(0d0)), dimension(0:nl-1) :: wr,wi !< eigenvalues - real(kind(0d0)), dimension(0:nl-1,0:nl-1) :: zr,zi !< eigenvectors + real(kind(0d0)), dimension(0:nl-1) :: wr,wi !< eigenvalues + real(kind(0d0)), dimension(0:nl-1,0:nl-1) :: zr,zi !< eigenvectors real(kind(0d0)), dimension(0:nl-1) :: vr,vi,vnr,vni !< most unstable eigenvector real(kind(0d0)), dimension(5,0:m,0:n,0:p) :: wave - real(kind(0d0)) :: alpha,beta,ang,shift - real(kind(0d0)) :: norm + real(kind(0d0)) :: alpha,beta,ang,shift + real(kind(0d0)) :: norm real(kind(0d0)) :: tr,ti,cr,ci !< temporary memory - integer idx + integer idx integer i,j,k - ! Find the most unstable eigenvalue and corresponding eigenvector + ! Find the most unstable eigenvalue and corresponding eigenvector k=0 do i=1,nl-1 if (wi(i) .gt. wi(k)) then @@ -486,7 +486,7 @@ contains vr = zr(:,k) vi = zi(:,k) - ! Normalize the eigenvector by its component with the largest modulus. + ! Normalize the eigenvector by its component with the largest modulus. norm = 0d0 do i=0,nl-1 if (dsqrt(vr(i)**2+vi(i)**2) .gt. norm) then @@ -503,12 +503,12 @@ contains vni(i) = ci end do - ! Generate an instability wave + ! Generate an instability wave do i=0,m do j=0,n do k=0,p if (beta .eq. 0) then - ang = alpha*x_cc(i) + ang = alpha*x_cc(i) else ang = alpha*x_cc(i)+beta*z_cc(k)+shift end if From 3c0d121d889024d2e1d53751a393a86cbdbb3720 Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Tue, 28 Mar 2023 22:18:54 -0700 Subject: [PATCH 08/17] make indent correct --- src/pre_process/m_assign_variables.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pre_process/m_assign_variables.f90 b/src/pre_process/m_assign_variables.f90 index 0b36a13f07..829b539de9 100644 --- a/src/pre_process/m_assign_variables.f90 +++ b/src/pre_process/m_assign_variables.f90 @@ -485,7 +485,7 @@ subroutine s_assign_patch_species_primitive_variables_bubbles(patch_id, j, k, l, + (1d0 - eta)*orig_prim_vf(i + cont_idx%end)) end do - ! Set streamwise velocity to hypertangent function of y + ! Set streamwise velocity to hypertangent function of y if (vel_profile) then q_prim_vf(1 + cont_idx%end)%sf(j, k, l) = & (eta*patch_icpp(patch_id)%vel(1)*tanh(y_cc(k)) & @@ -686,7 +686,7 @@ subroutine s_assign_patch_species_primitive_variables(patch_id, j, k, l, & + (1d0 - eta)*orig_prim_vf(i + cont_idx%end)) end do - ! Set streamwise velocity to hypertangent function of y + ! Set streamwise velocity to hypertangent function of y if (vel_profile) then q_prim_vf(1 + cont_idx%end)%sf(j, k, l) = & (eta*patch_icpp(patch_id)%vel(1)*tanh(y_cc(k)) & From 6f4ee0fb88f96304db08ecf352a19384f16ecb49 Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Tue, 28 Mar 2023 22:46:22 -0700 Subject: [PATCH 09/17] Add comments --- src/pre_process/m_initial_condition.fpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp index 32e5e3cbb0..4a3c346888 100644 --- a/src/pre_process/m_initial_condition.fpp +++ b/src/pre_process/m_initial_condition.fpp @@ -348,11 +348,11 @@ contains end subroutine s_superposition_instability_wave ! ---------------------- - !> This subroutine computes instability waves corresponding to the spatial - !! wavenumbers (alpha, beta) for x and z directions, - !! respectively. The eigenvalue problem is derived from the - !! linearized Euler equations with parallel mean flow - !! assumption (See Sandham 1989 PhD thesis for details). + !> This subroutine computes instability waves for a given set of spatial + !! wavenumbers (alpha, beta) in x and z directions. + !! The eigenvalue problem is derived from the linearized + !! Euler equations with parallel mean flow assumption + !! (See Sandham 1989 PhD thesis for details). subroutine s_instability_wave(alpha,beta,tr,ti,wave,shift) real(kind(0d0)),intent(in) :: alpha, beta !< spatial wavenumbers real(kind(0d0)),dimension(0:n) :: rho_mean, u_mean, t_mean !< mean profiles From c3f4dcd60c3d905e02e1460301044a9ee41e5f53 Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Mon, 3 Apr 2023 15:28:22 -0700 Subject: [PATCH 10/17] Update documentation --- docs/documentation/case.md | 32 ++++++++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index da6ee6d7f5..3736de182c 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -66,8 +66,9 @@ There are multiple sets of parameters that must be specified in the python input 6. [Formatted Database and Structure Parameters](#6-formatted-output) 7. [(Optional) Acoustic Source Parameters](#7-acoustic-source) 8. [(Optional) Ensemble-Averaged Bubble Model Parameters](#8-ensemble-averaged-bubble-model) +9. [(Optional) Velocity field parameters] (#9-velocity-field) -Items 7 and 8 are optional sets of parameters that activate the acoustic source model and ensemble-averaged bubble model, respectively. +Items 7, 8, and 9 are optional sets of parameters that activate the acoustic source model, ensemble-averaged bubble model, and initial velocity field modification, respectively. Definition of the parameters is described in the following subsections. ### 1. Runtime @@ -295,10 +296,10 @@ Note that `time_stepper` $=$ 3 specifies the total variation diminishing (TVD), | `format` | Integer | Output format. [1]: Silo-HDF5; [2] Binary | | `precision` | Integer | [1] Single; [2] Double | | `parallel_io` | Logical | Parallel I/O | -| `cons_vars_wrt` | Logical | Write conservative variables \| +| `cons_vars_wrt` | Logical | Write conservative variables | | `prim_vars_wrt` | Logical | Write primitive variables | | `fourier_decomp` | Logical | Apply a spatial Fourier decomposition to the output variables | -| `alpha_rho_wrt(i)` | Logical | Add the partial density of the fluid $i$ to the database \| +| `alpha_rho_wrt(i)` | Logical | Add the partial density of the fluid $i$ to the database | | `rho_wrt` | Logical | Add the mixture density to the database | | `mom_wrt(i)` | Logical | Add the $i$-direction momentum to the database | | `vel_wrt(i)` | Logical | Add the $i$-direction velocity to the database | @@ -307,11 +308,12 @@ Note that `time_stepper` $=$ 3 specifies the total variation diminishing (TVD), | `alpha_wrt(i)` | Logical | Add the volume fraction of fluid $i$ to the database | | `gamma_wrt` | Logical | Add the specific heat ratio function to the database | | `heat_ratio_wrt` | Logical | Add the specific heat ratio to the database | -| `pi_inf_wrt` | Logical | Add the liquid stiffness function to the database \| +| `pi_inf_wrt` | Logical | Add the liquid stiffness function to the database | | `pres_inf_wrt` | Logical | Add the liquid stiffness to the formatted database | | `c_wrt` | Logical | Add the sound speed to the database | | `omega_wrt(i)` | Logical | Add the $i$-direction vorticity to the database | | `schlieren_wrt` | Logical | Add the numerical schlieren to the database| +| `qm_wrt` | Logical | Add the Q-criterion to the database| | `fd_order` | Integer | Order of finite differences for computing the vorticity and the numerical Schlieren function [1,2,4] | | `schlieren_alpha(i)` | Real | Intensity of the numerical Schlieren computed via `alpha(i)` | | `probe_wrt` | Logical | Write the flow chosen probes data files for each time step | @@ -433,6 +435,28 @@ When `polytropic` is set `False`, the gas compression is modeled as non-polytrop `gamma_v`, `M_v`, `mu_v`, and `k_v` specify the specific heat ratio, molecular weight, viscosity, and thermal conductivity of a chosen component. Implementation of the parameterse into the model follow [Ando (2010)](references.md#Ando10). +### 9. Velocity field +| Parameter | Type | Description | +| ---: | :----: | :--- | +| 'perturb_flow' | Logical | Perturb the initlal velocity field by random noise | +| 'perturb_sph' | Logical | Perturb the initial partial density by random noise | +| 'perturb_sph_fluid' | Integer | Fluid component whose partial density to be perturbed | +| `vel_profile` | Logical | Set the mean streamwise velocity to hyperbolic tangent profile | +| `instability_wave` | Logical | Perturb the initial velocity field by instability waves | + +The table lists velocity field parameters. The parameters are optionally used to define initial velocity profiles and perturbations. + +- 'perturb_flow' activates the perturbation of initial velocity by random noise. + +- 'perturb_sph' activates the perturbation of intial partial density by random noise. + +- 'perturb_sph_fluid' specifies the fluid component whose the partial density to be perturbed. + +- 'vel_profile' activates setting the mean streamwise velocity to hyperbolic tangent profile. This option works only for 2D and 3D cases. + +- 'instability_wave' activates the perturbation of initial velocity by instability waves obtained from linear stability analysis for a mixing layer with hyperbolic tangent mean streamwise velocity profile. This option only works for 2D and 3D cases, together with 'vel_profile'=true. + + ## Enumerations ### Boundary conditions From 671cfaffb9c60f45b56cb6240d67c0bddf4ef742 Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Mon, 3 Apr 2023 15:43:09 -0700 Subject: [PATCH 11/17] Update documentation --- docs/documentation/case.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 3736de182c..2a16e17e8f 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -68,7 +68,7 @@ There are multiple sets of parameters that must be specified in the python input 8. [(Optional) Ensemble-Averaged Bubble Model Parameters](#8-ensemble-averaged-bubble-model) 9. [(Optional) Velocity field parameters] (#9-velocity-field) -Items 7, 8, and 9 are optional sets of parameters that activate the acoustic source model, ensemble-averaged bubble model, and initial velocity field modification, respectively. +Items 7, 8, and 9 are optional sets of parameters that activate the acoustic source model, ensemble-averaged bubble model, and initial velocity field setup, respectively. Definition of the parameters is described in the following subsections. ### 1. Runtime @@ -436,6 +436,7 @@ When `polytropic` is set `False`, the gas compression is modeled as non-polytrop Implementation of the parameterse into the model follow [Ando (2010)](references.md#Ando10). ### 9. Velocity field + | Parameter | Type | Description | | ---: | :----: | :--- | | 'perturb_flow' | Logical | Perturb the initlal velocity field by random noise | From 0638b53d9f170590a6a69613a5d4e633a801659b Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Mon, 3 Apr 2023 15:53:23 -0700 Subject: [PATCH 12/17] Update documentation --- docs/documentation/case.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 2a16e17e8f..8392bd6747 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -66,7 +66,7 @@ There are multiple sets of parameters that must be specified in the python input 6. [Formatted Database and Structure Parameters](#6-formatted-output) 7. [(Optional) Acoustic Source Parameters](#7-acoustic-source) 8. [(Optional) Ensemble-Averaged Bubble Model Parameters](#8-ensemble-averaged-bubble-model) -9. [(Optional) Velocity field parameters] (#9-velocity-field) +9. [(Optional) Velocity Field Parameters] (#9-velocity-field) Items 7, 8, and 9 are optional sets of parameters that activate the acoustic source model, ensemble-averaged bubble model, and initial velocity field setup, respectively. Definition of the parameters is described in the following subsections. @@ -435,7 +435,7 @@ When `polytropic` is set `False`, the gas compression is modeled as non-polytrop `gamma_v`, `M_v`, `mu_v`, and `k_v` specify the specific heat ratio, molecular weight, viscosity, and thermal conductivity of a chosen component. Implementation of the parameterse into the model follow [Ando (2010)](references.md#Ando10). -### 9. Velocity field +### 9. Velocity Field | Parameter | Type | Description | | ---: | :----: | :--- | From 4bd3fb823107a2d42551d46f7ae1e178f030b746 Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Mon, 3 Apr 2023 15:57:32 -0700 Subject: [PATCH 13/17] Update documentation --- docs/documentation/case.md | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 8392bd6747..324e92d6ea 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -66,7 +66,7 @@ There are multiple sets of parameters that must be specified in the python input 6. [Formatted Database and Structure Parameters](#6-formatted-output) 7. [(Optional) Acoustic Source Parameters](#7-acoustic-source) 8. [(Optional) Ensemble-Averaged Bubble Model Parameters](#8-ensemble-averaged-bubble-model) -9. [(Optional) Velocity Field Parameters] (#9-velocity-field) +9. [(Optional) Velocity Field Parameters](#9-velocity-field) Items 7, 8, and 9 are optional sets of parameters that activate the acoustic source model, ensemble-averaged bubble model, and initial velocity field setup, respectively. Definition of the parameters is described in the following subsections. @@ -439,23 +439,23 @@ Implementation of the parameterse into the model follow [Ando (2010)](references | Parameter | Type | Description | | ---: | :----: | :--- | -| 'perturb_flow' | Logical | Perturb the initlal velocity field by random noise | -| 'perturb_sph' | Logical | Perturb the initial partial density by random noise | -| 'perturb_sph_fluid' | Integer | Fluid component whose partial density to be perturbed | +| `perturb_flow` | Logical | Perturb the initlal velocity field by random noise | +| `perturb_sph` | Logical | Perturb the initial partial density by random noise | +| `perturb_sph_fluid` | Integer | Fluid component whose partial density to be perturbed | | `vel_profile` | Logical | Set the mean streamwise velocity to hyperbolic tangent profile | | `instability_wave` | Logical | Perturb the initial velocity field by instability waves | The table lists velocity field parameters. The parameters are optionally used to define initial velocity profiles and perturbations. -- 'perturb_flow' activates the perturbation of initial velocity by random noise. +- `perturb_flow` activates the perturbation of initial velocity by random noise. -- 'perturb_sph' activates the perturbation of intial partial density by random noise. +- `perturb_sph` activates the perturbation of intial partial density by random noise. -- 'perturb_sph_fluid' specifies the fluid component whose the partial density to be perturbed. +- `perturb_sph_fluid` specifies the fluid component whose the partial density to be perturbed. -- 'vel_profile' activates setting the mean streamwise velocity to hyperbolic tangent profile. This option works only for 2D and 3D cases. +- `vel_profile` activates setting the mean streamwise velocity to hyperbolic tangent profile. This option works only for 2D and 3D cases. -- 'instability_wave' activates the perturbation of initial velocity by instability waves obtained from linear stability analysis for a mixing layer with hyperbolic tangent mean streamwise velocity profile. This option only works for 2D and 3D cases, together with 'vel_profile'=true. +- `instability_wave` activates the perturbation of initial velocity by instability waves obtained from linear stability analysis for a mixing layer with hyperbolic tangent mean streamwise velocity profile. This option only works for 2D and 3D cases, together with 'vel_profile'=true. ## Enumerations From 8baac1384c8014ec1eb1b150ab50492ee420ece4 Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Mon, 3 Apr 2023 15:59:26 -0700 Subject: [PATCH 14/17] Update documentation --- docs/documentation/case.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 324e92d6ea..5fc25c255d 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -455,7 +455,7 @@ The table lists velocity field parameters. The parameters are optionally used to - `vel_profile` activates setting the mean streamwise velocity to hyperbolic tangent profile. This option works only for 2D and 3D cases. -- `instability_wave` activates the perturbation of initial velocity by instability waves obtained from linear stability analysis for a mixing layer with hyperbolic tangent mean streamwise velocity profile. This option only works for 2D and 3D cases, together with 'vel_profile'=true. +- `instability_wave` activates the perturbation of initial velocity by instability waves obtained from linear stability analysis for a mixing layer with hyperbolic tangent mean streamwise velocity profile. This option only works for 2D and 3D cases, together with `vel_profile`=true. ## Enumerations From d10fc234c706f7f99bc8ed3e3f479abee64daaed Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Mon, 3 Apr 2023 16:04:05 -0700 Subject: [PATCH 15/17] Update documentation --- docs/documentation/case.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 5fc25c255d..0311a5806e 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -66,7 +66,7 @@ There are multiple sets of parameters that must be specified in the python input 6. [Formatted Database and Structure Parameters](#6-formatted-output) 7. [(Optional) Acoustic Source Parameters](#7-acoustic-source) 8. [(Optional) Ensemble-Averaged Bubble Model Parameters](#8-ensemble-averaged-bubble-model) -9. [(Optional) Velocity Field Parameters](#9-velocity-field) +9. [(Optional) Velocity Field Setup Parameters](#9-velocity-field-setup) Items 7, 8, and 9 are optional sets of parameters that activate the acoustic source model, ensemble-averaged bubble model, and initial velocity field setup, respectively. Definition of the parameters is described in the following subsections. @@ -435,7 +435,7 @@ When `polytropic` is set `False`, the gas compression is modeled as non-polytrop `gamma_v`, `M_v`, `mu_v`, and `k_v` specify the specific heat ratio, molecular weight, viscosity, and thermal conductivity of a chosen component. Implementation of the parameterse into the model follow [Ando (2010)](references.md#Ando10). -### 9. Velocity Field +### 9. Velocity Field Setup | Parameter | Type | Description | | ---: | :----: | :--- | From 44e75ed82a558d25dfda06fa45544b3d6ca1979f Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Mon, 3 Apr 2023 16:08:35 -0700 Subject: [PATCH 16/17] Update documentation --- docs/documentation/case.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 0311a5806e..899e039a54 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -455,7 +455,7 @@ The table lists velocity field parameters. The parameters are optionally used to - `vel_profile` activates setting the mean streamwise velocity to hyperbolic tangent profile. This option works only for 2D and 3D cases. -- `instability_wave` activates the perturbation of initial velocity by instability waves obtained from linear stability analysis for a mixing layer with hyperbolic tangent mean streamwise velocity profile. This option only works for 2D and 3D cases, together with `vel_profile`=true. +- `instability_wave` activates the perturbation of initial velocity by instability waves obtained from linear stability analysis for a mixing layer with hyperbolic tangent mean streamwise velocity profile. This option only works for 2D and 3D cases, together with `vel_profile`=TRUE. ## Enumerations From e76c85da5f1392cb69c6a2e5773d9b12bb3a34ba Mon Sep 17 00:00:00 2001 From: lee-hyeoksu Date: Mon, 3 Apr 2023 16:31:57 -0700 Subject: [PATCH 17/17] Update example case file --- examples/3D_turb_mixing/case.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/3D_turb_mixing/case.py b/examples/3D_turb_mixing/case.py index acfa7daab9..4ca8c18778 100644 --- a/examples/3D_turb_mixing/case.py +++ b/examples/3D_turb_mixing/case.py @@ -22,9 +22,9 @@ Lz = 59.0 # Number of grid cells -Nx = 319 -Ny = 319 -Nz = 319 +Nx = 255 +Ny = 255 +Nz = 255 # Grid spacing dx = Lx/float(Nx)