TxMagneticPoint; generalized adj. RxMagneticPoint

CHANGELOG.rst

@@ -9,6 +9,14 @@ Changelog
 latest
 ------
 
+- ``electrodes``:
+
+  - New source ``TxMagneticPoint`` (requires ``discretize``; mainly used as
+    adjoint source for magnetic receivers; does not work in the presence of
+    magnetic permeabilities in the vincinity of the source).
+  - Both receivers (``Rx{Electric;Magnetic}Point``) can now produce their
+    proper adjoint (thanks to @sgkang!).
+
 - Changes in Simulation and parallel execution.
 
   - Parallel computation is not sharing the simulation any longer.
@@ -21,7 +29,6 @@ latest
   - Models and Fields return itself (not a copy) when the grid provided to
     ``interpolate_to_grid`` is the same as the current one.
 
-- WIP: Proper adjoint for magnetic receivers.
 
 
 v1.2.1: Remove optimize & bug fix

emg3d/__init__.py

@@ -16,8 +16,10 @@
 
 # Import most important functions and classes
 from emg3d.electrodes import (
-        TxElectricPoint, TxElectricDipole, TxMagneticDipole, TxElectricWire,
-        RxElectricPoint, RxMagneticPoint,
+    TxElectricPoint, TxMagneticPoint,
+    TxElectricDipole, TxMagneticDipole,
+    TxElectricWire,
+    RxElectricPoint, RxMagneticPoint,
 )
 from emg3d.fields import Field, get_source_field, get_magnetic_field
 from emg3d.io import save, load

emg3d/electrodes.py

@@ -24,10 +24,14 @@
 
 from emg3d import fields, utils
 
-__all__ = ['Wire', 'Point', 'Dipole', 'Source', 'TxElectricDipole',
-           'TxMagneticDipole', 'TxElectricWire', 'RxElectricPoint',
-           'RxMagneticPoint', 'rotation', 'point_to_dipole', 'dipole_to_point',
-           'point_to_square_loop']
+__all__ = [
+    'Wire', 'Point', 'Dipole', 'Source',
+    'TxElectricPoint', 'TxMagneticPoint',
+    'TxElectricDipole', 'TxMagneticDipole',
+    'TxElectricWire',
+    'RxElectricPoint', 'RxMagneticPoint',
+    'rotation', 'point_to_dipole', 'dipole_to_point', 'point_to_square_loop',
+]
 
 
 # BASE ELECTRODE TYPES
@@ -444,7 +448,33 @@ class TxElectricPoint(Source, Point):
     """
 
     def __init__(self, coordinates, strength=1.0):
-        """Initiate an electric dipole source."""
+        """Initiate an electric point source."""
+
+        super().__init__(coordinates=coordinates, strength=strength)
+
+
+@utils._known_class
+class TxMagneticPoint(Source, Point):
+    """Magnetic point source.
+
+    .. note::
+
+        The magnetic point source is not implemented for magnetic permeability.
+
+
+    Parameters
+    ----------
+    coordinates : array_like
+        Point coordinates in the format (x, y, z, azimuth, elevation).
+
+    strength : float, default: 1.0
+        Source strength (A). (Note that the source field is always an electric
+        field; in this case the electric field due to a magnetic point.)
+
+    """
+
+    def __init__(self, coordinates, strength=1.0):
+        """Initiate an magnetic point source."""
 
         super().__init__(coordinates=coordinates, strength=strength)
 
@@ -607,25 +637,13 @@ class RxElectricPoint(Receiver, Point):
         sources, such as is the case in a :class:`emg3d.surveys.Survey`.
 
     """
+    _adjoint_source = TxElectricPoint
 
     def __init__(self, coordinates, relative=False):
         """Initiate an electric point receiver."""
 
         super().__init__(coordinates=coordinates, relative=relative)
 
-    def adjoint_source(self, source, frequency, strength, grid):
-        """TODO document & test."""
-
-        # Get absolute coordinates as fct of source.
-        # (Only relevant in case of "relative" receivers.)
-        coords = self.coordinates_abs(source)
-
-        # Create source.
-        src = TxElectricPoint(coords, strength=strength)
-
-        # Return adjoint source.
-        return src.get_field(grid=grid, frequency=frequency)
-
 
 @utils._known_class
 class RxMagneticPoint(Receiver, Point):
@@ -645,44 +663,13 @@ class RxMagneticPoint(Receiver, Point):
         sources, such as is the case in a :class:`emg3d.surveys.Survey`.
 
     """
+    _adjoint_source = TxMagneticPoint
 
     def __init__(self, coordinates, relative=False):
         """Initiate a magnetic point receiver."""
 
         super().__init__(coordinates=coordinates, relative=relative)
 
-    @utils._requires('discretize')
-    def adjoint_source(self, source, frequency, strength, grid):
-        """TODO document & test; generalize."""
-
-        # Get absolute coordinates as fct of source.
-        # (Only relevant in case of "relative" receivers.)
-        coords = np.array(self.coordinates_abs(source))
-
-        # TODO Requires a generalization, but should be simple by
-        # combining x, y, z.
-        if (self.azimuth == 0) & (self.elevation == 0):
-            location_type = 'Fx'
-        elif (self.azimuth == 90) & (self.elevation == 0):
-            location_type = 'Fy'
-        elif (self.azimuth == 0) & (self.elevation == 90):
-            location_type = 'Fz'
-        else:
-            raise NotImplementedError
-
-        # Use of discretize for calculating the adjoint source.
-        # TODO implement so it is possible also without discretize.
-        C = grid.edge_curl
-        P = grid.get_interpolation_matrix(
-                coords[:3], location_type=location_type)
-
-        # h = C*e / (i*omega*mu)
-        # smu0 = i*omega*mu
-        h_deriv = -(C.T @ P.T).toarray().ravel() * strength
-
-        # Return adjoint source.
-        return fields.Field(grid=grid, data=h_deriv, frequency=frequency)
-
 
 # ROTATIONS AND CONVERSIONS
 def point_to_dipole(point, length, deg=True):

emg3d/fields.py

@@ -497,6 +497,8 @@ def get_source_field(grid, source, frequency, **kwargs):
     # Get vector field
     if isinstance(source, electrodes.TxElectricPoint):
         vfield = _point_vector(grid, source.coordinates)
+    elif isinstance(source, electrodes.TxMagneticPoint):
+        vfield = _point_vector_magnetic(grid, source.coordinates, frequency)
     else:
         vfield = _dipole_vector(grid, source.points)
 
@@ -739,6 +741,50 @@ def get_index_and_strength(ic, nc, csrc, cc):
     return vfield
 
 
+@utils._requires('discretize')
+def _point_vector_magnetic(grid, coordinates, frequency):
+    """Get magnetic point source using discretize functionality.
+
+
+    Parameters
+    ----------
+    grid : TensorMesh
+        The grid; a :class:`emg3d.meshes.TensorMesh` instance.
+
+    coordinates : array_like
+        Source coordinates in the format (x, y, z, azimuth, elevation).
+
+    frequency : float
+        Source frequency (Hz).
+
+
+    Returns
+    -------
+    vfield : Field
+        Source field, a :class:`emg3d.fields.Field` instance.
+
+    """
+
+    # Use discretize to get the source field for a magnetic point.
+    coords = np.array(coordinates)
+    rot = electrodes.rotation(coords[3], coords[4])
+    interp = (
+        rot[0]*grid.get_interpolation_matrix(coords[:3], 'faces_x') +
+        rot[1]*grid.get_interpolation_matrix(coords[:3], 'faces_y') +
+        rot[2]*grid.get_interpolation_matrix(coords[:3], 'faces_z')
+    )
+
+    # Compute the magnetic vector.
+    vfield = Field(grid, frequency=frequency)
+    vfield.field = -(grid.edge_curl.T @ interp.T).toarray().ravel()
+
+    # Divide by s*mu_0.
+    if frequency is not None:
+        vfield.field /= -vfield.smu0
+
+    return vfield
+
+
 def _dipole_vector(grid, points, decimals=9):
     """Get n-segment dipole source by distributing them to the relevant cells.
 

emg3d/simulations.py

@@ -763,9 +763,6 @@ def gradient(self):
             interpolation, so the residual should come from a linear
             interpolation.
 
-            Also, the adjoint test for magnetic receivers does not yet pass.
-            Electric receivers are good to go.
-
         .. note::
 
             The currently implemented gradient is only for isotropic models
@@ -1018,13 +1015,15 @@ def _get_rfield(self, source, frequency):
             if np.isnan(residual[i]):
                 continue
 
+            # Get absolute coordinates as fct of source.
+            # (Only relevant in case of "relative" receivers.)
+            coords = rec.coordinates_abs(self.survey.sources[source])
+
+            # Create adjoint source.
+            src = rec._adjoint_source(coords, strength=strength[i])
+
             # Get and add this source field.
-            rfield.field += rec.adjoint_source(
-                source=self.survey.sources[source],
-                frequency=freq,
-                strength=strength[i],
-                grid=grid,
-            ).field
+            rfield.field += src.get_field(grid=grid, frequency=freq).field
 
         return rfield
 

tests/test_electrodes.py

@@ -236,6 +236,22 @@ def test_tx_electric_dipole():
     assert "m; e2={1.0" in rep
 
 
+def test_tx_magnetic_point():
+    coo = (0, 0, 0, 45, 45)
+    s1a = electrodes.TxMagneticPoint(coo, strength=np.pi)
+    assert s1a.xtype == 'magnetic'
+    assert s1a._prefix == 'TxMP'
+    s1b = electrodes.TxMagneticPoint.from_dict(s1a.to_dict())
+    assert s1a == s1b
+    assert_allclose(s1b.coordinates, coo)
+    assert_allclose(s1b.points, np.atleast_2d(coo[:3]))
+    assert_allclose(s1b.strength, np.pi)
+
+    rep = s1a.__repr__()
+    assert "3.1 A" in rep
+    assert "=0.0 m, θ=45.0°" in rep
+
+
 def test_tx_magnetic_dipole():
     s4a = electrodes.TxMagneticDipole(
             (0, 0, 0, 45, 45), strength=np.pi, length=4)
@@ -300,6 +316,9 @@ def test_receiver():
 
 def test_rx_electric_point():
     r1a = electrodes.RxElectricPoint((1200, -56, 23.214, 368, 15), True)
+
+    assert r1a._adjoint_source == electrodes.TxElectricPoint
+
     assert r1a.xtype == 'electric'
     assert r1a._prefix == 'RxEP'
     assert r1a.relative is True
@@ -325,6 +344,9 @@ def test_rx_electric_point():
 
 def test_rx_magnetic_point():
     r1a = electrodes.RxMagneticPoint((-1200, 56, -23.214, 0, 90))
+
+    assert r1a._adjoint_source == electrodes.TxMagneticPoint
+
     assert r1a.xtype == 'magnetic'
     assert r1a._prefix == 'RxMP'
     assert r1a.relative is False