Actual source code: ex56.c

petsc-3.3-p7 2013-05-11
```  1: static char help[] =
2: "ex56: 3D, bi-linear quadrilateral (Q1), displacement finite element formulation\n\
3: of plain strain linear elasticity, that uses the GAMG PC.  E=1.0, nu=0.25.\n\
4: Unit square domain with Dirichelet boundary condition on the y=0 side only.\n\
5: Load of 1.0 in x direction on all nodes (not a true uniform load).\n\
6:   -ne <size>      : number of (square) quadrilateral elements in each dimention\n\
7:   -alpha <v>      : scaling of material coeficient in embedded circle\n\n";

9: #include <petscksp.h>
10: #include <assert.h>

13:
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:   Vec            xx,bb;
23:   KSP            ksp;
24:   MPI_Comm       wcomm;
25:   PetscMPIInt    npe,mype;
26:   PC pc;
27:   PetscScalar DD[24][24],DD2[24][24];
29:   PetscLogStage  stage[6];
30: #endif
31:   PetscScalar DD1[24][24];
32:   const PCType type;

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

39:   /* log */
41:   PetscLogStageRegister("Setup", &stage[0]);
42:   PetscLogStageRegister("Solve", &stage[1]);
43:   PetscLogStageRegister("2nd Setup", &stage[2]);
44:   PetscLogStageRegister("2nd Solve", &stage[3]);
45:   PetscLogStageRegister("3rd Setup", &stage[4]);
46:   PetscLogStageRegister("3rd Solve", &stage[5]);
47: #endif
48:
49:   PetscOptionsGetInt(PETSC_NULL,"-ne",&ne,PETSC_NULL);
50:   h = 1./ne; nn = ne+1;
51:   /* ne*ne; number of global elements */
52:   PetscOptionsGetReal(PETSC_NULL,"-alpha",&soft_alpha,PETSC_NULL);
53:   M = 3*nn*nn*nn; /* global number of equations */
54:   if(npe==2) {
55:     if(mype==1) m=0;
56:     else m = nn*nn*nn;
57:     npe = 1;
58:   }
59:   else {
60:     m = nn*nn*nn/npe;
61:     if(mype==npe-1) m = nn*nn*nn - (npe-1)*m;
62:   }
63:   m *= 3; /* number of equations local*/
64:   /* Setup solver, get PC type and pc */
65:   KSPCreate(PETSC_COMM_WORLD,&ksp);
66:   KSPSetType( ksp, KSPCG );
67:   KSPSetComputeSingularValues( ksp, PETSC_TRUE );
68:   KSPGetPC( ksp, &pc );
69:   PCSetType( pc, PCGAMG );  /* default */
70:   KSPSetFromOptions( ksp );
71:   PCGetType( pc, &type );

73:   {
74:     /* configureation */
75:     const PetscInt NP = (PetscInt)(pow((double)npe,1./3.) + .5);
76:     if(npe!=NP*NP*NP)SETERRQ1(wcomm,PETSC_ERR_ARG_WRONG, "npe=%d: npe^{1/3} must be integer",npe);
77:     if(nn!=NP*(nn/NP))SETERRQ1(wcomm,PETSC_ERR_ARG_WRONG, "-ne %d: (ne+1)%(npe^{1/3}) must equal zero",ne);
78:     const PetscInt ipx = mype%NP, ipy = (mype%(NP*NP))/NP, ipz = mype/(NP*NP);
79:     const PetscInt Ni0 = ipx*(nn/NP), Nj0 = ipy*(nn/NP), Nk0 = ipz*(nn/NP);
80:     const PetscInt Ni1 = Ni0 + (m>0 ? (nn/NP) : 0), Nj1 = Nj0 + (nn/NP), Nk1 = Nk0 + (nn/NP);
81:     const PetscInt NN = nn/NP, id0 = ipz*nn*nn*NN + ipy*nn*NN*NN + ipx*NN*NN*NN;
82:     PetscInt *d_nnz, *o_nnz,osz[4]={0,9,15,19},nbc;
83:     PetscScalar vv[24], v2[24];
84:
85:     /* count nnz */
86:     PetscMalloc( (m+1)*sizeof(PetscInt), &d_nnz );
87:     PetscMalloc( (m+1)*sizeof(PetscInt), &o_nnz );
88:     for(i=Ni0,ic=0;i<Ni1;i++){
89:       for(j=Nj0;j<Nj1;j++){
90:         for(k=Nk0;k<Nk1;k++){
91:           nbc = 0;
92:           if(i==Ni0 || i==Ni1-1)nbc++;
93:           if(j==Nj0 || j==Nj1-1)nbc++;
94:           if(k==Nk0 || k==Nk1-1)nbc++;
95:           for(jj=0;jj<3;jj++,ic++){
96:             d_nnz[ic] = 3*(27-osz[nbc]);
97:             o_nnz[ic] = 3*osz[nbc];
98:           }
99:         }
100:       }
101:     }
102:     assert(ic==m);
103:
104:     /* create stiffness matrix */
105:     if( strcmp(type, PCPROMETHEUS) == 0 ){
106:       /* prometheus needs BAIJ */
107:       MatCreateBAIJ(wcomm,3,m,m,M,M,27,PETSC_NULL,19,PETSC_NULL,&Amat);
108:     }
109:     else {
110:       MatCreate(wcomm,&Amat);
111:       MatSetSizes(Amat,m,m,M,M);
112:       MatSetBlockSize(Amat,3);
113:       MatSetType(Amat,MATAIJ);
114:       MatSeqAIJSetPreallocation(Amat,0,d_nnz);
115:       MatMPIAIJSetPreallocation(Amat,0,d_nnz,0,o_nnz);
116:     }
117:     PetscFree( d_nnz );
118:     PetscFree( o_nnz );

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

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

171:     PetscMalloc( (m+1)*sizeof(PetscReal), &coords );
172:     coords[m] = -99.0;

174:     /* forms the element stiffness for the Laplacian and coordinates */
175:     for(i=Ni0,ic=0,ii=0;i<Ni1;i++,ii++){
176:       for(j=Nj0,jj=0;j<Nj1;j++,jj++){
177:         for(k=Nk0,kk=0;k<Nk1;k++,kk++,ic++){

179:           /* coords */
180:           x = coords[3*ic] = h*(PetscReal)i;
181:           y = coords[3*ic+1] = h*(PetscReal)j;
182:           z = coords[3*ic+2] = h*(PetscReal)k;
183:           /* matrix */
184:           id = id0 + ii + NN*jj + NN*NN*kk;
185:
186:           if( i<ne && j<ne && k<ne) {
189:                                                (z-.5+h/2)*(z-.5+h/2));
190:             PetscReal alpha = 1.0;
191:             PetscInt jx,ix,idx[8] = { id, id+1, id+NN+1, id+NN,
192:                                       id        + NN*NN, id+1    + NN*NN,
193:                                       id+NN+1 + NN*NN, id+NN + NN*NN };

195:             /* correct indices */
196:             if(i==Ni1-1 && Ni1!=nn){
197:               idx[1] += NN*(NN*NN-1);
198:               idx[2] += NN*(NN*NN-1);
199:               idx[5] += NN*(NN*NN-1);
200:               idx[6] += NN*(NN*NN-1);
201:             }
202:             if(j==Nj1-1 && Nj1!=nn) {
203:               idx[2] += NN*NN*(nn-1);
204:               idx[3] += NN*NN*(nn-1);
205:               idx[6] += NN*NN*(nn-1);
206:               idx[7] += NN*NN*(nn-1);
207:             }
208:             if(k==Nk1-1 && Nk1!=nn) {
209:               idx[4] += NN*(nn*nn-NN*NN);
210:               idx[5] += NN*(nn*nn-NN*NN);
211:               idx[6] += NN*(nn*nn-NN*NN);
212:               idx[7] += NN*(nn*nn-NN*NN);
213:             }
214:
215:             if( radius < 0.25 ){
216:               alpha = soft_alpha;
217:             }
218:             for(ix=0;ix<24;ix++)for(jx=0;jx<24;jx++) DD[ix][jx] = alpha*DD1[ix][jx];
219:             if( k>0 ) {
221:
223:
224:             }
225:             else {
226:               /* a BC */
227:               for(ix=0;ix<24;ix++)for(jx=0;jx<24;jx++) DD[ix][jx] = alpha*DD2[ix][jx];
229:
231:             }
232:           }
233:         }
234:       }

236:     }
237:     MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
238:     MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
239:     VecAssemblyBegin(bb);
240:     VecAssemblyEnd(bb);
241:   }
242:
243:   if( !PETSC_TRUE ) {
244:     PetscViewer viewer;
245:     PetscViewerASCIIOpen(wcomm, "Amat.m", &viewer);
246:     PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB);
247:     MatView(Amat,viewer);
248:     PetscViewerDestroy( &viewer );
249:   }

251:   /* finish KSP/PC setup */
252:   KSPSetOperators( ksp, Amat, Amat, SAME_NONZERO_PATTERN );
253:   PCSetCoordinates( pc, 3, m/3, coords );

256:   PetscLogStagePush(stage[0]);
257: #endif

259:   /* PC setup basically */
260:   KSPSetUp( ksp );

263:   PetscLogStagePop();
264:   PetscLogStagePush(stage[1]);
265: #endif

267:   /* 1st solve */
268:   KSPSolve( ksp, bb, xx );

271:   PetscLogStagePop();
272: #endif

274:   /* 2nd solve */
275: /* #define TwoSolve */
276: #if defined(TwoSolve)
277:   {
278:     PetscReal emax, emin;
280:     PetscLogStagePush(stage[2]);
281: #endif
282:     /* PC setup basically */
283:     MatScale( Amat, 100000.0 );
284:     KSPSetOperators( ksp, Amat, Amat, SAME_NONZERO_PATTERN );
285:     KSPSetUp( ksp );

288:     PetscLogStagePop();
289:     PetscLogStagePush(stage[3]);
290: #endif
291:     KSPSolve( ksp, bb, xx );
292:     KSPComputeExtremeSingularValues( ksp, &emax, &emin );
293:
295:     PetscLogStagePop();
296:     PetscLogStagePush(stage[4]);
297: #endif
298:
299:     /* 3rd solve */
300:     MatScale( Amat, 100000.0 );
301:     KSPSetOperators( ksp, Amat, Amat, SAME_NONZERO_PATTERN );
302:     KSPSetUp( ksp );
303:
305:     PetscLogStagePop();
306:     PetscLogStagePush(stage[5]);
307: #endif
308:
309:     KSPSolve( ksp, bb, xx );
310:
312:     PetscLogStagePop();
313: #endif
314:
315:     PetscReal norm,norm2;
316:     /* PetscViewer viewer; */
317:     Vec res;
318:
319:     VecNorm( bb, NORM_2, &norm2 );
320:
321:     VecDuplicate( xx, &res );
322:     MatMult( Amat, xx, res );
323:     VecAXPY( bb, -1.0, res );
324:     VecDestroy( &res );
325:     VecNorm( bb, NORM_2, &norm );
326:     PetscPrintf(PETSC_COMM_WORLD,"[%d]%s |b-Ax|/|b|=%e, |b|=%e, emax=%e\n",0,__FUNCT__,norm/norm2,norm2,emax);
327:     /*PetscViewerASCIIOpen(wcomm, "residual.m", &viewer);
328:      PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
329:      VecView(bb,viewer);
330:      PetscViewerDestroy( &viewer );*/
331:
332:
333:     /* PetscViewerASCIIOpen(wcomm, "rhs.m", &viewer);   */
334:     /* PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB ); */
335:     /* CHKERRQ( ierr ); */
336:     /* VecView( bb,viewer );            */
337:     /* PetscViewerDestroy( &viewer );   */
338:
339:     /* PetscViewerASCIIOpen(wcomm, "solution.m", &viewer);   */
340:     /* PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB ); */
341:     /*  */
342:     /* VecView( xx, viewer );  */
343:     /* PetscViewerDestroy( &viewer );  */
344:   }
345: #endif

347:   /* Free work space */
348: #if !defined(foo)
349:   KSPDestroy(&ksp);
350: #endif
351:   VecDestroy(&xx);
352:   VecDestroy(&bb);
353:   MatDestroy(&Amat);
354:   PetscFree( coords );
355:
356:   PetscFinalize();
357:   return 0;
358: }

```