Open acse-yw11823 opened 2 months ago
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 #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.
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!
Sent from Outlook for iOShttps://aka.ms/o0ukef
From: Alexandre Olender @.> Sent: Thursday, July 4, 2024 3:31:33 PM To: NDF-Poli-USP/spyro @.> Cc: Wu, Yifan @.>; Author @.> Subject: Re: [NDF-Poli-USP/spyro] Output data in Spyro (Issue #118)
This email from @.*** originates from outside Imperial. Do not click on links and attachments unless you recognise the sender. If you trust the sender, add them to your safe senders listhttps://spam.ic.ac.uk/SpamConsole/Senders.aspx to disable email stamping for this address.
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, 2022https://gmd.copernicus.org/articles/15/8639/2022/. Figure 8 in the paper illustrates why it is best to avoid degree=1.
Before completing the major refactoring in PR #104https://github.com/NDF-Poli-USP/spyro/pull/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 tutorialhttps://slimgroup.github.io/Devito-Examples/tutorials/accuracy/, Devito also needs to do this when comparing with the forward analytical solution. I think they also do it internally when calculating gradients.
— Reply to this email directly, view it on GitHubhttps://github.com/NDF-Poli-USP/spyro/issues/118#issuecomment-2209134227, or unsubscribehttps://github.com/notifications/unsubscribe-auth/BDENNF6PIKDTCXGQIHVVSO3ZKVMELAVCNFSM6AAAAABKLFAUY2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEMBZGEZTIMRSG4. You are receiving this because you authored the thread.Message ID: @.***>
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.
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
Yes! I also notice the PDEs they are solving is a bit different ( putting the c^2 in different places)! I'm having the look now.
Ok. Will have a look on the delay and other parameters! Thanks!
From: Alexandre Olender @.> Sent: Thursday, July 4, 2024 3:47 PM To: NDF-Poli-USP/spyro @.> Cc: Wu, Yifan @.>; Author @.> Subject: Re: [NDF-Poli-USP/spyro] Output data in Spyro (Issue #118)
This email from @.*** originates from outside Imperial. Do not click on links and attachments unless you recognise the sender. If you trust the sender, add them to your safe senders listhttps://spam.ic.ac.uk/SpamConsole/Senders.aspx to disable email stamping for this address.
Some more observations:
Starting from line 181 in the Devito file operators.pyhttps://github.com/devitocodes/devito/blob/master/examples/seismic/acoustic/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.
— Reply to this email directly, view it on GitHubhttps://github.com/NDF-Poli-USP/spyro/issues/118#issuecomment-2209163386, or unsubscribehttps://github.com/notifications/unsubscribe-auth/BDENNF35UH2DI2GYI5WQPZLZKVOA5AVCNFSM6AAAAABKLFAUY2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEMBZGE3DGMZYGY. You are receiving this because you authored the thread.Message ID: @.***>
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!!!!
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:
Conversely, in Devito, the amplitude values are significantly larger, as illustrated below:
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:
For spyro:
Many thanks for your time!! Thanks for the help a lot!!!!