Fixed some precision things in vgrid that were causing edges to be duplicated. Added description of vgrid formats to SKILLS.md because they frequently cause confusion.

Author: Eli

SKILLS.md

@@ -0,0 +1,40 @@
+# Vertical Grid
+
+## vgrid.in Format vs schimpy Internal Representation (ivcor=1, LSC2)
+
+### On-disk format (vgrid.in, version 5.10)
+
+```
+1                                    ← ivcor (1 = LSC2 localized sigma)
+nvrt                                 ← max number of vertical levels across all nodes
+kbp(1) kbp(2) ... kbp(n_nodes)      ← 1-based bottom level index per node
+1  sigma(1,1) sigma(1,2) ...         ← level 1 values for all nodes (one row per level)
+2  sigma(2,1) sigma(2,2) ...
+...
+nvrt  sigma(nvrt,1) sigma(nvrt,2)...
+```
+
+- **Sigma range**: −1 (bottom) to 0 (surface), monotonically increasing with level index.
+- **Level ordering**: level 1 is the deepest possible; level `nvrt` is the surface.
+- **Per-node bottom**: `kbp(i)` (1-based) is the first valid level for node i. Levels below `kbp` are padded with **−9.0** (sentinel for "not applicable").
+- **Right-justified**: valid sigma values occupy levels `kbp` through `nvrt`; everything to the left is −9.0.
+- **Output precision**: `%14.6f` — only 6 decimal places. Consecutive sigma values that differ by less than 5×10⁻⁷ will round to the same value and cause SCHISM to abort ("Weird side (3)").
+
+### schimpy internal representation (`SchismLocalVerticalMesh`)
+
+- **Left-justified**: `sigma[i, 0:Ni]` holds the `Ni` valid values; `sigma[i, Ni:]` is `NaN`.
+- **kbps is 0-based**: `kbps[i] = nvrt − Ni`. File writes `kbps + 1`.
+- **Same sigma convention**: values are in [−1, 0], bottom-first to surface-last within the valid slice.
+- **Sentinel mapping**: on read, −9.0 → NaN. On write, NaN positions are filled with −9.0.
+
+### Pipeline internal convention (lsc2_v2.py)
+
+- **Sigma in [0, 1]**: 0 = surface, 1 = bottom. Left-justified: `sigma[i, 0] = 0`, `sigma[i, Ni-1] = 1`.
+- **Conversion to SCHISM format**: negate (`-sigma`) then `flip_sigma` (reverses non-NaN slice). Result: [−1, ..., 0, NaN, ...] — bottom-first, matching `SchismLocalVerticalMesh` expectations.
+- **Gotcha**: `flip_sigma` reverses the non-NaN portion. If the input is a palindrome (e.g., `[0, -1, 0]`), the flip is a no-op and the malformed array passes through unchanged. This happens when dry nodes (depth ≤ 0) have uninitialized interface values that clip to symmetric patterns.
+
+### Key invariants
+
+1. After the full pipeline, each node's sigma must be **strictly monotonically increasing** when written to disk.
+2. Monotonicity must survive the `%14.6f` rounding — consecutive sigma must differ by at least ~1×10⁻⁶.
+3. The z-coordinate check in SCHISM is at **side midpoints** (averaging two nodes), so both nodes of every edge must individually be monotone.

schimpy/lsc2_v2.py

@@ -575,6 +575,44 @@ def fix_sigma_pileups(
                 "depth_m": D,
             })
 
+    # ------------------------------------------------------------------
+    # Final pass: snap sigma to the output precision used by the vgrid
+    # writer (6 decimal places, i.e. %14.6f).  Without this, two sigma
+    # values that differ by less than 5e-7 survive the TOL=1e-10 check
+    # above but become identical after rounding in vgrid.in, which causes
+    # SCHISM to abort with "Weird side (3)".
+    # ------------------------------------------------------------------
+    OUTPUT_DECIMALS = 6
+    for i in range(n):
+        Ni = int(NL[i])
+        if Ni < 2:
+            continue
+        s = np.round(S[i, :Ni], OUTPUT_DECIMALS)
+        # Ensure endpoints are exact
+        s[0] = 0.0
+        s[-1] = 1.0
+        # Remove consecutive duplicates introduced by rounding
+        keep = np.r_[True, np.diff(s) > 0.0]
+        # Always keep the last entry (surface = 1.0)
+        keep[-1] = True
+        if not keep.all():
+            s = s[keep]
+            S[i, :] = np.nan
+            S[i, :len(s)] = s
+            NL[i] = len(s)
+            logs.append({
+                "node": i,
+                "x": float(mesh.nodes[i, 0]) if have_xy else np.nan,
+                "y": float(mesh.nodes[i, 1]) if have_xy else np.nan,
+                "Ni_before": Ni,
+                "Ni_after": len(s),
+                "dupes_removed": int(Ni - len(s)),
+                "bottom_drops": 0,
+                "dz1_after_m": np.nan,
+                "dz2_after_m": np.nan,
+                "depth_m": float(depth[i]),
+            })
+
     df = pd.DataFrame(logs, columns=[
         "node","x","y","Ni_before","Ni_after",
         "dupes_removed","bottom_drops","dz1_after_m","dz2_after_m","depth_m"
@@ -1333,7 +1371,14 @@ def _dbg_levels(tag):
     # Finalization: single forward guard + exact endpoints (+ strict nudge)
     for i in range(n):
         Ni = int(Nlevels[i]); D = depth[i]
-        if Ni < 2 or D <= 0.0:
+        if Ni < 2:
+            continue
+
+        # Dry / non-positive-depth nodes: assign uniform degenerate sigma
+        if D <= 0.0:
+            D_eff = max(D, 1e-6)
+            for k in range(Ni):
+                h[i, k] = D_eff * k / (Ni - 1)
             continue
 
         # Surface exactly 0