Skip to content

Commit 6e60f80

Browse files
CWillbergCWillberg
committed
code optimization
1 parent ab14261 commit 6e60f80

File tree

5 files changed

+193
-129
lines changed

5 files changed

+193
-129
lines changed

‎src/Core/BC_manager.jl‎

Lines changed: 102 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -182,28 +182,28 @@ function apply_bc_dirichlet(allowed_variables::Vector{String},
182182
if haskey(dof_mapping, bc["Coordinate"])
183183
@views field_to_apply_bc = field[bc["Node Set"],
184184
dof_mapping[bc["Coordinate"]]]
185-
eval_bc!(field_to_apply_bc,
186-
bc["Value"],
187-
coordinates[bc["Node Set"], :],
188-
time,
189-
step_time,
190-
dof,
191-
bc["Initial"],
192-
name)
185+
bc["Value"] = eval_bc!(field_to_apply_bc,
186+
bc["Value"],
187+
coordinates[bc["Node Set"], :],
188+
time,
189+
step_time,
190+
dof,
191+
bc["Initial"],
192+
name)
193193
else
194194
@error "Coordinate in boundary condition must be x,y or z."
195195
return nothing
196196
end
197197
else
198198
@views field_to_apply_bc = field[bc["Node Set"]]
199-
eval_bc!(field_to_apply_bc,
200-
bc["Value"],
201-
coordinates[bc["Node Set"], :],
202-
time,
203-
step_time,
204-
dof,
205-
bc["Initial"],
206-
name)
199+
bc["Value"] = eval_bc!(field_to_apply_bc,
200+
bc["Value"],
201+
coordinates[bc["Node Set"], :],
202+
time,
203+
step_time,
204+
dof,
205+
bc["Initial"],
206+
name)
207207
end
208208
end
209209
end
@@ -233,30 +233,30 @@ function apply_bc_neumann(bcs::Dict, time::Float64, step_time::Float64)
233233
if haskey(dof_mapping, bc["Coordinate"])
234234
@views field_to_apply_bc = field[bc["Node Set"],
235235
dof_mapping[bc["Coordinate"]]]
236-
eval_bc!(field_to_apply_bc,
237-
bc["Value"],
238-
coordinates[bc["Node Set"], :],
239-
time,
240-
step_time,
241-
dof,
242-
bc["Initial"],
243-
name,
244-
true)
236+
bc["Value"] = eval_bc!(field_to_apply_bc,
237+
bc["Value"],
238+
coordinates[bc["Node Set"], :],
239+
time,
240+
step_time,
241+
dof,
242+
bc["Initial"],
243+
name,
244+
true)
245245
else
246246
@error "Coordinate in boundary condition must be x,y or z."
247247
return nothing
248248
end
249249
else
250250
@views field_to_apply_bc = field[bc["Node Set"]]
251-
eval_bc!(field_to_apply_bc,
252-
bc["Value"],
253-
coordinates[bc["Node Set"], :],
254-
time,
255-
step_time,
256-
dof,
257-
bc["Initial"],
258-
name,
259-
true)
251+
bc["Value"] = eval_bc!(field_to_apply_bc,
252+
bc["Value"],
253+
coordinates[bc["Node Set"], :],
254+
time,
255+
step_time,
256+
dof,
257+
bc["Initial"],
258+
name,
259+
true)
260260
end
261261
end
262262
end
@@ -345,6 +345,7 @@ function eval_bc!(field_values::Union{SubArray,NodeScalarField{Float64},
345345
(coordinates[:, 1], coordinates[:, 2], coordinates[:, 3],
346346
time,
347347
step_time)...)
348+
bc = dynamic_bc_3D_func
348349
else
349350
func_args = [:x, :y, :t, :st]
350351
dynamic_func_expr = quote
@@ -355,6 +356,7 @@ function eval_bc!(field_values::Union{SubArray,NodeScalarField{Float64},
355356
(coordinates[:, 1], coordinates[:, 2],
356357
time,
357358
step_time)...)
359+
bc = dynamic_2D_bc_func
358360
end
359361

360362
if isnothing(value) || (initial && time != 0.0)
@@ -370,6 +372,72 @@ function eval_bc!(field_values::Union{SubArray,NodeScalarField{Float64},
370372
else
371373
copyto!(field_values, value)
372374
end
375+
return bc
376+
end
377+
378+
"""
379+
eval_bc!(field_values::Union{NodeScalarField{Float64},NodeScalarField{Int64}}, bc::Union{Float64,Float64,Int64,String}, coordinates::Matrix{Float64}, time::Float64, dof::Int64)
380+
Working with if-statements
381+
"if t>2 0 else 20 end"
382+
works for scalars. If you want to evaluate a vector, please use the Julia notation as input
383+
"ifelse.(x .> y, 10, 20)"
384+
"""
385+
function eval_bc!(field_values::Union{SubArray,NodeScalarField{Float64},
386+
NodeScalarField{Int64}},
387+
bc::Function,
388+
coordinates::Union{Matrix{Float64},Matrix{Int64}},
389+
time::Float64,
390+
step_time::Float64,
391+
dof::Int64,
392+
initial::Bool,
393+
name::String = "BC_1",
394+
neumann::Bool = false)
395+
# reason for global
396+
# https://stackoverflow.com/questions/60105828/julia-local-variable-not-defined-in-expression-eval
397+
# the yaml input allows multiple types. But for further use this input has to be a string
398+
399+
if length(coordinates) == 0
400+
# @warn "Ignoring boundary condition $name.\n No nodes found, check Input Deck and or Node Sets."
401+
return
402+
end
403+
404+
if dof > 2
405+
#func_args = [:x, :y, :z, :t, :st]
406+
#dynamic_func_expr = quote
407+
# ($(func_args...),) -> $bc_value
408+
#end
409+
#dynamic_bc_3D_func = Base.eval(@__MODULE__, dynamic_func_expr)
410+
411+
value = Base.invokelatest(bc,
412+
(coordinates[:, 1], coordinates[:, 2], coordinates[:, 3],
413+
time,
414+
step_time)...)
415+
else
416+
#func_args = [:x, :y, :t, :st]
417+
#dynamic_func_expr = quote
418+
# ($(func_args...),) -> $bc_value
419+
#end
420+
#dynamic_2D_bc_func = Base.eval(@__MODULE__, dynamic_func_expr)
421+
value = Base.invokelatest(bc,
422+
(coordinates[:, 1], coordinates[:, 2],
423+
time,
424+
step_time)...)
425+
end
426+
427+
if isnothing(value) || (initial && time != 0.0)
428+
return
429+
end
430+
431+
if value isa Number
432+
if neumann
433+
field_values .+= value
434+
else
435+
fill!(field_values, value)
436+
end
437+
else
438+
copyto!(field_values, value)
439+
end
440+
return bc
373441
end
374442

375443
end

‎src/Core/Solver/Matrix_Verlet.jl‎

Lines changed: 78 additions & 82 deletions
Original file line numberDiff line numberDiff line change
@@ -189,10 +189,10 @@ function run_solver(solver_options::Dict{Any,Any},
189189
bcs::Dict{Any,Any},
190190
outputs::Dict{Int64,Dict{}},
191191
result_files::Vector{Dict},
192-
synchronise_field,
193-
write_results,
194-
compute_parabolic_problems_before_model_evaluation,
195-
compute_parabolic_problems_after_model_evaluation,
192+
synchronise_field::Function,
193+
write_results::Function,
194+
compute_parabolic_problems_before_model_evaluation::Function,
195+
compute_parabolic_problems_after_model_evaluation::Function,
196196
silent::Bool)
197197
dt::Float64 = solver_options["dt"]
198198
nsteps::Int64 = solver_options["Number of Steps"]
@@ -203,17 +203,17 @@ function run_solver(solver_options::Dict{Any,Any},
203203
max_damage::Float64 = 0
204204
damage_init::Bool = false
205205

206-
density = Data_Manager.get_field("Density")
206+
density::NodeScalarField{Float64} = Data_Manager.get_field("Density")
207207

208-
coor = Data_Manager.get_field("Coordinates")
208+
coor::NodeVectorField{Float64} = Data_Manager.get_field("Coordinates")
209209

210-
active_list = Data_Manager.get_field("Active")
211-
volume = Data_Manager.get_field("Volume")
210+
active_list::NodeScalarField{Bool} = Data_Manager.get_field("Active")
211+
volume::NodeScalarField{Float64} = Data_Manager.get_field("Volume")
212212

213213
if "Material" in solver_options["Models"]
214-
external_forces = Data_Manager.get_field("External Forces")
215-
external_force_densities = Data_Manager.get_field("External Force Densities")
216-
a = Data_Manager.get_field("Acceleration")
214+
external_forces::NodeVectorField{Float64} = Data_Manager.get_field("External Forces")
215+
external_force_densities::NodeVectorField{Float64} = Data_Manager.get_field("External Force Densities")
216+
a::NodeVectorField{Float64} = Data_Manager.get_field("Acceleration")
217217
end
218218

219219
comm = Data_Manager.get_comm()
@@ -233,86 +233,81 @@ function run_solver(solver_options::Dict{Any,Any},
233233

234234
@timeit "Matrix Verlet" begin
235235
@inbounds @fastmath for idt in iter
236-
uN::Matrix{Float64} = Data_Manager.get_field("Displacements", "N")
237-
uNP1::Matrix{Float64} = Data_Manager.get_field("Displacements", "NP1")
238-
forces::Matrix{Float64} = Data_Manager.get_field("Forces", "NP1")
239-
deformed_coorNP1::Matrix{Float64} = Data_Manager.get_field("Deformed Coordinates",
240-
"NP1")
241-
vN::Matrix{Float64} = Data_Manager.get_field("Velocity", "N")
242-
vNP1::Matrix{Float64} = Data_Manager.get_field("Velocity", "NP1")
243-
force_densities_NP1::Matrix{Float64} = Data_Manager.get_field("Force Densities",
244-
"NP1")
245-
active_nodes::Vector{Int64} = Data_Manager.get_field("Active Nodes")
246-
active_nodes = find_active_nodes(active_list, active_nodes,
247-
1:Data_Manager.get_nnodes())
248-
apply_bc_dirichlet(["Displacements", "Forces", "Force Densities"],
249-
bcs, time,
250-
step_time)
251-
252-
if solver_options["Model Reduction"] != false
253-
active_nodes = master_nodes
254-
active_list .= false
255-
active_list[Data_Manager.get_reduced_model_pd()] .= true
256-
end
257-
258-
vNP1[active_nodes, :] = (1 - numerical_damping) .*
259-
vN[active_nodes, :] .+
260-
0.5 * dt .* a[active_nodes, :]
236+
@timeit "solver start" begin
237+
uN::Matrix{Float64} = Data_Manager.get_field("Displacements", "N")
238+
uNP1::Matrix{Float64} = Data_Manager.get_field("Displacements", "NP1")
239+
forces::Matrix{Float64} = Data_Manager.get_field("Forces", "NP1")
240+
deformed_coorNP1::Matrix{Float64} = Data_Manager.get_field("Deformed Coordinates",
241+
"NP1")
242+
vN::Matrix{Float64} = Data_Manager.get_field("Velocity", "N")
243+
vNP1::Matrix{Float64} = Data_Manager.get_field("Velocity", "NP1")
244+
force_densities_NP1::Matrix{Float64} = Data_Manager.get_field("Force Densities",
245+
"NP1")
246+
active_nodes::Vector{Int64} = Data_Manager.get_field("Active Nodes")
247+
@timeit "active nodes" active_nodes=find_active_nodes(active_list,
248+
active_nodes,
249+
1:Data_Manager.get_nnodes())
261250

262-
uNP1[active_nodes, :] = uN[active_nodes, :] .+
263-
dt .* vNP1[active_nodes, :]
251+
if solver_options["Model Reduction"] != false
252+
active_nodes = master_nodes
253+
active_list .= false
254+
active_list[Data_Manager.get_reduced_model_pd()] .= true
255+
end
256+
end
257+
@timeit "compute Velocity" begin
258+
@views vNP1[active_nodes, :] = (1 - numerical_damping) .*
259+
vN[active_nodes, :] .+
260+
0.5 * dt .* a[active_nodes, :]
264261

265-
apply_bc_dirichlet(["Displacements", "Temperature"],
266-
bcs,
267-
time,
268-
step_time) #-> Dirichlet
269-
@views deformed_coorNP1 .= coor .+ uNP1
262+
@views uNP1[active_nodes, :] = uN[active_nodes, :] .+
263+
dt .* vNP1[active_nodes, :]
264+
end
265+
@timeit "apply BC" apply_bc_dirichlet(["Displacements", "Temperature"],
266+
bcs,
267+
time,
268+
step_time) #-> Dirichlet
269+
@timeit "deformed coor" @views deformed_coorNP1 .= coor .+ uNP1
270270
# thermal model, usw.
271271

272272
#mul!(a[active_nodes], K[active_nodes, active_nodes], uNP1[active_nodes])
273273

274-
sa = size(a[active_nodes, :])
275-
276-
force_densities_NP1 .= 0.0
274+
@timeit "sa" sa=size(a[active_nodes, :])
277275

278-
compute_models(block_nodes,
279-
dt,
280-
time,
281-
solver_options["Models"],
282-
synchronise_field)
276+
@timeit "compute models" compute_models(block_nodes,
277+
dt,
278+
time,
279+
solver_options["Models"],
280+
synchronise_field)
283281
#force_densities_NP1[active_nodes, :] .*= volume[active_nodes]
284-
@timeit "Force matrix computations" begin
282+
@timeit "Force computations" begin
285283
# check if valid if volume is different
286284

287285
@views fNP1 = force_densities_NP1[active_nodes, :]
288286

289287
#-= f_int(K,
290288
# vec(uNP1[active_nodes,
291289
# :]), sa)
292-
f_int_inplace!(fNP1, temp, K, vec(uNP1[active_nodes, :]), sa)
293-
fNP1 .+= external_force_densities[active_nodes,
294-
:] .+
295-
external_forces[active_nodes,
296-
:] ./
297-
volume[active_nodes]
290+
@timeit "Force matrix computations" f_int_inplace!(fNP1, temp, K,
291+
vec(uNP1[active_nodes,
292+
:]), sa)
293+
@views fNP1 .+= external_force_densities[active_nodes,
294+
:] .+
295+
external_forces[active_nodes,
296+
:] ./ volume[active_nodes]
297+
298+
@views forces[active_nodes, :] .= force_densities_NP1[active_nodes, :] .*
299+
volume[active_nodes]
298300
end
299-
# @timeit "download_from_cores" Data_Manager.synch_manager(synchronise_field,
300-
# "download_from_cores")
301-
#println(a)
302-
#@views a[active_nodes, :] = external_force_densities[active_nodes, :] .+
303-
# external_forces[active_nodes, :] ./
304-
# volume[active_nodes]
305-
forces[active_nodes, :] .= force_densities_NP1[active_nodes, :] .*
306-
volume[active_nodes]
307301
@timeit "Accelaration computation" begin
308302
if solver_options["Model Reduction"] != false
309-
a[active_nodes, :] = reshape(M_fact \
310-
vec(force_densities_NP1[active_nodes, :] .*
311-
volume[active_nodes]),
312-
sa...)
303+
@views a[active_nodes, :] = reshape(M_fact \
304+
vec(force_densities_NP1[active_nodes,
305+
:] .*
306+
volume[active_nodes]),
307+
sa...)
313308
else
314-
a[active_nodes, :] = force_densities_NP1[active_nodes, :] ./
315-
(density[active_nodes])
309+
@views a[active_nodes, :] = force_densities_NP1[active_nodes, :] ./
310+
(density[active_nodes])
316311
end
317312
end
318313
@timeit "write_results" result_files=write_results(result_files, time,
@@ -325,19 +320,20 @@ function run_solver(solver_options::Dict{Any,Any},
325320
break
326321
end
327322
@timeit "switch_NP1_to_N" Data_Manager.switch_NP1_to_N()
323+
@timeit "final solver steps" begin
324+
time += dt
325+
step_time += dt
326+
Data_Manager.set_current_time(time)
328327

329-
time += dt
330-
step_time += dt
331-
Data_Manager.set_current_time(time)
328+
if idt % ceil(nsteps / 100) == 0
329+
@info "Step: $idt / $(nsteps+1) [$time s]"
330+
end
331+
if rank == 0 && !silent
332+
set_postfix(iter, t = @sprintf("%.4e", time))
333+
end
332334

333-
if idt % ceil(nsteps / 100) == 0
334-
@info "Step: $idt / $(nsteps+1) [$time s]"
335-
end
336-
if rank == 0 && !silent
337-
set_postfix(iter, t = @sprintf("%.4e", time))
335+
barrier(comm)
338336
end
339-
340-
barrier(comm)
341337
#a+= + fexternal / density
342338
end
343339
end

0 commit comments

Comments
 (0)