diff --git a/examples/MHDTests/FastRotor/FastRotor_schemes.yml b/examples/MHDTests/FastRotor/FastRotor_schemes.yml
index f4910bd237aac733d4520e35fb2c24dcb05dedc4..7885ced39a8332e961375ee8056c138d55e78ad1 100644
--- a/examples/MHDTests/FastRotor/FastRotor_schemes.yml
+++ b/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:
diff --git a/examples/MHDTests/FastRotor/makeIC_schemes.py b/examples/MHDTests/FastRotor/makeIC_schemes.py
index 503aa4e1b190f4af5cc08e53212ab85297b5d5b2..17d76aa7a41a529abee38fe535b063a2d9a81fbd 100644
--- a/examples/MHDTests/FastRotor/makeIC_schemes.py
+++ b/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
diff --git a/examples/MHDTests/FastRotor/plot_schemes.py b/examples/MHDTests/FastRotor/plot_schemes.py
index 029b29197d38161f15d1cc560ca33f427d1c94d2..047dd8369bb14fe3e5fe212bdee9da7d1724a248 100644
--- a/examples/MHDTests/FastRotor/plot_schemes.py
+++ b/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
diff --git a/examples/MHDTests/FastRotor/run_schemes.sh b/examples/MHDTests/FastRotor/run_schemes.sh
index 7955066cd13c2e64fb0e9b9b1bf2e22bb3546455..704fe27ddfc07ccceb0b746cda9eb9e5f6c8b120 100755
--- a/examples/MHDTests/FastRotor/run_schemes.sh
+++ b/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
diff --git a/src/mhd/VPotential/mhd.h b/src/mhd/VPotential/mhd.h
index 965e0a529bf2b48879e5beb7d51ba82a9c933ab7..1f79c0de38216c873463373e6d0177d0ad23e0c1 100644
--- a/src/mhd/VPotential/mhd.h
+++ b/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);
diff --git a/src/mhd/VPotential/mhd_iact.h b/src/mhd/VPotential/mhd_iact.h
index 38b0c71d648a632ca6708aef5782106f370ff784..d987db5e6d9bb90968f651ec95836a7ad1150307 100644
--- a/src/mhd/VPotential/mhd_iact.h
+++ b/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];
}
diff --git a/src/mhd/VPotential/mhd_struct.h b/src/mhd/VPotential/mhd_struct.h
index 558b3ef12d9e6d9e61f4e5548d1adb5c60132cf6..9cdbcd0d4406063bcf8032689187061d419d10be 100644
--- a/src/mhd/VPotential/mhd_struct.h
+++ b/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;
};
/**