# Gees - GPL Euler equation solver

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


## 2 Answers

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
• 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. – Azrael3000 Feb 4 '15 at 1:33
• Referenced your comment here: github.com/Azrael3000/gees/issues/1 and already made some progress – Azrael3000 Feb 4 '15 at 2:07
1. Isolate chunks of code that are independent from the rest in separate functions.
2. $$\\pi\$$ is a constant that should be declared as a parameter, where it can be initialized with a call to acos.
3. Use more descriptive self-documenting names for your variables
4. Do not use explicit numbers for i/o units. Use "newunit" instead
5. dble(n) is less cluttery than real(n,8)
6. Avoid magic numbers (e.g., "0.1": what's that?)
7. Reduce the lifespan of variables across the code
8. Do not repeat code
9. Use defensive programming against invalid input

I rewrote the main following most of these criteria, except for the name of some variables, which is not clear from the context what they exactly are. Such variables are apparently fundamental but have names a single-character long (u, v, p, ...), which is far from ideal. As a result of this treatment, the complexity of the code is drastically reduced. Among the other things, the Runge-Kutta propagator could be collapsed to a single step in a loop, instead of being repeated four times, with minimal variations that are easily reproduced by appropriate initialization and exit conditions.

For folks not familiar with the simple syntax of modern Fortran, I suggest quick modern Fortran tutorial.

module mod_constants
real(kind(1d0)), parameter :: PI = acos(-1.d0)
end module mod_constants

program gees
use, intrinsic :: iso_fortran_env, only : OUTPUT_UNIT
use mod_fluxcalc
use mod_constants
implicit none
integer        , parameter :: N_OUTPUT_LISTINGS = 10
real(kind(1d0)), dimension(:)  , allocatable :: p, v
real(kind(1d0)), dimension(:,:), allocatable :: u, f
real(kind(1d0))  :: Time, TimeStep, GridStepSize, SaveTimeInterval, SoundSpeed, DiscretizationSpeed
integer          :: nTimeSteps, iTime, nGridPoints, SaveTimeIterations, ListProgressIterations
logical          :: IsLastIteration, TimeToSave, ListProgress

call PrintCredits()
call FetchParameters( nGridPoints, nTimeSteps, TimeStep, SaveTimeInterval, SoundSpeed)

GridStepSize           = 1.d0 / dble( nGridPoints )
DiscretizationSpeed    = GridStepSize / TimeStep
SaveTimeIterations     = max( nint( SaveTimeInterval / TimeStep ), 1 )
ListProgressIterations = max( nTimeSteps / N_OUTPUT_LISTINGS, 1 )

allocate(p(-1:nGridPoints+1))
allocate(v,source=p)
allocate(f(nGridPoints,2))
allocate(u(-1:nGridPoints+1,2))

call SetInitialConditions( u, nGridPoints )

! set p,v from u and calcuate boundary values at bdry and ghost cells
call bcs        ( u, p, v, nGridPoints, SoundSpeed )
call SaveResults( u, p, v, nGridPoints )

Time = 0.d0
do iTime = 1, nTimeSteps
Time = Time + TimeStep
call RungeKuttaStepper( u, v, f, nGridPoints, SoundSpeed, DiscretizationSpeed )
IsLastIteration = iTime == nTimeSteps
TimeToSave      = mod( iTime, SaveTimeIterations ) == 0
if( TimeToSave .or. IsLastIteration ) call SaveResults( u, p, v, nGridPoints )
ListProgress =  mod( iTime, ListProgressIterations ) == 0
if( ListProgress ) write(OUTPUT_UNIT,"(a,i0,a)") ' Calculated ',Percentage(iTime,nTimeSteps),"%"
enddo

contains

pure integer function Percentage(i,n) result(iRes)
integer, intent(in) :: i, n
iRes = int(dble(i)/dble(n)*100.d0)
end function Percentage

subroutine PrintCredits()
use, intrinsic :: iso_fortran_env, only : OUTPUT_UNIT
character(len=*), parameter  :: FMT="(a)"
write(OUTPUT_UNIT,FMT) ' Welcome to GEES'
write(OUTPUT_UNIT,FMT) ' your friendly GPL Euler Equation solver'
write(OUTPUT_UNIT,FMT) ' written 2012 by Arno Mayrhofer (www.amconception.de)'
write(OUTPUT_UNIT,*)
end subroutine PrintCredits

subroutine RequestParameters( nGridPoints, nTimeSteps, TimeStep, SaveTimeInterval, SoundSpeed )
use, intrinsic :: iso_fortran_env, only : INPUT_UNIT, OUTPUT_UNIT
integer        , intent(out) :: nGridPoints, nTimeSteps
real(kind(1d0)), intent(out) :: TimeStep, SaveTimeInterval, SoundSpeed
character(len=*), parameter  :: FMT="(a)"
real(kind(1d0))              :: TimeEnd
real(kind(1d0)),  parameter  :: MINIMUM_TIME_STEP = 1.d-10
write(OUTPUT_UNIT,FMT,advance="no") ' Number of grid points:'
read (INPUT_UNIT,*) nGridPoints
write(OUTPUT_UNIT,FMT,advance="no") ' End time:'
read (INPUT_UNIT,*) TimeEnd
write(OUTPUT_UNIT,FMT,advance="no") ' Time-step size:'
read (INPUT_UNIT,*) TimeStep
write(OUTPUT_UNIT,FMT,advance="no") ' Output dt:'
read (INPUT_UNIT,*) SaveTimeInterval
write(OUTPUT_UNIT,FMT,advance="no") ' Speed of Sound:'
read (INPUT_UNIT,*) SoundSpeed
if( TimeStep < MINIMUM_TIME_STEP )then
write(OUTPUT_UNIT,"(a,e14.6)") ' Time step too small. Reset to ',MINIMUM_TIME_STEP
TimeStep = MINIMUM_TIME_STEP
endif
nTimeSteps = max(nint(TimeEnd/TimeStep),1)
end subroutine RequestParameters

subroutine SetInitialConditions( Density, n )
use mod_constants
real(kind(1d0)), intent(out) :: Density(:,:)
integer        , intent(in)  :: n
integer :: i
do i = 1, nGridPoints - 1
Density(i,1) = 1.d3 + cos( 2.d0 * PI * dble(i) / dble(nGridPoints) ) * 10d0 !rho
Density(i,2) = 0.d0 !u(i,1)*0.1d0*sin(2*pi*real(i,8)/real(nx,8)) !rho*v
enddo
end subroutine SetInitialConditions

subroutine RungeKuttaStepper( u, v, f, nGridPoints, SoundSpeed, DiscretizationSpeed )
real(kind(1d0)), intent(inout)     :: u(:,:), v(:), f(:,:)
integer        , intent(in)        :: nGridPoints
real(kind(1d0)), intent(in)        :: SoundSpeed, DiscretizationSpeed
integer        , parameter         :: N_RUNGE_KUTTA_STEPS = 4
logical        ,              save :: FIRST_CALL = .TRUE.
real(kind(1d0)), allocatable, save :: utild(:,:), uold(:,:)
integer                            :: i, iStep
if( FIRST_CALL )then
allocate(utild,source=u)
allocate(uold ,source=u)
FIRST_CALL = .FALSE.
endif
uold  = u
utild = u
do iStep = 1, N_RUNGE_KUTTA_STEPS
call fluxcalc(utild,v,f,nGridPoints,SoundSpeed)
do i=1,nGridPoints-1
utild(i,:) = - (f(i+1,:)-f(i,:)) / DiscretizationSpeed
u(i,:)     = u(i,:) + utild(i,:) / 6.d0
utild(i,:) = uold(i,:) + utild(i,:) * 0.5d0
enddo
if( iStep == N_RUNGE_KUTTA_STEPS )exit
call bcs(utild,p,v,nGridPoints,SoundSpeed)
enddo
call bcs(u,p,v,nGridPoints,SoundSpeed)
end subroutine RungeKuttaStepper

subroutine SaveResults( u, p, v, nGridPoints )
implicit none
integer        , intent(in) :: nGridPoints
real(kind(1d0)), intent(in) :: p(:), u(:,:), v(:)
integer, save               :: nCall = 0
character(len=13)           :: FileName
integer                     :: uid, iostat, i
character(len=100)          :: iomsg
write(FileName,'(i5.5,a8)') nCall,'.out.csv'
open(newunit = uid     , &
file    = FileName, &
status  ="unknown", &
iostat  = iostat  , &
iomsg   = iomsg   )
if( iostat /= 0 )then
write(ERROR_UNIT,*) trim(iomsg)
error stop
endif
write(uid,*) 'xpos,p,rho,v'
do i=-1,nGridPoints+1
write(uid,'(4(a1,e20.10))') ' ', dble(i)/dble(nGridPoints),&
',', p(i),',', u(i,1),',', v(i)
enddo
close(uid)
nCall = nCall + 1
end subroutine SaveResults

end program gees