Actual source code: ex54.c

petsc-dev 2014-04-18
Report Typos and Errors
  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[2*m];
 95:     /* forms the element stiffness for the Laplacian and coordinates */
 96:     for (Ii=Istart,ix=0; Ii<Iend; Ii++,ix++) {
 97:       j = Ii/(ne+1); i = Ii%(ne+1);
 98:       /* coords */
 99:       x            = h*(Ii % (ne+1)); y = h*(Ii/(ne+1));
100:       coords[2*ix] = x; coords[2*ix+1] = y;
101:       if (i<ne && j<ne) {
102:         PetscInt jj,ii,idx[4] = {Ii, Ii+1, Ii + (ne+1) + 1, Ii + (ne+1)};
103:         /* radius */
104:         PetscReal radius = PetscSqrtScalar((x-.5+h/2)*(x-.5+h/2) + (y-.5+h/2)*(y-.5+h/2));
105:         PetscReal alpha  = 1.0;
106:         if (radius < 0.25) alpha = soft_alpha;


109:         for (ii=0; ii<4; ii++) {
110:           for (jj=0; jj<4; jj++) DD[ii][jj] = alpha*DD1[ii][jj];
111:         }
112:         MatSetValues(Pmat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
113:         if (j>0) {
114:           MatSetValues(Amat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
115:         } else {
116:           /* a BC */
117:           for (ii=0;ii<4;ii++) {
118:             for (jj=0;jj<4;jj++) DD[ii][jj] = alpha*DD2[ii][jj];
119:           }
120:           MatSetValues(Amat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
121:         }
122:       }
123:       if (j>0) {
124:         PetscScalar v  = h*h;
125:         PetscInt    jj = Ii;
126:         VecSetValues(bb,1,&jj,&v,INSERT_VALUES);
127:       }
128:     }
129:     MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
130:     MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
131:     MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
132:     MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);
133:     VecAssemblyBegin(bb);
134:     VecAssemblyEnd(bb);

136:     /* Setup solver */
137:     KSPCreate(PETSC_COMM_WORLD,&ksp);
138:     KSPSetType(ksp, KSPCG);
139:     KSPGetPC(ksp,&pc);
140:     PCSetType(pc,PCGAMG);
141:     KSPSetFromOptions(ksp);

143:     /* PCGAMGSetType(pc,"agg"); */

145:     /* finish KSP/PC setup */
146:     KSPSetOperators(ksp, Amat, Amat);
147:     PCSetCoordinates(pc, 2, m, coords);
148:   }

150:   if (!PETSC_TRUE) {
151:     PetscViewer viewer;
152:     PetscViewerASCIIOpen(comm, "Amat.m", &viewer);
153:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
154:     MatView(Amat,viewer);
155:     PetscViewerDestroy(&viewer);
156:   }

158:   /* solve */
159: #if defined(PETSC_USE_LOG)
160:   PetscLogStageRegister("Solve", &stage);
161:   PetscLogStagePush(stage);
162: #endif
163:   VecSet(xx,.0);

165:   KSPSetUp(ksp);

167:   KSPSolve(ksp,bb,xx);

169: #if defined(PETSC_USE_LOG)
170:   PetscLogStagePop();
171: #endif

173:   KSPGetIterationNumber(ksp,&its);

175:   if (!PETSC_TRUE) {
176:     PetscReal   norm,norm2;
177:     PetscViewer viewer;
178:     Vec         res;
179:     PetscViewerASCIIOpen(comm, "rhs.m", &viewer);
180:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
181:     VecView(bb,viewer);
182:     PetscViewerDestroy(&viewer);
183:     VecNorm(bb, NORM_2, &norm2);

185:     PetscViewerASCIIOpen(comm, "solution.m", &viewer);
186:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
187:     VecView(xx,viewer);
188:     PetscViewerDestroy(&viewer);

190:     VecDuplicate(xx, &res);
191:     MatMult(Amat, xx, res);
192:     VecAXPY(bb, -1.0, res);
193:     VecDestroy(&res);
194:     VecNorm(bb,NORM_2,&norm);
195:     PetscPrintf(PETSC_COMM_WORLD,"[%d]%s |b-Ax|/|b|=%e, |b|=%e\n",0,__FUNCT__,norm/norm2,norm2);

197:     PetscViewerASCIIOpen(comm, "residual.m", &viewer);
198:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
199:     VecView(bb,viewer);
200:     PetscViewerDestroy(&viewer);
201:   }

203:   /* Free work space */
204:   KSPDestroy(&ksp);
205:   VecDestroy(&xx);
206:   VecDestroy(&bb);
207:   MatDestroy(&Amat);
208:   MatDestroy(&Pmat);

210:   PetscFinalize();
211:   return 0;
212: }