added NoCT for 1D

.gitignore

@@ -25,8 +25,12 @@ whislerwavelin/
 whislerwaveHLL/
 hallharrisres/
 hallharrisres2/
+hallharrisresCT/
+hallharrisresCT2/
+hallharrisresCT2long/
 subprojects
 
 debug/
 frames/
+framestest/
 pyMHD/pyMHD.cpython*

README.md

@@ -1,12 +1,5 @@
 ## TODO
-
-# By priority order in short term
-- Implement diagnostics and testing (the rough debugging is done).
-- Implement post-processing
-- Make the code more modular
-
-# And then
-- Extend the model (3D, Non-Ideal MHD, more sofisticated Riemann solvers, more sophisticated CT ...).
+- Extend the model (3D, more sofisticated Riemann solvers, more sophisticated CT ...).
 - Clean up the code.
 - Integrate in PHARE.
 - Add coupling (Fluid-Fluid, Fluid-Hybrid)

diagnostics/Hallharris/colormesh.py

@@ -12,7 +12,7 @@ def reshape_data(data, ny, nx):
 nx = 250
 ny = 250
 
-file_path = "hallharrisres/PRK2_15_0001697500.txt"
+file_path = "hallharrisresCT/PRK2_0_5002794389.txt"
 
 data = read_data(file_path)
 
@@ -22,8 +22,10 @@ def reshape_data(data, ny, nx):
 dy = 0.1
 
 # Extract Bx and By
+rho = reshaped_data[:, :, 0]
 Bx = reshaped_data[:, :, 4]
 By = reshaped_data[:, :, 5]
+Bz = reshaped_data[:, :, 6]
 
 # Calculate the derivatives
 dBy_dx = np.gradient(By, dx, axis=0)
@@ -32,24 +34,26 @@ def reshape_data(data, ny, nx):
 # Calculate Jz
 Jz = (dBy_dx - dBx_dy)
 
-plt.pcolormesh(Jz.T, cmap='coolwarm')  
-plt.title('Contour Plot of Jz at t=30')
+toPlot = Bz
+
+plt.pcolormesh(toPlot.T, cmap='coolwarm')  
+plt.title('Contour Plot of Bz at t=0.5')
 plt.xlabel('X')
 plt.ylabel('Y')
 plt.show()
 
-x_middle = nx // 2
+# x_middle = nx // 2
 
-Jz_cutx = Jz[x_middle, :]
+# Jz_cutx = Jz[x_middle, :]
 
-y_positions = np.arange(ny)
+# 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)
+# # 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()
+# plt.show()

diagnostics/Hallharris/colormesh2.py

@@ -0,0 +1,52 @@
+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_path1 = "hallharrisresCT/PRK2_0_5002794389.txt"
+file_path2 = "hallharrisresCT2/PRK2_0_5000859584.txt"
+# file_path1 = "hallharrisres2/PRK2_2_4174122015.txt"
+# file_path2 = "hallharrisresCT2long/PRK2_2_4039801735.txt"
+
+# Read and process data from the first file
+data1 = read_data(file_path1)
+reshaped_data1 = reshape_data(data1, ny, nx)
+
+# Extract Bz for the first file
+Bz1 = reshaped_data1[:, :, 6]
+
+# Read and process data from the second file
+data2 = read_data(file_path2)
+reshaped_data2 = reshape_data(data2, ny, nx)
+
+# Extract Bz for the second file
+Bz2 = reshaped_data2[:, :, 6]
+
+# Create subplots
+fig, axes = plt.subplots(1, 2, figsize=(14, 6))
+
+# Plot Bz for the first file
+im1 = axes[0].pcolormesh(Bz1.T, cmap='coolwarm')
+axes[0].set_title('Contour Plot of Bz at t=0.5 (NoCT)')
+axes[0].set_xlabel('X')
+axes[0].set_ylabel('Y')
+fig.colorbar(im1, ax=axes[0])
+
+# Plot Bz for the second file
+im2 = axes[1].pcolormesh(Bz2.T, cmap='coolwarm')
+axes[1].set_title('Contour Plot of Bz at t=0.5 (CT)')
+axes[1].set_xlabel('X')
+axes[1].set_ylabel('Y')
+fig.colorbar(im2, ax=axes[1])
+
+plt.tight_layout()
+plt.show()

diagnostics/Hallharris/colormesh_anim.py

@@ -24,7 +24,7 @@ def read_times(file_paths):
         times.append(time)
     return times
 
-results_dir = "hallharrisres/"
+results_dir = "hallharrisresCT2long/"
 file_paths = [results_dir + file for file in os.listdir(results_dir) if file.startswith("PRK2_") and file.endswith(".txt")]
 
 nx = 250
@@ -41,8 +41,10 @@ def read_times(file_paths):
 dy = 0.1
 
 # Extract Bx and By
+rho = [reshaped_data[i][:, :, 0] for i in range(len(data))]
 Bx = [reshaped_data[i][:, :, 4] for i in range(len(data))]
 By = [reshaped_data[i][:, :, 5] for i in range(len(data))]
+Bz = [reshaped_data[i][:, :, 6] for i in range(len(data))]
 
 # Calculate the derivatives
 dBy_dx = [np.gradient(By[i], dx, axis=0) for i in range(len(data))]
@@ -51,7 +53,7 @@ def read_times(file_paths):
 # Calculate Jz
 Jz = [(dBy_dx[i] - dBx_dy[i]) for i in range(len(data))]
 
-toPlot = Jz
+toPlot = Bz
 
 data_min = np.min(toPlot[-1])
 data_max = np.max(toPlot[-1])
@@ -70,7 +72,7 @@ def read_times(file_paths):
 
 def update(frame):
     im.set_array(toPlot[frame].T)
-    plt.title(f'Jz at t={times[frame]}')
+    plt.title(f'rho at t={times[frame]}')
     plt.savefig(f'{output_dir}/frame_{frame:04d}.png')
     return im,
 
@@ -83,4 +85,4 @@ def update(frame):
 plt.ylabel('Y')
 plt.show()
 
-#ffmpeg -r 10 -i frames/frame_%04d.png -vcodec mpeg4 -q:v 5 orszagtangP.mp4
+#ffmpeg -r 10 -i frames/frame_%04d.png -vcodec mpeg4 -q:v 5 hallharris.mp4

diagnostics/Hallharris/hallharris.py

@@ -13,7 +13,7 @@
 Dx = 0.1
 Dy = 0.1
 Dt = 0.0
-FinalTime = 30
+FinalTime = 0.5
 nghost = 2
 
 boundaryconditions = p.BoundaryConditions.Periodic
@@ -104,7 +104,7 @@ def P_(x, y):
 
 #############################################################################################################################################################################
 
-result_dir = 'hallharrisres2/'
+result_dir = 'hallharrisresCT2/'
 if os.path.exists(result_dir):
     shutil.rmtree(result_dir)
 

diagnostics/MHDrotor/MHDrotor.py

@@ -28,12 +28,11 @@
 
 ##############################################################################################################################################################################
 B0 = 5/(np.sqrt(4*np.pi))
-v0 = 1
+v0 = 2
 
 r0 = 0.1
 r1 = 0.115
 
-# Domain [-0.5, 0.5]*[-0.5, 0.5]
 def r(x, y):
     return np.sqrt((x-0.5)**2 + (y-0.5)**2)
 
@@ -51,14 +50,14 @@ def vx_(x, y):
     r_ = r(x, y)
     f_ = f(r_)
 
-    vx_values = np.where(r_ <= r0, -f_ * v0 * (y - 0.5) / r0, np.where(r_ < r1, -f_ * v0 * (y - 0.5) / r_, 0.0))
+    vx_values = np.where(r_ <= r0, - v0 * (y - 0.5) / r0, np.where(r_ < r1, -f_ * v0 * (y - 0.5) / r_, 0.0))
     return vx_values
 
 def vy_(x, y):
     r_ = r(x, y)
     f_ = f(r_)
 
-    vy_values = np.where(r_ <= r0, f_ * v0 * (x - 0.5) / r0, np.where(r_ < r1, f_ * v0 * (x - 0.5) / r_, 0.0))
+    vy_values = np.where(r_ <= r0, v0 * (x - 0.5) / r0, np.where(r_ < r1, f_ * v0 * (x - 0.5) / r_, 0.0))
     return vy_values
 
 def vz_(x, y):
@@ -74,7 +73,7 @@ def Bz_(x, y):
     return 0.0 
 
 def P_(x, y):
-    return 1.0
+    return 0.5
 
 
 x = np.arange(nx) * Dx + 0.5 * Dx

diagnostics/MHDrotor/colormesh.py

@@ -0,0 +1,27 @@
+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 = 100
+ny = 100
+
+file_path = "MHDrotorres/URK2_0_0.txt"
+
+data = read_data(file_path)
+
+reshaped_data = reshape_data(data, nx, ny)
+
+rho = reshaped_data[:, :, 0]
+
+plt.pcolormesh(rho.T, cmap='coolwarm')  
+plt.title('Contour Plot of rho at t=0')
+plt.xlabel('X')
+plt.ylabel('Y')
+plt.show()

diagnostics/MHDrotor/colormesh_anim.py

@@ -25,8 +25,8 @@ def read_times(file_paths):
 results_dir = "MHDrotorres/"
 file_paths = [results_dir + file for file in os.listdir(results_dir) if file.startswith("URK2_") and file.endswith(".txt")]
 
-nx = 200
-ny = 200
+nx = 100
+ny = 100
 
 data = [read_data(file_path) for file_path in file_paths]
 times = read_times(file_paths)

diagnostics/orszag-tang/colormesh_anim.py

@@ -2,6 +2,7 @@
 import matplotlib.pyplot as plt
 import os
 from matplotlib.colors import Normalize
+from matplotlib.animation import FuncAnimation
 import shutil
 
 # Function to read data from file
@@ -59,8 +60,7 @@ def update(frame):
     plt.savefig(f'{output_dir}/frame_{frame:04d}.png')
     return im,
 
-for frame in range(len(times)):
-    update(frame)
+ani = FuncAnimation(fig, update, frames=len(times), interval=100)
 
 plt.xlabel('X')
 plt.ylabel('Y')

diagnostics/orszag-tang/orszagtang.py

@@ -8,26 +8,26 @@
 
 #############################################################################################################################################################################
 
-nx = 128
-ny = 128
+nx = 1024
+ny = 1024
 Dx = 1/nx
 Dy = 1/ny
 Dt = 0.0
-FinalTime = 0.5
-nghost = 1
+FinalTime = 1
+nghost = 2
 
 boundaryconditions = p.BoundaryConditions.Periodic
 
 reconstruction = p.Reconstruction.Constant
 slopelimiter = p.Slope.VanLeer
 riemannsolver = p.RiemannSolver.Rusanov
 constainedtransport = p.CTMethod.Average
-timeintegrator = p.Integrator.EulerIntegrator
+timeintegrator = p.Integrator.TVDRK2Integrator
 
 dumpvariables = p.dumpVariables.Primitive
 dumpfrequency = 80
 
-OptionalPhysics = p.OptionalPhysics.HallResHyper
+OptionalPhysics = p.OptionalPhysics.Off
 
 ##############################################################################################################################################################################
 B0 = 1./(np.sqrt(4.*np.pi))

diagnostics/whistlerwave/fft.py

@@ -8,8 +8,13 @@
 nx = 128
 Dx = 0.1
 
-Dt = 0.0025
-finalTime = 1
+Dt = 0.002
+
+lx=nx*Dx
+m = 4
+kt=2*np.pi/lx * m
+w = (kt**2 /2) *(np.sqrt(1+4/kt**2) + 1)
+finalTime = 2*np.pi / w
 
 column_names = ['rho', 'vx', 'vy', 'vz', 'Bx', 'By', 'Bz', 'P']
 
@@ -54,21 +59,17 @@ def leftAlfvenIonCyclotron(k):
 kplus = k[:halfk]
 wplus = w[:halfw]
 
-"""modes = [4]
-kmodes = 2*np.pi*np.asarray([1/(nx*Dx) * m for m in modes])"""
+modes = [m]
+kmodes = 2*np.pi*np.asarray([1/(nx*Dx) * m for m in modes])
 
-"""
-fig,ax = plt.subplots()
-ax.plot(kplus, whistler(kplus))
-ax.plot(kplus, np.abs(fft(quantities['By'][-1,:] + 1j*fft(quantities['Bz'][-1,:]))[:halfk]))
-for km in kmodes:
-    ax.axvline(km, ls='--', color='k')
-    ax.plot(km, km**2, marker='o')
+B_combined = quantities['By'] + 1j * quantities['Bz']
 
-plt.show()
-"""
+# window_space = np.hanning(nx)
+# window_time = np.hanning(int(finalTime / Dt) + 1)
+# window = np.outer(window_time, window_space)
+
+# B_combined *= window
 
-B_combined = quantities['By'] + 1j * quantities['Bz']
 B_fft = fft2(B_combined)
 B_fft_abs = np.abs(B_fft)[:halfw,:halfk]
 
@@ -80,8 +81,8 @@ def leftAlfvenIonCyclotron(k):
 ax.plot(kplus, whistler(kplus), marker = '+', color='k')
 ax.plot(kplus, leftAlfvenIonCyclotron(kplus), marker = '*', color='g')
 ax.plot(kplus, kplus)
-"""for km in kmodes:
+for km in kmodes:
     ax.axvline(km, ls='--', color='k')
-    ax.plot(km, km**2, marker='o')"""
+    ax.plot(km, whistler(km), marker='o')
 
 plt.show()