diff --git a/src/FastMultipole.jl b/src/FastMultipole.jl index 73b15f6..f2491db 100644 --- a/src/FastMultipole.jl +++ b/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 diff --git a/src/b2m.jl b/src/b2m.jl index 939c316..09d7dc4 100644 --- a/src/b2m.jl +++ b/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 diff --git a/src/compatibility.jl b/src/compatibility.jl index 684a743..e3d4e74 100644 --- a/src/compatibility.jl +++ b/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 diff --git a/src/containers.jl b/src/containers.jl index b637b59..31396ae 100644 --- a/src/containers.jl +++ b/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 diff --git a/src/derivatives.jl b/src/derivatives.jl index b16f00e..301f864 100644 --- a/src/derivatives.jl +++ b/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 \ No newline at end of file diff --git a/src/direct.jl b/src/direct.jl index 439066f..2203707 100644 --- a/src/direct.jl +++ b/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 diff --git a/src/estimate_cost.jl b/src/estimate_cost.jl index fc26eab..645ab53 100644 --- a/src/estimate_cost.jl +++ b/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 diff --git a/src/fmm.jl b/src/fmm.jl index c290c0b..fbf2f6c 100644 --- a/src/fmm.jl +++ b/src/fmm.jl @@ -1,38 +1,38 @@ ##### ##### direct interactions ##### -function P2P!(target_system, target_branch::SingleBranch, source_system, source_branch::SingleBranch) - _direct!(target_system, target_branch.bodies_index, source_system, source_branch.bodies_index) +function P2P!(target_system, target_branch::SingleBranch, derivatives_switch, source_system, source_branch::SingleBranch) + _direct!(target_system, target_branch.bodies_index, derivatives_switch, source_system, source_branch.bodies_index) end -function P2P!(target_system, target_branch::SingleBranch, source_systems, source_branch::MultiBranch) +function P2P!(target_system, target_branch::SingleBranch, derivatives_switch, source_systems, source_branch::MultiBranch) for (source_system, source_bodies_index) in zip(source_systems, source_branch.bodies_index) - _direct!(target_system, target_branch.bodies_index, source_system, source_bodies_index) + _direct!(target_system, target_branch.bodies_index, derivatives_switch, source_system, source_bodies_index) end end -function P2P!(target_systems, target_branch::MultiBranch, source_system, source_branch::SingleBranch) - for (target_system, target_bodies_index) in zip(target_systems, target_branch.bodies_index) - _direct!(target_system, target_bodies_index, source_system, source_branch.bodies_index) +function P2P!(target_systems, target_branch::MultiBranch, derivatives_switches, source_system, source_branch::SingleBranch) + for (target_system, target_bodies_index, derivatives_switch) in zip(target_systems, target_branch.bodies_index, derivatives_switches) + _direct!(target_system, target_bodies_index, derivatives_switch, source_system, source_branch.bodies_index) end end -function P2P!(target_systems, target_branch::MultiBranch, source_systems, source_branch::MultiBranch) +function P2P!(target_systems, target_branch::MultiBranch, derivatives_switches, source_systems, source_branch::MultiBranch) for (source_system, source_bodies_index) in zip(source_systems, source_branch.bodies_index) - for (target_system, target_bodies_index) in zip(target_systems, target_branch.bodies_index) - _direct!(target_system, target_bodies_index, source_system, source_bodies_index) + for (target_system, target_bodies_index, derivatives_switch) in zip(target_systems, target_branch.bodies_index, derivatives_switches) + _direct!(target_system, target_bodies_index, derivatives_switch, source_system, source_bodies_index) end end end -function P2P!(target_systems, target_branch::MultiBranch, source_system, source_bodies_index::UnitRange) - for (target_system, target_bodies_index) in zip(target_systems, target_branch.bodies_index) - _direct!(target_system, target_bodies_index, source_system, source_bodies_index) +function P2P!(target_systems, target_branch::MultiBranch, derivatives_switches, source_system, source_bodies_index::UnitRange) + for (target_system, target_bodies_index, derivatives_switch) in zip(target_systems, target_branch.bodies_index, derivatives_switches) + _direct!(target_system, target_bodies_index, derivatives_switch, source_system, source_bodies_index) end end -function P2P!(target_system, target_branch::SingleBranch, source_system, source_bodies_index::UnitRange) - _direct!(target_system, target_branch.bodies_index, source_system, source_bodies_index) +function P2P!(target_system, target_branch::SingleBranch, derivatives_switch, source_system, source_bodies_index::UnitRange) + _direct!(target_system, target_branch.bodies_index, derivatives_switch, source_system, source_bodies_index) end ##### @@ -187,13 +187,13 @@ end ##### ##### horizontal pass ##### -function nearfield_singlethread!(target_system, target_branches, source_system, source_branches, direct_list) +function nearfield_singlethread!(target_system, target_branches, derivatives_switch, source_system, source_branches, direct_list) for (i_target, j_source) in direct_list - P2P!(target_system, target_branches[i_target], source_system, source_branches[j_source]) + P2P!(target_system, target_branches[i_target], derivatives_switch, source_system, source_branches[j_source]) end end -function nearfield_multithread!(target_system, target_branches, source_systems::Tuple, source_branches, direct_list) +function nearfield_multithread!(target_system, target_branches, derivatives_switch, source_systems::Tuple, source_branches, direct_list) ## load balance n_threads = Threads.nthreads() assignments = Vector{UnitRange{Int64}}(undef,n_threads) @@ -241,7 +241,7 @@ function nearfield_multithread!(target_system, target_branches, source_systems:: target_branch = target_branches[i_target] source_bodies_index = source_branches[j_source].bodies_index[i_source_system] Threads.lock(target_branch.lock) do - P2P!(target_system, target_branch, source_system, source_bodies_index) + P2P!(target_system, target_branch, derivatives_switch, source_system, source_bodies_index) end end end @@ -250,7 +250,7 @@ function nearfield_multithread!(target_system, target_branches, source_systems:: return nothing end -function nearfield_multithread!(target_system, target_branches, source_system, source_branches, direct_list) +function nearfield_multithread!(target_system, target_branches, derivatives_switch, source_system, source_branches, direct_list) ## load balance n_threads = Threads.nthreads() assignments = Vector{UnitRange{Int64}}(undef,n_threads) @@ -295,7 +295,7 @@ function nearfield_multithread!(target_system, target_branches, source_system, s for (i_target, j_source) in view(direct_list, assignment) target_branch = target_branches[i_target] Threads.lock(target_branch.lock) do - P2P!(target_system, target_branch, source_system, source_branches[j_source]) + P2P!(target_system, target_branch, derivatives_switch, source_system, source_branches[j_source]) end end end @@ -371,19 +371,20 @@ end potential_jacobian = zeros(float_type,3,4) potential_hessian = zeros(float_type,3,3,4) derivative_harmonics = zeros(expansion_type, 2, ((P+1) * (P+2)) >> 1) - derivative_harmonics_multipole_acceptance_criterion = zeros(expansion_type, 2, ((P+1) * (P+2)) >> 1) - derivative_harmonics_multipole_acceptance_criterion_2 = zeros(expansion_type, 2, ((P+1) * (P+2)) >> 1) + derivative_harmonics_theta = zeros(expansion_type, 2, ((P+1) * (P+2)) >> 1) + derivative_harmonics_theta_2 = zeros(expansion_type, 2, ((P+1) * (P+2)) >> 1) workspace = zeros(float_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 downward_pass_singlethread!(branches, systems, expansion_order::Val{P}) where P +function downward_pass_singlethread!(branches, systems, derivatives_switches, expansion_order::Val{P}) where P regular_harmonics = zeros(eltype(branches[1].multipole_expansion), 2, (P+1)*(P+1)) # L = zeros(eltype(branches[1].multipole_expansion),4) - vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_multipole_acceptance_criterion, derivative_harmonics_multipole_acceptance_criterion_2, workspace = preallocate_l2b(eltype(branches[1]), eltype(branches[1].multipole_expansion), expansion_order) + vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_theta, derivative_harmonics_theta_2, workspace = preallocate_l2b(eltype(branches[1]), eltype(branches[1].multipole_expansion), expansion_order) for branch in branches if branch.n_branches == 0 # leaf level - L2B!(systems, branch, expansion_order, vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_multipole_acceptance_criterion, derivative_harmonics_multipole_acceptance_criterion_2, workspace) + L2B!(systems, branch, derivatives_switches, expansion_order, vector_potential, potential_jacobian, potential_hessian, + derivative_harmonics, derivative_harmonics_theta, derivative_harmonics_theta_2, workspace) else for child_branch in view(branches,branch.branch_index) L2L!(branch, child_branch, regular_harmonics, branch.ML, expansion_order) @@ -419,7 +420,7 @@ function translate_locals_multithread!(branches, expansion_order::Val{P}, levels return nothing end -function local_2_body_multithread!(branches, systems, expansion_order, leaf_index) +function local_2_body_multithread!(branches, systems, derivatives_switches, expansion_order, leaf_index) # create assignments n_threads = Threads.nthreads() @@ -456,20 +457,20 @@ function local_2_body_multithread!(branches, systems, expansion_order, leaf_inde # spread remainder across rem chunks Threads.@threads for i_thread in eachindex(assignments) assignment = assignments[i_thread] - vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_multipole_acceptance_criterion, derivative_harmonics_multipole_acceptance_criterion_2, workspace = containers[i_thread] + vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_theta, derivative_harmonics_theta_2, workspace = containers[i_thread] for i_leaf in view(leaf_index,assignment) leaf = branches[i_leaf] - L2B!(systems, leaf, expansion_order, vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_multipole_acceptance_criterion, derivative_harmonics_multipole_acceptance_criterion_2, workspace) + L2B!(systems, leaf, derivatives_switches, expansion_order, vector_potential, potential_jacobian, potential_hessian, derivative_harmonics, derivative_harmonics_theta, derivative_harmonics_theta_2, workspace) end end end -function downward_pass_multithread!(branches, systems, expansion_order, levels_index, leaf_index) +function downward_pass_multithread!(branches, systems, derivatives_switch, expansion_order, levels_index, leaf_index) # m2m translation translate_locals_multithread!(branches, expansion_order, levels_index) # local to body interaction - local_2_body_multithread!(branches, systems, expansion_order, leaf_index) + local_2_body_multithread!(branches, systems, derivatives_switch, expansion_order, leaf_index) end ##### @@ -508,11 +509,27 @@ end ##### ##### running FMM ##### -function fmm!(tree::Tree, systems; multipole_acceptance_criterion=0.4, reset_tree=true, nearfield=true, farfield=true, self_induced=true, unsort_bodies=true) - fmm!(tree, systems, tree, systems; multipole_acceptance_criterion, reset_source_tree=reset_tree, reset_target_tree=false, nearfield=nearfield, farfield=farfield, self_induced, unsort_source_bodies=unsort_bodies, unsort_target_bodies=false) +function fmm!(tree::Tree, systems; + scalar_potential=true, vector_potential=true, velocity=true, velocity_gradient=true, + multipole_acceptance_criterion=0.4, reset_tree=true, + nearfield=true, farfield=true, self_induced=true, + unsort_bodies=true +) + fmm!(tree, systems, tree, systems; + scalar_potential, vector_potential, velocity, velocity_gradient, + multipole_acceptance_criterion, reset_source_tree=reset_tree, reset_target_tree=false, + nearfield, farfield, self_induced, + unsort_source_bodies=unsort_bodies, unsort_target_bodies=false + ) end -function fmm!(target_tree::Tree, target_systems, source_tree::Tree, source_systems; multipole_acceptance_criterion=0.4, reset_source_tree=true, reset_target_tree=true, nearfield=true, farfield=true, self_induced=true, unsort_source_bodies=true, unsort_target_bodies=true) +function fmm!(target_tree::Tree, target_systems, source_tree::Tree, source_systems; + scalar_potential=true, vector_potential=true, velocity=true, velocity_gradient=true, + multipole_acceptance_criterion=0.4, + reset_source_tree=true, reset_target_tree=true, + nearfield=true, farfield=true, self_induced=true, + unsort_source_bodies=true, unsort_target_bodies=true +) # reset multipole/local expansions reset_target_tree && (reset_expansions!(source_tree)) reset_source_tree && (reset_expansions!(source_tree)) @@ -520,17 +537,20 @@ function fmm!(target_tree::Tree, target_systems, source_tree::Tree, source_syste # create interaction lists m2l_list, direct_list = build_interaction_lists(target_tree.branches, source_tree.branches, multipole_acceptance_criterion, farfield, nearfield, self_induced) + # assemble derivatives switch + derivatives_switch = DerivativesSwitch(scalar_potential, vector_potential, velocity, velocity_gradient, target_systems) + # run FMM if Threads.nthreads() == 1 - nearfield && (nearfield_singlethread!(target_systems, target_tree.branches, source_systems, source_tree.branches, direct_list)) + nearfield && (nearfield_singlethread!(target_systems, target_tree.branches, derivatives_switch, source_systems, source_tree.branches, direct_list)) if farfield upward_pass_singlethread!(source_tree.branches, source_systems, source_tree.expansion_order) horizontal_pass_singlethread!(target_tree.branches, source_tree.branches, m2l_list, source_tree.expansion_order) - downward_pass_singlethread!(target_tree.branches, target_systems, target_tree.expansion_order) + downward_pass_singlethread!(target_tree.branches, target_systems, derivatives_switch, target_tree.expansion_order) end else # multithread # println("nearfield") - nearfield && length(direct_list) > 0 && (nearfield_multithread!(target_systems, target_tree.branches, source_systems, source_tree.branches, direct_list)) + nearfield && length(direct_list) > 0 && (nearfield_multithread!(target_systems, target_tree.branches, derivatives_switch, source_systems, source_tree.branches, direct_list)) # @time nearfield && length(direct_list) > 0 && (nearfield_multithread!(target_systems, target_tree.branches, source_systems, source_tree.branches, direct_list)) if farfield # println("upward pass") @@ -540,7 +560,7 @@ function fmm!(target_tree::Tree, target_systems, source_tree::Tree, source_syste length(m2l_list) > 0 && (horizontal_pass_multithread!(target_tree.branches, source_tree.branches, m2l_list, source_tree.expansion_order)) # @time length(m2l_list) > 0 && (horizontal_pass_multithread!(target_tree.branches, source_tree.branches, m2l_list, source_tree.expansion_order)) # println("downward pass") - downward_pass_multithread!(target_tree.branches, target_systems, target_tree.expansion_order, target_tree.levels_index, target_tree.leaf_index) + downward_pass_multithread!(target_tree.branches, target_systems, derivatives_switch, target_tree.expansion_order, target_tree.levels_index, target_tree.leaf_index) # @time downward_pass_multithread!(target_tree.branches, target_systems, target_tree.expansion_order, target_tree.levels_index, target_tree.leaf_index) end end @@ -550,12 +570,23 @@ function fmm!(target_tree::Tree, target_systems, source_tree::Tree, source_syste unsort_source_bodies && (unsort!(source_systems, source_tree)) end -function fmm!(systems; expansion_order=5, n_per_branch=50, multipole_acceptance_criterion=0.4, nearfield=true, farfield=true, self_induced=true, unsort_bodies=true, shrink_recenter=true, save_tree=false, save_name="tree") +function fmm!(systems; + scalar_potential=true, vector_potential=true, velocity=true, velocity_gradient=true, + expansion_order=5, n_per_branch=50, multipole_acceptance_criterion=0.4, + nearfield=true, farfield=true, self_induced=true, + unsort_bodies=true, shrink_recenter=true, + save_tree=false, save_name="tree" +) # create tree tree = Tree(systems; expansion_order, n_per_branch, shrink_recenter) # perform fmm - fmm!(tree, systems; multipole_acceptance_criterion, reset_tree=false, nearfield, farfield, self_induced, unsort_bodies) + fmm!(tree, systems; + scalar_potential, vector_potential, velocity, velocity_gradient, + multipole_acceptance_criterion, reset_tree=false, + nearfield, farfield, self_induced, + unsort_bodies + ) # visualize save_tree && (visualize(save_name, systems, tree)) @@ -596,16 +627,29 @@ Apply all interactions of `source_systems` acting on `target_systems` using the - `multipole_acceptance_criterion::Float64`: number between 0 and 1 (often denoted theta in [0,1]) controls the accuracy by determining the non-dimensional distance after which multipoles are used; 0 means an infinite distance (no error, high cost), and 1 means barely convergent (high error, low cost) """ -function fmm!(target_systems, source_systems; expansion_order=5, n_per_branch_source=50, n_per_branch_target=50, multipole_acceptance_criterion=0.4, nearfield=true, farfield=true, self_induced=true, unsort_source_bodies=true, unsort_target_bodies=true, source_shrink_recenter=true, target_shrink_recenter=true, save_tree=false, save_name="tree") +function fmm!(target_systems, source_systems; + scalar_potential=true, vector_potential=true, velocity=true, velocity_gradient=true, + expansion_order=5, n_per_branch_source=50, n_per_branch_target=50, multipole_acceptance_criterion=0.4, + nearfield=true, farfield=true, self_induced=true, + unsort_source_bodies=true, unsort_target_bodies=true, + source_shrink_recenter=true, target_shrink_recenter=true, + save_tree=false, save_name="tree" +) # check for duplicate systems target_systems = wrap_duplicates(target_systems, source_systems) # create trees - source_tree = Tree(source_systems; expansion_order=expansion_order, n_per_branch=n_per_branch_source, shrink_recenter=source_shrink_recenter) - target_tree = Tree(target_systems; expansion_order=expansion_order, n_per_branch=n_per_branch_target, shrink_recenter=target_shrink_recenter) + source_tree = Tree(source_systems; expansion_order, n_per_branch=n_per_branch_source, shrink_recenter=source_shrink_recenter) + target_tree = Tree(target_systems; expansion_order, n_per_branch=n_per_branch_target, shrink_recenter=target_shrink_recenter) # perform fmm - fmm!(target_tree, target_systems, source_tree, source_systems; multipole_acceptance_criterion=multipole_acceptance_criterion, reset_source_tree=false, reset_target_tree=false, nearfield=nearfield, farfield=farfield, self_induced, unsort_source_bodies=unsort_source_bodies, unsort_target_bodies=unsort_target_bodies) + fmm!(target_tree, target_systems, source_tree, source_systems; + scalar_potential, vector_potential, velocity, velocity_gradient, + multipole_acceptance_criterion, + reset_source_tree=false, reset_target_tree=false, + nearfield, farfield, self_induced, + unsort_source_bodies, unsort_target_bodies + ) # visualize save_tree && (visualize(save_name, systems, tree)) diff --git a/src/spherical.jl b/src/spherical.jl index dc476bd..6d13521 100644 --- a/src/spherical.jl +++ b/src/spherical.jl @@ -1,100 +1,111 @@ function cartesian_2_spherical!(dx; EPSILON=1e-10) r = sqrt(dx' * dx) - multipole_acceptance_criterion = r < EPSILON ? 0.0 : acos(dx[3] / r) + theta = r < EPSILON ? 0.0 : acos(dx[3] / r) phi = r < EPSILON ? 0.0 : atan(dx[2], dx[1]) - dx .= r, multipole_acceptance_criterion, phi + dx .= r, theta, phi end function cartesian_2_spherical(dx; EPSILON=1e-10) r = sqrt(dx' * dx) - multipole_acceptance_criterion = r < EPSILON ? 0.0 : acos(dx[3] / r) + theta = r < EPSILON ? 0.0 : acos(dx[3] / r) phi = r < EPSILON ? 0.0 : atan(dx[2], dx[1]) - dx = SVector{3}(r, multipole_acceptance_criterion, phi) + dx = SVector{3}(r, theta, phi) end function cartesian_2_spherical(x, y, z; EPSILON=1e-10) r = sqrt(x*x + y*y + z*z) - multipole_acceptance_criterion = r < EPSILON ? 0.0 : acos(z / r) + theta = r < EPSILON ? 0.0 : acos(z / r) phi = r < EPSILON ? 0.0 : atan(y, x) - return r, multipole_acceptance_criterion, phi + return r, theta, phi end -function spherical_2_cartesian!(potential_jacobian, potential_hessian, workspace, rho, multipole_acceptance_criterion, phi) +spherical_2_cartesian!(potential_jacobian, potential_hessian, workspace, rho, theta, phi, derivatives_switch::DerivativesSwitch{PS,false,false}) where PS = nothing + +function spherical_2_cartesian!(potential_jacobian, potential_hessian, workspace, rho, theta, phi, derivatives_switch::DerivativesSwitch{PS,VPS,VS,GS}) where {PS,VPS,VS,GS} a = 0 # get partial derivatives of the coordinates - s_multipole_acceptance_criterion, c_multipole_acceptance_criterion = sincos(multipole_acceptance_criterion) + s_theta, c_theta = sincos(theta) s_phi, c_phi = sincos(phi) drjdxi = @SMatrix [ - s_multipole_acceptance_criterion*c_phi c_multipole_acceptance_criterion*c_phi/rho -s_phi/rho/s_multipole_acceptance_criterion - s_multipole_acceptance_criterion * s_phi c_multipole_acceptance_criterion * s_phi / rho c_phi / rho / s_multipole_acceptance_criterion - c_multipole_acceptance_criterion -s_multipole_acceptance_criterion / rho 0.0 + s_theta*c_phi c_theta*c_phi/rho -s_phi/rho/s_theta + s_theta * s_phi c_theta * s_phi / rho c_phi / rho / s_theta + c_theta -s_theta / rho 0.0 ] - # convert Hessian to cartesian coordinates - workspace3x3 = view(workspace,:,1:3) - for ind in 1:4 - # potential_hessian[:,:,ind] .= drjdxi * potential_hessian[:,:,ind] * transpose(drjdxi) - mul!(workspace3x3, drjdxi, view(potential_hessian,:,:,ind)) - mul!(view(potential_hessian,:,:,ind), workspace3x3, transpose(drjdxi)) - end - for k_coord in 1:3 # loop over r, multipole_acceptance_criterion, and phi to save me some space on the stack - if k_coord == 1 # r coordinate - drkdxidxj = @SMatrix [ - (1-c_phi^2 * s_multipole_acceptance_criterion^2)/rho -s_multipole_acceptance_criterion^2*c_phi*s_phi/rho -s_multipole_acceptance_criterion*c_phi*c_multipole_acceptance_criterion/rho; - (-s_multipole_acceptance_criterion^2*c_phi*s_phi)/rho (1-s_multipole_acceptance_criterion^2*s_phi^2)/rho -s_multipole_acceptance_criterion*s_phi*c_multipole_acceptance_criterion/rho; - -s_multipole_acceptance_criterion*c_phi*c_multipole_acceptance_criterion/rho -s_multipole_acceptance_criterion*s_phi*c_multipole_acceptance_criterion/rho s_multipole_acceptance_criterion^2/rho - ] - elseif k_coord == 2 # multipole_acceptance_criterion coordinate - drkdxidxj = @SMatrix [ - c_multipole_acceptance_criterion/s_multipole_acceptance_criterion*(1-c_phi^2*(1+2*s_multipole_acceptance_criterion^2))/rho^2 -c_multipole_acceptance_criterion/s_multipole_acceptance_criterion*s_phi*c_phi*(1+2*s_multipole_acceptance_criterion^2)/rho^2 c_phi*(1-2*c_multipole_acceptance_criterion^2)/rho^2; - -c_multipole_acceptance_criterion/s_multipole_acceptance_criterion*s_phi*c_phi*(1+2*s_multipole_acceptance_criterion^2)/rho^2 c_multipole_acceptance_criterion/s_multipole_acceptance_criterion*(1-s_phi^2*(1+2*s_multipole_acceptance_criterion^2))/rho^2 (2*s_multipole_acceptance_criterion^2-1)/rho^2*s_phi; - c_phi*(1-2*c_multipole_acceptance_criterion^2)/rho^2 (2*s_multipole_acceptance_criterion^2-1)/rho^2*s_phi 2*s_multipole_acceptance_criterion*c_multipole_acceptance_criterion/rho^2 - ] - else # phi coordinate - drkdxidxj = @SMatrix [ - 2*c_phi*s_phi/rho^2/s_multipole_acceptance_criterion^2 (2*s_phi^2-1)/rho^2/s_multipole_acceptance_criterion^2 0; - (2*s_phi^2-1)/rho^2/s_multipole_acceptance_criterion^2 -2*s_phi*c_phi/rho^2/s_multipole_acceptance_criterion^2 0; - 0 0 0 - ] - end + + if GS + # convert Hessian to cartesian coordinates + workspace3x3 = view(workspace,:,1:3) for ind in 1:4 - view(potential_hessian,:,:,ind) .+= drkdxidxj * potential_jacobian[k_coord,ind] + # potential_hessian[:,:,ind] .= drjdxi * potential_hessian[:,:,ind] * transpose(drjdxi) + mul!(workspace3x3, drjdxi, view(potential_hessian,:,:,ind)) + mul!(view(potential_hessian,:,:,ind), workspace3x3, transpose(drjdxi)) + end + for k_coord in 1:3 # loop over r, theta, and phi to save me some space on the stack + if k_coord == 1 # r coordinate + drkdxidxj = @SMatrix [ + (1-c_phi^2 * s_theta^2)/rho -s_theta^2*c_phi*s_phi/rho -s_theta*c_phi*c_theta/rho; + (-s_theta^2*c_phi*s_phi)/rho (1-s_theta^2*s_phi^2)/rho -s_theta*s_phi*c_theta/rho; + -s_theta*c_phi*c_theta/rho -s_theta*s_phi*c_theta/rho s_theta^2/rho + ] + elseif k_coord == 2 # theta coordinate + drkdxidxj = @SMatrix [ + c_theta/s_theta*(1-c_phi^2*(1+2*s_theta^2))/rho^2 -c_theta/s_theta*s_phi*c_phi*(1+2*s_theta^2)/rho^2 c_phi*(1-2*c_theta^2)/rho^2; + -c_theta/s_theta*s_phi*c_phi*(1+2*s_theta^2)/rho^2 c_theta/s_theta*(1-s_phi^2*(1+2*s_theta^2))/rho^2 (2*s_theta^2-1)/rho^2*s_phi; + c_phi*(1-2*c_theta^2)/rho^2 (2*s_theta^2-1)/rho^2*s_phi 2*s_theta*c_theta/rho^2 + ] + else # phi coordinate + drkdxidxj = @SMatrix [ + 2*c_phi*s_phi/rho^2/s_theta^2 (2*s_phi^2-1)/rho^2/s_theta^2 0; + (2*s_phi^2-1)/rho^2/s_theta^2 -2*s_phi*c_phi/rho^2/s_theta^2 0; + 0 0 0 + ] + end + for ind in 1:4 + view(potential_hessian,:,:,ind) .+= drkdxidxj * potential_jacobian[k_coord,ind] + end end end - workspace .= potential_jacobian - # for some reason mul! allocates for nonsquare matrix matrix products - mul!(view(potential_jacobian,:,1),drjdxi,view(workspace,1:3,1)) - mul!(view(potential_jacobian,:,2:4),drjdxi,view(workspace,1:3,2:4)) + if VS + workspace .= potential_jacobian + # for some reason mul! allocates for nonsquare matrix matrix products + mul!(view(potential_jacobian,:,1),drjdxi,view(workspace,1:3,1)) + mul!(view(potential_jacobian,:,2:4),drjdxi,view(workspace,1:3,2:4)) + end return nothing end -function flatten_derivatives!(jacobian, hessian) - # velocity - jacobian[1,1] = -jacobian[1,1] + jacobian[2,4] - jacobian[3,3] - jacobian[2,1] = -jacobian[2,1] + jacobian[3,2] - jacobian[1,4] - jacobian[3,1] = -jacobian[3,1] + jacobian[1,3] - jacobian[2,2] +function flatten_derivatives!(jacobian, hessian, derivatives_switch::DerivativesSwitch{PS,VPS,VS,GS}) where {PS,VPS,VS,GS} + if VS + # velocity + jacobian[1,1] = -jacobian[1,1] + jacobian[2,4] - jacobian[3,3] + jacobian[2,1] = -jacobian[2,1] + jacobian[3,2] - jacobian[1,4] + jacobian[3,1] = -jacobian[3,1] + jacobian[1,3] - jacobian[2,2] + end - # velocity gradient - hessian[1,1,1] = -hessian[1,1,1]+hessian[2,1,4]-hessian[3,1,3] - hessian[2,1,1] = -hessian[2,1,1]+hessian[3,1,2]-hessian[1,1,4] - hessian[3,1,1] = -hessian[3,1,1]+hessian[1,1,3]-hessian[2,1,2] - hessian[1,2,1] = -hessian[1,2,1]+hessian[2,2,4]-hessian[3,2,3] - hessian[2,2,1] = -hessian[2,2,1]+hessian[3,2,2]-hessian[1,2,4] - hessian[3,2,1] = -hessian[3,2,1]+hessian[1,2,3]-hessian[2,2,2] - hessian[1,3,1] = -hessian[1,3,1]+hessian[2,3,4]-hessian[3,3,3] - hessian[2,3,1] = -hessian[2,3,1]+hessian[3,3,2]-hessian[1,3,4] - hessian[3,3,1] = -hessian[3,3,1]+hessian[1,3,3]-hessian[2,3,2] + if GS + # velocity gradient + hessian[1,1,1] = -hessian[1,1,1]+hessian[2,1,4]-hessian[3,1,3] + hessian[2,1,1] = -hessian[2,1,1]+hessian[3,1,2]-hessian[1,1,4] + hessian[3,1,1] = -hessian[3,1,1]+hessian[1,1,3]-hessian[2,1,2] + hessian[1,2,1] = -hessian[1,2,1]+hessian[2,2,4]-hessian[3,2,3] + hessian[2,2,1] = -hessian[2,2,1]+hessian[3,2,2]-hessian[1,2,4] + hessian[3,2,1] = -hessian[3,2,1]+hessian[1,2,3]-hessian[2,2,2] + hessian[1,3,1] = -hessian[1,3,1]+hessian[2,3,4]-hessian[3,3,3] + hessian[2,3,1] = -hessian[2,3,1]+hessian[3,3,2]-hessian[1,3,4] + hessian[3,3,1] = -hessian[3,3,1]+hessian[1,3,3]-hessian[2,3,2] + end end @inline odd_or_even(n::Int) = (n & 1) == 1 ? -1 : 1 @inline ipow2l(n::Int) = n >= 0 ? 1 : odd_or_even(n); -function regular_harmonic!(harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, rho, multipole_acceptance_criterion, phi::TF, P) where TF - y,x = sincos(multipole_acceptance_criterion) - invY = y == 0 ? 0 : 1 / y +function regular_harmonic!(harmonics, harmonics_theta, harmonics_theta_2, rho, theta, phi::TF, P) where TF + y,x = sincos(theta) + invY = iszero(y) ? 0 : 1 / y fact = 1.0 pl = 1.0 rhom = 1.0 @@ -113,12 +124,12 @@ function regular_harmonic!(harmonics, harmonics_multipole_acceptance_criterion, # harmonics[lml] = conj(harmonics[lpl]) p1 = p p = x * (2 * m + 1) * p1 - #harmonics_multipole_acceptance_criterion[lpl] = rhom * (p - (m + 1) * x * p1) * invY * eim - harmonics_multipole_acceptance_criterion[1,lpl] = rhom * (p - (m + 1) * x * p1) * invY * eim[1] - harmonics_multipole_acceptance_criterion[2,lpl] = rhom * (p - (m + 1) * x * p1) * invY * eim[2] - #harmonics_multipole_acceptance_criterion_2[lpl] = rhom * (-x * p + (-m + (m+1)^2 * x^2) * p1) * invY^2 * eim - harmonics_multipole_acceptance_criterion_2[1,lpl] = rhom * (-x * p + (-m + (m+1)^2 * x^2) * p1) * invY^2 * eim[1] - harmonics_multipole_acceptance_criterion_2[2,lpl] = rhom * (-x * p + (-m + (m+1)^2 * x^2) * p1) * invY^2 * eim[2] + #harmonics_theta[lpl] = rhom * (p - (m + 1) * x * p1) * invY * eim + harmonics_theta[1,lpl] = rhom * (p - (m + 1) * x * p1) * invY * eim[1] + harmonics_theta[2,lpl] = rhom * (p - (m + 1) * x * p1) * invY * eim[2] + #harmonics_theta_2[lpl] = rhom * (-x * p + (-m + (m+1)^2 * x^2) * p1) * invY^2 * eim + harmonics_theta_2[1,lpl] = rhom * (-x * p + (-m + (m+1)^2 * x^2) * p1) * invY^2 * eim[1] + harmonics_theta_2[2,lpl] = rhom * (-x * p + (-m + (m+1)^2 * x^2) * p1) * invY^2 * eim[2] rhom *= rho rhol = rhom for l=m+1:P @@ -133,12 +144,12 @@ function regular_harmonic!(harmonics, harmonics_multipole_acceptance_criterion, p2 = p1 p1 = p p = (x * (2 * l + 1) * p1 - (l + m) * p2) / (l - m + 1) - #harmonics_multipole_acceptance_criterion[lpm] = rhol * ((l - m + 1) * p - (l + 1) * x * p1) * invY * eim - harmonics_multipole_acceptance_criterion[1,lpm] = rhol * ((l - m + 1) * p - (l + 1) * x * p1) * invY * eim[1] - harmonics_multipole_acceptance_criterion[2,lpm] = rhol * ((l - m + 1) * p - (l + 1) * x * p1) * invY * eim[2] - #harmonics_multipole_acceptance_criterion_2[lpm] = rhol * ((m-l-1) * x * p + (m^2 - l*(l+1) + (l+1)^2 * x^2) * p1) * invY^2 * eim - harmonics_multipole_acceptance_criterion_2[1,lpm] = rhol * ((m-l-1) * x * p + (m^2 - l*(l+1) + (l+1)^2 * x^2) * p1) * invY^2 * eim[1] - harmonics_multipole_acceptance_criterion_2[2,lpm] = rhol * ((m-l-1) * x * p + (m^2 - l*(l+1) + (l+1)^2 * x^2) * p1) * invY^2 * eim[2] + #harmonics_theta[lpm] = rhol * ((l - m + 1) * p - (l + 1) * x * p1) * invY * eim + harmonics_theta[1,lpm] = rhol * ((l - m + 1) * p - (l + 1) * x * p1) * invY * eim[1] + harmonics_theta[2,lpm] = rhol * ((l - m + 1) * p - (l + 1) * x * p1) * invY * eim[2] + #harmonics_theta_2[lpm] = rhol * ((m-l-1) * x * p + (m^2 - l*(l+1) + (l+1)^2 * x^2) * p1) * invY^2 * eim + harmonics_theta_2[1,lpm] = rhol * ((m-l-1) * x * p + (m^2 - l*(l+1) + (l+1)^2 * x^2) * p1) * invY^2 * eim[1] + harmonics_theta_2[2,lpm] = rhol * ((m-l-1) * x * p + (m^2 - l*(l+1) + (l+1)^2 * x^2) * p1) * invY^2 * eim[2] rhol *= rho end rhom /= -(2 * m + 2) * (2 * m + 1) @@ -149,8 +160,8 @@ function regular_harmonic!(harmonics, harmonics_multipole_acceptance_criterion, end end -function regular_harmonic!(harmonics, rho, multipole_acceptance_criterion, phi::TF, P) where TF - y,x = sincos(multipole_acceptance_criterion) +function regular_harmonic!(harmonics, rho, theta, phi::TF, P) where TF + y,x = sincos(theta) fact = 1.0 pl = 1.0 rhom = 1.0 # rho^l / (l+m)! * (-1)^l @@ -195,8 +206,8 @@ function regular_harmonic!(harmonics, rho, multipole_acceptance_criterion, phi:: end end -function irregular_harmonic!(harmonics, rho, multipole_acceptance_criterion, phi::TF, P) where TF - y, x = sincos(multipole_acceptance_criterion) +function irregular_harmonic!(harmonics, rho, theta, phi::TF, P) where TF + y, x = sincos(theta) fact = 1 pl = 1 invR = -1.0 / rho @@ -254,8 +265,8 @@ end function M2B!(target_potential, target, center, irregular_harmonics, multipole_expansion, expansion_order) dx = target[1:3] - center - r, multipole_acceptance_criterion, phi = cartesian_2_spherical(dx) - irregular_harmonic!(irregular_harmonics, r, multipole_acceptance_criterion, phi, expansion_order) + r, theta, phi = cartesian_2_spherical(dx) + irregular_harmonic!(irregular_harmonics, r, theta, phi, expansion_order) d_potential = zeros(4) for l in 0:expansion_order for m in 0:l @@ -273,8 +284,8 @@ end function M2B!(target_potential, target, center, irregular_harmonics, multipole_expansion::Matrix{Complex{Float64}}, expansion_order) dx = target[1:3] - center - r, multipole_acceptance_criterion, phi = cartesian_2_spherical(dx) - irregular_harmonic!(irregular_harmonics, r, multipole_acceptance_criterion, phi, expansion_order) + r, theta, phi = cartesian_2_spherical(dx) + irregular_harmonic!(irregular_harmonics, r, theta, phi, expansion_order) d_potential = zeros(4) for l in 0:expansion_order for m in 0:l @@ -299,8 +310,8 @@ end function M2M!(branch, child, harmonics, M, expansion_order::Val{P}) where P # get distance vector dx, dy, dz = branch.center - child.center - r, multipole_acceptance_criterion, phi = cartesian_2_spherical(dx, dy, dz) - regular_harmonic!(harmonics, r, multipole_acceptance_criterion, phi, P) + r, theta, phi = cartesian_2_spherical(dx, dy, dz) + regular_harmonic!(harmonics, r, theta, phi, P) for j in 0:P # iterate over new Multipole coefficients B_j^k for k in 0:j @@ -385,23 +396,23 @@ end function M2L!(target_branch, source_branch, harmonics, L, expansion_order::Val{P}) where P dx, dy, dz = target_branch.center - source_branch.center - r, multipole_acceptance_criterion, phi = cartesian_2_spherical(dx, dy, dz) - irregular_harmonic!(harmonics, r, multipole_acceptance_criterion, phi, P<<1) + r, theta, phi = cartesian_2_spherical(dx, dy, dz) + irregular_harmonic!(harmonics, r, theta, phi, P<<1) M2L_loop!(target_branch.local_expansion, L, source_branch.multipole_expansion, harmonics, expansion_order) end function M2L!(target_branch, source_branch, expansion_order::Val{P}) where P dx, dy, dz = target_branch.center - source_branch.center - r, multipole_acceptance_criterion, phi = cartesian_2_spherical(dx, dy, dz) - irregular_harmonic!(target_branch.harmonics, r, multipole_acceptance_criterion, phi, P<<1) + r, theta, phi = cartesian_2_spherical(dx, dy, dz) + irregular_harmonic!(target_branch.harmonics, r, theta, phi, P<<1) M2L_loop!(target_branch.local_expansion, target_branch.ML, source_branch.multipole_expansion, target_branch.harmonics, expansion_order) end function B2L!(tree::Tree{<:Any,P}, i_branch, source_position, source_strength) where P branch = tree.branches[i_branch] irregular_harmonics = Matrix{eltype(branch.multipole_expansion)}(undef, 2, (P+1)^2) - r, multipole_acceptance_criterion, phi = cartesian_2_spherical(source_position - branch.center) - irregular_harmonic!(irregular_harmonics, r, multipole_acceptance_criterion, -phi, P) + r, theta, phi = cartesian_2_spherical(source_position - branch.center) + irregular_harmonic!(irregular_harmonics, r, theta, -phi, P) for l in 0:P for m in 0:l i_abb = (l * (l+1)) >> 1 + m + 1 @@ -417,8 +428,8 @@ end function L2L!(branch, child, regular_harmonics, L, expansion_order::Val{P}) where P dx, dy, dz = child.center - branch.center - r, multipole_acceptance_criterion, phi = cartesian_2_spherical(dx, dy, dz) - regular_harmonic!(regular_harmonics, r, multipole_acceptance_criterion, phi, P) + r, theta, phi = cartesian_2_spherical(dx, dy, dz) + regular_harmonic!(regular_harmonics, r, theta, phi, P) for j in 0:P for k in 0:j jks = (j * (j + 1)) >> 1 + k + 1 @@ -470,104 +481,127 @@ end # end # end -function L2B!(systems, branch::MultiBranch, expansion_order, vector_potential, potential_jacobian, potential_hessian, harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, workspace) - for (i_system, system) in enumerate(systems) - L2B!(system, branch.bodies_index[i_system], branch.local_expansion, expansion_order, branch.center, vector_potential, potential_jacobian, potential_hessian, harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, workspace) +function L2B!(systems, branch::MultiBranch, derivatives_switches, expansion_order, vector_potential, potential_jacobian, potential_hessian, harmonics, harmonics_theta, harmonics_theta_2, workspace) + for i in eachindex(systems) + L2B!(systems[i], branch.bodies_index[i], branch.local_expansion, derivatives_switches[i], expansion_order, branch.center, vector_potential, potential_jacobian, potential_hessian, harmonics, harmonics_theta, harmonics_theta_2, workspace) end end -function L2B!(system, branch::SingleBranch, expansion_order, vector_potential, potential_jacobian, potential_hessian, harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, workspace) - L2B!(system, branch.bodies_index, branch.local_expansion, expansion_order, branch.center, vector_potential, potential_jacobian, potential_hessian, harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, workspace) +function L2B!(system, branch::SingleBranch, derivatives_switch, expansion_order, vector_potential, potential_jacobian, potential_hessian, harmonics, harmonics_theta, harmonics_theta_2, workspace) + L2B!(system, branch.bodies_index, branch.local_expansion, derivatives_switch, expansion_order, branch.center, vector_potential, potential_jacobian, potential_hessian, harmonics, harmonics_theta, harmonics_theta_2, workspace) end -function L2B!(system, bodies_index, local_expansion, expansion_order, expansion_center, vector_potential, potential_jacobian, potential_hessian, harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, workspace) +function L2B!(system, bodies_index, local_expansion, derivatives_switch::DerivativesSwitch{PS,VPS,VS,GS}, expansion_order, expansion_center, vector_potential, potential_jacobian, potential_hessian, harmonics, harmonics_theta, harmonics_theta_2, workspace) where {PS,VPS,VS,GS} for i_body in bodies_index vector_potential .= zero(eltype(vector_potential)) potential_jacobian .= zero(eltype(potential_jacobian)) potential_hessian .= zero(eltype(potential_hessian)) body_position = system[i_body,POSITION] - scalar_potential = L2B_loop!(vector_potential, potential_jacobian, potential_hessian, body_position, expansion_center, local_expansion, harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, expansion_order, workspace) - system[i_body,SCALAR_POTENTIAL] += scalar_potential - # note: system[i,VECTOR_POTENTIAL], system[i,VELOCITY], and system[i,VELOCITY_GRADIENT] must be mutable - vpx, vpy, vpz = system[i_body,VECTOR_POTENTIAL] - system[i_body,VECTOR_POTENTIAL] = SVector{3}(vpx+vector_potential[1],vpy+vector_potential[2],vpz+vector_potential[3]) - vpx, vpy, vpz = system[i_body,VELOCITY] - system[i_body,VELOCITY] = SVector{3}(potential_jacobian[1,1]+vpx, potential_jacobian[2,1]+vpy, potential_jacobian[3,1]+vpz) - v1, v2, v3, v4, v5, v6, v7, v8, v9 = system[i_body,VELOCITY_GRADIENT] - system[i_body,VELOCITY_GRADIENT] = SMatrix{3,3}( - potential_hessian[1] + v1, - potential_hessian[2] + v2, - potential_hessian[3] + v3, - potential_hessian[4] + v4, - potential_hessian[5] + v5, - potential_hessian[6] + v6, - potential_hessian[7] + v7, - potential_hessian[8] + v8, - potential_hessian[9] + v9 - ) + scalar_potential = L2B_loop!(vector_potential, potential_jacobian, potential_hessian, body_position, expansion_center, local_expansion, harmonics, harmonics_theta, harmonics_theta_2, derivatives_switch, expansion_order, workspace) + if PS + system[i_body,SCALAR_POTENTIAL] += scalar_potential + end + if VPS + # note: system[i,VECTOR_POTENTIAL], system[i,VELOCITY], and system[i,VELOCITY_GRADIENT] must be mutable + vpx, vpy, vpz = system[i_body,VECTOR_POTENTIAL] + system[i_body,VECTOR_POTENTIAL] = SVector{3}(vpx+vector_potential[1],vpy+vector_potential[2],vpz+vector_potential[3]) + end + if VS + vpx, vpy, vpz = system[i_body,VELOCITY] + system[i_body,VELOCITY] = SVector{3}(potential_jacobian[1,1]+vpx, potential_jacobian[2,1]+vpy, potential_jacobian[3,1]+vpz) + end + if GS + v1, v2, v3, v4, v5, v6, v7, v8, v9 = system[i_body,VELOCITY_GRADIENT] + system[i_body,VELOCITY_GRADIENT] = SMatrix{3,3}( + potential_hessian[1] + v1, + potential_hessian[2] + v2, + potential_hessian[3] + v3, + potential_hessian[4] + v4, + potential_hessian[5] + v5, + potential_hessian[6] + v6, + potential_hessian[7] + v7, + potential_hessian[8] + v8, + potential_hessian[9] + v9 + ) + end end end -function L2B_loop!(vector_potential, potential_jacobian, potential_hessian, body_position, expansion_center, local_expansion, harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, expansion_order::Val{P}, workspace) where P +function L2B_loop!(vector_potential, potential_jacobian, potential_hessian, body_position, expansion_center, local_expansion, harmonics, harmonics_theta, harmonics_theta_2, derivatives_switch::DerivativesSwitch{PS,VPS,VS,GS}, expansion_order::Val{P}, workspace) where {PS,VPS,VS,GS,P} dx, dy, dz = body_position - expansion_center - r, multipole_acceptance_criterion, phi = cartesian_2_spherical(dx, dy, dz) - regular_harmonic!(harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, r, multipole_acceptance_criterion, phi, P) + r, theta, phi = cartesian_2_spherical(dx, dy, dz) + regular_harmonic!(harmonics, harmonics_theta, harmonics_theta_2, r, theta, phi, P) scalar_potential = zero(eltype(vector_potential)) for n in 0:P # nm = n * n + n + 1 # m = 0 nms = (n * (n+1)) >> 1 + 1 # m = 0 - #scalar_potential += real(local_expansion[1,nms] * harmonics[nms]) - scalar_potential += local_expansion[1,1,nms] * harmonics[1,nms] - local_expansion[2,1,nms] * harmonics[2,nms] - #vector_potential[1] += real(local_expansion[2,nms] * harmonics[nms]) - #vector_potential[2] += real(local_expansion[3,nms] * harmonics[nms]) - #vector_potential[3] += real(local_expansion[4,nms] * harmonics[nms]) - vector_potential[1] += local_expansion[1,2,nms]*harmonics[1,nms] - local_expansion[2,2,nms]*harmonics[2,nms] - vector_potential[2] += local_expansion[1,3,nms]*harmonics[1,nms] - local_expansion[2,3,nms]*harmonics[2,nms] - vector_potential[3] += local_expansion[1,4,nms]*harmonics[1,nms] - local_expansion[2,4,nms]*harmonics[2,nms] + if PS + #scalar_potential += real(local_expansion[1,nms] * harmonics[nms]) + scalar_potential += local_expansion[1,1,nms] * harmonics[1,nms] - local_expansion[2,1,nms] * harmonics[2,nms] + #vector_potential[1] += real(local_expansion[2,nms] * harmonics[nms]) + #vector_potential[2] += real(local_expansion[3,nms] * harmonics[nms]) + #vector_potential[3] += real(local_expansion[4,nms] * harmonics[nms]) + end + if VPS + vector_potential[1] += local_expansion[1,2,nms]*harmonics[1,nms] - local_expansion[2,2,nms]*harmonics[2,nms] + vector_potential[2] += local_expansion[1,3,nms]*harmonics[1,nms] - local_expansion[2,3,nms]*harmonics[2,nms] + vector_potential[3] += local_expansion[1,4,nms]*harmonics[1,nms] - local_expansion[2,4,nms]*harmonics[2,nms] + end for ind in 1:4 - # store derivatives of the potential in spherical coordinates here - #potential_jacobian[1,ind] += n/r * real(local_expansion[ind,nms] * harmonics[nms]) # dPsi/dr - #potential_jacobian[2,ind] += real(local_expansion[ind,nms] * harmonics_multipole_acceptance_criterion[nms]) # dPsi/dmultipole_acceptance_criterion - potential_jacobian[1,ind] += n/r * (local_expansion[1,ind,nms] * harmonics[1,nms] - local_expansion[2,ind,nms] * harmonics[2,nms]) # dPsi/dr - potential_jacobian[2,ind] += local_expansion[1,ind,nms] * harmonics_multipole_acceptance_criterion[1,nms] - local_expansion[2,ind,nms] * harmonics_multipole_acceptance_criterion[2,nms] # dPsi/dmultipole_acceptance_criterion - # dJ_potential[3,ind] += 0 # dPsi/dphi - potential_hessian[1,1,ind] += n * (n-1) / r^2 * (local_expansion[1,ind,nms] * harmonics[1,nms] - local_expansion[2,ind,nms] * harmonics[2,nms]) # d2Psi/dr2 - potential_hessian[2,1,ind] += n/r * (local_expansion[1,ind,nms] * harmonics_multipole_acceptance_criterion[1,nms] - local_expansion[2,ind,nms] * harmonics_multipole_acceptance_criterion[2,nms]) # d2Psi/dmultipole_acceptance_criterion dr - # potential_hessian[3,1,ind] += 0 # d2Psi/dphi dr - potential_hessian[1,2,ind] += n/r * (local_expansion[1,ind,nms] * harmonics_multipole_acceptance_criterion[1,nms] - local_expansion[2,ind,nms] * harmonics_multipole_acceptance_criterion[2,nms]) # d2Psi/dr dmultipole_acceptance_criterion - potential_hessian[2,2,ind] += local_expansion[1,ind,nms] * harmonics_multipole_acceptance_criterion_2[1,nms] - local_expansion[2,ind,nms] * harmonics_multipole_acceptance_criterion_2[2,nms] # d2Psi/dmultipole_acceptance_criterion2 - # potential_hessian[3,2,ind] += 0 # d2Psi/dphi dmultipole_acceptance_criterion - # potential_hessian[1,3,ind] += 0 # d2Psi/dr dphi - # potential_hessian[2,3,ind] += 0 # d2Psi/dmultipole_acceptance_criterion dphi - # potential_hessian[3,3,ind] += 0 # d2Psi/dphi2 - + if VS + # store derivatives of the potential in spherical coordinates here + #potential_jacobian[1,ind] += n/r * real(local_expansion[ind,nms] * harmonics[nms]) # dPsi/dr + #potential_jacobian[2,ind] += real(local_expansion[ind,nms] * harmonics_theta[nms]) # dPsi/dtheta + potential_jacobian[1,ind] += n/r * (local_expansion[1,ind,nms] * harmonics[1,nms] - local_expansion[2,ind,nms] * harmonics[2,nms]) # dPsi/dr + potential_jacobian[2,ind] += local_expansion[1,ind,nms] * harmonics_theta[1,nms] - local_expansion[2,ind,nms] * harmonics_theta[2,nms] # dPsi/dtheta + # dJ_potential[3,ind] += 0 # dPsi/dphi + end + if GS + potential_hessian[1,1,ind] += n * (n-1) / r^2 * (local_expansion[1,ind,nms] * harmonics[1,nms] - local_expansion[2,ind,nms] * harmonics[2,nms]) # d2Psi/dr2 + potential_hessian[2,1,ind] += n/r * (local_expansion[1,ind,nms] * harmonics_theta[1,nms] - local_expansion[2,ind,nms] * harmonics_theta[2,nms]) # d2Psi/dtheta dr + # potential_hessian[3,1,ind] += 0 # d2Psi/dphi dr + potential_hessian[1,2,ind] += n/r * (local_expansion[1,ind,nms] * harmonics_theta[1,nms] - local_expansion[2,ind,nms] * harmonics_theta[2,nms]) # d2Psi/dr dtheta + potential_hessian[2,2,ind] += local_expansion[1,ind,nms] * harmonics_theta_2[1,nms] - local_expansion[2,ind,nms] * harmonics_theta_2[2,nms] # d2Psi/dtheta2 + # potential_hessian[3,2,ind] += 0 # d2Psi/dphi dtheta + # potential_hessian[1,3,ind] += 0 # d2Psi/dr dphi + # potential_hessian[2,3,ind] += 0 # d2Psi/dtheta dphi + # potential_hessian[3,3,ind] += 0 # d2Psi/dphi2 + end end for m in 1:n # m > 0 # nm = n * n + n + m + 1 nms = (n * (n + 1)) >> 1 + m + 1 - scalar_potential += 2 * (local_expansion[1,1,nms] * harmonics[1,nms] - local_expansion[2,1,nms] * harmonics[2,nms]) - vector_potential[1] += 2 * (local_expansion[1,2,nms] * harmonics[1,nms] - local_expansion[2,2,nms] * harmonics[2,nms]) - vector_potential[2] += 2 * (local_expansion[1,3,nms] * harmonics[1,nms] - local_expansion[2,3,nms] * harmonics[2,nms]) - vector_potential[3] += 2 * (local_expansion[1,4,nms] * harmonics[1,nms] - local_expansion[2,4,nms] * harmonics[2,nms]) + if PS + scalar_potential += 2 * (local_expansion[1,1,nms] * harmonics[1,nms] - local_expansion[2,1,nms] * harmonics[2,nms]) + end + if VPS + vector_potential[1] += 2 * (local_expansion[1,2,nms] * harmonics[1,nms] - local_expansion[2,2,nms] * harmonics[2,nms]) + vector_potential[2] += 2 * (local_expansion[1,3,nms] * harmonics[1,nms] - local_expansion[2,3,nms] * harmonics[2,nms]) + vector_potential[3] += 2 * (local_expansion[1,4,nms] * harmonics[1,nms] - local_expansion[2,4,nms] * harmonics[2,nms]) + end for ind in 1:4 - # store derivatives of the potential in spherical harmonics here - potential_jacobian[1,ind] += 2 * n/r * (local_expansion[1,ind,nms] * harmonics[1,nms] - local_expansion[2,ind,nms] * harmonics[2,nms]) # dPsi/dr - potential_jacobian[2,ind] += 2 * (local_expansion[1,ind,nms] * harmonics_multipole_acceptance_criterion[1,nms] - local_expansion[2,ind,nms] * harmonics_multipole_acceptance_criterion[2,nms]) # dPsi/dmultipole_acceptance_criterion - #potential_jacobian[3,ind] += 2 * m * real(im * local_expansion[ind,nms] * harmonics[nms]) # dPsi/dphi - potential_jacobian[3,ind] += -2 * m * (local_expansion[1,ind,nms] * harmonics[2,nms] + local_expansion[2,ind,nms] * harmonics[1,nms]) # dPsi/dphi - potential_hessian[1,1,ind] += 2 * n * (n-1) / r^2 * (local_expansion[1,ind,nms] * harmonics[1,nms] - local_expansion[2,ind,nms] * harmonics[2,nms]) # d2Psi/dr2 - potential_hessian[2,1,ind] += 2 * n/r * (local_expansion[1,ind,nms] * harmonics_multipole_acceptance_criterion[1,nms] - local_expansion[2,ind,nms] * harmonics_multipole_acceptance_criterion[2,nms]) # d2Psi/dmultipole_acceptance_criterion dr - potential_hessian[3,1,ind] += -2 * n * m / r * (local_expansion[1,ind,nms] * harmonics[2,nms] + local_expansion[2,ind,nms] * harmonics[1,nms]) # d2Psi/dphi dr - potential_hessian[1,2,ind] += 2 * n/r * (local_expansion[1,ind,nms] * harmonics_multipole_acceptance_criterion[1,nms] - local_expansion[2,ind,nms] * harmonics_multipole_acceptance_criterion[2,nms]) # d2Psi/dr dmultipole_acceptance_criterion - potential_hessian[2,2,ind] += 2 * (local_expansion[1,ind,nms] * harmonics_multipole_acceptance_criterion_2[1,nms] - local_expansion[2,ind,nms] * harmonics_multipole_acceptance_criterion_2[2,nms]) # d2Psi/dmultipole_acceptance_criterion2 - potential_hessian[3,2,ind] += -2 * m * (local_expansion[1,ind,nms] * harmonics_multipole_acceptance_criterion[2,nms] + local_expansion[2,ind,nms] * harmonics_multipole_acceptance_criterion[1,nms]) # d2Psi/dphi dmultipole_acceptance_criterion - potential_hessian[1,3,ind] += -2 * n * m / r * (local_expansion[1,ind,nms] * harmonics[2,nms] + local_expansion[2,ind,nms] * harmonics[1,nms]) # d2Psi/dr dphi - potential_hessian[2,3,ind] += -2 * m * (local_expansion[1,ind,nms] * harmonics_multipole_acceptance_criterion[2,nms] + local_expansion[2,ind,nms] * harmonics_multipole_acceptance_criterion[1,nms]) # d2Psi/dmultipole_acceptance_criterion dphi - potential_hessian[3,3,ind] += 2 * -m^2 * (local_expansion[1,ind,nms] * harmonics[2,nms] - local_expansion[2,ind,nms] * harmonics[1,nms]) # d2Psi/dphi2 + if VS + # store derivatives of the potential in spherical harmonics here + potential_jacobian[1,ind] += 2 * n/r * (local_expansion[1,ind,nms] * harmonics[1,nms] - local_expansion[2,ind,nms] * harmonics[2,nms]) # dPsi/dr + potential_jacobian[2,ind] += 2 * (local_expansion[1,ind,nms] * harmonics_theta[1,nms] - local_expansion[2,ind,nms] * harmonics_theta[2,nms]) # dPsi/dtheta + #potential_jacobian[3,ind] += 2 * m * real(im * local_expansion[ind,nms] * harmonics[nms]) # dPsi/dphi + potential_jacobian[3,ind] += -2 * m * (local_expansion[1,ind,nms] * harmonics[2,nms] + local_expansion[2,ind,nms] * harmonics[1,nms]) # dPsi/dphi + end + if GS + potential_hessian[1,1,ind] += 2 * n * (n-1) / r^2 * (local_expansion[1,ind,nms] * harmonics[1,nms] - local_expansion[2,ind,nms] * harmonics[2,nms]) # d2Psi/dr2 + potential_hessian[2,1,ind] += 2 * n/r * (local_expansion[1,ind,nms] * harmonics_theta[1,nms] - local_expansion[2,ind,nms] * harmonics_theta[2,nms]) # d2Psi/dtheta dr + potential_hessian[3,1,ind] += -2 * n * m / r * (local_expansion[1,ind,nms] * harmonics[2,nms] + local_expansion[2,ind,nms] * harmonics[1,nms]) # d2Psi/dphi dr + potential_hessian[1,2,ind] += 2 * n/r * (local_expansion[1,ind,nms] * harmonics_theta[1,nms] - local_expansion[2,ind,nms] * harmonics_theta[2,nms]) # d2Psi/dr dtheta + potential_hessian[2,2,ind] += 2 * (local_expansion[1,ind,nms] * harmonics_theta_2[1,nms] - local_expansion[2,ind,nms] * harmonics_theta_2[2,nms]) # d2Psi/dtheta2 + potential_hessian[3,2,ind] += -2 * m * (local_expansion[1,ind,nms] * harmonics_theta[2,nms] + local_expansion[2,ind,nms] * harmonics_theta[1,nms]) # d2Psi/dphi dtheta + potential_hessian[1,3,ind] += -2 * n * m / r * (local_expansion[1,ind,nms] * harmonics[2,nms] + local_expansion[2,ind,nms] * harmonics[1,nms]) # d2Psi/dr dphi + potential_hessian[2,3,ind] += -2 * m * (local_expansion[1,ind,nms] * harmonics_theta[2,nms] + local_expansion[2,ind,nms] * harmonics_theta[1,nms]) # d2Psi/dtheta dphi + potential_hessian[3,3,ind] += 2 * -m^2 * (local_expansion[1,ind,nms] * harmonics[2,nms] - local_expansion[2,ind,nms] * harmonics[1,nms]) # d2Psi/dphi2 + end end end end - spherical_2_cartesian!(potential_jacobian, potential_hessian, workspace, r, multipole_acceptance_criterion, phi) - flatten_derivatives!(potential_jacobian, potential_hessian) # compute velocity and velocity gradient + spherical_2_cartesian!(potential_jacobian, potential_hessian, workspace, r, theta, phi, derivatives_switch) + flatten_derivatives!(potential_jacobian, potential_hessian, derivatives_switch) # compute velocity and velocity gradient return scalar_potential end diff --git a/src/tree.jl b/src/tree.jl index 6e43ccc..8310b52 100644 --- a/src/tree.jl +++ b/src/tree.jl @@ -659,6 +659,7 @@ function shrink_leaf!(branch, system) # recenter # turn this off for now- only shrink/expand radius # new_center = center_nonzero_radius(system, bodies_index) new_center = branch[].center + length(bodies_index) == 1 && (new_center += SVector{3}(sqrt(eps()),sqrt(eps()),sqrt(eps()))) # shrink radius new_radius = zero(branch[].radius) diff --git a/test/gravitational.jl b/test/gravitational.jl index 277b8d9..bef8d29 100644 --- a/test/gravitational.jl +++ b/test/gravitational.jl @@ -63,7 +63,7 @@ fmm.buffer_element(g::Gravitational) = (deepcopy(g.bodies[1]),zeros(eltype(g),52 fmm.B2M!(system::Gravitational, args...) = fmm.B2M!_sourcepoint(system, args...) -function fmm.direct!(target_system, target_index, source_system::Gravitational, source_index) +function fmm.direct!(target_system, target_index, derivatives_switch, source_system::Gravitational, source_index) # nbad = 0 for i_source in source_index source_x, source_y, source_z = source_system[i_source,fmm.POSITION] diff --git a/test/runtests.jl b/test/runtests.jl index 5182f0f..9470a36 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,8 +19,8 @@ test_dir = @__DIR__ include(joinpath(test_dir, "gravitational.jl")) fmm.M2M!(branch, child, harmonics, M, expansion_order::Int) = fmm.M2M!(branch, child, harmonics, M, Val(expansion_order)) -fmm.L2B_loop!(vector_potential, potential_jacobian, potential_hessian, body_position, expansion_center, local_expansion, harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, expansion_order::Int, workspace) = - fmm.L2B_loop!(vector_potential, potential_jacobian, potential_hessian, body_position, expansion_center, local_expansion, harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, Val(expansion_order), workspace) +fmm.L2B_loop!(vector_potential, potential_jacobian, potential_hessian, body_position, expansion_center, local_expansion, harmonics, harmonics_theta, harmonics_theta_2, expansion_order::Int, workspace) = + fmm.L2B_loop!(vector_potential, potential_jacobian, potential_hessian, body_position, expansion_center, local_expansion, harmonics, harmonics_theta, harmonics_theta_2, Val(expansion_order), workspace) fmm.M2L!(target_branch, source_branch, harmonics, L, expansion_order::Int) = fmm.M2L!(target_branch, source_branch, harmonics, L, Val(expansion_order)) fmm.upward_pass_singlethread!(branches, systems, expansion_order::Int) = fmm.upward_pass_singlethread!(branches, systems, Val(expansion_order)) fmm.M2L!(target_branch, source_branch, expansion_order::Int) = fmm.M2L!(target_branch, source_branch, Val(expansion_order)) @@ -160,12 +160,12 @@ end @testset "cartesian to spherical" begin # cartesian to spherical rho = 1.0 - multipole_acceptance_criterion = pi/4 + theta = pi/4 phi = pi/2 - that = [rho, multipole_acceptance_criterion, phi] - x = rho * sin(multipole_acceptance_criterion) * cos(phi) - y = rho * sin(multipole_acceptance_criterion) * sin(phi) - z = rho * cos(multipole_acceptance_criterion) + that = [rho, theta, phi] + x = rho * sin(theta) * cos(phi) + y = rho * sin(theta) * sin(phi) + z = rho * cos(theta) this = [x,y,z] fmm.cartesian_2_spherical!(this) for i in 1:3 @@ -178,31 +178,31 @@ end # once I begin using multiple kernels @testset "solid harmonics" begin -function Ylm(multipole_acceptance_criterion, phi, l, m) +function Ylm(theta, phi, l, m) lm! = sqrt(factorial(big(l-abs(m)))/ factorial(big(l+abs(m)))) - plm = Plm(cos(multipole_acceptance_criterion), l, abs(m)) + plm = Plm(cos(theta), l, abs(m)) eim = exp(im * m * phi) ylm = lm! * plm * eim end -function regular_harmonic_manual(rho, multipole_acceptance_criterion, phi, p) +function regular_harmonic_manual(rho, theta, phi, p) reg_harmonics = Vector{Complex{Float64}}(undef,(p+1)^2) i = 1 for l in 0:p for m in -l:l - reg_harmonics[i] = Ylm(multipole_acceptance_criterion, phi, l, m) * rho^l + reg_harmonics[i] = Ylm(theta, phi, l, m) * rho^l i+=1 end end return reg_harmonics end -function irregular_harmonic_manual(rho, multipole_acceptance_criterion, phi, p) +function irregular_harmonic_manual(rho, theta, phi, p) reg_harmonics = Vector{Complex{Float64}}(undef,(p+1)^2) i = 1 for l in 0:p for m in -l:l - reg_harmonics[i] = Ylm(multipole_acceptance_criterion, phi, l, m) / rho^(l+1) + reg_harmonics[i] = Ylm(theta, phi, l, m) / rho^(l+1) i+=1 end end @@ -217,8 +217,8 @@ P = 3 rh_man = regular_harmonic_manual(rho, alpha, beta, P) rh_fmm = zeros(Complex{Float64},length(rh_man)) -rh_fmm_multipole_acceptance_criterion = zeros(Complex{Float64},length(rh_man)) -fmm.regular_harmonic!(rh_fmm, rh_fmm_multipole_acceptance_criterion, rho, alpha, beta, P) +rh_fmm_theta = zeros(Complex{Float64},length(rh_man)) +fmm.regular_harmonic!(rh_fmm, rh_fmm_theta, rho, alpha, beta, P) for i in 1:length(rh_man) @test isapprox(rh_man[i], rh_fmm[i]; atol=1e-11) @@ -287,8 +287,8 @@ u_fmm = target_potential[1] dx = x_target - xs[1,:] u_check = ms[1] / sqrt(dx' * dx) -function Ylm(multipole_acceptance_criterion, phi, l, m) - ylm = sqrt(factorial(big(l-abs(m)))/ factorial(big(l+abs(m)))) * Plm(cos(multipole_acceptance_criterion), l, abs(m)) * exp(im * m * phi) +function Ylm(theta, phi, l, m) + ylm = sqrt(factorial(big(l-abs(m)))/ factorial(big(l+abs(m)))) * Plm(cos(theta), l, abs(m)) * exp(im * m * phi) end function evaluate_biot_savart(x_source, x_target, q_source, P) @@ -403,9 +403,9 @@ source_i = new_order_index[1] # just needs to be farther away than the target to dx_source = fmm.cartesian_2_spherical(elements[source_i,fmm.POSITION] - tree.branches[branch_i].center) dx_target = fmm.cartesian_2_spherical(elements[target_i,fmm.POSITION] - tree.branches[branch_i].center) -local_coefficients_multipole_acceptance_criterion = zeros(Float64, 2, ((expansion_order+1)*(expansion_order+2))>>1) +local_coefficients_theta = zeros(Float64, 2, ((expansion_order+1)*(expansion_order+2))>>1) local_coefficients_expanded = zeros(Float64, 2, (expansion_order+1)^2) -local_coefficients_expanded_multipole_acceptance_criterion = zeros(Float64, 2, (expansion_order+1)^2) +local_coefficients_expanded_theta = zeros(Float64, 2, (expansion_order+1)^2) fmm.irregular_harmonic!(local_coefficients_expanded, dx_source..., expansion_order) local_coefficients_expanded .*= ms[1] regular_harmonics_expanded = zeros(Float64, 2, (expansion_order+1)^2) @@ -414,10 +414,10 @@ fmm.regular_harmonic!(regular_harmonics_expanded, dx_target..., expansion_order) fmm.B2L!(tree, branch_i, elements[source_i,fmm.POSITION], elements.bodies[source_i].strength) harmonics = zeros(Float64, 2, (expansion_order+1)^2) -harmonics_multipole_acceptance_criterion = zeros(Float64, 2, (expansion_order+1)^2) -harmonics_multipole_acceptance_criterion_2 = zeros(Float64, 2, (expansion_order+1)^2) +harmonics_theta = zeros(Float64, 2, (expansion_order+1)^2) +harmonics_theta_2 = zeros(Float64, 2, (expansion_order+1)^2) workspace = zeros(3,4) -fmm.L2B!(elements, target_i:target_i, tree.branches[branch_i].local_expansion, expansion_order, tree.branches[branch_i].center, zeros(3), zeros(3,4), zeros(3,3,4), harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, workspace) +fmm.L2B!(elements, target_i:target_i, tree.branches[branch_i].local_expansion, fmm.DerivativesSwitch(), Val(expansion_order), tree.branches[branch_i].center, zeros(3), zeros(3,4), zeros(3,3,4), harmonics, harmonics_theta, harmonics_theta_2, workspace) u_fmm = elements.potential[1,target_i] dx_direct = xs[4,:] - xs[1,:] @@ -465,11 +465,11 @@ fmm.B2L!(tree, 2, elements[new_order_index[1],fmm.POSITION], elements.bodies[new # check L2P now: harmonics = zeros(Float64,2,(expansion_order+1)^2) -harmonics_multipole_acceptance_criterion = zeros(Float64,2,(expansion_order+1)^2) -harmonics_multipole_acceptance_criterion_2 = zeros(Float64,2,(expansion_order+1)^2) +harmonics_theta = zeros(Float64,2,(expansion_order+1)^2) +harmonics_theta_2 = zeros(Float64,2,(expansion_order+1)^2) workspace = zeros(3,4) spherical_potential = zeros(52) -fmm.L2B!(elements, new_order_index[5]:new_order_index[5], tree.branches[2].local_expansion, expansion_order, tree.branches[2].center, zeros(3), zeros(3,4), zeros(3,3,4), harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, workspace) +fmm.L2B!(elements, new_order_index[5]:new_order_index[5], tree.branches[2].local_expansion, fmm.DerivativesSwitch(), Val(expansion_order), tree.branches[2].center, zeros(3), zeros(3,4), zeros(3,3,4), harmonics, harmonics_theta, harmonics_theta_2, workspace) u_fmm_no_x = elements.potential[1,new_order_index[5]] elements.potential[1,new_order_index[5]] *= 0 @@ -482,7 +482,7 @@ fmm.irregular_harmonic!(local_coefficients_check, dx_check, dy_check, dz_check, local_coefficients_check .*= ms[1] # evaluate local expansion at mass 5 -fmm.L2B!((elements,), tree.branches[7], expansion_order, zeros(3), zeros(3,4), zeros(3,3,4), harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, workspace) +fmm.L2B!((elements,), tree.branches[7], (fmm.DerivativesSwitch(),), Val(expansion_order), zeros(3), zeros(3,4), zeros(3,3,4), harmonics, harmonics_theta, harmonics_theta_2, workspace) u_fmm = elements.potential[1,new_order_index[5]] dx_direct = elements[new_order_index[5],fmm.POSITION] - elements[new_order_index[1],fmm.POSITION] @@ -530,8 +530,8 @@ tree = fmm.Tree((elements,); expansion_order, n_per_branch=1, shrink_recenter=fa i_branch_multipole = 7 # mass 5 i_branch_local = 5 # mass 1 harmonics = zeros(Float64,2, (expansion_order+1)^2) -harmonics_multipole_acceptance_criterion = zeros(Float64,2, (expansion_order+1)^2) -harmonics_multipole_acceptance_criterion_2 = zeros(Float64,2, (expansion_order+1)^2) +harmonics_theta = zeros(Float64,2, (expansion_order+1)^2) +harmonics_theta_2 = zeros(Float64,2, (expansion_order+1)^2) workspace = zeros(3,4) fmm.B2M!(elements, tree.branches[i_branch_multipole], new_order_index[5]:new_order_index[5], harmonics, tree.expansion_order) @@ -551,7 +551,7 @@ fmm.B2M!(elements, tree.branches[i_branch_multipole], new_order_index[5]:new_ord m2l_harmonics = zeros(eltype(tree.branches[1].multipole_expansion), 2, (expansion_order<<1 + 1)*(expansion_order<<1 + 1)) L = zeros(eltype(tree.branches[1].local_expansion), 2, 4) fmm.M2L!(tree.branches[i_branch_local], tree.branches[i_branch_multipole], m2l_harmonics, L, expansion_order) -fmm.L2B!(elements, new_order_index[1]:new_order_index[1], tree.branches[i_branch_local].local_expansion, expansion_order, tree.branches[i_branch_local].center, zeros(3), zeros(3,4), zeros(3,3,4), harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, workspace) +fmm.L2B!(elements, new_order_index[1]:new_order_index[1], tree.branches[i_branch_local].local_expansion, fmm.DerivativesSwitch(), Val(expansion_order), tree.branches[i_branch_local].center, zeros(3), zeros(3,4), zeros(3,3,4), harmonics, harmonics_theta, harmonics_theta_2, workspace) u_fmm = elements.potential[1,new_order_index[1]] local_exp = tree.branches[i_branch_local].local_expansion[1] @@ -623,7 +623,7 @@ u_fmm_67 = mass_target_potential[1] # perform horizontal pass m2l_list, direct_list = fmm.build_interaction_lists(tree.branches, tree.branches, multipole_acceptance_criterion, true, true, true) -fmm.nearfield_singlethread!((elements,), tree.branches, (elements,), tree.branches, direct_list) +fmm.nearfield_singlethread!((elements,), tree.branches, (fmm.DerivativesSwitch(),), (elements,), tree.branches, direct_list) fmm.horizontal_pass_singlethread!(tree.branches, tree.branches, m2l_list, expansion_order) # consider the effect on branch 3 (mass 2) @@ -632,9 +632,9 @@ elements.potential[i_POTENTIAL,new_order_index[2]] .*= 0 # reset potential at ma # elements.direct!(elements.potential[i_POTENTIAL,new_order_index[2]], elements.bodies[i_POSITION,new_order_index[2]], elements.bodies[:,new_order_index[1]]) # elements.direct!(elements.potential[i_POTENTIAL,new_order_index[2]], elements.bodies[i_POSITION,new_order_index[2]], elements.bodies[:,new_order_index[2]]) # elements.direct!(elements.potential[i_POTENTIAL,new_order_index[2]], elements.bodies[i_POSITION,new_order_index[2]], elements.bodies[:,new_order_index[3]]) -fmm.P2P!((elements,), tree.branches[3], (elements,), tree.branches[3]) -fmm.P2P!((elements,), tree.branches[3], (elements,), tree.branches[4]) -fmm.P2P!((elements,), tree.branches[3], (elements,), tree.branches[5]) +fmm.P2P!((elements,), tree.branches[3], (fmm.DerivativesSwitch(),), (elements,), tree.branches[3]) +fmm.P2P!((elements,), tree.branches[3], (fmm.DerivativesSwitch(),), (elements,), tree.branches[4]) +fmm.P2P!((elements,), tree.branches[3], (fmm.DerivativesSwitch(),), (elements,), tree.branches[5]) u_fmm_123 = elements.potential[i_POTENTIAL[1],new_order_index[2]] dx_12 = elements[new_order_index[2],fmm.POSITION] - elements[new_order_index[1],fmm.POSITION] @@ -647,10 +647,10 @@ u_direct_123 = u_direct_12 + u_direct_22 + u_direct_32 # M2L is performed from branches 6, 7 to branch 3 (containing mass 2) harmonics = zeros(Float64,2, (expansion_order+1)^2) -harmonics_multipole_acceptance_criterion = zeros(Float64,2, (expansion_order+1)^2) -harmonics_multipole_acceptance_criterion_2 = zeros(Float64,2, (expansion_order+1)^2) +harmonics_theta = zeros(Float64,2, (expansion_order+1)^2) +harmonics_theta_2 = zeros(Float64,2, (expansion_order+1)^2) workspace = zeros(3,4) -fmm.L2B!((elements,), tree.branches[3], expansion_order, zeros(3), zeros(3,4), zeros(3,3,4), harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, workspace) +fmm.L2B!((elements,), tree.branches[3], (fmm.DerivativesSwitch(),), Val(expansion_order), zeros(3), zeros(3,4), zeros(3,3,4), harmonics, harmonics_theta, harmonics_theta_2, workspace) u_fmm_12345 = elements.potential[i_POTENTIAL[1],new_order_index[2]] dx_42 = elements[new_order_index[4],fmm.POSITION] - elements[new_order_index[2],fmm.POSITION] @@ -693,21 +693,21 @@ include("vortex.jl") """ 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 @@ -715,37 +715,37 @@ end function cartesian_2_spherical(x,y,z; vec=false) r = sqrt(x^2+y^2+z^2) - multipole_acceptance_criterion = acos(z/r) + theta = acos(z/r) phi = atan(y,x) - vec && return [r,multipole_acceptance_criterion,phi] - return r, multipole_acceptance_criterion, phi + vec && return [r,theta,phi] + return r, theta, phi end function d2rdx2_cart(x,y,z) - r, multipole_acceptance_criterion, phi = cartesian_2_spherical(x,y,z) - return d2rdx2(r,multipole_acceptance_criterion,phi) + r, theta, phi = cartesian_2_spherical(x,y,z) + return d2rdx2(r,theta,phi) 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 function drdx_cart(x,y,z) - r, multipole_acceptance_criterion, phi = cartesian_2_spherical(x,y,z) - return drdx(r, multipole_acceptance_criterion, phi) + r, theta, phi = cartesian_2_spherical(x,y,z) + return drdx(r, theta, phi) end function d2rdx2_fd(x,y,z) @@ -786,21 +786,21 @@ end """ 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 @@ -809,16 +809,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 @@ -826,24 +826,24 @@ end # """ # dr_i/dx_j # """ -# 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(multipole_acceptance_criterion)/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(theta)/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 -function derivatives_2_cartesian(first_derivatives, second_derivatives, r, multipole_acceptance_criterion, phi) +function derivatives_2_cartesian(first_derivatives, second_derivatives, r, theta, phi) @assert size(spherical_derivatives) == (3,3) - d2_unit_vec = d2rdx2(r, multipole_acceptance_criterion, phi) - d_unit_vec = drdx(r,multipole_acceptance_criterion,phi) + d2_unit_vec = d2rdx2(r, theta, phi) + d_unit_vec = drdx(r,theta,phi) cartesian_derivatives = zeros(3,3) for k in 1:3 cartesian_derivatives .+= d2_unit_vec[:,:,k] .* first_derivatives[k] @@ -859,23 +859,23 @@ end function simple(X) r = X[1] - multipole_acceptance_criterion = X[2] + theta = X[2] phi = X[3] - return 3*r^4*cos(multipole_acceptance_criterion) + 2*multipole_acceptance_criterion - sin(phi) + return 3*r^4*cos(theta) + 2*theta - sin(phi) end function cartesian_2_spherical(x,y,z; vec=false) r = sqrt(x^2+y^2+z^2) - multipole_acceptance_criterion = acos(z/r) + theta = acos(z/r) phi = atan(y/x) - vec && return [r,multipole_acceptance_criterion,phi] - return r, multipole_acceptance_criterion, phi + vec && return [r,theta,phi] + return r, theta, phi end -function spherical_2_cartesian(r,multipole_acceptance_criterion,phi; vec=false) - x = r*sin(multipole_acceptance_criterion)*cos(phi) - y = r*sin(multipole_acceptance_criterion)*sin(phi) - z = r*cos(multipole_acceptance_criterion) +function spherical_2_cartesian(r,theta,phi; vec=false) + x = r*sin(theta)*cos(phi) + y = r*sin(theta)*sin(phi) + z = r*cos(theta) vec && return [x,y,z] return x,y,z end @@ -885,20 +885,20 @@ function simple_cart(X) return simple(args) end -r0, multipole_acceptance_criterion0, phi0 = rand(3) +r0, theta0, phi0 = rand(3) -dsdr(r,multipole_acceptance_criterion,phi) = 12 * r^3 * cos(multipole_acceptance_criterion) -dsdt(r,multipole_acceptance_criterion,phi) = 3*r^4*-sin(multipole_acceptance_criterion) + 2 -dsdp(r,multipole_acceptance_criterion,phi) = -cos(phi) -d2sdr2(r,multipole_acceptance_criterion,phi) = 36 * r^2 * cos(multipole_acceptance_criterion) -d2sdrdt(r,multipole_acceptance_criterion,phi) = -12*r^3*sin(multipole_acceptance_criterion) -d2sdrdp(r,multipole_acceptance_criterion,phi) = 0.0 -d2sdt2(r,multipole_acceptance_criterion,phi) = -3*r^4*cos(multipole_acceptance_criterion) -d2sdtdp(r,multipole_acceptance_criterion,phi) = 0.0 -d2sdp2(r,multipole_acceptance_criterion,phi) = sin(phi) +dsdr(r,theta,phi) = 12 * r^3 * cos(theta) +dsdt(r,theta,phi) = 3*r^4*-sin(theta) + 2 +dsdp(r,theta,phi) = -cos(phi) +d2sdr2(r,theta,phi) = 36 * r^2 * cos(theta) +d2sdrdt(r,theta,phi) = -12*r^3*sin(theta) +d2sdrdp(r,theta,phi) = 0.0 +d2sdt2(r,theta,phi) = -3*r^4*cos(theta) +d2sdtdp(r,theta,phi) = 0.0 +d2sdp2(r,theta,phi) = sin(phi) -testd = cs_derivative(simple,1,Complex{Float64}[r0,multipole_acceptance_criterion0,phi0]) -@test isapprox(testd, dsdr(r0,multipole_acceptance_criterion0,phi0);atol=1e-12) +testd = cs_derivative(simple,1,Complex{Float64}[r0,theta0,phi0]) +@test isapprox(testd, dsdr(r0,theta0,phi0);atol=1e-12) function second_derivative(func,i,j,args; step=1e-8) first_func(args) = cs_derivative(func,i,args) @@ -917,11 +917,11 @@ function fd_hessian(func, args) end return simple_jacobian end -simple_jacobian = fd_hessian(simple,[r0,multipole_acceptance_criterion0,phi0]) +simple_jacobian = fd_hessian(simple,[r0,theta0,phi0]) simple_jacobian_anal = [ - d2sdr2(r0,multipole_acceptance_criterion0,phi0) d2sdrdt(r0,multipole_acceptance_criterion0,phi0) d2sdrdp(r0,multipole_acceptance_criterion0,phi0); - d2sdrdt(r0,multipole_acceptance_criterion0,phi0) d2sdt2(r0,multipole_acceptance_criterion0,phi0) d2sdtdp(r0,multipole_acceptance_criterion0,phi0); - d2sdrdp(r0,multipole_acceptance_criterion0,phi0) d2sdtdp(r0,multipole_acceptance_criterion0,phi0) d2sdp2(r0,multipole_acceptance_criterion0,phi0); + d2sdr2(r0,theta0,phi0) d2sdrdt(r0,theta0,phi0) d2sdrdp(r0,theta0,phi0); + d2sdrdt(r0,theta0,phi0) d2sdt2(r0,theta0,phi0) d2sdtdp(r0,theta0,phi0); + d2sdrdp(r0,theta0,phi0) d2sdtdp(r0,theta0,phi0) d2sdp2(r0,theta0,phi0); ] for i in 1:length(simple_jacobian) @@ -929,7 +929,7 @@ for i in 1:length(simple_jacobian) end # test chain rule -args = [r0,multipole_acceptance_criterion0,phi0] +args = [r0,theta0,phi0] spherical_grad = [cs_derivative(simple,i,convert(Vector{Complex{Float64}},args)) for i in 1:3] spherical_hessian = fd_hessian(simple, args) d2rkdxidxj = d2rdx2(args...) @@ -957,7 +957,7 @@ for i in 1:4 end workspace = zeros(3,4) -fmm.spherical_2_cartesian!(potential_jacobian, potential_hessian, workspace, r0, multipole_acceptance_criterion0, phi0) +fmm.spherical_2_cartesian!(potential_jacobian, potential_hessian, workspace, r0, theta0, phi0, fmm.DerivativesSwitch()) for i in 1:3 for j in 1:3 @@ -1031,11 +1031,11 @@ fmm.M2L!(tree.branches[2], tree.branches[3], m2l_harmonics, L, Val(expansion_ord fmm.M2L!(tree.branches[3], tree.branches[2], m2l_harmonics, L, Val(expansion_order)) # @show tree.branches[2].local_expansion harmonics = zeros(Float64,2, (expansion_order+1)^2) -harmonics_multipole_acceptance_criterion = zeros(Float64,2, (expansion_order+1)^2) -harmonics_multipole_acceptance_criterion_2 = zeros(Float64,2, (expansion_order+1)^2) +harmonics_theta = zeros(Float64,2, (expansion_order+1)^2) +harmonics_theta_2 = zeros(Float64,2, (expansion_order+1)^2) workspace = zeros(3,4) -fmm.L2B!((vortexparticles,), tree.branches[2], Val(expansion_order), zeros(3), zeros(3,4), zeros(3,3,4), harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, workspace) -fmm.L2B!((vortexparticles,), tree.branches[3], Val(expansion_order), zeros(3), zeros(3,4), zeros(3,3,4), harmonics, harmonics_multipole_acceptance_criterion, harmonics_multipole_acceptance_criterion_2, workspace) +fmm.L2B!((vortexparticles,), tree.branches[2], (fmm.DerivativesSwitch(),), Val(expansion_order), zeros(3), zeros(3,4), zeros(3,3,4), harmonics, harmonics_theta, harmonics_theta_2, workspace) +fmm.L2B!((vortexparticles,), tree.branches[3], (fmm.DerivativesSwitch(),), Val(expansion_order), zeros(3), zeros(3,4), zeros(3,3,4), harmonics, harmonics_theta, harmonics_theta_2, workspace) update_velocity_stretching!(vortexparticles) for i in 1:2 @@ -1420,7 +1420,7 @@ multipole_expansion_test = [3.15347369905358e-5 0.0 0.0 0.0; 0.0 0.0 0.0 0.0;;; expansion_order = Val{14}() fmm._B2M!_panel(multipole_expansion, qnm_prev, jnm_prev, inm_prev, R0, Ru, Rv, strength, normal, expansion_order, fmm.UniformSourcePanel()) for i in eachindex(multipole_expansion) - @test isapprox(multipole_expansion[i], multipole_expansion_test[i]; atol=2e-11) + @test isapprox(multipole_expansion[i], multipole_expansion_test[i]; atol=1e-12) end # single tri panel diff --git a/test/vortex.jl b/test/vortex.jl index 9cd37fb..6586432 100644 --- a/test/vortex.jl +++ b/test/vortex.jl @@ -73,7 +73,7 @@ Base.eltype(::VortexParticles{TF}) where TF = TF """ Classical formulation so far. """ -function fmm.direct!(target_system, target_index, source_system::VortexParticles, source_index) +function fmm.direct!(target_system, target_index, derivatives_switch, source_system::VortexParticles, source_index) jacobian = zeros(eltype(target_system),3,4) hessian = zeros(eltype(target_system),3,3,4) for j_source in source_index @@ -110,7 +110,7 @@ function fmm.direct!(target_system, target_index, source_system::VortexParticles hessian[3,2,i_POTENTIAL_VECTOR] += gamma_over_R * 3 * dx[2] * dx[3] # dydz hessian[3,3,i_POTENTIAL_VECTOR] += gamma_over_R .* (-dx[1]^2 - dx[2]^2 + 2 * dx[3]^2) # dz2 end - fmm.flatten_derivatives!(jacobian, hessian) + fmm.flatten_derivatives!(jacobian, hessian, derivatives_switch) target_system[i_target,fmm.VELOCITY] .+= view(jacobian,:,1) target_system[i_target,fmm.VELOCITY_GRADIENT] .+= view(hessian,:,:,1)