Actual source code: ex54.c
petsc-3.11.4 2019-09-28
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 coeficient in embedded circle\n\n";
6: #include <petscksp.h>
8: int main(int argc,char **args)
9: {
10: Mat Amat,Pmat;
12: PetscInt i,m,M,its,Istart,Iend,j,Ii,ix,ne=4;
13: PetscReal x,y,h;
14: Vec xx,bb;
15: KSP ksp;
16: PetscReal soft_alpha = 1.e-3;
17: MPI_Comm comm;
18: PetscMPIInt npe,mype;
19: PetscScalar DD[4][4],DD2[4][4];
20: #if defined(PETSC_USE_LOG)
21: PetscLogStage stage;
22: #endif
23: #define DIAG_S 0.0
24: PetscScalar DD1[4][4] = { {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} };
29: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
30: comm = PETSC_COMM_WORLD;
31: MPI_Comm_rank(comm, &mype);
32: MPI_Comm_size(comm, &npe);
33: PetscOptionsGetInt(NULL,NULL,"-ne",&ne,NULL);
34: h = 1./ne;
35: /* ne*ne; number of global elements */
36: PetscOptionsGetReal(NULL,NULL,"-alpha",&soft_alpha,NULL);
37: M = (ne+1)*(ne+1); /* global number of nodes */
39: /* create stiffness matrix (2) */
40: MatCreate(comm,&Amat);
41: MatSetSizes(Amat,PETSC_DECIDE,PETSC_DECIDE,M,M);
42: MatSetType(Amat,MATMPIAIJ);
43: MatSetFromOptions(Amat);
44: MatSeqAIJSetPreallocation(Amat,81,NULL);
45: MatMPIAIJSetPreallocation(Amat,81,NULL,57,NULL);
47: MatCreate(comm,&Pmat);
48: MatSetSizes(Pmat,PETSC_DECIDE,PETSC_DECIDE,M,M);
49: MatSetType(Pmat,MATMPIAIJ);
50: MatSetFromOptions(Pmat);
51: MatSeqAIJSetPreallocation(Pmat,81,NULL);
52: MatMPIAIJSetPreallocation(Pmat,81,NULL,57,NULL);
54: /* vectors */
55: MatCreateVecs(Amat,&bb,&xx);
56: VecSet(bb,.0);
57: /* generate element matrices -- see ex56.c on how to use different data set */
58: {
59: DD1[0][0] = 0.66666666666666663;
60: DD1[0][1] = -0.16666666666666669;
61: DD1[0][2] = -0.33333333333333343;
62: DD1[0][3] = -0.16666666666666666;
63: DD1[1][0] = -0.16666666666666669;
64: DD1[1][1] = 0.66666666666666663;
65: DD1[1][2] = -0.16666666666666666;
66: DD1[1][3] = -0.33333333333333343;
67: DD1[2][0] = -0.33333333333333343;
68: DD1[2][1] = -0.16666666666666666;
69: DD1[2][2] = 0.66666666666666663;
70: DD1[2][3] = -0.16666666666666663;
71: DD1[3][0] = -0.16666666666666666;
72: DD1[3][1] = -0.33333333333333343;
73: DD1[3][2] = -0.16666666666666663;
74: DD1[3][3] = 0.66666666666666663;
76: /* BC version of element */
77: for (i=0;i<4;i++) {
78: for (j=0;j<4;j++) {
79: if (i<2 || j < 2) {
80: if (i==j) DD2[i][j] = .1*DD1[i][j];
81: else DD2[i][j] = 0.0;
82: } else DD2[i][j] = DD1[i][j];
83: }
84: }
85: }
86: {
87: PetscReal *coords;
88: PC pc;
89: /* forms the element stiffness for the Laplacian and coordinates */
90: MatGetOwnershipRange(Amat,&Istart,&Iend);
91: m = Iend-Istart;
92: PetscMalloc1(2*m,&coords);
93: for (Ii=Istart,ix=0; Ii<Iend; Ii++,ix++) {
94: j = Ii/(ne+1); i = Ii%(ne+1);
95: /* coords */
96: x = h*(Ii % (ne+1)); y = h*(Ii/(ne+1));
97: coords[2*ix] = x; coords[2*ix+1] = y;
98: if (i<ne && j<ne) {
99: PetscInt jj,ii,idx[4];
100: /* radius */
101: PetscReal radius = PetscSqrtReal((x-.5+h/2)*(x-.5+h/2) + (y-.5+h/2)*(y-.5+h/2));
102: PetscReal alpha = 1.0;
103: idx[0] = Ii; idx[1] = Ii+1; idx[2] = Ii + (ne+1) + 1; idx[3] = Ii + (ne+1);
104: if (radius < 0.25) alpha = soft_alpha;
105: for (ii=0; ii<4; ii++) {
106: for (jj=0; jj<4; jj++) DD[ii][jj] = alpha*DD1[ii][jj];
107: }
108: MatSetValues(Pmat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
109: if (j>0) {
110: MatSetValues(Amat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
111: } else {
112: /* a BC */
113: for (ii=0;ii<4;ii++) {
114: for (jj=0;jj<4;jj++) DD[ii][jj] = alpha*DD2[ii][jj];
115: }
116: MatSetValues(Amat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
117: }
118: }
119: if (j>0) {
120: PetscScalar v = h*h;
121: PetscInt jj = Ii;
122: VecSetValues(bb,1,&jj,&v,INSERT_VALUES);
123: }
124: }
125: MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
126: MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
127: MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
128: MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);
129: VecAssemblyBegin(bb);
130: VecAssemblyEnd(bb);
132: /* Setup solver */
133: KSPCreate(PETSC_COMM_WORLD,&ksp);
134: KSPSetFromOptions(ksp);
136: /* finish KSP/PC setup */
137: KSPSetOperators(ksp, Amat, Amat);
139: KSPGetPC(ksp,&pc);
140: PCSetCoordinates(pc, 2, m, coords);
141: PetscFree(coords);
142: }
144: if (!PETSC_TRUE) {
145: PetscViewer viewer;
146: PetscViewerASCIIOpen(comm, "Amat.m", &viewer);
147: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
148: MatView(Amat,viewer);
149: PetscViewerPopFormat(viewer);
150: PetscViewerDestroy(&viewer);
151: }
153: /* solve */
154: #if defined(PETSC_USE_LOG)
155: PetscLogStageRegister("Solve", &stage);
156: PetscLogStagePush(stage);
157: #endif
158: VecSet(xx,.0);
160: KSPSetUp(ksp);
162: KSPSolve(ksp,bb,xx);
164: #if defined(PETSC_USE_LOG)
165: PetscLogStagePop();
166: #endif
168: KSPGetIterationNumber(ksp,&its);
170: if (!PETSC_TRUE) {
171: PetscReal norm,norm2;
172: PetscViewer viewer;
173: Vec res;
174: PetscViewerASCIIOpen(comm, "rhs.m", &viewer);
175: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
176: VecView(bb,viewer);
177: PetscViewerPopFormat(viewer);
178: PetscViewerDestroy(&viewer);
179: VecNorm(bb, NORM_2, &norm2);
181: PetscViewerASCIIOpen(comm, "solution.m", &viewer);
182: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
183: VecView(xx,viewer);
184: PetscViewerPopFormat(viewer);
185: PetscViewerDestroy(&viewer);
187: VecDuplicate(xx, &res);
188: MatMult(Amat, xx, res);
189: VecAXPY(bb, -1.0, res);
190: VecDestroy(&res);
191: VecNorm(bb,NORM_2,&norm);
192: PetscPrintf(PETSC_COMM_WORLD,"[%d]%s |b-Ax|/|b|=%e, |b|=%e\n",0,PETSC_FUNCTION_NAME,norm/norm2,norm2);
194: PetscViewerASCIIOpen(comm, "residual.m", &viewer);
195: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
196: VecView(bb,viewer);
197: PetscViewerPopFormat(viewer);
198: PetscViewerDestroy(&viewer);
199: }
201: /* Free work space */
202: KSPDestroy(&ksp);
203: VecDestroy(&xx);
204: VecDestroy(&bb);
205: MatDestroy(&Amat);
206: MatDestroy(&Pmat);
208: PetscFinalize();
209: return ierr;
210: }
214: /*TEST
216: build:
217: requires: !complex
219: test:
220: nsize: 4
221: args: -ne 49 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -ksp_converged_reason -mg_levels_esteig_ksp_type cg
223: test:
224: suffix: seqaijmkl
225: nsize: 4
226: requires: mkl_sparse
227: args: -ne 19 -alpha 1.e-3 -pc_type gamg -pc_gamg_agg_nsmooths 1 -mg_levels_ksp_max_it 3 -ksp_monitor -ksp_converged_reason -ksp_type cg -mat_seqaij_type seqaijmkl
229: test:
230: suffix: Classical
231: args: -ne 49 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_type classical -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -ksp_converged_reason -mg_levels_esteig_ksp_type cg
232: output_file: output/ex54_classical.out
234: test:
235: suffix: geo
236: nsize: 4
237: args: -ne 49 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_type geo -pc_gamg_coarse_eq_limit 200 -mg_levels_pc_type jacobi -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -ksp_monitor_short -mg_levels_esteig_ksp_type cg -mg_levels_ksp_max_it 3
238: requires: triangle
239: output_file: output/ex54_0.out
241: TEST*/