progress on hall debugging

.gitignore

@@ -22,6 +22,7 @@ harrisres/
 currentres/
 whislerwaveres/
 whislerwavelin/
+whislerwaveHLL/
 
 frames/
 pyMHD/pyMHD.cpython*

diagnostics/whislerwave/plotJ.py

@@ -0,0 +1,70 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from matplotlib.animation import FuncAnimation
+import os
+
+
+Dx = 0.2
+
+Jz = []
+Jy = []
+times = []
+results_dir = 'whislerwaveres/'
+
+def read_file(filename):
+    df = pd.read_csv(filename, delim_whitespace=True, header=None)
+    return df
+
+for filename in os.listdir(results_dir):
+    if filename.startswith("Jy_") and filename.endswith(".txt"):
+        time_str = filename.split('_')[1].split('.')[0]
+        time = float(time_str)
+        times.append(time)
+
+        df = read_file(os.path.join(results_dir, filename))
+        Jy.append(df.values.reshape((5,37)))
+
+    if filename.startswith("Jz_") and filename.endswith(".txt"):
+
+        df = read_file(os.path.join(results_dir, filename))
+        Jz.append(df.values.reshape((6,37)))
+
+x=Dx*np.arange(37)-2*Dx
+
+def update(frame):
+    plt.clf()
+    plt.plot(x, Jy[frame][2,:], 'b-', marker = 'o')
+    plt.plot(x, Jz[frame][2,:], 'r--',marker = '*')
+    plt.xlabel('x')
+    plt.grid(True)
+    #plt.yscale("log")
+    plt.tight_layout()
+
+    eps = 0.01
+    min_val = np.min(Jy) - eps
+    max_val = np.max(Jz) + eps
+
+    plt.ylim(min_val, max_val)
+
+
+fig = plt.figure(figsize=(8, 6))
+ani = FuncAnimation(fig, update, frames=len(times), interval=500)
+
+# Define a function for manually stepping through the frames
+def onclick(event):
+    if event.key == 'right':
+        ani.frame_seq = ani.new_frame_seq()
+        fig.canvas.draw()
+
+# Connect the key press event to the function
+fig.canvas.mpl_connect('key_press_event', onclick)
+
+plt.show()
+
+plt.plot(x, Jy[10][2,:], 'b-', marker = 'o')
+plt.plot(x, Jz[10][2,:], 'r--',marker = '*')
+plt.axvline(x[2], ls='--', color='k')
+plt.axvline(x[-3], ls='--', color='k')
+
+plt.show()

diagnostics/whislerwave/plot_anim.py

@@ -4,13 +4,11 @@
 from matplotlib.animation import FuncAnimation
 import os
 
-nx = 512
-Dx = 0.1
+nx = 128
+Dx = 0.05
 
-Dt = 0.002
 
-
-quantity_name = 'Bz'
+quantity_name = 'vz'
 fixed_index = 0
 ny = 1
 
@@ -20,7 +18,8 @@ def read_file(filename):
     df = pd.read_csv(filename, delim_whitespace=True, header=None, names=column_names)
     return df
 
-results_dir = 'whislerwaveres'
+results_dir = 'whislerwaveres/'
+results_dir2 = 'whislerwaveHLL/'
 
 quantities = {
     'rho': [],
@@ -34,11 +33,21 @@ def read_file(filename):
 }
 times = []
 
+quantitieshll = {
+    'rho': [],
+    'vx': [],
+    'vy': [],
+    'vz': [],
+    'Bx': [],
+    'By': [],
+    'Bz': [],
+    'P': []
+}
+timeshll = []
+
 for filename in os.listdir(results_dir):
-    if filename.startswith("URK2_") and filename.endswith(".txt"):
-        # Extract time from filename and format it properly
-        time_str = filename.split('_')[2].split('.')[0]
-        time_str = time_str.replace('_', '.')  # Replace underscore with dot
+    if filename.startswith("PRK2_") and filename.endswith(".txt"):
+        time_str = filename.split('_')[1]+'.'+filename.split('_')[2].split('.')[0]
         time = float(time_str)
         times.append(time)
 
@@ -50,28 +59,48 @@ def read_file(filename):
 for quantity in quantities.keys():
     quantities[quantity] = np.array(quantities[quantity])
 
-times = np.array(times)
+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)
+
+        df = read_file(os.path.join(results_dir, filename))
+
+        for quantity in quantitieshll.keys():
+            quantitieshll[quantity].append(df[quantity].values.reshape((ny, nx)))
+
+for quantity in quantitieshll.keys():
+    quantitieshll[quantity] = np.array(quantitieshll[quantity])
 
 x=Dx*np.arange(nx) + 0.5*Dx
 
+def update(frame):
+    lx = nx*Dx
+    m = 20#int(nx/4)
 
-def update(frame):    
+    k = 2 * np.pi / lx * m
+
+    expected_value = 1e-2 * np.cos(k * x + k**2 * times[frame] + 0.5488135)
 
     plt.clf()
-    plt.plot(x, quantities[quantity_name][frame, fixed_index, :], color='blue', markersize=3) # t,y,x
-    plt.title(f'{quantity_name} at y={fixed_index}, t={frame * Dt}')  # Format time to one decimal place
+    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, expected_value)
+    plt.title(f'{quantity_name} at t={times[frame]}')  # Format time to one decimal place
     plt.xlabel('x')
     plt.ylabel(quantity_name)
     plt.grid(True)
     #plt.yscale("log")
     plt.tight_layout()
 
-    eps = 0.01
+    eps = 0.001
     min_val = np.min(quantities[quantity_name][:, fixed_index, :]) - eps
     max_val = np.max(quantities[quantity_name][:, fixed_index, :]) + eps
 
     plt.ylim(min_val, max_val)
-
+    plt.axvline(x[1]-Dx/2, ls='--', color='k')
+    plt.axvline(x[-1]-Dx/2, ls='--', color='k')
 
 fig = plt.figure(figsize=(8, 6))
 ani = FuncAnimation(fig, update, frames=len(times), interval=100)
@@ -85,4 +114,16 @@ def onclick(event):
 # Connect the key press event to the function
 fig.canvas.mpl_connect('key_press_event', onclick)
 
-plt.show()
+plt.show()
+
+"""
+lx = nx*Dx
+m = 4#int(nx/4)
+
+k = 2 * np.pi / lx * m
+plt.plot(x, quantities[quantity_name][0, fixed_index, :], color='blue', marker = 'x', markersize=3) # t,y,x
+plt.plot(x, 1e-3 * np.cos(k * x + k**2 * times[0] + 0.5488135))
+plt.axvline(x[1]-Dx/2, ls='--', color='k')
+plt.axvline(x[-1]-Dx/2, ls='--', color='k')
+plt.show()
+"""

diagnostics/whislerwave/space_convergence/space_convergence_plot.py

@@ -0,0 +1,121 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import os
+
+results_dir = './space_results'
+
+nx_initial = 32
+ny = 1
+
+fixed_index = 0
+time_index = 1
+
+studied_quantity = 'By'
+
+quantities = {
+    'rho': {},
+    'rhovx': {},
+    'rhovy': {},
+    'rhovz': {},
+    'Bx': {},
+    'By': {},
+    'Bz': {},
+    'Etot': {},
+}
+n_errors = 5
+times = []
+errors = np.zeros(n_errors)
+nx_values = {}
+
+# List and sort the files in the results directory
+filenames = os.listdir(results_dir)
+filenames.sort(key=lambda x: (int(x.split('_')[0]), float(x.split('_')[2].split('.')[0].replace('_', '.'))))
+
+for filename in filenames:
+    stepDx_str = filename.split('_')[0]  # Extract the part before "URK2_"
+    stepDx = int(stepDx_str)
+
+    # Compute nx for this stepDx
+    nx = nx_initial * (2 ** stepDx)
+    nx_values[stepDx] = nx
+
+    # Create a new dictionary entry for this Dx if it doesn't exist
+    if stepDx not in quantities['rho']:
+        for key in quantities:
+            quantities[key][stepDx] = []
+
+    # Extract time from filename and format it properly
+    time_str = filename.split('_')[2]+'.'+filename.split('_')[3].split('.')[0] # Extract the part after "URK2_" (length of "URK2_" is 5)
+    time_str = time_str.replace('_', '.')  # Replace underscore with dot
+    time = float(time_str)
+    times.append(time)
+
+    df = pd.read_csv(os.path.join(results_dir, filename), delim_whitespace=True, header=None, names=['rho', 'rhovx', 'rhovy', 'rhovz', 'Bx', 'By', 'Bz', 'Etot'])
+    for quantity in quantities:
+        quantities[quantity][stepDx].append(df[quantity].values.reshape(nx))   # Store the entire data array for each quantity
+
+times = np.array(times)
+times = np.unique(times)
+
+Dx =np.asarray( [0.2 / (2 ** i) for i in range(len(nx_values))])
+
+# Function to calculate L2 norm of error
+def calculate_error(computed, expected, nx):
+    #return np.linalg.norm(computed -expected)
+    return np.sum(abs(computed - expected))/nx
+
+for stepDx in quantities[studied_quantity]:
+    # Get nx for this stepDx
+    nx = nx_values[stepDx]
+
+    x = Dx[stepDx] * np.arange(nx) + 0.5*Dx[stepDx]
+
+    lx = nx*Dx[stepDx]
+    m = int(nx/4)
+
+    k = 2 * np.pi / lx * m
+
+    expected_value = 1e-2 * np.cos(k * x - k**2 * times[time_index] + 0.5488135)
+
+    # Extract computed value for this time_index
+    computed_value = quantities[studied_quantity][stepDx][time_index]
+
+    #print(computed_value, expected_value, "#####################")
+
+    # Calculate error for this Dx
+    error = calculate_error(computed_value, expected_value , nx) #/ error0
+
+    errors[stepDx] = error
+
+# Assuming Dx and errors are your data arrays
+# Log-transform the data
+log_Dx = np.log(Dx)
+log_errors = np.log(errors)
+
+# Perform a linear fit to the log-log data
+coefficients = np.polyfit(log_Dx, log_errors, 1)
+
+# Extract the slope (order of accuracy) and intercept
+slope, intercept = coefficients
+
+# Generate the fitted line for plotting
+fitted_log_errors = np.polyval(coefficients, log_Dx)
+fitted_errors = np.exp(fitted_log_errors)
+
+# Plot the original data and the fitted line
+plt.figure(figsize=(12, 8))
+plt.loglog(Dx, errors, marker='o', label='Numerical Error')
+#plt.loglog(Dx,np.exp(log_Dx*slope + intercept), label="manual")
+plt.loglog(Dx, fitted_errors, linestyle='--', label=f'Fitted Line (slope={slope:.2f})')
+plt.xlabel('Dx')
+plt.ylabel('L2 Norm of Error')
+plt.title('L2 Norm of Error vs. Dx')
+plt.grid(True)
+plt.legend()
+plt.show()
+
+for stepDx in quantities[studied_quantity]:
+    computed_value = quantities[studied_quantity][stepDx][time_index]
+    plt.plot(Dx[stepDx]*np.arange(nx_values[stepDx]), computed_value, label=f"dx = {Dx[stepDx]}")
+plt.legend()