I have written a simple Euler solver for the 1D shock tube problem. Eventually, I plan to extend this code to solve the full 3D compressible Navier-Stokes equations. Therefore, I want to start with good programming practices as it will be more difficult to modify the code further down the road. The idea is to find good balance between readability and performance.
My thoughts so far:
- Should I use Fortran structures?
- Functions vs. subroutines: I have found (tested) that functions can be a bit less efficient than subroutine (~5%), but I do like them better as it is always clear to the reader which variable was modified upon a procedure call. The problem is down the road, 5% can mean an extra week of running time.
- Where should I store local variables? A few choices: declare local variables in procedures, store all of them in module declaration. I am not considering global variables - with the exception of
PARAMETERS
- as I do not consider them a good idea, even though I have read it can boost performance.
I am of course open to any other recommendation/advice.
Makefile:
FC = gfortran
FLAGS = -Wall -Wtabs
SPE = main.f90
SRCS = global_vars.f90 initProblem.f90 AUSMmethod.f90 WENOmethod.f90 file_io.f90 main.f90
SOBJ = $(SRCS:.f90=.o)
EXEC = $(SPE:.f90=)
all: $(EXEC)
touch $*.o $*.mod
$(EXEC): $(SOBJ)
$(FC) $(FLAGS) -o executable $^
%.o: %.f90
$(FC) $(FLAGS) -c $<
clean:
rm -rf *o *mod executable
Fortran files:
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! main.f90
PROGRAM main
USE global_vars, ONLY: sp,dp, Nx
USE initProblem, ONLY: get_Xgrid, set_initialConditions
USE AUSMmethod, ONLY: AUSMscheme
USE file_io, ONLY: write_solution
IMPLICIT NONE
REAL(dp), DIMENSION(:), ALLOCATABLE :: x ! x grid locations
REAL(dp), DIMENSION(:), ALLOCATABLE :: p ! Pressure
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ucon ! Conservative vector
! Initialize variables
ALLOCATE(x (0:Nx-1))
ALLOCATE(p (0:Nx-1))
ALLOCATE(Ucon(0:Nx-1,0:2))
CALL get_Xgrid(x)
CALL set_initialConditions(x,p,Ucon)
CALL AUSMscheme(p,Ucon)
CALL write_solution(x,p,Ucon(:,0),Ucon(:,1)/Ucon(:,0))
END PROGRAM main
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! global_vars.f90
MODULE global_vars
USE, INTRINSIC:: ISO_FORTRAN_ENV, sp => real32, dp => real64
IMPLICIT NONE
! Grid Parameters
INTEGER, PARAMETER :: Nx = 500 ! Number of points
REAL(dp), PARAMETER :: x_o = 0._dp ! Lower x-boundary
REAL(dp), PARAMETER :: x_f = 1._dp ! Upper x-boundary
REAL(dp), PARAMETER :: dx = (x_f - x_o)/(Nx-1) ! x grid spacing
! Physical Parameters
REAL(dp), PARAMETER :: mGamma = 1.4_dp ! Specific Heat ratio
REAL(dp), PARAMETER :: CFL = 0.5_dp ! CFL number for stability
REAL(dp), PARAMETER :: t_max = 0.1452_dp ! Maximum time at which simulation is
END MODULE global_vars
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! initProblem.f90
MODULE initProblem
USE global_vars, ONLY: dp, Nx, x_o, x_f, dx, mGamma
IMPLICIT NONE
CONTAINS
! Initialize computational grid
PURE SUBROUTINE get_Xgrid(x_grid)
REAL(dp), DIMENSION(0:), INTENT(OUT) :: x_grid
INTEGER :: ii
DO ii=0,Nx-1
x_grid(ii) = x_o + ii*dx
END DO
END SUBROUTINE get_Xgrid
! Setup initial conditions of the 1D shock tube problem
PURE SUBROUTINE set_initialConditions(x,p,Ucon)
REAL(dp), DIMENSION(0:), INTENT(IN) :: x
REAL(dp), DIMENSION(0:), INTENT(OUT) :: p
REAL(dp), DIMENSION(0:,0:), INTENT(OUT) :: Ucon
WHERE(x <= 0.5)
p = 1._dp
Ucon(:,0) = 1._dp
ELSE WHERE
p = 0.1_dp
Ucon(:,0) = 0.125_dp
END WHERE
Ucon(:,1) = 0._dp
Ucon(:,2) = p/(mGamma-1) + 0.5_dp*Ucon(:,1)**2/Ucon(:,0)
END SUBROUTINE set_initialConditions
END MODULE initProblem
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AUSMmethod.f90
MODULE AUSMmethod
USE global_vars, ONLY: dp, Nx, dx, mGamma, CFL, t_max
IMPLICIT NONE
CONTAINS
! Applies boundary coniditon on conservative vector
PURE SUBROUTINE apply_bc(Ucon)
REAL(dp), DIMENSION(0:,0:), INTENT(INOUT) :: Ucon
Ucon( 0,:) = Ucon(1,:)
Ucon(Nx-1,:) = Ucon(Nx-2,:)
END SUBROUTINE apply_bc
! Based on the sound and convective speed, returns the time step required for stability
PURE SUBROUTINE get_TimeStep(sound,speed,dt)
REAL(dp), DIMENSION(0:), INTENT(IN) :: sound
REAL(dp), DIMENSION(0:), INTENT(IN) :: speed
REAL(dp), INTENT(OUT) :: dt
dt = CFL*dx/MAXVAL(sound + ABS(speed))
END SUBROUTINE get_TimeStep
! Based on the pressure and density, returns the Mach number
PURE SUBROUTINE get_sound(pressure,density,sound)
REAL(dp), DIMENSION(0:), INTENT(IN) :: density
REAL(dp), DIMENSION(0:), INTENT(IN) :: pressure
REAL(dp), DIMENSION(0:), INTENT(OUT) :: sound
sound = SQRT(mGamma*pressure/density)
END SUBROUTINE get_sound
! Based on the pressure and density, returns the Mach number
PURE SUBROUTINE get_Mach(sound,velocity,Mach)
REAL(dp), DIMENSION(0:), INTENT(IN) :: sound
REAL(dp), DIMENSION(0:), INTENT(IN) :: velocity
REAL(dp), DIMENSION(0:), INTENT(OUT) :: Mach
Mach = velocity/sound
END SUBROUTINE get_Mach
! Based on the density, velocity and total energy, returns the pressure
PURE SUBROUTINE get_Pressure(density,velocity,totEnergy,pressure)
REAL(dp), DIMENSION(0:), INTENT(IN) :: density
REAL(dp), DIMENSION(0:), INTENT(IN) :: velocity
REAL(dp), DIMENSION(0:), INTENT(IN) :: totEnergy
REAL(dp), DIMENSION(0:), INTENT(OUT) :: pressure
pressure = (mGamma - 1) * (totEnergy - 0.5_dp*density*velocity**2)
END SUBROUTINE get_Pressure
! Based on the Mach number at neighboring points, returns the corresponding Mach number at the
! interface
PURE SUBROUTINE get_Mhalf(Mach,Mhalf)
REAL(dp), DIMENSION(0:), INTENT(IN) :: Mach
REAL(dp), DIMENSION(0:), INTENT(OUT) :: Mhalf
REAL(dp), DIMENSION(:), ALLOCATABLE :: M_plus
REAL(dp), DIMENSION(:), ALLOCATABLE :: M_minus
ALLOCATE(M_plus (0:Nx-1))
ALLOCATE(M_minus(0:Nx-1))
! Depending on the flow conditions (subsonic vs supersonic), use different formulas
WHERE(ABS(Mach) <= 1)
M_plus = 0.25_dp*(Mach + 1)**2
M_minus = -0.25_dp*(Mach - 1)**2
ELSE WHERE
M_plus = 0.50_dp*(Mach + ABS(Mach))
M_minus = 0.50_dp*(Mach - ABS(Mach))
END WHERE
Mhalf(0:Nx-2) = M_plus(0:Nx-2) + M_minus(1:Nx-1)
DEALLOCATE(M_plus)
DEALLOCATE(M_minus)
END SUBROUTINE get_Mhalf
! Based on the Mach number and pressure, returns the pressure flux
PURE SUBROUTINE get_Pflux(p,Mach,Pflux)
REAL(dp), DIMENSION(0:), INTENT(IN) :: p
REAL(dp), DIMENSION(0:), INTENT(IN) :: Mach
REAL(dp), DIMENSION(0:), INTENT(OUT) :: Pflux
REAL(dp), DIMENSION(:), ALLOCATABLE :: Pplus
REAL(dp), DIMENSION(:), ALLOCATABLE :: Pminus
ALLOCATE(Pplus (0:Nx-1))
ALLOCATE(Pminus(0:Nx-1))
! Depending on the flow conditions (subsonic vs supersonic), use different formulas
WHERE(ABS(Mach) <= 1)
Pplus = 0.25_dp* p * (Mach+1)**2 * (2-Mach)
Pminus = 0.25_dp* p * (Mach-1)**2 * (2+Mach)
ELSE WHERE
Pplus = 0.50_dp* p * (Mach+ABS(Mach))/Mach
Pminus = 0.50_dp* p * (Mach-ABS(Mach))/Mach
END WHERE
Pflux(0:Nx-2) = Pplus(0:Nx-2) + Pminus(1:Nx-1)
DEALLOCATE(Pplus )
DEALLOCATE(Pminus)
END SUBROUTINE get_Pflux
! Based on the sound speed, Mach number, pressure and conservative variables, returns the total
! x-directional flux
PURE SUBROUTINE get_AUSMflux(sound,Mach,p,Ucon,AUSMflux)
REAL(dp), DIMENSION(0:), INTENT(IN) :: sound
REAL(dp), DIMENSION(0:), INTENT(IN) :: Mach
REAL(dp), DIMENSION(0:), INTENT(IN) :: p
REAL(dp), DIMENSION(0:,0:), INTENT(IN) :: Ucon
REAL(dp), DIMENSION(0:,0:), INTENT(OUT) :: AUSMflux
REAL(dp), DIMENSION(:), ALLOCATABLE :: M_half
REAL(dp), DIMENSION(:), ALLOCATABLE :: Pflux
! Initialize variables
ALLOCATE(M_half(0:Nx-2))
ALLOCATE(Pflux (0:Nx-2))
! Calculate the speed of sound and modified Mach number
CALL get_Mhalf(Mach,M_half)
CALL get_Pflux(p,Mach,Pflux)
! Calculate Flux
WHERE(M_half >= 0)
AUSMflux(0:Nx-2,0) = M_half(0:Nx-2) * sound(0:Nx-2) * Ucon(0:Nx-2,0)
AUSMflux(0:Nx-2,1) = M_half(0:Nx-2) * sound(0:Nx-2) * Ucon(0:Nx-2,1) + Pflux(0:Nx-2)
AUSMflux(0:Nx-2,2) = M_half(0:Nx-2) * sound(0:Nx-2) *(Ucon(0:Nx-2,2) + p (0:Nx-2))
ELSE WHERE
AUSMflux(0:Nx-2,0) = M_half(0:Nx-2) * sound(1:Nx-1) * Ucon(1:Nx-1,0)
AUSMflux(0:Nx-2,1) = M_half(0:Nx-2) * sound(1:Nx-1) * Ucon(1:Nx-1,1) + Pflux(1:Nx-1)
AUSMflux(0:Nx-2,2) = M_half(0:Nx-2) * sound(1:Nx-1) *(Ucon(1:Nx-1,2) + p (1:Nx-1))
END WHERE
! Deallocate variables
DEALLOCATE(M_half)
DEALLOCATE(Pflux)
END SUBROUTINE get_AUSMflux
! Main subroutine: applies the AUSM scheme up to a given time
SUBROUTINE AUSMscheme(p,Ucon)
REAL(dp), DIMENSION(0:), INTENT(INOUT) :: p
REAL(dp), DIMENSION(0:,0:), INTENT(INOUT) :: Ucon
INTEGER :: tt
REAL(dp) :: time, dt
REAL(dp), DIMENSION(:), ALLOCATABLE :: Mach
REAL(dp), DIMENSION(:), ALLOCATABLE :: sound
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Uflux
ALLOCATE(Mach (0:Nx-1))
ALLOCATE(sound(0:Nx-1))
ALLOCATE(Uflux (0:Nx-2,0:2))
! Apply numerical scheme until tmax is reached
time = 0._dp; tt = 0;
DO WHILE(time <= t_max)
! Caluclate fluid variables
CALL get_Pressure(Ucon(:,0),Ucon(:,1)/Ucon(:,0),Ucon(:,2),p)
CALL get_sound(p,Ucon(:,0),sound)
CALL get_Mach (sound,Ucon(:,1)/Ucon(:,0),Mach)
! Calculate required time step depending on maximum wave speed
CALL get_TimeStep(sound,Ucon(:,1)/Ucon(:,0),dt)
! Get fluxes
CALL get_AUSMflux(sound,Mach,p,Ucon,Uflux)
! Calculate the the conservative vector for the new time step
Ucon(1:Nx-2,:) = Ucon(1:Nx-2,:) - dt/dx * (Uflux(1:Nx-2,:) - Uflux(0:Nx-3,:))
! Increment the time and keep the time step count
time = time + dt; tt = tt +1;
END DO
PRINT '(A12,F8.5)', "t = ", time
PRINT '(A12,I3)', "Time step = ", tt
END SUBROUTINE AUSMscheme
END MODULE AUSMmethod
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! file_io.f90
MODULE file_io
USE global_vars, ONLY: dp, Nx
IMPLICIT NONE
CONTAINS
! Write solution to output file
SUBROUTINE write_solution(x,var1,var2,var3,var4)
REAL(dp), DIMENSION(0:), INTENT(IN) :: x
REAL(dp), DIMENSION(0:), INTENT(IN) :: var1
REAL(dp), DIMENSION(0:), INTENT(IN) :: var2
REAL(dp), DIMENSION(0:), INTENT(IN) :: var3
REAL(dp), DIMENSION(0:), INTENT(IN), OPTIONAL :: var4
INTEGER :: ii
OPEN(UNIT=100,FILE='output/output.dat',STATUS='REPLACE')
IF(PRESENT(var4)) THEN
DO ii=0,Nx-1
WRITE(100,*) x(ii), var1(ii), var2(ii), var3(ii), var4(ii)
END DO
ELSE
DO ii=0,Nx-1
WRITE(100,*) x(ii), var1(ii), var2(ii), var3(ii)
END DO
END IF
CLOSE(100)
END SUBROUTINE write_solution
END MODULE file_io
Edit after answers' input:
- So capitalized letter in Fortran aren't the norm anymore? I didn't know that!
- I do plan to parallelize the code with either/both OpenMP, MPI
- @KyleKanos I was not sure of using local variables because that requires allocating/deallocating large arrays upon some subroutines' call (e.g.
M_plus
andM_minus
inget_Mhalf
), which seems quite inefficient to me. The alternative would be to have these variables global to the module only, but @haraldkl advises against it.