diff --git a/src/mhd/DirectInduction/mhd.h b/src/mhd/DirectInduction/mhd.h
index 97bbc566afdd3ee62206915bf6af501a296e7a99..66419c05b80d15d89bca96f81dc1362c1904eed6 100644
--- a/src/mhd/DirectInduction/mhd.h
+++ b/src/mhd/DirectInduction/mhd.h
@@ -292,7 +292,6 @@ __attribute__((always_inline)) INLINE static void mhd_prepare_gradient(
 
   p->mhd_data.Alfven_speed = mhd_get_comoving_Alfven_speed(p, mu_0);
 
-  p->force.balsara = 1.f;
 }
 
 /**
@@ -375,7 +374,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_comoving_magnetosonic_speed(p);
+  const float ch = 0.5f * p->viscosity.v_sig;
 
   /* Compute Dedner cleaning scalar time derivative. */
   const float hyp = hydro_props->mhd.hyp_dedner;
@@ -386,11 +385,14 @@ __attribute__((always_inline)) INLINE static float mhd_get_psi_over_ch_dt(
   const float div_v = a * a * hydro_get_div_v(p) - 3.0f * a * a * H;
   const float psi_over_ch = p->mhd_data.psi_over_ch;
 
+  const float cp = ch;
+  const float tau_inv = par * cp * h_inv; 
+
   const float hyperbolic_term =
       -hyp * a * a * a_factor_sound_speed * a_factor_sound_speed * ch * div_B;
   const float hyperbolic_divv_term = -hyp_divv * psi_over_ch * div_v;
   const float parabolic_term =
-      -par * a * a_factor_sound_speed * psi_over_ch * ch * h_inv;
+    - a * a_factor_sound_speed * psi_over_ch * tau_inv;
   const float Hubble_term = a * a * H * psi_over_ch;
 
   return hyperbolic_term + hyperbolic_divv_term + parabolic_term + Hubble_term;
diff --git a/src/mhd/DirectInduction/mhd_iact.h b/src/mhd/DirectInduction/mhd_iact.h
index 570a45db1e81c6012d5eff88995089c9ccfefcec..1200a6b488529101c121c721c422392e392752e3 100644
--- a/src/mhd/DirectInduction/mhd_iact.h
+++ b/src/mhd/DirectInduction/mhd_iact.h
@@ -326,7 +326,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
   /* B dot r. */
   const float Bri = Bi[0] * dx[0] + Bi[1] * dx[1] + Bi[2] * dx[2];
   const float Brj = Bj[0] * dx[0] + Bj[1] * dx[1] + Bj[2] * dx[2];
-
+  const float dBdr = Bri - Brj;
+  
   /* Compute gradient terms */
   const float over_rho2_i = 1.0f / (rhoi * rhoi) * f_ij;
   const float over_rho2_j = 1.0f / (rhoj * rhoj) * f_ji;
@@ -335,12 +336,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
   const float grad_term_i = over_rho2_i * wi_dr * r_inv;
   const float grad_term_j = over_rho2_j * wj_dr * r_inv;
 
-  float divB_i = Bri * grad_term_i;
-  divB_i += Brj * grad_term_j;
-  float divB_j = divB_i;
+  const float asym_grad_term_i = f_ij * wi_dr * r_inv / rhoi;
+  const float asym_grad_term_j = f_ji * wj_dr * r_inv / rhoj;
+  
+  const float divB_i = dBdr * asym_grad_term_i;
+  const float divB_j = dBdr * asym_grad_term_j;
 
-  pi->mhd_data.divB += rhoi * mj * divB_i;
-  pj->mhd_data.divB -= rhoj * mi * divB_j;
+  pi->mhd_data.divB -= mj * divB_i;
+  pj->mhd_data.divB -= mi * divB_j;
 
   /* SPH acceleration term in x direction, i_th particle */
   float sph_acc_term_i[3] = {0.f, 0.f, 0.f};
@@ -522,7 +525,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
 
   pi->u_dt -= 0.5f * mj * permeability_inv * art_diff_pref_i * dB_2;
   pj->u_dt -= 0.5f * mi * permeability_inv * art_diff_pref_j * dB_2;
-
+    
   /* Store AR terms */
   pi->mhd_data.B_over_rho_dt_AR[0] += mj * art_diff_pref_i * dB[0];
   pi->mhd_data.B_over_rho_dt_AR[1] += mj * art_diff_pref_i * dB[1];
@@ -537,21 +540,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
 
   /*Divergence diffusion */
 
-  const float vsig_Dedner_i = mhd_get_comoving_magnetosonic_speed(pi);
-  const float vsig_Dedner_j = mhd_get_comoving_magnetosonic_speed(pj);
+  const float vsig_Dedner_i = 0.5f * pi->viscosity.v_sig;
+  const float vsig_Dedner_j = 0.5f * pj->viscosity.v_sig;
 
-  const float delta_psi =
-      psi_over_ch_i * vsig_Dedner_i - psi_over_ch_j * vsig_Dedner_j;
-  const float grad_psi_i = -grad_term_i * delta_psi;
-  const float grad_psi_j = -grad_term_j * delta_psi;
+  float grad_psi = grad_term_i * psi_over_ch_i * vsig_Dedner_i;
+  grad_psi += grad_term_j * psi_over_ch_j * vsig_Dedner_j;
 
-  pi->mhd_data.B_over_rho_dt[0] -= mj * grad_psi_i * dx[0];
-  pi->mhd_data.B_over_rho_dt[1] -= mj * grad_psi_i * dx[1];
-  pi->mhd_data.B_over_rho_dt[2] -= mj * grad_psi_i * dx[2];
+  pi->mhd_data.B_over_rho_dt[0] -= mj * grad_psi * dx[0];
+  pi->mhd_data.B_over_rho_dt[1] -= mj * grad_psi * dx[1];
+  pi->mhd_data.B_over_rho_dt[2] -= mj * grad_psi * dx[2];
 
-  pj->mhd_data.B_over_rho_dt[0] -= mi * grad_psi_j * dx[0];
-  pj->mhd_data.B_over_rho_dt[1] -= mi * grad_psi_j * dx[1];
-  pj->mhd_data.B_over_rho_dt[2] -= mi * grad_psi_j * dx[2];
+  pj->mhd_data.B_over_rho_dt[0] += mi * grad_psi * dx[0];
+  pj->mhd_data.B_over_rho_dt[1] += mi * grad_psi * dx[1];
+  pj->mhd_data.B_over_rho_dt[2] += mi * grad_psi * dx[2];
 
   /* Save induction sources */
 
@@ -561,11 +562,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_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];
     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] += mj * art_diff_pref_i * dB[i];
+    pi->mhd_data.Adv_B_source[i] -= mj * grad_psi * dx[i];
+    pj->mhd_data.Adv_B_source[i] += mi * grad_psi * dx[i];
+    pi->mhd_data.Diff_B_source[i] += mj * dB_dt_pref_PR_i * wi_dr * dB[i];
+    pj->mhd_data.Diff_B_source[i] -= mi * dB_dt_pref_PR_j * wj_dr * 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];
     pj->mhd_data.Delta_B[i] -= mi * dB_dt_pref_Lap_j * wj_dr * dB[i];
@@ -653,7 +654,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
   /* B dot r. */
   const float Bri = Bi[0] * dx[0] + Bi[1] * dx[1] + Bi[2] * dx[2];
   const float Brj = Bj[0] * dx[0] + Bj[1] * dx[1] + Bj[2] * dx[2];
-
+  const float dBdr = Bri - Brj;
+  
   /* Compute gradient terms */
   const float over_rho2_i = 1.0f / (rhoi * rhoi) * f_ij;
   const float over_rho2_j = 1.0f / (rhoj * rhoj) * f_ji;
@@ -662,10 +664,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
   const float grad_term_i = over_rho2_i * wi_dr * r_inv;
   const float grad_term_j = over_rho2_j * wj_dr * r_inv;
 
-  float divB_i = Bri * grad_term_i;
-  divB_i += Brj * grad_term_j;
+  const float asym_grad_term_i = f_ij * wi_dr * r_inv / rhoi;
+  
+  const float divB_i = dBdr * asym_grad_term_i;
 
-  pi->mhd_data.divB += rhoi * mj * divB_i;
+  pi->mhd_data.divB -= mj * divB_i;
 
   /* SPH acceleration term in x direction, i_th particle */
   float sph_acc_term_i[3] = {0.f, 0.f, 0.f};
@@ -772,10 +775,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
   }
 
   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;
-
+  
   float dv_cross_dx[3];
   dv_cross_dx[0] = dv[1] * dx[2] - dv[2] * dx[1];
   dv_cross_dx[1] = dv[2] * dx[0] - dv[0] * dx[2];
@@ -785,7 +788,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
                           dv_cross_dx[1] * dv_cross_dx[1] +
                           dv_cross_dx[2] * dv_cross_dx[2];
   const float v_sig_B = sqrtf(v_sig_B_2) * r_inv;
-
+  
   const float art_diff_pref = 0.5f * art_diff_beta * v_sig_B *
                               (wi_dr * over_rho2_i + wj_dr * over_rho2_j);
 
@@ -794,7 +797,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
   pi->mhd_data.B_over_rho_dt[2] += mj * art_diff_pref * dB[2];
 
   pi->u_dt -= 0.5f * mj * permeability_inv * art_diff_pref * dB_2;
-
+    
   /* Store AR terms */
   pi->mhd_data.B_over_rho_dt_AR[0] += mj * art_diff_pref * dB[0];
   pi->mhd_data.B_over_rho_dt_AR[1] += mj * art_diff_pref * dB[1];
@@ -802,26 +805,25 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
 
   pi->mhd_data.u_dt_AR -= 0.5f * mj * permeability_inv * art_diff_pref * dB_2;
 
-  /*Divergence diffusion */
+  /* Divergence diffusion */
 
-  const float vsig_Dedner_i = mhd_get_comoving_magnetosonic_speed(pi);
-  const float vsig_Dedner_j = mhd_get_comoving_magnetosonic_speed(pj);
+  const float vsig_Dedner_i = 0.5f * pi->viscosity.v_sig;
+  const float vsig_Dedner_j = 0.5f * pj->viscosity.v_sig;
 
-  const float delta_psi =
-      psi_over_ch_i * vsig_Dedner_i - psi_over_ch_j * vsig_Dedner_j;
-  const float grad_psi_i = -grad_term_i * delta_psi;
+  float grad_psi = grad_term_i * psi_over_ch_i * vsig_Dedner_i;
+  grad_psi += grad_term_j * psi_over_ch_j * vsig_Dedner_j;
 
-  pi->mhd_data.B_over_rho_dt[0] -= mj * grad_psi_i * dx[0];
-  pi->mhd_data.B_over_rho_dt[1] -= mj * grad_psi_i * dx[1];
-  pi->mhd_data.B_over_rho_dt[2] -= mj * grad_psi_i * dx[2];
+  pi->mhd_data.B_over_rho_dt[0] -= mj * grad_psi * dx[0];
+  pi->mhd_data.B_over_rho_dt[1] -= mj * grad_psi * dx[1];
+  pi->mhd_data.B_over_rho_dt[2] -= mj * grad_psi * dx[2];
 
   /* Save induction sources */
   const float dB_dt_pref_Lap = 2.0f * r_inv / rhoj;
 
   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.Adv_B_source[i] -= mj * grad_psi * dx[i];
+    pi->mhd_data.Diff_B_source[i] += mj * dB_dt_pref_PR * wi_dr * 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];
   }