Actual source code: ex11.c

petsc-3.11.4 2019-09-28
Report Typos and Errors
  1: static char help[] = "Test DMStag ghosted boundaries in 2d\n\n";
  2: /* This solves a very contrived problem - the "pressure" terms are set to a constant function
  3:    and the "velocity" terms are just the sum of neighboring values of these, hence twice the
  4:    constant */
  5:  #include <petscdm.h>
  6:  #include <petscksp.h>
  7:  #include <petscdmstag.h>

  9: #define PRESSURE_CONST 2.0

 11: PetscErrorCode ApplyOperator(Mat,Vec,Vec);

 13: int main(int argc,char **argv)
 14: {
 15:   PetscErrorCode  ierr;
 16:   DM              dmSol;
 17:   Vec             sol,solRef,solRefLocal,rhs,rhsLocal;
 18:   Mat             A;
 19:   KSP             ksp;
 20:   PC              pc;
 21:   PetscInt        startx,starty,nx,ny,ex,ey,nExtrax,nExtray;
 22:   PetscInt        iux,iuy,ip;
 23:   PetscInt        dof0,dof1,dof2;
 24:   PetscScalar     ***arrSol,***arrRHS;

 26:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 27:   /* Note: these defaults are chosen to suit the problem. We allow adjusting
 28:      them to check that things still work when you add ununsed extra dof */
 29:   dof0 = 0;
 30:   dof1 = 1;
 31:   dof2 = 1;
 32:   PetscOptionsGetInt(NULL,NULL,"-dof2",&dof2,NULL);
 33:   DMStagCreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_GHOSTED,DM_BOUNDARY_GHOSTED,3,3,PETSC_DECIDE,PETSC_DECIDE,dof0,dof1,dof2,DMSTAG_STENCIL_BOX,1,NULL,NULL,&dmSol);
 34:   DMSetFromOptions(dmSol);
 35:   DMSetUp(dmSol);

 37:   /* Compute reference solution on the grid, using direct array access */
 38:   DMCreateGlobalVector(dmSol,&rhs);
 39:   DMCreateGlobalVector(dmSol,&solRef);
 40:   DMGetLocalVector(dmSol,&solRefLocal);
 41:   DMGetLocalVector(dmSol,&rhsLocal);
 42:   DMStagVecGetArrayDOF(dmSol,solRefLocal,&arrSol);

 44:   DMStagGetCorners(dmSol,&startx,&starty,NULL,&nx,&ny,NULL,&nExtrax,&nExtray,NULL);
 45:   DMStagVecGetArrayDOF(dmSol,rhsLocal,&arrRHS);

 47:   /* Get the correct entries for each of our variables in local element-wise storage */
 48:   DMStagGetLocationSlot(dmSol,DMSTAG_LEFT,0,&iux);
 49:   DMStagGetLocationSlot(dmSol,DMSTAG_DOWN,0,&iuy);
 50:   DMStagGetLocationSlot(dmSol,DMSTAG_ELEMENT,0,&ip);
 51:     for (ey=starty; ey<starty+ny+nExtray; ++ey) {
 52:       for (ex=startx; ex<startx+nx+nExtrax; ++ex) {
 53:         arrSol[ey][ex][iux] = 2*PRESSURE_CONST;
 54:         arrRHS[ey][ex][iux] = 0.0;
 55:         arrSol[ey][ex][iuy] = 2*PRESSURE_CONST;
 56:         arrRHS[ey][ex][iuy] = 0.0;
 57:         if (ex < startx+nx && ey < starty+ny) {
 58:           arrSol[ey][ex][ip] = PRESSURE_CONST;
 59:           arrRHS[ey][ex][ip] = PRESSURE_CONST;
 60:         }
 61:       }
 62:     }
 63:   DMStagVecRestoreArrayDOF(dmSol,rhsLocal,&arrRHS);
 64:   DMLocalToGlobalBegin(dmSol,rhsLocal,INSERT_VALUES,rhs);
 65:   DMLocalToGlobalEnd(dmSol,rhsLocal,INSERT_VALUES,rhs);
 66:   DMStagVecRestoreArrayDOF(dmSol,solRefLocal,&arrSol);
 67:   DMLocalToGlobalBegin(dmSol,solRefLocal,INSERT_VALUES,solRef);
 68:   DMLocalToGlobalEnd(dmSol,solRefLocal,INSERT_VALUES,solRef);
 69:   DMRestoreLocalVector(dmSol,&solRefLocal);
 70:   DMRestoreLocalVector(dmSol,&rhsLocal);

 72:   /* Matrix-free Operator */
 73:   DMSetMatType(dmSol,MATSHELL);
 74:   DMCreateMatrix(dmSol,&A);
 75:   MatShellSetOperation(A,MATOP_MULT,(void(*) (void)) ApplyOperator);

 77:   /* Solve */
 78:   DMCreateGlobalVector(dmSol,&sol);
 79:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 80:   KSPSetOperators(ksp,A,A);
 81:   KSPGetPC(ksp,&pc);
 82:   PCSetType(pc,PCNONE);
 83:   KSPSetFromOptions(ksp);
 84:   KSPSolve(ksp,rhs,sol);

 86:   /* Check Solution */
 87:   {
 88:     Vec       diff;
 89:     PetscReal normsolRef,errAbs,errRel;

 91:     VecDuplicate(sol,&diff);
 92:     VecCopy(sol,diff);
 93:     VecAXPY(diff,-1.0,solRef);
 94:     VecNorm(diff,NORM_2,&errAbs);
 95:     VecNorm(solRef,NORM_2,&normsolRef);
 96:     errRel = errAbs/normsolRef;
 97:     if (errAbs > 1e14 || errRel > 1e14) {
 98:       PetscPrintf(PetscObjectComm((PetscObject)dmSol),"Error (abs): %g\nError (rel): %g\n",errAbs,errRel);
 99:       PetscPrintf(PetscObjectComm((PetscObject)dmSol),"Non-zero error. Probable failure.\n");
100:     }
101:     VecDestroy(&diff);
102:   }

104:   KSPDestroy(&ksp);
105:   VecDestroy(&sol);
106:   VecDestroy(&solRef);
107:   VecDestroy(&rhs);
108:   MatDestroy(&A);
109:   DMDestroy(&dmSol);
110:   PetscFinalize();
111:   return ierr;
112: }

114: PetscErrorCode ApplyOperator(Mat A,Vec in,Vec out)
115: {
116:   PetscErrorCode    ierr;
117:   DM                dm;
118:   Vec               inLocal,outLocal;
119:   PetscScalar       ***arrIn;
120:   PetscScalar       ***arrOut;
121:   PetscInt          startx,starty,nx,ny,nExtrax,nExtray,ex,ey,idxP,idxUx,idxUy,startGhostx,startGhosty,nGhostx,nGhosty;
122:   PetscBool         isFirstx,isFirsty,isFirstz,isLastx,isLasty,isLastz;

125:   MatGetDM(A,&dm);
126:   DMGetLocalVector(dm,&inLocal);
127:   DMGetLocalVector(dm,&outLocal);
128:   DMGlobalToLocalBegin(dm,in,INSERT_VALUES,inLocal);
129:   DMGlobalToLocalEnd(dm,in,INSERT_VALUES,inLocal);
130:   DMStagGetCorners(dm,&startx,&starty,NULL,&nx,&ny,NULL,&nExtrax,&nExtray,NULL);
131:   DMStagGetGhostCorners(dm,&startGhostx,&startGhosty,NULL,&nGhostx,&nGhosty,NULL);
132:   DMStagVecGetArrayDOFRead(dm,inLocal,&arrIn);
133:   DMStagVecGetArrayDOF(dm,outLocal,&arrOut);
134:   DMStagGetLocationSlot(dm,DMSTAG_LEFT,0,&idxUx);
135:   DMStagGetLocationSlot(dm,DMSTAG_DOWN,0,&idxUy);
136:   DMStagGetLocationSlot(dm,DMSTAG_ELEMENT,0,&idxP);
137:   DMStagGetIsFirstRank(dm,&isFirstx,&isFirsty,&isFirstz);
138:   DMStagGetIsLastRank(dm,&isLastx,&isLasty,&isLastz);

140:   /* Set "pressures" on ghost boundaries by copying neighboring values*/
141:   if (isFirstx) {
142:       for (ey=starty; ey<starty+ny+nExtray; ++ey) {
143:         arrIn[ey][-1][idxP] = arrIn[ey][0][idxP];
144:       }
145:   }
146:   if (isLastx){
147:       for (ey=starty; ey<starty+ny+nExtray; ++ey) {
148:         arrIn[ey][startx + nx][idxP] = arrIn[ey][startx + nx - 1][idxP];
149:       }
150:   }
151:   if (isFirsty) {
152:       for (ex=startx; ex<startx+nx+nExtrax; ++ex) {
153:         arrIn[-1][ex][idxP] = arrIn[0][ex][idxP];
154:       }
155:   }
156:   if  (isLasty) {
157:       for (ex=startx; ex<startx+nx+nExtrax; ++ex) {
158:         arrIn[starty + ny][ex][idxP] = arrIn[starty + ny - 1][ex][idxP];
159:       }
160:   }

162:   /* Apply operator on physical points */
163:   for (ey=starty; ey<starty+ny+nExtray; ++ey) {
164:     for (ex=startx; ex<startx+nx+nExtrax; ++ex) {
165:       if (ex < startx+nx && ey < starty+ny) {/* Don't compute pressure outside domain */
166:         arrOut[ey][ex][idxP] = arrIn[ey][ex][idxP];
167:       }
168:       arrOut[ey][ex][idxUx] = arrIn[ey][ex][idxP] + arrIn[ey][ex-1][idxP] - arrIn[ey][ex][idxUx];
169:       arrOut[ey][ex][idxUy] = arrIn[ey][ex][idxP] + arrIn[ey-1][ex][idxP] - arrIn[ey][ex][idxUy];
170:     }
171:   }
172:   DMStagVecRestoreArrayDOFRead(dm,inLocal,&arrIn);
173:   DMStagVecRestoreArrayDOF(dm,outLocal,&arrOut);
174:   DMLocalToGlobalBegin(dm,outLocal,INSERT_VALUES,out);
175:   DMLocalToGlobalEnd(dm,outLocal,INSERT_VALUES,out);
176:   DMRestoreLocalVector(dm,&inLocal);
177:   DMRestoreLocalVector(dm,&outLocal);
178:   return(0);
179: }

181: /*TEST

183:    test:
184:       suffix: 1
185:       nsize: 1

187:    test:
188:       suffix: 2
189:       nsize: 4

191:    test:
192:       suffix: 3
193:       nsize: 1
194:       args: -stag_stencil_width 2

196:    test:
197:       suffix: 4
198:       nsize: 4
199:       args: -stag_grid_x 4 -stag_grid_y 5 -stag_stencil_width 2

201:    test:
202:       suffix: 5
203:       nsize: 4
204:       args: -stag_dof_0 3 -stag_dof_1 2 -stag_dof_2 4 -stag_stencil_width 3 -stag_grid_x 6 -stag_grid_y 6

206: TEST*/