added more correct density ratio

src/mhd/DirectInduction/mhd.h

@@ -371,9 +371,9 @@ __attribute__((always_inline)) INLINE static void mhd_end_gradient(
 
   Cos_Ind_Diff = (Adv_B_source[0]*Delta_B[0]+Adv_B_source[1]*Delta_B[1]+Adv_B_source[2]*Delta_B[2]);
   
-  p->mhd_data.eta_OWAR += 1.0f/OW * ( 0.5f * p->h * p->h / (absB+FLT_MIN)) * 0.5f * fmaxf(0.0f,1.0f-Cos_Ind_Diff) * Abs_Adv_B_source * (p->mass * kernel_root);
+  p->mhd_data.eta_OWAR += 1.0f/OW * ( 0.5f * p->h * p->h / (absB+FLT_MIN)) * 0.5f * fmaxf(0.0f,1.0f-Cos_Ind_Diff) * Abs_Adv_B_source * p->rho * (p->mass * kernel_root);
 
-  p->mhd_data.eta_OWAR *= pow_dimension(1.0f/p->h);
+  p->mhd_data.eta_OWAR *= pow_dimension(1.0f/p->h) / p->rho;
 
   if (p->mhd_data.eta_OWAR<0.0f){
     error(

src/mhd/DirectInduction/mhd_iact.h

@@ -183,7 +183,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_gradient(
   }
 
   Cos_Ind_Diffi = (Adv_B_sourcei[0]*Delta_Bi[0]+Adv_B_sourcei[1]*Delta_Bi[1]+Adv_B_sourcei[2]*Delta_Bi[2]);
-  pj->mhd_data.eta_OWAR += 1.0f/OW * ( 0.5f * hi * hi / (absBi+FLT_MIN)) * 0.5f * fmaxf(0.0f,1.0f-Cos_Ind_Diffi) * Abs_Adv_B_sourcei * (wj*mi);
+  pj->mhd_data.eta_OWAR += 1.0f/OW * ( 0.5f * hi * hi / (absBi+FLT_MIN)) * 0.5f * fmaxf(0.0f,1.0f-Cos_Ind_Diffi) * Abs_Adv_B_sourcei * rhoi * (wj*mi);
   if (pj->mhd_data.eta_OWAR<0.0f){
     error(
         "Error: incorrect OWAR "
@@ -210,7 +210,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_gradient(
   }
 
   Cos_Ind_Diffj = (Adv_B_sourcej[0]*Delta_Bj[0]+Adv_B_sourcej[1]*Delta_Bj[1]+Adv_B_sourcej[2]*Delta_Bj[2]);
-  pi->mhd_data.eta_OWAR += 1.0f/OW * ( 0.5f * hj * hj / (absBj+FLT_MIN)) * 0.5f * fmaxf(0.0f,1.0f-Cos_Ind_Diffj) * Abs_Adv_B_sourcej * (wi*mj);
+  pi->mhd_data.eta_OWAR += 1.0f/OW * ( 0.5f * hj * hj / (absBj+FLT_MIN)) * 0.5f * fmaxf(0.0f,1.0f-Cos_Ind_Diffj) * Abs_Adv_B_sourcej * rhoj * (wi*mj);
   if (pi->mhd_data.eta_OWAR<0.0f){
     error(
         "Error: incorrect OWAR "
@@ -336,7 +336,7 @@ runner_iact_nonsym_mhd_gradient(const float r2, const float dx[3],
 
   Cos_Ind_Diffj = (Adv_B_sourcej[0]*Delta_Bj[0]+Adv_B_sourcej[1]*Delta_Bj[1]+Adv_B_sourcej[2]*Delta_Bj[2]);
   
-  pi->mhd_data.eta_OWAR += 1.0f/OW * ( 0.5f * hj * hj / (absBj+FLT_MIN)) * 0.5f * fmaxf(0.0f,1.0f-Cos_Ind_Diffj) * Abs_Adv_B_sourcej * (wi*mj);
+  pi->mhd_data.eta_OWAR += 1.0f/OW * ( 0.5f * hj * hj / (absBj+FLT_MIN)) * 0.5f * fmaxf(0.0f,1.0f-Cos_Ind_Diffj) * Abs_Adv_B_sourcej * rhoj * (wi*mj);
   if (pi->mhd_data.eta_OWAR<0.0f){
     error(
         "Error: incorrect OWAR "