From 0b582e710af1e63bd07623ff1c73f6405240a486 Mon Sep 17 00:00:00 2001 From: UCaromel Date: Wed, 24 Jul 2024 19:03:29 +0200 Subject: [PATCH] progress on hall debugging --- .gitignore | 1 + diagnostics/whislerwave/plotJ.py | 70 ++++++++++ diagnostics/whislerwave/plot_anim.py | 75 ++++++++--- .../space_convergence_plot.py | 121 +++++++++++++++++ .../space_convergence/whislerwave.py | 122 ++++++++++++++++++ diagnostics/whislerwave/whislerwave.py | 25 ++-- src/AddGhostCells.hpp | 8 +- src/ConstainedTransport.cpp | 8 +- src/Interface.cpp | 45 +++++-- src/Interface.hpp | 2 + src/PhareMHD.cpp | 5 - src/RiemannSolver.cpp | 51 +++++++- src/TimeIntegrator.cpp | 11 ++ 13 files changed, 483 insertions(+), 61 deletions(-) create mode 100644 diagnostics/whislerwave/plotJ.py create mode 100644 diagnostics/whislerwave/space_convergence/space_convergence_plot.py create mode 100644 diagnostics/whislerwave/space_convergence/whislerwave.py diff --git a/.gitignore b/.gitignore index a55ef7d..548a9ce 100644 --- a/.gitignore +++ b/.gitignore @@ -22,6 +22,7 @@ harrisres/ currentres/ whislerwaveres/ whislerwavelin/ +whislerwaveHLL/ frames/ pyMHD/pyMHD.cpython* \ No newline at end of file diff --git a/diagnostics/whislerwave/plotJ.py b/diagnostics/whislerwave/plotJ.py new file mode 100644 index 0000000..1c830d9 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/diagnostics/whislerwave/plot_anim.py b/diagnostics/whislerwave/plot_anim.py index d37ed02..b3a7080 100644 --- a/diagnostics/whislerwave/plot_anim.py +++ b/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() \ No newline at end of file +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() +""" diff --git a/diagnostics/whislerwave/space_convergence/space_convergence_plot.py b/diagnostics/whislerwave/space_convergence/space_convergence_plot.py new file mode 100644 index 0000000..7861840 --- /dev/null +++ b/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() diff --git a/diagnostics/whislerwave/space_convergence/whislerwave.py b/diagnostics/whislerwave/space_convergence/whislerwave.py new file mode 100644 index 0000000..5bf4f92 --- /dev/null +++ b/diagnostics/whislerwave/space_convergence/whislerwave.py @@ -0,0 +1,122 @@ +import sys +sys.path.append('/home/caromel/Documents/PHAREMHD/pyMHD') + +import numpy as np +import pyMHD as p +import os +import shutil + +############################################################################################################################################################################# + +initial_nx = 32 +initial_Dx = 0.2 +ny = 1 +Dy = 1 +Dt = 1e-5 +FinalTime = 1 +nghost = 2 + +boundaryconditions = p.BoundaryConditions.Periodic + +reconstruction = p.Reconstruction.Linear +slopelimiter = p.Slope.VanLeer +riemannsolver = p.RiemannSolver.HLL +constainedtransport = p.CTMethod.Average +timeintegrator = p.Integrator.TVDRK2Integrator + +consts = p.Consts(sigmaCFL = 0.8, gam = 5/3, eta = 0.0, nu = 0.00) +physics = p.OptionalPhysics.HallResHyper + +dumpvariables = p.dumpVariables.Conservative +dumpfrequency = 1000 + +############################################################################################################################################################################## +def initialize_variables(nx, ny, Dx, Dy): + lx=nx*Dx + + np.random.seed(0) + + modes = [int(nx/4)] + phases = np.random.rand(len(modes)) + + def rho_(x, y): + return 1.0 + + def vx_(x, y): + return 0.0 + + def vy_(x, y): + ret = np.zeros((x.shape[0], y.shape[1])) + for m,phi in zip(modes, phases): + ret[:,:] += np.cos(2*np.pi*x/lx*m + phi)*0.01 + 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(2*np.pi*x/lx*m + phi)*0.01 + return ret + + def Bx_(x, y): + return 1.0 + + def By_(x, y): + ret = np.zeros((x.shape[0], y.shape[1])) + for m,phi in zip(modes, phases): + ret[:,:] += np.cos(2*np.pi*x/lx*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(2*np.pi*x/lx*m + phi)*0.01 + return ret + + def P_(x, y): + return 1e-3 + + + x = np.arange(nx) * Dx + 0.5 * Dx + y = np.arange(ny) * Dy + 0.5 * Dy + xf = np.arange(nx+1) * Dx + yf = np.arange(ny+1) * Dy + + xx, yy = np.meshgrid(x, y, indexing = 'ij') + xfx, yfx = np.meshgrid(xf, y, indexing = 'ij') + xfy, yfy = np.meshgrid(x, yf, indexing = 'ij') + + rho = np.full((nx, ny), rho_(xx, yy)).T + vx = np.full((nx, ny), vx_(xx, yy)).T + vy = np.full((nx, ny), vy_(xx, yy)).T + vz = np.full((nx, ny), vz_(xx, yy)).T + Bxf = np.full((nx + 1, ny), Bx_(xfx, yfx)).T + Byf = np.full((nx, ny + 1), By_(xfy, yfy)).T + Bz = np.full((nx, ny), Bz_(xx, yy)).T + P = np.full((nx, ny), P_(xx, yy)).T + return rho, vx, vy, vz, Bxf, Byf, Bz, P + +############################################################################################################################################################################ +results_dir = 'space_results/' + +if os.path.exists(results_dir): + shutil.rmtree(results_dir) +os.makedirs(results_dir, exist_ok=True) + +Dx = initial_Dx +nx = initial_nx +stepDx = 0 + +while Dx > initial_Dx / 32.0 and nx <= 512: + rho, vx, vy, vz, Bxf, Byf, Bz, P = initialize_variables(nx, ny, Dx, Dy) + P0cc = p.PrimitiveVariables(nx, ny) + P0cc.init(rho, vx, vy, vz, Bxf, Byf, Bz, P) + + current_results_dir = os.path.join(results_dir, f'{stepDx}_') + + p.PhareMHD(P0cc, current_results_dir, nghost, + boundaryconditions, reconstruction, slopelimiter, riemannsolver, constainedtransport, + timeintegrator, Dx, Dy, FinalTime, Dt, dumpfrequency=dumpfrequency, Consts = consts, OptionalPhysics = physics) + + stepDx += 1 + Dx /= 2.0 + nx *= 2 \ No newline at end of file diff --git a/diagnostics/whislerwave/whislerwave.py b/diagnostics/whislerwave/whislerwave.py index 5cf00ab..bad59cb 100644 --- a/diagnostics/whislerwave/whislerwave.py +++ b/diagnostics/whislerwave/whislerwave.py @@ -8,26 +8,27 @@ ############################################################################################################################################################################# -nx = 512 +nx = 128 ny = 1 -Dx = 0.1 +Dx = 0.05 Dy = 1 -Dt = 0.002 -FinalTime = 10 +Dt = 0.000625 +FinalTime = 0.5 nghost = 2 boundaryconditions = p.BoundaryConditions.Periodic -reconstruction = p.Reconstruction.Linear +reconstruction = p.Reconstruction.Constant + slopelimiter = p.Slope.VanLeer -riemannsolver = p.RiemannSolver.HLL +riemannsolver = p.RiemannSolver.Rusanov constainedtransport = p.CTMethod.Average timeintegrator = p.Integrator.TVDRK2Integrator consts = p.Consts(sigmaCFL = 0.8, gam = 5/3, eta = 0.0, nu = 0.000) physics = p.OptionalPhysics.HallResHyper -dumpvariables = p.dumpVariables.Conservative +dumpvariables = p.dumpVariables.Primitive dumpfrequency = 1 ############################################################################################################################################################################## @@ -36,7 +37,7 @@ np.random.seed(0) -modes = (64,128) +modes = [20]#[int(nx/4)] phases = np.random.rand(len(modes)) def rho_(x, y): @@ -48,13 +49,13 @@ def vx_(x, y): def vy_(x, y): ret = np.zeros((x.shape[0], y.shape[1])) for m,phi in zip(modes, phases): - ret[:,:] += np.cos(2*np.pi*x/lx*m + phi)*0.01 + ret[:,:] += np.cos(k*x*m + phi)*0.01 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(2*np.pi*x/lx*m + phi)*0.01 + ret[:,:] += np.sin(k*x*m + phi)*0.01 return ret def Bx_(x, y): @@ -63,13 +64,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(2*np.pi*x/lx*m + phi)*0.01 + 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(2*np.pi*x/lx*m + phi)*0.01 + ret[:,:] += np.sin(k*x*m + phi)*0.01 return ret def P_(x, y): diff --git a/src/AddGhostCells.hpp b/src/AddGhostCells.hpp index 01df281..568513c 100644 --- a/src/AddGhostCells.hpp +++ b/src/AddGhostCells.hpp @@ -161,8 +161,8 @@ void UpdateGhostJ(Variables& Vn, int nghost, BoundaryConditions bc) { for (int j = nghostJ; j < ny - nghostJ; j++) { for (int k = 1; k <= nghostJ; k++) { - Vn.Jy[j][nghostJ - k] = Vn.Jy[j][nx1 - k - nghostJ]; // Left side - Vn.Jy[j][nx1 - 1 + k - nghostJ] = Vn.Jy[j][k - 1 + nghostJ]; // Right side + Vn.Jy[j][nghostJ - k] = Vn.Jy[j][nx1 - k - nghostJ-1]; // Left side + Vn.Jy[j][nx1 - 1 + k - nghostJ] = Vn.Jy[j][k - 1 + nghostJ+1]; // Right side } } @@ -176,8 +176,8 @@ void UpdateGhostJ(Variables& Vn, int nghost, BoundaryConditions bc) { for (int j = nghostJ; j < ny1 - nghostJ; j++) { for (int k = 1; k <= nghostJ; k++) { - Vn.Jz[j][nghostJ - k] = Vn.Jz[j][nx1 - k - nghostJ]; // Left side - Vn.Jz[j][nx1 - 1 + k - nghostJ] = Vn.Jz[j][k - 1 + nghostJ]; // Right side + Vn.Jz[j][nghostJ - k] = Vn.Jz[j][nx1 - k - nghostJ-1]; // Left side + Vn.Jz[j][nx1 - 1 + k - nghostJ] = Vn.Jz[j][k - 1 + nghostJ+1]; // Right side } } diff --git a/src/ConstainedTransport.cpp b/src/ConstainedTransport.cpp index e72cd31..f86ea2e 100644 --- a/src/ConstainedTransport.cpp +++ b/src/ConstainedTransport.cpp @@ -41,7 +41,7 @@ std::vector> ConstrainedTransportAverage(const ConservativeV if (OptP == OptionalPhysics::HallResHyper){ // Edge-centered (reconstructions) - std::vector> rho(Cn.ny + 1 - 2 * nghost, std::vector(Cn.nx + 1 - 2 * nghost)); // used in hall effect. Since we only have une mass, m=1 -> n = rho + std::vector> rho(Cn.ny + 1 - 2 * nghost, std::vector(Cn.nx + 1 - 2 * nghost)); // used in hall effect. Since we only have one mass, m=1 -> n = rho std::vector> jx(Cn.ny + 1 - 2 * nghost, std::vector(Cn.nx + 1 - 2 * nghost)); std::vector> jy(Cn.ny + 1 - 2 * nghost, std::vector(Cn.nx + 1 - 2 * nghost)); @@ -56,9 +56,9 @@ std::vector> ConstrainedTransportAverage(const ConservativeV jx[j - nghost][i - nghost] = AVERAGEX(Cn.Jx, i1-1, j1); jy[j - nghost][i - nghost] = AVERAGEY(Cn.Jy, i1, j1-1); - Ez[j - nghost][i - nghost] += (1.0/rho[j - nghost][i - nghost]) * (jx[j - nghost][i - nghost] * By[j - nghost][i - nghost] + jy[j - nghost][i - nghost] * vy[j - nghost][i - nghost]); // Hall contribution - Ez[j - nghost][i - nghost] += pc.eta * (Cn.Jz[j1][i1]); // Resistivity - Ez[j - nghost][i - nghost] -= pc.nu * LAPLACIAN(Cn.Jz, i1, j1, Dx, Dy); // Hyper resistivity + Ez[j - nghost][i - nghost] += (1.0/rho[j - nghost][i - nghost]) * (jx[j - nghost][i - nghost] * By[j - nghost][i - nghost] - jy[j - nghost][i - nghost] * Bx[j - nghost][i - nghost]); // Hall contribution + //Ez[j - nghost][i - nghost] += pc.eta * (Cn.Jz[j1][i1]); // Resistivity + //Ez[j - nghost][i - nghost] -= pc.nu * LAPLACIAN(Cn.Jz, i1, j1, Dx, Dy); // Hyper resistivity } } } diff --git a/src/Interface.cpp b/src/Interface.cpp index 6cae46b..8e3efc8 100644 --- a/src/Interface.cpp +++ b/src/Interface.cpp @@ -31,13 +31,21 @@ 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 (dir == Dir::X) { + f.Bz += (1.0/u.rho) * (Jz * u.Bx - Jx * u.Bz); + //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 += pc.eta * (Jy * u.Bz - Jz * u.By); - f.P -= pc.nu * (LaplJy * u.Bz - LaplJz * u.By); + //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.Bz += - (1.0/u.rho) * (Jy * u.Bz - Jz * u.By); + //f.Bz += - pc.eta * Jx; + //f.Bz -= - pc.nu * LaplJx; + f.P += (1.0/u.rho) * ((u.Bx * Jx + u.By * Jy + u.Bz * Jz)*u.By - (u.Bx * u.Bx + u.By * u.By + u.Bz * u.Bz) * Jy); - f.P += pc.eta * (Jz * u.Bx - Jx * u.Bz); - f.P -= pc.nu * (LaplJz * u.Bx - LaplJx * u.Bz); + //f.P += pc.eta * (Jz * u.Bx - Jx * u.Bz); + //f.P -= pc.nu * (LaplJz * u.Bx - LaplJx * u.Bz); } } } @@ -147,6 +155,8 @@ std::pair, std::pair> ComputeRiemannJ( Interface::Interface() = default; Interface::Interface(const PrimitiveVariables& P_cc /* Assuming ghost cells are added */, int i /* (0 to nx) + nghost */, int j /* (0 to ny) + nghost */, double Dx, double Dy, int nghost, Reconstruction rec, Slope sl, OptionalPhysics OptP, Dir dir) { + // For riemann solver + OP = OptP; if (rec == Reconstruction::Constant) { if (dir == Dir::X) { @@ -234,15 +244,24 @@ Interface::Interface(const PrimitiveVariables& P_cc /* Assuming ghost cells are cfastyR = std::sqrt((c0R*c0R + caR*caR)*0.5 + (std::sqrt((c0R*c0R + caR*caR)*(c0R*c0R + caR*caR) - 4*c0R*c0R*cayR*cayR))*0.5); if (OptP == OptionalPhysics::HallResHyper) { - cwxL = std::abs(uL.Bx) * M_PI / (uL.rho * Dx); - cwyL = std::abs(uL.By) * M_PI / (uL.rho * Dy); - cwxR = std::abs(uL.Bx) * M_PI / (uR.rho * Dx); - cwyR = std::abs(uL.By) * M_PI / (uR.rho * Dy); - - cfastxL += cwxL; - cfastxR += cwxR; - cfastyL += cwyL; - cfastyR += cwyR; + double vwx = M_PI * (std::sqrt(1 + 0.25/(Dx*Dx)) + 0.5/Dx); + double vwy = M_PI * (std::sqrt(1 + 0.25/(Dy*Dy)) + 0.5/Dy); + cwxL = std::sqrt(uL.Bx*uL.Bx + uL.By*uL.By + uL.Bz*uL.Bz) / (uL.rho) * vwx; + 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; + + if(dir == Dir::X){ + SLb = std::min(uL.vx - cfastxL - cwxL, uR.vx - cfastxR - cwxR); + SRb = std::max(uL.vx + cfastxL + cwxL, uR.vx + cfastxR + cwxR); + // For rusanov : + Splusb = std::max(std::abs(uL.vx) + cfastxL + cwxL, std::abs(uR.vx) + cfastxR + cwxR); + }else if(dir == Dir::Y){ + SLb = std::min(uL.vy - cfastyL - cwyL, uR.vy - cfastyR - cwyR); + SRb = std::max(uL.vy + cfastyL + cwyL, uR.vy + cfastyR + cwyR); + // For rusanov : + Splusb = std::max(std::abs(uL.vy) + cfastyL + cwyL, std::abs(uR.vy) + cfastyR + cwyR); + } } // Wave speeds diff --git a/src/Interface.hpp b/src/Interface.hpp index 0c15af1..77b2969 100644 --- a/src/Interface.hpp +++ b/src/Interface.hpp @@ -26,8 +26,10 @@ class Interface { public: ReconstructedValues uL, uR, fL, fR; double SL, SR, Splus; + double SLb, SRb, Splusb; double cfastxL, cfastxR, cfastyL, cfastyR; double cwxL, cwxR, cwyL, cwyR; + OptionalPhysics OP; // For vector initialisation in UCT Interface(); diff --git a/src/PhareMHD.cpp b/src/PhareMHD.cpp index 88d32a5..490986d 100644 --- a/src/PhareMHD.cpp +++ b/src/PhareMHD.cpp @@ -6,11 +6,6 @@ void PhareMHD(const PrimitiveVariables& P0cc, std::string resultDir, int nghost, PhysicalConstants::getInstance().init(cst); PrimitiveVariables P0 = InitialiseGhostCells(P0cc, nghost, bc); // Because of the constructor, J is already the right size - - if (OptP == OptionalPhysics::HallResHyper){ - ComputeJ(P0, Dx, Dy, nghost); - UpdateGhostJ(P0, nghost, bc); - } ConservativeVariables U0(P0); diff --git a/src/RiemannSolver.cpp b/src/RiemannSolver.cpp index 7459b0f..02da2cf 100644 --- a/src/RiemannSolver.cpp +++ b/src/RiemannSolver.cpp @@ -1,17 +1,56 @@ #include "RiemannSolver.hpp" ReconstructedValues RusanovRiemannSolver(const Interface& inter){ - return 0.5 * (inter.fL + inter.fR) - 0.5 * inter.Splus * (inter.uR - inter.uL); + ReconstructedValues Frus; + Frus = 0.5 * (inter.fL + inter.fR) - 0.5 * inter.Splusb * (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); + Frus.By = 0.5 * (inter.fL.By + inter.fR.By) - 0.5 * inter.Splusb * (inter.uR.By - inter.uL.By); + Frus.Bz = 0.5 * (inter.fL.Bz + inter.fR.Bz) - 0.5 * inter.Splusb * (inter.uR.Bz - inter.uL.Bz); + 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)<<" "<= 0){ - return inter.fL; + ReconstructedValues fhll; + if(inter.SL > 0){ + fhll = inter.fL; } - else if((inter.SL <= 0) && (inter.SR >= 0)){ - return (inter.SR * inter.fL - inter.SL * inter.fR + inter.SL * inter.SR * (inter.uR - inter.uL)) / (inter.SR - inter.SL); + else if (inter.SR < 0){ + fhll = inter.fR; } else{ - return inter.fR; + fhll = (inter.SR * inter.fL - inter.SL * inter.fR + inter.SL * inter.SR * (inter.uR - inter.uL)) / (inter.SR - inter.SL); } + + + if (inter.OP == OptionalPhysics::HallResHyper){ + if(inter.SLb > 0){ + fhll.Bx = inter.fL.Bx; + fhll.By = inter.fL.By; + fhll.Bz = inter.fL.Bz; + fhll.P = inter.fL.P; + } + else if (inter.SRb < 0) { + fhll.Bx = inter.fR.Bx; + fhll.By = inter.fR.By; + fhll.Bz = inter.fR.Bz; + fhll.P = inter.fR.P; + } + else { + fhll.Bx = (inter.SRb * inter.fL.Bx - inter.SLb * inter.fR.Bx + inter.SLb * inter.SRb * (inter.uR.Bx - inter.uL.Bx)) / (inter.SRb - inter.SLb); + fhll.By = (inter.SRb * inter.fL.By - inter.SLb * inter.fR.By + inter.SLb * inter.SRb * (inter.uR.By - inter.uL.By)) / (inter.SRb - inter.SLb); + fhll.Bz = (inter.SRb * inter.fL.Bz - inter.SLb * inter.fR.Bz + inter.SLb * inter.SRb * (inter.uR.Bz - inter.uL.Bz)) / (inter.SRb - inter.SLb); + fhll.P = (inter.SRb * inter.fL.P - inter.SLb * inter.fR.P + inter.SLb * inter.SRb * (inter.uR.P - inter.uL.P)) / (inter.SRb - inter.SLb); + } + } + + return fhll; } \ No newline at end of file diff --git a/src/TimeIntegrator.cpp b/src/TimeIntegrator.cpp index b34bfe8..8fafbc4 100644 --- a/src/TimeIntegrator.cpp +++ b/src/TimeIntegrator.cpp @@ -19,6 +19,17 @@ ConservativeVariables Euler(ConservativeVariables& Un, double Dx, double Dy, dou if (OptP == OptionalPhysics::HallResHyper) { ComputeJ(Un, Dx, Dy, nghost); UpdateGhostJ(Un, nghost, bc); +/* static int i = 0; + std::ostringstream filenamex; + filenamex << "whislerwaveres/" <<"Jx_"<