Run bubble burst across all timesteps

%% Cell type:code id:0e02bb89-8003-47a2-839f-58536b067490 tags:

``` python

# Some parameters for the behavior of the extraction

# Specifies the spacing of the sampling grid.

lattice_parameter = 3.177

# The radius to use for a He splat

He_radius = 0.01

# Minimum number of He atoms that have to escape for a bubble to be

# considered burst.

burst_size = 25

# Minimum fraction of He atoms in a bubble that have to escape for

# the buble to be considered burst.

burst_fraction = 0.75

# Number of time steps to pass to check for escaped atoms.

time_to_escape = 10

# The maximum distance (squared) between bubble centers in two

# different timesteps for them to be considered the same bubble.

same_distance = 16

# Directory containing VTK files holding atoms.

#data_path = '/Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm/vtk-files'

data_path = '/Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files'

```

%% Cell type:code id:21fe4234-a349-402d-afdf-c9cb9a04f36a tags:

``` python

import sys

sys.path.append('/Users/4d5/builds/VTK/Release/lib/python3.11/site-packages')

import vtk

import numpy as np

from vtk.util.numpy_support import vtk_to_numpy

```

%% Cell type:code id:ed2a82cd-2c92-4e73-a56b-06fd2ffd931b tags:

``` python

class Bubble:

    '''Holds data for a bubble of atoms.

Bubble is constructed with a vtkDataSet that contains the atoms

of the bubble. Alternately, it can be constructed with a set of

classified atoms and the identifier of the atom in the `RegionId`

field (as comes from find_bubbles).'''

    def __init__(self, data, bubble_id = 0):

        if (bubble_id < 1):

            self.data = data

        else:

            threshold = vtk.vtkThreshold()

            threshold.SetInputData(data)

            threshold.SetInputArrayToProcess(

                0, 0, 0,

                'vtkDataObject::FIELD_ASSOCIATION_POINTS', 'RegionId')

            threshold.SetLowerThreshold(bubble_id - 0.1)

            threshold.SetUpperThreshold(bubble_id + 0.1)

            threshold.Update()

            self.data = threshold.GetOutput()

        self._center = None

    def get_number_of_atoms(self):

        return self.data.GetNumberOfPoints()

    def get_center(self):

        if not self._center:

            compute_center = vtk.vtkCenterOfMass()

            compute_center.SetInputData(self.data)

            compute_center.Update()

            self._center = compute_center.GetCenter()

        return self._center

```

%% Cell type:code id:122aa890-fdfe-4f24-9726-ec4fa0f9e8f6 tags:

``` python

class HeAtoms:

    '''Holds information about one timestep of He atoms.

This class is initialized by providing a VTK dataset representing

a LAMMPS atom file containing an `id` field giving unique identifiers

to atoms and a `type` field identifying the type of each atom

(`type` == 2 is He atoms). The data needs no cell data. The constructor

must also be given a timestep.

Use the `ReadVTKFile()` method to load the data from a VTK data set.

The `bubbles` field contains a list of `Bubble` objects holding the

extracted bubbles. The `timestep` field contains the time of the

atoms. The `filename` field provides the file the atoms came from,

when available.'''

    def ReadVTKFile(filename):

        import re

        timestep_match = re.search(r'([0-9]+)\.vtk', filename)

        timestep = int(timestep_match.group(1))

        reader = vtk.vtkDataSetReader()

        reader.SetFileName(filename)

        reader.Update()

        return HeAtoms(reader.GetOutput(), timestep, filename)

    def __init__(self, data, timestep, filename=None):

        self.timestep = timestep

        self.filename = filename

        # Many filters expect the data to have cells. Make one cell per point.

        tocloud = vtk.vtkConvertToPointCloud()

        tocloud.SetInputData(data)

        tocloud.SetCellGenerationMode(tocloud.VERTEX_CELLS)

        # Extract the He atoms.

        extract_he = vtk.vtkThreshold()

        extract_he.SetInputConnection(tocloud.GetOutputPort())

        extract_he.SetInputArrayToProcess(

            0, 0, 0, 'vtkDataObject::FIELD_ASSOCIATION_POINTS', 'type')

        extract_he.SetLowerThreshold(1.9)

        extract_he.SetUpperThreshold(2.1)

        # Create a density field of atoms by splatting the atoms into a grid.

        # Here we use `vtkGaussianSplatter`. If this is taking too long, it

        # is possible to speed things up using a `vtkFastSplatter`.

        bounds = data.GetBounds()

        sampling_resolution = (

            int((bounds[1] - bounds[0]) / (lattice_parameter / 2)),

            int((bounds[3] - bounds[2]) / (lattice_parameter / 2)),

            int((bounds[5] - bounds[4]) / (lattice_parameter / 2)))

        splatter = vtk.vtkGaussianSplatter()

        splatter.SetInputConnection(extract_he.GetOutputPort())

        splatter.SetSampleDimensions(sampling_resolution)

        splatter.SetRadius(He_radius)

        splatter.SetNormalWarping(0)

        splatter.SetAccumulationModeToSum()

        # Use an image connectivity filter to find regions with nearby atoms.

        # The filter will assign each connected region a unique ID and write

        # that to point data.

        connectivity = vtk.vtkImageConnectivityFilter()

        connectivity.SetInputConnection(splatter.GetOutputPort())

        connectivity.SetLabelScalarTypeToInt()

        connectivity.SetLabelModeToSizeRank()

        #connectivity.SetLabelModeToSizeRank()

        # Lowering threshold to include regions to prevent missing atoms at edge

        # of bubble.

        connectivity.SetScalarRange(0.1, 100000000000)

        # Frustratingly, there appears to be a bug in `vtkImageConnectivityFilter`

        # that does not properly pass the origin and spacing of the image, which

        # means that it will no longer align with the atom data. Restore it to

        # the proper values.

        splatter.Update()

        splatter_output = splatter.GetOutput()

        reset_position = vtk.vtkImageChangeInformation()

        reset_position.SetInputConnection(connectivity.GetOutputPort())

        reset_position.SetOutputOrigin(splatter_output.GetOrigin())

        reset_position.SetOutputSpacing(splatter_output.GetSpacing())

        # Sample the connectivity field back onto the atoms that generated them.

        # This will add a `RegionId` field to the atoms that identifes which

        # cluster each one belongs to.

        probe = vtk.vtkProbeFilter()

        probe.SetInputConnection(0, extract_he.GetOutputPort())

        probe.SetInputConnection(1, reset_position.GetOutputPort())

        probe.CategoricalDataOn()

        probe.PassPointArraysOn()

        connectivity.Update()

        bubble_sizes = connectivity.GetExtractedRegionSizes()

        probe.Update()

        self.atoms = probe.GetOutput()

        self.bubbles = []

        for bubble_index in range(1, bubble_sizes.GetNumberOfValues()):

            if bubble_sizes.GetValue(bubble_index) < burst_size:

                break

            self.bubbles.append(Bubble(self.atoms, bubble_index))

        for bubble_index in range(1, int(self.atoms.GetPointData().

                                         GetArray('RegionId').GetRange()[1])):

            bubble = Bubble(self.atoms, bubble_index)

            if bubble.get_number_of_atoms() >= burst_size:

                self.bubbles.append(bubble)

```

%% Cell type:code id:034539e1-5b21-4952-ae11-eca00cd35283 tags:

``` python

def find_burst_bubbles(reference_atoms, future_atoms):

    '''Given an `HeAtoms` object at a reference time and another such object

at a future timestep, determine which bubbles in the reference time

have burst. When a bubble bursts, the atoms become evacuated from their

cavity into a vacuum and disappear from the simulation. When a large

number of atoms in a bubble are not in the next timestep, that is

indicitive of the bubble bursting.

This function returns a list of `Bubble` objects that have been detected

as burst.'''

    burst_bubbles = []

    future_ids = \

        vtk_to_numpy(future_atoms.atoms.GetPointData().GetArray('id'))

    for bubble in reference_atoms.bubbles:

        reference_ids = \

            vtk_to_numpy(bubble.data.GetPointData().GetArray('id'))

        num_gone = np.isin(reference_ids, future_ids,

                           assume_unique=True, invert=True).sum()

        fraction = float(num_gone)/len(reference_ids)

        if (num_gone >= burst_size) and (fraction >= burst_fraction):

            print(num_gone, fraction)

            #print(num_gone, fraction)

            burst_bubbles.append(bubble)

    return burst_bubbles

```

%% Cell type:code id:0f551511-82e4-4e98-8557-6315b7beacb3 tags:

``` python

# Iterate over all the files in the data path and find potential bubble bursts.

import queue

bubble_queue = queue.Queue(time_to_escape - 1)

print('begin_time,end_time,center_x,center_y,center_z,size')

active_bursts = []

import glob

for filename in sorted(glob.glob(data_path + '/*-0890*.vtk')):

for filename in sorted(glob.glob(data_path + '/*.vtk')):

    next_atoms = HeAtoms.ReadVTKFile(filename)

    if bubble_queue.full():

        reference_atoms = bubble_queue.get_nowait()

        burst_bubbles = find_burst_bubbles(reference_atoms, next_atoms)

        for bubble in burst_bubbles:

            print(reference_atoms.timestep, bubble.get_center(), filename)

        for bubble in find_burst_bubbles(reference_atoms, next_atoms):

            is_active = False

            for active_burst in active_bursts:

                bubble_center = bubble.get_center()

                active_center = active_burst['bubble'].get_center()

                distsquare = (

                    (bubble_center[0]-active_center[0])**2 +

                    (bubble_center[1]-active_center[1])**2 +

                    (bubble_center[2]-active_center[2])**2

                )

                if distsquare < same_distance:

                    active_burst['endtime'] = next_atoms.timestep

                    is_active = True

            if not is_active:

                active_bursts.append({ 'bubble': bubble,

                                       'begintime' : reference_atoms.timestep,

                                       'endtime' : next_atoms.timestep })

    new_active = []

    for active_burst in active_bursts:

        if active_burst['endtime'] < next_atoms.timestep:

            center = active_burst['bubble'].get_center()

            print('{},{},{:0.1f},{:0.1f},{:0.1f},{}'.format(

                active_burst['begintime'],

                active_burst['endtime'],

                center[0], center[1], center[2],

                active_burst['bubble'].data.GetNumberOfPoints()))

        else:

            new_active.append(active_burst)

    active_bursts = new_active

    bubble_queue.put_nowait(next_atoms)

```

%% Output

    129 0.8164556962025317

    89018 (219.2584626403036, 128.282757288293, 4.523309201847974) /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-089027.vtk

    133 0.8417721518987342

    89019 (219.4156682461123, 128.26865285559546, 4.131269966311092) /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-089028.vtk

    134 0.8375

    89020 (219.19135427474976, 128.24084725379944, 3.909063812252134) /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-089029.vtk

    begin_time,end_time,center_x,center_y,center_x,size

    84241,84258,96.9,88.0,2.1,109

    85293,85302,170.9,113.4,2.6,29

    86015,86029,88.6,176.7,6.9,32

    86262,86271,47.9,2.6,3.6,29

    86264,86274,47.9,2.7,3.8,28

    86261,86277,47.7,244.9,2.2,32

    86267,86278,48.0,2.4,3.6,29

    86309,86325,61.8,217.1,1.0,56

    86537,86553,120.5,94.7,0.3,42

    86546,86561,188.8,132.3,3.5,60

    87201,87210,223.5,53.3,5.4,31

    87203,87214,223.4,53.2,5.4,30

    87418,87430,39.9,122.3,3.0,36

    87424,87433,40.1,122.9,3.3,27

    87473,87488,106.1,163.6,1.6,78

    89018,89029,219.3,128.3,4.5,158

    90220,90229,204.4,72.4,2.2,30

    90222,90232,204.5,72.6,2.2,30

    90226,90237,204.4,72.8,2.0,29

    90509,90518,10.3,86.5,2.6,100

    90512,90521,10.7,86.7,2.5,100

    91315,91324,37.1,138.5,1.7,35

    91317,91326,37.3,138.8,2.0,37

    91319,91328,37.1,138.7,1.5,36

    91828,91837,68.2,154.2,17.1,1480

    92003,92020,23.9,223.1,-0.0,119

    93481,93498,153.7,227.3,2.5,35

    94431,94440,45.6,63.9,14.9,1170

    95034,95050,180.9,131.6,0.3,116

    95299,95314,4.8,181.9,2.1,43

    96158,96167,116.9,149.2,-1.2,32

    96712,96723,229.0,79.9,-0.7,26

    96716,96728,228.9,80.1,-0.5,26

    96768,96777,44.8,192.5,-0.0,33

    98265,98281,56.0,178.7,-4.3,37

    98322,98338,178.4,190.1,18.6,1385

    99089,99106,29.1,124.6,3.8,51

    100375,100385,55.0,102.5,2.7,52

%% Cell type:code id:9dc7cc17-3cf3-4afd-a024-6ac6f7da8fbb tags:

%% Cell type:code id:60b38426-7bc1-4739-a365-6eb921df59d7 tags:

``` python

```