update box model

Updated PZ_box.jl to work with new box model interface

Author: johnryantaylor

Manifest.toml

@@ -1,9 +1,9 @@
 # This script runs the phytoplankton-zooplankton (PZ) model as a 0-D box model using OceanBioME
 
 # Import some required modules
+using Oceananigans
 using OceanBioME
 using OceanBioME: BoxModel
-import OceanBioME.BoxModels: update_boxmodel_state!
 using JLD2
 using Plots
 
@@ -18,46 +18,49 @@ biogeochemistry = PhytoplanktonZooplankton()
 #biogeochemistry = PhytoplanktonZooplankton(phytoplankton_growth_rate = 0.5)
 
 model = BoxModel(; biogeochemistry)
-model.Δt = 0.05
-model.stop_time = 50
 
 # Set the initial conditions
 set!(model, P = 0.1, Z = 0.1)
 
+simulation = Simulation(model; Δt = 0.05, stop_time = 50)
+
+simulation.output_writers[:fields] = JLD2Writer(model, model.fields; filename = "box_pz.jld2", schedule = IterationInterval(1), overwrite_existing = true)
+
+prog(sim) = @info "$(prettytime(time(sim))) in $(prettytime(simulation.run_wall_time))"
+
+simulation.callbacks[:progress] = Callback(prog, IterationInterval(1000))
+
 # ## Run the model and save the output every save_interval timesteps (should only take a few seconds)
-@info "Running the model..."
-run!(model, save_interval = 1, feedback_interval = Inf, save = SaveBoxModel("box_pz.jld2"))
+@info "Running the simulation..."
+run!(simulation)
 
 # Now, read the saved output and plot the results
 
-vars = (:P, :Z)
-file = jldopen("box_pz.jld2")
-times = parse.(Float64, keys(file["values"]))
-
-timeseries = NamedTuple{vars}(ntuple(t -> zeros(length(times)), length(vars)))
+# Now, read the data and plot the results. This is saved as a Field with metadata
+P_field = FieldTimeSeries("box_pz.jld2", "P")
+Z_field = FieldTimeSeries("box_pz.jld2", "Z")
 
-for (idx, time) in enumerate(times)
-    values = file["values/$time"]
-    for tracer in vars
-        getproperty(timeseries, tracer)[idx] = values[tracer]
-    end
-end
+# Extract the data at the single grid point corresponding to the box model
+P=P_field.data[1,1,1,:]
+Z=Z_field.data[1,1,1,:]
 
-close(file)
+# Extract the time from one of the variables
+t=P_field.times
 
-anim = @animate for i=1:5:length(times)
-  plt1 = plot(times[1:i], timeseries.P[1:i], linewidth = 2, xlabel = "time", ylabel = "P, Z", legend = :outertopleft, label = "P", linecolor = :blue);
-  scatter!(plt1, [times[i]], [timeseries.P[i]], markercolor = :blue, label = :none)
-  plot!(plt1, times[1:i], timeseries.Z[1:i], linewidth = 2, label="Z", linecolor = :red, xlimits = (0, maximum(times)), ylimits = (0, max(maximum(timeseries.P),maximum(timeseries.Z))));
-  scatter!(plt1, [times[i]], [timeseries.Z[i]], markercolor = :red, label = :none)
-  plt2 = plot(timeseries.P[1:i], timeseries.Z[1:i], linewidth = 2, xlabel="P", ylabel="Z", xlimits = (0, maximum(timeseries.P)), ylimits = (0,maximum(timeseries.Z)), linecolor = :black, legend = :none);
-  scatter!(plt2, [timeseries.P[i]], [timeseries.Z[i]], markercolor = :black, label = :none)
+# Make an animated plot of the evolution of the PZ model
+anim = @animate for i=1:5:length(P_field.times)
+  plt1 = plot(t[1:i], P[1:i], linewidth = 2, xlabel = "time", ylabel = "P, Z", legend = :outertopleft, label = "P", linecolor = :blue);
+  scatter!(plt1, [t[i]], [P[i]], markercolor = :blue, label = :none)
+  plot!(plt1, t[1:i], Z[1:i], linewidth = 2, label="Z", linecolor = :red, xlimits = (0, maximum(P_field.times)), ylimits = (0, max(maximum(P),maximum(Z))));
+  scatter!(plt1, [t[i]], [Z[i]], markercolor = :red, label = :none)
+  plt2 = plot(P[1:i], Z[1:i], linewidth = 2, xlabel="P", ylabel="Z", xlimits = (0, maximum(P)), ylimits = (0,maximum(Z)), linecolor = :black, legend = :none);
+  scatter!(plt2, [P[i]], [Z[i]], markercolor = :black, label = :none)
   plot(plt1, plt2, layout = (1,2))
 end
 
 mp4(anim, "PZ_box.mp4", fps = 20) # hide
 
-plt1 = plot(times, timeseries.P, linewidth = 2, xlabel = "time", ylabel = "P, Z", legend = :outertopleft, label = "P");
-plot!(plt1, times,timeseries.Z, linewidth = 2, label="Z");
-plt2 = plot(timeseries.P, timeseries.Z, linewidth = 2, xlabel="P", ylabel="Z", linecolor = :black, legend = :none);
+plt1 = plot(t, P, linewidth = 2, xlabel = "time", ylabel = "P, Z", legend = :outertopleft, label = "P");
+plot!(plt1, t, Z, linewidth = 2, label="Z");
+plt2 = plot(P, Z, linewidth = 2, xlabel="P", ylabel="Z", linecolor = :black, legend = :none);
 plot(plt1, plt2, layout = (1,2))