math part of b and dbdt corrected

[MasayukiMotoori]

Hi

I corrected potential bag about b and dbdt field calculation about whole space due to E-dipole.

Math for B : $\mathbf{h}(t) = \frac{Ids}{4 \pi r^3} \bigg ( \textrm{erf}(\theta r) - \frac{2}{\sqrt{\pi}} \theta r \, e^{-\theta^2 r^2} \bigg ) \big ( - z \ \mathbf{\hat y} + y \ \mathbf{\hat z} \big )$ Implemented:
$\mathbf{h}(t) = \frac{Ids}{4 \pi r^3} \bigg ( \textrm{erf}(\theta r) + \frac{2}{\sqrt{\pi}} \theta r \, e^{-\theta^2 r^2} \bigg ) \big ( - z \ \mathbf{\hat y} + y \ \mathbf{\hat z} \big )$

Math for dbdt: $\frac{\partial \mathbf{h}}{\partial t} = - \frac{2 \, \theta^5 Ids}{\pi^{3/2} \mu \sigma} e^{-\theta^2 r^2} \big ( - z \ \mathbf{\hat y} + y \ \mathbf{\hat z} \big )$ Implemented : $\frac{\partial \mathbf{h}}{\partial t} = - \frac{2 \, \theta^5 Ids}{\pi^{3/2} \mu \sigma} \big ( - z \ \mathbf{\hat y} + y \ \mathbf{\hat z} \big )$

I forked the repo and corrected math part, installed in environment and test it. I uploaded the note book about it in my git hub as well. https://github.com/MasayukiMotoori/EOSC556B/blob/main/90_Geoana_Check_revised.ipynb

[prisae]

Good catch!

I tried to re-run your notebook that you linked (https://github.com/MasayukiMotoori/EOSC556B/blob/main/90_Geoana_Check_revised.ipynb), but it gives me a 404 Error - maybe it is private?

I wanted to quick run it myself, and also check if, with some setting, I could improve the early empymod part.

[MasayukiMotoori]

Hi @prisae Thank you for your comment. I made the repository public, and it should be visible now.

[prisae]

I had a look at your example. Very early times is a challenge for the Fourier-Transform, and out of interest I looked at how much I can squeeze it. As you told me that you use empymod a lot, I post it here, as it might be useful for other computations too.

The main points (also as comments in the code):

• You can straight define a fullspace, no need to create a "pseudo"-fullspace with a halfspace, with source and receivers far away from the interface. This has the advantage that it can use the analytical solution in the frequency domain internally (see next point.)
• Setting xdirect=True can help, as the direct field is then computed analytically in the f-domain; in the case of a fullspace, the entire f-domain solution is computed analytically in the f-domain (rather than analytically in the wavenumber-domain, and transformed to the Fourier domain).
• As I mentioned already, very early times is hard for the Fourier Transform; sometimes you can achieve better results with FFTLog than with DLF. See https://empymod.emsig.xyz/en/stable/gallery/tdomain/step_and_impulse.html for a comparison of DLF, FFT, FFTLog, and QWE for time-domain responses.

I hope that helps! Feel free to reach out here or on Mattermost if there are things which I did not explain properly, or any other question you might have.

import empymod
import numpy as np
from scipy import special
import matplotlib.pyplot as plt
from scipy.constants import mu_0

def tem_edp_whl(sigma, t, mu, rec):
    """
    # whole space, transient, E-dipole.
    # E-dipole is in x direction
    """
    θ = np.sqrt((mu*sigma)/(4*t))
    r = np.linalg.norm(rec, 2)

    h = np.zeros((np.shape(t)[0], 3))
    hamp = 2/np.sqrt(np.pi)*θ*r*np.exp(-(r**2)*(θ**2)) + special.erfc(θ*r)
    hamp = 1-hamp
    hamp /= (4*np.pi*(r**2))
    h[:, 1] = -rec[2]/r*hamp
    h[:, 2] = -rec[1]/r*hamp

    dhdt = np.zeros((np.shape(t)[0], 3))
    dhamp = (θ**3)*r/(2*t*np.pi**1.5)*np.exp(-(θ**2)*(r**2))
    dhdt[:, 1] = -rec[2]/r*dhamp
    dhdt[:, 2] = -rec[1]/r*dhamp
    return h, dhdt

# Parameters
sigma = 3
off = 1.75
t = np.logspace(-7, -2, 51)

# empymod computation
model = {
    'src': [0., 0., 0.],
    'rec': [off, 0., 0.],
    # You can straight define a fullspace, no need to create a "pseudo"-fullspace
    # with a halfspace, with source and receivers far away from the interface.
    'depth': [],
    'res': 1/sigma,
    'freqtime': t,
    'verb': 1,
    'ab': 62,
    # Very early times is hard for the Fourier Transform; sometimes you can achieve better results
    # with FFTLog than with DLF. See https://empymod.emsig.xyz/en/stable/gallery/tdomain/step_and_impulse.html
    # for more examples
    'ft': 'fftlog',
    'ftarg': {'add_dec': [-5, 5], 'pts_per_dec': 20},
    # Setting xdirect=True can help, as direct field is then computed analytically in f-domain
    'xdirect': True,
}

empymod_b = -mu_0*empymod.dipole(signal=-1, **model)
empymod_dbdt = -mu_0*empymod.dipole(signal=0, **model)

# Hohmann computation
hohmann_b_tmp, hohmann_dbdt_tmp = tem_edp_whl(sigma, t, mu_0, [0., -off, 0.])
hohmann_b = mu_0*hohmann_b_tmp[:, 2]
hohmann_dbdt = mu_0*hohmann_dbdt_tmp[:, 2]

# Figure
fig, ax = plt.subplots(2, 1, sharex=True, constrained_layout=True)

ax[0].loglog(t, hohmann_b, "C3", label="Hohmann")
ax[0].loglog(t, empymod_b, "k--", label="empymod")
ax[0].set_title("b field ")

ax[1].loglog(t, hohmann_dbdt, "C3", label="Hohmann")
ax[1].loglog(t, empymod_dbdt, "k--", label="empymod")
ax[1].set_title("dbdt field")
ax[1].set_xlabel("time(sec)")

for a in ax:
    a.grid()
    a.legend()

[prisae]

Btw, to see the similarities to https://empymod.emsig.xyz/en/stable/gallery/tdomain/step_and_impulse.html better, you can simply change your loglog-commands to semilogx, which results in

[MasayukiMotoori]

@prisae Thank you for your looking at my notebook and showing me how to use empymod. Actually I used geoana to validate empymod before chargeability simulation. I didn't know how to use these empymod parameters. It is helpful for my research,

[lheagy]

Thanks so much for your contribution @MasayukiMotoori! And so cool to see the conversation here on simulating those early times @prisae -- thank you!!

[jcapriot]

Hey Thanks for catching this as well! It all looks good on my end now that you've corrected it, correctly matching equation 2.51 in ward in hohmann for the H-field, and dH/dT as well.

Looking even further through the code, the E-field seems off as well, I'll probably make a PR myself to get it all in shape! (Also I'll get the CI back up and running.)