tests/implementation for same number of bpc, still needs 3rd dimension handling

Author: ajozefiak

src/derivative_operators/derivative_operator_functions.jl

@@ -413,54 +413,68 @@ function LinearAlgebra.mul!(x_temp::AbstractArray{T,3}, A::AbstractDiffEqComposi
 
         # convolve boundary and interior points near boundary
         # partition operator indices along axis of differentiation
-        if pad[1] > 0 || pad[2] > 0
+        if pad[1] > 0 || pad[2] > 0 || pad[3] > 0
             ops_1 = Int64[]
             ops_1_max_bpc_idx = [0]
             ops_2 = Int64[]
             ops_2_max_bpc_idx = [0]
+            ops_3 = Int64[]
+            ops_3_max_bpc_idx = [0]
+
             for i in 1:length(opsA)
                 L = opsA[i]
                 if typeof(L).parameters[2] == 1
                     push!(ops_1,i)
                     if L.boundary_point_count == pad[1]
                         ops_1_max_bpc_idx[1] = i
                     end
-                else
+                elseif typeof(L).parameters[2] == 2
                     push!(ops_2,i)
                     if L.boundary_point_count == pad[2]
                         ops_2_max_bpc_idx[1]= i
                     end
+                else
+                    push!(ops_3,i)
+                    if L.boundary_point_count == pad[3]
+                        ops_3_max_bpc_idx[1]= i
+                    end
                 end
             end
 
             # need offsets since some axis may have ghost nodes and some may not
             offset_x = 0
             offset_y = 0
+            offset_z = 0
 
-            if length(ops_2) > 0
+            if length(ops_1) > 0
                 offset_x = 1
             end
-            if length(ops_1) > 0
+            if length(ops_2) > 0
                 offset_y = 1
             end
+            if length(ops_3) > 0
+                offset_z = 1
+            end
 
             # convolve boundaries and unaccounted for interior in axis 1
             if length(ops_1) > 0
                 for i in 1:size(x_temp)[2]
-                    convolve_BC_left!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[ops_1_max_bpc_idx...])
-                    convolve_BC_right!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[ops_1_max_bpc_idx...])
-                    if i <= pad[2] || i > size(x_temp)[2]-pad[2]
-                        convolve_interior!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[ops_1_max_bpc_idx...])
-                    end
+                    for j in 1:size(x_temp)[3]
+                        convolve_BC_left!(view(x_temp,:,i,j), view(M,:,i+offset_y,j+offset_z), opsA[ops_1_max_bpc_idx...])
+                        convolve_BC_right!(view(x_temp,:,i,j), view(M,:,i+offset_y,j+offset_z), opsA[ops_1_max_bpc_idx...])
+                        if i <= pad[2] || i > size(x_temp)[2]-pad[2] || j <= pad[3] || j > size(x_temp)[3]-pad[3]
+                            convolve_interior!(view(x_temp,:,i,j), view(M,:,i+offset_y,j+offset_z), opsA[ops_1_max_bpc_idx...])
+                        end
 
-                    for Lidx in ops_1
-                        if Lidx != ops_1_max_bpc_idx[1]
-                            convolve_BC_left!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[Lidx], overwrite = false)
-                            convolve_BC_right!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[Lidx], overwrite = false)
-                            if i <= pad[2] || i > size(x_temp)[2]-pad[2]
-                                convolve_interior!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[Lidx], overwrite = false)
-                            elseif pad[1] - opsA[Lidx].boundary_point_count > 0
-                                convolve_interior_add_range!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[Lidx], pad[1] - opsA[Lidx].boundary_point_count)
+                        for Lidx in ops_1
+                            if Lidx != ops_1_max_bpc_idx[1]
+                                convolve_BC_left!(view(x_temp,:,i,j), view(M,:,i+offset_y,j+offset_z), opsA[Lidx], overwrite = false)
+                                convolve_BC_right!(view(x_temp,:,i,j), view(M,:,i+offset_y,j+offset_z), opsA[Lidx], overwrite = false)
+                                if i <= pad[2] || i > size(x_temp)[2]-pad[2] || j <= pad[3] || j > size(x_temp)[3]-pad[3]
+                                    convolve_interior!(view(x_temp,:,i,j), view(M,:,i+offset_y,j+offset_z), opsA[Lidx], overwrite = false)
+                                elseif pad[1] - opsA[Lidx].boundary_point_count > 0 # TODO 475-477
+                                    convolve_interior_add_range!(view(x_temp,:,i), view(M,:,i+offset_x), opsA[Lidx], pad[1] - opsA[Lidx].boundary_point_count)
+                                end
                             end
                         end
                     end
@@ -469,30 +483,32 @@ function LinearAlgebra.mul!(x_temp::AbstractArray{T,3}, A::AbstractDiffEqComposi
             # convolve boundaries and unaccounted for interior in axis 2
             if length(ops_2) > 0
                 for i in 1:size(x_temp)[1]
-                    # in the case of no axis 1 operators, we need to overwrite x_temp
-                    if length(ops_1) == 0
-                        convolve_BC_left!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[ops_2_max_bpc_idx...])
-                        convolve_BC_right!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[ops_2_max_bpc_idx...])
-                        if i <= pad[1] || i > size(x_temp)[1]-pad[1]
-                            convolve_interior!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[ops_2_max_bpc_idx...])
-                        end
+                    for j in 1:size(x_temp)[3]
+                        # in the case of no axis 1 operators, we need to overwrite x_temp
+                        if length(ops_1) == 0
+                            convolve_BC_left!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[ops_2_max_bpc_idx...])
+                            convolve_BC_right!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[ops_2_max_bpc_idx...])
+                            if i <= pad[1] || i > size(x_temp)[1]-pad[1] #TODO 491-493
+                                convolve_interior!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[ops_2_max_bpc_idx...])
+                            end
 
-                    else
-                        convolve_BC_left!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[ops_2_max_bpc_idx...], overwrite = false)
-                        convolve_BC_right!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[ops_2_max_bpc_idx...], overwrite = false)
-                        if i <= pad[1] || i > size(x_temp)[1]-pad[1]
-                            convolve_interior!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[ops_2_max_bpc_idx...], overwrite = false)
-                        end
+                        else
+                            convolve_BC_left!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[ops_2_max_bpc_idx...], overwrite = false)
+                            convolve_BC_right!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[ops_2_max_bpc_idx...], overwrite = false)
+                            if i <= pad[1] || i > size(x_temp)[1]-pad[1] || j <= pad[3] || j > size(x_temp)[3]-pad[3]
+                                convolve_interior!(view(x_temp,i,:,j), view(M,i+offset_y,:,j+offset_z), opsA[ops_2_max_bpc_idx...], overwrite = false)
+                            end
 
-                    end
-                    for Lidx in ops_2
-                        if Lidx != ops_2_max_bpc_idx[1]
-                            convolve_BC_left!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[Lidx], overwrite = false)
-                            convolve_BC_right!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[Lidx], overwrite = false)
-                            if i <= pad[1] || i > size(x_temp)[1]-pad[1]
-                                convolve_interior!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[Lidx], overwrite = false)
-                            elseif pad[2] - opsA[Lidx].boundary_point_count > 0
-                                convolve_interior_add_range!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[Lidx], pad[2] - opsA[Lidx].boundary_point_count)
+                        end
+                        for Lidx in ops_2
+                            if Lidx != ops_2_max_bpc_idx[1]
+                                convolve_BC_left!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[Lidx], overwrite = false)
+                                convolve_BC_right!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[Lidx], overwrite = false)
+                                if i <= pad[1] || i > size(x_temp)[1]-pad[1] || j <= pad[3] || j > size(x_temp)[3]-pad[3]
+                                    convolve_interior!(view(x_temp,i,:,j), view(M,i+offset_x,:,j+offset_z), opsA[Lidx], overwrite = false)
+                                elseif pad[2] - opsA[Lidx].boundary_point_count > 0 #TODO 509-511
+                                    convolve_interior_add_range!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[Lidx], pad[2] - opsA[Lidx].boundary_point_count)
+                                end
                             end
                         end
                     end
@@ -502,7 +518,7 @@ function LinearAlgebra.mul!(x_temp::AbstractArray{T,3}, A::AbstractDiffEqComposi
 
         # Here we compute mul! (additively) for every operator in opsB
 
-        operating_dims = zeros(Int64,2)
+        operating_dims = zeros(Int64,3)
         # need to consider all dimensions and operators to determine the truncation
         # of M to x_temp
         for L in A.ops

test/2D_3D_fast_multiplication.jl

@@ -712,3 +712,75 @@ end
     @test M_temp ≈ ((Lx*M)[1:N,2:N+1,2:N+1] +(Ly*M)[2:N+1,1:N,2:N+1] + (Lxx*M)[1:N,2:N+1,2:N+1] +(Lyy*M)[2:N+1,1:N,2:N+1] + (Lz*M)[2:N+1,2:N+1,1:N] +(Lzz*M)[2:N+1,2:N+1,1:N])
 
 end
+
+@testset "3D Multiplication with identical bpc and dx = dy = dz = 1.0" begin
+
+    N = 100
+    M = zeros(N+2,N+2,N+2)
+    for i in 1:N+2
+        for j in 1:N+2
+            for k in 1:N+2
+                M[i,j,k] = cos(0.1i)+sin(0.1j) + exp(0.01k)
+            end
+        end
+    end
+
+    # Test a single axis, multiple operators: (Lxxx + Lxxxx)*M, dx = dy = dz = 1.0
+    # Lx3 and Lx4 have the same number of boundary points
+    Lx3 = CenteredDifference{1}(3,4,1.0,N)
+    Lx4 = CenteredDifference{1}(4,4,1.0,N)
+    A = Lx3 + Lx4
+
+    M_temp = zeros(N,N+2,N+2)
+    mul!(M_temp, A, M)
+
+    @test M_temp ≈ ((Lx3*M)+(Lx4*M))
+
+    # Test a single axis, multiple operators: (Lyyy + Lyyyy)*M, dx = dy = dz = 1.0, no coefficient
+    # Ly3 and Ly4 have the same number of boundary points
+    Ly3 = CenteredDifference{2}(3,4,1.0,N)
+    Ly4 = CenteredDifference{2}(4,4,1.0,N)
+    A = Ly3 + Ly4
+
+    M_temp = zeros(N+2,N,N+2)
+    mul!(M_temp, A, M)
+
+    @test M_temp ≈ ((Ly3*M)+(Ly4*M))
+
+    # Test a single axis, multiple operators: (Lzzz + Lzzzz)*M, dx = dy = dz = 1.0, no coefficient
+    # Lz3 and Lz4 have the same number of boundary points
+    Lz3 = CenteredDifference{3}(3,4,1.0,N)
+    Lz4 = CenteredDifference{3}(4,4,1.0,N)
+    A = Lz3 + Lz4
+
+    M_temp = zeros(N+2,N+2,N)
+    @test_broken mul!(M_temp, A, M)
+
+    @test_broken M_temp ≈ ((Lz3*M)+(Lz4*M))
+
+    # Test (Lxxxx + Lyyyy)*M, dx = 1.0, no coefficient, two boundary points on each axis
+
+    M_temp = zeros(N,N,N+2)
+    A = Lx4 + Ly4
+    mul!(M_temp, A, M)
+
+    @test M_temp ≈ ((Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])
+
+    # Test (Lxxx + Lyyy)*M, no coefficient. These operators have non-symmetric interior stencils
+    A = Lx3 + Ly3
+    M_temp = zeros(N,N,N+2)
+    mul!(M_temp, A, M)
+
+    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:])
+
+    # Test multiple operators on both axis: (Lxxx + Lyyy + Lxxxx + Lyyyy)*M, no coefficient
+    A = Lx3 + Ly3 + Lx4 + Ly4
+    M_temp = zeros(N,N,N+2)
+    mul!(M_temp, A, M)
+
+    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])
+
+    # TODO implement/test all combinations of x-z and y-z compositions.
+    # TODO implement/test combinations of x-y-z compositions
+
+end