Exact domain truncation for the Morse-Ingard equations¹¹1This work was supported by NSF SHF-1909176 and SHF-1911019.

Laser absorption spectroscopy is used for detecting trace amounts of gases in diverse application areas such as air quality monitoring, disease diagnosis, and manufacturing [3, 4, 5]. In photoacoustic spectroscopy, a laser is fired between the tines of a small quartz tuning fork, and in the presence of a particular gas, acoustic and thermal waves are generated. These waves, in turn, interact with the tuning fork to generate an electric signal via pyroelectric and piezoelectric effects. Two variants of these sensors are the so-called QEPAS (quartz-enhanced photoacoustic spectroscopy) and ROTADE (resonant optothermoacoustic detection) models [6, 7]. In QEPAS, the acoustic wave dominates the signal, while the thermal wave is more important in ROTADE. In many experimental configurations, both effects appear. While a full model of the sensor is an eventual goal, obtaining efficient and accurate models of the gas itself is a critical step.

Earlier work on modeling this problem [8, 9, 10] simplified the model to a single PDE that included an empirically-determined damping term to account for otherwise-neglected processes. This approach is only accurate in particular regimes, and the empirical corrections depend strongly on geometry as well as physical parameters, which limits the model’s utility in an optimal design context.

Hence, work began on coupled models including both thermal and acoustic effects. A finite element discretization of the coupled pressure-temperature system of Morse and Ingard [1] was first addressed in [11], where the difficulty of solving the linear system was noted. Kirby and Brennan gave a more rigorous treatment in [12], with analysis of the finite element error and preconditioner performance. Kaderli et al derived an analytical solution for the coupled system in idealized geometry in [13]. Their technique involves reformulating the system studied in [12] by an algebraic simplification that eliminates the temperature Laplacian from the pressure equation. In [14], this reformulation was seen to lose coercivity but still retain a Gårding-type inequality, leading to optimal-order finite element convergence theory and preconditioners. Work by Safin et al [15] began a more robust multi-physics study, coupling the Morse-Ingard equations for atmospheric pressure and temperature to heat conduction of the quartz tuning fork, although vibrational effects were still not considered. They also applied a perfectly-matched layer (PML) [16] to truncate the computational domain, and a Schwarz-type preconditioner that separates out the PML region was used to effectively reduce the cost of solving the linear system. They also include some favorable comparisons between the computational model and experimental data.

Previous numerical analysis of this problem in the cited literature has focused on volumetric discretizations based on finite elements. In [17], we derived a boundary integral formulation for a scattered-field form of the Morse-Ingard equations. As with other wave problems, this problem writes the solution as the sum of a Morse-Ingard solution that satisfies the forcing (evaluated by means of a fast volumetric convolution with a Green’s function) plus a field that satisfies Morse-Ingard with no volumetric forcing but Neumann data on the tuning fork such that the sum satisfies homogeneous boundary conditions. We then formulated a second-kind integral equation for the scattered field and approximated it with a boundary integral method. In this work we return to finite element discretization, but we make use of the results we obtained considering the integral form of the equations to make significant advances in imposing a far-field condition.

In [2], we developed a novel nonlocal boundary condition for truncating the domain of Helmholtz scattering problems. This condition, which uses Green’s representation of the solution on the artificial boundary to give a nonlocal Robin-type condition involving layer potentials, is exact – the solution of resulting BVP agrees exactly with the restriction of the solution of the original problem to the computational domain. The variational form of the problem inherits a Gårding-type inequality from the Helmholtz operator so that a Galerkin finite element method yields optimal asymptotic convergence rates. Moreover, empirical results suggest that the standard local operator with transmission boundary conditions serves as an excellent preconditioner. Hence, one only needs to apply the action of the resulting nonlocal operator, say, by a fast multipole method, in order to obtain fast GMRES convergence.

In this paper, extend this approach from the Helmholtz operator to the Morse-Ingard system, making using of several results developed in our boundary-integral formulation in [17]. In Section 2, we recall the Morse-Ingard equations. Then, Section 3 addresses far-field boundary conditions for the system and appropriate boundary conditions for domain truncation. By means of the transformation to a decoupled Helmholtz system, we are able to state an analogous far-field condition and associated transmission-type condition for the Morse-Ingard system. This allows a comparison to the ad hoc transmission boundary conditions used in [12, 14]. Moreover, we can derive an exact analog of the nonlocal Helmholtz boundary condition for Morse-Ingard. Although one may directly solve the decoupled Helmholtz equations rather than the coupled form of Morse-Ingard, formulating boundary conditions and directly simulating the coupled system serves several purposes. First, the pointwise transformation, thought it only involves a $2\times 2$ matrix, is quite ill-conditioned and seems to limit the accuracy we obtain on fine meshes for the decoupled form compared to the coupled system. Second, a more complete model of trace gas sensors [18] involves coupling Morse-Ingard to the tuning fork vibration, which in turn requires modeling the fluid flow. Domain truncation will still be required, but coupling of pressure and temperature to the fluid and tuning fork may limit the utility of the decoupled system. Additionally, as noted in [17], solving for the acoustic mode while neglecting the thermal mode turns out to be an effective approximation. After developing the boundary conditions in Section 3, we derive a finite element formulation for the Morse-Ingard system in Section 4. We discuss the structure of the linear system and approaches to preconditioning in Section 5 and we then provide some numerical in Section 6 before offering some final conclusions in Section 7.

The Morse-Ingard equations of thermoacoustics are a system of partial differential equations for the temperature and pressure of an excited gas. The model begins from a time-domain formulation. After assuming time-periodic forcing and performing nondimensionalization and some algebraic manipulations, we arrive at the form given in [13] and further analyzed in [14]:

Here, $T$ and $P$ are the non-dimensional temperature and pressure, respectively, within the gas. $S$ is a volumetric forcing function, modeling for example a laser pulse. $\gamma$ is the ratio of specific heat of the gas at constant pressure to that at constant volume. The dimensionless number $\mathcal{M}$ measures the ratio of the product of the characteristic thermal conduction scale and forcing frequency to sound speed, and $\Lambda$ does similarly for the viscous length scale. Typical values of parameters

are taken as in [13, 17].

We let $\Omega^{c}\subset\mathbb{R}^{d}$ (with $d=2,3$) be a bounded domain representing the tuning fork, and let its boundary be called $\Gamma$. The complement of $\overline{\Omega^{c}}$ will be the domain $\Omega$ on which we pose (1). On $\Gamma$, we impose homogeneous Neumann boundary conditions,

which posits that the tuning fork is thermally insulated from the gas, and that the tuning fork is sound-hard. More advanced models, in which the gas heats the tuning fork or the acoustic waves couple to tuning fork deformation, generalize this condition [15].

A suitable far-field condition is required to close the model, which requires some appropriate decay at infinity akin to the Sommerfeld radiation condition for the Helmholtz operator. Numerical methods based on volumetric discretization on a truncated domain have posed either some kind of transmission-type condition [12, 14] or perfectly-matched layers [15].

In [17], we gave a boundary integral method for Morse-Ingard based on a scattered-field formulation, which turns the volumetric inhomogeneity into an inhomogeneous Neumann condition on $\Gamma$. We discuss this in greater detail in Section 3.2. To arrive at the scattered-field formulation, we split the solution into incoming and scattered waves via

where $T^{i}$ and $P^{i}$ satisfy (1) with the given forcing function $S$ but have some inhomogeneous boundary conditions on $\Gamma$. Then, $T^{s}$ and $P^{s}$ are chosen to satisfy (1) with homogeneous forcing $S=0$ and such that the combined waves satisfy (3). The incoming waves $T^{i}$ and $P^{i}$ can be constructed by volumetric convolution of $S$ with a free-space Green’s function. With these in hand, their normal derivatives on $\Gamma$ can be computed, and the negative of these used as boundary conditions for $T^{s}$ and $P^{s}$. Consequently, we drop the superscripts ‘s’ for the scattered field and, for the rest of the paper, consider the system of PDE

together with boundary conditions

on $\Gamma$ together with an appropriate far field condition. Although we obtained satisfactory results with a boundary integral method in [17], we work with volumetric (finite element) discretizations here to chart a path towards modeling of additional volumetric physical phenomena without additional computational machinery in future work. Consequently, in this paper we are primarily interested in developing an analog of the nonlocal boundary condition developed in [2] for the Morse-Ingard system. To this end, we introduce a truncating boundary $\Sigma$ and define a domain $\Omega^{\prime}\subset\Omega$ to be that contained between $\Gamma$ and $\Sigma$, as shown in Figure 1.

In [17], we demonstrated that the Morse-Ingard system (1) could be decoupled into a pair of independent Helmholtz equations. After significant algebraic manipulations, we introduce modified material coefficients

as well as separate thermal and acoustic wave numbers

With the change of variables

the Morse-Ingard system (1) decouples into separate Helmholtz equations

where $a_{t}$ and $a_{p}$ are data-dependent constants. For parameter regimes of interest, $k_{t}$ has a large imaginary part, so that $V_{t}$ rapidly attenuates. On the other hand, $k_{p}$ has a very small imaginary part so that $V_{p}$ attenuates slowly. Typical values of these new parameters based on those from given above are

The same change of variables decouples the scattered-field formulation (5). In this case, we have that

and, by linearity, we take the normal derivative of (9) on $\Gamma$ to find boundary conditions

So, with $g_{T}$ and $g_{P}$ given a priori, the scattered field formulation of Morse-Ingard can be solved as a pair of decoupled Helmholtz scattering problems.

In this section, we give a finite element formulation of the Morse-Ingard equations (5) under various boundary conditions. First, we establish some notation. We let $L^{2}(\Omega^{\prime})$ denote the standard space of complex-valued functions with moduli square-integrable over $\Omega^{\prime}$ and $H^{k}(\Omega^{\prime})\subset L^{2}(\Omega)$ the subspace consisting of functions with weak derivatives up to and including order $k$ also lying in $L^{2}(\Omega^{\prime})$. For any Banach space $\mathcal{V}$, we let $\|\cdot\|_{\mathcal{V}}$ denote its norm, with the subscript typically omitted when $\mathcal{V}=L^{2}(\Omega^{\prime})$.

The space $L^{2}(\Omega^{\prime})$ is equipped with the standard inner product

and we also define the inner product over a portion of the boundary $\tilde{\Gamma}\subset\partial\Omega^{\prime}$ by

We partition $\Omega^{\prime}$ into a family of conforming, quasi-uniform triangulations [23] $\{\mathcal{T}_{h}\}_{h>0}$. Let $\mathcal{V}_{h}$ be the standard space of continuous piecewise polynomials of some degree $k\geq 1$ over $\mathcal{T}_{h}$. Since we are dealing with a system of two PDE, we define $\boldsymbol{\mathcal{V}}_{h}=\mathcal{V}_{h}\times\mathcal{V}_{h}$.

We multiply the equations of (5) by the conjugate of the test functions $v,w\in\mathcal{V}_{h}$, respectively and integrate by parts over $\Omega^{\prime}$. We let $U=(T,P)$ and $\Psi=(v,w)$. Applying the Neumann boundary conditions (6) on $\Gamma$ but not taking action yet on $\Sigma$, we have

We add these equations together and define $a_{0}:\mathbf{V}\times\mathbf{V}\rightarrow\mathbb{C}$ as consisting of the volumetric terms:

and $F_{0}$ involving those boundary terms on the right-hand side:

Then, we can write (41) as

At this point, we can close the system by selecting any of the boundary conditions discussed in Section 3 and substituting in the relevant expressions for $\tfrac{\partial T}{\partial n}$ and $\tfrac{\partial P}{\partial n}$. Following the general form of the local boundary condition (24), we define $a_{A}$ by

with

for the respective $A$ chosen, cf. (25) and (26). This leads to the variational problem of finding $U\in\mathbf{V}$ such that

for all $\Psi\in\mathbf{V}$.

We now consider variational problems corresponding to the exact, but nonlocal, boundary conditions. Recall that these boundary conditions are parametrized over the choice of $\boldsymbol{\sigma}=(\sigma_{t},\sigma_{p})$. Using boundary condition (38) in (41) motivates defining the bilinear form

and linear form

Then, we pose the variational problem of finding $U^{\boldsymbol{\sigma}}\in\boldsymbol{\mathcal{V}}$ such that

for all $\Psi\in\boldsymbol{\mathcal{V}}$. In fact, for any choice of $\boldsymbol{\sigma}$, $U^{\boldsymbol{\sigma}}$ is the solution to (5) with scattering boundary conditions (6) and the Sommerfeld-type far field condition (19).

A standard Galerkin discretization of this problem is obtained by restricting the test function $\Psi$ to $\boldsymbol{\mathcal{V}}_{h}$, seeking $U^{\boldsymbol{\sigma}}_{h}\in\boldsymbol{\mathcal{V}}_{h}$ such that

for all $\Psi_{h}\in\boldsymbol{\mathcal{V}}_{h}$.

Analyzing this discretization follows along the lines proposed in [2] for the scalar Helmholtz problem – one establishes a Gårding inequality for the variational form, by which existing theory [23] for Galerkin methods for elliptic operators provides solvability and optimal $H^{1}$ and $L^{2}$ and error estimates, subject to a sufficiently fine mesh. We have proven a Gårding-type inequality for the local form of Morse-Ingard in [14], and the same techniques used to handle the nonlocal terms for Helmholtz in [2] can be used for Morse-Ingard. Consequently,

A major feature of our nonlocal boundary condition is the opportunity for efficient solvers. For the Helmholtz problem in [2], we demonstrated empirically that preconditioning the entire operator with the local part led to very low GMRES iteration counts. This, of course, means the cost of inverting the local part of the operator drives the overall cost. It is well known that high wave numbers lead to notorious difficulty for iterative methods [24]. Fortunately, this is not the case for the parameter regime of interest for Morse-Ingard. In decoupled form (10), the thermal mode has a large imaginary part to complement the large real part, while the thermal mode has wave number approximately 1. Standard multigrid algorithms handle both of these situations effectively [25], so solving the problem in decoupled form proceeds along lines given in [2], followed by forming $T$ and $P$ from $V_{t}$ and $V_{p}$.

For the fully coupled system, we introduce a basis $\left\{\psi_{i}\right\}_{i=1}^{\dim\boldsymbol{\mathcal{V}}_{h}}$, and then we can write (50) as a linear system

where

and then $U^{\boldsymbol{\sigma}}=\sum_{i=1}^{\dim\boldsymbol{\mathcal{V}}_{h}}\mathbf{x}_{i}\psi_{i}$.

Following the partition of $a_{\boldsymbol{\sigma}}$ given in (47), we can write $A=A^{L}+A^{NL}$, where

corresponding to the local and nonlocal parts of the bilinear form.

We can effectively solve the linear using preconditioned GMRES [26], which is a parameter-free algorithm approximating the solution of a linear system $A\mathbf{x}=\mathbf{b}$ as the element of the Krylov subspace $\operatorname{span}\{A^{i}\mathbf{b}\}_{i=0}^{m}$ minimizing the equation residual. Building the subspace does not require access to entries of $A$, just the action of $A$ on vectors. Unlike conjugate gradients, GMRES is not restricted to operators that are symmetric and positive definite.

For most problems arising in the discretization of PDE, GMRES is most frequently used in conjunction with a preconditioner. Mathematically, we multiply the linear system through by some matrix $\widehat{P}^{-1}$:

and so the Krylov space then is $\operatorname{span}\{\left(\widehat{P}^{-1}A\right)^{i}\widehat{P}^{-1}\mathbf{b}\}_{i=0}^{m}$.

The overall performance of GMRES typically depends on two factors – the cost of building and applying the operators $\widehat{P}^{-1}$ and $A$, and the total number of iterations. One hopes to obtain a per-application cost that scales linearly (or log-linearly) with respect to the number of unknowns in the linear system, and a total number of GMRES iterations that is bounded independently of the number of unknowns. We think of $\widehat{P}^{-1}$ being an approximation to the inverse of some matrix $P$ that approximates $A$. As with the Helmholtz problem in [2], we will take $P=A^{L}$, the local part of the operator. Then, applying $\widehat{P}^{-1}$ might correspond to applying the inverse of $P$ by a sparse direct method, an application of some block preconditioner such as that analyzed in [14], or some other strategy.

As a partial justification of our choice of preconditioning matrix, when $\widehat{P}=A^{L}$ so that the inverse is applied exactly, we arrive at a preconditioned matrix of the form

Because $A^{NL}$ discretizes a compact operator (layer potential in weak form) and, moreover, $\left(A^{L}\right)^{-1}$ discretizes the inverse of an elliptic operator, the preconditioned matrix has the form of a (discretization of) a compact perturbation of the identity. We suggested obtaining such a form via preconditioning as a heuristic in [27]. Moret [28] gives rigorous GMRES convergence estimates for this situation once one establishes certain bounds on the operators. In practice, one might not apply the inverse of $A^{L}$ exactly. Supposing it is spectrally equivalent, one still obtains a compact perturbation of a bounded operator. Recently, Blechta [29] has extended Moret’s results to describe GMRES convergence for this setting as well.

In this section, we apply our various formulations of the Morse-Ingard equations. All of our numerical experiments are conducted using Firedrake [30]. We generate meshes with Gmsh [31], using its internal interface to OpenCascade for constructive solid geometry. Our Firedrake experiments make use of the recently-developed ExternalOperator capability, which is a foreign function interface allowing users to incorporate operators implemented by other packages into UFL expressions [32]. This technology will be documented elsewhere.

We use this to enable Firedrake to express layer potentials and evaluate them using the Pytential package. Pytential [33] is an open-source, MIT licensed software system for evaluating layer potentials from source geometry represented by unstructured meshes with high accuracy and near-optimal complexity. Pytential provides for the discretization of a source surface using tools for high-order accurate nonsingular quadrature [34, 35], its refinement according to accuracy requirements [36], and, finally, the evaluation of integral operators via quadrature by expansion (QBX) [37] and the associated GIGAQBX fast algorithm [38], with rigorous accuracy guarantees in two and three dimensions [39]. This fast algorithm, can, in turn make use of FMMLIB [40, 41] for the evaluation of translation operators in the moderate-frequency regime for the Helmholtz operator. In our two-dimensional experiments, we use an FMM order of 15, which provides sufficient accuracy for the accuracy of layer potential evaluation to not limit the overall accuracy obtained. While the integrals in our variational problem do not require the singular integral technology provided in QBX, it does provide robustness in the case of $\Sigma$ and $\Gamma$ are chosen to lie close together.

Our simulations are performed on a Debian Linux machine using dual Broadwell-EP 12-core processors running at 2.20GHz. The machine has 256GB of 2133GHz DDR4 RAM. In these simulations, Firedrake and PETSc were using a single core, Pytential evaluates the layer potentials in multicore mode through the PoCL implementation of OpenCL.

We have developed exact truncating boundary conditions for the Morse-Ingard equations. These boundary conditions use a Green’s formula representation of the solution in terms of layer potentials and work in general unstructured geometry. The action of the discrete operators may be evaluated efficiently using matrix-free finite elements and a fast multipole method for the layer potentials, and the linear system may be effectively preconditioned with the local part of the operator. Standard convergence theory holds for the Galerkin discretization, and the method gives good accuracy on small computational domains even with relatively coarse meshes.

In the future, we hope to pursue a rigorous suite of three-dimensional calculations, compute with iterative treatment for $A^{L}$, especially as ongoing Pytential improves its performance for three-dimensional problems. We also hope to study models in which the Morse-Ingard are equations are coupled to the tuning fork displacement.