================ Mesh Preparation ================ Cocoa uses unstructured triangular meshes stored in NetCDF format following UGRID conventions. This document describes the mesh format and how to prepare meshes for Cocoa simulations. Mesh Requirements ----------------- Cocoa meshes must satisfy the following requirements: - **Element Type**: Linear triangular elements only - **Connectivity**: Properly connected mesh with no isolated nodes - **Quality**: Reasonable element aspect ratios and low (<8) valency - **Numbering**: Consistent counter-clockwise node ordering - **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** directly read ADCIRC fort.14 files. Use the provided ``cocoa_mesh_tools.py`` script 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 (always 3 for triangles) * - ``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 (0-based node indices, start_index=0) **/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 ------------------- Every Cocoa mesh has a **mesh id**: a content-based fingerprint that identifies the mesh by its actual contents. The id is a short, filename-safe string. 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 * - ``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. - The domain-decomposition partition cache is keyed on the full id, so a partition file is provably paired 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.py`` computes, writes, and verifies mesh ids: .. code-block:: bash # Print the id and per-section digests python utils/cocoa_mesh_tools.py checksum mesh.nc # Store the id in the file's mesh_id global attribute python utils/cocoa_mesh_tools.py checksum mesh.nc --write # Recompute and compare to the stored id (non-zero exit on mismatch) python utils/cocoa_mesh_tools.py checksum mesh.nc --verify # Emit the full result as JSON python utils/cocoa_mesh_tools.py checksum mesh.nc --json The ``checksum`` code in ``cocoa_mesh_tools.py`` is the reference implementation; Cocoa's C++ implementation produces byte-identical 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. These allow different parameter values at each node, enabling more realistic representation of physical processes. 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 - Canopy coefficient for wind sheltering (0=full sheltering, 1=no effect) 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 **Configuration:** 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 controls the GWCE formulation weighting between the primitive continuity equation and the wave equation form: - tau0 = 0: Pure primitive continuity equation (mass conservation) - tau0 > 0: Wave equation form (improved stability) Typical values range from 0.001 to 0.03. Lower values are often used near open boundaries to reduce reflections. **Configuration:** 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 specifies the initial wet/dry state of each node, matching ADCIRC's ``startDry`` nodal attribute convention. This is useful for checkpoint restart simulations or when the initial water level would leave some areas dry. - 0 = Node starts wet (default) - 1 = Node starts dry **Behavior:** When this attribute is present in the mesh file, Cocoa will: 1. First check depth-based drying (nodes with depth + initial_zeta <= h0) 2. Then apply submergence state (nodes with value=1 are set dry) 3. Initialize element active status based on nodal states 4. Check for landlocked nodes (wet nodes with no active elements) 5. Set 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]``. - 12 components per node: one roughness length per 30-degree directional bin - Units: meters - Typical values: 0.001 (open water) to 1.0+ (dense forest/urban) - Default: 0.0 (disabled) **Configuration:** 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 specifies a per-node multiplicative canopy sheltering factor applied to wind velocity before stress computation. It represents the fraction of free-stream wind that reaches the water surface under forest canopy. - 1.0 = no sheltering (open area, default) - 0.0 = complete sheltering (no wind reaches surface) - Typical forested values: 0.3--0.7 **Configuration:** To use spatially-varying canopy coefficients 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.5 # 50% reduction everywhere Setting to ``1.0`` (default) disables the canopy reduction. .. _converting-from-adcirc: Converting from ADCIRC Format ----------------------------- Cocoa provides a Python script to convert ADCIRC fort.14 mesh files (and optionally fort.13 nodal attribute files) to Cocoa's NetCDF format. Script Location ^^^^^^^^^^^^^^^ The conversion script is located at: .. code-block:: text utils/cocoa_mesh_tools.py .. note:: ``cocoa_mesh_tools.py`` 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; and ``checksum`` computes, writes, or verifies a mesh's content id. Run ``cocoa_mesh_tools.py --help`` for the options of each. Usage ^^^^^ **Basic conversion (mesh only):** .. code-block:: bash python utils/cocoa_mesh_tools.py from_adcirc \ --mesh fort.14 \ --output mesh.nc **With nodal attributes:** .. code-block:: bash python utils/cocoa_mesh_tools.py from_adcirc \ --mesh fort.14 \ --attributes fort.13 \ --output mesh.nc **With SAL data:** .. code-block:: bash python utils/cocoa_mesh_tools.py from_adcirc \ --mesh fort.14 \ --attributes fort.13 \ --sal fort.24 \ --output mesh.nc **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. **Requirements:** The script requires Python 3 with the following packages: - ``netCDF4`` - ``numpy`` Install with: .. code-block:: bash pip install netCDF4 numpy What Gets Converted ^^^^^^^^^^^^^^^^^^^ The conversion script 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.) Adding Nodal Attributes ----------------------- You can add nodal attributes to an existing NetCDF mesh file using Python: .. code-block:: python import netCDF4 as nc import numpy as np # Open mesh file in append mode ds = nc.Dataset('mesh.nc', 'a') num_nodes = len(ds.dimensions['nnode']) # Add Manning's n attribute (spatially varying) if 'mannings_n_at_sea_floor' not in ds.variables: manning = ds.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) ds.variables['mannings_n_at_sea_floor'][:] = manning_values # Add tau0 attribute if 'primitive_weighting_in_continuity_equation' not in ds.variables: tau0 = ds.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 ds.variables['primitive_weighting_in_continuity_equation'][:] = 0.005 ds.close() 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 -------------------------------- Below is a complete CDL (NetCDF text format) example showing all variables and attributes 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.py"; // Added by `cocoa_mesh_tools.py checksum --write`; see Mesh Identification. :mesh_id = "cm1_P4b5652_Z4d037a_Ec05a73_Bef8a11_A65085e_Scec7ea"; }