add speedup test, initial transfer to gpu

README.md

@@ -7,7 +7,10 @@
 <!-- [![Publish to PyPI](https://github.com/jakep72/acoustipy/actions/workflows/python-publish.yml/badge.svg)](https://github.com/jakep72/acoustipy/actions/workflows/python-publish.yml) -->
 <!-- [![mkdocs-deploy](https://github.com/jakep72/acoustipy/actions/workflows/build_docs.yml/badge.svg)](https://github.com/jakep72/acoustipy/actions/workflows/build_docs.yml) -->
 
-# This branch is under active development and some functionality is currently broken.  The intent is to transition to pytorch from numpy to potentially utilize GPU processing in the transfer matrix operations and gain access to the gradients for parameter identification
+# This branch is under active development and some functionality is currently broken.  
+### The intent is to transition to pytorch from numpy to potentially utilize GPU processing in the transfer matrix operations and gain access to the gradients for parameter identification.  Ideally, no breaking changes to the current API will be introduced.
+### Preliminary testing indicates ~8x average speedup for transfer matrix calculations just from the transition from numpy to pytorch with no GPU acceleration.
+### 
 
 # acoustipy
 

speedup_test.py

@@ -0,0 +1,36 @@
+from src.acoustipy.TMM import AcousticTMM
+import time
+
+times = []
+# Create an AcousticTMM object, specifying a diffuse sound field at 20C
+for i in range(100):
+    s = time.time()
+    structure = AcousticTMM(incidence='Diffuse',air_temperature=20, device='cpu')
+
+    # Define the layers of the material using various models
+    layer1 = structure.Add_Resistive_Screen(thickness=1,flow_resistivity=100000,porosity=.86)
+    layer2 = structure.Add_DBM_Layer(thickness = 25.4,flow_resistivity=60000)
+    layer3 = structure.Add_Resistive_Screen(thickness = 1, flow_resistivity=500000,porosity=.75)
+
+    # Specify the material backing condition -- in this case a 400mm air gap
+    air = structure.Add_Air_Layer(thickness = 400)
+
+    # Build the total transfer matrix of the structure + air gap
+    transfer_matrix = structure.assemble_structure(layer1,layer2,layer3,air)
+
+    # Calculate the frequency dependent narrow band absorption coefficients
+    absorption = structure.absorption(transfer_matrix)
+
+    # Calculate the 3rd octave bands absorption coefficients
+    bands = structure.octave_bands(absorption)
+
+    # Calculate the four frequency average absorption
+    FFA = structure.FFA(bands)
+    t = time.time()-s
+    print(t)
+    times.append(t)
+
+print("Average Time: ", sum(times)/len(times))
+
+# Plot and display the narrow and 3rd band coefficients on the same figure
+# structure.plot_curve([absorption,bands],["absorption","third octave"])

src/acoustipy/TMM.py

@@ -4,6 +4,7 @@
 import pandas as pd
 import os
 import matplotlib.pyplot as plt
+import torch.backends.opt_einsum
 from acoustipy.Database import AcoustiBase
 
 
@@ -74,7 +75,8 @@ def __init__(self,
                  Cv:float = 0.717425,
                  viscosity:float = 1.825e-05,
                  Pr:float = .7157,
-                 P0:float = 101325
+                 P0:float = 101325,
+                 device:str = 'cuda:0'
                  ):
 
         #Define limits on the frequency range and angles
@@ -114,6 +116,7 @@ def __init__(self,
         self.viscosity = viscosity
         self.Pr = Pr
         self.P0 = P0
+        self.device = device
         self._custom_freq =  torch.arange(self.fmin,self.fmax+self.fs,self.fs)
         self.layers = []
 
@@ -220,8 +223,10 @@ def _create_layer_TM(self,
             Diffuse field --> 2 x 2 x len(frequency) x len(angles) array representing the transfer matrix
         
         """
+        kp = kp.to(self.device)
+        Zp = Zp.to(self.device)
         if self.incidence == "Normal":
-            TM = torch.zeros((2,2,len(self.frequency)),dtype = torch.complex64)
+            TM = torch.zeros((2,2,len(self.frequency)),dtype = torch.complex64).to(self.device)
             TM[0][0] = torch.cos(kp*thickness)
             TM[0][1] = 1j*Zp*torch.sin(kp*thickness)
             TM[1][0] = (1j/Zp)*torch.sin(kp*thickness)
@@ -231,12 +236,17 @@ def _create_layer_TM(self,
 
         elif self.incidence == "Diffuse":
             angles = torch.arange(self.angles[0],self.angles[1],self.angles[2])
-            vs = torch.sin(torch.deg2rad(angles))
-            TM = torch.zeros((2,2,len(self.frequency),len(angles)),dtype = torch.complex64)
-            kpx = torch.zeros((len(self.frequency),len(angles)),dtype = torch.complex64)
+            vs = torch.sin(torch.deg2rad(angles)).to(self.device)
+            TM = torch.zeros((2,2,len(self.frequency),len(angles)),dtype = torch.complex64).to(self.device)
+            kpx = torch.zeros((len(self.frequency),len(angles)),dtype = torch.complex64).to(self.device)
 
-            Kp = torch.einsum('ij,ij -> ij',torch.tile(kp[:,None], (1,len(angles))),torch.tile(kp[:,None], (1,len(angles))))
-            K0 = torch.einsum('ij,ij -> ij',torch.tile(self.k0[:,None],(1,len(angles))),torch.tile(self.k0[:,None],(1,len(angles))))
+            v1 = torch.tile(kp[:,None], (1,len(angles))).to(self.device)
+            v2 = torch.tile(kp[:,None], (1,len(angles))).to(self.device)
+            v3 = torch.tile(self.k0[:,None],(1,len(angles))).to(self.device)
+            v4 = torch.tile(self.k0[:,None],(1,len(angles))).to(self.device)
+
+            Kp = torch.einsum('ij,ij -> ij',v1,v2)
+            K0 = torch.einsum('ij,ij -> ij',v3,v4)
             VS = torch.einsum('ij,ij -> ij',vs[None,:],vs[None,:])
             kpx[:,:] = torch.sqrt((Kp)-(torch.einsum('ij,ij -> ij',K0,VS)))
 
@@ -1603,7 +1613,7 @@ def absorption(self,
             R = (Zst-self.Z0)/(Zst+self.Z0)
 
             A = 1-abs(R)**2
-            curve = torch.column_stack((self.frequency,A))
+            curve = torch.column_stack((self.frequency,A.to('cpu')))
 
             return (curve)
 
@@ -1612,7 +1622,7 @@ def absorption(self,
             angles = torch.arange(self.angles[0],self.angles[1],self.angles[2])
             v = torch.cos(torch.deg2rad(angles))
 
-            Zst = Tt[0][0][:][:]/ Tt[1][0][:][:]
+            Zst = (Tt[0][0][:][:]/ Tt[1][0][:][:]).to('cpu')
             R = ((Zst*v)-self.Z0)/((Zst*v)+self.Z0)
 
             a = 1-torch.abs(R)**2