Coagulation Basic 4: Particle Resolved¶

Introduction

In aerosol science, understanding particle-particle interactions is crucial for predicting the evolution of particle size distributions. One such interaction is coagulation, where two particles collide and merge into a larger particle. Accurately modeling coagulation at the level of individual particles is known as the particle-resolved method.

The particle-resolved method tracks each particle individually, considering its unique properties and interactions. This method provides the most detailed representation of aerosol dynamics, making it ideal for cases where precision is paramount, such as in cloud microphysics or laboratory-scale studies.

However, this approach is computationally intensive because it requires simulating every individual particle and its interactions. Unlike the super droplet method, which uses statistical representations to reduce computational load, the direct particle-resolved method does not aggregate particles into larger groups. Instead, every particle is treated independently, ensuring that every interaction is explicitly modeled.

This notebook provides a step-by-step guide to simulating coagulation using a pure particle-resolved approach, demonstrating how individual particles evolve over time without any simplifications or approximations in particle grouping.

Setup and Imports

We'll start by importing the necessary libraries and setting up the environment.

In [47]:

# %% particle resolved coagulation example
import numpy as np  # For numerical operations and array manipulations
import matplotlib.pyplot as plt  # For plotting graphs and visualizations

# Import specific modules from the particula package
from particula.dynamics.coagulation import brownian_kernel
from particula.dynamics.coagulation import particle_resolved_method

from particula.particles.properties.lognormal_size_distribution import (
    lognormal_sample_distribution,
)

# %% particle resolved coagulation example import numpy as np # For numerical operations and array manipulations import matplotlib.pyplot as plt # For plotting graphs and visualizations # Import specific modules from the particula package from particula.dynamics.coagulation import brownian_kernel from particula.dynamics.coagulation import particle_resolved_method from particula.particles.properties.lognormal_size_distribution import ( lognormal_sample_distribution, )

Generating Distribution¶

In this section, we generate a sample particle size distribution following a lognormal distribution. The lognormal distribution is commonly used in aerosol science to describe particle size distributions.

Coagulation Kernel

We also calculate the Brownian coagulation kernel for these particles, which quantifies the probability of coagulation between particles of different sizes.

Random seed

We set a random seed to ensure reproducibility of the results.

In [48]:

# lognormal spacing

radius_bins = np.logspace(
    -9, -6, num=20
)  # Define the radius bins for the distribution
mass_bins = (
    4 / 3 * np.pi * radius_bins**3 * 1000
)  # Calculate the mass of the particles in the bins

kernel = brownian_kernel.brownian_coagulation_kernel_via_system_state(
    radius_particle=radius_bins,
    mass_particle=mass_bins,
    temperature=298.15,
    pressure=101325,
)  # Calculate the Brownian coagulation kernel for the radius bins

random_generator = np.random.default_rng(12345)

# lognormal spacing radius_bins = np.logspace( -9, -6, num=20 ) # Define the radius bins for the distribution mass_bins = ( 4 / 3 * np.pi * radius_bins**3 * 1000 ) # Calculate the mass of the particles in the bins kernel = brownian_kernel.brownian_coagulation_kernel_via_system_state( radius_particle=radius_bins, mass_particle=mass_bins, temperature=298.15, pressure=101325, ) # Calculate the Brownian coagulation kernel for the radius bins random_generator = np.random.default_rng(12345)

Sampling the Particle Distribution

We then sample particles from the lognormal distribution. These particles will be sorted by size to prepare for the coagulation step.

In [49]:

# %% sample particle distribution
particle_radius = lognormal_sample_distribution(
    mode=np.array([1e-8, 1e-7]),
    geometric_standard_deviation=np.array([1.4, 1.4]),
    number_of_particles=np.array([5000, 1000]),
    number_of_samples=10_000,
)
# particle_radius = np.sort(particle_radius)
particles_original = particle_radius.copy()

# %% sample particle distribution particle_radius = lognormal_sample_distribution( mode=np.array([1e-8, 1e-7]), geometric_standard_deviation=np.array([1.4, 1.4]), number_of_particles=np.array([5000, 1000]), number_of_samples=10_000, ) # particle_radius = np.sort(particle_radius) particles_original = particle_radius.copy()

Coagulation Step¶

In the coagulation step, particles collide and merge over a given time step. The super droplet method efficiently simulates this process by adjusting the particle sizes and concentrations based on the calculated kernel and the specified volume and time step.

In [50]:

# %% Coagulation step

delta_t = 100  # time step in seconds
total_number_concentration = 1_000_000 * 1e6  # particles per m^3
total_number_tracked = len(particle_radius)
volume_sim = total_number_tracked / total_number_concentration

loss_gain_index = particle_resolved_method.particle_resolved_coagulation_step(
    particle_radius=particle_radius,
    kernel=kernel,
    kernel_radius=radius_bins,
    volume=volume_sim,
    time_step=delta_t,
    random_generator=random_generator,
)
particle_radius, gain, loss = particle_resolved_method.particle_resolved_update_step(
    particle_radius=particle_radius,
    gain=np.zeros_like(particle_radius),
    loss=np.zeros_like(particle_radius),
    small_index=loss_gain_index[:, 0],
    large_index=loss_gain_index[:, 1],
)

print(loss_gain_index.shape)
print(loss_gain_index)

# %% Coagulation step delta_t = 100 # time step in seconds total_number_concentration = 1_000_000 * 1e6 # particles per m^3 total_number_tracked = len(particle_radius) volume_sim = total_number_tracked / total_number_concentration loss_gain_index = particle_resolved_method.particle_resolved_coagulation_step( particle_radius=particle_radius, kernel=kernel, kernel_radius=radius_bins, volume=volume_sim, time_step=delta_t, random_generator=random_generator, ) particle_radius, gain, loss = particle_resolved_method.particle_resolved_update_step( particle_radius=particle_radius, gain=np.zeros_like(particle_radius), loss=np.zeros_like(particle_radius), small_index=loss_gain_index[:, 0], large_index=loss_gain_index[:, 1], ) print(loss_gain_index.shape) print(loss_gain_index)

(3039, 2)
[[ 759 9538]
 [3153 9522]
 [5383 9599]
 ...
 [9152 9417]
 [9173 9734]
 [8997 9213]]

Plotting the New Distribution

Finally, we plot the particle size distribution before and after coagulation. This visualization helps us understand the effect of the coagulation process on the particle size distribution.

In [51]:

# %% plot new distribution
fig, ax = plt.subplots()
ax.hist(
    particles_original, bins=100, histtype="step", color="black", density=True
)
ax.hist(
    particle_radius[particle_radius>0], bins=100, histtype="step", color="blue", density=True
)
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("Particle radius (m)")
ax.set_ylabel("Frequency")
plt.show()

# %% plot new distribution fig, ax = plt.subplots() ax.hist( particles_original, bins=100, histtype="step", color="black", density=True ) ax.hist( particle_radius[particle_radius>0], bins=100, histtype="step", color="blue", density=True ) ax.set_xscale("log") ax.set_yscale("log") ax.set_xlabel("Particle radius (m)") ax.set_ylabel("Frequency") plt.show()

Plotting the Loss and Gain of Particles

We also plot the loss and gain of particles due to coagulation. This visualization provides insights into the coagulation process and how it affects the particle population.

In [52]:

fig, ax = plt.subplots()
ax.hist(loss[loss>0], bins=100, histtype="step", color="red", density=True)
ax.hist(gain[gain>0], bins=100, histtype="step", color="green", density=True)
ax.set_xscale("log")
# ax.set_yscale("log")
ax.set_xlabel("Particle radius (m)")
ax.set_ylabel("Frequency")
plt.show()

fig, ax = plt.subplots() ax.hist(loss[loss>0], bins=100, histtype="step", color="red", density=True) ax.hist(gain[gain>0], bins=100, histtype="step", color="green", density=True) ax.set_xscale("log") # ax.set_yscale("log") ax.set_xlabel("Particle radius (m)") ax.set_ylabel("Frequency") plt.show()

Direct Time Stepping¶

With the first coagulation step completed, we can now proceed to the next time step. We repeat the coagulation process for the new particle distribution, updating the particle sizes and concentrations accordingly. This iterative approach allows us to simulate the evolution of the particle size distribution over time.

Here we use a simple for loop to perform multiple coagulation steps. In practice, more sophisticated time-stepping methods may be used to improve efficiency and accuracy.

In [53]:

# Initial distribution for the simulation
particles_0 = particles_original
particles_i = particles_0

# Define the time array for the simulation
time_array = np.linspace(
    start=0, stop=1000, num=100
)  # Time span of 1000 seconds
time_interval = (
    time_array[1] - time_array[0]
)  # Time interval between each step

# Array to store the distribution at each time step
particles_matrix = np.zeros([len(time_array), len(particles_0)])

# Simulate the coagulation process over time
for i, dpa in enumerate(time_array):
    if i > 0:

        loss_gain_index = particle_resolved_method.particle_resolved_coagulation_step(
            particle_radius=particles_i,
            kernel=kernel,
            kernel_radius=radius_bins,
            volume=volume_sim,
            time_step=time_interval,
            random_generator=random_generator,
        )
        particles_i, _, _ = particle_resolved_method.particle_resolved_update_step(
            particle_radius=particles_i,
            gain=np.zeros_like(particles_i),
            loss=np.zeros_like(particles_i),
            small_index=loss_gain_index[:, 0],
            large_index=loss_gain_index[:, 1],
        )

        # Ensure no negative concentrations (set to zero if less than zero)
        particles_i[particles_i < 0] = 0

    # Store the updated distribution for the current time step
    particles_matrix[i, :] = particles_i

# Initial distribution for the simulation particles_0 = particles_original particles_i = particles_0 # Define the time array for the simulation time_array = np.linspace( start=0, stop=1000, num=100 ) # Time span of 1000 seconds time_interval = ( time_array[1] - time_array[0] ) # Time interval between each step # Array to store the distribution at each time step particles_matrix = np.zeros([len(time_array), len(particles_0)]) # Simulate the coagulation process over time for i, dpa in enumerate(time_array): if i > 0: loss_gain_index = particle_resolved_method.particle_resolved_coagulation_step( particle_radius=particles_i, kernel=kernel, kernel_radius=radius_bins, volume=volume_sim, time_step=time_interval, random_generator=random_generator, ) particles_i, _, _ = particle_resolved_method.particle_resolved_update_step( particle_radius=particles_i, gain=np.zeros_like(particles_i), loss=np.zeros_like(particles_i), small_index=loss_gain_index[:, 0], large_index=loss_gain_index[:, 1], ) # Ensure no negative concentrations (set to zero if less than zero) particles_i[particles_i < 0] = 0 # Store the updated distribution for the current time step particles_matrix[i, :] = particles_i

Plotting the Final Distribution

Finally, we plot the final particle size distribution after multiple coagulation steps. This visualization shows how the particle size distribution evolves over time due to coagulation.

In [54]:

filtered = particles_matrix[-1, :]>0

fig, ax = plt.subplots()
ax.hist(
    particles_original, bins=100, histtype="step", color="black", density=True
)
ax.hist(
    particles_matrix[-1, filtered],
    bins=100,
    histtype="step",
    color="blue",
    density=True,
)
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("Particle radius (m)")
ax.set_ylabel("Frequency")
plt.show()

filtered = particles_matrix[-1, :]>0 fig, ax = plt.subplots() ax.hist( particles_original, bins=100, histtype="step", color="black", density=True ) ax.hist( particles_matrix[-1, filtered], bins=100, histtype="step", color="blue", density=True, ) ax.set_xscale("log") ax.set_yscale("log") ax.set_xlabel("Particle radius (m)") ax.set_ylabel("Frequency") plt.show()

In [55]:

# plot total number of particles
total_number = np.sum(particles_matrix>0, axis=1)

fig, ax = plt.subplots()
ax.plot(time_array, total_number)
ax.set_xlabel("Time (s)")
ax.set_ylabel("Number of particles")
plt.show()

# plot total number of particles total_number = np.sum(particles_matrix>0, axis=1) fig, ax = plt.subplots() ax.plot(time_array, total_number) ax.set_xlabel("Time (s)") ax.set_ylabel("Number of particles") plt.show()

In [56]:

# convert to concentration
total_concentration = total_number / volume_sim

fig, ax = plt.subplots()
ax.plot(time_array, total_concentration)
ax.set_xlabel("Time (s)")
ax.set_ylabel(r"Concentration $(m^{-3})$")
plt.show()

# convert to concentration total_concentration = total_number / volume_sim fig, ax = plt.subplots() ax.plot(time_array, total_concentration) ax.set_xlabel("Time (s)") ax.set_ylabel(r"Concentration $(m^{-3})$") plt.show()

Conclusion¶

This notebook demonstrated the use of the particle resolved coagulation method to simulate particle coagulation.

To be revise below¶

In [ ]:

"""coagulation rewirte"""

import numpy as np

from particula import particles
from particula.dynamics import coagulation
from particula.dynamics.coagulation import (
    particle_resolved_method,
    super_droplet_method,
)
from scipy.interpolate import RectBivariateSpline
from particula.util.input_handling import convert_units

# %%


particle_radius = particles.properties.lognormal_sample_distribution(
    mode=np.array([100e-9]),
    geometric_standard_deviation=np.array([1.3]),
    number_of_particles=np.array([1]),
    number_of_samples=100_000,
)

kernel_radius = np.logspace(-8, -6, 50)


# calc mass of particle from density 1000
mass_particles = 4 / 3 * np.pi * kernel_radius**3 * 1000
kernel = coagulation.brownian_coagulation_kernel_via_system_state(
    radius_particle=kernel_radius,
    mass_particle=mass_particles,
    temperature=300,
    pressure=101325,
)
# Step 4: Initialize a bivariate spline for interpolating kernel
# values between bin radii
interp_kernel = RectBivariateSpline(x=kernel_radius, y=kernel_radius, z=kernel)

random_generator = np.random.default_rng()

volume = 1 * convert_units("cm^3", "m^3")
time_step = 100
# %%

total_steps = 100
print(f"total time: {total_steps*time_step:.2f}")

gain_radius = np.zeros_like(particle_radius, dtype=np.float64)
loss_radius = np.zeros_like(particle_radius, dtype=np.float64)
for time_i in range(total_steps):
    number_in_bins, bin_indices = super_droplet_method.bin_particles(
        particle_radius, kernel_radius
    )
    # Step 3: Precompute unique bin pairs for efficient coagulation
    pair_indices = super_droplet_method.get_bin_pairs(bin_indices=bin_indices)

    small_index_total0 = np.array([], dtype=np.int64)
    large_index_total0 = np.array([], dtype=np.int64)

    # Iterate over each bin pair to calculate potential coagulation events
    for lower_bin, upper_bin in pair_indices:
        # get raidius indexes and filter out zeros
        small_indices = np.flatnonzero(
            (bin_indices == lower_bin) & (particle_radius > 0)
        )
        # filter small indices that are in small indices total
        small_indices = np.setdiff1d(small_indices, small_index_total0)

        large_indices = np.flatnonzero(
            (bin_indices == upper_bin) & (particle_radius > 0)
        )
        if np.size(small_indices) == 0 or np.size(large_indices) == 0:
            continue  # Skip to the next bin pair if no particles are present
        small_count = np.size(small_indices)
        large_count = np.size(large_indices)
        # print(f"small count: {small_count}, large count: {large_count}")

        # Retrieve the maximum kernel value for the current bin pair
        kernel_values = interp_kernel.ev(
            np.min(particle_radius[small_indices]),
            np.max(particle_radius[large_indices]),
        )
        # print(f"kernel max values: {kernel_values}")

        # Number of coagulation events
        if lower_bin == upper_bin:
            events = small_count * (large_count) / 2
        events = small_count * large_count

        # events = random_generator.poisson(events)
        tests = np.ceil(kernel_values * time_step * events / volume).astype(
            int
        )
        # print(f"events: {events}, tests: {tests}")

        if tests == 0:
            continue

        # Randomly select indices of particles involved in the coagulation
        small_replace = False if small_count > tests else True
        small_index = random_generator.choice(
            small_indices, tests, replace=small_replace
        )
        large_index = random_generator.choice(large_indices, tests)
        kernel_value = interp_kernel.ev(
            particle_radius[small_index], particle_radius[large_index]
        )
        # select diagonal
        if kernel_value.ndim > 1:
            kernel_value = np.diagonal(kernel_value)
        # print(f"kernel value: {kernel_value}")

        # Determine which coagulation events actually occur based on
        # interpolated kernel probabilities
        coagulation_probabilities = (
            kernel_value * time_step * events / (tests * volume)
        )
        # random number
        r = random_generator.uniform(size=tests)
        valid_indices = np.flatnonzero(r < coagulation_probabilities)
        # check if any valid indices are duplicate in small index
        # error of same small particle going to two different large particles
        _, unique_index = np.unique(
            small_index[valid_indices], return_index=True
        )
        small_index = small_index[valid_indices][unique_index]
        large_index = large_index[valid_indices][unique_index]

        # save the coagulation events
        small_index_total0 = np.append(small_index_total0, small_index)
        large_index_total0 = np.append(large_index_total0, large_index)

    commons, small_index_in_common, large_index_in_common = np.intersect1d(
        small_index_total0, large_index_total0, return_indices=True
    )
    # sort based on radius
    radius_argsort = np.argsort(particle_radius[commons])
    commons = commons[radius_argsort]
    small_index_in_common = small_index_in_common[radius_argsort]
    large_index_in_common = large_index_in_common[radius_argsort]

    # remap to largest particle in common
    for i, common in enumerate(commons):
        final_value = large_index_total0[small_index_in_common[i]]
        remap_index = np.flatnonzero(large_index_total0 == common)
        large_index_total0[remap_index] = final_value

    # update particle radius
    particle_radius, _, _ = (
        particle_resolved_method.particle_resolved_update_step(
            particle_radius=particle_radius,
            gain=gain_radius,
            loss=loss_radius,
            small_index=small_index_total0,
            large_index=large_index_total0,
        )
    )

print(f"concenration: {np.sum(number_in_bins) / volume:.4e}")

"""coagulation rewirte""" import numpy as np from particula import particles from particula.dynamics import coagulation from particula.dynamics.coagulation import ( particle_resolved_method, super_droplet_method, ) from scipy.interpolate import RectBivariateSpline from particula.util.input_handling import convert_units # %% particle_radius = particles.properties.lognormal_sample_distribution( mode=np.array([100e-9]), geometric_standard_deviation=np.array([1.3]), number_of_particles=np.array([1]), number_of_samples=100_000, ) kernel_radius = np.logspace(-8, -6, 50) # calc mass of particle from density 1000 mass_particles = 4 / 3 * np.pi * kernel_radius**3 * 1000 kernel = coagulation.brownian_coagulation_kernel_via_system_state( radius_particle=kernel_radius, mass_particle=mass_particles, temperature=300, pressure=101325, ) # Step 4: Initialize a bivariate spline for interpolating kernel # values between bin radii interp_kernel = RectBivariateSpline(x=kernel_radius, y=kernel_radius, z=kernel) random_generator = np.random.default_rng() volume = 1 * convert_units("cm^3", "m^3") time_step = 100 # %% total_steps = 100 print(f"total time: {total_steps*time_step:.2f}") gain_radius = np.zeros_like(particle_radius, dtype=np.float64) loss_radius = np.zeros_like(particle_radius, dtype=np.float64) for time_i in range(total_steps): number_in_bins, bin_indices = super_droplet_method.bin_particles( particle_radius, kernel_radius ) # Step 3: Precompute unique bin pairs for efficient coagulation pair_indices = super_droplet_method.get_bin_pairs(bin_indices=bin_indices) small_index_total0 = np.array([], dtype=np.int64) large_index_total0 = np.array([], dtype=np.int64) # Iterate over each bin pair to calculate potential coagulation events for lower_bin, upper_bin in pair_indices: # get raidius indexes and filter out zeros small_indices = np.flatnonzero( (bin_indices == lower_bin) & (particle_radius > 0) ) # filter small indices that are in small indices total small_indices = np.setdiff1d(small_indices, small_index_total0) large_indices = np.flatnonzero( (bin_indices == upper_bin) & (particle_radius > 0) ) if np.size(small_indices) == 0 or np.size(large_indices) == 0: continue # Skip to the next bin pair if no particles are present small_count = np.size(small_indices) large_count = np.size(large_indices) # print(f"small count: {small_count}, large count: {large_count}") # Retrieve the maximum kernel value for the current bin pair kernel_values = interp_kernel.ev( np.min(particle_radius[small_indices]), np.max(particle_radius[large_indices]), ) # print(f"kernel max values: {kernel_values}") # Number of coagulation events if lower_bin == upper_bin: events = small_count * (large_count) / 2 events = small_count * large_count # events = random_generator.poisson(events) tests = np.ceil(kernel_values * time_step * events / volume).astype( int ) # print(f"events: {events}, tests: {tests}") if tests == 0: continue # Randomly select indices of particles involved in the coagulation small_replace = False if small_count > tests else True small_index = random_generator.choice( small_indices, tests, replace=small_replace ) large_index = random_generator.choice(large_indices, tests) kernel_value = interp_kernel.ev( particle_radius[small_index], particle_radius[large_index] ) # select diagonal if kernel_value.ndim > 1: kernel_value = np.diagonal(kernel_value) # print(f"kernel value: {kernel_value}") # Determine which coagulation events actually occur based on # interpolated kernel probabilities coagulation_probabilities = ( kernel_value * time_step * events / (tests * volume) ) # random number r = random_generator.uniform(size=tests) valid_indices = np.flatnonzero(r < coagulation_probabilities) # check if any valid indices are duplicate in small index # error of same small particle going to two different large particles _, unique_index = np.unique( small_index[valid_indices], return_index=True ) small_index = small_index[valid_indices][unique_index] large_index = large_index[valid_indices][unique_index] # save the coagulation events small_index_total0 = np.append(small_index_total0, small_index) large_index_total0 = np.append(large_index_total0, large_index) commons, small_index_in_common, large_index_in_common = np.intersect1d( small_index_total0, large_index_total0, return_indices=True ) # sort based on radius radius_argsort = np.argsort(particle_radius[commons]) commons = commons[radius_argsort] small_index_in_common = small_index_in_common[radius_argsort] large_index_in_common = large_index_in_common[radius_argsort] # remap to largest particle in common for i, common in enumerate(commons): final_value = large_index_total0[small_index_in_common[i]] remap_index = np.flatnonzero(large_index_total0 == common) large_index_total0[remap_index] = final_value # update particle radius particle_radius, _, _ = ( particle_resolved_method.particle_resolved_update_step( particle_radius=particle_radius, gain=gain_radius, loss=loss_radius, small_index=small_index_total0, large_index=large_index_total0, ) ) print(f"concenration: {np.sum(number_in_bins) / volume:.4e}")