Corrected data passing in conversion constructor, added dumps and dat…

…a visualisation for hall mhd debugging

diagnostics/whislerwave/plotFlux.py

@@ -0,0 +1,94 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from matplotlib.animation import FuncAnimation
+import os
+
+nx = 128
+Dx = 0.05
+
+
+quantity_name = 'Bz'
+fixed_index = 0
+ny = 1
+
+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/'
+
+quantities = {
+    'rho': [],
+    'vx': [],
+    'vy': [],
+    'vz': [],
+    'Bx': [],
+    'By': [],
+    'Bz': [],
+    'P': []
+}
+times = []
+
+for filename in os.listdir(results_dir):
+    if filename.startswith("FluxDifx_") 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))
+
+        for quantity in quantities.keys():
+            quantities[quantity].append(df[quantity].values.reshape((ny, nx)))
+
+for quantity in quantities.keys():
+    quantities[quantity] = np.array(quantities[quantity])
+
+x=Dx*np.arange(nx) + 0.5*Dx
+
+def update(frame):
+
+    plt.clf()
+    plt.plot(x, quantities[quantity_name][frame, fixed_index, :], color='blue', marker = 'x', markersize=3) # t,y,x
+    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.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[2]-Dx/2, ls='--', color='k')
+    #plt.axvline(x[-2]-Dx/2, ls='--', color='k')
+
+fig = plt.figure(figsize=(8, 6))
+ani = FuncAnimation(fig, update, frames=len(times), interval=100)
+
+# 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()
+
+"""
+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/plotJ.py

@@ -4,8 +4,10 @@
 from matplotlib.animation import FuncAnimation
 import os
 
+nx = 128
+Dx = 0.05
 
-Dx = 0.2
+Dt = 0.000625
 
 Jz = []
 Jy = []
@@ -23,19 +25,30 @@ def read_file(filename):
         times.append(time)
 
         df = read_file(os.path.join(results_dir, filename))
-        Jy.append(df.values.reshape((5,37)))
+        Jy.append(df.values.reshape((5,133)))
 
     if filename.startswith("Jz_") and filename.endswith(".txt"):
 
         df = read_file(os.path.join(results_dir, filename))
-        Jz.append(df.values.reshape((6,37)))
+        Jz.append(df.values.reshape((6,133)))
 
-x=Dx*np.arange(37)-2*Dx
+x=Dx*np.arange(133)-2*Dx
 
 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, Jz[frame][2,:], 'r--',marker = '*')
     plt.xlabel('x')
     plt.grid(True)
     #plt.yscale("log")

diagnostics/whislerwave/plot_anim.py

@@ -8,7 +8,7 @@
 Dx = 0.05
 
 
-quantity_name = 'vz'
+quantity_name = 'Bz'
 fixed_index = 0
 ny = 1
 
@@ -33,7 +33,7 @@ def read_file(filename):
 }
 times = []
 
-quantitieshll = {
+"""quantitieshll = {
     'rho': [],
     'vx': [],
     'vy': [],
@@ -43,7 +43,7 @@ def read_file(filename):
     'Bz': [],
     'P': []
 }
-timeshll = []
+timeshll = []"""
 
 for filename in os.listdir(results_dir):
     if filename.startswith("PRK2_") and filename.endswith(".txt"):
@@ -59,33 +59,33 @@ def read_file(filename):
 for quantity in quantities.keys():
     quantities[quantity] = np.array(quantities[quantity])
 
-for filename in os.listdir(results_dir2):
+"""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))
+        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 quantitieshll.keys():
-    quantitieshll[quantity] = np.array(quantitieshll[quantity])
+    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)
+    m = 10#int(nx/4)
 
     k = 2 * np.pi / lx * m
 
-    expected_value = 1e-2 * np.cos(k * x + k**2 * times[frame] + 0.5488135)
+    expected_value = 1e-2 * np.sin(k * x + k**2 * 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, 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')
@@ -99,8 +99,8 @@ def update(frame):
     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')
+    #plt.axvline(x[2]-Dx/2, ls='--', color='k')
+    #plt.axvline(x[-2]-Dx/2, ls='--', color='k')
 
 fig = plt.figure(figsize=(8, 6))
 ani = FuncAnimation(fig, update, frames=len(times), interval=100)

diagnostics/whislerwave/whislerwave.py

@@ -14,7 +14,7 @@
 Dy = 1
 Dt = 0.000625
 FinalTime = 0.5
-nghost = 2
+nghost = 1
 
 boundaryconditions = p.BoundaryConditions.Periodic
 
@@ -23,7 +23,7 @@
 slopelimiter = p.Slope.VanLeer
 riemannsolver = p.RiemannSolver.Rusanov
 constainedtransport = p.CTMethod.Average
-timeintegrator = p.Integrator.TVDRK2Integrator
+timeintegrator = p.Integrator.EulerIntegrator
 
 consts = p.Consts(sigmaCFL = 0.8, gam = 5/3, eta = 0.0, nu = 0.000)
 physics = p.OptionalPhysics.HallResHyper
@@ -37,7 +37,7 @@
 
 np.random.seed(0)
 
-modes = [20]#[int(nx/4)]
+modes = [10]#[int(nx/4)]
 phases = np.random.rand(len(modes))
 
 def rho_(x, y):

src/AddGhostCells.hpp

@@ -137,7 +137,7 @@ void UpdateGhostJ(Variables& Vn, int nghost, BoundaryConditions bc) {
 
     if(bc == BoundaryConditions::Periodic){
         // Jx
-        for (int i = nghostJ; i < nx - nghostJ; i++) {
+        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

src/ComputeJ.hpp

@@ -25,7 +25,9 @@ void ComputeJ(Variables& Vn, double Dx, double Dy, int nghost) {
     {
         for (int i = nghost; i <= Vn.nx - nghost; ++i)
         {
-            Vn.Jz[j + 1][i + 1] = (Vn.Byf[j][i] - Vn.Byf[j][i - 1]) / (Dx) - (Vn.Bxf[j][i] - Vn.Bxf[j - 1][i]) / (Dy);
+            // In 1D (no CT), use cell centered Bx and By
+            Vn.Jz[j + 1][i + 1] = (Vn.By[j][i] - Vn.By[j][i - 1]) / (Dx);
+            //Vn.Jz[j + 1][i + 1] = (Vn.Byf[j][i] - Vn.Byf[j][i - 1]) / (Dx) - (Vn.Bxf[j][i] - Vn.Bxf[j - 1][i]) / (Dy);
         }
     }
 }

src/ConservativeVariables.cpp

@@ -71,14 +71,14 @@ ConservativeVariables::ConservativeVariables(const PrimitiveVariables& P_cc) : n
             Jy[j][i] = P_cc.Jy[j][i];
             Jz[j][i] = P_cc.Jz[j][i];
         }
-        Jy[j][nxJ] = P_cc.Jx[j][nxJ];
-        Jz[j][nxJ] = P_cc.Jx[j][nxJ];
+        Jy[j][nxJ] = P_cc.Jy[j][nxJ];
+        Jz[j][nxJ] = P_cc.Jz[j][nxJ];
     }
     for (int i = 0; i < nxJ; i++) {
         Jx[nyJ][i] = P_cc.Jx[nyJ][i];
-        Jz[nyJ][i] = P_cc.Jx[nyJ][i];
+        Jz[nyJ][i] = P_cc.Jz[nyJ][i];
     }
-    Jz[nyJ][nxJ] = P_cc.Jx[nyJ][nxJ];
+    Jz[nyJ][nxJ] = P_cc.Jz[nyJ][nxJ];
 }
 
 ConservativeVariables::~ConservativeVariables() = default;

src/GodunovFlux.cpp

@@ -1,5 +1,7 @@
 #include "GodunovFlux.hpp"
 
+#include "WrittingUtils.hpp"
+
 std::vector<ReconstructedValues> GodunovFluxX(const PrimitiveVariables& P_cc, double Dx, double Dy, int nghost, Reconstruction rec, Slope sl, Riemann rs, OptionalPhysics OptP){
     RiemannSolverFunction ChosenRiemannSolver = getRiemannSolver(rs);
     std::vector<ReconstructedValues> NumFluxx;
@@ -34,6 +36,15 @@ ConservativeVariables ComputeFluxDifferenceX(PrimitiveVariables& P_cc, double Dx
             FluxDifx.set(NumFluxx[(i+1)+(P_cc.nx + 1 - 2*nghost)*j] - NumFluxx[i+(P_cc.nx + 1 - 2*nghost)*j], i + nghost, j + nghost);
         }
     }
+
+
+    static int i = 0;
+    std::ostringstream fdifx;
+    fdifx << "whislerwaveres/FluxDifx_" << i <<".txt";
+    saveConcervativeVariables(FluxDifx, fdifx.str(), nghost);
+    i++;
+
+
     return FluxDifx;
 }
 

src/Interface.cpp

@@ -29,16 +29,22 @@ ReconstructedValues ComputeFluxVector(const ReconstructedValues& u, Dir dir) {
 }
 
 void AddNonIdealFlux(ReconstructedValues& f, const ReconstructedValues& u, double Jx, double Jy, double Jz, double LaplJx, double LaplJy, double LaplJz, OptionalPhysics OptP, Dir dir){
-    if (OptP == OptionalPhysics::HallResHyper) {
+    if (OptP == OptionalPhysics::HallRes) {
         if (dir == Dir::X) {
-            f.Bz += (1.0/u.rho) * (Jz * u.Bx - Jx * u.Bz);
+            f.By += - (1.0/u.rho) * (Jx * u.By - Jy * u.Bx) *1.0;
+            //f.By += - pc.eta * Jz;
+            //f.By -= - pc.nu * LaplJz;
+            f.Bz += (1.0/u.rho) * (Jz * u.Bx - Jx * u.Bz) *1.0;
             //f.Bz += pc.eta * Jy;
             //f.Bz -= pc.nu * LaplJy;
 
-            f.P += (1.0/u.rho) * ((u.Bx * Jx + u.By * Jy + u.Bz * Jz)*u.Bx - (u.Bx * u.Bx + u.By * u.By + u.Bz * u.Bz) * Jx);
+            f.P += (1.0/u.rho) * ((u.Bx * Jx + u.By * Jy + u.Bz * Jz)*u.Bx - (u.Bx * u.Bx + u.By * u.By + u.Bz * u.Bz) * Jx) *1.0;
             //f.P += pc.eta * (Jy * u.Bz - Jz * u.By);
             //f.P -= pc.nu * (LaplJy * u.Bz - LaplJz * u.By);
         } else if (dir == Dir::Y) {
+            f.Bx += (1.0/u.rho) * (Jx * u.By - Jy * u.Bx);
+            //f.Bx += pc.eta * Jz;
+            //f.Bx -= pc.nu * LaplJz;
             f.Bz += - (1.0/u.rho) * (Jy * u.Bz - Jz * u.By);
             //f.Bz += - pc.eta * Jx;
             //f.Bz -= - pc.nu * LaplJx;
@@ -250,6 +256,14 @@ Interface::Interface(const PrimitiveVariables& P_cc /* Assuming ghost cells are
         cwyL = std::sqrt(uL.Bx*uL.Bx + uL.By*uL.By + uL.Bz*uL.Bz) / (uL.rho) * vwy;
         cwxR = std::sqrt(uR.Bx*uR.Bx + uR.By*uR.By + uR.Bz*uR.Bz) / (uR.rho) * vwx;
         cwyR = std::sqrt(uR.Bx*uR.Bx + uR.By*uR.By + uR.Bz*uR.Bz) / (uR.rho) * vwy;
+        //double hallxL = std::sqrt(uL.Bx*uL.Bx + uL.By*uL.By + uL.Bz*uL.Bz) / (2.0*(uL.rho)*Dx);
+        //double hallyL = std::sqrt(uL.Bx*uL.Bx + uL.By*uL.By + uL.Bz*uL.Bz) / (2.0*(uL.rho)*Dy);
+        //double hallxR = std::sqrt(uR.Bx*uR.Bx + uR.By*uR.By + uR.Bz*uR.Bz) / (2.0*(uR.rho)*Dx);
+        //double hallyR = std::sqrt(uR.Bx*uR.Bx + uR.By*uR.By + uR.Bz*uR.Bz) / (2.0*(uR.rho)*Dy);
+        //cwxL = std::fabs(hallxL) + std::sqrt(hallxL*hallxL + caxL*caxL);
+        //cwyL = std::fabs(hallyL) + std::sqrt(hallyL*hallyL + cayL*cayL);
+        //cwxR = std::fabs(hallxR) + std::sqrt(hallxR*hallxR + caxR*caxR);
+        //cwyR = std::fabs(hallyR) + std::sqrt(hallyR*hallyR + cayR*cayR);
 
         if(dir == Dir::X){
             SLb = std::min(uL.vx - cfastxL - cwxL, uR.vx - cfastxR - cwxR);

src/PrimitiveVariables.cpp

@@ -73,14 +73,14 @@ PrimitiveVariables::PrimitiveVariables(const ConservativeVariables& C_cc) : nx(C
             Jy[j][i] = C_cc.Jy[j][i];
             Jz[j][i] = C_cc.Jz[j][i];
         }
-        Jy[j][nxJ] = C_cc.Jx[j][nxJ];
-        Jz[j][nxJ] = C_cc.Jx[j][nxJ];
+        Jy[j][nxJ] = C_cc.Jy[j][nxJ];
+        Jz[j][nxJ] = C_cc.Jz[j][nxJ];
     }
     for (int i = 0; i < nxJ; i++) {
         Jx[nyJ][i] = C_cc.Jx[nyJ][i];
-        Jz[nyJ][i] = C_cc.Jx[nyJ][i];
+        Jz[nyJ][i] = C_cc.Jz[nyJ][i];
     }
-    Jz[nyJ][nxJ] = C_cc.Jx[nyJ][nxJ];
+    Jz[nyJ][nxJ] = C_cc.Jz[nyJ][nxJ];
 
 }
 

src/RiemannSolver.cpp

@@ -2,7 +2,7 @@
 
 ReconstructedValues RusanovRiemannSolver(const Interface& inter){
     ReconstructedValues Frus;
-    Frus = 0.5 * (inter.fL + inter.fR) - 0.5 * inter.Splusb * (inter.uR - inter.uL);
+    Frus = 0.5 * (inter.fL + inter.fR) - 0.5 * inter.Splus * (inter.uR - inter.uL);
 
     if (inter.OP == OptionalPhysics::HallResHyper){
         Frus.Bx = 0.5 * (inter.fL.Bx + inter.fR.Bx) - 0.5 * inter.Splusb * (inter.uR.Bx - inter.uL.Bx);
@@ -11,10 +11,6 @@ ReconstructedValues RusanovRiemannSolver(const Interface& inter){
         Frus.P = 0.5 * (inter.fL.P + inter.fR.P) - 0.5 * inter.Splusb * (inter.uR.P - inter.uL.P);
     }
 
-    std::cout<<(inter.fL.vy + inter.fR.vy)<<" "<<inter.Splus<<
-    " "<<inter.Splusb<<" "<<(inter.uR.vy - inter.uL.vy)<<" "<<inter.SL<<
-    " "<<inter.SR<<" "<<inter.SLb<<" "<<inter.SRb<<std::endl;
-
     return Frus;
 }