diff --git a/examples/MHDTests/RobertsFlow/plot_snapshot_slice.py b/examples/MHDTests/RobertsFlow/plot_snapshot_slice.py
new file mode 100755
index 0000000000000000000000000000000000000000..f59ebed1399d925dfd2a64c866425a6afa2e6145
--- /dev/null
+++ b/examples/MHDTests/RobertsFlow/plot_snapshot_slice.py
@@ -0,0 +1,355 @@
+from swiftsimio import load
+from swiftsimio.visualisation.slice import slice_gas
+from swiftsimio.visualisation.rotation import rotation_matrix_from_vector
+import numpy as np
+import sys
+import h5py
+
+filename = sys.argv[1]
+
+slice_height = sys.argv[3]
+
+projection = sys.argv[4]
+
+prefered_color = 'magma'
+
+cpts = 100
+
+
+to_plot = 'errors' # 'B' or 'A' or 'errors'
+
+with h5py.File(filename, "r") as handle:
+ gamma = handle["HydroScheme"].attrs["Adiabatic index"][0]
+ boxsize = handle["Header"].attrs["BoxSize"][0]
+ t = handle["Header"].attrs["Time"][0]
+ git = handle["Code"].attrs["Git Revision"]
+ gitBranch = handle["Code"].attrs["Git Branch"]
+ scheme = handle["/HydroScheme"].attrs["Scheme"]
+ kernel = handle["/HydroScheme"].attrs["Kernel function"]
+ neighbours = handle["/HydroScheme"].attrs["Kernel target N_ngb"]
+ mu0 = handle["/PhysicalConstants/InternalUnits"].attrs["vacuum_permeability"]
+ #dedhyp = handle["/HydroScheme"].attrs["Dedner Hyperbolic Constant"]
+ #dedpar = handle["/HydroScheme"].attrs["Dedner Parabolic Constant"]
+ mhdflavour = handle["/HydroScheme"].attrs["MHD Flavour"]
+
+
+data = load(filename)
+
+# getting values from snapshots
+
+img_res = 512
+divreg = 1e-30
+above_noise=10
+reg_err=0.01
+err_upper_bnd=0.999
+
+if projection=='yz':
+ rotation_matrix = np.array([[0,1,0],[0,0,1],[1,0,0]])
+ for i in range(len(data.gas.coordinates)):
+ data.gas.coordinates[i] = np.matmul(rotation_matrix, data.gas.coordinates[i])
+ data.gas.velocities[i] = np.matmul(rotation_matrix, data.gas.velocities[i])
+ data.gas.magnetic_flux_densities[i] = np.matmul(rotation_matrix, data.gas.magnetic_flux_densities[i])
+ data.gas.magnetic_vector_potentials[i] =np.matmul(rotation_matrix, data.gas.magnetic_vector_potentials[i])
+if projection=='xz':
+ rotation_matrix = np.array([[1,0,0],[0,0,1],[0,1,0]])
+ for i in range(len(data.gas.coordinates)):
+ data.gas.coordinates[i] = np.matmul(rotation_matrix, data.gas.coordinates[i])
+ data.gas.velocities[i] = np.matmul(rotation_matrix, data.gas.velocities[i])
+ data.gas.magnetic_flux_densities[i] = np.matmul(rotation_matrix, data.gas.magnetic_flux_densities[i])
+ data.gas.magnetic_vector_potentials[i] =np.matmul(rotation_matrix, data.gas.magnetic_vector_potentials[i])
+
+def abs_vec(vec):
+ res = np.sqrt(vec[:, 0] ** 2 + vec[:, 1] ** 2 + vec[:, 2] ** 2)
+ return res
+def dot_vec(vec1,vec2):
+ res = vec1[:,0] * vec2[:,0] + vec1[:,1] * vec2[:,1] + vec1[:,2] * vec2[:,2]
+ return res
+def cross_vec(vec1,vec2):
+ res_vec = np.zeros((len(vec1),3))
+ res_vec[:,0] = vec1[:,1]*vec2[:,2]-vec1[:,2]*vec2[:,1]
+ res_vec[:,1] = vec1[:,2]*vec2[:,0]-vec1[:,0]*vec2[:,2]
+ res_vec[:,2] = vec1[:,0]*vec2[:,1]-vec1[:,1]*vec2[:,0]
+ return res_vec
+def rms_vec(vec):
+ res = np.sqrt(np.mean(vec[:, 0] ** 2 + vec[:, 1] ** 2 + vec[:, 2] ** 2))
+ return res
+# see available fields in snapshot
+#print(data.metadata.gas_properties.field_names)
+
+# Get physical quantities
+x = data.gas.coordinates[:, 0].value
+y = data.gas.coordinates[:, 1].value
+z = data.gas.coordinates[:, 2].value
+rho = data.gas.densities.value
+h = data.gas.smoothing_lengths.value
+v = data.gas.velocities.value
+P = data.gas.pressures.value
+B = data.gas.magnetic_flux_densities.value
+
+R0 = data.gas.r0.value
+R1 = data.gas.r1.value
+R2 = data.gas.r2.value
+R3 = data.gas.r3.value
+
+#Get RMS values
+Brms = rms_vec(B)
+if to_plot=='A':
+ A = data.gas.magnetic_vector_potentials.value
+ Arms = rms_vec(A)
+
+#Print Brms
+#print(f'Brms = {Brms}')
+
+def make_slice(key, slice_frac_z=float(slice_height)):
+ res = slice_gas(
+ data,
+ z_slice = slice_frac_z *data.metadata.boxsize[2],
+ resolution=img_res,
+ project=key,
+ parallel=True,
+ periodic=True,
+ )
+ return res
+
+masses = make_slice("masses").value
+dimy = len(masses)
+dimx = len(masses[0])
+mass_map = masses.flatten()
+mass_map = mass_map + divreg
+
+l = len(mass_map)
+
+def prepare_sliced_quantity(quantity, isvec=False):
+ if isvec:
+ data.gas.mass_weighted_temp_qx = data.gas.masses * quantity[:,0]
+ data.gas.mass_weighted_temp_qy = data.gas.masses * quantity[:,1]
+ data.gas.mass_weighted_temp_qz = data.gas.masses * quantity[:,2]
+ sliced_quantity = np.zeros((l, 3))
+ sliced_quantity[:,0] = make_slice("mass_weighted_temp_qx").value.flatten() / mass_map
+ sliced_quantity[:,1] = make_slice("mass_weighted_temp_qy").value.flatten() / mass_map
+ sliced_quantity[:,2] = make_slice("mass_weighted_temp_qz").value.flatten() / mass_map
+ else:
+ data.gas.mass_weighted_temp_q = data.gas.masses * quantity[:]
+ sliced_quantity = make_slice("mass_weighted_temp_q").value.flatten() / mass_map
+
+ return sliced_quantity
+
+
+#slicing scalars
+rho = prepare_sliced_quantity(rho)
+h = prepare_sliced_quantity(h)
+P = prepare_sliced_quantity(P)
+
+#slicing vectors
+v = prepare_sliced_quantity(v, isvec=True)
+B = prepare_sliced_quantity(B, isvec=True)
+if to_plot=='A':
+ A = prepare_sliced_quantity(A, isvec=True)
+
+#Get sliced error metrics
+R0 = prepare_sliced_quantity(R0)
+R1 = prepare_sliced_quantity(R1)
+R2 = prepare_sliced_quantity(R2)
+R3 = prepare_sliced_quantity(R3)
+
+#mask error metrics
+R0[R0<reg_err]=reg_err
+R1[R1<reg_err]=reg_err
+R2[R2<reg_err]=reg_err
+R3[R3<reg_err]=reg_err
+
+# plot everything
+
+from matplotlib.pyplot import imsave
+from matplotlib.colors import LogNorm
+import matplotlib.pyplot as plt
+from matplotlib import ticker
+
+new_x = np.linspace(0.0, data.metadata.boxsize.value[0], dimx) #np.linspace(np.min(x), np.max(x), dimx)
+new_y = np.linspace(0.0, data.metadata.boxsize.value[1], dimy) #np.linspace(np.min(y), np.max(y), dimy)
+
+def make_color_levels(cmin, cmax, c_res=10, log_sc=True):
+ cmax += divreg
+ if log_sc:
+ levmin = int(np.floor(np.log10(cmin)))
+ levmax = int(np.ceil(np.log10(cmax)))
+ levels = []
+ levels_short = []
+ for i in range(levmin, levmax):
+ levels_short += [10 ** i]
+ for j in range(int(c_res / 10), c_res):
+ levels += [(10 / c_res * j) * 10 ** i]
+
+ else:
+ levels = [cmin + (cmax - cmin) / c_res * k for k in range(c_res)]
+ levels_short = [cmin + (cmax - cmin) / c_res * k for k in range(0, c_res, 10)]
+ return levels, levels_short
+
+def make_density_plot(
+ Q, cmin, cmax, i, j, Q_name, c_res=10, log_sc=True, cmap="viridis"
+):
+ levels, levels_short = make_color_levels(cmin, cmax, c_res, log_sc)
+ if log_sc:
+ to_plot = ax[j].contourf(
+ new_x,
+ new_y,
+ Q.transpose(),
+ levels=np.array(levels),
+ locator=ticker.LogLocator(),
+ cmap=cmap,
+ extend="both",
+ )
+ else:
+ to_plot = ax[j].contourf(
+ new_x, new_y, Q.transpose(), levels=np.array(levels), cmap=cmap,
+ extend="both",
+ )
+
+ fig.colorbar(to_plot, ticks=levels_short)
+ ax[j].set_ylabel(Q_name)
+ ax[j].set_xlim(min(new_x),max(new_x))
+ ax[j].set_ylim(min(new_y),max(new_y))
+ return 0
+
+#Plot magnetic fields
+if to_plot=='A':
+ Ax = A[:,0]/Arms
+ Ay = A[:,1]/Arms
+ Az = A[:,2]/Arms
+
+ fig, ax = plt.subplots(1,3, sharex=True, figsize=(6*3, 5))
+ make_density_plot(
+ Ax.reshape((dimx, dimy)),
+ -1.0,
+ 1.0,
+ 0,
+ 0,
+ '$A_x$/$A_{rms}$',
+ c_res=cpts,
+ log_sc=False,
+ cmap=prefered_color
+ )
+ make_density_plot(
+ Ay.reshape((dimx, dimy)),
+ -1.0,
+ 1.0,
+ 0,
+ 1,
+ '$A_y$/$A_{rms}$',
+ c_res=cpts,
+ log_sc=False,
+ cmap=prefered_color
+ )
+ make_density_plot(
+ Az.reshape((dimx, dimy)),
+ -1.0,
+ 1.0,
+ 0,
+ 2,
+ '$A_z$/$A_{rms}$',
+ c_res=cpts,
+ log_sc=False,
+ cmap=prefered_color
+ )
+ ax[0].streamplot(new_x, new_y, np.transpose(Ax.reshape((dimx, dimy))), np.transpose(Ay.reshape((dimx, dimy))), color='w', density=2.0, linewidth=0.5, arrowsize=0.8)
+ ax[1].streamplot(new_x, new_y, np.transpose(Ax.reshape((dimx, dimy))), np.transpose(Ay.reshape((dimx, dimy))), color='w', density=2.0, linewidth=0.5, arrowsize=0.8)
+ ax[2].streamplot(new_x, new_y, np.transpose(Ax.reshape((dimx, dimy))), np.transpose(Ay.reshape((dimx, dimy))), color='w', density=2.0, linewidth=0.5, arrowsize=0.8)
+
+if to_plot=='B':
+ Bx = B[:,0]/Brms
+ By = B[:,1]/Brms
+ Bz = B[:,2]/Brms
+
+ fig, ax = plt.subplots(1,3, sharex=True, figsize=(6*3, 5))
+ make_density_plot(
+ Bx.reshape((dimx, dimy)),
+ -1.0,
+ 1.0,
+ 0,
+ 0,
+ '$B_x$/$B_{rms}$',
+ c_res=cpts,
+ log_sc=False,
+ cmap=prefered_color
+ )
+ make_density_plot(
+ By.reshape((dimx, dimy)),
+ -1.0,
+ 1.0,
+ 0,
+ 1,
+ '$B_y$/$B_{rms}$',
+ c_res=cpts,
+ log_sc=False,
+ cmap=prefered_color
+ )
+ make_density_plot(
+ Bz.reshape((dimx, dimy)),
+ -1.0,
+ 1.0,
+ 0,
+ 2,
+ '$B_z$/$B_{rms}$',
+ c_res=cpts,
+ log_sc=False,
+ cmap=prefered_color
+ )
+ ax[0].streamplot(new_x, new_y, np.transpose(Bx.reshape((dimx, dimy))), np.transpose(By.reshape((dimx, dimy))), color='w', density=2.0, linewidth=0.5, arrowsize=0.8)
+ ax[1].streamplot(new_x, new_y, np.transpose(Bx.reshape((dimx, dimy))), np.transpose(By.reshape((dimx, dimy))), color='w', density=2.0, linewidth=0.5, arrowsize=0.8)
+ ax[2].streamplot(new_x, new_y, np.transpose(Bx.reshape((dimx, dimy))), np.transpose(By.reshape((dimx, dimy))), color='w', density=2.0, linewidth=0.5, arrowsize=0.8)
+
+if to_plot=='errors':
+
+ fig, ax = plt.subplots(1,4, sharex=True, figsize=(24, 5)) #for 3 plts 18 for 4 plts use 24
+
+ make_density_plot(
+ R0.reshape((dimx, dimy)),
+ reg_err,
+ 1.0,
+ 0,
+ 0,
+ "$R_0$",
+ c_res=cpts,
+ #log_sc=False,
+ cmap=prefered_color
+ )
+ make_density_plot(
+ R1.reshape((dimx, dimy)),
+ reg_err,
+ 1.0,
+ 0,
+ 1,
+ "$R_1$",
+ c_res=cpts,
+ #log_sc=False,
+ cmap=prefered_color
+ )
+ make_density_plot(
+ R2.reshape((dimx, dimy)),
+ reg_err,
+ 1.0,
+ 0,
+ 2,
+ "$R_2$",
+ c_res=cpts,
+ #log_sc=False,
+ cmap=prefered_color
+ )
+ make_density_plot(
+ R3.reshape((dimx, dimy)),
+ reg_err,
+ 1.0,
+ 0,
+ 3,
+ "$R_3$",
+ c_res=cpts,
+ log_sc=False,
+ cmap=prefered_color
+ )
+
+ax[0].set_title(f"Nneigh={int(neighbours[0]):}, Npart={len(data.gas.coordinates):}")
+ax[1].set_title(f"t={t:.2e}")
+ax[2].set_title("$z_{slice}/L_{box}$="+slice_height)
+fig.tight_layout()
+
+plt.savefig(sys.argv[2], dpi=100)
diff --git a/src/mhd/DInduction/mhd_io.h b/src/mhd/DInduction/mhd_io.h
index f3a0911241fafcd11c78a5d2bca9b6c9db1ff5db..c5c69e7748187bd3f66911209c162f0bc8d23f49 100644
--- a/src/mhd/DInduction/mhd_io.h
+++ b/src/mhd/DInduction/mhd_io.h
@@ -45,6 +45,92 @@ INLINE static void convert_B(const struct engine* e, const struct part* p,
ret[2] = p->mhd_data.BPred[2];
}
+
+/* Calculation of error metrics*/
+
+INLINE static void calculate_R0(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+/* Calculating R0 error metric */
+ const float b = sqrt(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]);
+ ret[0] = fabs(p->mhd_data.divB * p->h / (b+1.e-30));
+
+/* Discarding noise*/
+ const float abs_SPH_divB_err = fabs(p->mhd_data.BPred[0] * p->mhd_data.mean_grad_SPH_err[0] +
+ p->mhd_data.BPred[1] * p->mhd_data.mean_grad_SPH_err[1] +
+ p->mhd_data.BPred[2] * p->mhd_data.mean_grad_SPH_err[2]);
+ if (fabs(p->mhd_data.divB)<10*abs_SPH_divB_err){
+ ret[0] = 0.f;
+ }
+}
+
+INLINE static void calculate_R1(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+/* Calculating R1 error metric */
+ const float b = sqrt(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 fmag = sqrt(p->mhd_data.tot_mag_F[0] * p->mhd_data.tot_mag_F[0] +
+ p->mhd_data.tot_mag_F[1] * p->mhd_data.tot_mag_F[1] +
+ p->mhd_data.tot_mag_F[2] * p->mhd_data.tot_mag_F[2]);
+ const float fmag_dot_b = p->mhd_data.BPred[0] * p->mhd_data.tot_mag_F[0] +
+ p->mhd_data.BPred[1] * p->mhd_data.tot_mag_F[1] +
+ p->mhd_data.BPred[2] * p->mhd_data.tot_mag_F[2];
+ ret[0] = fabs(fmag_dot_b / ( b * fmag + 1.e-30));
+
+/* Discarding noise */
+ const float ftot = sqrt(p->mhd_data.tot_F[0] * p->mhd_data.tot_F[0] +
+ p->mhd_data.tot_F[1] * p->mhd_data.tot_F[1] +
+ p->mhd_data.tot_F[2] * p->mhd_data.tot_F[2]);
+ const float RF = fmag / (ftot+fmag+1.e-30);
+ const float mag_force_magnitude_err = b*b/(p->rho * 1+1.e-30);
+
+ const float abs_SPH_Fmag_err = sqrt((fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[0])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[0]))*(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[0])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[0]))+(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[1])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[1]))*(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[1])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[1]))+(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[2])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[2]))*(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[2])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[2])));
+
+ if (fmag<10*abs_SPH_Fmag_err || RF<0.1){
+ ret[0] = 0.f;
+ }
+}
+
+INLINE static void calculate_R2(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+/* Calculating R2 error metric */
+ const float cb = sqrt(p->mhd_data.curlB[0] * p->mhd_data.curlB[0] +
+ p->mhd_data.curlB[1] * p->mhd_data.curlB[1] +
+ p->mhd_data.curlB[2] * p->mhd_data.curlB[2]);
+ ret[0] = fabs(p->mhd_data.divB / (cb+1.e-30));
+
+/* Discarding noise*/
+ const float abs_SPH_divB_err = fabs(p->mhd_data.BPred[0] * p->mhd_data.mean_grad_SPH_err[0] +
+ p->mhd_data.BPred[1] * p->mhd_data.mean_grad_SPH_err[1] +
+ p->mhd_data.BPred[2] * p->mhd_data.mean_grad_SPH_err[2]);
+ const float abs_SPH_curlB_err = sqrt((p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[2]-p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[1])*(p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[2]-p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[1])+(p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[0]-p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[2])*(p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[0]-p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[2])+(p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[1]-p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[0])*(p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[1]-p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[0]));
+
+ if (fabs(p->mhd_data.divB)<10*abs_SPH_divB_err || cb<10*abs_SPH_curlB_err){
+ ret[0] = 0.f;
+ }
+}
+
+INLINE static void calculate_R3(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+/* Calculating R3 error metric */
+ const float b = sqrt(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 cb = sqrt(p->mhd_data.curlB[0] * p->mhd_data.curlB[0] +
+ p->mhd_data.curlB[1] * p->mhd_data.curlB[1] +
+ p->mhd_data.curlB[2] * p->mhd_data.curlB[2]);
+ ret[0] = cb * p->h / (b+1.e-30);
+
+/* Discarding noise*/
+ const float abs_SPH_curlB_err = sqrt((p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[2]-p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[1])*(p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[2]-p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[1])+(p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[0]-p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[2])*(p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[0]-p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[2])+(p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[1]-p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[0])*(p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[1]-p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[0]));
+
+ if (cb<10*abs_SPH_curlB_err){
+ ret[0] = 0.f;
+ }
+}
+
/**
* @brief Specifies which particle fields to write to a dataset
*
@@ -71,33 +157,28 @@ INLINE static int mhd_write_particles(const struct part* parts,
UNIT_CONV_ELECTRIC_CHARGE_FIELD_STRENGTH,
mhd_comoving_factor + 1.f, parts, mhd_data.phi,
"Dedner scalars associated to the particles");
+
+ /* Error metrics */
+ list[3] = io_make_output_field_convert_part(
+ "R0", FLOAT, 1, 1,
+ 1, parts, xparts, calculate_R0,
+ "Classical error metric, indicates places with large divergence");
+ list[4] = io_make_output_field_convert_part(
+ "R1", FLOAT, 1, 1,
+ 1, parts, xparts, calculate_R1,
+ "Error metric, angle between B field and total Fmag. Indicates unpysical magnetic force");
+ list[5] = io_make_output_field_convert_part(
+ "R2", FLOAT, 1, 1,
+ 1, parts, xparts, calculate_R2,
+ "Error metric, ratio of divB and |curlB|. Estimates upper limit on B_monopole/B_physical");
+ list[6] = io_make_output_field_convert_part(
+ "R3", FLOAT, 1, 1,
+ 1, parts, xparts, calculate_R3,
+ "Error metric, shows relation of smoothing length to characteristic B gradient scale");
- /* Quantities for error monitoring. Check dimensionality and cosmo factors*/
- list[3] = io_make_output_field(
- "MagneticFluxCurl", FLOAT, 3, 1,
- 1, parts, mhd_data.curlB,
- "The curl of Magnetic flux densities of the particles");
-
- list[4] = io_make_output_field(
- "SPH_1", FLOAT, 1, 1,
- 1, parts, mhd_data.mean_SPH_err,
- "SPH smoothed 1 associated to a particle. Used to estimate SPH sum noise");
-
- list[5] = io_make_output_field(
- "SPH_grad_1", FLOAT, 3, 1,
- 1, parts, mhd_data.mean_grad_SPH_err,
- "SPH smoothed gradient of 1 associated to a particle. Used to estimate SPH sum noise of any gradient");
-
- list[6] = io_make_output_field(
- "Fmag", FLOAT, 3, 1,
- 1, parts, mhd_data.tot_mag_F,
- "Magnetic force acting on a particle with correcitons");
-
- list[6] = io_make_output_field(
- "Ftot", FLOAT, 3, 1,
- 1, parts, mhd_data.tot_F,
- "Total force acting on a particle");
return 7;
+
+
}
/**
diff --git a/src/mhd/DirectInduction/mhd_iact.h b/src/mhd/DirectInduction/mhd_iact.h
index d2d40c736bb1c78157ebf16328341a03be41cc1a..b06aa2dc04ff0d5313f83d16f39582c989975af5 100644
--- a/src/mhd/DirectInduction/mhd_iact.h
+++ b/src/mhd/DirectInduction/mhd_iact.h
@@ -478,8 +478,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
/* Save forces */
for (int k=1;k<3;k++){
- pi->mhd_data.tot_mag_F[k] = - mj * sph_acc_term_i[k];
- pj->mhd_data.tot_mag_F[k] = - mi * sph_acc_term_j[k];
+ pi->mhd_data.tot_mag_F[k] -= mj * sph_acc_term_i[k];
+ pj->mhd_data.tot_mag_F[k] -= mi * sph_acc_term_j[k];
pi->mhd_data.tot_F[k] = pi->a_hydro[k];
pj->mhd_data.tot_F[k] = pj->a_hydro[k];
}
@@ -815,7 +815,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
/* Save forces */
for (int k=1;k<3;k++){
- pi->mhd_data.tot_mag_F[k] = - mj * sph_acc_term_i[k];
+ pi->mhd_data.tot_mag_F[k] -= mj * sph_acc_term_i[k];
pi->mhd_data.tot_F[k] = pi->a_hydro[k];
}
diff --git a/src/mhd/DirectInduction/mhd_io.h b/src/mhd/DirectInduction/mhd_io.h
index 5b7d4fc67db12c6211a80d8623d2a3a1dd164ac8..60c150985246dfc200714a842454d9f40c4144ee 100644
--- a/src/mhd/DirectInduction/mhd_io.h
+++ b/src/mhd/DirectInduction/mhd_io.h
@@ -46,6 +46,93 @@ INLINE static void convert_B(const struct engine* e, const struct part* p,
ret[2] = xp->mhd_data.B_over_rho_full[2] * p->rho;
}
+INLINE static void calculate_R0(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+/* Calculating R0 error metric */
+ const float b = sqrt(xp->mhd_data.B_over_rho_full[0] * xp->mhd_data.B_over_rho_full[0] +
+ xp->mhd_data.B_over_rho_full[1] * xp->mhd_data.B_over_rho_full[1] +
+ xp->mhd_data.B_over_rho_full[2] * xp->mhd_data.B_over_rho_full[2]) * p->rho;
+ ret[0] = fabs(p->mhd_data.divB * p->h / (b+1.e-30));
+
+/* Discarding noise*/
+ const float abs_SPH_divB_err = fabs(xp->mhd_data.B_over_rho_full[0] * p->mhd_data.mean_grad_SPH_err[0] +
+ xp->mhd_data.B_over_rho_full[1] * p->mhd_data.mean_grad_SPH_err[1] +
+ xp->mhd_data.B_over_rho_full[2] * p->mhd_data.mean_grad_SPH_err[2]) * p->rho;
+ if (fabs(p->mhd_data.divB)<10*abs_SPH_divB_err){
+ ret[0] = 0.f;
+ }
+}
+
+
+INLINE static void calculate_R1(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+/* Calculating R1 error metric */
+ const float b = sqrt(xp->mhd_data.B_over_rho_full[0] * xp->mhd_data.B_over_rho_full[0] +
+ xp->mhd_data.B_over_rho_full[1] * xp->mhd_data.B_over_rho_full[1] +
+ xp->mhd_data.B_over_rho_full[2] * xp->mhd_data.B_over_rho_full[2]) * p->rho;
+ const float fmag = sqrt(p->mhd_data.tot_mag_F[0] * p->mhd_data.tot_mag_F[0] +
+ p->mhd_data.tot_mag_F[1] * p->mhd_data.tot_mag_F[1] +
+ p->mhd_data.tot_mag_F[2] * p->mhd_data.tot_mag_F[2]);
+ const float fmag_dot_b = (xp->mhd_data.B_over_rho_full[0] * p->mhd_data.tot_mag_F[0] +
+ xp->mhd_data.B_over_rho_full[1] * p->mhd_data.tot_mag_F[1] +
+ xp->mhd_data.B_over_rho_full[2] * p->mhd_data.tot_mag_F[2]) * p->rho;
+ ret[0] = fabs(fmag_dot_b / ( b * fmag + 1.e-30));
+
+/* Discarding noise */
+ const float ftot = sqrt(p->mhd_data.tot_F[0] * p->mhd_data.tot_F[0] +
+ p->mhd_data.tot_F[1] * p->mhd_data.tot_F[1] +
+ p->mhd_data.tot_F[2] * p->mhd_data.tot_F[2]);
+ const float RF = fmag / (ftot+fmag+1.e-30);
+ const float mag_force_magnitude_err = b*b/(p->rho * 1+1.e-30);
+
+ const float abs_SPH_Fmag_err = sqrt((fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[0])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[0]))*(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[0])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[0]))+(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[1])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[1]))*(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[1])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[1]))+(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[2])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[2]))*(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[2])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[2])));
+
+ if (fmag<10*abs_SPH_Fmag_err || RF<0.1){
+ ret[0] = 0.f;
+ }
+}
+
+
+INLINE static void calculate_R2(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+/* Calculating R2 error metric */
+ const float cb = sqrt(p->mhd_data.curl_B[0] * p->mhd_data.curl_B[0] +
+ p->mhd_data.curl_B[1] * p->mhd_data.curl_B[1] +
+ p->mhd_data.curl_B[2] * p->mhd_data.curl_B[2]);
+ ret[0] = fabs(p->mhd_data.divB / (cb+1.e-30));
+
+/* Discarding noise*/
+ const float abs_SPH_divB_err = fabs(xp->mhd_data.B_over_rho_full[0] * p->mhd_data.mean_grad_SPH_err[0] +
+ xp->mhd_data.B_over_rho_full[1] * p->mhd_data.mean_grad_SPH_err[1] +
+ xp->mhd_data.B_over_rho_full[2] * p->mhd_data.mean_grad_SPH_err[2]) * p->rho;
+
+ const float abs_SPH_curlB_err = sqrt((p->mhd_data.mean_grad_SPH_err[1] * xp->mhd_data.B_over_rho_full[2]-p->mhd_data.mean_grad_SPH_err[2] * xp->mhd_data.B_over_rho_full[1])*(p->mhd_data.mean_grad_SPH_err[1] * xp->mhd_data.B_over_rho_full[2]-p->mhd_data.mean_grad_SPH_err[2] * xp->mhd_data.B_over_rho_full[1])+(p->mhd_data.mean_grad_SPH_err[2] * xp->mhd_data.B_over_rho_full[0]-p->mhd_data.mean_grad_SPH_err[0] * xp->mhd_data.B_over_rho_full[2])*(p->mhd_data.mean_grad_SPH_err[2] * xp->mhd_data.B_over_rho_full[0]-p->mhd_data.mean_grad_SPH_err[0] * xp->mhd_data.B_over_rho_full[2])+(p->mhd_data.mean_grad_SPH_err[0] * xp->mhd_data.B_over_rho_full[1]-p->mhd_data.mean_grad_SPH_err[1] * xp->mhd_data.B_over_rho_full[0])*(p->mhd_data.mean_grad_SPH_err[0] * xp->mhd_data.B_over_rho_full[1]-p->mhd_data.mean_grad_SPH_err[1] * xp->mhd_data.B_over_rho_full[0])) * p->rho;
+
+ if (fabs(p->mhd_data.divB)<10*abs_SPH_divB_err || cb<10*abs_SPH_curlB_err){
+ ret[0] = 0.f;
+ }
+}
+
+
+INLINE static void calculate_R3(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+/* Calculating R3 error metric */
+ const float b = sqrt(xp->mhd_data.B_over_rho_full[0] * xp->mhd_data.B_over_rho_full[0] +
+ xp->mhd_data.B_over_rho_full[1] * xp->mhd_data.B_over_rho_full[1] +
+ xp->mhd_data.B_over_rho_full[2] * xp->mhd_data.B_over_rho_full[2]) * p->rho;
+ const float cb = sqrt(p->mhd_data.curl_B[0] * p->mhd_data.curl_B[0] +
+ p->mhd_data.curl_B[1] * p->mhd_data.curl_B[1] +
+ p->mhd_data.curl_B[2] * p->mhd_data.curl_B[2]);
+ ret[0] = cb * p->h / (b+1.e-30);
+
+/* Discarding noise*/
+ const float abs_SPH_curlB_err = sqrt((p->mhd_data.mean_grad_SPH_err[1] * xp->mhd_data.B_over_rho_full[2]-p->mhd_data.mean_grad_SPH_err[2] * xp->mhd_data.B_over_rho_full[1])*(p->mhd_data.mean_grad_SPH_err[1] * xp->mhd_data.B_over_rho_full[2]-p->mhd_data.mean_grad_SPH_err[2] * xp->mhd_data.B_over_rho_full[1])+(p->mhd_data.mean_grad_SPH_err[2] * xp->mhd_data.B_over_rho_full[0]-p->mhd_data.mean_grad_SPH_err[0] * xp->mhd_data.B_over_rho_full[2])*(p->mhd_data.mean_grad_SPH_err[2] * xp->mhd_data.B_over_rho_full[0]-p->mhd_data.mean_grad_SPH_err[0] * xp->mhd_data.B_over_rho_full[2])+(p->mhd_data.mean_grad_SPH_err[0] * xp->mhd_data.B_over_rho_full[1]-p->mhd_data.mean_grad_SPH_err[1] * xp->mhd_data.B_over_rho_full[0])*(p->mhd_data.mean_grad_SPH_err[0] * xp->mhd_data.B_over_rho_full[1]-p->mhd_data.mean_grad_SPH_err[1] * xp->mhd_data.B_over_rho_full[0])) * p->rho;
+
+ if (cb<10*abs_SPH_curlB_err){
+ ret[0] = 0.f;
+ }
+}
+
/**
* @brief Specifies which particle fields to write to a dataset
*
@@ -93,27 +180,23 @@ INLINE static int mhd_write_particles(const struct part* parts,
"AlphaAR", FLOAT, 1, UNIT_CONV_NO_UNITS, 1.f, parts, mhd_data.alpha_AR,
"Artificial resistivity switch of the particles");
- /* Quantities for error monitoring. Check dimensionality and cosmo factors*/
-
- list[7] = io_make_output_field(
- "SPH_1", FLOAT, 1, 1,
- 1, parts, mhd_data.mean_SPH_err,
- "SPH smoothed 1 associated to a particle. Used to estimate SPH sum noise");
-
- list[8] = io_make_output_field(
- "SPH_grad_1", FLOAT, 3, 1,
- 1, parts, mhd_data.mean_grad_SPH_err,
- "SPH smoothed gradient of 1 associated to a particle. Used to estimate SPH sum noise of any gradient");
-
- list[9] = io_make_output_field(
- "Fmag", FLOAT, 3, 1,
- 1, parts, mhd_data.tot_mag_F,
- "Magnetic force acting on a particle with correcitons");
-
- list[10] = io_make_output_field(
- "Ftot", FLOAT, 3, 1,
- 1, parts, mhd_data.tot_F,
- "Total force acting on a particle");
+/* Error metrics */
+ list[7] = io_make_output_field_convert_part(
+ "R0", FLOAT, 1, 1,
+ 1, parts, xparts, calculate_R0,
+ "Classical error metric, indicates places with large divergence");
+ list[8] = io_make_output_field_convert_part(
+ "R1", FLOAT, 1, 1,
+ 1, parts, xparts, calculate_R1,
+ "Error metric, angle between B field and total Fmag. Indicates unpysical magnetic force");
+ list[9] = io_make_output_field_convert_part(
+ "R2", FLOAT, 1, 1,
+ 1, parts, xparts, calculate_R2,
+ "Error metric, ratio of divB and |curlB|. Estimates upper limit on B_monopole/B_physical");
+ list[10] = io_make_output_field_convert_part(
+ "R3", FLOAT, 1, 1,
+ 1, parts, xparts, calculate_R3,
+ "Error metric, shows relation of smoothing length to characteristic B gradient scale");
return 11;
}
diff --git a/src/mhd/VPotential/mhd_io.h b/src/mhd/VPotential/mhd_io.h
index 6947dfb9254227492a3f568b11ccf0a2510477d2..005ceaa2f0b2d3d289b1518aace082d88b64ae90 100644
--- a/src/mhd/VPotential/mhd_io.h
+++ b/src/mhd/VPotential/mhd_io.h
@@ -48,6 +48,94 @@ INLINE static void convert_B(const struct engine* e, const struct part* p,
ret[2] = p->mhd_data.BPred[2];
}
+
+INLINE static void calculate_R0(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+/* Calculating R0 error metric */
+ const float b = sqrt(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]);
+ ret[0] = fabs(p->mhd_data.divB * p->h / (b+1.e-30));
+
+/* Discarding noise*/
+ const float abs_SPH_divB_err = fabs(p->mhd_data.BPred[0] * p->mhd_data.mean_grad_SPH_err[0] +
+ p->mhd_data.BPred[1] * p->mhd_data.mean_grad_SPH_err[1] +
+ p->mhd_data.BPred[2] * p->mhd_data.mean_grad_SPH_err[2]);
+ if (fabs(p->mhd_data.divB)<10*abs_SPH_divB_err){
+ ret[0] = 0.f;
+ }
+}
+
+
+INLINE static void calculate_R1(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+/* Calculating R1 error metric */
+ const float b = sqrt(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 fmag = sqrt(p->mhd_data.tot_mag_F[0] * p->mhd_data.tot_mag_F[0] +
+ p->mhd_data.tot_mag_F[1] * p->mhd_data.tot_mag_F[1] +
+ p->mhd_data.tot_mag_F[2] * p->mhd_data.tot_mag_F[2]);
+ const float fmag_dot_b = p->mhd_data.BPred[0] * p->mhd_data.tot_mag_F[0] +
+ p->mhd_data.BPred[1] * p->mhd_data.tot_mag_F[1] +
+ p->mhd_data.BPred[2] * p->mhd_data.tot_mag_F[2];
+ ret[0] = fabs(fmag_dot_b / ( b * fmag + 1.e-30));
+
+/* Discarding noise */
+ const float ftot = sqrt(p->mhd_data.tot_F[0] * p->mhd_data.tot_F[0] +
+ p->mhd_data.tot_F[1] * p->mhd_data.tot_F[1] +
+ p->mhd_data.tot_F[2] * p->mhd_data.tot_F[2]);
+ const float RF = fmag / (ftot+fmag+1.e-30);
+ const float mag_force_magnitude_err = b*b/(p->rho * 1+1.e-30);
+
+ const float abs_SPH_Fmag_err = sqrt((fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[0])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[0]))*(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[0])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[0]))+(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[1])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[1]))*(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[1])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[1]))+(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[2])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[2]))*(fabs(mag_force_magnitude_err * p->mhd_data.mean_grad_SPH_err[2])+fabs(p->mhd_data.mean_SPH_err * p->mhd_data.tot_mag_F[2])));
+
+ if (fmag<10*abs_SPH_Fmag_err || RF<0.1){
+ ret[0] = 0.f;
+ }
+}
+
+
+INLINE static void calculate_R2(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+/* Calculating R2 error metric */
+ const float cb = sqrt(p->mhd_data.curlB[0] * p->mhd_data.curlB[0] +
+ p->mhd_data.curlB[1] * p->mhd_data.curlB[1] +
+ p->mhd_data.curlB[2] * p->mhd_data.curlB[2]);
+ ret[0] = fabs(p->mhd_data.divB / (cb+1.e-30));
+
+/* Discarding noise*/
+ const float abs_SPH_divB_err = fabs(p->mhd_data.BPred[0] * p->mhd_data.mean_grad_SPH_err[0] +
+ p->mhd_data.BPred[1] * p->mhd_data.mean_grad_SPH_err[1] +
+ p->mhd_data.BPred[2] * p->mhd_data.mean_grad_SPH_err[2]);
+ const float abs_SPH_curlB_err = sqrt((p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[2]-p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[1])*(p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[2]-p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[1])+(p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[0]-p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[2])*(p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[0]-p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[2])+(p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[1]-p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[0])*(p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[1]-p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[0]));
+
+ if (fabs(p->mhd_data.divB)<10*abs_SPH_divB_err || cb<10*abs_SPH_curlB_err){
+ ret[0] = 0.f;
+ }
+}
+
+
+INLINE static void calculate_R3(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+/* Calculating R3 error metric */
+ const float b = sqrt(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 cb = sqrt(p->mhd_data.curlB[0] * p->mhd_data.curlB[0] +
+ p->mhd_data.curlB[1] * p->mhd_data.curlB[1] +
+ p->mhd_data.curlB[2] * p->mhd_data.curlB[2]);
+ ret[0] = cb * p->h / (b+1.e-30);
+
+/* Discarding noise*/
+ const float abs_SPH_curlB_err = sqrt((p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[2]-p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[1])*(p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[2]-p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[1])+(p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[0]-p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[2])*(p->mhd_data.mean_grad_SPH_err[2] * p->mhd_data.BPred[0]-p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[2])+(p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[1]-p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[0])*(p->mhd_data.mean_grad_SPH_err[0] * p->mhd_data.BPred[1]-p->mhd_data.mean_grad_SPH_err[1] * p->mhd_data.BPred[0]));
+
+ if (cb<10*abs_SPH_curlB_err){
+ ret[0] = 0.f;
+ }
+}
+
+
/**
* @brief Specifies which particle fields to write to a dataset
*
@@ -87,33 +175,25 @@ INLINE static int mhd_write_particles(const struct part* parts,
mhd_comoving_factor, parts, mhd_data.divA,
"Co-moving vector potential divergences of particles");
- /* Quantities for error monitoring. Check dimensionality and cosmo factors*/
- list[5] = io_make_output_field(
- "MagneticFluxCurl", FLOAT, 3, 1,
- 1, parts, mhd_data.curlB,
- "The curl of Magnetic flux densities of the particles");
-
- list[6] = io_make_output_field(
- "SPH_1", FLOAT, 1, 1,
- 1, parts, mhd_data.mean_SPH_err,
- "SPH smoothed 1 associated to a particle. Used to estimate SPH sum noise");
-
- list[7] = io_make_output_field(
- "SPH_grad_1", FLOAT, 3, 1,
- 1, parts, mhd_data.mean_grad_SPH_err,
- "SPH smoothed gradient of 1 associated to a particle. Used to estimate SPH sum noise of any gradient");
-
- list[8] = io_make_output_field(
- "Fmag", FLOAT, 3, 1,
- 1, parts, mhd_data.tot_mag_F,
- "Magnetic force acting on a particle with correcitons");
-
- list[9] = io_make_output_field(
- "Ftot", FLOAT, 3, 1,
- 1, parts, mhd_data.tot_F,
- "Total force acting on a particle");
-
- return 10;
+/* Error metrics */
+ list[5] = io_make_output_field_convert_part(
+ "R0", FLOAT, 1, 1,
+ 1, parts, xparts, calculate_R0,
+ "Classical error metric, indicates places with large divergence");
+ list[6] = io_make_output_field_convert_part(
+ "R1", FLOAT, 1, 1,
+ 1, parts, xparts, calculate_R1,
+ "Error metric, angle between B field and total Fmag. Indicates unpysical magnetic force");
+ list[7] = io_make_output_field_convert_part(
+ "R2", FLOAT, 1, 1,
+ 1, parts, xparts, calculate_R2,
+ "Error metric, ratio of divB and |curlB|. Estimates upper limit on B_monopole/B_physical");
+ list[8] = io_make_output_field_convert_part(
+ "R3", FLOAT, 1, 1,
+ 1, parts, xparts, calculate_R3,
+ "Error metric, shows relation of smoothing length to characteristic B gradient scale");
+
+ return 9;
}
/**