Skip to content

Commit

Permalink
Timing xdmf write for modal solver
Browse files Browse the repository at this point in the history 
Writing results to Xdmf in Modal analysis is unacceptable slow.
Finding reason and fixing.
  • Loading branch information
@ahojukka5
ahojukka5 committed Aug 20, 2017
1 parent 0ba3e5e commit fbe9060
Showing 1 changed file with 61 additions and 121 deletions.
182 changes: 61 additions & 121 deletions src/solvers_modal.jl
View file
Expand Up @@ -111,39 +111,10 @@ 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}()
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)
#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]
M_ = -C2[S,M]

Expand All @@ -154,52 +125,18 @@ function eliminate_boundary_conditions!(K_red::SparseMatrixCSC,
else
P = D_ \ M_
end
#@assert isdiag(D_)
info("Matrix P ready.")
Id = ones(ndim)
#Id[S] = 0
#Id[M] = 0
Q = spdiagm(Id)
Q[M,S] += P'
info("Matrix Q ready.")

# testing
#K_red_orig = copy(K_red)
#K_red[N,M] += K_red[N,S]*P
#K_red[M,N] += P'*K_red[S,N]
#K_red[M,M] += P'*K_red[S,S]*P
#K_red[S,:] = 0.0
#K_red[:,S] = 0.0

info("K transform")

K_red[:,:] = Q*K_red*Q'
K_red[S,:] = 0.0
K_red[:,S] = 0.0

#K_res = K_red - K_red_2
#SparseArrays.droptol!(K_res, 1.0e-9)

info("K transform ready")
#info("Create matrices, M")
#info("Sum matricse, M")
#M_red_orig = copy(M_red)
#M_red[N,M] += M_red[N,S]*P
#M_red[M,N] += P'*M_red[S,N]
#M_red[M,M] += P'*M_red[S,S]*P
#M_red[S,:] = 0.0
#M_red[:,S] = 0.0
info("M transform")
M_red[:,:] = Q*M_red*Q'
M_red[S,:] = 0.0
M_red[:,S] = 0.0
info("M transform ready")

#M_res = M_red - M_red_2
#SparseArrays.droptol!(M_res, 1.0e-9)
#info("Diff")
#println(K_res)
#println(M_res)

return true
end

Expand Down Expand Up @@ -366,24 +303,25 @@ function update_xdmf!(solver::Solver{Modal})
end

if maximum(abs.(imag(solver.properties.eigvals))) > 1.0e-9
error("Writing imaginary eigenvalues for Xdmf not supported.")
info("Writing imaginary eigenvalues for Xdmf not supported.")
return
end

xdmf = get(solver.xdmf)

# geometry
X_ = solver("geometry", solver.time)
node_ids = sort(collect(keys(X_)))
X = hcat([X_[nid] for nid in node_ids]...)
ndim, nnodes = size(X)
geom_type = (ndim == 2 ? "XY" : "XYZ")

# topology
nid_mapping = Dict(j=>i for (i, j) in enumerate(node_ids))
all_elements = get_all_elements(solver)
nelements = length(all_elements)
debug("Saving topology: $nelements elements total.")
element_types = unique(map(get_element_type, all_elements))
@timeit "extract geometry" begin
@timeit "fetch geometry" X_ = solver("geometry", solver.time)
@timeit "sort node_ids" node_ids = sort(collect(keys(X_)))
@timeit "create array" X = hcat([X_[nid] for nid in node_ids]...)
ndim, nnodes = size(X)
geom_type = (ndim == 2 ? "XY" : "XYZ")
end
@timeit "extract topology" begin
@timeit "create node mappign" nid_mapping = Dict(j=>i for (i, j) in enumerate(node_ids))
@timeit "get all elements" all_elements = get_all_elements(solver)
@timeit "find unique element types" element_types = unique(map(get_element_type, all_elements))
nelements = length(all_elements)
end

xdmf_element_mapping = Dict(
"Poi1" => "Polyvertex",
Expand All @@ -405,10 +343,10 @@ function update_xdmf!(solver::Solver{Modal})
# save modes

temporal_collection = get_temporal_collection(xdmf)

unknown_field_name = ucfirst(get_unknown_field_name(solver))
frames = []
for (j, eigval) in enumerate(real(solver.properties.eigvals))

@timeit "save modes" for (j, eigval) in enumerate(real(solver.properties.eigvals))
if eigval < 0.0
warn("negative real eigenvalue found, om2=$eigval, setting to zero.")
eigval = 0.0
Expand All @@ -421,53 +359,55 @@ function update_xdmf!(solver::Solver{Modal})
time = new_child(frame, "Time")
set_attribute(time, "Value", freq)

# add geometry
geometry = new_element("Geometry")
set_attribute(geometry, "Type", geom_type)
data_node_ids = new_dataitem(xdmf, "/Node IDs", node_ids)
data_geometry = new_dataitem(xdmf, "/Geometry", X)
add_child(geometry, data_geometry)
add_child(frame, geometry)
@timeit "add geometry" begin
geometry = new_element("Geometry")
set_attribute(geometry, "Type", geom_type)
data_node_ids = new_dataitem(xdmf, "/Node IDs", node_ids)
data_geometry = new_dataitem(xdmf, "/Geometry", X)
add_child(geometry, data_geometry)
add_child(frame, geometry)
end

# add topology
for element_type in element_types
elements = filter_by_element_type(element_type, all_elements)
nelements = length(elements)
debug("Xdmf save: $nelements elements of type $element_type")
sort!(elements, by=get_element_id)
element_ids = map(get_element_id, elements)
element_conn = map(element -> [nid_mapping[j]-1 for j in get_connectivity(element)], elements)
element_conn = hcat(element_conn...)
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_element("Topology")
set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code])
set_attribute(topology, "NumberOfElements", length(elements))
add_child(topology, dataitem)
add_child(frame, topology)
timeit("save topology of element type $element_type") do
@timeit "filter by element type" elements = filter_by_element_type(element_type, all_elements)
nelements = length(elements)
@timeit "sort elements" sort!(elements, by=get_element_id)
@timeit "get element ids" element_ids = map(get_element_id, elements)
@timeit "creaet element_conn" element_conn = map(element -> [nid_mapping[j]-1 for j in get_connectivity(element)], elements)
@timeit "hcat elcon" element_conn = hcat(element_conn...)
@timeit "element codd" element_code = split(string(element_type), ".")[end]
@timeit "new dataitem elids" dataitem = new_dataitem(xdmf, "/Topology/$element_code/Element IDs", element_ids)
@timeit "new dataitem elcon" dataitem = new_dataitem(xdmf, "/Topology/$element_code/Connectivity", element_conn)
@timeit "new element" topology = new_element("Topology")
set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code])
set_attribute(topology, "NumberOfElements", length(elements))
@timeit "add topo child" add_child(topology, dataitem)
@timeit "add topo frame" add_child(frame, topology)
end
end

mode = zeros(X)
mode_ = solver.properties.eigvecs[:,j]
mode_ = reshape(mode_, ndim, round(Int, length(mode_)/ndim))
for nid in node_ids
loc = nid_mapping[nid]
mode[:,loc] = mode_[:,nid]
@timeit "store eigenmode" begin
@timeit "init X" mode = zeros(X)
@timeit "get mode" mode_ = solver.properties.eigvecs[:,j]
@timeit "reshape mode" mode_ = reshape(mode_, ndim, round(Int, length(mode_)/ndim))
@timeit "pick by nid" for nid in node_ids
loc = nid_mapping[nid]
mode[:,loc] = mode_[:,nid]
end
field_type = ndim == 1 ? "Scalar" : "Vector"
field_center = "Node"
attribute = new_child(frame, "Attribute")
set_attribute(attribute, "Name", unknown_field_name)
set_attribute(attribute, "Center", field_center)
set_attribute(attribute, "AttributeType", field_type)
@timeit "new dataitem mode" dataitem = new_dataitem(xdmf, path, mode)
add_child(attribute, dataitem)
add_child(frame, attribute)
add_child(temporal_collection, frame)
end

field_type = ndim == 1 ? "Scalar" : "Vector"
field_center = "Node"
attribute = new_child(frame, "Attribute")
set_attribute(attribute, "Name", unknown_field_name)
set_attribute(attribute, "Center", field_center)
set_attribute(attribute, "AttributeType", field_type)
dataitem = new_dataitem(xdmf, path, mode)
add_child(attribute, dataitem)
add_child(frame, attribute)
add_child(temporal_collection, frame)
end

info("Saving Xdmf")

save!(xdmf)
end

0 comments on commit fbe9060

Please sign in to comment.