Numerical simulation of four-field extended magnetohydrodynamics in dynamically adaptive curvilinear coordinates via Newton–Krylov–Schwarz

Numerical simulations of the four-field extended magnetohydrodynamics (MHD) equations with hyper-resistivity terms present a difficult challenge because of demanding spatial resolution requirements. A time-dependent sequence of r-refinement adaptive grids obtained from solving a single Monge–Ampère (MA) equation addresses the high-resolution requirements near the x-point for numerical simulation of the magnetic reconnection problem. The MHD equations are transformed from Cartesian coordinates to solution-defined curvilinear coordinates. After the application of an implicit scheme to the time-dependent problem, the parallel Newton–Krylov–Schwarz (NKS) algorithm is used to solve the system at each time step. Convergence and accuracy studies show that the curvilinear solution requires less computational effort than a pure Cartesian treatment. This is due both to the more optimal placement of the grid points and to the improved convergence of the implicit solver, nonlinearly and linearly. The latter effect, which is significant (more than an order of magnitude in number of inner linear iterations for equivalent accuracy), does not yet seem to be widely appreciated.     

Particle-in-cell modeling of relativistic laser–plasma interaction with the adjustable-damping,direct implicit method

Implicit particle-in-cell codes offer advantages over their explicit counterparts in that they suffer weaker stability constraints on the need to resolve the higher frequency modes of the system. This feature may prove particularly valuable for modeling the interaction of high-intensity laser pulses with overcritical plasmas, in the case where the electrostatic modes in the denser regions are of negligible influence on the physical processes under study. To this goal, we have developed the new two-dimensional electromagnetic code ELIXIRS (standing for ELectromagnetic Implicit X-dimensional Iterative Relativistic Solver) based on the relativistic extension of the so-called Direct Implicit Method [D. Hewett, A.B. Langdon, Electromagnetic direct implicit plasma simulation, J. Comput. Phys. 72 (1987) 121–155]. Dissipation-free propagation of light waves into vacuum is achieved by an adjustable-damping electromagnetic solver. In the high-density case where the Debye length is not resolved, satisfactory energy conservation is ensured by the use of high-order weight factors. In this paper, we first derive the electromagnetic direct implicit method as a simplified Newton scheme. Its linear properties are then investigated through numerically solving the relation dispersions obtained for both light and plasma waves, accounting for finite space and time steps. Finally, our code is successfully benchmarked against explicit particle-in-cell simulations for two kinds of physical problems: plasma expansion into vacuum and relativistic laser–plasma interaction. In both cases, we will demonstrate the robustness of the implicit solver for crude discretizations, as well as the gains in efficiency which can be realized over standard explicit simulations.     

Influence of the inverse sheath on divertor plasma performance in tokamak edge plasma simulations

Recently, it was shown that strong electron thermionic emission from material walls could result in the formation of an “inverse sheath,” which prevents the flow of cold ions to the wall.^[1–3] Such regimes look very favourably from the point of view of plasma–material interactions at the edge of magnetic fusion devices, where the problem of the erosion of plasma-facing components under ion irradiation is one of the key issues for the development of future magnetic fusion reactors. However, it is not clear whether such regimes are compatible with edge plasma parameters and heat removal requirements in fusion reactors. To address the issue of practicality of the “inverse sheath” regime for edge tokamak plasma conditions, we perform a set of numerical simulations with 2D edge plasma transport code UEDGE^[4] for a DIII-D-like geometry and magnetic configuration. To describe both “standard” and “inverse sheath” conditions within the framework of the UEDGE code (which does not consider the sheath region per se), at the material surfaces, we apply effective boundary conditions that emulate both “standard” and “inverse sheath” regimes. We demonstrate that, for the same input parameters, spatial distributions of edge plasma parameters corresponding to detached divertor and “inverse sheath” regimes are similar, with only a few minor differences. We discuss the compatibility of “inverse sheath” regimes with core plasma parameters.     

Comparison benchmark between tokamak simulation code and TokSys for Chinese Fusion Engineering Test Reactor vertical displacement control design

Vertical displacement event(VDE) is a big challenge to the existing tokamak equipment and that being designed. As a Chinese next-step tokamak, the Chinese Fusion Engineering Test Reactor(CFETR) has to pay attention to the VDE study with full-fledged numerical codes during its conceptual design. The tokamak simulation code(TSC) is a free boundary time-dependent axisymmetric tokamak simulation code developed in PPPL, which advances the MHD equations describing the evolution of the plasma in a rectangular domain. The electromagnetic interactions between the surrounding conductor circuits and the plasma are solved self-consistently. The TokSys code is a generic modeling and simulation environment developed in GA. Its RZIP model treats the plasma as a fixed spatial distribution of currents which couple with the surrounding conductors through circuit equations. Both codes have been individually used for the VDE study on many tokamak devices, such as JT-60U, EAST, NSTX, DIII-D, and ITER. Considering the model differences, benchmark work is needed to answer whether they reproduce each other's results correctly. In this paper, the TSC and TokSys codes are used for analyzing the CFETR vertical instability passive and active controls design simultaneously. It is shown that with the same inputs, the results from these two codes conform with each other.     

Fully implicit 1D radiation hydrodynamics: Validation and verification

A fully implicit finite difference scheme has been developed to solve the hydrodynamic equations coupled with radiation transport. Solution of the time-dependent radiation transport equation is obtained using the discrete ordinates method and the energy flow into the Lagrangian meshes as a result of radiation interaction is fully accounted for. A tridiagonal matrix system is solved at each time step to determine the hydrodynamic variables implicitly. The results obtained from this fully implicit radiation hydrodynamics code in the planar geometry agrees well with the scaling law for radiation driven strong shock propagation in aluminium. For the point explosion problem the self similar solutions are compared with results for pure hydrodynamic case in spherical geometry. Results obtained when radiation interaction is also accounted agree with those of point explosion with heat conduction for lower input energies. Having, thus, benchmarked the code, self convergence of the method w.r.t. time step is studied in detail for both the planar and spherical problems. Spatial as well as temporal convergence rates are ?1 as expected from the difference forms of mass, momentum and energy conservation equations. This shows that the asymptotic convergence rate of the code is realized properly.     

Progress with the 5D full-F continuum gyrokinetic code COGENT

COGENT is an Eulerian gyrokinetic code being developed for edge plasma modelling. The code is distinguished by the use of a high-order finite-volume (conservative) discretization combined with mapped multi-block grid technology. Our recent work is focused on the development of a 5D full-F COGENT version. A numerical algorithm utilizing locally a field-aligned multi-block coordinate system is implemented to facilitate simulations of highly anisotropic microturbulence in the presence of a strong magnetic shear. In this approach, the toroidal direction is divided into blocks such that, within each block, the cells are field-aligned and a non-matching (non-conformal) grid interface is allowed at the block boundaries. Here we report on details of the numerical implementation and present preliminary results of verification studies performed for the case of the ion temperature gradient (ITG) instability in a sheared toroidal annulus geometry.     

An energy- and charge-conserving,implicit, electrostatic particle-in-cell algorithm

This paper discusses a novel fully implicit formulation for a one-dimensional electrostatic particle-in-cell (PIC) plasma simulation approach. Unlike earlier implicit electrostatic PIC approaches (which are based on a linearized Vlasov–Poisson formulation), ours is based on a nonlinearly converged Vlasov–Ampére (VA) model. By iterating particles and fields to a tight nonlinear convergence tolerance, the approach features superior stability and accuracy properties, avoiding most of the accuracy pitfalls in earlier implicit PIC implementations. In particular, the formulation is stable against temporal (Courant–Friedrichs–Lewy) and spatial (aliasing) instabilities. It is charge- and energy-conserving to numerical round-off for arbitrary implicit time steps (unlike the earlier “energy-conserving” explicit PIC formulation, which only conserves energy in the limit of arbitrarily small time steps). While momentum is not exactly conserved, errors are kept small by an adaptive particle sub-stepping orbit integrator, which is instrumental to prevent particle tunneling (a deleterious effect for long-term accuracy). The VA model is orbit-averaged along particle orbits to enforce an energy conservation theorem with particle sub-stepping. As a result, very large time steps, constrained only by the dynamical time scale of interest, are possible without accuracy loss. Algorithmically, the approach features a Jacobian-free Newton–Krylov solver. A main development in this study is the nonlinear elimination of the new-time particle variables (positions and velocities). Such nonlinear elimination, which we term particle enslavement, results in a nonlinear formulation with memory requirements comparable to those of a fluid computation, and affords us substantial freedom in regards to the particle orbit integrator. Numerical examples are presented that demonstrate the advertised properties of the scheme. In particular, long-time ion acoustic wave simulations show that numerical accuracy does not degrade even with very large implicit time steps, and that significant CPU gains are possible.     

The non-equilibrium Marshak wave problem

An analytic solution to a particular Marshak wave problem is given. The radiative transfer model used is the gray, non-equilibrium diffusion approximation which allows the radiation and material fields to be out of equilibrium. This solution should be useful as a reference problem for validating time-dependent radiative transfer computer codes, as well as investigating the convergence, as a function of space and time step size, for such codes. The coupling of the radiation field to the material field in a multigroup code, a difficult numerical problem, can also be tested against this solution. Typical numerical results are given for surface quantities, integral quantities, and the distribution of radiative energy and material temperature as a function of space and time.     

光码分多址系统中光学相关接收机判决阈值的分析

介绍了光码分多址系统中常用地址码（一维扩时码、二维码和三维码）的特点，并对它们各自的互相关均值和方差进行了理论分析。基于非相干光码分多址系统中光学相关接收机的基本原理，结合不同的用户地址码，对系统误码率性能进行了分析，得到了接收机最佳判决阈值与地址码基本特性参数和系统同时用户数间的关系。最后，给出了数值仿真结果。结果表明，对于采用特定地址码的光码分多址系统，只有选择合适的接收机判决阈值，系统的误码率性能才能达到最佳。研究结果对光码分多址系统中接收机判决阈值的选取具有一定的参考作用。     

Curvilinear finite-volume schemes using high-order compact interpolation

During the last years, the need of high fidelity simulations on complex geometries for aeroacoustics predictions has grown. Most of high fidelity numerical schemes, in terms of low dissipative and low dispersive effects, lie on finite-difference (FD) approach. But for industrial applications, FD schemes are less robust compared to finite-volume (FV) ones. Thus the present study focuses on the development of a new compact FV scheme for two- and three-dimensional applications.The proposed schemes are formulated in the physical space and not in the computational space as it is the case in most of the known works. Therefore, they are more appropriate for general grids. They are based on compact interpolation to approximate interface-averaged field values using known cell-averaged values. For each interface, the interpolation coefficients are determined by matching Taylor series expansions around the interface center. Two types of schemes can be distinguished. The first one uses only the curvilinear abscissa along a mesh line to derive a sixth-order compact interpolation formulae while the second, more general, uses coordinates in a spatial three-dimensional frame well chosen. This latter is formally sixth-order accurate in a preferred direction almost orthogonal to the interface and at most fourth-order accurate in transversal directions.For non-linear problems, different approaches can be used to keep the high-order scheme. However, in the present paper, a MUSCL-like formulation was sufficient to address the presented test cases.All schemes have been modified to treat multiblock and periodic interfaces in such a way that high-order accuracy, stability, good spectral resolution, conservativeness and low computational costs are guaranteed. This is a first step to insure good scalability of the schemes although parallel performances issues are not addressed. As high frequency waves, badly resolved, could be amplified and then destabilize the scheme, compact filtering operators have been used.Numerous test cases as the linear convection of a Gaussian wave, the convection of a Lamb–Oseen vortex and the diffraction of an acoustic wave on a plane have been realized to validate the schemes. The most efficient schemes are shown to be at least fifth-order accurate on linear and non-linear convection problems. They are also less dissipative and less dispersive on non-uniform curvilinear grids than schemes using implicit interpolation with constant coefficients of the same order on uniform cartesian grids.     

Operator Splitting Implicit Integration Factor Methods for Stiff Reaction-Diffusion-Advection Systems

For reaction-diffusion-advection equations, the stiffness from the reaction and diffusion terms often requires very restricted time step size, while the nonlinear advection term may lead to a sharp gradient in localized spatial regions. It is challenging to design numerical methods that can efficiently handle both difficulties. For reaction-diffusion systems with both stiff reaction and diffusion terms, implicit integration factor (IIF) method and its higher dimensional analog compact IIF (cIIF) serve as an efficient class of time-stepping methods, and their second order version is linearly unconditionally stable. For nonlinear hyperbolic equations, weighted essentially non-oscillatory (WENO) methods are a class of schemes with a uniformly high-order of accuracy in smooth regions of the solution, which can also resolve the sharp gradient in an accurate and essentially non-oscillatory fashion. In this paper, we couple IIF/cIIF with WENO methods using the operator splitting approach to solve reaction-diffusion-advection equations. In particular, we apply the IIF/cIIF method to the stiff reaction and diffusion terms and the WENO method to the advection term in two different splitting sequences. Calculation of local truncation error and direct numerical simulations for both splitting approaches show the second order accuracy of the splitting method, and linear stability analysis and direct comparison with other approaches reveals excellent efficiency and stability properties. Applications of the splitting approach to two biological systems demonstrate that the overall method is accurate and efficient, and the splitting sequence consisting of two reaction-diffusion steps is more desirable than the one consisting of two advection steps, because CWC exhibits better accuracy and stability.     

Periodic Filtering as a subgrid-scale model for LES of laminar separation bubble flows

Laminar separation bubbles develop over many blades and airfoils at moderate angles of attack and Reynolds numbers ranging from 10⁴ to 10⁵. More accurate simulation tools are necessary to enable higher fidelity design optimisation for these airfoils and blades as well as to test flow control schemes. Following previous investigators, an equivalent problem is formulated by imposing suitable boundary conditions for flow over a flat plate which allows to use a high accuracy spectral solver. Large eddy simulation (LES) of such a flow were performed at drastically reduced resolution to assess the accuracy of several LES modelling approaches: the explicit dynamic Smagorinsky model, implicit LES, and the truncated Navier–Stokes approach (TNS). To mimic dissipation that occurs in implicit LES, the solution on a coarse mesh is filtered at every time step and two different filter strengths are used. In the TNS approach, the solution is filtered periodically, every few hundred time steps. The performance of each approach is evaluated against benchmark direct numerical simulation (DNS) data focusing on pressure and skin friction distributions, which are critical to airfoil designers. TNS results confirm that periodic filtering can act as an apt substitute for explicit subgrid-scale models, whereas filtering at every time step demonstrates the dependence of implicit LES on details of numerics.     

求解Burgers方程的一种交替分段隐式算法

首先提出一个新的求解Burgers方程的差分格式,然后在此差分格式的基础上构造了便于并行计算的交替分段隐格式,并作了线性化稳定性分析.数值结果表明,本方法具有较高的精度,尤适于扩散项系数较小时的计算,且有效避免了数值结果的非物理振荡.     

XTOR-2F: A fully implicit Newton–Krylov solver applied to nonlinear 3D extended MHD in tokamaks

XTOR-2F solves a set of extended magnetohydrodynamic (MHD) equations in toroidal tokamak geometry. In the original XTOR code, the time stepping is handled by a semi-implicit method 1, 2 and 3. Moderate changes were necessary to transform it into a fully implicit one using the NITSOL library with Newton–Krylov methods of solution for nonlinear system of equations [4]. After addressing the sensitive issue of preconditioning and time step tuning, the performances of the semi-implicit and the implicit methods are compared for the nonlinear simulation of an internal kink mode test case within the framework of resistive MHD including anisotropic thermal transport. A convergence study comparing the semi-implicit and the implicit schemes is presented. Our main conclusion is that on one hand the Newton–Krylov implicit method, when applied to basic one fluid MHD is more computationally costly than the semi-implicit one by a factor 3 for a given numerical accuracy. But on the other hand, the implicit method allows to address challenging issues beyond MHD. By testing the Newton–Krylov method with diamagnetic modifications on the dynamics of the internal kink, some numerical issues, to be addressed further, are emphasized.     

A finite volume implicit time integration method for solving the equations of ideal magnetohydrodynamics for the hyperbolic divergence cleaning approach

A finite volume numerical technique is proposed to solve the compressible ideal MHD equations for steady and unsteady problems based on a quasi-Newton implicit time integration strategy. The solenoidal constraint is handled by a hyperbolic divergence cleaning approach allowing its satisfaction up to machine accuracy. The conservation of the magnetic flux is computed in a consistent way using the numerical flux of the finite volume discretization. For the unsteady problem, the time accuracy is obtained by a Newton subiteration at each physical timestep thereby converging the solenoidal constraint to steady state. We perform extensive numerical experiments to validate and demonstrate the capabilities of the proposed numerical technique.     

A mixed implicit–explicit finite difference scheme for heat transport in magnetised plasmas

An explicit/implicit domain decomposition method is applied to the time-dependent heat-conduction problem in a 2-d, strongly anisotropic medium (a magnetised plasma), using a formulation of the spatial derivatives which avoids the pollution of perpendicular by parallel heat fluxes. The time-stepping at the sub-domain boundaries is done using a DuFort–Frankel scheme, which leads to a time step limit given not by instabilities, but by the damping rate of numerical oscillations driven by inconsistencies in the formulation of initial conditions or the temporal variations in the true physical solution. These limitations can be minimized, however, by aligning the subdomain boundaries as much as possible with magnetic flux surfaces. The time step limit depends on the ratio of the implicit grid spacing to the distance between subdomain boundaries (DuFort–Frankel lines in 2-d, surfaces in 3-d).     

热辐射输运问题的高效蒙特卡罗模拟方法

惯性约束聚变研究中,热辐射光子在介质中的输运以及热辐射光子与介质的相互作用是重要研究课题,蒙特卡罗方法是该类问题的重要研究手段之一.隐式蒙特卡罗方法虽然能正确地模拟热辐射在介质中的输运过程,但当模拟重介质(材料的吸收系数大)问题时,该方法花费的计算时间将变得很长,导致模拟效率很低.本文以离散扩散蒙特卡罗方法为基础,开发了"离散扩散蒙特卡罗方法辐射输运模拟程序",可以较好地解决重介质区的计算效率问题,但是离散扩散蒙卡罗方法在模拟轻介质区时精度不够高.辐射输运问题中通常既有轻介质也有重介质,为了能同时解决蒙特卡罗方法模拟的效率和精度问题,本文研究了离散扩散蒙特卡罗方法与隐式蒙特卡罗方法相结合的模拟方法,并提出了新的扩散区与输运区界面处理方法,研制了混合蒙特卡罗方法的辐射输运模拟程序.典型辐射输运问题模拟显示:在模拟重介质问题时,该程序能大幅缩短模拟时间,且能取得与隐式蒙特卡罗方法一致的结果;在模拟轻重介质均存在的问题时,与隐式蒙特卡罗方法相比,混合蒙特卡罗方法的模拟精度与其相当且计算效率同样能够得到显著提升.     

EAST低杂波电流驱动的数值模拟及分析

介绍了波的射线轨迹程序C3PO和求解三维Fokker-Planck方程程序LUKE的基本计算过程及其物理机制.该程序采用的是反弹平均方法,考虑了相对论效应和捕获粒子效应.利用此程序模拟并分析了EAST偏滤器位形下等离子体电流对低杂波传播以及驱动效率的影响.计算结果表明,随着电流的增大,低杂波射线轨迹在极向的旋转加快,功...     

A second-order implicit particle mover with adjustable damping

A new algorithm for the calculation of particle trajectories is introduced. The algorithm combines second-order accuracy in the real frequency with third-order user-adjustable attenuation. It requires little storage of data from previous time levels. The method was designed for use in implicit particle-in-cell plasma simulation codes, and this application is treated in detail. It may also prove useful in other applications where one seeks to preserve the accuracy of low-frequency oscillations while rapidly damping under-resolved high frequency motions, e.g., solution of the field equations in electromagnetic particle codes. An explicit variation, wherein future quantities are obtained by extrapolation, can provide attenuation but not large-timestep stability.     

Numerical modeling of thermal force in a plasma for test-ion transport simulation based on Monte Carlo Binary Collision Model

We present a new numerical model of the thermal force in a plasma, based on the Monte Carlo Binary Collision Model (BCM) [T. Takizuka, H. Abe, J. Comput.Phys. 25 (1977) 205]. This model can be applied for the transport simulation of test ions. The model consists of two steps: (i) choosing a background plasma ion velocity from a distorted Maxwell distribution under the temperature gradient, and (ii) calculating a Coulomb collision between a test particle and the above chosen ion by using the BCM. For the step (i), we developed a velocity sampling method from a distorted Maxwellian, which enables the BCM to bring the thermal force on a test particle in the step (ii).A systematic series of simulations has been performed under various conditions to examine the model. The results of these simulations have been compared with the theoretical values, and it is shown that our model simulates the thermal force correctly for important characteristic features; dependences on the temperature gradient, the test particle velocity, and the background plasma density.