Remove allocations in the spline_interpolation function (#60)

* make the 1d spline_inteprolation functions allocation free

* reduce allocations for the 2D spline_interpolation routines

* small updates to the docs

* apply formatter

* fix wrong matrix for bilinear interpolation

Author: andrewwinters5000

Project.toml

@@ -7,6 +7,7 @@ version = "0.1.2-DEV"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [weakdeps]
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
@@ -19,4 +20,5 @@ Downloads = "1.6"
 LinearAlgebra = "1"
 Makie = "0.21, 0.22"
 SparseArrays = "1"
+StaticArrays = "1.9"
 julia = "1.10"

docs/src/trixi_jl_examples.md

@@ -52,9 +52,10 @@ and [B-spline function](https://trixi-framework.github.io/TrixiBottomTopography.
 In this case, a cubic B-spline interpolation function with free end condition is chosen.
 
 ```@example trixi1d
-# B-spline interpolation of the underlying data
-spline_struct = CubicBSpline(Rhine_data)
-spline_func(x) = spline_interpolation(spline_struct, x)
+# B-spline interpolation of the underlying data.
+# The type of this struct is fixed as `CubicBSpline`.
+const spline_struct = CubicBSpline(Rhine_data)
+spline_func(x::Float64) = spline_interpolation(spline_struct, x)
 ```
 
 Now that the B-spline interpolation function is determined, the one dimensional shallow water equations implemented in Trixi.jl can be defined by calling:
@@ -242,9 +243,10 @@ nothing #hide
 Using the data, a bicubic B-spline interpolation is performed on the data to define a bottom topography function.
 
 ```@example trixi2d
-# B-spline interpolation of the underlying data
-spline_struct = BicubicBSpline(Rhine_data)
-spline_func(x,y) = spline_interpolation(spline_struct, x, y)
+# B-spline interpolation of the underlying data.
+# The type of this struct is fixed as `BicubicBSpline`.
+const spline_struct = BicubicBSpline(Rhine_data)
+spline_func(x::Float64, y::Float64) = spline_interpolation(spline_struct, x, y)
 ```
 
 Then the two dimensional shallow water equations are defined, where the gravitational constant has been chosen to be `3.0` and the initial water height `55.0`. Afterwards, the initial condition is defined. Similar to the one dimensional case, in the center of the domain, a circular part with a diameter of `100.0` is chosen where the initial water height is chosen to be `10.0` units higher.

src/1D/spline_cache_1D.jl

@@ -74,8 +74,8 @@ function LinearBSpline(x::Vector, y::Vector)
 
     Delta = x[2] - x[1]
 
-    IP = [-1 1;
-          1 0]
+    IP = @SMatrix [-1 1;
+                   1 0]
 
     Q = y
 
@@ -256,10 +256,10 @@ function CubicBSpline(x::Vector, y::Vector; end_condition = "free", smoothing_fa
 
     n = length(x)
     P = vcat(0, y, 0)
-    IP = [-1 3 -3 1;
-          3 -6 3 0;
-          -3 0 3 0;
-          1 4 1 0]
+    IP = @SMatrix [-1 3 -3 1;
+                   3 -6 3 0;
+                   -3 0 3 0;
+                   1 4 1 0]
 
     # Free end condition
     if end_condition == "free"

src/1D/spline_methods_1D.jl

@@ -48,7 +48,15 @@ function spline_interpolation(b_spline::LinearBSpline, x::Number)
 
     kappa_i = (x - x_vec[i]) / Delta
 
-    c_i1 = [kappa_i, 1]' * IP * Q[i:(i + 1)]
+    # Allocation free version of the naive implementation
+    # c_i1 = [kappa_i, 1]' * IP * Q[i:(i + 1)]
+
+    # Note, IP has already been constructed as an `SMatrix`
+    # in the `LinearBSpline` constructor
+    kappa_vec = @SVector [kappa_i, 1]
+    Q_slice = @SVector [Q[i], Q[i + 1]]
+
+    c_i1 = kappa_vec' * IP * Q_slice
 
     return c_i1
 end
@@ -101,7 +109,15 @@ function spline_interpolation(b_spline::CubicBSpline, x::Number)
 
     kappa_i = (x - x_vec[i]) / Delta
 
-    c_i3 = 1 / 6 * [kappa_i^3, kappa_i^2, kappa_i, 1]' * IP * Q[i:(i + 3)]
+    # Allocation free version of the naive implementation
+    # c_i3 = 1 / 6 * [kappa_i^3, kappa_i^2, kappa_i, 1]' * IP * Q[i:(i + 3)]
+
+    # Note, IP has already been constructed as an `SMatrix`
+    # in the `CubicBSpline` constructor
+    kappa_vec = @SVector [kappa_i^3, kappa_i^2, kappa_i, 1]
+    Q_slice = @SVector [Q[i], Q[i + 1], Q[i + 2], Q[i + 3]]
+
+    c_i3 = (1 / 6) * (kappa_vec' * IP * Q_slice)
 
     return c_i3
 end

src/2D/spline_cache_2D.jl

@@ -97,8 +97,8 @@ function BilinearBSpline(x::Vector, y::Vector, z::Matrix)
     Delta = x[2] - x[1]
 
     P = vcat(reshape(z', (m * n, 1)))
-    IP = [-1 1;
-          1 0]
+    IP = @SMatrix [-1 1;
+                   1 0]
 
     Q = reshape(P, (n, m))
 
@@ -300,10 +300,10 @@ function BicubicBSpline(x::Vector, y::Vector, z::Matrix; end_condition = "free",
     inner_elmts = m * n
     P = vcat(reshape(z', (inner_elmts, 1)), zeros(boundary_elmts))
 
-    IP = [-1 3 -3 1;
-          3 -6 3 0;
-          -3 0 3 0;
-          1 4 1 0]
+    IP = @SMatrix [-1 3 -3 1;
+                   3 -6 3 0;
+                   -3 0 3 0;
+                   1 4 1 0]
 
     # Mapping matrix Phi
     Phi = spzeros((m + 2) * (n + 2), (m + 2) * (n + 2))

src/2D/spline_methods_2D.jl

@@ -47,9 +47,18 @@ function spline_interpolation(b_spline::BilinearBSpline, x::Number, y::Number)
     my = (x - x_vec[i]) / Delta
     ny = (y - y_vec[j]) / Delta
 
-    Q_temp = [Q[i, j:(j + 1)] Q[(i + 1), j:(j + 1)]]
+    # Allocation free version of the naive implementation
+    # Q_temp = [Q[i, j:(j + 1)] Q[(i + 1), j:(j + 1)]]
+    # c = [ny, 1]' * IP * Q_temp * IP' * [my, 1]
 
-    c = [ny, 1]' * IP * Q_temp * IP' * [my, 1]
+    # Note, IP has already been constructed as an `SMatrix`
+    # in the `BilinearBSpline` constructor
+    ny_vec = @SVector [ny, 1]
+    my_vec = @SVector [my, 1]
+    Q_temp = @SMatrix [Q[i, j] Q[i + 1, j];
+                       Q[i, j + 1] Q[i + 1, j + 1]]
+
+    c = ny_vec' * IP * Q_temp * IP' * my_vec
 
     return c
 end
@@ -113,10 +122,20 @@ function spline_interpolation(b_spline::BicubicBSpline, x::Number, y::Number)
     my = (x - x_vec[i]) / Delta
     ny = (y - y_vec[j]) / Delta
 
-    Q_temp = [Q[i, j:(j + 3)] Q[(i + 1), j:(j + 3)] Q[(i + 2), j:(j + 3)] Q[(i + 3),
-                                                                            j:(j + 3)]]
+    # Allocation free version of the naive implementation
+    # Q_temp = [Q[i, j:(j + 3)] Q[(i + 1), j:(j + 3)] Q[(i + 2), j:(j + 3)] Q[(i + 3), j:(j + 3)]]
+    # c = 1 / 36 * [ny^3, ny^2, ny, 1]' * IP * Q_temp * IP' * [my^3, my^2, my, 1]
+
+    # Note, IP has already been constructed as an `SMatrix`
+    # in the `BicubicBSpline` constructor
+    ny_vec = @SVector [ny^3, ny^2, ny, 1]
+    my_vec = @SVector [my^3, my^2, my, 1]
+    Q_temp = @SMatrix [Q[i, j] Q[i + 1, j] Q[i + 2, j] Q[i + 3, j];
+                       Q[i, j + 1] Q[i + 1, j + 1] Q[i + 2, j + 1] Q[i + 3, j + 1];
+                       Q[i, j + 2] Q[i + 1, j + 2] Q[i + 2, j + 2] Q[i + 3, j + 2];
+                       Q[i, j + 3] Q[i + 1, j + 3] Q[i + 2, j + 3] Q[i + 3, j + 3]]
 
-    c = 1 / 36 * [ny^3, ny^2, ny, 1]' * IP * Q_temp * IP' * [my^3, my^2, my, 1]
+    c = (1 / 36) * (ny_vec' * IP * Q_temp * IP' * my_vec)
 
     return c
 end

src/TrixiBottomTopography.jl

@@ -10,6 +10,7 @@ module TrixiBottomTopography
 # Include necessary packages
 using LinearAlgebra: norm, diagm, qr, Tridiagonal, SymTridiagonal
 using SparseArrays: sparse, spzeros
+using StaticArrays: SVector, @SVector, SMatrix, @SMatrix
 
 # Include one dimensional B-spline interpolation
 include("1D/spline_cache_1D.jl")