Figure out a mechanism for tagging faces in an Exodus mesh via DMPlex

[jeff-cohere]

We're interested in tagging sets of faces in our meshes so we can assign boundary conditions to them. Here's Matt's advice:

---

The way DMPlex handles this is to make a DMLabel object.

In ExodusII, you would mark the top of the domain with a "side set". Each side set has an id. When it is translated to DMPlex, all sides sets are placed in the "Face Sets" label. All the faces belonging to a side set are tagged with the side set id.

---

PETSc's Example 11 demonstrates how to use DMLabels with DMPlex to assign boundary conditions to faces. See also DMGetLabelByNum, DMGetLabel, DMGetNumLabels.

[knepley]

I think SNES ex12 is easier to understand, but lots of them do it. Is it clear how to make a sides set in ExodusII?

[jeff-cohere]

Thanks for the pointer! I've created a sides set in Exodus before, but it was a long time ago. I probably do still have the code that I used for it. If you happen to have anything handy, though, I'd welcome the info.

[knepley]

I have only ever done it in Cubit. If you are doing it programmatically, I would likely just make the DMPlex with the label I wanted and view it with the Exodus viewer that Blaise wrote. That should output the label correctly.

[jeff-cohere]

One example Exodus mesh in the PETSc repo that contains side sets is located in share/petsc/datafiles/meshes/blockcylinder-50.exo. This is a block-structured cylindrical mesh with the following side sets indicating the inner/outer (radial) and the top and bottom boundaries:

• InnerBoundary
• OuterBoundary
• TopBoundary
• BottomBoundary

This mesh is read in the tutorial example src/dm/label/tutorials/ex1.c, which writes all the labels discovered in the mesh to an output file.

Output for this tutorial example appears in src/dm/label/tutorials/output/ex1_0.out. Within this file, one can see the Face Sets label mentioned by Matt above, which includes faces that belong to all four of the cylinder's boundaries. We can also see individual labels for each of these boundaries:

petsc/src/dm/label/tutorials/output$ cat ex1_0.out
Number of labels: 12
Label 0: name: celltype
IS of values
  IS Object: 1 MPI processes
    type: general
  Number of indices in set 4
  0 0
  1 7
  2 4
  3 1
...

Label 3: name: Face Sets
IS of values
  IS Object: 1 MPI processes
    type: general
  Number of indices in set 4
  0 100
  1 101
  2 200
  3 201

Label 4: name: OuterBoundary
IS of values
  IS Object: 1 MPI processes
    type: general
  Number of indices in set 1
  0 100

Label 5: name: InnerBoundary
IS of values
  IS Object: 1 MPI processes
    type: general
  Number of indices in set 1
  0 101

Label 6: name: TopBoundary
IS of values
  IS Object: 1 MPI processes
    type: general
  Number of indices in set 1
  0 200

Label 7: name: BottomBoundary
IS of values
  IS Object: 1 MPI processes
    type: general
  Number of indices in set 1
  0 201
...

@bishtgautam

[bishtgautam]

@jeff-cohere , The .exo file is helpful and based on the ncdump, the mesh has four sidesets.

ncdump -h $PETSC_DIR/share/petsc/datafiles/meshes/blockcylinder-50.exo
netcdf blockcylinder-50 {
dimensions:
    len_string = 33 ;
    len_line = 81 ;
    four = 4 ;
    time_step = UNLIMITED ; // (0 currently)
    num_dim = 3 ;
    num_nodes = 125 ;
    num_elem = 56 ;
    num_el_blk = 1 ;
    num_qa_rec = 1 ;
    num_side_sets = 4 ;
    num_side_ss1 = 64 ;
    num_df_ss1 = 256 ;
    num_side_ss2 = 24 ;
    num_df_ss2 = 96 ;
    num_side_ss3 = 14 ;
    num_df_ss3 = 56 ;
    num_side_ss4 = 14 ;
    num_df_ss4 = 56 ;
    num_el_in_blk1 = 56 ;
    num_nod_per_el1 = 8 ;
variables:
    double time_whole(time_step) ;
    char qa_records(num_qa_rec, four, len_string) ;
    char coor_names(num_dim, len_string) ;
    char eb_names(num_el_blk, len_string) ;
    int ss_status(num_side_sets) ;
    int ss_prop1(num_side_sets) ;
        ss_prop1:name = "ID" ;
    char ss_names(num_side_sets, len_string) ;
    int elem_ss1(num_side_ss1) ;
    int side_ss1(num_side_ss1) ;
    double dist_fact_ss1(num_df_ss1) ;
    int elem_ss2(num_side_ss2) ;
    int side_ss2(num_side_ss2) ;
    double dist_fact_ss2(num_df_ss2) ;
    int elem_ss3(num_side_ss3) ;
    int side_ss3(num_side_ss3) ;
    double dist_fact_ss3(num_df_ss3) ;
    int elem_ss4(num_side_ss4) ;
    int side_ss4(num_side_ss4) ;
    double dist_fact_ss4(num_df_ss4) ;
    int elem_map(num_elem) ;
    int eb_status(num_el_blk) ;
    int eb_prop1(num_el_blk) ;
        eb_prop1:name = "ID" ;
    int connect1(num_el_in_blk1, num_nod_per_el1) ;
        connect1:elem_type = "HEX8" ;
    double coordx(num_nodes) ;
    double coordy(num_nodes) ;
    double coordz(num_nodes) ;

// global attributes:
        :api_version = 4.98f ;
        :version = 4.98f ;
        :floating_point_word_size = 8 ;
        :file_size = 1 ;
        :title = "cubit(blockcylinder-50.exo): 02/04/2013: 17:50:28" ;
}

The first sideset (*ss1) has following values:

[elem_ss1' side_ss1']
    6    3
    7    3
   10    2
   14    2
   20    3
   21    3
   24    2
   28    2
   34    3
   35    3
   38    2
   42    2
   48    3
   49    3
   52    2
   56    2
    8    2
    9    1
   10    1
   11    3
   22    2
   23    1
   24    1
   25    3
   36    2
   37    1
   38    1
   39    3
   50    2
   51    1
   52    1
   53    3
    1    3
    2    4
    8    1
   13    2
   15    3
   16    4
   22    1
   27    2
   29    3
   30    4
   36    1
   41    2
   43    3
   44    4
   50    1
   55    2
    4    2
    6    4
   12    3
   13    1
   18    2
   20    4
   26    3
   27    1
   32    2
   34    4
   40    3
   41    1
   46    2
   48    4
   54    3
   55    1

I'm assuming the elem_ss1 corresponds to the cell/control volume and side_ss1 corresponds to the face of each element. But, I'm don't know which face corresponds to the 3rd face of the 6th cell (first entries in elem_ss1' side_ss1'). The exodus manual in Figure 8/page 15 shows how vertices should be ordered to form different elements but doesn't tell us the order of faces.

In the TDycore, we need to

• read in all num_side_sets from .exo file, and
• map the side_ss* to the faces of each element.

@ghammond86 Do you know if the order of faces that we use in the unstructured grid format PFLOTRAN is the same as the Exodus II format?

cc: @knepley

[knepley]

Why write new code? Can't you just use the Plex I/O for Exodus?

[ghammond86]

@ghammond86https://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fghammond86&data=04%7C01%7Cglenn.hammond%40pnnl.gov%7Cdf2ba684b17f4dc3dc7908d90bef828e%7Cd6faa5f90ae240338c0130048a38deeb%7C0%7C1%7C637553945628722177%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=2N%2FCyxA4eyD%2FQEAODsBIg%2Boa671hbQ2ZkaW6cYe48w8%3D&reserved=0 Do you know if the order of faces that we use in the unstructured grid format PFLOTRAN is the same as the Exodus II format?

I believe so. That was the intent. Glenn

[bishtgautam]

@knepley We definitely don't want to write new code. My intent is to figure out how to correctly mark the faces in the .exo file.

[knepley]

@bishtgautam Cool. How about we try writing out the Plex+Vecs in the HDF5 file using the HDF5 viewer. Blaise has just worked on this and it should work. We could quickly find out if this will work for TDycore. Then all it would take is making the label and writing back out.

[jeff-cohere]

Hey @knepley ,

I'm happy to give this a try. Can you describe this a bit more? Do we want a separate program/utility that takes in an Exodus mesh without side sets, and adds side sets for specific faces, with specific label names? Or do you imagine having this be part of our "driver"?

[knepley]

Sure. I mistyped above. I meant Exodus viewer, not HDF5. We can also use that one. Here is how I would set things up. Make a small example with a VecView() and VecLoad(). We have a nice VecViewFromOptions(), but unfortunately I have not put in VecLoadFromOptions() yet. You make the kind of viewer you want, either Exodus or HDF5 (or ASCII...). We check that the Vec your write out is what you read in.