Abstract Operations on Discrete Data

[s-aucoin]

Hi,

I have a question about using abstract operations on discrete data, specifically FieldTimeSeries objects. I understand that you can use these operations (like the spatial derivatives) before running a simulation to calculate and store additional variables like the following:

u, v, w = model.velocities
         ξ = (∂y(w) - ∂z(v))                                  # (1/s) x-z vorticity
         η = (∂z(u) - ∂x(w))                                  # (1/s) y-z vorticity
         ζ = (∂x(v) - ∂y(u))                                  # (1/s) x-y vorticity
      vort = (; ξ, η, ζ)

simulation.output_writers[:slices] =
    JLD2OutputWriter(model, merge(model.velocities, model.tracers, vort),
                     filename = filepath,
                     schedule = TimeInterval(output_interval),
                     overwrite_existing = true)

Due mainly to storage concerns, I only save the basic variables in my output file, so my question is can the same thing (or an analogous operation) be done after the simulation with the FieldTimeSeries data? Something like:

time_series = (u = FieldTimeSeries(load_filepath, "u"),            # (m s⁻¹) x component of velocity
               v = FieldTimeSeries(load_filepath, "v"),            # (m s⁻¹) y component of velocity
               w = FieldTimeSeries(load_filepath, "w"))            # (m s⁻¹) z component of velocity

extra_time_series = (ξ = (∂y(time_series.w) - ∂z(time_series.v)),                        # (s⁻¹) x-z vorticity
                     η = (∂z(time_series.u) - ∂x(time_series.w)),                        # (s⁻¹) y-z vorticity
                     ζ = (∂x(time_series.v) - ∂y(time_series.u)))                        # (s⁻¹) x-y vorticity

new_time_series = merge(time_series, extra_time_series)

I know I could define a discrete derivative myself separately, but I would like to ensure consistency in the discretization, and ideally would like to include the whole domain (using the boundary conditions). Is this possible?

[glwagner]

Yes, except that you cannot perform AbstractOperations directly on FieldTimeSeries. Instead, you have to extract a Field at a particular time from the FieldTimeSeries and perform operations on that:

ut = FieldTimeSeries(filename, "u")
vt = FieldTimeSeries(filename, "v")

for n = 1:length(ut.times)
    un = ut[n]
    vn = vt[n]
    ζn = ∂x(vn) - ∂y(un)
end

Also note that you probably want to save your data using with_halo = true so that you have halo points (which will allow you to compute derivatives on the boundaries).

[glwagner]

A quick search shows that jupyter kernels can quit when they run out of memory, so that could make sense. Oddly, that must have to do with the memory used during compilation (not the memory consumed by the arrays), since otherwise it would not quit during compute!. I wonder if it works at the REPL.

[s-aucoin]

I tried running the script as a .jl file using the REPL, but it also crashed. It's strange because it doesn't seem like a very memory heavy task, and my computer has a decent amount (16 GB).

[glwagner]

How does the REPL crash? Do you mean that it throws an error?

[s-aucoin]

It does throw an error, but then immediately closes the REPL. Running it in VSCode does the same thing but also gives the message that it exits with exit code 1.

[glwagner]

I've never seen that before! Crazy!

I have seen Out of memory errors. However, another thing that can happen is that your CPU starts using "swap memory" (which is on disk and extremely slow). Then things will continue to run but probably would never finish because its so slow.