diff --git a/src/mhd/DirectInduction/mhd.h b/src/mhd/DirectInduction/mhd.h
index c93a97d60b79a6dc928e65d45440dfd3beaca15a..97bbc566afdd3ee62206915bf6af501a296e7a99 100644
--- a/src/mhd/DirectInduction/mhd.h
+++ b/src/mhd/DirectInduction/mhd.h
@@ -143,10 +143,10 @@ __attribute__((always_inline)) INLINE static float mhd_compute_timestep(
 }
 
 /**
- * @brief Compute magnetosonic speed
+ * @brief Compute Alfven speed
  */
-__attribute__((always_inline)) INLINE static float mhd_get_magnetosonic_speed(
-    const struct part *restrict p, const float a, const float mu_0) {
+__attribute__((always_inline)) INLINE static float
+mhd_get_comoving_Alfven_speed(const struct part *restrict p, const float mu_0) {
 
   /* Recover some data */
   const float rho = p->rho;
@@ -158,24 +158,37 @@ __attribute__((always_inline)) INLINE static float mhd_get_magnetosonic_speed(
   /* B squared */
   const float B2 = B[0] * B[0] + B[1] * B[1] + B[2] * B[2];
 
-  const float permeability_inv = 1.0f / mu_0;
+  /* Square of Alfven speed */
+  const float vA2 = B2 / (mu_0 * rho);
 
-  /* Compute effective sound speeds */
-  const float cs = p->force.soundspeed;
+  return sqrtf(vA2);
+}
+
+/**
+ * @brief Compute magnetosonic speed
+ */
+__attribute__((always_inline)) INLINE static float
+mhd_get_comoving_magnetosonic_speed(const struct part *restrict p) {
+
+  /* Compute square of fast magnetosonic speed */
+  const float cs = hydro_get_comoving_soundspeed(p);
   const float cs2 = cs * cs;
-  const float v_A2 = permeability_inv * B2 / rho;
-  const float c_ms2 = cs2 + v_A2;
 
-  return sqrtf(c_ms2);
+  const float vA = p->mhd_data.Alfven_speed;
+  const float vA2 = vA * vA;
+
+  const float cms2 = cs2 + vA2;
+
+  return sqrtf(cms2);
 }
 
 /**
  * @brief Compute fast magnetosonic wave phase veolcity
  */
 __attribute__((always_inline)) INLINE static float
-mhd_get_fast_magnetosonic_wave_phase_velocity(const float dx[3],
-                                              const struct part *restrict p,
-                                              const float a, const float mu_0) {
+mhd_get_comoving_fast_magnetosonic_wave_phase_velocity(
+    const float dx[3], const struct part *restrict p, const float a,
+    const float mu_0) {
 
   /* Get r and 1/r. */
   const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
@@ -196,7 +209,7 @@ mhd_get_fast_magnetosonic_wave_phase_velocity(const float dx[3],
   /* Compute effective sound speeds */
   const float cs = p->force.soundspeed;
   const float cs2 = cs * cs;
-  const float c_ms = mhd_get_magnetosonic_speed(p, a, mu_0);
+  const float c_ms = mhd_get_comoving_magnetosonic_speed(p);
   const float c_ms2 = c_ms * c_ms;
   const float projection_correction = c_ms2 * c_ms2 - 4.0f * permeability_inv *
                                                           cs2 * Br * r_inv *
@@ -226,8 +239,8 @@ __attribute__((always_inline)) INLINE static float mhd_signal_velocity(
     const struct part *restrict pj, const float mu_ij, const float beta,
     const float a, const float mu_0) {
 
-  const float v_sigi = mhd_get_magnetosonic_speed(pi, a, mu_0);
-  const float v_sigj = mhd_get_magnetosonic_speed(pj, a, mu_0);
+  const float v_sigi = mhd_get_comoving_magnetosonic_speed(pi);
+  const float v_sigj = mhd_get_comoving_magnetosonic_speed(pj);
 
   const float v_sig = v_sigi + v_sigj - beta * mu_ij;
 
@@ -274,7 +287,10 @@ __attribute__((always_inline)) INLINE static void mhd_end_density(
  */
 __attribute__((always_inline)) INLINE static void mhd_prepare_gradient(
     struct part *restrict p, struct xpart *restrict xp,
-    const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float mu_0) {
+
+  p->mhd_data.Alfven_speed = mhd_get_comoving_Alfven_speed(p, mu_0);
 
   p->force.balsara = 1.f;
 }
@@ -359,7 +375,7 @@ __attribute__((always_inline)) INLINE static float mhd_get_psi_over_ch_dt(
   const float h_inv = 1.0f / h;
 
   /* Compute Dedner cleaning speed. */
-  const float ch = mhd_get_magnetosonic_speed(p, a, mu_0);
+  const float ch = mhd_get_comoving_magnetosonic_speed(p);
 
   /* Compute Dedner cleaning scalar time derivative. */
   const float hyp = hydro_props->mhd.hyp_dedner;
@@ -470,7 +486,8 @@ __attribute__((always_inline)) INLINE static void mhd_reset_acceleration(
  * @param cosmo The cosmological model
  */
 __attribute__((always_inline)) INLINE static void mhd_reset_predicted_values(
-    struct part *p, const struct xpart *xp, const struct cosmology *cosmo) {
+    struct part *p, const struct xpart *xp, const struct cosmology *cosmo,
+    const float mu_0) {
 
   /* Re-set the predicted magnetic flux densities */
   p->mhd_data.B_over_rho[0] = xp->mhd_data.B_over_rho_full[0];
@@ -478,6 +495,8 @@ __attribute__((always_inline)) INLINE static void mhd_reset_predicted_values(
   p->mhd_data.B_over_rho[2] = xp->mhd_data.B_over_rho_full[2];
 
   p->mhd_data.psi_over_ch = xp->mhd_data.psi_over_ch_full;
+
+  p->mhd_data.Alfven_speed = mhd_get_comoving_Alfven_speed(p, mu_0);
 }
 
 /**
@@ -507,6 +526,8 @@ __attribute__((always_inline)) INLINE static void mhd_predict_extra(
   p->mhd_data.B_over_rho[2] += p->mhd_data.B_over_rho_dt[2] * dt_therm;
 
   p->mhd_data.psi_over_ch += p->mhd_data.psi_over_ch_dt * dt_therm;
+
+  p->mhd_data.Alfven_speed = mhd_get_comoving_Alfven_speed(p, mu_0);
 }
 
 /**
@@ -596,7 +617,7 @@ __attribute__((always_inline)) INLINE static void mhd_kick_extra(
  */
 __attribute__((always_inline)) INLINE static void mhd_convert_quantities(
     struct part *p, struct xpart *xp, const struct cosmology *cosmo,
-    const struct hydro_props *hydro_props) {
+    const struct hydro_props *hydro_props, const float mu_0) {
   /* Set Restitivity Eta */
   p->mhd_data.resistive_eta = hydro_props->mhd.mhd_eta;
   /* Set Monopole subtraction factor */
@@ -614,11 +635,16 @@ __attribute__((always_inline)) INLINE static void mhd_convert_quantities(
   p->mhd_data.B_over_rho[1] *= pow(cosmo->a, 1.5f * hydro_gamma);
   p->mhd_data.B_over_rho[2] *= pow(cosmo->a, 1.5f * hydro_gamma);
 
+  /* Instantiate full step magnetic field */
   xp->mhd_data.B_over_rho_full[0] = p->mhd_data.B_over_rho[0];
   xp->mhd_data.B_over_rho_full[1] = p->mhd_data.B_over_rho[1];
   xp->mhd_data.B_over_rho_full[2] = p->mhd_data.B_over_rho[2];
 
+  /* Instantiate full step magnetic Dedner scalar */
   xp->mhd_data.psi_over_ch_full = p->mhd_data.psi_over_ch;
+
+  /* Instantiate Alfven speed */
+  p->mhd_data.Alfven_speed = mhd_get_comoving_Alfven_speed(p, mu_0);
 }
 
 /**
diff --git a/src/mhd/DirectInduction/mhd_debug.h b/src/mhd/DirectInduction/mhd_debug.h
index cf12f96129f028bfa78a3c0ab5993a86a7b7c42f..92325a4c325e5419233bca4117b92e7e1fe21ba0 100644
--- a/src/mhd/DirectInduction/mhd_debug.h
+++ b/src/mhd/DirectInduction/mhd_debug.h
@@ -33,12 +33,12 @@ __attribute__((always_inline)) INLINE static void mhd_debug_particle(
       "B/rho=[%.3e,%.3e,%.3e], "
       "B/rho (pred)=[%.3e,%.3e,%.3e], "
       "d(B/rho)/dt=[%.3e,%.3e,%.3e], \n"
-      "divB=%.3e, v_fm=%.3e, psi/cs=%.3e, d(psi/c)/dt=%.3e\n",
+      "divB=%.3e, v_Alfven=%.3e, psi/cs=%.3e, d(psi/c)/dt=%.3e\n",
       xp->mhd_data.B_over_rho_full[0], xp->mhd_data.B_over_rho_full[1],
       xp->mhd_data.B_over_rho_full[2], p->mhd_data.B_over_rho[0],
       p->mhd_data.B_over_rho[1], p->mhd_data.B_over_rho[2],
       p->mhd_data.B_over_rho_dt[0], p->mhd_data.B_over_rho_dt[1],
-      p->mhd_data.B_over_rho_dt[2], p->mhd_data.divB, p->mhd_data.v_fm,
+      p->mhd_data.B_over_rho_dt[2], p->mhd_data.divB, p->mhd_data.Alfven_speed,
       p->mhd_data.psi_over_ch, p->mhd_data.psi_over_ch_dt);
 }
 
diff --git a/src/mhd/DirectInduction/mhd_iact.h b/src/mhd/DirectInduction/mhd_iact.h
index 123d019e2976dcfc7b917a5fc0205d2d543319a6..5e94f6b9c1ca6dbfddade4624b7d6119a54313a5 100644
--- a/src/mhd/DirectInduction/mhd_iact.h
+++ b/src/mhd/DirectInduction/mhd_iact.h
@@ -534,8 +534,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
 
   /*Divergence diffusion */
 
-  const float vsig_Dedner_i = mhd_get_magnetosonic_speed(pi, a, mu_0);
-  const float vsig_Dedner_j = mhd_get_magnetosonic_speed(pj, a, mu_0);
+  const float vsig_Dedner_i = mhd_get_comoving_magnetosonic_speed(pi);
+  const float vsig_Dedner_j = mhd_get_comoving_magnetosonic_speed(pj);
 
   const float delta_psi =
       psi_over_ch_i * vsig_Dedner_i - psi_over_ch_j * vsig_Dedner_j;
@@ -796,8 +796,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
 
   /*Divergence diffusion */
 
-  const float vsig_Dedner_i = mhd_get_magnetosonic_speed(pi, a, mu_0);
-  const float vsig_Dedner_j = mhd_get_magnetosonic_speed(pj, a, mu_0);
+  const float vsig_Dedner_i = mhd_get_comoving_magnetosonic_speed(pi);
+  const float vsig_Dedner_j = mhd_get_comoving_magnetosonic_speed(pj);
 
   const float delta_psi =
       psi_over_ch_i * vsig_Dedner_i - psi_over_ch_j * vsig_Dedner_j;
diff --git a/src/mhd/DirectInduction/mhd_io.h b/src/mhd/DirectInduction/mhd_io.h
index 37b6068bbabfacfe74683843e7912d9509b0aaa9..da3c5e223301ec7b77b1b4afd20b40ad8f90c0a0 100644
--- a/src/mhd/DirectInduction/mhd_io.h
+++ b/src/mhd/DirectInduction/mhd_io.h
@@ -394,11 +394,11 @@ INLINE static int mhd_write_particles(const struct part* parts,
       -1.5f * hydro_gamma - 1.f, parts, mhd_data.psi_over_ch,
       "Dedner scalars over cleaning speeds of the particles");
 
-  list[3] = io_make_output_field(
-      "DednerScalarsOverCleaningSpeedsdt", FLOAT, 1,
-      UNIT_CONV_MAGNETIC_FIELD_PER_TIME, 1.f, parts,
-      mhd_data.psi_over_ch_dt,
-      "Time derivative of Dedner scalars over cleaning speeds of the particles");
+  list[3] = io_make_output_field("DednerScalarsOverCleaningSpeedsdt", FLOAT, 1,
+                                 UNIT_CONV_MAGNETIC_FIELD_PER_TIME, 1.f, parts,
+                                 mhd_data.psi_over_ch_dt,
+                                 "Time derivative of Dedner scalars over "
+                                 "cleaning speeds of the particles");
 
   list[4] = io_make_output_field(
       "MagneticFluxDensitiesdt", FLOAT, 3, UNIT_CONV_MAGNETIC_FIELD_PER_TIME,
diff --git a/src/mhd/DirectInduction/mhd_struct.h b/src/mhd/DirectInduction/mhd_struct.h
index 67054a31afe71b8839d24249f289967318dc86ba..ed21b62111595afae2e834dd88fd1b5c92d6c2fb 100644
--- a/src/mhd/DirectInduction/mhd_struct.h
+++ b/src/mhd/DirectInduction/mhd_struct.h
@@ -24,37 +24,50 @@
  */
 struct mhd_part_data {
 
+  /* Predicted magnetic field over density */
   float B_over_rho[3];
 
-  float divB;
-
-  /*! dB Direct Induction */
+  /* Time derivative of magnetic field over density */
   float B_over_rho_dt[3];
 
-  float v_fm;
+  /* Alfven speed (=sqrt(B2/(mu_0 * rho))) of the particle drifted to the
+   * current time */
+  float Alfven_speed;
+
+  /* Divergence of the magnetic field */
+  float divB;
 
+  /* Curl of the magnetic field */
   float curl_B[3];
 
-  /* Resistive Eta */
-  float resistive_eta;
+  /* Tensile instability correction multiplicative prefactor */
+  float monopole_beta;
 
+  /* Artifical resistivity multiplicative prefactor */
+  float art_diff_beta;
+
+  /* Spatial gradient tensor of the magnetic field */
   float grad_B_tensor[3][3];
 
+  /* Artificial resistivity gradient based switch */
   float alpha_AR;
 
-  float psi_over_ch;
-
-  float psi_over_ch_dt;
+  /* Artificial resistivity contribution to the time derivative of magnetic
+   * field over density */
+  float B_over_rho_dt_AR[3];
 
-  /*! Monopole subtraction in Lorentz Force*/
-  float monopole_beta;
+  /* Artificial resistivity contribution to the time derivative of thermal
+   * energy */
+  float u_dt_AR;
 
-  /*! Artifical Diffusion */
-  float art_diff_beta;
+  /* Physical resistive parameter */
+  float resistive_eta;
 
-  float B_over_rho_dt_AR[3];
+  /* Predicted Dedner scalar over divergence cleaning speed */
+  float psi_over_ch;
 
-  float u_dt_AR;
+  /* Time derivative of Dedner scalar over divergence cleaning speed */
+  float psi_over_ch_dt;
 
   /* SPH <1> error */
   float mean_SPH_err;
@@ -64,10 +77,13 @@ struct mhd_part_data {
 
   /* Magnetic force */
   float tot_mag_F[3];
+
   /* B advection source */
   float Adv_B_source[3];
+
   /* B total diffusion source */
   float Diff_B_source[3];
+
   /* Laplacian B */
   float Delta_B[3];
 };
@@ -77,10 +93,10 @@ struct mhd_part_data {
  */
 struct mhd_xpart_data {
 
-  /*! Full Step Magnetic Field */
+  /* Full step magnetic field over density */
   float B_over_rho_full[3];
 
-  /*! Full Step Dedner Scalar */
+  /* Full step dedner scalar over divergence cleaning speed */
   float psi_over_ch_full;
 };
 
diff --git a/src/mhd/None/mhd.h b/src/mhd/None/mhd.h
index 6370f0a9e67bcea151c376e2d49db32657fb7de9..2aa42f7ead34c8553722d40956400bc348f658d9 100644
--- a/src/mhd/None/mhd.h
+++ b/src/mhd/None/mhd.h
@@ -191,7 +191,7 @@ __attribute__((always_inline)) INLINE static void mhd_end_density(
  */
 __attribute__((always_inline)) INLINE static void mhd_prepare_gradient(
     struct part *p, struct xpart *xp, const struct cosmology *cosmo,
-    const struct hydro_props *hydro_props) {}
+    const struct hydro_props *hydro_props, const float mu_0) {}
 
 /**
  * @brief Resets the variables that are required for a gradient calculation.
@@ -271,7 +271,8 @@ __attribute__((always_inline)) INLINE static void mhd_reset_acceleration(
  * @param cosmo The cosmological model
  */
 __attribute__((always_inline)) INLINE static void mhd_reset_predicted_values(
-    struct part *p, const struct xpart *xp, const struct cosmology *cosmo) {}
+    struct part *p, const struct xpart *xp, const struct cosmology *cosmo,
+    const float mu_0) {}
 
 /**
  * @brief Predict additional particle fields forward in time when drifting
@@ -348,7 +349,7 @@ __attribute__((always_inline)) INLINE static void mhd_kick_extra(
  */
 __attribute__((always_inline)) INLINE static void mhd_convert_quantities(
     struct part *p, struct xpart *xp, const struct cosmology *cosmo,
-    const struct hydro_props *hydro_props) {}
+    const struct hydro_props *hydro_props, const float mu_0) {}
 
 /**
  * @brief Initialises the particles for the first time
diff --git a/src/mhd/VPotential/mhd.h b/src/mhd/VPotential/mhd.h
index fd37d9a8496579a0ad9cc3282d245bb99ddccfb0..ad8aca99c7f9b7a0b558f617cb4bc819936509d5 100644
--- a/src/mhd/VPotential/mhd.h
+++ b/src/mhd/VPotential/mhd.h
@@ -342,7 +342,8 @@ __attribute__((always_inline)) INLINE static void mhd_end_density(
  */
 __attribute__((always_inline)) INLINE static void mhd_prepare_gradient(
     struct part *restrict p, struct xpart *restrict xp,
-    const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float mu_0) {
 
   p->force.balsara = 1.f;
 }
@@ -496,7 +497,8 @@ __attribute__((always_inline)) INLINE static void mhd_reset_acceleration(
  * @param cosmo The cosmological model
  */
 __attribute__((always_inline)) INLINE static void mhd_reset_predicted_values(
-    struct part *p, const struct xpart *xp, const struct cosmology *cosmo) {}
+    struct part *p, const struct xpart *xp, const struct cosmology *cosmo,
+    const float mu_0) {}
 
 /**
  * @brief Predict additional particle fields forward in time when drifting
@@ -594,7 +596,7 @@ __attribute__((always_inline)) INLINE static void mhd_kick_extra(
  */
 __attribute__((always_inline)) INLINE static void mhd_convert_quantities(
     struct part *p, struct xpart *xp, const struct cosmology *cosmo,
-    const struct hydro_props *hydro_props) {
+    const struct hydro_props *hydro_props, const float mu_0) {
 
   const float a_fact = pow(cosmo->a, -mhd_comoving_factor - 1.f);
   /* Set Restitivity Eta */
diff --git a/src/runner_ghost.c b/src/runner_ghost.c
index 81748e6bcb28678ba99fc7f796937db31253bbca..2a8d1313174eeb6e222fe5d3d961cd922da3bfed 100644
--- a/src/runner_ghost.c
+++ b/src/runner_ghost.c
@@ -1110,6 +1110,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
   const struct star_formation *star_formation = e->star_formation;
   const struct hydro_props *hydro_props = e->hydro_properties;
   const struct pressure_floor_props *pressure_floor = e->pressure_floor_props;
+  const float mu_0 = e->physical_constants->const_vacuum_permeability;
 
   const int with_cosmology = (e->policy & engine_policy_cosmology);
   const int with_rt = (e->policy & engine_policy_rt);
@@ -1270,7 +1271,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
 
             /* Compute variables required for the gradient loop */
             hydro_prepare_gradient(p, xp, cosmo, hydro_props, pressure_floor);
-            mhd_prepare_gradient(p, xp, cosmo, hydro_props);
+            mhd_prepare_gradient(p, xp, cosmo, hydro_props, mu_0);
 
             /* The particle gradient values are now set.  Do _NOT_
                try to read any particle density variables! */
@@ -1458,7 +1459,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
 
         /* Compute variables required for the gradient loop */
         hydro_prepare_gradient(p, xp, cosmo, hydro_props, pressure_floor);
-        mhd_prepare_gradient(p, xp, cosmo, hydro_props);
+        mhd_prepare_gradient(p, xp, cosmo, hydro_props, mu_0);
 
         /* The particle gradient values are now set.  Do _NOT_
            try to read any particle density variables! */
diff --git a/src/runner_time_integration.c b/src/runner_time_integration.c
index 4b2417af68e1319da6e21eff8b46525d48086b7b..5834451a6a2fd67563e3ca833c80594e39f8979f 100644
--- a/src/runner_time_integration.c
+++ b/src/runner_time_integration.c
@@ -363,6 +363,7 @@ void runner_do_kick2(struct runner *r, struct cell *c, const int timer) {
   const struct hydro_props *hydro_props = e->hydro_properties;
   const struct pressure_floor_props *pressure_floor = e->pressure_floor_props;
   const struct entropy_floor_properties *entropy_floor = e->entropy_floor;
+  const float mu_0 = e->physical_constants->const_vacuum_permeability;
   const int with_cosmology = (e->policy & engine_policy_cosmology);
   const int periodic = e->s->periodic;
   const int count = c->hydro.count;
@@ -461,7 +462,7 @@ void runner_do_kick2(struct runner *r, struct cell *c, const int timer) {
 
         /* Prepare the values to be drifted */
         hydro_reset_predicted_values(p, xp, cosmo, pressure_floor);
-        mhd_reset_predicted_values(p, xp, cosmo);
+        mhd_reset_predicted_values(p, xp, cosmo, mu_0);
       }
     }
 
diff --git a/src/space.c b/src/space.c
index bcd289547e09943529714a0d9be5fc6513b41d7d..e25b20c9f6924423063def6415b4278040ac2f77 100644
--- a/src/space.c
+++ b/src/space.c
@@ -785,6 +785,7 @@ void space_convert_quantities_mapper(void *restrict map_data, int count,
   const struct cosmology *cosmo = s->e->cosmology;
   const struct hydro_props *hydro_props = s->e->hydro_properties;
   const struct pressure_floor_props *floor = s->e->pressure_floor_props;
+  const float mu_0 = s->e->physical_constants->const_vacuum_permeability;
   struct part *restrict parts = (struct part *)map_data;
   const ptrdiff_t index = parts - s->parts;
   struct xpart *restrict xparts = s->xparts + index;
@@ -795,7 +796,7 @@ void space_convert_quantities_mapper(void *restrict map_data, int count,
     if (parts[k].time_bin <= num_time_bins) {
       hydro_convert_quantities(&parts[k], &xparts[k], cosmo, hydro_props,
                                floor);
-      mhd_convert_quantities(&parts[k], &xparts[k], cosmo, hydro_props);
+      mhd_convert_quantities(&parts[k], &xparts[k], cosmo, hydro_props, mu_0);
     }
   }
 }