Actual source code: ex8.c

petsc-3.11.1 2019-04-12
Report Typos and Errors
  1: static char help[] = "Test DMStag ghosted boundaries in 3d\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,startz,nx,ny,nz,ex,ey,ez,nExtrax,nExtray,nExtraz;
 22:   PetscInt        iux,iuy,iuz,ip;
 23:   PetscInt        dof0,dof1,dof2,dof3;
 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 = 0;
 31:   dof2 = 1;
 32:   dof3 = 1;
 33:   DMStagCreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_GHOSTED,DM_BOUNDARY_GHOSTED,DM_BOUNDARY_GHOSTED,3,3,3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof0,dof1,dof2,dof3,DMSTAG_STENCIL_BOX,1,NULL,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,&startz,&nx,&ny,&nz,&nExtrax,&nExtray,&nExtraz);
 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_BACK,0,&iuz);
 51:   DMStagGetLocationSlot(dmSol,DMSTAG_ELEMENT,0,&ip);
 52:   for (ez=startz; ez<startz+nz+nExtraz; ++ez) {
 53:     for (ey=starty; ey<starty+ny+nExtray; ++ey) {
 54:       for (ex=startx; ex<startx+nx+nExtrax; ++ex) {
 55:         arrSol[ez][ey][ex][iux] = 2*PRESSURE_CONST;
 56:         arrRHS[ez][ey][ex][iux] = 0.0;
 57:         arrSol[ez][ey][ex][iuy] = 2*PRESSURE_CONST;
 58:         arrRHS[ez][ey][ex][iuy] = 0.0;
 59:         arrSol[ez][ey][ex][iuz] = 2*PRESSURE_CONST;
 60:         arrRHS[ez][ey][ex][iuz] = 0.0;
 61:         if (ex < startx+nx && ey < starty+ny && ez < startz+nz) {
 62:           arrSol[ez][ey][ex][ip] = PRESSURE_CONST;
 63:           arrRHS[ez][ey][ex][ip] = PRESSURE_CONST;
 64:         }
 65:       }
 66:     }
 67:   }
 68:   DMStagVecRestoreArrayDOF(dmSol,rhsLocal,&arrRHS);
 69:   DMLocalToGlobalBegin(dmSol,rhsLocal,INSERT_VALUES,rhs);
 70:   DMLocalToGlobalEnd(dmSol,rhsLocal,INSERT_VALUES,rhs);
 71:   DMStagVecRestoreArrayDOF(dmSol,solRefLocal,&arrSol);
 72:   DMLocalToGlobalBegin(dmSol,solRefLocal,INSERT_VALUES,solRef);
 73:   DMLocalToGlobalEnd(dmSol,solRefLocal,INSERT_VALUES,solRef);
 74:   DMRestoreLocalVector(dmSol,&solRefLocal);
 75:   DMRestoreLocalVector(dmSol,&rhsLocal);

 77:   /* Matrix-free Operator */
 78:   DMSetMatType(dmSol,MATSHELL);
 79:   DMCreateMatrix(dmSol,&A);
 80:   MatShellSetOperation(A,MATOP_MULT,(void(*) (void)) ApplyOperator);

 82: #if 0
 83:   {
 84:     Mat Aex;
 85:     MatComputeExplicitOperator(A,&Aex);
 86:     MatView(Aex,0);
 87:     MatDestroy(&Aex);
 88:   }
 89: #endif

 91:   /* Solve */
 92:   DMCreateGlobalVector(dmSol,&sol);
 93:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 94:   KSPSetOperators(ksp,A,A);
 95:   KSPGetPC(ksp,&pc);
 96:   PCSetType(pc,PCNONE);
 97:   KSPSetFromOptions(ksp);
 98:   KSPSolve(ksp,rhs,sol);

100:   /* Check Solution */
101:   {
102:     Vec       diff;
103:     PetscReal normsolRef,errAbs,errRel;

105:     VecDuplicate(sol,&diff);
106:     VecCopy(sol,diff);
107:     VecAXPY(diff,-1.0,solRef);
108:     VecNorm(diff,NORM_2,&errAbs);
109:     VecNorm(solRef,NORM_2,&normsolRef);
110:     errRel = errAbs/normsolRef;
111:     if (errAbs > 1e14 || errRel > 1e14) {
112:       PetscPrintf(PetscObjectComm((PetscObject)dmSol),"Error (abs): %g\nError (rel): %g\n",errAbs,errRel);
113:       PetscPrintf(PetscObjectComm((PetscObject)dmSol),"Non-zero error. Probable failure.\n");
114:     }
115:     VecDestroy(&diff);
116:   }

118:   KSPDestroy(&ksp);
119:   VecDestroy(&sol);
120:   VecDestroy(&solRef);
121:   VecDestroy(&rhs);
122:   MatDestroy(&A);
123:   DMDestroy(&dmSol);
124:   PetscFinalize();
125:   return ierr;
126: }

128: PetscErrorCode ApplyOperator(Mat A,Vec in,Vec out)
129: {
130:   PetscErrorCode    ierr;
131:   DM                dm;
132:   Vec               inLocal,outLocal;
133:   PetscScalar       ****arrIn;
134:   PetscScalar       ****arrOut;
135:   PetscInt          startx,starty,startz,nx,ny,nz,nExtrax,nExtray,nExtraz,ex,ey,ez,idxP,idxUx,idxUy,idxUz,startGhostx,startGhosty,startGhostz,nGhostx,nGhosty,nGhostz;
136:   PetscBool         isFirstx,isFirsty,isFirstz,isLastx,isLasty,isLastz;

139:   MatGetDM(A,&dm);
140:   DMGetLocalVector(dm,&inLocal);
141:   DMGetLocalVector(dm,&outLocal);
142:   DMGlobalToLocalBegin(dm,in,INSERT_VALUES,inLocal);
143:   DMGlobalToLocalEnd(dm,in,INSERT_VALUES,inLocal);
144:   DMStagGetCorners(dm,&startx,&starty,&startz,&nx,&ny,&nz,&nExtrax,&nExtray,&nExtraz);
145:   DMStagGetGhostCorners(dm,&startGhostx,&startGhosty,&startGhostz,&nGhostx,&nGhosty,&nGhostz);
146:   DMStagVecGetArrayDOFRead(dm,inLocal,&arrIn);
147:   DMStagVecGetArrayDOF(dm,outLocal,&arrOut);
148:   DMStagGetLocationSlot(dm,DMSTAG_LEFT,0,&idxUx);
149:   DMStagGetLocationSlot(dm,DMSTAG_DOWN,0,&idxUy);
150:   DMStagGetLocationSlot(dm,DMSTAG_BACK,0,&idxUz);
151:   DMStagGetLocationSlot(dm,DMSTAG_ELEMENT,0,&idxP);
152:   DMStagGetIsFirstRank(dm,&isFirstx,&isFirsty,&isFirstz);
153:   DMStagGetIsLastRank(dm,&isLastx,&isLasty,&isLastz);

155:   /* Set "pressures" on ghost boundaries by copying neighboring values*/
156:   if (isFirstx) {
157:     for (ez=startz; ez<startz+nz+nExtraz; ++ez) {
158:       for (ey=starty; ey<starty+ny+nExtray; ++ey) {
159:         arrIn[ez][ey][-1][idxP] = arrIn[ez][ey][0][idxP];
160:       }
161:     }
162:   }
163:   if (isLastx){
164:     for (ez=startz; ez<startz+nz+nExtraz; ++ez) {
165:       for (ey=starty; ey<starty+ny+nExtray; ++ey) {
166:         arrIn[ez][ey][startx + nx][idxP] = arrIn[ez][ey][startx + nx - 1][idxP];
167:       }
168:     }
169:   }
170:   if (isFirsty) {
171:     for (ez=startz; ez<startz+nz+nExtraz; ++ez) {
172:       for (ex=startx; ex<startx+nx+nExtrax; ++ex) {
173:         arrIn[ez][-1][ex][idxP] = arrIn[ez][0][ex][idxP];
174:       }
175:     }
176:   }
177:   if  (isLasty) {
178:     for (ez=startz; ez<startz+nz+nExtraz; ++ez) {
179:       for (ex=startx; ex<startx+nx+nExtrax; ++ex) {
180:         arrIn[ez][starty + ny][ex][idxP] = arrIn[ez][starty + ny - 1][ex][idxP];
181:       }
182:     }
183:   }

185:   if (isFirstz) {
186:     for (ey=starty; ey<starty+ny+nExtray; ++ey) {
187:       for (ex=startx; ex<startx+nx+nExtrax; ++ex) {
188:         arrIn[-1][ey][ex][idxP] = arrIn[0][ey][ex][idxP];
189:       }
190:     }
191:   }
192:   if (isLastz) {
193:     for (ey=starty; ey<starty+ny+nExtray; ++ey) {
194:       for (ex=startx; ex<startx+nx+nExtrax; ++ex) {
195:         arrIn[startz + nz][ey][ex][idxP] = arrIn[startz + nz -1][ey][ex][idxP];
196:       }
197:     }
198:   }

200:   /* Apply operator on physical points */
201:   for (ez=startz; ez<startz+nz+nExtraz; ++ez) {
202:     for (ey=starty; ey<starty+ny+nExtray; ++ey) {
203:       for (ex=startx; ex<startx+nx+nExtrax; ++ex) {
204:         if (ex < startx+nx && ey < starty+ny && ez < startz+nz) {/* Don't compute pressure outside domain */
205:           arrOut[ez][ey][ex][idxP] = arrIn[ez][ey][ex][idxP];
206:         }
207:         arrOut[ez][ey][ex][idxUx] = arrIn[ez][ey][ex][idxP] + arrIn[ez][ey][ex-1][idxP] - arrIn[ez][ey][ex][idxUx];
208:         arrOut[ez][ey][ex][idxUy] = arrIn[ez][ey][ex][idxP] + arrIn[ez][ey-1][ex][idxP] - arrIn[ez][ey][ex][idxUy];
209:         arrOut[ez][ey][ex][idxUz] = arrIn[ez][ey][ex][idxP] + arrIn[ez-1][ey][ex][idxP] - arrIn[ez][ey][ex][idxUz];
210:       }
211:     }
212:   }
213:   DMStagVecRestoreArrayDOFRead(dm,inLocal,&arrIn);
214:   DMStagVecRestoreArrayDOF(dm,outLocal,&arrOut);
215:   DMLocalToGlobalBegin(dm,outLocal,INSERT_VALUES,out);
216:   DMLocalToGlobalEnd(dm,outLocal,INSERT_VALUES,out);
217:   DMRestoreLocalVector(dm,&inLocal);
218:   DMRestoreLocalVector(dm,&outLocal);
219:   return(0);
220: }

222: /*TEST

224:    test:
225:       suffix: 1
226:       nsize: 1

228:    test:
229:       suffix: 2
230:       nsize: 8

232:    test:
233:       suffix: 3
234:       nsize: 1
235:       args: -stag_stencil_width 2

237:    test:
238:       suffix: 4
239:       nsize: 8
240:       args: -stag_grid_x 4 -stag_grid_y 5 -stag_grid_z 4 -stag_stencil_width 2

242:    test:
243:       suffix: 5
244:       nsize: 8
245:       args: -stag_dof_0 3 -stag_dof_1 2 -stag_dof_2 4 -stag_dof_3 2 -stag_stencil_width 3 -stag_grid_x 6 -stag_grid_y 6 -stag_grid_z 6

247: TEST*/