Actual source code: ex54.c
petsc-3.5.4 2015-05-23
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>
10: int main(int argc,char **args)
11: {
12: Mat Amat,Pmat;
14: PetscInt i,m,M,its,Istart,Iend,j,Ii,ix,ne=4;
15: PetscReal x,y,h;
16: Vec xx,bb;
17: KSP ksp;
18: PetscReal soft_alpha = 1.e-3;
19: MPI_Comm comm;
20: PetscMPIInt npe,mype;
21: PC pc;
22: PetscScalar DD[4][4],DD2[4][4];
23: #if defined(PETSC_USE_LOG)
24: PetscLogStage stage;
25: #endif
26: #define DIAG_S 0.0
27: PetscScalar DD1[4][4] = { {5.0+DIAG_S, -2.0, -1.0, -2.0},
28: {-2.0, 5.0+DIAG_S, -2.0, -1.0},
29: {-1.0, -2.0, 5.0+DIAG_S, -2.0},
30: {-2.0, -1.0, -2.0, 5.0+DIAG_S} };
32: PetscInitialize(&argc,&args,(char*)0,help);
33: comm = PETSC_COMM_WORLD;
34: MPI_Comm_rank(comm, &mype);
35: MPI_Comm_size(comm, &npe);
36: PetscOptionsGetInt(NULL,"-ne",&ne,NULL);
37: h = 1./ne;
38: /* ne*ne; number of global elements */
39: PetscOptionsGetReal(NULL,"-alpha",&soft_alpha,NULL);
40: M = (ne+1)*(ne+1); /* global number of nodes */
41: /* create stiffness matrix */
42: MatCreateAIJ(comm,PETSC_DECIDE,PETSC_DECIDE,M,M,
43: 18,NULL,6,NULL,&Amat);
44: MatCreateAIJ(comm,PETSC_DECIDE,PETSC_DECIDE,M,M,
45: 18,NULL,6,NULL,&Pmat);
46: MatGetOwnershipRange(Amat,&Istart,&Iend);
47: m = Iend-Istart;
48: /* Generate vectors */
49: VecCreate(comm,&xx);
50: VecSetSizes(xx,m,M);
51: VecSetFromOptions(xx);
52: VecDuplicate(xx,&bb);
53: VecSet(bb,.0);
54: /* generate element matrices */
55: {
56: FILE *file;
57: char fname[] = "data/elem_2d_therm.txt";
58: file = fopen(fname, "r");
59: if (file == 0) {
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;
76: } else {
77: for (i=0;i<4;i++) {
78: for (j=0;j<4;j++) {
79: fscanf(file, "%le", &DD1[i][j]);
80: }
81: }
82: }
83: /* BC version of element */
84: for (i=0;i<4;i++) {
85: for (j=0;j<4;j++) {
86: if (i<2 || j < 2) {
87: if (i==j) DD2[i][j] = .1*DD1[i][j];
88: else DD2[i][j] = 0.0;
89: } else DD2[i][j] = DD1[i][j];
90: }
91: }
92: }
93: {
94: PetscReal *coords;
95: PetscMalloc1(2*m,&coords);
96: /* forms the element stiffness for the Laplacian and coordinates */
97: for (Ii=Istart,ix=0; Ii<Iend; Ii++,ix++) {
98: j = Ii/(ne+1); i = Ii%(ne+1);
99: /* coords */
100: x = h*(Ii % (ne+1)); y = h*(Ii/(ne+1));
101: coords[2*ix] = x; coords[2*ix+1] = y;
102: if (i<ne && j<ne) {
103: PetscInt jj,ii,idx[4];
104: /* radius */
105: PetscReal radius = PetscSqrtScalar((x-.5+h/2)*(x-.5+h/2) + (y-.5+h/2)*(y-.5+h/2));
106: PetscReal alpha = 1.0;
108: idx[0] = Ii; idx[1] = Ii+1; idx[2] = Ii + (ne+1) + 1; idx[3] = Ii + (ne+1);
109: 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: MatSetValues(Pmat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
116: if (j>0) {
117: 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: 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: VecSetValues(bb,1,&jj,&v,INSERT_VALUES);
130: }
131: }
132: MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
133: MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
134: MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
135: MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);
136: VecAssemblyBegin(bb);
137: VecAssemblyEnd(bb);
139: /* Setup solver */
140: KSPCreate(PETSC_COMM_WORLD,&ksp);
141: KSPSetType(ksp, KSPCG);
142: KSPGetPC(ksp,&pc);
143: PCSetType(pc,PCGAMG);
144: KSPSetFromOptions(ksp);
146: /* PCGAMGSetType(pc,"agg"); */
148: /* finish KSP/PC setup */
149: KSPSetOperators(ksp, Amat, Amat);
150: PCSetCoordinates(pc, 2, m, coords);
151: PetscFree(coords);
152: }
154: if (!PETSC_TRUE) {
155: PetscViewer viewer;
156: PetscViewerASCIIOpen(comm, "Amat.m", &viewer);
157: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
158: MatView(Amat,viewer);
159: PetscViewerDestroy(&viewer);
160: }
162: /* solve */
163: #if defined(PETSC_USE_LOG)
164: PetscLogStageRegister("Solve", &stage);
165: PetscLogStagePush(stage);
166: #endif
167: VecSet(xx,.0);
169: KSPSetUp(ksp);
171: KSPSolve(ksp,bb,xx);
173: #if defined(PETSC_USE_LOG)
174: PetscLogStagePop();
175: #endif
177: KSPGetIterationNumber(ksp,&its);
179: if (!PETSC_TRUE) {
180: PetscReal norm,norm2;
181: PetscViewer viewer;
182: Vec res;
183: PetscViewerASCIIOpen(comm, "rhs.m", &viewer);
184: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
185: VecView(bb,viewer);
186: PetscViewerDestroy(&viewer);
187: VecNorm(bb, NORM_2, &norm2);
189: PetscViewerASCIIOpen(comm, "solution.m", &viewer);
190: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
191: VecView(xx,viewer);
192: PetscViewerDestroy(&viewer);
194: VecDuplicate(xx, &res);
195: MatMult(Amat, xx, res);
196: VecAXPY(bb, -1.0, res);
197: VecDestroy(&res);
198: VecNorm(bb,NORM_2,&norm);
199: PetscPrintf(PETSC_COMM_WORLD,"[%d]%s |b-Ax|/|b|=%e, |b|=%e\n",0,__FUNCT__,norm/norm2,norm2);
201: PetscViewerASCIIOpen(comm, "residual.m", &viewer);
202: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
203: VecView(bb,viewer);
204: PetscViewerDestroy(&viewer);
205: }
207: /* Free work space */
208: KSPDestroy(&ksp);
209: VecDestroy(&xx);
210: VecDestroy(&bb);
211: MatDestroy(&Amat);
212: MatDestroy(&Pmat);
214: PetscFinalize();
215: return 0;
216: }