Detailed listing of the differences between the DHM and EDHM source codes.

## Abstract

The Diffusion Hydrodynamic Model or DHM is a coupled one- and two-dimensional (2-D) surface flow model based upon a diffusion formulation of the well-known Navier–Stokes equations, developed by research hydrologists of the USGS (United States Geological Survey) for use in modeling floodplains and dam-break situations. The Fortran 77 source code and various applications were published in 1987 by the USGS as a Technical Report authored by Hromadka and Yen. The DHM program led to the development of several subsequent computational programs such as the FLO-2D computational model and other similar programs. The original DHM program had a limit of applications to problems with no more than 250 nodes and modeling grids. That limitation was recently removed by a program version named EDHM (Extended DHM), which provides for 9999 nodes and grids. However, the computational code is preserved in order that the baseline code algorithmic procedures are untouched. In this paper, the DHM and EDHM are rigorously compared and examined to identify any variations between the two Fortran codes. It is concluded from this investigation that the two sets of algorithm codes are identical, and outcomes from either program are similar for appropriately sized applications.

### Keywords

- legacy Fortran codes
- computational economy
- large scale application
- overland flow
- flow through a constriction

## 1. Introduction

Legacy Fortran 77 codes that have been developed in the 1980s continue to have a wide spread audience (for both research and commercial applications) across all the Computational Fluid Dynamics disciplines. Their popularity can be largely attributed to the extensive validation that these models have been subjected to with analytical and experimental data for single and multi-dimensional flows. The trust that the audience have in the end results from these legacy codes, and their ability to meet the user goals are other driving factors for making them popular among the modeling community. Few of the legacy Fortran 77 codes, developed by various groups include MFIX [1] (Open source software for simulating Multiphase Flow with Interphase eXchanges), VOF 2D [2] (Two- dimensional, transient, free-surface incompressible fluid dynamics program), FUN3D [3] (fully unstructured-grid fluid dynamic simulations spanning incompressible flow to transonic flow), LAURA [4] (Langley Aerothermodynamic Upwind Relaxation Algorithm for structured, multi-block, computational aerothermodynamic simulations) and INS3D [5] (incompressible Navier–Stokes equations in three dimensions for steady and transient flows).

Although theoretically, these legacy codes are on a firm footing, computationally, they are uneconomical. When applied over a large scale application, characterized by thousands of nodes, these codes are constrained by varying degrees. The limitations arise partly from (a) the lack of object-oriented tools and the absence of abstract modeling capabilities in Fortran 77 and (b) the required CPU time when these models are applied across large domains. Balancing the accuracy of simulation with acceptable CPU is a crucial element that the current modelers are looking for. Capturing the physics of some flow phenomena necessitates that the equations be applied across the small spatial grid and temporal scales, which can be an issue in applying legacy codes. Since the codes were written (at that time) for the then available computational resources, applying them, as they are, for large scale domains might not be feasible either due to the large CPU time that they need to complete the simulation or because of the array limitations or modifications that need to be made so that the codes can be compiled using the currently available Fortran 77 compilers. Addressing these limitations will result in a ‘modernized’ version of the legacy codes. Modernization does not mean a better numerical formulation or a more accurate code. It only refers to a computationally efficient code with perhaps a better user interface for input and for visualizing the results through colorful multi-dimensional graphs and tables. In fact, using any modernized code without extensive benchmark testing can lead to erroneous solutions.

The above first limitation was addressed by researchers either by rewriting the entire code from scratch using an object-oriented programming language or by using an incremental approach in which the computationally intensive modules in legacy codes are identified and replaced with their computationally efficient counterparts. For instance, in an implicit finite difference or finite element formulation, where the system of equations are assembled in the form Ax = B, much of the computational time is spent on solving for ‘x’ [6]. While for the application audience, a solution module merely represents a means to an end of solving the flow equations, for a solver developer, the application is a source of sparse equations to be solved. Using an appropriate solver and a preconditioner can significantly cut down the simulation time, thus partly addressing the limitation.

The advent of high-performance computing tools and their availability to all the audience (based on commodity processors) saw the evolution of legacy Fortran codes to their full or semi-parallel versions, thus addressing the second limitation in legacy codes. There are many advantages to using parallel codes. For an application modeler using a well-written parallel code, the primary advantage is the reduction in the computational time. This reduction is on the order of the number of processors used in the communicator. Parallelization allows a program to be executed by as many processors available within the sub-complex simultaneously, thus facilitating a steep reduction in the computational time required for the solution to converge to the desired simulation time. This reduction enables the code to be applied over domains with millions of nodes or cells, to better capture the physics of the flow. Using standard parallel libraries that use MPI protocols like PETSc [7], ScaLAPACK [8], and BLAS [9], existing serial codes can be converted to parallel codes with reduced effort. Alternatively, highly optimized parallel versions of legacy codes can be written from scratch, which involves higher costs for developing and testing the code.

Diffusion Hydrodynamic Model (DHM) is a legacy model developed in the 1980s for USGS by the first author and his colleague [10]. The report and the Fortran 77 source code are also available at the DHM companion web page

The layout of this document is as follows. In Section 2, the flow equations and other salient characteristics in DHM are briefly described. Section 3 lists the computational limitations in DHM. In Section 4, the modifications done in DHM to arrive at EDHM are detailed. Performance tests to validate the reliability of EDHM are discussed in Section 5. Conclusions are presented in Section 6.

## 2. Overview of DHM

The two-dimensional flow continuity and momentum equations along the X and Y axis (assuming a constant fluid density without sources or sinks in the flow field and hydrostatic pressure distribution) can be written as [10].

in which

The local and convective acceleration terms can be grouped and Eqs. 2 and 3 are rewritten as

where

Eq. 5 can be rewritten in the general case as

where

The symbol S in Eq. 7 indicates the flow direction which makes an angle of

Two-dimensional DHM is formulated by substituting Eq. 8 into Eq. 1

If the momentum term groupings were retained, Eq. 9 can be written as

where

and

To maintain continuity in the discussion, while salient aspects of the DHM numerical algorithm are presented here, readers are referred to [10] for a detailed description of the numerical formulation and the input file format. The domain is divided into uniform square grids or cells. For each interior grid, its connectivity with adjacent grids along the North, East, South, and West directions is specified. For grids that are on the boundaries, ‘0’ is specified along the directions that do not have adjacent nodes. The flow equation is solved using the nodal domain integration method. Apart from the grid connectors, at the center of each cell, the required input variables are the roughness value, ground elevation, initial flow depth. The number of inflow cells and the inflow hydrograph at each of them need to be inputted. The number and the outflow boundary cell numbers should be identified. Since the formulation is explicit, the choice of time step (Δt) is limited by the Courant-Friedrich-Lewy stability condition. Starting from time = 0, the explicit solution is marched in the direction of time, until the required transient time level is reached. DHM gives the option of printing the output variables at any time level. The output at the center of each cell includes the flow depth, elevation, and flow velocities along the four directions.

## 3. Computational limitations in DHM

For modeling flows over large sizes domain, the two primary shortcomings in DHM are

The maximum number of cells (nodes) that can be accommodated in DHM is limited to 250

Inflow and outflow boundary nodes are limited to 10

Both these limitations were largely due to the computational resources that were available to the developers in the 1980’s. Application of DHM over large flow domains would require using a higher number of nodes in the computational domain, warranting modifications to DHM, as discussed next.

## 4. Features in extended DHM (EDHM)

The changes made in DHM can be grouped into three categories, as detailed below. No changes were made to the format or structure of variables input and output files.

### 4.1 Major enhancements

Increased the array size of variables (FP,FC) that stores the location of cells and the initial depth, average elevation, roughness coefficient values, velocity (VEL), maximum water depth, and the corresponding time at the cell (DMAX, TIMEX) from 250 to 9999.

Increased the array size of variables that store the stage curve data for the channel (NOSTA, STA) from 10 to 99.

Increased effective rainfall intensity data pairs from 10 to 99.

Increased the array size of variables relating to the inflow and outflow hydrograph nodes (KIN, KOUT), depth of the specified stage-discharge curve (HOUT), inflow boundary condition nodal points (KINP) and the inflow hydrograph details (HP) from 10 to 99

Increased of the array which stores the nodal points where outflow hydrographs are being printed (NODFX, NODCC) from 50 to 99

The changes in the above array sizes were done in the main code and the associated subroutines FLOODC, QFP, QFC, and CHANPL.

### 4.2 Minor enhancements

The two minor enhancements that were made in DHM code are (a) to accommodate the increased number of cells in EDHM, the fixed format output descriptor has been expanded by one digit and (b) to better align the variables in the output file, the inter variable spacing was decreased by one digit. A detailed listing of all the major and minor changes made in the DHM source code, along with the corresponding line numbers, is shown in Table 1.

DHM | EDHM | ||
---|---|---|---|

Line # | Content | Line # | Content |

12 | COMMON/BLK 1/FP(250,8),FC(250,6) | 12 | COMMON/BLK 1/FP(9999,8),FC(9999,6) |

14 | COMMON/BLK 2/KIN(10),H(10,15,2), KOUT(10),HOUT(10,15,3) | 14 | COMMON/BLK 2/KIN(99),H(99,15,2), KOUT(99),HOUT(99,15,3) |

16 | COMMON/BLK 3/NOSTA(10),STA(10,15,2),NODFX(50) | 16 | COMMON/BLK 3/NOSTA(99),STA(99,15,2),NODFX(99) |

18 | COMMON/BLK 4/DMAX(250,2),TIMEX(250,2) | 18 | COMMON/BLK 4/DMAX(9999,2),TIMEX(9999,2) |

20 | COMMON/BLK 5/KINP(10),HP(10,15,2) | 20 | COMMON/BLK 5/KINP(99),HP(99,15,2) |

26 | DIMENSION NODDC(50),VEL(250,4),R(10,2),Q(4) | 26 | DIMENSION NODDC(99),VEL(9999,4),R(99,2),Q(4) |

271 | FORMAT(10X,5I4,1X,F6.4,2X,F6.1,1X,F5.1) | 271 | FORMAT(10X,5I5,1X,F6.4,2X,F6.1,1X,F5.1) |

281 | FORMAT(/,10X,'INFLOW HYDROGRAPH AT NODE #’,I3,/, | 281 | FORMAT(/,10X,'INFLOW HYDROGRAPH AT NODE #’,I4,/, |

291 | FORMAT(10X,I3,1X,I3) | 291 | FORMAT(10X,I4,1X,I4) |

297 | FORMAT(10X,I3 ,2X,F5.4,1X,F7.1,1X,F7.1,1X,F7.1,5X,F7.1) | 297 | FORMAT(10X,I4 ,2X,F5.4,1X,F7.1,1X,F7.1,1X,F7.1,5X,F7.1) |

303 | FORMAT(10X,'OUTFLOW NODE # ‘,I3, | 303 | FORMAT(10X,'OUTFLOW NODE # ‘,I4, |

311 | FORMAT(/,10X,'STAGE CURVE AT NODE #’,I3,/, | 311 | FORMAT(/,10X,'STAGE CURVE AT NODE #’,I4,/, |

333 | FORMAT(10X,'INFLOW RATE AT NODE ‘,I3,’ IS EQUAL TO ‘,F10.2) | 333 | FORMAT(10X,'INFLOW RATE AT NODE ‘,I4,’ IS EQUAL TO ‘,F10.2) |

335 | FORMAT(/,5X,'NODE’,7X,10(I3,8X)) | 335 | FORMAT(/,5X,'NODE’,7X,10(I4,7X)) |

357 | FORMAT(10X,'OUTFLOW RATE AT NODE ‘,I3,’ IS EQUAL TO ‘,F10.2) | 357 | FORMAT(10X,'OUTFLOW RATE AT NODE ‘,I4,’ IS EQUAL TO ‘,F10.2) |

1275 | COMMON/BLK 1/FP(250,8),FC(250,6) | 1275 | COMMON/BLK 1/FP(9999,8),FC(9999,6) |

1277 | COMMON/BLK 2/KIN(10),H(10,15,2), KOUT(10),HOUT(10,15,3) | 1277 | COMMON/BLK 2/KIN(99),H(99,15,2), KOUT(99),HOUT(99,15,3) |

1279 | COMMON/BLK 3/NOSTA(10),STA(10,15,2),NODFX(50) | 1279 | COMMON/BLK 3/NOSTA(99),STA(99,15,2),NODFX(99) |

1281 | COMMON/BLK 4/DMAX(250,2),TIMEX(250,2) | 1281 | COMMON/BLK 4/DMAX(9999,2),TIMEX(9999,2) |

1567 | COMMON/BLK 1/FP(250,8),FC(250,6) | 1567 | COMMON/BLK 1/FP(9999,8),FC(9999,6) |

1647 | COMMON/BLK 1/FP(250,8),FC(250,6) | 1647 | COMMON/BLK 1/FP(9999,8),FC(9999,6) |

1727 | COMMON/BLK 1/FP(250,8),FC(250,6) | 1727 | COMMON/BLK 1/FP(9999,8),FC(9999,6) |

### 4.3 Compiler details

After reviewing the currently available compilers in Windows for Fortran 77 codes, we have chosen the Intel Fortran Compiler within the Microsoft Visual Studio integrated development environment (IDE) to make the enhancements in DHM and for generating the.EXE file. This interface is ideal to debug and execute Fortran 77 programs. The compiler can optimize the performance of source codes for Intel CPUs. It offers broad support for current and previous Fortran standards and also tools by which a robust, high-performance code can be created in serial and parallel environments. The Math Kernel Library (Intel MKL) and the Debugger tools in the compiler, creates a solid foundation for building robust, high-performance codes. The end executable file (.EXE) although optimized for Intel CPUs, can also run on x86 compatible CPUs such as those from AMD.

## 5. Application of EDHM

EDHM is qualitatively and quantitatively evaluated by comparing its results with the output from DHM, for two test simulations. The focus was to check if the EDHM solution resembles DHM output for these two cases for varying inflow and other model parameters.

*Case 1: Flow in a transition.*

Open channel flow through a linear contraction under the framework of two-dimensional flow is a common phenomenon and has drawn the attention of many experimental and numerical studies. The flow characteristics in the contraction depend on the Froude number at the upstream end. Flow in a contraction has acted as a benchmark simulation in investigations [12] that compared the performance of various CFD and hydraulic models. Figure 1 is the definition sketch of the test problem. The rectangular channel is 380 ft. long and 260 ft. wide. The constriction portion of the channel is 60 ft. x 60 ft. The channel length before and after constriction is 120 ft. and 200 ft., respectively. The cell size in the domain is 20 ft., and the total number of cells are 239. Figure 2 illustrates the cell numbers in the domain. The elevation of cells along the north and south boundaries was assigned a high value to physically denote that they are walls. The flow is confined within these boundaries. At the upstream end, cells 3-11 (nine cells) were specified with a constant inflow of 33.33 cfs. At the downstream end, a free outfall boundary was specified. The transient simulation was carried out until time = 1 hour.

Figure 3 plots the water depth profile for the two models along the channel centerline at time = 0.9 hours for a channel bottom Manning’s roughness value of 0.015. The close agreement of results gives confidence that the changes made to arrive at EDHM did not lead to different output data. The analysis over varying inflow discharge and bottom roughness values did not change the trend of the results. Figures 4 and 5 compare the depth profile for the roughness coefficient of 0.024 and 0.04. A similar agreement in output data (including velocity components) was observed across other longitudinal sections.

*Case 2: Overland flow.*

Overland flow over a hill slope generated by a rainfall event is characterized by varying hydraulic properties, roughness values, topography, and physical features in the domain. For predicting the hydraulic and hydrologic properties of flow, various models that solve a range of equations from a one-dimensional hydrodynamic equation for homogeneous place surfaces [13, 14] to 2D full non-linear shallow water Equations [15] have been applied. Figure 6 shows the flow domain. The number of cells in the domain is 56 and are 30 ft. in size. The roughness value ranged between 0.015 to 0.03Inflow hydrograph (Figure 7) was applied at cells 1 and 15. The transient model was run until time = 2 hours. The elevation drop along the north and south end of the domain is 26 ft., and the drop along the west and east boundaries is 5 ft. Cells 42 to 56 were specified as critical depth outflow nodes. Figures 8 and 9 compare the transient depth profiles from DHM and EDHM (time = 0 to 2 hours) at cells 36 and 26, respectively. The trend of the depth and velocity values at other cells also indicated the close agreement of flow data from the two models.

## 6. Conclusions

USGS Diffusion Hydrodynamic Model is a legacy Fortran 77 that has been widely applied for multiple one and two-dimensional flow scenarios. The computational limitations in the model which prevent its application over large domains have been addressed in this paper. The enhanced model can accommodate 9999 cells and 99 inflow and outflow nodes. Based on the analysis that was carried out and the close agreement of the results between the two models gives confidence in the reliability of the extended model. Current findings encourage future development of parallel EDHM in order to reduce the computational time.