xref: /petsc/src/dm/tutorials/ex25.c (revision 4e278199b78715991f5c71ebbd945c1489263e6c)
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 coordindates (which we will copy down to the small DA). */
27   ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,PETSC_DECIDE,PETSC_DECIDE,dof,1,0,0,&da);CHKERRQ(ierr);
28   ierr = DMSetFromOptions(da);CHKERRQ(ierr);
29   ierr = DMSetUp(da);CHKERRQ(ierr);
30   ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);
31   /* Just as a simple example we use the coordinates as the variables in the vectors we wish to examine. */
32   ierr = DMGetCoordinates(da,&xy);CHKERRQ(ierr);
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   ierr = VecView(xy,0);CHKERRQ(ierr);
35 
36   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr);
37   if (!rank) comm = MPI_COMM_SELF;
38   else comm = MPI_COMM_NULL;
39 
40   ierr = DMPatchZoom(da,lower,upper,comm,&sda, NULL,&sf);CHKERRQ(ierr);
41   if (!rank) {
42     ierr = DMCreateGlobalVector(sda,&sxy);CHKERRQ(ierr);
43   } else {
44     ierr = VecCreateSeq(PETSC_COMM_SELF,0,&sxy);
45   }
46   /*  A PetscSF can also be used as a VecScatter context */
47   ierr = VecScatterBegin(sf,xy,sxy,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
48   ierr = VecScatterEnd(sf,xy,sxy,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
49   /* Only rank == 0 has the entries of the patch, so run code only at that rank */
50   if (!rank) {
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     ierr = VecView(sxy,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
58     /* Compute some trivial statistic of the coordinates */
59     ierr = DMDAGetLocalInfo(sda,&info);CHKERRQ(ierr);
60     ierr = DMDAVecGetArray(sda,sxy,&vars);CHKERRQ(ierr);
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     ierr = PetscPrintf(PETSC_COMM_SELF,"The sum of the x coordinates is %g\n",(double)PetscRealPart(sum));CHKERRQ(ierr);
68     ierr = DMDAVecRestoreArray(sda,sxy,&vars);CHKERRQ(ierr);
69   }
70 
71   ierr = VecDestroy(&sxy);CHKERRQ(ierr);
72   ierr = PetscSFDestroy(&sf);CHKERRQ(ierr);
73   ierr = DMDestroy(&sda);CHKERRQ(ierr);
74   ierr = DMDestroy(&da);CHKERRQ(ierr);
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