Adding global edge and face ID

[karthichockalingam]

Hello,

Similar to global node and element numbers/IDs. We are interested in adding a new global edge and face numbering. This is particularly useful in preconditioning Nedelec and Raviart-Thomas elements using Hypre AMS (https://hypre.readthedocs.io/en/latest/solvers-ams.html) and ADS (https://hypre.readthedocs.io/en/latest/solvers-ads.html) preconditioners respectively.

Can you please give us some pointers on how we can create such a data structure in Libmesh? Perhaps something lightweight to begin with which doesn't mess with the I/O/.

Many thanks,
Karthik.

cc @SohailSTFC

[roystgnr]

If you're using what we call "complete order" elements, then every face has its own node and every edge has its own node, and all those have their own id. The numbering is sparse (because vertices and interior nodes still have their ids, but that still qualifies as "similar to global node and element numbers/IDs", which can also be sparse (because a mesh file had them so, or because there's been element coarsening/deletion without renumbering).

If you need a contiguous numbering (in either case)? A contiguous->sparse map (really a vector) is cheap, and easy to parallelize since you can choose to keep your own numbering contiguous on each processor (as we do with dof ids, though not node/elem ids). A sparse->contiguous map would actually have to be a map (or unordered_map or a cool_third_party_thing, whatever), so that's a bit less lightweight, but you could piggyback on existing libMesh overhead with MeshBase::add_node_integer() to store the contiguous id with the Node itself.

[SohailSTFC]

Hi Roy, Hope all is well, I wanted to ask if their is a preferential method for numbering the Edge and face Ids in a global way within LibMesh, To summarise:
-Is there a way to separate entities into Volumes, Faces, Edge, Vertices with unique contiguous yet separate numberings from the Nodal numberings within LibMesh.
-Is there a preferred communication handle (MPI or other) that could be used to communicate numberings of non-local sub-entities (Volumes, Faces, Edges, Nodes)
-Is there an existing framework-style for this such that it could be conveniently called from LibMesh as oppose to an externally hacked together code
Many Thanks!
SohailSTFC

[roystgnr]

Is there a way to separate entities into Volumes, Faces, Edge, Vertices with unique contiguous yet separate numberings from the Nodal numberings within LibMesh.

Yes; C++ is Turing-complete. 😉

Is there a preferred communication handle (MPI or other) that could be used

Every mesh will give you a mesh.comm() to use; this is just an MPI_COMM_WORLD wrapper by default but if a mesh has been created with a subcommunicator it will return that instead.

Is there an existing framework-style for this such that it could be conveniently called from LibMesh as oppose to an externally hacked together code

The Elem class has methods like is_vertex() etc. that can be called on local node numbers to determine which type of node each is; other than that you're doing something new here.

[karthichockalingam]

Hi Roy,

Thanks much for your response.

Q(1) We are thinking of using map for faces and edges to create a contiguous id. Would this map go in MeshBase or some other class?

Q(2) In which function would we place this map construction? It would have to be after the mesh is read/built.

Q(3) Do the processes need to communicate at all to create a secondary continuous ID for example for faces?

We are planning to use this map to build the Hypre AMS and ADS preconditioners in an existing example.

[SohailSTFC]

Hi Roy,
Also a further question, is there a place where you store the neighbouring processes in LibMesh (i.e processes that share nodes/an interface)?
Or more generally, as I am looping over all the elements to find their respective sub-entities (Nodes (which can be Verts, Edges, Faces, Vols) ) is there a global node-number-to-PID ownership that I can use?

[roystgnr]

Would this map go in MeshBase or some other class?

I know you intended that as an "xor" rather than an "or", but you're in danger of the answer being simply "no".

Since my answer to "will you remove this library feature to make maintenance easier" is almost always no, and yet the answer to "will indefinite maintenance of this library feature be hassle-free" is always no, we've got to be at least a bit wary when adding new features, and the closer the new feature has to go to the library core, the more wary we should be. Migrating working middleware-level or application-level features deeper can be persuasive; starting at the library level while still trying to figure out the very basics (e.g. should lookup really be O(log N), how do I find a processor id) is not. Let's see if we can come up with something that will work even before having to change the library to accommodate it?

I'm not even convinced we're not making a huge XY Problem mistake here, but if we want to put overhead in the library then that needs to be fixed. DoF ids are not node ids, so why does an algebraic solver care whether a subset of node ids are sparse? We can iterate through nodes and get a sequential numbering for a subset on the fly as we do so, so what do we need the reverse map for? If we're calling a custom solver rather than using the existing solver interfaces then shouldn't any solver-specific overhead be in that custom solver?

In which function would we place this map construction?

It would obviously be in its own function; I can't pretend that libMesh is as well-factored as it should be but I can at least try to avoid making things worse. ;-)

Where to call that function is a more interesting question, though. For the most general usability it should obviously be from MeshBase::prepare_for_use(), but most of our clever callback opportunities are at the System level rather than the Mesh level. System::user_initialization() (whether via an override or attaching an initialization function or Initialization object) would be appropriate here, since that'll be called in between mesh changes and solves, but we might want to add something like that at the MeshBase level too and use that instead in the long term ... or we might want to leave everything in your solver class and just let initialization and reinitialization there do it.

Do the processes need to communicate at all to create a secondary continuous ID for example for faces?

On a serialized mesh no communication is necessary, but even there it can speed things up. On a distributed mesh there's no getting around communication for this.

Do Hypre AMS and ADS work in parallel?

is there a place where you store the neighbouring processes in LibMesh (i.e processes that share nodes/an interface)?

No; we use a non-blocking exchange ("NBX" algorithm, in timpi/parallel_sync.h) to avoid needing to cache that.

is there a global node-number-to-PID ownership that I can use?

mesh.node_ref(n).processor_id() if you already know that n is ghosted on the processor you're querying from; if you don't know whether or not it's ghosted, mesh.query_node_ptr(n) returns either a ghosted node or a nullptr.

[roystgnr]

DoF ids are not node ids, so why does an algebraic solver care whether a subset of node ids are sparse?

Wait a minute. You're not trying to create low-order Raviart-Thomas and Nedelec (one) elements from scratch, are you? If all you need is a numbering of edge degrees of freedom in 2D or face degrees of freedom in 3D, we already do that whenever the user adds a RAVIART_THOMAS or NEDELEC_ONE element with p=1.

[karthichockalingam]

You're not trying to create low-order Raviart-Thomas and Nedelec (one) elements from scratch, are you?

Not really, we want to create a G-operator which is of size number of vertices (N) times the number of edges (E), which has two entree columns (denoting the orientation of an edge). This is used for the AMS preconditioner. For ADS, in addition to the G-operator, we have a C-operator denoting the orientation of a face. Hence an id to uniquely denote an edge and a face can be very handy. This helps determine the size of the operators.

What we find is when solving the A-formulation using Nedelec elements we need to deploy direct solvers. A preconditioned solution would allow us to avoid solving it directly.

Do Hypre AMS and ADS work in parallel?

Yes, they do.

Let’s see if we can come up with something that will work even before having to change the library to accommodate it

That is exactly, what we are doing now, computing the orientation on the fly directly on a standalone example. We can send you whenever we have a working version for you to have a look.

[karthichockalingam]

Hello @roystgnr

Here, EdgeLocalID[nodeID] returns the contiguous list of edge-node-ids, which is constructed by sequentially looping over the edges from one element to the next. I am working with 1st-order RT elements, with one dof per edge. My question is will the contiguous edge ID list match the dof list (assuming serial case for now without amr)? If not, how can I make them match?

  // Form the edge-node pairing
  for (const auto & elem : mesh.active_local_element_ptr_range()){
    for(unsigned int I=0; I<elem->n_nodes(); I++){
      unsigned int nodeID = elem->node_id(I);
      const Node & node = mesh.node_ref(nodeID);
      if( node.processor_id() == mesh.processor_id()){
        if(elem->is_edge(I)) 
           if( EdgeLocalID.find(nodeID) == EdgeLocalID.end() ){
                 EdgeLocalID[nodeID] = J;
              J++;}
      }
    }
  }

[roystgnr]

Unrelated aside: prefer elem->node_ref(I) over mesh.node_ref(nodeID) in cases where you have I - the latter may be a tiny bit slower in some cases and O(log Nnodes) instead of O(1) in others.

In serial, without AMR, with a single non-subdomain-restricted 1st-order RT variable, I think your EdgeLocalID will match the DoF ordering. It is not guaranteed to match the DoF ordering, even in that case, though; this makes the difference between "if a future libMesh change breaks my code I should report it as a bug so it can be fixed" and "if a future libMesh change breaks my code I should keep quiet so I can't be laughed at".

To guarantee a match (whether because you're going to parallel or AMR or subdomain restrictions or more variables or because you wisely don't trust my "I think" to be unimpeachably correct), query the DoF ordering. node->dof_number(system_number, variable_number, component_number) will give you the DoF id you want, and if you've just got one system with just one variable that's lowest order then it's just node->dof_number(0,0,0). Make sure to only query it on nodes where there is a DoF. In your case the is_edge(I) is sufficient; in general node->n_comp(system_number, variable_number) says how many DoFs ("components" here) belong to that variable on that node, and should exceed any component_number you try to query.

[karthichockalingam]

Hi Roy,

Thank you for your response. I have constructed the G operator for Hypre AMS and the plan is to use it to solve vector_fe_ex3 (curl curl problem). I made use of node->dof_number(0,0,0) to ensure the rows align with the dofs.

// Set discrete gradient
petscErr = PCHYPRESetDiscreteGradient(pc, G);

// Set vertex coordinates
petscErr = PCHYPRESetEdgeConstantVectors(pc, par_xvec, par_yvec, par_zvec);

The preconditioner was constructed using PETSc type PC. How best to pass the above preconditioner pc in vector_fe_ex3.C to the solve?

[karthichockalingam]

Hi Roy, I am posting a related question. Look's like we can get the reference to pc from
https://github.com/libMesh/libmesh/blob/devel/include/solvers/petsc_linear_solver.h#L230

I need to call the following PETSc functionality from Libmesh example
https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/pc/impls/hypre/hypre.c?ref_type=heads#L1640

For that, I would need to know how I get the reference to the system matrix A from libmesh?

Thank you!

[roystgnr]

Yeah; in that example the DiffSolver will be a NewtonSolver by default; you ought to be able to dynamic_cast to that, then get_linear_solver(); the resulting LinearSolver will be a PetscLinearSolver ... but we don't actually have a setter for the PC there, only a getter? You could get the ksp() and set your preconditioner with that, though.

The matrix is literally System::matrix; you can dynamic_cast that to a PetscMatrix to get an accessor for the underlying Mat.

[karthichockalingam]

It does look like we can set pc using pc() from PetscLinearSolver maybe because it was returning a const pointer?

NewtonSolver & newton = dynamic_cast<NewtonSolver &>(solver);
const PC & pc = dynamic_cast<PetscLinearSolver<double> &>(newton.get_linear_solver()).pc();
const KSP & ksp = dynamic_cast<PetscLinearSolver<double> &>(newton.get_linear_solver()).ksp();

But instead, you are proposing to use KSPGetPC to set pc?

Many thanks!

[roystgnr]

It's returning a PC, not a const PC &. You're allowed to take a const reference to a temporary, because otherwise you couldn't pass it to functions that take a const reference, but you're not allowed to take a non-const reference to it, because that rule catches errors like this one where you think you're going to change an lvalue somewhere but you've really only got your hands on a temporary copy.

That might be enough here, since under the hood a PC is just a handle itself to PETSc ... but I think KSPGetPC is the way to go.

[karthichockalingam]

I have a follow up question. Is the following error because I am missing a header in vector_fe_ex3.C?

  // The global system matrix
  SparseMatrix<Number> & sparse_matrix = system.get_system_matrix();

  PetscMatrix<Number> & petsc_matrix = dynamic_cast<PetscMatrix<Number> &>(sparse_matrix);

  Mat A = petsc_matrix.mat();

  PCHYPRESetPoissonMatrix_HYPRE(pc, A, false);

vector_fe_ex3.C:141:3: error: 'PCHYPRESetPoissonMatrix_HYPRE' was not declared in this scope; did you mean 'PCHYPRESetBetaPoissonMatrix'?
  141 |   PCHYPRESetPoissonMatrix_HYPRE(pc, A, false);
      |   ^~~~~~~~~~~~~~~~~~~~~~~~~~~~

[roystgnr]

Almost certainly. In general you don't want to rely on getting a necessary header indirectly; if you need Foo then add #include <foo.h>, don't rely on #include <bar.h> picking it up for you or picking up another header that picks it up for you.