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
  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,
 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, 
 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) -> 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    F = forward_model_batched(points,transducers=board,p_ref=p_ref, transducer_radius=transducer_radius, k=k)
 70    Fx, Fy, Fz = forward_model_grad(points,transducers=board,p_ref=p_ref, transducer_radius=transducer_radius, k=k)
 71    Fxx, Fyy, Fzz = forward_model_second_derivative_unmixed(points,transducers=board,p_ref=p_ref, transducer_radius=transducer_radius, k=k)
 72    Fxy, Fxz, Fyz = forward_model_second_derivative_mixed(points,transducers=board,p_ref=p_ref, transducer_radius=transducer_radius, k=k)
 73
 74    p   = (F@activations)
 75    Px  = (Fx@activations)
 76    Py  = (Fy@activations)
 77    Pz  = (Fz@activations)
 78    Pxx = (Fxx@activations)
 79    Pyy = (Fyy@activations)
 80    Pzz = (Fzz@activations)
 81    Pxy = (Fxy@activations)
 82    Pxz = (Fxz@activations)
 83    Pyz = (Fyz@activations)
 84
 85
 86    # grad_p = torch.stack([Px,Py,Pz], dim=2).squeeze(3)
 87    # grad_px = torch.stack([Pxx,Pxy,Pxz])
 88    # grad_py = torch.stack([Pxy,Pyy,Pyz])
 89    # grad_pz = torch.stack([Pxz,Pyz,Pzz])
 90
 91    p_term_x= (p*Px.conj() + p.conj()*Px) 
 92    p_term_y= (p*Py.conj() + p.conj()*Py) 
 93    p_term_z= (p*Pz.conj() + p.conj()*Pz) 
 94
 95    # px_term = Px*grad_px.conj() + Px.conj()*grad_px
 96    # py_term = Py*grad_py.conj() + Py.conj()*grad_py
 97    # pz_term = Pz*grad_pz.conj() + Pz.conj()*grad_pz
 98
 99    # K1 = V / (4*c.p_0*c.c_0**2)
100    # K2 = 3*V / (4*(2*c.f**2 * c.p_0))
101    K1, K2 = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density)
102
103    # grad_U = K1 * p_term - K2 * (px_term + py_term + pz_term)
104
105    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))
106    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))
107    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))
108
109    grad_U = torch.stack([grad_U_x, grad_U_y, grad_U_z])
110
111    force = -(grad_U).real.squeeze(3).permute(1,2,0)
112    
113    if return_components:
114        return force[:,:,0], force[:,:,1], force[:,:,2] 
115    else:
116        return force 
117
118    
119def get_force_axis(activations:Tensor, points:Tensor,board:Tensor|None=None, axis:int=2, transducer_radius=c.radius, k=c.k,
120                 medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p) -> Tensor:
121    '''
122    Returns the force in one axis on a particle using the analytical derivative of the Gor'kov potential and the piston model \n
123    Equivalent to `compute_force(activations, points,return_components=True)[axis]` \n 
124
125    :param activations: Transducer hologram
126    :param points: Points to propagate to
127    :param board: Transducers to use if `None` uses `acoustools.Utilities.TRANSDUCERS`
128    :param axis: Axis to take the force in
129    :return: force  
130    '''
131    if board is None:
132        board = TRANSDUCERS
133    forces = compute_force(activations, points,return_components=True, board=board, transducer_radius=transducer_radius, k=k, 
134                           medium_density=medium_density, medium_speed=medium_speed,particle_density=particle_density, particle_speed=particle_speed)
135    force = forces[axis]
136
137    return force
138
139
140def force_mesh(activations:Tensor, points:Tensor, norms:Tensor, areas:Tensor, board:Tensor, grad_function:FunctionType=forward_model_grad, 
141               grad_function_args:dict={}, F_fun:FunctionType|None=forward_model_batched, F_function_args:dict={},
142               F:Tensor|None=None, Ax:Tensor|None=None, Ay:Tensor|None=None,Az:Tensor|None=None,
143               use_momentum:bool=False, return_components:bool=False, p_ref=c.P_ref, transducer_radius=c.radius, k=c.k, 
144               medium_density = c.p_0, wave_speed = c.c_0, angular_frequency = c.angular_frequency ) -> Tensor:
145    '''
146    Returns the force on a mesh using a discritised version of Eq. 1 in `Acoustical boundary hologram for macroscopic rigid-body levitation`\n
147    :param activations: Transducer hologram
148    :param points: Points to propagate to
149    :param norms: The normals to the mesh faces
150    :param areas: The areas of the mesh points
151    :param board: Transducers to use 
152    :param grad_function: The function to use to compute the gradient of pressure
153    :param grad_function_args: The argument to pass to `grad_function`
154    :param F_fun: Function to compute F
155    :param F_function_args:Fucntion to compute Grad F
156    :param F: A precomputed forward propagation matrix, if `None` will be computed
157    :param Ax: The gradient of `F` wrt x, if `None` will be computed
158    :param Ay: The gradient of `F` wrt y, if `None` will be computed
159    :param Az: The gradient of `F` wrt z, if `None` will be computed
160    :param use_mpmentum: If true will add the term for momentum advection, for sound hard boundaries should be false
161    :param return_components: If True will return force, momentum_flux (force is still the total force)
162    :return: the force on each mesh element
163    '''
164
165    if F is None:
166        F = F_fun(points=points, transducers=board, p_ref=p_ref, transducer_radius=transducer_radius, k=k ,**F_function_args)
167    p = propagate(activations,points,board,A=F, p_ref=p_ref, transducer_radius=transducer_radius, k=k)
168    pressure_square = torch.abs(p)**2
169    pressure_time_average = 1/2 * pressure_square
170
171    # return pressure_time_average.expand((1,3,-1)), None 
172    if Ax is None or Ay is None or Az is None:
173        Ax, Ay, Az = grad_function(points=points, transducers=board, p_ref=p_ref, k=k, transducer_radius=transducer_radius, **grad_function_args)
174    
175    px = (Ax@activations).squeeze(2).unsqueeze(0)
176    py = (Ay@activations).squeeze(2).unsqueeze(0)
177    pz = (Az@activations).squeeze(2).unsqueeze(0)
178
179    grad  = torch.cat((px,py,pz),dim=1) 
180    velocity = grad /( 1j * medium_density * angular_frequency)
181
182    
183    k0 = 1/( medium_density * wave_speed**2)
184    velocity_time_average = 1/2 * torch.sum(velocity * velocity.conj().resolve_conj(), dim=1, keepdim=True).real 
185
186    # + velocity_time_average / velocity_time_average.max()
187    # pressure_square / pressure_square.max()
188
189    force = ( 0.5 * k0 * pressure_time_average - (medium_density / 2) * velocity_time_average) * norms
190
191    if use_momentum:        
192        momentum = medium_density/2 * (torch.sum(velocity * norms, dim=1, keepdim=True) * velocity.conj().resolve_conj()).real + 0j
193        
194        force += momentum 
195    else:
196        momentum = 0
197    
198    force *= -areas # *0.7
199    # force = torch.real(force) #Im(F) == 0 but needs to be complex till now for dtype compatability
200    # print(torch.sgn(torch.sgn(force) * torch.log(torch.abs(force))) == torch.sgn(force))
201
202    if return_components: 
203        return force, momentum
204    
205    return force
206
207def torque_mesh(activations:Tensor, points:Tensor, norms:Tensor, areas:Tensor, centre_of_mass:Tensor, board:Tensor,force:Tensor|None=None, 
208                grad_function:FunctionType=forward_model_grad,grad_function_args:dict={},F:Tensor|None=None, 
209                Ax:Tensor|None=None, Ay:Tensor|None=None,Az:Tensor|None=None, transducer_radius=c.radius, k=c.k, 
210               medium_density = c.p_0, wave_speed = c.c_0, angular_frequency = c.angular_frequency, p_ref=c.P_ref) -> Tensor:
211    '''
212    Returns the torque on a mesh using a discritised version of Eq. 1 in `Acoustical boundary hologram for macroscopic rigid-body levitation`\n
213    :param activations: Transducer hologram
214    :param points: Points to propagate to
215    :param norms: The normals to the mesh faces
216    :param areas: The areas of the mesh points
217    :param centre_of_mass: The position of the centre of mass of the mesh
218    :param board: Transducers to use 
219    :param force: Precomputed force on the mesh faces, if `None` will be computed
220    :param grad_function: The function to use to compute the gradient of pressure
221    :param grad_function_args: The argument to pass to `grad_function`
222    :param F: A precomputed forward propagation matrix, if `None` will be computed
223    :param Ax: The gradient of F wrt x, if `None` will be computed
224    :param Ay: The gradient of F wrt y, if `None` will be computed
225    :param Az: The gradient of F wrt z, if `None` will be computed
226    :return: the force on each mesh element
227    '''
228
229    if force is None:
230        force = force_mesh(activations, points, norms, areas, board,grad_function,grad_function_args,F=F, Ax=Ax, Ay=Ay, Az=Az, 
231                           p_ref=p_ref, k=k, wave_speed=wave_speed, medium_density=medium_density, transducer_radius=transducer_radius, angular_frequency=angular_frequency)
232    force = force.to(DTYPE)
233    
234    displacement = points - centre_of_mass
235    displacement = displacement.to(DTYPE)
236    torque = torch.linalg.cross(displacement,force,dim=1)
237
238    return torch.real(torque)