add LSRKStep logging test

test/unit_tests/logging/CMakeLists.txt

@@ -16,7 +16,7 @@
 set(unit_tests)
 
 if(BUILD_ARKODE)
-  list(APPEND unit_tests "test_logging_arkode_erkstep.cpp\;")
+  # ARKStep
   list(APPEND unit_tests "test_logging_arkode_arkstep.cpp\;0") # ERK
   list(APPEND unit_tests "test_logging_arkode_arkstep.cpp\;1 1 1") # DIRK Newton
                                                                    # + Dense
@@ -30,6 +30,11 @@ if(BUILD_ARKODE)
                                                                    # + GMRES
   list(APPEND unit_tests "test_logging_arkode_arkstep.cpp\;2 0") # ImEx
                                                                  # Fixed-point
+  # ERKStep
+  list(APPEND unit_tests "test_logging_arkode_erkstep.cpp\;")
+  # LSRKStep
+  list(APPEND unit_tests "test_logging_arkode_lsrkstep.cpp\;")
+  # MRIStep
   list(APPEND unit_tests "test_logging_arkode_mristep.cpp\;0") # Ex-MRI
   list(APPEND unit_tests "test_logging_arkode_mristep.cpp\;1 1 1") # Im-MRI
                                                                    # Newton +

test/unit_tests/logging/test_logging_arkode_lsrkstep.cpp

@@ -0,0 +1,144 @@
+/* -----------------------------------------------------------------------------
+ * Programmer(s): David J. Gardner @ LLNL
+ * -----------------------------------------------------------------------------
+ * SUNDIALS Copyright Start
+ * Copyright (c) 2002-2024, 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
+ * -----------------------------------------------------------------------------
+ * Test logging output in LSRKStep
+ * ---------------------------------------------------------------------------*/
+
+#include <cmath>
+#include <cstdio>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+
+// Include desired integrators, vectors, linear solvers, and nonlinear solvers
+#include "arkode/arkode_lsrkstep.h"
+#include "nvector/nvector_serial.h"
+#include "sundials/sundials_logger.h"
+
+#include "problems/prv.hpp"
+#include "sundials/sundials_nvector.h"
+#include "utilities/check_return.hpp"
+
+using namespace std;
+using namespace problems::prv;
+
+int main(int argc, char* argv[])
+{
+  cout << "Start LSRKStep Logging test" << endl;
+
+  // SUNDIALS context object for this simulation
+  sundials::Context sunctx;
+
+  // Use RKC (0) or RKL (1)
+  int method_type = 0;
+  if (argc > 1) { method_type = stoi(argv[1]); }
+
+  // Ensure logging output goes to stdout
+  SUNLogger logger;
+  int flag = SUNContext_GetLogger(sunctx, &logger);
+  if (check_flag(flag, "SUNContext_GetLogger")) { return 1; }
+
+  SUNLogger_SetErrorFilename(logger, "stdout");
+  SUNLogger_SetWarningFilename(logger, "stdout");
+  SUNLogger_SetInfoFilename(logger, "stdout");
+  SUNLogger_SetDebugFilename(logger, "stdout");
+
+  // Create initial condition
+  N_Vector y = N_VNew_Serial(1, sunctx);
+  if (check_ptr(y, "N_VNew_Serial")) { return 1; }
+  N_VConst(true_solution(zero), y);
+
+  // Create LSRKStep memory structure
+  void* arkode_mem = LSRKStepCreateSTS(ode_rhs, zero, y, sunctx);
+  if (check_ptr(arkode_mem, "LSRKStepCreate")) { return 1; }
+
+  // Select method
+  if (method_type == 0)
+  {
+    flag = LSRKStepSetSTSMethodByName(arkode_mem, "ARKODE_LSRK_RKC_2");
+    if (check_flag(flag, "LSRKStepSetSTSMethodByName")) { return 1; }
+  }
+  else if (method_type == 1)
+  {
+    flag = LSRKStepSetSTSMethodByName(arkode_mem, "ARKODE_LSRK_RKL_2");
+    if (check_flag(flag, "LSRKStepSetSTSMethodByName")) { return 1; }
+  }
+  else
+  {
+    cerr << "Invalid method option\n";
+    return 1;
+  }
+
+  flag = ARKodeSetUserData(arkode_mem, &prv_data);
+  if (check_flag(flag, "ARKodeSetUserData")) { return 1; }
+
+  // Relative and absolute tolerances
+  const sunrealtype rtol = SUN_RCONST(1.0e-6);
+  const sunrealtype atol = SUN_RCONST(1.0e-10);
+
+  flag = ARKodeSStolerances(arkode_mem, rtol, atol);
+  if (check_flag(flag, "ARKodeSStolerances")) { return 1; }
+
+  // Specify dominant eigenvalue function
+  flag = LSRKStepSetDomEigFn(arkode_mem, ode_dom_eig);
+  if (check_flag(flag, "LSRKStepSetDomEigFn")) { return 1; }
+
+  // Initial time and fist output time
+  const sunrealtype dtout = one; // output interval
+  const int nout          = 3;   // number of outputs
+  sunrealtype tret        = zero;
+  sunrealtype tout        = tret + dtout;
+
+  // Output initial contion
+  cout << scientific;
+  cout << setprecision(numeric_limits<sunrealtype>::digits10);
+  cout << "           t              ";
+  cout << "          y              ";
+  cout << "        y err            ";
+  for (int i = 0; i < 7; i++) { cout << "--------------"; }
+  cout << endl;
+
+  sunrealtype* y_data = N_VGetArrayPointer(y);
+
+  cout << setw(22) << tret << setw(25) << y_data[0] << setw(25)
+       << abs(y_data[0] - true_solution(tret)) << endl;
+
+  // Advance in time
+  for (int i = 0; i < nout; i++)
+  {
+    flag = ARKodeEvolve(arkode_mem, tout, y, &tret, ARK_ONE_STEP);
+    if (check_flag(flag, "ARKodeEvolve")) { return 1; }
+
+    cout << setw(22) << tret << setw(25) << y_data[0] << setw(25)
+         << abs(y_data[0] - true_solution(tret)) << endl;
+
+    // update output time
+    tout += dtout;
+  }
+  for (int i = 0; i < 7; i++) { cout << "--------------"; }
+  cout << endl;
+
+  // Print some final statistics
+  flag = ARKodePrintAllStats(arkode_mem, stdout, SUN_OUTPUTFORMAT_TABLE);
+  if (check_flag(flag, "ARKodePrintAllStats")) { return 1; }
+
+  // Clean up and return with successful completion
+  N_VDestroy(y);
+  ARKodeFree(&arkode_mem);
+
+  cout << "End LSRKStep Logging test" << endl;
+
+  return 0;
+}
+
+/*---- end of file ----*/

test/unit_tests/problems/kpr.hpp

@@ -11,7 +11,7 @@
  * SPDX-License-Identifier: BSD-3-Clause
  * SUNDIALS Copyright End
  * -----------------------------------------------------------------------------
- * Kvaerno-Prothero-Robinson ODE test problem:
+ * Kvaerno-Prothero-Robinson (KPR) ODE test problem:
  *
  *   [u]' = [ a  b ] [ (-1 + u^2 - r(t)) / (2u) ] + [ r'(t) / (2u) ]
  *   [v]    [ c  d ] [ (-2 + v^2 - s(t)) / (2v) ]   [ s'(t) / (2v) ]

test/unit_tests/problems/prv.hpp

@@ -0,0 +1,101 @@
+/* -----------------------------------------------------------------------------
+ * Programmer(s): David J. Gardner @ LLNL
+ * -----------------------------------------------------------------------------
+ * SUNDIALS Copyright Start
+ * Copyright (c) 2002-2024, 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
+ * -----------------------------------------------------------------------------
+ * Prothero-Robinson ODE test problem with a time-varying coefficient (PRV):
+ *
+ *   y' = L(t) (y - phi(t)) + phi'(t),
+ *
+ * where phi(t) = atan(t), phi'(t) = 1 / (1 + t^2), and the coefficient is
+ * given by
+ *
+ *   L(t) = lambda - alpha * cos((10 - t) / 10 * pi).
+ *
+ * This problem has analytical solution y(t) = atan(t).
+ *
+ * The stiffness of the problem depends on the value of L(t) where the lambda
+ * determines the center of the parameter and alpha the radius of the interval
+ * in which the stiffness parameter lies. For a well-posed problem, the values
+ * should be chosen such that L(t) is negative.
+ * ---------------------------------------------------------------------------*/
+
+#ifndef PRV_
+#define PRV_
+
+#include <cmath>
+#include <sundials/sundials_core.hpp>
+
+namespace problems {
+namespace prv {
+
+// Problem constants
+static const sunrealtype pi = std::acos(SUN_RCONST(-1.0));
+
+constexpr sunrealtype zero = SUN_RCONST(0.0);
+constexpr sunrealtype one  = SUN_RCONST(1.0);
+constexpr sunrealtype ten  = SUN_RCONST(10.0);
+
+// lambda and alpha
+sunrealtype prv_data[2] = {SUN_RCONST(-1000.0), SUN_RCONST(10.0)};
+
+// -----------------------------------------------------------------------------
+// Helper functions
+// -----------------------------------------------------------------------------
+
+// Compute L(t)
+inline sunrealtype l_coef(sunrealtype t, sunrealtype c[2])
+{
+  return c[0] - c[1] * std::cos((ten - t) / ten * pi);
+}
+
+// Compute phi(t)
+inline sunrealtype phi(sunrealtype t) { return std::atan(t); }
+
+// Compute phi'(t)
+inline sunrealtype phi_prime(sunrealtype t) { return one / (one + t * t); }
+
+// Compute the true solution
+inline sunrealtype true_solution(sunrealtype t)
+{
+  return phi(t);
+}
+
+// ODE RHS function
+inline int ode_rhs(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data)
+{
+  sunrealtype* u_data  = static_cast<sunrealtype*>(user_data);
+  sunrealtype* y_data  = N_VGetArrayPointer(y);
+  sunrealtype* yd_data = N_VGetArrayPointer(ydot);
+
+  yd_data[0] = l_coef(t, u_data) * (y_data[0] - phi(t)) + phi_prime(t);
+
+  return 0;
+}
+
+// Dominant eigenvalue function
+inline int ode_dom_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)
+{
+  sunrealtype* u_data = static_cast<sunrealtype*>(user_data);
+
+  *lambdaR = l_coef(t, u_data);
+  *lambdaI = zero;
+
+  return 0;
+}
+
+} // namespace prv
+} // namespace problems
+
+#endif