Actual source code: ex54.c
petsc-3.3-p7 2013-05-11
2: static char help[] = "Creates a matrix using simple quadirlateral finite elements, and uses it to test GAMG\n\
3: -ne <size> : problem size\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,bs,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 wcomm;
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: wcomm = PETSC_COMM_WORLD;
34: MPI_Comm_rank( wcomm, &mype );
35: MPI_Comm_size( wcomm, &npe );
36: PetscOptionsGetInt(PETSC_NULL,"-ne",&ne,PETSC_NULL);
37: h = 1./ne;
38: /* ne*ne; number of global elements */
39: PetscOptionsGetReal(PETSC_NULL,"-alpha",&soft_alpha,PETSC_NULL);
40: M = (ne+1)*(ne+1); /* global number of nodes */
41: /* create stiffness matrix */
42: MatCreateAIJ(wcomm,PETSC_DECIDE,PETSC_DECIDE,M,M,
43: 18,PETSC_NULL,6,PETSC_NULL,&Amat);
44: MatCreateAIJ(wcomm,PETSC_DECIDE,PETSC_DECIDE,M,M,
45: 18,PETSC_NULL,6,PETSC_NULL,&Pmat);
46: MatGetOwnershipRange(Amat,&Istart,&Iend);
47: m = Iend-Istart;
48: bs = 1;
49: /* Generate vectors */
50: VecCreate(wcomm,&xx);
51: VecSetSizes(xx,m,M);
52: VecSetFromOptions(xx);
53: VecDuplicate(xx,&bb);
54: VecSet(bb,.0);
55: /* generate element matrices */
56: {
57: FILE *file;
58: char fname[] = "data/elem_2d_therm.txt";
59: file = fopen(fname, "r");
60: if (file == 0) {
61: DD1[0][0] = 0.66666666666666663;
62: DD1[0][1] = -0.16666666666666669;
63: DD1[0][2] = -0.33333333333333343;
64: DD1[0][3] = -0.16666666666666666;
65: DD1[1][0] = -0.16666666666666669;
66: DD1[1][1] = 0.66666666666666663;
67: DD1[1][2] = -0.16666666666666666;
68: DD1[1][3] = -0.33333333333333343;
69: DD1[2][0] = -0.33333333333333343;
70: DD1[2][1] = -0.16666666666666666;
71: DD1[2][2] = 0.66666666666666663;
72: DD1[2][3] = -0.16666666666666663;
73: DD1[3][0] = -0.16666666666666666;
74: DD1[3][1] = -0.33333333333333343;
75: DD1[3][2] = -0.16666666666666663;
76: DD1[3][3] = 0.66666666666666663;
77: }
78: else {
79: for(i=0;i<4;i++)
80: for(j=0;j<4;j++)
81: fscanf(file, "%le", &DD1[i][j]);
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: PetscReal coords[2*m];
93: /* forms the element stiffness for the Laplacian and coordinates */
94: for (Ii=Istart,ix=0; Ii<Iend; Ii++,ix++) {
95: j = Ii/(ne+1); i = Ii%(ne+1);
96: /* coords */
97: x = h*(Ii % (ne+1)); y = h*(Ii/(ne+1));
98: coords[2*ix] = x; coords[2*ix+1] = y;
99: if( i<ne && j<ne ) {
100: PetscInt jj,ii,idx[4] = {Ii, Ii+1, Ii + (ne+1) + 1, Ii + (ne+1)};
101: /* radius */
102: PetscReal radius = PetscSqrtScalar( (x-.5+h/2)*(x-.5+h/2) + (y-.5+h/2)*(y-.5+h/2) );
103: PetscReal alpha = 1.0;
104: if( radius < 0.25 ){
105: alpha = soft_alpha;
106: }
108: for(ii=0;ii<4;ii++)for(jj=0;jj<4;jj++) DD[ii][jj] = alpha*DD1[ii][jj];
109: MatSetValues(Pmat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
110: if( j>0 ) {
111: MatSetValues(Amat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
112: }
113: else {
114: /* a BC */
115: for(ii=0;ii<4;ii++)for(jj=0;jj<4;jj++) DD[ii][jj] = alpha*DD2[ii][jj];
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: KSPSetType( ksp, KSPCG );
135: KSPGetPC(ksp,&pc);
136: PCSetType(pc,PCGAMG);
137: KSPSetFromOptions(ksp);
138:
139: /* finish KSP/PC setup */
140: KSPSetOperators( ksp, Amat, Amat, SAME_NONZERO_PATTERN );
141: PCSetCoordinates( pc, 2, m, coords );
142: }
144: if( !PETSC_TRUE ) {
145: PetscViewer viewer;
146: PetscViewerASCIIOpen(wcomm, "Amat.m", &viewer);
147: PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB);
148: MatView(Amat,viewer);
149: PetscViewerDestroy( &viewer );
150: }
152: /* solve */
153: #if defined(PETSC_USE_LOG)
154: PetscLogStageRegister("Solve", &stage);
155: PetscLogStagePush(stage);
156: #endif
157: VecSet(xx,.0);
158:
159: KSPSolve(ksp,bb,xx);
160:
161: #if defined(PETSC_USE_LOG)
162: PetscLogStagePop();
163: #endif
165: KSPGetIterationNumber(ksp,&its);
167: if( !PETSC_TRUE ) {
168: PetscReal norm,norm2;
169: PetscViewer viewer;
170: Vec res;
171: PetscViewerASCIIOpen(wcomm, "rhs.m", &viewer);
172: PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB);
173: VecView(bb,viewer);
174: PetscViewerDestroy( &viewer );
175: VecNorm( bb, NORM_2, &norm2 );
177: PetscViewerASCIIOpen(wcomm, "solution.m", &viewer);
178: PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB);
179: VecView(xx,viewer);
180: PetscViewerDestroy( &viewer );
182: VecDuplicate( xx, &res );
183: MatMult( Amat, xx, res );
184: VecAXPY( bb, -1.0, res );
185: VecDestroy( &res );
186: VecNorm(bb,NORM_2,&norm);
187: PetscPrintf(PETSC_COMM_WORLD,"[%d]%s |b-Ax|/|b|=%e, |b|=%e\n",0,__FUNCT__,norm/norm2,norm2);
189: PetscViewerASCIIOpen(wcomm, "residual.m", &viewer);
190: PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB);
191: VecView(bb,viewer);
192: PetscViewerDestroy( &viewer );
193: }
195: /* Free work space */
196: KSPDestroy(&ksp);
197: VecDestroy(&xx);
198: VecDestroy(&bb);
199: MatDestroy(&Amat);
200: MatDestroy(&Pmat);
202: PetscFinalize();
203: return 0;
204: }