Error from multiple receivers per source dipole

DC and IP
1

Okay, so..

I’m using the read_dcip_xyz() function for a 3D DCIP inversion and run into issues with the number of sources vs receivers generated.

The actual survey was completed using a Pole-Dipole array, where the B electrode is located at ‘infinity’, i.e., far away from the survey site.

When running dc_simulation.dpred(starting_conductivity_model), I get the following error:


File C:\anaconda\Lib\site-packages\simpeg\simulation.py:195, in BaseSimulation.dpred(self, m, f)
    193     for rx in src.receiver_list:
    194         src_rx_slice = survey_slices[src, rx]
--> 195         dpred[src_rx_slice] = mkvc(rx.eval(src, self.mesh, f))
    196 return mkvc(dpred)

ValueError: could not broadcast input array from shape (36,) into shape (6,)

Looking at that source, I see 6 associated receivers:

In [21]: src
Out[21]: Dipole(a: [4.9074050e+05 6.7348901e+06 7.3039350e+02]; b: [4.9224050e+05 6.7353901e+06 7.3039350e+02])

In [22]: src.receiver_list
Out[22]: 
Dipole(m: [4.9094050e+05 6.7349901e+06 7.3039350e+02]; n: [4.9084050e+05 6.7349901e+06 7.3039350e+02]),
Dipole(m: [4.9064050e+05 6.7349901e+06 7.3039350e+02]; n: [4.9064050e+05 6.7349901e+06 7.3039350e+02]),
Dipole(m: [4.9064050e+05 6.7349901e+06 7.3039350e+02]; n: [4.9054050e+05 6.7349901e+06 7.3039350e+02]),
Dipole(m: [4.9054050e+05 6.7349901e+06 7.3039350e+02]; n: [4.9054050e+05 6.7349901e+06 7.3039350e+02]),
Dipole(m: [4.9054050e+05 6.7349901e+06 7.3039350e+02]; n: [4.9044050e+05 6.7349901e+06 7.3039350e+02]),
Dipole(m: [4.9044050e+05 6.7349901e+06 7.3039350e+02]; n: [4.9044050e+05 6.7349901e+06 7.3039350e+02])

I can get around this issue by manually assigning each receiver to the appropriate source:

def flatten_survey(survey):
    new_sources = []
    for src in survey.source_list:
        for rx in src.receiver_list:
            # Each rx contain multiple MN pairs
            for i in range(rx.locations_m.shape[0]):
                # One measurement per receiver
                m = rx.locations_m[i : i+1]
                n = rx.locations_n[i : i+1]
                new_rx = dc.receivers.Dipole(locations_m=m, locations_n=n, data_type=rx.data_type)

                # Build a new source with just this one receiver
                new_sources.append(dc.sources.Dipole([new_rx], src.location_a, src.location_b))

    # Create a new survey from the flattened sources
    return dc.Survey(new_sources)

Turning:

In [1]: dc_simulation.survey
Out[1]: Survey(#sources: 167; #data: 4757)

Into:

In [1]: dc_simulation.survey
Out[1]: Survey(#sources: 4757; #data: 4757)

But I think this is having a significant impact on computation time.

If anyone has any tips, or if I’m approaching this the wrong, way, the help would be greatly appreciated.

2

Hi @UTMZ17N ! Thanks for opening this thread.

I understand the issue you are running into. I can confirm that the workaround would work, but yes, it will have some impact on performance since you’ll be storing repeated fields for each repeated source location.

Nonetheless, I wasn’t able to reproduce it through a simple example like this one:

import numpy as np
import simpeg
from simpeg.electromagnetics.static import resistivity as dc
import discretize


# Define mesh
# ===========
hx = [(1.0, 40)]
hz = [(0.5, 80)]
h = [hx, hx, hz]
mesh = discretize.TensorMesh(h=h, origin="CCC")

_, _, z = mesh.cell_centers.T
active_cells = z <= 0
n_active_cells = active_cells.sum()


# Define survey
# =============
electrodes = np.array(
    [
        [-5.0, 0.0, 0.0],
        [-3.0, 0.0, 0.0],
        [-1.0, 0.0, 0.0],
        [1.0, 0.0, 0.0],
        [3.0, 0.0, 0.0],
        [5.0, 0.0, 0.0],
        [7.0, 0.0, 0.0],
    ]
)

receivers = [
    dc.receivers.Dipole(
        locations_m=electrodes[i], locations_n=electrodes[i + 1], data_type="volt"
    )
    for i in range(len(electrodes) - 1)
]

source = dc.sources.Dipole(
    receiver_list=receivers,
    location_a=np.array([-10.0, 0.0, 0.0]),
    location_b=np.array([-8.0, 0.0, 0.0]),
)

print("Number of receivers:", len(receivers))
for rx in source.receiver_list:
    print(rx)

survey = dc.Survey([source])


# Define simulation
# =================
air_conductivity = 1e-8
mapping = simpeg.maps.InjectActiveCells(
    mesh=mesh, active_cells=active_cells, value_inactive=air_conductivity
)
solver = simpeg.utils.get_default_solver()
simulation = dc.Simulation3DNodal(
    mesh=mesh, sigmaMap=mapping, survey=survey, solver=solver
)

# Run simulation
# ==============
conductivity_model = 1e-3 * np.ones(n_active_cells)
dpred = simulation.dpred(conductivity_model)

Here I’m defining a dipole source with 6 dipole receivers, and running a nodal simulation (although it also works with the cell centered one). The simulation runs without falling into the error you are experiencing. I’m running this script with SimPEG v0.25.1.

Could you provide more information about the survey, source and receivers in the survey you read from the file? Which version of SimPEG are you using (check it out with import simpeg; simpeg.__version__) ? Could you share the output of len(src.receiver_list)?

On a side note, you mention that your survey should be a Pole-Dipole array, but in the code I see that your source is a Dipole object. I doubt this is related to your issue, but just be aware that in order to run a Pole-Dipole, you might want to use the Pole source instead.

3

Hi @UTMZ17N. One more thing. Could you share which version of Numpy are you using when getting that error?

4

Hi @santisoler ,

Before I made this post, I attempted the fix from “Possible bug in BaseRx.eval() when using resistivity”, but to no avail, leading me to my flattening solution.

Since the error only occurred when I used apparent resistivity, I changed my input data to volts and things ran no problem.

Doing some further digging, I noticed the modifications I made to the “receivers.py” file were not located where I expected them on my drive. Long story short, I had two version of simpeg installed on my pc and I modified a version I was not using.

I went into the correct “receivers.py” file, made the suggested change and presto! Problem solved.

I changed

  • return v / self.geometric_factor[src]

to be

  • return v / self.geometric_factor[src][:, None]
5

Perfect! Glad to hear that. Thanks for reporting back.