Kernel times 3D model should give dT measurements from SpecFEM

[sstaehler]

Should be done after we are sure that the kernels are similar to Dahlen kernels

[auerl]

What would be your preference as to how I implement this test? Here are three options:

1.) Add a "lateral_heterogeneity" module to the kernel code, that - similar to the equiv. AxiSEM module - reads point clouds, creates an kdtree and returns values via NN or inverse distance weighting interpolation. Then carry out MC integration like for the background model integration.

2.) Add the s362wmani and s40rts evaluation tools to the kernel code, such that one can get model values for the MC integration directly from the original model files without an initial kdtree/interpolation step.

3.) Do everything outside of the kernel code and, instead of MC integration do a simple local averaging of the model onto our inversion grid.

Option 1.) would probably be easiest to implement, and the tests could be reproduced by others, but it would be important to make sure that the tomography model is sampled fine enough. Option 2.) is most accurate, but we would need to package the kernel code with external code. Option 3.) is also easy to implement, but local averaging is perhaps not as good as a mc integration.

I would tend towards option 1 since it would allow that others could do similar and more detailed tests most easily. If one samples the tomography models fine enough there shouldn't be dramatic interpolation artifacts.

What do you think?

[sstaehler]

I tend to prefer a completely external solution. The code is already quite huge (34k lines) and if the external test is in a Python module, other users can run it as well. @martinvandriel, your opinion?

[auerl]

Ok! With an external solution it would however probably be more difficult (or at least more work) to perform the test with different meshes. I have routines to project models onto a voxel mesh, but doing the same for a tetrahedral mesh and your tomography basis would involve translating some parts of random_tetrahedron.f90 to python (or do you have them). Doing the test inside the kerner would allow reusing all these things. Or enough to do the tests for voxel meshes?

[martinvandriel]

I think it would probably be good to have a proper projection of models onto out grids to minimize errors especially for large elements. About the implementation, I am unsure, I can see both the argument of code complexity and reusability...

[auerl]

What about adding it in a new branch, and deciding later if it makes sense to adding it to the the main branch?

[sstaehler]

Sounds good. I have been convinced by your arguments pro adding it to the kerner code!

[auerl]

I have added mc integration of 3D heterogeneities and their multiplication with the kernel. @sstaehler: could you review how i implemented it? See https://github.com/auerl/kerner/commits/master

Here an example of how savani looks like when projected onto an (admittedly rather coarse) 5x5deg kernel grid:

[sstaehler]

@auerl: The implementation looks very nice! What are the results when you run it? Are the predicted dT approximately correct?

[auerl]

I am just comparing the predicted dt's for a coarse and a fine mesh and the model SMEAN. Interestingly, I get smaller values for the finer mesh... Maybe due to an error in computing the cell volumes (will check that). Generally, the values are not completely off. They are all negative (maybe due to the heterogeneity structure at this particular location), and differ for the different frequencies.

1.25 deg mesh: 90deg_P_40s : -0.23435E-01 s 90deg_P_30s : -0.13048E-01 s 90deg_P_20s : -0.73470E-02 s 90deg_P_15s : -0.17405E-02 s

5.0 deg mesh: 90deg_P_40s : -0.12113E+00 s 90deg_P_30s : -0.12387E+00 s 90deg_P_20s : -0.10615E+00 s 90deg_P_15s : -0.93786E-01 s

The phi-jumps in the kernel are due to the irregular parameterization in our voxel mesh towards the poles.

[auerl]

And here for a tretraheder mesh. The values are again different.

Tetrahedral (~300km) mesh: 90deg_P_40s : -0.33740E-01 s 90deg_P_30s : -0.30487E-01 s 90deg_P_20s : -0.31019E-01 s 90deg_P_15s : -0.27710E-01 s

[sstaehler]

Hmm... Interesting. Are the time windows selected well? By looking at the seismogram_* and seismogramcut* files?

[auerl]

Hm, good point... the default time windows are for a source at the surface, right? I am just running a test where I multiplied the anomalies with a factor of three...

[sstaehler]

I think so.

[auerl]

Even more important: For one and the same mesh, the values for different runs differ quite significantly.

This is for two runs and the 5deg grid:

30deg_P_40s : -0.76981E-01 s 30deg_P_30s : -0.60391E-01 s 30deg_P_20s : 0.25605E-01 s 30deg_P_15s : 0.41367E-01 s

30deg_P_40s : -0.21920E-01 s 30deg_P_30s : -0.25863E-01 s 30deg_P_20s : -0.54780E-01 s 30deg_P_15s : -0.23115E-01 s

I would guess that this is caused by some not converging grid cells with very high values which dominate the results, but I need to verify this.

[sstaehler]

What did you choose for Maximum number of iterations?

Also, if you run on a completely new setup, it makes sense, to run it once with a low allowed error and maximum number of iterations = 1. Then, look at the resutls and set the allowed error to a sensible value (1e-3*value in interesting regions and relative error to 1e-2) and the number of iterations to a high value (around 1e5). This should ensure convergence in the interesting cells, without wasting time elsewhere.

[auerl]

Here are some comparisons between predicted and measured (specfem) traveltimes.

The tomographic model is s362ani referenced against anisotropic PREM (the original reference model is a model called REF), the earthquake is the default Specfem Bolovia 1994 CMTSOLUTION, the Receiver locations correspond to the default Specfem stations. The comparisons are for the direct P-wave phase. The mesh is a 1x1 degree voxel mesh with 30 layers (so blocks of "roughly" 100x100x100 km, depending on where you are). v_p is simply scaled to v_s by a factor of 0.55

Remaining issues and open questions:

• To get more or less correct predictions (in terms of the right sign of the predicted delay times) and make the kernels look like published kernels, I needed to change the sign in calc_misfit_kernel (kernel.f90; ~line 350). I remember that we discussed about this before. There is an inconsistency between the formulars in Tarjes and other papers. I think in the rho, lambda, mu parameterization we are using, only the rho kernel needs a minus sign.
• The computations are quite costly with the convergence parameters that seem necessary (in terms of predicted dT values not varying to extremely from one run to another). Number of iterations is set to 1e5, absolute error to 1e-2, relative error to 1e-3. It takes ~5 hrs on Mönch (CSCS) for one source receiver combination (kernels for 4 dominant frequencies). Its a bit slower because I also need to project the model onto the inversion mesh (but this is not a significant fraction of the total runtime).
• There is an absolute time offset of around 0.5s between the predictons and the SEM measurements, when the source and the receiver elements are excluded, and maybe a little bit less offset here and there when the source and receiver elements are not masked. Could there be some tiny differences in how the reference model is defined in Specfem and AxiSEM (they look exactly the same, however)? Maybe a slight differences in how attenuation is implemented? Any other ideas where the time shift could be coming from? Some initial processing steps? STF deconvolution?
• There are also some relative differences between the predicted and the measured traveltimes, and I am not sure where they come from. Maybe its just due to the relatively coarse inversion mesh or some elements that don't converge when projecting the background model onto the mesh (at elements extending across background model discontinuities?)
• Neither the predicted nor the measured traveltime delays seem to really depend a lot on frequency, maybe this is due to the tomographic model being smooth?

I am currently rerunning Specfem in the same setup but s362ani multiplied by a factor of 2, and will recompute the kernel predictions subsequently.

Ciao,

Ludwig

Without source & receiver element

With source and receiver element

Example cut through projected tomographic model (this is really how the model looks like, seems pretty weird in the transition zone).

Example cut through a kernel (the kernel looks a bit weird at some locations because its orientation is completely arbitrary and sometimes you cut through the in different "shells" which makes the kernel appear jumpy in paraview).

[martinvandriel]

This looks great, I mean we are really starting to converge and have to work on the details. Some comments and questions:

• the 5 hours, is that node or core hours?
• For PREM in specfem, do you know if they have a proper implementation of the crust? This could be compatible with some of the systematic delay. Did you compare the specfem PREM synthetics with Instaseis PREM?
• Also, it would be interesting to compare this to ray theory / Dahlen Kernels, if we do better, we already have a major point!

Cheers, Martin

[auerl]

Hi Martin,

sorry, i missed to say that its 5 hours on 10 cores, so 50 CPU hours. Really a bit slow. But I should add that the mesh I use here consists of almost 1 Million elements - much more than my typical inversion meshes (I used a fine mesh in order not to introduce additional error due to using a coarse mesh).

In Specfem one can turn off CRUST2.0 (which is the default) by adding "_1Dcrust" to the name of the desired 3d mantel model. So in the s362ani case it would be "s362ani_prem_1Dcrust". Looking at the reference model velocities and discontinuity definitions everything seems consistent with axisem prem_ani, but I didn't compare with the Instaseis PREM synthetics (will do so)

As for the comparison with Dahlen kernels: I think for such a test we might need to add the old kernel code as a module to the new kernel code, so we can switch between computing kernel values at random points from the convolved wavefields or the Dahlen way. Would also allow for a better comparison bw the two types of kernels, since the quadrature would be performed in the exactly same way.

@kasra-hosseini, is there an installation of the old kernel code on our computer in ZH which I could look at?

Ciao,

Ludwig

[kasra-hosseini]

@auerl Thanks Ludwig for the info/figures/comparison with SPECFEM. They look great! As for old kernel code installation, I have these two directories which might be useful:

/home/khosseini/codes/kerner
/home/khosseini/codes/kerner_ludwig

Hope this helps, Kasra

[sstaehler]

@auerl that looks great, thanks a lot!

If the offset is consistently 0.5s, it will not be a problem for tomography, if we allow for event time corrections. Anyway, we should know, where it comes from.

This Bolivia event is very deep, right (I am in a train, can't check)? For the paper, we should have a second test with a normal event to make sure that there are no surprises if the event is in the crust. Can you do that with 2011 Virginia please?

The cost does not sound that outrageous to me since you had an extremely fine mesh. But I don't understand completely: Did you run the code separately for each station or did it take 50h*200 stations = 1e4 CPUh?

Did you already try to put the Dahlen code into the new one? To me, it sounds like a lot of work (e.g. coordinate transform from Cartesian to ray-centred) and I'm not sure about the reward. One advantage might be using Dahlen's ray theoretical solutions for the source element, but maybe that's not necessary for now.

Thanks again!

[auerl]

The 0.5s offset is not really consistent as far as i can tell now. Ok, right, I am going to repeat with a shallow earthquake. The Bolovian one is indeed quite special. Yes, the CPUhrs are not really representative, since I ran separately for each receiver. Still seemed quite extreme to me, but probably its much better than. It was hard to find a node-memory configuration to run the code for all stations in one large job.

I started looking at the Dahlen code and it really seems like a lot of work, mostly translating the old routines. The coordinate transform should already be done somewhere in the Dahlen code. It would be good to have for really fair comparisons (same integration method) in terms of how well the kernels predict measurements, but not sure if I it makes sense working on this....

Another thing I am not sure about is why the Kerner predictions seem to show some dispersion, while the pyffproc measurements from the specfem synthetics do hardly show any frenquency dependence. I'll make a few more plots on this, as soon as I have rerun it for the Virginia earthquake.

[sstaehler]

Hi @auerl, could you upload the results somewhere?

Ad comparisons to the Dahlen code: I do not think that it's really fairer to compare to a MC-Dahlen code. And also I am afraid that the CO-transform from cartesian to ray-centered is probably not in the Dahlen code yet, since it's only needed the other way around.

[auerl]

Hi @sstaehler, I have uploaded the files you would need for the final traveltime comparison to

/scratch2/auerl/TTcomparison_s362ani_bolivia_stf10s_withsrcelem_flippedPolarity_forsimon

I didn't want to upload the kerner rundirectories to the ETH machine, since you are currently running a test. Let me know when you are done.

The pyffproc measurements are stored in the files ffproc.*, they are combined with the kernel based predictions in files result_?? where ?? is the number of the frequency band (i only ran for 4 out of the 8 original frequency bands). Column 9 is the measurement, the rightmost column is the frequency (its the format of karins measurement code, so i guess you know it). Note that all jobs for receivers > Rec Nr. 84 didn't run through bc of exceeding their walltime (large epicentral distances). Also the Predictions for Pdiff don't seem to make sense at the moment, I don't know why. The kernels look normal, I would say. The plots above are only for direct P.

Here another comparison: It seems that the fit is best for the 16.35 second window with 0.133 Hz Center frequency, as said the pyffproc measurements on the specfem synthetics hardly show any dispersion, while the kernel predictions are frequency dependent. Does somebody have an Idea why this could be?

Here the center frequencies for the 4 windows 0.033 Hz 0.066 Hz 0.133 Hz 0.266 Hz

Do we agree about the sign issue in calc_misfit_kernel (kernel.f90; ~line 350). We discussed about this before. There is an inconsistency between the formulars in Tarjes and other papers. I think in the rho, lambda, mu parameterization we are using, only the rho kernel needs a minus sign. I needed to change the code accordingly to get reasonable results.

This example is now with the source and receiver element included, which slows things down a little bit, but also seems to improve the prediction (see above). However, these elements reach the iteration maximum (see picture below). How do we proceed with these elements? Does the two-fold speedup justify just leaving them out? Shouldn't one then mask at least some elements in a zone around the reciever/source, especially when the source is close to a element boundary, avoiding only the actual source/receiver element is maybe not sufficient?

Does the amplitude difference of the kernels make sense?

[auerl]

Here an example for an Pdiff kernel, looks normal, right?

[kasra-hosseini]

@auerl It definitely does! Great! Do you have this VTK file somewhere on ETH machine? Can I have the directory and the file name? What is the frequency again? It is for a deep source with around 130 degrees epicentral distance?

Thanks for that!

[auerl]

Does it make sense that it's amplitude is larger compared to the direct P kernel above? I wonder why the Pdiff kernels don't yet give correct predictions. This is for 0.133 Hz central frequency. It's not on the ETH machine yet, but I can upload it if you want?

[kasra-hosseini]

That would be great if you upload the file somewhere on the ETH machine...OK, I just meant the shape (just comparing to what we have seen so far for Pdiff kernels). I am not sure about the amplitude..look at Liu and Tromp: Fig 2(b) compared to Fig.(3).

[sstaehler]

@auerl Thanks for uploading!

So just to avoid further confusion: In the results_**\ files, the last column is the Kernel prediction?

And do I understand it correctly, that you had to run the code separately for each receiver and each frequency band?

Concerning the minus sign. I am convinced by the results, so let's add it. Did you already change it in sstaehler/kerner?

[sstaehler]

I have plotted the measured vs. predicted dT in a slightly different way that is maybe a bit easier to look at. When separating stations <100° from the others, the difference is quite striking. The measured travel-times seem okay, but the predicted ones go out of the roof. The kernel should not become larger that much for Pdiff. Since the seismogram is in the denominator of the kernel, the decrease in amplitude after 100° increases its amplitude. But this should be canceled by the lower amplitudes in the wavefield. Should...

The differences in band_01 can maybe be explained by the minimum frequency of the wavefield. Which filters did you take? The original Sigloch ones or the output of the kerner?

[sstaehler]

Just for completeness: @auerl's results for the same event and double the perturbations

Looks good!

[auerl]

@sstaehler: we figured that the process of time window selection in (py)ffproc is probably not ideal (as it was tailored for the dahlen kernels) and the windows do often not accomodate the entire waveform wiggle. This is worst for the pdiff arrivals, and might be the reason for the wrong predictions we see. Kasra is working on adapting the window selection to our new case.

[auerl]

Here the predictions with slightly more appropriate time windows, and some other tiny fixes such as that we now make the measurements on displacement seismograms and not velocity seismograms and that the exact same stf was used in both cases.

Unfortunately the pdiff predictions are still off (too high amplitudes), eventhough the timewindows do now more or less capture most of the pdiff energy. See example below. But the direct p predictions look a little bit better than before.

Do you think the bug related to the forward wavefield containing velocity traces (which might lead to kernels having slightly incorrect amplitudes) could have the strongest effects for pdiff and other late phases?

[sstaehler]

I am not sure, whether the velocity wave fields can cause such a big problem, but which version of the code did you use? The new one with displacement wave fields?

[auerl]

No these tests were done before these changes, i am now running another test. Ciao, L.

Von meinem iPhone gesendet

Am 22.08.2015 um 07:02 schrieb Simon Stähler notifications@github.com:

I am not sure, whether the velocity wave fields can cause such a big problem, but which version of the code did you use? The new one with displacement wave fields?

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