Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Output data in Spyro #118

Open
acse-yw11823 opened this issue Jul 4, 2024 · 6 comments
Open

Output data in Spyro #118

acse-yw11823 opened this issue Jul 4, 2024 · 6 comments

Comments

@acse-yw11823
Copy link

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 Spyro implement any form of automatic data scaling or normalization for the amplitude values in seismic simulations? How should I manually change the data scalling? 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))

Many thanks for your time!! Thanks for the help a lot!!!!
@Olender
Copy link
Collaborator

Olender commented Jul 4, 2024

I will review these files later to fully address your questions. However, I'm currently behind on my thesis, so it might take some time.

In the meantime, here are a few points I can already highlight:

  • For acoustic propagation problems, I strongly recommend avoiding degree=1. Ideally, use degree=4 for better runtime and reduced memory usage. Degree=1 results in significant error, even with refined meshes. For more details, please refer to this paper: GMD, 2022. Figure 8 in the paper illustrates why it is best to avoid degree=1.

  • Before completing the major refactoring in PR Major code refactoring #104, I compared results and some computational aspects with Devito. I recall having to scale the Devito results by dx^2 (with dx from Devito). As shown in this Devito tutorial, Devito also needs to do this when comparing with the forward analytical solution. I think they also do it internally when calculating gradients.

@acse-yw11823
Copy link
Author

acse-yw11823 commented Jul 4, 2024 via email

@Olender
Copy link
Collaborator

Olender commented Jul 4, 2024

Some more observations:

Starting from line 181 in the Devito file operators.py, it's clear that you also need to also consider m (which should be c^2 where the source is located) to properly scale the Devito forward results. Applying this adjustment to the latest result you showed me, 0.32277313 becomes 5.164370080000001e-06. While this is still not 6.19580433e-06, it is closer. The remaining difference is likely due to various factors such as different delays, errors related to meshing, and the degree of 1 (please go higher).

Additionally, it is important to note that our source and Devito's source probably have different starting delays for the Ricker wavelet. I recommend plotting the results from a single receiver to better understand the discrepancies, or examining where Devito defines their delay.

For degree and mesh sizes, please follow the specifications mentioned in the spyro paper.

If you are doing an in depth look in accuracy between devito and spyro I also suggest switching to the newest version of spyro.

@Olender
Copy link
Collaborator

Olender commented Jul 4, 2024

Thank you for your reply. I hope your thesis goes well. I will continue researching this issue in the meantime, and I will let you know if I figure it out. Thanks!

Thanks! It is almost finished

@acse-yw11823
Copy link
Author

acse-yw11823 commented Jul 4, 2024 via email

@acse-yw11823
Copy link
Author

Hello, I applied the adjustment to the latest result, but this is what I get

Data([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
       0.        ],
      [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
       0.        ],
      [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
       0.        ],
      ...,
      [0.00698016, 0.00686842, 0.00670907, ..., 0.00670909, 0.00686841,
       0.00698017],
      [0.00676926, 0.00664856, 0.00648524, ..., 0.00648525, 0.00664857,
       0.00676925],
      [0.00654491, 0.00641714, 0.0062522 , ..., 0.0062522 , 0.00641716,
       0.0065449 ]], dtype=float32)

I attached the code I adjusted here.

# Define a physical size
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.)  # What is the location of the top left corner. This is necessary to define
# the absolute location of the source and receivers

v = np.full(shape, 2.5, dtype=np.float32)  
center_x, center_y = shape[0] // 2, shape[1] // 2  
radius = 30 

for x in range(shape[0]):
    for y in range(shape[1]):
        if (x - center_x) ** 2 + (y - center_y) ** 2 <= radius ** 2:
            v[x, y] = 3.0


# With the velocity and model size defined, we can create the seismic model that
# encapsulates this properties. We also define the size of the absorbing layer as 10 grid points
model = Model(vp=v, origin=origin, shape=shape, spacing=spacing,
              space_order=2, nbl=40, bcs="damp")

t0 = 0.  # Simulation starts a t=0
tn = 1000.  # Simulation last 1 second (1000 ms)
dt = model.critical_dt  # Time step from model grid spacing
time_range = TimeAxis(start=t0, stop=tn, step=dt)

f0 = 0.010  # Source peak frequency is 10Hz (0.010 kHz)
src = RickerSource(name='src', grid=model.grid, f0=f0,
                   npoint=1, time_range=time_range)
# Set source coordinates
src.coordinates.data[0, 1] = np.array(model.domain_size[1]) * .5
src.coordinates.data[0, 0] = 20.  # Depth is 20

# Create symbol for 101 receivers
rec = Receiver(name='rec', grid=model.grid, npoint=101, time_range=time_range)
# Prescribe even spacing for receivers along the x-axis
rec.coordinates.data[:, 1] = np.linspace(0, model.domain_size[1], num=101)
rec.coordinates.data[:, 0] = 980.  # Depth is 20m

# In order to represent the wavefield u and the square slowness we need symbolic objects 
# corresponding to time-space-varying field (u, TimeFunction) and 
# space-varying field (m, Function)
# Define the wavefield with the size of the model and the time dimension
u = TimeFunction(name="u", grid=model.grid, time_order=2, space_order=2)

# We can now write the PDE
pde = u.dt2 - (1./model.m)*u.laplace + model.damp * u.dt * (1./model.m)

# This discrete PDE can be solved in a time-marching way updating u(t+dt) from the previous time step
# Devito as a shortcut for u(t+dt) which is u.forward. We can then rewrite the PDE as 
# a time marching updating equation known as a stencil using customized SymPy functions
stencil = Eq(u.forward, solve(pde, u.forward))
# Finally we define the source injection and receiver read function to generate the corresponding code
src_term = src.inject(field=u.forward, expr=src * dt**2 * model.m)
# Create interpolation expression for receivers
rec_term = rec.interpolate(expr=u.forward)

op = Operator([stencil] + src_term + rec_term, subs=model.spacing_map)
op(time=100, dt=model.critical_dt)
plot_image(u.data[0, :, :], cmap="seismic")

u.data[:] = 0.0
op(time=time_range.num-1, dt=model.critical_dt)
plot_shotrecord(rec.data, model, t0, tn)

Could you please tell me where my code might be incorrect, or if there is a simpler way to adjust it? Many thanks for your help!!!!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@Olender @acse-yw11823