diff --git a/CHANGELOG.md b/CHANGELOG.md index bd86096342..e363b3ca68 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -80,6 +80,11 @@ provided. Fixed the loading of ARKStep's default first order explicit method. +Fixed a bug in ARKODE when enabling rootfinding with fixed step sizes and the +initial value of the rootfinding function is zero. In this case, uninitialized +right-hand side data was used to compute a state value near the initial +condition to determine if any rootfinding functions are initially active. + Fixed a CMake bug regarding usage of missing "print_warning" macro that was only triggered when the deprecated `CUDA_ARCH` option was used. diff --git a/doc/shared/RecentChanges.rst b/doc/shared/RecentChanges.rst index c4cac39aa7..8aa4ea07d0 100644 --- a/doc/shared/RecentChanges.rst +++ b/doc/shared/RecentChanges.rst @@ -90,6 +90,11 @@ of ``0`` was provided. Fixed the loading of ARKStep's default first order explicit method. +Fixed a bug in ARKODE when enabling rootfinding with fixed step sizes and the +initial value of the rootfinding function is zero. In this case, uninitialized +right-hand side data was used to compute a state value near the initial +condition to determine if any rootfinding functions are initially active. + Fixed a CMake bug regarding usage of missing "print_warning" macro that was only triggered when the deprecated ``CUDA_ARCH`` option was used. diff --git a/src/arkode/arkode.c b/src/arkode/arkode.c index 099e2e361c..1a952f0bc6 100644 --- a/src/arkode/arkode.c +++ b/src/arkode/arkode.c @@ -2138,13 +2138,7 @@ int arkInitialSetup(ARKodeMem ark_mem, sunrealtype tout) if (ark_mem->root_mem->nrtfn > 0) { retval = arkRootCheck1((void*)ark_mem); - - if (retval == ARK_RTFUNC_FAIL) - { - arkProcessError(ark_mem, ARK_RTFUNC_FAIL, __LINE__, __func__, __FILE__, - MSG_ARK_RTFUNC_FAILED, ark_mem->tcur); - return (ARK_RTFUNC_FAIL); - } + if (retval != ARK_SUCCESS) { return retval; } } } diff --git a/src/arkode/arkode_root.c b/src/arkode/arkode_root.c index 05d0a4a08c..447830ec1e 100644 --- a/src/arkode/arkode_root.c +++ b/src/arkode/arkode_root.c @@ -434,7 +434,12 @@ int arkRootCheck1(void* arkode_mem) retval = rootmem->gfun(rootmem->tlo, ark_mem->yn, rootmem->glo, rootmem->root_data); rootmem->nge = 1; - if (retval != 0) { return (ARK_RTFUNC_FAIL); } + if (retval != 0) + { + arkProcessError(ark_mem, ARK_RTFUNC_FAIL, __LINE__, __func__, __FILE__, + MSG_ARK_RTFUNC_FAILED, ark_mem->tcur); + return (ARK_RTFUNC_FAIL); + } zroot = SUNFALSE; for (i = 0; i < rootmem->nrtfn; i++) @@ -447,6 +452,20 @@ int arkRootCheck1(void* arkode_mem) } if (!zroot) { return (ARK_SUCCESS); } + /* call full RHS if needed */ + if (!(ark_mem->fn_is_current)) + { + retval = ark_mem->step_fullrhs(ark_mem, ark_mem->tn, ark_mem->yn, + ark_mem->fn, ARK_FULLRHS_START); + if (retval) + { + arkProcessError(ark_mem, ARK_RHSFUNC_FAIL, __LINE__, __func__, + __FILE__, MSG_ARK_RHSFUNC_FAILED, ark_mem->tcur); + return ARK_RHSFUNC_FAIL; + } + ark_mem->fn_is_current = SUNTRUE; + } + /* Some g_i is zero at t0; look at g at t0+(small increment). */ hratio = SUNMAX(rootmem->ttol / SUNRabs(ark_mem->h), TENTH); smallh = hratio * ark_mem->h; @@ -454,7 +473,12 @@ int arkRootCheck1(void* arkode_mem) N_VLinearSum(ONE, ark_mem->yn, smallh, ark_mem->fn, ark_mem->ycur); retval = rootmem->gfun(tplus, ark_mem->ycur, rootmem->ghi, rootmem->root_data); rootmem->nge++; - if (retval != 0) { return (ARK_RTFUNC_FAIL); } + if (retval != 0) + { + arkProcessError(ark_mem, ARK_RTFUNC_FAIL, __LINE__, __func__, __FILE__, + MSG_ARK_RTFUNC_FAILED, ark_mem->tcur); + return (ARK_RTFUNC_FAIL); + } /* We check now only the components of g which were exactly 0.0 at t0 * to see if we can 'activate' them. */