post processing helpers

.gitignore

@@ -7,6 +7,7 @@ __pycache__
 *.gif
 *.mp4
 *.png
+*.pdf
 space_results/
 time_results/
 2Dwave/

diagnostics/Hallharris/colormesh.py

@@ -12,7 +12,7 @@ def reshape_data(data, ny, nx):
 nx = 250
 ny = 250
 
-file_path = "hallharrisresCT/PRK2_0_5002794389.txt"
+file_path = "hallharrisres2/PRK2_15_1429222576.txt"
 
 data = read_data(file_path)
 
@@ -34,12 +34,14 @@ def reshape_data(data, ny, nx):
 # Calculate Jz
 Jz = (dBy_dx - dBx_dy)
 
-toPlot = Bz
+toPlot = rho
 
 plt.pcolormesh(toPlot.T, cmap='coolwarm')  
-plt.title('Contour Plot of Bz at t=0.5')
+plt.colorbar()
+plt.title('Contour Plot of rho')
 plt.xlabel('X')
 plt.ylabel('Y')
+plt.savefig('hall_harris_rho.png')
 plt.show()
 
 # x_middle = nx // 2

diagnostics/Hallharris/colormesh_anim.py

@@ -24,7 +24,7 @@ def read_times(file_paths):
         times.append(time)
     return times
 
-results_dir = "hallharrisresCT2long/"
+results_dir = "hallharrisres2/"
 file_paths = [results_dir + file for file in os.listdir(results_dir) if file.startswith("PRK2_") and file.endswith(".txt")]
 
 nx = 250
@@ -53,7 +53,7 @@ def read_times(file_paths):
 # Calculate Jz
 Jz = [(dBy_dx[i] - dBx_dy[i]) for i in range(len(data))]
 
-toPlot = Bz
+toPlot = rho
 
 data_min = np.min(toPlot[-1])
 data_max = np.max(toPlot[-1])

diagnostics/MHDrotor/MHDrotor.py

@@ -24,6 +24,8 @@
 constainedtransport = p.CTMethod.Arithmetic
 timeintegrator = p.Integrator.TVDRK2Integrator
 
+dumpvariables = p.dumpVariables.Primitive
+
 consts = p.Consts(sigmaCFL = 0.2, gam = 1.4)
 
 ##############################################################################################################################################################################
@@ -73,7 +75,7 @@ def Bz_(x, y):
     return 0.0 
 
 def P_(x, y):
-    return 0.5
+    return 1
 
 
 x = np.arange(nx) * Dx + 0.5 * Dx
@@ -108,4 +110,4 @@ def P_(x, y):
 
 p.PhareMHD(P0cc, result_dir, nghost, 
            boundaryconditions, reconstruction, slopelimiter, riemannsolver, constainedtransport, timeintegrator,
-           Dx, Dy, FinalTime, Dt, Consts = consts)
+           Dx, Dy, FinalTime, Dt, Consts = consts, dumpvariables = dumpvariables)

diagnostics/MHDrotor/colormesh.py

@@ -12,13 +12,13 @@ def reshape_data(data, ny, nx):
 nx = 100
 ny = 100
 
-file_path = "MHDrotorres/URK2_0_0.txt"
+file_path = "MHDrotorres/PRK2_0_0.txt"
 
 data = read_data(file_path)
 
 reshaped_data = reshape_data(data, nx, ny)
 
-rho = reshaped_data[:, :, 0]
+rho = reshaped_data[:, :, 1]
 
 plt.pcolormesh(rho.T, cmap='coolwarm')  
 plt.title('Contour Plot of rho at t=0')

diagnostics/MHDrotor/colormesh_anim.py

@@ -2,6 +2,10 @@
 import matplotlib.pyplot as plt
 import os
 from matplotlib.animation import FuncAnimation, PillowWriter
+from matplotlib.colors import Normalize
+
+
+studied_index = 1
 
 # Function to read data from file
 def read_data(file_path):
@@ -23,7 +27,7 @@ def read_times(file_paths):
     return times
 
 results_dir = "MHDrotorres/"
-file_paths = [results_dir + file for file in os.listdir(results_dir) if file.startswith("URK2_") and file.endswith(".txt")]
+file_paths = [results_dir + file for file in os.listdir(results_dir) if file.startswith("PRK2_") and file.endswith(".txt")]
 
 nx = 100
 ny = 100
@@ -33,13 +37,18 @@ def read_times(file_paths):
 
 reshaped_data = [reshape_data(d, nx, ny) for d in data]
 
+data_min = np.min(reshaped_data[-1][:, :, studied_index])
+data_max = np.max(reshaped_data[-1][:, :, studied_index])
+
+Norm = Normalize(vmin=data_min, vmax=data_max)
+
 fig, ax = plt.subplots()
-im = ax.pcolormesh(reshaped_data[0][:, :, 0].T, cmap='coolwarm')
+im = ax.pcolormesh(reshaped_data[0][:, :, studied_index].T, cmap='coolwarm')
 ax.set_aspect('equal')
 fig.colorbar(im, ax=ax)
 
 def update(frame):
-    im.set_array(reshaped_data[frame][:, :, 0].T)
+    im.set_array(reshaped_data[frame][:, :, studied_index].T)
     plt.title(f'qty at t={times[frame]}')
     return im,
 

diagnostics/alfvenwave/X/space_convergence/space_convergence_plot.py

@@ -1,120 +1,107 @@
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
+from matplotlib.animation import FuncAnimation
 import os
+import shutil
 
-results_dir = './space_results'
-
-nx_initial = 50
+# Parameters
+nx = 50
 ny = 3
-
+Dt = 5e-4
+quantity_name = 'By'
 fixed_index = 0
-time_index = -1
+lx = 1
 
-studied_quantity = 'By'
+# Define column names
+column_names = ['rho', 'rhovx', 'rhovy', 'rhovz', 'Bx', 'By', 'Bz', 'Etot']
 
+def read_file(filename):
+    df = pd.read_csv(filename, delim_whitespace=True, header=None, names=column_names)
+    return df
+
+# Results directory
+results_dir = './space_results'
+
+# Dictionary to store quantities
 quantities = {
-    'rho': {},
-    'rhovx': {},
-    'rhovy': {},
-    'rhovz': {},
-    'Bx': {},
-    'By': {},
-    'Bz': {},
-    'Etot': {},
+    '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('_', '.'))))
+# Read all files in the results directory
+for filename in os.listdir(results_dir):
+    if filename.startswith("0_URK2_") and filename.endswith(".txt"):
+        # Extract time from filename and format it properly
+        time_str = filename.split('_')[2]+'.'+filename.split('_')[3].split('.')[0]
+        time = float(time_str)
+        times.append(time)
+
+        df = read_file(os.path.join(results_dir, filename))
 
-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)
+        for quantity in quantities.keys():
+            quantities[quantity].append(df[quantity].values.reshape((ny, nx)))
 
-    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((ny, nx)))   # Store the entire data array for each quantity
+# Convert lists to numpy arrays
+for quantity in quantities.keys():
+    quantities[quantity] = np.array(quantities[quantity])
 
 times = np.array(times)
-times = np.unique(times)
+times = times
+
+# Define x-axis
+Dx = lx / nx
+x = Dx * np.arange(nx) + 0.5 * Dx
+
+# Create output directory for frames
+output_dir = 'frames'
+if os.path.exists(output_dir):
+    shutil.rmtree(output_dir)
+os.makedirs(output_dir, exist_ok=True)
+
+# Define the update function for the animation
+def update(frame):
+    plt.clf()  # Clear the previous plot
+    t = times[frame]
+
+    # Plot the desired quantity
+    plt.plot(x, quantities[quantity_name][frame, fixed_index, :], color='blue')
+    plt.title(f'{quantity_name} at y={fixed_index}, t={times[frame]:.4f}')
+    plt.xlabel('x')
+    plt.ylabel(quantity_name)
+    plt.grid(True)
+    plt.tight_layout()
+
+    # Set y-limits based on data range
+    min_val = np.min(quantities[quantity_name][:, fixed_index, :])
+    max_val = np.max(quantities[quantity_name][:, fixed_index, :])
+    plt.ylim(min_val, max_val)
 
-Dx =np.asarray( [0.02 / (2 ** i) for i in range(len(nx_values))])
+    # Save the current frame
+    plt.savefig(f'{output_dir}/frame_{frame:04d}.png')
 
-# 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
+# Create the figure
+fig = plt.figure(figsize=(8, 6))
 
-for stepDx in quantities[studied_quantity]:
-    factor = 2 ** stepDx
+# Create animation
+ani = FuncAnimation(fig, update, frames=len(times), interval=100)
 
-    # Get nx for this stepDx
-    nx = nx_values[stepDx]
-
-    x = Dx[stepDx] * np.arange(nx) + 0.5 * Dx[stepDx]
+# Manually step through the frames with arrow keys
+def onclick(event):
+    if event.key == 'right':
+        ani.frame_seq = ani.new_frame_seq()
+        fig.canvas.draw()
 
-    #expected0 = 1e-6 * np.cos(2 * np.pi * (x - times[0] * (Dt)))
-    #computed0 = quantities[studied_quantity][stepDx][0][fixed_index, :]
-    #error0 = calculate_error(computed0, expected0, nx)
-
-    expected_value = 1e-6 * np.cos(2 * np.pi * (x - times[time_index]))
-
-    # Extract computed value for this time_index
-    computed_value = quantities[studied_quantity][stepDx][time_index][fixed_index, :]
-
-    # 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()
+# Connect key press event
+fig.canvas.mpl_connect('key_press_event', onclick)
 
-for stepDx in quantities[studied_quantity]:
-    computed_value = quantities[studied_quantity][stepDx][time_index][fixed_index, :]
-    plt.plot(Dx[stepDx]*np.arange(nx_values[stepDx]), computed_value, label=f"dx = {Dx[stepDx]}")
-plt.legend()    
+# Show plot (optional, can comment out if just saving frames)
+plt.show()

diagnostics/alfvenwave2D/alfvenwave2D.py

@@ -14,7 +14,7 @@
 Dx = 0.01
 Dy = 0.01
 Dt = 0.0
-FinalTime = 0.5
+FinalTime = 0.7
 nghost = 2
 
 boundaryconditions = p.BoundaryConditions.Periodic

diagnostics/alfvenwave2D/colormesh.py

@@ -37,7 +37,9 @@ def read_times(file_paths):
 qty = reshaped_data[time_index][:, :, 5]
 
 plt.pcolormesh(qty.T, cmap='coolwarm')  
-plt.title('Contour Plot of qty')
+plt.colorbar()
+plt.title('Contour Plot of By at t=0')
 plt.xlabel('X')
 plt.ylabel('Y')
+plt.savefig('2Dwave_0.png')
 plt.show()