Added 3rd order softened derivatives

src/gravity/Default/gravity.h

@@ -75,8 +75,10 @@ gravity_compute_timestep_self(const struct gpart* const gp,
 
   const float ac_inv = (ac2 > 0.f) ? 1.f / sqrtf(ac2) : FLT_MAX;
 
+  const float epsilon = gravity_get_softening(gp);
+
   /* Note that 0.66666667 = 2. (from Gadget) / 3. (Plummer softening) */
-  const float dt = sqrtf(0.66666667f * grav_props->eta * gp->epsilon * ac_inv);
+  const float dt = sqrtf(0.66666667f * grav_props->eta * epsilon * ac_inv);
 
   return dt;
 }
@@ -165,7 +167,7 @@ __attribute__((always_inline)) INLINE static void gravity_init_softening(
     struct gpart* gp, const struct gravity_props* grav_props) {
 
   /* Note 3 is the Plummer-equivalent correction */
-  gp->epsilon = 3.f * grav_props->epsilon;
+  gp->epsilon = grav_props->epsilon;
 }
 
 #endif /* SWIFT_DEFAULT_GRAVITY_H */

src/gravity/Default/gravity_iact.h

@@ -64,6 +64,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
 
     kernel_grav_eval(ui, &W);
 
+    if (W < 0) error("W < 0");
+
     /* Get softened gravity */
     fi = mj * hi_inv3 * W * f_lr;
   }
@@ -81,6 +83,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
 
     kernel_grav_eval(uj, &W);
 
+    if (W < 0) error("W < 0");
+
     /* Get softened gravity */
     fj = mi * hj_inv3 * W * f_lr;
   }
@@ -131,6 +135,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
 
     kernel_grav_eval(ui, &W);
 
+    if (W < 0) error("W < 0");
+
     /* Get softened gravity */
     f = mj * hi_inv3 * W * f_lr;
   }

src/gravity_properties.c

@@ -52,7 +52,7 @@ void gravity_props_init(struct gravity_props *p,
   p->theta_crit_inv = 1. / p->theta_crit;
 
   /* Softening lengths */
-  p->epsilon = parser_get_param_double(params, "Gravity:epsilon");
+  p->epsilon = 3. * parser_get_param_double(params, "Gravity:epsilon");
   p->epsilon2 = p->epsilon * p->epsilon;
   p->epsilon_inv = 1. / p->epsilon;
 }
@@ -66,7 +66,8 @@ void gravity_props_print(const struct gravity_props *p) {
 
   message("Self-gravity opening angle:  theta=%.4f", p->theta_crit);
 
-  message("Self-gravity softening:    epsilon=%.4f", p->epsilon);
+  message("Self-gravity softening:    epsilon=%.4f (Plummer equivalent: %.4f)",
+          p->epsilon, p->epsilon / 3.);
 
   if (p->a_smooth != gravity_props_default_a_smooth)
     message("Self-gravity MM smoothing-scale: a_smooth=%f", p->a_smooth);
@@ -81,6 +82,8 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
 
   io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta);
   io_write_attribute_f(h_grpgrav, "Softening length", p->epsilon);
+  io_write_attribute_f(h_grpgrav, "Softening length (Plummer equivalent)",
+                       p->epsilon / 3.);
   io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
   io_write_attribute_d(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
   io_write_attribute_f(h_grpgrav, "MM a_smooth", p->a_smooth);

src/gravity_softened_derivatives.h

@@ -123,8 +123,10 @@ __attribute__((always_inline)) INLINE static double D_soft_200(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
 
   const double u = r * eps_inv;
-  return r_x * r_x * eps_inv * eps_inv * D_soft_2(u) +
-         (r_y * r_y + r_z * r_z) * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv3 = eps_inv * eps_inv2;
+  const double eps_inv5 = eps_inv3 * eps_inv2;
+  return r_x * r_x * eps_inv5 * D_soft_2(u) - eps_inv3 * D_soft_1(u);
 }
 
 /**
@@ -140,8 +142,10 @@ __attribute__((always_inline)) INLINE static double D_soft_020(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
 
   const double u = r * eps_inv;
-  return r_y * r_y * eps_inv * eps_inv * D_soft_2(u) +
-         (r_x * r_x + r_z * r_z) * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv3 = eps_inv * eps_inv2;
+  const double eps_inv5 = eps_inv3 * eps_inv2;
+  return r_y * r_y * eps_inv5 * D_soft_2(u) - eps_inv3 * D_soft_1(u);
 }
 
 /**
@@ -157,8 +161,10 @@ __attribute__((always_inline)) INLINE static double D_soft_002(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
 
   const double u = r * eps_inv;
-  return r_z * r_z * eps_inv * eps_inv * D_soft_2(u) +
-         (r_x * r_x + r_y * r_y) * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv3 = eps_inv * eps_inv2;
+  const double eps_inv5 = eps_inv3 * eps_inv2;
+  return r_z * r_z * eps_inv5 * D_soft_2(u) - eps_inv3 * D_soft_1(u);
 }
 
 /**
@@ -175,8 +181,9 @@ __attribute__((always_inline)) INLINE static double D_soft_110(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
 
   const double u = r * eps_inv;
-  return r_x * r_y * eps_inv * eps_inv * D_soft_2(u) -
-         r_x * r_y * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
+  return r_x * r_y * eps_inv5 * D_soft_2(u);
 }
 
 /**
@@ -193,8 +200,9 @@ __attribute__((always_inline)) INLINE static double D_soft_101(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
 
   const double u = r * eps_inv;
-  return r_x * r_z * eps_inv * eps_inv * D_soft_2(u) -
-         r_x * r_z * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
+  return r_x * r_z * eps_inv5 * D_soft_2(u);
 }
 
 /**
@@ -211,8 +219,9 @@ __attribute__((always_inline)) INLINE static double D_soft_011(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
 
   const double u = r * eps_inv;
-  return r_y * r_z * eps_inv * eps_inv * D_soft_2(u) -
-         r_y * r_z * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
+  return r_y * r_z * eps_inv5 * D_soft_2(u);
 }
 
 /*************************/
@@ -232,7 +241,11 @@ __attribute__((always_inline)) INLINE static double D_soft_300(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
 
   const double u = r * eps_inv;
-  return u * 0;
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
+  const double eps_inv7 = eps_inv5 * eps_inv2;
+  return -r_x * r_x * r_x * eps_inv7 * D_soft_3(u) +
+         3. * r_x * eps_inv5 * D_soft_2(u);
 }
 
 /**
@@ -246,8 +259,13 @@ __attribute__((always_inline)) INLINE static double D_soft_300(
  */
 __attribute__((always_inline)) INLINE static double D_soft_030(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
+
   const double u = r * eps_inv;
-  return u * 0;
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
+  const double eps_inv7 = eps_inv5 * eps_inv2;
+  return -r_y * r_y * r_y * eps_inv7 * D_soft_3(u) +
+         3. * r_y * eps_inv5 * D_soft_2(u);
 }
 
 /**
@@ -261,8 +279,13 @@ __attribute__((always_inline)) INLINE static double D_soft_030(
  */
 __attribute__((always_inline)) INLINE static double D_soft_003(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
+
   const double u = r * eps_inv;
-  return u * 0;
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
+  const double eps_inv7 = eps_inv5 * eps_inv2;
+  return -r_z * r_z * r_z * eps_inv7 * D_soft_3(u) +
+         3. * r_z * eps_inv5 * D_soft_2(u);
 }
 
 /**
@@ -278,8 +301,13 @@ __attribute__((always_inline)) INLINE static double D_soft_003(
  */
 __attribute__((always_inline)) INLINE static double D_soft_210(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
+
   const double u = r * eps_inv;
-  return u * 0;
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
+  const double eps_inv7 = eps_inv5 * eps_inv2;
+  return -r_x * r_x * r_y * eps_inv7 * D_soft_3(u) +
+         r_y * eps_inv5 * D_soft_2(u);
 }
 
 /**
@@ -295,8 +323,13 @@ __attribute__((always_inline)) INLINE static double D_soft_210(
  */
 __attribute__((always_inline)) INLINE static double D_soft_201(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
+
   const double u = r * eps_inv;
-  return u * 0;
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
+  const double eps_inv7 = eps_inv5 * eps_inv2;
+  return -r_x * r_x * r_z * eps_inv7 * D_soft_3(u) +
+         r_z * eps_inv5 * D_soft_2(u);
 }
 
 /**
@@ -312,8 +345,13 @@ __attribute__((always_inline)) INLINE static double D_soft_201(
  */
 __attribute__((always_inline)) INLINE static double D_soft_120(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
+
   const double u = r * eps_inv;
-  return u * 0;
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
+  const double eps_inv7 = eps_inv5 * eps_inv2;
+  return -r_x * r_y * r_y * eps_inv7 * D_soft_3(u) +
+         r_x * eps_inv5 * D_soft_2(u);
 }
 
 /**
@@ -329,8 +367,13 @@ __attribute__((always_inline)) INLINE static double D_soft_120(
  */
 __attribute__((always_inline)) INLINE static double D_soft_021(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
+
   const double u = r * eps_inv;
-  return u * 0;
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
+  const double eps_inv7 = eps_inv5 * eps_inv2;
+  return -r_y * r_y * r_z * eps_inv7 * D_soft_3(u) +
+         r_z * eps_inv5 * D_soft_2(u);
 }
 
 /**
@@ -346,8 +389,13 @@ __attribute__((always_inline)) INLINE static double D_soft_021(
  */
 __attribute__((always_inline)) INLINE static double D_soft_102(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
+
   const double u = r * eps_inv;
-  return u * 0;
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
+  const double eps_inv7 = eps_inv5 * eps_inv2;
+  return -r_x * r_z * r_z * eps_inv7 * D_soft_3(u) +
+         r_x * eps_inv5 * D_soft_2(u);
 }
 
 /**
@@ -363,8 +411,13 @@ __attribute__((always_inline)) INLINE static double D_soft_102(
  */
 __attribute__((always_inline)) INLINE static double D_soft_012(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
+
   const double u = r * eps_inv;
-  return u * 0;
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
+  const double eps_inv7 = eps_inv5 * eps_inv2;
+  return -r_y * r_z * r_z * eps_inv7 * D_soft_3(u) +
+         r_y * eps_inv5 * D_soft_2(u);
 }
 
 /**
@@ -379,8 +432,12 @@ __attribute__((always_inline)) INLINE static double D_soft_012(
  */
 __attribute__((always_inline)) INLINE static double D_soft_111(
     double r_x, double r_y, double r_z, double r, double eps_inv) {
+
   const double u = r * eps_inv;
-  return u * 0;
+  const double eps_inv2 = eps_inv * eps_inv;
+  const double eps_inv4 = eps_inv2 * eps_inv2;
+  const double eps_inv7 = eps_inv4 * eps_inv2 * eps_inv;
+  return -r_x * r_y * r_z * eps_inv7 * D_soft_3(u);
 }
 
 #endif /* SWIFT_GRAVITY_SOFTENED_DERIVATIVE_H */

src/kernel_gravity.h

@@ -73,62 +73,44 @@ __attribute__((always_inline)) INLINE static void kernel_grav_eval_double(
 
 __attribute__((always_inline)) INLINE static double D_soft_0(double u) {
 
-  /* phi(u) = 3u^7 - 15u^6 + 28u^5 - 21u^4 + 7u^2 - 3 */
-  double phi = 3. * u - 15.;
-  phi = phi * u + 28.;
-  phi = phi * u - 21.;
+  /* phi(u) = -3u^7 + 15u^6 - 28u^5 + 21u^4 - 7u^2 + 3 */
+  double phi = -3. * u + 15.;
+  phi = phi * u - 28.;
+  phi = phi * u + 21.;
   phi = phi * u;
-  phi = phi * u + 7.;
+  phi = phi * u - 7.;
   phi = phi * u;
-  phi = phi * u - 3.;
+  phi = phi * u + 3.;
 
   return phi;
 }
 
 __attribute__((always_inline)) INLINE static double D_soft_1(double u) {
 
-  /* phi'(u) = 21u^6 - 90u^5 + 140u^4 - 84u^3 + 14u */
+  /* phi'(u)/u = 21u^5 - 90u^4 + 140u^3 - 84u^2 + 14 */
   double phi = 21. * u - 90.;
   phi = phi * u + 140.;
   phi = phi * u - 84.;
   phi = phi * u;
   phi = phi * u + 14.;
-  phi = phi * u;
 
   return phi;
 }
 
 __attribute__((always_inline)) INLINE static double D_soft_2(double u) {
 
-  /* phi''(u) = 126u^5 - 450u^4 + 560u^3 - 252u^2 + 14 */
-  double phi = 126. * u - 450.;
-  phi = phi * u + 560.;
-  phi = phi * u - 252.;
-  phi = phi * u;
-  phi = phi * u + 14.;
+  /* (phi'(u)/u)'/u = -105u^3 + 360u^2 - 420u + 168 */
+  double phi = -105. * u + 360.;
+  phi = phi * u - 420.;
+  phi = phi * u + 168.;
 
   return phi;
 }
 
 __attribute__((always_inline)) INLINE static double D_soft_3(double u) {
 
-  /* phi'''(u) = 630u^4 - 1800u^3 + 1680u^2 - 504u */
-  double phi = 630. * u - 1800.;
-  phi = phi * u + 1680.;
-  phi = phi * u - 504.;
-  phi = phi * u;
-
-  return phi;
-}
-
-__attribute__((always_inline)) INLINE static double D_soft_4(double u) {
-
-  /* phi''''(u) = 2520u^3 - 5400u^2 + 3360u - 504 */
-  double phi = 2520. * u - 5400.;
-  phi = phi * u + 3360.;
-  phi = phi * u - 504.;
-
-  return phi;
+  /* ((phi'(u)/u)'/u)'/u = 315u - 720 + 420/u */
+  return 315. * u - 720. + 420. / u;
 }
 
 #endif /* SWIFT_KERNEL_GRAVITY_H */

theory/Multipoles/potential_derivatives.tex

@@ -19,7 +19,7 @@ we can construct the higher order terms by successively applying the
 \begin{align}
 D_{100}(\mathbf{r}) = \frac{\partial}{\partial r_x} \phi (\mathbf{r},H) = 
 \left\lbrace\begin{array}{rcl}
-\frac{r_x}{H^3} \left(-21u^5 + 90u^4 - 140u^3 + 84u^2 - 14\right) & \mbox{if} & u < 1,\\
+-\frac{r_x}{H^3} \left(21u^5 - 90u^4 + 140u^3 - 84u^2 + 14\right) & \mbox{if} & u < 1,\\
 -\frac{r_x}{r^3} & \mbox{if} & u \geq 1, 
 \end{array}
 \right.\nonumber
@@ -28,8 +28,8 @@ D_{100}(\mathbf{r}) = \frac{\partial}{\partial r_x} \phi (\mathbf{r},H) =
 \begin{align}
 D_{200}(\mathbf{r}) = \frac{\partial^2}{\partial r_x^2} \phi (\mathbf{r},H) = 
 \left\lbrace\begin{array}{rcl}
-\frac{r_x^2}{H^5}\left(-105u^3+360u^2-420u+168\right) +
-\frac{1}{H^3} \left(-21u^5 + 90u^4 - 140u^3 + 84u^2 - 14\right) & \mbox{if} & u < 1,\\
+\frac{r_x^2}{H^5}\left(-105u^3+360u^2-420u+168\right) -
+\frac{1}{H^3} \left(21u^5 - 90u^4 + 140u^3 - 84u^2 + 14\right) & \mbox{if} & u < 1,\\
 3\frac{r_x^2}{r^5} - \frac{1}{r^3} & \mbox{if} & u \geq 1, 
 \end{array}
 \right.\nonumber
@@ -43,3 +43,33 @@ D_{110}(\mathbf{r}) = \frac{\partial^2}{\partial r_x\partial r_y} \phi (\mathbf{
 \end{array}
 \right.\nonumber
 \end{align}
+
+\begin{align}
+D_{300}(\mathbf{r}) = \frac{\partial^3}{\partial r_x^3} \phi (\mathbf{r},H) = 
+\left\lbrace\begin{array}{rcl}
+-\frac{r_x^3}{H^7} \left(315u - 720 + 420u^{-1}\right) +
+\frac{3r_x}{H^5}\left(-105u^3+360u^2-420u+168\right) & \mbox{if} & u < 1,\\
+-15\frac{r_x^3}{r^7} + 9 \frac{r_x}{r^5} & \mbox{if} & u \geq 1, 
+\end{array}
+\right.\nonumber
+\end{align}
+
+\begin{align}
+D_{210}(\mathbf{r}) = \frac{\partial^3}{\partial r_x^3} \phi (\mathbf{r},H) = 
+\left\lbrace\begin{array}{rcl}
+-\frac{r_x^2r_y}{H^7} \left(315u - 720 + 420u^{-1}\right) +
+\frac{r_y}{H^5}\left(-105u^3+360u^2-420u+168\right) & \mbox{if} & u < 1,\\
+-15\frac{r_x^2r_y}{r^7} + 3 \frac{r_y}{r^5} & \mbox{if} & u \geq 1, 
+\end{array}
+\right.\nonumber
+\end{align}
+
+
+\begin{align}
+D_{111}(\mathbf{r}) = \frac{\partial^3}{\partial r_x\partial r_y\partial r_z} \phi (\mathbf{r},H) = 
+\left\lbrace\begin{array}{rcl}
+-\frac{r_xr_yr_z}{H^7} \left(315u - 720 + 420u^{-1}\right) & \mbox{if} & u < 1,\\
+-15\frac{r_xr_yr_z}{r^7} & \mbox{if} & u \geq 1, 
+\end{array}
+\right.\nonumber
+\end{align}