==================
Momentum Equations
==================

Cocoa solves the non-conservative form of the 2D depth-integrated momentum
equations in spherical coordinates.

Governing Equations
-------------------

The momentum equations in spherical coordinates on an equidistant cylindrical
projection (:math:`\lambda, \phi` for longitude and latitude):

.. math::

   \begin{split}
   \frac{\partial U}{\partial t}
   &+ \frac{U}{R\cos\phi}\frac{\partial U}{\partial \lambda}
    + \frac{V}{R}\frac{\partial U}{\partial \phi}
    - f_{\text{eff}} V \\
   &= -\frac{g}{R\cos\phi}\frac{\partial \zeta}{\partial \lambda}
    + \frac{\tau_{sx} - \tau_{bx}}{\rho_0 H} + M_x
   \end{split}

.. math::

   \begin{split}
   \frac{\partial V}{\partial t}
   &+ \frac{U}{R\cos\phi}\frac{\partial V}{\partial \lambda}
    + \frac{V}{R}\frac{\partial V}{\partial \phi}
    + f_{\text{eff}} U \\
   &= -\frac{g}{R}\frac{\partial \zeta}{\partial \phi}
    + \frac{\tau_{sy} - \tau_{by}}{\rho_0 H} + M_y
   \end{split}

where :math:`M_x` and :math:`M_y` represent lateral stress terms, and the
effective Coriolis parameter includes a spherical metric correction:

.. math::

   f_{\text{eff}} = 2\Omega\sin\phi + \frac{U \tan\phi}{R}

The second term accounts for the convergence of meridians and is evaluated
using the velocity from the previous timestep.

Coordinate Systems
------------------

Cocoa supports multiple cylindrical map projections for spherical coordinates.
The default is the equidistant cylindrical (CPP) projection, corresponding to
ADCIRC's ICS=21 setting. See :doc:`coordinate_systems` for complete details
on available projections and their scale factors.

**Key points for momentum assembly:**

- The mesh precomputes shape function derivatives in projected coordinates
- During assembly, x-derivatives are multiplied by the momentum scale factor
  :math:`S_x^{mom}` (stored in ``elem_scale_momentum_x``)
- The y-derivatives require no correction for the equidistant cylindrical
  projection since latitude spacing is uniform
- The effective Coriolis parameter includes a spherical metric correction
  using :math:`\tan\phi/R` (see :ref:`equidistant-cylindrical`)

Time Discretization
-------------------

Bottom friction and Coriolis use Crank-Nicolson weighting (:math:`\theta = 0.5`).
Writing the x-momentum equation:

.. math::

   \begin{split}
   \frac{U^{n+1} - U^n}{\Delta t}
   - f_{\text{eff}}\left(\frac{V^{n+1} + V^n}{2}\right)
   &= -\tau_{xx}\left(\frac{U^{n+1} + U^n}{2}\right) \\
   &\quad - \tau_{xy}\left(\frac{V^{n+1} + V^n}{2}\right) + F_x
   \end{split}

Rearranging:

.. math::

   \begin{split}
   \left(1 + \frac{\Delta t}{2}\tau_{xx}\right) U^{n+1}
   &+ \frac{\Delta t}{2}\left(\tau_{xy} - f_{\text{eff}}\right) V^{n+1} \\
   &= \left(1 - \frac{\Delta t}{2}\tau_{xx}\right) U^n
    - \frac{\Delta t}{2}\left(\tau_{xy} - f_{\text{eff}}\right) V^n
    + \Delta t \cdot F_x
   \end{split}

Per-Node 2x2 System
----------------------------------

.. graphviz::
   :align: center
   :caption: Momentum solver: per-node 2x2 system assembly and solve

   digraph momentum {
     rankdir=LR;
     node [shape=box, style="filled,rounded", fontname="Helvetica", fontsize=10];
     edge [color="#555555"];
     bgcolor="transparent";

     // RHS contributions (inputs)
     subgraph cluster_rhs {
       label="RHS Contributions (element loop)";
       style="filled,rounded";
       fillcolor="#F5F5F5";
       color="#BDBDBD";
       fontname="Helvetica";
       fontsize=10;

       pres [label="Pressure\nGradient", fillcolor="#E3F2FD", color="#1E88E5"];
       adv [label="Advection", fillcolor="#E3F2FD", color="#1E88E5"];
       lat [label="Lateral\nStress", fillcolor="#E3F2FD", color="#1E88E5"];
       tipg [label="Tide Potential\nGradient", fillcolor="#E3F2FD", color="#1E88E5"];
       wind [label="Wind Stress\n(nodal)", fillcolor="#E3F2FD", color="#1E88E5"];
     }

     // LHS contributions
     subgraph cluster_lhs {
       label="LHS Matrix (per node)";
       style="filled,rounded";
       fillcolor="#F5F5F5";
       color="#BDBDBD";
       fontname="Helvetica";
       fontsize=10;

       fric [label="Bottom Friction\n(TKM diagonal)", fillcolor="#FFF3E0", color="#FB8C00"];
       cor [label="Coriolis\n(off-diagonal)", fillcolor="#FFF3E0", color="#FB8C00"];
     }

     // Assembly
     assemble [label="Assemble 2x2\nSystem per Node", fillcolor="#E8F5E9", color="#43A047"];

     // Solve
     solve [label="Cramer's Rule\nDirect Solve", fillcolor="#FFEBEE", color="#E53935"];

     // Output
     vel [label="U, V\n(velocity)", fillcolor="#F3E5F5", color="#8E24AA"];
     flux [label="Qx, Qy\n(flux = H * U)", fillcolor="#F3E5F5", color="#8E24AA"];

     pres -> assemble;
     adv -> assemble;
     lat -> assemble;
     tipg -> assemble;
     wind -> assemble;
     fric -> assemble;
     cor -> assemble;
     assemble -> solve -> vel;
     solve -> flux;
   }

The momentum solve at each node:

.. math::
   :label: momentum-2x2

   \begin{bmatrix}
   1 + \frac{\Delta t}{2}\tau_{xx} & \frac{\Delta t}{2}(\tau_{xy} - f_{\text{eff}}) \\[6pt]
   \frac{\Delta t}{2}(\tau_{xy} + f_{\text{eff}}) & 1 + \frac{\Delta t}{2}\tau_{yy}
   \end{bmatrix}
   \begin{bmatrix} U^{n+1} \\ V^{n+1} \end{bmatrix}
   =
   \begin{bmatrix} R_x \\ R_y \end{bmatrix}

The right-hand side:

.. math::

   \begin{split}
   R_x = \left(1 - \frac{\Delta t}{2}\tau_{xx}\right) U^n
       &- \frac{\Delta t}{2}(\tau_{xy} - f_{\text{eff}}) V^n \\
       &+ \Delta t \cdot F_x
   \end{split}

.. math::

   \begin{split}
   R_y = \left(1 - \frac{\Delta t}{2}\tau_{yy}\right) V^n
       &- \frac{\Delta t}{2}(\tau_{xy} + f_{\text{eff}}) U^n \\
       &+ \Delta t \cdot F_y
   \end{split}

Coriolis appears only in the off-diagonal entries of Equation :eq:`momentum-2x2`.
The system is solved via Cramer's rule, keeping all four matrix entries in GPU
registers.

Bottom Friction (TKM Tensor)
----------------------------

The TKM tensor stores the linearized friction coefficient. For Manning's n
(see :eq:`manning-friction`):

.. math::
   :label: tkm-cd

   C_d = \frac{g n^2}{H^{1/3}}

.. math::
   :label: tkm-tau

   \tau = \min\left(\frac{C_d |\mathbf{U}|}{H}, \frac{2}{\Delta t}\right)

The upper bound :math:`2/\Delta t` prevents the implicit friction term from
overwhelming the solution in very shallow water.

For isotropic friction, :math:`\tau_{xx} = \tau_{yy} = \tau` and
:math:`\tau_{xy} = 0`. When internal tide friction is enabled, its
contribution is added to the TKM tensor after Manning's friction is computed
(see :ref:`internal-tide-friction` in the governing equations). With tensor
internal tide friction, :math:`\tau_{xy} \neq 0` and the full
:math:`2 \times 2` system in Equation :eq:`momentum-2x2` is required.

These values are stored directly in ``PhysicsData``
and used as shown in the matrix above.

RHS Forcing Terms (Continuous)
------------------------------

The forcing terms :math:`F_x` and :math:`F_y` include:

**Pressure Gradient** (time-centered):

.. math::

   \begin{aligned}
   F_x^{pres} &= -\frac{g}{2} \alpha_L
   \left(\frac{\partial \zeta^{n-1}}{\partial x}
   + \frac{\partial \zeta^{n+1}}{\partial x}\right) \\[4pt]
   F_y^{pres} &= -\frac{g}{2} \alpha_L
   \left(\frac{\partial \zeta^{n-1}}{\partial y}
   + \frac{\partial \zeta^{n+1}}{\partial y}\right)
   \end{aligned}

where :math:`\alpha_L` is the slope limiter from wet/dry processing
(Equation :eq:`slope-limiter`) that reduces pressure gradients in shallow water.

**Tide Potential Gradient:**

.. math::

   \begin{aligned}
   F_x^{tip} &= \frac{g}{2} \alpha_L
   \left(\frac{\partial \eta_{tip}^{n-1}}{\partial x}
   + \frac{\partial \eta_{tip}^{n+1}}{\partial x}\right) \\[4pt]
   F_y^{tip} &= \frac{g}{2} \alpha_L
   \left(\frac{\partial \eta_{tip}^{n-1}}{\partial y}
   + \frac{\partial \eta_{tip}^{n+1}}{\partial y}\right)
   \end{aligned}

The tide potential acts as a body force opposing the free surface gradient.

**Advection** (non-conservative form, time-weighted):

.. math::

   F_x^{adv} = -\left(\bar{U} \frac{\partial U}{\partial x}
   + \bar{V} \frac{\partial U}{\partial y}\right), \quad
   F_y^{adv} = -\left(\bar{U} \frac{\partial V}{\partial x}
   + \bar{V} \frac{\partial V}{\partial y}\right)

Cocoa uses 2-level time centering with weights between :math:`n-1` and :math:`n`.

**Wind Stress** (nodal application):

.. math::

   F_x^{wind} = \frac{f_w \, \tau_{sx}}{\rho_0 H}, \quad
   F_y^{wind} = \frac{f_w \, \tau_{sy}}{\rho_0 H}

where :math:`f_w` is the depth-dependent wind stress limiter that smoothly
ramps wind forcing to zero near the drying threshold. Wind stress is computed
from 10-meter wind velocity using the Garratt drag law and applied per node
(not per element). See :doc:`meteorological_forcing` for the drag law
formulation, wind limiter definition, and discretized form.

**Lateral Stress Divergence:**

.. math::

   F_x^{ls} = \frac{1}{H}\left(\frac{\partial \sigma_{xx}}{\partial x}
   + \frac{\partial \sigma_{xy}}{\partial y}\right), \quad
   F_y^{ls} = \frac{1}{H}\left(\frac{\partial \sigma_{yx}}{\partial x}
   + \frac{\partial \sigma_{yy}}{\partial y}\right)

The stress tensor uses the Smagorinsky eddy viscosity model.

Discretized RHS Assembly
------------------------

The forcing terms are discretized using linear triangular finite elements.
For an element with nodes :math:`(1, 2, 3)`, area :math:`A_e`, and shape function
derivatives stored as :math:`FDX_i = 2 A_e \frac{\partial N_i}{\partial x}`.

**Pressure Gradient (Discretized):**

.. math::

   \begin{aligned}
   F_x^{pres,e} &= -\frac{g \Delta t}{4} \alpha_L
   \sum_{i=1}^{3} (\zeta_i^{n-1} + \zeta_i^{n+1}) \cdot FDX_i \\[4pt]
   F_y^{pres,e} &= -\frac{g \Delta t}{4} \alpha_L
   \sum_{i=1}^{3} (\zeta_i^{n-1} + \zeta_i^{n+1}) \cdot FDY_i
   \end{aligned}

Note: The :math:`\Delta t` factor is included because the RHS represents
:math:`\Delta t \cdot F`, not :math:`F` alone. The factor of :math:`1/4`
combined with :math:`FDX = 2A_e \partial N/\partial x` (and similarly
:math:`FDY = 2A_e \partial N/\partial y`) yields the continuous
:math:`g/2` coefficient for the time-centered average.

**Advection (Discretized):**

The velocity gradients are computed as:

.. math::

   \frac{\partial U}{\partial x} = \frac{1}{2A_e}
   \sum_{i=1}^{3} U_i \cdot FDX_i, \quad
   \frac{\partial V}{\partial x} = \frac{1}{2A_e}
   \sum_{i=1}^{3} V_i \cdot FDX_i

and similarly for the y-derivatives using :math:`FDY_i`.

The advection contribution uses element-averaged velocities:

.. math::

   \begin{aligned}
   F_x^{adv,e} &= -\Delta t \left(\bar{U} \frac{\partial U}{\partial x}
   + \bar{V} \frac{\partial U}{\partial y}\right) \\[4pt]
   F_y^{adv,e} &= -\Delta t \left(\bar{U} \frac{\partial V}{\partial x}
   + \bar{V} \frac{\partial V}{\partial y}\right)
   \end{aligned}

where :math:`\bar{U} = \frac{1}{3} \, \sum_{i=1}^{3} U_i`. Cocoa applies time-weighting
between levels :math:`n-1` and :math:`n` with coefficients matching ADCIRC's
2-level centered scheme.

**Lateral Stress (Discretized):**

The Smagorinsky eddy viscosity is computed per-element using the deviatoric strain rate
(matching ADCIRC's formulation):

.. math::

   \nu_e = C_s \cdot A_e \cdot |S|, \quad
   |S| = \sqrt{
   \left(\frac{\partial U}{\partial x}
   - \frac{\partial V}{\partial y}\right)^{\!2}
   + \left(\frac{\partial U}{\partial y}
   + \frac{\partial V}{\partial x}\right)^{\!2}
   }

The stress tensor uses per-node depth :math:`H_i`. Cocoa uses ADCIRC's
symmetric IBPSV formulation (``CME_LS_IBPSV``), so the off-diagonal shear
components are equal:

.. math::

   \sigma_{xx,i} = \nu_e H_i \frac{\partial U}{\partial x}, \quad
   \sigma_{yy,i} = \nu_e H_i \frac{\partial V}{\partial y}

.. math::

   \sigma_{xy,i} = \sigma_{yx,i}
   = \tfrac{1}{2} \nu_e H_i
   \left(\frac{\partial U}{\partial y} + \frac{\partial V}{\partial x}\right)

The lateral stress divergence contribution (symmetric IBPSV formulation):

.. math::

   F_x^{ls,i} = -\frac{\Delta t}{H_i} \cdot 1.5 \cdot
   (\sigma_{xx,i} \cdot FDX_i + \sigma_{xy,i} \cdot FDY_i)

.. math::

   F_y^{ls,i} = -\frac{\Delta t}{H_i} \cdot 1.5 \cdot
   (\sigma_{yx,i} \cdot FDX_i + \sigma_{yy,i} \cdot FDY_i)

**Element-to-Node Assembly:**

Element contributions are accumulated to nodes using atomic operations (shown for
x-momentum; y-momentum follows the same pattern with :math:`F_y` terms):

.. math::

   \begin{aligned}
   RHS_{x,i} &\mathrel{+}= F_x^{pres,e} + F_x^{adv,e} + F_x^{ls,i} \\
   RHS_{y,i} &\mathrel{+}= F_y^{pres,e} + F_y^{adv,e} + F_y^{ls,i}
   \end{aligned}

After assembly, the accumulated RHS is normalized by the node's control volume:

.. math::

   RHS_{x,i} \leftarrow \frac{RHS_{x,i}}{A_{total,i}}, \quad
   RHS_{y,i} \leftarrow \frac{RHS_{y,i}}{A_{total,i}}

where :math:`A_{total,i} = \sum_{e \ni i} A_e/3` is the sum of one-third of each
element area surrounding node :math:`i`

Land Boundary Conditions
------------------------

At land boundaries, Cocoa enforces a no-normal-flow condition:

.. math::

   n_x U + n_y V = 0

This is implemented by transforming the :math:`2 \times 2` system to tangential-normal coordinates:

1. Transform RHS: :math:`(RHS_t, RHS_n) = R \cdot (RHS_x, RHS_y)`
2. Transform matrix to rotated coordinates
3. Set :math:`RHS_n = 0` (no normal flow)
4. Solve transformed system
5. Transform back to :math:`(U, V)`

The transformation uses:

- :math:`\cos\theta = n_x` (CSII in ADCIRC notation)
- :math:`\sin\theta = n_y` (SIII in ADCIRC notation)

Flux Computation
----------------

After solving for velocity, the flux is computed:

.. math::

   Q_x = H \cdot U, \quad Q_y = H \cdot V

where :math:`H = h + \zeta` is the total water depth.

A minimum depth threshold is applied to avoid division issues in very shallow water.

Implementation Files
--------------------

- ``MomentumSolver.hpp``: High-level solver class
- ``MomentumSolveKernels.hpp``: Per-node :math:`2 \times 2` solve kernel
- ``MomentumRhsKernels.hpp``: RHS assembly kernels
- ``BottomFriction.hpp``: TKM tensor computation
- ``LateralStress.hpp``: Smagorinsky turbulence model