Preserve time stepping information in restart

NEWS.md

@@ -4,9 +4,16 @@ TrixiParticles.jl follows the interpretation of
 [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
 
+## Version 0.5.2
+
+### Features
+- Added the ability to restart simulations from VTK solution files generated by `SolutionSavingCallback`.
+  Users can now pass the `restart_with` keyword argument to `semidiscretize` (#1190).
+
 ## Version 0.5.1
 
 ### Features
+
 - Implemented stage-level coupling for split integration (#1049).
 - Added `MechanicalWorkCalculatorCallback` (#940).
 

examples/postprocessing/restart_poiseuille_flow_2d.jl

@@ -0,0 +1,31 @@
+# ==========================================================================================
+# Restart Example: Poiseuille Flow 2D
+#
+# This example demonstrates how to restart a simulation.
+# We first run a simulation of 2D Poiseuille flow up to t=0.3s, then restart from the
+# saved state and continue the simulation until t=0.6s.
+# ==========================================================================================
+using TrixiParticles
+
+trixi_include(@__MODULE__,
+              joinpath(examples_dir(), "fluid", "poiseuille_flow_2d.jl"),
+              tspan=(0.0, 0.3), sound_speed_factor=10, particle_spacing=4e-5)
+
+# Get latest iteration
+iter = saving_callback.affect!.affect!.latest_saved_iter
+
+restart_file_fluid = joinpath("out", "fluid_1_$iter.vtu")
+restart_file_open_boundary = joinpath("out", "open_boundary_1_$iter.vtu")
+restart_file_boundary = joinpath("out", "boundary_1_$iter.vtu")
+
+ode_restart = semidiscretize(semi, (0.3, 0.6);
+                             restart_with=(restart_file_fluid,
+                                           restart_file_open_boundary,
+                                           restart_file_boundary))
+
+saving_callback = SolutionSavingCallback(dt=0.02, prefix="restart")
+
+callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback())
+
+sol_restart = solve(ode_restart, RDPK3SpFSAL35(), abstol=1e-5, reltol=1e-3, dtmax=1e-2,
+                    save_everystep=false, callback=callbacks)

src/TrixiParticles.jl

@@ -62,6 +62,7 @@ include("general/semidiscretization.jl")
 include("general/gpu.jl")
 include("preprocessing/preprocessing.jl")
 include("io/io.jl")
+include("general/restart.jl")
 include("visualization/recipes_plots.jl")
 
 export Semidiscretization, semidiscretize, restart_with!

src/general/neighborhood_search.jl

@@ -261,6 +261,11 @@ end
 end
 
 # === Initialization ===
+
+function initialize_neighborhood_searches!(semi, u0_ode, restart_with::Nothing)
+    initialize_neighborhood_searches!(semi)
+end
+
 function initialize_neighborhood_searches!(semi)
     foreach_system(semi) do system
         foreach_system(semi) do neighbor

src/general/restart.jl

@@ -0,0 +1,139 @@
+# Update `v0_ode` and `u0_ode` with the data from the restart files
+function set_initial_conditions!(v0_ode, u0_ode, semi, restart_with::Tuple{Vararg{String}})
+    # Check number of systems
+    if length(semi.systems) != length(restart_with)
+        throw(ArgumentError("number of systems in `semi` does not match number of restart files provided " *
+                            "in `restart_with`"))
+    end
+
+    # Check that systems match
+    expected_system_names = system_names(semi.systems)
+    foreach_system(semi) do system
+        system_index = system_indices(system, semi)
+        filename = restart_with[system_index]
+        expected_system_name = expected_system_names[system_index]
+        if !occursin(expected_system_name, basename(splitext(filename)[1]))
+            throw(ArgumentError("Filename '$filename' for system $system_index does not contain expected name '$expected_system_name'. " *
+                                "Expected a VTK file for system of type $(nameof(typeof(system)))."))
+        end
+    end
+
+    # Set initial conditions
+    foreach_noalloc(semi.systems, restart_with) do (system, restart_file)
+        v0_system = wrap_v(v0_ode, system, semi)
+        u0_system = wrap_u(u0_ode, system, semi)
+
+        restart_data = vtk2trixi(restart_file; create_initial_condition=false)
+        v_restart = restart_v(system, restart_data)
+        u_restart = restart_u(system, restart_data)
+
+        v0_system .= Adapt.adapt(semi.parallelization_backend, v_restart)
+        u0_system .= Adapt.adapt(semi.parallelization_backend, u_restart)
+
+        restore_previous_state!(system, restart_file)
+    end
+end
+
+# Compute a new `tspan` based on the restart files
+function time_span(tspan, restart_with::Tuple{Vararg{String}})
+    # Read restart times from all files
+    restart_times = [vtk2trixi(file).time for file in restart_with]
+    t_restart = convert(eltype(tspan), first(restart_times))
+
+    # Check if all restart files have the same time
+    if !all(isapprox(t, t_restart) for t in restart_times)
+        throw(ArgumentError("all restart files must start from the same time"))
+    end
+
+    if !isapprox(tspan[1], t_restart)
+        @info "Adjusting initial time from $(tspan[1]) to restart time $t_restart"
+    end
+
+    return (t_restart, tspan[2])
+end
+
+function write_density_and_pressure!(v_restart, system, density_calculator,
+                                     pressure, density)
+    return v_restart
+end
+
+function write_density_and_pressure!(v_restart, system,
+                                     density_calculator::ContinuityDensity,
+                                     pressure, density)
+    v_restart[size(v_restart, 1), :] = density
+
+    return v_restart
+end
+
+function write_density_and_pressure!(v_restart, system::EntropicallyDampedSPHSystem,
+                                     density_calculator::ContinuityDensity,
+                                     pressure, density)
+    v_restart[size(v_restart, 1), :] = density
+    v_restart[size(v_restart, 1) - 1, :] = pressure
+
+    return v_restart
+end
+
+function write_density_and_pressure!(v_restart, system::EntropicallyDampedSPHSystem,
+                                     density_calculator::SummationDensity,
+                                     pressure, density)
+    v_restart[size(v_restart, 1) - 1, :] = pressure
+
+    return v_restart
+end
+
+restore_previous_state!(system, restart_file) = system
+
+function initialize_neighborhood_searches!(semi, u0_ode,
+                                           restart_with::Tuple{Vararg{String}})
+    foreach_system(semi) do system
+        foreach_system(semi) do neighbor
+            initialize_neighborhood_search!(semi, system, neighbor, u0_ode)
+        end
+    end
+
+    return semi
+end
+
+function initialize_neighborhood_search!(semi, system, neighbor, u0_ode)
+    # TODO Initialize after adapting to the GPU.
+    # Currently, this cannot use `semi.parallelization_backend`
+    # because data is still on the CPU.
+    PointNeighbors.initialize!(get_neighborhood_search(system, neighbor, semi),
+                               initial_restart_coordinates(system, u0_ode, semi),
+                               initial_restart_coordinates(neighbor, u0_ode, semi),
+                               eachindex_y=each_active_particle(neighbor),
+                               parallelization_backend=PolyesterBackend())
+    return semi
+end
+
+function initialize_neighborhood_search!(semi, system::TotalLagrangianSPHSystem,
+                                         neighbor::TotalLagrangianSPHSystem, u0_ode)
+    # For TLSPH, the self-interaction NHS is already initialized in the system constructor
+    return semi
+end
+
+function initial_restart_coordinates(system, u0_ode, semi)
+    # Transfer to CPU if data is on the GPU. Do nothing if already on CPU.
+    return transfer2cpu(wrap_u(u0_ode, system, semi))
+end
+
+function initial_restart_coordinates(system::Union{WallBoundarySystem,
+                                                   AbstractStructureSystem}, u0_ode, semi)
+    return initial_coordinates(system)
+end
+
+function initialize!(semi::Semidiscretization, restart_with::Tuple{Vararg{String}})
+    foreach_system(semi) do system
+        # Initialize this system
+        initialize_restart!(system, semi)
+    end
+
+    return semi
+end
+
+initialize_restart!(system, semi) = initialize!(system, semi)
+
+
+# TODO
+# Store controller values for an adaptive time stepping scheme

src/general/semidiscretization.jl

@@ -184,7 +184,7 @@ end
 end
 
 """
-    semidiscretize(semi, tspan; reset_threads=true)
+    semidiscretize(semi, tspan; reset_threads=true, restart_with=nothing)
 
 Create an `ODEProblem` from the semidiscretization with the specified `tspan`.
 
@@ -193,6 +193,14 @@ Create an `ODEProblem` from the semidiscretization with the specified `tspan`.
 - `tspan`: The time span over which the simulation will be run.
 
 # Keywords
+- `restart_with`: Can be used to restart the simulation from VTK solution files (see [`SolutionSavingCallback`](@ref)).
+  This can be either `nothing` (default, no restart) or a `Tuple` of filenames,
+  one for each system in the [`Semidiscretization`](@ref).
+  The order of the filenames must match the order of the systems in the [`Semidiscretization`](@ref).
+  Note that `semidiscretize` replaces the initial time (`tspan[1]`) with the timestamp read
+  from the VTK files. If the user-provided `tspan[1]` does not match the restart time,
+  it is adjusted and an info message is logged. If multiple files are provided, their
+  timestamps must match.
 - `reset_threads`: A boolean flag to reset Polyester.jl threads before the simulation (default: `true`).
   After an error within a threaded loop, threading might be disabled. Resetting the threads before the simulation
   ensures that threading is enabled again for the simulation.
@@ -219,9 +227,16 @@ timespan: (0.0, 1.0)
 u0: ([...], [...]) *this line is ignored by filter*
 ```
 """
-function semidiscretize(semi, tspan; reset_threads=true)
+function semidiscretize(semi, tspan; reset_threads=true, restart_with=nothing)
     (; systems) = semi
 
+    if restart_with isa String
+        restart_with = (restart_with,)
+    elseif !isnothing(restart_with) && !(restart_with isa NTuple{<:Any, String})
+        throw(ArgumentError("`restart_with` must be `nothing`, a string, or a tuple of strings, " *
+                            "got $(typeof(restart_with))"))
+    end
+
     # Check that all systems have the same eltype
     first_system = first(systems)
     if !all(system -> eltype(system) === eltype(first_system), systems)
@@ -267,14 +282,11 @@ function semidiscretize(semi, tspan; reset_threads=true)
     end
 
     # Set initial condition
-    foreach_system_wrapped(semi, v0_ode, u0_ode) do system, v0_system, u0_system
-        write_u0!(u0_system, system)
-        write_v0!(v0_system, system)
-    end
+    set_initial_conditions!(v0_ode, u0_ode, semi, restart_with)
 
     # TODO initialize after adapting to the GPU.
     # Requires https://github.com/trixi-framework/PointNeighbors.jl/pull/86.
-    initialize_neighborhood_searches!(semi)
+    initialize_neighborhood_searches!(semi, u0_ode, restart_with)
 
     if semi.parallelization_backend isa KernelAbstractions.GPU
         # Convert all arrays in the systems to the correct array type.
@@ -304,10 +316,7 @@ function semidiscretize(semi, tspan; reset_threads=true)
     end
 
     # Initialize all particle systems
-    foreach_system(semi_new) do system
-        # Initialize this system
-        initialize!(system, semi_new)
-    end
+    initialize!(semi_new, restart_with)
 
     # Reset callback flag that will be set by the `UpdateCallback`
     semi_new.update_callback_used[] = false
@@ -318,7 +327,27 @@ function semidiscretize(semi, tspan; reset_threads=true)
     # here, since we cannot change `p` from within the callback (only its contents).
     p = @NamedTuple{semi::typeof(semi_new), split_integration_data::Any}((semi_new,
                                                                           nothing))
-    return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode, tspan, p)
+
+    return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode,
+                               time_span(tspan, restart_with), p)
+end
+
+function set_initial_conditions!(v0_ode, u0_ode, semi, restart_with::Nothing)
+    foreach_system_wrapped(semi, v0_ode, u0_ode) do system, v0_system, u0_system
+        write_u0!(u0_system, system)
+        write_v0!(v0_system, system)
+    end
+end
+
+time_span(tspan, restart_with::Nothing) = tspan
+
+function initialize!(semi::Semidiscretization, restart_with::Nothing)
+    foreach_system(semi) do system
+        # Initialize this system
+        initialize!(system, semi)
+    end
+
+    return semi
 end
 
 """

src/io/read_vtk.jl

@@ -51,8 +51,19 @@ function vtk2trixi(file; element_type=nothing, coordinates_eltype=nothing,
     point_coords = ReadVTK.get_points(vtk_file)
 
     cELTYPE = isnothing(coordinates_eltype) ? eltype(point_coords) : coordinates_eltype
-    ELTYPE = isnothing(element_type) ?
-             eltype(first(ReadVTK.get_data(point_data["pressure"]))) : element_type
+
+    if !isnothing(element_type)
+        ELTYPE = element_type
+    else
+        # Try to get element type from pressure or density (whichever exists)
+        ELTYPE = if "pressure" in keys(point_data)
+            eltype(first(ReadVTK.get_data(point_data["pressure"])))
+        elseif "material_density" in keys(point_data)
+            eltype(first(ReadVTK.get_data(point_data["material_density"])))
+        else
+            error("neither 'pressure' nor 'material_density' field found in VTK file")
+        end
+    end
 
     results = Dict{Symbol, Any}()
 
@@ -87,7 +98,8 @@ function vtk2trixi(file; element_type=nothing, coordinates_eltype=nothing,
                                  results[:particle_spacing]
     results[:coordinates] = coordinates
     results[:time] = "time" in keys(field_data) ?
-                     first(ReadVTK.get_data(field_data["time"])) : zero(ELTYPE)
+                     convert(ELTYPE, first(ReadVTK.get_data(field_data["time"]))) :
+                     zero(ELTYPE)
 
     append!(used_keys, ["index", "ndims"])
     # Load any custom quantities

src/io/write_vtk.jl

@@ -447,6 +447,16 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySystem)
                       for particle in eachparticle(system)]
     vtk["pressure"] = [current_pressure(v, system, particle)
                        for particle in eachparticle(system)]
+    vtk["zone_id"] = [system.boundary_zone_indices[particle]
+                      for particle in eachparticle(system)]
+
+    if any(pm -> isa(pm, AbstractPressureModel), system.cache.pressure_reference_values)
+        for (i, pressure_model) in enumerate(system.cache.pressure_reference_values)
+            if pressure_model isa AbstractPressureModel
+                vtk["boundary_zone_pressure_$i"] = system.cache.pressure_reference_values[i].pressure[]
+            end
+        end
+    end
 
     if system.calculate_flow_rate
         Q_total = zero(eltype(system))