DLR-AMR / t8code

Parallel algorithms and data structures for tree-based adaptive mesh refinement (AMR) with arbitrary element shapes.
https://dlr-amr.github.io/t8code/
GNU General Public License v2.0
132 stars 52 forks source link

t8code equivalent function to p4est_nodes_new_local: rank-local node IDs #1169

Open ackirby88 opened 1 month ago

ackirby88 commented 1 month ago

Feature request

We'd like to create the t8code equivalent of p4est_nodes_new_local.

We were wondering if there are any immediate internal plans for building such a function?
Our team may attempt to build the function ourselves but before doing so, we wanted to check in to see if the team foresees any significant challenges in building the local node list following the p4est recipe.

Estimated priority Priority: high

holke commented 1 month ago

Hi Andrew! Currently, @cburstedde started implementing the full lnodes structure. The nodes_new_local is a subset of its functionality. As i understand the lnodes will take some time to develop and having the nodes_new_local quickly would be beneficial.

However, what you will most likely need is vertex neighbors and ghosts. Currently, we only have this for Quads/Hexes and single tree on an experimental branch by @lukasdreyer. @sandro-elsweijer will integrate these in the next weeks into main.

Getting them for multiple trees would require large changes in the cmesh. So for now it would be restricted to a single tree.

holke commented 1 month ago

Ping @Davknapp for organization.

holke commented 1 month ago

@ackirby88 it would help us to better understand your needs if you could specify the exact functionality that you need. (the p4est doc says 'Create node information' which is not very explicative)

ackirby88 commented 1 month ago

Thanks @holke for the response.
As you alluded to in your initial address, we are specifically looking for the nodes->num_owned_indeps count and nodes->local_nodes connectivity array across all trees and vertices on an MPI rank, (i.e., a unique vertex ID count and list).

It would be great to have ghost vertex IDs so we don't have to match them ourselves but understand it may not be easy or quick to develop that functionality. We should be able to handle those partition boundary cases in the interim.

For example, here is a use-case for our code using p4est on quads/hexes (for filling in wke->nnodes and wke->ndc8):

void setup_wake3d_allnodes(p4est_t *p4est,ctx_t *ctx){
    wake3d_t    *wke  = &ctx->wake3d;
    geometry_t  *geom = &ctx->geometry;

    p4est_topidx_t first_local_tree = p4est->first_local_tree;
    p4est_topidx_t last_local_tree = p4est->last_local_tree;
    p4est_tree_t *tree;
    sc_array_t *quadrants;
    p4est_quadrant_t *q;
    p4est_nodes_t *nodes;
    sc_array_t *indeps;
    dg::memory<char> ln;

    int lqid,jt,zz,corner;
    int numz_quadrants;
    int indep_node_id;
    int i,j;
    const int display = 0;

    /* get quadrant nodes for wake3d data */
    /*
     * NOTE:
     *  * nodes->local_nodes is a list of node indices (orders as local quadrants and corners)
     *    that locate than local node into indeps->array
     *  * nodes->local_nodes contains only grid support nodes (not high-order nodes)
     *  * indeps->array contains a list of UNIQUE local nodes
     */

    nodes = p4est_nodes_new(p4est,NULL);
    indeps = &nodes->indep_nodes;

    wke->nnodes = nodes->num_owned_indeps;
    wke->iblank_cell.malloc(p4est->local_num_quadrants,field_mask);
    wke->iblank.malloc(nodes->num_owned_indeps,field_mask);
    wke->ndc8.malloc(P4EST_CHILDREN*p4est->local_num_quadrants);
    wke->xgeom.malloc(3*wke->nnodes);

    /* local independent nodes tag */
    ln.malloc(wke->nnodes,0);

    /* note: order of nodes modified for Tioga */
    for (i = 0; i < nodes->num_local_quadrants; ++i) {
        const int ind = i*P4EST_CHILDREN;

        wke->ndc8[ind+0] = nodes->local_nodes[ind+0] + BASE;
        wke->ndc8[ind+1] = nodes->local_nodes[ind+1] + BASE;
        wke->ndc8[ind+2] = nodes->local_nodes[ind+3] + BASE; /* reversed order */
        wke->ndc8[ind+3] = nodes->local_nodes[ind+2] + BASE; /* reversed order */
#ifdef P4_TO_P8
        wke->ndc8[ind+4] = nodes->local_nodes[ind+4] + BASE;
        wke->ndc8[ind+5] = nodes->local_nodes[ind+5] + BASE;
        wke->ndc8[ind+6] = nodes->local_nodes[ind+7] + BASE; /* reversed order */
        wke->ndc8[ind+7] = nodes->local_nodes[ind+6] + BASE; /* reversed order */
#endif
    }

    /* unique grid node coordinates */
    for (jt = first_local_tree, lqid = 0; jt <= last_local_tree; ++jt) {
        tree = p4est_tree_array_index (p4est->trees, jt);
        quadrants = &tree->quadrants;
        numz_quadrants = quadrants->elem_count;

        /* loop through quadrants */
        for (zz = 0; zz < numz_quadrants; ++zz, ++lqid) {
            q = p4est_quadrant_array_index (quadrants, zz);
            quad_data_t *data = (quad_data_t *) q->p.user_data;

            /* loop through corners (real element vertices) */
            for (corner = 0; corner < P4EST_CHILDREN; ++corner) {
                /* get independent node */
                indep_node_id = nodes->local_nodes[lqid*P4EST_CHILDREN+corner];

                /* if already processed continue to next node */
                if(ln[indep_node_id]) continue;

                /* process corner geometry */
                fill_geom_corner(data,&wke->xgeom[3*indep_node_id],corner);

                /* mark the node as processed */
                ln[indep_node_id] = 1;
            }
        }
    }

    /* Print out node numbering */
    if (display) {
        for (i = 0; i < nodes->num_local_quadrants; ++i) {
            printf("------------------------------------\n");

            const int ind = i*P4EST_CHILDREN;
            for (j = 0; j < P4EST_CHILDREN; ++j) {
                Real x = wke->xgeom[3*nodes->local_nodes[ind+j]+0];
                Real y = wke->xgeom[3*nodes->local_nodes[ind+j]+1];
                Real z = wke->xgeom[3*nodes->local_nodes[ind+j]+2];
                printf("NODE [rank%d]: "
                       "Local quadrant %d; "
                       "Corner %d; "
                       "Node unique ID %d, "
                       "Node coords(x,y,z) [",
                       ctx->rank,i,j,
                       nodes->local_nodes[ind+j]);

                (x<0.0) ? printf("%f ",x):printf(" %f ",x);
                (y<0.0) ? printf("%f ",y):printf(" %f ",y);
                (z<0.0) ? printf("%f ",z):printf(" %f ",z);
                printf("]\n");
            }
        }
    }

  // more calculations ...
}
cburstedde commented 1 month ago

Thanks for bringing this up @ackirby88, and thanks @holke for mentioning our thoughts on nodes development. There is indeed a data structure that is defining a set of global nodes, its local views, and all sharer information between processes, in p4est_lnodes. There are routines for gather/allreduce type operations on nodes, where the local-to-global parallel sum of a typical finite element code is the main use case.

The lnodes are currently generated for tensor-product Lobatto points of arbitrary order and Taylor-Hood-type velocity nodes, but it can be repurposed for practically any design (if we relax the degree and vnodes members in their meaning). It is possible to fill it locally, defining all nodes as owned with no sharers, within the specification, optionally using p4est_iterate to correctly fill the face_code flag array that encodes the small-to-large topology of hanging mesh elements.

One new application of lnodes is to represent the corners of a conforming simplex mesh generated by p4est, defined by additional meta data in p4est_tnodes, which is on the alpha branch feature-triangle; see cburstedde/p4est#313.

I would hope that the lnodes and tnodes data structures are sufficiently general for your needs, and t8code (at least) conforming exports in general. This is a guess however, and it would be exciting to discuss further to meet everybody's needs. To represent a t8code mesh linked by corners, I would think that writing an algorithm that hashes all local element corners by their tree number and tree-reference (integer) coordinates and uses those to assign the local nodes would be a good first step. After that, nodes would need to be unified at tree boundaries. Unifying node numbers in parallel is an independent last step, which is described in the 2015 p4est paper by Isaac et. al. There are two different implementations to possibly adopt code from, namely @tisaac's original work in p4est_lnodes and the (currently unused, but working) code to populate nodes and sharer information in p4est_tnodes.

Curious to discuss further intergration and synergies!

cburstedde commented 1 month ago

There is some hashing code to look at in p4est_geometry, which has been added to compute the coordinates to an lnodes structure after the fact.

holke commented 1 month ago

We have a global cmesh vertex enumeration ready now (at https://github.com/DLR-AMR/t8code/tree/feature-cmesh_vertex_connectivity) and it should suffice to generate unique element node numbers across tree boundaries, since all corner-only neighbors will coincide with a tree corner...?

ackirby88 commented 1 month ago

@cburstedde We had the same thoughts regarding your point of building hashes out of the local tree ID and integer coordinates, then unifying across trees, then across boundaries. At the moment, we used only the physical spacial coordinates to build an integer hash (scaled by 1000 to reduce hash collisions), which has caused some robustness issues.

My plan is to extend the algorithm to utilize the local tree information. I know we'll need to be careful with the tree boundary unification process due to neighboring trees sharing refined face nodes. @holke Do to this fact, do you agree that the global cmesh vertex enumeration is not sufficient if it only provides info at tree corners?

@cburstedde Do you have example code in p4est for doing the rank-local tree boundary vertex unification?

holke commented 1 month ago

Hi Andrew, @lukasdreyer and I are on holidays for now and we have to come back to you on that later. Yes, I do agree that in 3D the cmesh vertex numbers might not enough due to the element edge neighbors.

cburstedde commented 1 month ago

@cburstedde We had the same thoughts regarding your point of building hashes out of the local tree ID and integer coordinates, then unifying across trees, then across boundaries. At the moment, we used only the physical spacial coordinates to build an integer hash (scaled by 1000 to reduce hash collisions), which has caused some robustness issues.

I see, and yes, pre-map hashing with tree number and tree reference coordinates may be a suitable way to ensure exact node identification.

@cburstedde Do you have example code in p4est for doing the rank-local tree boundary vertex unification?

Please check out p4est_geometry_coordinates_lnodes. If you're rolling your own algorithm, putting the results into a well defined p4est_lnodes struct would enable some synergies. I'd also be curious on your thoughts about the p4est_tnodes design.