Skip to content

Commit

Permalink
examples: Update reduced tti
Browse files Browse the repository at this point in the history
  • Loading branch information
georgebisbas committed Jul 15, 2024
1 parent df2811b commit 140637a
Show file tree
Hide file tree
Showing 2 changed files with 0 additions and 67 deletions.
1 change: 0 additions & 1 deletion devito/ir/xdsl_iet/cluster_to_ssa.py
Original file line number Diff line number Diff line change
Expand Up @@ -189,7 +189,6 @@ def _visit_math_nodes(self, dim: SteppingDimension, node: Expr,
if isinstance(node, Indexed):
# If we have a time function, we compute its time offset
if isinstance(node.function, TimeFunction):
import pdb;pdb.set_trace()
time_offset = (node.indices[dim] - dim) % node.function.time_size
elif isinstance(node.function, (Function, PointSource)):
# import pdb;pdb.set_trace()
Expand Down
66 changes: 0 additions & 66 deletions tests/test_xdsl_passes.py
Original file line number Diff line number Diff line change
Expand Up @@ -211,69 +211,3 @@ def test_unary(shape, steps):
orig_data = u.data_with_halo.copy()

assert np.isclose(xdsl_data, orig_data, rtol=1e-06).all()


def test_xdsl_sine():
import numpy as np
from devito import (Function, TimeFunction, cos, sin, solve,
Eq, Operator)
from examples.seismic import TimeAxis, RickerSource, Receiver, demo_model
from matplotlib import pyplot as plt

# We will start with the definitions of the grid and the physical parameters $v_{p}, \theta, \epsilon, \delta$. For simplicity, we won't use any absorbing boundary conditions to avoid reflections of outgoing waves at the boundaries of the computational domain, but we will have boundary conditions (zero Dirichlet) at $x=0,nx$ and $z=0,nz$ for the solution of the Poisson equation. We use a homogeneous model. The model is discretized with a grid of $101 \times 101$ and spacing of 10 m. The $v_{p}, \epsilon, \delta$ and $\theta$ parameters of this model are 3600 m∕s, 0.23, 0.17, and 45°, respectively.

# NBVAL_IGNORE_OUTPUT

shape = (101, 101) # 101x101 grid
spacing = (10., 10.) # spacing of 10 meters
nbl = 0 # number of pad points

model = demo_model('layers-tti', spacing=spacing, space_order=8,
shape=shape, nbl=nbl, nlayers=1)

# initialize Thomsem parameters to those used in Mu et al., (2020)
model.update('vp', np.ones(shape)*3.6) # km/s
model.update('epsilon', np.ones(shape)*0.23)
model.update('delta', np.ones(shape)*0.17)
model.update('theta', np.ones(shape)*(45.*(np.pi/180.))) # radians

# In cell below, symbols used in the PDE definition are obtained from the `model` object. Note that trigonometric functions proper of Devito are exploited.

# %%
# Get symbols from model
theta = model.theta
delta = model.delta
epsilon = model.epsilon
m = model.m

# Use trigonometric functions from Devito
costheta = cos(theta)

# Values used to compute the time sampling
epsilonmax = np.max(np.abs(epsilon.data[:]))
deltamax = np.max(np.abs(delta.data[:]))
etamax = max(epsilonmax, deltamax)
vmax = model._max_vp
max_cos_sin = np.amax(np.abs(np.cos(theta.data[:]) - np.sin(theta.data[:])))
dvalue = min(spacing)

# Compute the dt and set time range
t0 = 0. # Simulation time start
tn = 150. # Simulation time end (0.15 second = 150 msec)
dt = (dvalue/(np.pi*vmax))*np.sqrt(1/(1+etamax*(max_cos_sin)**2)) # eq. above (cell 3)
time_range = TimeAxis(start=t0, stop=tn, step=dt)
print("time_range; ", time_range)

# time stepping
p = TimeFunction(name="p", grid=model.grid, time_order=2, space_order=2)
q = Function(name="q", grid=model.grid, space_order=8)

# Main equations
term1_p = (1 + 2*delta*(costheta**2) + 2*epsilon*costheta**4)*q.dx4

# Create stencil and boundary condition expressions
x, z = model.grid.dimensions

optime = Operator([term1_p], opt='xdsl')
optime.apply(time_M=100, dt=dt)
import pdb;pdb.set_trace()

0 comments on commit 140637a

Please sign in to comment.