How to use the pymatsolver.MumpsSolver function in pymatsolver

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

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

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

        try:
            import pyMKL
            def pardisoSolverFunc(A):
                from pyMKL import pardisoSolver
                factor = pardisoSolver(A, mtype=2) #Real and SPD
                factor.run_pardiso(12) #analysis and factor
                return factor
            self.prb.Solver = SolverWrapD(pardisoSolverFunc)
        except ImportError:
            try:
                from pymatsolver import MumpsSolver
                self.prb.Solver = MumpsSolver
            except ImportError:
                self.prb.Solver = SolverLU

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

        self.prb.pair(survey)
        self.mesh = mesh

survey = MT.Survey(srcList)

    ## Setup the problem object
    sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
    if expMap:
        problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= np.log(sigma1d) )
        problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
        problem.curModel = np.log(sig)
    else:
        problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= sigma1d)
        problem.curModel = sig
    problem.pair(survey)
    problem.verbose = False
    try:
        from pymatsolver import MumpsSolver
        problem.Solver = MumpsSolver
    except:
        pass

    return (survey, problem)

problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= sigma1d)
        problem.curModel = sig
    problem.pair(survey)
    problem.verbose = False
    try:
        import pyMKL
        def pardisoSolverFunc(A):
            from pyMKL import pardisoSolver
            factor = pardisoSolver(A, mtype=6) #Complex and symmetric
            factor.run_pardiso(12) #analysis and factor
            return factor
        problem.Solver = SolverWrapD(pardisoSolverFunc)
    except ImportError:
        try:
            from pymatsolver import MumpsSolver
            problem.Solver = MumpsSolver
        except ImportError:
            pass

    return (survey, problem)

from scipy.sparse.linalg import dsolve
from SimPEG.FLOW import Richards
from SimPEG import SolverWrapD

try:
    import pyMKL
    def pardisoSolverFunc(A):
        from pyMKL import pardisoSolver
        factor = pardisoSolver(A, mtype=11) #Real and not symmetric
        factor.run_pardiso(12) #analysis and factor
        return factor
    Solver = SolverWrapD(pardisoSolverFunc)
except ImportError:
    try:
        from pymatsolver import MumpsSolver
        Solver = MumpsSolver
    except ImportError:
        pass


TOL = 1E-8

np.random.seed(0)


class TestModels(unittest.TestCase):

    def test_BaseHaverkamp_Theta(self):
        mesh = Mesh.TensorMesh([50])
        hav = Richards.Empirical._haverkamp_theta(mesh)
        m = np.random.randn(50)
        def wrapper(u):

mu[blkind] = mu_0*1.1
    offset = 0.
    frequency = np.r_[10., 100., 1000.]
    rx0 = EM.FDEM.Rx(np.array([[8., 0., 30.]]), 'bzr')
    rx1 = EM.FDEM.Rx(np.array([[8., 0., 30.]]), 'bzi')
    srcLists = []
    nfreq = frequency.size
    for ifreq in range(nfreq):
        src = EM.FDEM.Src.CircularLoop([rx0, rx1], frequency[ifreq], np.array([[0., 0., 30.]]), radius=5.)
        srcLists.append(src)
    survey = EM.FDEM.Survey(srcLists)
    iMap = Maps.IdentityMap(nP=int(mesh.nC))
    # Use PhysPropMap
    maps = [('sigma', iMap), ('mu', iMap)]
    prob = EM.FDEM.Problem_b(mesh, mapping=maps)
    prob.Solver = MumpsSolver
    survey.pair(prob)
    m = np.r_[sigma, mu]
    survey0 = EM.FDEM.Survey(srcLists)
    prob0 = EM.FDEM.Problem_b(mesh, mapping=maps)
    prob0.Solver = MumpsSolver
    survey0.pair(prob0)
    m = np.r_[sigma, mu]
    m0 = np.r_[sigma, np.ones(mesh.nC)*mu_0]
    m00 = np.r_[np.ones(mesh.nC)*1e-8, np.ones(mesh.nC)*mu_0]
    # Anomalous conductivity and susceptibility
    F = prob.fields(m)
    # Only anomalous conductivity
    F0 = prob.fields(m0)
    # Primary field
    F00 = prob.fields(m00)

# Import
import SimPEG as simpeg
from SimPEG import MT, SolverWrapD
import numpy as np
try:
    import pyMKL
    def pardisoSolverFunc(A):
        from pyMKL import pardisoSolver
        factor = pardisoSolver(A, mtype=6) #Complex and symmetric
        factor.run_pardiso(12) #analysis and factor
        return factor
    Solver = SolverWrapD(pardisoSolverFunc)
except ImportError:
    try:
        from pymatsolver import MumpsSolver
        Solver = MumpsSolver
    except ImportError:
        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]

src = DC.Src.Dipole([rx], np.r_[-200, 0, -12.5], np.r_[+200, 0, -12.5])
    survey = DC.Survey([src])
    problem = DC.Problem3D_CC(mesh)
    problem.pair(survey)
    try:
        import pyMKL
        def pardisoSolverFunc(A):
            from pyMKL import pardisoSolver
            factor = pardisoSolver(A, mtype=11) #Real and not symmetric
            factor.run_pardiso(12) #analysis and factor
            return factor
        problem.Solver = SolverWrapD(pardisoSolverFunc)
    except ImportError:
        try:
            from pymatsolver import MumpsSolver
            problem.Solver = MumpsSolver
        except ImportError:
            pass
    data = survey.dpred(sigma)

    def DChalf(srclocP, srclocN, rxloc, sigma, I=1.):
        rp = (srclocP.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
        rn = (srclocN.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
        rP = np.sqrt(((rxloc-rp)**2).sum(axis=1))
        rN = np.sqrt(((rxloc-rn)**2).sum(axis=1))
        return I/(sigma*2.*np.pi)*(1/rP-1/rN)

    data_anaP = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxP, sighalf)
    data_anaN = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxN, sighalf)
    data_ana = data_anaP-data_anaN
    Data_ana = data_ana.reshape((21, 21), order = 'F')
    Data = data.reshape((21, 21), order = 'F')

srcLists = []
    nfreq = frequency.size
    for ifreq in range(nfreq):
        src = EM.FDEM.Src.CircularLoop([rx0, rx1], frequency[ifreq], np.array([[0., 0., 30.]]), radius=5.)
        srcLists.append(src)
    survey = EM.FDEM.Survey(srcLists)
    iMap = Maps.IdentityMap(nP=int(mesh.nC))
    # Use PhysPropMap
    maps = [('sigma', iMap), ('mu', iMap)]
    prob = EM.FDEM.Problem_b(mesh, mapping=maps)
    prob.Solver = MumpsSolver
    survey.pair(prob)
    m = np.r_[sigma, mu]
    survey0 = EM.FDEM.Survey(srcLists)
    prob0 = EM.FDEM.Problem_b(mesh, mapping=maps)
    prob0.Solver = MumpsSolver
    survey0.pair(prob0)
    m = np.r_[sigma, mu]
    m0 = np.r_[sigma, np.ones(mesh.nC)*mu_0]
    m00 = np.r_[np.ones(mesh.nC)*1e-8, np.ones(mesh.nC)*mu_0]
    # Anomalous conductivity and susceptibility
    F = prob.fields(m)
    # Only anomalous conductivity
    F0 = prob.fields(m0)
    # Primary field
    F00 = prob.fields(m00)

    if plotIt:
        import matplotlib.pyplot as plt
        def vizfields(ifreq=0, primsec="secondary",realimag="real"):

            titles = ["F[$\sigma$, $\mu$]", "F[$\sigma$, $\mu_0$]", "F[$\sigma$, $\mu$]-F[$\sigma$, $\mu_0$]"]