Question on DOF ordering

[jwpeterson]

This question was emailed to me by Liu Haochen, I've copied it here so that others can see/comment. --John

Dear developers,

We'd like to consult the ordering of freedoms in a global matrix arising from multi-variable problems such as the Navier-Stokes equations, as addressed in ../examples/systems_of_equations/ex2.

1. If we use a linear approximation for velocity and pressure (P1-P1), is the global matrix ordered in a node-by-node form as follows?

2. What is the ordering of the global matrix if we use "SECOND" order for velocity (P2-P1)? Since there are more dofs for velocity than pressure.

3. For the example "../examples/systems_of_equations/ex2", we find that it fails when we use the option ' -sub_pc_factor_mat_ordering_type rcm ', therefore we are confused with the ordering.

We are looking forward to your reply. Thank you so much!

[jwpeterson]

By default the variable ordering in libMesh is so-called "variable major", so in your example all the u components would appear first, followed by all the v components, and finally all the p components. Please note however that the variable ordering can be changed via command line argument, so one should avoid writing code which assumes any particular ordering.

we find that it fails when we use the option ' -sub_pc_factor_mat_ordering_type rcm '

There could be many reasons for the failure, but I wouldn't initially suspect it is necessarily related to libMesh's DOF ordering. I don't think I've ever tried this particular PETSc option before, and we definitely don't have any testing for it.

[roystgnr]

That's right (about the P1-P1 case; I always thought "Taylor Hood" implied a velocity polynomial degree one larger than the pressure polynomial degree).

[LHappyCureall]

I'm really sorry to bother you again. For the P1-P1, because it is ordered node by node, we can use block size 3 for the matrix storage. But for P2-P1, we don't know how to set the block size. Would you please let us know about the exact ordering when using P2-P1 if using the "--node-major-dofs" option?

[jwpeterson]

Would you please let us know about the exact ordering when using P2-P1 if using the "--node-major-dofs" option?

I think it will be like @roystgnr said:

In a Taylor Hood (u,v,p) case where u and v match but p doesn't, IIRC the ordering will look something like [u0, v0, u1, v1, u2, v2, ... p0, p1, p2,...]

But you should probably write a small test code that loops over the nodes and prints

node->dof_number(sys_num, var_num, 0);

Note that not every node will have p DOFs in this case, and I think dof_number() throws an error if you ask for a DOF number that doesn't exist.

[roystgnr]

For P2-P1 there is no valid block size, even with --node-major-dofs. Vertices will have [ui, vi, pi]; edge nodes (and if you use Q2-Q1 quads or hexes too, interiors and face nodes) will have just [ui, vi]. The ordering is consistent, but some nodes will have DoFs that others won't. You can't even group multiple nodes into a single global block size: even if you're in serial and you ensure the node numbering is locally consistent with that, the node count won't be. E.g. a regular NxN grid of 2N² triangles will have N²+2N+1 vertices but 3N²+2N edge nodes; asymptotically you have 9 DoFs per P2-P1 triangle pair but there's remainders at borders.

[roystgnr]

dof_number() throws an error if you ask for a DOF number that doesn't exist.

It does so in METHOD=dbg or devel builds. In profiling or optimized builds (or any custom build without NDEBUG defined) asking for a DoF number that doesn't exist summons nasal demons.

If you're not sure a variable has any components on node, you can use node->n_comp(sys_num, var_num) first to check.