improved harris test and post treatment

.gitignore

@@ -24,6 +24,7 @@ whislerwaveres/
 whislerwavelin/
 whislerwaveHLL/
 hallharrisres/
+hallharrisres2/
 subprojects
 
 debug/

CMakeLists.txt

@@ -35,22 +35,4 @@ set(SOURCES
     src/RiemannSolver.cpp
     src/TimeIntegrator.cpp)
 
-
-
-pybind11_add_module(pyMHD pyMHD/pyMHD.cpp
-    src/ConservativeVariables.cpp 
-    src/ConstrainedTransport.cpp
-    src/EquationOfState.cpp
-    src/GodunovFlux.cpp
-    src/Interface.cpp
-    src/ModularityUtils.cpp
-    src/PhareMHD.cpp
-    src/PhysicalConstants.cpp
-    src/PrimitiveVariables.cpp
-    src/RiemannSolver.cpp
-    src/TimeIntegrator.cpp)
-
-
-
-
-
+pybind11_add_module(pyMHD pyMHD/pyMHD.cpp ${SOURCES})

diagnostics/Hallharris/colormesh.py

@@ -0,0 +1,55 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+def read_data(file_path):
+    data = np.loadtxt(file_path)
+    return data
+
+def reshape_data(data, ny, nx):
+    reshaped_data = data.reshape((ny, nx, -1), order='F')
+    return reshaped_data
+
+nx = 250
+ny = 250
+
+file_path = "hallharrisres/PRK2_15_0001697500.txt"
+
+data = read_data(file_path)
+
+reshaped_data = reshape_data(data, nx, ny)
+
+dx = 0.1
+dy = 0.1
+
+# Extract Bx and By
+Bx = reshaped_data[:, :, 4]
+By = reshaped_data[:, :, 5]
+
+# Calculate the derivatives
+dBy_dx = np.gradient(By, dx, axis=0)
+dBx_dy = np.gradient(Bx, dy, axis=1)
+
+# Calculate Jz
+Jz = (dBy_dx - dBx_dy)
+
+plt.pcolormesh(Jz.T, cmap='coolwarm')  
+plt.title('Contour Plot of Jz at t=30')
+plt.xlabel('X')
+plt.ylabel('Y')
+plt.show()
+
+x_middle = nx // 2
+
+Jz_cutx = Jz[x_middle, :]
+
+y_positions = np.arange(ny)
+
+# Create the plot
+plt.figure(figsize=(10, 6))
+plt.plot(y_positions, Jz_cutx, marker='o')
+plt.title(f'1D cut of Jz at X = {x_middle/nx}')
+plt.xlabel('y')
+plt.ylabel('Jz')
+plt.grid(True)
+
+plt.show()

diagnostics/Hallharris/colormesh_anim.py

@@ -30,21 +30,36 @@ def read_times(file_paths):
 nx = 250
 ny = 250
 
+studied_index = 0
+
 data = [read_data(file_path) for file_path in file_paths]
 times = read_times(file_paths)
 
 reshaped_data = [reshape_data(d, nx, ny) for d in data]
 
+dx = 0.1
+dy = 0.1
 
-studied_index = 0
+# Extract Bx and By
+Bx = [reshaped_data[i][:, :, 4] for i in range(len(data))]
+By = [reshaped_data[i][:, :, 5] for i in range(len(data))]
+
+# Calculate the derivatives
+dBy_dx = [np.gradient(By[i], dx, axis=0) for i in range(len(data))]
+dBx_dy = [np.gradient(Bx[i], dy, axis=1) for i in range(len(data))]
+
+# Calculate Jz
+Jz = [(dBy_dx[i] - dBx_dy[i]) for i in range(len(data))]
+
+toPlot = Jz
 
-data_min = np.min(reshaped_data[-1][:, :, studied_index])
-data_max = np.max(reshaped_data[-1][:, :, studied_index])
+data_min = np.min(toPlot[-1])
+data_max = np.max(toPlot[-1])
 
 Norm = Normalize(vmin=data_min, vmax=data_max)
 
 fig, ax = plt.subplots()
-im = ax.pcolormesh(reshaped_data[0][:, :, studied_index].T, cmap='coolwarm', norm=Norm) 
+im = ax.pcolormesh(toPlot[0].T, cmap='coolwarm', norm=Norm) 
 ax.set_aspect('equal')
 fig.colorbar(im, ax=ax)
 
@@ -54,18 +69,18 @@ def read_times(file_paths):
 os.makedirs(output_dir, exist_ok=True)
 
 def update(frame):
-    im.set_array(reshaped_data[frame][:, :, studied_index].T)
-    plt.title(f'rho at t={times[frame]}')
-    #plt.savefig(f'{output_dir}/frame_{frame:04d}.png')
+    im.set_array(toPlot[frame].T)
+    plt.title(f'Jz at t={times[frame]}')
+    plt.savefig(f'{output_dir}/frame_{frame:04d}.png')
     return im,
 
 ani = FuncAnimation(fig, update, frames=len(times), interval=100)
 
 #gif_writer = PillowWriter(fps=10)  # Adjust fps as needed
-#ani.save('hallharris.gif', writer=gif_writer)
+#ani.save('orszagtangP.gif', writer=gif_writer)
 
 plt.xlabel('X')
 plt.ylabel('Y')
 plt.show()
 
-#ffmpeg -r 10 -i frames/frame_%04d.png -vcodec mpeg4 -q:v 5 hallharris.mp4
+#ffmpeg -r 10 -i frames/frame_%04d.png -vcodec mpeg4 -q:v 5 orszagtangP.mp4

diagnostics/Hallharris/hallharris.py

@@ -27,7 +27,7 @@
 OptionalPhysics = p.OptionalPhysics.HallResHyper
 
 dumpvariables = p.dumpVariables.Primitive
-dumpfrequency = 10
+dumpfrequency = 1000
 
 ##############################################################################################################################################################################
 Lx = nx*Dx
@@ -104,7 +104,7 @@ def P_(x, y):
 
 #############################################################################################################################################################################
 
-result_dir = 'hallharrisres/'
+result_dir = 'hallharrisres2/'
 if os.path.exists(result_dir):
     shutil.rmtree(result_dir)
 

diagnostics/harris/colormesh_anim.py

@@ -27,24 +27,39 @@ def read_times(file_paths):
 results_dir = "harrisres/"
 file_paths = [results_dir + file for file in os.listdir(results_dir) if file.startswith("PRK2_") and file.endswith(".txt")]
 
-nx = 100
-ny = 100
+nx = 250
+ny = 250
+
+studied_index = 0
 
 data = [read_data(file_path) for file_path in file_paths]
 times = read_times(file_paths)
 
 reshaped_data = [reshape_data(d, nx, ny) for d in data]
 
+dx = 0.1
+dy = 0.1
 
-studied_index = 0
+# Extract Bx and By
+Bx = [reshaped_data[i][:, :, 4] for i in range(len(data))]
+By = [reshaped_data[i][:, :, 5] for i in range(len(data))]
+
+# Calculate the derivatives
+dBy_dx = [np.gradient(By[i], dx, axis=0) for i in range(len(data))]
+dBx_dy = [np.gradient(Bx[i], dy, axis=1) for i in range(len(data))]
+
+# Calculate Jz
+Jz = [(dBy_dx[i] - dBx_dy[i]) for i in range(len(data))]
+
+toPlot = Jz
 
-data_min = np.min(reshaped_data[-1][:, :, studied_index])
-data_max = np.max(reshaped_data[-1][:, :, studied_index])
+data_min = np.min(toPlot[-1])
+data_max = np.max(toPlot[-1])
 
-Norm = Normalize(vmin=0, vmax=2.5)
+Norm = Normalize(vmin=data_min, vmax=data_max)
 
 fig, ax = plt.subplots()
-im = ax.pcolormesh(reshaped_data[0][:, :, studied_index].T, cmap='coolwarm', norm=Norm) 
+im = ax.pcolormesh(toPlot[0].T, cmap='coolwarm', norm=Norm) 
 ax.set_aspect('equal')
 fig.colorbar(im, ax=ax)
 
@@ -54,8 +69,8 @@ def read_times(file_paths):
 os.makedirs(output_dir, exist_ok=True)
 
 def update(frame):
-    im.set_array(reshaped_data[frame][:, :, studied_index].T)
-    plt.title(f'rho at t={times[frame]}')
+    im.set_array(toPlot[frame].T)
+    plt.title(f'{toPlot} at t={times[frame]}')
     plt.savefig(f'{output_dir}/frame_{frame:04d}.png')
     return im,
 

diagnostics/harris/harris.py

@@ -8,12 +8,12 @@
 
 #############################################################################################################################################################################
 
-nx = 100
-ny = 100
-Dx = 0.2
-Dy = 0.2
+nx = 250
+ny = 250
+Dx = 0.1
+Dy = 0.1
 Dt = 0.0
-FinalTime = 5
+FinalTime = 30
 nghost = 2
 
 boundaryconditions = p.BoundaryConditions.Periodic
@@ -80,7 +80,7 @@ def Bz_(x, y):
 
 
 def P_(x, y):
-    return 5./(12.*np.pi)
+    return 1.0 - (Bx_(x, y) ** 2)/2.0
 
 x = np.arange(nx) * Dx + 0.5 * Dx
 y = np.arange(ny) * Dy + 0.5 * Dy