Repository: simpeg/simpeg

Issue:
BUG: fdem.Simulation1DLayered fails for more than one receiver
### Describe the issue:

I followed the tutorial https://simpeg.xyz/user-tutorials/fwd-fdem-1d, and tried to change the receiver from a single receiver to many receivers, as indicated in the comment for the receiver locations (`or (N, 3) numpy.ndarray`). However, making `N` anything bigger than one causes SimPEG to fail.

### Reproducible code example:

```python
# Run this example: https://simpeg.xyz/user-tutorials/fwd-fdem-1d
# changing the following line
receiver_locations = np.array([10.0, 0.0, 30.0])  # or (N, 3) numpy.ndarray
# to a multi-receiver line, e.g.
receiver_locations = np.array([[10.0, 0.0, 30.0], [11.0, 0.0, 30.0]])
```

### Error message:

```shell
# When dpred_conductivity = simulation_conductivity.dpred(conductivity_model) is called, the error is:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[11], line 1
----> 1 dpred_conductivity = simulation_conductivity.dpred(conductivity_model)
      2 dpred_log_resistivity = simulation_log_resistivity.dpred(log_resistivity_model)
      3 dpred_parametric = simulation_parametric.dpred(parametric_model)

File ~/miniforge3/lib/python3.12/site-packages/simpeg/electromagnetics/frequency_domain/simulation_1d.py:97, in Simulation1DLayered.dpred(self, m, f)
     91 """
     92 Return predicted data.
     93 Predicted data, (`_pred`) are computed when
     94 self.fields is called.
     95 """
     96 if f is None:
---> 97     f = self.fields(m)
     99 return f

File ~/miniforge3/lib/python3.12/site-packages/simpeg/electromagnetics/frequency_domain/simulation_1d.py:126, in Simulation1DLayered.fields(self, m)
    124 rTE = rTE[i_freq]
    125 rTE = np.take_along_axis(rTE, inv_lambs, axis=1)
--> 126 v = W @ ((C0s * rTE) @ self._fhtfilt.j0 + (C1s * rTE) @ self._fhtfilt.j1)
    128 return self._project_to_data(v)

File ~/miniforge3/lib/python3.12/site-packages/scipy/sparse/_base.py:669, in _spbase.__matmul__(self, other)
    666 if isscalarlike(other):
    667     raise ValueError("Scalar operands are not allowed, "
    668                      "use '*' instead")
--> 669 return self._matmul_dispatch(other)

File ~/miniforge3/lib/python3.12/site-packages/scipy/sparse/_base.py:598, in _spbase._matmul_dispatch(self, other)
    595 if other.ndim == 1 or other.ndim == 2 and other.shape[1] == 1:
    596     # dense row or column vector
    597     if other.shape != (N,) and other.shape != (N, 1):
--> 598         raise ValueError('dimension mismatch')
    600     result = self._matmul_vector(np.ravel(other))
    602     if isinstance(other, np.matrix):

ValueError: dimension mismatch
```

### Runtime information:

```shell
--------------------------------------------------------------------------------
  Date: Sun Mar 16 12:10:19 2025 CET

                OS : Linux (Ubuntu 24.04)
            CPU(s) : 16
           Machine : x86_64
      Architecture : 64bit
               RAM : 93.6 GiB
       Environment : Jupyter
       File system : ext4

  Python 3.12.8 | packaged by conda-forge | (main, Dec  5 2024, 14:24:40) [GCC
  13.3.0]

            simpeg : 0.23.0
        discretize : 0.11.1
       pymatsolver : 0.3.1
             numpy : 2.0.2
             scipy : 1.14.1
        matplotlib : 3.10.0
            geoana : 0.7.2
            libdlf : 0.3.0
             numba : 0.60.0
            pandas : 2.2.3
           IPython : 8.31.0
        ipywidgets : 8.1.5
   memory_profiler : 0.61.0
--------------------------------------------------------------------------------
```

### Context for the issue:

_No response_

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/simpeg/simpeg/pull/1637?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 86.58%. Comparing base [(`33c4221`)](https://app.codecov.io/gh/simpeg/simpeg/commit/33c4221d8f4d9c2e4e3e51b948ac755c4d3d78c4?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg) to head [(`8788c32`)](https://app.codecov.io/gh/simpeg/simpeg/commit/8788c325b5460b68e68bc91fbff61c485746567d?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg).

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #1637      +/-   ##
==========================================
+ Coverage   86.55%   86.58%   +0.03%     
==========================================
  Files         409      409              
  Lines       53118    53184      +66     
  Branches     5024     5028       +4     
==========================================
+ Hits        45975    46051      +76     
+ Misses       5724     5717       -7     
+ Partials     1419     1416       -3     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/simpeg/simpeg/pull/1637?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg).

<details><summary> :rocket: New features to boost your workflow: </summary>

- :snowflake: [Test Analytics](https://docs.codecov.com/docs/test-analytics): Detect flaky tests, report on failures, and find test suite problems.
</details>
jcapriot: @prisae I think this should fix the multiple locations on a single receiver.

Patch:
diff --git a/simpeg/electromagnetics/base_1d.py b/simpeg/electromagnetics/base_1d.py
index a871d56bf..afba17062 100644
--- a/simpeg/electromagnetics/base_1d.py
+++ b/simpeg/electromagnetics/base_1d.py
@@ -558,8 +558,9 @@ class BaseEM1DSimulation(BaseSimulation):
                 C1s.append(np.exp(-lambd * (z + h)[:, None]) * C1 / offsets[:, None])
                 lambs.append(lambd)
                 n_w_past += n_w
-                Is.append(np.ones(n_w, dtype=int) * i_count)
-                i_count += 1
+                for _ in range(rx.locations.shape[0]):
+                    Is.append(np.ones(n_w, dtype=int) * i_count)
+                    i_count += 1
 
         # Store these on the simulation for faster future executions
         self._lambs = np.vstack(lambs)
diff --git a/simpeg/electromagnetics/frequency_domain/simulation_1d.py b/simpeg/electromagnetics/frequency_domain/simulation_1d.py
index 7d880dd09..f7fe4d867 100644
--- a/simpeg/electromagnetics/frequency_domain/simulation_1d.py
+++ b/simpeg/electromagnetics/frequency_domain/simulation_1d.py
@@ -288,6 +288,10 @@ class Simulation1DLayered(BaseEM1DSimulation):
                     out[i_dat:i_dat_p1] = v_slice.real
                 elif rx.component == "imag":
                     out[i_dat:i_dat_p1] = v_slice.imag
+                else:
+                    raise NotImplementedError(
+                        f"receiver component {rx.component} not implemented."
+                    )
                 i_dat = i_dat_p1
                 i_v = i_v_p1
         return out


Test patch:
diff --git a/tests/em/em1d/test_EM1D_FD_fwd.py b/tests/em/em1d/test_EM1D_FD_fwd.py
index 8a62a1398..498602433 100644
--- a/tests/em/em1d/test_EM1D_FD_fwd.py
+++ b/tests/em/em1d/test_EM1D_FD_fwd.py
@@ -9,6 +9,7 @@ from geoana.em.fdem import (
     vertical_magnetic_field_horizontal_loop as mag_field,
 )
 import empymod
+import pytest
 
 
 class EM1D_FD_test_failures(unittest.TestCase):
@@ -548,5 +549,62 @@ class EM1D_FD_LineCurrentTest(unittest.TestCase):
         self.assertLess(err, 1e-4)
 
 
+@pytest.mark.parametrize(
+    "rx_class",
+    [fdem.receivers.PointMagneticField, fdem.receivers.PointMagneticFieldSecondary],
+)
+@pytest.mark.parametrize("n_locs1", [1, 4])
+@pytest.mark.parametrize("n_locs2", [1, 4])
+@pytest.mark.parametrize("orientation", ["x", "y", "z"])
+@pytest.mark.parametrize("component", ["real", "imag", "both"])
+def test_rx_loc_shapes(rx_class, n_locs1, n_locs2, orientation, component):
+    offsets = np.full(n_locs1, 100.0)
+    rx1_locs = np.pad(offsets[:, None], ((0, 0), (0, 2)), constant_values=0)
+    offsets = np.full(n_locs2, 100.0)
+    rx2_locs = np.pad(offsets[:, None], ((0, 0), (0, 2)), constant_values=0)
+
+    rx_list = [
+        rx_class(rx1_locs, orientation=orientation, component=component),
+        rx_class(rx2_locs, orientation=orientation, component=component),
+    ]
+    n_d = n_locs1 + n_locs2
+    if component == "both":
+        n_d *= 2
+
+    src = fdem.sources.MagDipole(rx_list, frequency=0.1)
+    srv = fdem.Survey(src)
+
+    sim = fdem.Simulation1DLayered(survey=srv, sigma=[1])
+    d = sim.dpred(None)
+
+    # assert the shape is correct
+    assert d.shape == (n_d,)
+
+    # every value should be the same...
+    d1 = d[srv.get_slice(src, rx_list[0])]
+    d2 = d[srv.get_slice(src, rx_list[1])]
+
+    if component == "both":
+        d1 = d1[::2] + 1j * d1[1::2]
+        d2 = d2[::2] + 1j * d2[1::2]
+    d = np.r_[d1, d2]
+    np.testing.assert_allclose(d, d[0], rtol=1e-12)
+
+    sim.sigmaMap = maps.IdentityMap(nP=1)
+    # make sure forming J works
+    J = sim.getJ(np.ones(1))["ds"]
+    assert J.shape == (n_d, 1)
+
+    # and all of its values are the same too:
+    j1 = J[srv.get_slice(src, rx_list[0]), 0]
+    j2 = J[srv.get_slice(src, rx_list[1]), 0]
+
+    if component == "both":
+        j1 = j1[::2] + 1j * j1[1::2]
+        j2 = j2[::2] + 1j * j2[1::2]
+    J = np.r_[j1, j2]
+    np.testing.assert_allclose(J, J[0], rtol=1e-12)
+
+
 if __name__ == "__main__":
     unittest.main()
diff --git a/tests/em/em1d/test_EM1D_TD_general_jac_layers.py b/tests/em/em1d/test_EM1D_TD_general_jac_layers.py
index f6818ccf8..0580d5407 100644
--- a/tests/em/em1d/test_EM1D_TD_general_jac_layers.py
+++ b/tests/em/em1d/test_EM1D_TD_general_jac_layers.py
@@ -3,6 +3,7 @@ from simpeg import maps
 from discretize import tests
 import numpy as np
 import simpeg.electromagnetics.time_domain as tdem
+import pytest
 
 
 class EM1D_TD_general_Jac_layers_ProblemTests(unittest.TestCase):
@@ -305,5 +306,59 @@ class EM1D_TD_LineCurrent_Jac_layers_ProblemTests(unittest.TestCase):
         self.assertTrue(passed)
 
 
+@pytest.mark.parametrize(
+    "rx_class",
+    [
+        tdem.receivers.PointMagneticField,
+        tdem.receivers.PointMagneticFluxDensity,
+        tdem.receivers.PointMagneticFluxTimeDerivative,
+    ],
+)
+@pytest.mark.parametrize("n_locs1", [1, 4])
+@pytest.mark.parametrize("n_locs2", [1, 4])
+@pytest.mark.parametrize("orientation", ["x", "y", "z"])
+@pytest.mark.parametrize(
+    "waveform", [tdem.sources.StepOffWaveform(), tdem.sources.RampOffWaveform(1e-6)]
+)
+@pytest.mark.parametrize("comparison", ["dpred", "J"])
+def test_rx_loc_shapes(rx_class, n_locs1, n_locs2, orientation, waveform, comparison):
+    offsets = np.full(n_locs1, 100.0)
+    rx1_locs = np.pad(offsets[:, None], ((0, 0), (0, 2)), constant_values=0)
+    offsets = np.full(n_locs2, 100.0)
+    rx2_locs = np.pad(offsets[:, None], ((0, 0), (0, 2)), constant_values=0)
+
+    times = [1e-5, 1e-4]
+    rx_list = [
+        rx_class(rx1_locs, times=times, orientation=orientation),
+        rx_class(rx2_locs, times=times, orientation=orientation),
+    ]
+    n_d = (n_locs1 + n_locs2) * len(times)
+
+    src = tdem.sources.MagDipole(rx_list, waveform=waveform)
+    srv = tdem.Survey(src)
+
+    sim = tdem.Simulation1DLayered(survey=srv, sigma=[1])
+    if comparison == "dpred":
+        d = sim.dpred(None)
+    else:
+        sim.sigmaMap = maps.IdentityMap(nP=1)
+        d = sim.getJ(np.ones(1))["ds"][:, 0]
+
+    # assert the shape is correct
+    assert d.shape == (n_d,)
+
+    # every pair of values should be the same...
+    d1 = d[srv.get_slice(src, rx_list[0])]
+    d2 = d[srv.get_slice(src, rx_list[1])]
+
+    # get the data into an n_loc * n_times shape
+    d = np.r_[
+        d1.reshape(2, -1).T,
+        d2.reshape(2, -1).T,
+    ]
+    d_compare = d[0] * np.ones_like(d)
+    np.testing.assert_equal(d, d_compare)
+
+
 if __name__ == "__main__":
     unittest.main()