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