Skip to content

Commit

Permalink
Merge pull request #176 from abhinavmuta/deltasph
Browse files Browse the repository at this point in the history 
Fix delta sph equations
  • Loading branch information
@prabhuramachandran
prabhuramachandran committed Feb 27, 2019
2 parents 6a4dcb6 + 8e12803 commit 0b273ee
Show file tree
Hide file tree
Showing 3 changed files with 135 additions and 91 deletions.
95 changes: 58 additions & 37 deletions pysph/sph/scheme.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -392,10 +392,14 @@ def get_equations(self):
UpdateSmoothingLengthFerrari
)
from pysph.sph.wc.basic import (ContinuityEquationDeltaSPH,
ContinuityEquationDeltaSPHPreStep,
MomentumEquationDeltaSPH)
from pysph.sph.basic_equations import \
(ContinuityEquation, SummationDensity, XSPHCorrection)
from pysph.sph.wc.viscosity import LaminarViscosity
from pysph.sph.wc.viscosity import (LaminarViscosity,
LaminarViscosityDeltaSPH)
from pysph.sph.wc.kernel_correction import (GradientCorrectionPreStep,
GradientCorrection)

equations = []
g1 = []
Expand Down Expand Up @@ -425,58 +429,69 @@ def get_equations(self):
dest=name, sources=None, rho0=self.rho0, c0=self.c0,
gamma=self.gamma
))

equations.append(Group(equations=g1, real=False))

if self.delta_sph and not self.summation_density:
eq2_pre = []
for name in self.fluids:
eq2_pre.append(
GradientCorrectionPreStep(dest=name, sources=[name],
dim=self.dim)
)
equations.append(Group(equations=eq2_pre, real=False))

eq2 = []
for name in self.fluids:
eq2.extend([
GradientCorrection(dest=name, sources=[name]),
ContinuityEquationDeltaSPHPreStep(
dest=name, sources=[name]
)])
equations.append(Group(equations=eq2))

g2 = []
for name in self.solids:
g2.append(ContinuityEquation(dest=name, sources=self.fluids))

for name in self.fluids:
if self.delta_sph:
other = all[:]
other.remove(name)
if not self.summation_density:
g2.append(
ContinuityEquation(dest=name, sources=all)
)
if self.delta_sph and not self.summation_density:
g2.append(
ContinuityEquationDeltaSPH(
dest=name, sources=[name], c0=self.c0,
delta=self.delta
)
)
if len(other) > 0:
g2.append(ContinuityEquation(dest=name, sources=other))
))
# This is required since MomentumEquation (ME) adds artificial
# viscosity (AV), so make alpha 0.0 for ME and enable delta sph AV.
alpha = 0.0 if self.delta_sph else self.alpha
g2.append(
MomentumEquation(
dest=name, sources=all, c0=self.c0,
alpha=alpha, beta=self.beta,
gx=self.gx, gy=self.gy, gz=self.gz,
tensile_correction=self.tensile_correction
))
if self.delta_sph:
g2.append(
MomentumEquationDeltaSPH(
dest=name, sources=[name], rho0=self.rho0, c0=self.c0,
alpha=self.alpha, gx=self.gx, gy=self.gy, gz=self.gz,
)
)
if len(other) > 0:
g2.append(
MomentumEquation(
dest=name, sources=other, c0=self.c0,
alpha=self.alpha, beta=self.beta,
gx=self.gx, gy=self.gy, gz=self.gz,
tensile_correction=self.tensile_correction
)
)

g2.append(XSPHCorrection(dest=name, sources=[name]))
else:
if not self.summation_density:
g2.append(ContinuityEquation(dest=name, sources=all))
g2.extend([
MomentumEquation(
dest=name, sources=all, alpha=self.alpha,
beta=self.beta, gx=self.gx, gy=self.gy, gz=self.gz,
c0=self.c0, tensile_correction=self.tensile_correction
),
XSPHCorrection(dest=name, sources=[name])
])
alpha=self.alpha
))
g2.append(XSPHCorrection(dest=name, sources=[name]))

if abs(self.nu) > 1e-14:
eq = LaminarViscosity(
dest=name, sources=all, nu=self.nu
)
if self.delta_sph:
eq = LaminarViscosityDeltaSPH(
dest=name, sources=all, dim=self.dim, rho0=self.rho0,
nu=self.nu
)
else:
eq = LaminarViscosity(
dest=name, sources=all, nu=self.nu
)
g2.insert(-1, eq)
equations.append(Group(equations=g2))

Expand All @@ -496,6 +511,12 @@ def setup_properties(self, particles, clean=True):
props = list(dummy.properties.keys())
output_props = ['x', 'y', 'z', 'u', 'v', 'w', 'rho', 'm', 'h',
'pid', 'gid', 'tag', 'p']
if self.delta_sph:
delta_sph_props = [
{'name': 'm_mat', 'stride': 9},
{'name': 'gradrho', 'stride': 3},
]
props += delta_sph_props
for pa in particles:
self._ensure_properties(pa, props, clean)
pa.set_output_arrays(output_props)
Expand Down
84 changes: 37 additions & 47 deletions pysph/sph/wc/basic.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -296,8 +296,8 @@ class MomentumEquationDeltaSPH(Equation):
pp 1468--1480.
"""
def __init__(self, dest, sources, rho0, c0, alpha=1.0,
gx=0.0, gy=0.0, gz=0.0):

def __init__(self, dest, sources, rho0, c0, alpha=1.0):
r"""
Parameters
----------
Expand All @@ -308,12 +308,6 @@ def __init__(self, dest, sources, rho0, c0, alpha=1.0,
alpha : float
coefficient used to control the intensity of the
diffusion of velocity
gx : float
body force per unit mass along the x-axis
gy : float
body force per unit mass along the y-axis
gz : float
body force per unit mass along the z-axis
Notes
-----
Expand All @@ -324,52 +318,55 @@ def __init__(self, dest, sources, rho0, c0, alpha=1.0,
"""

self.alpha = alpha
self.gx = gx
self.gy = gy
self.gz = gz

self.c0 = c0
self.rho0 = rho0

super(MomentumEquationDeltaSPH, self).__init__(dest, sources)

def initialize(self, d_idx, d_au, d_av, d_aw):
d_au[d_idx] = 0.0
d_av[d_idx] = 0.0
d_aw[d_idx] = 0.0

def loop(self, 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):

# src paricle volume mj/rhoj
Vj = s_m[s_idx]/s_rho[s_idx]

p_i = d_p[d_idx]
pj = s_p[s_idx]

# viscous contribution second part of eqn (5b) in [Maronne2011]
vijdotxij = VIJ[0]*XIJ[0] + VIJ[1]*XIJ[1] + VIJ[2]*XIJ[2]
piij = self.alpha * HIJ * self.c0 * self.rho0 * vijdotxij/(R2IJ + EPS)
fac = self.alpha * HIJ * self.c0 * self.rho0
piij = vijdotxij/(R2IJ + EPS)

# gradient and viscous terms eqn 5b in [Maronne2011]
tmp = -Vj/d_rho[d_idx] * (p_i + pj) + piij * Vj/d_rho[d_idx]
tmp = fac * piij * Vj/d_rho[d_idx]

# accelerations
d_au[d_idx] += tmp * DWIJ[0]
d_av[d_idx] += tmp * DWIJ[1]
d_aw[d_idx] += tmp * DWIJ[2]

def post_loop(self, d_idx, d_au, d_av, d_aw, d_dt_force):
d_au[d_idx] += self.gx
d_av[d_idx] += self.gy
d_aw[d_idx] += self.gz

acc2 = (d_au[d_idx]*d_au[d_idx] +
d_av[d_idx]*d_av[d_idx] +
d_aw[d_idx]*d_aw[d_idx])
class ContinuityEquationDeltaSPHPreStep(Equation):
r"""**Continuity equation with dissipative terms**
See :class:`pysph.sph.wc.basic.ContinuityEquationDeltaSPH`
The matrix :math:`L_a` is multiplied to :math:`\nabla W_{ij}` in the
:class:`pysph.sph.scheme.WCSPHScheme` class by using
:class:`pysph.sph.wc.kernel_correction.GradientCorrectionPreStep` and
:class:`pysph.sph.wc.kernel_correction.GradientCorrection`.
"""

# store the square of the max acceleration
d_dt_force[d_idx] = acc2
def initialize(self, d_idx, d_gradrho):
d_gradrho[d_idx*3 + 0] = 0.0
d_gradrho[d_idx*3 + 1] = 0.0
d_gradrho[d_idx*3 + 2] = 0.0

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

rhoi = d_rho[d_idx]
rhoj = s_rho[s_idx]
Vj = s_m[s_idx]/rhoj

# renormalized density of eqn (5a)
d_gradrho[d_idx*3 + 0] += (rhoj - rhoi) * DWIJ[0] * Vj
d_gradrho[d_idx*3 + 1] += (rhoj - rhoi) * DWIJ[1] * Vj
d_gradrho[d_idx*3 + 2] += (rhoj - rhoi) * DWIJ[2] * Vj


class ContinuityEquationDeltaSPH(Equation):
Expand All @@ -386,6 +383,7 @@ class ContinuityEquationDeltaSPH(Equation):
violent impact flows", Computer Methods in Applied Mechanics and
Engineering, 200 (2011), pp 1526--1542.
"""

def __init__(self, dest, sources, c0, delta=0.1):
r"""
Parameters
Expand All @@ -399,29 +397,21 @@ def __init__(self, dest, sources, c0, delta=0.1):
self.delta = delta
super(ContinuityEquationDeltaSPH, self).__init__(dest, sources)

def initialize(self, d_idx, d_arho):
d_arho[d_idx] = 0.0

def loop(self, d_idx, d_arho, s_idx, s_m, d_cs, s_cs, d_rho, s_rho,
DWIJ, VIJ, XIJ, RIJ, HIJ, EPS):

def loop(self, d_idx, d_arho, s_idx, s_m, d_rho, s_rho, DWIJ, XIJ, R2IJ,
HIJ, EPS, d_gradrho, s_gradrho):
rhoi = d_rho[d_idx]
rhoj = s_rho[s_idx]
Vj = s_m[s_idx]/rhoj

# v_{ij} \cdot \nabla W
vijdotdwij = DWIJ[0]*VIJ[0] + DWIJ[1]*VIJ[1] + DWIJ[2]*VIJ[2]

# eta_{ij} \cdot \nabla W
etadotdwij = XIJ[0]*DWIJ[0] + XIJ[1]*DWIJ[1] + XIJ[2]*DWIJ[2]
etadotdwij /= (RIJ + EPS)
fac = -2.0 * (rhoj - rhoi)/(R2IJ + EPS)
psix = fac*XIJ[0] - d_gradrho[d_idx*3 + 0] - s_gradrho[s_idx*3 + 0]
psiy = fac*XIJ[1] - d_gradrho[d_idx*3 + 1] - s_gradrho[s_idx*3 + 1]
psiz = fac*XIJ[2] - d_gradrho[d_idx*3 + 2] - s_gradrho[s_idx*3 + 2]

# celerity (sound speed)
cij = self.c0
psi_ij = self.delta * HIJ * cij * (rhoj - rhoi)
psidotdwij = psix*DWIJ[0] + psiy*DWIJ[1] + psiz*DWIJ[2]

# standard term with dissipative penalization eqn (5a)
d_arho[d_idx] += rhoi*vijdotdwij*Vj + psi_ij*etadotdwij*Vj
d_arho[d_idx] += self.delta * HIJ * self.c0 * psidotdwij * Vj


class UpdateSmoothingLengthFerrari(Equation):
Expand Down
47 changes: 40 additions & 7 deletions pysph/sph/wc/viscosity.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,12 @@

from pysph.sph.equation import Equation


class LaminarViscosity(Equation):
def __init__(self, dest, sources, nu, eta=0.01):
self.nu = nu
self.eta = eta
super(LaminarViscosity,self).__init__(dest, sources)
super(LaminarViscosity, self).__init__(dest, sources)

def loop(self, d_idx, s_idx, s_m, d_rho, s_rho,
d_au, d_av, d_aw, DWIJ, XIJ, VIJ, R2IJ, HIJ):
Expand All @@ -18,17 +19,18 @@ def loop(self, d_idx, s_idx, s_m, d_rho, s_rho,

mb = s_m[s_idx]

tmp = mb * 4 * self.nu * Fij/( (rhoa + rhob)*(R2IJ + self.eta*HIJ*HIJ) )
tmp = mb * 4 * self.nu * Fij/((rhoa + rhob)*(R2IJ + self.eta*HIJ*HIJ))

# accelerations
d_au[d_idx] += tmp * VIJ[0]
d_av[d_idx] += tmp * VIJ[1]
d_aw[d_idx] += tmp * VIJ[2]


class MonaghanSignalViscosityFluids(Equation):
def __init__(self, dest, sources, alpha, h):
self.alpha=0.125 * alpha * h
super(MonaghanSignalViscosityFluids,self).__init__(dest, sources)
self.alpha = 0.125 * alpha * h
super(MonaghanSignalViscosityFluids, self).__init__(dest, sources)

def loop(self, d_idx, s_idx, d_rho, s_rho, s_m,
d_au, d_av, d_aw, d_cs, s_cs,
Expand All @@ -44,7 +46,8 @@ def loop(self, d_idx, s_idx, d_rho, s_rho, s_m,

vabdotrab = VIJ[0]*XIJ[0] + VIJ[1]*XIJ[1] + VIJ[2]*XIJ[2]

force = -16 * nua * nub/(nua*rhoa + nub*rhob) * vabdotrab/(HIJ * (RIJ + 0.01*HIJ*HIJ))
eta = nua * nub/(nua*rhoa + nub*rhob)
force = -16 * eta * vabdotrab/(HIJ * (RIJ + 0.01*HIJ*HIJ))

d_au[d_idx] += -mb * force * DWIJ[0]
d_av[d_idx] += -mb * force * DWIJ[1]
Expand Down Expand Up @@ -99,10 +102,40 @@ def loop(self, d_idx, s_idx, d_m, s_m, d_rho, s_rho, d_h, s_h, d_cs, s_cs,
dot = VIJ[0]*XIJ[0] + VIJ[1]*XIJ[1] + VIJ[2]*XIJ[2]

# Pi_ab term. Eq. (8.9) in [JM05]
rhoa = d_rho[d_idx]; rhob = s_rho[s_idx]
piab = -s_m[s_idx] * self.factor*mua*mub/(rhoa*rhob*(mua + mub)) * (dot/(R2IJ + EPS))
rhoa = d_rho[d_idx]
rhob = s_rho[s_idx]
eta = mua*mub/(rhoa*rhob*(mua + mub))
piab = -s_m[s_idx] * self.factor*eta * (dot/(R2IJ + EPS))

# accelerations due to viscosity Eq. (8.2) in [JM05]
d_au[d_idx] += piab * DWIJ[0]
d_av[d_idx] += piab * DWIJ[1]
d_aw[d_idx] += piab * DWIJ[2]


class LaminarViscosityDeltaSPH(Equation):
r"""
See section 2 of the below reference
- P. 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.
"""
def __init__(self, dest, sources, dim, rho0, nu):
self.dim = dim
self.rho0 = rho0
self.nu = nu
super(LaminarViscosityDeltaSPH, self).__init__(dest, sources)

def loop(self, d_idx, s_idx, s_m, s_rho, d_rho, d_au, d_av, d_aw, HIJ,
DWIJ, R2IJ, EPS, VIJ, XIJ):
Vj = s_m[s_idx]/s_rho[s_idx]
rhoi = d_rho[d_idx]
vdotxij = VIJ[0]*XIJ[0] + VIJ[1]*XIJ[1] + VIJ[2]*XIJ[2]
piij = vdotxij/(R2IJ + EPS)

fac = 2 * (self.dim + 2) * self.nu * self.rho0 * piij * Vj / rhoi
d_au[d_idx] += fac*DWIJ[0]
d_av[d_idx] += fac*DWIJ[1]
d_aw[d_idx] += fac*DWIJ[2]

0 comments on commit 0b273ee

Please sign in to comment.