Re-normalize the frequency-domain response in case of a magnetic source #45

ruboerner · 2019-09-16T07:33:08Z

With these small changes it seems possible to calculate the time-domain response of a layered half-space due to a magnetic dipole source. The reason for the code change is that prior to the f->t transform, the normalization (by i \omega \mu_0) has to be undone.

…ce prior to transform to time-domain

prisae · 2019-09-16T09:50:12Z

Thanks Ralph-Uwe for your contribution! Improvements to the magnetic-side of things are very welcome, as I hardly ever use that part of the code.

Couple of questions:

Just in case: Do you know of any analytical time-domain solution we could use for comparison/tests? The analytical solutions currently implemented are only in the frequency domain (for magnetic sources).
I first thought instantly that this should live in model.tem. But model.tem is not aware of the source-type, so you put it into the right place.

Two tests fail currently. The first one is a regression test, so expected to fail. The second one, however, should not fail. I will have to look into that in more detail.

I'll leave some comments in the code too. If you could change that great. Or, if you don't mind, then I'll do it when I get around to adjust the tests (probably after my holidays, so give me two weeks....)

prisae

Thanks for the PR! As I already mentioned, if you know of an analytical solution that we could use for comparison that would be great!

prisae · 2019-09-16T09:58:48Z

empymod/model.py

+            # back into time-domain has to be carried out.
+            for kk in range(len(freq)):
+                EM[kk, :] *= 2j * np.pi * freq[kk] * np.pi * 4e-7
+


Would you mind changing

for kk in range(len(freq)): EM[kk, :] *= 2j * np.pi * freq[kk] * np.pi * 4e-7

to the more pythonic

for kk, kfreq in enumerate(freq): EM[kk, :] *= 2j * np.pi * kfreq * np.pi * 4e-7

here and in the other case too?

Thinking about it I think we should actually avoid the loop entirely:

EM *= 2j * np.pi * freq[:, None] * np.pi * 4e-7

(Relevant codacy feedback: https://app.codacy.com/manual/prisae/empymod/pullRequest?prid=4165816)

That's perfectly ok!

…me domain

ruboerner · 2019-09-16T11:24:20Z

A simple function for checking the impulse response due to a vertical magnetic dipole is:
########################################
import numpy as np
from scipy.constants import mu_0
from scipy.special import erf

def dhzdt(t, sigma, r):
r"""Return the impulse response (i.e., dHz(t)/dt) of a homogeneous isotropic half-space
of electrical conductivity sigma in z>=0 due to a current shut-off in a vertical
magnetic dipole source of unit moment for times t>0.

Source and receiver are both at the plane z=0 and separated by the distance r.
"""
theta = np.sqrt(mu_0 * sigma / 4 / t)
theta_r = theta * r
s = 9 * erf(theta_r)
s -= 2 * theta_r /np.sqrt(np.pi) * (9 + 6 * theta_r**2 + 4 * theta_r**4) * np.exp(-theta_r**2)
s /= (-2 * np.pi * mu_0 * sigma * r**5)
return s

#################################################

prisae · 2019-09-16T12:19:15Z

Fantastic, thanks Ralph-Uwe! Do you have a reference for the analytical solution?

I might get time to include it tomorrow and add the relevant tests, but more likely it will have to wait two weeks until I am back. Either way, I will make a new release with the fix pushing it to v1.9.0, including other recent changes I made.

Thanks again for the fix. And a big hooray, you are actually the first contributor (as far as code goes) to empymod! I will include you in the CREDITS file, and probably change the (c) from my name to something more generic like "(c) The empymod developers".

I hope you found empymod useful (despite the bug ;-) ).

ruboerner · 2019-09-16T12:25:55Z

Thanks Dieter,

Kudos for the empymod package!

Here's the reference:

@incollection{ward1988electromagnetic,
title={Electromagnetic theory for geophysical applications},
author={Ward, Stanley H and Hohmann, Gerald W},
booktitle={Electromagnetic Methods in Applied Geophysics: Voume 1, Theory},
pages={130--311},
year={1988},
publisher={Society of Exploration Geophysicists}
}

...have a nice holiday!

prisae · 2019-09-17T15:09:56Z

OK. I thought about it again, and talked it through with Evert Slob, who was one of the co-others of the paper on which the code is based (Hunziker et al., 2015, Geophysics).

So currently it returns the H-field, the field due to a magnetic dipole. Your PR would turn that into a B-field for a coil, which turns your PR in a feature request rather than a bug fix. So it would have to be implemented optionally, not exclusively. And also, if we implement it for source coils we should implement it for receiver coils too.

I'll have a think about it. Maybe an additional input parameter for dipole and bipole. Or a new function coil in addition to dipole and bipole.

Inputs/ideas/opinions welcome!

prisae · 2019-09-18T06:08:55Z

One possibility is to have for msrc/mrec in addition to False/True a third state, 'loop'. That would only work for bipole though, not for dipole, but I think that is fine. bipole is the more flexible function that can also do everything which dipole can...

Do you mind if I work directly on your PR, @ruboerner ?

ruboerner · 2019-09-18T06:20:53Z

One possibility is to have for msrc/mrec in addition to False/True a third state, 'loop'. That would only work for bipole though, not for dipole, but I think that is fine. bipole is the more flexible function that can also do everything which dipole can...

I think it might be better to introduce a new function for magnetic loops as "bipole" suggests that electric bipoles are used. This might just be a stripped-down version of bipole without e sources?

Do you mind if I work directly on your PR, @ruboerner ?

I'm fine with that.

prisae · 2019-09-18T06:26:35Z

Sounds good. You think just loop-loop, or both, loop-loop and loop-dipole? As far as I know, in reality the magnetic source is always a loop, whereas the receiver can be both.

As I mentioned, I'll be off from today for two weeks, but will have a look at this afterwards. Is it OK if I bother you for some testing/feedback once I have something? I personally have not much experience with magnetic sources, so it would be great and highly appreciated to have your input!

ruboerner · 2019-09-18T06:32:37Z

Sounds good. You think just loop-loop, or both, loop-loop and loop-dipole? As far as I know, in reality the magnetic source is always a loop, whereas the receiver can be both.

Yes, most land-based TEM systems use loops as sources and coils or small loops as receivers, which is exactly what you suggest.

As I mentioned, I'll be off from today for two weeks, but will have a look at this afterwards. Is it OK if I bother you for some testing/feedback once I have something? I personally have not much experience with magnetic sources, so it would be great and highly appreciated to have your input!

Yes, don't worry, it's absolutely ok. I'm glad to help.

Have a nice time off!

prisae · 2019-10-04T12:02:08Z

I hope to draft it early next week. I though of using bipole as a starting point, but with the following simplifications:

Magnetic loop sources, with azimuth and dip; unit length plus strength.
Magnetic loop and dipole receivers, with azimuth and dip; unit length.
Only time (or also frequency)?

Does that make sense?

Any clever idea for the function name? Simply loop comes to my mind, but there might be a better alternative.

prisae · 2019-10-10T13:51:14Z

It is necessary for frequency and time domain modelling, so I moved it to fem. Have a look at the new example here:
https://github.com/empymod/empymod-examples/blob/master/5b_TEM.ipynb
where I reproduced the results by Ward and Hohmann, 1988, Figures 4.2-4.5.

Now, given that, I actually think it is a bug in empymod.kernel. \zeta = s\mu = i\omega\mu; there must be a factor \zeta missing somewhere for magnetic sources...

I'll have to dig into that in more detail, so don't expect it to resolve anytime soon. I don't want to include a workaround in model that fixes a bug in kernel.

coveralls · 2019-10-10T13:52:04Z

Coverage decreased (-1.6%) to 98.38% when pulling 03d769f on ruboerner:master into 7f65e3a on empymod:master.

prisae · 2019-10-10T14:55:45Z

Note to self to not forget: it should not be mu_0, but the mu of the layer in which the source resides, and then differently for the different wavenumber-kernels depending if affected by vertical or horizontal mu, mu_v or mu_h (hence use mu=mu_0*mu_{r;s}^h or mu=mu_0*mu_{r;s}^v). Of course, if the source is in the air that all simplifies to mu_0 again. But the current implementation only holds for that special case, make it more general.

prisae · 2019-10-10T18:30:40Z

I am turning in circles. I take back what I said before. Again, it is the difference between a magnetic dipole and a loop source. So not a bug. Just confusing for someone never actually involved with magnetic sources. I am still confused with the units though... I will discuss this through and come back with a solution.

prisae · 2019-10-11T11:46:46Z

@ruboerner , do you know of any analytical solutions for magnetic source with an electric receiver? What leaves me a bit uneasy with implementing a fix for if msrc is that this also affects magnetic-source to electric-receiver cases, and I am not sure if they are affected in the same way...

prisae · 2019-10-11T13:25:39Z

Comparison with simpegEM1D confirms the problem...

prisae · 2019-10-14T15:59:15Z

Extensive discussion, amongst others with the authors of the original article and their code EMmod, confirm that everything is fine with the current implementation. All sources and receivers are assumed to be dipoles.

This changes for loops, and the solution therefore is really to have a loop-function. If the source is a loop, it needs the factor j\omega\mu; if the receiver is a loop, it needs another j\omega\mu.

prisae · 2019-10-15T14:31:18Z

Hi @ruboerner - This took more effort than I anticipated, but it is in. There is a new version empymod v1.10.0, which has a new function model.loop. It takes the same arguments as model.bipole, except that

mrec can be [True, False, 'loop']`.

Furthermore, the following parameters are removed:

msrc=True: fixed, cannot be changed;
srcpts=1: fixed, cannot be changed.

This means that it is a loop source, measured by

an electric or magnetic, arbitrary rotated, finite length bipole, or
a loop

In the end I worked on a new branch, as everything changed a bit (see #48). I will therefore close this PR. I mentioned your help in the Credits-file.

I would appreciate if you could test the new version and let me know if it works for you!

Thanks.

prisae · 2019-10-19T10:25:05Z

Here a graphical summary of it, with explanations:
https://nbviewer.jupyter.org/github/empymod/empymod-examples/blob/master/7b_Dipoles_and_Loops.ipynb

ruboerner · 2019-10-22T07:07:14Z

Hi @prisae , sorry for the late response. Big kudos for your effort you invested in implementing magnetic sources. This is indeed a huge improvement!

prisae · 2019-10-22T07:11:18Z

Thanks @ruboerner . Did you test it, does it work work for you as expected?

There is obviously much more that could be done in empymod for TEM-style things. Like allowing for more complicated waveforms, windowing as often used in TEM's, etc. But this is far off my current topic... But I'd be happy to include more things or assist in implementing them, if there is interest.

prisae · 2019-11-06T10:07:18Z

Instead of using empymod.loop, we can also actually model an electric loop.

I added examples here:
https://github.com/empymod/empymod-examples/blob/master/7b_Dipoles_and_Loops.ipynb
modelling a loop with either four el. point sources, or four electric finite length dipole sources.

prisae · 2019-11-22T14:57:05Z

New example which goes into the same direction (the above mentioned examples are in the new sphinx-gallery too):

Re-normalize the frequency-domain response in case of a magnetic sour…

96f46de

…ce prior to transform to time-domain

prisae requested changes Sep 16, 2019

View reviewed changes

Add an analytical solution for the vertical magnetic dipole in the ti…

f0c01f8

…me domain

prisae added the enhancement New feature or request label Sep 17, 2019

prisae added 2 commits October 7, 2019 17:21

Initiate loop, not yet thought through

52df94c

Move to fem; I think it is a bug in kernel

03d769f

prisae mentioned this pull request Oct 10, 2019

Add possibility to have loops as source and receiver. #46

Closed

prisae added bug Something isn't working and removed bug Something isn't working labels Oct 11, 2019

prisae closed this Oct 15, 2019

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Re-normalize the frequency-domain response in case of a magnetic source #45

Re-normalize the frequency-domain response in case of a magnetic source #45

ruboerner commented Sep 16, 2019

prisae commented Sep 16, 2019

prisae left a comment

prisae Sep 16, 2019

prisae Sep 16, 2019

ruboerner Sep 16, 2019

ruboerner commented Sep 16, 2019

prisae commented Sep 16, 2019

ruboerner commented Sep 16, 2019 •

edited

Loading

prisae commented Sep 17, 2019

prisae commented Sep 18, 2019

ruboerner commented Sep 18, 2019

prisae commented Sep 18, 2019

ruboerner commented Sep 18, 2019

prisae commented Oct 4, 2019

prisae commented Oct 10, 2019

coveralls commented Oct 10, 2019

prisae commented Oct 10, 2019

prisae commented Oct 10, 2019

prisae commented Oct 11, 2019

prisae commented Oct 11, 2019

prisae commented Oct 14, 2019

prisae commented Oct 15, 2019

prisae commented Oct 19, 2019

ruboerner commented Oct 22, 2019

prisae commented Oct 22, 2019

prisae commented Nov 6, 2019

prisae commented Nov 22, 2019 •

edited

Loading

Re-normalize the frequency-domain response in case of a magnetic source #45

Re-normalize the frequency-domain response in case of a magnetic source #45

Conversation

ruboerner commented Sep 16, 2019

prisae commented Sep 16, 2019

prisae left a comment

Choose a reason for hiding this comment

prisae Sep 16, 2019

Choose a reason for hiding this comment

prisae Sep 16, 2019

Choose a reason for hiding this comment

ruboerner Sep 16, 2019

Choose a reason for hiding this comment

ruboerner commented Sep 16, 2019

prisae commented Sep 16, 2019

ruboerner commented Sep 16, 2019 • edited Loading

prisae commented Sep 17, 2019

prisae commented Sep 18, 2019

ruboerner commented Sep 18, 2019

prisae commented Sep 18, 2019

ruboerner commented Sep 18, 2019

prisae commented Oct 4, 2019

prisae commented Oct 10, 2019

coveralls commented Oct 10, 2019

prisae commented Oct 10, 2019

prisae commented Oct 10, 2019

prisae commented Oct 11, 2019

prisae commented Oct 11, 2019

prisae commented Oct 14, 2019

prisae commented Oct 15, 2019

prisae commented Oct 19, 2019

ruboerner commented Oct 22, 2019

prisae commented Oct 22, 2019

prisae commented Nov 6, 2019

prisae commented Nov 22, 2019 • edited Loading

ruboerner commented Sep 16, 2019 •

edited

Loading

prisae commented Nov 22, 2019 •

edited

Loading