2D Seismic Wave Propagation

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 Sym

  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 Introductory Example tutorial on Lid Driven Cavity flow and have familiarized yourself with the basics of Modulus Sym. Also, see the 1D Wave Equation tutorial for more information on defining new differential equations, and solving time-dependent problems in Modulus Sym.

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

velocity_model.png

Fig. 60 Velocity model of the medium

ricker_source.png

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 Sym 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.

Note

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:

Copy
Copied!
            

# 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\):

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.

Copy
Copied!
            

# 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 Sym PDE and NN nodes for this problem:

Copy
Copied!
            

# 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"), 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 Sym. 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:

Copy
Copied!
            

# define geometry dLen = 2 # km rec = Rectangle((0, 0), (dLen, dLen))

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:

Copy
Copied!
            

# 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")

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 Sym validator as shown below.

Copy
Copied!
            

# 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 Sym 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.

Copy
Copied!
            

plotter=WavePlotter(),

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

simnet_vs_devito_combined.PNG

Fig. 62 Comparison of Modulus Sym results with Devito solution

Previous 1D Wave Equation
Next Coupled Spring Mass ODE System
© Copyright 2023, NVIDIA Modulus Team. Last updated on Jan 25, 2024.