Skip to content
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 definite column integrals #962

Merged
merged 3 commits into from
Oct 1, 2022
Merged

Conversation

charleskawczynski
Copy link
Member

@charleskawczynski charleskawczynski commented Sep 26, 2022

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. 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:
Screen Shot 2022-09-26 at 7 57 05 AM

The interface for the added function captures two cases:

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

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.

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:

@. 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,

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

src/Fields/integrals.jl Outdated Show resolved Hide resolved
src/Fields/integrals.jl Outdated Show resolved Hide resolved
src/Fields/integrals.jl Outdated Show resolved Hide resolved
src/Spaces/finitedifference.jl Outdated Show resolved Hide resolved
src/Fields/Fields.jl Outdated Show resolved Hide resolved
src/Fields/integrals.jl Outdated Show resolved Hide resolved
src/Fields/integrals.jl Outdated Show resolved Hide resolved
src/Fields/integrals.jl Outdated Show resolved Hide resolved
src/Fields/integrals.jl Outdated Show resolved Hide resolved
@simonbyrne
Copy link
Member

Also cross referencing #693

@charleskawczynski charleskawczynski force-pushed the ck/vert_col_integrals branch 2 times, most recently from bddec0e to 4ba8810 Compare September 27, 2022 21:50
@charleskawczynski charleskawczynski changed the title Add support for allocation-free definite column integrals Add support for definite column integrals Sep 27, 2022
@simonbyrne
Copy link
Member

bors r+

bors bot added a commit that referenced this pull request 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
Copy link
Contributor

bors bot commented Sep 28, 2022

Build failed:

Apply some suggestions

More suggested changes

Avoid parent in loop

Update src/Operators/integrals.jl

Use nlevels function
@charleskawczynski
Copy link
Member Author

bors r+

@bors bors bot merged commit 175ec42 into main Oct 1, 2022
@bors bors bot deleted the ck/vert_col_integrals branch October 1, 2022 04:03
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants