Acoustic simulation - phase information needed (using Python)

[Annalisa]

Hi!
How can I extract the phase value of the pressure p(f) on a sensor point directly by Python?
I don' t know if that's possible ( to directly extract the phase).

Alternatively, I can extract the real and imaginary part of the complex pressure on my n sensors and then derive the phase values. BUT I can do it, extracting in separate excel files each pressure value. Not easy to handle.
I would need to allocate the n complex pressure values in a single array and then, calculate the phase from these values. I am moving my first steps in Python, so I really appreciate some helps.

Thanks

[montanaro]

Hi,

For a point sensor p(f):

Under properties you get the option for Complex component and can select Phase (deg).

from Python:

import numpy
import s4l_v1.analysis as analysis
import s4l_v1.document as document
import s4l_v1.model as model
import s4l_v1.units as units
from s4l_v1 import ReleaseVersion
from s4l_v1 import Unit

# Adding a new SimulationExtractor
simulation = document.AllSimulations["Single Element Focused Transducer"]
simulation_extractor = simulation.Results()

# Adding a new AcousticSensorExtractor
acoustic_sensor_extractor = simulation_extractor["Focus  (SEFT)"]

# Get Point Sensor 
input = acoustic_sensor_extractor.Outputs["p(f)"]
# Get Point Sensor data (output is complex pressure)
output = input.Data.GetComponent(0)

# Post Process: Angle
output_angle = numpy.angle(output)

[montanaro]

For Field sensor p(x,y,z,f):

Easiest is to create a 'Calculator' and calculate it yourself, then right click field and select 'To Python..' to get the corresponding Python code

Select overall field -> p(x,y,z,f) -> Field Data Tools -> Calculator -> Value Location: Cell Center -> Scalar Variables: F (Re{p}) -> Scalar Variables: G (Im{p}) -> Expression: atan(G/F) (change it to read this or whatever suitable formula you wish)

Then you can use the viewers (note, careful with division by zero, you'll need to adjust the color bar accordingly, also you'll get 'wrapped' angles)

(In the next reply I'll show you how you can extract the field, manipulate it, and plot it back ... Very useful to plot arbitrary numpy arrays in S4L)

Python code:

import numpy
import s4l_v1.analysis as analysis
import s4l_v1.document as document
import s4l_v1.model as model
import s4l_v1.units as units
from s4l_v1 import ReleaseVersion
from s4l_v1 import Unit

# Creating the analysis pipeline
# Adding a new SimulationExtractor
simulation = document.AllSimulations["Single Element Focused Transducer"]
simulation_extractor = simulation.Results()

# Adding a new AcousticSensorExtractor
acoustic_sensor_extractor = simulation_extractor["Overall Field"]

# Adding a new FieldCalculator
inputs = [acoustic_sensor_extractor.Outputs["p(x,y,z,f)"]]
field_calculator = analysis.field.FieldCalculator(inputs=inputs)
field_calculator.Expression = u"atan(G/F)"
field_calculator.ValueLocation = field_calculator.ValueLocation.enum.CellCenter
field_calculator.UpdateAttributes()
document.AllAlgorithms.Add(field_calculator)

# Adding a new SliceFieldViewer
inputs = [field_calculator.Outputs["Result(x,y,z)"]]
slice_field_viewer = analysis.viewers.SliceFieldViewer(inputs=inputs)
slice_field_viewer.Data.Mode = slice_field_viewer.Data.Mode.enum.QuantityRealPart
slice_field_viewer.Slice.Plane = slice_field_viewer.Slice.Plane.enum.XZ
slice_field_viewer.Slice.Index = 155
slice_field_viewer.UpdateAttributes()
document.AllAlgorithms.Add(slice_field_viewer)

[montanaro]

To get the field as a numpy array:

import numpy
import s4l_v1.analysis as analysis
import s4l_v1.document as document
import s4l_v1.model as model
import s4l_v1.units as units
from s4l_v1 import ReleaseVersion
from s4l_v1 import Unit

# Creating the analysis pipeline
# Adding a new SimulationExtractor
simulation = document.AllSimulations["Single Element Focused Transducer"]
simulation_extractor = simulation.Results()

# Adding a new AcousticSensorExtractor
acoustic_sensor_extractor = simulation_extractor["Overall Field"]
document.AllAlgorithms.Add(acoustic_sensor_extractor)

# Adding a new SliceFieldViewer
field = acoustic_sensor_extractor.Outputs["p(x,y,z,f)"]
# Remember to update
# Otherwise this script will fail unless you have already extracted results in the GUI
field.Update()

field_np = field.Data.Field(0) # Field as flattened 1D array
# 3D Field dimensions (which might be different to grid dimensions)
# Data is often times defined at cell centers, Grid at cell nodes
dims = numpy.array(field.Data.Grid.Dimensions)
if numpy.prod(dims) != len(field_np):
	dims = dims-1
field_np = field_np.reshape(dims, order='F') # Field as 3D array

xaxis = field.Data.Grid.XAxis
yaxis = field.Data.Grid.YAxis
zaxis = field.Data.Grid.ZAxis

# You can now manipulate field_np as you wish
field_angle = numpy.unwrap(numpy.angle(field_np))

To plot your numpy array back into sim4life you create what is known as a 'Trivial Producer' which holds your numpy array which can be manipulated in the analysis tab:

data = field_angle

grid = analysis.core.RectilinearGrid()
grid.XAxis = xaxis
grid.YAxis = yaxis
grid.ZAxis = zaxis

# Careful, the new field is a field of Doubles, changes to ComplexFieldData
# if you need to work with complex numbers
new_field = analysis.core.DoubleFieldData()  
new_field.Grid = grid
if data.size == (len(grid.XAxis) * len(grid.YAxis) * len(grid.ZAxis)):
	new_field.ValueLocation = analysis.core.eValueLocation.kNode
elif data.size == (len(grid.XAxis)-1) * (len(grid.YAxis)-1) * (len(grid.ZAxis)-1):
	new_field.ValueLocation = analysis.core.eValueLocation.kCellCenter
else:
	print "ERROR: Grid and Data don't match"
	assert False
new_field.NumberOfComponents = 1
new_field.NumberOfSnapshots = 1
new_field.Quantity.Name = "My Field"

# Note: memory layout is such that x is fastest, z slowest dimension
values = data.ravel('F')	
values = values.astype(numpy.float64) # Change if working with ComplexFieldData
new_field.SetField( 0, values )
assert new_field.Check()

producer = analysis.core.TrivialProducer()
producer.Description = "My Trivial Producer"
producer.SetDataObject(new_field)

document.AllAlgorithms.Add(producer)

[Annalisa]

Thank you very much for the detailed instructions.
I tried the first code you gave me and it works for a single sensor as expected.
(just need to add input.Update()).
In my case I have several sensors so I need to use a for loop in order to extract the phase value on each sensor.
My code is:

from future import absolute_import
from future import print_function
import sys, os, math
import h5py
import s4l_v1.model as model
import s4l_v1.simulation.acoustic as acoustic
import s4l_v1.analysis as analysis
import s4l_v1.document as document
import scipy.io

import numpy
import s4l_v1.document as document
import s4l_v1.materials.database as database
import s4l_v1.model as model
import s4l_v1.simulation.acoustic as acoustic
import s4l_v1.units as units
from s4l_v1 import ReleaseVersion
from s4l_v1 import Unit

try:
# Define the version to use for default values
ReleaseVersion.set_active(ReleaseVersion.version5_2)

# Creating the analysis pipeline
# Adding a new SimulationExtractor
simulation = document.AllSimulations["My simulation"]
simulation_extractor = simulation.Results()


for i in range(0,128):

	# Adding a new AcousticSensorExtractor
	acoustic_sensor_extractor = simulation_extractor["Sensor" + str(i)]
	document.AllAlgorithms.Add(acoustic_sensor_extractor)
	input = acoustic_sensor_extractor.Outputs["p(f)"]
	input.Update()

	# Get Point Sensor data (output is complex pressure)
	output = input.Data.GetComponent(i)
	
	# Post Process: Angle
	output_angle = numpy.angle(output)

except Exception as exc:
import traceback
traceback.print_exc(exc)
Reset active version to default
ReleaseVersion.reset()
raise(exc)

When I run this code the console gives this message:
File "C:\Users\PC\Documents\Sim4Life\5.2\Scripts\temp\BP p(f) extractor.py", line 58
Reset active version to default
^
SyntaxError: invalid syntax

any idea?

Thanks

Annalisa

[montanaro]

"input.Update()"

• you're right. As a general rule, make sure to always have this (if results are extracted and updated in the GUI you might not get this error)

Not sure why you get this error but this line seems wrong:
GetComponent(i) ... Should probably be GetComponent(0)

Components are typically useful when you have vector fields, In which case GetComponent(0) probably corresponds to E_x, GetComponent(1) to E_y, etc.. They should typically be set to 0 since most fields have only 1 'component'. (not to be confused with snapshots which typically refer to time snapshots or frequency snapshots, where 0 is the main harmonic)

[Annalisa]

GetComponent(0), yes you are right. That was a typo.
The problem still persists.
How can I extract the phase on each sensors using a "for" loop?

thanks again

[montanaro]

Just tried it on a simple example and seemed to work ... Maybe get rid of

try:
# Define the version to use for default values
ReleaseVersion.set_active(ReleaseVersion.version5_2)

and

except Exception as exc:
import traceback
traceback.print_exc(exc)
Reset active version to default
ReleaseVersion.reset()
raise(exc)

(seems that 'reset active version ...' line should be a comment)

This is what I tried:

# -*- coding: utf-8 -*-
# Creating the analysis pipeline
# Adding a new SimulationExtractor
simulation = document.AllSimulations["Single Element Focused Transducer - Copy"]
simulation_extractor = simulation.Results()

for i in range(1,5):

	# Adding a new AcousticSensorExtractor
	acoustic_sensor_extractor = simulation_extractor["Focus " + str(i)]
	document.AllAlgorithms.Add(acoustic_sensor_extractor)
	input = acoustic_sensor_extractor.Outputs["p(f)"]
	input.Update()

	# Get Point Sensor data (output is complex pressure)
	output = input.Data.GetComponent(0)
	
	# Post Process: Angle
	output_angle = numpy.angle(output)
	print i, output_angle

and this was my output:

1 [ 1.69136047]
2 [ 1.69136047]
3 [ 1.69136047]
4 [ 1.69136047]

(my point sensors were all at the same location)

[Annalisa]

This post is deleted!

[Annalisa]

Now, I want assign each phase value to an element of my transducer.
But using my code:

for i in range(0,128):
	
	source_settings = acoustic.SourceSettings()
	source_settings.Name = "Source Settings" + str(i)
	source_settings.Amplitude = 100000.0, Unit("Pa")		
	source_settings.Phase = -output_angle[i], units.Degrees
	simulation.Add(source_settings, components[i])

does not work because in output_angle (evaluated with the previous code) I have only the last value (corresponding to sensor#127).
I would need to allocate all the phase values in a callable list (?) in order to be able to use the for loop.
Thank for further helps

[montanaro]

You need to make output_angle into an array and append the new values at the end every time (the i-th element will have the i-th phase you want)

output_angle = numpy.array([])
for i in range(0,128):
	# your code
	output_angle = numpy.append(output_angle, numpy.angle(output))