Fast SVD-Based Linear Elastic Eigenvalue Problem Solver for Band Structures of 3D Phononic Crystals

Abstract

In this article, a Fast Linear Elastic Eigenvalue Problem Solver (FLEEPS) is developed to calculate the band structures of three-dimensional (3D) isotropic phononic crystals (PnCs). In brief, FLEEPS solves in linear time complexity the smallest few eigenvalues and associated eigenvectors of the linear elastic eigenvalue problem originating from the finite difference discretization of the frequency-domain linear elastic wave equation. Notably, FLEEPS employs the weighted singular value decomposition based preconditioner to greatly improve the convergence rate of the conjugate gradient iteration, and uses the fast Fourier transform algorithm to accelerate this preconditioner times a vector, based on the structured decomposition of the dense unitary factor T of this preconditioner. Band structure calculations of several 3D isotropic PnCs are presented to showcase the capabilities of FLEEPS. The preliminary MATLAB implementation of FLEEPS is available at https://github.com/FAME-GPU/FLEEPS-MATLAB.

Appendices

The Component-Wise Formulation of (1)

In order to facilitate the reader to derive (4), we write down the component-wise velocity-stress formulation of (1) in detail as follows:

$$\begin{aligned} \imath \omega \rho v_x&= \partial _x \tau _{xx} + \partial _y \tau _{xy} + \partial _z \tau _{xz},\quad \imath \omega \rho v_y= \partial _x \tau _{xy} + \partial _y \tau _{yy} + \partial _z \tau _{yz},\\ \imath \omega \rho v_z&= \partial _x \tau _{xz} + \partial _y \tau _{yz} + \partial _z \tau _{zz},\\ \imath \omega \tau _{xx}&=(\lambda +2\mu )\partial _x v_x + \lambda \partial _y v_y +\lambda \partial _z v_z, \quad \imath \omega \tau _{yz}= \mu (\partial _z v_y + \partial _y v_z),\\ \imath \omega \tau _{yy}&= \lambda \partial _x v_x + (\lambda +2\mu )\partial _y v_y+\lambda \partial _z v_z,\quad \imath \omega \tau _{xz}= \mu (\partial _z v_x + \partial _x v_z),\\ \imath \omega \tau _{zz}&= \lambda \partial _x v_x +\lambda \partial _y v_y +(\lambda +2\mu ) \partial _z v_z,\quad \imath \omega \tau _{xy}= \mu (\partial _y v_x + \partial _x v_y). \end{aligned}$$

Therefore, in some sense, only six components of $\tau (\textbf{r})$, i.e.,  $\tau _{xx}(\textbf{r}),\tau _{yy}(\textbf{r}),\tau _{zz}(\textbf{r})$, $\tau _{yz}(\textbf{r})$, $\tau _{xz}(\textbf{r})$ and $\tau _{xy}(\textbf{r})$, are actually needed to solve (1).

Definitions of $J_2$ and $J_3$ for 14 Bravais Lattices

For all 14 Bravais lattices, we can compute the QR factorization with column pivoting of $[\varvec{a}_1,\varvec{a}_2,\varvec{a}_3]$ such that

$$\begin{aligned} {[}\varvec{a}_1,\varvec{a}_2,\varvec{a}_3]\Pi =QR_3, \end{aligned}$$

where $\Pi$ is a suitable permutation, Q is orthonormal and $R_3$ is upper triangular. Define $S:=\text{ diag }\big \{\text{ sign }(R_3(1,1)),\text{ sign }(R_3(2,2)),\text{ sign }(R_3(3,3))\big \}$ and

$$\begin{aligned} {[}\widetilde{\varvec{t}}_1,\widetilde{\varvec{t}}_2,\widetilde{\varvec{t}}_3]=SR_3, \end{aligned}$$

whereby three nonnegative scalars $\rho _i$ ($i = 1, 2, 3$) are defined as follows:

$$\begin{aligned} \rho _1&= {\left\{ \begin{array}{ll} 0, &{} \quad \text{ if } \widetilde{\varvec{t}}_2(1) \ge 0, \\ 1, &{} \quad \text{ if } \widetilde{\varvec{t}}_2(1)< 0; \end{array}\right. } \ \rho _2 = {\left\{ \begin{array}{ll} 0, &{} \quad \text{ if } \widetilde{\varvec{t}}_3(1) \ge 0, \\ 1, &{} \quad \text{ if } \widetilde{\varvec{t}}_3(1)< 0; \end{array}\right. } \ \rho _3 = {\left\{ \begin{array}{ll} 0, &{} \quad \text{ if } \widetilde{\varvec{t}}_3(2) \ge 0, \\ 1, &{} \quad \text{ if } \widetilde{\varvec{t}}_3(2) < 0. \end{array}\right. } \end{aligned}$$

Furthermore, we define

$$\begin{aligned}{} & {} m_1:= \texttt{round}\left( \frac{\widetilde{\varvec{t}}_2(1) + \rho _1 \widetilde{\varvec{t}}_1(1)}{h_x} \right) \le n_1,\ m_2:= \texttt{round} \left( \frac{\widetilde{\varvec{t}}_3(1)+\rho _2 \widetilde{\varvec{t}}_1(1)}{h_x} \right)< n_1,\\{} & {} m_3:= \texttt{round} \left( \frac{\widetilde{\varvec{t}}_3(2)+\rho _3 \widetilde{\varvec{t}}_2(2)}{h_y}\right) < n_2, \end{aligned}$$

where round denotes rounding to the nearest integer. Given a wave vector $\textbf{k}$, using these notations, matrices $J_2$ and $J_3$, with those in (10b) and (11b) for the FCC lattice as the special case, for all 14 Bravais lattices are defined as follows:

$$\begin{aligned} J_2=e^{\imath 2 \pi \textbf{k}\cdot (\rho _1 \widetilde{\varvec{t}}_1 + \widetilde{\varvec{t}}_2)} { \begin{bmatrix} &{} e^{-\imath 2\pi \textbf{k} \cdot \widetilde{\varvec{t}}_1} I_{m_1} \\ I_{n_1-m_1} &{} \end{bmatrix}}; \end{aligned}$$

• (i) for $\widetilde{\varvec{t}}_3(2) \ge 0$ and $m_{4} \equiv m_2 - m_1 \ge 0$,

  $$\begin{aligned} \hspace{-10mm} J_3 = { \begin{bmatrix} &{} e^{-\imath 2\pi \textbf{k}\cdot ({\widetilde{\varvec{t}}}_{2} + (\rho _{1} - \rho _{2}){\widetilde{\varvec{t}}}_{1})}I_{m_{3}}\otimes J_{2,m_4} \\ e^{\imath 2 \pi \rho _2 \textbf{k}\cdot \widetilde{\varvec{t}}_1} I_{n_{2} - m_{3}}\otimes J_{2,m_2} &{} \end{bmatrix}}; \end{aligned}$$

• (ii) for $\widetilde{\varvec{t}}_3(2) \ge 0$ and $m_{4} < 0$,

  $$\begin{aligned} J_3&={ \begin{bmatrix} &{} e^{-\imath 2\pi \textbf{k}\cdot ({\widetilde{\varvec{t}}}_{2} + (\rho _{1}-\rho _{2} - 1){\widetilde{\varvec{t}}}_{1})}I_{m_{3}}\otimes J_{2,n_1+m_4} \\ e^{\imath 2 \pi \rho _2 \textbf{k}\cdot \widetilde{\varvec{t}}_1} I_{n_{2} - m_{3}}\otimes J_{2,m_2} &{} \end{bmatrix}}; \end{aligned}$$

• (iii) for $\widetilde{\varvec{t}}_3(2) < 0$ and $m_{5} \equiv m_1 + m_2 \le n_{1}$,

  $$\begin{aligned} \hspace{-10mm} J_3 = { \begin{bmatrix} &{} e^{\imath 2 \pi \rho _2 \textbf{k}\cdot \widetilde{\varvec{t}}_1} I_{m_{3}}\otimes J_{2,m_2}\\ e^{\imath 2\pi \textbf{k}\cdot (\widetilde{\varvec{t}}_2+(\rho _2+\rho _1)\widetilde{\varvec{t}}_1)}I_{n_{2} - m_{3}}\otimes J_{2,m_5} &{} \end{bmatrix}}; \end{aligned}$$

• (iv) for $\widetilde{\varvec{t}}_3(2) < 0$ and $m_{5} > n_{1}$,

  $$\begin{aligned} J_3&={ \begin{bmatrix} &{} e^{\imath 2 \pi \rho _2 \textbf{k}\cdot \widetilde{\varvec{t}}_1} I_{m_{3}}\otimes J_{2,m_2}\\ e^{\imath 2\pi \textbf{k}\cdot (\widetilde{\varvec{t}}_2+(\rho _2+\rho _1-1)\widetilde{\varvec{t}}_1)}I_{n_{2} - m_{3}}\otimes J_{2,m_5-n_1} &{} \end{bmatrix}}, \end{aligned}$$

where $J_{2,m} \equiv { \begin{bmatrix} &{} e^{-\imath 2\pi \textbf{k}\cdot \widetilde{\varvec{t}}_1} I_{m} \\ I_{n_1-m} &{} \end{bmatrix}}$.

Some Details on Derivations of Lemmas 3, 4 and Theorem 5

To derive eigendecompositions in Lemmas 3, 4 and Theorem 5 in Sect. 4.1, here we provide some necessary proofs and matrix properties.

Lemma 18

${\begin{bmatrix} &{}\hspace{-0.5cm} e^{-\imath 2\pi \textbf{k}\cdot \varvec{a}_1} I_{q}\\ I_{n_1-q} &{} \end{bmatrix}}\!\!=\! (K_1^*)^q,q=1,\ldots ,n_1$, with $K_1$ defined in (9). Hence, $J_2$ defined in (10b) satisfies

$$\begin{aligned} J_2= J_2^{-*} = K_1^{-n_1/2}. \end{aligned}$$

(C1)

Proof

By applying mathematical induction with respect to q. $\square$

Corollary 19

$(I_{n_2}\!\otimes \! K_1)K_2\!=\!K_2(I_{n_2}\!\otimes \!K_1),C_1C_2=C_2C_1$ and $C_1C_2^*=C_2^*C_1$.

Proof

By direct verification and using (C1), we have $(I_{n_2}\!\otimes \! K_1)K_2\!=\!K_2(I_{n_2}\!\otimes \!K_1)$ and $(I_{n_2}\!\otimes \! K_1)K_2^*\!=\!K_2^*(I_{n_2}\!\otimes \!K_1)$, therefore, $C_1C_2=C_2C_1$ and $C_1C_2^*=C_2^*C_1$, using the definitions of $C_1$ and $C_2$. $\square$

Lemma 20

It holds that ${\begin{bmatrix} &{}\hspace{-1cm} e^{-\imath 2\pi \textbf{k}\cdot \varvec{a}_2} I_{q}\otimes J_2^*\\ I_{n_2-q}\otimes I_{n_1} &{} \end{bmatrix}}=(K_2^*)^q$, $q=1,\ldots ,n_2$, with $K_2$ and $J_2$ in (10b).

Proof

By applying mathematical induction with respect to q. $\square$

Corollary 21

With $K_2$ and $J_3$ in (10b) and (11b), respectively, we have

$$\begin{aligned} J_3 = K_2^{-n_2/3}(I_{n_2}\otimes K_1)^{-n_1/2}=(I_{n_2}\otimes K_1)^{-n_1/2}K_2^{-n_2/3}=J_3^{-*}. \end{aligned}$$

(C2)

Proof

$J_3 (I_{n_2}\!\otimes \!K_1^{n_1/2})\!=\! J_3(I_{n_2}\!\otimes \! J_2^*)\!=\!(I_{n_2}\!\otimes \! K_1^{n_1/2})\! J_3\! =\! (I_{n_2}\!\otimes \! J_2^*)\! J_3$, which further equals ${\begin{bmatrix} &{} \hspace{-1.2cm}e^{-\imath 2\pi \textbf{k}\cdot \varvec{a}_2}\!I_{n_2/3}\otimes J_2^*\\ I_{2n_2/3}\otimes I_{n_1} &{} \end{bmatrix}} = (K_2^*)^{n_2/3}$, by Lemma 20. $\square$

Corollary 22

$C_3$ commutes with $C_1, C_1^*, C_2$ and $C_2^*$.

Proof

By (C2), $J_3$ commutes with $I_{n_2}\otimes K_1, I_{n_2}\otimes K_1^*, K_2$ and $K_2^*$, hence, by direct verification, we have $C_3 C_\ell =C_\ell C_3$ and $C_3C_\ell ^*=C_\ell ^*C_3$, $\ell =1,2$. $\square$

Theorem 23

$\{C_\ell , C_\ell ^*\}_{\ell =1}^3$ is a set of commutative matrices.

Proof

By Corollary 19, Corollary 22 and the normality of $C_1,C_2$ and $C_3$. $\square$

Therefore, $C_\ell$ and $C_\ell ^*$, $\ell =1,2,3$, can be simultaneously diagonalized by a common unitary matrix. which is determined instantly.

Lemma 24

([30]) Given $p,q\in \mathbb {N}$, let $M(x)=\sum _{s=0}^{q-1}x^s M_s + x^q I_p$ with $M_s\in \mathbb {C}^{p\times p},s=0,1,\ldots ,q-1$, then $\det M(x) = \det (xI_{mq} - C_{BF})$, with

$$\begin{aligned} C_{BF}:= \begin{bmatrix} 0 &{} \quad I_p &{} \quad \cdots &{} \quad 0 \\ \vdots &{} \quad \vdots &{} \quad \ddots &{} \quad \vdots \\ 0 &{} \quad 0 &{} \quad \cdots &{} \quad I_p \\ -M_0 &{} \quad -M_1 &{} \quad \cdots &{} \quad -M_{q-1} \\ \end{bmatrix}. \end{aligned}$$

Moreover, if $\textbf{x}\in \mathbb {C}^p$ and $\beta _0\in \mathbb {C}$ satisfy $M(\beta _0)\textbf{x}=0$, then the eigenvector of $C_{BF}$ associated with the eigenvalue $\beta _0$ is $[\beta _0^0,\ldots ,\beta _0^{q-1}]^\top \otimes \textbf{x}$.

Proof of Lemma 3

By direct verification, and noting that $K_1$ is the companion matrix of the monic polynomial $(x^{n_1} - e^{\imath 2\pi \textbf{k}\cdot \varvec{a}_1})$. Moreover, for $i,i'\in \mathbb {N}_1$ but $i\ne i'$, we have $e^{\imath \theta _i}\ne e^{\imath \theta _{i'}}$, then $\textbf{x}_i^*\textbf{x}_{i'}=0$, since $K_1$ is normal. $\square$

Proof of Lemma 4

Note that $K_2$ defined in (10b) can be seen as the block companion matrix of the monic matrix polynomial $(x^{n_2}I_{n_1}-e^{\imath 2\pi \textbf{k}\cdot \varvec{a}_2}J_2)$. Hence, by (C1), Lemma 24 and Lemma 3, $\det (x^{n_2}I_{n_1}-e^{\imath 2\pi \textbf{k}\cdot \varvec{a}_2}J_2)=0$ if and only if $x^{n_2}= e^{\imath n_1\theta _i/2}$ for $i\in \mathbb {N}_1$, i.e.,  $x=e^{\imath \theta _{ji}}$ with $\theta _{ji}$ defined in (19a), and the associated eigenvector of $K_2$ is just $\textbf{y}_{ji}$ for $i\in \mathbb {N}_1, j\in \mathbb {N}_2$. Due to the orthogonality of $\textbf{x}_i$, we have $(\textbf{y}_{j'i'}\otimes \textbf{x}_{i'})^*(\textbf{y}_{ji}\otimes \textbf{x}_i) =0$ if $i'\ne i$. Moreover, for $j,j'\in \mathbb {N}_2$ but $j\ne j'$, we have $e^{\imath \theta _{ji}}\ne e^{\imath \theta _{j'i}}$, then $(\textbf{y}_{ji}\otimes \textbf{x}_{i})^*(\textbf{y}_{j'i}\otimes \textbf{x}_{i})=0$, since $K_2$ is normal. The eigen-decomposition of $K_3$ can be similarly found. $\square$

Proof of Theorem 5

By Lemma 4, T consists of n orthonormal columns, hence is unitary. Moreover, (21) can be directly verified using (9), (10a), (11a), Lemmas 3, 4 and the properties of the Kronecker product. $\square$

Computating the WSVD of $\Phi _m$

Here, we briefly describe the numerical procedures to compute the WSVD of $\Phi _m$ in (24b). Specifically, given a HPD matrix $\Psi \!\in \!\mathbb {C}^{6\times 6}$, after multiplying $\Phi _m$ with $\Psi ^{1/2}$, we compute the SVD of $\Phi _m\Psi ^{1/2}$, i.e., 

$$\begin{aligned} \Phi _m\Psi ^{1/2}\!=\!P_m\texttt{diag}(\sigma _{m,1},\sigma _{m,2},\sigma _{m,3})\widetilde{Q}_m^*,\;\text{ with }\; P_m^*P_m\!=\!I_3\!=\!\widetilde{Q}_m^*\widetilde{Q}_m, \end{aligned}$$

(D3)

then we compute $Q_m:=\Psi ^{-1/2}{\widetilde{Q}}_m$. Now, it is easy to see that $Q_m$ satisfies $Q_m^{*}\Psi Q_m=I_3$ and $\Phi _m=P_m\texttt{diag}(\sigma _{1,m},\sigma _{2,m},\sigma _{3,m})Q_m^{*}$. In other words, the WSVD (25) of $\Phi _m$ is obtained. In particular, if $\Psi$ is set to $(\lambda _0 {\textbf{1}}_{3\times 3}+2\mu _0I_3)\oplus \mu _0 I_3$ following (38b), then its HPD square root and inverse square root are

$$\begin{aligned} \Psi ^{1/2}&=\left( (\sqrt{3\lambda _0+2\mu _0}-\sqrt{2\mu _0}){\textbf{1}}_{3\times 3}+ I_3\sqrt{2\mu _0}\right) \oplus I_3\sqrt{\mu _0}, \end{aligned}$$

(D4a)

$$\begin{aligned} \Psi ^{-1/2}&=\left( (1/\sqrt{3\lambda _0+2\mu _0}-1/\sqrt{2\mu _0}){\textbf{1}}_{3\times 3}+I_3/\sqrt{2\mu _0} \right) \oplus I_3/\sqrt{\mu _0}, \end{aligned}$$

(D4b)

respectively.

Proof of Corollary 11

Recall that it is claimed in Corollary 11 that $\Psi$ in (38b) along with (36b) changes the RHS of (31) into (39), i.e., 

$$\begin{aligned} \kappa ({\mathcal {W}}^{-1/2}(\mathcal {V}_v\oplus \mathcal {V}_f)\mathcal W^{-1/2})\!=\!\frac{\max \!\left( \!\left\{ \frac{3\lambda _\ell +2\mu _\ell }{3\lambda _0+2\mu _0}\!\right\} _{\ell =1}^2\!\bigcup \!\left\{ \frac{\mu _\ell }{\mu _0}\!\right\} _{\ell =1}^2 \right) }{\min \!\left( \!\left\{ \frac{3\lambda _\ell +2\mu _\ell }{3\lambda _0+2\mu _0}\!\right\} _{\ell =1}^2\!\bigcup \!\left\{ \frac{\mu _\ell }{\mu _0}\!\right\} _{\ell =1}^2\right) }<\frac{\kappa (\mathcal {V}_v\oplus \mathcal {V}_f)}{2}, \end{aligned}$$

which is proved as follows.

Proof

Denote $\Psi _{3\times 3}:=\lambda _0 {\textbf{1}}_{3\times 3}+2\mu _0I_3$. Following (33), eigenvalues of $\mathcal W^{-1/2}(\mathcal {V}_v\oplus \mathcal {V}_f){\mathcal {W}}^{-1/2}$ are $\mu _1/\mu _0, \mu _2/\mu _0$ and union of eigenvalues $\beta _m$ of $\Psi _{3\times 3}^{-1/2}(L_{v,mm}{\textbf{1}}_{3\times 3}+2M_{v,mm} I_3)\Psi _{3\times 3}^{-1/2}$ for $m=1,\ldots ,n$. It is easy to see

$$\begin{aligned} 0&=\det \left( (L_{v,mm}{\textbf{1}}_{3\times 3}+2M_{v,mm} I_3)-\beta _m\Psi _{3\times 3}\right) \\&=4(M_{v,mm} - \beta _m\mu _0)^2(3L_{v,mm}+2M_{v,mm}-\beta _m(3\lambda _0+2\mu _0)), \end{aligned}$$

thus $\beta _m= \tfrac{3L_{v,mm}+2M_{v,mm}}{3\lambda _0+2\mu _0},\; \tfrac{M_{v,mm}}{\mu _0}$, where $L_{v,mm},M_{v,mm}$ can be either $\lambda _1,\mu _1$ or $\lambda _2,\mu _2$. Hence, the equality in (39) holds by definition.

Without loss of generality, let $\mu _1<\mu _2$, which, in view of (35), implies

$$\begin{aligned} \kappa (\mathcal {V}_v\oplus \mathcal {V}_f)\ge \tfrac{3\lambda _1+2\mu _1}{\mu _1}, \quad \kappa (\mathcal {V}_v\oplus \mathcal {V}_f)\ge \tfrac{3\lambda _2+2\mu _2}{\mu _1}. \end{aligned}$$

(E5)

If $\mu _0=\mu _1$ and $\lambda _0=\lambda _1$, we have $\mu _2 / \mu _0 = \mu _2/\mu _1>1$, $(3\lambda _2+2\mu _2)/(3\lambda _0+2\mu _0)=(3\lambda _2+2\mu _2)/(3\lambda _1+2\mu _1)$ and $\mu _1/\mu _0=1=(3\lambda _1+2\mu _1)/(3\lambda _0+2\mu _0)$. Then $\mu _\ell /\mu _0$ and $(3\lambda _\ell +2\mu _\ell )/(3\lambda _0+2\mu _0)$, $\ell =1,2$, can be arranged in three possible orders as follows:

• 1. $1\le (3\lambda _2+2\mu _2)/(3\lambda _1+2\mu _1)\le \mu _2/\mu _1$, which, with (E5), implies

  $$\begin{aligned} \kappa ({\mathcal {W}}^{-1/2}(\mathcal {V}_v\oplus \mathcal {V}_f)\mathcal W^{-1/2}) = \tfrac{\mu _2}{\mu _1}< \tfrac{1}{2}\left( \tfrac{2\mu _2}{\mu _1}+\tfrac{3\lambda _2}{\mu _1}\right) \le \tfrac{\kappa (\mathcal {V}_v\oplus \mathcal {V}_f)}{2}; \end{aligned}$$

• 2. $1<\mu _2/\mu _1\le (3\lambda _2+2\mu _2)/(3\lambda _1+2\mu _1)$, which, with (E5), implies

  $$\begin{aligned} \kappa ({\mathcal {W}}^{-1/2}(\mathcal {V}_v\oplus \mathcal {V}_f)\mathcal W^{-1/2})= \tfrac{3\lambda _2+2\mu _2}{3\lambda _1+2\mu _1}< \tfrac{1}{2}\tfrac{3\lambda _2+2\mu _2}{\mu _1}\le \tfrac{\kappa (\mathcal {V}_v\oplus \mathcal {V}_f)}{2}; \end{aligned}$$

• 3. $(3\lambda _2+2\mu _2)/(3\lambda _1+2\mu _1)\le 1<\mu _2/\mu _1$, which, with (E5), implies

  $$\begin{aligned} \kappa ({\mathcal {W}}^{-1/2}(\mathcal {V}_v\oplus \mathcal {V}_f){\mathcal {W}}^{-1/2})&= \tfrac{\mu _2}{\mu _1}\tfrac{3\lambda _1+2\mu _1}{3\lambda _2+2\mu _2}=\left( \tfrac{3\lambda _2+2\mu _2}{\mu _2}\right) ^{-1}\tfrac{3\lambda _1+2\mu _1}{\mu _1} \\&<\tfrac{1}{2}\tfrac{3\lambda _1+2\mu _1}{\mu _1}\le \tfrac{\kappa (\mathcal {V}_v\oplus \mathcal {V}_f)}{2}. \end{aligned}$$

On the other hand, if $\mu _0=\mu _2$ and $\lambda _0=\lambda _2$, we can similarly show that $\kappa (\mathcal A)<\kappa (\mathcal {V}_v\oplus \mathcal {V}_f)/2$. Therefore, the inequality in (39) holds. $\square$