High-Order Diffusion with MFEM + Ginkgo

Solves $-\nabla \cdot (\kappa \nabla u) = 1$ with a spatially varying coefficient $\kappa(x) = 1 + x_0^2$ using quartic ($p=4$) finite elements. Shows two integration approaches: manual CSR wrapping and MFEM’s built-in Ginkgo interface.

Integration Overview

High-order elements produce denser matrix rows (more non-zeros per row), making solver and preconditioner choice more impactful. This example demonstrates both integration paths:

1. Manual CSR wrapping — full control over Ginkgo’s solver factory, preconditioner, and stopping criteria

2. MFEM built-in (requires MFEM_USE_GINKGO=ON) — simpler but less configurable

Step-by-Step

Variable Coefficient

// A spatially varying diffusion coefficient kappa(x) = 1 + x^2.
double kappa_func(const Vector& x) { return 1.0 + x(0) * x(0); }

Assembly with High-Order Elements

The only change from linear elements is order = 4 in the H1_FECollection. The assembled matrix is larger and denser.

    // Essential boundary DOFs.
    Array<int> ess_tdof_list;
    if (mesh.bdr_attributes.Size()) {
        Array<int> ess_bdr(mesh.bdr_attributes.Max());
        ess_bdr = 1;
        fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
    }

    // RHS: f = 1.
    LinearForm b(&fespace);
    ConstantCoefficient one(1.0);
    b.AddDomainIntegrator(new DomainLFIntegrator(one));
    b.Assemble();

    GridFunction x(&fespace);
    x = 0.0;

    // Stiffness matrix with variable diffusion coefficient.
    FunctionCoefficient kappa(kappa_func);
    BilinearForm a(&fespace);
    a.AddDomainIntegrator(new DiffusionIntegrator(kappa));
    a.Assemble();

    SparseMatrix A;
    Vector B, X;
    a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
    std::cout << "System: " << A.Height() << " x " << A.Width()
              << ", nnz: " << A.NumNonZeroElems()
              << " (avg nnz/row: " << A.NumNonZeroElems() / A.Height() << ")"
              << std::endl;

Approach 1: Manual Ginkgo CSR Wrapping

// Approach 1: Manual CSR wrapping — full control over Ginkgo
// solver/preconditioner.
void solve_manual(SparseMatrix& A, Vector& B, Vector& X)
{
    A.Finalize();
    const int n = A.Height();
    const int nnz = A.NumNonZeroElems();

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

    auto gko_A = gko::share(gko::matrix::Csr<double, int>::create(
        exec,
        gko::dim<2>(static_cast<gko::size_type>(n),
                    static_cast<gko::size_type>(n)),
        gko::array<double>::view(exec, nnz, A.GetData()),
        gko::array<int>::view(exec, nnz, A.GetJ()),
        gko::array<int>::view(exec, n + 1, A.GetI())));

    auto gko_b = gko::matrix::Dense<double>::create(
        exec, gko::dim<2>(static_cast<gko::size_type>(n), 1),
        gko::array<double>::view(exec, n, B.GetData()), 1);
    auto gko_x = gko::matrix::Dense<double>::create(
        exec, gko::dim<2>(static_cast<gko::size_type>(n), 1),
        gko::array<double>::view(exec, n, X.GetData()), 1);

    auto solver =
        gko::solver::Cg<double>::build()
            .with_preconditioner(gko::preconditioner::Jacobi<double>::build())
            .with_criteria(
                gko::stop::Iteration::build().with_max_iters(1000u),
                gko::stop::ResidualNorm<double>::build().with_reduction_factor(
                    1e-12))
            .on(exec)
            ->generate(gko_A);

    solver->apply(gko_b, gko_x);
}

Approach 2: MFEM Built-in Ginkgo Interface

// Approach 2: MFEM's built-in Ginkgo wrappers — simpler, less control.
// Requires MFEM built with MFEM_USE_GINKGO=ON.
void solve_builtin(SparseMatrix& A, Vector& B, Vector& X)
{
    Ginkgo::GinkgoExecutor exec(Ginkgo::GinkgoExecutor::OMP);
    Ginkgo::CGSolver solver(exec);
    solver.SetRelTol(1e-12);
    solver.SetMaxIter(1000);
    solver.SetOperator(A);
    solver.Mult(B, X);
}

Building

cmake -S . -B build --preset default
cmake --build build
./build/high_order -o 4 -r 3

Expected Output

Options used:
   --order 4
   --refine 2
Order: 4
Elements: 1024
DOFs: 16641
System: 16641 x 16641, nnz: 591361 (avg nnz/row: 35)

Solving with manual Ginkgo CSR wrapping...
Converged.
Output written to solution.gf and mesh.mesh