Feature: Remove Implicit RHS Evaluation in Autonomous Problems With the Trivial Predictor

CHANGELOG.md

@@ -38,6 +38,12 @@ Added the following MRI coupling tables
 Added `ARKodeButcherTable_ERKIDToName` and `ARKodeButcherTable_DIRKIDToName` to
 convert a Butcher table ID to a string representation.
 
+Added the function ``ARKodeSetAutonomous`` in ARKODE to indicate that the
+implicit right-hand side function does not explicitly depend on time. When using
+the trivial predictor, an autonomous problem may reuse implicit function
+evaluations across stage solves to reduce the total number of function
+evaluations.
+
 Users may now disable interpolated output in ARKODE by passing `ARK_INTERP_NONE`
 to `ARKodeSetInterpolantType`. When interpolation is disabled, rootfinding is
 not supported, implicit methods must use the trivial predictor (the default

doc/arkode/guide/source/Usage/User_callable.rst

@@ -1769,6 +1769,7 @@ Optional input                                                  Function name
 ==============================================================  ======================================  ============
 Specify that the implicit RHS is linear                         :c:func:`ARKodeSetLinear`               ``SUNFALSE``
 Specify that the implicit RHS nonlinear                         :c:func:`ARKodeSetNonlinear`            ``SUNTRUE``
+Specify that the implicit RHS is autonomous                     :c:func:`ARKodeSetAutonomous`           ``SUNFALSE``
 Implicit predictor method                                       :c:func:`ARKodeSetPredictorMethod`      0
 User-provided implicit stage predictor                          :c:func:`ARKodeSetStagePredictFn`       ``NULL``
 RHS function for nonlinear system evaluations                   :c:func:`ARKodeSetNlsRhsFn`             ``NULL``
@@ -1843,6 +1844,58 @@ Specify if the implicit RHS is deduced after a nonlinear solve  :c:func:`ARKodeS
    .. versionadded:: x.y.z
 
 
+.. c:function:: int ARKodeSetAutonomous(void* arkode_mem, sunbooleantype autonomous)
+
+   Specifies that the implicit portion of the problem is autonomous i.e., does
+   not explicitly depend on time.
+
+   When using an implicit or ImEx method with the trivial predictor, this option
+   enables reusing the implicit right-hand side evaluation at the predicted
+   state across stage solves within a step. This reuse reduces the total number
+   of implicit RHS function evaluations.
+
+   :param arkode_mem: pointer to the ARKODE memory block.
+   :param autonomous: flag denoting if the implicit RHS function,
+                      :math:`f^I(t,y)`, is autonomous (``SUNTRUE``) or
+                      non-autonomous (``SUNFALSE``).
+
+   :retval ARK_SUCCESS: the function exited successfully.
+   :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL``.
+   :retval ARK_ILL_INPUT: an argument had an illegal value.
+   :retval ARK_STEPPER_UNSUPPORTED: implicit solvers are not supported by the
+                                    current time-stepping module.
+
+   .. warning::
+
+      Results may differ when enabling both :c:func:`ARKodeSetAutonomous` and
+      :c:func:`ARKodeSetDeduceImplicitRhs` with a stiffly accurate implicit
+      method and using the trivial predictor. The differences are due to reusing
+      the deduced implicit right-hand side (RHS) value in the initial nonlinear
+      residual computation rather than evaluating the implicit RHS function. The
+      significance of the difference will depend on how well the deduced RHS
+      approximates the RHS evaluated at the trivial predictor. This behavior can
+      be observed in ``examples/arkode/C_serial/ark_brusselator.c`` by comparing
+      the outputs with :c:func:`ARKodeSetAutonomous` enabled/disabled.
+
+      Similarly programs that assume the nonlinear residual will always call the
+      implicit RHS function will need to be updated to account for the RHS value
+      reuse when using :c:func:`ARKodeSetAutonomous`. For example,
+      ``examples/arkode/C_serial/ark_KrylovDemo_prec.c`` assumes that the
+      nonlinear residual will be called and will evaluate the implicit RHS
+      function before calling the preconditioner setup function. Based on this
+      assumption, this example code saves some computations in the RHS
+      evaluation for reuse in the preconditioner setup. However, when
+      :c:func:`ARKodeSetAutonomous` is enabled, the call to the nonlinear
+      residual before the preconditioner setup reuses a saved RHS evaluation and
+      the saved data is actually from an earlier RHS evaluation that is not
+      consistent with the state and RHS values passed to the preconditioner
+      setup function. For this example, the code should not save data in the RHS
+      evaluation but instead evaluate the necessary quantities within the
+      preconditioner setup function using the input values.
+
+   .. versionadded:: x.y.z
+
+
 .. c:function:: int ARKodeSetPredictorMethod(void* arkode_mem, int method)
 
    Specifies the method from :numref:`ARKODE.Mathematics.Predictors` to use

doc/shared/RecentChanges.rst

@@ -36,6 +36,12 @@ Added :c:func:`ARKodeButcherTable_ERKIDToName` and
 :c:func:`ARKodeButcherTable_DIRKIDToName` to convert a Butcher table ID to a
 string representation.
 
+Added the function :c:func:`ARKodeSetAutonomous` in ARKODE to indicate that the
+implicit right-hand side function does not explicitly depend on time. When using
+the trivial predictor, an autonomous problem may reuse implicit function
+evaluations across stage solves to reduce the total number of function
+evaluations.
+
 Users may now disable interpolated output in ARKODE by passing
 ``ARK_INTERP_NONE`` to :c:func:`ARKodeSetInterpolantType`. When interpolation is
 disabled, rootfinding is not supported, implicit methods must use the trivial

examples/arkode/CXX_serial/ark_analytic_sys.cpp

@@ -162,6 +162,7 @@ int main()
   // Initialize dense matrix data structure and solver
   A = SUNDenseMatrix(NEQ, NEQ, sunctx);
   if (check_flag((void*)A, "SUNDenseMatrix", 0)) { return 1; }
+
   LS = SUNLinSol_Dense(y, A, sunctx);
   if (check_flag((void*)LS, "SUNLinSol_Dense", 0)) { return 1; }
 
@@ -176,20 +177,25 @@ int main()
   flag = ARKodeSetUserData(arkode_mem,
                            (void*)&lamda); // Pass lamda to user functions
   if (check_flag(&flag, "ARKodeSetUserData", 1)) { return 1; }
+
   flag = ARKodeSStolerances(arkode_mem, reltol, abstol); // Specify tolerances
   if (check_flag(&flag, "ARKodeSStolerances", 1)) { return 1; }
 
   // Linear solver interface
   flag = ARKodeSetLinearSolver(arkode_mem, LS,
                                A); // Attach matrix and linear solver
   if (check_flag(&flag, "ARKodeSetLinearSolver", 1)) { return 1; }
+
   flag = ARKodeSetJacFn(arkode_mem, Jac); // Set Jacobian routine
   if (check_flag(&flag, "ARKodeSetJacFn", 1)) { return 1; }
 
   // Specify linearly implicit RHS, with non-time-dependent Jacobian
   flag = ARKodeSetLinear(arkode_mem, 0);
   if (check_flag(&flag, "ARKodeSetLinear", 1)) { return 1; }
 
+  flag = ARKodeSetAutonomous(arkode_mem, SUNTRUE);
+  if (check_flag(&flag, "ARKodeSetAutonomous", 1)) { return 1; }
+
   // Open output stream for results, output comment line
   FILE* UFID = fopen("solution.txt", "w");
   fprintf(UFID, "# t y1 y2 y3\n");

examples/arkode/CXX_serial/ark_analytic_sys.out

@@ -20,7 +20,7 @@ Analytical ODE test problem:
 
 Final Solver Statistics:
    Internal solver steps = 67 (attempted = 69)
-   Total RHS evals:  Fe = 0,  Fi = 693
+   Total RHS evals:  Fe = 0,  Fi = 348
    Total linear solver setups = 22
    Total RHS evals for setting up the linear system = 0
    Total number of Jacobian evaluations = 3

examples/arkode/CXX_serial/ark_pendulum.cpp

@@ -213,6 +213,9 @@ int main(int argc, char* argv[])
   flag = ARKodeSetNonlinConvCoef(arkode_mem, SUN_RCONST(0.01));
   if (check_flag(flag, "ARKodeSetNonlinConvCoef")) { return 1; }
 
+  flag = ARKodeSetAutonomous(arkode_mem, SUNTRUE);
+  if (check_flag(flag, "ARKodeSetAutonomous")) { return 1; }
+
   /* --------------- *
    * Advance in Time *
    * --------------- */

examples/arkode/CXX_serial/ark_pendulum.out

@@ -20,7 +20,7 @@ Nonlinear Pendulum problem:
 Final Solver Statistics:
   Internal solver steps = 8552 (attempted = 8562)
   Total number of error test failures = 10
-  Total RHS evals:  Fe = 0,  Fi = 52696
+  Total RHS evals:  Fe = 0,  Fi = 35561
   Total number of Newton iterations = 27018
   Total number of linear solver convergence failures = 11
   Total linear solver setups = 471

examples/arkode/C_openmp/ark_brusselator1D_omp.c

@@ -225,6 +225,7 @@ int main(int argc, char* argv[])
   /* Initialize matrix and linear solver data structures */
   A = SUNBandMatrix(NEQ, 4, 4, ctx);
   if (check_flag((void*)A, "SUNBandMatrix", 0)) { return 1; }
+
   LS = SUNLinSol_Band(y, A, ctx);
   if (check_flag((void*)LS, "SUNLinSol_Band", 0)) { return 1; }
 
@@ -239,16 +240,21 @@ int main(int argc, char* argv[])
   flag = ARKodeSetUserData(arkode_mem,
                            (void*)udata); /* Pass udata to user functions */
   if (check_flag(&flag, "ARKodeSetUserData", 1)) { return 1; }
+
   flag = ARKodeSStolerances(arkode_mem, reltol, abstol); /* Specify tolerances */
   if (check_flag(&flag, "ARKodeSStolerances", 1)) { return 1; }
 
   /* Linear solver specification */
   flag = ARKodeSetLinearSolver(arkode_mem, LS,
                                A); /* Attach matrix and linear solver */
   if (check_flag(&flag, "ARKodeSetLinearSolver", 1)) { return 1; }
+
   flag = ARKodeSetJacFn(arkode_mem, Jac); /* Set the Jacobian routine */
   if (check_flag(&flag, "ARKodeSetJacFn", 1)) { return 1; }
 
+  flag = ARKodeSetAutonomous(arkode_mem, SUNTRUE);
+  if (check_flag(&flag, "ARKodeSetAutonomous", 1)) { return 1; }
+
   /* output spatial mesh to disk */
   FID = fopen("bruss_mesh.txt", "w");
   for (i = 0; i < N; i++) { fprintf(FID, "  %.16" ESYM "\n", udata->dx * i); }

examples/arkode/C_openmp/ark_brusselator1D_omp.out

@@ -112,7 +112,7 @@
 
 Final Solver Statistics:
    Internal solver steps = 99 (attempted = 99)
-   Total RHS evals:  Fe = 0,  Fi = 1694
+   Total RHS evals:  Fe = 0,  Fi = 1186
    Total linear solver setups = 35
    Total RHS evals for setting up the linear system = 0
    Total number of Jacobian evaluations = 14

examples/arkode/C_serial/ark_KrylovDemo_prec.c

@@ -193,6 +193,13 @@ typedef struct
   int jxr[NGX], jyr[NGY];
   sunrealtype acoef[NS][NS], bcoef[NS], diff[NS];
   sunrealtype cox[NS], coy[NS], dx, dy, srur;
+  /* fsave contains saved rate data from the last RHS evaluation for reuse in
+     the preconditioner setup. This reuse relies on the nonlinear residual
+     evaluating the RHS function before calling the preconditioner setup
+     function. This assumption is incorrect when using ARKodeSetAutonomous which
+     enables reuse of RHS evaluations with the trivial predictor. To use
+     ARKodeSetAutonomous in this example, the saved values would need to be
+     recomputed in the preconditioner setup function using the input state. */
   sunrealtype fsave[NEQ];
   N_Vector tmp;
   N_Vector rewt;

examples/arkode/C_serial/ark_brusselator.c

@@ -171,27 +171,40 @@ int main(void)
   flag = ARKodeSetUserData(arkode_mem,
                            (void*)rdata); /* Pass rdata to user functions */
   if (check_flag(&flag, "ARKodeSetUserData", 1)) { return 1; }
+
   flag = ARKodeSStolerances(arkode_mem, reltol, abstol); /* Specify tolerances */
   if (check_flag(&flag, "ARKodeSStolerances", 1)) { return 1; }
+
   flag = ARKodeSetInterpolantType(arkode_mem,
                                   ARK_INTERP_LAGRANGE); /* Specify stiff interpolant */
   if (check_flag(&flag, "ARKodeSetInterpolantType", 1)) { return 1; }
+
   flag = ARKodeSetDeduceImplicitRhs(arkode_mem, 1); /* Avoid eval of f after stage */
   if (check_flag(&flag, "ARKodeSetDeduceImplicitRhs", 1)) { return 1; }
 
   /* Initialize dense matrix data structure and solver */
   A = SUNDenseMatrix(NEQ, NEQ, ctx);
   if (check_flag((void*)A, "SUNDenseMatrix", 0)) { return 1; }
+
   LS = SUNLinSol_Dense(y, A, ctx);
   if (check_flag((void*)LS, "SUNLinSol_Dense", 0)) { return 1; }
 
   /* Linear solver interface */
   flag = ARKodeSetLinearSolver(arkode_mem, LS,
                                A); /* Attach matrix and linear solver */
   if (check_flag(&flag, "ARKodeSetLinearSolver", 1)) { return 1; }
+
   flag = ARKodeSetJacFn(arkode_mem, Jac); /* Set Jacobian routine */
   if (check_flag(&flag, "ARKodeSetJacFn", 1)) { return 1; }
 
+  /* Signal that this problem does not explicitly depend on time. In this case
+     enabling/disabling this option will lead to slightly different results when
+     combined with ARKodeSetDeduceImplicitRhs because the implicit method is
+     stiffly accurate. The differences are due to reuse of the deduced implicit
+     RHS evaluation in the nonlinear residual. */
+  flag = ARKodeSetAutonomous(arkode_mem, SUNTRUE);
+  if (check_flag(&flag, "ARKodeSetAutonomous", 1)) { return 1; }
+
   /* Open output stream for results, output comment line */
   UFID = fopen("solution.txt", "w");
   fprintf(UFID, "# t u v w\n");

examples/arkode/C_serial/ark_brusselator.out

@@ -8,24 +8,24 @@ Brusselator ODE test problem:
    -------------------------------------------
     0.000000    1.200000    3.100000    3.000000
     1.000000    1.103849    3.013164    3.499981
-    2.000000    0.687998    3.521383    3.499988
-    3.000000    0.409468    4.277886    3.499993
-    4.000000    0.367888    4.942004    3.499994
-    5.000000    0.413854    5.510621    3.499993
-    6.000000    0.589240    5.855669    3.499990
-    7.000000    4.756630    0.735405    3.499917
-    8.000000    1.813466    1.575769    3.499968
-    9.000000    0.527904    2.807358    3.499991
-   10.000000    0.305600    3.657365    3.499995
+    2.000000    0.687997    3.521384    3.499988
+    3.000000    0.409467    4.277888    3.499993
+    4.000000    0.367889    4.942005    3.499994
+    5.000000    0.413860    5.510615    3.499993
+    6.000000    0.589251    5.855655    3.499990
+    7.000000    4.756472    0.735414    3.499917
+    8.000000    1.813403    1.575799    3.499968
+    9.000000    0.527886    2.807385    3.499991
+   10.000000    0.305598    3.657381    3.499995
    -------------------------------------------
 
 Final Solver Statistics:
-   Internal solver steps = 250 (attempted = 262)
-   Total RHS evals:  Fe = 0,  Fi = 3282
-   Total linear solver setups = 118
+   Internal solver steps = 252 (attempted = 262)
+   Total RHS evals:  Fe = 0,  Fi = 1910
+   Total linear solver setups = 111
    Total RHS evals for setting up the linear system = 0
-   Total number of Jacobian evaluations = 72
-   Total number of Newton iterations = 3279
-   Total number of nonlinear solver convergence failures = 71
-   Total number of error test failures = 11
+   Total number of Jacobian evaluations = 71
+   Total number of Newton iterations = 3286
+   Total number of nonlinear solver convergence failures = 70
+   Total number of error test failures = 9
    Total number of failed steps from solver failure = 1

examples/arkode/C_serial/ark_brusselator1D.c

@@ -216,16 +216,21 @@ int main(void)
   /* Initialize band matrix data structure and solver -- A will be factored, so set smu to ml+mu */
   A = SUNBandMatrix(NEQ, 4, 4, ctx);
   if (check_flag((void*)A, "SUNBandMatrix", 0)) { return 1; }
+
   LS = SUNLinSol_Band(y, A, ctx);
   if (check_flag((void*)LS, "SUNLinSol_Band", 0)) { return 1; }
 
   /* Linear solver interface */
   flag = ARKodeSetLinearSolver(arkode_mem, LS,
                                A); /* Attach matrix and linear solver */
   if (check_flag(&flag, "ARKodeSetLinearSolver", 1)) { return 1; }
+
   flag = ARKodeSetJacFn(arkode_mem, Jac); /* Set the Jacobian routine */
   if (check_flag(&flag, "ARKodeSetJacFn", 1)) { return 1; }
 
+  flag = ARKodeSetAutonomous(arkode_mem, SUNTRUE);
+  if (check_flag(&flag, "ARKodeSetAutonomous", 1)) { return 1; }
+
   /* output spatial mesh to disk */
   FID = fopen("bruss_mesh.txt", "w");
   for (i = 0; i < N; i++) { fprintf(FID, "  %.16" ESYM "\n", udata->dx * i); }

examples/arkode/C_serial/ark_brusselator1D.out

@@ -111,7 +111,7 @@
 
 Final Solver Statistics:
    Internal solver steps = 99 (attempted = 99)
-   Total RHS evals:  Fe = 0,  Fi = 1694
+   Total RHS evals:  Fe = 0,  Fi = 1186
    Total linear solver setups = 35
    Total RHS evals for setting up the linear system = 0
    Total number of Jacobian evaluations = 14

examples/arkode/C_serial/ark_brusselator1D_FEM_slu.c

@@ -296,6 +296,10 @@ int main(int argc, char* argv[])
   retval = ARKodeResStolerance(arkode_mem, abstol);
   if (check_retval(&retval, "ARKodeResStolerance", 1)) { return (1); }
 
+  /* Specify an autonomous problem */
+  retval = ARKodeSetAutonomous(arkode_mem, SUNTRUE);
+  if (check_retval(&retval, "ARKodeSetAutonomous", 1)) { return (1); }
+
   /* Initialize sparse matrix data structure and linear solvers (system and mass) */
   NNZ = 15 * NEQ;
 

examples/arkode/C_serial/ark_brusselator1D_FEM_slu.out

@@ -112,7 +112,7 @@
 
 Final Solver Statistics:
    Internal solver steps = 95 (attempted = 95)
-   Total RHS evals:  Fe = 0,  Fi = 1596
+   Total RHS evals:  Fe = 0,  Fi = 1111
    Total mass matrix setups = 1
    Total mass matrix solves = 193
    Total mass times evals = 1688

examples/arkode/C_serial/ark_brusselator1D_klu.c

@@ -233,15 +233,20 @@ int main(void)
   NNZ = 5 * NEQ;
   A   = SUNSparseMatrix(NEQ, NEQ, NNZ, CSC_MAT, ctx);
   if (check_flag((void*)A, "SUNSparseMatrix", 0)) { return 1; }
+
   LS = SUNLinSol_KLU(y, A, ctx);
   if (check_flag((void*)LS, "SUNLinSol_KLU", 0)) { return 1; }
 
   /* Attach the matrix, linear solver, and Jacobian construction routine to ARKODE */
   flag = ARKodeSetLinearSolver(arkode_mem, LS, A); /* Attach matrix and LS */
   if (check_flag(&flag, "ARKodeSetLinearSolver", 1)) { return 1; }
+
   flag = ARKodeSetJacFn(arkode_mem, Jac); /* Supply Jac routine */
   if (check_flag(&flag, "ARKodeSetJacFn", 1)) { return 1; }
 
+  flag = ARKodeSetAutonomous(arkode_mem, SUNTRUE);
+  if (check_flag(&flag, "ARKodeSetAutonomous", 1)) { return 1; }
+
   /* output spatial mesh to disk */
   FID = fopen("bruss_mesh.txt", "w");
   for (i = 0; i < N; i++) { fprintf(FID, "  %.16" ESYM "\n", udata->dx * i); }

examples/arkode/C_serial/ark_brusselator1D_klu.out

@@ -21,7 +21,7 @@
 
 Final Solver Statistics:
    Internal solver steps = 99 (attempted = 99)
-   Total RHS evals:  Fe = 0,  Fi = 1694
+   Total RHS evals:  Fe = 0,  Fi = 1186
    Total linear solver setups = 35
    Total number of Jacobian evaluations = 14
    Total number of nonlinear iterations = 1196

examples/arkode/C_serial/ark_brusselator_fp.c

@@ -196,12 +196,17 @@ int main(int argc, char* argv[])
   flag = ARKodeSetUserData(arkode_mem,
                            (void*)rdata); /* Pass rdata to user functions */
   if (check_flag(&flag, "ARKodeSetUserData", 1)) { return 1; }
+
   flag = ARKodeSStolerances(arkode_mem, reltol, abstol); /* Specify tolerances */
   if (check_flag(&flag, "ARKodeSStolerances", 1)) { return 1; }
+
   flag = ARKodeSetMaxNonlinIters(arkode_mem,
                                  maxcor); /* Increase default iterations */
   if (check_flag(&flag, "ARKodeSetMaxNonlinIters", 1)) { return 1; }
 
+  flag = ARKodeSetAutonomous(arkode_mem, SUNTRUE);
+  if (check_flag(&flag, "ARKodeSetAutonomous", 1)) { return 1; }
+
   /* Open output stream for results, output comment line */
   UFID = fopen("solution.txt", "w");
   fprintf(UFID, "# t u v w\n");

examples/arkode/C_serial/ark_brusselator_fp.out

@@ -20,7 +20,7 @@ Brusselator ODE test problem, fixed-point solver:
 
 Final Solver Statistics:
    Internal solver steps = 729 (attempted = 730)
-   Total RHS evals:  Fe = 4382,  Fi = 18792
+   Total RHS evals:  Fe = 4382,  Fi = 15142
    Total number of fixed-point iterations = 14410
    Total number of nonlinear solver convergence failures = 0
    Total number of error test failures = 1