Actual source code: ex12.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  1: static char help[] = "Test DMStag 2d star stencil\n\n";
  2:  #include <petscdm.h>
  3:  #include <petscdmstag.h>

  5: int main(int argc,char **argv)
  6: {
  7:   PetscErrorCode  ierr;
  8:   DM              dm;
  9:   Vec             vec,vecLocal1,vecLocal2;
 10:   PetscScalar     *a,***a1,***a2,expected,sum;
 11:   PetscInt        startx,starty,nx,ny,i,j,d,is,js,dof0,dof1,dof2,dofTotal,stencilWidth,ngx,ngy;
 12:   DMBoundaryType  boundaryTypex,boundaryTypey;
 13:   PetscMPIInt     rank;

 15:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 16:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 17:   dof0 = 1;
 18:   dof1 = 1;
 19:   dof2 = 1;
 20:   stencilWidth = 2;
 21:   DMStagCreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,4,4,PETSC_DECIDE,PETSC_DECIDE,dof0,dof1,dof2,DMSTAG_STENCIL_STAR,stencilWidth,NULL,NULL,&dm);
 22:   DMSetFromOptions(dm);
 23:   DMSetUp(dm);
 24:   DMStagGetDOF(dm,&dof0,&dof1,&dof2,NULL);
 25:   dofTotal = dof0 + 2*dof1 + dof2;
 26:   DMStagGetStencilWidth(dm,&stencilWidth);

 28:   DMCreateLocalVector(dm,&vecLocal1);
 29:   VecDuplicate(vecLocal1,&vecLocal2);

 31:   DMCreateGlobalVector(dm,&vec);
 32:   VecSet(vec,1.0);
 33:   VecSet(vecLocal1,0.0);
 34:   DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal1);
 35:   DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal1);

 37:   DMStagGetCorners(dm,&startx,&starty,NULL,&nx,&ny,NULL,NULL,NULL,NULL);
 38:   DMStagVecGetArrayDOFRead(dm,vecLocal1,&a1);
 39:   DMStagVecGetArrayDOF(dm,vecLocal2,&a2);
 40:   for (j=starty; j<starty + ny; ++j) {
 41:     for (i=startx; i<startx + nx; ++i) {
 42:       for (d=0; d<dofTotal; ++d) {
 43:         if (a1[j][i][d] != 1.0) {
 44:           PetscPrintf(PETSC_COMM_SELF,"[%d] Unexpected value %g (expecting %g)\n",rank,a1[j][i][d],1.0);
 45:         }
 46:         a2[j][i][d] = 0.0;
 47:         for (js = -stencilWidth; js <= stencilWidth; ++js) {
 48:           a2[j][i][d] += a1[j+js][i][d];
 49:         }
 50:         for (is = -stencilWidth; is <= stencilWidth; ++is) {
 51:           a2[j][i][d] += a1[j][i+is][d];
 52:         }
 53:         a2[j][i][d] -= a1[j][i][d];
 54:       }
 55:     }
 56:   }
 57:   DMStagVecRestoreArrayDOFRead(dm,vecLocal1,&a1);
 58:   DMStagVecRestoreArrayDOF(dm,vecLocal2,&a2);

 60:   DMLocalToGlobalBegin(dm,vecLocal2,INSERT_VALUES,vec);
 61:   DMLocalToGlobalEnd(dm,vecLocal2,INSERT_VALUES,vec);

 63:   /* For the all-periodic case, some additional checks */
 64:   DMStagGetBoundaryTypes(dm,&boundaryTypex,&boundaryTypey,NULL);
 65:   if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC) {

 67:     DMStagGetGhostCorners(dm,NULL,NULL,NULL,&ngx,&ngy,NULL);
 68:     expected = (ngx*ngy - 4*stencilWidth*stencilWidth)*dofTotal;
 69:     VecSum(vecLocal1,&sum);
 70:     if (sum != expected) {
 71:       PetscPrintf(PETSC_COMM_SELF,"[%d] Unexpected sum of local entries %g (expected %g)\n",rank,sum,expected);
 72:     }

 74:     VecGetArray(vec,&a);
 75:     expected = 1 + 4*stencilWidth;
 76:     for (i=0; i<ny*nx*dofTotal; ++i) {
 77:       if (a[i] != expected) {
 78:         PetscPrintf(PETSC_COMM_SELF,"[%d] Unexpected value %g (expecting %g)\n",rank,a[i],expected);
 79:       }
 80:     }
 81:     VecRestoreArray(vec,&a);
 82:   }

 84:   VecDestroy(&vec);
 85:   VecDestroy(&vecLocal1);
 86:   VecDestroy(&vecLocal2);
 87:   DMDestroy(&dm);
 88:   PetscFinalize();
 89:   return ierr;
 90: }

 92: /*TEST

 94:    test:
 95:       suffix: 1
 96:       nsize: 4
 97:       args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_stencil_width 1

 99:    test:
100:       suffix: 2
101:       nsize: 6
102:       args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_dof_0 2 -stag_grid_x 6

104:    test:
105:       suffix: 3
106:       nsize: 4
107:       args: -stag_dof_0 3 -stag_dof_1 2 -stag_dof_2 4 2 -stag_stencil_width 3 -stag_grid_x 6 -stag_grid_y 6

109:    test:
110:       suffix: 4
111:       nsize: 4
112:       args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_boundary_type_x ghosted

114:    test:
115:       suffix: 5
116:       nsize: 4
117:       args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_boundary_type_y ghosted

119:    test:
120:       suffix: 6
121:       nsize: 4
122:       args: -stag_stencil_width 1 -stag_grid_x 3 -stag_grid_y 2 -stag_boundary_type_x ghosted -stag_boundary_type_y ghosted

124:    test:
125:       suffix: 7
126:       nsize: 4
127:       args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_boundary_type_y ghosted

129:    test:
130:       suffix: 8
131:       nsize: 6
132:       args: -stag_stencil_width 1 -stag_grid_y 2 -stag_grid_x 19 -stag_boundary_type_y ghosted -stag_ranks_x 6
133: TEST*/