Spherical point source solver sometimes fails to converge

[jbmuir]

I'm trying to use the point source solver in spherical coordinates, and have found that for my problem and for some source locations, the algorithm fails to converge properly, instead giving an array of infinities. Its not clear to me what the bug is that is causing this. Changing the value of the phi coordinate slightly seems to stop this problem from happening for this particular example.

import pykonal
import numpy as np

lamax, lomin, depmin = 46.750, -113.0, -3.0
lanode, lonode, depnode = 32, 41, 164
ladx, lodx, depdx = 0.125, 0.125, 1.0

rmin = 6371 - depmin + 1 - depnode*depdx
thmin = np.deg2rad(90-lamax)
phmin = np.deg2rad(lomin)
rnode, thnode, phnode = depnode, lanode, lonode
rdx = depdx
phdx = np.deg2rad(lodx)
thdx = np.deg2rad(ladx)

model = 5*np.ones((rnode, thnode, phnode))

solver = pykonal.solver.PointSourceSolver(coord_sys="spherical")
solver.velocity.min_coords = rmin, thmin, phmin
solver.velocity.node_intervals = rdx, thdx, phdx
solver.velocity.npts = rnode, thnode, phnode
solver.velocity.values = model
solver.src_loc = (6373.149, 0.7995353303386024, -1.9405094222448556)
solver.solve()
solver.traveltime.values

[malcolmw]

Hi, @jbmuir, Thanks for opening this issue. This is likely due to an internal inconsistency in how I handle the range of phi (i.e. phi in [0, 2 pi) versus phi in [-pi, pi); see Issue #8). The most recent patch addressed this bug but seems to fail in your case. I will take a look and get back to you.

[malcolmw]

@jbmuir, it looks like this is actually a different issue than what I originally thought.

Just add solver.nrho *= 2 right above solver.solve(), and your code above should work.

The PointSourceSolver works by combining two computational grids: (a) a refined near-field grid and (b) a coarse far-field grid. The user defines the far-field grid as he would with the ordinary EikonalSolver, and the code automagically defines the near-field grid "appropriately" (I use quotations because this was obviously failing). The near-field grid was being defined such that no nodes from the far-field grid fell within the near-field grid, so the transition between grids failed. Increasing the size of the near-field grid (via the solver.nrho attribute) resolves the issue. This solution is a hack, but will keep you going for now. I will work on a proper patch when I have time.

Let me know if that fix works for you, and I will close this issue if it does.

And Happy New Year!

[jbmuir]

Hi Malcolm, looks like it works! I'm excited to try this code out (writing a quick julia wrapper for it with PyCall.jl and it seems to be a great solution for traveltime calculations)

[malcolmw]

Great! Let me know if you make the Julia wrapper public, and I will put a pointer to it in the documentation.

Closing this Issue as solved.

[jbmuir]

I have something working now; I don't think I will wrap the entire API just yet as I don't have a need for it and my PhD has to finish some time, but I'll add the necessary components below:

[jbmuir]

to go in the deps/ folder as build.jl

using Conda

Conda.add("cython>=0.29.14")
Conda.add("h5py")
Conda.add("numpy")
Conda.add("scipy")
pykonal_command_initial = `git clone https://github.com/malcolmw/pykonal`

if !ispath("pykonal")
    run(pykonal_command_initial)
else
    cd("pykonal")
    run(`git pull`)
    cd("..")
end

Conda.pip_interop(true)
cd("pykonal")
Conda.pip("install", ".")
cd("..")

[jbmuir]

and a function for running the point source solver (currently specialized for a cube surrounding Yellowstone, but easy enough to change)

struct Station{T<:Real}
    name::String
    latitude::T
    longitude::T
    elevation::T
    theta::T
    phi::T
    radius::T
end

Station(name, latitude::T, longitude::T, elevation::T) where T = Station(name, latitude, longitude, elevation, deg2rad(90-latitude), deg2rad(longitude), convert(T, 6371.0+elevation/1000.0))

using PyCall

const pykonal = pyimport("pykonal")

# constants describe specific study area, and conversion into spherical from LLA coordinates
const lamax, lomin, depmin = 46.750, -113.0, -3.0
const lanode, lonode, depnode = 32, 41, 164
const ladx, lodx, depdx = 0.125, 0.125, 1.0

const rmin = 6371 - depmin + 1 - depnode*depdx
const thmin = deg2rad(90-lamax)
const phmin = deg2rad(lomin)
const rnode, thnode, phnode = depnode, lanode, lonode
const rdx = depdx
const phdx = deg2rad(lodx)
const thdx = deg2rad(ladx)

function station_travel_times(station::Station{T}, model::Array{T,3}) where T
    solver = pykonal.solver.PointSourceSolver(coord_sys="spherical")
    solver.velocity.min_coords = rmin, thmin, phmin
    solver.velocity.node_intervals = rdx, thdx, phdx
    solver.velocity.npts = rnode, thnode, phnode
    solver.velocity.values = model
    solver.src_loc = station.radius, station.theta, station.phi
    solver.nrho = 2*solver.nrho
    solver.solve()
    return solver.traveltime.values
end

[malcolmw]

Awesome, thanks for sharing @jbmuir. I will likely include this as short tutorial or example.