revision

LLNL · Sep 19, 2024 · 7cd161f · 7cd161f
1 parent 90c4ea6
commit 7cd161f
Show file tree

Hide file tree

Showing 5 changed files with 24 additions and 31 deletions.
diff --git a/doc/arkode/guide/source/Usage/index.rst b/doc/arkode/guide/source/Usage/index.rst
@@ -77,6 +77,6 @@ ARKBBDPRE can only be used with NVECTOR_PARALLEL.
    Preconditioners
    ARKStep/index.rst
    ERKStep/index.rst
+   LSRKStep/index.rst
    SPRKStep/index.rst
    MRIStep/index.rst
-   LSRKStep/index.rst
diff --git a/examples/arkode/CXX_parallel/ark_heat2D_lsrk_p.cpp b/examples/arkode/CXX_parallel/ark_heat2D_lsrk_p.cpp
@@ -196,8 +196,9 @@ struct UserData
 static int f(sunrealtype t, N_Vector u, N_Vector f, void* user_data);
 
 // Spectral radius estimation routine
-static int eig(sunrealtype t, N_Vector y, sunrealtype* lambdaR,
-               sunrealtype* lambdaI, void* user_data);
+static int eig(sunrealtype t, N_Vector y, N_Vector fn, sunrealtype* lambdaR,
+               sunrealtype* lambdaI, void* user_data, N_Vector temp1,
+               N_Vector temp2, N_Vector temp3);
 
 // -----------------------------------------------------------------------------
 // Helper functions
@@ -1120,8 +1121,9 @@ static int WaitRecv(UserData* udata)
 }
 
 // Spectral radius estimation routine
-static int eig(sunrealtype t, N_Vector y, sunrealtype* lambdaR,
-               sunrealtype* lambdaI, void* user_data)
+static int eig(sunrealtype t, N_Vector y, N_Vector fn, sunrealtype* lambdaR,
+               sunrealtype* lambdaI, void* user_data, N_Vector temp1,
+               N_Vector temp2, N_Vector temp3)
 {
   // Access problem data
   UserData* udata = (UserData*)user_data;

diff --git a/examples/arkode/C_serial/lsrk_analytic.c b/examples/arkode/C_serial/lsrk_analytic.c
@@ -127,7 +127,7 @@ int main(void)
   if (check_flag(&flag, "LSRKStepSetDomEigFn", 1)) { return 1; }
 
   /* Specify after how many successful steps dom_eig is recomputed 
-     Note that nsteps = 0 refers to Constant Jacobian */
+     Note that nsteps = 0 refers to constant dominant eigenvalue */
   flag = LSRKStepSetDomEigFrequency(arkode_mem, 0);
   if (check_flag(&flag, "LSRKStepSetDomEigFrequency", 1)) { return 1; }
 

diff --git a/examples/arkode/C_serial/lsrk_analytic_varjac.c b/examples/arkode/C_serial/lsrk_analytic_varjac.c
@@ -157,7 +157,7 @@ int main(void)
   if (check_flag(&flag, "LSRKStepSetDomEigSafetyFactor", 1)) { return 1; }
 
   /* Specify the Runge--Kutta--Chebyshev LSRK method */
-  flag = LSRKStepSetMethod(arkode_mem, ARKODE_LSRK_RKC_2);
+  flag = LSRKStepSetMethod(arkode_mem, ARKODE_LSRK_RKL_2);
   if (check_flag(&flag, "LSRKStepSetMethod", 1)) { return 1; }
 
   /* Open output stream for results, output comment line */

diff --git a/src/arkode/arkode_lsrkstep.c b/src/arkode/arkode_lsrkstep.c
@@ -408,13 +408,11 @@ int lsrkStep_Init(ARKodeMem ark_mem, int init_type)
     ark_mem->e_data    = ark_mem;
   }
 
-  /* Set default parameters for the default stepper */
+  /* Set the default stepper */
   if (ark_mem->step == lsrkStep_TakeStepRKC)
   {
-    step_mem->LSRKmethod   = ARKODE_LSRK_RKC_2;
-    step_mem->nfusedopvecs = 5;
-    step_mem->q = ark_mem->hadapt_mem->q = 2;
-    step_mem->p = ark_mem->hadapt_mem->p = 2;
+    retval = LSRKStepSetMethod(arkode_mem, ARKODE_LSRK_RKC_2);
+    if (retval != ARK_SUCCESS) { return (retval); }
   }
 
   /* Allocate ARK RHS vector memory, update storage requirements */
@@ -595,21 +593,14 @@ int lsrkStep_TakeStepRKC(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
   }
 
   /* determine the number of required stages */
-  int ss;
-  for (ss = 1; ss <= step_mem->stagemaxlimit; ss++)
+  int ss = SUNRceil(SUNRsqrt(onep54 * SUNRabs(ark_mem->h) * step_mem->sprad));
+  step_mem->reqstages = SUNMAX(ss, 2);
+  if (step_mem->reqstages == step_mem->stagemaxlimit)
   {
-    if (SUNSQR(ss) >= (onep54 * SUNRabs(ark_mem->h) * step_mem->sprad))
-    {
-      step_mem->reqstages = SUNMAX(ss, 2);
-      break;
-    }
-    if (ss == step_mem->stagemaxlimit)
-    {
-      hmax         = SUN_RCONST(0.95) * SUNSQR(ss) / (onep54 * step_mem->sprad);
-      ark_mem->eta = hmax / ark_mem->h;
-      *nflagPtr    = ARK_RETRY_STEP;
-      return (ARK_RETRY_STEP);
-    }
+    hmax         = SUN_RCONST(0.95) * SUNSQR(ss) / (onep54 * step_mem->sprad);
+    ark_mem->eta = hmax / ark_mem->h;
+    *nflagPtr    = ARK_RETRY_STEP;
+    return (ARK_RETRY_STEP);
   }
 
   step_mem->stagemax = SUNMAX(step_mem->reqstages, step_mem->stagemax);
@@ -798,17 +789,17 @@ int lsrkStep_TakeStepRKL(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
   }
 
   /* determine the number of required stages */
-  int ss;
-  for (ss = 1; ss <= step_mem->stagemaxlimit; ss++)
+  int ss = floor(SUNRsqrt(TWO * SUNRabs(ark_mem->h) * step_mem->sprad + THREE) - ONE);
+  for (ss; ss <= step_mem->stagemaxlimit; ss++)
   {
-    if ((SUNSQR(ss) + ss - 2) >= 2 * (SUNRabs(ark_mem->h) * step_mem->sprad))
+    if ((SUNSQR(ss) + ss - TWO) >= TWO * (SUNRabs(ark_mem->h) * step_mem->sprad))
     {
       step_mem->reqstages = SUNMAX(ss, 2);
       break;
     }
     if (ss == step_mem->stagemaxlimit)
     {
-      hmax = SUN_RCONST(0.95) * (SUNSQR(ss) + ss - 2) / (2.0 * step_mem->sprad);
+      hmax = SUN_RCONST(0.95) * (SUNSQR(ss) + ss - TWO) / (TWO * step_mem->sprad);
       ark_mem->eta = hmax / ark_mem->h;
       *nflagPtr    = ARK_RETRY_STEP;
       return (ARK_RETRY_STEP);
@@ -1295,7 +1286,7 @@ int lsrkStep_TakeStepSSP104(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt
                    ark_mem->fn, ark_mem->tempv1);
     }
   }
-  N_VLinearSum(SUN_RCONST(1.0 / 25.0), ark_mem->tempv2, SUN_RCONST(9.0 / 25.0),
+  N_VLinearSum(SUN_RCONST(1.0) / SUN_RCONST(25.0), ark_mem->tempv2, SUN_RCONST(9.0) / SUN_RCONST(25.0),
                ark_mem->ycur, ark_mem->tempv2);
   N_VLinearSum(SUN_RCONST(15.0), ark_mem->tempv2, SUN_RCONST(-5.0),
                ark_mem->ycur, ark_mem->ycur);