3.1 Mathematical basics

To solve partial differential equations (PDEs), the py-pde package provides differential operators to express spatial derivatives. These operators are implemented using the finite difference method to support various boundary conditions. The time evolution of the PDE is then calculated using the method of lines by explicitly discretizing space using the grid classes. This reduces the PDEs to a set of ordinary differential equations, which can be solved using standard methods as described below.

3.1.1 Curvilinear coordinates

The package supports multiple curvilinear coordinate systems through pde.grids.coordinates. They allow to exploit symmetries present in physical systems. Consequently, many grids implemented in py-pde inherently assume symmetry of the described fields. However, a drawback of curvilinear coordinates are the fact that the basis vectors now depend on position, which makes tensor fields less intuitive and complicates the expression of differential operators. To avoid confusion, we here specify the used coordinate systems explicitly:

3.1.1.1 Polar coordinates

Polar coordinates describe points by a radius \(r\) and an angle \(\phi\) in a two-dimensional coordinates system. They are defined by the transformation

\[\begin{split}\begin{cases} x = r \cos(\phi) &\\ y = r \sin(\phi) & \end{cases} \text{for} \; r \in [0, \infty] \; \text{and} \; \phi \in [0, 2\pi)\end{split}\]

The associated symmetric grid PolarSymGrid assumes that fields only depend on the radial coordinate \(r\). Note that vector and tensor fields can still have components in the polar direction. In particular, vector fields still have two components: \(\vec v(r) = v_r(r) \vec e_r + v_\phi(r) \vec e_\phi\).

3.1.1.2 Spherical coordinates

Spherical coordinates describe points by a radius \(r\), an azimuthal angle \(\theta\), and a polar angle \(\phi\). The conversion to ordinary Cartesian coordinates reads

\[\begin{split}\begin{cases} x = r \sin(\theta) \cos(\phi) &\\ y = r \sin(\theta) \sin(\phi) &\\ z = r \cos(\theta) \end{cases} \text{for} \; r \in [0, \infty], \; \theta \in [0, \pi], \; \text{and} \; \phi \in [0, 2\pi)\end{split}\]

The associated symmetric grid SphericalSymGrid assumes that fields only depend on the radial coordinate \(r\). Note that vector and tensor fields can still have components in the two angular direction.

Warning

Not all results of differential operators on vectorial and tensorial fields can be expressed in terms of fields that only depend on the radial coordinate \(r\). In particular, the gradient of a vector field can only be calculated if the azimuthal component of the vector field vanishes. Similarly, the divergence of a tensor field can only be taken in special situations.

3.1.1.3 Cylindrical coordinates

Cylindrical coordinates describe points by a radius \(r\), an axial coordinate \(z\), and a polar angle \(\phi\). The conversion to ordinary Cartesian coordinates reads

\[\begin{split}\begin{cases} x = r \cos(\phi) &\\ y = r \sin(\phi) &\\ z = z \end{cases} \text{for} \; r \in [0, \infty], z \in \mathbb{R}, \; \text{and} \; \phi \in [0, 2\pi)\end{split}\]

The associated symmetric grid CylindricalSymGrid assumes that fields only depend on the coordinates \(r\) and \(z\). Vector and tensor fields still specify all components in the three-dimensional space.

Warning

The order of components in the vector and tensor fields defined on cylindrical grids is different than in ordinary math. While it is common to use \((r, \phi, z)\), we here use the order \((r, z, \phi)\). It might thus be best to access components by name instead of index.

3.1.2 Spatial discretization

The finite differences scheme used by py-pde is currently restricted to orthogonal coordinate systems with uniform discretization. Because of the orthogonality, each axis of the grid can be discretized independently. For simplicity, we only consider uniform grids, where the support points are spaced equidistantly along a given axis, i.e., the discretization \(\Delta x\) is constant. If a given axis covers values in a range \([x_\mathrm{min}, x_\mathrm{max}]\), a discretization with \(N\) support points can then be though of as covering the axis with \(N\) equal-sized boxes; see inset. Field values are then specified for each box, i.e., the support points lie at the centers of the box:

\[\begin{split} x_i &= x_\mathrm{min} + \left(i + \frac12\right) \Delta x \quad \text{for} \quad i = 0, \ldots, N - 1 \\ \Delta x &= \frac{x_\mathrm{max} - x_\mathrm{min}}{N}\end{split}\]

which is also indicated in the inset.

Differential operators are implemented using the usual second-order central differences. This requires introducing virtual support points at \(x_{-1}\) and \(x_N\), which can be determined from the boundary conditions at \(x=x_\mathrm{min}\) and \(x=x_\mathrm{max}\), respectively. The field classes automate this transparently. However, if you need more control over boundary conditions, you can access the full underlying data using the field._data_full, which will have \(N + 2\) entries along an axis that has \(N\) support points. In this case, the first and last entries (data_full[0] and data_full[N + 1]) denote the lower and upper virtual point, respectively. The actual field data can be obtained using data_full[1:-1] or the field.data attribute for convenience. Note that functions evaluating differential operators generally expect the full data as input while they return only valid data.

3.1.3 Temporal evolution

Once the fields have been discretized, the PDE reduces to a set of coupled ordinary differential equations (ODEs), which can be solved using standard methods. This reduction is also known as the method of lines. The py-pde package implements the simple Euler scheme and a more advanced Runge-Kutta scheme in the ExplicitSolver class. For the simple implementations of these explicit methods, the user typically specifies a fixed time step, although adaptive methods, which adjust the time step automatically, are also often used and available in the package. One problem with explicit solvers is that they require small time steps to stably evolve some PDEs; such PDEs are then often called ‘stiff’. Stiff PDEs can sometimes be solved more efficiently by using implicit methods. This package provides a simple implementation of the Backward Euler method in the ImplicitSolver class. Finally, more advanced methods are available by wrapping the scipy.integrate.solve_ivp() in the ScipySolver class.