From 2e086c43ab9f7c386e05a407fe8914ec7ee2f3e2 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 13 Oct 2016 01:00:08 +0300 Subject: [PATCH] improved eigenvalue solver code --- src/solvers_modal.jl | 148 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 127 insertions(+), 21 deletions(-) diff --git a/src/solvers_modal.jl b/src/solvers_modal.jl index 34e0c3e..6188ac7 100644 --- a/src/solvers_modal.jl +++ b/src/solvers_modal.jl @@ -23,24 +23,15 @@ function Modal(nev=10, which=:SM) solver = Modal(false, Vector(), Matrix(), nev, which) end -function call(solver::Solver{Modal}; show_info=true, debug=false, bc_invertible=false) +function call(solver::Solver{Modal}; show_info=true, debug=false, + bc_invertible=false, P=nothing, + empty_assemblies_before_solution=true) 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)) initialize!(solver) - # assemble all field problems - info("Assembling problems ...") - tic() - for problem in get_field_problems(solver) - assemble!(problem, solver.time) - assemble!(problem, solver.time, Val{:mass_matrix}) - end - for problem in get_boundary_problems(solver) - assemble!(problem, solver.time) - end - t1 = round(toq(), 2) - info("Assembled in $t1 seconds.") + assemble!(solver; with_mass_matrix=true) M, K, Kg, f = get_field_assembly(solver) Kb, C1, C2, D, fb, g = get_boundary_assembly(solver) K = K + Kb @@ -49,6 +40,14 @@ function call(solver::Solver{Modal}; show_info=true, debug=false, bc_invertible= K += Kg end + # free up some memory before solution + if empty_assemblies_before_solution + for problem in get_problems(solver) + empty!(problem.assembly) + end + gc() + end + @assert nnz(D) == 0 @assert C1 == C2 @@ -56,26 +55,38 @@ function call(solver::Solver{Modal}; show_info=true, debug=false, bc_invertible= nboundary_problems = length(get_boundary_problems(solver)) - if nboundary_problems != 0 + if !(P == nothing) + info("using custom P") + K_red = P'*K*P + M_red = P'*M*P + elseif nboundary_problems != 0 if bc_invertible + info("Invertible C, calculating P") P, h = create_projection(C1, g, Val{:invertible}) else + info("Contacts, calculate P") P, h = create_projection(C1, g) end + K_red = P'*K*P + M_red = P'*M*P else + info("No dirichlet boundaryes, P = I") P = speye(size(K, 1)) + K_red = K + M_red = M end - K_red = P'*K*P - M_red = P'*M*P - - # make sure matrices are symmetric - K_red = 1/2*(K_red + K_red') - M_red = 1/2*(M_red + M_red') - t1 = round(toq(), 2) info("Eliminated dirichlet boundaries in $t1 seconds.") + # make sure matrices are symmetric + info("Making matrices symmetric") + tic() + K_red = 1/2*(K_red + K_red') + M_red = 1/2*(M_red + M_red') + t1 = round(toq(), 2) + info("Finished in $t1 seconds.") + nz = get_nonzero_rows(K_red) ndofs = solver.ndofs props = solver.properties @@ -116,8 +127,10 @@ function call(solver::Solver{Modal}; show_info=true, debug=false, bc_invertible= rethrow() end + t1 = round(toq(), 2) info("Eigenvalues computed in $t1 seconds. Eigenvalues: $om2") + tic() props.eigvals = om2 props.eigvecs = zeros(ndofs, length(om2)) v = zeros(ndofs) @@ -128,6 +141,7 @@ function call(solver::Solver{Modal}; show_info=true, debug=false, bc_invertible= end t1 = round(toq(), 2) +#= for i=1:length(om2) freq = real(sqrt(om2[i])/(2.0*pi)) u = props.eigvecs[:,i] @@ -145,7 +159,99 @@ function call(solver::Solver{Modal}; show_info=true, debug=false, bc_invertible= end end end +=# return true end +function update_xdmf!(solver::Solver{Modal}; show_info=true) + xdmf = get(solver.xdmf) + temporal_collection = get_temporal_collection(xdmf) + + frame = new_element("Grid") + new_child(frame, "Time", Dict("Value" => solver.time)) + + # save geometry + X = solver("geometry", solver.time) + node_ids = sort(collect(keys(X))) + geometry = hcat([X[nid] for nid in node_ids]...) + ndim, nnodes = size(geometry) + geom_type = ndim == 2 ? "XY" : "XYZ" + dataitem = new_dataitem(xdmf, "/Node IDs", node_ids) + geom = new_child(frame, "Geometry", Dict("Type" => geom_type)) + dataitem = new_dataitem(xdmf, "/Geometry", geometry) + add_child(geom, dataitem) + + # save topology + all_elements = get_all_elements(solver) + nelements = length(all_elements) + element_types = unique(map(get_element_type, all_elements)) + + xdmf_element_mapping = Dict( + "Seg2" => "Polyline", + "Tri3" => "Triangle", + "Quad4" => "Quadrilateral", + "Tet4" => "Tetrahedron", + "Pyramid5" => "Pyramid", + "Wedge6" => "Wedge", + "Hex8" => "Hexahedron", + "Seg3" => "Edge_3", + "Tri6" => "Tri_6", + "Quad8" => "Quad_8", + "Tet10" => "Tet_10", + "Pyramid13" => "Pyramid_13", + "Wedge15" => "Wedge_15", + "Hex20" => "Hex_20") + + for element_type in element_types + elements = filter_by_element_type(element_type, all_elements) + sort!(elements, by=get_element_id) + element_ids = map(get_element_id, elements) + element_conn = map(get_connectivity, elements) + element_conn = transpose(hcat(element_conn...)) - 1 + element_code = split(string(element_type), ".")[end] + dataitem = new_dataitem(xdmf, "/Topology/$element_code/Element IDs", element_ids) + dataitem = new_dataitem(xdmf, "/Topology/$element_code/Connectivity", element_conn) + topology = new_child(frame, "Topology") + set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code]) + set_attribute(topology, "NumberOfElements", length(elements)) + add_child(topology, dataitem) + end + + # save solved fields + unknown_field_name = get_unknown_field_name(solver) + U = solver(unknown_field_name, solver.time) + node_ids2 = sort(collect(keys(U))) + + @assert node_ids == node_ids2 + ndim = length(U[first(node_ids)]) + field_type = ndim == 1 ? "Scalar" : "Vector" + field_center = "Node" + if ndim == 2 + for nid in node_ids + U[nid] = [U[nid]; 0.0] + end + ndim = 3 + end + U = hcat([U[nid] for nid in node_ids]...) + unknown_field_name = ucfirst(unknown_field_name) + time = solver.time + path = "" + if S == Nonlinear + iteration = solver.properties.iteration + path = "/Results/Time $time/Iteration $iteration/Nodal Fields/$unknown_field_name" + elseif S == Linear + path = "/Results/Time $time/Nodal Fields/$unknown_field_name" + end + dataitem = new_dataitem(xdmf, path, U) + attribute = new_child(frame, "Attribute") + set_attribute(attribute, "Name", unknown_field_name) + set_attribute(attribute, "Center", field_center) + set_attribute(attribute, "AttributeType", field_type) + add_child(attribute, dataitem) + if (S == Linear) || ((S == Nonlinear) && has_converged(solver)) + add_child(temporal_collection, frame) + end + save!(xdmf) +end +