Issue 1982 - First stiffness matrix computation by lsawade · Pull Request #1998 · PrincetonUniversity/SPECFEMPP

lsawade · 2026-07-17T12:05:43Z

Description

This implements the computation of the stiffness kernel in linear systems

Currently the elastic tag is somewhat hardcoded, which will be removed in subsequent PRs where we call the stiffness computation separately.

Issue Number

Closes #1982

Checklist

Please make sure to check developer documentation on specfem docs.

I ran the code through pre-commit to check style
THE DOCUMENTATION BUILDS WITHOUT WARNINGS/ERRORS
I have added labels to the PR (see right hand side of the PR page)
My code passes all the integration tests
I have added sufficient unittests to test my changes
I have added/updated documentation for the changes I am proposing
I have updated CMakeLists to ensure my code builds
My code builds across all platforms

Copilot

Pull request overview

Implements first pass of explicit stiffness-kernel extraction for SPECFEM++ linear systems (specfem::linear_system) by probing the existing matrix-free element operator to build dense per-element stiffness blocks, with scope validation and accompanying unit tests/docs.

Changes:

Add compute_element_stiffness (NGLL-dispatched) and validate_stiffness_scope for 3D elastic isotropic, attenuation-none, natural-boundary meshes.
Introduce a Kokkos probe kernel that applies gradient → stress → divergence to unit local DOFs to populate dense element blocks.
Add a unit test that checks symmetry/nullspace properties and matches against the production matrix-free stiffness operator; expose the new API in Sphinx/Doxygen docs.

Reviewed changes

Copilot reviewed 9 out of 9 changed files in this pull request and generated 2 comments.

Show a summary per file

File	Description
`tests/unit-tests/serial.cmake`	Adds the new `element_stiffness_tests` target to the serial unit-test suite.
`tests/unit-tests/linear_system/element_stiffness_tests.cpp`	New tests validating stiffness-block properties and agreement with matrix-free operator.
`docs/sections/api/specfem/linear_system/index.rst`	New API doc entry for `specfem::linear_system`.
`docs/sections/api/specfem/index.rst`	Adds `linear_system` to the top-level API index.
`core/specfem/linear_system/impl/stiffness_probe_kernel.hpp`	New Kokkos kernel to probe the element stiffness operator into dense blocks.
`core/specfem/linear_system/element_stiffness.hpp`	Declares local DOF ordering helper, scope validator, and stiffness-block computation APIs.
`core/specfem/linear_system/element_stiffness.tpp`	Implements the NGLL-specific stiffness probe execution and buffer validation.
`core/specfem/linear_system/element_stiffness.cpp`	Implements runtime NGLL dispatcher + scope validation and provides explicit instantiations.
`cmake/trilinos.cmake`	Adjusts macOS rpath handling for Trilinos shared libraries.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

+  const auto &element_types = assembly.element_types;
+  const auto elements = element_types.get_elements_on_host(Tags::medium_tag);
+
+  for (int i = 0; i < elements.size(); ++i) {


+    if(APPLE)
+        get_filename_component(SPECFEM_TRILINOS_LIBRARY_DIR
+            "${Trilinos_DIR}/../.." ABSOLUTE)
+        list(APPEND CMAKE_BUILD_RPATH "${SPECFEM_TRILINOS_LIBRARY_DIR}")
+        list(APPEND CMAKE_INSTALL_RPATH "${SPECFEM_TRILINOS_LIBRARY_DIR}")
+    endif()


icui

Can the static_assert by replaced by requires?

#1984) Foundations for the implicit Newmark solver: - Link Belos and Ifpack2 into specfem_linear_system (MueLu stays deferred: the float-only TROMP install does not instantiate Xpetra::Matrix<double>). linear_solver_smoke_test() gates the toolchain: a tridiagonal solve through Ifpack2-RILUK + Belos PseudoBlockGmres on type_real. - StiffnessScope: the strict scope validator gains a with_stacey opt-in -- the displacement probe runs at zero velocity, where the Stacey dashpot contributes nothing to K; callers must assemble C separately. - assemble_mass_vector(): lumped diagonal mass M as a Tpetra::Vector via the production initialize_mass_matrix path with dt = 0 (the Stacey (dt/2) C 1 lumped term is linear in dt and vanishes exactly). - DampingAssembler: Stacey damping C by probing the velocity path of the production stiffness kernel (u = 0, v = e_c at all points; C is block-diagonal so ncomp launches recover all blocks), assembled on a compact graph whose entries all exist in K's graph. Unit tests cover positivity and rho*V total mass of M, block structure / symmetry / PSD of C, agreement of C v with the production kernel for a random velocity, and the cross-path identity M(dt) - M(0) = (dt/2) C 1 through the explicit solver's mass path. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

specfem::solver::ImplicitNewmarkSolver, adjacent to time_marching: one Belos GMRES solve per time step on the assembled displacement-form operator A = M/(beta dt^2) + gamma/(beta dt) C + K, built from the linear_system assemblers (Stacey-tolerant stiffness scope, velocity-path damping probe, dt=0 lumped mass). A lives on K's static graph, is assembled and preconditioned once per time step size, and each step is a warm-started solve with an Ifpack2 RILUK right preconditioner. RILUK defaults to zero fill: spectral-element rows carry ~2000 nonzeros, so level-1 fill (~the pattern of A^2) explodes combinatorially in Ifpack2's symbolic setup. MueLu stays deferred (float-only cluster installs). Run with the dissipative Newmark preset (gamma > 1/2, beta = (gamma + 1/2)^2/4) and the velocity-based steady-state stopping criterion, the solver acts as a static solver: it recreates an explicit run driven to steady state in a few large steps. The criterion is velocity/acceleration-based because a nonzero-net-force source on a Stacey box converges to a constant-velocity drift, not a displacement fixed point. Tests: constructor scope rejections, exact operator identity against independently applied K/C/M, and a dynamic-equivalence integration test against the explicit solver on the natural-boundary fixture. The equivalence tolerances (0.15/0.35/0.55 relative L2 for u/v/a) are scheme-difference-bound: average acceleration elongates periods while central difference shortens them, and the measured gaps (0.10/0.27/0.48 at the fixture's coarse dt against a 1.25 Hz Ricker) are dominated by that dispersion, growing with derivative order. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

…1984) The ImplicitNewmark integration suite gets its own data tree mirroring Newmark/serial/dim3: each fixture carries the full provenance (fortran DATA + specfempp mesh files) and a database-only Snakefile that regenerates database.bin via the shared rule modules (cross-suite include of Newmark/specfem3d.smk; hoisting the modules is a follow-up). These are not trace-regression fixtures -- no traces/ and no tests.yaml entry; the reference is always a live explicit run. Steady-state test (the static solver of issue #1984): a Heaviside step force (hdur 4 s -- large enough that the implicit dt = 0.7 does not trigger the Heaviside hdur >= 5 dt clamp, which would silently give the two runs different forcings) on the Stacey halfspace. The implicit solver taking 30 dissipative steps (gamma 0.6) recreates the explicit solver taking 600 steps to the same final time T = 21 s within 5% relative L2 on the full displacement and velocity fields. Fields are compared at the same T because the box converges to a constant-velocity drift, not a displacement fixed point. The detector assertion is mechanism-only (fires at a coarse tolerance, reusing the solver instance and its assembled operators via the new set_steady_state_tolerance()); where a tight tolerance fires depends on the fixture's physical ring-down, which the solver now logs per step when the detector is armed. Solver: accept a non-Converged Belos return when the true residual b - A u meets the tolerance -- single-precision PseudoBlockGmres aborts on its internal loss-of-accuracy check while the solution is fine (observed: flagged at step 29 with true relative residual 2e-6 against a 1e-5 request). Also: Heaviside branch in the dim3 YAML source parser, dynamic-equivalence fixture copied into the suite tree, unit tests merged into one process (construction dominates), and the integration tests labeled `long` (exclude locally with ctest -LE long). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

@tparam

Addresses icui's review on #1998: the dim3 and attenuation-none constraints move from class-scope static_asserts to a C++20 requires clause on the template head, matching the existing pattern in medium/*/frechet_derivative.hpp. The assert message prose moves into the @tparam Tags documentation. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Issue 1982 - Phase 2, StiffnessAssembler

Rohit-Kakodkar

Minor changes.

Rohit-Kakodkar · 2026-08-19T16:56:57Z

+   * @return DofMap for the medium
+   */
+  template <specfem::element::medium_tag MediumTag>
+  static DofMap from_assembly(


Is this just a constructor?

Issue 1984 - Phase 1: Mass/damping assembly and Belos gate for the implicit solver

Issue 1984 - Phase 2: Implicit Newmark solver

Issue 1984 - Phase 3: Static-solver steady-state test and fixture suite

lsawade · 2026-09-01T19:38:09Z

Can the static_assert by replaced by requires?

address, @icui

codecov-commenter · 2026-09-01T22:26:30Z

Codecov Report

❌ Patch coverage is 65.31792% with 60 lines in your changes missing coverage. Please review.
✅ Project coverage is 60.99%. Comparing base (d0ee652) to head (1d560ab).

Files with missing lines	Patch %	Lines
core/specfem/linear_system/element_stiffness.cpp	41.37%	27 Missing and 7 partials ⚠️
core/specfem/linear_system/element_stiffness.tpp	40.54%	16 Missing and 6 partials ⚠️
core/specfem/source/dim3/source.tpp	0.00%	4 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##            devel    #1998      +/-   ##
==========================================
- Coverage   62.85%   60.99%   -1.87%     
==========================================
  Files         489      562      +73     
  Lines       20127    23456    +3329     
  Branches     2615     2922     +307     
==========================================
+ Hits        12651    14306    +1655     
- Misses       6621     8113    +1492     
- Partials      855     1037     +182

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

First prototype for the element stiffness matrix kernel

79c76ce

lsawade requested review from Rohit-Kakodkar, Copilot and icui July 17, 2026 12:05

lsawade added enhancement New feature or request new physics This labels the introduction of new physics to the code labels Jul 17, 2026

Copilot started reviewing on behalf of lsawade July 17, 2026 12:06 View session

Merge branch 'trilinos-integration' into issue-1982

4c3af48

Copilot AI reviewed Jul 17, 2026

View reviewed changes

Base automatically changed from trilinos-integration to devel July 17, 2026 13:09

lsawade added 3 commits July 17, 2026 16:26

StiffnessAssembler

4e918b3

Merge branch 'devel' into issue-1982

ae22584

Merge branch 'issue-1982' into issue-1982-phase-2

df0c6a5

Rohit-Kakodkar approved these changes Jul 28, 2026

View reviewed changes

icui approved these changes Jul 30, 2026

View reviewed changes

lsawade and others added 8 commits August 11, 2026 14:38

Merge branch 'devel' into issue-1982

bd08604

Updated the require statements for the functions as well.

9c15d00

Updated time limit

c195e87

Merge pull request #2002 from PrincetonUniversity/issue-1982-phase-2

1e74f69

Issue 1982 - Phase 2, StiffnessAssembler

Rohit-Kakodkar requested changes Aug 19, 2026

View reviewed changes

Rohit-Kakodkar reviewed Aug 19, 2026

View reviewed changes

Comment thread core/specfem/linear_system/dof_map.hpp

Rohit-Kakodkar and others added 4 commits August 31, 2026 09:29

Merge pull request #2023 from PrincetonUniversity/issue-1984-phase-1

ca2e114

Issue 1984 - Phase 1: Mass/damping assembly and Belos gate for the implicit solver

Merge pull request #2024 from PrincetonUniversity/issue-1984-phase-2

7c414c6

Issue 1984 - Phase 2: Implicit Newmark solver

Merge pull request #2025 from PrincetonUniversity/issue-1984

154b8ce

Issue 1984 - Phase 3: Static-solver steady-state test and fixture suite

Updated the Trilinos tests

97968cd

lsawade requested a review from Rohit-Kakodkar September 1, 2026 19:38

Merge branch 'devel' into issue-1982

1d560ab

Navigation Menu

Uh oh!

Issue 1982 - First stiffness matrix computation - #1998

Issue 1982 - First stiffness matrix computation#1998
lsawade wants to merge 18 commits into
develfrom
issue-1982

lsawade commented Jul 17, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

icui left a comment

Uh oh!

Rohit-Kakodkar left a comment

Uh oh!

Uh oh!

Rohit-Kakodkar Aug 19, 2026

Uh oh!

lsawade Sep 1, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

lsawade commented Sep 1, 2026

Uh oh!

codecov-commenter commented Sep 1, 2026 •

edited

Loading

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

Uh oh!

Conversation

lsawade commented Jul 17, 2026

Description

Issue Number

Checklist

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

icui left a comment

Choose a reason for hiding this comment

Uh oh!

Rohit-Kakodkar left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Rohit-Kakodkar Aug 19, 2026

Choose a reason for hiding this comment

Uh oh!

lsawade Sep 1, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

lsawade commented Sep 1, 2026

Uh oh!

codecov-commenter commented Sep 1, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

codecov-commenter commented Sep 1, 2026 •

edited

Loading