Generating Meshes for XGC Simulations ===================================== XGC Simulations Rely on Unstructured Triangular Configuration Space Meshes -------------------------------------------------------------------------- Most gyrokinetic codes use flux-coordinates, i.e., coordinates in which the field line pitch is constant on a flux-surface, because they are the natural choice for separating motion and spatial scales parallel and perpendicular to the magnetic field. Flux-coordinates enable simulations with low resolution along the magnetic field **B** for perturbations with low parallel wavenumber :math:`k_\parallel=\frac{\mathbf{B}}{|B|}\cdot\mathbf{k}`. **However**, the metric of flux-coordinates has singularities on the magnetic axis and the separatrix. This is a problem for whole-volume simulations. XGC computes particle trajectories in right-handed cylindrical coordinates :math:`(R,\varphi,Z)` to avoid coordinate singularities. This approach is good for whole-volume simulations from the magnetic axis to the material wall. In cylindrical coordinates, general unstructured meshes would require very high resolution in all three dimensions. But for simulation of tokamak and stellarator turbulence, which is characterized by low :math:`k_\parallel` and high perpendicular wavenumbers (:math:`k_\perp \rho_s \sim 1`), it is still possible to use reduced toroidal resolution similar to flux-coordinates. This is achieved by aligning the mesh vertices on a finite number of magnetic flux-surfaces placed such that they are approximately aligned with the magnetic field. See the following sections for more details. **References:** * `M. F. Adams, et al., J. Phys.: Conf. Ser. 180, 012036 (2009) `_ * `F. Zhang et al., Engineering with Computers 32, pp. 285-293 (2015) `_ * `U. Riaz et al., Comp. Phys. Comm. 295, 108982 (2023) `_ Configuration Space Discretization ---------------------------------- Toroidal Discretization ^^^^^^^^^^^^^^^^^^^^^^^ In the toroidal (:math:`\varphi`) direction, configuration space is discretized on a uniform toroidal angle grid of size :math:`N_\varphi`. That is, the topological torus is divided into :math:`N_\varphi` sections. .. image:: ./img/tor_dec.png :width: 600 :align: center :alt: Alternative text Radial-Poloidal Plane ^^^^^^^^^^^^^^^^^^^^^ On the radial-poloidal (:math:`\varphi=\text{const.}`) plane, configuration space is discretized using unstructured triangle meshes with linear finite elements. That is, there is an unstructured mesh on each boundary between adjacent toroidal sections. In axisymmetric magnetic confinement devices (i.e., tokamaks), those meshes are identical on each radial-poloidal plane. In stellarator geometry, each plane has a different mesh. However, the planar meshes are not entirely unstructured. Most vertices are placed on a finite number of flux-surfaces in an approximately field-aligned way. This means that a magnetic field-line starting from a vertex on plane :math:`\varphi_i=2\pi i/N_\varphi` intersects plane :math:`\varphi_{i+1}` at or very close to another vertex. The approximate field-alignment of the XGC meshes enables simulation of strongly anisotropic turbulence (low :math:`k_\parallel`, high :math:`k_\perp`) with low toroidal resolution similar to gyrokinetic codes using flux-coordinates, but without the associated coordinate singularities [see `R. Hager et al., Physics of Plasmas 29, 112308 (2022) `_]. .. figure:: ./img/field_aligned_mesh.png :width: 500 :align: center :alt: Alternative text Illustration of approximate field alignment Typical Mesh Size ^^^^^^^^^^^^^^^^^ * Mesh size in radial-poloidal plane: ~1-4 mm * Device size: :math:`\sim 1 \text{m} \rightarrow N_n \sim 10^{4}\text{-}10^{6}` vertices per radial-poloidal plane * :math:`\sim 10^4` marker particles per vertex * DIII-D size simulation * Axisymmetric (:math:`N_\varphi=1`): :math:`\sim 30,000-50,000` vertices, :math:`\sim 500\cdot 10^6` particles/species * Non-axisymmetric (:math:`N_\varphi \sim 32`): :math:`\sim 5 \cdot 10^6` vertices, :math:`\sim 50 \cdot 10^9` particles/species * ITER size simulation * Non-axisymmetric (:math:`N_\varphi=32`): :math:`\sim 32 \cdot 10^6` vertices, :math:`>320\cdot 10^9` particles/species .. figure:: ./img/mesh_image.png :width: 600 :align: center :alt: Alternative text C-Mod and D3D Meshes The Tokamak Modeling and Meshing Software (TOMMS) ------------------------------------------------- Meshes for XGC simulations are generated using the Tokamak Modeling and Meshing Software (`TOMMS `_) developed at the Scientific Computation Research Center (SCOREC) at Rensselaer Polytechnic Institute (RPI) in collaboration with `Simmetrix `_. The meshing tool is installed on the PPPL cluster Portal (``portal.pppl.gov``) and actively maintained. Installation on other resources is possible but requires a license for the Simmetrix products SimModSuite and SimModeler. TOMMS can be run as a command line tool or with a Graphical User Interface (GUI) available in Simmetrix SimModeler. TONMS 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 a tokamak. * Flux-surfaces with one element deep mesh (i.e. one element between two adjacent flux surface). * Unstructured meshes for undefined parts of the scrape-off layer and private flux regions. * The code supports tokamak geometries with multiple X-points. We rely on feedback from users to continually improve this tool. If you experience problems, contact Robert Hager (rhager@pppl.gov) or Seung-Hoe Ku (sku@pppl.gov) and/or open a Github issue at ``_. Environment Setup on Portal ^^^^^^^^^^^^^^^^^^^^^^^^^^^ TOMMS is available on ``portal.pppl.gov``. To set up your shell environment for running TOMMS, just run .. code-block:: bash > source /p/xgc/Software/bin/setup_tomms_meshgen If you maintain your own installation of TOMMS, remember to set the number of OpenMP threads to 1: .. code-block:: bash > export OMP_NUM_THREADS=1 To run the command-line version of TOMMS: .. code-block:: bash > tomms_meshgen For the GUI: .. code-block:: bash > simmodeler In SimModeler, select the "Fusion" tab. TOMMS Documentation and Support ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ * A user guide and description of every control parameter can be found on the following public pages: * `User Guide `_ * `Input Parameters `_ * Sample test cases with working input files (``xgc_mesh_in``) can be found in the following public directory: * `Meshing Examples `_ * The cases already contain a working input file (``xgc_mesh_in``). The user needs to set up the environment and call the TOMMS executable as described above. A Step-by-step Guide to XGC Tokamak Mesh Generation --------------------------------------------------- Collect the Required Input Data ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Certain input files are mandatory, others are optional. At the very least, the following two input files are required: * Magnetic equilibrium in * EFIT geqdsk format, or * XGC-eqd format: generated with ``efit2eqd`` included in the XGC repository; works only for cases without X-point, e.g. Cyclone). * The resolution of the EFIT geqdsk/eqd files should be at least 256x256. * Input file with meshing parameters (default name: ``xgc_mesh_in``). * In the SimModeler GUI, meshing parameters can be manually entered or read from file. Optional input files are: * Ion temperature profiles as function of the normalized poloidal magnetic flux :math:`\psi_N=\frac{\psi-\psi_{axis}}{\psi_{bd}-\psi_{axis}}`, where :math:`\psi_{bd}` is the magnetic flux on the boundary surface (separatrix in divertor geometry) * This is used for setting the radial (i.e., inter flux-surface) mesh size based on the local ion gyroradius. Format: .. code-block:: bash [Number of entries N] psi_1 T_i1 [in eV] psi_2 T_i2 [...] psi_N T_iN * Input file specifying the position (in terms of :math:`\psi_N`) of the desired flux-surfaces: .. code-block:: bash [Number of surfaces N] psi_1 psi_2 psi_3 [...] psi_N * A separate file specifying the flux-surfaces in two private flux regions can be used, if desired: .. code-block:: bash [Number of surfaces N] psi_1a psi_1b psi_2a psi_2b [...] psi_Na psi_Nb * Input file specifying the desired poloidal resolution in meter (:math:`m`) as function of :math:`\psi_N`. .. code-block:: bash [Number of entries N] psi_1 dx_1 psi_2 dx_2 [...] psi_N dx_N For best results and more economical simulation, it is recommend to always prepare input files for radial and poloidal resolution. In this case, it is useful to utilize the different (ion/impurity/electron) density, temperature and toroidal rotation profiles of the plasma to be simulated with XGC because they help with choosing the appropriate mesh resolution. This is a minimal example for a TOMMS input file: .. code-block:: fortran ! License file for Simmetrix software ! For Portal: /usr/pppl/Simmetrix/simmodsuite.lic SimLic [LICENSE FILE] ! ! Input EFIT file containing psi grid and wall boundary. eqdskFile [PATH TO EFIT geqdsk FILE] ! ! If you want to use a custom wall curve instead of the ! one included in the EFIT geqdsk file --> !limiterFile [PATH TO WALL CURVE FILE] ! ! If set to 1, this reverses the sign of Psi. ! Set to 1 if your geqdsk file has a maximum of psi_N ! on the magnetic axis. reverse_psi 0 ! ! Set to one to use a randomized starting position ! around the outer midplane for each flux-surface random_start 1 ! ! Set to 1 for advanced log output adv_log 1 ! ! Set to 1 if you want a mesh for only an annulus ! (excluding the core around the magnetic axis) finite_domain_on 0 ! ! Choose a way of uniformity in radial direction (1: psi, 2: sqrt(psi)) ! The recommended mode is -1 --> Read position of flux-surfaces ! from an input file --> see "inter_curve_spacing_file" ! and "inter_curve_spacing_private_file" (private flux regions) radial_uniform_meshing_unit -1 ! ! This is the desired spacing between psi curves, in meters ! The first three are optional with radial_uniform_meshing_unit=-1 !inter_curve_spacing 0.002 !inter_curve_spacing_max 0.03 !inter_curve_spacing_min 0.00001 inter_curve_spacing_axis_forcing 0.004 inter_curve_spacing_file [SURFACE FILE] ! ! Private region inter-curve spacing ! (optional but recommended) inter_curve_spacing_private_on 1 inter_curve_spacing_private_file [SURFACE FILE PRIVATE REGION] ! ! This is the desired spacing between points on a psi curve, in meters ! The options -1 and -2 both read the poloidal target spacing from a file. ! -1 will generate a field-aligned mesh for turbulence simulations ! -2 will generate a non-field-aligned mesh for axisymmetric simulations intra_curve_spacing_option -2 intra_curve_spacing_file [POLOIDAL SPACING FILE] ! ! For intra_curve_spacing_options=-2, the poloidal spacing is adjusted ! for flux-expansion, i.e., larger poloidal spacing where |grad(psi_N)| ! is smaller, smaller spacing where it's higher. intra_curve_spacing_propf_min 0.1 intra_curve_spacing_propf_max 4.0 ! For approximately constant intra_curve_spacing !intra_curve_spacing 0.001 ! ! Curve with psi value larger than it will be dropped drop_open_psi 0.090808868 drop_open_psi_hfs 0.050852966 ! ! Set to 1 if you want flux-surface/field-aligned mesh ! in the private flux-region fluxC_in_privateR_on 1 ! Set to 1 for field-aligned meshes with high resolution ! (i.e. for turbulence simulations) transitionFacesOn 0 ! ! Enforce minimum distance of vertices from the wall curve intra_curve_minlength_lastedge 0.001 ! ! Defines the mesh size at the boundary wall. meshSizeWall 2.5 ! ! Number of 2D planes around the torus, used in spacing points on a curve. ! Set >1 for turbulence simulations num_phi_planes 1 ! ! The following tolerances are useful to fine-tune the poloidal mesh simulation ! They basically determine the threshold (as a fraction/multiple of the desired ! poloidal resolution) at which the field-alignment is given up to prevent the mesh ! to become too fine. ! If possible, point spacing on a curve will fall within this fraction of the desired value. spacing_tolerance_optimal 0.7 ! ! Point spacing on a curve will always be within this fraction of the desired value. ! The acceptable range is calculated as 1/(1+tol) to (1+tol) so a value of 1.0 gives a range of .5 to 2.0 times the desired spacing. spacing_tolerance_absolute 0.7 ! !The output node and element files nodeFile [OUTPUT FILE NAME].node elmFile [OUTPUT FILE NAME].ele addinfoflxFile [OUTPUT FILE NAME].flx.aif More examples can be found at ``_. Prepare Your Plasma Profile Data ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A common practical issue in the preparation of XGC meshes and simulations based on experimental data is that the plasma profiles (densities :math:`n_s(\psi_N)`, temperatures :math:`T_s(\psi_N)`, toroidal rotation :math:`\omega_T(\psi_N)`) are provided only up to the separatrix, :math:`\psi_N \leq 1`. Extending the plasma profiles into the scrape-off layer is an important step of the meshing workflow because the extended profile data is used in setting the mesh resolution as well as in the actual XGC simulation later. A simple but effective approach is to fit an exponential decay to the experimental data: .. math:: y(\psi_N)=y_0 \exp[-\gamma (\psi_N-\psi_c)]+y_1. From the experimental data, we know :math:`y(\psi_c)` and :math:`y^\prime(\psi_c)`. Hence, :math:`y_0 = y(\psi_c)-y_1`, and :math:`\gamma=-\frac{y^\prime(\psi_c)}{y(\psi_c)-y_1}`, where :math:`y_1 \rho_i`. (:math:`\rho_s` is the main ion gyroradius.) * For turbulence simulations (total-f or delta-f), typically :math:`\Delta r \sim \rho_i`. (Meshes with global electron scale resolution :math:`\Delta r \sim \rho_e` are too large to be practical.) * The area of interest: * For simulations of core plasma, the edge can have lower resolution. * For simulations of edge plasma, the core can have lower resolution. * For the typical application to edge plasma turbulence, using lower resolution in the core, e.g., at :math:`\psi_N<0.7`, can lead to a reduction of the number of vertices by up to 50%. Hence, setting the radial resolution via input file ``radial_uniform_meshing_unit=-1`` is the recommended default option. * If you want to solve the for the electrostatic potential in the private flux region (``fluxC_in_privateR_on=1``), it is recommended to also set ``inter_curve_spacing_private_on=1`` to specify the flux-surfaces in the private region via input file. Otherwise, the flux-surfaces in the private regions will be generated at the same values of :math:`\psi_N` as in the confined region and the scrape-off layer. But in practice, the flux-surfaces in the private region should be denser in :math:`\psi_N` than in the scrape-off layer to adequately resolve the potential gradients perpendicular to the flux-surfaces. Too low resolution in the private regions often manifests in numerical instabilities on relatively long time scales. See the `TOMMS documentation `_ for more details. Choose the Poloidal Resolution ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The required poloidal resolution is determined by * Whether the mesh is for an axisymmetric (neoclassical) or non-axisymmetric (turbulence) simulation. * **Axisymmetric:** use ``intra_curve_spacing_option=-2`` to generate a non-field-aligned mesh. Set the base poloidal mesh size :math:`\Delta l_\theta` to about 1-3 times the radial resolution. The actual poloidal resolution is adjusted proportional to :math:`\frac{|\nabla(\psi)|(\theta)}{|\nabla(\psi)|(\theta=0)}`, where :math:`\theta=0` corresponds to the outer midplane. The upper and lower bounds of this scaling are set with the input parameters ``intra_curve_spacing_propf_min`` and ``intra_curve_spacing_propf_max``. The latter may need to be decreased if the meshing process fails. * **Non-axisymmetric:** use ``intra_curve_spacing_option=-1`` to generate a field-aligned mesh. Set :math:`\Delta l_\theta \approx \Delta r`. This is generally sufficient for electrostatic simulations. Electromagnetic simulations may require more fine-tuning as explained below. For some applications, it may make sense to use anisotropic radial-poloidal resolution. * The area of interest: use isotropic radial-poloidal resolution in the main area of interest. Elsewhere, the poloidal mesh size can be larger than the radial mesh size to limit the mesh size (i.e., the cost of the simulation). Choose the Toroidal Resolution ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ For axisymmetric simulations, the choice is simple. (There is none.) Just set ``num_phi_planes=1``. For non-axisymmetric simulations, the choice is somewhat tricky because the poloidal and (effective) toroidal resolutions are correlated. To explain this, we have to dive into the details of how TOMMS generates field aligned meshes. Generating a Field-aligned Mesh """"""""""""""""""""""""""""""" A field aligned mesh for axisymmetric magnetic fields (tokamaks) can be generated as follows: * For each flux-surface, one picks a starting position, which will be the first vertex on that surface in the mesh. * From the starting position, trace the magnetic field around the torus and record the :math:`(R,\varphi,Z)` position in intervals of :math:`\Delta \varphi=2\pi/\mathtt{num\_phi\_planes}` until one circuit around the magnetic axis is (almost) complete. The resulting set of coordinates represents a field-aligned discretization of a single field-line. * One could repeat this process for more starting positions on each flux-surface to generate a set of ``num_phi_planes`` planar, unstructured meshes that are field-aligned. But since the magnetic field is axisymmetric in tokamaks, it is sufficient to replicate the :math:`(R,Z)` pairs obtained from tracing a single field line on all ``num_phi_planes`` planes. The resulting mesh is field-aligned and identical on each radial-poloidal plane. From this discussion, it is obvious that in a strictly field-aligned mesh the number of radial-poloidal planes (i.e., the toroidal resolution) dictates the poloidal resolution, and that the input parameters for the poloidal and toroidal mesh resolutions can be conflicting. It is also obvious that control over the poloidal mesh size is limited by the field line pitch, roughly :math:`\Delta l_\theta \propto B_P/B_T`. When the poloidal magnetic field decreases, the poloidal mesh size becomes smaller (higher resolution). This is problematic in the vicinity of critical points with :math:`B_P=0` (magnetic axis, X-points). How TOMMS Resolves the Conflict between Poloidal and Toroidal Resolution """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" In practice, the notion that XGC meshes are approximately field aligned means that the poloidal mesh size constraint (``intra_curve_spacing``) has higher priority in the mesh generation than the field-alignment constraint. The targeted :math:`\Delta l_\theta^{(target)}` (``intra_curve_spacing``) together with a tolerance parameter :math:`\sigma` (``spacing_tolerance_absolute``) ensures that the effective poloidal mesh size satisfies the relation .. math:: \frac{\Delta l_\theta^{(target)}}{1+\sigma} \leq \Delta l_\theta^{(actual)} \leq (1+\sigma) \Delta l_\theta^{(target)} Whenever field-aligned meshing would result in a value of :math:`\Delta l_\theta^{(actual)}` outside those bounds, the toroidal resolution used in the mesh generation is (locally) increased or decreased by a factor of 2. Effectively, the toroidal resolution becomes a function of :math:`(R,Z)`. Relation to Flux-coordinates """""""""""""""""""""""""""" The field-aligned meshing process explained above is closely related to flux-coordinates. Following `R. Hager et al., Phys. Plasmas 29, 112308 (2022) `_, the straight field-line angle :math:`\theta^\ast` and the safety factor :math:`q` can be defined as .. math:: \theta^\ast = \frac{1}{q(\psi_N)} \int_0^{l_\theta} \frac{B_T}{R B_P} \mathrm{d}l_{\theta}^\prime,\\ q(\psi_N) = \frac{1}{2 \pi} \int_0^{l_{\theta,max}} \frac{B_T}{R B_P} \mathrm{d}l_{\theta}^\prime, where :math:`l_\theta` is the arc length along a flux-surface in the poloidal direction, and integration from :math:`0` to :math:`l_{\theta,max}` is over the whole flux-surface. (Refer to the paper cited above for more details of this definition.) The coordinate set :math:`(\psi_N,\varphi,\theta^\ast)` forms a right-handed flux-coordinate system. Hence, when tracing a magnetic field line from :math:`\varphi_i` to :math:`\varphi_{i+1}`, the corresponding change of the straight field-line angle :math:`\theta^\ast` is .. math:: \Delta \theta^\ast = \Delta \varphi/q(\psi_N) = \frac{2\pi}{N_\varphi q(\psi_N)}. Thus, the number of radial poloidal planes :math:`N_\varphi` (``num_phi_planes``) used for field-aligned meshing determines the poloidal mesh spacing in :math:`\theta^\ast`. Also, a strictly field-aligned mesh for :math:`N_\varphi` planes has uniform mesh spacing :math:`\Delta \theta^\ast`. The meshing process described above, in which :math:`N_\varphi` is adjusted dynamically implies that the resulting mesh will only have piecewise uniform :math:`\Delta \theta^\ast`, i.e., the mesh spacing :math:`\Delta \theta^\ast` on a given flux-surface can have discrete jumps. It is useful to relate this property to Fourier analysis. One can define a 2D Fourier decomposition in the two angle coordinates :math:`\varphi` and :math:`\theta^\ast`, .. math:: \lambda_{n,m}(\varphi,\theta^\ast) = \Lambda_{n,m} \exp[i(n\varphi - m \theta^\ast)], where :math:`n` and :math:`m` are the toroidal and poloidal mode numbers. Using the relation :math:`\Delta \theta^\ast = \Delta \varphi/q(\psi_N)` from above, the variation of :math:`\lambda_{n,m}` along a magnetic field line can be expressed as either :math:`\exp[i(n-m/q)\varphi]`, or :math:`\exp[i(nq-m)\theta^\ast]`. Both expressions show that Fourier modes with :math:`m=nq` are constant along the magnetic field, i.e, :math:`k_\parallel=0`. In addition, the first one shows that the parallel mode number, i.e., the number of parallel mode periods over one toroidal circuit along the magnetic field, is :math:`n_\parallel = |m/q-n| \approx R_0 k_\parallel`. The range of Fourier modes that are adequately resolved on the WHOLE flux-surface (the local resolution may be higher) is determined by: * The maximum of :math:`\Delta \theta^\ast` on the surface: :math:`m_{max}(\psi_N) = \frac{2\pi}{2 \max[\Delta \theta^\ast(\psi_N,\theta^\ast)]}`; hence, :math:`|n| \leq m_{max}(\psi_N)/q(\psi_N)`. * And, :math:`n_{\parallel,max} = \max(|m/q-n|) = \frac{\min[N_\varphi(\psi_N,\theta^\ast)]}{2}`. Conclusion: How to Choose the Toroidal Resolution """"""""""""""""""""""""""""""""""""""""""""""""" The observations in the previous sections can be condensed into the following recommendations for choosing ``num_phi_planes``: * If the goal is to resolve perturbations with :math:`|k_\theta \rho_i|<\alpha`, estimate the corresponding poloidal mode number :math:`m`, e.g., along the outer midplane. Using the safety factor profile, which is contained in the EFIT geqdsk file (or otherwhise known), :math:`n=m/q`, and :math:`N_\varphi \geq 2n`. * If the goal is to resolve modes up to a given toroidal mode number :math:`n`, use :math:`m \sim nq` and :math:`\Delta \theta^\ast < \pi/m`. Then estimate :math:`|\nabla(\theta^\ast)|` along the outer midplane to obtain an estimate for the required poloidal mesh spacing. * In any case, use :math:`k_\parallel \approx |m/q-n|/R_0` to adjust above estimates if a constraint on :math:`k_\parallel` is used. These recommendations can only yield approximations for the settings required for a particular simulation. After generating a mesh, it is recommended to calculate its :math:`\theta^\ast` and check whether the resolution is what is needed. If not, adjust the meshing parameters and iterate until the desired result is obtained. Please note that the number of radial-poloidal planes used in TOMMS (``num_phi_planes``) does not have to be identical to the number of planes in XGC (``sml_nphi_total``). However, for proper field-alignment, ``num_phi_planes`` must be divisible by the product of ``sml_nphi_total`` and ``sml_wedge_n``. Similarly, the effective upper bound on :math:`k_\parallel` in the XGC simulation (up to the intrinsic resolution limit of the mesh itself) is determined by the XGC parameters ``sml_nphi_total*sml_wedge_n``. Trouble-shooting ^^^^^^^^^^^^^^^^ The most common errors that would cause TOMMS to crash are: * Too large poloidal mesh size on the separatrix. When this error is triggered, it causes a crash very early in the execution of TOMMS because generating the separatrix curves is among the first tasks in the mesh generation. This error is easily resolved by decreasing the mesh size on the separatrix. The output at the time of the crash might look similar to this: .. code-block:: bash finite_domain_on = 0 correct to get a counter-clockwise wall curve Generating separatrix curves from 0 th x-point, startpoint size = 2 pushed 4 points from x-point :1.30398, -1.17765 pushed 4 points from x-point :1.2818, -1.18771 seg_front = 1 ||seg_back = 1 seg front index = 2 * DEBUG data for MATLAB, separatrix_100.txt is being written separatrix info: nxPts=2, xPts = 18, 906, curve ID = 100, typeId = 2, subtypeId = 1 convert XGC data starts # regions created btw core and wall = 2 List of physics regions = 7 3 0x2c760b0 has starting points = 54 0x2c765d0 has starting points = 92 Multiple open curves (2) with same psi value 1.00409 * A flux-surface that intersects the wall curve close to a model vertex, e.g. close to a corner of the wall curve, or almost tangentially. When this error is triggered, TOMMS prints an error message with the :math:`(R,Z)` coordinate of the offending position. To resolve this error, find out which of the flux-surfaces TOMMS generated is causing the error. Then, if using ``radial_uniform_meshing_unit=-1``, slightly change the :math:`\psi_N` value of that flux-surface and try running TOMMS again. Iterate until the error is resolved. * With ``intra_curve_spacing_option=-2`` (meshing for axisymmetric simulations), TOMMS can crash when the upper boundary for the dynamic poloidal spacing adjustment (``intra_curve_spacing_propf_max``) is too large for the given poloidal mesh size. Either decrease ``intra_curve_spacing_propf_max`` or the poloidal mesh size. Iterate until the error is resolved. 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. * Look for bad mesh quality. Example Workflow (IDL) ^^^^^^^^^^^^^^^^^^^^^^ .. code-block:: idl ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Set file names ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; fnefit='[PATH TO GEQDSK FILE]' fname='[SOME NAME]' ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Read EFIT geqdsk file (see XGC-Devel/utils/IDL/retrieve_efit.pro ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; retrieve_efit,fnefit,psirz,rgrid,zgrid,nrgrid,nzgrid,boundr,boundz,limitr,limitz, $ rmaxis,zmaxis,psi_xp,psiaxis,qpsi,ppsi,pprimepsi,fpsi,ffprimepsi, $ rzero,bzero ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Calculate psi_N, B_phi, B_pol, |B| along the outer midplane ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; nr=(size(rgrid))(1) psi_mid=dblarr(nr) for i=0,nr-1 do psi_mid(i)=spline(zgrid,redarr(psirz(i,*)),zmaxis) ;psi_mid=(psi_mid-psiaxis)/(psi_xp-psiaxis) psi0=spline(rgrid,psi_mid,rmaxis) ;; ;; Normalize such that psi=0 on the magnetic axis and ;; 1 on the separatrix ;; psi_mid=(psi_mid-psi0)/(psi_xp-psi0) bphi=fpsi/rgrid bpol=abs(deriv(rgrid,psi_mid*(psi_xp-psiaxis))/rgrid) btot=sqrt(bphi^2+bpol^2) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Set up fine psi-grid and interpolate magnetic field ;;; to that grid. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; nr2=1000 ;; rgrid2=dindgen(nr2)/(nr2-1)*(max(rgrid)-rmaxis)+rmaxis psi_mid=psi_mid>0 psi2=(spline(rgrid,psi_mid,rgrid2)) psi2(0)=0d0 bphi2=interpol(bphi,rgrid,rgrid2) bpol2=interpol(bpol,rgrid,rgrid2) btot2=interpol(btot,rgrid,rgrid2) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Repeat the calculation for the magnetic field and psi along ;;; the inner midplane ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; rgrid2a = dindgen(nr2)/(nr2-1)*(min(limitr)-rmaxis)+rmaxis ix = where(deriv(psi_mid) le 0) psi2a = reverse(spline(rgrid(ix),psi_mid(ix),reverse(rgrid2a))) bphi2a = interpol(bphi,rgrid,rgrid2a) bpol2a = interpol(bpol,rgrid,rgrid2a) btot2a = interpol(btot,rgrid,rgrid2a) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Read the experimental rofiles ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; read_file,'[DENSITY FILE]' ,2,den_in read_file,'[T_e FILE]' ,2,te_in read_file,'[T_i FILE]' ,2,ti_in read_file,'[OMEGA_TOR FILE]' ,2,om_in ;; Set normalization coefficients in case the ;; input profiles are not in m^-3, eV, or rad/s --> den_norm = 1d0 t_norm = 1d0 om_norm = 1d0 ;; Normalize and interpolate to common psi-grid psi_in = [den_in(*,0)] den_in = den_in(*,1)*den_norm ti_in = interpol(ti_in(*,1) ,ti_in (*,0),psi_in)*t_norm te_in = interpol(te_in(*,1) ,te_in (*,0),psi_in)*t_norm om_in = interpol(om_in(*,1) ,om_in (*,0),psi_in)*om_norm ;; npsi_in = n_elements(psi_in) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Extrapolate profiles into SOL ;;; (if necessary) by fitting an ;;; exponential decay to the input profiles ;;; such that value and gradient are mathed ;;; at the reference point (here: separatrix). ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; Find separatrix index ;; isep = round(interpol(dindgen(npsi_in),psi_in,1)) ;; ;; dden = (deriv(psi_in,den_in))(isep) denx = den_in(isep) dti = (deriv(psi_in,ti_in ))(isep) tix = ti_in(isep) dte = (deriv(psi_in,te_in ))(isep) tex = te_in(isep) dom = (deriv(psi_in,om_in ))(isep) omx = om_in(isep) ;; ;; Set upper bound for profile extrapolation ;; psi_max=1.2d0 ;; ;; Set up psi-grid in the SOL (1psi_max ;; i0 = max(where(psi2 lt psi_max)) den=spline(sqrt(psi_in2),den_in2,sqrt(psi2)) > 9d17 den(i0+1:*)=den(i0) te =spline(sqrt(psi_in2), te_in2,sqrt(psi2)) > 10 te(i0+1:*)=te(i0) ti =spline(sqrt(psi_in2), ti_in2,sqrt(psi2)) > 20 ti(i0+1:*)=ti(i0) om =spline(sqrt(psi_in2), om_in2,sqrt(psi2)) om(i0+1:*)=om(i0) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Some basic quantities that may be helpful ;;; for setting up the mesh ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; charge = 1.6022d-19 ;; elementary charge mi = 3.34d-27 ;; deuterium mass me = 9.109d-31 ;; electron mass pe = charge*den*(te) ;; electron pressure beta_e = 2*4d-7*!dpi*pe/btot2^2 ;; electron beta ;; ;; Safety factor ;; np = (size(psirz))(1) qs = spline(dindgen(np)/(np-1),qpsi,psi2) i0 = min(where(psi2 ge 1)) qs(i0:*) = qs(i0-1)*0.8 ;; heuristic extrapolation!!! ;; ;; Use inner midplane B for Alfven speed ;; (because it is much higher) ;; btot2b = interpol(btot2a,psi2a,psi2) va = btot2b/sqrt(4d-7*!dpi*den*mi) ;; Alfen speed vth = sqrt(2*ti*charge/mi) ;; thermal speed taua = rmaxis/va ;; Alfven time taut = 2*!dpi*rmaxis/vth ;; Ion toroidal transit time (thermal) taut_norm = 2*!dpi*rmaxis $ ;; Ion toroidal transit time (200 eV) *sqrt(mi/(2*200*charge))*(dblarr(nr2)+1d0) ;; ;; Calculate the ion gyroradius and the ;; density and temperature gradient scale lengths ;; rhoi = sqrt(mi*ti*charge)/(charge*btot2) ;; Outer midplane gyroradius rhoia = sqrt(mi*charge*interpol(ti,psi2,psi2a))/(charge*btot2a) rhoib = interpol(rhoia,psi2a,psi2) ;; Inner midplane gyroradius i0 = max(where(psi2 le max(psi2a))) rhoib(i0+1:*) = rhoib(i0) ;; simple extrapolation ln_1 = abs(deriv(rgrid2,den))/den ;; Inverse density gradient scale length lte_1 = abs(deriv(rgrid2,te))/te ;; Inv. T_e grad. scale length lti_1 = abs(deriv(rgrid2,ti))/ti ;; Inv. T_i grad. scale length pi = charge*den*ti pe = charge*den*te ptot = pi+pe lptot_1 = abs(deriv(rgrid2,ptot))/ptot ;; Inv. pressure grad. scale length ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Compare some useful time scales ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; loadct,3 sml_dt = 1d-4 plot,psi2,sml_dt*taut_norm,xra=[0.,1.1],xtitle='psi_n',ytitle='!3Timescales (ms)', $ th=3,font=0,posi=[0.25,0.15,0.95,0.85],xth=3,yth=3,/ylog,yra=[1d-9,1d-2] oplot,psi2,sml_dt*taut_norm*8,co=0,th=3,li=2 oplot,psi2,taut,co=200,th=3 oplot,psi2,taua,co=150,th=3 oplot,psi2,taua/5,co=150,th=3,li=2 legend,['!3Sim. time step','!38x sim. time step','!3Transit time','!3Alfven time','!30.2x Alfven time'], $ color=[0,0,200,150,150],linest=[0,2,0,0,2],psym=[0,0,0,0,0], th=3,chars=0.75,font=0 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Turbulence ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; Set up an array with the radial mesh size ;; relative to the ion gyroradius ;; ;; psi_N dr/rho_i ;; dat = [[0.00d0 , 1.00d0], $ [0.10d0 , 1.00d0], $ [0.20d0 , 1.00d0], $ [0.30d0 , 1.00d0], $ [0.40d0 , 1.00d0], $ [0.50d0 , 1.00d0], $ [0.60d0 , 1.00d0], $ [0.70d0 , 1.00d0], $ [0.80d0 , 1.00d0], $ [0.88d0 , 1.00d0], $ [0.90d0 , 0.95d0], $ [0.97d0 , 0.85d0], $ [1.00d0 , 0.80d0], $ [1.02d0 , 0.90d0], $ [1.07d0 , 2.00d0], $ [1.12d0 , 3.00d0], $ [1.20d0 , 4.00d0] $ ] psidum = transpose(dat(0,*)) dr1 = transpose(dat(1,*))*interpol(rhoi,psi2,psidum) dr1 = dr1>1d-3 rdum=interpol(rgrid2,sqrt(psi2),sqrt(psidum)) ;; ;; Radial resolution on fine psi-grid ;; dr = spline(sqrt(psidum),dr1,sqrt(psi2)) suffix='_turb' ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; XGCa Grid ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; psi_N dr/rho_i ;; dat = [[0.00d0 , 3.00d0], $ [0.10d0 , 3.00d0], $ [0.20d0 , 3.00d0], $ [0.30d0 , 3.00d0], $ [0.40d0 , 3.00d0], $ [0.50d0 , 3.00d0], $ [0.60d0 , 2.85d0], $ [0.70d0 , 2.60d0], $ [0.80d0 , 2.40d0], $ [0.88d0 , 2.15d0], $ [0.90d0 , 1.95d0], $ [0.97d0 , 1.00d0], $ [1.00d0 , 0.80d0], $ [1.02d0 , 1.00d0], $ [1.07d0 , 2.50d0], $ [1.12d0 , 4.00d0], $ [1.20d0 , 5.00d0] $ ] psidum = transpose(dat(0,*)) dr1 = transpose(dat(1,*))*interpol(rhoi,psi2,psidum) dr1 = dr1>1d-3 rdum=interpol(rgrid2,sqrt(psi2),sqrt(psidum)) ;; ;; Radial resolution on fine psi-grid ;; dr = spline(sqrt(psidum),dr1,sqrt(psi2)) suffix='_axisym' ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Visualize the radial resolution ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; plot,psi2,rhoi,xra=[0.,1.2d0],yra=[1d-3,1d-1],/ylog,xtitle='psi_n',ytitle='Length (m)',th=3, $ font=0,posi=[0.25,0.15,0.95,0.85] oplot,psi2,rhoi/2,li=2,th=3 oplot,psi2,rhoi*2,li=2,th=3 oplot,psi2,1d0/(2*ln_1),co=200,th=3 oplot,psi2,1d0/(2*lte_1),co=150,th=3 oplot,psi2,dr,li=1,th=2 oplot,[0.9,0.9],[1d-10,1d0],li=1,th=2 oplot,[1.0,1.0],[1d-10,1d0],li=1,th=2 legend,['rho_i','rho_i/2','L_n/2','L_Te/2','dr'],linestyle=[0,2,0,0,1],psym=[0,0,0,0,0], $ color=[0,0,200,150,0],th=3,chars=0.75,font=0,/bottom_leg,/left_leg ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Set up a set of flux-surfaces that ;;; has the targeted radial mesh size ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; Outer cutoff ;; psi_max=1.20D0 rmax=interpol(rgrid2,sqrt(psi2),sqrt(psi_max)) ;; ;; Step outward by dr(psi) starting from the separatris ;; (psi=1) to generate flux-surfaces in the SOL ;; psiout=1D0 rout=interpol(rgrid2,sqrt(psi2),1) drout=interpol(dr,sqrt(psi2),1) i=1 while (rout(i-1) lt rmax and psiout(i-1) lt psi_max) do begin $ rout=[rout,rout(i-1)+drout(i-1)] &$ psiout=[psiout,(interpol(sqrt(psi2),rgrid2,rout(i)))^2] &$ drout=[drout,interpol(dr,rgrid2,rout(i))] &$ i++ &$ endwhile nout=i ;; ;; Now step inward by dr(psi) starting from the separatrix ;; to generate the flux-surfaces inside the separatrix ;; dr_now = interpol(dr,sqrt(psi2),sqrt(psiout(0))) &$ while (rout(0) gt rmaxis+2*dr_now and psiout(0) gt 0D0) do begin $ dr_now = interpol(dr,sqrt(psi2),sqrt(psiout(0))) &$ rout=[rout(0)-dr_now,rout] &$ psiout=[(interpol(sqrt(psi2),rgrid2,rout(0)))^2,psiout] &$ drout=[interpol(dr,rgrid2,rout(0)),drout] &$ nout++ &$ endwhile ;; ;; Add the magnetic axis ;; psiout = [0d0,psiout] ;; psi-values of the flux-surfaces rout = [rmaxis,rout] ;; corresponding cylindrical R-coordinate drout = [rout(1)-rout(0),drout] ;; corresponding local radial mesh size nout++ ;; ;; Make sure that psi=1 is in the set ;; i0=where(abs(psiout-1) eq min(abs(psiout-1))) psiout(i0)=1D0 rout(i0)=interpol(rgrid2,sqrt(psi2),sqrt(psiout(i0))) drout(i0)=interpol(dr,sqrt(psi2),sqrt(psiout(i0))) ;; ;; Write input file for TOMMS ;; openw,u,fname+suffix+'_surf.dat',/get_lun printf,u,format='(i8)',nout for i=0,nout-1 do printf,u,format='(e22.12)',psiout(i) close,u & free_lun,u ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Set the poloidal target resolution ;;; Note that the following is subject ;;; to iteration until you get a mesh that ;;; fits your resolution needs and computing ;;; time constraints. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; pol_res_fac = dindgen(nout) ;; used for dl_theta=pol_res_fac*dr rhoiout = interpol(rhoi ,sqrt(psi2),sqrt(psiout)) ;; Outer midplane rho_i on psi_out rhoibout = interpol(rhoib,sqrt(psi2),sqrt(psiout)) ;; Inner midplane rho_i on psi_out isep = min(where(psiout ge 1)) qout = interpol(qpsi,dindgen(513)/512,psiout) ;; safety factor on psi_out grid qout(isep+1:*) = 0.5*(total(qout(isep-1:isep))) ;; heuristic extrapolation of q fpsiout = interpol(fpsi,dindgen(513)/512,psiout) ;; B_T*R on psi_out grid fpsiout(isep+1:*) = fpsiout(isep) ;;; ;;; Non-axisymmetric simulation ;;; ;;; Example 1: ;;; Aim at same k_max*rho_i in radial and poloidal direction at outer midplane and ;;; make sure that poloidal resolution supports the poloidal mode number associated ;;; with the outer midplane k_max*rho_i. ;;; dl_pol/dtheta^* scales like R^2*|B_pol| --> the required ratio between inboard and outboard dl_pol ;;; is (R^2*|B_pol|)_inner/(R^2 |B_pol|)_outer; in contrast, the gyroradius scales only ;;; like R_inner/R_outer, which is usually the weaker condition. ;;; ==> Choose dl_pol~dr_out --> same outer midplane k_max*rho_i radially and poloidally ;;; Choose poloidal spacing tolerance in mesh generator as f1/f2<1 ;;; ;;; This is dl_theta/dtheta* at the outer midplane: ;;; dl_dth = qout*(rout^2*interpol(bpol2,psi2,psiout)/abs(fpsiout)) ;;; ;;; Use this to estimate the poloidal mode number corresponding to ;;; a given mesh size ;;; m_target = 1d0*dl_dth/rhoiout ;;; ;;; And this is the corresponding toroidal mode number ;;; (use an upper limit) ;;; n_target = (m_target/qout)<384d0 ;;; ;;; Poloidal resolution to resolve m=m_target with ;;; 8 vertices per mode period: ;;; dl_out = dl_dth*(2*!dpi/(8*qout*n_target)) ;;; ;;; dl_pol/dtheta^* is proportional to R^2*|B_pol| ;;; inner midplane --> f1 ;;; outer midplane --> f2 ;;; ==> f1/f2 is the ratio of inner to outer poloidal mesh size ;;; of a field-aligned mesh ;;; This ratio can be useful in setting the poloidal target resolution ;;; and the intra-curve spacing tolerance ;;; f1 = interpol(rgrid2a^2*bpol2a,psi2a,psi2) f2 = rgrid2^2*bpol2 ;;; ;;; Can use some heuristic radial-poloidal anisotropy based on region of interest ;;; (larger anisotropy where we care less about accurate physics) ;;; pol_res_fac = 1d0+(1d0*0.5d0*(tanh(2*(psiout-1.08)/0.05d0)+1d0)+0.4*(1-0.5*(tanh(2*(psiout-0.3)/0.3)+1))) ;;; ;;; Optionally scale with f1/f2 ;;; pol_res_fac*= interpol(f1/f2,psi2,psiout) ;;; ;;; ;;; Example 2: ;;; ;;; Target: n=3, q_sep~7.5 ==> m_max~22 ;;; Aim at resolving poloidal mode numbers up to roughly m=25 ;;; ==> ~200 points per surface --> dl_pol~2*pi*r_min/200 npt = 200l pol_res_fac = (2*!dpi*(rout-rmaxis)/npt*interpol(f1/f2,psi2,psiout))/drout>0.5 pol_res_fac = (pol_res_fac<5) ;;; ;;; Example 3: axisymmetric simulation ;;; Aim at a given number of vertices per surface (function of psi_N) ;;; npt1 = 512 npt0 = 32 npt = 0.5*((npt1-npt0)*tanh(2*(psiout-0.5)/0.3)+(npt1+npt0)) pol_res_fac = (2*!dpi*(rout-rmaxis)/npt)/drout>1 pol_res_fac = (pol_res_fac<3) dl_out = drout*pol_res_fac ;;; ;;; Write input file for TOMMS ;;; openw,u,fname+suffix+'_dpol.dat',/get_lun printf,u,format='(i8)',nout ;;; Here, we add some logic to prevent outlandish values ;;; especially close to the magnetic axis and in the SOL for i=0,nout-1 do printf,u,format='(2(e22.12,1x))',psiout(i),max([1d-3,min([0.5*(rout(i)-rmaxis),pol_res_fac(i)*dl_out(i),1d-2])]) close,u & free_lun,u ;;; ;;; A rough estimate of the number of vertices with above mesh resolution ;;; np=lonarr(nout) npsum=np np(0)=1 npsum(0)=1 for i=1,nout-1 do begin $ np(i)=1.5*2*!dpi*(rout(i)-rout(0))/(drout(i)*pol_res_fac(i)) &$ npsum(i)=npsum(i-1)+np(i) &$ endfor ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Write XGC profile files ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; den_out = interpol(den,sqrt(psi2),sqrt(psiout)) ti_out = interpol(ti ,sqrt(psi2),sqrt(psiout)) te_out = interpol(te ,sqrt(psi2),sqrt(psiout)) om_out = interpol(om ,sqrt(psi2),sqrt(psiout)) openw,u,fname+'_den.dat',/get_lun printf,u,format='(i8)',nout for i=0,nout-1 do printf,u,format='(2(e22.12,1x))',psiout(i),den_out(i) printf,u,format='(i8)',-1 close,u & free_lun,u openw,u,fname+'_ti.dat',/get_lun printf,u,format='(i8)',nout for i=0,nout-1 do printf,u,format='(2(e22.12,1x))',psiout(i),ti_out(i) printf,u,format='(i8)',-1 close,u & free_lun,u openw,u,fname+'_te.dat',/get_lun printf,u,format='(i8)',nout for i=0,nout-1 do printf,u,format='(2(e22.12,1x))',psiout(i),te_out(i) printf,u,format='(i8)',-1 close,u & free_lun,u openw,u,fname+'_om.dat',/get_lun printf,u,format='(i8)',nout for i=0,nout-1 do printf,u,format='(2(e22.12,1x))',psiout(i),om_out(i) printf,u,format='(i8)',-1 close,u & free_lun,u ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; .. |xrArr| unicode:: U+27F9 .. Long Right Double Arrow .. |rarr| unicode:: U+2192 .. Short Right Arrow Generating Stellarator Meshes for XGC-S --------------------------------------- Meshes for stellarator simulations are generated differently than tokamak meshes because each plane has a different mesh in stellarator geometry. The stellarator mesh generator creates XGC-S mesh files from a VMEC NetCDF magnetic equilibrium file. The source files for the XGC-S meshing tool are located in `XGC-Devel/XGC-S/utils/VMEC2XGC2 `_. The tool relies on the `STELLOPT `_ package (used for stellarator optimization) which must be built first. This guide describes how to build and run VMEC2XGC2 on the Princeton University CPU cluster Stellar with the Slurm Workload Manager, but the procedure should be almost the same on other machines except for the modules to use. Building VMEC2XGC2 on Stellar ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ First note that VMEC2XGC2 is already installed on Stellar in ``/projects/XGC/STELLAR/Software/src/STELLOPT/VMEC2XGC2/Release/xvmec2xgc2`` so these instructions are merely to describe how to do it yourself. Also note that VMEC2XGC2 uses other modules than XGC-S on Stellar. **Download and install STELLOPT** * Pull the STELLOPT repository. .. code-block:: bash > cd /projects/XGC/STELLAR/Software/src/ > git clone git@github.com:PrincetonUniversity/STELLOPT.git * Load modules and set environment variables. STELLOPT needs ``STELLOPT_PATH`` set to the path of your directory where you have pulled the repository and ``MACHINE`` to the name of the machine to select the right makefile. The loaded modules, besides the compiler, are ``mpi``, ``fftw``, ``hdf5`` and ``netcdf``. .. code-block:: bash > source /projects/XGC/STELLAR/Software/bin/set_up_vmec2xgc2.stellar > echo $STELLOPT_PATH /projects/XGC/STELLAR/Software/src/STELLOPT > echo $MACHINE stellar * Build STELLOPT (this is likely to take some time because STELLOPT includes several codes). .. code-block:: bash > cd /projects/XGC/STELLAR/Software/src/STELLOPT > ./build_all -j16 **Install VMEC2XGC2 inside STELLOPT** * Copy VMEC2XGC2 to the STELLOPT directory and build. Note that one has to create symbolic links in the ``Debug`` and ``Release`` subfolders of VMEC2XGC2 to make the build work. .. code-block:: bash > cp -r /XGC-Devel/XGC-S/utils/VMEC2XGC2 /projects/XGC/STELLAR/Software/src/STELLOPT/ > cd /projects/XGC/STELLAR/Software/src/STELLOPT/VMEC2XGC2 > cd Debug > ln -s ../makevmec2xgc makevmec2xgc > cd ../Release > ln -s ../makevmec2xgc makevmec2xgc > cd .. > make clean_release * The VMEC2XGC2 executable is in ``/projects/XGC/STELLAR/Software/src/STELLOPT/VMEC2XGC2/Release/xvmec2xgc2``. Running VMEC2XGC2 on Stellar ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A VMEC NetCDF output file normally starts with ``wout_`` and ends with ``.nc``. Here we run the mesh generator on a W7-X equilibrium ``wout_w7x.nc``. * Go to a ``/scratch/`` location by choice and add the equilibrium file. .. code-block:: bash > cd /scratch/gpfs/xgc/STELLAR/xgc-s/vmec2xgc2_test > ls wout_w7x.nc * Start an interactive job. .. code-block:: bash > salloc --account=pppl --nodes=1 --qos=stellar-debug --time=00:30:00 * Load modules. .. code-block:: bash > source /projects/XGC/STELLAR/Software/bin/set_up_vmec2xgc2.stellar * To display all input options for xvmec2xgc2 run: .. code-block:: bash > srun --ntasks=1 /projects/XGC/STELLAR/Software/src/STELLOPT/VMEC2XGC2/Release/xvmec2xgc2 -help * Example to generate mesh files for XGC-S (in the current folder): .. code-block:: bash > srun --ntasks=1 /projects/XGC/STELLAR/Software/src/STELLOPT/VMEC2XGC2/Release/xvmec2xgc2 -vmec wout_w7x.nc -nrad 64 -spline 100 -nu 300 -nv 37 -out w7x_test -truncate 0.95 -thres1 0.2 -thres2 0.4 -extra 1 -align -1 -format 2 * Finish interactive job. .. code-block:: bash > exit * Now the following mesh files should be in the directory: * ``w7x_test_****.dat`` - magnetic field data files in cylindrical :math:`(R,Z)` coordinates in planes ``****``. * ``w7x_test_****.ele`` - mesh element files containing the triangle connectivity in planes ``****``. * ``w7x_test_****.node`` - mesh node files containing :math:`(R,Z)` coordinates of the mesh vertices in planes ``****``. * ``w7x_test.grid`` - flux-surface data file. * ``w7x_test.iota`` - another flux-surface data file which is read by XGC-S but not used at the moment. There can be other files generated, e.g. ``extra.debug*``, ``*.test``, ``poloidal_mesh.data``, ``spline_diff.dat``, ``w7x_test_****a.poly``, but these are not needed for XGC-S. **Explanation of common input flags to xvmec2xgc2:** .. list-table:: :widths: 35 25 100 :header-rows: 1 * - *Input flag* - Default - Description * - *-vmec* - None - Filename of the VMEC NetCDF magnetic equilibrium input file. * - *-nrad* - None - Radial resolution of the unstructured mesh, i.e. number of flux-surfaces (if -radr is not specified). The flux-surfaces are placed equidistantly in normalized minor radius :math:`\sqrt{\psi/\psi_a}`. Alternatively to -nrad, -input can be used to manually place flux-surfaces at specific radial locations. * - *-input* - None - This flag overrides -nrad and can be used when a non-equidistant radial resolution is wanted. It takes the name of a file containing the flux-surface locations. The file should contain a header line with the number of flux-surfaces and 1 dummy integer. The following lines should contain 4 columns: index of the flux-surface (starting from 1), normalized minor radius :math:`\sqrt{\psi/\psi_a}`, number of mesh vertices on the flux-surface per cross-section, rotational transform (iota) of the flux-surface. If -align -1 is used the iota value is ignored so dummy data can be specified in the last column. An example file is given below. * - *-spline* - 100 - Resolution on each toroidal cross-section :math:`(R,Z)` in cylindrical coordinates (if -rads not specified). * - *-nu* - 36 - Resolution of the unstructured mesh in poloidal direction (if -radu is not specified). This is the number of mesh vertices along the outermost flux-surface on each cross-section. * - *-nv* - 36 - Resolution of the unstructured mesh in toroidal direction per field period and for the field data in the cylindrical coordinate. The first cross-section is identical to the last due to periodicity, e.g. using 36 implies that there are 35 toroidal cross-sections. When running XGC-S ``sml_nphi_total`` must be set to match the value of this flag, see `XGC-S Input Parameters `__ for information. * - *-out* - xgc_mesh - 8 character string to use as prefix for the output files. * - *-truncate* - 1.0 - Normalized minor radius :math:`\sqrt{\psi/\psi_a}`, i.e. square root of normalized poloidal flux (label), of outermost flux-surface with mesh vertices (should be :math:`\leq 1` ). When running XGC-S ``sml_outpsi`` must be set to match the value of this flag, see `XGC-S Input Parameters `__ for information. * - *-thres1* - 0.3 - Specifies the connection scheme among mesh vertices together with -thres2. For normalized minor radius :math:`\sqrt{\psi/\psi_a}` 0.0 to -thres1 connection is based on poloidal angle (field-aligned). For -thres1 to -thres2 connection is field-aligned and partially optimized for triangle shape. For normalized minor radius -thres2 to 1.0 connection is optimized for triangle shape. * - *-thres2* - 0.9 - Specifies the connection scheme among mesh vertices together with -thres1. * - *-extra* - 2 - Extrapolation scheme (0 - 3). Normally use option 1. * - *-align* - 1 - Poloidal discretization scheme. -1: discretization in VMEC angle. 0: discretization in PEST angle. 1: discretization in PEST angle + safety factor (field-aligned). Option -1 is recommended, option 0 and 1 can sometimes experience issues with convergence (seen by messages "error in Newton method" in the output). * - *-format* - 0 - Output data file format. 0: binary without header. 1: binary with header. 2: ASCII with header. Recommended to use option 2. .. code-block:: bash 116 0 1 0.0000000000000000 1 0.85406845076578197 2 1.4456615074101151E-002 6 0.85408570690255092 3 2.8913196426081207E-002 18 0.85413776686616283 4 4.3369723606038038E-002 22 0.85422556061328181 5 5.7826164705973813E-002 30 0.85435081200137253 [...] i (psi/psi_a)^(1/2) N_{node} iota [...] 114 0.96835826327620989 642 0.93436583907214810 115 0.98102340081213069 614 0.93654868687994164 116 1.0000000000000000 570 0.93990568385888495 *Example of an -input file for xvmec2xgc2 used to manually specify flux-surfaces.* Calculation of spline interpolation coefficients for the field data ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ XGC-S also needs files containing spline coefficients for the field data. The make and source files for the tool to generate the spline coefficients files are located in `XGC-Devel/XGC-S/utils/gen_spline_coeff `_. The tool needs an installation of ``PSPLINE`` (on Stellar located in /projects/XGC/STELLAR/Software/install) and the makefile coef.cmp should be updated with the path to ``PSPLINE`` and the Fortran compiler on the system. * Build and run spline coefficients tool in the folder. The tool will ask for input/output filenames. (The tool only needs the ``w7x_test_****.dat`` files as input.) .. code-block:: bash > cd /scratch/gpfs/xgc/STELLAR/xgc-s/vmec2xgc2_test > cp XGC-Devel/XGC-S/utils/gen_spline_coeff/* /scratch/gpfs/xgc/STELLAR/xgc-s/vmec2xgc2_test/ > ./coef.cmp > ./make_coef.go input field file name (8 characters) w7x_test output coefficient file name (8 characters) w7x_test * Now the following spline files (in binary format) should be generated in the directory: * ``w7x_test_****.mag1`` - magnetic field data files containing :math:`B_R` coefficients in planes ``****``. * ``w7x_test_****.mag2`` - magnetic field data files containing :math:`B_\varphi` coefficients in planes ``****``. * ``w7x_test_****.mag3`` - magnetic field data files containing :math:`B_Z` coefficients in planes ``****``. * ``w7x_test_****.flux`` - normalized poloidal flux (label) :math:`\psi` coefficients in planes ``****``. Collect the output files from VMEC2XGC2 and the spline tool and put them in a subfolder of the ``equilibrium`` folder for the XGC-S simulation.