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