class PBFSimulator

class deeptime.data.PBFSimulator(domain_size: ndarray, initial_positions: ndarray, interaction_distance: float, n_jobs=None, n_solver_iterations: int = 5, gravity: float = 10.0, epsilon: float = 10.0, timestep: float = 0.016, rest_density: float = 1.0, tensile_instability_distance: float = 0.2, tensile_instability_k: float = 0.1)

A position based fluids simulator for two-dimensional systems. [1]

Its underlying principle is by definition of a rest density \(\rho_0\), which the particles in the system try to reach by a smoothed particle hydrodynamics style simulation [2] [3]. Up to numerics the simulation is deterministic given a set of parameters.

Parameters:

• domain_size ((2,) ndarray) – A 1-dimensional ndarray with two elements describing the extension of the simulation box.

• initial_positions ((N, 2) ndarray) – A 2-dimensional ndarray describing the \(N\) particles’ initial positions. This also fixes the number of particles. It is preserved during runs of the simulation.

• interaction_distance (float) – Interaction distance for particles, influences the cell size for the neighbor list.

• n_jobs (int or None, default=None) – Number threads to use for simulation.

• n_solver_iterations (int, default=5) – The number of solver iterations for particle position updates. The more, the slower the simulation but also the more accurate.

• gravity (float, default=10) – Gravity parameter which acts on particles’ velocities in each time step as constant force.

• epsilon (float, default=10) – Damping parameter. The higher, the more damping.

• timestep (float, default=0.016) – The timestep used for propagation.

• rest_density (float, defaul=1.) – The rest density \(\rho_0\).

• tensile_instability_distance (float, default=0.2) – Parameter responsible for surface-tension-like effects.

• tensile_instability_k (float, default=0.1) – Also controls surface-tension effects.

Notes

Each particle has positions \(p_1,\ldots,p_n \in\mathbb{R}^d\) and velocity \(v_1,\ldots,v_n\in\mathbb{R^d}\). For the local density a standard SPH estimator is used

\[\rho_i = \sum_i m_i W(p_i - p_j, h), \]

where \(m_i\) refers to the particle’s mass (here \(m_i \equiv 1\) for simplicity), \(W\) is a convolution kernel with interaction distance \(h\), i.e., particles with distance larger than \(h\) get assigned a weight of zero in the convolution operation.

The whole system has an underlying family of density constraints

\[C_i(p_1,\ldots,p_n) = \frac{\rho_i}{\rho_0} - 1. \]

The simulation algorithm tries to fulfill this constraint by updating each particle’s position, i.e., it tries to find \(\Delta p_k\) so that

\[C_k(p_1,\ldots,p_k + \Delta p_k,\ldots,p_n) = 0. \]

This is done by performing Newton’s optimization method along the constraint’s gradient. Because of potentially vanishing gradients, the optimization is regularized / dampened with a parameter \(\varepsilon\). In this implementation, larger values of \(\varepsilon\) lead to more regularization.

When a particle has only very few neighbors, it might not be possible to reach the rest density. This can lead to particle clumping (i.e., areas of negative pressures). One remedy is an artificial internal pressure term, here referred to as tensile instability. It modifies the weighting inside the convolution kernel itself. The default parameters tensile_instability_distance and tensile_instability_k are the suggested ones in the PBF publication.

References

[1]

Miles Macklin and Matthias Müller. Position based fluids. ACM Transactions on Graphics (TOG), 32(4):1–12, 2013.

[2]

Robert A Gingold and Joseph J Monaghan. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Monthly notices of the royal astronomical society, 181(3):375–389, 1977.

[3]

Leon B Lucy. A numerical approach to the testing of the fission hypothesis. The astronomical journal, 82:1013–1024, 1977.

Attributes

domain_size	The size of the domain.
n_particles	Number of particles in the simulation.

Methods

make_animation(trajectories[, stride, mode])	Creates a matplotlib animation object consisting of either scatter plots of contour plots for plotting particles and densities, respectively.
run(n_steps[, drift])	Performs n_steps many simulation steps and returns a trajectory of recorded particle positions.
simulate_oscillatory_force(n_oscillations, ...)	Runs 2*n_oscillations*n_steps many simulation steps in total by applying the drift alternatingly in positive and negative direction for n_steps each, n_oscillation times.
transform_to_density(trajectory[, n_grid_x, ...])	Transforms a two-dimensional PBF particle trajectory to a trajectory of densities by performing KDEs.

make_animation(trajectories, stride=1, mode='scatter', **kw)

Creates a matplotlib animation object consisting of either scatter plots of contour plots for plotting particles and densities, respectively.

Parameters:

• trajectories (list of ndarray) – Input trajectories. Must all have same shape.

• stride (int, default=1) – Apply stride to frames so that fewer data is presented.

• mode (str, default="scatter") – Aside from “scatter” this also supports “contourf” for densities.

• **kw –

  Optional keyword arguments. The following are supported:

  • ”n_grid_x”, required when in mode “contourf” (see transform_to_density())

  • ”n_grid_y”, required when in mode “contourf” (see transform_to_density())

  • ”figsize” controls figure size, defaults to (n_traj * 8, 6)

  • ”ncols” controls number of columns, defaults to n_traj

  • ”nrows” controls number of rows, defaults to 1

  • ”titles” optional list of strings for titles

Returns:

animation – Matplotlib animation object.

Return type:

FuncAnimation

run(n_steps: int, drift: float = 0.0)

Performs n_steps many simulation steps and returns a trajectory of recorded particle positions.

Parameters:

• n_steps (int) – Number of steps to run.

• drift (float, default=0) – Drift that is added to particles’ velocity in x direction.

Returns:

trajectory – Output trajectory.

Return type:

(n_steps, 2) ndarray

simulate_oscillatory_force(n_oscillations, n_steps, drift=0.3)

Runs 2*n_oscillations*n_steps many simulation steps in total by applying the drift alternatingly in positive and negative direction for n_steps each, n_oscillation times.

Parameters:

• n_oscillations (int) – Number of oscillations.

• n_steps (int) – Number of steps per run with a certain drift setting. Is then run again with the negated drift.

• drift (float) – The magnitude of the drift force into x direction.

Returns:

trajectory – Output trajectory.

Return type:

(T, n) ndarray

transform_to_density(trajectory, n_grid_x=20, n_grid_y=10, n_jobs=None)

Transforms a two-dimensional PBF particle trajectory to a trajectory of densities by performing KDEs. Since we have the prior knowledge that no particles get lost on the way, the densities are normalized frame-wise.

Parameters:

• trajectory ((T, n, 2) ndarray) – The input trajectory for n particles.

• n_grid_x (int, default=20) – Number of evaluation points of simulation box in x direction.

• n_grid_y (int, default=10) – Number of evaluation points of simulation box in y direction.

• n_jobs (int or None, default=None) – Number of jobs to use when transforming to densities.

Returns:

trajectory – Output trajectory

Return type:

(T, n_grid_x * n_grid_y) ndarray

property domain_size

The size of the domain.

property n_particles

Number of particles in the simulation.