src.acoustools.Gorkov 

  1import torch
  2from acoustools.Utilities import device, propagate, forward_model_batched, forward_model_grad, TRANSDUCERS, DTYPE
  3import acoustools.Constants as c
  4
  5from torch import Tensor
  6from types import FunctionType
  7
  8
  9def gorkov(activations: Tensor, points: Tensor,board:Tensor|None=None, axis:str="XYZ", V:float=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius, transducer_norms=None,
 10                        medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p, **params) -> Tensor:
 11    '''
 12    Use to compute the gorkov potential at a point. Alias for `src.acoustools.Gorkov.gorkov_analytical`
 13    '''
 14    return gorkov_analytical(activations, points, board, axis, V=V, p_ref=p_ref, k=k, transducer_radius=transducer_radius, transducer_norms=transducer_norms,
 15                             medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density,particle_speed=particle_speed , **params)
 16
 17def gorkov_analytical(activations: Tensor, points: Tensor,board:Tensor|None=None, axis:str="XYZ", transducer_norms = None,
 18                        V:float=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius, 
 19                        medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p, angular_frequency=c.angular_frequency , **params) -> Tensor:
 20    '''
 21    Computes the Gorkov potential using analytical derivative of the piston model \n
 22    :param activation: The transducer activations to use 
 23    :param points: The points to compute the potential at
 24    :param board: The transducer boards to use
 25    :param axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis
 26    :param V: particle Volume
 27    :return: gorkov potential at each point
 28    ```Python
 29    from acoustools.Utilities import create_points, add_lev_sig
 30    from acoustools.Solvers import wgs
 31    from acoustools.Gorkov import gorkov_analytical
 32
 33    N=1
 34    B=1
 35    points = create_points(N,B)
 36    x = wgs(points)
 37    x = add_lev_sig(x)
 38    
 39    U_a = gorkov_analytical(x,points)
 40
 41    print("Analytical",U_a.data.squeeze())
 42    ```
 43    '''
 44
 45    if board is None:
 46        board = TRANSDUCERS
 47
 48    Fx, Fy, Fz = forward_model_grad(points, transducers=board, p_ref=p_ref, k=k, transducer_radius=transducer_radius, transducer_norms=transducer_norms)
 49    F = forward_model_batched(points,board, p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms=transducer_norms)
 50    
 51    p = torch.abs(F@activations)**2
 52
 53    px = (Fx@activations) if 'X' in axis else 0
 54    py = (Fy@activations) if 'Y' in axis else 0
 55    pz = (Fz@activations) if 'Z' in axis else 0
 56
 57    grad  = torch.cat((px,py,pz),dim=2)
 58
 59    K1, K2 = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density, angular_frequency=angular_frequency)
 60    g = (torch.sum(torch.abs(grad)**2, dim=2, keepdim=True))
 61    
 62    # K1 = 1/4*V*(1/(c.c_0**2*c.p_0) - 1/(c.c_p**2*c.p_p))
 63    # K2 = 3/4 * V * ((c.p_0 - c.p_p) / (c.angular_frequency**2 * c.p_0 * (c.p_0 * 2*c.p_p)))
 64    
 65    U = K1*p - K2*g
 66    return U
 67
 68
 69def get_gorkov_constants(V=c.V, p_0 = c.p_0, p_p=c.p_p, c_0=c.c_0, c_p=c.c_p, angular_frequency=c.angular_frequency ):
 70    '''
 71    Returns K1 and K2 for use in Gorkov computations, Uses the form shown in `Holographic acoustic elements for manipulation of levitated objects` \n
 72    :param: V: Particle Volume
 73    :param p_0: Density of medium
 74    :param p_p: Density of particle
 75    :param c_0: Speed of sound in medium
 76    :param c_p: speed of sound in particle
 77    :param angular_frequency: The angular frequency
 78    :returns K1, K2:
 791
 80    '''
 81    #Derived Bk.3 Pg.91
 82
 83    K1 = 1/4*V*(1/(c_0**2*p_0) - 1/(c_p**2*p_p)) 
 84    K2 = 3/4 * V * ((p_p - p_0) / (angular_frequency**2 * p_0 * (p_0 + 2*p_p))) 
 85
 86
 87    # exit()
 88    # K1 = V / (4*p_0*c_0**2)
 89    # K2 = 3*V / (4*(2*f**2 * p_0))
 90    # return 1, 1
 91    return K1, K2
 92
 93def gorkov_autograd(activations:Tensor, points:Tensor, K1:float|None=None, K2:float|None=None, V:float=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius,
 94                     medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p,
 95                     retain_graph:bool=False,board:Tensor|None=None, transducer_norms = None,**params) -> Tensor:
 96    '''
 97    Computes the Gorkov potential using pytorch's autograd system\n
 98    :param activation: The transducer activations to use 
 99    :param points: The points to compute the potential at. if `None` will use `acoustools.Utilities.TRANSDUCERS`
100    :param K1: The value for K1 in the Gorkov equation, if `None` will use `c.V / (4*c.p_0*c.c_0**2)`
101    :param K2: The value for K2 in the Gorkov equation, if `None` will use `3*c.V / (4*(2*c.f**2 * c.p_0))`
102    :param board: The transducer boards to use
103    :param retain_graph: Value will be passed to autograd
104    :return: gorkov potential at each point
105
106    ```Python
107    from acoustools.Utilities import create_points, add_lev_sig
108    from acoustools.Solvers import wgs
109    from acoustools.Gorkov import gorkov_autograd
110
111    N=1
112    B=1
113    points = create_points(N,B)
114    x = wgs(points)
115    x = add_lev_sig(x)
116    
117    U_ag = gorkov_autograd(x,points)
118
119    print("Autograd", U_ag.data.squeeze())
120    ```
121    '''
122
123    if board is None:
124        board = TRANSDUCERS
125
126    var_points = torch.autograd.Variable(points.data, requires_grad=True).to(device).to(DTYPE)
127
128    B = points.shape[0]
129    N = points.shape[2]
130    
131    if len(activations.shape) < 3:
132        activations.unsqueeze_(0)    
133    
134    pressure = propagate(activations.to(DTYPE),var_points,board=board, p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms=transducer_norms)
135    pressure.backward(torch.ones((B,N))+0j, inputs=var_points, retain_graph=retain_graph)
136    grad_pos = var_points.grad
137
138    if K1 is None or K2 is None:
139        K1_, K2_ = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density)
140        if K1 is None:
141            K1 = K1_
142        if K2 is None:
143            K2 = K2_
144
145    # K1 = 1/4*V*(1/(c.c_0**2*c.p_0) - 1/(c.c_p**2*c.p_p))
146    # K2 = 3/4 * V * ((c.p_0 - c.p_p) / (c.angular_frequency**2 * c.p_0 * (c.p_0 * 2*c.p_p)))
147
148
149    gorkov = K1 * torch.abs(pressure) **2 - K2 * torch.sum((torch.abs(grad_pos)**2),1)
150    return gorkov
151
152
153def gorkov_fin_diff(activations: Tensor, points:Tensor, axis:str="XYZ", stepsize:float = 0.000135156253,K1:float|None=None, K2:float|None=None,
154                    prop_function:FunctionType=propagate,prop_fun_args:dict={}, board:Tensor|None=None, transducer_norms=None,
155                    V=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius,
156                     medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p,) -> Tensor:
157    '''
158    Computes the Gorkov potential using finite differences to compute derivatives \n
159    :param activation: The transducer activations to use 
160    :param points: The points to compute the potential at
161    :param axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis
162    :param stepsize: The distance aroud points to add, default 0.000135156253
163    :param K1: The value for K1 in the Gorkov equation, if `None` will use `c.V / (4*c.p_0*c.c_0**2)`
164    :param K2: The value for K2 in the Gorkov equation, if `None` will use `3*c.V / (4*(2*c.f**2 * c.p_0))`
165    :param prop_function: Function to use to compute pressure
166    :param prop_fun_args: Arguments to pass to `prop_function`
167    :param board: The transducer boards to use if `None` use `acoustools.Utilities.TRANSDUCERS`
168    :param V: Particle Volume
169    :return: gorkov potential at each point
170
171    ```Python
172    from acoustools.Utilities import create_points, add_lev_sig
173    from acoustools.Solvers import wgs
174    from acoustools.Gorkov import gorkov_fin_diff
175
176    N=1
177    B=1
178    points = create_points(N,B)
179    x = wgs(points)
180    x = add_lev_sig(x)
181    
182    U_fd = gorkov_fin_diff(x,points)
183
184    print("Finite Differences",U_fd.data.squeeze())
185    ```
186    '''
187    # torch.autograd.set_detect_anomaly(True)
188    if board is None:
189        board = TRANSDUCERS
190    B = points.shape[0]
191    D = len(axis)
192    N = points.shape[2]
193
194    
195    if len(activations.shape) < 3:
196        activations = torch.unsqueeze(activations,0).clone().to(device)
197
198    fin_diff_points = get_finite_diff_points_all_axis(points, axis, stepsize)
199
200    pressure_points = prop_function(activations, fin_diff_points,board=board,p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms = transducer_norms,**prop_fun_args)
201    # if len(pressure_points.shape)>1:
202    # pressure_points = torch.squeeze(pressure_points,2)
203
204    pressure = pressure_points[:,:N]
205    pressure_fin_diff = pressure_points[:,N:]
206
207    split = torch.reshape(pressure_fin_diff,(B,2, ((2*D))*N // 2))
208    
209    grad = (split[:,0,:] - split[:,1,:]) / (2*stepsize)
210    
211    grad = torch.reshape(grad,(B,D,N))
212    grad_abs_square = torch.pow(torch.abs(grad),2)
213    grad_term = torch.sum(grad_abs_square,dim=1)
214
215    if K1 is None or K2 is None:
216        K1_, K2_ = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density)
217        if K1 is None:
218            K1 = K1_
219        if K2 is None:
220            K2 = K2_
221    
222    # p_in =  torch.abs(pressure)
223    p_in = torch.sqrt(torch.real(pressure) **2 + torch.imag(pressure)**2)
224    if len(p_in.shape) > 2:
225        p_in.squeeze_(2)
226    # p_in = torch.squeeze(p_in,2)
227
228    # K1 = 1/4*V*(1/(c.c_0**2*c.p_0) - 1/(c.c_p**2*c.p_p))
229    # K2 = 3/4 * V * ((c.p_0 - c.p_p) / (c.angular_frequency**2 * c.p_0 * (c.p_0 * 2*c.p_p)))
230
231    U = K1 * p_in**2 - K2 *grad_term
232    
233    return U
234
235
236def get_finite_diff_points(points:Tensor , axis:Tensor, stepsize:float = 0.000135156253) -> Tensor:
237    '''
238    Gets points for finite difference calculations in one axis\n
239    :param points: Points around which to find surrounding points
240    :param axis: The axis to add points in
241    :param stepsize: The distance aroud points to add, default 0.000135156253
242    :return: points 
243    '''
244    #points = Bx3x4
245    points_h = points.clone()
246    points_neg_h = points.clone()
247    points_h[:,axis,:] = points_h[:,axis,:] + stepsize
248    points_neg_h[:,axis,:] = points_neg_h[:,axis,:] - stepsize
249
250    return points_h, points_neg_h
251
252def get_finite_diff_points_all_axis(points: Tensor,axis: str="XYZ", stepsize:float = 0.000135156253) -> Tensor:
253    '''
254    Gets points for finite difference calculations\\
255    :param points: Points around which to find surrounding points\\
256    :param axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis\\
257    :param stepsize: The distance aroud points to add, default 0.000135156253\\
258    :return: Points
259    '''
260    B = points.shape[0]
261    D = len(axis)
262    N = points.shape[2]
263    fin_diff_points=  torch.zeros((B,3,((2*D)+1)*N)).to(device).to(DTYPE)
264    fin_diff_points[:,:,:N] = points.clone()
265
266    i = 2
267    if "X" in axis:
268        points_h, points_neg_h = get_finite_diff_points(points, 0, stepsize)
269        fin_diff_points[:,:,N:i*N] = points_h
270        fin_diff_points[:,:,D*N+(i-1)*N:D*N+i*N] = points_neg_h
271
272        i += 1
273
274    
275    if "Y" in axis:
276        points_h, points_neg_h = get_finite_diff_points(points, 1, stepsize)
277        fin_diff_points[:,:,(i-1)*N:i*N] = points_h
278        fin_diff_points[:,:,D*N+(i-1)*N:D*N+i*N] = points_neg_h
279        i += 1
280    
281    if "Z" in axis:
282        points_h, points_neg_h = get_finite_diff_points(points, 2, stepsize)
283        fin_diff_points[:,:,(i-1)*N:i*N] = points_h
284        fin_diff_points[:,:,D*N+(i-1)*N:D*N+i*N] = points_neg_h
285        i += 1
286    
287    return fin_diff_points
def gorkov( activations: torch.Tensor, points: torch.Tensor, board: torch.Tensor | None = None, axis: str = 'XYZ', V: float = 4.188790204666667e-09, p_ref=3.4000000000000004, k=732.7329804081634, transducer_radius=0.0045, transducer_norms=None, medium_density=1.2, medium_speed=343, particle_density=29.36, particle_speed=1052, **params) -> torch.Tensor: 
10def gorkov(activations: Tensor, points: Tensor,board:Tensor|None=None, axis:str="XYZ", V:float=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius, transducer_norms=None,
11                        medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p, **params) -> Tensor:
12    '''
13    Use to compute the gorkov potential at a point. Alias for `src.acoustools.Gorkov.gorkov_analytical`
14    '''
15    return gorkov_analytical(activations, points, board, axis, V=V, p_ref=p_ref, k=k, transducer_radius=transducer_radius, transducer_norms=transducer_norms,
16                             medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density,particle_speed=particle_speed , **params)

Use to compute the gorkov potential at a point. Alias for src.acoustools.Gorkov.gorkov_analytical

def gorkov_analytical( activations: torch.Tensor, points: torch.Tensor, board: torch.Tensor | None = None, axis: str = 'XYZ', transducer_norms=None, V: float = 4.188790204666667e-09, p_ref=3.4000000000000004, k=732.7329804081634, transducer_radius=0.0045, medium_density=1.2, medium_speed=343, particle_density=29.36, particle_speed=1052, angular_frequency=251327.41228000002, **params) -> torch.Tensor: 
18def gorkov_analytical(activations: Tensor, points: Tensor,board:Tensor|None=None, axis:str="XYZ", transducer_norms = None,
19                        V:float=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius, 
20                        medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p, angular_frequency=c.angular_frequency , **params) -> Tensor:
21    '''
22    Computes the Gorkov potential using analytical derivative of the piston model \n
23    :param activation: The transducer activations to use 
24    :param points: The points to compute the potential at
25    :param board: The transducer boards to use
26    :param axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis
27    :param V: particle Volume
28    :return: gorkov potential at each point
29    ```Python
30    from acoustools.Utilities import create_points, add_lev_sig
31    from acoustools.Solvers import wgs
32    from acoustools.Gorkov import gorkov_analytical
33
34    N=1
35    B=1
36    points = create_points(N,B)
37    x = wgs(points)
38    x = add_lev_sig(x)
39    
40    U_a = gorkov_analytical(x,points)
41
42    print("Analytical",U_a.data.squeeze())
43    ```
44    '''
45
46    if board is None:
47        board = TRANSDUCERS
48
49    Fx, Fy, Fz = forward_model_grad(points, transducers=board, p_ref=p_ref, k=k, transducer_radius=transducer_radius, transducer_norms=transducer_norms)
50    F = forward_model_batched(points,board, p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms=transducer_norms)
51    
52    p = torch.abs(F@activations)**2
53
54    px = (Fx@activations) if 'X' in axis else 0
55    py = (Fy@activations) if 'Y' in axis else 0
56    pz = (Fz@activations) if 'Z' in axis else 0
57
58    grad  = torch.cat((px,py,pz),dim=2)
59
60    K1, K2 = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density, angular_frequency=angular_frequency)
61    g = (torch.sum(torch.abs(grad)**2, dim=2, keepdim=True))
62    
63    # K1 = 1/4*V*(1/(c.c_0**2*c.p_0) - 1/(c.c_p**2*c.p_p))
64    # K2 = 3/4 * V * ((c.p_0 - c.p_p) / (c.angular_frequency**2 * c.p_0 * (c.p_0 * 2*c.p_p)))
65    
66    U = K1*p - K2*g
67    return U

Computes the Gorkov potential using analytical derivative of the piston model

Parameters
  • activation: The transducer activations to use
  • points: The points to compute the potential at
  • board: The transducer boards to use
  • axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis
  • V: particle Volume
Returns

gorkov potential at each point

from acoustools.Utilities import create_points, add_lev_sig
from acoustools.Solvers import wgs
from acoustools.Gorkov import gorkov_analytical

N=1
B=1
points = create_points(N,B)
x = wgs(points)
x = add_lev_sig(x)

U_a = gorkov_analytical(x,points)

print("Analytical",U_a.data.squeeze())
def get_gorkov_constants( V=4.188790204666667e-09, p_0=1.2, p_p=29.36, c_0=343, c_p=1052, angular_frequency=251327.41228000002): 
70def get_gorkov_constants(V=c.V, p_0 = c.p_0, p_p=c.p_p, c_0=c.c_0, c_p=c.c_p, angular_frequency=c.angular_frequency ):
71    '''
72    Returns K1 and K2 for use in Gorkov computations, Uses the form shown in `Holographic acoustic elements for manipulation of levitated objects` \n
73    :param: V: Particle Volume
74    :param p_0: Density of medium
75    :param p_p: Density of particle
76    :param c_0: Speed of sound in medium
77    :param c_p: speed of sound in particle
78    :param angular_frequency: The angular frequency
79    :returns K1, K2:
801
81    '''
82    #Derived Bk.3 Pg.91
83
84    K1 = 1/4*V*(1/(c_0**2*p_0) - 1/(c_p**2*p_p)) 
85    K2 = 3/4 * V * ((p_p - p_0) / (angular_frequency**2 * p_0 * (p_0 + 2*p_p))) 
86
87
88    # exit()
89    # K1 = V / (4*p_0*c_0**2)
90    # K2 = 3*V / (4*(2*f**2 * p_0))
91    # return 1, 1
92    return K1, K2

Returns K1 and K2 for use in Gorkov computations, Uses the form shown in Holographic acoustic elements for manipulation of levitated objects 

:param: V: Particle Volume
:param p_0: Density of medium
:param p_p: Density of particle
:param c_0: Speed of sound in medium
:param c_p: speed of sound in particle
:param angular_frequency: The angular frequency
:returns K1, K2:

1

def gorkov_autograd( activations: torch.Tensor, points: torch.Tensor, K1: float | None = None, K2: float | None = None, V: float = 4.188790204666667e-09, p_ref=3.4000000000000004, k=732.7329804081634, transducer_radius=0.0045, medium_density=1.2, medium_speed=343, particle_density=29.36, particle_speed=1052, retain_graph: bool = False, board: torch.Tensor | None = None, transducer_norms=None, **params) -> torch.Tensor: 
 94def gorkov_autograd(activations:Tensor, points:Tensor, K1:float|None=None, K2:float|None=None, V:float=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius,
 95                     medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p,
 96                     retain_graph:bool=False,board:Tensor|None=None, transducer_norms = None,**params) -> Tensor:
 97    '''
 98    Computes the Gorkov potential using pytorch's autograd system\n
 99    :param activation: The transducer activations to use 
100    :param points: The points to compute the potential at. if `None` will use `acoustools.Utilities.TRANSDUCERS`
101    :param K1: The value for K1 in the Gorkov equation, if `None` will use `c.V / (4*c.p_0*c.c_0**2)`
102    :param K2: The value for K2 in the Gorkov equation, if `None` will use `3*c.V / (4*(2*c.f**2 * c.p_0))`
103    :param board: The transducer boards to use
104    :param retain_graph: Value will be passed to autograd
105    :return: gorkov potential at each point
106
107    ```Python
108    from acoustools.Utilities import create_points, add_lev_sig
109    from acoustools.Solvers import wgs
110    from acoustools.Gorkov import gorkov_autograd
111
112    N=1
113    B=1
114    points = create_points(N,B)
115    x = wgs(points)
116    x = add_lev_sig(x)
117    
118    U_ag = gorkov_autograd(x,points)
119
120    print("Autograd", U_ag.data.squeeze())
121    ```
122    '''
123
124    if board is None:
125        board = TRANSDUCERS
126
127    var_points = torch.autograd.Variable(points.data, requires_grad=True).to(device).to(DTYPE)
128
129    B = points.shape[0]
130    N = points.shape[2]
131    
132    if len(activations.shape) < 3:
133        activations.unsqueeze_(0)    
134    
135    pressure = propagate(activations.to(DTYPE),var_points,board=board, p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms=transducer_norms)
136    pressure.backward(torch.ones((B,N))+0j, inputs=var_points, retain_graph=retain_graph)
137    grad_pos = var_points.grad
138
139    if K1 is None or K2 is None:
140        K1_, K2_ = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density)
141        if K1 is None:
142            K1 = K1_
143        if K2 is None:
144            K2 = K2_
145
146    # K1 = 1/4*V*(1/(c.c_0**2*c.p_0) - 1/(c.c_p**2*c.p_p))
147    # K2 = 3/4 * V * ((c.p_0 - c.p_p) / (c.angular_frequency**2 * c.p_0 * (c.p_0 * 2*c.p_p)))
148
149
150    gorkov = K1 * torch.abs(pressure) **2 - K2 * torch.sum((torch.abs(grad_pos)**2),1)
151    return gorkov

Computes the Gorkov potential using pytorch's autograd system

Parameters
  • activation: The transducer activations to use
  • points: The points to compute the potential at. if None will use acoustools.Utilities.TRANSDUCERS
  • K1: The value for K1 in the Gorkov equation, if None will use c.V / (4*c.p_0*c.c_0**2)
  • K2: The value for K2 in the Gorkov equation, if None will use 3*c.V / (4*(2*c.f**2 * c.p_0))
  • board: The transducer boards to use
  • retain_graph: Value will be passed to autograd
Returns

gorkov potential at each point

from acoustools.Utilities import create_points, add_lev_sig
from acoustools.Solvers import wgs
from acoustools.Gorkov import gorkov_autograd

N=1
B=1
points = create_points(N,B)
x = wgs(points)
x = add_lev_sig(x)

U_ag = gorkov_autograd(x,points)

print("Autograd", U_ag.data.squeeze())
def gorkov_fin_diff( activations: torch.Tensor, points: torch.Tensor, axis: str = 'XYZ', stepsize: float = 0.000135156253, K1: float | None = None, K2: float | None = None, prop_function: function = <function propagate>, prop_fun_args: dict = {}, board: torch.Tensor | None = None, transducer_norms=None, V=4.188790204666667e-09, p_ref=3.4000000000000004, k=732.7329804081634, transducer_radius=0.0045, medium_density=1.2, medium_speed=343, particle_density=29.36, particle_speed=1052) -> torch.Tensor: 
154def gorkov_fin_diff(activations: Tensor, points:Tensor, axis:str="XYZ", stepsize:float = 0.000135156253,K1:float|None=None, K2:float|None=None,
155                    prop_function:FunctionType=propagate,prop_fun_args:dict={}, board:Tensor|None=None, transducer_norms=None,
156                    V=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius,
157                     medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p,) -> Tensor:
158    '''
159    Computes the Gorkov potential using finite differences to compute derivatives \n
160    :param activation: The transducer activations to use 
161    :param points: The points to compute the potential at
162    :param axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis
163    :param stepsize: The distance aroud points to add, default 0.000135156253
164    :param K1: The value for K1 in the Gorkov equation, if `None` will use `c.V / (4*c.p_0*c.c_0**2)`
165    :param K2: The value for K2 in the Gorkov equation, if `None` will use `3*c.V / (4*(2*c.f**2 * c.p_0))`
166    :param prop_function: Function to use to compute pressure
167    :param prop_fun_args: Arguments to pass to `prop_function`
168    :param board: The transducer boards to use if `None` use `acoustools.Utilities.TRANSDUCERS`
169    :param V: Particle Volume
170    :return: gorkov potential at each point
171
172    ```Python
173    from acoustools.Utilities import create_points, add_lev_sig
174    from acoustools.Solvers import wgs
175    from acoustools.Gorkov import gorkov_fin_diff
176
177    N=1
178    B=1
179    points = create_points(N,B)
180    x = wgs(points)
181    x = add_lev_sig(x)
182    
183    U_fd = gorkov_fin_diff(x,points)
184
185    print("Finite Differences",U_fd.data.squeeze())
186    ```
187    '''
188    # torch.autograd.set_detect_anomaly(True)
189    if board is None:
190        board = TRANSDUCERS
191    B = points.shape[0]
192    D = len(axis)
193    N = points.shape[2]
194
195    
196    if len(activations.shape) < 3:
197        activations = torch.unsqueeze(activations,0).clone().to(device)
198
199    fin_diff_points = get_finite_diff_points_all_axis(points, axis, stepsize)
200
201    pressure_points = prop_function(activations, fin_diff_points,board=board,p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms = transducer_norms,**prop_fun_args)
202    # if len(pressure_points.shape)>1:
203    # pressure_points = torch.squeeze(pressure_points,2)
204
205    pressure = pressure_points[:,:N]
206    pressure_fin_diff = pressure_points[:,N:]
207
208    split = torch.reshape(pressure_fin_diff,(B,2, ((2*D))*N // 2))
209    
210    grad = (split[:,0,:] - split[:,1,:]) / (2*stepsize)
211    
212    grad = torch.reshape(grad,(B,D,N))
213    grad_abs_square = torch.pow(torch.abs(grad),2)
214    grad_term = torch.sum(grad_abs_square,dim=1)
215
216    if K1 is None or K2 is None:
217        K1_, K2_ = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density)
218        if K1 is None:
219            K1 = K1_
220        if K2 is None:
221            K2 = K2_
222    
223    # p_in =  torch.abs(pressure)
224    p_in = torch.sqrt(torch.real(pressure) **2 + torch.imag(pressure)**2)
225    if len(p_in.shape) > 2:
226        p_in.squeeze_(2)
227    # p_in = torch.squeeze(p_in,2)
228
229    # K1 = 1/4*V*(1/(c.c_0**2*c.p_0) - 1/(c.c_p**2*c.p_p))
230    # K2 = 3/4 * V * ((c.p_0 - c.p_p) / (c.angular_frequency**2 * c.p_0 * (c.p_0 * 2*c.p_p)))
231
232    U = K1 * p_in**2 - K2 *grad_term
233    
234    return U

Computes the Gorkov potential using finite differences to compute derivatives

Parameters
  • activation: The transducer activations to use
  • points: The points to compute the potential at
  • axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis
  • stepsize: The distance aroud points to add, default 0.000135156253
  • K1: The value for K1 in the Gorkov equation, if None will use c.V / (4*c.p_0*c.c_0**2)
  • K2: The value for K2 in the Gorkov equation, if None will use 3*c.V / (4*(2*c.f**2 * c.p_0))
  • prop_function: Function to use to compute pressure
  • prop_fun_args: Arguments to pass to prop_function
  • board: The transducer boards to use if None use acoustools.Utilities.TRANSDUCERS
  • V: Particle Volume
Returns

gorkov potential at each point

from acoustools.Utilities import create_points, add_lev_sig
from acoustools.Solvers import wgs
from acoustools.Gorkov import gorkov_fin_diff

N=1
B=1
points = create_points(N,B)
x = wgs(points)
x = add_lev_sig(x)

U_fd = gorkov_fin_diff(x,points)

print("Finite Differences",U_fd.data.squeeze())
def get_finite_diff_points( points: torch.Tensor, axis: torch.Tensor, stepsize: float = 0.000135156253) -> torch.Tensor: 
237def get_finite_diff_points(points:Tensor , axis:Tensor, stepsize:float = 0.000135156253) -> Tensor:
238    '''
239    Gets points for finite difference calculations in one axis\n
240    :param points: Points around which to find surrounding points
241    :param axis: The axis to add points in
242    :param stepsize: The distance aroud points to add, default 0.000135156253
243    :return: points 
244    '''
245    #points = Bx3x4
246    points_h = points.clone()
247    points_neg_h = points.clone()
248    points_h[:,axis,:] = points_h[:,axis,:] + stepsize
249    points_neg_h[:,axis,:] = points_neg_h[:,axis,:] - stepsize
250
251    return points_h, points_neg_h

Gets points for finite difference calculations in one axis

Parameters
  • points: Points around which to find surrounding points
  • axis: The axis to add points in
  • stepsize: The distance aroud points to add, default 0.000135156253
Returns

points

def get_finite_diff_points_all_axis( points: torch.Tensor, axis: str = 'XYZ', stepsize: float = 0.000135156253) -> torch.Tensor: 
253def get_finite_diff_points_all_axis(points: Tensor,axis: str="XYZ", stepsize:float = 0.000135156253) -> Tensor:
254    '''
255    Gets points for finite difference calculations\\
256    :param points: Points around which to find surrounding points\\
257    :param axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis\\
258    :param stepsize: The distance aroud points to add, default 0.000135156253\\
259    :return: Points
260    '''
261    B = points.shape[0]
262    D = len(axis)
263    N = points.shape[2]
264    fin_diff_points=  torch.zeros((B,3,((2*D)+1)*N)).to(device).to(DTYPE)
265    fin_diff_points[:,:,:N] = points.clone()
266
267    i = 2
268    if "X" in axis:
269        points_h, points_neg_h = get_finite_diff_points(points, 0, stepsize)
270        fin_diff_points[:,:,N:i*N] = points_h
271        fin_diff_points[:,:,D*N+(i-1)*N:D*N+i*N] = points_neg_h
272
273        i += 1
274
275    
276    if "Y" in axis:
277        points_h, points_neg_h = get_finite_diff_points(points, 1, stepsize)
278        fin_diff_points[:,:,(i-1)*N:i*N] = points_h
279        fin_diff_points[:,:,D*N+(i-1)*N:D*N+i*N] = points_neg_h
280        i += 1
281    
282    if "Z" in axis:
283        points_h, points_neg_h = get_finite_diff_points(points, 2, stepsize)
284        fin_diff_points[:,:,(i-1)*N:i*N] = points_h
285        fin_diff_points[:,:,D*N+(i-1)*N:D*N+i*N] = points_neg_h
286        i += 1
287    
288    return fin_diff_points

Gets points for finite difference calculations\

Parameters
  • points: Points around which to find surrounding points\
  • axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis\
  • stepsize: The distance aroud points to add, default 0.000135156253\
Returns

Points