Python PETScIO.arrayToVec方法代码示例

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]
def RemoveRowCol(AA,bb,VelPres):
    As = AA.sparray()
    As.eliminate_zeros()
    Adelete = remove_ij(As,VelPres-1,VelPres-1)
    A = PETSc.Mat().createAIJ(size=Adelete.shape,csr=(Adelete.indptr, Adelete.indices, Adelete.data))
    # StoreMatrix(Adelete, "As")
    b = np.delete(bb,VelPres-1,0)
    zeros = 0*b
    bb= IO.arrayToVec(b)
    x = IO.arrayToVec(zeros)
    return A,bb,x

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]
def SchurPCD(Mass,L,F, backend):
    Mass = Mass.sparray()
    F = F.sparray()
    F = F + 1e-10*sp.identity(Mass.shape[0])
    F = PETSc.Mat().createAIJ(size=F.shape,csr=(F.indptr, F.indices, F.data))
    Mass.tocsc()
    Schur = sp.rand(Mass.shape[0], Mass.shape[0], density=0.00, format='csr')
    ksp = PETSc.KSP().create()
    pc = ksp.getPC()
    ksp.setOperators(F,F)
    ksp.setType('preonly')
    pc.setType('lu')
    OptDB = PETSc.Options()
    OptDB['pc_factor_shift_amount'] = "0.1"
    # OptDB['pc_factor_shift_type'] = 'POSITIVE_DEFINITE'
    OptDB['pc_factor_mat_ordering_type'] = 'amd'
    # OptDB['rtol']  = 1e-8
    # ksp.max_it = 5
    ksp.setFromOptions()
    for i in range(0,Mass.shape[0]):
        Col = Mass.getcol(i)
        Col = Col.toarray()
        Col = IO.arrayToVec(Col)
        u = Col.duplicate()
        ksp.solve(Col,u)
        C = u.duplicate()
        L.mult(u,C)
        # print C.array
        Schur[i,:] = C.array

    if backend == "PETSc":
        return PETSc.Mat().createAIJ(size=Schur.transpose().shape,csr=(Schur.transpose().indptr, Schur.transpose().indices, Schur.transpose().data))
    else:
        return Schur.transpose()

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]
def Assemble(W, NS, Maxwell, Couple, L_ns, L_m, RHSform, BC, Type, IterType):
    
    tic() 
    if Type == 'NonLinear':
        F = assemble(NS[0])
        BC[0].apply(F)
        F = F.sparray()
        if IterType == 'Full':
            C = assemble(Couple[0])
            C = BC[4]*C.sparray()*BC[3]
        else:
            C = None
        bu = assemble(L_ns-RHSform[0])
        bp = assemble(-RHSform[1])
        bb = assemble(L_m-RHSform[2])
        br = assemble(-RHSform[3])
        BC[0].apply(bu)
        BC[1].apply(bb)
        BC[2].apply(br)
        b = np.concatenate((bu.array(),bp.array(),bb.array(),br.array()),axis = 0)

        MO.StrTimePrint("MHD non-linear matrix assembled, time: ",toc())
        return [F, C],b
    elif Type == 'Linear':

        M = assemble(Maxwell[0])
        D = assemble(Maxwell[2])
        SS = assemble(Maxwell[3])

        B = assemble(NS[2])
        S = assemble(NS[3])


        SS = 0*SS
        BC[1].apply(M)
        BC[2].apply(SS)  

        
        B = B.sparray()*BC[3]
        S = S.sparray()

        M = M.sparray()
        D = BC[4]*D.sparray()*BC[5]
        SS = SS.sparray()
        
        MO.StrTimePrint("MHD linear matrix assembled, time: ",toc())
        return [B,M,D,S,SS]
    else:
        bu = assemble(L_ns-RHSform[0])
        bp = assemble(-RHSform[1])
        bb = assemble(L_m-RHSform[2])
        br = assemble(-RHSform[3])
        BC[0].apply(bu)
        BC[1].apply(bb)
        BC[2].apply(br)
        b = np.concatenate((bu.array(),bp.array(),bb.array(),br.array()),axis = 0)
        return IO.arrayToVec(b)

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]
def BCapply(V,BC,x,opt = "PETSc"):
    v = Function(V)
    v.vector()[:] = x.array
    BC.apply(v.vector())
    if opt == "PETSc":
        x = IO.arrayToVec(v.vector().array())
        return x
    else:
        return v

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]
def SystemAssemble(W,A,b,SetupType,IterType):
    tic()
    if SetupType == 'Matrix':
        for i in range(len(A)):
            if A[i] != None:
                A[i].eliminate_zeros()
        if IterType == 'Full':
            A = CP.Scipy2PETSc(bmat([[A[0],A[2].T,-A[1].T,None],
                      [A[2],A[5],None,None],
                      [A[1],None,A[3],A[4]],
                      [None,None,A[4].T,A[6]]]))
        else:
            A = CP.Scipy2PETSc(bmat([[A[0],A[2].T,None,None],
                      [A[2],A[5],None,None],
                      [None,None,A[3],A[4]],
                      [None,None,A[4].T,A[6]]]))
        b = IO.arrayToVec(b)
        MO.StrTimePrint("MHD system assemble, time: ",toc())
        return A,b
    else:
        for i in range(len(A)):
            if A[i] != None:
                A[i] = CP.Scipy2PETSc(A[i])
        if IterType == 'Full':
            P = PETSc.Mat().createPython([W[0].dim()+W[1].dim()+W[2].dim()+W[3].dim(),W[0].dim()+W[1].dim()+W[2].dim()+W[3].dim()])
            P.setType('python')
            p = MHDmulti.MHDmat(W,A)
            P.setPythonContext(p)
        else:
            MatFluid = PETSc.Mat().createPython([W[0].dim()+W[1].dim(), W[0].dim()+W[1].dim()])
            MatFluid.setType('python')
            pFluid = MHDmulti.MatFluid([W[0],W[1]],A)
            MatFluid.setPythonContext(pFluid)

            MatMag = PETSc.Mat().createPython([W[2].dim()+W[3].dim(), W[2].dim()+W[3].dim()])
            MatMag.setType('python')
            pMag = MHDmulti.MatMag([W[2],W[3]],A)
            MatMag.setPythonContext(pMag)
            P = [MatFluid,MatMag]
        b = IO.arrayToVec(b)
        MO.StrTimePrint("MHD mult-class setup, time: ",toc())
        return P,b

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]

#.........这里部分代码省略.........
                
                # OptDB = PETSc.Options()

                PP = PETSc.Mat().create()
                PP.setSizes([A.size[0], A.size[0]])

                #PP = PETSc.Mat().createPython([A.size[0], A.size[0]])
                PP.setType('python')
                pp = MHDmulti.P(FFSS,A,MatrixLinearFluids[1],MatrixLinearFluids[0],kspFp,HiptmairMatrices[6])
                
                PP.setPythonContext(pp)
                PP.setUp()
                kspInner.setOperators(PP,PP)
                kspInner.setFromOptions()
                scale = b.norm()
                b = b/scale
                del A
                stime = time.time()
                kspOuter.solve(b,u)
                Soltime = time.time()- stime
                NSits += kspOuter.its
                Mits += kspInner.its
                u = u*scale
                SolutionTime = SolutionTime +Soltime
                MO.PrintStr("Number of outer iterations ="+str(kspOuter.its),60,"+","\n\n","\n\n")
                MO.PrintStr("Number of inner iterations ="+str(kspInner.its),60,"+","\n\n","\n\n")
            u1, p1, b1, r1, eps= Iter.PicardToleranceDecouple(u,x,FSpaces,dim,"2",iter, SaddlePoint = "Yes")
            p1.vector()[:] += - assemble(p1*dx)/assemble(ones*dx)
            u_k.assign(u1)
            p_k.assign(p1)
            b_k.assign(b1)
            r_k.assign(r1)
            uOld= np.concatenate((u_k.vector().array(),b_k.vector().array(),p_k.vector().array(),r_k.vector().array()), axis=0)
            x = IO.arrayToVec(uOld)



        XX= np.concatenate((u_k.vector().array(),p_k.vector().array(),b_k.vector().array(),r_k.vector().array()), axis=0)
        SolTime[xx-1] = SolutionTime/iter
        NSave[xx-1] = (float(NSits)/iter)
        Mave[xx-1] = (float(Mits)/iter)
        iterations[xx-1] = iter
        TotalTime[xx-1] = time.time() - TotalStart
        dim = [Velocity.dim(), Pressure.dim(), Magnetic.dim(),Lagrange.dim()]
#
#        ExactSolution = [u0,p0,b0,r0]
#        errL2u[xx-1], errH1u[xx-1], errL2p[xx-1], errL2b[xx-1], errCurlb[xx-1], errL2r[xx-1], errH1r[xx-1] = Iter.Errors(XX,mesh,FSpaces,ExactSolution,order,dim, "DG")
#
#        if xx > 1:
#            l2uorder[xx-1] =  np.abs(np.log2(errL2u[xx-2]/errL2u[xx-1]))
#            H1uorder[xx-1] =  np.abs(np.log2(errH1u[xx-2]/errH1u[xx-1]))
#
#            l2porder[xx-1] =  np.abs(np.log2(errL2p[xx-2]/errL2p[xx-1]))
#
#            l2border[xx-1] =  np.abs(np.log2(errL2b[xx-2]/errL2b[xx-1]))
#            Curlborder[xx-1] =  np.abs(np.log2(errCurlb[xx-2]/errCurlb[xx-1]))
#
#            l2rorder[xx-1] =  np.abs(np.log2(errL2r[xx-2]/errL2r[xx-1]))
#            H1rorder[xx-1] =  np.abs(np.log2(errH1r[xx-2]/errH1r[xx-1]))




    import pandas as pd

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]
def u_prev(u,p,b,r):
    uOld = np.concatenate((u.vector().array(),p.vector().array(),b.vector().array(),r.vector().array()), axis=0)
    x = IO.arrayToVec(uOld)
    return x

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]
def FuncToPETSc(x):
    return IO.arrayToVec(x.vector().array())

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]
    # <codecell>

    C = sparse.csr_matrix((V.dim(),Q.dim()))
    (v) = TrialFunction(V)
    (u) = TestFunction(V)
    tic()
    for i in range(0,Q.dim()):
        uOut = Function(V)
        uu = Function(Q)
        x = M.getVecRight()
        zero = np.zeros((Q.dim(),1))[:,0]
        zero[i] = 1
        uu.vector()[:] = zero
        L = assemble(inner(u, grad(uu))*dx)
        rhs = IO.arrayToVec(B[:,i].toarray())
        ksp.solve(rhs,x)
    #     x = project(grad(uu),V)
        P = x.array
        uOut.vector()[:] = P
        low_values_indices = np.abs(P) < 1e-3
        P[low_values_indices] = 0
        P=np.around(P)
        pn = P.nonzero()[0]
        for j in range(0,len(pn)):
            C[pn[j],i] = P[pn[j]]
        del uu
    print toc()
    print C.todense()
    # C = sparse.csr_matrix((Magnetic.dim(),Lagrange.dim()))
    # Mmap = Magnetic.dofmap()

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]

#.........这里部分代码省略.........
#                            A,b = MHDsetup.SystemAssemble(FSpaces,A,b,SetupType,IterType)
#                            u = b.duplicate()
#                        else: 
#                            Fnline,b = MHDsetup.Assemble(W,ns,maxwell,CoupleTerm,Lns,Lmaxwell,RHSform,bcs+BC, "NonLinear",IterType)
#                            A = Fnlin+Alin
#                            A,b = MHDsetup.SystemAssemble(FSpaces,A,b,SetupType,IterType)
                        AA, bb = assemble_system(maxwell+ns+CoupleTerm, (Lmaxwell + Lns) - RHSform,  bcs)
                        A,b = CP.Assemble(AA,bb)
                    # if iter == 1:
                    MO.StrTimePrint("MHD total assemble, time: ", time.time()-AssembleTime)
                    
                    u = b.duplicate()
                    #A,Q
                    kspFp, Fp = PrecondSetup.FluidNonLinearSetup(Pressure, MU, u_k)
                    print "Inititial guess norm: ", u.norm()
                    if u.norm()>1e50:
                        iter = 10000
                        break
                    stime = time.time()
                    kspF = 0
                    u, mits,nsits = S.solve(A,b,u,params,W,'Direct',IterType,OuterTol,InnerTol,HiptmairMatrices,Hiptmairtol,KSPlinearfluids, Fp,kspF)
                    Soltime = time.time()- stime
                    Mits += mits
                    NSits += nsits
                    SolutionTime += Soltime
                    
                    u1, p1, b1, r1, eps= Iter.PicardToleranceDecouple(u,x,FSpaces,dim,"2",iter)
                    p1.vector()[:] += - assemble(p1*dx)/assemble(ones*dx)
                    u_k.assign(u1)
                    p_k.assign(p1)
                    b_k.assign(b1)
                    r_k.assign(r1)
                    uOld= np.concatenate((u_k.vector().array(),p_k.vector().array(),b_k.vector().array(),r_k.vector().array()), axis=0)
                    x = IO.arrayToVec(uOld)



                XX= np.concatenate((u_k.vector().array(),p_k.vector().array(),b_k.vector().array(),r_k.vector().array()), axis=0)

                iterations[xx-1,qq] = iter
                dim = [Velocity.dim(), Pressure.dim(), Magnetic.dim(),Lagrange.dim()]
    #
#        ExactSolution = [u0,p0,b0,r0]
#        errL2u[xx-1], errH1u[xx-1], errL2p[xx-1], errL2b[xx-1], errCurlb[xx-1], errL2r[xx-1], errH1r[xx-1] = Iter.Errors(XX,mesh,FSpaces,ExactSolution,order,dim, "DG")
#
#        if xx > 1:
#            l2uorder[xx-1] =  np.abs(np.log2(errL2u[xx-2]/errL2u[xx-1]))
#            H1uorder[xx-1] =  np.abs(np.log2(errH1u[xx-2]/errH1u[xx-1]))
#
#            l2porder[xx-1] =  np.abs(np.log2(errL2p[xx-2]/errL2p[xx-1]))
#
#            l2border[xx-1] =  np.abs(np.log2(errL2b[xx-2]/errL2b[xx-1]))
#            Curlborder[xx-1] =  np.abs(np.log2(errCurlb[xx-2]/errCurlb[xx-1]))
#
#            l2rorder[xx-1] =  np.abs(np.log2(errL2r[xx-2]/errL2r[xx-1]))
#            H1rorder[xx-1] =  np.abs(np.log2(errH1r[xx-2]/errH1r[xx-1]))
#
#
#
#
#    import pandas as pd
#
#
#
#    LatexTitles = ["l","DoFu","Dofp","V-L2","L2-order","V-H1","H1-order","P-L2","PL2-order"]
#    LatexValues = np.concatenate((level,Velocitydim,Pressuredim,errL2u,l2uorder,errH1u,H1uorder,errL2p,l2porder), axis=1)

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]
def PETScToFunc(V,x):
    v = Function(V)
    v.vector()[:] = x.array
    return x

def FuncToPETSc(x):
    return IO.arrayToVec(x.vector().array())


Hdiv = FunctionSpace(mesh,'BDM',order)

f = Expression(('(x[0])','(x[1])'))

Ft = interpolate(f,Magnetic)
x = IO.arrayToVec(Ft.vector().array())
print x.array
Pxx = Px.createVecRight()
Px.multTranspose(x,Pxx)
Pyy  = Py.createVecRight()
Py.multTranspose(x,Pyy)
PPP = CP.PETSc2Scipy(Px)
print (PPP*PPP.T).nnz
print (PPP*PPP.T).diagonal()

MO.StoreMatrix(PPP,"P")

P = np.concatenate((Pxx.array,Pyy.array),axis=1)
# print P
f = BCapply(Magnetic,BC[0],x,"dolfin")
fVec = interpolate(f,VecLagrange)

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]
    OptDB = PETSc.Options()
    ksp.setFromOptions()
    C = sparse.csr_matrix((V.dim(),Q.dim()))

    tic()
    for i in range(0,Q.dim()):
        uOut = Function(V)
        uu = Function(Q)
        x = M.getVecRight()
        zero = np.zeros((Q.dim(),1))[:,0]
        zero[i] = 1
        uu.vector()[:] = zero
        L = assemble(inner(u, grad(uu))*dx)
        # bcb.apply(L)
        rhs = IO.arrayToVec(L.array())
        ksp.solve(rhs,x)
    #     x = project(grad(uu),V)
        P = x.array
        uOut.vector()[:] = P
        low_values_indices = np.abs(P) < 1e-3
        P[low_values_indices] = 0
        #P=np.around(P)
        pn = P.nonzero()[0]
        for j in range(0,len(pn)):
            C[pn[j],i] = P[pn[j]]
        del uu
    print toc()

    pathToGrad = "/home/mwathen/Dropbox/MastersResearch/MHD/FEniCS/GradMatrices/"
    name = pathToGrad+"UnitSquareCrossed_m="+str(nn)

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]
def Stokes(V,Q,BC,f,mu,boundaries):

    parameters = CP.ParameterSetup()

    W = V*Q

    (u, p) = TrialFunctions(W)
    (v, q) = TestFunctions(W)

    def boundary(x, on_boundary):
        return on_boundary


    bcu1 = DirichletBC(W.sub(0),BC[0], boundaries,2)
    bcu2 = DirichletBC(W.sub(0),BC[1], boundaries,1)
    bcs = [bcu1,bcu2]
    u_k = Function(V)
    a11 = mu*inner(grad(v), grad(u))*dx
    a12 = div(v)*p*dx
    a21 = div(u)*q*dx
    L1  =  inner(v,f)*dx
    a = mu*a11-a12-a21
    i = p*q*dx

    tic()
    AA, bb = assemble_system(a, L1, bcs)
    # A = as_backend_type(AA).mat()
    A,b = CP.Assemble(AA,bb)
    print toc()
    b = bb.array()
    zeros = 0*b
    del bb
    bb = IO.arrayToVec(b)
    x = IO.arrayToVec(zeros)


    pp = inner(grad(v), grad(u))*dx + (1/mu)*p*q*dx
    PP, Pb = assemble_system(pp,L1,bcs)
    P = CP.Assemble(PP)


    u_is = PETSc.IS().createGeneral(range(V.dim()))
    p_is = PETSc.IS().createGeneral(range(V.dim(),V.dim()+Q.dim()))

    ksp = PETSc.KSP().create()
    ksp.setOperators(A,P)
    pc = ksp.getPC()
    pc.setType(pc.Type.FIELDSPLIT)
    fields = [ ("field1", u_is), ("field2", p_is)]
    pc.setFieldSplitIS(*fields)
    pc.setFieldSplitType(0)

    OptDB = PETSc.Options()

    OptDB["field_split_type"] = "multiplicative"

    OptDB["fieldsplit_field1_ksp_type"] = "preonly"
    OptDB["fieldsplit_field1_pc_type"] = "hypre"
    OptDB["fieldsplit_field2_ksp_type"] = "cg"
    OptDB["fieldsplit_field2_pc_type"] = "jacobi"
    OptDB["fieldsplit_field2_ksp_rtol"] = 1e-8



    ksp.setFromOptions()
    ksp.setTolerances(1e-8)


    ksp.solve(bb, x)


    X = IO.vecToArray(x)
    x = X[0:V.dim()]
    p =  X[V.dim():]
    # x =
    u = Function(V)
    u.vector()[:] = x

    p_k = Function(Q)
    n = p.shape
    p_k.vector()[:] = p


    u_k.assign(u)


    return u_k, p_k

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]
 def Amatvec(v):
     x = IO.arrayToVec(v)
     B = x.duplicate()
     AA.mult(x,B)
     b = IO.vecToArray(B)
     return b

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import arrayToVec [as 别名]

#.........这里部分代码省略.........
            reshist = {}
            def monitor(ksp, its, fgnorm):
                reshist[its] = fgnorm
                print its, "    OUTER:", fgnorm
            # ksp.setMonitor(monitor)
            ksp.max_it = 100
            ksp.setTolerances(OuterTol)

            W = Fspace
            FFSS = [W.sub(0), W.sub(1), W.sub(2), W.sub(3)]

            #pc.setPythonContext(MHDprec.BlkInvA(FFSS, kspF, KSPlinearfluids[0], KSPlinearfluids[1], Fp, HiptmairMatrices[
            #                    3], HiptmairMatrices[4], HiptmairMatrices[2], HiptmairMatrices[0], HiptmairMatrices[1], HiptmairMatrices[6], Hiptmairtol))
            pc.setPythonContext(MHDprec.ApproxInv(FFSS, kspF, KSPlinearfluids[0], KSPlinearfluids[1], Fp, HiptmairMatrices[
                                3], HiptmairMatrices[4], HiptmairMatrices[2], HiptmairMatrices[0], HiptmairMatrices[1], HiptmairMatrices[6], Hiptmairtol))

            # OptDB = PETSc.Options()
            # OptDB['pc_factor_mat_solver_package']  = "umfpack"
            # OptDB['pc_factor_mat_ordering_type']  = "rcm"
            ksp.setFromOptions()
            scale = b.norm()
            b = b/scale
            ksp.setOperators(A, A)
            del A
            ksp.solve(b, u)
            # Mits +=dodim
            u = u*scale
            MO.PrintStr("Number iterations = "+str(ksp.its),
                        60, "+", "\n\n", "\n\n")
            return u, ksp.its, 0

        IS = MO.IndexSet(Fspace, '2by2')
        M_is = IS[1]
        NS_is = IS[0]
        kspNS = PETSc.KSP().create()
        kspM = PETSc.KSP().create()
        kspNS.setTolerances(OuterTol)

        kspNS.setOperators(A[0])
        kspM.setOperators(A[1])
        # print P.symmetric
        if IterType == "MD":
            kspNS.setType('gmres')
            kspNS.max_it = 500

            pcNS = kspNS.getPC()
            pcNS.setType(PETSc.PC.Type.PYTHON)
            pcNS.setPythonContext(NSpreconditioner.NSPCD(MixedFunctionSpace(
                [Fspace.sub(0), Fspace.sub(1)]), kspF, KSPlinearfluids[0], KSPlinearfluids[1], Fp))
        elif IterType == "CD":
            kspNS.setType('minres')
            pcNS = kspNS.getPC()
            pcNS.setType(PETSc.PC.Type.PYTHON)
            Q = KSPlinearfluids[1].getOperators()[0]
            Q = 1./params[2]*Q
            KSPlinearfluids[1].setOperators(Q, Q)
            pcNS.setPythonContext(StokesPrecond.MHDApprox(MixedFunctionSpace(
                [Fspace.sub(0), Fspace.sub(1)]), kspF, KSPlinearfluids[1]))
        reshist = {}
        def monitor(ksp, its, fgnorm):
            reshist[its] = fgnorm
            print fgnorm
        # kspNS.setMonitor(monitor)

        uNS = u.getSubVector(NS_is)
        bNS = b.getSubVector(NS_is)
        # print kspNS.view()
        scale = bNS.norm()
        bNS = bNS/scale
        print bNS.norm()
        kspNS.solve(bNS, uNS)
        uNS = uNS*scale
        NSits = kspNS.its
        kspNS.destroy()
        # for line in reshist.values():
        #     print line
        kspM.setFromOptions()
        kspM.setType(kspM.Type.MINRES)
        kspM.setTolerances(InnerTol)
        pcM = kspM.getPC()
        pcM.setType(PETSc.PC.Type.PYTHON)
        pcM.setPythonContext(MP.Hiptmair(MixedFunctionSpace([Fspace.sub(2), Fspace.sub(3)]), HiptmairMatrices[3], HiptmairMatrices[
                             4], HiptmairMatrices[2], HiptmairMatrices[0], HiptmairMatrices[1], HiptmairMatrices[6], Hiptmairtol))

        uM = u.getSubVector(M_is)
        bM = b.getSubVector(M_is)
        scale = bM.norm()
        bM = bM/scale
        print bM.norm()
        kspM.solve(bM, uM)
        uM = uM*scale
        Mits = kspM.its
        kspM.destroy()
        u = IO.arrayToVec(np.concatenate([uNS.array, uM.array]))

        MO.PrintStr("Number of M iterations = " +
                    str(Mits), 60, "+", "\n\n", "\n\n")
        MO.PrintStr("Number of NS/S iterations = " +
                    str(NSits), 60, "+", "\n\n", "\n\n")
        return u, NSits, Mits