# 1D shock tube problem written in Fortran

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 and M_minus in get_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. ## 3 Answers Should I use Fortran structures? Best guess: probably yes. If you plan on making a career (academic or private), you will need to learn to design before building. For instance, do you want to include something like adaptive mesh refinement (AMR)? Or use (massive) parallelization? If either of those are true, it might be worth designing the code now to handle a Fortran type so that you can more easily pass information between processors or AMR levels. A code I used for my research for my doctorate stored information such as type data_struct integer, dimension(3,2) :: local_bounds, global_bounds integer, dimension(:), allocatable :: neigh_list real(dp) :: dt, dx(3) type(boundary), dimension(:), pointer :: boundaries real(dp), dimension(:,:,:,:), pointer :: q end type type(data_struct), pointer :: data  where local_bounds is the size (sans ghost cells) of q and global_bounds its location within the whole domain (see it's use for parallel sims?). The data type boundary would be a rank-2 array that holds the x-y boundary for z, the y-z for x, etc. There were a few other things added in there as well (e.g., qOld that held the previous time-step's solution in case the current step got screwed up somehow), but what's here shows the usefulness of a designing the type beforehand. Speaking of which, you don't actually do anything with boundaries in your code, despite apply_bc being there. Probably okay for the shock tube since you'll stop the simulation long before the shock reaches the boundary, but it is a necessity for a generalized hydrodynamics code. If you don't plan on doing (massive) parallelization or AMR, then it probably would be fine to just do as you've done & have some global variables. Functions vs. subroutines: I've always gone by the adage: a function returns a single variable; subroutines modify multiple variables.1 Obviously this can be abused because Fortran uses pass by reference, but I think that it should be reduced to as minimal as possible. If you are worried about speed, you may want to rethink those where constructs. Sure it's cleaner code with where, but do loops tend to be faster. It used to be that where was intended to parallelize code using SIMD instructions, but since then, compilers have been able to do the same vectorization with do loops. I don't know that anyone really uses where in more professional codes (it's not in the hydrocode I used), but it could be used. Where should I store local variables? I'm not sure why this is a question. If you have a variable that's only used in one procedure, it should be declared only in that procedure. If you have a variable that's common to a few procedures (e.g., half=0.5_dp), then you may want to put it as a global variable. I was not sure of using local variables because that requires allocating/deallocating large arrays upon some subroutines' call (e.g. M_plus and M_minus in get_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. Actually, another alternative would be to use inline functions (as I mention below) and use single variables, rather than arrays. Another alternative would be to have all the procedures under a contains block of the AUSMscheme procedure: subroutine AUSMscheme ... declare all variables here ... contains subroutine get_mach ... end subroutine get_mach end subroutine AUSMscheme  In this case, the variables declared in AUSMscheme are global to the containsed procedures, so you wouldn't have to be constantly allocating & deallocating the arrays used in the procedures. Other recommendations • I'm not much a fan of my code "yelling" at me (i.e., I don't like capitalized keywords). • I'm not sure why you're allocateing arrays when Nx is a parameter, just give it the dimensions at initialization. • On second thought, this probably would come in handy if Nx is a run-time value, passed in through a namelist, for example. • You'll probably want to look at other file format types if you want to go to 3D (e.g., XML, VTK, HDF5, netCDF, etc). The "gnuplot" style x,y1,y2,y3,... works for 1D, but for going to 2/3 D, you'll want something more robust & visualization friendly (paraview, visit, yt, etc). • If you want to use different integration methods (as evidenced by WENOmethod.f90 in the Makefile), you may want to investigate procedure pointers. • You probably should output some diagnostics at each time-step (current time, current dt, average energy, etc); it might seem wasteful to be calling print or write, but it's actually pretty useful for figuring out where your simulation died (if it died). • You should probably add some extra global variables like gamma1=mGamma-1, gamma2=1/gamma1; the compiler may be smart enough to decide that it can turn the division by mGamma-1 into a multiplication of the value 1/(mGamma-1) but expressing it this way ensures that it will. • If you do opt to write out do loops instead of the where constructs, some of the subroutines you've written can be turned into inline functions (e.g., get_sound, get_Mach) which can lead to reduced run times. 1. Note that I say variable here, rather than "value;" I think it's okay for functions to return a single array. • Thanks for your asnwer, please see edit. I was wondering about structure and procedure pointers and you've convinced me. Also, I thought the point of the WHERE construct was that it would be more efficient than the do loop. – solalito Mar 6 '16 at 8:47 • I thought you mentioned a hydrocode you wrote the first time I looked at your answer. Did you remove the github link? – solalito Mar 6 '16 at 12:10 • I had linked one in a comment but removed it because I later felt that it wasn't that good & didn't show what I thought it would. It's called spaceMARINE and is on github (on mobile so I can't link right now) – Kyle Kanos Mar 6 '16 at 12:13 • No need, already found it :) – solalito Mar 6 '16 at 12:19 Just a few additions to Kyle Kanos extensive answer. I like that you are using only for the uses. In my experience this helps the reader quite a bit to track where stuff is coming from. In addition you might also want to use the private statement and turn everything you want the module to export public explicitly. But be warned, that we saw trouble with the Intel compiler with long use list with only and private, where we sometimes needed to reorder the statements to avoid internal compiler errors. I think, more recent versions did never yield this problem anymore, so maybe it is not an issue anymore. One note on subroutines vs. functions: If you plan to incorporate OpenMP it may be easier to use subroutines instead of functions. In the type declaration I prefer the explicit kind= keyword, so instead of REAL(dp) I'd use real(kind=dp). Regarding comments: You might want to use some automatic code documentation tool like FORD, and format your comments accordingly. Regarding module variables: You are wise avoiding them, any global data might cause problems with optimization, especially if you go for share memory parallelization. Note, that your summation with time+dt may be dangerous if you are doing many timesteps. Also, the while loop is deprecated. All in all, it looks like pretty well written Fortran code, except for the all capital notation, which irritates me, just like Kyle Kanos. Update to question updates: • I'd say capitalization is not a norm anymore, to me it seems that only people not used to Fortran, or those who somehow think the language stopped at F77 are using it. • Of course, allocation and deallocation is a thing to avoid. However, global variables are not the only alternative, you could pass working arrays also through the procedure interfaces. Obviously this has the drawback to clutter the interfaces. A solution to this could be to use derived data types, where you allocate the necessary working arrays for the kernel once and then pass around the kernel datastructure, without the need to care about its internals anymore. It would then be easier to use dedicated working arrays for each thread for example, when you go for parallelization. Another alternative, that Kyle Kanos already pointed to would be to use automatic arrays, which would put the data on the stack and should be much faster than allocation/deallocation, but you need to ensure a sufficiently large stack space. • Thanks for your answer. Why is the time+dt summation a bad idea if I am doing too many steps? – solalito Mar 6 '16 at 8:45 • @solalito Eventually you will reach the limit of accuracy and time+dt == time, the addition will not make a difference anymore. I remember a question on SO, where this was the problem. It's unlikely that this will be an issue with your application scenario right here, but better to be aware of that. – haraldkl Mar 6 '16 at 9:26 A couple of quick notes on your makefile: • Add a standard to your FLAGS (e.g., -std=f95, -std=f2003, etc). Just adding this flag to the build scripts of existing projects frequently identifier possible areas of non portability. • Explicitly declare that the .mod files are a target, by expanding your genetic object code building rule to: %.o %.mod: %.f90$(FC) $(FLAGS) -c$< -o $(subst .f90,.o,$<)

• Explicitly state inter-module dependencies (note that the dependency on the source code file - file_io.f90 in the example below - is captured automatically by the generic building rule):

file_io.o: global_vars.mod