Initial discrete adjoint capability with ERK

.clangd

@@ -0,0 +1,2 @@
+CompileFlags:
+  CompilationDatabase: builddir/

doc/arkode/guide/source/Mathematics.rst

@@ -2126,3 +2126,65 @@ by :math:`\eta_\text{rf}`.
 
 For more information on utilizing relaxation Runge--Kutta methods, see
 :numref:`ARKODE.Usage.Relaxation`.
+
+
+.. _ARKODE.Mathematics.ASA:
+
+Adjoint Sensitivity Analysis
+============================
+
+Consider :eq:`ARKODE_IVP_simple_explicit`, but where the ODE also depends on some parameters
+:math:`p` (that is, we have :math:`f(t,y,p)`). Now, let :math:`g(y(t_f),p)` be a cost function for
+which we would like to compute the gradients or :math:`\partial g/\partial y(t_0)` and/or
+:math:`\partial g/\partial p`. The adjoint method is one approach to obtaining the gradients that is
+particularly efficient when there are relatively few cost functionals and a large number of parameters.
+With the adjoint method we solve the adjoint ODEs for :math:`\lambda(t) \in \mathbb{R}^N` and
+:math:`\mu(t) \in \mathbb{R}^{N_s}`:
+
+.. math::
+   \lambda'(t) &= -f_y^T(t, y, p) \lambda,\quad \lambda(t_F) = g_y^T(y(t_f), p) \\
+   \mu'(t) &= -f_p^T(t, y, p) \mu,\quad \mu(t_F) = g_p^T(y(t_f), p), \quad t_f \geq t \geq t_0. \\
+
+(For a detailed derivation see :cite:p:`hager2000runge,sanduDiscrete2006`). Here :math:`f_y \equiv
+\partial f/\partial y` is the Jacobian with respect to the dependent variable and :math:`f_p \equiv
+\partial f/\partial p` is the Jacobian with respect to the parameters. The ARKStep module in ARKODE
+provides adjoint sensitivity analysis based on the *discrete* formulation. Given a s-stage explicit
+Runge--Kutta method (as in :eq:`ARKODE_ERK`, but without the embedding), the discrete adjoint works
+out to be
+
+.. math::
+   \Lambda_i &= h_n f_y^T(t_{n,i}, z_i) \left(b_i \lambda_{n+1} + \sum_{j=i+1}^s a_{j,i}
+   \Lambda_j \right), \quad \quad i = s, \dots, 1,\\
+   \nu_i     &= h_n f_p^T(t_{n,i}, z_i, p) \left(b_i \lambda_{n+1} + \sum_{j=i}^{s} a_{ji} \nu_j \right) \\
+   \lambda_n &= \lambda_{n+1} + \sum_{j=1}^{s} \Lambda_j \\
+   \mu_i     &= \mu_{n+1} + \sum_{j=1}^{s} \nu_j.
+   :label: ARKODE_ERK_Adjoint
+
+The result of intergating the adjoint equations with the discrete adjoint of the ARKStep method is
+then
+
+.. math::
+   \frac{\partial g}{\partial y_{(t_0)}} = \lambda_i, \quad
+   \frac{\partial g}{\partial p} = \mu_i + \lambda_i \left(\frac{\partial y(t_0)}{\partial p} \right).
+
+
+Discrete vs. Continuous Adjoint Method
+--------------------------------------
+
+At a high level there are two adjoint method approaches that one can use. In the *continuous*
+approach, we derive the sensitivity equations directly from the model and then we integrate them
+with a time integration method. This is the approach implemented in the SUNDIALS :ref:`CVODES
+<CVODES.Mathematics.ASA>` and :ref:`IDAS <IDAS.Mathematics.ASA>` packages. In the *discrete*
+approach, the model equations are discretized with the time integration method first, and then we
+derive the adjoints of the discretized equations. It is understood that the continuous adjoint
+method can be problematic in the context of optimization problems because the continuous adjoint
+method provides an approximation to the gradient of a continuous cost function while the optimizer
+is expecting the gradient of the discrete cost function. The discrepancy means that the optimizer
+can fail to converge further once it is near a local minimum :cite:p:`giles2000introduction`. One
+the other hand, the discrete adjoint method provides the exact gradient of the discrete cost
+function allowing the optimizer to fully converge. Consequently, the discrete adjoint method is
+often preferable in optimization despite its own drawbacks -- such as its (relatively) increased
+memory usage and the possible introduction of unphysical computational modes
+:cite:p:`sirkes1997finite`. This is not to say that the discrete adjoint method is always the better
+choice over the continuous adjoint method in optimization. Practical considerations may lead one to
+choose the continuous approach.

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

@@ -867,29 +867,30 @@ Optional inputs for ARKODE
 
 .. cssclass:: table-bordered
 
-================================================  =======================================  =======================
-Optional input                                    Function name                            Default
-================================================  =======================================  =======================
-Return ARKODE parameters to their defaults        :c:func:`ARKodeSetDefaults`              internal
-Set integrator method order                       :c:func:`ARKodeSetOrder`                 4
-Set dense output interpolation type (SPRKStep)    :c:func:`ARKodeSetInterpolantType`       ``ARK_INTERP_LAGRANGE``
-Set dense output interpolation type (others)      :c:func:`ARKodeSetInterpolantType`       ``ARK_INTERP_HERMITE``
-Set dense output polynomial degree                :c:func:`ARKodeSetInterpolantDegree`     5
-Disable time step adaptivity (fixed-step mode)    :c:func:`ARKodeSetFixedStep`             disabled
-Supply an initial step size to attempt            :c:func:`ARKodeSetInitStep`              estimated
-Maximum no. of warnings for :math:`t_n+h = t_n`   :c:func:`ARKodeSetMaxHnilWarns`          10
-Maximum no. of internal steps before *tout*       :c:func:`ARKodeSetMaxNumSteps`           500
-Maximum absolute step size                        :c:func:`ARKodeSetMaxStep`               :math:`\infty`
-Minimum absolute step size                        :c:func:`ARKodeSetMinStep`               0.0
-Set a value for :math:`t_{stop}`                  :c:func:`ARKodeSetStopTime`              undefined
-Interpolate at :math:`t_{stop}`                   :c:func:`ARKodeSetInterpolateStopTime`   ``SUNFALSE``
-Disable the stop time                             :c:func:`ARKodeClearStopTime`            N/A
-Supply a pointer for user data                    :c:func:`ARKodeSetUserData`              ``NULL``
-Maximum no. of ARKODE error test failures         :c:func:`ARKodeSetMaxErrTestFails`       7
-Set inequality constraints on solution            :c:func:`ARKodeSetConstraints`           ``NULL``
-Set max number of constraint failures             :c:func:`ARKodeSetMaxNumConstrFails`     10
-================================================  =======================================  =======================
-
+=================================================  =======================================  =======================
+Optional input                                     Function name                            Default
+=================================================  =======================================  =======================
+Return ARKODE parameters to their defaults         :c:func:`ARKodeSetDefaults`              internal
+Set integrator method order                        :c:func:`ARKodeSetOrder`                 4
+Set dense output interpolation type (SPRKStep)     :c:func:`ARKodeSetInterpolantType`       ``ARK_INTERP_LAGRANGE``
+Set dense output interpolation type (others)       :c:func:`ARKodeSetInterpolantType`       ``ARK_INTERP_HERMITE``
+Set dense output polynomial degree                 :c:func:`ARKodeSetInterpolantDegree`     5
+Disable time step adaptivity (fixed-step mode)     :c:func:`ARKodeSetFixedStep`             disabled
+Supply an initial step size to attempt             :c:func:`ARKodeSetInitStep`              estimated
+Maximum no. of warnings for :math:`t_n+h = t_n`    :c:func:`ARKodeSetMaxHnilWarns`          10
+Maximum no. of internal steps before *tout*        :c:func:`ARKodeSetMaxNumSteps`           500
+Maximum absolute step size                         :c:func:`ARKodeSetMaxStep`               :math:`\infty`
+Minimum absolute step size                         :c:func:`ARKodeSetMinStep`               0.0
+Set a value for :math:`t_{stop}`                   :c:func:`ARKodeSetStopTime`              undefined
+Interpolate at :math:`t_{stop}`                    :c:func:`ARKodeSetInterpolateStopTime`   ``SUNFALSE``
+Disable the stop time                              :c:func:`ARKodeClearStopTime`            N/A
+Supply a pointer for user data                     :c:func:`ARKodeSetUserData`              ``NULL``
+Maximum no. of ARKODE error test failures          :c:func:`ARKodeSetMaxErrTestFails`       7
+Set inequality constraints on solution             :c:func:`ARKodeSetConstraints`           ``NULL``
+Set max number of constraint failures              :c:func:`ARKodeSetMaxNumConstrFails`     10
+Set the checkpointing scheme to use (for adjoint)  :c:func:`ARKodeSetCheckpointScheme`      ``NULL``
+Set the checkpointing step index (for adjoint)     :c:func:`ARKodeSetCheckpointIndex        0
+=================================================  =======================================  =======================
 
 
 
@@ -1397,6 +1398,34 @@ Set max number of constraint failures             :c:func:`ARKodeSetMaxNumConstr
    .. versionadded:: 6.1.0
 
 
+.. c:function:: int ARKodeSetCheckpointScheme(void* arkode_mem, SUNAdjointCheckpointScheme checkpoint_scheme)
+
+   Specifies the :c:type:`SUNAdjointCheckpointScheme` to use for saving states
+   during the forward integration, and loading states during backward integration
+   of an adjoint system.
+
+   :param arkode_mem: pointer to the ARKODE memory block.
+   :param checkpoint_scheme: the checkpoint scheme to use
+
+   :retval ARK_SUCCESS: the function exited successfully.
+   :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL``.
+
+   .. versionadded:: x.x.x
+
+.. c:function:: int ARKodeSetCheckpointIndex(void* arkode_mem, long int step_index)
+
+   Specifies the step index (that is step number) to insert the next checkpoint at.
+   This is incremented along with the step count, but it is useful to be able to reset
+   this index during recomputations of missing states during the backward adjoint integration.
+
+   :param arkode_mem: pointer to the ARKODE memory block.
+   :param checkpoint_scheme: the checkpoint scheme to use
+
+   :retval ARK_SUCCESS: the function exited successfully.
+   :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL``.
+
+   .. versionadded:: x.x.x
+
 
 .. _ARKODE.Usage.ARKodeAdaptivityInputTable:
 

doc/arkode/guide/source/index.rst

@@ -65,6 +65,7 @@ with support by the `US Department of Energy <http://www.doe.gov>`_,
    sunlinsol/index.rst
    sunnonlinsol/index.rst
    sunadaptcontroller/index.rst
+   sunadjoint/index.rst
    sunmemory/index.rst
    sundials/Install_link.rst
    Constants

doc/arkode/guide/source/sunadjoint/SUNAdjoint_links.rst

@@ -0,0 +1,14 @@
+.. ----------------------------------------------------------------
+   SUNDIALS Copyright Start
+   Copyright (c) 2002-2024, Lawrence Livermore National Security
+   and Southern Methodist University.
+   All rights reserved.
+
+   See the top-level LICENSE and NOTICE files for details.
+
+   SPDX-License-Identifier: BSD-3-Clause
+   SUNDIALS Copyright End
+   ----------------------------------------------------------------
+
+.. include:: ../../../shared/sunadjoint/SUNAdjointCheckpointScheme.rst
+.. include:: ../../../shared/sunadjoint/SUNAdjointStepper.rst

doc/arkode/guide/source/sunadjoint/index.rst

@@ -0,0 +1,31 @@
+..
+   ----------------------------------------------------------------
+   SUNDIALS Copyright Start
+   Copyright (c) 2002-2024, Lawrence Livermore National Security
+   and Southern Methodist University.
+   All rights reserved.
+
+   See the top-level LICENSE and NOTICE files for details.
+
+   SPDX-License-Identifier: BSD-3-Clause
+   SUNDIALS Copyright End
+   ----------------------------------------------------------------
+
+.. _SUNAdaptController:
+
+###########################
+Adjoint Sensitivty Analysis
+###########################
+
+.. versionadded:: x.x.x
+
+The ``SUNAdjoint`` API consists of a few customizable modules that provide a framework for adjoint
+sensitivity analysis (ASA). The API itself does not implement ASA, but it provides a common
+interface for ASA capabilties implemented in the SUNDIALS packages. Right now it supports :ref:`the
+ASA capabilties in ARKODE <ARKODE.Mathematics.ASA>`, while the ASA capabilities in :ref:`CVODES
+<CVODES.Mathematics.ASA>` and :ref:`IDAS <IDAS.Mathematics.ASA>` must be used directly.
+
+.. toctree::
+   :maxdepth: 1
+
+   SUNAdjoint_links.rst

doc/shared/sunadjoint/SUNAdjointCheckpointScheme.rst

@@ -0,0 +1,17 @@
+..
+   ----------------------------------------------------------------
+   SUNDIALS Copyright Start
+   Copyright (c) 2002-2024, Lawrence Livermore National Security
+   and Southern Methodist University.
+   All rights reserved.
+
+   See the top-level LICENSE and NOTICE files for details.
+
+   SPDX-License-Identifier: BSD-3-Clause
+   SUNDIALS Copyright End
+   ----------------------------------------------------------------
+
+.. _SUNAdjoint.CheckpointScheme:
+
+The SUNAdjointCheckpointScheme API
+==================================

doc/shared/sunadjoint/SUNAdjointStepper.rst

@@ -0,0 +1,17 @@
+..
+   ----------------------------------------------------------------
+   SUNDIALS Copyright Start
+   Copyright (c) 2002-2024, Lawrence Livermore National Security
+   and Southern Methodist University.
+   All rights reserved.
+
+   See the top-level LICENSE and NOTICE files for details.
+
+   SPDX-License-Identifier: BSD-3-Clause
+   SUNDIALS Copyright End
+   ----------------------------------------------------------------
+
+.. _SUNAdjoint.Solver:
+
+The SUNAdjointStepper API
+========================

doc/shared/sunadjoint/SUNAdjoint_Description.rst

@@ -0,0 +1,14 @@
+..
+   ----------------------------------------------------------------
+   SUNDIALS Copyright Start
+   Copyright (c) 2002-2024, Lawrence Livermore National Security
+   and Southern Methodist University.
+   All rights reserved.
+
+   See the top-level LICENSE and NOTICE files for details.
+
+   SPDX-License-Identifier: BSD-3-Clause
+   SUNDIALS Copyright End
+   ----------------------------------------------------------------
+
+.. _SUNAdjoint.Description:

doc/shared/sundials.bib

@@ -2369,3 +2369,48 @@ @article{edwards2014kokkos
   issn    = {0743-7315},
   doi     = {10.1016/j.jpdc.2014.07.003}
 }
+
+%
+% Discrete adjoints
+%
+
+@article{giles2000introduction,
+  title={An introduction to the adjoint approach to design},
+  author={Giles, Michael B and Pierce, Niles A},
+  journal={Flow, turbulence and combustion},
+  volume={65},
+  number={3},
+  pages={393--415},
+  year={2000},
+  publisher={Springer}
+}
+
+@article{sirkes1997finite,
+  title={Finite difference of adjoint or adjoint of finite difference?},
+  author={Sirkes, Ziv and Tziperman, Eli},
+  journal={Monthly weather review},
+  volume={125},
+  number={12},
+  pages={3373--3378},
+  year={1997}
+}
+
+@article{hager2000runge,
+  title={Runge-Kutta methods in optimal control and the transformed adjoint system},
+  author={Hager, William W},
+  journal={Numerische Mathematik},
+  volume={87},
+  pages={247--282},
+  year={2000},
+  publisher={Springer}
+}
+
+@article{sanduDiscrete2006,
+  year = {2006},
+  title = {{On the Properties of Runge-Kutta Discrete Adjoints}},
+  author = {Sandu, Adrian},
+  journal = {Lecture Notes in Computer Science},
+  issn = {0302-9743},
+  doi = {10.1007/11758549_76},
+  pages = {550--557}
+}

doc/superbuild/source/index.rst

@@ -174,6 +174,7 @@ SUNDIALS License and Notices
    sunlinsol/index.rst
    sunnonlinsol/index.rst
    sunadaptcontroller/index.rst
+   sunadjoint/index.rst
    sunmemory/index.rst
    History_link.rst
    Changelog_link.rst

doc/superbuild/source/sunadjoint/SUNAdjoint_links.rst

@@ -0,0 +1,14 @@
+.. ----------------------------------------------------------------
+   SUNDIALS Copyright Start
+   Copyright (c) 2002-2024, Lawrence Livermore National Security
+   and Southern Methodist University.
+   All rights reserved.
+
+   See the top-level LICENSE and NOTICE files for details.
+
+   SPDX-License-Identifier: BSD-3-Clause
+   SUNDIALS Copyright End
+   ----------------------------------------------------------------
+
+.. include:: ../../../shared/sunadjoint/SUNAdjointCheckpointScheme.rst
+.. include:: ../../../shared/sunadjoint/SUNAdjointStepper.rst

doc/superbuild/source/sunadjoint/index.rst

@@ -0,0 +1,31 @@
+..
+   ----------------------------------------------------------------
+   SUNDIALS Copyright Start
+   Copyright (c) 2002-2024, Lawrence Livermore National Security
+   and Southern Methodist University.
+   All rights reserved.
+
+   See the top-level LICENSE and NOTICE files for details.
+
+   SPDX-License-Identifier: BSD-3-Clause
+   SUNDIALS Copyright End
+   ----------------------------------------------------------------
+
+.. _SUNAdaptController:
+
+###########################
+Adjoint Sensitivty Analysis
+###########################
+
+.. versionadded:: x.x.x
+
+The ``SUNAdjoint`` API consists of a few customizable modules that provide a framework for adjoint
+sensitivity analysis (ASA). The API itself does not implement ASA, but it provides a common
+interface for ASA capabilties implemented in the SUNDIALS packages. Right now it supports :ref:`the
+ASA capabilties in ARKODE <ARKODE.Mathematics.ASA>`, while the ASA capabilities in :ref:`CVODES
+<CVODES.Mathematics.ASA>` and :ref:`IDAS <IDAS.Mathematics.ASA>` must be used directly.
+
+.. toctree::
+   :maxdepth: 1
+
+   SUNAdjoint_links.rst

examples/arkode/C_petsc/CMakeLists.txt

@@ -62,6 +62,7 @@ foreach(example_tuple ${ARKODE_examples})
     ANSWER_FILE ${example}.out
     EXAMPLE_TYPE ${example_type})
 
+
   # libraries to link against
   target_link_libraries(${example} ${SUNDIALS_LIBS})