Actual source code: ex56.c
petsc-3.5.4 2015-05-23
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: }