paudetseis / PyRaysum

Teleseismic body wave modeling through stacks of (dipping/anisotropic) layers
https://paudetseis.github.io/PyRaysum/
MIT License
42 stars 14 forks source link

Computations taking infinite time for some velocity models #20

Open pasansherath opened 1 year ago

pasansherath commented 1 year ago

Hi @wsja

For some velocity models, receiver function computations take an infinite time. When I use verbose=1, I see the following output, and stays the same forever. What could be happening here? I tried tweaking the velocity model in question, but same outcome. I am sorry for opening all these issues, but it seems the only way to resolve them. Thanks in advance!

 phase:          400
 layer: type:
      8     1
      7     1
      6     1
      5     1
      4     1
      3     1
      2     1
      1     1
      1     5
      2     5
      3     5
      4     6
      5     6
      6     6
      7     6
      7     3
      6     3
      5     3
      4     3
      3     2
      2     2
      1     2

 phase:          401
 layer: type:
      8     1
      7     1
      6     1
      5     1
      4     1
      3     1
      2     1
      1     1
      1     5
      2     5
      3     5
      4     6
      5     6
      6     6
      7     6
      7     3
      6     3
      5     3
      4     3
      3     3
      2     3
      1     3
wasjabloch commented 1 year ago

@pasansherath, no reason to apologize. Your comments are very helpful for making PyRaysum better and may be valuable for other users. Thanks for raising them.

I'll come back to you later this week. In the meantime, could you please post the parameters that you were using (model file, and the non-default arguments to Control and Geometry)? Perhaps a minimal working example?

wasjabloch commented 1 year ago

From your output, I understand that you have set mults=2 in combination with a model with 8 layers. I understand and expect that seemingly infinite run times result from this setup, as the code will compute a vast combination of reflected wave types. In the standard configuration, this would yield a segmentation fault. Please see curve Obs-M2-F-RF or Num-M2 of Figure 6 of the PyRaysum paper.

You have several options to proceed from here:

  1. Isolate the specific reflections you are interested in and define them using Control.set_phaselist(descriptors). This is advanced usage and I am happy to get you started with that.
  2. Try to reduce the complexity of your model, using fewer layers
  3. If you do not require dipping layers, you can try Telewavesim instead.
wasjabloch commented 11 months ago

To limit the number of phases, there is now a new method of Model in the dev branch: Model.get_phaselist() gives you fine-grained control over which phases to include in your calculation. This could be useful for your use case. I'd be happy to get some feedback, if you'd like to give it a try :smile:

pasansherath commented 11 months ago

Hi @wsja, Thanks very much for your suggestions, and I am sorry for the delay in responding to you.

I was able to resolve the issue I had by splitting a ~20 km thick layer into two identical ~10 km thick layers.

pasansherath commented 11 months ago

I also need to model all possible phases, therefore setting the phase list is not ideal.