Call compute! in constructor for ComputedField (#4531)

* Update computed_field.jl

* Update oceananigans.bib

* remove obsolete compute statement

* add `compute` kw and retain original behavior in output_construction

* fix test plus docstrings

* Update src/OutputWriters/output_construction.jl

* fix doctests in operations.md

* fix doctest in seawater_density

* no need to compute

* no need to compute!

* fix docstrings

---------

Co-authored-by: Navid C. Constantinou <[email protected]>
Co-authored-by: tomaschor <[email protected]>

Author: 3 people

docs/oceananigans.bib

@@ -983,18 +983,18 @@ @article{BouZeid05
 }
 
 @article{Lan2022,
-    author = {Lan, Rihui and Ju, Lili and Wanh, Zhu and Gunzburger, Max and Jones, Philip},
-    title = {High-order multirate explicit time-stepping schemes for the baroclinic-barotropic split dynamics in primitive equations},
-    journal = {Journal of Computational Physics},
-    volume = {457},
-    pages = {111050},
-    year = {2022},
-    doi = {https://doi.org/10.1016/j.jcp.2022.111050},
+	author = {Lan, Rihui and Ju, Lili and Wanh, Zhu and Gunzburger, Max and Jones, Philip},
+	title = {High-order multirate explicit time-stepping schemes for the baroclinic-barotropic split dynamics in primitive equations},
+	journal = {Journal of Computational Physics},
+	volume = {457},
+	pages = {111050},
+	year = {2022},
+  doi = {10.1016/j.jcp.2022.111050},
 }
 
 @inproceedings{Verstappen14,
   title={Numerical scale separation in large-eddy simulation},
-  author={Verstappen, RWCP and Rozema, W and Bae, HJ},
+  author={Verstappen, Roel and Rozema, Wybe and Bae, H. Jane},
   booktitle={Proceedings of the Summer Program},
   pages={417--426},
   year={2014}

docs/src/operations.md

@@ -245,26 +245,10 @@ set!(c, 1)
 ├── operand: Integral of BinaryOperation at (Center, Center, Center) over dims (1, 2)
 ├── status: time=0.0
 └── data: 1×1×11 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:8) with eltype Float64 with indices 1:1×1:1×-2:8
-    └── max=0.0, min=0.0, mean=0.0
+    └── max=2.55032e14, min=2.55032e14, mean=2.55032e14
 ```
 
 A few remarks: note that the `∫c` has locations `Nothing, Nothing, Center`; this is because we have integrated in the first two dimensions and thus it's `reduced over dims = (1, 2)`.
-Further note that `∫c` is full of zeros; its max, min, and mean values are all 0.
-No computation has been done yet.
-To compute `∫c`, we call [`compute!`](@ref),
-
-```jldoctest operations_avg_int
-compute!(∫c)
-
-# output
-1×1×5 Field{Nothing, Nothing, Center} reduced over dims = (1, 2) on LatitudeLongitudeGrid on CPU
-├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (1, 1, 5)
-├── grid: 60×10×5 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
-├── operand: Integral of BinaryOperation at (Center, Center, Center) over dims (1, 2)
-├── status: time=0.0
-└── data: 1×1×11 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:8) with eltype Float64 with indices 1:1×1:1×-2:8
-    └── max=2.55032e14, min=2.55032e14, mean=2.55032e14
-```
 
 Above we see that the max, min and mean of the field are all the same.
 Let's check that these values are what we expect:
@@ -296,15 +280,14 @@ conditional_∫c = Field(Integral(c, dims=(1, 2), condition=cond)) # integrate o
 ├── operand: Integral of ConditionalOperation of BinaryOperation at (Center, Center, Center) with condition cond (generic function with 1 method) over dims (1, 2)
 ├── status: time=0.0
 └── data: 1×1×11 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:8) with eltype Float64 with indices 1:1×1:1×-2:8
-    └── max=0.0, min=0.0, mean=0.0
+    └── max=1.27516e14, min=1.27516e14, mean=1.27516e14
 ```
 
 Above we have attached a condition to the operand.
 Now the operand is applied only when the condition is true.
-Let's compute and see if we get 1/4 of the area of the sphere
+Let's see if that is 1/4 of the area of the sphere
 
 ```jldoctest operations_avg_int
-compute!(conditional_∫c)
 conditional_∫c[1, 1, 1] ≈ π * grid.radius^2 # area of spherical zone with 0ᵒ ≤ φ ≤ 30ᵒ
 
 # output
@@ -319,7 +302,6 @@ cond_array[:, 1:5, :] .= false # set the first half of the latitude range to fal
 
 conditional_∫c = Field(Integral(c, dims=(1, 2), condition=cond_array))
 
-compute!(conditional_∫c)
 conditional_∫c[1, 1, 1] ≈ π * grid.radius^2 # area of spherical zone with 0ᵒ ≤ φ ≤ 30ᵒ
 
 # output

docs/src/quick_start.md

@@ -37,7 +37,6 @@ using CairoMakie
 
 u, v, w = model.velocities
 ζ = Field(∂x(v) - ∂y(u))
-compute!(ζ)
 
 heatmap(ζ, axis=(; aspect=1))
 ```
@@ -48,7 +47,6 @@ A few more time-steps, and it's starting to get a little diffuse!
 simulation.stop_iteration += 400
 run!(simulation)
 
-compute!(ζ)
 heatmap(ζ, axis=(; aspect=1))
 ```
 
@@ -87,7 +85,6 @@ run!(simulation)
 
 u, v, w = model.velocities
 ζ = Field(∂x(v) - ∂y(u))
-compute!(ζ)
 
 fig = Figure(size=(1200, 600))
 axζ = Axis(fig[1, 1], aspect=1, title="vorticity")

examples/horizontal_convection.jl

@@ -306,12 +306,9 @@ nothing #hide
 
 for n = 1:length(t)
     ke = Field(Integral(1/2 * s_timeseries[n]^2 / (Lx * H)))
-    compute!(ke)
     kinetic_energy[n] = ke[1, 1, 1]
 
     χ = Field(Integral(χ_timeseries[n] / (Lx * H)))
-    compute!(χ)
-
     Nu[n] = χ[1, 1, 1] / χ_diff
 end
 

examples/shallow_water_Bickley_jet.jl

@@ -96,9 +96,8 @@ uh, vh, h = model.solution
 u = uh / h
 v = vh / h
 
-## Build and compute mean vorticity discretely
+## Build mean vorticity discretely
 ω = Field(∂x(v) - ∂y(u))
-compute!(ω)
 
 ## Copy mean vorticity to a new field
 ωⁱ = Field((Face, Face, Nothing), model.grid)

src/AbstractOperations/computed_field.jl

@@ -17,6 +17,7 @@ const ComputedField = Field{<:Any, <:Any, <:Any, <:OperationOrFunctionField}
           data = nothing,
           indices = indices(operand),
           boundary_conditions = FieldBoundaryConditions(operand.grid, location(operand)),
+          compute = true,
           recompute_safely = true)
 
 Return a field `f` where `f.data` is computed from `f.operand` by calling `compute!(f)`.
@@ -29,16 +30,20 @@ Keyword arguments
 
 `boundary_conditions` (`FieldBoundaryConditions`): Boundary conditions for `f`.
 
-`recompute_safely` (`Bool`): whether or not to _always_ "recompute" `f` if `f` is
+`recompute_safely` (`Bool`): Whether or not to _always_ "recompute" `f` if `f` is
                              nested within another computation via an `AbstractOperation` or `FunctionField`.
-                             If `data` is not provided then `recompute_safely=false` and
+                             If `data` is not provided then `recompute_safely = false` and
                              recomputation is _avoided_. If `data` is provided, then
                              `recompute_safely = true` by default.
+
+`compute`: If `true`, `compute!` the `Field` during construction, otherwise if `false`, initialize with zeros.
+           Default: `true`.
 """
 function Field(operand::OperationOrFunctionField;
                data = nothing,
                indices = indices(operand),
                boundary_conditions = FieldBoundaryConditions(operand.grid, location(operand)),
+               compute = true,
                recompute_safely = true)
 
     grid = operand.grid
@@ -54,13 +59,20 @@ function Field(operand::OperationOrFunctionField;
 
     status = recompute_safely ? nothing : FieldStatus()
 
-    return Field(loc, grid, data, boundary_conditions, indices, operand, status)
+    computed_field = Field(loc, grid, data, boundary_conditions, indices, operand, status)
+
+    if compute
+        compute!(computed_field)
+    end
+
+    return computed_field
 end
 
 """
-    compute!(comp::ComputedField)
+    compute!(comp::ComputedField, time=nothing)
 
 Compute `comp.operand` and store the result in `comp.data`.
+If `time` then computation happens if `time != field.status.time`.
 """
 function compute!(comp::ComputedField, time=nothing)
     # First compute `dependencies`:

src/Fields/scans.jl

@@ -41,6 +41,7 @@ Base.summary(s::Scan) = string(summary(s.type), " ",
 function Field(scan::Scan;
                data = nothing,
                indices = indices(scan.operand),
+               compute = true,
                recompute_safely = true)
 
     operand = scan.operand
@@ -56,7 +57,13 @@ function Field(scan::Scan;
     boundary_conditions = FieldBoundaryConditions(grid, loc, indices)
     status = recompute_safely ? nothing : FieldStatus()
 
-    return Field(loc, grid, data, boundary_conditions, indices, scan, status)
+    scan_field = Field(loc, grid, data, boundary_conditions, indices, scan, status)
+
+    if compute
+         compute!(scan_field)
+    end
+
+    return scan_field
 end
 
 const ScannedComputedField = Field{<:Any, <:Any, <:Any, <:Scan}
@@ -105,7 +112,7 @@ Base.show(io::IO, s::Scan) =
 
 Return a `Reduction` of `operand` with `reduce!`, where `reduce!` can be called with
 
-```
+```julia
 reduce!(field, operand)
 ```
 
@@ -118,18 +125,12 @@ Example
 using Oceananigans
 
 Nx, Ny, Nz = 3, 3, 3
-
 grid = RectilinearGrid(size=(Nx, Ny, Nz), x=(0, 1), y=(0, 1), z=(0, 1),
                        topology=(Periodic, Periodic, Periodic))
 
 c = CenterField(grid)
-
 set!(c, (x, y, z) -> x + y + z)
-
 max_c² = Field(Reduction(maximum!, c^2, dims=3))
-
-compute!(max_c²)
-
 max_c²[1:Nx, 1:Ny]
 
 # output
@@ -151,7 +152,7 @@ location(r::Reduction) = reduced_location(location(r.operand); dims=r.dims)
 
 Return a `Accumulation` of `operand` with `accumulate!`, where `accumulate!` can be called with
 
-```
+```julia
 accumulate!(field, operand; dims)
 ```
 
@@ -164,18 +165,12 @@ Example
 using Oceananigans
 
 Nx, Ny, Nz = 3, 3, 3
-
 grid = RectilinearGrid(size=(Nx, Ny, Nz), x=(0, 1), y=(0, 1), z=(0, 1),
                        topology=(Periodic, Periodic, Periodic))
 
 c = CenterField(grid)
-
 set!(c, (x, y, z) -> x + y + z)
-
 cumsum_c² = Field(Accumulation(cumsum!, c^2, dims=3))
-
-compute!(cumsum_c²)
-
 cumsum_c²[1:Nx, 1:Ny, 1:Nz]
 
 # output

src/Models/seawater_density.jl

@@ -106,16 +106,6 @@ julia> density_field = Field(density_operation)
 │   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
 ├── operand: KernelFunctionOperation at (Center, Center, Center)
 ├── status: time=0.0
-└── data: 1×1×106 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:103) with eltype Float64 with indices 1:1×1:1×-2:103
-    └── max=0.0, min=0.0, mean=0.0
-
-julia> compute!(density_field)
-1×1×100 Field{Center, Center, Center} on RectilinearGrid on CPU
-├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
-├── boundary conditions: FieldBoundaryConditions
-│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
-├── operand: KernelFunctionOperation at (Center, Center, Center)
-├── status: time=0.0
 └── data: 1×1×106 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:103) with eltype Float64 with indices 1:1×1:1×-2:103
     └── max=1032.38, min=1027.73, mean=1030.06
 ```

src/OutputWriters/output_construction.jl

@@ -54,7 +54,11 @@ end
 
 function construct_output(user_output::Union{AbstractField, Reduction}, grid, user_indices, with_halos)
     indices = output_indices(user_output, grid, user_indices, with_halos)
-    return Field(user_output; indices)
+
+    # Don't compute AbstractOperations or Reductions
+    additional_kw = user_output isa Field ? NamedTuple() : (; compute=false)
+
+    return Field(user_output; indices, additional_kw...)
 end
 
 #####

test/test_field_scans.jl

@@ -376,10 +376,12 @@ interior_array(a, i, j, k) = Array(interior(a, i, j, k))
             c = CenterField(grid)
             set!(c, (x, y, z) -> x + y + z)
 
-            max_c² = Field(Reduction(maximum, c^2, dims=3))
+            max_c² = Field(Reduction(maximum!, c^2, dims=3))
             ∫max_c² = Integral(max_c², dims=(1, 2))
-            compute!(∫max_c²)
             @test ∫max_c² isa Reduction
+            ∫max_c²_field = Field(∫max_c²)
+            @test ∫max_c²_field isa Field
+            @test ∫max_c²_field.operand === ∫max_c²
 
             @info "  Testing conditional reductions of immersed Fields [$(typeof(arch))]"