output to a netcdf file

[miniufo]

Hi, I am using this great package to do some idealized simulation. Just wondering if it can output the states into a netcdf file so that I can do some offline Lagrangian tracking. I am not too farmiliar with Julia, so I would like to do this using Python with a netcdf file. The default output into a jld2 file is quite new to me.

[glwagner]

I believe you can read JLD2 with h5py, but not sure that works for you or not.

[miniufo]

Thanks @glwagner, I tried using Python's xarray to load jld2 (with engine='h5netcdf') file but no data available. So I just used NCDataset.jl to do the output.

[glwagner]

A contribution that adds this as a feature that can be used by others is welcome! The benefit is the usual open source concept that sharing code leads to higher quality more featureful code.

[miniufo]

Here is my code for output a netcdf file:

using NCDatasets

filepath = "."
filename = joinpath(filepath, "/home/qianyk/Output2D_8192.nc")

dsOut = NCDatasets.Dataset(filename, "c")

defDim(dsOut, "x", n)
defDim(dsOut, "y", n)
defDim(dsOut, "time", Inf)  # unlimited time

defVar(dsOut, "x", grid.x, ("x",))     # x
defVar(dsOut, "y", grid.y, ("y",))     # y
defVar(dsOut, "time", prec, ("time",)) # time

defVar(dsOut, "u"  , prec, ("x", "y", "time"), chunksizes=(n, n, 1), deflatelevel=5)
defVar(dsOut, "v"  , prec, ("x", "y", "time"), chunksizes=(n, n, 1), deflatelevel=5)
defVar(dsOut, "vor", prec, ("x", "y", "time"), chunksizes=(n, n, 1), deflatelevel=5)

TwoDNavierStokes.updatevars!(prob)
# put the initial state
dsOut["time"][1] = 0
dsOut["u"  ][:,:,1] = Array(prob.vars.u)
dsOut["v"  ][:,:,1] = Array(prob.vars.v)
dsOut["vor"][:,:,1] = Array(prob.vars.ζ)
NCDatasets.sync(dsOut)

startwalltime = time()

l = 2
for t in 1:nsteps
    stepforward!(prob, diags, nsubs)
    
    if t%nsubs == 0
        # update vars only when outputing???
        TwoDNavierStokes.updatevars!(prob)
        
        uu = Array(prob.vars.u)
        vv = Array(prob.vars.v)
        zz = Array(prob.vars.ζ)
        
        dsOut["time"][l] = t * dt
        dsOut["u"   ][:,:,l] = uu
        dsOut["v"   ][:,:,l] = vv
        dsOut["vor" ][:,:,l] = zz
        
        NCDatasets.sync(dsOut)
        l = l + 1

        if t%(2*nsubs) == 0
            # output some diagnostics to make sure overflow is not happening
            maxu = maximum(uu)
            cfl = prob.clock.dt * maximum([maxu / grid.dx, maximum(vv) / grid.dy])
            log = @sprintf("t: %05d, (%.3f %%), cfl: %.3f, Umax: %.3f, E: %.4f, Z: %.4f, walltime: %.2f min",
              t, (t/nsteps*100), cfl, maxu, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60)
            println(log)
        end
    end
end

close(dsOut)

I am not quite farmiliar to julia, so I don't think I can integrate this part to the package and make it re-usable for others. Just sharing this here if others find it useful.

[glwagner]

Perhaps you are not familiar with julia, but you're off to a great start with julia it looks like!

I think what I was suggesting is more or less the process of putting this code into the Output struct, but I'm not sure. Maybe @navidcy can provide some guidance for others who want to use netcdf.

[glwagner]

@miniufo do you mind if we keep this issue open until this feature is supported?

[miniufo]

No problem. I thought people here prefer the jld2 file. Actually, I've no idea which format is more efficient. With netcdf, I can compress the data, which is useful to my case and it is transparent to users. It would be great if this package can output a netcdf, providing only a filename, a variable list, and the frequency of output.

[glwagner]

both are useful. JLD2 is faster, lighter weight and works on more systems (does not depend on a C library which cannot always be installed eg on powerpc). However, NetCDF has advantages, the key one being compatibility with xarray. Our goal is always to support as much as possible, we do not limit the scope

[glwagner]

JLD2 supports compression by default: https://juliaio.github.io/JLD2.jl/dev/compression/

but I dont' know how this compares to NetCDF compression.

There's some discussion about compression on the (now old) Oceananigans PR

CliMA/Oceananigans.jl#3599

[miniufo]

Thanks @glwagner. That's quite useful. I've played with SpeedyWeather.jl and its output of a netcdf file is really small (surprising to me). I've no idea if @milankl has already used BitInformation.jl to make SpeedyWeather's output smaller or just relies on NCDataset.jl's compressibility (specifying deflatelevel). I tried with netcdf output here (using deflatelevel as shown above) but the compression ratio is not as high as that in SpeedyWeather, which is a bit strange to me.

[milankl]

Yes SpeedyWeather uses bitrounding and netcdfs default zlib compression (deflate). It doesn't actually analyse the bitinformation it uses some default keepbits per variable.

[glwagner]

So is the bitrounding the extra source of magic here?

[milankl]

Yes. You have to remove some insignificant mantissa bits to make the data compressible. Otherwise the lossless compressor has a very hard time to pack the file any smaller as you're hitting Shannon's source coding theorem.