Steady-State Poisson with SUNDIALS KINSOL + Ginkgo

Solves the nonlinear equation $-u’’(x) - e^{u(x)} = 0$ on $[0,1]$ with $u(0) = u(1) = 0$ using SUNDIALS KINSOL (Newton iteration). Ginkgo provides the linear solver for each Newton step via SUNDIALS’ native Ginkgo interface.

Integration Overview

SUNDIALS v6.4+ provides native Ginkgo wrappers:

• sundials::ginkgo::Matrix<GkoCsr> wraps a Ginkgo CSR matrix as SUNMatrix

• sundials::ginkgo::LinearSolver<Factory> wraps a Ginkgo solver factory as SUNLinearSolver

These plug directly into KINSOL’s KINSetLinearSolver() — no custom glue code.

Step-by-Step

Problem Definition

// 1D discretization: -u''(x) - exp(u(x)) = 0 on [0,1], u(0)=u(1)=0.
// Central difference: (-u_{i-1} + 2*u_i - u_{i+1})/h^2 - exp(u_i) = 0.
struct PoissonData {
    int N;     // interior grid points
    double h;  // grid spacing
};

Nonlinear Residual

// Computes the nonlinear residual F(u) = 0.
int residual(N_Vector u, N_Vector fval, void* user_data)
{
    auto* data = static_cast<PoissonData*>(user_data);
    const int N = data->N;
    const double h2 = data->h * data->h;

    sunrealtype* u_data = N_VGetArrayPointer(u);
    sunrealtype* f_data = N_VGetArrayPointer(fval);

    for (int i = 0; i < N; i++) {
        double u_left = (i > 0) ? u_data[i - 1] : 0.0;
        double u_right = (i < N - 1) ? u_data[i + 1] : 0.0;
        double u_center = u_data[i];

        // -u'' - exp(u) via central differences.
        f_data[i] =
            (-u_left + 2.0 * u_center - u_right) / h2 - std::exp(u_center);
    }
    return 0;
}

Jacobian Assembly

The Jacobian callback fills the Ginkgo CSR matrix that lives inside the SUNMatrix wrapper:

// Computes the Jacobian J = dF/du as a tridiagonal CSR matrix.
int jacobian(N_Vector u, N_Vector fu, SUNMatrix J, void* user_data, N_Vector,
             N_Vector)
{
    auto* data = static_cast<PoissonData*>(user_data);
    const int N = data->N;
    const double h2 = data->h * data->h;
    sunrealtype* u_data = N_VGetArrayPointer(u);

    // Access the Ginkgo CSR matrix inside the SUNMatrix wrapper.
    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();

    // Fill tridiagonal Jacobian.
    int idx = 0;
    for (int i = 0; i < N; i++) {
        row_ptrs[i] = idx;

        if (i > 0) {
            col_idxs[idx] = i - 1;
            values[idx] = -1.0 / h2;
            idx++;
        }

        col_idxs[idx] = i;
        values[idx] = 2.0 / h2 - std::exp(u_data[i]);
        idx++;

        if (i < N - 1) {
            col_idxs[idx] = i + 1;
            values[idx] = -1.0 / h2;
            idx++;
        }
    }
    row_ptrs[N] = idx;

    return 0;
}

Ginkgo Setup

Create the Ginkgo executor, CSR matrix, and solver factory, then wrap them for SUNDIALS:

    // Create Ginkgo executor and objects for SUNDIALS.
    auto exec = gko::ReferenceExecutor::create();

    // Tridiagonal matrix has at most 3*N - 2 nonzeros.
    const int nnz = 3 * N - 2;

    using GkoCsr = gko::matrix::Csr<sunrealtype, sunindextype>;
    auto gko_jacobian =
        gko::share(GkoCsr::create(exec,
                                  gko::dim<2>{static_cast<gko::size_type>(N),
                                              static_cast<gko::size_type>(N)},
                                  static_cast<gko::size_type>(nnz)));

    // Wrap Ginkgo CSR matrix as SUNMatrix.
    sundials::ginkgo::Matrix<GkoCsr> sunA(gko_jacobian, sunctx);

    // Build Ginkgo GMRES solver factory and wrap as SUNLinearSolver.
    // 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-10)
                               .with_max_iters(static_cast<gko::uint64>(200)))
            .on(exec));

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

KINSOL Solve

    // Create KINSOL solver and attach Ginkgo linear solver.
    void* kin_mem = KINCreate(sunctx);
    KINInit(kin_mem, residual, u);
    KINSetUserData(kin_mem, &data);
    KINSetLinearSolver(kin_mem, sunLS, sunA);
    KINSetJacFn(kin_mem, jacobian);
    KINSetFuncNormTol(kin_mem, 1e-10);

    // Solve F(u) = 0 using Newton's method with line search.
    int flag = KINSol(kin_mem, u, KIN_LINESEARCH, scale, scale);
    if (flag >= 0)
        std::cout << "KINSOL converged with flag = " << flag << std::endl;
    else
        std::cout << "KINSOL failed with flag = " << flag << std::endl;

Building

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

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

Expected Output

KINSOL converged with flag = 0

Solution at selected points:
  u(0.00990099) = 0.00539008
  u(0.108911) = 0.0537863
  u(0.207921) = 0.091845
  u(0.306931) = 0.119166
  u(0.405941) = 0.135454
  u(0.50495) = 0.140527
  u(0.60396) = 0.134328
  u(0.70297) = 0.116927
  u(0.80198) = 0.0885165
  u(0.90099) = 0.0494044