Restart from VTU

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.condition.index[] - 1
+
+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

@@ -61,6 +61,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/restart.jl

@@ -0,0 +1,119 @@
+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
+
+function time_span(tspan, restart_with::Tuple{Vararg{String}})
+    restart_data = vtk2trixi(first(restart_with))
+    t_restart = convert(eltype(tspan), restart_data.time)
+
+    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
+
+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)

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,10 @@ 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 has to be a tuple of filenames, one for each system in the [`Semidiscretization`](@ref).
+  The order of the filenames has to match the order of the systems in the [`Semidiscretization`](@ref).
+  If no restart is desired, use `nothing` (default).
 - `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,7 +223,7 @@ 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
 
     # Check that all systems have the same eltype
@@ -267,14 +271,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,15 +305,31 @@ 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
 
-    return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode, tspan, semi_new)
+    return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode,
+                               time_span(tspan, restart_with), semi_new)
+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
 
 """
@@ -348,6 +365,39 @@ function restart_with!(semi, sol; reset_threads=true)
     return semi
 end
 
+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
+            initialize_neighborhood_search!(semi, system, neighbor)
+        end
+    end
+
+    return semi
+end
+
+function initialize_neighborhood_search!(semi, system, neighbor)
+    # 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_coordinates(system),
+                               initial_coordinates(neighbor),
+                               eachindex_y=each_active_particle(neighbor),
+                               parallelization_backend=PolyesterBackend())
+
+    return semi
+end
+
+function initialize_neighborhood_search!(semi, system::TotalLagrangianSPHSystem,
+                                         neighbor::TotalLagrangianSPHSystem)
+    # For TLSPH, the self-interaction NHS is already initialized in the system constructor
+    return semi
+end
+
 # We have to pass `system` here for type stability,
 # since the type of `system` determines the return type.
 @inline function wrap_v(v_ode, system, semi)

src/io/write_vtk.jl

@@ -458,6 +458,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))

src/schemes/boundary/open_boundary/system.jl

@@ -133,6 +133,9 @@ function initialize!(system::OpenBoundarySystem, semi)
     return system
 end
 
+# Skip during restart, as boundary zone indices are updated in `restore_previous_state!`
+initialize_restart!(system::OpenBoundarySystem, semi) = system
+
 function create_cache_open_boundary(boundary_model, fluid_system, initial_condition,
                                     density_diffusion, calculate_flow_rate, boundary_zones)
     reference_values = map(bz -> bz.reference_values, boundary_zones)
@@ -251,6 +254,8 @@ end
     return system.shifting_technique
 end
 
+@inline density_calculator(system::OpenBoundarySystem) = nothing
+
 system_sound_speed(system::OpenBoundarySystem) = system_sound_speed(system.fluid_system)
 
 @inline hydrodynamic_mass(system::OpenBoundarySystem, particle) = system.mass[particle]
@@ -680,6 +685,66 @@ function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone,
     return system
 end
 
+function restart_u(system::OpenBoundarySystem, data)
+    coords_total = zeros(coordinates_eltype(system), u_nvariables(system),
+                         n_integrated_particles(system))
+    coords_total .= coordinates_eltype(system)(1e16)
+
+    coords_active = data.coordinates
+    for particle in axes(coords_active, 2)
+        for dim in 1:ndims(system)
+            coords_total[dim, particle] = coords_active[dim, particle]
+        end
+    end
+
+    system.buffer.active_particle .= false
+    system.buffer.active_particle[1:size(coords_active, 2)] .= true
+
+    update_system_buffer!(system.buffer)
+
+    return coords_total
+end
+
+function restart_v(system::OpenBoundarySystem, data)
+    v_total = zeros(eltype(system), v_nvariables(system),
+                    n_integrated_particles(system))
+
+    v_active = zeros(eltype(system), v_nvariables(system), size(data.velocity, 2))
+
+    v_active[1:ndims(system), :] = data.velocity
+    write_density_and_pressure!(v_active, system.fluid_system,
+                                density_calculator(system), data.pressure, data.density)
+
+    for particle in axes(v_active, 2)
+        for i in axes(v_active, 1)
+            v_total[i, particle] = v_active[i, particle]
+        end
+    end
+
+    return v_total
+end
+
+function restore_previous_state!(system::OpenBoundarySystem, file)
+    # We cannot simply use `update_boundary_zone_indices!` because rounding errors during file I/O
+    # may result in particles being located outside their intended boundary zone, even though they
+    # were written as active particles.
+    set_zero!(system.boundary_zone_indices)
+
+    values = vtk2trixi(file; create_initial_condition=false)
+    system.boundary_zone_indices[each_integrated_particle(system)] .= values.zone_id
+
+    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
+                pressure_model.pressure[] = values[Symbol(:boundary_zone_pressure_, i)]
+                pressure_model.flow_rate[] = values[Symbol(:Q_, i)]
+            end
+        end
+    end
+
+    return system
+end
+
 # Open-boundary interpolation should reconstruct the surrounding fluid-like velocity field.
 # Therefore, only actual fluid systems and other open-boundary particles contribute;
 # rigid bodies, walls, and other non-fluid systems are intentionally excluded.

src/schemes/boundary/wall_boundary/dummy_particles.jl

@@ -328,6 +328,10 @@ function Base.show(io::IO, model::BoundaryModelDummyParticles)
     print(io, ")")
 end
 
+@inline function density_calculator(model::BoundaryModelDummyParticles)
+    return model.density_calculator
+end
+
 # Set the initial pressure to zero for visualization
 function initial_boundary_pressure(initial_density, density_calculator, _)
     return zero(initial_density)