18. 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 \(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.
18.1. 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.
18.1.1. 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.
18.1.2. 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.
18.1.3. 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.
// 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 );
}
18.1.4. 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.
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);
}
}
18.1.5. 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.
// 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());
18.1.6. 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.
// 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.,
\(\bullet\) for nodes:
// 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
\(\bullet\) for edges:
// 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
\(\bullet\) for faces/sides:
// 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
18.1.7. 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.
\(\bullet\) nodal field registration:
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<double>(stk::topology::NODE_RANK,
"density", numStates));
// put it on this part
stk::mesh::put_field_on_mesh(*density_, *part, nullptr);
}
\(\bullet\) edge field registration:
void
LowMachEquationSystem::register_edge_fields(
stk::mesh::Part *part)
{
const int nDim = metaData_.spatial_dimension();
edgeAreaVec_
= &(metaData_.declare_field<double>(stk::topology::EDGE_RANK,
"edge_area_vector"));
stk::mesh::put_field_on_mesh(*edgeAreaVec_, *part, nDim, nullptr);
stk::io::set_field_output_type(*edgeAreaVec_,
stk::io::FieldOutputType::VECTOR_3D);
}
\(\bullet\) side/face field registration:
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<stk::topology::rank_t>(metaData_.side_rank());
GenericFieldType *wallFrictionVelocityBip
= &(metaData_.declare_field<double>
(sideRank, "wall_friction_velocity_bip"));
stk::mesh::put_field_on_mesh(*wallFrictionVelocityBip, *part, numIp,
nullptr);
}
}
18.1.8. 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,
VectorFieldType *velocityRTM
= metaData_.get_field<double>(stk::topology::NODE_RANK, "velocity");
ScalarFieldType *density
= metaData_.get_field<double>(stk::topology::NODE_RANK, "density");}
VectorFieldType *edgeAreaVec
= metaData_.get_field<double>(stk::topology::EDGE_RANK,
"edge_area_vector");
GenericFieldType *massFlowRate
= metaData_.get_field<double>(stk::topology::ELEMENT_RANK,
"mass_flow_rate_scs");
GenericFieldType *wallFrictionVelocityBip_
= metaData_.get_field<double>(metaData_.side_rank(),
"wall_friction_velocity_bip");
18.1.9. 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,
ScalarFieldType *density
= metaData_.get_field<double>(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.
\(\bullet\) nodal field data access:
// gather some data (density at state Np1) into a local workset pointer
p_density[ni] = *stk::mesh::field_data(densityNp1, node );
- \(\bullet\) edge field data access:
(from an edge bucket loop with the same selector as defined above)
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];
}
}
18.2. High Level Nalu-Wind Abstractions
18.2.1. 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.
18.2.2. 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.