University of Southern California Dissertations and Theses

Material, shape, and topology estimation in coupled structural acoustics: sensitivity-based and data-driven approaches

(USC Thesis Other)

Material, shape, and topology estimation in coupled structural acoustics: sensitivity-based and data-driven approaches

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content MATERIAL, SHAPE, AND TOPOLOGY ESTIMATION IN COUPLED
STRUCTURAL ACOUSTICS:
SENSITIVITY-BASED AND DATA-DRIVEN APPROACHES
by
Harisankar Ramaswamy
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(MECHANICAL ENGINEERING)
March 2025
Copyright 2024 Harisankar Ramaswamy

To the creator of the universe
ii

Acknowledgments
I would like to begin by offering my salutations to the Omnipotent, the Omniscient, and the
Omnipresent for Thy divine grace—I bow before Thee in humble reverence.
First and foremost, I want to express my heartfelt gratitude to my advisor, Prof. Assad
Oberai, whose infinite patience and guidance led me through the challenges of my PhD
journey. I am deeply grateful for his willingness to accept me as a doctoral student and for
ensuring that I received continual financial support. His kind and optimistic nature has been
an invaluable source of inspiration and motivation, especially during moments when I felt
overwhelmed. I believe that our time together has inspired in me a more positive outlook
on both my abilities and my work.
I also thank my doctoral committee members, Prof. Paul Plucinsky and Prof. Bora
Gencturk, for their generous commitment of time, as well as their valuable expertise and
thoughtful critiques. Their insights and guidance have been instrumental in shaping and
refining my dissertation. A special thanks to Prof. Chia Wei Hsu for his time and encouraging
advice as my qualifying examiner. He had graciously agreed to serve on my dissertation
committee but sadly departed this world before my defense.
I wish to thank Agnimitra, who helped me with my research, being a postdoc in the
Computation and Data Driven Discovery (CD3) group. I appreciate his generosity in sharing
his code and assisting with debugging challenges. I am thankful to Ken Foo for providing
the painstakingly curated datasets that were essential to my research.
I gratefully recognize the financial support from the U.S. Naval Research Laboratory and
the National Institutes of Health, whose grants funded various parts of my Ph.D. Additionally, I acknowledge the Center for Advanced Research Computing (CARC) at the University
of Southern California for providing the computing resources that contributed to the research
results reported in this thesis.
iii

My sincere thanks go to Justin, my lab mate during my initial years, for his help with
Linux, Vim, and high-performance computing. His companionship made my transition
into the research group much smoother. I also extend my appreciation to all my other
lab group members for their support and intellectually stimulating discussions—especially
Dhruv, Deep, Orazio, Javier, Iman, and Bozhou (our unofficial group member)—all of whom
collaborated with me at some stage of my time at USC.
I also wish to acknowledge everyone else who was part of my journey, including those I
may have unintentionally left out, for their invaluable time and interactions.
Finally, my deepest gratitude goes to my parents for their unwavering faith, motivation,
and support—without which I would not have been able to complete, or even seriously
consider pursuing, a Ph.D.
In this thesis, I used AI writing tools for grammatical corrections and style refinements of certain sentences
for the sake of clarity. All ideas, data, arguments, and structure remain my own.
iv

Table of Contents
Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii
Chapter 1:
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Design optimization problems . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Inverse problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 Aim and scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.5.1 Material optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.5.2 Shape optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.5.3 Topology estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.6 Contributions of this work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.6.1 Material estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.6.2 Shape estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.6.3 Topology estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Chapter 2:
Material Parameter Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.1.1 Non-dimensionalization . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.1.2 Weak form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.1.3 Sensitivity analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.1.4 Choice of functionals . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.5 Choice of regularization . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.2 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.2.1 Sound minimization . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.2 Vibration isolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
v

Chapter 3:
Shape Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.1 Problem definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.1.1 Weak form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.1.2 Shape sensitivity analysis . . . . . . . . . . . . . . . . . . . . . . . . 50
3.2 Numerical approximation of gradients . . . . . . . . . . . . . . . . . . . . . . 54
3.3 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.3.1 Benchmark problem: Hole in a plate . . . . . . . . . . . . . . . . . . 57
3.3.2 Sound abatement problem . . . . . . . . . . . . . . . . . . . . . . . . 59
Chapter 4:
Topology Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.1 Problem definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.2 Conditional generative problem . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.3 Score based diffusion models . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.3.1 Langevin dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.3.2 Denoising score matching . . . . . . . . . . . . . . . . . . . . . . . . 71
4.4 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.4.1 Workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.4.2 Model architecture and hyper-parameter . . . . . . . . . . . . . . . . 78
4.5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.5.1 Benchmark problems in elastography . . . . . . . . . . . . . . . . . . 82
4.5.2 Inverse acoustic scattering . . . . . . . . . . . . . . . . . . . . . . . . 98
Chapter 5:
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.2 Future directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
vi

List of Tables
2.1 Material properties for the sound minimization problem. . . . . . . . . . . . 29
2.2 Parameters for the vibration isolation problem. . . . . . . . . . . . . . . . . 39
4.1 Training and sampling hyper-parameters for the cSDMs . . . . . . . . . . . . 82
4.2 Random variables comprising the parametric prior for the synthetic quasistatic elastography problem. Note, U(a, b) denotes an uniform random variable between a and b . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.3 Root mean squared error (RMSE) between the pixel-wise posterior statistics
estimated using cSDM and MCS, for the test phantoms in Figure 4.7 . . . . 87
4.4 Root mean squared error (RMSE) between the pixel-wise posterior statistics
estimated using cSDM and MCS, averaged over ten test phantoms for the
synthetic quasi-static elastography example . . . . . . . . . . . . . . . . . . . 89
4.5 Root mean squared error (RMSE) between the pixel-wise posterior statistics
estimated using cSDM and MCS for misspecified values of ση . . . . . . . . . 90
4.6 Random variables comprising the parametric prior for the tumor spheroid
application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.7 Random variables comprising the parametric prior for the acoustic scattering 102
vii

List of Figures
1.1 Types of design optimization problems in structural acoustics. (a) Example
problem. (b) Material Optimization. (c) Shape Optimization. (d) Topology
Optimization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1 Problem definition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2 Schematic of the sound minimization problem. RoI is a circle with 0.2 m
radius shifted 0.4 m from the center. The width of control region is 0.15 m. . 28
2.3 Finite element mesh used to solve the sound minimization problem. . . . . . 29
2.4 Taylor test for the sound minimization problem with ψ = 1−0.1( x+1.25
2.5
). The
dashed lines represent slope = 2. . . . . . . . . . . . . . . . . . . . . . . . . 30
2.5 Variation of the QoI as a function of the regularization parameter for the
sound minimization problem with smooth material properties. . . . . . . . . 31
2.6 Variation of the objective function as a function of iteration number for the
sound minimization problem (CV: smoothly graded material; Bimat: bimaterial distribution). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.7 Distribution of displacement and pressure for the sound minimization problem; QoI: 2.20e-1 (initial), 4.47e-4 (final). . . . . . . . . . . . . . . . . . . . . 32
2.8 Optimal Young’s modulus distribution for the sound minimization problem
with smooth material properties. . . . . . . . . . . . . . . . . . . . . . . . . 33
2.9 Optimal distribution for three different starting guesses. . . . . . . . . . . . . 34
2.10 Displacement and pressure fields of the optimal solutions corresponding to the
three initial guesses (a), (b), and (c); the QoI values are 4.47e-4, 2.10e-4, and
1.35e-4 respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
viii

2.11 Variation of the QoI as a function of the regularization parameters in the case
of binary distribution after 150 iterations. . . . . . . . . . . . . . . . . . . . 36
2.12 Distribution of displacement and pressure for the sound minimization problem
at optimal material properties (bi-material case); QoI = 9.91e-5. . . . . . . . 37
2.13 Optimal Young’s modulus distribution for the sound minimization problem
with bi-material properties. . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.14 Setup for the vibration isolation problem. The location of source: 1.5 m from
center of shell; the thickness of RoI and the control region are 0.1 m and 0.3
m; the width of PML is 0.7 m. . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.15 Finite element mesh used to solve the vibration isolation problem. . . . . . . 40
2.16 Taylor test for the vibration isolation problem with ψ = 1 − 0.1( x+0.4
0.8
). . . . 40
2.17 Variation of the objective function as a function of iteration number for the vibration isolation problem (CV: smoothly graded material; Bimat: bi-material
distribution). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.18 Initial (pre-optimization) distribution of displacement and pressure fields for
the vibration isolation problem with smooth material properties; QoI = 3.48e-2. 42
2.19 Final (post-optimization) distribution of displacement and pressure fields for
the vibration isolation problem with smooth material properties; QoI = 6.99e-4. 43
2.20 Optimal distribution of Young’s modulus within the structure for the vibration
isolation problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.21 Final (post-optimization) distribution of displacement and pressure fields for
the vibration isolation problem with bi-material property distribution; QoI =
5.12e-4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.22 Optimal distribution of Young’s modulus within the structure for the vibration
isolation problem with a bi-material distribution. . . . . . . . . . . . . . . . 45
3.1 Illustration of the domain of optimization problem. The gray region illustrates
the solid, blue region illustrates the fluid domain, and the violet region is the
region of interest. Γϵ
I
is the small incremental deformation of ΓI . . . . . . . 47
3.2 Construction of auxiliary domain around the closed interface. . . . . . . . . . 55
ix

3.3 Shape optimization on an elastic plate. The auxiliary region is highlighted in
yellow on the left figure (initial configuration) and in green on the right figure
(optimized configuration). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.4 Schematic of the sound minimization problem. The RoI is a circle with 0.2 m
radius shifted 0.4 m from the center. The width of the auxiliary region is 0.2 m. 60
3.5 Computation of the extended normal field and subsequently the curvature field. 60
3.6 Variation of the objective function as a function of iteration. . . . . . . . . . 61
3.7 Optimal shape of the interface along with deformed mesh of the solid region. 61
3.8 Distribution of displacement and pressure for the sound minimization problem; QoI: 2.44e-1 (initial), 7.19e-4 (final). . . . . . . . . . . . . . . . . . . . . 62
4.1 Inverse acoustic scattering problem. . . . . . . . . . . . . . . . . . . . . . . . 63
4.2 Illustration of the forward and reverse SDE processes in diffusion models. . . 67
4.3 An example of 1-dimensional forward diffusion process which maps a complex multi-modal distribution to a Gaussian. Each colored cross represents a
sample, with identical colors indicating how a single point transforms under
increasing noise through the SDE. . . . . . . . . . . . . . . . . . . . . . . . . 67
4.4 A schematic illustrating the forward diffusion process in one dimension across
various time (or noise) levels σ(t). The original multi-modal distribution
PX(x) at t = 0 gradually transitions into smoother distributions as noise
increases. Eventually, at t = T, the distribution becomes close to a simple
Gaussian N (µX, 5
2
). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.5 Five typical realizations of X and Y sampled from the joint distribution for
the synthetic quasi-static elastography example and belong to the training
dataset. The first row shows the shear modulus field, while the second row
shows the corresponding noisy vertical displacement field used as measurements. In this figure, the standard deviation of the Gaussian noise ση = 0.025.
Moreover, realizations of X are Y are min-max normalized to [0, 1] scale . . 84
4.6 Visualization of the training process using a single realization sampled using
ALD. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
x

4.7 Posterior statistics – pixel-wise posterior mean and standard deviation – estimated using the conditional diffusion models (cSDMs) and Monte Carlo
simulation (MCS) for the synthetic quasi-static elastography example on two
test samples. In subfigures (a) and (b), the three rows show the results corresponding to three different levels of measurement noise ση but for the same
ground truth. All values are normalized to [0, 1] scale . . . . . . . . . . . . . 88
4.8 Posterior statistics – pixel-wise posterior mean and standard deviation – for
test sample 1 estimated using cSDMs trained assuming measurement noise
with standard deviation ση = 0.05. The reference MC statistics are also
estimated with the misspecified measurement noise model. All values are
normalized to [0, 1] scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.9 Schematic showing the experimental setup in the tumor spheroid application 91
4.10 Workflow for generating the training data in the tumor spheroid application 94
4.11 Five typical realizations of X and Y sampled from the joint distribution that
form the training dataset for the tumor spheroid experiment. In the first row,
we have the log-normalized Young’s module fields, and in the second row, the
corresponding noisy measurements are displayed. All values are normalized
to [0, 1] scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
4.12 Posterior statistics estimated for selected samples using the cSDM for the
tumor spheroid problem. The field of view in each image is 256×256 µm2
.
Note that all images are presented at their true scales. GT: Ground truth. . 97
4.13 Schematic of the data generation process. The top-middle diagram illustrates
the analysis domain with its surrounding PML region. At top-left is an example realization of X, and at top-right is the corresponding scattered pressure
field (magnitude) obtained via FEA. The bottom plots display the boundary
measurements, color-coded according to the three different measurement cases. 99
4.14 Five typical realizations of X and Y sampled from the joint distribution
for the acoustic scattering example. The top row shows the binary scatterer
topology field X, and the subsequent rows present the noisy back-scattered
pressure signals Y . Here, Y0, Y1 correspond to the real and imaginary components of the single-frequency case, and Y2, Y3 represent the multi-frequency
case. All data are min–max normalized to [0, 1] scale. . . . . . . . . . . . . . 103
4.15 Posterior statistics estimated for case 1, case 2, and case 3 using measurements
corresponding to the ground truth signal (top row), for four representative test
samples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
xi

5.1 Comparison between the the results obtained using the state-of-the-art method
and the cSDM approach for the tumor spheroid problem in Section 4.5.1.2.
Note, in the existing method, the inferred quantity (Young’s modulus) is directly estimated by computing the strain field; hence uncertainty measure is
unavailable. Algebraic results from Foo et al. . . . . . . . . . . . . . . . . . . 111
xii

Abstract
This thesis explores novel methodologies in the estimation of material properties, shape,
and topology of the structure for coupled structural acoustics systems by leveraging both
traditional optimization techniques and innovative conditional inference approaches. The
motivation for this work stems from the increasing demand for robust and cost-effective
techniques for design applications, where accurate modeling of the structural geometry along
with interface dynamics is critical. The techniques presented in work can be employed in
a variety of applications such as vibrations of submerged structures and fluid filled containers, acoustic radiation scattering, functional acoustic metamaterials design and acoustic
medical sensors. They can also be extended to other wave-propagation domains including
electromagnetic propagation and scattering. This thesis can be broadly divided into two
parts: Part 1 utilizes conventional deterministic approaches, including sensitivity analysis,
to solve material and shape optimization problems in structural acoustics; Part 2 employs
probabilistic methods and state-of-the-art generative machine learning models, specifically
diffusion models, to solve inverse problems in elastography and acoustic scattering.
The first part of this dissertation focuses on semi-analytical deterministic approaches for
solving material and shape optimization problems in coupled fluid-structure acoustics under
time-harmonic loads. The fluid and structural components are modeled using the Helmholtz
and Navier-Cauchy equations, respectively, with interface conditions explicitly enforced. We
introduce a computationally efficient and generalized formulation for computing adjointbased sensitivities, which drive gradient-based minimization algorithms for both material
and shape optimization. This approach is robust, computationally efficient, and seamlessly
integrates into conventional finite element method (FEM) frameworks.
In material optimization, we explore continuous and bi-material property distributions,
leveraging recent advancements in manufacturing functionally graded materials. These adxiii

vancements have expanded the design space, allowing for more sophisticated and feasible
material configurations.
For shape optimization, we explore parameterization-free methods, known as node-based
sensitivity approaches. These approaches encompass both level-set and artificial densitybased methods that utilize conventional meshing techniques. However, level-set methods
often face numerical stability and initialization challenges, particularly in 3D applications,
while density-based methods tend to lack the sharp interfaces required for applications where
precise interface definitions are critical, such as in airfoil design. Although Isogeometric
Analysis (IGA) is often the preferred method for shape optimization due to its integration of
Computer-Aided Design (CAD) and Finite Element Analysis (FEA), the industrial adoption
of such advanced techniques remains in its nascent stages. To address these limitations, we
propose a novel method that integrates smoothly within conventional frameworks. A key
challenge in node-based sensitivity approaches, particularly for coupled physics problems, is
computing the shape derivative, which involves complex terms such as surface curvature and
the gradient of interface-specific quantities. Additionally, mesh updates in each iteration
often lead to distortions and artifacts, a recurring issue in node-based methods. Although
numerous studies have addressed mesh updating challenges, effective computation of sensitivities for interface problems remain underexplored in the literature. To overcome this
challenge, we introduce a novel technique that constructs a smooth extension of the surface normal vector field across an auxiliary domain intersecting the interface, enabling us
to compute curvature effectively without additional parameterization. The efficacy of this
approach is demonstrated through its application in an interior noise minimization problem
in structural acoustics.
The second part of this dissertation transitions to a data-driven probabilistic approach
to address inverse problems, with a specific focus on inferring the topology of acoustic scatterers, a critical component in underwater sensing applications such as SONAR. By utilizing
simulation-based inference, we sample from posterior distributions linked to an unknown
xiv

forward model, employing state-of-the-art generative machine learning models, notably conditional diffusion models, to learn the structure of these distributions. Although generative
models are proficient at handling geometrically complex priors, they encounter challenges due
to the high dimensionality and nonlinearities inherent in the forward model. These complexities can restrict the comprehensive exploration of the distribution space, potentially resulting
in predictions confined to regions around local modes, a phenomenon commonly referred to
as mode collapse. To overcome these limitations, we introduce a framework employing conditional score-based diffusion models for posterior inference. These models excel at approximating the score function of a posterior distribution using samples from an extensive dataset
of pairwise data generated via finite element forward solvers. This setup demonstrates our
method’s capability to manage black-box forward models and intricate measurement noise
effectively. Uncertainty quantification is accomplished by computing posterior statistics from
numerous realizations sampled through Langevin dynamics driven by the learned score function. The efficacy of this approach is first validated through applications in elastography,
employing both synthetic and actual experimental data. We demonstrate the model’s ability
at handling diverse measurement modalities, sophisticated geometries, and nonlinear noise
models, by comparing their performance to traditional Monte Carlo and algebraic methods.
We subsequently apply this methodology to solve the inverse scattering problem for three
distinct sensor configurations and derive intuitive insights from them.
xv

Chapter 1
Introduction
If you want to find the secrets of the
universe, think in terms of energy,
frequency, and vibration.
Nikola Tesla
1.1 Motivation
We define structural acoustics as the study of propagation of mechanical waves in coupled
fluid-structure systems, typically involving elastic solids interacting with acoustic fluids. The
dynamics of these systems depends on the structural properties such as material composition
within the structure, the geometry of the structure, and the shape of the interface that
separates the two domains.
Design methodologies for structural acoustic systems involve accurately estimating these
properties. Therefore, techniques for solving such estimation problems have wide ranging
practical relevance. For instance, in acoustic devices such as hearing aids, earphones, and microphones, acoustic performance is particularly sensitive to mechanical design modifications,
as strong mechanical-acoustic coupling in compact enclosures can result in undesirable feedback noise. At the macro-scale level, the ability to design structures capable of minimizing
1

interior noise and reducing vibration amplitudes has significant implications in automotive,
aerospace, civil, and naval engineering. Notable examples include designing quiet cabin interiors [1, 2, 3] and developing structural foundations that improve resilience against seismic
activity [4, 5, 6, 7]. Recent advances in acoustic metamaterials have opened up additional
avenues of research, including seismic metamaterials [8], acoustic cloaks [9], acoustic superlenses [10, 11], and enhanced acoustic sensors [12]. Designing these systems typically involve
optimizing performance metrics such as sound pressure levels (SPL), the strength of backscattered signals, and displacement or vibration amplitudes. Furthermore, these methods
can effectively be used to solve inverse problems such as in ultrasound elastography [13] for
biomedical imaging, and in sensing applications like sonar [14].
Estimation problems in structural acoustics can be broadly categorized into three main
types: material, shape, and topology. To understand the difference between these three
types of property estimation categories, let us look at an example shown in Figure Figure 1.1. Consider a linear elastic continuum solid (grey) which is immersed in an acoustic
fluid. Plane waves are incident from the left and are transmitted through the solid structure.
The task here is to design the structure such that it efficiently concentrates the transmitted
acoustic waves at a designated target location. In practical applications, this target could
correspond to the placement of an acoustic sensor, positioned to maximize signal receptivity. This design problem can be formulated as an optimization problem with an objective
to maximize the acoustic signal amplitude at the target location. Alternatively, the same
setup can be employed to formulate an inverse problem, where the goal is to estimate unknown structural properties that closely replicate a desired acoustic field distribution after
incorporating relevant prior knowledge. Solutions to these design or inverse problems can be
obtained by manipulating the structure’s material property, interface shape, or its topology,
while ensuring compliance with constraints such as staying within a predefined bounding
box (depicted in black).
Next, we introduce three categories of property estimation from a functional optimization
2

Figure 1.1: Types of design optimization problems in structural acoustics. (a) Example
problem. (b) Material Optimization. (c) Shape Optimization. (d) Topology Optimization.
perspective; note that these categories also extend to inference-based approaches. The only
difference is that, in a functional optimization setting, the property to be estimated becomes
the design variable, whereas, in an inference-based setting, it becomes the inferred variable.
Material optimization: In the case of material optimization, we vary the material property (such as elastic modulus, density, etc.) within the structure to achieve the desired
behavior of the system. Further, we could enforce constraints on this distribution that require it to be smooth or piecewise constant (an assembly of multiple components) depending
on the designer’s requirements. Note that the shape of the structure remains unchanged
after optimization, as seen in the example Figure 1.1b). One implicit advantage of optimizing the material is that it does not involve altering the fluid-structure interface. As a
result, the shape of the interface can be determined by other important considerations, such
as the desire to optimize aerodynamic quantities like the lift to drag ratio in an aircraft.
Furthermore, the ability to combine the optimization of the material parameters along with
shape (or topology) optimization opens up broader avenues to modify the specific acoustic impedance of coupled elasto-acoustic systems to achieve design/performance goals as it
relates to its structural-acoustic response.
3

Shape optimization: The idea behind shape optimization (Figure 1.1c)) comes from
the fact that the shape of the interface plays a critical role in determining the coupling
behavior between the two media; thus, the behavior of the system can be modified by
altering the shape. Compared to material optimization, the complexity in shape optimization
arises from the fact that coupling interface is changing throughout the optimization process.
Like any multi-physics problem, there are certain coupling conditions that must be satisfied
at the interface between the two media. In structural acoustics, these conditions are the
continuity of displacements and traction. Therefore, the question of how to represent the
shape accurately while it deforms becomes important in shape optimization.
Topology optimization: Topology optimization is a more generalized version of shape
optimization, where the final design is allowed to have disconnected regions. In addition,
these regions may have fluid filled cavities present inside them as shown in Figure 1.1d).
Typically, topology optimization problem is solved with some constraint on the maximum
amount of solid material that can be present. Further, since the solution space is vast
allowing all kinds of complex geometries, this often leads to designs which are impractical to
manufacture. Therefore, the challenge of imposing the constraints on the feasibility of the
design space becomes crucial.
Remark: We use the term “design optimization” to describe the nature of the problem,
while “functional optimization” denotes a particular approach applied to solve either design
optimization problems or inverse problems.
1.2 Design optimization problems
There are various approaches one could take to solve these optimization problems. The most
standard approach is to compute the local gradients of the objective function with respect to
the design variables and use this information to drive the optimization algorithm iteratively.
4

Some of the widely used algorithms that employ this approach include Broyden Fletcher
Goldfarb Shanno (BFGS), Sequential Quadratic Programming (SQP), and the Method of
Moving Asymptotes (MMA). The convergence of these algorithms is dependent on the accuracy of the supplied gradients. Thus, the development of an efficient method to compute
the sensitivities of the objective function with respect to the design variable becomes crucial
for any optimization problem.
The finite difference method offers the simplest way to compute the gradients, however
the cost of computation of the gradient at each iteration is directly proportional to number
of design variables times the cost of one forward solve (i.e. solving the governing equations
for finding the state variables of the system). Therefore, this approach becomes impractical
for large DOFs system (can be larger than a million for fidelity simulations). The effective
way to compute gradients in such scenarios is the adjoint method, which just requires a
single forward solve and an adjoint solve for gradient computation. The formulation of the
adjoint equation however may require some “mathematical effort”, depending on the complexity of the constraints. Both the forward and the adjoint equations are then discretized
and solved numerically using the Finite Element Method (FEM) or the Boundary Element
Method (BEM). The BEM offers a way to handle infinite unbounded domain problems with
discretization applied only to the boundaries.
Then there are other approaches, which have been employed for structural optimization, that do not require the computation gradients. Heuristic algorithms, such as genetic
algorithm [15], particle swarm optimization [16], and simulated annealing [17], fall in this
category. They can search for solutions that are arbitrarily close to the global optimum,
however, they require exceedingly large number of objective function evaluations (forward
solves), depending on the size of the search space. This limitation makes them unattractive
for the purpose of optimizing continuum structures. Yet, another potential approach which
has recently gained attention, is the use of Deep Neural Networks (DNNs) for learning the
map from user inputs to the optimal solution. In the context of topology optimization,
5

several researchers have attempted to use these data-driven tools in the supervised learning
paradigm to accelerate the optimization process [18, 19, 20, 21]. The key feature of this
approach is that they (after training) directly predict the optimal topology without the need
for an iterative strategy (optimizer), when design parameters such as loading conditions [20]
or partially converged topologies [21] are feed as input. However, they all require an extensive dataset containing optimized solutions for various input conditions to begin with and
a generative neural networks architecture is used as black-box interpolation tool to find the
optimal solution.
1.3 Inverse problems
The mathematical formulation of inverse problems involves defining a forward model F :
X 7→ Y that describes how the input parameters of a system x ∈ X relate to its observable
outputs y ∈ Y . The challenge is then to invert this model, often using noisy estimates of the
output quantities y¯, to recover the input parameter estimate x
∗
. A typical inverse problem
in mechanics involve inferring the spatially varying constitutive material parameters of a
system, given measurements of its response to external mechanical loads. However, X could
also be tailored to represent domain geometries, external mechanical loads distributions,
initial conditions depending on the problem. Here, we refer to X as the inferred field and
Y as the measured field.
The deterministic approach to solving such problems typically involves formulating them
as optimization problem, where the objective is to minimize a loss function. This loss function
is defined as the distance between the predicted and observed outputs. A regularization term
is usually added to the loss function to mitigate the effects of noise and prevent over-fitting.
x
∗ = argmin
x∈X
∥F(x) − y¯∥
2
2 + R(x) (1.1)
However, using minimization approach to solve inverse problem has two major limitations:
6

• Non-uniqueness Inverse problems generally lack a unique solution, in other words they
are ill-posed. However, deterministic methods yield a single solution, likely converging
to a local minimum that is close to the initial guess. Hence, this approach limits the
exploration of the broader space of feasible solutions and may overlook other potentially
interesting solutions.
• Uncertainty quantification The deterministic approach does not inherently provide
a means to quantify uncertainty in the predictions. Without an estimate of the uncertainty, it is difficult to assess the reliability of the result, which is important in fields
like medical imaging where diagnostic decisions are based on the inferred parameters.
The probabilistic approach, on the other hand, offers a way to handle both the limitations
by estimating the probability distribution of X given Y , as opposed to a single realization
of X. The well-known Bayes’ theorem tell us that the conditional density PX|Y (x|y¯) is proportional to the prior distribution of X, denoted by PX(x), and the likelihood of obtaining
the measurement y¯ when x is assumed as the ground truth, which is denoted by PY |X(y¯|x).
PX|Y (x|y¯) = PY |X(y¯|x) · PX(x). (1.2)
The relationship between the measurements and the inferred field is captured by the likelihood term which has the forward map implicitly baked into it along with measurement noise
model. Let η be the random variable associated with the noise model, then
PY |X(y¯|x) = Pη|Y (η|F(x)). (1.3)
For the special case where additive Gaussian noise is used, y¯ = η +y. However, we can have
more complex, non-linear noise models too, as we will see in one of our examples problems.
Note that we abuse notations hereon by using y in place of y¯ for the purpose of enhancing
readability.
7

In real-life applications, the prior information (distribution) typically comes in the form of
data-samples, providing insights into the probable values or geometries that X may assume.
Generally, for physics-based inverse problems, a closed form mathematical expression for
PX(x) may not be feasible, As a result, generative neural network models are widely used
for this purpose.
1.4 Aim and scope
The goal of this research is to develop scalable and robust methods for solving material,
shape, and topology estimation problems in structural acoustics. Most estimation approaches
found in literature are tailored towards topology optimization and there are very few that
consider material or shape optimization. Moreover, they are implemented for 2D problems
and extending them to 3D is either cumbersome or requires a disproportionate increase in
computational cost. In this work. we explore strategies which are agnostic in terms of
geometric complexities (no additional effort from user) while keeping the DOF minimal as
possible. Furthermore, to the best of our knowledge, no prior studies have applied probabilistic techniques in this context. In this work, we address both design optimization and
inverse problems; we demonstrate our deterministic optimization methods for material and
shape estimation, and a probabilistic inference method for topology estimation.
For solving material and shape optimization, we take a semi-analytical approach which
includes the computation of variational adjoint sensitivities in conjunction with the FEM.
The sensitivity formulations are derived for a general case so that they are applicable for both
interior and exterior problems (unbounded fluid domain) with any given objective function
or boundary conditions.
For solving topology inference, we adopt a data-driven probabilistic approach focused on
inferring the posterior distribution of the scatterer topology in inverse scattering problem.
This posterior distribution is conditioned on the observed measurement of the back-scattered
8

pressure field. The topology is estimated by sampling the posterior distribution with the
help of a deep neural network model.
We limit the scope of this work to problems involving time-harmonic loading conditions
where the acoustic performance of the system is characterized by a single dominating frequency. This assumption simplifies the physics from a computational standpoint, enabling
one to rewrite the general wave equation which is a hyperbolic partial differential equation
(PDE) as an elliptic PDE in the frequency domain.
Remark. The methods proposed in this work can be easily extended to other steadystate multi-physics problems.
1.5 Literature review
1.5.1 Material optimization
There is relatively little work on the topic of material parameter optimization (MPO) on
coupled fluid-structure systems, in comparison to shape and topology optimization. Prior
works have been limited to modifying the stiffness properties (like the cross-section) of ‘thin’
structures that were modeled using beam [22] and shell elements [23, 24, 25, 26, 27], in order
to minimize the radiated sound from the structure. Perhaps this is because a few decades
ago altering the geometry of a structure was a more viable option compared to altering
its material property distribution within the structure’s volume. However, over the past
couple of decades, advances in materials science and manufacturing have led to fabrication
techniques that can offer greater control of material property distributions [28, 29]. This
includes complex distributions comprised of two constituent materials, or components with
spatially graded material properties. Therefore, it is natural to explore whether these design
spaces lead to better design [30].
9

1.5.2 Shape optimization
The conventional approaches for shape optimization rely on the parametrization of the
boundary using a polynomial or splines [31, 32]. This enables one to compute the curvature for sensitivity calculations and at the same time provides a reliable way to deforming
the boundary in a smooth fashion, without incurring unwanted boundary distortions. For
instance, this approach has been employed in [33] to find the optimal target shape for the
inverse scattering problem with a rigid scatter.
The authors in [34] have performed shape optimization to design an acoustic lens for
an ultrasound transducer. They have considered the mixed displacement/pressure (u/p)
formulation [35] to avoid remeshing the structure as the interface boundary deforms. The
mixed u/p formulation allows one to model both the media using the same set of governing
equations, thereby circumventing the need to impose interface conditions explicitly, at the
expense of increasing the DOF of the system. The material parameters that characterize the
PDE then take discrete values to distinguish between the two domains. For example, the
elastic shear modulus is set to zero for the fluid and has a finite value in the solid region.
They use a parameterized level-set function (has value < 0 inside the solid domain, 0 on
the boundary, > 0 outside) composed with a smoothened indicator function for mapping
the spatial variable to a smooth approximation of a binarized material field (a smooth field
enables the computation of the sensitivities using the adjoint method). Finally, the shape
optimization problem is solved by optimizing the level-set function parameters.
Another class of approaches involve the use of Iso-Geometric Analysis (IGA) proposed
by [36]. The core idea of IGA is to integrate the numerical analysis with the Computer
Aided Design (CAD) model by representing and discretizing the geometry accurately using
Non-Uniform Rational B-Splines (NURBS). Further, NURBS are also used as basis function
for performing FEA and reducing the computational time for remeshing. The IGABEM
method (IGA combined with BEM) have been used by [37, 38] for shape optimization in
10

exterior problems.
1.5.3 Topology estimation
Initial attempts to solve the topology optimization problem in structural acoustics involved
treating the solid as a rigid domain [39, 40] or as a stiff fluid with a large bulk modulus
[41, 42, 43]. In subsequent studies, the structure was modeled as an elastic solid with
displacement degrees of freedom. In most of these problems, the fluid was also treated as a
“solid” with negligible shear modulus using mixed u/p formulation. The interface was then
defined implicitly using a level-set function [44, 45] or a density function based approach
[46, 35, 47]. The approaches based on density function rely on the Solid Isotropic Material
Penalization (SIMP) technique to compute the gradients. The SIMP method and its variants
are popular density interpolation schemes used in bi-material topology optimization that
allows for a smooth transition of density (normalized) values between 0 and 1 (for stiff and
soft material respectively) to facilitate the evaluation of sensitivities with respect to the
density field (design variable). It used in conjunction with a volume constraint to penalize
intermediate densities. This is achieved by raising the design variable ρ to a power p > 1
when interpolating for the stiffness, E = E0 + (E1 −E0)ρ
p
. In case of level-set based schemes
in [44], the evolution of the level-set function is determined by the topological derivative.
The topological derivative [48] computes the change in the objective function with respect
to small change in topology due to the inclusion of an infinitesimal hole.
In some cases the interface was modeled explicitly by using the level-set function in
conjunction with adaptive mesh refinement to accommodate the changing interface boundary
[49]. Further, the shape derivative is used to update structural boundary which is followed
by reinitialization of the level set function. This leads to topological changes because two
holes can merge, or an existing hole can get filled-up; but new holes will not be added in the
structure.
In [50], the authors used a bi-directional evolutionary structural optimization (BESO)
11

on a finely refined mesh to impose explicit interface constraints. The design variable defined
on a fixed Cartesian grid is only allowed to have binary values. This way, the need to
update the mesh after each iteration is avoided. The downside is that evolving interface will
have corrugations which could degrade the analysis accuracy, unless the mesh size on the
design domain is sufficiently small compared to the wavelength. Therefore, this significantly
increase the computational cost for analyzing 3D geometries, especially for high frequency
applications.
Yet another approach for solving SAO, namely the CutFEM, was proposed by [51]. Here,
the interface is represented through a level-set function and the finite element problem is
solved on fixed rectangular mesh domain. The elements in the mesh that are cut by the
zero level-set are triangulated in a manner such the boundary (interface) is approximated
by a piece-wise linear representation. These cut elements are treated separately during FE
analysis, using special integration rules (on those elements) without modifying the original stiffness matrix. This way, the need for remeshing the entire domain with conforming
elements along the interface is avoided.
Finally, deep learning-based approaches have also been explored for topology representation. For example, specially designed convolutional neural networks [52, 53] have been used
to solve inverse scattering problems; they produce a single topology field that minimizes the
mismatch with the measured scattered field that is provided as input to the network. These
networks are trained on pairwise data, all generated through simulations.
1.6 Contributions of this work
1.6.1 Material estimation
In our work we describe a systematic and efficient approach to material parameter optimization (MPO) of continuum structures in order to achieve a desired structural acoustic goal. In
contrast to earlier works that use thin elements along with the BEM, we model the structure
12

as a linear elastic solid with the Navier-Cauchy equation and the fluid using the Helmholtz
equation. The BEM based approaches are best suited for geometries having small surface
to volume ratio (like in the case of exterior problems), whereas our approach can be used to
solve interior problems as well.
Salient features:
• We partition the structure into two sub-domains: one where the material properties
are fixed and cannot be altered, and another where they may be altered to achieve the
desired goal. This reflects the constraint that there may be regions within the structure
where the material properties are determined by other considerations, like strength.
• We allow the quantity of interest to be defined in the structural or the fluid domain.
Thus, using the same formulation, we solve a vibration isolation problem where the
goal is to reduce vibrations in a subdomain of the structure, and a noise abatement
problem where the goal is to minimize the acoustic amplitude in the fluid.
• We consider two formulations: one with a prior distribution that promotes smoothly
graded material properties and another that promotes binary (two-component) material property distribution.
• When solving the optimization problem, we represent both the state variables and
the material properties (same as design variable) using a finite element basis. The
optimal properties are determined by evaluating the gradient of the objective function
with respect to the nodal values of the material parameter using the adjoint equations.
This gradient is then used to drive a gradient-based minimization algorithm.
1.6.2 Shape estimation
We propose a novel method for computing the shape derivative with respect to the coupled
interface in SAO using the standard FEM solver. The shape derivative typically involves
13

a curvature term which is hard to compute, solely by using a conventional meshed. Prior
works have relied on polynomial interface parameterizations and level-set functions to provide
closed form expressions for the curvature, to mitigate this bottle neck. However, the use of a
parameterized models can limit the search space to a set (of shapes) that is within a particular
parametric family. Further, it is inconvenient to represent complex 3D surfaces using these
techniques (especially level-set functions). Although IGA offers a reliable solution to handle
all these limitations, existing industrial frameworks (computational) and most commercial
solvers (FEA) are yet to adopt IGA due to its tight coupling requirement between the CAD
and FE solvers. We offer an alternate approach that relies purely on traditional FEM,
starting with an independent mesh (geometry) provided as input.
Salient features:
• We solve the SAO problem by explicitly considering interface constraints. Additionally,
we represent the deformable interface solely using the mesh.
• The expression for the shape derivative is derived for the ”pre-descretized” problem,
and therefore is independent of the mesh. The problem is described in a generalized
manner, such that our formulation is applicable for both interior and exterior problems.
• The evaluation of shape derivative requires the computation of curvature on the interface as well as gradients of quantities that are defined on the boundary. We tackle
this problem by obtaining a smooth extension of these quantities on an auxiliary domain adjacent to the interface. This method is then validated by comparing with the
curvature values of a geometry for which the curvature is known (analytically) a priori.
• We demonstrate the efficacy of our approach by solving an interior noise abatement
problem. The shape derivatives computed are used to drive the interface boundary to
achieve an optimal (local) solution. The optimizer chosen for this is a simple steepest
descent algorithm with a constant step size.
14

• In each optimization step, we simultaneously update the entire mesh along with the
deforming interface boundary, to prevent the emergence irregularities (in the mesh).
The movement of the mesh is guided by the solution to an elasticity-type boundary
value problem defined over the entire domain.
1.6.3 Topology estimation
The main challenge in estimating the topology of coupled-physics systems is the representation of the interface boundary. This is made harder by the fact the topology of the structure
is allowed to vary, and therefore cannot be mapped to simple primitives. Broadly speaking, the methods by which this problem has been previously addressed in literature, can
be classified into two categories: one approach uses a level-set function and the other uses
density-based method (artificial material model). Choosing either one would usually lead to
a trade-off between re-meshing and accurate representation of the interface constraints with
significant increase in computational cost. The CutFEM [51], which was recently proposed,
provides a promising solution by combining the best of both approaches. The effectiveness of this deterministic method has been demonstrated for 2D problems. In contrast, our
approach leverages probabilistic deep learning techniques to infer the topology of the structure by modeling posterior distributions with diffusion models. The novelty of this method
lies in the integration of noise and measurement conditioned convolutional network (U-Net)
with diffusion models to enable sampling from the posterior distribution. This framework is
used to solve inverse problems in physics, particularly the notoriously challenging acoustic
scattering problem.
In this work, we consider rigid scatterers. However, we note that it can easily be extended
to more complex and realistic representations of the structure, including cases where it is
treated as an elastic solid.
The salient features of our approach are:
• We apply our method to inverse scattering problems involving noisy far-field scattered
15

signals measured under three distinct sensor arrangements. For each arrangement, we
estimate the scatterer topology and evaluate the associated uncertainty.
• We employ Langevin dynamics to sample from the posterior distribution. These samples (of the inferred topology) are generated by gradually modifying random samples
from a simple tractable distribution. The modification process is achieved through
Markov Chain Monte Carlo (MCMC), which is guided by a gradient term that encodes the posterior.
• The gradient term that is used to drive the dynamics is learned using a deep neural
network (specifically a U-Net). This network is trained on synthetically generated
joint data that includes both the inferred (estimated) fields and the corresponding
measurements.
• We validate our methodology by solving benchmark problems in elastography and
compare our results with those obtained through Monte Carlo (MC) simulations and
real-life experimental results.
• We show that our approach is capable of handling complex measurement noise models
and different modalities of input data.
16

Chapter 2
Material Parameter Estimation
2.1 Problem formulation
Consider the domain (see Figure 2.1) which consists of both solid (linear elastic) and fluid
regions, represented by Ωs and Ωf respectively, with ΓI as the interface. Let u be the
displacement field in the solid and p be the acoustic pressure field in the fluid due to an
excitation, pext, or a traction vector field h. The solid domain contains a region Ωc ⊂ Ωs,
where we are allowed to modify the material parameters. Our objective is to find the
optimum material parameter distribution E(x)(Young’s modulus) in Ωc that minimizes the
Quantity of Interest (QoI) F(u, p) defined over the region Ωr. Note that Ωr ⊂ Ωs ∪ Ωf can
be defined over the fluid and solid regions. For convenience, we express our parameter field
in terms of a design variable ψ(x) as E(ψ) = f(ψ)Eo, where 0 ≤ ψ ≤ 1; the reference value
Eo is typically the modulus of the surrounding solid region.
The loss functional J for our problem is defined as
Minimize
ψ(x), x∈Ωc
J(ψ) = F(u(ψ), p(ψ)) + R(ψ) (2.1)
where R(ψ) is the regularization term. The regularization term imposes prior information
about the solution. We note that we have expressed the displacement and the pressure fields
as functions of ψ since for a given instance of ψ, and hence the Lame parameters µ and λ,
17

they can be determined by solving the structural acoustics problem defined next.
Figure 2.1: Problem definition.
The equations of structural acoustics in the time domain are given by,
ρs
∂
2u
∂t2
− ∇ · σ = fu in Ωs (2.2)
ρf
∂
2p
∂t2
− κ∆p = 0 in Ωf (2.3)
u = 0 on ΓD (2.4)
σ · nˆ = h on ΓN (2.5)
p = pext on ΓDf (2.6)
∂p
∂nˆ
= 0 on ΓN f (2.7)
∂
2u
∂t2
· nˆ =
∇p
ρf
· nˆ on ΓI (2.8)
σ · nˆ = −pnˆ on ΓI (2.9)
Since the material (solid) is assumed to be linearly elastic, its stress-strain relationship is
given by
σ = µ(∇u + ∇u
T
) + λ(∇ · u)I (2.10)
18

Assuming a fixed Poisson’s ratio, we write λ and µ, the first and the second Lam´e parameters,
in terms of the design variable ψ,
λ(ψ) = νE(ψ)
(1 + ν)(1 − ndν)
, µ(ψ) = E(ψ)
2(1 + ν)
, (2.11)
where nd = 1 in the case of two-dimensional plane stress and nd = 2 in three dimensions.
In order to solve the minimization problem, we need to first compute the sensitivity
(gradient vector) of J with respect to ψ; this requires solving the system of partial differential
Eqs. (2.2) to (2.10) for different variations in parameter field ψ. Since we are interested in
the response at a fixed frequency, we transform these equations into the frequency domain,
by expanding any function of space and time as v(x, t) = v(x)e
iωt, where ω is the circular
frequency. Consequently, the state variables u and p and the material parameters µ, λ, and
κ take on spatially dependent complex values. The latter account for any viscous dissipation
that may be present in the structure or the fluid.
2.1.1 Non-dimensionalization
We begin our analysis by considering suitable scaling terms Uo, Po, λ0, and µ0 for the state
variables and the material parameters, respectively. This leads us to dimensionless fields
uˆ = u/Uo, ˆp = p/pext, λˆ = λ/λ0, and ˆµ = µ/µ0. While µ0, λ0 and pext appear naturally in
the problem definition, we choose Uo to be some value of the order of maximum deformation
in the structure. In our analysis, we set Uo =
Lpext
µ0
, where L is a measure of the spatial
extent of the structure.
On substituting these dimensionless quantities in Eqs. (2.2) to (2.9), and working with the
time-harmonic version of these equations, we obtain the following set of non-dimensionalized
19

equations:
−uˆ − α1∇ · ˆ σˆ = fˆu in Ωs (2.12)
α4pˆ+ ∆ˆ ˆ p = 0 in Ωf (2.13)
uˆ = 0 on ΓD (2.14)
σˆ · nˆ = hˆ on ΓN (2.15)
pˆ = 1 on ΓDf (2.16)
∇ˆ pˆ · nˆ = 0 on ΓN f (2.17)
uˆ · nˆ = α5∇ˆ pˆ · nˆ on ΓI (2.18)
σˆ · nˆ = −α6pˆnˆ on ΓI (2.19)
where fˆu = α3fu is the dimensionless body force, hˆ =
hL
µ0U o is the dimensionless traction,
and the dimensionless version for the stress in the structure is given by
σˆ =
n
µˆ(∇ˆ uˆ + ∇ˆ uˆ
T
) + α2λˆ(∇ · ˆ uˆ)I
o
. (2.20)
Further, the dimensionless parameters, αi
, are
α1 =
µ0
ρsω2L2
α2 =
λ0
µ0
α3 =
1
ρsω2Uo
α4 =
ρfω
2L
2
κ
α5 =
pext
ρfω2LUo
α6 =
pextL
µ0Uo
Note that in the equations above the gradient operator ∇ has been replaced by its nondimensional counterpart ∇ˆ = L∇ (the scale length for the geometric model is also L).
2.1.1.0.1 Remark In subsequent sections, we will work only with dimensionless quantities; therefore, in order to keep the notation simple, we drop the hat sign (ˆ·) used to denote
20

dimensionless quantities from here on.
2.1.2 Weak form
Consider the space of trial solutions, S, and the space of weighing solutions, V, given by
S =

(u, p)| ui ∈ H1
(Ωs),u = 0 on ΓD; p ∈ H1
(Ωf ), p = 1 on ΓDf
V =

(w, q)| wi ∈ H1
(Ωs), w = 0 on ΓD; q ∈ H1
(Ωf ), q = 0 on ΓDf
Here, H1
is the Sobolev space of square integrable functions with square integrable derivatives.
The equivalent weak statement for Eqs. (2.10) to (2.19) can be written as follows: find
(u, p) ∈ S, such that
Bs(w,u; ψ) + α1α6C(w, p) + 1
α5
C(u, q) + Bf (q, p) = ⟨w, fu⟩ + ⟨w,h⟩ΓN
, ∀(w, q) ∈ V
(2.21)
where the bilinear forms Bs(·, ·; ψ), Bf (·, ·),⟨·, ·⟩, and C(·, ·) are defined as
Bs(w,u; ψ) = Re
Z
Ωs

α1µ(ψ)(∇w + ∇wT
) : ∇u + α1α2λ(ψ)(∇ · w)(∇ · u) − w · u

dΩ
Bf (q, p) = Re
Z
Ωf

∇q · ∇p − α4qp
dΩ!
C(u, q) = Re
Z
ΓI
u · nˆ q dΓ
⟨w, fu⟩ = Re
Z
Ωs
w · fu dΩ
⟨w,h⟩ΓN
= Re
Z
ΓN
w · h dΓ
Note that we choose to work with the real component of the bilinear forms instead of the
imaginary component. This choice is arbitrary since either component yields the same set

of equations.
The finite element approximation to the weak form is obtained by introducing S
h ⊂ S
and V
h ⊂ V as finite dimensional approximations to the original spaces. This leads us to
the following finite-dimensional problem: find (u
h
, ph
) ∈ Sh
, such that ∀ (wh
, qh
) ∈ Vh
,
Bs(wh
,u
h
; ψ
h
) + α1α6C(wh
, ph
) + 1
α5
C(u
h
, qh
) + Bf (q
h
, ph
) =

wh
, f
h
u

+

wh
,h

ΓN
.
(2.22)
2.1.3 Sensitivity analysis
We would like to derive the following estimate: If we perturb design field, ψ, by a small
variation, δψ, how would the functional J(ψ) change. This variation is given by
δJ = DψJ · δψ, (2.23)
where DψJ · δψ is the Gateaux derivative of the functional, given by,
DψJ · δψ =
d
dϵJ(ψ + ϵδψ)

ϵ=0
. (2.24)
One approach to determining this variation would involve computing the corresponding
change in the displacement and pressure fields, and then using these in the expression for
the functional. This approach would require solving a structural acoustics problem for every
variation δψ that we consider. This makes it impractical where many such calculations are
needed. An alternate approach, which involves using the adjoint of the structural acoustics
problem, avoids this expense and is described next.
We begin with writing an “expanded’ functional, which has the same form as J, but
releases the constraint on the fields u and p. This is given by
J˜(u, p, ψ) = F(u, p) + R(ψ). (2.25)
22

Next, we define a Lagrangian that appends to this functional the constraint equations:
L(w,u, ψ) = J˜(u, p, ψ)+Bs(w,u; ψ)+α1α6C(w, p)+Bf (q, p)+ 1
α5
C(u, q)−⟨w, fu⟩−⟨w,h⟩ΓN
.
(2.26)
When the functions u and p are chosen such that they are solutions to the weak form, we
note that all terms except the one involving the functional vanish from the definition of the
Lagrangian, and L = J˜ = J. Further, under the same constraints,
δJ = δL. (2.27)
Recognizing the dependence of L on all its fields, we note that
δL = DψL · δψ + (DuL · δu + DpL · δp) + (DwL · δw + DqL · δq). (2.28)
Evaluating the terms in the second parenthesis on the right hand side, we observe
DwL · δw + DqL · δq = Bs(δw,u; ψ) + α1α6C(δw, p) − ⟨δw, fu⟩ − ⟨δw,h⟩ΓN
+Bf (δq, p) + 1
α5
C(u, δq). (2.29)
The terms on the right hand side above are identical to the weak form of the constraint
equations, and are therefore equal to zero for any δw and δq, since u and p are chosen to
satisfy this constraint. Consequently from Eq. (2.28),
δL = DψL · δψ + (DuL · δu + DpL · δp). (2.30)
We now evaluate the terms within the parenthesis in the equation above,
DuL · δu + DpL · δp = DuF · δu + Bs(w, δu; ψ) + 1
α5
C(δu, q)
+DpF · δp + Bf (q, δp) + α1α6C(w, δp). (2.31)
23

Motivated by our goal to obtain a simple expression for δL, we select w and q to ensure
that the terms above will evaluate to zero for all δu and δp. This leads us to the following
problem: find (w, q) ∈ V, such that ∀(δu, δp) ∈ V,
Bs(w, δu; ψ) + 1
α5
C(δu, q) + α1α6C(w, δp) + Bf (q, δp) = −DuF · δu − DpF · δp. (2.32)
With this choice, from Eq. (2.30), we have,
δL = DψL · δψ = DψBs(w,u; ψ) · δψ + DψR(ψ) · δψ. (2.33)
Finally, from Eq. (2.27) and Eq. (2.33) we conclude that
δJ = DψBs(w,u; ψ) · δψ + DψR(ψ) · δψ, (2.34)
where the pair (u, p) satisfy Eq. (3.10) and the pair (w, q) satisfy Eq. (3.17).
The development above leads us to an efficient means for computing the change in the
functional corresponding to a change in the design variable. The steps are
1. Obtain (u, p) by solving the weak form Eq. (3.10).
2. Obtain (w, q) by solving the adjoint problem Eq. (3.17).
3. Use these in the Eq. (2.34).
Gradient vector. In practice, we discretize the forward and the adjoint problems using
the finite element. We also represent the design variable, ψ, with its finite-dimensional
counterpart, ψ
h
, where
ψ
h =
X
N
A=1
ψAϕA(x), (2.35)
24

and ψi are the coefficients used in the approximation of ψ
h
. Using this basis, the variations
in functional J may be written as
δJ =
X
N
A=1
≡gA
z }| {
DψBs(w,u; ψ) · ϕA + DψR(ψ) · ϕA

δψA. (2.36)
This leads us to the gradient vector g ∈ R
N , which contains information about the change
in the functional corresponding to a change in the coefficients of the design field.
2.1.4 Choice of functionals
In this work we consider two types of objectives. This first involves minimizing the square
of the pressure magnitude in a specified region within the acoustic domain. We refer to this
type of problem as a sound minimization problem. In this case
F(u, p) = Z
Ωr
|p|
2
dΩ (2.37)
Further in the adjoint Eq. (3.17),
DuF · δu = 0 (2.38)
DpF · δp = 2 Re

Z
Ωr
pδp ¯ dΩ
. (2.39)
The second problem is a vibration isolation problem, where we aim to minimize the
square of the displacement amplitude within a subdomain of the structure. In this case,
F(u, p) = Z
Ωr
|u|
2
dΩ. (2.40)

Further in the adjoint Eq. (3.17),
DuF · δu = 2 Re

Z
Ωr
u¯ · δu dΩ
(2.41)
DpF · δp = 0. (2.42)
2.1.5 Choice of regularization
As mentioned earlier, we express our prior belief regarding the design field through the
regularization term. By selecting different forms for this term we are able to design materials
that display spatial characteristics that we deem useful. In particular, we consider two cases.
In the first instance, we consider materials that are smoothly graded spatially. In this case,
R(ψ) = β
Z
Ωc
|∇ψ|
2
dΩ, (2.43)
and the Lam´e parameters are expressed as
µ(x) = µ0(0.1 + 0.9ψ(x)) (2.44)
λ(x) = λ0(0.1 + 0.9ψ(x)). (2.45)
In the expression for the regularization term, β is the regularization parameter.
With these choices, the explicit expression for the gradient vector is given by
gA = 0.9 α1
Z
Ωc

µ0∇wh
: (∇u
h + ∇u
hT
) + α2λ0(∇ · wh
)(∇ · u
h
)

ϕA dΩ
+ 2β
Z
Ωc
∇ψ
h
· ∇ϕA dΩ.
(2.46)
For the second case, we desire that the optimal material distribution have a binary
structure. In other words, we desire that the optimal distribution be comprised of two
material types. One with Lam´e parameters given by (µ0, λ0), and another softer material
with Lam´e parameters (0.1 × µ0, 0.1 × λ0). To achieve this, we replace the H1 regularization

term in Eq. (2.25) with the total variation (TV ) functional. Thus we have,
R(ψ) = β1
Z
Ωc
p
(|∇ψ|
2 + c) dΩ + β2
Z
Ωc
ψ (1 − ψ) dΩ. (2.47)
The first term in Eq. (2.47) is the T V regularization, which promotes sparsity in the |∇ψ|,
and c is a small positive number to ensure that this term has continuous derivatives when
|∇ψ| = 0. The second term is added to penalize deviations from ψ = 0, and ψ = 1; β1 and
β2 are regularization parameters. With this choice for regularization, the gradient vector is
given by
gA = 0.9 α1
Z
Ωc

µ0∇wh
: (∇u
h + ∇u
hT ) + α2λ0(∇ · wh
)(∇ · u
h
)

ϕA dΩ
+β1
Z
Ωc
∇ψ
h
· ∇ϕA
(|∇ψh|
2 + c)
1/2
dΩ + β2
Z
Ωc
(1 − 2ψ
h
)ϕA dΩ. (2.48)
The T V term attempts to minimize peaks and troughs in the solution by flattening out
most regions. As its name suggests, it only penalizes the total difference between maximum
and minimum values and not the rate of increase. The second penalization term, which
corresponds to β2, forces our solution to take binary values.
2.2 Numerical examples
In this section we apply the algorithms developed in the previous sections to two structural
acoustics problems. The first, which we refer to as the sound minimization problem, involves
optimizing the elastic properties of a structure so as to minimize the pressure intensity within
the acoustic domain enclosed by the structure. The second, which is referred to as the
vibration isolation problem involves optimizing the elastic properties of a sub-region of a
structure so as to minimize the vibration of another distinct sub-region. In both cases we
consider algorithms that produce optimal material properties that are smooth, and those
that have a distinct bi-material distribution. We use the multiphenics library [54] to solve
2

the structural acoustics PDEs and the L-BFGS-B algorithm in SciPy’s optimization package
[55] to solve the minimization problem.
2.2.1 Sound minimization
In this example, we minimize the sound intensity inside a region in the fluid that is enclosed
within a solid structure. This is an example of an interior structural acoustics problem.
We introduce the sound source as a time-harmonic external traction on the structure. The
structure itself is a two-dimensional elliptical shell with an outer semi-major axis of 1.25 m
and an outer semi-minor axis of 1 m, and a thickness of 0.25 m. The model setup for the
optimization problem is illustrated in Figure 2.2. The computational model is discretized
with 28,230 linear triangular finite elements (ensuring at least 10 elements per wavelength
in the fluid), with a total of 50,788 degrees of freedom (see Figure 2.3).
Figure 2.2: Schematic of the sound minimization problem. RoI is a circle with 0.2 m radius
shifted 0.4 m from the center. The width of control region is 0.15 m.
Our goal is to find the optimal material property distribution within the control region
Ωc, which would minimize the magnitude squared of the pressure inside the region of interest
(RoI) Ωr. The functional for this problem is given in Eq. (2.37). The material properties and
the parameter values used for the optimization problem are listed in Table 2.1. Further, the
damping effects in all our simulations have been kept to a minimum by choosing arbitrary
28

Figure 2.3: Finite element mesh used to solve the sound minimization problem.
small values of order 10e-8 for the viscosity coefficients (the complex components of κ, µ,
and λ).
Now, before we attempt to solve the optimization problem, it is pertinent to verify our
numerical implementation of the gradient vector. For this, we rely on a commonly used test
called the Taylor test. Referring to Eq. (2.1), where J represents the cost functional that
depends on the design variable ψ, we expand J in a Taylor series, for a small increment hδψ,
and arrive at
R = J(ψ + hδψ) − J(ψ) − DψJ · hδψ ≈ O(h
2
). (2.49)
Table 2.1: Material properties for the sound minimization problem.
Parameter Value Units
Bulk modulus κ 2.15 × 109 N/m2
Fluid density ρf 997 kg/m3
Lame’s 1st coef. λ0 6.59 × 1010 N/m2
Lame’s 2nd coef. µ0 7.69 × 1010 N/m2
Solid density ρs 7800 kg/m3
Forcing freq. ω 8000 rad/s
Normal traction |h| 0.1 kN/m2
29

Figure 2.4: Taylor test for the sound minimization problem with ψ = 1 − 0.1( x+1.25
2.5
). The
dashed lines represent slope = 2.
This formula can be used to determine the accuracy of a method used to compute the
variation DψJ · hδψ. This is done by first fixing h and δψ, and then computing the three
terms on the right hand side of Eq. (2.49) to evaluate R. This is repeated for different
values of h and the slope of the R versus h plot is evaluated. If this slope is close to 2 as h
becomes small, the Taylor test is considered to be successful. As can be seen in Figure 2.4,
our implementation passes the Taylor test.
2.2.1.1 Smoothly graded material
In order to ensure a smooth material property distribution we use the H1
regularization
term Eq. (2.43) and the corresponding formula for the gradient vector Eq. (2.46) to drive
our optimization algorithm.
The selection of the regularization parameter β is guided by principles that are used
to solve inverse problems. In particular, we sweep over the values of the regularization
30

Figure 2.5: Variation of the QoI as a function of the regularization parameter for the sound
minimization problem with smooth material properties.
Figure 2.6: Variation of the objective function as a function of iteration number for the sound
minimization problem (CV: smoothly graded material; Bimat: bi-material distribution).
parameter, and for each value we solve the optimization problem and determine the QoI.
Thereafter we generate a plot of the QoI versus the regularization parameter. In this plot,
31

solutions obtained using a low regularization parameter correspond to material fields that
are not smooth but have a small value of the QoI; while those with large values correspond
to smooth fields that have a relatively large value of the QoI. Given this, we select the largest
value of the regularization parameter that provides a sufficiently small QoI as the “optimal”
value. This plot is shown in Figure 2.5; from this plot we conclude that the optimal value is
β = 10−3
.
Figure 2.7: Distribution of displacement and pressure for the sound minimization problem;
QoI: 2.20e-1 (initial), 4.47e-4 (final).
As described in the previous section, for each value of the regularization parameter, we
solve the minimization problem using a gradient-based (L-BFGS) algorithm. We begin with
32

a uniform guess for the material parameter distribution (ψ = 1) and monitor the drop in
the objective function as a function of iteration number (see Figure 2.6). When this drop
stagnates, or when the gradient falls below a certain threshold, we terminate the algorithm.
For both the smoothly graded material and the bi-material configurations (considered in the
following subsection) we find that 80 iterations are sufficient to achieve this drop. Also in
both cases we achieve a significant drop in the objective function (more than a factor of 100).
In Figure 2.7, we have compared the initial and post-optimization distribution of the
displacement within the structure and pressure fields in the fluid. In the initial configuration
we observe that the response of the structure to the external traction is such that it sets up
standing waves within the fluid, and that this wave pattern attains a node close to the RoI.
In the optimized configuration, the displacement pattern in the structure appears to have
smaller wavelengths, which in turn modifies the standing wave pattern within the fluid. In
particular, it leads to a pattern with an anti-node close to the RoI, which drastically lowers
the pressure intensity in this region.
Figure 2.8: Optimal Young’s modulus distribution for the sound minimization problem with
smooth material properties.
The optimized Young’s modulus distribution within the structure, shown in Figure 2.8,
reveals how the smaller wavelengths within the structure are obtained. They are obtained
by lowering the structural stiffness and hence the wave-speed within the structure, and by
33

introducing spatial heterogeneity. From this figure we also observe that as stipulated by the
regularization term, the optimal material property distribution is relatively smooth.
Figure 2.9: Optimal distribution for three different starting guesses.
As the optimization is non-convex, the solution obtained here is that of a local minimum.
Therefore, we investigate the influence of the starting guess by looking at two other initial
distributions. In Figure 2.9, we have shown the initial guess in the left column and the final,
or converged state, in the right column. In case (a) the initial guess is uniform, in case (b) it
is given by ψ = 1−0.25(1 +x), and in case (c) it is given by ψ = 1−0.25(1−x). We observe
that while case (a) and (c) yield similar results, the final state for case (b) is different, with
34

Figure 2.10: Displacement and pressure fields of the optimal solutions corresponding to
the three initial guesses (a), (b), and (c); the QoI values are 4.47e-4, 2.10e-4, and 1.35e-4
respectively.
three (instead of two) regions of lower modulus. Here we observe that case (c) performs the
best. Further, it is instructive to look at the pressure and displacement fields for all three
cases in the final state (see Figure 2.10). While there are minor differences, the overall trend
35

in all the cases is similar leading to a zone of destructive interference in the RoI. The fact
that different material property distributions lead to similar structural acoustic response is
not very surprising since the inverse problem of determining material properties from the
elastic response of the system is known to have multiple solutions [56, 57, 58].
2.2.1.2 Bi-material distribution
We now consider the case where the material heterogeneity is confined to a distribution
composed of two types of materials: a stiff material with parameters µ0 and λ0, and a softer
material with parameters that are a tenth of these values.
Figure 2.11: Variation of the QoI as a function of the regularization parameters in the case
of binary distribution after 150 iterations.
We achieve this by using the regularization term in Eq. (2.47) and the corresponding
36

gradient that appears in Eq. (2.48). The parameters β1 and β2 are chosen after performing
a standard grid search (see Figure 2.11). The most suitable choice achieves the best compromise between the QoI value and the binary characteristic of the material distribution.
The displacement and pressure fields obtained at the optimal distribution of the material
parameters are shown in Figure 2.12 (β1 = 10−4
, β2 = 10−2
, c = 10−3
). Once again we
observe that the algorithm has attained a small value of pressure intensity in the RoI by
manipulating the displacement in the structure so as to produce an anti-node close to the
RoI. Similar to the smoothly graded case, the displacement field in the structure has smaller
wavelengths when compared to the initial field (see Figure 2.7 a and b).
Figure 2.12: Distribution of displacement and pressure for the sound minimization problem
at optimal material properties (bi-material case); QoI = 9.91e-5.
In order to examine the origin of this displacement field, we have plotted the Young’s
modulus distribution at four distinct iterations of the optimization algorithm in Figure 2.13.
We observe that the final distribution comprises of three relatively small, roughly semielliptical regions of the softer material embedded in the stiffer material. It is remarkable
that these soft regions are effective in significantly reducing the QoI. It is also interesting
to observe that the algorithm proceeds by first generating a modulus distribution that is
smooth, and then with increasing iterations sharpens the material property distribution.
37

Figure 2.13: Optimal Young’s modulus distribution for the sound minimization problem
with bi-material properties.
2.2.2 Vibration isolation
We now consider an exterior structural acoustics problem where the structure is surrounded
by an unbounded fluid domain. The structure is an evacuated cylindrical shell, which is 0.6
m thick and has an outer radius of 1 m. A time-harmonic cylindrical source is placed in
the fluid domain in the vicinity of the shell (see Figure 2.14). The unbounded domain is
modeled by the use of a perfectly matched layer (PML) [59].
The structure is further divided into three circular regions. In the outer region the
material parameters are held fixed. The middle region is the control region, where the
material parameters are allowed to vary, and the inner region is where the QoI is defined.
The QoI is given in Eq. (2.40) and is a measure of the magnitude of displacement within the
38

Figure 2.14: Setup for the vibration isolation problem. The location of source: 1.5 m from
center of shell; the thickness of RoI and the control region are 0.1 m and 0.3 m; the width
of PML is 0.7 m.
RoI. The material parameters for the problem are reported in Table 2.2.
The finite element mesh used to solve the problem is shown in Figure 2.15. It comprises
of 46,477 linear triangular finite elements with 63,166 degrees of freedom. Due to symmetry
we only discretize one half of the original domain.
As in the previous problem, we verify the numerical implementation of the gradient by
performing the Taylor test, and the results are shown in Figure 2.16. Once again, we obtained
the h
2 dependence of the residual, which confirms that our implementation is correct.
Table 2.2: Parameters for the vibration isolation problem.
Parameter Value Units
Bulk modulus κ 1.28 × 105 N/m2
Fluid density ρf 1.25 kg/m3
Lame’s 1st coef. λ0 6.59 × 1010 N/m2
Lame’s 2nd coef. µ0 7.69 × 1010 N/m2
Solid density ρs 7800 kg/m3
Forcing freq. ω 4000 rad/s
Source pext 0.1 kN/m2
39

Figure 2.15: Finite element mesh used to solve the vibration isolation problem.
Figure 2.16: Taylor test for the vibration isolation problem with ψ = 1 − 0.1( x+0.4
0.8
).
2.2.2.1 Smoothly graded material
We begin with uniform structural properties (ψ = 1) and track the objective function as a
function of iteration number in Figure 2.17. For both the smoothly graded material and the
bi-material configurations (considered in the following subsection) we find that 100 iterations
are sufficient to achieve a converged state. Also, in both cases we achieve a significant drop
40

in the objective function (more than a factor of 50).
Figure 2.17: Variation of the objective function as a function of iteration number for the
vibration isolation problem (CV: smoothly graded material; Bimat: bi-material distribution).
In Figure 2.18, we have the plotted real and absolute values of the pressure in the fluid
and the displacement field within the structure in the initial configuration. We observe that
there is significant coupling between the fluid and the structure, and within the RoI (the
inner part of the structure) we observe displacement amplitudes that are around 1 unit. We
also observe that the structure scatters the waves generated by the acoustic source in the
fluid and that there is a significant shadow zone behind the structure.
In Figure 2.19, we have plotted the solution fields obtained after optimizing the elastic
properties within the control region of the structure. By comparing these images with Figure 2.18, we observe that the pre- and post-optimization pressure fields are very similar, and
there are no discernible differences between the two. On the other hand, the magnitude of
the displacement within the structure is significantly smaller. Thus the vibration minimization algorithm has achieved its goal. We also note that this decrease is not restricted to the
RoI, but is spread through the entire structure.
In Figure 2.20 we have plotted the optimized distribution of the Young’s modulus in the
41

Figure 2.18: Initial (pre-optimization) distribution of displacement and pressure fields for
the vibration isolation problem with smooth material properties; QoI = 3.48e-2.
structure. We observe that within the control region the Young’s modulus is significantly
smaller. It begins at about 0.5 near the part that faces the acoustic source and smoothly
drops to 0.1 on the other end. This smooth variation of the Young’s modulus indicates that
the selection of H1
regularization has achieved its goal.
It is interesting to consider how the lower displacement amplitude within the RoI was
attained in this problem. By comparing Figure 2.18 and Figure 2.19 we note that the displacement field in the entire structure was lowered through optimization. This indicates
that the optimized structure does not couple with the acoustic excitation as efficiently as
the structure in its initial configuration. Implying thereby that none of the natural frequencies of the optimized structure are close to the frequency of acoustic excitation. Thus the
optimization algorithm has achieved vibration isolation by moving the natural frequencies
42

Figure 2.19: Final (post-optimization) distribution of displacement and pressure fields for
the vibration isolation problem with smooth material properties; QoI = 6.99e-4.
Figure 2.20: Optimal distribution of Young’s modulus within the structure for the vibration
isolation problem.
of the structure away from the frequency of excitation.
43

2.2.2.2 Bi-material distribution
We consider the vibration isolation problem and require the material distribution to conform
to two likely values: one corresponding to a stiff material with parameters (µ0, λ0), and
another corresponding to a soft material with parameters that are tenth of these values.
This is achieved by replacing the H1
regularization with a TV term and another term that is
minimal at these values Eq. (2.47). The optimized pressure and displacement fields obtained
from this approach are shown in Figure 2.21. We observe that these fields are quite similar
to the fields obtained for the smoothly graded material, and when compared with the initial
state, the vibration amplitude is significantly lower.
Figure 2.21: Final (post-optimization) distribution of displacement and pressure fields for
the vibration isolation problem with bi-material property distribution; QoI = 5.12e-4.
In Figure 2.22 we have plotted the variation in the material property of the structure
as a function of iteration number during the optimization procedure. We note that the
44

optimization algorithm achieves its goal by first generating two regions with lower moduli,
and then sharpening the transition between these regions and the surrounding material, and
pushing the values in these regions toward the lower limit.
Figure 2.22: Optimal distribution of Young’s modulus within the structure for the vibration
isolation problem with a bi-material distribution.
Parts of this chapter are taken from the reference:
Harisankar Ramaswamy, Saikat Dey, and Assad A. Oberai. Material parameter optimization for interior
and exterior fluid-structure acoustic problems. International Journal for Numerical Methods in Engineering,
121(24):5568–5589, 2020.
45

Chapter 3
Shape Estimation
3.1 Problem definition
Let us consider a typical coupled structural acoustic system consisting of two regions: a
linear elastic solid domain Ωs and an acoustic fluid Ωf . The dynamics of these two regions
are coupled across an interface ΓI . A generalized pictorial depiction of the problem domain
is shown in Figure 3.1. Our goal is to optimize the shape of the interface such that a given
quantity of interest, defined by a functional F, is minimized. This functional is defined over
the region of interest Ωr, which can either lie entirely on the solid or fluid domain. We solve
this optimization problem incrementally by making infinitesimal deformations of O(ϵ) to the
interface boundary along the sensitivity direction, all the while maintaining the mechanical
equilibrium of the system.
The steady-state response of the system due to a time-harmonic excitation, is determined by solving the Navier-Cauchy equation in Ωs and the Helmholtz equation in Ωf . The
46

Figure 3.1: Illustration of the domain of optimization problem. The gray region illustrates
the solid, blue region illustrates the fluid domain, and the violet region is the region of
interest. Γϵ
I
is the small incremental deformation of ΓI
governing equations are as follows:
ρsω
2u + ∇ · σ = 0 in Ωs (3.1)
ρfω
2
κ
p + ∇2
p = 0 in Ωf (3.2)
u = 0 on ΓD (3.3)
σ · nˆ = h on ΓN (3.4)
p = pext on ΓDf (3.5)
∇p · nˆ = 0 on ΓN f (3.6)
ω
2u · nˆ =
∇p
ρf
· nˆ on ΓI (3.7)
σ · nˆ = −pnˆ on ΓI (3.8)
Here, the state variables u and p are the complex displacement and pressure fields; σ is
the Cauchy’s stress tensor which is given by the constitutive relation σ = µ(∇u + ∇u
T
) +
47

λ(∇ · u)I, with λ and µ being the first and the second Lam´e parameters; ρs and ρf are the
densities of solid and fluid at rest; either h or pext can act as the driving excitation, ω is the
angular frequency of the excitation; nˆ is the outward normal. Further, we allow the elastic
moduli to be complex to account for viscous dissipation.
The optimization problem is posed as
Minimize
ΓI
F(u) = Z
Ωr
f(u(ΓI )) dΩ. (3.9)
where, u satisfies Eqs. (3.1) to (3.8) and the parameters representing the shape of ΓI are
the control variables. Without loss of generality we assume Ωr ⊂ Ωs, therefore F here is
some arbitrary functional of the corresponding state variable u. Complementarily, in an
example problem to demonstrate our methodology, we solve the noise abatement problem
from Chapter 2 by considering Ωr to be inside Ωf . Further, the shape of ΓI is represented
by the position of mesh nodes on the interface alone and is agnostic to any other parametric
representation; this is the key feature of our method.
Remarks. In our analyses we make the following assumptions:
1. There can be multiple interfaces in the problem definition, however, each deformable
interface is considered to be a closed manifold (in the sense that it encloses a simply
connected domain in R
n
). Note that the enclosed region can either be solid (which is
suspended in a fluid) or fluid (which is contained within a structure).
2. The deformable interfaces do not have segments which are fixed. Else, special care is
required in imposing additional constrains at the junction points in order to maintain
regularity.
3. The interface is sufficiently smooth so that the curvatures, more specifically the tangential divergence of the outward-normal, can be computed at any point on the interface.
This assumption is critical for computing the shape sensitivity as we will see in Section 3.2
48

4. The region of interest is defined sufficiently far from the interface such that it does not
touch or cross the deforming interface during the optimization process, ΓI ∩ ∂Ωr = ∅.
3.1.1 Weak form
We start our analysis by writing the non-dimensionalized weak expression for our primal
problem Eqs. (3.1) to (3.8). For this, we scale the equations in the same manner as described
in Section 2.1.1.
Find (u, p) ∈ S, such that
Bs(w,u) + α
2
1Bf (q, p) + α1C(w, p) + α1C(u, q) = α1 ⟨w,h⟩ΓN
, ∀(w, q) ∈ V (3.10)
where the first two terms denote the contributions from the domain integrals on Ωs and Ωf
respectively, the next two correspond to the two interface conditions, and the right most is
the external traction force term. The expansion of each of the terms in Eq. (3.10) are as
follows:
Bs(w,u) = Re
Z
Ωs

α1(∇w + ∇wT
) : ∇u + α1α2(∇ · w)(∇ · u) − w · u

dΩ
Bf (q, p) = Re
Z
Ωf

∇q · ∇p − α4qp
dΩ!
C(u, q) = Re
Z
ΓI
u · nˆ q dΓ
⟨w,h⟩ΓN
= Re
Z
ΓN
w · h dΓ
Further, the set of trial solutions S and weighting solutions V are given by
S =

(u, p)| ui ∈ H1
(Ωs),u = 0 on ΓD; p ∈ H1
(Ωf ), p = 1 on ΓDf
V =

(w, q)| wi ∈ H1
(Ωs), w = 0 on ΓD; q ∈ H1
(Ωf ), q = 0 on ΓDf

3.1.2 Shape sensitivity analysis
For solving the optimization problem stated in Eq. (3.9), we have to compute the derivative
of F(u) with respect to the shape of the interface ΓI . Although ΓI does not appear explicitly
in F, the steady state solution u is dependent on ΓI . We follow the classical approach used
the shape optimization literature for computing the shape derivative.
Let ϕϵ be a linear map which takes any point from the current configuration of the
interface to its deformed configuration, Γϵ
I = ϕϵ(ΓI ). We define ϕϵ as
ϕϵ(x) = {xϵ
|xϵ = x + ϵvn(x), x ∈ ΓI} (3.11)
where vn = vnˆ is the instantaneous rate of deformation along the normal direction of the
interface. We only consider the deformation along the normal direction because the presence
of tangential component would merely redistribute the points within the deformed interface.
Note that this makes the map ϕϵ(v) : ΓI → Γ
ϵ
I unique for a given (ΓI , Γ
ϵ
I
).
Now,the shape derivative of F is defined as
DϕF · v = limϵ→0
F(u(ϕϵ(v))) − F(u(ϕ0(v)))
ϵ
(3.12)
where DϕF · v is the Gateaux derivative of F with respect to ϕϵ
in the direction of v.
Now, we need to compute DϕF, which would give us the steepest descent direction
to solve the minimization problem. One approach would be to numerically compute the
derivative by using the chain rule followed by finite difference approximation,
DϕFk = DϕF.vk = DuF(u(ΓI )) ·
du(ϕϵ(vk))
dϵ

ϵ=0
k ∈ N (3.13)
where (DϕFk)
K
k=1 is the K-dimensional vector equivalent of DϕF, K is the number of nodes
(of the mesh) in the discretized boundary, and vk is the perturbation of the k
th node keeping
50

the remaining nodes fixed. Consequently, the number of primal problems that needs to be
solved is proportional to K, which can be computationally costly when K is large.
In our analysis, we use an alternate approach based on the adjoint method to circumvent
the issue of having to explicitly compute the derivative of u(ϕϵ) with respect to ϕϵ
. In
this method, we treat the state variables (u, p) along with ϕϵ as independent variables, by
appending the governing equations to the cost functional F using Lagrangian multipliers.
We define the Lagrangian L as
L(w, q,u, p; ϕ) = F(u)+Bs(w,u; ϕ)+α
2
1Bf (q, p; ϕ)+α1C(w, p; ϕ)+α1C(u, q; ϕ)−α1 ⟨w,h⟩ΓN
(3.14)
Here, Bs(·, ·) and Bs(·, ·) are dependent on ϕ, since Ωs and Ωf are dependent on ϕ. Variables
(w, q) ∈ V are the Lagrangian multipliers corresponding to the constraint equations for
(u, p). Note that the terms on the right-hand side, barring the first term, gives us the weak
form Eq. (3.10); therefore, L is identically equal to F.
Thus,
δL = δF (3.15)
Next, we expand δL using the chain rule
δL(w, q,u, p; ϕ) = DϕL · δv + (DwL · δw + DqL · δq) + (DuL · δu + DpL · δp) (3.16)
The expression inside the first parenthesis is identical to the weak form Eq. (3.10) and hence
equates to zero ∀(δw, δq) ∈ V. In addition, we simplify the equation for δL further, by
setting the expression inside the second parenthesis also to zero. This gives rise to the
adjoint problem which is formally stated as follows:
51

Find (w, q) ∈ V such that ∀(δu, δp) ∈ V,
(DuF(u) · δu + Bs(w, δu; ϕ0) + α1C(δu, q; ϕ0)) + (α
2
1Bf (q, δp; ϕ0) + α1C(w, δp; ϕ0)) = 0
(3.17)
Now, with the choice of u, p, w, and q obtained by solving Eqs. (3.10) and (3.17), we
observe that the Eq. (3.14) reduces to
δL = DϕL · δv
= DϕBs(w,u; ϕ) · δv + α
2
1DϕBf (q, p; ϕ) · δv + α1DϕC(w, p; ϕ) · δv
+ α1DϕC(u, q; ϕ) · δv
(3.18)
Recognize that in the equation above, we need to compute the derivative of integrals whose
domain of integration is moving; so to proceed further, we rely on the following identities
derived from the Leibniz integral rule.
Proposition 1. Domain integrals with moving boundary.
F =
Z
Ωϵ
f dΩ
δF =
Z
Ω0
∂f
∂ϵ δϵ dΩ + Z
Γ0
fnˆ · vδϵ dΓ
(3.19)
Similarly,
Bs(w,u; ϕ) = Z
Ωϵ
s
bs(w,u) dΩ
DϕBs(w,u; ϕ) · δv =
Z
ΓI
bs(w,u)δv dΓ
(3.20)
52

Proposition 2. Boundary integrals on the moving boundary.
G =
Z
∂Ωϵ
g dΓ
=
Z
∂Ωϵ
gn · n dΓ
=
Z
Ωϵ
∇ · (gn) dΩ
=
Z
Ωϵ
{g∇ · n + n · ∇g} dΩ
δG =
Z
Ω0
∂(g∇ · n + n · ∇g)
∂ϵ δϵ dΩ + Z
Γ0
(g∇ · n + n · ∇g)nˆ · vδϵ dΓ
=
Z
Ω0
∇ · (
∂g
∂ϵ n) dΩ + Z
Γ0
(g∇ · n + n · ∇g)nˆ · vδϵ dΓ
=
Z
∂Ω
∂g
∂ϵ dΓ + Z
Γ0
(g∇ · n + n · ∇g)nˆ · vδϵ dΓ
(3.21)
Similarly,
C(w, p; ϕ) = Z
Γ
ϵ
I
c(w, p) dΓ
DϕC(w, p; ϕ) · δv =
Z
ΓI
{c(w, p)∇ · nˆ + nˆ · ∇c(w, p)} · δv dΓ
(3.22)
where, nˆ · ∇c(w, p) = (nˆ · ∇p)(w · nˆ) + p(nˆ · ∇(w · nˆ))
Finally, using the identities derived above, we obtain:
DϕF · δv = δL
=
Z
ΓI
bs(w,u)δv dΓ − α
2
1
Z
ΓI
bf (q, p)δv dΓ
+ α1
Z
ΓI
{∇c(u, q) · nˆ + c(u, q)∇ · nˆ}δv dΓ
+ α1
Z
ΓI
{∇c(w, p) · nˆ + c(w, p)∇ · nˆ}δv dΓ
(3.23)
where nˆ is the normal pointing outwards with respect to Ωs.
53

3.2 Numerical approximation of gradients
We observe that the expression for the shape derivative in Eq. (3.26) is of the form:
DϕF · δv =
Z
ΓI
G(u, w, p, q)δv dΓ (3.24)
The above equation is true for an arbitrary δv; thus we interpret G as the variational derivative of the functional F with respect to function v, in terms of the L
2
inner product.
δF(u; ϕ)
δv = G(u, w, p, q) (3.25)
where,
G(u, w, p, q) = bs(w,u) − α
2
1
bf (q, p) + α1{∇c(u, q) · nˆ + c(u, q)∇ · nˆ}
+ α1{∇c(w, p) · nˆ + c(w, p)∇ · nˆ}
(3.26)
The function G is defined everywhere on the interface ΓI . Hence, in order to solve the
minimization problem Eq. (3.9) through gradient descent, one has to make small deformations to the boundary along the negative direction of G, or v = −G. However, this poses
another challenge; it is difficult to numerically compute G due to the fact that it involves
computing the gradient ∇(·) and divergence ∇ · (·) of quantities which are defined only on
the boundary, like ∇ · nˆ and ∇(w · nˆ). Inaccurate estimates of these quantities can lead
to a highly irregular (distorted) boundaries (with kinks) within few iterations of boundary
movement. Therefore, it is crucial to device a clever way of computing G that is accurate,
efficient, and robust.
Now to solve this problem, we propose a novel method to construct a smooth extension
of the normal vector into an auxiliary domain adjacent to the interface. Once we have the
extended normal field, the computation of G becomes straightforward. The procedure to
54

construct such a field is described as follows:
Step 1: We start by considering an domain Ωa which extends to either side of the interface
as shown in Figure 3.2. The two boundaries of this region can be arbitrary and can
conform to the local mesh pattern (including zig-zag patterns), since ultimately we
are only interested in the values on the interface. Nonetheless, we ensure that the
interface approximately lies in the middle and the width of Ωa is approximately
10 times the local mesh size (or more); this would be enough to minimize spurious
effects from the boundary.
Figure 3.2: Construction of auxiliary domain around the closed interface.
Step 2: Next, we solve a simple heat conduction problem on Ωa with Dirichlet conditions
on its boundaries and also along the interface. This problem is given by
−kc∇ · ∇s = 0 in Ωa (3.27)
s = −s0 on Γ1 = ∂Ωa ∪ ωs (3.28)
s = s0 on Γ2 = ∂Ωa ∪ ωf (3.29)
s = 0 on ΓI (3.30)
where, s is the unknown field. Suitable values for parameters kc and s0 can be
55

chosen; we use kc = 1 and s0 = 10 in our work. The key idea behind enforcing
homogeneous Dirichlet condition on ΓI is to ensure that the gradient of s is in the
same direction as the interface normal, i.e. ∇s/|∇s| = nˆ. It can be shown that a
non-zero tangential component of ∇s would violate the constraint s(x) = 0 ∀x ∈
ΓI .
Step 3: Once we solve for s, we compute the normal field in a weak sense (since ∇s may
not be defined on ΓI ), using the equation:
nˆ =
∇s
√
∇s · ∇s + ϵ
in Ωa (3.31)
where ϵ << 1 is a small positive number to prevent division by zero.
Step 4: Finally, we compute the term ∇ · nˆ by solving the variational problem with an add
regularization term for better numerical performance: find k ∈ H1
(Ωa), such that
Z
Ωa

vk + β∇v · ∇k

dΩ = Z
Ωa
v∇ · nˆ dΩ ∀v ∈ H1
(Ωa) (3.32)
The β-term penalizes sharp changes in curvature and suppresses the formation of
dendrite-like projections; values of O(10−3
) are used for β.
It may be noted that, the regularizing formulation in Eq. (3.32) is used to evaluate all
the quantities that involve partial derivatives (such as ∇ · nˆ, nˆ · ∇p, ∇(w · nˆ)) in Eq. (3.26).
Nonetheless, the total time taken for solving all such extension/approximation problems on
Ωa is insignificant when compared to solving a single primal or adjoint problem.
The overall procedure for computing the gradient vector for shape optimization can be
summarized as follows:
1. Solve for the primary/state variables (u, p) using Eq. (3.10).
2. Solve for the adjoint variables (w, q) using Eq. (3.17).
5

3. Obtain a smooth extension of the interface normal on an auxiliary domain Ωa.
4. Compute the interface deformation rate/gradient vector using Eq. (3.26).
Remarks. The gradient vector thus computed using the approach described in this
section, is used to move the interface nodes in each optimizer iteration. However, moving
just the interface nodes will result in a poor mesh quality near the interface; therefore, it is
necessary to deform the entire mesh such that it conforms to the new interface boundary.
In this work, we solve an elasticity-type problem, where we compute the displacement of
internal nodes, given the displacement at the interface. This is described by the following
equations:
∇ · σ(v) = 0 in Ωs ∪ Ωf (3.33)
v = 0 on Γ\ΓI (3.34)
v = Gnˆ on ΓI (3.35)
where v is the nodal displacement field that smoothly decays to zero away from ΓI and σ is
the Cauchy stress tensor.
3.3 Numerical examples
We illustrate our shape optimization methodology through two problems. First, we tackle a
benchmark test in linear elasticity to validate our implementation against a known solution.
Subsequently, we address a sound-abatement problem in a structural acoustic scenario.
3.3.1 Benchmark problem: Hole in a plate
A commonly known problem in shape optimization involves reducing stress concentrations
in a structure such as around a hole in a plate. As described in [31], the goal is to find
the optimal shape of the hole that minimizes stress concentration under an area constraint.
57

In this case, the known solution is a circular hole of equivalent area (of the initial hole
geometry).
We consider a 3 × 3 square plate that initially contains a square hole with a diagonal
length of 1 unit. Around this hole, we define an auxiliary region whose outer boundary is
offset 1 unit from the hole boundary. We apply uniform unit loads in both the vertical and
horizontal directions with roller supports so that no translational degrees of freedom remain.
We assume unit Young’s modulus and a Poisson’s ratio of 0.3. The linear elasticity equations
are solved to compute the displacement field u in the structure.
To mitigate stress concentrations, it suffices to minimize the structure’s compliance.
Mathematically, the compliance F(u) can be written as
Minimize
ΓI
F(u) = Z
ΓN
u · f dΓ. (3.36)
where ΓI is the hole boundary along which shape changes occur, and ΓN is the Neumann
boundary where the external normal loads of magnitude f are applied. Intuitively, minimizing compliance encourages the structure to adopt a geometry that distributes stress more
uniformly, reducing peak stress concentrations near the hole.
Pre-optimization Post-optimization
Figure 3.3: Shape optimization on an elastic plate. The auxiliary region is highlighted in
yellow on the left figure (initial configuration) and in green on the right figure (optimized
configuration).
Figure 3.3 shows the initial geometry with the auxiliary region around the hole high58

lighted in yellow. The plate is meshed with refined elements near the corners. We compute
shape derivatives using our proposed method and then iteratively update the shape using a
Lagrange multiplier–based approach. Notably, we impose node displacement constraints on
the moving boundary in a weak manner and solve the node displacement vector for the entire
mesh at each iteration. Here, imposing a Neumann boundary condition serves to regularize
the perturbations and prevents the excessive shape distortions that would otherwise occur
near sharp corners.
The outcome of the optimization appears in Figure 3.3 (right), indicating the new hole
boundary. As expected, the optimized hole converges to a circular shape, consistent with
established results.
3.3.2 Sound abatement problem
Having verified our shape optimization procedure on a well-known benchmark, we next consider the sound minimization example problem of Chapter 2. In that problem, we minimize
the sound intensity inside a region in the fluid that is enclosed within a solid structure due
to oscillating traction loads acting on the structure’s exterior surface. Here, our goal is to
find the optimal shape of the interface, which would minimize the magnitude squared of the
pressure inside the region of interest Ωr. The loss functional F for this problem is given by,
Minimize
ΓI
F(p) = Z
Ωr
pp¯ dΩ. (3.37)
In Figure 3.4 we show a schematic of the problem geometry along with the constructed
auxiliary region. The material properties and the parameter values used for the optimization
problem are listed in Table 2.1.
The reliability of our optimization algorithm hinges on the accuracy of the gradient
computation which in turn depends on the curvature calculations. Therefore we validate the
curvature calculations for the initial elliptical interface geometry (see Figure 3.5) obtained
59

Figure 3.4: Schematic of the sound minimization problem. The RoI is a circle with 0.2 m
radius shifted 0.4 m from the center. The width of the auxiliary region is 0.2 m.
a. Temperature field b. Normal vector c. Curvature
Figure 3.5: Computation of the extended normal field and subsequently the curvature field.
using the approach with analytical calculations. The red dots on Figure 3.5(c) are the
points where we compare the interface curvature with the theoretical value. We observe a
percentage error of 0.7% at the top most point on the interface and 2.1% at the right most
point.
For solving the optimization problem we use a gradient descent algorithm, which uses a
constant step size, coded in python. We start with the initial shape/mesh, feed the gradient
vector computed using the method described in the previous section, and monitor the loss
function as well as the mesh deformation in each iteration. The optimization algorithm
is stopped when the loss function stagnates at certain value or when there is excessive
deformation in the shape. We found that 20 iterations were sufficient for the algorithm to
60

Figure 3.6: Variation of the objective function as a function of iteration.
reach convergence while achieving 2.5 orders of magnitude reduction in the loss function.
The lost history (as function of iteration) is shown in Figure 3.6 and the optimized interface
geometry is shown in Figure 3.7.
Figure 3.7: Optimal shape of the interface along with deformed mesh of the solid region.
In Figure 3.8, we have compared the initial and post-optimization distribution of the
displacement within the structure and pressure fields in the fluid. We observe that the
external excitation leads to standing waves in the system, as one would expect, since we
are looking at the response of a closed undamped system. The optimization algorithm then
61

a. Initial fields
b. Post-optimization fields
Figure 3.8: Distribution of displacement and pressure for the sound minimization problem;
QoI: 2.44e-1 (initial), 7.19e-4 (final).
modifies the interface such that it leads to a standing wave pattern with an anti-node close
to the RoI, thereby drastically lowering the pressure intensity in this region.
Hence, this example demonstrates how the shape of the interface can effectively influence
the coupling behavior and consequently the performance of the system.
62

Chapter 4
Topology Estimation
4.1 Problem definition
In this chapter, we shift focus from design optimization to a different class of problems,
namely inverse problems. Specifically, we address inverse scattering, which aims to infer the
internal composition or topology of unknown scatterers by observing how an incident wave
interacts with them. As shown in Figure 4.1, an incoming wave (e.g., acoustic) illuminates
(through sound) a target region (represented by a dashed circle), and the scattered waves
provide indirect information about the shape, arrangement, and material properties of the
scatterers. Here, we concentrate on topology estimation, having already covered material
Figure 4.1: Inverse acoustic scattering problem.
63

and shape estimation in previous chapters. Given that inverse problems are often ill-posed
and that measurements may be sparse or noisy, uncertainty undoubtedly becomes central
to the reconstruction. To manage this, we adopt a conditional generative framework where
we generate realizations (samples) from the posterior distribution, enabling us not only to
identify the most plausible topologies but also to quantify the associated uncertainty.
Although acoustics is our primary application domain, the methodologies presented here
apply broadly to other mechanical and optical systems, such as elastography and optical
coherence tomography, where waves or mechanical loads probe hidden structures. In fact,
we include additional benchmark problems from these fields to demonstrate the efficiency
and robustness of our approach.
4.2 Conditional generative problem
The conditional generative problem involves generating samples of X from the conditional
distribution PX|Y (x|y), given an instance of Y = y. To formally define the problem, let
us consider the random field variable X ∈ R
Nx to represent the field to be inferred and
Y ∈ R
Ny to denote the noisy measurements. First, we generate samples x
(i)
from the prior
distribution PX(x). In real-life applications, the prior information typically comes in the
form of data samples from past observations, expert knowledge, or informed guesses, providing insights into the probable values or distributions that X may assume. Next, we generate
corresponding samples y
(i) by applying a measurement operator M(·), which represents the
composition of the forward map and the assumed noise model. This is equivalent to generating samples from the likelihood PY |X(y|x). The resulting dataset S = {x
(i)
, y
(i)} is
effectively drawn from the joint distribution PX,Y . The goal, given this dataset S and a
specific instance Y = y, is to generate samples from the posterior distribution PX|Y (x|y).
This task is central to probabilistic inference problems in general.
To solve this problem of computing the posterior, one approach is to use conditional
64

generative adversarial networks (cGANs), as proposed in [60]. These models provide a fast
way to generate samples after training and are effective at sampling from high-dimensional
distributions. However, adversarial training poses challenges due to the lack of a well-defined
loss function to monitor convergence. Additionally, cGANs can suffer from a lack of diversity
in the generated samples, which may not fully represent the true posterior distribution.
Another limitation of GAN-based models is their susceptibility to catastrophic forgetting,
where previously learned tasks are forgotten while learning new ones [61]. Furthermore,
they are known to exhibit generalization issues, particularly when applied to datasets not
encountered during training.
A second approach involves using Markov Chain Monte Carlo (MCMC) methods, particularly the Metropolis-Hastings algorithm, which is used for sampling from high-dimensional
distributions. This method requires knowledge of a function proportional to the posterior
density, typically involving the computation of the likelihood probability and the prior probability. The prior distribution can also be learned using GAN-based methods [62], which
map realizations from a known latent distribution, such as a standard Gaussian, to samples
from the prior. However, traditional MCMC sampling is computationally expensive, requiring significant time to converge, and often struggles to fully explore the posterior space.
Moreover, this method requires explicit knowledge of the likelihood density function, which
may not always be readily available.
In this work, we rely on a third approach that strikes a balance between the two: score
based diffusion models (SDMs) [63]. Similarly to cGANs, we learn the posterior distribution
from a data set that contains samples from the joint distribution. However, instead of direct
sampling, we employ MCMC sampling similar to the second approach, but based on annealed
Langevin dynamics.
Remark Throughout this work, we restrict attention to two-dimensional settings, treating
both the measurements and the inferred fields as image arrays.
65

4.3 Score based diffusion models
The goal of score-based generative models is to sample from a target data distribution by
progressively transforming a sample from an easy-to-sample distribution, such as a highdimensional Gaussian. This process is based on the observation that sequentially corrupting
samples of any distribution with increasing levels of noise eventually produces samples drawn
from a Gaussian distribution. Score-based generative models aim to reverse this process,
enabling the generation of samples from the target distribution by effectively denoising a
Gaussian sample. This generative process is based on Langevin dynamics and is governed
by a stochastic differential equation (SDE):
dxt = α(t)∇x log P(x, t) dt +
p
α(t) dWt
, (4.1)
where the first term, α(t)∇x log P(x, t) dt, is referred to as the drift term, and P(x, t) represents density function of X as it evolves. Here, the gradient of the logarithm of the density
function is called the score of the density function; hence the name score-based diffusion.
Note that P(x, t) can either be a marginal or conditional probability density of X depending on the application; this would not change the dynamics. The second term, p
α(t) dWt
,
represents the stochastic component with dWt denoting a standard Wiener process (or Brownian motion). Figure 4.2 illustrates this process, where the reverse SDE transforms noise
x(t = 0) into data x(t = T), effectively approximating the original data distribution.
The noising-denoising process has a strong connection with the diffusion equation, which
is why these models are referred to as diffusion models. This connection, along with the
specifics of Langevin dynamics, is discussed in detail in Section 4.3.1. Following this, Section 4.3.2 delves into learning the score function, where we employ the denoising score
matching objective to train the network. Finally, guidelines for hyperparameter selection
and design choices for the score network are also discussed.
66

Figure 4.2: Illustration of the forward and reverse SDE processes in diffusion models.
Figure 4.3: An example of 1-dimensional forward diffusion process which maps a complex
multi-modal distribution to a Gaussian. Each colored cross represents a sample, with identical colors indicating how a single point transforms under increasing noise through the SDE.
4.3.1 Langevin dynamics
Langevin dynamics is a widely-used mathematical framework in molecular physics to model
the motion of particles in a medium, such as fluids. It describes the position of a particle
as a function of time under the influence of deterministic forces (e.g., conservative forces,
damping) and random forces (e.g., interparticle collisions). In the context of generative
models, we use Langevin dynamics for modeling the process of incrementally adding noise to
samples (forward process) and subsequently deriving the reverse process, which corresponds
to denoising.
We begin by looking at the forward process. The forward process incrementally corrupts
data samples with Gaussian noise such that samples from the original distribution progres67

sively transform into samples from a Gaussian distribution. This process is formulated as a
stochastic differential equation (SDE):
xt+∆t = xt +
p
γ¯t∆t w, 0 ≤ t ≤ T (4.2)
where xt
is the progression of the sample at time t and w ∼ N (0, 1) represents Gaussian
noise. Figure 4.3 illustrates the forward process where samples x
(i)
0
from the initial data distribution PX(x) (left) are finally transformed into a Gaussian distribution N (µx, σ2
) (right).
From another perspective, each SDE has a corresponding Fokker-Planck [64] equation that
describes the evolution of the probability density function over time. Consequently, the SDE
in Eq. (4.2) is equivalent to the following partial differential equation (PDE):
∂pX
∂t =
γ¯(t)
2
∇2
x
pX, (4.3)
where pX is the probability density function of xt
. This process begins with the initial
condition pX(x, 0) = PX(x)
This process evolves the data distribution PX(x) towards a Gaussian distribution N (µx, σ2
T 1)
centered around µx =
1
N
PN
i=1 x
(i)
0
as t → T for sufficiently large σT where,
σ
2
T =
Z T
0
γ¯(t)dt (4.4)
Figure 4.4 illustrates this transformation process for a one-dimensional example. Notably,
this forward SDE is equivalent to the diffusion equation, whose analytical solution can be
obtained using Green’s functions, which are Gaussian kernels.
To reconstruct the original data distribution from noise, we formulate the reverse process.
We start by writing the reverse-time probability density ¯p(x, t), which evolves backward in
68

Figure 4.4: A schematic illustrating the forward diffusion process in one dimension across
various time (or noise) levels σ(t). The original multi-modal distribution PX(x) at t = 0
gradually transitions into smoother distributions as noise increases. Eventually, at t = T,
the distribution becomes close to a simple Gaussian N (µX, 5
2
).
time. We have,
∂p¯X(x, t)
∂t = −
∂
∂tpX(x, T − t) (4.5)
p¯X(x, t) = pX(x, T − t) (4.6)
Substituting this into Eq. (4.3), the reverse diffusion equation becomes,
∂p¯X
∂t = −
γ¯(T − t)
2
∇2
p¯X (4.7)
p¯X(x, 0) = N (µx, σ2
T 1) (4.8)
where γ(t) = ¯γ(T − t); 0 ≤ t ≤ T. However, this negative diffusion equation does not lead
69

to stochastic dynamics, so we apply the following manipulation trick:
For some ε ≥ 0, we have
∂p¯X
∂t = −
γ(t)
2
∇2
p¯X ±
εγ(t)
2
∇2
p¯X (4.9)
Next, we regroup the terms after adding and subtracting the ε term.
∂p¯X
∂t = −
γ(t)
2
(1 + ε)∇2
p¯X +
εγ(t)
2
∇2
p¯X (4.10)
In the first term, we multiply and divide by ¯pX,
∂p¯X
∂t = −∇ ·
γ(t)
2
(1 + ε)
∇p¯X
p¯X
p¯X

+
εγ(t)
2
∇2
p¯X (4.11)
We have,
∂p¯X
∂t = −∇ ·
γ(t)
2
(1 + ε)∇ log |p¯X| p¯X

+
εγ(t)
2
∇2
p¯X (4.12)
Thus,
∂p¯X
∂t = −∇ ·
γ(t)
2
(1 + ε) s(x, t) ¯pX

+
εγ(t)
2
∇2
p¯X (4.13)
where is s(x, t) is the score function.
The equivalent stochastic process for Eq. (4.13) is given by,
dxt = γ(t)s(xt
, t) dt +
p
γ(t) dWt
, (4.14)
where ε = 1 is assumed, following [63]. Remark: for ε = 0 we have a what we call probabilistic
ODE. Remarkably, when ε = 0 the equation reduces to what is referred to as a probabilistic
ODE. The discretized version of this stochastic process forms the basis for annealed Langevin
dynamics, and is given by
xt+∆t = xt + γts(xt
; , t)∆t +
p
γt∆t w, 0 ≤ t ≤ T, (4.15)
70

To solve the SDE in Eq. (4.15), it is necessary to approximate the score of the posterior
as a function of time. However, this score inherently relies on the solution of the SDE itself.
Interestingly, the score function for the time-reversed problem is equivalent to the timereversed counterpart of the score function for the forward problem described in Eq. (4.3).
As the forward diffusion equation can be resolved by convolving the posterior with a timedependent Gaussian filter, this provides an effective method for approximating the score
function from data, a technique known as De-noising Score Matching.
4.3.2 Denoising score matching
We aim to approximate the score function of the posterior as a function of time since it
evolves during the forward diffusion (or equivalently the de-noising) process. Importantly,
we can substitute ¯pX(t) with PX|Y (x|y;t) in the annealed Langevin dynamics equation
given in Eq. (4.13), and the process remains valid. We are now interested in finding a neural
network surrogate s(x, y, t) that approximates the score function ∇x log PX|Y (x|y;t) across
the entire support of X, for any given Y and time instance t.
To train the score network s(x, y, t, θ), we minimize the Fisher divergence, which is
equivalent to minimizing the mean squared error (MSE) between true score and the predicted
score. The loss function is given by
L(θ) = Z
DY

1
T
Z T
0
Z
DX

s(x, y, t; θ) − ∇x log PX|Y (x|y;t)

2
PX|Y (x|y;t) dx dt
PY (y) dy
(4.16)
where θ represents the weights of the network.
One can show that the score matching term in Eq. (4.16) can be simplified under mild
regularity conditions [65, 66] and is equivalent to minimizing
L(θ) = Z
DY

1
T
Z T
0
Z
DX
∥s(x, y, t; θ) + ∇x · s(x, y, t; θ)∥
2
PX|Y (x|y;t) dx dt
PY (y) dy
(4.17)
71

However, note that the term ∇x · s(x, y, t; θ) is expensive to compute in high dimensions,
therefore we use an alternate approach proposed by [66], known as the denoising score
matching objective.
In denoising score matching the random variable X ∼ PX|Y (x|y;t) in Eq. (4.16) is
replaced with a perturbed version X˜
σ, which is obtained by adding homoscedastic Gaussian
noise to X ∼ PX|Y (x|y;t = 0) for a given Y . The σ-subscript represents the standard
deviation of the perturbation η.
x˜ = x + η (4.18)
where x˜ ∼ X˜σ and η ∼ N (0, σ2 1). Note that X˜σ and Y are conditionally independent for
given X. Equivalently, the probability distribution for X˜
σ can be obtained by convolving
the original distribution of X with a Gaussian kernel PX˜σ|X(x˜|x, σ(t)):
PX˜σ|Y
(x˜|y, t) = Z
DX
PX˜σ|X(x˜|x, σ(t))PX|Y (x|y) dx (4.19)
where PX˜σ|X(x˜|x, σ(t)) = N (x˜; x, σ(t)
2 1). Note that the use of Gaussian kernel is consistent with the solution to the forward diffusion process in Eq. (4.3) with initial condition
PX|Y (x|y). Thus the denoising score matching objective can be rewritten as
L(θ) = Z
DY
"
1
T
Z T
0
Z
DX˜σ

s(x˜, y, t; θ) − ∇x˜ log PX˜σ|Y
(x˜|y)

2
PX˜σ|Y
(x˜|y) dx˜ dt#
PY (y) dy
(4.20)
We can show, by referring to [67], that Eq. (4.20) can be further simplified to
L(θ) = Z
DY
"
1
T
Z T
0
Z
DX˜σ
Z
DX

s − ∇x˜ log PX˜σ|X

2
PX|Y PX˜σ|X dx dx˜ dt#
PY dy (4.21)
In the previous equation we have simplified function notations in order to make it look
concise.
Next, we rearrange the integral and move the terms PX|Y PY dx dy outside the square
72

bracket. Thus the loss term can be written as having two expectations: the outer one
computed over the joint distribution PXY and the inner one over the conditional PX˜σ|X.
L(θ) = EX,Y ∼PXY
1
T
Z T
0
EX˜σ∼PX˜σ|X

s(x˜, y, t; θ) − ∇x˜ log PX˜σ|X(x˜|x, σ(t))

2
dt
(4.22)
To evaluate the loss function, we need to know the explicit dependence of the integrand
on t. Referring back to Eq. (4.4), the dynamics of the diffusion process is determined by
the dependence of σ on t. To approximate the integral, we use quadrature points in time,
corresponding to L discrete noise levels, denoted as σ1 ≤ σj ≤ σL, where σj = σ(jT/L).
Further, σ(t) is chosen to be an exponential function of t; the rationale for selecting an
exponential function is to concentrate more quadrature points near smaller noise values,
thereby capturing the score function with higher precision in regions where it varies more
rapidly. Additionally, we multiply a weight function λ(t) = σ(t)
2
to the Fisher divergence, as
proposed in [65]. Interestingly, minimizing this weighted combination can be approximated
to maximizing likelihood [68]. Thus, the loss function becomes:
L(θ) ≈ EX,Y ∼PXY "
1
L
X
L
j=1
σ
2
j EX˜σ∼PX˜σ|X

s(x˜, y, σj
; θ) − ∇x˜ log PX˜σ|X(x˜|x, σj )

2
#
(4.23)
Note that we also changed the score function to have σj as a parameter so that we eliminate
the direct dependence of the loss function on t. Additionally, since the kernel function is
Gaussian, we can simplify the estimated score term as follows:
PX˜σ|X(x˜ | x) ∝ exp
−
∥x˜ − x∥
2
2
2σ
2

=⇒ ∇x log PX˜σ|X(x˜ | x) = x − x˜
σ
2
, (4.24)
Substituting this expression into Eq. (4.23), the loss function becomes:
L(θ) ≈ EX,Y ∼PXY "
1
L
X
L
j=1
σ
2
j EX˜σ∼PX˜σ|X

s(x˜, y, σj
; θ) −
x − x˜
σ
2
j

2
#
(4.25)
73

Further, the outer expectation is with respect to the joint distribution of X and Y . We
can easily generate samples from the joint by first sampling X from the prior PX, and
subsequently Y from the likelihood PY |X, which requires solving the forward model F(x).
For a generated x, we also generate L samples of x˜ that corresponding to each noise level.
Thus, the Monte Carlo approximation of the loss function is given by
L(θ) ≈
1
LNtrain
N
Xtrain
i=1
X
L
j=1
σ
2
j

s(x˜
(i)
j
, y
(i)
, σj
; θ) −
x
(i) − x˜
(i)
j
σ
2
j

2
(4.26)
where (x
(i)
, y(i)
) are the i
th realizations of (X,Y ) sampled from the joint PXY . The inner
expectation is approximated using a single realization x˜
(i)
j
from PX˜σ|X(x
(i)
, σj ) that is sampled on the fly during training. This is done to reduce training time and has been observed
to be sufficient in practice [68].
This formulation yields an equation similar to the denoising score matching objective
provided in [68]. The exact dependence of s on σ will be discussed in Section 4.4.2.
4.4 Methodology
In this work, to sample from a complex posterior distribution using diffusion-based methods,
we adopt a technique known as Annealed Langevin Dynamics (ALD). This approach can
be viewed as a discretized, quasi-static version of the reverse diffusion formulation discussed
in Section 4.3.1. It was originally implemented by Song et al. [65] when the stochastic
differential equation (SDE) framework for diffusion models had not yet risen to prominence;
nevertheless, ALD remains a practical choice due to its open-source implementation being
readily available and thoroughly tested for reproducibility.
Although the reverse diffusion SDE provides a strong theoretical backbone, we focus here
on how the annealed version of Langevin dynamics slightly deviates from that continuoustime framework. Instead of solving the time-dependent SDE continuously, ALD applies a
74

quasi-static approximation by using Markov chains to find the stationary distributions across
a discrete set of noise levels {σ1, σ2, ..., σL}. At each noise level σi
, we consider an SDE whose
Fokker–Planck description is given by:
∂p¯X
∂t = −∇ ·
γ(σi)
2
s(x, y, σi) ¯pX

+
γ(σi)
2
∇2
p¯X (4.27)
Further, it can be shown that the SDE derived from this PDE is ergodic with mild regularity
conditions; this means regardless of the initial distribution, the dynamics converge to the
invariant distribution (i.e., the equilibrium) associated with that particular noise setting. Notably, Eq. (4.27) closely resembles the reverse diffusion Eq. (4.13) introduced in Section 4.3.1,
but its dynamics are entirely driven by the ε-terms. Consequently, the steady-state solution
(RHS = 0) is achieved exactly when ∇x log ¯pX = s(x, y, σi). Thus, the invariant distribution is none other than the target distribution PX|Y (x|y) convolved with a Gaussian kernel
whose variance is proportional to σ
2
i
.
We setup a Markov chain to solve the ALD with the following update condition:
xj ← xj−1 +
γ(σi)
2
s(xj−1; y, σi)ϵ +
p
γ(σi)ϵ w, 1 ≤ i ≤ L, (4.28)
where ϵ (or ∆t) is the step size parameter and σi > σi+1. Importantly, xj
from the final
iteration at the noise level σi then becomes the initial condition for the next noise level σi+1.
By progressively decreasing the value of σi
, we “anneal” from a broad initial distribution
to one that closely approximates the true posterior. Thus, ALD effectively de-noises the
samples in stages.
4.4.1 Workflow
This section outlines the overall workflow for solving the conditional probabilistic inference
problem using score-based diffusion models. The complete methodology can be divided into
three main steps: (1) generation of a training dataset, (2) training of a score network, and (3)
75

posterior sampling executed via ALD. The final sampling procedure solves the inverse problem by drawing realizations from the posterior for a given measurement, and subsequently
computes posterior statistics for uncertainty quantification.
Data generation We assume access to a sufficient number of realizations of the random
variable X, which accurately reflects the characteristics of the underlying prior distribution.
These may be previously measured data or samples drawn from a parametric prior distribution designed (based on empirical evidence) to capture the variability in X. Alongside
each realization xi of X, we also draw a corresponding noise sample ηi
from a distribution
representing measurement uncertainty or model discrepancy. Let yi = M(F(xi), ηi) denote
the observed (or predicted) measurement. Note that F(·) can be a black-box solver and need
not be compatible with automatic differentiation. By pairing each xi with its corresponding
yi
, we construct the training dataset S. This dataset is then used for training the score
network.
Training score-network We train the score network on the dataset S by minimizing the
loss function:
LSDM(θ) = 1
LNtrain
N
Xtrain
i=1
X
L
j=1

σjs(x˜
(i)
j
, y
(i)
, σj
; θ) −
x
(i) − x˜
(i)
j
σj

2
(4.29)
which corresponds to the same objective introduced in Eq. (4.26). First, we determine the
hyper-parameters σ1, σL, and L (details on selecting these are provided in Section 4.4.2)
and then train the score network using a suitable optimization algorithm (Adam [69]). As
is standard in machine learning, we use mini-batches of size Nbatch rather than the entire
dataset Ntrain at each iteration, thus estimating LSDM(θ) on a smaller subset to reduce
computational cost.
Additionally, exponential moving averages of the weights are maintained to stabilize
training. Letting θ
(i) denote the parameters after the i
th iteration and θ
′ be an independent
76

variable storing a copy of these parameters, we update θ
′ via
θ
′ = mθ
′ + (1 − m)θ
(i)
, (4.30)
with m = 0.999 to suppress sudden fluctuations in the loss. At the end of training, we
use these averaged parameters as the optimal (final) parameters θ
∗
for the score network.
A sufficiently expressive score network trained using LSDM approximates the score function
across all noise levels. Information regarding the architecture of the score network is given in
Section 4.4.2. Once trained, the score network s(· , · , · ; θ
∗
) can be employed with Annealed
Langevin Dynamics to sample from the target posterior distribution.
Sampling from the posterior In this step, we sample from the target posterior distribution via the Annealed Langevin dynamics shown in Algorithm 1. First, we specify the
number of Langevin steps T and step size parameter ϵ (see Section 4.4.2), and use the trained
score network s(· , · , · ; θ
∗
) to generate new realizations of X for a given measurement yˆ. The
intermediate realizations of Xσ are generated by iteratively running a Markov chain at each
noise level σi
. Unlike in reverse diffusion approach which requires the starting distribution
to be uncorrelated Gaussian), here we may initialize x0 from a uniform distribution, progressively transforming the samples until they closely match the target posterior distribution.
We iterate through the noise levels in descending order, each time running a Markov chain
(i.e., Langevin updates) described by Eq. (4.28). The step size γi
is a function of the current
noise level and is known to be proportional to σ
2
i
for constant time increments based on
Eq. (4.4). We then rescale γi by dividing by σ
2
L
and also absorb the parameter ϵ into it,
yielding the relation
γi = ϵ
σ
2
i
σ
2
L
. (4.31)
In practice, an extra denoising step (Tweedie’s formula [70]) is added to the output of the
last iteration, which has been noted by [71] to improve FID scores [72].
77

Algorithm 1: Annealed Langevin dynamics
Input: Trained score network s(x, y, σ; θ), measurement yˆ, sampling steps T, step
size parameter ϵ, and perturbation kernel standard deviations {σi}
L
i=1
Initialize x0 such that its k
th component {x0}k ∼ U(0, 1) ;
for i = 1 to L do
Set γi = ϵσ2
i
/σ2
L
for j = 1 to T do
Sample z ∼ N (0, 1Nx)
xj = xj−1 + 0.5 γi s(xj−1, yˆ, σi
; θ) + √γiz
end
Set x0 = xT ;
end
Denoise xfinal = x0 + σ
2
L s(x0, yˆ, σL; θ) ;
Output: Realization xfinal
Furthermore, Algorithm 1 is straightforwardly parallelizable, allowing multiple realizations to be generated in tandem. We then use these realizations to estimate the posterior
mean (reconstruction) and posterior standard deviation (uncertainty measure). Moreover,
once the score network is trained, it can be reused for any new measurement y with no
further retraining, making the approach well-suited for high-throughput inverse problems.
4.4.2 Model architecture and hyper-parameter
In this section, we discuss the architecture of the network used to approximate the score
function, highlighting its special design features. We also address hyper-parameter selection,
covering both those parameters related to training the score network and those linked to the
sampling procedure for generating new realizations.
Network Architecture The design of the score networks used in this work follows the
recommendations provided by Song and Ermon [68]. For the implementation, we employ the
NCSNv2 architecture, an advanced version of the U-Net, which adopts a network structure
inspired by RefineNet [73], a model widely used in semantic segmentation tasks. The ELU
activation function is used throughout the network.
A key feature of the architecture is the incorporation of σj as a special input, which
78

normalizes the output of the network’s final layer. This works because the norm of the
score function is empirically noted to scale inversely proportionally to σi
. The inferred field
x and the measurement y are treated as images of 2-dimensional continuum fields that
are interpolated over a Cartesian grid; these are processed as distinct input channels of
the first convolutional layer of the U-Net. For cases involving multiple measurement types,
the number of input channels are extended accordingly. For example, real and imaginary
measurement components are treated as separate input channels along with the noisy x.
Hyper-parameter We follow the recommendations given by Song et al. [68], when choosing hyper-parameters for both training and sampling.
The key parameters for score estimation include the set of noise levels {σ1, . . . , σL} and
the number of noise scales L. The initial noise scale σ1 strongly influences sample diversity;
a large noise value produces a unimodal, nearly Gaussian distribution, which improves convergence when learning from any suitable initial distribution. In addition, regions of X with
sparse training data pose a risk of inaccurate score estimates, causing particles to become
trapped near spurious local modes. Therefore, a gradual decrease in the noise level helps in
the smooth transformation of particles (samples) during ALD sampling. However, an excessively large σ1 increases the total number of noise scales, making the MCMC process more
computationally expensive. In practice, we choose σ1 to be comparable to the maximum
Euclidean distance between all pairs of training data points. On the other end, we set σL
to be very small (e.g., 0.001–0.01), ensuring that the final approximation of the target distribution has minimal added noise. This choice is especially appropriate when the training
data X are normalized to the [0, 1] range.
We choose the number of noise scales L such that the spacing between consecutive noise
levels is sufficiently small for samples at σi to serve as effective initializers for the next
level σi+1. For this, first we set {σ1, σ2, . . . , σL} to form a geometric progression with ratio
β = σi/σi+1. Next, we tune β to ensure an overlap of 0.5 in cumulative density between
79

adjacent noise distributions:
Φ
p
2 Nx (β − 1) + 3 β

− Φ
p
2 Nx (β − 1) − 3 β

≈ 0.5, (4.32)
where Φ is the cumulative density function of the standard normal distribution.
In addition to the noise-related parameters, there are standard training hyper-parameters
for the score network, such as the learning rate, batch size, and number of training iterations.
These are typically chosen to prevent over-fitting, often by monitoring the loss curve on a
validation data set.
Coming to ALD, the associated hyperparameters are the number of Langevin steps T at
each noise level and the step size ϵ. Increasing the number of noise scales L inherently raises
the total number of Langevin steps (and thus the computational cost). Hence, given a fixed
budget, we typically choose T to be as large as resources allow. Finally, we pick ϵ so that
particles at noise level σj can adequately reach their steady-state distribution at noise level
σj+1 within T steps. According to Song et al. [68], one may ensure the following condition
holds:

1 −
ϵ
σ
2
L
2T


β
2 −
2 ϵ
σ
2
L − σ
2
L

1 −
ϵ
σ
2
L
+
2 ϵ
σ
2
L − σ
2
L

1 −
ϵ
σ
2
L



≈ 1. (4.33)
In practice, we perform a grid search over probable values of ϵ and select the one for which
the left-hand side of the above expression is closest to 1.
4.5 Results and Discussion
Here, we demonstrate our methodology for problems that require inferring the geometry
and material properties of mechanical structures ranging from rigid scatters to biological
specimens. While the main focus is on the inverse acoustic scattering problem, we first
80

test and validate our approach on elastography problems. With this, we establish that our
method is both robust and flexible in characterizing the posterior distribution of the inferred
field considering different measurement modalities. We consider three main inverse problems:
Prob1. Quasi-Static Elastography — This is synthetic example where we aim to infer the
Young’s modulus of a specimen containing a stiff inclusion within a softer background. The displacement field under quasi-static loading serves as the measurement, and the prior is relatively simple (e.g., a circular inclusion of fixed radius).
We verify our inference results against full Monte Carlo simulations.
Prob2. Optical Coherence Elastography of Tumor Spheroids — This problem involves inferring the Young’s modulus field within a hydrogel containing tumor spheroids.
The measurements come from optical coherence microscopy (OCM), where phase
differences during compression provide information about the underlying elasticity.
We work with the datasets, both synthetic and experimental, originally compiled
by Foo et al. [74]. The prior distribution used here is geometrically more complex
than the previous elastography example.
Prob3. Inverse Acoustic Scattering — Finally, we infer the topology of rigid elliptical scatterers in a fluid medium by measuring back-scattered pressure signal at sensor
arrays placed some distance away from the scatterers. In addition, we measure the
time harmonic response at multiple frequencies. For this problem, we only consider
synthetic measurements.
For all of these studies, finite element models serve as black-box forward solvers to generate synthetic training data. And the training dataset is min-max normalized to the range
[0, 1] before training. We implement the proposed algorithm using PyTorch [75]. We list
the training and sampling hyper-parameters for the various inverse problems in Table 4.1.
81

Table 4.1: Training and sampling hyper-parameters for the cSDMs
Hyper-parameters Inverse problem
Circular
phantom
(§ 4.5.1.1)
Tumor
spheroids
(§ 4.5.1.2)
Acoustic
scattering
(§ 4.5.2)
Dimensionality Nx 56 × 56 256 × 256 56 × 56
First noise level σ1 20 50 30
Final noise level σL 0.01 0.01 0.01
No. of noise levels L 128 256 160
No. of Langevin steps T 5 5 5
Step size parameter ϵ 5.7e-6 5.7e-6 5.7e-6
Learning rate 0.0001 0.0001 0.0001
Batch size 128 64 128
Total training epochs 100,000 300,000 360,000
4.5.1 Benchmark problems in elastography
In this section, we apply conditional score-based diffusion models (cSDMs) to two benchmark
problems in elastography. These problems were chosen to cover a wide range of scenarios
involving different types of measurements and noise models. Specifically, we consider:
1. Synthetic vs. Real Measurements: Section 4.5.1.1 deals with synthetic measurements
for quasi-static elastography. Section 4.5.1.2 shifts focus to real measurements obtained
from physical experiments. Note that even for experimental measurements, predictions
are made using a score network that is trained on synthetic data drawn from a modeled
joint distribution.
2. Noise and Modeling Errors: In Section 4.5.1.1, the measurement noise is approximated
by homoscedastic zero-mean Gaussian random variables, demonstrating the methodology’s effectiveness under relatively standard assumptions. In contrast, Section 4.5.1.2
addresses non-Gaussian, non-additive noise, confirming that the proposed approach
can accommodate more complex (and realistic) measurement error scenarios.
These examples help validate the proposed approach by comparing the solutions with
Monte Carlo simulation-based inference and real-life experimental results.
82

4.5.1.1 Quasi-static elastography
The first example concerns quasi-static elastography, a medical imaging technique in which
tissues undergoing deformations under an external load are imaged using ultrasound [58,
76]. The ultrasound images are used to determine the displacements of the heterogeneous
specimen, and the goal is to infer the spatially varying shear modulus. This study features a
tissue phantom embedded with a stiffer circular inclusion within a softer surrounding matrix,
such as the realizations shown in the first row of Figure 4.5. In elastography applications,
ultrasound can accurately measure the displacement component along the transducer axis,
which in this case is the vertical direction. Thus, the inverse problem comprises inferring the
spatial distribution of the shear modulus from noisy measurements of vertical displacements
of the phantom.
The underlying mechanics are governed by a linear elastic model described by the following equilibrium and constitutive equations:
∇ · σ = 0 and (4.34)
σ = µ(∇u+∇u
T
) + λ(∇ · u) 1 (4.35)
where σ represents the Cauchy stress tensor for a linear isotropic elastic solid, u denotes
the displacement field, and µ denotes the shear modulus. Moreover, we assume that the
specimen is incompressible and in plane stress condition which means λ = 2µ.
In this example, the dimensions of the specimen is 1×1 cm2
. To simulate compression,
we fix the bottom edge and impose a downwards vertical displacement of 0.01 cm on the top
edge. The bottom edge is traction-free in the horizontal direction, whereas the left and right
edges are kept traction-free in both directions. Additionally, we pin the bottom left corner to
prevent rigid body motion. A realization of X corresponds to the spatial distribution map of
the shear modulus of the specimen discretized over a 56×56 Cartesian grid. We formulate a
parametric prior that represents a uniformly stiff circular inclusion set against a homogeneous
83

0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
x
y
Figure 4.5: Five typical realizations of X and Y sampled from the joint distribution for
the synthetic quasi-static elastography example and belong to the training dataset. The
first row shows the shear modulus field, while the second row shows the corresponding noisy
vertical displacement field used as measurements. In this figure, the standard deviation of
the Gaussian noise ση = 0.025. Moreover, realizations of X are Y are min-max normalized
to [0, 1] scale
Table 4.2: Random variables comprising the parametric prior for the synthetic quasi-static
elastography problem. Note, U(a, b) denotes an uniform random variable between a and b
Random variable Distribution
Distance of the inclusion’s center from the left edge (mm) U(0.2, 0.8)
Distance of the inclusion’s center from the bottom edge (mm) U(0.2, 0.8)
background: the coordinates of the inclusion’s center are uniformly distributed over the
domain such that the inclusion fits within the boundaries. We intentionally work with a
simple prior, consisting of only two random variables, to facilitate Monte Carlo simulations
(MCS). See Table 4.2 for details regarding the prior distribution for this example. We fix
the inclusion’s radius to be 0.12 cm and choose the shear modulus of the inclusion and the
background to be 1.5 and 0.1 kPa, respectively.
The training dataset comprises 10,000 pairwise realizations of X and Y . For a given
realization x of X, we generate the corresponding measurement y, which is a realization of
Y , by utilizing x in a finite element model consisting of linear triangular elements. We use
FEniCS [77] to conduct the finite element analysis. We discretize the measurement over the
same 56×56 Cartesian grid as the shear modulus field. Therefore, nX = nY = 56×56 in this
84

example. Subsequently, we add measurement noise — homoscedastic Gaussian noise with
standard deviation ση ×umax truncated to within the [−3ση ×umax, 3ση ×umax] range, where
umax is the maximum vertical displacement across all specimens in the training dataset.
Figure 4.5 illustrates five randomly selected realizations from the training dataset. Note
that in the training data, both X and Y , are min-max normalized to the range [0, 1]. We
work with three values of ση ×umax corresponding to varying intensity of measurement noise
— 2.5%, 5%, and 10% of umax — to create three training data sets, which we then use to train
the score network. Note, we train a separate score network for each level of measurement
noise. Table 4.1 lists the training and sampling hyper-parameters we use for this problem.
After training the score network, we employ Annealed Langevin dynamics, as described in
Section 4.4, to sample the posterior distribution.
iteration
iteration
iteration iteration iteration
iteration iteration iteration Figure 4.6: Visualization of the training process using a single realization sampled using
ALD.
An interesting pattern is observed during the cSDM training; the model appears to
learn the underlying circular prior relatively quickly, by around the 10000th iteration, before
gradually refining the placement of the circles to match the full conditional distribution
85

over the next few tens of thousands of iterations. This two-phase learning behavior has
been consistently observed across all our examples, suggesting that early training focuses on
capturing the broad structure of the prior, while subsequent epochs refine how measurement
data Y affects the distribution of X, ultimately converging to a more accurate conditional
model.
We also estimate the posterior statistics using Monte Carlo simulation (MCS), which
we treat as the true statistics, to validate the results obtained using the cSDMs. We use
a sample size equal to 500,000 for the Monte Carlo simulations, which is significantly more
than the size of the training dataset. The simulation procedure is outlined as follows:
1. Compute the likelihood probability of the measurement y˜ ∈ Y given x
(i) ∀ ∈ X, where
1 ≤ i ≤ NMCS.
PY |X(y˜|x
(i)
) = 1
q
(2π)
kσ
2k
η
exp
−
1
2σ
2
η
∥y˜ − F(x
(i)
; 0)∥
2

, k = 562
(4.36)
2. Compute the posterior probability ∀ x
(i) ∈ X given y˜.
PX|Y (x
(i)
|y˜) = PY |X(y˜|x
(i)
)
PNMCS
i=1 PY |X(y˜|x(i)
)
(4.37)
Note that the prior probability density is assumed to be constant for this experiment.
3. Compute the expected value of f(X) by taking the weighted average over all x
(i)
.
E[f(X)] =
N
XMCS
i=1
PX|Y (x
(i)
|y˜)f(x
(i)
) (4.38)
Note that the measurement noise is additive in this case and the likelihood function
can be used to obtain the importance weights, which we subsequently self-normalize. We
compare the posterior statistics estimated using cSDM and MCS on two test samples that
were not part of the training dataset; Figure 4.7 shows the results. In Figure 4.7(a) and (b),
86

column 1 shows the test phantom while column 2 shows the corresponding measurement for
three different noise levels. Note that the ground truth is not readily discernible through
visual inspection of the measurements in column 2 of Figure 4.7.
Figure 4.7 also shows the pixel-wise posterior mean and standard deviation estimated
using the cSDM and Monte Carlo simulations. We generate 10000, 4000, and 1000 posterior
realizations to compute these statistics using the trained score network for ση =10%, 5%, and
2.5%, respectively. These sample sizes were roughly close to the effective sample size of the
Monte Carlo samples for the three cases. Figure 4.7 shows that the cSDM mean successfully
locates the inclusion at all measurement noise levels, similar to MCS. Furthermore, the posterior statistics estimated using cSDM are qualitatively similar to those obtained using MCS,
with the root mean squared errors (RMSE) listed in Table Table 4.3. Both these methods
predict larger uncertainty around the edges of the inclusion, and the uncertainty predictably
grows with the level of measurement noise. For the lowest level of measurement noise, the
reconstructed inclusion in the posterior mean is sharp as there is very little uncertainty,
both in spread and magnitude, surrounding the inclusion. The cSDM yields results similar
to Monte Carlo inference for other test phantoms as well. Table 4.4 tabulates the RMSE
between the posterior statistics estimated using cSDM and MCS across ten test phantoms.
The consistently low average RMSE values suggest that the proposed method is accurate
and robust to the amount of measurement noise.
Table 4.3: Root mean squared error (RMSE) between the pixel-wise posterior statistics
estimated using cSDM and MCS, for the test phantoms in Figure 4.7
RMSE in posterior statistics
Measurement
noise (ση)
Mean Std. Dev.
Test sample 1 Test sample 2 Test sample 1 Test sample 2
10.0% 0.036 0.027 0.038 0.029
5.0% 0.039 0.025 0.032 0.029
2.5% 0.031 0.022 0.028 0.027
87

0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
Ground truth Measurement Mean Std. dev.
Posterior statistics (cSDM)
Mean Std. dev.
Posterior statistics (MCS)
σ
η = 0.1
σ
η = 0.05
σ
=
η
0.025
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
Ground truth Measurement Mean Std. dev.
Posterior statistics (cSDM)
Mean Std. dev.
Posterior statistics (MCS)
σ
η = 0.1
σ
η = 0.05
σ
η = 0.025
(a) Test sample 1
(b) Test sample 2
Figure 4.7: Posterior statistics – pixel-wise posterior mean and standard deviation – estimated using the conditional diffusion models (cSDMs) and Monte Carlo simulation (MCS)
for the synthetic quasi-static elastography example on two test samples. In subfigures
(a) and (b), the three rows show the results corresponding to three different levels of measurement noise ση but for the same ground truth. All values are normalized to [0, 1] scale
88

Table 4.4: Root mean squared error (RMSE) between the pixel-wise posterior statistics
estimated using cSDM and MCS, averaged over ten test phantoms for the synthetic quasistatic elastography example
Average RMSE in posterior statistics
Measurement
noise (ση)
Mean Std. Dev.
10.0% 0.051 0.051
5.0% 0.031 0.032
2.5% 0.025 0.029
Effect of model misspecfication To understand the effects of a misspecified model for
PY |X, we perform inference for the case when the measurement noise’s standard deviation
ση = 0.025 and 0.1 with the score network trained assuming ση = 0.05. For test sample 1,
Figure 4.8 shows the posterior statistics estimated using the cSDM and compares it against
the MCS under the misspecified noise model. Comparing Tables 4.3 and 4.5, the RMSE in the
posterior statistics estimated using the cSDM increases when the noise model is misspecified.
In particular, Figure 4.8 and Table 4.5 show that the pixel-wise posterior mean is particularly
biased when the noise in the measurement has standard deviation ση = 0.1 but inference is
performed assuming ση = 0.05. These results highlight the need for studying the effects of
model misspecification on the performance of cSDMs.
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
Ground truth Measurement Mean Std. dev.
Posterior statistics (cSDM)
Mean Std. dev.
Posterior statistics (MCS)
σ
η = 0.1
σ
η = 0.025
Figure 4.8: Posterior statistics – pixel-wise posterior mean and standard deviation – for
test sample 1 estimated using cSDMs trained assuming measurement noise with standard
deviation ση = 0.05. The reference MC statistics are also estimated with the misspecified
measurement noise model. All values are normalized to [0, 1] scale
89

Table 4.5: Root mean squared error (RMSE) between the pixel-wise posterior statistics
estimated using cSDM and MCS for misspecified values of ση
Measurement
noise (ση)
Mean Std. Dev.
10.0% 0.069 0.038
2.5% 0.024 0.037
4.5.1.2 Optical coherence elastography of tumor spheroids
Optical coherence tomography (OCT) is a technique that leverages optics, analogous to
acoustics in ultrasound imaging, to produce high-resolution tomograms of tissue. A functional extension of OCT that enables imaging based on mechanical contrast is optical coherence elastography (OCE) [78], which involves capturing the phase data alongside intensity
measurements of the interference of light back-scattered from the specimen with light reflected from a reference mirror at varied axial depths. One variant of OCE is phase-sensitive
compression OCE [79], wherein the phase shifts observed between B-scans (cross-sectional
images) of the specimen as it undergoes compressive loading quantify the specimen’s deformation along the optical axis. The mechanical properties of the specimen can be inferred
from these phase shifts because they encode within them information about the displacement component along the optical axis. In this study, we consider a high resolution variant
of phase-sensitive compression OCE, termed mechano-microscopy [80, 81], which utilizes a
high-resolution variant of OCT known as optical coherence microscopy (OCM). A conditional
score-based diffusion model is used to solve the inverse problem of determining elasticity profiles, specifically the spatial distribution of the Young’s modulus X of tumor spheroids from
noisy phase difference measurements Y . The data for this experiment were obtained from
measurements and datasets originally compiled by [74].
We provide a brief summary of the experimental setup here (Figure 4.9) and the reader
is referred to [74] for more details. In the experiment, a custom-built spectral-domain OCM
system in dual-arm configuration images tumor spheroids embedded in hydrogel specimens
90

Objective
Coverslip
Silicone
Specimen
Actuator
To OCM
Figure 4.9: Schematic showing the experimental setup in the tumor spheroid application
from the top. Following the objective lens, the sample arm contains, from top to bottom: a
rigid glass coverslip acting as an imaging window, a 300 µm thick sample containing tumor
spheroids, a 1 mm thick compliant silicone layer, and a piezoelectric actuator to apply
compressive displacement. The purpose of the compliant layer is to provide surface stress
measurements for mechano-microscopy [80]; we do not utilize this data in this study. Two
volumetric images at different compression levels were acquired at each location to produce
a volumetric phase difference image.
The plane corresponding to the central cross-section of the tumor spheroid was identified
in the volumetric image and was cropped to a 256×256 µm2 window around the spheroid.
This two-dimensional phase difference image was used as measurement in solving the inverse
problem.
The specimens were prepared by encapsulating MCF7 (non-metastatic breast cancer)
cells within a hydrogel matrix composed of lyophilized gelatin methacryloyl (GelMA). The
cells were cultured until they formed spheroids. Each tumor spheroid contains multiple
cells, and each cell domain can be broadly classified into two components: the cytoplasm
and the nucleus (see Figure 4.10). We take inspiration from this morphology to develop
the parametric prior for this study. We refer interested readers to [82] for typical threedimensional immuno-fluorescence images of tumor spheroids taken at different time points.
The synthetic samples, which we generate from the parametric prior distribution, repre91

sent a spheroid within a hydrogel on a 256×256 µm2 domain (first column in Figure 4.11).
We model the spheroid as an ellipse, and straight lines subdivide the spheroid into cells.
Each cell contains a circle representing a nucleus, surrounded by the cytoplasm of the cell.
We treat the following parameters as truncated normal random variables: spheroid width,
spheroid aspect ratio, spheroid position, number of nuclei per unit area, nuclei radii, nuclei
elasticity, cytoplasm elasticity, and hydrogel elasticity. Table 4.6 provides details on these
parameter’s distributions, which we derive from preliminary analyses of spheroid OCM images. We determine the nuclei locations (defined by their center coordinates) using a custom
algorithm that iteratively accumulates nuclei coordinates by creating a list of valid coordinates for the current nucleus to be positioned, randomly selecting a coordinate from this
list for the current nucleus, removing the coordinates occupied by the current nucleus from
the list of valid coordinates for the next nucleus, and repeating until the required number of
nuclei for the spheroid is reached or there are no more valid coordinates for the next nucleus.
We constrain the sizes and positions of the nuclei to ensure a minimum spacing of 2 µm
between nuclei and between nuclei and the ellipse boundary. We subsequently determine the
cell boundaries using a power diagram, also known as a Laguerre–Voronoi diagram [83].
Table 4.6: Random variables comprising the parametric prior for the tumor spheroid application
Parameter Variable Minimum Maximum
Spheroid width (µm) N (105, 8.33) 80 130
Spheroid height to width ratio N (0.90, 0.0167) 0.85 0.95
Spheroid (center) horizontal position (µm) N (0, 13.33) -40 40
Spheroid (center) vertical position (µm) N (0, 13.33) -40 40
Number of nuclei per square area (mm−2
) N (4750, 250) 4000 5500
Radius of nuclei (µm) N (7.5, 0.5) 6 9
Young’s modulus of nuclei (kPa) N (32.5, 5.833) 15 50
Young’s modulus of cytoplasm (kPa) N (10, 1.67) 5 15
Young’s modulus of hydrogel (kPa) N (7.5, 0.833) 5 10
Displacement applied to top edge (µm) N (−2.2, 0.6) −0.4 −4
*In this example, the normal distributions were truncated at three standard deviations from the
mean, with truncation limits explicitly defined by the specified minimum and maximum values.
92

Figure 4.10 shows the steps for generating realizations of Y from realizations of X.
Given a realization x of the Young’s modulus field, we simulate uniaxial compression on
the specimen via finite element analysis (FEA) in ABAQUS [84] to obtain the vertical
displacement field (third column in Figure 4.10). To mitigate edge effects, we initially expand
each sample to an 800×800 µm2 domain by padding with hydrogel on all four sides. We
displace the top edge downward while the bottom-left corner remains fixed, and the bottom
edge is constrained vertically. We sample the displacement on the top edge from a truncated
normal distribution; see Table 4.6. The magnitude of the displacement of the top edge ensures
that the corresponding strain values are within the order of a few millistrains, consistent with
observations from mechano-microscopy experiments. For this analysis, given the low levels
of strain involved, we assume that the specimen is linear and elastic (Eq. (4.35)), under plane
strain equilibrium conditions (Eq. (4.34)), and nearly incompressible with Poisson’s ratio ν
uniformly set to 0.49. Under plain strain, the Lam´e parameters in Eq. (4.35) are expressed
in terms of Young’s modulus and Poisson’s ratio using the relations
λ(ψ) = νE(ψ)
(1 + ν)(1 − 2ν)
, µ(ψ) = E(ψ)
2(1 + ν)
. (4.39)
We develop a 2D CAD model for the specimen (second column in Figure 4.10), which we
subsequently use to develop the FEA mesh geometry consisting of quadrilateral elements
for the main region and triangular elements for the padded regions. Once the FEA is
complete, we obtain 2D arrays of the Young’s modulus E and vertical displacement u2 by
interpolating the main 256×256 µm2 domain window to a 256×256 Cartesian grid. We
convert the displacement fields to phase difference fields ∆ϕ using the following equation:
∆ϕ =
4πnru2
λ0
, (4.40)
where nr = 1.4 is the refractive index of the hydrogel, and λ0 = 800 nm is the central
wavelength of the OCM light source.
93

Prior realization x CAD model Vertical displacement Noisy phase difference y
FEA
Nuclei
Cytoplasm
Hydrogel
Figure 4.10: Workflow for generating the training data in the tumor spheroid application
Subsequently, we add non-Gaussian measurement noise ϕη to the phase difference field to
obtain a noisy measurement image (fourth column in Figure 4.10). The noise in the phase
difference varies with depth and is modeled to follow a probability law characterized by the
following density function [74]:
pΦη
(ϕη) = 1
2π
exp
−
k
2
2

+
r
1
8π
k exp
−
k
2
2
sin2 ϕη

cos ϕη

1 + erf
k cos ϕη
√
2
, (4.41)
where k is related to the variance of this distribution. The variance varies with depth
along the optical axis in a B-scan due to optical attenuation and the confocal function of
the objective lens, and was determined as a function of the depth by curve fitting based
on experimental calibration conducted with hydrogel specimens [74]. After the addition of
noise, the phase difference is wrapped to ensure it remains within the range (−π, π]. We
obtain the final wrapped phase difference image ∆ϕw using the following formula
∆ϕw = arg (exp [i(∆ϕ + ϕη)]). (4.42)
This wrapped and noisy phase difference field is the realization y of the measurements
corresponding to the realization x of the Young’s modulus field.
In this manner, we initially synthesize a synthetic dataset comprising 5000 pairs of X
and Y realizations. We augment this dataset sixfold by shifting the 256×256 µm2
image
window 20 µm in each horizontal and vertical direction (left, right, up, and down), as well as
94

0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
x
y
Figure 4.11: Five typical realizations of X and Y sampled from the joint distribution that
form the training dataset for the tumor spheroid experiment. In the first row, we have
the log-normalized Young’s module fields, and in the second row, the corresponding noisy
measurements are displayed. All values are normalized to [0, 1] scale
applying a horizontal flip to each image. We divide the resulting dataset into a training set of
24000 images and a test set of 6000 images. As a preprocessing step, we rescale the elasticity
images to enhance contrast and facilitate better training convergence. The rescaling involves
the following transformation
E
′ = log10
E
Ehydrogel
, (4.43)
where E
′
is the rescaled Young’s modulus. This transformation means that it will only be
possible to recover the Young’s modulus field up to the normalization constant Ehydrogel,
which is the Young’s modulus of the hydrogel. Figure 4.11 illustrates five randomly selected
samples from the dataset. We utilize the training data to train the score network. After
training the score network, we employ Annealed Langevin dynamics to sample elasticity
images from the posterior. See Table 4.1 for additional information regarding the training
and sampling hyperparameters.
Figure 4.12 shows realizations from the posterior, and its statistics, estimated using the
conditional score-based diffusion model on both synthetic and experimental data. The first
two columns correspond to samples from the test dataset previously unseen by the model,
while columns three and four depict predictions on experimental measurements. The first
95

row shows the measurements used for sampling using the trained score network. Rows two
to four show three realizations sampled by the model for the corresponding measurements.
Rows five and six show the pixel-wise posterior mean and standard deviation, respectively,
estimated using 200 realizations. For the synthetic samples, we compare the estimated
posterior means for the Young’s modulus distribution with the corresponding ground truth.
In the case of experimental data, comparisons were made against OCM images that were
co-registered during the acquisition of the phase difference images. Recall that the recovery
of the Young’s modulus is possible up to a multiplicative constant: the hydrogel modulus
Ehydrogel in this case. We assume Ehydrogel = 10 kPa for the test samples. For the experimental
samples, Ehydrogel was estimated from mechano-microscopy measurements [74].
From the results for the synthetic data, we note that each realization from the posterior
distribution correctly captures the volume and the location of the spheroid. Further, there are
significant differences between each realization in the modulus distribution predicted within
each spheroid. We can attribute these differences to the ill-posed nature of the inverse
problem [57, 85], which stipulates that there are multiple modulus distributions that are
consistent with a given noisy measured displacement field. The differences in the modulus
field within the spheroid lead to a fuzzy mean modulus field. Remarkably, when we compare
this field to the ground truth, we observe that the location of several nuclei, which correspond
to regions with elevated stiffness, is captured accurately. The uncertainty in the modulus
field within the spheroid directly translates to large values of standard deviation within the
spheroid. We also observe that regions of the largest standard deviation (uncertainty) are
associated with the boundaries of the nuclei.
From the results for the experimental data, we observe that the measured phase field
is more heterogeneous and noisy than its synthetic counterpart. Thus, we anticipate more
variability and uncertainty in the predictions, which is precisely what we observe when
examining the three realizations from the posterior distribution. This translates to a fuzzy
image for the mean stiffness field and large standard deviation values (when compared with
96

−π
0
π
rad
1
10
100
kPa
1
10
100
kPa
1
10
100
kPa
1
10
100
kPa
1
10
100
kPa
1
10
100
kPa
−π
0
π
rad
1
10
100
kPa
1
10
100
kPa
1
10
100
kPa
1
10
100
kPa
1
10
100
kPa
1
10
100
kPa
−π
0
π
rad
1
10
100
kPa
1
10
100
kPa
1
10
100
kPa
1
10
100
kPa
1
10
100
kPa
0
10
20
30
40
dB
−π
0
π
rad
1
10
100
kPa
1
10
100
kPa
1
10
100
kPa
1
10
100
kPa
1
10
100
kPa
0
10
20
30
40
dB
Test sample1 Test sample2 Sample1 Sample2
Synthetic data Experimental data
Measurement Instance 1 Instance 2 Instance 3 Mean Std. dev. GT/OCM
Realizations Posterior statistics Reference
Figure 4.12: Posterior statistics estimated for selected samples using the cSDM for the tumor
spheroid problem. The field of view in each image is 256×256 µm2
. Note that all images are
presented at their true scales. GT: Ground truth.
97

synthetic data) within the spheroid. The underlying ground truth image is unavailable in the
experimental case. Therefore, we cannot quantitatively evaluate the accuracy of the modulus
field. However, from the OCM images (last row of Figure 4.12), we can infer the location
and the size of the spheroid, and when we compare this with the mean modulus image we
conclude that the location and the size of the spheroid have been captured accurately.
We conclude this section by noting that although the initial tumor spheroid results are
encouraging, they point to several extensions that are currently underway. These include
improving the forward model for measurement noise so that the synthetic measurements
resemble their experimental counterparts more closely and performing inference on tissuemimicking phantoms of tumor spheroids for a more quantitative comparison.
4.5.2 Inverse acoustic scattering
Inverse acoustic scattering seeks to recover information about unknown scatterers (e.g., their
location, shape, or physical properties) by probing a medium with incident acoustic waves
and measuring the resulting scattered fields. This problem arises in numerous applications,
from biomedical imaging [86] to geophysical exploration [87], where one aims to infer subsurface features or hidden objects in a non-invasive manner.
Here, our focus is on inferring the two-dimensional topology of rigid scatterers within a
given domain, based on measurements of the scattered pressure field recorded by an array of
sensors placed around the domain boundary. We consider both single-frequency and multifrequency time-harmonic pressure waves as incident signals. Additionally, we compare the
performance of our reconstruction (topology inference) approach in two sensor-placement
scenarios: sensors arranged only in one direction versus sensors distributed all around the
domain boundary.
Although inverse scattering has been studied extensively, achieving accurate and stable
solutions remains challenging due to its highly ill-posed nature and the oscillatory behavior of
the underlying physics. Even obtaining reliable forward simulations and measurements can
98

be difficult. Recent machine learning (ML) based methods [52, 53, 88] have made progress in
solving inverse problems within acoustics, yet uncertainty quantification remains relatively
under-explored. In this work, we propose a framework to estimate the topology of rigid
scatterers while quantifying uncertainty in the reconstructed domain. Although our work
focuses on rigid structures for validation and demonstration purposes, the methodology can
be readily applied to more general structural-acoustic systems with deformable scatterers.
Gaussian
Noise
Real Imaginary
PML
Scatterer +
Frequency
Case 1
Single frequency Multiple frequency
Case 2 Case 3
Figure 4.13: Schematic of the data generation process. The top-middle diagram illustrates
the analysis domain with its surrounding PML region. At top-left is an example realization
of X, and at top-right is the corresponding scattered pressure field (magnitude) obtained
via FEA. The bottom plots display the boundary measurements, color-coded according to
the three different measurement cases.
As an example problem, let us consider a square domain of size 2 × 2 m2
(the blue region
in Figure 4.13), containing elliptical scatterers located sufficiently far from the boundaries
99

so as not to intersect them. A time-harmonic propagating plane wave,
pi = e
−ikx
, (4.44)
is incident from the right boundary, where k denotes the wavenumber. We denote the
scattered wavefield by ps = p − pi
, with p being the total field. Physically, the problem
is modeled as a constant-density time-harmonic acoustic wave in the frequency domain,
governed by the Helmholtz equation within the domain Ω.
On the rigid scatter boundary ΓI , we impose the homogeneous Neumann condition.
Outside Ω, we introduce a perfectly matched layer (PML) [59] of width 0.5 m (the red region
in Figure 4.13), denoted by ΩPML, to approximate a radiating (non-reflective) boundary
condition. Consequently, the scattered field ps satisfies
∇2
ps + k
2
ps = 0 on Ω, (4.45)
∇¯ 2
ps + k
2
ps = 0 on ΩPML, (4.46)
∇ps = −∇pi on ΓI , (4.47)
∇ps = 0 on Γ, (4.48)
where Γ is the outer boundary of ΩPML. The modified gradient ∇¯ in the PML region is
defined by
∇¯ =
1
1 + i
σx(x)
k
∂
∂x +
1
1 + i
σy(y)
k
∂
∂y (4.49)
where σx and σy are the damping profiles in the PML, designed to attenuate outgoing waves.
Specifically,
σx(x) =



c0(xmin − x)
2 x ≤ xmin
0 xmin ≤ x ≤ xmax
c0(x − xmax)
2 xmax ≤ x
(4.50)
100

σy(y) =



c0(ymin − y)
2
y ≤ ymin
0 ymin ≤ y ≤ ymax
c0(y − ymax)
2
ymax ≤ y
(4.51)
with {xmin, xmax, ymin, ymax} denoting the bounds of Ω and c0 = 24 m−1
. This construction
ensures that the waves in ΩPML are progressively damped, preventing reflections at the outer
boundary.
In this work, we assume a parametric prior for the scatterer geometry consisting of
one or two ellipses placed randomly within a square region centered inside Ω. The center
coordinates, aspect ratio, size, and orientation of each ellipse are treated as random variables,
all drawn from uniform distributions as specified in Table 4.7. Additionally, we let the
number of ellipses follow a binary random variable equal to 1 or 2 with a probability ratio
of 1:3, thereby capturing both single-ellipse and double-ellipse configurations, which may be
either intersecting or separate. To avoid meshing challenges at contact points or extremely
narrow gaps, we resample from the prior if two ellipses touch at the boundary in a way that
provides neither sufficient overlap nor sufficient distance. Each resulting scatterer topology
is represented by a binary field, 1 within rigid scatterers and 0 elsewhere, which is mapped
onto a 56 × 56 Cartesian grid to form the realizations of X. We employ a conforming mesh
refined near scatterer boundaries, ensuring accurate geometry representation and numerical
stability.
Once we generate the scatterer geometries, we solve the forward problem Eqs. (4.45)
to (4.48) using FEA to obtain the scattered pressure field throughout the domain. As
an example, the amplitude of the total pressure field for one such topology is shown in
Figure 4.13 (top right). We assume the incident pressure wave has unit magnitude, since the
inference task here is insensitive to absolute signal amplitude. The wavenumber k depends
101

Table 4.7: Random variables comprising the parametric prior for the acoustic scattering
Parameter Variable
X-coordinate (center) of ellipse-1 (m) U(−0.65, 0.65)
Y-coordinate (center) of ellipse-1 (m) U(−0.65, 0.65)
Aspect ratio of ellipse-1 U(0.5, 1.0)
Major axis of ellipse-1 (m) U(0.4, 0.6)
Orientation of ellipse-1 (rad) U(−1.57, 1.57)
X-coordinate (center) of ellipse-2 (m) U(−0.65, 0.65)
Y-coordinate (center) of ellipse-2 (m) U(−0.65, 0.65)
Aspect ratio of ellipse-2 U(0.5, 1.0)
Major axis of ellipse-2 (m) U(0.4, 0.6)
Orientation of ellipse-2 (rad) U(−1.57, 1.57)
on the sound frequency ω, the bulk modulus κ, and the density ρ of the acoustic medium:
k =
r
ρ ω2
κ
. (4.52)
Here, the acoustic medium is water with ρ = 997 kg/m3 and κ = 2.15 GPa. We use 56
frequencies uniformly spaced between 16 − 45 krad/s and solve the time-harmonic acoustic
scattered field at each frequency via FEA. For measurements, we place an array of 56 sensors
along the boundary of Ω on each side and record the time-harmonic back-scattered pressure
field. These recorded data are then organized into a 56 × 56 grid, forming a realization of Y
corresponding to the scatterer configuration sampled from X. We consider three different
measurement scenarios:
Case 1: Single-frequency (one boundary). We record only at the highest frequency, where
the 56 sensor readings form a single 56 × 1 array repeated across all rows of the
56 × 56 grid. Two channels capture the real and imaginary parts of the scattered
field (Y0,Y1 in Figure 4.13).
Case 2: Multi-frequency (one boundary). We repeat the procedure of Case 1 but now each
row of the 56 × 56 grid corresponds to one of the 56 frequencies and is gathered
from the same boundary, increasing from top to bottom (Y2,Y3 in Figure 4.13).
102

Case 3: Multi-frequency (four boundaries). Here, measurements are taken around all four
edges of Ω and at multiple frequencies. Consequently, Y includes eight channels
{Y0,Y1, . . . ,Y9}, consisting of two channels (real and imaginary) per boundary side.
To emulate real-world measurement conditions, we add Gaussian noise (truncated to [−3σ, 3σ]
range) with zero mean and standard deviation equal to 3% of the maximum range of the
recorded scattered field values for all frequencies.
0.0
0.5
1.0
0.26
0.48
0.70
0.29
0.50
0.71
0.31
0.51
0.72
0.28
0.50
0.71
0.0
0.5
1.0
0.27
0.48
0.68
0.29
0.50
0.70
0.27
0.48
0.69
0.29
0.48
0.67
0.0
0.5
1.0
0.36
0.49
0.62
0.36
0.49
0.62
0.37
0.50
0.62
0.39
0.51
0.63
0.0
0.5
1.0
0.20
0.48
0.76
0.19
0.47
0.76
0.29
0.49
0.69
0.29
0.55
0.81
0.0
0.5
1.0
0.25
0.50
0.75
0.24
0.49
0.74
0.22
0.47
0.72
0.25
0.50
0.75
Sample 1 Sample 2 Sample 3 Sample 4 Sample 5
X
Y
2
Y
3
Y
0
Y
1
Figure 4.14: Five typical realizations of X and Y sampled from the joint distribution for
the acoustic scattering example. The top row shows the binary scatterer topology field X,
and the subsequent rows present the noisy back-scattered pressure signals Y . Here, Y0, Y1
correspond to the real and imaginary components of the single-frequency case, and Y2, Y3
represent the multi-frequency case. All data are min–max normalized to [0, 1] scale.
103

We generated 15,000 pairs of X and Y , where X has one channel and Y has ten channels.
To expand this dataset without incurring additional forward simulation costs, we flipped X
horizontally around the mid-axis and applied the same flip operation to the corresponding
measurement arrays in Y . The resulting augmented dataset was split into three groups, each
intended for a specific measurement case. In the first group, there are three channels overall
(counting both X and Y ), the second group also has three channels, and the third group
has nine channels. Each group contains 30,000 samples, which we partitioned into 25,000
for training and 5,000 for testing. Few examples from the training dataset are provided in
Figure 4.14.
We used the same general network architecture for all three groups, with adjustments
only to the initial input layer so that it matches the respective channel count. Table 4.1
provides details on the training hyper-parameters, including batch size and number of epochs.
After training the score network for each group, we sampled from the resulting posterior
distributions using Langevin dynamics (parameters also detailed in Table 4.1) to compute
the posterior statistics for the three cases.
Figure 4.15 shows the predictions obtained from the trained score network for four randomly chosen test samples. The top row shows the ground truth configurations used to
generate the corresponding measurements. The subsequent rows display the posterior mean
and standard deviation computed from 1,000 realizations for each of the three measurement
cases using three separately trained score networks.
As we progress from Case 1 to Case 3, increasing the amount of measurement information
available to the network, the mean predictions match more closely with the ground truth, and
the posterior uncertainty consistently decreases. Notably, the network accurately identifies
the scatterer locations across a variety of topologies, including single, double, and intersecting
scatterers. A closer look at the standard deviation plots reveals that the uncertainty is
concentrated near the scatterer boundaries. This reflects the intrinsic difficulty in capturing
the exact interface profile, especially under sparse or directional measurements. We also
104

0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
Test sample 1 Test sample 2 Test sample 3 Test sample 4
Ground truth Mean Std. dev. Mean Std. dev. Mean Std. dev.
Case 1 Case 2 Case 3
Figure 4.15: Posterior statistics estimated for case 1, case 2, and case 3 using measurements
corresponding to the ground truth signal (top row), for four representative test samples.
105

observe that when a scatterer is closer to the sensor array, the interface is reconstructed
with higher accuracy and lower uncertainty. This is especially evident in Case 3, which has
sensor coverage on all four sides. Conversely, in scenarios where the scatterer is situated
away from the sensor’s field of view (for example, see Test sample 4, Case 1), the single
frequency back-scattered data yields a lower accuracy interface reconstruction compared to
samples where the scatterer is located near the mid horizontal axis.
Based on our results, we observe that the network effectively captures the underlying
posterior distribution consisting of different topologies for all measurement scenarios. Moreover, when a scatterer is closer to a sensor array, the interface appears sharper with lower
standard deviation, highlighting how proximity to measurement locations improves the fidelity of local reconstruction. These observations align with our intuition that more diverse
or comprehensive data (multiple frequencies, additional sensors) yield both higher accuracy
and lower uncertainty in inverse scattering problems.
This analysis demonstrates how our technique can be applied to a range of acoustic sensing tasks, particularly in structural acoustic settings where topology estimation is crucial.
In practice, predictive performance will depend on factors such as the number of frequencies
analyzed, the complexity of the prior, and the placement of sensors. Optimizing sensor arrangements to maximize reconstruction accuracy thus emerges as a promising direction for
future work. Subsequent research could target more complex topologies and incorporate material property distributions along with geometry, further broadening the scope and impact
of this technique.
Parts of this chapter are taken from the reference:
Agnimitra Dasgupta, Harisankar Ramaswamy, et al. Conditional score-based diffusion models for solving
inverse elasticity problems. In: Computer Methods in Applied Mechanics and Engineering 433 (2025), p.
117425.
106

Chapter 5
Conclusion
5.1 Summary
The main purpose of this thesis has been to explore and develop efficient computational
frameworks for solving a variety of coupled structural–acoustic problems, which we have
broadly classified under optimization and inverse problems. Our efforts include material,
shape, and topology estimation in both deterministic and probabilistic settings, ultimately
leading to large-scale design and predictive tools.
In Chapter 2 we have presented a first-principles-based formulation, and associated computational procedures, for the optimization of coupled linear elasto-acoustics systems applicable to both interior and exterior problems in structural-acoustics. While we specifically
address optimizing the material parameters of the elastic structure, the formulation may be
extended in a straightforward manner to also optimize the fluid properties. Our approach
utilizes an adjoint-based approach to compute the sensitivity of a prescribed quantity-ofinterest to the elastic moduli which is used to drive a gradient-based optimization algorithm.
We have demonstrated the effective use of our formulation for obtaining the spatial distribution of the elastic modulus in a structure to achieve optimality of interior and exterior
problems. We show that the formulation easily allows for spatially continuous or piecewise
constant distribution of the elastic modulus through the use of appropriate regularization
107

terms added to the loss function. This is significant since for many practical applications
one may be constrained to not change the modulus in portions of the structure and require
sharp discontinuities in modulus distribution.
In Chapter 3 we next developed a computational formulation for shape optimization,
again leveraging an adjoint-based approach. Notably, this method requires no additional
surface parameterization; instead, it operates directly on a conventional finite element mesh,
eliminating the need for specialized geometry-handling techniques. We validated the framework by applying it to a noise abatement problem, demonstrating that shape modifications
guided by adjoint sensitivities can effectively reduce unwanted acoustic response. Compared
to prior methods using parametric or isogeometric representations, our approach underscores
the ease of integration with a conventional mesh-based FEM workflow.
The main drawback of the proposed adjoint based approaches is that, due to the nonconvex nature of the optimization landscape, the algorithm may converge to a local optimum
unless guided by a well-chosen initial guess that is based on empirical knowledge. Moreover,
for large-scale, particularly three-dimensional acoustic problems, iterative forward solvers
may suffer from severe ill-conditioning; on the other hand, direct solvers demand substantial
memory creating a computational bottleneck. Thus, there is a need for more robust and
efficient strategies for numerically solving forward problems.
Finally in Chapter 4, we shifted focus to a probabilistic framework for inverse problems,
highlighting how conditional score based diffusion models (cSDM) can be leveraged to solve
Bayesian inference in high-dimensional setting. By training on synthetic data produced by
black-box forward solvers, these models accommodate complex, high-dimensional priors and
require no re-training when new measurements become available; only a sampling procedure
via Annealed Langevin dynamics is needed. This enables near-real-time inference and facilitates uncertainty quantification, even for problems involving substantial measurement noise.
We demonstrated the effectiveness of this framework on a range of inverse problems, which
includes problems in acoustic scattering and elastography. We validated the proposed ap108

proach on two synthetic and one real-life examples. The quality of inference was comparable
to Monte Carlo simulations and aligned well with intuitive expectations.
In the context of inverse scattering problems, the cSDM approach effectively infers the
topology of multiple scatterers in the domain of interest while accommodating arbitrary
interface boundaries, a common challenge in structural acoustic problems, thereby removing the need for methods like level-sets or explicit interface constraints. In addition, we
demonstrated the practical utility of the approach in a biomedical context by applying it
to elastography imaging of tumor spheroids and validating the results with experimental
measurements. Overall, we showed that conditional score-based diffusion models can tackle
high-dimensional inverse problems involving complex prior distributions, different types of
forward physics (treated as black-box models), complex measurement noise models, and various measurement modalities. Despite the promise shown by these results, several important
research directions remain; primarily, the need to address the computational costs associated with generating large training datasets via repeated simulations. The development of
data-efficient conditional diffusion models is therefore a key challenge for the wider adoption
of this technology in scientific and engineering applications.
5.2 Future directions
There are promising avenues for future work to extend the deterministic methods discussed
in this thesis. Currently, these methods focus on reconstructing specific parameters such as
stiffness; however, they can be adapted to estimate other relevant material properties, including density and Poisson’s ratio within the domain. Such enhancements would be particularly
useful in designing advanced materials, such as acoustic cloaks, where both geometric and
material optimization are essential. In addition, incorporating global optimization algorithms could help overcome local minima issues, ensuring more robust convergence. In the
case of shape optimization, from a computational geometry perspective, there is a clear need
109

for better mesh movement techniques that are capable of handling large interface boundary deformations such as affine transformations (e.g., translations and rotations) and other
complex boundary modifications, without resorting to expensive re-meshing strategies.
The cSDM approach presented in this thesis also leads to several future research directions. First, incorporating physics into the diffusion model either through modified loss
functions, specialized network architectures, or physics-informed sampling techniques, could
further improve accuracy of the prediction. Second, adopting multi-fidelity datasets, which
have lower accuracy but cheaper simulations that augment high-fidelity data, could significantly reduce training costs and enhance scalability to large, real-world problems. Third,
using methods such as masking-based training [89] could enable reconstructions from sparse
measurements; these are more common in practical scenarios than the full-field measurements
examined here. Fourth, developing more realistic noise models, possibly derived empirically
from experiments, could significantly improve the reliability of the inversion process, especially in scenarios where the noise is complex or non-Gaussian; this is because we observed
that the cSDM’s uncertainty predictions depend heavily on the accuracy of the underlying
noise model.
On a final note, the application of diffusion-based methods to biomedical imaging, particularly optical coherence elastography (OCE) [90], presents a highly promising research direction. By incorporating more sophisticated tumor prior models that are tailored to specific
tissue types and growth patterns, a significant improvement in elastography-based reconstructions can be achieved. This would help bridge the gap between a purely computational
technique and a real-world clinical procedure. Our approach has already demonstrated improved performance compared to the state-of-the-art algebraic method (see Figure 5.1) [74].
In the algebraic method, the strain is estimated from optical phase-difference images by
calculating gradients in the vertical direction and inferring the modulus through a reciprocal strain value based on a simplistic uniform-stress assumption. In contrast, our deep
learning–based generative model captures heterogeneous material property distributions and
110

Modulus Strain
State-of-the-art Our method
Figure 5.1: Comparison between the the results obtained using the state-of-the-art method
and the cSDM approach for the tumor spheroid problem in Section 4.5.1.2. Note, in the
existing method, the inferred quantity (Young’s modulus) is directly estimated by computing
the strain field; hence uncertainty measure is unavailable. Algebraic results from Foo et al.
provides uncertainty estimates, leading to a more accurate and comprehensive characterization. With the ability to map mechanical properties at µm − mm scales, OCE can facilitate
deeper insights into cellular-environment interactions. This makes OCE a powerful tool for
111

accurate and early disease detection, particularly in fields such as ophthalmology, oncology,
and cardiology. With the integration of diffusion-based models into the OCE, this thesis
highlights a crucial advance in delivering robust, high-resolution elastography reconstructions with explicit uncertainty quantification. Ultimately, this approach has the potential to
transform clinical practice by enabling more accurate disease detection and providing deeper
insights into tissue biomechanics in a wide range of medical applications.
112

References
[1] Semyung Wang and Jeawon Lee. “Acoustic design sensitivity analysis and optimization
for reduced exterior noise”. In: AIAA journal 39.4 (2001). doi:10.2514/2.1376, pp. 574–
580.
[2] George T Gebhardt. “Acoustical design features of Boeing model 727”. In: Journal of
Aircraft 2.4 (1965). doi:10.2514/3.43652, pp. 272–277.
[3] Jun Dong, Kyung K Choi, and Nam H Kim. “Design optimization for structuralacoustic problems using FEA-BEA with adjoint variable method”. In: Journal of Mechanical Design 126.3 (2004). doi:10.1115/1.1701879, pp. 527–533.
[4] Marco Miniaci et al. “Large scale mechanical metamaterials as seismic shields”. In:
New Journal of Physics 18.8 (2016). doi:10.1088/1367-2630/18/8/083041, p. 083041.
[5] Antonio Palermo et al. “Engineered metabarrier as shield from seismic surface waves”.
In: Scientific reports 6 (2016). doi:10.1038/srep39356, p. 39356.
[6] ZF Shi, ZB Cheng, and C Xiong. “A new seismic isolation method by using a periodic foundation”. In: Earth and Space 2010: Engineering, Science, Construction, and
Operations in Challenging Environments. doi:10.1061/41096(366)242. 2010, pp. 2586–
2594.
[7] Sebastian Kr¨odel, Nicolas Thom´e, and Chiara Daraio. “Wide band-gap seismic metastructures”. In: Extreme Mechanics Letters 4 (2015). doi:10.1016/j.eml.2015.05.004,
pp. 111–117.
[8] CW Lim et al. “From photonic crystals to seismic metamaterials: A review via phononic
crystals and acoustic metamaterials”. In: Archives of Computational Methods in Engineering (2021), pp. 1–62.
[9] Lucian Zigoneanu, Bogdan-Ioan Popa, and Steven A Cummer. “Three-dimensional
broadband omnidirectional acoustic ground cloak”. In: Nature materials 13.4 (2014),
pp. 352–355.
[10] Shu Zhang, Leilei Yin, and Nicholas Fang. “Focusing ultrasound with an acoustic
metamaterial network”. In: Physical review letters 102.19 (2009), p. 194301.
[11] Meng Chen et al. “Design of an acoustic superlens using single-phase metamaterials
with a star-shaped lattice structure”. In: Scientific reports 8.1 (2018), pp. 1–8.
113

[12] Yongyao Chen et al. “Enhanced acoustic sensing through wave compression and pressure amplification in anisotropic metamaterials”. In: Nature communications 5.1 (2014),
pp. 1–9.
[13] Rosa MS Sigrist et al. “Ultrasound elastography: review of techniques and clinical
applications”. In: Theranostics 7.5 (2017), p. 1303.
[14] CH Greene et al. “Acoustical detection of high-density krill demersal layers in the
submarine canyons off Georges Bank”. In: Science 241.4863 (1988), pp. 359–361.
[15] SY Wang, Kang Tai, and MY2222890 Wang. “An enhanced genetic algorithm for
structural topology optimization”. In: International Journal for Numerical Methods in
Engineering 65.1 (2006), pp. 18–44.
[16] Ahmed Mostafa Shaaban et al. “Shape optimization by conventional and extended isogeometric boundary element method with PSO for two-dimensional Helmholtz acoustic
problems”. In: Engineering Analysis with Boundary Elements 113 (2020), pp. 156–169.
[17] NP Garcia-Lopez et al. “A hybrid topology optimization methodology combining simulated annealing and SIMP”. In: Computers & structures 89.15-16 (2011), pp. 1512–
1522.
[18] Qiyin Lin et al. “Investigation into the topology optimization for conductive heat
transfer based on deep learning approach”. In: International Communications in Heat
and Mass Transfer 97 (2018), pp. 103–109.
[19] Dalei Wang et al. “A deep convolutional neural network for topology optimization with
perceptible generalization ability”. In: Engineering Optimization 54.6 (2022), pp. 973–
988.
[20] Erva Ulu, Rusheng Zhang, and Levent Burak Kara. “A data-driven investigation and
estimation of optimal topologies under variable loading configurations”. In: Computer
Methods in Biomechanics and Biomedical Engineering: Imaging & Visualization 4.2
(2016), pp. 61–72.
[21] Yonggyun Yu et al. “Deep learning for determining a near-optimal topological design
without any iteration”. In: Structural and Multidisciplinary Optimization 59.3 (2019),
pp. 787–799.
[22] Walter Jesus Paucar Casas and Renato Pavanello. “Optimization of fluid-structure systems by eigenvalues gap separation with sensitivity analysis”. In: Applied Mathematical
Modelling 42 (2017). doi:10.1016/j.apm.2016.10.031, pp. 269–289.
114

[23] Jianbin Du and Niels Olhoff. “Minimization of sound radiation from vibrating bimaterial structures using topology optimization”. In: Structural and Multidisciplinary
Optimization 33.4-5 (2007). doi:10.1007/s00158-006-0088-9, pp. 305–321.
[24] Xiaopeng Zhang and Zhan Kang. “Topology optimization of damping layers for minimizing sound radiation of shell structures”. In: Journal of Sound and Vibration 332.10
(2013). doi:10.1016/j.jsv.2012.12.022, pp. 2500–2519.
[25] Linyuan Shang and Guozhong Zhao. “Optimality criteria-based topology optimization of a bi-material model for acoustic–structural coupled systems”. In: Engineering
Optimization 48.6 (2016). doi:10.1080/0305215X.2015.1082351, pp. 1060–1079.
[26] Wenchang Zhao et al. “Minimization of sound radiation in fully coupled structuralacoustic systems using FEM-BEM based topology optimization”. In: Structural and
Multidisciplinary Optimization 58.1 (2018). doi:10.1007/s00158-017-1881-3, pp. 115–
128.
[27] Wenchang Zhao et al. “Topology optimization of exterior acoustic-structure interaction systems using the coupled FEM-BEM method”. In: International Journal for
Numerical Methods in Engineering 119.5 (2019). doi:10.1002/nme.6055, pp. 404–431.
[28] B Kieback, A Neubrand, and H Riedel. “Processing techniques for functionally graded
materials”. In: Materials Science and Engineering: A 362.1-2 (2003), pp. 81–106.
[29] Minoo Naebe and Kamyar Shirvanimoghaddam. “Functionally graded materials: A
review of fabrication and properties”. In: Applied materials today 5 (2016), pp. 223–
245.
[30] Tao Fu et al. “Vibro-acoustic characteristics of eccentrically stiffened functionally
graded material sandwich cylindrical shell under external mean fluid”. In: Applied
Mathematical Modelling 91 (2021), pp. 214–231.
[31] Raphael T Haftka and Ramana V Grandhi. “Structural shape optimization—a survey”.
In: Computer methods in applied mechanics and engineering 57.1 (1986), pp. 91–106.
[32] Jamshid A Samareh. “Survey of shape parameterization techniques for high-fidelity
multidisciplinary shape optimization”. In: AIAA journal 39.5 (2001), pp. 877–884.
[33] Gonzalo R Feij´oo, Assad A Oberai, and Peter M Pinsky. “An application of shape optimization in the solution of inverse acoustic scattering problems”. In: Inverse problems
20.1 (2003), p. 199.
[34] Gilles PL Thomas et al. “Parametric shape optimization of lens-focused piezoelectric ultrasound transducers”. In: IEEE transactions on ultrasonics, ferroelectrics, and
frequency control 65.5 (2018), pp. 844–850.
115

[35] Gil Ho Yoon, Jakob Søndergaard Jensen, and Ole Sigmund. “Topology optimization of acoustic–structure interaction problems using a mixed finite element formulation”. In: International journal for numerical methods in engineering 70.9 (2007).
doi:10.1002/nme.1900, pp. 1049–1075.
[36] J Austin Cottrell et al. “Isogeometric analysis of structural vibrations”. In: Computer
methods in applied mechanics and engineering 195.41-43 (2006), pp. 5257–5296.
[37] LL Chen et al. “Structural shape optimization of three dimensional acoustic problems with isogeometric boundary element methods”. In: Computer Methods in Applied
Mechanics and Engineering 355 (2019), pp. 926–951.
[38] Ahmed Mostafa Shaaban et al. “3D isogeometric boundary element analysis and structural shape optimization for Helmholtz acoustic scattering problems”. In: Computer
Methods in Applied Mechanics and Engineering 384 (2021), p. 113950.
[39] Eddie Wadbro and Martin Berggren. “Topology optimization of an acoustic horn”. In:
Computer methods in applied mechanics and engineering 196.1-3 (2006). doi:10.1016/
j.cma.2006.05.005, pp. 420–436.
[40] Leilei Chen et al. “An isogeometric approach of two dimensional acoustic design
sensitivity analysis and topology optimization analysis for absorbing material distribution”. In: Computer Methods in Applied Mechanics and Engineering 336 (2018).
doi:10.1016/j.cma.2018.03.025, pp. 507–532.
[41] Maria B D¨uhring, Jakob S Jensen, and Ole Sigmund. “Acoustic design by topology optimization”. In: Journal of sound and vibration 317.3-5 (2008). doi:10.1016/j.jsv.2008.
03.042, pp. 557–575.
[42] Jacob Andkjær and Ole Sigmund. “Topology optimized cloak for airborne sound”. In:
Journal of Vibration and Acoustics 135.4 (2013). doi:10.1115/1.4023828, p. 041011.
[43] Seongyeol Goo et al. “Topology optimization of bounded acoustic problems using the
hybrid finite element-wave based method”. In: Computer Methods in Applied Mechanics and Engineering 313 (2017). doi:10.1016/j.cma.2016.10.027, pp. 834–856.
[44] Yuki Noguchi et al. “Optimum design of an acoustic metamaterial with negative bulk
modulus in an acoustic-elastic coupled system using a level set–based topology optimization method”. In: International Journal for Numerical Methods in Engineering
113.8 (2018). doi:10.1002/nme.5616, pp. 1300–1339.
[45] Kohei Miyata et al. “Optimum design of a multi-functional acoustic metasurface using topology optimization based on Zwicker’s loudness model”. In: Computer Methods in Applied Mechanics and Engineering 331 (2018). doi:10.1016/j.cma.2017.11.017,
pp. 116–137.
116

[46] Niels Olhoff and Jianbin Du. “Generalized incremental frequency method for topological design of continuum structures for minimum dynamic compliance subject to forced
vibration at a prescribed low or high value of the excitation frequency”. In: Structural and Multidisciplinary Optimization 54.5 (2016). doi:10.1007/s00158-016-1574-3,
pp. 1113–1141.
[47] Junghwan Kook. “Evolutionary topology optimization for acoustic-structure interaction problems using a mixed u/p formulation”. In: Mechanics Based Design of Structures and Machines 47.3 (2019). doi:10.1080/15397734.2018.1557527, pp. 356–374.
[48] Antonio Andr´e Novotny et al. “Topological sensitivity analysis”. In: Computer methods
in applied mechanics and engineering 192.7-8 (2003), pp. 803–829.
[49] Lei Shu, Michael Yu Wang, and Zhengdong Ma. “Level set based topology optimization
of vibrating structures for coupled acoustic–structural dynamics”. In: Computers &
Structures 132 (2014). doi:10.1016/j.compstruc.2013.10.019, pp. 34–42.
[50] WM Vicente et al. “Topology optimization of frequency responses of fluid–structure interaction systems”. In: Finite Elements in Analysis and Design 98 (2015). doi:10.1016/
j.finel.2015.01.009, pp. 1–13.
[51] Cetin B Dilgen and Niels Aage. “Generalized shape optimization of transient vibroacoustic problems using cut elements”. In: International Journal for Numerical Methods
in Engineering 122.6 (2021), pp. 1578–1601.
[52] Yuehaw Khoo and Lexing Ying. “SwitchNet: a neural network model for forward and
inverse scattering problems”. In: SIAM Journal on Scientific Computing 41.5 (2019),
A3182–A3201.
[53] Yuwei Fan and Lexing Ying. “Solving inverse wave scattering with deep learning”. In:
arXiv preprint arXiv:1911.13202 (2019).
[54] Francesco Ballarin and Gianluigi Rozza. multiphenics - Easy prototyping of multiphysics problems in FEniCS. https://mathlab.sissa.it/multiphenics.
[55] Pauli Virtanen et al. “SciPy 1.0–Fundamental Algorithms for Scientific Computing in
Python”. In: arXiv e-prints, arXiv:1907.10121 (2019), arXiv:1907.10121. arXiv: 1907.
10121 [cs.MS].
[56] Joyce R McLaughlin and Jeong-Rock Yoon. “Unique identifiability of elastic parameters from time-dependent interior displacement measurement”. In: Inverse Problems
20.1 (2003). 10.1088/0266-5611/20/1/002, p. 25.
117

[57] Paul E Barbone and Nachiket H Gokhale. “Elastic modulus imaging: on the uniqueness
and nonuniqueness of the elastography inverse problem in two dimensions”. In: Inverse
problems 20.1 (2004). 10.1088/0266-5611/20/1/017, p. 283.
[58] Paul E Barbone and Assad A Oberai. “Elastic modulus imaging: some exact solutions
of the compressible elastography inverse problem”. In: Physics in Medicine & Biology
52.6 (2007). 10.1088/0031-9155/52/6/003, p. 1577.
[59] Steven G Johnson. “Notes on perfectly matched layers (PMLs)”. In: arXiv preprint
arXiv:2108.05348 (2021).
[60] Deep Ray et al. “The efficacy and generalizability of conditional GANs for posterior
inference in physics-based inverse problems”. In: Numerical Algebra, Control and Optimization 14.1 (2024), pp. 160–189. issn: 2155-3289. doi: 10.3934/naco.2022038.
url: https://www.aimsciences.org/article/id/639862c5b2114e413cb35cd4.
[61] Hoang Thanh-Tung and Truyen Tran. “Catastrophic forgetting and mode collapse
in GANs”. In: 2020 international joint conference on neural networks (ijcnn). IEEE.
2020, pp. 1–10.
[62] Dhruv Patel and Assad A Oberai. “Bayesian inference with generative adversarial
network priors”. In: arXiv preprint arXiv:1907.09987 (2019).
[63] Yang Song et al. “Score-based generative modeling through stochastic differential equations”. In: arXiv preprint arXiv:2011.13456 (2020).
[64] Mark Kac. “On distributions of certain Wiener functionals”. In: Transactions of the
American Mathematical Society 65.1 (1949), pp. 1–13.
[65] Yang Song and Stefano Ermon. “Generative modeling by estimating gradients of the
data distribution”. In: Advances in neural information processing systems 32 (2019).
[66] Aapo Hyv¨arinen and Peter Dayan. “Estimation of non-normalized statistical models
by score matching.” In: Journal of Machine Learning Research 6.4 (2005).
[67] Agnimitra Dasgupta et al. “Conditional score-based diffusion models for solving inverse
elasticity problems”. In: Computer Methods in Applied Mechanics and Engineering 433
(2025), p. 117425.
[68] Yang Song and Stefano Ermon. “Improved techniques for training score-based generative models”. In: Advances in neural information processing systems 33 (2020),
pp. 12438–12448.
[69] Diederik P Kingma and Jimmy Ba. “Adam: A method for stochastic optimization”.
In: arXiv preprint arXiv:1412.6980 (2014).
118

[70] Bradley Efron. “Tweedie’s formula and selection bias”. In: Journal of the American
Statistical Association 106.496 (2011), pp. 1602–1614.
[71] Alexia Jolicoeur-Martineau et al. “Adversarial score matching and improved sampling
for image generation”. In: arXiv preprint arXiv:2009.05475 (2020).
[72] Martin Heusel et al. “Gans trained by a two time-scale update rule converge to a local
nash equilibrium”. In: Advances in neural information processing systems 30 (2017).
[73] Guosheng Lin et al. “Refinenet: Multi-path refinement networks for high-resolution
semantic segmentation”. In: Proceedings of the IEEE conference on computer vision
and pattern recognition. 2017, pp. 1925–1934.
[74] Ken Y Foo et al. “Tumor spheroid elasticity estimation using mechano-microscopy
combined with a conditional generative adversarial network”. In: Computer Methods
and Programs in Biomedicine 255 (2024), p. 108362.
[75] Adam Paszke et al. “Pytorch: An imperative style, high-performance deep learning
library”. In: Advances in neural information processing systems 32 (2019).
[76] Paul E Barbone and Assad A Oberai. “A review of the mathematical and computational foundations of biomechanical imaging”. In: Computational Modeling in Biomechanics (2010), pp. 375–408.
[77] Anders Logg, Kent-Andre Mardal, Garth N. Wells, et al. Automated Solution of Differential Equations by the Finite Element Method. Springer, 2012. doi: 10.1007/978-
3-642-23099-8.
[78] B. F. Kennedy, K. M. Kennedy, and D. D. Sampson. “A Review of Optical Coherence
Elastography: Fundamentals, Techniques and Prospects”. In: IEEE Journal of Selected
Topics in Quantum Electronics 20.2 (2014), pp. 272–288. issn: 1077-260X. doi: 10.
1109/JSTQE.2013.2291445.
[79] Brendan F Kennedy et al. “Optical coherence micro-elastography: mechanical-contrast
imaging of tissue microstructure”. In: Biomedical optics express 5.7 (2014), pp. 2113–
2124.
[80] Alireza Mowla et al. “Subcellular Mechano-Microscopy: High Resolution Three Dimensional Elasticity Mapping Using Optical Coherence Microscopy”. In: Optics Letters
47.13 (2022), pp. 3303–3306. issn: 1539-4794. doi: 10.1364/OL.451681.
[81] Alireza Mowla et al. “Multimodal Mechano-Microscopy Reveals Mechanical Phenotypes of Breast Cancer Spheroids in Three Dimensions”. In: bioRxiv (2024), p. 2024.04.
05.588260. doi: 10.1101/2024.04.05.588260.
119

[82] Karl Olofsson et al. “Acoustic formation of multicellular tumor spheroids enabling onchip functional and structural imaging”. In: Lab on a Chip 18.16 (2018), pp. 2466–
2476.
[83] Hiroshi Imai, Masao Iri, and Kazuo Murota. “Voronoi Diagram in the Laguerre Geometry and Its Applications”. In: SIAM Journal on Computing 14.1 (1985), pp. 93–105.
doi: 10.1137/0214006.
[84] Dassault Syst`emes. Simulia 2022 User Assistance. 2022.
[85] Elizabete Rodrigues Ferreira, Assad A Oberai, and Paul E Barbone. “Uniqueness of
the elastography inverse problem for incompressible nonlinear planar hyperelasticity”.
In: Inverse problems 28.6 (2012), p. 065008.
[86] Hermann Schomberg. “An improved approach to reconstructive ultrasound tomography”. In: Journal of Physics D: Applied Physics 11.15 (1978), p. L181.
[87] Nicholas Rawlinson, S Pozgay, and S Fishwick. “Seismic tomography: A window into
deep Earth”. In: Physics of the Earth and Planetary Interiors 178.3-4 (2010), pp. 101–
135.
[88] Borong Zhang, Leonardo Zepeda-N´u˜nez, and Qin Li. “Solving the wide-band inverse
scattering problem via equivariant neural networks”. In: Journal of Computational and
Applied Mathematics 451 (2024), p. 116050.
[89] Kaiming He et al. “Masked autoencoders are scalable vision learners”. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. 2022,
pp. 16000–16009.
[90] Kelsey M Kennedy et al. “Quantitative micro-elastography: imaging of tissue elasticity
using compression optical coherence elastography”. In: Scientific reports 5.1 (2015),
p. 15538.
120

Abstract (if available)

Abstract This thesis explores novel methodologies in the estimation of material properties, shape, and topology of the structure for coupled structural acoustics systems by leveraging both traditional optimization techniques and innovative conditional inference approaches. The motivation for this work stems from the increasing demand for robust and cost-effective techniques for design applications, where accurate modeling of the structural geometry along with interface dynamics is critical. The techniques presented in work can be employed in a variety of applications such as vibrations of submerged structures and fluid filled containers, acoustic radiation scattering, functional acoustic metamaterials design and acoustic medical sensors. They can also be extended to other wave-propagation domains including electromagnetic propagation and scattering. This thesis can be broadly divided into two parts: Part 1 utilizes conventional deterministic approaches, including sensitivity analysis, to solve material and shape optimization problems in structural acoustics; Part 2 employs probabilistic methods and state-of-the-art generative machine learning models, specifically diffusion models, to solve inverse problems in elastography and acoustic scattering.

The first part of this dissertation focuses on semi-analytical deterministic approaches for solving material and shape optimization problems in coupled fluid-structure acoustics under time-harmonic loads. The fluid and structural components are modeled using the Helmholtz and Navier-Cauchy equations, respectively, with interface conditions explicitly enforced. We introduce a computationally efficient and generalized formulation for computing adjoint-based sensitivities, which drive gradient-based minimization algorithms for both material and shape optimization. This approach is robust, computationally efficient, and seamlessly integrates into conventional finite element method (FEM) frameworks.

In material optimization, we explore continuous and bi-material property distributions, leveraging recent advancements in manufacturing functionally graded materials. These advancements have expanded the design space, allowing for more sophisticated and feasible material configurations.

For shape optimization, we explore parameterization-free methods, known as node-based sensitivity approaches. These approaches encompass both level-set and artificial density-based methods that utilize conventional meshing techniques. However, level-set methods often face numerical stability and initialization challenges, particularly in 3D applications, while density-based methods tend to lack the sharp interfaces required for applications where precise interface definitions are critical, such as in airfoil design. Although Isogeometric Analysis (IGA) is often the preferred method for shape optimization due to its integration of Computer-Aided Design (CAD) and Finite Element Analysis (FEA), the industrial adoption of such advanced techniques remains in its nascent stages. To address these limitations, we propose a novel method that integrates smoothly within conventional frameworks. A key challenge in node-based sensitivity approaches, particularly for coupled physics problems, is computing the shape derivative, which involves complex terms such as surface curvature and the gradient of interface-specific quantities. Additionally, mesh updates in each iteration often lead to distortions and artifacts, a recurring issue in node-based methods. Although numerous studies have addressed mesh updating challenges, effective computation of sensitivities for interface problems remain underexplored in the literature. To overcome this challenge, we introduce a novel technique that constructs a smooth extension of the surface normal vector field across an auxiliary domain intersecting the interface, enabling us to compute curvature effectively without additional parameterization. The efficacy of this approach is demonstrated through its application in an interior noise minimization problem in structural acoustics.

The second part of this dissertation transitions to a data-driven probabilistic approach to address inverse problems, with a specific focus on inferring the topology of acoustic scatterers, a critical component in underwater sensing applications such as SONAR. By utilizing simulation-based inference, we sample from posterior distributions linked to an unknown forward model, employing state-of-the-art generative machine learning models, notably conditional diffusion models, to learn the structure of these distributions. Although generative models are proficient at handling geometrically complex priors, they encounter challenges due to the high dimensionality and nonlinearities inherent in the forward model. These complexities can restrict the comprehensive exploration of the distribution space, potentially resulting in predictions confined to regions around local modes, a phenomenon commonly referred to as mode collapse. To overcome these limitations, we introduce a framework employing conditional score-based diffusion models for posterior inference. These models excel at approximating the score function of a posterior distribution using samples from an extensive dataset of pairwise data generated via finite element forward solvers. This setup demonstrates our method's capability to manage black-box forward models and intricate measurement noise effectively. Uncertainty quantification is accomplished by computing posterior statistics from numerous realizations sampled through Langevin dynamics driven by the learned score function. The efficacy of this approach is first validated through applications in elastography, employing both synthetic and actual experimental data. We demonstrate the model’s ability at handling diverse measurement modalities, sophisticated geometries, and nonlinear noise models, by comparing their performance to traditional Monte Carlo and algebraic methods. We subsequently apply this methodology to solve the inverse scattering problem for three distinct sensor configurations and derive intuitive insights from them.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Physics-based data-driven inference

PDF

Efficient inverse analysis with dynamic and stochastic reductions for large-scale models of multi-component systems

PDF

Deep learning application for fiber-reinforced composites with uncertainty

PDF

Solution of inverse scattering problems via hybrid global and local optimization

PDF

Novel optimization tools for structured signals recovery: channels estimation and compressible signal recovery

PDF

Reconstruction and estimation of shear flows using physics-based and data-driven models

PDF

Robust optimization under exogenous and endogenous uncertainty: methodology, software, and applications

PDF

Data-driven optimization for indoor localization

PDF

Optimization methods and algorithms for constrained magnetic resonance imaging

PDF

Human motion data analysis and compression using graph based techniques

PDF

Parameter estimation to infer injector-producer relationships in oil fields: from hybrid constrained nonlinear optimization to inference in probabilistic graphical model

PDF

Model falsification: new insights, methods and applications

PDF

Impurity gas detection for spent nuclear fuel (SNF) canisters using ultrasonic sensing and deep learning

PDF

A high-fidelity geometric contact method for thin fabrics using shell finite elements in explicit and implicit time scheme

PDF

Data-driven multi-fidelity modeling for physical systems

PDF

Design optimization under uncertainty for rotor blades of horizontal axis wind turbines

PDF

Scalable exact inference in probabilistic graphical models on multi-core platforms

PDF

Damage identification in spent nuclear fuel canisters using dynamic modal analysis and machine learning

PDF

Understanding diffusion process: inference and theory

PDF

Advances in linguistic data-oriented uncertainty modeling, reasoning, and intelligent decision making

Asset Metadata

Creator Ramaswamy, Harisankar (author)

Core Title Material, shape, and topology estimation in coupled structural acoustics: sensitivity-based and data-driven approaches

School Andrew and Erna Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Mechanical Engineering

Degree Conferral Date 2025-05

Publication Date 04/07/2025

Defense Date 09/12/2024

Publisher University of Southern California (original), Los Angeles, California (original), University of Southern California. Libraries (digital)

Tag design optimization,diffusion models,inverse problem,OAI-PMH Harvest,probabilistic inference,sensitivity analysis,structural acoustics

Format theses (aat)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Oberai, Assad (committee chair), Plucinsky, Paul (committee member), Gencturk, Bora (committee member)

Creator Email hramaswa@usc.edu,hsankar94@gmail.com

Unique identifier UC11399K9VN

Identifier etd-RamaswamyH-13903.pdf (filename)

Legacy Identifier etd-RamaswamyH-13903

Document Type Dissertation

Format theses (aat)

Rights Ramaswamy, Harisankar

Internet Media Type application/pdf

Type texts

Source 20250409-usctheses-batch-1250 (batch), University of Southern California Dissertations and Theses (collection), University of Southern California (contributing entity)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu