SPH equations¶

Basic SPH Equations¶

class pysph.sph.basic_equations.BodyForce(dest, sources, fx=0.0, fy=0.0, fz=0.0)[source]¶

Bases: pysph.sph.equation.Equation

Add a body force to the particles:

\(\boldsymbol{f} = f_x, f_y, f_z\)

Parameters:	fx (float) – Body force per unit mass along the x-axis fy (float) – Body force per unit mass along the y-axis fz (float) – Body force per unit mass along the z-axis

initialize(d_idx, d_au, d_av, d_aw)[source]¶

loop(d_idx, d_au, d_av, d_aw)[source]¶

class pysph.sph.basic_equations.ContinuityEquation(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Density rate:

\(\frac{d\rho_a}{dt} = \sum_b m_b \boldsymbol{v}_{ab}\cdot \nabla_a W_{ab}\)

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_arho)[source]¶

loop(d_idx, d_arho, s_idx, s_m, DWIJ, VIJ)[source]¶

class pysph.sph.basic_equations.IsothermalEOS(dest, sources, rho0, c0, p0)[source]¶

Bases: pysph.sph.equation.Equation

Compute the pressure using the Isothermal equation of state:

\(p = p_0 + c_0^2(\rho_0 - \rho)\)

Parameters:	rho0 (float) – Reference density of the fluid (\(\rho_0\)) c0 (float) – Maximum speed of sound expected in the system (\(c0\)) p0 (float) – Reference pressure in the system (\(p0\))

loop(d_idx, d_rho, d_p)[source]¶

class pysph.sph.basic_equations.MonaghanArtificialViscosity(dest, sources, alpha=1.0, beta=1.0)[source]¶

Bases: pysph.sph.equation.Equation

Classical Monaghan style artificial viscosity [Monaghan2005]

\[\begin{split}\frac{d\mathbf{v}_{a}}{dt}&=&-\sum_{b}m_{b}\Pi_{ab}\nabla_{a}W_{ab}\end{split}\]

where

\[\begin{split}\Pi_{ab}=\begin{cases}\frac{-\alpha_{\pi}\bar{c}_{ab}\phi_{ab}+ \beta_{\pi}\phi_{ab}^{2}}{\bar{\rho}_{ab}}, & \mathbf{v}_{ab}\cdot \mathbf{r}_{ab}<0\\0, & \mathbf{v}_{ab}\cdot\mathbf{r}_{ab}\geq0 \end{cases}\end{split}\]

with

\[\begin{split}\phi_{ab}=\frac{h\mathbf{v}_{ab}\cdot\mathbf{r}_{ab}} {|\mathbf{r}_{ab}|^{2}+\epsilon^{2}}\\\end{split}\]\[\begin{split}\bar{c}_{ab}&=&\frac{c_{a}+c_{b}}{2}\\\end{split}\]\[\begin{split}\bar{\rho}_{ab}&=&\frac{\rho_{a}+\rho_{b}}{2}\end{split}\]

References

[Monaghan2005]

J. Monaghan, “Smoothed particle hydrodynamics”, Reports on Progress in Physics, 68 (2005), pp. 1703-1759.

Parameters:	alpha (float) – produces a shear and bulk viscosity beta (float) – used to handle high Mach number shocks

initialize(d_idx, d_au, d_av, d_aw)[source]¶

loop(d_idx, s_idx, d_rho, d_cs, d_au, d_av, d_aw, s_m, s_rho, s_cs, VIJ, XIJ, HIJ, R2IJ, RHOIJ1, EPS, DWIJ, DT_ADAPT)[source]¶

class pysph.sph.basic_equations.SummationDensity(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Good old Summation density:

\(\rho_a = \sum_b m_b W_{ab}\)

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_rho)[source]¶

loop(d_idx, d_rho, s_idx, s_m, WIJ)[source]¶

class pysph.sph.basic_equations.VelocityGradient2D(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Compute the SPH evaluation for the velocity gradient tensor in 2D.

The expression for the velocity gradient is:

\(\frac{\partial v^i}{\partial x^j} = \sum_{b}\frac{m_b}{\rho_b}(v_b - v_a)\frac{\partial W_{ab}}{\partial x_a^j}\)

Notes

The tensor properties are stored in the variables v_ij where ‘i’ refers to the velocity component and ‘j’ refers to the spatial component. Thus v_21 is \(\frac{\partial v}{\partial x}\)

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_v00, d_v01, d_v10, d_v11)[source]¶

loop(d_idx, s_idx, s_m, s_rho, d_v00, d_v01, d_v10, d_v11, DWIJ, VIJ)[source]¶

class pysph.sph.basic_equations.XSPHCorrection(dest, sources, eps=0.5)[source]¶

Bases: pysph.sph.equation.Equation

Position stepping with XSPH correction [Monaghan1992]

\[\frac{d\mathbf{r}_{a}}{dt}=\mathbf{\hat{v}}_{a}=\mathbf{v}_{a}- \epsilon\sum_{b}m_{b}\frac{\mathbf{v}_{ab}}{\bar{\rho}_{ab}}W_{ab}\]

References

[Monaghan1992]

J. Monaghan, Smoothed Particle Hydrodynamics, “Annual Review of Astronomy and Astrophysics”, 30 (1992), pp. 543-574.

Parameters:	eps (float) – \(\epsilon\) as in the above equation

Notes

This equation must be used to advect the particles. XSPH can be turned off by setting the parameter eps = 0.

initialize(d_idx, d_ax, d_ay, d_az)[source]¶

loop(s_idx, d_idx, s_m, d_ax, d_ay, d_az, WIJ, RHOIJ1, VIJ)[source]¶

post_loop(d_idx, d_ax, d_ay, d_az, d_u, d_v, d_w)[source]¶

class pysph.sph.basic_equations.XSPHCorrectionForLeapFrog(dest, sources, eps=0.5)[source]¶

Bases: pysph.sph.equation.Equation

The XSPH correction [Monaghan1992] alone. This is meant to be used with a leap-frog integrator which already considers the velocity of the particles. It simply computes the correction term and adds that to ax, ay, az.

\[\frac{d\mathbf{r}_{a}}{dt}=\mathbf{\hat{v}}_{a}= - \epsilon\sum_{b}m_{b}\frac{\mathbf{v}_{ab}}{\bar{\rho}_{ab}}W_{ab}\]

References

[Monaghan1992]

J. Monaghan, Smoothed Particle Hydrodynamics, “Annual Review of Astronomy and Astrophysics”, 30 (1992), pp. 543-574.

Parameters:	eps (float) – \(\epsilon\) as in the above equation

Notes

This equation must be used to advect the particles. XSPH can be turned off by setting the parameter eps = 0.

initialize(d_idx, d_ax, d_ay, d_az)[source]¶

loop(s_idx, d_idx, s_m, d_ax, d_ay, d_az, WIJ, RHOIJ1, VIJ)[source]¶

Basic WCSPH Equations¶

class pysph.sph.wc.basic.ContinuityEquationDeltaSPH(dest, sources, c0, delta=0.1)[source]¶

Bases: pysph.sph.equation.Equation

Continuity equation with dissipative terms

\(\frac{d\rho_a}{dt} = \sum_b \rho_a \frac{m_b}{\rho_b} \left( \boldsymbol{v}_{ab}\cdot \nabla_a W_{ab} + \delta \eta_{ab} \cdot \nabla_{a} W_{ab} (h_{ab}\frac{c_{ab}}{\rho_a}(\rho_b - \rho_a)) \right)\)

References

[Marrone2011]

S. Marrone et al., “delta-SPH model for simulating violent impact flows”, Computer Methods in Applied Mechanics and Engineering, 200 (2011), pp 1526–1542.

Parameters:	c0 (float) – reference speed of sound delta (float) – coefficient used to control the intensity of diffusion of density

initialize(d_idx, d_arho)[source]¶

loop(d_idx, d_arho, s_idx, s_m, d_cs, s_cs, d_rho, s_rho, DWIJ, VIJ, XIJ, RIJ, HIJ, EPS)[source]¶

class pysph.sph.wc.basic.MomentumEquation(dest, sources, c0, alpha=1.0, beta=1.0, gx=0.0, gy=0.0, gz=0.0, tensile_correction=False)[source]¶

Bases: pysph.sph.equation.Equation

Classic Monaghan Style Momentum Equation with Artificial Viscosity

\[\frac{d\mathbf{v}_{a}}{dt}=-\sum_{b}m_{b}\left(\frac{p_{b}} {\rho_{b}^{2}}+\frac{p_{a}}{\rho_{a}^{2}}+\Pi_{ab}\right) \nabla_{a}W_{ab}\]

where

\[\begin{split}\Pi_{ab}=\begin{cases} \frac{-\alpha\bar{c}_{ab}\mu_{ab}+\beta\mu_{ab}^{2}}{\bar{\rho}_{ab}} & \mathbf{v}_{ab}\cdot\mathbf{r}_{ab}<0;\\ 0 & \mathbf{v}_{ab}\cdot\mathbf{r}_{ab}\geq0; \end{cases}\end{split}\]

with

\[\begin{split}\mu_{ab}=\frac{h\mathbf{v}_{ab}\cdot\mathbf{r}_{ab}} {\mathbf{r}_{ab}^{2}+\eta^{2}}\\\end{split}\]\[\begin{split}\bar{c}_{ab} = \frac{c_a + c_b}{2}\\\end{split}\]\[\bar{\rho}_{ab} = \frac{\rho_a + \rho_b}{2}\]

References

[Monaghan1992]

J. Monaghan, Smoothed Particle Hydrodynamics, “Annual Review of Astronomy and Astrophysics”, 30 (1992), pp. 543-574.

Parameters:	c0 (float) – reference speed of sound alpha (float) – produces a shear and bulk viscosity beta (float) – used to handle high Mach number shocks gx (float) – body force per unit mass along the x-axis gy (float) – body force per unit mass along the y-axis gz (float) – body force per unit mass along the z-axis tensilte_correction (bool) – switch for tensile instability correction (Default: False)

initialize(d_idx, d_au, d_av, d_aw)[source]¶

loop(d_idx, s_idx, d_rho, d_cs, d_p, d_au, d_av, d_aw, s_m, s_rho, s_cs, s_p, VIJ, XIJ, HIJ, R2IJ, RHOIJ1, EPS, DWIJ, DT_ADAPT, WIJ, WDP)[source]¶

post_loop(d_idx, d_au, d_av, d_aw, DT_ADAPT)[source]¶

class pysph.sph.wc.basic.MomentumEquationDeltaSPH(dest, sources, rho0, c0, alpha=1.0, gx=0.0, gy=0.0, gz=0.0)[source]¶

Bases: pysph.sph.equation.Equation

Momentum equation defined in JOSEPHINE and the delta-SPH model

\[\frac{du_{i}}{dt}=-\frac{1}{\rho_{i}}\sum_{j}\left(p_{j}+p_{i}\right) \nabla_{i}W_{ij}V_{j}+\mathbf{g}_{i}+\alpha hc_{0}\rho_{0}\sum_{j} \pi_{ij}\nabla_{i}W_{ij}V_{j}\]

where

\[\pi_{ij}=\frac{\mathbf{u}_{ij}\cdot\mathbf{r}_{ij}} {|\mathbf{r}_{ij}|^{2}}\]

References

[Marrone2011]

S. Marrone et al., “delta-SPH model for simulating violent impact flows”, Computer Methods in Applied Mechanics and Engineering, 200 (2011), pp 1526–1542.

[Cherfils2012]

J. M. Cherfils et al., “JOSEPHINE: A parallel SPH code for free-surface flows”, Computer Physics Communications, 183 (2012), pp 1468–1480.

Parameters:	rho0 (float) – reference density c0 (float) – reference speed of sound alpha (float) – coefficient used to control the intensity of the diffusion of velocity gx (float) – body force per unit mass along the x-axis gy (float) – body force per unit mass along the y-axis gz (float) – body force per unit mass along the z-axis

Notes

Artificial viscosity is used in this momentum equation and is controlled by the parameter \(\alpha\). This form of the artificial viscosity is similar but not identical to the Monaghan-style artificial viscosity.

initialize(d_idx, d_au, d_av, d_aw)[source]¶

loop(d_idx, s_idx, d_rho, d_cs, d_p, d_au, d_av, d_aw, s_m, s_rho, s_cs, s_p, VIJ, XIJ, HIJ, R2IJ, RHOIJ1, EPS, WIJ, DWIJ)[source]¶

post_loop(d_idx, d_au, d_av, d_aw, DT_ADAPT)[source]¶

class pysph.sph.wc.basic.PressureGradientUsingNumberDensity(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Pressure gradient discretized using number density:

\[\frac{d \boldsymbol{v}_a}{dt} = -\frac{1}{m_a}\sum_b (\frac{p_a}{V_a^2} + \frac{p_b}{V_b^2})\nabla_a W_{ab}\]

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_au, d_av, d_aw)[source]¶

loop(d_idx, s_idx, d_m, d_rho, s_rho, d_au, d_av, d_aw, d_p, s_p, d_V, s_V, DWIJ)[source]¶

class pysph.sph.wc.basic.TaitEOS(dest, sources, rho0, c0, gamma, p0=0.0)[source]¶

Bases: pysph.sph.equation.Equation

Tait equation of state for water-like fluids

\(p_a = \frac{c_{0}^2\rho_0}{\gamma}\left( \left(\frac{\rho_a}{\rho_0}\right)^{\gamma} -1\right)\)

References

[Cole1948]

H. R. Cole, “Underwater Explosions”, Princeton University Press, 1948.

[Batchelor2002]

G. Batchelor, “An Introduction to Fluid Dynamics”, Cambridge University Press, 2002.

[Monaghan2005]

J. Monaghan, “Smoothed particle hydrodynamics”, Reports on Progress in Physics, 68 (2005), pp. 1703-1759.

Parameters:	rho0 (float) – reference density of fluid particles c0 (float) – maximum speed of sound expected in the system gamma (float) – constant p0 (float) – reference pressure in the system (defaults to zero).

Notes

The reference speed of sound, c0, is to be taken approximately as 10 times the maximum expected velocity in the system. The particle sound speed is given by the usual expression:

\(c_a = \sqrt{\frac{\partial p}{\partial \rho}}\)

loop(d_idx, d_rho, d_p, d_cs)[source]¶

class pysph.sph.wc.basic.TaitEOSHGCorrection(dest, sources, rho0, c0, gamma)[source]¶

Bases: pysph.sph.equation.Equation

Tait Equation of State with Hughes and Graham Correction

\[p_a = \frac{c_{0}^2\rho_0}{\gamma}\left( \left(\frac{\rho_a}{\rho_0}\right)^{\gamma} -1\right)\]

where

\[\begin{split}\rho_{a}=\begin{cases}\rho_{a} & \rho_{a}\geq\rho_{0}\\ \rho_{0} & \rho_{a}<\rho_{0}\end{cases}`\end{split}\]

References

[Hughes2010]

J. P. Hughes and D. I. Graham, “Comparison of incompressible and weakly-compressible SPH models for free-surface water flows”, Journal of Hydraulic Research, 48 (2010), pp. 105-117.

Parameters:	rho0 (float) – reference density c0 (float) – reference speed of sound gamma (float) – constant

Notes

The correction is to be applied on boundary particles and imposes a minimum value of the density (rho0) which is set upon instantiation. This correction avoids particle sticking behaviour at walls.

loop(d_idx, d_rho, d_p, d_cs)[source]¶

class pysph.sph.wc.basic.UpdateSmoothingLengthFerrari(dest, sources, dim, hdx)[source]¶

Bases: pysph.sph.equation.Equation

Update the particle smoothing lengths

\(h_a = hdx \left(\frac{m_a}{\rho_a}\right)^{\frac{1}{d}}\)

References

[Ferrari2009]

A. Ferrari et al., “A new 3D parallel SPH scheme for free surface flows”, Computers and Fluids, 38 (2009), pp. 1203–1217.

Parameters:	dim (float) – number of dimensions hdx (float) – scaling factor

Notes

Ideally, the kernel scaling factor should be determined from the kernel used based on a linear stability analysis. The default value of (hdx=1) reduces to the formulation suggested by Ferrari et al. who used a Cubic Spline kernel.

Typically, a change in the smoothing length should mean the neighbors are re-computed which in PySPH means the NNPS must be updated. This equation should therefore be placed as the last equation so that after the final corrector stage, the smoothing lengths are updated and the new NNPS data structure is computed.

Note however that since this is to be used with incompressible flow equations, the density variations are small and hence the smoothing lengths should also not vary too much.

loop(d_idx, d_rho, d_h, d_m)[source]¶

Transport Velocity Formulation¶

References

class pysph.sph.wc.transport_velocity.ContinuityEquation(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Conservation of mass equation

Eq (6) in [Adami2012]:

\[\frac{d\rho_a}{dt} = \rho_a \sum_b \frac{m_b}{\rho_b} \boldsymbol{v}_{ab} \cdot \nabla_a W_{ab}\]

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_arho)[source]¶

loop(d_idx, s_idx, d_arho, s_m, d_rho, s_rho, VIJ, DWIJ)[source]¶

class pysph.sph.wc.transport_velocity.MomentumEquationArtificialStress(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Artificial stress contribution to the Momentum Equation

\[\frac{d\boldsymbol{v}_a}{dt} = \frac{1}{m_a}\sum_b (V_a^2 + V_b^2)\left[ \frac{1}{2}(\boldsymbol{A}_a + \boldsymbol{A}_b) : \nabla_a W_{ab}\right]\]

where the artificial stress terms are given by:

\[ \boldsymbol{A} = \rho \boldsymbol{v} (\tilde{\boldsymbol{v}} - \boldsymbol{v})\]

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_au, d_av, d_aw)[source]¶

loop(d_idx, s_idx, d_rho, d_u, d_v, d_w, d_V, d_uhat, d_vhat, d_what, d_au, d_av, d_aw, d_m, s_rho, s_u, s_v, s_w, s_V, s_uhat, s_vhat, s_what, DWIJ)[source]¶

class pysph.sph.wc.transport_velocity.MomentumEquationArtificialViscosity(dest, sources, c0, alpha=0.1)[source]¶

Bases: pysph.sph.equation.Equation

Artificial viscosity for the momentum equation

Eq. (11) in [Adami2012]:

\[\frac{d \boldsymbol{v}_a}{dt} = -\sum_b m_b \alpha h_{ab} c_{ab} \frac{\boldsymbol{v}_{ab}\cdot \boldsymbol{r}_{ab}}{\rho_{ab}\left(|r_{ab}|^2 + \epsilon \right)}\nabla_a W_{ab}\]

where

\[\begin{split}\rho_{ab} = \frac{\rho_a + \rho_b}{2}\\\end{split}\]\[\begin{split}c_{ab} = \frac{c_a + c_b}{2}\\\end{split}\]\[h_{ab} = \frac{h_a + h_b}{2}\]

Parameters:	alpha (float) – constant c0 (float) – speed of sound

initialize(d_idx, d_au, d_av, d_aw)[source]¶

loop(d_idx, s_idx, s_m, d_au, d_av, d_aw, RHOIJ1, R2IJ, EPS, DWIJ, VIJ, XIJ, HIJ)[source]¶

class pysph.sph.wc.transport_velocity.MomentumEquationPressureGradient(dest, sources, pb, gx=0.0, gy=0.0, gz=0.0, tdamp=0.0)[source]¶

Bases: pysph.sph.equation.Equation

Momentum equation for the Transport Velocity Formulation: Pressure

Eq. (8) in [Adami2013]:

\[\frac{d \boldsymbol{v}_a}{dt} = \frac{1}{m_a}\sum_b (V_a^2 + V_b^2)\left[-\bar{p}_{ab}\nabla_a W_{ab} \right]\]

where

\[\bar{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b}\]

Parameters:	pb (float) – background pressure gx (float) – Body force per unit mass along the x-axis gy (float) – Body force per unit mass along the y-axis gz (float) – Body force per unit mass along the z-axis tdamp (float) – damping time

Notes

This equation should have the destination as fluid and sources as fluid and boundary particles.

This function also computes the contribution to the background pressure and accelerations due to a body force or gravity.

The body forces are damped according to Eq. (13) in [Adami2012] to avoid instantaneous accelerations. By default, damping is neglected.

initialize(d_idx, d_au, d_av, d_aw, d_auhat, d_avhat, d_awhat)[source]¶

loop(d_idx, s_idx, d_m, d_rho, s_rho, d_au, d_av, d_aw, d_p, s_p, d_auhat, d_avhat, d_awhat, d_V, s_V, DWIJ)[source]¶

post_loop(d_idx, d_au, d_av, d_aw, t)[source]¶

class pysph.sph.wc.transport_velocity.MomentumEquationViscosity(dest, sources, nu)[source]¶

Bases: pysph.sph.equation.Equation

Momentum equation for the Transport Velocity Formulation: Viscosity

Eq. (8) in [Adami2013]:

\[\frac{d \boldsymbol{v}_a}{dt} = \frac{1}{m_a}\sum_b (V_a^2 + V_b^2)\left[ \bar{\eta}_{ab}\hat{r}_{ab}\cdot \nabla_a W_{ab} \frac{\boldsymbol{v}_{ab}}{|\boldsymbol{r}_{ab}|}\right]\]

where

\[\bar{\eta}_{ab} = \frac{2\eta_a \eta_b}{\eta_a + \eta_b}\]

Parameters:	nu (float) – kinematic viscosity

initialize(d_idx, d_au, d_av, d_aw)[source]¶

loop(d_idx, s_idx, d_rho, s_rho, d_m, d_V, s_V, d_au, d_av, d_aw, R2IJ, EPS, DWIJ, VIJ, XIJ)[source]¶

class pysph.sph.wc.transport_velocity.SetWallVelocity(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Extrapolating the fluid velocity on to the wall

Eq. (22) in [Adami2012]:

\[\tilde{\boldsymbol{v}}_a = \frac{\sum_b\boldsymbol{v}_b W_{ab}} {\sum_b W_{ab}}\]

Notes

The destination particle array for this equation should define the filtered velocity variables \(uf, vf, wf\).

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_uf, d_vf, d_wf, d_wij)[source]¶

loop(d_idx, s_idx, d_uf, d_vf, d_wf, s_u, s_v, s_w, d_wij, WIJ)[source]¶

post_loop(d_uf, d_vf, d_wf, d_wij, d_idx, d_ug, d_vg, d_wg, d_u, d_v, d_w)[source]¶

class pysph.sph.wc.transport_velocity.SolidWallNoSlipBC(dest, sources, nu)[source]¶

Bases: pysph.sph.equation.Equation

Solid wall boundary condition [Adami2012]

This boundary condition is to be used with fixed ghost particles in SPH simulations and is formulated for the general case of moving boundaries.

The velocity and pressure of the fluid particles is extrapolated to the ghost particles and these values are used in the equations of motion.

No-penetration:

Ghost particles participate in the continuity and state equations with fluid particles. This means as fluid particles approach the wall, the pressure of the ghost particles increases to generate a repulsion force that prevents particle penetration.

No-slip:

Extrapolation is used to set the dummy velocity of the ghost particles for viscous interaction. First, the smoothed velocity field of the fluid phase is extrapolated to the wall particles:

\[\tilde{v}_a = \frac{\sum_b v_b W_{ab}}{\sum_b W_{ab}}\]

In the second step, for the viscous interaction in Eqs. (10) in [Adami2012] and Eq. (8) in [Adami2013], the velocity of the ghost particles is assigned as:

\[v_b = 2v_w -\tilde{v}_a,\]

where \(v_w\) is the prescribed wall velocity and \(v_b\) is the ghost particle in the interaction.

Parameters:	nu (float) – kinematic viscosity

Notes

For this equation the destination particle array should be the fluid and the source should be ghost or boundary particles. The boundary particles must define a prescribed velocity \(u_0, v_0, w_0\)

initialize(d_idx, d_au, d_av, d_aw)[source]¶

loop(d_idx, s_idx, d_m, d_rho, s_rho, d_V, s_V, d_u, d_v, d_w, d_au, d_av, d_aw, s_ug, s_vg, s_wg, DWIJ, R2IJ, EPS, XIJ)[source]¶

class pysph.sph.wc.transport_velocity.SolidWallPressureBC(dest, sources, rho0, p0, b=1.0, gx=0.0, gy=0.0, gz=0.0)[source]¶

Bases: pysph.sph.equation.Equation

Solid wall pressure boundary condition [Adami2012]

This boundary condition is to be used with fixed ghost particles in SPH simulations and is formulated for the general case of moving boundaries.

The velocity and pressure of the fluid particles is extrapolated to the ghost particles and these values are used in the equations of motion.

Pressure boundary condition:

The pressure of the ghost particle is also calculated from the fluid particle by interpolation using:

\[p_g = \frac{\sum_f p_f W_{gf} + \boldsymbol{g - a_g} \cdot \sum_f \rho_f \boldsymbol{r}_{gf}W_{gf}}{\sum_f W_{gf}},\]

where the subscripts g and f relate to the ghost and fluid particles respectively.

Density of the wall particle is then set using this pressure

\[\rho_w=\rho_0\left(\frac{p_w - \mathcal{X}}{p_0} + 1\right)^{\frac{1}{\gamma}}\]

Parameters:	rho0 (float) – reference density p0 (float) – reference pressure b (float) – constant (default 1.0) gx (float) – Body force per unit mass along the x-axis gy (float) – Body force per unit mass along the y-axis gz (float) – Body force per unit mass along the z-axis

Notes

For a two fluid system (boundary, fluid), this equation must be instantiated with boundary as the destination and fluid as the source.

The boundary particle array must additionally define a property \(wij\) for the denominator in Eq. (27) from [Adami2012]. This array sums the kernel terms from the ghost particle to the fluid particle.

initialize(d_idx, d_p, d_wij)[source]¶

loop(d_idx, s_idx, d_p, s_p, d_wij, s_rho, d_au, d_av, d_aw, WIJ, XIJ)[source]¶

post_loop(d_idx, d_wij, d_p, d_rho)[source]¶

class pysph.sph.wc.transport_velocity.StateEquation(dest, sources, p0, rho0, b=1.0)[source]¶

Bases: pysph.sph.equation.Equation

Generalized Weakly Compressible Equation of State

\[p_a = p_0\left[ \left(\frac{\rho}{\rho_0}\right)^\gamma - b \right] + \mathcal{X}\]

Notes

This is the generalized Tait’s equation of state and the suggested values in [Adami2013] are \(\mathcal{X} = 0\), \(\gamma=1\) and \(b = 1\).

The reference pressure \(p_0\) is calculated from the artificial sound speed and reference density:

\[p_0 = \frac{c^2\rho_0}{\gamma}\]

Parameters:	p0 (float) – reference pressure rho0 (float) – reference density b (float) – constant (default 1.0).

loop(d_idx, d_p, d_rho)[source]¶

class pysph.sph.wc.transport_velocity.SummationDensity(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Summation density with volume summation

In addition to the standard summation density, the number density for the particle is also computed. The number density is important for multi-phase flows to define a local particle volume independent of the material density.

\[\begin{split}\rho_a = \sum_b m_b W_{ab}\\\end{split}\]\[\mathcal{V}_a = \frac{1}{\sum_b W_{ab}}\]

Notes

For this equation, the destination particle array must define the variable V for particle volume.

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_V, d_rho)[source]¶

loop(d_idx, d_V, d_rho, d_m, WIJ)[source]¶

class pysph.sph.wc.transport_velocity.VolumeFromMassDensity(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Set the inverse volume using mass density

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

loop(d_idx, d_V, d_rho, d_m)[source]¶

class pysph.sph.wc.transport_velocity.VolumeSummation(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Number density for volume computation

See SummationDensity

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_V)[source]¶

loop(d_idx, d_V, WIJ)[source]¶

SPH Boundary Equations¶

class pysph.sph.boundary_equations.MonaghanBoundaryForce(dest, sources, deltap)[source]¶

Bases: pysph.sph.equation.Equation

loop(d_idx, s_idx, s_m, s_rho, d_m, d_cs, s_cs, d_h, s_tx, s_ty, s_tz, s_nx, s_ny, s_nz, d_au, d_av, d_aw, XIJ)[source]¶

class pysph.sph.boundary_equations.MonaghanKajtarBoundaryForce(dest, sources, K=None, beta=None, h=None)[source]¶

Bases: pysph.sph.equation.Equation

loop(d_idx, s_idx, d_m, s_m, d_au, d_av, d_aw, RIJ, R2IJ, XIJ)[source]¶

pysph.sph.boundary_equations.wendland_quintic(rij=1.0, h=1.0)[source]¶

Basic Equations for Solid Mechanics¶

References

class pysph.sph.solid_mech.basic.EnergyEquationWithStress2D(dest, sources, alpha=1.0, beta=1.0, eta=0.01)[source]¶

Bases: pysph.sph.equation.Equation

initialize(d_idx, d_ae)[source]¶

loop(d_idx, s_idx, s_m, d_rho, s_rho, d_p, s_p, d_cs, s_cs, d_ae, XIJ, VIJ, DWIJ, HIJ, R2IJ, RHOIJ1)[source]¶

post_loop(d_idx, d_rho, d_s00, d_s01, d_s11, s_s00, s_s01, s_s11, d_v00, d_v01, d_v10, d_v11, d_ae)[source]¶

class pysph.sph.solid_mech.basic.HookesDeviatoricStressRate2D(dest, sources, shear_mod)[source]¶

Bases: pysph.sph.equation.Equation

Rate of change of stress (2D)

\[\frac{dS^{ij}}{dt} = 2\mu\left(\epsilon^{ij} - \frac{1}{3}\delta^{ij} \epsilon^{ij}\right) + S^{ik}\Omega^{jk} + \Omega^{ik}S^{kj}\]

where

\[\begin{split}\epsilon^{ij} = \frac{1}{2}\left(\frac{\partial v^i}{\partial x^j} + \frac{\partial v^j}{\partial x^i}\right)\\\end{split}\]\[\Omega^{ij} = \frac{1}{2}\left(\frac{\partial v^i}{\partial x^j} - \frac{\partial v^j}{\partial x^i} \right)\]

Parameters:	shear_mod (float) – shear modulus (\(\mu\))

initialize(d_idx, d_as00, d_as01, d_as11)[source]¶

loop(d_idx, d_s00, d_s01, d_s11, d_v00, d_v01, d_v10, d_v11, d_as00, d_as01, d_as11)[source]¶

class pysph.sph.solid_mech.basic.MomentumEquationWithStress2D(dest, sources, wdeltap=-1, n=1)[source]¶

Bases: pysph.sph.equation.Equation

Momentum Equation with Artificial Stress

\[\frac{D\vec{v_a}^i}{Dt} = \sum_b m_b\left(\frac{\sigma_a^{ij}}{\rho_a^2} +\frac{\sigma_b^{ij}}{\rho_b^2} + R_{ab}^{ij}f^n \right)\nabla_a W_{ab}\]

where

\[\begin{split}f_{ab} = \frac{W(r_{ab})}{W(\Delta p)}\\\end{split}\]\[R_{ab}^{ij} = R_{a}^{ij} + R_{b}^{ij}\]

Parameters:	wdeltap (float) – evaluated value of \(W(\Delta p)\) n (float) – constant with_correction (bool) – switch for using tensile instability correction

initialize(d_idx, d_au, d_av)[source]¶

loop(d_idx, s_idx, d_rho, s_rho, s_m, d_p, s_p, d_s00, d_s01, d_s11, s_s00, s_s01, s_s11, d_r00, d_r01, d_r11, s_r00, s_r01, s_r11, d_au, d_av, WIJ, DWIJ)[source]¶

class pysph.sph.solid_mech.basic.MonaghanArtificialStress(dest, sources, eps=0.3)[source]¶

Bases: pysph.sph.equation.Equation

Artificial stress to remove tensile instability

The dispersion relations in [Gray2001] are used to determine the different components of \(R\).

Angle of rotation for particle \(a\)

\[\tan{2 \theta_a} = \frac{2\sigma_a^{xy}}{\sigma_a^{xx} - \sigma_a^{yy}}\]

In rotated frame, the new components of the stress tensor are

\[\begin{split}\bar{\sigma}_a^{xx} = \cos^2{\theta_a} \sigma_a^{xx} + 2\sin{\theta_a} \cos{\theta_a}\sigma_a^{xy} + \sin^2{\theta_a}\sigma_a^{yy}\\\end{split}\]\[\bar{\sigma}_a^{yy} = \sin^2{\theta_a} \sigma_a^{xx} + 2\sin{\theta_a} \cos{\theta_a}\sigma_a^{xy} + \cos^2{\theta_a}\sigma_a^{yy}\]

Components of \(R\) in rotated frame:

\[\begin{split}\bar{R}_{a}^{xx}=\begin{cases}-\epsilon\frac{\bar{\sigma}_{a}^{xx}} {\rho^{2}} & \bar{\sigma}_{a}^{xx}>0\\0 & \bar{\sigma}_{a}^{xx}\leq0 \end{cases}\\\end{split}\]\[\begin{split}\bar{R}_{a}^{yy}=\begin{cases}-\epsilon\frac{\bar{\sigma}_{a}^{yy}} {\rho^{2}} & \bar{\sigma}_{a}^{yy}>0\\0 & \bar{\sigma}_{a}^{yy}\leq0 \end{cases}\end{split}\]

Components of \(R\) in original frame:

\[\begin{split}R_a^{xx} = \cos^2{\theta_a} \bar{R}_a^{xx} + \sin^2{\theta_a} \bar{R}_a^{yy}\\\end{split}\]\[\begin{split}R_a^{yy} = \sin^2{\theta_a} \bar{R}_a^{xx} + \cos^2{\theta_a} \bar{R}_a^{yy}\\\end{split}\]\[R_a^{xy} = \sin{\theta_a} \cos{\theta_a}\left(\bar{R}_a^{xx} - \bar{R}_a^{yy}\right)\]

Parameters:	eps (float) – constant

loop(d_idx, d_rho, d_p, d_s00, d_s01, d_s02, d_s11, d_s12, d_s22, d_r00, d_r01, d_r02, d_r11, d_r12, d_r22)[source]¶

Compute the stress terms

Parameters:	d_sxx (DoubleArray) – Stress Tensor Deviatoric components (Symmetric) d_rxx (DoubleArray) – Artificial stress components (Symmetric)

Equations for the High Velocity Impact Problems¶

class pysph.sph.solid_mech.hvi.MieGruneisenEOS(dest, sources, gamma, r0, c0, S)[source]¶

Bases: pysph.sph.equation.Equation

loop(d_idx, d_p, d_rho, d_e)[source]¶

class pysph.sph.solid_mech.hvi.VonMisesPlasticity2D(dest, sources, flow_stress)[source]¶

Bases: pysph.sph.equation.Equation

loop(d_idx, d_s00, d_s01, d_s11)[source]¶

Gas Dynamics¶

Basic equations for Gas-dynamics

class pysph.sph.gas_dynamics.basic.ADKEAccelerations(dest, sources, alpha, beta, g1, g2, k, eps)[source]¶

Bases: pysph.sph.equation.Equation

initialize(d_idx, d_au, d_av, d_aw, d_ae)[source]¶

loop(d_idx, s_idx, d_au, d_av, d_aw, d_ae, d_p, s_p, d_rho, s_rho, d_m, s_m, d_cs, s_cs, s_e, d_e, s_h, d_h, s_div, d_div, DWIJ, HIJ, XIJ, VIJ, R2IJ, EPS, RHOIJ, RHOIJ1)[source]¶

class pysph.sph.gas_dynamics.basic.IdealGasEOS(dest, sources, gamma)[source]¶

Bases: pysph.sph.equation.Equation

loop(d_idx, d_p, d_rho, d_e, d_cs)[source]¶

class pysph.sph.gas_dynamics.basic.MPMAccelerations(dest, sources, beta=2.0, update_alpha1=False, update_alpha2=False, alpha1_min=0.1, alpha2_min=0.1, sigma=0.1)[source]¶

Bases: pysph.sph.equation.Equation

initialize(d_idx, d_au, d_av, d_aw, d_ae, d_am, d_aalpha1, d_aalpha2, d_del2e)[source]¶

loop(d_idx, s_idx, d_m, s_m, d_p, s_p, d_cs, s_cs, d_e, s_e, d_rho, s_rho, d_au, d_av, d_aw, d_ae, d_omega, s_omega, XIJ, VIJ, DWI, DWJ, DWIJ, HIJ, d_del2e, d_alpha1, s_alpha1, d_alpha2, s_alpha2, EPS, RIJ, R2IJ, RHOIJ, DT_ADAPT)[source]¶

post_loop(d_idx, d_h, d_cs, d_alpha1, d_aalpha1, d_div, d_del2e, d_e, d_alpha2, d_aalpha2)[source]¶

class pysph.sph.gas_dynamics.basic.Monaghan92Accelerations(dest, sources, alpha=1.0, beta=2.0)[source]¶

Bases: pysph.sph.equation.Equation

initialize(d_idx, d_au, d_av, d_aw, d_ae)[source]¶

loop(d_idx, s_idx, d_rho, s_rho, d_p, s_p, d_cs, s_cs, d_au, d_av, d_aw, d_ae, s_m, VIJ, DWIJ, XIJ, EPS, HIJ, R2IJ, RHOIJ1)[source]¶

class pysph.sph.gas_dynamics.basic.ScaleSmoothingLength(dest, sources, factor=2.0)[source]¶

Bases: pysph.sph.equation.Equation

loop(d_idx, d_h)[source]¶

class pysph.sph.gas_dynamics.basic.SummationDensity(dest, sources, dim, density_iterations=False, iterate_only_once=False, k=1.2, htol=1e-06)[source]¶

Bases: pysph.sph.equation.Equation

Summation density with iterative solution of the smoothing lengths.

Parameters:

density_iterations : bint

Flag to indicate density iterations are required.

iterate_only_once : bint

Flag to indicate if only one iteration is required

k : double

Kernel scaling factor

htol : double

Iteration tolerance

converged()[source]¶

initialize(d_idx, d_rho, d_div, d_grhox, d_grhoy, d_arho, d_dwdh)[source]¶

loop(d_idx, s_idx, d_rho, d_grhox, d_grhoy, d_arho, d_dwdh, s_m, d_converged, VIJ, WI, DWI, GHI)[source]¶

post_loop(d_idx, d_arho, d_rho, d_div, d_omega, d_dwdh, d_h0, d_h, d_m, d_ah, d_converged)[source]¶

class pysph.sph.gas_dynamics.basic.SummationDensityADKE(dest, sources, k=1.0, eps=0.0)[source]¶

Bases: pysph.sph.equation.Equation

References

initialize(d_idx, d_arho, d_rho, d_h, d_h0)[source]¶

loop(d_idx, d_rho, d_arho, s_idx, s_m, VIJ, DWI, WIJ)[source]¶

post_loop(d_idx, d_rho, d_arho, d_div, d_logrho)[source]¶

reduce(d_rho, d_h, d_h0, dst)[source]¶

class pysph.sph.gas_dynamics.basic.UpdateSmoothingLengthFromVolume(dest, sources, dim, k=1.2)[source]¶

Bases: pysph.sph.equation.Equation

loop(d_idx, d_m, d_rho, d_h)[source]¶

Surface tension¶

Implementation of the equations used for surface tension modelling, for example in KHI simulations. The references are as under:

• M. Shadloo, M. Yildiz, “Numerical modelling of Kelvin-Helmholtz isntability using smoothed particle hydrodynamics”, IJNME, 2011, 87, pp 988–1006 [SY11]
• Joseph P. Morris “Simulating surface tension with smoothed particle hydrodynamics”, JCP, 2000, 33, pp 333–353 [JM00]
• Adami et al. “A new surface-tension formulation for multi-phase SPH using a reproducing divergence approximation”, JCP 2010, 229, pp 5011–5021 [A10]

class pysph.sph.surface_tension.AdamiColorGradient(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Gradient of color Eq. (14) in [A10]

\[\nabla c_a = \frac{1}{V_a}\sum_b \left[V_a^2 + V_b^2 \right]\tilde{c}_{ab}\nabla_a W_{ab}\,,\]

where, the average \(\tilde{c}_{ab}\) is defined as

\[\tilde{c}_{ab} = \frac{\rho_b}{\rho_a + \rho_b}c_a + \frac{\rho_a}{\rho_a + \rho_b}c_b\]

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_cx, d_cy, d_cz, d_nx, d_ny, d_nz, d_ddelta, d_N)[source]¶

loop(d_idx, s_idx, d_V, s_V, d_rho, s_rho, d_cx, d_cy, d_cz, d_color, s_color, DWIJ)[source]¶

post_loop(d_idx, d_cx, d_cy, d_cz, d_nx, d_ny, d_nz, d_ddelta, d_N)[source]¶

class pysph.sph.surface_tension.AdamiReproducingDivergence(dest, sources, dim)[source]¶

Bases: pysph.sph.equation.Equation

Reproducing divergence approximation Eq. (20) in [A10] to compute the curvature

\[\nabla \cdot \boldsymbol{\phi}_a = d\frac{\sum_b \boldsymbol{\phi}_{ab}\cdot \nabla_a W_{ab}V_b}{\sum_b\boldsymbol{x}_{ab}\cdot \nabla_a W_{ab} V_b}\]

initialize(d_idx, d_kappa, d_wij_sum)[source]¶

loop(d_idx, s_idx, d_kappa, d_wij_sum, d_nx, d_ny, d_nz, s_nx, s_ny, s_nz, d_V, s_V, DWIJ, XIJ, RIJ, EPS)[source]¶

post_loop(d_idx, d_kappa, d_wij_sum)[source]¶

class pysph.sph.surface_tension.CSFSurfaceTensionForce(dest, sources, sigma=0.1)[source]¶

Bases: pysph.sph.equation.Equation

Acceleration due to surface tension force Eq. (25) in [JM00]:

Note that as per Eq. (17) in [JM00], the un-normalized normal is basically the gradient of the color function. The acceleration term therefore depends on the gradient of the color field.

initialize(d_idx, d_au, d_av, d_aw)[source]¶

loop(d_idx, d_au, d_av, d_aw, d_kappa, d_cx, d_cy, d_cz, d_rho)[source]¶

class pysph.sph.surface_tension.ColorGradientUsingNumberDensity(dest, sources, epsilon=1e-06)[source]¶

Bases: pysph.sph.equation.Equation

Gradient of the color function using Eq. (13) of [SY11]:

\[\nabla C_a = \sum_b \frac{2 C_b - C_a}{\psi_a + \psi_a} \nabla_{a} W_{ab}\]

Using the gradient of the color function, the normal and discretized dirac delta is calculated in the post loop.

Singularities are avoided as per the recommendation by [JM00] (see eqs 20 & 21) using the parameter \(\epsilon\)

initialize(d_idx, d_cx, d_cy, d_cz, d_nx, d_ny, d_nz, d_ddelta, d_N)[source]¶

loop(d_idx, s_idx, d_scolor, s_scolor, d_cx, d_cy, d_cz, d_V, s_V, DWIJ)[source]¶

post_loop(d_idx, d_cx, d_cy, d_cz, d_nx, d_ny, d_nz, d_N, d_ddelta)[source]¶

class pysph.sph.surface_tension.InterfaceCurvatureFromNumberDensity(dest, sources, with_morris_correction=True)[source]¶

Bases: pysph.sph.equation.Equation

Interface curvature using number density. Eq. (15) in [SY11]:

\[\kappa_a = \sum_b \frac{2.0}{\psi_a + \psi_b} \left(\boldsymbol{n_a} - \boldsymbol{n_b}\right) \cdot \nabla_a W_{ab}\]

initialize(d_idx, d_kappa, d_wij_sum)[source]¶

loop(d_idx, s_idx, d_kappa, d_nx, d_ny, d_nz, s_nx, s_ny, s_nz, d_V, s_V, d_N, s_N, d_wij_sum, s_rho, s_m, WIJ, DWIJ)[source]¶

post_loop(d_idx, d_wij_sum, d_nx, d_kappa)[source]¶

class pysph.sph.surface_tension.MorrisColorGradient(dest, sources, epsilon=1e-06)[source]¶

Bases: pysph.sph.equation.Equation

Gradient of the color function using Eq. (17) of [JM00]:

\[\nabla c_a = \sum_b \frac{m_b}{\rho_b}(c_b - c_a) \nabla_{a} W_{ab}\,,\]

where a smoothed representation of the color is used in the equation. Using the gradient of the color function, the normal and discretized dirac delta is calculated in the post loop.

Singularities are avoided as per the recommendation by [JM00] (see eqs 20 & 21) using the parameter \(\epsilon\)

initialize(d_idx, d_cx, d_cy, d_cz, d_nx, d_ny, d_nz, d_ddelta, d_N)[source]¶

loop(d_idx, s_idx, d_scolor, s_scolor, d_cx, d_cy, d_cz, s_m, s_rho, DWIJ)[source]¶

post_loop(d_idx, d_cx, d_cy, d_cz, d_nx, d_ny, d_nz, d_N, d_ddelta)[source]¶

class pysph.sph.surface_tension.SY11ColorGradient(dest, sources, epsilon=1e-06)[source]¶

Bases: pysph.sph.equation.Equation

Gradient of the color function using Eq. (13) of [SY11]:

\[\nabla C_a = \sum_b \frac{2 C_b - C_a}{\psi_a + \psi_a} \nabla_{a} W_{ab}\]

Using the gradient of the color function, the normal and discretized dirac delta is calculated in the post loop.

initialize(d_idx, d_cx, d_cy, d_cz, d_nx, d_ny, d_nz, d_ddelta, d_N)[source]¶

loop(d_idx, s_idx, d_color, s_color, d_cx, d_cy, d_cz, d_V, s_V, DWIJ)[source]¶

post_loop(d_idx, d_cx, d_cy, d_cz, d_nx, d_ny, d_nz, d_N, d_ddelta)[source]¶

class pysph.sph.surface_tension.SY11DiracDelta(dest, sources, epsilon=1e-06)[source]¶

Bases: pysph.sph.equation.Equation

Discretized dirac-delta for the SY11 formulation Eq. (14) in [SY11]

This is essentially the same as computing the color gradient, the only difference being that this might be called with a reduced smoothing length.

Note that the normals should be computed using the SY11ColorGradient equation. This function will effectively overwrite the color gradient.

initialize(d_idx, d_cx, d_cy, d_cz, d_ddelta)[source]¶

loop(d_idx, s_idx, d_color, s_color, d_cx, d_cy, d_cz, d_V, s_V, DWIJ)[source]¶

post_loop(d_idx, d_cx, d_cy, d_cz, d_nx, d_ny, d_nz, d_N, d_ddelta)[source]¶

class pysph.sph.surface_tension.ShadlooYildizSurfaceTensionForce(dest, sources, sigma=0.1)[source]¶

Bases: pysph.sph.equation.Equation

Acceleration due to surface tension force Eq. (7,9) in [SY11]:

where, \(\delta^s\) is the discretized dirac delta function, \(\boldsymbol{n}\) is the interface normal, \(\kappa\) is the discretized interface curvature and \(\sigma\) is the surface tension force constant.

initialize(d_idx, d_au, d_av, d_aw)[source]¶

loop(d_idx, d_au, d_av, d_aw, d_kappa, d_nx, d_ny, d_nz, d_m, d_rho, d_ddelta)[source]¶

class pysph.sph.surface_tension.SmoothedColor(dest, sources, smooth=False)[source]¶

Bases: pysph.sph.equation.Equation

Smoothed color function. Eq. (17) in [JM00]

\[c_a = \sum_b \frac{m_b}{\rho_b} c_b^i \nabla_a W_{ab}\,,\]

where, \(c_b^i\) is the color index associated with a particle.

initialize(d_idx, d_scolor)[source]¶

loop(d_idx, s_idx, d_rho, s_rho, s_m, s_color, d_scolor, WIJ)[source]¶

post_loop(d_idx, d_color, d_scolor)[source]¶

Implicit Incompressible SPH¶

The basic equations for the IISPH formulation of

M. Ihmsen, J. Cornelis, B. Solenthaler, C. Horvath, M. Teschner, “Implicit Incompressible SPH,” IEEE Transactions on Visualization and Computer Graphics, vol. 20, no. 3, pp. 426-435, March 2014. http://dx.doi.org/10.1109/TVCG.2013.105

class pysph.sph.iisph.AdvectionAcceleration(dest, sources, gx=0.0, gy=0.0, gz=0.0)[source]¶

Bases: pysph.sph.equation.Equation

initialize(d_idx, d_au, d_av, d_aw, d_uadv, d_vadv, d_wadv)[source]¶

loop()[source]¶

post_loop(d_idx, d_au, d_av, d_aw, d_uadv, d_vadv, d_wadv, d_u, d_v, d_w, dt=0.0)[source]¶

class pysph.sph.iisph.ComputeAII(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_aii)[source]¶

loop(d_idx, d_aii, d_dii0, d_dii1, d_dii2, d_m, d_rho, s_idx, s_m, s_rho, DWIJ, dt=0.0)[source]¶

class pysph.sph.iisph.ComputeAIIBoundary(dest, sources, rho0)[source]¶

Bases: pysph.sph.equation.Equation

This is important and not really discussed in the original IISPH paper.

loop(d_idx, d_aii, d_dii0, d_dii1, d_dii2, d_rho, s_idx, s_m, s_V, DWIJ, dt=0.0)[source]¶

class pysph.sph.iisph.ComputeDII(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_dii0, d_dii1, d_dii2)[source]¶

loop(d_idx, d_rho, d_dii0, d_dii1, d_dii2, s_idx, s_m, DWIJ, dt=0.0)[source]¶

class pysph.sph.iisph.ComputeDIIBoundary(dest, sources, rho0)[source]¶

Bases: pysph.sph.equation.Equation

loop(d_idx, d_dii0, d_dii1, d_dii2, d_rho, s_idx, s_m, s_V, DWIJ, dt=0.0)[source]¶

class pysph.sph.iisph.ComputeDIJPJ(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_dijpj0, d_dijpj1, d_dijpj2)[source]¶

loop(d_idx, d_dijpj0, d_dijpj1, d_dijpj2, s_idx, s_m, s_rho, s_piter, DWIJ, dt=0.0)[source]¶

class pysph.sph.iisph.ComputeRhoAdvection(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_rho_adv, d_rho, d_p0, d_p, d_piter, d_aii)[source]¶

loop(d_idx, d_rho, d_rho_adv, d_uadv, d_vadv, d_wadv, d_u, d_v, d_w, s_idx, s_m, s_uadv, s_vadv, s_wadv, DWIJ, dt=0.0)[source]¶

class pysph.sph.iisph.ComputeRhoBoundary(dest, sources, rho0)[source]¶

Bases: pysph.sph.equation.Equation

loop(d_idx, d_rho, d_rho_adv, d_uadv, d_vadv, d_wadv, s_idx, s_u, s_v, s_w, s_V, WIJ, DWIJ, dt=0.0)[source]¶

class pysph.sph.iisph.IISPHStep[source]¶

Bases: pysph.sph.integrator_step.IntegratorStep

A straightforward and simple integrator to be used for IISPH.

initialize()[source]¶

stage1(d_idx, d_x, d_y, d_z, d_u, d_v, d_w, d_uadv, d_vadv, d_wadv, d_au, d_av, d_aw, d_ax, d_ay, d_az, dt)[source]¶

class pysph.sph.iisph.NormalizedSummationDensity(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_rho, d_rho_adv, d_rho0, d_V)[source]¶

loop(d_idx, d_rho, d_rho_adv, d_V, s_idx, s_m, s_rho0, WIJ)[source]¶

post_loop(d_idx, d_rho, d_rho_adv)[source]¶

class pysph.sph.iisph.NumberDensity(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_V)[source]¶

loop(d_idx, d_V, WIJ)[source]¶

class pysph.sph.iisph.PressureForce(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_au, d_av, d_aw)[source]¶

loop(d_idx, d_rho, d_p, d_au, d_av, d_aw, s_idx, s_m, s_rho, s_p, DWIJ)[source]¶

post_loop(d_idx, d_au, d_av, d_aw, d_uadv, d_vadv, d_wadv, DT_ADAPT)[source]¶

class pysph.sph.iisph.PressureForceBoundary(dest, sources, rho0)[source]¶

Bases: pysph.sph.equation.Equation

loop(d_idx, d_rho, d_au, d_av, d_aw, d_p, s_idx, s_V, DWIJ)[source]¶

class pysph.sph.iisph.PressureSolve(dest, sources, rho0, omega=0.5, tolerance=0.01, debug=False)[source]¶

Bases: pysph.sph.equation.Equation

converged()[source]¶

initialize(d_idx, d_p, d_compression)[source]¶

loop(d_idx, d_p, d_piter, d_rho, d_m, d_dijpj0, d_dijpj1, d_dijpj2, s_idx, s_m, s_dii0, s_dii1, s_dii2, s_piter, s_dijpj0, s_dijpj1, s_dijpj2, DWIJ, dt=0.0)[source]¶

post_loop(d_idx, d_piter, d_p0, d_p, d_aii, d_rho_adv, d_rho, d_compression, dt=0.0)[source]¶

reduce(dst)[source]¶

class pysph.sph.iisph.PressureSolveBoundary(dest, sources, rho0)[source]¶

Bases: pysph.sph.equation.Equation

loop(d_idx, d_p, d_rho, d_dijpj0, d_dijpj1, d_dijpj2, s_idx, s_V, DWIJ)[source]¶

class pysph.sph.iisph.SummationDensity(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_rho)[source]¶

loop(d_idx, d_rho, s_idx, s_m, WIJ)[source]¶

class pysph.sph.iisph.SummationDensityBoundary(dest, sources, rho0)[source]¶

Bases: pysph.sph.equation.Equation

loop(d_idx, d_rho, s_idx, s_V, WIJ)[source]¶

class pysph.sph.iisph.ViscosityAcceleration(dest, sources, nu)[source]¶

Bases: pysph.sph.equation.Equation

loop(d_idx, d_au, d_av, d_aw, s_idx, s_m, EPS, VIJ, XIJ, RHOIJ1, R2IJ, DWIJ)[source]¶

class pysph.sph.iisph.ViscosityAccelerationBoundary(dest, sources, rho0, nu)[source]¶

Bases: pysph.sph.equation.Equation

The acceleration on the fluid due to a boundary.

loop(d_idx, d_au, d_av, d_aw, d_rho, s_idx, s_V, EPS, VIJ, XIJ, R2IJ, DWIJ)[source]¶

Rigid body motion¶

Rigid body related equations.

class pysph.sph.rigid_body.BodyForce(dest, sources, gx=0.0, gy=0.0, gz=0.0)[source]¶

Bases: pysph.sph.equation.Equation

initialize(d_idx, d_m, d_fx, d_fy, d_fz)[source]¶

class pysph.sph.rigid_body.EulerStepRigidBody[source]¶

Bases: pysph.sph.integrator_step.IntegratorStep

Fast but inaccurate integrator. Use this for testing

initialize()[source]¶

stage1(d_idx, d_u, d_v, d_w, d_x, d_y, d_z, d_omega, d_omega_dot, d_vc, d_ac, d_num_body, dt=0.0)[source]¶

class pysph.sph.rigid_body.NumberDensity(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_V)[source]¶

loop(d_idx, d_V, WIJ)[source]¶

class pysph.sph.rigid_body.PressureRigidBody(dest, sources, rho0)[source]¶

Bases: pysph.sph.equation.Equation

The pressure acceleration on the fluid/solid due to a boundary. Implemented from Akinci et al. http://dx.doi.org/10.1145/2185520.2185558

Use this with the fluid as a destination and body as source.

loop(d_idx, d_m, d_rho, d_au, d_av, d_aw, d_p, s_idx, s_V, s_fx, s_fy, s_fz, DWIJ)[source]¶

class pysph.sph.rigid_body.RK2StepRigidBody[source]¶

Bases: pysph.sph.integrator_step.IntegratorStep

initialize(d_idx, d_x, d_y, d_z, d_x0, d_y0, d_z0, d_omega, d_omega0, d_vc, d_vc0, d_num_body)[source]¶

stage1(d_idx, d_u, d_v, d_w, d_x, d_y, d_z, d_x0, d_y0, d_z0, d_omega, d_omega_dot, d_vc, d_ac, d_omega0, d_vc0, d_num_body, dt=0.0)[source]¶

stage2(d_idx, d_u, d_v, d_w, d_x, d_y, d_z, d_x0, d_y0, d_z0, d_omega, d_omega_dot, d_vc, d_ac, d_omega0, d_vc0, d_num_body, dt=0.0)[source]¶

class pysph.sph.rigid_body.RigidBodyCollision(dest, sources, k=1.0, d=1.0, eta=1.0, kt=1.0)[source]¶

Bases: pysph.sph.equation.Equation

This is inspired from http://http.developer.nvidia.com/GPUGems3/gpugems3_ch29.html and

BK Mishra’s article on DEM http://dx.doi.org/10.1016/S0301-7516(03)00032-2

A review of computer simulation of tumbling mills by the discrete element method: Part I - contact mechanics

Note that d is a factor multiplied with the “h” of the particle.

loop(d_idx, d_fx, d_fy, d_fz, d_h, d_total_mass, XIJ, RIJ, R2IJ, VIJ)[source]¶

class pysph.sph.rigid_body.RigidBodyMoments(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_m, d_x, d_y, d_z, d_fx, d_fy, d_fz, d_mi, d_num_body, d_body_id)[source]¶

reduce(dst)[source]¶

class pysph.sph.rigid_body.RigidBodyMotion(dest, sources)[source]¶

Bases: pysph.sph.equation.Equation

Parameters:	dest (str) – name of the destination particle array sources (list of str or None) – names of the source particle arrays

initialize(d_idx, d_x, d_y, d_z, d_u, d_v, d_w, d_cm, d_vc, d_ac, d_omega, d_body_id)[source]¶

class pysph.sph.rigid_body.SummationDensityRigidBody(dest, sources, rho0)[source]¶

Bases: pysph.sph.equation.Equation

initialize(d_idx, d_rho)[source]¶

loop(d_idx, d_rho, s_idx, s_V, WIJ)[source]¶

class pysph.sph.rigid_body.ViscosityRigidBody(dest, sources, rho0, nu)[source]¶

Bases: pysph.sph.equation.Equation

The viscous acceleration on the fluid/solid due to a boundary. Implemented from Akinci et al. http://dx.doi.org/10.1145/2185520.2185558

Use this with the fluid as a destination and body as source.

loop(d_idx, d_m, d_au, d_av, d_aw, d_rho, s_idx, s_V, s_fx, s_fy, s_fz, EPS, VIJ, XIJ, R2IJ, DWIJ)[source]¶

pysph.sph.rigid_body.declare(*args)[source]¶

pysph.sph.rigid_body.get_alpha_dot()[source]¶

Use sympy to perform most of the math and use the resulting formulae to calculate:

inv(I) ( au - w x (I w))

pysph.sph.rigid_body.get_torque()[source]¶

Use sympy to perform some simple math. R x F C_m x F w x r

pysph.sph.rigid_body.skew(vec)[source]¶

Miscellaneous¶

Functions to reduce array data in serial or parallel.

pysph.base.reduce_array.dummy_reduce_array(array, op='sum')[source]¶

Simply returns the array for the serial case.

pysph.base.reduce_array.mpi_reduce_array(array, op='sum')[source]¶

Reduce an array given an array and a suitable reduction operation.

Currently, only ‘sum’, ‘max’, ‘min’ and ‘prod’ are supported.

Parameters

• array: numpy.ndarray: Any numpy array (1D).
• op: str: reduction operation, one of (‘sum’, ‘prod’, ‘min’, ‘max’)

pysph.base.reduce_array.parallel_reduce_array(array, op='sum')¶

Reduce an array given an array and a suitable reduction operation.

Currently, only ‘sum’, ‘max’, ‘min’ and ‘prod’ are supported.

Parameters

• array: numpy.ndarray: Any numpy array (1D).
• op: str: reduction operation, one of (‘sum’, ‘prod’, ‘min’, ‘max’)

pysph.base.reduce_array.serial_reduce_array(array, op='sum')[source]¶

Reduce an array given an array and a suitable reduction operation.

Currently, only ‘sum’, ‘max’, ‘min’ and ‘prod’ are supported.

Parameters

• array: numpy.ndarray: Any numpy array (1D).
• op: str: reduction operation, one of (‘sum’, ‘prod’, ‘min’, ‘max’)

Group of equations¶

class pysph.sph.equation.Group(equations, real=True, update_nnps=False, iterate=False, max_iterations=1, min_iterations=0)[source]¶

Bases: object

A group of equations.

This class provides some support for the code generation for the collection of equations.

Constructor.

Parameters:

equations (list) – a list of equation objects. real (bool) – specifies if only non-remote/non-ghost particles should be operated on. update_nnps (bool) – specifies if the neighbors should be re-computed locally after this group iterate (bool) – specifies if the group should continue iterating until each equation’s “converged()” methods returns with a positive value. max_iterations (int) – specifies the maximum number of times this group should be iterated. min_iterations (int) – specifies the minimum number of times this group should be iterated.

Notes

When running simulations in parallel, one should typically run the summation density over all particles (both local and remote) in each processor. This is because we must update the pressure/density of the remote neighbors in the current processor. Otherwise the results can be incorrect with the remote particles having an older density. This is also the case for the TaitEOS. In these cases the group that computes the equation should set real to False.