add access to linear residual

examples/arkode/C_serial/ark_analytic.c

@@ -54,15 +54,21 @@ static int f(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);
 
-/* >>> Function to access delta <<< */
-int Access(long int step, int stage, int iter, N_Vector delta, void* user_data)
+/* >>> Function to access delta and linear residuals <<< */
+int Access(long int step, int stage, int iter, N_Vector delta,
+ N_Vector ls_res, N_Vector ls_res_r, void* user_data)
 {
  printf("Step %li = \n", step);
  printf("Stage %i = \n", stage);
  printf("Iter %i = \n", iter);
  printf("Delta:\n");
  N_VPrintFile(delta, stdout);
+ printf("LS Res:\n");
+ N_VPrintFile(ls_res, stdout);
+ printf("LS Res Relaxed:\n");
+ N_VPrintFile(ls_res_r, stdout);
  printf("\n");
+ return 0;
 }
 
 /* Private function to check function return values */
@@ -142,7 +148,7 @@ int main(void)
  if (check_flag(&flag, "ARKodeSetLinear", 1)) { return 1; }
 
  /* >>> Attach function to access delta <<< */
- flag = ARKStepSetAccessDeltaFn(arkode_mem, Access);
+ flag = ARKStepSetAccessFn(arkode_mem, Access);
  if (check_flag(&flag, "ARKStepSetAccessDeltaFn", 1)) { return 1; }
 
  /* Open output stream for results, output comment line */

include/arkode/arkode_arkstep.h

@@ -30,8 +30,9 @@ extern "C" {
 #endif
 
 /* Callback to access delta */
-typedef int (*ARKStepAccessDeltaFn)(long int step, int stage, int iter,
- N_Vector delta, void* user_data);
+typedef int (*ARKStepAccessFn)(long int step, int stage, int iter,
+ N_Vector delta, N_Vector ls_res,
+ N_Vector ls_res2, void* user_data);
 
 /* -----------------
  * ARKStep Constants
@@ -106,7 +107,7 @@ SUNDIALS_EXPORT int ARKStepCreateMRIStepInnerStepper(void* arkode_mem,
  MRIStepInnerStepper* stepper);
 
 SUNDIALS_EXPORT
-int ARKStepSetAccessDeltaFn(void* arkode_mem, ARKStepAccessDeltaFn access_fn);
+int ARKStepSetAccessFn(void* arkode_mem, ARKStepAccessFn access_fn);
 
 /* --------------------------------------------------------------------------
  * Deprecated Functions -- all are superseded by shared ARKODE-level routines

src/arkode/arkode_arkstep_impl.h

@@ -161,7 +161,7 @@ typedef struct ARKodeARKStepMemRec
  sunrealtype* stage_times; /* workspace for applying forcing */
  sunrealtype* stage_coefs; /* workspace for applying forcing */
 
- ARKStepAccessDeltaFn access_fn;
+ ARKStepAccessFn access_fn;
 
 }* ARKodeARKStepMem;
 

src/arkode/arkode_arkstep_io.c

@@ -2408,7 +2408,7 @@ int ARKStepGetNumRelaxSolveIters(void* arkode_mem, long int* iters)
  EOF
  ===============================================================*/
 
-int ARKStepSetAccessDeltaFn(void* arkode_mem, ARKStepAccessDeltaFn access_fn)
+int ARKStepSetAccessFn(void* arkode_mem, ARKStepAccessFn access_fn)
 {
  ARKodeMem ark_mem;
  ARKodeARKStepMem step_mem;

src/arkode/arkode_arkstep_nls.c

@@ -26,13 +26,19 @@
 
 int arkAccessDeltaFn(int iter, N_Vector delta, void* arkode_mem)
 {
+ int retval = 0;
  ARKodeMem ark_mem = (ARKodeMem)arkode_mem;
 
  ARKodeARKStepMem step_mem;
- int retval = arkStep_AccessStepMem(ark_mem, __func__, &step_mem);
+ retval = arkStep_AccessStepMem(ark_mem, __func__, &step_mem);
  if (retval != ARK_SUCCESS) { return retval; }
 
+ ARKLsMem arkls_mem;
+ retval = arkLs_AccessLMem(ark_mem, __func__, &arkls_mem);
+ if (retval != ARK_SUCCESS) { return (retval); }
+
  retval = step_mem->access_fn(ark_mem->nst, step_mem->istage, iter, delta,
+ arkls_mem->ytemp, arkls_mem->ytemp2,
  ark_mem->user_data);
  if (retval != ARK_SUCCESS) { return retval; }
 

src/arkode/arkode_ls.c

@@ -24,6 +24,8 @@
 
 #include "arkode_impl.h"
 #include "arkode_ls_impl.h"
+#include "sundials/sundials_matrix.h"
+#include "sundials/sundials_nvector.h"
 
 /* constants */
 #define MIN_INC_MULT SUN_RCONST(1000.0)
@@ -264,6 +266,15 @@ int ARKodeSetLinearSolver(void* arkode_mem, SUNLinearSolver LS, SUNMatrix A)
  return (ARKLS_MEM_FAIL);
  }
 
+ if (!arkAllocVec(ark_mem, ark_mem->tempv1, &(arkls_mem->ytemp2)))
+ {
+ arkProcessError(ark_mem, ARKLS_MEM_FAIL, __LINE__, __func__, __FILE__,
+ MSG_LS_MEM_FAIL);
+ free(arkls_mem);
+ arkls_mem = NULL;
+ return (ARKLS_MEM_FAIL);
+ }
+
  if (!arkAllocVec(ark_mem, ark_mem->tempv1, &(arkls_mem->x)))
  {
  arkProcessError(ark_mem, ARKLS_MEM_FAIL, __LINE__, __func__, __FILE__,
@@ -3383,9 +3394,12 @@ int arkLsSolve(ARKodeMem ark_mem, N_Vector b, sunrealtype tnow, N_Vector ynow,
  }
  }
 
- /* Call solver, and copy x to b */
+ /* Call solver */
  retval = SUNLinSolSolve(arkls_mem->LS, arkls_mem->A, arkls_mem->x, b, delta);
- N_VScale(ONE, arkls_mem->x, b);
+
+ /* compute the residual r = Ax - b */
+ SUNMatMatvec(arkls_mem->savedJ, arkls_mem->x, arkls_mem->ytemp);
+ N_VLinearSum(ONE, arkls_mem->ytemp, -ONE, b, arkls_mem->ytemp);
 
  /* If using a direct or matrix-iterative solver, scale the correction to
  account for change in gamma (this is only beneficial if M==I) */
@@ -3399,9 +3413,16 @@ int arkLsSolve(ARKodeMem ark_mem, N_Vector b, sunrealtype tnow, N_Vector ynow,
  __FILE__, "An error occurred in ark_step_getgammas");
  return (arkls_mem->last_flag);
  }
- if (gamrat != ONE) { N_VScale(TWO / (ONE + gamrat), b, b); }
+ if (gamrat != ONE) { N_VScale(TWO / (ONE + gamrat), arkls_mem->x, arkls_mem->x); }
  }
 
+ /* compute the relaxed residual r = Ax - b */
+ SUNMatMatvec(arkls_mem->savedJ, arkls_mem->x, arkls_mem->ytemp2);
+ N_VLinearSum(ONE, arkls_mem->ytemp2, -ONE, b, arkls_mem->ytemp2);
+
+ /* copy x to b */
+ N_VScale(ONE, arkls_mem->x, b);
+
  /* Retrieve statistics from iterative linear solvers */
  resnorm = ZERO;
  nli_inc = 0;
@@ -3498,6 +3519,11 @@ int arkLsFree(ARKodeMem ark_mem)
  N_VDestroy(arkls_mem->ytemp);
  arkls_mem->ytemp = NULL;
  }
+ if (arkls_mem->ytemp2)
+ {
+ N_VDestroy(arkls_mem->ytemp2);
+ arkls_mem->ytemp2 = NULL;
+ }
  if (arkls_mem->x)
  {
  N_VDestroy(arkls_mem->x);

src/arkode/arkode_ls_impl.h

@@ -68,6 +68,7 @@ typedef struct ARKLsMemRec
  SUNMatrix A; /* A = M - gamma * df/dy */
  SUNMatrix savedJ; /* savedJ = old Jacobian */
  N_Vector ytemp; /* temp vector passed to jtimes and psolve */
+ N_Vector ytemp2;
  N_Vector x; /* solution vector used by SUNLinearSolver */
  N_Vector ycur; /* ptr to current y vector in ARKLs solve */
  N_Vector fcur; /* ptr to current fcur = fI(tcur, ycur) */