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.

Sign up for GitHub

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

Jump to bottom

Run LowerTriangularMatrix and output examples while building the docs #453

Merged
merged 1 commit into from
Feb 22, 2024
Merged

Run LowerTriangularMatrix and output examples while building the docs #453

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
103 changes: 34 additions & 69 deletions docs/src/lowertriangularmatrices.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,117 +11,82 @@ the spherical harmonic coefficients (see [Spectral packing](@ref)).

## Creation of `LowerTriangularMatrix`

A `LowerTriangularMatrix` can be created using `zeros`,`ones`,`rand`, or `randn`
```julia
julia> using SpeedyWeather.LowerTriangularMatrices

julia> L = rand(LowerTriangularMatrix{Float32},5,5)
5×5 LowerTriangularMatrix{Float32}:
0.912744 0.0 0.0 0.0 0.0
0.0737592 0.230592 0.0 0.0 0.0
0.799679 0.0765255 0.888098 0.0 0.0
0.670835 0.997938 0.505276 0.492966 0.0
0.949321 0.193692 0.793623 0.152817 0.357968
A `LowerTriangularMatrix` can be created using `zeros`, `ones`, `rand`, or `randn`
```@repl LowerTriangularMatrices
using SpeedyWeather.LowerTriangularMatrices

L = rand(LowerTriangularMatrix{Float32},5,5)
```
or the `undef` initializor `LowerTriangularMatrix{Float32}(undef,3,3)`.
The element type is arbitrary though, you can use any type `T` too.

Alternatively, it can be created through conversion from `Matrix`, which drops the upper triangle
entries and sets them to zero
```julia
julia> M = rand(Float16,3,3)
3×3 Matrix{Float16}:
0.2222 0.694 0.3452
0.2158 0.04443 0.274
0.9746 0.793 0.6294

julia> LowerTriangularMatrix(M)
3×3 LowerTriangularMatrix{Float16}:
0.2222 0.0 0.0
0.2158 0.04443 0.0
0.9746 0.793 0.6294
```@repl LowerTriangularMatrices
M = rand(Float16,3,3)

L = LowerTriangularMatrix(M)
```

## Indexing `LowerTriangularMatrix`

`LowerTriangularMatrix` supports two types of indexing: 1) by denoting two indices, column and row `[l,m]`
or 2) by denoting a single index `[lm]`. The double index works as expected
```julia
julia> L
3×3 LowerTriangularMatrix{Float16}:
0.1499 0.0 0.0
0.1177 0.478 0.0
0.1709 0.756 0.3223

julia> L[2,2]
Float16(0.478)
```@repl LowerTriangularMatrices
L

L[2,2]
```
But the single index skips the zero entries in the upper triangle, i.e.
```julia
julia> L[4]
Float16(0.478)
```@repl LowerTriangularMatrices
L[4]
```
which, important, is different from single indices of an `AbstractMatrix`
```julia
julia> Matrix(L)[4]
Float16(0.0)
```@repl LowerTriangularMatrices
Matrix(L)[4]
```
In performance-critical code a single index should be used, as this directly maps
to the index of the underlying data vector. The double index is somewhat slower
as it first has to be converted to the corresponding single index.

Consequently, many loops in SpeedyWeather.jl are build with the following structure
```julia
```@repl LowerTriangularMatrices
n,m = size(L)

ij = 0
for j in 1:m
for i in j:n
ij += 1
L[ij] = i+j
end

for j in 1:m, i in j:n
ij += 1
L[ij] = i+j
end
```
which loops over all lower triangle entries of `L::LowerTriangularMatrix` and the single
index `ij` is simply counted up. However, one could also use `[i,j]` as indices in the
loop body or to perform any calculation (`i+j` here).
An iterator over all entries in the lower triangle can be created by
```julia
```@repl LowerTriangularMatrices
for ij in eachindex(L)
# do something
end
```
The `setindex!` functionality of matrixes will throw a `BoundsError` when trying to write
into the upper triangle of a `LowerTriangularMatrix`, for example
```julia
julia> L[2,1] = 0 # valid index
0
```@repl LowerTriangularMatrices
L[2,1] = 0 # valid index

julia> L[1,2] = 0 # invalid index in the upper triangle
ERROR: BoundsError: attempt to access 3×3 LowerTriangularMatrix{Float32} at index [1, 2]
L[1,2] = 0 # invalid index in the upper triangle
```

## Linear algebra with `LowerTriangularMatrix`

The [LowerTriangularMatrices](@ref lowertriangularmatrices) module's main purpose is not linear algebra, and it's
implementation may not be efficient, however, many operations work as expected
```julia
julia> L = rand(LowerTriangularMatrix{Float32},3,3)
3×3 LowerTriangularMatrix{Float32}:
0.57649 0.0 0.0
0.348685 0.875371 0.0
0.881923 0.850552 0.998306

julia> L + L
3×3 LowerTriangularMatrix{Float32}:
1.15298 0.0 0.0
0.697371 1.75074 0.0
1.76385 1.7011 1.99661

julia> L * L
3×3 Matrix{Float32}:
0.332341 0.0 0.0
0.506243 0.766275 0.0
1.68542 1.59366 0.996616
```@repl LowerTriangularMatrices
L = rand(LowerTriangularMatrix{Float32},3,3)

L + L

L * L
```
Note, however, that the latter includes a conversion to `Matrix`, which is true for many
operations, including `inv` or `\`. Hence when trying to do more sophisticated linear
Expand All @@ -132,4 +97,4 @@ back to normal matrix-land.

```@autodocs
Modules = [SpeedyWeather.LowerTriangularMatrices]
```
```
11 changes: 5 additions & 6 deletions docs/src/output.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,14 +94,13 @@ So for example `output_Grid=FullClenshawGrid` and `nlat_half=48` will always int
regular 192x95 longitude-latitude grid of 1.875˚ resolution, regardless the grid and resolution used
for the model integration.
```julia
julia> my_output_writer = OutputWriter(spectral_grid, ShallowWater, output_Grid=FullClenshawGrid, nlat_half=48)
my_output_writer = OutputWriter(spectral_grid, ShallowWater, output_Grid=FullClenshawGrid, nlat_half=48)
```
Note that by default the output is on the corresponding full of the grid used in the dynamical core
so that interpolation only happens at most in the zonal direction as they share the location of the
latitude rings. You can check this by
```julia
julia> RingGrids.full_grid(OctahedralGaussianGrid)
FullGaussianGrid
```@example netcdf
RingGrids.full_grid(OctahedralGaussianGrid)
```
So the corresponding full grid of an `OctahedralGaussianGrid` is the `FullGaussiangrid` and the same resolution
`nlat_half` is chosen by default in the output writer (which you can change though as shown above).
Expand Down Expand Up @@ -131,12 +130,12 @@ julia> path = pwd()
julia> my_output_writer = OutputWriter(spectral_grid, PrimitiveDry, path=path)
```
This folder must already exist. If you want to give your run a name/identification you can pass on `id`
```
```julia
julia> my_output_writer = OutputWriter(spectral_grid,PrimitiveDry,id="diffusion_test");
```
which will be used instead of a 4 digit number like 0001, 0002 which is automatically determined if
`id` is not provided. You will see the id of the run in the progress bar
```
```julia
Weather is speedy: run diffusion_test 100%|███████████████████████| Time: 0:00:12 (19.20 years/day)
```
and the run folder, here `run_diffusion_test`, is also named accordingly
Expand Down
Loading