Field conserving adaptive mesh refinement (AMR) scheme on massively parallel adaptive octree meshes

Adaptive mesh refinement (AMR) is widely used to efficiently resolve localized features in time-dependent partial differential equations (PDEs) by selectively refining and coarsening the mesh. However, in long-horizon simulations, repeated intergrid interpolations can introduce systematic drift in conserved quantities, especially for variational discretizations with continuous basis functions. While interpolation from parent-to-child during refinement in continuous Galerkin (CG) discretizations is naturally conservative, the standard injection-based child-to-parent coarsening interpolation is generally not. We propose a simple, scalable field-conserving coarsening operator for parallel, octree-based AMR. The method enforces discrete global conservation during coarsening by first computing field conserving coarse-element values at quadrature points and then recovering coarse nodal degrees of freedom via an $L^2$ projection (mass-matrix solve), which simultaneously controls the $L_2$ error. We evaluate the approach on mass-conserving phase-field models, including the Cahn–Hilliard and Cahn–Hilliard–Navier–Stokes systems, and compare against injection in terms of conservation error, solution quality, and computational cost.

💡 Research Summary

This paper addresses a subtle but critical source of error in continuous Galerkin (CG) finite‑element simulations that employ adaptive octree meshes. While refinement (parent‑to‑child) naturally conserves integral quantities because the new child nodes are obtained by interpolating the parent’s basis functions, the standard coarsening (child‑to‑parent) operation—typically implemented by simple injection—discards the degrees of freedom that lie on interior child nodes. Over many AMR cycles this leads to a systematic drift in globally conserved quantities such as mass, which is especially problematic for long‑time simulations of phase‑field models (Cahn–Hilliard, Cahn–Hilliard–Navier–Stokes) where exact mass conservation is a physical requirement.

The authors propose a field‑conserving coarsening operator that enforces discrete global conservation on a per‑element basis and then reconstructs the coarse‑mesh nodal values via an L² projection. The algorithm proceeds as follows:

1. Quadrature‑point conservation – For each coarse (parent) element, evaluate the fine‑mesh field at the quadrature points of all its child elements. The sum of the weighted fine‑mesh values must equal the weighted sum of the unknown coarse‑mesh values at the parent’s quadrature points (Equation 5). This guarantees that the integral of the field over the parent element is exactly preserved to machine precision.

2. L² projection (mass‑matrix solve) – The unknown coarse‑mesh nodal values are obtained by minimizing the L² error between the fine‑mesh representation (restricted to the parent domain) and the coarse‑mesh representation. This leads to a local linear system (M U = b), where (M) is the element‑wise mass matrix and (b) contains the fine‑mesh contributions projected onto the coarse basis functions. By choosing Gauss‑Legendre quadrature points that coincide with the Lagrange nodes of the basis, the mass matrix becomes diagonal ((M_{ij}=w_i\delta_{ij})), allowing the projection to be performed with a simple weighted average rather than a full linear solve.

3. Extension to higher order and dimensions – The method works for any polynomial order (p\ge 1) because the same projection framework applies; only the size of the local system grows as ((p+1)^d). In multiple dimensions the 1‑D restriction matrix is extended via tensor products (e.g., (R_{2D}=R_{1D}\otimes R_{1D})). When more quadrature points than the minimum are required, the mass matrix is no longer diagonal, but it remains a small, element‑local system that can be pre‑factorized.

The key advantage of this approach is that it requires only element‑local operations; no global communication or global constraint solving is needed. Consequently, the method scales almost perfectly on massively parallel machines. The authors demonstrate scalability up to 65 k cores, with the coarsening step consuming only 4–5 % of total runtime (compared to 3 % for the injection method).

Numerical experiments focus on two benchmark phase‑field problems:

• Cahn–Hilliard – 2‑D and 3‑D simulations of spinodal decomposition with periodic boundaries. Over 10 000 AMR cycles, the injection method accumulates a mass error of order (10^{-3}), whereas the proposed method maintains mass error at machine precision ((<10^{-12})). The free‑energy decay curves are indistinguishable, confirming that solution quality is unaffected.

• Cahn–Hilliard–Navier–Stokes – A two‑phase incompressible flow problem where the interface dynamics are coupled to the Navier–Stokes equations. Again, mass drift is eliminated by the new coarsening operator, while velocity fields, pressure, and interface shapes remain virtually identical to those obtained with injection.

Performance measurements show a modest overhead of roughly 5–10 % in the coarsening stage due to the extra quadrature evaluations and the local projection, but this is outweighed by the elimination of the systematic mass drift, which would otherwise require smaller time steps or more frequent re‑meshings to control.

Discussion highlights that the problem is intrinsic to CG discretizations: because nodal continuity couples neighboring elements, any local modification of nodal values to enforce conservation inevitably influences the solution globally. Therefore, a purely local “mass redistribution” strategy (as used in finite‑volume or DG methods) is insufficient. The proposed element‑wise L² projection provides a mathematically rigorous and computationally cheap way to reconcile local conservation with the global continuity constraints of CG.

Conclusion – The paper delivers a practical, scalable, and mathematically sound solution to the long‑standing issue of non‑conservative coarsening in CG‑based octree AMR. By guaranteeing per‑element integral conservation and reconstructing coarse nodal values through an L² projection, the method eliminates global mass drift without compromising accuracy or incurring prohibitive computational cost. Its applicability to arbitrary polynomial orders, higher‑dimensional meshes, and large‑scale parallel environments makes it a valuable addition to any CG‑AMR framework, particularly for multiphysics simulations where long‑time mass conservation is non‑negotiable.