Skip to content

Commit 5486242

Browse files
LasNikasmarcelschurersvchbLasNikas
authored
Restart from VTU (#1037)
* Update vtk2trixi function to return `NamedTuple` instead of `InitialCondition` * Refactor vtk2trixi function to accept custom fields * Enhance vtk2trixi function to support scalar and vector custom quantities and update tests * Fix typo in documentation for vtk2trixi function * Fix capitalization in comments * add ELTYPE * Fix formatting in documentation * add `coordinates_eltype` * Corrected failing doctest output * first prototype * write Q * gpu fix * add checks * add docs * Corrected failing doctest output * adapt pressure model * use `foreach_no_alloc` * fix unit tests * fix unit tests * fix unit tests * fix doc tests * better docs * use stored type * add docs * add tests * fix formatting * fix eltype * revise #959 * fix buffer * first poc * add show * remove mass * add preconditioning * fix inactive particles * fix zone IDs * add test * make it gpu compatible * add gpu tests * fix tests * add checks * add example * add docs * fix typo * implement suggestions * implement suggestions * implement suggestions * use strings instead of RestartConditin * add test mixed types * add NEWS entry * revise * fix doc test * rm test file * fix merge bugs * fix preconditioning * don't read IC * fix merge bug * fix NEWS entry * fix merge bug * fix formatting * better name * add more tests * revise example * implement suggestions * add structure example * add tests * fix description * fix tests * avoid outputs * add NEWS entry * implement suggestions * implement suggestions * fix gpu * apply formatter * fix * remove beam example * implement suggestions * undo density_calculator dispatch * use non-adaptive scheme * apply formatter * fix cfl * fix restart --------- Co-authored-by: Marcel Schurer <marcel.schurer@outlook.de> Co-authored-by: Sven Berger <berger.sven@gmail.com> Co-authored-by: LasNikas <niklas.nehe@web.de>
1 parent 2b76eaf commit 5486242

18 files changed

Lines changed: 573 additions & 15 deletions

File tree

NEWS.md

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,9 +4,16 @@ TrixiParticles.jl follows the interpretation of
44
[semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
55
used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
66

7+
## Version 0.5.2
8+
9+
### Features
10+
- Added the ability to restart simulations from VTK solution files generated by `SolutionSavingCallback`.
11+
Users can now pass the `restart_with` keyword argument to `semidiscretize` (#1190).
12+
713
## Version 0.5.1
814

915
### Features
16+
1017
- Implemented stage-level coupling for split integration (#1049).
1118
- Added `MechanicalWorkCalculatorCallback` (#940).
1219

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
# ==========================================================================================
2+
# Restart Example: Poiseuille Flow 2D
3+
#
4+
# This example demonstrates how to restart a simulation.
5+
# We first run a simulation of 2D Poiseuille flow up to t=0.3s, then restart from the
6+
# saved state and continue the simulation until t=0.6s.
7+
# ==========================================================================================
8+
using TrixiParticles
9+
10+
trixi_include(@__MODULE__,
11+
joinpath(examples_dir(), "fluid", "poiseuille_flow_2d.jl"),
12+
tspan=(0.0, 0.3), sound_speed_factor=10, particle_spacing=4e-5)
13+
14+
# Get latest iteration
15+
iter = saving_callback.affect!.affect!.latest_saved_iter
16+
17+
restart_file_fluid = joinpath("out", "fluid_1_$iter.vtu")
18+
restart_file_open_boundary = joinpath("out", "open_boundary_1_$iter.vtu")
19+
restart_file_boundary = joinpath("out", "boundary_1_$iter.vtu")
20+
21+
ode_restart = semidiscretize(semi, (0.3, 0.6);
22+
restart_with=(restart_file_fluid,
23+
restart_file_open_boundary,
24+
restart_file_boundary))
25+
26+
saving_callback = SolutionSavingCallback(dt=0.02, prefix="restart")
27+
28+
callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback())
29+
30+
sol_restart = solve(ode_restart, RDPK3SpFSAL35(), abstol=1e-5, reltol=1e-3, dtmax=1e-2,
31+
save_everystep=false, callback=callbacks)

src/TrixiParticles.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -62,6 +62,7 @@ include("general/semidiscretization.jl")
6262
include("general/gpu.jl")
6363
include("preprocessing/preprocessing.jl")
6464
include("io/io.jl")
65+
include("general/restart.jl")
6566
include("visualization/recipes_plots.jl")
6667

6768
export Semidiscretization, semidiscretize, restart_with!

src/general/neighborhood_search.jl

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -261,6 +261,11 @@ end
261261
end
262262

263263
# === Initialization ===
264+
265+
function initialize_neighborhood_searches!(semi, u0_ode, restart_with::Nothing)
266+
initialize_neighborhood_searches!(semi)
267+
end
268+
264269
function initialize_neighborhood_searches!(semi)
265270
foreach_system(semi) do system
266271
foreach_system(semi) do neighbor

src/general/restart.jl

Lines changed: 135 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,135 @@
1+
# Update `v0_ode` and `u0_ode` with the data from the restart files
2+
function set_initial_conditions!(v0_ode, u0_ode, semi, restart_with::Tuple{Vararg{String}})
3+
# Check number of systems
4+
if length(semi.systems) != length(restart_with)
5+
throw(ArgumentError("number of systems in `semi` does not match number of restart files provided " *
6+
"in `restart_with`"))
7+
end
8+
9+
# Check that systems match
10+
expected_system_names = system_names(semi.systems)
11+
foreach_system(semi) do system
12+
system_index = system_indices(system, semi)
13+
filename = restart_with[system_index]
14+
expected_system_name = expected_system_names[system_index]
15+
if !occursin(expected_system_name, basename(splitext(filename)[1]))
16+
throw(ArgumentError("Filename '$filename' for system $system_index does not contain expected name '$expected_system_name'. " *
17+
"Expected a VTK file for system of type $(nameof(typeof(system)))."))
18+
end
19+
end
20+
21+
# Set initial conditions
22+
foreach_noalloc(semi.systems, restart_with) do (system, restart_file)
23+
v0_system = wrap_v(v0_ode, system, semi)
24+
u0_system = wrap_u(u0_ode, system, semi)
25+
26+
restart_data = vtk2trixi(restart_file; create_initial_condition=false)
27+
v_restart = restart_v(system, restart_data)
28+
u_restart = restart_u(system, restart_data)
29+
30+
v0_system .= Adapt.adapt(semi.parallelization_backend, v_restart)
31+
u0_system .= Adapt.adapt(semi.parallelization_backend, u_restart)
32+
33+
restore_previous_state!(system, restart_file)
34+
end
35+
end
36+
37+
# Compute a new `tspan` based on the restart files
38+
function time_span(tspan, restart_with::Tuple{Vararg{String}})
39+
# Read restart times from all files
40+
restart_times = [vtk2trixi(file).time for file in restart_with]
41+
t_restart = convert(eltype(tspan), first(restart_times))
42+
43+
# Check if all restart files have the same time
44+
if !all(isapprox(t, t_restart) for t in restart_times)
45+
throw(ArgumentError("all restart files must start from the same time"))
46+
end
47+
48+
if !isapprox(tspan[1], t_restart)
49+
@info "Adjusting initial time from $(tspan[1]) to restart time $t_restart"
50+
end
51+
52+
return (t_restart, tspan[2])
53+
end
54+
55+
function write_density_and_pressure!(v_restart, system, density_calculator,
56+
pressure, density)
57+
return v_restart
58+
end
59+
60+
function write_density_and_pressure!(v_restart, system,
61+
density_calculator::ContinuityDensity,
62+
pressure, density)
63+
v_restart[size(v_restart, 1), :] = density
64+
65+
return v_restart
66+
end
67+
68+
function write_density_and_pressure!(v_restart, system::EntropicallyDampedSPHSystem,
69+
density_calculator::ContinuityDensity,
70+
pressure, density)
71+
v_restart[size(v_restart, 1), :] = density
72+
v_restart[size(v_restart, 1) - 1, :] = pressure
73+
74+
return v_restart
75+
end
76+
77+
function write_density_and_pressure!(v_restart, system::EntropicallyDampedSPHSystem,
78+
density_calculator::SummationDensity,
79+
pressure, density)
80+
v_restart[size(v_restart, 1) - 1, :] = pressure
81+
82+
return v_restart
83+
end
84+
85+
restore_previous_state!(system, restart_file) = system
86+
87+
function initialize_neighborhood_searches!(semi, u0_ode,
88+
restart_with::Tuple{Vararg{String}})
89+
foreach_system(semi) do system
90+
foreach_system(semi) do neighbor
91+
initialize_neighborhood_search!(semi, system, neighbor, u0_ode)
92+
end
93+
end
94+
95+
return semi
96+
end
97+
98+
function initialize_neighborhood_search!(semi, system, neighbor, u0_ode)
99+
# TODO Initialize after adapting to the GPU.
100+
# Currently, this cannot use `semi.parallelization_backend`
101+
# because data is still on the CPU.
102+
PointNeighbors.initialize!(get_neighborhood_search(system, neighbor, semi),
103+
initial_restart_coordinates(system, u0_ode, semi),
104+
initial_restart_coordinates(neighbor, u0_ode, semi),
105+
eachindex_y=each_active_particle(neighbor),
106+
parallelization_backend=PolyesterBackend())
107+
return semi
108+
end
109+
110+
function initialize_neighborhood_search!(semi, system::TotalLagrangianSPHSystem,
111+
neighbor::TotalLagrangianSPHSystem, u0_ode)
112+
# For TLSPH, the self-interaction NHS is already initialized in the system constructor
113+
return semi
114+
end
115+
116+
function initial_restart_coordinates(system, u0_ode, semi)
117+
# Transfer to CPU if data is on the GPU. Do nothing if already on CPU.
118+
return transfer2cpu(wrap_u(u0_ode, system, semi))
119+
end
120+
121+
function initial_restart_coordinates(system::Union{WallBoundarySystem,
122+
AbstractStructureSystem}, u0_ode, semi)
123+
return initial_coordinates(system)
124+
end
125+
126+
function initialize!(semi::Semidiscretization, restart_with::Tuple{Vararg{String}})
127+
foreach_system(semi) do system
128+
# Initialize this system
129+
initialize_restart!(system, semi)
130+
end
131+
132+
return semi
133+
end
134+
135+
initialize_restart!(system, semi) = initialize!(system, semi)

src/general/semidiscretization.jl

Lines changed: 41 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -184,7 +184,7 @@ end
184184
end
185185

186186
"""
187-
semidiscretize(semi, tspan; reset_threads=true)
187+
semidiscretize(semi, tspan; reset_threads=true, restart_with=nothing)
188188
189189
Create an `ODEProblem` from the semidiscretization with the specified `tspan`.
190190
@@ -193,6 +193,14 @@ Create an `ODEProblem` from the semidiscretization with the specified `tspan`.
193193
- `tspan`: The time span over which the simulation will be run.
194194
195195
# Keywords
196+
- `restart_with`: Can be used to restart the simulation from VTK solution files (see [`SolutionSavingCallback`](@ref)).
197+
This can be either `nothing` (default, no restart) or a `Tuple` of filenames,
198+
one for each system in the [`Semidiscretization`](@ref).
199+
The order of the filenames must match the order of the systems in the [`Semidiscretization`](@ref).
200+
Note that `semidiscretize` replaces the initial time (`tspan[1]`) with the timestamp read
201+
from the VTK files. If the user-provided `tspan[1]` does not match the restart time,
202+
it is adjusted and an info message is logged. If multiple files are provided, their
203+
timestamps must match.
196204
- `reset_threads`: A boolean flag to reset Polyester.jl threads before the simulation (default: `true`).
197205
After an error within a threaded loop, threading might be disabled. Resetting the threads before the simulation
198206
ensures that threading is enabled again for the simulation.
@@ -219,9 +227,16 @@ timespan: (0.0, 1.0)
219227
u0: ([...], [...]) *this line is ignored by filter*
220228
```
221229
"""
222-
function semidiscretize(semi, tspan; reset_threads=true)
230+
function semidiscretize(semi, tspan; reset_threads=true, restart_with=nothing)
223231
(; systems) = semi
224232

233+
if restart_with isa String
234+
restart_with = (restart_with,)
235+
elseif !isnothing(restart_with) && !(restart_with isa NTuple{<:Any, String})
236+
throw(ArgumentError("`restart_with` must be `nothing`, a string, or a tuple of strings, " *
237+
"got $(typeof(restart_with))"))
238+
end
239+
225240
# Check that all systems have the same eltype
226241
first_system = first(systems)
227242
if !all(system -> eltype(system) === eltype(first_system), systems)
@@ -267,14 +282,11 @@ function semidiscretize(semi, tspan; reset_threads=true)
267282
end
268283

269284
# Set initial condition
270-
foreach_system_wrapped(semi, v0_ode, u0_ode) do system, v0_system, u0_system
271-
write_u0!(u0_system, system)
272-
write_v0!(v0_system, system)
273-
end
285+
set_initial_conditions!(v0_ode, u0_ode, semi, restart_with)
274286

275287
# TODO initialize after adapting to the GPU.
276288
# Requires https://github.com/trixi-framework/PointNeighbors.jl/pull/86.
277-
initialize_neighborhood_searches!(semi)
289+
initialize_neighborhood_searches!(semi, u0_ode, restart_with)
278290

279291
if semi.parallelization_backend isa KernelAbstractions.GPU
280292
# Convert all arrays in the systems to the correct array type.
@@ -304,10 +316,7 @@ function semidiscretize(semi, tspan; reset_threads=true)
304316
end
305317

306318
# Initialize all particle systems
307-
foreach_system(semi_new) do system
308-
# Initialize this system
309-
initialize!(system, semi_new)
310-
end
319+
initialize!(semi_new, restart_with)
311320

312321
# Reset callback flag that will be set by the `UpdateCallback`
313322
semi_new.update_callback_used[] = false
@@ -318,7 +327,27 @@ function semidiscretize(semi, tspan; reset_threads=true)
318327
# here, since we cannot change `p` from within the callback (only its contents).
319328
p = @NamedTuple{semi::typeof(semi_new), split_integration_data::Any}((semi_new,
320329
nothing))
321-
return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode, tspan, p)
330+
331+
return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode,
332+
time_span(tspan, restart_with), p)
333+
end
334+
335+
function set_initial_conditions!(v0_ode, u0_ode, semi, restart_with::Nothing)
336+
foreach_system_wrapped(semi, v0_ode, u0_ode) do system, v0_system, u0_system
337+
write_u0!(u0_system, system)
338+
write_v0!(v0_system, system)
339+
end
340+
end
341+
342+
time_span(tspan, restart_with::Nothing) = tspan
343+
344+
function initialize!(semi::Semidiscretization, restart_with::Nothing)
345+
foreach_system(semi) do system
346+
# Initialize this system
347+
initialize!(system, semi)
348+
end
349+
350+
return semi
322351
end
323352

324353
"""

src/io/read_vtk.jl

Lines changed: 15 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -51,8 +51,19 @@ function vtk2trixi(file; element_type=nothing, coordinates_eltype=nothing,
5151
point_coords = ReadVTK.get_points(vtk_file)
5252

5353
cELTYPE = isnothing(coordinates_eltype) ? eltype(point_coords) : coordinates_eltype
54-
ELTYPE = isnothing(element_type) ?
55-
eltype(first(ReadVTK.get_data(point_data["pressure"]))) : element_type
54+
55+
if !isnothing(element_type)
56+
ELTYPE = element_type
57+
else
58+
# Try to get element type from pressure or density (whichever exists)
59+
ELTYPE = if "pressure" in keys(point_data)
60+
eltype(first(ReadVTK.get_data(point_data["pressure"])))
61+
elseif "material_density" in keys(point_data)
62+
eltype(first(ReadVTK.get_data(point_data["material_density"])))
63+
else
64+
error("neither 'pressure' nor 'material_density' field found in VTK file")
65+
end
66+
end
5667

5768
results = Dict{Symbol, Any}()
5869

@@ -87,7 +98,8 @@ function vtk2trixi(file; element_type=nothing, coordinates_eltype=nothing,
8798
results[:particle_spacing]
8899
results[:coordinates] = coordinates
89100
results[:time] = "time" in keys(field_data) ?
90-
first(ReadVTK.get_data(field_data["time"])) : zero(ELTYPE)
101+
convert(ELTYPE, first(ReadVTK.get_data(field_data["time"]))) :
102+
zero(ELTYPE)
91103

92104
append!(used_keys, ["index", "ndims"])
93105
# Load any custom quantities

src/io/write_vtk.jl

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -447,6 +447,16 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySystem)
447447
for particle in eachparticle(system)]
448448
vtk["pressure"] = [current_pressure(v, system, particle)
449449
for particle in eachparticle(system)]
450+
vtk["zone_id"] = [system.boundary_zone_indices[particle]
451+
for particle in eachparticle(system)]
452+
453+
if any(pm -> isa(pm, AbstractPressureModel), system.cache.pressure_reference_values)
454+
for (i, pressure_model) in enumerate(system.cache.pressure_reference_values)
455+
if pressure_model isa AbstractPressureModel
456+
vtk["boundary_zone_pressure_$i"] = system.cache.pressure_reference_values[i].pressure[]
457+
end
458+
end
459+
end
450460

451461
if system.calculate_flow_rate
452462
Q_total = zero(eltype(system))

0 commit comments

Comments
 (0)