# 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

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

  References

  Monaghan2005

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

  Parameters

    - alpha (float) – produces a shear and bulk viscosity

    - beta (float) – used to handle high Mach number shocks

  initialize(d_idx, d_au, d_av, d_aw)
  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}\]

  References

  Monaghan1992

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

  Parameters

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

  Notes

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

  initialize(d_idx, d_ax, d_ay, d_az)
  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}\]

  References

  Monaghan1992

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

  Parameters

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

  Notes

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

  initialize(d_idx, d_ax, d_ay, d_az)
  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

  - 16. Sun, A. Colagrossi, S. Marrone, A. Zhang

  “The 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)

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)

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

## 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

## Surface tension

## Implicit Incompressible SPH

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

  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.

    - pre – A callable which is passed no arguments that is called after the group is completed.

  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.