Isogeometric shape optimization for nonlinear ultrasound focusing

The goal of this work is to improve focusing of high-intensity ultrasound by modifying the geometry of acoustic lenses through shape optimization. The shape optimization problem is formulated by introducing a tracking-type cost functional to match a desired pressure distribution in the focal region. Westervelt's equation, a nonlinear acoustic wave equation, is used to model the pressure field. We apply the optimize first, then discretize approach, where we first rigorously compute the shape derivative of our cost functional. A gradient-based optimization algorithm is then developed within the concept of isogeometric analysis, where the geometry is exactly represented by splines at every gradient step and the same basis is used to approximate the equations. Numerical experiments in a 2D setting illustrate our findings.


Introduction
High-intensity focused ultrasound (HIFU) has become a powerful medical tool in recent years and a research area under active development. As a non-invasive therapeutic technology, it is used for breaking up kidney stones by lithotripsy [43,56,64] and it is currently involved in various stages of clinical trials for treatments of cancer in different organs, including the kidney, liver, prostate and brain [2,29,37,38,63].
Focusing in ultrasound applications is achieved geometrically, with lenses or spherical focusing arrays of piezoelectric transducers (see [36] and Figure 1 for geometry description). Such setups can naturally benefit from shape optimization: we can find the shape of the lens or the transducer surface which results in a desired pressure distribution around the focal point. The lens design has been quite often determined by heuristic laboratory-based approaches [46,49,65]. The goal of the present work is to improve the focusing of high-intensity ultrasound by employing a shape optimization approach to modify geometry of a focusing lens.
The use of high-power ultrasound sources enhances nonlinear effects in the wave propagation, inducing the formation of a sawtooth shape of the wave and generation of higher harmonics. When the focused acoustic wave propagates toward the focal point, its amplitude increases in a nonlinear fashion. Shape optimization for problems depending on linear partial differential equations is a well-explored topic; see, for instance, [24,66] for results concerned with analysis in the continuous setting and works [23,41,54], which tackle numerical aspects. Investigations of nonlinear problems are not as common, we refer, for example, to [21] for the treatment of an optimal design problem in nonlinear magnetostatics and [3] for optimization of instationary Navier-Stokes flow.
Results on gradient-based shape optimization for problems depending on second order hyperbolic equations are mainly restricted to the analysis with some numerical results available for linear problems. In [45], an optimal design problem governed by a linear damped one-dimensional wave equation was investigated. Results on optimal design for the support of the control for the 2D linear wave equation can be found in [44]. In [52], shape sensitivity analysis was provided for the problem of optimizing the shape of the acoustic lens; the case of ultrasound focusing by a concave array of transducers was treated analytically in [34]. First steps toward numerical shape optimization for linear ultrasound were made in [61]. Numerical treatment of the shape optimization problems for nonlinear ultrasound has to the authors' best knowledge been missing up to now.
Another goal of the present work is to investigate shape optimization subject to nonlinear waves within the concept of isogeometric analysis. Isogeometric analysis (IGA) was introduced by Hughes et al. in [9,27]. The idea behind IGA is that the non-uniform rational B-spline (NURBS) basis used to model geometry is also used as a basis for the desired solution. The idea of using spline functions for the approximation of partial differential equations can even be found in earlier works, see, e.g., [25] and references therein. However with IGA, geometry and analysis always go hand in hand, which is convenient especially for shape optimization. Numerical isogeometric shape optimization has been by now applied to a number of linear problems arising in fluid mechanics [53], electromagnetic scattering [51] and structural mechanics [6,39,62], while foundations of a mathematical framework for shape optimization problems within isogeometric analysis were laid out in [19,20]. In comparison to classical shape optimization methods, shape optimization based on the isogeometric approach offers several advantages. Even complex geometries can be represented precisely by multi-patch NURBS geometries and the output of the optimization can then directly be exported to a computer aided design (CAD) system. In addition, NURBS-bases can easily achieve higher order of global regularity than classical finite element bases.  The rest of the paper is organized as follows. In Section 2, we formulate the shape optimization problem together with its underlying nonlinear wave equation. Section 3 contains basics of shape calculus and provides a result on continuity of the state with respect to domain perturbations that will be needed to formally compute the derivative. Section 4 is then dedicated to the calculation of the shape derivative of the cost function. In Section 5, we give a brief overview of the isogeometric analysis and then introduce the numerical schemes we use to solve our state and adjoint problems. The shape optimization algorithm is presented in Section 6. Finally, in Section 7, we illustrate the efficiency of our algorithm through several numerical experiments.

Problem setting
To properly capture the nonlinear propagation behavior (see Figure 2) we have to employ a model of nonlinear acoustic for computing the pressure field. There is a number of so-called weakly nonlinear models in nonlinear acoustics obtained as a one-equation approximation of the compressible Navier-Stokes system. In the present work, we will employ a classical model, the Westervelt equation, to compute the acoustic pressure u:ü − c 2 ∆u − b∆u = 2k(u 2 + uü).
In the equation, c denotes the speed of sound in the medium and b the diffusivity. The strong b-damping in the model accounts for the energy losses that occur due to viscosity and thermal conductivity of the medium. The coefficient k is given by where is the mass density. Here the parameter of nonlinearity B/A accounts for the nonlinear pressure-density relation up to second order and indicates the strength of the nonlinearity of a medium. Values of B/A for different nonlinear fluids and gases are given, for example, in [59]. A detailed derivation of the model from the governing equations and an overview of well-established acoustic models can be found, for instance, in [10,36]. Typically in ultrasound applications the acoustic lens is immersed in an acoustic fluid. Let Ω ⊂ R d , d ∈ {1, 2, 3}, be a fixed bounded domain with Lipschitz boundary ∂Ω and Ω l its subdomain, representing the lens, such that Ω l ⊂ Ω and Ω l has Lipschitz boundary ∂Ω l (see Figure 3). We denote by Ω f = Ω \ Ω l the part of the domain representing the fluid region. Γ = ∂Ω l ∩ ∂Ω f then represents the interface between Ω l and Ω f .
The magneto-mechanical excitation will be modeled by Neumann boundary conditions on the part of the boundary ∂Ω denoted Γ n : where n represents the outer unit normal to ∂Ω. An important issue in acoustics is how to truncate an open domain problem to a bounded domain Ω ⊂ R d , d = 1, 2, 3, without causing spurious reflections of the waves. In this work we will employ the zero order Enquist-Majda absorbing boundary conditions [13]:u where Γ a = ∂Ω \ Γ n denotes the part of the boundary with absorbing conditions prescribed.
In what follows, we will denote by n l the unit outer normal to the lens domain Ω l . Restriction of a function v to Ω i will be denoted by v i , where i ∈ {l, f }, and v = v l − v f will stand for the jump over the interface Γ. The acoustic pressure can be calculated via the following initial boundary value problem: The coefficients are assumed to be piecewise constant functions, i.e. (2) where χ Ω l is the indicator function of Ω l . Throughout the paper, we assume that t ∈ [0, T ], with a finite time horizon T . In order to simplify the presentation and in accordance with ultrasound applications, we have set the initial conditions to zero. We will study the variational form of (1), given by {üφ + c 2 ∇u · ∇φ + b∇u · ∇φ − 2k(u 2 + uü)φ} dx ds for all test functions φ ∈ X = L 2 (0, T ; H 1 (Ω)) and (u,u)| t=0 = (0, 0).
2.1. On the well-posedness of (3). Westervelt's equation is a strongly damped quasilinear wave equation. An interesting feature of the model becomes more evident when rewritten as If u = 1/(2k), the equation degenerates, which means that special attention in the analysis has to be paid to obtaining a bound in the L ∞ norm for the pressure to keep it away from 1/(2k). Typically, this is resolved by a higher-regularity result and an embedding, e.g. H 2 (Ω) → L ∞ (Ω), together with appropriate energy estimates. For details we refer to various papers [32,33,34,47] dealing with well-posedness of this model with constant coefficients under different boundary conditions. Going forward, we assume that problem (3) is well-posed locally in time and space: Assumption 1. We assume that there existsm > 0 and a final time T > 0 such that under sufficiently regular boundary data g problem (3) has a unique solution .
Well-posedness of the problem (3) with constant coefficients c, b, k > 0 was shown to hold in [34,Theorem 2], with a solution u ∈ U ∩ L ∞ (0, T ; H 3/2+s (Ω)),s ∈ (0, 1/2), however well-posedness in the case of piecewise-constant coefficients is to the authors' best knowledge still an open problem. The main difficulty is that the H 3/2+s regularity in space, which is employed to avoid degeneracy, is too high to achieve when the normal derivative of the pressure is allowed to jump, as is the case in our problem.
Note that the last condition in (4) is imposed to ensure that the degeneracy of the model is avoided. In practice, this condition is easily fulfilled since the coefficient k is quite small. For instance, in water k ≈ 1.5 · 10 −9 s 2 m/kg, whereas the maximal pressure values in HIFU applications typically do not exceed 80 MPa [5,36,49]. This will be reflected in our numerical experiments as well.
2.2. The shape optimization problem. We will consider a shape optimization problem of the following form with U defined as in (4). For the objective function, we choose an L 2 -tracking type functional with D being a fixed bounded Lipschitz domain, D ⊂ Ω f , where the pressure distributions should match (see Figure 3). The constraint is given by We assume that the desired pressure u d ∈ L 2 (0, T ; L 2 (Ω)) ∩ H 1 (0, T ; (H 1 (Ω)) ), where (H 1 (Ω)) denotes the dual space of H 1 (Ω). Together with the assumed regularity of the state, we have u − u d ∈ L 2 (0, T ; L 2 (Ω)).

Remark 1.
From the point of view of applications it is crucial that the target pressure is reached around the focal point, which is why we demand that the acoustic pressure achieves prescribed values in a fixed subregion D. This means that our cost functional is compactly supported. The class of shape optimization problems with an L 2 (D)-tracking cost functional subject to Poisson equation with either Dirichlet or Neumann data are well-investigated and known to be ill-posed [14,15,22]. The problems were proven to be unstable by considering the shape Hessian and showing that it is not strictly H 1/2 -coercive at the critical domain. In our setting, however, the fixed domain D does not intersect the domain Ω l subject to optimization. The numerical experiments we conducted were shown to be effective without any regularization, which is why we do not introduce it here.
2.3. The adjoint problem. Since we will use the adjoint approach to compute the shape derivative, here we investigate the adjoint problem, given by the veryweak-in-time form The weak form of the adjoint problem is then given by for all ζ ∈ X = L 2 (0, T ; H 1 (Ω)), with (p,ṗ)| t=T = (0, 0).
Proposition 1. Let Assumption 1 on well-posedness of the state problem hold.
Then for sufficiently small final time T > 0 andm > 0, withm as in (4), problem (7) has a unique solution we can obtain a forward problem with initial conditions (p,ṗ)| t=0 = (0, 0). Well-posedness of a forward problem obtained in that way, but with coefficients assumed to be constant c, b, k > 0 was proven in [34,Corollary 7]. The proof is based on a fixed-point argument, together with Galerkin-Faedo approximations in space and lower and higher-order energy estimates. By closely inspecting the proof, we see that the well-posedness result still holds with our piecewise constant coefficients assumed to be as in (2).

Shape calculus
In order to define the shape derivative, we take the well-established approach of Murat and Simon [48], where changing domains are described by transformations of a reference domain. To this end, we introduce a fixed vector field Θ ∈ C 1,1 (Ω, R d ) with supp Θ ∩ (D ∪ ∂Ω) = ∅, and a family of transformations F τ : Ω → R d given by perturbations of the identity mapping There exists τ 0 > 0 such that for |τ | < τ 0 , F τ (Ω) = Ω and F τ is a C 1,1 diffeomorphism (cf. [30]). We define the perturbed lens as Ω l,τ = F τ (Ω l ) and the perturbed lens boundary as Then Γ τ is Lipschitz continuous (cf. [26,Theorem 4.1]) and Ω l,τ ⊂ Ω for |τ | < τ 0 . Due to assumptions on Θ, the boundary of Ω and subdomain D remain fixed as τ changes.
The Eulerian derivative of J at Ω l in the direction of the vector field Θ is defined as It is said that the functional J is shape differentiable at Ω l if dJ(u, Ω l )Θ exists for all Θ ∈ C 1,1 (Ω, R d ) and defines a continuous linear functional on C 1,1 (Ω, R d ).
Since the shape derivative defines a linear and bounded functional, according to the Riesz representation theorem it has a unique representative in a Hilbert space. This representative is commonly referred to as the shape gradient, i.e. the gradient is defined as the solution of (∇J, φ) H = dJ(u, Ω l )φ, ∀φ ∈ H.
This allows a certain freedom in choosing the Hilbert space H and in turn the descent direction. A comparison of L 2 and weighted H 1 scalar products defined on the boundary can be found in [12]. In the present work we will consider the common case H = L 2 (∂Ω l ).
For sufficiently small |τ |, we then bring back the transported solution to the domain with the fixed lens by introducing which will satisfy for all φ ∈ X, where we have employed the following notation with DF τ being the Jacobian of the transformation F τ . Note that we have made use of the fact that F τ vanishes on Γ n ∪ Γ a to end up with the boundary terms in (9).
In what is to follow we will rely on some useful properties of F τ , which we for the convenience of the reader recall here.
1. An auxiliary result on Hölder continuity. We will compute the shape derivative of our pressure-tracking shape functional by the so-called variational approach, introduced in [31]. This method avoids the condition of shape differentiability of the state variable by replacing it with Hölder continuity with respect to domain perturbations. Although the derivation of the volume representation of the shape derivative is analogous to the one presented in [52], we include it here for the sake of completeness. We remark however, that the computation of the boundary representation of the shape derivative differs from [52] due to different boundary conditions (influencing piece-wise regularity of the solution) and a lack of nonlinear damping in the present model, which allows for simplifications of the final derivative expression.
By employing the transformation of integrals formula it can be shown that u τ is uniformly bounded with respect to τ , i.e.
Theorem 3.2. Under Assumption 1, the solution of problem (3) is Hölder continuous with exponent 1/2 with respect to domain perturbations; i.e. it holds Proof. We begin by noting that the difference z τ = u τ − u satisfies , where the right hand side is given by Testing (13) by φ =ż τ χ [0,t] ∈ X then leads to In the last estimate, we have made use of the continuous embedding H 1 (Ω) → L 4 (Ω), cf. [11,Corollary 4.53], with the embedding constant C Ω H 1 ,L 4 . For further estimating the terms in the last line we employ (since z τ | t=0 =ż τ | t=0 = 0), to conclude that for sufficiently smallm and final time T it holds for some sufficiently large C > 0 that does not depend on τ . We can estimate the terms on the right hand side in (14) as follows whereā is defined as in (4) and α 2 as in (10). Recalling that u 0 = 0, we can further derive which, together with uniform boundedness of u τ in the sense of (11), leads to By combining (16) with estimate (14) and then employing the properties of the mapping F τ from Lemma 3.1, we first conclude that Secondly, we divide (14) by τ = 0, and then it remains to show that This now follows from estimate (16), Lemma 3.1 and (17).

Shape derivative of the cost function
Next, we will make use of the tools of shape calculus and the continuity result of the previous section to compute the shape derivative. We will begin by calculating the derivative of J in terms of the volume integral over the whole domain Ω, which is sometimes regarded as the weak shape derivative. Proposition 2. Let Assumption 1 on the well-posedness of the state problem hold and let u d ∈ L 2 (0, T ; L 2 (Ω)) ∩ H 1 (0, T ; (H 1 (Ω)) ). The volume representation of the shape derivative of J at Ω l in the direction of Θ ∈ C 1,1 (Ω, R d ), where supp Θ ∩ (D ∪ ∂Ω) = ∅, is given by Proof. By recalling the definition of the shape derivative and the fact that supp Θ ∩ D = ∅, together with the continuity result of Theorem 3.2, we arrive at Note that the last integral appears as the right hand side in the adjoint problem (7) when the test function is u τ −u. We employ again the notation z τ = u τ −u, recalling that z τ |t=0 =ż τ |t=0 = 0 and p |t=T =ṗ |t=T = 0. Testing the adjoint problem (7) with z τ ∈ L 2 (0, T ; H 1 (Ω)) and integrating by parts with respect to time leads to where in the last line we have then also used the weak form (13) satisfied by z τ . Due to the continuity result of Theorem 3.2 and our computations reduce to Thanks to Theorem 3.2 and Lemma 3.1, it is straightforward to verify that By employing Lemma 3.1 again, we can directly compute the limit on the right hand side of (19), which leads to the expression (2). 4.0.1. Imposing stronger regularity assumptions. According to the Hadamard-Zolésio Structure Theorem [66, Theorem 3.5], when J(u, Ω l ) and the domain are smooth enough, the shape derivative can be represented in the form To be able to express the derivative in terms of an integral over the boundary of the lens we have to impose stronger assumptions with respect to regularity of u, p, u d and domains Ω l , Ω f . From this point on let Ω l , Ω f be of class C 1,1 . Moreover, we introduce the following assumption: Assumption 2. We assume that the state and the adjoint variable are sufficiently regular in space on each of the subdomains, i.e. (20) u where M(Ω) is the linear space induced by the norm for somes ∈ (0, 1 2 ). Let us recall a useful identity that generalizes integration by parts in M(Ω), which we will employ to transform expression (18).
where n denotes the outer unit normal to ∂Ω.
We are now ready to state the main result.
where Γ denotes the interface between Ω l and Ω f .
Proof. Under Assumption 2, we can rewrite the expression (18) for the weak shape derivative as the sum of integrals over Ω l and Ω f and directly apply formula (21) to the terms containing ∇u · ∇p and ∇u · ∇p. The remaining terms containing div Θ can be dealt with by employing integration by parts with respect to time and space Here we have also made use of the fact that (u,u) |t=0 = (p,ṗ) |t=T = (0, 0). The last term in (23) can be further rewritten as In this way we transform (18) into The first two sums on the right-hand side vanish and we are left with It is possible to slightly simplify expression (24) Hence we have the following identity as well as By plugging these identities into (24) and making use of uuṗ = 0, the derivative reduces to dJ(u, Furthermore, due to conditions (25) on the interface and again the fact that (u,u) |t=0 = (0, 0) and (p, which helps us further to obtain Lastly, we can make use of the fact that again due to interface conditions (24) it holds which finally leads to (22).

Isogeometric analysis
In this section, we give a brief overview of basic concepts of isogeometric analysis which will be used to solve our state and adjoint problems and set the notation for future use. For an in-depth study of the isogeometric approach, we refer the reader to classical literature [9,57,60].  where l is the polynomial order and n is the number of basis functions used to construct the B-spline curve. As knot values are allowed to repeat, let us also introduce Z = {ζ 1 , . . . , ζ m } as a knot vector without any repetitions. Z partitions the parameter spaceΩ = [0, 1] into elements.
We will denote by B i,q , i = 1, . . . , d, B-Splines defined on Ξ recursivelŷ where in case one of the denominators is zero the contribution is assumed to be zero and omitted. The spline space can then be introduced as The NURBS geometry is then defined as Ω = G(Ω). We note that in isogeometric analysis one map G :Ω → Ω takes the patch from the parameter space to the physical space. However, from an implementational point of view, it still makes sense to speak about elements in which the domain gets divided through Z, since, for instance, quadrature nodes can then be defined in a uniform way for all the elements.

5.2.
Description of the computational lens/fluid domain in the multipatch setting. In our numerical experiments, we will assume that the problem is symmetric with respect to the y-axis through the center of the lens. Hence it suffices to only discretize half of the domain Ω and use symmetry boundary conditions at the y axis.
In practical examples computational domains are often described with several patches. In our case the full (half-)domain is decomposed into N ptc = 7 patches, Ω (j) , j = 1, . . . , N ptc , six of which represent water and one the lens region (see Figure 4). It is assumed that Ω = Nptc j=1 Ω (j) and Ω (j) are assumed to be disjoint. Every patch has one map G (j) , j ∈ {1, . . . , N ptc }, which maps the parametric domain to the physical one. Since the pressure u is continuous across the interface, but its normal derivative is not, we assume that the meshes and the control points coincide on the interface and impose strong C 0 -but not C 1 -continuity between patches. Furthermore, this assumption holds in every optimization step, i.e. for every updated geometry.

5.3.
Numerical treatment of the state problem. We will employ a semidiscrete method, where we discretize the space using a Galerkin finite element scheme and formulate the problem for continuous time. The main principle of isogeometric analysis is that the NURBS basis which we used to exactly model our domain also serves as the basis for the solution space of our numerical method. Following [1], in each patch Ω (j) we define a discrete space The discrete domain on the whole space will consist of continuous functions which, restricted on each domain, belong to V h . In other words, we seek an approximate solution u h such that Remark 2. One has to distinguish between a patch-wise and a global set of degrees of freedom. Within an individual patch the numerical solution has the form If two patches are "glued together", due to our assumptions on coinciding control points and meshes, it is possible to identify degrees of freedom with each other and also to add their respective ansatz functions.
For example, when using two rectangular patches of 3 × 3 degrees of freedom each and aligning them on one common boundary we would have to identify the three degrees of freedom of that boundary of patch 2 with the respective degrees of freedom of patch 1, giving a total of 15 instead of 9 + 9 = 18 global degrees of freedom. This approach ensures global continuity and hanging nodes are avoided. In this global set of ansatz functions one introduces a new (global) numbering of degrees of freedom and calls the set of all global ansatz functions (after identification on all interfaces) I such that the global numerical solution can be written in the following way: where the N i,q are the global ansatz functions (after identification).
Let us denote coefficients of a function u in the global NURBS basis by where time dependency is omitted for notational convenience. After discretizing both state and domain with NURBS, we arrive at the following semidiscrete formulation of our state problem: where mass, damping, and stiffness matrices are defined by 1.) Computing with M (j) , j = 1, ..., N ptc being the patch-wise mass matrices, defined by where DG (j) denotes the Jacobian, (DG (j) ) −T the inverse transposed of the patch-wise mappings G (j) and we have omitted polynomial degree q of basis functions to simplify notation.
2.) Secondly, identifying degrees of freedom at interfaces by adding the respective rows and columns in the global matrix M . For instance, at the interface between patch 1 and 2 we add the rows corresponding to interface degrees of freedom of patch 2 to the respective ones of patch 1 and then remove them from the global matrix. Proceeding analogously with the columns and repeating the process for all interfaces, we obtain M from M . This corresponds to transforming M via M = E T M E, where E is a rectangular matrix that only consists of ones and zeros.
We proceed in the same manner to obtain the damping matrix C and stiffness matrix K, with the only difference being their patch-wise definitions compared to the mass matrix, i.e.: Remark 3. Even though the isogeometric approach does not rely on individual elements and hence mass and stiffness matrices can analytically be calculated as above via integrals over the whole (parametric) domainΩ, it is still more convenient to adopt the notion of elements from classical FEA to the support boundaries of the individual NURBS basis functions, that are given by Z. This speeds up the computation as the integrals are only evaluated over parts of the domainΩ, where the respective basis functions do not vanish. In this sense one can also reformulate the definition of M (j) , C (j) and K (j) in terms of the assembly operator known from classical finite elements with n (j) e being the number of such elements per patch: In the same way also the tensor T is defined as the global structure after interface identification of degrees of freedom of the patch-wise tensors and the right-hand side vector and the matrices corresponding to the absorbing conditions are patch-wise given as which are then -as shown for M before -put together into global structures F, A 1 and A 2 by interface identification of degrees of freedom.

Time stepping scheme.
Due to the nonlinear nature of our problem, the higher harmonics of the fundamental frequency of the initial wave will grow with the distance from the source. In numerical algorithms cutting off higher harmonics is known to cause spurious oscillations at the shock peaks; this is often regarded as the Gibbs phenomenon. We tackle this issue by employing the generalized α-scheme for time stepping, introduced first by Chung and Hulbert in [7]. The generalized α-method allows to control the level of numerical dissipation of high frequencies with minimal unwanted damping in the low frequency range. In this way the high frequency modes which are insufficiently resolved by the spatial or the temporal discretization are eliminated.
Let 0 = t 0 < t 1 < . . . t m = T be a given uniform partition of the time domain and let us define the time step as ∆t = t n+1 − t n . The α-method applied to the semi-discrete equation (26) yields the modified equation Note that we have adopted a mid-point approach for the nonlinearities in (27). Morover, the classical Newmark relations [50] are used which will be numerically realized as a predictor-corrector algorithm. When α m = α f = 0, the time stepping reduces to the standard Newmark scheme. Relations among parameters α m , α f , β, and γ that guarantee stability and desired dissipation levels are given in [16] for systems with nonlinear internal forces involving u, but not its time derivatives. Adopting the strategy from [28,36], we define predictors for u n+1 andu n+1 as follows and the effective mass matrix which leads to the following system The system (28) is nonlinear, so we have to solve it by an iterative scheme. Following [36], we will employ a fixed-point iteration with respect to the second time derivative. Note that the effective mass matrix remains unchanged during the time stepping. An LU decomposition can therefore be performed once in the first time step, while in subsequent steps matrix inversion is avoided.
perform corrector step: Numerical treatment of the adjoint problem. For solving the adjoint problem, we will use the same NURBS space and the same time discretization as for the state solver. Let us denote by the coefficients of the target pressure u d in the NURBS base and by u d,n the coefficients at time step n. After discretizing also the adjoint state with NURBS, we obtain the semidiscrete formulation of our adjoint problem The matrices M, C, K and tensor T in (29) are defined in the same way as for the state problem. Note that the integral appearing on the right hand side in the weak formulation (7) of the adjoint problem only needs to be numerically evaluated on the focal region D. Therefore, in order to define the matrix M D , we first introduce where j = 1, ..., N ptc . Compared to M (j) , now we only integrate over the part of D where the pressure is tracked. This is expressed by the use of the characteristic function χ D taking the value 1 within D and 0 elsewhere. We define The matrix M D is then obtained from M D by again performing interface identification of degrees of freedom as it has been already done for the global mass, damping, and stiffness matrices.
Although the adjoint problem is linear, the modified mass matrix M − T (u) is time-dependent and thus would have to be inverted in every time step of our numerical solver. To avoid this computationally expensive step, the adjoint problem is also solved iteratively. The time stepping is performed backward in time by a Newmark predictor-corrector scheme to obtain p and thus p h as the solution. The Newmark parameters β p and γ p for the adjoint problem are chosen in such a way that second order of accuracy is achieved.
Algorithm 2: Solving the adjoint problem 1 t = T , n = number of time steps 2 p n = 0,ṗ n = 0,p n = 0 3 compute the effective mass matrix:

Shape optimization algorithm
Having seen how to numerically solve the state and the adjoint problem, we can now move to the general form of our optimization algorithm with step-size control. The optimality condition in line 3 will in practice be realized by testing if the norm of the shape gradient at the current step is sufficiently small, i.e. below a certain tolerance level TOL grad , relative to the gradient norm at the first step. We also introduce s max as the maximal number of gradient steps.
Lines 8 to 21 in Algorithm 4 describe the step size control, where the new geometry is accepted if it is better than the previous one in the sense of the value of the cost function. If not, the update is repeated with a smaller descent step. We introduce a lower step size bound TOL step so that the search for a smaller step can be stopped after a certain number of attempts when the cost and gradient norm don't decrease significantly anymore. However, if the newly obtained geometry is better than the previous one and there was no repetition of the previous step, we allow the current descent step to become larger. 6.1. Updating geometry. Let us reflect on step 7 of Algorithm 3 where the geometry is updated. Within the isogeometric approach shapes are easily changed by moving control points. In our solver, we restrict ourselves to a 2D setting, where control points corresponding to the lens boundary are treated as design variables. We denote by I ∂Ω l the set of indices of degrees of freedom contributing to the boundary of the lens, i.e.
In every optimization step, control points C i = (x i , y i ), i ∈ I ∂Ω l , are moved in the y-direction (while they remain fixed in the x-direction), according to where α is the current descent step and e y = (0, 1) T . Since we identify degrees of freedom on the interfaces, the fluid boundaries common with the lens ones are moved in the same way, whereas the others remain fixed.
In (30), we emphasized the dependence of the shape gradient on the solution of the adjoint problem. When u h is the exact solution of the state and p h of the adjoint problem, then dJ(u h , p h , Ω l )Θ = dJ(u, Ω l )Θ.
6.1.1. Extending the boundary changes to the interior. After updating the lens boundary, the internal control points of the lens and other patches are recomputed using the Coons patch approach [58].
Let us denote the updated control points by C new The updated boundary curves of, for instance, the lens patch are given by each mapping the parametric interval [0, 1] to the respective boundary in the physical space. A Coons surface is then computed as a combination of the following three surfaces: • Two ruled NURBS surfaces by linear interpolation of opposite boundary curves: • The bilinear NURBS surface interpolating between the corner points: The new patch mapping G new is then defined as For simplicity, we have omitted the index (j) of the current patch in notation above. For instance, G is used instead of G (j) to denote the current NURBS mapping. We note that the construction of the Coons patch can be motivated purely geometrically and is hence independent of the chosen coordinate frame.

Remark 4.
In the present work we do not employ tools to avoid mesh tangling, since in our numerical experiments this unwanted effect was never present. One possibility would be to use relative positioning when control points are moved only in the y-direction as suggested in [19].

Numerical experiments
We will now test our algorithm in two different scenarios. In the first one, we validate our code by creating artificial pressure data with a known goal shape for the lens and then trying to get back to that goal shape by starting the optimization process with a different lens geometry. In the second scenario, the goal shape will not be known and the target pressure distribution will be given analytically.
In both cases, we will perform experiments with linear B-splines (where IGA and classical FEM coincide) as well as quadratic splines. In all the experiments we will keep the lens corner point fixed, i.e. the shape will be pinned at the right corner.
The numerical results presented here have been obtained with the help of the GeoPDEs package in MATLAB [17].

Coefficient values.
For all the performed tests, the fluid is taken to be water. The values of physical parameters in typical media can be found, for instance, in [18,59]. The coefficients used in the experiments are provided in Table 1. In this way our waveform resembles the short tone bursts used in medical ultrasound practices, while still being continuous in time and space. We set the source amplitude value to be g 0 = 4 · 10 9 Pa and frequency f = 70kHz. The angular frequency is then given by ω = 2πf . Figure 5 shows the propagation of the signal in an exemplary lens setting. Toward the end of the geometry the steepening due to the nonlinear effects can be observed, which is also where the focal point is assumed to be and where high pressure values are required. Since the analytical solution of our problem is not known, we will generate artificial data for testing purposes. Using the same solver to generate data and then to invert the problem is known as committing the "inverse crime" [8], since it can cause unrealistically good solutions. To avoid the inverse crime, we produce the synthetic pressure data on a domain with lens shape Ω l by solving the forward problem on a staggered grid (compared to the grid that is used for the optimization process) and pollute it with random Gaussian noise of 2 % of the maximal pressure amplitude. Ω l is then taken to be the exact lens shape and the artificially produced data the target pressure distribution u d .

Domain measures.
In the forthcoming numerical experiments we will use several lens shapes as initial Ω 0 l and goal domains Ω l . Table 2 lists all the specific measurements.  7.3.2. Spatial discretization. Table 3 specifies the number of degrees of freedom used in the numerical experiments. Here ndof x denotes the number of degrees of freedom in x-direction and ndof y in y-direction of the computational domain. The overall number of degrees of freedom is denoted by ndof. Such a fine mesh is necessary to cope with the nonlinear effects.
In the focal region, it amounts to about 35 elements in y-direction per fundamental wavelength. With this refinement, we have 36 design control points for the upper as well as for the lower lens boundary.
7.3.3. Time stepping. The time resolution as well as the generalized α-parameters used when solving the state problem are given in Table 4. For the adjoint solver only the parameters γ p and β p are relevant.  This is a thinner shape than the initial one, so the upper boundary should move downward in order to reach its goal. The exact domain measurements are given in Table 2.
Let us introduce the discrete shape gradient ∇J(Ω s l ), where s denotes the current gradient step, as containing the discrete sensitivities of all boundary control points. We are interested in the relative change of the cost functional J(Ω s l )/J(Ω 0 l ) and the relative change of the norm of the shape gradient ∇J(Ω s l ) / ∇J(Ω 0 l ) over the course of optimization.
A scaling factor, depending on the computational domain and the norm of the shape gradient, has to be introduced to the descent step size to ensure that geometry changes lead to a feasible lens shape, i.e. a shape which is still contained within the domain Ω and does not intersect itself. Its order of magnitude is taken to be 10 −3 / ∇J(Ω s l ) . We take the lower bound for the descent step to be TOL step = 10 −8 . Figure 6 depicts the final lens shape Ω final l together with the initial shape Ω 0 l and the desired goal shape Ω l . In Figure 7, the evolution of the relative change of the cost over the course of optimization is plotted as well as the relative change of the norm of the shape gradient. The cost, starting from J(Ω 0 l ) = 2.2762 · 10 5 ,  corresponding to a decrease of around 98.8 %. Similarly, the norm of the shape gradient was reduced by 99.99 % of its initial value ∇J(Ω 0 l ) = 2.1886 · 10 7 . We observe that the cost function and the geometry change rapidly in the beginning and then very little closer to the optimal shape. Such behavior is common for gradient flows (see, for instance, [12]). Furthermore, it can be observed that the cost values don't change significantly after the first 10 steps, hence a termination criterion that would stop the algorithm sooner can be employed.
As a way of measuring errors in the lens shape, we introduce the mean L 2 shape error (cf. [55]): 7.3.5. Both lens boundaries moving. For our second experiment, we allow both the upper and lower lens boundary to move. We start with an initial lens geometry Ω 0 l = Ω l,both perturbed , where both lens boundaries are perturbed with respect to the goal shape. Unlike the previous case, now the initial shape has to expand in order the match the target lens. The exact domain measurements are again given in Table 2.
We will first, as before, perform tests where linear NURBS are used to obtain the synthetic data and then to solve the problem. Secondly, we will employ quadratic NURBS to obtain the artificial pressure data and then solve the shape optimization problem.
In Figures 9 -12, we show the same quantities of interest as in the first example.  They exhibit a similar behavior. The cost functional, its gradient as well as the shape errors decay rapidly within the first few steps, while almost no changes happen once the shape is close to its goal. When using the linear NURBS, the cost functional decreased by 97.85 % of its initial value, going from J(Ω 0 l ) = 1.532·10 5 to J(Ω final l ) = 3.298·10 3 . The values are very similar in the quadratic NURBS case (see Figure 10).   Figure 11. Norm of the shape gradient over the course of optimization. left: Linear NURBS. right: Quadratic NURBS. Local, but not global minimum. The third experiment with artificial data demonstrates that the final shape obtained with our algorithm does not necessarily have to be a global minimum. We start with the initial lens with a straight upper boundary Ω 0 l = Ω l,upper straight , and take the target lens to be Ω l = Ω l,both down , see Table 2. Ideally, both lens boundaries should move down during optimization in order to reach the target shape. However, it seems that our initial lens is already close to a local minimum. The final shape is similar to the initial one and represents a local, but not a global minimum. However, the decrease in cost is still significant, going from J(Ω 0 l ) = 1.388630 · 10 4 to J(Ω final l ) = 3.200507 · 10 3 , which amounts to a reduction of about 76.95 % of the initial value. The experiment was conducted with the linear NURBS. The lens shapes and the relative change of the cost function are given in Figure 13.  Figure 13. left: Initial, final and goal shape; the final shape is a local, but not a global minimum. right: Relative cost decay.

Test 2: Unknown solution.
We now proceed to the case where the goal shape is not known. We assume the desired focal point to be on the axis of symmetry with coordinates (0, y fp ) T and associate with it the desired pressure distribution u d .
Finding an adequate analytical expression for the nonlinear target pressure distribution is challenging. Ideally, u d should assign high pressure values to the region around the focal point and low values elsewhere. We will employ a stationary Gaussian centered around (0, y fp ) T with a scaling factor in the order of magnitude of the maximal pressure amplitude within the computational domain: Here Σ denotes the covariance matrix with parameters σ x = 0.02, σ y = 0.004. The y-coordinate of the focal point is given by y fp = 0.105 and the amplitude of the Gaussian by A = 6 · 10 7 . Having in mind a realistic shape for the focusing lens, we will include an additional stopping criterion in the algorithm, that will stop the optimization process if an update would lead to a shape being thinner than some tolerance. A constant thickness tolerance over the width of the lens would already fail in the first gradient step, since the thickness of the lenses is zero at the tip. To determine a lower thickness limit that varies over the width of the lens, we set up a reference lens, which represents the thinnest lens design that is possible to manufacture (see Figure 14). If the newly computed lower and upper lens boundary do not fulfill the minimal thickness condition at all control points, we reduce the descent step size until either the minimal thickness limit is matched again for all control points or until the step size undergoes a lower bound, in which case the algorithm terminates.
We will use a reference lens which has two arcs as upper and lower boundaries. The arcs have horizontal tangents at the point on the axis of symmetry and a width on the axis of 0.01, see Figure 14. The thickness constraint is given by The result of the optimization process with the given method and parameters is shown in Figure 15. The initial shape is Ω l,Gauss , specified in Table 2, the domain D over which the pressure is tracked was chosen to be D = [0, 0.015] × [0.1, 0.11]. Figure 15 shows the results of the optimization process with linear and quadratic NURBS, where we took the lower bound for the descent step size to be TOL step = 10 −3 . For the employed mesh sizes, we refer again to Table 3, where we note that the numbers of degrees of freedom are of the same order of magnitude. In both linear and quadratic NURBS cases, we notice a decrease of the cost functional by 0.89 % − 1.15 % of its initial value. The small decrease is likely caused by using a stationary Gaussian as the desired pressure distribution, which our moving pressure wave can hardly fit to. Future research should include other objective functions which would not require an analytically given target pressure distribution, for instance, maximizing acoustic pressure in a certain area. The resulting final shapes, although similar, still differ. We would like to compare which shape performs better. We observe that the quadratic NURBS give a better relative decay of the cost functional (see Figure 15). In order to test if the quadratic NURBS found a better shape, we will bring the two shapes into the same setting. First, we interpolate the final quadratic NURBS shape into the linear NURBS space and compute the cost values. Secondly, we reverse the process and represent the "linear" final shape with quadratic NURBS and compute the cost in this quadratic NURBS space. The resulting values of the cost functionals are listed in Table 5, confirming that the quadratic NURBS perform better in both cases by around 0.235 %. In Figure 16, snapshots of the pressure wave propagation with the final quadratic NURBS shape are shown. Finally, we take the optimal "quadratic" lens (interpolated into the linear NURBS space) as the initial geometry and restart the optimization algorithm in the linear NURBS setting. The results of this post-processing optimization, given in Figure  17, confirm that this lens is already a local minimum. The algorithm performs 3 steps, resulting in a cost reduction of only 0.0729 % of its starting value, with no visible changes to the shape of the lens. The steps that the algorithm needed to take are likely due to the approximation error made by reducing the degree of the NURBS curve.

Conclusion and outlook
We have developed a gradient-based algorithm to find a locally optimal shape of an acoustic lens in nonlinear ultrasound applications. The shape derivative of the pressure-tracking cost functional was first rigorously derived by employing a variational approach. The algorithm was then numerically realized within the isogeometric concept, where geometries can be exactly represented in every gradient step. Due to the nonlinear nature of the problem and high source frequencies, a very fine mesh and time discretization had to be employed, increasing the computational complexity. Numerical experiments show that the quadratic NURBS perform better with respect to the final lens shape than the linear ones, reinforcing the potential of isogeometric shape optimization.
There are many possible directions for further research. Combining linear and higher-order NURBS in optimization in a clever way could result in a more efficient algorithm. Future research should also include considerations of different objective functionals, where an analytical expression for the target pressure distribution would not be needed. Exploring other acoustic models for nonlinear ultrasound propagation is of interest as well. Finally, from the point of view of applications, it is important to numerically tackle a 3D setting with even higher source frequencies, which would require developing a high-performance computational framework.