Architecture

This document describes the software architecture of Cocoa.

High-Level Overview

Cocoa follows a modular architecture with clear separation of concerns, built on the Kokkos [Edwards2014] [Trott2022] performance portability framework within the Trilinos [Heroux2005] ecosystem:

Fig. 34 High-level architecture showing the simulation pipeline

Source Library Dependencies

The source libraries have a clear dependency hierarchy:

Fig. 35 Library dependencies (solid = internal, dashed = external)

Note

Namespace map. Cocoa::IO spans two homes: the cocoa_io library (YAML parsing via ConfigurationReader, NetCDF read/write) and cocoa_kernel’s own io/ and core/config/ directories (mesh I/O, checkpointing, output management, and YAML -> config translation). Cocoa::Types lives in one home, cocoa_compute_types (Vec2, Mat2, Mat3, MathUtils, CoriolisParameter, and the scalar/index type aliases in Ordinals.hpp). When hunting for a symbol in Cocoa::IO, check both locations rather than assuming a single library owns it.

Data Container Hierarchy

ModelFields is the centralized data container passed through the simulation:

Fig. 36 ModelFields data container hierarchy

GWCE Solver Selection

The continuity equation (GWCE) is solved by one of two interchangeable implementations behind GwceSolver, selected at construction time by numeric.solver in the YAML configuration:

• Lumped (lumped): a diagonal mass matrix, solved directly with no iteration. This is the explicit path.

• Consistent (consistent): a sparse, symmetric positive-definite matrix, solved iteratively with a preconditioned Belos conjugate-gradient solver. This is the semi-implicit path.

The momentum solve is always explicit regardless of which GWCE solver is selected; only the continuity equation has an implicit option. Both GWCE variants reassemble their LHS matrix every step, so there is no wet/dry change flag to thread between assembly and solve.

Directory Structure

src/
├── cocoa/                        # Main application
│   ├── cocoa.cpp                 # Entry point
│   └── CommandLineArgs.hpp       # CLI argument parsing
│
├── cocoa_datetime/               # Date/time library
│   ├── DateTime.hpp              # Date/time representation
│   └── TimeDelta.hpp             # Time duration
│
├── cocoa_compute_types/          # Small math value types (Cocoa::Types)
│   ├── Vec2.hpp                  # 2D vector
│   ├── Mat2.hpp                  # 2x2 matrix
│   ├── Mat3.hpp                  # 3x3 matrix
│   ├── MathUtils.hpp             # Math helpers
│   └── Ordinals.hpp              # Scalar/index types
│
├── cocoa_constants/              # Shared physical/numerical constants
│   ├── Constants.hpp             # Physical constants
│   ├── Defaults.hpp              # Default parameter values
│   └── Thresholds.hpp            # Centralized numerical thresholds
│
├── cocoa_io/                     # I/O utilities library
│   ├── ConfigurationReader.hpp   # YAML config parsing
│   ├── Logger.hpp                # Logging facade
│   ├── NetcdfReader.hpp          # NetCDF input
│   ├── NetcdfWriter.hpp          # NetCDF output
│   └── NetcdfCommon.hpp          # Shared NetCDF types
│
├── cocoa_meteo/                   # Meteorological I/O library
│   ├── core/                     # Base classes and data types
│   │   ├── MeteoReaderBase.hpp   # Abstract reader interface
│   │   ├── MeteoReaderConfig.hpp # Format-agnostic configuration
│   │   ├── MeteoFormat.hpp       # Format enum (CF/OWI ASCII/OWI NetCDF/GRIB)
│   │   ├── MeteoGrid.hpp         # Regular/irregular grid types
│   │   ├── MeteoField.hpp        # 2D field container
│   │   ├── MeteoFieldSnapshot.hpp # Time snapshot bundle
│   │   ├── GribConfig.hpp        # Numeric GRIB2 field identity + presets
│   │   ├── LongitudeUnwrap.hpp   # Antimeridian-safe longitude shift
│   │   ├── TimeAxisUnion.hpp     # Sorted-unique union of time axes
│   │   ├── CfTimeUtils.hpp       # CF time unit parsing
│   │   └── CfPacking.hpp         # CF scale/offset/fill-value unpacking
│   ├── readers/                  # Format-specific readers
│   │   ├── cf/CfNetcdfReader.hpp
│   │   ├── owi_ascii/{OwiAsciiReader,OwiAsciiDomain}.hpp
│   │   ├── owi_netcdf/{OwiNetcdfReader,OwiNetcdfDomain}.hpp
│   │   └── grib/                 # GRIB2 via ecCodes (cocoa_ENABLE_GRIB)
│   │       ├── GribReader.hpp    # One domain per configured weather product
│   │       ├── GribDomain.hpp    # Startup-built index + grid per priority level
│   │       ├── GribFileIndex.hpp # Message lookup by (product, valid time)
│   │       ├── GribGridBuilder.hpp # Grid + row-order normalization
│   │       ├── GribHandle.hpp    # RAII ecCodes handle
│   │       └── GribProducts.hpp  # Preset numeric identities (GFS, HRRR, ...)
│   ├── interpolation/            # Spatial interpolation
│   │   ├── MeteoSpatialInterpolator.hpp  # Single-domain bilinear weights
│   │   ├── MultiDomainInterpolator.hpp   # Multi-domain combine, winner mask, activation weights
│   │   ├── FieldInterpolation.hpp        # Apply precomputed weights to a field
│   │   └── InterpolationWeight.hpp       # SoA bilinear weight storage
│   ├── CompositeMeteoReader.hpp   # Composes multiple sources into one flat domain space
│   └── MeteoReaderFactory.hpp    # Reader instantiation (single reader or composite)
│
├── cocoa_vortex/                 # GAHM parametric hurricane vortex
│   ├── core/                     # ATCF track parsing + parameter table
│   │   ├── AtcfTrack.hpp         # Immutable parsed track (snapshots, isotachs)
│   │   ├── AtcfParser.hpp        # ATCF b-deck text parser
│   │   ├── GahmParameterTable.hpp # SoA device parameter table
│   │   └── GahmHostEvaluator.hpp # Host reference evaluator (shares physics/ with device)
│   ├── physics/                  # Shared host/device closed-form math
│   │   ├── GahmLimits.hpp        # Solver bracket/clamp limits
│   │   ├── GahmPhysics.hpp       # KOKKOS_INLINE_FUNCTION GAHM profile math
│   │   ├── PressureEstimateConfig.hpp # Wind-pressure relationship selection
│   │   └── WindPressure.hpp      # Missing-pressure fill from wind speed
│   └── solver/                   # Once-at-startup host solve
│       ├── GahmModelConfig.hpp   # GAHM2026 / background-wind model config
│       ├── GahmSolver.hpp        # Rmax/Bg fixed-point solver + outcome status
│       └── GahmPreprocessor.hpp  # Track -> parameter table preprocessing
│
└── cocoa_kernel/                 # Core computational library
    ├── core/                     # Domain-meaningful primitives
    │   ├── Execution.hpp         # Execution space selection
    │   ├── KokkosProfileRegion.hpp # Performance profiling
    │   ├── ModelConfiguration.hpp # Configuration struct
    │   ├── NarrowCast.hpp        # Checked numeric casts
    │   ├── algorithms/
    │   │   └── StreamCompactor.hpp # Wet-list stream compaction backends
    │   ├── config/               # Per-domain config sub-structs
    │   │   ├── FlowBoundaryConfig.hpp
    │   │   ├── ForcingConfig.hpp
    │   │   ├── GwceConfig.hpp
    │   │   ├── CheckpointConfig.hpp
    │   │   ├── CocoaConfigSchema.hpp # Accessor for the embedded config schema
    │   │   ├── cocoa_config.schema.yaml # Config schema (JSON Schema draft-07)
    │   │   ├── ForcingParsers.hpp    # Forcing-section parsers
    │   │   ├── MeshConfig.hpp
    │   │   ├── ModelConfigurationFactory.hpp # YAML -> config
    │   │   ├── OutputConfig.hpp
    │   │   ├── PhysicsConfig.hpp
    │   │   ├── SimulationConfig.hpp
    │   │   └── TidalConfig.hpp
    │   └── types/                # Generic infrastructure types
    │       ├── KokkosAliases.hpp # View type aliases
    │       ├── LinearAlgebraTypes.hpp # Trilinos type aliases
    │       ├── Precision.hpp     # Float/double storage precision
    │       ├── RingBuffer.hpp    # Ring buffer type
    │       └── TemporalField.hpp # Multi-level temporal storage
    │
    ├── data/                     # Field data structures
    │   ├── ModelFields.hpp       # Master data container
    │   ├── HydrodynamicState.hpp # Elevation/velocity state
    │   ├── PhysicsData.hpp       # Physics parameters
    │   └── WetDryData.hpp        # Wet/dry algorithm state
    │
    ├── geometry/                 # Mesh and geometry
    │   ├── Mesh.hpp              # Top-level mesh: HostMesh + DeviceMesh + dist context
    │   ├── HostMesh.hpp          # Host-side topology (nodes, elements, boundaries)
    │   ├── DeviceMesh.hpp        # Device-side FE cache (gradients, areas, CSR)
    │   ├── Element.hpp           # Triangle element
    │   ├── Node.hpp              # Mesh node
    │   ├── Point.hpp             # 2D/3D point
    │   ├── NeighborTable.hpp     # Mesh connectivity
    │   ├── BasisFunctions.hpp    # FE shape functions
    │   ├── MeshProjector.hpp # Global mesh operations
    │   ├── CoordinateRotation.hpp # Coordinate rotation
    │   ├── RotationData.hpp      # Rotation matrices
    │   ├── VelocityTransform.hpp # Velocity transformations
    │   ├── ProjectionScaleFactor.hpp  # Map projection scaling
    │   ├── ProjectionTransformer.hpp  # Coordinate transforms
    │   ├── NodalAttribute.hpp    # Nodal attribute type
    │   ├── NodalAttributeData.hpp     # Attribute storage
    │   ├── NodalAttributeRegistry.hpp # Attribute registry
    │   └── boundaries/           # Boundary data structures
    │       ├── BoundaryData.hpp
    │       ├── BoundaryRawData.hpp
    │       └── BoundaryType.hpp
    │
    ├── simulation/               # Simulation control
    │   ├── Simulation.hpp        # Main simulation driver
    │   ├── TimeStepper.hpp       # Time stepping logic
    │   ├── TimestepLogger.hpp    # Per-step status logging
    │   └── Diagnostics.hpp       # Solution monitoring
    │
    ├── numeric/                  # Numerical algorithms
    │   ├── terms/                 # Assembly adapters conforming to TermsConcepts.hpp
    │   │   ├── TermsConcepts.hpp  # GwcePhysicsTerm/GwceScalarTerm/MomentumElementTerm/MomentumNodalTerm
    │   │   ├── MomentumElementContext.hpp # Per-element input bundle for momentum terms
    │   │   ├── AtmPressureGradient.hpp    # Atmospheric pressure gradient (GWCE + momentum)
    │   │   ├── BottomFrictionTerm.hpp     # Bottom friction (GWCE)
    │   │   ├── Coriolis.hpp               # Coriolis with spherical correction (GWCE)
    │   │   ├── GwceAdvection.hpp          # Non-conservative advection (GWCE + momentum)
    │   │   ├── LateralStressTerm.hpp      # Lateral stress divergence (GWCE + momentum)
    │   │   ├── PressureGradientTerm.hpp   # Barotropic pressure gradient (momentum)
    │   │   ├── Tau0.hpp                   # Tau0 contribution (GWCE)
    │   │   ├── Tau0SpatialGradient.hpp    # Spatially varying tau0 gradient (GWCE)
    │   │   ├── TidePotential.hpp          # Tide potential gradient (GWCE)
    │   │   └── WindStressTerm.hpp         # Wind stress with depth-dependent limiter (GWCE)
    │   │
    │   ├── continuity/           # GWCE solver
    │   │   ├── GwceSolver.hpp            # Solver interface
    │   │   ├── GwceVectorAssembler.hpp   # RHS assembly
    │   │   ├── GwceVectorAssemblyKernels.hpp
    │   │   ├── GwcePreprocessingKernels.hpp
    │   │   ├── GwceCommonKernels.hpp     # Shared GWCE kernels
    │   │   ├── OpenBoundaryCoefficients.hpp
    │   │   ├── NodePositionSorter.hpp
    │   │   │
    │   │   ├── consistent/       # Implicit solver
    │   │   │   ├── GwceSolverConsistent.hpp
    │   │   │   ├── GwceMatrixAssemblerConsistent.hpp
    │   │   │   ├── GwceMatrixAssemblyConsistentKernels.hpp
    │   │   │   ├── ConjugateGradientSolver.hpp
    │   │   │   └── JacobiPreconditioner.hpp
    │   │   │
    │   │   └── lumped/           # Explicit solver
    │   │       ├── GwceSolverLumped.hpp
    │   │       ├── GwceMatrixAssemblerLumped.hpp
    │   │       └── GwceSolverLumpedKernels.hpp
    │   │
    │   ├── momentum/             # Momentum solver
    │   │   ├── MomentumSolver.hpp
    │   │   ├── MomentumSolveKernels.hpp
    │   │   └── MomentumRhsKernels.hpp
    │   │
    │   ├── boundary/             # Flux boundary conditions (QFORCE)
    │   │   ├── BoundaryConditionManager.hpp
    │   │   ├── BoundaryProcessor.hpp
    │   │   ├── BoundaryKernels.hpp
    │   │   ├── BoundaryState.hpp
    │   │   └── qforce/           # Per-type QFORCE structs + dispatch
    │   │
    │   └── wetdry/               # Wet/dry algorithm
    │       ├── WetDry.hpp
    │       └── WetDryKernels.hpp
    │
    ├── forcing/                  # Boundary and body forcing
    │   ├── ForcingManager.hpp    # Forcing orchestration
    │   ├── TanhRamp.hpp          # The one tanh ramp shape
    │   ├── Ramp.hpp              # Cold-start spinup ramp
    │   ├── WindowedRamp.hpp      # Activation/deactivation window fade
    │   ├── flow/                 # Flow boundary time series
    │   │   └── FlowTimeSeries.hpp
    │   ├── meteorological/       # Meteorological forcing
    │   │   ├── MeteoForcingProvider.hpp  # Ring buffer + interpolation orchestration (PIMPL)
    │   │   ├── MeteoForcingConfig.hpp    # Per-source + vortex configuration
    │   │   ├── MeteoSnapshotSource.hpp   # Rank-0 reader thread; union time axis + per-domain windows
    │   │   ├── MetSnapshot.hpp           # Device-side snapshot slot (base + moving-grid fields)
    │   │   ├── MeteoInterpolateKernel.hpp # Eulerian temporal blend + reduce-mask kernel
    │   │   ├── DeviceUpload.hpp          # Host->device vector/mask upload helpers
    │   │   ├── ForcingWindowTracker.hpp  # Per-source start/stop log transitions
    │   │   ├── VortexOverlay.hpp         # GAHM device table + per-step vortex blend
    │   │   ├── WindReduction.hpp         # Directional roughness + canopy correction
    │   │   ├── WindStress.hpp            # wind_stress<DragLaw>() kinematic stress
    │   │   └── DragLaw.hpp               # Wind drag formulations
    │   └── tide/                 # Tidal forcing
    │       ├── boundary/         # Tidal boundary conditions
    │       │   ├── TideBoundaryForcing.hpp
    │       │   └── TideBoundaryConstituent.hpp
    │       └── potential/        # Tide potential forcing
    │           ├── TidePotentialInterface.hpp
    │           ├── astronomical/     # Astronomical tide potential
    │           │   ├── AstronomicalTidePotential.hpp
    │           │   ├── AstronomicalTidePotentialAdapter.hpp
    │           │   ├── AstronomicParameters.hpp
    │           │   ├── AstronomicConstants.hpp
    │           │   ├── MoonPosition.hpp
    │           │   ├── SunPosition.hpp
    │           │   ├── MoonSunPositionCalculator.hpp
    │           │   └── SiderealTime.hpp
    │           └── harmonics/        # Harmonic tide potential
    │               ├── TidePotentialHarmonics.hpp
    │               └── TidePotentialConstituent.hpp
    │
    ├── physics/                  # Physics kernels
    │   ├── BottomFriction.hpp    # Manning's friction
    │   ├── LateralStress.hpp     # Lateral stress tensor
    │   ├── InlineLateralStress.hpp # Inline stress computation
    │   ├── PressureGradient.hpp  # Pressure gradient terms
    │   ├── NonConservativeAdvection.hpp # Advection terms
    │   └── Smagorinsky.hpp       # Smagorinsky turbulence
    │
    ├── distributed/              # MPI distributed computing
    │   ├── DistributedContext.hpp    # MPI context management
    │   ├── CommunicationManager.hpp  # Communication orchestration
    │   ├── GhostExchange.hpp         # Ghost node exchange
    │   ├── MapFactory.hpp            # Tpetra map creation
    │   └── partition/            # Mesh partitioning
    │       ├── MeshPartitioner.hpp   # Zoltan2 partitioning
    │       ├── PartitionCache.hpp    # Cached partitions
    │       └── PartitionInfo.hpp     # Partition metadata
    │
    └── io/                       # Kernel I/O (split by direction)
        ├── async/                # Background writer thread
        │   ├── OutputWriterThread.hpp  # Bounded writer thread pool
        │   └── OutputSnapshot.hpp      # Recycled write buffers
        ├── checkpoint/           # Checkpoint write / restart
        │   ├── CheckpointReader.hpp
        │   └── CheckpointWriter.hpp
        ├── input/                # Read-side
        │   ├── MeshReader.hpp
        │   ├── MeshGather.hpp
        │   ├── NodalAttributeReader.hpp
        │   ├── NodalAttributeInitializer.hpp
        │   └── SalDataReader.hpp
        └── output/               # Write-side
            ├── OutputManager.hpp
            ├── OutputFile.hpp
            ├── OutputFieldAccess.hpp      # Field view lookup for output
            ├── OutputVariable.hpp
            ├── OutputVariableRegistry.hpp
            ├── NetcdfFileSetup.hpp        # Shared UGRID file setup
            └── DistributedIO.hpp          # Parallel field gather

Namespace Organization

Fig. 37 Namespace hierarchy within the Cocoa project

Meteorological Forcing Pipeline

The meteorological forcing subsystem moves atmospheric data from one or more configured sources through a staged pipeline into device-resident Kokkos views consumed by the GWCE and momentum kernels. Each configured source (forcing.meteorological.domains entry, or the flat single-source form) gets its own reader, of any format; cocoa_meteo composes them into a single flat domain space, handles file reading, and performs the host-side spatial combine. The cocoa_kernel forcing subsystem manages the union time axis, temporal buffering, per-source wind reduction, activation fades, and device transfer.

Fig. 38 Meteorological forcing data pipeline

Pipeline Stages

1. Reader Construction (rank 0, at startup). MeteoReaderFactory::create_reader builds one reader per configured source (CfNetcdfReader, OwiAsciiReader, OwiNetcdfReader, or GribReader when built with ecCodes). A single source is used directly; two or more are wrapped by CompositeMeteoReader, which flattens child order then child-local domain order into one flat index space – this flat order IS the coverage priority order (list order = priority; later wins at overlap).

2. Dataset Metadata + Union Time Axis (rank 0). MeteoSnapshotSource::dataset_info() queries every domain’s own time axis and grid, builds the sorted-unique union axis, and computes each domain’s window on that axis (the union-index span between the domain’s own first and last sample). A background std::jthread then reads snapshots in increasing union-index order.

3. Per-Snapshot Read + Cadence Interpolation (rank 0, background thread). For each union-axis snapshot, a domain either serves an exact sample from its own time axis or, when the union time falls strictly between two of its own samples (a coarser-cadence source), linearly interpolates its fields – and, for a moving grid, its coordinates – between its bracketing snapshots. Regular-moving (storm) domains are read and cached the same way but excluded from the host spatial combine in step 5; they get a per-timestep Lagrangian overlay on device instead (step 8).

4. Broadcast. Dataset metadata is sent once; each snapshot’s grid-space fields are sent to every rank via MPI_Bcast, so only rank 0 touches the filesystem.

5. Host Spatial Combine. MultiDomainInterpolator::interpolate bilinearly interpolates every active, non-moving domain onto the local mesh nodes and combines them: the innermost (highest-index) domain covering a node wins, subject to that domain’s per-node activation-window weight (a WindowedRamp evaluated once per snapshot), which cross-fades the winning domain toward the underlying field instead of hard-overriding it. The combine also records, per node, which domain supplied the value (the winner mask), used to build the per-source wind-reduction mask.

6. Host-to-Device Transfer. The combined pressure/wind_u/wind_v node arrays, the per-node reduce mask, and each moving domain’s raw grid-space fields are uploaded into a device-resident MetSnapshot slot via Kokkos::deep_copy.

7. Temporal Interpolation + Reduction (device). At each simulation timestep, interpolate_met_fields blends the two bracketing ring-buffer slots using the interpolation weight \(\alpha\) (see Temporal Interpolation) and, where any source reduces, ramps the per-node wind-reduction mask across the same bracket so a node whose winning source changes mid-bracket has no discontinuity.

8. Lagrangian Moving-Domain Overlay (device). Each regular-moving (storm-following) domain is sampled directly from its own advected grid at the current simulation time and overlaid onto the base field, in domain order.

9. Vortex Overlay (device). VortexOverlay::apply blends each active GAHM storm’s analytic wind/pressure into the current fields, applying its own activation WindowedRamp and, optionally, the same land wind-reduction correction as a gridded source.

Ring Buffer

The MeteoForcingProvider maintains a fixed-size circular buffer of MetSnapshot slots on the device (default 8 slots on CPU backends, 32 on GPU backends; configurable via cocoa_METEO_BUFFER_SLOTS_CPU / cocoa_METEO_BUFFER_SLOTS_GPU at configure time). Each slot holds the spatially-combined base fields (pressure, wind_u, wind_v, reduce mask) plus the raw grid-space fields of every moving (storm) domain. The slots are allocated once at initialization and reused throughout the simulation.

At any given simulation time, two adjacent slots form the active bracket: the snapshot immediately before and after the current time. The temporal interpolation weight \(\alpha\) selects the blend between them. The remaining slots are pre-filled with upcoming snapshots so that advancing the bracket never requires a synchronous file read followed by a device transfer.

When the simulation time passes the end of the current bracket:

1. The oldest slot is recycled (the ring buffer pops its front).

2. The next unread snapshot is read, spatially combined, and uploaded into the newly available slot.

3. The bracket indices advance by one.

This design provides two benefits:

• Amortized transfer cost. Because upcoming snapshots are pre-loaded into free slots, the GPU is never stalled waiting for a single deep_copy at bracket boundaries.

• No redundant work. Once a snapshot is spatially combined and uploaded, it remains device-resident until the buffer rotates past it.

Meteorological Input Distribution

Rank 0 reads the meteorological files (through however many sources are configured) and broadcasts each snapshot’s dataset metadata and grid-space fields to all compute ranks via MPI_Bcast, pre-reading upcoming snapshots on a background thread so file access overlaps computation. Every rank performs its own host spatial combine and maintains its own ring buffer of device-side data for its local mesh partition.