Repository: simpeg/simpeg
Pull request #1288 (2023-09-20 23:18:49)

Triaxial magnetic gradient forward modelling
<!-- 
Thanks for contributing a pull request to SimPEG!
Remember to use a personal fork of SimPEG to propose changes.

Check out the stages of a pull request at
https://docs.simpeg.xyz/content/getting_started/contributing/pull-requests.html

Note that we are a team of volunteers and we appreciate your
patience during the review process.

Again, thanks for contributing!

Feel free to remove lines from this template that do not apply to you pull request.
-->

#### Summary
<!-- Add a summary of this Pull Request -->
This PR adds magnetic forward modelling for the gradient of the TMI, measured when several magnetometers are mounted on the same aircraft. A good example is [the Heli-GT system from SHA geophysics ](https://www.shageophysics.com/#js__scroll-to-section_heli-gt) but many airborne systems have several magnetometers mounted.

#### PR Checklist
* [x] If this is a work in progress PR, set as a Draft PR
* [x] Linted my code according to the [style guides](https://docs.simpeg.xyz/content/getting_started/contributing/code-style.html).
* [x] Added [tests](https://docs.simpeg.xyz/content/getting_started/practices.html#testing) to verify changes to the code.
* [x] Added necessary documentation to any new functions/classes following the
      expect [style](https://docs.simpeg.xyz/content/getting_started/practices.html#documentation).
* [x] Marked as ready for review (if this is was a draft PR), and converted 
      to a Pull Request
* [x] Tagged ``@simpeg/simpeg-developers`` when ready for review.

#### Reference issue
<!--Example: write "Closes #NNNN" to automatically close that issue on merge.-->

#### What does this implement/fix?
<!--Please explain your changes.-->

#### Additional information
<!--Any additional information you think is important.-->


<!--
Once all tests pass and the code has been reviewed and approved, it will be merged into main
-->
<img width="1275" alt="Screenshot 2023-09-19 at 7 38 12 PM" src="https://github.com/simpeg/simpeg/assets/97514898/22282b87-2f94-4996-a698-a08ce239d5b6">



thibaut-kobold: @domfournier @jcapriot : what would be the best way to test it? Would we need to update `geoana` prior to this PR or is there a path to keep this PR self-sufficient ?
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/simpeg/simpeg/pull/1288?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg) Report
> Merging [#1288](https://app.codecov.io/gh/simpeg/simpeg/pull/1288?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg) (80c8aff) into [main](https://app.codecov.io/gh/simpeg/simpeg/commit/b42d0d02bba0ef5a39c0ad011073cfd7963d0b7c?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg) (b42d0d0) will **decrease** coverage by `61.42%`.
> The diff coverage is `48.27%`.

```diff
@@             Coverage Diff             @@
##             main    #1288       +/-   ##
===========================================
- Coverage   82.36%   20.95%   -61.42%     
===========================================
  Files         164      164               
  Lines       25098    25111       +13     
===========================================
- Hits        20673     5263    -15410     
- Misses       4425    19848    +15423     
```


| [Files Changed](https://app.codecov.io/gh/simpeg/simpeg/pull/1288?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg) | Coverage Δ | |
|---|---|---|
| [SimPEG/potential\_fields/magnetics/receivers.py](https://app.codecov.io/gh/simpeg/simpeg/pull/1288?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg#diff-U2ltUEVHL3BvdGVudGlhbF9maWVsZHMvbWFnbmV0aWNzL3JlY2VpdmVycy5weQ==) | `100.00% <ø> (ø)` | |
| [SimPEG/potential\_fields/magnetics/simulation.py](https://app.codecov.io/gh/simpeg/simpeg/pull/1288?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg#diff-U2ltUEVHL3BvdGVudGlhbF9maWVsZHMvbWFnbmV0aWNzL3NpbXVsYXRpb24ucHk=) | `50.23% <48.27%> (-21.51%)` | :arrow_down: |

... and [146 files with indirect coverage changes](https://app.codecov.io/gh/simpeg/simpeg/pull/1288/indirect-changes?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg)


thibaut-kobold: @simpeg/simpeg-developers

Patch:
diff --git a/SimPEG/potential_fields/magnetics/receivers.py b/SimPEG/potential_fields/magnetics/receivers.py
index 8a290959c..6f885ad63 100644
--- a/SimPEG/potential_fields/magnetics/receivers.py
+++ b/SimPEG/potential_fields/magnetics/receivers.py
@@ -23,6 +23,9 @@ class Point(survey.BaseRx):
         - "byy"  --> y-derivative of the y-component
         - "byz"  --> z-derivative of the y-component (and visa versa)
         - "bzz"  --> z-derivative of the z-component
+        - "tmi_x"--> x-derivative of the total magnetic intensity data
+        - "tmi_y"--> y-derivative of the total magnetic intensity data
+        - "tmi_z"--> z-derivative of the total magnetic intensity data
 
     Notes
     -----
@@ -51,6 +54,9 @@ class Point(survey.BaseRx):
                     "by",
                     "bz",
                     "tmi",
+                    "tmi_x",
+                    "tmi_y",
+                    "tmi_z",
                 ],
             )
         self.components = components
diff --git a/SimPEG/potential_fields/magnetics/simulation.py b/SimPEG/potential_fields/magnetics/simulation.py
index 8d926d5e4..11fdfb5a7 100644
--- a/SimPEG/potential_fields/magnetics/simulation.py
+++ b/SimPEG/potential_fields/magnetics/simulation.py
@@ -1,27 +1,24 @@
 import numpy as np
 import scipy.sparse as sp
-from scipy.constants import mu_0
-
-from SimPEG import utils
-from ..base import BasePFSimulation, BaseEquivalentSourceLayerSimulation
-from ...base import BaseMagneticPDESimulation
-from .survey import Survey
-from .analytics import CongruousMagBC
-
-from SimPEG import Solver
-from SimPEG import props
-
-from SimPEG.utils import mkvc, mat_utils, sdiag
-from SimPEG.utils.code_utils import validate_string, deprecate_property, validate_type
 from geoana.kernels import (
-    prism_fzz,
-    prism_fzx,
-    prism_fzy,
-    prism_fzzz,
     prism_fxxy,
     prism_fxxz,
     prism_fxyz,
+    prism_fzx,
+    prism_fzy,
+    prism_fzz,
+    prism_fzzz,
 )
+from scipy.constants import mu_0
+
+from SimPEG import Solver, props, utils
+from SimPEG.utils import mat_utils, mkvc, sdiag
+from SimPEG.utils.code_utils import deprecate_property, validate_string, validate_type
+
+from ...base import BaseMagneticPDESimulation
+from ..base import BaseEquivalentSourceLayerSimulation, BasePFSimulation
+from .analytics import CongruousMagBC
+from .survey import Survey
 
 
 class Simulation3DIntegral(BasePFSimulation):
@@ -246,7 +243,7 @@ class Simulation3DIntegral(BasePFSimulation):
 
         components: list[str]
             List of magnetic components chosen from:
-            'bx', 'by', 'bz', 'bxx', 'bxy', 'bxz', 'byy', 'byz', 'bzz'
+            'tmi', 'bx', 'by', 'bz', 'bxx', 'bxy', 'bxz', 'byy', 'byz', 'bzz', 'tmi_x', 'tmi_y', 'tmi_z'
 
         OUTPUT:
         Tx = [Txx Txy Txz]
@@ -283,33 +280,33 @@ class Simulation3DIntegral(BasePFSimulation):
             #     # inside an active cell.
             #     node_evals["gzz"] = -node_evals["gxx"] - node_evals["gyy"]
 
-        if "bxx" in components:
+        if "bxx" in components or "tmi_x" in components:
             node_evals["gxxx"] = prism_fzzz(dy, dz, dx)
             node_evals["gxxy"] = prism_fxxy(dx, dy, dz)
             node_evals["gxxz"] = prism_fxxz(dx, dy, dz)
-        if "bxy" in components:
+        if "bxy" in components or "tmi_x" in components or "tmi_y" in components:
             if "gxxy" not in node_evals:
                 node_evals["gxxy"] = prism_fxxy(dx, dy, dz)
             node_evals["gyyx"] = prism_fxxz(dy, dz, dx)
             node_evals["gxyz"] = prism_fxyz(dx, dy, dz)
-        if "bxz" in components:
+        if "bxz" in components or "tmi_x" in components or "tmi_z" in components:
             if "gxxz" not in node_evals:
                 node_evals["gxxz"] = prism_fxxz(dx, dy, dz)
             if "gxyz" not in node_evals:
                 node_evals["gxyz"] = prism_fxyz(dx, dy, dz)
             node_evals["gzzx"] = prism_fxxy(dz, dx, dy)
-        if "byy" in components:
+        if "byy" in components or "tmi_y" in components:
             if "gyyx" not in node_evals:
                 node_evals["gyyx"] = prism_fxxz(dy, dz, dx)
             node_evals["gyyy"] = prism_fzzz(dz, dx, dy)
             node_evals["gyyz"] = prism_fxxy(dy, dz, dx)
-        if "byz" in components:
+        if "byz" in components or "tmi_y" in components or "tmi_z" in components:
             if "gxyz" not in node_evals:
                 node_evals["gxyz"] = prism_fxyz(dx, dy, dz)
             if "gyyz" not in node_evals:
                 node_evals["gyyz"] = prism_fxxy(dy, dz, dx)
             node_evals["gzzy"] = prism_fxxz(dz, dx, dy)
-        if "bzz" in components:
+        if "bzz" in components or "tmi_z" in components:
             if "gzzx" not in node_evals:
                 node_evals["gzzx"] = prism_fxxy(dz, dx, dy)
             if "gzzy" not in node_evals:
@@ -355,6 +352,57 @@ class Simulation3DIntegral(BasePFSimulation):
                     + tmi[1] * node_evals["gyz"]
                     + tmi[2] * node_evals["gzz"]
                 )
+            elif component == "tmi_x":
+                tmi = self.tmi_projection
+                vals_x = (
+                    tmi[0] * node_evals["gxxx"]
+                    + tmi[1] * node_evals["gxxy"]
+                    + tmi[2] * node_evals["gxxz"]
+                )
+                vals_y = (
+                    tmi[0] * node_evals["gxxy"]
+                    + tmi[1] * node_evals["gyyx"]
+                    + tmi[2] * node_evals["gxyz"]
+                )
+                vals_z = (
+                    tmi[0] * node_evals["gxxz"]
+                    + tmi[1] * node_evals["gxyz"]
+                    + tmi[2] * node_evals["gzzx"]
+                )
+            elif component == "tmi_y":
+                tmi = self.tmi_projection
+                vals_x = (
+                    tmi[0] * node_evals["gxxy"]
+                    + tmi[1] * node_evals["gyyx"]
+                    + tmi[2] * node_evals["gxyz"]
+                )
+                vals_y = (
+                    tmi[0] * node_evals["gyyx"]
+                    + tmi[1] * node_evals["gyyy"]
+                    + tmi[2] * node_evals["gyyz"]
+                )
+                vals_z = (
+                    tmi[0] * node_evals["gxyz"]
+                    + tmi[1] * node_evals["gyyz"]
+                    + tmi[2] * node_evals["gzzy"]
+                )
+            elif component == "tmi_z":
+                tmi = self.tmi_projection
+                vals_x = (
+                    tmi[0] * node_evals["gxxz"]
+                    + tmi[1] * node_evals["gxyz"]
+                    + tmi[2] * node_evals["gzzx"]
+                )
+                vals_y = (
+                    tmi[0] * node_evals["gxyz"]
+                    + tmi[1] * node_evals["gyyz"]
+                    + tmi[2] * node_evals["gzzy"]
+                )
+                vals_z = (
+                    tmi[0] * node_evals["gzzx"]
+                    + tmi[1] * node_evals["gzzy"]
+                    + tmi[2] * node_evals["gzzz"]
+                )
             elif component == "bxx":
                 vals_x = node_evals["gxxx"]
                 vals_y = node_evals["gxxy"]
@@ -949,11 +997,11 @@ def MagneticsDiffSecondaryInv(mesh, model, data, **kwargs):
 
     """
     from SimPEG import (
-        optimization,
-        regularization,
         directives,
-        objective_function,
         inversion,
+        objective_function,
+        optimization,
+        regularization,
     )
 
     prob = Simulation3DDifferential(mesh, survey=data, mu=model)


Test patch:
diff --git a/tests/pf/test_forward_Mag_Linear.py b/tests/pf/test_forward_Mag_Linear.py
index 3f412651d..662556b93 100644
--- a/tests/pf/test_forward_Mag_Linear.py
+++ b/tests/pf/test_forward_Mag_Linear.py
@@ -1,10 +1,12 @@
 import unittest
+
 import discretize
-from SimPEG import utils, maps
-from SimPEG.potential_fields import magnetics as mag
+import numpy as np
 from geoana.em.static import MagneticPrism
 from scipy.constants import mu_0
-import numpy as np
+
+from SimPEG import maps, utils
+from SimPEG.potential_fields import magnetics as mag
 
 
 def test_ana_mag_forward():
@@ -47,7 +49,7 @@ def test_ana_mag_forward():
     active_cells = model != 0.0
     model_reduced = model[active_cells]
 
-    # Create reduced identity map for Linear Pproblem
+    # Create reduced identity map for Linear Problem
     idenMap = maps.IdentityMap(nP=int(sum(active_cells)))
 
     # Create plane of observations
@@ -59,10 +61,12 @@ def test_ana_mag_forward():
     components = ["bx", "by", "bz", "tmi"]
 
     rxLoc = mag.Point(locXyz, components=components)
-    srcField = mag.SourceField([rxLoc], parameters=H0)
+    srcField = mag.UniformBackgroundField(
+        [rxLoc], amplitude=H0[0], inclination=H0[1], declination=H0[2]
+    )
     survey = mag.Survey(srcField)
 
-    # Creat reduced identity map for Linear Pproblem
+    # Creat reduced identity map for Linear Problem
     idenMap = maps.IdentityMap(nP=int(sum(active_cells)))
 
     sim = mag.Simulation3DIntegral(
@@ -101,6 +105,116 @@ def test_ana_mag_forward():
     np.testing.assert_allclose(d_t, d @ tmi)
 
 
+def test_ana_mag_tmi_grad_forward():
+    nx = 61
+    ny = 61
+
+    H0 = (50000.0, 60.0, 250.0)
+    b0 = mag.analytics.IDTtoxyz(-H0[1], H0[2], H0[0])
+    chi1 = 0.01
+    chi2 = 0.02
+
+    # Define a mesh
+    cs = 0.2
+    hxind = [(cs, 41)]
+    hyind = [(cs, 41)]
+    hzind = [(cs, 41)]
+    mesh = discretize.TensorMesh([hxind, hyind, hzind], "CCC")
+
+    # create a model of two blocks, 1 inside the other
+    block1 = np.array([[-1.5, 1.5], [-1.5, 1.5], [-1.5, 1.5]])
+    block2 = np.array([[-0.7, 0.7], [-0.7, 0.7], [-0.7, 0.7]])
+
+    def get_block_inds(grid, block):
+        return np.where(
+            (grid[:, 0] > block[0, 0])
+            & (grid[:, 0] < block[0, 1])
+            & (grid[:, 1] > block[1, 0])
+            & (grid[:, 1] < block[1, 1])
+            & (grid[:, 2] > block[2, 0])
+            & (grid[:, 2] < block[2, 1])
+        )
+
+    block1_inds = get_block_inds(mesh.cell_centers, block1)
+    block2_inds = get_block_inds(mesh.cell_centers, block2)
+
+    model = np.zeros(mesh.n_cells)
+    model[block1_inds] = chi1
+    model[block2_inds] = chi2
+
+    active_cells = model != 0.0
+    model_reduced = model[active_cells]
+
+    # Create reduced identity map for Linear Problem
+    idenMap = maps.IdentityMap(nP=int(sum(active_cells)))
+
+    # Create plane of observations
+    xr = np.linspace(-20, 20, nx)
+    dxr = xr[1] - xr[0]
+    yr = np.linspace(-20, 20, ny)
+    dyr = yr[1] - yr[0]
+    X, Y = np.meshgrid(xr, yr)
+    Z = np.ones_like(X) * 3.0
+    locXyz = np.c_[X.reshape(-1), Y.reshape(-1), Z.reshape(-1)]
+    components = ["tmi", "tmi_x", "tmi_y", "tmi_z"]
+
+    rxLoc = mag.Point(locXyz, components=components)
+    srcField = mag.UniformBackgroundField(
+        [rxLoc], amplitude=H0[0], inclination=H0[1], declination=H0[2]
+    )
+    survey = mag.Survey(srcField)
+
+    # Creat reduced identity map for Linear Problem
+    idenMap = maps.IdentityMap(nP=int(sum(active_cells)))
+
+    sim = mag.Simulation3DIntegral(
+        mesh,
+        survey=survey,
+        chiMap=idenMap,
+        ind_active=active_cells,
+        store_sensitivities="forward_only",
+        n_processes=None,
+    )
+
+    data = sim.dpred(model_reduced)
+    tmi = data[0::4]
+    d_x = data[1::4]
+    d_y = data[2::4]
+    d_z = data[3::4]
+
+    # Compute analytical response from magnetic prism
+    prism_1 = MagneticPrism(block1[:, 0], block1[:, 1], chi1 * b0 / mu_0)
+    prism_2 = MagneticPrism(block2[:, 0], block2[:, 1], -chi1 * b0 / mu_0)
+    prism_3 = MagneticPrism(block2[:, 0], block2[:, 1], chi2 * b0 / mu_0)
+
+    d = (
+        prism_1.magnetic_field_gradient(locXyz)
+        + prism_2.magnetic_field_gradient(locXyz)
+        + prism_3.magnetic_field_gradient(locXyz)
+    ) * mu_0
+    tmi_x = (d[:, 0, 0] * b0[0] + d[:, 0, 1] * b0[1] + d[:, 0, 2] * b0[2]) / H0[0]
+    tmi_y = (d[:, 1, 0] * b0[0] + d[:, 1, 1] * b0[1] + d[:, 1, 2] * b0[2]) / H0[0]
+    tmi_z = (d[:, 2, 0] * b0[0] + d[:, 2, 1] * b0[1] + d[:, 2, 2] * b0[2]) / H0[0]
+    np.testing.assert_allclose(d_x, tmi_x, rtol=1e-10, atol=1e-12)
+    np.testing.assert_allclose(d_y, tmi_y, rtol=1e-10, atol=1e-12)
+    np.testing.assert_allclose(d_z, tmi_z, rtol=1e-10, atol=1e-12)
+
+    # finite difference test y-grad
+    np.testing.assert_allclose(
+        np.diff(tmi.reshape(nx, ny, order="F")[:, ::2], axis=1) / (2 * dyr),
+        tmi_y.reshape(nx, ny, order="F")[:, 1::2],
+        atol=1.0,
+        rtol=1e-1,
+    )
+    # finite difference test x-grad
+    np.testing.assert_allclose(
+        np.diff(tmi.reshape(nx, ny, order="F")[::2, :], axis=0) / (2 * dxr),
+        tmi_x.reshape(nx, ny, order="F")[1::2, :],
+        atol=1.0,
+        rtol=1e-1,
+    )
+
+
 def test_ana_mag_grad_forward():
     nx = 5
     ny = 5
@@ -141,7 +255,7 @@ def test_ana_mag_grad_forward():
     active_cells = model != 0.0
     model_reduced = model[active_cells]
 
-    # Create reduced identity map for Linear Pproblem
+    # Create reduced identity map for Linear Problem
     idenMap = maps.IdentityMap(nP=int(sum(active_cells)))
 
     # Create plane of observations
@@ -153,10 +267,12 @@ def test_ana_mag_grad_forward():
     components = ["bxx", "bxy", "bxz", "byy", "byz", "bzz"]
 
     rxLoc = mag.Point(locXyz, components=components)
-    srcField = mag.SourceField([rxLoc], parameters=H0)
+    srcField = mag.UniformBackgroundField(
+        [rxLoc], amplitude=H0[0], inclination=H0[1], declination=H0[2]
+    )
     survey = mag.Survey(srcField)
 
-    # Creat reduced identity map for Linear Pproblem
+    # Creat reduced identity map for Linear Problem
     idenMap = maps.IdentityMap(nP=int(sum(active_cells)))
 
     sim = mag.Simulation3DIntegral(
@@ -245,10 +361,12 @@ def test_ana_mag_vec_forward():
     components = ["bx", "by", "bz", "tmi"]
 
     rxLoc = mag.Point(locXyz, components=components)
-    srcField = mag.SourceField([rxLoc], parameters=H0)
+    srcField = mag.UniformBackgroundField(
+        [rxLoc], amplitude=H0[0], inclination=H0[1], declination=H0[2]
+    )
     survey = mag.Survey(srcField)
 
-    # Create reduced identity map for Linear Pproblem
+    # Create reduced identity map for Linear Problem
     idenMap = maps.IdentityMap(nP=int(sum(active_cells)) * 3)
 
     sim = mag.Simulation3DIntegral(
@@ -331,10 +449,12 @@ def test_ana_mag_amp_forward():
     components = ["bx", "by", "bz"]
 
     rxLoc = mag.Point(locXyz, components=components)
-    srcField = mag.SourceField([rxLoc], parameters=H0)
+    srcField = mag.UniformBackgroundField(
+        [rxLoc], amplitude=H0[0], inclination=H0[1], declination=H0[2]
+    )
     survey = mag.Survey(srcField)
 
-    # Create reduced identity map for Linear Pproblem
+    # Create reduced identity map for Linear Problem
     idenMap = maps.IdentityMap(nP=int(sum(active_cells)) * 3)
 
     sim = mag.Simulation3DIntegral(