Package biossphere overview#

biosspheres is a python-based solver for Laplace and Helmholtz scattering by multiple disjoint spheres, utilizing spherical harmonic decomposition and local multiple trace formulations. It also provides routines for coupling the problems with ordinary differential equations in time.

It is divided into seven main modules, which are further subdivided:

  • formulations

  • helmholtz

  • laplace

  • miscella

  • quadratures

  • timecouplings

  • utils

In the following, a general background is given, with the main notations used along the code. After it, a small overview of each of the modules is made. The order of the sections are alphabetical with respect to the name of each module.

There are additional Jupyter notebooks that go into more detail of some functionalities. It is addressed when it is the case in this overview.

Geometry#

In most settings, we consider \(N\) disjoint spheres in \(\mathbb{R}^3\), with the following notation:

  • Position vectors of the sphere centers: \(\mathbf{p_j} \in \mathbb{R}^3\), with \(j \in \{1, ..., N\}\).

  • Radii: \(\tilde{r}_j \in \mathbb{R}^+\), with \(j \in \{1, ..., N\}\).

  • Interior of sphere \(j\): \(\Omega_j\), with \(j \in \{1, ..., N\}\).

  • Boundary of sphere \(j\): \(\Gamma_j := \partial \Omega_j\).

  • Exterior medium: \(\Omega_0\), defined as \(\mathbb{R}^3\) without the spheres and their boundaries.

When using \(N=1\), we sometimes use the sub-indexes \(e\) for exterior and \(i\) for interior.

Problems that can be solved#

The code was written with a focus of solving Helmholtz or Laplace problems with transmission boundary conditions, by using a system of boundary integral equations. However, the routines can be used to solve problems with other boundary conditions. Further details are given in the overviews of the modules laplace, helmoholtzand formulations.

An example volume problem that can be posed as a boundary integral system is:

We want to find \(u_j \in H^1_{loc}(\Omega_j)\), \(j \{0, ..., N \}\), such that they solve:

\[\text{div} (\sigma_j \nabla u_j) = 0 \quad \text{ in } \Omega_j, \text{ with } j \in \{0, ..., N \},\]
\[-\gamma_d^{0j} u_0 + \gamma_d^{j} u_j = \gamma_d^{0j} \phi_e \quad \text{ on } \Gamma_j, \text{ with } j \in \{1, ..., N \},\]
\[\sigma_0 \gamma_n^{0j} u_0 + \sigma_j \gamma_n^{j} u_j = - \sigma_0\gamma_n^{0j} \phi_e \quad \text{ on } \Gamma_j, \text{ with } j \in \{1, ..., N \},\]
\[\lim_{||\mathbf{x}||_2 \rightarrow \infty} u_0(\mathbf{x}) = \mathcal{O}\left(\frac{1}{||\mathbf{x}||_2}\right),\]

with

\[\gamma_d^{0j} u_0 := u_0|_{\Gamma_j},\]
\[\gamma_d^{j} u_j := u_j|_{\Gamma_j},\]
\[\gamma_n^{0j} u_0 := \nabla u_0|_{\Gamma_j} \cdot \widehat{n}_{0j},\]
\[\gamma_n^{j} u_j := \nabla u_j|_{\Gamma_j}\cdot \widehat{n}_{j},\]

where \(\widehat{n}_{j}\) is the exterior normal of \(\Omega_j\), with \(j\in \{ 1, ..., N\}\) and \(\widehat{n}_{0j}=-\widehat{n}_{j}\).

\(\sigma_j\), with \(j\in \{ 1, ..., N\}\), are positive parameters and \(\phi_e\) a given harmonic function, with domain in \(\Omega_0\).

Coordinate systems used#

Spherical coordinate system#

A vector is written as \(\mathbf{r}=\left(r,\varphi,\theta\right)^t\), with \(r \in [0,\infty)\), \(\varphi \in [0,2\pi)\) and \(\theta \in [0,\pi]\), which in Cartesian coordinates is equivalent to \(\mathbf{r}=r\left(\sin \theta \cos \varphi,\sin \theta \sin \varphi,\cos \theta\right)^t\). The unitary vectors of the spherical coordinate system can be written in Cartesian coordinates as

\[\begin{split}\begin{aligned} &\widehat{\mathbf{e}}_r= (\sin \theta \cos \varphi , \sin \theta \sin \varphi , \cos \theta )^t, \\ &\widehat{\mathbf{e}}_\theta=(\cos \theta \cos \varphi , \cos \theta \sin \varphi , -\sin \theta )^t, \\ &\widehat{\mathbf{e}}_{\varphi}=(-\sin \varphi , \cos \varphi , 0 )^t. \end{aligned}\end{split}\]

Also, the gradient operator in spherical coordinates:

\[\begin{aligned} \nabla f&=\frac{\partial f}{\partial r} \widehat{\mathbf{e}}_r+\frac{1}{r} \frac{\partial f}{\partial \theta} \widehat{\mathbf{e}}_\theta+\frac{1}{r \sin \theta} \frac{\partial f}{\partial \varphi} \widehat{\mathbf{e}}_{\varphi} \end{aligned}\]

Discretization#

Preliminaries#

First, we introduce the functions that form the trial and test functions that are assumed when computing the discretized versions of the boundary integral operators.

Associated Legendre functions#

\(P_l^m\) are the associated Legendre functions of degree \(l\) and order \(m\) defined as:

\[P_{l}^m\left(x\right) := (-1)^m \left( 1- x^2\right)^{\frac{m}{2}} \frac{d^m}{dx^m}P_l(x), \quad \text{with} \quad P_{l}\left(x\right) := \frac{1}{2^ll!}\frac{d^l}{dx^l}(x^2-1)^l.\]

Here, the term \((-1)^m\) is the Condon-Shortley phase factor.

Real spherical harmonics#

Real spherical harmonics of degree \(l\) and order \(m\) are defined using spherical coordinates:

\[\begin{split}\begin{aligned} Y_{l,m}\left(\theta,\varphi\right) &:= \sqrt{ (2-\delta_{m,0}) \frac{\left(2l+1\right)\left(l-m\right)!}{4 \pi \left(l+m\right)!}} P_l^{m} \left(\cos\theta\right) \cos m \varphi , \text{ and}\\ Y_{l,-m}\left(\theta,\varphi\right) &:= \sqrt{ (2-\delta_{m,0})\frac{\left(2l+1\right)\left(l-m\right)!}{4 \pi \left(l+m\right)!}} P_l^{m} \left(\cos\theta\right) \sin m \varphi , \end{aligned}\end{split}\]

with \(l\in \mathbb{N}_0\), \(m\in \mathbb{Z}\) such that \(0\leq m\leq l\). If \(m=0\), then \(\delta_{m,0}=1\), and it is zero otherwise.

Complex spherical harmonics#

Complex spherical harmonics of degree \(l\) and order \(m\) are defined using spherical coordinates:

\[Y_{l,m}\left(\theta,\varphi\right) := \sqrt{ \frac{\left(2l+1\right)\left(l-m\right)!}{4 \pi \left(l+m\right)!}} P_l^{m} \left(\cos\left(\theta\right)\right) e^{i m \varphi},\]
\[Y_{l,-m}\left(\theta,\varphi\right) := (-1)^m\overline{Y}_{l,m}\left(\theta,\varphi\right),\]

with \(l\in \mathbb{N}_0\), \(m\in \mathbb{Z}\) such that \(0\leq m\leq l\). If \(m=0\), then \(\delta_{m,0}=1\), and it is zero otherwise.

Orthonormality of spherical harmonics#

Spherical harmonics are dense in \(C(\mathbb{S}^2)\), with \(\mathbb{S}^2\) the surface of the unit sphere, and form a complete orthonormal system in \(L^2(\mathbb{S}^2)\) with respect to the internal product defined by:

\[\left( \psi , \xi \right)_{L^2(\mathbb{S}^2)} = \int_{0}^{2\pi}\int_{0}^{\pi} \psi\left(\theta,\varphi\right) \overline{\xi\left(\theta,\varphi\right) }\sin\left(\theta\right) d\theta d\varphi,\]

with \(\overline{\xi\left(\theta,\varphi\right)}\) the conjugate of \(\xi\left(\theta,\varphi\right)\).

They also are orthogonal in \(H^1(\mathbb{S}^2)\).

Spherical harmonic discretization#

Let be \(j\in \{1,...,N\}\). We define the reference system \(j\) as the one centered at \(\mathbf{p_j}\) with the same orientation that the reference system centered in the origin. Furthermore, we denote by \(Y_{l,m,j}\) the spherical harmonic \(Y_{l,m}\) centered in the origin of the reference system \(j\). Thus, if \(\left( r_j, \varphi_j, \theta_j \right)\) are the vector spherical coordinates of \(\mathbf{r_j}\) in the reference system \(j\), we have that \(Y_{l,m,j}\left(\mathbf{r}_j\right)=Y_{l,m}\left(\theta_j, \varphi_j\right)\).

For \(L \in \mathbb{N}_0\) and \(j\in \{1,...,N\}\), we define subspaces

\[\mathcal{Y}_L\left(\Gamma_j \right):= \text{span}\left\lbrace Y_{l,m,j}: l \in \mathbb{N}_0, m \in \mathbb{Z}, l \leq L, |m|\leq l \right\rbrace,\]

equipped with the \(L^2(\Gamma_j)\)-norm. Notice that the dimension of each subspace is \((L+1)^2\), and that the sequence of subspaces \(\lbrace \mathcal{Y}_L \left(\Gamma_j \right) \rbrace_{L \in \mathbb{N}_0}\) is dense in \(H^{\frac{1}{2}}(\Gamma_j)\) and in \(H^{-\frac{1}{2}}(\Gamma_j)\). This result justifies the discretization of all boundary Dirichlet and Neumann unknowns with spherical harmonics.

In the example given in the section “Problems that can be solved”, we can discretize all the traces using spherical harmonics:

For \(j \in \lbrace 1, ..., N \rbrace\), we write \(u^L_{d,0j}\), \(u^L_{n,0j}\), \(u^L_{d,j}\) and \(u^L_{n,j}\) in \(\mathcal{Y}_L(\Gamma_j)\) for the approximations of \(\gamma_d^{0j} u_0\), \(\gamma_n^{0j} u_0\), \(\gamma_d^{j} u_j\) and \(\gamma_n^{j} u_j\), respectively. They can be written as the following series expansions:

\[\begin{split}\begin{aligned} u^L_{d,0j}&=\sum_{l=0}^{L} \sum_{m=-l}^l u^{l,m}_{d,0j} Y_{l,m,j}, & u^L_{n,0j} &=\sum_{l=0}^{L} \sum_{m=-l}^l u^{l,m}_{n,0j} Y_{l,m,j},\\ u^L_{d,j}&=\sum_{l=0}^{L} \sum_{m=-l}^l u^{l,m}_{d,j} Y_{l,m,j},& u^L_{n,j} &=\sum_{l=0}^{L} \sum_{m=-l}^l u^{l,m}_{n,j} Y_{l,m,j}, \end{aligned}\end{split}\]

with \(u^{l,m}_{d,0j}\), \(u^{l,m}_{n,0j}\), \(u^{l,m}_{d,j}\), and \(u^{l,m}_{n,j}\) being constants.

Computational comment#

We mostly use the package pyshtools when computing the Associated Legendre Functions and when obtaining the coefficients of the expansion in spherical harmonics of a function when a mapping in a sphere is given.

1. formulations module#

The formulations module has the following submodules:

  • massmatrices: python script with routines to obtain discretization of mass matrices for one or several spheres when using spherical harmonics basis as test and trial functions in a boundary integral formulation setting. The routines only return the diagonal vector of the matrices, because the rest of the entries are zero.

  • mtf: module that is further subdivided. It has routines for solving Laplace and Helmholtz problems using the multiple traces formulation.

This structure was designed like this to allow the addition of other formulations by incorporating new modules. Since massmatrices can be used across several formulations it was placed in that position.

1.1 massmatrices#

Python script with routines to obtain discretization of mass matrices for one or several spheres when using spherical harmonics basis as test and trial functions in a boundary integral formulation setting. The routines only return the diagonal vector of the matrices, because the rest of the entries are zero.

1.2 mtf module#

The mtf folder has the following submodules:

  • mtf: python script with routines for obtaining linear operators arising in multiple traces formulations when using spherical harmonics basis as test and trial functions.

  • propertieschecks: python script with routines for checking properties that hold in a mtf setting. This helps the user to check the correctness of the obtained solutions.

  • reconstructions: python script with routines that use the representation formula to retrieve the values of the solution in volume points.

  • righthands: python script with routines for building righthand side of a mtf problem for some given external functions.

  • solvertemplates: python script with routines that take the parameters of a problem and return the solution. They build the discretization matrices, righthand sides and solve.

  • timecouplings: submodule with routines for coupling ODEs with a multitraces formulation.

In the jupyter notebook e1_laplace_transmission_mtf is presented how to solve a Laplace transmission problem with the multiple traces formulation, along with the small explanation of the formulation used.

1.2.6 Time couplings#

Examples:

  • e2_mtf_time_coupled_Kavian.ipynb

  • e3_mtf_time_coupled_e2_FitzHughNagumo.ipynb

2. helmholtz module#

The helmholtz module has two submodules:

  • selfinteractions: python script with routines to obtain the evaluation and testing of boundary integral operators with complex spherical harmonics in one sphere. The boundary integral operators have Helmholtz kernel. There are routines for the four operators needed for the Calderón projector. For details see the Jupyter notebook called 2_1_selfinteractions_helmholtz_overview.ipynb.

  • crossinteractions: python script with routines to obtain the evaluation and testing of boundary integral operators with complex spherical harmonics between two different spheres. The boundary integral operators have Helmholtz kernel. There are routines for the four operators needed for the Calderón projector. For details see the Jupyter notebook called 2_2_crossinteractions_helmholtz_overview.ipynb.

3. laplace module#

The laplace module has three submodules:

  • selfinteractions: python script with routines to obtain the evaluation and testing of boundary integral operators with real spherical harmonics in one sphere. The boundary integral operators have Laplace kernel. There are routines for the four operators needed for the Calderón projector. For details see the Jupyter notebook called 3_1_selfinteractions_laplace_overview.ipynb.

  • crossinteractions: python script with routines to obtain the evaluation and testing of boundary integral operators with real spherical harmonics between two different spheres. The boundary integral operators have Laplace kernel. There are routines for the four operators needed for the Calderón projector. For details see the Jupyter notebook called 3_2_crossinteractions_laplace_overview.ipynb.

  • drawings

References#

Boundary integral operators in spheres:

MTF:

Pyshtools:

Numpy:

Scipy:

Matplotlib: