Search

Equilibriums in magnetic confinement devices result from force balancing between the Lorentz force and the plasma pressure gradient. In an axisymmetric configuration like a tokamak, such an equilibrium is described by an elliptic equation for the poloidal magnetic flux, commonly known as the Grad--Shafranov equation. It is challenging to develop a scalable and accurate free-boundary Grad--Shafranov solver, since it is a fully nonlinear optimization problem that simultaneouly solves for the magnetic field coil current outside the plasma to control the plasma shape. In this work, we develop a Newton-based free-boundary Grad--Shafranov solver using adaptive finite elements and preconditioning strategies. The free-boundary interaction leads to the evaluation of a domain-dependent nonlinear form of which its contribution to the Jacobian matrix is achieved through shape calculus. The optimization problem aims to minimize the distance between the plasma boundary and specified control points while satisfying two non-trivial constraints, which correspond to the nonlinear finite element discretization of the Grad--Shafranov equation and a constraint on the total plasma current involving a nonlocal coupling term. The linear system is solved by a block factorization, and AMG is called for subblock elliptic operators. The unique contributions of this work include the treatment of a global constraint, preconditioning strategies, nonlocal reformulation, and the implementation of adaptive finite elements. It is found that the resulting Newton solver is robust, successfully reducing the nonlinear residual to 1e-6 and lower in a small handful of iterations while addressing the challenging case to find a Taylor state equilibrium where conventional Picard-based solvers fail to converge.

Given a data graph G, a source vertex u and a target vertex v of a reachability query, the reachability query is used to answer whether there exists a path from u to v in G. Reachability query processing is one of the fundamental operations in graph data management, which is widely used in biological networks, communication networks, and social networks to assist data analysis. The data graphs in practical applications usually contain information such as quantization weights associated with the structural relationships, in addition to the structural relationships between vertices. Thus, in addition to the traditional reachability relationships, users may want to further understand whether such reachability relationships satisfy specific constraints. In this paper, we study the problem of efficiently processing k -step reachability queries with weighted constraints in weighted graphs. The k -step weighted reachability query questions are used to answer the question of whether there exists a path from a source vertex u to a goal vertex v in a given weighted graph. If it exists, the path needs to satisfy 1) all edges in the path satisfy the given weight constraints, and 2) the length of the path does not exceed the given distance threshold k. To address the problem, firstly, WKRI index supporting k -step weighted reachability query processing and index construction methods based on efficient pruning strategies are proposed. Secondly, the idea of constructing index based on part of the vertexs is proposed to reduce the size of the index. We design and implement two optimized indexes GWKRI and LWKRI based on the vertex coverage set. Finally, experiments are conducted on several real datasets. The experimental results verify the efficiency of the method proposed in this paper in answering k -step weighted reachability queries.

A physics-informed neural network (PINN) is used to evaluate the fast ion distribution in the hot spot of an inertial confinement fusion target. The use of tailored input and output layers to the neural network is shown to enable a PINN to learn the parametric solution to the Vlasov-Fokker-Planck equation in the absence of any synthetic or experimental data. As an explicit demonstration of the approach, the specific problem of Knudsen layer fusion yield reduction is treated. Here, predictions from the Vlasov-Fokker-Planck PINN are used to provide a non-perturbative solution of the fast ion tail in the vicinity of the hot spot thus allowing the spatial profile of the fusion reactivity to be evaluated for a range of collisionalities and hot spot conditions. Excellent agreement is found between the predictions of the Vlasov-Fokker-Planck PINN and results from traditional numerical solvers with respect to both the energy and spatial distribution of fast ions and the fusion reactivity profile demonstrating that the Vlasov-Fokker-Planck PINN provides an accurate and efficient means of determining the impact of Knudsen layer yield reduction across a broad range of plasma conditions.

Stand-off runaway electron termination by injected tungsten particulates offers a plausible option in the toolbox of disruption mitigation. Tungsten is an attractive material choice for this application due to large electron stopping power and high melting point. To assess the feasibility of this scheme, we simulate runaway collisions with tungsten particulates using the MCNP program for incident runaway energies ranging from 1 to 10 MeV. We assess runaway termination from energetics and collisional kinematics perspectives. Energetically, the simulations show that 99% of runaway beam energy is removed by tungsten particulates on a timescale of 4-9 $\mu$s. Kinematically, the simulations show that 99% of runaways are terminated by absorption or backscattering on a timescale of 3-4 $\mu$s. By either metric, the runaway beam is effectively terminated before the onset of particulate melting. Furthermore, the simulations show that secondary radiation emission by tungsten particulates does not significantly impact the runaway termination efficacy of this scheme. Secondary radiation is emitted at lower particle energies than the incident runaways and with a broad angular distribution such that the majority of secondary electrons emitted will not experience efficient runaway re-acceleration. Overall, the stand-off runaway termination scheme is a promising concept for last-ditch runaway mitigation in ITER, SPARC, and other future burning-plasma tokamaks., Comment: Submitted to: Nuclear Fusion - 16 pages (4 supplementary), 11 figures (5 supplementary), 4 tables

The simulation of plasma physics is computationally expensive because the underlying physical system is of high dimensions, requiring three spatial dimensions and three velocity dimensions. One popular numerical approach is Particle-In-Cell (PIC) methods owing to its ease of implementation and favorable scalability in high-dimensional problems. An unfortunate drawback of the method is the introduction of statistical noise resulting from the use of finitely many particles. In this paper we examine the application of the Smoothness-Increasing Accuracy-Conserving (SIAC) family of convolution kernel filters as denoisers for moment data arising from PIC simulations. We show that SIAC filtering is a promising tool to denoise PIC data in the physical space as well as capture the appropriate scales in the Fourier space. Furthermore, we demonstrate how the application of the SIAC technique reduces the amount of information necessary in the computation of quantities of interest in plasma physics such as the Bohm speed.

In this work we consider a relativistic drift-kinetic model for runaway electrons along with a Fokker-Planck operator for small-angle Coulomb collisions, a radiation damping operator, and a secondary knock-on (Boltzmann) collision source. We develop a new scalable fully implicit solver utilizing finite volume and conservative finite difference schemes and dynamic mesh adaptivity. A new data management framework in the PETSc library based on the p4est library is developed to enable simulations with dynamic adaptive mesh refinement (AMR), distributed memory parallelization, and dynamic load balancing of computational work. This framework and the runaway electron solver building on the framework are able to dynamically capture both bulk Maxwellian at the low-energy region and a runaway tail at the high-energy region. To effectively capture features via the AMR algorithm, a new AMR indicator prediction strategy is proposed that is performed alongside the implicit time evolution of the solution. This strategy is complemented by the introduction of computationally cheap feature-based AMR indicators that are analyzed theoretically. Numerical results quantify the advantages of the prediction strategy in better capturing features compared with nonpredictive strategies; and we demonstrate trade-offs regarding computational costs. The robustness with respect to model parameters, algorithmic scalability, and parallel scalability are demonstrated through several benchmark problems including manufactured solutions and solutions of different physics models. We focus on demonstrating the advantages of using implicit time stepping and AMR for runaway electron simulations., Comment: accepted

Force-free plasmas are a good approximation where the plasma pressure is tiny compared with the magnetic pressure, which is the case during the cold vertical displacement event (VDE) of a major disruption in a tokamak. On time scales long compared with the transit time of Alfven waves, the evolution of a force-free plasma is most efficiently described by the quasi-static magnetohydrodynamic (MHD) model, which ignores the plasma inertia. Here we consider a regularized quasi-static MHD model for force-free plasmas in tokamak disruptions and propose a mimetic finite difference (MFD) algorithm. The full geometry of an ITER-like tokamak reactor is treated, with a blanket module region, a vacuum vessel region, and the plasma region. Specifically, we develop a parallel, fully implicit, and scalable MFD solver based on PETSc and its DMStag data structure for the discretization of the five-field quasi-static perpendicular plasma dynamics model on a 3D structured mesh. The MFD spatial discretization is coupled with a fully implicit DIRK scheme. The algorithm exactly preserves the divergence-free condition of the magnetic field under the resistive Ohm's law. The preconditioner employed is a four-level fieldsplit preconditioner, which is created by combining separate preconditioners for individual fields, that calls multigrid or direct solvers for sub-blocks or exact factorization on the separate fields. The numerical results confirm the divergence-free constraint is strongly satisfied and demonstrate the performance of the fieldsplit preconditioner and overall algorithm. The simulation of ITER VDE cases over the actual plasma current diffusion time is also presented., Comment: 43 pages

The edge plasma turbulence and transport dynamics, as well as the divertor power loads during the thermal quench phase of tokamak disruptions are numerically investigated with BOUT++'s flux-driven, six-field electromagnetic turbulence model. Here a transient yet intense particle and energy sources are applied at the pedestal top to mimic the plasma power drive at the edge induced by a core thermal collapse, which flattens core temperature profile. Interesting features such as surging of divertor heat load (up to 50 times), and broadening of heat flux width (up to 4 times) on the outer divertor target plate, are observed in the simulation, in qualitative agreement with experimental observations. The dramatic changes of divertor heat load and width are due to the enhanced plasma turbulence activities inside the separatrix. Two cross-field transport mechanisms, namely the $E\times B$ turbulent convection and the stochastic parallel advection/conduction, are identified to play important roles in this process. Firstly, elevated edge pressure gradient drives instabilities and subsequent turbulence in the entire pedestal region. The enhanced turbulence not only transports particles and energy radially across the separatrix via $E\times B$ convection which causes the initial divertor heat load burst, but also induces an amplified magnetic fluctuation $\tilde{B}$. Once the magnetic fluctuation is large enough to break the magnetic flux surface, magnetic flutter effect provides an additional radial transport channel. In the late stage of our simulation, $|\tilde{B}_r/B_0|$ reaches to $10^{-4}$ level that completely breaks magnetic flux surfaces such that stochastic field-lines are directly connecting pedestal top plasma to the divertor target plates or first wall, further contributing to the divertor heat flux width broadening., Comment: 20 pages, 12 figures

We present a novel solver technique for the anisotropic heat flux equation, aimed at the high level of anisotropy seen in magnetic confinement fusion plasmas. Such problems pose two major challenges: (i) discretization accuracy and (ii) efficient implicit linear solvers. We simultaneously address each of these challenges by constructing a new finite element discretization with excellent accuracy properties, tailored to a novel solver approach based on algebraic multigrid (AMG) methods designed for advective operators. We pose the problem in a mixed formulation, introducing the directional temperature gradient as an auxiliary variable. The temperature and auxiliary fields are discretized in a scalar discontinuous Galerkin space with upwinding principles used for discretizations of advection. We demonstrate the proposed discretization's superior accuracy over other discretizations of anisotropic heat flux, achieving error $1000\times$ smaller for anisotropy ratio of $10^9$, for $closed$ $field$ $lines$. The block matrix system is reordered and solved in an approach where the two advection operators are inverted using AMG solvers based on approximate ideal restriction (AIR), which is particularly efficient for upwind discontinuous Galerkin discretizations of advection. To ensure that the advection operators are non-singular, in this paper we restrict ourselves to considering open (acyclic) magnetic field lines for the linear solvers. We demonstrate fast convergence of the proposed iterative solver in highly anisotropic regimes where other diffusion-based AMG methods fail., Comment: 29 pages, 5 figures

In tokamak disruptions where the magnetic connection length becomes comparable to or even shorter than the plasma mean-free-path, parallel transport can dominate the energy loss and the thermal quench of the core plasma goes through four phases (stages) that have distinct temperature ranges and durations. The main temperature drop occurs while the core plasma remains nearly collisionless, with the parallel electron temperature $T_{e\parallel}$ dropping in time $t$ as $T_{e\parallel}\propto t^{-2}$ and a cooling time that scales with the ion sound wave transit time over the length of the open magnetic field line. These surprising physics scalings are the result of effective suppression of parallel electron thermal conduction in an otherwise bounded collisionless plasma, which is fundamentally different from what are known to date on electron thermal conduction along the magnetic field in a nearly collisionless plasma.

For plasma velocity space instabilities driven by particle distributions significantly deviated from a Maxwellian, weak collisions can damp the instabilities by an amount that is significantly beyond the collisional rate itself. This is attributed to the dual role of collisions that tend to relax the plasma distribution toward a Maxwellian and to suppress the linearly perturbed distribution function. The former effect can dominate in cases where the unstable non-Maxwellian distribution is driven by collisionless transport on a time scale much shorter than that of collisions, and the growth rate of the ideal instability has a sensitive dependence on the distribution function. The whistler instability driven by electrostatically trapped electrons is used as an example to elucidate such a strong collisional damping effect of plasma velocity space instabilities, which is confirmed by first-principles kinetic simulations.

Benign termination of mega-ampere (MA) level runaway current has been convincingly demonstrated in recent JET and DIII-D experiments, establishing it as a leading candidate for runaway mitigation on ITER. This comes in the form of a runaway flush by parallel streaming loss along stochastic magnetic field lines formed by global magnetohydrodynamic instabilities, which are found to correlate with a low-Z injection that purges the high-Z impurities from a post-thermal-quench plasma. Here we show the competing physics that govern the post-flush reconstitution of the runaway current in a ITER-like reactor where significantly higher current is expected. The trapped ``runaways'' are found to dominate the seeding for runaway reconstitution, and the incomplete purge of high-Z impurities helps drain the seed but produces a more efficient avalanche, two of which compete to produce a 2-3~MA step in current drop before runaway reconstitution of the plasma current.

A large family of plasmas has collisional mean-free-path much longer than the non-neutral sheath width, which scales with the plasma Debye length. The plasmas, particularly the electrons, assume strong temperature anisotropy in the sheath. The temperature in the sheath flow direction ($T_{e\parallel}$) is lower and drops towards the wall as a result of the decompressional cooling by the accelerating sheath flow. The electron temperature in the transverse direction of the flow field ($T_{e\perp}$) not only is higher but also spikes up in the sheath. This abnormal behavior of $T_{e\perp}$ spike is found to be the result of a negative gradient of the parallel heat flux of transverse degrees of freedom ($q_{es}$) in the sheath. The non-zero heat flux $q_{es}$ is induced by pitch-angle scattering of electrons via either their interaction with self-excited electromagnetic waves in a nearly collisionless plasma or Coulomb collision in a collisional plasma, or both in the intermediate regime of plasma collisionality.

In a post-thermal-quench plasma, mitigated or unmitigated, the plasma power balance is mostly between collisional or Ohmic heating and plasma radiative cooling. In a plasma of atomic mixture $\{n_\alpha\}$ with $\alpha$ labeling the atomic species, the power balance sets the plasma temperature, ion charge state distribution $\{n_\alpha^i\}$ with $i$ the charge number, and through the electron temperature $T_e$ and ion charge state distribution $\{n_\alpha^i\},$ the parallel electric field $E_\parallel.$ Since the threshold electric field for runaway avalanche growth $E_{av}$ is also set by the atomic mixture, ion charge state distribution and its derived quantity, the electron density $n_e,$ the plasma power balance between Ohmic heating and radiative cooling imposes a stringent constraint on the plasma regime for avoiding and minimizing runaways when a fusion-grade tokamak plasma is rapidly terminated.

A large class of Laboratory, Space, and Astrophysical plasmas is nearly collisionless. When a localized energy or particle sink, for example, in the form of a radiative cooling spot or a black hole, is introduced into such a plasma, it can trigger a plasma thermal collapse, also known as a thermal quench in tokamak fusion. Here we show that the electron thermal conduction in such a nearly collisionless plasma follows the convective energy transport scaling in itself or in its spatial gradient, due to the constraint of ambipolar transport. As the result, a robust cooling flow aggregates mass toward the cooling spot and the thermal collapse of the surrounding plasma takes the form of four propagating fronts that originate from the radiative cooling spot, along the magnetic field line in a magnetized plasma. The slowest one, which is responsible for deep cooling, is a shock front.

As one of the fundamental graph operations, reachability queries processing has been extensively studied during the past decades. Many approaches followed the line of designing 2-hop labels to make acceleration. Considering that the index size cannot be bounded when using all nodes to construct 2-hop labels, researchers proposed to use a part of important nodes to construct 2-hop labels (partial 2-hop labels) to cover as much reachability information as possible. Then, we may achieve better query performance with limited index size and index construction time. However, partial 2-hop labels do not always perform well on different graphs. In this paper, we focus on the problem of how to efficiently compute reachability ratio, such that to help users determine whether partial 2-hop labels should be used to answer reachability queries for the given graph. Intuitively, reachability ratio denotes the ratio of the number of reachable queries that can be answered by partial 2-hop labels over the total number of reachable queries involved in the given graph. We discuss the difficulties of reachability ratio computation, and propose an incremental-partition algorithm for reachability ratio computation. We show by rich experimental results that our algorithm can efficiently get the result of reachability ratio, and show how the overall query performance is affected by different partial 2-hop labels. Based on the experimental results, we give out our findings on whether partial 2-hop labels should be used to the given graph for reachability queries processing.

The plasma exit flow speed at the sheath entrance is constrained by the Bohm criterion. The so-called Bohm speed regulates the plasma particle and power exhaust fluxes to the wall, and it is commonly deployed as a boundary condition to exclude the sheath region in quasi-neutral plasma modeling. Here the Bohm criterion analysis is performed in the intermediate plasma regime away from the previously known limiting cases of adiabatic laws and the asymptotic limit of infinitesimal Debye length in a finite-size system, using the transport equations of an anisotropic plasma. The resulting Bohm speed has explicit dependence on local plasma heat flux, temperature isotropization, and thermal force. Comparison with kinetic simulations demonstrates its accuracy over the plasma-sheath transition region in which quasineutrality is weakly perturbed and Bohm criterion applies.

Reliable plasma transport modeling for magnetic confinement fusion depends on accurately resolving the ion charge state distribution and radiative power losses of the plasma. These quantities can be obtained from solutions of a collisional-radiative (CR) model at each time step within a plasma transport simulation. However, even compact, approximate CR models can be computationally onerous to evaluate, and in-situ evaluations of these models within a coupled plasma transport code can lead to a rigid bottleneck. A way to bypass this bottleneck is to deploy artificial neural network surrogates for rapid evaluations of the necessary plasma quantities. However, one issue with training an accurate artificial neural network surrogate is the reliance on a sufficiently large and representative data set for both training and validation, which can be time-consuming to generate. In this study we further explore a data-driven active learning and training scheme to allow autonomous adaptive sampling of the problem parameter space that ensures a sufficiently large and meaningful set of training data assembled for the surrogate training. Our numerical experiments show that in order to produce a comparably accurate CR surrogate, the proposed approach requires a total number of data samples that is an order-of-magnitude smaller than a conventional approach.

The magnetohydrodynamics (MHD) equations are continuum models used in the study of a wide range of plasma physics systems, including the evolution of complex plasma dynamics in tokamak disruptions. However, efficient numerical solution methods for MHD are extremely challenging due to disparate time and length scales, strong hyperbolic phenomena, and nonlinearity. Therefore the development of scalable, implicit MHD algorithms and high-resolution adaptive mesh refinement strategies is of considerable importance. In this work, we develop a high-order stabilized finite-element algorithm for the reduced visco-resistive MHD equations based on the MFEM finite element library (mfem.org). The scheme is fully implicit, solved with the Jacobian-free Newton-Krylov (JFNK) method with a physics-based preconditioning strategy. Our preconditioning strategy is a generalization of the physics-based preconditioning methods in [Chacon, et al, JCP 2002] to adaptive, stabilized finite elements. Algebraic multigrid methods are used to invert sub-block operators to achieve scalability. A parallel adaptive mesh refinement scheme with dynamic load-balancing is implemented to efficiently resolve the multi-scale spatial features of the system. Our implementation uses the MFEM framework, which provides arbitrary-order polynomials and flexible adaptive conforming and non-conforming meshes capabilities. Results demonstrate the accuracy, efficiency, and scalability of the implicit scheme in the presence of large scale disparity. The potential of the AMR approach is demonstrated on an island coalescence problem in the high Lundquist-number regime ($\ge 10^7$) with the successful resolution of plasmoid instabilities and thin current sheets., Comment: 41 pages, 21 figures