 |
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 Poisson solver on 2D polar mesh; uses FFT in theta and 2nd-order FD in r.
This module solves the Poisson equation \( -\nabla^2 \phi = \rho \) (permittivity of free space \( \epsilon_0 = 1 \)) on a 2D polar mesh \( (r,\theta) \in [r_\textrm{min},r_\textrm{max}]\times[0,2\pi) \). The Poisson equation in polar coordinates reads
\[ -\bigg( \frac{\partial^2}{\partial r^{2}} +\frac{1}{r}\frac{\partial}{\partial r} +\frac{1}{r^{2}}\frac{\partial^{2}}{\partial\theta^{2}} \bigg)\phi(r,\theta) = \rho(r,\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 \).
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{\partial^2}{\partial r^{2}} +\frac{1}{r}\frac{\partial}{\partial r} -\frac{k^{2}}{r^{2}} \bigg)\widehat{\phi}_k(r) = \widehat{\rho}_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 Poisson equation is differentiated to obtain an expression for \( \partial^3\widehat{\phi}_k / \partial r^3 \) that involves only 1st and 2nd derivatives,
\[ \frac{\partial^3\widehat\phi_k}{\partial r^3} = -\frac{1}{r} \frac{\partial^2\widehat\phi_k}{\partial r^2} +\frac{1+k^2}{r^2} \frac{\partial \widehat\phi_k}{\partial r } -\frac{2 k^2}{r^3} { \widehat\phi_k} - \frac{\partial \widehat\rho_k}{\partial r }, \]
which are then approximated with the same finite-difference formulas.
Derived types and interfaces | |
| type | sll_t_poisson_2d_polar |
| Class for the Poisson solver in polar coordinate. More... | |
Functions/Subroutines | |
| subroutine, public | sll_s_poisson_2d_polar_init (solver, rmin, rmax, nr, ntheta, bc_rmin, bc_rmax, rgrid) |
| sll_o_initialize the Poisson solver in polar coordinates More... | |
| subroutine, public | sll_s_poisson_2d_polar_solve (solver, rho, phi) |
| Solve the Poisson equation and get the electrostatic potential. More... | |
| subroutine, public | sll_s_poisson_2d_polar_free (solver) |
| Delete contents (local storage) of Poisson's solver. More... | |
| type(sll_t_poisson_2d_polar) function, pointer, public | sll_f_new_poisson_2d_polar (rmin, rmax, nr, ntheta, bc_r, rgrid) |
| OO interface: allocate pointer to Poisson solver and initialize it. More... | |
| subroutine | s_compute_phi_from_rho (poisson, phi, rho) |
| OO interface: solve Poisson's equation. More... | |
| subroutine | s_compute_e_from_rho (poisson, E1, E2, rho) |
| OO interface: solve Poisson's equation and compute E field. More... | |
| real(f64) function | f_l2norm_squared (poisson, coefs_dofs) |
| OO interface: compute L2 norm squared of something (...) More... | |
| subroutine | s_compute_rhs_from_function (poisson, func, coefs_dofs) |
| OO interface: project 2D function onto Finite Element space. More... | |
| subroutine | s_free (poisson) |
|
private |
OO interface: compute L2 norm squared of something (...)
Definition at line 650 of file sll_m_poisson_2d_polar.F90.
|
private |
OO interface: solve Poisson's equation and compute E field.
Definition at line 640 of file sll_m_poisson_2d_polar.F90.
|
private |
OO interface: solve Poisson's equation.
Definition at line 618 of file sll_m_poisson_2d_polar.F90.
|
private |
OO interface: project 2D function onto Finite Element space.
Definition at line 660 of file sll_m_poisson_2d_polar.F90.
|
private |
Definition at line 669 of file sll_m_poisson_2d_polar.F90.
| type(sll_t_poisson_2d_polar) function, pointer, public sll_m_poisson_2d_polar::sll_f_new_poisson_2d_polar | ( | real(f64), intent(in) | rmin, |
| real(f64), intent(in) | rmax, | ||
| integer(i32), intent(in) | nr, | ||
| integer(i32), intent(in) | ntheta, | ||
| integer(i32), dimension(2), intent(in) | bc_r, | ||
| real(f64), dimension(:), intent(in), optional | rgrid | ||
| ) |
OO interface: allocate pointer to Poisson solver and initialize it.
| [in] | nr | number of cells radial |
| [in] | ntheta | number of cells angular |
| [in] | bc_r | boundary conditions at [r_min,r_max] |
| [in] | rgrid | grid points along r |
Definition at line 586 of file sll_m_poisson_2d_polar.F90.
| subroutine, public sll_m_poisson_2d_polar::sll_s_poisson_2d_polar_free | ( | type(sll_t_poisson_2d_polar), intent(inout) | solver | ) |
Delete contents (local storage) of Poisson's solver.
Definition at line 568 of file sll_m_poisson_2d_polar.F90.
| subroutine, public sll_m_poisson_2d_polar::sll_s_poisson_2d_polar_init | ( | type(sll_t_poisson_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), optional | rgrid | ||
| ) |
sll_o_initialize the Poisson solver in polar coordinates
| [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] | rgrid | grid points along r |
Definition at line 269 of file sll_m_poisson_2d_polar.F90.
| subroutine, public sll_m_poisson_2d_polar::sll_s_poisson_2d_polar_solve | ( | type(sll_t_poisson_2d_polar), intent(inout) | solver, |
| real(f64), dimension(:, :), intent(in) | rho, | ||
| real(f64), dimension(:, :), intent(out) | phi | ||
| ) |
Solve the Poisson equation and get the electrostatic potential.
| [in,out] | solver | Solver object |
| [in] | rho | Charge density |
| [out] | phi | Potential |
Definition at line 480 of file sll_m_poisson_2d_polar.F90.
1.9.1