revision

Author: maggul

doc/arkode/guide/source/Mathematics.rst

@@ -780,6 +780,7 @@ The SSPRK methods in ARKODE use the following Shu--Osher representation :cite:p:
    y_{n+1} &= z_s.
    :label: ARKODE_SSP
 
+The coefficients of Shu--Osher representation is not uniquely determined by Butcher table :cite:p:`SR:02`.
 In particular, the methods SSP(s,2), SSP(s,3), and SSP(10,4) implemented herein and presented in 
 :cite:p:`K:08` have "almost" all zero coefficients appearing in :math:`\alpha_{i,i-1}` and 
 :math:`\beta_{i,i-1}`. This feature facilitates their implementation in a low-storage manner. The 

doc/shared/sundials.bib

@@ -2421,3 +2421,14 @@ @article{SO:88
 publisher={Elsevier},
 doi={10.1016/0021-9991(88)90177-5},
 author={Shu, Chi-Wang and Osher, Stanley}}
+
+@article{SR:02, 
+title={A new class of optimal high-order strong-stability-preserving time discretization methods}, 
+journal={SIAM Journal on Numerical Analysis}, 
+volume={40},
+number={2}, 
+pages={469--491}, 
+year={2002}, 
+publisher={SIAM}, 
+doi={10.1137/S0036142901389025},
+author={Spiteri, Raymond J and Ruuth, Steven J}}

examples/arkode/C_serial/ark_analytic_lsrk_varjac.c

@@ -28,7 +28,7 @@
  * the value of alpha determines the radius of the interval in which
  * the stiffness parameter lies.
  *
- * The value of "lambda - alpha*cos((10 - t)/10*pi)" should
+ * The value of lambda - alpha*cos((10 - t)/10*pi) should
  * be negative to result in a well-posed ODE; for values with magnitude
  * larger than 100 the problem becomes quite stiff.
  *
@@ -97,8 +97,8 @@ int main(void)
   sunindextype NEQ   = 1;                  /* number of dependent vars. */
   sunrealtype reltol = SUN_RCONST(1.0e-8); /* tolerances */
   sunrealtype abstol = SUN_RCONST(1.0e-8);
-  sunrealtype lambda = SUN_RCONST(-1.0e+6); /* stiffness parameter 1*/
-  sunrealtype alpha  = SUN_RCONST(1.0e+2);  /* stiffness parameter 2*/
+  sunrealtype lambda = SUN_RCONST(-1.0e+6); /* stiffness parameter 1 */
+  sunrealtype alpha  = SUN_RCONST(1.0e+2);  /* stiffness parameter 2 */
   sunrealtype UserData[2];
   UserData[0] = lambda;
   UserData[1] = alpha;
@@ -237,9 +237,9 @@ static int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data)
 
   /* fill in the RHS function: "N_VGetArrayPointer" accesses the 0th entry of ydot */
   N_VGetArrayPointer(ydot)[0] =
-    (lambda - alpha * COS((10 - t) / 10 * ACOS(-1))) * u +
+    (lambda - alpha * COS((SUN_RCONST(10.0) - t) / SUN_RCONST(10.0) * ACOS(SUN_RCONST(-1.0)))) * u +
     SUN_RCONST(1.0) / (SUN_RCONST(1.0) + t * t) -
-    (lambda - alpha * COS((10 - t) / 10 * ACOS(-1))) * ATAN(t);
+    (lambda - alpha * COS((SUN_RCONST(10.0) - t) / SUN_RCONST(10.0) * ACOS(SUN_RCONST(-1.0)))) * ATAN(t);
 
   return 0; /* return with success */
 }
@@ -252,8 +252,8 @@ static int dom_eig(sunrealtype t, N_Vector y, N_Vector fn, sunrealtype* lambdaR,
   sunrealtype* rdata = (sunrealtype*)user_data; /* cast user_data to sunrealtype */
   sunrealtype lambda = rdata[0]; /* set shortcut for stiffness parameter 1 */
   sunrealtype alpha  = rdata[1]; /* set shortcut for stiffness parameter 2 */
-  *lambdaR           = (lambda - alpha * COS((10 - t) / 10 *
-                                             ACOS(-1))); /* access current solution value */
+  *lambdaR           = (lambda - alpha * COS((SUN_RCONST(10.0) - t) / SUN_RCONST(10.0) *
+                                             ACOS(SUN_RCONST(-1.0)))); /* access current solution value */
   *lambdaI           = SUN_RCONST(0.0);
 
   return 0; /* return with success */