==================== Extending the Model ==================== This guide covers three extension points: nodal attributes (spatially varying parameters), physics contributions (terms in the GWCE and momentum equations), and meteorological file formats. .. contents:: On This Page :local: :depth: 2 .. _adding-a-nodal-attribute: Adding a Nodal Attribute ------------------------ Nodal attributes are per-node fields loaded from the mesh file or set to a constant value from configuration. Examples include Manning's n and the primitive weighting parameter tau0. All nodal attributes are accessed via the ``NodalAttributeId`` enum. Strings appear only at the I/O boundary (reading from NetCDF); the internal API is entirely enum-based for O(1) access with no string hashing at runtime. Step 1: Register the Enum ^^^^^^^^^^^^^^^^^^^^^^^^^ Add an entry to ``NodalAttributeId`` in ``geometry/NodalAttribute.hpp``: .. code-block:: cpp enum class NodalAttributeId : int { ManningsN = 0, Tau0 = 1, // ... the other existing attributes ... AdvectionState = 6, AppleSauce = 7, // New attribute }; /// Number of defined nodal attributes (update when adding new entries) static constexpr size_t NUM_NODAL_ATTRIBUTES = 8; static_assert(static_cast(NodalAttributeId::AppleSauce) + 1 == NUM_NODAL_ATTRIBUTES, "Update NUM_NODAL_ATTRIBUTES when adding new attributes"); Add the string mapping in ``attribute_name()`` in the same file: .. code-block:: cpp constexpr auto attribute_name(const NodalAttributeId id) -> const char* { switch (id) { case NodalAttributeId::ManningsN: return "mannings_n_at_sea_floor"; case NodalAttributeId::Tau0: return "primitive_weighting_in_continuity_equation"; // ... the other existing attributes ... case NodalAttributeId::AppleSauce: return "apple_sauce_coefficient"; } Kokkos::abort("attribute_name: unhandled NodalAttributeId"); return ""; } The switch has no ``default``, so ``-Wswitch`` flags a missing case at compile time. The string returned by ``attribute_name()`` must match the NetCDF variable name in the mesh file. Step 2: Initialize the Attribute ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Add an initialization method to ``io/input/NodalAttributeInitializer.cpp``. The existing ``initialize_scalar_attribute()`` helper handles the common pattern of reading from mesh or using a constant: .. code-block:: cpp void NodalAttributeInitializer::initialize_apple_sauce() { constexpr double default_value = 1.0; constexpr double valid_min = 0.0; constexpr double valid_max = 100.0; initialize_scalar_attribute( Geometry::NodalAttributeId::AppleSauce, m_config.physics.apple_sauce_coeff, default_value, "m", "apple sauce coefficient", valid_min, valid_max); } Call it from ``initialize()``: .. code-block:: cpp void NodalAttributeInitializer::initialize() { m_mesh->initialize_attributes(m_fields.num_nodes()); initialize_manning_n(); initialize_tau0(); // ... the other attributes ... initialize_advection_state(); initialize_apple_sauce(); // New initialize_initial_water_level(); initialize_wetdry_state(); } Step 3: Access in Kernels ^^^^^^^^^^^^^^^^^^^^^^^^^ Access the attribute in any kernel via the enum: .. code-block:: cpp const auto apple_sauce = mesh->attributes().get(Geometry::NodalAttributeId::AppleSauce); Kokkos::parallel_for(..., KOKKOS_LAMBDA(const int i) { const auto val = apple_sauce(i); // Works for constant or per-node }); ``get()`` returns a ``const NodalAttributeAccessor&`` from a pre-built cache array: a single array index, no lookup. If the attribute was not registered, ``get()`` throws. Step 4: Add Configuration ^^^^^^^^^^^^^^^^^^^^^^^^^ Add the configuration parameter to the appropriate config struct in ``core/config/PhysicsConfig.hpp`` and parse it in ``core/config/ModelConfigurationFactory.cpp``. A value that is either a constant or read from the mesh is a ``std::optional``: ``parse_mesh_or_numeric()`` turns the string ``"mesh"`` into ``std::nullopt`` (read per-node values from the mesh file) and a number into the constant. Add the key to ``core/config/cocoa_config.schema.yaml`` as well (``mesh_or_numeric``); the schema rejects keys it does not list. Step 5: Update Tests ^^^^^^^^^^^^^^^^^^^^ Register the new attribute in test fixtures that exercise the kernels using it: .. code-block:: cpp mesh->attributes().register_constant( NodalAttributeId::AppleSauce, 1.0, "m"); Vector (Multi-Component) Attributes ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Some attributes have multiple values per node (e.g., internal tide friction stores 3 tensor components). These use vector storage. Read 2D variables ``[num_nodes, num_components]`` with ``NodalAttributeReader::read_vector_attribute()`` or ``read_vector_attribute_distributed()``: .. code-block:: cpp auto attr = m_mesh->is_parallel() ? NodalAttributeReader::read_vector_attribute_distributed( reader, variable_name, m_mesh->partition(), 0.0, "1/s") : NodalAttributeReader::read_vector_attribute( reader, variable_name, num_nodes, 0.0, "1/s"); The two-argument accessor returns one component: .. code-block:: cpp const auto accessor = mesh->attributes().get(NodalAttributeId::MyVectorAttr); Kokkos::parallel_for(..., KOKKOS_LAMBDA(const int i) { const auto xx = accessor(i, 0); // First component const auto yy = accessor(i, 1); // Second component const auto xy = accessor(i, 2); // Third component }); If an attribute may be either scalar (1D in the mesh file) or vector (2D), detect the shape with ``reader.get_variable_dimensions()`` and register accordingly, as ``initialize_internal_tide_friction()`` does. Optional Attributes ^^^^^^^^^^^^^^^^^^^ An optional attribute applies only when the user enables it. Skip its registration when it is disabled, and check ``has()`` before accessing it: .. code-block:: cpp // In initialization: skip if disabled if (config_value == 0.0) { return; // Don't register -> has() returns false } // In kernel: check before applying if (!mesh->attributes().has(NodalAttributeId::MyOptionalAttr)) { return; } const auto accessor = mesh->attributes().get(NodalAttributeId::MyOptionalAttr); // ... apply See ``internal_tide_friction`` for a complete example of an optional vector attribute. Adding a Physics Contribution ----------------------------- Physics contributions add terms to the GWCE (continuity) and momentum equation solvers. The framework uses C++20 concepts to enforce compile-time contracts while keeping physics math out of the assembly kernel lambdas. Architecture ^^^^^^^^^^^^ GWCE contributions use nested ``Input`` structs constructed inline at the call site. Each contribution's ``Input`` holds only the data that contribution needs, and that data dies when ``compute()`` returns, which keeps peak register pressure on GPU down. Momentum contributions share one ``MomentumElementContextT`` struct (``MomentumElementContext`` is the triangle alias). Terms are written once for both element types: a term whose only dependence on the element type is its vertex count is templated on the element-type tag (``Geometry::Triangle``, ``Geometry::Quadrilateral``), as ``GwceTau0T`` is, and each element kernel instantiates it. .. graphviz:: :align: center :caption: Physics contribution architecture for GWCE and Momentum digraph physics_arch { rankdir=TB; node [shape=box, style="filled,rounded", fontname="Helvetica", fontsize=10]; edge [color="#555555"]; bgcolor="transparent"; kernel [label="Assembly Kernel Lambda\n(loads data, constructs Input,\ncalls contributions, assembles RHS)", fillcolor="#E0E0E0", color="#616161"]; subgraph cluster_gwce { label="GWCE Path"; style="filled,rounded"; fillcolor="#E3F2FD"; color="#1E88E5"; fontname="Helvetica"; ginput [label="Contribution::Input\n(per-contribution struct)", fillcolor="#BBDEFB", color="#1565C0"]; gcompute [label="static compute(Input)\n-> Types::Vec2 {jx, jy}", fillcolor="#BBDEFB", color="#1565C0"]; } subgraph cluster_mom { label="Momentum Path"; style="filled,rounded"; fillcolor="#FFF3E0"; color="#FB8C00"; fontname="Helvetica"; mctx [label="MomentumElementContextT\n(shared per-node data)", fillcolor="#FFE0B2", color="#E65100"]; mcompute [label="static compute(ctx)\n-> Types::Vec2 {x, y}\nor -> Array, N>", fillcolor="#FFE0B2", color="#E65100"]; } kernel -> ginput; kernel -> mctx; ginput -> gcompute [label="data dies\nafter call", fontsize=8, color="#43A047", fontcolor="#43A047"]; mctx -> mcompute; } Concepts """""""" Three C++20 concepts in ``numeric/cg/terms/TermsConcepts.hpp`` define the contracts. The momentum concepts are checked per element type ``Elem``: .. list-table:: :header-rows: 1 :widths: 30 40 30 * - Concept - Required Signature - Use Case * - ``GwcePhysicsTerm`` - ``static compute(const T::Input&)`` ``-> Types::Vec2`` - Returns (jx, jy) contribution to both x- and y-equation GWCE RHS * - ``MomentumElementTerm`` - ``static compute(const MomentumElementContextT&)`` ``-> Types::Vec2`` - Returns one (x, y) contribution that reaches every vertex of the element * - ``MomentumNodalTerm`` - ``static compute(const MomentumElementContextT&)`` ``-> Kokkos::Array, Elem::num_nodes>`` - Returns per-vertex (x, y) contributions to both x- and y-momentum Step-by-Step Guide ^^^^^^^^^^^^^^^^^^ .. graphviz:: :align: center :caption: Workflow for adding a new physics contribution digraph workflow { rankdir=TB; node [shape=box, style="filled,rounded", fontname="Helvetica", fontsize=10]; edge [color="#555555"]; bgcolor="transparent"; s1 [label="Step 1: Data Storage\nAdd fields to PhysicsData", fillcolor="#FFEBEE", color="#E53935"]; s2 [label="Step 2: GWCE Contribution\nCreate numeric/cg/terms/ header\nwith Input struct + compute()", fillcolor="#E3F2FD", color="#1E88E5"]; s3 [label="Step 3: Momentum Contribution\nCreate numeric/cg/terms/ header\nwith compute()", fillcolor="#E3F2FD", color="#1E88E5"]; s4 [label="Step 4: Wire GWCE Kernel\nAdd view to RhsNodalFields +\ncall in assemble_element_rhs", fillcolor="#FFF3E0", color="#FB8C00"]; s5 [label="Step 5: Wire Momentum Kernel\nAdd view to MomentumRhsFields +\nload context in MomentumRhsKernels", fillcolor="#FFF3E0", color="#FB8C00"]; s6 [label="Step 6: Populate Data\nNodal attribute, computed kernel,\nor external forcing", fillcolor="#E8F5E9", color="#43A047"]; s7 [label="Step 7: Ghost Sum\nsum_and_distribute_fields\n(if element-assembled)", fillcolor="#F3E5F5", color="#8E24AA"]; s8 [label="Step 8: Configuration + Build\nAdd config, parse YAML,\nupdate CMakeLists.txt", fillcolor="#E0E0E0", color="#616161"]; s1 -> s2 -> s3 -> s4 -> s5 -> s6 -> s7 -> s8; s2 -> s3 [style=invis, weight=10]; } The following steps use a hypothetical physics term called "AppleSauce" as an example. Step 1: Data Storage """""""""""""""""""" If the new term needs per-node fields (loaded from input or computed each timestep), add a sub-container to ``data/PhysicsData.hpp``: .. code-block:: cpp /// AppleSauce forcing sub-container struct AppleSauceFields { DoubleView coefficient; DoubleView intensity; }; Add a private member, allocate in the constructor, and expose with const/non-const accessors: .. code-block:: cpp // In the public section: [[nodiscard]] auto apple_sauce() -> AppleSauceFields& { return m_apple_sauce; } [[nodiscard]] auto apple_sauce() const -> const AppleSauceFields& { return m_apple_sauce; } // In the private section: AppleSauceFields m_apple_sauce; Follow the pattern of ``FrictionFields``, ``LateralStressFields``, and ``TidePotentialFields``. Step 2: GWCE Contribution """"""""""""""""""""""""" Create ``numeric/cg/terms/AppleSauce.hpp``. Each GWCE contribution nests an ``Input`` struct containing only the data it needs, and provides a static ``compute(const Input&)`` method: .. code-block:: cpp #pragma once #include #include "cocoa_kernel/geometry/ElementGather.hpp" #include "cocoa_kernel/geometry/ElementType.hpp" #include "cocoa_kernel/numeric/cg/terms/TermsConcepts.hpp" namespace Cocoa::Numeric::GwceRhsTerms { template struct GwceAppleSauceT { using Pack = Kokkos::Array; struct Input { Pack apple_sauce; // Add only the fields this contribution needs }; KOKKOS_INLINE_FUNCTION static auto compute(const Input& input) -> Types::Vec2 { const auto as_avg = Geometry::avg(input.apple_sauce); // Physics math goes here const auto jx = /* ... */; const auto jy = /* ... */; return {jx, jy}; } }; static_assert(GwcePhysicsTerm>); static_assert(GwcePhysicsTerm>); } // namespace Cocoa::Numeric::GwceRhsTerms The ``static_assert`` lines check the concept contract for each element type at compile time, so a mismatched signature fails there with a readable error rather than inside a template instantiation at the call site. A term whose quadrilateral form is a different expression from the triangle one (rather than the same expression over four vertices) gets its own struct beside the triangle one (see ``GwceTau0SpatialGradient`` and ``GwceQuadTau0SpatialGradient``). Step 3: Momentum Contribution """"""""""""""""""""""""""""" Add the momentum term to ``numeric/cg/terms/AppleSauce.hpp`` (alongside the GWCE term), in ``namespace Cocoa::Numeric::MomentumTerms``. Choose the concept based on whether the result is the same at every vertex or varies per vertex. **Uniform across vertices**: follow ``MomentumAdvection`` in ``GwceAdvection.hpp``: .. code-block:: cpp struct MomentumAppleSauce { template KOKKOS_INLINE_FUNCTION static auto compute( const MomentumElementContextT& ctx) -> Types::Vec2 { // Physics math: same value scattered to every vertex return {/* x */, /* y */}; } }; static_assert(MomentumElementTerm); static_assert( MomentumElementTerm); **Different per vertex**: follow ``MomentumLateralStress`` in ``LateralStressTerm.hpp``: .. code-block:: cpp struct MomentumAppleSauce { template KOKKOS_INLINE_FUNCTION static auto compute( const MomentumElementContextT& ctx) -> Kokkos::Array, Elem::num_nodes> { // Physics math: different value per vertex return {/* one Vec2 per vertex */}; } }; static_assert(MomentumNodalTerm); static_assert( MomentumNodalTerm); Step 4: Wire Into GWCE Kernel """"""""""""""""""""""""""""" In ``numeric/cg/continuity/GwceVectorAssemblyKernels.hpp``, the views the RHS assembly reads are gathered on the host into ``RhsNodalFields`` (built by ``rhs_nodal_fields()``) and captured by value; the per-element body is ``assemble_element_rhs``, which ``rhs_assembly_kernel`` runs over each type's wet list. Make three changes: **a) Include the header and add the view** to ``RhsNodalFields`` and its builder, alongside the other physics views: .. code-block:: cpp // At file top: #include "cocoa_kernel/numeric/cg/terms/AppleSauce.hpp" // In RhsNodalFields: Types::Scalar1DView apple_sauce; // In rhs_nodal_fields(), alongside the other physics views: .apple_sauce = physics.apple_sauce().coefficient, **b) Call the contribution** in ``assemble_element_rhs``, constructing the ``Input`` inline with aggregate initialization and accumulating into the local result: .. code-block:: cpp using GwceRhsTerms::GwceAppleSauceT; jxy += GwceAppleSauceT::compute( {.apple_sauce = Geometry::gather_element_nodes(fields.apple_sauce, nodes)}); The result accumulates into the ``Types::Vec2 jxy`` variable, and the ``Input`` temporary (with any data gathered into it) dies after ``compute()`` returns. Step 5: Wire Into Momentum Kernel """""""""""""""""""""""""""""""""" In ``numeric/cg/momentum/MomentumRhsKernels.cpp``, the same pattern: the views are gathered into ``MomentumRhsFields`` (built by ``momentum_rhs_fields()``), and the per-element body is ``accumulate_element_momentum_rhs``, which ``momentum_rhs_element_kernel`` launches. **a) Include and add the view** to ``MomentumRhsFields`` and its builder: .. code-block:: cpp #include "cocoa_kernel/numeric/cg/terms/AppleSauce.hpp" // In MomentumRhsFields: Types::Scalar1DView apple_sauce; **b) Load into context.** Add a ``Pack apple_sauce{};`` member to ``MomentumElementContextT`` (``numeric/cg/terms/MomentumElementContext.hpp``) and gather it with the others: .. code-block:: cpp ctx.apple_sauce = Geometry::gather_element_nodes(fields.apple_sauce, ctx.nm); **c) Call the contribution.** For ``MomentumElementTerm``, add to the base contribution: .. code-block:: cpp base_contrib += MomentumAppleSauce::compute(ctx); For ``MomentumNodalTerm``, add in the per-vertex scatter loops alongside lateral stress (each element type has its own loop): .. code-block:: cpp const auto apple = MomentumAppleSauce::compute(ctx); for (size_t i = 0; i < Elem::num_nodes; ++i) { const auto contrib = base_contrib + lateral[i] + baroclinic_contrib[i] + apple[i]; Kokkos::atomic_add(&rhs_x(ctx.nm[i]), contrib.x); Kokkos::atomic_add(&rhs_y(ctx.nm[i]), contrib.y); } Step 6: Populate the Data """"""""""""""""""""""""" How the data gets into ``AppleSauceFields`` depends on its nature: **Static nodal attribute** (loaded once from mesh file): Register a ``NodalAttributeId`` and use ``NodalAttributeInitializer`` as described in :ref:`adding-a-nodal-attribute` above. **Computed each timestep** (like lateral stress): Write a computation kernel that populates ``AppleSauceFields``. Call it from ``run_preprocessing`` (``numeric/cg/continuity/GwcePreprocessingKernels.hpp``), which ``GwceVectorAssembler::assemble`` runs before the RHS assembly: .. code-block:: cpp // In run_preprocessing: COCOA_PROFILE_REGION("Preprocessing"); initialize_rhs_fields(model_fields, rhs_1d, model_fields.num_nodes()); elemental_preprocessing_kernel(model_fields, config); elemental_preprocessing_kernel(model_fields, config); sum_lateral_stress(model_fields, comm_manager); apple_sauce_kernel(model_fields, config); // New kernel nodal_postprocessing_kernel(model_fields); **External forcing** (like wind or atmospheric pressure): Add a forcing source to ``ForcingManager`` that updates the fields each timestep using the appropriate ramp factor (``m_tidal_ramp`` for tidal forcing, ``m_meteorological_ramp`` for meteorological forcing). Step 7: Ghost Sum (if needed) """"""""""""""""""""""""""""" A ghost sum is only needed if the field is **assembled from element contributions via atomics**. Elements have no ghosts, so each rank holds only its own elements' share of a node on a partition cut, owned or ghost. Before any rank reads the field, sum every rank's share into the owner and copy the total back to the ghosts with ``CommunicationManager::sum_and_distribute_fields``. Copying the owner's value to the ghosts alone (``exchange_scalar_fields``) is not enough: the owner's value is itself only a partial sum. Fields that are directly loaded from nodal attributes or computed per-node (without element scatter) need neither. Follow ``sum_lateral_stress`` (``GwcePreprocessingKernels.cpp``), which sums the lateral stress components before ``nodal_postprocessing_kernel`` normalizes them: .. code-block:: cpp void sum_apple_sauce(Data::ModelFields& model_fields, Distributed::CommunicationManager& comm_manager) { std::array fields = { Types::StridedScalar1DView(model_fields.physics().apple_sauce().intensity) }; comm_manager.sum_and_distribute_fields(fields); // no-op in serial } Call it after the AppleSauce computation kernel and before anything that normalizes or reads the field. .. note:: Pass several fields in one ``std::array`` to combine them in a single round-trip pair rather than calling ``sum_and_distribute_field`` once per view. ``max_and_distribute_fields`` is the exact reduction for flags and bounds any rank may raise. Step 8: Configuration and Build """"""""""""""""""""""""""""""" **Configuration.** Add parameters to ``core/config/PhysicsConfig.hpp``: .. code-block:: cpp struct AppleSauceConfig { bool active{false}; double coefficient{1.0}; }; Add ``AppleSauceConfig apple_sauce;`` to ``PhysicsConfig``. Parse in ``core/config/ModelConfigurationFactory.cpp``, and add the keys to ``core/config/cocoa_config.schema.yaml``. **CMakeLists.txt.** Add new headers to ``COCOA_KERNEL_HEADERS`` in ``src/cocoa_kernel/CMakeLists.txt``: .. code-block:: text numeric/cg/terms/AppleSauce.hpp Reference: Existing Contributions ---------------------------------- .. list-table:: :header-rows: 1 :widths: 35 20 45 * - File / Struct - Concept - Description * - ``Coriolis.hpp`` / ``GwceCoriolis`` - ``GwcePhysicsTerm`` - Coriolis force (jx, jy) with spherical correction * - ``BottomFrictionTerm.hpp`` / ``GwceBottomFrictionT`` - ``GwcePhysicsTerm`` - TKM bottom friction tensor (jx, jy) * - ``Tau0.hpp`` / ``GwceTau0T`` - ``GwcePhysicsTerm`` - Spatially variable tau0 weighting (jx, jy) * - ``TidePotential.hpp`` / ``GwceTidePotentialT`` - ``GwcePhysicsTerm`` - Tide potential gradient (jx, jy) * - ``WindStressTerm.hpp`` / ``GwceWindStressT`` - ``GwcePhysicsTerm`` - Kinematic wind stress (jx, jy) via ``wind_stress()`` * - ``AtmPressureGradient.hpp`` / ``GwceAtmPressureGradientT`` - ``GwcePhysicsTerm`` - Atmospheric pressure gradient (jx, jy) * - ``BaroclinicGradientTerm.hpp`` / ``GwceBaroclinicGradientT`` - ``GwcePhysicsTerm`` - Baroclinic pressure gradient (jx, jy) * - ``Tau0SpatialGradient.hpp`` / ``GwceTau0SpatialGradient``, ``GwceQuadTau0SpatialGradient`` - none (``compute`` returns ``double``) - Spatially varying tau0 gradient, a scalar (triangle and quadrilateral forms) * - ``LateralStressTerm.hpp`` / ``GwceLateralStressT`` - ``GwcePhysicsTerm`` - Lateral stress divergence (jx, jy) * - ``GwceAdvection.hpp`` / ``GwceAdvectionT`` - ``GwcePhysicsTerm`` - Non-conservative advection (jx, jy) * - ``PressureGradientTerm.hpp`` / ``MomentumPressureGradient`` - ``MomentumElementTerm`` - Pressure gradient (x, y) (uniform per element) * - ``GwceAdvection.hpp`` / ``MomentumAdvection`` - ``MomentumElementTerm`` - Non-conservative advection (x, y) (uniform) * - ``AtmPressureGradient.hpp`` / ``MomentumAtmPressureGradient`` - ``MomentumElementTerm`` - Atmospheric pressure gradient (x, y) (uniform) * - ``BaroclinicGradientTerm.hpp`` / ``MomentumBaroclinicGradient`` - ``MomentumNodalTerm`` - Baroclinic pressure gradient (x, y) (per-node) * - ``LateralStressTerm.hpp`` / ``MomentumLateralStress`` - ``MomentumNodalTerm`` - Inline lateral stress (x, y) (per-node) .. _adding-meteo-format: Adding a Meteorological File Format ------------------------------------ The ``cocoa_meteo`` library provides a pluggable reader framework for meteorological forcing data. To add support for a new file format, implement a reader class that derives from ``MeteoReaderBase`` and register it in the factory. Multi-source composition (mixing formats, or several instances of the same format) is handled by ``CompositeMeteoReader``, the wrapper ``MeteoReaderFactory`` builds over two or more configured sources. A new format reader therefore composes with every existing one without any of the steps below touching the composite. Architecture Overview ^^^^^^^^^^^^^^^^^^^^^ .. graphviz:: :align: center :caption: cocoa_meteo reader architecture digraph meteo_arch { rankdir=TB; node [shape=box, style="filled,rounded", fontname="Helvetica", fontsize=10]; edge [color="#555555"]; bgcolor="transparent"; factory [label="MeteoReaderFactory\ncreate_reader(config)", fillcolor="#E0E0E0", color="#616161"]; base [label="MeteoReaderBase\n(abstract interface)", fillcolor="#E3F2FD", color="#1E88E5"]; cf [label="CfNetcdfReader", fillcolor="#BBDEFB", color="#1565C0"]; owi_a [label="OwiAsciiReader", fillcolor="#BBDEFB", color="#1565C0"]; owi_n [label="OwiNetcdfReader", fillcolor="#BBDEFB", color="#1565C0"]; grib [label="GribReader\n(cocoa_ENABLE_GRIB)", fillcolor="#BBDEFB", color="#1565C0"]; yours [label="YourFormatReader", fillcolor="#C8E6C9", color="#2E7D32"]; composite [label="CompositeMeteoReader\n(composes 2+ sources of\nany format, list order = priority)", fillcolor="#F5F5F5", color="#9E9E9E"]; factory -> base [label="dispatches on\nMeteoFormat, per source"]; base -> cf; base -> owi_a; base -> owi_n; base -> grib; base -> yours [style=dashed, label="new"]; factory -> composite [label="2+ sources", style=dashed]; composite -> base [style=dashed, label="wraps"]; } Step 1: Add the Format Enum ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Add an entry to ``MeteoFormat`` in ``cocoa_meteo/core/MeteoFormat.hpp``: .. code-block:: cpp enum class MeteoFormat { CF_NETCDF, OWI_ASCII, OWI_NETCDF, GRIB, YOUR_FORMAT }; Step 2: Create the Reader Class ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Create a new directory under ``cocoa_meteo/readers/your_format/`` and implement a reader that inherits from ``MeteoReaderBase``: .. code-block:: cpp #pragma once #include "cocoa_meteo/core/MeteoReaderBase.hpp" #include "cocoa_meteo/core/MeteoReaderConfig.hpp" namespace Cocoa::Meteo { class YourFormatReader : public MeteoReaderBase { public: // Construction opens the file(s), indexes snapshot positions, and fails // loud on any format error. There is no separate open() step, so a // constructed reader is always usable. explicit YourFormatReader(const MeteoReaderConfig& config); [[nodiscard]] auto num_domains() const -> size_t override; [[nodiscard]] auto num_snapshots() const -> size_t override; [[nodiscard]] auto grid_info(size_t domain_idx) const -> MeteoGrid override; [[nodiscard]] auto time_values() const -> const std::vector& override; [[nodiscard]] auto read_snapshot(size_t snapshot_idx, size_t domain_idx) -> MeteoFieldSnapshot override; private: // Format-specific internal state }; } // namespace Cocoa::Meteo The pure virtual methods define the contract: .. list-table:: :header-rows: 1 :widths: 30 70 * - Method - Description * - ``num_domains()`` - Number of grid domains (1 for single-domain formats) * - ``num_snapshots()`` - Total number of time snapshots in the file * - ``grid_info(domain_idx)`` - Return a ``MeteoGrid`` (``RegularGrid`` or ``IrregularGrid``) * - ``time_values()`` - Return snapshot times as seconds since Unix epoch * - ``read_snapshot(idx, domain)`` - Read and return a ``MeteoFieldSnapshot`` for the given time/domain The ``MeteoFieldSnapshot`` bundles the grid definition with the three forcing fields (pressure, wind_u, wind_v), and the sea-ice concentration where the source carries it, for one time step and one domain. ``read_snapshot`` indexes the domain's own time axis (``domain_time_values(domain)``) and must throw on an out-of-range index. ``MeteoReaderBase`` also declares six **optional** virtuals with defaults that are correct for most formats; override only what your format needs: .. list-table:: :header-rows: 1 :widths: 30 20 50 * - Method - Default - Override when... * - ``is_moving_domain(domain_idx)`` - ``false`` - the domain's grid coordinates change between snapshots (e.g. a storm-following nest); moving domains are re-sent with every snapshot broadcast and have their interpolation weights recomputed. * - ``domain_name(domain_idx)`` - ``""`` (empty) - the format has a natural identity for the domain (a filename, file pair, or group name) to show in start/stop log messages. * - ``domain_time_values(domain_idx)`` - the shared ``time_values()`` - domains within one reader instance do not share a single time axis (e.g. a nested grid active only during part of the run). * - ``domain_source(domain_idx)`` - ``0`` - never override this in a leaf reader. It exists so ``CompositeMeteoReader`` can report which configured source (child reader) owns a flattened domain, routing the per-source ``wind_reduction`` flag to the right nodes; a single-format reader always has exactly one source, so the default is always correct. * - ``domain_has_ice(domain_idx)`` - ``false`` - the format can carry sea-ice concentration; the domain must then carry it at every one of its snapshots. * - ``domain_has_met_fields(domain_idx)`` - ``true`` - the format can legally omit pressure and wind (an ice-only source, as ``CfNetcdfReader`` allows). Step 3: Register in the Factory ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Add a case to ``create_reader()`` in ``cocoa_meteo/MeteoReaderFactory.cpp``: .. code-block:: cpp case MeteoFormat::YOUR_FORMAT: return std::make_unique(config); Step 4: Add Configuration Parsing ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ In ``src/cocoa_kernel/core/config/ForcingParsers.cpp``, add the format's string spellings to the ``METEO_FORMATS`` lookup table that ``parse_meteo_format()`` reads (bump the table's size template argument): .. code-block:: cpp constexpr Core::EnumTable METEO_FORMATS{ {{{"cf_netcdf", Meteo::MeteoFormat::CF_NETCDF}, {"cf", Meteo::MeteoFormat::CF_NETCDF}, {"netcdf", Meteo::MeteoFormat::CF_NETCDF}, {"owi_ascii", Meteo::MeteoFormat::OWI_ASCII}, {"owi", Meteo::MeteoFormat::OWI_ASCII}, {"owi_netcdf", Meteo::MeteoFormat::OWI_NETCDF}, {"grib", Meteo::MeteoFormat::GRIB}, {"grib2", Meteo::MeteoFormat::GRIB}, {"your_format", Meteo::MeteoFormat::YOUR_FORMAT}}}}; This table is consulted once per configured source, so a run with several ``domains:`` entries can mix ``your_format`` with any other format for free. Add any format-specific fields to ``MeteoReaderConfig`` if needed (e.g., custom variable names or file organization parameters), and parse them alongside the format string in the same file. Step 5: Update CMakeLists.txt ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Add the new source files to ``src/cocoa_meteo/CMakeLists.txt``: .. code-block:: text readers/your_format/YourFormatReader.hpp readers/your_format/YourFormatReader.cpp Key Data Structures ^^^^^^^^^^^^^^^^^^^ **MeteoGrid**: grid geometry returned by ``grid_info()``: .. code-block:: cpp // Regular rectangular grid struct RegularGrid { double xll, yll; // Lower-left corner (longitude, latitude) double dx, dy; // Grid spacing [degrees] size_t nx, ny; // Grid dimensions bool global; // True if 360-degree longitude wrap }; // Irregular grid (explicit coordinates per point) struct IrregularGrid { size_t nx, ny; MeteoField longitude; // [ny * nx] coordinate field MeteoField latitude; // [ny * nx] coordinate field }; using MeteoGrid = std::variant; Use ``RegularGrid`` for structured data on a latitude-longitude grid. Use ``IrregularGrid`` when grid coordinates vary per point (e.g., moving vortex grids). **MeteoField**: 2D field storage ``[ny * nx]``. Data is stored in row-major order: ``data[j * nx + i]`` where ``i`` is the longitude index and ``j`` is the latitude index. This matches the NetCDF ``[lat, lon]`` dimension convention. **Scale factors:** The ``MeteoReaderConfig`` provides ``pressure_scale`` and ``wind_scale`` fields. Apply them in your reader's ``read_snapshot()`` before returning, as every existing reader does, by calling ``apply_scale_factors()`` from ``core/MeteoFieldSnapshot.hpp``. Nothing downstream scales the snapshot, so a reader that returns raw values yields unscaled winds and pressures. Reference: Existing Readers ^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. list-table:: :header-rows: 1 :widths: 30 15 55 * - Reader - Domains - Notes * - ``CfNetcdfReader`` - Single - CF-compliant NetCDF with customizable variable names. Auto-detects global grids. * - ``OwiAsciiReader`` - Multiple - Paired pressure/wind ASCII files with fixed-width headers. Indexes file positions on open for random-access reads. * - ``OwiNetcdfReader`` - Multiple - NetCDF group-based format with support for moving grids (per-timestep coordinates). * - ``GribReader`` - Multiple - GRIB2 via ecCodes, gated behind ``cocoa_ENABLE_GRIB``. One domain per configured weather product (priority level); supports moving (storm-following) nests. * - ``CompositeMeteoReader`` - Flattened - Not a format reader; wraps two or more configured sources of any format(s) into one flat domain space, list order = coverage priority.