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