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

Fast magnetosonic wave speed carried by particles now

parent dbf1c093
No related branches found
No related tags found
No related merge requests found
...@@ -191,12 +191,12 @@ __attribute__((always_inline)) INLINE static float mhd_get_comoving_magnetosonic ...@@ -191,12 +191,12 @@ __attribute__((always_inline)) INLINE static float mhd_get_comoving_magnetosonic
/* Get Alfven speed squared */ /* Get Alfven speed squared */
const float vA = p->mhd_data.vA; const float vA = p->mhd_data.vA;
const float v_A2 = vA * vA; const float vA2 = vA * vA;
/* Compute effective sound speeds */ /* Compute effective sound speeds */
const float c_ms2 = cs2 + v_A2; const float cfms2 = cs2 + vA2;
return sqrtf(c_ms2); return sqrtf(cfms2);
} }
/** /**
...@@ -226,7 +226,7 @@ mhd_get_fast_magnetosonic_wave_phase_velocity(const float dx[3], ...@@ -226,7 +226,7 @@ mhd_get_fast_magnetosonic_wave_phase_velocity(const float dx[3],
/* Compute effective sound speeds */ /* Compute effective sound speeds */
const float cs = p->force.soundspeed; const float cs = p->force.soundspeed;
const float cs2 = cs * cs; const float cs2 = cs * cs;
const float c_ms = mhd_get_comoving_magnetosonic_speed(p); const float c_ms = p->mhd_data.c_fms;
const float c_ms2 = c_ms * c_ms; const float c_ms2 = c_ms * c_ms;
const float projection_correction = c_ms2 * c_ms2 - 4.0f * permeability_inv * const float projection_correction = c_ms2 * c_ms2 - 4.0f * permeability_inv *
cs2 * Br * r_inv * cs2 * Br * r_inv *
...@@ -255,8 +255,8 @@ __attribute__((always_inline)) INLINE static float mhd_signal_velocity( ...@@ -255,8 +255,8 @@ __attribute__((always_inline)) INLINE static float mhd_signal_velocity(
const float dx[3], const struct part *restrict pi, const float dx[3], const struct part *restrict pi,
const struct part *restrict pj, const float mu_ij, const float beta) { const struct part *restrict pj, const float mu_ij, const float beta) {
const float v_sigi = mhd_get_comoving_magnetosonic_speed(pi); const float v_sigi = pi->mhd_data.c_fms;
const float v_sigj = mhd_get_comoving_magnetosonic_speed(pj); const float v_sigj = pj->mhd_data.c_fms;
const float v_sig = v_sigi + v_sigj - beta * mu_ij; const float v_sig = v_sigi + v_sigj - beta * mu_ij;
...@@ -306,7 +306,8 @@ __attribute__((always_inline)) INLINE static void mhd_prepare_gradient( ...@@ -306,7 +306,8 @@ __attribute__((always_inline)) INLINE static void mhd_prepare_gradient(
const struct cosmology *cosmo, const struct hydro_props *hydro_props, const float mu_0) { const struct cosmology *cosmo, const struct hydro_props *hydro_props, const float mu_0) {
p->mhd_data.vA = mhd_get_comoving_alfven_speed(p, mu_0); p->mhd_data.vA = mhd_get_comoving_alfven_speed(p, mu_0);
p->mhd_data.c_fms = mhd_get_comoving_magnetosonic_speed(p);
} }
/** /**
...@@ -390,7 +391,7 @@ __attribute__((always_inline)) INLINE static float mhd_get_psi_over_ch_dt( ...@@ -390,7 +391,7 @@ __attribute__((always_inline)) INLINE static float mhd_get_psi_over_ch_dt(
const float h_inv = 1.0f / h; const float h_inv = 1.0f / h;
/* Compute Dedner cleaning speed. */ /* Compute Dedner cleaning speed. */
const float ch = mhd_get_comoving_magnetosonic_speed(p); const float ch = p->mhd_data.c_fms;
/* Compute Dedner cleaning scalar time derivative. */ /* Compute Dedner cleaning scalar time derivative. */
const float hyp = hydro_props->mhd.hyp_dedner; const float hyp = hydro_props->mhd.hyp_dedner;
...@@ -516,7 +517,7 @@ __attribute__((always_inline)) INLINE static void mhd_reset_predicted_values( ...@@ -516,7 +517,7 @@ __attribute__((always_inline)) INLINE static void mhd_reset_predicted_values(
/* Set Alfven speed */ /* Set Alfven speed */
p->mhd_data.vA = mhd_get_comoving_alfven_speed(p, mu_0); p->mhd_data.vA = mhd_get_comoving_alfven_speed(p, mu_0);
p->mhd_data.c_fms = mhd_get_comoving_magnetosonic_speed(p);
} }
/** /**
...@@ -548,6 +549,7 @@ __attribute__((always_inline)) INLINE static void mhd_predict_extra( ...@@ -548,6 +549,7 @@ __attribute__((always_inline)) INLINE static void mhd_predict_extra(
p->mhd_data.psi_over_ch += p->mhd_data.psi_over_ch_dt * dt_therm; p->mhd_data.psi_over_ch += p->mhd_data.psi_over_ch_dt * dt_therm;
p->mhd_data.vA = mhd_get_comoving_alfven_speed(p, mu_0); p->mhd_data.vA = mhd_get_comoving_alfven_speed(p, mu_0);
p->mhd_data.c_fms = mhd_get_comoving_magnetosonic_speed(p);
} }
...@@ -647,6 +649,7 @@ __attribute__((always_inline)) INLINE static void mhd_convert_quantities( ...@@ -647,6 +649,7 @@ __attribute__((always_inline)) INLINE static void mhd_convert_quantities(
p->mhd_data.B_over_rho[2] *= pow(cosmo->a, 1.5f * hydro_gamma); p->mhd_data.B_over_rho[2] *= pow(cosmo->a, 1.5f * hydro_gamma);
p->mhd_data.vA = mhd_get_comoving_alfven_speed(p, mu_0); p->mhd_data.vA = mhd_get_comoving_alfven_speed(p, mu_0);
p->mhd_data.c_fms = mhd_get_comoving_magnetosonic_speed(p);
xp->mhd_data.B_over_rho_full[0] = p->mhd_data.B_over_rho[0]; xp->mhd_data.B_over_rho_full[0] = p->mhd_data.B_over_rho[0];
xp->mhd_data.B_over_rho_full[1] = p->mhd_data.B_over_rho[1]; xp->mhd_data.B_over_rho_full[1] = p->mhd_data.B_over_rho[1];
......
...@@ -33,12 +33,12 @@ __attribute__((always_inline)) INLINE static void mhd_debug_particle( ...@@ -33,12 +33,12 @@ __attribute__((always_inline)) INLINE static void mhd_debug_particle(
"B/rho=[%.3e,%.3e,%.3e], " "B/rho=[%.3e,%.3e,%.3e], "
"B/rho (pred)=[%.3e,%.3e,%.3e], " "B/rho (pred)=[%.3e,%.3e,%.3e], "
"d(B/rho)/dt=[%.3e,%.3e,%.3e], \n" "d(B/rho)/dt=[%.3e,%.3e,%.3e], \n"
"divB=%.3e, v_fm=%.3e, psi/cs=%.3e, d(psi/c)/dt=%.3e\n", "divB=%.3e, c_fms=%.3e, psi/ch=%.3e, d(psi/ch)/dt=%.3e\n",
xp->mhd_data.B_over_rho_full[0], xp->mhd_data.B_over_rho_full[1], xp->mhd_data.B_over_rho_full[0], xp->mhd_data.B_over_rho_full[1],
xp->mhd_data.B_over_rho_full[2], p->mhd_data.B_over_rho[0], xp->mhd_data.B_over_rho_full[2], p->mhd_data.B_over_rho[0],
p->mhd_data.B_over_rho[1], p->mhd_data.B_over_rho[2], p->mhd_data.B_over_rho[1], p->mhd_data.B_over_rho[2],
p->mhd_data.B_over_rho_dt[0], p->mhd_data.B_over_rho_dt[1], p->mhd_data.B_over_rho_dt[0], p->mhd_data.B_over_rho_dt[1],
p->mhd_data.B_over_rho_dt[2], p->mhd_data.divB, p->mhd_data.v_fm, p->mhd_data.B_over_rho_dt[2], p->mhd_data.divB, p->mhd_data.c_fms,
p->mhd_data.psi_over_ch, p->mhd_data.psi_over_ch_dt); p->mhd_data.psi_over_ch, p->mhd_data.psi_over_ch_dt);
} }
......
...@@ -625,8 +625,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force( ...@@ -625,8 +625,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
// const float vsig_Dedner_i = pi->viscosity.v_sig; // const float vsig_Dedner_i = pi->viscosity.v_sig;
// const float vsig_Dedner_j = pj->viscosity.v_sig; // const float vsig_Dedner_j = pj->viscosity.v_sig;
const float vsig_Dedner_i = mhd_get_comoving_magnetosonic_speed(pi); const float vsig_Dedner_i = pi->mhd_data.c_fms;
const float vsig_Dedner_j = mhd_get_comoving_magnetosonic_speed(pj); const float vsig_Dedner_j = pj->mhd_data.c_fms;
float grad_psi_i = float grad_psi_i =
over_rho2_i * psi_over_ch_i * vsig_Dedner_i * wi_dr * r_inv; over_rho2_i * psi_over_ch_i * vsig_Dedner_i * wi_dr * r_inv;
...@@ -949,8 +949,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force( ...@@ -949,8 +949,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
// const float vsig_Dedner_i = pi->viscosity.v_sig; // const float vsig_Dedner_i = pi->viscosity.v_sig;
// const float vsig_Dedner_j = pj->viscosity.v_sig; // const float vsig_Dedner_j = pj->viscosity.v_sig;
const float vsig_Dedner_i = mhd_get_comoving_magnetosonic_speed(pi); const float vsig_Dedner_i = pi->mhd_data.c_fms;
const float vsig_Dedner_j = mhd_get_comoving_magnetosonic_speed(pj); const float vsig_Dedner_j = pj->mhd_data.c_fms;
float grad_psi_i = float grad_psi_i =
over_rho2_i * psi_over_ch_i * vsig_Dedner_i * wi_dr * r_inv; over_rho2_i * psi_over_ch_i * vsig_Dedner_i * wi_dr * r_inv;
......
...@@ -31,7 +31,7 @@ struct mhd_part_data { ...@@ -31,7 +31,7 @@ struct mhd_part_data {
/*! dB Direct Induction */ /*! dB Direct Induction */
float B_over_rho_dt[3]; float B_over_rho_dt[3];
float v_fm; float c_fms;
float vA; float vA;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment