Results are influenced by compiler/library versions or machine

[auerl]

Upon Simons request I created a new issue concerning the observation that there are large differences to the kernel values, when the code is compiled with different compiler/netcdf (?) versions.

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:

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

I will definitely do that. I just need to finalize the DMT paper this week. Once I am done, this will be my first task. Until then (hopefully end of this week), I can send you the kernels for any src-rcv pairs. Does that make sense? or is it too late? If it is urgent, I can change my plan...

[tnissen]

@ALL https://github.com/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.

i have computed Dahlen kernels back in the day 10 years ago... we could
try to install that code together on Wednesday/Thursday this week? I
should be at ETH mid afternoon Wed, then again Thursday until lunch...
then next week Mon-Wed. From what I remember it was fairly easy to run
the code, and I certainly haven't changed anything in it!

On the other issues... let's talk Wed/Thur too.

@kasra-hosseini https://github.com/kasra-hosseini: Could you install
such a test-version for us? I think you are the only one of us who
knows the code.

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

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]

I already have a python wrapper for all the codes involved: raydata (for common corrections and ray tracing) + raymatrix that calculates the projected kernels on the inversion grid automatically. There is also a code that can create input files for testing purposes: for instance for one event and different stations (since the input files should be in a specific format). At the end, one can create the vtk file and etc out of the outputs as well. I will try to set it up tomorrow. If it does not work, maybe you can try it while Tarje is in Zurich.

[sstaehler]

@auerl
I am extremely confused by this observation and did not encounter it yet in this magnitude.

From my tests, kernels on a mesh with ~300 km cell size (which I assume was the size of your mesh, just from looking at it) should have a value of ~10^7 in the interesting parts. Which mesh did you use?
Did you multiply the kernels with the volume?

[auerl]

I didn't multiply with the volume here, such to look at the kernel values independently of the mesh. Here another example (this time on a triangular mesh) computed on the Zurich machine, again without multiplying with cell volume... seems like very small kernel values on the Zurich machine are independent of the mesh.

Could you perhaps run my inparam file in /scratch/auerl/axisem0215/kerner_gfortran with the wavefields in /scratch/auerl/axisem0215/kerner_gfortran/wavefield

I noticed that the background model velocities (when printing them in master_queue.f90) are on the order of ~10^11.... something must go wrong somewhere.

[kasra-hosseini]

Hi All,

It was much easier to set-up Dahlen kernels on ETH machine than I expected. Now, we have all the source codes here:

/home/khosseini/dahlen_kernels

If you are interested to know about the source files (like projection on inversion grid and etc), you can find them here:

/home/khosseini/dahlen_kernels/src_raydata_raymatrix

raydata_src: This code does the ray-tracing, calculate the hessian, common corrections and etc ---> at the end, it provides the input files for the next code:
raymatrix_src: in which the calculation of kernels on ray-centric coordinate as well as the projection on the inversion grid happen.

Here is a step-by-step tutorial for generating Dahlen kernels:

cd /home/khosseini/dahlen_kernels/test_maker

Edit test_maker.py as you want. It should be enough to change the parameters only in the INPUT section.

python test_maker.py

This will create 6666.6666.666.a directory in which all the required information for other codes are there.
In the next step, we change the input file of the python wrapper:

vim /home/khosseini/dahlen_kernels/inputs/input_raydata_raymatrix.ini

Based on what we have selected for the bands = ['band03'] in test_maker.py, you can change the input_raydata_raymatrix.ini accordingly (line 22).

And this is all:

./clean_up.sh
python raydata_raymatrix.py inputs/input_raydata_raymatrix.ini

@auerl @martinvandriel Can we install pyvtk on zurich machine?

I have disabled this functionality (creating VTK file) for now. If we have it at some points, the code generates VTK outputs automatically.

All in all the outputs are like this:
RESULTS: contains all the files generated by the codes.
VTK file: in case that we have pyvtk

Example:

[auerl]

Hi Kasra, thanks for setting up the code! I did a sudo pip install pyvtk so it should be installed now. Cheers, L.

[kasra-hosseini]

Great! I just changed the code as well to create the VTK file automatically. I tested it and it works as stated above. Please write me in case that there are some difficulties or etc...

[kasra-hosseini]

OK, as Simon and Ludwig suggested. I set-up a test as follow:

1. inversion grid: 400km everywhere + respecting the CMB (see below)
2. calculating the kernels for 60 degrees and 90 degrees epicentral distance for the following Gabor filters: 30.0sec, 21.2sec, 15.0sec, and 10.6sec.

Results
Please refer to:

/home/khosseini/dahlen_vs_mc

There we have the following dirs/files:
facets.mc_kernel_test and vertices.mc_kernel_test ---> the ones that I used for calculating kernels
60deg and 90deg directories that contain the VTK files for 4 bands. Note that band01 is the lowest frequency band.

@auerl @sstaehler Did I miss something?

Here is the inversion grid:

[kasra-hosseini]

Forgot to say:

event:
latitude: 90
longitude: 0
depth: 0

stations:
latitude: 30.0/0.0
longitude: 0.0/0.0
elevation: 0.0/0.0

[sstaehler]

okay, I'll put that into MC Kernel

[sstaehler]

@kasra-hosseini
What are the exact time windows for the kernels again?

[kasra-hosseini]

For the calculation of kernels with Dahlen method, we do not need to specify any time-window (in the sense that we do in MC). It does it by ray tracing and the filter durations as follow:
band01: 65.1sec
band02: 46.2sec
band03: 32.7sec
band04: 23.1sec

So you can assume the arrival time of P at 60 and 90 which are: 608.28 and 781.33 + the above filter durations.

[sstaehler]

@auerl could you check whether this issue still remains with versions after d9e8ee4
Since you were calculating lambda kernels and the parameter "lam" was mis-interpreted by the code, it might be that the code in earlier revisions just calculated a kernel from any component of the strain tensor, which might have been machine-dependent.

[auerl]

After d9e8ee4 the lambda kernels are the same on Mönch and our local machine. Here is how the kernels look like now. They are computed with the RELATIVE kernel option set to true and volume multiplication turned on. They are on the order of 10^7 and the polarity is flipped, wrt to how they plotted before. I think this issue can be closed.