Restart from VTU

NEWS.md

@@ -4,12 +4,12 @@ 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.4.4
+## Version 0.5
 
-### Features
-
-- Added `StateEquationAdaptiveCole` an adaptive sound speed version of the Cole state equation (#875)
+### API Changes
 
+- Return type of `vtk2trixi` changed to `NamedTuple` including an optional
+  `:initial_condition` field if `create_initial_condition=true` is passed. (#959)
 
 ## Version 0.4.3
 
@@ -34,7 +34,7 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 
 - Added GPU and FP32 support for DEM (#979).
 
-  
+
 ### Performance
 - Improved GPU performance with shifting up to a factor of 10x (#974, #993).
 
@@ -360,5 +360,3 @@ Features:
 #### TLSPH
 
 An implementation of TLSPH (Total Lagrangian Smoothed Particle Hydrodynamics) for solid bodies enabling FSI (Fluid Structure Interactions).
-
-

examples/postprocessing/restart_poiseuille_flow_2d.jl

@@ -0,0 +1,34 @@
+# ==========================================================================================
+# Restart Example
+#
+# TODO
+# ==========================================================================================
+using TrixiParticles
+
+# Run full simulation
+# trixi_include(@__MODULE__,
+#               joinpath(examples_dir(), "fluid", "poiseuille_flow_2d.jl"),
+#               tspan=(0.0, 0.6), sound_speed_factor=10, particle_spacing=4e-5)
+
+# Run half simulation
+trixi_include(@__MODULE__,
+              joinpath(examples_dir(), "fluid", "poiseuille_flow_2d.jl"),
+              tspan=(0.0, 0.3), sound_speed_factor=10, particle_spacing=4e-5)
+
+iter = round(Int, 0.3 / 0.02)
+
+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

@@ -59,6 +59,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/buffer.jl

@@ -40,10 +40,10 @@ end
 # Dispatch by system type to handle systems that provide a buffer.
 @inline buffer(system) = nothing
 
-@inline update_system_buffer!(buffer::Nothing, semi) = buffer
+@inline update_system_buffer!(buffer::Nothing) = buffer
 
 # TODO `resize` allocates. Find a non-allocating version
-@inline function update_system_buffer!(buffer::SystemBuffer, semi)
+@inline function update_system_buffer!(buffer::SystemBuffer)
     (; active_particle) = buffer
 
     # TODO: Parallelize (see https://github.com/trixi-framework/TrixiParticles.jl/issues/810)

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)
+
+        precondition_system!(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
+
+precondition_system!(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

@@ -302,7 +302,7 @@ end
 @inline foreach_system(f, systems) = foreach_noalloc(f, systems)
 
 """
-    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`.
 
@@ -311,6 +311,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.
@@ -337,7 +341,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
@@ -378,17 +382,11 @@ function semidiscretize(semi, tspan; reset_threads=true)
     end
 
     # Set initial condition
-    foreach_system(semi) do system
-        u0_system = wrap_u(u0_ode, system, semi)
-        v0_system = wrap_v(v0_ode, system, semi)
-
-        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 to the correct array type.
@@ -415,15 +413,34 @@ 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(semi) do system
+        v0_system = wrap_v(v0_ode, system, semi)
+        u0_system = wrap_u(u0_ode, system, semi)
+
+        write_v0!(v0_system, system)
+        write_u0!(u0_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
 
 """
@@ -461,6 +478,10 @@ 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
@@ -1186,7 +1207,7 @@ function check_configuration(system::ImplicitIncompressibleSPHSystem, systems, n
                 neighbor.boundary_model.density_calculator isa PressureBoundaries)
                 time_step_boundary = neighbor.boundary_model.density_calculator.time_step
                 omega_boundary = neighbor.boundary_model.density_calculator.omega
-                if !(time_step==time_step_boundary && omega==omega_boundary)
+                if !(time_step == time_step_boundary && omega == omega_boundary)
                     throw(ArgumentError("`PressureBoundaries` parameters have to be the same as the
                     `ImplicitIncompressibleSPHSystem` parameters"))
                 end