From 43cce3593a1a1296b03e4f4bd3e15d014335160c Mon Sep 17 00:00:00 2001 From: Marco Sertoli Date: Mon, 7 Aug 2023 16:07:02 +0100 Subject: [PATCH 01/16] Initial commit adding draft neutral beam class and function to compute fast ion slowing down and fast particle pressure --- indica/models/neutral_beam.py | 210 +++++++++++++++++++++ indica/operators/slowingdown.py | 317 ++++++++++++++++++++++++++++++++ 2 files changed, 527 insertions(+) create mode 100644 indica/models/neutral_beam.py create mode 100644 indica/operators/slowingdown.py diff --git a/indica/models/neutral_beam.py b/indica/models/neutral_beam.py new file mode 100644 index 00000000..0556f1e0 --- /dev/null +++ b/indica/models/neutral_beam.py @@ -0,0 +1,210 @@ +from typing import Tuple + +import numpy as np + +from ..converters.line_of_sight import LineOfSightTransform + + +analytical_beam_defaults = { + "element": "H", + "amu": int(1), + "energy": 25.0 * 1e3, + "power": 500 * 1e3, + "fractions": (0.7, 0.1, 0.2, 0.0), + "divergence": (14 * 1e-3, 14e-3), + "width": (0.025, 0.025), + "location": (-0.3446, -0.9387, 0.0), + "direction": (0.707, 0.707, 0.0), + "focus": 1.8, +} + + +class NeutralBeam: + def __init__(self, name: str, use_defaults=True, **kwargs): + """Set beam parameters, initialisation""" + self.name = name # Beam name + + # Use beam defaults + if use_defaults: + for (prop, default) in analytical_beam_defaults.items(): + setattr(self, prop, kwargs.get(prop, default)) + else: + for (prop, default) in analytical_beam_defaults.items(): + setattr(self, prop, None) + return + + # Set line_of_sight transform for centre of beam-line + self.set_los_transform() + + # Set Attenuator + self.attenuator = None + + # Print statements for debugging + print(f"Beam = {self.name}") + print(f"Energy = {self.energy} electron-volts") + print(f"Power = {self.power} watts") + print(f"Fractions (full, 1/2, 1/3, imp) = {self.fractions} %") + print(f"divergence (x, y) = {self.divergence} rad") + print(f"width (x, y) = {self.width} metres") + + def set_los_transform( + self, + machine_dimensions: Tuple[Tuple[float, float], Tuple[float, float]] = ( + (0.175, 1.0), + (-2.0, 2.0), + ), + ): + self.transform = LineOfSightTransform( + np.array([self.location[0]]), + np.array([self.location[1]]), + np.array([self.location[2]]), + np.array([self.direction[0]]), + np.array([self.direction[1]]), + np.array([self.direction[2]]), + name=f"{self.name}_los", + dl=0.01, + machine_dimensions=machine_dimensions, + ) + + def run_BBNBI(self): + print("Add code to run BBNBI") + + # def run_analytical_beam( + # self, + # x_dash: LabeledArray = DataArray(np.linspace(-0.25, 0.25, 101, dtype=float)), + # y_dash: LabeledArray = DataArray(np.linspace(-0.25, 0.25, 101, dtype=float)), + # ): + # """Analytical beam based on double gaussian formula""" + # + # # Calculate Attenuator + # x2 = self.transform.x2 + # attenuation_factor = np.ones_like( + # x2 + # ) # Replace with Attenuation Object in future, interact with plasma + # + # # Calculate beam velocity + # # v_beam = self.beam_velocity() + # # e = 1.602 * 1e-19 + # + # # Neutral beam + # n_x = 101 + # n_y = 101 + # n_z = len(attenuation_factor) + # nb_dash = np.zeros((n_z, n_x, n_y), dtype=float) + # # for i_z in range(n_z): + # # for i_y in range(n_y): + # # y_here = y_dash[i_y] + # # exp_factor = np.exp( + # # -(x_dash**2 / self.width[0]**2) - (y_here**2 / self.width[1]**2) + # # ) + # # nb_dash[i_z, :, i_y] = \ + # # self.power * attenuation_factor[i_z] * exp_factor / \ + # # (np.pi * self.energy * e * np.prod(self.width) * v_beam) + # + # # mesh-grid of beam cross section coordinates + # z_dash = x2 + # X_dash, Y_dash, Z_dash = np.meshgrid(x_dash, y_dash, z_dash) + # R_dash = np.sqrt(X_dash**2 + Y_dash**2) + # T_dash = np.arctan2(Y_dash, X_dash) + # + # x_transform = self.transform.x + # y_transform = self.transform.y + # z_transform = self.transform.z + # r_transform = np.sqrt(self.transform.x**2 + self.transform.y**2) + # theta_transform = np.arctan2( + # y_transform[1] - y_transform[0], x_transform[1] - x_transform[0] + # ) + # phi_transform = np.arctan2( + # z_transform[1] - z_transform[0], r_transform[1] - r_transform[0] + # ) + # theta_n = theta_transform + (np.pi / 2) + # phi_n = phi_transform + (np.pi / 2) + # + # xd = np.zeros_like(X_dash) + # yd = np.zeros_like(Y_dash) + # zd = np.zeros_like(Z_dash) + # for i_z in range(len(x2)): + # delta_X_dash = R_dash[:, :, i_z] * np.cos(T_dash[:, :, i_z]) + # delta_Y_dash = R_dash[:, :, i_z] * np.sin(T_dash[:, :, i_z]) + # + # xd[:, :, i_z] = x_transform[i_z].data + delta_X_dash * np.cos( + # theta_n.data + # ) # + X_dash[:, :, i_z]*np.tan(theta_n) + # yd[:, :, i_z] = y_transform[i_z].data + delta_X_dash * np.sin( + # theta_n.data + # ) # + Y_dash[:, :, i_z]*np.tan(theta_n) + # zd[:, :, i_z] = z_transform[i_z].data + delta_Y_dash * np.sin(phi_n.data) + # + # plt.figure() + # for i_z in range(len(x2)): + # plt.plot(xd[:, :, i_z].flatten(), yd[:, :, i_z].flatten(), "r.") + # plt.plot(self.transform.x, self.transform.y, "k") + # plt.axis("equal") + # plt.show(block=True) + # + # # print(X_dash) + # # print(np.shape(X_dash)) + # # print("aa" ** 2) + # + # # # Interpolate over tok-grid + # # x_tok = DataArray(np.linspace(-1.0, 1.0, 101, dtype=float)) + # # y_tok = DataArray(np.linspace(-1.0, 1.0, 101, dtype=float)) + # # z_tok = DataArray(np.linspace(-0.5, 0.5, 51, dtype=float)) + # # points = (z_tok.data, x_tok.data, y_tok.data) + # + # # X_tok, Y_tok, Z_tok = np.meshgrid(x_tok, y_tok, z_tok) + # # X_tok = X_tok.flatten() + # # Y_tok = Y_tok.flatten() + # # Z_tok = Z_tok.flatten() + # # nb_tok = np.zeros_like(X_tok) + # # for i in range(len(X_tok)): + # # print(f'{i}out of {len(X_tok)}') + # # point = np.array([Z_tok[i], X_tok[i], Y_tok[i]]) + # # nb_tok[i] = interpn(points, nb_dash, point, method='linear') + # + # # print("aa" ** 2) + # + # if True: + # + # plt.figure() + # plt.plot(x_dash, np.sum(nb_dash[0, :, :], axis=1)) + # + # plt.figure() + # plt.contour(x_dash, y_dash, nb_dash[0, :, :], 100) + # plt.xlabel("X (m)") + # plt.ylabel("Y (m)") + # plt.title("Initial beam cross section") + # plt.show(block=True) + + def beam_velocity(self): + return 4.38 * 1e5 * np.sqrt(self.energy * 1e-3 / float(self.amu)) + + def set_energy(self, energy: float): + self.energy = energy + + def set_power(self, power: float): + self.power = power + + def set_divergence(self, divergence: tuple): + self.divergence = divergence + + def set_fractions(self, fractions: tuple): + self.fractions = fractions + + def set_width(self, width: tuple): + self.width = width + + def set_element(self, element: str): + self.element = element + if self.element == "H": + self.amu = int(1) + elif self.element == "D": + self.amu = int(2) + else: + raise ValueError + + def set_location(self, location: tuple): + self.location = location + + def set_direction(self, direction: tuple): + self.direction = direction diff --git a/indica/operators/slowingdown.py b/indica/operators/slowingdown.py new file mode 100644 index 00000000..81977473 --- /dev/null +++ b/indica/operators/slowingdown.py @@ -0,0 +1,317 @@ +import matplotlib.pyplot as plt +import numpy as np +from scipy.special import erf +import scipy.interpolate as interp + +epsilon0 = 8.854e-12 +hbar = 1.054e-34 + + +def coulomb_log(ma, qa, va, mb, qb, nb, Tb): + """ + Calculate coulomb logarithm. + + ma: Test particle mass (kg) + qa: Test particle charge (C) + va: Test particle velocity (m/s) + mb: Background species mass (kg) + qb: Background species charge (C) + nb: Background species density (1/m3) + Tb: Background species temperature (1/m3) + """ + debyeLength = np.sqrt(epsilon0 / np.sum(nb * qb * qb / Tb)) + vbar = va * va + 2 * Tb / mb + mr = ma * mb / (ma + mb) + bcl = np.abs(qa * qb / (4 * np.pi * epsilon0 * mr * vbar)) + bqm = np.abs(hbar / (2 * mr * np.sqrt(vbar))) + + clog = np.zeros(len(mb)) + + for i in range(len(clog)): + if bcl[i] > bqm[i]: + clog[i] = np.log(debyeLength / bcl[i]) + else: + clog[i] = np.log(debyeLength / bqm[i]) + + return clog + + +def slowingdown(ma, qa, Ea, mb, qb, nb, Tb, Nmc=10, dt=1): + """ + Calculate slowing-down for a fast particle. Collision operator adopted from + Hirvijoki et al 2014 CPC. + + ma: Fast particle mass (kg) + qa: Fast particle charge (C) + Ea: Fast particle energy (J) + mb: Background species mass (kg) + qb: Background species charge (C) + nb: Background species density (1/m3) + Tb: Background species temperature (eV) + dt: Timestep for output (s) + """ + max_tslow = 0.1 + + tv = np.arange(0, max_tslow, dt) + nfast = np.zeros(len(tv)) + pressure = np.zeros(len(tv)) + + for i in range(Nmc): + va = np.sqrt(2 * Ea / ma) + vb = np.sqrt(2 * Tb / mb) + + t = 0 + dt = 1e-3 + + while va > np.sqrt(2 * Tb / ma): + t += dt + mu0 = ( + erf(va / vb) - 2 * va / vb * np.exp(-va / vb * va / vb) / np.sqrt(np.pi) + ) / (va / vb * va / vb) + mu1 = erf(va / vb) - 0.5 * mu0 + clog = coulomb_log(ma, qa, va, mb, qb, nb, Tb) + Cb = nb * qa * qa * qb * qb * clog / (4 * np.pi * epsilon0 * epsilon0) + F = np.sum(-(1 / mb + 1 / ma) * Cb * mu0 / (ma * vb * vb)) + Dpara = np.sum(Cb * mu0 / (2 * ma * ma * va)) + Dperp = np.sum(Cb * mu1 / (2 * ma * ma * va)) + + dW = np.random.normal() + + vav = np.array([va, 0, 0]) + vav[0] += F * dt + np.sqrt(2 * Dpara * dt) * dW + vav[1] += np.sqrt(2 * Dperp * dt) * dW + vav[2] += np.sqrt(2 * Dperp * dt) * dW + va = np.linalg.norm(vav) + + it = np.argmin(np.abs(tv - t)) + + nfast[it] += dt + pressure[it] += ma * va * va / 3 * dt + + return nfast / Nmc, pressure / Nmc + + +def simulate_slowingdown( + ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10 +): + """ + Calculate steady-state (time independent) fast ion density and pressure. Assumes zero orbit + width and no FLR. + + ne: Background electron density profile (1/m3) - profile on rho + Te: Background electron temperature profile (eV) + mass: Background species mass (kg) - main ion + charge: Background species charge (C) - main ion + E_fast: Fast ion energy (eV) - vectors of energy components [E, E/2, E/3] + source_fast: Fast ion source profile (1/s m3) - (rho, source) source of beam components [E, E/2, E/3] + mass_fast: Fast ion mass (kg) - main ion + charge_fast: Fast ion charge (C) - main ion + + Returns: + nfast: Fast ion density profile (1/m3) - profile on same radial grid as the sources in input + pressure: Fast ion pressure profile (Pa) + """ + N = len(ne) + + nfast = np.zeros(N) + pressure = np.zeros(N) + + for i in range(N): + for j in range(3): + nfast_frac, pressure_frac = slowingdown( + mass_fast, + charge_fast, + E_fast[j] * 1.602e-19, + np.array([9.109e-31, mass]), + np.array([-1.602e-19, charge]), + ne[i], + Te[i] * 1.602e-19, + Nmc, + ) + nfast[i] += source_fast[i, j] * nfast_frac + pressure[i] += source_fast[i, j] * pressure_frac + + return nfast, pressure + + +def simulate_slowingdown_timedep( + tv, ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10 +): + """ + Calculate steady-state fast ion density and pressure. Assumes zero orbit + width and no FLR. + + ne: Background electron density profile (1/m3) + Te: Background electron temperature profile (eV) + mass: Background species mass (kg) + charge: Background species charge (C) + E_fast: Fast ion energy (eV) + source_fast: Fast ion source profile (1/s m3) + mass_fast: Fast ion mass (kg) + charge_fast: Fast ion charge (C) + + Returns: + nfast: Fast ion density profile (1/m3) + pressure: Fast ion pressure profile (Pa) + """ + N = len(ne[0, :]) + Nt = len(tv) + + dt = tv[1] - tv[0] + + nfast = np.zeros((Nt, N)) + pressure = np.zeros((Nt, N)) + + for it in range(Nt): + for i in range(N): + for j in range(3): + nfast_frac, pressure_frac = slowingdown( + mass_fast, + charge_fast, + E_fast[j] * 1.602e-19, + np.array([9.109e-31, mass]), + np.array([-1.602e-19, charge]), + ne[it, i], + Te[it, i] * 1.602e-19, + Nmc, + dt, + ) + + nfast_frac = nfast_frac[0 : Nt - it] + nfast[it : it + len(nfast_frac), i] += ( + source_fast[it, i, j] * nfast_frac + ) + + pressure_frac = pressure_frac[0 : Nt - it] + pressure[it : it + len(pressure_frac), i] += ( + source_fast[it, i, j] * pressure_frac + ) + + return nfast, pressure + + +def suzuki(E, ne, Te, Anum): + """ + Calculate beam-stopping coefficient from Suzuki et al. + + E: Neutral energy/amu (keV) + ne: Electron density (1/m3) + Te: Electron temperature (eV) + """ + A = np.array( + [ + [ + -52.9, + -1.36, + 0.0719, + 0.0137, + 0.454, + 0.403, + -0.22, + 0.0666, + -0.0677, + -0.00148, + ], + [ + -67.9, + -1.22, + 0.0814, + 0.0139, + 0.454, + 0.465, + -0.273, + 0.0751, + -0.063, + -0.000508, + ], + [ + -74.2, + -1.18, + 0.0843, + 0.0139, + 0.453, + 0.491, + -0.294, + 0.0788, + -0.0612, + -0.000185, + ], + ] + ) + + sigmav = ( + ne + * A[Anum - 1, 0] + * 1.0e-16 + / E + * (1 + A[Anum - 1, 1] * np.log(E) + A[Anum - 1, 2] * np.log(E) * np.log(E)) + * ( + 1 + + (1 - np.exp(-A[Anum - 1, 3] * ne * 1e-19)) ** A[Anum - 1, 4] + * ( + A[Anum - 1, 5] + + A[Anum - 1, 6] * np.log(E) + + A[Anum - 1, 7] * np.log(E) * np.log(E) + ) + ) + * ( + 1 + + A[Anum - 1, 8] * np.log(Te * 1e-3) + + A[Anum - 1, 9] * np.log(Te * 1e-3) * np.log(Te * 1e-3) + ) + ) / ne + + return sigmav + + +def simulate_source( + rhov, ne, Te, Anum, Rv, zv, rho2d, vol, source, direction, energy_per_anum +): + """ + rhov - profiles rho values + ne (m**-3) - profiles + Te (eV) - profiles + Anum - mass number main ion + Rv - rho2d major radius grid + zv - rho2d z grid + rho2d - equilibrium rho 2D (R, z) + vol (m**3) - equilibrium_data["vjac"] * dpsin, not cumulative sum! + source (m) - position of the source in (x, y, z) + direction (m) - position of the source in (x, y, z) + energy_per_anum (keV) - 1 energy only, call this three times for each energy + + Returns + ------- + + """ + ds = 0.001 + s = 0 + x = source + rho = evaluate_rho(x, Rv, zv, rho2d) + + # step until in plasma + while rho > 1: + x += ds * direction + rho = evaluate_rho(x, Rv, zv, rho2d) + + source_fast = np.zeros(len(ne)) + weight = 1 + + attenuation = [] + + while rho <= 1: + x += ds * direction + s += ds + rho = evaluate_rho(x, Rv, zv, rho2d) + irho = np.argmin(np.abs(rho - rhov)) + rate = ne[irho] * 1e-4 * suzuki(energy_per_anum, ne[irho], Te[irho], Anum) + source_fast[irho] += weight * rate * ds + weight -= weight * rate * ds + attenuation.append([s, weight]) + + return source_fast / vol, np.array(attenuation) + + +def evaluate_rho(x, Rv, zv, rho2d): + frho = interp.interp2d(Rv, zv, rho2d, kind="cubic") + return frho(np.sqrt(x[0] ** 2 + x[1] ** 2), x[2])[0] From 3fb87369fdd68224b05bf6b11a6260ae9ce5c13e Mon Sep 17 00:00:00 2001 From: Marco Sertoli Date: Fri, 11 Aug 2023 13:53:34 +0100 Subject: [PATCH 02/16] Added dummy workflow for slowingdown.py from Jari --- indica/workflows/run_slowingdown_dummy.py | 40 +++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 indica/workflows/run_slowingdown_dummy.py diff --git a/indica/workflows/run_slowingdown_dummy.py b/indica/workflows/run_slowingdown_dummy.py new file mode 100644 index 00000000..3594c168 --- /dev/null +++ b/indica/workflows/run_slowingdown_dummy.py @@ -0,0 +1,40 @@ +from MDSplus import Tree +import scipy.interpolate as interp +import slowingdown + +amu2kg = 1.672e-27 +eV2J = 1.602e-19 + +rhov=np.linspace(0.05,0.95,10) +Te = np.array([3935.23408201, 3788.92702285, 3521.80277339, 3184.75402347, + 2818.53606058, 2440.74809717, 2041.54215192, 1590.61672625, + 1056.6714395 , 409.09603571]) +ne = np.array([7.13449907e+19, 6.99981218e+19, 6.76047271e+19, 6.46470768e+19, + 6.14714559e+19, 5.80791621e+19, 5.38962101e+19, 4.75178926e+19, + 3.68334779e+19, 1.85729312e+19]) +anum_plasma = 2 +znum_plasma = 1 +anum_beam = 2 +znum_beam = 1 +full_energy = 55e3 +energy_frac = [full_energy, full_energy / 2, full_energy / 3] +power = 1e6 +power_frac = [0.6, 0.3, 0.1] + +# Rv, zv, rho2d, vol from equilibrium + +location_hnbi = array([0.33704, 0.93884, 0. ]) # hnbi entrance +direction_hnbi = np.array([-0.704,-0.709,0.0]) + +location_rfx = array([-0.341, -0.940, 0.0]) # rfx entrance +direction_rfx = np.array([0.704, 0.709, 0.0]) + +source = np.zeros((len(rhov), 3)) + +for i in range(3): + source[:,i] = slowingdown.simulate_finite_source(rhov, ne, Te, anum_plasma, Rv, zv, rho2d, vol, location_hnbi, direction_hnbi, energy_frac[i], anum_beam, power,width=0.25,n=10) + +out = slowingdown.simulate_slowingdown(ne, Te, anum_plasma * amu2kg, znum_plasma * eV2J, energy_frac, source, anum_beam * amu2kg, znum_beam * eV2J, Nmc=10) + +plt.plot(rhov, out["nfast"]) +plt.plot(rhov,out["pressure"]) From c5df5059a2815e34a942950559fefbac186cb26e Mon Sep 17 00:00:00 2001 From: Marco Sertoli Date: Fri, 11 Aug 2023 13:55:15 +0100 Subject: [PATCH 03/16] Updated slowingdown.py to latest version from Jari --- indica/operators/slowingdown.py | 275 +++++++++++++++----------------- 1 file changed, 127 insertions(+), 148 deletions(-) diff --git a/indica/operators/slowingdown.py b/indica/operators/slowingdown.py index 81977473..b08896f1 100644 --- a/indica/operators/slowingdown.py +++ b/indica/operators/slowingdown.py @@ -19,7 +19,7 @@ def coulomb_log(ma, qa, va, mb, qb, nb, Tb): nb: Background species density (1/m3) Tb: Background species temperature (1/m3) """ - debyeLength = np.sqrt(epsilon0 / np.sum(nb * qb * qb / Tb)) + debyeLength = np.sqrt(epsilon0/np.sum(nb*qb*qb/Tb)) vbar = va * va + 2 * Tb / mb mr = ma * mb / (ma + mb) bcl = np.abs(qa * qb / (4 * np.pi * epsilon0 * mr * vbar)) @@ -29,14 +29,13 @@ def coulomb_log(ma, qa, va, mb, qb, nb, Tb): for i in range(len(clog)): if bcl[i] > bqm[i]: - clog[i] = np.log(debyeLength / bcl[i]) + clog[i] = np.log( debyeLength / bcl[i]) else: - clog[i] = np.log(debyeLength / bqm[i]) + clog[i] = np.log( debyeLength / bqm[i]) return clog - -def slowingdown(ma, qa, Ea, mb, qb, nb, Tb, Nmc=10, dt=1): +def slowingdown(ma, qa, Ea, mb, qb, nb, Tb, Nmc=10, dt=1e-3): """ Calculate slowing-down for a fast particle. Collision operator adopted from Hirvijoki et al 2014 CPC. @@ -65,37 +64,33 @@ def slowingdown(ma, qa, Ea, mb, qb, nb, Tb, Nmc=10, dt=1): while va > np.sqrt(2 * Tb / ma): t += dt - mu0 = ( - erf(va / vb) - 2 * va / vb * np.exp(-va / vb * va / vb) / np.sqrt(np.pi) - ) / (va / vb * va / vb) - mu1 = erf(va / vb) - 0.5 * mu0 + mu0 = (erf(va/vb) - 2*va/vb*np.exp(-va/vb*va/vb) / np.sqrt(np.pi)) / (va/vb*va/vb) + mu1 = erf(va/vb) - 0.5 * mu0 clog = coulomb_log(ma, qa, va, mb, qb, nb, Tb) Cb = nb * qa * qa * qb * qb * clog / (4 * np.pi * epsilon0 * epsilon0) - F = np.sum(-(1 / mb + 1 / ma) * Cb * mu0 / (ma * vb * vb)) + F = np.sum(-(1/mb + 1/ma) * Cb * mu0 / (ma * vb * vb)) Dpara = np.sum(Cb * mu0 / (2 * ma * ma * va)) Dperp = np.sum(Cb * mu1 / (2 * ma * ma * va)) dW = np.random.normal() vav = np.array([va, 0, 0]) - vav[0] += F * dt + np.sqrt(2 * Dpara * dt) * dW - vav[1] += np.sqrt(2 * Dperp * dt) * dW - vav[2] += np.sqrt(2 * Dperp * dt) * dW + vav[0] += F*dt + np.sqrt(2*Dpara*dt) * dW + vav[1] += np.sqrt(2*Dperp*dt) * dW + vav[2] += np.sqrt(2*Dperp*dt) * dW va = np.linalg.norm(vav) - it = np.argmin(np.abs(tv - t)) + it = np.argmin(np.abs(tv-t)) nfast[it] += dt pressure[it] += ma * va * va / 3 * dt - return nfast / Nmc, pressure / Nmc + return nfast/Nmc, pressure/Nmc -def simulate_slowingdown( - ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10 -): +def simulate_slowingdown(ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10): """ - Calculate steady-state (time independent) fast ion density and pressure. Assumes zero orbit + Calculate steady-state fast ion density and pressure. Assumes zero orbit width and no FLR. ne: Background electron density profile (1/m3) - profile on rho @@ -118,25 +113,16 @@ def simulate_slowingdown( for i in range(N): for j in range(3): - nfast_frac, pressure_frac = slowingdown( - mass_fast, - charge_fast, - E_fast[j] * 1.602e-19, - np.array([9.109e-31, mass]), - np.array([-1.602e-19, charge]), - ne[i], - Te[i] * 1.602e-19, - Nmc, - ) - nfast[i] += source_fast[i, j] * nfast_frac - pressure[i] += source_fast[i, j] * pressure_frac - - return nfast, pressure - - -def simulate_slowingdown_timedep( - tv, ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10 -): + nfast_frac, pressure_frac = slowingdown(mass_fast, charge_fast, E_fast[j]*1.602e-19, np.array([9.109e-31,mass]), np.array([-1.602e-19,charge]), ne[i], Te[i]*1.602e-19,Nmc) + nfast[i] += source_fast[i,j] * np.sum(nfast_frac) + pressure[i] += source_fast[i,j] * np.sum(pressure_frac) + + out = {"nfast": nfast, "pressure": pressure} + + return out + + +def simulate_slowingdown_timedep(tv, ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10): """ Calculate steady-state fast ion density and pressure. Assumes zero orbit width and no FLR. @@ -154,10 +140,10 @@ def simulate_slowingdown_timedep( nfast: Fast ion density profile (1/m3) pressure: Fast ion pressure profile (Pa) """ - N = len(ne[0, :]) + N = len(ne[0,:]) Nt = len(tv) - dt = tv[1] - tv[0] + dt = tv[1]-tv[0] nfast = np.zeros((Nt, N)) pressure = np.zeros((Nt, N)) @@ -165,29 +151,17 @@ def simulate_slowingdown_timedep( for it in range(Nt): for i in range(N): for j in range(3): - nfast_frac, pressure_frac = slowingdown( - mass_fast, - charge_fast, - E_fast[j] * 1.602e-19, - np.array([9.109e-31, mass]), - np.array([-1.602e-19, charge]), - ne[it, i], - Te[it, i] * 1.602e-19, - Nmc, - dt, - ) - - nfast_frac = nfast_frac[0 : Nt - it] - nfast[it : it + len(nfast_frac), i] += ( - source_fast[it, i, j] * nfast_frac - ) - - pressure_frac = pressure_frac[0 : Nt - it] - pressure[it : it + len(pressure_frac), i] += ( - source_fast[it, i, j] * pressure_frac - ) - - return nfast, pressure + nfast_frac, pressure_frac = slowingdown(mass_fast, charge_fast, E_fast[j]*1.602e-19, np.array([9.109e-31,mass]), np.array([-1.602e-19,charge]), ne[it,i], Te[it,i]*1.602e-19,Nmc,dt) + + nfast_frac = nfast_frac[0:Nt-it] + nfast[it:it+len(nfast_frac),i] += source_fast[it,i,j] * nfast_frac + + pressure_frac = pressure_frac[0:Nt-it] + pressure[it:it+len(pressure_frac),i] += source_fast[it,i,j] * pressure_frac + + out = {"nfast": nfast, "pressure": pressure} + + return out def suzuki(E, ne, Te, Anum): @@ -199,91 +173,83 @@ def suzuki(E, ne, Te, Anum): Te: Electron temperature (eV) """ A = np.array( - [ - [ - -52.9, - -1.36, - 0.0719, - 0.0137, - 0.454, - 0.403, - -0.22, - 0.0666, - -0.0677, - -0.00148, - ], - [ - -67.9, - -1.22, - 0.0814, - 0.0139, - 0.454, - 0.465, - -0.273, - 0.0751, - -0.063, - -0.000508, - ], - [ - -74.2, - -1.18, - 0.0843, - 0.0139, - 0.453, - 0.491, - -0.294, - 0.0788, - -0.0612, - -0.000185, - ], - ] - ) - - sigmav = ( - ne - * A[Anum - 1, 0] - * 1.0e-16 - / E - * (1 + A[Anum - 1, 1] * np.log(E) + A[Anum - 1, 2] * np.log(E) * np.log(E)) - * ( - 1 - + (1 - np.exp(-A[Anum - 1, 3] * ne * 1e-19)) ** A[Anum - 1, 4] - * ( - A[Anum - 1, 5] - + A[Anum - 1, 6] * np.log(E) - + A[Anum - 1, 7] * np.log(E) * np.log(E) - ) - ) - * ( - 1 - + A[Anum - 1, 8] * np.log(Te * 1e-3) - + A[Anum - 1, 9] * np.log(Te * 1e-3) * np.log(Te * 1e-3) + [[-52.9, -1.36, 0.0719, 0.0137, 0.454, 0.403, -0.22, 0.0666, -0.0677, + -0.00148], + [-67.9, -1.22, 0.0814, 0.0139, 0.454, 0.465, -0.273, 0.0751, -0.063, + -0.000508], + [-74.2, -1.18, 0.0843, 0.0139, 0.453, 0.491, -0.294, 0.0788, -0.0612, + -0.000185]] ) - ) / ne + + sigmav = (ne * A[Anum-1,0] * 1.e-16 / E + * (1 + A[Anum-1,1]*np.log(E) + A[Anum-1,2]*np.log(E)*np.log(E)) + * (1 + (1 - np.exp(-A[Anum-1,3]*ne*1e-19)) ** A[Anum-1,4] + * (A[Anum-1,5] + A[Anum-1,6]*np.log(E) + A[Anum-1,7]*np.log(E)*np.log(E))) + * (1 + A[Anum-1,8]*np.log(Te*1e-3) + A[Anum-1,9]*np.log(Te*1e-3)*np.log(Te*1e-3))) / ne return sigmav -def simulate_source( - rhov, ne, Te, Anum, Rv, zv, rho2d, vol, source, direction, energy_per_anum -): +def simulate_finite_source(rhov, ne, Te, Anum, Rv, zv, rho2d, vol, source, direction, energy, Anum_fast, power, width, n=5): + """ + Calculate fast ion source from beam ionisation. + + rhov: Rho axis for profiles + ne: Background electron density profile (1/m3) + Te: Background electron temperature profile (eV) + Anum: Background species mass number + Rv: R axis for 2d rho map (m) + zv: z axis for 2d rho map (m) + rho2d: 2d (R,z) rho map + vol: differential volume profile (m3) + source: beam source location (x,y,z) (m) + direction: beam direction (dx,dy,dz) unit vector + energy: Neutral energy (eV) + Anum_fast: Beam mass number + power: Beam power (W) + width: Width of circular beam (m) + n: number of MC pencil beams + + Returns: + source_fast: Fast ion source profile (1/m3/s) + """ + source_fast = np.zeros(len(rhov)) + + dirz = np.array([0, 0, 1]) + dirh = np.cross(direction, dirz) + + for i in range(n): + source_shift = (np.random.rand()*width - width/2)*dirh + (np.random.rand()*width - width/2)*dirz + source_temp = simulate_source(rhov, ne, Te, Anum, Rv, zv, rho2d, vol, source + source_shift, direction, energy, Anum_fast, power) + source_fast += source_temp +# attenuation = attenuation_temp + + return source_fast/n + + +def simulate_source(rhov, ne, Te, Anum, Rv, zv, rho2d, vol, source, direction, energy, Anum_fast, power): """ - rhov - profiles rho values - ne (m**-3) - profiles - Te (eV) - profiles - Anum - mass number main ion - Rv - rho2d major radius grid - zv - rho2d z grid - rho2d - equilibrium rho 2D (R, z) - vol (m**3) - equilibrium_data["vjac"] * dpsin, not cumulative sum! - source (m) - position of the source in (x, y, z) - direction (m) - position of the source in (x, y, z) - energy_per_anum (keV) - 1 energy only, call this three times for each energy - - Returns - ------- + Calculate fast ion source from beam ionisation. + rhov: Rho axis for profiles + ne: Background electron density profile (1/m3) + Te: Background electron temperature profile (eV) + Anum: Background species mass number + Rv: R axis for 2d rho map (m) + zv: z axis for 2d rho map (m) + rho2d: 2d (R,z) rho map + vol: differential volume profile (m3) + source: beam source location (x,y,z) (m) + direction: beam direction (dx,dy,dz) unit vector + energy: Neutral energy (eV) + Anum_fast: Beam mass number + power: Beam power (W) + + Returns: + source_fast: Fast ion source profile (1/m3/s) """ + energy_per_anum = energy/1e3 / Anum_fast + ds = 0.001 s = 0 x = source @@ -292,26 +258,39 @@ def simulate_source( # step until in plasma while rho > 1: x += ds * direction + s += ds rho = evaluate_rho(x, Rv, zv, rho2d) source_fast = np.zeros(len(ne)) weight = 1 - attenuation = [] + attenuation=[] while rho <= 1: x += ds * direction s += ds rho = evaluate_rho(x, Rv, zv, rho2d) - irho = np.argmin(np.abs(rho - rhov)) - rate = ne[irho] * 1e-4 * suzuki(energy_per_anum, ne[irho], Te[irho], Anum) - source_fast[irho] += weight * rate * ds - weight -= weight * rate * ds - attenuation.append([s, weight]) + irho = np.argmin(np.abs(rho-rhov)) + ne_x = np.interp(rho,rhov,ne) + Te_x = np.interp(rho,rhov,Te) + rate = ne_x * 1e-4 * suzuki(energy_per_anum, ne_x, Te_x, Anum) + source_fast[irho] += weight*rate*ds + weight-=weight*rate*ds + attenuation.append([s,weight]) - return source_fast / vol, np.array(attenuation) + source_fast *= power / (energy*1.602e-19) + + return source_fast/vol def evaluate_rho(x, Rv, zv, rho2d): frho = interp.interp2d(Rv, zv, rho2d, kind="cubic") - return frho(np.sqrt(x[0] ** 2 + x[1] ** 2), x[2])[0] + return frho(np.sqrt(x[0]**2 + x[1]**2), x[2])[0] + + +def parabolic_plasma(rhov, neave, T0, Sn, ST): + ne = 0.1e19 + (3e19 - 0.1e19) * (1 - rhov**2)**Sn + ne = ne * neave / np.trapz(ne, rhov)/(rhov[-1] - rhov[0]) + Te = 100 + (T0 - 100) * (1 - rhov**2)**ST + + return ne, Te From 806f97d7ae9bff6d681b6eb196d62dbc6d5f1e52 Mon Sep 17 00:00:00 2001 From: Marco Sertoli Date: Tue, 15 Aug 2023 09:19:21 +0100 Subject: [PATCH 04/16] Fist stab at implementing NBI object that runs slowingdown.py. Still haven't tried it. Some pre-commit things to solve. List of potential TDOOs to harmonize with Indica and make workflow clearer. --- indica/models/neutral_beam.py | 428 +++++++++++++--------- indica/operators/slowingdown.py | 232 +++++++++--- indica/workflows/run_slowingdown_dummy.py | 75 +++- 3 files changed, 480 insertions(+), 255 deletions(-) diff --git a/indica/models/neutral_beam.py b/indica/models/neutral_beam.py index 0556f1e0..888bf8d8 100644 --- a/indica/models/neutral_beam.py +++ b/indica/models/neutral_beam.py @@ -1,16 +1,41 @@ -from typing import Tuple - import numpy as np +from scipy import constants +from xarray import DataArray + +from indica.converters.line_of_sight import LineOfSightTransform +from indica.datatypes import ELEMENTS +from indica.equilibrium import Equilibrium +from indica.models.abstractdiagnostic import DiagnosticModel +from indica.models.plasma import example_run as example_plasma +from indica.numpy_typing import LabeledArray +from indica.operators.slowingdown import simulate_finite_source +from indica.operators.slowingdown import simulate_slowingdown -from ..converters.line_of_sight import LineOfSightTransform +AMU2KG = constants.m_p +EV2J = constants.e +LOCATION_HNBI = np.array([[0.33704, 0.93884, 0.0]]) +DIRECTION_HNBI = np.array([[-0.704, -0.709, 0.0]]) -analytical_beam_defaults = { - "element": "H", - "amu": int(1), - "energy": 25.0 * 1e3, - "power": 500 * 1e3, - "fractions": (0.7, 0.1, 0.2, 0.0), +LOCATION_RFX = np.array([[-0.341, -0.940, 0.0]]) +DIRECTION_RFX = np.array([[0.704, 0.709, 0.0]]) + +RFX_DEFAULTS = { + "element": "d", + "energy": 25.0e3, + "power": 0.5e6, + "power_frac": (0.7, 0.1, 0.2), # TODO: rename fractions + "divergence": (14 * 1e-3, 14e-3), + "width": (0.025, 0.025), + "location": (-0.3446, -0.9387, 0.0), + "direction": (0.707, 0.707, 0.0), + "focus": 1.8, +} +HNBI_DEFAULTS = { + "element": "d", + "energy": 55.0e3, + "power": 1.0e6, + "power_frac": (0.6, 0.3, 0.1), "divergence": (14 * 1e-3, 14e-3), "width": (0.025, 0.025), "location": (-0.3446, -0.9387, 0.0), @@ -19,192 +44,231 @@ } -class NeutralBeam: - def __init__(self, name: str, use_defaults=True, **kwargs): - """Set beam parameters, initialisation""" - self.name = name # Beam name +class NeutralBeam(DiagnosticModel): + def __init__( + self, + name: str, + instrument_method: str = "get_neutral_beam", + element: str = "d", + energy: float = 55.0e3, + power: float = 1.0e6, + power_frac: LabeledArray = [0.6, 0.3, 0.1], + divergence: LabeledArray = [14 * 1e-3, 14e-3], + width: LabeledArray = [0.025, 0.025], + focus: float = 1.8, + n_beamlets: int = 10, + n_mc: int = 10, + **kwargs, + ): + """ - # Use beam defaults - if use_defaults: - for (prop, default) in analytical_beam_defaults.items(): - setattr(self, prop, kwargs.get(prop, default)) - else: - for (prop, default) in analytical_beam_defaults.items(): - setattr(self, prop, None) - return + Parameters + ---------- + name + String identifier of the system + instrument_method + Corresponding method to read system's data from ST40 database + element + Beam gas + energy + Beam energy (V) (TODO: substitute with voltage?) + power + Beam power (W) + power_frac + Fraction of power in 1st, 2nd and 3rd energy + n_beamlets + Number of beamlets + n_mc + Number of MC-like samples + TODO: substitute energy with voltage (engineering parameter) + TODO: substitute power with current? + TODO: n_beamlets to be incorporated in LOS-transform as different LOSs + --> will substitute evaluate_rho in slowingdown.py + --> rho2d interpolation performed once! + --> this can become a version of VOS-implementation + TODO: rename anum, znum -> atomic_mass, atomic_number + + """ + + self.name = name + self.instrument_method = instrument_method + self.element = element + self.anum_beam = ELEMENTS[element.lower()][1] + self.znum_beam = ELEMENTS[element.lower()][0] + self.power = power + self.power_frac = np.array(power_frac) + self.energy = energy + self.energy_frac = energy / np.arange(1, np.size(power_frac) + 1) + self.focus = focus + self.width = np.array(width) + self.divergence = np.array(divergence) + self.n_beamlets = n_beamlets + self.n_mc = n_mc + + def get_beam_parameters(self): + beam_params = { + "element": self.element, + "anum_beam": self.anum_beam, + "znum_beam": self.znum_beam, + "power": self.power, + "power_frac": self.power_frac, + "energy": self.energy, + "energy_frac": self.energy_frac, + "focus": self.focus, + "divergence": self.divergence, + "width": self.width, + } - # Set line_of_sight transform for centre of beam-line - self.set_los_transform() + return beam_params - # Set Attenuator - self.attenuator = None + def set_beam_parameters(self, parameters: dict): + for k, v in parameters: + setattr(self, k, v) - # Print statements for debugging - print(f"Beam = {self.name}") - print(f"Energy = {self.energy} electron-volts") - print(f"Power = {self.power} watts") - print(f"Fractions (full, 1/2, 1/3, imp) = {self.fractions} %") - print(f"divergence (x, y) = {self.divergence} rad") - print(f"width (x, y) = {self.width} metres") + def _build_bckc_dictionary(self): + self.bckc = { + "fast_density": self.fast_density, + "fast_pressure_parallel": self.fast_pressure_parallel, + "fast_pressure_perpendicular": self.fast_pressure_perpendicular, + } - def set_los_transform( + def __call__( self, - machine_dimensions: Tuple[Tuple[float, float], Tuple[float, float]] = ( - (0.175, 1.0), - (-2.0, 2.0), - ), + Ne: DataArray = None, + Te: DataArray = None, + equilibrium: Equilibrium = None, + main_ion: str = None, + t: LabeledArray = None, + **kwargs, ): - self.transform = LineOfSightTransform( - np.array([self.location[0]]), - np.array([self.location[1]]), - np.array([self.location[2]]), - np.array([self.direction[0]]), - np.array([self.direction[1]]), - np.array([self.direction[2]]), - name=f"{self.name}_los", - dl=0.01, - machine_dimensions=machine_dimensions, + """ + Parameters + ---------- + Ne + Electron density profile (m**-3) + Te + Electron temperature profile (eV) + equilibrium + indica.equilibrium.Equilibrium object + main_ion + String identifier of plasma main ion + t + Desired time(s) of analysis + kwargs + ... + Returns + ------- + Dictionary with fast ion densities and pressure + + TODO: rename Rv, zv, rho_v --> R_equil, z_equil, rho_equil + TODO: Indica native rho starts at 0 - check rho_profile below is self-consistent + TODO: current implementation assumes pure plasma: expand to include impurities? + TODO: fast pressure currently split equally between parallel and perpendicular + """ + + if self.plasma is not None: + if t is None: + t = self.plasma.time_to_calculate + + Ne = self.plasma.electron_density.interp(t=t, method="nearest") + Te = self.plasma.electron_temperature.interp(t=t, method="nearest") + equilibrium = self.plasma.equilibrium + main_ion = self.plasma.main_ion + else: + if Ne is None or Te is None or equilibrium is None or main_ion is None: + raise ValueError("Give inputs or assign plasma class!") + + R_equil = np.array(equilibrium.rho.R) + z_equil = np.array(equilibrium.rho.z) + rho_equil = np.array(equilibrium.rho.sel(t=t, method="nearest")) + rho_profile = np.array(Ne.rho_poloidal) + rho_profile[0] = rho_profile[1] / 2.0 + vol_profile = np.array( + equilibrium.volume.sel(t=t, method="nearest") + .diff("rho_poloidal") + .interp(rho_poloidal=rho_profile) ) - def run_BBNBI(self): - print("Add code to run BBNBI") - - # def run_analytical_beam( - # self, - # x_dash: LabeledArray = DataArray(np.linspace(-0.25, 0.25, 101, dtype=float)), - # y_dash: LabeledArray = DataArray(np.linspace(-0.25, 0.25, 101, dtype=float)), - # ): - # """Analytical beam based on double gaussian formula""" - # - # # Calculate Attenuator - # x2 = self.transform.x2 - # attenuation_factor = np.ones_like( - # x2 - # ) # Replace with Attenuation Object in future, interact with plasma - # - # # Calculate beam velocity - # # v_beam = self.beam_velocity() - # # e = 1.602 * 1e-19 - # - # # Neutral beam - # n_x = 101 - # n_y = 101 - # n_z = len(attenuation_factor) - # nb_dash = np.zeros((n_z, n_x, n_y), dtype=float) - # # for i_z in range(n_z): - # # for i_y in range(n_y): - # # y_here = y_dash[i_y] - # # exp_factor = np.exp( - # # -(x_dash**2 / self.width[0]**2) - (y_here**2 / self.width[1]**2) - # # ) - # # nb_dash[i_z, :, i_y] = \ - # # self.power * attenuation_factor[i_z] * exp_factor / \ - # # (np.pi * self.energy * e * np.prod(self.width) * v_beam) - # - # # mesh-grid of beam cross section coordinates - # z_dash = x2 - # X_dash, Y_dash, Z_dash = np.meshgrid(x_dash, y_dash, z_dash) - # R_dash = np.sqrt(X_dash**2 + Y_dash**2) - # T_dash = np.arctan2(Y_dash, X_dash) - # - # x_transform = self.transform.x - # y_transform = self.transform.y - # z_transform = self.transform.z - # r_transform = np.sqrt(self.transform.x**2 + self.transform.y**2) - # theta_transform = np.arctan2( - # y_transform[1] - y_transform[0], x_transform[1] - x_transform[0] - # ) - # phi_transform = np.arctan2( - # z_transform[1] - z_transform[0], r_transform[1] - r_transform[0] - # ) - # theta_n = theta_transform + (np.pi / 2) - # phi_n = phi_transform + (np.pi / 2) - # - # xd = np.zeros_like(X_dash) - # yd = np.zeros_like(Y_dash) - # zd = np.zeros_like(Z_dash) - # for i_z in range(len(x2)): - # delta_X_dash = R_dash[:, :, i_z] * np.cos(T_dash[:, :, i_z]) - # delta_Y_dash = R_dash[:, :, i_z] * np.sin(T_dash[:, :, i_z]) - # - # xd[:, :, i_z] = x_transform[i_z].data + delta_X_dash * np.cos( - # theta_n.data - # ) # + X_dash[:, :, i_z]*np.tan(theta_n) - # yd[:, :, i_z] = y_transform[i_z].data + delta_X_dash * np.sin( - # theta_n.data - # ) # + Y_dash[:, :, i_z]*np.tan(theta_n) - # zd[:, :, i_z] = z_transform[i_z].data + delta_Y_dash * np.sin(phi_n.data) - # - # plt.figure() - # for i_z in range(len(x2)): - # plt.plot(xd[:, :, i_z].flatten(), yd[:, :, i_z].flatten(), "r.") - # plt.plot(self.transform.x, self.transform.y, "k") - # plt.axis("equal") - # plt.show(block=True) - # - # # print(X_dash) - # # print(np.shape(X_dash)) - # # print("aa" ** 2) - # - # # # Interpolate over tok-grid - # # x_tok = DataArray(np.linspace(-1.0, 1.0, 101, dtype=float)) - # # y_tok = DataArray(np.linspace(-1.0, 1.0, 101, dtype=float)) - # # z_tok = DataArray(np.linspace(-0.5, 0.5, 51, dtype=float)) - # # points = (z_tok.data, x_tok.data, y_tok.data) - # - # # X_tok, Y_tok, Z_tok = np.meshgrid(x_tok, y_tok, z_tok) - # # X_tok = X_tok.flatten() - # # Y_tok = Y_tok.flatten() - # # Z_tok = Z_tok.flatten() - # # nb_tok = np.zeros_like(X_tok) - # # for i in range(len(X_tok)): - # # print(f'{i}out of {len(X_tok)}') - # # point = np.array([Z_tok[i], X_tok[i], Y_tok[i]]) - # # nb_tok[i] = interpn(points, nb_dash, point, method='linear') - # - # # print("aa" ** 2) - # - # if True: - # - # plt.figure() - # plt.plot(x_dash, np.sum(nb_dash[0, :, :], axis=1)) - # - # plt.figure() - # plt.contour(x_dash, y_dash, nb_dash[0, :, :], 100) - # plt.xlabel("X (m)") - # plt.ylabel("Y (m)") - # plt.title("Initial beam cross section") - # plt.show(block=True) - - def beam_velocity(self): - return 4.38 * 1e5 * np.sqrt(self.energy * 1e-3 / float(self.amu)) - - def set_energy(self, energy: float): - self.energy = energy + width = self.width.mean() + anum_plasma = ELEMENTS[main_ion][1] + znum_plasma = ELEMENTS[main_ion][0] + source = np.zeros((len(rho_profile), len(self.energy_frac))) - def set_power(self, power: float): - self.power = power + for i in range(len(self.energy_frac)): + source[:, i] = simulate_finite_source( + rho_profile, + Ne, + Te, + anum_plasma, + R_equil, + z_equil, + rho_equil, + vol_profile, + self.los_transform.origin, + self.los_transform.direction, + self.energy_frac[i], + self.anum_beam, + self.power, + width=width, + n=self.n_beamlets, + ) - def set_divergence(self, divergence: tuple): - self.divergence = divergence + result = simulate_slowingdown( + Ne, + Te, + anum_plasma * AMU2KG, + znum_plasma * EV2J, + self.energy_frac, + source, + self.anum_beam * AMU2KG, + self.znum_beam * EV2J, + Nmc=self.n_mc, + ) - def set_fractions(self, fractions: tuple): - self.fractions = fractions + self.Ne = Ne + self.Te = Te + self.equilibrium = equilibrium + self.main_ion = main_ion + self.t = t - def set_width(self, width: tuple): - self.width = width + self.fast_density = result["nfast"] + self.fast_pressure_parallel = result["pressure"] / 2.0 + self.fast_pressure_perpendicular = result["pressure"] / 2.0 - def set_element(self, element: str): - self.element = element - if self.element == "H": - self.amu = int(1) - elif self.element == "D": - self.amu = int(2) - else: - raise ValueError + self._build_bckc_dictionary() + + return self.bckc + + +def example_run( + pulse: int = None, t: float = None, n_beamlets: int = 10, n_mc: int = 10 +): + plasma = example_plasma(pulse=pulse) + if t is None: + t = plasma.t.mean() + + beam_name = "hnbi1" + beam_params = HNBI_DEFAULTS + beam_params["n_beamlets"] = n_beamlets + beam_params["n_mc"] = n_mc - def set_location(self, location: tuple): - self.location = location + origin = LOCATION_HNBI + direction = DIRECTION_RFX + los_transform = LineOfSightTransform( + origin[:, 0], + origin[:, 1], + origin[:, 2], + direction[:, 0], + direction[:, 1], + direction[:, 2], + name=beam_name, + machine_dimensions=plasma.machine_dimensions, + passes=1, + ) + los_transform.set_equilibrium(plasma.equilibrium) + model = NeutralBeam(beam_name, **beam_params) + model.set_los_transform(los_transform) + model.set_plasma(plasma) + bckc = model(t=t) - def set_direction(self, direction: tuple): - self.direction = direction + return plasma, model, bckc diff --git a/indica/operators/slowingdown.py b/indica/operators/slowingdown.py index b08896f1..b0a6078e 100644 --- a/indica/operators/slowingdown.py +++ b/indica/operators/slowingdown.py @@ -1,7 +1,7 @@ import matplotlib.pyplot as plt import numpy as np -from scipy.special import erf import scipy.interpolate as interp +from scipy.special import erf epsilon0 = 8.854e-12 hbar = 1.054e-34 @@ -19,7 +19,7 @@ def coulomb_log(ma, qa, va, mb, qb, nb, Tb): nb: Background species density (1/m3) Tb: Background species temperature (1/m3) """ - debyeLength = np.sqrt(epsilon0/np.sum(nb*qb*qb/Tb)) + debyeLength = np.sqrt(epsilon0 / np.sum(nb * qb * qb / Tb)) vbar = va * va + 2 * Tb / mb mr = ma * mb / (ma + mb) bcl = np.abs(qa * qb / (4 * np.pi * epsilon0 * mr * vbar)) @@ -29,15 +29,16 @@ def coulomb_log(ma, qa, va, mb, qb, nb, Tb): for i in range(len(clog)): if bcl[i] > bqm[i]: - clog[i] = np.log( debyeLength / bcl[i]) + clog[i] = np.log(debyeLength / bcl[i]) else: - clog[i] = np.log( debyeLength / bqm[i]) + clog[i] = np.log(debyeLength / bqm[i]) return clog + def slowingdown(ma, qa, Ea, mb, qb, nb, Tb, Nmc=10, dt=1e-3): """ - Calculate slowing-down for a fast particle. Collision operator adopted from + Calculate slowing-down for a fast particle. Collision operator adopted from Hirvijoki et al 2014 CPC. ma: Fast particle mass (kg) @@ -64,31 +65,35 @@ def slowingdown(ma, qa, Ea, mb, qb, nb, Tb, Nmc=10, dt=1e-3): while va > np.sqrt(2 * Tb / ma): t += dt - mu0 = (erf(va/vb) - 2*va/vb*np.exp(-va/vb*va/vb) / np.sqrt(np.pi)) / (va/vb*va/vb) - mu1 = erf(va/vb) - 0.5 * mu0 + mu0 = ( + erf(va / vb) - 2 * va / vb * np.exp(-va / vb * va / vb) / np.sqrt(np.pi) + ) / (va / vb * va / vb) + mu1 = erf(va / vb) - 0.5 * mu0 clog = coulomb_log(ma, qa, va, mb, qb, nb, Tb) Cb = nb * qa * qa * qb * qb * clog / (4 * np.pi * epsilon0 * epsilon0) - F = np.sum(-(1/mb + 1/ma) * Cb * mu0 / (ma * vb * vb)) + F = np.sum(-(1 / mb + 1 / ma) * Cb * mu0 / (ma * vb * vb)) Dpara = np.sum(Cb * mu0 / (2 * ma * ma * va)) Dperp = np.sum(Cb * mu1 / (2 * ma * ma * va)) dW = np.random.normal() vav = np.array([va, 0, 0]) - vav[0] += F*dt + np.sqrt(2*Dpara*dt) * dW - vav[1] += np.sqrt(2*Dperp*dt) * dW - vav[2] += np.sqrt(2*Dperp*dt) * dW + vav[0] += F * dt + np.sqrt(2 * Dpara * dt) * dW + vav[1] += np.sqrt(2 * Dperp * dt) * dW + vav[2] += np.sqrt(2 * Dperp * dt) * dW va = np.linalg.norm(vav) - it = np.argmin(np.abs(tv-t)) + it = np.argmin(np.abs(tv - t)) nfast[it] += dt pressure[it] += ma * va * va / 3 * dt - return nfast/Nmc, pressure/Nmc + return nfast / Nmc, pressure / Nmc -def simulate_slowingdown(ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10): +def simulate_slowingdown( + ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10 +): """ Calculate steady-state fast ion density and pressure. Assumes zero orbit width and no FLR. @@ -113,16 +118,27 @@ def simulate_slowingdown(ne, Te, mass, charge, E_fast, source_fast, mass_fast, c for i in range(N): for j in range(3): - nfast_frac, pressure_frac = slowingdown(mass_fast, charge_fast, E_fast[j]*1.602e-19, np.array([9.109e-31,mass]), np.array([-1.602e-19,charge]), ne[i], Te[i]*1.602e-19,Nmc) - nfast[i] += source_fast[i,j] * np.sum(nfast_frac) - pressure[i] += source_fast[i,j] * np.sum(pressure_frac) + nfast_frac, pressure_frac = slowingdown( + mass_fast, + charge_fast, + E_fast[j] * 1.602e-19, + np.array([9.109e-31, mass]), + np.array([-1.602e-19, charge]), + ne[i], + Te[i] * 1.602e-19, + Nmc, + ) + nfast[i] += source_fast[i, j] * np.sum(nfast_frac) + pressure[i] += source_fast[i, j] * np.sum(pressure_frac) out = {"nfast": nfast, "pressure": pressure} return out -def simulate_slowingdown_timedep(tv, ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10): +def simulate_slowingdown_timedep( + tv, ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10 +): """ Calculate steady-state fast ion density and pressure. Assumes zero orbit width and no FLR. @@ -140,10 +156,10 @@ def simulate_slowingdown_timedep(tv, ne, Te, mass, charge, E_fast, source_fast, nfast: Fast ion density profile (1/m3) pressure: Fast ion pressure profile (Pa) """ - N = len(ne[0,:]) + N = len(ne[0, :]) Nt = len(tv) - dt = tv[1]-tv[0] + dt = tv[1] - tv[0] nfast = np.zeros((Nt, N)) pressure = np.zeros((Nt, N)) @@ -151,13 +167,27 @@ def simulate_slowingdown_timedep(tv, ne, Te, mass, charge, E_fast, source_fast, for it in range(Nt): for i in range(N): for j in range(3): - nfast_frac, pressure_frac = slowingdown(mass_fast, charge_fast, E_fast[j]*1.602e-19, np.array([9.109e-31,mass]), np.array([-1.602e-19,charge]), ne[it,i], Te[it,i]*1.602e-19,Nmc,dt) - - nfast_frac = nfast_frac[0:Nt-it] - nfast[it:it+len(nfast_frac),i] += source_fast[it,i,j] * nfast_frac - - pressure_frac = pressure_frac[0:Nt-it] - pressure[it:it+len(pressure_frac),i] += source_fast[it,i,j] * pressure_frac + nfast_frac, pressure_frac = slowingdown( + mass_fast, + charge_fast, + E_fast[j] * 1.602e-19, + np.array([9.109e-31, mass]), + np.array([-1.602e-19, charge]), + ne[it, i], + Te[it, i] * 1.602e-19, + Nmc, + dt, + ) + + nfast_frac = nfast_frac[0 : Nt - it] + nfast[it : it + len(nfast_frac), i] += ( + source_fast[it, i, j] * nfast_frac + ) + + pressure_frac = pressure_frac[0 : Nt - it] + pressure[it : it + len(pressure_frac), i] += ( + source_fast[it, i, j] * pressure_frac + ) out = {"nfast": nfast, "pressure": pressure} @@ -173,24 +203,88 @@ def suzuki(E, ne, Te, Anum): Te: Electron temperature (eV) """ A = np.array( - [[-52.9, -1.36, 0.0719, 0.0137, 0.454, 0.403, -0.22, 0.0666, -0.0677, - -0.00148], - [-67.9, -1.22, 0.0814, 0.0139, 0.454, 0.465, -0.273, 0.0751, -0.063, - -0.000508], - [-74.2, -1.18, 0.0843, 0.0139, 0.453, 0.491, -0.294, 0.0788, -0.0612, - -0.000185]] + [ + [ + -52.9, + -1.36, + 0.0719, + 0.0137, + 0.454, + 0.403, + -0.22, + 0.0666, + -0.0677, + -0.00148, + ], + [ + -67.9, + -1.22, + 0.0814, + 0.0139, + 0.454, + 0.465, + -0.273, + 0.0751, + -0.063, + -0.000508, + ], + [ + -74.2, + -1.18, + 0.0843, + 0.0139, + 0.453, + 0.491, + -0.294, + 0.0788, + -0.0612, + -0.000185, + ], + ] + ) + + sigmav = ( + ne + * A[Anum - 1, 0] + * 1.0e-16 + / E + * (1 + A[Anum - 1, 1] * np.log(E) + A[Anum - 1, 2] * np.log(E) * np.log(E)) + * ( + 1 + + (1 - np.exp(-A[Anum - 1, 3] * ne * 1e-19)) ** A[Anum - 1, 4] + * ( + A[Anum - 1, 5] + + A[Anum - 1, 6] * np.log(E) + + A[Anum - 1, 7] * np.log(E) * np.log(E) + ) ) - - sigmav = (ne * A[Anum-1,0] * 1.e-16 / E - * (1 + A[Anum-1,1]*np.log(E) + A[Anum-1,2]*np.log(E)*np.log(E)) - * (1 + (1 - np.exp(-A[Anum-1,3]*ne*1e-19)) ** A[Anum-1,4] - * (A[Anum-1,5] + A[Anum-1,6]*np.log(E) + A[Anum-1,7]*np.log(E)*np.log(E))) - * (1 + A[Anum-1,8]*np.log(Te*1e-3) + A[Anum-1,9]*np.log(Te*1e-3)*np.log(Te*1e-3))) / ne + * ( + 1 + + A[Anum - 1, 8] * np.log(Te * 1e-3) + + A[Anum - 1, 9] * np.log(Te * 1e-3) * np.log(Te * 1e-3) + ) + ) / ne return sigmav -def simulate_finite_source(rhov, ne, Te, Anum, Rv, zv, rho2d, vol, source, direction, energy, Anum_fast, power, width, n=5): +def simulate_finite_source( + rhov, + ne, + Te, + Anum, + Rv, + zv, + rho2d, + vol, + source, + direction, + energy, + Anum_fast, + power, + width, + n=5, +): """ Calculate fast ion source from beam ionisation. @@ -219,15 +313,33 @@ def simulate_finite_source(rhov, ne, Te, Anum, Rv, zv, rho2d, vol, source, direc dirh = np.cross(direction, dirz) for i in range(n): - source_shift = (np.random.rand()*width - width/2)*dirh + (np.random.rand()*width - width/2)*dirz - source_temp = simulate_source(rhov, ne, Te, Anum, Rv, zv, rho2d, vol, source + source_shift, direction, energy, Anum_fast, power) + source_shift = (np.random.rand() * width - width / 2) * dirh + ( + np.random.rand() * width - width / 2 + ) * dirz + source_temp = simulate_source( + rhov, + ne, + Te, + Anum, + Rv, + zv, + rho2d, + vol, + source + source_shift, + direction, + energy, + Anum_fast, + power, + ) source_fast += source_temp -# attenuation = attenuation_temp + # attenuation = attenuation_temp - return source_fast/n + return source_fast / n -def simulate_source(rhov, ne, Te, Anum, Rv, zv, rho2d, vol, source, direction, energy, Anum_fast, power): +def simulate_source( + rhov, ne, Te, Anum, Rv, zv, rho2d, vol, source, direction, energy, Anum_fast, power +): """ Calculate fast ion source from beam ionisation. @@ -248,7 +360,7 @@ def simulate_source(rhov, ne, Te, Anum, Rv, zv, rho2d, vol, source, direction, e Returns: source_fast: Fast ion source profile (1/m3/s) """ - energy_per_anum = energy/1e3 / Anum_fast + energy_per_anum = energy / 1e3 / Anum_fast ds = 0.001 s = 0 @@ -264,33 +376,35 @@ def simulate_source(rhov, ne, Te, Anum, Rv, zv, rho2d, vol, source, direction, e source_fast = np.zeros(len(ne)) weight = 1 - attenuation=[] + attenuation = [] while rho <= 1: x += ds * direction s += ds + print(x, ds) + rho = evaluate_rho(x, Rv, zv, rho2d) - irho = np.argmin(np.abs(rho-rhov)) - ne_x = np.interp(rho,rhov,ne) - Te_x = np.interp(rho,rhov,Te) + irho = np.argmin(np.abs(rho - rhov)) + ne_x = np.interp(rho, rhov, ne) + Te_x = np.interp(rho, rhov, Te) rate = ne_x * 1e-4 * suzuki(energy_per_anum, ne_x, Te_x, Anum) - source_fast[irho] += weight*rate*ds - weight-=weight*rate*ds - attenuation.append([s,weight]) + source_fast[irho] += weight * rate * ds + weight -= weight * rate * ds + attenuation.append([s, weight]) - source_fast *= power / (energy*1.602e-19) + source_fast *= power / (energy * 1.602e-19) - return source_fast/vol + return source_fast / vol def evaluate_rho(x, Rv, zv, rho2d): frho = interp.interp2d(Rv, zv, rho2d, kind="cubic") - return frho(np.sqrt(x[0]**2 + x[1]**2), x[2])[0] + return frho(np.sqrt(x[0] ** 2 + x[1] ** 2), x[2])[0] def parabolic_plasma(rhov, neave, T0, Sn, ST): - ne = 0.1e19 + (3e19 - 0.1e19) * (1 - rhov**2)**Sn - ne = ne * neave / np.trapz(ne, rhov)/(rhov[-1] - rhov[0]) - Te = 100 + (T0 - 100) * (1 - rhov**2)**ST + ne = 0.1e19 + (3e19 - 0.1e19) * (1 - rhov**2) ** Sn + ne = ne * neave / np.trapz(ne, rhov) / (rhov[-1] - rhov[0]) + Te = 100 + (T0 - 100) * (1 - rhov**2) ** ST return ne, Te diff --git a/indica/workflows/run_slowingdown_dummy.py b/indica/workflows/run_slowingdown_dummy.py index 3594c168..10837b5c 100644 --- a/indica/workflows/run_slowingdown_dummy.py +++ b/indica/workflows/run_slowingdown_dummy.py @@ -1,17 +1,38 @@ -from MDSplus import Tree import scipy.interpolate as interp import slowingdown amu2kg = 1.672e-27 eV2J = 1.602e-19 -rhov=np.linspace(0.05,0.95,10) -Te = np.array([3935.23408201, 3788.92702285, 3521.80277339, 3184.75402347, - 2818.53606058, 2440.74809717, 2041.54215192, 1590.61672625, - 1056.6714395 , 409.09603571]) -ne = np.array([7.13449907e+19, 6.99981218e+19, 6.76047271e+19, 6.46470768e+19, - 6.14714559e+19, 5.80791621e+19, 5.38962101e+19, 4.75178926e+19, - 3.68334779e+19, 1.85729312e+19]) +rhov = np.linspace(0.05, 0.95, 10) +Te = np.array( + [ + 3935.23408201, + 3788.92702285, + 3521.80277339, + 3184.75402347, + 2818.53606058, + 2440.74809717, + 2041.54215192, + 1590.61672625, + 1056.6714395, + 409.09603571, + ] +) +ne = np.array( + [ + 7.13449907e19, + 6.99981218e19, + 6.76047271e19, + 6.46470768e19, + 6.14714559e19, + 5.80791621e19, + 5.38962101e19, + 4.75178926e19, + 3.68334779e19, + 1.85729312e19, + ] +) anum_plasma = 2 znum_plasma = 1 anum_beam = 2 @@ -23,18 +44,44 @@ # Rv, zv, rho2d, vol from equilibrium -location_hnbi = array([0.33704, 0.93884, 0. ]) # hnbi entrance -direction_hnbi = np.array([-0.704,-0.709,0.0]) +location_hnbi = array([0.33704, 0.93884, 0.0]) # hnbi entrance +direction_hnbi = np.array([-0.704, -0.709, 0.0]) -location_rfx = array([-0.341, -0.940, 0.0]) # rfx entrance +location_rfx = array([-0.341, -0.940, 0.0]) # rfx entrance direction_rfx = np.array([0.704, 0.709, 0.0]) source = np.zeros((len(rhov), 3)) for i in range(3): - source[:,i] = slowingdown.simulate_finite_source(rhov, ne, Te, anum_plasma, Rv, zv, rho2d, vol, location_hnbi, direction_hnbi, energy_frac[i], anum_beam, power,width=0.25,n=10) + source[:, i] = slowingdown.simulate_finite_source( + rhov, + ne, + Te, + anum_plasma, + Rv, + zv, + rho2d, + vol, + location_hnbi, + direction_hnbi, + energy_frac[i], + anum_beam, + power, + width=0.25, + n=10, + ) -out = slowingdown.simulate_slowingdown(ne, Te, anum_plasma * amu2kg, znum_plasma * eV2J, energy_frac, source, anum_beam * amu2kg, znum_beam * eV2J, Nmc=10) +out = slowingdown.simulate_slowingdown( + ne, + Te, + anum_plasma * amu2kg, + znum_plasma * eV2J, + energy_frac, + source, + anum_beam * amu2kg, + znum_beam * eV2J, + Nmc=10, +) plt.plot(rhov, out["nfast"]) -plt.plot(rhov,out["pressure"]) +plt.plot(rhov, out["pressure"]) From 954f17b3407f98c634a888acbc68c9429b1c3414 Mon Sep 17 00:00:00 2001 From: Jari Varje Date: Mon, 11 Sep 2023 12:13:41 +0100 Subject: [PATCH 05/16] Separate beam LoS calculation --- indica/operators/slowingdown.py | 114 +++++++++++++++++++++++--------- 1 file changed, 84 insertions(+), 30 deletions(-) diff --git a/indica/operators/slowingdown.py b/indica/operators/slowingdown.py index b0a6078e..f1a49139 100644 --- a/indica/operators/slowingdown.py +++ b/indica/operators/slowingdown.py @@ -307,34 +307,29 @@ def simulate_finite_source( Returns: source_fast: Fast ion source profile (1/m3/s) """ - source_fast = np.zeros(len(rhov)) dirz = np.array([0, 0, 1]) dirh = np.cross(direction, dirz) + dist_los = [] + rho_los = [] + for i in range(n): source_shift = (np.random.rand() * width - width / 2) * dirh + ( np.random.rand() * width - width / 2 ) * dirz - source_temp = simulate_source( - rhov, - ne, - Te, - Anum, - Rv, - zv, - rho2d, - vol, - source + source_shift, - direction, - energy, - Anum_fast, - power, - ) - source_fast += source_temp - # attenuation = attenuation_temp - return source_fast / n + dist_los_temp, rho_los_temp = precalc_source_los(Rv, zv, rho2d, source + source_shift, direction) + + dist_los.append(dist_los_temp) + rho_los.append(rho_los_temp) + + source_fast = np.zeros(len(ne)) + + for i in range(n): + source_fast += simulate_source_los(rhov, ne, Te, Anum, dist_los[i], rho_los[i], energy, Anum_fast, power) + + return source_fast / n / vol def simulate_source( @@ -360,7 +355,33 @@ def simulate_source( Returns: source_fast: Fast ion source profile (1/m3/s) """ - energy_per_anum = energy / 1e3 / Anum_fast + + dist_los, rho_los = precalc_source_los(Rv, zv, rho2d, source, direction) + + source_fast = simulate_source_los(rhov, ne, Te, Anum, dist_los, rho_los, energy, Anum_fast, power) + + return source_fast / vol + + +def precalc_source_los( + Rv, zv, rho2d, source, direction +): + """ + Calculate fast ion source from beam ionisation. + + Rv: R axis for 2d rho map (m) + zv: z axis for 2d rho map (m) + rho2d: 2d (R,z) rho map + vol: differential volume profile (m3) + source: beam source location (x,y,z) (m) + direction: beam direction (dx,dy,dz) unit vector + energy: Neutral energy (eV) + Anum_fast: Beam mass number + power: Beam power (W) + + Returns: + source_fast: Fast ion source profile (1/m3/s) + """ ds = 0.001 s = 0 @@ -373,28 +394,61 @@ def simulate_source( s += ds rho = evaluate_rho(x, Rv, zv, rho2d) - source_fast = np.zeros(len(ne)) - weight = 1 + rho_los = [] + dist_los = [] - attenuation = [] + s = 0 while rho <= 1: + rho_los.append(rho) + dist_los.append(s) + x += ds * direction s += ds - print(x, ds) - rho = evaluate_rho(x, Rv, zv, rho2d) - irho = np.argmin(np.abs(rho - rhov)) - ne_x = np.interp(rho, rhov, ne) - Te_x = np.interp(rho, rhov, Te) + + return dist_los, rho_los + + +def simulate_source_los( + rhov, ne, Te, Anum, dist_los, rho_los, energy, Anum_fast, power +): + """ + Calculate fast ion source from beam ionisation for a single precalculated + LoS. + + rhov: Rho axis for profiles + ne: Background electron density profile (1/m3) + Te: Background electron temperature profile (eV) + Anum: Background species mass number + dist_los: Distance axis along the LoS + rho_los: Corresponding rho values anong the LoS + energy: Neutral energy (eV) + Anum_fast: Beam mass number + power: Beam power (W) + + Returns: + source_fast: Fast ion source profile (1/m3/s) + """ + energy_per_anum = energy / 1e3 / Anum_fast + + ds = dist_los[1] - dist_los[0] + s = 0 + + source_fast = np.zeros(len(ne)) + weight = 1 + + for i in range(len(rho_los)): + ne_x = np.interp(rho_los[i], rhov, ne) + Te_x = np.interp(rho_los[i], rhov, Te) rate = ne_x * 1e-4 * suzuki(energy_per_anum, ne_x, Te_x, Anum) + irho = np.argmin(np.abs(rho_los[i] - rhov)) source_fast[irho] += weight * rate * ds weight -= weight * rate * ds - attenuation.append([s, weight]) source_fast *= power / (energy * 1.602e-19) - return source_fast / vol + return source_fast def evaluate_rho(x, Rv, zv, rho2d): From dc78ebd0f30727866e2f653a340741dda7ee3866 Mon Sep 17 00:00:00 2001 From: Jari Varje Date: Tue, 12 Sep 2023 23:28:13 +0100 Subject: [PATCH 06/16] Add option to precalculate slowingdown --- indica/operators/slowingdown.py | 91 ++++++++++++++++++++++++++++++++- 1 file changed, 90 insertions(+), 1 deletion(-) diff --git a/indica/operators/slowingdown.py b/indica/operators/slowingdown.py index f1a49139..1acafcdb 100644 --- a/indica/operators/slowingdown.py +++ b/indica/operators/slowingdown.py @@ -6,6 +6,8 @@ epsilon0 = 8.854e-12 hbar = 1.054e-34 +max_tslow = 0.1 + def coulomb_log(ma, qa, va, mb, qb, nb, Tb): """ @@ -50,7 +52,6 @@ def slowingdown(ma, qa, Ea, mb, qb, nb, Tb, Nmc=10, dt=1e-3): Tb: Background species temperature (eV) dt: Timestep for output (s) """ - max_tslow = 0.1 tv = np.arange(0, max_tslow, dt) nfast = np.zeros(len(tv)) @@ -136,6 +137,31 @@ def simulate_slowingdown( return out +def simulate_slowingdown_precalc( + ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, precalc_data +): + N = len(ne) + + nfast = np.zeros(N) + pressure = np.zeros(N) + + for i in range(N): + for j in range(3): + nfast_frac = np.zeros(precalc_data["tv"].shape) + pressure_frac = np.zeros(precalc_data["tv"].shape) + + for k in range(len(nfast_frac)): + nfast_frac[k] = precalc_data["nfast"]([E_fast[j], ne[i], Te[i], precalc_data["tv"][k]]) + pressure_frac[k] = precalc_data["pressure"]([E_fast[j], ne[i], Te[i], precalc_data["tv"][k]]) + + nfast[i] += source_fast[i, j] * np.sum(nfast_frac) + pressure[i] += source_fast[i, j] * np.sum(pressure_frac) + + out = {"nfast": nfast, "pressure": pressure} + + return out + + def simulate_slowingdown_timedep( tv, ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10 ): @@ -194,6 +220,39 @@ def simulate_slowingdown_timedep( return out +def simulate_slowingdown_timedep_precalc( + tv, ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, precalc_data): + N = len(ne[0, :]) + Nt = len(tv) + + dt = tv[1] - tv[0] + + nfast = np.zeros((Nt, N)) + pressure = np.zeros((Nt, N)) + + for it in range(Nt): + for i in range(N): + for j in range(3): + nfast_frac = np.zeros(precalc_data["tv"].shape) + pressure_frac = np.zeros(precalc_data["tv"].shape) + for k in range(len(nfast_frac)): + nfast_frac[k] = precalc_data["nfast"]([E_fast[j], ne[i], Te[i], precalc_data["tv"][k]]) + pressure_frac[k] = precalc_data["pressure"]([E_fast[j], ne[i], Te[i], precalc_data["tv"][k]]) + + nfast_frac = nfast_frac[0 : Nt - it] + nfast[it : it + len(nfast_frac), i] += ( + source_fast[it, i, j] * nfast_frac + ) + + pressure_frac = pressure_frac[0 : Nt - it] + pressure[it : it + len(pressure_frac), i] += ( + source_fast[it, i, j] * pressure_frac + ) + + out = {"nfast": nfast, "pressure": pressure} + + return out + def suzuki(E, ne, Te, Anum): """ Calculate beam-stopping coefficient from Suzuki et al. @@ -462,3 +521,33 @@ def parabolic_plasma(rhov, neave, T0, Sn, ST): Te = 100 + (T0 - 100) * (1 - rhov**2) ** ST return ne, Te + + +def precalc_slowingdown(anum_plasma, anum_beam, Ev, nev, Tev,Nmc=10,dt=1e-3): + tv = np.arange(0, max_tslow, dt) + + out = dict() + out["Ev"] = Ev + out["nev"] = nev + out["Tev"] = Tev + out["tv"] = tv + nfast = np.zeros((len(Ev), len(nev), len(Tev), len(tv))) + pressure = np.zeros((len(Ev), len(nev), len(Tev), len(tv))) + + for k in range(len(Ev)): + for l in range(len(nev)): + for m in range(len(Tev)): + nfast_temp, pressure_temp = slowingdown(anum_beam * 1.661e-27, 1.602e-19, Ev[k]*1.602e-19, np.array([9.109e-31, anum_plasma * 1.661e-27]), np.array([-1.602e-19,1.602e-19]), nev[l], Tev[m]*1.602e-19, Nmc=Nmc, dt=dt) + nfast[k,l,m,:] = nfast_temp + pressure[k,l,m,:] = pressure_temp + + out["nfast"] = interp.RegularGridInterpolator( + (Ev, nev, Tev, tv), + nfast + ) + + out["pressure"] = interp.RegularGridInterpolator( + (Ev, nev, Tev, tv), + pressure + ) + return out From 8b63cc07b1b89f815557508567963dabbf5ef50c Mon Sep 17 00:00:00 2001 From: Jari Varje Date: Thu, 21 Sep 2023 10:50:05 +0100 Subject: [PATCH 07/16] Neutral beam model now running --- indica/models/neutral_beam.py | 32 +++-- indica/operators/slowingdown.py | 238 ++++++++++++++++---------------- 2 files changed, 139 insertions(+), 131 deletions(-) diff --git a/indica/models/neutral_beam.py b/indica/models/neutral_beam.py index 888bf8d8..599eddbc 100644 --- a/indica/models/neutral_beam.py +++ b/indica/models/neutral_beam.py @@ -8,8 +8,10 @@ from indica.models.abstractdiagnostic import DiagnosticModel from indica.models.plasma import example_run as example_plasma from indica.numpy_typing import LabeledArray +from indica.operators.slowingdown import precalc_finite_source_los from indica.operators.slowingdown import simulate_finite_source from indica.operators.slowingdown import simulate_slowingdown +from indica.operators.slowingdown import load_precalc_data AMU2KG = constants.m_p EV2J = constants.e @@ -37,7 +39,7 @@ "power": 1.0e6, "power_frac": (0.6, 0.3, 0.1), "divergence": (14 * 1e-3, 14e-3), - "width": (0.025, 0.025), + "width": (0.075, 0.075), "location": (-0.3446, -0.9387, 0.0), "direction": (0.707, 0.707, 0.0), "focus": 1.8, @@ -194,35 +196,39 @@ def __call__( znum_plasma = ELEMENTS[main_ion][0] source = np.zeros((len(rho_profile), len(self.energy_frac))) + precalc_los = precalc_finite_source_los( + R_equil, + z_equil, + rho_equil, + self.los_transform.origin[0], + self.los_transform.direction[0], + width=width, + n=self.n_beamlets + ) + for i in range(len(self.energy_frac)): source[:, i] = simulate_finite_source( rho_profile, Ne, Te, anum_plasma, - R_equil, - z_equil, - rho_equil, vol_profile, - self.los_transform.origin, - self.los_transform.direction, + precalc_los, self.energy_frac[i], self.anum_beam, - self.power, - width=width, - n=self.n_beamlets, + self.power_frac[i] * self.power, ) result = simulate_slowingdown( - Ne, - Te, + Ne.values, + Te.values, anum_plasma * AMU2KG, znum_plasma * EV2J, self.energy_frac, source, self.anum_beam * AMU2KG, self.znum_beam * EV2J, - Nmc=self.n_mc, + Nmc=self.n_mc ) self.Ne = Ne @@ -253,7 +259,7 @@ def example_run( beam_params["n_mc"] = n_mc origin = LOCATION_HNBI - direction = DIRECTION_RFX + direction = DIRECTION_HNBI los_transform = LineOfSightTransform( origin[:, 0], origin[:, 1], diff --git a/indica/operators/slowingdown.py b/indica/operators/slowingdown.py index 1acafcdb..6e20fb72 100644 --- a/indica/operators/slowingdown.py +++ b/indica/operators/slowingdown.py @@ -93,7 +93,8 @@ def slowingdown(ma, qa, Ea, mb, qb, nb, Tb, Nmc=10, dt=1e-3): def simulate_slowingdown( - ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10 + ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10, + precalc_data = None ): """ Calculate steady-state fast ion density and pressure. Assumes zero orbit @@ -119,40 +120,25 @@ def simulate_slowingdown( for i in range(N): for j in range(3): - nfast_frac, pressure_frac = slowingdown( - mass_fast, - charge_fast, - E_fast[j] * 1.602e-19, - np.array([9.109e-31, mass]), - np.array([-1.602e-19, charge]), - ne[i], - Te[i] * 1.602e-19, - Nmc, - ) - nfast[i] += source_fast[i, j] * np.sum(nfast_frac) - pressure[i] += source_fast[i, j] * np.sum(pressure_frac) - - out = {"nfast": nfast, "pressure": pressure} - - return out - - -def simulate_slowingdown_precalc( - ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, precalc_data -): - N = len(ne) + if precalc_data is not None: + nfast_frac = np.zeros(precalc_data["tv"].shape) + pressure_frac = np.zeros(precalc_data["tv"].shape) - nfast = np.zeros(N) - pressure = np.zeros(N) + for k in range(len(nfast_frac)): + nfast_frac[k] = precalc_data["nfast_interp"]([E_fast[j], ne[i], Te[i], precalc_data["tv"][k]]) + pressure_frac[k] = precalc_data["pressure_interp"]([E_fast[j], ne[i], Te[i], precalc_data["tv"][k]]) - for i in range(N): - for j in range(3): - nfast_frac = np.zeros(precalc_data["tv"].shape) - pressure_frac = np.zeros(precalc_data["tv"].shape) - - for k in range(len(nfast_frac)): - nfast_frac[k] = precalc_data["nfast"]([E_fast[j], ne[i], Te[i], precalc_data["tv"][k]]) - pressure_frac[k] = precalc_data["pressure"]([E_fast[j], ne[i], Te[i], precalc_data["tv"][k]]) + else: + nfast_frac, pressure_frac = slowingdown( + mass_fast, + charge_fast, + E_fast[j] * 1.602e-19, + np.array([9.109e-31, mass]), + np.array([-1.602e-19, charge]), + ne[i], + Te[i] * 1.602e-19, + Nmc, + ) nfast[i] += source_fast[i, j] * np.sum(nfast_frac) pressure[i] += source_fast[i, j] * np.sum(pressure_frac) @@ -163,7 +149,7 @@ def simulate_slowingdown_precalc( def simulate_slowingdown_timedep( - tv, ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10 + tv, ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, Nmc=10, precalc_data=None ): """ Calculate steady-state fast ion density and pressure. Assumes zero orbit @@ -193,17 +179,26 @@ def simulate_slowingdown_timedep( for it in range(Nt): for i in range(N): for j in range(3): - nfast_frac, pressure_frac = slowingdown( - mass_fast, - charge_fast, - E_fast[j] * 1.602e-19, - np.array([9.109e-31, mass]), - np.array([-1.602e-19, charge]), - ne[it, i], - Te[it, i] * 1.602e-19, - Nmc, - dt, - ) + if precalc_data is not None: + nfast_frac = np.zeros(precalc_data["tv"].shape) + pressure_frac = np.zeros(precalc_data["tv"].shape) + + for k in range(len(nfast_frac)): + nfast_frac[k] = precalc_data["nfast"]([E_fast[j], ne[i], Te[i], precalc_data["tv"][k]]) + pressure_frac[k] = precalc_data["pressure"]([E_fast[j], ne[i], Te[i], precalc_data["tv"][k]]) + + else: + nfast_frac, pressure_frac = slowingdown( + mass_fast, + charge_fast, + E_fast[j] * 1.602e-19, + np.array([9.109e-31, mass]), + np.array([-1.602e-19, charge]), + ne[it, i], + Te[it, i] * 1.602e-19, + Nmc, + dt, + ) nfast_frac = nfast_frac[0 : Nt - it] nfast[it : it + len(nfast_frac), i] += ( @@ -220,39 +215,6 @@ def simulate_slowingdown_timedep( return out -def simulate_slowingdown_timedep_precalc( - tv, ne, Te, mass, charge, E_fast, source_fast, mass_fast, charge_fast, precalc_data): - N = len(ne[0, :]) - Nt = len(tv) - - dt = tv[1] - tv[0] - - nfast = np.zeros((Nt, N)) - pressure = np.zeros((Nt, N)) - - for it in range(Nt): - for i in range(N): - for j in range(3): - nfast_frac = np.zeros(precalc_data["tv"].shape) - pressure_frac = np.zeros(precalc_data["tv"].shape) - for k in range(len(nfast_frac)): - nfast_frac[k] = precalc_data["nfast"]([E_fast[j], ne[i], Te[i], precalc_data["tv"][k]]) - pressure_frac[k] = precalc_data["pressure"]([E_fast[j], ne[i], Te[i], precalc_data["tv"][k]]) - - nfast_frac = nfast_frac[0 : Nt - it] - nfast[it : it + len(nfast_frac), i] += ( - source_fast[it, i, j] * nfast_frac - ) - - pressure_frac = pressure_frac[0 : Nt - it] - pressure[it : it + len(pressure_frac), i] += ( - source_fast[it, i, j] * pressure_frac - ) - - out = {"nfast": nfast, "pressure": pressure} - - return out - def suzuki(E, ne, Te, Anum): """ Calculate beam-stopping coefficient from Suzuki et al. @@ -332,17 +294,11 @@ def simulate_finite_source( ne, Te, Anum, - Rv, - zv, - rho2d, vol, - source, - direction, + precalc_los, energy, Anum_fast, - power, - width, - n=5, + power ): """ Calculate fast ion source from beam ionisation. @@ -367,28 +323,14 @@ def simulate_finite_source( source_fast: Fast ion source profile (1/m3/s) """ - dirz = np.array([0, 0, 1]) - dirh = np.cross(direction, dirz) - - dist_los = [] - rho_los = [] - - for i in range(n): - source_shift = (np.random.rand() * width - width / 2) * dirh + ( - np.random.rand() * width - width / 2 - ) * dirz - - dist_los_temp, rho_los_temp = precalc_source_los(Rv, zv, rho2d, source + source_shift, direction) - - dist_los.append(dist_los_temp) - rho_los.append(rho_los_temp) - source_fast = np.zeros(len(ne)) - for i in range(n): - source_fast += simulate_source_los(rhov, ne, Te, Anum, dist_los[i], rho_los[i], energy, Anum_fast, power) + n_los = len(precalc_los) + + for i_los in range(n_los): + source_fast += simulate_source_los(rhov, ne, Te, Anum, precalc_los[i_los], energy, Anum_fast, power) - return source_fast / n / vol + return source_fast / n_los / vol def simulate_source( @@ -415,13 +357,39 @@ def simulate_source( source_fast: Fast ion source profile (1/m3/s) """ - dist_los, rho_los = precalc_source_los(Rv, zv, rho2d, source, direction) + precalc_los = precalc_source_los(Rv, zv, rho2d, source, direction) - source_fast = simulate_source_los(rhov, ne, Te, Anum, dist_los, rho_los, energy, Anum_fast, power) + source_fast = simulate_source_los(rhov, ne, Te, Anum, precalc_los, energy, Anum_fast, power) return source_fast / vol +def precalc_finite_source_los( + Rv, + zv, + rho2d, + source, + direction, + width, + n=5 +): + dirz = np.array([0, 0, 1]) + dirh = np.cross(direction, dirz) + + precalc_los = [] + + for i in range(n): + source_shift = (np.random.rand() * width - width / 2) * dirh + ( + np.random.rand() * width - width / 2 + ) * dirz + + precalc_los_temp = precalc_source_los(Rv, zv, rho2d, source + source_shift, direction) + + precalc_los.append(precalc_los_temp) + + return precalc_los + + def precalc_source_los( Rv, zv, rho2d, source, direction ): @@ -466,11 +434,11 @@ def precalc_source_los( s += ds rho = evaluate_rho(x, Rv, zv, rho2d) - return dist_los, rho_los + return {"rho": rho_los, "dist": dist_los} def simulate_source_los( - rhov, ne, Te, Anum, dist_los, rho_los, energy, Anum_fast, power + rhov, ne, Te, Anum, precalc_los, energy, Anum_fast, power ): """ Calculate fast ion source from beam ionisation for a single precalculated @@ -491,17 +459,17 @@ def simulate_source_los( """ energy_per_anum = energy / 1e3 / Anum_fast - ds = dist_los[1] - dist_los[0] + ds = precalc_los["dist"][1] - precalc_los["dist"][0] s = 0 source_fast = np.zeros(len(ne)) weight = 1 - for i in range(len(rho_los)): - ne_x = np.interp(rho_los[i], rhov, ne) - Te_x = np.interp(rho_los[i], rhov, Te) + for i in range(len(precalc_los["rho"])): + ne_x = np.interp(precalc_los["rho"][i], rhov, ne) + Te_x = np.interp(precalc_los["rho"][i], rhov, Te) rate = ne_x * 1e-4 * suzuki(energy_per_anum, ne_x, Te_x, Anum) - irho = np.argmin(np.abs(rho_los[i] - rhov)) + irho = np.argmin(np.abs(precalc_los["rho"][i] - rhov)) source_fast[irho] += weight * rate * ds weight -= weight * rate * ds @@ -541,13 +509,47 @@ def precalc_slowingdown(anum_plasma, anum_beam, Ev, nev, Tev,Nmc=10,dt=1e-3): nfast[k,l,m,:] = nfast_temp pressure[k,l,m,:] = pressure_temp - out["nfast"] = interp.RegularGridInterpolator( + out["nfast"] = nfast + out["nfast_interp"] = interp.RegularGridInterpolator( (Ev, nev, Tev, tv), - nfast + nfast, + bounds_error = False, + fill_value = None ) - out["pressure"] = interp.RegularGridInterpolator( + out["pressure"] = pressure + out["pressure_interp"] = interp.RegularGridInterpolator( (Ev, nev, Tev, tv), - pressure + pressure, + bounds_error = False, + fill_value = None ) return out + +def save_precalc_data(fn, precalc_data): + out = dict() + out["Ev"] = precalc_data["Ev"] + out["nev"] = precalc_data["nev"] + out["Tev"] = precalc_data["Tev"] + out["tv"] = precalc_data["tv"] + out["nfast"] = precalc_data["nfast"] + out["pressure"] = precalc_data["pressure"] + np.save(fn, out) + +def load_precalc_data(fn): + data = np.load(fn, allow_pickle=True).item() + + data["pressure_interp"] = interp.RegularGridInterpolator( + (data["Ev"], data["nev"], data["Tev"], data["tv"]), + data["pressure"], + bounds_error = False, + fill_value = None + ) + data["nfast_interp"] = interp.RegularGridInterpolator( + (data["Ev"], data["nev"], data["Tev"], data["tv"]), + data["nfast"], + bounds_error = False, + fill_value = None + ) + + return data From 1ed82e90767ff874b5e4acd60f139249939afd6e Mon Sep 17 00:00:00 2001 From: Marco Sertoli Date: Thu, 21 Sep 2023 11:56:33 +0100 Subject: [PATCH 08/16] First commit new dev branch --- indica/converters/time.py | 2 +- indica/readers/read_st40.py | 10 +- indica/workflows/load_modelling_plasma.py | 122 ++++++++++++++++++---- 3 files changed, 108 insertions(+), 26 deletions(-) diff --git a/indica/converters/time.py b/indica/converters/time.py index cd4ced5f..7c35bd72 100644 --- a/indica/converters/time.py +++ b/indica/converters/time.py @@ -329,7 +329,7 @@ def bin_in_time_dt(tstart: float, tend: float, dt: float, data: DataArray) -> Da ------- : Array like the input, but binned along the time axis. - + TODO: add possibility of doing 50% overlap of time bins! """ check_bounds_bin(tstart, tend, dt, data) tlabels = get_tlabels_dt(tstart, tend, dt) diff --git a/indica/readers/read_st40.py b/indica/readers/read_st40.py index 34502099..29a3b10d 100644 --- a/indica/readers/read_st40.py +++ b/indica/readers/read_st40.py @@ -301,7 +301,7 @@ def plot_profile( def __call__( self, instruments: list = [], - revisions: list = None, + revisions: dict = None, map_raw: bool = False, tstart: float = None, tend: float = None, @@ -317,7 +317,7 @@ def __call__( if len(instruments) == 0: instruments = INSTRUMENTS if revisions is None: - revisions = [0] * len(instruments) + revisions = {instrument: 0 for instrument in instruments} if tstart is None: tstart = self.tstart if tend is None: @@ -327,13 +327,13 @@ def __call__( self.reset_data() self.get_equilibrium(R_shift=R_shift) - for i, instrument in enumerate(instruments): + for instrument in instruments: print(f"Reading {instrument}") if debug: - self.get_raw_data("", instrument, revisions[i]) + self.get_raw_data("", instrument, revisions[instrument]) else: try: - self.get_raw_data("", instrument, revisions[i]) + self.get_raw_data("", instrument, revisions[instrument]) except Exception as e: print(f"Error reading {instrument}: {e}") diff --git a/indica/workflows/load_modelling_plasma.py b/indica/workflows/load_modelling_plasma.py index 957577bc..38d06ecb 100644 --- a/indica/workflows/load_modelling_plasma.py +++ b/indica/workflows/load_modelling_plasma.py @@ -11,6 +11,7 @@ from indica.converters.line_of_sight import LineOfSightTransform from indica.converters.transect import TransectCoordinates from indica.equilibrium import Equilibrium +from indica.models.bolometer_camera import Bolometer from indica.models.charge_exchange import ChargeExchange from indica.models.diode_filters import BremsstrahlungDiode from indica.models.helike_spectroscopy import HelikeSpectrometer @@ -18,6 +19,7 @@ from indica.models.plasma import Plasma from indica.models.sxr_camera import SXRcamera from indica.models.thomson_scattering import ThomsonScattering +from indica.numpy_typing import RevisionLike from indica.readers.read_st40 import ReadST40 from indica.utilities import save_figure from indica.utilities import set_axis_sci @@ -37,10 +39,22 @@ "cxff_tws_c": ChargeExchange, "cxqf_tws_c": ChargeExchange, "brems": BremsstrahlungDiode, - "sxr_diode_1": SXRcamera, - "sxr_camera_4": SXRcamera, + "sxrc_xy1": Bolometer, + "sxrc_xy2": SXRcamera, + "sxr_spd": SXRcamera, + "blom_xy1": Bolometer, } -INSTRUMENTS = ["smmh1", "nirh1", "xrcs", "sxr_diode_1", "efit", "brems"] +INSTRUMENTS: list = [ + "smmh1", + "nirh1", + "xrcs", + "blom_xy1", + "sxrc_xy1", + "sxrc_xy2", + "efit", + "sxr_spd", +] +REVISIONS: dict = {instr: 0 for instr in INSTRUMENTS} FIG_PATH = f"/home/{getpass.getuser()}/figures/Indica/load_modelling_examples/" plt.ion() @@ -79,6 +93,7 @@ def plasma_code( n_rad = len(data[runs[0]]["ne"].rho_poloidal) main_ion = "h" impurities = ("ar", "c", "he") + impurities = ("c", "ar", "he") impurity_concentration = (0.001, 0.03, 0.01) _plasma = Plasma( tstart=tstart, @@ -112,11 +127,24 @@ def plasma_code( for element in _plasma.elements: _plasma.ion_temperature.loc[dict(element=element)] = Ti.values for i, impurity in enumerate(_plasma.impurities): - Nimp = ( - _data[f"niz{i+1}"].interp(rho_poloidal=_plasma.rho, t=_plasma.t) - * 1.0e19 - ) - _plasma.impurity_density.loc[dict(element=impurity)] = Nimp.values + if impurity == "c": + Nimp = ( + _data[f"niz{i+1}"].interp(rho_poloidal=_plasma.rho, t=_plasma.t) + * 1.0e19 + ) + _plasma.impurity_density.loc[dict(element=impurity)] = Nimp.values + elif impurity == "ar": + Nimp = ( + _data[f"niz{i}"].interp(rho_poloidal=_plasma.rho, t=_plasma.t) + * 1.0e19 + ) / 100 + _plasma.impurity_density.loc[dict(element=impurity)] = Nimp.values + else: + Nimp = ( + _data[f"niz{i+1}"].interp(rho_poloidal=_plasma.rho, t=_plasma.t) + * 1.0e19 + ) + _plasma.impurity_density.loc[dict(element=impurity)] = Nimp.values Nf = _data["nf"].interp(rho_poloidal=_plasma.rho, t=_plasma.t) * 1.0e19 _plasma.fast_density.values = Nf.values @@ -216,6 +244,7 @@ def example_run( pulse_code, pulse, equil, + equil_run, code, runs, comment, @@ -237,7 +266,8 @@ def example_run( if pulse is not None: print("Reading ST40 data") st40 = ReadST40(pulse, tstart, tend, dt=dt, tree="st40") - st40(instruments=INSTRUMENTS, map_diagnostics=False) + REVISIONS["efit"] = equil_run + st40(instruments=INSTRUMENTS, map_diagnostics=False, revisions=REVISIONS) if equil != code: equilibrium = {run: Equilibrium(st40.raw_data[equil]) for run in runs} @@ -283,7 +313,10 @@ def example_run( models[instrument].set_plasma(_plasma) - _bckc = models[instrument]() + if instrument == "xrcs": + _bckc = models[instrument](moment_analysis=True) + else: + _bckc = models[instrument]() bckc[run][instrument] = deepcopy(_bckc) if plot or save_fig: @@ -501,11 +534,13 @@ def plot_data_bckc_comparison( run_str += f"-{runs[1]}" fig_path = f"{FIG_PATH}{pulse_code}_{time}_{code}_{run_str}/" - norm = {} - norm["brems"] = True - norm["sxr_camera_4"] = True - norm["sxr_diode_1"] = True - norm["xrcs"] = True + norm: dict = {} + norm["xrcs"] = {} + norm["xrcs"]["spectra"] = True + # norm["brems"] = True + # norm["sxr_camera_4"] = True + norm["sxrc_xy2"] = {} + norm["sxrc_xy2"]["brightness"] = True y0 = {} y0["nirh1"] = True y0["smmh1"] = True @@ -590,10 +625,10 @@ def plot_data_bckc_comparison( _bckc = bckc[run][instrument][quantity].sel(t=tslice_binned) if instrument in norm.keys(): - _bckc = _bckc * _binned.max() / _bckc.max() - - if label is not None: - label += " (scaled)" + if quantity in norm[instrument].keys(): + _bckc = _bckc / _bckc.max() * _binned.max() + if label is not None: + label += " (scaled)" (_bckc).plot( label=label, @@ -626,6 +661,7 @@ def example_params(example: str, all_runs: bool = False): tend: float tplot: float runs = [f"RUN{run}" for run in (500 + np.arange(61, 77))] + equil_run: RevisionLike = 0 if example == "predictive": comment = "Tests using fixed-boundary predictive ASTRA" @@ -682,10 +718,56 @@ def example_params(example: str, all_runs: bool = False): tstart = 0.03 tend = 0.11 tplot = 0.1 + elif example == "aleksei_11228": + comment = "ASTRA using TS and invented Ti shapes" + pulse = 11228 + pulse_code = 13011228 + equil = "efit" # "astra" + equil_run = 1 + code = "astra" + if not all_runs: + runs = ["RUN610", "RUN611", "RUN612"] # C and Ar + runs = ["RUN623"] # C and Ar + tstart = 0.03 + tend = 0.11 + tplot = 0.08 + elif example == "alsu_11312": + comment = "ASTRA using TS and peaked Ti scaled to CXRS" + pulse = 11312 + pulse_code = 33011312 + equil = "astra" + code = "astra" + if not all_runs: + runs = ["RUN21"] + tstart = 0.065 + tend = 0.095 + tplot = 0.075 + elif example == "alsu_11314": + comment = "ASTRA using TS and peaked Ti scaled to CXRS" + pulse = 11314 + pulse_code = 33011314 + equil = "astra" + code = "astra" + if not all_runs: + runs = ["RUN12"] + tstart = 0.065 + tend = 0.095 + tplot = 0.075 + elif example == "alsu_11317": + comment = "ASTRA using TS and peaked Ti scaled to CXRS" + pulse = 11317 + pulse_code = 33011317 + equil = "astra" + code = "astra" + if not all_runs: + runs = ["RUN9"] + tstart = 0.065 + tend = 0.095 + tplot = 0.075 else: raise ValueError(f"No parameters for example {example}") - return pulse_code, pulse, equil, code, runs, comment, tstart, tend, tplot + return pulse_code, pulse, equil, equil_run, code, runs, comment, tstart, tend, tplot # def add_gacode_data( From 6af688e59735175095748e5cb61a006d4a33db14 Mon Sep 17 00:00:00 2001 From: Marco Sertoli Date: Mon, 2 Oct 2023 22:13:58 +0100 Subject: [PATCH 09/16] Fixed bug to read desired revision in ReadST40 wrapper that was affecting the assigned equilibrium. --- indica/readers/read_st40.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/indica/readers/read_st40.py b/indica/readers/read_st40.py index 29a3b10d..bcd2c577 100644 --- a/indica/readers/read_st40.py +++ b/indica/readers/read_st40.py @@ -318,6 +318,12 @@ def __call__( instruments = INSTRUMENTS if revisions is None: revisions = {instrument: 0 for instrument in instruments} + for instr in instruments: + if instr not in revisions.keys(): + revisions[instr] = 0 + if "efit" not in revisions: + revisions["efit"] = 0 + if tstart is None: tstart = self.tstart if tend is None: @@ -326,7 +332,7 @@ def __call__( dt = self.dt self.reset_data() - self.get_equilibrium(R_shift=R_shift) + self.get_equilibrium(R_shift=R_shift, revision=revisions["efit"]) for instrument in instruments: print(f"Reading {instrument}") if debug: From 7551fd32a327675795f56561d5d4731e19c6981e Mon Sep 17 00:00:00 2001 From: Marco Sertoli Date: Mon, 9 Oct 2023 09:07:20 +0100 Subject: [PATCH 10/16] Modified TS limits to physical range. Corrected bug in equilibrium reader time settings. --- indica/readers/read_st40.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/indica/readers/read_st40.py b/indica/readers/read_st40.py index bcd2c577..4589dd00 100644 --- a/indica/readers/read_st40.py +++ b/indica/readers/read_st40.py @@ -45,7 +45,7 @@ "sxrc_xy1": {"brightness": (0, np.inf)}, "sxrc_xy2": {"brightness": (0, np.inf)}, "blom_xy1": {"brightness": (0, np.inf)}, - "ts": {"te": (0, 10.0e3), "ne": (0, 1.0e21)}, + "ts": {"te": (0, np.inf), "ne": (0, np.inf)}, "pi": {"spectra": (0, np.inf)}, "tws_c": {"spectra": (0, np.inf)}, } @@ -91,8 +91,12 @@ def __init__( self.tend = tend self.dt = dt - self.reader = ST40Reader(pulse, tstart - 0.05, tend + 0.05, tree=tree) - self.reader_equil = ST40Reader(pulse, tstart - 0.1, tend + 0.1, tree=tree) + _tend = tend + dt * 2 + _tstart = tstart - dt * 2 + if _tstart < 0: + _tstart = 0.0 + self.reader = ST40Reader(pulse, _tstart, _tend, tree=tree) + self.reader_equil = ST40Reader(pulse, _tstart, _tend, tree=tree) self.equilibrium: Equilibrium self.raw_data: dict = {} From 65147e966e4d394439084cebf2ca21069f11a3a0 Mon Sep 17 00:00:00 2001 From: Marco Sertoli Date: Mon, 9 Oct 2023 09:07:53 +0100 Subject: [PATCH 11/16] Still having issues with TS MDS+ inconsistencies -> hardcoded x, y, z until this is solved. --- indica/readers/st40reader.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/indica/readers/st40reader.py b/indica/readers/st40reader.py index dc92d5ad..e4b2e7ec 100644 --- a/indica/readers/st40reader.py +++ b/indica/readers/st40reader.py @@ -3,6 +3,7 @@ """ +from copy import deepcopy from typing import Any from typing import Dict from typing import List @@ -1095,6 +1096,11 @@ def _get_thomson_scattering( y, y_path = self._get_signal(uid, instrument, ":y", revision) z, z_path = self._get_signal(uid, instrument, ":z", revision) R, R_path = self._get_signal(uid, instrument, ":R", revision) + # TODO: hardcoded correction to TS coordinates to be fixed in MDS+ + print("\n Hardcoded correction to TS coordinates to be fixed in MDS+ \n") + z = R * 0.0 + x = deepcopy(R) + y = 0 for q in quantities: qval, q_path = self._get_signal( From ef8512fde2577f4240bcca69d9577c76575edbec Mon Sep 17 00:00:00 2001 From: Marco Sertoli Date: Mon, 9 Oct 2023 09:08:28 +0100 Subject: [PATCH 12/16] Working implementation of simple spline fitting routine + example for TS data. --- indica/operators/spline_fit_easy.py | 90 +++++++++++++++++++---------- 1 file changed, 59 insertions(+), 31 deletions(-) diff --git a/indica/operators/spline_fit_easy.py b/indica/operators/spline_fit_easy.py index a012f5c0..d50dc463 100644 --- a/indica/operators/spline_fit_easy.py +++ b/indica/operators/spline_fit_easy.py @@ -99,63 +99,91 @@ def residuals(yknots): def spline_fit_ts( - pulse: int, - tstart: float = 0.0, - tend: float = 0.2, - plot: bool = False, + pulse: int = 11314, + tstart: float = 0.03, + tend: float = 0.1, + dt: float = 0.01, + quantity: str = "te", + R_shift: float = 0.0, + knots: list = None, + plot: bool = True, ): + st40 = ReadST40(pulse, tstart=tstart, tend=tend, dt=dt) + st40(["ts"], R_shift=R_shift) + + if quantity == "te" and knots is None: + knots = [0, 0.3, 0.6, 0.95, 1.1] + if quantity == "ne" and knots is None: + knots = [0, 0.3, 0.6, 0.8, 0.95, 1.1] + data_all = st40.raw_data["ts"][quantity] + t = data_all.t + transform = data_all.transform + transform.convert_to_rho_theta(t=data_all.t) - ST40 = ReadST40(pulse, tstart=tstart, tend=tend) - ST40(["ts"]) - - Te_data_all = ST40.binned_data["ts"]["te"] - t = ST40.binned_data["ts"]["te"].t - transform = ST40.binned_data["ts"]["te"].transform R = transform.R Rmag = transform.equilibrium.rmag.interp(t=t) # Fit all available TS data - ind = np.full_like(Te_data_all, True) + ind = np.full_like(data_all, True) rho = xr.where(ind, transform.rho, np.nan) - Te_data = xr.where(ind, Te_data_all, np.nan) - Te_err = xr.where(ind, Te_data_all.error, np.nan) - Te_fit = fit_profile( - rho, Te_data, Te_err, knots=[0, 0.3, 0.5, 0.75, 0.95, 1.05], virtual_knots=False - ) + data = xr.where(ind, data_all, np.nan) + err = xr.where(ind, data_all.error, np.nan) + fit = fit_profile(rho, data, err, knots=knots, virtual_knots=False) # Use only HFS channels ind = R <= Rmag - rho = xr.where(ind, transform.rho, np.nan) - Te_data = xr.where(ind, Te_data_all, np.nan) - Te_err = xr.where(ind, Te_data_all.error, np.nan) - Te_fit_hfs = fit_profile(rho, Te_data, Te_err, virtual_knots=True) + rho_hfs = xr.where(ind, transform.rho, np.nan) + data_hfs = xr.where(ind, data_all, np.nan) + err_hfs = xr.where(ind, data_all.error, np.nan) + fit_hfs = fit_profile(rho_hfs, data_hfs, err_hfs, knots=knots, virtual_knots=True) + + # Use only LFS channels + ind = R >= Rmag + rho_lfs = xr.where(ind, transform.rho, np.nan) + data_lfs = xr.where(ind, data_all, np.nan) + err_lfs = xr.where(ind, data_all.error, np.nan) + fit_lfs = fit_profile(rho_lfs, data_lfs, err_lfs, knots=knots, virtual_knots=True) if plot: - for t in Te_data_all.t: + for t in data_all.t: plt.ioff() plt.errorbar( - Te_data_all.transform.rho.sel(t=t), - Te_data_all.sel(t=t), - Te_data_all.error.sel(t=t), + rho_hfs.sel(t=t), + data_hfs.sel(t=t), + err_hfs.sel(t=t), marker="o", - label="data", + label="data HFS", + color="blue", + ) + plt.errorbar( + rho_lfs.sel(t=t), + data_lfs.sel(t=t), + err_lfs.sel(t=t), + marker="o", + label="data LFS", + color="red", + ) + fit.sel(t=t).plot( + linewidth=5, alpha=0.5, color="black", label="spline fit all" ) - Te_fit.sel(t=t).plot( - linewidth=5, alpha=0.5, color="orange", label="spline fit all" + fit_lfs.sel(t=t).plot( + linewidth=5, alpha=0.5, color="red", label="spline fit LFS" ) - Te_fit_hfs.sel(t=t).plot( - linewidth=5, alpha=0.5, color="red", label="spline fit HFS" + fit_hfs.sel(t=t).plot( + linewidth=5, alpha=0.5, color="blue", label="spline fit HFS" ) plt.legend() plt.show() - return Te_data_all, Te_fit + return data_all, fit if __name__ == "__main__": + plt.ioff() spline_fit_ts( - 10619, + 11314, tstart=0.0, tend=0.2, plot=True, ) + plt.show() From f50b613bbc02d7b521e493acc2131e4c92d6729e Mon Sep 17 00:00:00 2001 From: Marco Sertoli Date: Mon, 9 Oct 2023 09:09:00 +0100 Subject: [PATCH 13/16] Minor corrections to load_modelling.., but still too complex. --- indica/workflows/load_modelling_plasma.py | 38 +++++++++-------------- 1 file changed, 14 insertions(+), 24 deletions(-) diff --git a/indica/workflows/load_modelling_plasma.py b/indica/workflows/load_modelling_plasma.py index 38d06ecb..748c42f8 100644 --- a/indica/workflows/load_modelling_plasma.py +++ b/indica/workflows/load_modelling_plasma.py @@ -31,6 +31,7 @@ CMAP, COLORS = set_plot_colors() DIAGNOSTIC_MODELS = { + "smmh": Interferometry, "smmh1": Interferometry, "nirh1": Interferometry, "xrcs": HelikeSpectrometer, @@ -45,6 +46,7 @@ "blom_xy1": Bolometer, } INSTRUMENTS: list = [ + "smmh", "smmh1", "nirh1", "xrcs", @@ -93,7 +95,6 @@ def plasma_code( n_rad = len(data[runs[0]]["ne"].rho_poloidal) main_ion = "h" impurities = ("ar", "c", "he") - impurities = ("c", "ar", "he") impurity_concentration = (0.001, 0.03, 0.01) _plasma = Plasma( tstart=tstart, @@ -127,24 +128,11 @@ def plasma_code( for element in _plasma.elements: _plasma.ion_temperature.loc[dict(element=element)] = Ti.values for i, impurity in enumerate(_plasma.impurities): - if impurity == "c": - Nimp = ( - _data[f"niz{i+1}"].interp(rho_poloidal=_plasma.rho, t=_plasma.t) - * 1.0e19 - ) - _plasma.impurity_density.loc[dict(element=impurity)] = Nimp.values - elif impurity == "ar": - Nimp = ( - _data[f"niz{i}"].interp(rho_poloidal=_plasma.rho, t=_plasma.t) - * 1.0e19 - ) / 100 - _plasma.impurity_density.loc[dict(element=impurity)] = Nimp.values - else: - Nimp = ( - _data[f"niz{i+1}"].interp(rho_poloidal=_plasma.rho, t=_plasma.t) - * 1.0e19 - ) - _plasma.impurity_density.loc[dict(element=impurity)] = Nimp.values + Nimp = ( + _data[f"niz{i+1}"].interp(rho_poloidal=_plasma.rho, t=_plasma.t) + * 1.0e19 + ) + _plasma.impurity_density.loc[dict(element=impurity)] = Nimp.values Nf = _data["nf"].interp(rho_poloidal=_plasma.rho, t=_plasma.t) * 1.0e19 _plasma.fast_density.values = Nf.values @@ -229,9 +217,9 @@ def initialize_diagnostic_models( def example_run( dt: float = 0.01, verbose: bool = True, - example: str = "predictive", + example: str = "alsu_11314", all_runs: bool = False, - plot: bool = False, + plot: bool = True, save_fig: bool = False, fig_style="profiles", alpha: float = 1.0, @@ -539,8 +527,8 @@ def plot_data_bckc_comparison( norm["xrcs"]["spectra"] = True # norm["brems"] = True # norm["sxr_camera_4"] = True - norm["sxrc_xy2"] = {} - norm["sxrc_xy2"]["brightness"] = True + # norm["sxrc_xy2"] = {} + # norm["sxrc_xy2"]["brightness"] = True y0 = {} y0["nirh1"] = True y0["smmh1"] = True @@ -727,7 +715,7 @@ def example_params(example: str, all_runs: bool = False): code = "astra" if not all_runs: runs = ["RUN610", "RUN611", "RUN612"] # C and Ar - runs = ["RUN623"] # C and Ar + # runs = ["RUN623"] # C and Ar tstart = 0.03 tend = 0.11 tplot = 0.08 @@ -805,4 +793,6 @@ def example_params(example: str, all_runs: bool = False): if __name__ == "__main__": + plt.ioff() example_run() + plt.show() From d3ccd96a78b9bc088fbfe0d432c78b77acb10d05 Mon Sep 17 00:00:00 2001 From: Marco Sertoli Date: Mon, 9 Oct 2023 12:48:45 +0100 Subject: [PATCH 14/16] Added TODO notes for things that need to be changed. Issues should be raised once these CAN be changed. --- indica/readers/available_quantities.py | 55 ++++++-------------------- indica/readers/st40reader.py | 25 +++++------- 2 files changed, 23 insertions(+), 57 deletions(-) diff --git a/indica/readers/available_quantities.py b/indica/readers/available_quantities.py index bb5d9fb9..c0c61bf7 100644 --- a/indica/readers/available_quantities.py +++ b/indica/readers/available_quantities.py @@ -1,5 +1,6 @@ """ Quantities that can be read with the current abstract reader implementation +TODO: change the tuple to DataArray (long_name, units) - see examples in abstractreader """ from typing import Dict @@ -13,9 +14,7 @@ "te": ("temperature", "electrons"), "chi2": ("chi-squared", "fit"), }, - "get_spectrometer": { - "spectra": ("emission", "spectral"), - }, + "get_spectrometer": {"spectra": ("emission", "spectral"),}, "get_charge_exchange": { # "angf": ("angular_freq", "ion"), "vtor": ("linear_rotation", "ion"), @@ -23,9 +22,7 @@ "spectra": ("spectra", "experimental"), "fit": ("spectra", "fit"), }, - "get_bremsstrahlung_spectroscopy": { - "zeff": ("effective_charge", "plasma"), - }, + "get_bremsstrahlung_spectroscopy": {"zeff": ("effective_charge", "plasma"),}, "get_helike_spectroscopy": { "int_w": ("intensity", "spectral_line"), "int_k": ("intensity", "spectral_line"), @@ -37,12 +34,8 @@ "ti_z": ("temperature", "ions"), "spectra": ("spectra", "passive"), }, - "get_diode_filters": { - "brightness": ("luminous_flux", None), - }, - "get_interferometry": { - "ne": ("density", "electrons"), - }, + "get_diode_filters": {"brightness": ("luminous_flux", None),}, + "get_interferometry": {"ne": ("density", "electrons"),}, "get_equilibrium": { "f": ("f_value", "plasma"), "faxs": ("magnetic_flux_axis", "poloidal"), @@ -62,12 +55,8 @@ "wp": ("energy", "plasma"), "psin": ("poloidal_flux", "normalised"), }, - "get_cyclotron_emissions": { - "te": ("temperature", "electrons"), - }, - "get_radiation": { - "brightness": ("luminous_flux", None), - }, + "get_cyclotron_emissions": {"te": ("temperature", "electrons"),}, + "get_radiation": {"brightness": ("luminous_flux", None),}, "get_astra": { "f": ("f_value", "plasma"), "faxs": ("magnetic_flux_axis", "poloidal"), @@ -88,22 +77,10 @@ "rbnd": ("major_rad", "separatrix"), "zbnd": ("z", "separatrix"), "ipla": ("current", "plasma"), - "upl": ( - "voltage", - "loop", - ), # Loop voltage V - "wth": ( - "stored_energy", - "equilibrium", - ), - "wtherm": ( - "stored_energy", - "thermal", - ), - "wfast": ( - "stored_energy", - "fast", - ), # Thermal stored energy + "upl": ("voltage", "loop",), # Loop voltage V + "wth": ("stored_energy", "equilibrium",), + "wtherm": ("stored_energy", "thermal",), + "wfast": ("stored_energy", "fast",), # Thermal stored energy "j_bs": ("current_density", "bootstrap"), # Bootstrap current density,MA/m2 "j_nbi": ( "current_density", @@ -127,14 +104,8 @@ "electron", ), # Beam power density to electrons, MW/m3 "qnbi": ("power_density_nbi", "ion"), # Beam power density to ions, MW/m3 - "q_oh": ( - "power_density_ohm", - "total", - ), # Ohmic heating power profile, MW/m3 - "q_rf": ( - "power_density_rf", - "electron", - ), # RF power density to electron,MW/m3 + "q_oh": ("power_density_ohm", "total",), # Ohmic heating power profile, MW/m3 + "q_rf": ("power_density_rf", "electron",), # RF power density to electron,MW/m3 "sbm": ("particle_source", "nbi"), # Particle source from beam, 10^19/m^3/s "swall": ( "particle_source", diff --git a/indica/readers/st40reader.py b/indica/readers/st40reader.py index ff149fab..60299ad2 100644 --- a/indica/readers/st40reader.py +++ b/indica/readers/st40reader.py @@ -97,6 +97,7 @@ class ST40Reader(DataReader): "tws_c": "get_spectrometer", "ts": "get_thomson_scattering", } + #TODO: this will not be necessary once the MDS+ standardisation is complete UIDS_MDS = { "xrcs": "sxr", "princeton": "spectrom", @@ -297,7 +298,8 @@ class ST40Reader(DataReader): }, } - _IMPLEMENTATION_QUANTITIES = { # TODO: these will be different diagnostics!!!!!!! + # TODO: this can be deleted once MDS+ standardisation is complete + _IMPLEMENTATION_QUANTITIES = { "diode_arrays": { # GETTING THE DATA OF THE SXR CAMERA "sxr_camera_1": ("brightness", "total"), "sxr_camera_2": ("brightness", "50_Al_filtered"), @@ -306,6 +308,7 @@ class ST40Reader(DataReader): }, } + # TODO: this can be deleted once MDS+ standardisation is complete _RADIATION_RANGES = { "sxr_camera_1": (1, 20), "sxr_camera_2": (21, 40), @@ -769,6 +772,7 @@ def _get_charge_exchange( z, z_path = self._get_signal(uid, instrument, ":z", revision) R, R_path = self._get_signal(uid, instrument, ":R", revision) + # TODO: temporary fix until geometry sorted (especially pulse if statement..) try: location, location_path = self._get_signal( uid, instrument, ".geometry:location", revision @@ -780,7 +784,6 @@ def _get_charge_exchange( location = np.array([location]) direction = np.array([direction]) - # TODO: temporary fix until geometry sorted if location.shape[0] != x.shape[0]: if self.pulse > 10200: index = np.arange(18, 36) @@ -788,7 +791,6 @@ def _get_charge_exchange( index = np.arange(21, 36) location = location[index] direction = direction[index] - except TreeNNF: location = None direction = None @@ -813,10 +815,8 @@ def _get_charge_exchange( ) except TreeNNF: qval_err = np.full_like(qval, 0.0) - # q_path_err = "" dimensions, _ = self._get_signal_dims(q_path, len(qval.shape)) - results[q + "_records"] = q_path results[q] = qval results[f"{q}_error"] = qval_err @@ -869,13 +869,6 @@ def _get_spectrometer( location = np.array([location]) direction = np.array([direction]) - # if self.pulse > 10200: - # index = np.arange(18, 36) - # else: - # index = np.arange(21, 36) - # location = location[index] - # direction = direction[index] - for q in quantities: qval, q_path = self._get_signal( uid, @@ -893,7 +886,6 @@ def _get_spectrometer( ) except TreeNNF: qval_err = np.full_like(qval, 0.0) - # q_path_err = "" dimensions, _ = self._get_signal_dims(q_path, len(qval.shape)) @@ -954,6 +946,8 @@ def _get_diode_filters( _labels, _ = self._get_signal(uid, instrument, ":label", revision) if type(_labels[0]) == np.bytes_: labels = np.array([label.decode("UTF-8") for label in _labels]) + else: + labels = _labels results["times"] = times results["labels"] = labels @@ -1086,6 +1080,8 @@ def _get_thomson_scattering( revision = results["revision"] times, times_path = self._get_signal(uid, instrument, ":time", revision) + # TODO: hardcoded correction to TS coordinates to be fixed in MDS+ + print("\n Hardcoded correction to TS coordinates to be fixed in MDS+ \n") # location, location_path = self._get_signal( # uid, instrument, ".geometry:location", revision # ) @@ -1096,8 +1092,6 @@ def _get_thomson_scattering( y, y_path = self._get_signal(uid, instrument, ":y", revision) z, z_path = self._get_signal(uid, instrument, ":z", revision) R, R_path = self._get_signal(uid, instrument, ":R", revision) - # TODO: hardcoded correction to TS coordinates to be fixed in MDS+ - print("\n Hardcoded correction to TS coordinates to be fixed in MDS+ \n") z = R * 0.0 x = deepcopy(R) y = 0 @@ -1176,6 +1170,7 @@ def get_revision_name(self, revision) -> str: """Return string defining RUN## or BEST if revision = 0""" if type(revision) == int: + rev_str = "" if revision < 0: rev_str = "" elif revision == 0: From 738d75fb35dc4d2e31076b12ef750fe55db75e95 Mon Sep 17 00:00:00 2001 From: Marco Sertoli Date: Mon, 9 Oct 2023 12:49:14 +0100 Subject: [PATCH 15/16] Ran pre-commit --- indica/operators/spline_fit_easy.py | 25 +++++-------------------- 1 file changed, 5 insertions(+), 20 deletions(-) diff --git a/indica/operators/spline_fit_easy.py b/indica/operators/spline_fit_easy.py index d50dc463..b2c5b13a 100644 --- a/indica/operators/spline_fit_easy.py +++ b/indica/operators/spline_fit_easy.py @@ -22,14 +22,9 @@ def fit_profile( """Fit a profile""" def residuals(yknots): - spline = CubicSpline( - xknots, - yknots, - axis=0, - bc_type=bc_type, - ) + spline = CubicSpline(xknots, yknots, axis=0, bc_type=bc_type,) bckc = np.interp(x, xspl, spline(xspl)) - residuals = np.sqrt((y - bckc) ** 2 / err**2) + residuals = np.sqrt((y - bckc) ** 2 / err ** 2) return residuals yspl = xr.DataArray( @@ -81,12 +76,7 @@ def residuals(yknots): bounds=(lower_bound, upper_bound), verbose=verbose, ) - spline = CubicSpline( - xknots, - fit.x, - axis=0, - bc_type=bc_type, - ) + spline = CubicSpline(xknots, fit.x, axis=0, bc_type=bc_type,) _yspl = spline(xspl) except ValueError: _yspl = np.full_like(xspl, 0.0) @@ -112,7 +102,7 @@ def spline_fit_ts( st40(["ts"], R_shift=R_shift) if quantity == "te" and knots is None: - knots = [0, 0.3, 0.6, 0.95, 1.1] + knots = [0, 0.3, 0.6, 0.8, 1.1] if quantity == "ne" and knots is None: knots = [0, 0.3, 0.6, 0.8, 0.95, 1.1] data_all = st40.raw_data["ts"][quantity] @@ -180,10 +170,5 @@ def spline_fit_ts( if __name__ == "__main__": plt.ioff() - spline_fit_ts( - 11314, - tstart=0.0, - tend=0.2, - plot=True, - ) + spline_fit_ts(11089, quantity="ne") plt.show() From cb7908ae4395ac488798e1d1e3ec37fad428b688 Mon Sep 17 00:00:00 2001 From: Marco Sertoli Date: Wed, 11 Oct 2023 16:31:56 +0100 Subject: [PATCH 16/16] Minor mods and to decrease complexity and pass linting --- indica/operators/spline_fit_easy.py | 16 +- indica/readers/available_quantities.py | 54 +++- indica/readers/st40reader.py | 2 +- indica/workflows/load_modelling_plasma.py | 301 +++++++++++----------- 4 files changed, 211 insertions(+), 162 deletions(-) diff --git a/indica/operators/spline_fit_easy.py b/indica/operators/spline_fit_easy.py index b2c5b13a..d6549067 100644 --- a/indica/operators/spline_fit_easy.py +++ b/indica/operators/spline_fit_easy.py @@ -22,9 +22,14 @@ def fit_profile( """Fit a profile""" def residuals(yknots): - spline = CubicSpline(xknots, yknots, axis=0, bc_type=bc_type,) + spline = CubicSpline( + xknots, + yknots, + axis=0, + bc_type=bc_type, + ) bckc = np.interp(x, xspl, spline(xspl)) - residuals = np.sqrt((y - bckc) ** 2 / err ** 2) + residuals = np.sqrt((y - bckc) ** 2 / err**2) return residuals yspl = xr.DataArray( @@ -76,7 +81,12 @@ def residuals(yknots): bounds=(lower_bound, upper_bound), verbose=verbose, ) - spline = CubicSpline(xknots, fit.x, axis=0, bc_type=bc_type,) + spline = CubicSpline( + xknots, + fit.x, + axis=0, + bc_type=bc_type, + ) _yspl = spline(xspl) except ValueError: _yspl = np.full_like(xspl, 0.0) diff --git a/indica/readers/available_quantities.py b/indica/readers/available_quantities.py index c0c61bf7..133d892e 100644 --- a/indica/readers/available_quantities.py +++ b/indica/readers/available_quantities.py @@ -14,7 +14,9 @@ "te": ("temperature", "electrons"), "chi2": ("chi-squared", "fit"), }, - "get_spectrometer": {"spectra": ("emission", "spectral"),}, + "get_spectrometer": { + "spectra": ("emission", "spectral"), + }, "get_charge_exchange": { # "angf": ("angular_freq", "ion"), "vtor": ("linear_rotation", "ion"), @@ -22,7 +24,9 @@ "spectra": ("spectra", "experimental"), "fit": ("spectra", "fit"), }, - "get_bremsstrahlung_spectroscopy": {"zeff": ("effective_charge", "plasma"),}, + "get_bremsstrahlung_spectroscopy": { + "zeff": ("effective_charge", "plasma"), + }, "get_helike_spectroscopy": { "int_w": ("intensity", "spectral_line"), "int_k": ("intensity", "spectral_line"), @@ -34,8 +38,12 @@ "ti_z": ("temperature", "ions"), "spectra": ("spectra", "passive"), }, - "get_diode_filters": {"brightness": ("luminous_flux", None),}, - "get_interferometry": {"ne": ("density", "electrons"),}, + "get_diode_filters": { + "brightness": ("luminous_flux", None), + }, + "get_interferometry": { + "ne": ("density", "electrons"), + }, "get_equilibrium": { "f": ("f_value", "plasma"), "faxs": ("magnetic_flux_axis", "poloidal"), @@ -55,8 +63,12 @@ "wp": ("energy", "plasma"), "psin": ("poloidal_flux", "normalised"), }, - "get_cyclotron_emissions": {"te": ("temperature", "electrons"),}, - "get_radiation": {"brightness": ("luminous_flux", None),}, + "get_cyclotron_emissions": { + "te": ("temperature", "electrons"), + }, + "get_radiation": { + "brightness": ("luminous_flux", None), + }, "get_astra": { "f": ("f_value", "plasma"), "faxs": ("magnetic_flux_axis", "poloidal"), @@ -77,10 +89,22 @@ "rbnd": ("major_rad", "separatrix"), "zbnd": ("z", "separatrix"), "ipla": ("current", "plasma"), - "upl": ("voltage", "loop",), # Loop voltage V - "wth": ("stored_energy", "equilibrium",), - "wtherm": ("stored_energy", "thermal",), - "wfast": ("stored_energy", "fast",), # Thermal stored energy + "upl": ( + "voltage", + "loop", + ), # Loop voltage V + "wth": ( + "stored_energy", + "equilibrium", + ), + "wtherm": ( + "stored_energy", + "thermal", + ), + "wfast": ( + "stored_energy", + "fast", + ), # Thermal stored energy "j_bs": ("current_density", "bootstrap"), # Bootstrap current density,MA/m2 "j_nbi": ( "current_density", @@ -104,8 +128,14 @@ "electron", ), # Beam power density to electrons, MW/m3 "qnbi": ("power_density_nbi", "ion"), # Beam power density to ions, MW/m3 - "q_oh": ("power_density_ohm", "total",), # Ohmic heating power profile, MW/m3 - "q_rf": ("power_density_rf", "electron",), # RF power density to electron,MW/m3 + "q_oh": ( + "power_density_ohm", + "total", + ), # Ohmic heating power profile, MW/m3 + "q_rf": ( + "power_density_rf", + "electron", + ), # RF power density to electron,MW/m3 "sbm": ("particle_source", "nbi"), # Particle source from beam, 10^19/m^3/s "swall": ( "particle_source", diff --git a/indica/readers/st40reader.py b/indica/readers/st40reader.py index 60299ad2..92f567c2 100644 --- a/indica/readers/st40reader.py +++ b/indica/readers/st40reader.py @@ -97,7 +97,7 @@ class ST40Reader(DataReader): "tws_c": "get_spectrometer", "ts": "get_thomson_scattering", } - #TODO: this will not be necessary once the MDS+ standardisation is complete + # TODO: this will not be necessary once the MDS+ standardisation is complete UIDS_MDS = { "xrcs": "sxr", "princeton": "spectrom", diff --git a/indica/workflows/load_modelling_plasma.py b/indica/workflows/load_modelling_plasma.py index 748c42f8..d3167b0b 100644 --- a/indica/workflows/load_modelling_plasma.py +++ b/indica/workflows/load_modelling_plasma.py @@ -19,7 +19,6 @@ from indica.models.plasma import Plasma from indica.models.sxr_camera import SXRcamera from indica.models.thomson_scattering import ThomsonScattering -from indica.numpy_typing import RevisionLike from indica.readers.read_st40 import ReadST40 from indica.utilities import save_figure from indica.utilities import set_axis_sci @@ -541,39 +540,38 @@ def plot_data_bckc_comparison( for run in runs: bckc[run]["efit"] = {"wp": plasma[run].wp} - for instrument in st40.raw_data.keys(): - for quantity in bckc[run][instrument].keys(): - print(instrument) - print(f" {quantity}") - run = runs[0] - + instruments = st40.raw_data.keys() + for instrument in instruments: + quantities = bckc[runs[0]][instrument].keys() + for quantity in quantities: if ( quantity not in st40.binned_data[instrument].keys() or quantity not in st40.raw_data[instrument].keys() ): continue - plt.figure() + print(instrument) + print(f" {quantity}") _raw = st40.raw_data[instrument][quantity] _binned = st40.binned_data[instrument][quantity] - _bckc = bckc[run][instrument][quantity] + _bckc = bckc[runs[0]][instrument][quantity] + tslice = slice(_bckc.t.min().values, _bckc.t.max().values) + str_to_add = "" + tslice_raw = tslice + tslice_binned = tslice + if "error" not in _binned.attrs: _binned.attrs["error"] = xr.full_like(_binned, 0.0) if "stdev" not in _binned.attrs: _binned.attrs["stdev"] = xr.full_like(_binned, 0.0) - err = np.sqrt(_binned.error**2 + _binned.stdev**2) err = xr.where(err / _binned.values < 1.0, err, 0.0) - if len(_bckc.dims) > 1: + if len(_binned.dims) > 1: str_to_add = f" @ {time:.3f} s" tslice_binned = _binned.t.sel(t=time, method="nearest") tslice_raw = _raw.t.sel(t=time, method="nearest") - else: - str_to_add = "" - tslice_raw = tslice - tslice_binned = tslice _raw = _raw.sel(t=tslice_raw) _binned = _binned.sel(t=tslice_binned) @@ -585,6 +583,7 @@ def plot_data_bckc_comparison( _binned -= bgnd _raw -= bgnd + plt.figure() _raw.plot( label="Raw", color=COLORS["raw"], @@ -606,11 +605,8 @@ def plot_data_bckc_comparison( markersize=markersize, ) - label = None + label: str = "Model" for run in runs: - if run == runs[-1]: - label = "Model" - _bckc = bckc[run][instrument][quantity].sel(t=tslice_binned) if instrument in norm.keys(): if quantity in norm[instrument].keys(): @@ -624,138 +620,151 @@ def plot_data_bckc_comparison( linewidth=rcParams["lines.linewidth"] * 1.5, alpha=alpha, ) - set_axis_sci() - plt.title(f"{instrument.upper()} {quantity}" + str_to_add) - if instrument in y0.keys(): - plt.ylim( - 0, - ) + del label - if quantity == "spectra": - # TODO: wavelength axis is sorted from max to min... - plt.xlim(_bckc.wavelength.min(), _bckc.wavelength.max()) + set_axis_sci() + plt.title(f"{instrument.upper()} {quantity}" + str_to_add) + if instrument in y0.keys(): + plt.ylim( + 0, + ) + if quantity == "spectra": + plt.xlim(_bckc.wavelength.min(), _bckc.wavelength.max()) - plt.legend() - save_figure(fig_path, f"{instrument}_{quantity}", save_fig=save_fig) + plt.legend() + save_figure(fig_path, f"{instrument}_{quantity}", save_fig=save_fig) def example_params(example: str, all_runs: bool = False): - comment: str - pulse_code: int - pulse: int - equil: str - code: str - tstart: float - tend: float - tplot: float - runs = [f"RUN{run}" for run in (500 + np.arange(61, 77))] - equil_run: RevisionLike = 0 - - if example == "predictive": - comment = "Tests using fixed-boundary predictive ASTRA" - pulse_code = 13110009 - pulse = 10009 - equil = "astra" - code = "astra" - if not all_runs: - runs = ["RUN2621"] - tstart = 0.02 - tend = 0.08 - tplot = 0.06 - elif example == "interpretative_10009": - comment = "interpretative ASTRA using HDA/EFIT" - pulse_code = 13110009 - pulse = 10009 - equil = "efit" - code = "astra" - if not all_runs: - runs = ["RUN564"] - tstart = 0.03 - tend = 0.1 - tplot = 0.06 - elif example == "interpretative_9850": - comment = "ASTRA interpretative using HDA/EFIT" - pulse = 9850 - pulse_code = 13109850 - equil = "efit" - code = "astra" - if not all_runs: - runs = ["RUN564"] - tstart = 0.02 - tend = 0.1 - tplot = 0.08 - elif example == "interpretative_9229": - comment = "ASTRA interpretative using HDA/EFIT" - pulse = 9229 - pulse_code = 13109229 - equil = "efit" - code = "astra" - if not all_runs: - runs = ["RUN572"] - tstart = 0.03 - tend = 0.11 - tplot = 0.06 - elif example == "diverted": - comment = "predictive ASTRA using for diverted scenario" - pulse_code = 13000040 - pulse = 10009 - equil = "astra" - code = "astra" - if not all_runs: - runs = ["RUN292"] - tstart = 0.03 - tend = 0.11 - tplot = 0.1 - elif example == "aleksei_11228": - comment = "ASTRA using TS and invented Ti shapes" - pulse = 11228 - pulse_code = 13011228 - equil = "efit" # "astra" - equil_run = 1 - code = "astra" - if not all_runs: - runs = ["RUN610", "RUN611", "RUN612"] # C and Ar - # runs = ["RUN623"] # C and Ar - tstart = 0.03 - tend = 0.11 - tplot = 0.08 - elif example == "alsu_11312": - comment = "ASTRA using TS and peaked Ti scaled to CXRS" - pulse = 11312 - pulse_code = 33011312 - equil = "astra" - code = "astra" - if not all_runs: - runs = ["RUN21"] - tstart = 0.065 - tend = 0.095 - tplot = 0.075 - elif example == "alsu_11314": - comment = "ASTRA using TS and peaked Ti scaled to CXRS" - pulse = 11314 - pulse_code = 33011314 - equil = "astra" - code = "astra" - if not all_runs: - runs = ["RUN12"] - tstart = 0.065 - tend = 0.095 - tplot = 0.075 - elif example == "alsu_11317": - comment = "ASTRA using TS and peaked Ti scaled to CXRS" - pulse = 11317 - pulse_code = 33011317 - equil = "astra" - code = "astra" - if not all_runs: - runs = ["RUN9"] - tstart = 0.065 - tend = 0.095 - tplot = 0.075 - else: - raise ValueError(f"No parameters for example {example}") - - return pulse_code, pulse, equil, equil_run, code, runs, comment, tstart, tend, tplot + runs_all: list = [f"RUN{run}" for run in (500 + np.arange(61, 77))] + + params = { + "predictive": dict( + comment="Tests using fixed-boundary predictive ASTRA", + pulse_code=13110009, + pulse=10009, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN2621"], + tstart=0.02, + tend=0.08, + tplot=0.06, + ), + "interpretative_10009": dict( + comment="interpretative ASTRA using HDA/EFIT", + pulse_code=13110009, + pulse=10009, + equil="efit", + equil_run=0, + code="astra", + runs=["RUN564"], + tstart=0.03, + tend=0.1, + tplot=0.06, + ), + "interpretative_9850": dict( + comment="ASTRA interpretative using HDA/EFIT", + pulse=9850, + pulse_code=13109850, + equil="efit", + equil_run=0, + code="astra", + runs=["RUN564"], + tstart=0.02, + tend=0.1, + tplot=0.08, + ), + "interpretative_9229": dict( + comment="ASTRA interpretative using HDA/EFIT", + pulse=9229, + pulse_code=13109229, + equil="efit", + equil_run=0, + code="astra", + runs=["RUN572"], + tstart=0.03, + tend=0.11, + tplot=0.06, + ), + "diverted": dict( + comment="predictive ASTRA using for diverted scenario", + pulse_code=13000040, + pulse=10009, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN292"], + tstart=0.03, + tend=0.11, + tplot=0.1, + ), + "aleksei_11228": dict( + comment="ASTRA using TS and invented Ti shapes", + pulse=11228, + pulse_code=13011228, + equil="efit", # "astra" + equil_run=1, + code="astra", + runs=["RUN610", "RUN611", "RUN612"], + tstart=0.03, + tend=0.11, + tplot=0.08, + ), + "alsu_11312": dict( + comment="ASTRA using TS and peaked Ti scaled to CXRS", + pulse=11312, + pulse_code=33011312, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN21"], + tstart=0.065, + tend=0.095, + tplot=0.075, + ), + "alsu_11314": dict( + comment="ASTRA using TS and peaked Ti scaled to CXRS", + pulse=11314, + pulse_code=33011314, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN12"], + tstart=0.065, + tend=0.095, + tplot=0.075, + ), + "alsu_11317": dict( + comment="ASTRA using TS and peaked Ti scaled to CXRS", + pulse=11317, + pulse_code=33011317, + equil="astra", + equil_run=0, + code="astra", + runs=["RUN9"], + tstart=0.065, + tend=0.095, + tplot=0.075, + ), + } + + _params = params[example] + if all_runs: + _params["runs"] = runs_all + + return ( + _params["pulse_code"], + _params["pulse"], + _params["equil"], + _params["equil_run"], + _params["code"], + _params["runs"], + _params["comment"], + _params["tstart"], + _params["tend"], + _params["tplot"], + ) # def add_gacode_data(