sll_m_poisson_2d_polar_par Module Reference

Description

Parallel Poisson solver on 2D polar mesh; uses FFT in theta and 2nd-order FD in r.

Authors

Yaman Güçlü, IPP Garching

Edoardo Zoni, IPP Garching

Model equations

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 \).

Numerical methods

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.

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.

Parallelization

• The user must provide two compatible 2D layouts:
  • layout_a sequential in \( \theta \);
  • layout_r sequential in \( r \);
• Input and output data are given in layout_a;
• layout_r is only used internally;
• FFTs are computed on one row of layout_a;
• The linear solver (for each mode \( k \)) works on one column of layout_r;
• The maximum number of parallel solve is \( N_\theta \), where \( N_\theta\) is the number of cells along \( \theta \). Therefore, instead of working with \( N_\theta/2+1 \) Fourier modes, we work with \( N_\theta \) real modes. For this purpose, the r2r interface of FFTW is used: this corresponds to DFTs of real input and complex-Hermitian output in halfcomplex format.

Usage example

use sll_m_poisson_2d_polar_par, only: &
  sll_t_poisson_2d_polar_par,         &
  sll_s_poisson_2d_polar_par_init,    &
  sll_s_poisson_2d_polar_par_solve,   &
  sll_s_poisson_2d_polar_par_free
 
use sll_m_remapper, only: &
  sll_t_layout_2d
 
...
 
type(sll_t_poisson_2d_polar_par) :: solver
type(sll_t_layout_2d)            :: layout_a, layout_r
real(f64), allocatable           :: rho(:,:), phi(:,:)
real(f64)                        :: rmin, rmax
integer(i32)                     :: nr, ntheta, bc_rmin, bc_rmax
 
...
 
call sll_s_poisson_2d_polar_par_init ( solver, layout_r, layout_a, rmin, rmax, nr, ntheta, bc_rmin, bc_rmax )
call sll_s_poisson_2d_polar_par_solve( solver, rho, phi )
call sll_s_poisson_2d_polar_par_free ( solver )

Derived types and interfaces
type	sll_t_poisson_2d_polar_par
	Class for the Poisson solver in polar coordinate. More...

Functions/Subroutines
subroutine, public	sll_s_poisson_2d_polar_par_init (solver, layout_r, layout_a, 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_par_solve (solver, rho, phi)
	Solve the Poisson equation and get the electrostatic potential. More...
subroutine, public	sll_s_poisson_2d_polar_par_free (solver)
	Delete contents (local storage) of Poisson's solver. More...
subroutine	verify_argument_sizes_par (layout, array, array_name)
	Check if array sizes are compatible with the layout. More...

Function/Subroutine Documentation

sll_s_poisson_2d_polar_par_free()

subroutine, public sll_m_poisson_2d_polar_par::sll_s_poisson_2d_polar_par_free

(

type(sll_t_poisson_2d_polar_par), intent(inout)

solver

)

Delete contents (local storage) of Poisson's solver.

Definition at line 644 of file sll_m_poisson_2d_polar_par.F90.

sll_s_poisson_2d_polar_par_init()

subroutine, public sll_m_poisson_2d_polar_par::sll_s_poisson_2d_polar_par_init	(	type(sll_t_poisson_2d_polar_par), intent(out)	solver,
		type(sll_t_layout_2d), pointer	layout_r,
		type(sll_t_layout_2d), pointer	layout_a,
		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, target	rgrid
	)

sll_o_initialize the Poisson solver in polar coordinates

Parameters

[out]	solver	solver object
	layout_r	layout sequential in r direction
	layout_a	layout sequential in theta direction
[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 285 of file sll_m_poisson_2d_polar_par.F90.

Here is the call graph for this function:

sll_s_poisson_2d_polar_par_solve()

subroutine, public sll_m_poisson_2d_polar_par::sll_s_poisson_2d_polar_par_solve	(	type(sll_t_poisson_2d_polar_par), intent(inout)	solver,
		real(f64), dimension(:, :), intent(in)	rho,
		real(f64), dimension(:, :), intent(out), target	phi
	)

Solve the Poisson equation and get the electrostatic potential.

Parameters

[in,out]	solver	Solver object
[in]	rho	Charge density
[out]	phi	Potential

Definition at line 544 of file sll_m_poisson_2d_polar_par.F90.

Here is the call graph for this function:

verify_argument_sizes_par()

subroutine sll_m_poisson_2d_polar_par::verify_argument_sizes_par ( type(sll_t_layout_2d), pointer  layout, real(f64), dimension(:, :), intent(in)  array, character(len=*), intent(in)  array_name  )

private

Check if array sizes are compatible with the layout.

Definition at line 669 of file sll_m_poisson_2d_polar_par.F90.

Here is the caller graph for this function: