8
\$\begingroup\$

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.
\$\endgroup\$
5
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ 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. \$\endgroup\$ – solalito Mar 6 '16 at 8:47
  • \$\begingroup\$ I thought you mentioned a hydrocode you wrote the first time I looked at your answer. Did you remove the github link? \$\endgroup\$ – solalito Mar 6 '16 at 12:10
  • \$\begingroup\$ 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) \$\endgroup\$ – Kyle Kanos Mar 6 '16 at 12:13
  • \$\begingroup\$ No need, already found it :) \$\endgroup\$ – solalito Mar 6 '16 at 12:19
3
\$\begingroup\$

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.
\$\endgroup\$
  • \$\begingroup\$ Thanks for your answer. Why is the time+dt summation a bad idea if I am doing too many steps? \$\endgroup\$ – solalito Mar 6 '16 at 8:45
  • \$\begingroup\$ @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. \$\endgroup\$ – haraldkl Mar 6 '16 at 9:26
2
\$\begingroup\$

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
    
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.