================ Mesh Preparation ================ Cocoa uses unstructured meshes stored in NetCDF format following UGRID conventions. Mesh Requirements ----------------- Cocoa meshes must satisfy these requirements: - Element type: linear (P1) triangular elements, bilinear (Q1) quadrilateral elements, or a hybrid mesh mixing both in a single MPI partition. See :ref:`quad-elements` for the quadrilateral discretization. - Connectivity: a connected mesh with no isolated nodes - Quality: reasonable element aspect ratios and low (<8) valency; quadrilaterals must be convex and simple (non-self-intersecting) - Numbering: consistent counter-clockwise node ordering (Cocoa reorients clockwise elements automatically and logs the count it fixed at debug level) - Values: every node has finite coordinates and depth (no NaN, infinity, or fill value) and a latitude within 90 degrees of the equator, and every triangle has a non-zero area. Cocoa rejects a mesh that breaks these when it reads it, naming the file and the node or element. - Format: NetCDF with UGRID conventions Mesh Format ----------- Cocoa uses a NetCDF-4 format based on UGRID conventions. The format stores mesh geometry, boundaries, optional nodal attributes, and optional self-attraction and loading (SAL) data in a single file, organized into NetCDF groups. The file uses a grouped layout: - The root group holds only the UGRID mesh topology: the ``cocoa_mesh`` dummy variable, the node coordinates (``x``, ``y``), the bathymetry (``depth``), and the element connectivity (``mesh_face_nodes``). - A ``/boundaries`` group holds the boundary segment data (present only if the mesh has boundaries). - A ``/nodal_attributes`` group holds one variable per nodal attribute (present only if the mesh has nodal attributes). - A ``/self_attraction_loading`` group holds SAL tidal data (present only if the mesh carries SAL data). Keeping the core topology in the root group lets generic UGRID/CF tooling (QGIS, ParaView, ``xarray``, MDAL) discover and render the mesh, while the cocoa-specific boundary, attribute, and SAL blocks live in their own groups. .. note:: Cocoa does not read ADCIRC fort.14 files directly. Use the ``cocoa-mesh-tools`` command to convert ADCIRC meshes to Cocoa format. See :ref:`converting-from-adcirc` below. Dimensions ^^^^^^^^^^ .. list-table:: :header-rows: 1 :widths: 25 20 55 * - Dimension - Group - Description * - ``nnode`` - root - Number of mesh nodes * - ``nface`` - root - Number of mesh elements (faces/triangles) * - ``max_face_nodes`` - root - Maximum nodes per face: 3 for a triangle-only mesh, 4 for a mesh that encodes any quadrilateral (see :ref:`hybrid-connectivity` below) * - ``nboundary`` - /boundaries - Number of boundary segments * - ``nboundary_node`` - /boundaries - Total number of boundary nodes across all segments * - ``nsal_constituents`` - /self_attraction_loading - Number of tidal constituents in the SAL data * - ``sal_name_len`` - /self_attraction_loading - Character length of SAL constituent names Required Variables ^^^^^^^^^^^^^^^^^^ **Root group (mesh topology):** .. list-table:: :header-rows: 1 :widths: 20 15 65 * - Variable - Dimensions - Description * - ``cocoa_mesh`` - (scalar) - UGRID mesh topology variable with UGRID attributes * - ``x`` - (nnode) - Node x-coordinates (longitude for spherical) [degrees_east] * - ``y`` - (nnode) - Node y-coordinates (latitude for spherical) [degrees_north] * - ``depth`` - (nnode) - Bathymetric depth at each node [m, positive down] * - ``mesh_face_nodes`` - (nface, max_face_nodes) - Element connectivity; see :ref:`hybrid-connectivity` for the ``max_face_nodes=4`` encoding and below for ``start_index`` .. _hybrid-connectivity: Triangle/Quadrilateral Connectivity Encoding ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ``mesh_face_nodes`` honors the variable's ``start_index`` attribute, which must be exactly ``0`` (UGRID's own default, and every fixture shipped here) or ``1`` (the convention of files numbered from 1); Cocoa rejects any other value at startup. Fill entries are recognized before the shift is applied and are never shifted: an entry equal to the declared ``_FillValue`` (or any negative value, when no ``_FillValue`` is declared) is a sentinel, not a node id. A triangle-only mesh writes ``mesh_face_nodes`` with ``max_face_nodes=3`` and every row fully populated. A mesh with any quadrilateral element writes ``max_face_nodes=4``: every row is either a triangle, using its first 3 columns with the 4th column carrying an unused ("fill") node index, or a quadrilateral, using all 4 columns. - Fill (unused) slots must be trailing: a triangle row's 4th column is unused, never its 1st, 2nd, or 3rd. - The unused-slot value is the variable's declared ``_FillValue`` attribute when present, or any negative value when the attribute is absent (the plain UGRID convention). Cocoa reads whichever the file declares; it does not require any particular fill value. - Column count must be exactly 3 or 4; any other width is rejected at read time. - Element order in the file is preserved end to end. Cocoa classifies each row into a triangle or a quadrilateral list once at read time, and file order stays the external element id space: it is what the NetCDF output connectivity, the partition cache, and global element ids all use, in serial and under MPI alike, and checkpoint files carry the same file-order ``element`` table at the same 3- or 4-column width. - The mesh id, checkpoints, and partition caches are the same regardless of encoding: the id hashes a canonical topology (fill stripped, an encoding-independent 3- or 4-wide array) rather than the raw file bytes, so a triangle-only mesh gets the identical id whether it is stored with ``max_face_nodes=3`` or ``4``, and anything keyed on the id (references, partition caches, checkpoints) is shared across the two encodings. A mesh with real quadrilaterals gets a different id from its triangle-split equivalent, since they describe different topology. - A quadrilateral's 4 nodes must be distinct and in range, and the resulting quadrilateral must be convex and simple (not self-intersecting); Cocoa reorients a clockwise quadrilateral automatically, the same as it does for triangles, but does not otherwise repair a malformed one; it fails at startup naming the file element index (both 0-based and 1-based) and the offending node ids. A triangle-only mesh encoded with ``max_face_nodes=4`` (every row a triangle, column 4 unused) reads and runs identically to the same mesh encoded with ``max_face_nodes=3``; the encoding only sets the connectivity-table width. **/boundaries group:** This group is present only when the mesh has boundaries. It owns the ``nboundary`` and ``nboundary_node`` dimensions. The first four variables are always present; the remaining variables are optional and are written only when they carry real (non-fill) data. For example, a mesh with no weirs/levees/pipes omits ``boundary_pair_node`` and the barrier/pipe variables entirely rather than storing arrays full of ``-9999``. .. list-table:: :header-rows: 1 :widths: 38 22 40 * - Variable - Dimensions - Description * - ``boundary_types`` - (nboundary) - ADCIRC boundary type code (-1=open, 0/1=land, etc.) * - ``boundary_start_index`` - (nboundary) - Start index into boundary_nodes for each segment * - ``boundary_end_index`` - (nboundary) - End index into boundary_nodes for each segment * - ``boundary_nodes`` - (nboundary_node) - Node indices for all boundary nodes * - ``boundary_pair_node`` (optional) - (nboundary_node) - Paired node for internal barriers (weir/levee) * - ``boundary_elevation`` (optional) - (nboundary_node) - Barrier crest elevation [m] * - ``boundary_supercritical_coefficient`` (optional) - (nboundary_node) - Supercritical flow coefficient * - ``boundary_subcritical_coefficient`` (optional) - (nboundary_node) - Subcritical flow coefficient * - ``boundary_pipe_coefficient`` (optional) - (nboundary_node) - Cross-barrier pipe flow coefficient * - ``boundary_pipe_height`` (optional) - (nboundary_node) - Pipe centerline height [m] * - ``boundary_pipe_diameter`` (optional) - (nboundary_node) - Pipe diameter [m] **/nodal_attributes group:** This group is present only when the mesh has nodal attributes. It contains one variable per attribute, dimensioned ``(nnode)`` for scalar attributes or ``(nnode, na_length_N)`` for multi-value attributes (for example, ``surface_directional_effective_roughness_length`` is ``(nnode, ndir)`` with ``ndir=12``). The on-disk variable order is irrelevant. The available attributes are described in :ref:`nodal-attributes` below. **/self_attraction_loading group:** This group is present only when the mesh carries SAL tidal data. It owns the ``nsal_constituents`` and ``sal_name_len`` dimensions. .. list-table:: :header-rows: 1 :widths: 30 30 40 * - Variable - Dimensions - Description * - ``sal_constituent_names`` - (nsal_constituents, sal_name_len) - Tidal constituent names (char) * - ``sal_frequency`` - (nsal_constituents) - Constituent angular frequency [rad/s, double] * - ``sal_amplitude`` - (nnode, nsal_constituents) - SAL amplitude per node and constituent [m, float] * - ``sal_phase`` - (nnode, nsal_constituents) - SAL phase per node and constituent [degrees, float] .. _mesh-identification: Mesh Identification ------------------- Every Cocoa mesh has a mesh id: a short, filename-safe fingerprint computed from its contents. The same mesh always produces the same id, regardless of the machine it is computed on. Any change to the node coordinates, bathymetry, connectivity, boundaries, nodal attributes, or SAL data changes the id. Format ^^^^^^ A mesh id has the form:: cm1_P<6hex>_Z<6hex>_E<6hex>_B<6hex>_A<6hex>_S<6hex> for example:: cm1_P4b5652_Z4d037a_Ec05a73_Bef8a11_A65085e_Scec7ea The ``cm1`` prefix is the scheme/version tag. Each of the six tagged sections is a BLAKE2b digest (3 bytes rendered as 6 hex characters) over a canonical serialization of one part of the mesh: .. list-table:: :header-rows: 1 :widths: 15 85 * - Section - Covers * - ``P`` - Points: node ``x`` / ``y`` coordinates * - ``Z`` - Depth (bathymetry) * - ``E`` - Elements: face-node connectivity, in the canonical (encoding-independent) form; see :ref:`hybrid-connectivity` * - ``B`` - Boundary block * - ``A`` - Nodal attributes * - ``S`` - Self-attraction and loading Order independence ^^^^^^^^^^^^^^^^^^ Nodal-attribute names and SAL constituents are sorted before hashing. Two files that differ only in the on-disk order of their attributes or constituents therefore produce the same id. Where it appears ^^^^^^^^^^^^^^^^ The mesh id is used throughout Cocoa for provenance and consistency checking: - Cocoa computes the id when it reads the mesh and logs it at startup. - The id is stamped into every output NetCDF file as the ``mesh_id`` global attribute, so any result traces back to the exact mesh that produced it. - Checkpoints carry the id, and a restart compares it with the mesh it runs on (see :doc:`checkpoint`). - The domain-decomposition partition cache is keyed on the full id, so a partition file is only reused with the mesh it was built from. - Cocoa warns if a file's stored ``mesh_id`` does not match the id recomputed from its contents. Computing and verifying the id ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The ``checksum`` subcommand of ``cocoa-mesh-tools`` computes, writes, and verifies mesh ids: .. code-block:: bash # Print the id and per-section digests cocoa-mesh-tools checksum mesh.nc # Store the id in the file's mesh_id global attribute cocoa-mesh-tools checksum mesh.nc --write # Recompute and compare to the stored id (non-zero exit on mismatch) cocoa-mesh-tools checksum mesh.nc --verify # Emit the full result as JSON cocoa-mesh-tools checksum mesh.nc --json The ``checksum`` code in the ``cocoa_mesh_tools`` package is the reference implementation; Cocoa's C++ implementation produces the same ids. .. note:: The mesh-identifier concept was proposed by Alex Crosby (Oceanweather, Inc.). .. _nodal-attributes: Nodal Attributes ---------------- Cocoa supports spatially-varying physics parameters through nodal attributes stored in the mesh file. Available Nodal Attributes ^^^^^^^^^^^^^^^^^^^^^^^^^^ .. list-table:: :header-rows: 1 :widths: 35 15 15 35 * - Variable Name - Units - Default - Description * - ``mannings_n_at_sea_floor`` - s/m^(1/3) - 0.025 - Manning's n roughness coefficient for bottom friction * - ``primitive_weighting_in_continuity_equation`` - 1/s - 0.03 - GWCE weighting parameter (tau0) * - ``surface_submergence_state`` - (unitless) - 0.0 - Initial wet/dry state (0=wet, 1=dry, matches ADCIRC startDry) * - ``surface_directional_effective_roughness_length`` - m - 0.0 - Directional land roughness lengths for wind reduction (12 values per node) * - ``surface_canopy_coefficient`` - (unitless) - 1.0 - Binary canopy mask for wind sheltering (0=full sheltering, 1=no effect) * - ``elemental_slope_limiter`` - m/m - (disabled) - Per-node elevation-gradient threshold for the selective post-continuity elevation limiter (see :doc:`elemental_slope_limiter`). Read only when ``physics.elemental_slope_limiter: mesh``. Values must be non-negative; zero is the always-active sentinel. * - ``advection_state`` - (unitless) - 1.0 - Binary mask switching the advective terms on (1) or off (0) at each node. Read only when ``physics.advection: mesh``; fractional values fail loud. An element advects only when all of its corner nodes do, so a node set to 0 disables advection in every element touching it. .. note:: ADCIRC's fort.13 ``advection_state`` is a different quantity: it stores a per-node depth threshold, and advection is active in an element only where every corner satisfies ``depth >= threshold``. Cocoa stores the resolved flag instead, so the comparison happens once at mesh-conversion time rather than every element every timestep. ``cocoa-mesh-tools from_adcirc`` performs that conversion against the mesh bathymetry when it writes the nodal attributes; a fort.13 copied verbatim would present depths where a mask is expected and the binary validator would reject the mesh. Manning's n (mannings_n_at_sea_floor) """"""""""""""""""""""""""""""""""""" The Manning's n coefficient controls bottom friction. Spatially-varying values allow representation of different bottom types: - Sandy bottoms: 0.015-0.025 - Gravel: 0.025-0.035 - Vegetated areas: 0.035-0.070 - Channels with vegetation: 0.05-0.15 To use spatially varying Manning's n from the mesh file: .. code-block:: yaml physics: manning_n: mesh # Read from mesh file To use a constant value instead: .. code-block:: yaml physics: manning_n: 0.025 # Constant value for entire domain Tau0 (primitive_weighting_in_continuity_equation) """"""""""""""""""""""""""""""""""""""""""""""""" The tau0 parameter weights the GWCE between the primitive continuity equation and the wave equation form (see :doc:`/theory/gwce`): - Larger tau0: stronger weighting toward the primitive continuity equation - Smaller tau0: closer to the pure wave equation Typical values range from 0.001 to 0.03. Lower values are often used near open boundaries to reduce reflections. To use spatially varying tau0 from the mesh file: .. code-block:: yaml physics: tau0: mesh # Read from mesh file To use a constant value instead: .. code-block:: yaml physics: tau0: 0.005 # Constant value for entire domain Surface Submergence State (surface_submergence_state) """"""""""""""""""""""""""""""""""""""""""""""""""""" This attribute sets the initial wet/dry state of each node, matching ADCIRC's ``startDry`` nodal attribute convention. Use it for checkpoint restart simulations or when the initial water level would leave some areas dry. - 0 = Node starts wet (default) - 1 = Node starts dry When this attribute is present in the mesh file, Cocoa: 1. Checks depth-based drying (nodes with depth + initial_zeta <= h0) 2. Applies the submergence state (nodes with value=1 are set dry) 3. Initializes element active status based on nodal states 4. Checks for landlocked nodes (wet nodes with no active elements) 5. Sets dry node elevations to minimum depth (zeta = h0 - depth) If this attribute is absent, only depth-based and landlocked checks apply. Directional Effective Roughness (surface_directional_effective_roughness_length) """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" This attribute provides land surface roughness lengths in 12 directional bins (30-degree spacing, starting from East and proceeding counter-clockwise) for each node. It is used by the directional wind reduction algorithm to reduce wind velocity over land based on the upwind roughness characteristics. This attribute is stored as a 2D variable with shape ``[num_nodes, 12]``. - Units: meters - Typical values: 0.001 (open water) to 1.0+ (dense forest/urban) - Default: 0.0 (disabled) To use directional roughness from the mesh file: .. code-block:: yaml physics: surface_directional_roughness: mesh # Read 12-directional roughness from mesh Setting to ``0`` (default) disables directional wind reduction entirely. .. note:: A constant scalar value is not meaningful for this attribute because it requires 12 directional components per node. If a non-zero constant is specified in the configuration, it will be ignored with a warning. See :ref:`wind-reduction-config` for details on the wind reduction algorithm. Surface Canopy Coefficient (surface_canopy_coefficient) """"""""""""""""""""""""""""""""""""""""""""""""""""""" This attribute is a per-node binary canopy mask applied to the wind before stress computation. It marks whether free-stream wind reaches the water surface under forest canopy. - ``1`` = no sheltering (open area, default) - ``0`` = complete sheltering (no wind, hence no stress, at that node) Values must be exactly ``0`` or ``1``; a fractional value fails loud at startup. The per-source wind reduction folds the canopy factor into the wind before the quadratic drag law, which is exact only for a 0/1 mask. To use a spatially varying canopy mask from the mesh file: .. code-block:: yaml physics: surface_canopy_coefficient: mesh # Read from mesh file To use a constant value instead: .. code-block:: yaml physics: surface_canopy_coefficient: 0 # fully sheltered everywhere Setting to ``1`` (default) disables the canopy reduction. .. _converting-from-adcirc: Converting from ADCIRC Format ----------------------------- The ``cocoa-mesh-tools`` command converts ADCIRC fort.14 mesh files (and optionally fort.13 nodal attribute and fort.24 SAL files) to Cocoa's NetCDF format. It is part of the ``cocoa-tools`` Python package (``pip install utils/cocoa_tools``; installed in the containers); see :doc:`cocoa_tools`. .. note:: ``cocoa-mesh-tools`` is a multi-command tool. ``from_adcirc`` (shown below) builds a grouped Cocoa mesh from ADCIRC files; ``group_variables`` re-groups an existing mesh netCDF; ``reorder`` reorders an existing mesh for GPU performance (see :ref:`element-ordering`); and ``checksum`` computes, writes, or verifies a mesh's content id. Run ``cocoa-mesh-tools --help`` for the options of each. .. note:: ``from_adcirc`` does not support quadrilateral (SWAN-style ``NHY=4``) elements: it rejects a fort.14 with any such element instead of truncating it to a triangle. A mesh with quadrilaterals must be written directly as Cocoa NetCDF with ``max_face_nodes=4`` (see :ref:`hybrid-connectivity` above). Usage ^^^^^ **Basic conversion (mesh only):** .. code-block:: bash cocoa-mesh-tools from_adcirc \ --mesh fort.14 \ --output mesh.nc **With nodal attributes:** .. code-block:: bash cocoa-mesh-tools from_adcirc \ --mesh fort.14 \ --attributes fort.13 \ --output mesh.nc **With SAL data:** .. code-block:: bash cocoa-mesh-tools from_adcirc \ --mesh fort.14 \ --attributes fort.13 \ --sal fort.24 \ --output mesh.nc **With elements ordered for GPU performance** (recommended for production meshes; see :ref:`element-ordering`): .. code-block:: bash cocoa-mesh-tools from_adcirc \ --mesh fort.14 \ --attributes fort.13 \ --output mesh.nc \ --reorder ``--reorder all`` also renumbers the nodes first (see :ref:`node-renumbering`). **Command-line options:** .. list-table:: :header-rows: 1 :widths: 20 15 65 * - Option - Required - Description * - ``--mesh`` - Yes - Path to ADCIRC fort.14 mesh file * - ``--attributes`` - No - Path to ADCIRC fort.13 nodal attributes file * - ``--output`` - Yes - Path for output NetCDF file * - ``--sal`` - No - Path to ADCIRC SAL file (fort.24 ASCII or fort.24.nc NetCDF). Embeds per-node SAL amplitude and phase data into the mesh file. * - ``--reorder [elements-only|all]`` - No - Reorder the mesh after conversion (see :ref:`element-ordering`). ``--reorder`` alone, or ``--reorder elements-only``, sorts the elements by their lowest-numbered node: element numbering then differs from the fort.14, while nodes, boundaries, nodal attributes and SAL do not change. ``--reorder all`` first renumbers the nodes in reverse Cuthill-McKee order (see :ref:`node-renumbering`), so node numbering differs from the fort.14 and fort.13 as well. What Gets Converted ^^^^^^^^^^^^^^^^^^^ The conversion transfers: - Mesh geometry: node coordinates (x, y) and bathymetry (depth) - Element connectivity: triangle definitions - All boundary types: open boundaries, land boundaries, weirs, levees, etc. - Boundary attributes: crest elevations, flow coefficients, pipe parameters - Nodal attributes: all attributes from fort.13 (Manning's n, tau0, etc.) .. _element-ordering: Element Ordering ---------------- On a GPU, Cocoa processes elements in groups of 32 consecutive elements (one warp), and every element reads data stored at its nodes: depth, elevation, velocity, scale factors. When the 32 elements of a group share nearby node numbers, those reads land in a few contiguous memory blocks. When they are scattered across the mesh, each read is a separate memory transaction, and the element kernels, which are a large part of every time step, spend most of their time waiting on memory. Nodes in ADCIRC meshes are usually numbered to keep neighbors close (bandwidth or RCM ordering), but elements often are not: they arrive in whatever order the mesh generator or later editing left them. On the global STOFS mesh (12.8 million nodes), 32 consecutive elements touch nodes spread over a median of 26,000 node numbers, and the worst tenth over 180,000. Sorting the elements by their lowest-numbered node cuts that to 12,000 and 24,000, and on one V100 GPU it made every element kernel 1.5 to 2.4 times faster and the whole time step 16% faster, with the same results. Under MPI the effect is larger still, because each rank inherits its own slice of the file's element order: on two GPUs, the rank with the less favorable slice ran its element kernels 1.3 times slower per element, and the other rank waited for it. To sort the elements, convert with ``--reorder`` (above), or reorder an existing Cocoa mesh: .. code-block:: bash cocoa-mesh-tools reorder \ --mesh mesh.nc \ --output mesh_reordered.nc Both report how far 32 consecutive elements reach before and after, for example:: element order sorted by lowest node node numbers spanned by each 32 consecutive elements: median 26,482 -> 12,381, 90th percentile 179,588 -> 24,470 A large drop means the mesh benefits; little or no change means it was already well ordered. Every ``from_adcirc`` and ``group_variables`` run ends with a summary of the mesh it wrote (node and element counts, extent, depth range, boundaries by type, nodal attributes, SAL constituents, ordering, mesh id). When the elements were left unsorted and sorting them would shrink that span appreciably, it also prints a warning with the span it would reach and the command that sorts them:: warning: elements are not in GPU order. A block of 32 consecutive elements spans 40,510 node numbers (median; 118,966 at the 90th percentile), where sorted elements would span 4,208 (7,923), and every element kernel on the GPU slows with the spread. Convert with --reorder, or run: cocoa-mesh-tools reorder --mesh mesh.nc --output ``group_variables`` also warns about any source variable that fits none of the layout's groups, which it does not copy. Both commands write every nodal attribute they are given, but mark the ones Cocoa does not read (for example ``average_horizontal_eddy_viscosity_in_sea_water_wrt_depth`` or ``sea_surface_height_above_geoid`` from an ADCIRC fort.13) and warn about them: they have no effect on a simulation, yet they are part of the mesh id. Reordering changes only the order of the elements. Nodes, coordinates, depth, boundaries (including weirs and levees), nodal attributes and self-attraction and loading data are copied unchanged, and any other per-element variable is reordered with the elements. Because element numbering is part of the mesh identity, the reordered mesh has a new mesh id (only its topology section changes, and a stored ``mesh_id`` attribute is updated). As a consequence: - Partition cache files built for the original mesh do not apply; Cocoa builds a new one on the first parallel run. - Checkpoints written on the original mesh cannot be restarted on the reordered one; Cocoa refuses them at startup. Reorder before a simulation series, not in the middle of one. - Results agree with the original mesh to round-off, not exactly: the order in which element contributions are summed changes. Reordering rewrites variables inside the new file, and HDF5 does not reclaim the space the old data held, so the file is rewritten into a fresh copy as a last step. The report's ``storage`` line shows the size before and after. ``reorder`` writes a new file and refuses to overwrite its input. It also refuses a mesh that stores element numbers as values (UGRID ``face_face_connectivity`` and similar), which Cocoa meshes do not carry. .. _node-renumbering: Node Renumbering (optional) ^^^^^^^^^^^^^^^^^^^^^^^^^^^ Sorting elements by their lowest-numbered node only helps as much as the node numbering keeps neighbours close. ADCIRC meshes usually arrive bandwidth- or RCM-ordered, but a mesh that has been edited (nodes appended where it was refined, or removed where it was coarsened) or written by a generator that does not order its nodes can be far from it. For those, renumber the nodes in reverse Cuthill-McKee (RCM) order before the element sort: .. code-block:: bash cocoa-mesh-tools reorder --renumber-nodes \ --mesh mesh.nc \ --output mesh_reordered.nc or convert with ``from_adcirc --reorder all``. The report adds the node bandwidth and profile before and after, for example:: node order renumbered, reverse Cuthill-McKee bandwidth 101,782 -> 7,140, profile 1.617e+10 -> 1.426e+10 Each connected part of the mesh is ordered from every candidate start of a George-Liu pseudo-peripheral search, keeping the ordering with the smallest profile; on large meshes the start alone moves the bandwidth by tens of percent. Node renumbering is off by default because node numbers are what other files and workflows refer to. Every per-node variable (coordinates, depth, nodal attributes, SAL) is permuted, and the connectivity, boundary node lists and internal weir pair nodes are remapped, so the mesh itself is unchanged; but station files, partition caches, checkpoints, and anything outside the mesh file that names nodes by number no longer match. A mesh carrying other node connectivity (UGRID ``edge_node_connectivity`` or ``boundary_node_connectivity``) is refused, since those would be left stale. A mesh whose nodes are already well ordered gains little and can lose a little element locality, so compare the report against a plain ``reorder``. Adding Nodal Attributes ----------------------- You can add nodal attributes to an existing NetCDF mesh file using Python. Write them into the ``/nodal_attributes`` group: when that group exists, Cocoa reads nodal attributes only from it, and the mesh id covers only the attributes in it. .. code-block:: python import netCDF4 as nc import numpy as np # Open mesh file in append mode ds = nc.Dataset('mesh.nc', 'a') if 'nodal_attributes' in ds.groups: attrs = ds.groups['nodal_attributes'] else: attrs = ds.createGroup('nodal_attributes') # Add Manning's n attribute (spatially varying) if 'mannings_n_at_sea_floor' not in attrs.variables: manning = attrs.createVariable('mannings_n_at_sea_floor', 'f8', ('nnode',), fill_value=-9999.0, zlib=True, complevel=2) manning.long_name = "Manning's n at sea floor" manning.units = "s/m^(1/3)" # Set values (example: higher in shallow areas) depth = ds.variables['depth'][:] manning_values = np.where(depth < 10, 0.035, 0.025) attrs.variables['mannings_n_at_sea_floor'][:] = manning_values # Add tau0 attribute if 'primitive_weighting_in_continuity_equation' not in attrs.variables: tau0 = attrs.createVariable('primitive_weighting_in_continuity_equation', 'f8', ('nnode',), fill_value=-9999.0, zlib=True, complevel=2) tau0.long_name = "tau0 weighting parameter" tau0.units = "1/s" # Set constant tau0 attrs.variables['primitive_weighting_in_continuity_equation'][:] = 0.005 ds.close() Adding or changing an attribute changes the mesh id: existing partition caches no longer match the mesh, and a restart from an earlier checkpoint logs a warning naming the changed section. If the file stores a ``mesh_id`` attribute, rerun ``cocoa-mesh-tools checksum mesh.nc --write`` to update it; otherwise Cocoa warns at startup that the stored id does not match the file's contents. Mesh Generation Tools --------------------- Recommended tools for generating ADCIRC-format meshes (which can then be converted to Cocoa format): - SMS (Surface-water Modeling System) - OceanMesh2D (MATLAB/Octave mesh generator) - QGIS with mesh plugins - Triangle (open-source mesh generator) - Gmsh (open-source mesh generator) Complete NetCDF Format Reference -------------------------------- A complete CDL (NetCDF text format) example with every variable and attribute in the Cocoa mesh format: .. code-block:: text netcdf cocoa_mesh { dimensions: nnode = 8303; nface = 14761; max_face_nodes = 3; variables: // Mesh topology variable (UGRID convention) int cocoa_mesh; cocoa_mesh:cf_role = "mesh_topology"; cocoa_mesh:long_name = "Topology data of 2D unstructured mesh"; cocoa_mesh:topology_dimension = 2; cocoa_mesh:node_coordinates = "x y"; cocoa_mesh:face_node_connectivity = "mesh_face_nodes"; cocoa_mesh:face_dimension = "nface"; // Node coordinates double x(nnode); x:long_name = "x coordinate of mesh nodes"; x:standard_name = "longitude"; x:units = "degrees_east"; x:axis = "X"; x:_FillValue = -9999.0; double y(nnode); y:long_name = "y coordinate of mesh nodes"; y:standard_name = "latitude"; y:units = "degrees_north"; y:axis = "Y"; y:_FillValue = -9999.0; // Bathymetry double depth(nnode); depth:long_name = "z coordinate of mesh nodes"; depth:standard_name = "z"; depth:units = "meters"; depth:positive = "down"; depth:axis = "Z"; depth:_FillValue = -9999.0; // Element connectivity (0-based) int mesh_face_nodes(nface, max_face_nodes); mesh_face_nodes:long_name = "nodes of faces"; mesh_face_nodes:cf_role = "face_node_connectivity"; mesh_face_nodes:start_index = 0; mesh_face_nodes:_FillValue = -9999; // Boundary segment data (present only if the mesh has boundaries) group: boundaries { dimensions: nboundary = 45; nboundary_node = 1972; variables: int boundary_types(nboundary); boundary_types:long_name = "adcirc boundary type code"; boundary_types:standard_name = "boundary_types"; boundary_types:_FillValue = -9999; int boundary_start_index(nboundary); boundary_start_index:long_name = "start index of boundary nodes"; boundary_start_index:_FillValue = -9999; int boundary_end_index(nboundary); boundary_end_index:long_name = "end index of boundary nodes"; boundary_end_index:_FillValue = -9999; int boundary_nodes(nboundary_node); boundary_nodes:long_name = "boundary nodes data"; boundary_nodes:_FillValue = -9999; // The following are written only when they carry non-fill data int boundary_pair_node(nboundary_node); boundary_pair_node:long_name = "boundary pair node data"; boundary_pair_node:_FillValue = -9999; double boundary_elevation(nboundary_node); boundary_elevation:long_name = "boundary elevation above datum"; boundary_elevation:_FillValue = -9999.0; double boundary_supercritical_coefficient(nboundary_node); boundary_supercritical_coefficient:long_name = "boundary supercritical coefficient"; boundary_supercritical_coefficient:_FillValue = -9999.0; double boundary_subcritical_coefficient(nboundary_node); boundary_subcritical_coefficient:long_name = "boundary subcritical coefficient"; boundary_subcritical_coefficient:_FillValue = -9999.0; double boundary_pipe_coefficient(nboundary_node); boundary_pipe_coefficient:long_name = "boundary pipe coefficient"; boundary_pipe_coefficient:_FillValue = -9999.0; double boundary_pipe_height(nboundary_node); boundary_pipe_height:long_name = "boundary pipe height"; boundary_pipe_height:_FillValue = -9999.0; double boundary_pipe_diameter(nboundary_node); boundary_pipe_diameter:long_name = "boundary pipe diameter"; boundary_pipe_diameter:_FillValue = -9999.0; } // group boundaries // Nodal attributes (present only if the mesh has nodal attributes) group: nodal_attributes { dimensions: // Multi-value attributes declare a per-attribute na_length_ dimension na_length_12 = 12; variables: double mannings_n_at_sea_floor(nnode); mannings_n_at_sea_floor:long_name = "Manning's n at sea floor"; mannings_n_at_sea_floor:units = "s/m^(1/3)"; mannings_n_at_sea_floor:_FillValue = -9999.0; double primitive_weighting_in_continuity_equation(nnode); primitive_weighting_in_continuity_equation:long_name = "tau0 weighting"; primitive_weighting_in_continuity_equation:units = "1/s"; primitive_weighting_in_continuity_equation:_FillValue = -9999.0; double surface_submergence_state(nnode); surface_submergence_state:long_name = "Surface submergence state"; surface_submergence_state:_FillValue = -9999.0; double surface_directional_effective_roughness_length(nnode, na_length_12); surface_directional_effective_roughness_length:long_name = "Directional effective roughness length"; surface_directional_effective_roughness_length:units = "m"; surface_directional_effective_roughness_length:_FillValue = -9999.0; // na_length_12 = 12 (30-degree bins starting from East, counter-clockwise) double surface_canopy_coefficient(nnode); surface_canopy_coefficient:long_name = "Surface canopy coefficient"; surface_canopy_coefficient:units = "1"; surface_canopy_coefficient:_FillValue = -9999.0; } // group nodal_attributes // Self-attraction and loading data (present only if the mesh carries SAL) group: self_attraction_loading { dimensions: nsal_constituents = 8; sal_name_len = 10; variables: char sal_constituent_names(nsal_constituents, sal_name_len); sal_constituent_names:long_name = "SAL tidal constituent names"; double sal_frequency(nsal_constituents); sal_frequency:long_name = "SAL constituent angular frequency"; sal_frequency:units = "rad/s"; float sal_amplitude(nnode, nsal_constituents); sal_amplitude:long_name = "SAL amplitude"; sal_amplitude:units = "m"; float sal_phase(nnode, nsal_constituents); sal_phase:long_name = "SAL phase"; sal_phase:units = "degrees"; } // group self_attraction_loading // Global attributes :Conventions = "UGRID-1.0"; :title = "Cocoa mesh file"; :source = "Generated by cocoa-mesh-tools"; // Added by `cocoa-mesh-tools checksum --write`; see Mesh Identification. :mesh_id = "cm1_P4b5652_Z4d037a_Ec05a73_Bef8a11_A65085e_Scec7ea"; }