Pass the gravity properties to the M2L kernel. Better documentation.

src/gravity_properties.c

@@ -48,7 +48,9 @@ void gravity_props_init(struct gravity_props *p,
   p->eta = parser_get_param_float(params, "Gravity:eta");
 
   /* Softening lengths */
-  p->epsilon = parser_get_param_float(params, "Gravity:epsilon");
+  p->epsilon = parser_get_param_double(params, "Gravity:epsilon");
+  p->epsilon2 = p->epsilon * p->epsilon;
+  p->epsilon_inv = 1. / p->epsilon;
 }
 
 void gravity_props_print(const struct gravity_props *p) {

src/gravity_properties.h

@@ -38,11 +38,17 @@ struct gravity_props {
   float a_smooth;
   float r_cut;
 
-  /* Time integration parameters */
+  /*! Time integration dimensionless multiplier */
   float eta;
 
-  /* Softening lengths */
-  float epsilon;
+  /*! Softening length */
+  double epsilon;
+
+  /*! Square of softening length */
+  double epsilon2;
+
+  /*! Inverse of softening length */
+  double epsilon_inv;
 };
 
 void gravity_props_print(const struct gravity_props *p);

src/multipole.h

@@ -33,6 +33,7 @@
 #include "const.h"
 #include "error.h"
 #include "gravity_derivatives.h"
+#include "gravity_properties.h"
 #include "inline.h"
 #include "kernel_gravity.h"
 #include "minmax.h"
@@ -164,20 +165,23 @@ struct multipole {
  */
 struct gravity_tensors {
 
-  /*! Linking pointer for "memory management". */
-  struct gravity_tensors *next;
+  union {
 
-  /*! Centre of mass of the matter dsitribution */
-  double CoM[3];
+    /*! Linking pointer for "memory management". */
+    struct gravity_tensors *next;
 
-  /*! The actual content */
-  struct {
+    /*! The actual content */
+    struct {
 
-    /*! Multipole mass */
-    struct multipole m_pole;
+      /*! Centre of mass of the matter dsitribution */
+      double CoM[3];
 
-    /*! Field tensor for the potential */
-    struct grav_tensor pot;
+      /*! Multipole mass */
+      struct multipole m_pole;
+
+      /*! Field tensor for the potential */
+      struct grav_tensor pot;
+    };
   };
 } SWIFT_STRUCT_ALIGN;
 
@@ -206,6 +210,11 @@ INLINE static void gravity_drift(struct gravity_tensors *m, double dt) {
   m->CoM[2] += m->m_pole.vel[2] * dt;
 }
 
+/**
+ * @brief Zeroes all the fields of a field tensor
+ *
+ * @param l The field tensor.
+ */
 INLINE static void gravity_field_tensors_init(struct grav_tensor *l) {
 
   bzero(l, sizeof(struct grav_tensor));
@@ -354,6 +363,11 @@ INLINE static void gravity_field_tensors_print(const struct grav_tensor *l) {
 #endif
 }
 
+/**
+ * @brief Zeroes all the fields of a multipole.
+ *
+ * @param m The multipole
+ */
 INLINE static void gravity_multipole_init(struct multipole *m) {
 
   bzero(m, sizeof(struct multipole));
@@ -1461,11 +1475,13 @@ INLINE static void gravity_M2M(struct multipole *m_a,
  * @param m_a The multipole creating the field.
  * @param pos_b The position of the field tensor.
  * @param pos_a The position of the multipole.
+ * @param props The #gravity_props of this calculation.
  * @param periodic Is the calculation periodic ?
  */
 INLINE static void gravity_M2L(struct grav_tensor *l_b,
                                const struct multipole *m_a,
                                const double pos_b[3], const double pos_a[3],
+                               const struct gravity_props *props,
                                int periodic) {
 
   double dx, dy, dz;

src/runner_doiact_grav.h

@@ -93,6 +93,7 @@ void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
                            struct cell *restrict cj) {
 
   const struct engine *e = r->e;
+  const struct gravity_props *props = e->gravity_properties;
   const int periodic = e->s->periodic;
   const struct multipole *multi_j = &cj->multipole->m_pole;
   // const float a_smooth = e->gravity_properties->a_smooth;
@@ -114,7 +115,7 @@ void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
 
   /* Let's interact at this level */
   gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
-              cj->multipole->CoM, periodic * 0);
+              cj->multipole->CoM, props, periodic * 0);
 
   TIMER_TOC(timer_dopair_grav_mm);
 }