Actual source code: ex56.c
1: static char help[] = "3D, tri-quadratic hexahedra (Q1), displacement finite element formulation\n\
2: of linear elasticity. E=1.0, nu=1/3.\n\
3: Unit cube domain with Dirichlet boundary\n\n";
5: #include <petscdmplex.h>
6: #include <petscsnes.h>
7: #include <petscds.h>
8: #include <petscdmforest.h>
10: static PetscReal s_soft_alpha=1.e-3;
11: static PetscReal s_mu=0.4;
12: static PetscReal s_lambda=0.4;
14: static void f0_bd_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
15: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
16: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
17: PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
18: {
19: f0[0] = 1; /* x direction pull */
20: f0[1] = -x[2]; /* add a twist around x-axis */
21: f0[2] = x[1];
22: }
24: static void f1_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
25: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
26: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
27: PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
28: {
29: const PetscInt Ncomp = dim;
30: PetscInt comp, d;
31: for (comp = 0; comp < Ncomp; ++comp) {
32: for (d = 0; d < dim; ++d) {
33: f1[comp*dim+d] = 0.0;
34: }
35: }
36: }
38: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
39: static void f1_u_3d_alpha(PetscInt dim, PetscInt Nf, PetscInt NfAux,
40: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
41: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
42: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
43: {
44: PetscReal trace,mu=s_mu,lambda=s_lambda,rad;
45: PetscInt i,j;
46: for (i=0,rad=0.;i<dim;i++) {
47: PetscReal t=x[i];
48: rad += t*t;
49: }
50: rad = PetscSqrtReal(rad);
51: if (rad>0.25) {
52: mu *= s_soft_alpha;
53: lambda *= s_soft_alpha; /* we could keep the bulk the same like rubberish */
54: }
55: for (i=0,trace=0; i < dim; ++i) {
56: trace += PetscRealPart(u_x[i*dim+i]);
57: }
58: for (i=0; i < dim; ++i) {
59: for (j=0; j < dim; ++j) {
60: f1[i*dim+j] = mu*(u_x[i*dim+j]+u_x[j*dim+i]);
61: }
62: f1[i*dim+i] += lambda*trace;
63: }
64: }
66: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
67: static void f1_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
68: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
69: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
70: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
71: {
72: PetscReal trace,mu=s_mu,lambda=s_lambda;
73: PetscInt i,j;
74: for (i=0,trace=0; i < dim; ++i) {
75: trace += PetscRealPart(u_x[i*dim+i]);
76: }
77: for (i=0; i < dim; ++i) {
78: for (j=0; j < dim; ++j) {
79: f1[i*dim+j] = mu*(u_x[i*dim+j]+u_x[j*dim+i]);
80: }
81: f1[i*dim+i] += lambda*trace;
82: }
83: }
85: /* 3D elasticity */
86: #define IDX(ii,jj,kk,ll) (27*ii+9*jj+3*kk+ll)
88: void g3_uu_3d_private( PetscScalar g3[], const PetscReal mu, const PetscReal lambda)
89: {
90: if (1) {
91: g3[0] += lambda;
92: g3[0] += mu;
93: g3[0] += mu;
94: g3[4] += lambda;
95: g3[8] += lambda;
96: g3[10] += mu;
97: g3[12] += mu;
98: g3[20] += mu;
99: g3[24] += mu;
100: g3[28] += mu;
101: g3[30] += mu;
102: g3[36] += lambda;
103: g3[40] += lambda;
104: g3[40] += mu;
105: g3[40] += mu;
106: g3[44] += lambda;
107: g3[50] += mu;
108: g3[52] += mu;
109: g3[56] += mu;
110: g3[60] += mu;
111: g3[68] += mu;
112: g3[70] += mu;
113: g3[72] += lambda;
114: g3[76] += lambda;
115: g3[80] += lambda;
116: g3[80] += mu;
117: g3[80] += mu;
118: } else {
119: int i,j,k,l;
120: static int cc=-1;
121: cc++;
122: for (i = 0; i < 3; ++i) {
123: for (j = 0; j < 3; ++j) {
124: for (k = 0; k < 3; ++k) {
125: for (l = 0; l < 3; ++l) {
126: if (k==l && i==j) g3[IDX(i,j,k,l)] += lambda;
127: if (i==k && j==l) g3[IDX(i,j,k,l)] += mu;
128: if (i==l && j==k) g3[IDX(i,j,k,l)] += mu;
129: if (k==l && i==j && !cc) (void) PetscPrintf(PETSC_COMM_WORLD,"g3[%d] += lambda;\n",IDX(i,j,k,l));
130: if (i==k && j==l && !cc) (void) PetscPrintf(PETSC_COMM_WORLD,"g3[%d] += mu;\n",IDX(i,j,k,l));
131: if (i==l && j==k && !cc) (void) PetscPrintf(PETSC_COMM_WORLD,"g3[%d] += mu;\n",IDX(i,j,k,l));
132: }
133: }
134: }
135: }
136: }
137: }
139: static void g3_uu_3d_alpha(PetscInt dim, PetscInt Nf, PetscInt NfAux,
140: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
141: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
142: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
143: {
144: PetscReal mu=s_mu, lambda=s_lambda,rad;
145: PetscInt i;
146: for (i=0,rad=0.;i<dim;i++) {
147: PetscReal t=x[i];
148: rad += t*t;
149: }
150: rad = PetscSqrtReal(rad);
151: if (rad>0.25) {
152: mu *= s_soft_alpha;
153: lambda *= s_soft_alpha; /* we could keep the bulk the same like rubberish */
154: }
155: g3_uu_3d_private(g3,mu,lambda);
156: }
158: static void g3_uu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
159: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
160: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
161: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
162: {
163: g3_uu_3d_private(g3,s_mu,s_lambda);
164: }
166: static void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
167: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
168: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
169: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
170: {
171: const PetscInt Ncomp = dim;
172: PetscInt comp;
174: for (comp = 0; comp < Ncomp; ++comp) f0[comp] = 0.0;
175: }
177: /* PI_i (x_i^4 - x_i^2) */
178: static void f0_u_x4(PetscInt dim, PetscInt Nf, PetscInt NfAux,
179: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
180: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
181: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
182: {
183: const PetscInt Ncomp = dim;
184: PetscInt comp,i;
186: for (comp = 0; comp < Ncomp; ++comp) {
187: f0[comp] = 1e5;
188: for (i = 0; i < Ncomp; ++i) {
189: f0[comp] *= /* (comp+1)* */(x[i]*x[i]*x[i]*x[i] - x[i]*x[i]); /* assumes (0,1]^D domain */
190: }
191: }
192: }
194: PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
195: {
196: const PetscInt Ncomp = dim;
197: PetscInt comp;
199: for (comp = 0; comp < Ncomp; ++comp) u[comp] = 0;
200: return 0;
201: }
203: int main(int argc,char **args)
204: {
205: Mat Amat;
206: PetscErrorCode ierr;
207: SNES snes;
208: KSP ksp;
209: MPI_Comm comm;
210: PetscMPIInt rank;
211: #if defined(PETSC_USE_LOG)
212: PetscLogStage stage[17];
213: #endif
214: PetscBool test_nonzero_cols = PETSC_FALSE,use_nearnullspace = PETSC_TRUE,attach_nearnullspace = PETSC_FALSE;
215: Vec xx,bb;
216: PetscInt iter,i,N,dim = 3,cells[3] = {1,1,1},max_conv_its,local_sizes[7],run_type = 1;
217: DM dm,distdm,basedm;
218: PetscBool flg;
219: char convType[256];
220: PetscReal Lx,mdisp[10],err[10];
221: const char * const options[10] = {"-ex56_dm_refine 0",
222: "-ex56_dm_refine 1",
223: "-ex56_dm_refine 2",
224: "-ex56_dm_refine 3",
225: "-ex56_dm_refine 4",
226: "-ex56_dm_refine 5",
227: "-ex56_dm_refine 6",
228: "-ex56_dm_refine 7",
229: "-ex56_dm_refine 8",
230: "-ex56_dm_refine 9"};
232: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
233: comm = PETSC_COMM_WORLD;
234: MPI_Comm_rank(comm, &rank);
235: /* options */
236: PetscOptionsBegin(comm,NULL,"3D bilinear Q1 elasticity options","");
237: {
238: i = 3;
239: PetscOptionsIntArray("-cells", "Number of (flux tube) processor in each dimension", "ex56.c", cells, &i, NULL);
241: Lx = 1.; /* or ne for rod */
242: max_conv_its = 3;
243: PetscOptionsInt("-max_conv_its","Number of iterations in convergence study","",max_conv_its,&max_conv_its,NULL);
244: if (max_conv_its<=0 || max_conv_its>7) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_USER, "Bad number of iterations for convergence test (%D)",max_conv_its);
245: PetscOptionsReal("-lx","Length of domain","",Lx,&Lx,NULL);
246: PetscOptionsReal("-alpha","material coefficient inside circle","",s_soft_alpha,&s_soft_alpha,NULL);
247: PetscOptionsBool("-test_nonzero_cols","nonzero test","",test_nonzero_cols,&test_nonzero_cols,NULL);
248: PetscOptionsBool("-use_mat_nearnullspace","MatNearNullSpace API test","",use_nearnullspace,&use_nearnullspace,NULL);
249: PetscOptionsBool("-attach_mat_nearnullspace","MatNearNullSpace API test (via MatSetNearNullSpace)","",attach_nearnullspace,&attach_nearnullspace,NULL);
250: PetscOptionsInt("-run_type","0: twisting load on cantalever, 1: 3rd order accurate convergence test","",run_type,&run_type,NULL);
251: }
252: PetscOptionsEnd();
253: PetscLogStageRegister("Mesh Setup", &stage[16]);
254: for (iter=0 ; iter<max_conv_its ; iter++) {
255: char str[] = "Solve 0";
256: str[6] += iter;
257: PetscLogStageRegister(str, &stage[iter]);
258: }
259: /* create DM, Plex calls DMSetup */
260: PetscLogStagePush(stage[16]);
261: DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, cells, NULL, NULL, NULL, PETSC_TRUE, &dm);
262: {
263: DMLabel label;
264: IS is;
265: DMCreateLabel(dm, "boundary");
266: DMGetLabel(dm, "boundary", &label);
267: DMPlexMarkBoundaryFaces(dm, 1, label);
268: if (!run_type) {
269: DMGetStratumIS(dm, "boundary", 1, &is);
270: DMCreateLabel(dm,"Faces");
271: if (is) {
272: PetscInt d, f, Nf;
273: const PetscInt *faces;
274: PetscInt csize;
275: PetscSection cs;
276: Vec coordinates ;
277: DM cdm;
278: ISGetLocalSize(is, &Nf);
279: ISGetIndices(is, &faces);
280: DMGetCoordinatesLocal(dm, &coordinates);
281: DMGetCoordinateDM(dm, &cdm);
282: DMGetLocalSection(cdm, &cs);
283: /* Check for each boundary face if any component of its centroid is either 0.0 or 1.0 */
284: for (f = 0; f < Nf; ++f) {
285: PetscReal faceCoord;
286: PetscInt b,v;
287: PetscScalar *coords = NULL;
288: PetscInt Nv;
289: DMPlexVecGetClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
290: Nv = csize/dim; /* Calculate mean coordinate vector */
291: for (d = 0; d < dim; ++d) {
292: faceCoord = 0.0;
293: for (v = 0; v < Nv; ++v) faceCoord += PetscRealPart(coords[v*dim+d]);
294: faceCoord /= Nv;
295: for (b = 0; b < 2; ++b) {
296: if (PetscAbs(faceCoord - b) < PETSC_SMALL) { /* domain have not been set yet, still [0,1]^3 */
297: DMSetLabelValue(dm, "Faces", faces[f], d*2+b+1);
298: }
299: }
300: }
301: DMPlexVecRestoreClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
302: }
303: ISRestoreIndices(is, &faces);
304: }
305: ISDestroy(&is);
306: DMGetLabel(dm, "Faces", &label);
307: DMPlexLabelComplete(dm, label);
308: }
309: }
310: {
311: PetscInt dimEmbed, i;
312: PetscInt nCoords;
313: PetscScalar *coords,bounds[] = {0,1,-.5,.5,-.5,.5,}; /* x_min,x_max,y_min,y_max */
314: Vec coordinates;
315: bounds[1] = Lx;
316: if (run_type==1) {
317: for (i = 0; i < 2*dim; i++) bounds[i] = (i%2) ? 1 : 0;
318: }
319: DMGetCoordinatesLocal(dm,&coordinates);
320: DMGetCoordinateDim(dm,&dimEmbed);
321: if (dimEmbed != dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"dimEmbed != dim %D",dimEmbed);
322: VecGetLocalSize(coordinates,&nCoords);
323: if (nCoords % dimEmbed) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Coordinate vector the wrong size");
324: VecGetArray(coordinates,&coords);
325: for (i = 0; i < nCoords; i += dimEmbed) {
326: PetscInt j;
327: PetscScalar *coord = &coords[i];
328: for (j = 0; j < dimEmbed; j++) {
329: coord[j] = bounds[2 * j] + coord[j] * (bounds[2 * j + 1] - bounds[2 * j]);
330: }
331: }
332: VecRestoreArray(coordinates,&coords);
333: DMSetCoordinatesLocal(dm,coordinates);
334: }
336: /* convert to p4est, and distribute */
338: PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
339: PetscOptionsFList("-dm_type","Convert DMPlex to another format (should not be Plex!)","ex56.c",DMList,DMPLEX,convType,256,&flg);
340: PetscOptionsEnd();
341: if (flg) {
342: DM newdm;
343: DMConvert(dm,convType,&newdm);
344: if (newdm) {
345: const char *prefix;
346: PetscBool isForest;
347: PetscObjectGetOptionsPrefix((PetscObject)dm,&prefix);
348: PetscObjectSetOptionsPrefix((PetscObject)newdm,prefix);
349: DMIsForest(newdm,&isForest);
350: if (!isForest) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER, "Converted to non Forest?");
351: DMDestroy(&dm);
352: dm = newdm;
353: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER, "Convert failed?");
354: } else {
355: PetscPartitioner part;
356: /* Plex Distribute mesh over processes */
357: DMPlexGetPartitioner(dm,&part);
358: PetscPartitionerSetFromOptions(part);
359: DMPlexDistribute(dm, 0, NULL, &distdm);
360: if (distdm) {
361: const char *prefix;
362: PetscObjectGetOptionsPrefix((PetscObject)dm,&prefix);
363: PetscObjectSetOptionsPrefix((PetscObject)distdm,prefix);
364: DMDestroy(&dm);
365: dm = distdm;
366: }
367: }
368: PetscLogStagePop();
369: basedm = dm; dm = NULL;
371: for (iter=0 ; iter<max_conv_its ; iter++) {
372: PetscLogStagePush(stage[16]);
373: /* make new DM */
374: DMClone(basedm, &dm);
375: PetscObjectSetOptionsPrefix((PetscObject) dm, "ex56_");
376: PetscObjectSetName( (PetscObject)dm,"Mesh");
377: if (max_conv_its > 1) {
378: /* If max_conv_its == 1, then we are not doing a convergence study. */
379: PetscOptionsInsertString(NULL,options[iter]);
380: }
381: DMSetFromOptions(dm); /* refinement done here in Plex, p4est */
382: /* snes */
383: SNESCreate(comm, &snes);
384: SNESSetDM(snes, dm);
385: /* fem */
386: {
387: const PetscInt Ncomp = dim;
388: const PetscInt components[] = {0,1,2};
389: const PetscInt Nfid = 1, Npid = 1;
390: const PetscInt fid[] = {1}; /* The fixed faces (x=0) */
391: const PetscInt pid[] = {2}; /* The faces with loading (x=L_x) */
392: PetscFE fe;
393: PetscDS prob;
394: DMLabel label;
395: DM cdm = dm;
397: PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, dim, PETSC_FALSE, NULL, PETSC_DECIDE, &fe); /* elasticity */
398: PetscObjectSetName((PetscObject) fe, "deformation");
399: /* FEM prob */
400: DMSetField(dm, 0, NULL, (PetscObject) fe);
401: DMCreateDS(dm);
402: DMGetDS(dm, &prob);
403: /* setup problem */
404: if (run_type==1) {
405: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu_3d);
406: PetscDSSetResidual(prob, 0, f0_u_x4, f1_u_3d);
407: } else {
408: PetscWeakForm wf;
409: PetscInt bd, i;
411: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu_3d_alpha);
412: PetscDSSetResidual(prob, 0, f0_u, f1_u_3d_alpha);
414: DMGetLabel(dm, "Faces", &label);
415: DMAddBoundary(dm, DM_BC_NATURAL, "traction", label, Npid, pid, 0, Ncomp, components, NULL, NULL, NULL, &bd);
416: PetscDSGetBoundary(prob, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
417: for (i = 0; i < Npid; ++i) {PetscWeakFormSetIndexBdResidual(wf, label, pid[i], 0, 0, 0, f0_bd_u_3d, 0, f1_bd_u);}
418: }
419: /* bcs */
420: if (run_type==1) {
421: PetscInt id = 1;
422: DMGetLabel(dm, "boundary", &label);
423: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) zero, NULL, NULL, NULL);
424: } else {
425: DMGetLabel(dm, "Faces", &label);
426: DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed", label, Nfid, fid, 0, Ncomp, components, (void (*)(void)) zero, NULL, NULL, NULL);
427: }
428: while (cdm) {
429: DMCopyDisc(dm, cdm);
430: DMGetCoarseDM(cdm, &cdm);
431: }
432: PetscFEDestroy(&fe);
433: }
434: /* vecs & mat */
435: DMCreateGlobalVector(dm,&xx);
436: VecDuplicate(xx, &bb);
437: PetscObjectSetName((PetscObject) bb, "b");
438: PetscObjectSetName((PetscObject) xx, "u");
439: DMCreateMatrix(dm, &Amat);
440: MatSetOption(Amat,MAT_SYMMETRIC,PETSC_TRUE); /* Some matrix kernels can take advantage of symmetry if we set this. */
441: MatSetOption(Amat,MAT_SYMMETRY_ETERNAL,PETSC_TRUE); /* Inform PETSc that Amat is always symmetric, so info set above isn't lost. */
442: MatSetBlockSize(Amat,3);
443: MatSetOption(Amat,MAT_SPD,PETSC_TRUE);
444: VecGetSize(bb,&N);
445: local_sizes[iter] = N;
446: PetscInfo2(snes,"%D global equations, %D vertices\n",N,N/dim);
447: if ((use_nearnullspace || attach_nearnullspace) && N/dim > 1) {
448: /* Set up the near null space (a.k.a. rigid body modes) that will be used by the multigrid preconditioner */
449: DM subdm;
450: MatNullSpace nearNullSpace;
451: PetscInt fields = 0;
452: PetscObject deformation;
453: DMCreateSubDM(dm, 1, &fields, NULL, &subdm);
454: DMPlexCreateRigidBody(subdm, 0, &nearNullSpace);
455: DMGetField(dm, 0, NULL, &deformation);
456: PetscObjectCompose(deformation, "nearnullspace", (PetscObject) nearNullSpace);
457: DMDestroy(&subdm);
458: if (attach_nearnullspace) {
459: MatSetNearNullSpace(Amat,nearNullSpace);
460: }
461: MatNullSpaceDestroy(&nearNullSpace); /* created by DM and destroyed by Mat */
462: }
463: DMPlexSetSNESLocalFEM(dm,NULL,NULL,NULL);
464: SNESSetJacobian(snes, Amat, Amat, NULL, NULL);
465: SNESSetFromOptions(snes);
466: DMSetUp(dm);
467: PetscLogStagePop();
468: PetscLogStagePush(stage[16]);
469: /* ksp */
470: SNESGetKSP(snes, &ksp);
471: KSPSetComputeSingularValues(ksp,PETSC_TRUE);
472: /* test BCs */
473: VecZeroEntries(xx);
474: if (test_nonzero_cols) {
475: if (rank == 0) {
476: VecSetValue(xx,0,1.0,INSERT_VALUES);
477: }
478: VecAssemblyBegin(xx);
479: VecAssemblyEnd(xx);
480: }
481: VecZeroEntries(bb);
482: VecGetSize(bb,&i);
483: local_sizes[iter] = i;
484: PetscInfo2(snes,"%D equations in vector, %D vertices\n",i,i/dim);
485: PetscLogStagePop();
486: /* solve */
487: PetscLogStagePush(stage[iter]);
488: SNESSolve(snes, bb, xx);
489: PetscLogStagePop();
490: VecNorm(xx,NORM_INFINITY,&mdisp[iter]);
491: DMViewFromOptions(dm, NULL, "-dm_view");
492: {
493: PetscViewer viewer = NULL;
494: PetscViewerFormat fmt;
495: PetscOptionsGetViewer(comm,NULL,"ex56_","-vec_view",&viewer,&fmt,&flg);
496: if (flg) {
497: PetscViewerPushFormat(viewer,fmt);
498: VecView(xx,viewer);
499: VecView(bb,viewer);
500: PetscViewerPopFormat(viewer);
501: }
502: PetscViewerDestroy(&viewer);
503: }
504: /* Free work space */
505: DMDestroy(&dm);
506: SNESDestroy(&snes);
507: VecDestroy(&xx);
508: VecDestroy(&bb);
509: MatDestroy(&Amat);
510: }
511: DMDestroy(&basedm);
512: if (run_type==1) err[0] = 59.975208 - mdisp[0]; /* error with what I think is the exact solution */
513: else err[0] = 171.038 - mdisp[0];
514: for (iter=1 ; iter<max_conv_its ; iter++) {
515: if (run_type==1) err[iter] = 59.975208 - mdisp[iter];
516: else err[iter] = 171.038 - mdisp[iter];
517: PetscPrintf(PETSC_COMM_WORLD,"[%d] %D) N=%12D, max displ=%9.7e, disp diff=%9.2e, error=%4.3e, rate=%3.2g\n",rank,iter,local_sizes[iter],(double)mdisp[iter],
518: (double)(mdisp[iter]-mdisp[iter-1]),(double)err[iter],(double)(PetscLogReal(err[iter-1]/err[iter])/PetscLogReal(2.)));
519: }
521: PetscFinalize();
522: return ierr;
523: }
525: /*TEST
527: test:
528: suffix: 0
529: nsize: 4
530: requires: !single
531: args: -cells 2,2,1 -max_conv_its 2 -petscspace_degree 2 -snes_max_it 2 -ksp_max_it 100 -ksp_type cg -ksp_rtol 1.e-10 -ksp_norm_type unpreconditioned -snes_rtol 1.e-10 -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 10 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -ksp_converged_reason -snes_monitor_short -ksp_monitor_short -snes_converged_reason -use_mat_nearnullspace true -mg_levels_ksp_max_it 2 -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.1 -mg_levels_pc_type jacobi -petscpartitioner_type simple -matptap_via scalable -ex56_dm_view
532: timeoutfactor: 2
534: # HYPRE PtAP broken with complex numbers
535: test:
536: suffix: hypre
537: requires: hypre !single !complex !defined(PETSC_HAVE_HYPRE_DEVICE)
538: nsize: 4
539: args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -pc_type hypre -pc_hypre_type boomeramg -pc_hypre_boomeramg_no_CF true -pc_hypre_boomeramg_agg_nl 1 -pc_hypre_boomeramg_coarsen_type HMIS -pc_hypre_boomeramg_interp_type ext+i -ksp_converged_reason -use_mat_nearnullspace true -petscpartitioner_type simple
541: test:
542: suffix: ml
543: requires: ml !single
544: nsize: 4
545: args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_converged_reason -ksp_rtol 1.e-8 -pc_type ml -mg_levels_ksp_type chebyshev -mg_levels_ksp_max_it 3 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -mg_levels_pc_type sor -petscpartitioner_type simple -use_mat_nearnullspace
547: test:
548: suffix: hpddm
549: requires: hpddm slepc !single defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
550: nsize: 4
551: args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type fgmres -ksp_monitor_short -ksp_converged_reason -ksp_rtol 1.e-8 -pc_type hpddm -petscpartitioner_type simple -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 6 -pc_hpddm_coarse_p 1 -pc_hpddm_coarse_pc_type svd
553: test:
554: suffix: repart
555: nsize: 4
556: requires: parmetis !single
557: args: -cells 8,2,2 -max_conv_its 1 -petscspace_degree 2 -snes_max_it 4 -ksp_max_it 100 -ksp_type cg -ksp_rtol 1.e-2 -ksp_norm_type unpreconditioned -snes_rtol 1.e-3 -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -use_mat_nearnullspace true -mg_levels_ksp_max_it 2 -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -mg_levels_pc_type jacobi -pc_gamg_mat_partitioning_type parmetis -pc_gamg_repartition true -snes_converged_reason -pc_gamg_process_eq_limit 20 -pc_gamg_coarse_eq_limit 10 -ksp_converged_reason -snes_converged_reason -pc_gamg_reuse_interpolation true
559: test:
560: suffix: bddc
561: nsize: 4
562: requires: !single
563: args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -ksp_converged_reason -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type {{sbaij baij aij}} -pc_type bddc
565: testset:
566: nsize: 4
567: requires: !single
568: args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-10 -ksp_converged_reason -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type aij -pc_type bddc -attach_mat_nearnullspace {{0 1}separate output}
569: test:
570: suffix: bddc_approx_gamg
571: args: -pc_bddc_switch_static -prefix_push pc_bddc_dirichlet_ -approximate -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -prefix_pop -prefix_push pc_bddc_neumann_ -approximate -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 10 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -prefix_pop
572: # HYPRE PtAP broken with complex numbers
573: test:
574: requires: hypre !complex !defined(PETSC_HAVE_HYPRE_DEVICE)
575: suffix: bddc_approx_hypre
576: args: -pc_bddc_switch_static -prefix_push pc_bddc_dirichlet_ -pc_type hypre -pc_hypre_boomeramg_no_CF true -pc_hypre_boomeramg_strong_threshold 0.75 -pc_hypre_boomeramg_agg_nl 1 -pc_hypre_boomeramg_coarsen_type HMIS -pc_hypre_boomeramg_interp_type ext+i -prefix_pop -prefix_push pc_bddc_neumann_ -pc_type hypre -pc_hypre_boomeramg_no_CF true -pc_hypre_boomeramg_strong_threshold 0.75 -pc_hypre_boomeramg_agg_nl 1 -pc_hypre_boomeramg_coarsen_type HMIS -pc_hypre_boomeramg_interp_type ext+i -prefix_pop
577: test:
578: requires: ml
579: suffix: bddc_approx_ml
580: args: -pc_bddc_switch_static -prefix_push pc_bddc_dirichlet_ -approximate -pc_type ml -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -prefix_pop -prefix_push pc_bddc_neumann_ -approximate -pc_type ml -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -prefix_pop
582: test:
583: suffix: fetidp
584: nsize: 4
585: requires: !single
586: args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type fetidp -fetidp_ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -ksp_converged_reason -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type {{sbaij baij aij}}
588: test:
589: suffix: bddc_elast
590: nsize: 4
591: requires: !single
592: args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -ksp_converged_reason -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type sbaij -pc_type bddc -pc_bddc_monolithic -attach_mat_nearnullspace
594: test:
595: suffix: fetidp_elast
596: nsize: 4
597: requires: !single
598: args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -petscspace_degree 2 -ksp_type fetidp -fetidp_ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -ksp_converged_reason -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type sbaij -fetidp_bddc_pc_bddc_monolithic -attach_mat_nearnullspace
600: testset:
601: nsize: 4
602: requires: !single
603: args: -cells 2,2,1 -max_conv_its 2 -petscspace_degree 2 -snes_max_it 2 -ksp_max_it 100 -ksp_type cg -ksp_rtol 1.e-10 -ksp_norm_type unpreconditioned -snes_rtol 1.e-10 -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 10 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -use_mat_nearnullspace true -mg_levels_ksp_max_it 2 -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -mg_levels_pc_type jacobi -ksp_monitor_short -ksp_converged_reason -snes_converged_reason -snes_monitor_short -ex56_dm_view -petscpartitioner_type simple -pc_gamg_process_eq_limit 20
604: output_file: output/ex56_cuda.out
606: test:
607: suffix: cuda
608: requires: cuda
609: args: -ex56_dm_mat_type aijcusparse -ex56_dm_vec_type cuda
611: test:
612: suffix: viennacl
613: requires: viennacl
614: args: -ex56_dm_mat_type aijviennacl -ex56_dm_vec_type viennacl
616: test:
617: suffix: kokkos
618: requires: kokkos_kernels
619: args: -ex56_dm_mat_type aijkokkos -ex56_dm_vec_type kokkos
620: # Don't run AIJMKL caes with complex scalars because of convergence issues.
621: # Note that we need to test both single and multiple MPI rank cases, because these use different sparse MKL routines to implement the PtAP operation.
622: test:
623: suffix: seqaijmkl
624: nsize: 1
625: requires: defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE) !single !complex
626: args: -cells 2,2,1 -max_conv_its 2 -petscspace_degree 2 -snes_max_it 2 -ksp_max_it 100 -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned -snes_rtol 1.e-10 -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 1000 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -ksp_converged_reason -snes_monitor_short -ksp_monitor_short -snes_converged_reason -use_mat_nearnullspace true -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.1 -mg_levels_pc_type jacobi -petscpartitioner_type simple -mat_block_size 3 -ex56_dm_view -run_type 1 -mat_seqaij_type seqaijmkl
627: timeoutfactor: 2
629: test:
630: suffix: mpiaijmkl
631: nsize: 2
632: requires: defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE) !single !complex
633: args: -cells 2,2,1 -max_conv_its 2 -petscspace_degree 2 -snes_max_it 2 -ksp_max_it 100 -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned -snes_rtol 1.e-10 -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 1000 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -ksp_converged_reason -snes_monitor_short -ksp_monitor_short -snes_converged_reason -use_mat_nearnullspace true -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.1 -mg_levels_pc_type jacobi -petscpartitioner_type simple -mat_block_size 3 -ex56_dm_view -run_type 1 -mat_seqaij_type seqaijmkl
634: timeoutfactor: 2
636: TEST*/