diff --git a/src/mhd/DirectInduction/mhd.h b/src/mhd/DirectInduction/mhd.h
index 00826cdc2552a4986d637a444347f58e45a3f498..6807133e2456428ab7be89a86884c25d6c76cd9b 100644
--- a/src/mhd/DirectInduction/mhd.h
+++ b/src/mhd/DirectInduction/mhd.h
@@ -133,26 +133,13 @@ __attribute__((always_inline)) INLINE static float mhd_compute_timestep(
const struct hydro_props *hydro_properties, const struct cosmology *cosmo,
const float mu_0) {
- /* Dt from 1/DivOperator(Alfven speed) */
-
- const float divB = p->mhd_data.divB;
-
- const float dt_B_factor = fabsf(divB);
-
- const float dt_B_derivatives =
- dt_B_factor != 0.f
- ? hydro_properties->CFL_condition * cosmo->a /
- cosmo->a_factor_sound_speed *
- sqrtf(p->rho * mu_0 / (dt_B_factor * dt_B_factor))
- : FLT_MAX;
-
const float dt_eta = p->mhd_data.resistive_eta != 0.f
? hydro_properties->CFL_condition * cosmo->a *
cosmo->a * p->h * p->h /
p->mhd_data.resistive_eta
: FLT_MAX;
- return fminf(dt_B_derivatives, dt_eta);
+ return dt_eta;
}
/**
@@ -312,7 +299,6 @@ __attribute__((always_inline)) INLINE static void mhd_reset_gradient(
struct part *p) {
/* Zero the fields updated by the mhd gradient loop */
- p->mhd_data.divB = 0.0f;
p->mhd_data.curl_B[0] = 0.0f;
p->mhd_data.curl_B[1] = 0.0f;
p->mhd_data.curl_B[2] = 0.0f;
@@ -465,6 +451,8 @@ __attribute__((always_inline)) INLINE static void mhd_reset_acceleration(
struct part *restrict p) {
/* Zero the fields updated by the mhd force loop */
+ p->mhd_data.divB = 0.0f;
+
p->mhd_data.B_over_rho_dt[0] = 0.0f;
p->mhd_data.B_over_rho_dt[1] = 0.0f;
p->mhd_data.B_over_rho_dt[2] = 0.0f;
@@ -502,6 +490,8 @@ __attribute__((always_inline)) INLINE static void mhd_reset_predicted_values(
p->mhd_data.B_over_rho[0] = xp->mhd_data.B_over_rho_full[0];
p->mhd_data.B_over_rho[1] = xp->mhd_data.B_over_rho_full[1];
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;
}
/**
@@ -549,8 +539,11 @@ __attribute__((always_inline)) INLINE static void mhd_end_force(
struct part *p, const struct cosmology *cosmo,
const struct hydro_props *hydro_props, const float mu_0) {
- /* Hubble expansion contribution to induction equation */
+ /* Get time derivative of Dedner scalar */
+ p->mhd_data.psi_over_ch_dt = mhd_get_psi_over_ch_dt(
+ p, cosmo->a, cosmo->a_factor_sound_speed, cosmo->H, hydro_props, mu_0);
+ /* Hubble expansion contribution to induction equation */
const float Hubble_induction_pref =
cosmo->a * cosmo->a * cosmo->H * (1.5f * hydro_gamma - 2.f);
p->mhd_data.B_over_rho_dt[0] +=
@@ -592,6 +585,8 @@ __attribute__((always_inline)) INLINE static void mhd_kick_extra(
xp->mhd_data.B_over_rho_full[0] = xp->mhd_data.B_over_rho_full[0] + delta_Bx;
xp->mhd_data.B_over_rho_full[1] = xp->mhd_data.B_over_rho_full[1] + delta_By;
xp->mhd_data.B_over_rho_full[2] = xp->mhd_data.B_over_rho_full[2] + delta_Bz;
+
+ xp->mhd_data.psi_over_ch_full += p->mhd_data.psi_over_ch_dt * dt_therm;
}
/**
diff --git a/src/mhd/DirectInduction/mhd_iact.h b/src/mhd/DirectInduction/mhd_iact.h
index 104a492e662b406d81392c4efb398ab081bceb5e..b89deab5c2fecb7a7b120e697dfe763fe93fc022 100644
--- a/src/mhd/DirectInduction/mhd_iact.h
+++ b/src/mhd/DirectInduction/mhd_iact.h
@@ -119,10 +119,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_gradient(
const float f_ij = 1.f - pi->force.f / mj;
const float f_ji = 1.f - pj->force.f / mi;
- /* 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];
-
/* dB cross r */
float dB_cross_dx[3];
dB_cross_dx[0] = dB[1] * dx[2] - dB[2] * dx[1];
@@ -133,12 +129,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_gradient(
const float over_rho_i = 1.0f / rhoi * f_ij;
const float over_rho_j = 1.0f / rhoj * f_ji;
- /* Calculate monopole term */
- float divB_i = -over_rho_i * (Bri - Brj) * wi_dr * r_inv;
- float divB_j = -over_rho_j * (Bri - Brj) * wj_dr * r_inv;
- pi->mhd_data.divB += mj * divB_i;
- pj->mhd_data.divB += mi * divB_j;
-
/* Calculate curl */
pi->mhd_data.curl_B[0] += mj * over_rho_i * wi_dr * r_inv * dB_cross_dx[0];
pi->mhd_data.curl_B[1] += mj * over_rho_i * wi_dr * r_inv * dB_cross_dx[1];
@@ -234,10 +224,6 @@ runner_iact_nonsym_mhd_gradient(const float r2, const float dx[3],
const float f_ij = 1.f - pi->force.f / mj;
// const float f_ji = 1.f - pj->force.f / mi;
- /* 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]);
-
/* dB cross r */
float dB_cross_dx[3];
dB_cross_dx[0] = dB[1] * dx[2] - dB[2] * dx[1];
@@ -247,10 +233,6 @@ runner_iact_nonsym_mhd_gradient(const float r2, const float dx[3],
/* Compute gradient terms */
const float over_rho_i = 1.0f / rhoi * f_ij;
- /* Calculate monopole term */
- float divB_i = -over_rho_i * (Bri - Brj) * wi_dr * r_inv;
- pi->mhd_data.divB += mj * divB_i;
-
/* Calculate curl */
pi->mhd_data.curl_B[0] += mj * over_rho_i * wi_dr * r_inv * dB_cross_dx[0];
pi->mhd_data.curl_B[1] += mj * over_rho_i * wi_dr * r_inv * dB_cross_dx[1];
@@ -319,9 +301,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
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 normBi = sqrtf(B2i);
- const float normBj = sqrtf(B2j);
-
/*
float curlBi[3];
float curlBj[3];
@@ -373,6 +352,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
const float over_rho2_i = 1.0f / (rhoi * rhoi) * f_ij;
const float over_rho2_j = 1.0f / (rhoj * rhoj) * f_ji;
+ /* Compute symmetric div(B) */
+ 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;
+
+ pi->mhd_data.divB += rhoi * mj * divB_i;
+ pj->mhd_data.divB -= rhoj * mi * divB_j;
+
/* SPH acceleration term in x direction, i_th particle */
float sph_acc_term_i[3] = {0.f, 0.f, 0.f};
@@ -621,35 +611,21 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
/*Divergence diffusion */
- // const float vsig_Dedner_i = pi->viscosity.v_sig;
- // const float vsig_Dedner_j = pj->viscosity.v_sig;
-
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);
- float grad_psi_i =
- over_rho2_i * psi_over_ch_i * vsig_Dedner_i * wi_dr * r_inv;
- grad_psi_i += over_rho2_j * psi_over_ch_j * vsig_Dedner_j * wj_dr * r_inv;
- float grad_psi_j = grad_psi_i;
-
- const float psi_over_ch_i_inv =
- psi_over_ch_i != 0.f ? 1.f / psi_over_ch_i : 0.;
- const float psi_over_ch_j_inv =
- psi_over_ch_j != 0.f ? 1.f / psi_over_ch_j : 0.;
-
- const float corr_ratio_i = fabsf(normBi * psi_over_ch_i_inv);
- const float corr_ratio_j = fabsf(normBj * psi_over_ch_j_inv);
+ 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;
- const float Qi = corr_ratio_i < 2 ? 0.5 * corr_ratio_i : 1.0f;
- const float Qj = corr_ratio_j < 2 ? 0.5 * corr_ratio_j : 1.0f;
+ 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 * Qi * grad_psi_i * dx[0];
- pi->mhd_data.B_over_rho_dt[1] -= mj * Qi * grad_psi_i * dx[1];
- pi->mhd_data.B_over_rho_dt[2] -= mj * Qi * grad_psi_i * dx[2];
-
- pj->mhd_data.B_over_rho_dt[0] += mi * Qj * grad_psi_j * dx[0];
- pj->mhd_data.B_over_rho_dt[1] += mi * Qj * grad_psi_j * dx[1];
- pj->mhd_data.B_over_rho_dt[2] += mi * Qj * grad_psi_j * 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];
/* Save induction sources */
@@ -659,8 +635,8 @@ __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 * Qi * grad_psi_i * dx[i];
- pj->mhd_data.Adv_B_source[i] += mi * Qj * grad_psi_j * dx[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] += 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];
@@ -716,8 +692,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
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 normBi = sqrtf(B2i);
-
/*
float curlBi[3];
float curlBj[3];
@@ -769,6 +743,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
const float over_rho2_i = 1.0f / (rhoi * rhoi) * f_ij;
const float over_rho2_j = 1.0f / (rhoj * rhoj) * f_ji;
+ /* Compute symmetric div(B) */
+ 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;
+
+ pi->mhd_data.divB += rhoi * mj * divB_i;
+
/* SPH acceleration term in x direction, i_th particle */
float sph_acc_term_i[3] = {0.f, 0.f, 0.f};
@@ -945,33 +928,23 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
/*Divergence diffusion */
- // const float vsig_Dedner_i = pi->viscosity.v_sig;
- // const float vsig_Dedner_j = pj->viscosity.v_sig;
-
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);
- float grad_psi_i =
- over_rho2_i * psi_over_ch_i * vsig_Dedner_i * wi_dr * r_inv;
- grad_psi_i += over_rho2_j * psi_over_ch_j * vsig_Dedner_j * wj_dr * r_inv;
-
- const float psi_over_ch_i_inv =
- psi_over_ch_i != 0.f ? 1.f / psi_over_ch_i : 0.;
-
- const float corr_ratio_i = fabsf(normBi * psi_over_ch_i_inv);
-
- const float Qi = corr_ratio_i < 2 ? 0.5 * corr_ratio_i : 1.0f;
+ 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;
- pi->mhd_data.B_over_rho_dt[0] -= mj * Qi * grad_psi_i * dx[0];
- pi->mhd_data.B_over_rho_dt[1] -= mj * Qi * grad_psi_i * dx[1];
- pi->mhd_data.B_over_rho_dt[2] -= mj * Qi * grad_psi_i * dx[2];
+ 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];
/* 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 * Qi * grad_psi_i * dx[i];
+ pi->mhd_data.Adv_B_source[i] -= mj * grad_psi_i * 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];