From a637736f068240aa02a5488ba33786e76f6df8de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Fri, 25 Oct 2024 11:01:36 +0200 Subject: [PATCH] Add possibility for el. current density in bipole (#220) --- CHANGELOG.rst | 8 ++++++++ empymod/model.py | 31 +++++++++++++++++++++++++++---- tests/test_model.py | 22 ++++++++++++++++++++++ 3 files changed, 57 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 9ca4145..b670ce8 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -10,6 +10,14 @@ v2.3.x """""" +v2.4.0 El. current density +-------------------------- + +**2024-10-25** + +The ``mrec`` keyword in ``empymod.bipole`` can now be set to ``j``, in which case the electric current density (A/m2) is returned. + + v2.3.3 pyproject.toml --------------------- diff --git a/empymod/model.py b/empymod/model.py index 551f52b..48ff55a 100644 --- a/empymod/model.py +++ b/empymod/model.py @@ -151,8 +151,13 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, #mpermH = #mpermV = #res. If mpermH is provided but not mpermV, isotropic behaviour is assumed. - msrc, mrec : bool, default: False + msrc, mrec : bool or string, default: False If True, source/receiver (msrc/mrec) is magnetic, else electric. + The receiver can also be set to `mrec='j'`, in which case the electric + current density is returned (with the approximation that the current + density is proportional to the electric field, J=σE). Only implemented + for isotropic resistivities and electric permittivities at receiver + level. srcpts, recpts : int, default: 1 Number of integration points for bipole source/receiver: @@ -209,7 +214,7 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, - If > 0, interpolation is used. - `diff_quad`: criteria when to swap to QUAD (only relevant if - pts_per_dec=-1) (default: 100) + pts_per_dec>0) (default: 100) - `a`: lower limit for QUAD (default: first interval from QWE) - `b`: upper limit for QUAD (default: last interval from QWE) - `limit`: limit for quad (default: maxint) @@ -310,6 +315,7 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, - If rec is electric, returns E [V/m]. - If rec is magnetic, returns H [A/m]. + - If rec is `j`, returns J [A/m2]. EMArray is a subclassed ndarray with `.pha` and `.amp` attributes (only relevant for frequency-domain data). @@ -422,6 +428,13 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, src, nsrc, nsrcz, srcdipole = check_bipole(src, 'src') rec, nrec, nrecz, recdipole = check_bipole(rec, 'rec') + # Check if receiver is a `j`. + if mrec == 'j': + rec_j = True + mrec = False + else: + rec_j = False + # === 3. EM-FIELD CALCULATION ============ # Pre-allocate output EM array @@ -487,6 +500,11 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, # Get layer number in which rec resides lrec, zrec = get_layer_nr(tirec, depth) + # Check eta at receiver level (only isotropic implemented). + if rec_j and verb > 0 and etaH[0, lrec] != etaV[0, lrec]: + print("* WARNING :: `etaH != etaV` at receiver level, " + "only `etaH` considered for e-current density.") + # Gather variables finp = (off, angle, zsrc, zrec, lsrc, lrec, depth, freq, etaH, etaV, zetaH, zetaV, xdirect, isfullspace, ht, @@ -528,6 +546,10 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, src_rec_w *= np.tile(rec_w, isrc) sEM *= src_rec_w + # Multiply with eta of the rec-layer if ecurrent. + if rec_j: + sEM *= etaH[:, lrec, None] + # Add this src-rec signal if nrec == nrecz: if nsrc == nsrcz: # Case 1: Looped over each src and each rec @@ -898,7 +920,8 @@ def loop(src, rec, depth, res, freqtime, signal=None, aniso=None, epermH=None, - True: Magnetic dipole receiver; - False: Electric dipole receiver; - - 'loop': Magnetic receiver consisting of an electric-wire loop. + - 'loop': Magnetic receiver consisting of an electric-wire loop. Only + implemented for isotropic magnetic permeability at loop levels. recpts : int, default: 1 Number of integration points for bipole receiver: @@ -1119,7 +1142,7 @@ def loop(src, rec, depth, res, freqtime, signal=None, aniso=None, epermH=None, # Get layer number in which rec resides lrec, zrec = get_layer_nr(tirec, depth) - # Check mu at receiver level. + # Check mu at receiver level (only isotropic implemented). if rec_loop and verb > 0 and mpermH[lrec] != mpermV[lrec]: print("* WARNING :: `mpermH != mpermV` at receiver level, " "only `mpermH` considered for loop factor.") diff --git a/tests/test_model.py b/tests/test_model.py index 99c503d..3a077d6 100644 --- a/tests/test_model.py +++ b/tests/test_model.py @@ -588,6 +588,28 @@ def test_shape(self): assert_allclose(out[:, 0], b) assert_allclose(out[:, 1], d) + def test_j(self, capsys): + # Compare a electric * sigma with j + inp = { + 'src': [-25, 25, -25, 25, 100, 170.7107], + 'rec': [8000, 200, 300, 0, 0], + 'depth': [0, 250], + 'res': [1e20, 0.3, 5], + 'freqtime': 1, + } + efield = bipole(**inp) + ecurr = bipole(mrec='j', **inp) + assert_allclose(efield/5, ecurr) + + # Check warning + inp['aniso'] = [1, 1, 2] + efield = bipole(**inp) + _, _ = capsys.readouterr() # Empty it + ecurr = bipole(mrec='j', **inp) + out, _ = capsys.readouterr() + assert_allclose(efield/5, ecurr) + assert "* WARNING :: `etaH != etaV` at receiver level, " in out + def test_dipole(): # As this is a subset of bipole, just run two tests to ensure