From cd0d0623859cf89fc0238da1f03843c6a986ae3a Mon Sep 17 00:00:00 2001 From: David Gardner Date: Mon, 17 Jul 2023 11:15:56 -0700 Subject: [PATCH 1/2] Relaxation Runge-Kutta Methods (#298) Add support for explicit, implicit, and IMEX relaxation Runge-Kutta methods in ERKStep and ARKStep --------- Co-authored-by: Daniel R. Reynolds Co-authored-by: Cody Balos --- CHANGELOG.md | 3 + doc/arkode/guide/source/Constants.rst | 21 + doc/arkode/guide/source/Introduction.rst | 5 + doc/arkode/guide/source/Mathematics.rst | 65 ++ .../Usage/ARKStep_c_interface/Relaxation.rst | 346 +++++++ .../Usage/ARKStep_c_interface/index.rst | 1 + .../Usage/ERKStep_c_interface/Relaxation.rst | 335 +++++++ .../Usage/ERKStep_c_interface/index.rst | 1 + .../guide/source/Usage/User_supplied.rst | 49 + doc/shared/sundials.bib | 91 ++ examples/arkode/CXX_serial/CMakeLists.txt | 7 + examples/arkode/CXX_serial/ark_pendulum.cpp | 433 +++++++++ examples/arkode/CXX_serial/ark_pendulum.out | 35 + .../C_parallel/ark_diurnal_kry_bbd_p.out | 4 +- .../arkode/C_parallel/ark_diurnal_kry_p.out | 2 +- .../arkode/C_parhyp/ark_diurnal_kry_ph.out | 2 +- examples/arkode/C_serial/CMakeLists.txt | 42 +- .../arkode/C_serial/ark_KrylovDemo_prec.out | 8 +- .../arkode/C_serial/ark_KrylovDemo_prec_1.out | 8 +- .../arkode/C_serial/ark_KrylovDemo_prec_2.out | 8 +- .../C_serial/ark_conserved_exp_entropy_ark.c | 528 +++++++++++ .../ark_conserved_exp_entropy_ark_1_0.out | 25 + .../ark_conserved_exp_entropy_ark_1_1.out | 53 ++ .../C_serial/ark_conserved_exp_entropy_erk.c | 436 +++++++++ .../ark_conserved_exp_entropy_erk_1.out | 25 + .../C_serial/ark_dissipated_exp_entropy.c | 478 ++++++++++ .../ark_dissipated_exp_entropy_1_0.out | 26 + .../ark_dissipated_exp_entropy_1_1.out | 59 ++ examples/templates/cmakelists_serial_C_ex.in | 3 + examples/utilities/example_utilities.hpp | 28 +- examples/utilities/plot_data_time_series.py | 208 +++++ include/arkode/arkode.h | 19 + include/arkode/arkode_arkstep.h | 30 + include/arkode/arkode_erkstep.h | 29 + scripts/arkode | 13 + src/arkode/CMakeLists.txt | 1 + src/arkode/arkode.c | 43 +- src/arkode/arkode_arkstep.c | 298 ++++++ src/arkode/arkode_arkstep_impl.h | 6 + src/arkode/arkode_arkstep_io.c | 80 +- src/arkode/arkode_erkstep.c | 158 ++++ src/arkode/arkode_erkstep_impl.h | 6 + src/arkode/arkode_erkstep_io.c | 80 ++ src/arkode/arkode_impl.h | 27 +- src/arkode/arkode_io.c | 8 + src/arkode/arkode_relaxation.c | 850 ++++++++++++++++++ src/arkode/arkode_relaxation_impl.h | 140 +++ src/arkode/arkode_root.c | 2 - src/arkode/arkode_types_impl.h | 23 + src/arkode/fmod/farkode_arkstep_mod.c | 228 +++++ src/arkode/fmod/farkode_arkstep_mod.f90 | 426 +++++++++ src/arkode/fmod/farkode_erkstep_mod.c | 228 +++++ src/arkode/fmod/farkode_erkstep_mod.f90 | 424 +++++++++ src/arkode/fmod/farkode_mod.f90 | 22 +- src/arkode/fmod/farkode_mristep_mod.c | 14 + src/arkode/fmod/farkode_mristep_mod.f90 | 26 + src/cvode/fmod/fcvode_mod.c | 14 + src/cvode/fmod/fcvode_mod.f90 | 26 + src/cvodes/fmod/fcvodes_mod.c | 14 + src/cvodes/fmod/fcvodes_mod.f90 | 26 + test/answers | 2 +- 61 files changed, 6538 insertions(+), 60 deletions(-) create mode 100644 doc/arkode/guide/source/Usage/ARKStep_c_interface/Relaxation.rst create mode 100644 doc/arkode/guide/source/Usage/ERKStep_c_interface/Relaxation.rst create mode 100644 examples/arkode/CXX_serial/ark_pendulum.cpp create mode 100644 examples/arkode/CXX_serial/ark_pendulum.out create mode 100644 examples/arkode/C_serial/ark_conserved_exp_entropy_ark.c create mode 100644 examples/arkode/C_serial/ark_conserved_exp_entropy_ark_1_0.out create mode 100644 examples/arkode/C_serial/ark_conserved_exp_entropy_ark_1_1.out create mode 100644 examples/arkode/C_serial/ark_conserved_exp_entropy_erk.c create mode 100644 examples/arkode/C_serial/ark_conserved_exp_entropy_erk_1.out create mode 100644 examples/arkode/C_serial/ark_dissipated_exp_entropy.c create mode 100644 examples/arkode/C_serial/ark_dissipated_exp_entropy_1_0.out create mode 100644 examples/arkode/C_serial/ark_dissipated_exp_entropy_1_1.out create mode 100755 examples/utilities/plot_data_time_series.py create mode 100644 src/arkode/arkode_relaxation.c create mode 100644 src/arkode/arkode_relaxation_impl.h create mode 100644 src/arkode/arkode_types_impl.h diff --git a/CHANGELOG.md b/CHANGELOG.md index 3d7f302794..e5d4be0291 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,9 @@ ## Changes to SUNDIALS in release 6.6.0 +Added support for relaxation Runge-Kutta methods to ERKStep and ARKStep in +ARKODE. + Added the second order IMEX method from Giraldo, Kelly, and Constantinescu 2013 as the default second order IMEX method in ARKStep. The explicit table is given by `ARKODE_ARK2_ERK_3_1_2` and the implicit table by `ARKODE_ARK2_DIRK_3_1_2`. diff --git a/doc/arkode/guide/source/Constants.rst b/doc/arkode/guide/source/Constants.rst index 8558965ee4..fd4e8f6777 100644 --- a/doc/arkode/guide/source/Constants.rst +++ b/doc/arkode/guide/source/Constants.rst @@ -62,6 +62,16 @@ contains the ARKODE output constants. +---------------------------------------------+-----------------------------------------------------------+ | | +---------------------------------------------+-----------------------------------------------------------+ + | **Relaxtion module input constants** | + +---------------------------------------------+-----------------------------------------------------------+ + | :index:`ARK_RELAX_BRENT` | Specifies Brent's method as the relaxation nonlinear | + | | solver. | + +---------------------------------------------+-----------------------------------------------------------+ + | :index:`ARK_RELAX_NEWTON` | Specifies Newton's method as the relaxation nonlinear | + | | solver. | + +---------------------------------------------+-----------------------------------------------------------+ + | | + +---------------------------------------------+-----------------------------------------------------------+ | **Explicit Butcher table specification** | +---------------------------------------------+-----------------------------------------------------------+ | :index:`ARKODE_HEUN_EULER_2_1_2` | Use the Heun-Euler-2-1-2 ERK method. | @@ -377,6 +387,17 @@ contains the ARKODE output constants. +-------------------------------------+------+------------------------------------------------------------+ | :index:`ARK_INVALID_TABLE` | -41 | An invalid Butcher or MRI table was encountered. | +-------------------------------------+------+------------------------------------------------------------+ + | :index:`ARK_CONTEXT_ERR` | -42 | An error occurred with the SUNDIALS context object | + +-------------------------------------+------+------------------------------------------------------------+ + | :index:`ARK_RELAX_FAIL` | -43 | An error occurred in computing the relaxation parameter | + +-------------------------------------+------+------------------------------------------------------------+ + | :index:`ARK_RELAX_MEM_FAIL` | -44 | The relaxation memory structure is ``NULL`` | + +-------------------------------------+------+------------------------------------------------------------+ + | :index:`ARK_RELAX_FUNC_FAIL` | -45 | The relaxation function returned an unrecoverable error | + +-------------------------------------+------+------------------------------------------------------------+ + | :index:`ARK_RELAX_JAC_FAIL` | -46 | The relaxation Jacobian function returned an unrecoverable | + | | | error | + +-------------------------------------+------+------------------------------------------------------------+ | :index:`ARK_UNRECOGNIZED_ERROR` | -99 | An unknown error was encountered. | +-------------------------------------+------+------------------------------------------------------------+ | | diff --git a/doc/arkode/guide/source/Introduction.rst b/doc/arkode/guide/source/Introduction.rst index 4703edde37..accfa96d36 100644 --- a/doc/arkode/guide/source/Introduction.rst +++ b/doc/arkode/guide/source/Introduction.rst @@ -121,6 +121,11 @@ Changes from previous versions Changes in v5.6.0 ----------------- + +Added support for relaxation Runge-Kutta methods in ERKStep and ARKStep, see +:numref:`ARKODE.Mathematics.Relaxation`, :numref:`ARKODE.Usage.ERKStep.Relaxation`, +and :numref:`ARKODE.Usage.ARKStep.Relaxation` for more information. + Added the second order IMEX method from :cite:p:`giraldo2013implicit` as the default second order IMEX method in ARKStep. The explicit table is given by ``ARKODE_ARK2_ERK_3_1_2`` (see :numref:`Butcher.ARK2_ERK`) and the implicit diff --git a/doc/arkode/guide/source/Mathematics.rst b/doc/arkode/guide/source/Mathematics.rst index 38c36b0c24..fbe18f3e40 100644 --- a/doc/arkode/guide/source/Mathematics.rst +++ b/doc/arkode/guide/source/Mathematics.rst @@ -2074,3 +2074,68 @@ step attempt, or fails with the minimum step size, then the integration is halte and an error is returned. In this case the user may need to employ other strategies as discussed in :numref:`ARKODE.Usage.ARKStep.Tolerances` and :numref:`ARKODE.Usage.ERKStep.Tolerances` to satisfy the inequality constraints. + +.. _ARKODE.Mathematics.Relaxation: + +Relaxation Methods +================== + +When the solution of :eq:`ARKODE_IVP` is conservative or dissipative with +respect to a smooth *convex* function :math:`\xi(y(t))`, it is desirable to have +the numerical method preserve these properties. That is +:math:`\xi(y_n) = \xi(y_{n-1}) = \ldots = \xi(y_{0})` for conservative systems +and :math:`\xi(y_n) \leq \xi(y_{n-1})` for dissipative systems. For examples +of such problems, see the references below and the citations there in. + +For such problems, ARKODE supports relaxation methods +:cite:p:`ketcheson2019relaxation, kang2022entropy, ranocha2020relaxation, ranocha2020hamiltonian` +applied to ERK, DIRK, or ARK methods to ensure dissipation or preservation of +the global function. The relaxed solution is given by + +.. math:: + y_r = y_{n-1} + r d = r y_n + (1 - r) y_{n - 1} + :label: ARKODE_RELAX_SOL + +where :math:`d` is the update to :math:`y_n` (i.e., +:math:`h_n \sum_{i=1}^{s}(b^E_i \hat{f}^E_i + b^I_i \hat{f}^I_i)` for ARKStep +and :math:`h_n \sum_{i=1}^{s} b_i f_i` for ERKStep) and :math:`r` is the +relaxation factor selected to ensure conservation or dissipation. Given an ERK, +DIRK, or ARK method of at least second order with non-negative solution weights +(i.e., :math:`b_i \geq 0` for ERKStep or :math:`b^E_i \geq 0` and +:math:`b^I_i \geq 0` for ARKStep), the factor :math:`r` is computed by solving +the auxiliary scalar nonlinear system + +.. math:: + F(r) = \xi(y_{n-1} + r d) - \xi(y_{n-1}) - r e = 0 + :label: ARKODE_RELAX_NLS + +at the end of each time step. The estimated change in :math:`\xi` is given by +:math:`e \equiv h_n \sum_{i=1}^{s} \langle \xi'(z_i), b^E_i f^E_i + b^I_i f^I_i \rangle` +where :math:`\xi'` is the Jacobian of :math:`\xi`. + +Two iterative methods are provided for solving :eq:`ARKODE_RELAX_NLS`, Newton's +method and Brent's method. When using Newton's method (the default), the +iteration is halted either when the residual tolerance is met, +:math:`F(r^{(k)}) < \epsilon_{\mathrm{relax\_res}}`, or when the difference +between successive iterates satisfies the relative and absolute tolerances, +:math:`|\delta_r^{(k)}| = |r^{(k)} - r^{(k-1)}| < \epsilon_{\mathrm{relax\_rtol}} |r^{(k-1)}| + \epsilon_{\mathrm{relax\_atol}}`. +Brent's method applies the same residual tolerance check and additionally halts +when the bisection update satisfies the relative and absolute tolerances, +:math:`|0.5 (r_c - r^{k})| < \epsilon_{\mathrm{relax\_rtol}} |r^{(k)}| + 0.5 \epsilon_{\mathrm{relax\_atol}}` +where :math:`r_c` and :math:`r^{(k)}` bound the root. + +If the nonlinear solve fails to meet the specified tolerances within the maximum +allowed number of iterations, the step size is reduced by the factor +:math:`\eta_\mathrm{rf}` (default 0.25) and the step is repeated. Additionally, +the solution of :eq:`ARKODE_RELAX_NLS` should be +:math:`r = 1 + \mathcal{O}(h_n^{q - 1})` for a method of order :math:`q` +:cite:p:`ranocha2020relaxation`. As such, limits are imposed on the range of +relaxation values allowed (i.e., limiting the maximum change in step size due to +relaxation). A relaxation value greater than :math:`r_\text{max}` (default 1.2) +or less than :math:`r_\text{min}` (default 0.8), is considered as a failed +relaxation application and the step will is repeated with the step size reduced +by :math:`\eta_\text{rf}`. + +For more information on utilizing relaxation Runge--Kutta methods, see +:numref:`ARKODE.Usage.ERKStep.Relaxation` and +:numref:`ARKODE.Usage.ARKStep.Relaxation`. diff --git a/doc/arkode/guide/source/Usage/ARKStep_c_interface/Relaxation.rst b/doc/arkode/guide/source/Usage/ARKStep_c_interface/Relaxation.rst new file mode 100644 index 0000000000..6d2027265d --- /dev/null +++ b/doc/arkode/guide/source/Usage/ARKStep_c_interface/Relaxation.rst @@ -0,0 +1,346 @@ +.. ----------------------------------------------------------------------------- + Programmer(s): David J. Gardner @ LLNL + ----------------------------------------------------------------------------- + SUNDIALS Copyright Start + Copyright (c) 2002-2022, 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 + ----------------------------------------------------------------------------- + +.. _ARKODE.Usage.ARKStep.Relaxation: + +Relaxation Methods +================== + +This section describes user-callable functions for applying relaxation methods +with ARKStep. For more information on relaxation Runge--Kutta methods see +:numref:`ARKODE.Mathematics.Relaxation`. + +.. note:: + + Relaxation support as not been evaluated with non-identity mass matrices. + While this usage mode is supported, feedback from users who explore this + combination would be appreciated. + +Enabling or Disabling Relaxation +-------------------------------- + +.. c:function:: int ARKStepSetRelaxFn(void* arkode_mem, ARKRelaxFn rfn, ARKRelaxJacFn rjac) + + Attaches the user supplied functions for evaluating the relaxation function + (``rfn``) and its Jacobian (``rjac``). + + Both ``rfn`` and ``rjac`` are required and an error will be returned if only + one of the functions is ``NULL``. If both ``rfn`` and ``rjac`` are ``NULL``, + relaxation is disabled. + + With DIRK and IMEX-ARK methods or when a fixed mass matrix is present, + applying relaxation requires allocating :math:`s` additional state vectors + (where :math:`s` is the number of stages in the method). + + :param arkode_mem: the ARKStep memory structure + :param rfn: the user-defined function to compute the relaxation function + :math:`\xi(y)` + :param rjac: the user-defined function to compute the relaxation Jacobian + :math:`\xi'(y)` + + :retval ARK_SUCCESS: the function exited successfully + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_ILL_INPUT: an invalid input combination was provided (see the + output error message for more details) + :retval ARK_MEM_FAIL: a memory allocation failed + + .. warning:: + + Applying relaxation requires using a method of at least second order with + :math:`b^E_i \geq 0` and :math:`b^I_i \geq 0`. If these conditions are not + satisfied, :c:func:`ARKStepEvolve` will return with an error during + initialization. + + .. note:: + + When combined with fixed time step sizes, ARKStep will attempt each step + using the specified step size. If the step is successful, relaxation will + be applied, effectively modifying the step size for the current step. If + the step fails or applying relaxation fails, :c:func:`ARKStepEvolve` will + return with an error. + + .. versionadded:: 5.6.0 + +Optional Input Functions +------------------------ + +This section describes optional input functions used to control applying +relaxation. + +.. c:function:: int ARKStepSetRelaxEtaFail(void* arkode_mem, sunrealtype eta_rf) + + Sets the step size reduction factor applied after a failed relaxation + application. + + The default value is 0.25. Input values :math:`\leq 0` or :math:`\geq 1` will + result in the default value being used. + + :param arkode_mem: the ARKStep memory structure + :param eta_rf: the step size reduction factor + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ARKStepSetRelaxLowerBound(void* arkode_mem, sunrealtype lower) + + Sets the smallest acceptable value for the relaxation parameter. + + Values smaller than the lower bound will result in a failed relaxation + application and the step will be repeated with a smaller step size + (determined by :c:func:`ARKStepSetRelaxEtaFail`). + + The default value is 0.8. Input values :math:`\leq 0` or :math:`\geq 1` will + result in the default value being used. + + :param arkode_mem: the ARKStep memory structure + :param lower: the relaxation parameter lower bound + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ARKStepSetRelaxUpperBound(void* arkode_mem, sunrealtype upper) + + Sets the largest acceptable value for the relaxation parameter. + + Values larger than the upper bound will result in a failed relaxation + application and the step will be repeated with a smaller step size + (determined by :c:func:`ARKStepSetRelaxEtaFail`). + + The default value is 1.2. Input values :math:`\leq 1` will result in the + default value being used. + + :param arkode_mem: the ARKStep memory structure + :param upper: the relaxation parameter upper bound + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ARKStepSetRelaxMaxFails(void* arkode_mem, int max_fails) + + Sets the maximum number of times applying relaxation can fail within a step + attempt before the integration is halted with an error. + + The default value is 10. Input values :math:`\leq 0` will result in the + default value being used. + + :param arkode_mem: the ARKStep memory structure + :param max_fails: the maximum number of failed relaxation applications + allowed in a step + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ARKStepSetRelaxMaxIters(void* arkode_mem, int max_iters) + + Sets the maximum number of nonlinear iterations allowed when solving for the + relaxation parameter. + + If the maximum number of iterations is reached before meeting the solve + tolerance (determined by :c:func:`ARKStepSetRelaxResTol` and + :c:func:`ARKStepSetRelaxTol`), the step will be repeated with a smaller + step size (determined by :c:func:`ARKStepSetRelaxEtaFail`). + + The default value is 10. Input values :math:`\leq 0` will result in the + default value being used. + + :param arkode_mem: the ARKStep memory structure + :param max_iters: the maximum number of solver iterations allowed + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ARKStepSetRelaxSolver(void* arkode_mem, ARKRelaxSolver solver) + + Sets the nonlinear solver method used to compute the relaxation parameter. + + The default value is ``ARK_RELAX_NEWTON``. + + :param arkode_mem: the ARKStep memory structure + :param solver: the nonlinear solver to use: ``ARK_RELAX_BRENT`` or + ``ARK_RELAX_NEWTON`` + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + :retval ARK_ILL_INPUT: an invalid solver option was provided + + .. versionadded:: 5.6.0 + +.. c:function:: int ARKStepSetRelaxResTol(void* arkode_mem, sunrealtype res_tol) + + Sets the nonlinear solver residual tolerance to use when solving + :eq:`ARKODE_RELAX_NLS`. + + If the residual or iteration update tolerance (see + :c:func:`ARKStepSetRelaxMaxIter`) is not reached within the maximum number of + iterations (determined by :c:func:`ARKStepSetRelaxMaxIters`), the step will + be repeated with a smaller step size (determined by + :c:func:`ARKStepSetRelaxEtaFail`). + + The default value is :math:`4 \epsilon` where :math:`\epsilon` is + floating-point precision. Input values :math:`\leq 0` will result in the + default value being used. + + :param arkode_mem: the ARKStep memory structure + :param res_tol: the nonlinear solver residual tolerance to use + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ARKStepSetRelaxTol(void* arkode_mem, sunrealtype rel_tol, sunrealtype abs_tol) + + Sets the nonlinear solver relative and absolute tolerance on changes in + :math:`r` iterates when solving :eq:`ARKODE_RELAX_NLS`. + + If the residual (see :c:func:`ARKStepSetRelaxResTol`) or iterate update + tolerance is not reached within the maximum number of iterations (determined + by :c:func:`ARKStepSetRelaxMaxIters`), the step will be repeated with a + smaller step size (determined by :c:func:`ARKStepSetRelaxEtaFail`). + + The default relative and absolute tolerances are :math:`4 \epsilon` and + :math:`10^{-14}`, respectively, where :math:`\epsilon` is floating-point + precision. Input values :math:`\leq 0` will result in the default value being + used. + + :param arkode_mem: the ARKStep memory structure + :param rel_tol: the nonlinear solver relative solution tolerance to use + :param abs_tol: the nonlinear solver absolute solution tolerance to use + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +Optional Output Functions +------------------------- + +This section describes optional output functions used to retrieve information +about the performance of the relaxation method. + +.. c:function:: int ARKStepGetNumRelaxFnEvals(void* arkode_mem, long int* r_evals) + + Get the number of times the user's relaxation function was evaluated. + + :param arkode_mem: the ARKStep memory structure + :param r_evals: the number of relaxation function evaluations + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ARKStepGetNumRelaxJacEvals(void* arkode_mem, long int* J_evals) + + Get the number of times the user's relaxation Jacobian was evaluated. + + :param arkode_mem: the ARKStep memory structure + :param J_evals: the number of relaxation Jacobian evaluations + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ARKStepGetNumRelaxFails(void* arkode_mem, long int* fails) + + Get the total number of times applying relaxation failed. + + The counter includes the sum of the number of nonlinear solver failures + (see :c:func:`ARKStepGetNumRelaxSolveFails`) and the number of failures due + an unacceptable relaxation value (see :c:func:`ARKStepSetRelaxBoundFactor`). + + :param arkode_mem: the ARKStep memory structure + :param fails: the total number of failed relaxation attempts + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + + +.. c:function:: int ARKStepGetNumRelaxBoundFails(void* arkode_mem, long int* fails) + + Get the number of times the relaxation parameter was deemed unacceptable. + + :param arkode_mem: the ARKStep memory structure + :param fails: the number of failures due to an unacceptable relaxation + parameter value + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ARKStepGetNumRelaxSolveFails(void* arkode_mem, long int* fails) + + Get the number of times the relaxation parameter nonlinear solver failed. + + :param arkode_mem: the ARKStep memory structure + :param fails: the number of relaxation nonlinear solver failures + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ARKStepGetNumRelaxSolveIters(void* arkode_mem, long int* iters) + + Get the number of relaxation parameter nonlinear solver iterations. + + :param arkode_mem: the ARKStep memory structure + :param iters: the number of relaxation nonlinear solver iterations + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 diff --git a/doc/arkode/guide/source/Usage/ARKStep_c_interface/index.rst b/doc/arkode/guide/source/Usage/ARKStep_c_interface/index.rst index 1d2c076a2c..2cd52cd7fd 100644 --- a/doc/arkode/guide/source/Usage/ARKStep_c_interface/index.rst +++ b/doc/arkode/guide/source/Usage/ARKStep_c_interface/index.rst @@ -59,5 +59,6 @@ detailed in the following subsections. Skeleton User_callable + Relaxation Preconditioners XBraid diff --git a/doc/arkode/guide/source/Usage/ERKStep_c_interface/Relaxation.rst b/doc/arkode/guide/source/Usage/ERKStep_c_interface/Relaxation.rst new file mode 100644 index 0000000000..ab45d7f54c --- /dev/null +++ b/doc/arkode/guide/source/Usage/ERKStep_c_interface/Relaxation.rst @@ -0,0 +1,335 @@ +.. ----------------------------------------------------------------------------- + Programmer(s): David J. Gardner @ LLNL + ----------------------------------------------------------------------------- + SUNDIALS Copyright Start + Copyright (c) 2002-2022, 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 + ----------------------------------------------------------------------------- + +.. _ARKODE.Usage.ERKStep.Relaxation: + +Relaxation Methods +================== + +This section describes user-callable functions for applying relaxation methods +with ERKStep. For more information on relaxation Runge--Kutta methods see +:numref:`ARKODE.Mathematics.Relaxation`. + +Enabling or Disabling Relaxation +-------------------------------- + +.. c:function:: int ERKStepSetRelaxFn(void* arkode_mem, ARKRelaxFn rfn, ARKRelaxJacFn rjac) + + Attaches the user supplied functions for evaluating the relaxation function + (``rfn``) and its Jacobian (``rjac``). + + Both ``rfn`` and ``rjac`` are required and an error will be returned if only + one of the functions is ``NULL``. If both ``rfn`` and ``rjac`` are ``NULL``, + relaxation is disabled. + + :param arkode_mem: the ERKStep memory structure + :param rfn: the user-defined function to compute the relaxation function + :math:`\xi(y)` + :param rjac: the user-defined function to compute the relaxation Jacobian + :math:`\xi'(y)` + + :retval ARK_SUCCESS: the function exited successfully + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_ILL_INPUT: an invalid input combination was provided (see the + output error message for more details) + :retval ARK_MEM_FAIL: a memory allocation failed + + .. warning:: + + Applying relaxation requires using a method of at least second order with + :math:`b_i \geq 0`. If these conditions are not satisfied, + :c:func:`ERKStepEvolve` will return with an error during initialization. + + .. note:: + + When combined with fixed time step sizes, ERKStep will attempt each step + using the specified step size. If the step is successful, relaxation will + be applied, effectively modifying the step size for the current step. If + the step fails or applying relaxation fails, :c:func:`ERKStepEvolve` will + return with an error. + + .. versionadded:: 5.6.0 + +Optional Input Functions +------------------------ + +This section describes optional input functions used to control applying +relaxation. + +.. c:function:: int ERKStepSetRelaxEtaFail(void* arkode_mem, sunrealtype eta_rf) + + Sets the step size reduction factor applied after a failed relaxation + application. + + The default value is 0.25. Input values :math:`\leq 0` or :math:`\geq 1` will + result in the default value being used. + + :param arkode_mem: the ERKStep memory structure + :param eta_rf: the step size reduction factor + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ERKStepSetRelaxLowerBound(void* arkode_mem, sunrealtype lower) + + Sets the smallest acceptable value for the relaxation parameter. + + Values smaller than the lower bound will result in a failed relaxation + application and the step will be repeated with a smaller step size + (determined by :c:func:`ERKStepSetRelaxEtaFail`). + + The default value is 0.8. Input values :math:`\leq 0` or :math:`\geq 1` will + result in the default value being used. + + :param arkode_mem: the ERKStep memory structure + :param lower: the relaxation parameter lower bound + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ERKStepSetRelaxUpperBound(void* arkode_mem, sunrealtype upper) + + Sets the largest acceptable value for the relaxation parameter. + + Values larger than the upper bound will result in a failed relaxation + application and the step will be repeated with a smaller step size + (determined by :c:func:`ERKStepSetRelaxEtaFail`). + + The default value is 1.2. Input values :math:`\leq 1` will result in the + default value being used. + + :param arkode_mem: the ERKStep memory structure + :param upper: the relaxation parameter upper bound + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ERKStepSetRelaxMaxFails(void* arkode_mem, int max_fails) + + Sets the maximum number of times applying relaxation can fail within a step + attempt before the integration is halted with an error. + + The default value is 10. Input values :math:`\leq 0` will result in the + default value being used. + + :param arkode_mem: the ERKStep memory structure + :param max_fails: the maximum number of failed relaxation applications + allowed in a step + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ERKStepSetRelaxMaxIters(void* arkode_mem, int max_iters) + + Sets the maximum number of nonlinear iterations allowed when solving for the + relaxation parameter. + + If the maximum number of iterations is reached before meeting the solve + tolerance (determined by :c:func:`ERKStepSetRelaxResTol` and + :c:func:`ERKStepSetRelaxTol`), the step will be repeated with a smaller + step size (determined by :c:func:`ERKStepSetRelaxEtaFail`). + + The default value is 10. Input values :math:`\leq 0` will result in the + default value being used. + + :param arkode_mem: the ERKStep memory structure + :param max_iters: the maximum number of solver iterations allowed + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ERKStepSetRelaxSolver(void* arkode_mem, ARKRelaxSolver solver) + + Sets the nonlinear solver method used to compute the relaxation parameter. + + The default value is ``ARK_RELAX_NEWTON``. + + :param arkode_mem: the ERKStep memory structure + :param solver: the nonlinear solver to use: ``ARK_RELAX_BRENT`` or + ``ARK_RELAX_NEWTON`` + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + :retval ARK_ILL_INPUT: an invalid solver option was provided + + .. versionadded:: 5.6.0 + +.. c:function:: int ERKStepSetRelaxResTol(void* arkode_mem, sunrealtype res_tol) + + Sets the nonlinear solver residual tolerance to use when solving + :eq:`ARKODE_RELAX_NLS`. + + If the residual or solution tolerance (see :c:func:`ERKStepSetRelaxMaxIter`) + is not reached within the maximum number of iterations (determined by + :c:func:`ERKStepSetRelaxMaxIters`), the step will be repeated with a smaller + step size (determined by :c:func:`ERKStepSetRelaxEtaFail`). + + The default value is :math:`4 \epsilon` where :math:`\epsilon` is + floating-point precision. Input values :math:`\leq 0` will result in the + default value being used. + + :param arkode_mem: the ERKStep memory structure + :param res_tol: the nonlinear solver residual tolerance to use + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ERKStepSetRelaxTol(void* arkode_mem, sunrealtype rel_tol, sunrealtype abs_tol) + + Sets the nonlinear solver relative and absolute tolerance on changes in + :math:`r` when solving :eq:`ARKODE_RELAX_NLS`. + + + If the residual (see :c:func:`ERKStepSetRelaxResTol`) or solution tolerance + is not reached within the maximum number of iterations (determined by + :c:func:`ERKStepSetRelaxMaxIters`), the step will be repeated with a smaller + step size (determined by :c:func:`ERKStepSetRelaxEtaFail`). + + The default relative and absolute tolerances are :math:`4 \epsilon` and + :math:`10^{-14}`, respectively, where :math:`\epsilon` is floating-point + precision. Input values :math:`\leq 0` will result in the default value being + used. + + :param arkode_mem: the ERKStep memory structure + :param rel_tol: the nonlinear solver relative solution tolerance to use + :param abs_tol: the nonlinear solver absolute solution tolerance to use + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +Optional Output Functions +------------------------- + +This section describes optional output functions used to retrieve information +about the performance of the relaxation method. + +.. c:function:: int ERKStepGetNumRelaxFnEvals(void* arkode_mem, long int* r_evals) + + Get the number of times the user's relaxation function was evaluated. + + :param arkode_mem: the ERKStep memory structure + :param r_evals: the number of relaxation function evaluations + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ERKStepGetNumRelaxJacEvals(void* arkode_mem, long int* J_evals) + + Get the number of times the user's relaxation Jacobian was evaluated. + + :param arkode_mem: the ERKStep memory structure + :param J_evals: the number of relaxation Jacobian evaluations + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ERKStepGetNumRelaxFails(void* arkode_mem, long int* fails) + + Get the total number of times applying relaxation failed. + + The counter includes the sum of the number of nonlinear solver failures + (see :c:func:`ERKStepGetNumRelaxSolveFails`) and the number of failures due + an unacceptable relaxation value (see :c:func:`ERKStepSetRelaxBoundFactor`). + + :param arkode_mem: the ERKStep memory structure + :param fails: the total number of failed relaxation attempts + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + + +.. c:function:: int ERKStepGetNumRelaxBoundFails(void* arkode_mem, long int* fails) + + Get the number of times the relaxation parameter was deemed unacceptable. + + :param arkode_mem: the ERKStep memory structure + :param fails: the number of failures due to an unacceptable relaxation + parameter value + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ERKStepGetNumRelaxSolveFails(void* arkode_mem, long int* fails) + + Get the number of times the relaxation parameter nonlinear solver failed. + + :param arkode_mem: the ERKStep memory structure + :param fails: the number of relaxation nonlinear solver failures + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 + +.. c:function:: int ERKStepGetNumRelaxSolveIters(void* arkode_mem, long int* iters) + + Get the number of relaxation parameter nonlinear solver iterations. + + :param arkode_mem: the ERKStep memory structure + :param iters: the number of relaxation nonlinear solver iterations + + :retval ARK_SUCCESS: the value was successfully set + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL`` + :retval ARK_RELAX_MEM_NULL: the internal relaxation memory structure was + ``NULL`` + + .. versionadded:: 5.6.0 diff --git a/doc/arkode/guide/source/Usage/ERKStep_c_interface/index.rst b/doc/arkode/guide/source/Usage/ERKStep_c_interface/index.rst index 760c38e58d..184ea401e2 100644 --- a/doc/arkode/guide/source/Usage/ERKStep_c_interface/index.rst +++ b/doc/arkode/guide/source/Usage/ERKStep_c_interface/index.rst @@ -42,3 +42,4 @@ detailed in the following sub-sections. Skeleton User_callable + Relaxation diff --git a/doc/arkode/guide/source/Usage/User_supplied.rst b/doc/arkode/guide/source/Usage/User_supplied.rst index 21ed62475e..9a5c83be79 100644 --- a/doc/arkode/guide/source/Usage/User_supplied.rst +++ b/doc/arkode/guide/source/Usage/User_supplied.rst @@ -63,6 +63,9 @@ The user-supplied functions for ARKODE consist of: by the outer integrator to the inner integrator, or state data supplied by the inner integrator to the outer integrator. +* if relaxation is enabled (optional), a function that evaluates the + conservative or dissipative function :math:`\xi(y(t))` (required) and a + function to evaluate its Jacobian :math:`\xi'(y(t))` (required). .. _ARKODE.Usage.ODERHS: @@ -1114,3 +1117,49 @@ outer integrator for the outer integration. **Notes:** In a heterogeneous computing environment if any data copies between the host and device vector data are necessary, this is where that should occur. + + +.. _ARKODE.Usage.RelaxFn: + +Relaxation function +------------------- + +.. c:type:: int (*ARKRelaxFn)(N_Vector y, realtype* r, void* user_data) + + When applying relaxation, an :c:func:`ARKRelaxFn` function is required to + compute the conservative or dissipative function :math:`\xi(y)`. + + **Arguments:** + * *y* -- the current value of the dependent variable vector. + * *r* -- the value of :math:`\xi(y)`. + * *user_data* -- the ``user_data`` pointer that was passed to + :c:func:`ARKStepSetUserData`. + + **Return value:** + An :c:func:`ARKRelaxFn` function should return 0 if successful, a positive + value if a recoverable error occurred, or a negative value if an + unrecoverable error occurred. If a recoverable error occurs, the step size + will be reduced and the step repeated. + +.. _ARKODE.Usage.RelaxJacFn: + +Relaxation Jacobian function +---------------------------- + +.. c:type:: int (*ARKRelaxJacFn)(N_Vector y, N_Vector J, void* user_data); + + When applying relaxation, an :c:func:`ARKRelaxJacFn` function is required to + compute the Jacobian :math:`\xi'(y)` of the :c:func:`ARKRelaxFn` + :math:`\xi(y)`. + + **Arguments:** + * *y* -- the current value of the dependent variable vector. + * *J* -- the Jacobian vector :math:`\xi'(y)`. + * *user_data* -- the ``user_data`` pointer that was passed to + :c:func:`ARKStepSetUserData`. + + **Return value:** + An :c:func:`ARKRelaxJacFn` function should return 0 if successful, a + positive value if a recoverable error occurred, or a negative value if an + unrecoverable error occurred. If a recoverable error occurs, the step size + will be reduced and the step repeated. diff --git a/doc/shared/sundials.bib b/doc/shared/sundials.bib index d62f0d9c22..9577ea3e62 100644 --- a/doc/shared/sundials.bib +++ b/doc/shared/sundials.bib @@ -2002,6 +2002,97 @@ @techreport{Zon:63 year = {1963} } +% +% Relaxation Runge-Kutta Methods +% + +% IMEX and Multirate Relaxation Methods +@article{kang2022entropy, + title = {{Entropy--Preserving and Entropy--Stable Relaxation IMEX and Multirate Time--Stepping Methods}}, + author = {Kang, Shinhoo and Constantinescu, Emil M}, + journal = {Journal of Scientific Computing}, + volume = {93}, + number = {1}, + pages = {1--31}, + year = {2022}, + publisher = {Springer}, + doi = {10.1007/s10915-022-01982-w} +} + +@article{ketcheson2019relaxation, + title = {{Relaxation Runge--Kutta methods: Conservation and stability for inner-product norms}}, + author = {Ketcheson, David I}, + journal = {SIAM Journal on Numerical Analysis}, + volume = {57}, + number = {6}, + pages = {2850--2870}, + year = {2019}, + publisher = {SIAM}, + doi = {10.1137/19M1263662} +} + +% Local Relaxation Runge-Kutta +@article{ranocha2020fully, + title = {{Fully discrete explicit locally entropy-stable schemes for the compressible Euler and Navier--Stokes equations}}, + author = {Ranocha, Hendrik and Dalcin, Lisandro and Parsani, Matteo}, + journal = {Computers \& Mathematics with Applications}, + volume = {80}, + number = {5}, + pages = {1343--1359}, + year = {2020}, + publisher = {Elsevier}, + doi = {10.1016/j.camwa.2020.06.016} +} + +% Relaxation in Multistep and General Linear Methods +@article{ranocha2020general, + title = {General relaxation methods for initial-value problems with application to multistep schemes}, + author = {Ranocha, Hendrik and L{\'o}czi, Lajos and Ketcheson, David I}, + journal = {Numerische Mathematik}, + volume = {146}, + number = {4}, + pages = {875--906}, + year = {2020}, + publisher = {Springer}, + doi = {10.1007/s00211-020-01158-4} +} + +% Global Relaxation Runge-Kutta +@article{ranocha2020relaxation, + title = {{Relaxation Runge--Kutta Methods: Fully Discrete Explicit Entropy-Stable Schemes for the Compressible Euler and Navier--Stokes Equations}}, + author = {Ranocha, Hendrik and Sayyari, Mohammed and Dalcin, Lisandro and Parsani, Matteo and Ketcheson, David I}, + journal = {SIAM Journal on Scientific Computing}, + volume = {42}, + number = {2}, + pages = {A612--A638}, + year = {2020}, + publisher = {SIAM}, + doi = {10.1137/19M1263480} +} + +% Relaxation Hamiltontian +@article{ranocha2020hamiltonian, + title = {{Relaxation Runge--Kutta methods for Hamiltonian problems}}, + author = {Ranocha, Hendrik and Ketcheson, David I}, + journal = {Journal of Scientific Computing}, + volume = {84}, + number = {1}, + pages = {1--27}, + year = {2020}, + publisher = {Springer}, + doi = {10.1007/s10915-020-01277-y} +} + +@article{rogowski2022performance, + title = {{Performance analysis of relaxation Runge--Kutta methods}}, + author = {Rogowski, Marcin and Dalcin, Lisandro and Parsani, Matteo and Keyes, David E}, + journal = {The International Journal of High Performance Computing Applications}, + pages = {10943420221085947}, + year = {2022}, + publisher = {SAGE Publications Sage UK: London, England}, + doi = {10.1177/10943420221085947} +} + % % Kokkos % diff --git a/examples/arkode/CXX_serial/CMakeLists.txt b/examples/arkode/CXX_serial/CMakeLists.txt index f92b4008e4..bded241a5f 100644 --- a/examples/arkode/CXX_serial/CMakeLists.txt +++ b/examples/arkode/CXX_serial/CMakeLists.txt @@ -34,6 +34,7 @@ set(ARKODE_examples "ark_kpr_Mt.cpp\;2 4 0 -10 0\;develop" "ark_kpr_Mt.cpp\;0 4 0 -10 1 10 1\;develop" "ark_kpr_Mt.cpp\;0 4 0 -10 0 10 1\;develop" + "ark_pendulum.cpp\;\;develop" ) # Header files to install @@ -70,6 +71,11 @@ foreach(example_tuple ${ARKODE_examples}) target_compile_features(${example_target} PRIVATE cxx_std_14) + # directories to include + target_include_directories(${example_target} + PRIVATE + ${PROJECT_SOURCE_DIR}/examples/utilities) + # libraries to link against target_link_libraries(${example_target} sundials_arkode @@ -122,6 +128,7 @@ if(EXAMPLES_INSTALL) README ${ARKODE_extras} ${ARKODE_headers} + "${PROJECT_SOURCE_DIR}/examples/utilities/example_utilities.hpp" TEST_INSTALL CXX_serial ) diff --git a/examples/arkode/CXX_serial/ark_pendulum.cpp b/examples/arkode/CXX_serial/ark_pendulum.cpp new file mode 100644 index 0000000000..d56ee752b8 --- /dev/null +++ b/examples/arkode/CXX_serial/ark_pendulum.cpp @@ -0,0 +1,433 @@ +/* ----------------------------------------------------------------------------- + * Programmer(s): David J. Gardner @ LLNL + * ----------------------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2002-2022, 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 + * ----------------------------------------------------------------------------- + * This example problem is adapted from: + * + * H. Ranocha, M. Sayyari, L. Dalcin, M. Parsani, and D.I. Ketcheson, + * "Relaxation Runge-Kutta Methods: Fully-Discrete Explicit Entropy-Stable + * Schemes for the Compressible Euler and Navier-Stokes Equations," SIAM Journal + * on Scientific Computing, 42(2), 2020, https://doi.org/10.1137/19M1263480. + * ----------------------------------------------------------------------------- + * This example evolves system y = [u, v]^T + * + * du/dt = -sin(v) + * dv/dt = u + * + * for t in the interval [0, 10] with the initial condition y0 = [1.5 1.0]^T. + * The conserved energy and its Jacobian for the system are + * + * e(y) = 0.5 u^2 - cos(v) and e'(y) = [u, sin(v)]^T + * + * The problem is advanced in time with an explicit or implicit relaxed + * Runge-Kutta method to ensure conservation of the energy. + * ---------------------------------------------------------------------------*/ + +// Standard headers +#include +#include +#include +#include +#include +#include + +// Common utility functions +#include + +// SUNDIALS headers +#include +#include +#include +#include + +#include "arkode/arkode_butcher.h" +#include "sundials/sundials_types.h" + +/* ----------------------- * + * User-supplied functions * + * ----------------------- */ + +// ODE RHS function +int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data); + +// ODE RHS Jacobian function +int Jac(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, + N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); + +// Energy function +int Eng(N_Vector y, sunrealtype* e, void* user_data); + +// Energy Jacobian function +int JacEng(N_Vector y, N_Vector J, void* user_data); + +/* ------------ * + * Main Program * + * ------------ */ + +int main(int argc, char* argv[]) +{ + // Create the SUNDIALS context object for this simulation + sundials::Context ctx; + + // Initial and final times + sunrealtype t0 = SUN_RCONST(0.0); + sunrealtype tf = SUN_RCONST(10.0); + + // Relative and absolute tolerances + sunrealtype reltol = SUN_RCONST(1.0e-6); + sunrealtype abstol = SUN_RCONST(1.0e-10); + + // Command line options + int relax = 1; // enable relaxation + int implicit = 1; // implicit + sunrealtype fixed_h = SUN_RCONST(0.0); // adaptive + + /* -------------------- * + * Output Problem Setup * + * -------------------- */ + + if (argc > 1) relax = atoi(argv[1]); + if (argc > 2) implicit = atoi(argv[2]); + if (argc > 3) fixed_h = atof(argv[3]); + + std::cout << "Nonlinear Pendulum problem:\n"; + if (implicit) { std::cout << " method = DIRK\n"; } + else { std::cout << " method = ERK\n"; } + if (relax) { std::cout << " relaxation = ON\n"; } + else { std::cout << " relaxation = OFF\n"; } + std::cout << " reltol = " << reltol << "\n"; + std::cout << " abstol = " << abstol << "\n"; + if (fixed_h > SUN_RCONST(0.0)) + { + std::cout << " fixed h = " << fixed_h << "\n"; + } + std::cout << std::endl; + + /* ------------ * + * Setup ARKODE * + * ------------ */ + + // Create serial vector and set the initial condition values + N_Vector y = N_VNew_Serial(2, ctx); + if (check_ptr(y, "N_VNew_Serial")) return 1; + + sunrealtype* ydata = N_VGetArrayPointer(y); + if (check_ptr(ydata, "N_VGetArrayPointer")) return 1; + + ydata[0] = SUN_RCONST(1.5); + ydata[1] = SUN_RCONST(1.0); + + // Initial energy + sunrealtype eng0; + int flag = Eng(y, &eng0, nullptr); + if (check_flag(flag, "Eng")) return 1; + + // Initialize ARKStep + void* arkode_mem = nullptr; + if (implicit) { arkode_mem = ARKStepCreate(nullptr, f, t0, y, ctx); } + else { arkode_mem = ARKStepCreate(f, nullptr, t0, y, ctx); } + if (check_ptr(arkode_mem, "ARKStepCreate")) return 1; + + // Specify tolerances + flag = ARKStepSStolerances(arkode_mem, reltol, abstol); + if (check_flag(flag, "ARKStepSStolerances")) return 1; + + if (relax) + { + // Enable relaxation methods + flag = ARKStepSetRelaxFn(arkode_mem, Eng, JacEng); + if (check_flag(flag, "ARKStepSetRelaxFn")) return 1; + } + + SUNMatrix A = nullptr; + SUNLinearSolver LS = nullptr; + + if (implicit) + { + // Create dense matrix and linear solver + A = SUNDenseMatrix(2, 2, ctx); + if (check_ptr(A, "SUNDenseMatrix")) return 1; + + LS = SUNLinSol_Dense(y, A, ctx); + if (check_ptr(LS, "SUNLinSol_Dense")) return 1; + + // Attach the matrix and linear solver + flag = ARKStepSetLinearSolver(arkode_mem, LS, A); + if (check_flag(flag, "ARKStepSetLinearSolver")) return 1; + + // Set Jacobian routine + flag = ARKStepSetJacFn(arkode_mem, Jac); + if (check_flag(flag, "ARKStepSetJacFn")) return 1; + + if (fixed_h > SUN_RCONST(0.0)) + { + // 3rd-order SDIRK method of Norsett + ARKodeButcherTable B = ARKodeButcherTable_Alloc(2, SUNFALSE); + + const sunrealtype gamma = ((SUN_RCONST(3.0) + std::sqrt(SUN_RCONST(3.0))) + / SUN_RCONST(6.0)); + + B->A[0][0] = gamma; + B->A[1][0] = SUN_RCONST(1.0) - SUN_RCONST(2.0) * gamma; + B->A[1][1] = gamma; + + B->c[0] = gamma; + B->c[1] = SUN_RCONST(1.0) - gamma; + + B->b[0] = SUN_RCONST(0.5); + B->b[1] = SUN_RCONST(0.5); + + B->q = 3; + B->p = 0; + + flag = ARKStepSetTables(arkode_mem, 3, 0, B, nullptr); + if (check_flag(flag, "ARKStepSetTables")) return 1; + + ARKodeButcherTable_Free(B); + } + else + { + // Select a Butcher table with non-negative b values + flag = ARKStepSetTableName(arkode_mem, "ARKODE_SDIRK_2_1_2", + "ARKODE_ERK_NONE"); + if (check_flag(flag, "ARKStepSetTableName")) return 1; + } + } + + if (fixed_h > SUN_RCONST(0.0)) + { + flag = ARKStepSetFixedStep(arkode_mem, fixed_h); + if (check_flag(flag, "ARKStepSetFixedStep")) return 1; + } + + /* --------------- * + * Advance in Time * + * --------------- */ + + // Initial time + sunrealtype t = t0; + + // Output the initial condition and energy + int swidth = 8; + int rwidth = std::numeric_limits::digits10 + 8; + + std::ofstream outfile("ark_pendulum.txt"); + outfile << "# vars: t u v energy energy_err\n"; + outfile << std::scientific; + outfile << std::setprecision(std::numeric_limits::digits10); + outfile << t << " " << ydata[0] << " " << ydata[1] << " " << eng0 << " " + << SUN_RCONST(0.0) << std::endl; + + std::cout << std::setw(swidth) << "step" << std::setw(rwidth) << "t" + << std::setw(rwidth) << "u" << std::setw(rwidth) << "v" + << std::setw(rwidth) << "e" << std::setw(rwidth) << "e err" + << std::endl; + for (int i = 0; i < swidth + 5 * rwidth; i++) std::cout << "-"; + + std::cout << std::endl; + std::cout << std::scientific; + std::cout << std::setprecision(std::numeric_limits::digits10); + std::cout << std::setw(swidth) << 0 << std::setw(rwidth) << t + << std::setw(rwidth) << ydata[0] << std::setw(rwidth) << ydata[1] + << std::setw(rwidth) << eng0 << std::setw(rwidth) + << SUN_RCONST(0.0); + std::cout << std::endl; + + while (t < tf) + { + // Evolve in time + flag = ARKStepEvolve(arkode_mem, tf, y, &t, ARK_ONE_STEP); + if (check_flag(flag, "ARKStepEvolve")) break; + + // Output solution and errors + sunrealtype eng; + flag = Eng(y, &eng, nullptr); + if (check_flag(flag, "Eng")) return 1; + + sunrealtype eng_chng = eng - eng0; + + /* Output to the screen periodically */ + long int nst; + flag = ARKStepGetNumSteps(arkode_mem, &nst); + check_flag(flag, "ARKStepGetNumSteps"); + + if (nst % 1000 == 0) + { + std::cout << std::setw(swidth) << nst << std::setw(rwidth) << t + << std::setw(rwidth) << ydata[0] << std::setw(rwidth) + << ydata[1] << std::setw(rwidth) << eng << std::setw(rwidth) + << eng_chng << std::endl; + } + + /* Write all steps to file */ + outfile << t << " " << ydata[0] << " " << ydata[1] << " " << eng << " " + << eng_chng << std::endl; + } + + for (int i = 0; i < swidth + 5 * rwidth; i++) std::cout << "-"; + std::cout << std::endl; + outfile.close(); + + /* ------------ * + * Output Stats * + * ------------ */ + + // ARKODE statistics + long int nst, nst_a, netf, nfe, nfi; + + // Get final statistics on how the solve progressed + flag = ARKStepGetNumSteps(arkode_mem, &nst); + check_flag(flag, "ARKStepGetNumSteps"); + + flag = ARKStepGetNumStepAttempts(arkode_mem, &nst_a); + check_flag(flag, "ARKStepGetNumStepAttempts"); + + flag = ARKStepGetNumErrTestFails(arkode_mem, &netf); + check_flag(flag, "ARKStepGetNumErrTestFails"); + + flag = ARKStepGetNumRhsEvals(arkode_mem, &nfe, &nfi); + check_flag(flag, "ARKStepGetNumRhsEvals"); + + std::cout << std::endl; + std::cout << "Final Solver Statistics:\n"; + std::cout << " Internal solver steps = " << nst << " (attempted = " << nst_a + << ")\n"; + std::cout << " Total number of error test failures = " << netf << "\n"; + std::cout << " Total RHS evals: Fe = " << nfe << ", Fi = " << nfi << "\n"; + + if (implicit) + { + long int nsetups, nje, nfeLS, nni, ncfn; + + flag = ARKStepGetNumNonlinSolvIters(arkode_mem, &nni); + check_flag(flag, "ARKStepGetNumNonlinSolvIters"); + + flag = ARKStepGetNumNonlinSolvConvFails(arkode_mem, &ncfn); + check_flag(flag, "ARKStepGetNumNonlinSolvConvFails"); + + flag = ARKStepGetNumLinSolvSetups(arkode_mem, &nsetups); + check_flag(flag, "ARKStepGetNumLinSolvSetups"); + + flag = ARKStepGetNumJacEvals(arkode_mem, &nje); + check_flag(flag, "ARKStepGetNumJacEvals"); + + flag = ARKStepGetNumLinRhsEvals(arkode_mem, &nfeLS); + check_flag(flag, "ARKStepGetNumLinRhsEvals"); + + std::cout << " Total number of Newton iterations = " << nni << "\n"; + std::cout << " Total number of linear solver convergence failures = " << ncfn + << "\n"; + std::cout << " Total linear solver setups = " << nsetups << "\n"; + std::cout << " Total number of Jacobian evaluations = " << nje << "\n"; + std::cout << " Total RHS evals for setting up the linear system = " << nfeLS + << "\n"; + } + + if (relax) + { + long int nre, nrje, nrf, nrbf, nrnlsi, nrnlsf; + + flag = ARKStepGetNumRelaxFnEvals(arkode_mem, &nre); + check_flag(flag, "ARKStepGetNumRelaxFnEvals"); + + flag = ARKStepGetNumRelaxJacEvals(arkode_mem, &nrje); + check_flag(flag, "ARKStepGetNumRelaxJacEvals"); + + flag = ARKStepGetNumRelaxFails(arkode_mem, &nrf); + check_flag(flag, "ARKStepGetNumRelaxFails"); + + flag = ARKStepGetNumRelaxBoundFails(arkode_mem, &nrbf); + check_flag(flag, "ARKStepGetNumRelaxBoundFails"); + + flag = ARKStepGetNumRelaxSolveFails(arkode_mem, &nrnlsf); + check_flag(flag, "ARKStepGetNumRelaxSolveFails"); + + flag = ARKStepGetNumRelaxSolveIters(arkode_mem, &nrnlsi); + check_flag(flag, "ARKStepGetNumRelaxSolveIters"); + + std::cout << " Total Relaxation Fn evals = " << nre << "\n"; + std::cout << " Total Relaxation Jac evals = " << nrje << "\n"; + std::cout << " Total Relaxation fails = " << nrf << "\n"; + std::cout << " Total Relaxation bound fails = " << nrbf << "\n"; + std::cout << " Total Relaxation NLS fails = " << nrnlsf << "\n"; + std::cout << " Total Relaxation NLS iters = " << nrnlsi << "\n"; + } + std::cout << "\n"; + + /* -------- * + * Clean up * + * -------- */ + + // Free ARKStep integrator and SUNDIALS objects + ARKStepFree(&arkode_mem); + SUNLinSolFree(LS); + SUNMatDestroy(A); + N_VDestroy(y); + + return flag; +} + +/* ----------------------- * + * User-supplied functions * + * ----------------------- */ + +// ODE RHS function f(t,y). +int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* fdata = N_VGetArrayPointer(ydot); + + fdata[0] = -std::sin(ydata[1]); + fdata[1] = ydata[0]; + + return 0; +} + +// ODE RHS Jacobian function J(t,y) = df/dy. +int Jac(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, + N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* Jdata = SUNDenseMatrix_Data(J); + + // column 0 + Jdata[0] = SUN_RCONST(0.0); + Jdata[1] = SUN_RCONST(1.0); + + // column 1 + Jdata[2] = -std::cos(ydata[1]); + Jdata[3] = SUN_RCONST(0.0); + + return 0; +} + +// Energy function e(y) +int Eng(N_Vector y, sunrealtype* e, void* user_data) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + + *e = SUN_RCONST(0.5) * ydata[0] * ydata[0] - std::cos(ydata[1]); + + return 0; +} + +// Energy function Jacobian Je(y) = de/dy +int JacEng(N_Vector y, N_Vector J, void* user_data) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* jdata = N_VGetArrayPointer(J); + + jdata[0] = ydata[0]; + jdata[1] = std::sin(ydata[1]); + + return 0; +} diff --git a/examples/arkode/CXX_serial/ark_pendulum.out b/examples/arkode/CXX_serial/ark_pendulum.out new file mode 100644 index 0000000000..dca1a58ccb --- /dev/null +++ b/examples/arkode/CXX_serial/ark_pendulum.out @@ -0,0 +1,35 @@ +Nonlinear Pendulum problem: + method = DIRK + relaxation = ON + reltol = 1e-06 + abstol = 1e-10 + + step t u v e e err +--------------------------------------------------------------------------------------------------------------------------- + 0 0.000000000000000e+00 1.500000000000000e+00 1.000000000000000e+00 5.846976941318602e-01 0.000000000000000e+00 + 1000 1.226511879534299e+00 3.458488320982976e-01 2.123384535074106e+00 5.846976941318555e-01 -4.773959005888173e-15 + 2000 2.398784336118560e+00 -6.416442371616857e-01 1.959343228667432e+00 5.846976941318377e-01 -2.253752739989068e-14 + 3000 3.618483293504332e+00 -1.718124366999066e+00 4.706406225355917e-01 5.846976941318435e-01 -1.676436767183986e-14 + 4000 4.580865391430372e+00 -1.411700239944070e+00 -1.146421559917496e+00 5.846976941318369e-01 -2.331468351712829e-14 + 5000 5.877850661141863e+00 -2.026377125271683e-01 -2.170219947005763e+00 5.846976941318408e-01 -1.942890293094024e-14 + 6000 7.079795692175875e+00 8.320973704310155e-01 -1.811622127832077e+00 5.846976941318234e-01 -3.685940441755520e-14 + 7000 8.313428698722841e+00 1.777791699609239e+00 -9.412012907795819e-02 5.846976941318288e-01 -3.141931159689193e-14 + 8000 9.289846347799413e+00 1.192220279192671e+00 1.444463648145965e+00 5.846976941318289e-01 -3.130828929442941e-14 +--------------------------------------------------------------------------------------------------------------------------- + +Final Solver Statistics: + Internal solver steps = 8518 (attempted = 8528) + Total number of error test failures = 10 + Total RHS evals: Fe = 0, Fi = 50909 + Total number of Newton iterations = 25332 + Total number of linear solver convergence failures = 5 + Total linear solver setups = 459 + Total number of Jacobian evaluations = 146 + Total RHS evals for setting up the linear system = 0 + Total Relaxation Fn evals = 25494 + Total Relaxation Jac evals = 25494 + Total Relaxation fails = 0 + Total Relaxation bound fails = 0 + Total Relaxation NLS fails = 0 + Total Relaxation NLS iters = 8438 + diff --git a/examples/arkode/C_parallel/ark_diurnal_kry_bbd_p.out b/examples/arkode/C_parallel/ark_diurnal_kry_bbd_p.out index d3e94e7725..5d53186f79 100644 --- a/examples/arkode/C_parallel/ark_diurnal_kry_bbd_p.out +++ b/examples/arkode/C_parallel/ark_diurnal_kry_bbd_p.out @@ -58,7 +58,7 @@ At top right: c1, c2 = -4.190e-08 4.163e+11 Final Statistics: -lenrw = 3901 leniw = 265 +lenrw = 3901 leniw = 267 lenrwls = 2455 leniwls = 126 nst = 7252 nfe = 0 nfe = 76109 nfels = 82856 @@ -127,7 +127,7 @@ At top right: c1, c2 = -6.607e-08 4.163e+11 Final Statistics: -lenrw = 3901 leniw = 270 +lenrw = 3901 leniw = 272 lenrwls = 2455 leniwls = 126 nst = 7364 nfe = 0 nfe = 77457 nfels = 101131 diff --git a/examples/arkode/C_parallel/ark_diurnal_kry_p.out b/examples/arkode/C_parallel/ark_diurnal_kry_p.out index a3b8276583..97f07c3221 100644 --- a/examples/arkode/C_parallel/ark_diurnal_kry_p.out +++ b/examples/arkode/C_parallel/ark_diurnal_kry_p.out @@ -52,7 +52,7 @@ At top right: c1, c2 = 1.567e-26 4.163e+11 Final Statistics: -lenrw = 3301 leniw = 241 +lenrw = 3301 leniw = 243 lenrwls = 2455 leniwls = 126 nst = 7103 nfe = 0 nfi = 74602 nfels = 87206 diff --git a/examples/arkode/C_parhyp/ark_diurnal_kry_ph.out b/examples/arkode/C_parhyp/ark_diurnal_kry_ph.out index d88203c2b1..411f0a09fc 100644 --- a/examples/arkode/C_parhyp/ark_diurnal_kry_ph.out +++ b/examples/arkode/C_parhyp/ark_diurnal_kry_ph.out @@ -52,7 +52,7 @@ At top right: c1, c2 = -4.442e-20 4.163e+11 Final Statistics: -lenrw = 3301 leniw = 241 +lenrw = 3301 leniw = 243 lenrwls = 2455 leniwls = 126 nst = 7237 nfe = 0 nfi = 75881 nfels = 88695 diff --git a/examples/arkode/C_serial/CMakeLists.txt b/examples/arkode/C_serial/CMakeLists.txt index 5269485065..a9ea064c0c 100644 --- a/examples/arkode/C_serial/CMakeLists.txt +++ b/examples/arkode/C_serial/CMakeLists.txt @@ -26,6 +26,11 @@ set(ARKODE_examples "ark_brusselator\;\;develop" "ark_brusselator_fp\;\;exclude-single" "ark_brusselator1D\;\;exclude-single" + "ark_conserved_exp_entropy_ark\;1 0\;develop" + "ark_conserved_exp_entropy_ark\;1 1\;develop" + "ark_conserved_exp_entropy_erk\;1\;develop" + "ark_dissipated_exp_entropy\;1 0\;develop" + "ark_dissipated_exp_entropy\;1 1\;develop" "ark_heat1D\;\;develop" "ark_heat1D_adapt\;\;develop" "ark_KrylovDemo_prec\;\;exclude-single" @@ -90,14 +95,6 @@ set(ARKODE_extras ark_robertson_stats.csv ) -# Specify libraries to link against -set(ARKODE_LIB sundials_arkode) -set(NVECS_LIB sundials_nvecserial) - -# Set-up linker flags and link libraries -set(SUNDIALS_LIBS ${ARKODE_LIB} ${NVECS_LIB} ${EXE_EXTRA_LINK_LIBS}) - - # Add the build and install targets for each ARKODE example foreach(example_tuple ${ARKODE_examples}) @@ -114,7 +111,10 @@ foreach(example_tuple ${ARKODE_examples}) set_target_properties(${example} PROPERTIES FOLDER "Examples") # libraries to link against - target_link_libraries(${example} ${SUNDIALS_LIBS}) + target_link_libraries(${example} + sundials_arkode + sundials_nvecserial + ${EXE_EXTRA_LINK_LIBS}) endif() # check if example args are provided and set the test name @@ -143,12 +143,6 @@ endforeach(example_tuple ${ARKODE_examples}) # Add the build and install targets for each KLU example (if needed) if(BUILD_SUNLINSOL_KLU) - # Sundials KLU linear solver module - set(SUNLINSOLKLU_LIBS sundials_sunlinsolklu) - - # KLU libraries - list(APPEND SUNLINSOLKLU_LIBS) - foreach(example_tuple ${ARKODE_examples_KLU}) # parse the example tuple @@ -167,7 +161,11 @@ if(BUILD_SUNLINSOL_KLU) EXAMPLE_TYPE ${example_type}) # libraries to link against - target_link_libraries(${example} ${SUNDIALS_LIBS} ${SUNLINSOLKLU_LIBS}) + target_link_libraries(${example} + sundials_arkode + sundials_nvecserial + sundials_sunlinsolklu + ${EXE_EXTRA_LINK_LIBS}) # install example source and out files if(EXAMPLES_INSTALL) @@ -183,12 +181,6 @@ endif() # Add the build and install targets for each SuperLU_MT example (if needed) if(BUILD_SUNLINSOL_SUPERLUMT) - # Sundials SuperLU_MT linear solver module - set(SUNLINSOLSLUMT_LIBS sundials_sunlinsolsuperlumt) - - # SuperLU_MT libraries - list(APPEND SUNLINSOLSLUMT_LIBS ${SUPERLUMT_LIBRARIES}) - foreach(example_tuple ${ARKODE_examples_SUPERLUMT}) # parse the example tuple @@ -207,7 +199,11 @@ if(BUILD_SUNLINSOL_SUPERLUMT) EXAMPLE_TYPE ${example_type}) # libraries to link against - target_link_libraries(${example} ${SUNDIALS_LIBS} ${SUNLINSOLSLUMT_LIBS}) + target_link_libraries(${example} + sundials_arkode + sundials_nvecserial + sundials_sunlinsolsuperlumt + ${EXE_EXTRA_LINK_LIBS}) # install example source and out files if(EXAMPLES_INSTALL) diff --git a/examples/arkode/C_serial/ark_KrylovDemo_prec.out b/examples/arkode/C_serial/ark_KrylovDemo_prec.out index e6390da296..3a63a46785 100644 --- a/examples/arkode/C_serial/ark_KrylovDemo_prec.out +++ b/examples/arkode/C_serial/ark_KrylovDemo_prec.out @@ -419,7 +419,7 @@ Species 6 Final statistics for this run: ARKStep real workspace length = 3557 - ARKStep integer workspace length = 129 + ARKStep integer workspace length = 131 ARKLS real workspace length = 2647 ARKLS integer workspace length = 42 Number of steps = 266 @@ -488,7 +488,7 @@ t = 1.00e+01 nst = 251 nfe = 0 nfi = 4101 nni = 2746 hu = 7.09e-02 Final statistics for this run: ARKStep real workspace length = 3557 - ARKStep integer workspace length = 134 + ARKStep integer workspace length = 136 ARKLS real workspace length = 2647 ARKLS integer workspace length = 42 Number of steps = 251 @@ -557,7 +557,7 @@ t = 1.00e+01 nst = 359 nfe = 0 nfi = 5722 nni = 3907 hu = 1.35e-01 Final statistics for this run: ARKStep real workspace length = 3557 - ARKStep integer workspace length = 139 + ARKStep integer workspace length = 141 ARKLS real workspace length = 2647 ARKLS integer workspace length = 42 Number of steps = 359 @@ -626,7 +626,7 @@ t = 1.00e+01 nst = 359 nfe = 0 nfi = 5728 nni = 3912 hu = 1.31e-01 Final statistics for this run: ARKStep real workspace length = 3557 - ARKStep integer workspace length = 144 + ARKStep integer workspace length = 146 ARKLS real workspace length = 2647 ARKLS integer workspace length = 42 Number of steps = 359 diff --git a/examples/arkode/C_serial/ark_KrylovDemo_prec_1.out b/examples/arkode/C_serial/ark_KrylovDemo_prec_1.out index e6390da296..58b34e703d 100644 --- a/examples/arkode/C_serial/ark_KrylovDemo_prec_1.out +++ b/examples/arkode/C_serial/ark_KrylovDemo_prec_1.out @@ -419,7 +419,7 @@ Species 6 Final statistics for this run: ARKStep real workspace length = 3557 - ARKStep integer workspace length = 129 + ARKStep integer workspace length = 131 ARKLS real workspace length = 2647 ARKLS integer workspace length = 42 Number of steps = 266 @@ -488,7 +488,7 @@ t = 1.00e+01 nst = 251 nfe = 0 nfi = 4101 nni = 2746 hu = 7.09e-02 Final statistics for this run: ARKStep real workspace length = 3557 - ARKStep integer workspace length = 134 + ARKStep integer workspace length = 136 ARKLS real workspace length = 2647 ARKLS integer workspace length = 42 Number of steps = 251 @@ -557,7 +557,7 @@ t = 1.00e+01 nst = 359 nfe = 0 nfi = 5722 nni = 3907 hu = 1.35e-01 Final statistics for this run: ARKStep real workspace length = 3557 - ARKStep integer workspace length = 139 + ARKStep integer workspace length = 141 ARKLS real workspace length = 2647 ARKLS integer workspace length = 42 Number of steps = 359 @@ -626,7 +626,7 @@ t = 1.00e+01 nst = 359 nfe = 0 nfi = 5728 nni = 3912 hu = 1.31e-01 Final statistics for this run: ARKStep real workspace length = 3557 - ARKStep integer workspace length = 144 + ARKStep integer workspace length = 146 ARKLS real workspace length = 2647 ARKLS integer workspace length = 42 Number of steps = 359 diff --git a/examples/arkode/C_serial/ark_KrylovDemo_prec_2.out b/examples/arkode/C_serial/ark_KrylovDemo_prec_2.out index e6390da296..1af6fb8046 100644 --- a/examples/arkode/C_serial/ark_KrylovDemo_prec_2.out +++ b/examples/arkode/C_serial/ark_KrylovDemo_prec_2.out @@ -419,7 +419,7 @@ Species 6 Final statistics for this run: ARKStep real workspace length = 3557 - ARKStep integer workspace length = 129 + ARKStep integer workspace length = 131 ARKLS real workspace length = 2647 ARKLS integer workspace length = 42 Number of steps = 266 @@ -488,7 +488,7 @@ t = 1.00e+01 nst = 251 nfe = 0 nfi = 4101 nni = 2746 hu = 7.09e-02 Final statistics for this run: ARKStep real workspace length = 3557 - ARKStep integer workspace length = 134 + ARKStep integer workspace length = 136 ARKLS real workspace length = 2647 ARKLS integer workspace length = 42 Number of steps = 251 @@ -557,7 +557,7 @@ t = 1.00e+01 nst = 359 nfe = 0 nfi = 5722 nni = 3907 hu = 1.35e-01 Final statistics for this run: ARKStep real workspace length = 3557 - ARKStep integer workspace length = 139 + ARKStep integer workspace length = 141 ARKLS real workspace length = 2647 ARKLS integer workspace length = 42 Number of steps = 359 @@ -626,7 +626,7 @@ t = 1.00e+01 nst = 359 nfe = 0 nfi = 5728 nni = 3912 hu = 1.31e-01 Final statistics for this run: ARKStep real workspace length = 3557 - ARKStep integer workspace length = 144 + ARKStep integer workspace length = 146 ARKLS real workspace length = 2647 ARKLS integer workspace length = 42 Number of steps = 359 diff --git a/examples/arkode/C_serial/ark_conserved_exp_entropy_ark.c b/examples/arkode/C_serial/ark_conserved_exp_entropy_ark.c new file mode 100644 index 0000000000..9064c56987 --- /dev/null +++ b/examples/arkode/C_serial/ark_conserved_exp_entropy_ark.c @@ -0,0 +1,528 @@ +/* ----------------------------------------------------------------------------- + * Programmer(s): David J. Gardner @ LLNL + * ----------------------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2002-2022, 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 + * ----------------------------------------------------------------------------- + * This example problem is adapted from: + * + * H. Ranocha, M. Sayyari, L. Dalcin, M. Parsani, and D.I. Ketcheson, + * "Relaxation Runge-Kutta Methods: Fully-Discrete Explicit Entropy-Stable + * Schemes for the Compressible Euler and Navier-Stokes Equations," SIAM Journal + * on Scientific Computing, 42(2), 2020, https://doi.org/10.1137/19M1263480. + * ----------------------------------------------------------------------------- + * This example evolves system + * + * du/dt = -exp(v) + * dv/dt = exp(u) + * + * for t in the interval [0, 5] with the initial condition + * + * u(0) = 1.0 + * v(0) = 0.5 + * + * The system has the analytic solution + * + * u = log(e + e^(3/2)) - log(b) + * v = log(a * e^(a * t)) - log(b) + * + * where log is the natural logarithm, a = sqrt(e) + e, and + * b = sqrt(e) + e^(a * t). + * + * The conserved exponential entropy for the system is given by + * ent(u,v) = exp(u) + exp(v) with the Jacobian + * ent'(u,v) = [ de/du de/dv ]^T = [ exp(u) exp(v) ]^T. + * + * The problem is advanced in time with an explicit or implicit relaxed + * Runge-Kutta method from ARKStep to ensure conservation of the entropy. + * ---------------------------------------------------------------------------*/ + +/* Header files */ +#include +#include + +/* SUNDIALS headers */ +#include +#include +#include +#include + +/* Value of the natural number e */ +#define EVAL 2.718281828459045235360287471352662497757247093699959574966 + +/* Convince macros for calling precision-specific math functions */ +#if defined(SUNDIALS_DOUBLE_PRECISION) +#define EXP(x) (exp((x))) +#define SQRT(x) (sqrt((x))) +#define LOG(x) (log((x))) +#elif defined(SUNDIALS_SINGLE_PRECISION) +#define EXP(x) (expf((x))) +#define SQRT(x) (sqrtf((x))) +#define LOG(x) (logf((x))) +#elif defined(SUNDIALS_EXTENDED_PRECISION) +#define EXP(x) (expl((x))) +#define SQRT(x) (sqrtl((x))) +#define LOG(x) (logl((x))) +#endif + +/* Convince macros for using precision-specific format specifiers */ +#if defined(SUNDIALS_EXTENDED_PRECISION) +#define GSYM "Lg" +#define ESYM "Le" +#define FSYM "Lf" +#else +#define GSYM "g" +#define ESYM "e" +#define FSYM "f" +#endif + +/* ----------------------- * + * User-supplied functions * + * ----------------------- */ + +/* ODE RHS function */ +int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data); + +/* ODE RHS Jacobian function */ +int Jac(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, + N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); + +/* Entropy function */ +int Ent(N_Vector y, sunrealtype* e, void* user_data); + +/* Entropy Jacobian function */ +int JacEnt(N_Vector y, N_Vector J, void* user_data); + +/* ----------------- * + * Utility functions * + * ----------------- */ + +/* Analytic solution */ +int ans(sunrealtype t, N_Vector y); + +/* Check for an unrecoverable (negative) return flag from a SUNDIALS function */ +int check_flag(int flag, const char* funcname); + +/* Check if a function returned a NULL pointer */ +int check_ptr(void* ptr, const char* funcname); + +/* ------------ * + * Main Program * + * ------------ */ + +int main(int argc, char* argv[]) +{ + /* Error-checking flag */ + int flag; + + /* Initial and final times */ + sunrealtype t0 = SUN_RCONST(0.0); + sunrealtype tf = SUN_RCONST(5.0); + + /* Relative and absolute tolerances */ + sunrealtype reltol = SUN_RCONST(1.0e-6); + sunrealtype abstol = SUN_RCONST(1.0e-10); + + /* SUNDIALS context, vector, matrix, and linear solver objects */ + SUNContext ctx = NULL; + N_Vector y = NULL; + N_Vector ytrue = NULL; + SUNMatrix A = NULL; + SUNLinearSolver LS = NULL; + + /* Pointer to vector data array */ + sunrealtype* ydata; + sunrealtype* ytdata; + + /* Initial, current, and change in entropy value */ + sunrealtype ent0, ent, ent_err; + + /* Solution errors */ + sunrealtype u_err, v_err; + + /* ARKODE memory structure */ + void* arkode_mem = NULL; + + /* ARKODE statistics */ + long int nst, nst_a, nfe, nfi; + long int nrf, nrbf, nre, nrje, nrnlsi, nrnlsf; + long int nsetups, nje, nfeLS, nni, ncfn, netf; + + /* Output time and file */ + sunrealtype t; + FILE* UFID; + + /* Command line options */ + int relax = 1; /* enable relaxation */ + int implicit = 1; /* implicit */ + sunrealtype fixed_h = SUN_RCONST(0.0); /* adaptive stepping */ + + /* -------------------- * + * Output Problem Setup * + * -------------------- */ + + if (argc > 1) relax = atoi(argv[1]); + if (argc > 2) implicit = atoi(argv[2]); + if (argc > 3) fixed_h = atof(argv[3]); + + printf("\nConserved Exponential Entropy problem:\n"); + if (implicit) + { + printf(" method = DIRK\n"); + } + else + { + printf(" method = ERK\n"); + } + printf(" reltol = %.1" ESYM "\n", reltol); + printf(" abstol = %.1" ESYM "\n", abstol); + if (fixed_h > SUN_RCONST(0.0)) + { + printf(" fixed h = %.1" ESYM "\n", fixed_h); + } + if (relax) + { + printf(" relaxation = ON\n"); + } + else + { + printf(" relaxation = OFF\n"); + } + printf("\n"); + + /* ------------ * + * Setup ARKODE * + * ------------ */ + + /* Create the SUNDIALS context object for this simulation */ + flag = SUNContext_Create(NULL, &ctx); + if (check_flag(flag, "SUNContext_Create")) return 1; + + /* Create serial vector and set the initial condition values */ + y = N_VNew_Serial(2, ctx); + if (check_ptr(y, "N_VNew_Serial")) return 1; + + ydata = N_VGetArrayPointer(y); + if (check_ptr(ydata, "N_VGetArrayPointer")) return 1; + + ydata[0] = SUN_RCONST(1.0); + ydata[1] = SUN_RCONST(0.5); + + ytrue = N_VClone(y); + if (check_ptr(ytrue, "N_VClone")) return 1; + + ytdata = N_VGetArrayPointer(ytrue); + if (check_ptr(ytdata, "N_VGetArrayPointer")) return 1; + + /* Initialize ARKStep */ + if (implicit) + { + arkode_mem = ARKStepCreate(NULL, f, t0, y, ctx); + } + else + { + arkode_mem = ARKStepCreate(f, NULL, t0, y, ctx); + } + if (check_ptr(arkode_mem, "ARKStepCreate")) return 1; + + /* Specify tolerances */ + flag = ARKStepSStolerances(arkode_mem, reltol, abstol); + if (check_flag(flag, "ARKStepSStolerances")) return 1; + + if (relax) + { + /* Enable relaxation methods */ + flag = ARKStepSetRelaxFn(arkode_mem, Ent, JacEnt); + if (check_flag(flag, "ARKStepSetRelaxFn")) return 1; + } + + if (implicit) + { + /* Create dense matrix and linear solver */ + A = SUNDenseMatrix(2, 2, ctx); + if (check_ptr(A, "SUNDenseMatrix")) return 1; + + LS = SUNLinSol_Dense(y, A, ctx); + if (check_ptr(LS, "SUNLinSol_Dense")) return 1; + + /* Attach the matrix and linear solver */ + flag = ARKStepSetLinearSolver(arkode_mem, LS, A); + if (check_flag(flag, "ARKStepSetLinearSolver")) return 1; + + /* Set Jacobian routine */ + flag = ARKStepSetJacFn(arkode_mem, Jac); + if (check_flag(flag, "ARKStepSetJacFn")) return 1; + + /* Select a Butcher table with non-negative b values */ + flag = ARKStepSetTableName(arkode_mem, "ARKODE_ARK2_DIRK_3_1_2", + "ARKODE_ERK_NONE"); + if (check_flag(flag, "ARKStepSetTableName")) return 1; + } + + if (fixed_h > SUN_RCONST(0.0)) + { + flag = ARKStepSetFixedStep(arkode_mem, fixed_h); + if (check_flag(flag, "ARKStepSetFixedStep")) return 1; + } + + /* Open output stream for results, output comment line */ + UFID = fopen("ark_conserved_exp_entropy_ark.txt", "w"); + fprintf(UFID, "# vars: t u v entropy u_err v_err entropy_error\n"); + + /* --------------- * + * Advance in Time * + * --------------- */ + + /* Initial time */ + t = t0; + + /* Output the initial condition and entropy */ + flag = Ent(y, &ent0, NULL); + if (check_flag(flag, "Ent")) return 1; + + fprintf(UFID, + "%23.16" ESYM " %23.16" ESYM " %23.16" ESYM " %23.16" ESYM + " %23.16" ESYM " %23.16" ESYM " %23.16" ESYM "\n", + t0, ydata[0], ydata[1], ent0, SUN_RCONST(0.0), SUN_RCONST(0.0), + SUN_RCONST(0.0)); + + printf(" step t u v e " + " delta e\n"); + printf(" --------------------------------------------------------------------" + "-----------\n"); + printf("%5d %14.6" ESYM " %14.6" ESYM " %14.6" ESYM " %14.6" ESYM " %14.6" ESYM + "\n", + 0, t, ydata[0], ydata[1], ent0, SUN_RCONST(0.0)); + + while (t < tf) + { + /* Evolve in time */ + flag = ARKStepEvolve(arkode_mem, tf, y, &t, ARK_ONE_STEP); + if (check_flag(flag, "ARKStepEvolve")) break; + + /* Output solution and errors */ + flag = Ent(y, &ent, NULL); + if (check_flag(flag, "Ent")) return 1; + + flag = ans(t, ytrue); + if (check_flag(flag, "ans")) return 1; + + ent_err = ent - ent0; + u_err = ydata[0] - ytdata[0]; + v_err = ydata[1] - ytdata[1]; + + /* Output to the screen periodically */ + flag = ARKStepGetNumSteps(arkode_mem, &nst); + check_flag(flag, "ARKStepGetNumSteps"); + + if (nst % 40 == 0) + { + printf("%5ld %14.6" ESYM " %14.6" ESYM " %14.6" ESYM " %14.6" ESYM " %14.6" ESYM + "\n", + nst, t, ydata[0], ydata[1], ent, ent_err); + } + + /* Write all steps to file */ + fprintf(UFID, + "%23.16" ESYM " %23.16" ESYM " %23.16" ESYM " %23.16" ESYM + " %23.16" ESYM " %23.16" ESYM " %23.16" ESYM "\n", + t, ydata[0], ydata[1], ent, u_err, v_err, ent_err); + } + + printf(" --------------------------------------------------------------------" + "-----------\n"); + fclose(UFID); + + /* ------------ * + * Output Stats * + * ------------ */ + + /* Get final statistics on how the solve progressed */ + flag = ARKStepGetNumSteps(arkode_mem, &nst); + check_flag(flag, "ARKStepGetNumSteps"); + + flag = ARKStepGetNumStepAttempts(arkode_mem, &nst_a); + check_flag(flag, "ARKStepGetNumStepAttempts"); + + flag = ARKStepGetNumErrTestFails(arkode_mem, &netf); + check_flag(flag, "ARKStepGetNumErrTestFails"); + + flag = ARKStepGetNumRhsEvals(arkode_mem, &nfe, &nfi); + check_flag(flag, "ARKStepGetNumRhsEvals"); + + printf("\nFinal Solver Statistics:\n"); + printf(" Internal solver steps = %li (attempted = %li)\n", nst, nst_a); + printf(" Total number of error test failures = %li\n", netf); + printf(" Total RHS evals: Fe = %li, Fi = %li\n", nfe, nfi); + + if (implicit) + { + flag = ARKStepGetNumNonlinSolvIters(arkode_mem, &nni); + check_flag(flag, "ARKStepGetNumNonlinSolvIters"); + + flag = ARKStepGetNumNonlinSolvConvFails(arkode_mem, &ncfn); + check_flag(flag, "ARKStepGetNumNonlinSolvConvFails"); + + flag = ARKStepGetNumLinSolvSetups(arkode_mem, &nsetups); + check_flag(flag, "ARKStepGetNumLinSolvSetups"); + + flag = ARKStepGetNumJacEvals(arkode_mem, &nje); + check_flag(flag, "ARKStepGetNumJacEvals"); + + flag = ARKStepGetNumLinRhsEvals(arkode_mem, &nfeLS); + check_flag(flag, "ARKStepGetNumLinRhsEvals"); + + printf(" Total number of Newton iterations = %li\n", nni); + printf(" Total number of linear solver convergence failures = %li\n", ncfn); + printf(" Total linear solver setups = %li\n", nsetups); + printf(" Total number of Jacobian evaluations = %li\n", nje); + printf(" Total RHS evals for setting up the linear system = %li\n", nfeLS); + } + + if (relax) + { + flag = ARKStepGetNumRelaxFnEvals(arkode_mem, &nre); + check_flag(flag, "ARKStepGetNumRelaxFnEvals"); + + flag = ARKStepGetNumRelaxJacEvals(arkode_mem, &nrje); + check_flag(flag, "ARKStepGetNumRelaxJacEvals"); + + flag = ARKStepGetNumRelaxFails(arkode_mem, &nrf); + check_flag(flag, "ARKStepGetNumRelaxFails"); + + flag = ARKStepGetNumRelaxBoundFails(arkode_mem, &nrbf); + check_flag(flag, "ARKStepGetNumRelaxBoundFails"); + + flag = ARKStepGetNumRelaxSolveFails(arkode_mem, &nrnlsf); + check_flag(flag, "ARKStepGetNumRelaxSolveFails"); + + flag = ARKStepGetNumRelaxSolveIters(arkode_mem, &nrnlsi); + check_flag(flag, "ARKStepGetNumRelaxSolveIters"); + + printf(" Total Relaxation Fn evals = %li\n", nre); + printf(" Total Relaxation Jac evals = %li\n", nrje); + printf(" Total Relaxation fails = %li\n", nrf); + printf(" Total Relaxation bound fails = %li\n", nrbf); + printf(" Total Relaxation NLS fails = %li\n", nrnlsf); + printf(" Total Relaxation NLS iters = %li\n", nrnlsi); + } + printf("\n"); + + /* -------- * + * Clean up * + * -------- */ + + /* Free ARKStep integrator and SUNDIALS objects */ + ARKStepFree(&arkode_mem); + SUNLinSolFree(LS); + SUNMatDestroy(A); + N_VDestroy(y); + N_VDestroy(ytrue); + SUNContext_Free(&ctx); + + return flag; +} + +/* ----------------------- * + * User-supplied functions * + * ----------------------- */ + +/* ODE RHS function f(t,y). */ +int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* fdata = N_VGetArrayPointer(ydot); + + fdata[0] = -EXP(ydata[1]); + fdata[1] = EXP(ydata[0]); + + return 0; +} + +/* ODE RHS Jacobian function J(t,y) = df/dy. */ +int Jac(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, + N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* Jdata = SUNDenseMatrix_Data(J); + + /* column 0 */ + Jdata[0] = SUN_RCONST(0.0); + Jdata[1] = EXP(ydata[0]); + + /* column 1 */ + Jdata[2] = -EXP(ydata[1]); + Jdata[3] = SUN_RCONST(0.0); + + return 0; +} + +/* Entropy function e(y) */ +int Ent(N_Vector y, sunrealtype* e, void* user_data) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + + *e = EXP(ydata[0]) + EXP(ydata[1]); + + return 0; +} + +/* Entropy function Jacobian Je(y) = de/dy */ +int JacEnt(N_Vector y, N_Vector J, void* user_data) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* jdata = N_VGetArrayPointer(J); + + jdata[0] = EXP(ydata[0]); + jdata[1] = EXP(ydata[1]); + + return 0; +} + +/* ----------------- * + * Utility functions * + * ----------------- */ + +/* Analytic solution */ +int ans(sunrealtype t, N_Vector y) +{ + sunrealtype a, b; + sunrealtype* ydata = N_VGetArrayPointer(y); + + a = SQRT(EVAL) + EVAL; + b = SQRT(EVAL) + EXP(a * t); + + ydata[0] = LOG(EVAL + EXP(SUN_RCONST(1.5))) - LOG(b); + ydata[1] = LOG(a * EXP(a * t)) - LOG(b); + + return 0; +} + +/* Check for an unrecoverable (negative) return flag from a SUNDIALS function */ +int check_flag(int flag, const char* funcname) +{ + if (flag < 0) + { + fprintf(stderr, "ERROR: %s() returned %d\n", funcname, flag); + return 1; + } + return 0; +} + +/* Check if a function returned a NULL pointer */ +int check_ptr(void* ptr, const char* funcname) +{ + if (!ptr) + { + fprintf(stderr, "ERROR: %s() returned NULL\n", funcname); + return 1; + } + return 0; +} diff --git a/examples/arkode/C_serial/ark_conserved_exp_entropy_ark_1_0.out b/examples/arkode/C_serial/ark_conserved_exp_entropy_ark_1_0.out new file mode 100644 index 0000000000..05f6b3f4ce --- /dev/null +++ b/examples/arkode/C_serial/ark_conserved_exp_entropy_ark_1_0.out @@ -0,0 +1,25 @@ + +Conserved Exponential Entropy problem: + method = ERK + reltol = 1.0e-06 + abstol = 1.0e-10 + relaxation = ON + + step t u v e delta e + ------------------------------------------------------------------------------- + 0 0.000000e+00 1.000000e+00 5.000000e-01 4.367003e+00 0.000000e+00 + 40 6.911091e-01 -1.121528e+00 1.396547e+00 4.367003e+00 -1.776357e-15 + 80 1.946408e+00 -6.526227e+00 1.473742e+00 4.367003e+00 4.440892e-15 + ------------------------------------------------------------------------------- + +Final Solver Statistics: + Internal solver steps = 99 (attempted = 124) + Total number of error test failures = 0 + Total RHS evals: Fe = 722, Fi = 0 + Total Relaxation Fn evals = 622 + Total Relaxation Jac evals = 1019 + Total Relaxation fails = 25 + Total Relaxation bound fails = 0 + Total Relaxation NLS fails = 25 + Total Relaxation NLS iters = 399 + diff --git a/examples/arkode/C_serial/ark_conserved_exp_entropy_ark_1_1.out b/examples/arkode/C_serial/ark_conserved_exp_entropy_ark_1_1.out new file mode 100644 index 0000000000..77742dc2bd --- /dev/null +++ b/examples/arkode/C_serial/ark_conserved_exp_entropy_ark_1_1.out @@ -0,0 +1,53 @@ + +Conserved Exponential Entropy problem: + method = DIRK + reltol = 1.0e-06 + abstol = 1.0e-10 + relaxation = ON + + step t u v e delta e + ------------------------------------------------------------------------------- + 0 0.000000e+00 1.000000e+00 5.000000e-01 4.367003e+00 0.000000e+00 + 40 4.754706e-02 9.164641e-01 6.241023e-01 4.367003e+00 8.881784e-16 + 80 8.098591e-02 8.514226e-01 7.050883e-01 4.367003e+00 0.000000e+00 + 120 1.285350e-01 7.498036e-01 8.111164e-01 4.367003e+00 1.776357e-15 + 160 1.711613e-01 6.495637e-01 8.970255e-01 4.367003e+00 8.881784e-16 + 200 2.117650e-01 5.461535e-01 9.709320e-01 4.367003e+00 5.329071e-15 + 240 2.453290e-01 4.549914e-01 1.026344e+00 4.367003e+00 0.000000e+00 + 280 2.845614e-01 3.421715e-01 1.084852e+00 4.367003e+00 1.776357e-15 + 320 3.187813e-01 2.385199e-01 1.130639e+00 4.367003e+00 1.776357e-15 + 360 3.470183e-01 1.495118e-01 1.164942e+00 4.367003e+00 4.440892e-15 + 400 3.679487e-01 8.161060e-02 1.188444e+00 4.367003e+00 1.776357e-15 + 440 3.849834e-01 2.519301e-02 1.206417e+00 4.367003e+00 7.105427e-15 + 480 3.919101e-01 1.965927e-03 1.213438e+00 4.367003e+00 3.552714e-15 + 520 3.949458e-01 -8.265002e-03 1.216465e+00 4.367003e+00 2.664535e-15 + 560 4.054406e-01 -4.387006e-02 1.226690e+00 4.367003e+00 -1.776357e-15 + 600 4.224627e-01 -1.023806e-01 1.242516e+00 4.367003e+00 0.000000e+00 + 640 4.483804e-01 -1.931947e-01 1.264884e+00 4.367003e+00 -1.776357e-15 + 680 4.858087e-01 -3.277677e-01 1.293760e+00 4.367003e+00 1.776357e-15 + 720 5.292378e-01 -4.885049e-01 1.322678e+00 4.367003e+00 -4.440892e-15 + 760 5.861685e-01 -7.057075e-01 1.354092e+00 4.367003e+00 -6.217249e-15 + 800 6.559300e-01 -9.802131e-01 1.384235e+00 4.367003e+00 -2.664535e-15 + 840 7.351394e-01 -1.300674e+00 1.409682e+00 4.367003e+00 8.881784e-16 + 880 8.653344e-01 -1.841818e+00 1.437100e+00 4.367003e+00 -8.881784e-16 + 920 1.014019e+00 -2.473635e+00 1.454590e+00 4.367003e+00 8.881784e-16 + 960 1.342023e+00 -3.891227e+00 1.469390e+00 4.367003e+00 3.552714e-15 + 1000 2.028294e+00 -6.883726e+00 1.473842e+00 4.367003e+00 5.329071e-15 + ------------------------------------------------------------------------------- + +Final Solver Statistics: + Internal solver steps = 1024 (attempted = 1376) + Total number of error test failures = 86 + Total RHS evals: Fe = 0, Fi = 9934 + Total number of Newton iterations = 5803 + Total number of linear solver convergence failures = 37 + Total linear solver setups = 844 + Total number of Jacobian evaluations = 304 + Total RHS evals for setting up the linear system = 0 + Total Relaxation Fn evals = 7023 + Total Relaxation Jac evals = 8651 + Total Relaxation fails = 266 + Total Relaxation bound fails = 14 + Total Relaxation NLS fails = 252 + Total Relaxation NLS iters = 4523 + diff --git a/examples/arkode/C_serial/ark_conserved_exp_entropy_erk.c b/examples/arkode/C_serial/ark_conserved_exp_entropy_erk.c new file mode 100644 index 0000000000..a13add2ac3 --- /dev/null +++ b/examples/arkode/C_serial/ark_conserved_exp_entropy_erk.c @@ -0,0 +1,436 @@ +/* ----------------------------------------------------------------------------- + * Programmer(s): David J. Gardner @ LLNL + * ----------------------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2002-2022, 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 + * ----------------------------------------------------------------------------- + * This example problem is adapted from: + * + * H. Ranocha, M. Sayyari, L. Dalcin, M. Parsani, and D.I. Ketcheson, + * "Relaxation Runge-Kutta Methods: Fully-Discrete Explicit Entropy-Stable + * Schemes for the Compressible Euler and Navier-Stokes Equations," SIAM Journal + * on Scientific Computing, 42(2), 2020, https://doi.org/10.1137/19M1263480. + * ----------------------------------------------------------------------------- + * This example evolves system + * + * du/dt = -exp(v) + * dv/dt = exp(u) + * + * for t in the interval [0, 5] with the initial condition + * + * u(0) = 1.0 + * v(0) = 0.5 + * + * The system has the analytic solution + * + * u = log(e + e^(3/2)) - log(b) + * v = log(a * e^(a * t)) - log(b) + * + * where log is the natural logarithm, a = sqrt(e) + e, and + * b = sqrt(e) + e^(a * t). + * + * The conserved exponential entropy for the system is given by + * ent(u,v) = exp(u) + exp(v) with the Jacobian + * ent'(u,v) = [ de/du de/dv ]^T = [ exp(u) exp(v) ]^T. + * + * The problem is advanced in time with an explicit relaxed + * Runge-Kutta method from ERKStep to ensure conservation of the entropy. + * ---------------------------------------------------------------------------*/ + +/* Header files */ +#include +#include + +/* SUNDIALS headers */ +#include +#include + +/* Value of the natural number e */ +#define EVAL 2.718281828459045235360287471352662497757247093699959574966 + +/* Convince macros for calling precision-specific math functions */ +#if defined(SUNDIALS_DOUBLE_PRECISION) +#define EXP(x) (exp((x))) +#define SQRT(x) (sqrt((x))) +#define LOG(x) (log((x))) +#elif defined(SUNDIALS_SINGLE_PRECISION) +#define EXP(x) (expf((x))) +#define SQRT(x) (sqrtf((x))) +#define LOG(x) (logf((x))) +#elif defined(SUNDIALS_EXTENDED_PRECISION) +#define EXP(x) (expl((x))) +#define SQRT(x) (sqrtl((x))) +#define LOG(x) (logl((x))) +#endif + +/* Convince macros for using precision-specific format specifiers */ +#if defined(SUNDIALS_EXTENDED_PRECISION) +#define GSYM "Lg" +#define ESYM "Le" +#define FSYM "Lf" +#else +#define GSYM "g" +#define ESYM "e" +#define FSYM "f" +#endif + +/* ----------------------- * + * User-supplied functions * + * ----------------------- */ + +/* ODE RHS function */ +int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data); + +/* Entropy function */ +int Ent(N_Vector y, sunrealtype* e, void* user_data); + +/* Entropy Jacobian function */ +int JacEnt(N_Vector y, N_Vector J, void* user_data); + +/* ----------------- * + * Utility functions * + * ----------------- */ + +/* Analytic solution */ +int ans(sunrealtype t, N_Vector y); + +/* Check for an unrecoverable (negative) return flag from a SUNDIALS function */ +int check_flag(int flag, const char* funcname); + +/* Check if a function returned a NULL pointer */ +int check_ptr(void* ptr, const char* funcname); + +/* ------------ * + * Main Program * + * ------------ */ + +int main(int argc, char* argv[]) +{ + /* Error-checking flag */ + int flag; + + /* Initial and final times */ + sunrealtype t0 = SUN_RCONST(0.0); + sunrealtype tf = SUN_RCONST(5.0); + + /* Relative and absolute tolerances */ + sunrealtype reltol = SUN_RCONST(1.0e-6); + sunrealtype abstol = SUN_RCONST(1.0e-10); + + /* SUNDIALS context, vector, matrix, and linear solver objects */ + SUNContext ctx = NULL; + N_Vector y = NULL; + N_Vector ytrue = NULL; + + /* Pointer to vector data array */ + sunrealtype* ydata; + sunrealtype* ytdata; + + /* Initial, current, and change in entropy value */ + sunrealtype ent0, ent, ent_err; + + /* Solution errors */ + sunrealtype u_err, v_err; + + /* ARKODE memory structure */ + void* arkode_mem = NULL; + + /* ARKODE statistics */ + long int nst, nst_a, nfe; + long int nrf, nrbf, nre, nrje, nrnlsi, nrnlsf, netf; + + /* Output time and file */ + sunrealtype t; + FILE* UFID; + + /* Command line options */ + int relax = 1; /* enable relaxation */ + sunrealtype fixed_h = SUN_RCONST(0.0); /* adaptive stepping */ + + /* -------------------- * + * Output Problem Setup * + * -------------------- */ + + if (argc > 1) relax = atoi(argv[1]); + if (argc > 2) fixed_h = atof(argv[2]); + + printf("\nConserved Exponential Entropy problem:\n"); + printf(" method = ERK\n"); + printf(" reltol = %.1" ESYM "\n", reltol); + printf(" abstol = %.1" ESYM "\n", abstol); + if (fixed_h > SUN_RCONST(0.0)) + { + printf(" fixed h = %.1" ESYM "\n", fixed_h); + } + if (relax) + { + printf(" relaxation = ON\n"); + } + else + { + printf(" relaxation = OFF\n"); + } + printf("\n"); + + /* ------------ * + * Setup ARKODE * + * ------------ */ + + /* Create the SUNDIALS context object for this simulation */ + flag = SUNContext_Create(NULL, &ctx); + if (check_flag(flag, "SUNContext_Create")) return 1; + + /* Create serial vector and set the initial condition values */ + y = N_VNew_Serial(2, ctx); + if (check_ptr(y, "N_VNew_Serial")) return 1; + + ydata = N_VGetArrayPointer(y); + if (check_ptr(ydata, "N_VGetArrayPointer")) return 1; + + ydata[0] = SUN_RCONST(1.0); + ydata[1] = SUN_RCONST(0.5); + + ytrue = N_VClone(y); + if (check_ptr(ytrue, "N_VClone")) return 1; + + ytdata = N_VGetArrayPointer(ytrue); + if (check_ptr(ytdata, "N_VGetArrayPointer")) return 1; + + /* Initialize ERKStep */ + arkode_mem = ERKStepCreate(f, t0, y, ctx); + if (check_ptr(arkode_mem, "ERKStepCreate")) return 1; + + /* Specify tolerances */ + flag = ERKStepSStolerances(arkode_mem, reltol, abstol); + if (check_flag(flag, "ERKStepSStolerances")) return 1; + + if (relax) + { + /* Enable relaxation methods */ + flag = ERKStepSetRelaxFn(arkode_mem, Ent, JacEnt); + if (check_flag(flag, "ERKStepSetRelaxFn")) return 1; + } + + if (fixed_h > SUN_RCONST(0.0)) + { + flag = ERKStepSetFixedStep(arkode_mem, fixed_h); + if (check_flag(flag, "ERKStepSetFixedStep")) return 1; + } + + /* Open output stream for results, output comment line */ + UFID = fopen("ark_conserved_exp_entropy_erk.txt", "w"); + fprintf(UFID, "# vars: t u v entropy u_err v_err entropy_error\n"); + + /* --------------- * + * Advance in Time * + * --------------- */ + + /* Initial time */ + t = t0; + + /* Output the initial condition and entropy */ + flag = Ent(y, &ent0, NULL); + if (check_flag(flag, "Ent")) return 1; + + fprintf(UFID, + "%23.16" ESYM " %23.16" ESYM " %23.16" ESYM " %23.16" ESYM + " %23.16" ESYM " %23.16" ESYM " %23.16" ESYM "\n", + t0, ydata[0], ydata[1], ent0, SUN_RCONST(0.0), SUN_RCONST(0.0), + SUN_RCONST(0.0)); + + printf(" step t u v e " + " delta e\n"); + printf(" --------------------------------------------------------------------" + "-----------\n"); + printf("%5d %14.6" ESYM " %14.6" ESYM " %14.6" ESYM " %14.6" ESYM " %14.6" ESYM + "\n", + 0, t, ydata[0], ydata[1], ent0, SUN_RCONST(0.0)); + + while (t < tf) + { + /* Evolve in time */ + flag = ERKStepEvolve(arkode_mem, tf, y, &t, ARK_ONE_STEP); + if (check_flag(flag, "ERKStepEvolve")) break; + + /* Output solution and errors */ + flag = Ent(y, &ent, NULL); + if (check_flag(flag, "Ent")) return 1; + + flag = ans(t, ytrue); + if (check_flag(flag, "ans")) return 1; + + ent_err = ent - ent0; + u_err = ydata[0] - ytdata[0]; + v_err = ydata[1] - ytdata[1]; + + /* Output to the screen periodically */ + flag = ERKStepGetNumSteps(arkode_mem, &nst); + check_flag(flag, "ERKStepGetNumSteps"); + + if (nst % 40 == 0) + { + printf("%5ld %14.6" ESYM " %14.6" ESYM " %14.6" ESYM " %14.6" ESYM " %14.6" ESYM + "\n", + nst, t, ydata[0], ydata[1], ent, ent_err); + } + + /* Write all steps to file */ + fprintf(UFID, + "%23.16" ESYM " %23.16" ESYM " %23.16" ESYM " %23.16" ESYM + " %23.16" ESYM " %23.16" ESYM " %23.16" ESYM "\n", + t, ydata[0], ydata[1], ent, u_err, v_err, ent_err); + } + + printf(" --------------------------------------------------------------------" + "-----------\n"); + fclose(UFID); + + /* ------------ * + * Output Stats * + * ------------ */ + + /* Get final statistics on how the solve progressed */ + flag = ERKStepGetNumSteps(arkode_mem, &nst); + check_flag(flag, "ERKStepGetNumSteps"); + + flag = ERKStepGetNumStepAttempts(arkode_mem, &nst_a); + check_flag(flag, "ERKStepGetNumStepAttempts"); + + flag = ERKStepGetNumErrTestFails(arkode_mem, &netf); + check_flag(flag, "ERKStepGetNumErrTestFails"); + + flag = ERKStepGetNumRhsEvals(arkode_mem, &nfe); + check_flag(flag, "ERKStepGetNumRhsEvals"); + + printf("\nFinal Solver Statistics:\n"); + printf(" Internal solver steps = %li (attempted = %li)\n", nst, nst_a); + printf(" Total number of error test failures = %li\n", netf); + printf(" Total RHS evals = %li\n", nfe); + + if (relax) + { + flag = ERKStepGetNumRelaxFnEvals(arkode_mem, &nre); + check_flag(flag, "ERKStepGetNumRelaxFnEvals"); + + flag = ERKStepGetNumRelaxJacEvals(arkode_mem, &nrje); + check_flag(flag, "ERKStepGetNumRelaxJacEvals"); + + flag = ERKStepGetNumRelaxFails(arkode_mem, &nrf); + check_flag(flag, "ERKStepGetNumRelaxFails"); + + flag = ERKStepGetNumRelaxBoundFails(arkode_mem, &nrbf); + check_flag(flag, "ERKStepGetNumRelaxBoundFails"); + + flag = ERKStepGetNumRelaxSolveFails(arkode_mem, &nrnlsf); + check_flag(flag, "ERKStepGetNumRelaxSolveFails"); + + flag = ERKStepGetNumRelaxSolveIters(arkode_mem, &nrnlsi); + check_flag(flag, "ERKStepGetNumRelaxSolveIters"); + + printf(" Total Relaxation Fn evals = %li\n", nre); + printf(" Total Relaxation Jac evals = %li\n", nrje); + printf(" Total Relaxation fails = %li\n", nrf); + printf(" Total Relaxation bound fails = %li\n", nrbf); + printf(" Total Relaxation NLS fails = %li\n", nrnlsf); + printf(" Total Relaxation NLS iters = %li\n", nrnlsi); + } + printf("\n"); + + /* -------- * + * Clean up * + * -------- */ + + /* Free ERKStep integrator and SUNDIALS objects */ + ERKStepFree(&arkode_mem); + N_VDestroy(y); + N_VDestroy(ytrue); + SUNContext_Free(&ctx); + + return flag; +} + +/* ----------------------- * + * User-supplied functions * + * ----------------------- */ + +/* ODE RHS function f(t,y). */ +int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* fdata = N_VGetArrayPointer(ydot); + + fdata[0] = -EXP(ydata[1]); + fdata[1] = EXP(ydata[0]); + + return 0; +} + +/* Entropy function e(y) */ +int Ent(N_Vector y, sunrealtype* e, void* user_data) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + + *e = EXP(ydata[0]) + EXP(ydata[1]); + + return 0; +} + +/* Entropy function Jacobian Je(y) = de/dy */ +int JacEnt(N_Vector y, N_Vector J, void* user_data) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* jdata = N_VGetArrayPointer(J); + + jdata[0] = EXP(ydata[0]); + jdata[1] = EXP(ydata[1]); + + return 0; +} + +/* ----------------- * + * Utility functions * + * ----------------- */ + +/* Analytic solution */ +int ans(sunrealtype t, N_Vector y) +{ + sunrealtype a, b; + sunrealtype* ydata = N_VGetArrayPointer(y); + + a = SQRT(EVAL) + EVAL; + b = SQRT(EVAL) + EXP(a * t); + + ydata[0] = LOG(EVAL + EXP(SUN_RCONST(1.5))) - LOG(b); + ydata[1] = LOG(a * EXP(a * t)) - LOG(b); + + return 0; +} + +/* Check for an unrecoverable (negative) return flag from a SUNDIALS function */ +int check_flag(int flag, const char* funcname) +{ + if (flag < 0) + { + fprintf(stderr, "ERROR: %s() returned %d\n", funcname, flag); + return 1; + } + return 0; +} + +/* Check if a function returned a NULL pointer */ +int check_ptr(void* ptr, const char* funcname) +{ + if (!ptr) + { + fprintf(stderr, "ERROR: %s() returned NULL\n", funcname); + return 1; + } + return 0; +} diff --git a/examples/arkode/C_serial/ark_conserved_exp_entropy_erk_1.out b/examples/arkode/C_serial/ark_conserved_exp_entropy_erk_1.out new file mode 100644 index 0000000000..9e1649f48e --- /dev/null +++ b/examples/arkode/C_serial/ark_conserved_exp_entropy_erk_1.out @@ -0,0 +1,25 @@ + +Conserved Exponential Entropy problem: + method = ERK + reltol = 1.0e-06 + abstol = 1.0e-10 + relaxation = ON + + step t u v e delta e + ------------------------------------------------------------------------------- + 0 0.000000e+00 1.000000e+00 5.000000e-01 4.367003e+00 0.000000e+00 + 40 6.759227e-01 -1.060311e+00 1.391445e+00 4.367003e+00 -1.776357e-15 + 80 1.694830e+00 -5.428257e+00 1.473071e+00 4.367003e+00 -8.881784e-16 + ------------------------------------------------------------------------------- + +Final Solver Statistics: + Internal solver steps = 95 (attempted = 121) + Total number of error test failures = 2 + Total RHS evals = 582 + Total Relaxation Fn evals = 608 + Total Relaxation Jac evals = 995 + Total Relaxation fails = 24 + Total Relaxation bound fails = 0 + Total Relaxation NLS fails = 24 + Total Relaxation NLS iters = 390 + diff --git a/examples/arkode/C_serial/ark_dissipated_exp_entropy.c b/examples/arkode/C_serial/ark_dissipated_exp_entropy.c new file mode 100644 index 0000000000..76da27a0c1 --- /dev/null +++ b/examples/arkode/C_serial/ark_dissipated_exp_entropy.c @@ -0,0 +1,478 @@ +/* ----------------------------------------------------------------------------- + * Programmer(s): David J. Gardner @ LLNL + * ----------------------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2002-2022, 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 + * ----------------------------------------------------------------------------- + * This example problem is adapted from: + * + * H. Ranocha, M. Sayyari, L. Dalcin, M. Parsani, and D.I. Ketcheson, + * "Relaxation Runge-Kutta Methods: Fully-Discrete Explicit Entropy-Stable + * Schemes for the Compressible Euler and Navier-Stokes Equations," SIAM Journal + * on Scientific Computing, 42(2), 2020, https://doi.org/10.1137/19M1263480. + * ----------------------------------------------------------------------------- + * This example evolves the equation du/dt = -exp(u) for t in the interval + * [0, 5] with the initial condition u(0) = 0.5. The equation has the analytic + * solution u(t) = -log(e^{-0.5} + t) and the dissipated exponential entropy is + * given by ent(u) = exp(u) with Jacobian ent'(u) = de/du = exp(u). + * + * The problem is advanced in time with an explicit or implicit relaxed + * Runge-Kutta method to ensure dissipation of the entropy. + * ---------------------------------------------------------------------------*/ + +/* Header files */ +#include +#include + +/* SUNDIALS headers */ +#include +#include +#include +#include + +/* Convince macros for calling precision-specific math functions */ +#if defined(SUNDIALS_DOUBLE_PRECISION) +#define EXP(x) (exp((x))) +#define SQRT(x) (sqrt((x))) +#define LOG(x) (log((x))) +#elif defined(SUNDIALS_SINGLE_PRECISION) +#define EXP(x) (expf((x))) +#define SQRT(x) (sqrtf((x))) +#define LOG(x) (logf((x))) +#elif defined(SUNDIALS_EXTENDED_PRECISION) +#define EXP(x) (expl((x))) +#define SQRT(x) (sqrtl((x))) +#define LOG(x) (logl((x))) +#endif + +/* Convince macros for using precision-specific format specifiers */ +#if defined(SUNDIALS_EXTENDED_PRECISION) +#define GSYM "Lg" +#define ESYM "Le" +#define FSYM "Lf" +#else +#define GSYM "g" +#define ESYM "e" +#define FSYM "f" +#endif + +/* ----------------------- * + * User-supplied functions * + * ----------------------- */ + +/* ODE RHS function */ +int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data); + +/* ODE RHS Jacobian function */ +int Jac(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, + N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); + +/* Entropy function */ +int Ent(N_Vector y, sunrealtype* e, void* user_data); + +/* Entropy Jacobian function */ +int JacEnt(N_Vector y, N_Vector J, void* user_data); + +/* ----------------- * + * Utility functions * + * ----------------- */ + +/* Analytic solution */ +int ans(sunrealtype t, N_Vector y); + +/* Check for an unrecoverable (negative) return flag from a SUNDIALS function */ +int check_flag(int flag, const char* funcname); + +/* Check if a function returned a NULL pointer */ +int check_ptr(void* ptr, const char* funcname); + +/* ------------ * + * Main Program * + * ------------ */ + +int main(int argc, char* argv[]) +{ + /* Error-checking flag */ + int flag; + + /* Initial and final times */ + sunrealtype t0 = SUN_RCONST(0.0); + sunrealtype tf = SUN_RCONST(5.0); + + /* Relative and absolute tolerances */ + sunrealtype reltol = SUN_RCONST(1.0e-6); + sunrealtype abstol = SUN_RCONST(1.0e-10); + + /* SUNDIALS context, vector, matrix, and linear solver objects */ + SUNContext ctx = NULL; + N_Vector y = NULL; + N_Vector ytrue = NULL; + SUNMatrix A = NULL; + SUNLinearSolver LS = NULL; + + /* Pointer to vector data array */ + sunrealtype* ydata; + sunrealtype* ytdata; + + /* Initial, current, and change in entropy value */ + sunrealtype ent0, ent, delta_ent; + + /* Solution errors */ + sunrealtype u_err; + + /* ARKODE memory structure */ + void* arkode_mem = NULL; + + /* ARKODE statistics */ + long int nst, nst_a, nfe, nfi; + long int nrf, nrbf, nre, nrje, nrnlsi, nrnlsf; + long int nsetups, nje, nfeLS, nni, ncfn, netf; + + /* Output time and file */ + sunrealtype t; + FILE* UFID; + + /* Command line options */ + int relax = 1; /* enable relaxation */ + int implicit = 1; /* implicit */ + sunrealtype fixed_h = SUN_RCONST(0.0); /* adaptive stepping */ + + /* -------------------- * + * Output Problem Setup * + * -------------------- */ + + if (argc > 1) relax = atoi(argv[1]); + if (argc > 2) implicit = atoi(argv[2]); + if (argc > 3) fixed_h = atof(argv[3]); + + printf("\nDissipated Exponential Entropy problem:\n"); + if (implicit) + { + printf(" method = DIRK\n"); + } + else + { + printf(" method = ERK\n"); + } + printf(" reltol = %.1" ESYM "\n", reltol); + printf(" abstol = %.1" ESYM "\n", abstol); + if (fixed_h > SUN_RCONST(0.0)) + { + printf(" fixed h = %.1" ESYM "\n", fixed_h); + } + if (relax) + { + printf(" relaxation = ON\n"); + } + else + { + printf(" relaxation = OFF\n"); + } + printf("\n"); + + /* ------------ * + * Setup ARKODE * + * ------------ */ + + /* Create the SUNDIALS context object for this simulation */ + flag = SUNContext_Create(NULL, &ctx); + if (check_flag(flag, "SUNContext_Create")) return 1; + + /* Create serial vector and set the initial condition values */ + y = N_VNew_Serial(2, ctx); + if (check_ptr(y, "N_VNew_Serial")) return 1; + + ydata = N_VGetArrayPointer(y); + if (check_ptr(ydata, "N_VGetArrayPointer")) return 1; + + ydata[0] = SUN_RCONST(1.0); + ydata[1] = SUN_RCONST(0.5); + + ytrue = N_VClone(y); + if (check_ptr(ytrue, "N_VClone")) return 1; + + ytdata = N_VGetArrayPointer(ytrue); + if (check_ptr(ytdata, "N_VGetArrayPointer")) return 1; + + /* Initialize ARKStep */ + if (implicit) + { + arkode_mem = ARKStepCreate(NULL, f, t0, y, ctx); + } + else + { + arkode_mem = ARKStepCreate(f, NULL, t0, y, ctx); + } + if (check_ptr(arkode_mem, "ARKStepCreate")) return 1; + + /* Specify tolerances */ + flag = ARKStepSStolerances(arkode_mem, reltol, abstol); + if (check_flag(flag, "ARKStepSStolerances")) return 1; + + if (relax) + { + /* Enable relaxation methods */ + flag = ARKStepSetRelaxFn(arkode_mem, Ent, JacEnt); + if (check_flag(flag, "ARKStepSetRelaxFn")) return 1; + } + + if (implicit) + { + /* Create dense matrix and linear solver */ + A = SUNDenseMatrix(2, 2, ctx); + if (check_ptr(A, "SUNDenseMatrix")) return 1; + + LS = SUNLinSol_Dense(y, A, ctx); + if (check_ptr(LS, "SUNLinSol_Dense")) return 1; + + /* Attach the matrix and linear solver */ + flag = ARKStepSetLinearSolver(arkode_mem, LS, A); + if (check_flag(flag, "ARKStepSetLinearSolver")) return 1; + + /* Set Jacobian routine */ + flag = ARKStepSetJacFn(arkode_mem, Jac); + if (check_flag(flag, "ARKStepSetJacFn")) return 1; + + /* Select a Butcher table with non-negative b values */ + flag = ARKStepSetTableName(arkode_mem, "ARKODE_ARK2_DIRK_3_1_2", + "ARKODE_ERK_NONE"); + if (check_flag(flag, "ARKStepSetTableName")) return 1; + } + + if (fixed_h > SUN_RCONST(0.0)) + { + flag = ARKStepSetFixedStep(arkode_mem, fixed_h); + if (check_flag(flag, "ARKStepSetFixedStep")) return 1; + } + + /* Open output stream for results, output comment line */ + UFID = fopen("ark_dissipated_exp_entropy.txt", "w"); + fprintf(UFID, "# vars: t u entropy u_err delta_entropy\n"); + + /* --------------- * + * Advance in Time * + * --------------- */ + + /* Initial time */ + t = t0; + + /* Output the initial condition and entropy */ + flag = Ent(y, &ent0, NULL); + if (check_flag(flag, "Ent")) return 1; + + fprintf(UFID, + "%23.16" ESYM " %23.16" ESYM " %23.16" ESYM " %23.16" ESYM " %23.16" ESYM "\n", + t0, ydata[0], ent0, SUN_RCONST(0.0), SUN_RCONST(0.0)); + + printf(" step t u e u_err delta e\n"); + printf(" --------------------------------------------------------------------" + "-----------\n"); + printf("%5d %14.6" ESYM " %14.6" ESYM " %14.6" ESYM " %14.6" ESYM " %14.6" ESYM "\n", + 0, t, ydata[0], ent0, SUN_RCONST(0.0), SUN_RCONST(0.0)); + + while (t < tf) + { + /* Evolve in time */ + flag = ARKStepEvolve(arkode_mem, tf, y, &t, ARK_ONE_STEP); + if (check_flag(flag, "ARKStepEvolve")) break; + + /* Output solution and errors */ + flag = Ent(y, &ent, NULL); + if (check_flag(flag, "Ent")) return 1; + + flag = ans(t, ytrue); + if (check_flag(flag, "ans")) return 1; + + delta_ent = ent - ent0; + u_err = ydata[0] - ytdata[0]; + + /* Output to the screen periodically */ + flag = ARKStepGetNumSteps(arkode_mem, &nst); + check_flag(flag, "ARKStepGetNumSteps"); + + if (nst % 40 == 0) + { + printf("%5ld %14.6" ESYM " %14.6" ESYM " %14.6" ESYM " %14.6" ESYM " %14.6" ESYM "\n", + nst, t, ydata[0], ent, u_err, delta_ent); + } + + /* Write all steps to file */ + fprintf(UFID, + "%23.16" ESYM " %23.16" ESYM " %23.16" ESYM " %23.16" ESYM " %23.16" ESYM "\n", + t, ydata[0], ent, u_err, delta_ent); + } + + printf(" --------------------------------------------------------------------" + "-----------\n"); + fclose(UFID); + + /* ------------ * + * Output Stats * + * ------------ */ + + /* Get final statistics on how the solve progressed */ + flag = ARKStepGetNumSteps(arkode_mem, &nst); + check_flag(flag, "ARKStepGetNumSteps"); + + flag = ARKStepGetNumStepAttempts(arkode_mem, &nst_a); + check_flag(flag, "ARKStepGetNumStepAttempts"); + + flag = ARKStepGetNumErrTestFails(arkode_mem, &netf); + check_flag(flag, "ARKStepGetNumErrTestFails"); + + flag = ARKStepGetNumRhsEvals(arkode_mem, &nfe, &nfi); + check_flag(flag, "ARKStepGetNumRhsEvals"); + + printf("\nFinal Solver Statistics:\n"); + printf(" Internal solver steps = %li (attempted = %li)\n", nst, nst_a); + printf(" Total number of error test failures = %li\n", netf); + printf(" Total RHS evals: Fe = %li, Fi = %li\n", nfe, nfi); + + if (implicit) + { + flag = ARKStepGetNumNonlinSolvIters(arkode_mem, &nni); + check_flag(flag, "ARKStepGetNumNonlinSolvIters"); + + flag = ARKStepGetNumNonlinSolvConvFails(arkode_mem, &ncfn); + check_flag(flag, "ARKStepGetNumNonlinSolvConvFails"); + + flag = ARKStepGetNumLinSolvSetups(arkode_mem, &nsetups); + check_flag(flag, "ARKStepGetNumLinSolvSetups"); + + flag = ARKStepGetNumJacEvals(arkode_mem, &nje); + check_flag(flag, "ARKStepGetNumJacEvals"); + + flag = ARKStepGetNumLinRhsEvals(arkode_mem, &nfeLS); + check_flag(flag, "ARKStepGetNumLinRhsEvals"); + + printf(" Total number of Newton iterations = %li\n", nni); + printf(" Total number of linear solver convergence failures = %li\n", ncfn); + printf(" Total linear solver setups = %li\n", nsetups); + printf(" Total number of Jacobian evaluations = %li\n", nje); + printf(" Total RHS evals for setting up the linear system = %li\n", nfeLS); + } + + if (relax) + { + flag = ARKStepGetNumRelaxFnEvals(arkode_mem, &nre); + check_flag(flag, "ARKStepGetNumRelaxFnEvals"); + + flag = ARKStepGetNumRelaxJacEvals(arkode_mem, &nrje); + check_flag(flag, "ARKStepGetNumRelaxJacEvals"); + + flag = ARKStepGetNumRelaxFails(arkode_mem, &nrf); + check_flag(flag, "ARKStepGetNumRelaxFails"); + + flag = ARKStepGetNumRelaxBoundFails(arkode_mem, &nrbf); + check_flag(flag, "ARKStepGetNumRelaxBoundFails"); + + flag = ARKStepGetNumRelaxSolveFails(arkode_mem, &nrnlsf); + check_flag(flag, "ARKStepGetNumRelaxSolveFails"); + + flag = ARKStepGetNumRelaxSolveIters(arkode_mem, &nrnlsi); + check_flag(flag, "ARKStepGetNumRelaxSolveIters"); + + printf(" Total Relaxation Fn evals = %li\n", nre); + printf(" Total Relaxation Jac evals = %li\n", nrje); + printf(" Total Relaxation fails = %li\n", nrf); + printf(" Total Relaxation bound fails = %li\n", nrbf); + printf(" Total Relaxation NLS fails = %li\n", nrnlsf); + printf(" Total Relaxation NLS iters = %li\n", nrnlsi); + } + printf("\n"); + + /* -------- * + * Clean up * + * -------- */ + + /* Free ARKStep integrator and SUNDIALS objects */ + ARKStepFree(&arkode_mem); + SUNLinSolFree(LS); + SUNMatDestroy(A); + N_VDestroy(y); + N_VDestroy(ytrue); + SUNContext_Free(&ctx); + + return flag; +} + +/* ----------------------- * + * User-supplied functions * + * ----------------------- */ + +/* ODE RHS function f(t,y). */ +int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* fdata = N_VGetArrayPointer(ydot); + fdata[0] = -EXP(ydata[0]); + return 0; +} + +/* ODE RHS Jacobian function J(t,y) = df/dy. */ +int Jac(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, + N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* Jdata = SUNDenseMatrix_Data(J); + Jdata[0] = -EXP(ydata[0]); + return 0; +} + +/* Entropy function e(y) */ +int Ent(N_Vector y, sunrealtype* e, void* user_data) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + *e = EXP(ydata[0]); + return 0; +} + +/* Entropy function Jacobian Je(y) = de/dy */ +int JacEnt(N_Vector y, N_Vector J, void* user_data) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* jdata = N_VGetArrayPointer(J); + jdata[0] = EXP(ydata[0]); + return 0; +} + +/* ----------------- * + * Utility functions * + * ----------------- */ + +/* Analytic solution */ +int ans(sunrealtype t, N_Vector y) +{ + sunrealtype* ydata = N_VGetArrayPointer(y); + ydata[0] = LOG(EXP(SUN_RCONST(-0.5)) + t); + return 0; +} + +/* Check for an unrecoverable (negative) return flag from a SUNDIALS function */ +int check_flag(int flag, const char* funcname) +{ + if (flag < 0) + { + fprintf(stderr, "ERROR: %s() returned %d\n", funcname, flag); + return 1; + } + return 0; +} + +/* Check if a function returned a NULL pointer */ +int check_ptr(void* ptr, const char* funcname) +{ + if (!ptr) + { + fprintf(stderr, "ERROR: %s() returned NULL\n", funcname); + return 1; + } + return 0; +} diff --git a/examples/arkode/C_serial/ark_dissipated_exp_entropy_1_0.out b/examples/arkode/C_serial/ark_dissipated_exp_entropy_1_0.out new file mode 100644 index 0000000000..90df8380be --- /dev/null +++ b/examples/arkode/C_serial/ark_dissipated_exp_entropy_1_0.out @@ -0,0 +1,26 @@ + +Dissipated Exponential Entropy problem: + method = ERK + reltol = 1.0e-06 + abstol = 1.0e-10 + relaxation = ON + + step t u e u_err delta e + ------------------------------------------------------------------------------- + 0 0.000000e+00 1.000000e+00 2.718282e+00 0.000000e+00 0.000000e+00 + 40 4.675332e-01 1.798296e-01 1.197013e+00 1.083801e-01 -1.521269e+00 + 80 1.147092e+00 -4.153968e-01 6.600783e-01 -9.770808e-01 -2.058204e+00 + 120 3.654055e+00 -1.391763e+00 2.486366e-01 -2.841170e+00 -2.469645e+00 + ------------------------------------------------------------------------------- + +Final Solver Statistics: + Internal solver steps = 130 (attempted = 130) + Total number of error test failures = 0 + Total RHS evals: Fe = 783, Fi = 0 + Total Relaxation Fn evals = 392 + Total Relaxation Jac evals = 782 + Total Relaxation fails = 0 + Total Relaxation bound fails = 0 + Total Relaxation NLS fails = 0 + Total Relaxation NLS iters = 132 + diff --git a/examples/arkode/C_serial/ark_dissipated_exp_entropy_1_1.out b/examples/arkode/C_serial/ark_dissipated_exp_entropy_1_1.out new file mode 100644 index 0000000000..49be4a5967 --- /dev/null +++ b/examples/arkode/C_serial/ark_dissipated_exp_entropy_1_1.out @@ -0,0 +1,59 @@ + +Dissipated Exponential Entropy problem: + method = DIRK + reltol = 1.0e-06 + abstol = 1.0e-10 + relaxation = ON + + step t u e u_err delta e + ------------------------------------------------------------------------------- + 0 0.000000e+00 1.000000e+00 2.718282e+00 0.000000e+00 0.000000e+00 + 40 5.532834e-02 8.598921e-01 2.362906e+00 1.272595e+00 -3.553760e-01 + 80 1.097061e-01 7.390122e-01 2.093866e+00 1.072757e+00 -6.244157e-01 + 120 1.621192e-01 6.348810e-01 1.886798e+00 8.980007e-01 -8.314842e-01 + 160 2.128767e-01 5.434247e-01 1.721894e+00 7.425986e-01 -9.963881e-01 + 200 2.816359e-01 4.315291e-01 1.539610e+00 5.501251e-01 -1.178672e+00 + 240 3.484648e-01 3.335947e-01 1.395977e+00 3.796434e-01 -1.322305e+00 + 280 4.133160e-01 2.469303e-01 1.280090e+00 2.272780e-01 -1.438192e+00 + 320 4.764552e-01 1.692067e-01 1.184365e+00 8.948479e-02 -1.533917e+00 + 360 5.342836e-01 1.029604e-01 1.108447e+00 -2.878187e-02 -1.609834e+00 + 400 5.805889e-01 5.290718e-02 1.054332e+00 -1.186227e-01 -1.663950e+00 + 440 6.131438e-01 1.915938e-02 1.019344e+00 -1.794246e-01 -1.698938e+00 + 480 6.301872e-01 1.935497e-03 1.001937e+00 -2.105255e-01 -1.716344e+00 + 520 6.350100e-01 -2.884950e-03 9.971192e-01 -2.192380e-01 -1.721163e+00 + 560 6.468882e-01 -1.465940e-02 9.854475e-01 -2.405343e-01 -1.732834e+00 + 600 6.693569e-01 -3.655950e-02 9.641007e-01 -2.802015e-01 -1.754181e+00 + 640 7.021244e-01 -6.766195e-02 9.345763e-01 -3.366619e-01 -1.783705e+00 + 680 7.450858e-01 -1.070275e-01 8.985009e-01 -4.083288e-01 -1.819781e+00 + 720 8.015384e-01 -1.565058e-01 8.551265e-01 -4.987252e-01 -1.863155e+00 + 760 8.754911e-01 -2.178256e-01 8.042657e-01 -6.112329e-01 -1.914016e+00 + 800 9.627556e-01 -2.856560e-01 7.515211e-01 -7.362769e-01 -1.966761e+00 + 840 1.068726e+00 -3.622828e-01 6.960855e-01 -8.782493e-01 -2.022196e+00 + 880 1.205356e+00 -4.531342e-01 6.356328e-01 -1.047503e+00 -2.082649e+00 + 920 1.360691e+00 -5.472947e-01 5.785127e-01 -1.223917e+00 -2.139769e+00 + 960 1.555080e+00 -6.538651e-01 5.200319e-01 -1.424719e+00 -2.198250e+00 + 1000 1.799105e+00 -7.733362e-01 4.614709e-01 -1.651150e+00 -2.256811e+00 + 1040 2.083363e+00 -8.965948e-01 4.079565e-01 -1.886097e+00 -2.310325e+00 + 1080 2.428941e+00 -1.028483e+00 3.575490e-01 -2.138850e+00 -2.360733e+00 + 1120 2.880022e+00 -1.178008e+00 3.078913e-01 -2.426922e+00 -2.410391e+00 + 1160 3.423853e+00 -1.332822e+00 2.637319e-01 -2.726684e+00 -2.454550e+00 + 1200 4.082008e+00 -1.492878e+00 2.247249e-01 -3.037999e+00 -2.493557e+00 + 1240 4.890808e+00 -1.659881e+00 1.901617e-01 -3.364145e+00 -2.528120e+00 + ------------------------------------------------------------------------------- + +Final Solver Statistics: + Internal solver steps = 1244 (attempted = 1246) + Total number of error test failures = 2 + Total RHS evals: Fe = 0, Fi = 9633 + Total number of Newton iterations = 5892 + Total number of linear solver convergence failures = 81 + Total linear solver setups = 149 + Total number of Jacobian evaluations = 93 + Total RHS evals for setting up the linear system = 0 + Total Relaxation Fn evals = 3773 + Total Relaxation Jac evals = 5019 + Total Relaxation fails = 0 + Total Relaxation bound fails = 0 + Total Relaxation NLS fails = 0 + Total Relaxation NLS iters = 1281 + diff --git a/examples/templates/cmakelists_serial_C_ex.in b/examples/templates/cmakelists_serial_C_ex.in index 804cda9968..33c620e643 100644 --- a/examples/templates/cmakelists_serial_C_ex.in +++ b/examples/templates/cmakelists_serial_C_ex.in @@ -83,6 +83,9 @@ set(SUNDIALS_LIBRARIES ${SUNDIALS_MANYVEC_LIB} ${SUNDIALS_EXTRA_LIBS}) +# Include the current directory +include_directories(.) + # ------------------------------------------------------------------------------ # Set the names of the examples to be built and their dependencies diff --git a/examples/utilities/example_utilities.hpp b/examples/utilities/example_utilities.hpp index a2f4d2ffc1..2dcc7b5364 100644 --- a/examples/utilities/example_utilities.hpp +++ b/examples/utilities/example_utilities.hpp @@ -11,14 +11,15 @@ * SPDX-License-Identifier: BSD-3-Clause * SUNDIALS Copyright End * ----------------------------------------------------------------------------- - * Utility functions for C++ examples + * This header file should *NOT* be included in user codes and exists *ONLY* to + * reduce duplicate utility functions across example programs. * ---------------------------------------------------------------------------*/ #include #include #include -// Check function return flag +// Check for an unrecoverable (negative) return value from a SUNDIALS function int check_flag(const int flag, const std::string funcname) { if (!flag) return 0; @@ -35,17 +36,30 @@ int check_ptr(const void* ptr, const std::string funcname) return 1; } -inline void find_arg(std::vector& args, const std::string key, sunrealtype& dest) +// Functions for parsing vectors of command line inputs +inline void find_arg(std::vector& args, const std::string key, float& dest) { auto it = std::find(args.begin(), args.end(), key); if (it != args.end()) { -#if defined(SUNDIALS_SINGLE_PRECISION) dest = stof(*(it + 1)); -#elif defined(SUNDIALS_DOUBLE_PRECISION) + args.erase(it, it + 2); + } +} + +inline void find_arg(std::vector& args, const std::string key, double& dest) +{ + auto it = std::find(args.begin(), args.end(), key); + if (it != args.end()) { dest = stod(*(it + 1)); -#elif defined(SUNDIALS_EXTENDED_PRECISION) + args.erase(it, it + 2); + } +} + +inline void find_arg(std::vector& args, const std::string key, long double& dest) +{ + auto it = std::find(args.begin(), args.end(), key); + if (it != args.end()) { dest = stold(*(it + 1)); -#endif args.erase(it, it + 2); } } diff --git a/examples/utilities/plot_data_time_series.py b/examples/utilities/plot_data_time_series.py new file mode 100755 index 0000000000..46c5d1e831 --- /dev/null +++ b/examples/utilities/plot_data_time_series.py @@ -0,0 +1,208 @@ +#!/usr/bin/env python3 +# ----------------------------------------------------------------------------- +# Programmer(s): David J. Gardner @ LLNL +# ----------------------------------------------------------------------------- +# SUNDIALS Copyright Start +# Copyright (c) 2002-2022, 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 +# ----------------------------------------------------------------------------- +# Script to plot and compare time series data from multiple files +# +# The input data files must contain an M x (N + 1) matrix of values where the +# first column of each row is the output time and the remaining values are the +# N output quantities. Additionally the data files must contain header comment +# starting with "vars:" that lists the name for each data column i.e., +# +# # vars: time quantity_0 quantity_1 . . . quantity_N-1 +# t_0 q_0,0 q_0,1 q_0,2 . . . q_0,N-1 +# t_1 q_1,0 q_1,1 q_1,2 . . . q_1,N-1 +# . . . . . . +# . . . . . . +# . . . . . . +# t_M-1 q_M-1,0 q_M-1,1 q_M-1,2 . . . q_M-1,N-1 +# +# where t_m is the m-th output time and q_m,n is the n-th quantity at the m-th +# output time. +# ----------------------------------------------------------------------------- + +# ----------------------------------------------------------------------------- +# main routine +# ----------------------------------------------------------------------------- +def main(): + + import argparse + import matplotlib.pyplot as plt + import numpy as np + import shlex + + parser = argparse.ArgumentParser(description='''Plot data files''') + + parser.add_argument('quantity', type=str, + help='''Quantity to plot''') + + parser.add_argument('datafiles', type=str, nargs='+', + help='''Data files to plot''') + + # Plot display options + + parser.add_argument('--save', action='store_true', + help='''Save figure to file''') + + parser.add_argument('--labels', type=str, nargs='+', + help='''Data file labels for plot legend''') + + parser.add_argument('--title', type=str, + help='''Plot title''') + + parser.add_argument('--xlabel', type=str, + help='''x-axis label''') + + parser.add_argument('--ylabel', type=str, + help='''y-axis label''') + + parser.add_argument('--grid', action='store_true', + help='''Add grid to plot''') + + # Axis scaling + logscale = parser.add_mutually_exclusive_group() + + logscale.add_argument('--logx', action='store_true', + help='''Plot with log scale x-axis''') + + logscale.add_argument('--logy', action='store_true', + help='''Plot with log scale y-axis''') + + logscale.add_argument('--loglog', action='store_true', + help='''Use log scale x and y axes''') + + # Debugging options + + parser.add_argument('--debug', action='store_true', + help='Enable debugging') + + # Parse command line args + args = parser.parse_args() + + # Create figure and axes + fig, ax = plt.subplots() + + for i, datafile in enumerate(args.datafiles): + + quantities = None + with open(datafile) as fn: + # read the file line by line + for line in fn: + # skip empty lines + if not line.strip(): + continue + # exit after reading initial comment lines + if "#" not in line: + break + # split line into a list + text = shlex.split(line) + # extract quantity names + if "vars:" in line: + quantities = text[2:] + continue + + # Check inputs + if quantities is None: + print("ERROR: quantity names not provided") + sys.exit() + + if args.quantity not in quantities: + print("ERROR: quantity not found") + print(f" Possible values: {quantities}") + sys.exit() + + # Get data column index + idx = quantities.index(args.quantity) + + # Load data + data = np.loadtxt(datafile, dtype=np.double) + + if args.debug: + print(np.shape(data)) + print(data) + + # Extract t and q data + tdata = data[:,0] # first column has t values + qdata = data[:,idx] # remaining columns have q values + + # line colors: matplotlib.org/stable/tutorials/colors/colormaps.html + # and colorbrewer2.org) + if len(args.datafiles) < 22: + colors = ["#d62728", "#1f77b4", "#2ca02c", "#9467bd", "#ff7f0e", + "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf", + "#000000", "#ff9896", "#aec7e8", "#98df8a", "#c5b0d5", + "#ffbb78", "#c49c94", "#f7b6d2", "#c7c7c7", "#dbdb8d", + "#9edae5"] + else: + print("ERROR: ncols > ncolors") + sys.exit() + + # Set plot label for legend + if (args.labels): + label=args.labels[i] + else: + label=None + + if args.logx or args.logy or args.loglog: + ax.plot(tdata, np.abs(qdata), label=label, + color=colors[i]) + else: + ax.plot(tdata, qdata, label=label, + color=colors[i]) + + # Change axis scale + if args.logx: + ax.set_xscale("log") + elif args.logy: + ax.set_yscale("log") + elif args.loglog: + ax.set_xscale("log") + ax.set_yscale("log") + + # Add title + if args.title: + plt.title(args.title) + + # Add x-axis label + if args.xlabel: + plt.xlabel(args.xlabel) + else: + plt.xlabel("time") + + # Add y-axis label + if args.ylabel: + plt.ylabel(args.ylabel) + else: + plt.ylabel(args.quantity.replace("_"," ")); + + # Add legend + if args.labels: + ax.legend(bbox_to_anchor=(1, 0.95), loc="upper right") + + # Add grid + if args.grid: + plt.grid() + + # Save plot to file + if args.save: + plt.savefig("fig.pdf") + else: + plt.show() + +# ----------------------------------------------------------------------------- +# run the main routine +# ----------------------------------------------------------------------------- + +if __name__ == '__main__': + import sys + sys.exit(main()) diff --git a/include/arkode/arkode.h b/include/arkode/arkode.h index be8bd7a543..3b12589182 100644 --- a/include/arkode/arkode.h +++ b/include/arkode/arkode.h @@ -132,6 +132,12 @@ extern "C" { #define ARK_CONTEXT_ERR -42 +#define ARK_RELAX_FAIL -43 +#define ARK_RELAX_MEM_NULL -44 +#define ARK_RELAX_FUNC_FAIL -45 +#define ARK_RELAX_JAC_FAIL -46 + + #define ARK_UNRECOGNIZED_ERROR -99 /* ------------------------------ @@ -170,12 +176,25 @@ typedef int (*ARKPostProcessFn)(realtype t, N_Vector y, typedef int (*ARKStagePredictFn)(realtype t, N_Vector zpred, void *user_data); +typedef int (*ARKRelaxFn)(N_Vector y, realtype* r, void* user_data); + +typedef int (*ARKRelaxJacFn)(N_Vector y, N_Vector J, void* user_data); + /* -------------------------- * MRIStep Inner Stepper Type * -------------------------- */ typedef _SUNDIALS_STRUCT_ _MRIStepInnerStepper *MRIStepInnerStepper; +/* -------------------------- + * Relaxation Solver Options + * -------------------------- */ + +typedef enum { + ARK_RELAX_BRENT, + ARK_RELAX_NEWTON +} ARKRelaxSolver; + #ifdef __cplusplus } #endif diff --git a/include/arkode/arkode_arkstep.h b/include/arkode/arkode_arkstep.h index d293f7611c..a4951bb888 100644 --- a/include/arkode/arkode_arkstep.h +++ b/include/arkode/arkode_arkstep.h @@ -475,6 +475,36 @@ SUNDIALS_EXPORT void ARKStepPrintMem(void* arkode_mem, FILE* outfile); SUNDIALS_EXPORT int ARKStepCreateMRIStepInnerStepper(void *arkode_mem, MRIStepInnerStepper *stepper); +/* Relaxation functions */ +SUNDIALS_EXPORT int ARKStepSetRelaxFn(void* arkode_mem, ARKRelaxFn rfn, + ARKRelaxJacFn rjac); +SUNDIALS_EXPORT int ARKStepSetRelaxEtaFail(void* arkode_mem, + sunrealtype eta_rf); +SUNDIALS_EXPORT int ARKStepSetRelaxLowerBound(void* arkode_mem, + sunrealtype lower); +SUNDIALS_EXPORT int ARKStepSetRelaxMaxFails(void* arkode_mem, int max_fails); +SUNDIALS_EXPORT int ARKStepSetRelaxMaxIters(void* arkode_mem, int max_iters); +SUNDIALS_EXPORT int ARKStepSetRelaxSolver(void* arkode_mem, + ARKRelaxSolver solver); +SUNDIALS_EXPORT int ARKStepSetRelaxResTol(void* arkode_mem, + sunrealtype res_tol); +SUNDIALS_EXPORT int ARKStepSetRelaxTol(void* arkode_mem, sunrealtype rel_tol, + sunrealtype abs_tol); +SUNDIALS_EXPORT int ARKStepSetRelaxUpperBound(void* arkode_mem, + sunrealtype upper); +SUNDIALS_EXPORT int ARKStepGetNumRelaxFnEvals(void* arkode_mem, + long int* r_evals); +SUNDIALS_EXPORT int ARKStepGetNumRelaxJacEvals(void* arkode_mem, + long int* J_evals); +SUNDIALS_EXPORT int ARKStepGetNumRelaxFails(void* arkode_mem, + long int* relax_fails); +SUNDIALS_EXPORT int ARKStepGetNumRelaxBoundFails(void* arkode_mem, + long int* fails); +SUNDIALS_EXPORT int ARKStepGetNumRelaxSolveFails(void* arkode_mem, + long int* fails); +SUNDIALS_EXPORT int ARKStepGetNumRelaxSolveIters(void* arkode_mem, + long int* iters); + #ifdef __cplusplus } #endif diff --git a/include/arkode/arkode_erkstep.h b/include/arkode/arkode_erkstep.h index a9f6f2d13e..9f5b56f072 100644 --- a/include/arkode/arkode_erkstep.h +++ b/include/arkode/arkode_erkstep.h @@ -262,6 +262,35 @@ SUNDIALS_EXPORT void ERKStepFree(void **arkode_mem); /* Output the ERKStep memory structure (useful when debugging) */ SUNDIALS_EXPORT void ERKStepPrintMem(void* arkode_mem, FILE* outfile); +/* Relaxation functions */ +SUNDIALS_EXPORT int ERKStepSetRelaxFn(void* arkode_mem, ARKRelaxFn rfn, + ARKRelaxJacFn rjac); +SUNDIALS_EXPORT int ERKStepSetRelaxEtaFail(void* arkode_mem, + sunrealtype eta_rf); +SUNDIALS_EXPORT int ERKStepSetRelaxLowerBound(void* arkode_mem, + sunrealtype lower); +SUNDIALS_EXPORT int ERKStepSetRelaxMaxFails(void* arkode_mem, int max_fails); +SUNDIALS_EXPORT int ERKStepSetRelaxMaxIters(void* arkode_mem, int max_iters); +SUNDIALS_EXPORT int ERKStepSetRelaxSolver(void* arkode_mem, + ARKRelaxSolver solver); +SUNDIALS_EXPORT int ERKStepSetRelaxResTol(void* arkode_mem, + sunrealtype res_tol); +SUNDIALS_EXPORT int ERKStepSetRelaxTol(void* arkode_mem, sunrealtype rel_tol, + sunrealtype abs_tol); +SUNDIALS_EXPORT int ERKStepSetRelaxUpperBound(void* arkode_mem, + sunrealtype upper); +SUNDIALS_EXPORT int ERKStepGetNumRelaxFnEvals(void* arkode_mem, + long int* r_evals); +SUNDIALS_EXPORT int ERKStepGetNumRelaxJacEvals(void* arkode_mem, + long int* J_evals); +SUNDIALS_EXPORT int ERKStepGetNumRelaxFails(void* arkode_mem, + long int* relax_fails); +SUNDIALS_EXPORT int ERKStepGetNumRelaxBoundFails(void* arkode_mem, + long int* fails); +SUNDIALS_EXPORT int ERKStepGetNumRelaxSolveFails(void* arkode_mem, + long int* fails); +SUNDIALS_EXPORT int ERKStepGetNumRelaxSolveIters(void* arkode_mem, + long int* iters); #ifdef __cplusplus } diff --git a/scripts/arkode b/scripts/arkode index f87b467d4a..d6ab949254 100755 --- a/scripts/arkode +++ b/scripts/arkode @@ -83,8 +83,11 @@ $tar $tarfile $distrobase/src/arkode/arkode_mristep.c $tar $tarfile $distrobase/src/arkode/arkode_mristep_impl.h $tar $tarfile $distrobase/src/arkode/arkode_mristep_io.c $tar $tarfile $distrobase/src/arkode/arkode_mristep_nls.c +$tar $tarfile $distrobase/src/arkode/arkode_relaxation.c +$tar $tarfile $distrobase/src/arkode/arkode_relaxation_impl.h $tar $tarfile $distrobase/src/arkode/arkode_root.c $tar $tarfile $distrobase/src/arkode/arkode_root_impl.h +$tar $tarfile $distrobase/src/arkode/arkode_types_impl.h $tar $tarfile $distrobase/src/arkode/xbraid/CMakeLists.txt $tar $tarfile $distrobase/src/arkode/xbraid/arkode_xbraid.c @@ -184,6 +187,14 @@ $tar $tarfile $distrobase/examples/arkode/C_serial/ark_brusselator_mri.c $tar $tarfile $distrobase/examples/arkode/C_serial/ark_brusselator_mri.out $tar $tarfile $distrobase/examples/arkode/C_serial/ark_brusselator_1D_mri.c $tar $tarfile $distrobase/examples/arkode/C_serial/ark_brusselator_1D_mri.out +$tar $tarfile $distrobase/examples/arkode/C_serial/ark_conserved_exp_entropy_ark.c +$tar $tarfile $distrobase/examples/arkode/C_serial/ark_conserved_exp_entropy_ark_1_0.out +$tar $tarfile $distrobase/examples/arkode/C_serial/ark_conserved_exp_entropy_ark_1_1.out +$tar $tarfile $distrobase/examples/arkode/C_serial/ark_conserved_exp_entropy_erk.c +$tar $tarfile $distrobase/examples/arkode/C_serial/ark_conserved_exp_entropy_erk_1.out +$tar $tarfile $distrobase/examples/arkode/C_serial/ark_dissipated_exp_entropy.c +$tar $tarfile $distrobase/examples/arkode/C_serial/ark_dissipated_exp_entropy_1_0.out +$tar $tarfile $distrobase/examples/arkode/C_serial/ark_dissipated_exp_entropy_1_1.out $tar $tarfile $distrobase/examples/arkode/C_serial/ark_heat1D.c $tar $tarfile $distrobase/examples/arkode/C_serial/ark_heat1D.out $tar $tarfile $distrobase/examples/arkode/C_serial/ark_heat1D_adapt.c @@ -284,6 +295,8 @@ $tar $tarfile $distrobase/examples/arkode/CXX_serial/ark_kpr_Mt_2_-5_0_-10.out $tar $tarfile $distrobase/examples/arkode/CXX_serial/ark_kpr_Mt_2_8_0_-10.out $tar $tarfile $distrobase/examples/arkode/CXX_serial/ark_kpr_Mt_0_4_0_-10_0_10_1.out $tar $tarfile $distrobase/examples/arkode/CXX_serial/ark_kpr_Mt_0_4_0_-10_1_10_1.out +$tar $tarfile $distrobase/examples/arkode/CXX_serial/ark_pendulum.cpp +$tar $tarfile $distrobase/examples/arkode/CXX_serial/ark_pendulum.out $tar $tarfile $distrobase/examples/arkode/CXX_serial/plot_heat2D.py $tar $tarfile $distrobase/examples/arkode/CXX_serial/plot_sol.py diff --git a/src/arkode/CMakeLists.txt b/src/arkode/CMakeLists.txt index 530ac098cd..e6b9964a4e 100644 --- a/src/arkode/CMakeLists.txt +++ b/src/arkode/CMakeLists.txt @@ -36,6 +36,7 @@ set(arkode_SOURCES arkode_mristep_io.c arkode_mristep_nls.c arkode_mristep.c + arkode_relaxation.c arkode_root.c arkode.c ) diff --git a/src/arkode/arkode.c b/src/arkode/arkode.c index 724a5651c4..13e0b4d4fa 100644 --- a/src/arkode/arkode.c +++ b/src/arkode/arkode.c @@ -94,13 +94,17 @@ ARKodeMem arkCreate(SUNContext sunctx) ark_mem->constraintsSet = SUNFALSE; ark_mem->constraints = NULL; + /* Initialize relaxation variables */ + ark_mem->relax_enabled = SUNFALSE; + ark_mem->relax_mem = NULL; + /* Initialize diagnostics reporting variables */ ark_mem->report = SUNFALSE; ark_mem->diagfp = NULL; /* Initialize lrw and liw */ ark_mem->lrw = 18; - ark_mem->liw = 39; /* fcn/data ptr, int, long int, sunindextype, booleantype */ + ark_mem->liw = 41; /* fcn/data ptr, int, long int, sunindextype, booleantype */ /* No mallocs have been done yet */ ark_mem->VabstolMallocDone = SUNFALSE; @@ -642,6 +646,7 @@ int arkEvolve(ARKodeMem ark_mem, realtype tout, N_Vector yout, booleantype inactive_roots; realtype dsm; int nflag, attempts, ncf, nef, constrfails; + int relax_fails; /* Check and process inputs */ @@ -825,6 +830,7 @@ int arkEvolve(ARKodeMem ark_mem, realtype tout, N_Vector yout, /* Looping point for step attempts */ dsm = ZERO; attempts = ncf = nef = constrfails = ark_mem->last_kflag = 0; + relax_fails = 0; nflag = FIRST_CALL; for(;;) { @@ -850,6 +856,17 @@ int arkEvolve(ARKodeMem ark_mem, realtype tout, N_Vector yout, kflag = arkCheckConvergence(ark_mem, &nflag, &ncf); if (kflag < 0) break; + /* Perform relaxation: + - computes relaxation parameter + - on success, updates ycur, h, and dsm + - on recoverable failure, updates eta and signals to retry step + - on fatal error, returns negative error flag */ + if (ark_mem->relax_enabled && (kflag == ARK_SUCCESS)) + { + kflag = arkRelax(ark_mem, &relax_fails, &dsm, &nflag); + if (kflag < 0) break; + } + /* perform constraint-handling (if selected, and if solver check passed) */ if (ark_mem->constraintsSet && (kflag == ARK_SUCCESS)) { kflag = arkCheckConstraints(ark_mem, &constrfails, &nflag); @@ -1098,6 +1115,13 @@ void arkFree(void **arkode_mem) ark_mem->root_mem = NULL; } + /* free the relaxation module */ + if (ark_mem->relax_mem) + { + (void) arkRelaxDestroy(ark_mem->relax_mem); + ark_mem->relax_mem = NULL; + } + free(*arkode_mem); *arkode_mem = NULL; } @@ -2518,6 +2542,23 @@ int arkHandleFailure(ARKodeMem ark_mem, int flag) arkProcessError(ark_mem, ARK_INVALID_TABLE, "ARKODE", "ARKODE", "ARKODE was provided an invalid method table"); break; + case ARK_RELAX_FAIL: + arkProcessError(ark_mem, ARK_RELAX_FAIL, "ARKODE", "ARKODE", + "At t = %Lg the relaxation module failed", + (long double) ark_mem->tcur); + break; + case ARK_RELAX_MEM_NULL: + arkProcessError(ark_mem, ARK_RELAX_MEM_NULL, "ARKODE", "ARKODE", + "The ARKODE relaxation module memory is NULL"); + break; + case ARK_RELAX_FUNC_FAIL: + arkProcessError(ark_mem, ARK_RELAX_FUNC_FAIL, "ARKODE", "ARKODE", + "The relaxation function failed unrecoverably"); + break; + case ARK_RELAX_JAC_FAIL: + arkProcessError(ark_mem, ARK_RELAX_JAC_FAIL, "ARKODE", "ARKODE", + "The relaxation Jacobian failed unrecoverably"); + break; default: /* This return should never happen */ arkProcessError(ark_mem, ARK_UNRECOGNIZED_ERROR, "ARKODE", "ARKODE", diff --git a/src/arkode/arkode_arkstep.c b/src/arkode/arkode_arkstep.c index 8b5bde3a1b..e1245168e0 100644 --- a/src/arkode/arkode_arkstep.c +++ b/src/arkode/arkode_arkstep.c @@ -691,6 +691,15 @@ void ARKStepFree(void **arkode_mem) ark_mem->liw -= step_mem->stages; } + /* free stage vectors */ + if (step_mem->z != NULL) { + for(j=0; jstages; j++) + arkFreeVec(ark_mem, &step_mem->z[j]); + free(step_mem->z); + step_mem->z = NULL; + ark_mem->liw -= step_mem->stages; + } + /* free the reusable arrays for fused vector interface */ if (step_mem->cvals != NULL) { free(step_mem->cvals); @@ -1139,6 +1148,14 @@ int arkStep_Init(void* arkode_mem, int init_type) return(ARK_ILL_INPUT); } + /* Relaxation is incompatible with implicit RHS deduction */ + if (ark_mem->relax_enabled && step_mem->implicit && step_mem->deduce_rhs) + { + arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKODE::ARKStep", + "arkStep_Init", "Relaxation cannot be performed when deducing implicit RHS values"); + return ARK_ILL_INPUT; + } + /* Allocate ARK RHS vector memory, update storage requirements */ /* Allocate Fe[0] ... Fe[stages-1] if needed */ if (step_mem->explicit) { @@ -1162,6 +1179,20 @@ int arkStep_Init(void* arkode_mem, int init_type) ark_mem->liw += step_mem->stages; /* pointers */ } + /* Allocate stage storage for relaxation with implicit/IMEX methods or if a + fixed mass matrix is present (since we store f(t,y) not M^{-1} f(t,y)) */ + if (ark_mem->relax_enabled && (step_mem->implicit || + step_mem->mass_type == MASS_FIXED)) + { + if (step_mem->z == NULL) + step_mem->z = (N_Vector *) calloc(step_mem->stages, sizeof(N_Vector)); + for (j = 0; j < step_mem->stages; j++) { + if (!arkAllocVec(ark_mem, ark_mem->ewt, &(step_mem->z[j]))) + return(ARK_MEM_FAIL); + } + ark_mem->liw += step_mem->stages; /* pointers */ + } + /* Allocate reusable arrays for fused vector operations */ step_mem->nfusedopvecs = 2 * step_mem->stages + 2 + step_mem->nforcing; if (step_mem->cvals == NULL) { @@ -1725,6 +1756,14 @@ int arkStep_TakeStep_Z(void* arkode_mem, realtype *dsmPtr, int *nflagPtr) } /* successful stage solve */ + + /* store stage (if necessary for relaxation) */ + if (ark_mem->relax_enabled && (step_mem->implicit || + step_mem->mass_type == MASS_FIXED)) + { + N_VScale(ONE, ark_mem->ycur, step_mem->z[is]); + } + /* store implicit RHS (value in Fi[is] is from preceding nonlinear iteration) */ if (step_mem->implicit) { @@ -2164,6 +2203,48 @@ int arkStep_CheckButcherTables(ARKodeMem ark_mem) } } + /* Check if the method is compatible with relaxation */ + if (ark_mem->relax_enabled) + { + if (step_mem->q < 2) + { + arkProcessError(ark_mem, ARK_INVALID_TABLE, "ARKODE::ARKStep", + "arkStep_CheckButcherTables", + "The Butcher table(s) must be at least second order!"); + return ARK_INVALID_TABLE; + } + + if (step_mem->explicit) + { + /* Check if all b values are positive */ + for (i = 0; i < step_mem->stages; i++) + { + if (step_mem->Be->b[i] < ZERO) + { + arkProcessError(ark_mem, ARK_INVALID_TABLE, "ARKODE::ARKStep", + "arkStep_CheckButcherTables", + "The explicit Butcher table has a negative b value!"); + return ARK_INVALID_TABLE; + } + } + } + + if (step_mem->implicit) + { + /* Check if all b values are positive */ + for (i = 0; i < step_mem->stages; i++) + { + if (step_mem->Bi->b[i] < ZERO) + { + arkProcessError(ark_mem, ARK_INVALID_TABLE, "ARKODE::ARKStep", + "arkStep_CheckButcherTables", + "The implicit Butcher table has a negative b value!"); + return ARK_INVALID_TABLE; + } + } + } + } + return(ARK_SUCCESS); } @@ -2971,6 +3052,223 @@ int arkStep_SetInnerForcing(void* arkode_mem, realtype tshift, realtype tscale, return(0); } +/* ----------------------------------------------------------------------------- + * arkStep_RelaxDeltaY + * + * Computes the RK update to yn for use in relaxation methods + * ---------------------------------------------------------------------------*/ + +int arkStep_RelaxDeltaY(ARKodeMem ark_mem, N_Vector delta_y) +{ + int i, nvec, retval; + realtype* cvals; + N_Vector* Xvecs; + ARKodeARKStepMem step_mem; + + /* Access the stepper memory structure */ + if (!(ark_mem->step_mem)) + { + arkProcessError(ark_mem, ARK_MEM_NULL, "ARKODE::ARKStep", + "arkStep_RelaxDeltaY", MSG_ARKSTEP_NO_MEM); + return ARK_MEM_NULL; + } + step_mem = (ARKodeARKStepMem)(ark_mem->step_mem); + + /* Set arrays for fused vector operation */ + cvals = step_mem->cvals; + Xvecs = step_mem->Xvecs; + + nvec = 0; + for (i = 0; i < step_mem->stages; i++) + { + /* Explicit pieces */ + if (step_mem->explicit) + { + cvals[nvec] = ark_mem->h * step_mem->Be->b[i]; + Xvecs[nvec] = step_mem->Fe[i]; + nvec++; + } + /* Implicit pieces */ + if (step_mem->implicit) + { + cvals[nvec] = ark_mem->h * step_mem->Bi->b[i]; + Xvecs[nvec] = step_mem->Fi[i]; + nvec++; + } + } + + /* Compute time step update (delta_y) */ + retval = N_VLinearCombination(nvec, cvals, Xvecs, delta_y); + if (retval) return ARK_VECTOROP_ERR; + + if (step_mem->mass_type == MASS_FIXED) + { + /* Solve to compute update M^{-1} h * sum_j Ae[i,j] Fe[j] + Ai[i,j] Fi[j] */ + retval = step_mem->msolve((void *) ark_mem, delta_y, step_mem->nlscoef); + if (retval) { return ARK_MASSSOLVE_FAIL; } + } + + return ARK_SUCCESS; +} + +/* ----------------------------------------------------------------------------- + * arkStep_RelaxDeltaE + * + * Computes the change in the relaxation functions for use in relaxation methods + * delta_e = h * sum_i b_i * + * + * With implicit and IMEX methods it is necessary to store the method stages + * (or compute the delta_e estimate along the way) to avoid inconsistencies + * between z_i, F(z_i), and J_relax(z_i) that arise from reconstructing stages + * from stored RHS values like with ERK methods. As such the take step function + * stores the stages along the way but only when there is an implicit RHS. When + * a fixed mass matrix is present the stages are also stored to avoid additional + * mass matrix solves in reconstructing the stages for an ERK method. + * + * Future: when ARKStep can exploit the structure of low storage methods, it may + * be necessary to compute the delta_e estimate along the way with explicit + * methods to avoid storing additional RHS or stage values. + * ---------------------------------------------------------------------------*/ + +int arkStep_RelaxDeltaE(ARKodeMem ark_mem, ARKRelaxJacFn relax_jac_fn, + long int* num_relax_jac_evals, + sunrealtype* delta_e_out) +{ + int i, j, nvec, retval; + sunrealtype* cvals; + N_Vector* Xvecs; + ARKodeARKStepMem step_mem; + N_Vector z_stage = ark_mem->tempv2; + N_Vector J_relax = ark_mem->tempv3; + N_Vector rhs_tmp = NULL; + sunrealtype bi = ONE; + + /* Access the stepper memory structure */ + if (!(ark_mem->step_mem)) + { + arkProcessError(ark_mem, ARK_MEM_NULL, "ARKODE::ARKStep", + "arkStep_RelaxDeltaE", MSG_ARKSTEP_NO_MEM); + return ARK_MEM_NULL; + } + step_mem = (ARKodeARKStepMem)(ark_mem->step_mem); + + /* Initialize output */ + *delta_e_out = ZERO; + + /* Set arrays for fused vector operation */ + cvals = step_mem->cvals; + Xvecs = step_mem->Xvecs; + + for (i = 0; i < step_mem->stages; i++) + { + if (step_mem->implicit || step_mem->mass_type == MASS_FIXED) + { + /* Use stored stages */ + z_stage = step_mem->z[i]; + } + else + { + /* Reconstruct explicit stages */ + nvec = 0; + + cvals[nvec] = ONE; + Xvecs[nvec] = ark_mem->yn; + nvec++; + + for (j = 0; j < i; j++) + { + cvals[nvec] = ark_mem->h * step_mem->Be->A[i][j]; + Xvecs[nvec] = step_mem->Fe[j]; + nvec++; + } + + retval = N_VLinearCombination(nvec, cvals, Xvecs, z_stage); + if (retval) return ARK_VECTOROP_ERR; + } + + /* Evaluate the Jacobian at z_i */ + retval = relax_jac_fn(z_stage, J_relax, ark_mem->user_data); + (*num_relax_jac_evals)++; + if (retval < 0) { return ARK_RELAX_JAC_FAIL; } + if (retval > 0) { return ARK_RELAX_JAC_RECV; } + + /* Reset temporary RHS alias */ + rhs_tmp = z_stage; + + /* Compute delta_e = h * sum_i b_i * */ + if (step_mem->explicit && step_mem->implicit) + { + N_VLinearSum(step_mem->Be->b[i], step_mem->Fe[i], + step_mem->Bi->b[i], step_mem->Fi[i], + rhs_tmp); + bi = ONE; + } + else if (step_mem->explicit) + { + if (step_mem->mass_type == MASS_FIXED) + { + N_VScale(ONE, step_mem->Fe[i], rhs_tmp); + } + else + { + rhs_tmp = step_mem->Fe[i]; + } + bi = step_mem->Be->b[i]; + } + else + { + if (step_mem->mass_type == MASS_FIXED) + { + N_VScale(ONE, step_mem->Fi[i], rhs_tmp); + } + else + { + rhs_tmp = step_mem->Fi[i]; + } + bi = step_mem->Bi->b[i]; + } + + if (step_mem->mass_type == MASS_FIXED) + { + retval = step_mem->msolve((void *) ark_mem, rhs_tmp, + step_mem->nlscoef); + if (retval) { return ARK_MASSSOLVE_FAIL; } + } + + /* Update estimate of relaxation function change */ + if (J_relax->ops->nvdotprodlocal && J_relax->ops->nvdotprodmultiallreduce) + { + *delta_e_out += bi * N_VDotProdLocal(J_relax, rhs_tmp); + } + else + { + *delta_e_out += bi * N_VDotProd(J_relax, rhs_tmp); + } + } + + if (J_relax->ops->nvdotprodlocal && J_relax->ops->nvdotprodmultiallreduce) + { + retval = N_VDotProdMultiAllReduce(1, J_relax, delta_e_out); + if (retval) { return ARK_VECTOROP_ERR; } + } + + *delta_e_out *= ark_mem->h; + + return ARK_SUCCESS; +} + +/* ----------------------------------------------------------------------------- + * arkStep_GetOrder + * + * Returns the method order + * ---------------------------------------------------------------------------*/ + +int arkStep_GetOrder(ARKodeMem ark_mem) +{ + ARKodeARKStepMem step_mem = (ARKodeARKStepMem)(ark_mem->step_mem); + return step_mem->q; +} + /*=============================================================== EOF ===============================================================*/ diff --git a/src/arkode/arkode_arkstep_impl.h b/src/arkode/arkode_arkstep_impl.h index 307f0f8555..91ff0e6489 100644 --- a/src/arkode/arkode_arkstep_impl.h +++ b/src/arkode/arkode_arkstep_impl.h @@ -77,6 +77,7 @@ typedef struct ARKodeARKStepMemRec { /* ARK method storage and parameters */ N_Vector *Fe; /* explicit RHS at each stage */ N_Vector *Fi; /* implicit RHS at each stage */ + N_Vector *z; /* stages (for relaxation) */ N_Vector sdata; /* old stage data in residual */ N_Vector zpred; /* predicted stage solution */ N_Vector zcor; /* stage correction */ @@ -229,6 +230,11 @@ int arkStep_MRIStepInnerFullRhs(MRIStepInnerStepper stepper, realtype t, int arkStep_MRIStepInnerReset(MRIStepInnerStepper stepper, realtype tR, N_Vector yR); +/* private functions for relaxation */ +int arkStep_RelaxDeltaY(ARKodeMem ark_mem, N_Vector delta_y); +int arkStep_RelaxDeltaE(ARKodeMem ark_mem, ARKRelaxJacFn relax_jac_fn, + long int* relax_jac_fn_evals, sunrealtype* delta_e_out); +int arkStep_GetOrder(ARKodeMem ark_mem); /*=============================================================== Reusable ARKStep Error Messages diff --git a/src/arkode/arkode_arkstep_io.c b/src/arkode/arkode_arkstep_io.c index 93c99e6014..50f9182a75 100644 --- a/src/arkode/arkode_arkstep_io.c +++ b/src/arkode/arkode_arkstep_io.c @@ -28,7 +28,6 @@ #include #include - /*=============================================================== ARKStep Optional input functions (wrappers for generic ARKODE utility routines). All are documented in arkode_io.c. @@ -263,7 +262,86 @@ int ARKStepGetLastMassFlag(void *arkode_mem, long int *flag) { char *ARKStepGetLinReturnFlagName(long int flag) { return(arkLSGetReturnFlagName(flag)); } +/* ----------------------------------------------------------------------------- + * Wrappers for the ARKODE relaxation module + * ---------------------------------------------------------------------------*/ + +int ARKStepSetRelaxFn(void* arkode_mem, ARKRelaxFn rfn, ARKRelaxJacFn rjac) +{ + return arkRelaxCreate(arkode_mem, rfn, rjac, arkStep_RelaxDeltaY, + arkStep_RelaxDeltaE, arkStep_GetOrder); +} + +int ARKStepSetRelaxEtaFail(void* arkode_mem, sunrealtype eta_rf) +{ + return arkRelaxSetEtaFail(arkode_mem, eta_rf); +} + +int ARKStepSetRelaxLowerBound(void* arkode_mem, sunrealtype lower) +{ + return arkRelaxSetLowerBound(arkode_mem, lower); +} +int ARKStepSetRelaxMaxFails(void* arkode_mem, int max_fails) +{ + return arkRelaxSetMaxFails(arkode_mem, max_fails); +} + +int ARKStepSetRelaxMaxIters(void* arkode_mem, int max_iters) +{ + return arkRelaxSetMaxIters(arkode_mem, max_iters); +} + +int ARKStepSetRelaxSolver(void* arkode_mem, ARKRelaxSolver solver) +{ + return arkRelaxSetSolver(arkode_mem, solver); +} + +int ARKStepSetRelaxResTol(void* arkode_mem, sunrealtype res_tol) +{ + return arkRelaxSetResTol(arkode_mem, res_tol); +} + +int ARKStepSetRelaxTol(void* arkode_mem, sunrealtype rel_tol, + sunrealtype abs_tol) +{ + return arkRelaxSetTol(arkode_mem, rel_tol, abs_tol); +} + +int ARKStepSetRelaxUpperBound(void* arkode_mem, sunrealtype upper) +{ + return arkRelaxSetUpperBound(arkode_mem, upper); +} + +int ARKStepGetNumRelaxFnEvals(void* arkode_mem, long int* r_evals) +{ + return arkRelaxGetNumRelaxFnEvals(arkode_mem, r_evals); +} + +int ARKStepGetNumRelaxJacEvals(void* arkode_mem, long int* J_evals) +{ + return arkRelaxGetNumRelaxJacEvals(arkode_mem, J_evals); +} + +int ARKStepGetNumRelaxFails(void* arkode_mem, long int* relax_fails) +{ + return arkRelaxGetNumRelaxFails(arkode_mem, relax_fails); +} + +int ARKStepGetNumRelaxBoundFails(void* arkode_mem, long int* fails) +{ + return arkRelaxGetNumRelaxBoundFails(arkode_mem, fails); +} + +int ARKStepGetNumRelaxSolveFails(void* arkode_mem, long int* fails) +{ + return arkRelaxGetNumRelaxSolveFails(arkode_mem, fails); +} + +int ARKStepGetNumRelaxSolveIters(void* arkode_mem, long int* iters) +{ + return arkRelaxGetNumRelaxSolveIters(arkode_mem, iters); +} /*=============================================================== ARKStep optional input functions -- stepper-specific diff --git a/src/arkode/arkode_erkstep.c b/src/arkode/arkode_erkstep.c index 0d9e699e0d..f21aa2602b 100644 --- a/src/arkode/arkode_erkstep.c +++ b/src/arkode/arkode_erkstep.c @@ -1061,6 +1061,29 @@ int erkStep_CheckButcherTable(ARKodeMem ark_mem) return(ARK_INVALID_TABLE); } + /* check if all b values are positive for relaxation */ + if (ark_mem->relax_enabled) + { + if (step_mem->q < 2) + { + arkProcessError(ark_mem, ARK_INVALID_TABLE, "ARKODE::ERKStep", + "erkStep_CheckButcherTables", + "The Butcher table must be at least second order!"); + return ARK_INVALID_TABLE; + } + + for (i = 0; i < step_mem->stages; i++) + { + if (step_mem->B->b[i] < ZERO) + { + arkProcessError(ark_mem, ARK_INVALID_TABLE, "ARKODE::ERKStep", + "erkStep_CheckButcherTable", + "The Butcher table has a negative b value!"); + return ARK_INVALID_TABLE; + } + } + } + return(ARK_SUCCESS); } @@ -1146,6 +1169,141 @@ int erkStep_ComputeSolutions(ARKodeMem ark_mem, realtype *dsmPtr) return(ARK_SUCCESS); } +/* ----------------------------------------------------------------------------- + * erkStep_RelaxDeltaY + * + * Computes the RK update to yn for use in relaxation methods + * ---------------------------------------------------------------------------*/ + +int erkStep_RelaxDeltaY(ARKodeMem ark_mem, N_Vector delta_y) +{ + int i, nvec, retval; + realtype* cvals; + N_Vector* Xvecs; + ARKodeERKStepMem step_mem; + + /* Access the stepper memory structure */ + if (!(ark_mem->step_mem)) + { + arkProcessError(ark_mem, ARK_MEM_NULL, "ARKODE::ERKStep", + "erkStep_RelaxDeltaY", MSG_ERKSTEP_NO_MEM); + return ARK_MEM_NULL; + } + step_mem = (ARKodeERKStepMem)(ark_mem->step_mem); + + /* Set arrays for fused vector operation */ + cvals = step_mem->cvals; + Xvecs = step_mem->Xvecs; + + nvec = 0; + for (i = 0; i < step_mem->stages; i++) + { + cvals[nvec] = ark_mem->h * step_mem->B->b[i]; + Xvecs[nvec] = step_mem->F[i]; + nvec++; + } + + /* Compute time step update (delta_y) */ + retval = N_VLinearCombination(nvec, cvals, Xvecs, delta_y); + if (retval) return ARK_VECTOROP_ERR; + + return ARK_SUCCESS; +} + +/* ----------------------------------------------------------------------------- + * erkStep_RelaxDeltaE + * + * Computes the change in the relaxation functions for use in relaxation methods + * delta_e = h * sum_i b_i * + * ---------------------------------------------------------------------------*/ + +int erkStep_RelaxDeltaE(ARKodeMem ark_mem, ARKRelaxJacFn relax_jac_fn, + long int* num_relax_jac_evals, sunrealtype* delta_e_out) +{ + int i, j, nvec, retval; + realtype* cvals; + N_Vector* Xvecs; + ARKodeERKStepMem step_mem; + N_Vector z_stage = ark_mem->tempv2; + N_Vector J_relax = ark_mem->tempv3; + + /* Access the stepper memory structure */ + if (!(ark_mem->step_mem)) + { + arkProcessError(ark_mem, ARK_MEM_NULL, "ARKODE::ERKStep", + "erkStep_RelaxDeltaE", MSG_ERKSTEP_NO_MEM); + return ARK_MEM_NULL; + } + step_mem = (ARKodeERKStepMem)(ark_mem->step_mem); + + /* Initialize output */ + *delta_e_out = ZERO; + + /* Set arrays for fused vector operation */ + cvals = step_mem->cvals; + Xvecs = step_mem->Xvecs; + + for (i = 0; i < step_mem->stages; i++) + { + /* Construct stages z[i] = y_n + h * sum_j Ae[i,j] Fe[j] + Ai[i,j] Fi[j] */ + nvec = 0; + + cvals[nvec] = ONE; + Xvecs[nvec] = ark_mem->yn; + nvec++; + + for (j = 0; j < i; j++) + { + cvals[nvec] = ark_mem->h * step_mem->B->A[i][j]; + Xvecs[nvec] = step_mem->F[j]; + nvec++; + } + + retval = N_VLinearCombination(nvec, cvals, Xvecs, z_stage); + if (retval) return ARK_VECTOROP_ERR; + + /* Evaluate the Jacobian at z_i */ + retval = relax_jac_fn(z_stage, J_relax, ark_mem->user_data); + (*num_relax_jac_evals)++; + if (retval < 0) { return ARK_RELAX_JAC_FAIL; } + if (retval > 0) { return ARK_RELAX_JAC_RECV; } + + /* Update estimates */ + if (J_relax->ops->nvdotprodlocal && J_relax->ops->nvdotprodmultiallreduce) + { + *delta_e_out += step_mem->B->b[i] * N_VDotProdLocal(J_relax, + step_mem->F[i]); + } + else + { + *delta_e_out += step_mem->B->b[i] * N_VDotProd(J_relax, + step_mem->F[i]); + } + } + + if (J_relax->ops->nvdotprodlocal && J_relax->ops->nvdotprodmultiallreduce) + { + retval = N_VDotProdMultiAllReduce(1, J_relax, delta_e_out); + if (retval) { return ARK_VECTOROP_ERR; } + } + + *delta_e_out *= ark_mem->h; + + return ARK_SUCCESS; +} + +/* ----------------------------------------------------------------------------- + * erkStep_GetOrder + * + * Returns the method order + * ---------------------------------------------------------------------------*/ + +int erkStep_GetOrder(ARKodeMem ark_mem) +{ + ARKodeERKStepMem step_mem = (ARKodeERKStepMem)(ark_mem->step_mem); + return step_mem->q; +} + /*=============================================================== EOF ===============================================================*/ diff --git a/src/arkode/arkode_erkstep_impl.h b/src/arkode/arkode_erkstep_impl.h index b8dec90c5d..f3684d9d39 100644 --- a/src/arkode/arkode_erkstep_impl.h +++ b/src/arkode/arkode_erkstep_impl.h @@ -83,6 +83,12 @@ int erkStep_SetButcherTable(ARKodeMem ark_mem); int erkStep_CheckButcherTable(ARKodeMem ark_mem); int erkStep_ComputeSolutions(ARKodeMem ark_mem, realtype *dsm); +/* private functions for relaxation */ +int erkStep_RelaxDeltaY(ARKodeMem ark_mem, N_Vector delta_y); +int erkStep_RelaxDeltaE(ARKodeMem ark_mem, ARKRelaxJacFn relax_jac_fn, + long int* relax_jac_fn_evals, sunrealtype* delta_e_out); +int erkStep_GetOrder(ARKodeMem ark_mem); + /*=============================================================== Reusable ERKStep Error Messages ===============================================================*/ diff --git a/src/arkode/arkode_erkstep_io.c b/src/arkode/arkode_erkstep_io.c index d6cdb2fa35..34bb7d46c2 100644 --- a/src/arkode/arkode_erkstep_io.c +++ b/src/arkode/arkode_erkstep_io.c @@ -153,6 +153,86 @@ int ERKStepGetUserData(void *arkode_mem, void** user_data) { char *ERKStepGetReturnFlagName(long int flag) { return(arkGetReturnFlagName(flag)); } +/* ----------------------------------------------------------------------------- + * Wrappers for the ARKODE relaxation module + * ---------------------------------------------------------------------------*/ + +int ERKStepSetRelaxFn(void* arkode_mem, ARKRelaxFn rfn, ARKRelaxJacFn rjac) +{ + return arkRelaxCreate(arkode_mem, rfn, rjac, erkStep_RelaxDeltaY, + erkStep_RelaxDeltaE, erkStep_GetOrder); +} + +int ERKStepSetRelaxEtaFail(void* arkode_mem, sunrealtype eta_rf) +{ + return arkRelaxSetEtaFail(arkode_mem, eta_rf); +} + +int ERKStepSetRelaxLowerBound(void* arkode_mem, sunrealtype lower) +{ + return arkRelaxSetLowerBound(arkode_mem, lower); +} + +int ERKStepSetRelaxMaxFails(void* arkode_mem, int max_fails) +{ + return arkRelaxSetMaxFails(arkode_mem, max_fails); +} + +int ERKStepSetRelaxMaxIters(void* arkode_mem, int max_iters) +{ + return arkRelaxSetMaxIters(arkode_mem, max_iters); +} + +int ERKStepSetRelaxSolver(void* arkode_mem, ARKRelaxSolver solver) +{ + return arkRelaxSetSolver(arkode_mem, solver); +} + +int ERKStepSetRelaxResTol(void* arkode_mem, sunrealtype res_tol) +{ + return arkRelaxSetResTol(arkode_mem, res_tol); +} + +int ERKStepSetRelaxTol(void* arkode_mem, sunrealtype rel_tol, + sunrealtype abs_tol) +{ + return arkRelaxSetTol(arkode_mem, rel_tol, abs_tol); +} + +int ERKStepSetRelaxUpperBound(void* arkode_mem, sunrealtype upper) +{ + return arkRelaxSetUpperBound(arkode_mem, upper); +} + +int ERKStepGetNumRelaxFnEvals(void* arkode_mem, long int* r_evals) +{ + return arkRelaxGetNumRelaxFnEvals(arkode_mem, r_evals); +} + +int ERKStepGetNumRelaxJacEvals(void* arkode_mem, long int* J_evals) +{ + return arkRelaxGetNumRelaxJacEvals(arkode_mem, J_evals); +} + +int ERKStepGetNumRelaxFails(void* arkode_mem, long int* relax_fails) +{ + return arkRelaxGetNumRelaxFails(arkode_mem, relax_fails); +} + +int ERKStepGetNumRelaxBoundFails(void* arkode_mem, long int* fails) +{ + return arkRelaxGetNumRelaxBoundFails(arkode_mem, fails); +} + +int ERKStepGetNumRelaxSolveFails(void* arkode_mem, long int* fails) +{ + return arkRelaxGetNumRelaxSolveFails(arkode_mem, fails); +} + +int ERKStepGetNumRelaxSolveIters(void* arkode_mem, long int* iters) +{ + return arkRelaxGetNumRelaxSolveIters(arkode_mem, iters); +} /*=============================================================== ERKStep optional input functions -- stepper-specific diff --git a/src/arkode/arkode_impl.h b/src/arkode/arkode_impl.h index 2b2350d64f..f280538ea3 100644 --- a/src/arkode/arkode_impl.h +++ b/src/arkode/arkode_impl.h @@ -18,14 +18,18 @@ #define _ARKODE_IMPL_H #include + #include #include #include #include -#include "arkode_adapt_impl.h" -#include "arkode_root_impl.h" #include #include + +#include "arkode_types_impl.h" +#include "arkode_adapt_impl.h" +#include "arkode_relaxation_impl.h" +#include "arkode_root_impl.h" #include "sundials_context_impl.h" #include "sundials_logger_impl.h" @@ -48,6 +52,14 @@ extern "C" { #define ARK_PROFILER ark_mem->sunctx->profiler #define ARK_LOGGER ark_mem->sunctx->logger +/*=============================================================== + MACROS + ===============================================================*/ + +/* TODO(DJG): replace with signbit when C99+ is required */ +#define DIFFERENT_SIGN(a,b) ( ( (a) < 0 && (b) > 0 ) || ( (a) > 0 && (b) < 0 ) ) +#define SAME_SIGN(a,b) ( ( (a) > 0 && (b) > 0 ) || ( (a) < 0 && (b) < 0 ) ) + /*=============================================================== ARKODE Private Constants ===============================================================*/ @@ -67,6 +79,7 @@ extern "C" { #define HALF RCONST(0.5) /* real 0.5 */ #define ONE RCONST(1.0) /* real 1.0 */ #define TWO RCONST(2.0) /* real 2.0 */ +#define THREE RCONST(3.0) /* real 3.0 */ #define FOUR RCONST(4.0) /* real 4.0 */ #define FIVE RCONST(5.0) /* real 5.0 */ @@ -272,8 +285,8 @@ typedef struct ARKodeMassMemRec { The type ARKodeMem is type pointer to struct ARKodeMemRec. This structure contains fields to keep track of problem state. ---------------------------------------------------------------*/ -typedef struct ARKodeMemRec { - +struct ARKodeMemRec +{ SUNContext sunctx; realtype uround; /* machine unit roundoff */ @@ -405,6 +418,10 @@ typedef struct ARKodeMemRec { /* Rootfinding Data */ ARKodeRootMem root_mem; /* root-finding structure */ + /* Relaxation Data */ + sunbooleantype relax_enabled; /* is relaxation enabled? */ + ARKodeRelaxMem relax_mem; /* relaxation data structure */ + /* User-supplied step solution post-processing function */ ARKPostProcessFn ProcessStep; void* ps_data; /* pointer to user_data */ @@ -423,7 +440,7 @@ typedef struct ARKodeMemRec { force_pass is true and is used by the XBraid interface to determine if a time step passed or failed the time step error test. */ -} *ARKodeMem; +}; diff --git a/src/arkode/arkode_io.c b/src/arkode/arkode_io.c index b4fda7d21c..e60b76bf98 100644 --- a/src/arkode/arkode_io.c +++ b/src/arkode/arkode_io.c @@ -1717,6 +1717,7 @@ int arkGetUserData(void *arkode_mem, void** user_data) int arkPrintAllStats(void *arkode_mem, FILE *outfile, SUNOutputFormat fmt) { + int retval; ARKodeMem ark_mem; ARKodeRootMem ark_root_mem; @@ -1771,6 +1772,13 @@ int arkPrintAllStats(void *arkode_mem, FILE *outfile, SUNOutputFormat fmt) return(ARK_ILL_INPUT); } + /* Print relaxation stats */ + if (ark_mem->relax_enabled) + { + retval = arkRelaxPrintAllStats(arkode_mem, outfile, fmt); + if (retval != ARK_SUCCESS) return(retval); + } + return(ARK_SUCCESS); } diff --git a/src/arkode/arkode_relaxation.c b/src/arkode/arkode_relaxation.c new file mode 100644 index 0000000000..cc9fb5d0b8 --- /dev/null +++ b/src/arkode/arkode_relaxation.c @@ -0,0 +1,850 @@ +/* ----------------------------------------------------------------------------- + * Programmer(s): David J. Gardner @ LLNL + * ----------------------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2002-2022, 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 + * ----------------------------------------------------------------------------- + * This is the implementation file for ARKODE's relaxation (in time) + * functionality + * + * Temporary vectors utilized in the functions below: + * tempv2 - holds delta_y, the update direction vector + * tempv3 - holds y_relax, the relaxed solution vector + * tempv4 - holds J_relax, the Jacobian of the relaxation function + * ---------------------------------------------------------------------------*/ + +#include +#include +#include +#include + +#include "arkode_impl.h" +#include "arkode_relaxation_impl.h" +#include "sundials/sundials_types.h" + +/* ============================================================================= + * Private Functions + * ===========================================================================*/ + +/* Access the ARKODE and relaxation memory structures */ +static int arkRelaxAccessMem(void* arkode_mem, const char* fname, + ARKodeMem* ark_mem, ARKodeRelaxMem* relax_mem) +{ + if (!arkode_mem) + { + arkProcessError(NULL, ARK_MEM_NULL, "ARKODE", fname, MSG_ARK_NO_MEM); + return ARK_MEM_NULL; + } + *ark_mem = (ARKodeMem)arkode_mem; + + if (!((*ark_mem)->relax_mem)) + { + arkProcessError(*ark_mem, ARK_RELAX_MEM_NULL, "ARKODE", fname, + MSG_RELAX_MEM_NULL); + return ARK_RELAX_MEM_NULL; + } + *relax_mem = (ARKodeRelaxMem)((*ark_mem)->relax_mem); + + return ARK_SUCCESS; +} + +/* Evaluates the relaxation residual function */ +static int arkRelaxResidual(sunrealtype relax_param, sunrealtype* relax_res, + ARKodeMem ark_mem) +{ + int retval; + sunrealtype e_old = ark_mem->relax_mem->e_old; + sunrealtype delta_e = ark_mem->relax_mem->delta_e; + N_Vector delta_y = ark_mem->tempv2; + N_Vector y_relax = ark_mem->tempv3; + void* user_data = ark_mem->user_data; + + /* y_relax = y_n + r * delta_y */ + N_VLinearSum(ONE, ark_mem->yn, relax_param, delta_y, y_relax); + + /* Evaluate entropy function */ + retval = ark_mem->relax_mem->relax_fn(y_relax, relax_res, user_data); + ark_mem->relax_mem->num_relax_fn_evals++; + if (retval < 0) { return ARK_RELAX_FUNC_FAIL; } + if (retval > 0) { return ARK_RELAX_FUNC_RECV; } + + /* Compute relaxation residual */ + *relax_res = *relax_res - e_old - relax_param * delta_e; + + return ARK_SUCCESS; +} + +/* Evaluates the Jacobian of the relaxation residual function */ +static int arkRelaxResidualJacobian(sunrealtype relax_param, + sunrealtype* relax_jac, ARKodeMem ark_mem) +{ + int retval; + N_Vector delta_y = ark_mem->tempv2; + N_Vector y_relax = ark_mem->tempv3; + N_Vector J_relax = ark_mem->tempv4; + sunrealtype delta_e = ark_mem->relax_mem->delta_e; + void* user_data = ark_mem->user_data; + + /* y_relax = y_n + r * delta_y */ + N_VLinearSum(ONE, ark_mem->yn, relax_param, delta_y, y_relax); + + /* Evaluate Jacobian of entropy functions */ + retval = ark_mem->relax_mem->relax_jac_fn(y_relax, J_relax, user_data); + ark_mem->relax_mem->num_relax_jac_evals++; + if (retval < 0) { return ARK_RELAX_JAC_FAIL; } + if (retval > 0) { return ARK_RELAX_JAC_RECV; } + + /* Compute relaxation residual Jacobian */ + *relax_jac = N_VDotProd(delta_y, J_relax); + *relax_jac -= delta_e; + + return ARK_SUCCESS; +} + +/* Solve the relaxation residual equation using Newton's method */ +static int arkRelaxNewtonSolve(ARKodeMem ark_mem) +{ + int i, retval; + sunrealtype tol, delta; + ARKodeRelaxMem relax_mem = ark_mem->relax_mem; + + for (i = 0; i < ark_mem->relax_mem->max_iters; i++) + { + /* Compute the current residual */ + retval = arkRelaxResidual(relax_mem->relax_param, &(relax_mem->res), + ark_mem); + if (retval) return retval; + + /* Check for convergence */ + if (SUNRabs(relax_mem->res) < relax_mem->res_tol) { return ARK_SUCCESS; } + + /* Compute Jacobian */ + retval = arkRelaxResidualJacobian(relax_mem->relax_param, &(relax_mem->jac), + ark_mem); + if (retval) return retval; + + /* Update step length tolerance and solution */ + tol = (relax_mem->rel_tol * SUNRabs(relax_mem->relax_param) + + relax_mem->abs_tol); + + delta = relax_mem->res / relax_mem->jac; + relax_mem->relax_param -= delta; + + /* Update cumulative iteration count */ + relax_mem->nls_iters++; + + /* Check for small update */ + if (SUNRabs(delta) < tol) { return ARK_SUCCESS; } + } + + return ARK_RELAX_SOLVE_RECV; +} + +/* Solve the relaxation residual equation using Brent's method */ +static int arkRelaxBrentSolve(ARKodeMem ark_mem) +{ + int i, retval; + sunrealtype xa, fa; /* previous solution and function value */ + sunrealtype xb, fb; /* current solution and function value */ + sunrealtype xc, fc; /* together brac and curr bracket zero */ + sunrealtype xm; /* midpoint between brac and curr */ + sunrealtype old_update; /* previous iteration update */ + sunrealtype new_update; /* new iteration update */ + sunrealtype tol; /* iteration tolerance */ + sunrealtype pt, qt, rt, st; /* temporary values */ + + ARKodeRelaxMem relax_mem = ark_mem->relax_mem; + + /* Compute interval that brackets the root */ + xa = SUN_RCONST(0.9) * relax_mem->relax_param; + xb = SUN_RCONST(1.1) * relax_mem->relax_param; + + for (i = 0; i < 10; i++) + { + /* Compute relaxation residual */ + retval = arkRelaxResidual(xa, &fa, ark_mem); + ark_mem->relax_mem->num_relax_fn_evals++; + if (retval < 0) { return ARK_RELAX_FUNC_FAIL; } + if (retval > 0) { return ARK_RELAX_FUNC_RECV; } + + /* Check if we got lucky */ + if (SUNRabs(fa) < relax_mem->res_tol) + { + relax_mem->res = fa; + relax_mem->relax_param = xa; + return ARK_SUCCESS; + } + + if (fa < ZERO) break; + + fb = fa; + xb = xa; + xa *= SUN_RCONST(0.9); + } + if (fa > ZERO) { return ARK_RELAX_SOLVE_RECV; } + + for (i = 0; i < 10; i++) + { + /* Compute relaxation residual */ + retval = arkRelaxResidual(xb, &fb, ark_mem); + ark_mem->relax_mem->num_relax_fn_evals++; + if (retval < 0) { return ARK_RELAX_FUNC_FAIL; } + if (retval > 0) { return ARK_RELAX_FUNC_RECV; } + + /* Check if we got lucky */ + if (SUNRabs(fb) < relax_mem->res_tol) + { + relax_mem->res = fb; + relax_mem->relax_param = xb; + return ARK_SUCCESS; + } + + if (fb > ZERO) break; + + fa = fb; + xa = xb; + xb *= SUN_RCONST(1.1); + } + if (fb < ZERO) { return ARK_RELAX_SOLVE_RECV; } + + /* Initialize values bracketing values to lower bound and updates */ + xc = xa; + fc = fa; + + old_update = ZERO; + new_update = ZERO; + + /* Find root */ + for (i = 0; i < ark_mem->relax_mem->max_iters; i++) + { + /* Ensure xc and xb bracket zero */ + if (SAME_SIGN(fc,fb)) + { + xc = xa; + fc = fa; + old_update = new_update = xb - xa; + } + + /* Ensure xb is closer to zero than xc */ + if (SUNRabs(fb) > SUNRabs(fc)) + { + xa = xb; + xb = xc; + xc = xa; + + fa = fb; + fb = fc; + fc = fa; + } + + /* Update tolerance */ + tol = relax_mem->rel_tol * SUNRabs(xb) + HALF * relax_mem->abs_tol; + + /* Compute midpoint for bisection */ + xm = HALF * (xc - xb); + + /* Check for convergence */ + if (SUNRabs(xm) < tol || SUNRabs(fb) < relax_mem->res_tol) + { + relax_mem->res = fb; + relax_mem->relax_param = xb; + return ARK_SUCCESS; + } + + /* Compute iteration update */ + if (SUNRabs(old_update) >= tol && SUNRabs(fb) < SUNRabs(fa)) + { + /* Converging sufficiently fast, interpolate solution */ + st = fb / fa; + + if (xa == xc) + { + /* Two unique values available, try linear interpolant (secant) */ + pt = TWO * xm * st ; + qt = ONE - st; + } + else + { + /* Three unique values available, try inverse quadratic interpolant */ + qt = fa / fc; + rt = fb / fc; + pt = st * (TWO * xm * qt * (qt - rt) - (xb - xa) * (rt - ONE)); + qt = (qt - ONE) * (rt - ONE) * (st - ONE); + } + + /* Ensure updates produce values within [xc, xb] or [xb, xc] */ + if (pt > ZERO) { qt = -qt; } + else { pt = -pt; } + + /* Check if interpolant is acceptable, otherwise use bisection */ + st = THREE * xm * qt - SUNRabs(tol * qt); + rt = SUNRabs(old_update * qt); + + if (TWO * pt < SUNMIN(st, rt)) + { + old_update = new_update; + new_update = pt / qt; + } + else + { + new_update = xm; + old_update = xm; + } + } + else + { + /* Converging too slowly, use bisection */ + new_update = xm; + old_update = xm; + } + + /* Update solution */ + xa = xb; + fa = fb; + + /* If update is small, use tolerance in bisection direction */ + if (SUNRabs(new_update) > tol) + { + xb += new_update; + } + else + { + /* TODO(DJG): Replace with copysign when C99+ required */ + if (xm > ZERO) { xb += tol; } + else { xb -= tol; } + } + + /* Compute relaxation residual */ + retval = arkRelaxResidual(xb, &fb, ark_mem); + ark_mem->relax_mem->num_relax_fn_evals++; + if (retval < 0) { return ARK_RELAX_FUNC_FAIL; } + if (retval > 0) { return ARK_RELAX_FUNC_RECV; } + } + + return ARK_RELAX_SOLVE_RECV; +} + +/* Compute and apply relaxation parameter */ +int arkRelaxSolve(ARKodeMem ark_mem, ARKodeRelaxMem relax_mem, + sunrealtype* relax_val_out) +{ + int retval; + + /* Get the change in entropy (uses temp vectors 2 and 3) */ + retval = relax_mem->delta_e_fn(ark_mem, + relax_mem->relax_jac_fn, + &(relax_mem->num_relax_jac_evals), + &(relax_mem->delta_e)); + if (retval) return retval; + + /* Get the change in state (delta_y = tempv2) */ + retval = relax_mem->delta_y_fn(ark_mem, ark_mem->tempv2); + if (retval) return retval; + + /* Store the current relaxation function value */ + retval = relax_mem->relax_fn(ark_mem->yn, &(relax_mem->e_old), + ark_mem->user_data); + relax_mem->num_relax_fn_evals++; + if (retval < 0) { return ARK_RELAX_FUNC_FAIL; } + if (retval > 0) { return ARK_RELAX_FUNC_RECV; } + + /* Initial guess for relaxation parameter */ + relax_mem->relax_param = relax_mem->relax_param_prev; + + switch(relax_mem->solver) + { + case(ARK_RELAX_BRENT): + retval = arkRelaxBrentSolve(ark_mem); + break; + case(ARK_RELAX_NEWTON): + retval = arkRelaxNewtonSolve(ark_mem); + break; + default: + return ARK_ILL_INPUT; + break; + } + + /* Check for solver failure */ + if (retval) + { + relax_mem->nls_fails++; + return retval; + } + + /* Check for bad relaxation value */ + if (ark_mem->relax_mem->relax_param < relax_mem->lower_bound || + ark_mem->relax_mem->relax_param > relax_mem->upper_bound) + { + relax_mem->bound_fails++; + return ARK_RELAX_SOLVE_RECV; + } + + /* Save parameter for next initial guess */ + relax_mem->relax_param_prev = relax_mem->relax_param; + + /* Return relaxation value */ + *relax_val_out = ark_mem->relax_mem->relax_param; + + return ARK_SUCCESS; +} + +/* ============================================================================= + * User Functions + * ===========================================================================*/ + +/* ----------------------------------------------------------------------------- + * Set functions + * ---------------------------------------------------------------------------*/ + +int arkRelaxSetEtaFail(void* arkode_mem, sunrealtype eta_fail) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxSetEtaFail", &ark_mem, + &relax_mem); + if (retval) return retval; + + if (eta_fail > ZERO && eta_fail < ONE) { relax_mem->eta_fail = eta_fail; } + else { relax_mem->eta_fail = ARK_RELAX_DEFAULT_ETA_FAIL; } + + return ARK_SUCCESS; +} + +int arkRelaxSetLowerBound(void* arkode_mem, sunrealtype lower) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxSetLowerBound", &ark_mem, + &relax_mem); + if (retval) return retval; + + if (lower > ZERO && lower < ONE) { relax_mem->lower_bound = lower; } + else { relax_mem->lower_bound = ARK_RELAX_DEFAULT_LOWER_BOUND; } + + return ARK_SUCCESS; +} + +int arkRelaxSetMaxFails(void* arkode_mem, int max_fails) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxSetMaxFails", &ark_mem, + &relax_mem); + if (retval) return retval; + + if (max_fails > 0) { relax_mem->max_fails = max_fails; } + else { relax_mem->max_fails = ARK_RELAX_DEFAULT_MAX_FAILS; } + + return ARK_SUCCESS; +} + +int arkRelaxSetMaxIters(void* arkode_mem, int max_iters) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxSetMaxIters", &ark_mem, + &relax_mem); + if (retval) return retval; + + if (max_iters > 0) { relax_mem->max_iters = max_iters; } + else { relax_mem->max_iters = ARK_RELAX_DEFAULT_MAX_ITERS; } + + return ARK_SUCCESS; +} + +int arkRelaxSetSolver(void* arkode_mem, ARKRelaxSolver solver) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxSetSolver", &ark_mem, + &relax_mem); + if (retval) return retval; + + if (solver != ARK_RELAX_BRENT && solver != ARK_RELAX_NEWTON) + { + arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKODE", "arkRelaxSetSolver", + "An invalid relaxation solver option was provided."); + return ARK_ILL_INPUT; + } + + relax_mem->solver = solver; + + return ARK_SUCCESS; +} + +int arkRelaxSetResTol(void* arkode_mem, sunrealtype res_tol) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxSetResTol", &ark_mem, + &relax_mem); + if (retval) return retval; + + if (res_tol > ZERO) { relax_mem->res_tol = res_tol; } + else { relax_mem->res_tol = ARK_RELAX_DEFAULT_RES_TOL; } + + return ARK_SUCCESS; +} + +int arkRelaxSetTol(void* arkode_mem, sunrealtype rel_tol, sunrealtype abs_tol) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxSetTol", &ark_mem, + &relax_mem); + if (retval) return retval; + + if (rel_tol > ZERO) { relax_mem->rel_tol = rel_tol; } + else { relax_mem->rel_tol = ARK_RELAX_DEFAULT_REL_TOL; } + + if (abs_tol > ZERO) { relax_mem->abs_tol = abs_tol; } + else { relax_mem->abs_tol = ARK_RELAX_DEFAULT_ABS_TOL; } + + return ARK_SUCCESS; +} + +int arkRelaxSetUpperBound(void* arkode_mem, sunrealtype upper) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxSetUpperBound", &ark_mem, + &relax_mem); + if (retval) return retval; + + if (upper > ONE) { relax_mem->upper_bound = upper; } + else { relax_mem->upper_bound = ARK_RELAX_DEFAULT_UPPER_BOUND; } + + return ARK_SUCCESS; +} + +/* ----------------------------------------------------------------------------- + * Get functions + * ---------------------------------------------------------------------------*/ + +int arkRelaxGetNumRelaxFnEvals(void* arkode_mem, long int* r_evals) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxGetNumRelaxFnEvals", &ark_mem, + &relax_mem); + if (retval) return retval; + + *r_evals = relax_mem->num_relax_fn_evals; + + return ARK_SUCCESS; +} + +int arkRelaxGetNumRelaxJacEvals(void* arkode_mem, long int* J_evals) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxGetNumRelaxJacEvals", &ark_mem, + &relax_mem); + if (retval) return retval; + + *J_evals = relax_mem->num_relax_jac_evals; + + return ARK_SUCCESS; +} + +int arkRelaxGetNumRelaxFails(void* arkode_mem, long int* relax_fails) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxGetNumRelaxFails", &ark_mem, + &relax_mem); + if (retval) return retval; + + *relax_fails = relax_mem->num_fails; + + return ARK_SUCCESS; +} + +int arkRelaxGetNumRelaxSolveFails(void* arkode_mem, long int* fails) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxGetNumRelaxSolveFails", + &ark_mem, &relax_mem); + if (retval) return retval; + + *fails = relax_mem->nls_fails; + + return ARK_SUCCESS; +} + +int arkRelaxGetNumRelaxBoundFails(void* arkode_mem, long int* fails) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxGetNumRelaxBoundFails", + &ark_mem, &relax_mem); + if (retval) return retval; + + *fails = relax_mem->bound_fails; + + return ARK_SUCCESS; +} + +int arkRelaxGetNumRelaxSolveIters(void* arkode_mem, long int* iters) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxGetNumRelaxSolveIters", + &ark_mem, &relax_mem); + if (retval) return retval; + + *iters = relax_mem->nls_iters; + + return ARK_SUCCESS; +} + +int arkRelaxPrintAllStats(void* arkode_mem, FILE* outfile, SUNOutputFormat fmt) +{ + int retval; + ARKodeMem ark_mem; + ARKodeRelaxMem relax_mem; + + retval = arkRelaxAccessMem(arkode_mem, "arkRelaxPrintAllStats", &ark_mem, + &relax_mem); + if (retval) return retval; + + switch(fmt) + { + case SUN_OUTPUTFORMAT_TABLE: + fprintf(outfile, "Relax fn evals = %ld\n", + relax_mem->num_relax_fn_evals); + fprintf(outfile, "Relax Jac evals = %ld\n", + relax_mem->num_relax_jac_evals); + fprintf(outfile, "Relax fails = %ld\n", + relax_mem->num_fails); + fprintf(outfile, "Relax bound fails = %ld\n", + relax_mem->bound_fails); + fprintf(outfile, "Relax NLS iters = %ld\n", + relax_mem->nls_iters); + fprintf(outfile, "Relax NLS fails = %ld\n", + relax_mem->nls_fails); + break; + case SUN_OUTPUTFORMAT_CSV: + fprintf(outfile, ",Relax fn evals,%ld", relax_mem->num_relax_fn_evals); + fprintf(outfile, ",Relax Jac evals,%ld", relax_mem->num_relax_jac_evals); + fprintf(outfile, ",Relax fails,%ld", relax_mem->num_fails); + fprintf(outfile, ",Relax bound fails,%ld", relax_mem->bound_fails); + fprintf(outfile, ",Relax NLS iters,%ld", relax_mem->nls_iters); + fprintf(outfile, ",Relax NLS fails,%ld", relax_mem->nls_fails); + break; + default: + arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKODE", "arkRelaxPrintAllStats", + "Invalid formatting option."); + return ARK_ILL_INPUT; + } + + return ARK_SUCCESS; +} + +/* ============================================================================= + * Driver and Stepper Functions + * ===========================================================================*/ + +/* Constructor called by stepper */ +int arkRelaxCreate(void* arkode_mem, ARKRelaxFn relax_fn, + ARKRelaxJacFn relax_jac_fn, ARKRelaxDeltaYFn delta_y_fn, + ARKRelaxDeltaEFn delta_e_fn, ARKRelaxGetOrderFn get_order_fn) +{ + ARKodeMem ark_mem; + + /* Check inputs */ + if (!arkode_mem) + { + arkProcessError(NULL, ARK_MEM_NULL, "ARKODE", "arkRelaxCreate", + MSG_ARK_NO_MEM); + return ARK_MEM_NULL; + } + ark_mem = (ARKodeMem)arkode_mem; + + /* Disable relaxation if both user inputs are NULL */ + if (!relax_fn && !relax_jac_fn) + { + ark_mem->relax_enabled = SUNFALSE; + return ARK_SUCCESS; + } + + /* Ensure both the relaxation function and Jacobian are provided */ + if (!relax_fn) + { + arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKODE", "arkRelaxCreate", + "The relaxation function is NULL."); + return ARK_ILL_INPUT; + } + + if (!relax_jac_fn) + { + arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKODE", "arkRelaxCreate", + "The relaxation Jacobian function is NULL."); + return ARK_ILL_INPUT; + } + + /* Ensure stepper supplied inputs are provided */ + if (!delta_y_fn || !delta_e_fn || !get_order_fn) + { + arkProcessError(ark_mem, ARK_ILL_INPUT, "ARKODE", "arkRelaxCreate", + "The Delta y, Delta e, or get order function is NULL."); + return ARK_ILL_INPUT; + } + + /* Allocate and initialize relaxation memory structure */ + if (!(ark_mem->relax_mem)) + { + ark_mem->relax_mem = (ARKodeRelaxMem)malloc(sizeof(*(ark_mem->relax_mem))); + if (!(ark_mem->relax_mem)) return ARK_MEM_FAIL; + memset(ark_mem->relax_mem, 0, sizeof(struct ARKodeRelaxMemRec)); + + /* Set defaults */ + ark_mem->relax_mem->max_fails = ARK_RELAX_DEFAULT_MAX_FAILS; + ark_mem->relax_mem->lower_bound = ARK_RELAX_DEFAULT_LOWER_BOUND; + ark_mem->relax_mem->upper_bound = ARK_RELAX_DEFAULT_UPPER_BOUND; + ark_mem->relax_mem->eta_fail = ARK_RELAX_DEFAULT_ETA_FAIL; + ark_mem->relax_mem->solver = ARK_RELAX_NEWTON; + ark_mem->relax_mem->res_tol = ARK_RELAX_DEFAULT_RES_TOL; + ark_mem->relax_mem->rel_tol = ARK_RELAX_DEFAULT_REL_TOL; + ark_mem->relax_mem->abs_tol = ARK_RELAX_DEFAULT_ABS_TOL; + ark_mem->relax_mem->max_iters = ARK_RELAX_DEFAULT_MAX_ITERS; + + /* Initialize values */ + ark_mem->relax_mem->relax_param_prev = ONE; + + /* Update workspace sizes */ + ark_mem->lrw += 12; + ark_mem->liw += 14; + } + + /* Set function pointers */ + ark_mem->relax_mem->relax_fn = relax_fn; + ark_mem->relax_mem->relax_jac_fn = relax_jac_fn; + ark_mem->relax_mem->delta_y_fn = delta_y_fn; + ark_mem->relax_mem->delta_e_fn = delta_e_fn; + ark_mem->relax_mem->get_order_fn = get_order_fn; + + /* Enable relaxation */ + ark_mem->relax_enabled = SUNTRUE; + + return ARK_SUCCESS; +} + +/* Destructor called by driver */ +int arkRelaxDestroy(ARKodeRelaxMem relax_mem) +{ + if (!relax_mem) return ARK_SUCCESS; + + /* Free structure */ + free(relax_mem); + + return ARK_SUCCESS; +} + +/* Compute and apply relaxation, called by driver */ +int arkRelax(ARKodeMem ark_mem, int* relax_fails, realtype* dsm_inout, + int* nflag_out) +{ + int retval; + sunrealtype relax_val; + ARKodeRelaxMem relax_mem = ark_mem->relax_mem; + + /* Get the relaxation memory structure */ + if (!relax_mem) + { + arkProcessError(ark_mem, ARK_RELAX_MEM_NULL, "ARKODE", + "arkRelax", MSG_RELAX_MEM_NULL); + return ARK_RELAX_MEM_NULL; + } + + /* Compute the relaxation parameter */ + retval = arkRelaxSolve(ark_mem, relax_mem, &relax_val); + if (retval < 0) return retval; + if (retval > 0) + { + /* Update failure counts */ + relax_mem->num_fails++; + (*relax_fails)++; + + /* Check for max fails in a step */ + if (*relax_fails == relax_mem->max_fails) { return ARK_RELAX_FAIL; } + + /* Return with an error if |h| == hmin */ + if (SUNRabs(ark_mem->h) <= ark_mem->hmin * ONEPSM) + { + return ARK_RELAX_FAIL; + } + + /* Return with error if using fixed step sizes */ + if (ark_mem->fixedstep) { return(ARK_RELAX_FAIL); } + + /* Cut step size and try again */ + ark_mem->eta = relax_mem->eta_fail; + +#if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_INFO + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_INFO, + "ARKODE::arkStep_TakeStep_Z", "relaxation", + "relaxation failed"); +#endif + + return TRY_AGAIN; + } + + /* Update step size and error estimate */ + ark_mem->h *= relax_val; + *dsm_inout *= SUNRpowerI(relax_val, relax_mem->get_order_fn(ark_mem)); + + /* Relax solution */ + N_VLinearSum(relax_val, ark_mem->ycur, (ONE - relax_val), ark_mem->yn, + ark_mem->ycur); + +#if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_INFO + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_INFO, + "ARKODE::arkStep_TakeStep_Z", "relaxation", + "relaxation parameter = %"RSYM", relaxed h = %"RSYM + ", relaxed error = %"RSYM, + relax_val, ark_mem->h, *dsm_inout); +#endif + + return ARK_SUCCESS; +} + +/* ============================================================================= + * EOF + * ===========================================================================*/ diff --git a/src/arkode/arkode_relaxation_impl.h b/src/arkode/arkode_relaxation_impl.h new file mode 100644 index 0000000000..9a6e5b4c38 --- /dev/null +++ b/src/arkode/arkode_relaxation_impl.h @@ -0,0 +1,140 @@ +/* ----------------------------------------------------------------------------- + * Programmer(s): David J. Gardner @ LLNL + * ----------------------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2002-2022, 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 + * ----------------------------------------------------------------------------- + * Implementation header file for ARKODE's relaxation (in time) functionality. + * ---------------------------------------------------------------------------*/ + +#ifndef _ARKODE_RELAX_IMPL_H +#define _ARKODE_RELAX_IMPL_H + +#include + +#include +#include +#include + +#include "arkode_types_impl.h" + +/* ----------------------------------------------------------------------------- + * Relaxation Constants + * ---------------------------------------------------------------------------*/ + +#define ARK_RELAX_DEFAULT_MAX_FAILS 10 +#define ARK_RELAX_DEFAULT_RES_TOL (4 * SUN_UNIT_ROUNDOFF) +#define ARK_RELAX_DEFAULT_REL_TOL (4 * SUN_UNIT_ROUNDOFF) +#define ARK_RELAX_DEFAULT_ABS_TOL SUN_RCONST(1.0e-14) +#define ARK_RELAX_DEFAULT_MAX_ITERS 10 +#define ARK_RELAX_DEFAULT_LOWER_BOUND SUN_RCONST(0.8) +#define ARK_RELAX_DEFAULT_UPPER_BOUND SUN_RCONST(1.2) +#define ARK_RELAX_DEFAULT_ETA_FAIL SUN_RCONST(0.25) + +/* ----------------------------------------------------------------------------- + * Relaxation Private Return Values (see arkode/arkode.h for public values) + * ---------------------------------------------------------------------------*/ + +#define ARK_RELAX_FUNC_RECV 1; +#define ARK_RELAX_JAC_RECV 2; +#define ARK_RELAX_SOLVE_RECV 3; + +/* ----------------------------------------------------------------------------- + * Stepper Supplied Relaxation Functions + * ---------------------------------------------------------------------------*/ + +/* Compute the change in state for the current step y_new = y_old + delta_y */ +typedef int (*ARKRelaxDeltaYFn)(ARKodeMem ark_mem, N_Vector delta_y); + +/* Compute the estimated change in entropy for this step delta_e */ +typedef int (*ARKRelaxDeltaEFn)(ARKodeMem ark_mem, ARKRelaxJacFn relax_jac_fn, + long int* evals_out, sunrealtype* delta_e_out); + +/* Get the method order */ +typedef int (*ARKRelaxGetOrderFn)(ARKodeMem ark_mem); + +/* ----------------------------------------------------------------------------- + * Relaxation Data Structure + * ---------------------------------------------------------------------------*/ + +struct ARKodeRelaxMemRec +{ + /* user-supplied and stepper supplied functions */ + ARKRelaxFn relax_fn; /* user relaxation function ("entropy") */ + ARKRelaxJacFn relax_jac_fn; /* user relaxation Jacobian */ + ARKRelaxDeltaYFn delta_y_fn; /* get delta y from stepper */ + ARKRelaxDeltaEFn delta_e_fn; /* get delta entropy from stepper */ + ARKRelaxGetOrderFn get_order_fn; /* get the method order */ + + /* relaxation variables */ + int max_fails; /* max allowed relax fails in a step */ + long int num_relax_fn_evals; /* counter for total function evals */ + long int num_relax_jac_evals; /* counter for total jacobian evals */ + long int num_fails; /* counter for total relaxation fails */ + sunrealtype e_old; /* entropy at start of step y(t_{n-1}) */ + sunrealtype delta_e; /* change in entropy */ + sunrealtype res; /* relaxation residual value */ + sunrealtype jac; /* relaxation Jacobian value */ + sunrealtype relax_param; /* current relaxation parameter value */ + sunrealtype relax_param_prev; /* previous relaxation parameter value */ + sunrealtype lower_bound; /* smallest allowed relaxation value */ + sunrealtype upper_bound; /* largest allowed relaxation value */ + sunrealtype eta_fail; /* failed relaxation step size factor */ + + /* nonlinear solver settings */ + ARKRelaxSolver solver; /* choice of relaxation solver */ + sunrealtype res_tol; /* nonlinear residual solve tolerance */ + sunrealtype rel_tol; /* nonlinear iterate relative tolerance */ + sunrealtype abs_tol; /* nonlinear iterate absolute tolerance */ + int max_iters; /* nonlinear solve max iterations */ + long int nls_iters; /* total nonlinear iterations */ + long int nls_fails; /* number of nonlinear solver fails */ + long int bound_fails; /* number of relax param bound fails */ +}; + +/* ----------------------------------------------------------------------------- + * Relaxation Functions + * ---------------------------------------------------------------------------*/ + +/* Driver and Stepper Functions */ +int arkRelaxCreate(void* arkode_mem, ARKRelaxFn relax_fn, + ARKRelaxJacFn relax_jac_fn, ARKRelaxDeltaYFn delta_y_fn, + ARKRelaxDeltaEFn delta_e_fn, + ARKRelaxGetOrderFn get_order_fn); +int arkRelaxDestroy(ARKodeRelaxMem relax_mem); +int arkRelax(ARKodeMem ark_mem, int* relax_fails, sunrealtype* dsm_inout, + int* nflag_out); + +/* User Functions */ +int arkRelaxSetEtaFail(void* arkode_mem, sunrealtype eta_fail); +int arkRelaxSetLowerBound(void* arkode_mem, sunrealtype lower); +int arkRelaxSetMaxFails(void* arkode_mem, int max_fails); +int arkRelaxSetMaxIters(void* arkode_mem, int max_iters); +int arkRelaxSetSolver(void* arkode_mem, ARKRelaxSolver solver); +int arkRelaxSetResTol(void* arkode_mem, sunrealtype res_tol); +int arkRelaxSetTol(void* arkode_mem, sunrealtype rel_tol, sunrealtype abs_tol); +int arkRelaxSetUpperBound(void* arkode_mem, sunrealtype upper); + +int arkRelaxGetNumRelaxFnEvals(void* arkode_mem, long int* r_evals); +int arkRelaxGetNumRelaxJacEvals(void* arkode_mem, long int* j_evals); +int arkRelaxGetNumRelaxFails(void* arkode_mem, long int* relax_fails); +int arkRelaxGetNumRelaxBoundFails(void* arkode_mem, long int* fails); +int arkRelaxGetNumRelaxSolveFails(void* arkode_mem, long int* fails); +int arkRelaxGetNumRelaxSolveIters(void* arkode_mem, long int* iters); + +int arkRelaxPrintAllStats(void* arkode_mem, FILE* outfile, SUNOutputFormat fmt); + +/* ----------------------------------------------------------------------------- + * Error Messages + * ---------------------------------------------------------------------------*/ + +#define MSG_RELAX_MEM_NULL "Relaxation memory is NULL." + +#endif diff --git a/src/arkode/arkode_root.c b/src/arkode/arkode_root.c index 6acdb578ce..7c0b49df81 100644 --- a/src/arkode/arkode_root.c +++ b/src/arkode/arkode_root.c @@ -526,8 +526,6 @@ int arkRootCheck3(void* arkode_mem) return(RTFOUND); } -#define DIFFERENT_SIGN(a,b) ( ( (a) < 0 && (b) > 0 ) || ( (a) > 0 && (b) < 0 ) ) - /*--------------------------------------------------------------- arkRootfind diff --git a/src/arkode/arkode_types_impl.h b/src/arkode/arkode_types_impl.h new file mode 100644 index 0000000000..ec9cbef30d --- /dev/null +++ b/src/arkode/arkode_types_impl.h @@ -0,0 +1,23 @@ +/* ----------------------------------------------------------------------------- + * Programmer(s): David J. Gardner @ LLNL + * ----------------------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2002-2022, 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 + * ----------------------------------------------------------------------------- + * Implementation header file with forward declarations of ARKODE types. + * ---------------------------------------------------------------------------*/ + +#ifndef _ARKODE_TYPES_IMPL_H +#define _ARKODE_TYPES_IMPL_H + +typedef struct ARKodeMemRec* ARKodeMem; +typedef struct ARKodeRelaxMemRec* ARKodeRelaxMem; + +#endif diff --git a/src/arkode/fmod/farkode_arkstep_mod.c b/src/arkode/fmod/farkode_arkstep_mod.c index 543f624102..2ae77231a3 100644 --- a/src/arkode/fmod/farkode_arkstep_mod.c +++ b/src/arkode/fmod/farkode_arkstep_mod.c @@ -1089,6 +1089,20 @@ SWIGEXPORT int _wrap_FARKStepSetMaxStep(void *farg1, double const *farg2) { } +SWIGEXPORT int _wrap_FARKStepSetInterpolateStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)ARKStepSetInterpolateStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FARKStepSetStopTime(void *farg1, double const *farg2) { int fresult ; void *arg1 = (void *) 0 ; @@ -2445,4 +2459,218 @@ SWIGEXPORT int _wrap_FARKStepCreateMRIStepInnerStepper(void *farg1, void *farg2) } +SWIGEXPORT int _wrap_FARKStepSetRelaxFn(void *farg1, ARKRelaxFn farg2, ARKRelaxJacFn farg3) { + int fresult ; + void *arg1 = (void *) 0 ; + ARKRelaxFn arg2 = (ARKRelaxFn) 0 ; + ARKRelaxJacFn arg3 = (ARKRelaxJacFn) 0 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (ARKRelaxFn)(farg2); + arg3 = (ARKRelaxJacFn)(farg3); + result = (int)ARKStepSetRelaxFn(arg1,arg2,arg3); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FARKStepSetRelaxEtaFail(void *farg1, double const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + sunrealtype arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (sunrealtype)(*farg2); + result = (int)ARKStepSetRelaxEtaFail(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FARKStepSetRelaxLowerBound(void *farg1, double const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + sunrealtype arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (sunrealtype)(*farg2); + result = (int)ARKStepSetRelaxLowerBound(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FARKStepSetRelaxMaxFails(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)ARKStepSetRelaxMaxFails(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FARKStepSetRelaxMaxIters(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)ARKStepSetRelaxMaxIters(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FARKStepSetRelaxSolver(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + ARKRelaxSolver arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (ARKRelaxSolver)(*farg2); + result = (int)ARKStepSetRelaxSolver(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FARKStepSetRelaxResTol(void *farg1, double const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + sunrealtype arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (sunrealtype)(*farg2); + result = (int)ARKStepSetRelaxResTol(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FARKStepSetRelaxTol(void *farg1, double const *farg2, double const *farg3) { + int fresult ; + void *arg1 = (void *) 0 ; + sunrealtype arg2 ; + sunrealtype arg3 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (sunrealtype)(*farg2); + arg3 = (sunrealtype)(*farg3); + result = (int)ARKStepSetRelaxTol(arg1,arg2,arg3); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FARKStepSetRelaxUpperBound(void *farg1, double const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + sunrealtype arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (sunrealtype)(*farg2); + result = (int)ARKStepSetRelaxUpperBound(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FARKStepGetNumRelaxFnEvals(void *farg1, long *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + long *arg2 = (long *) 0 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (long *)(farg2); + result = (int)ARKStepGetNumRelaxFnEvals(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FARKStepGetNumRelaxJacEvals(void *farg1, long *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + long *arg2 = (long *) 0 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (long *)(farg2); + result = (int)ARKStepGetNumRelaxJacEvals(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FARKStepGetNumRelaxFails(void *farg1, long *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + long *arg2 = (long *) 0 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (long *)(farg2); + result = (int)ARKStepGetNumRelaxFails(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FARKStepGetNumRelaxBoundFails(void *farg1, long *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + long *arg2 = (long *) 0 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (long *)(farg2); + result = (int)ARKStepGetNumRelaxBoundFails(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FARKStepGetNumRelaxSolveFails(void *farg1, long *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + long *arg2 = (long *) 0 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (long *)(farg2); + result = (int)ARKStepGetNumRelaxSolveFails(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FARKStepGetNumRelaxSolveIters(void *farg1, long *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + long *arg2 = (long *) 0 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (long *)(farg2); + result = (int)ARKStepGetNumRelaxSolveIters(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + diff --git a/src/arkode/fmod/farkode_arkstep_mod.f90 b/src/arkode/fmod/farkode_arkstep_mod.f90 index 5002a7b1c4..c52006c929 100644 --- a/src/arkode/fmod/farkode_arkstep_mod.f90 +++ b/src/arkode/fmod/farkode_arkstep_mod.f90 @@ -49,9 +49,11 @@ module farkode_arkstep_mod integer(C_INT), parameter, public :: ARKSTEP_DEFAULT_DIRK_3 = ARKODE_ARK324L2SA_DIRK_4_2_3 integer(C_INT), parameter, public :: ARKSTEP_DEFAULT_DIRK_4 = ARKODE_SDIRK_5_3_4 integer(C_INT), parameter, public :: ARKSTEP_DEFAULT_DIRK_5 = ARKODE_ARK548L2SA_DIRK_8_4_5 + integer(C_INT), parameter, public :: ARKSTEP_DEFAULT_ARK_ETABLE_2 = ARKODE_ARK2_ERK_3_1_2 integer(C_INT), parameter, public :: ARKSTEP_DEFAULT_ARK_ETABLE_3 = ARKODE_ARK324L2SA_ERK_4_2_3 integer(C_INT), parameter, public :: ARKSTEP_DEFAULT_ARK_ETABLE_4 = ARKODE_ARK436L2SA_ERK_6_3_4 integer(C_INT), parameter, public :: ARKSTEP_DEFAULT_ARK_ETABLE_5 = ARKODE_ARK548L2SA_ERK_8_4_5 + integer(C_INT), parameter, public :: ARKSTEP_DEFAULT_ARK_ITABLE_2 = ARKODE_ARK2_DIRK_3_1_2 integer(C_INT), parameter, public :: ARKSTEP_DEFAULT_ARK_ITABLE_3 = ARKODE_ARK324L2SA_DIRK_4_2_3 integer(C_INT), parameter, public :: ARKSTEP_DEFAULT_ARK_ITABLE_4 = ARKODE_ARK436L2SA_DIRK_6_3_4 integer(C_INT), parameter, public :: ARKSTEP_DEFAULT_ARK_ITABLE_5 = ARKODE_ARK548L2SA_DIRK_8_4_5 @@ -117,6 +119,7 @@ module farkode_arkstep_mod public :: FARKStepSetInitStep public :: FARKStepSetMinStep public :: FARKStepSetMaxStep + public :: FARKStepSetInterpolateStopTime public :: FARKStepSetStopTime public :: FARKStepClearStopTime public :: FARKStepSetFixedStep @@ -210,6 +213,21 @@ module farkode_arkstep_mod public :: FARKStepFree public :: FARKStepPrintMem public :: FARKStepCreateMRIStepInnerStepper + public :: FARKStepSetRelaxFn + public :: FARKStepSetRelaxEtaFail + public :: FARKStepSetRelaxLowerBound + public :: FARKStepSetRelaxMaxFails + public :: FARKStepSetRelaxMaxIters + public :: FARKStepSetRelaxSolver + public :: FARKStepSetRelaxResTol + public :: FARKStepSetRelaxTol + public :: FARKStepSetRelaxUpperBound + public :: FARKStepGetNumRelaxFnEvals + public :: FARKStepGetNumRelaxJacEvals + public :: FARKStepGetNumRelaxFails + public :: FARKStepGetNumRelaxBoundFails + public :: FARKStepGetNumRelaxSolveFails + public :: FARKStepGetNumRelaxSolveIters ! WRAPPER DECLARATIONS interface @@ -758,6 +776,15 @@ function swigc_FARKStepSetMaxStep(farg1, farg2) & integer(C_INT) :: fresult end function +function swigc_FARKStepSetInterpolateStopTime(farg1, farg2) & +bind(C, name="_wrap_FARKStepSetInterpolateStopTime") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +integer(C_INT), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + function swigc_FARKStepSetStopTime(farg1, farg2) & bind(C, name="_wrap_FARKStepSetStopTime") & result(fresult) @@ -1628,6 +1655,143 @@ function swigc_FARKStepCreateMRIStepInnerStepper(farg1, farg2) & integer(C_INT) :: fresult end function +function swigc_FARKStepSetRelaxFn(farg1, farg2, farg3) & +bind(C, name="_wrap_FARKStepSetRelaxFn") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_FUNPTR), value :: farg2 +type(C_FUNPTR), value :: farg3 +integer(C_INT) :: fresult +end function + +function swigc_FARKStepSetRelaxEtaFail(farg1, farg2) & +bind(C, name="_wrap_FARKStepSetRelaxEtaFail") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FARKStepSetRelaxLowerBound(farg1, farg2) & +bind(C, name="_wrap_FARKStepSetRelaxLowerBound") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FARKStepSetRelaxMaxFails(farg1, farg2) & +bind(C, name="_wrap_FARKStepSetRelaxMaxFails") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +integer(C_INT), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FARKStepSetRelaxMaxIters(farg1, farg2) & +bind(C, name="_wrap_FARKStepSetRelaxMaxIters") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +integer(C_INT), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FARKStepSetRelaxSolver(farg1, farg2) & +bind(C, name="_wrap_FARKStepSetRelaxSolver") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +integer(C_INT), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FARKStepSetRelaxResTol(farg1, farg2) & +bind(C, name="_wrap_FARKStepSetRelaxResTol") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FARKStepSetRelaxTol(farg1, farg2, farg3) & +bind(C, name="_wrap_FARKStepSetRelaxTol") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +real(C_DOUBLE), intent(in) :: farg3 +integer(C_INT) :: fresult +end function + +function swigc_FARKStepSetRelaxUpperBound(farg1, farg2) & +bind(C, name="_wrap_FARKStepSetRelaxUpperBound") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FARKStepGetNumRelaxFnEvals(farg1, farg2) & +bind(C, name="_wrap_FARKStepGetNumRelaxFnEvals") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FARKStepGetNumRelaxJacEvals(farg1, farg2) & +bind(C, name="_wrap_FARKStepGetNumRelaxJacEvals") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FARKStepGetNumRelaxFails(farg1, farg2) & +bind(C, name="_wrap_FARKStepGetNumRelaxFails") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FARKStepGetNumRelaxBoundFails(farg1, farg2) & +bind(C, name="_wrap_FARKStepGetNumRelaxBoundFails") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FARKStepGetNumRelaxSolveFails(farg1, farg2) & +bind(C, name="_wrap_FARKStepGetNumRelaxSolveFails") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FARKStepGetNumRelaxSolveIters(farg1, farg2) & +bind(C, name="_wrap_FARKStepGetNumRelaxSolveIters") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +integer(C_INT) :: fresult +end function + end interface @@ -2647,6 +2811,22 @@ function FARKStepSetMaxStep(arkode_mem, hmax) & swig_result = fresult end function +function FARKStepSetInterpolateStopTime(arkode_mem, interp) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_INT), intent(in) :: interp +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = arkode_mem +farg2 = interp +fresult = swigc_FARKStepSetInterpolateStopTime(farg1, farg2) +swig_result = fresult +end function + function FARKStepSetStopTime(arkode_mem, tstop) & result(swig_result) use, intrinsic :: ISO_C_BINDING @@ -4234,5 +4414,251 @@ function FARKStepCreateMRIStepInnerStepper(arkode_mem, stepper) & swig_result = fresult end function +function FARKStepSetRelaxFn(arkode_mem, rfn, rjac) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +type(C_FUNPTR), intent(in), value :: rfn +type(C_FUNPTR), intent(in), value :: rjac +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_FUNPTR) :: farg2 +type(C_FUNPTR) :: farg3 + +farg1 = arkode_mem +farg2 = rfn +farg3 = rjac +fresult = swigc_FARKStepSetRelaxFn(farg1, farg2, farg3) +swig_result = fresult +end function + +function FARKStepSetRelaxEtaFail(arkode_mem, eta_rf) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +real(C_DOUBLE), intent(in) :: eta_rf +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 + +farg1 = arkode_mem +farg2 = eta_rf +fresult = swigc_FARKStepSetRelaxEtaFail(farg1, farg2) +swig_result = fresult +end function + +function FARKStepSetRelaxLowerBound(arkode_mem, lower) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +real(C_DOUBLE), intent(in) :: lower +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 + +farg1 = arkode_mem +farg2 = lower +fresult = swigc_FARKStepSetRelaxLowerBound(farg1, farg2) +swig_result = fresult +end function + +function FARKStepSetRelaxMaxFails(arkode_mem, max_fails) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_INT), intent(in) :: max_fails +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = arkode_mem +farg2 = max_fails +fresult = swigc_FARKStepSetRelaxMaxFails(farg1, farg2) +swig_result = fresult +end function + +function FARKStepSetRelaxMaxIters(arkode_mem, max_iters) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_INT), intent(in) :: max_iters +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = arkode_mem +farg2 = max_iters +fresult = swigc_FARKStepSetRelaxMaxIters(farg1, farg2) +swig_result = fresult +end function + +function FARKStepSetRelaxSolver(arkode_mem, solver) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(ARKRelaxSolver), intent(in) :: solver +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = arkode_mem +farg2 = solver +fresult = swigc_FARKStepSetRelaxSolver(farg1, farg2) +swig_result = fresult +end function + +function FARKStepSetRelaxResTol(arkode_mem, res_tol) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +real(C_DOUBLE), intent(in) :: res_tol +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 + +farg1 = arkode_mem +farg2 = res_tol +fresult = swigc_FARKStepSetRelaxResTol(farg1, farg2) +swig_result = fresult +end function + +function FARKStepSetRelaxTol(arkode_mem, rel_tol, abs_tol) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +real(C_DOUBLE), intent(in) :: rel_tol +real(C_DOUBLE), intent(in) :: abs_tol +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 +real(C_DOUBLE) :: farg3 + +farg1 = arkode_mem +farg2 = rel_tol +farg3 = abs_tol +fresult = swigc_FARKStepSetRelaxTol(farg1, farg2, farg3) +swig_result = fresult +end function + +function FARKStepSetRelaxUpperBound(arkode_mem, upper) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +real(C_DOUBLE), intent(in) :: upper +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 + +farg1 = arkode_mem +farg2 = upper +fresult = swigc_FARKStepSetRelaxUpperBound(farg1, farg2) +swig_result = fresult +end function + +function FARKStepGetNumRelaxFnEvals(arkode_mem, r_evals) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_LONG), dimension(*), target, intent(inout) :: r_evals +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 + +farg1 = arkode_mem +farg2 = c_loc(r_evals(1)) +fresult = swigc_FARKStepGetNumRelaxFnEvals(farg1, farg2) +swig_result = fresult +end function + +function FARKStepGetNumRelaxJacEvals(arkode_mem, j_evals) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_LONG), dimension(*), target, intent(inout) :: j_evals +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 + +farg1 = arkode_mem +farg2 = c_loc(j_evals(1)) +fresult = swigc_FARKStepGetNumRelaxJacEvals(farg1, farg2) +swig_result = fresult +end function + +function FARKStepGetNumRelaxFails(arkode_mem, relax_fails) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_LONG), dimension(*), target, intent(inout) :: relax_fails +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 + +farg1 = arkode_mem +farg2 = c_loc(relax_fails(1)) +fresult = swigc_FARKStepGetNumRelaxFails(farg1, farg2) +swig_result = fresult +end function + +function FARKStepGetNumRelaxBoundFails(arkode_mem, fails) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_LONG), dimension(*), target, intent(inout) :: fails +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 + +farg1 = arkode_mem +farg2 = c_loc(fails(1)) +fresult = swigc_FARKStepGetNumRelaxBoundFails(farg1, farg2) +swig_result = fresult +end function + +function FARKStepGetNumRelaxSolveFails(arkode_mem, fails) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_LONG), dimension(*), target, intent(inout) :: fails +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 + +farg1 = arkode_mem +farg2 = c_loc(fails(1)) +fresult = swigc_FARKStepGetNumRelaxSolveFails(farg1, farg2) +swig_result = fresult +end function + +function FARKStepGetNumRelaxSolveIters(arkode_mem, iters) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_LONG), dimension(*), target, intent(inout) :: iters +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 + +farg1 = arkode_mem +farg2 = c_loc(iters(1)) +fresult = swigc_FARKStepGetNumRelaxSolveIters(farg1, farg2) +swig_result = fresult +end function + end module diff --git a/src/arkode/fmod/farkode_erkstep_mod.c b/src/arkode/fmod/farkode_erkstep_mod.c index d990659fb9..46836ba780 100644 --- a/src/arkode/fmod/farkode_erkstep_mod.c +++ b/src/arkode/fmod/farkode_erkstep_mod.c @@ -757,6 +757,20 @@ SWIGEXPORT int _wrap_FERKStepSetMaxStep(void *farg1, double const *farg2) { } +SWIGEXPORT int _wrap_FERKStepSetInterpolateStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)ERKStepSetInterpolateStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FERKStepSetStopTime(void *farg1, double const *farg2) { int fresult ; void *arg1 = (void *) 0 ; @@ -1348,4 +1362,218 @@ SWIGEXPORT void _wrap_FERKStepPrintMem(void *farg1, void *farg2) { } +SWIGEXPORT int _wrap_FERKStepSetRelaxFn(void *farg1, ARKRelaxFn farg2, ARKRelaxJacFn farg3) { + int fresult ; + void *arg1 = (void *) 0 ; + ARKRelaxFn arg2 = (ARKRelaxFn) 0 ; + ARKRelaxJacFn arg3 = (ARKRelaxJacFn) 0 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (ARKRelaxFn)(farg2); + arg3 = (ARKRelaxJacFn)(farg3); + result = (int)ERKStepSetRelaxFn(arg1,arg2,arg3); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FERKStepSetRelaxEtaFail(void *farg1, double const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + sunrealtype arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (sunrealtype)(*farg2); + result = (int)ERKStepSetRelaxEtaFail(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FERKStepSetRelaxLowerBound(void *farg1, double const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + sunrealtype arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (sunrealtype)(*farg2); + result = (int)ERKStepSetRelaxLowerBound(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FERKStepSetRelaxMaxFails(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)ERKStepSetRelaxMaxFails(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FERKStepSetRelaxMaxIters(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)ERKStepSetRelaxMaxIters(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FERKStepSetRelaxSolver(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + ARKRelaxSolver arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (ARKRelaxSolver)(*farg2); + result = (int)ERKStepSetRelaxSolver(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FERKStepSetRelaxResTol(void *farg1, double const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + sunrealtype arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (sunrealtype)(*farg2); + result = (int)ERKStepSetRelaxResTol(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FERKStepSetRelaxTol(void *farg1, double const *farg2, double const *farg3) { + int fresult ; + void *arg1 = (void *) 0 ; + sunrealtype arg2 ; + sunrealtype arg3 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (sunrealtype)(*farg2); + arg3 = (sunrealtype)(*farg3); + result = (int)ERKStepSetRelaxTol(arg1,arg2,arg3); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FERKStepSetRelaxUpperBound(void *farg1, double const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + sunrealtype arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (sunrealtype)(*farg2); + result = (int)ERKStepSetRelaxUpperBound(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FERKStepGetNumRelaxFnEvals(void *farg1, long *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + long *arg2 = (long *) 0 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (long *)(farg2); + result = (int)ERKStepGetNumRelaxFnEvals(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FERKStepGetNumRelaxJacEvals(void *farg1, long *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + long *arg2 = (long *) 0 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (long *)(farg2); + result = (int)ERKStepGetNumRelaxJacEvals(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FERKStepGetNumRelaxFails(void *farg1, long *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + long *arg2 = (long *) 0 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (long *)(farg2); + result = (int)ERKStepGetNumRelaxFails(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FERKStepGetNumRelaxBoundFails(void *farg1, long *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + long *arg2 = (long *) 0 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (long *)(farg2); + result = (int)ERKStepGetNumRelaxBoundFails(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FERKStepGetNumRelaxSolveFails(void *farg1, long *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + long *arg2 = (long *) 0 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (long *)(farg2); + result = (int)ERKStepGetNumRelaxSolveFails(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FERKStepGetNumRelaxSolveIters(void *farg1, long *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + long *arg2 = (long *) 0 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (long *)(farg2); + result = (int)ERKStepGetNumRelaxSolveIters(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + diff --git a/src/arkode/fmod/farkode_erkstep_mod.f90 b/src/arkode/fmod/farkode_erkstep_mod.f90 index 7e2b05ae23..3857ac1a7d 100644 --- a/src/arkode/fmod/farkode_erkstep_mod.f90 +++ b/src/arkode/fmod/farkode_erkstep_mod.f90 @@ -84,6 +84,7 @@ module farkode_erkstep_mod public :: FERKStepSetInitStep public :: FERKStepSetMinStep public :: FERKStepSetMaxStep + public :: FERKStepSetInterpolateStopTime public :: FERKStepSetStopTime public :: FERKStepClearStopTime public :: FERKStepSetFixedStep @@ -125,6 +126,21 @@ module farkode_erkstep_mod public :: FERKStepGetStepStats public :: FERKStepFree public :: FERKStepPrintMem + public :: FERKStepSetRelaxFn + public :: FERKStepSetRelaxEtaFail + public :: FERKStepSetRelaxLowerBound + public :: FERKStepSetRelaxMaxFails + public :: FERKStepSetRelaxMaxIters + public :: FERKStepSetRelaxSolver + public :: FERKStepSetRelaxResTol + public :: FERKStepSetRelaxTol + public :: FERKStepSetRelaxUpperBound + public :: FERKStepGetNumRelaxFnEvals + public :: FERKStepGetNumRelaxJacEvals + public :: FERKStepGetNumRelaxFails + public :: FERKStepGetNumRelaxBoundFails + public :: FERKStepGetNumRelaxSolveFails + public :: FERKStepGetNumRelaxSolveIters ! WRAPPER DECLARATIONS interface @@ -461,6 +477,15 @@ function swigc_FERKStepSetMaxStep(farg1, farg2) & integer(C_INT) :: fresult end function +function swigc_FERKStepSetInterpolateStopTime(farg1, farg2) & +bind(C, name="_wrap_FERKStepSetInterpolateStopTime") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +integer(C_INT), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + function swigc_FERKStepSetStopTime(farg1, farg2) & bind(C, name="_wrap_FERKStepSetStopTime") & result(fresult) @@ -844,6 +869,143 @@ subroutine swigc_FERKStepPrintMem(farg1, farg2) & type(C_PTR), value :: farg2 end subroutine +function swigc_FERKStepSetRelaxFn(farg1, farg2, farg3) & +bind(C, name="_wrap_FERKStepSetRelaxFn") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_FUNPTR), value :: farg2 +type(C_FUNPTR), value :: farg3 +integer(C_INT) :: fresult +end function + +function swigc_FERKStepSetRelaxEtaFail(farg1, farg2) & +bind(C, name="_wrap_FERKStepSetRelaxEtaFail") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FERKStepSetRelaxLowerBound(farg1, farg2) & +bind(C, name="_wrap_FERKStepSetRelaxLowerBound") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FERKStepSetRelaxMaxFails(farg1, farg2) & +bind(C, name="_wrap_FERKStepSetRelaxMaxFails") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +integer(C_INT), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FERKStepSetRelaxMaxIters(farg1, farg2) & +bind(C, name="_wrap_FERKStepSetRelaxMaxIters") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +integer(C_INT), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FERKStepSetRelaxSolver(farg1, farg2) & +bind(C, name="_wrap_FERKStepSetRelaxSolver") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +integer(C_INT), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FERKStepSetRelaxResTol(farg1, farg2) & +bind(C, name="_wrap_FERKStepSetRelaxResTol") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FERKStepSetRelaxTol(farg1, farg2, farg3) & +bind(C, name="_wrap_FERKStepSetRelaxTol") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +real(C_DOUBLE), intent(in) :: farg3 +integer(C_INT) :: fresult +end function + +function swigc_FERKStepSetRelaxUpperBound(farg1, farg2) & +bind(C, name="_wrap_FERKStepSetRelaxUpperBound") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FERKStepGetNumRelaxFnEvals(farg1, farg2) & +bind(C, name="_wrap_FERKStepGetNumRelaxFnEvals") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FERKStepGetNumRelaxJacEvals(farg1, farg2) & +bind(C, name="_wrap_FERKStepGetNumRelaxJacEvals") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FERKStepGetNumRelaxFails(farg1, farg2) & +bind(C, name="_wrap_FERKStepGetNumRelaxFails") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FERKStepGetNumRelaxBoundFails(farg1, farg2) & +bind(C, name="_wrap_FERKStepGetNumRelaxBoundFails") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FERKStepGetNumRelaxSolveFails(farg1, farg2) & +bind(C, name="_wrap_FERKStepGetNumRelaxSolveFails") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +integer(C_INT) :: fresult +end function + +function swigc_FERKStepGetNumRelaxSolveIters(farg1, farg2) & +bind(C, name="_wrap_FERKStepGetNumRelaxSolveIters") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +integer(C_INT) :: fresult +end function + end interface @@ -1479,6 +1641,22 @@ function FERKStepSetMaxStep(arkode_mem, hmax) & swig_result = fresult end function +function FERKStepSetInterpolateStopTime(arkode_mem, interp) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_INT), intent(in) :: interp +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = arkode_mem +farg2 = interp +fresult = swigc_FERKStepSetInterpolateStopTime(farg1, farg2) +swig_result = fresult +end function + function FERKStepSetStopTime(arkode_mem, tstop) & result(swig_result) use, intrinsic :: ISO_C_BINDING @@ -2178,5 +2356,251 @@ subroutine FERKStepPrintMem(arkode_mem, outfile) call swigc_FERKStepPrintMem(farg1, farg2) end subroutine +function FERKStepSetRelaxFn(arkode_mem, rfn, rjac) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +type(C_FUNPTR), intent(in), value :: rfn +type(C_FUNPTR), intent(in), value :: rjac +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_FUNPTR) :: farg2 +type(C_FUNPTR) :: farg3 + +farg1 = arkode_mem +farg2 = rfn +farg3 = rjac +fresult = swigc_FERKStepSetRelaxFn(farg1, farg2, farg3) +swig_result = fresult +end function + +function FERKStepSetRelaxEtaFail(arkode_mem, eta_rf) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +real(C_DOUBLE), intent(in) :: eta_rf +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 + +farg1 = arkode_mem +farg2 = eta_rf +fresult = swigc_FERKStepSetRelaxEtaFail(farg1, farg2) +swig_result = fresult +end function + +function FERKStepSetRelaxLowerBound(arkode_mem, lower) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +real(C_DOUBLE), intent(in) :: lower +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 + +farg1 = arkode_mem +farg2 = lower +fresult = swigc_FERKStepSetRelaxLowerBound(farg1, farg2) +swig_result = fresult +end function + +function FERKStepSetRelaxMaxFails(arkode_mem, max_fails) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_INT), intent(in) :: max_fails +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = arkode_mem +farg2 = max_fails +fresult = swigc_FERKStepSetRelaxMaxFails(farg1, farg2) +swig_result = fresult +end function + +function FERKStepSetRelaxMaxIters(arkode_mem, max_iters) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_INT), intent(in) :: max_iters +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = arkode_mem +farg2 = max_iters +fresult = swigc_FERKStepSetRelaxMaxIters(farg1, farg2) +swig_result = fresult +end function + +function FERKStepSetRelaxSolver(arkode_mem, solver) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(ARKRelaxSolver), intent(in) :: solver +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = arkode_mem +farg2 = solver +fresult = swigc_FERKStepSetRelaxSolver(farg1, farg2) +swig_result = fresult +end function + +function FERKStepSetRelaxResTol(arkode_mem, res_tol) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +real(C_DOUBLE), intent(in) :: res_tol +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 + +farg1 = arkode_mem +farg2 = res_tol +fresult = swigc_FERKStepSetRelaxResTol(farg1, farg2) +swig_result = fresult +end function + +function FERKStepSetRelaxTol(arkode_mem, rel_tol, abs_tol) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +real(C_DOUBLE), intent(in) :: rel_tol +real(C_DOUBLE), intent(in) :: abs_tol +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 +real(C_DOUBLE) :: farg3 + +farg1 = arkode_mem +farg2 = rel_tol +farg3 = abs_tol +fresult = swigc_FERKStepSetRelaxTol(farg1, farg2, farg3) +swig_result = fresult +end function + +function FERKStepSetRelaxUpperBound(arkode_mem, upper) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +real(C_DOUBLE), intent(in) :: upper +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 + +farg1 = arkode_mem +farg2 = upper +fresult = swigc_FERKStepSetRelaxUpperBound(farg1, farg2) +swig_result = fresult +end function + +function FERKStepGetNumRelaxFnEvals(arkode_mem, r_evals) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_LONG), dimension(*), target, intent(inout) :: r_evals +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 + +farg1 = arkode_mem +farg2 = c_loc(r_evals(1)) +fresult = swigc_FERKStepGetNumRelaxFnEvals(farg1, farg2) +swig_result = fresult +end function + +function FERKStepGetNumRelaxJacEvals(arkode_mem, j_evals) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_LONG), dimension(*), target, intent(inout) :: j_evals +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 + +farg1 = arkode_mem +farg2 = c_loc(j_evals(1)) +fresult = swigc_FERKStepGetNumRelaxJacEvals(farg1, farg2) +swig_result = fresult +end function + +function FERKStepGetNumRelaxFails(arkode_mem, relax_fails) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_LONG), dimension(*), target, intent(inout) :: relax_fails +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 + +farg1 = arkode_mem +farg2 = c_loc(relax_fails(1)) +fresult = swigc_FERKStepGetNumRelaxFails(farg1, farg2) +swig_result = fresult +end function + +function FERKStepGetNumRelaxBoundFails(arkode_mem, fails) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_LONG), dimension(*), target, intent(inout) :: fails +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 + +farg1 = arkode_mem +farg2 = c_loc(fails(1)) +fresult = swigc_FERKStepGetNumRelaxBoundFails(farg1, farg2) +swig_result = fresult +end function + +function FERKStepGetNumRelaxSolveFails(arkode_mem, fails) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_LONG), dimension(*), target, intent(inout) :: fails +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 + +farg1 = arkode_mem +farg2 = c_loc(fails(1)) +fresult = swigc_FERKStepGetNumRelaxSolveFails(farg1, farg2) +swig_result = fresult +end function + +function FERKStepGetNumRelaxSolveIters(arkode_mem, iters) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_LONG), dimension(*), target, intent(inout) :: iters +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 + +farg1 = arkode_mem +farg2 = c_loc(iters(1)) +fresult = swigc_FERKStepGetNumRelaxSolveIters(farg1, farg2) +swig_result = fresult +end function + end module diff --git a/src/arkode/fmod/farkode_mod.f90 b/src/arkode/fmod/farkode_mod.f90 index dd41dfce77..b0f26eda24 100644 --- a/src/arkode/fmod/farkode_mod.f90 +++ b/src/arkode/fmod/farkode_mod.f90 @@ -100,7 +100,18 @@ module farkode_mod integer(C_INT), parameter, public :: ARK_INTERP_FAIL = -40_C_INT integer(C_INT), parameter, public :: ARK_INVALID_TABLE = -41_C_INT integer(C_INT), parameter, public :: ARK_CONTEXT_ERR = -42_C_INT + integer(C_INT), parameter, public :: ARK_RELAX_FAIL = -43_C_INT + integer(C_INT), parameter, public :: ARK_RELAX_MEM_NULL = -44_C_INT + integer(C_INT), parameter, public :: ARK_RELAX_FUNC_FAIL = -45_C_INT + integer(C_INT), parameter, public :: ARK_RELAX_JAC_FAIL = -46_C_INT integer(C_INT), parameter, public :: ARK_UNRECOGNIZED_ERROR = -99_C_INT + ! typedef enum ARKRelaxSolver + enum, bind(c) + enumerator :: ARK_RELAX_BRENT + enumerator :: ARK_RELAX_NEWTON + end enum + integer, parameter, public :: ARKRelaxSolver = kind(ARK_RELAX_BRENT) + public :: ARK_RELAX_BRENT, ARK_RELAX_NEWTON public :: FARKBandPrecInit public :: FARKBandPrecGetWorkSpace public :: FARKBandPrecGetNumRhsEvals @@ -192,7 +203,8 @@ module farkode_mod enumerator :: ARKODE_ESDIRK437L2SA_7_3_4 enumerator :: ARKODE_ESDIRK547L2SA_7_4_5 enumerator :: ARKODE_ESDIRK547L2SA2_7_4_5 - enumerator :: ARKODE_MAX_DIRK_NUM = ARKODE_ESDIRK547L2SA2_7_4_5 + enumerator :: ARKODE_ARK2_DIRK_3_1_2 + enumerator :: ARKODE_MAX_DIRK_NUM = ARKODE_ARK2_DIRK_3_1_2 end enum integer, parameter, public :: ARKODE_DIRKTableID = kind(ARKODE_DIRK_NONE) public :: ARKODE_DIRK_NONE, ARKODE_MIN_DIRK_NUM, ARKODE_SDIRK_2_1_2, ARKODE_BILLINGTON_3_3_2, ARKODE_TRBDF2_3_3_2, & @@ -200,7 +212,8 @@ module farkode_mod ARKODE_KVAERNO_5_3_4, ARKODE_ARK436L2SA_DIRK_6_3_4, ARKODE_KVAERNO_7_4_5, ARKODE_ARK548L2SA_DIRK_8_4_5, & ARKODE_ARK437L2SA_DIRK_7_3_4, ARKODE_ARK548L2SAb_DIRK_8_4_5, ARKODE_ESDIRK324L2SA_4_2_3, ARKODE_ESDIRK325L2SA_5_2_3, & ARKODE_ESDIRK32I5L2SA_5_2_3, ARKODE_ESDIRK436L2SA_6_3_4, ARKODE_ESDIRK43I6L2SA_6_3_4, ARKODE_QESDIRK436L2SA_6_3_4, & - ARKODE_ESDIRK437L2SA_7_3_4, ARKODE_ESDIRK547L2SA_7_4_5, ARKODE_ESDIRK547L2SA2_7_4_5, ARKODE_MAX_DIRK_NUM + ARKODE_ESDIRK437L2SA_7_3_4, ARKODE_ESDIRK547L2SA_7_4_5, ARKODE_ESDIRK547L2SA2_7_4_5, ARKODE_ARK2_DIRK_3_1_2, & + ARKODE_MAX_DIRK_NUM public :: FARKodeButcherTable_LoadDIRK type, bind(C) :: SwigArrayWrapper type(C_PTR), public :: data = C_NULL_PTR @@ -243,14 +256,15 @@ module farkode_mod enumerator :: ARKODE_KNOTH_WOLKE_3_3 enumerator :: ARKODE_ARK437L2SA_ERK_7_3_4 enumerator :: ARKODE_ARK548L2SAb_ERK_8_4_5 - enumerator :: ARKODE_MAX_ERK_NUM = ARKODE_ARK548L2SAb_ERK_8_4_5 + enumerator :: ARKODE_ARK2_ERK_3_1_2 + enumerator :: ARKODE_MAX_ERK_NUM = ARKODE_ARK2_ERK_3_1_2 end enum integer, parameter, public :: ARKODE_ERKTableID = kind(ARKODE_ERK_NONE) public :: ARKODE_ERK_NONE, ARKODE_MIN_ERK_NUM, ARKODE_HEUN_EULER_2_1_2, ARKODE_BOGACKI_SHAMPINE_4_2_3, & ARKODE_ARK324L2SA_ERK_4_2_3, ARKODE_ZONNEVELD_5_3_4, ARKODE_ARK436L2SA_ERK_6_3_4, ARKODE_SAYFY_ABURUB_6_3_4, & ARKODE_CASH_KARP_6_4_5, ARKODE_FEHLBERG_6_4_5, ARKODE_DORMAND_PRINCE_7_4_5, ARKODE_ARK548L2SA_ERK_8_4_5, & ARKODE_VERNER_8_5_6, ARKODE_FEHLBERG_13_7_8, ARKODE_KNOTH_WOLKE_3_3, ARKODE_ARK437L2SA_ERK_7_3_4, & - ARKODE_ARK548L2SAb_ERK_8_4_5, ARKODE_MAX_ERK_NUM + ARKODE_ARK548L2SAb_ERK_8_4_5, ARKODE_ARK2_ERK_3_1_2, ARKODE_MAX_ERK_NUM public :: FARKodeButcherTable_LoadERK public :: FARKodeButcherTable_LoadERKByName integer(C_INT), parameter, public :: ARKLS_SUCCESS = 0_C_INT diff --git a/src/arkode/fmod/farkode_mristep_mod.c b/src/arkode/fmod/farkode_mristep_mod.c index fcd9d093b5..8ba7ad5296 100644 --- a/src/arkode/fmod/farkode_mristep_mod.c +++ b/src/arkode/fmod/farkode_mristep_mod.c @@ -1061,6 +1061,20 @@ SWIGEXPORT int _wrap_FMRIStepSetStopTime(void *farg1, double const *farg2) { } +SWIGEXPORT int _wrap_FMRIStepSetInterpolateStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)MRIStepSetInterpolateStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FMRIStepClearStopTime(void *farg1) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/arkode/fmod/farkode_mristep_mod.f90 b/src/arkode/fmod/farkode_mristep_mod.f90 index 2c0b32a367..b80ca727f8 100644 --- a/src/arkode/fmod/farkode_mristep_mod.f90 +++ b/src/arkode/fmod/farkode_mristep_mod.f90 @@ -160,6 +160,7 @@ module farkode_mristep_mod public :: FMRIStepSetNonlinConvCoef public :: FMRIStepSetMaxHnilWarns public :: FMRIStepSetStopTime + public :: FMRIStepSetInterpolateStopTime public :: FMRIStepClearStopTime public :: FMRIStepSetFixedStep public :: FMRIStepSetRootDirection @@ -734,6 +735,15 @@ function swigc_FMRIStepSetStopTime(farg1, farg2) & integer(C_INT) :: fresult end function +function swigc_FMRIStepSetInterpolateStopTime(farg1, farg2) & +bind(C, name="_wrap_FMRIStepSetInterpolateStopTime") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +integer(C_INT), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + function swigc_FMRIStepClearStopTime(farg1) & bind(C, name="_wrap_FMRIStepClearStopTime") & result(fresult) @@ -2313,6 +2323,22 @@ function FMRIStepSetStopTime(arkode_mem, tstop) & swig_result = fresult end function +function FMRIStepSetInterpolateStopTime(arkode_mem, interp) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_INT), intent(in) :: interp +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = arkode_mem +farg2 = interp +fresult = swigc_FMRIStepSetInterpolateStopTime(farg1, farg2) +swig_result = fresult +end function + function FMRIStepClearStopTime(arkode_mem) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/src/cvode/fmod/fcvode_mod.c b/src/cvode/fmod/fcvode_mod.c index db7e75ce06..5b63d2b605 100644 --- a/src/cvode/fmod/fcvode_mod.c +++ b/src/cvode/fmod/fcvode_mod.c @@ -628,6 +628,20 @@ SWIGEXPORT int _wrap_FCVodeSetStopTime(void *farg1, double const *farg2) { } +SWIGEXPORT int _wrap_FCVodeSetInterpolateStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)CVodeSetInterpolateStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FCVodeClearStopTime(void *farg1) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/cvode/fmod/fcvode_mod.f90 b/src/cvode/fmod/fcvode_mod.f90 index d107e94c9a..738f334a3f 100644 --- a/src/cvode/fmod/fcvode_mod.f90 +++ b/src/cvode/fmod/fcvode_mod.f90 @@ -103,6 +103,7 @@ module fcvode_mod public :: FCVodeSetNonlinearSolver public :: FCVodeSetStabLimDet public :: FCVodeSetStopTime + public :: FCVodeSetInterpolateStopTime public :: FCVodeClearStopTime public :: FCVodeSetUseIntegratorFusedKernels public :: FCVodeSetUserData @@ -471,6 +472,15 @@ function swigc_FCVodeSetStopTime(farg1, farg2) & integer(C_INT) :: fresult end function +function swigc_FCVodeSetInterpolateStopTime(farg1, farg2) & +bind(C, name="_wrap_FCVodeSetInterpolateStopTime") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +integer(C_INT), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + function swigc_FCVodeClearStopTime(farg1) & bind(C, name="_wrap_FCVodeClearStopTime") & result(fresult) @@ -1823,6 +1833,22 @@ function FCVodeSetStopTime(cvode_mem, tstop) & swig_result = fresult end function +function FCVodeSetInterpolateStopTime(cvode_mem, interp) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: cvode_mem +integer(C_INT), intent(in) :: interp +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = cvode_mem +farg2 = interp +fresult = swigc_FCVodeSetInterpolateStopTime(farg1, farg2) +swig_result = fresult +end function + function FCVodeClearStopTime(cvode_mem) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/src/cvodes/fmod/fcvodes_mod.c b/src/cvodes/fmod/fcvodes_mod.c index 1fa4e6bd9f..7c5b85e0af 100644 --- a/src/cvodes/fmod/fcvodes_mod.c +++ b/src/cvodes/fmod/fcvodes_mod.c @@ -702,6 +702,20 @@ SWIGEXPORT int _wrap_FCVodeSetStopTime(void *farg1, double const *farg2) { } +SWIGEXPORT int _wrap_FCVodeSetInterpolateStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)CVodeSetInterpolateStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FCVodeClearStopTime(void *farg1) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/cvodes/fmod/fcvodes_mod.f90 b/src/cvodes/fmod/fcvodes_mod.f90 index 39d8868e2f..505ed54f8a 100644 --- a/src/cvodes/fmod/fcvodes_mod.f90 +++ b/src/cvodes/fmod/fcvodes_mod.f90 @@ -133,6 +133,7 @@ module fcvodes_mod public :: FCVodeSetNonlinearSolver public :: FCVodeSetStabLimDet public :: FCVodeSetStopTime + public :: FCVodeSetInterpolateStopTime public :: FCVodeClearStopTime public :: FCVodeSetUserData public :: FCVodeSetEtaFixedStepBounds @@ -645,6 +646,15 @@ function swigc_FCVodeSetStopTime(farg1, farg2) & integer(C_INT) :: fresult end function +function swigc_FCVodeSetInterpolateStopTime(farg1, farg2) & +bind(C, name="_wrap_FCVodeSetInterpolateStopTime") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +integer(C_INT), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + function swigc_FCVodeClearStopTime(farg1) & bind(C, name="_wrap_FCVodeClearStopTime") & result(fresult) @@ -3229,6 +3239,22 @@ function FCVodeSetStopTime(cvode_mem, tstop) & swig_result = fresult end function +function FCVodeSetInterpolateStopTime(cvode_mem, interp) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: cvode_mem +integer(C_INT), intent(in) :: interp +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = cvode_mem +farg2 = interp +fresult = swigc_FCVodeSetInterpolateStopTime(farg1, farg2) +swig_result = fresult +end function + function FCVodeClearStopTime(cvode_mem) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/test/answers b/test/answers index 96d6e170c1..97011ebcf0 160000 --- a/test/answers +++ b/test/answers @@ -1 +1 @@ -Subproject commit 96d6e170c15f997d1e9062d4e6478e618d3f30ca +Subproject commit 97011ebcf0322a6372bbff86acf7a95c5546169c From 3247363e60c7cd8de5a6db1c0482328bcf2d1821 Mon Sep 17 00:00:00 2001 From: Shahbaj Sohal <89656941+shahbajsohal@users.noreply.github.com> Date: Mon, 17 Jul 2023 16:28:12 -0700 Subject: [PATCH 2/2] Bugfix - Creating pcg parallel directory (#297) There was an error where even when MPI was enabled the `../examples/sunlinsol/pcg/parallel` directory was not being built because it was not being added as a subdirectory in `../examples/sunlinsol/CMakeLists.txt`. --------- Co-authored-by: shahbajsohal Co-authored-by: Balos, Cody, J Co-authored-by: David J. Gardner --- examples/sunlinsol/CMakeLists.txt | 1 + examples/sunlinsol/pcg/parallel/CMakeLists.txt | 11 ++++++++++- .../pcg/parallel/test_sunlinsol_pcg_parallel.c | 11 +++++------ scripts/shared | 4 ++-- 4 files changed, 18 insertions(+), 9 deletions(-) diff --git a/examples/sunlinsol/CMakeLists.txt b/examples/sunlinsol/CMakeLists.txt index 3628997f36..e36b239234 100644 --- a/examples/sunlinsol/CMakeLists.txt +++ b/examples/sunlinsol/CMakeLists.txt @@ -39,6 +39,7 @@ endif() target_link_libraries(test_sunlinsol_obj PRIVATE sundials_sunlinsoldense) if(ENABLE_MPI AND MPI_C_FOUND) + add_subdirectory(pcg/parallel) add_subdirectory(spgmr/parallel) add_subdirectory(spfgmr/parallel) add_subdirectory(spbcgs/parallel) diff --git a/examples/sunlinsol/pcg/parallel/CMakeLists.txt b/examples/sunlinsol/pcg/parallel/CMakeLists.txt index e6a130c84d..8b298db30e 100644 --- a/examples/sunlinsol/pcg/parallel/CMakeLists.txt +++ b/examples/sunlinsol/pcg/parallel/CMakeLists.txt @@ -14,12 +14,21 @@ # CMakeLists.txt file for sunlinsol PCG examples # --------------------------------------------------------------- +# Set tolerance for linear solver test based on Sundials precision +if(SUNDIALS_PRECISION MATCHES "SINGLE") + set(TOL "1e-2") +elseif(SUNDIALS_PRECISION MATCHES "DOUBLE") + set(TOL "1e-10") +else() + set(TOL "1e-13") +endif() + # Example lists are tuples "name\;args\;nodes\;tasks\;type" where the # type is develop for examples excluded from 'make test' in releases # Examples using the SUNDIALS PCG linear solver set(sunlinsol_pcg_examples - "test_sunlinsol_pcg_parallel\;100 100 1e-3 0\;1\;4\;" + "test_sunlinsol_pcg_parallel\;100 200 ${TOL} 0\;1\;4\;" ) # Dependencies for nvector examples diff --git a/examples/sunlinsol/pcg/parallel/test_sunlinsol_pcg_parallel.c b/examples/sunlinsol/pcg/parallel/test_sunlinsol_pcg_parallel.c index 0963050fb7..644079a94d 100644 --- a/examples/sunlinsol/pcg/parallel/test_sunlinsol_pcg_parallel.c +++ b/examples/sunlinsol/pcg/parallel/test_sunlinsol_pcg_parallel.c @@ -103,7 +103,6 @@ int main(int argc, char *argv[]) { int fails=0; /* counter for test failures */ int passfail=0; /* overall passfail flag */ - int mpierr; /* mpi error flag */ SUNLinearSolver LS; /* linear solver object */ N_Vector xhat, x, b; /* test vectors */ UserData ProbData; /* problem data structure */ @@ -158,8 +157,8 @@ int main(int argc, char *argv[]) if (ProbData.myid == 0) { printf("\nPCG linear solver test:\n"); printf(" nprocs = %i\n", ProbData.nprocs); - printf(" local/global problem sizes = %ld/%ld\n", ProbData.Nloc, - ProbData.nprocs * ProbData.Nloc, sunctx); + printf(" local/global problem sizes = %ld/%ld\n", (long int) ProbData.Nloc, + (long int) ProbData.nprocs * ProbData.Nloc); printf(" Maximum Krylov subspace dimension = %i\n", maxl); printf(" Solver Tolerance = %g\n", tol); printf(" timing output flag = %i\n\n", print_timing); @@ -338,7 +337,7 @@ int main(int argc, char *argv[]) } /* check if any other process failed */ - mpierr = MPI_Allreduce(&passfail, &fails, 1, MPI_INT, MPI_MAX, ProbData.comm); + (void) MPI_Allreduce(&passfail, &fails, 1, MPI_INT, MPI_MAX, ProbData.comm); /* Free solver and vectors */ SUNLinSolFree(LS); @@ -489,8 +488,8 @@ static int check_flag(void *flagvalue, const char *funcname, int opt) * --------------------------------------------------------------------*/ int check_vector(N_Vector X, N_Vector Y, realtype tol) { - int failure = 0; - long int i; + int failure = 0; + sunindextype i; realtype *Xdata, *Ydata, maxerr; Xdata = N_VGetArrayPointer(X); diff --git a/scripts/shared b/scripts/shared index b9a7da6458..28ecda2d92 100755 --- a/scripts/shared +++ b/scripts/shared @@ -716,8 +716,8 @@ $tar $tarfile $distrobase/examples/sunlinsol/magmadense/CMakeLists.txt $tar $tarfile $distrobase/examples/sunlinsol/onemkldense/CMakeLists.txt $tar $tarfile $distrobase/examples/sunlinsol/onemkldense/test_sunlinsol_onemkldense.cpp -#$tar $tarfile $distrobase/examples/sunlinsol/pcg/parallel/test_sunlinsol_pcg_parallel.c -#$tar $tarfile $distrobase/examples/sunlinsol/pcg/parallel/CMakeLists.txt +$tar $tarfile $distrobase/examples/sunlinsol/pcg/parallel/test_sunlinsol_pcg_parallel.c +$tar $tarfile $distrobase/examples/sunlinsol/pcg/parallel/CMakeLists.txt $tar $tarfile $distrobase/examples/sunlinsol/pcg/serial/test_sunlinsol_pcg_serial.c $tar $tarfile $distrobase/examples/sunlinsol/pcg/serial/test_fsunlinsol_pcg_mod_serial.f90