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>
12: #define ADD_STAGES
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];
28: #if defined(PETSC_USE_LOG) && defined(ADD_STAGES)
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 */
40: #if defined(PETSC_USE_LOG) && defined(ADD_STAGES)
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) {
187: /* radius */
188: PetscReal radius = PetscSqrtScalar((x-.5+h/2)*(x-.5+h/2)+(y-.5+h/2)*(y-.5+h/2)+
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 ) {
220: MatSetValuesBlocked(Amat,8,idx,8,idx,(const PetscScalar*)DD,ADD_VALUES);
221:
222: VecSetValuesBlocked(bb,8,idx,(const PetscScalar*)vv,ADD_VALUES);
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];
228: MatSetValuesBlocked(Amat,8,idx,8,idx,(const PetscScalar*)DD,ADD_VALUES);
229:
230: VecSetValuesBlocked(bb,8,idx,(const PetscScalar*)v2,ADD_VALUES);
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 );
255: #if defined(PETSC_USE_LOG) && defined(ADD_STAGES)
256: PetscLogStagePush(stage[0]);
257: #endif
259: /* PC setup basically */
260: KSPSetUp( ksp );
262: #if defined(PETSC_USE_LOG) && defined(ADD_STAGES)
263: PetscLogStagePop();
264: PetscLogStagePush(stage[1]);
265: #endif
267: /* 1st solve */
268: KSPSolve( ksp, bb, xx );
270: #if defined(PETSC_USE_LOG) && defined(ADD_STAGES)
271: PetscLogStagePop();
272: #endif
274: /* 2nd solve */
275: /* #define TwoSolve */
276: #if defined(TwoSolve)
277: {
278: PetscReal emax, emin;
279: #if defined(PETSC_USE_LOG) && defined(ADD_STAGES)
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 );
287: #if defined(PETSC_USE_LOG) && defined(ADD_STAGES)
288: PetscLogStagePop();
289: PetscLogStagePush(stage[3]);
290: #endif
291: KSPSolve( ksp, bb, xx );
292: KSPComputeExtremeSingularValues( ksp, &emax, &emin );
293:
294: #if defined(PETSC_USE_LOG) && defined(ADD_STAGES)
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:
304: #if defined(PETSC_USE_LOG) && defined(ADD_STAGES)
305: PetscLogStagePop();
306: PetscLogStagePush(stage[5]);
307: #endif
308:
309: KSPSolve( ksp, bb, xx );
310:
311: #if defined(PETSC_USE_LOG) && defined(ADD_STAGES)
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: }