Actual source code: ex5.c

petsc-3.11.4 2019-09-28
Report Typos and Errors
  1: static char help[] = "Test DMStag ghosted boundaries in 1d\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 LEFT    DMSTAG_LEFT
 10: #define RIGHT   DMSTAG_RIGHT
 11: #define ELEMENT DMSTAG_ELEMENT

 13: #define PRESSURE_CONST 2.0

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

 17: int main(int argc,char **argv)
 18: {
 19:   PetscErrorCode  ierr;
 20:   DM              dmSol;
 21:   Vec             sol,solRef,solRefLocal,rhs,rhsLocal;
 22:   Mat             A;
 23:   KSP             ksp;
 24:   PC              pc;
 25:   PetscInt        start,n,e,nExtra;
 26:   PetscInt        iu,ip;
 27:   PetscScalar     **arrSol,**arrRHS;

 29:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 30:   DMStagCreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_GHOSTED,3,1,1,DMSTAG_STENCIL_BOX,1,NULL,&dmSol);
 31:   DMSetFromOptions(dmSol);
 32:   DMSetUp(dmSol);

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

 41:   DMStagGetCorners(dmSol,&start,NULL,NULL,&n,NULL,NULL,&nExtra,NULL,NULL);
 42:   DMStagVecGetArrayDOF(dmSol,rhsLocal,&arrRHS);

 44:   /* Get the correct entries for each of our variables in local element-wise storage */
 45:   DMStagGetLocationSlot(dmSol,LEFT,0,&iu);
 46:   DMStagGetLocationSlot(dmSol,ELEMENT,0,&ip);
 47:   for (e=start; e<start+n+nExtra; ++e) {
 48:     {
 49:       arrSol[e][iu] = 2*PRESSURE_CONST;
 50:       arrRHS[e][iu] = 0.0;
 51:     }
 52:     if (e < start+n) {
 53:       arrSol[e][ip] = PRESSURE_CONST;
 54:       arrRHS[e][ip] = PRESSURE_CONST;
 55:     }
 56:   }
 57:   DMStagVecRestoreArrayDOF(dmSol,rhsLocal,&arrRHS);
 58:   DMLocalToGlobalBegin(dmSol,rhsLocal,INSERT_VALUES,rhs);
 59:   DMLocalToGlobalEnd(dmSol,rhsLocal,INSERT_VALUES,rhs);
 60:   DMStagVecRestoreArrayDOF(dmSol,solRefLocal,&arrSol);
 61:   DMLocalToGlobalBegin(dmSol,solRefLocal,INSERT_VALUES,solRef);
 62:   DMLocalToGlobalEnd(dmSol,solRefLocal,INSERT_VALUES,solRef);
 63:   DMRestoreLocalVector(dmSol,&solRefLocal);
 64:   DMRestoreLocalVector(dmSol,&rhsLocal);

 66:   /* Matrix-free Operator */
 67:   DMSetMatType(dmSol,MATSHELL);
 68:   DMCreateMatrix(dmSol,&A);
 69:   MatShellSetOperation(A,MATOP_MULT,(void(*) (void)) ApplyOperator);

 71:   /* Solve */
 72:   DMCreateGlobalVector(dmSol,&sol);
 73:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 74:   KSPSetOperators(ksp,A,A);
 75:   KSPGetPC(ksp,&pc);
 76:   PCSetType(pc,PCNONE);
 77:   KSPSetFromOptions(ksp);
 78:   KSPSolve(ksp,rhs,sol);

 80:   /* Check Solution */
 81:   {
 82:     Vec       diff;
 83:     PetscReal normsolRef,errAbs,errRel;

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

 98:   KSPDestroy(&ksp);
 99:   VecDestroy(&sol);
100:   VecDestroy(&solRef);
101:   VecDestroy(&rhs);
102:   MatDestroy(&A);
103:   DMDestroy(&dmSol);
104:   PetscFinalize();
105:   return ierr;
106: }

108: PetscErrorCode ApplyOperator(Mat A,Vec in,Vec out)
109: {
110:   PetscErrorCode    ierr;
111:   DM                dm;
112:   Vec               inLocal,outLocal;
113:   PetscScalar       **arrIn;
114:   PetscScalar       **arrOut;
115:   PetscInt          start,n,nExtra,ex,idxP,idxU,startGhost,nGhost;
116:   DMBoundaryType    boundaryType;
117:   PetscBool         isFirst,isLast;

120:   MatGetDM(A,&dm);
121:   DMStagGetBoundaryTypes(dm,&boundaryType,NULL,NULL);
122:   if (boundaryType != DM_BOUNDARY_GHOSTED) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"Ghosted boundaries required");
123:   DMGetLocalVector(dm,&inLocal);
124:   DMGetLocalVector(dm,&outLocal);
125:   DMGlobalToLocalBegin(dm,in,INSERT_VALUES,inLocal);
126:   DMGlobalToLocalEnd(dm,in,INSERT_VALUES,inLocal);
127:   DMStagGetCorners(dm,&start,NULL,NULL,&n,NULL,NULL,&nExtra,NULL,NULL);
128:   DMStagGetGhostCorners(dm,&startGhost,NULL,NULL,&nGhost,NULL,NULL);
129:   DMStagVecGetArrayDOFRead(dm,inLocal,&arrIn);
130:   DMStagVecGetArrayDOF(dm,outLocal,&arrOut);
131:   DMStagGetLocationSlot(dm,LEFT,0,&idxU);
132:   DMStagGetLocationSlot(dm,ELEMENT,0,&idxP);
133:   DMStagGetIsFirstRank(dm,&isFirst,NULL,NULL);
134:   DMStagGetIsLastRank(dm,&isLast,NULL,NULL);

136:   /* Set "pressures" on ghost boundaries by copying neighboring values*/
137:   if (isFirst) {
138:     arrIn[-1][idxP] = arrIn[0][idxP];
139:   }
140:   if (isLast){
141:     arrIn[start + n][idxP] = arrIn[start + n - 1][idxP];
142:   }

144:   /* Apply operator on physical points */
145:   for (ex=start; ex<start + n + nExtra; ++ex) {
146:     if (ex < start + n) { /* Don't compute pressure outside domain */
147:       arrOut[ex][idxP] = arrIn[ex][idxP];
148:     }
149:     arrOut[ex][idxU] = arrIn[ex][idxP] + arrIn[ex-1][idxP] - arrIn[ex][idxU];
150:   }
151:   DMStagVecRestoreArrayDOFRead(dm,inLocal,&arrIn);
152:   DMStagVecRestoreArrayDOF(dm,outLocal,&arrOut);
153:   DMLocalToGlobalBegin(dm,outLocal,INSERT_VALUES,out);
154:   DMLocalToGlobalEnd(dm,outLocal,INSERT_VALUES,out);
155:   DMRestoreLocalVector(dm,&inLocal);
156:   DMRestoreLocalVector(dm,&outLocal);
157:   return(0);
158: }

160: /*TEST

162:    test:
163:       suffix: 1
164:       nsize: 1

166:    test:
167:       suffix: 2
168:       nsize: 2

170:    test:
171:       suffix: 3
172:       nsize: 3
173:       args: -stag_grid_x 19

175:    test:
176:       suffix: 4
177:       nsize: 5
178:       args: -stag_grid_x 21 -stag_stencil_width 2

180: TEST*/