 |
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
|
|
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
|
Serial quasi-neutrality solver on 2D polar mesh; uses FFT in theta and 2nd-order FD in r.
This module solves the quasi-neutrality equation on a 2D polar mesh \( (r,\theta) \in [r_\textrm{min},r_\textrm{max}]\times[0,2\pi) \). The general quasi-neutrality equation is of the form
\[ -\nabla_\perp\cdot \bigg[\frac{\rho_{m0}}{\epsilon_0 B^2}\nabla_\perp\phi\bigg] +\frac{1}{\lambda_D^2}\big[\phi-\chi\langle\phi\rangle_f\big] = \frac{1}{\epsilon_0}\rho_{c1}, \]
where \( \phi \) is the electrostatic potential, \( \langle\phi\rangle_f \) is the flux-surface average of \( \phi \), \( \rho_{c1} \) is the charge perturbation density due to the kinetic species, \( \rho_{m0} \) is the equilibrium mass density, \( B \) is the intensity of the equilibrium background magnetic field, \( \epsilon_0 \) is the permittivity of free space and \( \lambda_D \) is the electron Debye length
\[ \lambda_D\equiv\sqrt{\frac{\epsilon_0\kappa T_e}{q_e^{2}n_{e0}}}. \]
The equation in polar coordinates reads
\[ -\bigg[ \frac{g}{r}\frac{\partial}{\partial r} +\frac{\partial}{\partial r}\bigg(g\frac{\partial}{\partial r}\bigg) +\frac{g}{r^2}\frac{\partial^2}{\partial\theta^2} \bigg]\phi(r,\theta) +\frac{1}{\lambda_D^2}\big[\phi(r,\theta)-\chi\langle\phi\rangle_\theta(r)\big] = \frac{1}{\epsilon_0}\rho_{c1}(r,\theta), \]
where \( g\equiv\rho_{m0}/\epsilon_0 B^{2} \) is assumed to be independent of \( \theta \) and the flux-surface average is replaced by a simple average over \( \theta \).
The boundary conditions (BCs) on \( \phi(r,\theta) \) are set as follows. \( \phi(r,\theta) \) is \( 2\pi\)-periodic along \( \theta \) and the BCs along \( r \) can be chosen among the following types ( \( \overline{r} = r_\textrm{min} \) or \( \overline{r} = r_\textrm{max} \)):
Option name | Description | Condition on solution | Condition on Fourier modes |
|---|---|---|---|
sll_p_dirichlet | Homogeneous Dirichlet | \( \phi(\overline{r},\theta)\equiv 0\) | \( \widehat{\phi}_k(\overline{r})= 0\quad \forall k \) |
sll_p_neumann | Homogeneous Neumann | \( \partial_r\phi(\overline{r},\theta)\equiv 0\) | \( \partial_r\widehat{\phi}_k(\overline{r})=0\quad \forall k \) |
sll_p_neumann_mode_0 | Homogeneous Neumann mode 0 | \( \partial_\theta\phi(\overline{r},\theta)\equiv 0 \) and \( \int_0^{2\pi} \partial_r\phi(\overline{r},\theta) d\theta = 0 \) | \( \partial_r\widehat{\phi}_0(\overline{r})=0 \) and \( \widehat{\phi}_k(\overline{r})=0 \) for \( k\neq 0 \) |
sll_p_polar_origin | \( \overline{r}=r_\textrm{min}=0 \) is center of circular disk | \( \phi(-\Delta r/2,\theta) \equiv \phi(\Delta r/2,\theta+\pi) \) | \( \widehat{\phi}_k(-\Delta r/2) = (-1)^k \widehat{\phi}_k(\Delta r/2)\quad \forall k \) |
Please note that the option sll_p_polar_origin can only be used at \( \overline{r} = r_\textrm{min} \) when one sets \( r_\textrm{min} = 0 \).
The following arguments are given by the user at initialization: \( \rho_{m0} \), \( B \), \( \lambda_D \), \( \epsilon_0 \) and the additional parameter \( \chi \) (default is \( \chi=1 \)). If \( \lambda_D \) is not given to the solver, we assume \( 1/\lambda_D^2=0 \): the electrons form a kinetic species and their contribution goes into \( \rho_{c1} \). If \( \epsilon_0 \) is not given to the solver, we assume \( \epsilon_0=8.854187817\times10^{-12} [F/m]\) according to the parameter sll_p_epsilon_0 in module sll_m_constants.
This module uses fast Fourier transforms (FFTs) in \( \theta \) and a 2nd-order finite-difference method in \( r \). Thanks to the linearity of the differential operator and the periodicity of the domain, a discrete Fourier transform (DFT) in \( \theta \) is applied to both sides of the above elliptic PDE. Then, each Fourier coefficient \( \widehat{\phi}_k(r) \) solves an independent 1D boundary value problem on \( [r_\textrm{min},r_\textrm{max}] \):
\[ -\bigg[ \frac{g}{r}\frac{\partial}{\partial r} +\frac{\partial}{\partial r}\bigg(g\frac{\partial}{\partial r}\bigg) -\frac{k^2}{r^2}g \bigg]\widehat{\phi}_k(r) +\frac{1}{\lambda_D^2}(1-\chi\,\delta_{k0})\widehat{\phi}_k(r) = \frac{1}{\epsilon_0}\widehat{\rho}_{c1,k}(r), \]
For each mode \( k \), the resulting ODE is solved with a 2nd-order finite-difference collocation method.
Since \( \phi(r,\theta) \) is real, we have \( \widehat{\phi}_{-k}=\overline{\widehat{\phi}_k} \) for each \( k=-N_\theta/2,-N_\theta/2+1,\dots,-1,0,1,\dots,N_\theta/2-1,N_\theta/2\), where \( N_\theta\) is the number of cells along \( \theta \). Therefore, it is enough to consider only \( N_\theta/2+1 \) coefficients \( \widehat{\phi}_0, \dots, \widehat{\phi}_{N_\theta/2} \). For this purpose, the r2c and c2r interfaces of FFTW are used for the forward and backward FFTs, respectively.
If the user does not provide a radial grid, a collocation grid with uniform spacing \( \Delta r \) is calculated from the input parameters \( r_\textrm{min} \), \( r_\textrm{max} \) and \(N_r\) (number of radial cells). In the case of a circular annulus ( \( r_\textrm{min} > 0\)), the collocation grid has \( N_r+1 \) points, with first point \( r_1 = r_\textrm{min} \) and last point \( r_{N_r+1} = r_\textrm{max} \). In the case of a circular disc ( \( r_\textrm{min} = 0 \)), the collocation grid has only \( N_r \) points, with first point \( r_1 = \Delta r/2 \) and last point \( r_{N_r} = r_\textrm{max} \).
For a given mode \(k\), at each interior point in the collocation grid ( \( r_i \neq r_\textrm{min}, r_\textrm{max} \)) the equation above is approximated using centered finite-differences on a 3-point stencil, valid on a general non-uniform grid with \( h_i^{+}=r_{i+1}-r_i \) and \( h_i^{-}=r_{i}-r_{i-1} \):
\[ \frac{\partial}{\partial r} \widehat\phi_k(r_i) = \frac{1}{h_i^+ + h_i^-} \left[ -\frac{h_i^+}{h_i^-}\, \widehat{\phi}_k(r_{i-1}) +\left( \frac{h_i^+}{h_i^-} - \frac{h_i^-}{h_i^+} \right)\, \widehat{\phi}_k(r_i) +\frac{h_i^-}{h_i^+}\, \widehat{\phi}_k(r_{i+1}) \right] \,+\, \left\{ -\frac{\partial^3}{\partial r^3} \widehat\phi_k(r_i) \frac{h_i^+ h_i^-}{6} \right\} \,+\, \textit{H.O.T.}, \]
\[ \frac{\partial^2}{\partial r^2} \widehat\phi_k(r_i) = \frac{1}{h_i^+ h_i^-} \left[ \frac{2h_i^+}{h_i^++h_i^-}\, \widehat{\phi}_k(r_{i-1}) -2 \widehat{\phi}_k(r_i) +\frac{2h_i^-}{h_i^++h_i^-}\, \widehat{\phi}_k(r_{i+1}) \right] \,+\, \left\{ -\frac{\partial^3}{\partial r^3} \widehat\phi_k(r_i) \frac{h_i^+ - h_i^-}{3} \, -\frac{\partial^4}{\partial r^4} \widehat\phi_k(r_i) \frac{(h_i^+)^2 - h_i^+ h_i^- +(h_i^-)^2}{12} \right\} \,+\, \textit{H.O.T.}. \]
On the right-hand side of each equation, the first term is the finite-difference approximation used in the code, the terms in curly braces represent the leading truncation error, and \( \textit{H.O.T} \) means "higher order terms". The coefficients that multiply \(\widehat{\phi}_k(r_{i-1})\), \(\widehat{\phi}_k(r_i)\) and \(\widehat{\phi}_k(r_{i+1})\) are the so-called "finite-difference coefficients". Although both formulas are exact for a parabolic profile, for a general profile the finite-difference approximation to the second derivative is only 1st-order accurate when \( h_i^+ \neq h_i^- \). In the future the finite-difference coefficients for the second derivative can be modified to cancel the 1st-order error term: the Fourier-transformed quasi-neutrality equation can be differentiated to obtain an expression for \( \partial^3\widehat{\phi}_k / \partial r^3 \) that involves only 1st and 2nd derivatives, which are then approximated with the same finite-difference formulas.
Derived types and interfaces | |
| type | sll_t_qn_solver_2d_polar |
| Class for 2D Poisson solver in polar coordinates. More... | |
Functions/Subroutines | |
| subroutine, public | sll_s_qn_solver_2d_polar_init (solver, rmin, rmax, nr, ntheta, bc_rmin, bc_rmax, rho_m0, b_magn, lambda, use_zonal_flow, epsilon_0, rgrid) |
| Initialize the solver. More... | |
| subroutine, public | sll_s_qn_solver_2d_polar_solve (solver, rho, phi) |
| Solve the quasi-neutrality equation and get the electrostatic potential. More... | |
| subroutine, public | sll_s_qn_solver_2d_polar_free (solver) |
| Delete contents (local storage) of quasi-neutrality solver. More... | |
| subroutine, public sll_m_qn_solver_2d_polar::sll_s_qn_solver_2d_polar_free | ( | type(sll_t_qn_solver_2d_polar), intent(inout) | solver | ) |
Delete contents (local storage) of quasi-neutrality solver.
Definition at line 640 of file sll_m_qn_solver_2d_polar.F90.
| subroutine, public sll_m_qn_solver_2d_polar::sll_s_qn_solver_2d_polar_init | ( | type(sll_t_qn_solver_2d_polar), intent(out) | solver, |
| real(f64), intent(in) | rmin, | ||
| real(f64), intent(in) | rmax, | ||
| integer(i32), intent(in) | nr, | ||
| integer(i32), intent(in) | ntheta, | ||
| integer(i32), intent(in) | bc_rmin, | ||
| integer(i32), intent(in) | bc_rmax, | ||
| real(f64), dimension(:), intent(in) | rho_m0, | ||
| real(f64), dimension(:), intent(in) | b_magn, | ||
| real(f64), dimension(:), intent(in), optional | lambda, | ||
| logical, intent(in), optional | use_zonal_flow, | ||
| real(f64), intent(in), optional | epsilon_0, | ||
| real(f64), dimension(:), intent(in), optional, target | rgrid | ||
| ) |
Initialize the solver.
| [out] | solver | solver object |
| [in] | nr | number of cells radial |
| [in] | ntheta | number of cells angular |
| [in] | bc_rmin | boundary condition at r_min |
| [in] | bc_rmax | boundary condition at r_max |
| [in] | rho_m0 | radial profile: total mass density of equilibrium |
| [in] | b_magn | radial profile: intensity of magnetic field |
| [in] | lambda | radial profile: electron Debye length |
| [in] | use_zonal_flow | if .false. set flux average to zero |
| [in] | epsilon_0 | override default: vacuum permittivity |
| [in] | rgrid | grid points along r |
Definition at line 285 of file sll_m_qn_solver_2d_polar.F90.
| subroutine, public sll_m_qn_solver_2d_polar::sll_s_qn_solver_2d_polar_solve | ( | type(sll_t_qn_solver_2d_polar), intent(inout) | solver, |
| real(f64), dimension(:, :), intent(in) | rho, | ||
| real(f64), dimension(:, :), intent(out) | phi | ||
| ) |
Solve the quasi-neutrality equation and get the electrostatic potential.
| [in,out] | solver | Solver object |
| [in] | rho | Charge density |
| [out] | phi | Potential |
Definition at line 549 of file sll_m_qn_solver_2d_polar.F90.
1.9.1