Adaptive Refinement with deal.II + Ginkgo#

Solves the Poisson equation $-\Delta u = 1$ on a disk with adaptive mesh refinement (AMR). After each refinement cycle, the mesh changes, the system is re-assembled, and a new Ginkgo solver is generated from the updated matrix.

Based on deal.II tutorial step-6.

Integration Overview#

The key difference from the basic Poisson example is that AMR changes the matrix dimensions and sparsity pattern every cycle. The Ginkgo solver (and its preconditioner) must be re-created from scratch after each refinement. This example shows the pattern for solver re-generation.

Step-by-Step#

CSR Extraction Helper#

A reusable function that extracts CSR arrays from a deal.II SparseMatrix:

// Extracts CSR arrays from a deal.II SparseMatrix for Ginkgo consumption.
void extract_csr(const SparseMatrix<double>& matrix, std::vector<int>& row_ptrs,
                 std::vector<int>& col_idxs, std::vector<double>& values)
{
    const auto n = matrix.m();
    const auto nnz = matrix.n_nonzero_elements();

    row_ptrs.resize(n + 1);
    col_idxs.resize(nnz);
    values.resize(nnz);

    row_ptrs[0] = 0;
    for (unsigned int row = 0; row < n; ++row) {
        int idx = row_ptrs[row];
        for (auto it = matrix.begin(row); it != matrix.end(row); ++it) {
            col_idxs[idx] = static_cast<int>(it->column());
            values[idx] = it->value();
            ++idx;
        }
        row_ptrs[row + 1] = idx;
    }
}

Solve Function with Ginkgo#

Each cycle calls this function, which creates a fresh Ginkgo solver. After solving, constraints.distribute() enforces hanging node constraints.

// Builds a new Ginkgo solver from the current system matrix and solves.
// A new solver must be generated each cycle because the matrix changes
// after refinement (different size, different sparsity pattern).
void solve_with_ginkgo(const SparseMatrix<double>& system_matrix,
                       const Vector<double>& system_rhs,
                       Vector<double>& solution,
                       const AffineConstraints<double>& constraints)
{
    auto exec = gko::ReferenceExecutor::create();

    std::vector<int> row_ptrs;
    std::vector<int> col_idxs;
    std::vector<double> values;
    extract_csr(system_matrix, row_ptrs, col_idxs, values);

    const auto n = system_matrix.m();
    const auto nnz = system_matrix.n_nonzero_elements();

    auto mtx = gko::share(gko::matrix::Csr<double, int>::create(
        exec, gko::dim<2>(n, n),
        gko::array<double>::view(exec, nnz, values.data()),
        gko::array<int>::view(exec, nnz, col_idxs.data()),
        gko::array<int>::view(exec, n + 1, row_ptrs.data())));

    auto b = gko::matrix::Dense<double>::create(
        exec, gko::dim<2>(n, 1),
        gko::array<double>::view(exec, n,
                                 const_cast<double*>(system_rhs.begin())),
        1);
    auto x = gko::matrix::Dense<double>::create(
        exec, gko::dim<2>(n, 1),
        gko::array<double>::view(exec, n, solution.begin()), 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(mtx);

    solver->apply(b, x);

    // Enforce hanging node constraints on the solution.
    constraints.distribute(solution);
}

Building#

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

Expected Output#

Cycle 0:
  Cells: 20  DOFs: 89
  Solved.
Cycle 1:
  Cells: 44  DOFs: 209
  Solved.
...
Cycle 7:
  Cells: 2513  DOFs: 12725
  Solved.

Each cycle produces a solution-N.vtk file for visualization.