first working hall mhd commit

diagnostics/orszag-tang/orszagtang.py

@@ -14,19 +14,21 @@
 Dy = 1/ny
 Dt = 0.0
 FinalTime = 0.5
-nghost = 2
+nghost = 1
 
 boundaryconditions = p.BoundaryConditions.Periodic
 
-reconstruction = p.Reconstruction.Linear
+reconstruction = p.Reconstruction.Constant
 slopelimiter = p.Slope.VanLeer
 riemannsolver = p.RiemannSolver.Rusanov
 constainedtransport = p.CTMethod.Average
-timeintegrator = p.Integrator.TVDRK2Integrator
+timeintegrator = p.Integrator.EulerIntegrator
 
 dumpvariables = p.dumpVariables.Primitive
 dumpfrequency = 80
 
+OptionalPhysics = p.OptionalPhysics.HallResHyper
+
 ##############################################################################################################################################################################
 B0 = 1./(np.sqrt(4.*np.pi))
 
@@ -85,4 +87,4 @@ def P_(x, y):
 
 p.PhareMHD(P0cc, result_dir, nghost, 
            boundaryconditions, reconstruction, slopelimiter, riemannsolver, constainedtransport, timeintegrator,
-           Dx, Dy, FinalTime, Dt, dumpvariables = dumpvariables, dumpfrequency = dumpfrequency)
+           Dx, Dy, FinalTime, Dt, dumpvariables = dumpvariables, dumpfrequency = dumpfrequency, OptionalPhysics = OptionalPhysics)

diagnostics/orszag-tang/plotJ.py

@@ -0,0 +1,95 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from matplotlib.animation import FuncAnimation
+import os
+
+nx = 128
+ny = 128
+Dx = 0.01
+Dy = 0.01
+
+nghost = 1
+nghostJ = 2
+
+Dt = 0.000625
+
+Jx = []
+Jy = []
+Jz = []
+times = []
+results_dir = 'orszagtangCTAverage/'
+
+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("Jx_") 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))
+        Jx.append(df.values.reshape((ny + 2*nghostJ + 1, nx + 2*nghostJ)))
+
+    if filename.startswith("Jy_") and filename.endswith(".txt"):
+
+        df = read_file(os.path.join(results_dir, filename))
+        Jy.append(df.values.reshape((ny + 2*nghostJ, nx + 2*nghostJ + 1)))
+
+    if filename.startswith("Jz_") and filename.endswith(".txt"):
+
+        df = read_file(os.path.join(results_dir, filename))
+        Jz.append(df.values.reshape((ny + 2*nghostJ + 1, nx + 2*nghostJ + 1)))
+
+x=Dx*np.arange(nx + 2*nghostJ + 1)-nghostJ*Dx
+y=Dy*np.arange(ny + 2*nghostJ)-nghostJ*Dy
+
+def update(frame):
+
+    m = 10
+    lx = nx*Dx
+    k = 2*np.pi/lx * m
+    #expectedJy = -k * np.cos(k*x + k**2 * times[frame]*Dt + 0.5488135)*0.01
+    #expectedJz = -k * np.sin(k*x + k**2 * times[frame]*Dt + 0.5488135)*0.01
+
+    plt.clf()
+
+    #plt.plot(x, expectedJy, 'y-', marker = '+')
+    #plt.plot(x, expectedJz, 'g-', marker = 'x')
+
+    plt.plot(x, Jy[frame][2,:], 'b-', marker = 'o')
+    #plt.plot(x, Jz[frame][2,:], 'r--',marker = '*')
+
+    #plt.plot(x, Jy2[frame][2,:], 'y-', marker = 'o')
+    #plt.plot(x, Jz2[frame][2,:], 'g--',marker = '*')
+
+    plt.xlabel('x')
+    plt.grid(True)
+    #plt.yscale("log")
+    plt.tight_layout()
+    plt.axvline(y[2], ls='--', color='k')
+    plt.axvline(y[-3], ls='--', color='k')
+
+
+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/fft.py → diagnostics/whistlerwave/fft.py

@@ -5,10 +5,10 @@
 import pandas as pd
 from matplotlib.colors import LogNorm, Normalize, PowerNorm
 
-nx = 512
+nx = 128
 Dx = 0.1
 
-Dt = 0.002
+Dt = 0.002411
 finalTime = 10
 
 column_names = ['rho', 'vx', 'vy', 'vz', 'Bx', 'By', 'Bz', 'P']
@@ -31,7 +31,7 @@ def read_file(filename):
 }
 
 for filename in os.listdir(results_dir):
-    if filename.startswith("URK2_") and filename.endswith(".txt"):
+    if filename.startswith("PRK2_") and filename.endswith(".txt"):
         df = read_file(os.path.join(results_dir, filename))
 
         for quantity in quantities.keys():
@@ -41,7 +41,10 @@ def read_file(filename):
     quantities[quantity] = np.array(quantities[quantity])
 
 def whistler(k):
-    return k**2
+    return (k**2 /2) * (np.sqrt(1+4/k**2) + 1)
+
+def leftAlfvenIonCyclotron(k):
+    return (k**2 /2) * (np.sqrt(1+4/k**2) - 1)
 
 k = 2*np.pi*fftfreq(nx, Dx)
 w = 2*np.pi*fftfreq(int(finalTime/Dt) + 1, Dt)
@@ -51,8 +54,8 @@ def whistler(k):
 kplus = k[:halfk]
 wplus = w[:halfw]
 
-modes = (1,2,4)
-kmodes = 2*np.pi*np.asarray([1/(nx*Dx) * m for m in modes])
+"""modes = [4]
+kmodes = 2*np.pi*np.asarray([1/(nx*Dx) * m for m in modes])"""
 
 """
 fig,ax = plt.subplots()
@@ -72,12 +75,13 @@ def whistler(k):
 Xx, Yy = np.meshgrid(kplus, wplus)
 
 fig, ax = plt.subplots()
-pcm = ax.pcolormesh(Xx, Yy, B_fft_abs, cmap='plasma', shading='auto', norm=LogNorm(vmin=B_fft_abs.min(), vmax=B_fft_abs.max()))
+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:
+"""for km in kmodes:
     ax.axvline(km, ls='--', color='k')
-    ax.plot(km, km**2, marker='o')
+    ax.plot(km, km**2, marker='o')"""
 
 plt.show()

diagnostics/whislerwave/plotJ.py → diagnostics/whistlerwave/plotJ.py

@@ -9,8 +9,9 @@
 
 Dt = 0.000625
 
-Jz = []
+
 Jy = []
+Jz = []
 times = []
 results_dir = 'whislerwaveres/'
 
@@ -25,14 +26,14 @@ def read_file(filename):
         times.append(time)
 
         df = read_file(os.path.join(results_dir, filename))
-        Jy.append(df.values.reshape((5,133)))
+        Jy.append(df.values.reshape((7,135)))
 
     if filename.startswith("Jz_") and filename.endswith(".txt"):
 
         df = read_file(os.path.join(results_dir, filename))
-        Jz.append(df.values.reshape((6,133)))
+        Jz.append(df.values.reshape((8,135)))
 
-x=Dx*np.arange(133)-2*Dx
+x=Dx*np.arange(135)-2*Dx
 
 def update(frame):
 
@@ -44,21 +45,27 @@ def update(frame):
 
     plt.clf()
 
-    plt.plot(x, expectedJy, 'y-', marker = '+')
+    #plt.plot(x, expectedJy, 'y-', marker = '+')
     #plt.plot(x, expectedJz, 'g-', marker = 'x')
 
-    plt.plot(x, Jy[frame][2,:], 'b-', marker = 'o')
+    plt.plot(x, Jy[frame][3,:], 'b-', marker = 'o')
     #plt.plot(x, Jz[frame][2,:], 'r--',marker = '*')
+
+    #plt.plot(x, Jy2[frame][2,:], 'y-', marker = 'o')
+    #plt.plot(x, Jz2[frame][2,:], 'g--',marker = '*')
+
     plt.xlabel('x')
     plt.grid(True)
     #plt.yscale("log")
     plt.tight_layout()
 
-    eps = 0.01
+    eps = 1e-7
     min_val = np.min(Jy) - eps
     max_val = np.max(Jz) + eps
 
     plt.ylim(min_val, max_val)
+    plt.axvline(x[2], ls='--', color='k')
+    plt.axvline(x[-3], ls='--', color='k')
 
 
 fig = plt.figure(figsize=(8, 6))
@@ -75,9 +82,9 @@ def onclick(event):
 
 plt.show()
 
-plt.plot(x, Jy[10][2,:], 'b-', marker = 'o')
+"""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()
+plt.show()"""

diagnostics/whislerwave/plot_anim.py → diagnostics/whistlerwave/plot_anim.py

@@ -3,12 +3,13 @@
 import matplotlib.pyplot as plt
 from matplotlib.animation import FuncAnimation
 import os
+import shutil
 
 nx = 128
-Dx = 0.05
+Dx = 0.1
 
 
-quantity_name = 'Bz'
+quantity_name = 'By'
 fixed_index = 0
 ny = 1
 
@@ -75,13 +76,20 @@ def read_file(filename):
 
 x=Dx*np.arange(nx) + 0.5*Dx
 
+output_dir = 'frames'
+if os.path.exists(output_dir):
+    shutil.rmtree(output_dir)
+os.makedirs(output_dir, exist_ok=True)
+
 def update(frame):
     lx = nx*Dx
-    m = 10#int(nx/4)
+    m = 4#int(nx/4)
 
     k = 2 * np.pi / lx * m
+
+    w = (k**2 /2) *(np.sqrt(1+4/k**2) + 1)
 
-    expected_value = 1e-2 * np.sin(k * x + k**2 * times[frame] + 0.5488135)
+    expected_value = 1e-7 * np.cos(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
@@ -94,11 +102,12 @@ def update(frame):
     #plt.yscale("log")
     plt.tight_layout()
 
-    eps = 0.001
+    eps = 1e-7
     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.savefig(f'{output_dir}/frame_{frame:04d}.png')
     #plt.axvline(x[2]-Dx/2, ls='--', color='k')
     #plt.axvline(x[-2]-Dx/2, ls='--', color='k')
 

diagnostics/whislerwave/whislerwave.py → diagnostics/whistlerwave/whislerwave.py

@@ -10,15 +10,15 @@
 
 nx = 128
 ny = 1
-Dx = 0.05
+Dx = 0.1
 Dy = 1
-Dt = 0.000625
+Dt = 0.00
 FinalTime = 0.5
-nghost = 1
+nghost = 2
 
 boundaryconditions = p.BoundaryConditions.Periodic
 
-reconstruction = p.Reconstruction.Constant
+reconstruction = p.Reconstruction.Linear
 
 slopelimiter = p.Slope.VanLeer
 riemannsolver = p.RiemannSolver.Rusanov
@@ -37,7 +37,7 @@
 
 np.random.seed(0)
 
-modes = [10]#[int(nx/4)]
+modes = [4]#[int(nx/4)]
 phases = np.random.rand(len(modes))
 
 def rho_(x, y):
@@ -48,14 +48,16 @@ def vx_(x, y):
 
 def vy_(x, y):
     ret = np.zeros((x.shape[0], y.shape[1]))
+    w = (k**2 /2) *(np.sqrt(1+4/k**2) + 1)
     for m,phi in zip(modes, phases):
-        ret[:,:] += np.cos(k*x*m + phi)*0.01
+        ret[:,:] += -np.cos(k*x*m + phi)*1e-7*k
     return ret
 
 def vz_(x, y):
     ret = np.zeros((x.shape[0], y.shape[1]))
+    w = (k**2 /2) *(np.sqrt(1+4/k**2) + 1)
     for m,phi in zip(modes, phases):
-        ret[:,:] += np.sin(k*x*m + phi)*0.01
+        ret[:,:] += np.sin(k*x*m + phi)*1e-7*k
     return ret
 
 def Bx_(x, y):
@@ -64,17 +66,17 @@ 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)*0.01
+        ret[:,:] += np.cos(k*x*m + phi)*1e-7
     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)*0.01
+        ret[:,:] += -np.sin(k*x*m + phi)*1e-7
     return ret
 
 def P_(x, y):
-    return 1e-3
+    return 1.0
 
 x = np.arange(nx) * Dx + 0.5 * Dx 
 y = np.arange(ny) * Dy + 0.5 * Dy 

src/AddGhostCells.hpp

@@ -139,8 +139,8 @@ void UpdateGhostJ(Variables& Vn, int nghost, BoundaryConditions bc) {
         // Jx
         for (int i = nghostJ; i < nx - nghostJ; i++) { // Ghost cells still not good for Jx, dir y and ZeroGradient (maybe)
             for (int k = 1; k <= nghostJ; k++) {
-                Vn.Jx[nghostJ - k][i] = Vn.Jx[ny1 - k - nghostJ][i]; // Bottom side
-                Vn.Jx[ny1 - 1 + k - nghostJ][i] = Vn.Jx[k - 1 + nghostJ][i]; // Top side
+                Vn.Jx[nghostJ - k][i] = Vn.Jx[ny1 - k - nghostJ-1][i]; // Bottom side
+                Vn.Jx[ny1 - 1 + k - nghostJ][i] = Vn.Jx[k - 1 + nghostJ+1][i]; // Top side
             }
         }