xref: /petsc/src/dm/tutorials/ex25.c (revision 40badf4fbc550ac1f60bd080eaff6de6d55b946d)
1 
2 static char help[] = "Takes a patch of a large DMDA vector to one process.\n\n";
3 
4 #include <petscdm.h>
5 #include <petscdmda.h>
6 #include <petscdmpatch.h>
7 #include <petscsf.h>
8 
9 typedef struct {
10   PetscScalar x,y;
11 } Field;
12 
13 int main(int argc,char **argv)
14 {
15   Vec            xy,sxy;
16   DM             da,sda = NULL;
17   PetscSF        sf;
18   PetscErrorCode ierr;
19   PetscInt       m = 10, n = 10, dof = 2;
20   MatStencil     lower = {0,3,2,0}, upper = {0,7,8,0}; /* These are in the order of the z, y, x, logical coordinates, the fourth entry is ignored */
21   MPI_Comm       comm;
22   PetscMPIInt    rank;
23 
24   ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
25 
26   /* create the large DMDA and set coordinates (which we will copy down to the small DA). */
27   CHKERRQ(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,PETSC_DECIDE,PETSC_DECIDE,dof,1,0,0,&da));
28   CHKERRQ(DMSetFromOptions(da));
29   CHKERRQ(DMSetUp(da));
30   CHKERRQ(DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0));
31   /* Just as a simple example we use the coordinates as the variables in the vectors we wish to examine. */
32   CHKERRQ(DMGetCoordinates(da,&xy));
33   /* The vector entries are displayed in the "natural" ordering on the two dimensional grid; interlaced x and y with with the x variable increasing more rapidly than the y */
34   CHKERRQ(VecView(xy,0));
35 
36   CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
37   if (rank == 0) comm = MPI_COMM_SELF;
38   else comm = MPI_COMM_NULL;
39 
40   CHKERRQ(DMPatchZoom(da,lower,upper,comm,&sda, NULL,&sf));
41   if (rank == 0) {
42     CHKERRQ(DMCreateGlobalVector(sda,&sxy));
43   } else {
44     CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,0,&sxy));
45   }
46   /*  A PetscSF can also be used as a VecScatter context */
47   CHKERRQ(VecScatterBegin(sf,xy,sxy,INSERT_VALUES,SCATTER_FORWARD));
48   CHKERRQ(VecScatterEnd(sf,xy,sxy,INSERT_VALUES,SCATTER_FORWARD));
49   /* Only rank == 0 has the entries of the patch, so run code only at that rank */
50   if (rank == 0) {
51     Field         **vars;
52     DMDALocalInfo info;
53     PetscInt      i,j;
54     PetscScalar   sum = 0;
55 
56     /* The vector entries of the patch are displayed in the "natural" ordering on the two grid; interlaced x and y with with the x variable increasing more rapidly */
57     CHKERRQ(VecView(sxy,PETSC_VIEWER_STDOUT_SELF));
58     /* Compute some trivial statistic of the coordinates */
59     CHKERRQ(DMDAGetLocalInfo(sda,&info));
60     CHKERRQ(DMDAVecGetArray(sda,sxy,&vars));
61     /* Loop over the patch of the entire domain */
62     for (j=info.ys; j<info.ys+info.ym; j++) {
63       for (i=info.xs; i<info.xs+info.xm; i++) {
64         sum += vars[j][i].x;
65       }
66     }
67     CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"The sum of the x coordinates is %g\n",(double)PetscRealPart(sum)));
68     CHKERRQ(DMDAVecRestoreArray(sda,sxy,&vars));
69   }
70 
71   CHKERRQ(VecDestroy(&sxy));
72   CHKERRQ(PetscSFDestroy(&sf));
73   CHKERRQ(DMDestroy(&sda));
74   CHKERRQ(DMDestroy(&da));
75   ierr = PetscFinalize();
76   return ierr;
77 }
78 
79 /*TEST
80 
81    test:
82 
83    test:
84      nsize: 4
85      suffix: 2
86 
87 TEST*/
88