Data type in Devito

[acse-yw11823]

I concurrently use Spyro and Devito tools to process the most simple velocity model. In Spyro, I obtained minimal data values in true_d.data, as shown below:

array([[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       ...,
       [6.05248748e-06, 6.05709782e-06, 5.79909842e-06, ...,
        5.79183276e-06, 6.12123652e-06, 6.19926250e-06],
       [6.06644009e-06, 6.01030630e-06, 5.70104761e-06, ...,
        5.68715204e-06, 6.06524176e-06, 6.20124016e-06],
       [6.07264028e-06, 5.95774122e-06, 5.59937139e-06, ...,
        5.57923501e-06, 6.00384318e-06, 6.19580433e-06]])

Conversely, in Devito, the amplitude values are significantly larger, as illustrated below:

Data([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
       0.        ],
      [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
       0.        ],
      [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
       0.        ],
      ...,
      [0.32997185, 0.32673766, 0.32016267, ..., 0.31328549, 0.31975521,
       0.32302558],
      [0.32989321, 0.32637454, 0.31955268, ..., 0.3127043 , 0.31940411,
       0.32293749],
      [0.32973529, 0.32594303, 0.31887859, ..., 0.31206041, 0.31898795,
       0.32277313]])

How should I interpret these differences? Does Devito implement any form of automatic data scaling or normalization for the amplitude values in seismic simulations? I will also attach the relevant sections of code from both Devito and Spyro for your reference.

For devito:

# Define true and initial model
shape = (101, 101)  # Number of grid point (nx, nz)
spacing = (10., 10.)  # Grid spacing in m. The domain size is now 1km by 1km
origin = (0., 0.)  # Need origin to define relative source and receiver locations

#Set the number of boundary layers
nbl = 40

#True Vp model
model = demo_model('circle-isotropic', vp_circle=3.0, vp_background=2.5,
                    origin=origin, shape=shape, spacing=spacing, nbl=nbl, dt=0.5, dtype=np.float64)

#Smooth/Initial Vp model
model0 = demo_model('circle-isotropic', vp_circle=2.5, vp_background=2.5,  
                     origin=origin, shape=shape, spacing=spacing, nbl=nbl,
                    grid = model.grid, dt=0.5, dtype=np.float64)

# Define acquisition geometry: 
# nshots = 9  # Number of shots to used to generate the gradient
nreceivers = 101  # Number of receiver locations per shot 
# fwi_iterations = 5  # Number of outer FWI iterations

# Note, distances below correspond to meters.
# Position the source:
src_coordinates = np.empty((1, 2))
src_coordinates[0, 1] = np.array(model.domain_size[1]) * .5
src_coordinates[0, 0] = 20.

# Initialize receivers for synthetic and imaging data
rec_coordinates = np.empty((nreceivers, 2))
rec_coordinates[:, 1] = np.linspace(0, model.domain_size[0], num=nreceivers)
rec_coordinates[:, 0] = 980.

# Create the Geometry
# Along with storing the source and receiver locations defined above this will work out the points in time at
# at which we'll need the source signal and compute them according to the definiction of the Ricker-wavelet
t0 = 0.
tn = 1000. # 1000ms
f0 = 0.010 # 10Hz (=0.01 kHz)
geometry = AcquisitionGeometry(model, rec_coordinates, src_coordinates, t0, tn, f0=f0, src_type='Ricker')

solver = AcousticWaveSolver(model, geometry, space_order=4)

# Compute synthetic data with forward operator and calculate the total running time
t1=time.time()
true_d,a, b = solver.forward(vp=model.vp)
end_time = time.time()
running_time = end_time - t1
# Print the running time
print("Running Time - devito: {:.2f} seconds".format(running_time))

# Compute initial data with forward operator 
smooth_d, _, _ = solver.forward(vp=model0.vp)

For spyro:

sou_pos = [(0.02, 0.5)]
rec_pos = spyro.create_transect((0.98, 0.0), (0.98, 1.0), 101)

model = {}
model["opts"] = {
    "method": "KMV",  # either CG or KMV
    "quadrature": "KMV",  # Equi or KMV
    "degree": 1,  # p order
    "dimension": 2,  # dimension
    # "regularization": False,  # regularization is on?
    # "gamma": 1.0,  # regularization parameter
}
model["parallelism"] = {
    "type": "spatial", 
}
model["mesh"] = {
    "Lz": 1.0,  # depth in km - always positive
    "Lx": 1.0,  # width in km - always positive
    "Ly": 0.0,  # thickness in km - always positive
    "meshfile": "not_used.msh", 
    "initmodel": "not_used.hdf5", 
    "truemodel": "not_used.hdf5", 
}
model["BCs"] = {
    "status": True,  # True or false
    "outer_bc": "non-reflective",  # None or non-reflective (outer boundary condition)
    "damping_type": "polynomial",  # polynomial, hyperbolic, shifted_hyperbolic
    "exponent": 2,  # damping layer has a exponent variation
    "cmax": 4.5,  # maximum acoustic wave velocity in PML - km/s
    "R": 1e-6,  # theoretical reflection coefficient
    "lz": 0.4,  # thickness of the PML in the z-direction (km) - always positive
    "lx": 0.4,  # thickness of the PML in the x-direction (km) - always positive
    "ly": 0.0,  # thickness of the PML in the y-direction (km) - always positive
}
model["acquisition"] = {
    "source_type": "Ricker", 
    "source_pos": sou_pos, 
    "frequency": 10.0, 
    # "delay": 0.1, 
    "delay": 0.1, 
    "receiver_locations": rec_pos, 
}
model["timeaxis"] = {
    "t0": 0.0,  # Initial time for event
    "tf": 1.00,  # Final time for event
    "dt": 0.0005, 
    "amplitude": 1,  # the Ricker has an amplitude of 1.
    "nspool": 100,  # how frequently to output solution to pvds
    "fspool": 100,  # how frequently to save solution to RAM
    "skip": 4,
}

# mesh = UnitSquareMesh(100, 100, 1.0, 1.0)
mesh = RectangleMesh(101, 101, 1.0, 1.0)
# mesh.coordinates.dat.data[:, 0] -= 0.25 将x坐标向左平移0.25个单位
# mesh.coordinates.dat.data[:, 1] -= 1.25 将y坐标向下平移1.25个单位

comm = spyro.utils.mpi_init(model)
element = spyro.domains.space.FE_method(mesh, model["opts"]["method"], model["opts"]["degree"])
V = FunctionSpace(mesh, element)

# Create a simple two-layer seismic velocity model `vp`.
x, y = SpatialCoordinate(mesh)
radius = 0.3
center_x, center_y = 0.5, 0.5
condition = (x - center_x)**2 + (y - center_y)**2 < radius**2
velocity = conditional(condition, 3.0, 2.5)
vp = Function(V, name="velocity").interpolate(velocity)
File("true_velocity_model_circle.pvd").write(vp)

sources = spyro.Sources(model, mesh, V, comm)
receivers = spyro.Receivers(model, mesh, V, comm)

d0=model["timeaxis"]["t0"]
dt=model["timeaxis"]["dt"]
tf=model["timeaxis"]["tf"]
freq=model["acquisition"]["frequency"]
wavelet = spyro.full_ricker_wavelet( dt= dt, tf = tf, freq= freq,)

# Calculate running time
start_time = time.time()
p_field, p_at_recv = spyro.solvers.forward(model, mesh, comm, vp, sources, wavelet, receivers)
end_time = time.time()
running_time = end_time - start_time
print("Running Time: {:.2f} seconds".format(running_time))

[mloubout]

Well I don't know what spyro does for the source so I'm not sure I can tell you the difference. But I can already see that you use an odd set of units with spyro (time in seconds, but velocity in m/ms ) so the scaling will change.

I would probably look at normalized data to avoid trying to infer the differences in choices you made inside Spyro with Devito. THe other option is to make the acoustic solver by hand with devito so you can choose your source injection expression how you want it.

[acse-yw11823]

Hello, I noticed that the difference might be due to the fact that the parameter m we are considering is different between Devito and Spyro. In Devito, the PDE related to m (which should be c^2where the source is located) in the operator is:

pde = model.m * u.dt2 - u.laplace + model.damp * u.dt 

However, in Spyro, the PDE is:

pde = u.dt2 - (1./model.m) * u.laplace + model.damp * u.dt

Is there an easy way to adjust this in Devito, for example, when I call the simple circle model, so I can compare the results between Devito and Spyro? Many thanks for the help!!

[mloubout]

So two things:

• To be clear, these are example built with devito with some choices of parameterizations. If you want to implement different PDEs, you can use them as a skeleton, but you will need to implement them yourself,
• These two are still equivalent given proper source scaling, so it really depends what you did in spyro (you also didn't divide damp for some reason)

[acse-yw11823]

Thanks for the relpy. I am having the look now.