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: PetscLogStage stage[2];
23: PetscScalar DD1[8][8] = {
24: {5.333333333333333E-01, 2.0000E-01, -3.333333333333333E-01, 0.0000E+00, -2.666666666666667E-01, -2.0000E-01, 6.666666666666667E-02, 0.0000E-00 },
25: {2.0000E-01, 5.333333333333333E-01, 0.0000E-00, 6.666666666666667E-02, -2.0000E-01, -2.666666666666667E-01, 0.0000E-00, -3.333333333333333E-01},
26: {-3.333333333333333E-01, 0.0000E-00, 5.333333333333333E-01, -2.0000E-01, 6.666666666666667E-02, 0.0000E-00, -2.666666666666667E-01, 2.0000E-01 },
27: {0.0000E+00, 6.666666666666667E-02, -2.0000E-01, 5.333333333333333E-01, 0.0000E-00, -3.333333333333333E-01, 2.0000E-01, -2.666666666666667E-01},
28: {-2.666666666666667E-01, -2.0000E-01, 6.666666666666667E-02, 0.0000E-00, 5.333333333333333E-01, 2.0000E-01, -3.333333333333333E-01, 0.0000E+00 },
29: {-2.0000E-01, -2.666666666666667E-01, 0.0000E-00, -3.333333333333333E-01, 2.0000E-01, 5.333333333333333E-01, 0.0000E-00, 6.666666666666667E-02 },
30: {6.666666666666667E-02, 0.0000E-00, -2.666666666666667E-01, 2.0000E-01, -3.333333333333333E-01, 0.0000E-00, 5.333333333333333E-01, -2.0000E-01 },
31: {0.0000E-00, -3.333333333333333E-01, 2.0000E-01, -2.666666666666667E-01, 0.0000E-00, 6.666666666666667E-02, -2.0000E-01, 5.333333333333333E-01 }
32: };
34: PetscFunctionBeginUser;
35: PetscCall(PetscInitialize(&argc, &args, NULL, help));
36: comm = PETSC_COMM_WORLD;
37: PetscCallMPI(MPI_Comm_rank(comm, &mype));
38: PetscCallMPI(MPI_Comm_size(comm, &npe));
39: PetscCall(PetscOptionsGetInt(NULL, NULL, "-ne", &ne, NULL));
40: h = 1. / ne;
41: /* ne*ne; number of global elements */
42: PetscCall(PetscOptionsGetReal(NULL, NULL, "-alpha", &soft_alpha, NULL));
43: PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_coordinates", &use_coords, NULL));
44: M = 2 * (ne + 1) * (ne + 1); /* global number of equations */
45: m = (ne + 1) * (ne + 1) / npe;
46: if (mype == npe - 1) m = (ne + 1) * (ne + 1) - (npe - 1) * m;
47: m *= 2;
48: /* create stiffness matrix */
49: PetscCall(MatCreate(comm, &Amat));
50: PetscCall(MatSetSizes(Amat, m, m, M, M));
51: PetscCall(MatSetType(Amat, MATAIJ));
52: PetscCall(MatSetOption(Amat, MAT_SPD, PETSC_TRUE));
53: PetscCall(MatSetOption(Amat, MAT_SPD_ETERNAL, PETSC_TRUE));
54: PetscCall(MatSetFromOptions(Amat));
55: PetscCall(MatSetBlockSize(Amat, 2));
56: PetscCall(MatSeqAIJSetPreallocation(Amat, 18, NULL));
57: PetscCall(MatMPIAIJSetPreallocation(Amat, 18, NULL, 18, NULL));
58: #if defined(PETSC_HAVE_HYPRE)
59: PetscCall(MatHYPRESetPreallocation(Amat, 18, NULL, 18, NULL));
60: #endif
62: PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
63: PetscCheck(m == Iend - Istart, PETSC_COMM_SELF, PETSC_ERR_PLIB, "m %" PetscInt_FMT " does not equal Iend %" PetscInt_FMT " - Istart %" PetscInt_FMT, m, Iend, Istart);
64: /* Generate vectors */
65: PetscCall(MatCreateVecs(Amat, &xx, &bb));
66: PetscCall(VecSet(bb, .0));
67: /* generate element matrices -- see ex56.c on how to use different data set */
68: {
69: DD[0][0] = 0.53333333333333321;
70: DD[0][1] = 0.20000000000000001;
71: DD[0][2] = -0.33333333333333331;
72: DD[0][3] = 0.0000000000000000;
73: DD[0][4] = -0.26666666666666666;
74: DD[0][5] = -0.20000000000000001;
75: DD[0][6] = 6.66666666666666796E-002;
76: DD[0][7] = 6.93889390390722838E-018;
77: DD[1][0] = 0.20000000000000001;
78: DD[1][1] = 0.53333333333333333;
79: DD[1][2] = 7.80625564189563192E-018;
80: DD[1][3] = 6.66666666666666935E-002;
81: DD[1][4] = -0.20000000000000001;
82: DD[1][5] = -0.26666666666666666;
83: DD[1][6] = -3.46944695195361419E-018;
84: DD[1][7] = -0.33333333333333331;
85: DD[2][0] = -0.33333333333333331;
86: DD[2][1] = 1.12757025938492461E-017;
87: DD[2][2] = 0.53333333333333333;
88: DD[2][3] = -0.20000000000000001;
89: DD[2][4] = 6.66666666666666935E-002;
90: DD[2][5] = -6.93889390390722838E-018;
91: DD[2][6] = -0.26666666666666666;
92: DD[2][7] = 0.19999999999999998;
93: DD[3][0] = 0.0000000000000000;
94: DD[3][1] = 6.66666666666666935E-002;
95: DD[3][2] = -0.20000000000000001;
96: DD[3][3] = 0.53333333333333333;
97: DD[3][4] = 4.33680868994201774E-018;
98: DD[3][5] = -0.33333333333333331;
99: DD[3][6] = 0.20000000000000001;
100: DD[3][7] = -0.26666666666666666;
101: DD[4][0] = -0.26666666666666666;
102: DD[4][1] = -0.20000000000000001;
103: DD[4][2] = 6.66666666666666935E-002;
104: DD[4][3] = 8.67361737988403547E-019;
105: DD[4][4] = 0.53333333333333333;
106: DD[4][5] = 0.19999999999999998;
107: DD[4][6] = -0.33333333333333331;
108: DD[4][7] = -3.46944695195361419E-018;
109: DD[5][0] = -0.20000000000000001;
110: DD[5][1] = -0.26666666666666666;
111: DD[5][2] = -1.04083408558608426E-017;
112: DD[5][3] = -0.33333333333333331;
113: DD[5][4] = 0.19999999999999998;
114: DD[5][5] = 0.53333333333333333;
115: DD[5][6] = 6.93889390390722838E-018;
116: DD[5][7] = 6.66666666666666519E-002;
117: DD[6][0] = 6.66666666666666796E-002;
118: DD[6][1] = -6.93889390390722838E-018;
119: DD[6][2] = -0.26666666666666666;
120: DD[6][3] = 0.19999999999999998;
121: DD[6][4] = -0.33333333333333331;
122: DD[6][5] = 6.93889390390722838E-018;
123: DD[6][6] = 0.53333333333333321;
124: DD[6][7] = -0.20000000000000001;
125: DD[7][0] = 6.93889390390722838E-018;
126: DD[7][1] = -0.33333333333333331;
127: DD[7][2] = 0.19999999999999998;
128: DD[7][3] = -0.26666666666666666;
129: DD[7][4] = 0.0000000000000000;
130: DD[7][5] = 6.66666666666666519E-002;
131: DD[7][6] = -0.20000000000000001;
132: DD[7][7] = 0.53333333333333321;
134: /* BC version of element */
135: for (i = 0; i < 8; i++) {
136: for (j = 0; j < 8; j++) {
137: if (i < 4 || j < 4) {
138: if (i == j) DD2[i][j] = .1 * DD1[i][j];
139: else DD2[i][j] = 0.0;
140: } else DD2[i][j] = DD1[i][j];
141: }
142: }
143: }
144: {
145: PetscReal *coords;
146: PetscCall(PetscMalloc1(m, &coords));
147: /* forms the element stiffness and coordinates */
148: for (Ii = Istart / 2, ix = 0; Ii < Iend / 2; Ii++, ix++) {
149: j = Ii / (ne + 1);
150: i = Ii % (ne + 1);
151: /* coords */
152: x = h * (Ii % (ne + 1));
153: y = h * (Ii / (ne + 1));
154: coords[2 * ix] = x;
155: coords[2 * ix + 1] = y;
156: if (i < ne && j < ne) {
157: PetscInt jj, ii, idx[4];
158: /* radius */
159: PetscReal radius = PetscSqrtReal((x - .5 + h / 2) * (x - .5 + h / 2) + (y - .5 + h / 2) * (y - .5 + h / 2));
160: PetscReal alpha = 1.0;
161: if (radius < 0.25) alpha = soft_alpha;
163: idx[0] = Ii;
164: idx[1] = Ii + 1;
165: idx[2] = Ii + (ne + 1) + 1;
166: idx[3] = Ii + (ne + 1);
167: for (ii = 0; ii < 8; ii++) {
168: for (jj = 0; jj < 8; jj++) DD[ii][jj] = alpha * DD1[ii][jj];
169: }
170: if (j > 0) {
171: PetscCall(MatSetValuesBlocked(Amat, 4, idx, 4, idx, (const PetscScalar *)DD, ADD_VALUES));
172: } else {
173: /* a BC */
174: for (ii = 0; ii < 8; ii++) {
175: for (jj = 0; jj < 8; jj++) DD[ii][jj] = alpha * DD2[ii][jj];
176: }
177: PetscCall(MatSetValuesBlocked(Amat, 4, idx, 4, idx, (const PetscScalar *)DD, ADD_VALUES));
178: }
179: }
180: if (j > 0) {
181: PetscScalar v = h * h;
182: PetscInt jj = 2 * Ii; /* load in x direction */
183: PetscCall(VecSetValues(bb, 1, &jj, &v, INSERT_VALUES));
184: }
185: }
186: PetscCall(MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY));
187: PetscCall(MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY));
188: PetscCall(VecAssemblyBegin(bb));
189: PetscCall(VecAssemblyEnd(bb));
191: /* Setup solver */
192: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
193: PetscCall(KSPSetFromOptions(ksp));
195: /* finish KSP/PC setup */
196: PetscCall(KSPSetOperators(ksp, Amat, Amat));
197: if (use_coords) {
198: PC pc;
200: PetscCall(KSPGetPC(ksp, &pc));
201: PetscCall(PCSetCoordinates(pc, 2, m / 2, coords));
202: }
203: PetscCall(PetscFree(coords));
204: }
206: if (!PETSC_TRUE) {
207: PetscViewer viewer;
208: PetscCall(PetscViewerASCIIOpen(comm, "Amat.m", &viewer));
209: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
210: PetscCall(MatView(Amat, viewer));
211: PetscCall(PetscViewerPopFormat(viewer));
212: PetscCall(PetscViewerDestroy(&viewer));
213: }
215: /* solve */
216: PetscCall(PetscLogStageRegister("Setup", &stage[0]));
217: PetscCall(PetscLogStageRegister("Solve", &stage[1]));
218: PetscCall(PetscLogStagePush(stage[0]));
219: PetscCall(KSPSetUp(ksp));
220: PetscCall(PetscLogStagePop());
222: PetscCall(VecSet(xx, .0));
224: PetscCall(PetscLogStagePush(stage[1]));
225: PetscCall(KSPSolve(ksp, bb, xx));
226: PetscCall(PetscLogStagePop());
228: PetscCall(KSPGetIterationNumber(ksp, &its));
230: if (0) {
231: PetscReal norm, norm2;
232: PetscViewer viewer;
233: Vec res;
235: PetscCall(PetscObjectGetComm((PetscObject)bb, &comm));
236: PetscCall(VecNorm(bb, NORM_2, &norm2));
238: PetscCall(VecDuplicate(xx, &res));
239: PetscCall(MatMult(Amat, xx, res));
240: PetscCall(VecAXPY(bb, -1.0, res));
241: PetscCall(VecDestroy(&res));
242: PetscCall(VecNorm(bb, NORM_2, &norm));
243: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%d]%s |b-Ax|/|b|=%e, |b|=%e\n", 0, PETSC_FUNCTION_NAME, (double)(norm / norm2), (double)norm2));
244: PetscCall(PetscViewerASCIIOpen(comm, "residual.m", &viewer));
245: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
246: PetscCall(VecView(bb, viewer));
247: PetscCall(PetscViewerPopFormat(viewer));
248: PetscCall(PetscViewerDestroy(&viewer));
249: }
251: /* Free work space */
252: PetscCall(KSPDestroy(&ksp));
253: PetscCall(VecDestroy(&xx));
254: PetscCall(VecDestroy(&bb));
255: PetscCall(MatDestroy(&Amat));
257: PetscCall(PetscFinalize());
258: return 0;
259: }
261: /*TEST
263: test:
264: suffix: 1
265: nsize: 4
266: 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
267: output_file: output/ex55_sa.out
269: test:
270: suffix: Classical
271: nsize: 4
272: requires: !complex
273: 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 -pc_gamg_mat_coarsen_misk_distance 2 -pc_gamg_threshold 0
274: output_file: output/ex55_classical.out
275: filter: sed -e "s/Linear solve converged due to CONVERGED_RTOL iterations 17/Linear solve converged due to CONVERGED_RTOL iterations 18/g"
277: test:
278: suffix: NC
279: nsize: 4
280: 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
281: filter: sed -e "s/Linear solve converged due to CONVERGED_RTOL iterations 12/Linear solve converged due to CONVERGED_RTOL iterations 11/g"
283: test:
284: suffix: geo
285: nsize: 4
286: 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
287: output_file: output/ex55_0.out
288: requires: triangle
290: test:
291: suffix: hypre
292: nsize: 4
293: requires: hypre !complex !defined(PETSC_HAVE_HYPRE_DEVICE)
294: args: -ne 29 -alpha 1.e-3 -ksp_type cg -pc_type hypre -pc_hypre_type boomeramg -ksp_monitor_short
296: test:
297: suffix: hypre_ilu
298: nsize: 4
299: requires: hypre !complex !defined(PETSC_HAVE_HYPRE_DEVICE)
300: args: -ne 29 -alpha 1.e-3 -pc_type hypre -pc_hypre_type ilu -ksp_monitor_short
302: test:
303: suffix: kok
304: nsize: 4
305: requires: kokkos_kernels
306: args: -ne 29 -ksp_type cg -pc_type pbjacobi -mat_type aijkokkos -ksp_converged_reason
307: filter: grep -v CONVERGED_RTOL
308: output_file: output/empty.out
310: # command line options match GPU defaults
311: test:
312: suffix: hypre_device
313: nsize: 4
314: requires: hypre !complex
315: 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
317: TEST*/