Model times kernel is not travel time.

sstaehler commented 9 years ago

The most basic test, implemented in master_queue.f90 fails

auerl commented 9 years ago

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):

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?
Are there unrealistically large values close to the source and the receiver, which may obscure the integrated kernel values? Looking at the K_x array it seems that there are some elements with extremely large values, which may completely dominate the summation with the background model.
Do the kernels look the same for more "complicated" (i.e. ones where the source is not at the northpole) source-receiver geometries?
Do the kernels visually look similar to published kernels?
Do the kernel computations work for all possible phases?

kasra-hosseini commented 9 years ago

Just some quick thoughts:

Are the kernel values (after dividing by the element volumes) independent of which mesh is used?

I think so! at least by assuming a very fine initial mesh, it should be possible to project it onto other inversion grids. I am saying very fine because otherwise we have not sampled the domain correctly.

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

I assume that the best is to check the integral of that (as we have already talked about) and compare it with the travel time...but I have also seen the same values. They actually report it like that to be independent of the inversion grid (if I am not mistaken).

Are the kernel values independent of the source magnitude used in the fwd simulation?

Yes! This is why it is divided by the length of Vp (last term in eq 5.13)

martinvandriel commented 9 years ago

@auerl:

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.

@kasra-hosseini:

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.

kasra-hosseini commented 9 years ago

@martinvandriel You are right, I was just thinking in this way:

\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?

martinvandriel commented 9 years ago

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.

kasra-hosseini commented 9 years ago

Agree! I remember that derivation from Warnemuende...that would be indeed the best.

sstaehler commented 9 years ago

@auerl

Are there unrealistically large values close to the source and the receiver, which may obscure the integrated kernel values? Looking at the K_x array it seems that there are some elements with extremely large values, which may completely dominate the summation with the background model.

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.

@martinvandriel

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.

I agree. The code can calculate volumetric kernels on any mesh, so we might also do that with tetrahedrons. But either should work.

sstaehler commented 9 years ago

@auerl I might as well answer all your questions. The answers are about the current state of the code (3c937554cc34029d73c43297add95baf6b09ab17):

Are the kernel values (after dividing by the element volumes) independent of which mesh is used?

Difficult to say. Should test that

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

See Martin's answer on Kx vs Aij. The only plots I know are of Kx.

Are the kernel values independent of the source magnitude used in the fwd simulation?

They should be and as I see it, they are.

Do the kernels look the same for more "complicated" (i.e. ones where the source is not at the northpole) source-receiver geometries?

Yes, that's fine.

Do the kernels visually look similar to published kernels?

Yes, they do

Do the kernel computations work for all possible phases?

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 $\int K_X(t) dt$. Could be tested with the Parseval-theorem based integration in frequency domain.

So to summarize @auerl, @kasra-hosseini : Tasks for the inclined would be:

[ ] Test more complicated phases
[ ] Test different meshes. To remove the multiplication by volume, set the volume term in l. 367 in slave.f90 https://github.com/sstaehler/kerner/blob/master/slave.f90#L367 to 1:

    call int_kernel(ibasisfunc)%initialize_montecarlo(parameters%nkernel, &
    1, &
    parameters%allowed_error, &
    parameters%allowed_relative_error)

martinvandriel commented 9 years ago

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.

auerl commented 9 years ago

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.

kernel_5deg_lugano

kernel_5deg_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:

kernel_1 25deg_lugano

kernel_tet_lugano

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:

velocity_integration

@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.

kasra-hosseini commented 9 years ago

@auerl Sure, I will try to do it this weekend. Is that OK?

@martinvandriel I am sure I saw it somewhere in Nolet's book, but unfortunately I do not have the book now. I just wrote it down and you can find the scan version here...It is for ray-based tomography.

Kernel

martinvandriel commented 9 years ago

@kasra-hosseini In the first equation, you only integrate the kernel itself, not the kernel times the background model? Where is the K in the following?

tnissen commented 9 years ago

Aren't the attached couple of lines based on born scattering sufficient, at least intuitively?

-------- Original message -------- From: Martin van Driel notifications@github.com Date: 12/03/2015 13:53 (GMT+00:00) To: sstaehler/kerner kerner@noreply.github.com Subject: Re: [kerner] Model times kernel is not travel time. (#21)

@kasra-hosseini In the first equation, you only integrate the kernel itself, not the kernel times the background model? Where is the K in the following?

— Reply to this email directly or view it on GitHub.

martinvandriel commented 9 years ago

@tnissen attachment?

tnissen commented 9 years ago

Attachment didn't seem to work... Here: http://seis.earth.ox.ac.uk/20150312_141808.jpg

On 12/03/2015 14:23, Martin van Driel wrote:

@tnissen https://github.com/tnissen attachment?

— Reply to this email directly or view it on GitHub https://github.com/sstaehler/kerner/issues/21#issuecomment-78488659.

Tarje

<>--<>--<>--<>--<>--<> Dept. of Earth Sciences Oxford University South Parks Road Oxford OX1 3AN; UK tel: +44 1865 282149 fax: +44 1865 272072 web: seis.earth.ox.ac.uk http://seis.earth.ox.ac.uk <>--<>--<>--<>--<>--<>

kasra-hosseini commented 9 years ago

@martinvandriel Well the very first equation (at top) that I wrote was really only defining the question, not carefully written. However, about your second question: I think the whole kernel in ray theory is:

-\int_{ray}{1./V}ds

Not sure actually...

martinvandriel commented 9 years ago

@tnissen well now we have a - b = c - d. That is obviously true if a = c and b = d, but we need to show that this is the only solution.

kasra-hosseini commented 9 years ago

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.

martinvandriel commented 9 years ago

I think you are now confusing the integral over time and volume. The one over volume should be the relevant one here.

kasra-hosseini commented 9 years ago

You are right! similar to -\int_{ray}{1./V}ds = -T in ray theory, the relevant one should be the spatial integration.

martinvandriel commented 9 years ago

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...

tnissen commented 9 years ago

I would say Born scattering validates that... ?

On 12/03/2015 14:36, Martin van Driel wrote:

@tnissen https://github.com/tnissen well now we have a - b = c - d. That is obviously true it a = c and b = d, but we need to show that this is the only solution.

— Reply to this email directly or view it on GitHub https://github.com/sstaehler/kerner/issues/21#issuecomment-78491468.

Tarje

<>--<>--<>--<>--<>--<> Dept. of Earth Sciences Oxford University South Parks Road Oxford OX1 3AN; UK tel: +44 1865 282149 fax: +44 1865 272072 web: seis.earth.ox.ac.uk http://seis.earth.ox.ac.uk <>--<>--<>--<>--<>--<>

martinvandriel commented 9 years ago

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?

auerl commented 9 years ago

... 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]).

eqtarje2

eqtarje1

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].

trompeq

Does somebody have an idea?

tnissen commented 9 years ago

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 -------- From: auerl notifications@github.com Date: 13/03/2015 12:26 (GMT+00:00) To: sstaehler/kerner kerner@noreply.github.com Cc: tnissen tarje@alumni.princeton.edu Subject: Re: [kerner] Model times kernel is not travel time. (#21)

... 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].

[Uploading trompeq.png . . .]()

Does somebody have an idea?

— Reply to this email directly or view it on GitHub.

auerl commented 9 years ago

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

auerl commented 9 years ago

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.

tnissen commented 9 years ago

Yes, that's how I had it in the original kernel code as well - it was indeed related to the units of the backward wavefield.

The 7 orders of magnitude w.r.t. Dahlen... could that have to do with Mij units between cgs to SI? My recollection of the Dahlen 2000 paper is that the kernel represents an explosion (isotropic radiation), but I don't recall the incorporation of the moment. Is the backward wavefield for a delta function of the correct integral equal 1? At least in an earlier version we multiplied it with a large number to avoid having all too small numbers in the wavefield. I bet you take care of all this, just checking.

On 17/03/2015 15:25, auerl wrote:

I went through the units of our implementation and everything seems consistent.That 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). Having the units of all the involved quantities, I'll next check whether they are all approximately on the right order.

image1 https://cloud.githubusercontent.com/assets/5452469/6690557/84e18a72-ccc1-11e4-852b-ff3fadc6240f.jpeg

image2 https://cloud.githubusercontent.com/assets/5452469/6690559/890e7dda-ccc1-11e4-9348-229032c72631.jpeg

— Reply to this email directly or view it on GitHub https://github.com/sstaehler/kerner/issues/21#issuecomment-82411358.

Tarje

<>--<>--<>--<>--<>--<> Dept. of Earth Sciences Oxford University South Parks Road Oxford OX1 3AN; UK tel: +44 1865 282149 fax: +44 1865 272072 web: seis.earth.ox.ac.uk http://seis.earth.ox.ac.uk <>--<>--<>--<>--<>--<>

sstaehler commented 9 years ago

The scalar moment should not affect the kernel value, since the convolved wavefield is divided by the squared seismogram in 5.13. I'll run a test to see whether this is the case in our code.

On the amplitude of the backward field: Good question, I'll have a look

sstaehler commented 9 years ago

I close this discussion now. The underlying assumption seems to have been wrong.

seismology / mc_kernel

Model times kernel is not travel time. #21