Computational design of acoustic materials using an adaptive optimization algorithm

TL;DR

This paper presents an adaptive optimization algorithm for designing acoustic materials with specific properties, utilizing local mesh refinement and symmetry considerations to improve efficiency and accuracy.

Contribution

It introduces a novel adaptive mesh refinement strategy combined with interpolation for the efficient design of symmetric acoustic structures.

Findings

Efficient construction of symmetric acoustic materials demonstrated in 2D.

Adaptive refinement improves design accuracy and computational efficiency.

Method can be extended to arbitrary geometries and higher dimensions.

Abstract

We consider the problem of the construction of the acoustic structure of arbitrary geometry with prescribed desired properties. We use optimization approach for the solution of this problem and minimize the Tikhonov functional on adaptively refined meshes. These meshes are refined locally only in places where the acoustic structure should be designed. Our special symmetric mesh refinement strategy together with interpolation procedure allows the construction of the symmetric acoustic material with prescribed properties. Efficiency of the presented adaptive optimization algorithm is illustrated on the construction of the symmetric acoustic material in two dimensions.

Equations

$$\{\begin{array}{ll}\tilde{c}(x)\frac{{\partial }^{2}u}{\partial t^2}-△u=0& \text{mbox}in{\mathbb{R}}^2\times (0,\infty ),\\ u(x,0)=f_0(x),u_t(x,0)=0& \text{mbox}inD.\end{array}$$

$$f_0\in H^1(D),\tilde{c}(x)\in C^2(D).$$

$$\begin{array}{rl}{H}_u^1(D_T)& :=\{w\in H^1(D_T):w(\cdot ,0)=0\},\\ H_{\lambda }^1(D_T)& :=\{w\in H^1(D_T):w(\cdot ,T)=0\},\\ U^1& =H_u^1(D_T)\times H_{\lambda }^1(D_T)\times C\left(\overline{D}\right).\end{array}$$

$$\{\begin{array}{ll}\tilde{c}\frac{{\partial }^{2}u}{\partial t^2}-△u=0& \text{mbox}in{D}_T,\\ u(x,0)=f_0(x),u_t(x,0)=0& \text{mbox}inD,\\ {\partial }_{n}u=p(x,t)& \text{mbox}on{S}_{1,1},\\ {\partial }_{n}u=-{\partial }_{t}u& \text{mbox}on{S}_{1,2}\cup S_2,\\ {\partial }_{n}u=0& \text{mbox}on{S}_3\cup S_4.\end{array}$$

$$\{\begin{array}{ll}\tilde{c}\left(x\right)\in \left[1,M\right],M=const.>0,& \text{ for }x\in G_1,\\ \tilde{c}(x)=1& \text{ for }x\in D_{FDM}\cup G_2.\end{array}$$

$$u\left(x,t\right)=\tilde{u}\left(x,t\right),\forall \left(x,t\right)\in S_T.$$

$$F(u,\tilde{c})=\frac{1}{2}{\int }_{S_T}(u-\tilde{u}{)}^2{z}_{\delta }(t)dSdt+\frac{1}{2}\gamma {\int }_{G_1}(\tilde{c}-{\tilde{c}}_0{)}^{2}dx,$$

$$\begin{array}{rl}L(v)& =F(u,\tilde{c})-{\int }_{D_T}\tilde{c}\frac{\partial \lambda }{\partial t}\frac{\partial u}{\partial t}dxdt+{\int }_{D_T}(\nabla u)(\nabla \lambda )dxdt\\ & -{\int }_{S_{1,1}}\lambda p(x,t)dSdt+{\int }_{S_{1,2}\cup S_2}\lambda {\partial }_{t}udSdt,\end{array}$$

$$L^{\prime }(v;\bar{v})=0,$$

$$L^{\prime }(v;\bar{v})=\frac{\partial L}{\partial \lambda }(v)(\bar{\lambda })+\frac{\partial L}{\partial u}(v)(\bar{u})+\frac{\partial L}{\partial \tilde{c}}(v)(\overline{\stackrel{~}{c}})=0,$$

$$\begin{array}{rl}0& =\frac{\partial L}{\partial \lambda }(v)(\bar{\lambda })=\\ & -{\int }_{D_T}\tilde{c}\frac{\partial \bar{\lambda }}{\partial t}\frac{\partial u}{\partial t}dxdt+{\int }_{D_T}(\nabla u)(\nabla \bar{\lambda })dxdt\\ & -{\int }_{S_{1,1}}\bar{\lambda }p(x,t)dSdt\\ & +{\int }_{S_{1,2}\cup S_2}\bar{\lambda }{\partial }_{t}udSdt,\forall \bar{\lambda }\in H_{\lambda }^1(D_T),\end{array}$$

$$\begin{array}{rl}0& =\frac{\partial L}{\partial u}(v)(\bar{u})=\\ & {\int }_{S_T}(u-\tilde{u})\bar{u}{z}_{\delta }dSdt-{\int }_D\tilde{c}\frac{\partial \lambda }{\partial t}(x,0)\bar{u}(x,0)dx\\ & -{\int }_{S_{1,2}\cup S_2}\frac{\partial \lambda }{\partial t}\bar{u}dSdt\\ & -{\int }_{D_T}\tilde{c}\frac{\partial \lambda }{\partial t}\frac{\partial \bar{u}}{\partial t}dxdt\\ & +{\int }_{D_T}(\nabla \lambda )(\nabla \bar{u})dxdt,\forall \bar{u}\in H_u^1(D_T),\end{array}$$

$$\begin{array}{rl}0& =\frac{\partial L}{\partial \tilde{c}}(v)(\overline{\stackrel{~}{c}})=-{\int }_{D_T}\frac{\partial \lambda }{\partial t}\frac{\partial u}{\partial t}\overline{\stackrel{~}{c}}dxdt\\ & +\gamma {\int }_{G_1}(\tilde{c}-{\tilde{c}}_0)\overline{\stackrel{~}{c}}dx,x\in D.\end{array}$$

$$\{\begin{array}{ll}\tilde{c}\frac{{\partial }^2\lambda }{\partial t^2}-△\lambda =-(u-\tilde{u})z_{\delta }& x\in S_T,\\ \lambda (\cdot ,T)=\frac{\partial \lambda }{\partial t}(\cdot ,T)=0,& \\ {\partial }_n\lambda ={\partial }_t\lambda & \text{mbox}on{S}_{1,2}\cup S_2,\\ {\partial }_n\lambda =0& \text{mbox}on{S}_3\cup S_4\cup S_{1,1}.\end{array}$$

$$F(u(\tilde{c}),\tilde{c})=L(v(\tilde{c})),$$

$$\begin{array}{rl}{F}^{\prime }(\tilde{c}):=& F^{\prime }(u(\tilde{c}),\tilde{c})=\frac{\partial F}{\partial \tilde{c}}(u(\tilde{c}),\tilde{c})=\frac{\partial L}{\partial \tilde{c}}(v(\tilde{c})).\end{array}$$

$$\begin{array}{rl}{F}^{\prime }(\tilde{c})(x)& :=F^{\prime }(u(\tilde{c}),\tilde{c})(x)=\\ & -{\int }_0^T\frac{\partial \lambda (\tilde{c})}{\partial t}\frac{\partial u(\tilde{c})}{\partial t}(x,t)dt+\gamma (\tilde{c}-{\tilde{c}}_0)(x).\end{array}$$

$$\begin{array}{rl}((a,b){)}_{D_T}& ={\int }_D{\int }_0^{T}abdxdt,\Vert a{\Vert }_{L_2(D_T)}^2=((a,a){)}_{D_T},\\ (a,b{)}_D& ={\int }_{D}abdx,\Vert a{\Vert }_{L_2(D)}^2=(a,a{)}_D.\end{array}$$

$$\begin{array}{rl}{\Vert \sqrt{\tilde{c}}{\partial }_{t}u(x,t)\Vert }_{L_2\left(D\right)}^2& +{\Vert \nabla u(x,t)\Vert }_{L_2\left(D\right)}^2\\ & \le A\left[{\Vert p(x,t)\Vert }_{L_2\left(D_t\right)}^2+{\Vert \nabla f_0\Vert }_{L_2\left(D\right)}^2\right].\end{array}$$

$$\begin{array}{r}{\Vert \sqrt{\tilde{c}}{\partial }_t\lambda (x,t)\Vert }_{L_2\left(D\right)}^2+{\Vert \nabla \lambda (x,t)\Vert }_{L_2\left(D\right)}^2\le B{\Vert (u-\tilde{u})z_{\delta }\Vert }_{L_2\left(D_{t_a}\right)}^2.\end{array}$$

$$\begin{array}{rl}& \underset{t}{\overset{T}{\int }}\underset{D}{\int }2\tilde{c}{\partial }_{tt}\lambda {\partial }_t\lambda dxd\tau -\underset{t}{\overset{T}{\int }}\underset{D}{\int }2\nabla \cdot (\nabla \lambda ){\partial }_t\lambda dxd\tau \\ & =-2\underset{t}{\overset{T}{\int }}\underset{{\partial }_{1}D\cup {\partial }_{2}D}{\int }(u-\tilde{u})z_{\delta }{\partial }_t\lambda dSd\tau .\end{array}$$

$$\begin{array}{r}\underset{t}{\overset{T}{\int }}\underset{D}{\int }{\partial }_t\left(\tilde{c}{\partial }_t{\lambda }^2\right)dxd\tau =-\underset{D}{\int }\left(\tilde{c}{\partial }_t{\lambda }^2\right)\left(x,t\right)dx.\end{array}$$

$$\begin{array}{rl}& 2\underset{t}{\overset{T}{\int }}\underset{D}{\int }\nabla \cdot \left(\nabla \lambda \right){\partial }_t\lambda dxd\tau =2\underset{t}{\overset{T}{\int }}\underset{\partial D}{\int }\left({\partial }_t\lambda \right){\partial }_n\lambda dSd\tau \\ & -2\underset{t}{\overset{T}{\int }}\underset{D}{\int }\left(\nabla \lambda \right)\left(\nabla {\partial }_t\lambda \right)dxd\tau \\ & =2\underset{t}{\overset{T}{\int }}\underset{{\partial }_{1}D\cup {\partial }_{2}D}{\int }{\left({\partial }_t\lambda \right)}^{2}dSd\tau -\underset{t}{\overset{T}{\int }}\underset{D}{\int }{\partial }_t|\nabla \lambda {|}^{2}dxd\tau .\end{array}$$

$$\begin{array}{rl}& \underset{t}{\overset{T}{\int }}\underset{D}{\int }{\partial }_t|\nabla \lambda {|}^{2}dxd\tau =\underset{D}{\int }|\nabla \lambda {|}^2\left(x,T\right)dx-\underset{D}{\int }|\nabla \lambda {|}^2\left(x,t\right)dx\\ & =-\underset{D}{\int }|\nabla \lambda {|}^2\left(x,t\right)dx.\end{array}$$

$$\begin{array}{rl}& -\underset{D}{\int }\left(\tilde{c}{\partial }_t{\lambda }^2\right)\left(x,t\right)dx-\underset{D}{\int }|\nabla \lambda {|}^2\left(x,t\right)dx\\ & =2\left(\underset{t}{\overset{T}{\int }}\underset{{\partial }_{1}D\cup {\partial }_{2}D}{\int }{\left({\partial }_t\lambda \right)}^2-(u-\tilde{u})z_{\delta }{\partial }_t\lambda \right)dSd\tau .\end{array}$$

$$\begin{array}{rl}& \underset{D}{\int }\left(\tilde{c}{\partial }_t{\lambda }^2\right)\left(x,t\right)dx+\underset{D}{\int }|\nabla \lambda {|}^2\left(x,t\right)dx\\ & =2\left(\underset{t}{\overset{T}{\int }}\underset{{\partial }_{1}D\cup {\partial }_{2}D}{\int }(u-\tilde{u})z_{\delta }{\partial }_t\lambda -{\left({\partial }_t\lambda \right)}^2\right)dSd\tau .\end{array}$$

$$\begin{array}{rl}& \underset{D}{\int }\left(\tilde{c}{\partial }_t{\lambda }^2+{\left|\nabla \lambda \right|}^2\right)(x,t)dx\\ & \le B\underset{t}{\overset{T}{\int }}\underset{{\partial }_{1}D\cup {\partial }_{2}D}{\int }|(\tilde{u}-u)z_{\delta }{|}^2(x,\tau )dSd\tau ,\end{array}$$

$$\begin{array}{rl}{W}_h^u:=& \{w\in H_u^1(Q_T):w{|}_{K\times J}\in P_1(K)\times P_1(J),\\ & \forall K\in K_h,\forall J\in J_{\tau }\},\end{array}$$

$$\begin{array}{rl}{W}_h^{\lambda }:=& \{w\in H_{\lambda }^1(Q_T):w{|}_{K\times J}\in P_1(K)\times P_1(J),\\ & \forall K\in K_h,\forall J\in J_{\tau }\}.\end{array}$$

$$C_h:=\{u\in L_2(D):u{|}_K\in P_0(K),\forall K\in K_h\},$$

Full text

1 Introduction

In this work we present a new adaptive optimization algorithm which can construct acoustic materials with arbitrary geometry from desired scattering parameters. We formulate our problem as a Coefficient Inverse Problem (CIP), and our goal is to determine an unknown spatially distributed wave speed of the acoustic wave equation from boundary measurements on the adaptively refined meshes. To solve our CIP, we minimize the Tikhonov functional in order to find the wave speed distribution inside designed domain which satisfies prescribed scattering properties. In the case of numerical simulations of Section 7 we formulate these properties as obtaining as small as possible reflections from the designed structure. For minimization of the Tikhonov functional we use Lagrangian approach and search for a stationary point of it on the adaptively refined meshes. Compared with other works on this subject [2, 4, 6] we need to refine mesh locally only inside the known geometry. For construction of a new mesh we use symmetric mesh refinement strategy combined with the interpolation procedure over the neighboring vertices for every element in the mesh. This allows us finally to get acoustic material of the symmetric structure.

To construct the desired acoustic structure we formulate an adaptive optimization algorithm which includes solution of the forward and adjoint problems for the acoustic wave equation. The domain decomposition finite element/finite difference (FE/FD) method of [3] is used for the computational solution of these problems. This method is implemented efficiently using the software packages WavES [19] and PETSc [16]. In the theoretical part of this work we present proof of the energy estimate for a hyperbolic equation with one unknown function - the wave speed- and different boundary conditions for the case of our domain decomposition. We illustrate efficiency of the proposed method in numerical examples on the construction of new acoustic material in two dimensions. The goal of our numerical simulations is to reconstruct the wave speed function of the hyperbolic equation from single observations of the solution of this equation in space and time which gives us as small reflections as possible. We note that the domain decomposition approach in this case is particularly feasible for implementing of absorbing boundary conditions [10].

Developed in this work adaptive optimization method can be used in construction and design of new materials including nano-materials with so-called cloaking properties, see [13, 14, 15]. To obtain cloaking structures in all these works are used methods of transformational optics which are based on the accordance between material parameters and coordinate transformations. In the current work we propose to use an adaptive optimization algorithm which is an alternative approach for the construction of an approximate cloaking. Depending on applications, this method can be used alone or as a compliment to the method of transformational optics. Advantage of a new technique compared to the transformational optics is fast construction of any material of arbitrary geometry with desired symmetric structure of any size. This structure is not dependent on the coordinate transformation and can be adapted to desired properties of the physical material. The mesh size of the symmetric structure can be defined as a parameter in the adaptive mesh refinement procedure used in the optimization algorithm. Thus, the new algorithm allows efficiently compute a new material of any symmetric structure with desired properties. A first version of a such algorithm was presented in [8] for design of a nanophotonic structure.

The paper is organized as follows. In Section 2 we present statements of the forward and inverse problems and in Section 3 we describe the Lagrangian approach for solution of our CIP. Stability estimates for the solution of forward and adjoint problems are given in Section 4. In Section 5 we present the domain decomposition FEM/FDM to solve the minimization problem of Section 3, and in Section 6 we present an adaptive conjugate gradient algorithm for the solution of our CIP. Finally, in our concluding Section 7 we demonstrate efficiency of the adaptive optimization algorithm identifying the wave speed function in two dimensions to construct material of symmetric structure which produce as small reflections as possible.

2 Statement of the forward and inverse problems

Let $x=(x_{1},x_{2})$ denote a point in $\mathbb{R}^{2}$ in an unbounded domain $D$. We model the wave propagation by the following Cauchy problem for the scalar wave equation:

[TABLE]

Here, $u$ is the total wave pressure generated by the plane wave $p(t)$ which is incident at $x_{1}=x_{0}$ and propagates along $x_{2}$ axis, $\tilde{c}(x)=\frac{1}{c(x)^{2}}$ is the isotropic function with the spatially distributed wave speed $c(x)$.

Let now $D\subset\mathbb{R}^{2}$ be a bounded domain with the boundary $\partial D$. We use the notation $D_{T}:=D\times(0,T),\partial D_{T}:=\partial D\times(0,T),T>0$ and assume that

[TABLE]

For computational solution of (1) we use the domain decomposition finite element/finite difference (FE/FD) method of [3] which was applied for the solution of different coefficient inverse problems for the acoustic wave equation in works [2, 3, 4, 5]. To apply method of [3] we decompose $D$ into two regions $D_{FEM}$ and $D_{FDM}$ such that the whole domain $D=D_{FEM}\cup D_{FDM}$, see Figure 1. In $D_{FEM}$ we use the finite element method (FEM), and in $D_{FDM}$ we will use the Finite Difference Method (FDM), see details in [3]. Furthermore, we decompose the domain $D_{FEM}$ into three regions $G_{0},G_{1},G_{2}$ such that $D_{FEM}=G_{0}\cup G_{1}\cup G_{2}$, where $G_{0}$ is the innermost subdomain with the boundary $\partial G_{0}$, $G_{1}$ is the subdomain where we want to design the acoustic material, and $G_{2}$ is the outermost subdomain, see Figure 1-b).

Let the boundary $\partial D$ be decomposed as $\partial D=\partial_{1}D\cup\partial_{2}D\cup\partial_{3}D$ where $\partial_{1}D$ and $\partial_{2}D$ are top and bottom sides of the domain $D$, respectively, and $\partial_{3}D$ is the union of left and right sides of this domain. At $S_{T}:=(\partial_{1}D\cup\partial_{2}D)\times(0,T)$ we have time-dependent observations. We define $S_{1}=\partial_{1}D\times(0,T)$, $S_{1,1}=\partial_{1}D\times(0,t_{1}]$, $S_{1,2}=\partial_{1}D\times(t_{1},T)$, $S_{2}=\partial_{2}D\times(0,T)$ and $S_{3}=\partial_{3}D\times(0,T)$, $S_{4}=\partial G_{0}\times(0,T)$.

We also introduce the following spaces of real valued functions

[TABLE]

In our computations we have used the following model problem

[TABLE]

In (4) we use the first order absorbing boundary conditions [10] and $p(x,t)\in L_{2}(S_{1,1})$. We note that these conditions are exact in the case of computations of Section 7, since in our computations we initialize the plane wave orthogonal to the domain of propagation.

We choose the coefficient $\tilde{c}(x)$ in (4) such that

[TABLE]

We consider the following inverse problem

Inverse Problem (IP)

*Let the coefficient $\tilde{c}\left(x\right)$ in the problem (4) satisfies conditions (5) and assume that $\tilde{c}\left(x\right)$ is unknown in the domain $G_{1}$. Determine the function $\tilde{c}\left(x\right)$ in (4) for $x\in G_{1}$ * assuming that the following function $\widetilde{u}\left(x,t\right)$ is known

[TABLE]

3 Optimization method

In this section we present the reconstruction method to solve inverse problem IP. This method is based on the finding of the stationary point of the following Tikhonov functional

[TABLE]

where $u$ satisfies the equations (4), $\tilde{c}_{0}$ is the initial guess for $\tilde{c}$ (see details about choice of this guess in Section 7 and [6, 7]), $\widetilde{u}$ is the observed field at $S_{T}$, $\gamma>0$ is the regularization parameter and $z_{\delta}(t)$ is the compatibility function in time and can be chosen as in [6].

To find minimum of (7) we use the Lagrangian approach [2, 4, 6] and define the following Lagrangian in the week form

[TABLE]

where $v=(u,\lambda,\tilde{c})\in U^{1}$, and search for a stationary point with respect to $v$ satisfying for all $\bar{v}=(\bar{u},\bar{\lambda},\bar{\tilde{c}})\in U^{1}$

[TABLE]

where $L^{\prime}(v;\cdot)$ is the Jacobian of $L$ at $v$.

In order to find the Fréchet derivative (9) of the Lagrangian (8) we consider $L(v+\bar{v})-L(v)~{}\forall\bar{v}\in U^{1}$ and single out the linear part of the obtained expression with respect to $\bar{v}$. When we derive the Fréchet derivative we assume that in the Lagrangian (8) functions in $v=(u,\lambda,\tilde{c})\in U^{1}$ can be varied independent on each others. We note that by doing so we get the same Fréchet derivative of the Lagrangian (8) as by assuming that functions $u$ and $\lambda$ are dependent on the coefficient $\tilde{c}$, see details in Chapter 4 of [7]. Similar to [2, 3, 4] we use conditions $\lambda\left(x,T\right)=\partial_{t}\lambda\left(x,T\right)=0$ and imply such conditions on the function $\lambda$ to deduce that $L\left(u,\lambda,\tilde{c}\right):=L\left(v\right)=F\left(u,\tilde{c}\right).$ We also use conditions (5) on $\partial D$, together with initial and boundary conditions of (4) to get that for all $\bar{v}\in U^{1}$ we have

[TABLE]

or

[TABLE]

[TABLE]

[TABLE]

We observe that (11) is the weak formulation of the state equation (4) and (12) is the weak formulation of the following adjoint problem

[TABLE]

We define by $u(\tilde{c}),\lambda(\tilde{c})$ exact solutions of the forward and adjoint problems, respectively, for the known function $\tilde{c}$. Then using the fact that exact solutions $u(\tilde{c}),\lambda(\tilde{c})$ are sufficiently stable (see Chapter 5 of book [12] for details), we get from (8)

[TABLE]

and the Fréchet derivative of the Tikhonov functional can be written as

[TABLE]

Inserting (13) into (16), we get the following space-dependent function:

[TABLE]

4 Stability estimates

The stability estimate for the forward problem (4) follows from the stability estimate of [3] and can be derived using the technique of [12]. For analysis we first introduce the $L_{2}$ inner product and the norm over $D_{T}$ and $D$, correspondingly, as

[TABLE]

Theorem

*Assume that the condition (5) for the function $\tilde{c}(x)$ holds. Let $D\subset\mathbb{R}^{n},n=2,3,$ be a bounded domain with a piecewise smooth boundary $\partial D$. * For any $t\in\left(0,T\right)$ we define $D_{t}=\partial_{1}D\times\left(0,t_{1}\right).$ * Assume that there exists a solution $u$ of the problem (4). Then $u\in H^{1}(D_{T})$ is unique and there exists a positive constant $A=A(\|\tilde{c}\|_{D},t)$ * such that the following energy estimate is true for all $t\in(0,T]$

[TABLE]

Proof.

A proof of this theorem follows from the stability estimate given in [3]. $\square$

The stability result for the adjoint problem is obtained similarly as for the forward problem, the only difference is in the integration in time $(t,T)$.

Theorem

*Assume that the condition (5) for the function $\tilde{c}(x)$ holds. Let $D\subset\mathbb{R}^{n},n=2,3$ be a bounded domain with a piecewise smooth boundary $\partial D$. * For any $t\in\left(0,T\right)$ we define by $D_{t_{a}}=(\partial_{1}D\cup\partial_{2}D)\times\left(t,T\right).$ * Assume that there exists a solution $\lambda$ of the problem (14) and a solution $u$ of the problem (4). Then $\lambda\in H^{1}(D_{T})$ is unique and there exists a positive constant $B=B(\|\tilde{c}\|_{D},t)$ * such that the following energy estimate is true for all $t\in(0,T]$

[TABLE]

Proof.

We multiply the equation in (14) by $2\partial_{t}\lambda$ and integrate over $D\times\left(t,T\right)$ to get

[TABLE]

Next, we integrate by parts in time the first term of (20) and noting zero initial condition in (14), we have

[TABLE]

Next, we integrate by parts in space the second term of (20). From (5) it follows that $\tilde{c}=1$ on $\partial D$. Thus, using (5) and absorbing boundary condition in (14), we get

[TABLE]

Integrating last term of (22) in time and using initial conditions of the equation (14), we obtain

[TABLE]

We insert (21)-(23) in (20) to get

[TABLE]

The equation above can be rewritten as

[TABLE]

Young’s inequality applied to (25) directly leads to

[TABLE]

with a constant $B=0.5$ which is the desired result.

$\square$

5 The finite element method in $D_{FEM}$

As was mentioned above for the numerical solution of (4) we use the domain decomposition FE/FD method of [3]. Similarly with this work, in our computations we decompose the finite difference domain $D_{FDM}$ into squares, and the finite element domain $D_{FEM}$ - into triangles. In $D_{FDM}$ we use the standard finite difference discretization of the equation (4) and obtain an explicit scheme as in [3].

For the finite element discretization of $D_{FEM}$ we define a partition $K_{h}=\{K\}$ which consists of triangles. We define by $h$ the mesh function as $h|_{K}=h_{K}$, where $h_{K}$ is the local diameter of the element $K$, and assume the minimal angle condition on the $K_{h}$ [9]. Let $J_{\tau}=\left\{J\right\}$ be a partition of the time interval $(0,\,T)$ into subintervals $J=(t_{k-1},\,t_{k}]$ of uniform length $\tau=t_{k}-t_{k-1}$.

To solve the state problem (4) and the adjoint problem (14) we define the finite element spaces, $W_{h}^{u}\subset H_{u}^{1}\left(Q_{T}\right)$ and $W_{h}^{\lambda}\subset H_{\lambda}^{1}\left(Q_{T}\right)$. First, we introduce the finite element trial space $W_{h}^{u}$

[TABLE]

where $P_{1}(K)$ and $P_{1}(J)$ denote the set of linear functions on $K$ and $J$, respectively. We also introduce the finite element test space $W_{h}^{\lambda}$ as

[TABLE]

To approximate the function $\tilde{c}$, we use the space of piecewise constant functions $C_{h}\subset L_{2}\left(D\right)$,

[TABLE]

where $P_{0}(K)$ is the set of constant functions on $K$.

Setting $V_{h}=W_{h}^{u}\times W_{h}^{\lambda}\times C_{h}$, the finite element method for (9) now reads: Find $v_{h}\in V_{h}$, such that

[TABLE]

To find approximate solution $v_{h}\in V_{h}$ we need to solve the forward problem (4), the adjoint problem (14) and then find the discrete gradient $L^{\prime}_{\tilde{c}}(v_{h})$. For the fully discrete schemes of these equations we refer to [3].

6 Adaptive conjugate gradient algorithm

To compute minimum of the functional (7) we use the adaptive conjugate gradient method (ACGM). The regularization parameter $\gamma$ in ACGM is computed iteratively via rules of [1]. For the local mesh refinement we use a posteriori error estimate of [2, 4] which means that the finite element mesh in $D_{FEM}$ should be locally refined where the maximum norm of the Fréchet derivative of the Lagrangian with respect to the coefficient is large. However, since our goal is to design material inside the known domain $G_{1}$, we refine mesh only inside this domain.

Now we define

[TABLE]

where $\tilde{c}_{h}^{m}$ is approximation of the function $\tilde{c}_{h}$ on the iteration step $m$ in AGCM, $E_{h}\left(x,t,\tilde{c}_{h}^{m}\right),\lambda_{h}\left(x,t,\tilde{c}_{h}^{m}\right)$ are computed by solving the state problem (4) and the adjoint problem (14), respectively, with $\tilde{c}:=\tilde{c}_{h}^{m}$. In our computations of section 7 we use the following algorithm.

Algorithm   (AGCM)

• •

Step 0. Set number of mesh refinements $j:=0$. Choose initial mesh $K_{h}^{j}$ in $D_{FEM}$ and time partition $J_{\tau}^{j}$ of the time interval $\left(0,T\right)$ as described in section 5. Start with the initial approximation $\tilde{c}_{h}^{0}=\tilde{c}_{0}$ at $K_{h}^{0}$ and compute the sequences of $\tilde{c}_{h}^{m}$ via the following steps:

• •

Step 1. Compute solutions $E_{h}\left(x,t,\tilde{c}_{h}^{m}\right)$ and $\lambda_{h}\left(x,t,\tilde{c}_{h}^{m}\right)$ of state (4) and adjoint (14) problems, respectively, on $K_{h}^{j}$ and $J_{\tau}^{j}$.

• •

Step 2. Update the coefficient $\tilde{c}_{h}:=\tilde{c}_{h}^{m+1}$ on $K_{h}^{j}$ (only inside the discretized domain $G_{1}$) and $J_{\tau}^{j}$ using the conjugate gradient method

[TABLE]

where

[TABLE]

with

[TABLE]

where $d^{0}(x)=-g^{0}(x)$. In (28) the step size $\alpha$ in the gradient update is computed as

[TABLE]

and the regularization parameter $\gamma^{m}$ at iteration $m$ is computed iteratively accordingly to [1] as

[TABLE]

• •

Step 3. Stop computing $\tilde{c}_{h}^{m}$ and obtain the function $\tilde{c}_{h}$ at $M=m$ if either $\|g^{m}\|_{L_{2}(D_{FEM})}\leq\theta$ or norms $\|g^{m}\|_{L_{2}(D_{FEM})}$ are stabilized. Here $\theta$ is the tolerance in updates $m$ of gradient method. Otherwise set $m:=m+1$ and go to step 1.

• •

Step 4. Refine the mesh $K_{h}^{j}$ inside $G_{1}$ using symmetric mesh refinement procedure, for example, as shown in Figure 5.

• •

Step 5. Set $j:=j+1$ and construct a new mesh $K_{h}^{j}$ in $D_{FEM}$ and a new partition $J_{\tau}^{j}$ of the time interval $\left(0,\,T\right)$ with the new time step $\tau$ which should be chosen correspondingly to the CFL condition of [20].

• •

Step 6. Interpolate the approximation $\tilde{c}_{h}$ computed on the step 3, from every element $K^{j-1}$ on the previous space mesh $K_{h}^{j-1}$ to the new elements $K^{j}$ in the mesh $K_{h}^{j}$, and obtain the initial guess $\tilde{c}_{0}$ on a new mesh. Set $m=1$ and return to step 1.

• •

Step 7. Stop refinements of $K_{h}^{j}$ and $J_{\tau}^{j}$ if norms defined in step 3 either increase or stabilize, compared to the previous space mesh.

7 Numerical Studies

The goal of this section is to present possibility of the computational design of an acoustic structure with the property to generate as small reflections as possible. This problem is equivalent to IP. Thus, we will reconstruct a function $\tilde{c}(x)$ inside a domain $G_{1}$ using the ACGM algorithm of section 6. We assume, that this function is known inside $D_{FDM}\cup G_{2}$ and is set to be $\tilde{c}(x)=1$.

Our computational geometry $D$ is split into two geometries $D_{FEM}$ and $D_{FDM}$ as described in section 2, see Figure 1. We denote by $\partial D_{FEM}$ the outer boundary of $D_{FEM}$ and by $\partial D_{FDM}$ the inner boundary of $D_{FDM}$. We set the dimensionless computational domain $D$ as

[TABLE]

and the domain $D_{FEM}$ as

[TABLE]

The spatial mesh in $D_{FEM}$ and in $D_{FDM}$ consists of triangles and squares, respectively. We choose the initial mesh size $h=0.02$ in $D=D_{FEM}\cup D_{FDM}$, as well as in the contiguous regions between FE/FD domains. We also decompose the domain $\mathrm{D}_{FEM}$ into three different domains $G_{0},G_{1},G_{2}$ such that $\mathrm{D}_{FEM}=G_{0}\cup G_{1}\cup G_{2}$ which are intersecting only by their boundaries, see Figure 1. The goal of our numerical tests is to reconstruct the function $\tilde{c}$ of the domain $G_{1}$ of Figure 1 which produces as small reflections as possible.

We initialize a plane wave in $D$ in time $T=[0,2.0]$ such that

[TABLE]

As for the forward problem in $D_{FDM}$ we solve the problem (4) choosing $\tilde{c}=1$, and in $D_{FEM}$ we solve

[TABLE]

Here, $\partial D_{{FDM}_{I}}$ denotes internal structured nodes of $D_{FDM}$ which have the same coordinates as structured nodes at the boundary $\partial D_{FEM}$, see details in [3]. We note, that we use the boundary condition $\partial_{n}u=0$ on $\partial G_{0}$ which implies that waves are not penetrated into $G_{0}$.

We also note that in $D_{FDM}$ the adjoint problem will be the following wave equation with $\tilde{c}(x)=1$ for $x\in D_{FDM}$:

[TABLE]

which we solve using finite difference method. In $D_{FEM}$ we solve the problem

[TABLE]

using finite element method. Here, $\partial D_{{FEM}_{I}}$ denotes internal structured nodes of $D_{FEM}$ lying on the inner boundary $\partial D_{FDM}$ of $D_{FDM}$, see details in [3] for the exchange procedure between FE/FD solutions.

As initial guess $\tilde{c}_{0}(x)$ we take different constant values of the function $\tilde{c}(x)$ inside domain of $G_{1}$ of Figure 1 on the coarse non-refined mesh, and we take $\tilde{c}(x)=1.0$ everywhere else in $D$. We choose three different constant values of $\tilde{c}_{0}(x)=\{1.5,2.0,2.5\}$ inside $G_{1}$. We define that the minimal and maximal values of the function $\tilde{c}(x)$ belong to the following set $M_{\tilde{c}}$ of admissible parameters

[TABLE]

The time step is chosen to be $\tau=0.002$ which satisfies the CFL condition [20].

7.1 Reconstructions

We generate data at the observation points at $S_{T}$ by solving the forward problem (4) in the time interval $t=[0,2.0]$, with function $p(t)$ given by (31) and for different values of $\omega=\{40,60,80,100\}$. To generate non-reflected data $\tilde{u}$ at $S_{T}$ we take the function $\tilde{c}(x)=1$ for all $x$ in $D$ and solve the problem (4) with a plane wave (31) and $\omega=\{40,60,80,100\}$. We regularize the solution of the inverse problem by starting computations with regularization parameter $\gamma=0.01$ in (7) and then updating this parameter iteratively in ACGM by formula (30). Computing the regularization parameter in this way is optimal for our problem. We refer to [11] for different techniques for choice of a regularization parameter.

Figure 2 shows real part of the Fourier transform of the time-dependent solution $u(x,t)$ of (4) when the initial guess for $\tilde{c}$ was $\tilde{c}_{0}=1.5$ in all points of $G_{1}$ (left figures), and after application of the adaptive optimization algorithm on three times refined mesh in $G_{1}$ (right figures) for different values of $\omega$ in (31). All right figures in Figure 2 show significant reduction of backscattered reflections for all tested frequencies compared with left figures.

Figures 3, 4 present reconstructions of $\tilde{c}$ which we have obtained on three time adaptively refined mesh inside the domain $G_{1}$ for different values of $\omega$ in (31). We note that different initial guesses $\tilde{c}_{0}$ in (7) produce different symmetric structures inside $G_{1}$ with different values of the function $\tilde{c}(x)$, compare reconstructions presented on Figures 5. Left images of Figure 5 present reconstructions obtained in ACGM when the optimized function $\tilde{c}$, obtained on a coarse mesh, is sequentially interpolated on the one, two and three times refined mesh. Then this interpolated function is taken as an initial guess $\tilde{c}_{0}$ in (7) and optimized further to get reconstruction on the third refined mesh. Right images of Figure 5 are obtained after direct application of the adaptive algorithm of Section 6. Optimized values of $\tilde{c}(x)$ obtained on Figures 3–5 can be of physical interest since they present symmetric structured domains with almost the same material in every structured layer.

Acknowledgments

The research of L.B. is supported by the sabbatical programme at the Faculty of Science, University of Gothenburg. The research of E.S. is supported by the Ministry of Education and Science of the Russian Federation, Project No. 1.894.2017/$\Pi$.