Visualizing Diffusion: A Step-by-Step Guide with Python

Diffusion is a fundamental process in nature, governing phenomena from the spread of pollutants in the air to the mixing of substances at the molecular level. This process is driven by the random motion of particles and can be mathematically described by Fick's laws of diffusion. Visualizing diffusion can provide deeper insights into these mechanisms. In this blog, we'll explore how to create an animation of diffusion using Python, leveraging libraries such as NumPy and Matplotlib.

The Concept of Diffusion

Diffusion is the net movement of particles from regions of higher concentration to regions of lower concentration, resulting from random thermal motion. This process continues until a uniform concentration is achieved. The rate of diffusion can be quantified by the diffusion coefficient, which depends on factors such as temperature, the medium in which diffusion occurs, and the nature of the diffusing particles.

Why Use Agent-Based Simulation Rather Than Fick's Laws?

Fick's laws of diffusion provide a macroscopic view of diffusion, describing how concentration changes over time and space in a continuous medium. They are powerful tools for predicting diffusion rates and profiles in many scenarios. However, agent-based simulations offer several advantages that make them more suitable for certain types of analyses:

Microscopic Detail:

Agent-Based Simulation: Models individual particles (agents) and their interactions, capturing the stochastic nature of diffusion at a microscopic level.

Fick's Laws: Provide a deterministic, averaged view of particle concentration, lacking the granularity of individual particle behavior.

Stochastic Effects:

Agent-Based Simulation: Naturally incorporates randomness and fluctuations in particle movement, which are especially important in systems with low particle numbers or where local interactions are critical.

Fick's Laws: Assume a continuous, smooth concentration field, which may not accurately reflect systems with significant stochastic variations.

Complex Boundaries and Heterogeneities:

Agent-Based Simulation: Easily accommodates complex geometries, boundaries, and heterogeneous environments, where diffusion properties may vary spatially.

Fick's Laws: Require complex boundary conditions and may not handle heterogeneities as intuitively or accurately.

Interacting Particles:

Agent-Based Simulation: Can simulate interactions between particles, such as reactions or binding events, providing a richer model for systems where particle interactions are significant.

Fick's Laws: Primarily focus on diffusion alone, requiring additional equations and complexity to model interactions.

Setting Up the Environment

To simulate and visualize diffusion, we need to set up our environment with the necessary Python libraries

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation        

We'll use NumPy for efficient numerical operations and Matplotlib for creating the animation.

The DiffusionAnimation Class

The core of our simulation is the DiffusionAnimation class, which handles the initialization and updating of the particle grid. This class encapsulates the entire diffusion process.

class DiffusionAnimation:
    def __init__(self, size, central_density, diffusion_rate=0.1):
        self.size = size
        self.central_density = central_density  # Initial particle density in the center
        self.diffusion_rate = diffusion_rate  # Probability of particle diffusion
        self.grid = np.zeros((size, size))
        self.initialize_central_density()        

Initializing the Central Density

In the initialize_central_density method, we define a central region in the grid and populate it with particles:

    def initialize_central_density(self):
        central_region_size = self.size // 4  # Size of the central region
        start = self.size // 2 - central_region_size // 2
        end = start + central_region_size
        # Fill the central region with particles
        self.grid[start:end, start:end] = self.central_density
        

Central Region Size Calculation: central_region_size = self.size // 4: Defines the size of the central region as one-fourth of the total grid size. Example: For a grid size of 100, the central region will be 25x25.

Determining Indices: start = self.size // 2 - central_region_size // 2: Calculates the starting index to center the region. end = start + central_region_size: Determines the ending index. For a 100x100 grid, with central_region_size of 25, start will be 38 and end will be 63.

Initializing Grid Values: self.grid[start:end, start:end] = self.central_density: Sets the particle density in the central region to the specified value (e.g., 0.8).

Updating the Grid

The update method simulates the diffusion process by allowing particles to move to adjacent cells based on a probabilistic model.

    def update(self, frame):
        new_grid = self.grid.copy()
        for i in range(self.size):
            for j in range(self.size):
                # Random diffusion to adjacent cells
                if self.grid[i, j] > 0:
                    for di, dj in [(0, 1), (0, -1), (1, 0), (-1, 0)]:
                        if np.random.rand() < self.diffusion_rate:
                            ni, nj = (i + di) % self.size, (j + dj) % self.size
                            new_grid[ni, nj] += self.grid[i, j] * self.diffusion_rate
                            new_grid[i, j] -= self.grid[i, j] * self.diffusion_rate
        self.grid[:] = new_grid
        self.mat.set_data(self.grid)        

new_grid = self.grid.copy(): Creates a copy to ensure updates are based on the original state. Nested loops iterate over each cell in the grid.

Diffusion Condition: if self.grid[i, j] > 0: Ensures only cells with particles attempt diffusion.

Adjacent Cells Diffusion: for di, dj in [(0, 1), (0, -1), (1, 0), (-1, 0)] Loops through adjacent cells (right, left, down, up). if np.random.rand() < self.diffusion_rate: Applies a probabilistic condition for diffusion. (i + di) % self.size, (j + dj) % self.size: Ensures indices wrap around grid boundaries. Ie it doesnt move beyond the boundaries , but wrap back to the center.

Animating the Diffusion

The animate method sets up and runs the animation, providing a visual representation of the diffusion process.

    def animate(self):
        fig, ax = plt.subplots()
        self.mat = ax.matshow(self.grid, cmap='hot', vmin=0, vmax=1)
        plt.colorbar(self.mat, label='Particle Density')
        ani = FuncAnimation(fig, self.update, frames=100, interval=100)
        return ani        

Figure and Axes Creation: fig, ax = plt.subplots(): Initializes the plotting area.

Initial Grid Display: self.mat = ax.matshow(self.grid, cmap='hot', vmin=0, vmax=1): Displays the initial grid as a heatmap.

Colorbar Addition: plt.colorbar(self.mat, label='Particle Density'): Adds a colorbar to represent particle densities.

Animation Setup: ani = FuncAnimation(fig, self.update, frames=100, interval=100): Creates the animation by repeatedly calling the update method. frames=100: Specifies the number of frames. interval=100: Sets the interval between frames (100 ms). Returning the Animation: return ani: Returns the animation object.

Running the Simulation

To visualize the diffusion process, create an instance of the DiffusionAnimation class and call the animate method:

# Usage
size = 100
central_density = 0.8  # High density of particles in the center
diffusion_rate = 0.25  # Probability of particle diffusion to an adjacent cell
animation = DiffusionAnimation(size, central_density, diffusion_rate)
ani = animation.animate()
plt.show()        

Results

The attached figures represent a two-dimensional simulation of particle diffusion, generated using a discrete agent-based. This simulation visualizes the evolution of particle density over a grid as particles diffuse from a central high-density region to the surrounding low-density areas.

Initial State

In the first image, we observe the initial state of the simulation:Grid Dimensions: The grid is a 100x100 matrix, initialized to contain particle densities ranging from 0 to 1. Central Density Initialization: The central 25x25 region (from indices 38 to 62) is initialized with a uniform particle density of 0.8. This setup creates a high-density region in the middle of the grid, which will serve as the starting point for the diffusion process. Color Map: The color map used is 'hot', where colors range from black (density = 0) to red, orange, and yellow (density = 1). The color bar on the right side of the figure indicates the particle density scale.

Evolved State

In the second image, we observe the state of the grid after several iterations of the diffusion process:

Diffusion Process: The diffusion is modeled as a random walk, where each particle has a certain probability (diffusion_rate = 0.25) of moving to an adjacent cell in each iteration. This probabilistic movement is governed by a diffusion coefficient, reflecting the random thermal motion of particles. Resultant Density Distribution: Over time, particles spread from the high-density central region to adjacent cells, leading to a more even distribution. The density in the central region decreases, while surrounding areas see an increase in particle density. Visual Changes: The central high-density region becomes more dispersed, and the overall grid shows a gradient from higher density at the center to lower density at the edges. This indicates the particles' tendency to move from high to low concentration regions, consistent with Fick's first law of diffusion.

Conclusions

The agent-based approach to simulating diffusion provides a rich, detailed, and intuitive way to explore the dynamics of particle movement. It complements the macroscopic view provided by Fick's laws, offering insights into the microscopic and stochastic nature of diffusion. By running the animation, you can visually observe how particles diffuse from a high-density central region to the rest of the grid, providing a dynamic and intuitive understanding of the diffusion process.