-
Notifications
You must be signed in to change notification settings - Fork 8
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.
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
Model times kernel is not travel time. #21
Comments
A basic question regarding the units of the kernels we compute: Our kernels are essentially given by eqn. 5.13 in Tarjes thesis, right? So their unit should be s^2 / m . When multiplied with the background model what should come out is the absolute traveltime. But is this consistent with eqn. 5.12 in Tarjes thesis which relates differential traveltimes with (percentual) model perturbations via (relative?) kernels. I probably don't quite understand the difference between absolute and relative kernels, can sb explain? Here are some more basic things that can/should probably be checked to make sure that our kernels are plausible (not all directly related to this issue reported by Simon):
|
Just some quick thoughts: Are the kernel values (after dividing by the element volumes) independent of which mesh is used?
In which range do the kernel values lie? Are they similar to published kernels? In the literature the kernels are usually plotted on scales b/w -5 to 5 * 10-6 s/km3
Are the kernel values independent of the source magnitude used in the fwd simulation?
|
I think in 5.12 and 5.13 the perturbations are absolute, as this comes from 5.3. To understand the integral with the model, it is maybe useful to think about the analogy with rays. However, I have not seen a rigorous proof and I could not make it up myself in 5 minutes. When you come across it, please let me know. Kernel values: I guess the difference is, whether we talk about the kernel as a function of space K(x) or the values Aij in the matrix. The first ones should be independent of the mesh, the latter depend on the normalization chosen such that it can be multiplied with the modelvector. |
@martinvandriel \int_{Volume}KdV ---> that should be in second (arrival time) which is A (as you said) in dT = Ac If we only calculate K, it should be then sec/km3; However, we dont do this. So if we go for fine inversion grid and then divide it by volume of each (sort of similar to numerical integration), we should be able to say it is actually K (ad hoc), and can be projected on any other inversion grid (which is relatively larger than the initial one). Does that make sense? |
I think to just look at the kernel values it would be easier to just use an orthonormal basis like Ludwigs voxels, because in this case the relation between the kernel and the matrix is trivial, see the draft of the kernel paper, eq 7. We could then just plot the kernels with the tools we have and look at the values. |
Agree! I remember that derivation from Warnemuende...that would be indeed the best. |
Can you quantify that? The values at source and receiver should be considerably larger, since all energy has to go through these cells, which gives them a large influence on the result.
I agree. The code can calculate volumetric kernels on any mesh, so we might also do that with tetrahedrons. But either should work. |
@auerl
Difficult to say. Should test that
See Martin's answer on Kx vs Aij. The only plots I know are of Kx.
They should be and as I see it, they are.
Yes, that's fine.
Yes, they do
I tested it a bit, but that could be done more rigorously. It did not really work well for more complicated phases with the time-domain integration for So to summarize @auerl, @kasra-hosseini :
call int_kernel(ibasisfunc)%initialize_montecarlo(parameters%nkernel, &
1, &
parameters%allowed_error, &
parameters%allowed_relative_error) |
We (Ludwig, Dr. Boehm, me) tried to write down why the kernel times the model should be absolute traveltime and I think we did not quite get there. Does anyone of you have notes about it? If not, do you know someone who has them and we could ask? Without understanding this properly, there is IMHO no point in using this as a test. |
Hi all, I just compared actual kernel values computed at cscs and our machine in Zurich, and surprisingly there is quite a difference. See the figures below. No idea what this could be. The compiler is gfortran in both cases (4.8.1 in lugano and 4.8.2 in zurich), wavefields are the same, version of the code is the same, inparam file is the same, netcdf version is 4.3.1 in lugano and 4.1.3 in zurich. Maybe an issue related to using different (incompatible) versions of netcdf to write / read the netcdf wavefields. When running the code in Zurich also the projected background model velocities are 6 (!) orders too large. From what I see, the kernel values (if not multiplied with the volume of the grid cells) does not seem to be influenced by the chosen mesh), which is good: Another thing that I noticed is that at places where the reference model has jumps (CMB, surface, probably also the internal discontinuities) the default MC-integration parameters that are set for projecting the model onto the inversion grid don't seem to ensure convergence, since there is lateral variations. See the following screenshot: @ALL: I think, what would be extremely helpful would be to have a running version of the Princeton's group (Dahlen) kernel code on our machine in Zurich, to which everybody has access. Would allow us to output values at various locations for similar kernels, and help to get a feeling for how large they need to be. Also, it would be good to see how exactly the test of multiplying kernels with the background model to get absolute traveltimes is implemented. @kasra-hosseini: Could you install such a test-version for us? I think you are the only one of us who knows the code. |
@auerl @martinvandriel |
@kasra-hosseini |
Aren't the attached couple of lines based on born scattering sufficient, at least intuitively? -------- Original message -------- @kasra-hosseini — |
@tnissen attachment? |
Attachment didn't seem to work... Here: On 12/03/2015 14:23, Martin van Driel wrote:
Tarje <>--<>--<>--<>--<>--<> |
@martinvandriel -\int_{ray}{1./V}ds Not sure actually... |
@tnissen |
Sorry, I meant: -1/V is the kernel (ray theory)...without integral. Just some thoughts: In all sensitivity kernels: K_rho, K_lambda and K_mu we have: -\int_{0}^{t} something dt ---> that something will be normalized and if we consider it is constant (just approximate), we will end up with: -\int_{0}^{t}dt which is -T (travel time). Does that make sense? Or I am making myself more and more confused :D. |
I think you are now confusing the integral over time and volume. The one over volume should be the relevant one here. |
You are right! similar to -\int_{ray}{1./V}ds = -T in ray theory, the relevant one should be the spatial integration. |
My line of thought at the moment is to use travelime Kernels with respect to slowness instead of velocity or moduli. For these it is straightforward to show the property we are looking for. The missing step is now to relate these to the other kernels in first order. I'll try to write it down properly... |
I would say Born scattering validates that... ? On 12/03/2015 14:36, Martin van Driel wrote:
Tarje <>--<>--<>--<>--<>--<> |
The more I think about it, the less clear it becomes: what is T in a finite frequency sense? delta T is clear, as it comes from comparing two signals, but what is the absolute traveltime T? |
... another thing we are struggling with right now concerns the units of our kernels: When looking at eqns. 5.12 and 5.13 in Tarjes thesis (on which our implementation is based on), there doesn't seem to be correspondence between the left and the right side in 5.12 in terms of the physical units. Can somebody explain what I am missing here? In 5.12 \delta T is given in [s] and \delta rho is given in [kg/m^3] so the kernels as defined in 5.13 ought to have the unit [s/kg]. However, if I see it correctly, 5.13 has the unit [s^2/m] (since vp is the velocity seismogram, the K inside the time integral is in [s] and the normalization factor is in [s/m^2]. Usually the kernels in all kernel galleries are plotted in the unit [s/m^3](or [km^3]). I have the same problem in the banana donought section of Tromps paper on "Seismic tomography, adjoint methods, time reversal and banana-donought kernels" (2004, GJI). The following equations don't seem to make sense. In eqn. 46 [s] is equated with [kgm^2/s]. Does somebody have an idea? |
I will be in Zurich next week (wed afternoon/thur morning), let's sit together on this, i can't beforehand in more detail. When writing my thesis I spent quite some time exactly on this...as I recall (almost 10 years ago!) it is related to the fact that the backward/Adjoint wavefield has weird units... definitely not those of a seismic wavefield -------- Original message -------- ... another thing we are struggling with right now concerns the units of our kernels: When looking at eqns. 5.12 and 5.13 in Tarjes thesis (on which our implementation is based on), there doesn't seem to be correspondence between the left and the right side in 5.12 in terms of the physical units. Can somebody explain what I am missing here? In 5.12 \delta T is given in [s] and \delta rho is given in [kg/m^3] so the kernels as defined in 5.13 ought to have the unit [s/kg]. However, if I see it correctly, 5.13 has the unit [s^2/m] (since vp is the velocity seismogram, the K inside the time integral is in [s] and the normalization factor is in [s/m^2]. Usually the kernels in all kernel galleries are plotted in the unit s/m^3. I have the same problem in the banana donought section of Tromps paper on "Seismic tomography, adjoint methods, time reversal and banana-donought kernels" (2004, GJI). The following equations don't seem to make sense. In eqn. 46 [s] is equated with [kgm^2/s]. Does somebody have an idea? — |
Hi all, It seems that what I missed before was that the delta function in the force term has a unit. See Carl Tapes thesis at page 184: "In practice, the adjoint sources have a m^−3 or m^−2 quantity as well, due to the 3D or 2D volume for the delta function, δ(x), that is applied at the source location" or wikipedia "Eine direkte Folgerung aus der Skalierungseigenschaft ist die Dimension bzw. Maßeinheit der Delta-Distribution. Sie entspricht genau der reziproken Dimension ihres Arguments". Accounting for the additional 1/m^3 the units in the equations in Tromps paper make sense and one gets the [s^2/kg m] = [1/N] as the unit of the adjoint wavefield, if I see it correctly. Greets, Ludwig |
I went through the units of our implementation and everything seems consistent. The kernels, as we have them right now, are for absolute perturbations, and not relative ones. For absolute ones we just need to add an extra multiplication with the background model (see 2nd page). Should we introduce an option in the inparam file for that? Relative ones are probably more important. Having the units of all the involved quantities, I'll now check whether they are all approximately on the right order. |
Yes, that's how I had it in the original kernel code as well - it was The 7 orders of magnitude w.r.t. Dahlen... could that have to do with On 17/03/2015 15:25, auerl wrote:
Tarje <>--<>--<>--<>--<>--<> |
The scalar moment should not affect the kernel value, since the convolved wavefield is divided by the squared seismogram in 5.13. On the amplitude of the backward field: Good question, I'll have a look |
I close this discussion now. The underlying assumption seems to have been wrong. |
- Takes some time for very large meshes and does not seem to make much sense anyway (see #21)
The most basic test, implemented in master_queue.f90 fails
The text was updated successfully, but these errors were encountered: