! ! This fortran 90 program is used in MATH 365 - Computational ! Fluid Dynamics at James Madison University. This program was ! written by James Sochacki of the Department of Mathematics ! at James Madison University as part of the NSF Grant - A Collaborative ! Compuatational Science Program between JMU and North Carolina Central ! University. His info is ! ! Jim Sochacki ! Department of Mathematics ! James Madison University ! Harrisonburg, VA 22807 ! jim@math.jmu.edu program euler_1d_lf ! This code calculates the leap frog approximations for the ! density, velocity and energy to the one dimensional Euler equations of ! gas dynamics. The equations are those given in the module intro_cfd.mod implicit none real :: rho_left,rho_right,u_left,u_right,E_left,E_right real :: s_u,s_rho,s_E real :: dx,dt,L,x,gamma,time real,dimension(:),allocatable :: m0,m1,m2 real,dimension(:),allocatable :: rho0,rho1,rho2 real,dimension(:),allocatable :: E0,E1,E2 integer :: count,interval,tstep,jmax integer :: output_type,i,j character*6 type_d,type_v,type_E character*11 filename print*,' This code computes the density, velocity and energy for ' print*,' the one dimensional nonlinear Euler equations using ' print*,' leap frog. ' print*,' Give the gas constant, gamma ' read*,gamma print*,' Give u_left the value of u(x,0) at -infinity ' read*,u_left print*,' Give u_right the value of u(x,0) at -infinity ' read*,u_right print*,' Give s_u the "slope" of u(0,0) ' read*,s_u print*,' Give rho_left the value of rho(x,0) at -infinity ' read*,rho_left print*,' Give rho_right the value of rho(x,0) at -infinity ' read*,rho_right print*,' Give s_rho the "slope" of rho(0,0) ' read*,s_rho print*,' Give E_left the value of E(x,0) at -infinity ' read*,E_left print*,' Give E_right the value of E(x,0) at -infinity ' read*,E_right print*,' Give s_E the "slope" of E(0,0) ' read*,s_E print*,' Give the time to run the model ' read*,time print*,' Give the time step ' read*, dt print*,' Give L ' read*,L print*,' Give the grid size ' read*, dx print*,' Type a 0 for density output, a 1 for velocity output and ' print*,' a 2 for energy output and a 3 or larger for all to be ' print*,' outputted.' read*,output_type if (output_type==0 .or. output_type>=3) then type_d = 'densty' end if if (output_type==1 .or. output_type>=3) then type_v = 'velcty' end if if (output_type>=2) then type_E = 'energy' end if print*,' Give an integer indicating how often a snapshot ' print*,' of the profile should be taken. ' read*,interval ! jmax is the number of grid points in the model jmax = 2*L/dx tstep=time/dt ! Set up the memory allocate(m0(0:jmax+1)) allocate(m1(0:jmax+1)) allocate(m2(0:jmax+1)) allocate(rho0(0:jmax+1)) allocate(rho1(0:jmax+1)) allocate(rho2(0:jmax+1)) allocate(E0(0:jmax+1)) allocate(E1(0:jmax+1)) allocate(E2(0:jmax+1)) x = -L ! ! Set up the initial values for the leap frog scheme. ! do j = 0,jmax+1 x = x+dx rho0(j) = p(rho_left,rho_right,s_rho,x) m0(j) = p(rho_left,rho_right,s_rho,x)*p(u_left,u_right,s_u,x) E0(j) = p(E_left,E_right,s_E,x) end do E1(0) = E0(0) E1(jmax+1) = E0(jmax+1) rho1(0) = rho0(0) rho1(jmax+1) = rho0(jmax+1) m1(0) = m0(0) m1(jmax+1) = m0(jmax+1) ! ! Start the calculations. ! tstep is the number of time steps ! We initialize the first step in leap frog using forward euler in time do j = 1,jmax m1(j) = m0(j)+(dt/(2*dx))*((1-gamma)*(E0(j+1)-E0(j-1))+ & ((gamma-3)/2)*(m0(j+1)**2/rho0(j+1)-m0(j-1)**2/rho0(j-1))) rho1(j) = rho0(j)-(dt/(2*dx))*(m0(j+1)-m0(j-1)) E1(j) = E0(j)+(dt/(2*dx))*(-gamma*(E0(j+1)*m0(j+1)/rho0(j+1)- & E0(j-1)*m0(j-1)/rho0(j-1))+ & ((gamma-1)/2)* & (m0(j+1)**3/rho0(j+1)**2-m0(j-1)**3/rho0(j-1)**2)) end do count = 0 i = 0 if (output_type==0 .or. output_type>=3) then write (filename,'(a6,i5.5)') type_d,i open (6,file=filename,status='new') do j=1,jmax write (6,'(f14.8)') rho2(j) end do close (6) end if if (output_type==1 .or. output_type>=3) then write (filename,'(a6,i5.5)') type_v,i open (6,file=filename,status='new') do j=1,jmax write (6,'(f14.8)') m2(j)/rho2(j) end do close (6) end if if (output_type>=2) then write (filename,'(a6,i5.5)') type_E,i open (6,file=filename,status='new') do j=1,jmax write (6,'(f14.8)') E2(j) end do close (6) end if do i=1,tstep count=count+1 if (count==interval) print*,' step = ',i do j = 1 , jmax m2(j) = m0(j)+(dt/dx)*((1-gamma)*(E1(j+1)-E1(j-1))+ & ((gamma-3)/2)*(m1(j+1)**2/rho1(j+1)-m1(j-1)**2/rho1(j-1))) rho2(j) = rho0(j)-(dt/dx)*(m1(j+1)-m1(j-1)) E2(j) = E0(j)+(dt/dx)*(-gamma*(E1(j+1)*m1(j+1)/rho1(j+1)- & E1(j-1)*m1(j-1)/rho1(j-1))+ & ((gamma-1)/2)* & (m1(j+1)**3/rho1(j+1)**2-m1(j-1)**3/rho1(j-1)**2)) end do ! Create the snapshots. if (count==interval)then count = 0 if (output_type==0 .or. output_type>=3) then write (filename,'(a6,i5.5)') type_d,i open (6,file=filename,status='new') do j=1,jmax write (6,'(f14.8)') rho2(j) end do close (6) end if if (output_type==1 .or. output_type>=3) then write (filename,'(a6,i5.5)') type_v,i open (6,file=filename,status='new') do j=1,jmax write (6,'(f14.8)') m2(j)/rho2(j) end do close (6) end if if (output_type>=2) then write (filename,'(a6,i5.5)') type_E,i open (6,file=filename,status='new') do j=1,jmax write (6,'(f14.8)') E2(j) end do close (6) end if end if ! ! Update the calculations in time ! do j = 1 , jmax m0(j)=m1(j) rho0(j)=rho1(j) E0(j)=E1(j) m1(j)=m2(j) rho1(j)=rho2(j) E1(j)=E2(j) end do end do if (output_type==0 .or. output_type>=3) then print*,' the snapshots are in file: ',type_d end if if (output_type==1 .or. output_type>=3) then print*,' the snapshots are in file: ',type_v end if if (output_type>=2) then print*,' the snapshots are in file: ',type_E end if contains function p(a_l,a_r,s,x) real :: x,a_l,a_r,s,p,q,r r = (a_l+a_r)/2; q = (a_r-a_l)/(4*atan(1.0)); p = q*atan(s*x)+r return end function end program euler_1d_lf