Add docstrings for the public mass matrix helpers

Also mark __check_dae_adaptivity as public, since it is imported by
BoundaryValueDiffEqMIRK and BoundaryValueDiffEqFIRK and the
ExplicitImports QA check requires explicit imports to be public, and
render all six docstrings on the internal interfaces docs page as
required by the API docs QA check.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Author: SebastianM-C

docs/src/devdocs/internal_interfaces.md

@@ -22,9 +22,13 @@ BoundaryValueDiffEqCore.__FastShortcutNonlinearPolyalg
 BoundaryValueDiffEqCore.__Fix3
 BoundaryValueDiffEqCore.__add_singular_term!
 BoundaryValueDiffEqCore.__any_sparse_ad
+BoundaryValueDiffEqCore.__apply_algebraic_constraint!
+BoundaryValueDiffEqCore.__apply_algebraic_constraint_oop!
+BoundaryValueDiffEqCore.__apply_mass_matrix!
 BoundaryValueDiffEqCore.__build_cost
 BoundaryValueDiffEqCore.__build_solution
 BoundaryValueDiffEqCore.__cache_trait
+BoundaryValueDiffEqCore.__check_dae_adaptivity
 BoundaryValueDiffEqCore.__concrete_kwargs
 BoundaryValueDiffEqCore.__concrete_solve_algorithm
 BoundaryValueDiffEqCore.__construct_internal_problem
@@ -36,6 +40,7 @@ BoundaryValueDiffEqCore.__extract_mesh
 BoundaryValueDiffEqCore.__extract_problem_details
 BoundaryValueDiffEqCore.__extract_u0
 BoundaryValueDiffEqCore.__flatten_initial_guess
+BoundaryValueDiffEqCore.__get_algebraic_indices
 BoundaryValueDiffEqCore.__get_bcresid_prototype
 BoundaryValueDiffEqCore.__get_non_sparse_ad
 BoundaryValueDiffEqCore.__has_initial_guess
@@ -52,6 +57,7 @@ BoundaryValueDiffEqCore.__needs_diffcache
 BoundaryValueDiffEqCore.__resize!
 BoundaryValueDiffEqCore.__restructure_sol
 BoundaryValueDiffEqCore.__split_kwargs
+BoundaryValueDiffEqCore.__subtract_mass_stage!
 BoundaryValueDiffEqCore.__tunable_part
 BoundaryValueDiffEqCore.__use_both_error_control
 BoundaryValueDiffEqCore.__vec

lib/BoundaryValueDiffEqCore/src/BoundaryValueDiffEqCore.jl

@@ -82,6 +82,7 @@ export BVPVerbosity, _process_verbose_param, DEFAULT_VERBOSE
     __vec_f, __vec_f!, __vec_so_bc, __vec_so_bc!, _sparse_like,
     __apply_mass_matrix!, __get_algebraic_indices, __subtract_mass_stage!,
     __apply_algebraic_constraint!, __apply_algebraic_constraint_oop!,
+    __check_dae_adaptivity,
     concrete_jacobian_algorithm, diff!, eval_bc_residual, eval_bc_residual!,
     get_dense_ad, interval, nodual_value, recursive_flatten, recursive_flatten!,
     recursive_flatten_twopoint!, recursive_unflatten!, safe_similar, _unwrap_val

lib/BoundaryValueDiffEqCore/src/utils.jl

@@ -940,18 +940,29 @@ end
     return nothing
 end
 
+"""
+    __apply_mass_matrix!(residᵢ, mass_matrix, tmp)
+
+Overwrite `residᵢ` with `mass_matrix * residᵢ` in place, using `tmp` as workspace. Used in
+defect estimation so that the defect compares `M * z′` against `f(z)` for mass-matrix
+problems. No-op for `UniformScaling` mass matrices, so plain ODEs pay no cost.
+"""
 @inline function __apply_mass_matrix!(residᵢ, mass_matrix::UniformScaling, tmp)
-    # Identity matrix - no modification needed
     return nothing
 end
 
 @inline function __apply_mass_matrix!(residᵢ, mass_matrix::AbstractMatrix, tmp)
-    # Apply M * residᵢ, using tmp as workspace
     mul!(tmp, mass_matrix, residᵢ)
     copyto!(residᵢ, tmp)
     return nothing
 end
 
+"""
+    __get_algebraic_indices(mass_matrix)
+
+Return the indices of the zero rows of `mass_matrix`, i.e. the algebraic equations of an
+index-1 DAE, or `nothing` when there are none (always for `UniformScaling`).
+"""
 @inline function __get_algebraic_indices(mass_matrix::UniformScaling)
     return nothing
 end
@@ -961,6 +972,13 @@ end
     return isempty(indices) ? nothing : indices
 end
 
+"""
+    __subtract_mass_stage!(res, mass_matrix, K_r, tmp)
+
+Subtract `mass_matrix * K_r` from `res` in place, using `tmp` as workspace. Turns the
+collocation stage residual `f(...) - K_r` into `f(...) - M * K_r`, i.e. collocation of
+`M * u′ = f(u, p, t)`.
+"""
 @inline function __subtract_mass_stage!(res, ::UniformScaling, K_r, tmp)
     res .-= K_r
     return nothing
@@ -972,6 +990,14 @@ end
     return nothing
 end
 
+"""
+    __apply_algebraic_constraint!(residᵢ, algebraic_indices, f!, yᵢ₊₁, p, t, tmp)
+
+Overwrite the algebraic rows of the continuity residual `residᵢ` with the DAE constraint
+residual evaluated at the right mesh point `yᵢ₊₁`, following the unprojected collocation
+approach for index-1 DAEs of Ascher & Spiteri (1994). `f!` is evaluated in place into
+`tmp`. No-op when `algebraic_indices === nothing`.
+"""
 @inline function __apply_algebraic_constraint!(
         residᵢ, ::Nothing, f!, yᵢ₊₁, p, t, tmp
     )
@@ -986,6 +1012,11 @@ end
     return nothing
 end
 
+"""
+    __apply_algebraic_constraint_oop!(residᵢ, algebraic_indices, f, yᵢ₊₁, p, t)
+
+Same as [`__apply_algebraic_constraint!`](@ref) for an out-of-place `f`.
+"""
 @inline function __apply_algebraic_constraint_oop!(
         residᵢ, ::Nothing, f, yᵢ₊₁, p, t
     )
@@ -1000,6 +1031,13 @@ end
     return nothing
 end
 
+"""
+    __check_dae_adaptivity(algebraic_indices, adaptive)
+
+Throw an `ArgumentError` when mesh adaptivity is requested for a DAE problem
+(`algebraic_indices !== nothing`): the collocation interpolant is inaccurate for algebraic
+variables, so defect-based mesh refinement cannot converge.
+"""
 @inline __check_dae_adaptivity(::Nothing, adaptive::Bool) = nothing
 
 @inline function __check_dae_adaptivity(::Vector{Int}, adaptive::Bool)