Actual source code: ex8.c
petsc-3.11.3 2019-06-26
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*/