Python PETScIO类代码示例

def PicardToleranceDecouple(x,U,FSpaces,dim,NormType,iter):
    X = IO.vecToArray(x)
    uu = X[0:dim[0]]
    pp = X[dim[0]:dim[0]+dim[1]-1]
    bb = X[dim[0]+dim[1]-1:dim[0]+dim[1]+dim[2]-1]
    rr = X[dim[0]+dim[1]+dim[2]-1:]

    u = Function(FSpaces[0])
    u.vector()[:] = uu
    diffu = u.vector().array()

    p = Function(FSpaces[1])
    n = pp.shape
    p.vector()[:] = np.insert(pp,n,0)
    # ones = Function(FSpaces[1])
    # ones.vector()[:]=(0*ones.vector().array()+1)
    # pp = Function(FSpaces[1])
    # p.vector()[:] = p.vector().array()- assemble(p*dx)/assemble(ones*dx)

    b = Function(FSpaces[2])
    b.vector()[:] = bb
    diffb = b.vector().array()

    r = Function(FSpaces[3])
    print r.vector().array().shape
    print rr.shape
    r.vector()[:] = rr


    if (NormType == '2'):
        epsu = splin.norm(diffu)/sqrt(dim[0])
        epsp = splin.norm(p.vector().array())/sqrt(dim[1])
        epsb = splin.norm(diffb)/sqrt(dim[2])
        epsr = splin.norm(r.vector().array())/sqrt(dim[3])
    elif (NormType == 'inf'):
        epsu = splin.norm(diffu, ord=np.Inf)
        epsp = splin.norm(p.vector().array(),ord=np.inf)
        epsb = splin.norm(diffb, ord=np.Inf)
        epsr = splin.norm(r.vector().array(),ord=np.inf)
    else:
        print "NormType must be 2 or inf"
        quit()


    RHS = IO.vecToArray(U)
    u.vector()[:] = uu + RHS[0:dim[0]]
    p.vector()[:] = p.vector().array() + U.array[dim[0]:dim[0]+dim[1]]
    b.vector()[:] = bb + RHS[dim[0]+dim[1]:dim[0]+dim[1]+dim[2]]
    r.vector()[:] = rr + RHS[dim[0]+dim[1]+dim[2]:]




    print 'iter=%d:\n u-norm=%g   p-norm=%g  \n b-norm=%g   r-norm=%g' % (iter, epsu,epsp,epsb,epsr), '\n\n\n'
    return u,p,b,r,epsu+epsp+epsb+epsr

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

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()

def PicardTolerance(x,u_k,b_k,FSpaces,dim,NormType,iter):
    X = IO.vecToArray(x)
    uu = X[0:dim[0]]
    bb = X[dim[0]+dim[1]:dim[0]+dim[1]+dim[2]]

    u = Function(FSpaces[0])
    u.vector()[:] = u.vector()[:] + uu
    diffu = u.vector().array() - u_k.vector().array()

    b = Function(FSpaces[2])
    b.vector()[:] = b.vector()[:] + bb
    diffb = b.vector().array() - b_k.vector().array()
    if (NormType == '2'):
        epsu = splin.norm(diffu)/sqrt(dim[0])
        epsb = splin.norm(diffb)/sqrt(dim[0])
    elif (NormType == 'inf'):
        epsu = splin.norm(diffu, ord=np.Inf)
        epsb = splin.norm(diffb, ord=np.Inf)
    else:
        print "NormType must be 2 or inf"
        quit()

    print 'iter=%d: u-norm=%g   b-norm=%g ' % (iter, epsu,epsb)
    u_k.assign(u)
    b_k.assign(b)


    return u_k,b_k,epsu,epsb

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)

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

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

    maxiter = 100        # max no of iterations allowed

    while eps > tol and iter < maxiter:
            iter += 1
            x = Function(W)

            uu = Function(W)
            tic()
            AA, bb = assemble_system(a, L1, bcs)
            A = as_backend_type(AA).mat()
            print toc()


            b = bb.array()
            zeros = 0*b
            bb = IO.arrayToVec(b)
            x = IO.arrayToVec(zeros)

            tic()
            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,A)
            ksp.setTolerances(1e-6)
            ksp.setType('gmres')
            pc = ksp.getPC()
            pc.setType(pc.Type.FIELDSPLIT)
            fields = [ ("field1", u_is), ("field2", p_is)]
            pc.setFieldSplitIS(*fields)
            pc.setFieldSplitType(0)

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

    f = -MU*Laplacian+gradPres



    n = FacetNormal(mesh)
    h = CellSize(mesh)
    h_avg =avg(h)
    d = 0
    u_k,p_k = common.Stokes(V,Q,u0,f,[1,1,MU])

    uOld = np.concatenate((u_k.vector().array(),p_k.vector().array()), axis=0)
    r = IO.arrayToVec(uOld)

    a11 = MU*inner(grad(v), grad(u))*dx + inner((grad(u)*u_k),v)*dx + (1/2)*div(u_k)*inner(u,v)*dx- (1/2)*inner(u_k,n)*inner(u,v)*ds
    a12 = div(v)*p*dx
    a21 = div(u)*q*dx
    L1  = inner(v, f)*dx
    a = a11-a12-a21


    r11 = MU*inner(grad(v), grad(u_k))*dx + inner((grad(u_k)*u_k),v)*dx + (1/2)*div(u_k)*inner(u_k,v)*dx- (1/2)*inner(u_k,n)*inner(u_k,v)*ds
    r12 = div(v)*p_k*dx
    r21 = div(u_k)*q*dx
    RHSform = r11-r12-r21


    p11 = inner(u,v)*dx

def foo():
    m = 6


    errL2u =np.zeros((m-1,1))
    errH1u =np.zeros((m-1,1))
    errL2p =np.zeros((m-1,1))
    errL2b =np.zeros((m-1,1))
    errCurlb =np.zeros((m-1,1))
    errL2r =np.zeros((m-1,1))
    errH1r =np.zeros((m-1,1))



    l2uorder =  np.zeros((m-1,1))
    H1uorder =np.zeros((m-1,1))
    l2porder =  np.zeros((m-1,1))
    l2border =  np.zeros((m-1,1))
    Curlborder =np.zeros((m-1,1))
    l2rorder =  np.zeros((m-1,1))
    H1rorder = np.zeros((m-1,1))

    NN = np.zeros((m-1,1))
    DoF = np.zeros((m-1,1))
    Velocitydim = np.zeros((m-1,1))
    Magneticdim = np.zeros((m-1,1))
    Pressuredim = np.zeros((m-1,1))
    Lagrangedim = np.zeros((m-1,1))
    Wdim = np.zeros((m-1,1))
    iterations = np.zeros((m-1,1))
    SolTime = np.zeros((m-1,1))
    udiv = np.zeros((m-1,1))
    MU = np.zeros((m-1,1))
    level = np.zeros((m-1,1))
    NSave = np.zeros((m-1,1))
    Mave = np.zeros((m-1,1))
    TotalTime = np.zeros((m-1,1))
    
    nn = 2

    dim = 2
    ShowResultPlots = 'yes'
    split = 'Linear'

    MU[0]= 1e0
    for xx in xrange(1,m):
        print xx
        level[xx-1] = xx+ 0
        nn = 2**(level[xx-1])



        # Create mesh and define function space
        nn = int(nn)
        NN[xx-1] = nn/2
        # parameters["form_compiler"]["quadrature_degree"] = 6
        # parameters = CP.ParameterSetup()
        mesh = UnitCubeMesh(nn,nn,nn)

        order = 2
        parameters['reorder_dofs_serial'] = False
        Velocity = VectorFunctionSpace(mesh, "CG", order)
        Pressure = FunctionSpace(mesh, "CG", order-1)
        Magnetic = FunctionSpace(mesh, "N1curl", order)
        Lagrange = FunctionSpace(mesh, "CG", order)
        W = MixedFunctionSpace([Velocity, Pressure, Magnetic,Lagrange])
        # W = Velocity*Pressure*Magnetic*Lagrange
        Velocitydim[xx-1] = Velocity.dim()
        Pressuredim[xx-1] = Pressure.dim()
        Magneticdim[xx-1] = Magnetic.dim()
        Lagrangedim[xx-1] = Lagrange.dim()
        Wdim[xx-1] = W.dim()
        print "\n\nW:  ",Wdim[xx-1],"Velocity:  ",Velocitydim[xx-1],"Pressure:  ",Pressuredim[xx-1],"Magnetic:  ",Magneticdim[xx-1],"Lagrange:  ",Lagrangedim[xx-1],"\n\n"
        dim = [Velocity.dim(), Pressure.dim(), Magnetic.dim(), Lagrange.dim()]


        def boundary(x, on_boundary):
            return on_boundary

        u0, p0,b0, r0, Laplacian, Advection, gradPres,CurlCurl, gradR, NS_Couple, M_Couple = ExactSol.MHD3D(4,1)


        bcu = DirichletBC(Velocity,u0, boundary)
        bcb = DirichletBC(Magnetic,b0, boundary)
        bcr = DirichletBC(Lagrange,r0, boundary)

        # bc = [u0,p0,b0,r0]
        bcs = [bcu,bcb,bcr]
        FSpaces = [Velocity,Pressure,Magnetic,Lagrange]


        (u, b, p, r) = TrialFunctions(W)
        (v, c, q, s) = TestFunctions(W)
        kappa = 1.0
        Mu_m =10.0
        MU = 1.0/1
        IterType = 'Full'
        Split = "No"
        Saddle = "No"
        Stokes = "No"
#.........这里部分代码省略.........

        del A,P


    if (Solving == 'Iterative' or Solving == 'Direct'):
        
        ue = u0
        pe = p0

        Ve = FunctionSpace(mesh,"N1curl",6)
        u = interpolate(ue,Ve)
        Qe = FunctionSpace(mesh,"CG",6)
        p = interpolate(pe,Qe)



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

        pp = X[V.dim():]
        pa = Function(Q)

        pa.vector()[:] = pp

        parameters["form_compiler"]["quadrature_degree"] = 16

        ErrorB = Function(V)
        ErrorR = Function(Q)

        ErrorB = u-ua

    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)

def solve(A, b, u, params, Fspace, SolveType, IterType, OuterTol, InnerTol, HiptmairMatrices, Hiptmairtol, KSPlinearfluids, Fp, kspF):

    if SolveType == "Direct":
        ksp = PETSc.KSP()
        ksp.create(comm=PETSc.COMM_WORLD)
        pc = ksp.getPC()
        ksp.setType('preonly')
        pc.setType('lu')
        OptDB = PETSc.Options()
        OptDB['pc_factor_mat_solver_package'] = "pastix"
        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
    elif SolveType == "Direct-class":
        ksp = PETSc.KSP()
        ksp.create(comm=PETSc.COMM_WORLD)
        pc = ksp.getPC()
        ksp.setType('gmres')
        pc.setType('none')
        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

    else:

        # u = b.duplicate()
        if IterType == "Full":
            ksp = PETSc.KSP()
            ksp.create(comm=PETSc.COMM_WORLD)
            pc = ksp.getPC()
            ksp.setType('gmres')
            pc.setType('python')

            OptDB = PETSc.Options()
            OptDB['ksp_gmres_restart'] = 100
            # FSpace = [Velocity,Magnetic,Pressure,Lagrange]
            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)
#.........这里部分代码省略.........

    # <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()

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