Using ChatGPT to Understand what I programmed in 1988 for SWMM4

Introduction

Background for this: I recently had to revisit the numerical solution for SWMM4 Runoff as it is part of our software programs XPSWMM and ICM InfoWorks SWMM Hydrology. Looking at the code without any comments made it a bit hard to understand. I used GPT4 to analyze the Fortran code. It found it to be a form of Richardson's equation and explained nicely how the subroutine works.

What does the Fortran Code Look like?

It is subroutine overlnd.for written in 1988. The structure is simple, but there are no comments other than my name. I guess I thought that was important at the time. The whole purpose of the routine is to solve WCC*DO**1.6666667, or the Runoff equation for SWMM4.

      SUBROUTINE OVERLND(KK,WC1,R1,WD,DEL,RDELT)
C     RUNOFF BLOCK
C     CALLED BY WSHED NEAR LINE 397
C=======================================================================
C      NONLINEAR RESERVOIR ROUTING FOR SURFACE RUNOFF
C                CREATED MARCH, 1988 BY BOB DICKINSON
C                UPDATED OCTOBER, 1988 BY BOB DICKINSON
C=======================================================================
      INCLUDE 'TAPES.INC'
      INCLUDE 'NEW88.INC'
      DIMENSION NSEQ(4),YZ(4),DY(4),DZ(4)
      DATA NSEQ/1,2,4,8/,NUSE/4/
C=======================================================================
      KLND(KK) = KLND(KK) + 1
      DO 100 J = 1,NUSE
      DOI      = WD
      DELT     = RDELT/FLOAT(NSEQ(J))
      RSTAR    = R1*DELT
      WCC      = WC1*DELT
      DO 400 K = 1,NSEQ(J)
      DELR     = 1.0E-10
      DO 300 I = 1,21
      DO       =  DOI + DELR/2.0
      IF(DO.LE.0.0) THEN
                    TEST = -RSTAR/WCC
                    IF(TEST.GT.0.0) THEN
                                    DO  = TEST**0.666667
                                    DEL = (DO - DOI)/2.0
                                    ELSE
                                    DEL = -DOI
                                    ENDIF
                    GO TO 400
                    ENDIF
      F   = DELR  - RSTAR - WCC*DO**1.6666667
      DF  =  1.0  - 0.8333333*WCC*DO**0.6666667
      DEL = DELR  - F/DF
      IF(I.NE.1.AND.ABS(DEL-DELR).LE.ABS(0.001*DELR)) GO TO 400
  300 DELR  = DEL
  400 DOI   = DOI + DEL
      YZ(J) = DOI - WD
      IF(J.EQ.1) DY(1) = YZ(1)
      IF(J.EQ.2) THEN
                 DY(2) = (4.0*YZ(2) - YZ(1))/3.0
                 IF(ABS(DY(2)-DY(1)).LE.0.001*ABS(DY(1))) THEN
                                                    DEL     = DY(2)
                                                    LND(KK) = LND(KK)+3
                                                    RETURN
                                                    ENDIF
                 ENDIF
      IF(J.EQ.3) THEN
                 DY(3) = (4.0*YZ(3)  - YZ(2))/3.0
                 DZ(3) = (16.0*DY(3) - DY(2))/15.0
                 IF(ABS(DZ(3)-DY(3)).LE.0.001*ABS(DZ(3))) THEN
                                                    DEL     = DZ(3)
                                                    LND(KK) = LND(KK)+7
                                                    RETURN
                                                    ENDIF
                 ENDIF
      IF(J.EQ.4) THEN
                 DY(4)   = (4.0*YZ(4)  - YZ(3))/3.0
                 DZ(4)   = (16.0*DY(4) - DY(3))/15.0
                 DZ(4)   = (64.0*DZ(4) - DZ(3))/63.0
                 DEL     = DZ(4)
                 LND(KK) = LND(KK) + 15
                 ENDIF
 100  CONTINUE
      RETURN
      END        

How does ChatGPT explain the SWMM4 Runoff code?

It is a Richardson extrapolation. While the code doesn’t correspond to a standard, widely-known numerical method like Runge-Kutta or Euler's method, it seems to be a customized implementation that combines elements of Richardson extrapolation for error control and an iterative approach for solving the SWMM4 nonlinear runoff equation.

The code snippet you provided appears to implement a numerical method for solving differential equations. Based on the structure and approach of the code, it seems to be utilizing a form of Richardson Extrapolation combined with an iterative solution approach, likely for a type of root-finding or differential equation problem.

Here’s a breakdown of what the method seems to be doing:
Richardson Extrapolation: This is a technique used to increase the accuracy of an approximation by combining approximations with different step sizes. In your code, the NSEQ array and its use in scaling DELT suggest that calculations are being performed at different step sizes (1, 2, 4, 8), and the results are combined to improve accuracy. This is a hallmark of Richardson Extrapolation.

Iterative Solution: The code shows a loop (DO 300) indicating an iterative method. This loop seems to be refining the value of DELR, which suggests an iterative process to find a solution (likely a root of an equation, given the structure of the code).

Convergence Check: The code checks for convergence (ABS(DEL-DELR).LE.ABS(0.001*DELR)) to determine if the solution has stabilized within a certain tolerance. This is common in numerical methods where iterative refinement is used.

Sequential Refinement: The outer loop (DO 100) and the checks (IF(J.EQ.2), etc.) indicate that the solution is being refined sequentially, using the results of smaller step sizes to inform and improve the accuracy of larger step sizes.
Error Estimation and Adaptation: The method appears to estimate the error (DY(2), DZ(3), DZ(4)) and adapts the calculation based on this estimation. This is typical in methods that aim to control the error of a numerical solution.

In conclusion, while the code doesn’t correspond to a standard, widely-known numerical method like Runge-Kutta or Euler's method, it seems to be a customized implementation that combines elements of Richardson Extrapolation for error control and an iterative approach for solving a specific type of problem, possibly related to root finding or differential equations. The specific context of the problem it’s designed to solve would provide more insights into the nuances of the method.        

The given code snippet appears to be implementing a Richardson extrapolation process within a numerical computation context. This is a Fortran program, which is often used for scientific computing due to its efficiency with array operations and numerical tasks. Let's break down how Richardson extrapolation is being applied here, step by step.

Solution By Steps
Step 1: Initializing Variables and Arrays
NSEQ array represents a sequence of step sizes, typically halved at each step (1, 2, 4, 8). This is crucial for Richardson extrapolation, which requires approximations at different step sizes.
YZ, DY, and DZ arrays are used to store intermediate and final results of the computations.
NUSE is the number of extrapolations to use.

Step 2: Iterative Computation
The outer loop (DO 100 J = 1, NUSE) iterates over the different step sizes.
For each step size, the code calculates a new value (DOI = DOI + DEL) based on some function, which seems to involve a root-solving method (like Newton-Raphson, judging from the update DEL = DELR - F/DF).

Step 3: Richardson Extrapolation Application
The Richardson extrapolation is applied in the nested IF statements within the loop.

For J=2, it performs a first-order extrapolation: DY(2) = (4.0*YZ(2) - YZ(1))/3.0. This formula is an application of the Richardson extrapolation where a higher order of accuracy is achieved by combining results from two different step sizes.

For J=3 and J=4, the code further refines the approximation, increasing the order of accuracy each time. This is done by combining results from previous steps and the current step, using coefficients derived from Richardson extrapolation theory.

Step 4: Convergence Check and Update
After each extrapolation, there is a check to see if the current approximation is sufficiently close to the previous one (ABS(DZ(3)-DY(3)).LE.0.001*ABS(DZ(3))). If it is, the loop exits, and the current approximation is considered accurate enough.

The LND(KK) variable is updated based on the level of accuracy achieved.
Final AnswerThis Fortran code applies Richardson extrapolation in a numerical computation, using a sequence of decreasing step sizes to iteratively refine an approximation. It starts with a basic approximation and improves its accuracy through successive applications of extrapolation formulas, each eliminating more of the error terms, until a convergence criterion is met.

Key ConceptRichardson extrapolation in numerical methods.
Key Concept ExplanationRichardson extrapolation is a technique for increasing the accuracy of approximations in numerical analysis by combining estimates from different step sizes. It leverages the idea that the error in numerical approximations can often be expressed as a power series of the step size. By calculating estimates with progressively smaller step sizes and combining these estimates in a specific way, Richardson extrapolation systematically eliminates higher-order error terms, leading to a more accurate result.        

How is this solved in the SWMM5 C Code?

It is solved as a 5th-order RK. Interestingly enough, the code is from Numerical Recipes in C. I took Richardson's Extrapolation code from the book Numerical Recipes in Fortran or a different edition of the book

//------------------------------------------------------------------------
//   odesolve.c
//
//   Fifth-order Runge-Kutta integration with adaptive step size control
//   based on code from Numerical Recipes in C (Cambridge University
//   Press, 1992).
//
//   Date:     11/15/06
//   Author:   L. Rossman
//----------        

Closing Note: Thank you so much for journeying with me through this content. This space is reserved for future updates and insights. Your engagement and time are truly appreciated. Until next time! You can also see my past articles on LinkedIn (91 in 2023).