How can I prevent or mitigate this issue (I have already tried adjusting the time step, grid spacing)?

How can I prevent or mitigate this issue? I have already tried adjusting the time step, grid spacing.

import numpy as np from scipy.sparse import diags from scipy.sparse.linalg import spsolve # --------------------------------------------------------------- # Simulation setup # --------------------------------------------------------------- dt = 0.1/nu_ei0 total_time = 30/nu_ei0 num_steps = int(total_time / dt) n_e = density_profile # initial condition Te_int = temperature_profile # initial temperature profile # iteration control max_iter = 1000 alpha = 0.01 Te_history = np.zeros((num_steps+1, nx)) Qe_history_current = np.zeros((num_steps, nx-1)) Qe_history_previous = np.zeros((num_steps, nx-1)) Te_history[0, :] = Te_int.copy() for n in range(num_steps): print(f"\n================ TIME STEP {n+1}/{num_steps} ================") if n == 0: Te_n = Te_int.copy() else: Te_n = Te.copy() Te_k = Te_n.copy() divQ_SH_prev = np.zeros_like(Te_n) S_nl_corr_prev = np.zeros_like(Te_n) # ----------------------------------------------------------- # Inner nonlinear iteration loop # ----------------------------------------------------------- for k in range(1, max_iter + 1): #print(f" Iteration {k} ...") # STEP 1: compute Te_k Te_face = (Te_k[:-1] + Te_k[1:])/2 kappa_i = ka * Te_face**(2.5) k_right = np.zeros(nx) k_left = np.zeros(nx) k_right[:-1] = kappa_i k_left[1:] = k_right[:-1] n_interior = nx - 2 rhs = np.zeros(n_interior) a = np.zeros(nx) b = np.zeros(nx) c = np.zeros(nx) for i in range(1, nx-1): a[i] = -k_left[i] / dx**2 b[i] = (k_left[i] + k_right[i]) / dx**2 + (1.5 * n_e[i] / dt) c[i] = -k_right[i] / dx**2 # reflective BC corrections --- b[1] += a[1]; a[1] = 0.0 b[nx-2] += c[nx-2]; c[nx-2] = 0.0 # diagonals for sparse matrix lower_diag = np.zeros(n_interior-1) main_diag = np.zeros(n_interior) upper_diag = np.zeros(n_interior-1) for j, i in enumerate(range(1, nx-1)): if i > 1: lower_diag[j-1] = a[i] main_diag[j] = b[i] if i < nx-2: upper_diag[j] = c[i] rhs[j] = (1.5 * n_e[i] * Te_n[i] / dt) + S_nl_corr_prev[i] # build sparse matrix directly A_sparse = diags([lower_diag, main_diag, upper_diag], [-1, 0, 1], format="csc") # print(A_sparse) Te_new_bulk = spsolve(A_sparse, rhs) # reflective boundary for temperature --- Te_new = np.zeros(nx) Te_new[1:-1] = Te_new_bulk Te_new[0] = Te_new[1] Te_new[-1] = Te_new[-2] #Te_newkk = 2 * Te_new[:-1] * Te_new[1:] / (Te_new[:-1] + Te_new[1:]) Te_new_face = (Te_new[:-1] + Te_new[1:])/2 # STEP 2: compute ∇·Q_SH kappa_face = ka * Te_new_face**2.5 # Cao flux_face = - kappa_face * (Te_new[1:] - Te_new[:-1]) / dx #plt.plot(kappa_face1) #plt.show() divQ_SH_now = np.zeros_like(Te_new) divQ_SH_now[1:-1] = (flux_face[1:] - flux_face[:-1]) / dx divQ_SH_now[0] = divQ_SH_now[1] divQ_SH_now[-1] = divQ_SH_now[-2] # STEP 3: convergence check (skip the k=1 interation) if k > 1: flux_diff = np.abs(divQ_SH_now - divQ_SH_prev) tol_array = alpha * (1.5 * n_e * Te_new / dt) if np.all(flux_diff <= tol_array): print(f" ✅ Converged at iteration {k}") break # STEP 4 & 5: compute nonlocal heat flux and its divergence Q_nl, divQ_nl = Limiter_solver(Te_new) # STEP 7: update nonlocal correction S_nl_corr = divQ_SH_now - divQ_nl # update for next iteration divQ_SH_prev = divQ_SH_now.copy() S_nl_corr_prev = S_nl_corr.copy() #Te_k = 0.5 * Te_k + 0.5 * Te_new Te_k = Te_new.copy() # update temperature for next time step Te_history[n+1, :] = Te_new.copy() Te = Te_new.copy() Qe_history_previous[n, :] = Q_nl.copy() Q_nlt, divQ_nlt = Limiter_solver(Te) Qe_history_current[n, :] = Q_nlt.copy() #plt.plot(Te) # optional: diagnostic #total_energy = np.trapz(1.5 * n_e * Te, x) #print(f" Total energy = {total_energy:.3e} J/m²") print(f" Max Te = {Te.max():.3f} eV, Min Te = {Te.min():.3f} eV") print("\nSimulation complete.") 

import numpy as np from scipy.sparse import diags from scipy.sparse.linalg import spsolve from scipy.integrate import quad # --------------------------- # Spitzer-Härm flux # ---------------------------  import numpy as np import matplotlib.pyplot as plt from math import sqrt, log, pi, exp # Constants c0 = 299792458 # the vacuum light velocity in unit of m/s me = 9.10938356e-31 # the electron mass in kg epsilon0 = 8.8541878128e-12 # the dielectric constant in vacuum (C^2⋅kg^−1⋅m^−3⋅s^2) qe = 1.60217662e-19 # the electron charge in unit of coulomb (C) K = 1.783e-33*c0*c0 # in unit of kg*(m/s)^2/(keV), i.e., no units kB = 1.380649e-23 # in the unit of J/K mu00 = 1.25663706127e-6 # vacuum permeability in the unit of (N⋅(C/s)^−2) # OR mu0 = 1/(c0**2*epsilon0) # the relation connected by the light velocity Te_central = 1000 # the temperature in the central wall, in the unit of eV Te_per = 4000 Ne = 1.5e20 # the central plasma density 1/m^3 omegap = sqrt(Ne*qe**2/(epsilon0*me)) # plasma frequency (1/s) kp = omegap/c0 # plasma wavenumber p = Ne*kB*Te_central*11604.51812 # the thermal pressure (J/m^3) in kept as constant, 1eV = 11604.51812 K TeRatio = 2 # the temperature ratio between hot and cold Teh = 2*Te_central/(1+TeRatio)*TeRatio # the temperature in the hot wall, in the unit of eV Tel = 2*Te_central/(1+TeRatio) # the temperature in the cold wall, in the unit of eV vthMean_SI = sqrt(Te_central/1000*1*K/me) # thermal velocity in the hot wall, in the unit of SI (m/s) vthh_SI = sqrt(Teh/1000*1*K/me) # thermal velocity in the hot wall, in the unit of SI (m/s) vthl_SI = sqrt(Tel/1000*1*K/me) # thermal velocity in the cold wall, in the unit of SI (m/s) Z = 1 # ion charge state LambdaCoulomb1 = 23 - log((Ne*1.e-6)**(1/2)*Z*(Te_central)**(-3/2)) # Ti*me/mi < Te < 10z^2 eV LambdaCoulomb2 = 24 - log((Ne*1.e-6)**(1/2)*(Te_central)**(-1)) # Ti*me/mi < 10z^2 eV < Te if Te_central <= 10*Z**2: CoulombLogarithm = LambdaCoulomb1 else: CoulombLogarithm = LambdaCoulomb2 nu_ei0 = Z*Ne*qe**4*CoulombLogarithm/(4*pi*epsilon0**2*me**2)/(vthMean_SI**3) mfp0 = vthMean_SI/nu_ei0 # centrel mean free path in m LL = 160*mfp0 # simulation box length Xmin = -LL/2 Xmax = LL/2 ll = 8.44*mfp0 # Create spatial grid num_points = 800 central = int(num_points / 2) # Ensures an integer value x = np.linspace(Xmin, Xmax, num_points) nx = len(x) dx = x[1]-x[0] temperature_profile = Te_central + Te_per*np.exp(-x**2/ll**2) # temperature profile in eV thermal_velocity_profile = np.sqrt(temperature_profile / 1000 * 1 * K / me) # thermal velocity profile in m/s #density_profile = Te_central*Ne/temperature_profile # electron density 1/m^3 density_profile = np.full(len(x), Ne) # electron density 1/m^3 Z_factor = (Z + 0.24) / (Z + 4.2) cont_nu_ei = Z*density_profile*qe**4* CoulombLogarithm/(4*pi*epsilon0**2*me**2) nu_ei_local = cont_nu_ei/(thermal_velocity_profile**3) mfp_profile = thermal_velocity_profile/nu_ei_local # mean free path in m (SI units) mfp_sherlock = mfp_profile*np.sqrt(Z)*np.sqrt(Z_factor) # mean free path in m (SI units) lambda_lmv = 32 * np.sqrt(Z + 1) * mfp_profile dT_dx = np.gradient(temperature_profile, x) const_factor = 128 / (3 * np.pi) * Z_factor * (3 * np.pi) / np.sqrt(2 * np.pi) kappa = const_factor * 4*pi*epsilon0**2*me**2 / (Z*qe**4*CoulombLogarithm) * np.sqrt(temperature_profile/1000*1*K/me)**5 ka = const_factor * 4*pi*epsilon0**2*me**2 / (Z*qe**4*CoulombLogarithm) * np.sqrt(1/1000*1*K/me)**5 q_sh = -kappa * dT_dx q_free = density_profile * temperature_profile * thermal_velocity_profile def SH_solver(Te): Te_center = (Te[:-1] + Te[1:]) /2 kappa_face = ka * Te_center**2.5 # Cao Q_nl = -kappa_face * (Te[1:] - Te[:-1]) / dx * 0.01 divQ_nl = np.zeros_like(Te) divQ_nl[1:-1] = (Q_nl[1:] - Q_nl[:-1]) / dx divQ_nl[0] = divQ_nl[1] divQ_nl[-1] = divQ_nl[-2] return Q_nl, divQ_nl import numpy as np from scipy.sparse import diags from scipy.sparse.linalg import spsolve # --------------------------------------------------------------- # Simulation setup # --------------------------------------------------------------- dt = 0.1/nu_ei0 total_time = 30/nu_ei0 num_steps = int(total_time / dt) n_e = density_profile # initial condition Te_int = temperature_profile # initial temperature profile # iteration control max_iter = 1000 alpha = 0.01 Te_history = np.zeros((num_steps+1, nx)) Qe_history_current = np.zeros((num_steps, nx-1)) Qe_history_previous = np.zeros((num_steps, nx-1)) Te_history[0, :] = Te_int.copy() for n in range(num_steps): print(f"\n================ TIME STEP {n+1}/{num_steps} ================") if n == 0: Te_n = Te_int.copy() else: Te_n = Te.copy() Te_k = Te_n.copy() divQ_SH_prev = np.zeros_like(Te_n) S_nl_corr_prev = np.zeros_like(Te_n) # ----------------------------------------------------------- # Inner nonlinear iteration loop # ----------------------------------------------------------- for k in range(1, max_iter + 1): #print(f" Iteration {k} ...") # STEP 1: compute Te_k Te_face = (Te_k[:-1] + Te_k[1:])/2 kappa_i = ka * Te_face**(2.5) k_right = np.zeros(nx) k_left = np.zeros(nx) k_right[:-1] = kappa_i k_left[1:] = k_right[:-1] n_interior = nx - 2 rhs = np.zeros(n_interior) a = np.zeros(nx) b = np.zeros(nx) c = np.zeros(nx) for i in range(1, nx-1): a[i] = -k_left[i] / dx**2 b[i] = (k_left[i] + k_right[i]) / dx**2 + (1.5 * n_e[i] / dt) c[i] = -k_right[i] / dx**2 # reflective BC corrections --- b[1] += a[1]; a[1] = 0.0 b[nx-2] += c[nx-2]; c[nx-2] = 0.0 # diagonals for sparse matrix lower_diag = np.zeros(n_interior-1) main_diag = np.zeros(n_interior) upper_diag = np.zeros(n_interior-1) for j, i in enumerate(range(1, nx-1)): if i > 1: lower_diag[j-1] = a[i] main_diag[j] = b[i] if i < nx-2: upper_diag[j] = c[i] rhs[j] = (1.5 * n_e[i] * Te_n[i] / dt) + S_nl_corr_prev[i] # build sparse matrix directly A_sparse = diags([lower_diag, main_diag, upper_diag], [-1, 0, 1], format="csc") # print(A_sparse) Te_new_bulk = spsolve(A_sparse, rhs) # reflective boundary for temperature --- Te_new = np.zeros(nx) Te_new[1:-1] = Te_new_bulk Te_new[0] = Te_new[1] Te_new[-1] = Te_new[-2] #Te_newkk = 2 * Te_new[:-1] * Te_new[1:] / (Te_new[:-1] + Te_new[1:]) Te_new_face = (Te_new[:-1] + Te_new[1:])/2 # STEP 2: compute ∇·Q_SH kappa_face = ka * Te_new_face**2.5 # Cao flux_face = - kappa_face * (Te_new[1:] - Te_new[:-1]) / dx #plt.plot(kappa_face1) #plt.show() divQ_SH_now = np.zeros_like(Te_new) divQ_SH_now[1:-1] = (flux_face[1:] - flux_face[:-1]) / dx divQ_SH_now[0] = divQ_SH_now[1] divQ_SH_now[-1] = divQ_SH_now[-2] # STEP 3: convergence check (skip the k=1 interation) if k > 1: flux_diff = np.abs(divQ_SH_now - divQ_SH_prev) tol_array = alpha * (1.5 * n_e * Te_new / dt) if np.all(flux_diff <= tol_array): print(f" ✅ Converged at iteration {k}") break # STEP 4 & 5: compute nonlocal heat flux and its divergence Q_nl, divQ_nl = SH_solver(Te_new) # STEP 7: update nonlocal correction S_nl_corr = divQ_SH_now - divQ_nl # update for next iteration divQ_SH_prev = divQ_SH_now.copy() S_nl_corr_prev = S_nl_corr.copy() #Te_k = 0.5 * Te_k + 0.5 * Te_new Te_k = Te_new.copy() # update temperature for next time step Te_history[n+1, :] = Te_new.copy() Te = Te_new.copy() Qe_history_previous[n, :] = Q_nl.copy() Q_nlt, divQ_nlt = SH_solver(Te) Qe_history_current[n, :] = Q_nlt.copy() #plt.plot(Te) # optional: diagnostic #total_energy = np.trapz(1.5 * n_e * Te, x) #print(f" Total energy = {total_energy:.3e} J/m²") print(f" Max Te = {Te.max():.3f} eV, Min Te = {Te.min():.3f} eV") print("\nSimulation complete.")