Vecview Petsc

[Lalo Scalf]

Vectorsare used to store solutions and right-hand sides for linearsystems of the form Ax=b.This tutorial is an introduction to vector objects (Vec) in PETSc andis based on PETSc Vec tutorial ex1and PETSc Vec tutorial ex3.This tutorial demonstrates how to:

This program creates two PETSc vectors and performs basic operations on them.To utilize the vector class, we need to include petscvec.h whichautomatically includes header files for basic PETSc routines, index sets,and viewers.MPI is set up using PetscInitialize which invokes MPI_Initand reads the command line for PETSc runtime options. In this case, we give thevariable n a default value of 20 with the option of specifying it at runtime.

PETSc supports sequential and parallel vectors. The type can be set in the code ordetermined at runtime.Below, we let the vector type be specified at runtime using command line arguments -vec_type seq or -vec_type mpi.The vector type is then set by VecSetFromOptions().Alternately, sequential vectors can be created with VecCreateSeq() and parallelwith VecCreateMPI().

We set the size of a vector by passing the global size n and PETSC_DECIDE for thelocal size to VecSetSizes(). PETSc will choose a reasonable partitioning trying to put nearly anequal number of elements on each processor.Alternately, we can pass the local size and PETSC_DECIDE and PETSc willcompute the global size. We can duplicate a vector to create another vector,which will have the same format and partitioning as the first.

We can set all entries in a vector to a single value with VecSet().Elements of a vector can be set with the VecSetValues() function.Parallel vectors in PETSc are partitioned into continguous chunksacross processors. We determine which chunks of the vector are locallyowned and set each element using global indexing.A process can set an element in any location, even elements it does not own.Contributions to a vector can either be summed (ADD_VALUES) or replaced (INSERT_VALUES).We assemble a vector in a two-step process. Computation can be done whilemessages are in transition by placing code between these two statements.

PETSc supports many operations on vectors and vector routines such as dot products andnorms. Some of these are demonstrated here.We can visualize vectors and check our results by printing to the consol usingVecView() and PETSc print function PetscPrintf().For large vectors, a PetscViewer can be useful to create Matlab or HDF5 compatible files.

We use the command line to run PETSc programs. Navigate to ex1.cusing the cd command followed by the path to the file. PETSc makesuse of make to build tests and tutorials, so compiling ex1.c isstraightforward and we get an executable called ex1. We run thisprogram in parallel with two MPI ranks. Additionaly, we set the vectorsize and type with command line options.

Vectors (denoted by Vec) are used to store discrete PDE solutions, right-hand sides forlinear systems, etc. Users can create and manipulate entries in vectors directly with a basic, low-level interface orthey can use the PETSc DM objects to connect actions on vectors to the type of discretization and grid that they areworking with. These higher-level interfaces handle much of the details of the interactions with vectors and hence, are preferredin most situations. This chapter is organized as follows:

PETSc provides many ways to create vectors. The most basic, where the user is responsible for managing theparallel distribution of the vector entries, and a variety of higher-level approaches, based on DM, for classes of problems suchas structured grids, staggered grids, unstructured grids, networks, and particles.

which automatically generates the appropriate vector type (sequential orparallel) over all processes in comm. The option -vec_type can be used in conjunction withVecSetFromOptions() to specify the use of a particular type of vector. For example, for NVIDIA GPU CUDA, use cuda.The GPU-based vectors allowone to set values on either the CPU or GPU but do their computations on the GPU.

We emphasize that all processes in comm must call the vectorcreation routines since these routines are collective on allprocesses in the communicator. If you are unfamiliar with MPIcommunicators, see the discussion in Writing PETSc Programs. In addition, if a sequence of creation routines isused, they must be called in the same order for each process in thecommunicator.

For applications running in parallel that involve multi-dimensional structured grids, unstructured grids, networks, etc, it is cumbersomefor users to explicitly manage the needed local and global sizes of the vectors. Hence, PETSc provides a powerful abstractobject called the DM to help manage the vectors and matrices needed for such applications. Parallel vectors can be created easily with

The DM object, see DMDA - Creating vectors for structured grids, DMSTAG - Creating vectors for staggered grids, and DMPlex: Unstructured Grids for more details on DM for structured grids, staggeredstructured grids, and for unstructured grids,manages creating the correctly sized parallel vectors efficiently. One controls the type of vector that DM creates by calling

Each DM type is suitable for a family of problems. The first of these, DMDAare intended for use with logically structured rectangular gridswhen communication of nonlocal data is needed before certain localcomputations can occur. DMDA is designed only forthe case in which data can be thought of as being stored in a standardmultidimensional array; thus, DMDA are not intended forparallelizing unstructured grid problems, etc.

The arguments M and N indicate the global numbers of grid pointsin each direction, while m and n denote the process partition ineach direction; m*n must equal the number of processes in the MPIcommunicator, comm. Instead of specifying the process layout, onemay use PETSC_DECIDE for m and n so that PETSc willselect the partition. The type of periodicity of the arrayis specified by xperiod and yperiod, which can beDM_BOUNDARY_NONE (no periodicity), DM_BOUNDARY_PERIODIC(periodic in that direction), DM_BOUNDARY_TWIST (periodic in thatdirection, but identified in reverse order), DM_BOUNDARY_GHOSTED ,or DM_BOUNDARY_MIRROR. The argument dof indicates the number ofdegrees of freedom at each array point, and s is the stencil width(i.e., the width of the ghost point region). The optional arrays lxand ly may contain the number of nodes along the x and y axis foreach cell, i.e. the dimension of lx is m and the dimension ofly is n; alternately, NULL may be passed in.

Two types of DMDA communication data structures can becreated, as specified by st. Star-type stencils that radiate outwardonly in the coordinate directions are indicated byDMDA_STENCIL_STAR, while box-type stencils are specified byDMDA_STENCIL_BOX. For example, for the two-dimensional case,DMDA_STENCIL_STAR with width 1 corresponds to the standard 5-pointstencil, while DMDA_STENCIL_BOX with width 1 denotes the standard9-point stencil. In both instances, the ghost points are identical, theonly difference being that with star-type stencils, certain ghost pointsare ignored, substantially decreasing the number of messages sent. Notethat the DMDA_STENCIL_STAR stencils can save interprocesscommunication in two and three dimensions.

These DMDA stencils have nothing directly to do with a specific finitedifference stencil one might choose to use for discretization; theyonly ensure that the correct values are in place for the application of auser-defined finite difference stencil (or any other discretizationtechnique).

For structured grids with staggered data (living on elements, faces, edges,and/or vertices), the DMSTAG object is available. It behaves muchlike DMDA.See DMSTAG: Staggered, Structured Grid for discussion of creating vectors with DMSTAG.

To print the vector to the screen, one can use the viewerPETSC_VIEWER_STDOUT_WORLD, which ensures that parallel vectors areprinted correctly to stdout. To display the vector in an X-window,one can use the default X-windows viewer PETSC_VIEWER_DRAW_WORLD, orone can create a viewer with the routine PetscViewerDrawOpen(). Avariety of viewers are discussed further inViewers: Looking at PETSc Objects.

This routine creates an array of pointers to vectors. The two routinesare useful because they allow one to write library code that doesnot depend on the particular format of the vectors being used. Instead,the subroutines can automatically create work vectors based onthe specified existing vector. As discussed inDuplicating Multiple Vectors, the Fortran interface forVecDuplicateVecs() differs slightly.

any number of times on any or all of the processes. The argument ngives the number of components being set in this insertion. The integerarray indices contains the global component indices, andvalues is the array of values to be inserted at those global component index locations. Any process can setany vector components; PETSc ensures that they are automaticallystored in the correct location. Once all of the values have beeninserted with VecSetValues(), one must call

Again, one must call the assembly routines VecAssemblyBegin() andVecAssemblyEnd() after all of the values have been added. Note thataddition and insertion calls to VecSetValues() cannot be mixed.Instead, one must add and insert vector elements in phases, withintervening calls to the assembly routines. This phased assemblyprocedure overcomes the nondeterministic behavior that would occur iftwo different processes generated values for the same location, with oneprocess adding while the other is inserting its value. (In this case, theaddition and insertion actions could be performed in either order, thusresulting in different values at the particular location. Since PETScdoes not allow the simultaneous use of INSERT_VALUES andADD_VALUES this nondeterministic behavior will not occur in PETSc.)

For vectors obtained with DMCreateGlobalVector(), on can use VecSetValuesLocal() to set values intoa global vector but using the local (ghosted) vector indexing of the vector entries. See also Local to global mappingsthat allows one to provide arbitrary local-to-global mapping when not working with a DM.

3a8082e126