Code Abstractions ----------------- The Nalu-Wind code base is a c++ code-base that significantly leverages the Sierra Toolkit and Trilinos infrastructure. This section is designed to provide a high level overview of the underlying abstractions that the code base exercises. For more detailed code information, the developer is referred to the Trilinos project (github.com). In the sections that follow, only a high level overview is provided. The Nalu-Wind code base emerged as a small testbed unit test to evaluate the STK infrastructure. Interestingly, the first “algorithm” implementation was a simple :math:`L_2` projected nodal gradient. This effort involved reading in a mesh, registering a nodal (vector) field, iterating elements and exposed surfaces to assemble the projected nodal gradient to the nodes of the mesh (in parallel). When evaluating kokkos, this algorithm was also used to learn about the parallel NGP abstraction provided. Sierra Toolkit Abstractions +++++++++++++++++++++++++++ Consider a typical mesh that consists of nodes, sides of elements and elements. Such a mesh, when using the Exodus standard, will liekly be represented by a collection of “element blocks”, “sidesets” and, possibly, “nodesets”. The definition of the mesh (generated by the user through commercial meshing packages such as pointwise or ICM-CFD) will provide the required spatial definitions of the volume physics and the required boundary conditions. An element block is a homegeneous collection of elements of the same underlying topology, e.g., HEXAHEDRAL-8. A sideset is a set of exposed element faces on which a boundary condition is to be applied. Finally, a nodeset is a collection of nodes. In general, nodesets are possibly output entities as the code does not exercise enforcing physics or boundary conditions on nodesets. Although Nalu-Wind supports an edge-based scheme, an edge, which is an entity connecting two nodes, is not part of the Exodus standard and must be generated within the STK infrastructure. Therefore, a particular discretization choice may require ``stk::mesh::Entity`` types of element, face (or side), edge and node. Once the mesh is read in, a variety of routine operations are generally required. For example, a low-Mach physics equation set may want to be applied to ``block_1`` while inflow, open, symmetry, periodic and wall boundary conditions can be applied to a variety of sidesets. For example, ``surface_1`` might be of an “inflow” type. Therefore, the high level set of requirements on a mesh infrastructure might be to allow one to iterate parts of the mesh and, in the end, assemble a quantity to a nodal or elemental field. Meta and Bulk Data ~~~~~~~~~~~~~~~~~~ Meta and Bulk data are simply STK containers. MetaData is used to extract parts, extract ownership status, determine the side rank, field declaration, etc. BulkData is used to extract buckets, extract upward and downward connectivities and determine node count for a given entity. Parallel Rules ~~~~~~~~~~~~~~ In STK, elements are locally owned by a single rank. Elements may be ghosted to other parallel ranks through STK custom ghosting. Exposed faces are locally owned by the lower parallel rank. Nodes are also locally owned by the lower parallel rank and can also be shared by all parallel ranks touching them. Edges and internal faces (element:face:element connectivity) have the same rule of locally owned/shared and can also be ghosted. Again, edges and internal faces must be created by existing STK methods should the physics algorithm require them. In Nalu-Wind, the choice of element-based or edge-based is determined within the input file. Connectivity ~~~~~~~~~~~~ In an unstructured mesh, connectivity must be built from the mesh and can not be assumed to follow an assumed “i-j-k” data layout, i.e., structured. In general, one speaks of downward and upward relationships between the underlying entities. For example, if one has a particular element, one might like to extract all of the nodes connected to the element. Likewise, this represents a common opporation for faces and edges. Such examples are those in which downward relationships are required. However, one might also have a node and want to extract all of the connected elements to this node (consider some sort of patch recovery algorithm). STK provides the ability to extract such connectivities. In general, full downward and upward connectivities are created. For example, consider an example in which one has a pointer to an element and wants to extract the nodes of this element. At this point, the developer has not been exposed to abstractions such as buckets, selectors, etc. As such, this is a very high level overview with more details to come in subsequent sections. Therefore, the scope below is to assume that from an element-k of a “bucket”, b[k] (which is a collection of homogeneous RANK-ed entities) we will extract the nodes of this element using the STK bulk data. .. code-block:: c++ // extract element from this bucket stk::mesh::Entity elem = b[k]; // extract node relationship from bulk data stk::mesh::Entity const * node_rels = bulkData_.begin_nodes(elem); int num_nodes = bulkData_.num_nodes(elem); // iterate nodes for ( int ni = 0; ni < num_nodes; ++ni ) { stk::mesh::Entity node = node_rels[ni]; // set connected nodes connected_nodes[ni] = node; // gather some data, e.g., density at state Np1, // into a local workset pointer to a std::vector p_density[ni] = *stk::mesh::field_data(densityNp1, node ); } Parts ~~~~~ As noted before, a ``stk::mesh::Part`` is simply an abstraction that describes a set of mesh entities. If one has the name of the part from the mesh data base, one may extract the part. Once the part is in hand, one may iterate the underlying set of entities, walk relations, assemble data, etc. The following example simply extracts a part for each vector of names that lives in the vector ``targetNames`` and provides this part to all of the underlying equations that have been created for purposes of nodal field registration. Parts of the mesh that are not included within the ``targetNames`` vector would not be included in the field registration and, as such, if this missing part was used to extract the data, an error would occur. .. code-block:: c++ for ( size_t itarget = 0; itarget < targetNames.size(); ++itarget ) { stk::mesh::Part *targetPart = metaData_.get_part(targetNames[itarget]); // check for a good part if ( NULL == targetPart ) { throw std::runtime_error("Trouble with part " + targetNames[itarget]); } else { EquationSystemVector::iterator ii; for( ii=equationSystemVector_.begin(); ii!=equationSystemVector_.end(); ++ii ) (*ii)->register_nodal_fields(targetPart); } } Selectors ~~~~~~~~~ In order to arrive at the precise parts of the mesh and entities on which one desires to operate, one needs to “select” what is useful. The STK selector infrastructure provides this. In the following example, it is desired to obtain a selector that contains all of the parts of interest to a physics algorithm that are locally owned and active. .. code-block:: c++ // define the selector; locally owned, the parts I have served up and active stk::mesh::Selector s_locally_owned_union = metaData_.locally_owned_part() & stk::mesh::selectUnion(partVec_) & !(realm_.get_inactive_selector()); Buckets ~~~~~~~ Once a selector is defined (as above) an abstraction to provide access to the type of data can be defined. In STK, the mechanism to iterate entities on the mesh is through the ``stk::mesh::bucket`` interface. A bucket is a homogeneous collection of ``stk::mesh::Entity``. In the below example, the selector is used to define the bucket of entities that are provided to the developer. .. code-block:: c++ // given the defined selector, extract the buckets of type ``element'' stk::mesh::BucketVector const& elem_buckets = bulkData_.get_buckets( stk::topology::ELEMENT_RANK, s_locally_owned_union ); // loop over the vector of buckets for ( stk::mesh::BucketVector::const_iterator ib = elem_buckets.begin(); ib != elem_buckets.end() ; ++ib ) { stk::mesh::Bucket & b = **ib ; const stk::mesh::Bucket::size_type length = b.size(); // extract master element (homogeneous over buckets) MasterElement *meSCS = sierra::nalu::get_surface_master_element(b.topology()); for ( stk::mesh::Bucket::size_type k = 0 ; k < length ; ++k ) { // extract element from this bucket stk::mesh::Entity elem = b[k]; // etc... } } The look-and-feel for nodes, edges, face/sides is the same, e.g., :math:`\bullet` for nodes: .. code-block:: c++ // given the defined selector, extract the buckets of type ``node'' stk::mesh::BucketVector const& node_buckets = bulkData_.get_buckets( stk::topology::NODE_RANK, s_locally_owned_union ); // loop over the vector of buckets :math:`\bullet` for edges: .. code-block:: c++ // given the defined selector, extract the buckets of type ``edge'' stk::mesh::BucketVector const& edge_buckets = bulkData_.get_buckets( stk::topology::EDGE_RANK, s_locally_owned_union ); // loop over the vector of buckets :math:`\bullet` for faces/sides: .. code-block:: c++ // given the defined selector, extract the buckets of type ``face/side'' stk::mesh::BucketVector const& face_buckets = bulkData_.get_buckets( metaData_.side_rank(), s_locally_owned_union ); // loop over the vector of buckets Field Data Registration ~~~~~~~~~~~~~~~~~~~~~~~ Given a part, we would like to declare the field and put the field on the part of interest. The developer can register fields of type elemental, nodal, face and edge of desired size. :math:`\bullet` nodal field registration: .. code-block:: c++ void LowMachEquationSystem::register_nodal_fields( stk::mesh::Part *part) { // how many states? BDF2 requires Np1, N and Nm1 const int numStates = realm_.number_of_states(); // declare it density_ = &(metaData_.declare_field(stk::topology::NODE_RANK, "density", numStates)); // put it on this part stk::mesh::put_field(*density_, *part); } :math:`\bullet` edge field registration: .. code-block:: c++ void LowMachEquationSystem::register_edge_fields( stk::mesh::Part *part) { const int nDim = metaData_.spatial_dimension(); edgeAreaVec_ = &(metaData_.declare_field(stk::topology::EDGE_RANK, "edge_area_vector")); stk::mesh::put_field(*edgeAreaVec_, *part, nDim); } :math:`\bullet` side/face field registration: .. code-block:: c++ void MomentumEquationSystem::register_wall_bc( stk::mesh::Part *part, const stk::topology &theTopo, const WallBoundaryConditionData &wallBCData) { // Dirichlet or wall function bc if ( wallFunctionApproach ) { stk::topology::rank_t sideRank = static_cast(metaData_.side_rank()); GenericFieldType *wallFrictionVelocityBip = &(metaData_.declare_field (sideRank, "wall_friction_velocity_bip")); stk::mesh::put_field(*wallFrictionVelocityBip, *part, numIp); } } Field Data Access ~~~~~~~~~~~~~~~~~ Once we have the field registered and put on a part of the mesh, we can extract the field data anytime that we have the entity in hand. In the example below, we extract nodal field data and load a workset field. To obtain a pointer for a field that was put on a node, edge face/side or element field, the string name used for declaration is used in addition to the field template type, .. code-block:: c++ VectorFieldType *velocityRTM = metaData_.get_field(stk::topology::NODE_RANK, "velocity"); ScalarFieldType *density = metaData_.get_field(stk::topology::NODE_RANK, "density");} VectorFieldType *edgeAreaVec = metaData_.get_field(stk::topology::EDGE_RANK, "edge_area_vector"); GenericFieldType *massFlowRate = metaData_.get_field(stk::topology::ELEMENT_RANK, "mass_flow_rate_scs"); GenericFieldType *wallFrictionVelocityBip_ = metaData_.get_field(metaData_.side_rank(), "wall_friction_velocity_bip"); State ~~~~~ For fields that require state, the field should have been declared with the proper number of states (see field declaration section). Once the field pointer is in hand, the specific field with state is easily extracted, .. code-block:: c++ ScalarFieldType *density = metaData_.get_field(stk::topology::NODE_RANK, "density"); densityNm1_ = &(density->field_of_state(stk::mesh::StateNM1)); densityN_ = &(density->field_of_state(stk::mesh::StateN)); densityNp1_ = &(density->field_of_state(stk::mesh::StateNP1)); With the field pointer already in hand, obtaining the particular data is field the field data method. :math:`\bullet` nodal field data access: .. code-block:: c++ // gather some data (density at state Np1) into a local workset pointer p_density[ni] = *stk::mesh::field_data(densityNp1, node ); :math:`\bullet` edge field data access: (from an edge bucket loop with the same selector as defined above) .. code-block:: c++ stk::mesh::BucketVector const& edge_buckets = bulkData_.get_buckets( stk::topology::EDGE_RANK, s_locally_owned_union ); for ( stk::mesh::BucketVector::const_iterator ib = edge_buckets.begin(); ib != edge_buckets.end() ; ++ib ) { stk::mesh::Bucket & b = **ib ; const stk::mesh::Bucket::size_type length = b.size(); // pointer to edge area vector and mdot (all of the buckets) const double * av = stk::mesh::field_data(*edgeAreaVec_, b); const double * mdot = stk::mesh::field_data(*massFlowRate_, b); for ( stk::mesh::Bucket::size_type k = 0 ; k < length ; ++k ) { // copy edge area vector to a pointer for ( int j = 0; j < nDim; ++j ) p_areaVec[j] = av[k*nDim+j]; // save off mass flow rate for this edge const double tmdot = mdot[k]; } } High Level Nalu-Wind Abstractions +++++++++++++++++++++++++++++++++ Realm ~~~~~ A realm holds a particular physics set, e.g., low-Mach fluids. Realms are coupled loosely through a transfer operation. For example, one might have a turbulent fluids realm, a thermal heat conduction realm and a PMR realm. The realm also holds a BulkData and MetaData since a realm requires fields and parts to solve the desired physics set. EquationSystem ~~~~~~~~~~~~~~ An equation system holds the set of PDEs of interest. As Nalu-Wind uses a pressure projection scheme with split PDE systems, the pre-defined systems are, LowMach, MixtureFraction, Enthalpy, TurbKineticEnergy, etc. New monolithic equation system can be easily created and plugged into the set of all equation systems. In general, the creation of each equation system is of arbitrary order, however, in some cases fields required for MixtureFraction, e.g., ``mass_flow_rate`` might have only been registered on the low-Mach equation system. As such, if MixtureFraction is created before LowMachEOS, an error might be noted. This can be easily resolved by cleaning the code base such that each equation system is “autonomous”. Each equation system has a set of virtual methods expected to be implemented. These include, however, are not limited to registration of nodal fields, edge fields, boundary conditions of fixed type, e.g., wall, inflow, symmetry, etc.