From 7ede4ee861a6ad862ea5598d229af37ccfa290e0 Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Sun, 14 Jun 2026 14:33:44 -0400 Subject: [PATCH 1/2] fixing bug where if we assembled stiffness and mass the cscnzval which is really the nzval in the CSC matrix gets overwritten. This was likely only an issue for csc/csr matrices on the CPU. Need to add options to not explicitly store this data when the backend is anything but CPU. --- src/AppTools.jl | 101 +++++++++++++----------- src/Parameters.jl | 30 +++++++ src/assemblers/Assemblers.jl | 29 +++++-- src/assemblers/Diagonal.jl | 53 ++++++++++--- src/assemblers/SparseMatrixAssembler.jl | 9 ++- src/assemblers/SparsityPatterns.jl | 77 +++++++++--------- src/assemblers/Vector.jl | 1 - src/bcs/BoundaryConditions.jl | 40 +++++++--- 8 files changed, 220 insertions(+), 120 deletions(-) diff --git a/src/AppTools.jl b/src/AppTools.jl index 7d726cf5..8d34ebc4 100644 --- a/src/AppTools.jl +++ b/src/AppTools.jl @@ -15,6 +15,7 @@ import ..NeumannBC import ..PeriodicBC import ..RobinBC import ..Source +import ..TimeStepper import ..UnstructuredMesh import ..element_blocks import ..element_ids @@ -285,6 +286,7 @@ print_dict(l::LogFile, d::Dict{String, String}) = print_dict(l, d) # Input file ####################################################### struct FunctionSettings{N, T <: Number} + dict::Dict{String, Any} scalar_expr_funcs::Dict{String, ScalarExpressionFunction{T}} vector_expr_funcs::Dict{String, VectorExpressionFunction{N, T}} @@ -321,12 +323,15 @@ struct FunctionSettings{N, T <: Number} @assert false "Unsupported function type $type" end end + else + func_settings = Dict{String, Any}() end - return new{N, T}(scalar_functions, vector_functions) + return new{N, T}(func_settings, scalar_functions, vector_functions) end end struct BCSettings{N, T <: Number} + dict::Dict{String, Any} dirichlet::Vector{DirichletBC{ScalarExpressionFunction{T}}} neumann::Vector{NeumannBC{VectorExpressionFunction{N, T}}} periodic::Vector{PeriodicBC{ScalarExpressionFunction{T}}} @@ -435,17 +440,10 @@ struct BCSettings{N, T <: Number} end end - new{N, T}(dbcs, nbcs, pbcs, rbcs, srcs) + new{N, T}(bc_settings, dbcs, nbcs, pbcs, rbcs, srcs) end end - -# struct FunctionSpaceSettings -# field_type::String -# polynomial_type::String -# end - -# TODO currently only supports blocks struct ICSettings{T <: Number} ics::Vector{InitialCondition{ScalarExpressionFunction{T}}} @@ -489,13 +487,26 @@ struct ICSettings{T <: Number} end end +struct MaterialSettings + materials::Dict{String, Any} + + function MaterialSettings(log_file, parser) + if haskey(parser, "materials") + materials = InputFileParser.get_nested_block(parser, "materials") + else + materials = Dict{String, Any}() + end + new(materials) + end +end + struct MeshSettings dimension::Int file_path::String file_type::String - function MeshSettings(log_file, data) - mesh_settings = data["mesh"]::Dict{String, Any} + function MeshSettings(log_file, parser) + mesh_settings = InputFileParser.get_nested_block(parser, "mesh") dimension = mesh_settings["dimension"]::Int file_path = mesh_settings["file path"]::String file_type = lowercase(mesh_settings["file type"]::String) @@ -503,18 +514,38 @@ struct MeshSettings end end -# struct VariableSettings -# function_space::String -# name::String -# variable_type::String -# end +struct TimeSettings{T <: Number} + dict::Dict{String, Any} + time::Union{Nothing, TimeStepper{T}} + + function TimeSettings(log_file, parser) + if haskey(parser, "time") + time_settings = InputFileParser.get_nested_block(parser, "time") + end_time = time_settings["end time"]::Float64 + start_time = time_settings["start time"]::Float64 + if haskey(time_settings, "number of time steps") + nt = time_settings["number of time steps"]::Int + Δt = (end_time - start_time) / nt + elseif haskey(time_settings, "time step") + Δt = time_settings["time step"]::Float64 + else + @assert false "Requires either one of \"number of time steps\" or \"time step\" commands." + end + new{Float64}(time_settings, TimeStepper(start_time, end_time, start_time, Δt)) + else + new{Float64}(Dict{String, Any}(), nothing) + end + end +end struct InputSettings{N, T <: Number} bcs::BCSettings{N, T} functions::FunctionSettings{N, T} ics::ICSettings{T} + materials::MaterialSettings mesh::MeshSettings parser::InputFileParser.Parser + time::TimeSettings{T} end function InputSettings{N}(cli_args::CLIArgParser, log_file::LogFile, ::Type{T} = Float64) where {N, T <: Number} @@ -532,35 +563,11 @@ function InputSettings{N}(cli_args::CLIArgParser, log_file::LogFile, ::Type{T} = functions = FunctionSettings{N, T}(log_file, parser) bcs = BCSettings{N, T}(log_file, parser, functions) ics = ICSettings{T}(log_file, parser, functions) + materials = MaterialSettings(log_file, parser) mesh = MeshSettings(log_file, parser) - return InputSettings{N, T}(bcs, functions, ics, mesh, parser) -end - -# function _parse_function_space_settings(log_file, data) -# settings = data["function_spaces"]::Dict{String, Any} -# fspaces = Dict{String, FunctionSpaceSettings}() -# for (k, v) in settings -# name = k::String -# temp = v::Dict{String, Any} -# field_type = temp["field_type"]::String -# polynomial_type = temp["polynomial_type"]::String -# fspaces[name] = FunctionSpaceSettings(field_type, polynomial_type) -# end -# return fspaces -# end - -# function _parse_variable_settings(log_file, data) -# var_settings = data["variables"]::Dict{String, Any} -# vars = Dict{String, VariableSettings}() -# for (k, v) in var_settings -# name = k::String -# temp = v::Dict{String, Any} -# fspace = temp["function_space"]::String -# type = temp["type"]::String -# vars[name] = VariableSettings(fspace, name, type) -# end -# return vars -# end + time = TimeSettings(log_file, parser) + return InputSettings{N, T}(bcs, functions, ics, materials, mesh, parser, time) +end ####################################################### # MeshIO strongly typed helpers @@ -656,6 +663,7 @@ end struct Simulation{D, N, T <: Number, IO, Mesh} dbcs::Vector{DirichletBC{ScalarExpressionFunction{T}}} ics::Vector{InitialCondition{ScalarExpressionFunction{T}}} + input_settings::InputSettings{N, T} log_file::LogFile{IO} mesh::Mesh nbcs::Vector{NeumannBC{VectorExpressionFunction{N, T}}} @@ -676,12 +684,13 @@ struct Simulation{D, N, T <: Number, IO, Mesh} # println(log_file.io, fspace::FunctionSpace) # end if length(settings.ics.ics) > 1 - T = eltype(settings.ics[1]) + T = eltype(settings.ics.ics[1]) else T = Float64 end new{D, N, T, IO, typeof(mesh)}( - settings.bcs.dirichlet, settings.ics.ics, log_file, mesh, + settings.bcs.dirichlet, settings.ics.ics, + settings, log_file, mesh, settings.bcs.neumann, settings.bcs.periodic, settings.bcs.robin, settings.bcs.source ) diff --git a/src/Parameters.jl b/src/Parameters.jl index a49ec035..b6c316a9 100644 --- a/src/Parameters.jl +++ b/src/Parameters.jl @@ -253,6 +253,36 @@ struct TypeStableParameters{ physics, props, state_old, state_new, coords, field, field_old, hvp_scratch_field ) end + + function TypeStableParameters{SF, VF}( + mesh, assembler, physics, props, state_old, state_new, + ics, dbcs, nbcs, pbcs, srcs, times + ) where {SF, VF} + dof = assembler.dof + ND = size(dof, 1) + ics = InitialConditions{SF}(mesh, dof, ics) + dbcs = DirichletBCs{SF}(mesh, dof, dbcs) + nbcs = NeumannBCs{VF}(mesh, dof, nbcs) + pbcs = PeriodicBCs{SF}(mesh, dof, pbcs) + srcs = Sources{VF}(mesh, dof, srcs) + + coords = mesh.nodal_coords + field = create_field(assembler) + field_old = create_field(assembler) + hvp_scratch_field = create_field(assembler) + + # update assembler, where should this really live? + update_dofs!(assembler, dbcs, pbcs) + + new{ + SF, VF, Int, Float64, Vector{Int}, Vector{Float64}, Matrix{SVector{ND, Float64}}, + typeof(physics), typeof(props), typeof(mesh.nodal_coords), typeof(field) + }( + ics, dbcs, nbcs, pbcs, srcs, + times, + physics, props, state_old, state_new, coords, field, field_old, hvp_scratch_field + ) + end end function create_parameters( diff --git a/src/assemblers/Assemblers.jl b/src/assemblers/Assemblers.jl index dd009d58..92e0bbfc 100644 --- a/src/assemblers/Assemblers.jl +++ b/src/assemblers/Assemblers.jl @@ -250,14 +250,31 @@ function _quadrature_level_state(state::AbstractArray{<:Number, 3}, q::Int, e::I return state_q end -function _sparse_matrix(asm::AbstractAssembler, storage) +function _sparse_matrix_mass(asm::AbstractAssembler, coo_storage) type = _sparse_matrix_type(asm) if isa(type, COOMatrix) - return _coo_matrix(asm.matrix_pattern, asm.dof, storage) + return _coo_matrix(asm.matrix_pattern, asm.dof, coo_storage) elseif isa(type, CSCMatrix) - return _csc_matrix(asm.matrix_pattern, asm.dof, storage) + csc_storage = asm.matrix_pattern.cscnzval_mass + return _csc_matrix(asm.matrix_pattern, asm.dof, coo_storage, csc_storage) elseif isa(type, CSRMatrix) - return _csr_matrix(asm.matrix_pattern, asm.dof, storage) + csc_storage = asm.matrix_pattern.cscnzval_mass + return _csr_matrix(asm.matrix_pattern, asm.dof, coo_storage, csc_storage) + else + @assert false "Should never happen. Handle this case better. Type is $(_sparse_matrix_type(asm))" + end +end + +function _sparse_matrix_stiffness(asm::AbstractAssembler, coo_storage) + type = _sparse_matrix_type(asm) + if isa(type, COOMatrix) + return _coo_matrix(asm.matrix_pattern, asm.dof, coo_storage) + elseif isa(type, CSCMatrix) + csc_storage = asm.matrix_pattern.cscnzval_stiffness + return _csc_matrix(asm.matrix_pattern, asm.dof, coo_storage, csc_storage) + elseif isa(type, CSRMatrix) + csc_storage = asm.matrix_pattern.cscnzval_stiffness + return _csr_matrix(asm.matrix_pattern, asm.dof, coo_storage, csc_storage) else @assert false "Should never happen. Handle this case better. Type is $(_sparse_matrix_type(asm))" end @@ -314,7 +331,7 @@ function mass(asm::AbstractAssembler) return _zero_sparse_matrix(asm) end - M = _sparse_matrix(asm, asm.mass_storage) + M = _sparse_matrix_mass(asm, asm.mass_storage) if _is_condensed(asm.dof) _adjust_matrix_entries_for_constraints!(M, asm.constraint_storage) @@ -361,7 +378,7 @@ function stiffness(asm::AbstractAssembler) return _zero_sparse_matrix(asm) end - K = _sparse_matrix(asm, asm.stiffness_storage) + K = _sparse_matrix_stiffness(asm, asm.stiffness_storage) if _is_condensed(asm.dof) _adjust_matrix_entries_for_constraints!(K, asm.constraint_storage) diff --git a/src/assemblers/Diagonal.jl b/src/assemblers/Diagonal.jl index 087e687b..91649d2a 100644 --- a/src/assemblers/Diagonal.jl +++ b/src/assemblers/Diagonal.jl @@ -14,11 +14,26 @@ gives the true diagonal, whereas the row-sum approximation `K·1` can be zero at interior nodes of uniform meshes. """ function assemble_diagonal!( - assembler, func::F, Uu, p + assembler::AbstractAssembler, func::F, u, p ) where F <: Function - storage = assembler.residual_storage + assemble_diagonal!( + assembler.residual_storage, + assembler.vector_pattern, assembler.dof, + func, u, p; + use_inplace_methods = _use_inplace_methods(assembler) + ) +end + +# function assemble_diagonal!( +# assembler, func::F, Uu, p +# ) where F <: Function + # storage = assembler.residual_storage +function assemble_diagonal!( + storage, pattern, dof, func, Uu, p; + use_inplace_methods::Bool = false +) fill!(storage, zero(eltype(storage))) - dof = assembler.dof + # dof = assembler.dof fspace = function_space(dof) X = coordinates(p) t = current_time(p) @@ -29,16 +44,28 @@ function assemble_diagonal!( return_type = AssembledDiagonal() conns = fspace.elem_conns foreach_block(fspace, p) do physics, props, ref_fe, b - _assemble_block!( - storage, - conns.data, conns.offsets[b], - func, - physics, ref_fe, - X, t, Δt, - U, U_old, - block_view(p.state_old, b), block_view(p.state_new, b), props, - return_type - ) + if use_inplace_methods + _assemble_block!( + storage, + func, + physics, + t, Δt, + props, + block_view(p.state_old, b), block_view(p.state_new, b), + conns.data, conns.offsets[b], ref_fe, X, U, U_old + ) + else + _assemble_block!( + storage, + conns.data, conns.offsets[b], + func, + physics, ref_fe, + X, t, Δt, + U, U_old, + block_view(p.state_old, b), block_view(p.state_new, b), props, + return_type + ) + end end return nothing end diff --git a/src/assemblers/SparseMatrixAssembler.jl b/src/assemblers/SparseMatrixAssembler.jl index 14ff32a9..3c91336f 100644 --- a/src/assemblers/SparseMatrixAssembler.jl +++ b/src/assemblers/SparseMatrixAssembler.jl @@ -140,15 +140,16 @@ function _empty_matrix_pattern(dof::DofManager) similar(dof.unknown_dofs, 0), # Is similar(dof.unknown_dofs, 0), # Js similar(dof.unknown_dofs, 0), # unknown_dofs - Int[], # block_start_indices - [0], # max_entries (single zero — see create_assembler_cache) + Int[], # block_start_indices + [0], # max_entries (single zero — see create_assembler_cache) similar(dof.unknown_dofs, 0), # klasttouch similar(dof.unknown_dofs, 0), # csrrowptr similar(dof.unknown_dofs, 0), # csrcolval - Float64[], # csrnzval + Float64[], # csrnzval similar(dof.unknown_dofs, 0), # csccolptr similar(dof.unknown_dofs, 0), # cscrowval - Float64[], # cscnzval + Float64[], # cscnzval_mass + Float64[], # cscnzval_stiffness similar(dof.unknown_dofs, 0) # permutation ) end diff --git a/src/assemblers/SparsityPatterns.jl b/src/assemblers/SparsityPatterns.jl index a11b8b98..9425cc59 100644 --- a/src/assemblers/SparsityPatterns.jl +++ b/src/assemblers/SparsityPatterns.jl @@ -44,7 +44,8 @@ struct SparseMatrixPattern{ # additional cache arrays csccolptr::I cscrowval::I - cscnzval::R + cscnzval_mass::R + cscnzval_stiffness::R # trying something new permutation::I end @@ -89,9 +90,10 @@ function SparseMatrixPattern(dof::DofManager) csrcolval = zeros(Int64, length(Is)) csrnzval = zeros(Float64, length(Is)) - csccolptr = Vector{Int64}(undef, 0) - cscrowval = Vector{Int64}(undef, 0) - cscnzval = Vector{Float64}(undef, 0) + csccolptr = Vector{Int64}(undef, 0) + cscrowval = Vector{Int64}(undef, 0) + cscnzval_mass = Vector{Float64}(undef, 0) + cscnzval_stiffness = Vector{Float64}(undef, 0) # set permutation — pack (row, col) into one Int64 key so sortperm takes the # integer radix-sort path instead of an O(N log N) tuple comparison sort. @@ -107,7 +109,7 @@ function SparseMatrixPattern(dof::DofManager) # cache arrays klasttouch, csrrowptr, csrcolval, csrnzval, # additional cache arrays - csccolptr, cscrowval, cscnzval, + csccolptr, cscrowval, cscnzval_mass, cscnzval_stiffness, permutation ) @@ -115,25 +117,19 @@ function SparseMatrixPattern(dof::DofManager) end function Adapt.adapt_structure(to, asm::SparseMatrixPattern) - Is = adapt(to, asm.Is) - Js = adapt(to, asm.Js) - unknown_dofs = adapt(to, asm.unknown_dofs) - # - klasttouch = adapt(to, asm.klasttouch) - csrrowptr = adapt(to, asm.csrrowptr) - csrcolval = adapt(to, asm.csrcolval) - csrnzval = adapt(to, asm.csrnzval) - # - csccolptr = adapt(to, asm.csccolptr) - cscrowval = adapt(to, asm.cscrowval) - cscnzval = adapt(to, asm.cscnzval) - perm = adapt(to, asm.permutation) return SparseMatrixPattern( - Is, Js, - unknown_dofs, + adapt(to, asm.Is), + adapt(to, asm.Js), + adapt(to, asm.unknown_dofs), asm.block_start_indices, asm.max_entries, - klasttouch, csrrowptr, csrcolval, csrnzval, - csccolptr, cscrowval, cscnzval, perm + adapt(to, asm.klasttouch), + adapt(to, asm.csrrowptr), + adapt(to, asm.csrcolval), + adapt(to, asm.csrnzval), + adapt(to, asm.csccolptr), + adapt(to, asm.cscrowval), + adapt(to, asm.cscnzval_mass), adapt(to, asm.cscnzval_stiffness), + adapt(to, asm.permutation) ) end @@ -141,12 +137,12 @@ function KA.get_backend(pattern::SparseMatrixPattern) return KA.get_backend(pattern.Is) end -function SparseArrays.sparse!(pattern::SparseMatrixPattern, storage) +function _sparse_cpu!(pattern::SparseMatrixPattern, coo_storage, csc_storage) return @views SparseArrays.sparse!( - pattern.Is, pattern.Js, storage[pattern.unknown_dofs], + pattern.Is, pattern.Js, coo_storage[pattern.unknown_dofs], length(pattern.klasttouch), length(pattern.klasttouch), +, pattern.klasttouch, pattern.csrrowptr, pattern.csrcolval, pattern.csrnzval, - pattern.csccolptr, pattern.cscrowval, pattern.cscnzval + pattern.csccolptr, pattern.cscrowval, csc_storage ) end @@ -228,8 +224,6 @@ function _update_dofs!(pattern::SparseMatrixPattern, dof, dirichlet_dofs, period resize!(pattern.klasttouch, n_total_dofs) resize!(pattern.csrrowptr, n_total_dofs + 1) # TODO Not sure about below 2 sizes - # resize!(assembler.pattern.csrcolval, length(assembler.pattern.Is)) - # resize!(assembler.pattern.csrnzval, length(assembler.pattern.Is)) resize!(pattern.csrcolval, length(pattern.Is)) resize!(pattern.csrnzval, length(pattern.Is)) @@ -302,39 +296,44 @@ function _coo_matrix_constructor(backend::KA.Backend) @assert false "Need to implement for $backend" end -function _csc_matrix(pattern::SparseMatrixPattern, dof, storage) +function _csc_matrix(pattern::SparseMatrixPattern, dof, coo_storage, csc_storage) backend = KA.get_backend(pattern) - return _csc_matrix(backend, pattern, dof, storage) + return _csc_matrix(backend, pattern, dof, coo_storage, csc_storage) end -function _csc_matrix(backend::KA.Backend, pattern::SparseMatrixPattern, dof, storage) - coo = _coo_matrix(backend, pattern, dof, storage) +function _csc_matrix(backend::KA.Backend, pattern::SparseMatrixPattern, dof, coo_storage, csc_storage) + coo = _coo_matrix(backend, pattern, dof, coo_storage) constructor = _csc_matrix_constructor(backend) return constructor(coo) end -function _csc_matrix(::KA.CPU, pattern::SparseMatrixPattern, dof, storage) - return SparseArrays.sparse!(pattern, storage) +function _csc_matrix(::KA.CPU, pattern::SparseMatrixPattern, dof, coo_storage, csc_storage) + return @views SparseArrays.sparse!( + pattern.Is, pattern.Js, coo_storage[pattern.unknown_dofs], + length(pattern.klasttouch), length(pattern.klasttouch), +, pattern.klasttouch, + pattern.csrrowptr, pattern.csrcolval, pattern.csrnzval, + pattern.csccolptr, pattern.cscrowval, csc_storage + ) end function _csc_matrix_constructor(backend::KA.Backend) @assert false "Need to implement for $backend" end -function _csr_matrix(pattern::SparseMatrixPattern, dof, storage) +function _csr_matrix(pattern::SparseMatrixPattern, dof, coo_storage, csc_storage) backend = KA.get_backend(pattern) - return _csr_matrix(backend, pattern, dof, storage) + return _csr_matrix(backend, pattern, dof, coo_storage, csc_storage) end -function _csr_matrix(backend::KA.Backend, pattern::SparseMatrixPattern, dof, storage) - coo = _coo_matrix(backend, pattern, dof, storage) +function _csr_matrix(backend::KA.Backend, pattern::SparseMatrixPattern, dof, coo_storage, csc_storage) + coo = _coo_matrix(backend, pattern, dof, coo_storage) constructor = _csr_matrix_constructor(backend) return constructor(coo) end # TODO eventually write kernels to do this without going to csc first -function _csr_matrix(backend::KA.CPU, pattern::SparseMatrixPattern, dof, storage) - csc = _csc_matrix(backend, pattern, dof, storage) +function _csr_matrix(backend::KA.CPU, pattern::SparseMatrixPattern, dof, coo_storage, csc_storage) + csc = _csc_matrix(backend, pattern, dof, coo_storage, csc_storage) return SparseMatrixCSR(csc) end diff --git a/src/assemblers/Vector.jl b/src/assemblers/Vector.jl index fb0fc682..1913dee4 100644 --- a/src/assemblers/Vector.jl +++ b/src/assemblers/Vector.jl @@ -10,7 +10,6 @@ function assemble_vector!( storage = assembler.residual_storage end assemble_vector!( - # assembler.residual_storage, storage, assembler.vector_pattern, assembler.dof, func, Uu, p; diff --git a/src/bcs/BoundaryConditions.jl b/src/bcs/BoundaryConditions.jl index 6f356018..c9d9726c 100644 --- a/src/bcs/BoundaryConditions.jl +++ b/src/bcs/BoundaryConditions.jl @@ -16,19 +16,37 @@ struct VariableNameNotFoundError <: Exception end _var_not_found_err() = throw(VariableNameNotFoundError()) +# function _dof_index_from_var_name(dof, var_name) +# # dbc/nbc or robin for scalar case +# dof_index = findfirst(x -> x == var_name, names(dof.var)) +# if dof_index === nothing +# # try a vector/tensor neumann/robin bc type +# dof_index = findfirst(x -> occursin("$(var_name)_", x), names(dof.var)) + +# # if we still haven't found anything, this variable likely doesn't exist +# if dof_index === nothing +# _var_not_found_err() +# end +# end +# return dof_index +# end function _dof_index_from_var_name(dof, var_name) - # dbc/nbc or robin for scalar case - dof_index = findfirst(x -> x == var_name, names(dof.var)) - if dof_index === nothing - # try a vector/tensor neumann/robin bc type - dof_index = findfirst(x -> occursin("$(var_name)_", x), names(dof.var)) - - # if we still haven't found anything, this variable likely doesn't exist - if dof_index === nothing - _var_not_found_err() - end + var_names = names(dof.var) + + for i = 1:length(var_names) + if var_names[i] == var_name + return i + end end - return dof_index + + prefix = var_name * "_" + for i = 1:length(var_names) + if startswith(var_names[i], prefix) + return i + end + end + + _var_not_found_err() end function _unique_sort_perm(array::AbstractArray{T, 1}) where T <: Number From cc1fa14637d2c64cb3e04fa2a4693106ba313eae Mon Sep 17 00:00:00 2001 From: "Craig M. Hamel" Date: Sun, 14 Jun 2026 15:58:47 -0400 Subject: [PATCH 2/2] removing dud function. --- src/assemblers/SparsityPatterns.jl | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/assemblers/SparsityPatterns.jl b/src/assemblers/SparsityPatterns.jl index 9425cc59..565eecb9 100644 --- a/src/assemblers/SparsityPatterns.jl +++ b/src/assemblers/SparsityPatterns.jl @@ -137,15 +137,6 @@ function KA.get_backend(pattern::SparseMatrixPattern) return KA.get_backend(pattern.Is) end -function _sparse_cpu!(pattern::SparseMatrixPattern, coo_storage, csc_storage) - return @views SparseArrays.sparse!( - pattern.Is, pattern.Js, coo_storage[pattern.unknown_dofs], - length(pattern.klasttouch), length(pattern.klasttouch), +, pattern.klasttouch, - pattern.csrrowptr, pattern.csrcolval, pattern.csrnzval, - pattern.csccolptr, pattern.cscrowval, csc_storage - ) -end - function block_view(storage::AbstractVector, pattern::SparseMatrixPattern, b::Int) @assert b > 0 && b <= length(pattern.block_start_indices) if b == length(pattern.block_start_indices) || length(pattern.block_start_indices) == 1