In :
%%javascript
MathJax.Hub.Config({
TeX: { equationNumbers: { autoNumber: "AMS" } }
});

In :
from IPython.display import HTML

HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();$('div.prompt').hide();
} else {
$('div.input').show();$('div.prompt').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle); </script> <form action="javascript:code_toggle()"><input type="submit" value="Code Toggle"></form>''')  Out: In : from IPython.display import HTML HTML(''' <a href="https://github.com/usnistgov/pfhub/raw/master/benchmarks/benchmark8.ipynb" download> <button type="submit">Download Notebook</button> </a> ''')  Out: # Benchmark Problem 8: Homogeneous Nucleation¶ In : from IPython.display import HTML HTML('''{% include jupyter_benchmark_table.html num="" revision=1 %}''')  Out: ##### Revision History See the full benchmark table and table key See the journal publication entitled "Phase Field Benchmark Problems for Nucleation" for more details about the benchmark problems. Furthermore, read the extended essay for a discussion about the need for benchmark problems. ## Overview¶ Phase field modeling of nucleation has a long history and is covered in a number of reviews (Gránásy 2002, Castro 2003, Simmons 2004, Gránásy 2007, Warren 2009, Heo 2014, Gránásy 2019). The problem formulation of crystallization of an ideal pure liquid cooled below its melting point starts with homogeneous nucleation, a process in which the internal fluctuations of the undercooled liquid lead to the formation of crystal-like seeds able to grow to macroscopic sizes. The nucleation can be assisted by the presence of surfaces (container walls, foreign particles, etc.), in which case the process is termed heterogeneous nucleation. We note that homogeneous nucleation is an idealized formulation, and it is unlikely homogeneous nucleation occurs, due to impurities present in experimental apparatus. However, creating benchmark problems for homogeneous nucleation is still needed since the focus in nucleation research often lies on the underlying homogeneous nucleation, as it is the basis for advancing theoretical approaches to the much more complex phenomena governing heterogeneous nucleation. This benchmark (benchmark 8) targets homogeneous nucleation while benchmark 9 targets heterogeneous nucleation. There are two main modeling approaches to introduce nuclei into a metastable system: the Langevin noise method (Kubo 1966, Gránásy 2002) and the explicit nucleation method (Simmons 2000, Shi 2019) . We focus on the explicit method, where the critical nucleus size and the nucleation energy are determined by the classical nucleation theory. We break down our consideration of homogeneous nucleation into three parts. Part (a) considers the simple case of single-seed homogeneous nucleation. We explore how the particle size influences the evolution of the nucleus when the thermodynamic driving force is close to the critical value where we can observe whether the particle grows, dissolves or remains stationary. In parts (b) and (c) we consider multiple-seed homogeneous nucleation, with nuclei appearing at fixed time$t=0$or at random times distributed uniformly, respectively, with the latter part illustrating a scenario with constant nucleation rate. The Benchmark Problem probes the differences between transformation kinetics of the two cases and summarizes them using the Avrami plots based on the Johnson-Mehl-Avrami-Kolmogorov (JMAK) theory (Johnson 1939, Avrami 1939, Avrami 1940, Avrami 1941, Kolmogorov 1937). ## Model Formulation¶ ### The phase field model in classical nucleation theory¶ For this benchmark problems we use the simplest possible phase field model with a single non-conserved phase field$\phi$, which describes an isothermal pure substance with one liquid ($\phi=0$) and one solid ($\phi=1$) phase. The free energy of this system is \begin{equation} \mathcal{F}(\phi)=\int \left[\frac{\epsilon^2}{2}(\nabla \phi)^2+wg(\phi)-\Delta f p(\phi) \right] d{V}, \label{eqn:Fphi} \end{equation} set$g(\phi)=\phi^2(1-\phi)^2$, a simple double well function with minima at$\phi=0$and$\phi=1$,$ w$controls the double-well barrier height,$\epsilon^2$is the gradient energy coefficient, let$p(\phi)=\phi^3(10-15\phi+6\phi^2)$, which ensures that$p(0)=p'(0)=p'(1)=0$and$p(1)=1$, and$\Delta f$is the driving force for solidification at the simulation temperature ($\Delta f$is positive below the melting point). The time evolution of$\phi$is given by the Allen-Cahn equation (Allen 1979) \begin{equation} \frac{\partial \phi}{\partial t}=-M\frac{\delta \mathcal{F}}{\delta \phi}=M\left[\epsilon^2\nabla^2\phi-wg'(\phi)+\Delta fp'(\phi)\right] \label{eqn:dphidt} \end{equation} where$M$is the mobility parameter. We will restrict the problem to two dimensions (2D). For a planar interface, one can show that the equilibrium ($\Delta f$=0) solid-liquid interface profile with the interface centered at$x=x_0$is given by \begin{equation} \phi(x)=\frac{1}{2}\left[1-\tanh\left(\frac{x-x_0}{\sqrt{2}\ell}\right)\right], \label{eqn:phir} \end{equation} where$x$is the perpendicular distance from the interface. We use this expression as the initial condition when we introduce a nucleus. We can also obtain the width$\ell$and the free energy of the interface,$\gamma$, as \begin{equation} \ell=\sqrt{\frac{\epsilon^2}{w}} \label{eqn:l} \end{equation} and \begin{equation} \gamma=\frac{\sqrt{\epsilon^2w}}{3\sqrt{2}}. \label{eqn:gamma} \end{equation} Choosing a characteristic length scale$\xi=\ell$and a characteristic time scale$\tau=1/(Mw)$, we can obtain a nondimensional form of the problem (the nondimensional quantities are denoted by tildes), \begin{equation} \tilde{\mathcal{F}}(\phi)=\int \left[\frac{1}{2}(\tilde{\nabla} \phi)^2+g(\phi)- \tilde{\Delta f} p(\phi) \right] d{\tilde{V}} \label{eqn:tFphi} \end{equation} and \begin{equation} \frac{\partial \phi}{\partial \tilde{t}}=\tilde{\nabla}^2\phi-g'(\phi)+\tilde{\Delta f}p'(\phi), \label{eqn:dphidtt} \end{equation} with$\tilde{\Delta f}=\Delta f/w$. ### The properties of the classical nucleus¶ The explicit nucleation method that we employ introduces nuclei into a metastable system, with their radii and nucleation barrier determined from classical nucleation theory. Classical nucleation theory views crystallite fluctuations appearing in the undercooled liquid as small spherical domains of the bulk crystalline phase bounded by a mathematically sharp solid-liquid interface. For a 2D system, the free energy of a circular solid particle of radius$r$is \begin{equation} \Delta G(r) = 2 \pi r \gamma - \pi r^2 \Delta f \label{eqn:deltaG} \end{equation} where$\Delta f$is the nucleation driving force that is used in our phase field model, and$\gamma$is the free energy of the interface. The free energy of the particle is a balance between the energy cost in forming the solid-liquid interface, and the free energy from the driving force which is released when the crystalline particle forms. Once the rate of change of free energy with respect to particle size becomes negative, the particle can grow. Taking the derivative of$\Delta G(r)$we get the rate of change of free energy with respect to radius as \begin{equation} \frac{d\Delta G}{dr}=2\pi \gamma - 2\pi r\Delta f. \label{eqn:ddeltaGdr} \end{equation} By setting$d\Delta G/dr$to zero, we obtain the critical radius$r^*$as \begin{equation} r^* = \frac{\gamma}{\Delta f}. \label{eqn:rc} \end{equation} The corresponding critical nucleation free energy is \begin{equation} \Delta G^* = \frac{\pi \gamma^2}{\Delta f}. \label{eqn:deltaGc} \end{equation} Using the same units as before, the nondimensional forms of these quantities are \begin{equation} \tilde{r^*} = \frac{1}{3\sqrt{2}}\frac{1}{\tilde{\Delta f}} \label{eqn:trc} \end{equation} and \begin{equation} \tilde{\Delta G^*} = \frac{\pi}{18}\frac{1}{\tilde{\Delta f}} \label{eqn:tdeltaGc}, \end{equation} respectively. ### Avrami plots¶ Avrami plots describe how solids transform from one phase to another at constant temperature. In particular, they can describe the kinetics of nucleation. The Avrami plots come from the JMAK theory (Johnson 1939, Avrami 1939, Avrami 1940, Avrami 1941, Kolmogorov 1937), which makes a number of assumptions and simplifications: (i) nucleation occurs randomly and homogeneously over the entire un-transformed portion of the material, (ii) the growth rate is constant and does not depend on the extent of transformation, (iii) the particles have convex shape with the same orientation, and (iv) the size of the system is infinite (both in space and time). The basis of the Avrami plots is the transformed fraction (solid fraction)$Y(t)$vs. time. According to the JMAK theory, \begin{equation} Y(t)=1-\exp (-Kt^n), \label{eqn:Xt} \end{equation} where$K$is a constant depending on the nucleation and growth rates,$n=d+1$for continuous nucleation and$n=d$if nucleation happens only at$t=0$, and$d$is the number of spatial dimensions. If we plot$\log(-\log(1-Y))$vs.$\log(t)$, then for the JMAK kinetics (Eq. \ref{eqn:Xt}) we get a straight line with slope$n$. ## Benchmark 8 Specification¶ From here on, we will use only the nondimensional forms of the phase field (Eqs. \ref{eqn:tFphi}, \ref{eqn:dphidtt}, and \ref{eqn:phir}) and nucleation (Eqs. \ref{eqn:trc} and \ref{eqn:tdeltaGc}) equations, but we will drop the tildes. ### Part (a): single seed¶ In this problem we examine the morphology change of the nucleus for different initial radii. We consider a 2D simulation domain of size$100 \times 100$units centered at$x=y=0$with Neumann boundary conditions. The driving force is set to$\Delta f = \sqrt{2}/30$, which corresponds to a critical radius of$r^*=5$(Eq. \ref{eqn:trc}). Next, we place a circular seed at the center of the domain. Incorporating a nucleus into a diffuse interface phase field model leads to a small offset from a classical, sharp-interface model of a nucleus. To account for this, we incorporate a diffuse-interface seed using the profile$\phi(r)$given by modifying Eq. \ref{eqn:phir} to 2D where$r=\sqrt{x^2+y^2}$, with a radius$r_0=r^*$. This seed is the diffuse-interface approximation of the classical sharp interface nucleus corresponding to the given$\Delta f$, and therefore it should be fairly close to (an unstable) equilibrium. Part (a) is then defined by the Allen-Cahn equation \begin{equation} \frac{\partial \phi}{\partial {t}}={\nabla}^2\phi-g'(\phi)+{\Delta f}p'(\phi) \label{eqn:dphidtt_2} \end{equation} for the phase field$\phi$, and the three different initial conditions: in the first case with the seed radius$r_0$corresponding exactly to the critical one (critical nucleus), and in the other two cases, slightly below and above the critical radius (subcritical and supercritical nuclei): \begin{equation} \phi(r)=\frac{1}{2}\left[1-\tanh\left( \frac{r-r_0}{\sqrt{2}} \right)\right], \label{eqn:phitr_2} \end{equation} where$r_0=\{r^*, 0.99r^*, 1.01r^*\}$. Figure 1 shows the computational domain with an initial seed of radius$r=r^*$. The time evolution of the system is then followed for times up to$t = 200$units, and the solid fraction$Y$and the total free energy$F$are plotted as functions of time. Finally, the problem also includes a convergence test with respect to mesh (spatial resolution). The closer the initial radius$r$is to the critical radius$r^*$, the more sensitive the numerical integration will be to round-off errors. Therefore, a convergence test is included for$r_0=1.01r^*$using successive runs each halving the average spatial mesh size$\Delta x$between consecutive runs. At least 3 runs are required to estimate the rate of convergence based on the$L^2$error with the finest mesh being used as a gold standard. The simulations should be run to$t=200$and the final phase field (at$t=200$) for each successive run should be submitted. To submit results to the PFHub website, see the submission guidelines below. Illustration of the$100 \times 100$computational domain in 2D. ## Part (b): multiple initial seeds¶ The second part of the homogeneous nucleation benchmark problem focuses on the kinetics of nucleation using Avrami plots and compares them to JMAK theory. To obtain reasonable statistics for these two parts, the simulation volume needs to be larger so it can encompass a larger number of smaller nuclei. We therefore increase the domain size to$500 \times 500$units of length, and the driving force to$\Delta f = 1/(3 \sqrt{2})$, which corresponds to a critical radius of$r^*=1$. Random initial positions$\mathbf{r}_i$of 25 supercritical seeds$i, i=1,\ldots,25$are generated with$r_0=2.2r^*$drawn from a uniform distribution on the$500\times500$domain. The distribution of the phase field$\phi$is the sum of the phase fields$\phi_i$with profiles \begin{equation} \phi_i(r)=\frac{1}{2}\left[1-\tanh\left( \frac{|\mathbf{r}-\mathbf{r}_i|-2.2r^*}{\sqrt{2}} \right)\right]. \label{eqn:phitr_3} \end{equation} from the different seeds$i$. Overlaps between different seeds are managed by setting$\phi=1$in all regions where the sums of$\phi_i >1$. After placing initial seeds and adjusting the initial phase field$\phi$so that$\phi\leq1$everywhere, the simulation is run up to total time$t=200$, at which time the whole domain is transformed to a solid ($\phi=1$). This is repeated four times, each time with different random seeds, and the total free energy, time evolution of the solid fraction$Y$, discrete particle count$N$(the number of disjoint regions with$\phi=1$) are plotted for the five simulations, and Avrami plots are generated from the five simulations; from the Avrami plot the exponent$n$is estimated and compared to the JMAK theory. Both Part (b) and Part (c) involve a random number generator. Because different executions will in general generate different random numbers, specific individual solutions will in general not be reproduced. However, statistics such as the average slope of Avrami plots can be meaningfully compared. To submit results to the PFHub website, see the submission guidelines below. ## Part (c): multiple seeds at random times¶ For the third part of the homogeneous benchmark problem, the computational domain is first expanded to$1000 \times1000$. Instead of inserting nuclei at fixed time$t = 0$as in part two of the problem, in this case 100 random nucleation times$t_i$are generated,$i=1,\ldots,100$, drawn from a uniform distribution in the interval$t_i \in [0, 600)$with centers$\mathbf{r}_i$drawn from a uniform distribution on the$1000\times1000$domain. We decrease the driving force to$\Delta f=1/(6\sqrt{2})$, which corresponds to a critical radius or$r^*=2$, but keep the initial radii of the inserted nuclei the same as those in the second part of the problem \begin{equation} \phi_i(r)=\frac{1}{2}\left[1-\tanh\left( \frac{|\mathbf{r}-\mathbf{r}_i|-1.1r^*}{\sqrt{2}} \right)\right]. \end{equation} Again,$\phi$is set to unity in regions of overlaps of nuclei. The simulations are run up to times$t=600$, and repeated for four more random initial conditions. The total free energy plot the total free energy, the time evolution of the solid fraction$Y$and the discrete particle count$N$, are plotted as functions of time, and an Avrami plot is generated for the five runs. To submit results to the PFHub website, see the submission guidelines below. ## Submission Guidelines¶ Please use the upload form to upload your results. In addition to that specified, further data to upload can include a Youtube video, images of the nuclei at different times, or the entire phase field variable at different times. This data can be uploaded directly to the website or stored at a secondary location and only the link to the data provided. This data is not required, but will help others view your work. ### Part (a)¶ Part (a) requires running three simulations with a single seed with different values of$r_0$using the model formulation defined above. Optionally, a further two simulations can be run to check convergence, but this is not required. The three required simulations and two optional simulations are 1.$r_0 = 0.99 r^*$,${\Delta x} = {\Delta x}_0$2.$r_0 = r^*$,${\Delta x} = {\Delta x}_0$3.$r_0 = 1.01 r^*$,${\Delta x} = {\Delta x}_0$4.$r_0 = 1.01 r^*$,${\Delta x} = {\Delta x}_1$(optional) 5.$r_0 = 1.01 r^*$,${\Delta x} = {\Delta x}_2$(optional) where${\Delta x}_0, {\Delta x}_1$and${\Delta x}_2$represent simulations at different mesh sizes.$\Delta x$refinement should be selected so that successively refined meshes have roughly$4\times$as many elements. When uploading the results, the$\Delta x$values should be provided in the datafile alongside the values for$\phi$at each location in the domain. Each simulation should be run to$t=200$and the following data should be collected. • The solid fraction in the domain (the integral$\int_V \phi dV / \int_V dV$) during the transient (not necessarily at every time step but for ~100 data points). • The free energy integrated over the whole domain,$\mathcal{F}$, during the transient. • The phase field,$\phi$, at each grid point in the domain at$t=200$. • The grid spacing,$\Delta x$at each grid point in the domain at$t=200$if submitting the optional simulations. One PFHub upload is required for all the simulations. The solid fraction and free energy data can be included in the same CSV file or in separate data files, but either way still require 2 entries in the "Data Files" section on the upload form. The entries should be labeled as solid_fraction_{i} or free_energy_{i} depending on the simulation number (1 to 5). The CSV (JSON or TSV are also fine) files should have a format similar to the following. time,fraction,energy 0.000000000000000000e+00,8.214174263100689974e-03,3.579809478015085311e+00 1.000000000000000056e-01,8.205027223394931180e-03,3.577696619091880859e+00 2.000000000000000111e-01,8.197361418448430651e-03,3.576379039721957032e+00 ... The names of the columns are arbitrary as long as they are given in "Name of x-axis column" and "Name of y-axis column" in the form. The phase field data should be named phase_field_{i} in the "Short name of data" entry and have the following format x,y,phi,delta_x 0.0,1.0,0.5,0.44 0.0,1.5,0.49,0.44 ... where x and y is a position in the domain and phi is the value at that position. delta_x is the nominal grid spacing in that region given by the coordinates. This may be constant across the domain for gridded data, but vary for adaptive meshing. delta_x is only a representative value of the grid spacing in that region. The "3D" radio button should be selected for the phase_field_{i} entry. The grid spacing data should be named grid_spacing_{i} and have exactly the same data as above. Again the same file can be used, but two "Data Files" entries are required for both phase field and grid spacing. The "3D" radio button should be selected for the grid_spacing_{i} entry. The grid_spacing_{i} entries are not required, but if included will show convergence data for simulations 3, 4 and 5. ### Part (b)¶ Part (b) requires running 5 simulations with multiple seeds randomly positioned at$t=0$, see the defintion above. The simulations are run to$t=200$. The following data should be collected. • The solid fraction in the domain (the integral$\int_V \phi dV / \int_V dV$) during the transient (not necessarily at every time step but for ~100 data points). • The free energy integrated over the whole domain,$\mathcal{F}$, during the transient. • The discrete particle count,$N$, during the transient. One PFHub upload is required with five data files, but with a total of 15 entries in the "Data Files" section. The solid fraction, free energy and particle count can be included in the same data file for each simulation. The entries should be labeled as solid_fraction_{i}, free_energy_{i} or particle_count_{i} depending on the simulation number (1 to 5). The CSV (JSON or TSV are also fine) files should have a format similar to the following. time,fraction,energy,particle_count 0.000000000000000000e+00,8.214174263100689974e-03,3.579809478015085311e+00,24 1.000000000000000056e-01,8.205027223394931180e-03,3.577696619091880859e+00,24 2.000000000000000111e-01,8.197361418448430651e-03,3.576379039721957032e+00,23 ... ### Part (c)¶ Part (c) requires running 5 simulations with seeds randomly appearing between$t=0$and$t=600$, see the definition above. The simulations are run to$t=600$. The following data should be collected. • The solid fraction in the domain (the integral$\int_V \phi dV / \int_V dV$) during the transient (not necessarily at every time step but for ~100 data points). • The free energy integrated over the whole domain,$\mathcal{F}$, during the transient. • The discrete particle count,$N\$, during the transient.

One PFHub upload is required with five data files, but with a total of 15 entries in the "Data Files" section. The solid fraction, free energy and particle count can be included in the same data file for each simulation. The entries should be labeled as solid_fraction_{i}, free_energy_{i} or particle_count_{i} depending on the simulation number (1 to 5). The CSV (JSON or TSV are also fine) files should have a format similar to the following.

time,fraction,energy,particle_count
0.000000000000000000e+00,8.214174263100689974e-03,3.579809478015085311e+00,0
1.000000000000000056e-01,8.205027223394931180e-03,3.577696619091880859e+00,1
2.000000000000000111e-01,8.197361418448430651e-03,3.576379039721957032e+00,1
...

Please use the upload form to upload your results.

## Results¶

Results from this benchmark problem are displayed on the simulation result page for different codes.

[^][^]Gránásy, L and and Börzsönyi, T and Pusztai, T. 2002. Nucleation and bulk crystallization in binary phase field theory. URL

[^]Castro, Mario. 2003. Phase-field approach to heterogeneous nucleation. URL

[^]Simmons, J P and Wen, Y and Shen, C and Wang, Y. 2004. Microstructural development involving nucleation and growth phenomena simulated with the Phase Field method. URL

[^]Gránásy, L and T and Pusztai, T and Saylor, D and Warren, J A. 2007. Phase field theory of heterogeneous crystal nucleation. URL

[^]Warren, J A and Pusztai, T and Gránásy, L K and Gránásy, L. 2009. Phase field approach to heterogeneous crystal nucleation in alloys. URL

[^]Heo, T W and Chen, L. 2014. Phase-field modeling of nucleation in solid-state phase transformations. URL

[^]Gránásy, L and Tóth, G and Warren, J A and Podmaniczky, F and Tegze, G and Rátkai, L and Pusztai, T. 2019. Phase-field modeling of crystal nucleation in undercooled liquids -- A review. URL

[^]Kubo, R. 1966. The fluctuation-dissipation theorem. URL

[^]Simmons, J P and Shen, C and Wang, Y. 2000. Phase field modeling of simultaneous nucleation and growth by explicitly incorporating nucleation events. URL

[^]Shi, R and Khairallah, S and Heo, T W and Rolchigo, M and McKeown, J T and Matthews, M J. 2019. Integrated simulation framework for additively manufactured Ti-6Al-4V: melt pool dynamics, microstructure, solid-state phase transformation, and microelastic response. URL

[^][^]Johnson, William A and Mehl, R F. 1939. Reaction kinetics in processes of nucleation and growth. URL

[^][^]Avrami, Melvin. 1939. Kinetics of phase change. I General theory. URL

[^][^]Avrami, Melvin. 1940. Kinetics of phase change. II transformation-time relations for random distribution of nuclei. URL

[^][^]Avrami, Melvin. 1941. Granulation, phase change, and microstructure kinetics of phase change. III. URL

[^][^]Kolmogorov, A N. 1937. On the Statistical Theory of the Crystallization of Metals. URL

[^]Allen, Samuel M and Cahn, John W. 1979. A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. URL

## Appendix¶

The following code is given as a resource for those implementing the benchmark to help count the particles in Part (b) and Part (c). One possibility as outlined here is to use scipy.ndimage.label. This returns a tuple with a labeled array of unique features as the first item and a count of features as the second item. In the case of the phase field the values will require thresholding to 0 or 1 before using scipy.ndimage.label.

import numpy as np
import scipy.ndimage
import matplotlib.pyplot as plt

def make_particles(shape, n_particles, max_radius, seed=99):
"""Make an array with random particles

Args:
shape: shape of array
n_particles: number of particles in arrau
max_radius: maximum possible radius of particles
seed: random seed

Returns: numpy array mask with particles
"""
np.random.seed(seed)

def xy_grid(shape):
"""Generate a grid of x, y values in a single array

The returned array has shape (len(shape),) + shape + (1,). The last axis is
for the operations with particles.
"""
linspace = lambda x: np.linspace(0, x - 1, x) + 0.5
return np.array(np.meshgrid(*map(linspace, shape), indexing="ij"))[..., None]

# first axis is dimensions
# second and third axes are shape of domain
# last axis is the number of particles.
centers = (np.random.random((len(shape), n_particles)) * np.array(shape)[:, None])[:, None, None]
radii = (np.random.random(n_particles) * max_radius)[None, None]
return numpy.any(((xy_grid(shape) - centers)**2).sum(0) < radii**2, axis=-1)

arr = make_particles((100, 200), 20, 10)
plt.imshow(arr)
print('Actual number of particles:', scipy.ndimage.label(arr))

In :
import numpy as np
import scipy.ndimage
import matplotlib.pyplot as plt

def make_particles(shape, n_particles, max_radius, seed=99):
"""Make an array with random particles

Args:
shape: shape of array
n_particles: number of particles in arrau
max_radius: maximum possible radius of particles
seed: random seed

Returns: numpy array mask with particles
"""
np.random.seed(seed)

def xy_grid(shape):
"""Generate a grid of x, y values in a single array

The returned array has shape (len(shape),) + shape + (1,). The last axis is
for the operations with particles.
"""
linspace = lambda x: np.linspace(0, x - 1, x) + 0.5
return np.array(np.meshgrid(*map(linspace, shape), indexing="ij"))[..., None]

# first axis is dimensions
# second and third axes are shape of domain
# last axis is the number of particles.
centers = (np.random.random((len(shape), n_particles)) * np.array(shape)[:, None])[:, None, None]
radii = (np.random.random(n_particles) * max_radius)[None, None]
return numpy.any(((xy_grid(shape) - centers)**2).sum(0) < radii**2, axis=-1)

arr = make_particles((100, 200), 20, 10)
plt.imshow(arr)
print('Actual number of particles:', scipy.ndimage.label(arr))

Actual number of particles: 14 Upload Benchmark Results