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-:math:`k_{\vert\vert}` modes.
- **HOWEVER:** flux-coordinates have singularity on the magnetic axis and the separatrix |rarr| 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-:math:`k_{\vert\vert}` modes? |rarr| grids are **“semi-structured”** |rarr| Most vertices are aligned on a finite number of flux-surfaces and placed such that they are aligned with the magnetic field.
Toroidal Decomposition
----------------------
.. image:: ./img/tor_dec.png
:width: 600
:align: center
:alt: Alternative text
Radial-Poloidal Mesh
--------------------
* “Unstructured” triangular mesh, linear finite elements
* Typical size for XGC turbulence simulations
* :math:`\text{100,000 – 1,000,000 vertices per poloidal plane}`
* :math:`\text{N}_{\varphi} \sim \text{16 - 128}`
* :math:`\text{Velocity space grid:}\text{N}_{v\vert\vert} \times \text{N}_{v\perp} \sim 2000`
* :math:`10,000\; \text{particles/vertex}` |rarr| :math:`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 :math:`f_\text{g}` data reside on MPI rank corresponding to their position in configuration space
* XGC_core/search.F90 |rarr| type grid_type
.. figure:: ./img/mesh_image.png
:width: 600
:align: center
:alt: Alternative text
C-Mod and D3d Meshes
Typical Mesh Size
-----------------
* Mesh size in radial-poloidal plane: ~1-4 mm
* Device size: ~m |rarr| ~ :math:`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.
.. figure:: ./img/shape_function.png
:width: 400
:align: center
:alt: 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 :math:`\varphi_\text{i}` and :math:`\varphi_\text{i+1}`.
* Project along **B** to :math:`\varphi_\text{i+1/2}`.
* Gyroaverage for ions |rarr| E_rho_ff.
* Collect data from all planes (subroutine gather_field_info) |rarr| E_phi_ff.
.. figure:: ./img/field_solve.png
:width: 400
:align: center
:alt: 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 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.
Setting up for Meshing
----------------------
* Open the portal by ssh connection to portal.pppl.gov.
* 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 :math:`\psi`.
* Input file specifying the position (normalized :math:`\psi`-values) of flux-surfaces.
* Input file specifying the desired poloidal resolution as function of normalized :math:`\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 "Setting up for Meshing" 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 :math:`\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 :math:`\text{L}_{\nabla}/2`
* :math:`\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 :math:`\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 (:math:`\text{N}_\text{w}` =1) usually require :math:`\text{N}_{\varphi}^\text{XGC}` =32-128.
* For a wedge simulation with :math:`\text{N}_\text{w}` =3, you would use :math:`\text{N}_{\varphi}^\text{mesh}` =16-32.
* The number of planes used to generate the mesh should be :math:`\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 |rarr| sometimes you need to adjust the positions of the innermost surfaces.
* Where flux-surfaces come close to the wall curve |rarr| adjust flux-surface position if necessary.
* Around the X-point |rarr| adjust mesh resolution around the X-point if necessary.
* Look for bad mesh quality.
.. |xrArr| unicode:: U+27F9 .. Long Right Double Arrow
.. |rarr| unicode:: U+2192 .. Short Right Arrow