Actual source code: ex56.c

  1: static char help[] = "3D, tri-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 coefficient in embedded circle\n\n";

  8: #include <petscksp.h>

 10: static PetscBool      log_stages = PETSC_TRUE;
 11: static PetscErrorCode MaybeLogStagePush(PetscLogStage stage)
 12: {
 13:   return log_stages ? PetscLogStagePush(stage) : 0;
 14: }
 15: static PetscErrorCode MaybeLogStagePop()
 16: {
 17:   return log_stages ? PetscLogStagePop() : 0;
 18: }
 19: PetscErrorCode elem_3d_elast_v_25(PetscScalar *);

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

 36:   PetscInitialize(&argc, &args, (char *)0, help);
 37:   comm = PETSC_COMM_WORLD;
 38:   MPI_Comm_rank(comm, &mype);
 39:   MPI_Comm_size(comm, &npe);

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

 55:   if (log_stages) {
 56:     PetscLogStageRegister("Setup", &stage[0]);
 57:     PetscLogStageRegister("Solve", &stage[1]);
 58:     PetscLogStageRegister("2nd Setup", &stage[2]);
 59:     PetscLogStageRegister("2nd Solve", &stage[3]);
 60:     PetscLogStageRegister("3rd Setup", &stage[4]);
 61:     PetscLogStageRegister("3rd Solve", &stage[5]);
 62:   } else {
 63:     for (i = 0; i < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(stage); i++) stage[i] = -1;
 64:   }

 66:   h  = 1. / ne;
 67:   nn = ne + 1;
 68:   /* ne*ne; number of global elements */
 69:   M = 3 * nn * nn * nn; /* global number of equations */
 70:   if (npe == 2) {
 71:     if (mype == 1) m = 0;
 72:     else m = nn * nn * nn;
 73:     npe = 1;
 74:   } else {
 75:     m = nn * nn * nn / npe;
 76:     if (mype == npe - 1) m = nn * nn * nn - (npe - 1) * m;
 77:   }
 78:   m *= 3; /* number of equations local*/
 79:   /* Setup solver */
 80:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 81:   KSPSetComputeSingularValues(ksp, PETSC_TRUE);
 82:   KSPSetFromOptions(ksp);
 83:   {
 84:     /* configuration */
 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];

 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:     }

114:     /* create stiffness matrix */
115:     MatCreate(comm, &Amat);
116:     MatSetSizes(Amat, m, m, M, M);
117:     if (!test_late_bs) MatSetBlockSize(Amat, 3);
118:     MatSetType(Amat, MATAIJ);
119:     MatSetOption(Amat, MAT_SPD, PETSC_TRUE);
120:     MatSetOption(Amat, MAT_SPD_ETERNAL, PETSC_TRUE);
121:     MatSetFromOptions(Amat);
122:     MatSeqAIJSetPreallocation(Amat, 0, d_nnz);
123:     MatMPIAIJSetPreallocation(Amat, 0, d_nnz, 0, o_nnz);

125:     PetscFree(d_nnz);
126:     PetscFree(o_nnz);
127:     MatCreateVecs(Amat, &bb, &xx);

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

132:     /* generate element matrices */
133:     {
134:       PetscBool hasData = PETSC_TRUE;
135:       if (!hasData) {
136:         PetscPrintf(PETSC_COMM_WORLD, "\t No data is provided\n");
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:       } else {
144:         /* Get array data */
145:         elem_3d_elast_v_25((PetscScalar *)DD1);
146:       }

148:       /* BC version of element */
149:       for (i = 0; i < 24; i++) {
150:         for (j = 0; j < 24; j++) {
151:           if (i < 12 || (j < 12 && !test_nonzero_cols)) {
152:             if (i == j) DD2[i][j] = 0.1 * DD1[i][j];
153:             else DD2[i][j] = 0.0;
154:           } else DD2[i][j] = DD1[i][j];
155:         }
156:       }
157:       /* element residual/load vector */
158:       for (i = 0; i < 24; i++) {
159:         if (i % 3 == 0) vv[i] = h * h;
160:         else if (i % 3 == 1) vv[i] = 2.0 * h * h;
161:         else vv[i] = .0;
162:       }
163:       for (i = 0; i < 24; i++) {
164:         if (i % 3 == 0 && i >= 12) v2[i] = h * h;
165:         else if (i % 3 == 1 && i >= 12) v2[i] = 2.0 * h * h;
166:         else v2[i] = .0;
167:       }
168:     }

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

173:     /* forms the element stiffness and coordinates */
174:     for (i = Ni0, ic = 0, ii = 0; i < Ni1; i++, ii++) {
175:       for (j = Nj0, jj = 0; j < Nj1; j++, jj++) {
176:         for (k = Nk0, kk = 0; k < Nk1; k++, kk++, ic++) {
177:           /* coords */
178:           x = coords[3 * ic] = h * (PetscReal)i;
179:           y = coords[3 * ic + 1] = h * (PetscReal)j;
180:           z = coords[3 * ic + 2] = h * (PetscReal)k;
181:           /* matrix */
182:           id = id0 + ii + NN * jj + NN * NN * kk;
183:           if (i < ne && j < ne && k < ne) {
184:             /* radius */
185:             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));
186:             PetscReal alpha  = 1.0;
187:             PetscInt  jx, ix, idx[8], idx3[24];
188:             idx[0] = id;
189:             idx[1] = id + 1;
190:             idx[2] = id + NN + 1;
191:             idx[3] = id + NN;
192:             idx[4] = id + NN * NN;
193:             idx[5] = id + 1 + NN * NN;
194:             idx[6] = id + NN + 1 + NN * NN;
195:             idx[7] = id + NN + NN * NN;

197:             /* correct indices */
198:             if (i == Ni1 - 1 && Ni1 != nn) {
199:               idx[1] += NN * (NN * NN - 1);
200:               idx[2] += NN * (NN * NN - 1);
201:               idx[5] += NN * (NN * NN - 1);
202:               idx[6] += NN * (NN * NN - 1);
203:             }
204:             if (j == Nj1 - 1 && Nj1 != nn) {
205:               idx[2] += NN * NN * (nn - 1);
206:               idx[3] += NN * NN * (nn - 1);
207:               idx[6] += NN * NN * (nn - 1);
208:               idx[7] += NN * NN * (nn - 1);
209:             }
210:             if (k == Nk1 - 1 && Nk1 != nn) {
211:               idx[4] += NN * (nn * nn - NN * NN);
212:               idx[5] += NN * (nn * nn - NN * NN);
213:               idx[6] += NN * (nn * nn - NN * NN);
214:               idx[7] += NN * (nn * nn - NN * NN);
215:             }

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

219:             for (ix = 0; ix < 24; ix++) {
220:               for (jx = 0; jx < 24; jx++) DD[ix][jx] = alpha * DD1[ix][jx];
221:             }
222:             if (k > 0) {
223:               if (!test_late_bs) {
224:                 MatSetValuesBlocked(Amat, 8, idx, 8, idx, (const PetscScalar *)DD, ADD_VALUES);
225:                 VecSetValuesBlocked(bb, 8, idx, (const PetscScalar *)vv, ADD_VALUES);
226:               } else {
227:                 for (ix = 0; ix < 8; ix++) {
228:                   idx3[3 * ix]     = 3 * idx[ix];
229:                   idx3[3 * ix + 1] = 3 * idx[ix] + 1;
230:                   idx3[3 * ix + 2] = 3 * idx[ix] + 2;
231:                 }
232:                 MatSetValues(Amat, 24, idx3, 24, idx3, (const PetscScalar *)DD, ADD_VALUES);
233:                 VecSetValues(bb, 24, idx3, (const PetscScalar *)vv, ADD_VALUES);
234:               }
235:             } else {
236:               /* a BC */
237:               for (ix = 0; ix < 24; ix++) {
238:                 for (jx = 0; jx < 24; jx++) DD[ix][jx] = alpha * DD2[ix][jx];
239:               }
240:               if (!test_late_bs) {
241:                 MatSetValuesBlocked(Amat, 8, idx, 8, idx, (const PetscScalar *)DD, ADD_VALUES);
242:                 VecSetValuesBlocked(bb, 8, idx, (const PetscScalar *)v2, ADD_VALUES);
243:               } else {
244:                 for (ix = 0; ix < 8; ix++) {
245:                   idx3[3 * ix]     = 3 * idx[ix];
246:                   idx3[3 * ix + 1] = 3 * idx[ix] + 1;
247:                   idx3[3 * ix + 2] = 3 * idx[ix] + 2;
248:                 }
249:                 MatSetValues(Amat, 24, idx3, 24, idx3, (const PetscScalar *)DD, ADD_VALUES);
250:                 VecSetValues(bb, 24, idx3, (const PetscScalar *)v2, ADD_VALUES);
251:               }
252:             }
253:           }
254:         }
255:       }
256:     }
257:     MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
258:     MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
259:     VecAssemblyBegin(bb);
260:     VecAssemblyEnd(bb);
261:   }
262:   MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
263:   MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
264:   VecAssemblyBegin(bb);
265:   VecAssemblyEnd(bb);
266:   if (test_late_bs) {
267:     VecSetBlockSize(xx, 3);
268:     VecSetBlockSize(bb, 3);
269:     MatSetBlockSize(Amat, 3);
270:   }

272:   if (!PETSC_TRUE) {
273:     PetscViewer viewer;
274:     PetscViewerASCIIOpen(comm, "Amat.m", &viewer);
275:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
276:     MatView(Amat, viewer);
277:     PetscViewerPopFormat(viewer);
278:     PetscViewerDestroy(&viewer);
279:   }

281:   /* finish KSP/PC setup */
282:   KSPSetOperators(ksp, Amat, Amat);
283:   if (use_nearnullspace) {
284:     MatNullSpace matnull;
285:     Vec          vec_coords;
286:     PetscScalar *c;

288:     VecCreate(MPI_COMM_WORLD, &vec_coords);
289:     VecSetBlockSize(vec_coords, 3);
290:     VecSetSizes(vec_coords, m, PETSC_DECIDE);
291:     VecSetUp(vec_coords);
292:     VecGetArray(vec_coords, &c);
293:     for (i = 0; i < m; i++) c[i] = coords[i]; /* Copy since Scalar type might be Complex */
294:     VecRestoreArray(vec_coords, &c);
295:     MatNullSpaceCreateRigidBody(vec_coords, &matnull);
296:     MatSetNearNullSpace(Amat, matnull);
297:     MatNullSpaceDestroy(&matnull);
298:     VecDestroy(&vec_coords);
299:   } else {
300:     PC pc;
301:     KSPGetPC(ksp, &pc);
302:     PCSetCoordinates(pc, 3, m / 3, coords);
303:   }

305:   MaybeLogStagePush(stage[0]);

307:   /* PC setup basically */
308:   KSPSetUp(ksp);

310:   MaybeLogStagePop();
311:   MaybeLogStagePush(stage[1]);

313:   /* test BCs */
314:   if (test_nonzero_cols) {
315:     VecZeroEntries(xx);
316:     if (mype == 0) VecSetValue(xx, 0, 1.0, INSERT_VALUES);
317:     VecAssemblyBegin(xx);
318:     VecAssemblyEnd(xx);
319:     KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
320:   }

322:   /* 1st solve */
323:   KSPSolve(ksp, bb, xx);

325:   MaybeLogStagePop();

327:   /* 2nd solve */
328:   if (two_solves) {
329:     PetscReal emax, emin;
330:     PetscReal norm, norm2;
331:     Vec       res;

333:     MaybeLogStagePush(stage[2]);
334:     /* PC setup basically */
335:     MatScale(Amat, 100000.0);
336:     KSPSetOperators(ksp, Amat, Amat);
337:     KSPSetUp(ksp);

339:     MaybeLogStagePop();
340:     MaybeLogStagePush(stage[3]);
341:     KSPSolve(ksp, bb, xx);
342:     KSPComputeExtremeSingularValues(ksp, &emax, &emin);

344:     MaybeLogStagePop();
345:     MaybeLogStagePush(stage[4]);

347:     MaybeLogStagePop();
348:     MaybeLogStagePush(stage[5]);

350:     /* 3rd solve */
351:     KSPSolve(ksp, bb, xx);

353:     MaybeLogStagePop();

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

357:     VecDuplicate(xx, &res);
358:     MatMult(Amat, xx, res);
359:     VecAXPY(bb, -1.0, res);
360:     VecDestroy(&res);
361:     VecNorm(bb, NORM_2, &norm);
362:     PetscPrintf(PETSC_COMM_WORLD, "[%d]%s |b-Ax|/|b|=%e, |b|=%e, emax=%e\n", 0, PETSC_FUNCTION_NAME, (double)(norm / norm2), (double)norm2, (double)emax);
363:   }

365:   /* Free work space */
366:   KSPDestroy(&ksp);
367:   VecDestroy(&xx);
368:   VecDestroy(&bb);
369:   MatDestroy(&Amat);
370:   PetscFree(coords);

372:   PetscFinalize();
373:   return 0;
374: }

376: /* Data was previously provided in the file data/elem_3d_elast_v_25.tx */
377: PetscErrorCode elem_3d_elast_v_25(PetscScalar *dd)
378: {
379:   PetscScalar DD[] = {
380:     0.18981481481481474,       5.27777777777777568E-002,  5.27777777777777568E-002,  -5.64814814814814659E-002, -1.38888888888889072E-002, -1.38888888888889089E-002, -8.24074074074073876E-002, -5.27777777777777429E-002, 1.38888888888888725E-002,
381:     4.90740740740740339E-002,  1.38888888888889124E-002,  4.72222222222222071E-002,  4.90740740740740339E-002,  4.72222222222221932E-002,  1.38888888888888968E-002,  -8.24074074074073876E-002, 1.38888888888888673E-002,  -5.27777777777777429E-002,
382:     -7.87037037037036785E-002, -4.72222222222221932E-002, -4.72222222222222071E-002, 1.20370370370370180E-002,  -1.38888888888888742E-002, -1.38888888888888829E-002, 5.27777777777777568E-002,  0.18981481481481474,       5.27777777777777568E-002,
383:     1.38888888888889124E-002,  4.90740740740740269E-002,  4.72222222222221932E-002,  -5.27777777777777637E-002, -8.24074074074073876E-002, 1.38888888888888725E-002,  -1.38888888888889037E-002, -5.64814814814814728E-002, -1.38888888888888985E-002,
384:     4.72222222222221932E-002,  4.90740740740740478E-002,  1.38888888888888968E-002,  -1.38888888888888673E-002, 1.20370370370370058E-002,  -1.38888888888888742E-002, -4.72222222222221932E-002, -7.87037037037036785E-002, -4.72222222222222002E-002,
385:     1.38888888888888742E-002,  -8.24074074074073598E-002, -5.27777777777777568E-002, 5.27777777777777568E-002,  5.27777777777777568E-002,  0.18981481481481474,       1.38888888888889055E-002,  4.72222222222222002E-002,  4.90740740740740269E-002,
386:     -1.38888888888888829E-002, -1.38888888888888829E-002, 1.20370370370370180E-002,  4.72222222222222002E-002,  1.38888888888888985E-002,  4.90740740740740339E-002,  -1.38888888888888985E-002, -1.38888888888888968E-002, -5.64814814814814520E-002,
387:     -5.27777777777777568E-002, 1.38888888888888777E-002,  -8.24074074074073876E-002, -4.72222222222222002E-002, -4.72222222222221932E-002, -7.87037037037036646E-002, 1.38888888888888794E-002,  -5.27777777777777568E-002, -8.24074074074073598E-002,
388:     -5.64814814814814659E-002, 1.38888888888889124E-002,  1.38888888888889055E-002,  0.18981481481481474,       -5.27777777777777568E-002, -5.27777777777777499E-002, 4.90740740740740269E-002,  -1.38888888888889072E-002, -4.72222222222221932E-002,
389:     -8.24074074074073876E-002, 5.27777777777777568E-002,  -1.38888888888888812E-002, -8.24074074074073876E-002, -1.38888888888888742E-002, 5.27777777777777499E-002,  4.90740740740740269E-002,  -4.72222222222221863E-002, -1.38888888888889089E-002,
390:     1.20370370370370162E-002,  1.38888888888888673E-002,  1.38888888888888742E-002,  -7.87037037037036785E-002, 4.72222222222222002E-002,  4.72222222222222071E-002,  -1.38888888888889072E-002, 4.90740740740740269E-002,  4.72222222222222002E-002,
391:     -5.27777777777777568E-002, 0.18981481481481480,       5.27777777777777568E-002,  1.38888888888889020E-002,  -5.64814814814814728E-002, -1.38888888888888951E-002, 5.27777777777777637E-002,  -8.24074074074073876E-002, 1.38888888888888881E-002,
392:     1.38888888888888742E-002,  1.20370370370370232E-002,  -1.38888888888888812E-002, -4.72222222222221863E-002, 4.90740740740740339E-002,  1.38888888888888933E-002,  -1.38888888888888812E-002, -8.24074074074073876E-002, -5.27777777777777568E-002,
393:     4.72222222222222071E-002,  -7.87037037037036924E-002, -4.72222222222222140E-002, -1.38888888888889089E-002, 4.72222222222221932E-002,  4.90740740740740269E-002,  -5.27777777777777499E-002, 5.27777777777777568E-002,  0.18981481481481477,
394:     -4.72222222222222071E-002, 1.38888888888888968E-002,  4.90740740740740131E-002,  1.38888888888888812E-002,  -1.38888888888888708E-002, 1.20370370370370267E-002,  5.27777777777777568E-002,  1.38888888888888812E-002,  -8.24074074074073876E-002,
395:     1.38888888888889124E-002,  -1.38888888888889055E-002, -5.64814814814814589E-002, -1.38888888888888812E-002, -5.27777777777777568E-002, -8.24074074074073737E-002, 4.72222222222222002E-002,  -4.72222222222222002E-002, -7.87037037037036924E-002,
396:     -8.24074074074073876E-002, -5.27777777777777637E-002, -1.38888888888888829E-002, 4.90740740740740269E-002,  1.38888888888889020E-002,  -4.72222222222222071E-002, 0.18981481481481480,       5.27777777777777637E-002,  -5.27777777777777637E-002,
397:     -5.64814814814814728E-002, -1.38888888888889037E-002, 1.38888888888888951E-002,  -7.87037037037036785E-002, -4.72222222222222002E-002, 4.72222222222221932E-002,  1.20370370370370128E-002,  -1.38888888888888725E-002, 1.38888888888888812E-002,
398:     4.90740740740740408E-002,  4.72222222222222002E-002,  -1.38888888888888951E-002, -8.24074074074073876E-002, 1.38888888888888812E-002,  5.27777777777777637E-002,  -5.27777777777777429E-002, -8.24074074074073876E-002, -1.38888888888888829E-002,
399:     -1.38888888888889072E-002, -5.64814814814814728E-002, 1.38888888888888968E-002,  5.27777777777777637E-002,  0.18981481481481480,       -5.27777777777777568E-002, 1.38888888888888916E-002,  4.90740740740740339E-002,  -4.72222222222222210E-002,
400:     -4.72222222222221932E-002, -7.87037037037036924E-002, 4.72222222222222002E-002,  1.38888888888888742E-002,  -8.24074074074073876E-002, 5.27777777777777429E-002,  4.72222222222222002E-002,  4.90740740740740269E-002,  -1.38888888888888951E-002,
401:     -1.38888888888888846E-002, 1.20370370370370267E-002,  1.38888888888888916E-002,  1.38888888888888725E-002,  1.38888888888888725E-002,  1.20370370370370180E-002,  -4.72222222222221932E-002, -1.38888888888888951E-002, 4.90740740740740131E-002,
402:     -5.27777777777777637E-002, -5.27777777777777568E-002, 0.18981481481481480,       -1.38888888888888968E-002, -4.72222222222221932E-002, 4.90740740740740339E-002,  4.72222222222221932E-002,  4.72222222222222071E-002,  -7.87037037037036646E-002,
403:     -1.38888888888888742E-002, 5.27777777777777499E-002,  -8.24074074074073737E-002, 1.38888888888888933E-002,  1.38888888888889020E-002,  -5.64814814814814589E-002, 5.27777777777777568E-002,  -1.38888888888888794E-002, -8.24074074074073876E-002,
404:     4.90740740740740339E-002,  -1.38888888888889037E-002, 4.72222222222222002E-002,  -8.24074074074073876E-002, 5.27777777777777637E-002,  1.38888888888888812E-002,  -5.64814814814814728E-002, 1.38888888888888916E-002,  -1.38888888888888968E-002,
405:     0.18981481481481480,       -5.27777777777777499E-002, 5.27777777777777707E-002,  1.20370370370370180E-002,  1.38888888888888812E-002,  -1.38888888888888812E-002, -7.87037037037036785E-002, 4.72222222222222002E-002,  -4.72222222222222071E-002,
406:     -8.24074074074073876E-002, -1.38888888888888742E-002, -5.27777777777777568E-002, 4.90740740740740616E-002,  -4.72222222222222002E-002, 1.38888888888888846E-002,  1.38888888888889124E-002,  -5.64814814814814728E-002, 1.38888888888888985E-002,
407:     5.27777777777777568E-002,  -8.24074074074073876E-002, -1.38888888888888708E-002, -1.38888888888889037E-002, 4.90740740740740339E-002,  -4.72222222222221932E-002, -5.27777777777777499E-002, 0.18981481481481480,       -5.27777777777777568E-002,
408:     -1.38888888888888673E-002, -8.24074074074073598E-002, 5.27777777777777429E-002,  4.72222222222222002E-002,  -7.87037037037036785E-002, 4.72222222222222002E-002,  1.38888888888888708E-002,  1.20370370370370128E-002,  1.38888888888888760E-002,
409:     -4.72222222222222002E-002, 4.90740740740740478E-002,  -1.38888888888888951E-002, 4.72222222222222071E-002,  -1.38888888888888985E-002, 4.90740740740740339E-002,  -1.38888888888888812E-002, 1.38888888888888881E-002,  1.20370370370370267E-002,
410:     1.38888888888888951E-002,  -4.72222222222222210E-002, 4.90740740740740339E-002,  5.27777777777777707E-002,  -5.27777777777777568E-002, 0.18981481481481477,       1.38888888888888829E-002,  5.27777777777777707E-002,  -8.24074074074073598E-002,
411:     -4.72222222222222140E-002, 4.72222222222222140E-002,  -7.87037037037036646E-002, -5.27777777777777707E-002, -1.38888888888888829E-002, -8.24074074074073876E-002, -1.38888888888888881E-002, 1.38888888888888881E-002,  -5.64814814814814589E-002,
412:     4.90740740740740339E-002,  4.72222222222221932E-002,  -1.38888888888888985E-002, -8.24074074074073876E-002, 1.38888888888888742E-002,  5.27777777777777568E-002,  -7.87037037037036785E-002, -4.72222222222221932E-002, 4.72222222222221932E-002,
413:     1.20370370370370180E-002,  -1.38888888888888673E-002, 1.38888888888888829E-002,  0.18981481481481469,       5.27777777777777429E-002,  -5.27777777777777429E-002, -5.64814814814814659E-002, -1.38888888888889055E-002, 1.38888888888889055E-002,
414:     -8.24074074074074153E-002, -5.27777777777777429E-002, -1.38888888888888760E-002, 4.90740740740740408E-002,  1.38888888888888968E-002,  -4.72222222222222071E-002, 4.72222222222221932E-002,  4.90740740740740478E-002,  -1.38888888888888968E-002,
415:     -1.38888888888888742E-002, 1.20370370370370232E-002,  1.38888888888888812E-002,  -4.72222222222222002E-002, -7.87037037037036924E-002, 4.72222222222222071E-002,  1.38888888888888812E-002,  -8.24074074074073598E-002, 5.27777777777777707E-002,
416:     5.27777777777777429E-002,  0.18981481481481477,       -5.27777777777777499E-002, 1.38888888888889107E-002,  4.90740740740740478E-002,  -4.72222222222221932E-002, -5.27777777777777568E-002, -8.24074074074074153E-002, -1.38888888888888812E-002,
417:     -1.38888888888888846E-002, -5.64814814814814659E-002, 1.38888888888888812E-002,  1.38888888888888968E-002,  1.38888888888888968E-002,  -5.64814814814814520E-002, 5.27777777777777499E-002,  -1.38888888888888812E-002, -8.24074074074073876E-002,
418:     4.72222222222221932E-002,  4.72222222222222002E-002,  -7.87037037037036646E-002, -1.38888888888888812E-002, 5.27777777777777429E-002,  -8.24074074074073598E-002, -5.27777777777777429E-002, -5.27777777777777499E-002, 0.18981481481481474,
419:     -1.38888888888888985E-002, -4.72222222222221863E-002, 4.90740740740740339E-002,  1.38888888888888829E-002,  1.38888888888888777E-002,  1.20370370370370249E-002,  -4.72222222222222002E-002, -1.38888888888888933E-002, 4.90740740740740339E-002,
420:     -8.24074074074073876E-002, -1.38888888888888673E-002, -5.27777777777777568E-002, 4.90740740740740269E-002,  -4.72222222222221863E-002, 1.38888888888889124E-002,  1.20370370370370128E-002,  1.38888888888888742E-002,  -1.38888888888888742E-002,
421:     -7.87037037037036785E-002, 4.72222222222222002E-002,  -4.72222222222222140E-002, -5.64814814814814659E-002, 1.38888888888889107E-002,  -1.38888888888888985E-002, 0.18981481481481474,       -5.27777777777777499E-002, 5.27777777777777499E-002,
422:     4.90740740740740339E-002,  -1.38888888888889055E-002, 4.72222222222221932E-002,  -8.24074074074074153E-002, 5.27777777777777499E-002,  1.38888888888888829E-002,  1.38888888888888673E-002,  1.20370370370370058E-002,  1.38888888888888777E-002,
423:     -4.72222222222221863E-002, 4.90740740740740339E-002,  -1.38888888888889055E-002, -1.38888888888888725E-002, -8.24074074074073876E-002, 5.27777777777777499E-002,  4.72222222222222002E-002,  -7.87037037037036785E-002, 4.72222222222222140E-002,
424:     -1.38888888888889055E-002, 4.90740740740740478E-002,  -4.72222222222221863E-002, -5.27777777777777499E-002, 0.18981481481481469,       -5.27777777777777499E-002, 1.38888888888889072E-002,  -5.64814814814814659E-002, 1.38888888888889003E-002,
425:     5.27777777777777429E-002,  -8.24074074074074153E-002, -1.38888888888888812E-002, -5.27777777777777429E-002, -1.38888888888888742E-002, -8.24074074074073876E-002, -1.38888888888889089E-002, 1.38888888888888933E-002,  -5.64814814814814589E-002,
426:     1.38888888888888812E-002,  5.27777777777777429E-002,  -8.24074074074073737E-002, -4.72222222222222071E-002, 4.72222222222222002E-002,  -7.87037037037036646E-002, 1.38888888888889055E-002,  -4.72222222222221932E-002, 4.90740740740740339E-002,
427:     5.27777777777777499E-002,  -5.27777777777777499E-002, 0.18981481481481474,       4.72222222222222002E-002,  -1.38888888888888985E-002, 4.90740740740740339E-002,  -1.38888888888888846E-002, 1.38888888888888812E-002,  1.20370370370370284E-002,
428:     -7.87037037037036785E-002, -4.72222222222221932E-002, -4.72222222222222002E-002, 1.20370370370370162E-002,  -1.38888888888888812E-002, -1.38888888888888812E-002, 4.90740740740740408E-002,  4.72222222222222002E-002,  1.38888888888888933E-002,
429:     -8.24074074074073876E-002, 1.38888888888888708E-002,  -5.27777777777777707E-002, -8.24074074074074153E-002, -5.27777777777777568E-002, 1.38888888888888829E-002,  4.90740740740740339E-002,  1.38888888888889072E-002,  4.72222222222222002E-002,
430:     0.18981481481481477,       5.27777777777777429E-002,  5.27777777777777568E-002,  -5.64814814814814659E-002, -1.38888888888888846E-002, -1.38888888888888881E-002, -4.72222222222221932E-002, -7.87037037037036785E-002, -4.72222222222221932E-002,
431:     1.38888888888888673E-002,  -8.24074074074073876E-002, -5.27777777777777568E-002, 4.72222222222222002E-002,  4.90740740740740269E-002,  1.38888888888889020E-002,  -1.38888888888888742E-002, 1.20370370370370128E-002,  -1.38888888888888829E-002,
432:     -5.27777777777777429E-002, -8.24074074074074153E-002, 1.38888888888888777E-002,  -1.38888888888889055E-002, -5.64814814814814659E-002, -1.38888888888888985E-002, 5.27777777777777429E-002,  0.18981481481481469,       5.27777777777777429E-002,
433:     1.38888888888888933E-002,  4.90740740740740339E-002,  4.72222222222222071E-002,  -4.72222222222222071E-002, -4.72222222222222002E-002, -7.87037037037036646E-002, 1.38888888888888742E-002,  -5.27777777777777568E-002, -8.24074074074073737E-002,
434:     -1.38888888888888951E-002, -1.38888888888888951E-002, -5.64814814814814589E-002, -5.27777777777777568E-002, 1.38888888888888760E-002,  -8.24074074074073876E-002, -1.38888888888888760E-002, -1.38888888888888812E-002, 1.20370370370370249E-002,
435:     4.72222222222221932E-002,  1.38888888888889003E-002,  4.90740740740740339E-002,  5.27777777777777568E-002,  5.27777777777777429E-002,  0.18981481481481474,       1.38888888888888933E-002,  4.72222222222222071E-002,  4.90740740740740339E-002,
436:     1.20370370370370180E-002,  1.38888888888888742E-002,  1.38888888888888794E-002,  -7.87037037037036785E-002, 4.72222222222222071E-002,  4.72222222222222002E-002,  -8.24074074074073876E-002, -1.38888888888888846E-002, 5.27777777777777568E-002,
437:     4.90740740740740616E-002,  -4.72222222222222002E-002, -1.38888888888888881E-002, 4.90740740740740408E-002,  -1.38888888888888846E-002, -4.72222222222222002E-002, -8.24074074074074153E-002, 5.27777777777777429E-002,  -1.38888888888888846E-002,
438:     -5.64814814814814659E-002, 1.38888888888888933E-002,  1.38888888888888933E-002,  0.18981481481481477,       -5.27777777777777568E-002, -5.27777777777777637E-002, -1.38888888888888742E-002, -8.24074074074073598E-002, -5.27777777777777568E-002,
439:     4.72222222222222002E-002,  -7.87037037037036924E-002, -4.72222222222222002E-002, 1.38888888888888812E-002,  1.20370370370370267E-002,  -1.38888888888888794E-002, -4.72222222222222002E-002, 4.90740740740740478E-002,  1.38888888888888881E-002,
440:     1.38888888888888968E-002,  -5.64814814814814659E-002, -1.38888888888888933E-002, 5.27777777777777499E-002,  -8.24074074074074153E-002, 1.38888888888888812E-002,  -1.38888888888888846E-002, 4.90740740740740339E-002,  4.72222222222222071E-002,
441:     -5.27777777777777568E-002, 0.18981481481481477,       5.27777777777777637E-002,  -1.38888888888888829E-002, -5.27777777777777568E-002, -8.24074074074073598E-002, 4.72222222222222071E-002,  -4.72222222222222140E-002, -7.87037037037036924E-002,
442:     5.27777777777777637E-002,  1.38888888888888916E-002,  -8.24074074074073876E-002, 1.38888888888888846E-002,  -1.38888888888888951E-002, -5.64814814814814589E-002, -4.72222222222222071E-002, 1.38888888888888812E-002,  4.90740740740740339E-002,
443:     1.38888888888888829E-002,  -1.38888888888888812E-002, 1.20370370370370284E-002,  -1.38888888888888881E-002, 4.72222222222222071E-002,  4.90740740740740339E-002,  -5.27777777777777637E-002, 5.27777777777777637E-002,  0.18981481481481477,
444:   };

447:   PetscArraycpy(dd, DD, 576);
448:   return 0;
449: }

451: /*TEST

453:    test:
454:       nsize: 8
455:       args: -ne 13 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_agg_nsmooths 1 -pc_gamg_reuse_interpolation true -two_solves -ksp_converged_reason -use_mat_nearnullspace -mg_levels_ksp_max_it 2 -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.05 -pc_gamg_esteig_ksp_max_it 10 -pc_gamg_threshold -0.01 -pc_gamg_coarse_eq_limit 200 -pc_gamg_process_eq_limit 30 -pc_gamg_repartition false -pc_mg_cycle_type v -pc_gamg_use_parallel_coarse_grid_solver -mg_coarse_pc_type jacobi -mg_coarse_ksp_type cg -ksp_monitor_short -pc_gamg_rank_reduction_factors 2,2
456:       filter: grep -v variant

458:    test:
459:       suffix: 2
460:       nsize: 1
461:       args: -ne 11 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_agg_nsmooths 1 -pc_gamg_reuse_interpolation true -two_solves -ksp_converged_reason -use_mat_nearnullspace  -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.05 -pc_gamg_esteig_ksp_max_it 10 -pc_gamg_asm_use_agg true -mg_levels_sub_pc_type lu -mg_levels_pc_asm_overlap 0 -pc_gamg_threshold -0.01 -pc_gamg_coarse_eq_limit 200 -pc_gamg_process_eq_limit 30 -pc_gamg_repartition false -pc_mg_cycle_type v -pc_gamg_use_parallel_coarse_grid_solver -mg_coarse_pc_type jacobi -mg_coarse_ksp_type cg
462:       filter: grep -v variant

464:    test:
465:       suffix: latebs
466:       filter: grep -v variant
467:       nsize: 8
468:       args: -test_late_bs 0 -ne 9 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_agg_nsmooths 1 -pc_gamg_reuse_interpolation true -two_solves -ksp_converged_reason -use_mat_nearnullspace  -mg_levels_ksp_max_it 2 -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.05 -pc_gamg_esteig_ksp_max_it 10 -pc_gamg_threshold -0.01 -pc_gamg_coarse_eq_limit 200 -pc_gamg_process_eq_limit 30 -pc_gamg_repartition false -pc_mg_cycle_type v -pc_gamg_use_parallel_coarse_grid_solver -mg_coarse_pc_type jacobi -mg_coarse_ksp_type cg -ksp_monitor_short -ksp_view

470:    test:
471:       suffix: latebs-2
472:       filter: grep -v variant
473:       nsize: 8
474:       args: -test_late_bs -ne 9 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_agg_nsmooths 1 -pc_gamg_reuse_interpolation true -two_solves -ksp_converged_reason -use_mat_nearnullspace  -mg_levels_ksp_max_it 2 -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.05 -pc_gamg_esteig_ksp_max_it 10 -pc_gamg_threshold -0.01 -pc_gamg_coarse_eq_limit 200 -pc_gamg_process_eq_limit 30 -pc_gamg_repartition false -pc_mg_cycle_type v -pc_gamg_use_parallel_coarse_grid_solver -mg_coarse_pc_type jacobi -mg_coarse_ksp_type cg -ksp_monitor_short -ksp_view

476:    test:
477:       suffix: ml
478:       nsize: 8
479:       args: -ne 9 -alpha 1.e-3 -ksp_type cg -pc_type ml -mg_levels_ksp_type chebyshev -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -mg_levels_pc_type sor -ksp_monitor_short
480:       requires: ml

482:    test:
483:       suffix: nns
484:       args: -ne 9 -alpha 1.e-3 -ksp_converged_reason -ksp_type cg -ksp_max_it 50 -pc_type gamg -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10 -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 1000 -mg_levels_ksp_type chebyshev -mg_levels_pc_type sor -pc_gamg_reuse_interpolation true -two_solves -use_mat_nearnullspace -pc_gamg_use_sa_esteig 0 -mg_levels_esteig_ksp_max_it 10

486:    test:
487:       suffix: nns_telescope
488:       nsize: 2
489:       args: -use_mat_nearnullspace -ksp_monitor_short -pc_type telescope -pc_telescope_reduction_factor 2 -telescope_pc_type gamg -telescope_pc_gamg_esteig_ksp_type cg -telescope_pc_gamg_esteig_ksp_max_it 10

491:    test:
492:       suffix: nns_gdsw
493:       filter: grep -v "variant HERMITIAN"
494:       nsize: 8
495:       args: -use_mat_nearnullspace -ksp_monitor_short -pc_type mg -pc_mg_levels 2 -pc_mg_adapt_interp_coarse_space gdsw -pc_mg_galerkin -mg_levels_pc_type bjacobi -ne 3 -ksp_view

497:    test:
498:       suffix: seqaijmkl
499:       nsize: 8
500:       requires: mkl_sparse
501:       args: -ne 9 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_agg_nsmooths 1 -pc_gamg_reuse_interpolation true -two_solves -ksp_converged_reason -use_mat_nearnullspace -mg_levels_ksp_max_it 2 -mg_levels_ksp_type chebyshev -mg_levels_pc_type jacobi -mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.05 -pc_gamg_esteig_ksp_max_it 10 -pc_gamg_threshold 0.01 -pc_gamg_coarse_eq_limit 2000 -pc_gamg_process_eq_limit 200 -pc_gamg_repartition false -pc_mg_cycle_type v -ksp_monitor_short -mat_seqaij_type seqaijmkl

503: TEST*/