Add operator_matrix to MatrixFields, along with tests and docs #1399

dennisYatunin · 2023-07-28T01:03:07Z

Purpose

Second PR of #1230. Refactors what is currently in src/operators/operator2stencil.jl.

Content

Adds the function operator_matrix(op), which will replace Operator2Stencil(op). This function has a detailed docstring and throws descriptive errors for operators without well-defined operator matrices.
Improves the usability of operator matrices. With this new interface, users will no longer need to specify intermediate fields to compute operator matrices. For example, with our pre-existing code, the matrix of the interpolation operator (a bidiagonal matrix whose entries are -1/2 and +1/2) needs to be computed with @. matrix_field = interp_matrix(ones_field), where ones_field is a field filled with the number 1 that is used to infer the space and entry type of the matrix. Now, this matrix can just be computed with @. matrix_field = interp_matrix(). In order to make this work, interp_matrix is now defined as a "lazy operator". When a broadcast expression containing lazy operators is evaluated, each lazy operator is replaced with an actual operator, and it is given one or more fields as input arguments. In this case, interp_matrix is given the local geometry field as an input argument, and this field is used to infer the space and entry type of the operator matrix.
This usability improvement slightly changes the computation of derivative matrices. With our pre-existing code, @. op(func(field)) is equivalent to @. op_matrix(ones_field) ⋅ func(field), and the derivative of this expression with respect to field can be specified as @. op_matrix(func'(field)), where func' is the derivative of the point-wise function func. With this new interface, @. op(func(field)) is equivalent to @. op_matrix() ⋅ func(field), and the derivative can be specified as @. op_matrix() ⋅ DiagonalMatrixRow(func'(field)). Similarly, the derivative of @. op2(func2(op1(func1(field)))) with respect to field is op2_matrix(func2'(op1(func1(field)))) ⋅ op1_matrix(func1'(field)) with our pre-existing code and op2_matrix() ⋅ DiagonalMatrixRow(func2'(op1(func1(field)))) ⋅ op1_matrix() ⋅ DiagonalMatrixRow(func1'(field)) with the new interface. Although the new interface leads to longer derivative expressions, those expressions are more similar to how the chain rule is usually written out, and they can be debugged/analyzed more incrementally.
Adds support for computing operator matrices of multi-argument operators. For example, if op is the Upwind3rdOrderBiasedProductC2F operator, then @. op(velocity_field, tracer_field) is equivalent to @. op_matrix(velocity_field) ⋅ tracer_field. The implementation is similar to that of single-argument operators, except that it does not require the use of "lazy operators" (since there is already a field being passed to the operator matrix, the local geometry field can be obtained from that field during the evaluation of Base.Broadcast.broadcasted).
Adds support for computing operator matrices of operators with Extrapolate boundary conditions. These boundary conditions cause the matrices to have larger bandwidths than other boundary conditions.
Tests the operator_matrix function with every valid combination of finite difference operators and boundary conditions. The tests check for correctness, type stability, and lack of allocations. The tests are run on both CPUs and GPUs. In addition, the tests print out how the performance of @. op_matrix() ⋅ field compares to the performance of @. op(field); the two expressions are similarly fast on GPUs (between -70% and +40% relative change in speed), though the operator matrix expressions tend to be slower on CPUs.
Tests a few more complicated broadcast expressions involving products and linear combinations of operator matrices. These tests indicate that operator matrices are similarly performant to regular matrix fields.
Modifies test_field_broadcast so that it also tests whether get_result generates the same result as set_result!.
Modifies the * method for BandMatrixRow so that matrix fields can be scaled by vectors/covectors in addition to numbers. This simplifies a few of the complicated broadcast tests.
Fixes a typo (Geometery) in the method of stencil_left_boundary for GradientF2C.
Adds a missing method to rpromote_type that was preventing empty matrix rows from being constructed.
Modifies the Base.Broadcast.broadcasted method for FiniteDifferenceOperator and SpectralElementOperator so that lazy operators can work correctly (the original versions of these methods would always overwrite the LazyOperatorStyle with StencilStyle or SpectralStyle, respectively).
Unfortunately, this PR also adds 2 method invalidations. These are due to a new definition of broadcasted for lazy operators:
```
Base.Broadcast.broadcasted(::LazyOperatorStyle, f::F, args...) where {F} =
    LazyOperatorBroadcasted(f, args)
```
As explained in an accompanying code comment, removing this method requires modifying several other method definitions, and one of these modifications adds 11 invalidations. So, there doesn't seem to be a good way to avoid these 2 invalidations.

dennisYatunin · 2023-07-28T18:43:29Z

bors try

bors · 2023-07-28T19:09:41Z

try

Build failed:

buildkite/climacore-ci

dennisYatunin · 2023-07-31T17:24:48Z

bors try

bors · 2023-07-31T17:58:56Z

try

Build failed:

buildkite/climacore-ci

dennisYatunin · 2023-07-31T18:24:04Z

bors try

bors · 2023-07-31T19:08:30Z

try

Build failed:

test-os (ubuntu-20.04)

charleskawczynski · 2023-07-31T20:27:44Z

test/MatrixFields/matrix_field_broadcasting.jl

-    end
-end
+
+include("matrix_field_test_utils.jl")


dennisYatunin · 2023-07-31T21:26:51Z

bors try

bors · 2023-07-31T22:58:24Z

try

Build failed:

test-os (ubuntu-20.04)

src/MatrixFields/lazy_operators.jl

dennisYatunin · 2023-08-09T20:41:11Z

This PR has gotten too long, so I'm going to peel off a few smaller PRs before flagging this for a final review.

1423: Clean up MatrixFields tests r=dennisYatunin a=dennisYatunin Peeled off from #1399. Cleans up `test/MatrixFields` by creating `matrix_field_test_utils.jl` and moving all shared testing code there. Simplifies the matrix broadcast tests by removing a lot of variables that were being passed around as function arguments (e.g., replaces anonymous functions of the form `(args...) -> f(args...)` with `() -> f(args...)`, since both give the same benchmarking results). Sets the default values of `time_ratio_limit` and `max_eps_error_limit` to 10, which should help us avoid random CI failures. - [x] Code follows the [style guidelines](https://clima.github.io/ClimateMachine.jl/latest/DevDocs/CodeStyle/) OR N/A. - [x] Unit tests are included OR N/A. - [x] Code is exercised in an integration test OR N/A. - [x] Documentation has been added/updated OR N/A. Co-authored-by: Dennis Yatunin <dyatun@gmail.com>

1424: Fix matrix multiplication at boundaries r=dennisYatunin a=dennisYatunin Peeled off from #1399. Fixes a bug in matrix-matrix multiplication (i.e., the `MultiplyColumnwiseBandMatrixField()` operator with two matrix fields as inputs) that was causing entries from outside the second input matrix to be used in the computation of entries outside the product matrix. These entries are used for padding, allowing us to store the nonzero entries of band matrices in `Field`s with rectangular parent arrays, so using them during the computation of matrix-matrix products leads to an unnecessary performance loss. Moreover, it breaks the assumption that these outside entries have no effect on program behavior. In addition to fixing the bug, this PR - Updates the documentation of `MultiplyColumnwiseBandMatrixField` to reflect the change. - Adds a unit test which ensures that outside entries are never used during matrix-matrix multiplication. - Simplifies some of the code for `MultiplyColumnwiseBandMatrixField`, making sure that it reflects the updated documentation as closely as possible. - Introduces the `column_axes(matrix_field)` utility function, which can be used to obtain the space that corresponds to the columns of `matrix_field` (as opposed to `axes(matrix_field)`, which corresponds to the rows). This function simplifies the computation of interior/boundary indices during matrix multiplication. - Fixes a type instability in `field2arrays` and updates the corresponding unit test. This also decreases the memory allocated by the new unit test for outside entries. - Modifies the `show` method for matrix `Field`s so that it also displays the matrix shape. Co-authored-by: Dennis Yatunin <dyatun@gmail.com>

dennisYatunin · 2023-08-11T18:48:42Z

bors try

bors · 2023-08-11T20:54:12Z

try

Build failed:

buildkite/climacore-ci

dennisYatunin · 2023-08-11T21:28:50Z

bors try

bors · 2023-08-12T00:42:34Z

try

Build succeeded!

The publicly hosted instance of bors-ng is deprecated and will go away soon.

If you want to self-host your own instance, instructions are here.
For more help, visit the forum.

If you want to switch to GitHub's built-in merge queue, visit their help page.

dennisYatunin · 2023-08-14T16:17:39Z

There was one unit test whose reference function was taking too long to compile on Julia 1.9. I've disabled the reference function for now (so the test isn't going to check for correctness, but it will still check for allocations and type instabilities), though we should probably figure out how to re-enable it with reasonable compilation time in the future.

dennisYatunin · 2023-08-14T16:17:46Z

bors r+

bors · 2023-08-14T17:44:26Z

Build succeeded!

The publicly hosted instance of bors-ng is deprecated and will go away soon.

If you want to self-host your own instance, instructions are here.
For more help, visit the forum.

If you want to switch to GitHub's built-in merge queue, visit their help page.

dennisYatunin force-pushed the dy/operator_matrices branch 3 times, most recently from 96d568f to 5faf79d Compare July 28, 2023 18:37

dennisYatunin requested review from charleskawczynski and simonbyrne July 28, 2023 18:38

bors bot added a commit that referenced this pull request Jul 28, 2023

Try #1399:

b9352b1

dennisYatunin force-pushed the dy/operator_matrices branch from 5faf79d to 3ee210d Compare July 28, 2023 21:03

dennisYatunin mentioned this pull request Jul 28, 2023

Implicit Solver Interface Updates #1230

Closed

7 tasks

dennisYatunin force-pushed the dy/operator_matrices branch 5 times, most recently from a4371c6 to 9e8a6d3 Compare July 31, 2023 17:23

bors bot added a commit that referenced this pull request Jul 31, 2023

Try #1399:

ebf5ad6

dennisYatunin force-pushed the dy/operator_matrices branch from 9e8a6d3 to bde9f87 Compare July 31, 2023 18:23

bors bot added a commit that referenced this pull request Jul 31, 2023

Try #1399:

be33dc1

charleskawczynski reviewed Jul 31, 2023

View reviewed changes

test/MatrixFields/matrix_field_broadcasting.jl Outdated

end

end

include("matrix_field_test_utils.jl")

Copy link

Member

charleskawczynski Jul 31, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

❤️

dennisYatunin force-pushed the dy/operator_matrices branch from bde9f87 to 403cb67 Compare July 31, 2023 21:26

dennisYatunin marked this pull request as ready for review July 31, 2023 21:26

bors bot added a commit that referenced this pull request Jul 31, 2023

Try #1399:

51b75c3

dennisYatunin force-pushed the dy/operator_matrices branch from 403cb67 to 4b747af Compare July 31, 2023 23:22

charleskawczynski reviewed Jul 31, 2023

View reviewed changes

src/MatrixFields/lazy_operators.jl Outdated Show resolved Hide resolved

dennisYatunin force-pushed the dy/operator_matrices branch 2 times, most recently from 5d5aedc to 3003cb4 Compare August 9, 2023 20:20

dennisYatunin mentioned this pull request Aug 9, 2023

Clean up MatrixFields tests #1423

Merged

4 tasks

dennisYatunin force-pushed the dy/operator_matrices branch from 3003cb4 to a95e919 Compare August 9, 2023 23:34

dennisYatunin mentioned this pull request Aug 10, 2023

Fix matrix multiplication at boundaries #1424

Merged

dennisYatunin force-pushed the dy/operator_matrices branch 3 times, most recently from 8b40c6b to 70eba2e Compare August 11, 2023 18:48

bors bot added a commit that referenced this pull request Aug 11, 2023

Try #1399:

72b31dc

dennisYatunin requested a review from charleskawczynski August 11, 2023 18:57

bors bot added a commit that referenced this pull request Aug 11, 2023

Try #1399:

f01b0df

Add operator_matrix to MatrixFields, along with tests and docs

d284e69

dennisYatunin force-pushed the dy/operator_matrices branch from ce23cf0 to d284e69 Compare August 14, 2023 16:11

bors bot merged commit ca3864c into main Aug 14, 2023

bors bot deleted the dy/operator_matrices branch August 14, 2023 17:44

dennisYatunin mentioned this pull request Oct 17, 2023

Update to new implicit solver interface CliMA/ClimaAtmos.jl#2044

Merged

1 task

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add operator_matrix to MatrixFields, along with tests and docs #1399

Add operator_matrix to MatrixFields, along with tests and docs #1399

dennisYatunin commented Jul 28, 2023 •

edited

Loading

dennisYatunin commented Jul 28, 2023

bors bot commented Jul 28, 2023

dennisYatunin commented Jul 31, 2023

bors bot commented Jul 31, 2023

dennisYatunin commented Jul 31, 2023

bors bot commented Jul 31, 2023

charleskawczynski Jul 31, 2023

dennisYatunin commented Jul 31, 2023

bors bot commented Jul 31, 2023

dennisYatunin commented Aug 9, 2023

dennisYatunin commented Aug 11, 2023

bors bot commented Aug 11, 2023

dennisYatunin commented Aug 11, 2023

bors bot commented Aug 12, 2023

dennisYatunin commented Aug 14, 2023

dennisYatunin commented Aug 14, 2023

bors bot commented Aug 14, 2023

Add operator_matrix to MatrixFields, along with tests and docs #1399

Add operator_matrix to MatrixFields, along with tests and docs #1399

Conversation

dennisYatunin commented Jul 28, 2023 • edited Loading

Purpose

Content

dennisYatunin commented Jul 28, 2023

bors bot commented Jul 28, 2023

try

dennisYatunin commented Jul 31, 2023

bors bot commented Jul 31, 2023

try

dennisYatunin commented Jul 31, 2023

bors bot commented Jul 31, 2023

try

charleskawczynski Jul 31, 2023

Choose a reason for hiding this comment

dennisYatunin commented Jul 31, 2023

bors bot commented Jul 31, 2023

try

dennisYatunin commented Aug 9, 2023

dennisYatunin commented Aug 11, 2023

bors bot commented Aug 11, 2023

try

dennisYatunin commented Aug 11, 2023

bors bot commented Aug 12, 2023

try

dennisYatunin commented Aug 14, 2023

dennisYatunin commented Aug 14, 2023

bors bot commented Aug 14, 2023

dennisYatunin commented Jul 28, 2023 •

edited

Loading