Fix numerical issues in int/runge_kutta.c

int/runge_kutta.c
- Initialize Hacc, ErrOld, NewtonIncrementOld
- Add error checks to prevent div-by-zero errors
- Now set *IERR = Code in Rk_ErrorMsg (was Code = *IERR)
- Fixed missing decimal point for rkb[2] Radau2A coefficient
- Added missing rkT[0][0] and rkt[0][1] coefficients

CHANGELOG.md
- Updated accordingly

Signed-off-by: Bob Yantosca <[email protected]>

Author: yantosca

CHANGELOG.md

@@ -26,6 +26,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Updated rules to ignore files in `.gitignore` and updated comments accordingly
 - Fixed a bug that prevented `.ci-pipelines/ci-cleanup-script.sh` from removing KPP-generated files for MCM mechanisms
 - Fixed typo in error message in `int/rosenbrock_autoreduce.f90`
+- Fixed several issues in `int/runge_kutta.c` that were causing numerical instability or non-deterministic results
 
 ### Removed
 - Removed C-I tests on Microsoft Azure Dev Pipelines

int/runge_kutta.c

@@ -670,6 +670,9 @@ void RK_Integrator( int N,
        SkipJac,
        NewtonDone,
        SkipLU;
+   /* Initialize variables to prevent use of uninitialized values */
+   Hacc = ZERO;
+   ErrOld = ZERO;
 
 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 
@@ -769,6 +772,7 @@ Tloop:   while ( (Tend-*T)*Tdirection - Roundoff > ZERO ) {
       /*~~~> Initializations for Newton iteration */
       NewtonDone = 0;
       Fac = (KPP_REAL)0.5; /* Step reduction if too many iterations */
+      NewtonIncrementOld = Roundoff; /* Initialize to prevent division by zero */
 
       /* for NewtonLoop */
       for (NewtonIter = 1; NewtonIter <= NewtonMaxit; NewtonIter++)
@@ -791,8 +795,12 @@ Tloop:   while ( (Tend-*T)*Tdirection - Roundoff > ZERO ) {
 	   Theta = ABS(ThetaMin);
 	   NewtonRate = (KPP_REAL)2.0;
 	}
-	else {
+	   else {
 	   /*printf("Entering else from NewtonIter=%d == 1\n", NewtonIter);*/
+	   /* Safety check: prevent division by zero */
+	   if (NewtonIncrementOld <= Roundoff) {
+		NewtonIncrementOld = MAX(NewtonIncrement, Roundoff);
+	   }
 	   Theta = NewtonIncrement / NewtonIncrementOld;
 	   /*printf("Theta=%f\n", Theta);*/
 	   if (Theta < (KPP_REAL)0.99) {
@@ -874,6 +882,7 @@ Tloop:   while ( (Tend-*T)*Tdirection - Roundoff > ZERO ) {
       	/*~~~> Initializations for Newton iteration */
       	NewtonDone = 0;
       	Fac = (KPP_REAL)0.5; /* Step reduction factor if too many iterations */
+      	NewtonIncrementOld = Roundoff; /* Initialize to prevent division by zero */
 
       	/* for SDNewtonLoop */
       	for (NewtonIter = 1; NewtonIter <= NewtonMaxit; NewtonIter++)
@@ -884,6 +893,11 @@ Tloop:   while ( (Tend-*T)*Tdirection - Roundoff > ZERO ) {
 	   FUN_CHEM(*T+H,TMP,DZ4);	/* DZ4 <- Fun(Y+Z4) */
 	   ISTATUS[Nfun]++;
 	   /* DZ4[1,N] = (D[1, N]-Z4[1,N])*(rkGamma/H) + DZ4[1,N]; */
+	   /* Safety check: H should not be zero here, but add protection */
+	   if (ABS(H) <= Roundoff) {
+		*IERR = -10;
+		return;
+	   }
 	   WAXPY(N,-ONE*rkGamma/H,Z4,1,DZ4,1);
 	   WAXPY(N,rkGamma/H,G,1,DZ4,1);
 
@@ -902,6 +916,10 @@ Tloop:   while ( (Tend-*T)*Tdirection - Roundoff > ZERO ) {
 	   }
 	   else {
 	      /*printf("Entering else from NewtonIter=%d == 1\n", NewtonIter);*/
+	      /* Safety check: prevent division by zero */
+	      if (NewtonIncrementOld <= Roundoff) {
+		NewtonIncrementOld = MAX(NewtonIncrement, Roundoff);
+	      }
 	      ThetaSD = NewtonIncrement / NewtonIncrementOld;
 	      if (ThetaSD < (KPP_REAL)0.99) {
 		/*printf( "Entering if Theta=%g < 0.99\n", Theta );*/
@@ -1001,7 +1019,13 @@ Tloop:   while ( (Tend-*T)*Tdirection - Roundoff > ZERO ) {
       	if (Gustafsson == 1) {
 	   /*~~~> Predictive controller of Gustafsson */
 	   if (ISTATUS[Nacc] > 1) {
-	      FacGus = FacSafe*(H/Hacc)*pow(Err*Err/ErrOld,(KPP_REAL)(-0.25));
+	      /* Safety check: prevent division by zero or near-zero values */
+	      if (ABS(Hacc) > ((KPP_REAL)10.0)*WLAMCH('E') && 
+		  ABS(ErrOld) > ((KPP_REAL)10.0)*WLAMCH('E')) {
+		  FacGus = FacSafe*(H/Hacc)*pow(Err*Err/ErrOld,(KPP_REAL)(-0.25));
+	      } else {
+		  FacGus = Fac;
+	      }
 	      FacGus = MIN(FacMax,MAX(FacMin,FacGus));
 	      Fac = MIN(Fac,FacGus);	
 	      Hnew = Fac*H;
@@ -1067,7 +1091,7 @@ Tloop:   while ( (Tend-*T)*Tdirection - Roundoff > ZERO ) {
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 void RK_ErrorMsg(int Code, KPP_REAL T, KPP_REAL H, int* IERR)
 {
-   Code = *IERR;
+   *IERR = Code;
    printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
    printf("\nForced exit from RungeKutta due to the following error:\n");
 
@@ -1177,6 +1201,10 @@ void RK_Interpolate(char action[], int N, KPP_REAL H, KPP_REAL Hold, KPP_REAL Z1
    /* Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3 */
    if ( strncmp(action, "make", 4 ) ) {
 	den = (rkC[2]-rkC[1])*(rkC[1]-rkC[0])*(rkC[0]-rkC[2]);
+	/* Safety check: prevent division by zero or near-zero denominator */
+	if (ABS(den) <= ((KPP_REAL)10.0)*WLAMCH('E')) {
+		return; /* Denominator is too small, cannot interpolate */
+	}
 	for (i = 0; i < N; i++)
 	{
 	   CONT[i][0] = (-rkC[2]*rkC[2]*rkC[1]*Z1[i]
@@ -1195,6 +1223,10 @@ void RK_Interpolate(char action[], int N, KPP_REAL H, KPP_REAL Hold, KPP_REAL Z1
    }
    /* Evaluate quadratic polynomial */
    else if ( strncmp( action, "eval", 4 ) ) {
+        /* Safety check: prevent division by zero or near-zero Hold */
+        if (ABS(Hold) <= ((KPP_REAL)10.0)*WLAMCH('E')) {
+		return; /* Hold is too small, cannot interpolate */
+	}
         r  = H / Hold;
 	x1 = ONE + rkC[0]*r;
 	x2 = ONE + rkC[1]*r;
@@ -1258,6 +1290,11 @@ void RK_Decomp(int N, KPP_REAL H, KPP_REAL FJAC[], KPP_REAL E1[],
    KPP_REAL Alpha, Beta, Gamma;
    int i,j;
 
+   /* Safety check: prevent division by zero or near-zero H */
+   if (ABS(H) <= ((KPP_REAL)10.0)*WLAMCH('E')) {
+	*ISING = 1;
+	return;
+   }
    Gamma = rkGamma / H;
    Alpha = rkAlpha / H;
    Beta  = rkBeta / H;
@@ -1318,6 +1355,10 @@ void RK_Solve(int N, KPP_REAL H, KPP_REAL E1[], int IP1[], KPP_REAL E2R[],
    int i;
 
    /* Z <- h^{-1) T^{-1) A^{-1) x Z */
+   /* Safety check: prevent division by zero or near-zero H */
+   if (ABS(H) <= ((KPP_REAL)10.0)*WLAMCH('E')) {
+	return; /* H is too small, skip this operation */
+   }
    for (i = 0; i < N; i++)
    {
 	x1 = R1[i]/H;
@@ -1368,6 +1409,11 @@ void RK_ErrorEstimate(int N, KPP_REAL H, KPP_REAL T, KPP_REAL Y[],
    KPP_REAL HrkE1,HrkE2,HrkE3;
    int i;
 
+   /* Safety check: prevent division by zero or near-zero H */
+   if (ABS(H) <= ((KPP_REAL)10.0)*WLAMCH('E')) {
+	*Err = (KPP_REAL)1.0e10; /* Large error if H is too small */
+	return;
+   }
    HrkE1  = rkE[0]/H;
    HrkE2  = rkE[1]/H;
    HrkE3  = rkE[2]/H;
@@ -1458,7 +1504,7 @@ void Radau2A_Coefficients()
 
    rkB[0] = (KPP_REAL)(3.764030627004672750500754423692808e-01);
    rkB[1] = (KPP_REAL)(5.124858261884216138388134465196080e-01);
-   rkB[2] = (KPP_REAL)(1111111111111111111111111111111111e-01);
+   rkB[2] = (KPP_REAL)(1.111111111111111111111111111111111e-01);
 
    rkC[0] = (KPP_REAL)(1.550510257216821901802715925294109e-01);
    rkC[1] = (KPP_REAL)(6.449489742783178098197284074705891e-01);
@@ -1636,8 +1682,8 @@ void Lobatto3C_Coefficients ()
    rkAlpha = (KPP_REAL)1.687091590520766641994055533117359;
    rkBeta  = (KPP_REAL)2.508731754924880510838743672432351;
 
-   rkT[0][1] = ONE;
-   rkT[0][1] = ONE;
+   rkT[0][0] = (KPP_REAL).4554100411010284672111720348287483;
+   rkT[0][1] = (KPP_REAL)(-.6027050205505142336055860174143743);
    rkT[0][2] = ZERO;
    rkT[1][0] = (KPP_REAL).4554100411010284672111720348287483;
    rkT[1][1] = (KPP_REAL)(-.6027050205505142336055860174143743);