Feature: Remove Implicit RHS Evaluation in Autonomous Problems With t…

…he Trivial Predictor (#495)

Add the function ``ARKodeSetAutonomous`` 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 reducing the total number of function
evaluations.

---------

Co-authored-by: Daniel R. Reynolds <reynolds@smu.edu>
Co-authored-by: Cody Balos <balos1@llnl.gov>

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