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¶
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
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.
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 , XGC_core/poisson_extra.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.
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 portal.pppl.gov.
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 (rhager@pppl.gov).
Meshing problems are discussed in a regular conference call with RPI and Simmetrix.
Environment Setup¶
Meshing tools are available on portal.pppl.gov.
For the command line tool:
source /p/xgc/Software/bin/setup_tomms_meshgen
tomms_meshgen
For the GUI:
source /p/xgc/Software/bin/setup_tomms_meshgen
simmodeler
If you experience a license error:
module unload rlm
module load rlm
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¶
Sample test cases with working input files(xgc_mesh_in) can be found in the following public directory.
The cases already contain the working input file (xgc_mesh_in). The user needs to set the environment and call the executable (As described in the section “Environment Setup” above).
A user guide and description of every control parameter can be found in the following documents.
Contact Robert Hager (rhager@pppl.gov) or Seung-Hoe Ku (sku@pppl.gov) in case of any issues.
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.