Handling Large Field Data Results

[dbsim4]

Hello,

I'm working on Acoustic simulations, which require a small grid step and thus very large grid size.

The results file is ~40GB, probably due to the "Overall Field" sensor. I'm working on a laptop with 32GB RAM, and it takes a very long time to open/view the field data (Slice Viewer), and frequently crashes. The point sensors are not a problem.

I only want the animations of the pressure field in 2 planes (centered on the transducer, in 2 different orientations). Is there any way to only store this data, and avoid a huge results file?

Or a better Analysis workflow that avoids loading the whole file into RAM? I've tried Crop and Extract Slice, but it still takes an extremely long time to load.

I'd be fine with a Python script, or something that just extracts the 2 necessary slices over time & saves them. Or if there's a way to only record snapshots of the 2 slices instead of the entire 3D field...

Secondary question - we have an HPC cluster that we've been thinking of migrating the Ares solver to, but I'd probably still be working with the modeling & analysis remotely on my laptop. Will handling large results files still be an issue even with an HPC backend (or is there a recommended workflow)?

[LJ]

I don't think I have the exact answer to your question, but this might be useful.
I find that running my simulation using the python script helps avoid issues such as the GUI crashing (I assume because the entire field doesn't actually have to be loaded for me to use parts of it? Not sure). I tend to automate the process to extract the information I need (using masks or slices) in a loop that saves to an hdf5 file at each step. This helps with memory issues, and I can easily and quickly access just the information I need in the analysis step.

[dbsim4]

Thank you for the suggestion! I'd like to learn more about the Python interface, if you know any resources. (I'm very comfortable in Python but very new to Sim4Life). Do you run the scripts from the Python console in the GUI, or in a separate shell (bash/powershell/etc.)?

I'm having some luck with HDF Viewer and Python h5py library as suggested in this answer: https://forum.zmt.swiss/topic/265/s4l-is-unable-to-read-and-load-my-project/9

It seems like it can parse the metadata of the file nicely without loading everything into RAM all at once, and shows a nice tree/hierarchy of objects.

I've found the Overall Field p(x,y,z,t) snapshots I'm looking for, now I need to find my transducer center & point sensor coordinates with respect to grid indexes (I think?) in order to pull the correct slices (centered on the transducer).

[montanaro]

You can also create a wireframe around your volume of interest (Don't know if this works for a 2D slice / wireframe) and add that to your Field Sensors - then in PostPro you can extract only that. I don't know if you can disable recording from the Overall Field Sensor, but if you're lucky you can either delete the Overall Field Sensor or maybe there's an option to disable it in the Field Sensor part.

As for Python, the easiest way to get started is to create your simulation and then right click on the simulation and select 'To Python..' this automatically generates a script to recreate the entire simulation via Python.

Additionally you can also look here for more info, but I'd start with the 'To Python..' option and then maybe expand it with what you see in this script:

"C:\Users\Public\Documents\Sim4Life\6.2\Documentation\Tutorials\Tutorial Scripts\tutorial_acoustic_seft.py"

[montanaro]

Ah, you probably also want scripts to deal with post processing and data extraction ... The AnalyzeSimulation() in the tutorial script shows you how to do that

Here's a short script to go through all simulations with results and then take the absolute value of pressure and normalize by the max:

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

for sim in document.AllSimulations:
	# Creating the analysis pipeline
	# Adding a new SimulationExtractor
	if sim.HasResults() == False:
		continue
		
	simulation_extractor = sim.Results()

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

	# Get dimensions
	acoustic_sensor_extractor.Outputs["p(x,y,z,f)"].Update()
	
	dims = acoustic_sensor_extractor.Outputs["p(x,y,z,f)"].Data.Grid.Dimensions
	max = abs(acoustic_sensor_extractor.Outputs["p(x,y,z,f)"].Data.Field(0)).max()

	# 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"sqrt(F^2+G^2)/" + str(max)
	field_calculator.ValueLocation = field_calculator.ValueLocation.enum.CellCenter
	field_calculator.UpdateAttributes()
	document.AllAlgorithms.Add(field_calculator)

[dbsim4]

Thank you, this is all very helpful!

I'm getting comfortable with HDF and h5py; the one thing I'm still getting hung up on is finding the indexes/coordinates of objects in "grid/voxel space" vs. the "global model millimeter coordinate space."

In the <...>Input.h5 file,
I found the 'axis_x', 'axis_y', 'axis_z'. It seems like the axis_ arrays might correspond to the millimeter scale for each voxel/grid index, is that correct? But the axis_ arrays are 1 entry longer in each dimension - is this the "lines between the voxels"?

I also found '_Object/node_list' for each Point Sensor, which is one "large" integer. I thought this might be a linear index into the 3D voxel array, so I tried to use numpy.unravel_index() to get the x,y,z voxel indexes, and then use those indexes into the axis_x/y/z arrays to get millimeter coordinates.

This gets some results, but they don't seem to make sense with respect to the geometry. The x value seems reasonable but the y/z values don't... am I interpreting this correctly, or is there any documentation/reference for how to get the transducer/sensor coordinates in the output grid?

The end goal is to just save snapshots of single slices, on center with the transducer/sensors, for later animation. It's working now, but I'm finding the slice indexes manually by scrolling through them.

[montanaro]

This is due to the way that Sim4life defines the grid and voxels

A grid is made of nodes which make cells / voxels. The results (Pressure) values are defined at cell/voxel centers.

Therefore, if your x_axis says 0.0, 1.0, 2.0 that means your results will only have 2 values which will be located at x = 0.5 and 1.5

hint: pressure_x_axis = (x_axis[1:]+x_axis[:-1])/2
Then you'd have each axis aligned with your results

While you can work with h5py, what about trying my earlier suggestion?

"You can also create a wireframe around your volume of interest (Don't know if this works for a 2D slice / wireframe) and add that to your Field Sensors - then in PostPro you can extract only that." If you can't make a 2D wireframe, then you can just make it very thin. Although I guess you are probably working with h5py since you probably don't want to re run your simulation

[montanaro]

Once you have whatever you need as a numpy array I suggest you use one of these functions to visualize it within Sim4Life

import s4l_v1 as s4l
import s4l_v1.document as document
import s4l_v1.analysis as analysis
import s4l_v1.model as model

def visualizeArray(data, x=None,y=None,z=None, unit_name="", unit="", name="", scaling = 1.0, visualize_max = False, scalarrange = None, visualize_isosurface = False, debug = False):
	"""
	Create field in postpro to visualize data, if no axis provided
	then they are automatically assigned based on number of elements in data
	visualizeArray(data, x=None,y=None,z=None,unit_name="", unit="", name="", scaling = 0.001, goToMax = False)
	data - numpy array - 3D numpy array
	x - numpy array - 1D numpy array of x axis coordinates
	y - numpy array - 1D numpy array of y axis coordinates
	z - numpy array - 1D numpy array of z axis coordinates
	unit_name - string - unit name to be displayed
	unit - string - unit to be displayed
	name - string - name of field in postpro
	scaling - float - scaling factor for axes (default 0.001)
	goToMax - boolean - visualize slice field viewer of maximum in array
	visualize_slice
	scalarrange
	visualize_isosurface
	isosurface_value
	debug
	"""
	
	if x is None:
		x = np.arange(data.shape[0]+1)*.001
	if y is None:
		y = np.arange(data.shape[1]+1)*.001
	if z is None:
		z = np.arange(data.shape[2]+1)*.001

	grid = analysis.core.RectilinearGrid()
	grid.XAxis = np.array(x)*scaling
	grid.YAxis = np.array(y)*scaling
	grid.ZAxis = np.array(z)*scaling

	field = analysis.core.DoubleFieldData()
	field.Grid = grid
	if data.size == (len(x) * len(y) * len(z)):
		field.ValueLocation = analysis.core.eValueLocation.kNode
	elif data.size == (len(x)-1) * (len(y)-1) * (len(z)-1):
		field.ValueLocation = analysis.core.eValueLocation.kCellCenter
	else:
		print "ERROR: Grid and Data don't match"
		return False
	field.NumberOfComponents = 1
	field.NumberOfSnapshots = 1
	field.Quantity.Unit = s4l.Unit(unit)
	field.Quantity.Name = unit_name

	# Note: memory layout is such that x is fastest, z slowest dimension
	values = data.ravel('F')	
	values = values.astype(np.float64)
	field.SetField( 0, values )
	assert field.Check()

	producer = analysis.core.TrivialProducer()
	if name != "":
		producer.Description = name
	producer.SetDataObject(field)

	if visualize_max:
		sfv = analysis.viewers.SliceFieldViewer()
		sfv.Inputs[0].Connect( producer.Outputs[0] )
		sfv.Slice.Plane = sfv.Slice.Plane.enum.XY
		sfv.Update(0)
		sfv.GotoMaxSlice()
		sfv.Update(0)
		document.AllAlgorithms.Add(sfv)

		sfv = analysis.viewers.SliceFieldViewer()
		sfv.Inputs[0].Connect( producer.Outputs[0] )
		sfv.Slice.Plane = sfv.Slice.Plane.enum.YZ
		sfv.Update(0)
		sfv.GotoMaxSlice()
		sfv.Update(0)
		document.AllAlgorithms.Add(sfv)
		
		sfv = analysis.viewers.SliceFieldViewer()
		sfv.Inputs[0].Connect( producer.Outputs[0] )
		sfv.Slice.Plane = sfv.Slice.Plane.enum.XZ
		sfv.Update(0)
		sfv.GotoMaxSlice()
		sfv.Update(0)
		document.AllAlgorithms.Add(sfv)
	
	if visualize_isosurface:
		iso_surface_viewer = analysis.viewers.IsoSurfaceViewer()
		iso_surface_viewer.Inputs[0].Connect( producer.Outputs[0] )
		iso_surface_viewer.Data.Mode = iso_surface_viewer.Data.Mode.enum.QuantityRealModulus #H CHECK - maybe should just use default (delete this line)
		iso_surface_viewer.Visualization.ScalarBarVisible = False
		iso_surface_viewer.UpdateAttributes()
		iso_surface_viewer.Update(0)
		document.AllAlgorithms.Add(iso_surface_viewer)

	document.AllAlgorithms.Add(producer)
	
	return producer

def visualizeComplexArray(data, x=None,y=None,z=None, unit_name="", unit="", name="", scaling = 1.0, visualize_max = False, scalarrange = None, visualize_isosurface = False, debug = False):
	"""
	Create field in postpro to visualize data, if no axis provided
	then they are automatically assigned based on number of elements in data
	visualizeArray(data, x=None,y=None,z=None,unit_name="", unit="", name="", scaling = 0.001, goToMax = False)
	data - numpy array - 3D numpy array
	x - numpy array - 1D numpy array of x axis coordinates
	y - numpy array - 1D numpy array of y axis coordinates
	z - numpy array - 1D numpy array of z axis coordinates
	unit_name - string - unit name to be displayed
	unit - string - unit to be displayed
	name - string - name of field in postpro
	scaling - float - scaling factor for axes (default 0.001)
	goToMax - boolean - visualize slice field viewer of maximum in array
	visualize_slice
	scalarrange
	visualize_isosurface
	isosurface_value
	debug
	"""
	
	if x is None:
		x = np.arange(data.shape[0]+1)*.001
	if y is None:
		y = np.arange(data.shape[1]+1)*.001
	if z is None:
		z = np.arange(data.shape[2]+1)*.001

	grid = analysis.core.RectilinearGrid()
	grid.XAxis = np.array(x)*scaling
	grid.YAxis = np.array(y)*scaling
	grid.ZAxis = np.array(z)*scaling

	field = analysis.core.ComplexDoubleFieldData()
	field.Grid = grid
	if data.size == (len(x) * len(y) * len(z)):
		field.ValueLocation = analysis.core.eValueLocation.kNode
	elif data.size == (len(x)-1) * (len(y)-1) * (len(z)-1):
		field.ValueLocation = analysis.core.eValueLocation.kCellCenter
	else:
		print "ERROR: Grid and Data don't match"
		return False
	field.NumberOfComponents = 1
	field.NumberOfSnapshots = 1
	field.Quantity.Unit = s4l.Unit(unit)
	field.Quantity.Name = unit_name

	# Note: memory layout is such that x is fastest, z slowest dimension
	values = data.ravel('F')	
	#values = values.astype(np.complex64)
	field.SetField( 0, values )
	assert field.Check()

	producer = analysis.core.TrivialProducer()
	if name != "":
		producer.Description = name
	producer.SetDataObject(field)

	if visualize_max:
		sfv = analysis.viewers.SliceFieldViewer()
		sfv.Inputs[0].Connect( producer.Outputs[0] )
		sfv.Slice.Plane = sfv.Slice.Plane.enum.XY
		sfv.Update(0)
		sfv.GotoMaxSlice()
		sfv.Update(0)
		document.AllAlgorithms.Add(sfv)

		sfv = analysis.viewers.SliceFieldViewer()
		sfv.Inputs[0].Connect( producer.Outputs[0] )
		sfv.Slice.Plane = sfv.Slice.Plane.enum.YZ
		sfv.Update(0)
		sfv.GotoMaxSlice()
		sfv.Update(0)
		document.AllAlgorithms.Add(sfv)
		
		sfv = analysis.viewers.SliceFieldViewer()
		sfv.Inputs[0].Connect( producer.Outputs[0] )
		sfv.Slice.Plane = sfv.Slice.Plane.enum.XZ
		sfv.Update(0)
		sfv.GotoMaxSlice()
		sfv.Update(0)
		document.AllAlgorithms.Add(sfv)
	
	if visualize_isosurface:
		iso_surface_viewer = analysis.viewers.IsoSurfaceViewer()
		iso_surface_viewer.Inputs[0].Connect( producer.Outputs[0] )
		iso_surface_viewer.Data.Mode = iso_surface_viewer.Data.Mode.enum.QuantityRealModulus #H CHECK - maybe should just use default (delete this line)
		iso_surface_viewer.Visualization.ScalarBarVisible = False
		iso_surface_viewer.UpdateAttributes()
		iso_surface_viewer.Update(0)
		document.AllAlgorithms.Add(iso_surface_viewer)

	document.AllAlgorithms.Add(producer)
	
	return producer

def visualize2DArray(data, x=None,y=None,z=None, unit_name="", unit="", name="", scaling = 1., visualize_slice = False, scalarrange = None, visualize_isosurface = False, isosurface_value = None, debug = False):
	"""
	Create field in postpro to visualize data, if no axis provided
	then they are automatically assigned based on number of elements in data
	visualize2DArray(data, x=None,y=None,z=None, unit_name="", unit="", name="", scaling = 0.001, visualize = False, scalarrange = None)
	data - numpy array - 3D numpy array
	x - numpy array - 1D numpy array of x axis coordinates
	y - numpy array - 1D numpy array of y axis coordinates
	z - numpy array - 1D numpy array of z axis coordinates
	unit_name - string - unit name to be displayed
	unit - string - unit to be displayed
	name - string - name of field in postpro
	scaling - float - scaling factor for axes (default 0.001)
	visualize - bool - automatically extract field
	"""
	
	from numpy import newaxis

	# Deal with dimension of data
	if data.ndim == 2:
		data = data[:,:,newaxis]
	elif data.ndim < 2 or data.ndim > 3:
		print "Data Dimensions Error"
		return
		
	# Deal with scalar axis by turning into array
	if np.isscalar(x):
		x = [x]
	elif np.isscalar(y):
		y = [y]
	elif np.isscalar(z):
		z = [z]
	
	#If axes are not set, then make axes
	if x is None:
		x = np.arange(data.shape[0]+1)*.001 
	if y is None:
		y = np.arange(data.shape[1]+1)*.001
	if z is None:
		z = np.arange(data.shape[2]+1)*.001

	# Deal with monotonically decreasing axes
	if np.all(np.diff(x) < 0) and np.size(x) > 1:
		x = x[::-1]
		data = data[::-1]
		if debug == True:
			print "Warning: Monotonically decreasing x axes"
	if np.all(np.diff(y) < 0) and np.size(y) > 1:
		y = y[::-1]
		data = data[:,::-1]
		if debug == True:
			print "Warning: Monotonically decreasing y axes"
	if np.all(np.diff(z) < 0) and np.size(z) > 1:
		z = z[::-1]
		data = data[:,:,::-1]
		if debug == True:
			print "Warning: Monotonically decreasing z axes"
	
	grid = analysis.core.RectilinearGrid()
	grid.XAxis = np.array(x)*scaling
	grid.YAxis = np.array(y)*scaling
	grid.ZAxis = np.array(z)*scaling

	field = analysis.core.DoubleFieldData()
	field.Grid = grid
	if data.size == (len(x) * len(y) * len(z)):
		field.ValueLocation = analysis.core.eValueLocation.kNode
	elif data.size == (len(x)-1) * (len(y)-1) * (len(z)-1):
		field.ValueLocation = analysis.core.eValueLocation.kCellCenter
	else:
		print "ERROR: Grid and Data don't match"
		return
	field.NumberOfComponents = 1
	field.NumberOfSnapshots = 1
	field.Quantity.Unit = s4l.Unit(unit)
	field.Quantity.Name = unit_name

	# Note: memory layout is such that x is fastest, z slowest dimension
	values = data.ravel('F')	
	values = values.astype(np.float64)
	field.SetField( 0, values )
	assert field.Check()

	producer = analysis.core.TrivialProducer()
	if name != "":
		producer.Description = name
	producer.SetDataObject(field)
	
	# Adding a SliceSurfaceViewer
	if visualize_slice:
		sfv = analysis.viewers.SliceFieldViewer()
		sfv.Inputs[0].Connect( producer.Outputs[0] )
		if len(x) == 1:
			sfv.Slice.Plane = sfv.Slice.Plane.enum.YZ
		elif len(y) == 1:
			sfv.Slice.Plane = sfv.Slice.Plane.enum.XZ
		elif len(z) == 1:
			sfv.Slice.Plane = sfv.Slice.Plane.enum.XY
		sfv.Visualization.ScalarBarVisible = False
		if scalarrange != None:
			sfv.ScalarRange = scalarrange
		sfv.Update(0)
		document.AllAlgorithms.Add(sfv)
		
	# Adding a IsoSurfaceViewer
	if visualize_isosurface: 
		iso_surface_viewer = analysis.viewers.IsoSurfaceViewer()
		iso_surface_viewer.Inputs[0].Connect( producer.Outputs[0] )
		if isosurface_value is not None:
			iso_surface_viewer.IsoValues = (isosurface_value,)
		iso_surface_viewer.Data.Mode = iso_surface_viewer.Data.Mode.enum.QuantityRealModulus #H CHECK - maybe should just use default (delete this line)
		iso_surface_viewer.Visualization.ScalarBarVisible = False
		iso_surface_viewer.UpdateAttributes()
		iso_surface_viewer.Update()
		if scalarrange != None:
			iso_surface_viewer.ScalarRange = scalarrange
		iso_surface_viewer.UpdateAttributes()
		iso_surface_viewer.Update()
		document.AllAlgorithms.Add(iso_surface_viewer)

	document.AllAlgorithms.Add(producer)

	return producer

to use it then you would just need to do:

visualize2DArray(my_data, x_axis, y_axis, z_coordinate)

and then it should show up in the Postpro / Analysis tab for you to work this