Gees - GPL Euler equation solver

Question

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

Answer 1

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

• Azrael3000: 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: Referenced your comment here: github.com/Azrael3000/gees/issues/1 and already made some progress

Answer 2

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