revert to Fast Rotor

examples/MHDTests/FastRotor/FastRotor_schemes.yml

@@ -9,7 +9,7 @@ InternalUnitSystem:
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   0.3   # The end time of the simulation (in internal units).
+  time_end:   0.1   # The end time of the simulation (in internal units).
   dt_min:     1e-9  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-3  # The maximal time-step size of the simulation (in internal units).
 
@@ -26,17 +26,20 @@ Statistics:
 
 # Parameters for the hydrodynamics scheme
 SPH:
-  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel$
+  resolution_eta:        1.7 #1.595 #1.487    # 1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel$
   CFL_condition:         0.075    # Courant-Friedrich-Levy condition for time integration.
+#  viscosity_alpha: 10.0
+#  viscosity_alpha_max: 20
+
 
 # Parameters for MHD schemes
 MHD:
   monopole_subtraction:   1.0
-  artificial_diffusion:   1.0
+  artificial_diffusion:   0.0
   hyperbolic_dedner:      1.0
   hyperbolic_dedner_divv: 0.5
   parabolic_dedner:       1.0
-  resistive_eta:          0.001
+  resistive_eta:          0.01
 
 # Physical parameters
 PhysicalConstants:

examples/MHDTests/FastRotor/makeIC_schemes.py

@@ -27,7 +27,7 @@ unit_cell = glass["/PartType0/Coordinates"][:, :]
 h_unit_cell = glass["/PartType0/SmoothingLength"][:]
 
 N_unit_cell = len(h_unit_cell)
-times = 2
+times = 4
 
 ratio = np.cbrt(rho_in_0 / rho_out_0)
 
@@ -36,6 +36,7 @@ ratio = np.cbrt(rho_in_0 / rho_out_0)
 cx_out = times
 cy_out = times
 cz_out = 1
+cz_out = times
 
 cx_in = int(np.ceil(ratio * cx_out))
 cy_in = int(np.ceil(ratio * cy_out))
@@ -144,6 +145,7 @@ u[N_out_f:] = P_0 / (rho_in_0 * (gamma - 1))
 
 v[N_out_f:, 0] = -rot[N_out_f:] * omega_0 * np.sin(theta[N_out_f:])
 v[N_out_f:, 1] = rot[N_out_f:] * omega_0 * np.cos(theta[N_out_f:])
+v[N_out_f:, 2] = 0.0 
 
 # B[:, 0] = B_0
 
@@ -164,7 +166,7 @@ hh[N:] = h
 h = hh
 vel = np.zeros((N2, 3))
 vel[:N, :] = v[:, :]
-vel[N:, :] = v[:, :]
+vel[N:, :] = -v[:, :]
 v = vel
 ids = np.linspace(1, N2, N2)
 m = np.ones(N2) * rho_out_0 * vol / N_out

examples/MHDTests/FastRotor/plot_schemes.py

@@ -78,12 +78,14 @@ for ii in range(nini, nfin):
 
     B = data.gas.magnetic_flux_densities
     divB = data.gas.magnetic_divergences
+    #divB = data.gas.vector_potential_divergences
+    #divB = data.gas.velocities[:,2]
     P_mag = (B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2) / 2
     h = data.gas.smoothing_lengths
 
     normB = np.sqrt(B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2)
 
-    data.gas.DivB_error = np.maximum(h * abs(divB) / normB, 1e-10)
+    data.gas.DivB_error =  np.maximum(h * abs(divB) / normB, 1e-10)
 
     # Then create a mass-weighted B error dataset
     data.gas.mass_weighted_magnetic_divB_error = data.gas.masses * data.gas.DivB_error

examples/MHDTests/FastRotor/run_schemes.sh

@@ -53,7 +53,8 @@ case $WHAT in
 	  SCHEME_Nr=4
 	;;
 	all)
-	  SCHEME_Nr=( 0 1 2 3 4 )
+#	  SCHEME_Nr=( 0 1 2 3 4 )
+	  SCHEME_Nr=( 0 1 2 )
 	;;
 	*)
 	  echo $WHAT" wrong scheme"
@@ -80,14 +81,14 @@ do
       cur_dir=`pwd`
       cd ../../../../
       pwd
-      ./TestAllMHD.sh $IDD "--with-adiabatic-index=7/5"
+      ./TestAllMHD.sh $IDD "--with-adiabatic-index=7/5 --with-kernel=wendland-C2 --disable-hand-vec"
       cd $cur_dir
       cp ../../../../sw_$ID .
    fi
    cat <<-EOF > ./run.sh
 	#!/bin/bash
 	# Run SWIFT
-	./sw_$ID --hydro --threads=4 ../FastRotor_schemes.yml 2>&1 > out.log 
+	./sw_$ID --hydro --threads=12 ../FastRotor_schemes.yml 2>&1 > out.log 
 	
 	# Plot the evolution
 	python3 ../plot_schemes.py 0 61 2>&1 > plot.log

src/mhd/VPotential/mhd.h

@@ -127,22 +127,15 @@ __attribute__((always_inline)) INLINE static float mhd_compute_timestep(
     const struct hydro_props *hydro_properties, const struct cosmology *cosmo,
     const float mu_0) {
 
-  //const float afac_divB = pow(cosmo->a, -mhd_comoving_factor - 0.5f); //(cosmo->a / cosmo->a_factor_sound_speed);
   const float afac_divB = (cosmo->a / cosmo->a_factor_sound_speed);
   const float afac_resistive = cosmo->a * cosmo->a;
   /* Dt from 1/DivOperator(Alfven speed) */
 
-  float dt_divA =
-      p->mhd_data.divA != 0.0f
-          ? afac_divB * hydro_properties->CFL_condition * p->h *
-                sqrtf(p->rho / (p->mhd_data.divA * p->mhd_data.divA) * mu_0)
-          : FLT_MAX;
   float dt_divB =
       p->mhd_data.divB != 0.0f
           ? afac_divB * hydro_properties->CFL_condition *
                 sqrtf(p->rho / (p->mhd_data.divB * p->mhd_data.divB) * mu_0)
           : FLT_MAX;
-  dt_divB = fminf(dt_divA, dt_divB);
   const float resistive_eta = p->mhd_data.resistive_eta;
   const float dt_eta = resistive_eta != 0.0f
                            ? afac_resistive * hydro_properties->CFL_condition *
@@ -244,9 +237,46 @@ __attribute__((always_inline)) INLINE static float mhd_signal_velocity(
       mhd_get_fast_magnetosonic_wave_phase_velocity(dx, pj, a, mu_0);
   */
 
-  const float v_sig = v_sigi + v_sigj - beta * mu_ij;
+  const float v_sig = v_sigi + v_sigj - beta/2.f * mu_ij;
 
   return v_sig;
+/*  const float a_fac =
+      2.f * mhd_comoving_factor + 3.f + 3.f * (hydro_gamma - 1.f);
+  const float ci = pi->force.soundspeed;
+  const float cj = pj->force.soundspeed;
+  const float r2 = (dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
+  const float r_inv = r2 ? 1.f / sqrt(r2) : 0.0f;
+
+  const float b2_i = (pi->mhd_data.BPred[0] * pi->mhd_data.BPred[0] +
+                      pi->mhd_data.BPred[1] * pi->mhd_data.BPred[1] +
+                      pi->mhd_data.BPred[2] * pi->mhd_data.BPred[2]);
+  const float b2_j = (pj->mhd_data.BPred[0] * pj->mhd_data.BPred[0] +
+                      pj->mhd_data.BPred[1] * pj->mhd_data.BPred[1] +
+                      pj->mhd_data.BPred[2] * pj->mhd_data.BPred[2]);
+  const float vcsa2_i = ci * ci + pow(a, a_fac) * b2_i / pi->rho * 0.5 / mu_0;
+  const float vcsa2_j = cj * cj + pow(a, a_fac) * b2_j / pj->rho * 0.5 / mu_0;
+  float Bpro2_i =
+      (pi->mhd_data.BPred[0] * dx[0] + pi->mhd_data.BPred[1] * dx[1] +
+       pi->mhd_data.BPred[2] * dx[2]) *
+      r_inv;
+  Bpro2_i *= Bpro2_i;
+  float mag_speed_i = sqrtf(
+      0.5 * (vcsa2_i +
+             sqrtf(max((vcsa2_i * vcsa2_i - 4.f * ci * ci * pow(a, a_fac) *
+                                                Bpro2_i / pi->rho * 0.5 / mu_0),
+                       0.f))));
+  float Bpro2_j =
+      (pj->mhd_data.BPred[0] * dx[0] + pj->mhd_data.BPred[1] * dx[1] +
+       pj->mhd_data.BPred[2] * dx[2]) *
+      r_inv;
+  Bpro2_j *= Bpro2_j;
+  float mag_speed_j = sqrtf(
+      0.5 * (vcsa2_j +
+             sqrtf(max((vcsa2_j * vcsa2_j - 4.f * cj * cj * pow(a, a_fac) *
+                                                Bpro2_j / pj->rho * 0.5 / mu_0),
+                       0.f))));
+
+  return (mag_speed_i + mag_speed_j - beta  * mu_ij);*/
 }
 
 /**
@@ -259,10 +289,7 @@ __attribute__((always_inline)) INLINE static float mhd_signal_velocity(
  * @param p The particle to act upon
  */
 __attribute__((always_inline)) INLINE static void mhd_init_part(
-    struct part *p) 
- 
-    {
-
+    struct part *p) {
       p->mhd_data.divA = 0.f;
       // XXX todo, really is not the predicted variable will be the full step
       p->mhd_data.BPred[0] = 0.f;
@@ -288,7 +315,7 @@ __attribute__((always_inline)) INLINE static void mhd_end_density(
   const float h_inv_dim_plus_one =
       pow_dimension_plus_one(1.f / p->h); /*1/h^(d+1) */
   const float rho_inv = 1.f / p->rho;
-  //p->mhd_data.divA *= h_inv_dim_plus_one * rho_inv;
+  p->mhd_data.divA *= h_inv_dim_plus_one * rho_inv;
   for (int i = 0; i < 3; i++)
     p->mhd_data.BPred[i] *= h_inv_dim_plus_one * rho_inv;
 }
@@ -321,7 +348,6 @@ __attribute__((always_inline)) INLINE static void mhd_prepare_gradient(
  */
 __attribute__((always_inline)) INLINE static void mhd_reset_gradient(
     struct part *p) {
-  p->mhd_data.divA = 0.f;
   p->mhd_data.divB = 0.f;
 
   p->mhd_data.BSmooth[0] = 0.f;
@@ -395,21 +421,23 @@ __attribute__((always_inline)) INLINE static void mhd_part_has_no_neighbours(
  */
 __attribute__((always_inline)) INLINE static float hydro_get_dGau_dt(
     const struct part *restrict p, const float Gauge, const struct cosmology *c)
-    //,  const float mu_0) 
     {
 
   //const float v_sig = mhd_get_magnetosonic_speed(p, c->a, mu_0);
   //const float v_sig = p->viscosity.v_sig;
   const float v_sig = hydro_get_signal_velocity(p);
+  //const float v_sig = p->mhd_data.MaxVsig;
   const float afac1 = pow(c->a_factor_sound_speed, 2.f) * c->a * c->a;
   const float afac2 = c->a_factor_sound_speed * c->a;
   // Everything is with a2 for the time integration
   /* Hyperbolic term */
-  const float Source_Term = 1.f * afac1 * p->mhd_data.divA * (v_sig * v_sig);
+  const float Source_Term = 0.1 * afac1 * p->mhd_data.divA * (v_sig * v_sig);
+  //const float Source_Term = 0.1 * afac1 * p->mhd_data.divB * p->h * (v_sig * v_sig);
+  //const float Source_Term = 0.1 * afac1 * p->mhd_data.divA * p->h * v_sig;
   /* Parabolic evolution term */
-  const float Damping_Term = 4.f * afac2 * v_sig * Gauge / p->h;
+  const float Damping_Term = 2.f * afac2 * v_sig * Gauge / p->h;
   /* Density change term */
-  const float DivV_Term =  hydro_get_div_v(p) * Gauge * c->a * c->a ;
+  const float DivV_Term = 0.f * hydro_get_div_v(p) * Gauge * c->a * c->a ;
   /* Cosmological term */
   const float Hubble_Term = (2.f + mhd_comoving_factor + 3.f) * c->H * c->a * c->a * Gauge;
 
@@ -438,33 +466,25 @@ __attribute__((always_inline)) INLINE static void mhd_prepare_force(
     const struct hydro_props *hydro_props, const float dt_alpha,
     const float mu_0) {
 
-//  const float mu_0_1 = 1.f / mu_0;
-//  const float pressure = hydro_get_comoving_pressure(p);
-//  const float b2 = (p->mhd_data.BPred[0] * p->mhd_data.BPred[0] +
-//                    p->mhd_data.BPred[1] * p->mhd_data.BPred[1] +
-//                    p->mhd_data.BPred[2] * p->mhd_data.BPred[2]);
+  const float mu_0_1 = 1.f / mu_0;
+  const float pressure = hydro_get_comoving_pressure(p);
+  const float b2 = (p->mhd_data.BPred[0] * p->mhd_data.BPred[0] +
+                    p->mhd_data.BPred[1] * p->mhd_data.BPred[1] +
+                    p->mhd_data.BPred[2] * p->mhd_data.BPred[2]);
   /* Estimation of the tensile instability due divB */
-//  p->mhd_data.Q0 = pressure / (b2 / 2.0f * mu_0_1);  // Plasma Beta
-//  p->mhd_data.Q0 =
-//      p->mhd_data.Q0 < 10.0f ? 1.0f : 0.0f;  // No correction if not magnetized
+  p->mhd_data.Q0 = pressure / (b2 / 2.0f * mu_0_1);  // Plasma Beta
+  p->mhd_data.Q0 =
+      p->mhd_data.Q0 < 10.0f ? 1.0f : 0.0f;  // No correction if not magnetized
   /* divB contribution */
-//  const float ACC_corr = fabs(p->mhd_data.divB * sqrt(b2) * mu_0_1);
+  const float ACC_corr = fabs(p->mhd_data.divB * sqrt(b2) * mu_0_1);
   // this should go with a /p->h, but I
   // take simplify becasue of ACC_mhd also.
   /* isotropic magnetic presure */
   // add the correct hydro acceleration?
-//  const float ACC_mhd = (b2 / p->h) * mu_0_1;
+  const float ACC_mhd = (b2 / p->h) * mu_0_1;
   /* Re normalize the correction in the momentum from the DivB errors*/
-//  p->mhd_data.Q0 =
-//      ACC_corr > ACC_mhd ? p->mhd_data.Q0 * ACC_mhd / ACC_corr : p->mhd_data.Q0;
-  //float delim = 1.f;
-  //float change_Gau =
-  //    hydro_get_dGau_dt(p, xp->mhd_data.Gau_full, cosmo) * dt_alpha;
-  //delim = fabs(change_Gau/xp->mhd_data.Gau_full) < 0.25f ?  1.f : 0.25f * fabs(xp->mhd_data.Gau_full / change_Gau);
-  
-  //xp->mhd_data.Gau_full += change_Gau * delim;
-  
-  //p->mhd_data.Gau = xp->mhd_data.Gau_full;
+  p->mhd_data.Q0 =
+      ACC_corr > ACC_mhd ? p->mhd_data.Q0 * ACC_mhd / ACC_corr : p->mhd_data.Q0;
 }
 
 /**
@@ -509,9 +529,6 @@ __attribute__((always_inline)) INLINE static void mhd_reset_predicted_values(
   p->mhd_data.APred[1] = xp->mhd_data.Afull[1];
   p->mhd_data.APred[2] = xp->mhd_data.Afull[2];
   p->mhd_data.Gau = xp->mhd_data.Gau_full;
-  ///xp->mhd_data.Afull[0] = p->mhd_data.APred[0] ;
-  ///xp->mhd_data.Afull[1] = p->mhd_data.APred[1] ;
-  ///xp->mhd_data.Afull[2] = p->mhd_data.APred[2] ;
 }
 
 /**
@@ -539,14 +556,10 @@ __attribute__((always_inline)) INLINE static void mhd_predict_extra(
   p->mhd_data.APred[0] += p->mhd_data.dAdt[0] * dt_therm;
   p->mhd_data.APred[1] += p->mhd_data.dAdt[1] * dt_therm;
   p->mhd_data.APred[2] += p->mhd_data.dAdt[2] * dt_therm;
-  ///xp->mhd_data.Afull[0] += p->mhd_data.dAdt[0] * dt_therm;
-  ///xp->mhd_data.Afull[1] += p->mhd_data.dAdt[1] * dt_therm;
-  ///xp->mhd_data.Afull[2] += p->mhd_data.dAdt[2] * dt_therm;
-  float delim = 1.f;
+
   float change_Gau =
-      hydro_get_dGau_dt(p, xp->mhd_data.Gau_full, cosmo) * dt_therm;
-  //float delim = fabs(change_Gau/p->mhd_data.Gau) < 0.25f ?  1.f : 0.25f * fabs(p->mhd_data.Gau / change_Gau);
-  p->mhd_data.Gau += change_Gau * delim;
+      hydro_get_dGau_dt(p, p->mhd_data.Gau, cosmo) * dt_therm;
+  p->mhd_data.Gau += change_Gau;
 }
 
 /**
@@ -593,29 +606,11 @@ __attribute__((always_inline)) INLINE static void mhd_kick_extra(
     const struct cosmology *cosmo, const struct hydro_props *hydro_props,
     const struct entropy_floor_properties *floor_props) {
   // Slope delimiter
-  /*const float dAdt_2 = 
-	  p->mhd_data.dAdt[0] * p->mhd_data.dAdt[0] +
-	  p->mhd_data.dAdt[1] * p->mhd_data.dAdt[1] +
-	  p->mhd_data.dAdt[2] * p->mhd_data.dAdt[2] ;
-  const float Afull_2 = 
-	  xp->mhd_data.Afull[0] * xp->mhd_data.Afull[0] +
-	  xp->mhd_data.Afull[1] * xp->mhd_data.Afull[1] +
-	  xp->mhd_data.Afull[2] * xp->mhd_data.Afull[2];
-  float delim = dAdt_2 * dt_therm * dt_therm / Afull_2 < 1.0f  ? 1.f : sqrt(Afull_2 / dAdt_2) / dt_therm ;
-  */
-  float delim = 1.f;
-  xp->mhd_data.Afull[0] += p->mhd_data.dAdt[0] * dt_therm * delim;
-  xp->mhd_data.Afull[1] += p->mhd_data.dAdt[1] * dt_therm * delim;
-  xp->mhd_data.Afull[2] += p->mhd_data.dAdt[2] * dt_therm * delim;
-  ///p->mhd_data.APred[0] += p->mhd_data.dAdt[0] * dt_therm * delim;
-  ///p->mhd_data.APred[1] += p->mhd_data.dAdt[1] * dt_therm * delim;
-  ///p->mhd_data.APred[2] += p->mhd_data.dAdt[2] * dt_therm * delim;
-
-  float change_Gau =
-      hydro_get_dGau_dt(p, xp->mhd_data.Gau_full, cosmo) * dt_therm;
-  //delim = fabs(change_Gau/xp->mhd_data.Gau_full) < 0.25f ?  1.f : 0.25f * fabs(xp->mhd_data.Gau_full / change_Gau);
-  xp->mhd_data.Gau_full += change_Gau * delim;
-   
+  xp->mhd_data.Afull[0] += p->mhd_data.dAdt[0] * dt_therm;
+  xp->mhd_data.Afull[1] += p->mhd_data.dAdt[1] * dt_therm;
+  xp->mhd_data.Afull[2] += p->mhd_data.dAdt[2] * dt_therm;
+  float change_Gau = hydro_get_dGau_dt(p, p->mhd_data.Gau, cosmo) * dt_therm;
+  xp->mhd_data.Gau_full += change_Gau;
 }
 
 /**
@@ -648,9 +643,9 @@ __attribute__((always_inline)) INLINE static void mhd_convert_quantities(
   xp->mhd_data.Afull[1] = p->mhd_data.APred[1];
   xp->mhd_data.Afull[2] = p->mhd_data.APred[2];
   
-  xp->mhd_data.Gau_full = 0.f;
-  p->mhd_data.Gau = 0.f;
-
+  p->mhd_data.dAdt[0] = 0.0;
+  p->mhd_data.dAdt[1] = 0.0;
+  p->mhd_data.dAdt[2] = 0.0;
 }
 
 /**
@@ -666,6 +661,12 @@ __attribute__((always_inline)) INLINE static void mhd_convert_quantities(
 __attribute__((always_inline)) INLINE static void mhd_first_init_part(
     struct part *restrict p, struct xpart *restrict xp,
     const struct mhd_global_data *mhd_data, const double Lsize) {
+  xp->mhd_data.Gau_full = 0.0f;
+  p->mhd_data.divB = 0.0f;
+  
+  p->mhd_data.dAdt[0] = 0.0;
+  p->mhd_data.dAdt[1] = 0.0;
+  p->mhd_data.dAdt[2] = 0.0;
 
   mhd_reset_acceleration(p);
   mhd_init_part(p);

src/mhd/VPotential/mhd_iact.h

@@ -61,22 +61,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_density(
    for (int i = 0; i < 3; ++i)
       dA[i] = pi->mhd_data.APred[i] - pj->mhd_data.APred[i];
    
-   //const double dAdr = dA[0] * dx[0] + dA[1] * dx[1] + dA[2] * dx[2];
-
-   //pi->mhd_data.divA -= faci * dAdr;
-   //pj->mhd_data.divA -= facj * dAdr;
-   /////
-   // bi = dj ak - dk aj
-   // bj = dk ai - di ak
-   // bk = di aj - dj ai
-   //
-   for (int i = 0; i < 3; ++i) {
-      pi->mhd_data.BPred[i] += faci * (dA[(i + 1) % 3] * dx[(i + 2) % 3] -
+   const double dAdr = dA[0] * dx[0] + dA[1] * dx[1] + dA[2] * dx[2];
+  pi->mhd_data.divA -= faci * dAdr;
+  pj->mhd_data.divA -= facj * dAdr;
+  /////
+  // bi = dj ak - dk aj
+  // bj = dk ai - di ak
+  // bk = di aj - dj ai
+  //
+  for (int i = 0; i < 3; ++i) {
+    pi->mhd_data.BPred[i] += faci * (dA[(i + 1) % 3] * dx[(i + 2) % 3] -
                                dA[(i + 2) % 3] * dx[(i + 1) % 3]);
-      pj->mhd_data.BPred[i] += facj * (dA[(i + 1) % 3] * dx[(i + 2) % 3] -
+    pj->mhd_data.BPred[i] += facj * (dA[(i + 1) % 3] * dx[(i + 2) % 3] -
                                dA[(i + 2) % 3] * dx[(i + 1) % 3]);
-   }
-   return;
+  }
 }
 
 
@@ -113,13 +111,12 @@ runner_iact_nonsym_mhd_density(const float r2, const float dx[3],
    double dA[3];
    for (int i = 0; i < 3; ++i)
       dA[i] = pi->mhd_data.APred[i] - pj->mhd_data.APred[i];
-   //const double dAdr = dA[0] * dx[0] + dA[1] * dx[1] + dA[2] * dx[2];
-   //pi->mhd_data.divA -= faci * dAdr;
+   const double dAdr = dA[0] * dx[0] + dA[1] * dx[1] + dA[2] * dx[2];
+   pi->mhd_data.divA -= faci * dAdr;
   
    for (int i = 0; i < 3; ++i)
       pi->mhd_data.BPred[i] += faci * (dA[(i + 1) % 3] * dx[(i + 2) % 3] -
                              dA[(i + 2) % 3] * dx[(i + 1) % 3]);
-   return;
 }
 
 /**
@@ -158,12 +155,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_gradient(
   const float rhoj = pj->rho;
 
   float Bi[3], Bj[3];
-  float Ai[3], Aj[3];
   for (int i = 0; i < 3; ++i) {
     Bi[i] = pi->mhd_data.BPred[i];
     Bj[i] = pj->mhd_data.BPred[i];
-    Ai[i] = pi->mhd_data.APred[i];
-    Aj[i] = pj->mhd_data.APred[i];
   }
 
   /* Get the kernel for hi. */
@@ -187,9 +181,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_gradient(
   /* 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];
-  /* A dot r. */
-  const float Ari = Ai[0] * dx[0] + Ai[1] * dx[1] + Ai[2] * dx[2];
-  const float Arj = Aj[0] * dx[0] + Aj[1] * dx[1] + Aj[2] * dx[2];
   /* Compute gradient terms */
   const float over_rho_i = 1.0f / rhoi * f_ij;
   const float over_rho_j = 1.0f / rhoj * f_ji;
@@ -199,11 +190,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_gradient(
   float B_mon_j = -over_rho_j * (Bri - Brj) * wj_dr * r_inv;
   pi->mhd_data.divB += mj * B_mon_i;
   pj->mhd_data.divB += mi * B_mon_j;
-  /* Calculate divergence term */
-  float A_mon_i = -over_rho_i * (Ari - Arj) * wi_dr * r_inv;
-  float A_mon_j = -over_rho_j * (Ari - Arj) * wj_dr * r_inv;
-  pi->mhd_data.divA += mj * A_mon_i;
-  pj->mhd_data.divA += mi * A_mon_j;
 
   /* dB */
   float dB[3];
@@ -227,8 +213,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_gradient(
   for (int k = 0; k < 3; k++) {
     pi->mhd_data.curlB[k] += mj * over_rho_i * wi_dr * r_inv * dB_cross_dx[k];
     pj->mhd_data.curlB[k] += mi * over_rho_j * wj_dr * r_inv * dB_cross_dx[k];
-    //pi->mhd_data.BSmooth[k] += mj * over_rho_i * wi_dr * r_inv * dA_cross_dx[k];
-    //pj->mhd_data.BSmooth[k] += mi * over_rho_j * wj_dr * r_inv * dA_cross_dx[k];
     pi->mhd_data.BSmooth[k] += mj * wi * Bi[k];
     pj->mhd_data.BSmooth[k] += mi * wj * Bj[k];
   }
@@ -283,12 +267,9 @@ runner_iact_nonsym_mhd_gradient(const float r2, const float dx[3],
   const float rhoi = pi->rho;
 
   float Bi[3], Bj[3];
-  float Ai[3], Aj[3];
   for (int i = 0; i < 3; ++i) {
     Bi[i] = pi->mhd_data.BPred[i];
     Bj[i] = pj->mhd_data.BPred[i];
-    Ai[i] = pi->mhd_data.APred[i];
-    Aj[i] = pj->mhd_data.APred[i];
   }
 
   /* Get the kernel for hi. */
@@ -304,18 +285,13 @@ runner_iact_nonsym_mhd_gradient(const float r2, const float dx[3],
   /* 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]);
-  /* A dot r. */
-  const float Ari = Ai[0] * dx[0] + Ai[1] * dx[1] + Ai[2] * dx[2];
-  const float Arj = Aj[0] * dx[0] + Aj[1] * dx[1] + Aj[2] * dx[2];
+
   /* Compute gradient terms */
   const float over_rho_i = 1.0f / rhoi * f_ij;
 
   /* Calculate divergence term */
   float B_mon_i = -over_rho_i * (Bri - Brj) * wi_dr * r_inv;
   pi->mhd_data.divB += mj * B_mon_i;
-  /* Calculate divergence term */
-  float A_mon_i = -over_rho_i * (Ari - Arj) * wi_dr * r_inv;
-  pi->mhd_data.divA += mj * A_mon_i;
 
   /* dB */
   float dB[3];
@@ -338,7 +314,6 @@ runner_iact_nonsym_mhd_gradient(const float r2, const float dx[3],
   /* Calculate Curl */
   for (int k = 0; k < 3; k++) {
     pi->mhd_data.curlB[k] += mj * over_rho_i * wi_dr * r_inv * dB_cross_dx[k];
-    //pi->mhd_data.BSmooth[k] += mj * over_rho_i * wi_dr * r_inv * dA_cross_dx[k];
     pi->mhd_data.BSmooth[k] += mj * wi * Bi[k];
   }
   pi->mhd_data.Q0 += pj->mass * wi;
@@ -376,8 +351,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
   /* Recover some data */
   const float mi = pi->mass;
   const float mj = pj->mass;
-  const float Pi = pi->force.pressure;
-  const float Pj = pj->force.pressure;
+//  const float Pi = pi->force.pressure;
+//  const float Pj = pj->force.pressure;
   const float rhoi = pi->rho;
   const float rhoj = pj->rho;
 
@@ -414,8 +389,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
     Bi[i] = pi->mhd_data.BPred[i];
     Bj[i] = pj->mhd_data.BPred[i];
   }
-  const float B2i = Bi[0] * Bi[0] + Bi[1] * Bi[1] + Bi[2] * Bi[2];
-  const float B2j = Bj[0] * Bj[0] + Bj[1] * Bj[1] + Bj[2] * Bj[2];
+//  const float B2i = Bi[0] * Bi[0] + Bi[1] * Bi[1] + Bi[2] * Bi[2];
+//  const float B2j = Bj[0] * Bj[0] + Bj[1] * Bj[1] + Bj[2] * Bj[2];
 
 
   ///////////////////////////// FORCE MAXWELL TENSOR
@@ -428,14 +403,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
     mm_i[j][j] -= 0.5 * (Bi[0] * Bi[0] + Bi[1] * Bi[1] + Bi[2] * Bi[2]);
     mm_j[j][j] -= 0.5 * (Bj[0] * Bj[0] + Bj[1] * Bj[1] + Bj[2] * Bj[2]);
   }
-  const float plasma_beta_i = B2i != 0.0f ? 2.0f * mu_0 * Pi / B2i : FLT_MAX;
+/*  const float plasma_beta_i = B2i != 0.0f ? 2.0f * mu_0 * Pi / B2i : FLT_MAX;
   const float plasma_beta_j = B2j != 0.0f ? 2.0f * mu_0 * Pj / B2j : FLT_MAX;
 
   const float scale_i = 0.125f * (10.0f - plasma_beta_i);
   const float scale_j = 0.125f * (10.0f - plasma_beta_j);
 
   pi->mhd_data.Q0 = fmaxf(0.0f, fminf(scale_i, 1.0f));
-  pj->mhd_data.Q0 = fmaxf(0.0f, fminf(scale_j, 1.0f));
+  pj->mhd_data.Q0 = fmaxf(0.0f, fminf(scale_j, 1.0f));*/
   //////////////////////////// Apply to the Force and DIVB TERM SUBTRACTION
   for (int i = 0; i < 3; i++)
     for (int j = 0; j < 3; j++) {
@@ -477,23 +452,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
   dv[1] = pi->v[1] - pj->v[1];
   dv[2] = pi->v[2] - pj->v[2];
  
-  const float sourceAi = dv[0] * pi->mhd_data.APred[0] +
-                         dv[1] * pi->mhd_data.APred[1] +
-                         dv[2] * pi->mhd_data.APred[2];
-  const float sourceAj = dv[0] * pj->mhd_data.APred[0] +
-                         dv[1] * pj->mhd_data.APred[1] +
-                         dv[2] * pj->mhd_data.APred[2];
-
-//  float SAi = sourceAi + a * a * (pi->mhd_data.Gau - pj->mhd_data.Gau);
-//  float SAj = sourceAj + a * a * (pi->mhd_data.Gau - pj->mhd_data.Gau);
+  const float sAi = dv[0] * pi->mhd_data.APred[0] + dv[1] * pi->mhd_data.APred[1] + dv[2] * pi->mhd_data.APred[2];
+  const float sAj = dv[0] * pj->mhd_data.APred[0] + dv[1] * pj->mhd_data.APred[1] + dv[2] * pj->mhd_data.APred[2];
+//  const float sAi = -(dA[0] * pi->v[0] + dA[1] * pi->v[1] + dA[2] * pi->v[2]);
+//  const float sAj = -(dA[0] * pj->v[0] + dA[1] * pj->v[1] + dA[2] * pj->v[2]);
+  const float sourceAi = sAi + a * a * (pi->mhd_data.Gau - pj->mhd_data.Gau);
+  const float sourceAj = sAj + a * a * (pi->mhd_data.Gau - pj->mhd_data.Gau);
 
   for (int i = 0; i < 3; i++) {
-//    pi->mhd_data.dAdt[i] += mj * mag_VPIndi * SAi * dx[i];
-//    pj->mhd_data.dAdt[i] += mi * mag_VPIndj * SAj * dx[i];
     pi->mhd_data.dAdt[i] += mj * mag_VPIndi * sourceAi * dx[i];
     pj->mhd_data.dAdt[i] += mi * mag_VPIndj * sourceAj * dx[i];
-    pi->mhd_data.dAdt[i] += mj * a * a * (mag_VPIndi * pi->mhd_data.Gau + mag_VPIndj * pj->mhd_data.Gau) * dx[i];
-    pj->mhd_data.dAdt[i] += mi * a * a * (mag_VPIndi * pi->mhd_data.Gau + mag_VPIndj * pj->mhd_data.Gau) * dx[i];
+//    pi->mhd_data.dAdt[i] += mj * dA[i] * mag_VPIndi * (dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2]);
+//    pj->mhd_data.dAdt[i] += mi * dA[i] * mag_VPIndj * (dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2]);
   }
   /// DISSSIPATION
   const float mag_Disi =
@@ -511,8 +481,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
   for (int i = 0; i < 3; i++) {
     pi->mhd_data.Adv_A_source[i] += mj * mag_VPIndi * sourceAi * dx[i];
     pj->mhd_data.Adv_A_source[i] += mi * mag_VPIndj * sourceAj * dx[i];
-    pi->mhd_data.Adv_A_source[i] += mj * a * a * (mag_VPIndi * pi->mhd_data.Gau + mag_VPIndj * pj->mhd_data.Gau) * dx[i];
-    pj->mhd_data.Adv_A_source[i] += mi * a * a * (mag_VPIndi * pi->mhd_data.Gau + mag_VPIndj * pj->mhd_data.Gau) * dx[i];
     pi->mhd_data.Diff_A_source[i] += mj * 8.0 * pi->mhd_data.resistive_eta * mag_Disi * dA[i];
     pj->mhd_data.Diff_A_source[i] -= mi * 8.0 * pj->mhd_data.resistive_eta * mag_Disj * dA[i];
     pi->mhd_data.Delta_A[i] += mj * 8.0 * mag_Disi * dA[i];
@@ -549,7 +517,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
   const float mi = pi->mass;
   const float mj = pj->mass;
 
-  const float Pi = pi->force.pressure;
+// const float Pi = pi->force.pressure;
   const float rhoi = pi->rho;
   const float rhoj = pj->rho;
 
@@ -590,8 +558,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
     Bi[i] = pi->mhd_data.BPred[i];
     Bj[i] = pj->mhd_data.BPred[i];
   }
-  const float B2i = Bi[0] * Bi[0] + Bi[1] * Bi[1] + Bi[2] * Bi[2];
-
+//  const float B2i = Bi[0] * Bi[0] + Bi[1] * Bi[1] + Bi[2] * Bi[2];
+//
 
   ///////////////////////////// FORCE MAXWELL TENSOR
   for (int i = 0; i < 3; i++)
@@ -603,9 +571,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
     mm_i[j][j] -= 0.5 * (Bi[0] * Bi[0] + Bi[1] * Bi[1] + Bi[2] * Bi[2]);
     mm_j[j][j] -= 0.5 * (Bj[0] * Bj[0] + Bj[1] * Bj[1] + Bj[2] * Bj[2]);
   }
-  const float plasma_beta_i = B2i != 0.0f ? 2.0f * mu_0 * Pi / B2i : FLT_MAX;
-  const float scale_i = 0.125f * (10.0f - plasma_beta_i);
-  pi->mhd_data.Q0 = fmaxf(0.0f, fminf(scale_i, 1.0f));
+//  const float plasma_beta_i = B2i != 0.0f ? 2.0f * mu_0 * Pi / B2i : FLT_MAX;
+//  const float scale_i = 0.125f * (10.0f - plasma_beta_i);
+//  pi->mhd_data.Q0 = fmaxf(0.0f, fminf(scale_i, 1.0f));
   //////////////////////////// Apply to the Force and DIVB TERM SUBTRACTION
   for (int i = 0; i < 3; i++)
     for (int j = 0; j < 3; j++) {
@@ -628,7 +596,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
 
   /////////////////////////// VP INDUCTION
   const float mag_VPIndi = wi_dr * r_inv / rhoi;
-  const float mag_VPIndj = wj_dr * r_inv / rhoj;
+//  const float mag_VPIndj = wj_dr * r_inv / rhoj;
   // Normal Gauge
   double dA[3];
   for (int i = 0; i < 3; i++)
@@ -638,15 +606,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
   dv[1] = pi->v[1] - pj->v[1];
   dv[2] = pi->v[2] - pj->v[2];
   
-  const float sourceAi = dv[0] * pi->mhd_data.APred[0] +
-                         dv[1] * pi->mhd_data.APred[1] +
-                         dv[2] * pi->mhd_data.APred[2];
-
- // float SAi = sourceAi + a * a * (pi->mhd_data.Gau - pj->mhd_data.Gau);
-  for (int i = 0; i < 3; i++)
+  const float sAi = dv[0] * pi->mhd_data.APred[0] + dv[1] * pi->mhd_data.APred[1] +
+                    dv[2] * pi->mhd_data.APred[2];
+//  const float sAj = dv[0] * pj->mhd_data.APred[0] + dv[1] * pj->mhd_data.APred[1] + dv[2] * pj->mhd_data.APred[2];
+//  const float sAi = -(dA[0] * pi->v[0] + dA[1] * pi->v[1] + dA[2] * pi->v[2]);
+  const float sourceAi = sAi + a * a * (pi->mhd_data.Gau - pj->mhd_data.Gau);
+//  const float sourceAj =  sAj + a * a * (pi->mhd_data.Gau - pj->mhd_data.Gau);
+  for (int i = 0; i < 3; i++){
     pi->mhd_data.dAdt[i] += mj * mag_VPIndi * sourceAi * dx[i];
-  for (int i = 0; i < 3; i++)
-    pi->mhd_data.dAdt[i] += mj * a * a * (mag_VPIndi * pi->mhd_data.Gau + mag_VPIndj * pj->mhd_data.Gau) * dx[i];
+//    pi->mhd_data.dAdt[i] += mj * dA[i] * mag_VPIndi * (dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2]);
+  }
   /// DISSSIPATION
   const float mag_Disi =
       (wi_dr + wj_dr) / 2.f * r_inv * rhoi / (rho_ij * rho_ij);
@@ -656,7 +625,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
  /* Save induction sources */
   for (int i = 0; i < 3; i++) {
     pi->mhd_data.Adv_A_source[i] += mj * mag_VPIndi * sourceAi * dx[i];
-    pi->mhd_data.Adv_A_source[i] += mj * a * a * (mag_VPIndi * pi->mhd_data.Gau + mag_VPIndj * pj->mhd_data.Gau) * dx[i];
     pi->mhd_data.Diff_A_source[i] += mj * 8.0 * pi->mhd_data.resistive_eta * mag_Disi * dA[i];
     pi->mhd_data.Delta_A[i] += mj * 8.0 * mag_Disi * dA[i];
   }

src/mhd/VPotential/mhd_struct.h

@@ -54,6 +54,7 @@ struct mhd_part_data {
   float Diff_A_source[3];
   /* Laplacian A */
   float Delta_A[3];
+  float MaxVsig;
 };
 
 /**