diff --git a/doc/arkode/guide/source/Mathematics.rst b/doc/arkode/guide/source/Mathematics.rst index d2e5ac0b06..7222bee571 100644 --- a/doc/arkode/guide/source/Mathematics.rst +++ b/doc/arkode/guide/source/Mathematics.rst @@ -641,9 +641,11 @@ algorithm for a single step: diagonally-implicit, or additive Runge--Kutta stage update, .. math:: - z_i - \theta_{i,i} h^S f^I(t_{n-1}+c_i^S h^S, z_i) = a_i. + z_i - \theta_{i,i} h^S f^I(t_{n,i}^S, z_i) = a_i. :label: MRI_implicit_solve + where :math:`t_{n,j}^S = t_{n-1}*h^S c^S_j`. + #. Set :math:`y_{n} = z_{s}`. #. If the method has an embedding, compute the embedded solution, @@ -682,19 +684,19 @@ MIS, MRI-GARK, and IMEX-MRI-GARK Methods The methods in IMEX-MRI-GARK family, which includes MIS and MRI-GARK methods, are defined by a vector of slow stage time abscissae, :math:`c^S \in -\mathbb{R}^{s+1}`, and a set of coupling tensors, -:math:`\Omega\in\mathbb{R}^{(s+1)\times(s+1)\times k}` and -:math:`\Gamma\in\mathbb{R}^{(s+1)\times(s+1)\times k}`, that specify the +\mathbb{R}^{s}`, and a set of coupling tensors, +:math:`\Omega\in\mathbb{R}^{(s+1)\times s \times k}` and +:math:`\Gamma\in\mathbb{R}^{(s+1)\times s \times k}`, that specify the slow-to-fast coupling for the explicit and implicit components, respectively. The fast stage IVPs, :eq:`MRI_fast_IVP`, are evolved over non-overlapping -intervals :math:`[t_{0,i},t_{F,i}] = [t_{n-1}+c_{i-1}h^S, t_{n-1}+c_ih^S]` with +intervals :math:`[t_{0,i},t_{F,i}] = [t_{n,i-1}^S, t_{n,i}^S]` with the initial condition :math:`v_{0,i}=z_{i-1}`. The fast IVP forcing function is given by .. math:: - r_i(t) = \frac{1}{\Delta c_i^S} \sum\limits_{j=1}^{i-1} \omega_{i,j}(\tau) f^E(t_{n,j}^I, z_j) - + \frac{1}{\Delta c_i^S} \sum\limits_{j=1}^i \gamma_{i,j}(\tau) f^I(t_{n,j}^I, z_j) + r_i(t) = \frac{1}{\Delta c_i^S} \sum\limits_{j=1}^{i-1} \omega_{i,j}(\tau) f^E(t_{n,j}^S, z_j) + + \frac{1}{\Delta c_i^S} \sum\limits_{j=1}^i \gamma_{i,j}(\tau) f^I(t_{n,j}^S, z_j) where :math:`\Delta c_i^S=\left(c^S_i - c^S_{i-1}\right)`, :math:`\tau = (t - t_{n,i-1}^S)/(h^S \Delta c_i^S)` is the normalized time, the coefficients @@ -721,7 +723,7 @@ stage is computed as :label: ARKODE_MRI_delta_c_zero Similarly, the embedded solution IVP, :eq:`MRI_embedding_fast_IVP`, is evolved -over the interval :math:`[\tilde{t}_{0},\tilde{t}_{F}] = [t_{n-1}+c_{s-1}h^S, t_{n}]` +over the interval :math:`[\tilde{t}_{0},\tilde{t}_{F}] = [t_{n,s-1}^S, t_{n}]` with the initial condition :math:`\tilde{v}_0=z_{s-1}`. As with standard ARK and DIRK methods, implicitness at the slow time scale is @@ -746,26 +748,87 @@ The IMEX-MRI-SR family of methods perform *both* the fast IVP evolution, :eq:`MRI_fast_IVP` or :eq:`MRI_embedding_fast_IVP`, *and* stage update, :eq:`MRI_implicit_solve` or :eq:`MRI_embedding_implicit_solve`, in every stage (but these methods typically have far fewer stages than implicit MRI-GARK or -IMEX-MRI-GARK methods). +IMEX-MRI-GARK methods). These methods are defined by a vector of slow stage +time abscissae :math:`c^S \in \mathbb{R}^{s}`, a set of coupling tensors +:math:`\Omega\in\mathbb{R}^{(s+1)\times s\times k}`, and a Butcher table of +slow-implicit coefficients, :math:`\Gamma\in\mathbb{R}^{(s+1) \times s}`. + +The fast stage IVPs, :eq:`MRI_fast_IVP`, are evolved on overlapping +intervals :math:`[t_{0,i},t_{F,i}] = [t_{n-1}, t_{n,i}^S]` with +the initial condition :math:`v_{0,i}=y_{n-1}`. The fast IVP forcing function is +given by + +.. math:: + r_i(t) = \frac{1}{c_i^S} \sum\limits_{j=1}^{i-1} \omega_{i,j}(\tau) \left( f^E(t_{n,j}^S, z_j) + f^I(t_{n,j}^S, z_j)\right), + :label: IMEXMRISR_forcing + +where :math:`\tau = (t - t_n)/(h^S c_i^S)` is the normalized time, and the coefficients +:math:`\omega_{i,j}` are polynomials in time of degree :math:`k` that are also given by +:eq:`ARKODE_MRI_coupling`. The solution of these fast IVPs defines an intermediate stage +solution, :math:`\tilde{z}_i`. + +The implicit solve that follows each fast IVP must solve the algebraic equation for :math:`z_i` + +.. math:: + z_i = \tilde{z}_i + h^S \sum_{j=1}^{i} \gamma_{i,j} f^I(t_{n,j}^S, z_j). + :label: ARKODE_MRISR_implicit + +We note that IMEX-MRI-SR methods are solve-decoupled by construction, and thus the structure +of a given stage never needs to be deduced based on :math:`\Delta c_i^S`. However, ARKODE +still checks the value of :math:`\gamma_{i,i}`, since if it zero then the stage update +equation :eq:`ARKODE_MRISR_implicit` simplifies to a simple explicit Runge--Kutta-like stage +update. + +The overall time step solution is given by the final internal stage solution, +i.e., :math:`y_{n} = z_{s}`. The embedded solution is computing using the above +algorithm for stage index :math:`s+1`, under the definition that :math:`c_{s+1}^S=1` +(and thus the fast IVP portion is evolved over the full time step, +:math:`[t_{0,i},t_{F,i}] = [t_{n-1}, t_{n}]`). -These methods evolve the fast IVPs, :eq:`MRI_fast_IVP` and -:eq:`MRI_embedding_fast_IVP`, on overlapping time intervals -:math:`[t_{0,i},t_{F,i}] = [t_{n-1}, t_{n-1}+c_i h^S]` and -:math:`[\tilde{t}_{0},\tilde{t}_{F}] = [t_{n-1}, t_{n}]`, using initial -conditions :math:`v_{0,i}=y_{n-1}` and :math:`\tilde{v}_0=y_{n-1}`. MERK Methods ------------ -MERK family of methods never include the stage updates, :eq:`MRI_implicit_solve` -or :eq:`MRI_embedding_implicit_solve`. +The MERK family of methods are only defined for multirate applications that +are explicit at the slow time scale, i.e., :math:`f^I=0`, but otherwise they are +nearly identical to IMEX-MRI-SR methods. Specifically, like IMEX-MRI-SR methods, +these evolve the fast IVPs +:eq:`MRI_fast_IVP` and :eq:`MRI_embedding_fast_IVP` over the intervals +:math:`[t_{0,i},t_{F,i}] = [t_{n-1}, t_{n,i}^S]` and +:math:`[t_{0,i},t_{F,i}] = [t_{n-1}, t_{n}]`, respectively, and begin +with the initial condition :math:`v_{0,i}=y_{n-1}`. Furthermore, the fast IVP +forcing functions are given by :eq:`IMEXMRISR_forcing` with :math:`f^I=0`. +As MERK-based applications lack the implicit slow operator, they do not require +the solution of implicit algebraic equations. + +However, unlike other MRI families, MERK methods were designed to admit a useful +efficiency improvement. Since each fast IVP begins with the same initial condition, +:math:`v_{0,i}=y_{n-1}`, if multiple stages share the same forcing function +:math:`r_i(t)`, then they may be "grouped" together. This is achieved by sorting the +final IVP solution time for each stage, :math:`t_{n,i}^S`, and then evolving the inner +solver to each of these stage times in order, storing the corresponding inner solver +solutions at these times as the stages :math:`z_i`. For example, the +:index:`ARKODE_MERK54` method involves 11 stages, that may be combined into 5 distinct +groups. The fourth group contains stages 7, 8, 9, and the embedding, corresponding to +the :math:`c_i^S` values :math:`7/10`, :math:`1/2`, :math:`2/3`, and :math:`1`. +Sorting these, a single fast IVP for this group must be evolved over the interval +:math:`[t_{0,i},t_{F,i}] = [t_{n-1}, t_{n}]`, first pausing at :math:`t_{n-1}+\frac12 h^S` +to store :math:`z_8`, then pausing at :math:`t_{n-1}+\frac{2}{3} h^S` to store +:math:`z_9`, then pausing at :math:`t_{n-1}+\frac{7}{10} h^S` to store :math:`z_7`, +and finally finishing the IVP solve to :math:`t_{n-1}+h^S` to obtain :math:`\tilde{y}_n`. + +.. note:: + + Although all MERK methods were derived in :cite:p:`Luan:20` under an assumption that + the fast function :math:`f^F(t,y)` is linear in :math:`y`, in :cite:p:`Fish:24` it + was proven that MERK methods also satisfy all nonlinear order conditions up through + their linear order. The lone exception is :index:`ARKODE_MERK54`, where it was only + proven to satisfy all nonlinear conditions up to order 4 (since :cite:p:`Fish:24` did + not establish the formulas for the order 5 conditions). All our numerical tests to + date have shown :index:`ARKODE_MERK54` to achieve fifth order for nonlinear problems, + and so we conjecture that it also satisfies the nonlinear fifth order conditions. -These methods evolve the fast IVPs, :eq:`MRI_fast_IVP` and -:eq:`MRI_embedding_fast_IVP`, on overlapping time intervals -:math:`[t_{0,i},t_{F,i}] = [t_{n-1}, t_{n-1}+c_i h^S]` and -:math:`[\tilde{t}_{0},\tilde{t}_{F}] = [t_{n-1}, t_{n}]`, using initial -conditions :math:`v_{0,i}=y_{n-1}` and :math:`\tilde{v}_0=y_{n-1}`. .. _ARKODE.Mathematics.Error.Norm: @@ -1063,17 +1126,25 @@ that it may have used a potentially different relative solution tolerance, estimates using either a "maximum accumulation" strategy, .. math:: - \varepsilon^F_{max} = \text{RTOL}^F \max_{m\in M} \|\varepsilon^F_{n,m}\|_{WRMS}, + \varepsilon^F_{max} = \text{RTOL}^F \max_{m\in \mathcal{S}} \|\varepsilon^F_{n,m}\|_{WRMS}, + :label: maximum_accumulation + +an "additive accumulation" strategy, + +.. math:: + \varepsilon^F_{sum} = \text{RTOL}^F \sum_{m\in \mathcal{S}} \|\varepsilon^F_{n,m}\|_{WRMS}, :label: maximum_accumulation or using an "averaged accumulation" strategy, .. math:: - \varepsilon^F_{avg} = \frac{\text{RTOL}^F}{|M|} \sum_{m\in M} \|\varepsilon^F_{n,m}\|_{WRMS}, + \varepsilon^F_{avg} = \frac{\text{RTOL}^F}{\Delta t_{\mathcal{S}}} \sum_{m\in \mathcal{S}} h_{n,m} \|\varepsilon^F_{n,m}\|_{WRMS}, :label: average_accumulation -where the norm is taken using the tolerance-informed error-weight vector. In -both cases, the sum or the maximum is taken over the set of all steps :math:`M` +where :math:`h_{n,m}` is the step size that gave rise to :math:`\varepsilon^F_{n,m}`, +:math:`\Delta t_{\mathcal{S}}` denotes the elapsed time over which :math:`\mathcal{S}` +is taken, and the norms are taken using the tolerance-informed error-weight vector. In +each case, the sum or the maximum is taken over the set of all steps :math:`\mathcal{S}` since the fast error accumulator has been reset. @@ -1152,7 +1223,7 @@ initial state, where we again satisfy the inequality with a value of :math:`\frac12` to obtain .. math:: - |h_0| = \frac{1}{2\left\| f(t_0,y_0)\right\|}, + |h_0| = \frac{1}{2\left\| f(t_0,y_0)\right\|}. :label: H0_TSExp0 @@ -1565,7 +1636,7 @@ but may be modified by the user. When using the dense and band SUNMatrix objects for the linear systems :eq:`ARKODE_modified_Newton_system`, the Jacobian :math:`J` may be supplied -by a user routine, or approximated internally by finite-differences. +by a user routine, or approximated internally with finite-differences. In the case of differencing, we use the standard approximation .. math:: diff --git a/doc/arkode/guide/source/Usage/MRIStep/Custom_Inner_Stepper/Description.rst b/doc/arkode/guide/source/Usage/MRIStep/Custom_Inner_Stepper/Description.rst index a2daa41e4b..4076cb5b0e 100644 --- a/doc/arkode/guide/source/Usage/MRIStep/Custom_Inner_Stepper/Description.rst +++ b/doc/arkode/guide/source/Usage/MRIStep/Custom_Inner_Stepper/Description.rst @@ -237,7 +237,7 @@ Setting Member Functions :retval ARK_SUCCESS: if successful :retval ARK_ILL_INPUT: if the stepper is ``NULL`` - + .. versionadded: x.y.z @@ -251,7 +251,7 @@ Setting Member Functions :retval ARK_SUCCESS: if successful :retval ARK_ILL_INPUT: if the stepper is ``NULL`` - + .. versionadded: x.y.z @@ -265,7 +265,7 @@ Setting Member Functions :retval ARK_SUCCESS: if successful :retval ARK_ILL_INPUT: if the stepper is ``NULL`` - + .. versionadded: x.y.z @@ -416,6 +416,7 @@ following member functions: **Example codes:** * ``examples/arkode/CXX_parallel/ark_diffusion_reaction_p.cpp`` + .. versionchanged:: v5.7.0 Supplying a full right-hand side function was made optional. @@ -424,7 +425,7 @@ following member functions: This function resets the inner (fast) stepper state to the provided independent variable value and dependent variable vector. - + If provided, the :c:type:`MRIStepInnerResetFn` function will be called *before* a call to :c:type:`MRIStepInnerEvolveFn` when the state was updated at the slow timescale. @@ -439,16 +440,17 @@ following member functions: value if a recoverable error occurred, or a negative value if it failed unrecoverably. - .. note:: - - **Example codes:** * ``examples/arkode/CXX_parallel/ark_diffusion_reaction_p.cpp`` .. c:type:: int (*MRIStepInnerGetAccumulatedError)(MRIStepInnerStepper stepper, sunrealtype* accum_error) - This function returns an estimate of the accumulated solution error arising from the inner stepper. + This function returns an estimate of the accumulated solution error arising from the + inner stepper. Both the :c:type:`MRIStepInnerGetAccumulatedError` and + :c:type:`MRIStepInnerResetAccumulatedError` functions should be provided, or not; if only + one is provided then MRIStep will disable multirate temporal adaptivity and call neither. + **Arguments:** * *stepper* -- the inner stepper object. @@ -461,20 +463,19 @@ following member functions: .. note:: This function is only called when multirate temporal adaptivity has been enabled, - using a :c:type:`SUNAdaptController` module having type :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_TOL`. + using a :c:type:`SUNAdaptController` module having type :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_H_TOL`. If provided, the :c:type:`MRIStepInnerGetAccumulatedError` function will always be called *after* a preceding call to the :c:type:`MRIStepInnerResetAccumulatedError` function. - Both the :c:type:`MRIStepInnerGetAccumulatedError` and - :c:type:`MRIStepInnerResetAccumulatedError` functions should be provided, or not; if only - one is provided then MRIStep will disable multirate temporal adaptivity and call neither. - .. c:type:: int (*MRIStepInnerResetAccumulatedError)(MRIStepInnerStepper stepper) This function resets the inner stepper's accumulated solution error to zero. + This function performs a different role within MRIStep than the + :c:type:`MRIStepInnerResetFn`, and thus an implementation should make no + assumptions about the frequency/ordering of calls to either. **Arguments:** * *stepper* -- the inner stepper object. @@ -486,7 +487,7 @@ following member functions: .. note:: This function is only called when multirate temporal adaptivity has been enabled, - using a :c:type:`SUNAdaptController` module having type :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_TOL`. + using a :c:type:`SUNAdaptController` module having type :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_H_TOL`. The :c:type:`MRIStepInnerResetAccumulatedError` function will always be called *before* any calls to the :c:type:`MRIStepInnerGetAccumulatedError` function. @@ -495,14 +496,13 @@ following member functions: :c:type:`MRIStepInnerResetAccumulatedError` functions should be provided, or not; if only one is provided then MRIStep will disable multirate temporal adaptivity and call neither. - This function performs a different role within MRIStep than the - :c:type:`MRIStepInnerResetFn`, and thus an implementation should make no - assumptions about the frequency/ordering of calls to either. - .. c:type:: int (*MRIStepInnerSetRTol)(MRIStepInnerStepper stepper, sunrealtype rtol) - This function accepts a relative tolerance for the inner stepper to use in its upcoming adaptive solve. + This function accepts a relative tolerance for the inner stepper to use in its + upcoming adaptive solve. It is assumed that if the inner stepper supports absolute + tolerances as well, then these have been set up directly by the user to indicate the + "noise" level for solution components. **Arguments:** * *stepper* -- the inner stepper object. @@ -515,8 +515,4 @@ following member functions: .. note:: This function is only called when multirate temporal adaptivity has been enabled - using a :c:type:`SUNAdaptController` module having type :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_TOL`. - - It is assumed that if the inner stepper supports absolute tolerances as well, then - these have been set up directly by the user to indicate the "noise" level for - solution components. + using a :c:type:`SUNAdaptController` module having type :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_H_TOL`. diff --git a/doc/arkode/guide/source/Usage/MRIStep/MRIStepCoupling.rst b/doc/arkode/guide/source/Usage/MRIStep/MRIStepCoupling.rst index bda7590bfc..46590c75a7 100644 --- a/doc/arkode/guide/source/Usage/MRIStep/MRIStepCoupling.rst +++ b/doc/arkode/guide/source/Usage/MRIStep/MRIStepCoupling.rst @@ -273,8 +273,6 @@ are defined ``arkode/arkode_mristep.h``. A ``NULL`` pointer if ``stages`` was invalid, an allocation error occurred, or the input data arrays are inconsistent with the method type. - .. note:: - .. c:function:: MRIStepCoupling MRIStepCoupling_MIStoMRI(ARKodeButcherTable B, int q, int p) @@ -425,19 +423,23 @@ with values specified for each method below (e.g., ``ARKODE_MIS_KW3``). Notes regarding the above table: -#. The default method for each order when using fixed step sizes is marked with an asterisk (:math:`^*`). +#. The default method for each order when using fixed step sizes is marked with an + asterisk (:math:`^*`). -#. The default method for each order when using adaptive time stepping is marked with a circle (:math:`^\circ`). +#. The default method for each order when using adaptive time stepping is marked + with a circle (:math:`^\circ`). -#. The "Slow RHS Calls" column corresponds to the number of calls to the slow right-hand side function, - :math:`f^E`, per time step. +#. The "Slow RHS Calls" column corresponds to the number of calls to the slow + right-hand side function, :math:`f^E`, per time step. -#. Note A: although all MERK methods were derived in :cite:p:`Luan:20` under an assumption that the fast time - scale is linear in the solution, in :cite:p:`Fish:24` it was proven that they also satisfy all nonlinear - order conditions up through their linear order. The lone exception is MERK54, where it was only proven to - satisfy all nonlinear conditions up to order 4, since :cite:p:`Fish:24` did not establish the formulas for - the order 5 conditions. All our numerical tests to date have shown MERK54 to achieve fifth order for - nonlinear problems, and so we conjecture that it also satisfies the nonlinear fifth order conditions. +#. Note A: although all MERK methods were derived in :cite:p:`Luan:20` under an assumption + that the fast function :math:`f^F(t,y)` is linear in :math:`y`, in :cite:p:`Fish:24` it + was proven that MERK methods also satisfy all nonlinear order conditions up through + their linear order. The lone exception is :index:`ARKODE_MERK54`, where it was only + proven to satisfy all nonlinear conditions up to order 4 (since :cite:p:`Fish:24` did + not establish the formulas for the order 5 conditions). All our numerical tests to + date have shown :index:`ARKODE_MERK54` to achieve fifth order for nonlinear problems, + and so we conjecture that it also satisfies the nonlinear fifth order conditions. .. table:: Diagonally-implicit, solve-decoupled MRI-GARK coupling tables. The default diff --git a/doc/arkode/guide/source/Usage/MRIStep/Skeleton.rst b/doc/arkode/guide/source/Usage/MRIStep/Skeleton.rst index 1a7bf46c42..c8e239086c 100644 --- a/doc/arkode/guide/source/Usage/MRIStep/Skeleton.rst +++ b/doc/arkode/guide/source/Usage/MRIStep/Skeleton.rst @@ -85,11 +85,16 @@ unchanged from the skeleton program presented in inputs to :c:func:`MRIStepCreate` is the ``MRIStepInnerStepper`` object for solving the fast (inner) IVP created in the previous step. -#. Set the slow step size - - Call :c:func:`ARKodeSetFixedStep()` on the MRIStep object to specify the +#. If using fixed step sizes, then set the slow step size by calling + :c:func:`ARKodeSetFixedStep()` on the MRIStep object to specify the slow time step size. + If using adaptive slow steps, then specify the desired integration tolerances + as normal. By default, MRIStep will use a "decoupled" (see + :numref:`ARKODE.Mathematics.MultirateControllers`) I controller (see + :numref:`SUNAdaptController.Soderlind`), Alternately, create and attach a + multirate temporal controller (see :numref:`SUNAdaptController.MRIHTol`). + #. Create and configure implicit solvers (*as appropriate*) Specifically, if MRIStep is configured with an implicit slow right-hand side diff --git a/doc/arkode/guide/source/Usage/User_callable.rst b/doc/arkode/guide/source/Usage/User_callable.rst index 5c394b0de1..e07b08d4e5 100644 --- a/doc/arkode/guide/source/Usage/User_callable.rst +++ b/doc/arkode/guide/source/Usage/User_callable.rst @@ -1452,7 +1452,7 @@ Reset accumulated error :c:func:`ARKodeReset Not all time-stepping modules are compatible with all types of :c:type:`SUNAdaptController` objects. While all steppers that support temporal adaptivity support controllers with :c:type:`SUNAdaptController_Type` type ``SUN_ADAPTCONTROLLER_H``, only MRIStep supports - inputs with type ``SUN_ADAPTCONTROLLER_MRI_TOL``. + inputs with type ``SUN_ADAPTCONTROLLER_MRI_H_TOL``. .. versionadded:: 6.1.0 @@ -1758,23 +1758,25 @@ Reset accumulated error :c:func:`ARKodeReset The following routines are used to control algorithms that ARKODE can use to estimate -the accumulated temporal error over multiple time steps. For time-stepping modules that -compute both a solution and embedding, :math:`y_n` and :math:`\tilde{y}_n`, these may be -combined to create a vector-valued local temporal error estimate for the current internal -step, :math:`y_n - \tilde{y}_n`. These local errors may be accumulated by ARKODE in a -variety of ways, as determined by the enumerated type :c:enum:`ARKAccumError`. In each -of the cases below, the accumulation is taken over all steps since the most recent -call to either :c:func:`ARKodeSetAccumulatedErrorType` or +the accumulated temporal error over multiple time steps. While these may be informational +for users on their applications, this functionality is required when using multirate +temporal adaptivity in MRIStep via the :ref:`SUNAdaptController_MRIHTol ` +module. For time-stepping modules that compute both a solution and embedding, :math:`y_n` +and :math:`\tilde{y}_n`, these may be combined to create a vector-valued local temporal error +estimate for the current internal step, :math:`y_n - \tilde{y}_n`. These local errors may be +accumulated by ARKODE in a variety of ways, as determined by the enumerated type +:c:enum:`ARKAccumError`. In each of the cases below, the accumulation is taken over all steps +since the most recent call to either :c:func:`ARKodeSetAccumulatedErrorType` or :c:func:`ARKodeResetAccumulatedError`. Below the set :math:`\mathcal{S}` contains the indices of the steps since the last call to either of the aforementioned functions. -The norm is taken using the tolerance-informed -error-weight vector (see :c:func:`ARKodeGetErrWeights`), and ``reltol`` is the -user-specified relative solution tolerance. +The norm is taken using the tolerance-informed error-weight vector (see +:c:func:`ARKodeGetErrWeights`), and ``reltol`` is the user-specified relative solution +tolerance. .. c:enum:: ARKAccumError The type of error accumulation that ARKODE should use. - + .. versionadded:: x.y.z .. c:enumerator:: ARK_ACCUMERROR_NONE @@ -1783,21 +1785,29 @@ user-specified relative solution tolerance. .. c:enumerator:: ARK_ACCUMERROR_MAX - Computes :math:`\text{reltol} \max_{i \in \mathcal{S}} \|y_i - \tilde{y}_i\|_{WRMS}` + Computes :math:`\text{reltol} \max\limits_{i \in \mathcal{S}} \|y_i - \tilde{y}_i\|_{WRMS}` .. c:enumerator:: ARK_ACCUMERROR_SUM - Computes :math:`\text{reltol} \sum_{i \in \mathcal{S}} \|y_i - \tilde{y}_i\|_{WRMS}` + Computes :math:`\text{reltol} \sum\limits_{i \in \mathcal{S}} \|y_i - \tilde{y}_i\|_{WRMS}` .. c:enumerator:: ARK_ACCUMERROR_AVG - Computes :math:`\frac{\text{reltol}}{|\mathcal{S}|} \sum_{i \in \mathcal{S}} \|y_i - \tilde{y}_i\|_{WRMS}`. + Computes :math:`\frac{\text{reltol}}{\Delta t_{\mathcal{S}}} \sum\limits_{i \in \mathcal{S}} h_i \|y_i - \tilde{y}_i\|_{WRMS}`, + where :math:`h_i` is the step size used when computing :math:`y_i`, and + :math:`\Delta t_{\mathcal{S}}` denotes the elapsed time over which + :math:`\mathcal{S}` is taken. .. c:function:: int ARKodeSetAccumulatedErrorType(void* arkode_mem, ARKAccumError accum_type) Sets the strategy to use for accumulating a temporal error estimate - over multiple time steps. + over multiple time steps. By default, ARKODE will not accumulate any + local error estimates (i.e., the default *accum_type* is ``ARK_ACCUMERROR_NONE``). + + A non-default error accumulation strategy can be disabled by calling + :c:func:`ARKodeSetAccumulatedErrorType` with the argument ``ARK_ACCUMERROR_NONE``. + :param arkode_mem: pointer to the ARKODE memory block. :param accum_type: accumulation strategy. @@ -1807,14 +1817,6 @@ user-specified relative solution tolerance. :retval ARK_STEPPER_UNSUPPORTED: temporal error estimation is not supported by the current time-stepping module. - .. note:: - - By default, ARKODE will not accumulate any local error estimates (i.e., - the default *accum_type* is ``ARK_ACCUMERROR_NONE``). - - A non-default error accumulation strategy can be disabled by calling - :c:func:`ARKodeSetAccumulatedErrorType` with the argument ``ARK_ACCUMERROR_NONE``. - .. versionadded:: x.y.z diff --git a/doc/shared/sunadaptcontroller/SUNAdaptController_Description.rst b/doc/shared/sunadaptcontroller/SUNAdaptController_Description.rst index fa7e048017..3407edc9f4 100644 --- a/doc/shared/sunadaptcontroller/SUNAdaptController_Description.rst +++ b/doc/shared/sunadaptcontroller/SUNAdaptController_Description.rst @@ -86,7 +86,7 @@ The virtual table structure is defined as .. c:member:: SUNErrCode (*estimatesteptol)(SUNAdaptController C, sunrealtype H, sunrealtype tolfac, int P, sunrealtype DSM, sunrealtype dsm, sunrealtype* Hnew, sunrealtype* tolfacnew) The function implementing :c:func:`SUNAdaptController_EstimateStepTol` - + .. versionadded:: x.y.z .. c:member:: SUNErrCode (*reset)(SUNAdaptController C) @@ -111,7 +111,7 @@ The virtual table structure is defined as .. c:member:: SUNErrCode (*updatemritol)(SUNAdaptController C, sunrealtype H, sunrealtype tolfac, sunrealtype DSM, sunrealtype dsm) - The function implementing :c:func:`SUNAdaptController_UpdateMRITol` + The function implementing :c:func:`SUNAdaptController_UpdateMRIHTol` .. versionadded:: x.y.z @@ -144,11 +144,11 @@ following set of SUNAdaptController types: Controls a single-rate step size. -.. c:enumerator:: SUN_ADAPTCONTROLLER_MRI_TOL +.. c:enumerator:: SUN_ADAPTCONTROLLER_MRI_H_TOL Controls both a slow time step and a tolerance factor to apply on the next-faster time scale within a multirate simulation that has an arbitrary number of time scales. - + .. versionadded:: x.y.z @@ -191,11 +191,6 @@ note these requirements below. Additionally, we note the behavior of the base SU :param C: the :c:type:`SUNAdaptController` object. :return: :c:type:`SUNErrCode` indicating success or failure. - Usage: - - .. code-block:: c - - retval = SUNAdaptController_DestroyEmpty(C); .. c:function:: SUNAdaptController_Type SUNAdaptController_GetType(SUNAdaptController C) @@ -205,11 +200,6 @@ note these requirements below. Additionally, we note the behavior of the base SU :param C: the :c:type:`SUNAdaptController` object. :return: :c:type:`SUNAdaptController_Type` type identifier. - Usage: - - .. code-block:: c - - SUNAdaptController_Type id = SUNAdaptController_GetType(C); .. c:function:: SUNErrCode SUNAdaptController_Destroy(SUNAdaptController C) @@ -223,11 +213,6 @@ note these requirements below. Additionally, we note the behavior of the base SU :param C: the :c:type:`SUNAdaptController` object. :return: :c:type:`SUNErrCode` indicating success or failure. - Usage: - - .. code-block:: c - - retval = SUNAdaptController_Destroy(C); .. c:function:: SUNErrCode SUNAdaptController_EstimateStep(SUNAdaptController C, sunrealtype h, int p, sunrealtype dsm, sunrealtype* hnew) @@ -242,20 +227,15 @@ note these requirements below. Additionally, we note the behavior of the base SU :param hnew: (output) the estimated step size. :return: :c:type:`SUNErrCode` indicating success or failure. - Usage: - - .. code-block:: c - - retval = SUNAdaptController_EstimateStep(C, hcur, p, dsm, &hnew); .. c:function:: SUNErrCode SUNAdaptController_EstimateStepTol(SUNAdaptController C, sunrealtype H, sunrealtype tolfac, int P, sunrealtype DSM, sunrealtype dsm, sunrealtype* Hnew, sunrealtype* tolfacnew) Estimates a slow step size and a fast tolerance multiplication factor for two adjacent time scales within a multirate application. - - This routine is required for controllers of type :c:enumerator`SUN_ADAPTCONTROLLER_MRI_TOL`. + + This routine is required for controllers of type :c:enumerator`SUN_ADAPTCONTROLLER_MRI_H_TOL`. If the current time scale has relative tolerance ``rtol``, then the - next-faster time scale will be called with relative tolerance ``tolfac * rtol``. + next-faster time scale will be called with relative tolerance ``tolfac * rtol``. If this is not provided by the implementation, the base class method will set ``*Hnew = H`` and ``*tolfacnew = tolfac`` and return. @@ -270,11 +250,7 @@ note these requirements below. Additionally, we note the behavior of the base SU :return: :c:type:`SUNErrCode` indicating success or failure. .. versionadded:: x.y.z - Usage: - - .. code-block:: c - retval = SUNAdaptController_EstimateStepTol(C, Hcur, tolfac, P, DSM, dsm, &Hnew, &tolfacnew); .. c:function:: SUNErrCode SUNAdaptController_Reset(SUNAdaptController C) @@ -283,7 +259,7 @@ note these requirements below. Additionally, we note the behavior of the base SU :param C: the :c:type:`SUNAdaptController` object. :return: :c:type:`SUNErrCode` indicating success or failure. - + .. versionadded:: x.y.z @@ -293,7 +269,7 @@ note these requirements below. Additionally, we note the behavior of the base SU :param C: the :c:type:`SUNAdaptController` object. :return: :c:type:`SUNErrCode` indicating success or failure. - + .. versionadded:: x.y.z @@ -332,9 +308,9 @@ note these requirements below. Additionally, we note the behavior of the base SU :return: :c:type:`SUNErrCode` indicating success or failure. -.. c:function:: SUNErrCode SUNAdaptController_UpdateMRITol(SUNAdaptController C, sunrealtype H, sunrealtype tolfac, sunrealtype DSM, sunrealtype dsm) +.. c:function:: SUNErrCode SUNAdaptController_UpdateMRIHTol(SUNAdaptController C, sunrealtype H, sunrealtype tolfac, sunrealtype DSM, sunrealtype dsm) - Notifies a controller of type :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_TOL` that a successful time step + Notifies a controller of type :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_H_TOL` that a successful time step was taken with slow stepsize ``H`` and fast relative tolerance factor ``tolfac``, and that the step had slow and fast local error factors ``DSM`` and ``dsm``, indicating that these can be saved for subsequent controller functions. This is typically relevant for controllers that @@ -346,7 +322,7 @@ note these requirements below. Additionally, we note the behavior of the base SU :param DSM: the successful slow temporal error estimate. :param dsm: the successful fast temporal error estimate. :return: :c:type:`SUNErrCode` indicating success or failure. - + .. versionadded:: x.y.z diff --git a/doc/shared/sunadaptcontroller/SUNAdaptController_MRIHTol.rst b/doc/shared/sunadaptcontroller/SUNAdaptController_MRIHTol.rst index db8136668b..1d7c8c9f8c 100644 --- a/doc/shared/sunadaptcontroller/SUNAdaptController_MRIHTol.rst +++ b/doc/shared/sunadaptcontroller/SUNAdaptController_MRIHTol.rst @@ -55,8 +55,8 @@ results from two types of error: :label: inner_solver_assumption where :math:`\text{RTOL}_n^F = \text{RTOL}^S \text{tolfac}_n^F`, - :math:`\text{RTOL}^S` is the relative tolerance that was supplied to the - current time scale solver, and where + the relative tolerance that was supplied to the current time scale + solver is :math:`\text{RTOL}^S`, and :math:`\kappa(t_n) = c(t_n) h_n^S \text{RTOL}^S` is independent of the relative tolerance factor, :math:`\text{tolfac}_n^F`. @@ -93,38 +93,39 @@ Implementation -------------- The SUNAdaptController_MRIHTol controller is implemented as a derived -SUNAdaptController class, and defines its *content* field as: +:c:type:`SUNAdaptController` class, and its *content* field is defined by +the :c:struct:`SUNAdaptControllerContent_MRIHTol_` structure: -.. code-block:: c +.. c:struct:: SUNAdaptControllerContent_MRIHTol_ - struct SUNAdaptControllerContent_MRIHTol_ - { - SUNAdaptController HControl; - SUNAdaptController TolControl; - sunrealtype inner_max_relch; - sunrealtype inner_min_tolfac; - sunrealtype inner_max_tolfac; - }; + The member data structure for an MRIHTol controller -These entries of the *content* field contain the following information: + .. c:member:: SUNAdaptController HControl -* ``HControl`` - single time-scale SUNAdaptController object to adapt - the current step size, :math:`h^S_n`. + A single time-scale controller to adapt the current step size, :math:`h^S_n`. -* ``TolControl`` - single time-scale SUNAdaptController object to adapt - the inner solver relative tolerance factor, :math:`\text{RTOL}^F_n`. + .. c:member:: SUNAdaptController TolControl -* ``inner_max_relch`` - the parameter :math:`relch_{\text{max}}` above. + A single time-scale controller to adapt the inner solver relative tolerance + factor, :math:`\text{reltol}^F_n`. -* ``inner_min_tolfac`` - the parameter :math:`\text{tolfac}_{min}` above. + .. c:member:: sunrealtype inner_max_relch -* ``inner_max_tolfac`` - the parameter :math:`\text{tolfac}_{max}` above. + The parameter :math:`relch_{\text{max}}` above. + + .. c:member:: sunrealtype inner_min_tolfac + + The parameter :math:`\text{tolfac}_{min}` above. + + .. c:member:: sunrealtype inner_max_tolfac + + The parameter :math:`\text{tolfac}_{max}` above. The header file to be included when using this module is ``sunadaptcontroller/sunadaptcontroller_mrihtol.h``. The SUNAdaptController_MRIHTol class provides implementations of all operations -relevant to a :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_TOL` controller listed in +relevant to a :c:enumerator:`SUN_ADAPTCONTROLLER_MRI_H_TOL` controller listed in :numref:`SUNAdaptController.Description.operations`. This class also provides the following additional user-callable routines: @@ -152,3 +153,60 @@ also provides the following additional user-callable routines: :param inner_min_tolfac: the parameter :math:`\text{tolfac}_{min}`. :param inner_max_tolfac: the parameter :math:`\text{tolfac}_{max}`. :returns: :c:type:`SUNErrCode` indicating success or failure. + + +Usage +----- + +Since this adaptivity controller is constructed using multiple single-rate adaptivity +controllers, there are a few steps required when setting this up in an application +(the steps below in *italics* correspond to the surrounding steps described in the +:ref:`MRIStep usage skeleton `. + +#. *Create an inner stepper object to solve the fast (inner) IVP* + +#. Configure the inner stepper to use temporal adaptivity. For example, when using + an ARKODE inner stepper and the :c:func:`ARKodeCreateMRIStepInnerStepper` + function, then either use its default adaptivity approach or supply a + single-rate SUNAdaptController object, e.g. + + .. code:: C + + void* inner_arkode_mem = ERKStepCreate(f_f, T0, y, sunctx); + MRIStepInnerStepper inner_stepper = nullptr; + retval = ARKodeCreateMRIStepInnerStepper(inner_arkode_mem, &inner_stepper); + SUNAdaptController fcontrol = SUNAdaptController_PID(sunctx); + retval = ARKodeSetAdaptController(inner_arkode_mem, fcontrol); + +#. If using an ARKODE inner stepper, then set the desired temporal error accumulation + estimation strategy via a call to :c:func:`ARKodeSetAccumulatedErrorType`, e.g., + + .. code:: C + + retval = ARKodeSetAccumulatedErrorType(inner_arkode_mem, ARK_ACCUMERROR_MAX); + +#. *Create an MRIStep object for the slow (outer) integration* + +#. Create single-rate controllers for both the slow step size and inner solver + tolerance, e.g., + + .. code:: C + + SUNAdaptController scontrol_H = SUNAdaptController_PI(sunctx); + SUNAdaptController scontrol_Tol = SUNAdaptController_I(sunctx); + +#. Create the multirate controller object, e.g., + + .. code:: C + + SUNAdaptController scontrol = SUNAdaptController_MRIHTol(scontrol_H, scontrol_Tol, sunctx); + +#. Attach the multirate controller object to MRIStep, e.g., + + .. code:: C + + retval = ARKodeSetAdaptController(arkode_mem, scontrol); + +An example showing the above steps is provided in +``examples/arkode/CXX_serial/ark_kpr_nestedmri.cpp``, where multirate controller objects +are used for both the slow and intermediate time scales in a 3-time-scale simulation. diff --git a/examples/arkode/CXX_serial/ark_kpr_nestedmri.cpp b/examples/arkode/CXX_serial/ark_kpr_nestedmri.cpp index 42568dba60..56d4e68a1f 100644 --- a/examples/arkode/CXX_serial/ark_kpr_nestedmri.cpp +++ b/examples/arkode/CXX_serial/ark_kpr_nestedmri.cpp @@ -934,11 +934,11 @@ int main(int argc, char* argv[]) // Main time-stepping loop: calls ARKodeEvolve to perform the // integration, then prints results. Stops when the final time // has been reached - sunrealtype t = T0; - sunrealtype t2 = T0; - sunrealtype dTout = (Tf - T0) / Nt; - sunrealtype tout = T0 + dTout; - sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype t = T0; + sunrealtype t2 = T0; + sunrealtype dTout = (Tf - T0) / Nt; + sunrealtype tout = T0 + dTout; + sunrealtype* ydata = N_VGetArrayPointer(y); sunrealtype* yrefdata = N_VGetArrayPointer(yref); sunrealtype u, v, w, uerr, verr, werr, uerrtot, verrtot, werrtot, errtot, accuracy; @@ -949,9 +949,9 @@ int main(int argc, char* argv[]) "---------------------------------------------------------------------" "-------\n"); printf(" %10.6" FSYM " %10.6" FSYM " %10.6" FSYM " %10.6" FSYM " %.2" ESYM - while (Tf - t > SUN_RCONST(1.0e-8)) + " %.2" ESYM " %.2" ESYM "\n", t, ydata[0], ydata[1], ydata[2], uerr, verr, werr); - while (Tf - t > SUN_RCONST(1.0e-8)) + while (Tf - t > SUN_RCONST(1.0e-8)) { // reset reference solver so that it begins with identical state retval = ARKodeReset(arkode_ref, t, y); @@ -985,15 +985,12 @@ int main(int argc, char* argv[]) verrtot += verr * verr; werrtot += werr * werr; errtot += uerr * uerr + verr * verr + werr * werr; - accuracy = - std::max(accuracy, - uerr / std::abs(opts.atol + opts.rtol * NV_Ith_S(yref, 0))); - accuracy = - std::max(accuracy, - verr / std::abs(opts.atol + opts.rtol * NV_Ith_S(yref, 1))); - accuracy = - std::max(accuracy, - werr / std::abs(opts.atol + opts.rtol * NV_Ith_S(yref, 2))); + accuracy = std::max(accuracy, + uerr / std::abs(opts.atol + opts.rtol * yrefdata[0])); + accuracy = std::max(accuracy, + verr / std::abs(opts.atol + opts.rtol * yrefdata[1])); + accuracy = std::max(accuracy, + werr / std::abs(opts.atol + opts.rtol * yrefdata[2])); // Periodically output current results to screen if (t >= tout) @@ -1021,8 +1018,10 @@ int main(int argc, char* argv[]) check_flag(retval, "ARKodeGetNumStepAttempts"); retval = ARKodeGetNumErrTestFails(arkode_mem, &netfs); check_flag(retval, "ARKodeGetNumErrTestFails"); - retval = MRIStepGetNumRhsEvals(arkode_mem, &nfse, &nfsi); - check_flag(retval, "MRIStepGetNumRhsEvals"); + retval = ARKodeGetNumRhsEvals(arkode_mem, 0, &nfse); + check_flag(retval, "ARKodeGetNumRhsEvals"); + retval = ARKodeGetNumRhsEvals(arkode_mem, 1, &nfsi); + check_flag(retval, "ARKodeGetNumRhsEvals"); // Get some intermediate integrator statistics long int nstm, nattm, netfm, nfme, nfmi; @@ -1032,8 +1031,10 @@ int main(int argc, char* argv[]) check_flag(retval, "ARKodeGetNumStepAttempts"); retval = ARKodeGetNumErrTestFails(mid_arkode_mem, &netfm); check_flag(retval, "ARKodeGetNumErrTestFails"); - retval = MRIStepGetNumRhsEvals(mid_arkode_mem, &nfme, &nfmi); - check_flag(retval, "MRIStepGetNumRhsEvals"); + retval = ARKodeGetNumRhsEvals(mid_arkode_mem, 0, &nfme); + check_flag(retval, "ARKodeGetNumRhsEvals"); + retval = ARKodeGetNumRhsEvals(mid_arkode_mem, 1, &nfmi); + check_flag(retval, "ARKodeGetNumRhsEvals"); // Get some fast integrator statistics long int nstf, nattf, netff, nff; @@ -1043,14 +1044,14 @@ int main(int argc, char* argv[]) check_flag(retval, "ARKodeGetNumStepAttempts"); retval = ARKodeGetNumErrTestFails(inner_arkode_mem, &netff); check_flag(retval, "ARKodeGetNumErrTestFails"); - retval = ERKStepGetNumRhsEvals(inner_arkode_mem, &nff); - check_flag(retval, "ERKStepGetNumRhsEvals"); + retval = ARKodeGetNumRhsEvals(inner_arkode_mem, 0, &nff); + check_flag(retval, "ARKodeGetNumRhsEvals"); // Print some final statistics - uerrtot = std::sqrt(uerrtot / nsts); - verrtot = std::sqrt(verrtot / nsts); - werrtot = std::sqrt(werrtot / nsts); - errtot = std::sqrt(errtot / nsts / 3); + uerrtot = std::sqrt(uerrtot / (sunrealtype)nsts); + verrtot = std::sqrt(verrtot / (sunrealtype)nsts); + werrtot = std::sqrt(werrtot / (sunrealtype)nsts); + errtot = std::sqrt(errtot / SUN_RCONST(3.0) / (sunrealtype)nsts); std::cout << "\nFinal Solver Statistics:\n"; std::cout << " Slow steps = " << nsts << " (attempts = " << natts << ", fails = " << netfs << ")\n"; @@ -1074,10 +1075,12 @@ int main(int argc, char* argv[]) retval = ARKodeGetNumJacEvals(arkode_mem, &njes); check_flag(retval, "ARKodeGetNumJacEvals"); std::cout << " Slow Newton iters = " << nnis << std::endl; - std::cout << " Slow Newton iters/attempt = " << nnis/natts << std::endl; + std::cout << " Slow Newton iters/attempt = " + << (sunrealtype)nnis / (sunrealtype)natts << std::endl; std::cout << " Slow Newton conv fails = " << nncs << std::endl; std::cout << " Slow Jacobian evals = " << njes << std::endl; - std::cout << " Slow Jacobian evals/Newton = " << njes/nnis << std::endl; + std::cout << " Slow Jacobian evals/Newton = " + << (sunrealtype)njes / (sunrealtype)nnis << std::endl; } // Get/print intermediate integrator implicit solver statistics @@ -1089,10 +1092,12 @@ int main(int argc, char* argv[]) retval = ARKodeGetNumJacEvals(mid_arkode_mem, &njem); check_flag(retval, "ARKodeGetNumJacEvals"); std::cout << " Intermediate Newton iters = " << nnim << std::endl; - std::cout << " Intermediate Newton iters/attempt = " << nnim/nattm << std::endl; + std::cout << " Intermediate Newton iters/attempt = " + << (sunrealtype)nnim / (sunrealtype)nattm << std::endl; std::cout << " Intermediate Newton conv fails = " << nncm << std::endl; std::cout << " Intermediate Jacobian evals = " << njem << std::endl; - std::cout << " Intermediate Jacobian evals/Newton = " << njem/nnim << std::endl; + std::cout << " Intermediate Jacobian evals/Newton = " + << (sunrealtype)njem / (sunrealtype)nnim << std::endl; } // Clean up and return @@ -1144,14 +1149,12 @@ static int fn(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) // [ G e e ] [(u^2-p-2)/(2u)] + [ pdot(t)/(2u) ] // [ e al be ] [(v^2-q-2)/(2v)] [ qdot(t)/(2v) ] // [ e -be al ] [(w^2-r-2)/(2w)] [ rdot(t)/(2w) ] - tmp1 = (-TWO + u * u - p(t, *opts)) / (TWO * u); - tmp2 = (-TWO + v * v - q(t, *opts)) / (TWO * v); - tmp3 = (-TWO + w * w - r(t, *opts)) / (TWO * w); + tmp1 = (-TWO + u * u - p(t, *opts)) / (TWO * u); + tmp2 = (-TWO + v * v - q(t, *opts)) / (TWO * v); + tmp3 = (-TWO + w * w - r(t, *opts)) / (TWO * w); ydotdata[0] = G * tmp1 + e * tmp2 + e * tmp3 + pdot(t, *opts) / (TWO * u); - ydotdata[1] = e * tmp1 + al * tmp2 + be * tmp3 + - qdot(t, *opts) / (TWO * v); - ydotdata[2] = e * tmp1 - be * tmp2 + al * tmp3 + - rdot(t, *opts) / (TWO * w); + ydotdata[1] = e * tmp1 + al * tmp2 + be * tmp3 + qdot(t, *opts) / (TWO * v); + ydotdata[2] = e * tmp1 - be * tmp2 + al * tmp3 + rdot(t, *opts) / (TWO * w); // Return with success return 0; @@ -1175,13 +1178,12 @@ static int ff(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) // [ 0 0 0 ] [(u^2-p-2)/(2u)] + [ 0 ] // [ 0 0 0 ] [(v^2-q-2)/(2v)] [ 0 ] // [ e -be al ] [(w^2-r-2)/(2w)] [ rdot(t)/(2w) ] - tmp1 = (-TWO + u * u - p(t, *opts)) / (TWO * u); - tmp2 = (-TWO + v * v - q(t, *opts)) / (TWO * v); - tmp3 = (-TWO + w * w - r(t, *opts)) / (TWO * w); + tmp1 = (-TWO + u * u - p(t, *opts)) / (TWO * u); + tmp2 = (-TWO + v * v - q(t, *opts)) / (TWO * v); + tmp3 = (-TWO + w * w - r(t, *opts)) / (TWO * w); ydotdata[0] = ZERO; ydotdata[1] = ZERO; - ydotdata[2] = e * tmp1 - be * tmp2 + al * tmp3 + - rdot(t, *opts) / (TWO * w); + ydotdata[2] = e * tmp1 - be * tmp2 + al * tmp3 + rdot(t, *opts) / (TWO * w); // Return with success return 0; @@ -1205,12 +1207,11 @@ static int fm(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) // [ 0 0 0 ] [(u^2-p-2)/(2u)] + [ 0 ] // [ e al be ] [(v^2-q-2)/(2v)] [ qdot(t)/(2v) ] // [ 0 0 0 ] [(w^2-r-2)/(2w)] [ 0 ] - tmp1 = (-TWO + u * u - p(t, *opts)) / (TWO * u); - tmp2 = (-TWO + v * v - q(t, *opts)) / (TWO * v); - tmp3 = (-TWO + w * w - r(t, *opts)) / (TWO * w); + tmp1 = (-TWO + u * u - p(t, *opts)) / (TWO * u); + tmp2 = (-TWO + v * v - q(t, *opts)) / (TWO * v); + tmp3 = (-TWO + w * w - r(t, *opts)) / (TWO * w); ydotdata[0] = ZERO; - ydotdata[1] = e * tmp1 + al * tmp2 + be * tmp3 + - qdot(t, *opts) / (TWO * v); + ydotdata[1] = e * tmp1 + al * tmp2 + be * tmp3 + qdot(t, *opts) / (TWO * v); ydotdata[2] = ZERO; return 0; @@ -1281,9 +1282,9 @@ static int fs(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) // [ G e e ] [(u^2-p-2)/(2u)] + [ pdot(t)/(2u) ] // [ 0 0 0 ] [(v^2-q-2)/(2v)] [ 0 ] // [ 0 0 0 ] [(w^2-r-2)/(2w)] [ 0 ] - tmp1 = (-TWO + u * u - p(t, *opts)) / (TWO * u); - tmp2 = (-TWO + v * v - q(t, *opts)) / (TWO * v); - tmp3 = (-TWO + w * w - r(t, *opts)) / (TWO * w); + tmp1 = (-TWO + u * u - p(t, *opts)) / (TWO * u); + tmp2 = (-TWO + v * v - q(t, *opts)) / (TWO * v); + tmp3 = (-TWO + w * w - r(t, *opts)) / (TWO * w); ydotdata[0] = G * tmp1 + e * tmp2 + e * tmp3 + pdot(t, *opts) / (TWO * u); ydotdata[1] = ZERO; ydotdata[2] = ZERO; @@ -1294,7 +1295,7 @@ static int fs(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) // fse routine to compute the explicit slow portion of the ODE RHS. static int fse(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { - Options* opts. = static_cast(user_data); + Options* opts = static_cast(user_data); sunrealtype* ydata = N_VGetArrayPointer(y); sunrealtype* ydotdata = N_VGetArrayPointer(ydot); const sunrealtype u = ydata[0]; @@ -1342,11 +1343,11 @@ static int fsi(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int Jm(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { - Options* opts = static_cast(user_data); - sunrealtype* ydata = N_VGetArrayPointer(y); - const sunrealtype u = ydata[0]; - const sunrealtype v = ydata[1]; - const sunrealtype w = ydata[2]; + Options* opts = static_cast(user_data); + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; + const sunrealtype w = ydata[2]; sunrealtype t11, t22, t33; // fill in the Jacobian: @@ -1374,11 +1375,11 @@ static int Jm(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, static int Jmi(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { - Options* opts = static_cast(user_data); - sunrealtype* ydata = N_VGetArrayPointer(y); - const sunrealtype u = ydata[0]; - const sunrealtype v = ydata[1]; - const sunrealtype w = ydata[2]; + Options* opts = static_cast(user_data); + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; + const sunrealtype w = ydata[2]; sunrealtype t11, t22, t33; // fill in the Jacobian: @@ -1406,11 +1407,11 @@ static int Jmi(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, static int Js(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { - Options* opts = static_cast(user_data); - sunrealtype* ydata = N_VGetArrayPointer(y); - const sunrealtype u = ydata[0]; - const sunrealtype v = ydata[1]; - const sunrealtype w = ydata[2]; + Options* opts = static_cast(user_data); + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; + const sunrealtype w = ydata[2]; sunrealtype t11, t22, t33; // fill in the Jacobian: @@ -1438,11 +1439,11 @@ static int Js(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, static int Jsi(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { - Options* opts = static_cast(user_data); - sunrealtype* ydata = N_VGetArrayPointer(y); - const sunrealtype u = ydata[0]; - const sunrealtype v = ydata[1]; - const sunrealtype w = ydata[2]; + Options* opts = static_cast(user_data); + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; + const sunrealtype w = ydata[2]; sunrealtype t11, t22, t33; // fill in the Jacobian: diff --git a/examples/arkode/C_serial/CMakeLists.txt b/examples/arkode/C_serial/CMakeLists.txt index 74e5665ea9..4efb840356 100644 --- a/examples/arkode/C_serial/CMakeLists.txt +++ b/examples/arkode/C_serial/CMakeLists.txt @@ -61,7 +61,7 @@ set(ARKODE_examples "ark_kepler\;\;develop" "ark_kpr_mri\;0 1 0.005\;develop" "ark_kpr_mri\;1 0 0.01\;develop" - "ark_kpr_mri\;1 1 0.002\;develop" + "ark_kpr_mri\;1 1 0.002\;exclude-single" "ark_kpr_mri\;2 4 0.002\;develop" "ark_kpr_mri\;3 2 0.001\;develop" "ark_kpr_mri\;4 3 0.001\;develop" @@ -73,7 +73,7 @@ set(ARKODE_examples "ark_kpr_mri\;10 4 0.001 -100 100 0.5 1\;exclude-single" "ark_kpr_mri\;11 2 0.001\;develop" "ark_kpr_mri\;12 3 0.005\;develop" - "ark_kpr_mri\;13 4 0.01\;develop" + "ark_kpr_mri\;13 4 0.01\;exclude-single" "ark_KrylovDemo_prec\;\;exclude-single" "ark_KrylovDemo_prec\;1\;exclude-single" "ark_KrylovDemo_prec\;2\;exclude-single" diff --git a/examples/arkode/C_serial/ark_kpr_mri.c b/examples/arkode/C_serial/ark_kpr_mri.c index d9e96e6a1f..6edfb14cee 100644 --- a/examples/arkode/C_serial/ark_kpr_mri.c +++ b/examples/arkode/C_serial/ark_kpr_mri.c @@ -190,8 +190,8 @@ int main(int argc, char* argv[]) printf(" ark_kpr_mri slow_type fast_type h G w e deduce_rhs"); return (-1); } - slow_type = atoi(argv[1]); - fast_type = atoi(argv[2]); + slow_type = atoi(argv[1]); + fast_type = atoi(argv[2]); if (argc > 3) { hs = SUNStrToReal(argv[3]); } if (argc > 4) { G = SUNStrToReal(argv[4]); } if (argc > 5) { w = SUNStrToReal(argv[5]); } @@ -212,7 +212,7 @@ int main(int argc, char* argv[]) } if ((fast_type < 0) || (fast_type > 5)) { - printf("ERROR: fast_type be an integer in [0,5] \n"); + printf("ERROR: fast_type be an integer in [0,5] \n"); return (-1); } if ((slow_type == 0) && (fast_type == 0)) @@ -220,7 +220,7 @@ int main(int argc, char* argv[]) printf("ERROR: at least one of slow_type and fast_type must be nonzero\n"); return (-1); } - if (((slow_type == 9) || (slow_type == 10)) && (fast_type == 0)) + if ((slow_type >= 9) && (fast_type == 0)) { printf("ERROR: example not configured for ImEx slow solver with no fast " "solver\n"); @@ -432,8 +432,7 @@ int main(int argc, char* argv[]) /* Set Butcher table for fast integrator */ switch (fast_type) { - case (0): - break; + case (0): break; case (1): B = ARKodeButcherTable_Alloc(3, SUNFALSE); if (check_retval((void*)B, "ARKodeButcherTable_Alloc", 0)) { return 1; } diff --git a/include/arkode/arkode_mristep.h b/include/arkode/arkode_mristep.h index cd972e94ff..39be4a4407 100644 --- a/include/arkode/arkode_mristep.h +++ b/include/arkode/arkode_mristep.h @@ -22,6 +22,7 @@ #include #include #include +#include #ifdef __cplusplus /* wrapper to enable C++ usage */ extern "C" { @@ -118,7 +119,6 @@ typedef int (*MRIStepInnerGetAccumulatedError)(MRIStepInnerStepper stepper, typedef int (*MRIStepInnerResetAccumulatedError)(MRIStepInnerStepper stepper); - typedef int (*MRIStepInnerSetRTol)(MRIStepInnerStepper stepper, sunrealtype rtol); /*--------------------------------------------------------------- @@ -127,13 +127,13 @@ typedef int (*MRIStepInnerSetRTol)(MRIStepInnerStepper stepper, sunrealtype rtol struct MRIStepCouplingMem { MRISTEP_METHOD_TYPE type; /* flag to encode the MRI method type */ - int nmat; /* number of MRI coupling matrices */ - int stages; /* size of coupling matrices ((stages+1) * stages) */ - int q; /* method order of accuracy */ - int p; /* embedding order of accuracy */ - sunrealtype* c; /* stage abscissae */ - sunrealtype*** W; /* explicit coupling matrices [nmat][stages+1][stages] */ - sunrealtype*** G; /* implicit coupling matrices [nmat][stages+1][stages] */ + int nmat; /* number of MRI coupling matrices */ + int stages; /* size of coupling matrices ((stages+1) * stages) */ + int q; /* method order of accuracy */ + int p; /* embedding order of accuracy */ + sunrealtype* c; /* stage abscissae */ + sunrealtype*** W; /* explicit coupling matrices [nmat][stages+1][stages] */ + sunrealtype*** G; /* implicit coupling matrices [nmat][stages+1][stages] */ int ngroup; /* number of stage groups (MERK-specific) */ int** group; /* stages to integrate together (MERK-specific) */ diff --git a/include/sunadaptcontroller/sunadaptcontroller_mrihtol.h b/include/sunadaptcontroller/sunadaptcontroller_mrihtol.h index e53f17327f..9f100353f1 100644 --- a/include/sunadaptcontroller/sunadaptcontroller_mrihtol.h +++ b/include/sunadaptcontroller/sunadaptcontroller_mrihtol.h @@ -44,9 +44,9 @@ typedef struct SUNAdaptControllerContent_MRIHTol_* SUNAdaptControllerContent_MRI * ------------------ */ SUNDIALS_EXPORT -SUNAdaptController SUNAdaptController_MRIHTol(SUNAdaptController HControl, - SUNAdaptController TolControl, - SUNContext sunctx); +SUNAdaptController SUNAdaptController_MRIHTol(SUNAdaptController HControl, + SUNAdaptController TolControl, + SUNContext sunctx); SUNDIALS_EXPORT SUNErrCode SUNAdaptController_SetParams_MRIHTol(SUNAdaptController C, sunrealtype inner_max_relch, @@ -74,9 +74,9 @@ SUNDIALS_EXPORT int SUNAdaptController_SetErrorBias_MRIHTol(SUNAdaptController C, sunrealtype bias); SUNDIALS_EXPORT -int SUNAdaptController_UpdateMRITol_MRIHTol(SUNAdaptController C, sunrealtype H, - sunrealtype tolfac, sunrealtype DSM, - sunrealtype dsm); +int SUNAdaptController_UpdateMRIHTol_MRIHTol(SUNAdaptController C, + sunrealtype H, sunrealtype tolfac, + sunrealtype DSM, sunrealtype dsm); SUNDIALS_EXPORT int SUNAdaptController_Space_MRIHTol(SUNAdaptController C, long int* lenrw, long int* leniw); diff --git a/include/sundials/sundials_adaptcontroller.h b/include/sundials/sundials_adaptcontroller.h index 4809ea64eb..4d5c9bf9ca 100644 --- a/include/sundials/sundials_adaptcontroller.h +++ b/include/sundials/sundials_adaptcontroller.h @@ -41,7 +41,7 @@ typedef enum { SUN_ADAPTCONTROLLER_NONE, SUN_ADAPTCONTROLLER_H, - SUN_ADAPTCONTROLLER_MRI_TOL + SUN_ADAPTCONTROLLER_MRI_H_TOL } SUNAdaptController_Type; /* ----------------------------------------------------------------- @@ -64,7 +64,7 @@ struct _generic_SUNAdaptController_Ops SUNErrCode (*estimatestep)(SUNAdaptController C, sunrealtype h, int p, sunrealtype dsm, sunrealtype* hnew); - /* REQUIRED for controllers of SUN_ADAPTCONTROLLER_MRI_TOL type. */ + /* REQUIRED for controllers of SUN_ADAPTCONTROLLER_MRI_H_TOL type. */ SUNErrCode (*estimatesteptol)(SUNAdaptController C, sunrealtype H, sunrealtype tolfac, int P, sunrealtype DSM, sunrealtype dsm, sunrealtype* Hnew, @@ -77,9 +77,9 @@ struct _generic_SUNAdaptController_Ops SUNErrCode (*write)(SUNAdaptController C, FILE* fptr); SUNErrCode (*seterrorbias)(SUNAdaptController C, sunrealtype bias); SUNErrCode (*updateh)(SUNAdaptController C, sunrealtype h, sunrealtype dsm); - SUNErrCode (*updatemritol)(SUNAdaptController C, sunrealtype H, - sunrealtype tolfac, sunrealtype DSM, - sunrealtype dsm); + SUNErrCode (*updatemrihtol)(SUNAdaptController C, sunrealtype H, + sunrealtype tolfac, sunrealtype DSM, + sunrealtype dsm); SUNErrCode (*space)(SUNAdaptController C, long int* lenrw, long int* leniw); }; @@ -174,9 +174,9 @@ SUNErrCode SUNAdaptController_UpdateH(SUNAdaptController C, sunrealtype h, DSM and dsm, indicating that the step size, tolerance factor, or local error factors can be saved for subsequent controller functions. */ SUNDIALS_EXPORT -SUNErrCode SUNAdaptController_UpdateMRITol(SUNAdaptController C, sunrealtype H, - sunrealtype tolfac, sunrealtype DSM, - sunrealtype dsm); +SUNErrCode SUNAdaptController_UpdateMRIHTol(SUNAdaptController C, sunrealtype H, + sunrealtype tolfac, sunrealtype DSM, + sunrealtype dsm); /* Function to return the memory requirements of the controller object. */ SUNDIALS_EXPORT diff --git a/src/arkode/arkode.c b/src/arkode/arkode.c index fcdaca944e..c55c6bde33 100644 --- a/src/arkode/arkode.c +++ b/src/arkode/arkode.c @@ -815,7 +815,7 @@ int ARKodeEvolve(void* arkode_mem, sunrealtype tout, N_Vector yout, /* Check for too much accuracy requested */ nrm = N_VWrmsNorm(ark_mem->yn, ark_mem->ewt); ark_mem->tolsf = ark_mem->uround * nrm; - if (ark_mem->tolsf > ONE) + if (ark_mem->tolsf > ONE && !ark_mem->fixedstep) { arkProcessError(ark_mem, ARK_TOO_MUCH_ACC, __LINE__, __func__, __FILE__, MSG_ARK_TOO_MUCH_ACC, ark_mem->tcur); @@ -1373,7 +1373,7 @@ int ARKodeCreateMRIStepInnerStepper(void* inner_arkode_mem, ark_mem = (ARKodeMem)inner_arkode_mem; /* return with an error if the ARKODE solver does not support forcing */ - if (!ark_mem->step_supports_forcing) + if (ark_mem->step_setforcing == NULL) { arkProcessError(ark_mem, ARK_STEPPER_UNSUPPORTED, __LINE__, __func__, __FILE__, "time-stepping module does not support forcing"); @@ -1851,6 +1851,14 @@ int arkInitialSetup(ARKodeMem ark_mem, sunrealtype tout) sunrealtype tout_hin, rh, htmp; sunbooleantype conOK; + /* Check that user has supplied an initial step size if fixedstep mode is on */ + if ((ark_mem->fixedstep) && (ark_mem->hin == ZERO)) + { + arkProcessError(ark_mem, ARK_ILL_INPUT, __LINE__, __func__, __FILE__, + "Fixed step mode enabled, but no step size set"); + return (ARK_ILL_INPUT); + } + /* If using a built-in routine for error/residual weights with abstol==0, ensure that N_VMin is available */ if ((!ark_mem->user_efun) && (ark_mem->atolmin0) && (!ark_mem->yn->ops->nvmin)) @@ -1867,6 +1875,30 @@ int arkInitialSetup(ARKodeMem ark_mem, sunrealtype tout) return (ARK_ILL_INPUT); } + /* Test input tstop for legality (correct direction of integration) */ + if (ark_mem->tstopset) + { + htmp = (ark_mem->h == ZERO) ? tout - ark_mem->tcur : ark_mem->h; + if ((ark_mem->tstop - ark_mem->tcur) * htmp <= ZERO) + { + arkProcessError(ark_mem, ARK_ILL_INPUT, __LINE__, __func__, __FILE__, + MSG_ARK_BAD_TSTOP, ark_mem->tstop, ark_mem->tcur); + return (ARK_ILL_INPUT); + } + } + + /* Check to see if y0 satisfies constraints */ + if (ark_mem->constraintsSet) + { + conOK = N_VConstrMask(ark_mem->constraints, ark_mem->yn, ark_mem->tempv1); + if (!conOK) + { + arkProcessError(ark_mem, ARK_ILL_INPUT, __LINE__, __func__, __FILE__, + MSG_ARK_Y0_FAIL_CONSTR); + return (ARK_ILL_INPUT); + } + } + /* Load initial error weights */ retval = ark_mem->efun(ark_mem->yn, ark_mem->ewt, ark_mem->e_data); if (retval != 0) @@ -1923,38 +1955,6 @@ int arkInitialSetup(ARKodeMem ark_mem, sunrealtype tout) } } - /* Check that user has supplied an initial step size if fixedstep mode is on */ - if ((ark_mem->fixedstep) && (ark_mem->hin == ZERO)) - { - arkProcessError(ark_mem, ARK_ILL_INPUT, __LINE__, __func__, __FILE__, - "Fixed step mode enabled, but no step size set"); - return (ARK_ILL_INPUT); - } - - /* Test input tstop for legality (correct direction of integration) */ - if (ark_mem->tstopset) - { - htmp = (ark_mem->h == ZERO) ? tout - ark_mem->tcur : ark_mem->h; - if ((ark_mem->tstop - ark_mem->tcur) * htmp <= ZERO) - { - arkProcessError(ark_mem, ARK_ILL_INPUT, __LINE__, __func__, __FILE__, - MSG_ARK_BAD_TSTOP, ark_mem->tstop, ark_mem->tcur); - return (ARK_ILL_INPUT); - } - } - - /* Check to see if y0 satisfies constraints */ - if (ark_mem->constraintsSet) - { - conOK = N_VConstrMask(ark_mem->constraints, ark_mem->yn, ark_mem->tempv1); - if (!conOK) - { - arkProcessError(ark_mem, ARK_ILL_INPUT, __LINE__, __func__, __FILE__, - MSG_ARK_Y0_FAIL_CONSTR); - return (ARK_ILL_INPUT); - } - } - /* Create default Hermite interpolation module (if needed) */ if (ark_mem->interp_type != ARK_INTERP_NONE && !(ark_mem->interp)) { @@ -2594,7 +2594,11 @@ int arkCompleteStep(ARKodeMem ark_mem, sunrealtype dsm) { ark_mem->AccumError = SUNMAX(dsm, ark_mem->AccumError); } - else { ark_mem->AccumError += dsm; } + else if (ark_mem->AccumErrorType == ARK_ACCUMERROR_SUM) + { + ark_mem->AccumError += dsm; + } + else /* ARK_ACCUMERROR_AVG */ { ark_mem->AccumError += (dsm * ark_mem->h); } } /* apply user-supplied step postprocessing function (if supplied) */ diff --git a/src/arkode/arkode_arkstep.c b/src/arkode/arkode_arkstep.c index a2d25adebe..dfd941c6e9 100644 --- a/src/arkode/arkode_arkstep.c +++ b/src/arkode/arkode_arkstep.c @@ -142,7 +142,6 @@ void* ARKStepCreate(ARKRhsFn fe, ARKRhsFn fi, sunrealtype t0, N_Vector y0, ark_mem->step_setforcing = arkStep_SetInnerForcing; ark_mem->step_supports_adaptive = SUNTRUE; ark_mem->step_supports_implicit = SUNTRUE; - ark_mem->step_supports_forcing = SUNTRUE; ark_mem->step_supports_massmatrix = SUNTRUE; ark_mem->step_supports_relaxation = SUNTRUE; ark_mem->step_mem = (void*)step_mem; @@ -1807,20 +1806,20 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::arkStep_TakeStep_Z", "explicit stage", - "z[%i] =", 0); + "z_%i(:) =", 0); N_VPrintFile(ark_mem->ycur, ARK_LOGGER->debug_fp); if (step_mem->implicit) { SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::arkStep_TakeStep_Z", "implicit RHS", - "Fi[%i] =", 0); + "Fi_%i(:) =", 0); N_VPrintFile(step_mem->Fi[0], ARK_LOGGER->debug_fp); } if (step_mem->explicit) { SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::arkStep_TakeStep_Z", "explicit RHS", - "Fe[%i] =", 0); + "Fe_%i(:) =", 0); N_VPrintFile(step_mem->Fe[0], ARK_LOGGER->debug_fp); } #endif @@ -1888,7 +1887,8 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::arkStep_TakeStep_Z", "predictor", "zpred =", ""); + "ARKODE::arkStep_TakeStep_Z", "predictor", + "zpred(:) =", ""); N_VPrintFile(step_mem->zpred, ARK_LOGGER->debug_fp); #endif @@ -1898,7 +1898,8 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::arkStep_TakeStep_Z", "rhs data", "sdata =", ""); + "ARKODE::arkStep_TakeStep_Z", "rhs data", + "sdata(:) =", ""); N_VPrintFile(step_mem->sdata, ARK_LOGGER->debug_fp); #endif @@ -1913,7 +1914,7 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::arkStep_TakeStep_Z", "implicit stage", - "z[%i] =", is); + "z_%i(:) =", is); N_VPrintFile(ark_mem->ycur, ARK_LOGGER->debug_fp); #endif @@ -1937,7 +1938,7 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::arkStep_TakeStep_Z", "explicit stage", - "z[%i] =", is); + "z_%i(:) =", is); N_VPrintFile(ark_mem->ycur, ARK_LOGGER->debug_fp); #endif } @@ -1982,7 +1983,7 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::arkStep_TakeStep_Z", "implicit RHS", - "Fi[%i] =", is); + "Fi_%i(:) =", is); N_VPrintFile(step_mem->Fi[is], ARK_LOGGER->debug_fp); #endif @@ -2000,7 +2001,7 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::arkStep_TakeStep_Z", "explicit RHS", - "Fe[%i] =", is); + "Fe_%i(:) =", is); N_VPrintFile(step_mem->Fe[is], ARK_LOGGER->debug_fp); #endif @@ -2019,7 +2020,7 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::arkStep_TakeStep_Z", "M^{-1} implicit RHS", - "Fi[%i] =", is); + "Fi_%i(:) =", is); N_VPrintFile(step_mem->Fi[is], ARK_LOGGER->debug_fp); #endif if (*nflagPtr != ARK_SUCCESS) { return (TRY_AGAIN); } @@ -2031,7 +2032,7 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::arkStep_TakeStep_Z", "M^{-1} explicit RHS", - "Fe[%i] =", is); + "Fe_%i(:) =", is); N_VPrintFile(step_mem->Fe[is], ARK_LOGGER->debug_fp); #endif if (*nflagPtr != ARK_SUCCESS) { return (TRY_AGAIN); } @@ -2053,16 +2054,15 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::arkStep_TakeStep_Z", - "updated solution", "ycur =", ""); + "updated solution", "ycur(:) =", ""); N_VPrintFile(ark_mem->ycur, ARK_LOGGER->debug_fp); #endif #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_INFO SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_INFO, "ARKODE::arkStep_TakeStep_Z", "end-step", - "step = %li, t = %" RSYM ", h = %" RSYM ", dsm = %" RSYM - ", nflag = %d", - ark_mem->nst, ark_mem->tn, ark_mem->h, *dsmPtr, *nflagPtr); + "step = %li, h = %" RSYM ", dsm = %" RSYM ", nflag = %d", + ark_mem->nst, ark_mem->h, *dsmPtr, *nflagPtr); #endif return (ARK_SUCCESS); diff --git a/src/arkode/arkode_arkstep_io.c b/src/arkode/arkode_arkstep_io.c index 3ea4a5ce47..7c07248265 100644 --- a/src/arkode/arkode_arkstep_io.c +++ b/src/arkode/arkode_arkstep_io.c @@ -757,7 +757,7 @@ int arkStep_SetDefaults(ARKodeMem ark_mem) if (ark_mem->hadapt_mem->hcontroller == NULL) { arkProcessError(ark_mem, ARK_MEM_FAIL, __LINE__, __func__, __FILE__, - "SUNAdaptControllerPID allocation failure"); + "SUNAdaptController_PID allocation failure"); return (ARK_MEM_FAIL); } ark_mem->hadapt_mem->owncontroller = SUNTRUE; diff --git a/src/arkode/arkode_arkstep_nls.c b/src/arkode/arkode_arkstep_nls.c index 5f1811b156..939a82ad02 100644 --- a/src/arkode/arkode_arkstep_nls.c +++ b/src/arkode/arkode_arkstep_nls.c @@ -425,7 +425,7 @@ int arkStep_Nls(ARKodeMem ark_mem, int nflag) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::arkStep_Nls", - "correction", "zcor =", ""); + "correction", "zcor(:) =", ""); N_VPrintFile(step_mem->zcor, ARK_LOGGER->debug_fp); #endif diff --git a/src/arkode/arkode_erkstep.c b/src/arkode/arkode_erkstep.c index 873af9b888..3030ddcad2 100644 --- a/src/arkode/arkode_erkstep.c +++ b/src/arkode/arkode_erkstep.c @@ -107,7 +107,6 @@ void* ERKStepCreate(ARKRhsFn f, sunrealtype t0, N_Vector y0, SUNContext sunctx) ark_mem->step_setforcing = erkStep_SetInnerForcing; ark_mem->step_supports_adaptive = SUNTRUE; ark_mem->step_supports_relaxation = SUNTRUE; - ark_mem->step_supports_forcing = SUNTRUE; ark_mem->step_mem = (void*)step_mem; /* Set default values for optional inputs */ @@ -691,10 +690,10 @@ int erkStep_TakeStep(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::erkStep_TakeStep", - "stage", "z[0] =", ""); + "stage", "z_0(:) =", ""); N_VPrintFile(ark_mem->ycur, ARK_LOGGER->debug_fp); SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::erkStep_TakeStep", - "stage RHS", "F[0] =", ""); + "stage RHS", "F_0(:) =", ""); N_VPrintFile(step_mem->F[0], ARK_LOGGER->debug_fp); #endif @@ -770,7 +769,7 @@ int erkStep_TakeStep(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::erkStep_TakeStep", "stage RHS", "F[%i] =", is); + "ARKODE::erkStep_TakeStep", "stage RHS", "F_%i(:) =", is); N_VPrintFile(step_mem->F[is], ARK_LOGGER->debug_fp); #endif @@ -782,15 +781,14 @@ int erkStep_TakeStep(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::erkStep_TakeStep", - "updated solution", "ycur =", ""); + "updated solution", "ycur(:) =", ""); N_VPrintFile(ark_mem->ycur, ARK_LOGGER->debug_fp); #endif #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::erkStep_TakeStep", - "end-step", - "step = %li, t = %" RSYM ", h = %" RSYM ", dsm = %" RSYM, - ark_mem->nst, ark_mem->tn, ark_mem->h, *dsmPtr); + "error-test", "step = %li, h = %" RSYM ", dsm = %" RSYM, + ark_mem->nst, ark_mem->h, *dsmPtr); #endif return (ARK_SUCCESS); diff --git a/src/arkode/arkode_erkstep_io.c b/src/arkode/arkode_erkstep_io.c index d86a87317a..586d3bf580 100644 --- a/src/arkode/arkode_erkstep_io.c +++ b/src/arkode/arkode_erkstep_io.c @@ -326,7 +326,7 @@ int erkStep_SetDefaults(ARKodeMem ark_mem) if (ark_mem->hadapt_mem->hcontroller == NULL) { arkProcessError(ark_mem, ARK_MEM_FAIL, __LINE__, __func__, __FILE__, - "SUNAdaptControllerPI allocation failure"); + "SUNAdaptController_PI allocation failure"); return (ARK_MEM_FAIL); } ark_mem->hadapt_mem->owncontroller = SUNTRUE; diff --git a/src/arkode/arkode_impl.h b/src/arkode/arkode_impl.h index 1568ccb56f..029ca962a6 100644 --- a/src/arkode/arkode_impl.h +++ b/src/arkode/arkode_impl.h @@ -463,7 +463,6 @@ struct ARKodeMemRec ARKMassMultFn step_mmult; /* Time stepper module -- forcing */ - sunbooleantype step_supports_forcing; ARKTimestepSetForcingFn step_setforcing; /* N_Vector storage */ @@ -540,9 +539,9 @@ struct ARKodeMemRec sunrealtype terr; /* error in tn for compensated sums */ sunrealtype hold; /* last successful h value used */ sunrealtype tolsf; /* tolerance scale factor (suggestion to user) */ - ARKAccumError AccumErrorType; /* accumulated error estimation type */ - long int AccumErrorStep; /* time step of last accumulated error reset */ - sunrealtype AccumError; /* accumulated error estimate */ + ARKAccumError AccumErrorType; /* accumulated error estimation type */ + sunrealtype AccumErrorStart; /* time of last accumulated error reset */ + sunrealtype AccumError; /* accumulated error estimate */ sunbooleantype VabstolMallocDone; sunbooleantype VRabstolMallocDone; sunbooleantype MallocDone; diff --git a/src/arkode/arkode_io.c b/src/arkode/arkode_io.c index dfa7281a9e..d22ff863a4 100644 --- a/src/arkode/arkode_io.c +++ b/src/arkode/arkode_io.c @@ -1992,8 +1992,8 @@ int ARKodeResetAccumulatedError(void* arkode_mem) } /* Reset value and counter, and return */ - ark_mem->AccumErrorStep = ark_mem->nst; - ark_mem->AccumError = ZERO; + ark_mem->AccumErrorStart = ark_mem->tn; + ark_mem->AccumError = ZERO; return (ARK_SUCCESS); } @@ -2447,7 +2447,7 @@ int ARKodeGetAccumulatedError(void* arkode_mem, sunrealtype* accum_error) /* Get number of steps since last accumulated error reset (set floor of 1 to safeguard against division-by-zero) */ - long int steps = SUNMAX(1, ark_mem->nst - ark_mem->AccumErrorStep); + sunrealtype time_interval = ark_mem->tcur - ark_mem->AccumErrorStart; /* Fill output based on error accumulation type */ if (ark_mem->AccumErrorType == ARK_ACCUMERROR_MAX) @@ -2460,7 +2460,7 @@ int ARKodeGetAccumulatedError(void* arkode_mem, sunrealtype* accum_error) } else if (ark_mem->AccumErrorType == ARK_ACCUMERROR_AVG) { - *accum_error = ark_mem->AccumError * ark_mem->reltol / ((sunrealtype)steps); + *accum_error = ark_mem->AccumError * ark_mem->reltol / time_interval; } else { diff --git a/src/arkode/arkode_mri_tables.c b/src/arkode/arkode_mri_tables.c index 11d8e14b17..ecae4773f3 100644 --- a/src/arkode/arkode_mri_tables.c +++ b/src/arkode/arkode_mri_tables.c @@ -69,7 +69,8 @@ MRIStepCoupling MRIStepCoupling_LoadTableByName(const char* method) /*--------------------------------------------------------------- Routine to allocate an empty MRIStepCoupling structure ---------------------------------------------------------------*/ -MRIStepCoupling MRIStepCoupling_Alloc(int nmat, int stages, MRISTEP_METHOD_TYPE type) +MRIStepCoupling MRIStepCoupling_Alloc(int nmat, int stages, + MRISTEP_METHOD_TYPE type) { int i, j; sunbooleantype hasOmegas, hasGammas; @@ -894,8 +895,7 @@ int mriStepCoupling_GetStageMap(MRIStepCoupling MRIC, int* stage_map, /* Number of stage RHS vectors active */ *nstages_active = MRIC->stages; - /* Check if a stage corresponds to a column of zeros for all coupling - * matrices by computing the column sums */ + /* Create an identity map (all columns are non-zero) */ for (j = 0; j < MRIC->stages; j++) { stage_map[j] = j; } return (ARK_SUCCESS); } diff --git a/src/arkode/arkode_mristep.c b/src/arkode/arkode_mristep.c index 8a8fbea714..fc42735c13 100644 --- a/src/arkode/arkode_mristep.c +++ b/src/arkode/arkode_mristep.c @@ -143,7 +143,6 @@ void* MRIStepCreate(ARKRhsFn fse, ARKRhsFn fsi, sunrealtype t0, N_Vector y0, ark_mem->step_setforcing = mriStep_SetInnerForcing; ark_mem->step_supports_adaptive = SUNTRUE; ark_mem->step_supports_implicit = SUNTRUE; - ark_mem->step_supports_forcing = SUNTRUE; ark_mem->step_mem = (void*)step_mem; /* Set default values for optional inputs */ @@ -886,7 +885,7 @@ int mriStep_GetGammas(ARKodeMem ark_mem, sunrealtype* gamma, sunrealtype* gamrat With initialization types FIRST_INIT this routine: - sets/checks the coefficient tables to be used - - allocates any internal memory that depend on the MRI method + - allocates any internal memory that depends on the MRI method structure or solver options With other initialization types, this routine does nothing. @@ -952,6 +951,10 @@ int mriStep_Init(ARKodeMem ark_mem, int init_type) case MRISTEP_IMEX: ark_mem->step = mriStep_TakeStepMRIGARK; break; case MRISTEP_MERK: ark_mem->step = mriStep_TakeStepMERK; break; case MRISTEP_MRISR: ark_mem->step = mriStep_TakeStepMRISR; break; + default: + arkProcessError(ark_mem, ARK_ILL_INPUT, __LINE__, __func__, __FILE__, + "Unknown method type"); + return (ARK_ILL_INPUT); } /* Retrieve/store method and embedding orders now that tables are finalized */ @@ -962,7 +965,7 @@ int mriStep_Init(ARKodeMem ark_mem, int init_type) /* Ensure that if adaptivity or error accumulation is enabled, then method includes embedding coefficients */ if ((!ark_mem->fixedstep || (ark_mem->AccumErrorType != ARK_ACCUMERROR_NONE)) && - (step_mem->p == 0)) + (step_mem->p <= 0)) { arkProcessError(ark_mem, ARK_ILL_INPUT, __LINE__, __func__, __FILE__, "Temporal error estimation cannot be performed without embedding coefficients"); @@ -1247,9 +1250,8 @@ int mriStep_Init(ARKodeMem ark_mem, int init_type) if (ark_mem->fixedstep) { - /* Non-adaptive controller: user must have supplied initial step - size, and indicated fixed time stepping */ - if ((ark_mem->hin == ZERO) || (ark_mem->fixedstep == SUNFALSE)) + /* Fixed step sizes: user must supply the initial step size */ + if (ark_mem->hin == ZERO) { arkProcessError(ark_mem, ARK_ILL_INPUT, __LINE__, __func__, __FILE__, "Timestep adaptivity disabled, but missing user-defined fixed stepsize"); @@ -1259,7 +1261,7 @@ int mriStep_Init(ARKodeMem ark_mem, int init_type) else { /* ensure that a compatible adaptivity controller is provided */ - if ((adapt_type != SUN_ADAPTCONTROLLER_MRI_TOL) && + if ((adapt_type != SUN_ADAPTCONTROLLER_MRI_H_TOL) && (adapt_type != SUN_ADAPTCONTROLLER_H)) { arkProcessError(ark_mem, ARK_ILL_INPUT, __LINE__, __func__, __FILE__, @@ -1270,7 +1272,7 @@ int mriStep_Init(ARKodeMem ark_mem, int init_type) /* Controller provides adaptivity (at least at the slow time scale): - verify that the MRI method includes an embedding, and - estimate initial slow step size (store in ark_mem->hin) */ - if (step_mem->MRIC->p == 0) + if (step_mem->MRIC->p <= 0) { arkProcessError(ark_mem, ARK_ILL_INPUT, __LINE__, __func__, __FILE__, "Timestep adaptivity enabled, but non-embedded MRI table specified"); @@ -1297,13 +1299,13 @@ int mriStep_Init(ARKodeMem ark_mem, int init_type) } /* Perform additional setup for (H,tol) controller */ - if (adapt_type == SUN_ADAPTCONTROLLER_MRI_TOL) + if (adapt_type == SUN_ADAPTCONTROLLER_MRI_H_TOL) { /* Verify that adaptivity type is supported by inner stepper */ if (!mriStepInnerStepper_SupportsRTolAdaptivity(step_mem->stepper)) { arkProcessError(ark_mem, ARK_ILL_INPUT, __LINE__, __func__, - __FILE__, "MRI-TOL SUNAdaptController provided, but unsupported by inner stepper"); + __FILE__, "MRI H-TOL SUNAdaptController provided, but unsupported by inner stepper"); return (ARK_ILL_INPUT); } @@ -1692,7 +1694,6 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt /* initial time and state for step */ N_VScale(ONE, ark_mem->yn, ark_mem->ycur); - t0 = ark_mem->tn; /* determine whether embedding stage is needed */ do_embedding = !ark_mem->fixedstep || @@ -1702,11 +1703,11 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt and send appropriate control parameter to the fast integrator */ adapt_type = SUNAdaptController_GetType(ark_mem->hadapt_mem->hcontroller); need_inner_dsm = SUNFALSE; - if (adapt_type == SUN_ADAPTCONTROLLER_MRI_TOL) + if (adapt_type == SUN_ADAPTCONTROLLER_MRI_H_TOL) { need_inner_dsm = SUNTRUE; step_mem->inner_dsm = ZERO; - retval = mriStepInnerStepper_ResetError(step_mem->stepper); + retval = mriStepInnerStepper_ResetAccumulatedError(step_mem->stepper); if (retval != ARK_SUCCESS) { arkProcessError(ark_mem, ARK_INNERSTEP_FAIL, __LINE__, __func__, __FILE__, @@ -1728,7 +1729,8 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt to the beginning of this step (in case of recomputation) */ if (!ark_mem->fixedstep) { - retval = mriStepInnerStepper_Reset(step_mem->stepper, t0, ark_mem->yn); + retval = mriStepInnerStepper_Reset(step_mem->stepper, ark_mem->tn, + ark_mem->yn); if (retval != ARK_SUCCESS) { arkProcessError(ark_mem, ARK_INNERSTEP_FAIL, __LINE__, __func__, __FILE__, @@ -1758,7 +1760,7 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt /* compute the explicit component */ if (step_mem->explicit_rhs) { - retval = step_mem->fse(t0, ark_mem->yn, step_mem->Fse[0], + retval = step_mem->fse(ark_mem->tn, ark_mem->yn, step_mem->Fse[0], ark_mem->user_data); step_mem->nfse++; if (retval) { return ARK_RHSFUNC_FAIL; } @@ -1767,7 +1769,7 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt /* compute the implicit component */ if (step_mem->implicit_rhs) { - retval = step_mem->fsi(t0, ark_mem->yn, step_mem->Fsi[0], + retval = step_mem->fsi(ark_mem->tn, ark_mem->yn, step_mem->Fsi[0], ark_mem->user_data); step_mem->nfsi++; if (retval) { return ARK_RHSFUNC_FAIL; } @@ -1779,7 +1781,7 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt step_mem->cvals[0] = ONE; step_mem->Xvecs[0] = step_mem->Fse[0]; nvec = 1; - mriStep_ApplyForcing(step_mem, t0, ONE, &nvec); + mriStep_ApplyForcing(step_mem, ark_mem->tn, ONE, &nvec); N_VLinearCombination(nvec, step_mem->cvals, step_mem->Xvecs, step_mem->Fse[0]); } @@ -1788,7 +1790,7 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt step_mem->cvals[0] = ONE; step_mem->Xvecs[0] = step_mem->Fsi[0]; nvec = 1; - mriStep_ApplyForcing(step_mem, t0, ONE, &nvec); + mriStep_ApplyForcing(step_mem, ark_mem->tn, ONE, &nvec); N_VLinearCombination(nvec, step_mem->cvals, step_mem->Xvecs, step_mem->Fsi[0]); } @@ -1798,27 +1800,24 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt #ifdef SUNDIALS_DEBUG printf(" MRIStep step %li, stage 0, h = %" RSYM ", t_n = %" RSYM "\n", - ark_mem->nst, ark_mem->h, t0); + ark_mem->nst, ark_mem->h, ark_mem->tn); #endif #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRIGARK", "slow stage", - "z[0] =", ""); - N_VPrintFile(ark_mem->yn, ARK_LOGGER->debug_fp); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "slow stage", "z_0(:) =", ""); + N_VPrintFile(ark_mem->ycur, ARK_LOGGER->debug_fp); if (step_mem->explicit_rhs) { - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRIGARK", "slow explicit RHS", - "Fse[0] =", ""); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "slow explicit RHS", "Fse_0(:) =", ""); N_VPrintFile(step_mem->Fse[0], ARK_LOGGER->debug_fp); } if (step_mem->implicit_rhs) { - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRIGARK", "slow implicit RHS", - "Fsi[0] =", ""); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "slow implicit RHS", "Fsi_0(:) =", ""); N_VPrintFile(step_mem->Fsi[0], ARK_LOGGER->debug_fp); } #endif @@ -1836,7 +1835,7 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt /* Solver diagnostics reporting */ #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRIGARK", "start-stage", + "ARKODE::mriStep_TakeStep", "start-stage", "step = %li, stage = %i, stage type = %d, h = %" RSYM ", tcur = %" RSYM, ark_mem->nst, is, step_mem->stagetypes[is], ark_mem->h, @@ -1850,11 +1849,7 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt { case (MRISTAGE_ERK_FAST): retval = mriStep_ComputeInnerForcing(ark_mem, step_mem, is, t0, tf); - if (retval != ARK_SUCCESS) - { - *nflagPtr = CONV_FAIL; - break; - } + if (retval != ARK_SUCCESS) { return retval; } retval = mriStep_StageERKFast(ark_mem, step_mem, t0, tf, ark_mem->ycur, ark_mem->tempv2, SUNFALSE, need_inner_dsm); if (retval != ARK_SUCCESS) { *nflagPtr = CONV_FAIL; } @@ -1874,8 +1869,7 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRIGARK", "slow stage", - "z[%i] =", is); + "ARKODE::mriStep_TakeStep", "slow stage", "z_%i(:) =", is); N_VPrintFile(ark_mem->ycur, ARK_LOGGER->debug_fp); #endif @@ -1883,7 +1877,8 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt if ((ark_mem->ProcessStage != NULL) && (step_mem->stagetypes[is] != MRISTAGE_STIFF_ACC)) { - retval = ark_mem->ProcessStage(tf, ark_mem->ycur, ark_mem->user_data); + retval = ark_mem->ProcessStage(ark_mem->tcur, ark_mem->ycur, + ark_mem->user_data); if (retval != 0) { return (ARK_POSTPROCESS_STAGE_FAIL); } } @@ -1922,7 +1917,7 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt /* store explicit slow rhs */ if (step_mem->explicit_rhs) { - retval = step_mem->fse(tf, ark_mem->ycur, + retval = step_mem->fse(ark_mem->tcur, ark_mem->ycur, step_mem->Fse[step_mem->stage_map[is]], ark_mem->user_data); step_mem->nfse++; @@ -1942,8 +1937,8 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRIGARK", - "slow explicit RHS", "Fse[%i] =", is); + "ARKODE::mriStep_TakeStep", "slow explicit RHS", + "Fse_%i(:) =", is); N_VPrintFile(step_mem->Fse[step_mem->stage_map[is]], ARK_LOGGER->debug_fp); #endif @@ -1955,7 +1950,7 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt if (!step_mem->deduce_rhs || (step_mem->stagetypes[is] != MRISTAGE_DIRK_NOFAST)) { - retval = step_mem->fsi(tf, ark_mem->ycur, + retval = step_mem->fsi(ark_mem->tcur, ark_mem->ycur, step_mem->Fsi[step_mem->stage_map[is]], ark_mem->user_data); step_mem->nfsi++; @@ -1983,8 +1978,8 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRIGARK", - "slow implicit RHS", "Fsi[%i] =", is); + "ARKODE::mriStep_TakeStep", "slow implicit RHS", + "Fsi_%i(:) =", is); N_VPrintFile(step_mem->Fsi[step_mem->stage_map[is]], ARK_LOGGER->debug_fp); #endif @@ -2001,9 +1996,8 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt } /* loop over stages */ #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRIGARK", "updated solution", - "ycur =", ""); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "updated solution", "ycur(:) =", ""); N_VPrintFile(ark_mem->ycur, ARK_LOGGER->debug_fp); #endif @@ -2028,7 +2022,7 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt /* Solver diagnostics reporting */ #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRIGARK", "embedding-stage", + "ARKODE::mriStep_TakeStep", "embedding-stage", "step = %li, stage = %i, stage type = %d, h = %" RSYM ", tcur = %" RSYM, ark_mem->nst, step_mem->stages, @@ -2044,11 +2038,7 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt case (MRISTAGE_ERK_FAST): retval = mriStep_ComputeInnerForcing(ark_mem, step_mem, step_mem->stages, t0, tf); - if (retval != ARK_SUCCESS) - { - *nflagPtr = CONV_FAIL; - break; - } + if (retval != ARK_SUCCESS) { return retval; } retval = mriStep_StageERKFast(ark_mem, step_mem, t0, tf, ark_mem->ycur, ark_mem->tempv2, SUNTRUE, SUNFALSE); if (retval != ARK_SUCCESS) { *nflagPtr = CONV_FAIL; } @@ -2068,9 +2058,8 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt if (retval != ARK_SUCCESS) { return (retval); } #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRIGARK", "embedded solution", - "ytilde =", ""); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "embedded solution", "ytilde(:) =", ""); N_VPrintFile(ark_mem->ycur, ARK_LOGGER->debug_fp); #endif @@ -2087,10 +2076,9 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt /* Solver diagnostics reporting */ #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_DEBUG - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRIGARK", "end-step", - "step = %li, t = %" RSYM ", h = %" RSYM ", dsm = %" RSYM, - ark_mem->nst, ark_mem->tn, ark_mem->h, *dsmPtr); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "error-test", "step = %li, h = %" RSYM ", dsm = %" RSYM, + ark_mem->nst, ark_mem->h, *dsmPtr); #endif return (ARK_SUCCESS); @@ -2168,11 +2156,11 @@ int mriStep_TakeStepMRISR(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) and send appropriate control parameter to the fast integrator */ adapt_type = SUNAdaptController_GetType(ark_mem->hadapt_mem->hcontroller); need_inner_dsm = SUNFALSE; - if (adapt_type == SUN_ADAPTCONTROLLER_MRI_TOL) + if (adapt_type == SUN_ADAPTCONTROLLER_MRI_H_TOL) { need_inner_dsm = SUNTRUE; step_mem->inner_dsm = ZERO; - retval = mriStepInnerStepper_ResetError(step_mem->stepper); + retval = mriStepInnerStepper_ResetAccumulatedError(step_mem->stepper); if (retval != ARK_SUCCESS) { arkProcessError(ark_mem, ARK_INNERSTEP_FAIL, __LINE__, __func__, __FILE__, @@ -2264,23 +2252,21 @@ int mriStep_TakeStepMRISR(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #endif #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRISR", "slow stage", "z[0] =", ""); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "slow stage", "z_0(:) =", ""); N_VPrintFile(ark_mem->yn, ARK_LOGGER->debug_fp); if (step_mem->explicit_rhs) { - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRISR", "slow explicit RHS", - "Fse[0] =", ""); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "slow explicit RHS", "Fse_0(:) =", ""); N_VPrintFile(step_mem->Fse[0], ARK_LOGGER->debug_fp); } if (step_mem->implicit_rhs) { - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRISR", "slow implicit RHS", - "Fsi[0] =", ""); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "slow implicit RHS", "Fsi_0(:) =", ""); N_VPrintFile(step_mem->Fsi[0], ARK_LOGGER->debug_fp); } #endif @@ -2293,10 +2279,9 @@ int mriStep_TakeStepMRISR(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) { /* Solver diagnostics reporting */ #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_DEBUG - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRISR", "start-stage", - "step = %li, stage = %i, h = %" RSYM, ark_mem->nst, - stage, ark_mem->h); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "start-stage", "step = %li, stage = %i, h = %" RSYM, + ark_mem->nst, stage, ark_mem->h); #endif /* Determine if this is an "embedding" or "solution" stage */ @@ -2326,11 +2311,7 @@ int mriStep_TakeStepMRISR(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) step_mem->explicit_rhs = SUNTRUE; retval = mriStep_ComputeInnerForcing(ark_mem, step_mem, stage, ark_mem->tn, ark_mem->tn + cstage * ark_mem->h); - if (retval != ARK_SUCCESS) - { - *nflagPtr = CONV_FAIL; - break; - } + if (retval != ARK_SUCCESS) { return retval; } step_mem->implicit_rhs = store_imprhs; step_mem->explicit_rhs = store_exprhs; @@ -2377,8 +2358,8 @@ int mriStep_TakeStepMRISR(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRISR", "predictor", - "zpred =", ""); + "ARKODE::mriStep_TakeStep", "predictor", + "zpred(:) =", ""); N_VPrintFile(step_mem->zpred, ARK_LOGGER->debug_fp); #endif @@ -2398,8 +2379,8 @@ int mriStep_TakeStepMRISR(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRISR", "rhs data", - "sdata =", ""); + "ARKODE::mriStep_TakeStep", "rhs data", + "sdata(:) =", ""); N_VPrintFile(step_mem->sdata, ARK_LOGGER->debug_fp); #endif @@ -2431,9 +2412,8 @@ int mriStep_TakeStepMRISR(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) } #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRISR", "slow stage", - "z[%i] =", stage); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "slow stage", "z_%i(:) =", stage); N_VPrintFile(ark_mem->ycur, ARK_LOGGER->debug_fp); #endif @@ -2479,8 +2459,8 @@ int mriStep_TakeStepMRISR(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRISR", "slow explicit RHS", - "Fse[%i] =", stage); + "ARKODE::mriStep_TakeStep", "slow explicit RHS", + "Fse_%i(:) =", stage); N_VPrintFile(step_mem->Fse[stage], ARK_LOGGER->debug_fp); #endif } @@ -2517,8 +2497,8 @@ int mriStep_TakeStepMRISR(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRISR", "slow implicit RHS", - "Fsi[%i] =", stage); + "ARKODE::mriStep_TakeStep", "slow implicit RHS", + "Fsi_%i(:) =", stage); N_VPrintFile(step_mem->Fsi[stage], ARK_LOGGER->debug_fp); #endif } @@ -2537,9 +2517,8 @@ int mriStep_TakeStepMRISR(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) } /* loop over stages */ #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRISR", "updated solution", - "ycur =", ""); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "updated solution", "ycur(:) =", ""); N_VPrintFile(ark_mem->ycur, ARK_LOGGER->debug_fp); #endif @@ -2555,8 +2534,8 @@ int mriStep_TakeStepMRISR(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) /* Solver diagnostics reporting */ #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_DEBUG - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMRISR", "end-step", + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "end-step", "step = %li, t = %" RSYM ", h = %" RSYM ", dsm = %" RSYM, ark_mem->nst, ark_mem->tn, ark_mem->h, *dsmPtr); #endif @@ -2638,11 +2617,11 @@ int mriStep_TakeStepMERK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) and send appropriate control parameter to the fast integrator */ adapt_type = SUNAdaptController_GetType(ark_mem->hadapt_mem->hcontroller); need_inner_dsm = SUNFALSE; - if (adapt_type == SUN_ADAPTCONTROLLER_MRI_TOL) + if (adapt_type == SUN_ADAPTCONTROLLER_MRI_H_TOL) { need_inner_dsm = SUNTRUE; step_mem->inner_dsm = ZERO; - retval = mriStepInnerStepper_ResetError(step_mem->stepper); + retval = mriStepInnerStepper_ResetAccumulatedError(step_mem->stepper); if (retval != ARK_SUCCESS) { arkProcessError(ark_mem, ARK_INNERSTEP_FAIL, __LINE__, __func__, __FILE__, @@ -2702,17 +2681,13 @@ int mriStep_TakeStepMERK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #endif #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMERK", "slow stage", "z[0] =", ""); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "slow stage", "z_0(:) =", ""); N_VPrintFile(ark_mem->yn, ARK_LOGGER->debug_fp); - if (step_mem->explicit_rhs) - { - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMERK", "slow explicit RHS", - "Fse[0] =", ""); - N_VPrintFile(step_mem->Fse[0], ARK_LOGGER->debug_fp); - } + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "slow explicit RHS", "Fse_0(:) =", ""); + N_VPrintFile(step_mem->Fse[0], ARK_LOGGER->debug_fp); #endif /* The first stage is the previous time-step solution, so its RHS @@ -2723,10 +2698,9 @@ int mriStep_TakeStepMERK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) { /* Solver diagnostics reporting */ #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_DEBUG - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMERK", "start-stage", - "step = %li, stage group = %i, h = %" RSYM, ark_mem->nst, - ig, ark_mem->h); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "start-stage", "step = %li, stage group = %i, h = %" RSYM, + ark_mem->nst, ig, ark_mem->h); #endif /* Set up fast RHS for this stage group */ @@ -2792,8 +2766,8 @@ int mriStep_TakeStepMERK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMERK", "slow stage", - "z[%i] =", stage); + "ARKODE::mriStep_TakeStep", "slow stage", + "z_%i(:) =", stage); N_VPrintFile(ark_mem->ycur, ARK_LOGGER->debug_fp); #endif @@ -2837,8 +2811,8 @@ int mriStep_TakeStepMERK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMERK", "slow explicit RHS", - "Fse[%i] =", is); + "ARKODE::mriStep_TakeStep", "slow explicit RHS", + "Fse_%i(:) =", is); N_VPrintFile(step_mem->Fse[stage], ARK_LOGGER->debug_fp); #endif } @@ -2851,9 +2825,8 @@ int mriStep_TakeStepMERK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) } /* loop over stage groups */ #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMERK", "updated solution", - "ycur =", ""); + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "updated solution", "ycur(:) =", ""); N_VPrintFile(ark_mem->ycur, ARK_LOGGER->debug_fp); #endif @@ -2867,8 +2840,8 @@ int mriStep_TakeStepMERK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) /* Solver diagnostics reporting */ #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_DEBUG - SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, - "ARKODE::mriStep_TakeStepMERK", "end-step", + SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_TakeStep", + "end-step", "step = %li, t = %" RSYM ", h = %" RSYM ", dsm = %" RSYM, ark_mem->nst, ark_mem->tn, ark_mem->h, *dsmPtr); #endif @@ -3357,10 +3330,10 @@ int mriStep_StageERKFast(ARKodeMem ark_mem, ARKodeMRIStepMem step_mem, if (get_inner_dsm) { /* if the fast integrator uses adaptive steps, retrieve the error estimate */ - if (adapt_type == SUN_ADAPTCONTROLLER_MRI_TOL) + if (adapt_type == SUN_ADAPTCONTROLLER_MRI_H_TOL) { - retval = mriStepInnerStepper_GetError(step_mem->stepper, - &(step_mem->inner_dsm)); + retval = mriStepInnerStepper_GetAccumulatedError(step_mem->stepper, + &(step_mem->inner_dsm)); if (retval != ARK_SUCCESS) { arkProcessError(ark_mem, ARK_INNERSTEP_FAIL, __LINE__, __func__, @@ -3489,7 +3462,7 @@ int mriStep_StageDIRKNoFast(ARKodeMem ark_mem, ARKodeMRIStepMem step_mem, #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_StageDIRKNoFast", "predictor", - "zpred =", ""); + "zpred(:) =", ""); N_VPrintFile(step_mem->zpred, ARK_LOGGER->debug_fp); #endif @@ -3505,7 +3478,7 @@ int mriStep_StageDIRKNoFast(ARKodeMem ark_mem, ARKodeMRIStepMem step_mem, #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_StageDIRKNoFast", "rhs data", - "sdata =", ""); + "sdata(:) =", ""); N_VPrintFile(step_mem->sdata, ARK_LOGGER->debug_fp); #endif @@ -3678,7 +3651,7 @@ int mriStep_ComputeInnerForcing(SUNDIALS_MAYBE_UNUSED ARKodeMem ark_mem, { SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::mriStep_ComputeInnerForcing", "forcing", - "forcing[%i] =", k); + "forcing_%i(:) =", k); N_VPrintFile(step_mem->stepper->forcing[k], ARK_LOGGER->debug_fp); } #endif @@ -3958,7 +3931,7 @@ int mriStep_StageSetup(ARKodeMem ark_mem) ARKTimestepFullRHSFn. This is only used to determine an initial slow time-step size to use when one is not specified by the user (i.e., mode should correspond with - ARK_FULLRHS_OTHER. + ARK_FULLRHS_START. ---------------------------------------------------------------*/ int mriStep_SlowRHS(ARKodeMem ark_mem, sunrealtype t, N_Vector y, N_Vector f, SUNDIALS_MAYBE_UNUSED int mode) @@ -4030,38 +4003,6 @@ int mriStep_SlowRHS(ARKodeMem ark_mem, sunrealtype t, N_Vector y, N_Vector f, return (ARK_SUCCESS); } -/*--------------------------------------------------------------- - mriStep_FastRHS: - - This is just a wrapper to call the fast RHS function, - f(t,y) = ff(t,y), with API matching ARKTimestepFullRHSFn. This - is only used to determine an initial fast time-step size to use - when one is not specified by the user and H-H MRI time step - adaptivity is enabled. - ---------------------------------------------------------------*/ -int mriStep_FastRHS(ARKodeMem ark_mem, sunrealtype t, N_Vector y, N_Vector f, - SUNDIALS_MAYBE_UNUSED int mode) -{ - ARKodeMRIStepMem step_mem; - int retval; - - /* access ARKodeMRIStepMem structure */ - retval = mriStep_AccessStepMem(ark_mem, __func__, &step_mem); - if (retval != ARK_SUCCESS) { return (retval); } - - /* call ff */ - retval = mriStepInnerStepper_FullRhs(step_mem->stepper, t, y, f, - ARK_FULLRHS_OTHER); - if (retval != ARK_SUCCESS) - { - arkProcessError(ark_mem, ARK_RHSFUNC_FAIL, __LINE__, __func__, __FILE__, - MSG_ARK_RHSFUNC_FAILED, t); - return (ARK_RHSFUNC_FAIL); - } - - return (ARK_SUCCESS); -} - /*--------------------------------------------------------------- mriStep_Hin @@ -4390,7 +4331,7 @@ int mriStepInnerStepper_Evolve(MRIStepInnerStepper stepper, sunrealtype t0, if (stepper->ops == NULL) { return ARK_ILL_INPUT; } if (stepper->ops->evolve == NULL) { return ARK_ILL_INPUT; } -#if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_INFO +#if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_DEBUG SUNLogger_QueueMsg(stepper->sunctx->logger, SUN_LOGLEVEL_INFO, "ARKODE::mriStepInnerStepper_Evolve", "start-inner-evolve", "t0 = %" RSYM ", tout = %" RSYM, t0, tout); @@ -4398,7 +4339,7 @@ int mriStepInnerStepper_Evolve(MRIStepInnerStepper stepper, sunrealtype t0, stepper->last_flag = stepper->ops->evolve(stepper, t0, tout, y); -#if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_INFO +#if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_DEBUG SUNLogger_QueueMsg(stepper->sunctx->logger, SUN_LOGLEVEL_INFO, "ARKODE::mriStepInnerStepper_Evolve", "end-inner-evolve", "flag = %i", stepper->last_flag); @@ -4447,8 +4388,8 @@ int mriStepInnerStepper_Reset(MRIStepInnerStepper stepper, sunrealtype tR, } /* Gets the inner (fast) stepper accumulated error */ -int mriStepInnerStepper_GetError(MRIStepInnerStepper stepper, - sunrealtype* accum_error) +int mriStepInnerStepper_GetAccumulatedError(MRIStepInnerStepper stepper, + sunrealtype* accum_error) { if (stepper == NULL) { return ARK_ILL_INPUT; } if (stepper->ops == NULL) { return ARK_ILL_INPUT; } @@ -4462,7 +4403,7 @@ int mriStepInnerStepper_GetError(MRIStepInnerStepper stepper, } /* Resets the inner (fast) stepper accumulated error */ -int mriStepInnerStepper_ResetError(MRIStepInnerStepper stepper) +int mriStepInnerStepper_ResetAccumulatedError(MRIStepInnerStepper stepper) { if (stepper == NULL) { return ARK_ILL_INPUT; } if (stepper->ops == NULL) { return ARK_ILL_INPUT; } @@ -4671,7 +4612,7 @@ void mriStep_ApplyForcing(ARKodeMRIStepMem step_mem, sunrealtype t, Sets an array of coefficient vectors for a time-dependent external polynomial forcing term in the ODE RHS i.e., y' = f(t,y) + p(t). This function is primarily intended for using MRIStep as an inner integrator within another - [outer] instance of MRIStep, where this instance is is used to solve a + [outer] instance of MRIStep, where this instance is used to solve a modified ODE at a fast time scale. The polynomial is of the form p(t) = sum_{i = 0}^{nvecs - 1} forcing[i] * ((t - tshift) / (tscale))^i diff --git a/src/arkode/arkode_mristep_controller.c b/src/arkode/arkode_mristep_controller.c index dc71c69762..f0e399b58d 100644 --- a/src/arkode/arkode_mristep_controller.c +++ b/src/arkode/arkode_mristep_controller.c @@ -41,7 +41,7 @@ SUNAdaptController SUNAdaptController_MRIStep(ARKodeMem ark_mem, /* Return with failure if input controller is NULL or has unsupported type */ if (CMRI == NULL) { return (NULL); } - if (SUNAdaptController_GetType(CMRI) != SUN_ADAPTCONTROLLER_MRI_TOL) + if (SUNAdaptController_GetType(CMRI) != SUN_ADAPTCONTROLLER_MRI_H_TOL) { return (NULL); } @@ -111,9 +111,9 @@ SUNErrCode SUNAdaptController_UpdateH_MRIStep(SUNAdaptController C, if ((ark_mem == NULL) || (step_mem == NULL)) { return SUN_ERR_MEM_FAIL; } /* Update MRI controller */ - SUNErrCode retval = SUNAdaptController_UpdateMRITol(MRICONTROL_C(C), H, - step_mem->inner_rtol_factor, - DSM, step_mem->inner_dsm); + SUNErrCode retval = SUNAdaptController_UpdateMRIHTol(MRICONTROL_C(C), H, + step_mem->inner_rtol_factor, + DSM, step_mem->inner_dsm); if (retval != SUN_SUCCESS) { return (retval); } /* Update inner controller parameter to most-recent prediction */ diff --git a/src/arkode/arkode_mristep_impl.h b/src/arkode/arkode_mristep_impl.h index b84a36232c..8a39fd3887 100644 --- a/src/arkode/arkode_mristep_impl.h +++ b/src/arkode/arkode_mristep_impl.h @@ -109,7 +109,7 @@ typedef struct ARKodeMRIStepMemRec sunrealtype nlscoef; /* coefficient in nonlin. convergence test */ int msbp; /* positive => max # steps between lsetup - negative => call at each Newton iter */ + negative => call at each Newton iter */ long int nstlp; /* step number of last setup call */ int maxcor; /* max num iterations for solving the @@ -288,8 +288,6 @@ int mriStep_NlsInit(ARKodeMem ark_mem); int mriStep_Nls(ARKodeMem ark_mem, int nflag); int mriStep_SlowRHS(ARKodeMem ark_mem, sunrealtype t, N_Vector y, N_Vector f, int mode); -int mriStep_FastRHS(ARKodeMem ark_mem, sunrealtype t, N_Vector y, N_Vector f, - int mode); int mriStep_Hin(ARKodeMem ark_mem, sunrealtype tcur, sunrealtype tout, N_Vector fcur, sunrealtype* h); void mriStep_ApplyForcing(ARKodeMRIStepMem step_mem, sunrealtype t, @@ -314,9 +312,9 @@ int mriStepInnerStepper_FullRhs(MRIStepInnerStepper stepper, sunrealtype t, N_Vector y, N_Vector f, int mode); int mriStepInnerStepper_Reset(MRIStepInnerStepper stepper, sunrealtype tR, N_Vector yR); -int mriStepInnerStepper_GetError(MRIStepInnerStepper stepper, - sunrealtype* accum_error); -int mriStepInnerStepper_ResetError(MRIStepInnerStepper stepper); +int mriStepInnerStepper_GetAccumulatedError(MRIStepInnerStepper stepper, + sunrealtype* accum_error); +int mriStepInnerStepper_ResetAccumulatedError(MRIStepInnerStepper stepper); int mriStepInnerStepper_SetRTol(MRIStepInnerStepper stepper, sunrealtype rtol); int mriStepInnerStepper_AllocVecs(MRIStepInnerStepper stepper, int count, N_Vector tmpl); @@ -344,7 +342,7 @@ int mriStep_RKCoeffs(MRIStepCoupling MRIC, int is, int* stage_map, rate context, leveraging MRIStep-specific knowledge of the slow+fast time scale relationship to CALL multi-rate controller functions (e.g., EstimateMRISteps, EstimateStepTol, UpdateMRIH, - and UpdateMRITol) provided by the underlying multi-rate + and UpdateMRIHTol) provided by the underlying multi-rate controller. ===============================================================*/ diff --git a/src/arkode/arkode_mristep_io.c b/src/arkode/arkode_mristep_io.c index c99e7ff0af..ed4eab3d5c 100644 --- a/src/arkode/arkode_mristep_io.c +++ b/src/arkode/arkode_mristep_io.c @@ -242,7 +242,7 @@ int mriStep_SetAdaptController(ARKodeMem ark_mem, SUNAdaptController C) SUNAdaptController_Type ctype = SUNAdaptController_GetType(C); /* If this does not have MRI type, then just pass to ARKODE */ - if (ctype != SUN_ADAPTCONTROLLER_MRI_TOL) + if (ctype != SUN_ADAPTCONTROLLER_MRI_H_TOL) { return (arkReplaceAdaptController(ark_mem, C)); } @@ -286,7 +286,8 @@ int mriStep_SetUserData(ARKodeMem ark_mem, void* user_data) int mriStep_SetDefaults(ARKodeMem ark_mem) { ARKodeMRIStepMem step_mem; - sunindextype lenrw, leniw; + sunindextype Clenrw, Cleniw; + long int lenrw, leniw; int retval; /* access ARKodeMRIStepMem structure */ @@ -319,12 +320,47 @@ int mriStep_SetDefaults(ARKodeMem ark_mem) /* Remove pre-existing coupling table */ if (step_mem->MRIC) { - MRIStepCoupling_Space(step_mem->MRIC, &leniw, &lenrw); - ark_mem->lrw -= lenrw; - ark_mem->liw -= leniw; + MRIStepCoupling_Space(step_mem->MRIC, &Cleniw, &Clenrw); + ark_mem->lrw -= Clenrw; + ark_mem->liw -= Cleniw; MRIStepCoupling_Free(step_mem->MRIC); } step_mem->MRIC = NULL; + + /* Remove pre-existing SUNAdaptController object, and replace with "I" */ + if (ark_mem->hadapt_mem->owncontroller) + { + retval = SUNAdaptController_Space(ark_mem->hadapt_mem->hcontroller, &lenrw, + &leniw); + if (retval == SUN_SUCCESS) + { + ark_mem->liw -= leniw; + ark_mem->lrw -= lenrw; + } + retval = SUNAdaptController_Destroy(ark_mem->hadapt_mem->hcontroller); + ark_mem->hadapt_mem->owncontroller = SUNFALSE; + if (retval != SUN_SUCCESS) + { + arkProcessError(ark_mem, ARK_MEM_FAIL, __LINE__, __func__, __FILE__, + "SUNAdaptController_Destroy failure"); + return (ARK_MEM_FAIL); + } + } + ark_mem->hadapt_mem->hcontroller = SUNAdaptController_I(ark_mem->sunctx); + if (ark_mem->hadapt_mem->hcontroller == NULL) + { + arkProcessError(ark_mem, ARK_MEM_FAIL, __LINE__, __func__, __FILE__, + "SUNAdaptController_I allocation failure"); + return (ARK_MEM_FAIL); + } + ark_mem->hadapt_mem->owncontroller = SUNTRUE; + retval = SUNAdaptController_Space(ark_mem->hadapt_mem->hcontroller, &lenrw, + &leniw); + if (retval == SUN_SUCCESS) + { + ark_mem->liw += leniw; + ark_mem->lrw += lenrw; + } return (ARK_SUCCESS); } diff --git a/src/arkode/arkode_relaxation.c b/src/arkode/arkode_relaxation.c index e2307448b9..c911077580 100644 --- a/src/arkode/arkode_relaxation.c +++ b/src/arkode/arkode_relaxation.c @@ -357,7 +357,7 @@ static int arkRelaxSolve(ARKodeMem ark_mem, ARKodeRelaxMem relax_mem, #ifdef SUNDIALS_LOGGING_EXTRA_DEBUG SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_DEBUG, "ARKODE::arkRelaxSolve", - "compute delta y", "delta_y =", ""); + "compute delta y", "delta_y(:) =", ""); N_VPrintFile(ark_mem->tempv2, ARK_LOGGER->debug_fp); #endif diff --git a/src/sunadaptcontroller/mrihtol/fmod_int32/CMakeLists.txt b/src/sunadaptcontroller/mrihtol/fmod_int32/CMakeLists.txt index ecd215bc10..befc6cb2ed 100644 --- a/src/sunadaptcontroller/mrihtol/fmod_int32/CMakeLists.txt +++ b/src/sunadaptcontroller/mrihtol/fmod_int32/CMakeLists.txt @@ -16,7 +16,6 @@ sundials_add_f2003_library( sundials_fsunadaptcontrollermrihtol_mod SOURCES fsunadaptcontroller_mrihtol_mod.f90 fsunadaptcontroller_mrihtol_mod.c LINK_LIBRARIES PUBLIC sundials_fcore_mod - OBJECT_LIBRARIES OUTPUT_NAME sundials_fsunadaptcontrollermrihtol_mod OBJECT_LIB_ONLY) message(STATUS "Added SUNAdaptController_MRIHTol F2003 interface") diff --git a/src/sunadaptcontroller/mrihtol/fmod_int32/fsunadaptcontroller_mrihtol_mod.c b/src/sunadaptcontroller/mrihtol/fmod_int32/fsunadaptcontroller_mrihtol_mod.c index dcb3ecddda..e8a0775fe4 100644 --- a/src/sunadaptcontroller/mrihtol/fmod_int32/fsunadaptcontroller_mrihtol_mod.c +++ b/src/sunadaptcontroller/mrihtol/fmod_int32/fsunadaptcontroller_mrihtol_mod.c @@ -596,7 +596,7 @@ SWIGEXPORT int _wrap_FSUNAdaptController_SetErrorBias_MRIHTol(SUNAdaptController } -SWIGEXPORT int _wrap_FSUNAdaptController_UpdateMRITol_MRIHTol(SUNAdaptController farg1, double const *farg2, double const *farg3, double const *farg4, double const *farg5) { +SWIGEXPORT int _wrap_FSUNAdaptController_UpdateMRIHTol_MRIHTol(SUNAdaptController farg1, double const *farg2, double const *farg3, double const *farg4, double const *farg5) { int fresult ; SUNAdaptController arg1 = (SUNAdaptController) 0 ; sunrealtype arg2 ; @@ -610,7 +610,7 @@ SWIGEXPORT int _wrap_FSUNAdaptController_UpdateMRITol_MRIHTol(SUNAdaptController arg3 = (sunrealtype)(*farg3); arg4 = (sunrealtype)(*farg4); arg5 = (sunrealtype)(*farg5); - result = (int)SUNAdaptController_UpdateMRITol_MRIHTol(arg1,arg2,arg3,arg4,arg5); + result = (int)SUNAdaptController_UpdateMRIHTol_MRIHTol(arg1,arg2,arg3,arg4,arg5); fresult = (int)(result); return fresult; } diff --git a/src/sunadaptcontroller/mrihtol/fmod_int32/fsunadaptcontroller_mrihtol_mod.f90 b/src/sunadaptcontroller/mrihtol/fmod_int32/fsunadaptcontroller_mrihtol_mod.f90 index c08850f180..d333a6dd9c 100644 --- a/src/sunadaptcontroller/mrihtol/fmod_int32/fsunadaptcontroller_mrihtol_mod.f90 +++ b/src/sunadaptcontroller/mrihtol/fmod_int32/fsunadaptcontroller_mrihtol_mod.f90 @@ -64,7 +64,7 @@ module fsunadaptcontroller_mrihtol_mod public :: FSUNAdaptController_SetDefaults_MRIHTol public :: FSUNAdaptController_Write_MRIHTol public :: FSUNAdaptController_SetErrorBias_MRIHTol - public :: FSUNAdaptController_UpdateMRITol_MRIHTol + public :: FSUNAdaptController_UpdateMRIHTol_MRIHTol public :: FSUNAdaptController_Space_MRIHTol ! WRAPPER DECLARATIONS @@ -271,8 +271,8 @@ function swigc_FSUNAdaptController_SetErrorBias_MRIHTol(farg1, farg2) & integer(C_INT) :: fresult end function -function swigc_FSUNAdaptController_UpdateMRITol_MRIHTol(farg1, farg2, farg3, farg4, farg5) & -bind(C, name="_wrap_FSUNAdaptController_UpdateMRITol_MRIHTol") & +function swigc_FSUNAdaptController_UpdateMRIHTol_MRIHTol(farg1, farg2, farg3, farg4, farg5) & +bind(C, name="_wrap_FSUNAdaptController_UpdateMRIHTol_MRIHTol") & result(fresult) use, intrinsic :: ISO_C_BINDING type(C_PTR), value :: farg1 @@ -632,7 +632,7 @@ function FSUNAdaptController_SetErrorBias_MRIHTol(c, bias) & swig_result = fresult end function -function FSUNAdaptController_UpdateMRITol_MRIHTol(c, h, tolfac, dsm, dsm4) & +function FSUNAdaptController_UpdateMRIHTol_MRIHTol(c, h, tolfac, dsm, dsm4) & result(swig_result) use, intrinsic :: ISO_C_BINDING integer(C_INT) :: swig_result @@ -653,7 +653,7 @@ function FSUNAdaptController_UpdateMRITol_MRIHTol(c, h, tolfac, dsm, dsm4) & farg3 = tolfac farg4 = dsm farg5 = dsm4 -fresult = swigc_FSUNAdaptController_UpdateMRITol_MRIHTol(farg1, farg2, farg3, farg4, farg5) +fresult = swigc_FSUNAdaptController_UpdateMRIHTol_MRIHTol(farg1, farg2, farg3, farg4, farg5) swig_result = fresult end function diff --git a/src/sunadaptcontroller/mrihtol/fmod_int64/CMakeLists.txt b/src/sunadaptcontroller/mrihtol/fmod_int64/CMakeLists.txt index ecd215bc10..befc6cb2ed 100644 --- a/src/sunadaptcontroller/mrihtol/fmod_int64/CMakeLists.txt +++ b/src/sunadaptcontroller/mrihtol/fmod_int64/CMakeLists.txt @@ -16,7 +16,6 @@ sundials_add_f2003_library( sundials_fsunadaptcontrollermrihtol_mod SOURCES fsunadaptcontroller_mrihtol_mod.f90 fsunadaptcontroller_mrihtol_mod.c LINK_LIBRARIES PUBLIC sundials_fcore_mod - OBJECT_LIBRARIES OUTPUT_NAME sundials_fsunadaptcontrollermrihtol_mod OBJECT_LIB_ONLY) message(STATUS "Added SUNAdaptController_MRIHTol F2003 interface") diff --git a/src/sunadaptcontroller/mrihtol/fmod_int64/fsunadaptcontroller_mrihtol_mod.c b/src/sunadaptcontroller/mrihtol/fmod_int64/fsunadaptcontroller_mrihtol_mod.c index dcb3ecddda..e8a0775fe4 100644 --- a/src/sunadaptcontroller/mrihtol/fmod_int64/fsunadaptcontroller_mrihtol_mod.c +++ b/src/sunadaptcontroller/mrihtol/fmod_int64/fsunadaptcontroller_mrihtol_mod.c @@ -596,7 +596,7 @@ SWIGEXPORT int _wrap_FSUNAdaptController_SetErrorBias_MRIHTol(SUNAdaptController } -SWIGEXPORT int _wrap_FSUNAdaptController_UpdateMRITol_MRIHTol(SUNAdaptController farg1, double const *farg2, double const *farg3, double const *farg4, double const *farg5) { +SWIGEXPORT int _wrap_FSUNAdaptController_UpdateMRIHTol_MRIHTol(SUNAdaptController farg1, double const *farg2, double const *farg3, double const *farg4, double const *farg5) { int fresult ; SUNAdaptController arg1 = (SUNAdaptController) 0 ; sunrealtype arg2 ; @@ -610,7 +610,7 @@ SWIGEXPORT int _wrap_FSUNAdaptController_UpdateMRITol_MRIHTol(SUNAdaptController arg3 = (sunrealtype)(*farg3); arg4 = (sunrealtype)(*farg4); arg5 = (sunrealtype)(*farg5); - result = (int)SUNAdaptController_UpdateMRITol_MRIHTol(arg1,arg2,arg3,arg4,arg5); + result = (int)SUNAdaptController_UpdateMRIHTol_MRIHTol(arg1,arg2,arg3,arg4,arg5); fresult = (int)(result); return fresult; } diff --git a/src/sunadaptcontroller/mrihtol/fmod_int64/fsunadaptcontroller_mrihtol_mod.f90 b/src/sunadaptcontroller/mrihtol/fmod_int64/fsunadaptcontroller_mrihtol_mod.f90 index c08850f180..d333a6dd9c 100644 --- a/src/sunadaptcontroller/mrihtol/fmod_int64/fsunadaptcontroller_mrihtol_mod.f90 +++ b/src/sunadaptcontroller/mrihtol/fmod_int64/fsunadaptcontroller_mrihtol_mod.f90 @@ -64,7 +64,7 @@ module fsunadaptcontroller_mrihtol_mod public :: FSUNAdaptController_SetDefaults_MRIHTol public :: FSUNAdaptController_Write_MRIHTol public :: FSUNAdaptController_SetErrorBias_MRIHTol - public :: FSUNAdaptController_UpdateMRITol_MRIHTol + public :: FSUNAdaptController_UpdateMRIHTol_MRIHTol public :: FSUNAdaptController_Space_MRIHTol ! WRAPPER DECLARATIONS @@ -271,8 +271,8 @@ function swigc_FSUNAdaptController_SetErrorBias_MRIHTol(farg1, farg2) & integer(C_INT) :: fresult end function -function swigc_FSUNAdaptController_UpdateMRITol_MRIHTol(farg1, farg2, farg3, farg4, farg5) & -bind(C, name="_wrap_FSUNAdaptController_UpdateMRITol_MRIHTol") & +function swigc_FSUNAdaptController_UpdateMRIHTol_MRIHTol(farg1, farg2, farg3, farg4, farg5) & +bind(C, name="_wrap_FSUNAdaptController_UpdateMRIHTol_MRIHTol") & result(fresult) use, intrinsic :: ISO_C_BINDING type(C_PTR), value :: farg1 @@ -632,7 +632,7 @@ function FSUNAdaptController_SetErrorBias_MRIHTol(c, bias) & swig_result = fresult end function -function FSUNAdaptController_UpdateMRITol_MRIHTol(c, h, tolfac, dsm, dsm4) & +function FSUNAdaptController_UpdateMRIHTol_MRIHTol(c, h, tolfac, dsm, dsm4) & result(swig_result) use, intrinsic :: ISO_C_BINDING integer(C_INT) :: swig_result @@ -653,7 +653,7 @@ function FSUNAdaptController_UpdateMRITol_MRIHTol(c, h, tolfac, dsm, dsm4) & farg3 = tolfac farg4 = dsm farg5 = dsm4 -fresult = swigc_FSUNAdaptController_UpdateMRITol_MRIHTol(farg1, farg2, farg3, farg4, farg5) +fresult = swigc_FSUNAdaptController_UpdateMRIHTol_MRIHTol(farg1, farg2, farg3, farg4, farg5) swig_result = fresult end function diff --git a/src/sunadaptcontroller/mrihtol/sunadaptcontroller_mrihtol.c b/src/sunadaptcontroller/mrihtol/sunadaptcontroller_mrihtol.c index ae8716258a..202779f7ef 100644 --- a/src/sunadaptcontroller/mrihtol/sunadaptcontroller_mrihtol.c +++ b/src/sunadaptcontroller/mrihtol/sunadaptcontroller_mrihtol.c @@ -82,7 +82,7 @@ SUNAdaptController SUNAdaptController_MRIHTol(SUNAdaptController HControl, C->ops->setdefaults = SUNAdaptController_SetDefaults_MRIHTol; C->ops->write = SUNAdaptController_Write_MRIHTol; C->ops->seterrorbias = SUNAdaptController_SetErrorBias_MRIHTol; - C->ops->updatemritol = SUNAdaptController_UpdateMRITol_MRIHTol; + C->ops->updatemrihtol = SUNAdaptController_UpdateMRIHTol_MRIHTol; C->ops->space = SUNAdaptController_Space_MRIHTol; /* Create content */ @@ -102,7 +102,7 @@ SUNAdaptController SUNAdaptController_MRIHTol(SUNAdaptController HControl, /* Attach content */ C->content = content; - return (C); + return C; } /* ----------------------------------------------------------------- @@ -144,7 +144,7 @@ SUNAdaptController SUNAdaptController_GetFastController_MRIHTol(SUNAdaptControll SUNAdaptController_Type SUNAdaptController_GetType_MRIHTol( SUNDIALS_MAYBE_UNUSED SUNAdaptController C) { - return SUN_ADAPTCONTROLLER_MRI_TOL; + return SUN_ADAPTCONTROLLER_MRI_H_TOL; } SUNErrCode SUNAdaptController_EstimateStepTol_MRIHTol( @@ -227,11 +227,11 @@ SUNErrCode SUNAdaptController_SetErrorBias_MRIHTol(SUNAdaptController C, return SUN_SUCCESS; } -SUNErrCode SUNAdaptController_UpdateMRITol_MRIHTol(SUNAdaptController C, - sunrealtype H, - sunrealtype tolfac, - sunrealtype DSM, - sunrealtype dsm) +SUNErrCode SUNAdaptController_UpdateMRIHTol_MRIHTol(SUNAdaptController C, + sunrealtype H, + sunrealtype tolfac, + sunrealtype DSM, + sunrealtype dsm) { SUNFunctionBegin(C->sunctx); SUNCheckCall(SUNAdaptController_UpdateH(MRIHTOL_CSLOW(C), H, DSM)); diff --git a/src/sundials/fmod_int32/fsundials_core_mod.c b/src/sundials/fmod_int32/fsundials_core_mod.c index 8590361c6a..d24e5f2f12 100644 --- a/src/sundials/fmod_int32/fsundials_core_mod.c +++ b/src/sundials/fmod_int32/fsundials_core_mod.c @@ -2651,7 +2651,7 @@ SWIGEXPORT int _wrap_FSUNAdaptController_UpdateH(SUNAdaptController farg1, doubl } -SWIGEXPORT int _wrap_FSUNAdaptController_UpdateMRITol(SUNAdaptController farg1, double const *farg2, double const *farg3, double const *farg4, double const *farg5) { +SWIGEXPORT int _wrap_FSUNAdaptController_UpdateMRIHTol(SUNAdaptController farg1, double const *farg2, double const *farg3, double const *farg4, double const *farg5) { int fresult ; SUNAdaptController arg1 = (SUNAdaptController) 0 ; sunrealtype arg2 ; @@ -2665,7 +2665,7 @@ SWIGEXPORT int _wrap_FSUNAdaptController_UpdateMRITol(SUNAdaptController farg1, arg3 = (sunrealtype)(*farg3); arg4 = (sunrealtype)(*farg4); arg5 = (sunrealtype)(*farg5); - result = (SUNErrCode)SUNAdaptController_UpdateMRITol(arg1,arg2,arg3,arg4,arg5); + result = (SUNErrCode)SUNAdaptController_UpdateMRIHTol(arg1,arg2,arg3,arg4,arg5); fresult = (SUNErrCode)(result); return fresult; } diff --git a/src/sundials/fmod_int32/fsundials_core_mod.f90 b/src/sundials/fmod_int32/fsundials_core_mod.f90 index ff9ffaa91c..60b6da39ce 100644 --- a/src/sundials/fmod_int32/fsundials_core_mod.f90 +++ b/src/sundials/fmod_int32/fsundials_core_mod.f90 @@ -511,10 +511,10 @@ module fsundials_core_mod enum, bind(c) enumerator :: SUN_ADAPTCONTROLLER_NONE enumerator :: SUN_ADAPTCONTROLLER_H - enumerator :: SUN_ADAPTCONTROLLER_MRI_TOL + enumerator :: SUN_ADAPTCONTROLLER_MRI_H_TOL end enum integer, parameter, public :: SUNAdaptController_Type = kind(SUN_ADAPTCONTROLLER_NONE) - public :: SUN_ADAPTCONTROLLER_NONE, SUN_ADAPTCONTROLLER_H, SUN_ADAPTCONTROLLER_MRI_TOL + public :: SUN_ADAPTCONTROLLER_NONE, SUN_ADAPTCONTROLLER_H, SUN_ADAPTCONTROLLER_MRI_H_TOL ! struct struct _generic_SUNAdaptController_Ops type, bind(C), public :: SUNAdaptController_Ops type(C_FUNPTR), public :: gettype @@ -526,7 +526,7 @@ module fsundials_core_mod type(C_FUNPTR), public :: write type(C_FUNPTR), public :: seterrorbias type(C_FUNPTR), public :: updateh - type(C_FUNPTR), public :: updatemritol + type(C_FUNPTR), public :: updatemrihtol type(C_FUNPTR), public :: space end type SUNAdaptController_Ops ! struct struct _generic_SUNAdaptController @@ -546,7 +546,7 @@ module fsundials_core_mod public :: FSUNAdaptController_Write public :: FSUNAdaptController_SetErrorBias public :: FSUNAdaptController_UpdateH - public :: FSUNAdaptController_UpdateMRITol + public :: FSUNAdaptController_UpdateMRIHTol public :: FSUNAdaptController_Space ! WRAPPER DECLARATIONS @@ -2051,8 +2051,8 @@ function swigc_FSUNAdaptController_UpdateH(farg1, farg2, farg3) & integer(C_INT) :: fresult end function -function swigc_FSUNAdaptController_UpdateMRITol(farg1, farg2, farg3, farg4, farg5) & -bind(C, name="_wrap_FSUNAdaptController_UpdateMRITol") & +function swigc_FSUNAdaptController_UpdateMRIHTol(farg1, farg2, farg3, farg4, farg5) & +bind(C, name="_wrap_FSUNAdaptController_UpdateMRIHTol") & result(fresult) use, intrinsic :: ISO_C_BINDING type(C_PTR), value :: farg1 @@ -4829,7 +4829,7 @@ function FSUNAdaptController_UpdateH(c, h, dsm) & swig_result = fresult end function -function FSUNAdaptController_UpdateMRITol(c, h, tolfac, dsm, dsm4) & +function FSUNAdaptController_UpdateMRIHTol(c, h, tolfac, dsm, dsm4) & result(swig_result) use, intrinsic :: ISO_C_BINDING integer(C_INT) :: swig_result @@ -4850,7 +4850,7 @@ function FSUNAdaptController_UpdateMRITol(c, h, tolfac, dsm, dsm4) & farg3 = tolfac farg4 = dsm farg5 = dsm4 -fresult = swigc_FSUNAdaptController_UpdateMRITol(farg1, farg2, farg3, farg4, farg5) +fresult = swigc_FSUNAdaptController_UpdateMRIHTol(farg1, farg2, farg3, farg4, farg5) swig_result = fresult end function diff --git a/src/sundials/fmod_int64/fsundials_core_mod.c b/src/sundials/fmod_int64/fsundials_core_mod.c index 98cb592e42..ac00ce7477 100644 --- a/src/sundials/fmod_int64/fsundials_core_mod.c +++ b/src/sundials/fmod_int64/fsundials_core_mod.c @@ -2651,7 +2651,7 @@ SWIGEXPORT int _wrap_FSUNAdaptController_UpdateH(SUNAdaptController farg1, doubl } -SWIGEXPORT int _wrap_FSUNAdaptController_UpdateMRITol(SUNAdaptController farg1, double const *farg2, double const *farg3, double const *farg4, double const *farg5) { +SWIGEXPORT int _wrap_FSUNAdaptController_UpdateMRIHTol(SUNAdaptController farg1, double const *farg2, double const *farg3, double const *farg4, double const *farg5) { int fresult ; SUNAdaptController arg1 = (SUNAdaptController) 0 ; sunrealtype arg2 ; @@ -2665,7 +2665,7 @@ SWIGEXPORT int _wrap_FSUNAdaptController_UpdateMRITol(SUNAdaptController farg1, arg3 = (sunrealtype)(*farg3); arg4 = (sunrealtype)(*farg4); arg5 = (sunrealtype)(*farg5); - result = (SUNErrCode)SUNAdaptController_UpdateMRITol(arg1,arg2,arg3,arg4,arg5); + result = (SUNErrCode)SUNAdaptController_UpdateMRIHTol(arg1,arg2,arg3,arg4,arg5); fresult = (SUNErrCode)(result); return fresult; } diff --git a/src/sundials/fmod_int64/fsundials_core_mod.f90 b/src/sundials/fmod_int64/fsundials_core_mod.f90 index 117f0567a3..a157807b05 100644 --- a/src/sundials/fmod_int64/fsundials_core_mod.f90 +++ b/src/sundials/fmod_int64/fsundials_core_mod.f90 @@ -511,10 +511,10 @@ module fsundials_core_mod enum, bind(c) enumerator :: SUN_ADAPTCONTROLLER_NONE enumerator :: SUN_ADAPTCONTROLLER_H - enumerator :: SUN_ADAPTCONTROLLER_MRI_TOL + enumerator :: SUN_ADAPTCONTROLLER_MRI_H_TOL end enum integer, parameter, public :: SUNAdaptController_Type = kind(SUN_ADAPTCONTROLLER_NONE) - public :: SUN_ADAPTCONTROLLER_NONE, SUN_ADAPTCONTROLLER_H, SUN_ADAPTCONTROLLER_MRI_TOL + public :: SUN_ADAPTCONTROLLER_NONE, SUN_ADAPTCONTROLLER_H, SUN_ADAPTCONTROLLER_MRI_H_TOL ! struct struct _generic_SUNAdaptController_Ops type, bind(C), public :: SUNAdaptController_Ops type(C_FUNPTR), public :: gettype @@ -526,7 +526,7 @@ module fsundials_core_mod type(C_FUNPTR), public :: write type(C_FUNPTR), public :: seterrorbias type(C_FUNPTR), public :: updateh - type(C_FUNPTR), public :: updatemritol + type(C_FUNPTR), public :: updatemrihtol type(C_FUNPTR), public :: space end type SUNAdaptController_Ops ! struct struct _generic_SUNAdaptController @@ -546,7 +546,7 @@ module fsundials_core_mod public :: FSUNAdaptController_Write public :: FSUNAdaptController_SetErrorBias public :: FSUNAdaptController_UpdateH - public :: FSUNAdaptController_UpdateMRITol + public :: FSUNAdaptController_UpdateMRIHTol public :: FSUNAdaptController_Space ! WRAPPER DECLARATIONS @@ -2051,8 +2051,8 @@ function swigc_FSUNAdaptController_UpdateH(farg1, farg2, farg3) & integer(C_INT) :: fresult end function -function swigc_FSUNAdaptController_UpdateMRITol(farg1, farg2, farg3, farg4, farg5) & -bind(C, name="_wrap_FSUNAdaptController_UpdateMRITol") & +function swigc_FSUNAdaptController_UpdateMRIHTol(farg1, farg2, farg3, farg4, farg5) & +bind(C, name="_wrap_FSUNAdaptController_UpdateMRIHTol") & result(fresult) use, intrinsic :: ISO_C_BINDING type(C_PTR), value :: farg1 @@ -4829,7 +4829,7 @@ function FSUNAdaptController_UpdateH(c, h, dsm) & swig_result = fresult end function -function FSUNAdaptController_UpdateMRITol(c, h, tolfac, dsm, dsm4) & +function FSUNAdaptController_UpdateMRIHTol(c, h, tolfac, dsm, dsm4) & result(swig_result) use, intrinsic :: ISO_C_BINDING integer(C_INT) :: swig_result @@ -4850,7 +4850,7 @@ function FSUNAdaptController_UpdateMRITol(c, h, tolfac, dsm, dsm4) & farg3 = tolfac farg4 = dsm farg5 = dsm4 -fresult = swigc_FSUNAdaptController_UpdateMRITol(farg1, farg2, farg3, farg4, farg5) +fresult = swigc_FSUNAdaptController_UpdateMRIHTol(farg1, farg2, farg3, farg4, farg5) swig_result = fresult end function diff --git a/src/sundials/sundials_adaptcontroller.c b/src/sundials/sundials_adaptcontroller.c index e841b08a60..4641f1751c 100644 --- a/src/sundials/sundials_adaptcontroller.c +++ b/src/sundials/sundials_adaptcontroller.c @@ -55,7 +55,7 @@ SUNAdaptController SUNAdaptController_NewEmpty(SUNContext sunctx) ops->write = NULL; ops->seterrorbias = NULL; ops->updateh = NULL; - ops->updatemritol = NULL; + ops->updatemrihtol = NULL; ops->space = NULL; /* attach ops and initialize content to NULL */ @@ -206,16 +206,16 @@ SUNErrCode SUNAdaptController_UpdateH(SUNAdaptController C, sunrealtype h, return (ier); } -SUNErrCode SUNAdaptController_UpdateMRITol(SUNAdaptController C, sunrealtype H, - sunrealtype tolfac, sunrealtype DSM, - sunrealtype dsm) +SUNErrCode SUNAdaptController_UpdateMRIHTol(SUNAdaptController C, sunrealtype H, + sunrealtype tolfac, sunrealtype DSM, + sunrealtype dsm) { SUNErrCode ier = SUN_SUCCESS; if (C == NULL) { return SUN_ERR_ARG_CORRUPT; } SUNFunctionBegin(C->sunctx); - if (C->ops->updatemritol) + if (C->ops->updatemrihtol) { - ier = C->ops->updatemritol(C, H, tolfac, DSM, dsm); + ier = C->ops->updatemrihtol(C, H, tolfac, DSM, dsm); } return (ier); } diff --git a/test/answers b/test/answers index c79677ea5e..cbf335f7e2 160000 --- a/test/answers +++ b/test/answers @@ -1 +1 @@ -Subproject commit c79677ea5ecd18b5a53456f7ec1314de1c637d6d +Subproject commit cbf335f7e244de7c12a4dbd638b4e0f02811a4a3 diff --git a/test/unit_tests/arkode/CXX_serial/ark_test_accumerror_brusselator.cpp b/test/unit_tests/arkode/CXX_serial/ark_test_accumerror_brusselator.cpp index 7869048a90..730e7f2a84 100644 --- a/test/unit_tests/arkode/CXX_serial/ark_test_accumerror_brusselator.cpp +++ b/test/unit_tests/arkode/CXX_serial/ark_test_accumerror_brusselator.cpp @@ -152,11 +152,12 @@ int main(int argc, char* argv[]) N_Vector* yref = NULL; // empty vectors for storing reference solution void* arkode_mem = NULL; // empty ARKStep memory structure void* arkode_ref = NULL; // empty ARKStep memory structure for reference solution - SUNMatrix A = NULL; // empty matrix for solver - SUNLinearSolver LS = NULL; // empty linear solver object - UserData udata; // user-data structure - udata.ep = SUN_RCONST(0.0004); // stiffness parameter - udata.Npart = 20; // partition size + SUNMatrix A = NULL; // empty matrix for solver + SUNLinearSolver LS = NULL; // empty linear solver object + UserData udata; // user-data structure + sunrealtype* ydata = NULL; + udata.ep = SUN_RCONST(0.0004); // stiffness parameter + udata.Npart = 20; // partition size // // Initialization @@ -256,12 +257,14 @@ int main(int argc, char* argv[]) if (check_retval((void*)y, "N_VNew_Serial", 0)) return 1; yref = N_VCloneVectorArray(udata.Npart + 1, y); if (check_retval((void*)yref, "N_VNew_Serial", 0)) return 1; + ydata = N_VGetArrayPointer(y); + if (check_retval((void*)ydata, "N_VGetArrayPointer", 0)) return 1; // Generate reference solution - NV_Ith_S(y, 0) = u0; - NV_Ith_S(y, 1) = v0; - NV_Ith_S(y, 2) = w0; - arkode_ref = ARKStepCreate(fn, NULL, T0, y, ctx); + ydata[0] = u0; + ydata[1] = v0; + ydata[2] = w0; + arkode_ref = ARKStepCreate(fn, NULL, T0, y, ctx); if (check_retval((void*)arkode_ref, "ARKStepCreate", 0)) return 1; retval = ARKodeSetUserData(arkode_ref, (void*)&udata); if (check_retval(&retval, "ARKodeSetUserData", 1)) return 1; @@ -285,9 +288,9 @@ int main(int argc, char* argv[]) ARKodeFree(&arkode_ref); // Set up ARKStep integrator - NV_Ith_S(y, 0) = u0; - NV_Ith_S(y, 1) = v0; - NV_Ith_S(y, 2) = w0; + ydata[0] = u0; + ydata[1] = v0; + ydata[2] = w0; if (rk_type == 0) { // DIRK method arkode_mem = ARKStepCreate(NULL, fn, T0, y, ctx); @@ -347,15 +350,17 @@ int main(int argc, char* argv[]) static int fn(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { - UserData* udata = (UserData*)user_data; - sunrealtype u = NV_Ith_S(y, 0); // access solution values - sunrealtype v = NV_Ith_S(y, 1); - sunrealtype w = NV_Ith_S(y, 2); + UserData* udata = (UserData*)user_data; + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; // access solution values + const sunrealtype v = ydata[1]; + const sunrealtype w = ydata[2]; // fill in the RHS function - NV_Ith_S(ydot, 0) = udata->a - (w + ONE) * u + v * u * u; - NV_Ith_S(ydot, 1) = w * u - v * u * u; - NV_Ith_S(ydot, 2) = (udata->b - w) / udata->ep - w * u; + dydata[0] = udata->a - (w + ONE) * u + v * u * u; + dydata[1] = w * u - v * u * u; + dydata[2] = (udata->b - w) / udata->ep - w * u; // Return with success return 0; @@ -364,10 +369,11 @@ static int fn(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int Jac(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { - UserData* udata = (UserData*)user_data; - sunrealtype u = NV_Ith_S(y, 0); // access solution values - sunrealtype v = NV_Ith_S(y, 1); - sunrealtype w = NV_Ith_S(y, 2); + UserData* udata = (UserData*)user_data; + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; // access solution values + const sunrealtype v = ydata[1]; + const sunrealtype w = ydata[2]; // fill in the Jacobian SM_ELEMENT_D(J, 0, 0) = -(w + ONE) + TWO * u * v; @@ -451,15 +457,14 @@ static int adaptive_run(void* arkode_mem, N_Vector y, sunrealtype T0, if (check_retval(&retval, "ARKodeGetNumSteps", 1)) break; // Compute/print solution error - sunrealtype uref = NV_Ith_S(yref[ipart + 1], 0); - sunrealtype vref = NV_Ith_S(yref[ipart + 1], 1); - sunrealtype wref = NV_Ith_S(yref[ipart + 1], 2); - sunrealtype udsm = abs(NV_Ith_S(y, 0) - uref) / - (abstol + rtols[irtol] * abs(uref)); - sunrealtype vdsm = abs(NV_Ith_S(y, 1) - vref) / - (abstol + rtols[irtol] * abs(vref)); - sunrealtype wdsm = abs(NV_Ith_S(y, 2) - wref) / - (abstol + rtols[irtol] * abs(wref)); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* yrefdata = N_VGetArrayPointer(yref[ipart + 1]); + const sunrealtype udsm = abs(ydata[0] - yrefdata[0]) / + (abstol + rtols[irtol] * abs(yrefdata[0])); + const sunrealtype vdsm = abs(ydata[1] - yrefdata[1]) / + (abstol + rtols[irtol] * abs(yrefdata[1])); + const sunrealtype wdsm = abs(ydata[2] - yrefdata[2]) / + (abstol + rtols[irtol] * abs(yrefdata[2])); dsm[ipart] = rtols[irtol] * sqrt((udsm * udsm + vdsm * vdsm + wdsm * wdsm) / SUN_RCONST(3.0)); @@ -550,12 +555,14 @@ static int fixed_run(void* arkode_mem, N_Vector y, sunrealtype T0, sunrealtype T if (check_retval(&retval, "ARKodeGetNumSteps", 1)) break; // Compute/print solution error - sunrealtype udsm = abs(NV_Ith_S(y, 0) - NV_Ith_S(yref[ipart + 1], 0)) / - (abstol + reltol * abs(NV_Ith_S(yref[ipart + 1], 0))); - sunrealtype vdsm = abs(NV_Ith_S(y, 1) - NV_Ith_S(yref[ipart + 1], 1)) / - (abstol + reltol * abs(NV_Ith_S(yref[ipart + 1], 1))); - sunrealtype wdsm = abs(NV_Ith_S(y, 2) - NV_Ith_S(yref[ipart + 1], 2)) / - (abstol + reltol * abs(NV_Ith_S(yref[ipart + 1], 2))); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* yrefdata = N_VGetArrayPointer(yref[ipart + 1]); + const sunrealtype udsm = abs(ydata[0] - yrefdata[0]) / + (abstol + reltol * abs(yrefdata[0])); + const sunrealtype vdsm = abs(ydata[1] - yrefdata[1]) / + (abstol + reltol * abs(yrefdata[1])); + const sunrealtype wdsm = abs(ydata[2] - yrefdata[2]) / + (abstol + reltol * abs(yrefdata[2])); dsm[ipart] = reltol * sqrt((udsm * udsm + vdsm * vdsm + wdsm * wdsm) / SUN_RCONST(3.0)); cout << " h " << hvals[ih] << " rk_type " << rk_type << " order " @@ -631,12 +638,14 @@ static int fixed_run(void* arkode_mem, N_Vector y, sunrealtype T0, sunrealtype T dsm_est[ipart] = reltol * N_VWrmsNorm(y2, ewt); Nsteps[ipart] += nsteps2; - sunrealtype udsm = abs(NV_Ith_S(y, 0) - NV_Ith_S(yref[ipart + 1], 0)) / - (abstol + reltol * abs(NV_Ith_S(yref[ipart + 1], 0))); - sunrealtype vdsm = abs(NV_Ith_S(y, 1) - NV_Ith_S(yref[ipart + 1], 1)) / - (abstol + reltol * abs(NV_Ith_S(yref[ipart + 1], 1))); - sunrealtype wdsm = abs(NV_Ith_S(y, 2) - NV_Ith_S(yref[ipart + 1], 2)) / - (abstol + reltol * abs(NV_Ith_S(yref[ipart + 1], 2))); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* yrefdata = N_VGetArrayPointer(yref[ipart + 1]); + const sunrealtype udsm = abs(ydata[0] - yrefdata[0]) / + (abstol + reltol * abs(yrefdata[0])); + const sunrealtype vdsm = abs(ydata[1] - yrefdata[1]) / + (abstol + reltol * abs(yrefdata[1])); + const sunrealtype wdsm = abs(ydata[2] - yrefdata[2]) / + (abstol + reltol * abs(yrefdata[2])); dsm[ipart] = reltol * sqrt((udsm * udsm + vdsm * vdsm + wdsm * wdsm) / SUN_RCONST(3.0)); cout << " h " << hvals[ih] << " rk_type " << rk_type << " order " diff --git a/test/unit_tests/arkode/CXX_serial/ark_test_accumerror_kpr.cpp b/test/unit_tests/arkode/CXX_serial/ark_test_accumerror_kpr.cpp index b5bf8f48ae..3707bff97c 100644 --- a/test/unit_tests/arkode/CXX_serial/ark_test_accumerror_kpr.cpp +++ b/test/unit_tests/arkode/CXX_serial/ark_test_accumerror_kpr.cpp @@ -280,11 +280,7 @@ int main(int argc, char* argv[]) } // Clean up and return - if (rk_type == 0) - { // Free integrator memory - ARKodeFree(&arkode_mem); - } - else { ARKodeFree(&arkode_mem); } + ARKodeFree(&arkode_mem); if (LS != NULL) SUNLinSolFree(LS); // free system linear solver if (A != NULL) SUNMatDestroy(A); // free system matrix N_VDestroy(y); // Free y vector @@ -298,17 +294,19 @@ int main(int argc, char* argv[]) static int fn(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { UserData* udata = (UserData*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; sunrealtype tmp1, tmp2; // fill in the RHS function: // [G e]*[(-2+u^2-p(t))/(2*u)] + [pdot(t)/(2u)] // [e -1] [(-2+v^2-s(t))/(2*v)] [qdot(t)/(2v)] - tmp1 = (-TWO + u * u - p(t)) / (TWO * u); - tmp2 = (-TWO + v * v - q(t, *udata)) / (TWO * v); - NV_Ith_S(ydot, 0) = udata->G * tmp1 + udata->e * tmp2 + pdot(t) / (TWO * u); - NV_Ith_S(ydot, 1) = udata->e * tmp1 - tmp2 + qdot(t, *udata) / (TWO * v); + tmp1 = (-TWO + u * u - p(t)) / (TWO * u); + tmp2 = (-TWO + v * v - q(t, *udata)) / (TWO * v); + dydata[0] = udata->G * tmp1 + udata->e * tmp2 + pdot(t) / (TWO * u); + dydata[1] = udata->e * tmp1 - tmp2 + qdot(t, *udata) / (TWO * v); // Return with success return 0; @@ -318,8 +316,9 @@ static int Jn(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { UserData* udata = (UserData*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; sunrealtype t11, t22; // fill in the Jacobian: @@ -355,6 +354,7 @@ static int adaptive_run(void* arkode_mem, N_Vector y, sunrealtype T0, vector dsm(udata.Npart); vector dsm_est(udata.Npart); vector Nsteps(udata.Npart); + sunrealtype* ydata = N_VGetArrayPointer(y); // Loop over tolerances cout << "\nAdaptive-step runs:\n"; @@ -416,9 +416,9 @@ static int adaptive_run(void* arkode_mem, N_Vector y, sunrealtype T0, } // Compute/print solution error - sunrealtype udsm = abs(NV_Ith_S(y, 0) - utrue(t)) / + sunrealtype udsm = abs(ydata[0] - utrue(t)) / (abstol + rtols[irtol] * abs(utrue(t))); - sunrealtype vdsm = abs(NV_Ith_S(y, 1) - vtrue(t, udata)) / + sunrealtype vdsm = abs(ydata[1] - vtrue(t, udata)) / (abstol + rtols[irtol] * abs(vtrue(t, udata))); dsm[ipart] = rtols[irtol] * sqrt(0.5 * (udsm * udsm + vdsm * vdsm)); cout << " rtol " << rtols[irtol] << " rk_type " << rk_type @@ -455,6 +455,7 @@ static int fixed_run(void* arkode_mem, N_Vector y, sunrealtype T0, vector dsm(udata.Npart); vector dsm_est(udata.Npart); vector Nsteps(udata.Npart); + sunrealtype* ydata = N_VGetArrayPointer(y); // Loop over step sizes cout << "\nFixed-step runs:\n"; @@ -526,9 +527,9 @@ static int fixed_run(void* arkode_mem, N_Vector y, sunrealtype T0, } // Compute/print solution error - sunrealtype udsm = abs(NV_Ith_S(y, 0) - utrue(t)) / + sunrealtype udsm = abs(ydata[0] - utrue(t)) / (abstol + reltol * abs(utrue(t))); - sunrealtype vdsm = abs(NV_Ith_S(y, 1) - vtrue(t, udata)) / + sunrealtype vdsm = abs(ydata[1] - vtrue(t, udata)) / (abstol + reltol * abs(vtrue(t, udata))); dsm[ipart] = reltol * sqrt(0.5 * (udsm * udsm + vdsm * vdsm)); cout << " h " << hvals[ih] << " rk_type " << rk_type << " order " @@ -620,9 +621,9 @@ static int fixed_run(void* arkode_mem, N_Vector y, sunrealtype T0, N_VLinearSum(ONE, y2, -ONE, y, y2); dsm_est[ipart] = reltol * N_VWrmsNorm(y2, ewt); Nsteps[ipart] += nsteps2; - sunrealtype udsm = abs(NV_Ith_S(y, 0) - utrue(t)) / + sunrealtype udsm = abs(ydata[0] - utrue(t)) / (abstol + reltol * abs(utrue(t))); - sunrealtype vdsm = abs(NV_Ith_S(y, 1) - vtrue(t, udata)) / + sunrealtype vdsm = abs(ydata[1] - vtrue(t, udata)) / (abstol + reltol * abs(vtrue(t, udata))); dsm[ipart] = reltol * sqrt(0.5 * (udsm * udsm + vdsm * vdsm)); cout << " h " << hvals[ih] << " rk_type " << rk_type << " order " @@ -662,8 +663,9 @@ static sunrealtype vtrue(sunrealtype t, UserData& udata) static int Ytrue(sunrealtype t, N_Vector y, UserData& udata) { - NV_Ith_S(y, 0) = utrue(t); - NV_Ith_S(y, 1) = vtrue(t, udata); + sunrealtype* ydata = N_VGetArrayPointer(y); + ydata[0] = utrue(t); + ydata[1] = vtrue(t, udata); return (0); } diff --git a/test/unit_tests/arkode/CXX_serial/ark_test_brusselator_mriadapt.cpp b/test/unit_tests/arkode/CXX_serial/ark_test_brusselator_mriadapt.cpp index 57714351d6..63e2b4b39b 100644 --- a/test/unit_tests/arkode/CXX_serial/ark_test_brusselator_mriadapt.cpp +++ b/test/unit_tests/arkode/CXX_serial/ark_test_brusselator_mriadapt.cpp @@ -265,9 +265,11 @@ int main(int argc, char* argv[]) if (check_ptr((void*)yref, "N_VNew_Serial")) return 1; // Set initial conditions - NV_Ith_S(y, 0) = SUN_RCONST(1.2); - NV_Ith_S(y, 1) = SUN_RCONST(3.1); - NV_Ith_S(y, 2) = SUN_RCONST(3.0); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* yrefdata = N_VGetArrayPointer(yref); + ydata[0] = SUN_RCONST(1.2); + ydata[1] = SUN_RCONST(3.1); + ydata[2] = SUN_RCONST(3.0); N_VScale(ONE, y, yref); // Create and configure reference solver object @@ -575,7 +577,7 @@ int main(int argc, char* argv[]) scontrol_Tol = SUNAdaptController_ImExGus(sunctx); if (check_ptr((void*)scontrol_Tol, "SUNAdaptController_ImExGus (slow Tol)")) return 1; - scontrol = SUNAdaptController_MRIHTol(scontrol_H, scontrol_Tolsunctx); + scontrol = SUNAdaptController_MRIHTol(scontrol_H, scontrol_Tol, sunctx); if (check_ptr((void*)scontrol, "SUNAdaptController_MRIHTol")) return 1; if (std::min(opts.htol_relch, std::min(opts.htol_minfac, opts.htol_maxfac)) > -1) @@ -672,17 +674,14 @@ int main(int argc, char* argv[]) "----\n"); printf(" %10.6" FSYM " %10.6" FSYM " %10.6" FSYM " %10.6" FSYM " %.1" ESYM " %.1" ESYM " %.1" ESYM "\n", - t, NV_Ith_S(y, 0), NV_Ith_S(y, 1), NV_Ith_S(y, 2), uerr, verr, werr); - int Nout = 0; - while (Tf - t > 1.0e-8) + t, ydata[0], ydata[1], ydata[2], uerr, verr, werr); + while (Tf - t > SUN_RCONST(1.0e-8)) { // reset reference solver so that it begins with identical state retval = ARKodeReset(arkode_ref, t, y); // evolve solution in one-step mode - retval = ARKodeSetStopTime(arkode_mem, tout); - if (check_flag(retval, "ARKodeSetStopTime")) return 1; - retval = ARKodeEvolve(arkode_mem, tout, y, &t, ARK_ONE_STEP); + retval = ARKodeEvolve(arkode_mem, Tf, y, &t, ARK_ONE_STEP); if (retval < 0) { printf("ARKodeEvolve error (%i)\n", retval); @@ -692,7 +691,7 @@ int main(int argc, char* argv[]) // evolve reference solver to same time in "normal" mode retval = ARKodeSetStopTime(arkode_ref, t); if (check_flag(retval, "ARKodeSetStopTime")) return 1; - retval = ARKodeEvolve(arkode_ref, t, yref, &t2, ARK_NORMAL); + retval = ARKodeEvolve(arkode_ref, Tf, yref, &t2, ARK_NORMAL); if (retval < 0) { printf("ARKodeEvolve reference solution error (%i)\n", retval); @@ -700,23 +699,22 @@ int main(int argc, char* argv[]) } // access/print solution and error - u = NV_Ith_S(y, 0); - v = NV_Ith_S(y, 1); - w = NV_Ith_S(y, 2); - uerr = SUNRabs(NV_Ith_S(yref, 0) - u); - verr = SUNRabs(NV_Ith_S(yref, 1) - v); - werr = SUNRabs(NV_Ith_S(yref, 2) - w); + u = ydata[0]; + v = ydata[1]; + w = ydata[2]; + uerr = std::abs(yrefdata[0] - u); + verr = std::abs(yrefdata[1] - v); + werr = std::abs(yrefdata[2] - w); uerrtot += uerr * uerr; verrtot += verr * verr; werrtot += werr * werr; errtot += uerr * uerr + verr * verr + werr * werr; - accuracy = std::max(accuracy, uerr / SUNRabs(opts.atol + - opts.rtol * NV_Ith_S(yref, 0))); - accuracy = std::max(accuracy, verr / SUNRabs(opts.atol + - opts.rtol * NV_Ith_S(yref, 1))); - accuracy = std::max(accuracy, werr / SUNRabs(opts.atol + - opts.rtol * NV_Ith_S(yref, 2))); - Nout++; + accuracy = std::max(accuracy, + uerr / std::abs(opts.atol + opts.rtol * yrefdata[0])); + accuracy = std::max(accuracy, + verr / std::abs(opts.atol + opts.rtol * yrefdata[1])); + accuracy = std::max(accuracy, + werr / std::abs(opts.atol + opts.rtol * yrefdata[2])); // Periodically output current results to screen if (t >= tout) @@ -728,10 +726,6 @@ int main(int argc, char* argv[]) t, u, v, w, uerr, verr, werr); } } - uerrtot = SUNRsqrt(uerrtot / Nt); - verrtot = SUNRsqrt(verrtot / Nt); - werrtot = SUNRsqrt(werrtot / Nt); - errtot = SUNRsqrt(errtot / Nt / 3); printf(" " "---------------------------------------------------------------------" "----\n"); @@ -748,8 +742,10 @@ int main(int argc, char* argv[]) check_flag(retval, "ARKodeGetNumStepAttempts"); retval = ARKodeGetNumErrTestFails(arkode_mem, &netfs); check_flag(retval, "ARKodeGetNumErrTestFails"); - retval = MRIStepGetNumRhsEvals(arkode_mem, &nfse, &nfsi); - check_flag(retval, "MRIStepGetNumRhsEvals"); + retval = ARKodeGetNumRhsEvals(arkode_mem, 0, &nfse); + check_flag(retval, "ARKodeGetNumRhsEvals"); + retval = ARKodeGetNumRhsEvals(arkode_mem, 1, &nfsi); + check_flag(retval, "ARKodeGetNumRhsEvals"); // Get some fast integrator statistics long int nstf, nattf, netff, nff; @@ -759,10 +755,14 @@ int main(int argc, char* argv[]) check_flag(retval, "ARKodeGetNumStepAttempts"); retval = ARKodeGetNumErrTestFails(inner_arkode_mem, &netff); check_flag(retval, "ARKodeGetNumErrTestFails"); - retval = ERKStepGetNumRhsEvals(inner_arkode_mem, &nff); - check_flag(retval, "ERKStepGetNumRhsEvals"); + retval = ARKodeGetNumRhsEvals(inner_arkode_mem, 0, &nff); + check_flag(retval, "ARKodeGetNumRhsEvals"); // Print some final statistics + uerrtot = std::sqrt(uerrtot / (sunrealtype)nsts); + verrtot = std::sqrt(verrtot / (sunrealtype)nsts); + werrtot = std::sqrt(werrtot / (sunrealtype)nsts); + errtot = std::sqrt(errtot / SUN_RCONST(3.0) / (sunrealtype)nsts); std::cout << "\nFinal Solver Statistics:\n"; std::cout << " Slow steps = " << nsts << " (attempts = " << natts << ", fails = " << netfs << ")\n"; @@ -814,12 +814,14 @@ int main(int argc, char* argv[]) static int ff(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { Options* opts = (Options*)user_data; - const sunrealtype w = NV_Ith_S(y, 2); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype w = ydata[2]; // fill in the RHS function: - NV_Ith_S(ydot, 0) = ZERO; - NV_Ith_S(ydot, 1) = ZERO; - NV_Ith_S(ydot, 2) = (opts->b - w) / opts->ep; + dydata[0] = ZERO; + dydata[1] = ZERO; + dydata[2] = (opts->b - w) / opts->ep; // Return with success return 0; @@ -829,14 +831,16 @@ static int ff(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int fs(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { Options* opts = (Options*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); - const sunrealtype w = NV_Ith_S(y, 2); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; + const sunrealtype w = ydata[2]; // fill in the RHS function: - NV_Ith_S(ydot, 0) = opts->a - (w + ONE) * u + v * u * u; - NV_Ith_S(ydot, 1) = w * u - v * u * u; - NV_Ith_S(ydot, 2) = -w * u; + dydata[0] = opts->a - (w + ONE) * u + v * u * u; + dydata[1] = w * u - v * u * u; + dydata[2] = -w * u; // Return with success return 0; @@ -846,13 +850,15 @@ static int fs(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int fse(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { Options* opts = (Options*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; // fill in the RHS function: - NV_Ith_S(ydot, 0) = opts->a + v * u * u; - NV_Ith_S(ydot, 1) = -v * u * u; - NV_Ith_S(ydot, 2) = ZERO; + dydata[0] = opts->a + v * u * u; + dydata[1] = -v * u * u; + dydata[2] = ZERO; // Return with success return 0; @@ -861,13 +867,15 @@ static int fse(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) // fsi routine to compute the slow portion of the ODE RHS.(currently same as fse) static int fsi(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype w = NV_Ith_S(y, 2); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; + const sunrealtype w = ydata[2]; // fill in the RHS function: - NV_Ith_S(ydot, 0) = -(w + ONE) * u; - NV_Ith_S(ydot, 1) = w * u; - NV_Ith_S(ydot, 2) = -w * u; + dydata[0] = -(w + ONE) * u; + dydata[1] = w * u; + dydata[2] = -w * u; // Return with success return 0; @@ -876,14 +884,16 @@ static int fsi(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int fn(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { Options* opts = (Options*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); - const sunrealtype w = NV_Ith_S(y, 2); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; + const sunrealtype w = ydata[2]; // fill in the RHS function: - NV_Ith_S(ydot, 0) = opts->a - (w + ONE) * u + v * u * u; - NV_Ith_S(ydot, 1) = w * u - v * u * u; - NV_Ith_S(ydot, 2) = (opts->b - w) / opts->ep - w * u; + dydata[0] = opts->a - (w + ONE) * u + v * u * u; + dydata[1] = w * u - v * u * u; + dydata[2] = (opts->b - w) / opts->ep - w * u; // Return with success return 0; @@ -898,9 +908,10 @@ static int f0(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int Js(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); - const sunrealtype w = NV_Ith_S(y, 2); + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; + const sunrealtype w = ydata[2]; // fill in the Jacobian: SM_ELEMENT_D(J, 0, 0) = -(w + ONE) + TWO * u * v; @@ -922,8 +933,9 @@ static int Js(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, static int Jsi(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype w = NV_Ith_S(y, 2); + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; + const sunrealtype w = ydata[2]; // fill in the Jacobian: SM_ELEMENT_D(J, 0, 0) = -(w + ONE); @@ -946,9 +958,10 @@ static int Jn(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { Options* opts = (Options*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); - const sunrealtype w = NV_Ith_S(y, 2); + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; + const sunrealtype w = ydata[2]; // fill in the Jacobian: SM_ELEMENT_D(J, 0, 0) = -(w + ONE) + TWO * u * v; diff --git a/test/unit_tests/arkode/CXX_serial/ark_test_kpr_mriadapt.cpp b/test/unit_tests/arkode/CXX_serial/ark_test_kpr_mriadapt.cpp index 5611780f19..e6eb0ac097 100644 --- a/test/unit_tests/arkode/CXX_serial/ark_test_kpr_mriadapt.cpp +++ b/test/unit_tests/arkode/CXX_serial/ark_test_kpr_mriadapt.cpp @@ -271,6 +271,10 @@ int main(int argc, char* argv[]) if (check_ptr((void*)y, "N_VNew_Serial")) return 1; N_Vector yref = N_VNew_Serial(NEQ, refctx); if (check_ptr((void*)yref, "N_VNew_Serial")) return 1; + sunrealtype* ydata = N_VGetArrayPointer(y); + if (check_ptr((void*)ydata, "N_VGetArrayPointer")) return 1; + sunrealtype* yrefdata = N_VGetArrayPointer(yref); + if (check_ptr((void*)yrefdata, "N_VGetArrayPointer")) return 1; // Set initial conditions retval = Ytrue(T0, y, &opts); @@ -678,17 +682,14 @@ int main(int argc, char* argv[]) printf(" ------------------------------------------------------\n"); printf(" %10.6" FSYM " %10.6" FSYM " %10.6" FSYM " %.2" ESYM " %.2" ESYM "\n", - t, NV_Ith_S(y, 0), NV_Ith_S(y, 1), uerr, verr); - int Nout = 0; - while (Tf - t > 1.0e-8) + t, ydata[0], ydata[1], uerr, verr); + while (Tf - t > SUN_RCONST(1.0e-8)) { // reset reference solver so that it begins with identical state retval = ARKodeReset(arkode_ref, t, y); // evolve solution in one-step mode - retval = ARKodeSetStopTime(arkode_mem, tout); - if (check_flag(retval, "ARKodeSetStopTime")) return 1; - retval = ARKodeEvolve(arkode_mem, tout, y, &t, ARK_ONE_STEP); + retval = ARKodeEvolve(arkode_mem, Tf, y, &t, ARK_ONE_STEP); if (retval < 0) { printf("ARKodeEvolve error (%i)\n", retval); @@ -698,7 +699,7 @@ int main(int argc, char* argv[]) // evolve reference solver to same time in "normal" mode retval = ARKodeSetStopTime(arkode_ref, t); if (check_flag(retval, "ARKodeSetStopTime")) return 1; - retval = ARKodeEvolve(arkode_ref, t, yref, &t2, ARK_NORMAL); + retval = ARKodeEvolve(arkode_ref, Tf, yref, &t2, ARK_NORMAL); if (retval < 0) { printf("ARKodeEvolve reference solution error (%i)\n", retval); @@ -706,18 +707,17 @@ int main(int argc, char* argv[]) } // access/print solution and error - u = NV_Ith_S(y, 0); - v = NV_Ith_S(y, 1); - uerr = SUNRabs(NV_Ith_S(yref, 0) - u); - verr = SUNRabs(NV_Ith_S(yref, 1) - v); + u = ydata[0]; + v = ydata[1]; + uerr = std::abs(yrefdata[0] - u); + verr = std::abs(yrefdata[1] - v); uerrtot += uerr * uerr; verrtot += verr * verr; errtot += uerr * uerr + verr * verr; - accuracy = std::max(accuracy, uerr / SUNRabs(opts.atol + - opts.rtol * NV_Ith_S(yref, 0))); - accuracy = std::max(accuracy, verr / SUNRabs(opts.atol + - opts.rtol * NV_Ith_S(yref, 1))); - Nout++; + accuracy = std::max(accuracy, + uerr / std::abs(opts.atol + opts.rtol * yrefdata[0])); + accuracy = std::max(accuracy, + verr / std::abs(opts.atol + opts.rtol * yrefdata[1])); // Periodically output current results to screen if (t >= tout) @@ -729,9 +729,6 @@ int main(int argc, char* argv[]) t, u, v, uerr, verr); } } - uerrtot = SUNRsqrt(uerrtot / Nt); - verrtot = SUNRsqrt(verrtot / Nt); - errtot = SUNRsqrt(errtot / Nt / 2); printf(" ------------------------------------------------------\n"); // @@ -746,8 +743,10 @@ int main(int argc, char* argv[]) check_flag(retval, "ARKodeGetNumStepAttempts"); retval = ARKodeGetNumErrTestFails(arkode_mem, &netfs); check_flag(retval, "ARKodeGetNumErrTestFails"); - retval = MRIStepGetNumRhsEvals(arkode_mem, &nfse, &nfsi); - check_flag(retval, "MRIStepGetNumRhsEvals"); + retval = ARKodeGetNumRhsEvals(arkode_mem, 0, &nfse); + check_flag(retval, "ARKodeGetNumRhsEvals"); + retval = ARKodeGetNumRhsEvals(arkode_mem, 1, &nfsi); + check_flag(retval, "ARKodeGetNumRhsEvals"); // Get some fast integrator statistics long int nstf, nattf, netff, nff; @@ -757,10 +756,13 @@ int main(int argc, char* argv[]) check_flag(retval, "ARKodeGetNumStepAttempts"); retval = ARKodeGetNumErrTestFails(inner_arkode_mem, &netff); check_flag(retval, "ARKodeGetNumErrTestFails"); - retval = ERKStepGetNumRhsEvals(inner_arkode_mem, &nff); - check_flag(retval, "ERKStepGetNumRhsEvals"); + retval = ARKodeGetNumRhsEvals(inner_arkode_mem, 0, &nff); + check_flag(retval, "ARKodeGetNumRhsEvals"); // Print some final statistics + uerrtot = std::sqrt(uerrtot / (sunrealtype)nsts); + verrtot = std::sqrt(verrtot / (sunrealtype)nsts); + errtot = std::sqrt(errtot / SUN_RCONST(2.0) / (sunrealtype)nsts); std::cout << "\nFinal Solver Statistics:\n"; std::cout << " Slow steps = " << nsts << " (attempts = " << natts << ", fails = " << netfs << ")\n"; @@ -811,17 +813,19 @@ int main(int argc, char* argv[]) static int ff(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { Options* opts = (Options*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; sunrealtype tmp1, tmp2; // fill in the RHS function: // [0 0]*[(-2+u^2-r(t))/(2*u)] + [ 0 ] // [e -1] [(-2+v^2-s(t))/(2*v)] [sdot(t)/(2v)] - tmp1 = (-TWO + u * u - r(t, opts)) / (TWO * u); - tmp2 = (-TWO + v * v - s(t, opts)) / (TWO * v); - NV_Ith_S(ydot, 0) = ZERO; - NV_Ith_S(ydot, 1) = opts->e * tmp1 - tmp2 + sdot(t, opts) / (TWO * v); + tmp1 = (-TWO + u * u - r(t, opts)) / (TWO * u); + tmp2 = (-TWO + v * v - s(t, opts)) / (TWO * v); + dydata[0] = ZERO; + dydata[1] = opts->e * tmp1 - tmp2 + sdot(t, opts) / (TWO * v); // Return with success return 0; @@ -831,17 +835,19 @@ static int ff(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int fs(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { Options* opts = (Options*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; sunrealtype tmp1, tmp2; // fill in the RHS function: // [G e]*[(-2+u^2-r(t))/(2*u)] + [rdot(t)/(2u)] // [0 0] [(-2+v^2-s(t))/(2*v)] [ 0 ] - tmp1 = (-TWO + u * u - r(t, opts)) / (TWO * u); - tmp2 = (-TWO + v * v - s(t, opts)) / (TWO * v); - NV_Ith_S(ydot, 0) = opts->G * tmp1 + opts->e * tmp2 + rdot(t, opts) / (TWO * u); - NV_Ith_S(ydot, 1) = ZERO; + tmp1 = (-TWO + u * u - r(t, opts)) / (TWO * u); + tmp2 = (-TWO + v * v - s(t, opts)) / (TWO * v); + dydata[0] = opts->G * tmp1 + opts->e * tmp2 + rdot(t, opts) / (TWO * u); + dydata[1] = ZERO; // Return with success return 0; @@ -851,13 +857,15 @@ static int fs(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int fse(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { Options* opts = (Options*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; // fill in the slow explicit RHS function: // [rdot(t)/(2u)] // [ 0 ] - NV_Ith_S(ydot, 0) = rdot(t, opts) / (TWO * u); - NV_Ith_S(ydot, 1) = ZERO; + dydata[0] = rdot(t, opts) / (TWO * u); + dydata[1] = ZERO; // Return with success return 0; @@ -867,17 +875,19 @@ static int fse(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int fsi(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { Options* opts = (Options*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; sunrealtype tmp1, tmp2; // fill in the slow implicit RHS function: // [G e]*[(-2+u^2-r(t))/(2*u)] // [0 0] [(-2+v^2-s(t))/(2*v)] - tmp1 = (-TWO + u * u - r(t, opts)) / (TWO * u); - tmp2 = (-TWO + v * v - s(t, opts)) / (TWO * v); - NV_Ith_S(ydot, 0) = opts->G * tmp1 + opts->e * tmp2; - NV_Ith_S(ydot, 1) = ZERO; + tmp1 = (-TWO + u * u - r(t, opts)) / (TWO * u); + tmp2 = (-TWO + v * v - s(t, opts)) / (TWO * v); + dydata[0] = opts->G * tmp1 + opts->e * tmp2; + dydata[1] = ZERO; // Return with success return 0; @@ -886,17 +896,19 @@ static int fsi(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int fn(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { Options* opts = (Options*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; sunrealtype tmp1, tmp2; // fill in the RHS function: // [G e]*[(-2+u^2-r(t))/(2*u)] + [rdot(t)/(2u)] // [e -1] [(-2+v^2-s(t))/(2*v)] [sdot(t)/(2v)] - tmp1 = (-TWO + u * u - r(t, opts)) / (TWO * u); - tmp2 = (-TWO + v * v - s(t, opts)) / (TWO * v); - NV_Ith_S(ydot, 0) = opts->G * tmp1 + opts->e * tmp2 + rdot(t, opts) / (TWO * u); - NV_Ith_S(ydot, 1) = opts->e * tmp1 - tmp2 + sdot(t, opts) / (TWO * v); + tmp1 = (-TWO + u * u - r(t, opts)) / (TWO * u); + tmp2 = (-TWO + v * v - s(t, opts)) / (TWO * v); + dydata[0] = opts->G * tmp1 + opts->e * tmp2 + rdot(t, opts) / (TWO * u); + dydata[1] = opts->e * tmp1 - tmp2 + sdot(t, opts) / (TWO * v); // Return with success return 0; @@ -912,8 +924,9 @@ static int Js(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { Options* opts = (Options*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; sunrealtype t11, t22; // fill in the Jacobian: @@ -934,8 +947,9 @@ static int Jsi(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { Options* opts = (Options*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; sunrealtype t11, t22; // fill in the Jacobian: @@ -956,8 +970,9 @@ static int Jn(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { Options* opts = (Options*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; sunrealtype t11, t22; // fill in the Jacobian: @@ -997,18 +1012,19 @@ static sunrealtype sdot(sunrealtype t, Options* opts) static sunrealtype utrue(sunrealtype t, Options* opts) { - return (SUNRsqrt(TWO + r(t, opts))); + return (std::sqrt(TWO + r(t, opts))); } static sunrealtype vtrue(sunrealtype t, Options* opts) { - return (SUNRsqrt(TWO + s(t, opts))); + return (std::sqrt(TWO + s(t, opts))); } static int Ytrue(sunrealtype t, N_Vector y, Options* opts) { - NV_Ith_S(y, 0) = utrue(t, opts); - NV_Ith_S(y, 1) = vtrue(t, opts); + sunrealtype* ydata = N_VGetArrayPointer(y); + ydata[0] = utrue(t, opts); + ydata[1] = vtrue(t, opts); return (0); } diff --git a/test/unit_tests/arkode/CXX_serial/ark_test_slowerror_brusselator.cpp b/test/unit_tests/arkode/CXX_serial/ark_test_slowerror_brusselator.cpp index cc0b1b5a11..8b74a377fe 100644 --- a/test/unit_tests/arkode/CXX_serial/ark_test_slowerror_brusselator.cpp +++ b/test/unit_tests/arkode/CXX_serial/ark_test_slowerror_brusselator.cpp @@ -252,11 +252,13 @@ int main(int argc, char* argv[]) if (check_retval((void*)y, "N_VNew_Serial", 0)) return 1; N_Vector* yref = N_VCloneVectorArray(udata.Npart + 1, y); if (check_retval((void*)yref, "N_VNew_Serial", 0)) return 1; + sunrealtype* ydata = N_VGetArrayPointer(y); + if (check_retval((void*)ydata, "N_VGetArrayPointer", 0)) return 1; + ydata[0] = u0; + ydata[1] = v0; + ydata[2] = w0; // Generate reference solution - NV_Ith_S(y, 0) = u0; - NV_Ith_S(y, 1) = v0; - NV_Ith_S(y, 2) = w0; void* arkode_ref = ERKStepCreate(fn, T0, y, ctx); if (check_retval((void*)arkode_ref, "ERKStepCreate", 0)) return 1; retval = ARKodeSetUserData(arkode_ref, (void*)&udata); @@ -280,9 +282,9 @@ int main(int argc, char* argv[]) } // Set up fast ERKStep integrator as fifth-order adaptive method - NV_Ith_S(y, 0) = u0; - NV_Ith_S(y, 1) = v0; - NV_Ith_S(y, 2) = w0; + ydata[0] = u0; + ydata[1] = v0; + ydata[2] = w0; void* inner_arkode_mem = ERKStepCreate(f0, T0, y, ctx); if (check_retval((void*)inner_arkode_mem, "ERKStepCreate", 0)) return 1; retval = ARKodeSetOrder(inner_arkode_mem, 5); @@ -374,15 +376,17 @@ static int f0(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int fn(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { - UserData* udata = (UserData*)user_data; - sunrealtype u = NV_Ith_S(y, 0); // access solution values - sunrealtype v = NV_Ith_S(y, 1); - sunrealtype w = NV_Ith_S(y, 2); + UserData* udata = (UserData*)user_data; + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; // access solution values + const sunrealtype v = ydata[1]; + const sunrealtype w = ydata[2]; // fill in the RHS function - NV_Ith_S(ydot, 0) = udata->a - (w + ONE) * u + v * u * u; - NV_Ith_S(ydot, 1) = w * u - v * u * u; - NV_Ith_S(ydot, 2) = (udata->b - w) / udata->ep - w * u; + dydata[0] = udata->a - (w + ONE) * u + v * u * u; + dydata[1] = w * u - v * u * u; + dydata[2] = (udata->b - w) / udata->ep - w * u; // Return with success return 0; @@ -390,14 +394,16 @@ static int fn(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int fi(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { - UserData* udata = (UserData*)user_data; - sunrealtype u = NV_Ith_S(y, 0); // access solution values - sunrealtype w = NV_Ith_S(y, 2); + UserData* udata = (UserData*)user_data; + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; // access solution values + const sunrealtype w = ydata[2]; // fill in the RHS function - NV_Ith_S(ydot, 0) = ZERO; - NV_Ith_S(ydot, 1) = ZERO; - NV_Ith_S(ydot, 2) = (udata->b - w) / udata->ep - w * u; + dydata[0] = ZERO; + dydata[1] = ZERO; + dydata[2] = (udata->b - w) / udata->ep - w * u; // Return with success return 0; @@ -405,15 +411,17 @@ static int fi(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int fe(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { - UserData* udata = (UserData*)user_data; - sunrealtype u = NV_Ith_S(y, 0); // access solution values - sunrealtype v = NV_Ith_S(y, 1); - sunrealtype w = NV_Ith_S(y, 2); + UserData* udata = (UserData*)user_data; + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; // access solution values + const sunrealtype v = ydata[1]; + const sunrealtype w = ydata[2]; // fill in the RHS function - NV_Ith_S(ydot, 0) = udata->a - (w + ONE) * u + v * u * u; - NV_Ith_S(ydot, 1) = w * u - v * u * u; - NV_Ith_S(ydot, 2) = ZERO; + dydata[0] = udata->a - (w + ONE) * u + v * u * u; + dydata[1] = w * u - v * u * u; + dydata[2] = ZERO; // Return with success return 0; @@ -422,10 +430,11 @@ static int fe(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int Jn(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { - UserData* udata = (UserData*)user_data; - sunrealtype u = NV_Ith_S(y, 0); // access solution values - sunrealtype v = NV_Ith_S(y, 1); - sunrealtype w = NV_Ith_S(y, 2); + UserData* udata = (UserData*)user_data; + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; // access solution values + const sunrealtype v = ydata[1]; + const sunrealtype w = ydata[2]; // fill in the Jacobian SM_ELEMENT_D(J, 0, 0) = -(w + ONE) + TWO * u * v; @@ -447,9 +456,10 @@ static int Jn(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, static int Ji(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { - UserData* udata = (UserData*)user_data; - sunrealtype u = NV_Ith_S(y, 0); // access solution values - sunrealtype w = NV_Ith_S(y, 2); + UserData* udata = (UserData*)user_data; + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; // access solution values + const sunrealtype w = ydata[2]; // fill in the Jacobian SM_ELEMENT_D(J, 0, 0) = ZERO; @@ -481,10 +491,12 @@ static int run_test(void* mristep_mem, void* arkode_ref, N_Vector y, int retval; sunrealtype hpart = (Tf - T0) / udata.Npart; sunrealtype t, t2; - N_Vector y2 = N_VClone(y); - N_Vector ele = N_VClone(y); - N_Vector ewt = N_VClone(y); - N_Vector vtemp = N_VClone(y); + N_Vector y2 = N_VClone(y); + N_Vector ele = N_VClone(y); + N_Vector ewt = N_VClone(y); + N_Vector vtemp = N_VClone(y); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* y2data = N_VGetArrayPointer(y2); // Set storage for errors vector> dsm(Hvals.size(), @@ -525,12 +537,12 @@ static int run_test(void* mristep_mem, void* arkode_ref, N_Vector y, dsm_est[iH][ipart] = N_VWrmsNorm(ewt, ele); // Compute/print solution error - sunrealtype udsm = abs(NV_Ith_S(y, 0) - NV_Ith_S(y2, 0)) / - (abstol + reltol * abs(NV_Ith_S(y2, 0))); - sunrealtype vdsm = abs(NV_Ith_S(y, 1) - NV_Ith_S(y2, 1)) / - (abstol + reltol * abs(NV_Ith_S(y2, 1))); - sunrealtype wdsm = abs(NV_Ith_S(y, 2) - NV_Ith_S(y2, 2)) / - (abstol + reltol * abs(NV_Ith_S(y2, 2))); + sunrealtype udsm = abs(ydata[0] - y2data[0]) / + (abstol + reltol * abs(y2data[0])); + sunrealtype vdsm = abs(ydata[1] - y2data[1]) / + (abstol + reltol * abs(y2data[1])); + sunrealtype wdsm = abs(ydata[2] - y2data[2]) / + (abstol + reltol * abs(y2data[2])); dsm[iH][ipart] = sqrt((udsm * udsm + vdsm * vdsm + wdsm * wdsm) / SUN_RCONST(3.0)); cout << " H " << Hvals[iH] << " method " << method << " t " << t diff --git a/test/unit_tests/arkode/CXX_serial/ark_test_slowerror_kpr.cpp b/test/unit_tests/arkode/CXX_serial/ark_test_slowerror_kpr.cpp index ef9c5e97e6..583be495ef 100644 --- a/test/unit_tests/arkode/CXX_serial/ark_test_slowerror_kpr.cpp +++ b/test/unit_tests/arkode/CXX_serial/ark_test_slowerror_kpr.cpp @@ -309,17 +309,19 @@ static int f0(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int fn(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { UserData* udata = (UserData*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; sunrealtype tmp1, tmp2; // fill in the RHS function: // [G e]*[(-2+u^2-p(t))/(2*u)] + [pdot(t)/(2u)] // [e -1] [(-2+v^2-s(t))/(2*v)] [qdot(t)/(2v)] - tmp1 = (-TWO + u * u - p(t)) / (TWO * u); - tmp2 = (-TWO + v * v - q(t, *udata)) / (TWO * v); - NV_Ith_S(ydot, 0) = udata->G * tmp1 + udata->e * tmp2 + pdot(t) / (TWO * u); - NV_Ith_S(ydot, 1) = udata->e * tmp1 - tmp2 + qdot(t, *udata) / (TWO * v); + tmp1 = (-TWO + u * u - p(t)) / (TWO * u); + tmp2 = (-TWO + v * v - q(t, *udata)) / (TWO * v); + dydata[0] = udata->G * tmp1 + udata->e * tmp2 + pdot(t) / (TWO * u); + dydata[1] = udata->e * tmp1 - tmp2 + qdot(t, *udata) / (TWO * v); // Return with success return 0; @@ -328,17 +330,19 @@ static int fn(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int fi(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { UserData* udata = (UserData*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; sunrealtype tmp1, tmp2; // fill in the RHS function: // [G e]*[(-2+u^2-p(t))/(2*u)] + [pdot(t)/(2u)] // [0 0] [(-2+v^2-s(t))/(2*v)] [ 0 ] - tmp1 = (-TWO + u * u - p(t)) / (TWO * u); - tmp2 = (-TWO + v * v - q(t, *udata)) / (TWO * v); - NV_Ith_S(ydot, 0) = udata->G * tmp1 + udata->e * tmp2 + pdot(t) / (TWO * u); - NV_Ith_S(ydot, 1) = ZERO; + tmp1 = (-TWO + u * u - p(t)) / (TWO * u); + tmp2 = (-TWO + v * v - q(t, *udata)) / (TWO * v); + dydata[0] = udata->G * tmp1 + udata->e * tmp2 + pdot(t) / (TWO * u); + dydata[1] = ZERO; // Return with success return 0; @@ -347,17 +351,19 @@ static int fi(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int fe(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { UserData* udata = (UserData*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* dydata = N_VGetArrayPointer(ydot); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; sunrealtype tmp1, tmp2; // fill in the RHS function: // [0 0]*[(-2+u^2-p(t))/(2*u)] + [ 0 ] // [e -1] [(-2+v^2-s(t))/(2*v)] [qdot(t)/(2v)] - tmp1 = (-TWO + u * u - p(t)) / (TWO * u); - tmp2 = (-TWO + v * v - q(t, *udata)) / (TWO * v); - NV_Ith_S(ydot, 0) = ZERO; - NV_Ith_S(ydot, 1) = udata->e * tmp1 - tmp2 + qdot(t, *udata) / (TWO * v); + tmp1 = (-TWO + u * u - p(t)) / (TWO * u); + tmp2 = (-TWO + v * v - q(t, *udata)) / (TWO * v); + dydata[0] = ZERO; + dydata[1] = udata->e * tmp1 - tmp2 + qdot(t, *udata) / (TWO * v); // Return with success return 0; @@ -367,8 +373,9 @@ static int Jn(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { UserData* udata = (UserData*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; sunrealtype t11, t22; // fill in the Jacobian: @@ -389,8 +396,9 @@ static int Ji(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { UserData* udata = (UserData*)user_data; - const sunrealtype u = NV_Ith_S(y, 0); - const sunrealtype v = NV_Ith_S(y, 1); + sunrealtype* ydata = N_VGetArrayPointer(y); + const sunrealtype u = ydata[0]; + const sunrealtype v = ydata[1]; // fill in the Jacobian: // [G/2 + (G*(1+p(t))-pdot(t))/(2*u^2) e/2+e*(2+s(t))/(2*v^2)] @@ -418,9 +426,10 @@ static int run_test(void* mristep_mem, N_Vector y, sunrealtype T0, int retval; sunrealtype hpart = (Tf - T0) / udata.Npart; sunrealtype t; - N_Vector ele = N_VClone(y); - N_Vector ewt = N_VClone(y); - N_Vector vtemp = N_VClone(y); + N_Vector ele = N_VClone(y); + N_Vector ewt = N_VClone(y); + N_Vector vtemp = N_VClone(y); + sunrealtype* ydata = N_VGetArrayPointer(y); // Set storage for errors vector> dsm(Hvals.size(), @@ -455,9 +464,9 @@ static int run_test(void* mristep_mem, N_Vector y, sunrealtype T0, dsm_est[iH][ipart] = N_VWrmsNorm(ewt, ele); // Compute/print solution error - sunrealtype udsm = abs(NV_Ith_S(y, 0) - utrue(t)) / + sunrealtype udsm = abs(ydata[0] - utrue(t)) / (abstol + reltol * abs(utrue(t))); - sunrealtype vdsm = abs(NV_Ith_S(y, 1) - vtrue(t, udata)) / + sunrealtype vdsm = abs(ydata[1] - vtrue(t, udata)) / (abstol + reltol * abs(vtrue(t, udata))); dsm[iH][ipart] = sqrt(0.5 * (udsm * udsm + vdsm * vdsm)); cout << " H " << Hvals[iH] << " method " << method << " t " << t @@ -505,8 +514,9 @@ static sunrealtype vtrue(sunrealtype t, UserData& udata) static int Ytrue(sunrealtype t, N_Vector y, UserData& udata) { - NV_Ith_S(y, 0) = utrue(t); - NV_Ith_S(y, 1) = vtrue(t, udata); + sunrealtype* ydata = N_VGetArrayPointer(y); + ydata[0] = utrue(t); + ydata[1] = vtrue(t, udata); return (0); } diff --git a/test/unit_tests/arkode/CXX_serial/ark_test_slowerror_polynomial.cpp b/test/unit_tests/arkode/CXX_serial/ark_test_slowerror_polynomial.cpp index 1f236da7d6..96541382b2 100644 --- a/test/unit_tests/arkode/CXX_serial/ark_test_slowerror_polynomial.cpp +++ b/test/unit_tests/arkode/CXX_serial/ark_test_slowerror_polynomial.cpp @@ -240,8 +240,9 @@ static int f0(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) static int fn(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) { // fill in the RHS function and return with success - UserData* udata = (UserData*)user_data; - NV_Ith_S(ydot, 0) = udata->a * t * t + udata->b * t + udata->c; + UserData* udata = (UserData*)user_data; + sunrealtype* dydata = N_VGetArrayPointer(ydot); + dydata[0] = udata->a * t * t + udata->b * t + udata->c; return 0; } @@ -264,9 +265,11 @@ static int run_test(void* mristep_mem, N_Vector y, sunrealtype T0, // Reused variables int retval; sunrealtype t; - N_Vector vtemp = N_VClone(y); - N_Vector ele = N_VClone(y); - N_Vector ewt = N_VClone(y); + N_Vector vtemp = N_VClone(y); + N_Vector ele = N_VClone(y); + N_Vector ewt = N_VClone(y); + sunrealtype* ydata = N_VGetArrayPointer(y); + sunrealtype* vdata = N_VGetArrayPointer(vtemp); // Set storage for errors vector dsm(Hvals.size(), ZERO); @@ -298,18 +301,17 @@ static int run_test(void* mristep_mem, N_Vector y, sunrealtype T0, // Compute/print solution error retval = Ytrue(t, vtemp, udata); if (check_retval(&retval, "Ytrue", 1)) return 1; - dsm[iH] = abs(NV_Ith_S(y, 0) - NV_Ith_S(vtemp, 0)) / - (abstol + reltol * abs(NV_Ith_S(vtemp, 0))); + dsm[iH] = abs(ydata[0] - vdata[0]) / (abstol + reltol * abs(vdata[0])); if (method == "ARKODE_MRI_GARK_ERK22a") { printf(" H %.5f dsm %.2e dsm_est %.2e dsm_anal %.2e " " dsm_est_anal %.2e\n", Hvals[iH], dsm[iH], dsm_est[iH], Hvals[iH] * Hvals[iH] * Hvals[iH] * abs(udata.a / 12.0) / - (abstol + reltol * abs(NV_Ith_S(vtemp, 0))), + (abstol + reltol * abs(vdata[0])), Hvals[iH] * Hvals[iH] * abs(udata.a * Hvals[iH] / 4.0 + udata.b / 2.0) / - (abstol + reltol * abs(NV_Ith_S(vtemp, 0)))); + (abstol + reltol * abs(vdata[0]))); } else if (method == "ARKODE_MRI_GARK_IRK21a") { @@ -317,10 +319,10 @@ static int run_test(void* mristep_mem, N_Vector y, sunrealtype T0, " dsm_est_anal %.2e\n", Hvals[iH], dsm[iH], dsm_est[iH], Hvals[iH] * Hvals[iH] * Hvals[iH] * abs(udata.a / 6.0) / - (abstol + reltol * abs(NV_Ith_S(vtemp, 0))), + (abstol + reltol * abs(vdata[0])), Hvals[iH] * Hvals[iH] * abs(udata.a * Hvals[iH] / 2.0 + udata.b / 2.0) / - (abstol + reltol * abs(NV_Ith_S(vtemp, 0)))); + (abstol + reltol * abs(vdata[0]))); } else if (method == "ARKODE_IMEX_MRI_SR21") { @@ -328,11 +330,11 @@ static int run_test(void* mristep_mem, N_Vector y, sunrealtype T0, " dsm_est_anal %.2e\n", Hvals[iH], dsm[iH], dsm_est[iH], Hvals[iH] * Hvals[iH] * Hvals[iH] * abs(udata.a * 3137.0 / 50370.0) / - (abstol + reltol * abs(NV_Ith_S(vtemp, 0))), + (abstol + reltol * abs(vdata[0])), Hvals[iH] * Hvals[iH] * abs(udata.a * Hvals[iH] * 20191.0 / 755550.0 - udata.b * 19.0 / 30.0) / - (abstol + reltol * abs(NV_Ith_S(vtemp, 0)))); + (abstol + reltol * abs(vdata[0]))); } else { @@ -349,8 +351,9 @@ static int run_test(void* mristep_mem, N_Vector y, sunrealtype T0, static int Ytrue(sunrealtype t, N_Vector y, UserData& udata) { - NV_Ith_S(y, 0) = udata.a / SUN_RCONST(3.0) * t * t * t + - udata.b / SUN_RCONST(2.0) * t * t + udata.c * t + ONE; + sunrealtype* ydata = N_VGetArrayPointer(y); + ydata[0] = udata.a / SUN_RCONST(3.0) * t * t * t + + udata.b / SUN_RCONST(2.0) * t * t + udata.c * t + ONE; return (0); } diff --git a/test/unit_tests/arkode/C_serial/ark_test_erkstepsetforcing.c b/test/unit_tests/arkode/C_serial/ark_test_erkstepsetforcing.c deleted file mode 100644 index 1a40f17082..0000000000 --- a/test/unit_tests/arkode/C_serial/ark_test_erkstepsetforcing.c +++ /dev/null @@ -1,386 +0,0 @@ -/* ----------------------------------------------------------------------------- - * Programmer(s): Daniel R. Reynolds @ SMU - * ----------------------------------------------------------------------------- - * SUNDIALS Copyright Start - * Copyright (c) 2002-2023, 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 - * ----------------------------------------------------------------------------- - * Unit test for erkStep_SetInnerForcing - * - * erkStep_SetInnerForcing is used to provide a polynomial forcing term in - * ERKStep when it is used as the inner integrator under MRIStep. To check that - * the forcing is computed and applied correctly we integrate an ODE in time - * with ERKStep where to RHS consists of only the polynomial forcing term. The - * solution should be exact when the method order is greater than or equal to - * the polynomial order. - * ---------------------------------------------------------------------------*/ - -#include -#include -#include - -#include "arkode/arkode_erkstep.h" -#include "arkode/arkode_erkstep_impl.h" -#include "nvector/nvector_serial.h" -#include "sundials/sundials_math.h" -#include "sundials/sundials_types.h" - -#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 Called by the Solver */ -static int f(realtype t, N_Vector y, N_Vector ydot, void* user_data); - -/* Private function to check function return values */ -static int check_flag(void* flagvalue, const char* funcname, int opt); - -/* Private function to check computed solution */ -static int compute_ans(realtype t, realtype tshift, realtype tscale, - N_Vector* forcing, int order, N_Vector ans); -static int compute_error(N_Vector y, N_Vector ans, N_Vector tmp, realtype rtol, - realtype atol); - -/* Main Program */ -int main(int argc, char* argv[]) -{ - SUNContext sunctx = NULL; - - /* default input values */ - sunindextype NEQ = 1; /* number of dependent vars. */ - int order = 3; /* order of polynomial forcing */ - realtype T0 = RCONST(0.0); /* initial time */ - realtype Tf = RCONST(1.0); /* final time */ - realtype tshift = T0; /* time shift for normalization */ - realtype tscale = Tf; /* time scale for normalization */ - - /* tolerances */ - realtype reltol = SUNRsqrt(UNIT_ROUNDOFF); - realtype abstol = SUNRsqrt(UNIT_ROUNDOFF) / 100; - - /* general problem variables */ - int flag; /* reusable error-checking flag */ - N_Vector y = NULL; /* vector for storing the computed solution */ - N_Vector ans = NULL; /* vector for storing the true solution */ - N_Vector tmp = NULL; /* temporary workspace vector */ - N_Vector* forcing = NULL; /* array of forcing vectors */ - void* arkode_mem = NULL; /* ARKode memory structure */ - int i, j; /* loop counters */ - realtype tret; /* integrator return time */ - realtype* data; /* array for accessing vector data */ - long int nst, nst_a; /* number of integrator steps */ - long int mxsteps = 100000; /* max steps before output */ - - /* check inputs */ - if (argc > 1) - { - NEQ = (sunindextype)atol(argv[1]); - if (NEQ <= 0) - { - printf("ERROR: The problem size must be a positive integer\n"); - return (1); - } - } - - if (argc > 2) - { - order = (int)atol(argv[2]); - if (order < 0) - { - printf("ERROR: The polynomial order must be a non-negative integer\n"); - return (1); - } - } - - if (argc > 3) - { - if (argc < 5) - { - printf("ERROR: Both the initial and final time are required\n"); - return (1); - } - T0 = (realtype)atof(argv[3]); - Tf = (realtype)atof(argv[4]); - if (SUNRabs(T0) >= SUNRabs(Tf)) - { - printf("ERROR: |T0| must be less than |Tf|\n"); - return (1); - } - } - - if (argc > 5) - { - if (argc < 7) - { - printf("ERROR: Both tshift and tscale are required\n"); - return (1); - } - tshift = (realtype)atof(argv[5]); - tscale = (realtype)atof(argv[6]); - if (SUNRabs(tscale) < TINY) - { - printf("ERROR: |tscale| must be greater than %" GSYM "\n", TINY); - return (1); - } - } - - /* Output test setup */ - printf("\nerkStep_SetInnerForcing unit test:\n"); - printf(" NEQ = %li\n", (long int)NEQ); - printf(" order = %i\n", order); - printf(" t0 = %.1" GSYM "\n", T0); - printf(" tf = %.1" GSYM "\n", Tf); - printf(" tshift = %.1" GSYM "\n", tshift); - printf(" tscale = %.1" GSYM "\n", tscale); - printf(" reltol = %.1" ESYM "\n", reltol); - printf(" abstol = %.1" ESYM "\n\n", abstol); - - /* Create the SUNDIALS context object for this simulation. */ - SUNContext_Create(NULL, &sunctx); - - /* Create solution vector and initialize to zero */ - y = N_VNew_Serial(NEQ, sunctx); - if (check_flag((void*)y, "N_VNew_Serial", 0)) return 1; - - ans = N_VClone(y); - if (check_flag((void*)ans, "N_VClone", 0)) return 1; - - tmp = N_VClone(y); - if (check_flag((void*)tmp, "N_VClone", 0)) return 1; - - /* allocate vector array for polynomial forcing */ - forcing = N_VCloneVectorArray(order + 1, y); - if (check_flag((void*)forcing, "N_VCloneVectorArray", 0)) return 1; - for (i = 0; i < order + 1; i++) - { - if (check_flag((void*)forcing[i], "N_VCloneVectorArray", 0)) return 1; - } - - /* fill forcing vectors with random data */ - for (i = 0; i < order + 1; i++) - { - data = N_VGetArrayPointer(forcing[i]); - for (j = 0; j < NEQ; j++) - { - data[j] = (realtype)rand() / (realtype)RAND_MAX; - } - } - - /* compute the true solution */ - flag = compute_ans(Tf, tshift, tscale, forcing, order, ans); - if (check_flag(&flag, "compute_ans", 1)) return 1; - - printf("True solution:\n"); - N_VPrint_Serial(ans); - - /* --------------------------------------------------------------------------- - * explicit test - * -------------------------------------------------------------------------*/ - - /* initialize solution vector */ - flag = compute_ans(T0, tshift, tscale, forcing, order, y); - if (check_flag(&flag, "compute_ans", 1)) return 1; - - /* Create ARKode mem structure */ - arkode_mem = ERKStepCreate(f, T0, y, sunctx); - if (check_flag((void*)arkode_mem, "ERKStepCreate", 0)) return 1; - - /* Specify tolerances */ - flag = ARKodeSStolerances(arkode_mem, reltol, abstol); - if (check_flag(&flag, "ARKodeSStolerances", 1)) return 1; - - /* Set stop time */ - flag = ARKodeSetStopTime(arkode_mem, Tf); - if (check_flag(&flag, "ARKodeSetStopTime", 1)) return 1; - - /* Set max steps before output */ - flag = ARKodeSetMaxNumSteps(arkode_mem, mxsteps); - if (check_flag(&flag, "ARKodeSetMaxNumSteps", 1)) return 1; - - /* Set forcing */ - flag = erkStep_SetInnerForcing(arkode_mem, tshift, tscale, forcing, order + 1); - if (check_flag(&flag, "erkStep_SetInnerForcing", 1)) return 1; - - /* Integrate the problem */ - flag = ARKodeEvolve(arkode_mem, Tf, y, &tret, ARK_NORMAL); - - /* check for errors */ - if (flag < 0) - { - fprintf(stderr, "ARKodeEvolve failure, flag = %d\n", flag); - return 1; - } - - /* get some integrator stats */ - flag = ARKodeGetNumSteps(arkode_mem, &nst); - check_flag(&flag, "ARKodeGetNumSteps", 1); - - flag = ARKodeGetNumStepAttempts(arkode_mem, &nst_a); - check_flag(&flag, "ARKodeGetNumStepAttempts", 1); - - printf("Stats:\n"); - printf("Steps = %li (attempted = %li)\n\n", nst, nst_a); - - /* Free integrator memory */ - ARKodeFree(&arkode_mem); - arkode_mem = NULL; - - /* print solution */ - printf("Explicit solution:\n"); - N_VPrint_Serial(y); - - /* check the solution error */ - flag = compute_error(y, ans, tmp, reltol, abstol); - if (flag != 0) return (1); - - /* --------------------------------------------------------------------------- - * Clean up and return - * -------------------------------------------------------------------------*/ - - N_VDestroy(y); /* Free vectors */ - N_VDestroy(ans); - N_VDestroy(tmp); - - N_VDestroyVectorArray(forcing, order + 1); - - SUNContext_Free(&sunctx); - return flag; -} - -/*------------------------------- - * Functions called by the solver - *-------------------------------*/ - -/* ODE RHS function */ -static int f(realtype t, N_Vector y, N_Vector ydot, void* user_data) -{ - N_VConst(ZERO, ydot); - return 0; -} - -/*------------------------------- - * Private helper functions - *-------------------------------*/ - -/* Check function return value... - opt == 0 means SUNDIALS function allocates memory so check if - returned NULL pointer - opt == 1 means SUNDIALS function returns a flag so check if - flag < 0 - opt == 2 means function allocates memory so check if returned - NULL pointer -*/ -static int check_flag(void* flagvalue, const char* funcname, int opt) -{ - int* errflag; - - /* Check if SUNDIALS function returned NULL pointer - no memory allocated */ - if (opt == 0 && flagvalue == NULL) - { - fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n", - funcname); - return 1; - } - - /* Check if flag < 0 */ - else if (opt == 1) - { - errflag = (int*)flagvalue; - if (*errflag < 0) - { - fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n", - funcname, *errflag); - return 1; - } - } - - /* Check if function returned NULL pointer - no memory allocated */ - else if (opt == 2 && flagvalue == NULL) - { - fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n", - funcname); - return 1; - } - - return 0; -} - -/* computed the true solution at a given time */ -static int compute_ans(realtype t, realtype tshift, realtype tscale, - N_Vector* forcing, int order, N_Vector ans) -{ - int i; - realtype tau; - N_Vector* vecs; - realtype* vals; - - vals = NULL; - vals = (realtype*)calloc(order + 1, sizeof(realtype)); - if (vals == NULL) return (1); - - vecs = NULL; - vecs = (N_Vector*)calloc(order + 1, sizeof(N_Vector)); - if (vecs == NULL) return (1); - - /* compute normalized time */ - tau = (t - tshift) / tscale; - - /* compute true solution */ - for (i = 0; i < order + 1; i++) - { - vals[i] = ((SUNRpowerI(tau, i + 1)) / (i + 1)) * tscale; - vecs[i] = forcing[i]; - } - N_VLinearCombination(order + 1, vals, vecs, ans); - - free(vals); - free(vecs); - - return (0); -} - -/* compure the weighted max norm of the difference of two vectors */ -static int compute_error(N_Vector y, N_Vector ans, N_Vector tmp, realtype rtol, - realtype atol) -{ - int status; /* success (0) or failure (1) flag */ - realtype error; - - /* compute the error in y */ - N_VLinearSum(ONE, y, -ONE, ans, y); - - /* compute error weights */ - N_VAbs(ans, tmp); - N_VScale(rtol, tmp, tmp); - N_VAddConst(tmp, atol, tmp); - N_VInv(tmp, tmp); - - /* compute weighted max norm */ - N_VProd(tmp, y, y); - error = N_VMaxNorm(y); - - /* is the solution within the tolerances? */ - status = (error < ONE) ? 0 : 1; - - if (status) - fprintf(stdout, "ERROR: Test failed with wmax error = %" GSYM "\n\n", error); - else - fprintf(stdout, "SUCCESS: Test passed with wmax error = %" GSYM "\n\n", - error); - - return (status); -} - -/*---- end of file ----*/