added arithmetic constrained transport, tested

.gitignore

@@ -11,23 +11,14 @@ space_results/
 time_results/
 2Dwave/
 orszagtangres/
-orszagtangCTAverage/
-orszagtangUCTHLL/
-orszagtangCTContact/
-orszagtang1024
-orszagtang1024r2
+orszagtang*/
 MHDrotorres/
 shockres/
 harrisres/
 currentres/
-whislerwaveres/
-whislerwavelin/
-whislerwaveHLL/
-hallharrisres/
-hallharrisres2/
-hallharrisresCT/
-hallharrisresCT2/
-hallharrisresCT2long/
+whislerwaveres*/
+ideal*/
+hallharrisres*/
 subprojects
 
 debug/

diagnostics/Hallharris/colormesh2.py

@@ -13,7 +13,7 @@ def reshape_data(data, ny, nx):
 ny = 250
 
 file_path1 = "hallharrisresCT/PRK2_0_5002794389.txt"
-file_path2 = "hallharrisresCT2/PRK2_0_5000859584.txt"
+file_path2 = "hallharrisresCT2/PRK2_0_5001043321.txt"
 # file_path1 = "hallharrisres2/PRK2_2_4174122015.txt"
 # file_path2 = "hallharrisresCT2long/PRK2_2_4039801735.txt"
 

diagnostics/Hallharris/hallharris.py

@@ -21,7 +21,7 @@
 reconstruction = p.Reconstruction.Linear
 slopelimiter = p.Slope.MinMod
 riemannsolver = p.RiemannSolver.Rusanov
-constainedtransport = p.CTMethod.Average
+constainedtransport = p.CTMethod.Arithmetic
 timeintegrator = p.Integrator.TVDRK2Integrator
 
 OptionalPhysics = p.OptionalPhysics.HallResHyper

diagnostics/MHDrotor/MHDrotor.py

@@ -21,7 +21,7 @@
 reconstruction = p.Reconstruction.Constant
 slopelimiter = p.Slope.MinMod
 riemannsolver = p.RiemannSolver.Rusanov
-constainedtransport = p.CTMethod.Average
+constainedtransport = p.CTMethod.Arithmetic
 timeintegrator = p.Integrator.TVDRK2Integrator
 
 consts = p.Consts(sigmaCFL = 0.2, gam = 1.4)

diagnostics/alfvenwave2D/alfvenwave2D.py

@@ -1,6 +1,6 @@
 import sys
 
-# sys.path.append('/home/caromel/Documents/PHAREMHD/pyMHD')
+sys.path.append('/home/caromel/Documents/PHAREMHD/pyMHD')
 
 import numpy as np
 import pyMHD as p
@@ -22,7 +22,7 @@
 reconstruction = p.Reconstruction.Linear
 slopelimiter = p.Slope.VanLeer
 riemannsolver = p.RiemannSolver.Rusanov
-constainedtransport = p.CTMethod.Average
+constainedtransport = p.CTMethod.Arithmetic
 timeintegrator = p.Integrator.TVDRK2Integrator
 
 dumpvariables = p.dumpVariables.Conservative

diagnostics/currentsheet/colormesh_anim.py

@@ -29,14 +29,13 @@ def read_times(file_paths):
 nx = 200
 ny = 200
 
+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]
 
-
-studied_index = 0
-
 data_min = np.min(reshaped_data[-1][:, :, studied_index])
 data_max = np.max(reshaped_data[-1][:, :, studied_index])
 

diagnostics/currentsheet/current.py

@@ -18,17 +18,17 @@
 
 boundaryconditions = p.BoundaryConditions.Periodic
 
-reconstruction = p.Reconstruction.Constant
+reconstruction = p.Reconstruction.Linear
 slopelimiter = p.Slope.VanLeer
 riemannsolver = p.RiemannSolver.Rusanov
-constainedtransport = p.CTMethod.Average
+constainedtransport = p.CTMethod.Arithmetic
 timeintegrator = p.Integrator.TVDRK2Integrator
 
 dumpvariables = p.dumpVariables.Primitive
 dumpfrequency = 5
 
 ##############################################################################################################################################################################
-u0 = 0.001
+u0 = 0.1
 B0 = 1
 
 def rho_(x, y):

diagnostics/orszag-tang/orszagtang.py

@@ -8,8 +8,8 @@
 
 #############################################################################################################################################################################
 
-nx = 1024
-ny = 1024
+nx = 128
+ny = 128
 Dx = 1/nx
 Dy = 1/ny
 Dt = 0.0
@@ -18,10 +18,10 @@
 
 boundaryconditions = p.BoundaryConditions.Periodic
 
-reconstruction = p.Reconstruction.Constant
+reconstruction = p.Reconstruction.Linear
 slopelimiter = p.Slope.VanLeer
 riemannsolver = p.RiemannSolver.Rusanov
-constainedtransport = p.CTMethod.Average
+constainedtransport = p.CTMethod.Arithmetic
 timeintegrator = p.Integrator.TVDRK2Integrator
 
 dumpvariables = p.dumpVariables.Primitive

diagnostics/whistlerwave/fft.py

@@ -6,23 +6,27 @@
 from matplotlib.colors import LogNorm, Normalize, PowerNorm
 
 nx = 128
-Dx = 0.1
+Dx = 0.8
 
-Dt = 0.002
+Dt = 0.077
+
+modes = [1,2,4,8]
 
 lx=nx*Dx
-m = 4
+
+m = 1
+
 kt=2*np.pi/lx * m
-w = (kt**2 /2) *(np.sqrt(1+4/kt**2) + 1)
-finalTime = 2*np.pi / w
+w = kt#(kt**2 /2) *(np.sqrt(1+4/kt**2) + 1)
+finalTime = 2*np.pi / w *10
 
 column_names = ['rho', 'vx', 'vy', 'vz', 'Bx', 'By', 'Bz', 'P']
 
 def read_file(filename):
     df = pd.read_csv(filename, delim_whitespace=True, header=None, names=column_names)
     return df
 
-results_dir = 'whislerwaveres'
+results_dir = 'whislerwaveres2'
 
 quantities = {
     'rho': [],
@@ -59,10 +63,9 @@ def leftAlfvenIonCyclotron(k):
 kplus = k[:halfk]
 wplus = w[:halfw]
 
-modes = [m]
 kmodes = 2*np.pi*np.asarray([1/(nx*Dx) * m for m in modes])
 
-B_combined = quantities['By'] + 1j * quantities['Bz']
+B_combined = quantities['Bz']# + 1j * quantities['Bz']
 
 # window_space = np.hanning(nx)
 # window_time = np.hanning(int(finalTime / Dt) + 1)
@@ -75,14 +78,31 @@ def leftAlfvenIonCyclotron(k):
 
 Xx, Yy = np.meshgrid(kplus, wplus)
 
+def find_peak_frequency(k_index, fft_data):
+    frequency_slice = fft_data[:, k_index]
+    peak_index = np.argmax(frequency_slice)
+    return wplus[peak_index]
+
 fig, ax = plt.subplots()
 pcm = ax.pcolormesh(Xx, Yy, B_fft_abs, cmap='plasma', shading='auto', norm=PowerNorm(gamma=0.3, vmin=B_fft_abs.min(), vmax=B_fft_abs.max()))
 fig.colorbar(pcm, ax=ax, label='Magnitude of FFT(By + i*Bz)')
 ax.plot(kplus, whistler(kplus), marker = '+', color='k')
 ax.plot(kplus, leftAlfvenIonCyclotron(kplus), marker = '*', color='g')
 ax.plot(kplus, kplus)
-for km in kmodes:
-    ax.axvline(km, ls='--', color='k')
-    ax.plot(km, whistler(km), marker='o')
-
-plt.show()
+excited_modes = []
+for m in modes:
+    km = 2 * np.pi / lx * m
+    k_index = np.argmin(np.abs(kplus - km))
+    peak_w = find_peak_frequency(k_index, B_fft_abs)
+    ax.plot(kplus[k_index], peak_w, 'ro', markersize=10)
+    ax.text(kplus[k_index], peak_w, f'm={m}', fontsize=8, ha='right', va='bottom')
+    excited_modes.append((kplus[k_index], peak_w))
+
+ax.set_xlabel('k')
+ax.set_ylabel('ω')
+ax.set_title('2D FFT with Excited Modes')
+plt.show()
+
+print("Excited modes (k, ω):")
+for m, (k, w) in zip(modes, excited_modes):
+    print(f"m = {m}: ({k:.4f}, {w:.4f})")

diagnostics/whistlerwave/plot_anim.py

@@ -6,10 +6,10 @@
 import shutil
 
 nx = 128
-Dx = 0.1
+Dx = 0.05
 
 
-quantity_name = 'By'
+quantity_name = 'Bz'
 fixed_index = 0
 ny = 1
 
@@ -19,8 +19,8 @@ def read_file(filename):
     df = pd.read_csv(filename, delim_whitespace=True, header=None, names=column_names)
     return df
 
-results_dir = 'whislerwaveres/'
-results_dir2 = 'whislerwaveHLL/'
+results_dir = 'AICres2/'
+results_dir2 = 'whislerwaveres2/'
 
 quantities = {
     'rho': [],
@@ -34,7 +34,7 @@ def read_file(filename):
 }
 times = []
 
-"""quantitieshll = {
+quantities2 = {
     'rho': [],
     'vx': [],
     'vy': [],
@@ -44,7 +44,7 @@ def read_file(filename):
     'Bz': [],
     'P': []
 }
-timeshll = []"""
+times2 = []
 
 for filename in os.listdir(results_dir):
     if filename.startswith("PRK2_") and filename.endswith(".txt"):
@@ -60,19 +60,19 @@ def read_file(filename):
 for quantity in quantities.keys():
     quantities[quantity] = np.array(quantities[quantity])
 
-"""for filename in os.listdir(results_dir2):
-    if filename.startswith("PRK2_") and filename.endswith(".txt"):
-        time_str = filename.split('_')[1]+'.'+filename.split('_')[2].split('.')[0]
-        time = float(time_str)
-        timeshll.append(time)
+# for filename in os.listdir(results_dir2):
+#     if filename.startswith("PRK2_") and filename.endswith(".txt"):
+#         time_str = filename.split('_')[1]+'.'+filename.split('_')[2].split('.')[0]
+#         time = float(time_str)
+#         times2.append(time)
 
-        df = read_file(os.path.join(results_dir2, filename))
+#         df = read_file(os.path.join(results_dir2, filename))
 
-        for quantity in quantitieshll.keys():
-            quantitieshll[quantity].append(df[quantity].values.reshape((ny, nx)))
+#         for quantity in quantities2.keys():
+#             quantities2[quantity].append(df[quantity].values.reshape((ny, nx)))
 
-for quantity in quantitieshll.keys():
-    quantitieshll[quantity] = np.array(quantitieshll[quantity])"""
+# for quantity in quantities2.keys():
+#     quantities2[quantity] = np.array(quantities2[quantity])
 
 x=Dx*np.arange(nx) + 0.5*Dx
 
@@ -83,17 +83,17 @@ def read_file(filename):
 
 def update(frame):
     lx = nx*Dx
-    m = 4#int(nx/4)
+    m = 12#int(nx/4)
 
     k = 2 * np.pi / lx * m
 
     w = (k**2 /2) *(np.sqrt(1+4/k**2) + 1)
 
-    expected_value = 1e-7 * np.cos(k * x - w * times[frame] + 0.5488135)
+    expected_value = -1 * np.sin(k * x - w * times[frame] + 0.5488135)
 
     plt.clf()
     plt.plot(x, quantities[quantity_name][frame, fixed_index, :], color='blue', marker = 'x', markersize=3) # t,y,x
-    #plt.plot(x, quantitieshll[quantity_name][frame, fixed_index, :], color='green', marker = 'o', markersize=3) # t,y,x
+    #plt.plot(x, quantities2[quantity_name][frame, fixed_index, :], color='green', marker = 'o', markersize=3) # t,y,x
     plt.plot(x, expected_value)
     plt.title(f'{quantity_name} at t={times[frame]}')  # Format time to one decimal place
     plt.xlabel('x')

diagnostics/whistlerwave/whislerwaveX.py

@@ -10,9 +10,9 @@
 
 nx = 128
 ny = 1
-Dx = 0.1
+Dx = 0.8
 Dy = 1
-Dt = 0.002
+Dt = 0.077
 
 nghost = 2
 
@@ -22,7 +22,7 @@
 
 slopelimiter = p.Slope.VanLeer
 riemannsolver = p.RiemannSolver.Rusanov
-constainedtransport = p.CTMethod.Average
+constainedtransport = p.CTMethod.NoCT
 timeintegrator = p.Integrator.EulerIntegrator
 
 consts = p.Consts(sigmaCFL = 0.8, gam = 5/3, eta = 0.0, nu = 0.000)
@@ -35,15 +35,15 @@
 lx=nx*Dx
 k=2*np.pi/lx
 
-m = 4
+m = 1
 
 kt=2*np.pi/lx * m
 w = (kt**2 /2) *(np.sqrt(1+4/kt**2) + 1)
-FinalTime = 2*np.pi / w
+FinalTime = 2*np.pi / w *10
 
 np.random.seed(0)
 
-modes = [m]#[int(nx/4)]
+modes = [1,2,4,8]
 phases = np.random.rand(len(modes))
 
 def rho_(x, y):
@@ -56,13 +56,13 @@ def vy_(x, y):
     ret = np.zeros((x.shape[0], y.shape[1]))
 
     for m,phi in zip(modes, phases):
-        ret[:,:] += -np.cos(k*x*m + phi)*1e-7*k
+        ret[:,:] += -np.cos(k*x*m + phi)*0.01*k
     return ret
 
 def vz_(x, y):
     ret = np.zeros((x.shape[0], y.shape[1]))
     for m,phi in zip(modes, phases):
-        ret[:,:] += np.sin(k*x*m + phi)*1e-7*k
+        ret[:,:] += np.sin(k*x*m + phi)*0.01*k
     return ret
 
 def Bx_(x, y):
@@ -71,13 +71,13 @@ def Bx_(x, y):
 def By_(x, y):
     ret = np.zeros((x.shape[0], y.shape[1]))
     for m,phi in zip(modes, phases):
-        ret[:,:] += np.cos(k*x*m + phi)*1e-7
+        ret[:,:] += np.cos(k*x*m + phi)*0.01
     return ret
 
 def Bz_(x, y):
     ret = np.zeros((x.shape[0], y.shape[1]))
     for m,phi in zip(modes, phases):
-        ret[:,:] += -np.sin(k*x*m + phi)*1e-7
+        ret[:,:] += -np.sin(k*x*m + phi)*0.01
     return ret
 
 def P_(x, y):
@@ -103,7 +103,7 @@ def P_(x, y):
 
 #############################################################################################################################################################################
 
-result_dir = 'whislerwaveres/'
+result_dir = 'idealres/'
 if os.path.exists(result_dir):
     shutil.rmtree(result_dir)