diff --git a/src/assemblers/matrix_free/eigensolve.jl b/src/assemblers/matrix_free/eigensolve.jl index f373eb7..edde69c 100644 --- a/src/assemblers/matrix_free/eigensolve.jl +++ b/src/assemblers/matrix_free/eigensolve.jl @@ -1,5 +1,5 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +# SPDX-FileCopyrightText: 2015-2026 Jukka Aho +# SPDX-License-Identifier: MIT import IterativeSolvers import LinearOperators @@ -226,18 +226,22 @@ end # --------------------------------------------------------------------------- """ - solve_eigenproblem(cache, asm, kernel, mesh; + solve_eigenproblem(cache, asm, mesh; nev = 1, tol = 1e-8, maxiter = 200, p = nothing, verbose = false, dirichlet = nothing, mpc = nothing, shift = 0.0) -> (λ, V) + solve_eigenproblem(cache, asm, kernel, mesh; …) Convenience wrapper around `lowest_eigenpairs`: assembles matrix-free `K` and `M` operators (via `apply_K!` / `apply_M!`) and runs subspace iteration to extract the lowest `nev` generalized eigenpairs of -`K φ = λ M φ`. +`K φ = λ M φ`. Volume kernels are read from `cache.kernel_column`. -`dirichlet` and `mpc` are forwarded to `matrix_free_op` so the +The four-argument form ignores `kernel` (backward compatibility; emits +`Base.depwarn` once per session, same as the matrix-free operator overloads). + +`dirichlet` and `mpc` are forwarded to [`MatrixFreeOperator`](@ref) so the constrained operator is solved directly. `shift` adds `σ M` to `K` internally and subtracts `σ` from the returned eigenvalues — useful for problems with rigid-body / null-space modes (free-free elasticity, @@ -245,22 +249,23 @@ unconstrained heat) where the unshifted `K` is singular and the inner CG cannot invert it. A shift slightly larger than the smallest non-trivial eigenvalue is sufficient. """ -function solve_eigenproblem(cache::DOFBasedCOOCache, - asm::DOFBasedCOOAssembler, - kernel::AbstractKernel, - mesh::AbstractMesh; - nev::Int = 1, - tol::Real = 1e-8, - maxiter::Int = 200, - p::Union{Nothing,Int} = nothing, - verbose::Bool = false, - dirichlet = nothing, - mpc = nothing, - shift::Real = 0.0) +function solve_eigenproblem( + cache::DOFBasedCOOCache, + asm::DOFBasedCOOAssembler, + mesh::AbstractMesh; + nev::Int = 1, + tol::Real = 1e-8, + maxiter::Int = 200, + p::Union{Nothing,Int} = nothing, + verbose::Bool = false, + dirichlet = nothing, + mpc = nothing, + shift::Real = 0.0, +) n = cache.ndofs - op_K_base = MatrixFreeOperator(cache, asm, kernel, mesh; + op_K_base = MatrixFreeOperator(cache, asm, mesh; dirichlet = dirichlet, mpc = mpc) - op_M = MatrixFreeMassOperator(cache, asm, kernel, mesh) + op_M = MatrixFreeMassOperator(cache, asm, mesh) # Optional shift: K_shift = K + σ M ⇒ λ_shift = λ + σ. op_K = if shift == 0.0 @@ -286,3 +291,14 @@ function solve_eigenproblem(cache::DOFBasedCOOCache, end return (λ, V) end + +@inline function solve_eigenproblem( + cache::DOFBasedCOOCache, + asm::DOFBasedCOOAssembler, + ::AbstractKernel, + mesh::AbstractMesh; + kwargs..., +) + _depwarn_redundant_kernel_arg!(:solve_eigenproblem) + return solve_eigenproblem(cache, asm, mesh; kwargs...) +end diff --git a/src/assemblers/matrix_free/preconditioners.jl b/src/assemblers/matrix_free/preconditioners.jl index ba46e65..7006c5d 100644 --- a/src/assemblers/matrix_free/preconditioners.jl +++ b/src/assemblers/matrix_free/preconditioners.jl @@ -1,5 +1,5 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +# SPDX-FileCopyrightText: 2015-2026 Jukka Aho +# SPDX-License-Identifier: MIT """ Lightweight, matrix-free preconditioners for Krylov solves built on top @@ -32,8 +32,8 @@ the `A` block when that slice is SPD enough for an incomplete factorisation. using IterativeSolvers, LinearOperators, JuliaFEM bc = PenaltyDirichlet(fixed_dofs, vals; penalty = 1e8) -op = matrix_free_op(cache, asm, kernel, mesh; dirichlet = bc) -P = JacobiPreconditioner(cache, asm, kernel, mesh; dirichlet = bc) +op = matrix_free_op(cache, asm, mesh; dirichlet = bc) +P = JacobiPreconditioner(cache, asm, mesh; dirichlet = bc) linop = LinearOperator(Float64, n, n, true, true, op) @@ -57,15 +57,14 @@ using SparseArrays: SparseMatrixCSC, sparse, nnz, rowvals, nonzeros, getcolptr Diagonal preconditioner with stored *inverse* diagonal so that `ldiv!(P, x)` is a single multiplication per entry. Direct constructor; -prefer the `(cache, asm, kernel, mesh; dirichlet)` factory for the -common matrix-free path. +prefer the `(cache, asm, mesh; dirichlet)` factory for the common matrix-free path. """ struct JacobiPreconditioner inv_diag::Vector{Float64} end """ - JacobiPreconditioner(cache, asm, kernel, mesh; dirichlet = nothing, mpc = nothing) + JacobiPreconditioner(cache, asm, mesh; dirichlet = nothing, mpc = nothing) Build a Jacobi preconditioner from the matrix-free operator. Internally: @@ -78,16 +77,21 @@ Build a Jacobi preconditioner from the matrix-free operator. Internally: 3. invert entry-wise (zeros are kept as `1.0` so the preconditioner stays well-defined; in practice every DOF row of `K` has a positive diagonal in any well-posed FEM problem). + + JacobiPreconditioner(cache, asm, kernel, mesh; …) + +Four-argument backward-compatibility overload: `kernel` is ignored; volume +physics come from `cache.kernel_column` only (same as [`compute_diagonal!`](@ref)). +Emits `Base.depwarn` once per session (shared with other redundant-kernel overloads). """ function JacobiPreconditioner(cache::DOFBasedCOOCache, asm::DOFBasedCOOAssembler, - kernel::AbstractKernel, mesh::AbstractMesh; dirichlet::Union{AbstractDirichletConstraint, Nothing} = nothing, mpc::Union{AbstractMultipointConstraint, Nothing} = nothing) n = cache.ndofs d = zeros(Float64, n) - compute_diagonal!(d, cache, asm, kernel, mesh) + compute_diagonal!(d, cache, asm, mesh) if dirichlet !== nothing apply_constraint_diag!(d, dirichlet) end @@ -102,6 +106,15 @@ function JacobiPreconditioner(cache::DOFBasedCOOCache, return JacobiPreconditioner(inv_d) end +@inline function JacobiPreconditioner(cache::DOFBasedCOOCache, + asm::DOFBasedCOOAssembler, + ::AbstractKernel, + mesh::AbstractMesh; + kwargs...) + _depwarn_redundant_kernel_arg!(:JacobiPreconditioner) + return JacobiPreconditioner(cache, asm, mesh; kwargs...) +end + # `ldiv!(y, P, x)` is the contract IterativeSolvers.cg! / Krylov.jl call # for left-preconditioning: compute `y = P^{-1} * x`. For a diagonal # `P`, this is just an entry-wise multiplication by the stored inverse. @@ -143,7 +156,6 @@ end compute_block_diagonal!(blocks::Array{Float64,3}, cache::DOFBasedCOOCache, asm::DOFBasedCOOAssembler, - kernel::AbstractKernel, mesh::AbstractMesh) -> blocks Assemble the `N × N` block-diagonal of the stiffness matrix into @@ -153,14 +165,18 @@ for vector fields). Same DOF-by-element traversal as `assemble!` but only the entries with both DOFs in the same block contribute. Allocation-free after warmup. + + compute_block_diagonal!(blocks, cache, asm, kernel, mesh) + +Five-argument overload: `kernel` is ignored; volume kernels are read from +`cache.kernel_column`. Emits `Base.depwarn` once per session when used. """ function compute_block_diagonal!( blocks::Array{Float64,3}, - cache::DOFBasedCOOCache{T,B,IPS,E,GC,Buf,FieldType,StateType}, + cache::DOFBasedCOOCache{T,B,IPS,E,GC,Buf,FieldType,StateType,KS}, asm::DOFBasedCOOAssembler, - kernel::AbstractKernel, mesh::AbstractMesh, -) where {T,B,IPS,E<:AbstractElement,GC,Buf,FieldType,StateType} +) where {T,B,IPS,E<:AbstractElement,GC,Buf,FieldType,StateType,KS} N, N2, n_blocks = size(blocks) @assert N == N2 ( "compute_block_diagonal!: blocks must be (N,N,n_blocks); got $(size(blocks))") @@ -169,7 +185,7 @@ function compute_block_diagonal!( "cache.ndofs = $(cache.ndofs)") fill!(blocks, 0.0) - _prepare_caches!(cache, kernel, mesh) + _prepare_caches!(cache, mesh) elements = cache.elements element_caches = cache.element_caches @@ -180,6 +196,7 @@ function compute_block_diagonal!( ndofs_elem = length(layout) @inbounds for elem_idx in 1:length(elements) + k_e = kernel_at(cache, elem_idx) ec = element_caches[elem_idx] gc = geometry_caches[elem_idx] qp = view(qp_buffers, :, elem_idx) @@ -200,7 +217,7 @@ function compute_block_diagonal!( comp_j = mod(dof_j - 1, N) + 1 entry_j = layout[lj] - K_ij = evaluate_entry(kernel, gc, qp, entry_i, entry_j, elem_idx) + K_ij = evaluate_entry(k_e, gc, qp, entry_i, entry_j, elem_idx) blocks[comp_i, comp_j, blk_i] += K_ij end end @@ -209,6 +226,17 @@ function compute_block_diagonal!( return blocks end +@inline function compute_block_diagonal!( + blocks::Array{Float64,3}, + cache::DOFBasedCOOCache{T,B,IPS,E,GC,Buf,FieldType,StateType,KS}, + asm::DOFBasedCOOAssembler, + ::AbstractKernel, + mesh::AbstractMesh, +) where {T,B,IPS,E<:AbstractElement,GC,Buf,FieldType,StateType,KS} + _depwarn_redundant_kernel_arg!(:compute_block_diagonal!) + return compute_block_diagonal!(blocks, cache, asm, mesh) +end + """ apply_constraint_block_diag!(blocks::Array{Float64,3}, c) -> blocks @@ -276,7 +304,7 @@ struct BlockJacobiPreconditioner{N} end """ - BlockJacobiPreconditioner{N}(cache, asm, kernel, mesh; dirichlet = nothing) + BlockJacobiPreconditioner{N}(cache, asm, mesh; dirichlet = nothing) Build a block-Jacobi preconditioner from the matrix-free operator. Internally: @@ -288,11 +316,14 @@ Internally: singular blocks, fall back to identity to keep the preconditioner well-defined (a singular nodal block typically signals an ill-posed problem, so this is intentionally permissive). + + BlockJacobiPreconditioner{N}(cache, asm, kernel, mesh; …) + +Five-argument overload: `kernel` is ignored (emits `Base.depwarn` once per session). """ function BlockJacobiPreconditioner{N}( cache::DOFBasedCOOCache, asm::DOFBasedCOOAssembler, - kernel::AbstractKernel, mesh::AbstractMesh; dirichlet::Union{AbstractDirichletConstraint, Nothing} = nothing, ) where {N} @@ -302,7 +333,7 @@ function BlockJacobiPreconditioner{N}( "BlockJacobiPreconditioner{N=$N}: ndofs = $(cache.ndofs) is not divisible by N") blocks = zeros(Float64, N, N, n_blocks) - compute_block_diagonal!(blocks, cache, asm, kernel, mesh) + compute_block_diagonal!(blocks, cache, asm, mesh) if dirichlet !== nothing apply_constraint_block_diag!(blocks, dirichlet) end @@ -326,6 +357,17 @@ function BlockJacobiPreconditioner{N}( return BlockJacobiPreconditioner{N}(inv_blocks) end +@inline function BlockJacobiPreconditioner{N}( + cache::DOFBasedCOOCache, + asm::DOFBasedCOOAssembler, + ::AbstractKernel, + mesh::AbstractMesh; + kwargs..., +) where {N} + _depwarn_redundant_kernel_arg!(:BlockJacobiPreconditioner) + return BlockJacobiPreconditioner{N}(cache, asm, mesh; kwargs...) +end + # ldiv!(y, P, x): y = P^{-1} x, block by block. function ldiv!(y::AbstractVector{Float64}, P::BlockJacobiPreconditioner{N}, @@ -504,7 +546,7 @@ end Lower-triangular Incomplete-Cholesky preconditioner. Build it with either `ICholPreconditioner(K::SparseMatrixCSC)` (factor an existing -SPD matrix) or `ICholPreconditioner(cache, asm, kernel, mesh; dirichlet)` +SPD matrix) or `ICholPreconditioner(cache, asm, mesh; dirichlet)` (assemble + factor in one shot for the matrix-free workflow). `ldiv!(P, x)` solves `L L^T y = x` via two sparse triangular solves @@ -583,25 +625,28 @@ function ICholPreconditioner(K::SparseMatrixCSC{Float64,Int}) end """ - ICholPreconditioner(cache::DOFBasedCOOCache, asm, kernel, mesh; + ICholPreconditioner(cache::DOFBasedCOOCache, asm, mesh; dirichlet = nothing) -> ICholPreconditioner Convenience constructor for the matrix-free workflow: - 1. assemble K via `assemble!(cache, asm, kernel, mesh)` and - extract it with `extract_system(cache)`, + 1. assemble K via `assemble!(cache, asm, mesh)` and extract it with + `extract_system(cache)`, 2. apply the Dirichlet constraint (if given) to `K`, 3. build `IC(0)` of the modified `K`. This *does* materialise the sparse `K` once — the price of IC(0). The returned preconditioner then plugs into the matrix-free Krylov solve. + + ICholPreconditioner(cache, asm, kernel, mesh; …) + +Four-argument overload: `kernel` is ignored (same as [`assemble!`](@ref)). """ function ICholPreconditioner(cache::DOFBasedCOOCache, asm::DOFBasedCOOAssembler, - kernel::AbstractKernel, mesh::AbstractMesh; dirichlet::Union{AbstractDirichletConstraint, Nothing} = nothing) - assemble!(cache, asm, kernel, mesh) + assemble!(cache, asm, mesh) K, _ = extract_system(cache) if dirichlet !== nothing apply_constraint!(K, dirichlet) @@ -609,6 +654,15 @@ function ICholPreconditioner(cache::DOFBasedCOOCache, return ICholPreconditioner(K) end +@inline function ICholPreconditioner(cache::DOFBasedCOOCache, + asm::DOFBasedCOOAssembler, + ::AbstractKernel, + mesh::AbstractMesh; + kwargs...) + _depwarn_redundant_kernel_arg!(:ICholPreconditioner) + return ICholPreconditioner(cache, asm, mesh; kwargs...) +end + # Crout left-looking IC(0) on CSC storage. Returns the factor `L` # such that `L * L^T ≈ K` and `nnz(L) == nnz(tril(K))`. function _ic0_factor(K::SparseMatrixCSC{Float64,Int})