Love wave

Model

We assume the following subsurface velocity and density distribution with 4 layers:

Layer #	thickness (m)	$V_s(km{s}^{-1})$	$\rho (gc{m}^{-3})$
1	8000	4	2
2	4000	3.5	2
3	8000	4.2	2
4	$\mathrm{\infty }$	7.65	2

The LWModel is defined as :

vs = [3.2, 4.39731, 4.40192, 4.40653, 4.41113, 4.41574]
ρ = [2.6, 3.38014, 3.37797, 3.37579, 3.37362, 3.37145]
h = [20.0, 20.0, 20.0, 20.0, 20.0] .* 1e3
m = LWModel(vs, h, ρ)

1D Love Wave Model
┌───────┬────────┬──────────┬───────────┬──────────┐
│ Layer │     vₛ │        h │         z │        ρ │
│       │ [km/s] │      [m] │       [m] │ [kg⋅m⁻³] │
├───────┼────────┼──────────┼───────────┼──────────┤
│     1 │   3.20 │ 20000.00 │  20000.00 │      2.6 │
│     2 │   4.40 │ 20000.00 │  40000.00 │     3.38 │
│     3 │   4.40 │ 20000.00 │  60000.00 │    3.378 │
│     4 │   4.41 │ 20000.00 │  80000.00 │    3.376 │
│     5 │   4.41 │ 20000.00 │ 100000.00 │    3.374 │
│     6 │   4.42 │        ∞ │         ∞ │    3.371 │
└───────┴────────┴──────────┴───────────┴──────────┘

Code for this figure

fig = Figure()
ax = Axis(fig[1, 1]; xlabel="Vs (km/s)", ylabel="depth (m)", backgroundcolor=(
    :magenta, 0.05))

plot_model!(ax, m)

Warn

Always use Float64 or Float32 types while defining the vectors for velocities, densities and thickness. Using Int will throw an InexactError, e.g. : InexactError: Int64(4193.453970907305)

Response

To obtain the response, that is the dispersion curve c. We first need to define the frequencies :

freq = exp10.(-2:0.1:1)
T = inv.(freq)

31-element Vector{Float64}:
 100.0
  79.43282347242813
  63.09573444801933
  50.11872336272722
  39.810717055349734
  31.622776601683796
  25.118864315095795
  19.9526231496888
  15.848931924611133
  12.589254117941675
   ⋮
   0.6309573444801932
   0.5011872336272724
   0.39810717055349726
   0.31622776601683794
   0.251188643150958
   0.199526231496888
   0.15848931924611132
   0.1258925411794167
   0.1

and then call the forward dispatch and get SurfaceWaveResponse:

resp = forward(m, T)

Code for this figure

fig = Figure() #size = (800, 600))
ax = Axis(fig[1, 1]; backgroundcolor=(:magenta, 0.05), aspect=AxisAspect(1))

plot_response!([ax], T, resp; plt_type=:scatter)

Note

By default, the forward response calculates the phase velocity at fundamental mode. This information is passed through params.

default_params(LWModel)

(mode = 0, dc = 0.01, dt = 0.001, type = Val{:phase}(), c_low = 3.0, c_high = 6.0)

Various symbols are defined as :

• mode : fundamental mode corresponds to 0, higher modes correspond to 1,2,...

• dc : step size used to obtain solution of the propagater matrix

• dt : step size used to obtain the derivative for group velocity, unused for phase velocity

• Val(:phase) : to calculate phase velocity, pass Val(:group) to calculate group velocity

In-place operations

Mutating forward calls are also supported. This leads to no extra allocations while calculating the forward responses. This also speeds up the performance, though the calculations in the well-optimized forward calls are way more expensive than allocations to get a significant boost here. We now call the in-place variant forward! and provide a SurfaceWaveResponse variable to be overwritten :

forward!(resp, m, T)

We can compare the allocations made in the two calls:

@allocated forward!(resp, m, T)

0

@allocated forward(m, T)

320

Benchmarks

While the exact runtimes will vary across processors, the runtimes increase linearly with number of layers.

Benchmark code

n_layers = Int.(exp2.(2:1:6))
freq = exp10.(-1:0.2:2)
T = 2π .* freq

btime_iip = Float64.(zero(n_layers))
btime_oop = Float64.(zero(n_layers))
for i in eachindex(n_layers)
    vs = 3 .+ 0.005 .* randn(n_layers[i])
    ρ = 2 .* ones(n_layers[i])
    h = 100 .* rand(n_layers[i]-1)
    m = LWModel(vs, h, ρ)
    resp = forward(m, T)
    time_oop = @timed begin
        for i in 1:1000
            forward(m, T)
        end
    end
    btime_oop[i] = time_oop.time/1e3
    forward!(resp, m, T)
    time_iip = @timed begin
        for i in 1:1000
            forward!(resp, m, T)
        end
    end
    btime_iip[i] = time_iip.time/1e3
end

fig = Figure()
ax = Axis(fig[1, 1]; xscale=log2, backgroundcolor=(:magenta, 0.05),
    xlabel="no. of layers", ylabel="time (ms)")
lines!(n_layers, btime_oop .* 1e3; label="out-of place forward calls", color=:tomato)
lines!(n_layers, btime_iip .* 1e3; label="in place forward calls",
    linestyle=:dash, color=:steelblue3)
Legend(fig[1, 2], ax)

┌ Warning: Assignment to `vs` in soft scope is ambiguous because a global variable by the same name exists: `vs` will be treated as a new local. Disambiguate by using `local vs` to suppress this warning or `global vs` to assign to the existing global variable.
└ @ love.md:138
┌ Warning: Assignment to `ρ` in soft scope is ambiguous because a global variable by the same name exists: `ρ` will be treated as a new local. Disambiguate by using `local ρ` to suppress this warning or `global ρ` to assign to the existing global variable.
└ @ love.md:139
┌ Warning: Assignment to `h` in soft scope is ambiguous because a global variable by the same name exists: `h` will be treated as a new local. Disambiguate by using `local h` to suppress this warning or `global h` to assign to the existing global variable.
└ @ love.md:140
┌ Warning: Assignment to `m` in soft scope is ambiguous because a global variable by the same name exists: `m` will be treated as a new local. Disambiguate by using `local m` to suppress this warning or `global m` to assign to the existing global variable.
└ @ love.md:141
┌ Warning: Assignment to `resp` in soft scope is ambiguous because a global variable by the same name exists: `resp` will be treated as a new local. Disambiguate by using `local resp` to suppress this warning or `global resp` to assign to the existing global variable.
└ @ love.md:142

Reproducibility

The above benchmarks were obtained on AMD EPYC 9V74 80-Core Processor