4
\$\begingroup\$

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
\$\endgroup\$
4
\$\begingroup\$

A few points I'd do different:

  • use the name of the procedures/modules also in their end statement
  • not writing enddo and endif in a single word, but seperately. There are newer language constructs, where writing them together is not allowed, and separating them all is more consistent
  • as mentioned in a comment above, use selected_real_kind instead of a hardwired 8
  • use >,<, == instead of .gt., .lt. and .eq.
  • provide some explanation on routines arguments and the purpose of each routine
  • indent module routines
  • use spaces around the condition of if-clauses
  • use the private statement in the module and explicitly mark visible entities with the public keyword
  • turn the initialization and the time loop into subroutines each
\$\endgroup\$
  • \$\begingroup\$ Thank you for your comments. I recently got some feedback via email and hope to incorporate also your suggestions. The code is now on github so if you want to contribute feel free to send pull requests. \$\endgroup\$ – Azrael3000 Feb 4 '15 at 1:33
  • \$\begingroup\$ Referenced your comment here: github.com/Azrael3000/gees/issues/1 and already made some progress \$\endgroup\$ – Azrael3000 Feb 4 '15 at 2:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.