diff --git a/src/mhd/DirectInduction/mhd.h b/src/mhd/DirectInduction/mhd.h
index df5e9ed867616f778325c9238c7725e83a787d68..6b16f6b4eea79b0da73812c822c8ef5b47c547ce 100644
--- a/src/mhd/DirectInduction/mhd.h
+++ b/src/mhd/DirectInduction/mhd.h
@@ -382,6 +382,8 @@ __attribute__((always_inline)) INLINE static void mhd_prepare_force(
const float h = p->h;
const float rho = p->rho;
+ const float P = p->force.pressure;
+
float B[3];
B[0] = p->mhd_data.B_over_rho[0] * rho;
B[1] = p->mhd_data.B_over_rho[1] * rho;
@@ -401,6 +403,12 @@ __attribute__((always_inline)) INLINE static void mhd_prepare_force(
p->mhd_data.alpha_AR =
normB ? fminf(1.0f, h * sqrtf(grad_B_mean_square) / normB) : 0.0f;
+
+ const float plasma_beta = B2 != 0.0f ? 2.0f * mu_0 * P / B2 : FLT_MAX;
+ const float scale = 0.125f * (10.0f - plasma_beta);
+ const float tensile_correction_scale = fmaxf(0.0f, fminf(scale, 1.0f));
+
+ p->mhd_data.TIC_beta = p->mhd_data.monopole_beta * ensile_correction_scale;
}
/**
diff --git a/src/mhd/DirectInduction/mhd_iact.h b/src/mhd/DirectInduction/mhd_iact.h
index 0af0d5538428fa87337de3248bf179993e2307d9..e325f7f022808489dcf37192f8c3001c466c04c6 100644
--- a/src/mhd/DirectInduction/mhd_iact.h
+++ b/src/mhd/DirectInduction/mhd_iact.h
@@ -298,8 +298,6 @@ __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 rhoi = pi->rho;
const float rhoj = pj->rho;
@@ -429,46 +427,39 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
/* Manifestly *NOT* symmetric in i <-> j */
// const float monopole_beta = hydro_props->mhd.monopole_subtraction;
- const float monopole_beta = pi->mhd_data.monopole_beta;
- 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 TIC_beta_i = pi->mhd_data.TIC_beta;
+ const float TIC_beta_j = pj->mhd_data.TIC_beta;
- const float scale_i = 0.125f * (10.0f - plasma_beta_i);
- const float scale_j = 0.125f * (10.0f - plasma_beta_j);
+ sph_acc_term_i[0] += TIC_beta_i * over_rho2_i * wi_dr * permeability_inv *
+ Bri * r_inv * Bi[0];
+ sph_acc_term_i[0] += TIC_beta_i * over_rho2_j * wj_dr * permeability_inv *
+ Brj * r_inv * Bi[0];
- const float tensile_correction_scale_i = fmaxf(0.0f, fminf(scale_i, 1.0f));
- const float tensile_correction_scale_j = fmaxf(0.0f, fminf(scale_j, 1.0f));
+ sph_acc_term_i[1] += TIC_beta_i * over_rho2_i * wi_dr * permeability_inv *
+ Bri * r_inv * Bi[1];
+ sph_acc_term_i[1] += TIC_beta_i * over_rho2_j * wj_dr * permeability_inv *
+ Brj * r_inv * Bi[1];
- sph_acc_term_i[0] += monopole_beta * over_rho2_i * wi_dr * permeability_inv *
- Bri * r_inv * Bi[0] * tensile_correction_scale_i;
- sph_acc_term_i[0] += monopole_beta * over_rho2_j * wj_dr * permeability_inv *
- Brj * r_inv * Bi[0] * tensile_correction_scale_i;
+ sph_acc_term_i[2] += TIC_beta_i * over_rho2_i * wi_dr * permeability_inv *
+ Bri * r_inv * Bi[2];
+ sph_acc_term_i[2] += TIC_beta_i * over_rho2_j * wj_dr * permeability_inv *
+ Brj * r_inv * Bi[2];
- sph_acc_term_i[1] += monopole_beta * over_rho2_i * wi_dr * permeability_inv *
- Bri * r_inv * Bi[1] * tensile_correction_scale_i;
- sph_acc_term_i[1] += monopole_beta * over_rho2_j * wj_dr * permeability_inv *
- Brj * r_inv * Bi[1] * tensile_correction_scale_i;
+ sph_acc_term_j[0] -= TIC_beta_j * over_rho2_i * wi_dr * permeability_inv *
+ Bri * r_inv * Bj[0];
+ sph_acc_term_j[0] -= TIC_beta_j * over_rho2_j * wj_dr * permeability_inv *
+ Brj * r_inv * Bj[0];
- sph_acc_term_i[2] += monopole_beta * over_rho2_i * wi_dr * permeability_inv *
- Bri * r_inv * Bi[2] * tensile_correction_scale_i;
- sph_acc_term_i[2] += monopole_beta * over_rho2_j * wj_dr * permeability_inv *
- Brj * r_inv * Bi[2] * tensile_correction_scale_i;
+ sph_acc_term_j[1] -= TIC_beta_j * over_rho2_i * wi_dr * permeability_inv *
+ Bri * r_inv * Bj[1];
+ sph_acc_term_j[1] -= TIC_beta_j * over_rho2_j * wj_dr * permeability_inv *
+ Brj * r_inv * Bj[1];
- sph_acc_term_j[0] -= monopole_beta * over_rho2_i * wi_dr * permeability_inv *
- Bri * r_inv * Bj[0] * tensile_correction_scale_j;
- sph_acc_term_j[0] -= monopole_beta * over_rho2_j * wj_dr * permeability_inv *
- Brj * r_inv * Bj[0] * tensile_correction_scale_j;
-
- sph_acc_term_j[1] -= monopole_beta * over_rho2_i * wi_dr * permeability_inv *
- Bri * r_inv * Bj[1] * tensile_correction_scale_j;
- sph_acc_term_j[1] -= monopole_beta * over_rho2_j * wj_dr * permeability_inv *
- Brj * r_inv * Bj[1] * tensile_correction_scale_j;
-
- sph_acc_term_j[2] -= monopole_beta * over_rho2_i * wi_dr * permeability_inv *
- Bri * r_inv * Bj[2] * tensile_correction_scale_j;
- sph_acc_term_j[2] -= monopole_beta * over_rho2_j * wj_dr * permeability_inv *
- Brj * r_inv * Bj[2] * tensile_correction_scale_j;
+ sph_acc_term_j[2] -= TIC_beta_j * over_rho2_i * wi_dr * permeability_inv *
+ Bri * r_inv * Bj[2];
+ sph_acc_term_j[2] -= TIC_beta_j * over_rho2_j * wj_dr * permeability_inv *
+ Brj * r_inv * Bj[2];
/* Use the force Luke ! */
pi->a_hydro[0] -= mj * sph_acc_term_i[0];
@@ -666,7 +657,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
/* Recover some data */
const float mi = pi->mass;
const float mj = pj->mass;
- const float Pi = pi->force.pressure;
const float rhoi = pi->rho;
const float rhoj = pj->rho;
@@ -789,26 +779,23 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
/* Manifestly *NOT* symmetric in i <-> j */
// const float monopole_beta = hydro_props->mhd.monopole_subtraction;
- const float monopole_beta = pi->mhd_data.monopole_beta;
-
- 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);
- const float tensile_correction_scale_i = fmaxf(0.0f, fminf(scale_i, 1.0f));
-
- sph_acc_term_i[0] += monopole_beta * over_rho2_i * wi_dr * permeability_inv *
- Bri * r_inv * Bi[0] * tensile_correction_scale_i;
- sph_acc_term_i[0] += monopole_beta * over_rho2_j * wj_dr * permeability_inv *
- Brj * r_inv * Bi[0] * tensile_correction_scale_i;
-
- sph_acc_term_i[1] += monopole_beta * over_rho2_i * wi_dr * permeability_inv *
- Bri * r_inv * Bi[1] * tensile_correction_scale_i;
- sph_acc_term_i[1] += monopole_beta * over_rho2_j * wj_dr * permeability_inv *
- Brj * r_inv * Bi[1] * tensile_correction_scale_i;
-
- sph_acc_term_i[2] += monopole_beta * over_rho2_i * wi_dr * permeability_inv *
- Bri * r_inv * Bi[2] * tensile_correction_scale_i;
- sph_acc_term_i[2] += monopole_beta * over_rho2_j * wj_dr * permeability_inv *
- Brj * r_inv * Bi[2] * tensile_correction_scale_i;
+
+ const float TIC_beta_i = pi->mhd_data.TIC_beta;
+
+ sph_acc_term_i[0] += TIC_beta_i * over_rho2_i * wi_dr * permeability_inv *
+ Bri * r_inv * Bi[0];
+ sph_acc_term_i[0] += TIC_beta_i * over_rho2_j * wj_dr * permeability_inv *
+ Brj * r_inv * Bi[0];
+
+ sph_acc_term_i[1] += TIC_beta_i * over_rho2_i * wi_dr * permeability_inv *
+ Bri * r_inv * Bi[1];
+ sph_acc_term_i[1] += TIC_beta_i * over_rho2_j * wj_dr * permeability_inv *
+ Brj * r_inv * Bi[1];
+
+ sph_acc_term_i[2] += TIC_beta_i * over_rho2_i * wi_dr * permeability_inv *
+ Bri * r_inv * Bi[2];
+ sph_acc_term_i[2] += TIC_beta_i * over_rho2_j * wj_dr * permeability_inv *
+ Brj * r_inv * Bi[2];
/* Use the force Luke ! */
pi->a_hydro[0] -= mj * sph_acc_term_i[0];
pi->a_hydro[1] -= mj * sph_acc_term_i[1];
diff --git a/src/mhd/DirectInduction/mhd_io.h b/src/mhd/DirectInduction/mhd_io.h
index 2cef812932b403ce886e37ba7c1f6e0415248ac7..29ad782da5b2b505332160df7e7c94d930be7496 100644
--- a/src/mhd/DirectInduction/mhd_io.h
+++ b/src/mhd/DirectInduction/mhd_io.h
@@ -291,36 +291,40 @@ INLINE static int mhd_write_particles(const struct part* parts,
"The curl of Magnetic flux densities of the particles");
list[6] = io_make_output_field(
+ "TIC_beta", FLOAT, 1, UNIT_CONV_NO_UNITS, 1.f, parts, mhd_data.TIC_beta,
+ "Tensile instability correction term coupling strength");
+
+ list[7] = io_make_output_field(
"AlphaAR", FLOAT, 1, UNIT_CONV_NO_UNITS, 1.f, parts, mhd_data.alpha_AR,
"Artificial resistivity switch of the particles");
/* Error metrics */
- list[7] = io_make_output_field_convert_part(
+ list[8] = io_make_output_field_convert_part(
"R0", FLOAT, 1, UNIT_CONV_NO_UNITS, 0, parts, xparts, calculate_R0,
"Classical error metric, indicates places with large divergence. "
"Sensetivity to particle noise depends on signal_to_noise parameter, "
"default is 10 (if 1 - weak noise filtering, if 100 - strong noise "
"filtering)");
- list[8] = io_make_output_field_convert_part(
+ list[9] = io_make_output_field_convert_part(
"R1", FLOAT, 1, UNIT_CONV_NO_UNITS, 0, parts, xparts, calculate_R1,
"Error metric, angle between B field and total Fmag. Indicates unpysical "
"magnetic force. Sensetivity to particle noise depends on "
"signal_to_noise parameter, default is 10 (if 1 - weak noise filtering, "
"if 100 - strong noise filtering)");
- list[9] = io_make_output_field_convert_part(
+ list[10] = io_make_output_field_convert_part(
"R2", FLOAT, 1, UNIT_CONV_NO_UNITS, 0, parts, xparts, calculate_R2,
"Error metric, ratio of divB and |curlB|. Estimates upper limit on "
"B_monopole/B_physical. Sensetivity to particle noise depends on "
"signal_to_noise parameter, default is 10 (if 1 - weak noise filtering, "
"if 100 - strong noise filtering)");
- list[10] = io_make_output_field_convert_part(
+ list[11] = io_make_output_field_convert_part(
"R3", FLOAT, 1, UNIT_CONV_NO_UNITS, 0, parts, xparts, calculate_R3,
"Error metric, shows relation of smoothing length to characteristic B "
"gradient scale. Sensetivity to particle noise depends on "
"signal_to_noise parameter, default is 10 (if 1 - weak noise filtering, "
"if 100 - strong noise filtering)");
- return 11;
+ return 12;
}
/**
diff --git a/src/mhd/DirectInduction/mhd_struct.h b/src/mhd/DirectInduction/mhd_struct.h
index 350b4af294d3d9951cd4d7abd5330b5a7c3d1fa7..6dffc5bde2fa292e0346ab9a85330598a43b9721 100644
--- a/src/mhd/DirectInduction/mhd_struct.h
+++ b/src/mhd/DirectInduction/mhd_struct.h
@@ -46,10 +46,13 @@ struct mhd_part_data {
float psi_over_ch_dt;
- /*! Monopole subtraction in Lorentz Force*/
+ /* Monopole subtraction in Lorentz Force */
float monopole_beta;
- /*! Artifical Diffusion */
+ /* Tensile instability correction coupling strength */
+ float TIC_beta;
+
+ /* Artifical Diffusion */
float art_diff_beta;
/* SPH <1> error */