Skip to content
  • Search
Skins
  • Light
  • Cerulean
  • Cosmo
  • Flatly
  • Journal
  • Litera
  • Lumen
  • Lux
  • Materia
  • Minty
  • Morph
  • Pulse
  • Sandstone
  • Simplex
  • Sketchy
  • Spacelab
  • United
  • Yeti
  • Zephyr
  • Dark
  • Cyborg
  • Darkly
  • Quartz
  • Slate
  • Solar
  • Superhero
  • Vapor

  • Default (No Skin)
  • No Skin
Collapse

ZMT zurich med tech

  1. Home
  2. Sim4Life
  3. Analysis & Postprocessing
  4. Handling Large Field Data Results

Handling Large Field Data Results

Scheduled Pinned Locked Moved Analysis & Postprocessing
analysismemory acoustic
8 Posts 3 Posters 1.2k Views
  • Oldest to Newest
  • Newest to Oldest
  • Most Votes
Reply
  • Reply as topic
Log in to reply
This topic has been deleted. Only users with topic management privileges can see it.
  • D Offline
    D Offline
    dbsim4
    wrote on last edited by
    #1

    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)?

    1 Reply Last reply
    0
    • L Offline
      L Offline
      LJ
      wrote on last edited by LJ
      #2

      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.

      1 Reply Last reply
      0
      • D Offline
        D Offline
        dbsim4
        wrote on last edited by
        #3

        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).

        1 Reply Last reply
        0
        • M Offline
          M Offline
          montanaro
          wrote on last edited by montanaro
          #4

          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"

          1 Reply Last reply
          0
          • M Offline
            M Offline
            montanaro
            wrote on last edited by
            #5

            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)
            
            1 Reply Last reply
            0
            • D Offline
              D Offline
              dbsim4
              wrote on last edited by dbsim4
              #6

              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.

              1 Reply Last reply
              0
              • M Offline
                M Offline
                montanaro
                wrote on last edited by
                #7

                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

                1 Reply Last reply
                0
                • M Offline
                  M Offline
                  montanaro
                  wrote on last edited by montanaro
                  #8

                  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

                  1 Reply Last reply
                  0
                  Reply
                  • Reply as topic
                  Log in to reply
                  • Oldest to Newest
                  • Newest to Oldest
                  • Most Votes


                  • Login

                  • Don't have an account? Register

                  • Login or register to search.
                  • First post
                    Last post
                  0
                  • Search