Actual source code: ex55.c

  1: static char help[] = "2D, bi-linear quadrilateral (Q1), displacement finite element formulation\n\
  2: of plain strain 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 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: int main(int argc, char **args)
 11: {
 12:   Mat         Amat;
 13:   PetscInt    i, m, M, its, Istart, Iend, j, Ii, ix, ne = 4;
 14:   PetscReal   x, y, h;
 15:   Vec         xx, bb;
 16:   KSP         ksp;
 17:   PetscReal   soft_alpha = 1.e-3;
 18:   MPI_Comm    comm;
 19:   PetscBool   use_coords = PETSC_FALSE;
 20:   PetscMPIInt npe, mype;
 21:   PetscScalar DD[8][8], DD2[8][8];
 22: #if defined(PETSC_USE_LOG)
 23:   PetscLogStage stage[2];
 24: #endif
 25:   PetscScalar DD1[8][8] = {
 26:     {5.333333333333333E-01,  2.0000E-01,             -3.333333333333333E-01, 0.0000E+00,             -2.666666666666667E-01, -2.0000E-01,            6.666666666666667E-02,  0.0000E-00            },
 27:     {2.0000E-01,             5.333333333333333E-01,  0.0000E-00,             6.666666666666667E-02,  -2.0000E-01,            -2.666666666666667E-01, 0.0000E-00,             -3.333333333333333E-01},
 28:     {-3.333333333333333E-01, 0.0000E-00,             5.333333333333333E-01,  -2.0000E-01,            6.666666666666667E-02,  0.0000E-00,             -2.666666666666667E-01, 2.0000E-01            },
 29:     {0.0000E+00,             6.666666666666667E-02,  -2.0000E-01,            5.333333333333333E-01,  0.0000E-00,             -3.333333333333333E-01, 2.0000E-01,             -2.666666666666667E-01},
 30:     {-2.666666666666667E-01, -2.0000E-01,            6.666666666666667E-02,  0.0000E-00,             5.333333333333333E-01,  2.0000E-01,             -3.333333333333333E-01, 0.0000E+00            },
 31:     {-2.0000E-01,            -2.666666666666667E-01, 0.0000E-00,             -3.333333333333333E-01, 2.0000E-01,             5.333333333333333E-01,  0.0000E-00,             6.666666666666667E-02 },
 32:     {6.666666666666667E-02,  0.0000E-00,             -2.666666666666667E-01, 2.0000E-01,             -3.333333333333333E-01, 0.0000E-00,             5.333333333333333E-01,  -2.0000E-01           },
 33:     {0.0000E-00,             -3.333333333333333E-01, 2.0000E-01,             -2.666666666666667E-01, 0.0000E-00,             6.666666666666667E-02,  -2.0000E-01,            5.333333333333333E-01 }
 34:   };

 37:   PetscInitialize(&argc, &args, (char *)0, help);
 38:   comm = PETSC_COMM_WORLD;
 39:   MPI_Comm_rank(comm, &mype);
 40:   MPI_Comm_size(comm, &npe);
 41:   PetscOptionsGetInt(NULL, NULL, "-ne", &ne, NULL);
 42:   h = 1. / ne;
 43:   /* ne*ne; number of global elements */
 44:   PetscOptionsGetReal(NULL, NULL, "-alpha", &soft_alpha, NULL);
 45:   PetscOptionsGetBool(NULL, NULL, "-use_coordinates", &use_coords, NULL);
 46:   M = 2 * (ne + 1) * (ne + 1); /* global number of equations */
 47:   m = (ne + 1) * (ne + 1) / npe;
 48:   if (mype == npe - 1) m = (ne + 1) * (ne + 1) - (npe - 1) * m;
 49:   m *= 2;
 50:   /* create stiffness matrix */
 51:   MatCreate(comm, &Amat);
 52:   MatSetSizes(Amat, m, m, M, M);
 53:   MatSetType(Amat, MATAIJ);
 54:   MatSetOption(Amat, MAT_SPD, PETSC_TRUE);
 55:   MatSetOption(Amat, MAT_SPD_ETERNAL, PETSC_TRUE);
 56:   MatSetFromOptions(Amat);
 57:   MatSetBlockSize(Amat, 2);
 58:   MatSeqAIJSetPreallocation(Amat, 18, NULL);
 59:   MatMPIAIJSetPreallocation(Amat, 18, NULL, 18, NULL);
 60: #if defined(PETSC_HAVE_HYPRE)
 61:   MatHYPRESetPreallocation(Amat, 18, NULL, 18, NULL);
 62: #endif

 64:   MatGetOwnershipRange(Amat, &Istart, &Iend);
 66:   /* Generate vectors */
 67:   MatCreateVecs(Amat, &xx, &bb);
 68:   VecSet(bb, .0);
 69:   /* generate element matrices -- see ex56.c on how to use different data set */
 70:   {
 71:     DD[0][0] = 0.53333333333333321;
 72:     DD[0][1] = 0.20000000000000001;
 73:     DD[0][2] = -0.33333333333333331;
 74:     DD[0][3] = 0.0000000000000000;
 75:     DD[0][4] = -0.26666666666666666;
 76:     DD[0][5] = -0.20000000000000001;
 77:     DD[0][6] = 6.66666666666666796E-002;
 78:     DD[0][7] = 6.93889390390722838E-018;
 79:     DD[1][0] = 0.20000000000000001;
 80:     DD[1][1] = 0.53333333333333333;
 81:     DD[1][2] = 7.80625564189563192E-018;
 82:     DD[1][3] = 6.66666666666666935E-002;
 83:     DD[1][4] = -0.20000000000000001;
 84:     DD[1][5] = -0.26666666666666666;
 85:     DD[1][6] = -3.46944695195361419E-018;
 86:     DD[1][7] = -0.33333333333333331;
 87:     DD[2][0] = -0.33333333333333331;
 88:     DD[2][1] = 1.12757025938492461E-017;
 89:     DD[2][2] = 0.53333333333333333;
 90:     DD[2][3] = -0.20000000000000001;
 91:     DD[2][4] = 6.66666666666666935E-002;
 92:     DD[2][5] = -6.93889390390722838E-018;
 93:     DD[2][6] = -0.26666666666666666;
 94:     DD[2][7] = 0.19999999999999998;
 95:     DD[3][0] = 0.0000000000000000;
 96:     DD[3][1] = 6.66666666666666935E-002;
 97:     DD[3][2] = -0.20000000000000001;
 98:     DD[3][3] = 0.53333333333333333;
 99:     DD[3][4] = 4.33680868994201774E-018;
100:     DD[3][5] = -0.33333333333333331;
101:     DD[3][6] = 0.20000000000000001;
102:     DD[3][7] = -0.26666666666666666;
103:     DD[4][0] = -0.26666666666666666;
104:     DD[4][1] = -0.20000000000000001;
105:     DD[4][2] = 6.66666666666666935E-002;
106:     DD[4][3] = 8.67361737988403547E-019;
107:     DD[4][4] = 0.53333333333333333;
108:     DD[4][5] = 0.19999999999999998;
109:     DD[4][6] = -0.33333333333333331;
110:     DD[4][7] = -3.46944695195361419E-018;
111:     DD[5][0] = -0.20000000000000001;
112:     DD[5][1] = -0.26666666666666666;
113:     DD[5][2] = -1.04083408558608426E-017;
114:     DD[5][3] = -0.33333333333333331;
115:     DD[5][4] = 0.19999999999999998;
116:     DD[5][5] = 0.53333333333333333;
117:     DD[5][6] = 6.93889390390722838E-018;
118:     DD[5][7] = 6.66666666666666519E-002;
119:     DD[6][0] = 6.66666666666666796E-002;
120:     DD[6][1] = -6.93889390390722838E-018;
121:     DD[6][2] = -0.26666666666666666;
122:     DD[6][3] = 0.19999999999999998;
123:     DD[6][4] = -0.33333333333333331;
124:     DD[6][5] = 6.93889390390722838E-018;
125:     DD[6][6] = 0.53333333333333321;
126:     DD[6][7] = -0.20000000000000001;
127:     DD[7][0] = 6.93889390390722838E-018;
128:     DD[7][1] = -0.33333333333333331;
129:     DD[7][2] = 0.19999999999999998;
130:     DD[7][3] = -0.26666666666666666;
131:     DD[7][4] = 0.0000000000000000;
132:     DD[7][5] = 6.66666666666666519E-002;
133:     DD[7][6] = -0.20000000000000001;
134:     DD[7][7] = 0.53333333333333321;

136:     /* BC version of element */
137:     for (i = 0; i < 8; i++) {
138:       for (j = 0; j < 8; j++) {
139:         if (i < 4 || j < 4) {
140:           if (i == j) DD2[i][j] = .1 * DD1[i][j];
141:           else DD2[i][j] = 0.0;
142:         } else DD2[i][j] = DD1[i][j];
143:       }
144:     }
145:   }
146:   {
147:     PetscReal *coords;
148:     PetscMalloc1(m, &coords);
149:     /* forms the element stiffness and coordinates */
150:     for (Ii = Istart / 2, ix = 0; Ii < Iend / 2; Ii++, ix++) {
151:       j = Ii / (ne + 1);
152:       i = Ii % (ne + 1);
153:       /* coords */
154:       x                  = h * (Ii % (ne + 1));
155:       y                  = h * (Ii / (ne + 1));
156:       coords[2 * ix]     = x;
157:       coords[2 * ix + 1] = y;
158:       if (i < ne && j < ne) {
159:         PetscInt jj, ii, idx[4];
160:         /* radius */
161:         PetscReal radius = PetscSqrtReal((x - .5 + h / 2) * (x - .5 + h / 2) + (y - .5 + h / 2) * (y - .5 + h / 2));
162:         PetscReal alpha  = 1.0;
163:         if (radius < 0.25) alpha = soft_alpha;

165:         idx[0] = Ii;
166:         idx[1] = Ii + 1;
167:         idx[2] = Ii + (ne + 1) + 1;
168:         idx[3] = Ii + (ne + 1);
169:         for (ii = 0; ii < 8; ii++) {
170:           for (jj = 0; jj < 8; jj++) DD[ii][jj] = alpha * DD1[ii][jj];
171:         }
172:         if (j > 0) {
173:           MatSetValuesBlocked(Amat, 4, idx, 4, idx, (const PetscScalar *)DD, ADD_VALUES);
174:         } else {
175:           /* a BC */
176:           for (ii = 0; ii < 8; ii++) {
177:             for (jj = 0; jj < 8; jj++) DD[ii][jj] = alpha * DD2[ii][jj];
178:           }
179:           MatSetValuesBlocked(Amat, 4, idx, 4, idx, (const PetscScalar *)DD, ADD_VALUES);
180:         }
181:       }
182:       if (j > 0) {
183:         PetscScalar v  = h * h;
184:         PetscInt    jj = 2 * Ii; /* load in x direction */
185:         VecSetValues(bb, 1, &jj, &v, INSERT_VALUES);
186:       }
187:     }
188:     MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
189:     MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
190:     VecAssemblyBegin(bb);
191:     VecAssemblyEnd(bb);

193:     /* Setup solver */
194:     KSPCreate(PETSC_COMM_WORLD, &ksp);
195:     KSPSetFromOptions(ksp);

197:     /* finish KSP/PC setup */
198:     KSPSetOperators(ksp, Amat, Amat);
199:     if (use_coords) {
200:       PC pc;

202:       KSPGetPC(ksp, &pc);
203:       PCSetCoordinates(pc, 2, m / 2, coords);
204:     }
205:     PetscFree(coords);
206:   }

208:   if (!PETSC_TRUE) {
209:     PetscViewer viewer;
210:     PetscViewerASCIIOpen(comm, "Amat.m", &viewer);
211:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
212:     MatView(Amat, viewer);
213:     PetscViewerPopFormat(viewer);
214:     PetscViewerDestroy(&viewer);
215:   }

217:   /* solve */
218: #if defined(PETSC_USE_LOG)
219:   PetscLogStageRegister("Setup", &stage[0]);
220:   PetscLogStageRegister("Solve", &stage[1]);
221:   PetscLogStagePush(stage[0]);
222: #endif
223:   KSPSetUp(ksp);
224: #if defined(PETSC_USE_LOG)
225:   PetscLogStagePop();
226: #endif

228:   VecSet(xx, .0);

230: #if defined(PETSC_USE_LOG)
231:   PetscLogStagePush(stage[1]);
232: #endif
233:   KSPSolve(ksp, bb, xx);
234: #if defined(PETSC_USE_LOG)
235:   PetscLogStagePop();
236: #endif

238:   KSPGetIterationNumber(ksp, &its);

240:   if (0) {
241:     PetscReal   norm, norm2;
242:     PetscViewer viewer;
243:     Vec         res;

245:     PetscObjectGetComm((PetscObject)bb, &comm);
246:     VecNorm(bb, NORM_2, &norm2);

248:     VecDuplicate(xx, &res);
249:     MatMult(Amat, xx, res);
250:     VecAXPY(bb, -1.0, res);
251:     VecDestroy(&res);
252:     VecNorm(bb, NORM_2, &norm);
253:     PetscPrintf(PETSC_COMM_WORLD, "[%d]%s |b-Ax|/|b|=%e, |b|=%e\n", 0, PETSC_FUNCTION_NAME, (double)(norm / norm2), (double)norm2);
254:     PetscViewerASCIIOpen(comm, "residual.m", &viewer);
255:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
256:     VecView(bb, viewer);
257:     PetscViewerPopFormat(viewer);
258:     PetscViewerDestroy(&viewer);
259:   }

261:   /* Free work space */
262:   KSPDestroy(&ksp);
263:   VecDestroy(&xx);
264:   VecDestroy(&bb);
265:   MatDestroy(&Amat);

267:   PetscFinalize();
268:   return 0;
269: }

271: /*TEST

273:    test:
274:       suffix: 1
275:       nsize: 4
276:       args: -ne 29 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -use_coordinates -ksp_converged_reason -pc_gamg_esteig_ksp_max_it 5 -ksp_rtol 1.e-3 -ksp_monitor_short -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.2 -pc_gamg_aggressive_coarsening 0
277:       output_file: output/ex55_sa.out

279:    test:
280:       suffix: Classical
281:       nsize: 4
282:       args: -ne 29 -alpha 1.e-3 -ksp_type gmres -pc_type gamg -pc_gamg_type classical -mg_levels_ksp_max_it 5 -ksp_converged_reason -ksp_rtol 1e-3 -mat_coarsen_misk_distance 2 -pc_gamg_threshold 0
283:       output_file: output/ex55_classical.out

285:    test:
286:       suffix: NC
287:       nsize: 4
288:       args: -ne 29 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -ksp_converged_reason -pc_gamg_esteig_ksp_max_it 10 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.2 -pc_gamg_aggressive_coarsening 0 -pc_gamg_threshold 0

290:    test:
291:       suffix: geo
292:       nsize: 4
293:       args: -ne 29 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_type geo -use_coordinates -ksp_monitor_short -ksp_type cg -ksp_norm_type unpreconditioned  -mg_levels_ksp_max_it 3 -pc_gamg_threshold 0
294:       output_file: output/ex55_0.out
295:       requires: triangle

297:    test:
298:       suffix: hypre
299:       nsize: 4
300:       requires: hypre !complex !defined(PETSC_HAVE_HYPRE_DEVICE)
301:       args: -ne 29 -alpha 1.e-3 -ksp_type cg -pc_type hypre -pc_hypre_type boomeramg -ksp_monitor_short

303:    # command line options match GPU defaults
304:    test:
305:       suffix: hypre_device
306:       nsize: 4
307:       requires: hypre !complex
308:       args: -mat_type hypre -ksp_view -ne 29 -alpha 1.e-3 -ksp_type cg -pc_type hypre -pc_hypre_type boomeramg -ksp_monitor_short -pc_hypre_boomeramg_relax_type_all l1scaled-Jacobi -pc_hypre_boomeramg_interp_type ext+i -pc_hypre_boomeramg_coarsen_type PMIS -pc_hypre_boomeramg_no_CF -pc_mg_galerkin_mat_product_algorithm hypre

310: TEST*/