==================== 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 ^^^^^^^^^^^^^^^^^^ .. code-block:: cpp namespace Cocoa::Types { struct LinearAlgebraTypes { // Node type wraps Kokkos execution space for Tpetra using ExecutionNode = Tpetra::KokkosCompat::KokkosDeviceWrapperNode; // Core Tpetra types parameterized by ordinal types and node using Map = Tpetra::Map; using Vector = Tpetra::Vector; using MultiVector = Tpetra::MultiVector; using Import = Tpetra::Import; using Export = Tpetra::Export; // Sparse matrix types using CrsGraph = Tpetra::CrsGraph; using CrsMatrix = Tpetra::CrsMatrix; // Operator interface for solvers using Operator = Tpetra::Operator; }; } LinearSolverTypes ^^^^^^^^^^^^^^^^^ .. code-block:: cpp 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; using SolverManager = Belos::SolverManager; }; } 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: .. code-block:: cpp // 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()) { 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: .. code-block:: cpp // 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(fields, mesh, rows, ...); assemble_gwce_consistent_matrix_kernel(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``: .. code-block:: cpp 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: .. code-block:: cpp 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: .. code-block:: cpp 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 ^^^^^^^^^^^^ .. code-block:: cpp // 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 factory; auto solver = factory.create(solver_name(), params); Jacobi Preconditioner ^^^^^^^^^^^^^^^^^^^^^ From ``JacobiPreconditioner.hpp``: .. code-block:: cpp // Create Ifpack2 preconditioner via the static factory auto prec = Ifpack2::Factory::create("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 ^^^^^^^ .. code-block:: cpp // 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: .. code-block:: cpp // 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``): .. code-block:: 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``: .. code-block:: cpp 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