Skip to content

Commit 191c161

Browse files
CWillbergCWillberg
committed
printing works
1 parent 89faca3 commit 191c161

File tree

2 files changed

+692
-121
lines changed

2 files changed

+692
-121
lines changed

‎src/Core/Solver/Matrix_linear_static.jl‎

Lines changed: 232 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,27 @@ using ..Model_Factory.Pre_Calculation.Bond_Deformation
3232
export init_solver
3333
export run_solver
3434

35+
mutable struct DisplacementSolverCache
36+
# Workspace
37+
F_modified::Vector{Float64}
38+
bc_mask::BitVector
39+
temp::Vector{Float64}
40+
u_free::Vector{Float64}
41+
42+
# LU cache
43+
K_free_lu::Union{Nothing,Any}
44+
last_non_BCs::Vector{Int}
45+
46+
# Active DOFs cache
47+
K_active_cached::Union{Nothing,AbstractMatrix{Float64}}
48+
last_active_dofs::Vector{Int}
49+
50+
function DisplacementSolverCache(n_total::Int)
51+
new(Float64[], BitVector(undef, n_total), Float64[], Float64[],
52+
nothing, Int[], nothing, Int[])
53+
end
54+
end
55+
3556
"""
3657
compute_thermodynamic_critical_time_step(nodes::AbstractVector{Int64}, lambda::Float64, Cv::Float64)
3758
@@ -124,7 +145,8 @@ function init_solver(solver_options::Dict{Any,Any},
124145
model_param = Data_Manager.get_properties(block, "Material Model")
125146
Correspondence_matrix_based.init_model(nodes, model_param, block)
126147
end
127-
Correspondence_matrix_based.init_matrix()
148+
149+
@timeit "init matrix" Correspondence_matrix_based.init_matrix()
128150

129151
solver_options["Matrix update"] = get(params["Linear Static Matrix Based"],
130152
"Matrix update", false)
@@ -206,9 +228,7 @@ function run_solver(solver_options::Dict{Any,Any},
206228
dof = Data_Manager.get_dof()
207229
nnodes = Data_Manager.get_nnodes()
208230

209-
volume = Data_Manager.get_field("Volume")
210231
delta_u = Data_Manager.get_field("Delta Displacements")
211-
external_force_densities = Data_Manager.get_field("External Force Densities")
212232

213233
matrix_update::Bool = solver_options["Matrix update"]
214234

@@ -220,7 +240,7 @@ function run_solver(solver_options::Dict{Any,Any},
220240

221241
displacement_solver_cache = DisplacementSolverCache(nnodes * dof)
222242

223-
# Bestimme ob alle Knoten aktiv sind
243+
# Determine if all nodes are active
224244
all_nodes_active = length(active_nodes) == nnodes
225245

226246
non_BCs_global = Data_Manager.get_bc_free_dof()
@@ -267,7 +287,7 @@ function run_solver(solver_options::Dict{Any,Any},
267287
active_list = Data_Manager.get_field("Active")
268288
active_nodes = find_active_nodes(active_list, active_nodes, nodes)
269289

270-
# test status
290+
# Test current status
271291
current_all_active = length(active_nodes) == nnodes
272292

273293
if !current_all_active
@@ -276,26 +296,69 @@ function run_solver(solver_options::Dict{Any,Any},
276296
K = Data_Manager.get_stiffness_matrix()
277297
end
278298

279-
# Lokale Indizierung für reduzierte Matrix
280-
active_dofs_local,
281-
local_to_global = get_active_dof(active_nodes, dof,
282-
nnodes)
283-
active_non_BCs_local = filter_and_map_bc(non_BCs_global, local_to_global)
299+
# Build global DOF indices for active nodes
300+
# Using block-style indexing: global_idx = (d - 1) * nnodes + node
301+
active_dofs_global = Vector{Int}(undef, length(active_nodes) * dof)
302+
idx = 1
303+
for d in 1:dof
304+
for node in active_nodes
305+
active_dofs_global[idx] = (d - 1) * nnodes + node
306+
idx += 1
307+
end
308+
end
309+
sort!(active_dofs_global)
310+
311+
# Filter non_BCs to active DOFs
312+
active_non_BCs_global = intersect(non_BCs_global, active_dofs_global)
284313

285314
@timeit "compute_displacements (partial)" begin
286-
@views compute_displacements!(K,
287-
active_dofs_local,
288-
active_non_BCs_local,
289-
uNP1[active_nodes, :],
290-
force_densities_NP1[active_nodes, :],
291-
external_force_densities[active_nodes, :],
292-
displacement_solver_cache)
293-
for iID in active_nodes
294-
@. @views forces[iID, :] = force_densities_NP1[iID, :] * volume[iID]
315+
if length(active_non_BCs_global) > 0
316+
# Create views to active node data in matrix form
317+
u_active_nodes = @view uNP1[active_nodes, :]
318+
F_int_active_nodes = @view force_densities_NP1[active_nodes, :]
319+
F_ext_active_nodes = @view external_force_densities[active_nodes, :]
320+
321+
# Build mapping from global to active indices
322+
global_to_active = Dict{Int,Int}()
323+
sizehint!(global_to_active, length(active_dofs_global))
324+
for (active_idx, global_idx) in enumerate(active_dofs_global)
325+
global_to_active[global_idx] = active_idx
326+
end
327+
328+
# Map BC indices to active system
329+
active_non_BCs_local = Vector{Int}(undef,
330+
length(active_non_BCs_global))
331+
@inbounds for (i, g) in enumerate(active_non_BCs_global)
332+
active_non_BCs_local[i] = global_to_active[g]
333+
end
334+
335+
# Extract submatrix for active DOFs
336+
K_active = K[active_dofs_global, active_dofs_global]
337+
338+
# Solve system using same structure as reference
339+
# but with mapping for active nodes
340+
compute_displacements_active_subset!(K_active,
341+
active_dofs_global,
342+
active_non_BCs_local,
343+
u_active_nodes,
344+
F_int_active_nodes,
345+
F_ext_active_nodes,
346+
volume[active_nodes],
347+
active_nodes,
348+
nnodes,
349+
dof,
350+
displacement_solver_cache)
351+
end
352+
353+
# Update forces for active nodes
354+
@inbounds for iID in active_nodes
355+
for d in 1:dof
356+
forces[iID, d] = force_densities_NP1[iID, d] * volume[iID]
357+
end
295358
end
296359
end
297360
else
298-
# all nodes are active
361+
# All nodes are active
299362
if matrix_update
300363
@timeit "update matrix (all active)" begin
301364
compute_matrix(active_nodes)
@@ -313,6 +376,13 @@ function run_solver(solver_options::Dict{Any,Any},
313376
force_densities_NP1,
314377
external_force_densities,
315378
displacement_solver_cache)
379+
380+
# Update forces for all nodes
381+
@inbounds for iID in nodes
382+
for d in 1:dof
383+
forces[iID, d] = force_densities_NP1[iID, d] * volume[iID]
384+
end
385+
end
316386
end
317387
end
318388

@@ -324,9 +394,7 @@ function run_solver(solver_options::Dict{Any,Any},
324394

325395
compute_parabolic_problems_after_model_evaluation(active_nodes, solver_options,
326396
dt)
327-
for iID in nodes
328-
@. @views forces[iID, :] = force_densities_NP1[iID, :] * volume[iID]
329-
end
397+
330398
@. velocities = (uNP1 - uN) / dt
331399
@. @views deformed_coorNP1 = coor + uNP1
332400

@@ -345,7 +413,7 @@ function run_solver(solver_options::Dict{Any,Any},
345413
Data_Manager.set_current_time(time)
346414

347415
if idt % ceil(nsteps / 100) == 0
348-
@info "Step: $idt / $(nsteps+1) [$time s]"
416+
@info "Step: $idt / $(nsteps) [$time s]"
349417
end
350418
if rank == 0 && !silent
351419
set_postfix(iter, t = @sprintf("%.4e", time))
@@ -358,6 +426,145 @@ function run_solver(solver_options::Dict{Any,Any},
358426
return result_files
359427
end
360428

429+
# New function for active subset that mimics reference implementation
430+
function compute_displacements_active_subset!(K_active::AbstractMatrix{Float64},
431+
active_dofs_global::AbstractVector{Int},
432+
active_non_BCs::AbstractVector{Int},
433+
u_active::AbstractMatrix{Float64},
434+
F_int_active::AbstractMatrix{Float64},
435+
F_ext_active::AbstractMatrix{Float64},
436+
volume_active::AbstractVector{Float64},
437+
active_nodes::AbstractVector{Int},
438+
nnodes::Int,
439+
dof::Int,
440+
solver::DisplacementSolverCache)
441+
isempty(active_non_BCs) && return nothing
442+
443+
n_active_nodes = length(active_nodes)
444+
n_active_dofs = n_active_nodes * dof
445+
446+
# Vectorize matrices (same as reference)
447+
u_vec = vec(u_active)
448+
F_int_vec = vec(F_int_active)
449+
F_ext_vec = vec(F_ext_active)
450+
451+
n_free = length(active_non_BCs)
452+
has_BCs = n_free < n_active_dofs
453+
454+
# Resize buffers
455+
length(solver.F_modified) != n_free && resize!(solver.F_modified, n_free)
456+
length(solver.bc_mask) != n_active_dofs && resize!(solver.bc_mask, n_active_dofs)
457+
458+
# Modify force vector to account for prescribed displacements (same as reference)
459+
if has_BCs
460+
fill!(solver.bc_mask, true)
461+
@inbounds for i in active_non_BCs
462+
solver.bc_mask[i] = false
463+
end
464+
465+
K_sub = K_active[active_non_BCs, solver.bc_mask]
466+
u_bc = @view u_vec[solver.bc_mask]
467+
468+
length(solver.temp) != n_free && resize!(solver.temp, n_free)
469+
mul!(solver.temp, K_sub, u_bc)
470+
471+
@inbounds for (idx, i) in enumerate(active_non_BCs)
472+
F_int_vec[i] += solver.temp[idx]
473+
end
474+
end
475+
476+
# Build modified force vector (same as reference)
477+
@inbounds for (idx, i) in enumerate(active_non_BCs)
478+
solver.F_modified[idx] = F_ext_vec[i] - F_int_vec[i]
479+
end
480+
481+
# LU factorization and solve (same as reference)
482+
@timeit "LU factorization" begin
483+
K_free = K_active[active_non_BCs, active_non_BCs]
484+
try
485+
solver.K_free_lu = lu(K_free)
486+
catch e
487+
@error "LU factorization failed for reduced system"
488+
@error "Matrix size: $(size(K_free))"
489+
@error "Rank: $(rank(K_free))"
490+
rethrow(e)
491+
end
492+
end
493+
494+
length(solver.u_free) != n_free && resize!(solver.u_free, n_free)
495+
496+
@timeit "ldiv" ldiv!(solver.u_free, solver.K_free_lu, solver.F_modified)
497+
498+
# Write solution back (same as reference)
499+
@inbounds for (idx, i) in enumerate(active_non_BCs)
500+
u_vec[i] = solver.u_free[idx]
501+
end
502+
503+
return nothing
504+
end
505+
506+
# Neue Hilfsfunktion für reduziertes System
507+
function compute_displacements_reduced!(K_active::AbstractMatrix{Float64},
508+
active_non_BCS::AbstractVector{Int},
509+
u_active::AbstractVector{Float64},
510+
F_int_active::AbstractVector{Float64},
511+
F_ext_active::AbstractVector{Float64},
512+
solver::DisplacementSolverCache)
513+
isempty(active_non_BCS) && return nothing
514+
515+
n_total = length(u_active)
516+
n_free = length(active_non_BCS)
517+
518+
length(solver.F_modified) != n_free && resize!(solver.F_modified, n_free)
519+
length(solver.bc_mask) != n_total && resize!(solver.bc_mask, n_total)
520+
521+
has_BCs = n_free < n_total
522+
523+
if has_BCs
524+
fill!(solver.bc_mask, true)
525+
@inbounds for i in active_non_BCS
526+
solver.bc_mask[i] = false
527+
end
528+
529+
K_sub = K_active[active_non_BCS, solver.bc_mask]
530+
u_bc = @view u_active[solver.bc_mask]
531+
532+
length(solver.temp) != n_free && resize!(solver.temp, n_free)
533+
mul!(solver.temp, K_sub, u_bc)
534+
535+
@inbounds for (idx, i) in enumerate(active_non_BCS)
536+
F_int_active[i] += solver.temp[idx]
537+
end
538+
end
539+
540+
@inbounds for (idx, i) in enumerate(active_non_BCS)
541+
solver.F_modified[idx] = F_ext_active[i] - F_int_active[i]
542+
end
543+
544+
# LU factorization
545+
@timeit "LU factorization" begin
546+
K_free = K_active[active_non_BCS, active_non_BCS]
547+
try
548+
solver.K_free_lu = lu(K_free)
549+
catch e
550+
@error "LU factorization failed: $e"
551+
@error "Matrix size: $(size(K_free))"
552+
@error "Condition number: $(cond(K_free))"
553+
rethrow(e)
554+
end
555+
end
556+
557+
length(solver.u_free) != n_free && resize!(solver.u_free, n_free)
558+
559+
@timeit "ldiv" ldiv!(solver.u_free, solver.K_free_lu, solver.F_modified)
560+
561+
@inbounds for (idx, i) in enumerate(active_non_BCS)
562+
u_active[i] = solver.u_free[idx]
563+
end
564+
565+
return nothing
566+
end
567+
361568
function get_active_dof(active_nodes::AbstractVector{Int64}, dof::Int64, nnodes::Int64)
362569
n_active = length(active_nodes)
363570

@@ -394,27 +601,6 @@ function compute_matrix(nodes::AbstractVector{Int64})
394601
Correspondence_matrix_based.compute_model(nodes)
395602
end
396603

397-
mutable struct DisplacementSolverCache
398-
# Workspace
399-
F_modified::Vector{Float64}
400-
bc_mask::BitVector
401-
temp::Vector{Float64}
402-
u_free::Vector{Float64}
403-
404-
# LU cache
405-
K_free_lu::Union{Nothing,Any}
406-
last_non_BCs::Vector{Int}
407-
408-
# Active DOFs cache
409-
K_active_cached::Union{Nothing,AbstractMatrix{Float64}}
410-
last_active_dofs::Vector{Int}
411-
412-
function DisplacementSolverCache(n_total::Int)
413-
new(Float64[], BitVector(undef, n_total), Float64[], Float64[],
414-
nothing, Int[], nothing, Int[])
415-
end
416-
end
417-
418604
"""
419605
compute_displacements!(K, non_BCs, u, F, F_temp, K_reduced, lu_fact, temp)
420606
@@ -486,6 +672,7 @@ function compute_displacements!(K::AbstractMatrix{Float64},
486672
if solver.K_free_lu === nothing || solver.last_non_BCs != active_non_BCS
487673
@timeit "LU factorization" begin
488674
K_free = K_active[active_non_BCS, active_non_BCS]
675+
489676
try
490677
solver.K_free_lu = lu(K_free)
491678
catch e

0 commit comments

Comments
 (0)