Actual source code: ex54.c
1: static char help[] = "Creates a matrix from quadrilateral finite elements in 2D, Laplacian \n\
2: -ne <size> : problem size in number of elements (eg, -ne 31 gives 32^2 grid)\n\
3: -alpha <v> : scaling of material coefficient in embedded circle\n\n";
5: #include <petscksp.h>
7: int main(int argc, char **args)
8: {
9: Mat Amat, Pmat;
10: PetscInt i, m, M, its, Istart, Iend, j, Ii, ix, ne = 4;
11: PetscReal x, y, h;
12: Vec xx, bb;
13: KSP ksp;
14: PetscReal soft_alpha = 1.e-3;
15: MPI_Comm comm;
16: PetscMPIInt npe, mype;
17: PetscScalar DD[4][4], DD2[4][4];
18: PetscLogStage stage;
19: #define DIAG_S 0.0
20: PetscScalar DD1[4][4] = {
21: {5.0 + DIAG_S, -2.0, -1.0, -2.0 },
22: {-2.0, 5.0 + DIAG_S, -2.0, -1.0 },
23: {-1.0, -2.0, 5.0 + DIAG_S, -2.0 },
24: {-2.0, -1.0, -2.0, 5.0 + DIAG_S}
25: };
27: PetscFunctionBeginUser;
28: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
29: comm = PETSC_COMM_WORLD;
30: PetscCallMPI(MPI_Comm_rank(comm, &mype));
31: PetscCallMPI(MPI_Comm_size(comm, &npe));
32: PetscCall(PetscOptionsGetInt(NULL, NULL, "-ne", &ne, NULL));
33: h = 1. / ne;
34: /* ne*ne; number of global elements */
35: PetscCall(PetscOptionsGetReal(NULL, NULL, "-alpha", &soft_alpha, NULL));
36: M = (ne + 1) * (ne + 1); /* global number of nodes */
38: /* create stiffness matrix (2) */
39: PetscCall(MatCreate(comm, &Amat));
40: PetscCall(MatSetSizes(Amat, PETSC_DECIDE, PETSC_DECIDE, M, M));
41: PetscCall(MatSetType(Amat, MATAIJ));
42: PetscCall(MatSetOption(Amat, MAT_SPD, PETSC_TRUE));
43: PetscCall(MatSetOption(Amat, MAT_SPD_ETERNAL, PETSC_TRUE));
44: PetscCall(MatSetFromOptions(Amat));
45: PetscCall(MatSeqAIJSetPreallocation(Amat, 81, NULL));
46: PetscCall(MatMPIAIJSetPreallocation(Amat, 81, NULL, 57, NULL));
48: PetscCall(MatCreate(comm, &Pmat));
49: PetscCall(MatSetSizes(Pmat, PETSC_DECIDE, PETSC_DECIDE, M, M));
50: PetscCall(MatSetType(Pmat, MATMPIAIJ));
51: PetscCall(MatSetFromOptions(Pmat));
52: PetscCall(MatSeqAIJSetPreallocation(Pmat, 81, NULL));
53: PetscCall(MatMPIAIJSetPreallocation(Pmat, 81, NULL, 57, NULL));
55: /* vectors */
56: PetscCall(MatCreateVecs(Amat, &bb, &xx));
57: PetscCall(VecSet(bb, .0));
58: /* generate element matrices -- see ex56.c on how to use different data set */
59: {
60: DD1[0][0] = 0.66666666666666663;
61: DD1[0][1] = -0.16666666666666669;
62: DD1[0][2] = -0.33333333333333343;
63: DD1[0][3] = -0.16666666666666666;
64: DD1[1][0] = -0.16666666666666669;
65: DD1[1][1] = 0.66666666666666663;
66: DD1[1][2] = -0.16666666666666666;
67: DD1[1][3] = -0.33333333333333343;
68: DD1[2][0] = -0.33333333333333343;
69: DD1[2][1] = -0.16666666666666666;
70: DD1[2][2] = 0.66666666666666663;
71: DD1[2][3] = -0.16666666666666663;
72: DD1[3][0] = -0.16666666666666666;
73: DD1[3][1] = -0.33333333333333343;
74: DD1[3][2] = -0.16666666666666663;
75: DD1[3][3] = 0.66666666666666663;
77: /* BC version of element */
78: for (i = 0; i < 4; i++) {
79: for (j = 0; j < 4; j++) {
80: if (i < 2 || j < 2) {
81: if (i == j) DD2[i][j] = .1 * DD1[i][j];
82: else DD2[i][j] = 0.0;
83: } else DD2[i][j] = DD1[i][j];
84: }
85: }
86: }
87: {
88: PetscReal *coords;
89: PC pc;
90: /* forms the element stiffness for the Laplacian and coordinates */
91: PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
92: m = Iend - Istart;
93: PetscCall(PetscMalloc1(2 * m, &coords));
94: for (Ii = Istart, ix = 0; Ii < Iend; Ii++, ix++) {
95: j = Ii / (ne + 1);
96: i = Ii % (ne + 1);
97: /* coords */
98: x = h * (Ii % (ne + 1));
99: y = h * (Ii / (ne + 1));
100: coords[2 * ix] = x;
101: coords[2 * ix + 1] = y;
102: if (i < ne && j < ne) {
103: PetscInt jj, ii, idx[4];
104: /* radius */
105: PetscReal radius = PetscSqrtReal((x - .5 + h / 2) * (x - .5 + h / 2) + (y - .5 + h / 2) * (y - .5 + h / 2));
106: PetscReal alpha = 1.0;
107: idx[0] = Ii;
108: idx[1] = Ii + 1;
109: idx[2] = Ii + (ne + 1) + 1;
110: idx[3] = Ii + (ne + 1);
111: if (radius < 0.25) alpha = soft_alpha;
112: for (ii = 0; ii < 4; ii++) {
113: for (jj = 0; jj < 4; jj++) DD[ii][jj] = alpha * DD1[ii][jj];
114: }
115: PetscCall(MatSetValues(Pmat, 4, idx, 4, idx, (const PetscScalar *)DD, ADD_VALUES));
116: if (j > 0) {
117: PetscCall(MatSetValues(Amat, 4, idx, 4, idx, (const PetscScalar *)DD, ADD_VALUES));
118: } else {
119: /* a BC */
120: for (ii = 0; ii < 4; ii++) {
121: for (jj = 0; jj < 4; jj++) DD[ii][jj] = alpha * DD2[ii][jj];
122: }
123: PetscCall(MatSetValues(Amat, 4, idx, 4, idx, (const PetscScalar *)DD, ADD_VALUES));
124: }
125: }
126: if (j > 0) {
127: PetscScalar v = h * h;
128: PetscInt jj = Ii;
129: PetscCall(VecSetValues(bb, 1, &jj, &v, INSERT_VALUES));
130: }
131: }
132: PetscCall(MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY));
133: PetscCall(MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY));
134: PetscCall(MatAssemblyBegin(Pmat, MAT_FINAL_ASSEMBLY));
135: PetscCall(MatAssemblyEnd(Pmat, MAT_FINAL_ASSEMBLY));
136: PetscCall(VecAssemblyBegin(bb));
137: PetscCall(VecAssemblyEnd(bb));
139: /* Setup solver */
140: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
141: PetscCall(KSPSetFromOptions(ksp));
143: /* finish KSP/PC setup */
144: PetscCall(KSPSetOperators(ksp, Amat, Amat));
146: PetscCall(KSPGetPC(ksp, &pc));
147: PetscCall(PCSetCoordinates(pc, 2, m, coords));
148: PetscCall(PetscFree(coords));
149: }
151: if (!PETSC_TRUE) {
152: PetscViewer viewer;
153: PetscCall(PetscViewerASCIIOpen(comm, "Amat.m", &viewer));
154: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
155: PetscCall(MatView(Amat, viewer));
156: PetscCall(PetscViewerPopFormat(viewer));
157: PetscCall(PetscViewerDestroy(&viewer));
158: }
160: /* solve */
161: PetscCall(PetscLogStageRegister("Solve", &stage));
162: PetscCall(PetscLogStagePush(stage));
163: PetscCall(VecSet(xx, .0));
165: PetscCall(KSPSetUp(ksp));
167: PetscCall(KSPSolve(ksp, bb, xx));
169: PetscCall(PetscLogStagePop());
171: PetscCall(KSPGetIterationNumber(ksp, &its));
173: if (!PETSC_TRUE) {
174: PetscReal norm, norm2;
175: PetscViewer viewer;
176: Vec res;
177: PetscCall(PetscViewerASCIIOpen(comm, "rhs.m", &viewer));
178: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
179: PetscCall(VecView(bb, viewer));
180: PetscCall(PetscViewerPopFormat(viewer));
181: PetscCall(PetscViewerDestroy(&viewer));
182: PetscCall(VecNorm(bb, NORM_2, &norm2));
184: PetscCall(PetscViewerASCIIOpen(comm, "solution.m", &viewer));
185: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
186: PetscCall(VecView(xx, viewer));
187: PetscCall(PetscViewerPopFormat(viewer));
188: PetscCall(PetscViewerDestroy(&viewer));
190: PetscCall(VecDuplicate(xx, &res));
191: PetscCall(MatMult(Amat, xx, res));
192: PetscCall(VecAXPY(bb, -1.0, res));
193: PetscCall(VecDestroy(&res));
194: PetscCall(VecNorm(bb, NORM_2, &norm));
195: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%d]%s |b-Ax|/|b|=%e, |b|=%e\n", 0, PETSC_FUNCTION_NAME, (double)(norm / norm2), (double)norm2));
197: PetscCall(PetscViewerASCIIOpen(comm, "residual.m", &viewer));
198: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
199: PetscCall(VecView(bb, viewer));
200: PetscCall(PetscViewerPopFormat(viewer));
201: PetscCall(PetscViewerDestroy(&viewer));
202: }
204: /* Free work space */
205: PetscCall(KSPDestroy(&ksp));
206: PetscCall(VecDestroy(&xx));
207: PetscCall(VecDestroy(&bb));
208: PetscCall(MatDestroy(&Amat));
209: PetscCall(MatDestroy(&Pmat));
211: PetscCall(PetscFinalize());
212: return 0;
213: }
215: /*TEST
217: build:
218: requires: !complex
220: test:
221: nsize: 4
222: args: -ne 19 -alpha 1.e-3 -ksp_type cg -pc_type gamg -mg_levels_ksp_max_it 2 -ksp_monitor -ksp_converged_reason -pc_gamg_esteig_ksp_max_it 5 -pc_gamg_esteig_ksp_type cg -mg_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -pc_gamg_aggressive_coarsening 0
224: test:
225: suffix: seqaijmkl
226: nsize: 4
227: requires: mkl_sparse
228: args: -ne 19 -alpha 1.e-3 -ksp_type cg -pc_type gamg -mg_levels_ksp_max_it 2 -ksp_monitor -ksp_converged_reason -pc_gamg_esteig_ksp_max_it 5 -pc_gamg_esteig_ksp_type cg -mg_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -mat_seqaij_type seqaijmkl -pc_gamg_aggressive_coarsening 0
230: test:
231: suffix: Classical
232: args: -ne 49 -alpha 1.e-3 -ksp_type cg -pc_type gamg -mg_levels_ksp_max_it 2 -pc_gamg_type classical -ksp_monitor -ksp_converged_reason -mg_levels_esteig_ksp_type cg -mg_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -mat_coarsen_type mis
233: output_file: output/ex54_classical.out
235: test:
236: suffix: geo
237: nsize: 4
238: args: -ne 49 -alpha 1.e-3 -ksp_type cg -pc_type gamg -mg_levels_ksp_max_it 4 -pc_gamg_type geo -pc_gamg_coarse_eq_limit 200 -mg_levels_esteig_ksp_type cg -mg_levels_esteig_ksp_max_it 10 -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.05 -ksp_monitor_short -ksp_converged_reason -ksp_rtol 1e-3 -ksp_norm_type unpreconditioned
239: requires: triangle
240: output_file: output/ex54_0.out
242: test:
243: requires: !single !__float128
244: nsize: 4
245: suffix: hem
246: args: -ne 39 -ksp_type cg -pc_type gamg -pc_gamg_type agg -ksp_rtol 1e-4 -pc_gamg_square_graph 0 -ksp_monitor_short -ksp_norm_type unpreconditioned -mat_coarsen_type hem -ksp_monitor_short -pc_gamg_aggressive_coarsening 0
247: TEST*/