This tutorial extends the previous 1D wave equation example and solve a 2D seismic wave propagation problem commonly used in seismic surveying. In this tutorial, you will learn the following:

How to solve a 2D time-dependent problem in Modulus

How to define an open boundary condition using custom equations

How to present a variable velocity model as an additional network

This tutorial assumes that you have completed the Introductory Example tutorial on Lid Driven Cavity flow and have familiarized yourself with the basics of Modulus. Also, see the 1D Wave Equation tutorial for more information on defining new differential equations, and solving time-dependent problems in Modulus.

The acoustic wave equation for the square slowness, defined as \(m=1/c^2\) where \(c\) is the speed of sound (velocity) of a given physical medium with constant density, and a source \(q\) is given by:

(140)\[u_{tt} = c^2u_{xx} + c^2u_{yy} + q \quad \mbox{ in } \Omega\]

Where \(u(\mathbf{x},t)\) represents the pressure response (known as the “wavefield”) at location vector \(\mathbf{x}\) and time \(t\) in an acoustic medium. Despite its linearity, the wave equation is notoriously challenging to solve in complex media, because the dynamics of the wavefield at the interfaces of the media can be highly complex, with multiple types of waves with large range of amplitudes and frequencies interfering simultaneously.

In this tutorial, you will solve the 2D acoustic wave equation with a single Ricker Source in a layered velocity model, 1.0 \(km/s\) at the top layer and 2.0 \(km/s\) the bottom (Fig. 60). Sources in seismic surveys are positioned at a single or a few physical locations where artificial pressure is injected into the domain you want to model. In the case of land survey, it is usually dynamite blowing up at a given location, or a vibroseis (a vibrating plate generating continuous sound waves). For a marine survey, the source is an air gun sending a bubble of compressed air into the water that will expand and generate a seismic wave.

This problem, uses a domain size 2 \(km\) x 2 \(km\), and a single source is located at the center of the domain. The source signature is modelled using a Ricker wavelet, illustrated in Fig. 61, with a peak wavelet frequency of \(f_0=15 Hz\).

Fig. 60 *Velocity model of the medium*

Fig. 61 *Ricker source signature*

The problem uses wavefield data at time steps (150 \(ms\) - 300 \(ms\)) generated from finite difference simulations, using Devito, as constraints for the temporal boundary conditions, and train Modulus to produce wavefields at later time steps (300 \(ms\) – 1000 \(ms\)).

The setup for this problem is similar to tutorial 1D Wave Equation. So, only the main highlights of this problem is discussed here.

The full python script for this problem can be found at `examples/seismic_wave/wave_2d.py`

A second-order PDE requires strict BC on both the initial wavefield and its derivatives for its solution to be unique. In the field, the seismic wave propagates in every direction to an “infinite” distance. In the finite computational domain, Absorbing BC (ABC) or Perfectly Matched Layers (PML) are artificial boundary conditions that are typically applied to approximate an infinite media by damping and absorbing the waves at the edges of the domain, in order to avoid reflections from the boundary. However, a NN solver is meshless – it is not suitable to implement ABC or PML. To enable a wave to leave the computational domain and travel undisturbed through the boundaries, you can apply an open boundary condition, also called a radiation condition, which imposes the first-order PDEs at the boundaries:

(141)\[\begin{split}\begin{split}
\frac{\partial u}{\partial t} - c\frac{\partial u}{\partial x} &= 0,\quad \mbox{ at } x = 0 \\
\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} &= 0,\quad \mbox{ at } x = dLen \\
\frac{\partial u}{\partial t} - c\frac{\partial u}{\partial y} &= 0,\quad \mbox{ at } y = 0 \\
\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial y} &= 0,\quad \mbox{ at } y = dLen \\
\end{split}\end{split}\]

Previous tutorials have described how to define custom PDEs. Similarly
here you create a class `OpenBoundary`

that inherits from `PDES`

:

```
```*# define open boundary conditions*
class OpenBoundary(PDE):
"""
Open boundary condition for wave problems
Ref: http://hplgit.github.io/wavebc/doc/pub/._wavebc_cyborg002.html
Parameters
==========
u : str
The dependent variable.
c : float, Sympy Symbol/Expr, str
Wave speed coefficient. If `c` is a str then it is
converted to Sympy Function of form 'c(x,y,z,t)'.
If 'c' is a Sympy Symbol or Expression then this
is substituted into the equation.
dim : int
Dimension of the wave equation (1, 2, or 3). Default is 2.
time : bool
If time-dependent equations or not. Default is True.
"""
name = "OpenBoundary"
def __init__(self, u="u", c="c", dim=3, time=True):
*# set params*
self.u = u
self.dim = dim
self.time = time
*# coordinates*
x, y, z = Symbol("x"), Symbol("y"), Symbol("z")
*# normal*
normal_x, normal_y, normal_z = (
Symbol("normal_x"),
Symbol("normal_y"),
Symbol("normal_z"),
)
*# time*
t = Symbol("t")
*# make input variables*
input_variables = {"x": x, "y": y, "z": z, "t": t}
if self.dim == 1:
input_variables.pop("y")
input_variables.pop("z")
elif self.dim == 2:
input_variables.pop("z")
if not self.time:
input_variables.pop("t")
*# Scalar function*
assert type(u) == str, "u needs to be string"
u = Function(u)(*input_variables)
*# wave speed coefficient*
if type(c) is str:
c = Function(c)(*input_variables)
elif type(c) in [float, int]:
c = Number(c)
*# set equations*
self.equations = {}
self.equations["open_boundary"] = (
u.diff(t)
+ normal_x * c * u.diff(x)
+ normal_y * c * u.diff(y)
+ normal_z * c * u.diff(z)
)

In the 1D Wave Equation tutorial, the velocity (\(c\)) of the physical medium is constant. In this problem, you have a velocity model that varies with locations \(x\). Use a \(tanh\) function form to represent the velocity \(c\) as a function of \(x\) and \(y\):

A \(tanh\) function was used to avoid the sharp discontinuity at the interface. Such smoothing of sharp boundaries helps the convergence of the Neural network solver.

```
```*# define target velocity model*
*# 2.0 km/s at the bottom and 1.0 km/s at the top using tanh function*
mesh_x, mesh_y = np.meshgrid(
np.linspace(0, 2, 512), np.linspace(0, 2, 512), indexing="ij"
)
wave_speed_invar = {}
wave_speed_invar["x"] = np.expand_dims(mesh_x.flatten(), axis=-1)
wave_speed_invar["y"] = np.expand_dims(mesh_y.flatten(), axis=-1)
wave_speed_outvar = {}
wave_speed_outvar["c"] = np.tanh(80 * (wave_speed_invar["y"] - 1.0)) / 2 + 1.5

First, define the Modulus PDE and NN nodes for this problem:

```
```*# define networks and nodes*
wave_net = instantiate_arch(
input_keys=[Key("x"), Key("y"), Key("t")],
output_keys=[Key("u")],
cfg=cfg.arch.fully_connected,
)
speed_net = instantiate_arch(
input_keys=[Key("x"), Key("y")],
output_keys=[Key("c")],
cfg=cfg.arch.fully_connected,
)
nodes = (
we.make_nodes(detach_names=["c"])
+ ob.make_nodes(detach_names=["c"])
+ [
wave_net.make_node(name="wave_network"),
speed_net.make_node(name="speed_network"),
]
)

The PDE nodes are created using both the `OpenBoundary`

PDE defined above and the `WaveEquation`

PDE defined within Modulus. Since this is a time-dependent problem, use the `time=True`

arguments.

Here two neural network nodes are defined where `wave_net`

is used to learn the solution to the wave equation and `speed_net`

is used to learn the velocity model given the training data (`wave_speed_invar`

and `wave_speed_outvar`

) defined above.

The 2D geometry of the computational domain as well as the time range of the solution is defined:

```
```*# define geometry*
dLen = 2 *# km*
rec = Rectangle((0, 0), (dLen, dLen))
*# define sympy domain variables*
x, y, t = Symbol("x"), Symbol("y"), Symbol("t")
*# define time range*
time_length = 1
time_range = {t: (0.15, time_length)}

A total of four different constraints is sufficient to define this problem. The first constraint is a simple supervised constraint which matches the velocity model to the training data above. The second is another supervised constraint, constraining the wavefield solution to match the numerical wavefield generated using Devito at 4 starting time steps. The third constraint ensures the `WaveEquation`

PDE is honoured in the interior of the domain, whilst the fourth constraint ensures the `OpenBoundary`

PDE is honoured along the boundaries of the domain:

```
```*# make domain*
domain = Domain()
*# add velocity constraint*
velocity = PointwiseConstraint.from_numpy(
nodes=nodes, invar=wave_speed_invar, outvar=wave_speed_outvar, batch_size=1024
)
domain.add_constraint(velocity, "Velocity")
*# add initial timesteps constraints*
batch_size = 1024
for i, ms in enumerate(np.linspace(150, 300, 4)):
timestep_invar, timestep_outvar = read_wf_data(ms, dLen)
lambda_weighting = {}
lambda_weighting["u"] = np.full_like(timestep_invar["x"], 10.0 / batch_size)
timestep = PointwiseConstraint.from_numpy(
nodes,
timestep_invar,
timestep_outvar,
batch_size,
lambda_weighting=lambda_weighting,
)
domain.add_constraint(timestep, f"BC{i:04d}")
*# add interior constraint*
interior = PointwiseInteriorConstraint(
nodes=nodes,
geometry=rec,
outvar={"wave_equation": 0},
batch_size=4096,
bounds={x: (0, dLen), y: (0, dLen)},
lambda_weighting={"wave_equation": 0.0001},
parameterization=time_range,
)
domain.add_constraint(interior, "Interior")
*# add open boundary constraint*
edges = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=rec,
outvar={"open_boundary": 0},
batch_size=1024,
lambda_weighting={"open_boundary": 0.01 * time_length},
parameterization=time_range,
)
domain.add_constraint(edges, "Edges")
*# add validators*
for i, ms in enumerate(np.linspace(350, 950, 13)):
val_invar, val_true_outvar = read_wf_data(ms, dLen)
validator = PointwiseValidator(
nodes=nodes,
invar=val_invar,
true_outvar=val_true_outvar,
batch_size=1024,
plotter=WavePlotter(),
)
domain.add_validator(validator, f"VAL_{i:04d}")
validator = PointwiseValidator(
nodes=nodes,
invar=wave_speed_invar,
true_outvar=wave_speed_outvar,
batch_size=1024,
plotter=WavePlotter(),
)
domain.add_validator(validator, "Velocity")

The supervised constraints are added using the `PointwiseConstraint.from_numpy(...)`

constructor.

Finally, you can validate your results using 13 later time steps from Devito, defined within a Modulus validator as shown below.

```
```*# add validators*
for i, ms in enumerate(np.linspace(350, 950, 13)):
val_invar, val_true_outvar = read_wf_data(ms, dLen)
validator = PointwiseValidator(
nodes=nodes,
invar=val_invar,
true_outvar=val_true_outvar,
batch_size=1024,
plotter=WavePlotter(),
)
domain.add_validator(validator, f"VAL_{i:04d}")
validator = PointwiseValidator(
nodes=nodes,
invar=wave_speed_invar,
true_outvar=wave_speed_outvar,
batch_size=1024,
plotter=WavePlotter(),
)
domain.add_validator(validator, "Velocity")

You can define a custom tensorboard plotter within `examples/seismic_wave/wave_2d.py`

to add validation plots to TensorBoard, see TensorBoard in Modulus for more details on how to do this.

Now that you have defined the model, its constraints and validation data, you can form a solver to train your model. All of the hyperparameters used for the example can be found in the project’s Hydra configuration file.

```
```slv = Solver(cfg, domain)
slv.solve()

The training results can be viewed in TensorBoard. It can be seen that the Modulus results are noticeably better than Devito, predicting the wavefield with much less boundary reflections, especially at later time steps.

Fig. 62 *Comparison of Modulus results with Devito solution*