NumPy and HDF5 output

Ragged dispersion rows

Each frequency and azimuth can produce a different number of valid modes. SpecSWD represents phase and group results with independent CSR containers:

  • indptr: uint64, shape (nsolve + 1,);

  • mode_index: int32, shape (nentry,);

  • values: complex128, shape (nentry, ncomponent); and

  • component_names: labels for the final dimension.

Solve rows are frequency-major and azimuth-minor:

row = frequency_index * table.azimuth_deg.size + azimuth_index
start, stop = csr.indptr[row:row + 2]
modes = csr.mode_index[start:stop]
values = csr.values[start:stop]

An empty row has equal adjacent offsets. Independent phase and group CSR tables allow their mode counts and mode selections to differ at the same solve point.

Python API

SpecWorkSpace.compute_dispersion returns NumPy-backed data directly:

import numpy as np

table = ws.compute_dispersion(
    frequencies_hz=np.geomspace(0.02, 1.0, 50),
    azimuths_deg=np.array([0.0, 30.0, 60.0]),
    include_group=True,
)
table.to_hdf5("dispersion.h5")

restored = DispersionTable.from_hdf5("dispersion.h5")
modes, phase = restored.phase_row(4, 1)
group_modes, group_values = restored.group_row(4, 1)

Use DispersionTable.from_rows to pack existing variable-length NumPy arrays:

from specd import DispersionTable

table = DispersionTable.from_rows(
    frequency_hz=np.array([0.05, 0.1, 0.2]),
    phase_rows=[
        np.array([3.0 + 0.01j, 3.4 + 0.02j]),
        np.array([3.1 + 0.01j]),
        np.empty(0, dtype=np.complex128),
    ],
    phase_mode_rows=[np.array([0, 2]), np.array([0]), np.empty(0, dtype=np.int32)],
)

The complete command-line result is available through SolverOutput:

from specd import SolverOutput

result = SolverOutput.from_hdf5("out/specswd.h5")
solve_row = result.dispersion.solve_row(0, 0)
mode_start = int(result.dispersion.phase_velocity.indptr[solve_row])

depth = result.sem_depth_row(solve_row)
eigenfunction = result.eigenfunctions.row(mode_start)
elastic_kernel = result.elastic_kernels.velocity[mode_start]

eigenfunction has shape (n_sem_point, ncomponent). A kernel table has shape (nphase_entry, nparameter, n_model_depth). Attenuating runs also provide propagation_q_inverse with the same shape. All public numeric containers are contiguous NumPy arrays; xarray and NumPy object arrays are not used. For VTI Rayleigh and Scholte waves, horizontal and vertical displacement components use real polarization amplitudes. The physical quadrature phase of the vertical component is implicit rather than stored as a factor of i.

Command-line output

The solvers write one self-describing file and replace it on each run:

surflove model.txt 0.02 1.0 50 0
surfrayl model.txt 0.02 1.0 50 1
surfani model.txt 0.02 1.0 50 30.0 1

The last value selects phase kernels (0) or group kernels (1). Anisotropic group velocity is saved with x and y components. Its group kernel is the radial derivative. Each command writes out/specswd.h5.

File schema

The implemented schema is specswd-dispersion, version 1:

/
  attrs: schema, schema_version, layout, complex_encoding
  attrs for complete output: content, wave_type, attenuation, kernel_type
  /solutions
    frequency_hz                 float64 (nfreq), units="Hz"
    azimuth_deg                  float64 (nazimuth), units="degree"
    /phase
      indptr                     uint64 (nsolve+1)
      mode_index                 int32 (nphase)
      component_names            ["phase"]
      values                     complex (nphase,1), units="km/s"
    /group                       optional, independent CSR
      indptr                     uint64 (nsolve+1)
      mode_index                 int32 (ngroup)
      component_names            ["radial"] or ["x","y"]
      values                     complex (ngroup,ncomponent), units="km/s"
  /model                         complete output only
    depth                        float64 (n_model_depth), units="km"
  /mesh
    depth_indptr                 uint64 (nsolve+1)
    depth_values                 float64 (n_sem_point_total), units="km"
  /eigenfunctions
    indptr                       uint64 (nphase+1)
    component_names              wave-dependent component labels
    values                       complex (n_eigen_point_total,ncomponent)
  /kernels
    attrs: observable, entry_alignment="solutions/phase"
    /elastic                     optional
      parameter_names
      velocity                   float64 (nphase,nparameter,n_model_depth)
      propagation_q_inverse      optional, same shape
    /acoustic                    optional, same datasets

depth_indptr handles a different SEM mesh at every solve point. eigenfunctions/indptr handles the corresponding different eigenfunction length for every phase entry. Projected kernels use the fixed input model depth grid and therefore remain dense after the ragged phase-entry dimension has been packed.

HDF5 has no portable native complex type. The file stores complex values as a compound with float64 members r and i. Python exposes this representation as numpy.complex128; C++ converts it to and from specswd::Complex.

C++ API

IO is declared in shared/io.hpp and is part of the shared target:

#include "shared/io.hpp"

using specswd::Complex;
using specswd::io::ComplexCSR;
using specswd::io::DispersionTable;

DispersionTable table;
table.frequency_hz = {0.05, 0.1, 0.2};
table.azimuth_deg = {0.0};
table.phase_velocity = ComplexCSR::from_rows(
    {{{3.0, 0.01}, {3.4, 0.02}}, {{3.1, 0.01}}, {}},
    {"phase"},
    {{0, 2}, {0}, {}}
);

specswd::io::write_hdf5("dispersion.h5", table);
auto restored = specswd::io::read_hdf5("dispersion.h5");
auto [begin, end] = restored.phase_velocity.row_range(1);

SolverOutput, RaggedComplexField, and KernelTable provide the complete output form used by the executables. The implementation is in src/shared/io.cpp and calls the HDF5 C API directly. Linking shared propagates the required HDF5 include and library settings.

Both implementations validate offsets, shapes, mode indices, coordinate counts, units, schema metadata, eigenfunction-to-depth alignment, and kernel dimensions before accepting or writing data.