search

CliMA / Oceananigans.jl

🌊 Julia software for fast, friendly, flexible, ocean-flavored fluid dynamics on CPUs and GPUs
https://clima.github.io/OceananigansDocumentation/stable
MIT License
985 stars 193 forks source link

Cannot output Windowed fields with together with standard fields #3835

Open simone-silvestri opened 2 weeks ago

simone-silvestri commented 2 weeks ago

MWE 

grid = RectilinearGrid(size = (5, 5, 5), extent = (1, 1, 1))
model = HydrostaticFreeSurfaceModel(; grid)
c = CenterField(grid; indices = (:, :, grid.Nz))
d = CenterField(grid)
JLD2OutputWriter(model, (; c, d), filename = "test", schedule = IterationInterval(1))

fails with 

julia> JLD2OutputWriter(model, (; c, d), filename = "test", schedule = IterationInterval(1))
ERROR: ArgumentError: view indices (1:5, 1:5, 1:5) do not intersect field indices (Colon(), Colon(), 5:5)
Stacktrace:
 [1] view(f::Field{…}, i::UnitRange{…}, j::UnitRange{…}, k::UnitRange{…})
   @ Oceananigans.Fields ~/development/Oceananigans.jl/src/Fields/field.jl:319
 [2] Field
   @ ~/development/Oceananigans.jl/src/Fields/field.jl:184 [inlined]
 [3] construct_output(user_output::Field{…}, grid::RectilinearGrid{…}, user_indices::Tuple{…}, with_halos::Bool)
   @ Oceananigans.OutputWriters ~/development/Oceananigans.jl/src/OutputWriters/output_construction.jl:46
 [4] (::Oceananigans.OutputWriters.var"#28#29"{Tuple{…}, Bool, HydrostaticFreeSurfaceModel{…}})(name::Symbol)
   @ Oceananigans.OutputWriters ./none:0
 [5] iterate
   @ ./generator.jl:47 [inlined]
 [6] merge(a::@NamedTuple{}, itr::Base.Generator{Tuple{…}, Oceananigans.OutputWriters.var"#28#29"{…}})
   @ Base ./namedtuple.jl:361
 [7] NamedTuple(itr::Base.Generator{Tuple{…}, Oceananigans.OutputWriters.var"#28#29"{…}})
   @ Base ./namedtuple.jl:151
 [8] JLD2OutputWriter(model::HydrostaticFreeSurfaceModel{…}, outputs::@NamedTuple{…}; filename::String, schedule::IterationInterval, dir::String, indices::Tuple{…}, with_halos::Bool, array_type::Type, file_splitting::Oceananigans.OutputWriters.NoFileSplitting, overwrite_existing::Bool, init::typeof(Oceananigans.OutputWriters.noinit), including::Vector{…}, verbose::Bool, part::Int64, jld2_kw::Dict{…})
   @ Oceananigans.OutputWriters ~/development/Oceananigans.jl/src/OutputWriters/jld2_output_writer.jl:185
 [9] top-level scope
   @ REPL[8]:1
Some type information was truncated. Use `show(err)` to see complete types.

It looks like there are problems also if trying to save the windowed field by itself

julia> JLD2OutputWriter(model, (; c), filename = "test1", schedule = IterationInterval(1))
ERROR: ArgumentError: view indices (1:5, 1:5, 1:5) do not intersect field indices (Colon(), Colon(), 5:5)
Stacktrace:
 [1] view(f::Field{…}, i::UnitRange{…}, j::UnitRange{…}, k::UnitRange{…})
   @ Oceananigans.Fields ~/development/Oceananigans.jl/src/Fields/field.jl:319
 [2] Field
   @ ~/development/Oceananigans.jl/src/Fields/field.jl:184 [inlined]
 [3] construct_output(user_output::Field{…}, grid::RectilinearGrid{…}, user_indices::Tuple{…}, with_halos::Bool)
   @ Oceananigans.OutputWriters ~/development/Oceananigans.jl/src/OutputWriters/output_construction.jl:46
 [4] (::Oceananigans.OutputWriters.var"#28#29"{Tuple{…}, Bool, HydrostaticFreeSurfaceModel{…}})(name::Symbol)
   @ Oceananigans.OutputWriters ./none:0
 [5] iterate
   @ ./generator.jl:47 [inlined]
 [6] merge(a::@NamedTuple{}, itr::Base.Generator{Tuple{…}, Oceananigans.OutputWriters.var"#28#29"{…}})
   @ Base ./namedtuple.jl:361
 [7] NamedTuple(itr::Base.Generator{Tuple{Symbol}, Oceananigans.OutputWriters.var"#28#29"{Tuple{…}, Bool, HydrostaticFreeSurfaceModel{…}}})
   @ Base ./namedtuple.jl:151
 [8] JLD2OutputWriter(model::HydrostaticFreeSurfaceModel{…}, outputs::@NamedTuple{…}; filename::String, schedule::IterationInterval, dir::String, indices::Tuple{…}, with_halos::Bool, array_type::Type, file_splitting::Oceananigans.OutputWriters.NoFileSplitting, overwrite_existing::Bool, init::typeof(Oceananigans.OutputWriters.noinit), including::Vector{…}, verbose::Bool, part::Int64, jld2_kw::Dict{…})
   @ Oceananigans.OutputWriters ~/development/Oceananigans.jl/src/OutputWriters/jld2_output_writer.jl:185
 [9] top-level scope
   @ REPL[11]:1
Some type information was truncated. Use `show(err)` to see complete types.

However, by specifying the indices it works

julia> JLD2OutputWriter(model, (; c), filename = "test1", schedule = IterationInterval(1), indices = (:, :, grid.Nz))
JLD2OutputWriter scheduled on IterationInterval(1):
├── filepath: ./test1.jld2
├── 1 outputs: c
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 20.0 KiB
glwagner commented 2 weeks ago

What happens if you use with_halos=true? This might be a problem of trying to restrict to the interior.

simone-silvestri commented 2 weeks ago

Indeed:

julia> JLD2OutputWriter(model, (; c), filename = "test3", schedule = IterationInterval(1), with_halos = true)
JLD2OutputWriter scheduled on IterationInterval(1):
├── filepath: ./test3.jld2
├── 1 outputs: c
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 20.3 KiB

julia> JLD2OutputWriter(model, (; c, d), filename = "test4", schedule = IterationInterval(1), with_halos = true)
JLD2OutputWriter scheduled on IterationInterval(1):
├── filepath: ./test4.jld2
├── 2 outputs: (c, d)
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 25.9 KiB

I ll take a look at the restrict_to_interior function

glwagner commented 2 weeks ago

It's trying to use the wrong indices:

ERROR: ArgumentError: view indices (1:5, 1:5, 1:5) do not intersect field indices (Colon(), Colon(), 5:5)

It should be trying to view with indices (1:5, 1:5, 5:5).