Fix Issue #263

examples/n_body/n_body_solar_system.jl

@@ -39,7 +39,7 @@ tspan = (0.0, 10year)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=100000)
-saving_callback = SolutionSavingCallback(dt=10day)
+saving_callback = SolutionSavingCallback(dt=10day, max_coordinates=Inf)
 
 callbacks = CallbackSet(info_callback, saving_callback)
 

src/callbacks/solution_saving.jl

@@ -2,7 +2,7 @@
     SolutionSavingCallback(; interval::Integer=0, dt=0.0,
                            save_initial_solution=true,
                            save_final_solution=true,
-                           output_directory="out", append_timestamp=false,
+                           output_directory="out", append_timestamp=false, max_coordinates=2^15,
                            custom_quantities...)
 
 Callback to save the current numerical solution in VTK format in regular intervals.
@@ -28,6 +28,7 @@ To ignore a custom quantity for a specific system, return `nothing`.
 - `custom_quantities...`:       Additional user-defined quantities.
 - `write_meta_data`:            Write meta data.
 - `verbose=false`:              Print to standard IO when a file is written.
+- `max_coordinates=2^15`        Particles with absolute coordinates larger than this will be clipped.
 
 # Examples
 ```julia
@@ -57,6 +58,7 @@ struct SolutionSavingCallback{I, CQ}
     verbose               :: Bool
     output_directory      :: String
     prefix                :: String
+    max_coordinates       :: Float64
     custom_quantities     :: CQ
     latest_saved_iter     :: Vector{Int}
 end
@@ -66,7 +68,7 @@ function SolutionSavingCallback(; interval::Integer=0, dt=0.0,
                                 save_final_solution=true,
                                 output_directory="out", append_timestamp=false,
                                 prefix="", verbose=false, write_meta_data=true,
-                                custom_quantities...)
+                                max_coordinates=Float64(2^15), custom_quantities...)
     if dt > 0 && interval > 0
         throw(ArgumentError("Setting both interval and dt is not supported!"))
     end
@@ -82,7 +84,8 @@ function SolutionSavingCallback(; interval::Integer=0, dt=0.0,
     solution_callback = SolutionSavingCallback(interval,
                                                save_initial_solution, save_final_solution,
                                                write_meta_data, verbose, output_directory,
-                                               prefix, custom_quantities, [-1])
+                                               prefix, max_coordinates, custom_quantities,
+                                               [-1])
 
     if dt > 0
         # Add a `tstop` every `dt`, and save the final solution.
@@ -137,7 +140,7 @@ end
 # affect!
 function (solution_callback::SolutionSavingCallback)(integrator)
     (; interval, output_directory, custom_quantities, write_meta_data,
-    verbose, prefix, latest_saved_iter) = solution_callback
+    verbose, prefix, latest_saved_iter, max_coordinates) = solution_callback
 
     vu_ode = integrator.u
     semi = integrator.p
@@ -162,6 +165,7 @@ function (solution_callback::SolutionSavingCallback)(integrator)
                                                     output_directory=output_directory,
                                                     prefix=prefix,
                                                     write_meta_data=write_meta_data,
+                                                    max_coordinates=max_coordinates,
                                                     custom_quantities...)
 
     # Tell OrdinaryDiffEq that u has not been modified

src/visualization/write2vtk.jl

@@ -15,6 +15,7 @@ Convert Trixi simulation data to VTK format.
 - `prefix`:               Prefix for output files. Defaults to an empty string.
 - `write_meta_data`:      Write meta data.
 - `custom_quantities...`: Additional custom quantities to include in the VTK output. TODO.
+- `max_coordinates=Inf`   Particles with absolute coordinates larger than this will be clipped.
 
 
 # Example
@@ -24,7 +25,7 @@ trixi2vtk(sol.u[end], semi, 0.0, iter=1, output_directory="output", prefix="solu
 TODO: example for custom_quantities
 """
 function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix="",
-                   write_meta_data=true, custom_quantities...)
+                   write_meta_data=true, max_coordinates=Inf, custom_quantities...)
     (; systems) = semi
     v_ode, u_ode = vu_ode.x
 
@@ -45,13 +46,15 @@ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix
         trixi2vtk(v, u, t, system, periodic_box;
                   output_directory=output_directory,
                   system_name=filenames[system_index], iter=iter, prefix=prefix,
-                  write_meta_data=write_meta_data, custom_quantities...)
+                  write_meta_data=write_meta_data, max_coordinates=max_coordinates,
+                  custom_quantities...)
     end
 end
 
 # Convert data for a single TrixiParticle system to VTK format
 function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix="",
                    iter=nothing, system_name=vtkname(system), write_meta_data=true,
+                   max_coordinates=Inf,
                    custom_quantities...)
     mkpath(output_directory)
 
@@ -71,6 +74,13 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix
     points = periodic_coords(current_coordinates(u, system), periodic_box)
     cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in axes(points, 2)]
 
+    if abs(maximum(points)) > max_coordinates || abs(minimum(points)) > max_coordinates
+        println("Warning: At least one particle's absolute coordinates exceed `max_coordinates` and have been clipped")
+        for i in eachindex(points)
+            points[i] = clamp(points[i], -max_coordinates, max_coordinates)
+        end
+    end
+
     vtk_grid(file, points, cells) do vtk
         write2vtk!(vtk, v, u, t, system, write_meta_data=write_meta_data)
 
@@ -107,10 +117,12 @@ end
 Convert coordinate data to VTK format.
 
 # Arguments
-- `coordinates`:                 Coordinates to be saved.
-- `output_directory` (optional): Output directory path. Defaults to `"out"`.
-- `prefix` (optional):           Prefix for the output file. Defaults to an empty string.
-- `filename` (optional):         Name of the output file. Defaults to `"coordinates"`.
+- `coordinates`: Coordinates to be saved.
+
+# Keywords
+- `output_directory="out"`: Output directory path.
+- `prefix=""`:              Prefix for the output file.
+- `filename="coordinates"`: Name of the output file.
 
 # Returns
 - `file::AbstractString`: Path to the generated VTK file.