How to use pymatsolver - 10 common examples

eta[blk_ind] = 0.1

        nElecs = 5
        x_temp = np.linspace(-250, 250, nElecs)
        aSpacing = x_temp[1]-x_temp[0]
        y_temp = 0.
        xyz = Utils.ndgrid(x_temp, np.r_[y_temp], np.r_[0.])
        srcList = DC.Utils.WennerSrcList(nElecs,aSpacing)
        survey = DC.SurveyIP(srcList)
        imap   = Maps.IdentityMap(mesh)
        problem = DC.ProblemIP(mesh, sigma=sigma, mapping= imap)
        problem.pair(survey)

        try:
            from pymatsolver import MumpsSolver
            problem.Solver = MumpsSolver
        except ImportError, e:
            problem.Solver = SolverLU

        mSynth = eta
        survey.makeSyntheticData(mSynth)

        # Now set up the problem to do some minimization
        dmis = DataMisfit.l2_DataMisfit(survey)
        reg = Regularization.Tikhonov(mesh)
        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
        inv = Inversion.BaseInversion(invProb)

        self.inv = inv
        self.reg = reg
        self.p = problem

srcs = []
    for ii in range(nSrc):
        rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20 + ii), rxType) for rxType in rxTypes.split(',')]
        srcs += [EM.TDEM.SrcTDEM_VMD_MVP(rxs,np.array([0., 0., 0.]))]

    survey = EM.TDEM.SurveyTDEM(srcs)

    prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
    # prb.timeSteps = [1e-5]
    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
    # prb.timeSteps = [(1e-05, 100)]

    try:
        from pymatsolver import PardisoSolver
        prb.Solver = PardisoSolver
    except ImportError:
        prb.Solver  = SolverLU

    sigma = np.ones(mesh.nCz)*1e-8
    sigma[mesh.vectorCCz<0] = 1e-1
    sigma = np.log(sigma[active])

    prb.pair(survey)
    return prb, mesh, sigma

def get_prob(mesh, mapping, formulation):
    prb = getattr(EM.TDEM, 'Problem3D_{}'.format(formulation))(
        mesh, sigmaMap=mapping
    )
    prb.timeSteps = [(1e-3, 5), (1e-4, 5), (5e-5, 10), (5e-5, 10), (1e-4, 10)]
    prb.Solver = Solver
    return prb

# Setup the problem object
    sigma1d = M.r(sigBG, 'CC', 'CC', 'M')[0, 0, :]

    if expMap:
        problem = Problem3D_ePrimSec(M, sigmaPrimary=np.log(sigma1d))
        problem.sigmaMap = simpeg.Maps.ExpMap(problem.mesh)
        problem.model = np.log(sig)
    else:
        problem = Problem3D_ePrimSec(M, sigmaPrimary=sigma1d)
        problem.sigmaMap = simpeg.Maps.IdentityMap(problem.mesh)
        problem.model = sig
    problem.pair(survey)
    problem.verbose = False
    try:
        from pymatsolver import Pardiso
        problem.Solver = Pardiso
    except:
        pass

    return (survey, problem)

#     freq=freq
            # ), # less accurate
        ]

        survey = EM.FDEM.Survey(SrcList)

        sig = 1e-1
        sigma = np.ones(mesh.nC)*sig
        sigma[mesh.gridCC[:, 2] > 0] = 1e-8

        prb = EM.FDEM.Problem3D_b(mesh, sigma=sigma)
        prb.pair(survey)

        try:
            from pymatsolver import Pardiso
            prb.Solver = Pardiso
        except ImportError:
            prb.Solver = SolverLU

        self.prb = prb
        self.mesh = mesh
        self.sig = sig

        print(' starting solve ...')
        u = self.prb.fields()
        print(' ... done')
        self.u = u

orientation=orientation))

        # primary
        self.primaryProblem = FDEM.Problem3D_j(meshp, sigmaMap=primaryMapping)
        self.primaryProblem.solver = Solver
        s_e = np.zeros(meshp.nF)
        inds = meshp.nFx + Utils.closestPoints(meshp, src_loc, gridLoc='Fz')
        s_e[inds] = 1./csz
        primarySrc = FDEM.Src.RawVec_e(
            self.rxlist, freq=freq, s_e=s_e/meshp.area
        )
        self.primarySurvey = FDEM.Survey([primarySrc])

        # Secondary Problem
        self.secondaryProblem = FDEM.Problem3D_e(meshs, sigmaMap=mapping)
        self.secondaryProblem.Solver = Solver
        self.secondarySrc = FDEM.Src.PrimSecMappedSigma(
                self.rxlist, freq, self.primaryProblem,
                self.primarySurvey, primaryMap2Meshs
        )
        self.secondarySurvey = FDEM.Survey([self.secondarySrc])
        self.secondaryProblem.pair(self.secondarySurvey)

        # Full 3D problem to compare with
        self.problem3D = FDEM.Problem3D_e(meshs, sigmaMap=mapping)
        self.problem3D.Solver = Solver
        s_e3D = np.zeros(meshs.nE)
        inds = (meshs.nEx + meshs.nEy +
                Utils.closestPoints(meshs, src_loc, gridLoc='Ez'))
        s_e3D[inds] = [1./(len(inds))] * len(inds)
        self.problem3D.model = model
        src3D = FDEM.Src.RawVec_e(self.rxlist, freq=freq, s_e=s_e3D)

)

    if srctype == "MagDipole":
        src = EM.TDEM.Src.MagDipole(
            [rx], waveform=EM.TDEM.Src.StepOffWaveform(),
            loc=np.array([0., 0., 0.])
        )
    elif srctype == "CircularLoop":
        src = EM.TDEM.Src.CircularLoop(
            [rx], waveform=EM.TDEM.Src.StepOffWaveform(),
            loc=np.array([0., 0., 0.]), radius=0.1
        )

    survey = EM.TDEM.Survey([src])
    prb = EM.TDEM.Problem3D_b(mesh, sigmaMap=mapping)
    prb.Solver = Solver

    prb.timeSteps = [(1e-06, 40), (5e-06, 40), (1e-05, 40), (5e-05, 40),
                     (0.0001, 40), (0.0005, 40)]

    sigma = np.ones(mesh.nCz)*1e-8
    sigma[active] = sig_half
    sigma = np.log(sigma[active])
    prb.pair(survey)
    if srctype == "MagDipole":
        bz_ana = mu_0*EM.Analytics.hzAnalyticDipoleT(rx.locs[0][0]+1e-3,
                                                     rx.times, sig_half)
    elif srctype == "CircularLoop":
        bz_ana = mu_0*EM.Analytics.hzAnalyticDipoleT(13, rx.times, sig_half)

    bz_calc = survey.dpred(sigma)
    ind = np.logical_and(rx.times > bounds[0], rx.times < bounds[1])

rxtimes = np.logspace(-4, -3, 20)
    rx = getattr(EM.TDEM.Rx, 'Point_{}'.format(rxcomp[:-1]))(
        locs=rxlocs, times=rxtimes, orientation=rxcomp[-1]
    )
    Aloc = np.r_[-10., 0., 0.]
    Bloc = np.r_[10., 0., 0.]
    srcloc = np.vstack((Aloc, Bloc))

    src = EM.TDEM.Src.LineCurrent([rx], loc=srcloc, waveform = EM.TDEM.Src.StepOffWaveform())
    survey = EM.TDEM.Survey([src])

    prb = getattr(EM.TDEM, 'Problem3D_{}'.format(prbtype))(mesh, sigmaMap=mapping)

    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]

    prb.Solver = Solver

    m = (
        np.log(1e-1)*np.ones(prb.sigmaMap.nP) +
        1e-3*np.random.randn(prb.sigmaMap.nP)
    )

    prb.pair(survey)
    mesh = mesh

    return prb, m, mesh

def setUpClass(self):

        print('\n------- Testing Primary Secondary Source HJ -> EB --------\n')
        # receivers
        self.rxlist = []
        for rxtype in ['b', 'e']:
            rx = getattr(FDEM.Rx, 'Point_{}'.format(rxtype))
            for orientation in ['x', 'y', 'z']:
                for comp in ['real', 'imag']:
                    self.rxlist.append(rx(rx_locs, component=comp,
                                       orientation=orientation))

        # primary
        self.primaryProblem = FDEM.Problem3D_j(meshp, sigmaMap=primaryMapping)
        self.primaryProblem.solver = Solver
        s_e = np.zeros(meshp.nF)
        inds = meshp.nFx + Utils.closestPoints(meshp, src_loc, gridLoc='Fz')
        s_e[inds] = 1./csz
        primarySrc = FDEM.Src.RawVec_e(
            self.rxlist, freq=freq, s_e=s_e/meshp.area
        )
        self.primarySurvey = FDEM.Survey([primarySrc])

        # Secondary Problem
        self.secondaryProblem = FDEM.Problem3D_e(meshs, sigmaMap=mapping)
        self.secondaryProblem.Solver = Solver
        self.secondarySrc = FDEM.Src.PrimSecMappedSigma(
                self.rxlist, freq, self.primaryProblem,
                self.primarySurvey, primaryMap2Meshs
        )
        self.secondarySurvey = FDEM.Survey([self.secondarySrc])

sigma0 = sigma*(1-eta)

        nElecs = 11
        x_temp = np.linspace(-100, 100, nElecs)
        aSpacing = x_temp[1]-x_temp[0]
        y_temp = 0.
        xyz = Utils.ndgrid(x_temp, np.r_[y_temp], np.r_[0.])
        srcList = DC.Utils.WennerSrcList(nElecs,aSpacing)
        survey = DC.SurveyDC(srcList)

        imap   = Maps.IdentityMap(mesh)
        problem = DC.ProblemDC_CC(mesh, mapping=imap)

        try:
            from pymatsolver import MumpsSolver
            solver = MumpsSolver
        except ImportError, e:
            solver = SolverLU

        problem.Solver = solver
        problem.pair(survey)

        phi0 = survey.dpred(sigma0)
        phiInf = survey.dpred(sigmaInf)

        phiIP_true = phi0-phiInf

        surveyIP = DC.SurveyIP(srcList)
        problemIP = DC.ProblemIP(mesh, sigma=sigma)
        problemIP.pair(surveyIP)

        problemIP.Solver = solver