As a little helper I recently had to write a code that solves the 1-D Euler equations. As it serves my purpose well I though others could make use of it as well. The homepage of the code can be found here, including a download containing the Makefile.
This version of the code is in this Github repository.
File gees.f90:
program gees use mod_fluxcalc implicit none real(8), dimension(:),allocatable :: p,v real(8), dimension(:,:), allocatable :: u,f, utild, uold real(8) :: temps, dt, tend, dx, odt, pi, c0 integer :: i, nt, it, nx, io, lout character(len=13) :: fname write(*,*) 'Welcome to GEES' write(*,*) 'your friendly GPL Euler Equation solver' write(*,*) 'written 2012 by Arno Mayrhofer (www.amconception.de)' write(*,*) write(*,*) 'Number of grid points:' read(*,*) nx write(*,*) nx write(*,*) 'End time:' read(*,*) tend write(*,*) tend write(*,*) 'Time-step size' read(*,*) dt write(*,*) dt write(*,*) 'Output dt:' read(*,*) odt write(*,*) odt write(*,*) 'Speed of sound:' read(*,*) c0 write(*,*) c0 ! number of time-steps nt = int(tend/dt+1d-14) ! grid size dx = 1D0/real(nx,8) pi = acos(-1d0) ! list ouput index lout = 1 allocate(p(-1:nx+1),v(-1:nx+1),f(nx,2),u(-1:nx+1,2), & utild(-1:nx+1,2), uold(-1:nx+1,2)) ! init do i=1,nx-1 u(i,1) = 1d3+cos(2*pi*real(i,8)/real(nx,8))*10d0 !rho u(i,2) = 0d0 !u(i,1)*0.1d0*sin(2*pi*real(i,8)/real(nx,8)) !rho*v enddo ! set p,v from u and calcuate boundary values at bdry and ghost cells call bcs(u,p,v,nx,c0) ! file output index io = 0 ! simulation time temps = 0d0 ! output write(fname,'(i5.5,a8)') io,'.out.csv' open(file=fname,unit=800) write(800,*) 'xpos,p,rho,v' do i=-1,nx+1 write(800,'(4(a1,e20.10))') ' ', real(i,8)/real(nx,8),',', p(i),',', u(i,1),',', v(i) enddo close(800) io = io+1 ! loop over all timesteps do it=1,nt ! list output if(real(nt*lout)*0.1d0.lt.real(it))then write(*,*) 'Calculated ', int(real(lout)*10.), '%' lout = lout + 1 endif temps = temps + dt uold = u ! First Runge-Kutta step ! calculate flux at mid points using u call fluxcalc(u,v,f,nx,c0) do i=1,nx-1 ! calc k1 = -dt/dx(f(u,i+1/2) - f(u,i-1/2)) utild(i,:) = -dt/dx*(f(i+1,:)-f(i,:)) ! u^n = uold^n + 1/6 k1 u(i,:) = uold(i,:) + utild(i,:)/6d0 ! utild = uold^n + 1/2 k_1 utild(i,:) = uold(i,:) + utild(i,:)*0.5d0 enddo ! calculate p,v + bcs for uold + 1/2 k1 = utild call bcs(utild,p,v,nx,c0) ! Second Runge-Kutta step ! calculate flux at mid points using utild call fluxcalc(utild,v,f,nx,c0) do i=1,nx-1 ! calc k2 = -dt/dx(f(utild,i+1/2) - f(utild,i-1/2)) utild(i,:) = -dt/dx*(f(i+1,:)-f(i,:)) ! u^n = u^n + 1/6 k2 u(i,:) = u(i,:) + utild(i,:)/6d0 ! utild = uold^n + 1/2 k_2 utild(i,:) = uold(i,:) + utild(i,:)*0.5d0 enddo ! calculate p,v + bcs for uold + 1/2 k2 = utild call bcs(utild,p,v,nx,c0) ! Third Runge-Kutta step ! calculate flux at mid points using utild call fluxcalc(utild,v,f,nx,c0) do i=1,nx-1 ! calc k3 = -dt/dx(f(utild,i+1/2) - f(utild,i-1/2)) utild(i,:) = -dt/dx*(f(i+1,:)-f(i,:)) ! u^n = u^n + 1/6 k3 u(i,:) = u(i,:) + utild(i,:)/6d0 ! utild = uold^n + k_3 utild(i,:) = uold(i,:) + utild(i,:) enddo ! calculate p,v + bcs for uold + k3 = utild call bcs(utild,p,v,nx,c0) ! Fourth Runge-Kutta step ! calculate flux at mid points using utild call fluxcalc(utild,v,f,nx,c0) do i=1,nx-1 ! calc k4 = -dt/dx(f(utild,i+1/2) - f(utild,i-1/2)) utild(i,:) = -dt/dx*(f(i+1,:)-f(i,:)) ! u^n = u^n + 1/6 k4 u(i,:) = u(i,:) + utild(i,:)/6d0 enddo ! calculate p,v + bcs for (uold + k1 + k2 + k3 + k4)/6 = u call bcs(u,p,v,nx,c0) ! output if(abs(temps-odt*real(io,8)).lt.abs(temps+dt-odt*real(io,8)).or.it.eq.nt)then write(fname,'(i5.5,a8)') io,'.out.csv' open(file=fname,unit=800) write(800,*) 'xpos,p,rho,v' do i=-1,nx+1 write(800,'(4(a1,e20.10))') ' ', real(i,8)/real(nx,8),',', p(i),',', u(i,1),',', v(i) enddo close(800) io = io + 1 endif enddo end program gees fluxcalc.f90:
module mod_fluxcalc contains subroutine fluxcalc(u, v, f, nx, c0) implicit none integer, intent(in) :: nx real(8), dimension(-1:nx+1), intent(inout) :: v real(8), dimension(-1:nx+1,2), intent(inout) :: u real(8), dimension(nx,2), intent(inout) :: f real(8), intent(in) :: c0 integer :: i, j real(8), dimension(2) :: ur, ul real(8) :: a, pl, pr, nom, denom, ri, rim1 real(8), dimension(4) :: lam ! calculate flux at midpoints, f(i) = f_{i-1/2} do i=1,nx do j=1,2 ! calculate r_{i} = \frac{u_i-u_{i-1}}{u_{i+1}-u_i} nom = (u(i ,j)-u(i-1,j)) denom = (u(i+1,j)-u(i ,j)) ! make sure division by 0 does not happen if(abs(nom).lt.1d-14)then ! nom = 0 nom = 0d0 denom = 1d0 elseif(nom.gt.1d-14.and.abs(denom).lt.1d-14)then ! nom > 0 => r = \inf nom = 1d14 denom = 1d0 elseif(nom.lt.-1d-14.and.abs(denom).lt.1d-14)then ! nom < 0 => r = 0 nom = -1d14 denom = 1d0 endif ri = nom/denom ! calculate r_{i-1} = \frac{u_{i-1}-u_{i-2}}{u_i-u_{i-1}} nom = (u(i-1,j)-u(i-2,j)) denom = (u(i ,j)-u(i-1,j)) ! make sure division by 0 does not happen if(abs(nom).lt.1d-14)then nom = 0d0 denom = 1d0 elseif(nom.gt.1d-14.and.abs(denom).lt.1d-14)then nom = 1d14 denom = 1d0 elseif(nom.lt.-1d-14.and.abs(denom).lt.1d-14)then nom = -1d14 denom = 1d0 endif rim1 = nom/denom ! u^l_{i-1/2} = u_{i-1} + 0.5*phi(r_{i-1})*(u_i-u_{i-1}) ul(j) = u(i-1,j)+0.5d0*phi(rim1)*(u(i ,j)-u(i-1,j)) ! u^r_{i-1/2} = u_i + 0.5*phi(r_i)*(u_{i+1}-u_i) ur(j) = u(i,j) -0.5d0*phi(ri )*(u(i+1,j)-u(i ,j)) enddo ! calculate eigenvalues of \frac{\partial F}{\parial u} at u_{i-1} lam(1) = ev(v,u,c0,i-1, 1d0,nx) lam(2) = ev(v,u,c0,i-1,-1d0,nx) ! calculate eigenvalues of \frac{\partial F}{\parial u} at u_i lam(3) = ev(v,u,c0,i , 1d0,nx) lam(4) = ev(v,u,c0,i ,-1d0,nx) ! max spectral radius (= max eigenvalue of dF/du) of flux Jacobians (u_i, u_{i-1}) a = maxval(abs(lam),dim=1) ! calculate pressure via equation of state: ! p = \frac{rho_0 c0^2}{xi}*((\frac{\rho}{\rho_0})^xi-1), (xi=7, rho_0=1d3) pr = 1d3*c0**2/7d0*((ur(1)/1d3)**7d0-1d0) pl = 1d3*c0**2/7d0*((ul(1)/1d3)**7d0-1d0) ! calculate flux based on the Kurganov and Tadmor central scheme ! F_{i-1/2} = 0.5*(F(u^r_{i-1/2})+F(u^l_{i-1/2}) - a*(u^r_{i-1/2} - u^l_{i-1/2})) ! F_1 = rho * v = u_2 f(i,1) = 0.5d0*(ur(2)+ul(2)-a*(ur(1)-ul(1))) ! F_2 = p + rho * v**2 = p + \frac{u_2**2}{u_1} f(i,2) = 0.5d0*(pr+ur(2)**2/ur(1)+pl+ul(2)**2/ul(1)-a*(ur(2)-ul(2))) enddo end subroutine ! flux limiter real(8) function phi(r) implicit none real(8), intent(in) :: r phi = 0d0 ! ospre flux limiter phi(r) = \frac{1.5*(r^2+r)}{r^2+r+1} if(r.gt.0d0)then phi = 1.5d0*(r**2+r)/(r**2+r+1d0) endif ! van leer ! phi = (r+abs(r))/(1d0+abs(r)) end function ! eigenvalue calc real(8) function ev(v,u,c0,i,sgn,nx) implicit none real(8), dimension(-1:nx+1), intent(inout) :: v real(8), dimension(-1:nx+1,2), intent(inout) :: u integer, intent(in) :: i, nx real(8), intent(in) :: sgn, c0 ! calculate root of characteristic equation ! \lambda = 0.5*(3*v \pm sqrt{5v^2+4c0^2(\frac{\rho}{\rho_0})^{xi-1}}) ev = 0.5d0*(3d0*v(i)+sgn*sqrt(5d0*v(i)**2+4d0*c0**2*(u(i,1)/1d3)**6)) return end function subroutine bcs(u,p,v,nx,c0) implicit none integer, intent(in) :: nx real(8), intent(in) :: c0 real(8), dimension(-1:nx+1), intent(inout) :: p,v real(8), dimension(-1:nx+1,2), intent(inout) :: u integer :: i ! calculate velocity and pressure do i=1,nx-1 ! v = u_2 / u_1 v(i) = u(i,2)/u(i,1) ! p = \frac{rho_0 c0^2}{xi}*((\frac{\rho}{\rho_0})^xi-1), (xi=7, rho_0=1d3) p(i) = 1d3*c0**2/7d0*((u(i,1)/1d3)**7d0-1d0) enddo ! calculate boundary conditions ! note: for periodic boundary conditions set ! f(0) = f(nx-1), f(-1) = f(nx-2), f(nx) = f(1), f(nx+1) = f(2) \forall f ! for rho: d\rho/dn = 0 using second order extrapolation u(0,1) = (4d0*u(1,1)-u(2,1))/3d0 u(-1,1) = u(1,1) u(nx,1) = (4d0*u(nx-1,1)-u(nx-2,1))/3d0 u(nx+1,1) = u(nx-1,1) ! for p: dp/dn = 0 (thus can use EOS) p(0) = 1d3*c0**2/7d0*((u(0,1)/1d3)**7d0-1d0) p(-1) = 1d3*c0**2/7d0*((u(-1,1)/1d3)**7d0-1d0) p(nx) = 1d3*c0**2/7d0*((u(nx,1)/1d3)**7d0-1d0) p(nx+1) = 1d3*c0**2/7d0*((u(nx+1,1)/1d3)**7d0-1d0) ! for v: v = 0 using second order extrapolation v(0) = 0d0 v(-1) = v(2)-3d0*v(1) v(nx) = 0d0 v(nx+1) = v(nx-2)-3d0*v(nx-1) ! for rho*v u(0,2) = u(0,1)*v(0) u(-1,2) = u(-1,1)*v(-1) u(nx,2) = u(nx,1)*v(nx) u(nx+1,2) = u(nx+1,1)*v(nx+1) end subroutine end module The code is covered in comments as it should allow the beginner to understand the inner workings fairly easily. If you have any comments on what could be done differently let me know. I hope this code is bug free, but in case there is still one in the hiding tell me and I'll fix it.
Implementation details: Time-stepping: 4th order explicit Runge Kutta Flux calculation: MUSCL scheme (Kurganov and Tadmor central scheme) Flux limiter: Ospre Boundary conditions: Second order polynomial extrapolation Equation of state: Tait
