diff --git a/CHANGELOG.md b/CHANGELOG.md index 5f40a1099d..43a8409da7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,13 @@ functions to queue and flush log messages. Updated `examples/cvode/petsc/cv_petsc_ex7.c` to support PETSc 3.25.0. +Added the function `ARKodeSkipAdaptStopTime` to specify that +stop-time-limited steps should be disregarded when selecting step sizes +for time step adaptivity. Added the functions `CVodeSkipAdaptStopTime` +and `IDASkipAdaptStopTime` to specify that stop-time-limited steps should +be disregarded when adapting the step size and method order for CVODE(S) +and IDA(S), respectively. + ### Bug Fixes Fixed memory leaks in CVODES, IDAS, and KINSOL in the unlikely event of a failed @@ -111,6 +118,37 @@ Fixed a bug where passing an empty string to corresponding logging stream ([Issue #844](https://github.com/llnl/sundials/issues/844)). +Fixed a bug in logging output from ARKODE, where for some time stepping modules, +the current "time" output in the logger was incorrect. + +Fixed a bug in the ARKODE discrete adjoint checkpointing where an incorrect +state would be stored on the first step if the output vector passed to +`ARKodeEvolve` did not contain the initial condition on the first call. + +Fixed a bug in MRIStep when using a custom inner integrator that relies on the +input state being the initial condition for the fast integration rather than +retaining the result from the last inner integration or most recent reset call +and the output vector passed to `ARKodeEvolve` does not contain the initial +condition on the first call or the last returned solution on subsequent calls. + +Removed an extraneous copy of the output vector in each step with SplittingStep. + +Added a missing call to `SUNNonlinSolSetup` in MRIStep when using an +IMEX-MRI-SR method. + +Fixed a potential bug in LSRKStep's `ARKODE_LSRK_SSP_S_3` method, where a real +number was used instead of an integer, potentially resulting in a rounding error. + +Fixed a bug in LSRKStep where an incorrect state vector could be passed to a +user-supplied dominant eigenvalue function on the first step unless the output +vector passed to `ARKodeEvolve` contained the initial condition and when an +eigenvalue estimate is requested on the first step in a subsequent call to +`ARKodeEvolve` unless the output vector passed contained the most recently returned +solution. + +Fixed a bug in MRIStep for estimating the first "slow" time step in an adaptive +multirate calculation. + ### Deprecation Notices The `CVodeSetMonitorFn` and `CVodeSetMonitorFrequency` functions have been @@ -120,7 +158,7 @@ Several CMake options have been deprecated in favor of namespaced versions prefixed with `SUNDIALS_` to avoid naming collisions in applications that include SUNDIALS directly within their CMake builds. Additionally, a consistent naming convention (`SUNDIALS_ENABLE`) is now used for all boolean options. The -table below lists the old CMake option names and the new replacements. +Removed an extraneous copy of the output vector in each step with SplittingStep. | Old Option | New Option | |-----------------------------------------|------------------------------------------------| diff --git a/bindings/sundials4py/arkode/arkode_generated.hpp b/bindings/sundials4py/arkode/arkode_generated.hpp index 26cb224ad6..4c457cf204 100644 --- a/bindings/sundials4py/arkode/arkode_generated.hpp +++ b/bindings/sundials4py/arkode/arkode_generated.hpp @@ -203,6 +203,9 @@ m.def("ARKodeSetStopTime", ARKodeSetStopTime, nb::arg("arkode_mem"), m.def("ARKodeClearStopTime", ARKodeClearStopTime, nb::arg("arkode_mem")); +m.def("ARKodeSkipAdaptStopTime", ARKodeSkipAdaptStopTime, nb::arg("arkode_mem"), + nb::arg("skipadapttstop")); + m.def("ARKodeSetFixedStep", ARKodeSetFixedStep, nb::arg("arkode_mem"), nb::arg("hfixed")); diff --git a/bindings/sundials4py/cvodes/cvodes_generated.hpp b/bindings/sundials4py/cvodes/cvodes_generated.hpp index dbbde6f0ad..540de53e48 100644 --- a/bindings/sundials4py/cvodes/cvodes_generated.hpp +++ b/bindings/sundials4py/cvodes/cvodes_generated.hpp @@ -166,6 +166,9 @@ m.def("CVodeSetInterpolateStopTime", CVodeSetInterpolateStopTime, m.def("CVodeClearStopTime", CVodeClearStopTime, nb::arg("cvode_mem")); +m.def("CVodeSkipAdaptStopTime", CVodeSkipAdaptStopTime, nb::arg("cvode_mem"), + nb::arg("skipadapttstop")); + m.def("CVodeSetEtaFixedStepBounds", CVodeSetEtaFixedStepBounds, nb::arg("cvode_mem"), nb::arg("eta_min_fx"), nb::arg("eta_max_fx")); diff --git a/bindings/sundials4py/idas/idas_generated.hpp b/bindings/sundials4py/idas/idas_generated.hpp index 02902af5a2..317e2f559e 100644 --- a/bindings/sundials4py/idas/idas_generated.hpp +++ b/bindings/sundials4py/idas/idas_generated.hpp @@ -117,6 +117,9 @@ m.def("IDASetStopTime", IDASetStopTime, nb::arg("ida_mem"), nb::arg("tstop")); m.def("IDAClearStopTime", IDAClearStopTime, nb::arg("ida_mem")); +m.def("IDASkipAdaptStopTime", IDASkipAdaptStopTime, nb::arg("ida_mem"), + nb::arg("skipadapttstop")); + m.def("IDASetMaxErrTestFails", IDASetMaxErrTestFails, nb::arg("ida_mem"), nb::arg("maxnef")); diff --git a/doc/arkode/guide/source/Usage/User_callable.rst b/doc/arkode/guide/source/Usage/User_callable.rst index ea67827470..090e5ab28f 100644 --- a/doc/arkode/guide/source/Usage/User_callable.rst +++ b/doc/arkode/guide/source/Usage/User_callable.rst @@ -905,6 +905,7 @@ Minimum absolute step size :c:func:`ARKodeSetMinStep` Set a value for :math:`t_{stop}` :c:func:`ARKodeSetStopTime` undefined Interpolate at :math:`t_{stop}` :c:func:`ARKodeSetInterpolateStopTime` ``SUNFALSE`` Disable the stop time :c:func:`ARKodeClearStopTime` N/A +Disregard stop time limited steps in adaptivity :c:func:`ARKodeSkipAdaptStopTime` ``SUNFALSE`` Supply a pointer for user data :c:func:`ARKodeSetUserData` ``NULL`` Maximum no. of ARKODE error test failures :c:func:`ARKodeSetMaxErrTestFails` 7 Set inequality constraints on solution :c:func:`ARKodeSetConstraints` ``NULL`` @@ -1466,6 +1467,30 @@ Use compensated summation for accumulating time :c:func:`ARKodeSetUseCompensa .. versionadded:: 6.1.0 +.. c:function:: int ARKodeSkipAdaptStopTime(void* arkode_mem, sunbooleantype skip) + + Specifies whether stop-time-limited steps should be disregarded + when selecting step sizes for time step adaptivity. + + :param arkode_mem: pointer to the ARKODE memory block. + :param skip: flag indicating to disregard (1) or retain (0) stop time + limited steps from the temporal adaptivity algorithm. + + :retval ARK_SUCCESS: the function exited successfully. + :retval ARK_MEM_NULL: ``arkode_mem`` was ``NULL``. + + .. note:: + + The default behavior is to use all successful time steps + (including stop-time-limited steps) when determining an + adapted time step size. + + This routine will be called by :c:func:`ARKodeSetOptions` + when using the key "arkid.skip_adapt_stop_time". + + .. versionadded:: x.y.z + + .. c:function:: int ARKodeSetUserData(void* arkode_mem, void* user_data) Specifies the user data block *user_data* and diff --git a/doc/cvode/guide/source/Usage/index.rst b/doc/cvode/guide/source/Usage/index.rst index 2929ab59c4..95bf791f8b 100644 --- a/doc/cvode/guide/source/Usage/index.rst +++ b/doc/cvode/guide/source/Usage/index.rst @@ -864,6 +864,9 @@ Main solver optional input functions +-------------------------------+---------------------------------------------+----------------+ | Disable the stop time | :c:func:`CVodeClearStopTime` | N/A | +-------------------------------+---------------------------------------------+----------------+ + | Disregard stop time limited | :c:func:`CVodeSkipAdaptStopTime` | ``SUNFALSE`` | + | steps in adaptivity | | | + +-------------------------------+---------------------------------------------+----------------+ | Maximum no. of error test | :c:func:`CVodeSetMaxErrTestFails` | 7 | | failures | | | +-------------------------------+---------------------------------------------+----------------+ @@ -1207,6 +1210,28 @@ Main solver optional input functions .. versionadded:: 6.5.1 +.. c:function:: int CVodeSkipAdaptStopTime(void* cvode_mem, sunbooleantype skip) + + Specifies whether stop-time-limited steps should be disregarded + when adapting step sizes and method order. + + **Arguments:** + * ``cvode_mem`` -- pointer to the CVODE memory block. + * ``skip`` -- flag indicating to disregard (1) or retain (0) stop time limited steps from the temporal adaptivity algorithm. + + **Return value:** + * ``CV_SUCCESS`` if successful + * ``CV_MEM_NULL`` if the CVODE memory is ``NULL`` + + **Notes:** + The default behavior is to use all successful time steps + (including stop-time-limited steps) when performing adaptivity. + + This routine will be called by :c:func:`CVodeSetOptions` + when using the key "cvid.skip_adapt_stop_time". + + .. versionadded:: x.y.z + .. c:function:: int CVodeSetMaxErrTestFails(void* cvode_mem, int maxnef) The function ``CVodeSetMaxErrTestFails`` specifies the maximum number of error test failures permitted in attempting one step. diff --git a/doc/cvodes/guide/source/Usage/SIM.rst b/doc/cvodes/guide/source/Usage/SIM.rst index 5ee88ad335..d71366789a 100644 --- a/doc/cvodes/guide/source/Usage/SIM.rst +++ b/doc/cvodes/guide/source/Usage/SIM.rst @@ -873,6 +873,9 @@ Main solver optional input functions +-------------------------------+---------------------------------------------+----------------+ | Disable the stop time | :c:func:`CVodeClearStopTime` | N/A | +-------------------------------+---------------------------------------------+----------------+ + | Disregard stop time limited | :c:func:`CVodeSkipAdaptStopTime` | ``SUNFALSE`` | + | steps in adaptivity | | | + +-------------------------------+---------------------------------------------+----------------+ | Maximum no. of error test | :c:func:`CVodeSetMaxErrTestFails` | 7 | | failures | | | +-------------------------------+---------------------------------------------+----------------+ @@ -1213,6 +1216,28 @@ Main solver optional input functions .. versionadded:: 6.5.1 +.. c:function:: int CVodeSkipAdaptStopTime(void* cvode_mem, sunbooleantype skip) + + Specifies whether stop-time-limited steps should be disregarded + when adapting step sizes and method order. + + **Arguments:** + * ``cvode_mem`` -- pointer to the CVODE memory block. + * ``skip`` -- flag indicating to disregard (1) or retain (0) stop time limited steps from the temporal adaptivity algorithm. + + **Return value:** + * ``CV_SUCCESS`` if successful + * ``CV_MEM_NULL`` if the CVODE memory is ``NULL`` + + **Notes:** + The default behavior is to use all successful time steps + (including stop-time-limited steps) when performing adaptivity. + + This routine will be called by :c:func:`CVodeSetOptions` + when using the key "cvid.skip_adapt_stop_time". + + .. versionadded:: x.y.z + .. c:function:: int CVodeSetMaxErrTestFails(void* cvode_mem, int maxnef) The function ``CVodeSetMaxErrTestFails`` specifies the maximum number of error test failures permitted in attempting one step. diff --git a/doc/ida/guide/source/Usage/index.rst b/doc/ida/guide/source/Usage/index.rst index 2f3b992c04..002b90be9e 100644 --- a/doc/ida/guide/source/Usage/index.rst +++ b/doc/ida/guide/source/Usage/index.rst @@ -931,6 +931,8 @@ Main solver optional input functions +--------------------------------------------------------------------+---------------------------------------+----------------+ | Disable the stop time | :c:func:`IDAClearStopTime` | N/A | +--------------------------------------------------------------------+---------------------------------------+----------------+ + | Disregard stop time limited steps in adaptivity | :c:func:`IDASkipAdaptStopTime` | ``SUNFALSE`` | + +--------------------------------------------------------------------+---------------------------------------+----------------+ | Maximum no. of error test failures | :c:func:`IDASetMaxErrTestFails` | 10 | +--------------------------------------------------------------------+---------------------------------------+----------------+ | Suppress alg. vars. from error test | :c:func:`IDASetSuppressAlg` | ``SUNFALSE`` | @@ -1189,6 +1191,28 @@ Main solver optional input functions .. versionadded:: 6.5.1 +.. c:function:: int IDASkipAdaptStopTime(void* ida_mem, sunbooleantype skip) + + Specifies whether stop-time-limited steps should be disregarded + when adapting step sizes and method order. + + **Arguments:** + * ``ida_mem`` -- pointer to the IDA memory block. + * ``skip`` -- flag indicating to disregard (1) or retain (0) stop time limited steps from the temporal adaptivity algorithm. + + **Return value:** + * ``IDA_SUCCESS`` if successful + * ``IDA_MEM_NULL`` if the IDA memory is ``NULL`` + + **Notes:** + The default behavior is to use all successful time steps + (including stop-time-limited steps) when performing adaptivity. + + This routine will be called by :c:func:`IDASetOptions` + when using the key "idaid.skip_adapt_stop_time". + + .. versionadded:: x.y.z + .. c:function:: int IDASetMaxErrTestFails(void * ida_mem, int maxnef) The function ``IDASetMaxErrTestFails`` specifies the maximum number of error diff --git a/doc/idas/guide/source/Usage/SIM.rst b/doc/idas/guide/source/Usage/SIM.rst index 1630210b12..1fbb8f898c 100644 --- a/doc/idas/guide/source/Usage/SIM.rst +++ b/doc/idas/guide/source/Usage/SIM.rst @@ -936,6 +936,8 @@ Main solver optional input functions +--------------------------------------------------------------------+---------------------------------------+----------------+ | Disable the stop time | :c:func:`IDAClearStopTime` | N/A | +--------------------------------------------------------------------+---------------------------------------+----------------+ + | Disregard stop time limited steps in adaptivity | :c:func:`IDASkipAdaptStopTime` | ``SUNFALSE`` | + +--------------------------------------------------------------------+---------------------------------------+----------------+ | Maximum no. of error test failures | :c:func:`IDASetMaxErrTestFails` | 10 | +--------------------------------------------------------------------+---------------------------------------+----------------+ | Suppress alg. vars. from error test | :c:func:`IDASetSuppressAlg` | ``SUNFALSE`` | @@ -1194,6 +1196,28 @@ Main solver optional input functions .. versionadded:: 6.5.1 +.. c:function:: int IDASkipAdaptStopTime(void* ida_mem, sunbooleantype skip) + + Specifies whether stop-time-limited steps should be disregarded + when adapting step sizes and method order. + + **Arguments:** + * ``ida_mem`` -- pointer to the IDA memory block. + * ``skip`` -- flag indicating to disregard (1) or retain (0) stop time limited steps from the temporal adaptivity algorithm. + + **Return value:** + * ``IDA_SUCCESS`` if successful + * ``IDA_MEM_NULL`` if the IDA memory is ``NULL`` + + **Notes:** + The default behavior is to use all successful time steps + (including stop-time-limited steps) when performing adaptivity. + + This routine will be called by :c:func:`IDASetOptions` + when using the key "idaid.skip_adapt_stop_time". + + .. versionadded:: x.y.z + .. c:function:: int IDASetMaxErrTestFails(void * ida_mem, int maxnef) The function :c:func:`IDASetMaxErrTestFails` specifies the maximum number of error diff --git a/doc/shared/RecentChanges.rst b/doc/shared/RecentChanges.rst index d1855658a3..8340be1afb 100644 --- a/doc/shared/RecentChanges.rst +++ b/doc/shared/RecentChanges.rst @@ -10,6 +10,13 @@ user-defined functions to queue and flush log messages. Updated ``examples/cvode/petsc/cv_petsc_ex7.c`` to support PETSc 3.25.0. +Added the function :c:func:`ARKodeSkipAdaptStopTime` to specify that +stop-time-limited steps should be disregarded when selecting step sizes +for time step adaptivity. Added the functions +:c:func:`CVodeSkipAdaptStopTime` and :c:func:`IDASkipAdaptStopTime` +to specify that stop-time-limited steps should be disregarded when +adapting the step size and method order for CVODE(S) and IDA(S), respectively. + **Bug Fixes** Fixed memory leaks in CVODES, IDAS, and KINSOL in the unlikely event of a failed diff --git a/examples/arkode/CXX_serial/CMakeLists.txt b/examples/arkode/CXX_serial/CMakeLists.txt index 0c6a0b6279..a6a4b4ef3a 100644 --- a/examples/arkode/CXX_serial/CMakeLists.txt +++ b/examples/arkode/CXX_serial/CMakeLists.txt @@ -39,6 +39,7 @@ set(ARKODE_examples "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_kpr_nestedmri.cpp\;\;exclude-single" + "ark_kpr_nestedmri.cpp\;arkfast.skip_adapt_stop_time 1 arkmod.skip_adapt_stop_time 1\;exclude-single" "ark_pendulum.cpp\;\;develop") # Header files to install diff --git a/examples/arkode/CXX_serial/ark_kpr_nestedmri.cpp b/examples/arkode/CXX_serial/ark_kpr_nestedmri.cpp index c574f5814f..27cc64edba 100644 --- a/examples/arkode/CXX_serial/ark_kpr_nestedmri.cpp +++ b/examples/arkode/CXX_serial/ark_kpr_nestedmri.cpp @@ -451,6 +451,8 @@ int main(int argc, char* argv[]) if (check_flag(retval, "ARKodeSetMaxNumSteps")) return 1; retval = ARKodeSetUserData(inner_arkode_mem, (void*)&opts); if (check_flag(retval, "ARKodeSetUserData")) return 1; + retval = ARKodeSetOptions(inner_arkode_mem, "arkfast", nullptr, argc, argv); + if (check_flag(retval, "ARKodeSetOptions")) return 1; // Create inner stepper MRIStepInnerStepper inner_stepper = nullptr; @@ -873,6 +875,8 @@ int main(int argc, char* argv[]) retval = ARKodeSetFixedStep(mid_arkode_mem, opts.hm); if (check_flag(retval, "ARKodeSetFixedStep")) return 1; } + retval = ARKodeSetOptions(mid_arkode_mem, "arkmid", nullptr, argc, argv); + if (check_flag(retval, "ARKodeSetOptions")) return 1; // Create intermediate stepper MRIStepInnerStepper intermediate_stepper = nullptr; @@ -931,6 +935,8 @@ int main(int argc, char* argv[]) retval = ARKodeSetFixedStep(arkode_mem, opts.hs); if (check_flag(retval, "ARKodeSetFixedStep")) return 1; } + retval = ARKodeSetOptions(arkode_mem, "arkslow", nullptr, argc, argv); + if (check_flag(retval, "ARKodeSetOptions")) return 1; // // Integrate ODE diff --git a/examples/arkode/CXX_serial/ark_kpr_nestedmri_arkfast.skip_adapt_stop_time_1_arkmod.skip_adapt_stop_time_1.out b/examples/arkode/CXX_serial/ark_kpr_nestedmri_arkfast.skip_adapt_stop_time_1_arkmod.skip_adapt_stop_time_1.out new file mode 100644 index 0000000000..1d1d669c09 --- /dev/null +++ b/examples/arkode/CXX_serial/ark_kpr_nestedmri_arkfast.skip_adapt_stop_time_1_arkmod.skip_adapt_stop_time_1.out @@ -0,0 +1,51 @@ + +Adaptive nested multirate nonlinear Kvaerno-Prothero-Robinson test problem: + time domain: (0,5] + G = -10 + e = 0.5 + al = -1 + be = 1 + om = 50 + + Slow integrator: ARKODE_MRI_GARK_ERK45a (explicit) + + Intermediate integrator: ARKODE_MRI_GARK_ERK45a (explicit) + MRI-HTOL controller (using I for H) based on order of MRI embedding + rtol = 0.0001, atol = 1e-11 + fast error accumulation strategy = 0 + + Fast order 4 + I controller for fast time scale, based on order of RK embedding + fast_rtol = 0.0001, atol = 1e-11 + t u v w uerr verr werr + ---------------------------------------------------------------------------- + 0.000000 1.581139 1.732051 1.732051 0.00e+00 0.00e+00 0.00e+00 + 0.250000 1.575568 1.691653 0.982638 8.18e-06 4.44e-07 9.29e-05 + 0.500000 1.560231 1.085823 1.209956 1.49e-05 8.07e-06 4.30e-04 + 0.750000 1.536102 1.469088 1.256444 3.00e-05 3.10e-05 2.84e-03 + 1.000000 1.504205 1.644139 1.534738 5.98e-05 8.65e-06 4.35e-04 + 1.250000 1.465597 1.089258 1.075900 2.87e-05 3.16e-06 1.15e-04 + 1.500000 1.422645 1.424070 1.332193 2.34e-06 7.84e-09 3.97e-08 + 1.750000 1.377092 1.701792 1.324408 4.49e-06 2.63e-05 3.93e-04 + 2.000000 1.332370 1.534071 0.886999 1.10e-10 3.34e-13 6.01e-12 + 2.250000 1.290314 1.296504 1.093453 5.11e-06 1.44e-05 3.29e-03 + 2.500000 1.254200 1.020992 0.935999 1.43e-08 5.85e-10 1.37e-07 + 2.750000 1.229247 1.082903 1.598229 8.41e-06 2.44e-06 5.59e-04 + 3.000000 1.216144 1.091610 1.130624 2.39e-05 4.74e-06 1.40e-03 + 3.250000 1.216297 1.240783 1.639081 3.44e-06 8.55e-07 6.50e-06 + 3.500000 1.230701 1.413878 1.540425 6.51e-05 1.73e-04 2.06e-03 + 3.750000 1.254951 1.365115 1.401758 2.29e-05 7.27e-06 2.13e-04 + 4.000000 1.287747 0.974938 1.690643 3.14e-06 1.77e-05 2.69e-03 + 4.250000 1.328130 1.675308 0.966312 2.42e-06 7.03e-09 1.21e-05 + 4.500000 1.371964 1.605167 1.479730 5.60e-07 4.61e-08 2.75e-08 + 4.750000 1.417588 1.514725 1.160359 8.37e-06 3.40e-07 4.86e-06 + 5.000000 1.461167 1.473354 1.644126 6.42e-07 1.35e-06 9.85e-05 + ---------------------------------------------------------------------------- + +Final Solver Statistics: + Slow steps = 239 (attempts = 309, fails = 70, innerfails = 0) + Intermediate steps = 2820 (attempts = 2955, fails = 135, innerfails = 0) + Fast steps = 87774 (attempts = 122683, fails = 34909) + u error = 9.45885e-05, v error = 0.000106298, total error = 0.00175566 + Relative accuracy = 94.0089 + Total RHS evals: Fse = 1477, Fsi = 0, Fme = 14883, Fmi = 0, Ff = 496888 diff --git a/examples/arkode/C_serial/ark_analytic.out b/examples/arkode/C_serial/ark_analytic.out index f822f3862f..7e9c0055ff 100644 --- a/examples/arkode/C_serial/ark_analytic.out +++ b/examples/arkode/C_serial/ark_analytic.out @@ -29,6 +29,7 @@ Soderlind SUNAdaptController module: previous-previous step = 1 firststeps = 0 historysize = 0 + Skip tstop-limited steps from affecting temporal adaptivity = 0 Maximum number of error test failures = 7 Maximum number of convergence test failures = 10 ARKStep time step module parameters: diff --git a/examples/arkode/C_serial/ark_analytic_arkode.scalar_tolerances_1e-6_1e-8_arkode.table_names_ARKODE_ESDIRK547L2SA_7_4_5_ARKODE_ERK_NONE.out b/examples/arkode/C_serial/ark_analytic_arkode.scalar_tolerances_1e-6_1e-8_arkode.table_names_ARKODE_ESDIRK547L2SA_7_4_5_ARKODE_ERK_NONE.out index 78af88676a..136fe90a32 100644 --- a/examples/arkode/C_serial/ark_analytic_arkode.scalar_tolerances_1e-6_1e-8_arkode.table_names_ARKODE_ESDIRK547L2SA_7_4_5_ARKODE_ERK_NONE.out +++ b/examples/arkode/C_serial/ark_analytic_arkode.scalar_tolerances_1e-6_1e-8_arkode.table_names_ARKODE_ESDIRK547L2SA_7_4_5_ARKODE_ERK_NONE.out @@ -29,6 +29,7 @@ Soderlind SUNAdaptController module: previous-previous step = 1 firststeps = 0 historysize = 0 + Skip tstop-limited steps from affecting temporal adaptivity = 0 Maximum number of error test failures = 7 Maximum number of convergence test failures = 10 ARKStep time step module parameters: diff --git a/examples/cvode/CXX_serial/cv_kpr.cpp b/examples/cvode/CXX_serial/cv_kpr.cpp index ff0aad9f79..b8ada93a92 100644 --- a/examples/cvode/CXX_serial/cv_kpr.cpp +++ b/examples/cvode/CXX_serial/cv_kpr.cpp @@ -140,6 +140,12 @@ int main(int argc, char* argv[]) // Advance in time for (int i = 0; i < opts.nout; i++) { + if (opts.use_tstop) + { + flag = CVodeSetStopTime(cvode_mem, tout); + if (check_flag(flag, "CVodeSetStopTime")) { return 1; } + } + flag = CVode(cvode_mem, tout, y, &tret, CV_NORMAL); if (check_flag(flag, "CVode")) { return 1; } diff --git a/examples/cvode/CXX_serial/cv_kpr.hpp b/examples/cvode/CXX_serial/cv_kpr.hpp index ef5dd77d84..272e70c11d 100644 --- a/examples/cvode/CXX_serial/cv_kpr.hpp +++ b/examples/cvode/CXX_serial/cv_kpr.hpp @@ -54,6 +54,9 @@ struct Options // Finite difference Jacobian bool fd_jac = false; + // Enforce stop time exactly (instead of interpolating) + bool use_tstop = false; + // Output options sunrealtype dtout = ONE; // output interval int nout = 10; // number of outputs @@ -99,6 +102,7 @@ static void InputHelp() std::cout << " --fdjac : finite-difference Jacobian\n"; std::cout << " --dtout : output interval\n"; std::cout << " --nout : number of outputs\n"; + std::cout << " --tstop : enforce stop time exactly\n"; } static int ReadInputs(std::vector& args, Options& opts, @@ -116,6 +120,7 @@ static int ReadInputs(std::vector& args, Options& opts, find_arg(args, "--fdjac", opts.fd_jac); find_arg(args, "--dtout", opts.dtout); find_arg(args, "--nout", opts.nout); + find_arg(args, "--tstop", opts.use_tstop); return 0; } diff --git a/include/arkode/arkode.h b/include/arkode/arkode.h index 6fae986edc..5adf808d3c 100644 --- a/include/arkode/arkode.h +++ b/include/arkode/arkode.h @@ -289,6 +289,8 @@ SUNDIALS_EXPORT int ARKodeSetInterpolateStopTime(void* arkode_mem, sunbooleantype interp); SUNDIALS_EXPORT int ARKodeSetStopTime(void* arkode_mem, sunrealtype tstop); SUNDIALS_EXPORT int ARKodeClearStopTime(void* arkode_mem); +SUNDIALS_EXPORT int ARKodeSkipAdaptStopTime(void* arkode_mem, + sunbooleantype skipadapttstop); SUNDIALS_EXPORT int ARKodeSetFixedStep(void* arkode_mem, sunrealtype hfixed); SUNDIALS_EXPORT int ARKodeSetStepDirection(void* arkode_mem, sunrealtype stepdir); SUNDIALS_EXPORT int ARKodeSetUserData(void* arkode_mem, void* user_data); diff --git a/include/cvode/cvode.h b/include/cvode/cvode.h index e08f462e0a..624bfcf7c5 100644 --- a/include/cvode/cvode.h +++ b/include/cvode/cvode.h @@ -155,6 +155,8 @@ SUNDIALS_EXPORT int CVodeSetStopTime(void* cvode_mem, sunrealtype tstop); SUNDIALS_EXPORT int CVodeSetInterpolateStopTime(void* cvode_mem, sunbooleantype interp); SUNDIALS_EXPORT int CVodeClearStopTime(void* cvode_mem); +SUNDIALS_EXPORT int CVodeSkipAdaptStopTime(void* cvode_mem, + sunbooleantype skipadapttstop); SUNDIALS_EXPORT int CVodeSetUseIntegratorFusedKernels(void* cvode_mem, sunbooleantype onoff); SUNDIALS_EXPORT int CVodeSetUserData(void* cvode_mem, void* user_data); diff --git a/include/cvodes/cvodes.h b/include/cvodes/cvodes.h index f5e7621a85..92d05ff6f5 100644 --- a/include/cvodes/cvodes.h +++ b/include/cvodes/cvodes.h @@ -224,6 +224,8 @@ SUNDIALS_EXPORT int CVodeSetStopTime(void* cvode_mem, sunrealtype tstop); SUNDIALS_EXPORT int CVodeSetInterpolateStopTime(void* cvode_mem, sunbooleantype interp); SUNDIALS_EXPORT int CVodeClearStopTime(void* cvode_mem); +SUNDIALS_EXPORT int CVodeSkipAdaptStopTime(void* cvode_mem, + sunbooleantype skipadapttstop); SUNDIALS_EXPORT int CVodeSetUserData(void* cvode_mem, void* user_data); /* Optional step adaptivity input functions */ diff --git a/include/ida/ida.h b/include/ida/ida.h index ab5a4bcff8..c217189a59 100644 --- a/include/ida/ida.h +++ b/include/ida/ida.h @@ -146,6 +146,8 @@ SUNDIALS_EXPORT int IDASetMaxStep(void* ida_mem, sunrealtype hmax); SUNDIALS_EXPORT int IDASetMinStep(void* ida_mem, sunrealtype hmin); SUNDIALS_EXPORT int IDASetStopTime(void* ida_mem, sunrealtype tstop); SUNDIALS_EXPORT int IDAClearStopTime(void* ida_mem); +SUNDIALS_EXPORT int IDASkipAdaptStopTime(void* ida_mem, + sunbooleantype skipadapttstop); SUNDIALS_EXPORT int IDASetMaxErrTestFails(void* ida_mem, int maxnef); SUNDIALS_EXPORT int IDASetSuppressAlg(void* ida_mem, sunbooleantype suppressalg); SUNDIALS_EXPORT int IDASetId(void* ida_mem, N_Vector id); diff --git a/include/idas/idas.h b/include/idas/idas.h index 310777d624..7768465550 100644 --- a/include/idas/idas.h +++ b/include/idas/idas.h @@ -210,6 +210,8 @@ SUNDIALS_EXPORT int IDASetMaxStep(void* ida_mem, sunrealtype hmax); SUNDIALS_EXPORT int IDASetMinStep(void* ida_mem, sunrealtype hmin); SUNDIALS_EXPORT int IDASetStopTime(void* ida_mem, sunrealtype tstop); SUNDIALS_EXPORT int IDAClearStopTime(void* ida_mem); +SUNDIALS_EXPORT int IDASkipAdaptStopTime(void* ida_mem, + sunbooleantype skipadapttstop); SUNDIALS_EXPORT int IDASetMaxErrTestFails(void* ida_mem, int maxnef); SUNDIALS_EXPORT int IDASetSuppressAlg(void* ida_mem, sunbooleantype suppressalg); SUNDIALS_EXPORT int IDASetId(void* ida_mem, N_Vector id); diff --git a/src/arkode/arkode.c b/src/arkode/arkode.c index 453814a25e..d56231257f 100644 --- a/src/arkode/arkode.c +++ b/src/arkode/arkode.c @@ -128,6 +128,7 @@ int ARKodeResize(void* arkode_mem, N_Vector y0, sunrealtype hscale, ark_mem->hprime *= hscale; /* If next step would overtake tstop, adjust stepsize */ + ark_mem->tstoplimited = SUNFALSE; if (ark_mem->tstopset) { if ((ark_mem->tcur + ark_mem->hprime - ark_mem->tstop) * ark_mem->hprime > @@ -136,6 +137,7 @@ int ARKodeResize(void* arkode_mem, N_Vector y0, sunrealtype hscale, ark_mem->hprime = (ark_mem->tstop - ark_mem->tcur) * (ONE - FOUR * ark_mem->uround); ark_mem->eta = ark_mem->hprime / ark_mem->h; + if (ark_mem->skipadapttstop) { ark_mem->tstoplimited = SUNTRUE; } } } } @@ -866,8 +868,8 @@ int ARKodeEvolve(void* arkode_mem, sunrealtype tout, N_Vector yout, /* Update parameter for upcoming step size */ if (ark_mem->hprime != ark_mem->h) { - ark_mem->h = ark_mem->h * ark_mem->eta; - ark_mem->next_h = ark_mem->h; + ark_mem->h = ark_mem->h * ark_mem->eta; + if (!ark_mem->skipadapttstop) ark_mem->next_h = ark_mem->h; } if (ark_mem->fixedstep) { @@ -1093,6 +1095,7 @@ int ARKodeEvolve(void* arkode_mem, sunrealtype tout, N_Vector yout, } /* Check if tn is at tstop or near tstop */ + ark_mem->tstoplimited = SUNFALSE; if (ark_mem->tstopset) { troundoff = FUZZ_FACTOR * ark_mem->uround * @@ -1129,6 +1132,7 @@ int ARKodeEvolve(void* arkode_mem, sunrealtype tout, N_Vector yout, ark_mem->hprime = (ark_mem->tstop - ark_mem->tcur) * (ONE - FOUR * ark_mem->uround); ark_mem->eta = ark_mem->hprime / ark_mem->h; + if (ark_mem->skipadapttstop) ark_mem->tstoplimited = SUNTRUE; } } @@ -1344,6 +1348,8 @@ void ARKodePrintMem(void* arkode_mem, FILE* outfile) fprintf(outfile, "user_efun = %i\n", ark_mem->user_efun); fprintf(outfile, "tstopset = %i\n", ark_mem->tstopset); fprintf(outfile, "tstopinterp = %i\n", ark_mem->tstopinterp); + fprintf(outfile, "tstoplimited = %i\n", ark_mem->tstoplimited); + fprintf(outfile, "skipadapttstop = %i\n", ark_mem->skipadapttstop); fprintf(outfile, "tstop = " SUN_FORMAT_G "\n", ark_mem->tstop); fprintf(outfile, "VabstolMallocDone = %i\n", ark_mem->VabstolMallocDone); fprintf(outfile, "MallocDone = %i\n", ark_mem->MallocDone); @@ -2306,6 +2312,7 @@ int arkInitialSetup(ARKodeMem ark_mem, sunrealtype tout) ark_mem->h0u = ark_mem->h; ark_mem->eta = ONE; ark_mem->hprime = ark_mem->h; + ark_mem->hold = ark_mem->h; } else { @@ -2812,8 +2819,9 @@ int arkCompleteStep(ARKodeMem ark_mem, sunrealtype dsm) N_VScale(ONE, ark_mem->ycur, ark_mem->yn); ark_mem->fn_is_current = SUNFALSE; - /* Notify time step controller object of successful step */ - if (ark_mem->hadapt_mem->hcontroller) + /* Notify time step controller object of successful step + (skip this if the previous step was stop-time-limited) */ + if (ark_mem->hadapt_mem->hcontroller && !ark_mem->tstoplimited) { retval = SUNAdaptController_UpdateH(ark_mem->hadapt_mem->hcontroller, ark_mem->h, dsm); @@ -3451,6 +3459,14 @@ int arkCheckTemporalError(ARKodeMem ark_mem, int* nflagPtr, int* nefPtr, } hadapt_mem = ark_mem->hadapt_mem; + /* If a stop-time-limited step will succeed, set eta so that the step size + reverts to the last "full size" successful step for the next attempt */ + if (ark_mem->tstoplimited && dsm <= ONE) + { + ark_mem->eta = ark_mem->hold / ark_mem->h; + return (ARK_SUCCESS); + } + /* consider change of step size for next step attempt (may be larger/smaller than current step, depending on dsm) */ ttmp = (dsm <= ONE) ? ark_mem->tn + ark_mem->h : ark_mem->tn; diff --git a/src/arkode/arkode_cli.c b/src/arkode/arkode_cli.c index 91aa2c78f4..ad49060e41 100644 --- a/src/arkode/arkode_cli.c +++ b/src/arkode/arkode_cli.c @@ -84,6 +84,7 @@ static int arkSetFromCommandLine(void* arkode_mem, const char* arkid, int argc, {"max_nonlin_iters", ARKodeSetMaxNonlinIters}, {"max_hnil_warns", ARKodeSetMaxHnilWarns}, {"interpolate_stop_time", ARKodeSetInterpolateStopTime}, + {"skip_adapt_stop_time", ARKodeSkipAdaptStopTime}, {"max_num_constr_fails", ARKodeSetMaxNumConstrFails}, {"adaptivity_adjustment", ARKodeSetAdaptivityAdjustment}, {"small_num_efails", ARKodeSetSmallNumEFails}, diff --git a/src/arkode/arkode_impl.h b/src/arkode/arkode_impl.h index 86cedbf3ae..e785d593af 100644 --- a/src/arkode/arkode_impl.h +++ b/src/arkode/arkode_impl.h @@ -502,6 +502,8 @@ struct ARKodeMemRec /* Tstop information */ sunbooleantype tstopset; sunbooleantype tstopinterp; + sunbooleantype tstoplimited; + sunbooleantype skipadapttstop; sunrealtype tstop; /* Time step data */ diff --git a/src/arkode/arkode_io.c b/src/arkode/arkode_io.c index 2a179b0b95..bf887ad897 100644 --- a/src/arkode/arkode_io.c +++ b/src/arkode/arkode_io.c @@ -80,12 +80,14 @@ int ARKodeSetDefaults(void* arkode_mem) ark_mem->maxncf = MAXNCF; /* max convergence fails */ ark_mem->maxconstrfails = MAXCONSTRFAILS; /* max number of constraint fails */ ark_mem->preallocated = SUNFALSE; /* data was not preallocated */ - ark_mem->hin = ZERO; /* determine initial step on-the-fly */ - ark_mem->hmin = ZERO; /* no minimum step size */ - ark_mem->hmax_inv = ZERO; /* no maximum step size */ - ark_mem->tstopset = SUNFALSE; /* no stop time set */ - ark_mem->tstopinterp = SUNFALSE; /* copy at stop time */ - ark_mem->tstop = ZERO; /* no fixed stop time */ + ark_mem->hin = ZERO; /* determine initial step on-the-fly */ + ark_mem->hmin = ZERO; /* no minimum step size */ + ark_mem->hmax_inv = ZERO; /* no maximum step size */ + ark_mem->tstopset = SUNFALSE; /* no stop time set */ + ark_mem->tstopinterp = SUNFALSE; /* copy at stop time */ + ark_mem->tstoplimited = SUNFALSE; /* tstop did not limit last step */ + ark_mem->skipadapttstop = SUNFALSE; /* tstop-limited steps can affect adaptivity */ + ark_mem->tstop = ZERO; /* no fixed stop time */ ark_mem->hadapt_mem->etamx1 = ETAMX1; /* max change on first step */ ark_mem->hadapt_mem->etamxf = ETAMXF; /* max change on error-failed step */ ark_mem->hadapt_mem->etamin = ETAMIN; /* min bound on time step reduction */ @@ -1272,6 +1274,27 @@ int ARKodeClearStopTime(void* arkode_mem) return (ARK_SUCCESS); } +/*--------------------------------------------------------------- + ARKodeSkipAdaptStopTime: + + Specifies whether stop-time-limited steps should be disregarded + when selecting step sizes for time step adaptivity. + ---------------------------------------------------------------*/ +int ARKodeSkipAdaptStopTime(void* arkode_mem, sunbooleantype skip) +{ + ARKodeMem ark_mem; + if (arkode_mem == NULL) + { + arkProcessError(NULL, ARK_MEM_NULL, __LINE__, __func__, __FILE__, + MSG_ARK_NO_MEM); + return (ARK_MEM_NULL); + } + ark_mem = (ARKodeMem)arkode_mem; + ark_mem->skipadapttstop = skip; + + return (ARK_SUCCESS); +} + /*--------------------------------------------------------------- ARKodeSetFixedStep: @@ -3494,6 +3517,8 @@ int ARKodeWriteParameters(void* arkode_mem, FILE* fp) { (void)SUNAdaptController_Write(ark_mem->hadapt_mem->hcontroller, fp); } + fprintf(fp, " Skip tstop-limited steps from affecting temporal adaptivity = %i\n", + ark_mem->skipadapttstop); fprintf(fp, " Maximum number of error test failures = %i\n", ark_mem->maxnef); fprintf(fp, " Maximum number of convergence test failures = %i\n", diff --git a/src/arkode/fmod_int32/farkode_mod.c b/src/arkode/fmod_int32/farkode_mod.c index 16c35a0e87..916426739f 100644 --- a/src/arkode/fmod_int32/farkode_mod.c +++ b/src/arkode/fmod_int32/farkode_mod.c @@ -620,6 +620,20 @@ SWIGEXPORT int _wrap_FARKodeClearStopTime(void *farg1) { } +SWIGEXPORT int _wrap_FARKodeSkipAdaptStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)ARKodeSkipAdaptStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FARKodeSetFixedStep(void *farg1, double const *farg2) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/arkode/fmod_int32/farkode_mod.f90 b/src/arkode/fmod_int32/farkode_mod.f90 index 6869accd7a..462eeb08f5 100644 --- a/src/arkode/fmod_int32/farkode_mod.f90 +++ b/src/arkode/fmod_int32/farkode_mod.f90 @@ -147,6 +147,7 @@ module farkode_mod public :: FARKodeSetInterpolateStopTime public :: FARKodeSetStopTime public :: FARKodeClearStopTime + public :: FARKodeSkipAdaptStopTime public :: FARKodeSetFixedStep public :: FARKodeSetStepDirection public :: FARKodeSetUserData @@ -697,6 +698,15 @@ function swigc_FARKodeClearStopTime(farg1) & integer(C_INT) :: fresult end function +function swigc_FARKodeSkipAdaptStopTime(farg1, farg2) & +bind(C, name="_wrap_FARKodeSkipAdaptStopTime") & +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_FARKodeSetFixedStep(farg1, farg2) & bind(C, name="_wrap_FARKodeSetFixedStep") & result(fresult) @@ -2935,6 +2945,22 @@ function FARKodeClearStopTime(arkode_mem) & swig_result = fresult end function +function FARKodeSkipAdaptStopTime(arkode_mem, skipadapttstop) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_INT), intent(in) :: skipadapttstop +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = arkode_mem +farg2 = skipadapttstop +fresult = swigc_FARKodeSkipAdaptStopTime(farg1, farg2) +swig_result = fresult +end function + function FARKodeSetFixedStep(arkode_mem, hfixed) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/src/arkode/fmod_int64/farkode_mod.c b/src/arkode/fmod_int64/farkode_mod.c index 904917178c..9e1bb6e0dc 100644 --- a/src/arkode/fmod_int64/farkode_mod.c +++ b/src/arkode/fmod_int64/farkode_mod.c @@ -620,6 +620,20 @@ SWIGEXPORT int _wrap_FARKodeClearStopTime(void *farg1) { } +SWIGEXPORT int _wrap_FARKodeSkipAdaptStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)ARKodeSkipAdaptStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FARKodeSetFixedStep(void *farg1, double const *farg2) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/arkode/fmod_int64/farkode_mod.f90 b/src/arkode/fmod_int64/farkode_mod.f90 index 0eb2e94994..1dc2368461 100644 --- a/src/arkode/fmod_int64/farkode_mod.f90 +++ b/src/arkode/fmod_int64/farkode_mod.f90 @@ -147,6 +147,7 @@ module farkode_mod public :: FARKodeSetInterpolateStopTime public :: FARKodeSetStopTime public :: FARKodeClearStopTime + public :: FARKodeSkipAdaptStopTime public :: FARKodeSetFixedStep public :: FARKodeSetStepDirection public :: FARKodeSetUserData @@ -697,6 +698,15 @@ function swigc_FARKodeClearStopTime(farg1) & integer(C_INT) :: fresult end function +function swigc_FARKodeSkipAdaptStopTime(farg1, farg2) & +bind(C, name="_wrap_FARKodeSkipAdaptStopTime") & +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_FARKodeSetFixedStep(farg1, farg2) & bind(C, name="_wrap_FARKodeSetFixedStep") & result(fresult) @@ -2935,6 +2945,22 @@ function FARKodeClearStopTime(arkode_mem) & swig_result = fresult end function +function FARKodeSkipAdaptStopTime(arkode_mem, skipadapttstop) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_INT), intent(in) :: skipadapttstop +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = arkode_mem +farg2 = skipadapttstop +fresult = swigc_FARKodeSkipAdaptStopTime(farg1, farg2) +swig_result = fresult +end function + function FARKodeSetFixedStep(arkode_mem, hfixed) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/src/cvode/cvode.c b/src/cvode/cvode.c index 7c07dd65ac..bb867238c7 100644 --- a/src/cvode/cvode.c +++ b/src/cvode/cvode.c @@ -306,6 +306,8 @@ void* CVodeCreate(int lmm, SUNContext sunctx) cv_mem->cv_small_nef = SMALL_NEF_DEFAULT; cv_mem->cv_tstopset = SUNFALSE; cv_mem->cv_tstopinterp = SUNFALSE; + cv_mem->cv_tstoplimited = SUNFALSE; + cv_mem->cv_skipadapttstop = SUNFALSE; cv_mem->cv_maxnef = MXNEF; cv_mem->cv_maxncf = MXNCF; cv_mem->cv_nlscoef = CORTES; @@ -1209,10 +1211,17 @@ int CVode(void* cvode_mem, sunrealtype tout, N_Vector yout, sunrealtype* tret, /* Check for approach to tstop */ + cv_mem->cv_tstoplimited = SUNFALSE; if (cv_mem->cv_tstopset) { if ((cv_mem->cv_tn + cv_mem->cv_h - cv_mem->cv_tstop) * cv_mem->cv_h > ZERO) { + if (cv_mem->cv_skipadapttstop) + { + cv_mem->cv_tstoplimited = SUNTRUE; + cv_mem->cv_hsave = cv_mem->cv_h; + cv_mem->cv_qsave = cv_mem->cv_q; + } cv_mem->cv_h = (cv_mem->cv_tstop - cv_mem->cv_tn) * (ONE - FOUR * cv_mem->cv_uround); } @@ -1367,6 +1376,7 @@ int CVode(void* cvode_mem, sunrealtype tout, N_Vector yout, sunrealtype* tret, } /* end of root stop check */ /* Test for tn at tstop or near tstop */ + cv_mem->cv_tstoplimited = SUNFALSE; if (cv_mem->cv_tstopset) { /* Test for tn at tstop */ @@ -1399,6 +1409,12 @@ int CVode(void* cvode_mem, sunrealtype tout, N_Vector yout, sunrealtype* tret, cv_mem->cv_h > ZERO) { + if (cv_mem->cv_skipadapttstop) + { + cv_mem->cv_tstoplimited = SUNTRUE; + cv_mem->cv_hsave = cv_mem->cv_hprime; + cv_mem->cv_qsave = cv_mem->cv_q; + } cv_mem->cv_hprime = (cv_mem->cv_tstop - cv_mem->cv_tn) * (ONE - FOUR * cv_mem->cv_uround); cv_mem->cv_eta = cv_mem->cv_hprime / cv_mem->cv_h; @@ -1593,6 +1609,7 @@ int CVode(void* cvode_mem, sunrealtype tout, N_Vector yout, sunrealtype* tret, } /* Check if tn is at tstop or near tstop */ + cv_mem->cv_tstoplimited = SUNFALSE; if (cv_mem->cv_tstopset) { troundoff = FUZZ_FACTOR * cv_mem->cv_uround * @@ -1621,6 +1638,12 @@ int CVode(void* cvode_mem, sunrealtype tout, N_Vector yout, sunrealtype* tret, cv_mem->cv_h > ZERO) { + if (cv_mem->cv_skipadapttstop) + { + cv_mem->cv_tstoplimited = SUNTRUE; + cv_mem->cv_hsave = cv_mem->cv_hprime; + cv_mem->cv_qsave = cv_mem->cv_q; + } cv_mem->cv_hprime = (cv_mem->cv_tstop - cv_mem->cv_tn) * (ONE - FOUR * cv_mem->cv_uround); cv_mem->cv_eta = cv_mem->cv_hprime / cv_mem->cv_h; @@ -3627,6 +3650,14 @@ static void cvPrepareNextStep(CVodeMem cv_mem, sunrealtype dsm) cv_mem->cv_hprime = cv_mem->cv_h; cv_mem->cv_eta = ONE; } + else if (cv_mem->cv_skipadapttstop && cv_mem->cv_tstoplimited) + { + /* If the current step was limited by tstop, set the upcoming step size + and order to match the values just prior to the tstop-limited step */ + cv_mem->cv_qprime = cv_mem->cv_qsave; + cv_mem->cv_hprime = cv_mem->cv_hsave; + cv_mem->cv_eta = cv_mem->cv_hsave / cv_mem->cv_h; + } else { /* etaq is the ratio of new to old h at the current order */ diff --git a/src/cvode/cvode_cli.c b/src/cvode/cvode_cli.c index cf42ef4080..53dcdf8a48 100644 --- a/src/cvode/cvode_cli.c +++ b/src/cvode/cvode_cli.c @@ -81,6 +81,7 @@ static int cvSetFromCommandLine(void* cvode_mem, const char* cvid, int argc, {"max_order", CVodeSetMaxOrd}, {"stab_lim_det", CVodeSetStabLimDet}, {"interpolate_stop_time", CVodeSetInterpolateStopTime}, + {"skip_adapt_stop_time", CVodeSkipAdaptStopTime}, {"use_integrator_fused_kernels", CVodeSetUseIntegratorFusedKernels}, {"num_fails_eta_max_err_fail", CVodeSetNumFailsEtaMaxErrFail}, {"linear_solution_scaling", CVodeSetLinearSolutionScaling}, diff --git a/src/cvode/cvode_impl.h b/src/cvode/cvode_impl.h index 73cf16b9b8..1f8ef6b657 100644 --- a/src/cvode/cvode_impl.h +++ b/src/cvode/cvode_impl.h @@ -250,6 +250,10 @@ typedef struct CVodeMemRec sunbooleantype cv_tstopset; sunbooleantype cv_tstopinterp; + sunbooleantype cv_tstoplimited; + sunbooleantype cv_skipadapttstop; + int cv_qsave; + sunrealtype cv_hsave; sunrealtype cv_tstop; /*--------- diff --git a/src/cvode/cvode_io.c b/src/cvode/cvode_io.c index 7a99daeb79..ebde9e47bb 100644 --- a/src/cvode/cvode_io.c +++ b/src/cvode/cvode_io.c @@ -736,6 +736,28 @@ int CVodeClearStopTime(void* cvode_mem) return (CV_SUCCESS); } +/* + * CVodeSkipAdaptStopTime + * + * Specifies whether stop-time-limited steps should be disregarded + * when performing step size and order adaptivity. + */ + +int CVodeSkipAdaptStopTime(void* cvode_mem, sunbooleantype skip) +{ + CVodeMem cv_mem; + + if (cvode_mem == NULL) + { + cvProcessError(NULL, CV_MEM_NULL, __LINE__, __func__, __FILE__, MSGCV_NO_MEM); + return (CV_MEM_NULL); + } + cv_mem = (CVodeMem)cvode_mem; + cv_mem->cv_skipadapttstop = skip; + + return (CV_SUCCESS); +} + /* * CVodeSetMaxErrTestFails * diff --git a/src/cvode/fmod_int32/fcvode_mod.c b/src/cvode/fmod_int32/fcvode_mod.c index 572c165c62..91ab09e508 100644 --- a/src/cvode/fmod_int32/fcvode_mod.c +++ b/src/cvode/fmod_int32/fcvode_mod.c @@ -663,6 +663,20 @@ SWIGEXPORT int _wrap_FCVodeClearStopTime(void *farg1) { } +SWIGEXPORT int _wrap_FCVodeSkipAdaptStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)CVodeSkipAdaptStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FCVodeSetUseIntegratorFusedKernels(void *farg1, int const *farg2) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/cvode/fmod_int32/fcvode_mod.f90 b/src/cvode/fmod_int32/fcvode_mod.f90 index fd03a4021d..36574792a5 100644 --- a/src/cvode/fmod_int32/fcvode_mod.f90 +++ b/src/cvode/fmod_int32/fcvode_mod.f90 @@ -95,6 +95,7 @@ module fcvode_mod public :: FCVodeSetStopTime public :: FCVodeSetInterpolateStopTime public :: FCVodeClearStopTime + public :: FCVodeSkipAdaptStopTime public :: FCVodeSetUseIntegratorFusedKernels public :: FCVodeSetUserData public :: FCVodeSetEtaFixedStepBounds @@ -484,6 +485,15 @@ function swigc_FCVodeClearStopTime(farg1) & integer(C_INT) :: fresult end function +function swigc_FCVodeSkipAdaptStopTime(farg1, farg2) & +bind(C, name="_wrap_FCVodeSkipAdaptStopTime") & +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_FCVodeSetUseIntegratorFusedKernels(farg1, farg2) & bind(C, name="_wrap_FCVodeSetUseIntegratorFusedKernels") & result(fresult) @@ -1884,6 +1894,22 @@ function FCVodeClearStopTime(cvode_mem) & swig_result = fresult end function +function FCVodeSkipAdaptStopTime(cvode_mem, skipadapttstop) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: cvode_mem +integer(C_INT), intent(in) :: skipadapttstop +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = cvode_mem +farg2 = skipadapttstop +fresult = swigc_FCVodeSkipAdaptStopTime(farg1, farg2) +swig_result = fresult +end function + function FCVodeSetUseIntegratorFusedKernels(cvode_mem, onoff) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/src/cvode/fmod_int64/fcvode_mod.c b/src/cvode/fmod_int64/fcvode_mod.c index 790f1aa937..c3d54c7f7f 100644 --- a/src/cvode/fmod_int64/fcvode_mod.c +++ b/src/cvode/fmod_int64/fcvode_mod.c @@ -663,6 +663,20 @@ SWIGEXPORT int _wrap_FCVodeClearStopTime(void *farg1) { } +SWIGEXPORT int _wrap_FCVodeSkipAdaptStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)CVodeSkipAdaptStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FCVodeSetUseIntegratorFusedKernels(void *farg1, int const *farg2) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/cvode/fmod_int64/fcvode_mod.f90 b/src/cvode/fmod_int64/fcvode_mod.f90 index e9cc5a6e03..55d34e7147 100644 --- a/src/cvode/fmod_int64/fcvode_mod.f90 +++ b/src/cvode/fmod_int64/fcvode_mod.f90 @@ -95,6 +95,7 @@ module fcvode_mod public :: FCVodeSetStopTime public :: FCVodeSetInterpolateStopTime public :: FCVodeClearStopTime + public :: FCVodeSkipAdaptStopTime public :: FCVodeSetUseIntegratorFusedKernels public :: FCVodeSetUserData public :: FCVodeSetEtaFixedStepBounds @@ -484,6 +485,15 @@ function swigc_FCVodeClearStopTime(farg1) & integer(C_INT) :: fresult end function +function swigc_FCVodeSkipAdaptStopTime(farg1, farg2) & +bind(C, name="_wrap_FCVodeSkipAdaptStopTime") & +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_FCVodeSetUseIntegratorFusedKernels(farg1, farg2) & bind(C, name="_wrap_FCVodeSetUseIntegratorFusedKernels") & result(fresult) @@ -1884,6 +1894,22 @@ function FCVodeClearStopTime(cvode_mem) & swig_result = fresult end function +function FCVodeSkipAdaptStopTime(cvode_mem, skipadapttstop) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: cvode_mem +integer(C_INT), intent(in) :: skipadapttstop +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = cvode_mem +farg2 = skipadapttstop +fresult = swigc_FCVodeSkipAdaptStopTime(farg1, farg2) +swig_result = fresult +end function + function FCVodeSetUseIntegratorFusedKernels(cvode_mem, onoff) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/src/cvodes/cvodes.c b/src/cvodes/cvodes.c index 1f188d72d7..0636cb14dc 100644 --- a/src/cvodes/cvodes.c +++ b/src/cvodes/cvodes.c @@ -544,6 +544,8 @@ void* CVodeCreate(int lmm, SUNContext sunctx) cv_mem->cv_small_nef = SMALL_NEF_DEFAULT; cv_mem->cv_tstopset = SUNFALSE; cv_mem->cv_tstopinterp = SUNFALSE; + cv_mem->cv_tstoplimited = SUNFALSE; + cv_mem->cv_skipadapttstop = SUNFALSE; cv_mem->cv_maxnef = MXNEF; cv_mem->cv_maxncf = MXNCF; cv_mem->cv_nlscoef = CORTES; @@ -3143,10 +3145,17 @@ int CVode(void* cvode_mem, sunrealtype tout, N_Vector yout, sunrealtype* tret, /* Check for approach to tstop */ + cv_mem->cv_tstoplimited = SUNFALSE; if (cv_mem->cv_tstopset) { if ((cv_mem->cv_tn + cv_mem->cv_h - cv_mem->cv_tstop) * cv_mem->cv_h > ZERO) { + if (cv_mem->cv_skipadapttstop) + { + cv_mem->cv_tstoplimited = SUNTRUE; + cv_mem->cv_hsave = cv_mem->cv_h; + cv_mem->cv_qsave = cv_mem->cv_q; + } cv_mem->cv_h = (cv_mem->cv_tstop - cv_mem->cv_tn) * (ONE - FOUR * cv_mem->cv_uround); } @@ -3343,6 +3352,7 @@ int CVode(void* cvode_mem, sunrealtype tout, N_Vector yout, sunrealtype* tret, } /* end of root stop check */ /* Test for tn at tstop or near tstop */ + cv_mem->cv_tstoplimited = SUNFALSE; if (cv_mem->cv_tstopset) { /* Test for tn at tstop */ @@ -3375,6 +3385,12 @@ int CVode(void* cvode_mem, sunrealtype tout, N_Vector yout, sunrealtype* tret, cv_mem->cv_h > ZERO) { + if (cv_mem->cv_skipadapttstop) + { + cv_mem->cv_tstoplimited = SUNTRUE; + cv_mem->cv_hsave = cv_mem->cv_hprime; + cv_mem->cv_qsave = cv_mem->cv_q; + } cv_mem->cv_hprime = (cv_mem->cv_tstop - cv_mem->cv_tn) * (ONE - FOUR * cv_mem->cv_uround); cv_mem->cv_eta = cv_mem->cv_hprime / cv_mem->cv_h; @@ -3622,6 +3638,7 @@ int CVode(void* cvode_mem, sunrealtype tout, N_Vector yout, sunrealtype* tret, } /* Check if tn is at tstop or near tstop */ + cv_mem->cv_tstoplimited = SUNFALSE; if (cv_mem->cv_tstopset) { troundoff = FUZZ_FACTOR * cv_mem->cv_uround * @@ -3650,6 +3667,12 @@ int CVode(void* cvode_mem, sunrealtype tout, N_Vector yout, sunrealtype* tret, cv_mem->cv_h > ZERO) { + if (cv_mem->cv_skipadapttstop) + { + cv_mem->cv_tstoplimited = SUNTRUE; + cv_mem->cv_hsave = cv_mem->cv_hprime; + cv_mem->cv_qsave = cv_mem->cv_q; + } cv_mem->cv_hprime = (cv_mem->cv_tstop - cv_mem->cv_tn) * (ONE - FOUR * cv_mem->cv_uround); cv_mem->cv_eta = cv_mem->cv_hprime / cv_mem->cv_h; @@ -7928,6 +7951,14 @@ static void cvPrepareNextStep(CVodeMem cv_mem, sunrealtype dsm) cv_mem->cv_hprime = cv_mem->cv_h; cv_mem->cv_eta = ONE; } + else if (cv_mem->cv_skipadapttstop && cv_mem->cv_tstoplimited) + { + /* If the current step was limited by tstop, set the upcoming step size + and order to match the values just prior to the tstop-limited step */ + cv_mem->cv_qprime = cv_mem->cv_qsave; + cv_mem->cv_hprime = cv_mem->cv_hsave; + cv_mem->cv_eta = cv_mem->cv_hsave / cv_mem->cv_h; + } else { /* etaq is the ratio of new to old h at the current order */ diff --git a/src/cvodes/cvodes_cli.c b/src/cvodes/cvodes_cli.c index 9524ea6869..2051cec61b 100644 --- a/src/cvodes/cvodes_cli.c +++ b/src/cvodes/cvodes_cli.c @@ -81,6 +81,7 @@ static int cvSetFromCommandLine(void* cvode_mem, const char* cvid, int argc, {"max_order", CVodeSetMaxOrd}, {"stab_lim_det", CVodeSetStabLimDet}, {"interpolate_stop_time", CVodeSetInterpolateStopTime}, + {"skip_adapt_stop_time", CVodeSkipAdaptStopTime}, {"num_fails_eta_max_err_fail", CVodeSetNumFailsEtaMaxErrFail}, {"quad_err_con", CVodeSetQuadErrCon}, {"sens_err_con", CVodeSetSensErrCon}, diff --git a/src/cvodes/cvodes_impl.h b/src/cvodes/cvodes_impl.h index 5f26db2302..6baeba67d9 100644 --- a/src/cvodes/cvodes_impl.h +++ b/src/cvodes/cvodes_impl.h @@ -380,6 +380,10 @@ typedef struct CVodeMemRec sunbooleantype cv_tstopset; sunbooleantype cv_tstopinterp; + sunbooleantype cv_tstoplimited; + sunbooleantype cv_skipadapttstop; + int cv_qsave; + sunrealtype cv_hsave; sunrealtype cv_tstop; /*--------- diff --git a/src/cvodes/cvodes_io.c b/src/cvodes/cvodes_io.c index dd3792b77a..76a9962708 100644 --- a/src/cvodes/cvodes_io.c +++ b/src/cvodes/cvodes_io.c @@ -738,6 +738,28 @@ int CVodeClearStopTime(void* cvode_mem) return (CV_SUCCESS); } +/* + * CVodeSkipAdaptStopTime + * + * Specifies whether stop-time-limited steps should be disregarded + * when performing step size and order adaptivity. + */ + +int CVodeSkipAdaptStopTime(void* cvode_mem, sunbooleantype skip) +{ + CVodeMem cv_mem; + + if (cvode_mem == NULL) + { + cvProcessError(NULL, CV_MEM_NULL, __LINE__, __func__, __FILE__, MSGCV_NO_MEM); + return (CV_MEM_NULL); + } + cv_mem = (CVodeMem)cvode_mem; + cv_mem->cv_skipadapttstop = skip; + + return (CV_SUCCESS); +} + /* * CVodeSetMaxErrTestFails * diff --git a/src/cvodes/fmod_int32/fcvodes_mod.c b/src/cvodes/fmod_int32/fcvodes_mod.c index 0974c4ebeb..60ed1a1043 100644 --- a/src/cvodes/fmod_int32/fcvodes_mod.c +++ b/src/cvodes/fmod_int32/fcvodes_mod.c @@ -737,6 +737,20 @@ SWIGEXPORT int _wrap_FCVodeClearStopTime(void *farg1) { } +SWIGEXPORT int _wrap_FCVodeSkipAdaptStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)CVodeSkipAdaptStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FCVodeSetUserData(void *farg1, void *farg2) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/cvodes/fmod_int32/fcvodes_mod.f90 b/src/cvodes/fmod_int32/fcvodes_mod.f90 index 90e4ff52a7..d0c4f96b78 100644 --- a/src/cvodes/fmod_int32/fcvodes_mod.f90 +++ b/src/cvodes/fmod_int32/fcvodes_mod.f90 @@ -126,6 +126,7 @@ module fcvodes_mod public :: FCVodeSetStopTime public :: FCVodeSetInterpolateStopTime public :: FCVodeClearStopTime + public :: FCVodeSkipAdaptStopTime public :: FCVodeSetUserData public :: FCVodeSetEtaFixedStepBounds public :: FCVodeSetEtaMaxFirstStep @@ -660,6 +661,15 @@ function swigc_FCVodeClearStopTime(farg1) & integer(C_INT) :: fresult end function +function swigc_FCVodeSkipAdaptStopTime(farg1, farg2) & +bind(C, name="_wrap_FCVodeSkipAdaptStopTime") & +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_FCVodeSetUserData(farg1, farg2) & bind(C, name="_wrap_FCVodeSetUserData") & result(fresult) @@ -3302,6 +3312,22 @@ function FCVodeClearStopTime(cvode_mem) & swig_result = fresult end function +function FCVodeSkipAdaptStopTime(cvode_mem, skipadapttstop) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: cvode_mem +integer(C_INT), intent(in) :: skipadapttstop +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = cvode_mem +farg2 = skipadapttstop +fresult = swigc_FCVodeSkipAdaptStopTime(farg1, farg2) +swig_result = fresult +end function + function FCVodeSetUserData(cvode_mem, user_data) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/src/cvodes/fmod_int64/fcvodes_mod.c b/src/cvodes/fmod_int64/fcvodes_mod.c index 54a1a268f6..db85a4d738 100644 --- a/src/cvodes/fmod_int64/fcvodes_mod.c +++ b/src/cvodes/fmod_int64/fcvodes_mod.c @@ -737,6 +737,20 @@ SWIGEXPORT int _wrap_FCVodeClearStopTime(void *farg1) { } +SWIGEXPORT int _wrap_FCVodeSkipAdaptStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)CVodeSkipAdaptStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FCVodeSetUserData(void *farg1, void *farg2) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/cvodes/fmod_int64/fcvodes_mod.f90 b/src/cvodes/fmod_int64/fcvodes_mod.f90 index 2396671fb1..2d3b6261c0 100644 --- a/src/cvodes/fmod_int64/fcvodes_mod.f90 +++ b/src/cvodes/fmod_int64/fcvodes_mod.f90 @@ -126,6 +126,7 @@ module fcvodes_mod public :: FCVodeSetStopTime public :: FCVodeSetInterpolateStopTime public :: FCVodeClearStopTime + public :: FCVodeSkipAdaptStopTime public :: FCVodeSetUserData public :: FCVodeSetEtaFixedStepBounds public :: FCVodeSetEtaMaxFirstStep @@ -660,6 +661,15 @@ function swigc_FCVodeClearStopTime(farg1) & integer(C_INT) :: fresult end function +function swigc_FCVodeSkipAdaptStopTime(farg1, farg2) & +bind(C, name="_wrap_FCVodeSkipAdaptStopTime") & +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_FCVodeSetUserData(farg1, farg2) & bind(C, name="_wrap_FCVodeSetUserData") & result(fresult) @@ -3302,6 +3312,22 @@ function FCVodeClearStopTime(cvode_mem) & swig_result = fresult end function +function FCVodeSkipAdaptStopTime(cvode_mem, skipadapttstop) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: cvode_mem +integer(C_INT), intent(in) :: skipadapttstop +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = cvode_mem +farg2 = skipadapttstop +fresult = swigc_FCVodeSkipAdaptStopTime(farg1, farg2) +swig_result = fresult +end function + function FCVodeSetUserData(cvode_mem, user_data) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/src/ida/fmod_int32/fida_mod.c b/src/ida/fmod_int32/fida_mod.c index 06f332f7db..40a2c086b7 100644 --- a/src/ida/fmod_int32/fida_mod.c +++ b/src/ida/fmod_int32/fida_mod.c @@ -572,6 +572,20 @@ SWIGEXPORT int _wrap_FIDAClearStopTime(void *farg1) { } +SWIGEXPORT int _wrap_FIDASkipAdaptStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)IDASkipAdaptStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FIDASetMaxErrTestFails(void *farg1, int const *farg2) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/ida/fmod_int32/fida_mod.f90 b/src/ida/fmod_int32/fida_mod.f90 index 5731f02b26..e46f921fe8 100644 --- a/src/ida/fmod_int32/fida_mod.f90 +++ b/src/ida/fmod_int32/fida_mod.f90 @@ -88,6 +88,7 @@ module fida_mod public :: FIDASetMinStep public :: FIDASetStopTime public :: FIDAClearStopTime + public :: FIDASkipAdaptStopTime public :: FIDASetMaxErrTestFails public :: FIDASetSuppressAlg public :: FIDASetId @@ -402,6 +403,15 @@ function swigc_FIDAClearStopTime(farg1) & integer(C_INT) :: fresult end function +function swigc_FIDASkipAdaptStopTime(farg1, farg2) & +bind(C, name="_wrap_FIDASkipAdaptStopTime") & +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_FIDASetMaxErrTestFails(farg1, farg2) & bind(C, name="_wrap_FIDASetMaxErrTestFails") & result(fresult) @@ -1594,6 +1604,22 @@ function FIDAClearStopTime(ida_mem) & swig_result = fresult end function +function FIDASkipAdaptStopTime(ida_mem, skipadapttstop) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: ida_mem +integer(C_INT), intent(in) :: skipadapttstop +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = ida_mem +farg2 = skipadapttstop +fresult = swigc_FIDASkipAdaptStopTime(farg1, farg2) +swig_result = fresult +end function + function FIDASetMaxErrTestFails(ida_mem, maxnef) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/src/ida/fmod_int64/fida_mod.c b/src/ida/fmod_int64/fida_mod.c index 1c3dc00058..cd1461897f 100644 --- a/src/ida/fmod_int64/fida_mod.c +++ b/src/ida/fmod_int64/fida_mod.c @@ -572,6 +572,20 @@ SWIGEXPORT int _wrap_FIDAClearStopTime(void *farg1) { } +SWIGEXPORT int _wrap_FIDASkipAdaptStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)IDASkipAdaptStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FIDASetMaxErrTestFails(void *farg1, int const *farg2) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/ida/fmod_int64/fida_mod.f90 b/src/ida/fmod_int64/fida_mod.f90 index 8a25821956..44761a711f 100644 --- a/src/ida/fmod_int64/fida_mod.f90 +++ b/src/ida/fmod_int64/fida_mod.f90 @@ -88,6 +88,7 @@ module fida_mod public :: FIDASetMinStep public :: FIDASetStopTime public :: FIDAClearStopTime + public :: FIDASkipAdaptStopTime public :: FIDASetMaxErrTestFails public :: FIDASetSuppressAlg public :: FIDASetId @@ -402,6 +403,15 @@ function swigc_FIDAClearStopTime(farg1) & integer(C_INT) :: fresult end function +function swigc_FIDASkipAdaptStopTime(farg1, farg2) & +bind(C, name="_wrap_FIDASkipAdaptStopTime") & +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_FIDASetMaxErrTestFails(farg1, farg2) & bind(C, name="_wrap_FIDASetMaxErrTestFails") & result(fresult) @@ -1594,6 +1604,22 @@ function FIDAClearStopTime(ida_mem) & swig_result = fresult end function +function FIDASkipAdaptStopTime(ida_mem, skipadapttstop) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: ida_mem +integer(C_INT), intent(in) :: skipadapttstop +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = ida_mem +farg2 = skipadapttstop +fresult = swigc_FIDASkipAdaptStopTime(farg1, farg2) +swig_result = fresult +end function + function FIDASetMaxErrTestFails(ida_mem, maxnef) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/src/ida/ida.c b/src/ida/ida.c index bed14c9140..83ad81ca9c 100644 --- a/src/ida/ida.c +++ b/src/ida/ida.c @@ -307,32 +307,34 @@ void* IDACreate(SUNContext sunctx) IDA_mem->ida_uround = SUN_UNIT_ROUNDOFF; /* Set default values for integrator optional inputs */ - IDA_mem->ida_res = NULL; - IDA_mem->ida_user_data = NULL; - IDA_mem->ida_itol = IDA_NN; - IDA_mem->ida_atolmin0 = SUNTRUE; - IDA_mem->ida_user_efun = SUNFALSE; - IDA_mem->ida_efun = NULL; - IDA_mem->ida_edata = NULL; - IDA_mem->ida_maxord = MAXORD_DEFAULT; - IDA_mem->ida_mxstep = MXSTEP_DEFAULT; - IDA_mem->ida_hmax_inv = HMAX_INV_DEFAULT; - IDA_mem->ida_hmin = HMIN_DEFAULT; - IDA_mem->ida_eta_max_fx = ETA_MAX_FX_DEFAULT; - IDA_mem->ida_eta_min_fx = ETA_MIN_FX_DEFAULT; - IDA_mem->ida_eta_max = ETA_MAX_DEFAULT; - IDA_mem->ida_eta_low = ETA_LOW_DEFAULT; - IDA_mem->ida_eta_min = ETA_MIN_DEFAULT; - IDA_mem->ida_eta_min_ef = ETA_MIN_EF_DEFAULT; - IDA_mem->ida_eta_cf = ETA_CF_DEFAULT; - IDA_mem->ida_hin = ZERO; - IDA_mem->ida_epcon = EPCON; - IDA_mem->ida_maxnef = MXNEF; - IDA_mem->ida_maxncf = MXNCF; - IDA_mem->ida_suppressalg = SUNFALSE; - IDA_mem->ida_id = NULL; - IDA_mem->ida_tstopset = SUNFALSE; - IDA_mem->ida_dcj = DCJ_DEFAULT; + IDA_mem->ida_res = NULL; + IDA_mem->ida_user_data = NULL; + IDA_mem->ida_itol = IDA_NN; + IDA_mem->ida_atolmin0 = SUNTRUE; + IDA_mem->ida_user_efun = SUNFALSE; + IDA_mem->ida_efun = NULL; + IDA_mem->ida_edata = NULL; + IDA_mem->ida_maxord = MAXORD_DEFAULT; + IDA_mem->ida_mxstep = MXSTEP_DEFAULT; + IDA_mem->ida_hmax_inv = HMAX_INV_DEFAULT; + IDA_mem->ida_hmin = HMIN_DEFAULT; + IDA_mem->ida_eta_max_fx = ETA_MAX_FX_DEFAULT; + IDA_mem->ida_eta_min_fx = ETA_MIN_FX_DEFAULT; + IDA_mem->ida_eta_max = ETA_MAX_DEFAULT; + IDA_mem->ida_eta_low = ETA_LOW_DEFAULT; + IDA_mem->ida_eta_min = ETA_MIN_DEFAULT; + IDA_mem->ida_eta_min_ef = ETA_MIN_EF_DEFAULT; + IDA_mem->ida_eta_cf = ETA_CF_DEFAULT; + IDA_mem->ida_hin = ZERO; + IDA_mem->ida_epcon = EPCON; + IDA_mem->ida_maxnef = MXNEF; + IDA_mem->ida_maxncf = MXNCF; + IDA_mem->ida_suppressalg = SUNFALSE; + IDA_mem->ida_id = NULL; + IDA_mem->ida_tstopset = SUNFALSE; + IDA_mem->ida_tstoplimited = SUNFALSE; + IDA_mem->ida_skipadapttstop = SUNFALSE; + IDA_mem->ida_dcj = DCJ_DEFAULT; /* Initialize inequality constraint variables */ IDA_mem->ida_constraints = NULL; @@ -1204,6 +1206,7 @@ int IDASolve(void* ida_mem, sunrealtype tout, sunrealtype* tret, N_Vector yret, /* Check for approach to tstop */ + IDA_mem->ida_tstoplimited = SUNFALSE; if (IDA_mem->ida_tstopset) { if ((IDA_mem->ida_tstop - IDA_mem->ida_tn) * IDA_mem->ida_hh <= ZERO) @@ -1217,6 +1220,12 @@ int IDASolve(void* ida_mem, sunrealtype tout, sunrealtype* tret, N_Vector yret, IDA_mem->ida_hh > ZERO) { + if (IDA_mem->ida_skipadapttstop) + { + IDA_mem->ida_tstoplimited = SUNTRUE; + IDA_mem->ida_hsave = IDA_mem->ida_hh; + IDA_mem->ida_ksave = IDA_mem->ida_kk; + } IDA_mem->ida_hh = (IDA_mem->ida_tstop - IDA_mem->ida_tn) * (ONE - FOUR * IDA_mem->ida_uround); } @@ -2182,6 +2191,7 @@ static int IDAStopTest1(IDAMem IDA_mem, sunrealtype tout, sunrealtype* tret, int ier; sunrealtype troundoff; + IDA_mem->ida_tstoplimited = SUNFALSE; if (IDA_mem->ida_tstopset) { /* Test for tn past tstop */ @@ -2219,6 +2229,12 @@ static int IDAStopTest1(IDAMem IDA_mem, sunrealtype tout, sunrealtype* tret, IDA_mem->ida_hh > ZERO) { + if (IDA_mem->ida_skipadapttstop) + { + IDA_mem->ida_tstoplimited = SUNTRUE; + IDA_mem->ida_hsave = IDA_mem->ida_hh; + IDA_mem->ida_ksave = IDA_mem->ida_kk; + } IDA_mem->ida_hh = (IDA_mem->ida_tstop - IDA_mem->ida_tn) * (ONE - FOUR * IDA_mem->ida_uround); } @@ -2292,6 +2308,7 @@ static int IDAStopTest2(IDAMem IDA_mem, sunrealtype tout, sunrealtype* tret, /* int ier; */ sunrealtype troundoff; + IDA_mem->ida_tstoplimited = SUNFALSE; if (IDA_mem->ida_tstopset) { troundoff = HUNDRED * IDA_mem->ida_uround * @@ -2315,6 +2332,12 @@ static int IDAStopTest2(IDAMem IDA_mem, sunrealtype tout, sunrealtype* tret, IDA_mem->ida_hh > ZERO) { + if (IDA_mem->ida_skipadapttstop) + { + IDA_mem->ida_tstoplimited = SUNTRUE; + IDA_mem->ida_hsave = IDA_mem->ida_hh; + IDA_mem->ida_ksave = IDA_mem->ida_kk; + } IDA_mem->ida_hh = (IDA_mem->ida_tstop - IDA_mem->ida_tn) * (ONE - FOUR * IDA_mem->ida_uround); } @@ -3270,7 +3293,11 @@ static void IDACompleteStep(IDAMem IDA_mem, sunrealtype err_k, sunrealtype err_k stepsize and order are set by the usual local error algorithm. Note that, after the first step, the order is not increased, as not all - of the necessary information is available yet. */ + of the necessary information is available yet. + + Also, if the current step size was reduced due to reaching a requested stop + time, reset the next step order to equal the order before the stop-time-reduced + step. */ if (IDA_mem->ida_phase == 0) { @@ -3288,7 +3315,12 @@ static void IDACompleteStep(IDAMem IDA_mem, sunrealtype err_k, sunrealtype err_k /* Set action = LOWER/MAINTAIN/RAISE to specify order decision */ - if (IDA_mem->ida_knew == IDA_mem->ida_kk - 1) + if (IDA_mem->ida_skipadapttstop && IDA_mem->ida_tstoplimited) + { + /* Returning order to the value used prior to the stop-time-reduced step. */ + action = MAINTAIN; + } + else if (IDA_mem->ida_knew == IDA_mem->ida_kk - 1) { /* Already decided to reduce the order */ action = LOWER; @@ -3355,6 +3387,9 @@ static void IDACompleteStep(IDAMem IDA_mem, sunrealtype err_k, sunrealtype err_k else { err_knew = err_k; } /* Compute tmp = tentative ratio hnew/hh from error norm estimate. + 0. If the current step size was reduced due to reaching a requested + stop time, set eta to return the step size to the size just preceding + the stop-time-reduced step. 1. If eta >= eta_max_fx (default = 2), increase hh to at most eta_max (default = 2) i.e., double the step size 2. If eta <= eta_min_fx (default = 1), reduce hh to between eta_min @@ -3364,7 +3399,12 @@ static void IDACompleteStep(IDAMem IDA_mem, sunrealtype err_k, sunrealtype err_k IDA_mem->ida_eta = ONE; tmp = SUNRpowerR(TWO * err_knew + PT0001, -ONE / (IDA_mem->ida_kk + 1)); - if (tmp >= IDA_mem->ida_eta_max_fx) + if (IDA_mem->ida_skipadapttstop && IDA_mem->ida_tstoplimited) + { + /* Returning step size to the value used prior to the stop-time-reduced step. */ + IDA_mem->ida_eta = IDA_mem->ida_hsave / IDA_mem->ida_hh; + } + else if (tmp >= IDA_mem->ida_eta_max_fx) { /* Enforce max growth factor bound and max step size */ IDA_mem->ida_eta = SUNMIN(tmp, IDA_mem->ida_eta_max); diff --git a/src/ida/ida_cli.c b/src/ida/ida_cli.c index 06f4b8a74f..5603721248 100644 --- a/src/ida/ida_cli.c +++ b/src/ida/ida_cli.c @@ -84,7 +84,8 @@ static int idaSetFromCommandLine(void* ida_mem, const char* idaid, int argc, {"max_conv_fails", IDASetMaxConvFails}, {"max_nonlin_iters", IDASetMaxNonlinIters}, {"linear_solution_scaling", IDASetLinearSolutionScaling}, - {"max_num_constraint_fails", IDASetMaxNumConstraintFails}}; + {"max_num_constraint_fails", IDASetMaxNumConstraintFails}, + {"skip_adapt_stop_time", IDASkipAdaptStopTime}}; static const int num_int_keys = sizeof(int_pairs) / sizeof(*int_pairs); static const struct sunKeyLongPair long_pairs[] = { diff --git a/src/ida/ida_impl.h b/src/ida/ida_impl.h index 3215f345ad..07426fb2dd 100644 --- a/src/ida/ida_impl.h +++ b/src/ida/ida_impl.h @@ -171,6 +171,10 @@ typedef struct IDAMemRec /* Tstop information */ sunbooleantype ida_tstopset; + sunbooleantype ida_tstoplimited; + sunbooleantype ida_skipadapttstop; + int ida_ksave; + sunrealtype ida_hsave; sunrealtype ida_tstop; /* Step Data */ diff --git a/src/ida/ida_io.c b/src/ida/ida_io.c index 4c8f357365..6f40a10358 100644 --- a/src/ida/ida_io.c +++ b/src/ida/ida_io.c @@ -413,6 +413,23 @@ int IDAClearStopTime(void* ida_mem) /*-----------------------------------------------------------------*/ +int IDASkipAdaptStopTime(void* ida_mem, sunbooleantype skip) +{ + IDAMem IDA_mem; + + if (ida_mem == NULL) + { + IDAProcessError(NULL, IDA_MEM_NULL, __LINE__, __func__, __FILE__, MSG_NO_MEM); + return (IDA_MEM_NULL); + } + IDA_mem = (IDAMem)ida_mem; + IDA_mem->ida_skipadapttstop = skip; + + return (IDA_SUCCESS); +} + +/*-----------------------------------------------------------------*/ + int IDASetNonlinConvCoef(void* ida_mem, sunrealtype epcon) { IDAMem IDA_mem; diff --git a/src/idas/fmod_int32/fidas_mod.c b/src/idas/fmod_int32/fidas_mod.c index d2dc442bab..f6af554629 100644 --- a/src/idas/fmod_int32/fidas_mod.c +++ b/src/idas/fmod_int32/fidas_mod.c @@ -647,6 +647,20 @@ SWIGEXPORT int _wrap_FIDAClearStopTime(void *farg1) { } +SWIGEXPORT int _wrap_FIDASkipAdaptStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)IDASkipAdaptStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FIDASetMaxErrTestFails(void *farg1, int const *farg2) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/idas/fmod_int32/fidas_mod.f90 b/src/idas/fmod_int32/fidas_mod.f90 index 4b6f5c20c1..a5941db94d 100644 --- a/src/idas/fmod_int32/fidas_mod.f90 +++ b/src/idas/fmod_int32/fidas_mod.f90 @@ -113,6 +113,7 @@ module fidas_mod public :: FIDASetMinStep public :: FIDASetStopTime public :: FIDAClearStopTime + public :: FIDASkipAdaptStopTime public :: FIDASetMaxErrTestFails public :: FIDASetSuppressAlg public :: FIDASetId @@ -575,6 +576,15 @@ function swigc_FIDAClearStopTime(farg1) & integer(C_INT) :: fresult end function +function swigc_FIDASkipAdaptStopTime(farg1, farg2) & +bind(C, name="_wrap_FIDASkipAdaptStopTime") & +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_FIDASetMaxErrTestFails(farg1, farg2) & bind(C, name="_wrap_FIDASetMaxErrTestFails") & result(fresult) @@ -3054,6 +3064,22 @@ function FIDAClearStopTime(ida_mem) & swig_result = fresult end function +function FIDASkipAdaptStopTime(ida_mem, skipadapttstop) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: ida_mem +integer(C_INT), intent(in) :: skipadapttstop +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = ida_mem +farg2 = skipadapttstop +fresult = swigc_FIDASkipAdaptStopTime(farg1, farg2) +swig_result = fresult +end function + function FIDASetMaxErrTestFails(ida_mem, maxnef) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/src/idas/fmod_int64/fidas_mod.c b/src/idas/fmod_int64/fidas_mod.c index 6becf9f0c2..4e15d5e299 100644 --- a/src/idas/fmod_int64/fidas_mod.c +++ b/src/idas/fmod_int64/fidas_mod.c @@ -647,6 +647,20 @@ SWIGEXPORT int _wrap_FIDAClearStopTime(void *farg1) { } +SWIGEXPORT int _wrap_FIDASkipAdaptStopTime(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)IDASkipAdaptStopTime(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FIDASetMaxErrTestFails(void *farg1, int const *farg2) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/idas/fmod_int64/fidas_mod.f90 b/src/idas/fmod_int64/fidas_mod.f90 index ac98811f2b..3b5ac17cf9 100644 --- a/src/idas/fmod_int64/fidas_mod.f90 +++ b/src/idas/fmod_int64/fidas_mod.f90 @@ -113,6 +113,7 @@ module fidas_mod public :: FIDASetMinStep public :: FIDASetStopTime public :: FIDAClearStopTime + public :: FIDASkipAdaptStopTime public :: FIDASetMaxErrTestFails public :: FIDASetSuppressAlg public :: FIDASetId @@ -575,6 +576,15 @@ function swigc_FIDAClearStopTime(farg1) & integer(C_INT) :: fresult end function +function swigc_FIDASkipAdaptStopTime(farg1, farg2) & +bind(C, name="_wrap_FIDASkipAdaptStopTime") & +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_FIDASetMaxErrTestFails(farg1, farg2) & bind(C, name="_wrap_FIDASetMaxErrTestFails") & result(fresult) @@ -3054,6 +3064,22 @@ function FIDAClearStopTime(ida_mem) & swig_result = fresult end function +function FIDASkipAdaptStopTime(ida_mem, skipadapttstop) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: ida_mem +integer(C_INT), intent(in) :: skipadapttstop +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = ida_mem +farg2 = skipadapttstop +fresult = swigc_FIDASkipAdaptStopTime(farg1, farg2) +swig_result = fresult +end function + function FIDASetMaxErrTestFails(ida_mem, maxnef) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/src/idas/idas.c b/src/idas/idas.c index 429262f3e9..b259aba9c1 100644 --- a/src/idas/idas.c +++ b/src/idas/idas.c @@ -437,32 +437,34 @@ void* IDACreate(SUNContext sunctx) IDA_mem->ida_uround = SUN_UNIT_ROUNDOFF; /* Set default values for integrator optional inputs */ - IDA_mem->ida_res = NULL; - IDA_mem->ida_user_data = NULL; - IDA_mem->ida_itol = IDA_NN; - IDA_mem->ida_atolmin0 = SUNTRUE; - IDA_mem->ida_user_efun = SUNFALSE; - IDA_mem->ida_efun = NULL; - IDA_mem->ida_edata = NULL; - IDA_mem->ida_maxord = MAXORD_DEFAULT; - IDA_mem->ida_mxstep = MXSTEP_DEFAULT; - IDA_mem->ida_hmax_inv = HMAX_INV_DEFAULT; - IDA_mem->ida_hmin = HMIN_DEFAULT; - IDA_mem->ida_eta_max_fx = ETA_MAX_FX_DEFAULT; - IDA_mem->ida_eta_min_fx = ETA_MIN_FX_DEFAULT; - IDA_mem->ida_eta_max = ETA_MAX_DEFAULT; - IDA_mem->ida_eta_low = ETA_LOW_DEFAULT; - IDA_mem->ida_eta_min = ETA_MIN_DEFAULT; - IDA_mem->ida_eta_min_ef = ETA_MIN_EF_DEFAULT; - IDA_mem->ida_eta_cf = ETA_CF_DEFAULT; - IDA_mem->ida_hin = ZERO; - IDA_mem->ida_epcon = EPCON; - IDA_mem->ida_maxnef = MXNEF; - IDA_mem->ida_maxncf = MXNCF; - IDA_mem->ida_suppressalg = SUNFALSE; - IDA_mem->ida_id = NULL; - IDA_mem->ida_tstopset = SUNFALSE; - IDA_mem->ida_dcj = DCJ_DEFAULT; + IDA_mem->ida_res = NULL; + IDA_mem->ida_user_data = NULL; + IDA_mem->ida_itol = IDA_NN; + IDA_mem->ida_atolmin0 = SUNTRUE; + IDA_mem->ida_user_efun = SUNFALSE; + IDA_mem->ida_efun = NULL; + IDA_mem->ida_edata = NULL; + IDA_mem->ida_maxord = MAXORD_DEFAULT; + IDA_mem->ida_mxstep = MXSTEP_DEFAULT; + IDA_mem->ida_hmax_inv = HMAX_INV_DEFAULT; + IDA_mem->ida_hmin = HMIN_DEFAULT; + IDA_mem->ida_eta_max_fx = ETA_MAX_FX_DEFAULT; + IDA_mem->ida_eta_min_fx = ETA_MIN_FX_DEFAULT; + IDA_mem->ida_eta_max = ETA_MAX_DEFAULT; + IDA_mem->ida_eta_low = ETA_LOW_DEFAULT; + IDA_mem->ida_eta_min = ETA_MIN_DEFAULT; + IDA_mem->ida_eta_min_ef = ETA_MIN_EF_DEFAULT; + IDA_mem->ida_eta_cf = ETA_CF_DEFAULT; + IDA_mem->ida_hin = ZERO; + IDA_mem->ida_epcon = EPCON; + IDA_mem->ida_maxnef = MXNEF; + IDA_mem->ida_maxncf = MXNCF; + IDA_mem->ida_suppressalg = SUNFALSE; + IDA_mem->ida_id = NULL; + IDA_mem->ida_tstopset = SUNFALSE; + IDA_mem->ida_tstoplimited = SUNFALSE; + IDA_mem->ida_skipadapttstop = SUNFALSE; + IDA_mem->ida_dcj = DCJ_DEFAULT; /* Initialize inequality constraint variables */ IDA_mem->ida_constraints = NULL; @@ -2692,6 +2694,7 @@ int IDASolve(void* ida_mem, sunrealtype tout, sunrealtype* tret, N_Vector yret, /* Check for approach to tstop */ + IDA_mem->ida_tstoplimited = SUNFALSE; if (IDA_mem->ida_tstopset) { if ((IDA_mem->ida_tstop - IDA_mem->ida_tn) * IDA_mem->ida_hh <= ZERO) @@ -2705,6 +2708,12 @@ int IDASolve(void* ida_mem, sunrealtype tout, sunrealtype* tret, N_Vector yret, IDA_mem->ida_hh > ZERO) { + if (IDA_mem->ida_skipadapttstop) + { + IDA_mem->ida_tstoplimited = SUNTRUE; + IDA_mem->ida_hsave = IDA_mem->ida_hh; + IDA_mem->ida_ksave = IDA_mem->ida_kk; + } IDA_mem->ida_hh = (IDA_mem->ida_tstop - IDA_mem->ida_tn) * (ONE - FOUR * IDA_mem->ida_uround); } @@ -5557,6 +5566,7 @@ static int IDAStopTest1(IDAMem IDA_mem, sunrealtype tout, sunrealtype* tret, int ier; sunrealtype troundoff; + IDA_mem->ida_tstoplimited = SUNFALSE; if (IDA_mem->ida_tstopset) { /* Test for tn past tstop */ @@ -5594,6 +5604,12 @@ static int IDAStopTest1(IDAMem IDA_mem, sunrealtype tout, sunrealtype* tret, IDA_mem->ida_hh > ZERO) { + if (IDA_mem->ida_skipadapttstop) + { + IDA_mem->ida_tstoplimited = SUNTRUE; + IDA_mem->ida_hsave = IDA_mem->ida_hh; + IDA_mem->ida_ksave = IDA_mem->ida_kk; + } IDA_mem->ida_hh = (IDA_mem->ida_tstop - IDA_mem->ida_tn) * (ONE - FOUR * IDA_mem->ida_uround); } @@ -5667,6 +5683,7 @@ static int IDAStopTest2(IDAMem IDA_mem, sunrealtype tout, sunrealtype* tret, /* int ier; */ sunrealtype troundoff; + IDA_mem->ida_tstoplimited = SUNFALSE; if (IDA_mem->ida_tstopset) { troundoff = HUNDRED * IDA_mem->ida_uround * @@ -5690,6 +5707,12 @@ static int IDAStopTest2(IDAMem IDA_mem, sunrealtype tout, sunrealtype* tret, IDA_mem->ida_hh > ZERO) { + if (IDA_mem->ida_skipadapttstop) + { + IDA_mem->ida_tstoplimited = SUNTRUE; + IDA_mem->ida_hsave = IDA_mem->ida_hh; + IDA_mem->ida_ksave = IDA_mem->ida_kk; + } IDA_mem->ida_hh = (IDA_mem->ida_tstop - IDA_mem->ida_tn) * (ONE - FOUR * IDA_mem->ida_uround); } @@ -7499,7 +7522,11 @@ static void IDACompleteStep(IDAMem IDA_mem, sunrealtype err_k, sunrealtype err_k stepsize and order are set by the usual local error algorithm. Note that, after the first step, the order is not increased, as not all - of the necessary information is available yet. */ + of the necessary information is available yet. + + Also, if the current step size was reduced due to reaching a requested stop + time, reset the next step order to equal the order before the stop-time-reduced + step. */ if (IDA_mem->ida_phase == 0) { @@ -7517,7 +7544,12 @@ static void IDACompleteStep(IDAMem IDA_mem, sunrealtype err_k, sunrealtype err_k /* Set action = LOWER/MAINTAIN/RAISE to specify order decision */ - if (IDA_mem->ida_knew == IDA_mem->ida_kk - 1) + if (IDA_mem->ida_skipadapttstop && IDA_mem->ida_tstoplimited) + { + /* Returning order to the value used prior to the stop-time-reduced step. */ + action = MAINTAIN; + } + else if (IDA_mem->ida_knew == IDA_mem->ida_kk - 1) { /* Already decided to reduce the order */ action = LOWER; @@ -7617,6 +7649,9 @@ static void IDACompleteStep(IDAMem IDA_mem, sunrealtype err_k, sunrealtype err_k else { err_knew = err_k; } /* Compute tmp = tentative ratio hnew/hh from error norm estimate. + 0. If the current step size was reduced due to reaching a requested + stop time, set eta to return the step size to the size just preceding + the stop-time-reduced step. 1. If eta >= eta_max_fx (default = 2), increase hh to at most eta_max (default = 2) i.e., double the step size 2. If eta <= eta_min_fx (default = 1), reduce hh to between eta_min @@ -7626,7 +7661,12 @@ static void IDACompleteStep(IDAMem IDA_mem, sunrealtype err_k, sunrealtype err_k IDA_mem->ida_eta = ONE; tmp = SUNRpowerR(TWO * err_knew + PT0001, -ONE / (IDA_mem->ida_kk + 1)); - if (tmp >= IDA_mem->ida_eta_max_fx) + if (IDA_mem->ida_skipadapttstop && IDA_mem->ida_tstoplimited) + { + /* Returning step size to the value used prior to the stop-time-reduced step. */ + IDA_mem->ida_eta = IDA_mem->ida_hsave / IDA_mem->ida_hh; + } + else if (tmp >= IDA_mem->ida_eta_max_fx) { /* Enforce max growth factor bound and max step size */ IDA_mem->ida_eta = SUNMIN(tmp, IDA_mem->ida_eta_max); diff --git a/src/idas/idas_cli.c b/src/idas/idas_cli.c index ddee644944..48ac5f41a9 100644 --- a/src/idas/idas_cli.c +++ b/src/idas/idas_cli.c @@ -88,7 +88,8 @@ static int idaSetFromCommandLine(void* ida_mem, const char* idaid, int argc, {"sens_max_nonlin_iters", IDASetSensMaxNonlinIters}, {"quad_sens_err_con", IDASetQuadSensErrCon}, {"linear_solution_scaling", IDASetLinearSolutionScaling}, - {"max_num_constraint_fails", IDASetMaxNumConstraintFails}}; + {"max_num_constraint_fails", IDASetMaxNumConstraintFails}, + {"skip_adapt_stop_time", IDASkipAdaptStopTime}}; static const int num_int_keys = sizeof(int_pairs) / sizeof(*int_pairs); static const struct sunKeyLongPair long_pairs[] = { diff --git a/src/idas/idas_impl.h b/src/idas/idas_impl.h index 2853c5cb02..d12290e79a 100644 --- a/src/idas/idas_impl.h +++ b/src/idas/idas_impl.h @@ -295,6 +295,10 @@ typedef struct IDAMemRec /* Tstop information */ sunbooleantype ida_tstopset; + sunbooleantype ida_tstoplimited; + sunbooleantype ida_skipadapttstop; + int ida_ksave; + sunrealtype ida_hsave; sunrealtype ida_tstop; /* Step Data */ diff --git a/src/idas/idas_io.c b/src/idas/idas_io.c index dce0a076f0..5b2981dd42 100644 --- a/src/idas/idas_io.c +++ b/src/idas/idas_io.c @@ -412,6 +412,23 @@ int IDAClearStopTime(void* ida_mem) /*-----------------------------------------------------------------*/ +int IDASkipAdaptStopTime(void* ida_mem, sunbooleantype skip) +{ + IDAMem IDA_mem; + + if (ida_mem == NULL) + { + IDAProcessError(NULL, IDA_MEM_NULL, __LINE__, __func__, __FILE__, MSG_NO_MEM); + return (IDA_MEM_NULL); + } + IDA_mem = (IDAMem)ida_mem; + IDA_mem->ida_skipadapttstop = skip; + + return (IDA_SUCCESS); +} + +/*-----------------------------------------------------------------*/ + int IDASetNonlinConvCoef(void* ida_mem, sunrealtype epcon) { IDAMem IDA_mem; diff --git a/test/answers b/test/answers index 49ef12abaf..9dde9a1354 160000 --- a/test/answers +++ b/test/answers @@ -1 +1 @@ -Subproject commit 49ef12abaf84ae79a8d6f6c1abea53a06e0e9771 +Subproject commit 9dde9a135477860c2f7c8a98d0f77b1a983434d5