Trilinos Integration

Cocoa uses Trilinos for distributed linear algebra and iterative solvers.

Trilinos Type Aliases

core/types/LinearAlgebraTypes.hpp defines the Trilinos types on Cocoa’s Kokkos execution space:

LinearAlgebraTypes

namespace Cocoa::Types {

struct LinearAlgebraTypes {
  // Node type wraps Kokkos execution space for Tpetra
  using ExecutionNode =
      Tpetra::KokkosCompat::KokkosDeviceWrapperNode<CocoaExecutionSpace>;

  // Core Tpetra types parameterized by ordinal types and node
  using Map = Tpetra::Map<LocalOrdinal, GlobalOrdinal, ExecutionNode>;
  using Vector = Tpetra::Vector<double, LocalOrdinal, GlobalOrdinal, ExecutionNode>;
  using MultiVector = Tpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, ExecutionNode>;
  using Import = Tpetra::Import<LocalOrdinal, GlobalOrdinal, ExecutionNode>;
  using Export = Tpetra::Export<LocalOrdinal, GlobalOrdinal, ExecutionNode>;

  // Sparse matrix types
  using CrsGraph = Tpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, ExecutionNode>;
  using CrsMatrix = Tpetra::CrsMatrix<double, LocalOrdinal, GlobalOrdinal, ExecutionNode>;

  // Operator interface for solvers
  using Operator = Tpetra::Operator<double, LocalOrdinal, GlobalOrdinal, ExecutionNode>;
};

}

LinearSolverTypes

namespace Cocoa::Types {

struct LinearSolverTypes {
  // Belos template parameters
  using SC = double;
  using MV = LinearAlgebraTypes::MultiVector;
  using OP = LinearAlgebraTypes::Operator;

  // Ifpack2 preconditioner
  using PreconditionerType = Ifpack2::Preconditioner<
      LinearAlgebraTypes::CrsMatrix::scalar_type,
      LinearAlgebraTypes::CrsMatrix::local_ordinal_type,
      LinearAlgebraTypes::CrsMatrix::global_ordinal_type,
      LinearAlgebraTypes::CrsMatrix::node_type>;

  // Belos solver types
  using LinearProblem = Belos::LinearProblem<SC, MV, OP>;
  using SolverManager = Belos::SolverManager<SC, MV, OP>;
};

}

GWCE Matrix Assembly

The GWCE matrix is a plain Tpetra::CrsMatrix of the owned rows. Its pattern is fixed for the run, so it is fill-completed once at construction and never reopened: each step overwrites its values in place. Reopening a matrix (resumeFill) discards the state Tpetra builds for the solve’s sparse products, which then has to be rebuilt, with a device allocation and a synchronizing free, on every step.

Graph Construction

From GwceSolverConsistent.cpp. Each rank inserts its elements’ blocks into the rows of every local node, owned and ghost, and exports them to the owned rows, so a row on a partition boundary gets the union of every rank’s columns:

// Local indices; the row bound is the largest node valence plus one
const auto row_bound = mesh.device().max_node_valence() + 1;
LinearAlgebraTypes::CrsGraph local_blocks(overlap_map, overlap_map, row_bound);
for (const auto& element : mesh.host().elements<Triangle>()) {
  insert_dense_block(local_blocks, element.nodes());
}
// ... quadrilaterals the same way
local_blocks.fillComplete(overlap_map, overlap_map);

auto graph = Teuchos::rcp(
    new LinearAlgebraTypes::CrsGraph(owned_map, column_map, row_bound));
graph->doExport(local_blocks,
                LinearAlgebraTypes::Export(overlap_map, owned_map),
                Tpetra::INSERT);
graph->fillComplete(owned_map, owned_map);

Matrix Assembly Pattern

Elements have no ghosts, so a rank’s elements contribute to its owned rows and to its ghost nodes’ rows. The owned rows are the matrix; the ghost rows are a small CSR held by SharedRowExchange. AssemblyRows presents the two to the element kernels as one target with one column numbering (the local node index), so the kernels never ask which half a row is in:

// GwceMatrixAssemblerConsistent::assemble
Kokkos::deep_copy(space,
    matrix->getLocalValuesDevice(Tpetra::Access::OverwriteAll), 0.0);
const AssemblyRows rows{.owned = matrix->getLocalMatrixDevice(),
                        .ghost = shared_rows.ghost_rows()};
Kokkos::deep_copy(space, rows.ghost.values, 0.0);

// Element-parallel scatter over each type's wet list: a triangle adds each
// of its rows with rows.sumIntoValues(...) (atomic); a quadrilateral adds
// through slots found once at construction (ElementAssemblyMap).
assemble_gwce_consistent_matrix_kernel<Triangle>(fields, mesh, rows, ...);
assemble_gwce_consistent_matrix_kernel<Quadrilateral>(fields, mesh, rows, ...);

// Move open-boundary columns out, on the rank that computed each entry.
eliminate_bc_columns(rows, obc_coef);

// Send the ghost rows to their owners and add them in.
shared_rows.add_into(*matrix);

// Owned rows are now complete: compute EP, then pin BC and dry rows.
compute_elevation_penalty(matrix, diagonal_offsets, num_owned_nodes, comm, ep);
apply_dirichlet_bc_kernel(matrix, diagonal_offsets, obc_coef, ep.value,
                          dry_rows, num_owned_nodes);

The exchange plan is built once. Every ghost-row entry is named by the global id row gid * num_global_nodes + column gid and moved by a Tpetra::Export with ADD from the senders’ overlapping map to the owners’ one-to-one map; the owner adds each received value at its CSR position, found at construction. Serially there are no ghost rows and the exchange does nothing.

Boundary Condition Application

Each row’s diagonal position is found once, at construction, with CrsGraph::getLocalDiagOffsets (DiagonalOffsets.hpp); construction fails if an owned row has no diagonal, which means a node in no element. Elevation-specified (Dirichlet) rows, and rows no wet element reached, are pinned to the elevation penalty in one pass over the owned rows, in apply_dirichlet_bc_kernel:

Kokkos::parallel_for("Apply BC and Fix Dry Nodes",
    Types::DefaultRangePolicy(0, num_owned_rows),
    KOKKOS_LAMBDA(const int row) {
      const bool dry = row_is_dry(local_matrix, diagonal_offsets, row, eps);
      dry_rows(row) = dry ? 1 : 0;  // the RHS reads the same verdict
      if (dry or bc_flags(row) == 1) {
        // Zero the row, then write EP on its diagonal
        pin_row_to_diagonal(local_matrix.row(row),
                            diagonal_index(diagonal_offsets, row),
                            ep_value());
      }
    });

The boundary columns of every other row were already moved into OpenBoundaryCoefficients by eliminate_bc_columns, before the shared-row exchange, so the right-hand side can apply them.

Linear Solver Configuration

Cocoa solves the GWCE system with a Jacobi-preconditioned Conjugate Gradient solver, configured in ConjugateGradientSolver.hpp. The default implementation is standard Belos CG (CgSolverType::BelosCG). Two Tpetra-native implementations are selectable via CgSolverType: TpetraSingleReduce (single-reduce CG) and TpetraCgPipeline (pipelined CG). They need fewer global reductions per iteration, but the single-reduce recurrence was measured to take more iterations per solve and to be slower overall, which is why Belos CG is the default.

Solver Configuration

The SolverConfig struct controls solver behavior:

struct SolverConfig {
  CgSolverType solver_type{CgSolverType::BelosCG};  // CG implementation
  double tolerance{1e-5};            // Convergence tolerance
  int max_iterations{50};            // Maximum CG iterations
  int verbosity{0};                  // 0=silent, 1=summary, 2=detailed
  ResidualScaling residual_scaling{ResidualScaling::NormOfInitialResidual};  // Convergence scaling
};

Residual Scaling Options

The ResidualScaling enum controls how residual norms are scaled for convergence testing:

enum class ResidualScaling {
  NormOfInitialResidual,  // ||r|| / ||r_0|| < tol (default)
  NormOfPrecInitRes,      // ||M^{-1}r|| / ||M^{-1}r_0|| < tol
  None                    // ||r|| < tol (absolute)
};

NormOfInitialResidual is the default; at similar tolerances it gives iteration counts comparable to ADCIRC’s ITPACK JCG solver.

Solver Setup

// Create linear problem
auto problem = Teuchos::rcp(new LinearSolverTypes::LinearProblem(
    matrix, solution, rhs));

// Configure solver parameters
auto params = Teuchos::rcp(new Teuchos::ParameterList());
params->set("Maximum Iterations", config.max_iterations);  // Default: 50
params->set("Convergence Tolerance", config.tolerance);    // Default: 1e-5
params->set("Implicit Residual Scaling",
            "Norm of Initial Residual");                   // Default scaling

// Verbosity (Belos CG): errors always, then more with each level
int verbosity = Belos::Errors;
if (config.verbosity >= 1) {
  verbosity |= Belos::Warnings | Belos::FinalSummary;
}
if (config.verbosity >= 2) {
  verbosity |= Belos::StatusTestDetails;
}
params->set("Verbosity", verbosity);

// Create the solver: "CG" for BelosCG, or a Tpetra-native CG
Belos::SolverFactory<double, MultiVector, Operator> factory;
auto solver = factory.create(solver_name(), params);

Jacobi Preconditioner

From JacobiPreconditioner.hpp:

// Create Ifpack2 preconditioner via the static factory
auto prec = Ifpack2::Factory::create<Matrix>("RELAXATION", matrix);

// Configure as Jacobi
const auto prec_params = Teuchos::rcp(new Teuchos::ParameterList());
prec_params->set("relaxation: type", "Jacobi");
prec_params->set("relaxation: sweeps", 1);
prec_params->set("relaxation: damping factor", 1.0);
prec->setParameters(*prec_params);

// Initialize sparsity pattern
prec->initialize();

// Compute values (must be called when matrix changes)
prec->compute();

// Set on linear problem
problem->setRightPrec(prec);

Solving

// Set up the problem
problem->setProblem();

// Solve
auto result = solver->solve();

if (result != Belos::Converged) {
  // Handle non-convergence
}

Vector Operations

Tpetra vectors integrate with Kokkos views for device-side access:

// Create vector
auto vec = Teuchos::rcp(new LinearAlgebraTypes::Vector(map));

// Get device view for parallel operations
auto view = vec->getLocalViewDevice(Tpetra::Access::ReadWrite);

// Use in Kokkos kernel
Kokkos::parallel_for("Set values",
    Types::DefaultRangePolicy(0, num_nodes),
    KOKKOS_LAMBDA(int i) {
      view(i, 0) = initial_value(i);
    });

Lumped Diagonal Solve:

The lumped GWCE solver has the same shape. It takes the local device views and divides in one pass, in diagonal_solve_kernel (numeric/cg/continuity/lumped/GwceSolverLumpedKernels.cpp):

auto sln_1d = Kokkos::subview(sln->getLocalViewDevice(Tpetra::Access::ReadWrite),
                              Kokkos::ALL(), 0);
// ... rhs_1d, diag_1d likewise, ReadOnly ...
Kokkos::parallel_for(
    "GWCE Lumped Diagonal Solve",
    Types::DefaultRangePolicy(0, num_local_nodes),
    KOKKOS_LAMBDA(const int node) {
      sln_1d(node) = rhs_1d(node) / diag_1d(node);
    });

The division is unguarded: the lumped diagonal is the nodal mass, which the assembly leaves strictly positive. Only local rows are touched, so the pass needs no communication.

Elevation Penalty Computation

The elevation penalty (EP) maintains numerical conditioning for BC rows. It is the root mean square of the owned rows’ diagonals, read through the diagonal offsets, from compute_elevation_penalty:

auto local_matrix = matrix->getLocalMatrixDevice();

// Owned rows only, so no node is counted twice across ranks
Kokkos::parallel_reduce("DiagonalRMS",
    Types::DefaultRangePolicy(0, num_owned_nodes),
    KOKKOS_LAMBDA(const int row, double& local_sum) {
      const double diagonal = local_matrix.row(row).value(
          diagonal_index(diagonal_offsets, row));
      local_sum += diagonal * diagonal;
    },
    penalty.sum_sq);  // device-resident scalar

// sqrt(sum / count): partitioned runs sum both over all ranks first
finalize_elevation_penalty(num_owned_nodes, comm, penalty);

Implementation Files

• core/types/LinearAlgebraTypes.hpp: Type aliases

• numeric/cg/continuity/GwceSolver.hpp: Main GWCE solver class

• numeric/cg/continuity/consistent/ConjugateGradientSolver.hpp: Belos CG wrapper

• numeric/cg/continuity/consistent/JacobiPreconditioner.hpp: Ifpack2 Jacobi wrapper

• numeric/cg/continuity/consistent/GwceSolverConsistent.cpp: Graph and matrix construction

• numeric/cg/continuity/consistent/GwceMatrixAssemblerConsistent.hpp: Per-step matrix assembly

• numeric/cg/continuity/consistent/GwceMatrixAssemblyConsistentKernels.hpp: Matrix assembly kernels

• numeric/cg/continuity/consistent/SharedRowExchange.hpp: Ghost rows summed into their owners

• numeric/cg/continuity/ElementAssemblyMap.hpp: AssemblyRows and the quadrilateral slot map

• numeric/cg/continuity/DiagonalOffsets.hpp: Row diagonal positions