Merge pull request #48 from sbrus89/sbrus89/omega/mesh-class

Add Mesh class

components/omega/doc/design/HorzMeshClass.md

@@ -1,5 +1,5 @@
 (omega-design-horz-mesh)=
-# Horizontal Mesh 
+# Horizontal Mesh
 
 ## 1 Overview
 
@@ -13,7 +13,7 @@ The OMEGA mesh information should be compatible with the [MPAS Mesh Specificatio
 
 ### 2.2 Requirement: Functionality is needed to read the mesh on the host and transfer relevant data to the device for computation
 
-Not all mesh information is required in computing the tendency terms on the device, e.g. lonCell, latCell, etc. 
+Not all mesh information is required in computing the tendency terms on the device, e.g. lonCell, latCell, etc.
 However, other arrays will need to be allocated and copied to the device for use in tendency computation.
 The mesh class will explicitly include host and device YAKL arrays for each variable.
 A class method will be included to copy the relevant mesh information to the device.
@@ -26,13 +26,13 @@ Although the existing MPAS Mesh spec uses a one-based mesh numbering, zero-based
 
 The mesh class will reference the partitioned connectivity arrays created by the Decomp class.
 
-### 2.5 Requirement: Mesh variables will be associated with metadata to describe data 
+### 2.5 Requirement: Mesh variables will be associated with metadata to describe data
 
 Following the Metadata and IO designs, the YAKL arrays for the mesh variables will be associated with information about the represented values.
 
 ### 2.6 Requirement: I/O to obtain mesh data
 
-The Mesh class will have a method to read in the mesh information not obtained by the Decomp class. 
+The Mesh class will have a method to read in the mesh information not obtained by the Decomp class.
 
 ### 2.7 Desired: Ability to support multiple independent mesh objects
 
@@ -42,13 +42,13 @@ Additionally, this flexibility can be used to support separate domain decomposit
 ### 2.8 Desired: OMEGA can read in a reduced number of mesh variables and compute the remaining array information online.
 
 Many of the mesh variables are not independent, e.g.  areaCell, weightsOnEdge, etc., and can be computed from a reduced set of mesh variables.
-This functionality could be used to reduce mesh/restart file size for high resolution meshes.  
+This functionality could be used to reduce mesh/restart file size for high resolution meshes.
 Building in this flexibility would allow for all mesh-related calculations to be available within the code base instead of spreading them amongst various utility programs in MPAS-Tools.
 This will improve the ability to maintain and unit test the mesh calculations in MPAS-Tools.
 The functions used to compute the dependent mesh information can also be used to verify any mesh information that is provided in the mesh input stream.
 
 As standard practice, all necessary internally computed mesh information will be output in a single file for post-processing purposes.
-A checksumming strategy will be implemented to avoid situations where simulation data and mesh information are mismatched during post processing. 
+A checksumming strategy will be implemented to avoid situations where simulation data and mesh information are mismatched during post processing.
 
 Where appropriate, some additional derived quantities (e.g. reciprocals) will also be included to improve the performance of device-side calculations.
 
@@ -78,7 +78,7 @@ public:
   Array1DR8 AreaCell;
   ArrayHost1DR8 AreaCellH;
 
-  Array2DI4 CellsOnCell; 
+  Array2DI4 CellsOnCell;
   ArrayHost2DI4 CellsOnCellH;
 
 }
@@ -106,7 +106,7 @@ A destructor will be available to release memory.
 #### 4.2.2 Read
 The mesh class requires a method to read in all other available mesh information that has been provided in the mesh file, but has not been initialized by the decomposition. This will be a private method called by the constructor.
 
-#### 4.2.3 Compute 
+#### 4.2.3 Compute
 The compute method will be a private method called by the constructor. It will be resonsible for calculating any dependent mesh information that is not provided in the mesh input file.
 
 #### 4.2.4 Device copy creation
@@ -123,6 +123,6 @@ The metadata associated with each mesh variable will be registred within the I/O
 
 ## 5 Verification and Testing
 
-### 5.1 Test mesh compute routines 
+### 5.1 Test mesh compute routines
 
 The sample domain used for the Decomp test will be used to test obtaining the correct local values and the mesh computation routines.

components/omega/doc/devGuide/HorzMesh.md

@@ -0,0 +1,55 @@
+(omega-dev-horz-mesh)=
+
+## Horizontal Mesh
+
+The OMEGA horizontal mesh uses the [MPAS Mesh
+Specification](https://mpas-dev.github.io/files/documents/MPAS-MeshSpec.pdf).
+A mesh object is created by the `init` method, which assumes that `Decomp` has
+already been initialized.
+```c++
+Err = OMEGA::Decomp::init();
+Err = OMEGA::HorzMesh::init();
+```
+The constructor replicates the subdomain mesh cell/edge/vertex counts and
+connectivity information from Decomp so this information can be passed among the
+computational routines, alongside the other local mesh information.  It then
+creates several parallel I/O decompositions and reads in the remaining subdomain
+mesh information.  Finally, any mesh information needed on the device is copied
+from the host to a device YAKL array. Arrays such as the coordinate variables,
+which are not involved in tendency calculations, are not transferred to the
+device. These tasks are organized into several private methods. Eventually,
+dependent mesh variables will be computed from the minimum set of required mesh
+information.
+
+After initialization, the default mesh object can be retrieved via:
+```
+OMEGA::HorzMesh *HMesh = OMEGA::HorzMesh::getDefault();
+```
+Once this retrieval has been performed, public member variables can be accessed
+using:
+```
+OMEGA::I4 NCellsOwned = HMesh->NCellsOwned;
+OMEGA::Array1DR8 AreaCell = HMesh->AreaCell;
+```
+
+The HorzMesh is meant to be a container that allows the mesh information to be
+passed to the PDE solver routines:
+```
+void computeFluxTendency(OMEGA::HorzMesh *HMesh, ...) {
+yakl::c::parallel_for(yakl::c::Bounds<2>(HMesh->NCellsOwned,HMesh->MaxEdges),
+                      YAKL_LAMBDA (int Cell, int Edge) {
+  if (Edge < HMesh->NEdgesOnCell(Cell)) {
+      Var(Cell) = Var(Cell) + Flux(Cell,Edge);
+  }
+}
+```
+
+For member variables that are host arrays, variable names are appended with an
+`H`.  Array variable names not ending in `H` are device arrays.  The copy from
+host to device array is performed in the constructor via:
+```c++
+AreaCell = AreaCellH.createDeviceCopy();
+```
+
+The device arrays are deallocated by the `HorzMesh::clear()` method, which is
+necessary before calling `yakl::finalize`.

components/omega/doc/index.md

@@ -30,6 +30,7 @@ userGuide/Logging
 userGuide/Decomp
 userGuide/IO
 userGuide/Halo
+userGuide/HorzMesh
 ```
 
 ```{toctree}
@@ -50,6 +51,7 @@ devGuide/Logging
 devGuide/Decomp
 devGuide/IO
 devGuide/Halo
+devGuide/HorzMesh
 ```
 
 ```{toctree}

components/omega/doc/userGuide/HorzMesh.md

@@ -0,0 +1,43 @@
+(omega-user-horz-mesh)=
+
+## Horizontal Mesh
+
+OMEGA uses the MPAS mesh specification found
+[here](https://mpas-dev.github.io/files/documents/MPAS-MeshSpec.pdf). The names
+of the mesh variables have been retained, with the caveat that they now begin
+with a capital letter.
+
+The Mesh class is meant to be a container for all mesh variables local to a
+decomposed sub-domain that can be easily passed among the dycore routines. It
+depends on a given [Decomp](#omega-user-decomp) and reproduces the
+cell/edge/vertex totals and connectivity information from that class. The Mesh
+class also creates parallel I/O decompositions that are used to read in the
+additional mesh variables, which are not required for Decomp.
+
+Currently, the Mesh class reads in all variables from the MPAS mesh
+specification except those read by [Decomp](#omega-dev-decomp).
+This includes the following variables:
+
+| Variable Name | Description | Units |
+| ------------- | ----------- | ----- |
+| XCell, YCell, ZCell | Cartesian coordinates of cell centers | m |
+| XEdge, YEdge, ZEdge | Cartesian coordinates of edge centers | m |
+| XVertex, YVertex, ZVertex | Cartesian coordinates of vertices | m |
+| BottomDepth | Depth of the bottom of the ocean at cell centers | m |
+| FCell, FEdge, FVertex | Coriolis parameter at cell centers/edges/vertices | radians/s |
+| LonCell, LatCell | Longitude/latitude coordinates of cell centers | radians |
+| LonEdge, LatEdge | Longitude/latitude coordinates of edge centers | radians |
+| LonVertex, LatVertex | Longitude/latitude coordinates of vertices | radians |
+| AreaCell | Area of each cell | m^2 |
+| AreaTriangle | Area of each triangle in the dual grid | m^2 |
+| KiteAreasOnVertex | Area of the portions of each dual cell that are part of each cellsOnVertex | m^2 |
+| DvEdge | Length of each edge, computed as the distance between verticesOnEdge | m |
+| DcEdge | Length of each edge, computed as the distance between CellsOnEdge | m |
+| AngleEdge | Angle the edge normal makes with local eastward direction | radians |
+| MeshDensity | Value of density function used to generate a particular mesh at cell centers | - |
+| WeightsOnEdge | Reconstruction weights associated with each of the edgesOnEdge | - |
+
+In the future, the Mesh class will optionally compute the mesh variables that
+are dependent on the Cartesian mesh coordinates internally.
+This includes the various areas, lengths, angles, and weights needed for the
+TRiSK discretization (e.g. rows 5-11 in the table above).

components/omega/src/base/Decomp.cpp

@@ -352,12 +352,12 @@ int Decomp::init() {
 // NPart partitions of the mesh.
 
 Decomp::Decomp(
-    const std::string &Name,        //< [in] Name for new decomposition
-    const MachEnv *InEnv,           //< [in] MachEnv for the new partition
-    I4 NParts,                      //< [in] num of partitions for new decomp
-    PartMethod Method,              //< [in] method for partitioning
-    I4 InHaloWidth,                 //< [in] width of halo in new decomp
-    const std::string &MeshFileName //< [in] name of file with mesh info
+    const std::string &Name,         //< [in] Name for new decomposition
+    const MachEnv *InEnv,            //< [in] MachEnv for the new partition
+    I4 NParts,                       //< [in] num of partitions for new decomp
+    PartMethod Method,               //< [in] method for partitioning
+    I4 InHaloWidth,                  //< [in] width of halo in new decomp
+    const std::string &MeshFileName_ //< [in] name of file with mesh info
 ) {
 
    int Err = 0; // internal error code
@@ -371,7 +371,8 @@ Decomp::Decomp(
 
    // Open the mesh file for reading (assume IO has already been initialized)
    int FileID;
-   Err = IO::openFile(FileID, MeshFileName, IO::ModeRead);
+   MeshFileName = MeshFileName_;
+   Err          = IO::openFile(FileID, MeshFileName, IO::ModeRead);
    if (Err != 0)
       LOG_CRITICAL("Decomp: error opening mesh file");
 

components/omega/src/base/Decomp.h

@@ -135,6 +135,8 @@ class Decomp {
    // Since these are used frequently, we make them public to reduce the
    // number of retrievals required.
 
+   std::string MeshFileName; ///< The name of the file with mesh info
+
    // Sizes and global IDs
    // Note that all sizes are actual counts (1-based) so that loop extents
    // should always use the 0:NCellsXX-1 form.
@@ -226,7 +228,7 @@ class Decomp {
           I4 NParts,               ///< [in] num of partitions for new decomp
           PartMethod Method,       ///< [in] method for partitioning
           I4 InHaloWidth,          ///< [in] width of halo in new decomp
-          const std::string &MeshFileName ///< [in] name of file with mesh info
+          const std::string &MeshFileName_ ///< [in] name of file with mesh info
    );
 
    /// Destructor - deallocates all memory and deletes a Decomp.