Extending the Model

This guide covers how to extend Cocoa with new model capabilities, including nodal attributes (spatially varying parameters) and physics contributions (terms in the GWCE and momentum equations).

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:

enum class NodalAttributeId : int {
  ManningsN = 0,
  Tau0 = 1,
  AppleSauce = 2,  // New attribute
};

/// Number of defined nodal attributes (update when adding new entries)
static constexpr size_t NUM_NODAL_ATTRIBUTES = 3;
static_assert(static_cast<int>(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:

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";
    case NodalAttributeId::AppleSauce:
      return "apple_sauce_coefficient";
    default:
      return "";
  }
}

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:

void NodalAttributeInitializer::initialize_apple_sauce() {
  constexpr auto default_value = 1.0_s;
  constexpr auto valid_min = 0.0_s;
  constexpr auto valid_max = 100.0_s;

  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():

void NodalAttributeInitializer::initialize() {
  m_mesh->initialize_attributes(m_fields.num_nodes());
  initialize_manning_n();
  initialize_tau0();
  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:

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
});

The get() method returns a const NodalAttributeAccessor& from a pre-built cache array — pure array indexing with zero overhead. If the attribute was not registered, get() throws with a clear error message.

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. Use the mesh sentinel value (Sentinels::MESH_ATTRIBUTE) to indicate the value should be read from the mesh file rather than using a constant.

Step 5: Update Tests

Register the new attribute in test fixtures that exercise the kernels using it:

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:

Reading: Use NodalAttributeReader::read_vector_attribute() or read_vector_attribute_distributed() for 2D variables [num_nodes, num_components]:

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");

Access: Use the two-argument accessor to get a specific component:

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
});

Auto-detecting dimensionality: If an attribute may be either scalar (1D in the mesh file) or vector (2D), use reader.get_variable_dimensions() to detect the shape and register accordingly. See initialize_internal_tide_friction() for an example.

Optional Attributes

Some attributes are optional — they should only be applied when enabled by the user. For these, skip registration when disabled and check has() before accessing:

// 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

The physics contribution system is built around two patterns:

GWCE contributions use nested Input structs constructed inline at the call site. Each contribution’s Input holds only the data that contribution needs. Data loaded for one contribution’s Input dies after compute() returns, reducing peak register pressure on GPU.

Momentum contributions use a shared MomentumElementContext struct.

Fig. 36 Physics contribution architecture for GWCE and Momentum

Concepts

Four C++20 concepts in numeric/terms/TermsConcepts.hpp define the contracts:

Concept	Required Signature	Use Case
GwcePhysicsTerm	static compute(const T::Input&) -> Types::Vec2<double>	Returns (jx, jy) contribution to both x- and y-equation GWCE RHS
GwceScalarTerm	static compute(const T::Input&) -> double	Returns a scalar GWCE RHS contribution (uniform across all nodes)
MomentumElementTerm	static compute(const MomentumElementContext&) -> Types::Vec2<double>	Returns uniform (x, y) contribution to both x- and y-momentum for all 3 nodes
MomentumNodalTerm	static compute_nodal(const MomentumElementContext&) -> Kokkos::Array<Types::Vec2<double>, 3>	Returns per-node (x, y) contributions to both x- and y-momentum

Step-by-Step Guide

Fig. 37 Workflow for adding a new physics contribution

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:

/// AppleSauce forcing sub-container
struct AppleSauceFields {
  ScalarView coefficient;
  ScalarView intensity;
};

Add a private member, allocate in the constructor, and expose with const/non-const accessors:

// 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/terms/AppleSauce.hpp. Each GWCE contribution nests an Input struct containing only the data it needs, and provides a static compute(const Input&) method:

#pragma once
#include <Kokkos_Core.hpp>
#include "numeric/terms/TermsConcepts.hpp"

namespace Cocoa::Numeric::GwceRhsTerms {

struct GwceAppleSauce {
  struct Input {
    Kokkos::Array<Types::Scalar, 3> apple_sauce;
    // Add only the fields this contribution needs
  };

  KOKKOS_INLINE_FUNCTION
  static Types::Vec2<double> compute(const Input& input) {
    const auto as_avg = (input.apple_sauce[0] + input.apple_sauce[1] +
                         input.apple_sauce[2]) / 3.0_s;
    // Physics math goes here
    const auto jx = /* ... */;
    const auto jy = /* ... */;
    return {jx, jy};
  }
};

static_assert(GwcePhysicsTerm<GwceAppleSauce>);

}  // namespace Cocoa::Numeric::GwceRhsTerms

The static_assert verifies the concept contract at compile time. If your method signature doesn’t match, you get a clear error instead of a cryptic template failure at the call site.

The nested Input pattern ensures that data loaded for this contribution dies when compute() returns, freeing registers for subsequent contributions.

Step 3: Momentum Contribution

Add the momentum term to numeric/terms/AppleSauce.hpp (alongside the GWCE term). Choose the concept based on whether the result is uniform across all 3 nodes or varies per node.

Uniform across nodes – follow MomentumAdvection.hpp:

struct MomentumAppleSauce {
  KOKKOS_INLINE_FUNCTION
  static Types::Vec2<double> compute(const MomentumElementContext& ctx) {
    // Physics math: same value scattered to all 3 nodes
    return {/* x */, /* y */};
  }
};
static_assert(MomentumElementTerm<MomentumAppleSauce>);

Different per node – follow MomentumLateralStress.hpp:

struct MomentumAppleSauce {
  KOKKOS_INLINE_FUNCTION
  static Kokkos::Array<Types::Vec2<double>, 3> compute_nodal(
      const MomentumElementContext& ctx) {
    // Physics math: different value per node
    return {/* node 0 */, /* node 1 */, /* node 2 */};
  }
};
static_assert(MomentumNodalTerm<MomentumAppleSauce>);

Step 4: Wire Into GWCE Kernel

In numeric/continuity/GwceVectorAssemblyKernels.hpp, make two changes to main_rhs_assembly_kernel:

a) Include the header and capture the view at the top of the function, alongside the other physics view captures:

// At file top:
#include "numeric/terms/AppleSauce.hpp"

// Inside the function, alongside other view captures:
const auto apple_sauce = physics.apple_sauce().coefficient;

b) Call the contribution inside the lambda, constructing the Input inline with aggregate initialization and accumulating into the local result:

jxy += GwceAppleSauce::compute(
    {load_nodes(apple_sauce)});

Add the using declaration at the top of the lambda:

using GwceRhsTerms::GwceAppleSauce;

That’s it–one line to capture the view, one line to call. The Input temporary (and any data loaded into it) dies after compute() returns. The result is accumulated into a Types::Vec2<double> jxy variable.

Step 5: Wire Into Momentum Kernel

In numeric/momentum/MomentumRhsKernels.cpp, same three-part pattern in momentum_rhs_element_kernel:

a) Include and capture:

#include "numeric/terms/AppleSauce.hpp"

// Capture the view
const auto apple_sauce = physics.apple_sauce().coefficient;

b) Load into context:

ctx.apple_sauce = {apple_sauce(ctx.nm[0]), apple_sauce(ctx.nm[1]),
                   apple_sauce(ctx.nm[2])};

c) Call the contribution. For MomentumElementTerm, add to the base contribution:

const auto apple = MomentumAppleSauce::compute(ctx);
const auto base_contrib = pressure + advect + apple;

For MomentumNodalTerm, add in the per-node loop alongside lateral stress:

const auto apple = MomentumAppleSauce::compute_nodal(ctx);
for (int i = 0; i < nodes_per_triangle; ++i) {
  const auto contrib = base_contrib + lateral[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 Adding a Nodal Attribute above.

Computed each timestep (like lateral stress):

Write a computation kernel that populates AppleSauceFields. Call it from GwceVectorAssembler::assemble during the preprocessing phase:

// In GwceVectorAssembler::assemble_rhs, preprocessing block:
{
  const auto preprocess_region =
      Core::KokkosProfileRegion("GWCE::Preprocessing");
  nodal_preprocessing_kernel(model_fields);
  elemental_preprocessing_kernel(model_fields, config);
  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 Exchange (if needed)

Ghost exchange is only needed if the field is assembled from element contributions via atomics. In that case, ghost nodes only see contributions from local elements and miss contributions from elements on their owning rank.

Fields that are directly loaded from nodal attributes or computed per-node (without element scatter) do not need ghost exchange.

Add a method to CommunicationManager:

// In CommunicationManager.hpp (public):
void post_apple_sauce_exchange(Data::PhysicsData& physics);

// In CommunicationManager.cpp:
void CommunicationManager::post_apple_sauce_exchange(
    Data::PhysicsData& physics) {
  if (!m_is_parallel) return;
  std::array fields = {
    Types::StridedScalar1DView(physics.apple_sauce().intensity)
  };
  exchange_scalar_fields(fields);
}

Call from TimeStepper after the AppleSauce computation kernel, before the GWCE/momentum assembly that reads the data.

Note

Use the batched exchange_scalar_fields with std::array to exchange multiple fields in a single MPI round-trip rather than calling exchange_scalar_field once per view.

Step 8: Configuration and Build

Configuration – Add parameters to core/config/PhysicsConfig.hpp:

struct AppleSauceConfig {
  bool active{false};
  Types::Scalar coefficient{1.0};
};

Add AppleSauceConfig apple_sauce; to PhysicsConfig. Parse in core/config/ModelConfigurationFactory.cpp.

CMakeLists.txt – Add new headers to COCOA_KERNEL_HEADERS in src/cocoa_kernel/CMakeLists.txt:

numeric/terms/AppleSauce.hpp

Reference: Existing Contributions

File / Struct	Concept	Description
Coriolis.hpp – GwceCoriolis	GwcePhysicsTerm	Coriolis force (jx, jy) with spherical correction
BottomFrictionTerm.hpp – GwceBottomFriction	GwcePhysicsTerm	TKM bottom friction tensor (jx, jy)
Tau0.hpp – GwceTau0	GwcePhysicsTerm	Spatially variable tau0 weighting (jx, jy)
TidePotential.hpp – GwceTidePotential	GwcePhysicsTerm	Tide potential gradient (jx, jy)
WindStressTerm.hpp – GwceWindStress	GwcePhysicsTerm	Kinematic wind stress (jx, jy) via wind_stress<DragLaw>()
AtmPressureGradient.hpp – GwceAtmPressureGradient	GwcePhysicsTerm	Atmospheric pressure gradient (jx, jy)
LateralStressTerm.hpp – GwceLateralStress	GwcePhysicsTerm	Lateral stress divergence (jx, jy)
GwceAdvection.hpp – GwceAdvection	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)
LateralStressTerm.hpp – MomentumLateralStress	MomentumNodalTerm	Inline lateral stress (x, y) (per-node)

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.

Architecture Overview

Fig. 38 cocoa_meteo reader architecture

Step 1: Add the Format Enum

Add an entry to MeteoFormat in cocoa_meteo/core/MeteoFormat.hpp:

enum class MeteoFormat { CF_NETCDF, OWI_ASCII, OWI_NETCDF, 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:

#pragma once
#include "cocoa_meteo/core/MeteoReaderBase.hpp"
#include "cocoa_meteo/core/MeteoReaderConcept.hpp"
#include "cocoa_meteo/core/MeteoReaderConfig.hpp"

namespace Cocoa::Meteo {

class YourFormatReader : public MeteoReaderBase {
 public:
  explicit YourFormatReader(const MeteoReaderConfig& config);

  void open() override;
  void close() override;
  [[nodiscard]] auto is_open() const -> bool override;
  [[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<double>& override;
  [[nodiscard]] auto read_snapshot(size_t snapshot_idx,
                                   size_t domain_idx)
      -> MeteoFieldSnapshot override;

 private:
  // Format-specific internal state
};

// Verify concept compliance at compile time
static_assert(MeteoReader<YourFormatReader>);

}  // namespace Cocoa::Meteo

The pure virtual methods define the contract:

Method	Description
open()	Open file(s) and index snapshot positions
close()	Release file handles
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) for one time step and one domain.

Step 3: Register in the Factory

Add a case to create_reader() in cocoa_meteo/MeteoReaderFactory.cpp:

case MeteoFormat::YOUR_FORMAT:
  return std::make_unique<YourFormatReader>(config);

Step 4: Add Configuration Parsing

In cocoa_core/config/ModelConfigurationFactory.cpp, add the string-to-enum mapping for the new format:

if (format_str == "your_format") {
  reader_config.format = Meteo::MeteoFormat::YOUR_FORMAT;
}

Add any format-specific fields to MeteoReaderConfig if needed (e.g., custom variable names or file organization parameters).

Step 5: Update CMakeLists.txt

Add the new source files to src/cocoa_meteo/CMakeLists.txt:

readers/your_format/YourFormatReader.hpp
readers/your_format/YourFormatReader.cpp

Key Data Structures

MeteoGrid — Grid geometry returned by grid_info():

// 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<RegularGrid, IrregularGrid>;

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 these in your reader’s read_snapshot() or let the caller apply them — the convention is to return raw values and let the interpolation pipeline handle scaling.

Reference: Existing Readers

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).