tests for same bpc fleshed out and all are passing

Author: ajozefiak

src/derivative_operators/derivative_operator_functions.jl

@@ -488,8 +488,8 @@ function LinearAlgebra.mul!(x_temp::AbstractArray{T,3}, A::AbstractDiffEqComposi
                         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...])
+                            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[ops_2_max_bpc_idx...])
                             end
 
                         else
@@ -514,6 +514,40 @@ function LinearAlgebra.mul!(x_temp::AbstractArray{T,3}, A::AbstractDiffEqComposi
                     end
                 end
             end
+            # convolve boundaries and unaccounted for interior in axis 3
+            if length(ops_3) > 0
+                for i in 1:size(x_temp)[1]
+                    for j in 1:size(x_temp)[2]
+                        # in the case of no axis 1 and 2 operators, we need to overwrite x_temp
+                        if length(ops_1) == 0 && length(ops_2) == 0
+                            convolve_BC_left!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[ops_3_max_bpc_idx...])
+                            convolve_BC_right!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[ops_3_max_bpc_idx...])
+                            if i <= pad[1] || i > size(x_temp)[1]-pad[1] #TODO 525-527
+                                convolve_interior!(view(x_temp,i,:), view(M,i+offset_y,:), opsA[ops_3_max_bpc_idx...])
+                            end
+
+                        else
+                            convolve_BC_left!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[ops_3_max_bpc_idx...], overwrite = false)
+                            convolve_BC_right!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[ops_3_max_bpc_idx...], overwrite = false)
+                            if i <= pad[1] || i > size(x_temp)[1]-pad[1] || j <= pad[2] || j > size(x_temp)[2]-pad[2]
+                                convolve_interior!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[ops_3_max_bpc_idx...], overwrite = false)
+                            end
+
+                        end
+                        for Lidx in ops_3
+                            if Lidx != ops_3_max_bpc_idx[1]
+                                convolve_BC_left!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[Lidx], overwrite = false)
+                                convolve_BC_right!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[Lidx], overwrite = false)
+                                if i <= pad[1] || i > size(x_temp)[1]-pad[1] || j <= pad[2] || j > size(x_temp)[2]-pad[2]
+                                    convolve_interior!(view(x_temp,i,j,:), view(M,i+offset_x,j+offset_y,:), opsA[Lidx], overwrite = false)
+                                elseif pad[3] - opsA[Lidx].boundary_point_count > 0 #TODO 543-545
+                                    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
+                end
+            end
         end
 
         # Here we compute mul! (additively) for every operator in opsB

test/2D_3D_fast_multiplication.jl

@@ -754,12 +754,11 @@ end
     A = Lz3 + Lz4
 
     M_temp = zeros(N+2,N+2,N)
-    @test_broken mul!(M_temp, A, M)
+    mul!(M_temp, A, M)
 
-    @test_broken M_temp ≈ ((Lz3*M)+(Lz4*M))
+    @test 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)
@@ -780,7 +779,67 @@ end
 
     @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
+    # Test (Lxxxx + Lzzzz)*M, dx = 1.0, no coefficient, two boundary points on each axis
+    M_temp = zeros(N,N+2,N)
+    A = Lx4 + Lz4
+    mul!(M_temp, A, M)
+
+    @test M_temp ≈ ((Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])
+
+    # Test (Lxxx + Lzzz)*M, no coefficient. These operators have non-symmetric interior stencils
+    A = Lx3 + Lz3
+    M_temp = zeros(N,N+2,N)
+    mul!(M_temp, A, M)
+
+    @test M_temp ≈ ((Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N])
+
+    # Test multiple operators on both axis: (Lxxx + Lzzz + Lxxxx + Lzzzz)*M, no coefficient
+    A = Lx3 + Lz3 + Lx4 + Lz4
+    M_temp = zeros(N,N+2,N)
+    mul!(M_temp, A, M)
+
+    @test M_temp ≈ ((Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N] + (Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])
+
+    # Test (Lyyyy + Lzzzz)*M, dx = 1.0, no coefficient, two boundary points on each axis
+    M_temp = zeros(N+2,N,N)
+    A = Ly4 + Lz4
+    mul!(M_temp, A, M)
+
+    @test M_temp ≈ ((Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])
+
+    # Test (Lyyy + Lzzz)*M, no coefficient. These operators have non-symmetric interior stencils
+    A = Ly3 + Lz3
+    M_temp = zeros(N+2,N,N)
+    mul!(M_temp, A, M)
+
+    @test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N])
+
+    # Test multiple operators on both axis: (Lyyy + Lzzz + Lyyyy + Lzzzz)*M, no coefficient
+    A = Ly3 + Lz3 + Ly4 + Lz4
+    M_temp = zeros(N+2,N,N)
+    mul!(M_temp, A, M)
+
+    @test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N] + (Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])
+
+    # Test a single operator on each axis: (Lx3 + Ly3 + Lz3)*M, no coefficient
+    A = Lx3 + Ly3 + Lz3
+    M_temp = zeros(N,N,N)
+    mul!(M_temp, A, M)
+
+    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1,2:N+1] + (Ly3*M)[2:N+1,1:N,2:N+1] +(Lz3*M)[2:N+1,2:N+1,1:N])
+
+    # Test a single operator on each axis: (Lx4 + Ly4 + Lz4)*M, no coefficient
+    A = Lx4 + Ly4 + Lz4
+    M_temp = zeros(N,N,N)
+    mul!(M_temp, A, M)
+
+    @test M_temp ≈ ((Lx4*M)[1:N,2:N+1,2:N+1] + (Ly4*M)[2:N+1,1:N,2:N+1] +(Lz4*M)[2:N+1,2:N+1,1:N])
+
+    # Test multiple operators on each axis: (Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4)*M, no coefficient
+    A = Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4
+    M_temp = zeros(N,N,N)
+    mul!(M_temp, A, M)
+
+    @test M_temp ≈ ((Lx3*M)[1:N,2:Ns+1,2:N+1] + (Ly3*M)[2:N+1,1:N,2:N+1] +(Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] + (Ly4*M)[2:N+1,1:N,2:N+1] +(Lz4*M)[2:N+1,2:N+1,1:N])
 
 end