How to use the petsc4py.PETSc.Vec function in petsc4py

# --- Interpolate (or map) in parallel the fluid interface loads on the solid interface ---
	#self.MappingMatrix.transpose()
        if FSI_config['MATCHING_MESH'] == 'NO' and (FSI_config['MESH_INTERP_METHOD'] == 'RBF' or FSI_config['MESH_INTERP_METHOD'] == 'TPS'):
          if self.have_MPI == True:
            gamma_array_LoadX = PETSc.Vec().create(self.comm)
            gamma_array_LoadY = PETSc.Vec().create(self.comm)
            gamma_array_LoadZ = PETSc.Vec().create(self.comm)
            gamma_array_LoadX.setType('mpi')
            gamma_array_LoadY.setType('mpi')
            gamma_array_LoadZ.setType('mpi')
            KSP_solver = PETSc.KSP().create(self.comm)
          else:
            gamma_array_LoadX = PETSc.Vec().create()
            gamma_array_LoadY = PETSc.Vec().create()
            gamma_array_LoadZ = PETSc.Vec().create()
            gamma_array_LoadX.setType('seq')
            gamma_array_LoadY.setType('seq')
            gamma_array_LoadZ.setType('seq')
            KSP_solver = PETSc.KSP().create()
          gamma_array_LoadX.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF)
          gamma_array_LoadY.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF)
          gamma_array_LoadZ.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF)
          gamma_array_LoadX.set(0.0)
          gamma_array_LoadY.set(0.0)
          gamma_array_LoadZ.set(0.0)
          KSP_solver.setType('fgmres')
          KSP_solver.getPC().setType('jacobi')
          KSP_solver.setOperators(self.MappingMatrixA_T)
          KSP_solver.setFromOptions()
          self.MappingMatrixB_T.mult(self.fluidLoads_array_X, gamma_array_LoadX)
          self.MappingMatrixB_T.mult(self.fluidLoads_array_Y, gamma_array_LoadY)

sizes = []
        for isurf in xrange(self.nSurf):
            sizes.append([self.surfs[isurf].Nu,self.surfs[isurf].Nv])
        # end for

        Npts, g_index,l_index = self.calcGlobalNumbering(sizes)
        
        self._initJacobian(Npts,Nctl)
        if not self.NO_PRINT:
            print '------------- Fitting Surfaces Globally ------------------'
            print 'Npts:',Npts
            print 'Nctl:',Nctl

        if USE_PETSC:
            pts = PETSc.Vec()
            pts = pts.createSeq(Npts*3)

            X = PETSc.Vec()
            X.createSeq(Nctl*3)

            PETSC_INSERT_MODE = PETSc.InsertMode.ADD_VALUES

        else:
            pts = zeros(Npts*3)
            X = zeros(Nctl*3)
        # end if 

        # Now Fill up the pt list
        
        for ii in xrange(len(g_index)):
            isurf = g_index[ii][0][0]

""" Build the linear algebra operators needed to compute the point gauges.

        The operators are all local since the point gauge measurements are calculated locally.
        """

        for field, field_id in zip(self.fields, self.field_ids):
            m = PETSc.Mat().create(PETSc.COMM_SELF)
            m.setSizes([len(self.measuredQuantities[field]),
                        self.u[field_id].femSpace.dim])
            m.setType('aij')
            m.setUp()
            # matrices are a list in same order as fields
            self.pointGaugeMats.append(m)
            # dofs are a list in same order as fields as well
            dofs = self.u[field_id].dof
            dofsVec = PETSc.Vec().createWithArray(dofs, comm=PETSc.COMM_SELF)
            self.dofsVecs.append(dofsVec)

        for field, field_id, m in zip(self.fields, self.field_ids, self.pointGaugeMats):
            # get the FiniteElementFunction object for this quantity
            femFun = self.u[field_id]
            for quantity_id, quantity in enumerate(self.measuredQuantities[field]):
                location, node = quantity
                logEvent("Gauge for: %s at %e %e %e is on local operator row %d" % (field, location[0], location[1],
                                                                      location[2], quantity_id), 3)
                self.buildQuantityRow(m, femFun, quantity_id, quantity)
            pointGaugesVec = PETSc.Vec().create(comm=PETSc.COMM_SELF)
            pointGaugesVec.setSizes(len(self.measuredQuantities[field]))
            pointGaugesVec.setUp()
            self.pointGaugeVecs.append(pointGaugesVec)

        for m in self.pointGaugeMats:

self.coef = []
        # Now Fill up the self.coef list:
        for ii in xrange(len(self.g_index)):
            isurf = self.g_index[ii][0][0]
            i = self.g_index[ii][0][1]
            j = self.g_index[ii][0][2]
            self.coef.append( self.surfs[isurf].coef[i,j])
        # end for
            
        # Finally turn self.coef into a complex array
        self.coef = array(self.coef,'D')

        # Create a PETSc vector of the global coefficients
        if USE_PETSC:
            self.petsc_coef = PETSc.Vec()
            self.petsc_coef.createSeq(3*self.Ncoef)
            self.petsc_coef[:] = self.coef.flatten().astype('d')
            self.petsc_coef.assemble()
        # end
        return

self._initialize_b_x()

        # PETSc
        self.ksp.setOperators(self.petsc_A)
        self.ksp.setFromOptions()

        self.ksp.solve(self.petsc_b, self.petsc_x)

        # Convert solution vector from PETSc.Vec instance
        # to a numpy array
        self.solution = PETSc.Vec.getArray(self.petsc_x)

        # Destroy petsc solver, coefficients matrix, rhs, and solution vectors
        PETSc.KSP.destroy(self.ksp)
        PETSc.Mat.destroy(self.petsc_A)
        PETSc.Vec.destroy(self.petsc_b)
        PETSc.Vec.destroy(self.petsc_x)

        return(self.solution)

x_values = np.array(x_values)

    x = bb.copy()
    x.zeroEntries()
    x.setValues(rows, x_values)

    # Apply to monolithic
    len(rows) and AA.zeroRowsColumns(rows, diag=diag_val, x=x, b=bb)

    blocks = []
    for first, last in zip(offsets[:-1], offsets[1:]):
        blocks.append(PETSc.IS().createStride(last-first, first, 1))

    # Reasamble
    comm = mpi_comm_world()
    b = [PETSc.Vec().createWithArray(bb.getValues(block), comm=comm) for block in blocks]
    for bi in b:
        bi.assemblyBegin()
        bi.assemblyEnd()
    b = block_vec(map(PETScVector, b))

    # PETSc 3.7.x
    try:
        AA.getSubMatrix
        A = [[PETScMatrix(AA.getSubMatrix(block_row, block_col)) for block_col in blocks]
             for block_row in blocks]
    # NOTE: 3.8+ does not ahve getSubMatrix, there is getLocalSubMatrix
    # but cannot get that to work - petsc error 73. So for now everything
    # is done with scipy    
    except AttributeError:
        AA = csr_matrix(AA.getValuesCSR()[::-1], shape=AA.size)

def _vector_load(directory, filename, mpi_comm):
        if _file_exists(directory, filename + ".dat", mpi_comm):
            viewer = PETSc.Viewer().createBinary(os.path.join(str(directory), filename + ".dat"), "r", mpi_comm)
            return PETSc.Vec().load(viewer)
        else:
            raise OSError

wing_box.writeTecplot('./c172_layout.dat')

pre_con = 'mg'

if pre_con == 'asm':
    nmesh = 1
    structure, tacs = wing_box.finalize( nlevels = nmesh )
    mesh = structure.getMesh()

    tacs_assembler = TACS.TACS_PETScAssembler_Real( tacs[0] )

    mat   = PETSc.Mat()
    pcmat = PETSc.Mat()
    res = PETSc.Vec()
    vec = PETSc.Vec()

    tacs_assembler.createMat( mat )
    tacs_assembler.createMat( pcmat )
    tacs_assembler.createVec( vec )
    tacs_assembler.createVec( res )

    loadCase =0
    res.zeroEntries()        
    tacs_assembler.addDirectLoads( loadCase, res )
    res.scale(-1.0)
    tacs_assembler.assembleRes( loadCase, res )

    paramSet = TACSAnalysis.ParameterSet()
    
    if MPI != None and MPI.COMM_WORLD.size > 1:
        # Global preconditioner options

ub.array = 1*np.ones(n)
    tao = PETSc.TAO().create(tao_comm)
    tao.setType(params['localopt_method'])

    if params['localopt_method'] == 'pounders':
        f = PETSc.Vec().create(tao_comm)
        f.setSizes(m)
        f.setFromOptions()

        delta_0 = params['delta_0_mult']*np.min([np.min(ub.array-x.array), np.min(x.array-lb.array)])

        PETSc.Options().setValue('-tao_pounders_delta',str(delta_0))
        PETSc.Options().setValue('-pounders_subsolver_tao_type','bqpip')
        tao.setSeparableObjective(lambda tao, x, f: pounders_obj_func(tao, x, f, Run_H), f)
    elif params['localopt_method'] == 'blmvm':
        g = PETSc.Vec().create(tao_comm)
        g.setSizes(n)
        g.setFromOptions()
        tao.setObjectiveGradient(lambda tao, x, g: blmvm_obj_func(tao, x, g, Run_H))

    # Set everything for tao before solving
    PETSc.Options().setValue('-tao_max_funcs',str(total_pts_in_run+1))
    tao.setFromOptions()
    tao.setVariableBounds((lb,ub))
    # tao.setObjectiveTolerances(fatol=params['fatol'], frtol=params['frtol'])
    # tao.setGradientTolerances(grtol=params['grtol'], gatol=params['gatol'])
    tao.setTolerances(grtol=params['grtol'], gatol=params['gatol'])
    tao.setInitial(x)

    tao.solve(x)

    x_opt = tao.getSolution().getArray()

col_lgmap = PETSc.LGMap().create(map(int, col_lgmap), comm=comm)
        mat.setLGMap(row_lgmap, col_lgmap)

        mat.assemblyBegin()
        for row, columns, values in entries[index]:
            mat.setValues([row], columns, values, PETSc.InsertMode.INSERT_VALUES)
        mat.assemblyEnd()

        matT = PETSc.Mat()
        mat.transpose(matT)
        matT.setLGMap(col_lgmap, row_lgmap)

        mats.append(mat)
        matsT.append(matT)

    rhs = PETSc.Vec().createWithArray(np.zeros(len(relations)))

    return map(lambda x: PETScMatrix(x), mats),\
           map(lambda x: PETScMatrix(x), matsT),\
           PETScVector(rhs)