Enable easy restart (#87)

* Start writing restart function

* Finish restart framework

* Update restart sim

* Improve re-start documentation

* Output writer constructor comments added

* Give start_tstep default value

* Add documentation that ghost floes are saved

documentation.md

@@ -569,23 +569,25 @@ Note that other than `∆tout`, all values have default values so you only need
 
 #### FloeOutputWriter
 The floe output writer allows you to save floe values into a JLD2 file. This will give you the ability to easily analyze floe fields across timesteps. A `FloeOutputWriter` has four field:
-- `outputs`, which specifies which floe fields should be included
 - `Δtout`, which specifies the umber of timesteps between floe outputs starting from the first timestep
+- `outputs`, which specifies which floe fields should be included
 - `filename`, including the path to the file, to save the file to
 - `overwrite` boolean that specifies whether a file with the same filename should be overwritter or if an error should be thrown during creation.
 
 You can create an `FloeOutputWriter` directly as a struct or with the following function call:
 ```julia
 floewriter = FloeOutputWriter(
- outputs,
- Δtout,
+ Δtout;
+ outputs = collect(fieldnames(Floe)),
  dir = ".",
  filename = "floes.jld2",
  overwrite = false,
  jld2_kw = Dict{Symbol, Any}(),
 )
 ```
-The `outputs` field takes in a list of symbols corresponding to floe fields. For example, if you want the floe output writer to output the floes centroid and coordinates then `outputs = [:centroid, :coords]`. If you want all floe fields then you can simply admit the outputs field altogether and all floe fields will be output. Note that other than `outputs` and `∆tout`, all values have default values so you only need to pass in arguments that you wish to change from these values. Furthermore, all inputs but ∆tout are keyword arguments, so you must use the keyword when passing in new values.
+The `outputs` field takes in a list of symbols corresponding to floe fields. For example, if you want the floe output writer to output the floes centroid and coordinates then `outputs = [:centroid, :coords]`. If you want all floe fields then you can simply omit the outputs field all together and all floe fields will be output. Note that other than `∆tout`, all values have default values so you only need to pass in arguments that you wish to change from these values. Furthermore, all inputs but ∆tout are keyword arguments, so you must use the keyword when passing in new values.
+
+Note that if you have Periodic calls, and thus ghost floes in your simulation, these will also be saved by the `FloeOutputWriter`. If you want to exclude these floes from your analysis or when otherwise using the `FloeOutputWriter` output, you can do so by only including floes with a `ghost_id = 0`.
 
 #### GridOutputWriter
 The grid output writer allows you to floe values averaged onto a course grid to a NetCDF file. This will give you the ability to easily analyze floe characteristics on a grid. A `GridOutputWriter` has eight field:
@@ -602,10 +604,10 @@ You can create an `GridOutputWriter` directly as a struct or with the following
 ```julia
 gridwriter = GridOutputWriter(
  FT,
- outputs,
  Δtout,
  grid,
  dims;
+ outputs = collect(get_known_grid_outputs()),
  dir = ".",
  filename = "gridded_data.nc",
  overwrite = false,
@@ -614,7 +616,7 @@ gridwriter = GridOutputWriter(
 ```
 The `outputs` field takes in a list of symbols. To see all possible outputs, call the `get_known_grid_outputs()` function. For example, if you want the grid output writer to output the floe masses and areas averaged on the grid then `outputs = [:mass_grid, :area_grid]`. If you want all possible fields then you can simply admit the outputs field altogether and all grid fields will be output. The `grid` field is the simulation grid, and then `dims` field specifies the dimensions of the grid you would like the output calculated on.
 
-Note that other than `outputs`, `∆tout`, `grid`, and `dims`, all values have default values so you only need to pass in arguments that you wish to change from these values. Furthermore, all inputs but ∆tout are keyword arguments, so you must use the keyword when passing in new values.
+Note that other than `∆tout`, `grid`, and `dims`, all values have default values so you only need to pass in arguments that you wish to change from these values. Furthermore, all inputs but ∆tout are keyword arguments, so you must use the keyword when passing in new values.
 
 #### OutputWriters
 Once you have created all of the types of output writers you need, you must combine them into one `OutputWriters` object that will be a simulation field. 
@@ -691,6 +693,20 @@ timestep_sim!(
 
 Note that we are working on a more elegant solution to coupling with Oceananigans and CliMA and this page will be updated once that is in place. 
 
+If you run your simulation in multiple parts and need to re-start your simulation from files, the `restart!` function will be a good place to start. However, note that it is quite simple and users may need to write their own restart function if they want any complex behavior. 
+
+The provided `restart!` function takes in the output file from both an `InitialStateOutputWriter` and a `CheckpointOutputWriter` to restart the simulation. In addition to providing these two files, the user must also provide the number of timesteps to run the next part of the simulation for (`new_nΔt`) and new output writers. The user also has an option to specify a non-zero starting timestep for the simulation using the keyword argument `start_tstep`.
+
+```
+restart!(
+ initial_state_fn,
+ checkpointer_fn,
+ new_nΔt,
+ new_output_writers;
+ start_tstep = 0,
+)
+```
+
 ### Plotting
 
 If your simulation has both a `FloeOutputWriter` and an `InitialStateOutputWriter`, you can use the built in plotting function to make an MP4 file with each frame as a timestep saved by the `FloeOutputWriter`. You do this as follows:

examples/restart_sim.jl

@@ -0,0 +1,113 @@
+using JLD2, Subzero, Random, Statistics
+
+const FT = Float64
+const Δt = 10
+const nΔt = 5000
+const n_part_sim = 3
+const nfloes = 20
+const L = 1e5
+const Δgrid = 1e4
+const hmean = 2
+const concentration = 0.7
+const uomax = 2
+
+dirs = [joinpath("output/restart_sim", "run_" * string(i)) for i in 1:n_part_sim]
+
+# Build initial model / simulation
+grid = RegRectilinearGrid(
+ (0.0, L),
+ (0.0, L),
+ Δgrid,
+ Δgrid,
+)
+ngrid = Int(L/Δgrid) + 1
+ygrid = range(0,L,ngrid)
+
+uoprofile = @. uomax * (1 - abs(1 - 2 * ygrid/L))
+uvels_ocean = repeat(
+ uoprofile,
+ outer = (1, ngrid),
+)
+ocean = Ocean(
+ uvels_ocean',
+ zeros(grid.Nx + 1, grid.Ny + 1),
+ zeros(grid.Nx + 1, grid.Ny + 1),
+)
+
+atmos = Atmos(FT, grid, 0.0, 0.0, 0.0)
+
+nboundary = PeriodicBoundary(North, grid)
+sboundary = PeriodicBoundary(South, grid)
+eboundary = PeriodicBoundary(East, grid)
+wboundary = PeriodicBoundary(West, grid)
+domain = Domain(nboundary, sboundary, eboundary, wboundary)
+
+floe_settings = FloeSettings(subfloe_point_generator = SubGridPointsGenerator(grid, 2))
+floe_arr = initialize_floe_field(
+ FT,
+ nfloes,
+ [concentration],
+ domain,
+ hmean,
+ 0;
+ rng = Xoshiro(1),
+ floe_settings = floe_settings
+)
+
+model = Model(grid, ocean, atmos, domain, floe_arr)
+
+modulus = 1.5e3*(mean(sqrt.(floe_arr.area)) + minimum(sqrt.(floe_arr.area)))
+consts = Constants(E = modulus, f = 0, turnθ = 0)
+
+initwriter = InitialStateOutputWriter(dir = dirs[1], overwrite = true)
+checkpointer = CheckpointOutputWriter(
+ 250,
+ dir = dirs[1],
+ filename = "checkpoint.jld2",
+ overwrite = true,
+ jld2_kw = Dict{Symbol, Any}(),
+)
+floewriter = FloeOutputWriter(50, dir = dirs[1], overwrite = true)
+writers = OutputWriters(initwriter, floewriter, checkpointer)
+
+simulation = Simulation(
+ model = model,
+ consts = consts,
+ Δt = Δt,
+ nΔt = nΔt,
+ verbose = true,
+ writers = writers,
+ rng = Xoshiro(1),
+)
+
+# Run the first part of the simulation
+run!(simulation)
+
+# Run remaining parts of the simulation
+for i in 2:n_part_sim
+ # Build new output writers
+ global initwriter = InitialStateOutputWriter(initwriter; dir = dirs[i])
+ global floewriter = FloeOutputWriter(floewriter; dir = dirs[i])
+ writers = if i < n_part_sim
+ global checkpointer = CheckpointOutputWriter(checkpointer; dir = dirs[i])
+ OutputWriters(initwriter, floewriter, checkpointer)
+ else
+ OutputWriters(initwriter, floewriter)
+ end
+ # Run next part of the simulation
+ Subzero.restart!(
+ dirs[i-1] * "/initial_state.jld2",
+ dirs[i-1] * "/checkpoint.jld2",
+ nΔt, writers; start_tstep = nΔt * (i - 1))
+end
+
+# Plot all simulation parts
+for i in 1:n_part_sim
+ plot_sim(
+ joinpath(dirs[i], "floes.jld2"),
+ joinpath(dirs[i], "initial_state.jld2"),
+ Δt,
+ joinpath(dirs[i], "shear_flow.mp4"),
+ )
+ println("Simulation part $i plotted.")
+end