SPH equations

class pysph.sph.equation.Equation(dest, sources)

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

converged()

Return > 0 to indicate converged iterations and < 0 otherwise.

Basic SPH Equations

References

[Monaghan1992]

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

[Monaghan2005]

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

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

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)

loop(d_idx, d_au, d_av, d_aw)

class pysph.sph.basic_equations.ContinuityEquation(dest, sources)

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)

loop(d_idx, d_arho, s_idx, s_m, DWIJ, VIJ)

class pysph.sph.basic_equations.IsothermalEOS(dest, sources, rho0, c0, p0)

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)

class pysph.sph.basic_equations.MonaghanArtificialViscosity(dest, sources, alpha=1.0, beta=1.0)

Classical Monaghan style artificial viscosity [Monaghan2005]

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

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{align}\begin{aligned}\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}\\\bar{\rho}_{ab}&=&\frac{\rho_{a}+\rho_{b}}{2}\end{aligned}\end{align} \]

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)

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)

class pysph.sph.basic_equations.SummationDensity(dest, sources)

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)

loop(d_idx, d_rho, s_idx, s_m, WIJ)

class pysph.sph.basic_equations.VelocityGradient2D(dest, sources)

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_10 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)

loop(d_idx, s_idx, s_m, s_rho, d_v00, d_v01, d_v10, d_v11, DWIJ, VIJ)

class pysph.sph.basic_equations.VelocityGradient3D(dest, sources)

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_v02, d_v10, d_v11, d_v12, d_v20, d_v21, d_v22)

loop(d_idx, s_idx, s_m, s_rho, d_v00, d_v01, d_v02, d_v10, d_v11, d_v12, d_v20, d_v21, d_v22, DWIJ, VIJ)

class pysph.sph.basic_equations.XSPHCorrection(dest, sources, eps=0.5)

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}\]

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)

loop(s_idx, d_idx, s_m, d_ax, d_ay, d_az, WIJ, RHOIJ1, VIJ)

post_loop(d_idx, d_ax, d_ay, d_az, d_u, d_v, d_w)

class pysph.sph.basic_equations.XSPHCorrectionForLeapFrog(dest, sources, eps=0.5)

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}\]

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)

loop(s_idx, d_idx, s_m, d_ax, d_ay, d_az, WIJ, RHOIJ1, VIJ)

Basic WCSPH Equations

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

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

loop(d_idx, d_arho, s_idx, s_m, d_rho, s_rho, DWIJ, XIJ, R2IJ, HIJ, EPS, d_gradrho, s_gradrho)

class pysph.sph.wc.basic.ContinuityEquationDeltaSPHPreStep(dest, sources)

Continuity equation with dissipative terms See pysph.sph.wc.basic.ContinuityEquationDeltaSPH The matrix \(L_a\) is multiplied to \(\nabla W_{ij}\) in the pysph.sph.scheme.WCSPHScheme class by using pysph.sph.wc.kernel_correction.GradientCorrectionPreStep and pysph.sph.wc.kernel_correction.GradientCorrection.

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_gradrho)

loop(d_idx, d_arho, s_idx, d_rho, s_rho, s_m, d_gradrho, DWIJ)

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)

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{align}\begin{aligned}\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}\end{aligned}\end{align} \]

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, d_dt_cfl)

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, WIJ, WDP, d_dt_cfl)

post_loop(d_idx, d_au, d_av, d_aw, d_dt_force)

class pysph.sph.wc.basic.MomentumEquationDeltaSPH(dest, sources, rho0, c0, alpha=1.0)

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

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.

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)

class pysph.sph.wc.basic.PressureGradientUsingNumberDensity(dest, sources)

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)

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)

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

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)

class pysph.sph.wc.basic.TaitEOSHGCorrection(dest, sources, rho0, c0, gamma)

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)

class pysph.sph.wc.basic.UpdateSmoothingLengthFerrari(dest, sources, dim, hdx)

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)

Viscosity functions

class pysph.sph.wc.viscosity.ClearyArtificialViscosity(dest, sources, dim, alpha=1.0)

Artificial viscosity proposed By P. Cleary:

\[\mathcal{Pi}_{ab} = -\frac{16}{\mu_a \mu_b}{\rho_a \rho_b (\mu_a + \mu_b)}\left( \frac{\boldsymbol{v}_{ab} \cdot \boldsymbol{r}_{ab}}{\boldsymbol{r}_{ab}^2 + \epsilon} \right),\]

where the viscosity is determined from the parameter \(\alpha\) as

\[\mu_a = \frac{1}{8}\alpha h_a c_a \rho_a\]

This equation is described in the 2005 review paper by Monaghan

• J. J. Monaghan, “Smoothed Particle Hydrodynamics”, Reports on Progress in Physics, 2005, 68, pp 1703–1759 [JM05]

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)

loop(d_idx, s_idx, d_m, s_m, d_rho, s_rho, d_h, s_h, d_cs, s_cs, d_au, d_av, d_aw, XIJ, VIJ, R2IJ, EPS, DWIJ)

class pysph.sph.wc.viscosity.LaminarViscosity(dest, sources, nu, eta=0.01)

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

loop(d_idx, s_idx, s_m, d_rho, s_rho, d_au, d_av, d_aw, DWIJ, XIJ, VIJ, R2IJ, HIJ)

class pysph.sph.wc.viscosity.LaminarViscosityDeltaSPH(dest, sources, dim, rho0, nu)

See section 2 of the below reference [Sun2017]

References

[Sun2017]

P. Sun, A. Colagrossi, S. Marrone, A. Zhang “The $delta$ plus-SPH model: simple procedures for a further improvement of the SPH scheme”, Computer Methods in Applied Mechanics and Engineering 315 (2017), pp. 25-49.

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

loop(d_idx, s_idx, s_m, s_rho, d_rho, d_au, d_av, d_aw, HIJ, DWIJ, R2IJ, EPS, VIJ, XIJ)

class pysph.sph.wc.viscosity.MonaghanSignalViscosityFluids(dest, sources, alpha, h)

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

loop(d_idx, s_idx, d_rho, s_rho, s_m, d_au, d_av, d_aw, d_cs, s_cs, RIJ, HIJ, VIJ, XIJ, DWIJ)

Transport Velocity Formulation

References

[Adami2012] (1,2)

S. Adami et. al “A generalized wall boundary condition for smoothed particle hydrodynamics”, Journal of Computational Physics (2012), pp. 7057–7075.

[Adami2013]

S. Adami et. al “A transport-velocity formulation for smoothed particle hydrodynamics”, Journal of Computational Physics (2013), pp. 292–307.

class pysph.sph.wc.transport_velocity.ContinuityEquation(dest, sources)

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)

loop(d_idx, s_idx, d_arho, s_m, s_rho, d_rho, VIJ, DWIJ)

class pysph.sph.wc.transport_velocity.ContinuitySolid(dest, sources)

Continuity equation for the solid’s ghost particles.

The key difference is that we use the ghost velocity ug, and not the particle velocity u.

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

loop(d_idx, s_idx, d_rho, d_u, d_v, d_w, d_arho, s_m, s_rho, s_ug, s_vg, s_wg, DWIJ)

class pysph.sph.wc.transport_velocity.MomentumEquationArtificialStress(dest, sources)

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)

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)

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

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{align}\begin{aligned}\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}\end{aligned}\end{align} \]

Parameters:

• alpha (float) – constant

• c0 (float) – speed of sound

initialize(d_idx, d_au, d_av, d_aw)

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

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

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)

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)

post_loop(d_idx, d_au, d_av, d_aw, t)

class pysph.sph.wc.transport_velocity.MomentumEquationViscosity(dest, sources, nu)

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)

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)

class pysph.sph.wc.transport_velocity.SetWallVelocity(dest, sources)

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)

loop(d_idx, s_idx, d_uf, d_vf, d_wf, s_u, s_v, s_w, d_wij, WIJ)

post_loop(d_uf, d_vf, d_wf, d_wij, d_idx, d_ug, d_vg, d_wg, d_u, d_v, d_w)

class pysph.sph.wc.transport_velocity.SolidWallNoSlipBC(dest, sources, nu)

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)

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)

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

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)

loop(d_idx, s_idx, d_p, s_p, d_wij, s_rho, d_au, d_av, d_aw, WIJ, XIJ)

post_loop(d_idx, d_wij, d_p, d_rho)

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

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)

class pysph.sph.wc.transport_velocity.SummationDensity(dest, sources)

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{align}\begin{aligned}\begin{split}\rho_a = \sum_b m_b W_{ab}\\\end{split}\\\mathcal{V}_a = \frac{1}{\sum_b W_{ab}}\end{aligned}\end{align} \]

Notes

Note that in the pysph implementation, V is the inverse volume of a particle, i.e. the equation computes V as follows:

\[\mathcal{V}_a = \sum_b W_{ab}\]

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)

loop(d_idx, d_V, d_rho, d_m, WIJ)

class pysph.sph.wc.transport_velocity.VolumeFromMassDensity(dest, sources)

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)

class pysph.sph.wc.transport_velocity.VolumeSummation(dest, sources)

Number density for volume computation

See SummationDensity

Note that the quantity V is really \(\sigma\) of the original paper, i.e. inverse of the 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)

loop(d_idx, d_V, WIJ)

Generalized Transport Velocity Formulation

Some notes on the paper,

• In the viscosity term of equation (17) a factor of ‘2’ is missing.

• A negative sign is missing from equation (22) i.e, either put a negative sign in equation (22) or at the integrator step equation(25).

• The Solid Mechanics Equations are not tested.

References

[ZhangHuAdams2017]

Chi Zhang, Xiangyu Y. Hu, Nikolaus A. Adams “A generalized transport-velocity formulation for smoothed particle hydrodynamics”, Journal of Computational Physics 237 (2017), pp. 216–232.

class pysph.sph.wc.gtvf.ContinuityEquationGTVF(dest, sources)

Evolution of density

From [ZhangHuAdams2017], equation (12),

\[\frac{\tilde{d} \rho_i}{dt} = \rho_i \sum_j \frac{m_j}{\rho_j} \nabla W_{ij} \cdot \tilde{\boldsymbol{v}}_{ij}\]

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

initialize(d_arho, d_idx)

loop(d_idx, s_idx, s_m, d_rho, s_rho, d_uhat, d_vhat, d_what, s_uhat, s_vhat, s_what, d_arho, DWIJ)

class pysph.sph.wc.gtvf.CorrectDensity(dest, sources)

Density correction

From [ZhangHuAdams2017], equation (13),

\[\rho_i = \frac{\sum_j m_j W_{ij}} {\min(1, \sum_j \frac{m_j}{\rho_j^{*}} W_{ij})}\]

where,

\[\rho_j^{*} = \text{density before this correction is applied.}\]

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_rho0, d_rhodiv)

loop(d_idx, s_idx, d_rho, d_rhodiv, s_m, WIJ, s_rho0)

post_loop(d_idx, d_rho, d_rhodiv)

class pysph.sph.wc.gtvf.DeviatoricStressRate(dest, sources, dim, G)

Stress rate for solids

From [ZhangHuAdams2017], equation (5),

\[\frac{d \boldsymbol{\sigma}'}{dt} = 2 G (\boldsymbol{\epsilon} - \frac{1}{3} \text{Tr}(\boldsymbol{\epsilon})\textbf{I}) + \boldsymbol{\sigma}' \cdot \boldsymbol{\Omega}^{T} + \boldsymbol{\Omega} \cdot \boldsymbol{\sigma}'\]

where,

\[\boldsymbol{\Omega_{i/j}} = \frac{1}{2} \left(\nabla \otimes \boldsymbol{v}_{i/j} - (\nabla \otimes \boldsymbol{v}_{i/j})^{T}\right)\]

\[\boldsymbol{\epsilon_{i/j}} = \frac{1}{2} \left(\nabla \otimes \boldsymbol{v}_{i/j} + (\nabla \otimes \boldsymbol{v}_{i/j})^{T}\right)\]

see the class VelocityGradient for \(\nabla \otimes \boldsymbol{v}_i\)

Parameters:

• dim (int) – Dimensionality of the problem.

• G (float) – value of shear modulus

initialize(d_idx, d_sigma, d_asigma, d_gradvhat)

class pysph.sph.wc.gtvf.GTVFIntegrator(**kw)

Pass fluid names and suitable IntegratorStep instances.

For example:

>>> integrator = Integrator(fluid=WCSPHStep(), solid=WCSPHStep())

where “fluid” and “solid” are the names of the particle arrays.

one_timestep(t, dt)

User written function that actually does one timestep.

This function is used in the high-performance Cython implementation. The assumptions one may make are the following:

• t and dt are passed.

• the following methods are available:

  • self.initialize()

  • self.stage1(), self.stage2() etc. depending on the number of stages available.

  • self.compute_accelerations(index=0, update_nnps=True)

  • self.do_post_stage(stage_dt, stage_count_from_1)

  • self.update_domain()

Please see any of the concrete implementations of the Integrator class to study. By default the Integrator implements a predict-evaluate-correct method, the same as PECIntegrator.

class pysph.sph.wc.gtvf.GTVFScheme(fluids, solids, dim, rho0, c0, nu, h0, pref, gx=0.0, gy=0.0, gz=0.0, b=1.0, alpha=0.0)

Parameters:

• fluids (list) – List of names of fluid particle arrays.

• solids (list) – List of names of solid particle arrays.

• dim (int) – Dimensionality of the problem.

• rho0 (float) – Reference density.

• c0 (float) – Reference speed of sound.

• nu (float) – Real viscosity of the fluid.

• h0 (float) – Reference smoothing length.

• pref (float) – reference pressure for rate of change of transport velocity.

• gx (float) – Body force acceleration components in x direction.

• gy (float) – Body force acceleration components in y direction.

• gz (float) – Body force acceleration components in z direction.

• b (float) – constant for the equation of state.

configure_solver(kernel=None, integrator_cls=None, extra_steppers=None, **kw)

Configure the solver to be generated.

Parameters:

• kernel (Kernel instance.) – Kernel to use, if none is passed a default one is used.

• integrator_cls (pysph.sph.integrator.Integrator) – Integrator class to use, use sensible default if none is passed.

• extra_steppers (dict) – Additional integration stepper instances as a dict.

• **kw (extra arguments) – Any additional keyword args are passed to the solver instance.

get_equations()

setup_properties(particles, clean=True)

Setup the particle arrays so they have the right set of properties for this scheme.

Parameters:

• particles (list) – List of particle arrays.

• clean (bool) – If True, removes any unnecessary properties.

class pysph.sph.wc.gtvf.GTVFStep

stage1(d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, d_uhat, d_vhat, d_what, d_auhat, d_avhat, d_awhat, dt)

stage2(d_idx, d_uhat, d_vhat, d_what, d_x, d_y, d_z, d_rho, d_arho, d_sigma, d_asigma, dt)

stage3(d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, dt)

class pysph.sph.wc.gtvf.MomentumEquationArtificialStress(dest, sources, dim)

Momentum equation Artificial stress for solids

See the class MomentumEquationPressureGradient for details.

Parameters:

dim (int) – Dimensionality of the problem.

loop(d_idx, s_idx, d_rho, s_rho, d_u, d_v, d_w, d_uhat, d_vhat, d_what, s_u, s_v, s_w, s_uhat, s_vhat, s_what, d_au, d_av, d_aw, s_m, DWIJ)

class pysph.sph.wc.gtvf.MomentumEquationArtificialStressSolid(dest, sources, dim)

Momentum equation Artificial stress for solids

See the class MomentumEquationPressureGradient for details.

Parameters:

dim (int) – Dimensionality of the problem.

loop(d_idx, s_idx, d_sigma, s_sigma, d_au, d_av, d_aw, s_m, DWIJ)

class pysph.sph.wc.gtvf.MomentumEquationPressureGradient(dest, sources, pref, gx=0.0, gy=0.0, gz=0.0)

Momentum Equation

From [ZhangHuAdams2017], equation (17),

\[\frac{\tilde{d} \boldsymbol{v}_i}{dt} = - \sum_j m_j \nabla W_{ij} \cdot \left[\left(\frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right)\textbf{I} - \left(\frac{\boldsymbol{A_i}}{\rho_i^2} + \frac{\boldsymbol{A_j}}{\rho_j^2} \right)\right] + \sum_j \frac{\eta_{ij}\boldsymbol{v}_{ij}}{\rho_i \rho_j r_{ij}} \nabla W_{ij} \cdot \boldsymbol{x}_{ij}\]

where,

\[\boldsymbol{A_{i/j}} = \rho_{i/j} \boldsymbol{v}_{i/j} \otimes (\tilde{\boldsymbol{v}}_{i/j} - \boldsymbol{v}_{i/j})\]

\[\eta_{ij} = \frac{2\eta_i \eta_j}{\eta_i + \eta_j}\]

\[\eta_{i/j} = \rho_{i/j} \nu\]

for solids, replace \(\boldsymbol{A}_{i/j}\) with \(\boldsymbol{\sigma}'_{i/j}\).

The rate of change of transport velocity is given by,

\[(\frac{d\boldsymbol{v}_i}{dt})_c = -p_i^0 \sum_j \frac{m_j} {\rho_i^2} \nabla \tilde{W}_{ij}\]

where,

\[\tilde{W}_{ij} = W(\boldsymbol{x}_ij, \tilde{0.5 h_{ij}})\]

\[p_i^0 = \min(10|p_i|, p_{ref})\]

Notes:

A negative sign in \((\frac{d\boldsymbol{v}_i}{dt})_c\) is missing in the paper [ZhangHuAdams2017].

Parameters:

• pref (float) – reference 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

initialize(d_idx, d_au, d_av, d_aw, d_auhat, d_avhat, d_awhat, d_p0, d_p)

loop(d_rho, s_rho, d_idx, s_idx, d_p, s_p, s_m, d_au, d_av, d_aw, DWIJ, d_p0, d_auhat, d_avhat, d_awhat, XIJ, RIJ, SPH_KERNEL, HIJ)

class pysph.sph.wc.gtvf.MomentumEquationViscosity(dest, sources, nu)

Momentum equation Artificial stress for solids

See the class MomentumEquationPressureGradient for details.

Notes:

A factor of ‘2’ is missing in the viscosity equation given by [ZhangHuAdams2017].

Parameters:

nu (float) – viscosity of the fluid.

loop(d_idx, s_idx, d_rho, s_rho, s_m, d_au, d_av, d_aw, VIJ, R2IJ, EPS, DWIJ, XIJ)

class pysph.sph.wc.gtvf.VelocityGradient(dest, sources, dim)

Gradient of velocity vector

\[(\nabla \otimes \tilde{\boldsymbol{v}})_i = \sum_j \frac{m_j} {\rho_j} \tilde{\boldsymbol{v}}_{ij} \otimes \nabla W_{ij}\]

Parameters:

dim (int) – Dimensionality of the problem.

initialize(d_idx, d_gradvhat)

loop(s_idx, d_idx, s_m, d_uhat, d_vhat, d_what, s_uhat, s_vhat, s_what, s_rho, d_gradvhat, DWIJ)

pysph.sph.wc.gtvf.get_particle_array_gtvf(constants=None, **props)

class pysph.sph.wc.density_correction.MLSFirstOrder2D(dest, sources)

Moving Least Squares density reinitialization This is a first order density reinitialization

\[W_{ab}^{MLS} = \beta\left(\mathbf{r_{a}}\right)\cdot\left( \mathbf{r}_a - \mathbf{r}_b\right)W_{ab}\]

\[\beta\left(\mathbf{r_{a}}\right) = A^{-1} \left[1 0 0\right]^{T}\]

where

\[A = \sum_{b}W_{ab}\tilde{A}\frac{m_{b}}{\rho_{b}}\]

\[\tilde{A} = pp^{T}\]

where

\[p = \left[1 x_{a}-x_{b} y_{a}-y_{b}\right]^{T}\]

\[\rho_{a} = \sum_{b} \m_{b}W_{ab}^{MLS}\]

References

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_rhotmp)

loop_all(d_idx, d_rho, d_x, d_y, s_x, s_y, d_h, s_h, s_m, s_rhotmp, SPH_KERNEL, NBRS, N_NBRS)

class pysph.sph.wc.density_correction.MLSFirstOrder3D(dest, sources)

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_rhotmp)

loop_all(d_idx, d_rho, d_x, d_y, d_z, s_x, s_y, s_z, d_h, s_h, s_m, s_rhotmp, SPH_KERNEL, NBRS, N_NBRS)

class pysph.sph.wc.density_correction.ShepardFilter(dest, sources)

Shepard Filter density reinitialization This is a zeroth order density reinitialization

\[\tilde{W_{ab}} = \frac{W_{ab}}{\sum_{b} W_{ab}\frac{m_{b}} {\rho_{b}}}\]

\[\rho_{a} = \sum_{b} \m_{b}\tilde{W_{ab}}\]

References

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_rhotmp)

loop_all(d_idx, d_rho, d_x, d_y, d_z, s_m, s_rhotmp, s_x, s_y, s_z, d_h, s_h, SPH_KERNEL, NBRS, N_NBRS)

Kernel Corrections

These are the equations for the kernel corrections that are mentioned in the paper by Bonet and Lok [BonetLok1999].

References

[BonetLok1999]

Bonet, J. and Lok T.-S.L. (1999) Variational and Momentum Preservation Aspects of Smoothed Particle Hydrodynamic Formulations.

class pysph.sph.wc.kernel_correction.GradientCorrection(dest, sources, dim=2, tol=0.1)

Kernel Gradient Correction

From [BonetLok1999], equations (42) and (45)

\[\nabla \tilde{W}_{ab} = L_{a}\nabla W_{ab}\]

\[L_{a} = \left(\sum \frac{m_{b}}{\rho_{b}} \nabla W_{ab} \mathbf{\otimes}x_{ba} \right)^{-1}\]

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_m_mat, DWIJ, HIJ)

class pysph.sph.wc.kernel_correction.GradientCorrectionPreStep(dest, sources, dim=2)

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_mat)

loop_all(d_idx, d_m_mat, s_m, s_rho, d_x, d_y, d_z, d_h, s_x, s_y, s_z, s_h, SPH_KERNEL, NBRS, N_NBRS)

class pysph.sph.wc.kernel_correction.KernelCorrection(dest, sources)

Kernel Correction

From [BonetLok1999], equation (53):

\[\mathbf{f}_{a} = \frac{\sum_{b}\frac{m_{b}}{\rho_{b}} \mathbf{f}_{b}W_{ab}}{\sum_{b}\frac{m_{b}}{\rho_{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_cwij)

loop(d_idx, s_idx, d_cwij, s_m, s_rho, WIJ)

class pysph.sph.wc.kernel_correction.MixedGradientCorrection(dest, sources, dim=2, tol=0.1)

Mixed Kernel Gradient Correction

This is as per [BonetLok1999]. See the MixedKernelCorrectionPreStep for the equations.

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_m_mat, d_dw_gamma, d_cwij, DWIJ, HIJ)

class pysph.sph.wc.kernel_correction.MixedKernelCorrectionPreStep(dest, sources, dim=2)

Mixed Kernel Correction

From [BonetLok1999], equations (54), (57) and (58)

\[\tilde{W}_{ab} = \frac{W_{ab}}{\sum_{b} V_{b}W_{ab}}\]

\[\nabla \tilde{W}_{ab} = L_{a}\nabla \bar{W}_{ab}\]

where,

\[L_{a} = \left(\sum_{b} V_{b} \nabla \bar{W}_{ab} \mathbf{\otimes}x_{ba} \right)^{-1}\]

\[\nabla \bar{W}_{ab} = \frac{\nabla W_{ab} - \gamma} {\sum_{b} V_{b}W_{ab}}\]

\[\gamma = \frac{\sum_{b} V_{b}\nabla W_{ab}} {\sum_{b} V_{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_m_mat)

loop_all(d_idx, d_x, d_y, d_z, d_h, s_x, s_y, s_z, s_h, SPH_KERNEL, N_NBRS, NBRS, d_m_mat, s_m, s_rho, d_cwij, d_dw_gamma)

CRKSPH corrections

These are equations for the basic kernel corrections in [CRKSPH2017].

References

[CRKSPH2017]

Nicholas Frontiere, Cody D. Raskin, J. Michael Owen “CRKSPH - A Conservative Reproducing Kernel Smoothed Particle Hydrodynamics Scheme”, Journal of Computational Physics 332 (2017), pp. 160–209.

class pysph.sph.wc.crksph.CRKSPH(dest, sources, dim=2, tol=0.5)

Conservative Reproducing Kernel SPH

Equations from the paper [CRKSPH2017].

\[W_{ij}^{R} = A_{i}\left(1+B_{i}^{\alpha}x_{ij}^{\alpha} \right)W_{ij}\]

\[\partial_{\gamma}W_{ij}^{R} = A_{i}\left(1+B_{i}^{\alpha} x_{ij}^{\alpha}\right)\partial_{\gamma}W_{ij} + \partial_{\gamma}A_{i}\left(1+B_{i}^{\alpha}x_{ij}^{\alpha} \right)W_{ij} + A_{i}\left(\partial_{\gamma}B_{i}^{\alpha} x_{ij}^{\alpha} + B_{i}^{\gamma}\right)W_{ij}\]

\[\nabla\tilde{W}_{ij} = 0.5 * \left(\nabla W_{ij}^{R}-\nabla W_{ji}^{R} \right)\]

where,

\[A_{i} = \left[m_{0} - \left(m_{2}^{-1}\right)^{\alpha \beta} m_1^{\beta}m_1^{\alpha}\right]^{-1}\]

\[B_{i}^{\alpha} = -\left(m_{2}^{-1}\right)^{\alpha \beta} m_{1}^{\beta}\]

\[\partial_{\gamma}A_{i} = -A_{i}^{2}\left(\partial_{\gamma} m_{0}-\left(m_{2}^{-1}\right)^{\alpha \beta}\left( m_{1}^{\beta}\partial_{\gamma}m_{1}^{\alpha} + \partial_{\gamma}m_{1}^{\beta}m_{1}^{\alpha}\right) + \left(m_{2}^{-1}\right)^{\alpha \phi}\partial_{\gamma} m_{2}^{\phi \psi}\left(m_{2}^{-1}\right)^{\psi \beta} m_{1}^{\beta}m_{1}^{\alpha} \right)\]

\[\partial_{\gamma}B_{i}^{\alpha} = -\left(m_{2}^{-1}\right)^ {\alpha \beta}\partial_{\gamma}m_{1}^{\beta} + \left(m_{2}^{-1}\right)^ {\alpha \phi}\partial_{\gamma}m_{2}^{\phi \psi}\left(m_{2}^ {-1}\right)^{\psi \beta}m_{1}^{\beta}\]

\[m_{0} = \sum_{j}V_{j}W_{ij}\]

\[m_{1}^{\alpha} = \sum_{j}V_{j}x_{ij}^{\alpha}W_{ij}\]

\[m_{2}^{\alpha \beta} = \sum_{j}V_{j}x_{ij}^{\alpha} x_{ij}^{\beta}W_{ij}\]

\[\partial_{\gamma}m_{0} = \sum_{j}V_{j}\partial_{\gamma} W_{ij}\]

\[\partial_{\gamma}m_{1}^{\alpha} = \sum_{j}V_{j}\left[ x_{ij}^{\alpha}\partial_{\gamma}W_{ij}+\delta^ {\alpha \gamma}W_{ij} \right]\]

\[\partial_{\gamma}m_{2}^{\alpha \beta} = \sum_{j}V_{j}\left[ x_{ij}^{\alpha}x_{ij}^{\beta}\partial_{\gamma}W_{ij} + \left(x_{ij}^{\alpha}\delta^{\beta \gamma} + x_{ij}^{\beta} \delta^{\alpha \gamma} \right)W_{ij} \right]\]

Parameters:

• dim (int) – Dimensionality of the problem.

• tol (float) – Tolerence value to decide std or corrected kernel

loop(d_idx, s_idx, d_ai, d_gradai, d_cwij, d_bi, d_gradbi, WIJ, DWIJ, XIJ, HIJ)

class pysph.sph.wc.crksph.CRKSPHIntegrator(**kw)

Pass fluid names and suitable IntegratorStep instances.

For example:

>>> integrator = Integrator(fluid=WCSPHStep(), solid=WCSPHStep())

where “fluid” and “solid” are the names of the particle arrays.

one_timestep(t, dt)

User written function that actually does one timestep.

This function is used in the high-performance Cython implementation. The assumptions one may make are the following:

• t and dt are passed.

• the following methods are available:

  • self.initialize()

  • self.stage1(), self.stage2() etc. depending on the number of stages available.

  • self.compute_accelerations(index=0, update_nnps=True)

  • self.do_post_stage(stage_dt, stage_count_from_1)

  • self.update_domain()

Please see any of the concrete implementations of the Integrator class to study. By default the Integrator implements a predict-evaluate-correct method, the same as PECIntegrator.

class pysph.sph.wc.crksph.CRKSPHPreStep(dest, sources, dim=2)

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

loop_all(d_idx, d_x, d_y, d_z, d_h, s_x, s_y, s_z, s_h, s_m, s_rho, SPH_KERNEL, NBRS, N_NBRS, d_ai, d_gradai, d_bi, s_V, d_gradbi)

class pysph.sph.wc.crksph.CRKSPHScheme(fluids, dim, rho0, c0, nu, h0, p0, gx=0.0, gy=0.0, gz=0.0, cl=2, cq=1, gamma=7.0, eta_crit=0.3, eta_fold=0.2, tol=0.5, has_ghosts=False)

Parameters:

• fluids (list) – a list with names of fluid particle arrays

• solids (list) – a list with names of solid (or boundary) particle arrays

configure_solver(kernel=None, integrator_cls=None, extra_steppers=None, **kw)

Configure the solver to be generated.

Parameters:

• kernel (Kernel instance.) – Kernel to use, if none is passed a default one is used.

• integrator_cls (pysph.sph.integrator.Integrator) – Integrator class to use, use sensible default if none is passed.

• extra_steppers (dict) – Additional integration stepper instances as a dict.

• **kw (extra arguments) – Any additional keyword args are passed to the solver instance.

get_equations()

setup_properties(particles, clean=True)

Setup the particle arrays so they have the right set of properties for this scheme.

Parameters:

• particles (list) – List of particle arrays.

• clean (bool) – If True, removes any unnecessary properties.

class pysph.sph.wc.crksph.CRKSPHStep

stage1(d_idx, d_u, d_v, d_w, d_u0, d_v0, d_w0)

stage2(d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, dt)

stage3(d_idx, d_e, d_ae, d_u, d_v, d_w, d_u0, d_v0, d_w0, d_x, d_y, d_z, dt)

class pysph.sph.wc.crksph.CRKSPHSymmetric(dest, sources, dim=2, tol=0.5)

Conservative Reproducing Kernel SPH

This is symmetric and will only work for particles of the same array.

Equations from the paper [CRKSPH2017].

\[W_{ij}^{R} = A_{i}\left(1+B_{i}^{\alpha}x_{ij}^{\alpha} \right)W_{ij}\]

\[\partial_{\gamma}W_{ij}^{R} = A_{i}\left(1+B_{i}^{\alpha} x_{ij}^{\alpha}\right)\partial_{\gamma}W_{ij} + \partial_{\gamma}A_{i}\left(1+B_{i}^{\alpha}x_{ij}^{\alpha} \right)W_{ij} + A_{i}\left(\partial_{\gamma}B_{i}^{\alpha} x_{ij}^{\alpha} + B_{i}^{\gamma}\right)W_{ij}\]

\[\nabla\tilde{W}_{ij} = 0.5 * \left(\nabla W_{ij}^{R}-\nabla W_{ji}^{R} \right)\]

where,

\[A_{i} = \left[m_{0} - \left(m_{2}^{-1}\right)^{\alpha \beta} m1_{\beta}m1_{\alpha}\right]^{-1}\]

\[B_{i}^{\alpha} = -\left(m_{2}^{-1}\right)^{\alpha \beta} m_{1}^{\beta}\]

\[\partial_{\gamma}A_{i} = -A_{i}^{2}\left(\partial_{\gamma} m_{0}-\left[m_{2}^{-1}\right]^{\alpha \beta}\left[ m_{1}^{\beta}\partial_{\gamma}m_{1}^{\beta}m_{1}^{\alpha} + \partial_{\gamma}m_{1}^{\alpha}m_{1}^{\beta}\right] + \left[m_{2}^{-1}\right]^{\alpha \phi}\partial_{\gamma} m_{2}^{\phi \psi}\left[m_{2}^{-1}\right]^{\psi \beta} m_{1}^{\beta}m_{1}^{\alpha} \right)\]

\[\partial_{\gamma}B_{i}^{\alpha} = -\left(m_{2}^{-1}\right)^ {\alpha \beta}\partial_{\gamma}m_{1}^{\beta} + \left(m_{2}^{-1}\right)^ {\alpha \phi}\partial_{\gamma}m_{2}^{\phi \psi}\left(m_{2}^ {-1}\right)^{\psi \beta}m_{1}^{\beta}\]

\[m_{0} = \sum_{j}V_{j}W_{ij}\]

\[m_{1}^{\alpha} = \sum_{j}V_{j}x_{ij}^{\alpha}W_{ij}\]

\[m_{2}^{\alpha \beta} = \sum_{j}V_{j}x_{ij}^{\alpha} x_{ij}^{\beta}W_{ij}\]

\[\partial_{\gamma}m_{0} = \sum_{j}V_{j}\partial_{\gamma} W_{ij}\]

\[\partial_{\gamma}m_{1}^{\alpha} = \sum_{j}V_{j}\left[ x_{ij}^{\alpha}\partial_{\gamma}W_{ij}+\delta^ {\alpha \gamma}W_{ij} \right]\]

\[\partial_{\gamma}m_{2}^{\alpha \beta} = \sum_{j}V_{j}\left[ x_{ij}^{\alpha}x_{ij}^{\beta}\partial_{\gamma}W_{ij} + \left(x_{ij}^{\alpha}\delta^{\beta \gamma} + x_{ij}^{\beta} \delta^{\alpha \gamma} \right)W_{ij} \right]\]

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

loop(d_idx, s_idx, d_ai, d_gradai, d_cwij, d_bi, d_gradbi, s_ai, s_gradai, s_bi, s_gradbi, d_h, s_h, WIJ, DWIJ, XIJ, HIJ, RIJ, DWI, DWJ, WI, WJ, SPH_KERNEL)

class pysph.sph.wc.crksph.CRKSPHUpdateGhostProps(dest, sources=None, dim=2)

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_orig_idx, d_p, d_cs, d_V, d_ai, d_bi, d_tag, d_gradbi, d_gradai, d_rho, d_u0, d_v0, d_w0, d_gradv, d_u, d_v, d_w)

class pysph.sph.wc.crksph.EnergyEquation(dest, sources, dim, gamma, gx=0.0, gy=0.0, gz=0.0, cl=2, cq=1, eta_crit=0.5, eta_fold=0.2, tol=0.5)

Energy Equation

From [CRKSPH2017], equation (66):

\[\Delta u_{ij} = \frac{f_{ij}}{2}\left(v_j^{\alpha}(t) + v_j^{\alpha}(t + \Delta t) - v_i^{\alpha}(t) - v_i^{\alpha}(t + \Delta t)\right) \frac{Dv_{ij}^{\alpha}}{Dt}\]

\[\begin{split}f_{ij} = \begin{cases} 1/2 &|s_i - s_j| = 0,\\ s_{\min} / (s_{\min} + s_{\max}) &\Delta u_{ij}\times(s_i - s_j) > 0\\ s_{\max} / (s_{\min} + s_{\max}) &\Delta u_{ij}\times(s_i - s_j) < 0\\ \end{cases}\end{split}\]

\[s_{\min} = \min(|s_i|, |s_j|)\]

\[s_{\max} = \max(|s_i|, |s_j|)\]

\[s_{i/j} = \frac{p_{i/j}}{\rho_{i/j}^\gamma}\]

see MomentumEquation for \(\frac{Dv_{ij}^{\alpha}}{Dt}\)

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_ae)

loop(d_idx, s_idx, d_au, d_av, d_aw, d_ae, s_au, s_av, s_aw, d_u0, d_v0, d_w0, s_u0, s_v0, s_w0, d_u, d_v, d_w, s_u, s_v, s_w, d_p, d_rho, s_p, s_rho, d_m, d_V, s_V, d_cs, s_cs, d_h, s_h, XIJ, d_gradv, s_gradv, EPS, DWIJ)

class pysph.sph.wc.crksph.MomentumEquation(dest, sources, dim, gx=0.0, gy=0.0, gz=0.0, cl=2, cq=1, eta_crit=0.3, eta_fold=0.2, tol=0.5)

Momentum Equation

From [CRKSPH2017], equation (64):

\[\frac{Dv_{i}^{\alpha}}{Dt} = -\frac{1}{2 m_i}\sum_{j} V_i V_j (P_i + P_j + Q_i + Q_j) (\partial_{\alpha}W_{ij}^R - \partial_{\alpha} W_{ji}^R)\]

where,

\[V_{i/j} = \text{dest/source particle number density}\]

\[P_{i/j} = \text{dest/source particle pressure}\]

\[Q_i = \rho_{i} (-C_{l} c_{i} \mu_{i} + C_{q} \mu_{i}^{2})\]

\[\mu_i = \min \left(0, \frac{\hat{v}_{ij} \eta_{i}^{\alpha}}{\eta_{i}^{\alpha}\eta_{i}^{\alpha} + \epsilon^{2}}\right)\]

\[\hat{v}_{ij}^{\alpha} = v_{i}^{\alpha} - v_{j}^{\alpha} - \frac{\phi_{ij}}{2}\left(\partial_{\beta} v_i^{\alpha} + \partial_{\beta}v_j^{\alpha}\right) x_{ij}^{\beta}\]

\[\begin{split}\phi_{ij} = \max \left[0, \min \left[1, \frac{4r_{ij}}{(1 + r_{ij})^2}\right]\right] \times \begin{cases} \exp{\left[-\left((\eta_{ij} - \eta_{crit})/\eta_{fold}\right)^2\right]}, &\eta_{ij} < \eta_{crit} \\ 1, & \eta_{ij} >= \eta_{crit} \end{cases}\end{split}\]

\[\eta_{ij} = \min(\eta_i, \eta_j)\]

\[\eta_{i/j} = (x_{ij}^{\alpha} x_{ij}^{\alpha})^{1/2} / h_{i/j}\]

\[r_{ij} = \frac{\partial_{\beta} v_i^{\alpha} x_{ij}^{\alpha} x_{ij}^{\beta}}{\partial_{\beta} v_j^{\alpha}x_{ij}^{\alpha} x_{ij}^{\beta}}\]

\[\partial_{\beta} v_i^{\alpha} = -\sum_j V_j v_{ij}^{\alpha} \partial_{\beta} W_{ij}^R\]

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)

loop(d_idx, s_idx, d_m, s_m, d_rho, s_rho, d_p, s_p, d_cs, s_cs, d_u, d_v, d_w, s_u, s_v, s_w, d_gradv, s_gradv, d_h, s_h, s_ai, s_bi, s_gradai, s_gradbi, d_au, d_av, d_aw, d_V, s_V, XIJ, VIJ, DWIJ, EPS, HIJ, WIJ, DWI, DWJ)

class pysph.sph.wc.crksph.NumberDensity(dest, sources)

Number Density

From [CRKSPH2017], equation (75):

\[V_{i}^{-1} = \sum_{j} W_{i}\]

Note that the quantity V is the inverse of particle volume, so when using in the equation use \(\frac{1}{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_V)

loop(d_idx, d_V, WI)

class pysph.sph.wc.crksph.SpeedOfSound(dest, sources=None, gamma=7.0)

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

initialize(d_cs, d_idx, d_p, d_rho)

class pysph.sph.wc.crksph.StateEquation(dest, sources, gamma)

State Equation

State equation for ideal gas, from [CRKSPH2017] equation (77):

\[p_i = (\gamma - 1)\rho_{i} u_i\]

where, \(u_i\) is the specific thermal energy.

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_p, d_rho, d_e)

class pysph.sph.wc.crksph.SummationDensityCRKSPH(dest, sources)

Summation Density CRKSPH

From [CRKSPH2017], equation (76):

\[\rho_{i} = \frac{\sum_j m_{ij} V_j W_{ij}^R} {\sum_{j} V_{j}^{2} W_{ij}^R}\]

where,

\[\begin{split}mij = \begin{cases} m_j, &i \text{ and } j \text{ are same material} \\ m_i, &i \text{ and } j \text{ are different material} \end{cases}\end{split}\]

Note that in this equation we are just using \(m_{ij}\) to be \(m_i\) as the mass remains constant through out the simulation.

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_rhofac)

loop(d_idx, s_idx, d_m, d_rho, d_rhofac, s_V, WIJ, d_cwij)

post_loop(d_idx, d_rho, d_rhofac)

class pysph.sph.wc.crksph.VelocityGradient(dest, sources, dim)

Velocity Gradient

From [CRKSPH2017], equation (74)

\[\partial_{\beta} v_i^{\alpha} = -\sum_j V_j v_{ij}^{\alpha} \partial_{\beta} W_{ij}^R\]

Parameters:

dim (int) – Dimensionality of the Problem.

initialize(d_idx, d_gradv)

loop(d_idx, s_idx, s_V, d_gradv, d_h, s_h, XIJ, DWIJ, VIJ, DWI)

pysph.sph.wc.crksph.get_particle_array_crksph(constants=None, **props)

Predictive-Corrective Incompressible SPH (PCISPH)

References

[SolPaj2009]

B. Solenthaler, R. Pajarola “Predictive-Corrective Incompressible SPH”, ACM Trans. Graph 28 (2009), pp. 1–6.

class pysph.sph.wc.pcisph.ComputePressure(dest, sources, rho0)

Compute Pressure

Compute pressure iteratively maintaining density within a given tolerance.

\[p_i += \delta \rho^{*}_{{err}_i}\]

where,

\[\rho_{err_i} = \rho_i^{*} - \rho_0\]

\[\delta = \frac{-1}{\beta (-\sum_j \nabla W_{ij} \cdot \sum_j \nabla W_{ij} - \sum_j \nabla W_{ij} \nabla W_{ij})}\]

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_dw, d_dwij2)

loop(d_idx, d_dw, d_dwij2, DWIJ)

post_loop(d_idx, d_dw, d_m, dt, d_dwij2, d_p, d_rho)

class pysph.sph.wc.pcisph.MomentumEquationPressureGradient(dest, sources, rho0, tolerance, debug)

Momentum Equation pressure gradient

Standard WCSPH pressure gradient,

\[\frac{d\mathbf{v}}{dt} = - \sum_j m_j \left(\frac{p_i}{\rho_i^2} + \frac{p_i}{\rho_i^2}\right) \nabla W(x_{ij}, h)\]

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

converged()

Return > 0 to indicate converged iterations and < 0 otherwise.

loop(d_idx, s_idx, d_p, s_p, d_rho, s_rho, s_m, d_aup, d_avp, d_awp, DWIJ)

reduce(dst, t, dt)

class pysph.sph.wc.pcisph.MomentumEquationViscosity(dest, sources, nu=0.0, gx=0.0, gy=0.0, gz=0.0)

Momentum Equation Viscosity

See “pysph.sph.wc.viscosity.LaminarViscocity”

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)

loop(d_idx, s_idx, s_m, d_rho, s_rho, d_au, d_av, d_aw, DWIJ, XIJ, VIJ, R2IJ, EPS)

post_loop(d_idx, d_au, d_av, d_aw, d_u, d_v, d_w, d_p, d_aup, d_avp, d_awp, dt)

class pysph.sph.wc.pcisph.PCISPHIntegrator(**kw)

Pass fluid names and suitable IntegratorStep instances.

For example:

>>> integrator = Integrator(fluid=WCSPHStep(), solid=WCSPHStep())

where “fluid” and “solid” are the names of the particle arrays.

initial_acceleration(t, dt)

Compute the initial accelerations if needed before the iterations start.

The default implementation only does this for the first acceleration evaluator. So if you have multiple evaluators, you must override this method in a subclass.

one_timestep(t, dt)

User written function that actually does one timestep.

This function is used in the high-performance Cython implementation. The assumptions one may make are the following:

• t and dt are passed.

• the following methods are available:

  • self.initialize()

  • self.stage1(), self.stage2() etc. depending on the number of stages available.

  • self.compute_accelerations(index=0, update_nnps=True)

  • self.do_post_stage(stage_dt, stage_count_from_1)

  • self.update_domain()

Please see any of the concrete implementations of the Integrator class to study. By default the Integrator implements a predict-evaluate-correct method, the same as PECIntegrator.

class pysph.sph.wc.pcisph.PCISPHScheme(fluids, dim, rho0, nu, gx=0.0, gy=0.0, gz=0.0, tolerance=0.1, debug=False, show_itercount=False)

Parameters:

• fluids (list) – List of names of fluid particle arrays.

• solids (list) – List of names of solid particle arrays (or boundaries).

• dim (int) – Dimensionality of the problem.

add_user_options(group)

configure_solver(kernel=None, integrator_cls=None, extra_steppers=None, **kw)

Configure the solver to be generated.

Parameters:

• kernel (Kernel instance.) – Kernel to use, if none is passed a default one is used.

• integrator_cls (pysph.sph.integrator.Integrator) – Integrator class to use, use sensible default if none is passed.

• extra_steppers (dict) – Additional integration stepper instances as a dict.

• **kw (extra arguments) – Any additional keyword args are passed to the solver instance.

consume_user_options(options)

get_equations()

setup_properties(particles, clean=True)

Setup the particle arrays so they have the right set of properties for this scheme.

Parameters:

• particles (list) – List of particle arrays.

• clean (bool) – If True, removes any unnecessary properties.

class pysph.sph.wc.pcisph.PCISPHStep(show_itercount=False)

initialize(d_idx, d_u, d_v, d_w, d_u0, d_v0, d_w0, d_x, d_y, d_z, d_x0, d_y0, d_z0, d_rho, d_rho0)

py_stage1(dst, t, dt)

stage1(d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, d_x, d_y, d_z, d_aup, d_avp, d_awp, d_u0, d_v0, d_w0, d_x0, d_y0, d_z0, dt)

class pysph.sph.wc.pcisph.Predict(dest, sources)

Predict velocity and position

\[\mathbf{v}^{*}(t+1) = \mathbf{v}(t) + dt \left(\frac{d \mathbf{v}_{visc, g}(t)}{dt} + \frac{d \mathbf{v}_{p} (t)}{dt} \right)\]

\[\mathbf{x}^{*}(t+1) = \mathbf{x}(t) + dt * \mathbf{v}(t+1)\]

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_u, d_v, d_w, d_aup, d_avp, d_awp, d_x, d_y, d_z, d_au, d_av, d_aw, d_u0, d_v0, d_w0, d_x0, d_y0, d_z0, dt)

pysph.sph.wc.pcisph.get_particle_array_pcisph(constants=None, **props)

SPH Boundary Equations

class pysph.sph.boundary_equations.MonaghanBoundaryForce(dest, sources, deltap)

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

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)

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

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

loop(d_idx, s_idx, d_m, s_m, d_au, d_av, d_aw, RIJ, R2IJ, XIJ)

pysph.sph.boundary_equations.wendland_quintic(rij=1.0, h=1.0)

Basic Equations for Solid Mechanics

References

[Gray2001]

J. P. Gray et al., “SPH elastic dynamics”, Computer Methods in Applied Mechanics and Engineering, 190 (2001), pp 6641 - 6662.

class pysph.sph.solid_mech.basic.ElasticSolidsScheme(elastic_solids, solids, dim, artificial_stress_eps=0.3, xsph_eps=0.5, alpha=1.0, beta=1.0)

Parameters:

• fluids (list) – List of names of fluid particle arrays.

• solids (list) – List of names of solid particle arrays (or boundaries).

• dim (int) – Dimensionality of the problem.

configure_solver(kernel=None, integrator_cls=None, extra_steppers=None, **kw)

Configure the solver to be generated.

Parameters:

• kernel (Kernel instance.) – Kernel to use, if none is passed a default one is used.

• integrator_cls (pysph.sph.integrator.Integrator) – Integrator class to use, use sensible default if none is passed.

• extra_steppers (dict) – Additional integration stepper instances as a dict.

• **kw (extra arguments) – Any additional keyword args are passed to the solver instance.

get_equations()

class pysph.sph.solid_mech.basic.EnergyEquationWithStress(dest, sources, alpha=1.0, beta=1.0, eta=0.01)

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_ae)

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)

post_loop(d_idx, d_rho, d_s00, d_s01, d_s02, d_s11, d_s12, d_s22, d_v00, d_v01, d_v02, d_v10, d_v11, d_v12, d_v20, d_v21, d_v22, d_ae)

class pysph.sph.solid_mech.basic.HookesDeviatoricStressRate(dest, sources)

Rate of change of stress

\[\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{align}\begin{aligned}\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)\end{aligned}\end{align} \]

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_as00, d_as01, d_as02, d_as11, d_as12, d_as22)

loop(d_idx, d_s00, d_s01, d_s02, d_s11, d_s12, d_s22, d_v00, d_v01, d_v02, d_v10, d_v11, d_v12, d_v20, d_v21, d_v22, d_as00, d_as01, d_as02, d_as11, d_as12, d_as22, d_G)

class pysph.sph.solid_mech.basic.IsothermalEOS(dest, sources)

Compute the pressure using the Isothermal equation of state:

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

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_rho, d_p, d_c0_ref, d_rho_ref)

class pysph.sph.solid_mech.basic.MomentumEquationWithStress(dest, sources)

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{align}\begin{aligned}\begin{split}f_{ab} = \frac{W(r_{ab})}{W(\Delta p)}\\\end{split}\\R_{ab}^{ij} = R_{a}^{ij} + R_{b}^{ij}\end{aligned}\end{align} \]

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)

loop(d_idx, s_idx, d_rho, s_rho, s_m, d_p, s_p, d_s00, d_s01, d_s02, d_s11, d_s12, d_s22, s_s00, s_s01, s_s02, s_s11, s_s12, s_s22, d_r00, d_r01, d_r02, d_r11, d_r12, d_r22, s_r00, s_r01, s_r02, s_r11, s_r12, s_r22, d_au, d_av, d_aw, d_wdeltap, d_n, WIJ, DWIJ)

class pysph.sph.solid_mech.basic.MonaghanArtificialStress(dest, sources, eps=0.3)

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{align}\begin{aligned}\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}\end{aligned}\end{align} \]

Components of \(R\) in rotated frame:

\[ \begin{align}\begin{aligned}\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}\end{aligned}\end{align} \]

Components of \(R\) in original frame:

\[ \begin{align}\begin{aligned}\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)\end{aligned}\end{align} \]

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)

Compute the stress terms

Parameters:

• d_sxx (DoubleArray) – Stress Tensor Deviatoric components (Symmetric)

• d_rxx (DoubleArray) – Artificial stress components (Symmetric)

pysph.sph.solid_mech.basic.get_bulk_mod(G, nu)

Get the bulk modulus from shear modulus and Poisson ratio

pysph.sph.solid_mech.basic.get_particle_array_elastic_dynamics(constants=None, **props)

Return a particle array for the Standard SPH formulation of solids.

Parameters:

• constants (dict) – Dictionary of constants

• props (dict) – Additional keywords passed are set as the property arrays.

See also

get_particle_array

pysph.sph.solid_mech.basic.get_shear_modulus(E, nu)

pysph.sph.solid_mech.basic.get_speed_of_sound(E, nu, rho0)

Equations for the High Velocity Impact Problems

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

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_p, d_rho, d_e)

class pysph.sph.solid_mech.hvi.StiffenedGasEOS(dest, sources, gamma, r0, c0)

Stiffened-gas EOS from “A Free Lagrange Augmented Godunov Method for the Simulation of Elastic-Plastic Solids”, B. P. Howell and G. J. Ball, JCP (2002). http://dx.doi.org/10.1006/jcph.2001.6931

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_e, d_rho, d_p, d_cs)

class pysph.sph.solid_mech.hvi.VonMisesPlasticity2D(dest, sources, flow_stress)

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_s00, d_s01, d_s02, d_s11, d_s12, d_s22)

Gas Dynamics

Basic equations for Gas-dynamics

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

ADKE as discussed in the reference [KP14].

References

[KP14] (1,2,3)

A comparison of SPH schemes for the compressible Euler equations, 2014, Journal of Computational Physics, 256, pp 308 – 333 (http://dx.doi.org/10.1016/j.jcp.2013.08.060)

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, d_ae)

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)

class pysph.sph.gas_dynamics.basic.ADKEUpdateGhostProps(dest, sources=None, dim=2)

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_orig_idx, d_p, d_cs, d_tag, d_rho)

class pysph.sph.gas_dynamics.basic.IdealGasEOS(dest, sources, gamma)

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_p, d_rho, d_e, d_cs)

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)

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, d_ae, d_am, d_aalpha1, d_aalpha2, d_del2e, d_dt_cfl)

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, d_dt_cfl)

post_loop(d_idx, d_h, d_cs, d_alpha1, d_aalpha1, d_div, d_del2e, d_e, d_alpha2, d_aalpha2)

class pysph.sph.gas_dynamics.basic.MPMUpdateGhostProps(dest, sources=None, dim=2)

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_orig_idx, d_p, d_cs, d_tag)

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

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, d_ae)

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)

class pysph.sph.gas_dynamics.basic.ScaleSmoothingLength(dest, sources, factor=2.0)

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_h)

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

Summation density with iterative solution of the smoothing lengths.

Parameters:

density_iterationsbint

Flag to indicate density iterations are required.

iterate_only_oncebint

Flag to indicate if only one iteration is required

kdouble

Kernel scaling factor

htoldouble

Iteration tolerance

converged()

Return > 0 to indicate converged iterations and < 0 otherwise.

initialize(d_idx, d_rho, d_div, d_grhox, d_grhoy, d_grhoz, d_arho, d_dwdh)

loop(d_idx, s_idx, d_rho, d_grhox, d_grhoy, d_grhoz, d_arho, d_dwdh, s_m, d_converged, VIJ, WI, DWI, GHI)

post_loop(d_idx, d_arho, d_rho, d_div, d_omega, d_dwdh, d_h0, d_h, d_m, d_ah, d_converged)

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

References

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, d_rho, d_h, d_h0)

loop(d_idx, d_rho, d_arho, s_idx, s_m, VIJ, DWI, WIJ)

post_loop(d_idx, d_rho, d_arho, d_div, d_logrho)

reduce(dst, t, dt)

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

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_m, d_rho, d_h)

Boundary equations for Gas-dynamics

class pysph.sph.gas_dynamics.boundary_equations.WallBoundary(dest, sources)

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_p, d_rho, d_e, d_m, d_cs, d_div, d_h, d_htmp, d_h0, d_u, d_v, d_w, d_wij)

loop(d_idx, s_idx, d_p, d_rho, d_e, d_m, d_cs, d_div, d_h, d_u, d_v, d_w, d_wij, d_htmp, s_p, s_rho, s_e, s_m, s_cs, s_h, s_div, s_u, s_v, s_w, WI)

post_loop(d_idx, d_p, d_rho, d_e, d_m, d_cs, d_div, d_h, d_u, d_v, d_w, d_wij, d_htmp)

Surface tension

Implementation of the equations used for surface tension modelling, for example in KHI simulations. The references are [SY11], [JM00], [A10], [XA06].

References

[SY11]

M. Shadloo, M. Yildiz, “Numerical modelling of Kelvin-Helmholtz instability using smoothed particle hydrodynamics”, IJNME, 2011, 87, pp 988–1006

[JM00]

Joseph P. Morris “Simulating surface tension with smoothed particle hydrodynamics”, JCP, 2000, 33, pp 333–353

[A10]

Adami et al. “A new surface-tension formulation for multi-phase SPH using a reproducing divergence approximation”, JCP 2010, 229, pp 5011–5021

[XA06]

X.Y.Hu, N.A. Adams. “A multi-phase SPH method for macroscopic and mesoscopic flows”, JCP 2006, 213, pp 844-861 [XA06]

class pysph.sph.surface_tension.AdamiColorGradient(dest, sources)

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)

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

post_loop(d_idx, d_cx, d_cy, d_cz, d_h, d_nx, d_ny, d_nz, d_ddelta, d_N)

class pysph.sph.surface_tension.AdamiReproducingDivergence(dest, sources, dim)

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}\]

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_kappa, d_wij_sum)

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, d_N, s_N, d_color, s_color)

post_loop(d_idx, d_kappa, d_wij_sum)

class pysph.sph.surface_tension.CSFSurfaceTensionForce(dest, sources, sigma=0.1)

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.

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)

loop(d_idx, d_au, d_av, d_aw, d_kappa, d_cx, d_cy, d_cz, d_rho)

class pysph.sph.surface_tension.CSFSurfaceTensionForceAdami(dest, sources, sigma)

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)

post_loop(d_idx, d_au, d_av, d_aw, d_kappa, d_cx, d_cy, d_cz, d_m, d_alpha, d_rho)

class pysph.sph.surface_tension.ColorGradientAdami(dest, sources)

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)

loop(d_idx, d_cx, d_cy, d_cz, d_V, s_V, d_color, s_color, DWIJ, s_idx)

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

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\)

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)

loop(d_idx, s_idx, d_scolor, s_scolor, d_cx, d_cy, d_cz, d_V, s_V, DWIJ)

post_loop(d_idx, d_cx, d_cy, d_cz, d_nx, d_ny, d_nz, d_N, d_ddelta)

class pysph.sph.surface_tension.ConstructStressMatrix(dest, sources, sigma, d=2)

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

initialize(d_pi00, d_pi01, d_pi02, d_pi10, d_pi11, d_pi12, d_pi20, d_pi21, d_pi22, d_cx, d_cy, d_cz, d_idx, d_N)

class pysph.sph.surface_tension.InterfaceCurvatureFromDensity(dest, sources, with_morris_correction=True)

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_kappa, d_wij_sum)

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)

post_loop(d_idx, d_wij_sum, d_nx, d_kappa)

class pysph.sph.surface_tension.InterfaceCurvatureFromNumberDensity(dest, sources, with_morris_correction=True)

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}\]

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_kappa, d_wij_sum)

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)

post_loop(d_idx, d_wij_sum, d_nx, d_kappa)

class pysph.sph.surface_tension.MomentumEquationPressureGradientAdami(dest, sources, gx=0.0, gy=0.0, gz=0.0)

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)

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)

post_loop(d_idx, d_au, d_av, d_aw)

class pysph.sph.surface_tension.MomentumEquationPressureGradientHuAdams(dest, sources, gx=0.0, gy=0.0, gz=0.0)

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

initialize(d_au, d_av, d_aw, d_idx)

loop(d_idx, d_V, d_au, d_av, d_aw, s_V, d_p, s_p, DWIJ, s_idx, d_m)

post_loop(d_idx, d_au, d_av, d_aw)

class pysph.sph.surface_tension.MomentumEquationPressureGradientMorris(dest, sources)

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)

loop(d_idx, s_idx, d_au, d_av, d_aw, s_m, d_p, s_p, DWIJ, d_rho, s_rho)

class pysph.sph.surface_tension.MomentumEquationViscosityAdami(dest, sources)

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

initialize(d_au, d_av, d_aw, d_idx)

loop(d_idx, d_V, d_au, d_av, d_aw, s_V, d_p, s_p, DWIJ, s_idx, d_m, R2IJ, XIJ, EPS, VIJ, d_nu, s_nu)

class pysph.sph.surface_tension.MomentumEquationViscosityMorris(dest, sources, eta=0.01)

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

loop(d_idx, s_idx, d_au, d_av, d_aw, s_m, d_nu, s_nu, d_rho, s_rho, DWIJ, R2IJ, VIJ, HIJ, XIJ)

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

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\)

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)

loop(d_idx, s_idx, d_scolor, s_scolor, d_cx, d_cy, d_cz, s_m, s_rho, DWIJ)

post_loop(d_idx, d_cx, d_cy, d_cz, d_nx, d_ny, d_nz, d_N, d_ddelta)

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

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.

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)

loop(d_idx, s_idx, d_color, s_color, d_cx, d_cy, d_cz, d_V, s_V, DWIJ)

post_loop(d_idx, d_cx, d_cy, d_cz, d_nx, d_ny, d_nz, d_N, d_ddelta)

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

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.

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_ddelta)

loop(d_idx, s_idx, d_color, s_color, d_cx, d_cy, d_cz, d_V, s_V, DWIJ)

post_loop(d_idx, d_cx, d_cy, d_cz, d_nx, d_ny, d_nz, d_N, d_ddelta)

class pysph.sph.surface_tension.ShadlooViscosity(dest, sources, alpha)

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)

loop(d_idx, d_au, d_av, d_aw, d_h, s_idx, s_h, d_cs, s_cs, d_rho, s_rho, VIJ, XIJ, d_V, s_V, R2IJ, EPS, DWIJ)

class pysph.sph.surface_tension.ShadlooYildizSurfaceTensionForce(dest, sources, sigma=0.1)

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.

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)

loop(d_idx, d_au, d_av, d_aw, d_kappa, d_nx, d_ny, d_nz, d_m, d_rho, d_ddelta)

class pysph.sph.surface_tension.SmoothedColor(dest, sources)

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.

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_scolor)

loop(d_idx, s_idx, d_rho, s_rho, s_m, s_color, d_scolor, WIJ)

class pysph.sph.surface_tension.SolidWallPressureBCnoDensity(dest, sources)

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_p, d_wij)

loop(d_idx, s_idx, d_p, s_p, d_wij, s_rho, d_au, d_av, d_aw, WIJ, XIJ)

post_loop(d_idx, d_wij, d_p, d_rho)

class pysph.sph.surface_tension.SummationDensitySourceMass(dest, sources)

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)

loop(d_idx, d_V, d_rho, d_m, WIJ, s_m, s_idx)

post_loop(d_idx, d_V, d_rho, d_m)

class pysph.sph.surface_tension.SurfaceForceAdami(dest, sources)

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

initialize(d_au, d_av, d_idx)

loop(d_au, d_av, d_aw, d_idx, d_m, DWIJ, d_pi00, d_pi01, d_pi02, d_pi10, d_pi11, d_pi12, d_pi20, d_pi21, d_pi22, s_pi00, s_pi01, s_pi02, s_pi10, s_pi11, s_pi12, s_pi20, s_pi21, s_pi22, d_V, s_V, s_idx)

pysph.sph.surface_tension.get_surface_tension_equations(fluids, solids, scheme, rho0, p0, c0, b, factor1, factor2, nu, sigma, d, epsilon, gamma, real=False)

This function returns the required equations for the multiphase formulation taking inputs of the fluid particles array, solid particles array, the scheme to be used and other physical parameters

Parameters:

• fluids (list) – List of names of fluid particle arrays

• solids (list) – List of names of solid particle arrays

• scheme (string) –

  The scheme with which the equations are to be setup. Supported Schemes:

  1. TVF scheme with Morris’ surface tension. String to be used: “tvf”

  2. Adami’s surface tension implementation which doesn’t involve calculation of curvature. String to be used: “adami_stress”

  3. Adami’s surface tension implementation which involves calculation of curvature. String to be used: “adami”

  4. Shadloo Yildiz surface tension formulation. String to be used: “shadloo”

  5. Morris’ surface tension formulation. This is the default scheme which will be used if none of the above strings are input as scheme.

• rho0 (float) – The reference density of the medium (Currently multiple reference densities for different particles is not supported)

• p0 (float) – The background pressure of the medium(Currently multiple background pressures for different particles is not supported)

• c0 (float) – The speed of sound of the medium(Currently multiple speeds of sounds for different particles is not supported)

• b (float) – The b parameter of the generalized Tait Equation of State. Refer to the Tait Equation’s documentation for reference

• factor1 (float) – The factor for scaling of smoothing length for calculation of interface curvature number for shadloo’s scheme

• factor2 (float) – The factor for scaling back of smoothing length for calculation of forces after calculating the interface curvature number in shadloo’s scheme

• nu (float) – The kinematic viscosity of the medium

• sigma (float) – The surface tension of the system

• d (int) – The number of dimensions of the problem in the cartesian space

• epsilon (float) – Put this option false if the equations are supposed to be evaluated for the ghost particles, else keep it True

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)

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, d_uadv, d_vadv, d_wadv)

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)

class pysph.sph.iisph.ComputeAII(dest, sources)

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)

loop(d_idx, d_aii, d_dii0, d_dii1, d_dii2, d_m, d_rho, s_idx, s_m, s_rho, DWIJ)

class pysph.sph.iisph.ComputeAIIBoundary(dest, sources, rho0)

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

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_m, d_aii, d_dii0, d_dii1, d_dii2, d_rho, s_idx, s_m, s_V, DWIJ)

class pysph.sph.iisph.ComputeDII(dest, sources)

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)

loop(d_idx, d_rho, d_dii0, d_dii1, d_dii2, s_idx, s_m, DWIJ)

class pysph.sph.iisph.ComputeDIIBoundary(dest, sources, rho0)

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_dii0, d_dii1, d_dii2, d_rho, s_idx, s_m, s_V, DWIJ)

class pysph.sph.iisph.ComputeDIJPJ(dest, sources)

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)

loop(d_idx, d_dijpj0, d_dijpj1, d_dijpj2, s_idx, s_m, s_rho, s_piter, DWIJ)

class pysph.sph.iisph.ComputeRhoAdvection(dest, sources)

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)

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)

class pysph.sph.iisph.ComputeRhoBoundary(dest, sources, rho0)

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_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)

class pysph.sph.iisph.IISPHScheme(fluids, solids, dim, rho0, nu=0.0, gx=0.0, gy=0.0, gz=0.0, omega=0.5, tolerance=0.01, debug=False, has_ghosts=False)

The IISPH scheme

Parameters:

• fluids (list(str)) – List of names of fluid particle arrays.

• solids (list(str)) – List of names of solid particle arrays.

• dim (int) – Dimensionality of the problem.

• rho0 (float) – Density of fluid.

• nu (float) – Kinematic viscosity.

• gx (float) – Componenents of body acceleration (gravity, external forcing etc.)

• gy (float) – Componenents of body acceleration (gravity, external forcing etc.)

• gz (float) – Componenents of body acceleration (gravity, external forcing etc.)

• omega (float) – Relaxation parameter for relaxed-Jacobi iterations.

• tolerance (float) – Tolerance for the convergence of pressure iterations as a fraction.

• debug (bool) – Produce some debugging output on iterations.

• has_ghosts (bool) – The problem has ghost particles so add equations for those.

add_user_options(group)

configure_solver(kernel=None, integrator_cls=None, extra_steppers=None, **kw)

Configure the solver to be generated.

This is to be called before get_solver is called.

Parameters:

• dim (int) – Number of dimensions.

• kernel (Kernel instance.) – Kernel to use, if none is passed a default one is used.

• integrator_cls (pysph.sph.integrator.Integrator) – Integrator class to use, use sensible default if none is passed.

• extra_steppers (dict) – Additional integration stepper instances as a dict.

• **kw (extra arguments) – Any additional keyword args are passed to the solver instance.

consume_user_options(options)

get_equations()

setup_properties(particles, clean=True)

Setup the particle arrays so they have the right set of properties for this scheme.

Parameters:

• particles (list) – List of particle arrays.

• clean (bool) – If True, removes any unnecessary properties.

class pysph.sph.iisph.IISPHStep

A straightforward and simple integrator to be used for IISPH.

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)

class pysph.sph.iisph.NormalizedSummationDensity(dest, sources)

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)

loop(d_idx, d_rho, d_rho_adv, d_V, s_idx, s_m, s_rho0, WIJ)

post_loop(d_idx, d_rho, d_rho_adv)

class pysph.sph.iisph.NumberDensity(dest, sources)

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)

loop(d_idx, d_V, WIJ)

class pysph.sph.iisph.PressureForce(dest, sources)

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)

loop(d_idx, d_rho, d_p, d_au, d_av, d_aw, s_idx, s_m, s_rho, s_p, DWIJ)

post_loop(d_idx, d_au, d_av, d_aw, d_uadv, d_vadv, d_wadv, d_dt_cfl, d_dt_force)

class pysph.sph.iisph.PressureForceBoundary(dest, sources, rho0)

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_rho, d_au, d_av, d_aw, d_p, s_idx, s_V, DWIJ)

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

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

converged()

Return > 0 to indicate converged iterations and < 0 otherwise.

initialize(d_idx, d_p, d_compression)

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)

post_loop(d_idx, d_piter, d_p0, d_p, d_aii, d_rho_adv, d_rho, d_compression, dt=0.0)

reduce(dst, t, dt)

class pysph.sph.iisph.PressureSolveBoundary(dest, sources, rho0)

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_p, d_rho, d_dijpj0, d_dijpj1, d_dijpj2, s_idx, s_V, DWIJ)

class pysph.sph.iisph.SummationDensity(dest, sources)

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)

loop(d_idx, d_rho, s_idx, s_m, WIJ)

class pysph.sph.iisph.SummationDensityBoundary(dest, sources, rho0)

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_rho, s_idx, s_V, WIJ)

class pysph.sph.iisph.UpdateGhostPressure(dest, sources)

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_tag, d_orig_idx, d_p, d_piter)

class pysph.sph.iisph.UpdateGhostProps(dest, sources=None)

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_tag, d_orig_idx, d_dijpj0, d_dijpj1, d_dijpj2, d_dii0, d_dii1, d_dii2, d_piter)

class pysph.sph.iisph.ViscosityAcceleration(dest, sources, nu)

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_au, d_av, d_aw, s_idx, s_m, EPS, VIJ, XIJ, RHOIJ1, R2IJ, DWIJ)

class pysph.sph.iisph.ViscosityAccelerationBoundary(dest, sources, rho0, nu)

The acceleration on the fluid due to a boundary.

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_au, d_av, d_aw, d_rho, s_idx, s_V, EPS, VIJ, XIJ, R2IJ, DWIJ)

Hopkins’ ‘Traditional’ SPH (TSPH)

References

[Hopkins2013] (1,2)

Hopkins, Philip F. “A General Class of Lagrangian Smoothed Particle Hydrodynamics Methods and Implications for Fluid Mixing Problems.” Monthly Notices of the Royal Astronomical Society 428, no. 4 (February 1, 2013): 2840–56. https://doi.org/10.1093/mnras/sts210.

[Hopkins2015] (1,2,3,4,5,6,7,8,9,10)

Hopkins, Philip F. “A New Class of Accurate, Mesh-Free Hydrodynamic Simulation Methods.” Monthly Notices of the Royal Astronomical Society 450, no. 1 (June 11, 2015): 53–110. https://doi.org/10.1093/mnras/stv195.

class pysph.sph.gas_dynamics.tsph.TSPHScheme(fluids, solids, dim, gamma, hfact, beta=2.0, fkern=1.0, max_density_iterations=250, alphamax=1.0, density_iteration_tolerance=0.001, has_ghosts=False)

Density-energy formulation [Hopkins2013] including Balsara’s artificial viscocity switch with modifications, as presented in Appendix F1 of [Hopkins2015] .

Notes

Is this exactly in accordance with what is proposed in [Hopkins2015] ?

Not quite.

What is different then?

1. Adapting smoothing length using MPM [KP14] procedure from SummationDensity. In this, calculation of grad-h terms are changed to that specified for this scheme.

2. Using the PEC integrator step. No individual adaptive time-stepping.

3. Using Gaussian Kernel by default instead of Cubic Spline with radius scale 1.

Tip: Reduce the number of points if particle penetration is encountered. This has to be done while running gas_dynamics.wc_blastwave and gas_dynamics.robert.

Parameters:

• fluids (list) – List of names of fluid particle arrays.

• solids (list) – List of names of solid particle arrays (or boundaries), currently not supported

• dim (int) – Dimensionality of the problem.

• gamma (float) – \(\gamma\) for Equation of state.

• hfact (float) – \(h_{fact}\) for smoothing length adaptivity, also referred to as kernel_factor in other gas dynamics schemes.

• beta (float, optional) – \(\beta\) for artificial viscosity, by default 2.0

• fkern (float, optional) – \(f_{kern}\), Factor to scale smoothing length for equivalence with classic kernel when using kernel with altered radius_scale is being used, by default 1.

• max_density_iterations (int, optional) – Maximum number of iterations to run for one density step, by default 250.

• density_iteration_tolerance (float, optional) – Maximum difference allowed in two successive density iterations, by default 1e-3

• has_ghosts (bool, optional) – if ghost particles (either mirror or periodic) is used, by default False

• alphamax (float, optional) – \(\alpha_{av}\) for artificial viscosity switch, by default 1.0

add_user_options(group)

consume_user_options(options)

configure_solver(kernel=None, integrator_cls=None, extra_steppers=None, **kw)

Configure the solver to be generated.

Parameters:

• kernel (Kernel instance.) – Kernel to use, if none is passed a default one is used.

• integrator_cls (pysph.sph.integrator.Integrator) – Integrator class to use, use sensible default if none is passed.

• extra_steppers (dict) – Additional integration stepper instances as a dict.

• **kw (extra arguments) – Any additional keyword args are passed to the solver instance.

get_equations()

setup_properties(particles, clean=True)

Setup the particle arrays so they have the right set of properties for this scheme.

Parameters:

• particles (list) – List of particle arrays.

• clean (bool) – If True, removes any unnecessary properties.

class pysph.sph.gas_dynamics.tsph.SummationDensity(dest, sources, dim, density_iterations=False, iterate_only_once=False, hfact=1.2, htol=1e-06)

SummationDensity modified to use number density for calculation of grad-h terms.

Ref. Appendix F1 [Hopkins2015]

initialize(d_idx, d_rho, d_arho, d_drhosumdh, d_n, d_dndh, d_prevn, d_prevdndh, d_prevdrhosumdh, d_an)

loop(d_idx, s_idx, d_rho, d_arho, d_drhosumdh, s_m, VIJ, WI, DWI, GHI, d_n, d_dndh, d_h, d_prevn, d_prevdndh, d_prevdrhosumdh, d_an)

post_loop(d_idx, d_h0, d_h, d_ah, d_converged, d_n, d_dndh, d_an)

converged()

Return > 0 to indicate converged iterations and < 0 otherwise.

class pysph.sph.gas_dynamics.tsph.IdealGasEOS(dest, sources, gamma)

IdealGasEOS modified to avoid repeated calculations using loop(). Doing the same using post_loop().

post_loop(d_idx, d_p, d_rho, d_e, d_cs)

class pysph.sph.gas_dynamics.tsph.VelocityGradDivC1(dest, sources, dim)

First Order consistent velocity gradient and divergence

initialize(d_gradv, d_idx, d_invtt, d_divv)

loop(d_idx, d_invtt, s_m, s_idx, VIJ, DWI, XIJ, d_gradv)

post_loop(d_idx, d_gradv, d_invtt, d_divv)

class pysph.sph.gas_dynamics.tsph.BalsaraSwitch(dest, sources, alphaav, fkern)

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

post_loop(d_h, d_idx, d_cs, d_divv, d_gradv, d_alpha)

class pysph.sph.gas_dynamics.tsph.MomentumAndEnergy(dest, sources, dim, fkern, beta=2.0)

TSPH Momentum and Energy Equations with artificial viscosity.

Possible typo in that has been taken care of:

Instead of Equation F3 [Hopkins2015] for evolution of total energy sans artificial viscosity and artificial conductivity,

\[\frac{\mathrm{d} E_{i}}{\mathrm{~d} t}=\boldsymbol{v}_{i} \cdot \frac{\mathrm{d} \boldsymbol{P}_{i}}{\mathrm{~d} t}- \sum_{j} m_{i} m_{j}\left(\boldsymbol{v}_{i}- \boldsymbol{v}_{j}\right) \cdot\left[\frac{P_{i}} {\bar{\rho}_{i}^{2}} f_{i, j} \nabla_{i} W_{i j}\left(h_{i}\right)\right],\]

it should have been,

\[\frac{\mathrm{d} E_{i}}{\mathrm{~d} t}=\boldsymbol{v}_{i} \cdot \frac{\mathrm{d} \boldsymbol{P}_{i}}{\mathrm{~d} t}+ \sum_{j} m_{i} m_{j}\left(\boldsymbol{v}_{i}- \boldsymbol{v}_{j}\right) \cdot\left[\frac{P_{i}} {\bar{\rho}_{i}^{2}} f_{i, j} \nabla_{i} W_{i j}\left(h_{i}\right)\right].\]

Specific thermal energy, \(u\), would therefore be evolved using,

\[\frac{\mathrm{d} u_{i}}{\mathrm{~d} t}= \sum_{j} m_{j}\left(\boldsymbol{v}_{i}- \boldsymbol{v}_{j}\right) \cdot\left[\frac{P_{i}} {\bar{\rho}_{i}^{2}} f_{i, j} \nabla_{i} W_{i j}\left(h_{i}\right)\right]\]

initialize(d_idx, d_au, d_av, d_aw, d_ae)

loop(d_idx, s_idx, d_m, s_m, d_p, s_p, d_cs, s_cs, d_rho, s_rho, d_au, d_av, d_aw, d_ae, XIJ, VIJ, DWI, DWJ, HIJ, d_alpha, s_alpha, R2IJ, RHOIJ1, d_h, d_dndh, d_n, d_drhosumdh, s_h, s_dndh, s_n, s_drhosumdh)

class pysph.sph.gas_dynamics.tsph.WallBoundary(dest, sources)

WallBoundary modified for TSPH.

Most importantly, mass of the boundary particle should never be zero since it appears in denominator of fij. This has been addressed.

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_p, d_rho, d_e, d_m, d_cs, d_h, d_htmp, d_h0, d_u, d_v, d_w, d_wij, d_n, d_dndh, d_drhosumdh, d_divv, d_m0)

loop(d_idx, s_idx, d_p, d_rho, d_e, d_m, d_cs, d_divv, d_h, d_u, d_v, d_w, d_wij, d_htmp, s_p, s_rho, s_e, s_m, s_cs, s_h, s_divv, s_u, s_v, s_w, WI, s_n, d_n, s_dndh, d_dndh, d_drhosumdh, s_drhosumdh)

post_loop(d_idx, d_p, d_rho, d_e, d_m, d_cs, d_divv, d_h, d_u, d_v, d_w, d_wij, d_htmp, d_n, d_dndh, d_drhosumdh, d_m0)

class pysph.sph.gas_dynamics.tsph.UpdateGhostProps(dest, sources=None, dim=2)

MPMUpdateGhostProps modified for TSPH

initialize(d_idx, d_orig_idx, d_p, d_tag, d_h, d_rho, d_dndh, d_psumdh, d_n)

Hopkins’ ‘Modern’ SPH (PSPH)

References

[CullenDehnen2010] (1,2)

Cullen, Lee, and Walter Dehnen. “Inviscid Smoothed Particle Hydrodynamics: Inviscid Smoothed Particle Hydrodynamics.” Monthly Notices of the Royal Astronomical Society 408, no. 2 (October 21, 2010): 669–83. https://doi.org/10.1111/j.1365-2966.2010.17158.x.

[ReadHayfield2012]

Read, J. I., and T. Hayfield. “SPHS: Smoothed Particle Hydrodynamics with a Higher Order Dissipation Switch: SPH with a Higher Order Dissipation Switch.” Monthly Notices of the Royal Astronomical Society 422, no. 4 (June 1, 2012): 3037–55. https://doi.org/10.1111/j.1365-2966.2012.20819.x.

class pysph.sph.gas_dynamics.psph.PSPHScheme(fluids, solids, dim, gamma, hfact, betab=2.0, fkern=1.0, max_density_iterations=250, alphac=0.25, density_iteration_tolerance=0.001, has_ghosts=False, alphamin=0.02, alphamax=2.0, betac=0.7, betad=0.05, betaxi=1.0)

Pressure-energy formulation [Hopkins2013] including Cullen-Dehnen artificial viscocity switch [CullenDehnen2010] with modifications, as presented in Appendix F2 of [Hopkins2015] .

Notes

Is this exactly in accordance with what is proposed in [Hopkins2015] ?

Not quite.

What is different then?

1. Adapting smoothing length using MPM [KP14] procedure from SummationDensity. In this, calculation of grad-h terms are changed to that specified for this scheme.

2. Using the PEC integrator step. No individual adaptive time-stepping.

3. Using Gaussian Kernel by default instead of Cubic Spline with radius scale 1.

Tip: Reduce the number of points if particle penetration is encountered. This has to be done while running gas_dynamics.wc_blastwave and gas_dynamics.robert.

Parameters:

• fluids (list) – List of names of fluid particle arrays.

• solids (list) – List of names of solid particle arrays (or boundaries), currently not supported

• dim (int) – Dimensionality of the problem.

• gamma (float) – \(\gamma\) for Equation of state.

• hfact (float) – \(h_{fact}\) for smoothing length adaptivity, also referred to as kernel_factor in other gas dynamics schemes.

• betab (float, optional) – \(\beta_b\) for artificial viscosity, by default 2.0

• fkern (float, optional) – \(f_{kern}\), Factor to scale smoothing length for equivalence with classic kernel when using kernel with altered radius_scale is being used, by default 1.

• max_density_iterations (int, optional) – Maximum number of iterations to run for one density step, by default 250.

• alphac (float, optional) – \(\alpha_c\) for artificial conductivity, by default 0.25

• density_iteration_tolerance (float, optional) – Maximum difference allowed in two successive density iterations, by default 1e-3

• has_ghosts (bool, optional) – if ghost particles (either mirror or periodic) is used, by default False

• alphamin (float, optional) – \(\alpha_{min}\) for artificial viscosity switch, by default 0.02

• alphamax (float, optional) – \(\alpha_{max}\) for artificial viscosity switch, by default 2.0

• betac (float, optional) – \(\beta_c\) for artificial viscosity switch, by default 0.7

• betad (float, optional) – \(\beta_d\) for artificial viscosity switch, by default 0.05

• betaxi (float, optional) – \(\beta_{\xi}\) for artificial viscosity switch, by default 1.0

add_user_options(group)

consume_user_options(options)

configure_solver(kernel=None, integrator_cls=None, extra_steppers=None, **kw)

Configure the solver to be generated.

Parameters:

• kernel (Kernel instance.) – Kernel to use, if none is passed a default one is used.

• integrator_cls (pysph.sph.integrator.Integrator) – Integrator class to use, use sensible default if none is passed.

• extra_steppers (dict) – Additional integration stepper instances as a dict.

• **kw (extra arguments) – Any additional keyword args are passed to the solver instance.

get_equations()

setup_properties(particles, clean=True)

Setup the particle arrays so they have the right set of properties for this scheme.

Parameters:

• particles (list) – List of particle arrays.

• clean (bool) – If True, removes any unnecessary properties.

class pysph.sph.gas_dynamics.psph.PSPHSummationDensityAndPressure(dest, sources, dim, gamma, density_iterations=False, iterate_only_once=False, hfact=1.2, htol=1e-06)

SummationDensity modified to use number density for calculation of grad-h terms and to calculate pressure and speed of sound as well.

Ref. Appendix F2 [Hopkins2015]

Parameters:

• density_iterations (bint, optional) – Flag to indicate density iterations are required, by default False

• iterate_only_once (bint, optional) –

  Flag to indicate if only one iteration is required,

  by default False

• hfact (float, optional) – \(h_{fact}\), by default 1.2

• htol (double, optional) – Iteration tolerance, by default 1e-6

initialize(d_idx, d_rho, d_arho, d_n, d_dndh, d_prevn, d_prevdndh, d_p, d_dpsumdh, d_dprevpsumdh, d_an)

loop(d_idx, s_idx, d_rho, d_arho, s_m, VIJ, WI, DWI, GHI, d_n, d_dndh, d_h, d_prevn, d_prevdndh, s_e, d_p, d_dpsumdh, d_e, d_an)

post_loop(d_idx, d_rho, d_h0, d_h, d_ah, d_converged, d_cs, d_p, d_n, d_dndh, d_an)

converged()

Return > 0 to indicate converged iterations and < 0 otherwise.

class pysph.sph.gas_dynamics.psph.GradientKinsfolkC1(dest, sources, dim)

First order consistent,

• Velocity gradient, grad(v)

• Acceleration gradient, grad(a)

• Velocity divergence, div(v)

• Velocity divergence rate, d (div(v)) / dt

• Traceless symmetric strain rate, S

• trace(dot(S,transpose(S)))

Ref. Appendix B [CullenDehnen2010]

initialize(d_gradv, d_idx, d_invtt, d_divv, d_grada, d_adivv, d_trssdsst)

loop(d_idx, d_invtt, s_m, s_idx, VIJ, DWI, XIJ, d_gradv, d_grada, d_au, s_au, d_av, s_av, d_aw, s_aw)

post_loop(d_idx, d_gradv, d_invtt, d_divv, d_grada, d_adivv, d_ss, d_trssdsst)

class pysph.sph.gas_dynamics.psph.LimiterAndAlphas(dest, sources, alphamin=0.02, alphamax=2.0, betac=0.7, betad=0.05, betaxi=1.0, fkern=1.0)

Cullen Dehnen’s limiter for artificial viscosity modified by Hopkins.

Ref. Appendix F2 [Hopkins2015]

initialize(d_idx, d_xi)

loop(d_idx, s_idx, s_m, d_xi, s_divv, WI)

post_loop(d_idx, d_xi, d_rho, d_h, d_adivv, d_cs, d_alpha0, d_vsig, dt, d_divv, d_trssdsst, d_alpha)

class pysph.sph.gas_dynamics.psph.MomentumAndEnergy(dest, sources, dim, fkern, gamma, betab=2.0, alphac=0.25)

PSPH Momentum and Energy Equations with artificial viscosity and artificial conductivity.

Possible typos in that have been taken care of,

1. Instead of Equation F15 [Hopkins2015] for evolution of total energy sans artificial viscosity and artificial conductivity,

   \[\frac{\mathrm{d} E_{i}}{\mathrm{~d} t}= \boldsymbol{v}_{i} \cdot \frac{\mathrm{d} \boldsymbol{P}_{i}}{\mathrm{~d} t}- \sum_{j=1}^{N}(\gamma-1)^{2} m_{i} m_{j} u_{i} u_{j} \frac{f_{i j}}{\bar{P}_{i}}\left(\boldsymbol{v}_{i}- \boldsymbol{v}_{j}\right) \cdot \nabla_{i} W_{i j} \left(h_{i}\right),\]

   it should have been,

   \[\frac{\mathrm{d} E_{i}}{\mathrm{~d} t}= \boldsymbol{v}_{i} \cdot \frac{\mathrm{d} \boldsymbol{P}_{i}}{\mathrm{~d} t}+ \sum_{j=1}^{N}(\gamma-1)^{2} m_{i} m_{j} u_{i} u_{j} \frac{f_{i j}}{\bar{P}_{i}}\left(\boldsymbol{v}_{i}- \boldsymbol{v}_{j}\right) \cdot \nabla_{i} W_{i j} \left(h_{i}\right).\]

   Specific thermal energy, \(u\), would therefore be evolved using,

   \[\frac{\mathrm{d} u_{i}}{\mathrm{~d} t}= \sum_{j=1}^{N}(\gamma-1)^{2} m_{j} u_{i} u_{j} \frac{f_{i j}}{\bar{P}_{i}}\left(\boldsymbol{v}_{i}- \boldsymbol{v}_{j}\right) \cdot \nabla_{i} W_{i j} \left(h_{i}\right).\]

2. Equation F18 [Hopkins2015] for contribution of artificial viscosity to the evolution of total energy is,

   \[\frac{\mathrm{d} E_{i}}{\mathrm{~d} t}= \alpha_{\mathrm{C}} \sum_{j} m_{i} m_{j} \alpha_{i j} \tilde{v}_{s}\left(u_{i}- u_{j}\right) \times \frac{\left|P_{i}-P_{j}\right|}{P_{i}+ P_{j}} \frac{\nabla_{i} W_{i j}\left(h_{i}\right)+ \nabla_{i} W_{i j}\left(h_{j}\right)}{\bar{\rho}_{i}+ \bar{\rho}_{j}} .\]

   Carefully comparing with [ReadHayfield2012] and [KP14], specific thermal energy, \(u\), should be evolved using,

   \[\begin{split}\frac{\mathrm{d} u_{i}}{\mathrm{~d} t}= \alpha_{\mathrm{C}} \sum_{j} & m_{j} \alpha_{i j} \tilde{v}_{s}\left(u_{i}- u_{j}\right) \frac{\left|P_{i}-P_{j}\right|}{P_{i}+ P_{j}} \\ & \frac{\nabla_{i} W_{i j}\left(h_{i}\right)+ \nabla_{i} W_{i j}\left(h_{j}\right)}{\bar{\rho}_{i}+ \bar{\rho}_{j}} \cdot \frac{\left(\boldsymbol{x}_{i}- \boldsymbol{x}_{j}\right)}{\left|\boldsymbol{x}_{i}- \boldsymbol{x}_{j}\right|}\end{split}\]

initialize(d_idx, d_au, d_av, d_aw, d_ae)

loop(d_idx, s_idx, d_m, s_m, d_p, s_p, d_cs, s_cs, d_au, d_av, d_aw, d_ae, XIJ, VIJ, DWI, DWJ, d_alpha, s_alpha, RIJ, d_h, d_dndh, d_n, s_h, s_dndh, s_n, d_e, s_e, d_dpsumdh, s_dpsumdh, RHOIJ1)

class pysph.sph.gas_dynamics.psph.WallBoundary(dest, sources)

WallBoundary modified for PSPH

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_p, d_rho, d_e, d_m, d_cs, d_h, d_htmp, d_h0, d_u, d_v, d_w, d_wij, d_n, d_dndh, d_dpsumdh, d_m0)

loop(d_idx, s_idx, d_p, d_rho, d_e, d_m, d_cs, d_u, d_v, d_w, d_wij, d_htmp, s_p, s_rho, s_e, s_m, s_cs, s_h, s_u, s_v, s_w, WI, s_n, d_n, d_dndh, s_dndh, d_dpsumdh, s_dpsumdh)

post_loop(d_idx, d_p, d_rho, d_e, d_m, d_cs, d_h, d_u, d_v, d_w, d_wij, d_htmp, d_dndh, d_dpsumdh, d_n, d_m0)

class pysph.sph.gas_dynamics.psph.UpdateGhostProps(dest, sources=None, dim=2)

MPMUpdateGhostProps modified for PSPH

initialize(d_idx, d_orig_idx, d_p, d_tag, d_h, d_rho, d_dndh, d_psumdh, d_n)

MAGMA2

References

[Rosswog2009]

Rosswog, Stephan. “Astrophysical smooth particle hydrodynamics.” New Astronomy Reviews 53, no. 4-6 (2009): 78-104. https://doi.org/10.1016/j.newar.2009.08.007

[Rosswog2015]

Rosswog, Stephan. “Boosting the accuracy of SPH techniques: Newtonian and special-relativistic tests.” Monthly Notices of the Royal Astronomical Society 448, no. 4 (2015): 3628-3664. https://doi.org/10.1093/mnras/stv225.

[Rosswog2020a]

Rosswog, Stephan. “A simple, entropy-based dissipation trigger for SPH.” The Astrophysical Journal 898, no. 1 (2020): 60. https://doi.org/10.3847/1538-4357/ab9a2e.

[Rosswog2020b] (1,2,3,4)

Rosswog, Stephan. “The Lagrangian hydrodynamics code MAGMA2.” Monthly Notices of the Royal Astronomical Society 498, no. 3 (2020): 4230-4255. https://doi.org/10.1093/mnras/staa2591.

class pysph.sph.gas_dynamics.magma2.IncreaseSmoothingLength(dest, sources)

Increase smoothing length by 10%.

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_h)

class pysph.sph.gas_dynamics.magma2.UpdateSmoothingLength(dest, sources, ndes)

Sorts neighbours based on distance and uses the distance of nearest \((n_{des} + 1)^{th}\) particle to set the smoothing length. Here, \(n_{des}\) is the desired number of neighbours to be in the kernel support of each particle.

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

loop_all(d_idx, d_x, d_y, d_z, d_h, s_x, s_y, s_z, NBRS, N_NBRS, SPH_KERNEL)

class pysph.sph.gas_dynamics.magma2.SummationDensityMPMStyle(dest, sources, dim, density_iterations=False, iterate_only_once=False, hfact=1.2, htol=1e-06)

SummationDensity modified to use number density and without grad-h terms.

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_arho, d_n, d_dndh, d_prevn, d_prevdndh, d_an)

loop(d_idx, s_idx, d_rho, d_arho, s_m, VIJ, WI, DWI, GHI, d_n, d_dndh, d_an)

post_loop(d_idx, d_h0, d_h, d_ah, d_converged, d_n, d_dndh, d_an)

converged()

Return > 0 to indicate converged iterations and < 0 otherwise.

class pysph.sph.gas_dynamics.magma2.IdealGasEOS(dest, sources, gamma)

IdealGasEOS modified to avoid repeated calculations using loop(). Doing the same using post_loop().

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

post_loop(d_idx, d_p, d_rho, d_e, d_cs)

class pysph.sph.gas_dynamics.magma2.AuxiliaryGradient(dest, sources, dim)

Auxiliary first gradient calculated using analytical gradient of kernel and without using density.

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

initialize(d_dvaux, d_idx, d_invdm, d_deaux)

loop(d_idx, VIJ, XIJ, d_invdm, DWI, d_dvaux, s_m, s_idx, d_deaux, d_e, s_e)

post_loop(d_idx, d_invdm, d_dvaux, d_deaux)

class pysph.sph.gas_dynamics.magma2.CorrectionMatrix(dest, sources, dim)

Correction matrix, C, that accounts for the local particle distribution and used in calculation of gradients without using analytical derivatives of kernel.

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

initialize(d_cm, d_idx)

loop(d_idx, s_m, s_idx, XIJ, s_rho, d_cm, WI)

post_loop(d_idx, d_cm)

class pysph.sph.gas_dynamics.magma2.FirstGradient(dest, sources, dim)

First gradient and divergence calculated using matrix inversion formulation without analytical derivative of the kernel.

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

initialize(d_dv, d_idx, d_divv, d_de)

loop(d_idx, VIJ, XIJ, d_dv, WI, s_m, s_rho, s_idx, d_e, s_e, d_de)

post_loop(d_idx, d_dv, d_divv, d_cm, d_de)

class pysph.sph.gas_dynamics.magma2.SecondGradient(dest, sources, dim)

Second gradient calculated from auxiliary gradient using matrix inversion formulation without analytical derivative of the kernel.

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

initialize(d_ddv, d_idx, d_dde)

loop(d_idx, XIJ, d_dvaux, s_dvaux, WI, d_ddv, s_m, s_rho, s_idx, s_deaux, d_deaux, d_dde)

post_loop(d_idx, d_cm, d_ddv, d_dde)

class pysph.sph.gas_dynamics.magma2.EntropyBasedDissipationTrigger(dest, sources, alphamax, alphamin, fkern, l0, l1, gamma)

Simple, entropy-based dissipation trigger from [Rosswog2020a]

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

post_loop(d_h, d_idx, d_cs, d_alpha, d_s, d_p, d_rho, dt, d_aalpha)

class pysph.sph.gas_dynamics.magma2.WallBoundary(dest, sources, dim)

WallBoundary modified for MAGMA2.

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_p, d_rho, d_e, d_m, d_cs, d_h, d_htmp, d_h0, d_u, d_v, d_w, d_wij, d_n, d_dndh, d_divv, d_alpha, d_ddv, d_dv, d_de, d_cm, d_dde, d_rho0)

loop(d_idx, s_idx, d_p, d_rho, d_e, d_m, d_cs, d_divv, d_u, d_v, d_w, d_wij, d_htmp, s_p, s_rho, s_e, s_m, s_cs, s_h, s_divv, s_u, s_v, s_w, WI, s_n, d_n, s_dndh, d_dndh, d_alpha, s_alpha, d_de, s_de, d_dv, d_cm, d_dde, s_dv, s_cm, s_dde, s_ddv, d_ddv)

post_loop(d_idx, d_p, d_rho, d_e, d_m, d_cs, d_divv, d_h, d_u, d_v, d_w, d_wij, d_htmp, d_n, d_dndh, d_de, d_dv, d_cm, d_dde, d_ddv, d_rho0)

class pysph.sph.gas_dynamics.magma2.MomentumAndEnergyStdGrad(dest, sources, dim, fkern, eta_crit=0.3, eta_fold=0.2, beta=2.0, alphac=0.05, eps=0.01)

Standard Gradient formulation (stdGrad) momentum and energy equations with artificial viscosity and artificial conductivity from [Rosswog2020b]

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

loop(d_idx, s_idx, s_m, d_p, s_p, d_cs, s_cs, d_rho, s_rho, d_au, d_av, d_aw, d_ae, XIJ, VIJ, d_alpha, s_alpha, d_ddv, s_ddv, RHOIJ1, d_h, s_h, DWI, DWJ, d_dv, s_dv, d_de, s_de, d_dde, s_dde, d_e, s_e)

class pysph.sph.gas_dynamics.magma2.MomentumAndEnergyMI1(dest, sources, dim, fkern, eta_crit=0.3, eta_fold=0.2, beta=2.0, alphac=0.05, eps=0.01)

Matrix inversion formulation 1 (MI1) momentum and energy equations with artificial viscosity and artificial conductivity from [Rosswog2020b]

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

loop(d_idx, s_idx, s_m, d_p, s_p, d_cs, s_cs, d_rho, s_rho, d_au, d_av, d_aw, d_ae, XIJ, VIJ, d_alpha, s_alpha, d_ddv, s_ddv, RHOIJ1, d_h, s_h, d_cm, s_cm, WI, WJ, d_dv, s_dv, d_de, s_de, d_dde, s_dde, d_e, s_e)

class pysph.sph.gas_dynamics.magma2.MomentumAndEnergyMI2(dest, sources, dim, fkern, eta_crit=0.3, eta_fold=0.2, beta=2.0, alphac=0.05, eps=0.01)

Matrix inversion formulation 2 (MI2) momentum and energy equations with artificial viscosity and artificial conductivity from [Rosswog2020b]

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

loop(d_idx, s_idx, s_m, d_p, s_p, d_cs, s_cs, d_rho, s_rho, d_au, d_av, d_aw, d_ae, XIJ, VIJ, d_alpha, s_alpha, d_ddv, s_ddv, RHOIJ1, d_h, s_h, d_cm, s_cm, WI, WJ, d_dv, s_dv, d_de, s_de, d_dde, s_dde, d_e, s_e)

class pysph.sph.gas_dynamics.magma2.EvaluateTildeMu(dest, sources, dim)

Find \(\tilde{\mu}\) to calculate time step.

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_tilmu)

loop(d_tilmu, d_idx, d_h, VIJ, XIJ, R2IJ)

class pysph.sph.gas_dynamics.magma2.SettleByArtificialPressure(dest, sources, xi=0.5, fkern=1.0)

Equation 40 of [Rosswog2020b] . Combined with an equation to update density (and smoothing length, if preferred), this equation can be evaluated through SPHEvaluator to settle the particles and obtain an initial distribution.

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

initialize(d_deltax, d_deltay, d_deltaz, d_idx, d_n, d_pouerr)

loop(d_rho, d_idx, d_rhodes, s_rho, s_rhodes, s_idx, d_deltax, d_deltay, d_deltaz, DWI, d_n, WI, s_m, d_pouerr)

post_loop(d_deltax, d_deltay, d_deltaz, d_idx, d_h, d_m, d_pouerr, d_rhodes, d_n, d_x, d_y, d_z)

Rigid body motion

Rigid body related equations.

class pysph.sph.rigid_body.AkinciRigidFluidCoupling(dest, sources, fluid_rho=1000)

Force between a solid sphere and a SPH fluid particle. This is implemented using Akinci’s[1] force and additional force from solid bodies pressure which is implemented by Liu[2]

[1]’Versatile Rigid-Fluid Coupling for Incompressible SPH’

URL: https://graphics.ethz.ch/~sobarbar/papers/Sol12/Sol12.pdf

[2]A 3D Simulation of a Moving Solid in Viscous Free-Surface Flows by Coupling SPH and DEM

https://doi.org/10.1155/2017/3174904

Note: Here forces for both the phases are added at once.

Please make sure that this force is applied only once for both the particle properties.

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_m, d_rho, d_au, d_av, d_aw, d_p, s_idx, s_V, s_fx, s_fy, s_fz, DWIJ, s_m, s_p, s_rho)

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

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_fx, d_fy, d_fz)

class pysph.sph.rigid_body.EulerStepRigidBody

Fast but inaccurate integrator. Use this for testing

initialize()

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)

class pysph.sph.rigid_body.LiuFluidForce(dest, sources)

Force between a solid sphere and a SPH fluid particle. This is implemented using Akinci’s[1] force and additional force from solid bodies pressure which is implemented by Liu[2]

[1]’Versatile Rigid-Fluid Coupling for Incompressible SPH’

URL: https://graphics.ethz.ch/~sobarbar/papers/Sol12/Sol12.pdf

[2]A 3D Simulation of a Moving Solid in Viscous Free-Surface Flows by Coupling SPH and DEM

https://doi.org/10.1155/2017/3174904

Note: Here forces for both the phases are added at once.

Please make sure that this force is applied only once for both the particle properties.

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_m, d_rho, d_au, d_av, d_aw, d_p, s_idx, s_V, s_fx, s_fy, s_fz, DWIJ, s_m, s_p, s_rho)

class pysph.sph.rigid_body.NumberDensity(dest, sources)

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)

loop(d_idx, d_V, WIJ)

class pysph.sph.rigid_body.PressureRigidBody(dest, sources, rho0)

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.

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_m, d_rho, d_au, d_av, d_aw, d_p, s_idx, s_V, s_fx, s_fy, s_fz, DWIJ)

class pysph.sph.rigid_body.RK2StepRigidBody

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)

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)

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)

class pysph.sph.rigid_body.RigidBodyCollision(dest, sources, kn=1000.0, mu=0.5, en=0.8)

Force between two spheres is implemented using DEM contact force law.

Refer https://doi.org/10.1016/j.powtec.2011.09.019 for more information.

Open-source MFIX-DEM software for gas–solids flows: Part I—Verification studies .

Initialise the required coefficients for force calculation.

Keyword arguments: kn – Normal spring stiffness (default 1e3) mu – friction coefficient (default 0.5) en – coefficient of restitution (0.8)

Given these coefficients, tangential spring stiffness, normal and tangential damping coefficient are calculated by default.

loop(d_idx, d_fx, d_fy, d_fz, d_h, d_total_mass, d_rad_s, d_tang_disp_x, d_tang_disp_y, d_tang_disp_z, d_tang_velocity_x, d_tang_velocity_y, d_tang_velocity_z, s_idx, s_rad_s, XIJ, RIJ, R2IJ, VIJ)

class pysph.sph.rigid_body.RigidBodyForceGPUGems(dest, sources, k=1.0, d=1.0, eta=1.0, kt=1.0)

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)

class pysph.sph.rigid_body.RigidBodyMoments(dest, sources)

Parameters:

• dest (str) – name of the destination particle array

• sources (list of str or None) – names of the source particle arrays

reduce(dst, t, dt)

class pysph.sph.rigid_body.RigidBodyMotion(dest, sources)

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)

class pysph.sph.rigid_body.RigidBodyWallCollision(dest, sources, kn=1000.0, mu=0.5, en=0.8)

Force between sphere and a wall is implemented using DEM contact force law.

Refer https://doi.org/10.1016/j.powtec.2011.09.019 for more information.

Open-source MFIX-DEM software for gas–solids flows: Part I—Verification studies .

Initialise the required coefficients for force calculation.

Keyword arguments: kn – Normal spring stiffness (default 1e3) mu – friction coefficient (default 0.5) en – coefficient of restitution (0.8)

Given these coefficients, tangential spring stiffness, normal and tangential damping coefficient are calculated by default.

loop(d_idx, d_fx, d_fy, d_fz, d_h, d_total_mass, d_rad_s, d_tang_disp_x, d_tang_disp_y, d_tang_disp_z, d_tang_velocity_x, d_tang_velocity_y, d_tang_velocity_z, s_idx, XIJ, RIJ, R2IJ, VIJ, s_nx, s_ny, s_nz)

class pysph.sph.rigid_body.SummationDensityBoundary(dest, sources, fluid_rho=1000.0)

Equation to find the density of the fluid particle due to any boundary or a rigid body

\(\rho_a = \sum_b {\rho}_fluid V_b W_{ab}\)

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_rho, s_idx, s_m, s_V, WIJ)

class pysph.sph.rigid_body.SummationDensityRigidBody(dest, sources, rho0)

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)

loop(d_idx, d_rho, s_idx, s_V, WIJ)

class pysph.sph.rigid_body.ViscosityRigidBody(dest, sources, rho0, nu)

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.

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_m, d_au, d_av, d_aw, d_rho, s_idx, s_V, s_fx, s_fy, s_fz, EPS, VIJ, XIJ, R2IJ, DWIJ)

pysph.sph.rigid_body.declare(*args)

pysph.sph.rigid_body.get_alpha_dot()

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()

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

pysph.sph.rigid_body.skew(vec)

Miscellaneous

Functions for advection

class pysph.sph.misc.advection.Advect(dest, sources)

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_ax, d_ay, d_u, d_v)

class pysph.sph.misc.advection.MixingVelocityUpdate(dest, sources, T)

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_u, d_v, d_u0, d_v0, t=0.1)

Functions to reduce array data in serial or parallel.

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

Simply returns the array for the serial case.

pysph.base.reduce_array.mpi_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.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')

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, pre=None, post=None, condition=None, start_idx=0, stop_idx=None, name=None)

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.

• pre (callable) – A callable which is passed no arguments that is called before anything in the group is executed.

• post (callable) – A callable which is passed no arguments that is called after the group is completed.

• condition (callable) – A callable that is passed (t, dt). If this callable returns True, the group is executed, otherwise it is not. If condition is None, the group is always executed. Note that this should work even if the group has many destination arrays.

• start_idx (int or str) – Start looping from this destination index. Starts from the given number if an integer is passed. If a string is look for a property/constant and use its first value as the loop count.

• stop_idx (int or str) – Loop up to this destination index instead of over all possible values. Defaults to all particles. Ends at the given number if an integer is passed. If a string is passed, look for a property/constant and use its first value as the loop count. Note that this works like a range stop parameter so the last value is not included.

• name (str) – The passed string is used to name the Group in the profiling info csv file to make it easy to read. If a string is not passed it defaults to the name ‘Group’.

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.

class pysph.sph.equation.MultiStageEquations(groups)

A class that allows a user to specify different equations for different stages.

The object doesn’t do much, except contain the different collections of equations.

Parameters:

groups (list/tuple) – A list/tuple of list of groups/equations, one for each stage.