From 24d16ddc9ef1c873beed2e599a19e3088751a729 Mon Sep 17 00:00:00 2001
From: Federico Stasyszyn <fstasyszyn@unc.edu.ar>
Date: Wed, 29 Nov 2023 15:39:08 +0100
Subject: [PATCH] formating script
---
.../HydroTests/SedovBlast_3D/plotSolution.py | 10 +-
examples/MHDTests/CosmoVolume/makeIC.py | 68 +++---
examples/MHDTests/CosmoVolume/plotRhoB.py | 14 +-
examples/MHDTests/CosmoVolume/plot_all.py | 168 +++++++++------
examples/MHDTests/FastRotor/makeIC.py | 60 +++---
examples/MHDTests/FastRotor/plot_all.py | 199 ++++++++++++------
examples/MHDTests/MagneticBlastWave/makeIC.py | 64 +++---
.../MHDTests/MagneticBlastWave/plot_all.py | 180 ++++++++++------
src/mhd/DInduction/mhd.h | 7 +-
src/mhd/DirectInduction/mhd.h | 8 +-
src/mhd/DirectInduction/mhd_iact.h | 4 +-
src/mhd/VPotential/mhd.h | 72 ++++---
src/mhd/VPotential/mhd_iact.h | 15 +-
src/mhd/VPotential/mhd_io.h | 11 +-
src/units.c | 2 +-
15 files changed, 518 insertions(+), 364 deletions(-)
diff --git a/examples/HydroTests/SedovBlast_3D/plotSolution.py b/examples/HydroTests/SedovBlast_3D/plotSolution.py
index 2dd850cb41..d5fb1712a8 100644
--- a/examples/HydroTests/SedovBlast_3D/plotSolution.py
+++ b/examples/HydroTests/SedovBlast_3D/plotSolution.py
@@ -393,13 +393,13 @@ tight_layout()
savefig("Sedov.png")
data = np.loadtxt("statistics.txt")
-t = data[:,1]
-E_kin = data[:,13]
-E_int = data[:,14]
+t = data[:, 1]
+E_kin = data[:, 13]
+E_int = data[:, 14]
E_tot = E_kin + E_int
print(E_tot)
figure()
plot(t, E_tot)
-plot(t, E_int, '--')
-plot(t, E_kin, ':')
+plot(t, E_int, "--")
+plot(t, E_kin, ":")
savefig("Energy.png")
diff --git a/examples/MHDTests/CosmoVolume/makeIC.py b/examples/MHDTests/CosmoVolume/makeIC.py
index 3297bcdcbd..b336e7f6e6 100644
--- a/examples/MHDTests/CosmoVolume/makeIC.py
+++ b/examples/MHDTests/CosmoVolume/makeIC.py
@@ -24,62 +24,62 @@ head = infile["/Header"]
pos_in = infile["/PartType0/Coordinates"][:, :]
BoxSize = head.attrs["BoxSize"]
-afact = head.attrs["Time"]
+afact = head.attrs["Time"]
N_in = head.attrs["NumPart_Total"][0]
-print("BoxSize [in Comoving]: ",BoxSize)
-print("Total Gas Part:",N_in)
-print("Expansion Factor:",afact)
+print("BoxSize [in Comoving]: ", BoxSize)
+print("Total Gas Part:", N_in)
+print("Expansion Factor:", afact)
infile.close()
#### physical B field in Gauss #####
-Bini = 1E-12
+Bini = 1e-12
-B = np.zeros((N_in,3))
-A = np.zeros((N_in,3))
+B = np.zeros((N_in, 3))
+A = np.zeros((N_in, 3))
wavelen = BoxSize / 10.0
wavenum = 2.0 * np.pi / wavelen
-Aini = Bini / wavenum * afact;
+Aini = Bini / wavenum * afact
-#print(wavelen,wavenum)
+# print(wavelen,wavenum)
#####################################################
-#### Positions in Comoving ####
+#### Positions in Comoving ####
#### Other quatities in Physical in the IC files ####
-A[:,0] = Aini * (np.sin(pos_in[:,2]*wavenum) + np.cos(pos_in[:,1]*wavenum))
-A[:,1] = Aini * (np.sin(pos_in[:,0]*wavenum) + np.cos(pos_in[:,2]*wavenum))
-A[:,2] = Aini * (np.sin(pos_in[:,1]*wavenum) + np.cos(pos_in[:,0]*wavenum))
-B[:,0] = Bini * (np.sin(pos_in[:,2]*wavenum) + np.cos(pos_in[:,1]*wavenum))
-B[:,1] = Bini * (np.sin(pos_in[:,0]*wavenum) + np.cos(pos_in[:,2]*wavenum))
-B[:,2] = Bini * (np.sin(pos_in[:,1]*wavenum) + np.cos(pos_in[:,0]*wavenum))
+A[:, 0] = Aini * (np.sin(pos_in[:, 2] * wavenum) + np.cos(pos_in[:, 1] * wavenum))
+A[:, 1] = Aini * (np.sin(pos_in[:, 0] * wavenum) + np.cos(pos_in[:, 2] * wavenum))
+A[:, 2] = Aini * (np.sin(pos_in[:, 1] * wavenum) + np.cos(pos_in[:, 0] * wavenum))
+B[:, 0] = Bini * (np.sin(pos_in[:, 2] * wavenum) + np.cos(pos_in[:, 1] * wavenum))
+B[:, 1] = Bini * (np.sin(pos_in[:, 0] * wavenum) + np.cos(pos_in[:, 2] * wavenum))
+B[:, 2] = Bini * (np.sin(pos_in[:, 1] * wavenum) + np.cos(pos_in[:, 0] * wavenum))
-print("VPotencial min/max: ",min(A[:,0]),max(A[:,0]))
-print("Bfield min/max",min(B[:,0]),max(B[:,0]))
+print("VPotencial min/max: ", min(A[:, 0]), max(A[:, 0]))
+print("Bfield min/max", min(B[:, 0]), max(B[:, 0]))
#### Constant field in X direction ####
-#B[:,0] = np.where(pos_in[:,1] < BoxSize * 0.5, Bini, -Bini)
-#B[:,1] = 0.0
-#B[:,2] = 0.0
-#A[:,0] = 0.0
-#A[:,1] = 0.0
-#A[:,2] = np.where(pos_in[:,1] < BoxSize *0.5, Bini * pos_in[:,1], Bini * (BoxSize-pos_in[:,1]))
+# B[:,0] = np.where(pos_in[:,1] < BoxSize * 0.5, Bini, -Bini)
+# B[:,1] = 0.0
+# B[:,2] = 0.0
+# A[:,0] = 0.0
+# A[:,1] = 0.0
+# A[:,2] = np.where(pos_in[:,1] < BoxSize *0.5, Bini * pos_in[:,1], Bini * (BoxSize-pos_in[:,1]))
-#print(min(A[:,2]),max(A[:,2]))
-#print(min(B[:,0]),max(B[:,0]))
+# print(min(A[:,2]),max(A[:,2]))
+# print(min(B[:,0]),max(B[:,0]))
-os.system("cp "+fileInputName+" "+fileOutputName)
+os.system("cp " + fileInputName + " " + fileOutputName)
# File
fileOutput = h5py.File(fileOutputName, "a")
# Units
-#print(list(fileOutput.keys()))
-#grpu = fileOutput.require_group(["/Units"])
+# print(list(fileOutput.keys()))
+# grpu = fileOutput.require_group(["/Units"])
grpu = fileOutput["/Units"]
-#print(grpu.attrs["Unit length in cgs (U_L)"])
-#print(grpu.attrs["Unit mass in cgs (U_M)"])
-#print(grpu.attrs["Unit time in cgs (U_t)"])
-#print(grpu.attrs["Unit current in cgs (U_I)"])
-#print(grpu.attrs["Unit temperature in cgs (U_T)"])
+# print(grpu.attrs["Unit length in cgs (U_L)"])
+# print(grpu.attrs["Unit mass in cgs (U_M)"])
+# print(grpu.attrs["Unit time in cgs (U_t)"])
+# print(grpu.attrs["Unit current in cgs (U_I)"])
+# print(grpu.attrs["Unit temperature in cgs (U_T)"])
grpu.attrs["Unit current in cgs (U_I)"] = 4.788e6 # amperes to be Gauss
# Particle group
diff --git a/examples/MHDTests/CosmoVolume/plotRhoB.py b/examples/MHDTests/CosmoVolume/plotRhoB.py
index 67742724b5..a8539732f9 100644
--- a/examples/MHDTests/CosmoVolume/plotRhoB.py
+++ b/examples/MHDTests/CosmoVolume/plotRhoB.py
@@ -31,7 +31,7 @@ import numpy as np
import h5py
import sys
-#plt.style.use("../../../../tools/stylesheets/mnras.mplstyle")
+# plt.style.use("../../../../tools/stylesheets/mnras.mplstyle")
snap = int(sys.argv[1])
@@ -83,14 +83,14 @@ log_B_max = -6
bins_x = np.linspace(log_rho_min, log_rho_max, 54)
bins_y = np.linspace(log_B_min, log_B_max, 54)
-#H, _, _ = np.histogram2d(log_rho, log_B, bins=[bins_x, bins_y], normed=True)
+# H, _, _ = np.histogram2d(log_rho, log_B, bins=[bins_x, bins_y], normed=True)
H, _, _ = np.histogram2d(log_rho, log_B, bins=[bins_x, bins_y])
# Plot the interesting quantities
plt.figure()
-#plt.pcolormesh(bins_x, bins_y, np.log10(H).T)
+# plt.pcolormesh(bins_x, bins_y, np.log10(H).T)
plt.pcolormesh(bins_x, bins_y, H.T)
plt.text(-5, 8.0, "$z=%.2f$" % z)
@@ -99,15 +99,15 @@ plt.xticks(
[-5, -4, -3, -2, -1, 0, 1, 2, 3],
["", "$10^{-4}$", "", "$10^{-2}$", "", "$10^0$", "", "$10^2$", ""],
)
-#plt.yticks(
+# plt.yticks(
# [2, 3, 4, 5, 6, 7, 8], ["$10^{2}$", "", "$10^{4}$", "", "$10^{6}$", "", "$10^8$"]
-#)
-plt.plot([-2, 3], [-13, 2/3*(3+2)-13], lw=1.5, alpha=0.7,color="red")
+# )
+plt.plot([-2, 3], [-13, 2 / 3 * (3 + 2) - 13], lw=1.5, alpha=0.7, color="red")
plt.xlabel("${\\rm Physical~Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
plt.ylabel("${\\rm Magnetic~Field}~B[{\\rm G}]$", labelpad=0)
plt.xlim(log_rho_min, log_rho_max)
plt.ylim(log_B_min, log_B_max)
-#plt.tight_layout()
+# plt.tight_layout()
plt.savefig("rhoB_%04d.png" % snap, dpi=200)
diff --git a/examples/MHDTests/CosmoVolume/plot_all.py b/examples/MHDTests/CosmoVolume/plot_all.py
index 6d8a49e6f3..dcd17322f3 100644
--- a/examples/MHDTests/CosmoVolume/plot_all.py
+++ b/examples/MHDTests/CosmoVolume/plot_all.py
@@ -9,6 +9,7 @@ from matplotlib.colors import LogNorm
from matplotlib.colors import Normalize
from mpl_toolkits.axes_grid1 import make_axes_locatable
+
def set_colorbar(ax, im):
"""
Adapt the colorbar a bit for axis object <ax> and
@@ -19,29 +20,30 @@ def set_colorbar(ax, im):
plt.colorbar(im, cax=cax)
return
+
filename_base = "snapshots/snap_"
-nini=int(sys.argv[1])
-nfin=int(sys.argv[2])
-#for ii in range(61):
-for ii in range(nini,nfin):
+nini = int(sys.argv[1])
+nfin = int(sys.argv[2])
+# for ii in range(61):
+for ii in range(nini, nfin):
print(ii)
filename = filename_base + str(ii).zfill(4) + ".hdf5"
- #mask=sw_mask(filename)
- #bs=mask.metadata.boxsize
- #print(bs)
- #mask.constrain_spatial([[0.,1.0],[0.,1.0],[0.,0.5]])
- #mask.constrain_spatial([[0.*bs[0],bs[0]],[0.*bs[0], 0.5*bs[1]],[0.*bs[2],bs[2]]])
- #data = load(filename,mask=mask)
+ # mask=sw_mask(filename)
+ # bs=mask.metadata.boxsize
+ # print(bs)
+ # mask.constrain_spatial([[0.,1.0],[0.,1.0],[0.,0.5]])
+ # mask.constrain_spatial([[0.*bs[0],bs[0]],[0.*bs[0], 0.5*bs[1]],[0.*bs[2],bs[2]]])
+ # data = load(filename,mask=mask)
data = load(filename)
- #print(data.metadata.gas_properties.field_names)
+ # print(data.metadata.gas_properties.field_names)
boxsize = data.metadata.boxsize
extent = [0.0, boxsize[0], 0.0, boxsize[1]]
-
+
gas_gamma = data.metadata.gas_gamma
- print("Gas Gamma:",gas_gamma)
+ print("Gas Gamma:", gas_gamma)
mhdflavour = data.metadata.hydro_scheme["MHD Flavour"]
mhd_scheme = data.metadata.hydro_scheme["MHD Scheme"]
@@ -51,7 +53,7 @@ for ii in range(nini,nfin):
scheme = data.metadata.hydro_scheme["Scheme"]
kernel = data.metadata.hydro_scheme["Kernel function"]
neighbours = data.metadata.hydro_scheme["Kernel target N_ngb"]
-
+
try:
dedhyp = data.metadata.hydro_scheme["Dedner Hyperbolic Constant"]
dedpar = data.metadata.hydro_scheme["Dedner Parabolic Constant"]
@@ -61,35 +63,37 @@ for ii in range(nini,nfin):
try:
deddivV = data.metadata.hydro_scheme["Dedner Hyperbolic div(v) Constant"]
- tensile = data.metadata.hydro_scheme["MHD Tensile Instability Correction Prefactor"]
+ tensile = data.metadata.hydro_scheme[
+ "MHD Tensile Instability Correction Prefactor"
+ ]
artdiff = data.metadata.hydro_scheme["Artificial Diffusion Constant"]
except:
deddivV = 0.0
artdiff = 0.0
tensile = 1.0
-
+
# First create a mass-weighted temperature dataset
-
+
B = data.gas.magnetic_flux_densities
divB = data.gas.magnetic_divergences
- P_mag = (B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2)
+ P_mag = B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2
P_mag = np.sqrt(P_mag)
h = data.gas.smoothing_lengths
- #A = data.gas.magnetic_vector_potentials
-
+ # A = data.gas.magnetic_vector_potentials
+
normB = np.sqrt(B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2)
- print("Bfield Min/Max: [",min(normB),"/",max(normB),"]")
-
+ print("Bfield Min/Max: [", min(normB), "/", max(normB), "]")
+
DivB_error = np.maximum(h * abs(divB) / normB, 1e-10)
- mmasses=data.gas.masses
+ mmasses = data.gas.masses
+
+ pressure = data.gas.pressures
- pressure=data.gas.pressures
-
# Then create a mass-weighted B error dataset
- #cdata.gas.mass_weighted_magnetic_divB_error = mmasses * DivB_error
- data.gas.mass_weighted_magnetic_divB_error = DivB_error #/ mmasses
-
+ # cdata.gas.mass_weighted_magnetic_divB_error = mmasses * DivB_error
+ data.gas.mass_weighted_magnetic_divB_error = DivB_error # / mmasses
+
# Then create a mass-weighted B pressure dataset
data.gas.mass_weighted_magnetic_pressures = mmasses * P_mag
@@ -104,95 +108,116 @@ for ii in range(nini,nfin):
data.gas.mass_weighted_speeds = data.gas.masses * (
v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2
)
-
+
# Then create a mass-weighted densities dataset
data.gas.mass_weighted_densities = data.gas.masses * data.gas.densities
# Map in mass per area
mass_map = project_gas(data, resolution=1024, project="masses", parallel=True)
-
+
# Map in density per area
mw_density_map = project_gas(
data, resolution=1024, project="mass_weighted_densities", parallel=True
)
# Map in magnetism squared times mass per area
mw_magnetic_pressure_map = project_gas(
- data, resolution=1024, project="mass_weighted_magnetic_pressures", parallel=True)
+ data, resolution=1024, project="mass_weighted_magnetic_pressures", parallel=True
+ )
# Map in pressure times mass per area
mw_pressure_map = project_gas(
- data, resolution=1024, project="mass_weighted_pressures", parallel=True, )
+ data, resolution=1024, project="mass_weighted_pressures", parallel=True
+ )
# Map in speed squared times mass per area
mw_speeds_map = project_gas(
- data, resolution=1024, project="mass_weighted_speeds", parallel=True, )
+ data, resolution=1024, project="mass_weighted_speeds", parallel=True
+ )
# Map in divB error times mass per area
mw_ErrDivB_map = project_gas(
- data, resolution=1024, project="mass_weighted_magnetic_divB_error", parallel=True, )
-
+ data,
+ resolution=1024,
+ project="mass_weighted_magnetic_divB_error",
+ parallel=True,
+ )
+
# Map in Plasma Beta times mass per area
plasma_beta_map = project_gas(
- data, resolution=1024, project="plasma_beta", parallel=True, )
+ data, resolution=1024, project="plasma_beta", parallel=True
+ )
rho_map = mw_density_map / mass_map
magnetic_pressure_map = mw_magnetic_pressure_map / mass_map
speed_map = mw_speeds_map / mass_map
pressure_map = mw_pressure_map / mass_map
- # ErrDivB_map = np.maximum( float(mw_ErrDivB_map[:]), 1E-10) / mass_map
+ # ErrDivB_map = np.maximum( float(mw_ErrDivB_map[:]), 1E-10) / mass_map
ErrDivB_map = mw_ErrDivB_map * mass_map / rho_map
- #plasma_beta_map
+ # plasma_beta_map
- #fig = plt.figure(figsize=(12, 11), dpi=100)
+ # fig = plt.figure(figsize=(12, 11), dpi=100)
fig = plt.figure(figsize=(12, 8), dpi=100)
-
+
ax1 = fig.add_subplot(231)
- im1 = ax1.imshow(rho_map.T, origin="lower", extent=extent, cmap="inferno", norm=LogNorm())
+ im1 = ax1.imshow(
+ rho_map.T, origin="lower", extent=extent, cmap="inferno", norm=LogNorm()
+ )
ax1.set_title("Density")
set_colorbar(ax1, im1)
-
+
ax2 = fig.add_subplot(232)
- #im2 = ax2.imshow(magnetic_pressure_map.T, origin="lower", extent=extent, cmap="plasma", norm=Normalize(vmax=120,vmin=30))
- im2 = ax2.imshow(magnetic_pressure_map.T, origin="lower", extent=extent, cmap="plasma", norm=LogNorm())
+ # im2 = ax2.imshow(magnetic_pressure_map.T, origin="lower", extent=extent, cmap="plasma", norm=Normalize(vmax=120,vmin=30))
+ im2 = ax2.imshow(
+ magnetic_pressure_map.T,
+ origin="lower",
+ extent=extent,
+ cmap="plasma",
+ norm=LogNorm(),
+ )
ax2.set_title("Magnetic field [G]")
set_colorbar(ax2, im2)
-
+
ax3 = fig.add_subplot(233)
- im3 = ax3.imshow(speed_map.T, origin="lower", extent=extent, cmap="cividis", norm=LogNorm())
+ im3 = ax3.imshow(
+ speed_map.T, origin="lower", extent=extent, cmap="cividis", norm=LogNorm()
+ )
ax3.set_title("Speed")
set_colorbar(ax3, im3)
-
+
ax4 = fig.add_subplot(234)
- im4 = ax4.imshow(pressure_map.T, origin="lower", extent=extent, cmap="viridis", norm=LogNorm())
+ im4 = ax4.imshow(
+ pressure_map.T, origin="lower", extent=extent, cmap="viridis", norm=LogNorm()
+ )
ax4.set_title("Internal Pressure")
set_colorbar(ax4, im4)
-
- #ax5 = fig.add_subplot(235)
- #im5 = ax5.imshow(plasma_beta_map.T, origin="lower", extent=extent, cmap="magma", norm=LogNorm())
- #ax5.set_title("plasma beta")
- #set_colorbar(ax5, im5)
-
+
+ # ax5 = fig.add_subplot(235)
+ # im5 = ax5.imshow(plasma_beta_map.T, origin="lower", extent=extent, cmap="magma", norm=LogNorm())
+ # ax5.set_title("plasma beta")
+ # set_colorbar(ax5, im5)
+
ax5 = fig.add_subplot(235)
- #im5 = ax5.imshow(ErrDivB_map.T, origin="lower", extent=extent, cmap="gray", norm=LogNorm(vmax=1,vmin=0.01))
- im5 = ax5.imshow(ErrDivB_map.T, origin="lower", extent=extent, cmap="gray", norm=LogNorm())
+ # im5 = ax5.imshow(ErrDivB_map.T, origin="lower", extent=extent, cmap="gray", norm=LogNorm(vmax=1,vmin=0.01))
+ im5 = ax5.imshow(
+ ErrDivB_map.T, origin="lower", extent=extent, cmap="gray", norm=LogNorm()
+ )
ax5.set_title("Err(DivB)")
set_colorbar(ax5, im5)
for ax in [ax1, ax2, ax3, ax4, ax5]:
ax.set_xlabel("x ")
ax.set_ylabel("y ")
-
+
ax6 = fig.add_subplot(236)
-
+
text_fontsize = 8
ax6.text(
- 0.1,
- 0.9,
- "Cosmo run $t=%.2f$" % data.metadata.time,
- fontsize=text_fontsize,
+ 0.1, 0.9, "Cosmo run $t=%.2f$" % data.metadata.time, fontsize=text_fontsize
)
ax6.text(0.1, 0.85, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
- ax6.text(0.1, 0.8, "$Branch$ %s" % gitBranch.decode("utf-8"), fontsize=text_fontsize)
+ ax6.text(
+ 0.1, 0.8, "$Branch$ %s" % gitBranch.decode("utf-8"), fontsize=text_fontsize
+ )
ax6.text(0.1, 0.75, scheme.decode("utf-8"), fontsize=text_fontsize)
ax6.text(0.1, 0.7, kernel.decode("utf-8"), fontsize=text_fontsize)
ax6.text(0.1, 0.65, "$%.2f$ neighbours" % (neighbours), fontsize=text_fontsize)
@@ -204,21 +229,28 @@ for ii in range(nini,nfin):
)
ax6.text(0.1, 0.5, "$Resitivity_\\eta:%.4f$ " % (mhdeta), fontsize=text_fontsize)
ax6.text(0.1, 0.45, "$Dedner Parameters: $", fontsize=text_fontsize)
- ax6.text(0.1, 0.4, "$[hyp, par, div] [:%.3f,%.3f,%.3f]$ " % (dedhyp,dedpar,deddivV), fontsize=text_fontsize)
+ ax6.text(
+ 0.1,
+ 0.4,
+ "$[hyp, par, div] [:%.3f,%.3f,%.3f]$ " % (dedhyp, dedpar, deddivV),
+ fontsize=text_fontsize,
+ )
ax6.text(0.1, 0.35, "$Tensile Prefactor:%.4f$ " % (tensile), fontsize=text_fontsize)
ax6.text(0.1, 0.3, "$Art. Diffusion:%.4f$ " % (artdiff), fontsize=text_fontsize)
- ax6.tick_params(left=False, right=False, labelleft=False, labelbottom=False, bottom=False)
+ ax6.tick_params(
+ left=False, right=False, labelleft=False, labelbottom=False, bottom=False
+ )
plt.tight_layout()
plt.savefig(filename_base + str(ii).zfill(4) + ".jpg")
plt.close()
- #from matplotlib.pyplot import imsave
+ # from matplotlib.pyplot import imsave
# Normalize and save
- #imsave(
+ # imsave(
# input_filename_base + str(ii).zfill(4) + ".png",
# np.rot90(density_map.value),
# cmap="jet",
# vmin=0.1,
# vmax=2.4,
- #)
+ # )
diff --git a/examples/MHDTests/FastRotor/makeIC.py b/examples/MHDTests/FastRotor/makeIC.py
index 22f4daafaa..58fefc4d30 100644
--- a/examples/MHDTests/FastRotor/makeIC.py
+++ b/examples/MHDTests/FastRotor/makeIC.py
@@ -134,8 +134,8 @@ rot = np.sqrt(
) # + (pos[:,2]-lz_c)**2)
theta = np.arctan2((pos[:, 1] - ly_c), (pos[:, 0] - lx_c))
v = np.zeros((N, 3))
-#B = np.zeros((N, 3))
-#A = np.zeros((N, 3))
+# B = np.zeros((N, 3))
+# A = np.zeros((N, 3))
ids = np.linspace(1, N, N)
m = np.ones(N) * rho_out_0 * vol / N_out
u = np.ones(N)
@@ -145,49 +145,49 @@ 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:])
-#B[:, 0] = B_0
+# B[:, 0] = B_0
###---------------------------###
-N2=int(2*N)
-p=np.zeros((N2, 3))
-p[:N,0]=pos[:,0]
-p[N:,0]=pos[:,0]
-p[:N,1]=pos[:,1]
-p[N:,1]=pos[:,1]+1.0
-p[:N,2]=pos[:,2]
-p[N:,2]=pos[:,2]
-pos=p
-hh =np.zeros(N2)
-hh[:N]=h
-hh[N:]=h
-h=hh
-vel=np.zeros((N2,3))
-vel[:N,:]=v[:,:]
-vel[N: ,:]=v[:,:]
-v=vel
+N2 = int(2 * N)
+p = np.zeros((N2, 3))
+p[:N, 0] = pos[:, 0]
+p[N:, 0] = pos[:, 0]
+p[:N, 1] = pos[:, 1]
+p[N:, 1] = pos[:, 1] + 1.0
+p[:N, 2] = pos[:, 2]
+p[N:, 2] = pos[:, 2]
+pos = p
+hh = np.zeros(N2)
+hh[:N] = h
+hh[N:] = h
+h = hh
+vel = np.zeros((N2, 3))
+vel[:N, :] = v[:, :]
+vel[N:, :] = v[:, :]
+v = vel
ids = np.linspace(1, N2, N2)
m = np.ones(N2) * rho_out_0 * vol / N_out
-uu =np.zeros(N2)
-uu[:N]=u
-uu[N:]=u
-u=uu
+uu = np.zeros(N2)
+uu[:N] = u
+uu[N:] = u
+u = uu
B = np.zeros((N2, 3))
A = np.zeros((N2, 3))
B[:N, 0] = B_0
B[N:, 0] = -B_0
-A[:N, 2] = B_0 * pos[:N,1]
-A[N:, 2] = B_0 * (2.0-pos[N:,1])
+A[:N, 2] = B_0 * pos[:N, 1]
+A[N:, 2] = B_0 * (2.0 - pos[N:, 1])
# File
fileOutput = h5py.File(fileOutputName, "w")
# Header
grp = fileOutput.create_group("/Header")
-grp.attrs["BoxSize"] = [lx, 2.* ly, lz] #####
-grp.attrs["NumPart_Total"] = [2*N, 0, 0, 0, 0, 0]
+grp.attrs["BoxSize"] = [lx, 2.0 * ly, lz] #####
+grp.attrs["NumPart_Total"] = [2 * N, 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_ThisFile"] = [2*N, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [2 * N, 0, 0, 0, 0, 0]
grp.attrs["Time"] = 0.0
grp.attrs["NumFileOutputsPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
@@ -211,7 +211,7 @@ grp.create_dataset("SmoothingLength", data=h, dtype="f")
grp.create_dataset("InternalEnergy", data=u, dtype="f")
grp.create_dataset("ParticleIDs", data=ids, dtype="L")
grp.create_dataset("MagneticFluxDensities", data=B, dtype="f")
-grp.create_dataset("MagneticVectorPotentials", data = A, dtype = 'f')
+grp.create_dataset("MagneticVectorPotentials", data=A, dtype="f")
# grp.create_dataset("EPalpha", data = epa, dtype = 'f')
# grp.create_dataset("EPbeta" , data = epb, dtype = 'f')
diff --git a/examples/MHDTests/FastRotor/plot_all.py b/examples/MHDTests/FastRotor/plot_all.py
index 62466c2896..029b29197d 100644
--- a/examples/MHDTests/FastRotor/plot_all.py
+++ b/examples/MHDTests/FastRotor/plot_all.py
@@ -4,13 +4,15 @@ import numpy as np
import sys
import os
import matplotlib
-#matplotlib.use("pdf")
+
+# matplotlib.use("pdf")
import matplotlib.pyplot as plt
from matplotlib.pyplot import imsave
from matplotlib.colors import LogNorm
from matplotlib.colors import Normalize
from mpl_toolkits.axes_grid1 import make_axes_locatable
+
def set_colorbar(ax, im):
"""
Adapt the colorbar a bit for axis object <ax> and
@@ -21,26 +23,27 @@ def set_colorbar(ax, im):
plt.colorbar(im, cax=cax)
return
+
filename_base = "FastRotor_LR_"
-nini=int(sys.argv[1])
-nfin=int(sys.argv[2])
-#for ii in range(61):
-for ii in range(nini,nfin):
+nini = int(sys.argv[1])
+nfin = int(sys.argv[2])
+# for ii in range(61):
+for ii in range(nini, nfin):
print(ii)
filename = filename_base + str(ii).zfill(4) + ".hdf5"
data = load(filename)
- #print(data.metadata.gas_properties.field_names)
+ # print(data.metadata.gas_properties.field_names)
boxsize = data.metadata.boxsize
- extent = [0, 1.0 , 0, 1.0]
- # cut the domian in half
- data.metadata.boxsize*=[1.0,0.5,1.0]
-
+ extent = [0, 1.0, 0, 1.0]
+ # cut the domian in half
+ data.metadata.boxsize *= [1.0, 0.5, 1.0]
+
gas_gamma = data.metadata.gas_gamma
- print("Gas Gamma:",gas_gamma)
- if gas_gamma != 7./5.:
+ print("Gas Gamma:", gas_gamma)
+ if gas_gamma != 7.0 / 5.0:
print("WRONG GAS GAMMA")
exit()
@@ -52,7 +55,7 @@ for ii in range(nini,nfin):
scheme = data.metadata.hydro_scheme["Scheme"]
kernel = data.metadata.hydro_scheme["Kernel function"]
neighbours = data.metadata.hydro_scheme["Kernel target N_ngb"]
-
+
try:
dedhyp = data.metadata.hydro_scheme["Dedner Hyperbolic Constant"]
dedpar = data.metadata.hydro_scheme["Dedner Parabolic Constant"]
@@ -62,27 +65,29 @@ for ii in range(nini,nfin):
try:
deddivV = data.metadata.hydro_scheme["Dedner Hyperbolic div(v) Constant"]
- tensile = data.metadata.hydro_scheme["MHD Tensile Instability Correction Prefactor"]
+ tensile = data.metadata.hydro_scheme[
+ "MHD Tensile Instability Correction Prefactor"
+ ]
artdiff = data.metadata.hydro_scheme["Artificial Diffusion Constant"]
except:
deddivV = 0.0
artdiff = 0.0
tensile = 1.0
-
+
# First create a mass-weighted temperature dataset
-
+
B = data.gas.magnetic_flux_densities
divB = data.gas.magnetic_divergences
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)
-
+
# Then create a mass-weighted B error dataset
data.gas.mass_weighted_magnetic_divB_error = data.gas.masses * data.gas.DivB_error
-
+
# Then create a mass-weighted B pressure dataset
data.gas.mass_weighted_magnetic_pressures = data.gas.masses * P_mag
@@ -101,92 +106,151 @@ for ii in range(nini,nfin):
data.gas.mass_weighted_densities = data.gas.masses * data.gas.densities
# Map in mass per area
- mass_map = slice_gas(data, z_slice=0.5 * data.metadata.boxsize[2], resolution=1024, project="masses", parallel=True )
-
+ mass_map = slice_gas(
+ data,
+ z_slice=0.5 * data.metadata.boxsize[2],
+ resolution=1024,
+ project="masses",
+ parallel=True,
+ )
+
# Map in density per area
rho_map = slice_gas(
- data, z_slice=0.5 * data.metadata.boxsize[2], resolution=1024,
- project="mass_weighted_densities", parallel=True)
+ data,
+ z_slice=0.5 * data.metadata.boxsize[2],
+ resolution=1024,
+ project="mass_weighted_densities",
+ parallel=True,
+ )
# Map in magnetism squared times mass per area
mw_magnetic_pressure_map = slice_gas(
- data, z_slice=0.5 * data.metadata.boxsize[2], resolution=1024,
- project="mass_weighted_magnetic_pressures", parallel=True)
+ data,
+ z_slice=0.5 * data.metadata.boxsize[2],
+ resolution=1024,
+ project="mass_weighted_magnetic_pressures",
+ parallel=True,
+ )
# Map in pressure times mass per area
mw_pressure_map = slice_gas(
- data, z_slice=0.5 * data.metadata.boxsize[2], resolution=1024,
- project="mass_weighted_pressures", parallel=True, )
+ data,
+ z_slice=0.5 * data.metadata.boxsize[2],
+ resolution=1024,
+ project="mass_weighted_pressures",
+ parallel=True,
+ )
# Map in speed squared times mass per area
mw_speeds_map = slice_gas(
- data, z_slice=0.5 * data.metadata.boxsize[2], resolution=1024,
- project="mass_weighted_speeds", parallel=True, )
-
+ data,
+ z_slice=0.5 * data.metadata.boxsize[2],
+ resolution=1024,
+ project="mass_weighted_speeds",
+ parallel=True,
+ )
+
# Map in divB error times mass per area
mw_ErrDivB_map = slice_gas(
- data, z_slice=0.5 * data.metadata.boxsize[2], resolution=1024,
- project="mass_weighted_magnetic_divB_error", parallel=True, )
-
+ data,
+ z_slice=0.5 * data.metadata.boxsize[2],
+ resolution=1024,
+ project="mass_weighted_magnetic_divB_error",
+ parallel=True,
+ )
+
# Map in Plasma Beta times mass per area
plasma_beta_map = slice_gas(
- data, z_slice=0.5 * data.metadata.boxsize[2], resolution=1024,
- project="plasma_beta", parallel=True, )
+ data,
+ z_slice=0.5 * data.metadata.boxsize[2],
+ resolution=1024,
+ project="plasma_beta",
+ parallel=True,
+ )
rho_map = rho_map / mass_map
magnetic_pressure_map = mw_magnetic_pressure_map / mass_map
speed_map = mw_speeds_map / mass_map
pressure_map = mw_pressure_map / mass_map
ErrDivB_map = mw_ErrDivB_map / mass_map
- #plasma_beta_map
+ # plasma_beta_map
- #fig = plt.figure(figsize=(12, 11), dpi=100)
+ # fig = plt.figure(figsize=(12, 11), dpi=100)
fig = plt.figure(figsize=(12, 8), dpi=100)
-
+
ax1 = fig.add_subplot(231)
- im1 = ax1.imshow(rho_map.T, origin="lower", extent=extent, cmap="inferno", norm=LogNorm(vmax=14,vmin=1))
+ im1 = ax1.imshow(
+ rho_map.T,
+ origin="lower",
+ extent=extent,
+ cmap="inferno",
+ norm=LogNorm(vmax=14, vmin=1),
+ )
ax1.set_title("Density")
set_colorbar(ax1, im1)
-
+
ax2 = fig.add_subplot(232)
- im2 = ax2.imshow(magnetic_pressure_map.T, origin="lower", extent=extent, cmap="plasma", norm=Normalize(vmax=3,vmin=0.1))
+ im2 = ax2.imshow(
+ magnetic_pressure_map.T,
+ origin="lower",
+ extent=extent,
+ cmap="plasma",
+ norm=Normalize(vmax=3, vmin=0.1),
+ )
ax2.set_title("Magnetic Pressure")
set_colorbar(ax2, im2)
-
+
ax3 = fig.add_subplot(233)
- im3 = ax3.imshow(speed_map.T, origin="lower", extent=extent, cmap="cividis", norm=Normalize(vmax=1.6,vmin=0.1))
+ im3 = ax3.imshow(
+ speed_map.T,
+ origin="lower",
+ extent=extent,
+ cmap="cividis",
+ norm=Normalize(vmax=1.6, vmin=0.1),
+ )
ax3.set_title("v^2")
set_colorbar(ax3, im3)
-
+
ax4 = fig.add_subplot(234)
- im4 = ax4.imshow(pressure_map.T, origin="lower", extent=extent, cmap="viridis", norm=Normalize(vmax=2.1,vmin=0.1))
+ im4 = ax4.imshow(
+ pressure_map.T,
+ origin="lower",
+ extent=extent,
+ cmap="viridis",
+ norm=Normalize(vmax=2.1, vmin=0.1),
+ )
ax4.set_title("Internal Pressure")
set_colorbar(ax4, im4)
-
- #ax5 = fig.add_subplot(235)
- #im5 = ax5.imshow(plasma_beta_map.T, origin="lower", extent=extent, cmap="magma", norm=LogNorm())
- #ax5.set_title("plasma beta")
- #set_colorbar(ax5, im5)
-
+
+ # ax5 = fig.add_subplot(235)
+ # im5 = ax5.imshow(plasma_beta_map.T, origin="lower", extent=extent, cmap="magma", norm=LogNorm())
+ # ax5.set_title("plasma beta")
+ # set_colorbar(ax5, im5)
+
ax5 = fig.add_subplot(235)
- im5 = ax5.imshow(ErrDivB_map.T, origin="lower", extent=extent, cmap="gray", norm=LogNorm(vmax=1,vmin=0.001))
+ im5 = ax5.imshow(
+ ErrDivB_map.T,
+ origin="lower",
+ extent=extent,
+ cmap="gray",
+ norm=LogNorm(vmax=1, vmin=0.001),
+ )
ax5.set_title("Err(DivB)")
set_colorbar(ax5, im5)
for ax in [ax1, ax2, ax3, ax4, ax5]:
ax.set_xlabel("x ")
ax.set_ylabel("y ")
-
+
ax6 = fig.add_subplot(236)
-
+
text_fontsize = 8
ax6.text(
- 0.1,
- 0.9,
- "Fast Rotor $t=%.2f$" % data.metadata.time,
- fontsize=text_fontsize,
+ 0.1, 0.9, "Fast Rotor $t=%.2f$" % data.metadata.time, fontsize=text_fontsize
)
ax6.text(0.1, 0.85, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
- ax6.text(0.1, 0.8, "$Branch$ %s" % gitBranch.decode("utf-8"), fontsize=text_fontsize)
+ ax6.text(
+ 0.1, 0.8, "$Branch$ %s" % gitBranch.decode("utf-8"), fontsize=text_fontsize
+ )
ax6.text(0.1, 0.75, scheme.decode("utf-8"), fontsize=text_fontsize)
ax6.text(0.1, 0.7, kernel.decode("utf-8"), fontsize=text_fontsize)
ax6.text(0.1, 0.65, "$%.2f$ neighbours" % (neighbours), fontsize=text_fontsize)
@@ -198,12 +262,19 @@ for ii in range(nini,nfin):
)
ax6.text(0.1, 0.5, "$Resitivity_\\eta:%.4f$ " % (mhdeta), fontsize=text_fontsize)
ax6.text(0.1, 0.45, "$Dedner Parameters: $", fontsize=text_fontsize)
- ax6.text(0.1, 0.4, "$[hyp, par, div] [:%.3f,%.3f,%.3f]$ " % (dedhyp,dedpar,deddivV), fontsize=text_fontsize)
+ ax6.text(
+ 0.1,
+ 0.4,
+ "$[hyp, par, div] [:%.3f,%.3f,%.3f]$ " % (dedhyp, dedpar, deddivV),
+ fontsize=text_fontsize,
+ )
ax6.text(0.1, 0.35, "$Tensile Prefactor:%.4f$ " % (tensile), fontsize=text_fontsize)
ax6.text(0.1, 0.3, "$Art. Diffusion:%.4f$ " % (artdiff), fontsize=text_fontsize)
- ax6.tick_params(left=False, right=False, labelleft=False, labelbottom=False, bottom=False)
- #ax6.set_xlabel("")
- #ax6.plot(frameon=False)
+ ax6.tick_params(
+ left=False, right=False, labelleft=False, labelbottom=False, bottom=False
+ )
+ # ax6.set_xlabel("")
+ # ax6.plot(frameon=False)
plt.tight_layout()
plt.savefig(filename_base + str(ii).zfill(4) + ".jpg")
diff --git a/examples/MHDTests/MagneticBlastWave/makeIC.py b/examples/MHDTests/MagneticBlastWave/makeIC.py
index c9fa05b9c1..99874e572d 100644
--- a/examples/MHDTests/MagneticBlastWave/makeIC.py
+++ b/examples/MHDTests/MagneticBlastWave/makeIC.py
@@ -13,7 +13,7 @@ r_in = 0.125
rho_0 = 1.0
P_in_0 = 100.0
P_out_0 = 1
-B_0 = 10.0
+B_0 = 10.0
gamma = 5.0 / 3.0
fileOutputName = "MagneticBlastWave_LR.hdf5"
@@ -70,47 +70,47 @@ vol = lx * ly * lz
###---------------------------###
rot = np.sqrt((pos[:, 0] - lx_c) ** 2 + (pos[:, 1] - ly_c) ** 2)
-#v = np.zeros((N, 3))
-#B = np.zeros((N, 3))
-#ids = np.linspace(1, N, N)
-#m = np.ones(N) * rho_0 * vol / N
+# v = np.zeros((N, 3))
+# B = np.zeros((N, 3))
+# ids = np.linspace(1, N, N)
+# m = np.ones(N) * rho_0 * vol / N
u = np.ones(N)
u[rot < r_in] = P_in_0 / (rho_0 * (gamma - 1))
u[rot >= r_in] = P_out_0 / (rho_0 * (gamma - 1))
-#B[:, 0] = B_0
+# B[:, 0] = B_0
-print("Max Pos: [",max(pos[:,0]),max(pos[:,1]),max(pos[:,2]),"]")
-print("Min Pos: [",min(pos[:,0]),min(pos[:,1]),min(pos[:,2]),"]")
+print("Max Pos: [", max(pos[:, 0]), max(pos[:, 1]), max(pos[:, 2]), "]")
+print("Min Pos: [", min(pos[:, 0]), min(pos[:, 1]), min(pos[:, 2]), "]")
###---------------------------###
-N2=int(2*N)
-p=np.zeros((N2, 3))
-p[:N,0]=pos[:,0]
-p[N:,0]=pos[:,0]
-p[:N,1]=pos[:,1]
-p[N:,1]=pos[:,1]+ly
-p[:N,2]=pos[:,2]
-p[N:,2]=pos[:,2]
-pos=p
-hh =np.zeros(N2)
-hh[:N]=h
-hh[N:]=h
-h=hh
-v=np.zeros((N2,3))
+N2 = int(2 * N)
+p = np.zeros((N2, 3))
+p[:N, 0] = pos[:, 0]
+p[N:, 0] = pos[:, 0]
+p[:N, 1] = pos[:, 1]
+p[N:, 1] = pos[:, 1] + ly
+p[:N, 2] = pos[:, 2]
+p[N:, 2] = pos[:, 2]
+pos = p
+hh = np.zeros(N2)
+hh[:N] = h
+hh[N:] = h
+h = hh
+v = np.zeros((N2, 3))
ids = np.linspace(1, N2, N2)
m = np.ones(N2) * rho_0 * vol / N
-uu =np.zeros(N2)
-uu[:N]=u
-uu[N:]=u
-u=uu
+uu = np.zeros(N2)
+uu[:N] = u
+uu[N:] = u
+u = uu
B = np.zeros((N2, 3))
A = np.zeros((N2, 3))
B[:N, 0] = B_0
B[N:, 0] = -B_0
-A[:N, 2] = B_0 * pos[:N,1]
-A[N:, 2] = B_0 * (2*ly-pos[N:,1])
+A[:N, 2] = B_0 * pos[:N, 1]
+A[N:, 2] = B_0 * (2 * ly - pos[N:, 1])
-print("Max Pos: [",max(pos[:,0]),max(pos[:,1]),max(pos[:,2]),"]")
-print("Min Pos: [",min(pos[:,0]),min(pos[:,1]),min(pos[:,2]),"]")
+print("Max Pos: [", max(pos[:, 0]), max(pos[:, 1]), max(pos[:, 2]), "]")
+print("Min Pos: [", min(pos[:, 0]), min(pos[:, 1]), min(pos[:, 2]), "]")
# File
@@ -118,7 +118,7 @@ fileOutput = h5py.File(fileOutputName, "w")
# Header
grp = fileOutput.create_group("/Header")
-grp.attrs["BoxSize"] = [lx, 2. * ly, lz] #####
+grp.attrs["BoxSize"] = [lx, 2.0 * ly, lz] #####
grp.attrs["NumPart_Total"] = [N2, 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [N2, 0, 0, 0, 0, 0]
@@ -145,6 +145,6 @@ grp.create_dataset("SmoothingLength", data=h, dtype="f")
grp.create_dataset("InternalEnergy", data=u, dtype="f")
grp.create_dataset("ParticleIDs", data=ids, dtype="L")
grp.create_dataset("MagneticFluxDensities", data=B, dtype="f")
-grp.create_dataset("MagneticVectorPotentials", data = A, dtype = 'f')
+grp.create_dataset("MagneticVectorPotentials", data=A, dtype="f")
fileOutput.close()
diff --git a/examples/MHDTests/MagneticBlastWave/plot_all.py b/examples/MHDTests/MagneticBlastWave/plot_all.py
index 3705048ff1..ef47af117a 100644
--- a/examples/MHDTests/MagneticBlastWave/plot_all.py
+++ b/examples/MHDTests/MagneticBlastWave/plot_all.py
@@ -9,6 +9,7 @@ from matplotlib.colors import LogNorm
from matplotlib.colors import Normalize
from mpl_toolkits.axes_grid1 import make_axes_locatable
+
def set_colorbar(ax, im):
"""
Adapt the colorbar a bit for axis object <ax> and
@@ -19,31 +20,32 @@ def set_colorbar(ax, im):
plt.colorbar(im, cax=cax)
return
+
filename_base = "MagneticBlastWave_LR_"
-nini=int(sys.argv[1])
-nfin=int(sys.argv[2])
-#for ii in range(61):
-for ii in range(nini,nfin):
+nini = int(sys.argv[1])
+nfin = int(sys.argv[2])
+# for ii in range(61):
+for ii in range(nini, nfin):
print(ii)
filename = filename_base + str(ii).zfill(4) + ".hdf5"
- #mask=sw_mask(filename)
- #bs=mask.metadata.boxsize
- #print(bs)
- #mask.constrain_spatial([[0.,1.0],[0.,1.0],[0.,0.5]])
- #mask.constrain_spatial([[0.*bs[0],bs[0]],[0.*bs[0], 0.5*bs[1]],[0.*bs[2],bs[2]]])
- #data = load(filename,mask=mask)
+ # mask=sw_mask(filename)
+ # bs=mask.metadata.boxsize
+ # print(bs)
+ # mask.constrain_spatial([[0.,1.0],[0.,1.0],[0.,0.5]])
+ # mask.constrain_spatial([[0.*bs[0],bs[0]],[0.*bs[0], 0.5*bs[1]],[0.*bs[2],bs[2]]])
+ # data = load(filename,mask=mask)
data = load(filename)
- #print(data.metadata.gas_properties.field_names)
- #boxsize = data.metadata.boxsize
+ # print(data.metadata.gas_properties.field_names)
+ # boxsize = data.metadata.boxsize
extent = [0.0, 1.0, 0.0, 1.0]
- data.metadata.boxsize*=[1.0,0.5,1.0]
-
+ data.metadata.boxsize *= [1.0, 0.5, 1.0]
+
gas_gamma = data.metadata.gas_gamma
- print("Gas Gamma:",gas_gamma)
- if gas_gamma != 5./3.:
+ print("Gas Gamma:", gas_gamma)
+ if gas_gamma != 5.0 / 3.0:
print("WRONG GAS GAMMA")
exit()
@@ -55,7 +57,7 @@ for ii in range(nini,nfin):
scheme = data.metadata.hydro_scheme["Scheme"]
kernel = data.metadata.hydro_scheme["Kernel function"]
neighbours = data.metadata.hydro_scheme["Kernel target N_ngb"]
-
+
try:
dedhyp = data.metadata.hydro_scheme["Dedner Hyperbolic Constant"]
dedpar = data.metadata.hydro_scheme["Dedner Parabolic Constant"]
@@ -65,32 +67,34 @@ for ii in range(nini,nfin):
try:
deddivV = data.metadata.hydro_scheme["Dedner Hyperbolic div(v) Constant"]
- tensile = data.metadata.hydro_scheme["MHD Tensile Instability Correction Prefactor"]
+ tensile = data.metadata.hydro_scheme[
+ "MHD Tensile Instability Correction Prefactor"
+ ]
artdiff = data.metadata.hydro_scheme["Artificial Diffusion Constant"]
except:
deddivV = 0.0
artdiff = 0.0
tensile = 1.0
-
+
# First create a mass-weighted temperature dataset
-
+
B = data.gas.magnetic_flux_densities
divB = data.gas.magnetic_divergences
P_mag = (B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2) / 2
h = data.gas.smoothing_lengths
- #A = data.gas.magnetic_vector_potentials
-
+ # A = data.gas.magnetic_vector_potentials
+
normB = np.sqrt(B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2)
-
+
DivB_error = np.maximum(h * abs(divB) / normB, 1e-10)
- mmasses=data.gas.masses
+ mmasses = data.gas.masses
+
+ pressure = data.gas.pressures
- pressure=data.gas.pressures
-
# Then create a mass-weighted B error dataset
data.gas.mass_weighted_magnetic_divB_error = mmasses * DivB_error
-
+
# Then create a mass-weighted B pressure dataset
data.gas.mass_weighted_magnetic_pressures = mmasses * P_mag
@@ -105,93 +109,130 @@ for ii in range(nini,nfin):
data.gas.mass_weighted_speeds = data.gas.masses * (
v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2
)
-
+
# Then create a mass-weighted densities dataset
data.gas.mass_weighted_densities = data.gas.masses * data.gas.densities
# Map in mass per area
mass_map = project_gas(data, resolution=1024, project="masses", parallel=True)
-
+
# Map in density per area
mw_density_map = project_gas(
data, resolution=1024, project="mass_weighted_densities", parallel=True
)
# Map in magnetism squared times mass per area
mw_magnetic_pressure_map = project_gas(
- data, resolution=1024, project="mass_weighted_magnetic_pressures", parallel=True)
+ data, resolution=1024, project="mass_weighted_magnetic_pressures", parallel=True
+ )
# Map in pressure times mass per area
mw_pressure_map = project_gas(
- data, resolution=1024, project="mass_weighted_pressures", parallel=True, )
+ data, resolution=1024, project="mass_weighted_pressures", parallel=True
+ )
# Map in speed squared times mass per area
mw_speeds_map = project_gas(
- data, resolution=1024, project="mass_weighted_speeds", parallel=True, )
+ data, resolution=1024, project="mass_weighted_speeds", parallel=True
+ )
# Map in divB error times mass per area
mw_ErrDivB_map = project_gas(
- data, resolution=1024, project="mass_weighted_magnetic_divB_error", parallel=True, )
-
+ data,
+ resolution=1024,
+ project="mass_weighted_magnetic_divB_error",
+ parallel=True,
+ )
+
# Map in Plasma Beta times mass per area
plasma_beta_map = project_gas(
- data, resolution=1024, project="plasma_beta", parallel=True, )
+ data, resolution=1024, project="plasma_beta", parallel=True
+ )
rho_map = mw_density_map / mass_map
magnetic_pressure_map = mw_magnetic_pressure_map / mass_map
speed_map = mw_speeds_map / mass_map
pressure_map = mw_pressure_map / mass_map
- ErrDivB_map = np.maximum(mw_ErrDivB_map / mass_map,1E-10)
- #plasma_beta_map
+ ErrDivB_map = np.maximum(mw_ErrDivB_map / mass_map, 1e-10)
+ # plasma_beta_map
- #fig = plt.figure(figsize=(12, 11), dpi=100)
+ # fig = plt.figure(figsize=(12, 11), dpi=100)
fig = plt.figure(figsize=(12, 8), dpi=100)
-
+
ax1 = fig.add_subplot(231)
- im1 = ax1.imshow(rho_map.T, origin="lower", extent=extent, cmap="inferno", norm=Normalize(vmax=3,vmin=0))
+ im1 = ax1.imshow(
+ rho_map.T,
+ origin="lower",
+ extent=extent,
+ cmap="inferno",
+ norm=Normalize(vmax=3, vmin=0),
+ )
ax1.set_title("Density")
set_colorbar(ax1, im1)
-
+
ax2 = fig.add_subplot(232)
- #im2 = ax2.imshow(magnetic_pressure_map.T, origin="lower", extent=extent, cmap="plasma", norm=Normalize(vmax=120,vmin=30))
- im2 = ax2.imshow(magnetic_pressure_map.T, origin="lower", extent=extent, cmap="plasma", norm=Normalize(vmax=65,vmin=25))
+ # im2 = ax2.imshow(magnetic_pressure_map.T, origin="lower", extent=extent, cmap="plasma", norm=Normalize(vmax=120,vmin=30))
+ im2 = ax2.imshow(
+ magnetic_pressure_map.T,
+ origin="lower",
+ extent=extent,
+ cmap="plasma",
+ norm=Normalize(vmax=65, vmin=25),
+ )
ax2.set_title("Magnetic Pressure")
set_colorbar(ax2, im2)
-
+
ax3 = fig.add_subplot(233)
- im3 = ax3.imshow(speed_map.T, origin="lower", extent=extent, cmap="cividis", norm=Normalize(vmax=30,vmin=0.))
+ im3 = ax3.imshow(
+ speed_map.T,
+ origin="lower",
+ extent=extent,
+ cmap="cividis",
+ norm=Normalize(vmax=30, vmin=0.0),
+ )
ax3.set_title("Speed")
set_colorbar(ax3, im3)
-
+
ax4 = fig.add_subplot(234)
- im4 = ax4.imshow(pressure_map.T, origin="lower", extent=extent, cmap="viridis", norm=Normalize(vmax=50,vmin=0.))
+ im4 = ax4.imshow(
+ pressure_map.T,
+ origin="lower",
+ extent=extent,
+ cmap="viridis",
+ norm=Normalize(vmax=50, vmin=0.0),
+ )
ax4.set_title("Internal Pressure")
set_colorbar(ax4, im4)
-
- #ax5 = fig.add_subplot(235)
- #im5 = ax5.imshow(plasma_beta_map.T, origin="lower", extent=extent, cmap="magma", norm=LogNorm())
- #ax5.set_title("plasma beta")
- #set_colorbar(ax5, im5)
-
+
+ # ax5 = fig.add_subplot(235)
+ # im5 = ax5.imshow(plasma_beta_map.T, origin="lower", extent=extent, cmap="magma", norm=LogNorm())
+ # ax5.set_title("plasma beta")
+ # set_colorbar(ax5, im5)
+
ax5 = fig.add_subplot(235)
- im5 = ax5.imshow(ErrDivB_map.T, origin="lower", extent=extent, cmap="gray", norm=LogNorm(vmax=1,vmin=0.001))
+ im5 = ax5.imshow(
+ ErrDivB_map.T,
+ origin="lower",
+ extent=extent,
+ cmap="gray",
+ norm=LogNorm(vmax=1, vmin=0.001),
+ )
ax5.set_title("Err(DivB)")
set_colorbar(ax5, im5)
for ax in [ax1, ax2, ax3, ax4, ax5]:
ax.set_xlabel("x ")
ax.set_ylabel("y ")
-
+
ax6 = fig.add_subplot(236)
-
+
text_fontsize = 8
ax6.text(
- 0.1,
- 0.9,
- "Blast Wave $t=%.2f$" % data.metadata.time,
- fontsize=text_fontsize,
+ 0.1, 0.9, "Blast Wave $t=%.2f$" % data.metadata.time, fontsize=text_fontsize
)
ax6.text(0.1, 0.85, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
- ax6.text(0.1, 0.8, "$Branch$ %s" % gitBranch.decode("utf-8"), fontsize=text_fontsize)
+ ax6.text(
+ 0.1, 0.8, "$Branch$ %s" % gitBranch.decode("utf-8"), fontsize=text_fontsize
+ )
ax6.text(0.1, 0.75, scheme.decode("utf-8"), fontsize=text_fontsize)
ax6.text(0.1, 0.7, kernel.decode("utf-8"), fontsize=text_fontsize)
ax6.text(0.1, 0.65, "$%.2f$ neighbours" % (neighbours), fontsize=text_fontsize)
@@ -203,21 +244,28 @@ for ii in range(nini,nfin):
)
ax6.text(0.1, 0.5, "$Resitivity_\\eta:%.4f$ " % (mhdeta), fontsize=text_fontsize)
ax6.text(0.1, 0.45, "$Dedner Parameters: $", fontsize=text_fontsize)
- ax6.text(0.1, 0.4, "$[hyp, par, div] [:%.3f,%.3f,%.3f]$ " % (dedhyp,dedpar,deddivV), fontsize=text_fontsize)
+ ax6.text(
+ 0.1,
+ 0.4,
+ "$[hyp, par, div] [:%.3f,%.3f,%.3f]$ " % (dedhyp, dedpar, deddivV),
+ fontsize=text_fontsize,
+ )
ax6.text(0.1, 0.35, "$Tensile Prefactor:%.4f$ " % (tensile), fontsize=text_fontsize)
ax6.text(0.1, 0.3, "$Art. Diffusion:%.4f$ " % (artdiff), fontsize=text_fontsize)
- ax6.tick_params(left=False, right=False, labelleft=False, labelbottom=False, bottom=False)
+ ax6.tick_params(
+ left=False, right=False, labelleft=False, labelbottom=False, bottom=False
+ )
plt.tight_layout()
plt.savefig(filename_base + str(ii).zfill(4) + ".jpg")
plt.close()
- #from matplotlib.pyplot import imsave
+ # from matplotlib.pyplot import imsave
# Normalize and save
- #imsave(
+ # imsave(
# input_filename_base + str(ii).zfill(4) + ".png",
# np.rot90(density_map.value),
# cmap="jet",
# vmin=0.1,
# vmax=2.4,
- #)
+ # )
diff --git a/src/mhd/DInduction/mhd.h b/src/mhd/DInduction/mhd.h
index 62c51f14cf..511d2f7e5e 100644
--- a/src/mhd/DInduction/mhd.h
+++ b/src/mhd/DInduction/mhd.h
@@ -264,8 +264,8 @@ __attribute__((always_inline)) INLINE static float hydro_get_dphi_dt(
const struct part *restrict p, const float hyp, const float par,
const struct cosmology *c, const float mu0) {
- //const float v_sig = hydro_get_signal_velocity(p);
- const float v_sig = mhd_get_magnetosonic_speed(p,c->a,mu0);
+ // const float v_sig = hydro_get_signal_velocity(p);
+ const float v_sig = mhd_get_magnetosonic_speed(p, c->a, mu0);
// const float div_v = hydro_get_div_v(p);
const float afac1 = pow(c->a, 2.f * c->a_factor_sound_speed);
const float afac2 = pow(c->a, 1.f + c->a_factor_sound_speed);
@@ -468,7 +468,7 @@ __attribute__((always_inline)) INLINE static void mhd_predict_extra(
p->mhd_data.BPred[0] += p->mhd_data.dBdt[0] * dt_therm;
p->mhd_data.BPred[1] += p->mhd_data.dBdt[1] * dt_therm;
p->mhd_data.BPred[2] += p->mhd_data.dBdt[2] * dt_therm;
-
+
const float hyp = hydro_props->mhd.hyp_dedner;
const float par = hydro_props->mhd.par_dedner;
p->mhd_data.phi += hydro_get_dphi_dt(p, hyp, par, cosmo, mu_0) * dt_therm;
@@ -524,7 +524,6 @@ __attribute__((always_inline)) INLINE static void mhd_kick_extra(
xp->mhd_data.Bfld_full[0] += p->mhd_data.dBdt[0] * dt_therm;
xp->mhd_data.Bfld_full[1] += p->mhd_data.dBdt[1] * dt_therm;
xp->mhd_data.Bfld_full[2] += p->mhd_data.dBdt[2] * dt_therm;
-
}
/**
diff --git a/src/mhd/DirectInduction/mhd.h b/src/mhd/DirectInduction/mhd.h
index a69beaa534..774c029c67 100644
--- a/src/mhd/DirectInduction/mhd.h
+++ b/src/mhd/DirectInduction/mhd.h
@@ -232,12 +232,10 @@ __attribute__((always_inline)) INLINE static float mhd_signal_velocity(
const struct part *restrict pj, const float mu_ij, const float beta,
const float a, const float mu_0) {
- const float v_sigi =
- mhd_get_magnetosonic_speed(pi, a, mu_0);
- const float v_sigj =
- mhd_get_magnetosonic_speed(pj, a, mu_0);
+ const float v_sigi = mhd_get_magnetosonic_speed(pi, a, mu_0);
+ const float v_sigj = mhd_get_magnetosonic_speed(pj, a, mu_0);
- /*
+ /*
const float v_sigi =
mhd_get_fast_magnetosonic_wave_phase_velocity(dx, pi, a, mu_0);
const float v_sigj =
diff --git a/src/mhd/DirectInduction/mhd_iact.h b/src/mhd/DirectInduction/mhd_iact.h
index ec101145d9..9530cab1ce 100644
--- a/src/mhd/DirectInduction/mhd_iact.h
+++ b/src/mhd/DirectInduction/mhd_iact.h
@@ -606,7 +606,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
const float corr_ratio_i = fabsf(normBi * psi_over_ch_i_inv);
const float corr_ratio_j = fabsf(normBj * psi_over_ch_j_inv);
-
+
const float Qi = corr_ratio_i < 2 ? 0.5 * corr_ratio_i : 1.0f;
const float Qj = corr_ratio_j < 2 ? 0.5 * corr_ratio_j : 1.0f;
@@ -897,7 +897,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
psi_over_ch_i != 0.f ? 1.f / psi_over_ch_i : 0.;
const float corr_ratio_i = fabsf(normBi * psi_over_ch_i_inv);
-
+
const float Qi = corr_ratio_i < 2 ? 0.5 * corr_ratio_i : 1.0f;
pi->mhd_data.B_over_rho_dt[0] -= mj * Qi * grad_psi_i * dx[0];
diff --git a/src/mhd/VPotential/mhd.h b/src/mhd/VPotential/mhd.h
index 208be1925b..4c5381ece0 100644
--- a/src/mhd/VPotential/mhd.h
+++ b/src/mhd/VPotential/mhd.h
@@ -121,7 +121,8 @@ __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 = cosmo->a/cosmo->a_factor_sound_speed ;//pow(cosmo->a, -mhd_comoving_factor - 0.5f);
+ // const float afac_divB = cosmo->a/cosmo->a_factor_sound_speed
+ // ;//pow(cosmo->a, -mhd_comoving_factor - 0.5f);
const float afac_divB = pow(cosmo->a, -mhd_comoving_factor - 0.5f);
const float afac_resistive = cosmo->a * cosmo->a;
@@ -132,10 +133,11 @@ __attribute__((always_inline)) INLINE static float mhd_compute_timestep(
? afac_divB * hydro_properties->CFL_condition *
sqrtf(p->rho / (p->mhd_data.divB * p->mhd_data.divB) * mu_0)
: FLT_MAX;
- //dt_divB *= afac_divB;
- //const float resistive_eta = p->mhd_data.resistive_eta + (mhd_comoving_factor + 2.f) * afac_resistive * p->h * p->h * cosmo->H;
+ // dt_divB *= afac_divB;
+ // const float resistive_eta = p->mhd_data.resistive_eta +
+ // (mhd_comoving_factor + 2.f) * afac_resistive * p->h * p->h * cosmo->H;
const float resistive_eta = p->mhd_data.resistive_eta;
- //const float resistive_eta = p->mhd_data.resistive_eta;
+ // const float resistive_eta = p->mhd_data.resistive_eta;
const float dt_eta = resistive_eta != 0.0f
? afac_resistive * hydro_properties->CFL_condition *
p->h * p->h / resistive_eta * 0.5
@@ -273,21 +275,22 @@ __attribute__((always_inline)) INLINE static float mhd_signal_velocity(
* @param Gauge Gauge
*/
__attribute__((always_inline)) INLINE static float hydro_get_dGau_dt(
- const struct part *restrict p, const float Gauge,
- const struct cosmology *c, const float mu0) {
+ const struct part *restrict p, const float Gauge, const struct cosmology *c,
+ const float mu0) {
- //const float v_sig = hydro_get_signal_velocity(p);
- const float v_sig = mhd_get_magnetosonic_speed(p,c->a,mu0);
- //const float afac1 = pow(c->a, 2.f * mhd_comoving_factor - 1.f);
- //const float afac1 = pow(c->a, 2.f * mhd_comoving_factor - 2.f);
+ // const float v_sig = hydro_get_signal_velocity(p);
+ const float v_sig = mhd_get_magnetosonic_speed(p, c->a, mu0);
+ // const float afac1 = pow(c->a, 2.f * mhd_comoving_factor - 1.f);
+ // const float afac1 = pow(c->a, 2.f * mhd_comoving_factor - 2.f);
const float afac1 = pow(c->a, 2.f * c->a_factor_sound_speed);
- //const float afac2 = pow(c->a, mhd_comoving_factor + 1.f);
- //const float afac2 = pow(c->a, 2.f * mhd_comoving_factor + 1.5f);
+ // const float afac2 = pow(c->a, mhd_comoving_factor + 1.f);
+ // const float afac2 = pow(c->a, 2.f * mhd_comoving_factor + 1.5f);
const float afac2 = pow(c->a, c->a_factor_sound_speed + 1.f);
return (-p->mhd_data.divA * v_sig * v_sig * 0.1 * afac1 -
- 1.0f * v_sig * Gauge / p->h * afac2
- - (2.f + mhd_comoving_factor) * c->H * Gauge) * c->a * c->a;
+ 1.0f * v_sig * Gauge / p->h * afac2 -
+ (2.f + mhd_comoving_factor) * c->H * Gauge) *
+ c->a * c->a;
}
/**
@@ -322,14 +325,14 @@ __attribute__((always_inline)) INLINE static void mhd_init_part(
__attribute__((always_inline)) INLINE static void mhd_end_density(
struct part *p, const struct cosmology *cosmo) {
- const float h_inv_dim_plus_one = pow_dimension(1.f / p->h) /p->h; /*1/h^(d+1) */
+ const float h_inv_dim_plus_one =
+ pow_dimension(1.f / p->h) / 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;
for (int i = 0; i < 3; i++)
p->mhd_data.BPred[i] *= h_inv_dim_plus_one * rho_inv;
}
-
/**
* @brief Prepare a particle for the gradient calculation.
*
@@ -433,16 +436,17 @@ __attribute__((always_inline)) INLINE static void mhd_prepare_force(
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);
- // this should go with a /p->h, but I
- // take simplify becasue of ACC_mhd also.
+ // 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;
- // p->mhd_data.Q0 = 0.0f;
+ // p->mhd_data.Q0 =
+ // ACC_corr > ACC_mhd ? p->mhd_data.Q0 * ACC_mhd / ACC_corr :
+ // p->mhd_data.Q0;
+ // p->mhd_data.Q0 = 0.0f;
}
/**
@@ -470,8 +474,7 @@ __attribute__((always_inline)) INLINE static void mhd_reset_acceleration(
* @param cosmo The cosmological model
*/
__attribute__((always_inline)) INLINE static void mhd_reset_predicted_values(
- struct part *p, const struct xpart *xp, const struct cosmology *cosmo) {
-}
+ struct part *p, const struct xpart *xp, const struct cosmology *cosmo) {}
/**
* @brief Predict additional particle fields forward in time when drifting
@@ -499,7 +502,8 @@ __attribute__((always_inline)) INLINE static void mhd_predict_extra(
p->mhd_data.APred[1] += p->mhd_data.dAdt[1] * dt_therm;
p->mhd_data.APred[2] += p->mhd_data.dAdt[2] * dt_therm;
- float change_Gau = hydro_get_dGau_dt(p, p->mhd_data.Gau,cosmo,mu_0) * dt_therm;
+ float change_Gau =
+ hydro_get_dGau_dt(p, p->mhd_data.Gau, cosmo, mu_0) * dt_therm;
p->mhd_data.Gau += change_Gau;
}
@@ -516,9 +520,9 @@ __attribute__((always_inline)) INLINE static void mhd_predict_extra(
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void mhd_end_force(
- struct part *p, const struct cosmology *cosmo,
+ struct part *p, const struct cosmology *cosmo,
const struct hydro_props *hydro_props, const float mu_0) {
-
+
float a_fac = (1.f + mhd_comoving_factor) * cosmo->a * cosmo->a * cosmo->H;
p->mhd_data.dAdt[0] -= a_fac * p->mhd_data.APred[0];
p->mhd_data.dAdt[1] -= a_fac * p->mhd_data.APred[1];
@@ -545,10 +549,9 @@ __attribute__((always_inline)) INLINE static void mhd_kick_extra(
struct part *p, struct xpart *xp, const float dt_therm, const float dt_grav,
const float dt_hydro, const float dt_kick_corr,
const struct cosmology *cosmo, const struct hydro_props *hydro_props,
- const struct entropy_floor_properties *floor_props) {
-}
+ const struct entropy_floor_properties *floor_props) {}
- /**
+/**
* @brief Converts MHD quantities of a particle at the start of a run
*
* This function is called once at the end of the engine_init_particle()
@@ -568,10 +571,9 @@ __attribute__((always_inline)) INLINE static void mhd_convert_quantities(
/* Set Restitivity Eta */
p->mhd_data.resistive_eta = hydro_props->mhd.mhd_eta;
- p->mhd_data.APred[0] *= pow(cosmo->a, -mhd_comoving_factor -1.f);
- p->mhd_data.APred[1] *= pow(cosmo->a, -mhd_comoving_factor -1.f);
- p->mhd_data.APred[2] *= pow(cosmo->a, -mhd_comoving_factor -1.f);
-
+ p->mhd_data.APred[0] *= pow(cosmo->a, -mhd_comoving_factor - 1.f);
+ p->mhd_data.APred[1] *= pow(cosmo->a, -mhd_comoving_factor - 1.f);
+ p->mhd_data.APred[2] *= pow(cosmo->a, -mhd_comoving_factor - 1.f);
}
/**
diff --git a/src/mhd/VPotential/mhd_iact.h b/src/mhd/VPotential/mhd_iact.h
index c577322e3b..0988871c7e 100644
--- a/src/mhd/VPotential/mhd_iact.h
+++ b/src/mhd/VPotential/mhd_iact.h
@@ -395,8 +395,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
SAi = SourceAi + a * a * (pi->mhd_data.Gau - pj->mhd_data.Gau);
SAj = SourceAj + a * a * (pi->mhd_data.Gau - pj->mhd_data.Gau);
- //SAi = SourceAi - (pi->mhd_data.Gau - pj->mhd_data.Gau);
- //SAj = SourceAj - (pi->mhd_data.Gau - pj->mhd_data.Gau);
+ // SAi = SourceAi - (pi->mhd_data.Gau - pj->mhd_data.Gau);
+ // SAj = SourceAj - (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];
@@ -410,10 +410,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_mhd_force(
for (int i = 0; i < 3; i++) {
pi->mhd_data.dAdt[i] +=
mj * 8.0 * (pi->mhd_data.resistive_eta) * mag_Disi * dA[i];
- //mj * 8.0 * (pi->mhd_data.resistive_eta + ( mhd_comoving_factor + 2.f ) * a * a * hi * hi * H) * mag_Disi * dA[i];
+ // mj * 8.0 * (pi->mhd_data.resistive_eta + ( mhd_comoving_factor + 2.f ) *
+ // a * a * hi * hi * H) * mag_Disi * dA[i];
pj->mhd_data.dAdt[i] -=
mi * 8.0 * (pj->mhd_data.resistive_eta) * mag_Disj * dA[i];
- //mi * 8.0 * (pj->mhd_data.resistive_eta + ( mhd_comoving_factor + 2.f ) * a * a * hj * hj * H) * mag_Disj * dA[i];
+ // mi * 8.0 * (pj->mhd_data.resistive_eta + ( mhd_comoving_factor + 2.f ) *
+ // a * a * hj * hj * H) * mag_Disj * dA[i];
}
}
@@ -524,7 +526,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
// float SAi = SourceAi + sourceAi;
float SAi = (SourceAi + sourceAi) / 2.f;
- SAi = SourceAi + a * a * (pi->mhd_data.Gau - pj->mhd_data.Gau);
+ SAi = SourceAi + 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];
/// DISSSIPATION
@@ -533,6 +535,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_mhd_force(
for (int i = 0; i < 3; i++)
pi->mhd_data.dAdt[i] +=
mj * 8.0 * (pi->mhd_data.resistive_eta) * mag_Disi * dA[i];
- //mj * 8.0 * (pi->mhd_data.resistive_eta + 0.f * ( mhd_comoving_factor + 2.f ) * a * a * hi * hi * H) * mag_Disi * dA[i];
+ // mj * 8.0 * (pi->mhd_data.resistive_eta + 0.f * ( mhd_comoving_factor + 2.f
+ // ) * a * a * hi * hi * H) * mag_Disi * dA[i];
}
#endif /* SWIFT_VECTOR_POTENTIAL_MHD_H */
diff --git a/src/mhd/VPotential/mhd_io.h b/src/mhd/VPotential/mhd_io.h
index 438dd13fbf..5edd4913b3 100644
--- a/src/mhd/VPotential/mhd_io.h
+++ b/src/mhd/VPotential/mhd_io.h
@@ -76,11 +76,12 @@ INLINE static int mhd_write_particles(const struct part* parts,
parts, mhd_data.APred,
"Co-moving Magnetic Vector Potentials of the particles");
- list[3] = io_make_output_field("VectorPotentialScalarGauges", FLOAT, 1,
- UNIT_CONV_MAGNETIC_FIELD_VECTOR_POTENTIAL_GAUGE,
- mhd_comoving_factor + 2.f, parts, mhd_data.Gau,
- "Co-coving gauge scalar associated to the "
- "magnetic vector potentials per particle");
+ list[3] =
+ io_make_output_field("VectorPotentialScalarGauges", FLOAT, 1,
+ UNIT_CONV_MAGNETIC_FIELD_VECTOR_POTENTIAL_GAUGE,
+ mhd_comoving_factor + 2.f, parts, mhd_data.Gau,
+ "Co-coving gauge scalar associated to the "
+ "magnetic vector potentials per particle");
list[4] = io_make_output_field(
"VectorPotentialDivergences", FLOAT, 1, UNIT_CONV_MAGNETIC_FIELD,
diff --git a/src/units.c b/src/units.c
index f63cd8930f..773163bc7b 100644
--- a/src/units.c
+++ b/src/units.c
@@ -406,7 +406,7 @@ void units_get_base_unit_exponents_array(float baseUnitsExp[5],
baseUnitsExp[UNIT_TIME] = -2.f;
baseUnitsExp[UNIT_CURRENT] = -1.f;
break;
-
+
case UNIT_CONV_MAGNETIC_FIELD_VECTOR_POTENTIAL_GAUGE:
baseUnitsExp[UNIT_MASS] = 1.f;
baseUnitsExp[UNIT_LENGTH] = 2.f;
--
GitLab