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 %%CODEBLOCK_0%%lt;               

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: w

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

StackExchange Code Review Q#121874, answer score: 6