From 1f776f994c5c1f435f47992fd49f5144617ddaf4 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 22 Dec 2016 10:06:18 +0200 Subject: [PATCH] possibility to solve positive semi-definite problems. --- src/solvers_modal.jl | 116 +++++++++++++++++++++++++------------------ 1 file changed, 67 insertions(+), 49 deletions(-) diff --git a/src/solvers_modal.jl b/src/solvers_modal.jl index b1f80da..970216a 100644 --- a/src/solvers_modal.jl +++ b/src/solvers_modal.jl @@ -68,32 +68,9 @@ function calc_projection(problem::Problem{Mortar}, ndim::Int) @assert C1 == C2 #@assert problem.properties.dual_basis == true @assert problem.properties.adjust == false - - # determine master and slave dofs - dim = get_unknown_field_dimension(problem) - M = Set{Int64}() - S = Set{Int64}() - slave_elements = get_slave_elements(problem) - master_elements = setdiff(get_elements(problem), slave_elements) - - for element in slave_elements - for j in get_connectivity(element) - for i=1:dim - push!(S, dim*(j-1)+i) - end - end - end - - for element in master_elements - for j in get_connectivity(element) - for i=1:dim - push!(M, dim*(j-1)+i) - end - end - end - - S = sort(collect(S)) - M = sort(collect(M)) + + S = get_nonzero_rows(C2) + M = setdiff(get_nonzero_columns(C2), S) # Construct matrix P = D^-1*M D_ = C2[S,S] @@ -133,6 +110,7 @@ function eliminate_boundary_conditions!(K_red::SparseMatrixCSC, info("Eliminating mesh tie constraint $(problem.name) using static condensation") + #= # determine master and slave dofs dim = get_unknown_field_dimension(problem) M = Set{Int64}() @@ -158,8 +136,11 @@ function eliminate_boundary_conditions!(K_red::SparseMatrixCSC, S = sort(collect(S)) M = sort(collect(M)) - N = setdiff(get_nonzero_rows(K_red), union(S, M)) - info("#S = $(length(S)), #M = $(length(M)), #N = $(length(N))") + =# + S = get_nonzero_rows(C2) + M = setdiff(get_nonzero_columns(C2), S) + #N = setdiff(get_nonzero_rows(K_red), get_nonzero_columns(C2)) + info("# slave dofs = $(length(S)), # master dofs = $(length(M))") # Construct matrix P = D^-1*M D_ = C2[S,S] @@ -221,13 +202,21 @@ function eliminate_boundary_conditions!(K_red::SparseMatrixCSC, return true end +""" +Parameters +---------- + +sigma + Shift stiffness matrix by adding diagonal term, i.e. K_shifted = K + sigma*I +""" function (solver::Solver{Modal})(; show_info=true, debug=false, bc_invertible=false, P=nothing, symmetric=true, - empty_assemblies_before_solution=true, dense=false) - show_info && info(repeat("-", 80)) - show_info && info("Starting natural frequency solver") - show_info && info("Increment time t=$(round(solver.time, 3))") - show_info && info(repeat("-", 80)) + empty_assemblies_before_solution=true, dense=false, + info_matrices=false, sigma=0.0) + info(repeat("-", 80)) + info("Starting natural frequency solver") + info("Increment time t=$(round(solver.time, 3))") + info(repeat("-", 80)) initialize!(solver) assemble!(solver; with_mass_matrix=true) M, K, Kg, f = get_field_assembly(solver) @@ -250,6 +239,7 @@ function (solver::Solver{Modal})(; show_info=true, debug=false, t1 = round(toq(), 2) info("Transform ready in $t1 seconds.") elseif nboundary_problems != 0 + tic() info("Eliminate boundary conditions from system.") for boundary_problem in get_boundary_problems(solver) eliminate_boundary_conditions!(K_red, M_red, boundary_problem, dim) @@ -273,6 +263,11 @@ function (solver::Solver{Modal})(; show_info=true, debug=false, nz = get_nonzero_rows(K_red) K_red = K_red[nz,nz] M_red = M_red[nz,nz] + + if sigma != 0.0 + info("Adding diagonal term $sigma to stiffness matrix") + end + ndofs = solver.ndofs props = solver.properties @@ -287,39 +282,62 @@ function (solver::Solver{Modal})(; show_info=true, debug=false, tic() - om2 = nothing - X = nothing - if symmetric K_red = 1/2*(K_red + transpose(K_red)) M_red = 1/2*(M_red + transpose(M_red)) end - info("is K symmetric? ", issymmetric(K_red)) - info("is M symmetric? ", issymmetric(M_red)) - info("is K positive definite? ", isposdef(K_red)) - info("is M positive definite? ", isposdef(M_red)) + if info_matrices + info("is K symmetric? ", issymmetric(K_red)) + info("is M symmetric? ", issymmetric(M_red)) + info("is K positive definite? ", isposdef(K_red)) + info("is M positive definite? ", isposdef(M_red)) + end if dense K_red = full(K_red) M_red = full(M_red) end + + om2 = nothing + X = nothing + passed = false + try - om2, X = eigs(K_red, M_red; nev=props.nev, which=props.which) + om2, X = eigs(K_red + sigma*I, M_red; nev=props.nev, which=props.which) + passed = true catch - info("failed to calculate eigenvalues") - info("reduced system") - info("is K symmetric? ", issym(K_red)) - info("is M symmetric? ", issym(M_red)) + info("failed to calculate eigenvalues for problem. Maybe stiffness matrix is not positive definite, checking...") + info("is K symmetric? ", issymmetric(K_red)) + info("is M symmetric? ", issymmetric(M_red)) info("is K positive definite? ", isposdef(K_red)) info("is M positive definite? ", isposdef(M_red)) - dump(full(K_red[1:10,1:10])) + info("Probably the reason is that stiffness matrix is not positive definite and Cholesky factorization is failing.") + info("To work around this problem, use arguments `sigma = ` when calling solver, i.e.") + info("solver(; sigma=1.0e-9") + info("Be aware that using sigma shifts eigenvalues up and a bit different results can be expected.") if size(K_red, 1) < 2000 + info("stiffness matrix is small, using dense eigenvalue solver to check eigenvalues ...") om2 = eigvals(full(K_red)) - info("om2 = $om2") + info("squared eigenvalues om2 = $om2") + end + if sigma != 0.0 + info("sigma is manually set and did not work, giving up, try increase sigma.") + rethrow() + end + end + + if !passed + sigma = 1.0e-9 + info("Calculation of eigenvalues failed, trying again using sigma value sigma=$sigma") + try + om2, X = eigs(K_red + sigma*I, M_red; nev=props.nev, which=props.which) + passed = true + catch + info("Failed to calculate eigenvalues with sigma=$sigma, manually set sigma to something larger and try again.") + rethrow() end - rethrow() end t1 = round(toq(), 2) @@ -335,7 +353,7 @@ function (solver::Solver{Modal})(; show_info=true, debug=false, S, M, P = calc_projection(problem, dim) # FIXME: store projection to boundary problem, i.e. # update!(problem, "master-slave projection", time => P) - # us = P*um + # us = P*um props.eigvecs[S,i] = P*props.eigvecs[M,i] end end