from pysph.sph.equation import Equation from pysph.sph.integrator_step import IntegratorStep from math import sqrt class SummationDensity(Equation): def initialize(self, d_idx, d_V, d_rho): d_V[d_idx] = 0.0 d_rho[d_idx] = 0.0 def loop(self, d_idx, d_V, WIJ): d_V[d_idx] += WIJ def post_loop(self, d_idx, d_rho, d_V, d_m, d_max_rho, d_min_rho): d_rho[d_idx] = d_m[d_idx] * d_V[d_idx] if d_max_rho[d_idx] < d_rho[d_idx]: d_max_rho[d_idx] = d_rho[d_idx] elif d_min_rho[d_idx] > d_rho[d_idx]: d_min_rho[d_idx] = d_rho[d_idx] class StateEquation(Equation): def __init__(self, dest, sources, c_l, c_g, rho_l, rho_g, p0): self.c_l = c_l self.c_g = c_g self.rho_l = rho_l self.rho_g = rho_g self.p0 = p0 super(StateEquation, self).__init__(dest, sources) def loop(self, d_idx, d_p, d_rho, d_phase): if d_phase[d_idx] == 0.0: d_p[d_idx] = self.c_g * self.c_g * (d_rho[d_idx] - self.rho_g) + self.p0 else: d_p[d_idx] = self.c_l * self.c_l * (d_rho[d_idx] - self.rho_l) + self.p0 def post_loop(self, d_idx, d_p, d_min_p): if d_min_p[d_idx] > d_p[d_idx]: d_min_p[d_idx] = d_p[d_idx] class AdamiColorGradient(Equation): def initialize(self, d_idx, d_cx, d_cy, d_cz, d_nx, d_ny, d_nz, d_N, d_ddelta): d_cx[d_idx] = 0.0 d_cy[d_idx] = 0.0 d_cz[d_idx] = 0.0 d_nx[d_idx] = 0.0 d_ny[d_idx] = 0.0 d_nz[d_idx] = 0.0 # reliability indicator for normals d_N[d_idx] = 0.0 # Discretized dirac-delta d_ddelta[d_idx] = 0.0 def loop(self, d_idx, s_idx, d_V, s_V, d_rho, s_rho, d_cx, d_cy, d_cz, d_phase, s_phase, DWIJ): if d_phase[d_idx] != s_phase[s_idx]: V_ij = 1.0 / (d_V[d_idx] * d_V[d_idx]) + 1.0 / (s_V[s_idx] * s_V[s_idx]) rho_ij = d_rho[d_idx] / (d_rho[d_idx] + s_rho[s_idx]) d_cx[d_idx] += d_V[d_idx] * V_ij * rho_ij * DWIJ[0] d_cy[d_idx] += d_V[d_idx] * V_ij * rho_ij * DWIJ[1] d_cz[d_idx] += d_V[d_idx] * V_ij * rho_ij * DWIJ[2] def post_loop(self, d_idx, d_cx, d_cy, d_cz, d_h, d_nx, d_ny, d_nz, d_ddelta, d_N): # absolute value of the color gradient mod_gradc2 = d_cx[d_idx] * d_cx[d_idx] + d_cy[d_idx] * d_cy[d_idx] + d_cz[d_idx] * d_cz[d_idx] # avoid sqrt computations on non-interface particles h2 = d_h[d_idx] * d_h[d_idx] if mod_gradc2 > 1e-12 / h2: # this normal is reliable in the sense of [JM00] d_N[d_idx] = 1.0 # compute the normals mod_gradc = sqrt(mod_gradc2) d_nx[d_idx] = d_cx[d_idx] / mod_gradc d_ny[d_idx] = d_cy[d_idx] / mod_gradc d_nz[d_idx] = d_cz[d_idx] / mod_gradc # discretized dirac delta d_ddelta[d_idx] = mod_gradc else: d_cx[d_idx] = 0.0 d_cy[d_idx] = 0.0 d_cz[d_idx] = 0.0 class AdamiReproducingDivergence(Equation): def __init__(self, dest, sources, dim): self.dim = dim super(AdamiReproducingDivergence, self).__init__(dest, sources) def initialize(self, d_idx, d_kappa, d_wij_sum): d_kappa[d_idx] = 0.0 d_wij_sum[d_idx] = 0.0 def loop(self, d_idx, s_idx, d_kappa, d_wij_sum, d_nx, d_ny, d_nz, s_nx, s_ny, s_nz, s_V, d_phase, s_phase, DWIJ, XIJ): if d_phase[d_idx] != s_phase[s_idx]: phi_ij = -1.0 else: phi_ij = 1.0 # dot product in the numerator of Eq. (20) nijdotdwij = (d_nx[d_idx] - phi_ij * s_nx[s_idx]) * DWIJ[0] + (d_ny[d_idx] - phi_ij * s_ny[s_idx]) * DWIJ[1] + (d_nz[d_idx] - phi_ij * s_nz[s_idx]) * DWIJ[2] # dot product in the denominator of Eq. (20) xijdotdwij = XIJ[0] * DWIJ[0] + XIJ[1] * DWIJ[1] + XIJ[2] * DWIJ[2] # accumulate the contributions d_kappa[d_idx] += nijdotdwij / s_V[s_idx] d_wij_sum[d_idx] += xijdotdwij / s_V[s_idx] def post_loop(self, d_idx, d_kappa, d_wij_sum): # normalize the curvature estimate if d_wij_sum[d_idx] > 1e-12: d_kappa[d_idx] /= d_wij_sum[d_idx] d_kappa[d_idx] *= -self.dim class MomentumEquationPressureGradientAdami(Equation): 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_V, s_V, d_au, d_av, d_aw, d_p, s_p, d_m, d_rho, s_rho, DWIJ): p_ij = (d_rho[d_idx] * s_p[s_idx] + s_rho[s_idx] * d_p[d_idx]) / (d_rho[d_idx] + s_rho[s_idx]) V_ij = 1.0 / (d_V[d_idx] * d_V[d_idx]) + 1.0 / (s_V[s_idx] * s_V[s_idx]) d_au[d_idx] += -p_ij * V_ij * DWIJ[0] / d_m[d_idx] d_av[d_idx] += -p_ij * V_ij * DWIJ[1] / d_m[d_idx] d_aw[d_idx] += -p_ij * V_ij * DWIJ[2] / d_m[d_idx] class MomentumEquationViscosityAdami(Equation): 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_V, s_V, d_au, d_av, d_aw, d_m, d_mu, s_mu, DWIJ, R2IJ, XIJ, EPS, VIJ): mu_ij = 2.0 * d_mu[d_idx] * s_mu[s_idx] / (d_mu[d_idx] + s_mu[s_idx]) V_ij = 1.0 / (d_V[d_idx] * d_V[d_idx]) + 1.0 / (s_V[s_idx] * s_V[s_idx]) dwijdotrij = (DWIJ[0] * XIJ[0] + DWIJ[1] * XIJ[1] + DWIJ[2] * XIJ[2]) / (R2IJ + 0.01 * EPS) factor = mu_ij * V_ij * dwijdotrij / d_m[d_idx] d_au[d_idx] += factor * VIJ[0] d_av[d_idx] += factor * VIJ[1] d_aw[d_idx] += factor * VIJ[2] class CSFSurfaceTensionForceAdami(Equation): 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 post_loop(self, d_idx, d_au, d_av, d_aw, d_kappa, d_cx, d_cy, d_cz, d_sigma, d_rho): d_au[d_idx] += -d_sigma[d_idx] * d_kappa[d_idx] * d_cx[d_idx] / d_rho[d_idx] d_av[d_idx] += -d_sigma[d_idx] * d_kappa[d_idx] * d_cy[d_idx] / d_rho[d_idx] d_aw[d_idx] += -d_sigma[d_idx] * d_kappa[d_idx] * d_cz[d_idx] / d_rho[d_idx] class InterfaceSharpnessForce(Equation): 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_V, s_V, d_au, d_av, d_aw, d_p, s_p, d_m, d_rho, s_rho, d_phase, s_phase, DWIJ): if d_phase[d_idx] != s_phase[s_idx]: p_ij = (d_rho[d_idx] * s_p[s_idx] + s_rho[s_idx] * d_p[d_idx]) / (d_rho[d_idx] + s_rho[s_idx]) V_ij = 1.0 / (d_V[d_idx] * d_V[d_idx]) + 1.0 / (s_V[s_idx] * s_V[s_idx]) d_au[d_idx] += -0.07 * p_ij * V_ij * DWIJ[0] / d_m[d_idx] d_av[d_idx] += -0.07 * p_ij * V_ij * DWIJ[1] / d_m[d_idx] d_aw[d_idx] += -0.07 * p_ij * V_ij * DWIJ[2] / d_m[d_idx] class DropletAcceleration(Equation): def post_loop(self, d_idx, d_av, d_phase, t): if t > 0.35 and t < 0.4: if d_phase[d_idx] == 1.0: d_av[d_idx] += 200.0 elif d_phase[d_idx] == 2.0: d_av[d_idx] -= 200.0 class AdamiVerletStep(IntegratorStep): def initialize(self): pass def stage1(self, d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, d_x, d_y, d_z, dt): dtb2 = 0.5 * dt # velocity predictor eqn (14) d_u[d_idx] += dtb2 * d_au[d_idx] d_v[d_idx] += dtb2 * d_av[d_idx] d_w[d_idx] += dtb2 * d_aw[d_idx] # position predictor eqn (15) d_x[d_idx] += dtb2 * d_u[d_idx] d_y[d_idx] += dtb2 * d_v[d_idx] d_z[d_idx] += dtb2 * d_w[d_idx] def stage2(self, d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, d_x, d_y, d_z, d_vmag2, d_max_vmag2, dt): dtb2 = 0.5 * dt # position corrector eqn (17) d_x[d_idx] += dtb2 * d_u[d_idx] d_y[d_idx] += dtb2 * d_v[d_idx] d_z[d_idx] += dtb2 * d_w[d_idx] # velocity corrector eqn (18) d_u[d_idx] += dtb2 * d_au[d_idx] d_v[d_idx] += dtb2 * d_av[d_idx] d_w[d_idx] += dtb2 * d_aw[d_idx] d_vmag2[d_idx] = d_u[d_idx] * d_u[d_idx] + d_v[d_idx] * d_v[d_idx] + d_w[d_idx] * d_w[d_idx] if d_max_vmag2[d_idx] < d_vmag2[d_idx]: d_max_vmag2[d_idx] = d_vmag2[d_idx]