1 / 20

Two Fortran codes are provided for you: Explicit:

Two Fortran codes are provided for you: Explicit: FTBS (forward in time, backward in space) , name of code is ftbs_quasi1d.f Implicit: LBI (Linearized Block Implicit method), name of folder containing files is quasi1d_lbi. Quasi1d_lbi folder contains the following files:

ide
Download Presentation

Two Fortran codes are provided for you: Explicit:

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Two Fortran codes are provided for you: Explicit: FTBS (forward in time, backward in space), name of code is ftbs_quasi1d.f Implicit: LBI (Linearized Block Implicit method), name of folder containing files is quasi1d_lbi

  2. Quasi1d_lbi folder contains the following files: input.f - sets all input parameters, geometry, and initial profiles. Main.f - main code with set up of block matrices and calls to LU solvers. Gambc.f, gauss.f, gsolve.f, vambc.f, vsolve.f - BLAS routines for LU decomposition. Fort.3 - input data file. Fort.7 - re-start data file. Header, makefile.

  3. subroutine input include 'header’ character*10 ans,file1 dimension ro(0:n+1),ri(0:n+1),area1(0:n+1) dimension rho1(0:190),u1(0:190),t1(0:190) pi = 4.0d0 * datan(1.0d0) cc refl = 13 cm. and xthroat = 3 cm.c refl = 13.0 xthroat = 3.0/refl cc Normalize dx with respect to lref so that dx is non-dimensionalc dx = 1.0 / float(n+1) dx2 = dx * dx Input.f

  4. do i = 0,n+1 x(i) = float(i) * dx zetax(i) = 1.0 end do do i = 0,n+1 if(x(i).le.xthroat)then area(i) = 4.0 - (2.0*x(i)*refl/3.0) else area(i) = 2.0 + (6.0*(x(i)-xthroat)*refl/10.0) end if end do do i = 0,n+1 write(15,12)i,sngl(x(i)*refl),sngl(area(i)) end do 12 format(i4,2x,2(f10.6,3x)) Input.f

  5. gamma = 5.0/3.0d0 rgas = 1.0d0 / gamma cp = 1.0d0 / (gamma - 1.0d0) read(3,*)nlast read(3,*)tstag read(3,*)pstag read(3,*)pexit read(3,*)tin read(3,*)uin read(3,'(a8)')timestep read(3,*)qin write(6,*)'do you want to read input data ?' read(5,'(a10)')ans if (ans .eq. 'y') then write(6,*)'enter name of input file' read(5,'(a10)')file1 open(unit=7,file=file1,status='unknown') Input.f

  6. do i = 0,n+1 read(7,*)rho(i),u(i),t(i) end do close(unit=7) write(6,*)'reading input data finished' else c open(unit=7,file='initial.d',status='unknown') do i = 0,n+1 t(i) = tstag - (tstag - tin) , * (x(i) / x(n+1)) u(i) = 0.1 + (uin - 0.1) * (x(i) / x(n+1))c u(i) = 0.001d0 + (uin - 0.001d0) c , * (x(i) - x(0)) / (x(n+1) - x(0)) rho(i) = 1.0 - (1.0 - 0.2) , * (x(i) / x(n+1)) write(7,*)t(i),u(i),x(i) end do end ifc close(unit=7) Input.f

  7. do i = 0,n+1 if (x(i) .ge. 2.0d0 .and. x(i) .le. 2.5d0) then c q(i) = qin q(i) = 0.0d0 else q(i) = 0.0d0 end if end do print*,'input',tstag,pstag,tin,u(0)/sqrt(t(0)) return end Input.f

  8. 2000 /no. of time steps 1.0d0 /tstag 1.0d0 /pstag 0.00018 /pexit 0.06d0 /tin 1.7d0 /uin varying 0.0d0 Fort.3

  9. Main.f include 'header' c... mu_art is the artificial viscosity real*8 mu_art c... the arrays d, a, b, and f correspond exactly to the quantities c... D (diagonal block), A (Above diagonal block), B (Below diagonal c... block), and F (right hand side column vector) described in c... lectures 10 and 11 dimension d(nblk,nblk,n),a(nblk,nblk,n), , b(nblk,nblk,n),f(nblk,n) dimension work(nblk,nblk) dimension rho_nm1(0:n+1),u_nm1(0:n+1),t_nm1(0:n+1) call input c... the following dt is up to you to change. dt = 1.0d-3 c... time stepping begins

  10. Main.f do nt = 1,nlast c... the following coefficient of artificial viscosity, mu_art_, is c... also up to you to change. do i = 1,n if(u(i)/sqrt(t(i)).le.1.3)then mu_art = 1.0d0 * rho(i) * dabs(u(i)) * (x(i+1) - x(i)) else cc coefficient changed to 0.1 after 4000 time steps with dt=1.d-3, c and changed to 1.0d-5 after 6000 time steps with same dt c mu_art = 1.0d-5 * rho(i) * dabs(u(i)) * (x(i+1) - x(i)) end if c... Throughout this code, you may ignore the variable zetax(i) - it c... is simply a parameter that enables use of a non-uniform grid, i.e c... unequal deltax. You may ignore it since it is set to 1 in the c... subroutine input.f

  11. Main.f f(1,i) = ((-0.5d0 / (area(i) * dx) * zetax(i)) , * ((rho(i+1) * u(i+1) * area(i+1)) , - (rho(i-1) * u(i-1) * area(i-1))) , + mu_art / dx2 , * (rho(i+1) - 2.0d0 * rho(i) + rho(i-1))) f(2,i) = ((-0.5d0 / (area(i) * dx) * zetax(i)) , * ((rho(i+1) * u(i+1) * area(i+1)) * u(i+1) , - (rho(i-1) * u(i-1) * area(i-1)) * u(i-1)) , - rgas * 0.5d0 / dx * (rho(i+1) * t(i+1) , - rho(i-1) * t(i-1)) * zetax(i) , + mu_art / dx2 * (u(i+1) - 2.0d0 * u(i) + u(i-1))) f(3,i) = ((-0.5d0 / (area(i) * dx) * zetax(i)) * gamma , * ((rho(i+1) * u(i+1) * area(i+1)) * (t(i+1) , + (gamma - 1.0d0) * 0.5d0 * u(i+1) * u(i+1)) , - (rho(i-1) * u(i-1) * area(i-1)) * (t(i-1) , + (gamma - 1.0d0) * 0.5d0 * u(i-1) * u(i-1))) , + mu_art / dx2 * (t(i+1) - 2.0d0 * t(i) + t(i-1)) , + q(i))

  12. Main.f d(1,1,i) = 1.0d0 + dt * mu_art / dx2 d(1,2,i) = 0.0d0 d(1,3,i) = 0.0d0 a(1,1,i) = 0.25d0 * dt / dx * u(i+1) * area(i+1) / area(i) , * zetax(i) , - dt * mu_art / dx2 * 0.5d0 a(1,2,i) = 0.25d0 * dt / dx * rho(i+1) * area(i+1) / area(i) , * zetax(i) a(1,3,i) = 0.0d0 b(1,1,i) = -0.25d0 * dt / dx * u(i-1) * area(i-1) / area(i) , * zetax(i) , - dt * mu_art / dx2 * 0.5d0 b(1,2,i) = -0.25d0 * dt / dx * rho(i-1) * area(i-1) / area(i) , * zetax(i) b(1,3,i) = 0.0d0 f(1,i) = dt * f(1,i)

  13. Main.f d(2,1,i) = u(i) d(2,2,i) = rho(i) , + dt * mu_art / dx2 d(2,3,i) = 0.0d0 a(2,1,i) = 0.25d0 * dt / dx * u(i+1) * u(i+1) , * area(i+1) / area(i) * zetax(i) , + rgas * 0.25d0 * dt / dx * t(i+1) * zetax(i) a(2,2,i) = 0.5d0 * dt / dx * rho(i+1) * u(i+1) , * area(i+1) / area(i) * zetax(i) , - dt * mu_art / dx2 * 0.5d0 a(2,3,i) = rgas * 0.25d0 * dt / dx * rho(i+1) * zetax(i) a(2,3,i) = rgas * 0.25d0 * dt / dx * rho(i+1) * zetax(i) b(2,1,i) = -0.25d0 * dt / dx * u(i-1) * u(i-1) , * area(i-1) / area(i) * zetax(i) , - rgas * 0.25d0 * dt / dx * t(i-1) * zetax(i) b(2,2,i) = -0.5d0 * dt / dx * rho(i-1) * u(i-1) , * area(i-1) / area(i) * zetax(i) , - dt * mu_art / dx2 * 0.5d0 b(2,3,i) = -rgas * 0.25d0 * dt / dx * rho(i-1) * zetax(i) f(2,i) = dt * f(2,i)

  14. Main.f d(3,1,i) = t(i) + gamma * (gamma - 1.0d0) , * 0.5d0 * u(i) * u(i) d(3,2,i) = rho(i) * gamma * (gamma - 1.0d0) * u(i) d(3,3,i) = rho(i) , + dt * mu_art / dx2 a(3,1,i) = 0.25d0 * gamma * dt / dx * u(i+1) , * area(i+1) / area(i) * (t(i+1) , + 0.5d0 * (gamma - 1.0d0) * u(i+1) * u(i+1)) * zetax(i) a(3,2,i) = 0.25d0 * gamma * dt / dx * area(i+1) / area(i) , * rho(i+1) * (t(i+1) + (gamma - 1.0d0) * 1.5d0 , * u(i+1) * u(i+1)) * zetax(i) a(3,3,i) = 0.25d0 * gamma * dt / dx * u(i+1) * rho(i+1) , * area(i+1) / area(i) * zetax(i) , - dt * mu_art / dx2 * 0.5d0 b(3,1,i) = -0.25d0 * gamma * dt / dx * u(i-1) , * area(i-1) / area(i) * (t(i-1) , + 0.5d0 * (gamma - 1.0d0) * u(i-1) * u(i-1)) * zetax(i) b(3,2,i) = -0.25d0 * gamma * dt / dx * area(i-1) / area(i) , * rho(i-1) * (t(i-1) + (gamma - 1.0d0) * 1.5d0 , * u(i-1) * u(i-1)) * zetax(i) b(3,3,i) = -0.25d0 * gamma * dt / dx * u(i-1) * rho(i-1) , * area(i-1) / area(i) * zetax(i) , - dt * mu_art / dx2 * 0.5d0 f(3,i) = dt * f(3,i) end do

  15. Main.f c... boundary conditions c... subsonic inlet c... eliminate density and temperature in favor of velocity cin1 = 1.0d0 / (gamma - 1.0d0) cin2 = (2.0d0 - gamma) / (gamma - 1.0d0) c... The following formula for eliminating B1 by modifying D1 and A1 c... is the same as that given in slide #21 of lecture 11. Note that c... the variable cp is NOT the specific heat. Here, because of the c... non-dimensionalization, cp=1/(gamma-1). Also comparing with the c... formula on slide 21, lecture 11, Pstag & Tstag don't appear c... because they are set equal to 1. do ib = 1,nblk b(ib,2,1) = b(ib,2,1) - u(0) / cp , * (b(ib,3,1) + b(ib,1,1) * cin1 , * t(0) ** cin2) end do c... extrapolate for velocity do ib = 1,nblk d(ib,2,1) = d(ib,2,1) + b(ib,2,1) , * (1.0d0 + (x(1) - x(0)) * zetax(1) / dx) a(ib,2,1) = a(ib,2,1) - b(ib,2,1) , * (x(1) - x(0)) * zetax(1) / dx end do

  16. Main.f c... exit boundary condition. The following lines of code take care c... of automatically choosing the correct BC. By computing the Mach c... number at the exit plane, it is determined if the flow is supersonic c... or subsonic. If supersonic, extrapolation boundary conditions are c... implemented. If subsonic, a pressure specified in fort.3 and read c... in the subroutine input is used in concert with implicit extrapolation c... for u and T or u and rho. exitma = u(n+1) / dsqrt(t(n+1)) if (exitma .gt. 1.0d0) then subsonic = 0.0d0 supsonic = 1.0d0 else subsonic = 1.0d0 supsonic = 0.0d0 end if

  17. Main.f do ib = 1,nblk a(ib,3,n) = a(ib,3,n) - subsonic * a(ib,1,n) , * rho(n+1) / t(n+1) d(ib,1,n) = d(ib,1,n) + supsonic * a(ib,1,n) , * (1.0d0 + (x(n+1) - x(n)) * zetax(n) / dx) b(ib,1,n) = b(ib,1,n) - supsonic * a(ib,1,n) , * (x(n+1) - x(n)) * zetax(n) / dx d(ib,2,n) = d(ib,2,n) + a(ib,2,n) , * (1.0d0 + (x(n+1) - x(n)) * zetax(n) / dx) b(ib,2,n) = b(ib,2,n) - a(ib,2,n) , * (x(n+1) - x(n)) * zetax(n) / dx d(ib,3,n) = d(ib,3,n) + a(ib,3,n) , * (1.0d0 + (x(n+1) - x(n)) * zetax(n) / dx) b(ib,3,n) = b(ib,3,n) - a(ib,3,n) , * (x(n+1) - x(n)) * zetax(n) / dx end do

  18. Main.f c... Now that the block tri-diagonal matrix is set up, BLAS routines are used for c… direct solution by LU decomposition. To minimize storage requirements, the c… solution psi vector is returned as the f vector. call gauss(d(1,1,1),nblk) call gsolve(d(1,1,1),a(1,1,1),nblk) do i = 2,n-1 call gambc(d(1,1,i),b(1,1,i),a(1,1,i-1),work,nblk) call gauss(d(1,1,i),nblk) call gsolve(d(1,1,i),a(1,1,i),nblk) end do call gambc(d(1,1,n),b(1,1,n),a(1,1,n-1),work,nblk) call gauss(d(1,1,n),nblk) c... backward substitution call vsolve(d(1,1,1),f(1,1),nblk) do i = 2,n call vambc(f(1,i),b(1,1,i),f(1,i-1),work,nblk) call vsolve(d(1,1,i),f(1,i),nblk) end do do i = n-1,1,-1 call vambc(f(1,i),a(1,1,i),f(1,i+1),work,nblk) end do

  19. Main.f c... update variables do i = 1,n rho(i) = rho(i) + f(1,i) u(i) = u(i) + f(2,i) t(i) = t(i) + f(3,i) end do c... After updating all the interior points, the boundary points are updated: u(0) = (1.0d0 + (x(1) - x(0)) * zetax(1) / dx) * u(1) - (x(1) - x(0)) * zetax(1) / dx * u(2) t(0) = tstag - 0.5d0 / cp * u(0) * u(0) rho(0) = t(0) ** (1.0d0 / (gamma - 1.0d0)) u(n+1) = (1.0d0 + (x(n+1) - x(n)) * zetax(n) / dx) , * u(n) - (x(n+1) - x(n)) * zetax(n) / dx , * u(n-1) t(n+1) = (1.0d0 + (x(n+1) - x(n)) * zetax(n) / dx) , * t(n) - (x(n+1) - x(n)) * zetax(n) / dx , * t(n-1) rho(n+1) = supsonic * , ((1.0d0 + (x(n+1) - x(n)) * zetax(n) / dx) , * rho(n) - (x(n+1) - x(n)) * zetax(n) / dx , * rho(n-1)) , + subsonic * pexit / t(n+1) if(nt/500*500.eq.nt)write(6,*)nt , ,fmax,dfdtmax,sngl(u(n+1)/dsqrt(t(n+1)))

  20. Main.f c... time stepping ends end do12 format(i3,3x,3(e12.5,2x)) open(unit=7,file='cold.d',status='unknown') do i = 0,n+1 write(7,*)rho(i),u(i),t(i) end do close(unit=7) do i = 0,n+1 write(11,*)i,sngl(rho(i)),u(i)/dsqrt(t(i)),sngl(t(i)) end do 11 stop end

More Related