Bugs in SWMS_2D and SWMS_3D

[mecej4]

I am posting here to report some bugs in the old software packages SWMS_2D and SWMS_3D, the source code of which is available at this site (www.pc-progress.com), or from https://github.com/usda-ars-ussl/legacy . These old packages, I realize, are no longer maintained and may see little use these days. However, a recent post by a user on the Intel Fortran Compiler user forum prompted me to look into the source codes of these packages, and I found that both contain a number of bugs. These bugs were not serious when the old Microsoft Powerstation Fortran compiler (FPS4) was used, and floating point calculations were performed using x87 FPU instructions. With modern compilers, however, when building 64-bit applications or using SSE2 and later FPU instructions, serious errors can occur.

I have not seen the source code of your more recent software, such as Hydrus_3D, but it is worthwhile to check if vestiges of the bugs in SWMS have been carried over into the new software.

With the Intel Fortran compilers of the era 2010-2020, there were compiler flags available to build 32-bit applications that behaved similarly to those compiled using FPS4. Those features are going away. Using current Ifort and Gfortran compilers, the user who posted to the Intel forum reported seeing lots of NaNs in the output of SWMS_3D. With SWMS_2D, I found, there were no NaNs but EXEs compiled with different compilers gave substantially different results on some example problems, and errors of this type can easily go unnoticed, whereas NaNs do not.

There are three categories of bugs.

1. Uninitialized variables. The FPS4 compiler gave all variables the SAVE attribute and initialized variables to zero. Few current compilers do that nowadays, and compiler flags such as /Qzero and /Qsave are like sledge-hammers. I'll be happy to provide a list of these errors -- there are just a handful, and the fixes are easy, and I will make them available if there is interest.

2. Insufficient floating-point precision. The x87 instructions generated by FPS4 often resulted in intermediate results having higher precision than one would expect from 4-byte REAL variables, because the intermediate results were held in the 80-bit FPU registers. With SSE2 and AVX instructions, this does not happen, and the resulting loss of precision can be catastrophic. This is the cause of the NaNs and unreliable results. There are similar differences relating to denormal numbers and gradual underflow to zero.

Here are some examples, followed by the fix for the problem:

Last lines of console output from SWMS_2D, Example.2, run using PGI 15.3 for Windows:

Without the fix: 273.0000 3 2794 -2.543E+01 4.432E+01 1.241E+00 -101. -190. 3.
With the fix: 273.0000 4 2733 -2.543E+01 4.434E+01 1.499E+00 -98. -180. 11.

Same, using Gfortran 10.2 on Cygwin-64, Windows 10:

Without the fix: 273.0000 4 2752 -2.543E+01 4.434E+01 1.566E+00 -98. -177. 13.
With the fix: 273.0000 4 2733 -2.543E+01 4.434E+01 1.499E+00 -98. -180. 11.

Fortunately, there is a simple fix.

In watflow2.f, subroutine RESET, add the declaration: double precision ae,qn,amul,bmul,fmul,xmul,cone,Bi,Ci
In watflow3.f, subroutine RESET, add the declaration: double precision det,ve,QN,caxx,cayy,cazz,caxy,caxz,cayz,
+ Amul,Bmul,Fmul,VolR,conE,BetaE,SinkE,hNewE

3. In subroutine WEFACT, the expression 1/tanh(x) - 1/x is evaluated in single precision. For small values of x, such as 4E-4, in single precision this is evaluated by Gfortran compiled code as zero, the true value being about 1E-4. Amazingly, in a few instances the evaluated value is 1 instead of 0 or the true value! A simple fix is to switch over to the approximation (1-x*x/15)*x/3 when x is smaller than, say, 0.1, or to evaluate the sensitive expression using double precision.

I'll be happy to provide any further details needed or clarifications, and look forward to your responses.

[Jirka]

Thanks for your comments on this legacy software. Indeed, we had to do a lot of the modifications that you list (e.g., initializations, double precision, etc) when we moved from the MS Fortran PowerStation to Digital Fortran, and later Intel Fortran.

Do you have a reference for this approximation: (1-x*x/15)*x/3 ?

I would be happy to look at further details if you could email them to me to Jiri.Simunek@ucr.edu.

Thank you. Jirka

[mecej4]

I appreciate the prompt response.

I have sent you email with the URL to a shared folder in the cloud that contains the corrected source files for SWMS_2D and SWMS_3D.

The polynomial approximation to the expression 1/tanh(x) - 1/x was obtained using the Taylor expansion at x = 0:

Using Maple mathematical software, for example,

R:\SWMS3D>d:\MATH\Maple17\bin.X86_64_WINDOWS\cmaple.exe
|\^/| Maple 17 (X86 64 WINDOWS)
._|\| |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2013
\ MAPLE / All rights reserved. Maple is a trademark of
<____ ____> Waterloo Maple Inc.
| Type ? for help.
> series(1/tanh(x)-1/x,x=0,7);
1/3 x - 1/45 x^3 + O(x^5 )

[Jirka]

Thanks. J.

[mecej4]

The legacy package CHAIN_2D shares uses some of the same subroutines as SWMS_2D, so it contains similar bugs.

The loss of precision in subroutine RESET and in the calculation of 1/tanh(x) - 1/x in subroutine WEFACT affects CHAIN_2D, as well. For Example.2, compare the last line of the console output obtained with different compilers.

As is, only A,B, B1 double precision in RESET:

Code: Select all

Gfortran 9.3       273.0000  6  1051  -2.543E+01  4.434E+01  1.560E+00  -102.  -200.    12.
Intel 9.1 X64      273.0000  6  1051  -2.543E+01  4.434E+01  1.559E+00  -102.  -200.    12.
Intel 19.1U4 X64   273.0000  6  1052  -2.543E+01  4.434E+01  1.560E+00  -102.  -200.    12.
PGI 15.3 X64       273.0000  6  1083  -2.543E+01  4.432E+01  1.235E+00  -106.  -219.     2.

With ae,qn,amul,bmul,fmul,xmul,cone,Bi,Ci,Bj,Cj,Bk,Ck promoted to double precision in RESET:

Code: Select all

Gfortran 9.3       273.0000  6  1067  -2.543E+01  4.434E+01  1.493E+00  -103.  -204.    10.
Intel 9.1 X64      273.0000  6  1067  -2.543E+01  4.434E+01  1.493E+00  -103.  -204.    10.
Intel 19.1U4 X64   273.0000  6  1067  -2.543E+01  4.434E+01  1.493E+00  -103.  -204.    10.
PGI 15.3 X64       273.0000  6  1067  -2.543E+01  4.434E+01  1.493E+00  -103.  -204.    10.

Note the significant differences in the 4th column from the right end.

I have added a Zip file with my modified sources of CHAIN_2D to the same shared folder in the cloud.