Stiff ODE System with SUNDIALS CVODE + Ginkgo#

Solves the Robertson chemical kinetics problem — a classic stiff ODE test — using SUNDIALS CVODE (BDF method) with Ginkgo providing the linear solver for the implicit Newton steps.

$$\frac{du_1}{dt} = -0.04 u_1 + 10^4 u_2 u_3$$ $$\frac{du_2}{dt} = 0.04 u_1 - 10^4 u_2 u_3 - 3 \times 10^7 u_2^2$$ $$\frac{du_3}{dt} = 3 \times 10^7 u_2^2$$

Initial conditions: $u(0) = (1, 0, 0)$. Integrated to $t = 4 \times 10^7$.

Integration Overview#

Same native SUNDIALS-Ginkgo interface as the Poisson example:

sundials::ginkgo::Matrix<GkoCsr> wraps the Jacobian
sundials::ginkgo::LinearSolver<Factory> wraps the Ginkgo GMRES solver

CVODE calls Ginkgo’s GMRES at every implicit time step to solve the Newton system $(I - \gamma J) \delta = -F$.

Step-by-Step#

Right-Hand Side#

// Right-hand side: Robertson chemical kinetics.
int rhs(sunrealtype t, N_Vector u, N_Vector udot, void*)
{
    sunrealtype* y = N_VGetArrayPointer(u);
    sunrealtype* yd = N_VGetArrayPointer(udot);

    yd[0] = -0.04 * y[0] + 1.0e4 * y[1] * y[2];
    yd[1] = 0.04 * y[0] - 1.0e4 * y[1] * y[2] - 3.0e7 * y[1] * y[1];
    yd[2] = 3.0e7 * y[1] * y[1];

    return 0;
}

Jacobian#

The 3x3 analytical Jacobian is stored as CSR (9 entries):

// Analytical Jacobian df/du (3x3 dense stored as CSR).
int jac(sunrealtype t, N_Vector u, N_Vector fu, SUNMatrix J, void*, N_Vector,
        N_Vector, N_Vector)
{
    sunrealtype* y = N_VGetArrayPointer(u);

    using GkoCsr = gko::matrix::Csr<sunrealtype, sunindextype>;
    auto* gko_wrapper =
        static_cast<sundials::ginkgo::Matrix<GkoCsr>*>(J->content);
    auto gko_J = gko_wrapper->GkoMtx();

    auto* row_ptrs = gko_J->get_row_ptrs();
    auto* col_idxs = gko_J->get_col_idxs();
    auto* values = gko_J->get_values();

    // 3x3 dense matrix stored in CSR with 9 entries.
    // Row 0: df0/du0, df0/du1, df0/du2
    row_ptrs[0] = 0;
    col_idxs[0] = 0;
    values[0] = -0.04;
    col_idxs[1] = 1;
    values[1] = 1.0e4 * y[2];
    col_idxs[2] = 2;
    values[2] = 1.0e4 * y[1];

    // Row 1: df1/du0, df1/du1, df1/du2
    row_ptrs[1] = 3;
    col_idxs[3] = 0;
    values[3] = 0.04;
    col_idxs[4] = 1;
    values[4] = -1.0e4 * y[2] - 6.0e7 * y[1];
    col_idxs[5] = 2;
    values[5] = -1.0e4 * y[1];

    // Row 2: df2/du0, df2/du1, df2/du2
    row_ptrs[2] = 6;
    col_idxs[6] = 0;
    values[6] = 0.0;
    col_idxs[7] = 1;
    values[7] = 6.0e7 * y[1];
    col_idxs[8] = 2;
    values[8] = 0.0;

    row_ptrs[3] = 9;

    return 0;
}

Ginkgo Setup#

    auto exec = gko::ReferenceExecutor::create();

    // 3x3 dense stored as CSR (9 entries).
    using GkoCsr = gko::matrix::Csr<sunrealtype, sunindextype>;
    auto gko_jacobian = gko::share(GkoCsr::create(exec, gko::dim<2>{3, 3}, 9));

    sundials::ginkgo::Matrix<GkoCsr> sunA(gko_jacobian, sunctx);

    // Use Ginkgo GMRES as the linear solver for Newton steps.
    // Use SUNDIALS' DefaultStop criterion so that SUNDIALS can control
    // the tolerance and properly extract residual norms after each solve.
    using GkoSolver = gko::solver::Gmres<sunrealtype>;
    using DefaultStop = sundials::ginkgo::DefaultStop;
    auto solver_factory = gko::share(
        GkoSolver::build()
            .with_criteria(DefaultStop::build()
                               .with_reduction_factor(1e-14)
                               .with_max_iters(static_cast<gko::uint64>(50)))
            .on(exec));

    sundials::ginkgo::LinearSolver<GkoSolver, GkoCsr> sunLS(solver_factory,
                                                            sunctx);

Time Integration#

    // Create CVODE with BDF method (stiff problems).
    void* cvode_mem = CVodeCreate(CV_BDF, sunctx);
    CVodeInit(cvode_mem, rhs, t0, y);
    CVodeSStolerances(cvode_mem, 1e-8, 1e-10);

    // Attach Ginkgo linear solver.
    CVodeSetLinearSolver(cvode_mem, sunLS, sunA);
    CVodeSetJacFn(cvode_mem, jac);

    // Time integration loop.
    std::cout << std::scientific << std::setprecision(6);
    std::cout << "  t            u1          u2           u3" << std::endl;
    std::cout << "  " << t0 << "  " << y_data[0] << "  " << y_data[1] << "  "
              << y_data[2] << std::endl;

    sunrealtype t = t0;
    sunrealtype tout = 0.4;
    while (tout <= tf) {
        int flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
        if (flag < 0) {
            std::cerr << "CVODE error at t = " << t << std::endl;
            break;
        }

        std::cout << "  " << t << "  " << y_data[0] << "  " << y_data[1] << "  "
                  << y_data[2] << std::endl;

        tout *= 10.0;
    }

Building#

Requires SUNDIALS v6.4+ built with -DENABLE_GINKGO=ON.

cmake -S . -B build --preset default
cmake --build build
./build/ode_system

Expected Output#

  t            u1          u2           u3
000000e+00  1.000000e+00  0.000000e+00  0.000000e+00
000000e-01  9.851721e-01  3.386395e-05  1.479402e-02
000000e+00  9.055187e-01  2.240476e-05  9.445890e-02
000000e+01  7.158271e-01  9.185536e-06  2.841637e-01
000000e+02  4.505186e-01  3.222901e-06  5.494781e-01
000000e+03  1.832022e-01  8.942367e-07  8.167969e-01
000000e+04  3.898339e-02  1.621769e-07  9.610164e-01
000000e+05  4.938276e-03  1.984995e-08  9.950617e-01
000000e+06  5.168100e-04  2.068296e-09  9.994832e-01
000000e+07  5.203083e-05  2.081340e-10  9.999480e-01