Generating Meshes for XGC Simulations

XGC Simulations Rely on Unstructured Triangular Configuration Space Grids

  • Most gyrokinetic codes use flux-coordinates

    • Natural choice for separating motion parallel and perpendicular to the magnetic field.

    • Field-aligned coordinates enable low resolution along B for low-\(k_{\vert\vert}\) modes.

    • HOWEVER: flux-coordinates have singularity on the magnetic axis and the separatrix → problematic for whole-device simulation.

  • XGC computes particle trajectories in cylindrical coordinates to avoid coordinate singularities.

    • Good for whole-device simulation.

    • HOWEVER: Does this require high toroidal resolution even for low-\(k_{\vert\vert}\) modes? → grids are “semi-structured” → Most vertices are aligned on a finite number of flux-surfaces and placed such that they are aligned with the magnetic field.

Toroidal Decomposition

Alternative text

Radial-Poloidal Mesh

  • “Unstructured” triangular mesh, linear finite elements

  • Typical size for XGC turbulence simulations

    • \(\text{100,000 – 1,000,000 vertices per poloidal plane}\)

    • \(\text{N}_{\varphi} \sim \text{16 - 128}\)

    • \(\text{Velocity space grid:}\text{N}_{v\vert\vert} \times \text{N}_{v\perp} \sim 2000\)

    • \(10,000\; \text{particles/vertex}\)\(O(10^{10}\text{-}10^{12})\;\text{particles}\)

  • Planar mesh is divided into patches (guided by physics); each patch is assigned to one MPI rank

  • Data locality: Particle and \(f_\text{g}\) data reside on MPI rank corresponding to their position in configuration space

  • XGC_core/search.F90 → type grid_type

Alternative text

C-Mod and D3d Meshes

Typical Mesh Size

  • Mesh size in radial-poloidal plane: ~1-4 mm

  • Device size: ~m → ~ \(10^{4}\text{-}10^{6}\) vertices

  • ~104 marker particles per vertex

  • DIII-D size simulation

    • XGCa: ~30-50k vertices, ~500M particles (~1024 KNL nodes, 2.5 TFLOPS, ~1 day)

    • XGC1: ~2.5M vertices, ~25B particles (~2048 KNL nodes, 5 TFLOPS, ~2 days)

  • ITER size simulation

    • XGC1: ~18M+ vertices, ~180B particles (full-scale Titan, 27 TFLOPS, several days)

Particle Shape Function

  • For calculation of fluid moments on a (configuration/velocity/phase space) grid (e.g. density, temperature), particles are interpolated to the grid with a set of interpolation coefficients determined by shape functions Γ(z).

  • XGC uses linear shape functions in all dimensions.

Alternative text

Sketch of mesh cells in XGC1 (field-aligned) and XGCa (axisymmetric toroidal) and linear shape function in the toroidal direction

Field Solve and Electric Field Calculation

  • XGC1/es_poisson.F90 , XGC1/poisson_extra_xgc1.F90

  • Field solvers on each poloidal plane are independent linear finite element solvers.

  • After solving for the electrostatic potential, data is sent/received with MPI such that each process has data from 4 planes.

  • This data is used to compute the electric field E at \(\varphi_\text{i}\) and \(\varphi_\text{i+1}\).

  • Project along B to \(\varphi_\text{i+1/2}\).

  • Gyroaverage for ions → E_rho_ff.

  • Collect data from all planes (subroutine gather_field_info) → E_phi_ff.

Alternative text

The RPI/Simmetrix Meshing Tool for XGC

  • The RPI group (Prof. Mark Shephard) of our SciDAC project together with Simmetrix developed a meshing tool for XGC.

  • The meshing tool is available on

    • Can be run as command line tool.

    • A Graphical User Interface (GUI) is available via Simmetrix SimModeler.

  • The meshing tool produces high quality meshes and has the following key features

    • Triangular 2D mesh of a poloidal plane with field-following vertices for the user-defined region of Tokamak.

    • Flux surface curves with one element deep mesh (i.e. one element between two adjacent flux surface).

    • Unstructured meshes for undefined part of scrape-off-layer regions.

    • The code supports 0,1 and 2 X-point geometries.

  • We rely on feedback from users to continually improve this tool.

  • If you experience problems, talk to Robert Hager (

  • Meshing problems are discussed in a regular conference call with RPI and Simmetrix.

Setting up for Meshing

  • Open the portal by ssh connection to

  • Load the following modules (unload conflicting modules first, if necessary)(Updated: Dec 3rd, 2020)

    • module load intel/2019.u3 openmpi/4.0.3 hdf5-parallel/1.10.5 cmake

    • module load simmodsuite/15.0-200714 simmodeler/9.0-200714-dev

    • Sometimes, you need to reload the license server (if you experience a license error)

      • module unload rlm

      • module load rlm

    • Export OMP_NUM_THREADS=1

    • Executable: /p/epsi/rpi/install/centos7-intel2019u3-openmpi4.0.1/15.0-200714/bin/xgc_meshgen

  • Start the GUI by calling “simmodeler”.

Requirements for Meshing

  • Magnetic equilibrium in EFIT geqdsk format or XGC-eqd format (works only for cases without X-point, e.g. Cyclone).

  • Input file with meshing parameters.

    • Meshing parameters can be manually entered in the GUI or read from file.

  • Optional:

    • Density and temperature profiles as function of normalized \(\psi\).

    • Input file specifying the position (normalized \(\psi\)-values) of flux-surfaces.

    • Input file specifying the desired poloidal resolution as function of normalized \(\psi\).

  • For best results and more economical simulation, it is recommend to always prepare input files for radial and poloidal resolution.

Help on Meshing Parameters

How to Determine the Radial-Poloidal Resolution

  • First, think about your requirements for radial and poloidal resolution

    • Resolution affects the cost of your simulation (10k particles per configuration space vertex).

    • Is this a test or a production run?

    • What physics do you need to resolve (neoclassical, ITG, TEM, ETG,..)?

  • How do you want to specify mesh resolution?

    • Constant resolution everywhere?

      • Radial_uniform_meshing_unit=[1,2], intra_curve_spacing_option=1

    • Resolution is function of normalized \(\psi\)

      • Radial_uniform_meshing_unit=[-1], intra_curve_spacing_option=[-1 (XGC1),-2 (XGCa)]

  • Does the mesh have to be field-line following (XGC1, XGCa+RMP) or not (XGCa)?

  • The main factors that dictate the required resolution are

    • Density/Temperature gradient scale length

    • Ion gyroradius (electron gyroradius for ETG simulations)

  • The geqdsk/eqd files with the magnetic equilibrium and the input profiles are sufficient to calculate those characteristic lengths

  • Compute the gyroradius and the gradient scale lengths along the outer midplane

  • Radial resolution should be

    • less than \(\text{L}_{\nabla}/2\)

    • \(\text{~}\rho_\text{i}\) (thermal) if turbulence is to be resolved

  • Once the radial resolution is determined, you can compute the radial locations (normalized pol. Flux) of flux-surfaces

  • For turbulence simulations, the radial and poloidal resolution should be similar

  • For neoclassical simulations, the poloidal resolution can be larger than the radial resolution (by ~3-5x, depending on when the meshing code starts to fail due to invalid elements)

Toroidal Resolution

  • In case of field-aligned meshes, the number of poloidal planes \(\text{N}_{\varphi}^\text{mesh}\) has to be given in the input file.

  • This number refers to the full torus even if you run XGC only for a wedge.

  • Full-torus simulations (\(\text{N}_\text{w}\) =1) usually require \(\text{N}_{\varphi}^\text{XGC}\) =32-128.

  • For a wedge simulation with \(\text{N}_\text{w}\) =3, you would use \(\text{N}_{\varphi}^\text{mesh}\) =16-32.

  • The number of planes used to generate the mesh should be \(\text{N}_{\varphi}^\text{mesh}= \text{n} \times \text{N}_\text{w} \times \text{N}_{\varphi}^\text{XGC}\) , where n is a positive integer, typically between 2 and 8.

Check your Mesh

  • After the meshing tool finishes, visually check the final mesh, especially.

    • Around the magnetic axis → sometimes you need to adjust the positions of the innermost surfaces.

    • Where flux-surfaces come close to the wall curve → adjust flux-surface position if necessary.

    • Around the X-point → adjust mesh resolution around the X-point if necessary.

    • Look for bad mesh quality.