import numpy as np
from adios2 import FileReader
[docs]
class AnalyticDiffusionProfiles:
"""
Creates diffusion profile data using analytic functions.
Attributes:
density : A 3D array (shape (n_species, n_samples, n_surf)) of the diffusion density.
flow : A 3D array (shape (n_species, n_samples, n_surf) of the parallel mean flow.
temp : A 3D array (shape (n_species, n_samples, n_surf)) of temperature (in eV).
n_samples : An integer (n_samples) representing the number of time samples.
n_species : An integer (default 2) representing the number of particle species.
n_surf : An integer (n_surf) representing the number of ψ-surfaces.
psi : A 1D array (shape (n_surf,)) containing the poloidal magnetic flux values.
steps : A 1D array (shape (n_samples,)) of time step indices.
time : A 1D array (shape (n_samples,)) of simulation times (in seconds).
"""
def __init__(self, psi, time, S, Vol, grad_avg, psi_norm):
"""
Alternative constructor that generates a time series using analytic formulas.
Inputs:
psi : 1D numpy array of psi values (assumed monotonically increasing).
time : 1D numpy array of time values.
S : 1D numpy array of the surface area of the flux-surfaces.
V : 1D numpy array of V(psi) at the psi grid points.
grad_avg : <|grad(psi)|>
psi_norm : Psi value on the separatrix.
The initial profiles (at t = time[0]) for density (n), flow (u), and temperature (T)
are defined by:
Y(psi, t=0) = Y₀ exp[-w*k * tanh((psi - ψ₀) / w)]
with:
w = 0.3, k = 2.0, ψ₀ = 0.5.
For n, Y₀ = n₀ = 1e19;
for u, Y₀ = u₀ = 1e4;
for T, Y₀ = T₀ = 1e3.
The transport coefficients are defined as:
Y(psi) = Y₀ exp[-((psi - ψ₀) / w_coef)²]
with ψ₀ = 0.5, w_coef = 0.2; for D, eta, chi use Y₀ = 1.0, and for vₚ, use Y₀ = -0.2.
The transport equations are integrated in time via forward Euler:
Density:
∂ₜ n = 1/Vprime * d/dψ { S*(D dn/dψ - vₚ n) }
Flow:
∂ₜ u = 1/Vprime * d/dψ { S*(eta du/dψ) } + (D+eta)/n * (dn/dψ)*(du/dψ) - vₚ * (du/dψ)
Temperature:
∂ₜ T = (2/3)/Vprime * d/dψ { S*(chi dT/dψ) }
+ (2T/(3n))/Vprime * d/dψ { S*(D dn/dψ - vₚ n) }
+ (5D+2chi)/(3n)*(dn/dψ)*(dT/dψ)
+ (2*eta/3)*(du/dψ)² - (5/3)*vₚ*(dT/dψ)
Here we assume S(psi) = 1 and m = e = 1.
"""
# Number of psi grid points and number of time steps
self.time = np.array(time,dtype=np.float64)
self.psi = np.array(psi,dtype=np.float64)
self.n_surf = len(psi)
self.n_samples = len(time)
self.steps = 1
self.n_species = 2
self.steps = np.arange(self.n_samples)
self.density = np.zeros([self.n_species, self.n_samples, self.n_surf], order='C', dtype=np.float64)
self.flow = np.zeros([self.n_species, self.n_samples, self.n_surf], order='C', dtype=np.float64)
self.temp = np.zeros([self.n_species, self.n_samples, self.n_surf], order='C', dtype=np.float64)
nt = self.n_samples
n_psi = self.n_surf
# Parameters for the profiles at t=0
w = 0.3
k = 2.0
psi0 = 0.5
n0 = 1e19
u0 = 1e4
T0 = 1e3
# Function for initial profile for n, u, or T.
def Y_profile(Y0, psi, w, k, psi0):
return Y0 * np.exp(-w * k * np.tanh((psi - psi0) / w))
# Initial profiles (at t = time[0])
n0_profile = Y_profile(n0, self.psi/psi_norm, w, k, psi0)
u0_profile = Y_profile(u0, self.psi/psi_norm, w, k, psi0)
T0_profile = Y_profile(T0, self.psi/psi_norm, w, k, psi0)
# Parameters for transport coefficients
w_coef = 0.2
def trans_coeff(Y0, psi, psi0, w_coef):
return Y0 * np.exp(-((psi - psi0)/w_coef)**2)
D = trans_coeff( 1.0, self.psi/psi_norm, psi0, w_coef)
vp = trans_coeff(-0.2, self.psi/psi_norm, psi0+0.1, w_coef) # note the negative for v_p
eta = trans_coeff( 1.0, self.psi/psi_norm, psi0, w_coef)
chi = trans_coeff( 1.0, self.psi/psi_norm, psi0, w_coef)
# Helper: centered finite difference for a nonuniform psi grid.
def dfdpsi(f, psi):
N = len(psi)
df = np.zeros_like(f)
# Interior points: centered difference
for i in range(1, N-1):
dpsi = psi[i+1] - psi[i-1]
df[i] = (f[i+1] - f[i-1]) / dpsi
# Forward difference at the first point
df[0] = (f[1] - f[0]) / (psi[1] - psi[0])
# Backward difference at the last point
df[-1] = (f[-1] - f[-2]) / (psi[-1] - psi[-2])
return df
Vprime = dfdpsi(Vol,self.psi)
# Preallocate time series arrays (each: nt x n_psi)
n_series = np.zeros((nt, n_psi), order='C', dtype=np.float64)
u_series = np.zeros((nt, n_psi), order='C', dtype=np.float64)
T_series = np.zeros((nt, n_psi), order='C', dtype=np.float64)
# Set initial profiles at time[0]
n_series[0, :] = n0_profile
u_series[0, :] = u0_profile
T_series[0, :] = T0_profile
m = 3.34e-27
e_val = 1.6022e-19
# Time integration using a simple forward Euler method
for t_idx in range(nt - 1):
dt = self.time[t_idx+1] - self.time[t_idx]
# CURRENT profiles
n_curr = n_series[t_idx, :]
u_curr = u_series[t_idx, :]
T_curr = T_series[t_idx, :]
# Compute spatial derivatives needed for the transport equations:
dn_dpsi = dfdpsi(n_curr, self.psi) * grad_avg
du_dpsi = dfdpsi(u_curr, self.psi) * grad_avg
dT_dpsi = dfdpsi(T_curr, self.psi) * grad_avg
# ---------------------------
# Density Equation:
# Flux: F_n = D * (dn/dpsi) - vp * n.
F_n = S*(D * dn_dpsi - vp * n_curr)
dF_n_dpsi = dfdpsi(F_n, psi)
# ∂t n = 1/Vprime * dF_n/dpsi (S is 1)
dt_n = dF_n_dpsi / Vprime
# ---------------------------
# Flow Equation:
# Flux: F_u = eta * (du/dpsi)
F_u = S * eta * du_dpsi
dF_u_dpsi = dfdpsi(F_u, psi)
term_flow = (D + eta) / n_curr * dn_dpsi * du_dpsi # pointwise multiplication
dt_u = dF_u_dpsi / Vprime + term_flow - vp * du_dpsi
# ---------------------------
# Temperature Equation:
# Term 1: (2/3)/Vprime * d/dpsi { chi*dT/dpsi }
F_T1 = S * chi * dT_dpsi
dF_T1_dpsi = dfdpsi(F_T1, psi)
term_T1 = (2/3) * dF_T1_dpsi / Vprime
# Term 2: (2T/(3n))/Vprime * d/dpsi { D*dn/dpsi - vp*n }
F_T2 = S * (D * dn_dpsi - vp * n_curr)
dF_T2_dpsi = dfdpsi(F_T2, psi)
term_T2 = (2 * T_curr / (3 * n_curr)) * dF_T2_dpsi / Vprime
# Term 3: (5D+2chi)/(3n)*(dn/dpsi)*(dT/dpsi)
term_T3 = ((5 * D + 2 * chi) / (3 * n_curr)) * dn_dpsi * dT_dpsi
# Term 4: (2*eta/3)*(du/dpsi)^2 -- note that m/e = 1.
term_T4 = (2 * eta * m / (3 * e_val)) * (du_dpsi**2)
# Term 5: - (5/3) * vp * dT/dpsi
term_T5 = - (5/3) * vp * dT_dpsi
dt_T = term_T1 + term_T2 + term_T3 + term_T4 + term_T5
# Update profiles with forward Euler:
n_series[t_idx+1, :] = n_curr + dt * dt_n
u_series[t_idx+1, :] = u_curr + dt * dt_u
T_series[t_idx+1, :] = T_curr + dt * dt_T
for isp in range(self.n_species):
self.density[isp, :, :] = n_series
self.flow[isp, :, :] = u_series
self.temp[isp, :, :] = T_series
def __repr__(self):
return (f"DiffusionProfiles(n_species={self.n_species}, n_samples={self.n_samples}, "
f"n_surf={self.n_surf}, density.shape={self.density.shape}, "
f"psi.shape={self.psi.shape})")