Actual source code: ex56.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  1: static char help[] = "3D, bi-linear quadrilateral (Q1), displacement finite element formulation\n\
  2: of linear elasticity.  E=1.0, nu=0.25.\n\
  3: Unit square domain with Dirichelet boundary condition on the y=0 side only.\n\
  4: Load of 1.0 in x + 2y direction on all nodes (not a true uniform load).\n\
  5:   -ne <size>      : number of (square) quadrilateral elements in each dimension\n\
  6:   -alpha <v>      : scaling of material coeficient in embedded circle\n\n";

  8: #include <petscksp.h>

 10: static PetscBool log_stages = PETSC_TRUE;
 11: static PetscErrorCode MaybeLogStagePush(PetscLogStage stage) { return log_stages ? PetscLogStagePush(stage) : 0; }
 12: static PetscErrorCode MaybeLogStagePop() { return log_stages ? PetscLogStagePop() : 0; }

 16: int main(int argc,char **args)
 17: {
 18:   Mat            Amat;
 20:   PetscInt       m,nn,M,Istart,Iend,i,j,k,ii,jj,kk,ic,ne=4,id;
 21:   PetscReal      x,y,z,h,*coords,soft_alpha=1.e-3;
 22:   PetscBool      two_solves=PETSC_FALSE,test_nonzero_cols=PETSC_FALSE,use_nearnullspace=PETSC_FALSE;
 23:   Vec            xx,bb;
 24:   KSP            ksp;
 25:   MPI_Comm       comm;
 26:   PetscMPIInt    npe,mype;
 27:   PC             pc;
 28:   PetscScalar    DD[24][24],DD2[24][24];
 29:   PetscLogStage  stage[6];
 30:   PetscScalar    DD1[24][24];
 31:   PCType         type;

 33:   PetscInitialize(&argc,&args,(char*)0,help);
 34:   comm = PETSC_COMM_WORLD;
 35:   MPI_Comm_rank(comm, &mype);
 36:   MPI_Comm_size(comm, &npe);

 38:   PetscOptionsBegin(comm,NULL,"3D bilinear Q1 elasticity options","");
 39:   {
 40:     char nestring[256];
 41:     PetscSNPrintf(nestring,sizeof nestring,"number of elements in each direction, ne+1 must be a multiple of %D (sizes^{1/3})",(PetscInt)(PetscPowReal((PetscReal)npe,1./3.) + .5));
 42:     PetscOptionsInt("-ne",nestring,"",ne,&ne,NULL);
 43:     PetscOptionsBool("-log_stages","Log stages of solve separately","",log_stages,&log_stages,NULL);
 44:     PetscOptionsReal("-alpha","material coefficient inside circle","",soft_alpha,&soft_alpha,NULL);
 45:     PetscOptionsBool("-two_solves","solve additional variant of the problem","",two_solves,&two_solves,NULL);
 46:     PetscOptionsBool("-test_nonzero_cols","nonzero test","",test_nonzero_cols,&test_nonzero_cols,NULL);
 47:     PetscOptionsBool("-use_mat_nearnullspace","MatNearNullSpace API test","",use_nearnullspace,&use_nearnullspace,NULL);
 48:   }
 49:   PetscOptionsEnd();

 51:   if (log_stages) {
 52:     PetscLogStageRegister("Setup", &stage[0]);
 53:     PetscLogStageRegister("Solve", &stage[1]);
 54:     PetscLogStageRegister("2nd Setup", &stage[2]);
 55:     PetscLogStageRegister("2nd Solve", &stage[3]);
 56:     PetscLogStageRegister("3rd Setup", &stage[4]);
 57:     PetscLogStageRegister("3rd Solve", &stage[5]);
 58:   } else {
 59:     for (i=0; i<(PetscInt)(sizeof(stage)/sizeof(stage[0])); i++) stage[i] = -1;
 60:   }

 62:   h = 1./ne; nn = ne+1;
 63:   /* ne*ne; number of global elements */
 64:   M = 3*nn*nn*nn; /* global number of equations */
 65:   if (npe==2) {
 66:     if (mype==1) m=0;
 67:     else m = nn*nn*nn;
 68:     npe = 1;
 69:   } else {
 70:     m = nn*nn*nn/npe;
 71:     if (mype==npe-1) m = nn*nn*nn - (npe-1)*m;
 72:   }
 73:   m *= 3; /* number of equations local*/
 74:   /* Setup solver, get PC type and pc */
 75:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 76:   KSPSetType(ksp, KSPCG);
 77:   KSPSetComputeSingularValues(ksp, PETSC_TRUE);
 78:   KSPGetPC(ksp, &pc);
 79:   PCSetType(pc, PCGAMG); /* default */
 80:   KSPSetFromOptions(ksp);
 81:   PCGetType(pc, &type);

 83:   {
 84:     /* configureation */
 85:     const PetscInt NP = (PetscInt)(PetscPowReal((PetscReal)npe,1./3.) + .5);
 86:     const PetscInt ipx = mype%NP, ipy = (mype%(NP*NP))/NP, ipz = mype/(NP*NP);
 87:     const PetscInt Ni0 = ipx*(nn/NP), Nj0 = ipy*(nn/NP), Nk0 = ipz*(nn/NP);
 88:     const PetscInt Ni1 = Ni0 + (m>0 ? (nn/NP) : 0), Nj1 = Nj0 + (nn/NP), Nk1 = Nk0 + (nn/NP);
 89:     const PetscInt NN  = nn/NP, id0 = ipz*nn*nn*NN + ipy*nn*NN*NN + ipx*NN*NN*NN;
 90:     PetscInt       *d_nnz, *o_nnz,osz[4]={0,9,15,19},nbc;
 91:     PetscScalar    vv[24], v2[24];
 92:     if (npe!=NP*NP*NP) SETERRQ1(comm,PETSC_ERR_ARG_WRONG, "npe=%d: npe^{1/3} must be integer",npe);
 93:     if (nn!=NP*(nn/NP)) SETERRQ1(comm,PETSC_ERR_ARG_WRONG, "-ne %d: (ne+1)%(npe^{1/3}) must equal zero",ne);

 95:     /* count nnz */
 96:     PetscMalloc1((m+1), &d_nnz);
 97:     PetscMalloc1((m+1), &o_nnz);
 98:     for (i=Ni0,ic=0; i<Ni1; i++) {
 99:       for (j=Nj0; j<Nj1; j++) {
100:         for (k=Nk0; k<Nk1; k++) {
101:           nbc = 0;
102:           if (i==Ni0 || i==Ni1-1) nbc++;
103:           if (j==Nj0 || j==Nj1-1) nbc++;
104:           if (k==Nk0 || k==Nk1-1) nbc++;
105:           for (jj=0; jj<3; jj++,ic++) {
106:             d_nnz[ic] = 3*(27-osz[nbc]);
107:             o_nnz[ic] = 3*osz[nbc];
108:           }
109:         }
110:       }
111:     }
112:     if (ic != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ic %D does not equal m %D",ic,m);

114:     /* create stiffness matrix */
115:     MatCreate(comm,&Amat);
116:     MatSetSizes(Amat,m,m,M,M);
117:     MatSetBlockSize(Amat,3);
118:     MatSetType(Amat,MATAIJ);
119:     MatSeqAIJSetPreallocation(Amat,0,d_nnz);
120:     MatMPIAIJSetPreallocation(Amat,0,d_nnz,0,o_nnz);

122:     PetscFree(d_nnz);
123:     PetscFree(o_nnz);

125:     MatGetOwnershipRange(Amat,&Istart,&Iend);

127:     if (m != Iend - Istart) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"m %D does not equal Iend %D - Istart %D",m,Iend,Istart);
128:     /* Generate vectors */
129:     VecCreate(comm,&xx);
130:     VecSetSizes(xx,m,M);
131:     VecSetBlockSize(xx,3);
132:     VecSetFromOptions(xx);
133:     VecDuplicate(xx,&bb);
134:     VecSet(bb,.0);
135:     /* generate element matrices */
136:     {
137:       FILE *file;
138:       char fname[] = "data/elem_3d_elast_v_25.txt";
139:       file = fopen(fname, "r");
140:       if (file == 0) {
141:         PetscPrintf(PETSC_COMM_WORLD,"\t%s failed to open input file '%s'\n",__FUNCT__,fname);
142:         for (i=0; i<24; i++) {
143:           for (j=0; j<24; j++) {
144:             if (i==j) DD1[i][j] = 1.0;
145:             else DD1[i][j] = -.25;
146:           }
147:         }
148:       } else {
149:         double dd;
150:         for (i=0; i<24; i++) {
151:           for (j=0; j<24; j++) {
152:             fscanf(file, "%le", &dd);
153:             DD1[i][j] = dd;
154:           }
155:         }
156:       }
157:       fclose(file);
158:       /* BC version of element */
159:       for (i=0; i<24; i++) {
160:         for (j=0; j<24; j++) {
161:           if (i<12 || (j < 12 && !test_nonzero_cols)) {
162:             if (i==j) DD2[i][j] = 0.1*DD1[i][j];
163:             else DD2[i][j] = 0.0;
164:           } else DD2[i][j] = DD1[i][j];
165:         }
166:       }
167:       /* element residual/load vector */
168:       for (i=0; i<24; i++) {
169:         if (i%3==0) vv[i] = h*h;
170:         else if (i%3==1) vv[i] = 2.0*h*h;
171:         else vv[i] = .0;
172:       }
173:       for (i=0; i<24; i++) {
174:         if (i%3==0 && i>=12) v2[i] = h*h;
175:         else if (i%3==1 && i>=12) v2[i] = 2.0*h*h;
176:         else v2[i] = .0;
177:       }
178:     }

180:     PetscMalloc1((m+1), &coords);
181:     coords[m] = -99.0;

183:     /* forms the element stiffness and coordinates */
184:     for (i=Ni0,ic=0,ii=0; i<Ni1; i++,ii++) {
185:       for (j=Nj0,jj=0; j<Nj1; j++,jj++) {
186:         for (k=Nk0,kk=0; k<Nk1; k++,kk++,ic++) {

188:           /* coords */
189:           x = coords[3*ic] = h*(PetscReal)i;
190:           y = coords[3*ic+1] = h*(PetscReal)j;
191:           z = coords[3*ic+2] = h*(PetscReal)k;
192:           /* matrix */
193:           id = id0 + ii + NN*jj + NN*NN*kk;

195:           if (i<ne && j<ne && k<ne) {
196:             /* radius */
197:             PetscReal radius = PetscSqrtReal((x-.5+h/2)*(x-.5+h/2)+(y-.5+h/2)*(y-.5+h/2)+(z-.5+h/2)*(z-.5+h/2));
198:             PetscReal alpha = 1.0;
199:             PetscInt  jx,ix,idx[8] = { id, id+1, id+NN+1, id+NN,
200:                                        id        + NN*NN, id+1    + NN*NN,
201:                                        id+NN+1 + NN*NN, id+NN + NN*NN };

203:             /* correct indices */
204:             if (i==Ni1-1 && Ni1!=nn) {
205:               idx[1] += NN*(NN*NN-1);
206:               idx[2] += NN*(NN*NN-1);
207:               idx[5] += NN*(NN*NN-1);
208:               idx[6] += NN*(NN*NN-1);
209:             }
210:             if (j==Nj1-1 && Nj1!=nn) {
211:               idx[2] += NN*NN*(nn-1);
212:               idx[3] += NN*NN*(nn-1);
213:               idx[6] += NN*NN*(nn-1);
214:               idx[7] += NN*NN*(nn-1);
215:             }
216:             if (k==Nk1-1 && Nk1!=nn) {
217:               idx[4] += NN*(nn*nn-NN*NN);
218:               idx[5] += NN*(nn*nn-NN*NN);
219:               idx[6] += NN*(nn*nn-NN*NN);
220:               idx[7] += NN*(nn*nn-NN*NN);
221:             }

223:             if (radius < 0.25) alpha = soft_alpha;

225:             for (ix=0; ix<24; ix++) {
226:               for (jx=0;jx<24;jx++) DD[ix][jx] = alpha*DD1[ix][jx];
227:             }
228:             if (k>0) {
229:               MatSetValuesBlocked(Amat,8,idx,8,idx,(const PetscScalar*)DD,ADD_VALUES);
230:               VecSetValuesBlocked(bb,8,idx,(const PetscScalar*)vv,ADD_VALUES);
231:             } else {
232:               /* a BC */
233:               for (ix=0;ix<24;ix++) {
234:                 for (jx=0;jx<24;jx++) DD[ix][jx] = alpha*DD2[ix][jx];
235:               }
236:               MatSetValuesBlocked(Amat,8,idx,8,idx,(const PetscScalar*)DD,ADD_VALUES);
237:               VecSetValuesBlocked(bb,8,idx,(const PetscScalar*)v2,ADD_VALUES);
238:             }
239:           }
240:         }
241:       }

243:     }
244:     MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
245:     MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
246:     VecAssemblyBegin(bb);
247:     VecAssemblyEnd(bb);
248:   }

250:   if (!PETSC_TRUE) {
251:     PetscViewer viewer;
252:     PetscViewerASCIIOpen(comm, "Amat.m", &viewer);
253:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
254:     MatView(Amat,viewer);
255:     PetscViewerDestroy(&viewer);
256:   }

258:   /* finish KSP/PC setup */
259:   KSPSetOperators(ksp, Amat, Amat);
260:   if (use_nearnullspace) {
261:     MatNullSpace matnull;
262:     Vec          vec_coords;
263:     PetscScalar  *c;

265:     VecCreate(MPI_COMM_WORLD,&vec_coords);
266:     VecSetBlockSize(vec_coords,3);
267:     VecSetSizes(vec_coords,m,PETSC_DECIDE);
268:     VecSetUp(vec_coords);
269:     VecGetArray(vec_coords,&c);
270:     for (i=0; i<m; i++) c[i] = coords[i]; /* Copy since Scalar type might be Complex */
271:     VecRestoreArray(vec_coords,&c);
272:     MatNullSpaceCreateRigidBody(vec_coords,&matnull);
273:     MatSetNearNullSpace(Amat,matnull);
274:     MatNullSpaceDestroy(&matnull);
275:     VecDestroy(&vec_coords);
276:   } else {
277:     PCSetCoordinates(pc, 3, m/3, coords);
278:   }

280:   MaybeLogStagePush(stage[0]);

282:   /* PC setup basically */
283:   KSPSetUp(ksp);

285:   MaybeLogStagePop();
286:   MaybeLogStagePush(stage[1]);

288:   /* test BCs */
289:   if (test_nonzero_cols) {
290:     VecZeroEntries(xx);
291:     if (mype==0) VecSetValue(xx,0,1.0,INSERT_VALUES);
292:     VecAssemblyBegin(xx);
293:     VecAssemblyEnd(xx);
294:     KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
295:   }

297:   /* 1st solve */
298:   KSPSolve(ksp, bb, xx);

300:   MaybeLogStagePop();

302:   /* 2nd solve */
303:   if (two_solves) {
304:     PetscReal emax, emin;
305:     PetscReal norm,norm2;
306:     Vec       res;

308:     MaybeLogStagePush(stage[2]);
309:     /* PC setup basically */
310:     MatScale(Amat, 100000.0);
311:     KSPSetOperators(ksp, Amat, Amat);
312:     KSPSetUp(ksp);

314:     MaybeLogStagePop();
315:     MaybeLogStagePush(stage[3]);
316:     KSPSolve(ksp, bb, xx);
317:     KSPComputeExtremeSingularValues(ksp, &emax, &emin);

319:     MaybeLogStagePop();
320:     MaybeLogStagePush(stage[4]);

322:     /* 3rd solve */
323:     MatScale(Amat, 100000.0);
324:     KSPSetOperators(ksp, Amat, Amat);
325:     KSPSetUp(ksp);

327:     MaybeLogStagePop();
328:     MaybeLogStagePush(stage[5]);

330:     KSPSolve(ksp, bb, xx);

332:     MaybeLogStagePop();


335:     VecNorm(bb, NORM_2, &norm2);

337:     VecDuplicate(xx, &res);
338:     MatMult(Amat, xx, res);
339:     VecAXPY(bb, -1.0, res);
340:     VecDestroy(&res);
341:     VecNorm(bb, NORM_2, &norm);
342:     PetscPrintf(PETSC_COMM_WORLD,"[%d]%s |b-Ax|/|b|=%e, |b|=%e, emax=%e\n",0,__FUNCT__,(double)(norm/norm2),(double)norm2,(double)emax);
343:     /*PetscViewerASCIIOpen(comm, "residual.m", &viewer);
344:      PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
345:      VecView(bb,viewer);
346:      PetscViewerDestroy(&viewer);*/


349:     /* PetscViewerASCIIOpen(comm, "rhs.m", &viewer); */
350:     /* PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB); */
351:     /*  */
352:     /* VecView(bb,viewer); */
353:     /* PetscViewerDestroy(&viewer); */

355:     /* PetscViewerASCIIOpen(comm, "solution.m", &viewer); */
356:     /* PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB); */
357:     /*  */
358:     /* VecView(xx, viewer); */
359:     /* PetscViewerDestroy(&viewer); */
360:   }

362:   /* Free work space */
363:   KSPDestroy(&ksp);
364:   VecDestroy(&xx);
365:   VecDestroy(&bb);
366:   MatDestroy(&Amat);
367:   PetscFree(coords);

369:   PetscFinalize();
370:   return 0;
371: }