Skip to content
Snippets Groups Projects
Commit 11a3ab93 authored by Orestis Karapiperis's avatar Orestis Karapiperis
Browse files

Added tensile instability correction term coupling strength to snapshots

parent 6cf519a2
No related branches found
No related tags found
1 merge request!1896Karapiperis/tic in snapshots
Hide whitespace changes
Inline Side-by-side
src/mhd/DirectInduction/mhd.h
+ 8
− 0
View file @ 11a3ab93
...@@ -382,6 +382,8 @@ __attribute__((always_inline)) INLINE static void mhd_prepare_force( ...@@ -382,6 +382,8 @@ __attribute__((always_inline)) INLINE static void mhd_prepare_force(
const float h = p->h; const float h = p->h;
const float rho = p->rho; const float rho = p->rho;
const float P = p->force.pressure;
float B[3]; float B[3];
B[0] = p->mhd_data.B_over_rho[0] * rho; B[0] = p->mhd_data.B_over_rho[0] * rho;
B[1] = p->mhd_data.B_over_rho[1] * rho; B[1] = p->mhd_data.B_over_rho[1] * rho;
...@@ -401,6 +403,12 @@ __attribute__((always_inline)) INLINE static void mhd_prepare_force( ...@@ -401,6 +403,12 @@ __attribute__((always_inline)) INLINE static void mhd_prepare_force(
p->mhd_data.alpha_AR = p->mhd_data.alpha_AR =
normB ? fminf(1.0f, h * sqrtf(grad_B_mean_square) / normB) : 0.0f; 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;
} }
/** /**
......
src/mhd/DirectInduction/mhd_iact.h
+ 43
− 56
View file @ 11a3ab93
...@@ -298,8 +298,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force( ...@@ -298,8 +298,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
/* Recover some data */ /* Recover some data */
const float mi = pi->mass; const float mi = pi->mass;
const float mj = pj->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 rhoi = pi->rho;
const float rhoj = pj->rho; const float rhoj = pj->rho;
...@@ -429,46 +427,39 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force( ...@@ -429,46 +427,39 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
/* Manifestly *NOT* symmetric in i <-> j */ /* Manifestly *NOT* symmetric in i <-> j */
// const float monopole_beta = hydro_props->mhd.monopole_subtraction; // 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 TIC_beta_i = pi->mhd_data.TIC_beta;
const float plasma_beta_j = B2j != 0.0f ? 2.0f * mu_0 * Pj / B2j : FLT_MAX; const float TIC_beta_j = pj->mhd_data.TIC_beta;
const float scale_i = 0.125f * (10.0f - plasma_beta_i); sph_acc_term_i[0] += TIC_beta_i * over_rho2_i * wi_dr * permeability_inv *
const float scale_j = 0.125f * (10.0f - plasma_beta_j); 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)); sph_acc_term_i[1] += TIC_beta_i * over_rho2_i * wi_dr * permeability_inv *
const float tensile_correction_scale_j = fmaxf(0.0f, fminf(scale_j, 1.0f)); 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 * sph_acc_term_i[2] += TIC_beta_i * over_rho2_i * wi_dr * permeability_inv *
Bri * r_inv * Bi[0] * tensile_correction_scale_i; Bri * r_inv * Bi[2];
sph_acc_term_i[0] += monopole_beta * over_rho2_j * wj_dr * permeability_inv * sph_acc_term_i[2] += TIC_beta_i * over_rho2_j * wj_dr * permeability_inv *
Brj * r_inv * Bi[0] * tensile_correction_scale_i; Brj * r_inv * Bi[2];
sph_acc_term_i[1] += monopole_beta * over_rho2_i * wi_dr * permeability_inv * sph_acc_term_j[0] -= TIC_beta_j * over_rho2_i * wi_dr * permeability_inv *
Bri * r_inv * Bi[1] * tensile_correction_scale_i; Bri * r_inv * Bj[0];
sph_acc_term_i[1] += monopole_beta * over_rho2_j * wj_dr * permeability_inv * sph_acc_term_j[0] -= TIC_beta_j * over_rho2_j * wj_dr * permeability_inv *
Brj * r_inv * Bi[1] * tensile_correction_scale_i; Brj * r_inv * Bj[0];
sph_acc_term_i[2] += monopole_beta * over_rho2_i * wi_dr * permeability_inv * sph_acc_term_j[1] -= TIC_beta_j * over_rho2_i * wi_dr * permeability_inv *
Bri * r_inv * Bi[2] * tensile_correction_scale_i; Bri * r_inv * Bj[1];
sph_acc_term_i[2] += monopole_beta * over_rho2_j * wj_dr * permeability_inv * sph_acc_term_j[1] -= TIC_beta_j * over_rho2_j * wj_dr * permeability_inv *
Brj * r_inv * Bi[2] * tensile_correction_scale_i; Brj * r_inv * Bj[1];
sph_acc_term_j[0] -= monopole_beta * over_rho2_i * wi_dr * permeability_inv * sph_acc_term_j[2] -= TIC_beta_j * over_rho2_i * wi_dr * permeability_inv *
Bri * r_inv * Bj[0] * tensile_correction_scale_j; Bri * r_inv * Bj[2];
sph_acc_term_j[0] -= monopole_beta * over_rho2_j * wj_dr * permeability_inv * sph_acc_term_j[2] -= TIC_beta_j * over_rho2_j * wj_dr * permeability_inv *
Brj * r_inv * Bj[0] * tensile_correction_scale_j; Brj * r_inv * Bj[2];
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;
/* Use the force Luke ! */ /* Use the force Luke ! */
pi->a_hydro[0] -= mj * sph_acc_term_i[0]; 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( ...@@ -666,7 +657,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
/* Recover some data */ /* Recover some data */
const float mi = pi->mass; const float mi = pi->mass;
const float mj = pj->mass; const float mj = pj->mass;
const float Pi = pi->force.pressure;
const float rhoi = pi->rho; const float rhoi = pi->rho;
const float rhoj = pj->rho; const float rhoj = pj->rho;
...@@ -789,26 +779,23 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force( ...@@ -789,26 +779,23 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
/* Manifestly *NOT* symmetric in i <-> j */ /* Manifestly *NOT* symmetric in i <-> j */
// const float monopole_beta = hydro_props->mhd.monopole_subtraction; // const float monopole_beta = hydro_props->mhd.monopole_subtraction;
const float monopole_beta = pi->mhd_data.monopole_beta;
const float TIC_beta_i = pi->mhd_data.TIC_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); sph_acc_term_i[0] += TIC_beta_i * over_rho2_i * wi_dr * permeability_inv *
const float tensile_correction_scale_i = fmaxf(0.0f, fminf(scale_i, 1.0f)); Bri * r_inv * Bi[0];
sph_acc_term_i[0] += TIC_beta_i * over_rho2_j * wj_dr * permeability_inv *
sph_acc_term_i[0] += monopole_beta * over_rho2_i * wi_dr * permeability_inv * Brj * r_inv * Bi[0];
Bri * r_inv * Bi[0] * tensile_correction_scale_i;
sph_acc_term_i[0] += monopole_beta * over_rho2_j * wj_dr * permeability_inv * sph_acc_term_i[1] += TIC_beta_i * over_rho2_i * wi_dr * permeability_inv *
Brj * r_inv * Bi[0] * tensile_correction_scale_i; Bri * r_inv * Bi[1];
sph_acc_term_i[1] += TIC_beta_i * over_rho2_j * wj_dr * permeability_inv *
sph_acc_term_i[1] += monopole_beta * over_rho2_i * wi_dr * permeability_inv * Brj * r_inv * Bi[1];
Bri * r_inv * Bi[1] * tensile_correction_scale_i;
sph_acc_term_i[1] += monopole_beta * over_rho2_j * wj_dr * permeability_inv * sph_acc_term_i[2] += TIC_beta_i * over_rho2_i * wi_dr * permeability_inv *
Brj * r_inv * Bi[1] * tensile_correction_scale_i; Bri * r_inv * Bi[2];
sph_acc_term_i[2] += TIC_beta_i * over_rho2_j * wj_dr * permeability_inv *
sph_acc_term_i[2] += monopole_beta * over_rho2_i * wi_dr * permeability_inv * Brj * r_inv * Bi[2];
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;
/* Use the force Luke ! */ /* Use the force Luke ! */
pi->a_hydro[0] -= mj * sph_acc_term_i[0]; pi->a_hydro[0] -= mj * sph_acc_term_i[0];
pi->a_hydro[1] -= mj * sph_acc_term_i[1]; pi->a_hydro[1] -= mj * sph_acc_term_i[1];
......
src/mhd/DirectInduction/mhd_io.h
+ 9
− 5
View file @ 11a3ab93
...@@ -291,36 +291,40 @@ INLINE static int mhd_write_particles(const struct part* parts, ...@@ -291,36 +291,40 @@ INLINE static int mhd_write_particles(const struct part* parts,
"The curl of Magnetic flux densities of the particles"); "The curl of Magnetic flux densities of the particles");
list[6] = io_make_output_field( 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, "AlphaAR", FLOAT, 1, UNIT_CONV_NO_UNITS, 1.f, parts, mhd_data.alpha_AR,
"Artificial resistivity switch of the particles"); "Artificial resistivity switch of the particles");
/* Error metrics */ /* 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, "R0", FLOAT, 1, UNIT_CONV_NO_UNITS, 0, parts, xparts, calculate_R0,
"Classical error metric, indicates places with large divergence. " "Classical error metric, indicates places with large divergence. "
"Sensetivity to particle noise depends on signal_to_noise parameter, " "Sensetivity to particle noise depends on signal_to_noise parameter, "
"default is 10 (if 1 - weak noise filtering, if 100 - strong noise " "default is 10 (if 1 - weak noise filtering, if 100 - strong noise "
"filtering)"); "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, "R1", FLOAT, 1, UNIT_CONV_NO_UNITS, 0, parts, xparts, calculate_R1,
"Error metric, angle between B field and total Fmag. Indicates unpysical " "Error metric, angle between B field and total Fmag. Indicates unpysical "
"magnetic force. Sensetivity to particle noise depends on " "magnetic force. Sensetivity to particle noise depends on "
"signal_to_noise parameter, default is 10 (if 1 - weak noise filtering, " "signal_to_noise parameter, default is 10 (if 1 - weak noise filtering, "
"if 100 - strong 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, "R2", FLOAT, 1, UNIT_CONV_NO_UNITS, 0, parts, xparts, calculate_R2,
"Error metric, ratio of divB and |curlB|. Estimates upper limit on " "Error metric, ratio of divB and |curlB|. Estimates upper limit on "
"B_monopole/B_physical. Sensetivity to particle noise depends on " "B_monopole/B_physical. Sensetivity to particle noise depends on "
"signal_to_noise parameter, default is 10 (if 1 - weak noise filtering, " "signal_to_noise parameter, default is 10 (if 1 - weak noise filtering, "
"if 100 - strong 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, "R3", FLOAT, 1, UNIT_CONV_NO_UNITS, 0, parts, xparts, calculate_R3,
"Error metric, shows relation of smoothing length to characteristic B " "Error metric, shows relation of smoothing length to characteristic B "
"gradient scale. Sensetivity to particle noise depends on " "gradient scale. Sensetivity to particle noise depends on "
"signal_to_noise parameter, default is 10 (if 1 - weak noise filtering, " "signal_to_noise parameter, default is 10 (if 1 - weak noise filtering, "
"if 100 - strong noise filtering)"); "if 100 - strong noise filtering)");
return 11; return 12;
} }
/** /**
......
src/mhd/DirectInduction/mhd_struct.h
+ 5
− 2
View file @ 11a3ab93
...@@ -46,10 +46,13 @@ struct mhd_part_data { ...@@ -46,10 +46,13 @@ struct mhd_part_data {
float psi_over_ch_dt; float psi_over_ch_dt;
/*! Monopole subtraction in Lorentz Force*/ /* Monopole subtraction in Lorentz Force */
float monopole_beta; float monopole_beta;
/*! Artifical Diffusion */ /* Tensile instability correction coupling strength */
float TIC_beta;
/* Artifical Diffusion */
float art_diff_beta; float art_diff_beta;
/* SPH <1> error */ /* SPH <1> error */
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment