Welcome to ALPINE documentation¶
What is ALPINE?¶
ALPINE is a data and visualization software project under development as part of the Department of Energy’s Exascale Computing Project (ECP). ALPINE will deliver exascale-ready in situ data analysis and visualization capabilities to ECP science applications.
ALPINE capabilities include infrastructure, analysis and visualization algorithms, and compression. This documentation will focus on ALPINE aglorithms. Links are provided to access information and documentation about ALPINE infrastructure and zfp compression.
This is a work in progress. New algorithms and documentation will be added as needed in the future. ECP science applications interested in using ALPINE capabilities can reach out to us via ECP Confluence.
All of our software is open source and available at the links noted within the documentation.
ALPINE Algorithms¶
As simulations move to exascale machines, I/O and storage become a major limitation to scientific knowledge discovery. Exascale system concurrency is expected to grow by five or six orders of magnitude, yet I/O bandwidth is only expected to grow by two orders of magnitude. On exascale systems, the increasing I/O bottleneck will make it necessary for most simulation data to be analyzed in situ, or on the supercomputer while the simulation is running. Furthermore, to meet data bandwidth constraints, it will be necessary to sharply reduce the volume of data moved on the machine and especially the data that are exported to persistent storage. ALPINE provides a suite of in situ algorithms that will address these concerns by using adaptable data-driven techniques to downsample data, identify important features in the data, and move analysis of the data into the simulation process.
The current list of ALPINE algorithms is below. This documentation is a work in progress and some algorithms may not yet be fully documented.. _label_introduction:. ECP partner applications can find more information and contacts on the ECP ALPINE Confluence page.
In Situ Statistical Feature Detection and Post Hoc Tracking Algorithm¶
Overview¶
In Situ Statistical Feature Detection Algorithm¶
The Statistical feature detection algorithm is an in situ feature detection algorithm that processes three-dimensional (3D) particle fields in situ and transforms the data into a feature similarity field, which is stored in the disks for further post hoc analysis. Even though our current version of the algorithm works on a particle field, the algorithm can be easily applied to any regular-grid scalar data with minor modifications. By analyzing data in situ and detecting user-interested feature regions, the algorithm outputs a statistically summarized data set, which is significantly smaller in size compared to the raw particle data and such summarized data can be analyzed interactively in post hoc analysis for further feature analysis. This algorithm follows the feature-driven data reduction paradigm to achieve significant data reduction while preserving important information so that post hoc analysis can be done on the reduced data promptly.
This algorithm works on an unstructured particle field and a feature is represented as a statistical probability distribution. Representing the feature in the form of distribution allows the application scientists to specify a descriptor of the features of interest without needing to precisely define it. In many application domains, a precise description of a feature is not readily available due to the complexity of the scientific data and so a statistical technique such as this is a promising solution for feature detection. An interactive user interface can be used where the users can move a cube object freely inside the data and put it in a region where they are interested. Next, a distribution representation (currently Gaussian distribution is used, but any other distribution model can be used) is created from the data points within that selected cube region and is used as the target feature descriptor that we want to find in the data. For our current target application, the algorithm takes a particle field as input and first transforms it into a regular grid particle density field, which is a scalar field and is used as an intermediate representation. Then the density field is passed through a 3D super voxel generating algorithm, called Simple Linear Iterative Clustering (SLIC) [2], that produces super voxels from the particle density field. Now, a (Gaussian) distribution is modeled for each super voxel using the density values from respective super voxels. Finally, we use a distribution similarity measure to compute the statistical similarity between each super voxel distribution and the user-provided target feature distribution. Using these distribution similarity values for each super voxel, the algorithm finally generates a new feature similarity scalar field where each data point in the regular grid indicates its similarity with the user-given target feature. This feature similarity field is stored into disks at each time step for further detailed interactive post hoc analysis. A detailed version of this algorithm can be found in [1].
This statistical feature detection algorithm is made available as a VTK-based ParaView/Catalyst filter for the user community. We have also developed a VTKm-based filter of this algorithm, which will be made available soon and the VTKm version of the algorithm will be using accelerators to improve computational performance.
Post Hoc Feature Tracking Algorithm¶
To enable post hoc analysis of the in situ detected features, we have developed a Python-based feature tracking algorithm. This algorithm allows extraction and tracking of user-interested features so that the temporal feature dynamics can be studied in detail. The feature tracking algorithm takes the in situ generated feature similarity fields as input and then can track a user-specified feature throughout its lifetime. The tracking algorithm also detects merge and split events for the target feature so that the interactivity of the features can be studied. The analysis results are visualized in the form of a Cinema database [3] and use an interactive CinemaExplorer interface. Also, we have developed a VTK-based visualization tool that allows visual analysis of tracking results directly in 3D space.
To solve the correct correspondence between features in consecutive time steps, we use an overlap-based testing method. First, the features are segmented using a suitable feature similarity threshold, and then a connected component algorithm is applied so that each segment gets a unique feature id assigned to it. Next, given a target feature, our algorithm tracks it through space and time as long as the feature lasts within the simulation bound. To estimate the overlap amount, instead of just checking the intersection between objects, we compute the Dice similarity index which gives the amount of intersection between two features. The Dice index is computed directly in 3D space so that the correct overlap is detected. Formally given two sets of points A and B, their Dice similarity index DI(A, B) is measured as:
The Dice index value can be interpreted as matching confidence and when an abruptly low value is found, such a time step is flagged for further investigation. In our method, we track feature volumes to explore potential merge and split events. During tracking, a sudden rise in the feature volume indicates a merge event while a significant drop in feature would indicate splitting. From the tracking results, we also estimate the feature rise velocity, an important feature characteristic of interest to the domain experts. The feature rise velocity is computed from differences of the centroids of matched features from the consecutive time steps.
Getting Started¶
In Situ Statistical Feature Detection Algorithm¶
The statistical distribution-based in situ feature detection algorithm is available as a VTK filter and is deployed in situ through ParaView/Catalyst pipeline. The name of the VTK filter is “vtkPFeatureAnalysis”. This filter generates a histogram-based density field from the particle data using distributed data communication. Then SLIC super voxel algorithm is applied to generate the super voxels (clusters) from the particle density field and finally, the statistical similarity field is computed and stored into the disks for post hoc analysis. The output field can be readily visualized in ParaView/VisIT.
This filter expects two input parameters from the users and can be specified in the Catalyst python script:
- Size of the 3D super voxel/clusters (“ClusterBlockSize”) - The expected X, Y, Z size of the super voxels to be generated.
- Parameters for user-specified target feature distribution (“FeatureGaussian”) - Currently takes mean and standard deviation values for the target feature distribution.
Currently, the output dimensions of the feature similarity scalar field are set to 128X16X128. In our upcoming version, we plan to make it an input parameter.
Post Hoc Feature Tracking Algorithm¶
The feature tracking code currently takes the in situ generated feature similarity fields as input. The script needs to be pointed to the directory that contains such outputs in VTI format. Next, the tracking script also expects a feature id and a time step as input which will indicate which feature to be tracked. We also provide two separate python scripts. The first script is used to segment all the features from the VTI files and generate new VTI output files where each feature is segmented and isolated for all the time steps. The second script computes the feature attributes such as feature volume, feature mass, etc., for each segmented feature so that they can be visualized and analyzed while exploring the tracking results.
The tracking script also produces a Cinema database output of the tracked feature that contains images of the tracked feature for all the tracked time steps. The images are precomputed using a Python-based script that uses ParaView to produce images of the rendered features.
Example Result - In situ bubble detection and post hoc Tracking in MFIX-Exa simulation¶
In Situ Statistical Feature Detection Algorithm¶
The following Fig. 1 shows a representative result of this algorithm when applied on a time step of MFIX-Exa particle field. The step-by-step in situ processing of the algorithm is shown. The final bubble similarity field is stored in the disks for post hoc bubble dynamics analysis.
Fig. 1 Steps of the in situ statistical feature detection algorithm.
Post Hoc Feature Tracking Algorithm¶
Fig. 2 Cinema database showing all the extracted bubbles with their different attributes using a Parallel Coordinates plot.
To visualize the bubble tracking results and the bubble dynamics, we have developed a Cinema-based post hoc feature tracking visualization workflow that uses a customized CinemaExplorer (https://github.com/cinemascience/cinema_explorer) interface to allow the application experts to explore the tracking results. We first create a Cinema database consisting of all the bubble features. Several bubble attributes such as bubble volume, location, aspect ratio, etc. can be studied to get a general overview of all the bubbles in the simulation data. Fig. 2 shows an example of the CinemaExplorer interface where several bubbles are filtered from the Parallel Coordinates plot (PCP) at the top panel and the filtered bubbles are shown at the bottom.
Fig. 3 Cinema database showing tracking result of a single bubble over time.
In Fig. 3, we show bubble tracking results of a specific bubble, tracked both forward and backward in time so that the complete time evolution of the bubble can be explored. After tracking, we compute bubble velocity and Dice similarity index and add them to the bubble tracking attribute list so that users can filter bubbles based on such attributes. We find that the bubbles typically merge/split as they move through space and time. The bubbles generally rise upward in the domain and merge/split events can be visualized using a scatterplot between bubble volume and time. A sudden increase in bubble volume in consecutive time steps would indicate a bubble merge and a sudden decrease in bubble volume would indicate a bubble split.
Fig. 4 A scatterplot between bubble volume (Y-axis) and time steps (X-axis). The bubble suddenly increases from time step 81 to 82 which indicates a bubble merge event has happened.
In Fig. 4 we show a scatterplot using the CinemaExplorer interface (X-axis time steps, Y-axis bubble volume) where we see that at time step 82, the bubble volume increases significantly, indicating a bubble merge event.
Fig. 5 3D bubble tracking visualization.
To investigate the bubble tracking results and bubble shapes directly in the 3D domain, we have also developed a Python and VTK based visualization interface (shown in Fig. 5) where the users can interactively visualize bubble tracking results. Together with the CinemaExplorer and 3D bubble tracking interfaces, the application experts can study the time-varying bubble dynamics flexibly and interactively.
Use Case Example¶
To demonstrate the use of statistical distribution-based feature detection VTK filter, here is an example Python-based Catalyst script file:
#--------------------------------------------------------------
# Global timestep output options
timeStepToStartOutputAt=0
forceOutputAtFirstCall=False
# Global screenshot output options
imageFileNamePadding=0
rescale_lookuptable=False
# Whether or not to request specific arrays from the adaptor.
requestSpecificArrays=False
# a root directory under which all Catalyst output goes
rootDirectory=''
# makes a cinema D index table
make_cinema_table=False
#--------------------------------------------------------------
# Code generated from cpstate.py to create the CoProcessor.
# paraview version 5.8.1-1-g40b41376db
#--------------------------------------------------------------
from paraview.simple import *
from paraview import coprocessing
# ----------------------- CoProcessor definition -----------------------
def CreateCoProcessor():
def _CreatePipeline(coprocessor, datadescription):
class Pipeline:
# state file generated using paraview version 5.8.1-1-g40b41376db
# ----------------------------------------------------------------
# setup views used in the visualization
# ----------------------------------------------------------------
# trace generated using paraview version 5.8.1-1-g40b41376db
#### disable automatic camera reset on 'Show'
paraview.simple._DisableFirstRenderCameraReset()
# get the material library
materialLibrary1 = GetMaterialLibrary()
# Create a new 'Render View'
renderView1 = CreateView('RenderView')
renderView1.ViewSize = [716, 352]
renderView1.AxesGrid = 'GridAxes3DActor'
renderView1.CenterOfRotation = [0.0020009127283896633, 0.001999172356439815, 0.0019999934867559444]
renderView1.StereoType = 'Crystal Eyes'
renderView1.CameraPosition = [0.036644586455975216, 0.21705065997450937, 0.0662849960547136]
renderView1.CameraFocalPoint = [0.054171492461448235, -0.2147265944903157, 0.06031747111213411]
renderView1.CameraViewUp = [0.9991579942858331, 0.04046464155513918, 0.0067760784031188894]
renderView1.CameraParallelScale = 0.11422577498907434
renderView1.BackEnd = 'OSPRay raycaster'
renderView1.OSPRayMaterialLibrary = materialLibrary1
# register the view with coprocessor
# and provide it with information such as the filename to use,
# how frequently to write the images, etc.
coprocessor.RegisterView(renderView1,
filename='RenderView1_%t.png', freq=1, fittoscreen=0,
magnification=1, width=716, height=352, cinema={}, compression=5)
renderView1.ViewTime = datadescription.GetTime()
SetActiveView(None)
# ----------------------------------------------------------------
# setup view layouts
# ----------------------------------------------------------------
# create new layout object 'Layout #1'
layout1 = CreateLayout(name='Layout #1')
layout1.AssignView(0, renderView1)
# ----------------------------------------------------------------
# restore active view
SetActiveView(renderView1)
# ----------------------------------------------------------------
# ----------------------------------------------------------------
# setup the data processing pipelines
# ----------------------------------------------------------------
# create a new 'AMReX/BoxLib Particles Reader'
# create a producer from a simulation input
dEM07_plt00050 = coprocessor.CreateProducer(datadescription, 'inputparticles')
# create a new 'Feature Analysis' filter
featureAnalysis1 = FeatureAnalysis(Input=dEM07_plt00050)
featureAnalysis1.ClusterBlockSize = [3, 3, 3]
featureAnalysis1.FeatureGaussian = [2.0, 10.0]
# create a new 'Resample To Image'
resampleToImage1 = ResampleToImage(Input=featureAnalysis1)
resampleToImage1.SamplingBounds = [1.9417382859558333e-06, 0.003984164024624505,
1.748457857490837e-08, 0.003865203601139793, 2.523888219899246e-06, 0.003966030690060129]
# ----------------------------------------------------------------
# setup the visualization in view 'renderView1'
# ----------------------------------------------------------------
# show data from resampleToImage1
resampleToImage1Display = Show(resampleToImage1, renderView1, 'UniformGridRepresentation')
# get color transfer function/color map for 'similarity'
similarityLUT = GetColorTransferFunction('similarity')
similarityLUT.ApplyPreset('Cold and Hot', True) ####
similarityLUT.InvertTransferFunction() ####
similarityLUT.ScalarRangeInitialized = 1.0
# get opacity transfer function/opacity map for 'similarity'
similarityPWF = GetOpacityTransferFunction('similarity')i
similarityPWF.Points = [0.0, 0.0, 0.5, 0.0, 0.37078651785850525, 0.02139037474989891,
0.5, 0.0, 0.7387640476226807, 0.04812834411859512, 0.5, 0.0, 0.8904494643211365, 0.27272728085517883, 0.5, 0.0, 1.0, 1.0, 0.5, 0.0]
##similarityPWF.Points = [0.0, 1.0, 0.5, 0.0, 0.25070422887802124, 0.3375000059604645,
0.5, 0.0, 0.490140825510025, 0.0, 0.5, 0.0, 1.0, 0.0, 0.5, 0.0]
similarityPWF.ScalarRangeInitialized = 1
# trace defaults for the display properties.
resampleToImage1Display.Representation = 'Volume'
resampleToImage1Display.ColorArrayName = ['POINTS', 'similarity']
resampleToImage1Display.LookupTable = similarityLUT
resampleToImage1Display.OSPRayScaleArray = 'similarity'
resampleToImage1Display.OSPRayScaleFunction = 'PiecewiseFunction'
resampleToImage1Display.SelectOrientationVectors = 'None'
resampleToImage1Display.ScaleFactor = 0.0003982222286338549
resampleToImage1Display.SelectScaleArray = 'None'
resampleToImage1Display.GlyphType = 'Arrow'
resampleToImage1Display.GlyphTableIndexArray = 'None'
resampleToImage1Display.GaussianRadius = 1.9911111431692744e-05
resampleToImage1Display.SetScaleArray = ['POINTS', 'similarity']
resampleToImage1Display.ScaleTransferFunction = 'PiecewiseFunction'
resampleToImage1Display.OpacityArray = ['POINTS', 'similarity']
resampleToImage1Display.OpacityTransferFunction = 'PiecewiseFunction'
resampleToImage1Display.DataAxesGrid = 'GridAxesRepresentation'
resampleToImage1Display.PolarAxes = 'PolarAxesRepresentation'
resampleToImage1Display.ScalarOpacityUnitDistance = 6.888499664772514e-05
resampleToImage1Display.ScalarOpacityFunction = similarityPWF
resampleToImage1Display.SliceFunction = 'Plane'
resampleToImage1Display.Slice = 49
# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction'
resampleToImage1Display.ScaleTransferFunction.Points = [0.0, 0.0, 0.5, 0.0, 0.9267358779907227, 1.0, 0.5, 0.0]
# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction'
resampleToImage1Display.OpacityTransferFunction.Points = [0.0, 0.0, 0.5, 0.0, 0.9267358779907227, 1.0, 0.5, 0.0]
# init the 'Plane' selected for 'SliceFunction'
resampleToImage1Display.SliceFunction.Origin = [0.0019930528814552304, 0.0019326105428591842, 0.0019842772891400145]
# setup the color legend parameters for each legend in this view
# get color legend/bar for similarityLUT in view renderView1
similarityLUTColorBar = GetScalarBar(similarityLUT, renderView1)
similarityLUTColorBar.Title = 'similarity'
similarityLUTColorBar.ComponentTitle = ''
similarityLUTColorBar.Orientation = 'Vertical' ##
similarityLUTColorBar.Position = [0.862998295615945, 0.13773147335584612] ##
similarityLUTColorBar.ScalarBarLength = 0.33000000000000007 ##
# set color bar visibility
similarityLUTColorBar.Visibility = 1
# show color legend
resampleToImage1Display.SetScalarBarVisibility(renderView1, True)
# ----------------------------------------------------------------
# setup color maps and opacity mapes used in the visualization
# note: the Get..() functions create a new object, if needed
# ----------------------------------------------------------------
# ----------------------------------------------------------------
# finally, restore active source
SetActiveSource(resampleToImage1)
# ----------------------------------------------------------------
# Now any catalyst writers
xMLPImageDataWriter1 = servermanager.writers.XMLPImageDataWriter(Input=resampleToImage1)
coprocessor.RegisterWriter(xMLPImageDataWriter1, filename='ResampleToImage1_%t.pvti',
freq=1, paddingamount=0, DataMode='Appended', HeaderType='UInt64',
EncodeAppendedData=False, CompressorType='None', CompressionLevel='6')
return Pipeline()
class CoProcessor(coprocessing.CoProcessor):
def CreatePipeline(self, datadescription):
self.Pipeline = _CreatePipeline(self, datadescription)
coprocessor = CoProcessor()
# these are the frequencies at which the coprocessor updates.
freqs = {'inputparticles': [1]}
coprocessor.SetUpdateFrequencies(freqs)
if requestSpecificArrays:
arrays = [['cpu', 0], ['density', 0], ['dragx', 0], ['dragy', 0],
['dragz', 0], ['id', 0], ['mass', 0], ['omegax', 0], ['omegay', 0],
['omegaz', 0], ['omoi', 0], ['phase', 0], ['radius', 0],
['state', 0], ['velx', 0], ['vely', 0], ['velz', 0], ['volume', 0]]
coprocessor.SetRequestedArrays('inputparticles', arrays)
coprocessor.SetInitialOutputOptions(timeStepToStartOutputAt,forceOutputAtFirstCall)
if rootDirectory:
coprocessor.SetRootDirectory(rootDirectory)
if make_cinema_table:
coprocessor.EnableCinemaDTable()
return coprocessor
#--------------------------------------------------------------
# Global variable that will hold the pipeline for each timestep
# Creating the CoProcessor object, doesn't actually create the ParaView pipeline.
# It will be automatically setup when coprocessor.UpdateProducers() is called the
# first time.
coprocessor = CreateCoProcessor()
#--------------------------------------------------------------
# Enable Live-Visualizaton with ParaView and the update frequency
coprocessor.EnableLiveVisualization(False, 1)
# ---------------------- Data Selection method ----------------------
def RequestDataDescription(datadescription):
"Callback to populate the request for current timestep"
global coprocessor
# setup requests for all inputs based on the requirements of the
# pipeline.
coprocessor.LoadRequestedData(datadescription)
# ------------------------ Processing method ------------------------
def DoCoProcessing(datadescription):
"Callback to do co-processing for current timestep"
global coprocessor
# Update the coprocessor by providing it the newly generated simulation data.
# If the pipeline hasn't been setup yet, this will setup the pipeline.
coprocessor.UpdateProducers(datadescription)
# Write output data, if appropriate.
coprocessor.WriteData(datadescription);
# Write image capture (Last arg: rescale lookup table), if appropriate.
coprocessor.WriteImages(datadescription, rescale_lookuptable=rescale_lookuptable,
image_quality=0, padding_amount=imageFileNamePadding)
# Live Visualization, if enabled.
coprocessor.DoLiveVisualization(datadescription, "localhost", 22222)
Unit Testing¶
Unit test cases for the VTK filter of the feature detection algorithm has been developed and is made available as a VTK test. The test case runs the feature detection algorithm on a known MFIX-Exa test data set and several estimated output values are checked to ensure the correctness of the test. The users can find the test case under the name TestPFeatureAnalysisFilter in the appropriate VTK repository.
A similar unit test case for the VTKm filter of the algorithm will be made available along with the release of the VTKm filter.
Performance¶
The VTK filter of the feature detection algorithm currently runs on CPUs in the distributed environment. The integration pipeline has been tested at OLCF Summit. We have done an initial performance study to evaluate the computational time taken by our algorithm when deployed as an in situ VTK filter and run on CPUs in Summit using ParaView Catalyst. To improve the performance of the algorithm, we added OpenMP acceleration to our VTK filter and improved the computation timings. Since our algorithm has a component that executes on a single node, we found that by using OpenMP, we can accelerate the code. In the Fig. 6, we show the performance plot, where the X-axis shows number of OpenMP threads used and Y-axis shows time taken. As we can see, with increased number of threads, the computation time goes down and our filter takes less time to execute in situ.
Fig. 6 Performance plot of the VTK filter of the statistical feature detection algorithm. We used OpenMP acceleration to obtain higher performance.
We will soon make the VTKm-based filter of the feature detection algorithm available for the users, which will be using accelerators/GPUs and could further improve performance over the current VTK filter.
Developers¶
- Soumya Dutta (sdutta@lanl.gov)
- Li-Ta Lo
- Dan Lipsa
- Patrick O’Leary
- Berk Geveci
References¶
- Soumya Dutta, Terece L. Turton, David Rogers, Jordan M. Musser, James Ahrens, and Ann S. Almgren, In Situ Feature Analysis for Large-scale Multiphase Flow Simulations, Journal of Computational Science (Elsivier), 2022.
- Soumya Dutta, Dan Lipsa, Terece L. Turton, Berk Geveci, and James Ahrens, In Situ Analysis and Visualization of Extreme-Scale Particle Simulations, WOIV: 6th International Workshop on In Situ Visualization, 2022 (co-located with ISC High Performance Conference).
- Soumya Dutta, Jonathan Woodring, Han-Wei Shen, Jen-Ping Chen, and James Ahrens, “Homogeneity Guided Probabilistic Data Summaries for Analysis and Visualization of Large-Scale Data Sets”, 2017 IEEE Pacific Visualization Symposium (PacificVis), Seoul, Korea (South), 2017, pp. 111-120.
- Radhakrishna Achanta, Appu Shaji, Kevin Smith, Aurelien Lucchi, Pascal Fua, and Sabine Süsstrunk, “SLIC Superpixels Compared to State-of-the-art Superpixel Methods”, in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 34, no. 11, pp. 2274-2282, Nov. 2012, doi: 10.1109/TPAMI.2012.120.
- James Ahrens, Sebastien Jourdain, Patrick O’Leary, John Patchett, David Rogers, and Mark Petersen, “An Image-Based Approach to Extreme Scale in Situ Visualization and Analysis,” SC ‘14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 2014, pp. 424-434, doi: 10.1109/SC.2014.40.
Adaptive Sampling¶
Overview¶
Sampling is an in situ data reduction approach for the scalar datasets generated by the scientific simulations. This generic feature-driven data reduction method is applicable for any regular-grid datasets at this stage and is available through Ascent as a VTK-m/VTK-h filter. This generic sampling method essentially analyzes the scalar data distribution and automatically assigns importance to the scalars.
Features of the data generally span much fewer data points compared to the background (“non-interesting”) regions of the data. Utilizing this idea, our sampling approach assigns more importance to the low probability scalars and less importance to the highly frequent scalars. In addition to utilizing the global data distribution, importance can also be assigned based on the local smoothness or homogeneity of data. In other words, to capture the shapes of important data-features in a better way, we can assign more importance to the samples around the feature boundaries (i.e. high-gradient regions) by using the joint distribution of the gradient magnitudes and the scalar values. This feature can be turned-on in our sampling filter by setting “use_gradient” parameter as described below.
Here is a brief outline of the two falvors of sampling algorithms available in our filter:
- Value-based Sampling: In this approach, to achieve sampling, first, the histogram of the dataset is created. Histograms, used as a representation of the data distribution, are used for identifying and assigning importance to the individual scalars. The goal is to pick an equal number of samples from each bin. Therefore, in this sampling approach the probabilty of picking low-frequency values is higher than picking high-frequency values. More details of the approach can be found in this paper: “In Situ Data-Driven Adaptive Sampling for Large-scale Simulation Data Summarization”, Ayan Biswas, Soumya Dutta, Jesus Pulido, and James Ahrens, In Situ Infrastructures for Enabling Extreme-scale Analysis and Visualization (ISAV 2018), co-located with Supercomputing 2018.
- Value and Gradient-based Sampling: The pervious approach first decides how many samples to pick from each bin of the scalar value distribution and then selects that many samples at random. Value and Gradient-based sampling is an extension of the pervious approach where gradient magnitude information is used to decide how to select the samples from each scalar value bin. For this, the joint histogram of the scalar values and their gradient maginitudes is first calculated. For each scalar value bin, while selecting the samples more importance is assigned to the samples with high corresponding gradient magnitude. More details of the approach can be found in this paper: “Probabilistic Data-Driven Sampling viaMulti-Criteria Importance Analysis”, Ayan Biswas, Soumya Dutta, Earl Lawrence, John Patchett, Jon C. Calhoun, and James Ahrens, IEEE Transcations of Visualization and Computer Graphics (TVCG), 2020 (doi: 10.1109/TVCG.2020.3006426)
The sampling method operates by creating a point representation of the selected data points at the in situ processing phase. These data points can then be restored back to its original size/resolution by an inverse operation or these samples can be directly visualized as a preview of the dataset using visualization tools (such as Paraview).
Getting Started¶
The distribution-driven sampling approach is now available through Ascent and is implemented using VTK-m and VTK-h filters. VTK-h is needed for creating a distributed histogram and making it available to all the individual processors where the VTK-m filter is executed. Due to the dependence on data histogram, the name of the filter is “HistSampling”. Since the light-weight in situ library Ascent can be connected to simulation code, the sampling algorithm can be directly executed in situ on the supercomputers.
This easy to use sampling filter expects only three parameters from the users:
- Scalar field name (“field”) - the name of the scalar field to be sampled.
- Rate of sampling (“sample_rate”) - the amount of data to be stored after sampling. If the total bandwidth is m and the dataset size for the current time step is n, then this value can be set to n/m or less. The default is 0.1, i.e., 10% of the original data will be stored.
- Number of bins for the histogram (“bins”) - number of bins for creating the data histogram. The default is 128.
- Whether to use the gradient information during sampling (“use_gradient”) ? The default is false
These parameters can be specified by assigning value via an ascent_actions.json file.
Use Case Examples¶
To demonstrate the use of histogram-based sampling filter, here is an example ascent_actions.json file
[
{
"action": "add_pipelines",
"pipelines":
{
"pl1":
{
"f1":
{
"type": "histsampling",
"params":
{
"field": "density",
"sample_rate": "0.05",
"use_gradient": "true",
"bins": "64"
}
}
}
}
},
{
"action": "execute"
},
{
"action": "reset"
}
]
In the use case above, we have selected parameters that are appropriate to down-sample a scalar field named “density”.
- “field”: We choose the scalar field “density” that is generated by the scientific simulation.
- “sample_rate”: The example sampling rate is equal to 0.05, i.e., 5% of the original data will be stored at each time step.
- “bins”: In this example we used 64 bins to create the histogram that will be used for sampling the scalar field.
- “use_gradient”: In this example we are using the Value and Gradient-based sampling approach because the flag is set to true.
Performance¶
To improve overall performance of the sampling operation, this filter is implemented in vtkm/vtkh and is capable of utilizing hardware accelerators. Further, creation of distributed histogram is additive across processors that also helps improve scalability.
Lagrangian Analysis¶
Overview¶
Lagrangian analysis is an in situ data reduction operator used for time-dependent vector field data generated by a simulation code. With the objective of storing/representing fluid dynamics data in its Lagrangian representation, the Lagrangian analysis functionality is implemented as a VTK-m filter. The filter operates by placing seeds and calculating the corresponding particle trajectories in the flow volume. These particle trajectories encode the underlying behavior of the flow field. Calculating and extracting a Lagrangian representation of a flow field offers significantly improved accuracy-storage propositions for time-dependent flow visualization compared to the traditional (Eulerian) method. Thus, the Lagrangian analysis filter enables data reduction of large vector fields while maintaining high data integrity.
Lagrangian analysis operates as a two-phase approach. The first phase is performed in situ and involves calculating and extracting flow field information. The second phase is performed post hoc and involves interpolating the extracted information to perform time-dependent flow visualization.
To calculate a Lagrangian representation in situ, we uniformly seed particles in the volume and advect the particles for nonoverlapping pre-determined intervals of time. The calculated particle trajectories are referred to as basis flows. The extracted basis flows are accurate given they use the complete spatial and temporal resolution of the simulation vector field during advection in situ. For each basis flow, the start location and end location of the particle trajectory is known.
The stored basis flow information can be interpolated post hoc to perform exploratory time-dependent flow visualization. Complete pathlines can be calculated by concatenating results from batches of basis flows.
Getting Started¶
Lagrangian analysis functionality is accessible via Ascent. In situ extraction of time-dependent vector field data is implemented using VTK-m, VTK-h filters and a VTK-m particle advection worklet. The Ascent library is light-weight infrastructure that is directly integrated with the simulation code. Depending on the particular use case (for example, operating on a single node), Lagrangian analysis can also be directly performed using the VTK-m filter.
To perform in situ extraction, several parameters can be specified to the Lagrangian analysis filter. These include:
- Vector field name (“field”) - the name of the vector field whose Lagrangian representation is to be calculated
- Step size (“step_size”)- the unit of time a particle advances between each cycle
- Write frequency (“write_frequency”) - the number of cycles between writing basis flow information to disk, i.e., the interval length
- Custom seeding resolution (“cust_res”) - 0 or 1. The default value 0 results in a 1:1 seed resolution, i.e., number of seeds to grid points. The value of 1 allows the user to specify a custom seed resolution. This is achieved by specifying the number of seeds to grid points along each axis:
- Reduction factor in X direction (“x_res”)
- Reduction factor in Y direction (“y_res”)
- Reduction factor in Z direction (“z_res”)
These parameters can be specified by manually setting the values using filter object functions or by passing values via an ascent_actions.json file.
Use Case Examples¶
To demonstrate the use of Lagrangian analysis, we use the Cloverleaf3D proxy application that is available as an example integration with the Ascent library.
Below is an example ascent_actions.json file that can be used to extract time-dependent vector fields in a Lagrangian representation from the Cloverleaf3D. The pipeline of operations to be performed contains the name of the filter - “lagrangian”, followed by the specific parameters required to configure the filter. For the Cloverleaf3D proxy application, a vector field named “velocity” is generated. For other simulation codes, the respective vector field name will need to be used. We can specify the step size used for advection (in this example we use “0.02”). If cust_res is set to 0, then the number of particles used is the same as the number of grid points, i.e., a 1:1 configuration. In our example, we set cust_res to 1, indicating we do not want to use the default seed resolution. We set x_res, y_res and z_res to 2, i.e., a particle is placed for every other grid point in each axis direction. Thus, for this custom seed resolution we are using 1:8 particles to grid points.
[
{
"action": "add_pipelines",
"pipelines":
{
"pl1":
{
"f1":
{
"type": "lagrangian",
"params":
{
"field": "velocity",
"step_size": 0.02,
"write_frequency": 25,
"cust_res": 1,
"x_res": 2,
"y_res": 2,
"z_res": 2
}
}
}
}
},
{
"action": "execute"
},
{
"action": "reset"
}
]
Additionally, in the event the user does not want to use an ascent_actions.json file, the code can be directly instrumented. The C++ code snippet below demonstrates the same Lagrangian analysis filter configuration as the instance above. Further, this example utilizes Ascent relay extract functionality to store the calculated basis flow data.
Ascent ascent;
Node ascent_opts;
ascent_opts["runtime/type"] = "ascent";
ascent.open(ascent_opts);
conduit::Node mesh_data;
// Populate mesh_data;
conduit::Node extracts;
extracts["e1/type"] = "relay";
extracts["e1/pipeline"] = "pl1";
extracts["e1/params/path"] = output_file_path;
conduit::Node pipelines;
// pipeline 1
pipelines["pl1/f1/type"] = "lagrangian";
// filter knobs
conduit::Node &lagrangian_params = pipelines["pl1/f1/params"];
lagrangian_params["field"] = "velocity";
lagrangian_params["step_size"] = 0.02;
lagrangian_params["write_frequency"] = 25;
lagrangian_params["cust_res"] = 1;
lagrangian_params["x_res"] = 2;
lagrangian_params["y_res"] = 2;
lagrangian_params["z_res"] = 2;
conduit::Node actions;
// add the pipeline
conduit::Node &add_pipelines = actions.append();
add_pipelines["action"] = "add_pipelines";
add_pipelines["pipelines"] = pipelines;
// add the extracts
conduit::Node &add_extracts = actions.append();
add_extracts["action"] = "add_extracts";
add_extracts["extracts"] = extracts;
// execute
conduit::Node &execute = actions.append();
execute["action"] = "execute";
// reset
conduit::Node &reset = actions.append();
reset["action"] = "reset";
ascent.publish(mesh_data);
ascent.execute(actions);
ascent.close();
Parameter Value Selection
In the use cases above, we have selected parameters that are appropriate to extract time-dependent vector field information from the Cloverleaf3D proxy application.
- “field”: We choose the vector field “velocity” that is generated by the Cloverleaf3D application.
- “step_size”: The example step size we use is equal to the duration that the simulation advances every cycle, i.e., one cycle corresponds to the simulation advancing by 0.02 units of time.
- “write_frequency”: Depending on how often the user wants to save data, an appropriate value for write frequency should be selected. Lagrangian basis flows represent an interval of time. In the example above, we select a value of 25. For the Cloverleaf3D application this would be considered a short/medium interval length. If the simulation executes for approximately 200 cycles - then 8 sets of basis flows, each representing an interval of 25 cycles, would be stored.
- “cust_res”: Custom resolution needs to be enabled (value = 1) in order to achieve a data reduction. By default, the Lagrangian filter will use a 1:1 configuration, i.e., it will use as many particles as grid points.
- “x_res”, “y_res”, “z_res”: These parameters control the degree on data reduction. We recommend using either values of 2, i.e., every other grid point in each direction or 1:8 data reduction, or 3, i.e., every third grid point in each direction or 1:27 data reduction. Setting each parameter to 4 would result in a 1:64 data reduction and could reduce data integrity.
Performance¶
To improve overall performance of the Lagrangian analysis the filter particle advection worklet is capable of utilizing hardware accelerators for improved performance. Further, we adopt a communication-free model and eliminate integration of basis flow trajectories across node boundaries to improve scalability.
Post Hoc Reconstruction¶
To utilize the extracted Lagrangian representation, a reconstruction filter is available in VTK-h and accessible through Ascent.
Required Input Format Information: The VTK-h filter currently accepts “.vtk” structured grid files with a “displacement” vector field and “valid” scalar field defined at each grid point. Further, the current setup requires the files follow a particular naming format in order to support loading data generated across multiple ranks and time-stamps. Files should be named “Lagrangian_Structured_RANK_INTERVAL.vtk”, where RANK specified the rank that generated the file and INTERVAL specified the interval the file encoded. For example, for particle trajectories calculated starting at cycle 0 and stored at cycle 10, and generated by rank 5, the file would be named Lagrangian_Structured_5_10.vtk. For the next interval, i.e., 10 to 20, the file would be named Lagrangian_Structured_5_20.vtk.
The second required information is a file containing seed locations. Currently, this information is read in from a text file following the format “X1 Y1 Z1” on each line, where each line specifies the location of one seed particle.
conduit::Node pipelines;
// pipeline 1
pipelines["pl1/f1/type"] = "lagrangian_interpolation";
// filter knobs
conduit::Node &lagrangian_params = pipelines["pl1/f1/params"];
lagrangian_params["field"] = "braid";
lagrangian_params["radius"] = 0.1;
lagrangian_params["num_seeds"] = 1000;
lagrangian_params["interval"] = 10;
lagrangian_params["start_cycle"] = 10;
lagrangian_params["end_cycle"] = 100;
lagrangian_params["seed_path"] = "path/to/seed/file/Seed.txt";
lagrangian_params["basis_path"] = "path/to/extracted/data/";
lagrangian_params["output_path"] = "path/to/output/data/Pathlines_"; // Name of output file begins with "Pathlines_"
conduit::Node actions;
// add the pipeline
conduit::Node &add_pipelines = actions.append();
add_pipelines["action"] = "add_pipelines";
add_pipelines["pipelines"] = pipelines;
//
// Run Ascent
//
Ascent ascent;
Node ascent_opts;
ascent_opts["runtime/type"] = "ascent";
ascent_opts["mpi_comm"] = MPI_Comm_c2f(comm);
ascent.open(ascent_opts);
// create dummy mesh using conduit blueprint
Node n_mesh;
conduit::blueprint::mesh::examples::braid("hexs",
2,
2,
2,
n_mesh);
// publish mesh to ascent
ascent.publish(n_mesh);
actions.append()["action"] = "execute";
Node &reset = actions.append();
reset["action"] = "reset";
// execute
ascent.execute(actions);
ascent.close();
Parameter Information The post hoc reconstruction routine accepts multiple parameters. We provide information for each of these below.
- “radius”: Each grid point with a valid = 0 is reconstructed using a scattered point interpolation. The radius determines which neighboring points to use for reconstruction.
- “num_seeds”: Specifies the number of seed particles to be reconstructed and specified in the seed file.
- “interval”: Temporal subsampling rate at which files were stored to disk.
- “start_cycle”: This parameter specifies the first interval that should be loaded from disk and used for reconstruction. Will typically be equal to interval, unless user wants to reconstruct pathlines starting at a different time in the simulation.
- “end_cycle”: This parameter specifies the last interval that should be reconstructed.
- “seed_path”: Path to the seed file.
- “basis_path”: Path to the extracted Lagrangian data.
- “output_path”: Specified the path to store calculated particle trajectories as well as allows specifying the output file name.
Topological Analysis¶
Overview¶
Isosurfaces often have a useful interpretation in the application domain making them an ubiquitous visual and data analysis technique, such as finding the molecular boundaries in chemical simulations. Topological analysis identifies relevant isosurfaces in data, and the contour tree summarizes the change in isocontours as isovalue changes, with individual connected isosurface components. Both of these techniques enable automatic data selection based on finding relevant contours. Furthermore, contour tree simplification can identify the most relevant topological features, thereby facilitating automatic data reduction by focusing analysis on the most important contours.
Background¶
The following section provides background information on the contour tree, simplifying the contour tree and its use in selecting relevant contours. While it is useful to understand this background, it is not strictly required to use contour tree based isovalue selection and you may skip to Getting Started if you are not interested in this background.
Contour Trees¶
Given a function from the simulation domain to a range of scalar values, a common way to visualize that function is to pick an isovalue and extract isoliness (2D) or isosurfaces (3D), i.e., lines/surfaces that connect all locations in space where the function assumes the isovalue. When considering isolines/isosurfaces, one property often of particular interest is how many connected components comprise the isocontour/isosurface. The contour tree, described in the following, provides this information for all scalar values in the range. To make this description more succinct, we use the term contour to refer to a single connected component of the isosurface.
Example of a contour tree for a simple 3D toy data set to illustrate concepts.
Fig. 7 shows the contour tree for a simple toy data set to illustrate concepts. Consider the evolution of the contours shown in the figure. As the isovalue increases, first the gay contour and later the red contour appear and subsequently merge into a single contour (colored gray in the image). For a higher isovalue, this contour splits into two contours (colored gray and blue in the image. Each of these contours splits a second time, the gray contour into two contours colored gray and green in the image and the blue contour into two contours colored blue and orange in the figure. Finally, all contours disappear at slightly different values.
The contour tree tracks this evolution of contours. Common graph layouts of the contour tree choose node height (y-coordinate) to correspond to the isovalue at which the corresponding event occurs. In the image, all edges in the contour tree share the color of the contours they represent. Degree one nodes correspond to appearing and disappearing contours and degree three (or higher) nodes correspond to contours merging or splitting. We note that the contour tree does not only encode how many contours exist for a given isovalue, but also their relationship to each other, i.e., it encodes the identify of which contours are involved in a merge or split event.
Contour Tree Simplification¶
Despite describing high-level isosurface behavior, contour trees can quickly become very complex.
Despite providing a high-level summary over isosurface behavior, the contour tree can quickly become very complex, even for simple data sets, as shown in Fig. 8. However, since the contour tree defines a relationship between contours, it is possible to further simplify and reduce it to the most relevant contours.
The contour tree from Fig. 8 can be simplified using various simplification metrics.
Fig. 9 shows the contour tree from Fig. 8 using various simplification metrics. A simplification metric is a means to rank the importance of individual contours and “prune” away unimportant branches of the contour tree. Three common choices for the simplification metric are persistence, volume and hypervolume. Persistence is the height of a branch in the contour tree, corresponding to how long this contour lives as a distinct entity as the isovalue is changed. Volume corresponds to contour volume, approximated as the number of mesh points it comprises (for a regular, rectilinear grid) and hypervolume corresponds to the integrated value of the scalar variable inside a contour, approximated as the Riemann sum. Based on these simplification metrics it is possible to rank all leaf branches in the contour tree and simplify it. At each node we compare the “weights” of the two branches leading to it and “prune” the the branch with the lower weight. For example, if we consider persistence in the contour tree of numref:ctfig, branches would be pruned in the order green, orange, red. Two types of simplification are common: (i) simplify until all remaining branches have a weight above a given threshold or (ii) simplify until only a specified number of “most important” branches remain.
Isovalue Selection¶
Our current approach for selecting the n most important isovalues proceeds as
follows. We simplify the contour tree down to the n+1 most important
branches. Subsequently, we pick a value that is one 
away from
the value at the branch point. This ensures that the contour has the maximum
possible size. Please note that we select the most relevant contours and then
use their isovalue to extract an entire isosurface and not just the contour. We
are currently working on a filter that will extract just the relevant contours.
Getting Started¶
Contour tree-based isosurface selection is available via Ascen. However,
currently Ascent needs to be built with special options (enabling building vtk-m
with MPI support) to enable this capability. To build Ascent with contour tree
support, follow the instructions for building Ascent via uberenv
but ensure that vtk-m is built with MPI support by adding --spec="%gcc
^vtkm+mpi" to the invocation on uberenv.py, e.g., to build for Linux:
python scripts/uberenv/uberenv.py --prefix uberenv_libs \
--spec "%GCC ^vtk-m+mpi"
When Ascent is built with MPI support, contour tree-based isosurface selection
is available through Ascent’s contour filter. Use the levels parameter
to select the desired number of isosurface levels and set the parameter
use_contour_tree to true to enable isovalue selection via contour tree
(instead of using equidistant isovalues. (Note: Currently the parameters for
contour selection are hard-coded and use volume-based simplification for single
compute node runs and persistence-based simplification for MPI runs. We are
working on making isovalue selection more customizable and are looking for
collaborators using this method to help us identify useful parametrizations and
parameter choices. Furthermore, our current approach may result in one isovalue
being chosen multiple times, e.g., for highly symmetric datasets. Thus, the
resulting isosurface will be extracted for up to levels isovalues, but
possibly for fewer. We are working on resolving this issue.)
Use Case Examples¶
The following use case examples plots the scalar variable named Ex using
equidistant 15 isovalues and saves the results as levels_tttt.png (with
tttt being the time step) and using 15 isovalues selected based on the contour tree
and saves the result as smart_ttt.png.
[
{
"action": "add_pipelines",
"pipelines":
{
"p1":
{
"f1":
{
"type" : "contour",
"params" :
{
"field" : "Ex",
"levels": 15
}
}
},
"p2":
{
"f1":
{
"type" : "contour",
"params" :
{
"field" : "Ex",
"levels": 15,
"use_contour_tree": "true"
}
}
}
}
},
{
"action": "add_scenes",
"scenes":
{
"s1":
{
"image_prefix": "levels_%04d",
"plots":
{
"p1":
{
"type": "pseudocolor",
"pipeline": "p1",
"field": "Ex"
}
}
},
"s2":
{
"image_prefix": "smart_%04d",
"plots":
{
"p1":
{
"type": "pseudocolor",
"pipeline": "p2",
"field": "Ex"
}
}
}
}
},
]
Performance¶
To improve the performance of contour tree calculation can use multithreading on compute nodes or accelerators such as GPUs (via vtk-m). We are working on making this functionality more easily accesible on machines like Summit.
Developers¶
- David Camp (LBNL)
- Hamish Carr (University of Leeds)
- Oliver Rübel (LBNL)
- Gunther H. Weber (LBNL)
Rotation Invariant Pattern Detection¶
Overview¶
Pattern detection can be used to identify known features in a simulation in situ to reduce the amount of data needing to be written to disk. For simulations where physically meaningful patterns are already known, the orientation of the pattern may not be known a priori. Pattern detection can be unnecessarily slowed if the pattern detection algorithm must search for all possible rotated copies of a pattern template. Therefore, rotation invariance is a critical requirement. Moment invariants can achieve rotation invariance without the need for point to point correlations, which are difficult to generate in smooth fields. For an introduction to moment invariants, we recommend:
Flusser, J., Suk, T., & Zitová, B. (2016). 2D and 3D Image Analysis by Moments. John Wiley & Sons.
ALPINE has implemented two VTK filters that together are able to perform rotation invariant pattern detection. The algorithm upon which the moment invariants pattern detection is based can be found in:
Bujack, R., & Hagen, H. (2017). Moment Invariants for Multi-Dimensional Data. In Modeling, Analysis, and Visualization of Anisotropy (pp. 43-64). Springer, Cham.
Algorithm Description¶
The input to pattern detection algorithm consists of three pieces which must match in data type and dimensionality:
- A 2D or 3D vtkImageData dataset of scalars, vectors, or matrices which will be searched for:
- A vtkImageData pattern
- A vtkImageData grid that defines the subset of (1) that is searched for the pattern; note that this may be the full dataset.
The first VTK filter, vtkComputeMments computes the moments while the second filter, vtkMomentsInvariants, performs the normalization based on the given pattern and computes the similarity. The architecture of the filter, with inputs and outputs, can be found in the following figure.
Extensions¶
The MomentInvariants module contains several extra algorithms and helper classes.
vtkMomentsHelper : a class that provides functions for the moments computation that
will be needed by vtkComputeMoments and vtkMomentInvariants.
vtkMomentsTensor : a class that provides the functionality to treat tensors of arbitrary
dimension and rank. It supports addition, outer product,
and contractions.
vtkSimilarityBalls : a filter that takes the similarity field produced by
vtkMomentInvariants and computes the local maxima in
space plus scale and produces the output
localMaxSimilarity that contains the similarity value
together with the corresponding radius at the maxima.
All other points are zero.
For further visualization, vtkSimilarityBalls also produces two output fields that encode the radius through drawing a solid ball or a hollow sphere around those locations. The second input, i.e. the grid, steers the resolution of the balls. It is helpful if its extent is a multiple of the first input’s. Then, the circles are centered nicely. The spheres/circles are useful for 2D visualizations, as they can be laid over a visualization of the field. The balls are good for 3D volume rendering or steering of the seeding of visualization elements.
The 2D visualzation is described in:
Bujack, R., Hotz, I., Scheuermann, G., & Hitzer, E. (2015). Moment invariants for 2D flow fields via normalization in detail. IEEE transactions on visualization and computer graphics, 21(8), 916-929
and the 3D counterpart in:
Bujack, R., Kasten, J., Hotz, I., Scheuermann, G., & Hitzer, E. (2015, April). Moment invariants for 3D flow fields via normalization. In Visualization Symposium (PacificVis), 2015 IEEE Pacific (pp. 9-16). IEEE.
A schematic overview of the use of vtkSimilarityBalls with example images is given in the following Figure.
vtkReconstructFromMoments : a filter that takes the momentData, as produced by
vtkComputeMoments or vtkMomentInvariants, and a grid.
It reconstructs the function from the moments (similar to
reconstructing a function from the coefficients of a Taylor
series). For the reconstruction, we first orthonormalize the
moments. Then, we multiply the coefficients with their
corresponding basis function and add them up.
Use Case Example - paraview-vis.py¶
The pattern detection algorithm can be demonstrated using ALPINE’s Ascent infrastructure and its built-in example integration, the Cloverleaf3D proxy application (http://uk-mac.github.io/CloverLeaf3D/).
In this case, the ascent_actions.json (below) points to a python script, paraview-vis.py that describes the visualization:
[
{
"action": "add_extracts",
"extracts":
{
"e1":
{
"type": "python",
"params":
{
"file": "paraview-vis.py"
}
}
}
},
{
"action": "execute"
},
{
"action": "reset"
}
]
An example python script and pattern can be found in: Instructions for In Situ ParaView Vis using Ascent Extract Interface https://github.com/danlipsa/ascent/tree/moment-invariants/src/examples/paraview-vis.
paraview-vis-cloverleaf3d-momentinvariants.py : example script calling the moments invariant
algorithm; make a symbolic link pointing
paraview-vis.py to this script.
expandingVortex.vti : example vtkImageData pattern for CloverLeaf3D
The pattern detection workflow was run in Ascent through ParaView. The images show the output of the algorithm for a vortex pattern for a single timestep of Cloverleaf running in parallel. On the left, we show the pattern visualized through streamlines. In the center, the 3D similarity output is volume rendered with red corresponding to high similarity to the pattern. On the right, we put a slice with line integral convolution (LIC) through the 3D data at the location of the strongest matches to verify the result.
Use Case Example - Bubble-finding¶
One example of using rotational invariant pattern detection is for data reduction in an MFIX-Exa bubbling bed simulation. The pattern that is used for the search is a simple density boundary – particles on one side, no particles on the other (left-hand image). The middle image shows the original dataset while the right-hand image shows the bubbles found by the rotational invariant pattern detection algorithm saving only 5% of the original data.
Repository Information¶
The moment invariants pattern detection code is found within the Kitware GitLab:
VTK: https://gitlab.kitware.com/vtk/vtk (dc2d04cdd3167d0a0aa95bc3efffa13f26c98516)
MomentInvariants: https://gitlab.kitware.com/vtk/MomentInvariants (df81d17f941989d9becdbcf10413e53af7a7ab10) Includes unit testing instructions.
Instructions for In Situ ParaView Vis using Ascent Extract Interface: https://github.com/danlipsa/ascent/tree/moment-invariants/src/examples/paraview-vis
The moment invariant pattern detection algorithm was developed by Roxana Bujack and Karen Tsai at Los Alamos National Laboratory and Dan Lipsa at Kitware, Inc.
Task Based Feature Detection¶
Overview¶
Segmentations of the domain derived from merge trees have been shown to efficiently encode threshold-based features. For example, segmented merge trees can be used to extract extinction regions defined as areas of high scalar dissipation in turbulent combustion simulations or also burning cells in turbulent combustion, eddies in the oceans, and bubbles in Raleigh-Taylor instabilities. Segmented merge-trees provide two key advantages over traditional threshold-based segmentation techniques: 1) they efficiently encode a wide range of possible segmentations, and 2) they allow for the selection of localized thresholds.
We implemented our feature extraction algorithm using a task based approach where we define the algorithm using an abstract task graphs that can be execute on multiple runtimes (using the BabelFlow library). The computation is composed of four types of tasks: local computation, join, correction and segmentation. The local computation at the leaf nodes (of the task graph) receive the input data and produces two outputs: a local tree and a boundary tree. These are sent to the correction and join tasks, respectively. All but the leaf tasks of the reduction perform the join of two (or more) boundary trees and send the other boundary tree to the next join and an augmented boundary tree to as many correction stages as there are leaves in the subtree of the join. The correction uses the augmented boundary tree and the local tree to update the local tree and sends it to the next correction stage. Once all joins have been passed to a correction, each local tree is passed to a final segmentation.
Getting Started¶
The task based feature segmentation is usable through Ascent when built via Spack/uberenv enabling the variant “+babelflow”. The filter is called “bflow_pmt” and it accepts the following set of parameters:
- task: the name of the algorithm to execute, now only “bflow_pmt” is currently supported (default = “bflow_pmt”)
- field: name of the input field (scalar and float64)
- fanin: defines the reduce/broadcast factor to use (e.g., fanin=2 will merge two merge trees at a time), it is preferrable to use 2 for 2D datasets and 8 for 3D datasets
- threshold: the threshold for the merge tree construction, all value below this thershold will be discarded
- gen_segment: generate (1) or not (0) a new field containing the segmentation named “segment”
Some optional parameters are:
- in_ghosts: number of ghost cells used by the input domain decomposition (default = [1,1,1,1,1,1])
- ugrid_select: if dataset contains multiple uniform grids, select the index of the grid to use for the analysis (default = 0)
Task-based compositing
Another feature of BabelFlow is task-based, parallel compositing of images. This feature is available as Ascent extract called “bflow_comp”. The extract is designed to be used with a Devil Ray filter called “dray_project_colors_2d”. This filter renders images of local data and passes them down the Ascent’s pipeline as Blueprint uniform mesh data with two fields. One field, called “color_field”, is the pixel data and the other one, called “depth_field”, is the depth data. This designs allows us to add future renderers and support other potential compositing patterns. The compositing extract accepts the following set of parameters:
- pipeline: the name of the pipeline on which to apply this extract
- color_field: the name of the field with the pixel data
- depth_field: the name of the field with depth data
- image_prefix: the prefix for the output image name
- compositing: specifies which compositing technique to use, 0 means single reduce tree, 1 means binary swap, and 2 means use the radix-k algorithm
- fanin: defines the reduce factor to use when image fragments are reduced into one image in the end of the Radix-k exchange pattern; for higher process counts it is recommended to use higher fanin such as 4 or 8 to reduce the number of levels in the reduce tree
- radices: array of radices for radix-k algorithm
Relevance computation
One extension on top of the feature segmentation filter is computing the relevance field. Relevance, which is a value in the range [0.0, 1.0], is computed from an existing merge tree and aims to rank points according to their relative distance in function space from the most dominant point in their neighborhood. Essentially, it reduces the potetially wide dynamic range of the original data into a fixed range that allows highlighting both local and global features.
The dataflow graph for the relevance computations is constructed by composing the PMT dataflow with a dataflow to compute the global merge tree and a layer to compute the relevance field. For the relevance computation the user needs to use the same “bflow_pmt” type as for the plain PMT filter, but there is an additional parameter called “rel_field” that needs to be switched on.
Use Case Examples¶
Input and output
Limitations: This filter currently supports only uniform grids and scalar fields (with datatype float64). The filter can generate a new field named “segment” containing a new uniform grid (equivalent to the input field’s one) containing the segmentation values.
Example
Here is an actions file which takes as input a scalar field “mag” and perform a topological segmentation using a threshold (e.g., 0.0025). The segmentation is stored by the filter in a new variable named “segment”, in the following example pipeline this variable saved by relay extract to a blueprint file.
-
action: "add_pipelines"
pipelines:
pl2:
f1:
type: "bflow_pmt"
params:
field: "mag"
fanin: 2
threshold: 0.0025
gen_segment: 1
- # rel_field: 1 # Optional parameter to enable relevance computation
action: “add_extracts” extracts:
- e1:
type: “relay” pipeline: “pl2” params:
path: “seg” protocol: “blueprint/mesh/hdf5” fields:
- “segment”
The equivaled pipeline in C++:
Ascent a;
conduit::Node ascent_opt;
ascent_opt["mpi_comm"] = MPI_Comm_c2f(MPI_COMM_WORLD);
ascent_opt["runtime/type"] = "ascent";
a.open(ascent_opt);
Node pipelines;
pipelines["pl1/f1/type"] = "bflow_pmt";
pipelines["pl1/f1/params/field"] = "mag";
pipelines["pl1/f1/params/fanin"] = 2;
pipelines["pl1/f1/params/threshold"] = 0.0025;
pipelines["pl1/f1/params/gen_segment"] = 1;
Node extract;
extract["e1/type"] = "relay";
extract["e1/pipeline"] = "pl1";
extract["e1/params/path"] = "seg";
extract["e1/params/protocol"] = "blueprint/mesh/hdf5";
extract["e1/params/fields"].append() = "segment";
Node actions;
Node &add_extract = actions.append();
add_extract["action"] = "add_extracts";
add_extract["extracts"] = extract;
actions.append()["action"] = "execute";
Node &add_pipelines = actions.append();
add_pipelines["action"] = "add_pipelines";
add_pipelines["pipelines"] = pipelines;
a.execute(actions);
a.close();
Task-based compositing
Below is the actions file for the compositing pipeline with Devil Ray renderer and BabelFlow compositing:
-
action: "add_pipelines"
pipelines:
pl1:
f1:
type: "dray_project_colors_2d"
params:
field: "density"
color_table:
name: "cool2warm"
camera:
azimuth: -30
elevation: 35
-
action: "add_extracts"
extracts:
e1:
type: "bflow_comp"
pipeline: pl1
params:
color_field: "colors"
depth_field: "depth"
image_prefix: "comp_img"
fanin: 2
compositing: 2
radices: [2, 4]
The equivaled pipeline in C++:
Ascent a;
conduit::Node ascent_opt;
ascent_opt["mpi_comm"] = MPI_Comm_c2f(MPI_COMM_WORLD);
ascent_opt["runtime/type"] = "ascent";
a.open(ascent_opt);
std::vector<int64_t> radices({2, 4});
Node pipelines;
pipelines["pl1/f1/type"] = "dray_project_colors_2d";
pipelines["pl1/f1/params/field"] = "density";
pipelines["pl1/f1/params/color_table/name"] = "cool2warm";
pipelines["pl1/f1/params/camera/azimuth"] = -30;
pipelines["pl1/f1/params/camera/elevation"] = 35;
Node extract;
extract["e1/type"] = "bflow_comp";
extract["e1/pipeline"] = "pl1";
extract["e1/params/color_field"] = "colors";
extract["e1/params/depth_field"] = "depth";
extract["e1/params/image_prefix"] = "comp_img";
extract["e1/params/fanin"] = 2;
extract["e1/params/compositing"] = 2;
extract["e1/params/radices"].set_int64_vector(radices);
Node actions;
Node &add_extract = actions.append();
add_extract["action"] = "add_extracts";
add_extract["extracts"] = extract;
actions.append()["action"] = "execute";
Node &add_pipelines = actions.append();
add_pipelines["action"] = "add_pipelines";
add_pipelines["pipelines"] = pipelines;
a.execute(actions);
a.close();
Performance¶
Our implementation using a task based abstraction layer (BabelFlow) allows for easy portability of analytics over multiple runtimes (i.e., the current implementation uses MPI). To avoid large global communications and improve overall performance reduction and broadcast communications are defined in the task graph as k-way reduction/broadcast trees (where k in the number of input/output nodes at each level of the tree).
Developers¶
- Peer-Timo Bremer (LLNL)
- Xuan Huang (University of Utah)
- Steve Petruzza (University of Utah)
- Sergei Shudler (LLNL)
ALPINE Infrastructure¶
ALPINE is delivering its capabilities through multiple in situ infrastructures and post-processing applications. Detailed information can be found by clicking on the included links.
- ParaView and Catalyst, ParaView’s in situ library.
- VisIt and LibSim, VisIt’s in situ library.
- Ascent, a flyweight in situ infrastructure.
zfp Compression¶
zpf is a floating-point array primitive using very high-speed, lossy (but optionally error-bounded) compression to significantly reduce data volumes. zfp reduces I/O time and off-line storage requirements by 1-2 orders of magnitude depending on accuracy requirements, as dictated by user-set error tolerances. Unique among data compressors, zfp also supports constant-time read/write random access to individual array elements from compressed storage. zfp’s compressed arrays can often replace conventional arrays in existing applications with minimal code changes, allowing, for example, the user to store tables of floating-point data in compressed form that otherwise would not fit in memory. When used in numerical computations, zfp arrays provide a fine-grained knob on precision while achieving accuracy comparable to IEEE floating point at half the storage, reducing both memory usage and bandwidth.
More information can be found on zfp compression on the LLNL zfp website or zfp GitHub page.
The ALPINE Team¶
ALPINE is a collaboration between teams at the following institutions:
- Los Alamos National Laboratory
- Lawrence Livermore National Laboratory
- Lawrence Berkeley National Laboratory
- Oak Ridge National Laboratory
- Sandia National Laboratories
- Argonne National Laboratory
- Kitware
- University of Oregon
- University of Utah
- University of Leeds