-
Notifications
You must be signed in to change notification settings - Fork 9
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add support for vertical integrals #748
Labels
enhancement
New feature or request
Comments
As discussed on slack, I think definite integrals are well-defined: you sum multiplying by the dxdxi[3,3] term, and return a field on the horizontal space. For indefinite integrals, I think we should return a face-aligned field? I think we also want an option for upward or downward. |
Duplicate of #693. |
bors bot
added a commit
that referenced
this issue
Sep 28, 2022
962: Add support for definite column integrals r=simonbyrne 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. Co-authored-by: Charles Kawczynski <[email protected]> Co-authored-by: Simon Byrne <[email protected]>
bors bot
added a commit
that referenced
this issue
Oct 1, 2022
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]>
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Is your feature request related to a problem? Please describe.
We need vertical integrals in various places, including TurbulenceConvection and ClimaAtmos. The current
sum
is just a volume integral which is equivalent for a single column of area 1, but definitely not what we want to use for things like liquid water path.Please link to any related issues/PRs.
Describe the solution you'd like
Define both
vertical∫
that computes a vertical definite integralvertical_indefinite∫
that computes a vertical indefinite integralDescribe alternatives you've considered
A clear and concise description of any alternative solutions or features you've considered.
Additional context
Add any other context or screenshots about the feature request here.
The text was updated successfully, but these errors were encountered: