Implement AgFEM on Octree Distributed Meshes

[ericneiva]

LAST UPDATE APR 17 2026

Hi all,

Here some minimal notes about this development. I might extend it next week.

Next steps (for now):

• @amartinhuertas confirms the ghost layer can be extended beyond first layer
  • https://github.com/cburstedde/p4est/blob/master/src/p4est_ghost.h#L408
• @amartinhuertas creates a model M in GridapP4est.jl with >1 ghost layers [WIP here]
• @JordiManyer checks that M runs through Gridap until assembly
  • For the current set of tests in GridapP4est, the code seems to run consistently with M from the beginning to end, including FE assembly and solution of the linear system.
  • We can now create cartesian models with arbitrary number of ghost layers using the following branch of GridapDistributed.
• @JordiManyer and @ericneiva adapt the distributed cell aggregation scheme to the scenario > 1 ghost layers
  • Working on this branch checking out expand-ghost from GridapDistributed.jl
  • @ericneiva designs a test on uniform meshes where a "compliant" path can only be found if minimum distance to centroid criterion is communicated: 2D test here
  • @ericneiva updates the sequential cell aggregation such that it provides the "objective function" to the distributed counterpart. See 1866a38.
  • @ericneiva Adds warning when the root cell of a last-layer ghost cell is an owned cell
  • @ericneiva Brief performance comparison w.r.t. old algorithm
    • @ericneiva Updates results with minimum wall times
    • @JordiManyer Optimise old aggregation to recover weak scalability
    • @ericneiva reruns performance comparison of dist-aggregation on uniform in MNV
    • (to be decided) Performance comparison including times to set up dofs
• @JordiManyer reworks the distribution of DoF ownership on uniform meshes according to his notes
  • [Pre-requisite] @ericneiva deals with the local distribution of DoFs on the patch-conforming uniform mesh (i.e. the sequential part of the task). See 7fd7a6c.
• @amartinhuertas implements support for the creation of a FESpace out of a Triangulation (currently in a GridapP4est.jl branch: https://github.com/gridap/GridapP4est.jl/tree/fe_space_on_triangulation)
• @ericneiva updates closest root cell criterion (relative distance) of Definition 2.2 in https://arxiv.org/pdf/2006.05373
• @janmodderman tests sequential AgFEM on octree meshes
  • https://github.com/janmodderman/GridapEmbedded.jl/tree/hydro_code
  • Convergence tests on Fichera corner PB or similar
    • h-adaptive refinements driven with a suitable error estimator
  • First with CutFEM, then with AgFEM
• distributed h-AgFEM
  • @amartinhuertas Updates GridapP4est to the latest Gridap release
  • @amartinhuertas Update set up of conforming space in GridapP4est.jl with rework of FESpaceWithLinearConstraints
  • @ericneiva provides h-AgFEM distribution of dofs on the "owned" dofs (raw constraint arrays)
    • Rewrite assuming raw constraint arrays as input
      • Attempted in 007d987 need LinearConstraintsMap which relies on FESpaceWithLinearConstraints data layout
      • Documented pitfall in uncoupling from internal data layout of FESpaceWithLinearConstraints
    • Rewrite using the constraints branch
    • Performance optimisations
      • @ericneiva use local dictionaries to compress repeated constraining dofs of hanging dofs
      • @ericneiva checks that OctreeDistributedDiscreteModels "have affine maps"
      • Precompute root-reindex cell maps on a bounding box of the aggregates
  • @JordiManyer distributed dof setup of h-AgFEM extending the algs developed for the uniform case
  • NOTE that, strictly, we only need to complete/make-consistent across processors the ill-posed free dof constraints, because the ill-posed hanging dofs can be locally resolved once the ill-posed free dofs are consistently resolved
  • TO EXPLORE An abstract framework for composing linear constraints
• @oriolcg Complete weighted-aggregated SBM formulation with the DG-like terms on cut/cut and cut/interior facets
• @janmodderman update Fichera code and implement the Li and Bettess refinement criterion

I will keep this issue updated.

FYI @oriolcg @santiagobadia

Let me know if I miss something important or feel free to comment it below.

Thanks :)

[amartinhuertas]

Related to #106

[amartinhuertas]

@ericneiva ... I have updated the steps above with some work in progress I have been doing.

One step that we need and we are currently missing is the ability to determine how many layers of ghost cells we will need. I guess that this should be relatively easy to infer as part of the distributed algorithm that builds the aggregates in parallel, which I believe already exists, right? Can you take care of that part? Thanks!

[JordiManyer]

I think it should be set by the user at the beginning, and we limit the aggregation to local cells. I dont know if we want the number of ghost layers to be dynamic. Can it be dynamic from p4est point of view?

I would instead create some kind of indicator within the aggregation algorithm that outputs a warning if the overlapping is not big enough.

[amartinhuertas]

I think it should be set by the user at the beginning, and we limit the aggregation to local cells.

Not sure if I understand the second part of the sentence. I guess that you mean that there should be enough ghost layers so that the aggregates can be built locally. But in the case that there are not enough layers, then the code crashes or something like that?

I dont know if we want the number of ghost layers to be dynamic. Can it be dynamic from p4est point of view?

Do you mean that the number of ghost layers grows as the aggregation algorithm moves forward? I think this could be a possibility, but an expensive one, as I believe we would need to reconstruct the model from scratch each time that we need to add a ghost layer (or at least this would be the naive way of implementing this, right now I am not sure if the way we have partition the functionality into functions and data structures would allow as to grow the ghost layers without recomputing the model from scratch). From the p4est point of view is doable, but not sure form the GridapP4est' data structures point of view.

I would instead create some kind of indicator within the aggregation algorithm that outputs a warning if the overlapping is not big enough.

I would need to refresh the algorithm the builds the aggregates with distributed data structures to be more certain about an answer to either this possibility or any of the other ones. Perhaps @ericneiva can also give us here some insights?

[JordiManyer]

Do you mean that the number of ghost layers grows as the aggregation algorithm moves forward? I think this could be a possibility, but an expensive one

Then I would definitely do what I said above.

But in the case that there are not enough layers, then the code crashes or something like that?

I don't know if we need to crash, a warning is probably enough. At the end of the day, if the aggregation is not good enough the code will still work, you will just get crappy convergence rates. So we should allow the user to under-aggregate if they want to.

But as I said above, we should have some sort of indicator that outputs a warning whenever the algorithm feels there is not enough ghost layers to aggregate cells in a robust way.

[ericneiva]

Dear @amartinhuertas and @JordiManyer,

Thanks a lot for your work and comments.

how many layers of ghost cells we will need

I agree with @JordiManyer, this should be a user-prescribed value.

Nevertheless, I would say the iter value here is an upper-bound of the amount of ghost layers needed.

This leads me to the question, @amartinhuertas:

If the user sets, e.g., num_ghost_layers = 5 and when the loop breaks (i.e., all_aggregated) iter = 2, would you still avoid reconstructing the model (with num_ghost_layers = 2)?

indicator that outputs a warning whenever the algorithm feels there is not enough ghost layers

Yes, we must add this, recommending to consider first (adaptive) refinements and, if difficult, to increase the number of ghost layers.

It should be relatively easy to add the indicator in the cell aggregation scheme. For instance, we could query if the global cell id of the root cell is outside the global indices of the (owned+ghost) local cells. I'll see if there is a better way.

if the aggregation is not good enough the code will still work, you will just get crappy convergence rates.

Or worse, I would not be surprised if the solver fails.

Can you take care of that part?

Of course, @amartinhuertas, happy to take care of adapting the distributed cell aggregation scheme to reflect that num_ghost_layers >= 1 + add the indicator. Tasklist updated.

[amartinhuertas]

Hi @ericneiva,

Let me answer your questions below, and also ask you further questions :-)

Nevertheless, I would say the iter value here is an upper-bound of the amount of ghost layers needed.

A couple of follow-up questions: (1) does this algorithm work verbatim for non-conforming meshes? (I guess so, but just double checking). (2) As far as I can see, at each iteration of this algorithm there is local work followed by nearest neighbour exchanges (assuming a single layer of ghost cells). Do we actually want to do this when we have more than one layer of ghost cells? or otherwise build the aggregates using just local information and avoid any communication?

If the user sets, e.g., num_ghost_layers = 5 and when the loop breaks (i.e., all_aggregated) iter = 2, would you still avoid reconstructing the model (with num_ghost_layers = 2)?

That's a good question. I think that, IDEALLY, once we have built the aggregates, we should try to keep the number of remote cells stored locally as small as possible, i.e., to the minimum required to resolve the constraints locally. This might imply to reduce the number of ghost cells even below those that are required whenever num_ghost_layers=2, otherwise, we would be bringing with each nearest neighbour comm (e.g., during a nonlinear iteration) much more information than we actually need. How easy that would be to implement, not sure at the moment. Perhaps as a first approach we can reconstruct the model with num_ghost_layers=2 and see whether that would be acceptable. Overall, having multiple layers of ghost cells concerns me that we may end up with an unbearable amount of extra space and extra communication overhead.

Yes, we must add this, recommending to consider first (adaptive) refinements and, if difficult, to increase the number of ghost layers.

Not sure what do you mean here by "adaptive refinements".

[JordiManyer]

build the aggregates using just local information and avoid any communication

isn't this the point of expanding the ghost layers? I think we should just limit ourselves to local information. If more layers are needed, we just tell the user.

otherwise, we would be bringing with each nearest neighbour comm (e.g., during a nonlinear iteration) much more information than we actually need

I'm not that sure about that. We can always use Triangulations and FESpaces that are limited to the cells touched by the aggregates. I'm not sure we communicate that much using the model itself.

[amartinhuertas]

We can always use Triangulations and FESpaces that are limited to the cells touched by the aggregates.

That would be ideal. However, at the moment, GridapP4est does not support the creation of a FESpace out of a Triangulation. You need to pass the model to the FESpace constructor. I guess it could be possible to extend it such that you can pass a Triangulation with the minimum possible set of ghost cells, but this is not possible at the moment.

[amartinhuertas]

isn't this the point of expanding the ghost layers? I think we should just limit ourselves to local information. If more layers are needed, we just tell the user.

When we talked in our last meeting, I had in my mind that having additional ghost layers was about having a local data structure where we can run the current algorithm that builds the hanging DoF constraints locally verbatim. I was not actually thinking in the construction of the aggregates themselves. But I think it also may make sense to leverage the extra ghost cells to carry the construction of the aggregates in a purely local fashion without communication (beyond any that was required to bring the extra ghost layers).

[JordiManyer]

I was not actually thinking in the construction of the aggregates themselves. But I think it also may make sense to leverage the extra ghost cells to carry the construction of the aggregates in a purely local fashion without communication (beyond any that was required to bring the extra ghost layers).

Yes, this is what was on my mind when I proposed it.

[ericneiva]

Dear @amartinhuertas and @JordiManyer,

build the aggregates using just local information and avoid any communication

The aggregation scheme constructs facet-connected cell paths (via cut or interior facets) from root cells to cut cells. The paths are incremented at every iteration of the scheme, until all cut cells have been visited.

A cell path is incremented by searching for non-visited cut cells among the facet neighbours of visited cells. If we build the aggregates with just local information, I would say the facets on the ghost boundary will not provide the full neighbour list.

GridapEmbedded.jl/src/Distributed/DistributedAgFEM.jl

Line 258 in b596a9a

neigh_cells = getindex!(c2,face_to_cells,face)

All this to say: if we do this, we can no longer guarantee that the scheme is invariant to the partition of the model. I would say it will work anyway, but I am not sure if this is acceptable or if there are other consequences.

Let me know what you think.

(1) does this algorithm work verbatim for non-conforming meshes?

If this

GridapEmbedded.jl/src/Distributed/DistributedAgFEM.jl

Line 258 in b596a9a

neigh_cells = getindex!(c2,face_to_cells,face)

provides both the conforming and non-conforming cells around, yes.

Not sure what do you mean here by "adaptive refinements".

I mean refining the mesh either uniformly or with h-adaptive refinements. If you cannot complete the aggregation in say num_ghost_layers = 4, this could imply your mesh size is too large at the boundary of the geometry. See Figure 5 of https://arxiv.org/pdf/2006.05373 -> here, you would rather have a smaller mesh size on the narrow tongue than a larger ghost layer.

Thanks, again!

[amartinhuertas]

If this GridapEmbedded.jl/src/Distributed/DistributedAgFEM.jl
provides both the conforming and non-conforming cells around, yes.

It does not. In the case of non-conforming meshes, we would need to augment face_to_cells with the information stored within NonConformingGlue, in particular, in the following field: https://github.com/gridap/GridapP4est.jl/blob/d61269a6a9c2eaec0a4ef63a7dbc1b5bbc265430/src/OctreeDistributedDiscreteModels.jl#L10