From 9683416ac6dab7b6d80c8e16894cddfeafaaa8b9 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Fri, 20 Oct 2023 10:55:38 +1100 Subject: [PATCH 01/10] VTEM pass through tx_rx_dz --- .../VTEM/DarlingParoo/McMC/01_read_data.jl | 1 + .../gradientbased/01_read_data.jl | 1 + .../synth/gradientbased/02_set_options.jl | 2 +- src/VTEM1DInversion.jl | 33 ++++++++++++------- 4 files changed, 24 insertions(+), 13 deletions(-) diff --git a/examples/VTEM/DarlingParoo/McMC/01_read_data.jl b/examples/VTEM/DarlingParoo/McMC/01_read_data.jl index d784a1fa..d57f8edd 100644 --- a/examples/VTEM/DarlingParoo/McMC/01_read_data.jl +++ b/examples/VTEM/DarlingParoo/McMC/01_read_data.jl @@ -15,6 +15,7 @@ d = [177,221] soundings = transD_GP.VTEM1DInversion.read_survey_files(; X, Y, Z, fid, linenum, frame_height, d, fname_dat, fname_specs_halt, + tx_rx_dz_pass_through = 0.01, # GA-AEM Z up, +ve is rx over tx startfrom = 1, skipevery = 1, dotillsounding = 2, diff --git a/examples/VTEM/DarlingParoo/gradientbased/01_read_data.jl b/examples/VTEM/DarlingParoo/gradientbased/01_read_data.jl index f6711da6..d0377e5f 100644 --- a/examples/VTEM/DarlingParoo/gradientbased/01_read_data.jl +++ b/examples/VTEM/DarlingParoo/gradientbased/01_read_data.jl @@ -15,6 +15,7 @@ d = [177,221] soundings = transD_GP.VTEM1DInversion.read_survey_files(; X, Y, Z, fid, linenum, frame_height, d, fname_dat, fname_specs_halt, + tx_rx_dz_pass_through = 0.01, # GA-AEM Z up, +ve is rx over tx startfrom = 1, skipevery = 50, dotillsounding = nothing, diff --git a/examples/VTEM/synth/gradientbased/02_set_options.jl b/examples/VTEM/synth/gradientbased/02_set_options.jl index 1c7aadaa..c34b649b 100644 --- a/examples/VTEM/synth/gradientbased/02_set_options.jl +++ b/examples/VTEM/synth/gradientbased/02_set_options.jl @@ -7,4 +7,4 @@ regtype = :R2 # :R1 is first difference regularization, :R2 is two succesive fir lo, hi = -3, 1 # min, max log10 conductivity λ²min, λ²max=-0.5, 7 # min, max of regularization parameter (Tikhonov parameter) ntries = 8 # how many max tries between λ²min, λ²max -β² = 0. # what fraction of roughness penalty to use to enforce a penalty in deviations from background \ No newline at end of file +β² = 0.001 # what fraction of roughness penalty to use to enforce a penalty in deviations from background \ No newline at end of file diff --git a/src/VTEM1DInversion.jl b/src/VTEM1DInversion.jl index ec73ae75..126949a3 100644 --- a/src/VTEM1DInversion.jl +++ b/src/VTEM1DInversion.jl @@ -34,6 +34,7 @@ function dBzdt(;times = [1.], ramp = [0. 1], rTx = 10., zTx = -30., + zRx = -30.01, # now Z down for use useML = false, z = [-1.], ρ = [-1.], @@ -53,7 +54,7 @@ function dBzdt(;times = [1.], F = AEM_VMD_HMD.HFieldDHT(; times, ramp, nmax, zTx, rTx, - rRx = 0., zRx = zTx-0.01, + rRx = 0., zRx, calcjacobian, nfreqsperdecade, ntimesperdecade, modelprimary ) @@ -100,6 +101,7 @@ mutable struct VTEMsoundingData <: Sounding fid :: Float64 linenum :: Int zTx :: Float64 + zRx :: Float64 rTx :: Float64 lowpassfcs :: Array{Float64, 1} times :: Array{Float64, 1} @@ -109,14 +111,15 @@ mutable struct VTEMsoundingData <: Sounding end returnforwrite(s::VTEMsoundingData) = [s.X, s.Y, s.Z, s.fid, - s.linenum, s.zTx, s.rTx] + s.linenum, s.zTx, s.zRx, s.rTx] function getndata(S::VTEMsoundingData) getndata(S.data)[1] end function VTEMsoundingData(;rRx=nothing, zRx=nothing, zTx=12., - rTx=-12.,lowpassfcs=[], + rTx=-12.,lowpassfcs=[], + tx_rx_dz = -0.01, # not reading but using format Z down -ve is rx above tx times=[1., 2.], ramp=[1 2; 3 4], noise=[1.], data=[1.], sounding_string="sounding", X=nothing, Y=nothing, Z=nothing, @@ -124,13 +127,13 @@ function VTEMsoundingData(;rRx=nothing, zRx=nothing, zTx=12., @assert rTx > 0 @assert zTx < 0 # my coordinate system z down isnothing(rRx) && (rRx = 0.) - isnothing(zRx) && (zRx = zTx-0.1) # place receiver just above tx centre + isnothing(zRx) && (zRx = zTx + tx_rx_dz) !isempty(lowpassfcs) && @assert all(lowpassfcs .> 0) @assert all(diff(times) .>0 ) @assert all(diff(ramp[:,1]) .>0 ) @assert all((noise .>0) .| isnan.(noise)) @assert length(data) == length(noise) - VTEMsoundingData(sounding_string, X, Y, Z, fid, linenum, zTx, rTx, + VTEMsoundingData(sounding_string, X, Y, Z, fid, linenum, zTx, zRx, rTx, lowpassfcs, times, ramp, noise, data) end @@ -152,6 +155,7 @@ function read_survey_files(; Z = -1, fid = -1, linenum = -1, + tx_rx_dz_pass_through = 0.01, # Z up GA-AEM reading convention +ve is rx above tx nanchar = "*") @assert frame_height > 0 @@ -171,8 +175,10 @@ function read_survey_files(; fiducial = soundings[:,fid] whichline = soundings[:,linenum] d = soundings[:,d[1]:d[2]] - zTx = -soundings[:,frame_height] # my coordinate system - + zTx = soundings[:,frame_height] # read in Z up + zRx = -(zTx .+ tx_rx_dz_pass_through) # my coordinate system Z down + zTx = -zTx # my coordinate system Z down + @info "reading $fname_specs_halt" include(fname_specs_halt) @assert size(d, 2) == length(times) @@ -181,13 +187,13 @@ function read_survey_files(; d[:] .*= units σ = sqrt.((multnoise*d).^2 .+ (σ_halt').^2) - makeqcplots && plotsoundingdata(d, σ, times, zTx; figsize, fontsize) + makeqcplots && plotsoundingdata(d, σ, times, zTx, zRx; figsize, fontsize) nsoundings = size(d, 1) s_array = Array{VTEMsoundingData, 1}(undef, nsoundings) fracdone = 0 for is in 1:nsoundings l, f = Int(whichline[is]), fiducial[is] - s_array[is] = VTEMsoundingData(;zTx=zTx[is], rTx, + s_array[is] = VTEMsoundingData(;zTx=zTx[is], zRx=zRx[is], rTx, times, ramp, noise=σ[is,:], data=d[is,:], sounding_string="sounding_$(l)_$f", X=easting[is], Y=northing[is], Z=topo[is], fid=f, @@ -201,7 +207,7 @@ function read_survey_files(; return s_array end -function plotsoundingdata(d, σ, times, zTx; figsize=(8,4), fontsize=1) +function plotsoundingdata(d, σ, times, zTx, zRx; figsize=(8,4), fontsize=1) f, ax = plt.subplots(2, 2, figsize=figsize, gridspec_kw=Dict("width_ratios" => [1,0.01])) nsoundings = size(d, 1) plot_d = permutedims(d) @@ -216,9 +222,11 @@ function plotsoundingdata(d, σ, times, zTx; figsize=(8,4), fontsize=1) axx.semilogy(mean(σ./abs.(d), dims=1)[:], times, "r") axx.semilogy(mean(σ./abs.(d), dims=1)[:], times, "--w") axx.set_xlabel("avg noise fraction") + ax[2].plot(1:nsoundings, zRx, label="zRx") ax[2].plot(1:nsoundings, zTx, label="zTx") ax[2].set_xlabel("sounding #") ax[2].set_ylabel("height m") + ax[2].invert_yaxis() ax[2].sharex(ax[1]) ax[2].legend() ax[2,2].axis("off") @@ -379,7 +387,7 @@ function makeoperator(sounding::VTEMsoundingData; ρ[z.>=zstart] .= ρbg aem = dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary, times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade, - rTx=sounding.rTx, zTx=sounding.zTx, z, ρ, calcjacobian, useML, showgates=plotfield) + rTx=sounding.rTx, zTx=sounding.zTx, zRx=sounding.zRx, z, ρ, calcjacobian, useML, showgates=plotfield) plotfield && plotmodelfield!(aem, log10.(ρ[2:end])) aem, zall, znall, zboundaries end @@ -390,7 +398,8 @@ function makeoperator(aem::dBzdt, sounding::VTEMsoundingData) modelprimary = aem.F.useprimary === 1. ? true : false dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary, times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade, - rTx=sounding.rTx, zTx=sounding.zTx, z=copy(aem.z), ρ=copy(aem.ρ), + rTx=sounding.rTx, zTx=sounding.zTx, zRx=sounding.zRx, + z=copy(aem.z), ρ=copy(aem.ρ), aem.F.calcjacobian, aem.useML, showgates=false) end From 91d02e04c35c2aa9e894fbec725dddf1f866dc73 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Sun, 22 Oct 2023 14:38:11 +1100 Subject: [PATCH 02/10] on ramp conv works with low pass --- .../gradientbased/electronics_halt.jl | 2 +- .../synth/forwardtests/02_compare_ross.jl | 4 ++-- .../synth/waveletapprox/electronics_halt.jl | 2 +- src/AEM_VMD_HMD.jl | 22 ++++++++++++++----- src/GradientInversion.jl | 14 ++++++++++-- src/VTEM1DInversion.jl | 13 ++++++----- 6 files changed, 40 insertions(+), 17 deletions(-) diff --git a/examples/VTEM/DarlingParoo/gradientbased/electronics_halt.jl b/examples/VTEM/DarlingParoo/gradientbased/electronics_halt.jl index 7b9a00c2..80accdd8 100644 --- a/examples/VTEM/DarlingParoo/gradientbased/electronics_halt.jl +++ b/examples/VTEM/DarlingParoo/gradientbased/electronics_halt.jl @@ -50,5 +50,5 @@ times = vec(10 .^mean(log10.([ 0.0100851 0.0113498 ]), dims=2)) include("hedgehog.jl") - +lowpassfcs = [] σ_halt = vec([0.238388 0.175366 0.126135 0.091645 0.071083 0.043692 0.030069 0.019933 0.015175 0.010694 0.008489 0.006127 0.005354 0.005058 0.003871 0.003524 0.002759 0.002185 0.001792 0.001505 0.001298 0.001112 0.000946 0.000805 0.000689 0.000620 0.000546 0.000465 0.000428 0.000399 0.000352 0.000333 0.000317 0.000293 0.000302 0.000334 0.000406 0.000465 0.000458 0.000356 0.000292 0.000193 0.000285 0.000195 0.000177]) \ No newline at end of file diff --git a/examples/VTEM/synth/forwardtests/02_compare_ross.jl b/examples/VTEM/synth/forwardtests/02_compare_ross.jl index bac90612..0f41c0a0 100644 --- a/examples/VTEM/synth/forwardtests/02_compare_ross.jl +++ b/examples/VTEM/synth/forwardtests/02_compare_ross.jl @@ -14,10 +14,10 @@ rRx = 0.0 zRx = -30.0 zTx = -30.0 freqlow = 1e-3 -freqhigh = 1e6 +freqhigh = 5e5 nkᵣeval = 50 ntimesperdecade = 10 -nfreqsperdecade = 6 +nfreqsperdecade = 5 calcjacobian = true include("../waveletapprox/electronics_halt.jl") ## LM operator diff --git a/examples/VTEM/synth/waveletapprox/electronics_halt.jl b/examples/VTEM/synth/waveletapprox/electronics_halt.jl index 7b9a00c2..80accdd8 100644 --- a/examples/VTEM/synth/waveletapprox/electronics_halt.jl +++ b/examples/VTEM/synth/waveletapprox/electronics_halt.jl @@ -50,5 +50,5 @@ times = vec(10 .^mean(log10.([ 0.0100851 0.0113498 ]), dims=2)) include("hedgehog.jl") - +lowpassfcs = [] σ_halt = vec([0.238388 0.175366 0.126135 0.091645 0.071083 0.043692 0.030069 0.019933 0.015175 0.010694 0.008489 0.006127 0.005354 0.005058 0.003871 0.003524 0.002759 0.002185 0.001792 0.001505 0.001298 0.001112 0.000946 0.000805 0.000689 0.000620 0.000546 0.000465 0.000428 0.000399 0.000352 0.000333 0.000317 0.000293 0.000302 0.000334 0.000406 0.000465 0.000458 0.000356 0.000292 0.000193 0.000285 0.000195 0.000177]) \ No newline at end of file diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl index d72f088f..31ce383b 100644 --- a/src/AEM_VMD_HMD.jl +++ b/src/AEM_VMD_HMD.jl @@ -103,8 +103,18 @@ function HFieldDHT(; ϵᵢ = similar(pz) rTE = zeros(length(pz)-1) rTM = similar(rTE) - if freqhigh < 3/minimum(times) - freqhigh = 3/minimum(times) + + mintime = minimum(times) + ## all this below for pesky zero time at start instead of shutoff ... + @assert mintime > 0 # else interpolation in log10 trouble + if doconvramp + if isapprox(ramp[1,1], 0, atol=1e-8) + mintime = ramp[2,1] # should be around 1.28e-5 to 2.5e-5 for HeliTEM + # @warn "automatically setting impulse response mintime to $mintime" + end + end + if freqhigh < 3/mintime + freqhigh = 3/mintime end if freqlow > 0.33/maximum(times) freqlow = 0.33/maximum(times) @@ -115,8 +125,8 @@ function HFieldDHT(; interpkᵣ = 10 .^range(minimum(log10.(Filter_base))-0.5, maximum(log10.(Filter_base))+0.5, length = nkᵣeval) J0_kernel_h, J1_kernel_h, J0_kernel_v, J1_kernel_v = map(x->zeros(ComplexF64, length(interpkᵣ), length(freqs)), 1:4) J01kernelhold = zeros(ComplexF64, length(Filter_base)) - log10ω = log10.(2*pi*freqs) - interptimes = 10 .^(minimum(log10.(times))-1:1/ntimesperdecade:maximum(log10.(times))+1) + log10ω = log10.(2*pi*freqs) + interptimes = 10 .^(log10(mintime)-1:1/ntimesperdecade:maximum(log10.(times))+1) HFD_z = zeros(ComplexF64, length(freqs)) # space domain fields in freq HFD_r = zeros(ComplexF64, length(freqs)) # space domain fields in freq HFD_az = zeros(ComplexF64, length(freqs)) # space domain fields in freq @@ -623,7 +633,7 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu dI = F.ramp[iramp+1,2] - F.ramp[iramp,2] dIdt = dI/dt - if rta > F.times[itime] + if rta >= F.times[itime] # geq instead of eq as we could have an unlcky time break end if rtb > F.times[itime] # end in this interval @@ -631,7 +641,7 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu end ta = F.times[itime]-rta - tb = max(F.times[itime]-rtb, 1e-8) # rtb > rta, so make sure this is not zero... + tb = max(F.times[itime]-rtb, 1e-8) # rtb > rta, so make sure this is not zero because integ is in log10... a, b = log10(ta), log10(tb) x, w = F.quadnodes, F.quadweights F.dBzdt[itime] += (b-a)/2*dot(getrampresponse((b-a)/2*x .+ (a+b)/2, splz), w)*dIdt diff --git a/src/GradientInversion.jl b/src/GradientInversion.jl index a7edb481..087e10d8 100644 --- a/src/GradientInversion.jl +++ b/src/GradientInversion.jl @@ -261,7 +261,12 @@ function gradientinv( m::AbstractVector, istep += 1 istep > nstepsmax && break end - isa(io, Nothing) || begin @info "Finished "*fname; close(io) end + if !isnothing(io) + @info "Finished "*fname + close(io) + nothing + return + end return map(x->x[1:istep], (mnew, χ², λ², oidx)) end @@ -334,7 +339,12 @@ function gradientinv( m::AbstractVector, istep += 1 istep > nstepsmax && break end - isa(io, Nothing) || begin @info "Finished "*fname; close(io) end + if !isnothing(io) + @info "Finished "*fname + close(io) + nothing + return + end # nuisances mess up occam so this is a kludge as sometimes a smoothest root is not found (istep > nstepsmax) && (istep -= 1) return map(x->x[1:istep], (mnew, nunew, χ², χ²nu, λ², oidx)) diff --git a/src/VTEM1DInversion.jl b/src/VTEM1DInversion.jl index 126949a3..0d3c09e8 100644 --- a/src/VTEM1DInversion.jl +++ b/src/VTEM1DInversion.jl @@ -47,6 +47,9 @@ function dBzdt(;times = [1.], σ = zeros(0), showgates = false, modelprimary = false, + lowpassfcs = [], + freqlow = 1e-3, + freqhigh = 5e5, ) @assert size(σ) == size(d) @@ -54,7 +57,8 @@ function dBzdt(;times = [1.], F = AEM_VMD_HMD.HFieldDHT(; times, ramp, nmax, zTx, rTx, - rRx = 0., zRx, + rRx = 0., zRx, + lowpassfcs, freqlow, freqhigh, calcjacobian, nfreqsperdecade, ntimesperdecade, modelprimary ) @@ -186,7 +190,6 @@ function read_survey_files(; σ_halt[:] .*= units d[:] .*= units σ = sqrt.((multnoise*d).^2 .+ (σ_halt').^2) - makeqcplots && plotsoundingdata(d, σ, times, zTx, zRx; figsize, fontsize) nsoundings = size(d, 1) s_array = Array{VTEMsoundingData, 1}(undef, nsoundings) @@ -194,7 +197,7 @@ function read_survey_files(; for is in 1:nsoundings l, f = Int(whichline[is]), fiducial[is] s_array[is] = VTEMsoundingData(;zTx=zTx[is], zRx=zRx[is], rTx, - times, ramp, noise=σ[is,:], data=d[is,:], + times, ramp, noise=σ[is,:], data=d[is,:], lowpassfcs, sounding_string="sounding_$(l)_$f", X=easting[is], Y=northing[is], Z=topo[is], fid=f, linenum=l) @@ -386,7 +389,7 @@ function makeoperator(sounding::VTEMsoundingData; z, ρ, = makezρ(zboundaries; zfixed, ρfixed) ρ[z.>=zstart] .= ρbg aem = dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary, - times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade, + times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade, lowpassfcs=sounding.lowpassfcs, rTx=sounding.rTx, zTx=sounding.zTx, zRx=sounding.zRx, z, ρ, calcjacobian, useML, showgates=plotfield) plotfield && plotmodelfield!(aem, log10.(ρ[2:end])) aem, zall, znall, zboundaries @@ -396,7 +399,7 @@ function makeoperator(aem::dBzdt, sounding::VTEMsoundingData) ntimesperdecade = gettimesperdec(aem.F.interptimes) nfreqsperdecade = gettimesperdec(aem.F.freqs) modelprimary = aem.F.useprimary === 1. ? true : false - dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary, + dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary, lowpassfcs=sounding.lowpassfcs, times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade, rTx=sounding.rTx, zTx=sounding.zTx, zRx=sounding.zRx, z=copy(aem.z), ρ=copy(aem.ρ), From 3397e45d8791a0a35b7d33ff76656db3eb8b9394 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Fri, 27 Oct 2023 15:43:43 +1100 Subject: [PATCH 03/10] Xcite works but check: -what is min time to go back to, ie., ta, tb can't be less than mintime as otherwise interpolation at those times is nonsensical --- .../synth/forwardtests/02_compare_ross.jl | 4 +- src/AEM_VMD_HMD.jl | 43 +++++++++++++++---- src/TEMPEST1DInversion.jl | 2 +- src/VTEM1DInversion.jl | 4 +- test/verify_skytem_td.jl | 2 +- 5 files changed, 40 insertions(+), 15 deletions(-) diff --git a/examples/VTEM/synth/forwardtests/02_compare_ross.jl b/examples/VTEM/synth/forwardtests/02_compare_ross.jl index 0f41c0a0..079efc49 100644 --- a/examples/VTEM/synth/forwardtests/02_compare_ross.jl +++ b/examples/VTEM/synth/forwardtests/02_compare_ross.jl @@ -13,8 +13,8 @@ z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed) rRx = 0.0 zRx = -30.0 zTx = -30.0 -freqlow = 1e-3 -freqhigh = 5e5 +freqlow = 1e-4 +freqhigh = 1e6 nkᵣeval = 50 ntimesperdecade = 10 nfreqsperdecade = 5 diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl index 31ce383b..cea6f9e9 100644 --- a/src/AEM_VMD_HMD.jl +++ b/src/AEM_VMD_HMD.jl @@ -21,6 +21,7 @@ mutable struct HFieldDHT <: HField ramp :: Array{Float64, 2} log10ω :: Array{Float64, 1} interptimes :: Array{Float64, 1} + minresptime HFD_z :: Array{ComplexF64, 1} HFD_r :: Array{ComplexF64, 1} HFD_az :: Array{ComplexF64, 1} @@ -92,6 +93,7 @@ function HFieldDHT(; getazimH = false, freqlow = 1e-4, freqhigh = 1e6, + minresptime = 1.e-6, # I think responses earlier than this are unstable calcjacobian = false ) @assert all(freqs .> 0.) @@ -104,29 +106,32 @@ function HFieldDHT(; rTE = zeros(length(pz)-1) rTM = similar(rTE) - mintime = minimum(times) ## all this below for pesky zero time at start instead of shutoff ... + mintime = minimum(times) + if freqhigh < 3/mintime + freqhigh = 3/mintime + end @assert mintime > 0 # else interpolation in log10 trouble if doconvramp + checkrampformintime(times, ramp, minresptime) if isapprox(ramp[1,1], 0, atol=1e-8) mintime = ramp[2,1] # should be around 1.28e-5 to 2.5e-5 for HeliTEM - # @warn "automatically setting impulse response mintime to $mintime" + # @warn "automatically setting step response mintime to $mintime" end + mintime = 10^(log10(mintime) - 1) # go back a decade further than asked for + mintime = max(mintime, minresptime) # we can't handle responses shorter than minresptime -- unstable IDFT end - if freqhigh < 3/mintime - freqhigh = 3/mintime - end + interptimes = 10 .^(log10(mintime):1/ntimesperdecade:maximum(log10.(times))+1) if freqlow > 0.33/maximum(times) freqlow = 0.33/maximum(times) end if isempty(freqs) - freqs = 10 .^(log10(freqlow):1/nfreqsperdecade:log10(freqhigh)) + freqs = 10 .^(log10(freqlow):1/nfreqsperdecade:log10(freqhigh)) end interpkᵣ = 10 .^range(minimum(log10.(Filter_base))-0.5, maximum(log10.(Filter_base))+0.5, length = nkᵣeval) J0_kernel_h, J1_kernel_h, J0_kernel_v, J1_kernel_v = map(x->zeros(ComplexF64, length(interpkᵣ), length(freqs)), 1:4) J01kernelhold = zeros(ComplexF64, length(Filter_base)) log10ω = log10.(2*pi*freqs) - interptimes = 10 .^(log10(mintime)-1:1/ntimesperdecade:maximum(log10.(times))+1) HFD_z = zeros(ComplexF64, length(freqs)) # space domain fields in freq HFD_r = zeros(ComplexF64, length(freqs)) # space domain fields in freq HFD_az = zeros(ComplexF64, length(freqs)) # space domain fields in freq @@ -171,7 +176,7 @@ function HFieldDHT(; # Jtemp = calcjacobian ? zeros(ComplexF64, nmax) : zeros(0) useprimary = modelprimary ? one(Float64) : zero(Float64) - HFieldDHT(thickness, pz, ϵᵢ, zintfc, rTE, rTM, zRx, zTx, rTx, rRx, freqs, times, ramp, log10ω, interptimes, + HFieldDHT(thickness, pz, ϵᵢ, zintfc, rTE, rTM, zRx, zTx, rTx, rRx, freqs, times, ramp, log10ω, interptimes, minresptime, HFD_z, HFD_r, HFD_az, HFD_z_interp, HFD_r_interp, HFD_az_interp, HTD_z_interp, HTD_r_interp, HTD_az_interp, dBzdt, dBrdt, dBazdt, J0_kernel_h, J1_kernel_h, J0_kernel_v, J1_kernel_v, J01kernelhold, lowpassfcs, quadnodes, quadweights, preallocate_ω_Hsc(interptimes, lowpassfcs)..., rxwithinloop, provideddt, doconvramp, useprimary, @@ -180,6 +185,26 @@ function HFieldDHT(; HFD_r_J, HTD_r_J_interp, dBrdt_J) end +function checkrampformintime(times, ramp, minresptime) + # this checks we don't have ultra small tobs - ramp_time_a + for itime = 1:length(times) + for iramp = 1:size(ramp,1)-1 + rta, rtb = ramp[iramp,1], ramp[iramp+1,1] + if rta >= times[itime] # geq instead of gt as we could have an unlcky time + break + end + if rtb > times[itime] # end in this interval + rtb = times[itime] + end + # just so we know, rta < rtb and rta < t so rta < rtb <= t + ta = times[itime]-rta + tb = max(times[itime]-rtb, minresptime) # rtb > rta, so make sure this is not zero because integ is in log10... + @assert ta>tb # else we're in trouble + end + end + nothing +end + #update geometry and dependent parameters - necessary for adjusting geometry #e.g. inverting nuisance parameters function update_ZR!(F::HFieldDHT, zTx, zRx, rTx, rRx; onlyprimary=false) @@ -641,7 +666,7 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu end ta = F.times[itime]-rta - tb = max(F.times[itime]-rtb, 1e-8) # rtb > rta, so make sure this is not zero because integ is in log10... + tb = max(F.times[itime]-rtb, F.minresptime)# rtb > rta, so make sure this is not zero because integ is in log10... a, b = log10(ta), log10(tb) x, w = F.quadnodes, F.quadweights F.dBzdt[itime] += (b-a)/2*dot(getrampresponse((b-a)/2*x .+ (a+b)/2, splz), w)*dIdt diff --git a/src/TEMPEST1DInversion.jl b/src/TEMPEST1DInversion.jl index bce1e390..2341b9dd 100644 --- a/src/TEMPEST1DInversion.jl +++ b/src/TEMPEST1DInversion.jl @@ -111,7 +111,7 @@ function Bfield(; F = AEM_VMD_HMD.HFieldDHT(;nmax, calcjacobian, ntimesperdecade = ntimesperdecade, nfreqsperdecade = nfreqsperdecade, - freqlow=1e-5, + # freqlow =1e-5, nkᵣeval = nkᵣeval, times = times, ramp = ramp, diff --git a/src/VTEM1DInversion.jl b/src/VTEM1DInversion.jl index 0d3c09e8..acce20af 100644 --- a/src/VTEM1DInversion.jl +++ b/src/VTEM1DInversion.jl @@ -48,8 +48,8 @@ function dBzdt(;times = [1.], showgates = false, modelprimary = false, lowpassfcs = [], - freqlow = 1e-3, - freqhigh = 5e5, + freqlow = 1e-4, + freqhigh = 1e6, ) @assert size(σ) == size(d) diff --git a/test/verify_skytem_td.jl b/test/verify_skytem_td.jl index 268fbb39..b15b1b8f 100644 --- a/test/verify_skytem_td.jl +++ b/test/verify_skytem_td.jl @@ -3,7 +3,7 @@ include("skytem_response.jl") ## modeling parameters ntimesperdecade = 10 nfreqsperdecade = 5 -freqlow = 1e-3 +freqlow = 1e-4 ## LM operator Flm = transD_GP.AEM_VMD_HMD.HFieldDHT( ntimesperdecade = ntimesperdecade, From d4780356a53672065ef8eba350c06a85742e0ed0 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Fri, 27 Oct 2023 20:16:54 +1100 Subject: [PATCH 04/10] I think fixed the early time stuff for ontime gates --- src/AEM_VMD_HMD.jl | 44 ++++++++++++++++++++++++++-------------- test/verify_skytem_td.jl | 2 +- test/verify_vmd_td.jl | 37 +++++++++++++++++---------------- 3 files changed, 49 insertions(+), 34 deletions(-) diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl index cea6f9e9..0cab257e 100644 --- a/src/AEM_VMD_HMD.jl +++ b/src/AEM_VMD_HMD.jl @@ -80,7 +80,7 @@ function HFieldDHT(; zTx = -35.0, zRx = -37.5, times = 10 .^LinRange(-6, -1, 50), - ramp = ones(10, 10), + ramp = ones(10, 2), nfreqsperdecade = 5, ntimesperdecade = 10, glegintegorder = 5, @@ -107,21 +107,17 @@ function HFieldDHT(; rTM = similar(rTE) ## all this below for pesky zero time at start instead of shutoff ... - mintime = minimum(times) - if freqhigh < 3/mintime - freqhigh = 3/mintime - end + mintime, maxtime = extrema(times) @assert mintime > 0 # else interpolation in log10 trouble + @show mintime = 10^(log10(mintime) - 1) # go back a decade further than asked for + @show maxtime = 10^(log10(maxtime) + 1) # go ahead a decade further if doconvramp - checkrampformintime(times, ramp, minresptime) - if isapprox(ramp[1,1], 0, atol=1e-8) - mintime = ramp[2,1] # should be around 1.28e-5 to 2.5e-5 for HeliTEM - # @warn "automatically setting step response mintime to $mintime" - end - mintime = 10^(log10(mintime) - 1) # go back a decade further than asked for - mintime = max(mintime, minresptime) # we can't handle responses shorter than minresptime -- unstable IDFT + mintime, maxtime = checkrampformintime(times, ramp, minresptime, mintime, maxtime) + end + interptimes = 10 .^(log10(mintime) : 1/ntimesperdecade : log10(maxtime)) + if freqhigh < 3/mintime + freqhigh = 3/mintime end - interptimes = 10 .^(log10(mintime):1/ntimesperdecade:maximum(log10.(times))+1) if freqlow > 0.33/maximum(times) freqlow = 0.33/maximum(times) end @@ -185,8 +181,10 @@ function HFieldDHT(; HFD_r_J, HTD_r_J_interp, dBrdt_J) end -function checkrampformintime(times, ramp, minresptime) +function checkrampformintime(times, ramp, minresptime, mintime, maxtime) # this checks we don't have ultra small tobs - ramp_time_a + minta = Inf + maxta = -Inf for itime = 1:length(times) for iramp = 1:size(ramp,1)-1 rta, rtb = ramp[iramp,1], ramp[iramp+1,1] @@ -200,9 +198,25 @@ function checkrampformintime(times, ramp, minresptime) ta = times[itime]-rta tb = max(times[itime]-rtb, minresptime) # rtb > rta, so make sure this is not zero because integ is in log10... @assert ta>tb # else we're in trouble + if ta < minta + minta = ta + end + if ta > maxta + maxta = ta + end end + # if minta < mintime + mintime = max(0.5minta, minresptime) + # though I believe mintime = minta always is safe. + # end + # if maxta > maxtime + # too high a time requires very low freqs + # maxtime = maxta # should work but doesn't with SkyTEM / VTEM checks + maxtime = min(3maxta, maxtime) + # end end - nothing + # @info minta, maxta + mintime, maxtime end #update geometry and dependent parameters - necessary for adjusting geometry diff --git a/test/verify_skytem_td.jl b/test/verify_skytem_td.jl index b15b1b8f..391f79ec 100644 --- a/test/verify_skytem_td.jl +++ b/test/verify_skytem_td.jl @@ -16,7 +16,7 @@ Flm = transD_GP.AEM_VMD_HMD.HFieldDHT( rTx = rTx, zRx = zRxLM, freqlow = freqlow - ) + ); ## HM operator Fhm = transD_GP.AEM_VMD_HMD.HFieldDHT( ntimesperdecade = ntimesperdecade, diff --git a/test/verify_vmd_td.jl b/test/verify_vmd_td.jl index 7b0cdee0..a8c5dc31 100644 --- a/test/verify_vmd_td.jl +++ b/test/verify_vmd_td.jl @@ -1,4 +1,4 @@ -using PyPlot, Revise, HiQGA.transD_GP, Random +using PyPlot, HiQGA.transD_GP, Random ## set up nfreqsperdecade = 20 ntimesperdecade = 20 @@ -8,30 +8,31 @@ doconvramp = false nkᵣeval = 200 times = 10 .^LinRange(-5,-1.1, 60) lowpassfcs = [1e6] -freqhigh = 1e6 -freqlow = 1e-8 +freqhigh = 1e5 +freqlow = 1e-7 ## model zfixed = [-1e5, 0] rho = [1e12, 100] nmax = 200 ## geometry rRx = 100. -zRx = 0. +zRx = -0.01 zTx = 0. -modelprimary = true ## -F = transD_GP.AEM_VMD_HMD.HFieldDHT(freqlow = freqlow, - freqhigh = freqhigh, - zTx = zTx, - rRx = rRx, - times = times, - nfreqsperdecade = nfreqsperdecade, - ntimesperdecade = ntimesperdecade, - zRx = zRx, - nkᵣeval = nkᵣeval, - modelprimary = modelprimary, - provideddt = provideddt, - lowpassfcs = lowpassfcs) +F = transD_GP.AEM_VMD_HMD.HFieldDHT(; + freqlow = freqlow, + freqhigh = freqhigh, + zTx = zTx, + rRx = rRx, + times = times, + doconvramp, + nfreqsperdecade = nfreqsperdecade, + ntimesperdecade = ntimesperdecade, + zRx = zRx, + nkᵣeval = nkᵣeval, + modelprimary = modelprimary, + provideddt = provideddt, + lowpassfcs = lowpassfcs); transD_GP.AEM_VMD_HMD.getfieldTD!(F, zfixed, rho) ## plot FD figure() @@ -60,7 +61,7 @@ title("Compare with W&H Fig 4.4") plt.tight_layout() ## Radial fields figure(figsize=(4,7)) -F.useprimary = 0. +# F.useprimary = 0. F.getradialH = true F.provideddt = true transD_GP.AEM_VMD_HMD.getfieldTD!(F, zfixed, rho) From 2cbd5ba550e49375ffd91c7ebf78f14dcab73da8 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Fri, 27 Oct 2023 20:19:54 +1100 Subject: [PATCH 05/10] remove debugs --- src/AEM_VMD_HMD.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl index 0cab257e..b68cbe5b 100644 --- a/src/AEM_VMD_HMD.jl +++ b/src/AEM_VMD_HMD.jl @@ -109,8 +109,8 @@ function HFieldDHT(; ## all this below for pesky zero time at start instead of shutoff ... mintime, maxtime = extrema(times) @assert mintime > 0 # else interpolation in log10 trouble - @show mintime = 10^(log10(mintime) - 1) # go back a decade further than asked for - @show maxtime = 10^(log10(maxtime) + 1) # go ahead a decade further + mintime = 10^(log10(mintime) - 1) # go back a decade further than asked for + maxtime = 10^(log10(maxtime) + 1) # go ahead a decade further if doconvramp mintime, maxtime = checkrampformintime(times, ramp, minresptime, mintime, maxtime) end From 904b351ae4d7ad24c19daf68602facc42a1e9c4e Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Fri, 27 Oct 2023 23:04:37 +1100 Subject: [PATCH 06/10] mintimes, maxtimes solidly fixed .. -needed to move maxtimes = 1.5maxta to ensure that we reach maxta in interptimes -also cleaner to do min, max outside loop in check --- src/AEM_VMD_HMD.jl | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl index b68cbe5b..5dcf00a3 100644 --- a/src/AEM_VMD_HMD.jl +++ b/src/AEM_VMD_HMD.jl @@ -205,17 +205,15 @@ function checkrampformintime(times, ramp, minresptime, mintime, maxtime) maxta = ta end end - # if minta < mintime - mintime = max(0.5minta, minresptime) - # though I believe mintime = minta always is safe. - # end - # if maxta > maxtime - # too high a time requires very low freqs - # maxtime = maxta # should work but doesn't with SkyTEM / VTEM checks - maxtime = min(3maxta, maxtime) - # end end - # @info minta, maxta + # if minta < mintime + mintime = max(0.5minta, minresptime) + # though I believe mintime = minta always is safe. + # end + # if maxta > maxtime + maxtime = min(1.5maxta, maxtime) + # 1.5 to make sure we clear maxta in the interptimes + # end mintime, maxtime end From 4f78d7b18b3d92769465deaaa3f6d916820ad4f7 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Sat, 28 Oct 2023 20:03:38 +1100 Subject: [PATCH 07/10] refreshing some old tests --- test/verify_vmd_loop_fd.jl | 13 +++++++++---- test/verify_vmd_td.jl | 10 +++++----- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/test/verify_vmd_loop_fd.jl b/test/verify_vmd_loop_fd.jl index 33dbc310..54a471ed 100644 --- a/test/verify_vmd_loop_fd.jl +++ b/test/verify_vmd_loop_fd.jl @@ -13,12 +13,14 @@ rRx = 100. zRx = -0.011 zTx = -0.01 nkᵣeval = 200 +doconvramp = false ## VMD modelprimary = true -Fvmd = transD_GP.AEM_VMD_HMD.HFieldDHT( +Fvmd = transD_GP.AEM_VMD_HMD.HFieldDHT(; nmax = nmax, zTx = zTx, rRx = rRx, + doconvramp, freqs = freqs, zRx = zRx, nkᵣeval = nkᵣeval, @@ -26,12 +28,13 @@ Fvmd = transD_GP.AEM_VMD_HMD.HFieldDHT( transD_GP.AEM_VMD_HMD.getfieldFD!(Fvmd, zfixed, rho) ## now use a tiny loop radius - VMD approximation is worse as radius gets larger rTx = 0.001 -Floop = transD_GP.AEM_VMD_HMD.HFieldDHT( +Floop = transD_GP.AEM_VMD_HMD.HFieldDHT(; nmax = nmax, zTx = zTx, rRx = rRx, rTx = rTx, freqs = freqs, + doconvramp, zRx = zRx, nkᵣeval = nkᵣeval, modelprimary = modelprimary) @@ -46,12 +49,13 @@ ax[1].legend() ax[1].set_title("Compare with W&H Fig 4.2") ## Compare H radial getradialH = true -Fvmd = transD_GP.AEM_VMD_HMD.HFieldDHT( +Fvmd = transD_GP.AEM_VMD_HMD.HFieldDHT(; nmax = nmax, zTx = 0, rRx = rRx, freqs = freqs, zRx = 0, + doconvramp, nkᵣeval = nkᵣeval, modelprimary = false, # says W&H for within plane of observation !! getradialH = getradialH) @@ -70,13 +74,14 @@ zTx = -0.01 zRx = -0.02 rRx = 0.001 rTx = 50.0 -Floopin = transD_GP.AEM_VMD_HMD.HFieldDHT( +Floopin = transD_GP.AEM_VMD_HMD.HFieldDHT(; nmax = nmax, zTx = zTx, rRx = rRx, rTx = rTx, freqs = freqs, zRx = zRx, + doconvramp, nkᵣeval = nkᵣeval, modelprimary = modelprimary) transD_GP.AEM_VMD_HMD.getfieldFD!(Floopin, zfixed, rho) diff --git a/test/verify_vmd_td.jl b/test/verify_vmd_td.jl index a8c5dc31..2a82483f 100644 --- a/test/verify_vmd_td.jl +++ b/test/verify_vmd_td.jl @@ -1,13 +1,13 @@ using PyPlot, HiQGA.transD_GP, Random ## set up -nfreqsperdecade = 20 -ntimesperdecade = 20 +nfreqsperdecade = 10 +ntimesperdecade = 10 # only for interpolation in plots modelprimary = true provideddt = true doconvramp = false -nkᵣeval = 200 -times = 10 .^LinRange(-5,-1.1, 60) -lowpassfcs = [1e6] +nkᵣeval = 120 # for HT +times = 10 .^LinRange(-5,-1.1, 5) +lowpassfcs = [] freqhigh = 1e5 freqlow = 1e-7 ## model From 203423a62bbb032516bd906c4a8d4bfc80965022 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Sat, 28 Oct 2023 20:22:42 +1100 Subject: [PATCH 08/10] refreshing some old tests --- test/verify_loop_td.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/verify_loop_td.jl b/test/verify_loop_td.jl index 4b459405..fe34a74f 100644 --- a/test/verify_loop_td.jl +++ b/test/verify_loop_td.jl @@ -5,15 +5,15 @@ rho = [1e20, 1, ] nmax = 200 ## geometry rTx = 15 -rRx = 0.001 -zRx = -0.001 -zTx = -0.002 +rRx = 0.0 +zRx = -0.01 +zTx = -0.005 doconvramp = false modelprimary = true nkᵣeval = 200 ntimesperdecade = 15 nfreqsperdecade = 15 -times = 10 .^LinRange(-6,1,40) +times = 10 .^LinRange(-5,0,5) ## build operator F = transD_GP.AEM_VMD_HMD.HFieldDHT( ntimesperdecade = ntimesperdecade, From e7305d16ee773d132617e5a08d156c21bd04b102 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Sat, 28 Oct 2023 21:17:14 +1100 Subject: [PATCH 09/10] test refreshed, mintime in convcheck removed input --- src/AEM_VMD_HMD.jl | 6 ++---- test/verify_tempest_geometry.jl | 4 ++-- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl index 5dcf00a3..3d229fbf 100644 --- a/src/AEM_VMD_HMD.jl +++ b/src/AEM_VMD_HMD.jl @@ -112,7 +112,7 @@ function HFieldDHT(; mintime = 10^(log10(mintime) - 1) # go back a decade further than asked for maxtime = 10^(log10(maxtime) + 1) # go ahead a decade further if doconvramp - mintime, maxtime = checkrampformintime(times, ramp, minresptime, mintime, maxtime) + mintime, maxtime = checkrampformintime(times, ramp, minresptime, maxtime) end interptimes = 10 .^(log10(mintime) : 1/ntimesperdecade : log10(maxtime)) if freqhigh < 3/mintime @@ -181,7 +181,7 @@ function HFieldDHT(; HFD_r_J, HTD_r_J_interp, dBrdt_J) end -function checkrampformintime(times, ramp, minresptime, mintime, maxtime) +function checkrampformintime(times, ramp, minresptime, maxtime) # this checks we don't have ultra small tobs - ramp_time_a minta = Inf maxta = -Inf @@ -206,10 +206,8 @@ function checkrampformintime(times, ramp, minresptime, mintime, maxtime) end end end - # if minta < mintime mintime = max(0.5minta, minresptime) # though I believe mintime = minta always is safe. - # end # if maxta > maxtime maxtime = min(1.5maxta, maxtime) # 1.5 to make sure we clear maxta in the interptimes diff --git a/test/verify_tempest_geometry.jl b/test/verify_tempest_geometry.jl index 8ecb2113..9c56cde8 100644 --- a/test/verify_tempest_geometry.jl +++ b/test/verify_tempest_geometry.jl @@ -8,7 +8,7 @@ nmax = 100 zstart = 0. extendfrac, dz = 1.06, 2. -zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=40, showplot=true) +zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=40, showplot=false) z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed) @info z @@ -66,7 +66,7 @@ tempest = transD_GP.TEMPEST1DInversion.Bfield( addprimary = true #this ensures that the geometry update actually changes everything that needs to be ) # plot before testing -transD_GP.TEMPEST1DInversion.plotmodelfield!(tempest,z,ρ) +transD_GP.TEMPEST1DInversion.plotmodelfield!(tempest,log10.(ρ[2:end])) μ = transD_GP.TEMPEST1DInversion.μ₀ Bx, By, Bz = tempest.Hx[1]*μ*1e15, tempest.Hy[1]*μ*1e15, tempest.Hz[1]*μ*1e15 From 1ef305d50c23bd0fdae1f99ec775747cb2a23b41 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Mon, 30 Oct 2023 11:08:52 +1100 Subject: [PATCH 10/10] version bump --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e1484fd7..57a5a0b6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "HiQGA" uuid = "01574e87-63b6-4466-925a-334cf7e21b6b" authors = ["richardt94 <29562583+richardt94@users.noreply.github.com>"] -version = "0.3.12" +version = "0.3.13" [deps] Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"