add derivatives switch

src/FastMultipole.jl

@@ -17,6 +17,6 @@ for file in ["containers", "complex", "derivatives", "element", "tree", "direct"
     include(file*".jl")
 end
 
-export fmm!, direct!, Tree, SortWrapper, ProbeSystem, reset!
+export fmm!, direct!, Tree, SortWrapper, ProbeSystem, reset!, DerivativesSwitch
 
 end # module

src/b2m.jl

@@ -10,8 +10,8 @@
         z = b_z - c_z
         y = b_y - c_y
         qx, qy, qz = system[i_body,VECTOR_STRENGTH]
-        r, multipole_acceptance_criterion, phi = cartesian_2_spherical(x,y,z)
-        regular_harmonic!(harmonics, r ,multipole_acceptance_criterion, -phi, P) # Ylm^* -> -dx[3]
+        r, theta, phi = cartesian_2_spherical(x,y,z)
+        regular_harmonic!(harmonics, r ,theta, -phi, P) # Ylm^* -> -dx[3]
         # update values
         for l in 0:P
             for m in 0:l
@@ -39,8 +39,8 @@ end
         z = b_z - c_z
         y = b_y - c_y
         q = system[i_body,SCALAR_STRENGTH]
-        r, multipole_acceptance_criterion, phi = cartesian_2_spherical(x,y,z)
-        regular_harmonic!(harmonics, r, multipole_acceptance_criterion, -phi, P) # Ylm^* -> -dx[3]
+        r, theta, phi = cartesian_2_spherical(x,y,z)
+        regular_harmonic!(harmonics, r, theta, -phi, P) # Ylm^* -> -dx[3]
         # update values
         for l in 0:P
             for m in 0:l

src/compatibility.jl

@@ -109,7 +109,7 @@ function B2M!(system, branch, bodies_index, harmonics, expansion_order)
     return nothing
 end
 
-function direct!(target_system, target_index, source_system, source_index)
+function direct!(target_system, target_index, derivatives_switch, source_system, source_index)
     @warn "direct! function not overloaded for type $(typeof(source_system)); interaction ignored"
     return nothing
 end

src/containers.jl

@@ -49,6 +49,24 @@ struct UniformSourcePanel <: AbstractKernel end
 struct UniformNormalDipolePanel <: AbstractKernel end
 struct UniformSourceNormalDipolePanel <: AbstractKernel end
 
+#####
+##### dispatch convenience functions to determine which derivatives are desired
+#####
+struct DerivativesSwitch{PS,VPS,VS,GS} end
+
+function DerivativesSwitch(scalar_potential::Bool, vector_potential::Bool, velocity::Bool, velocity_gradient::Bool, target_systems::Tuple)
+    return Tuple(DerivativesSwitch{scalar_potential, vector_potential, velocity, velocity_gradient}() for _ in target_systems)
+end
+
+function DerivativesSwitch(scalar_potential::Bool, vector_potential::Bool, velocity::Bool, velocity_gradient::Bool, target_system)
+    return DerivativesSwitch{scalar_potential, vector_potential, velocity, velocity_gradient}()
+end
+
+function DerivativesSwitch(scalar_potential, vector_potential, velocity, velocity_gradient)
+    return Tuple(DerivativesSwitch(ps,vps,vs,gs) for (ps,vps,vs,gs) in zip(scalar_potential, vector_potential, velocity, velocity_gradient))
+end
+
+DerivativesSwitch() = DerivativesSwitch{true, true, true, true}()
 
 #####
 ##### cost parameters

src/derivatives.jl

@@ -1,21 +1,21 @@
 """
 dr_k/dx_idx_j
 """
-function d2rdx2(r, multipole_acceptance_criterion, phi)
+function d2rdx2(r, theta, phi)
     derivatives = zeros(3,3,3)
     derivatives[:,:,1] .= [
-        (1-cos(phi)^2 * sin(multipole_acceptance_criterion)^2)/r -sin(multipole_acceptance_criterion)^2*cos(phi)*sin(phi)/r -sin(multipole_acceptance_criterion)*cos(phi)*cos(multipole_acceptance_criterion)/r;
-        (-sin(multipole_acceptance_criterion)^2*cos(phi)*sin(phi))/r (1-sin(multipole_acceptance_criterion)^2*sin(phi)^2)/r -sin(multipole_acceptance_criterion)*sin(phi)*cos(multipole_acceptance_criterion)/r;
-        -sin(multipole_acceptance_criterion)*cos(phi)*cos(multipole_acceptance_criterion)/r -sin(multipole_acceptance_criterion)*sin(phi)*cos(multipole_acceptance_criterion)/r sin(multipole_acceptance_criterion)^2/r
+        (1-cos(phi)^2 * sin(theta)^2)/r -sin(theta)^2*cos(phi)*sin(phi)/r -sin(theta)*cos(phi)*cos(theta)/r;
+        (-sin(theta)^2*cos(phi)*sin(phi))/r (1-sin(theta)^2*sin(phi)^2)/r -sin(theta)*sin(phi)*cos(theta)/r;
+        -sin(theta)*cos(phi)*cos(theta)/r -sin(theta)*sin(phi)*cos(theta)/r sin(theta)^2/r
     ]
     derivatives[:,:,2] .= [
-        cos(multipole_acceptance_criterion)/sin(multipole_acceptance_criterion)*(1-cos(phi)^2*(1+2*sin(multipole_acceptance_criterion)^2))/r^2 -cos(multipole_acceptance_criterion)/sin(multipole_acceptance_criterion)*sin(phi)*cos(phi)*(1+2*sin(multipole_acceptance_criterion)^2)/r^2 cos(phi)*(1-2*cos(multipole_acceptance_criterion)^2)/r^2;
-        -cos(multipole_acceptance_criterion)/sin(multipole_acceptance_criterion)*sin(phi)*cos(phi)*(1+2*sin(multipole_acceptance_criterion)^2)/r^2 cos(multipole_acceptance_criterion)/sin(multipole_acceptance_criterion)*(1-sin(phi)^2*(1+2*sin(multipole_acceptance_criterion)^2))/r^2 (2*sin(multipole_acceptance_criterion)^2-1)/r^2*sin(phi);
-        cos(phi)*(1-2*cos(multipole_acceptance_criterion)^2)/r^2 (2*sin(multipole_acceptance_criterion)^2-1)/r^2*sin(phi) 2*sin(multipole_acceptance_criterion)*cos(multipole_acceptance_criterion)/r^2
+        cos(theta)/sin(theta)*(1-cos(phi)^2*(1+2*sin(theta)^2))/r^2 -cos(theta)/sin(theta)*sin(phi)*cos(phi)*(1+2*sin(theta)^2)/r^2 cos(phi)*(1-2*cos(theta)^2)/r^2;
+        -cos(theta)/sin(theta)*sin(phi)*cos(phi)*(1+2*sin(theta)^2)/r^2 cos(theta)/sin(theta)*(1-sin(phi)^2*(1+2*sin(theta)^2))/r^2 (2*sin(theta)^2-1)/r^2*sin(phi);
+        cos(phi)*(1-2*cos(theta)^2)/r^2 (2*sin(theta)^2-1)/r^2*sin(phi) 2*sin(theta)*cos(theta)/r^2
     ]
     derivatives[:,:,3] .= [
-        2*cos(phi)*sin(phi)/r^2/sin(multipole_acceptance_criterion)^2 (2*sin(phi)^2-1)/r^2/sin(multipole_acceptance_criterion)^2 0;
-        (2*sin(phi)^2-1)/r^2/sin(multipole_acceptance_criterion)^2 -2*sin(phi)*cos(phi)/r^2/sin(multipole_acceptance_criterion)^2 0;
+        2*cos(phi)*sin(phi)/r^2/sin(theta)^2 (2*sin(phi)^2-1)/r^2/sin(theta)^2 0;
+        (2*sin(phi)^2-1)/r^2/sin(theta)^2 -2*sin(phi)*cos(phi)/r^2/sin(theta)^2 0;
         0 0 0
     ]
     return derivatives
@@ -24,16 +24,16 @@ end
 """
 dr_j/dx_i
 """
-function drdx(r,multipole_acceptance_criterion,phi)
+function drdx(r,theta,phi)
     derivatives = [
-        sin(multipole_acceptance_criterion)*cos(phi) cos(multipole_acceptance_criterion)*cos(phi)/r -sin(phi)/r/sin(multipole_acceptance_criterion);
-        sin(multipole_acceptance_criterion)*sin(phi) cos(multipole_acceptance_criterion)*sin(phi)/r cos(phi)/r/sin(multipole_acceptance_criterion);
-        cos(multipole_acceptance_criterion) -sin(multipole_acceptance_criterion)/r 0
+        sin(theta)*cos(phi) cos(theta)*cos(phi)/r -sin(phi)/r/sin(theta);
+        sin(theta)*sin(phi) cos(theta)*sin(phi)/r cos(phi)/r/sin(theta);
+        cos(theta) -sin(theta)/r 0
     ]
     # derivatives = [
-    #     sin(multipole_acceptance_criterion)*cos(phi) sin(multipole_acceptance_criterion)*sin(phi) cos(multipole_acceptance_criterion);
-    #     cos(multipole_acceptance_criterion)*cos(phi)/r cos(multipole_acceptance_criterion)*sin(phi)/r -sin(multipole_acceptance_criterion)/r;
-    #     -sin(phi)/r/sin(multipole_acceptance_criterion) cos(phi)/r/sin(multipole_acceptance_criterion) 0
+    #     sin(theta)*cos(phi) sin(theta)*sin(phi) cos(theta);
+    #     cos(theta)*cos(phi)/r cos(theta)*sin(phi)/r -sin(theta)/r;
+    #     -sin(phi)/r/sin(theta) cos(phi)/r/sin(theta) 0
     # ]
     return derivatives
 end

src/direct.jl

@@ -22,58 +22,58 @@
 
 Direct calculation of induced potential (no FMM acceleration).
 """
-function direct!(systems::Tuple)
+function direct!(systems::Tuple; derivatives_switches=DerivativesSwitch(true,true,true,true,systems))
     for source_system in systems
-        for target_system in systems
-            _direct!(target_system, 1:get_n_bodies(target_system), source_system, 1:get_n_bodies(source_system))
+        for (target_system, derivatives_switch) in zip(systems, derivatives_switches)
+            _direct!(target_system, 1:get_n_bodies(target_system), derivatives_switch, source_system, 1:get_n_bodies(source_system))
         end
     end
 end
 
-function direct!(system)
-    direct!(system, system)
+function direct!(system; derivatives_switch=DerivativesSwitch{true,true,true,true}())
+    direct!(system, system; derivatives_switch)
 end
 
-@inline function direct!(target_system, source_system)
-    _direct!(target_system, 1:get_n_bodies(target_system), source_system, 1:get_n_bodies(source_system))
+@inline function direct!(target_system, source_system; derivatives_switch=DerivativesSwitch{true,true,true,true}())
+    _direct!(target_system, 1:get_n_bodies(target_system), derivatives_switch, source_system, 1:get_n_bodies(source_system))
 end
 
-function direct!(target_systems::Tuple, source_systems::Tuple)
+function direct!(target_systems::Tuple, source_systems::Tuple; derivatives_switches=Tuple(DerivativesSwitch{true,true,true,true}() for _ in target_systems))
     for source_system in source_systems
-        for target_system in target_systems
-            direct!(target_system, source_system)
+        for (target_system, derivatives_switch) in zip(target_systems, derivatives_switches)
+            direct!(target_system, source_system; derivatives_switch)
         end
     end
 end
 
-function direct!(target_systems::Tuple, source_system)
-    for target_system in target_systems
-        direct!(target_system, source_system)
+function direct!(target_systems::Tuple, source_system; derivatives_switches=Tuple(DerivativesSwitch{true,true,true,true}() for _ in target_systems))
+    for (target_system, derivatives_switch) in zip(target_systems, derivatives_switches)
+        direct!(target_system, source_system; derivatives_switch)
     end
 end
 
-function direct!(target_systems, source_system::Tuple)
+function direct!(target_systems, source_system::Tuple; derivatives_switch=DerivativesSwitch{true,true,true,true}())
     for source_system in source_systems
-        direct!(target_system, source_system)
+        direct!(target_system, source_system; derivatives_switch)
     end
 end
 
 #####
 ##### private methods for dispatch
 #####
-@inline function _direct!(target_system, target_bodies_index, source_system, source_bodies_index)
-    direct!(target_system, target_bodies_index, source_system, source_bodies_index)
+@inline function _direct!(target_system, target_bodies_index, derivatives_switch, source_system, source_bodies_index)
+    direct!(target_system, target_bodies_index, derivatives_switch, source_system, source_bodies_index)
 end
 
-@inline function _direct!(target_system::SortWrapper, target_bodies_index, source_system, source_bodies_index)
-    direct!(target_system.system, target_system.index[target_bodies_index], source_system, source_bodies_index)
+@inline function _direct!(target_system::SortWrapper, target_bodies_index, derivatives_switch, source_system, source_bodies_index)
+    direct!(target_system.system, target_system.index[target_bodies_index], derivatives_switch, source_system, source_bodies_index)
 end
 
-@inline function _direct!(target_system, target_bodies_index, source_system::SortWrapper, source_bodies_index)
-    direct!(target_system, target_bodies_index, source_system.system, source_system.index[source_bodies_index])
+@inline function _direct!(target_system, target_bodies_index, derivatives_switch, source_system::SortWrapper, source_bodies_index)
+    direct!(target_system, target_bodies_index, source_system.system, derivatives_switch, source_system.index[source_bodies_index])
 end
 
-@inline function _direct!(target_system::SortWrapper, target_bodies_index, source_system::SortWrapper, source_bodies_index)
-    direct!(target_system.system, target_system.index[target_bodies_index], source_system.system, source_system.index[source_bodies_index])
+@inline function _direct!(target_system::SortWrapper, target_bodies_index, derivatives_switch, source_system::SortWrapper, source_bodies_index)
+    direct!(target_system.system, target_system.index[target_bodies_index], derivatives_switch, source_system.system, source_system.index[source_bodies_index])
 end
 

src/estimate_cost.jl

@@ -41,10 +41,10 @@ function allocate_l2b(expansion_order, type=Float64)
     potential_jacobian = zeros(type,3,4)
     potential_hessian = zeros(type,3,3,4)
     derivative_harmonics = zeros(Complex{type}, ((expansion_order+1) * (expansion_order+2)) >> 1)
-    derivative_harmonics_multipole_acceptance_criterion = zeros(Complex{type}, ((expansion_order+1) * (expansion_order+2)) >> 1)
-    derivative_harmonics_multipole_acceptance_criterion_2 = zeros(Complex{type}, ((expansion_order+1) * (expansion_order+2)) >> 1)
+    derivative_harmonics_theta = zeros(Complex{type}, ((expansion_order+1) * (expansion_order+2)) >> 1)
+    derivative_harmonics_theta_2 = zeros(Complex{type}, ((expansion_order+1) * (expansion_order+2)) >> 1)
     workspace = zeros(type,3,4)
-    return vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_multipole_acceptance_criterion, derivative_harmonics_multipole_acceptance_criterion_2, workspace
+    return vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_theta, derivative_harmonics_theta_2, workspace
 end
 
 function get_tau_b2m(system, branch, b2m_n_bodies, harmonics_m2m_l2l, expansion_order)
@@ -74,15 +74,15 @@ function estimate_tau_fmm(expansion_order, type=Float64)
     # harmonics_l2l, L = allocate_m2m_l2l(expansion_order, type)
     # alloc_l2l = @belapsed harmonics_l2l, L = allocate_m2m_l2l(expansion_order, type)
     # tau_l2l = @belapsed L2L!(branch1, branch2, harmonics_l2l, L, expansion_order)
-    vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_multipole_acceptance_criterion, derivative_harmonics_multipole_acceptance_criterion_2, workspace = allocate_l2b(expansion_order, type)
-    alloc_l2b = @elapsed vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_multipole_acceptance_criterion, derivative_harmonics_multipole_acceptance_criterion_2, workspace = allocate_l2b(expansion_order, type)
-    alloc_l2b = @elapsed vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_multipole_acceptance_criterion, derivative_harmonics_multipole_acceptance_criterion_2, workspace = allocate_l2b(expansion_order, type)
+    vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_theta, derivative_harmonics_theta_2, workspace = allocate_l2b(expansion_order, type)
+    alloc_l2b = @elapsed vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_theta, derivative_harmonics_theta_2, workspace = allocate_l2b(expansion_order, type)
+    alloc_l2b = @elapsed vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_theta, derivative_harmonics_theta_2, workspace = allocate_l2b(expansion_order, type)
     body_position = rand(SVector{3,type})
     expansion_center = 10 * rand(SVector{3,type})
     tau_l2b = @elapsed L2B_loop!(vector_potential, potential_jacobian, potential_hessian, body_position, expansion_center, 
-                branch2.local_expansion, derivative_harmonics, derivative_harmonics_multipole_acceptance_criterion, derivative_harmonics_multipole_acceptance_criterion_2, expansion_order, workspace)
+                branch2.local_expansion, derivative_harmonics, derivative_harmonics_theta, derivative_harmonics_theta_2, expansion_order, workspace)
     tau_l2b = @elapsed L2B_loop!(vector_potential, potential_jacobian, potential_hessian, body_position, expansion_center, 
-                branch2.local_expansion, derivative_harmonics, derivative_harmonics_multipole_acceptance_criterion, derivative_harmonics_multipole_acceptance_criterion_2, expansion_order, workspace)
+                branch2.local_expansion, derivative_harmonics, derivative_harmonics_theta, derivative_harmonics_theta_2, expansion_order, workspace)
     return alloc_m2m_l2l, tau_m2m_l2l, alloc_m2l, tau_m2l, alloc_l2b, tau_l2b
 end