Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
962: Add support for definite column integrals r=charleskawczynski a=charleskawczynski This PR adds `column_integral_definite!`, a non-allocating definite integral for columns. The idea is that it will replace the existing (allocating) definite integral [here](https://github.com/CliMA/ClimaAtmos.jl/blob/44cee55def2a51433c5dcdcdae010b7777cd741b/examples/hybrid/sphere/baroclinic_wave_utilities.jl#L167-L186). However, this implementation has several advantages: - It's compatible with `bycolumn`, so the computation can occur in parallel - It's allocation-free - It's _tested_ on single column and sphere configurations Looking at the flame graph for allocations (using the Julia 1.8 Allocs module with `sample_rate=1`), this is responsible for most of the remaining allocations in ClimaAtmos: <img width="1993" alt="Screen Shot 2022-09-26 at 7 57 05 AM" src="https://user-images.githubusercontent.com/1880641/192353757-03101c41-2c0b-4ccb-a8b9-43faa78680f8.png"> The interface for the added function captures two cases: ```julia function column_integral_definite!(col∫field::Field, field::Field) bycolumn(axes(field)) do colidx column_integral_definite!(col∫field[colidx], field[colidx]) nothing end return nothing end ``` and ```julia function column_integral_definite!( col∫field::PointField, field::ColumnField, ) `@inbounds` col∫field[] = column_integral_definite(field) return nothing end ``` A step towards closing #943, #748, [CA#686](CliMA/ClimaAtmos.jl#686). ## A note on an alternative approach Ideally, we would write this function as `column_integral_definite!(fn, args...)` where we might be able to write a broadcast statement like: ```julia `@.` f2D = column_integral_definite(f3D) do z f3D.x * z^2 end ``` however, this would require us to define broadcasting between planar and 3D domains, which is not trivial (or maybe not possible) because of ambiguities. The ambiguities arise because (2D, 3D) broadcasting may want different things in different cases, for example: - (f2D, f3D) -> f3D: mul full 3D field by planar surface value - (f2D, f3D) -> f2D: perform reduction over z-coordinate to yield 2D field The situation is similar when thinking about what happens when we make views. For example, ```julia Fields.bycolumn(axes(f3D)) do colidx `@.` f2D[colidx] = column_integral_definite(f3D[colidx]) do z f3D[colidx].x * z^2 end end ``` Now, we have to define how `DataF` data layouts broadcast with `VF`. Again, we have two cases: - (f0D, f1D) -> f1D: mul full 1D field by 0D field - (f0D, f1D) -> f0D: perform reduction over z-coordinate to yield 0D field My vote/preference is to support the first cases (which is partially supported already) and write custom functions (e.g., reductions) that operate on single fields for the second case. A step towards closing [CA#686](CliMA/ClimaAtmos.jl#686) Co-authored-by: Charles Kawczynski <[email protected]>
- Loading branch information
175ec42There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@JuliaRegistrator register
175ec42There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Registration pull request created: JuliaRegistries/General/69332
After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.
This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via: