diff --git a/src/hydro/SPHENIX/hydro.h b/src/hydro/SPHENIX/hydro.h
index 283be043ed93e8aff95727d1f0ae4f6e7ff333b6..c7781df12e30442d546bc872b168b3dc6acb2c01 100644
--- a/src/hydro/SPHENIX/hydro.h
+++ b/src/hydro/SPHENIX/hydro.h
@@ -628,7 +628,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
 
   /* Finish matrix and volume computations for FVPM Radiative Transfer */
   fvpm_compute_volume_and_matrix(p, h_inv_dim);
-  
+
 #ifdef SWIFT_HYDRO_DENSITY_CHECKS
   p->n_density += kernel_root;
   p->n_density *= h_inv_dim;
@@ -694,7 +694,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
       grad_h_term = common_factor * p->density.rho_dh / (1.f + grad_W_term);
     }
   }
-  
+
   /* Update variables. */
   p->force.f = grad_h_term;
   p->force.pressure = pressure_including_floor;
@@ -718,7 +718,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_gradient(
   p->viscosity.rot_v[2] = 0.f;
 
   p->viscosity.div_v = 0.f;
-  
+
   p->viscosity.v_sig = 2.f * p->force.soundspeed;
   p->force.alpha_visc_max_ngb = p->viscosity.alpha;
 }
@@ -734,14 +734,15 @@ __attribute__((always_inline)) INLINE static void hydro_reset_gradient(
  * @param p The particle to act upon.
  */
 __attribute__((always_inline)) INLINE static void hydro_end_gradient(
-    struct part *p, const struct cosmology *cosmo, const struct pressure_floor_props *pressure_floor) {
+    struct part *p, const struct cosmology *cosmo,
+    const struct pressure_floor_props *pressure_floor) {
 
   /* Some smoothing length multiples. */
   const float h = p->h;
   const float h_inv = 1.0f / h;                       /* 1/h */
   const float h_inv_dim = pow_dimension(h_inv);       /* 1/h^d */
   const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
-  
+
   const float rho_inv = 1.f / p->rho;
   const float a_inv2 = cosmo->a2_inv;
 
@@ -753,11 +754,11 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
   /* Finish calculation of the velocity divergence */
   p->viscosity.div_v *= h_inv_dim_plus_one * rho_inv * a_inv2;
   p->viscosity.div_v += cosmo->H * hydro_dimension;
-  
+
   /* Include the extra factors in the del^2 u */
 
   p->diffusion.laplace_u *= 2.f * h_inv_dim_plus_one;
-  
+
 #ifdef SWIFT_HYDRO_DENSITY_CHECKS
   p->n_gradient += kernel_root;
 #endif
@@ -842,7 +843,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
       gas_soundspeed_from_pressure(p->rho, pressure_including_floor);
   const float soundspeed_physical = soundspeed * cosmo->a_factor_sound_speed;
   const float fac_B = cosmo->a_factor_Balsara_eps;
-  
+
   const float sound_crossing_time_inverse =
       soundspeed_physical * kernel_support_physical_inv;
 
@@ -853,18 +854,19 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   /* Compute the norm of div v */
   const float abs_div_v = fabsf(p->viscosity.div_v);
-  
+
   /* Get the squares of the quantities necessary for the Balsara-like switch */
   const float fac_B_2 = fac_B * fac_B;
   const float curl_v_2 = curl_v * curl_v;
   const float abs_div_v_2 = abs_div_v * abs_div_v;
   const float soundspeed_2 = soundspeed * soundspeed;
   const float h_2 = p->h * p->h;
-  
-  /* Compute the Balsara-like switch (Price et a. (2018), eq. 47; a simplified version of eq. 18 in Cullen & Dehnen (2012)) */
-  const float balsara =
-      abs_div_v_2 / (abs_div_v_2 + curl_v_2 + 0.0001f * soundspeed_2 * fac_B_2 / h_2);
-  
+
+  /* Compute the Balsara-like switch (Price et a. (2018), eq. 47; a simplified
+   * version of eq. 18 in Cullen & Dehnen (2012)) */
+  const float balsara = abs_div_v_2 / (abs_div_v_2 + curl_v_2 +
+                                       0.0001f * soundspeed_2 * fac_B_2 / h_2);
+
   /* Construct time differential of div.v implicitly following the ANARCHY spec
    */
   const float div_v_dt =
@@ -877,8 +879,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   /* Source term is only activated if flow is converging (i.e. in the pre-
    * shock region) */
   const float S = p->viscosity.div_v < 0.f
-                      ? balsara * kernel_support_physical * kernel_support_physical *
-                            max(0.f, -1.f * div_v_dt)
+                      ? balsara * kernel_support_physical *
+                            kernel_support_physical * max(0.f, -1.f * div_v_dt)
                       : 0.f;
 
   /* We want the decay to occur based on the thermodynamic properties
diff --git a/src/hydro/SPHENIX/hydro_iact.h b/src/hydro/SPHENIX/hydro_iact.h
index 4a96f22c5123d2f23b7e13bf482f89d242929ba7..f5ea496a0eb614fd14c689efe7502a6838be5864 100644
--- a/src/hydro/SPHENIX/hydro_iact.h
+++ b/src/hydro/SPHENIX/hydro_iact.h
@@ -173,7 +173,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient(
    * and all of it's neighbours */
 
   float dv[3], curlvr[3];
-  
+
   const float r = sqrtf(r2);
   const float r_inv = r ? 1.0f / r : 0.0f;
 
@@ -188,11 +188,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient(
 
   kernel_deval(ui, &wi, &wi_dx);
   kernel_deval(uj, &wj, &wj_dx);
-  
+
   /* Variable smoothing length term */
   const float f_ij = 1.f - pi->force.f / mj;
   const float f_ji = 1.f - pj->force.f / mi;
-  
+
   /* Cosmology terms for the signal velocity */
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
@@ -201,8 +201,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient(
   dv[0] = pi->v[0] - pj->v[0];
   dv[1] = pi->v[1] - pj->v[1];
   dv[2] = pi->v[2] - pj->v[2];
-  const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];  
-  
+  const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
+
   /* Add Hubble flow */
 
   const float dvdr_Hubble = dvdr + a2_Hubble * r2;
@@ -213,7 +213,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient(
   /* Signal velocity */
   const float new_v_sig =
       signal_velocity(dx, pi, pj, mu_ij, const_viscosity_beta, a, mu_0);
-  
+
   /* Update if we need to */
   pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig);
   pj->viscosity.v_sig = max(pj->viscosity.v_sig, new_v_sig);
@@ -238,7 +238,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient(
   pj->viscosity.rot_v[0] += facj * curlvr[0];
   pj->viscosity.rot_v[1] += facj * curlvr[1];
   pj->viscosity.rot_v[2] += facj * curlvr[2];
-  
+
   /* Calculate Del^2 u for the thermal diffusion coefficient. */
   const float delta_u_factor = (pi->u - pj->u) * r_inv;
   pi->diffusion.laplace_u += pj->mass * delta_u_factor * wi_dx / pj->rho;
@@ -297,12 +297,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
   float wi, wi_dx;
 
   const float ui = r / hi;
-  
+
   kernel_deval(ui, &wi, &wi_dx);
 
   /* Variable smoothing length term */
   const float f_ij = 1.f - pi->force.f / mj;
-  
+
   /* Cosmology terms for the signal velocity */
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
@@ -312,7 +312,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
   dv[1] = pi->v[1] - pj->v[1];
   dv[2] = pi->v[2] - pj->v[2];
   const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
-  
+
   /* Add Hubble flow */
 
   const float dvdr_Hubble = dvdr + a2_Hubble * r2;
@@ -340,7 +340,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
   pi->viscosity.rot_v[0] += faci * curlvr[0];
   pi->viscosity.rot_v[1] += faci * curlvr[1];
   pi->viscosity.rot_v[2] += faci * curlvr[2];
-  
+
   /* Calculate Del^2 u for the thermal diffusion coefficient. */
   const float delta_u_factor = (pi->u - pj->u) * r_inv;
   pi->diffusion.laplace_u += pj->mass * delta_u_factor * wi_dx / pj->rho;
@@ -431,8 +431,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   /* Construct the full viscosity term */
   const float rho_ij = rhoi + rhoj;
   const float alpha = pi->viscosity.alpha + pj->viscosity.alpha;
-  const float visc =
-      -0.5f * alpha * v_sig * mu_ij / rho_ij;
+  const float visc = -0.5f * alpha * v_sig * mu_ij / rho_ij;
 
   /* Convolve with the kernel */
   const float visc_acc_term =
@@ -579,8 +578,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   /* Construct the full viscosity term */
   const float rho_ij = rhoi + rhoj;
   const float alpha = pi->viscosity.alpha + pj->viscosity.alpha;
-  const float visc =
-      -0.5f * alpha * v_sig * mu_ij / rho_ij;
+  const float visc = -0.5f * alpha * v_sig * mu_ij / rho_ij;
 
   /* Convolve with the kernel */
   const float visc_acc_term =
diff --git a/src/hydro/SPHENIX/hydro_io.h b/src/hydro/SPHENIX/hydro_io.h
index 211ef38e9968a2737d3aa3d743b7701c27922fa5..f4c3757eb8f725923c10a8a3c5953b54d71218c3 100644
--- a/src/hydro/SPHENIX/hydro_io.h
+++ b/src/hydro/SPHENIX/hydro_io.h
@@ -235,8 +235,7 @@ INLINE static void hydro_write_particles(const struct part* parts,
 
   list[9] = io_make_output_field_convert_part(
       "ViscosityParameters", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, parts, xparts,
-      convert_viscosity,
-      "Visosity coefficient (alpha_visc) of the particles");
+      convert_viscosity, "Visosity coefficient (alpha_visc) of the particles");
 
   list[10] = io_make_output_field_convert_part(
       "DiffusionParameters", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, parts, xparts,
diff --git a/src/hydro/SPHENIX/hydro_part.h b/src/hydro/SPHENIX/hydro_part.h
index 859176e34823793d2e54c2646e494e5c71351bb5..67f83e187f1da44eb24e9939e486834dcfd243fe 100644
--- a/src/hydro/SPHENIX/hydro_part.h
+++ b/src/hydro/SPHENIX/hydro_part.h
@@ -141,7 +141,7 @@ struct part {
 
     /*! Particle velocity curl. */
     float rot_v[3];
-    
+
     /*! Particle velocity divergence from previous step */
     float div_v_previous_step;
 
@@ -212,7 +212,7 @@ struct part {
     float alpha_visc_max_ngb;
 
   } force;
-  
+
   /*! Additional data used for adaptive softening */
   struct adaptive_softening_part_data adaptive_softening_data;
 
diff --git a/src/mhd/DirectInduction/mhd.h b/src/mhd/DirectInduction/mhd.h
index a4bddf0b75eb6e4f79c36e7a7abfdaa9e332d582..71b3f8ec96282de08380ee342bab25247b6f825e 100644
--- a/src/mhd/DirectInduction/mhd.h
+++ b/src/mhd/DirectInduction/mhd.h
@@ -291,7 +291,6 @@ __attribute__((always_inline)) INLINE static void mhd_prepare_gradient(
     const float mu_0) {
 
   p->mhd_data.Alfven_speed = mhd_get_comoving_Alfven_speed(p, mu_0);
-
 }
 
 /**
diff --git a/src/mhd/DirectInduction/mhd_iact.h b/src/mhd/DirectInduction/mhd_iact.h
index 570a45db1e81c6012d5eff88995089c9ccfefcec..a41639c2f49270a620898c1ffe5b1a81d4908aa9 100644
--- a/src/mhd/DirectInduction/mhd_iact.h
+++ b/src/mhd/DirectInduction/mhd_iact.h
@@ -481,18 +481,22 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
   const float resistive_eta_j = pj->mhd_data.resistive_eta;
 
   const float rho_term_PR = 1.0f / (rhoi * rhoj);
-  const float grad_term_PR = f_ij * wi_dr + f_ji * wj_dr; 
+  const float grad_term_PR = f_ij * wi_dr + f_ji * wj_dr;
 
   const float dB_dt_pref_PR = rho_term_PR * grad_term_PR * r_inv;
 
   for (int k = 0; k < 3; k++) {
-    pi->mhd_data.B_over_rho_dt[k] += resistive_eta_i * mj * dB_dt_pref_PR * dB[k];
-    pj->mhd_data.B_over_rho_dt[k] -= resistive_eta_j * mi * dB_dt_pref_PR * dB[k];
+    pi->mhd_data.B_over_rho_dt[k] +=
+        resistive_eta_i * mj * dB_dt_pref_PR * dB[k];
+    pj->mhd_data.B_over_rho_dt[k] -=
+        resistive_eta_j * mi * dB_dt_pref_PR * dB[k];
   }
 
-  pi->u_dt -= 0.5f * permeability_inv * resistive_eta_i * mj * dB_dt_pref_PR * dB_2;  
-  pj->u_dt -= 0.5f * permeability_inv * resistive_eta_j * mi * dB_dt_pref_PR * dB_2;  
-  
+  pi->u_dt -=
+      0.5f * permeability_inv * resistive_eta_i * mj * dB_dt_pref_PR * dB_2;
+  pj->u_dt -=
+      0.5f * permeability_inv * resistive_eta_j * mi * dB_dt_pref_PR * dB_2;
+
   /* Artificial resistivity */
   const float art_diff_beta_i = pi->mhd_data.art_diff_beta;
   const float art_diff_beta_j = pj->mhd_data.art_diff_beta;
@@ -563,8 +567,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
     pj->mhd_data.Adv_B_source[i] += mi * dB_dt_pref_j * dB_dt_j[i];
     pi->mhd_data.Adv_B_source[i] -= mj * grad_psi_i * dx[i];
     pj->mhd_data.Adv_B_source[i] -= mi * grad_psi_j * dx[i];
-    pi->mhd_data.Diff_B_source[i] += resistive_eta_i * mj * dB_dt_pref_PR * dB[i];
-    pj->mhd_data.Diff_B_source[i] -= resistive_eta_j * mi * dB_dt_pref_PR * dB[i];
+    pi->mhd_data.Diff_B_source[i] +=
+        resistive_eta_i * mj * dB_dt_pref_PR * dB[i];
+    pj->mhd_data.Diff_B_source[i] -=
+        resistive_eta_j * mi * dB_dt_pref_PR * dB[i];
     pi->mhd_data.Diff_B_source[i] += mj * art_diff_pref_i * dB[i];
     pj->mhd_data.Diff_B_source[i] -= mi * art_diff_pref_j * dB[i];
     pi->mhd_data.Delta_B[i] += mj * dB_dt_pref_Lap_i * wi_dr * dB[i];
@@ -763,16 +769,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
   const float resistive_eta_i = pi->mhd_data.resistive_eta;
 
   const float rho_term_PR = 1.0f / (rhoi * rhoj);
-  const float grad_term_PR = f_ij * wi_dr + f_ji * wj_dr; 
+  const float grad_term_PR = f_ij * wi_dr + f_ji * wj_dr;
 
   const float dB_dt_pref_PR = rho_term_PR * grad_term_PR * r_inv;
 
   for (int k = 0; k < 3; k++) {
-    pi->mhd_data.B_over_rho_dt[k] += resistive_eta_i * mj * dB_dt_pref_PR * dB[k];
+    pi->mhd_data.B_over_rho_dt[k] +=
+        resistive_eta_i * mj * dB_dt_pref_PR * dB[k];
   }
 
-  pi->u_dt -= 0.5f * permeability_inv * resistive_eta_i * mj * dB_dt_pref_PR * dB_2;  
-  
+  pi->u_dt -=
+      0.5f * permeability_inv * resistive_eta_i * mj * dB_dt_pref_PR * dB_2;
+
   /* Artificial resistivity */
   const float art_diff_beta = pi->mhd_data.art_diff_beta;
 
@@ -821,7 +829,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
   for (int i = 0; i < 3; i++) {
     pi->mhd_data.Adv_B_source[i] += mj * dB_dt_pref_i * dB_dt_i[i];
     pi->mhd_data.Adv_B_source[i] -= mj * grad_psi_i * dx[i];
-    pi->mhd_data.Diff_B_source[i] += resistive_eta_i * mj * dB_dt_pref_PR * dB[i];
+    pi->mhd_data.Diff_B_source[i] +=
+        resistive_eta_i * mj * dB_dt_pref_PR * dB[i];
     pi->mhd_data.Diff_B_source[i] += mj * art_diff_pref * dB[i];
     pi->mhd_data.Delta_B[i] += mj * dB_dt_pref_Lap * wi_dr * dB[i];
   }
diff --git a/src/statistics.c b/src/statistics.c
index ebe79b615d1f20d87108ef382c0a96df8c18f9d3..2c1c1c343870c287ed081c839f4101d09f69429f 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -259,7 +259,8 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
     stats.divB_error += mhd_get_divB_error(p, xp);
 
     /* Get maximal div B error */
-    stats.divB_error_max = fmaxf(stats.divB_error_max, mhd_get_divB_error(p, xp));
+    stats.divB_error_max =
+        fmaxf(stats.divB_error_max, mhd_get_divB_error(p, xp));
 
     /* Collect square of magnetic field vector norm */
     stats.Brms += mhd_get_Bms(p, xp);
@@ -899,7 +900,8 @@ void stats_write_file_header(FILE *file, const struct unit_system *restrict us,
       file,
       "#%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s "
       "%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s "
-      "%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s \n",
+      "%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s "
+      "\n",
       "(0)", "(1)", "(2)", "(3)", "(4)", "(5)", "(6)", "(7)", "(8)", "(9)",
       "(10)", "(11)", "(12)", "(13)", "(14)", "(15)", "(16)", "(17)", "(18)",
       "(19)", "(20)", "(21)", "(22)", "(23)", "(24)", "(25)", "(26)", "(27)",
@@ -909,15 +911,17 @@ void stats_write_file_header(FILE *file, const struct unit_system *restrict us,
       file,
       "#%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s "
       "%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s "
-      "%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s \n",
+      "%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s "
+      "\n",
       "Step", "Time", "a", "z", "Total mass", "Gas mass", "DM mass",
       "Sink mass", "Star mass", "BH mass", "Gas Z mass", "Star Z mass",
       "BH Z mass", "Kin. Energy", "Int. Energy", "Pot. energy", "Rad. energy",
       "Gas Entropy", "CoM x", "CoM y", "CoM z", "Mom. x", "Mom. y", "Mom. z",
       "Ang. mom. x", "Ang. mom. y", "Ang. mom. z", "BH acc. rate",
       "BH acc. mass", "BH sub. mass", "Gas H mass", "Gas H2 mass",
-      "Gas HI mass", "Gas He mass", "Mag. Energy", "DivB err", "Max divB err", "Cr. Helicity",
-      "Mag. Helicity", "RMS B field", "BH bol. lum.", "BH jet power");
+      "Gas HI mass", "Gas He mass", "Mag. Energy", "DivB err", "Max divB err",
+      "Cr. Helicity", "Mag. Helicity", "RMS B field", "BH bol. lum.",
+      "BH jet power");
 
   fflush(file);
 }
@@ -944,7 +948,7 @@ void stats_write_to_file(FILE *file, const struct statistics *stats,
       file,
       " %14d %14e %14.7f %14.7f %14e %14e %14e %14e %14e %14e %14e %14e %14e "
       "%14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e "
-      "%14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e " 
+      "%14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e "
       "%14e\n",
       step, time, a, z, stats->total_mass, stats->gas_mass, stats->dm_mass,
       stats->sink_mass, stats->star_mass, stats->bh_mass, stats->gas_Z_mass,
@@ -954,8 +958,8 @@ void stats_write_to_file(FILE *file, const struct statistics *stats,
       stats->mom[1], stats->mom[2], stats->ang_mom[0], stats->ang_mom[1],
       stats->ang_mom[2], stats->bh_accretion_rate, stats->bh_accreted_mass,
       stats->bh_subgrid_mass, stats->gas_H_mass, stats->gas_H2_mass,
-      stats->gas_HI_mass, stats->gas_He_mass, stats->E_mag, stats->divB_error, stats->divB_error_max,
-      stats->H_cross, stats->H_mag, stats->Brms,
+      stats->gas_HI_mass, stats->gas_He_mass, stats->E_mag, stats->divB_error,
+      stats->divB_error_max, stats->H_cross, stats->H_mag, stats->Brms,
       stats->bh_bolometric_luminosity, stats->bh_jet_power);
 
   fflush(file);