src.acoustools.Force

  1from acoustools.Gorkov import gorkov_fin_diff, get_finite_diff_points_all_axis, get_gorkov_constants
  2from acoustools.Utilities import forward_model_batched, forward_model_grad, forward_model_second_derivative_unmixed, forward_model_second_derivative_mixed, TRANSDUCERS, propagate, DTYPE, device
  3import acoustools.Constants as c
  4
  5import torch
  6from torch import Tensor
  7from types import FunctionType
  8
  9
 10def force_fin_diff(activations:Tensor, points:Tensor, axis:str="XYZ", stepsize:float= 0.000135156253,K1:float|None=None, 
 11                   K2:float|None=None,U_function:FunctionType=gorkov_fin_diff,U_fun_args:dict={}, board:Tensor|None=None, V=c.V, p_ref=c.P_ref,
 12                    k=c.k, transducer_radius = c.radius, transducer_norms = None,
 13                    medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p) -> Tensor:
 14    '''
 15    Returns the force on a particle using finite differences to approximate the derivative of the gor'kov potential\n
 16    :param activations: Transducer hologram
 17    :param points: Points to propagate to
 18    :param axis: string containing `X`, `Y` or `Z` defining the axis to take into account eg `XYZ` considers all 3 axes and `YZ` considers only the y and z-axes
 19    :param stepsize: stepsize to use for finite differences 
 20    :param K1: Value for K1 to be used in the gor'kov computation, see `Holographic acoustic elements for manipulation of levitated objects` for more information
 21    :param K2: Value for K1 to be used in the gor'kov computation, see `Holographic acoustic elements for manipulation of levitated objects` for more information
 22    :param U_function: The function used to compute the gor'kov potential
 23    :param U_fun_args: arguments for `U_function` 
 24    :param board: Transducers to use, if `None` uses `acoustools.Utilities.TRANSDUCERS`
 25    :parm V: Particle volume
 26    :return: Force
 27    '''
 28    B = points.shape[0]
 29    D = len(axis)
 30    N = points.shape[2]
 31
 32    if board is None:
 33        board = TRANSDUCERS
 34
 35    fin_diff_points = get_finite_diff_points_all_axis(points, axis, stepsize)
 36    
 37    U_points = U_function(activations, fin_diff_points, axis=axis, stepsize=stepsize/10 ,K1=K1,K2=K2,**U_fun_args, board=board,V=V, 
 38                            p_ref=p_ref, k=k, transducer_radius=transducer_radius, transducer_norms=transducer_norms,
 39                            medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density,particle_speed=particle_speed)
 40    U_grads = U_points[:,N:]
 41    split = torch.reshape(U_grads,(B,2,-1))
 42    # print(split)
 43    # print((split[:,0,:] - split[:,1,:]))
 44
 45    # print()
 46
 47    F = -1* (split[:,0,:] - split[:,1,:]) / (2*stepsize)
 48    F = F.reshape(B,3,N).permute(0,2,1)
 49    return F
 50
 51def compute_force(activations:Tensor, points:Tensor,board:Tensor|None=None,return_components:bool=False, V=c.V, p_ref=c.P_ref, 
 52                  transducer_radius=c.radius, k=c.k,
 53                 medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p, transducer_norms=None) -> Tensor | tuple[Tensor, Tensor, Tensor]:
 54    '''
 55    Returns the force on a particle using the analytical derivative of the Gor'kov potential and the piston model\n
 56    :param activations: Transducer hologram
 57    :param points: Points to propagate to
 58    :param board: Transducers to use, if `None` uses `acoustools.Utilities.TRANSDUCERS`
 59    :param return_components: If true returns force as one tensor otherwise returns Fx, Fy, Fz
 60    :param V: Particle volume
 61    :return: force  
 62    '''
 63
 64    #Bk.2 Pg.319
 65
 66    if board is None:
 67        board = TRANSDUCERS
 68
 69    if transducer_norms is None:
 70        transducer_norms = (torch.zeros_like(board) + torch.tensor([0,0,1], device=device)) * torch.sign(board[:,2].real).unsqueeze(1).to(DTYPE)
 71    
 72    F = forward_model_batched(points,transducers=board,p_ref=p_ref, transducer_radius=transducer_radius, k=k, norms=transducer_norms)
 73    Fx, Fy, Fz = forward_model_grad(points,transducers=board,p_ref=p_ref, transducer_radius=transducer_radius, k=k, transducer_norms=transducer_norms)
 74    Fxx, Fyy, Fzz = forward_model_second_derivative_unmixed(points,transducers=board,p_ref=p_ref, transducer_radius=transducer_radius, k=k,transducer_norms=transducer_norms)
 75    Fxy, Fxz, Fyz = forward_model_second_derivative_mixed(points,transducers=board,p_ref=p_ref, transducer_radius=transducer_radius, k=k,transducer_norms=transducer_norms)
 76
 77    p   = (F@activations)
 78    Px  = (Fx@activations)
 79    Py  = (Fy@activations)
 80    Pz  = (Fz@activations)
 81    Pxx = (Fxx@activations)
 82    Pyy = (Fyy@activations)
 83    Pzz = (Fzz@activations)
 84    Pxy = (Fxy@activations)
 85    Pxz = (Fxz@activations)
 86    Pyz = (Fyz@activations)
 87
 88
 89    # grad_p = torch.stack([Px,Py,Pz], dim=2).squeeze(3)
 90    # grad_px = torch.stack([Pxx,Pxy,Pxz])
 91    # grad_py = torch.stack([Pxy,Pyy,Pyz])
 92    # grad_pz = torch.stack([Pxz,Pyz,Pzz])
 93
 94    p_term_x= (p*Px.conj() + p.conj()*Px) 
 95    p_term_y= (p*Py.conj() + p.conj()*Py) 
 96    p_term_z= (p*Pz.conj() + p.conj()*Pz) 
 97
 98    # px_term = Px*grad_px.conj() + Px.conj()*grad_px
 99    # py_term = Py*grad_py.conj() + Py.conj()*grad_py
100    # pz_term = Pz*grad_pz.conj() + Pz.conj()*grad_pz
101
102    # K1 = V / (4*c.p_0*c.c_0**2)
103    # K2 = 3*V / (4*(2*c.f**2 * c.p_0))
104    K1, K2 = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density)
105
106    # grad_U = K1 * p_term - K2 * (px_term + py_term + pz_term)
107
108    grad_U_x = K1 * p_term_x - K2 * ((Pxx*Px.conj() + Px*Pxx.conj()) + (Pxy*Py.conj() + Py.conj()*Pxy) + (Pxz*Pz.conj() + Pxz.conj()*Pz))
109    grad_U_y = K1 * p_term_y - K2 * ((Pxy*Px.conj() + Px*Pxy.conj()) + (Pyy*Py.conj() + Py.conj()*Pyy) + (Pyz*Pz.conj() + Pyz.conj()*Pz))
110    grad_U_z = K1 * p_term_z - K2 * ((Pxz*Px.conj() + Px*Pxz.conj()) + (Pyz*Py.conj() + Py.conj()*Pyz) + (Pzz*Pz.conj() + Pxz.conj()*Pz))
111
112    grad_U = torch.stack([grad_U_x, grad_U_y, grad_U_z])
113
114    force = -(grad_U).real.squeeze(3).permute(1,2,0)
115    
116    if return_components:
117        return force[:,:,0], force[:,:,1], force[:,:,2] 
118    else:
119        return force 
120
121    
122def get_force_axis(activations:Tensor, points:Tensor,board:Tensor|None=None, axis:int=2, transducer_radius=c.radius, k=c.k, transducer_norms=None,
123                 medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p) -> Tensor:
124    '''
125    Returns the force in one axis on a particle using the analytical derivative of the Gor'kov potential and the piston model \n
126    Equivalent to `compute_force(activations, points,return_components=True)[axis]` \n 
127
128    :param activations: Transducer hologram
129    :param points: Points to propagate to
130    :param board: Transducers to use if `None` uses `acoustools.Utilities.TRANSDUCERS`
131    :param axis: Axis to take the force in
132    :return: force  
133    '''
134    if board is None:
135        board = TRANSDUCERS
136    forces = compute_force(activations, points,return_components=True, board=board, transducer_radius=transducer_radius, k=k, transducer_norms=transducer_norms,
137                           medium_density=medium_density, medium_speed=medium_speed,particle_density=particle_density, particle_speed=particle_speed)
138    force = forces[axis]
139
140    return force
141
142
143def force_mesh(activations:Tensor, points:Tensor, norms:Tensor, areas:Tensor, board:Tensor, grad_function:FunctionType=forward_model_grad, 
144               grad_function_args:dict={}, F_fun:FunctionType|None=forward_model_batched, F_function_args:dict={},
145               F:Tensor|None=None, Ax:Tensor|None=None, Ay:Tensor|None=None,Az:Tensor|None=None,
146               use_momentum:bool=False, return_components:bool=False, p_ref=c.P_ref, transducer_radius=c.radius, 
147               transducer_norms=None,k=c.k, 
148               medium_density = c.p_0, wave_speed = c.c_0, angular_frequency = c.angular_frequency ) -> Tensor:
149    '''
150    Returns the force on a mesh using a discritised version of Eq. 1 in `Acoustical boundary hologram for macroscopic rigid-body levitation`\n
151    :param activations: Transducer hologram
152    :param points: Points to propagate to
153    :param norms: The normals to the mesh faces
154    :param areas: The areas of the mesh points
155    :param board: Transducers to use 
156    :param grad_function: The function to use to compute the gradient of pressure
157    :param grad_function_args: The argument to pass to `grad_function`
158    :param F_fun: Function to compute F
159    :param F_function_args:Fucntion to compute Grad F
160    :param F: A precomputed forward propagation matrix, if `None` will be computed
161    :param Ax: The gradient of `F` wrt x, if `None` will be computed
162    :param Ay: The gradient of `F` wrt y, if `None` will be computed
163    :param Az: The gradient of `F` wrt z, if `None` will be computed
164    :param use_mpmentum: If true will add the term for momentum advection, for sound hard boundaries should be false
165    :param return_components: If True will return force, momentum_flux (force is still the total force)
166    :return: the force on each mesh element
167    '''
168
169    if F is None:
170        F = F_fun(points=points, transducers=board, p_ref=p_ref, transducer_radius=transducer_radius,norms=transducer_norms, k=k ,**F_function_args)
171    p = propagate(activations,points,board,A=F, p_ref=p_ref, transducer_radius=transducer_radius, k=k)
172    pressure_square = torch.abs(p)**2
173    pressure_time_average = 1/2 * pressure_square
174
175    # return pressure_time_average.expand((1,3,-1)), None 
176    if Ax is None or Ay is None or Az is None:
177        Ax, Ay, Az = grad_function(points=points, transducers=board, p_ref=p_ref, k=k, transducer_radius=transducer_radius, transducer_norms=transducer_norms, **grad_function_args)
178    
179    px = (Ax@activations).squeeze(2).unsqueeze(0)
180    py = (Ay@activations).squeeze(2).unsqueeze(0)
181    pz = (Az@activations).squeeze(2).unsqueeze(0)
182
183    grad  = torch.cat((px,py,pz),dim=1) 
184    velocity = grad /( 1j * medium_density * angular_frequency)
185
186    
187    k0 = 1/( medium_density * wave_speed**2)
188    velocity_time_average = 1/2 * torch.sum(velocity * velocity.conj().resolve_conj(), dim=1, keepdim=True).real 
189
190    # + velocity_time_average / velocity_time_average.max()
191    # pressure_square / pressure_square.max()
192
193    force = ( 0.5 * k0 * pressure_time_average - (medium_density / 2) * velocity_time_average) * norms
194
195    if use_momentum:        
196        momentum = medium_density/2 * (torch.sum(velocity * norms, dim=1, keepdim=True) * velocity.conj().resolve_conj()).real + 0j
197        
198        force += momentum 
199    else:
200        momentum = 0
201    
202    force *= -areas # *0.7
203    # force = torch.real(force) #Im(F) == 0 but needs to be complex till now for dtype compatability
204    # print(torch.sgn(torch.sgn(force) * torch.log(torch.abs(force))) == torch.sgn(force))
205
206    if return_components: 
207        return force, momentum
208    
209    return force
210
211def torque_mesh(activations:Tensor, points:Tensor, norms:Tensor, areas:Tensor, centre_of_mass:Tensor, board:Tensor,force:Tensor|None=None, 
212                grad_function:FunctionType=forward_model_grad,grad_function_args:dict={},F:Tensor|None=None, transducer_norms=None,
213                Ax:Tensor|None=None, Ay:Tensor|None=None,Az:Tensor|None=None, transducer_radius=c.radius, k=c.k, 
214               medium_density = c.p_0, wave_speed = c.c_0, angular_frequency = c.angular_frequency, p_ref=c.P_ref) -> Tensor:
215    '''
216    Returns the torque on a mesh using a discritised version of Eq. 1 in `Acoustical boundary hologram for macroscopic rigid-body levitation`\n
217    :param activations: Transducer hologram
218    :param points: Points to propagate to
219    :param norms: The normals to the mesh faces
220    :param areas: The areas of the mesh points
221    :param centre_of_mass: The position of the centre of mass of the mesh
222    :param board: Transducers to use 
223    :param force: Precomputed force on the mesh faces, if `None` will be computed
224    :param grad_function: The function to use to compute the gradient of pressure
225    :param grad_function_args: The argument to pass to `grad_function`
226    :param F: A precomputed forward propagation matrix, if `None` will be computed
227    :param Ax: The gradient of F wrt x, if `None` will be computed
228    :param Ay: The gradient of F wrt y, if `None` will be computed
229    :param Az: The gradient of F wrt z, if `None` will be computed
230    :return: the force on each mesh element
231    '''
232
233    if force is None:
234        force = force_mesh(activations, points, norms, areas, board,grad_function,grad_function_args,F=F, Ax=Ax, Ay=Ay, Az=Az, transducer_norms=transducer_norms,
235                           p_ref=p_ref, k=k, wave_speed=wave_speed, medium_density=medium_density, transducer_radius=transducer_radius, angular_frequency=angular_frequency)
236    force = force.to(DTYPE)
237    
238    displacement = points - centre_of_mass
239    displacement = displacement.to(DTYPE)
240    torque = torch.linalg.cross(displacement,force,dim=1)
241
242    return torch.real(torque)