2D Seismic Wave Propagation¶

Introduction¶

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:

1. How to solve a 2d time-dependent problem in Modulus

2. How to define an open boundary condition using custom equations

3. How to present a variable velocity model as an additional network

Note

This tutorial assumes that you have completed the Lid Driven Cavity Background 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.

Problem Description¶

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:

(130)$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. 42). 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. 43, with a peak wavelet frequency of $$f_0=15 Hz$$.

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$$).

Problem Setup¶

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

Note

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

Defining Boundary Conditions¶

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, apply an open boundary condition, also called a radiation condition, which imposes the first-order PDEs at the boundaries:

(131)$\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(PDES):
"""
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)
)

Variable Velocity Model¶

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$$:

Note

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

Creating PDE and Neural Network Nodes¶

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

# define PDEs
we = WaveEquation(u="u", c="c", dim=2, time=True)
ob = OpenBoundary(u="u", c="c", dim=2, time=True)

# 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", jit=cfg.jit),
speed_net.make_node(name="speed_network", jit=cfg.jit),
]
)

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.

Creating Geometry¶

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)}

Add four different constraints to your domain for 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 timesteps. 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()

velocity = PointwiseConstraint.from_numpy(
nodes=nodes, invar=wave_speed_invar, outvar=wave_speed_outvar, batch_size=1024
)

batch_size = 1024
for i, ms in enumerate(np.linspace(150, 300, 4)):
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,
)

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},
param_ranges=time_range,
)

edges = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=rec,
outvar={"open_boundary": 0},
batch_size=1024,
lambda_weighting={"open_boundary": 0.01 * time_length},
param_ranges=time_range,
)

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

Validation¶

Finally, validate your results using 13 later timesteps from Devito, defined within a Modulus validator as shown below.

for i, ms in enumerate(np.linspace(350, 950, 13)):
validator = PointwiseValidator(
val_invar, val_true_outvar, nodes, batch_size=1024, plotter=WavePlotter()
)
validator = PointwiseValidator(
wave_speed_invar,
wave_speed_outvar,
nodes,
batch_size=1024,
plotter=WavePlotter(),
)