MT: 3D: ForwardΒΆ

Forward model 3D MT data.

(Source code)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
# Test script to use SimPEG.MT platform to forward model synthetic data.

# Import
import SimPEG as simpeg
from SimPEG.EM import NSEM
import numpy as np
try:
    from pymatsolver import Pardiso as Solver
except:
    from SimPEG import Solver


def run(plotIt=True, nFreq=1):
    """
        MT: 3D: Forward
        ===============

        Forward model 3D MT data.

    """

    # Make a mesh
    M = simpeg.Mesh.TensorMesh(
        [
            [(100, 5, -1.5), (100., 10), (100, 5, 1.5)],
            [(100, 5, -1.5), (100., 10), (100, 5, 1.5)],
            [(100, 5, +1.6), (100., 10), (100, 3, 2)]
        ], x0=['C', 'C', -3529.5360]
    )
    # Setup the model
    conds = [1e-2, 1]
    sig = simpeg.Utils.ModelBuilder.defineBlock(
        M.gridCC, [-1000, -1000, -400], [1000, 1000, -200], conds
    )
    sig[M.gridCC[:, 2] > 0] = 1e-8
    sig[M.gridCC[:, 2] < -600] = 1e-1
    sigBG = np.zeros(M.nC) + conds[0]
    sigBG[M.gridCC[:, 2] > 0] = 1e-8

    # Setup the the survey object
    # Receiver locations
    rx_x, rx_y = np.meshgrid(np.arange(-500, 501, 50), np.arange(-500, 501, 50))
    rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x, 2), simpeg.Utils.mkvc(rx_y, 2), np.zeros((np.prod(rx_x.shape), 1))))
    # Make a receiver list
    rxList = []
    for loc in rx_loc:
        # NOTE: loc has to be a (1, 3) np.ndarray otherwise errors accure
        for rx_orientation in ['xx', 'xy', 'yx', 'yy']:
            rxList.append(NSEM.Rx.Point_impedance3D(simpeg.mkvc(loc, 2).T, rx_orientation, 'real'))
            rxList.append(NSEM.Rx.Point_impedance3D(simpeg.mkvc(loc, 2).T, rx_orientation, 'imag'))
        for rx_orientation in ['zx', 'zy']:
            rxList.append(NSEM.Rx.Point_tipper3D(simpeg.mkvc(loc, 2).T, rx_orientation, 'real'))
            rxList.append(NSEM.Rx.Point_tipper3D(simpeg.mkvc(loc, 2).T, rx_orientation, 'imag'))
    # Source list
    srcList = [
        NSEM.Src.Planewave_xy_1Dprimary(rxList, freq)
        for freq in np.logspace(3, -3, nFreq)
    ]
    # Survey MT
    survey = NSEM.Survey(srcList)

    # Setup the problem object
    problem = NSEM.Problem3D_ePrimSec(M, sigma=sig, sigmaPrimary=sigBG)

    problem.pair(survey)
    problem.Solver = Solver

    # Calculate the data
    fields = problem.fields()
    dataVec = survey.eval(fields)

    # Make the data
    mtData = NSEM.Data(survey, dataVec)

    # Add plots
    if plotIt:
        pass

if __name__ == '__main__':
    run()