Plotting the FDEM forward modelling data

Electromagnetics
1

Hello!
Thank you so much for the opportunity to use the SimPEG and ask questions!
I do 3D forward modelling for FDEM.

survey = FDEM.Survey(dg_p)

problem = FDEM.Problem3D_h(
mesh,
sigmaMap=Maps.IdentityMap(mesh),
Solver=Solver
)
problem.pair(survey)
t0 = time.time()
ModellingData = problem.fields(model)
h0 = ModellingData[dg_p, ‘h’]

Finally, I have h for each edge of my mesh.
Can I split them to hx, hy and hz?
How can I plot the data? I need maps for different depths.
Could you please help me?
Thank you in advance!

2 Likes
2

Welcome @Tatiana!

Here is a time-domain example where we set up plotting for the currents, magnetic fields and data measured at the surface. There are two functions in there that you could adapt for your purposes. The plot_currents function grabs the electric field (on edges) and plots them. In your example, the magnetic field h is on edges, so this function is very similar to what you need. Instead of the lines

S = np.kron(np.ones(3), sigma)
j = Utils.sdiag(S) * mesh.aveE2CCV * fields[srcList[iSrc], 'e', itime]

you can just do

hplot = mesh.aveE2CCV * h0

The operator mesh.aveE2CCV averages a vector (V) from edges (E) to cell-centers. Then we have the function

mesh.plotSlice(
    hplot, normal='Z', vType='CCv'
)

to plot a vector field (amplitudes are colors and arrows indicate the vector direction.

I hope this is a helpful starting point!

3

Dear @lheagy,
Thank you so much!
It is very useful for me!
Best regards,
Tatiana

1 Like
4

Hello Iheagy!

Could you please give me an example of the code for FDEM with horizontal current line?
I have plotted my results but it is quite strange for Hy component. I guess that the problem is connected with my source. I supposed I did some mistake. This is my code:3D_forward_modelling_line_in_mine.py (7.6 KB)

5

Hi @Tatiana, apologies for my very(!) late reply. It looks to me that in the code, the y-indices for the horizontal line current are not specified:

dgh_indx = ((mesh.gridFx[:, 0] > (A_x-csx*0.5)) &
                (mesh.gridFx[:, 0] <= (B_x+csx*0.5)))
dgh_indz = ((mesh.gridFx[:, 2] > (AB_z-csz*0.5)) & (mesh.gridFx[:, 2] <= (AB_z+csz*0.5)))
dgh_ind = dgh_indx & dgh_indz

without those, the line current would be defined everywhere in the plane that is defined by your x, z locations. So if you add a indy term, that should help

6

Thank you! Yes, it was my mistake