Actual source code: ex55.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  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 coeficient in embedded circle\n\n";

  8:  #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:   PetscBool      use_coords = PETSC_FALSE;
 21:   PetscMPIInt    npe,mype;
 22:   PetscScalar    DD[8][8],DD2[8][8];
 23: #if defined(PETSC_USE_LOG)
 24:   PetscLogStage stage[2];
 25: #endif
 26:   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 },
 27:                              {2.0000E-01,  5.333333333333333E-01,  0.0000E-00,  6.666666666666667E-02, -2.0000E-01, -2.666666666666667E-01, 0.0000E-00, -3.333333333333333E-01 },
 28:                              {-3.333333333333333E-01,  0.0000E-00,  5.333333333333333E-01, -2.0000E-01,  6.666666666666667E-02, 0.0000E-00, -2.666666666666667E-01,  2.0000E-01 },
 29:                              {0.0000E+00,  6.666666666666667E-02, -2.0000E-01,  5.333333333333333E-01,  0.0000E-00, -3.333333333333333E-01, 2.0000E-01, -2.666666666666667E-01 },
 30:                              {-2.666666666666667E-01, -2.0000E-01,  6.666666666666667E-02,  0.0000E-00,  5.333333333333333E-01,  2.0000E-01, -3.333333333333333E-01,  0.0000E+00 },
 31:                              {-2.0000E-01, -2.666666666666667E-01, 0.0000E-00, -3.333333333333333E-01,  2.0000E-01,  5.333333333333333E-01, 0.0000E-00,  6.666666666666667E-02 },
 32:                              {6.666666666666667E-02, 0.0000E-00, -2.666666666666667E-01,  2.0000E-01, -3.333333333333333E-01,  0.0000E-00, 5.333333333333333E-01, -2.0000E-01 },
 33:                              {0.0000E-00, -3.333333333333333E-01,  2.0000E-01, -2.666666666666667E-01, 0.0000E-00,  6.666666666666667E-02, -2.0000E-01,  5.333333333333333E-01 } };


 36:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 37:   comm = PETSC_COMM_WORLD;
 38:   MPI_Comm_rank(comm, &mype);
 39:   MPI_Comm_size(comm, &npe);
 40:   PetscOptionsGetInt(NULL,NULL,"-ne",&ne,NULL);
 41:   h     = 1./ne;
 42:   /* ne*ne; number of global elements */
 43:   PetscOptionsGetReal(NULL,NULL,"-alpha",&soft_alpha,NULL);
 44:   PetscOptionsGetBool(NULL,NULL,"-use_coordinates",&use_coords,NULL);
 45:   M    = 2*(ne+1)*(ne+1); /* global number of equations */
 46:   m    = (ne+1)*(ne+1)/npe;
 47:   if (mype==npe-1) m = (ne+1)*(ne+1) - (npe-1)*m;
 48:   m *= 2;
 49:   /* create stiffness matrix */
 50:   MatCreate(comm,&Amat);
 51:   MatSetSizes(Amat,m,m,M,M);
 52:   MatSetBlockSize(Amat,2);
 53:   MatSetType(Amat,MATAIJ);
 54:   MatSetFromOptions(Amat);
 55:   MatSeqAIJSetPreallocation(Amat,18,NULL);
 56:   MatMPIAIJSetPreallocation(Amat,18,NULL,18,NULL);

 58:   MatCreate(comm,&Pmat);
 59:   MatSetSizes(Pmat,m,m,M,M);
 60:   MatSetBlockSize(Pmat,2);
 61:   MatSetType(Pmat,MATAIJ);
 62:   MatSetFromOptions(Pmat);
 63:   MatSeqAIJSetPreallocation(Pmat,18,NULL);
 64:   MatMPIAIJSetPreallocation(Pmat,18,NULL,12,NULL);

 66:   MatGetOwnershipRange(Amat,&Istart,&Iend);
 67:   if (m != Iend - Istart) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"m %D does not equal Iend %D - Istart %D",m,Iend,Istart);
 68:   /* Generate vectors */
 69:   VecCreate(comm,&xx);
 70:   VecSetSizes(xx,m,M);
 71:   VecSetFromOptions(xx);
 72:   VecDuplicate(xx,&bb);
 73:   VecSet(bb,.0);
 74:   /* generate element matrices -- see ex56.c on how to use different data set */
 75:   {
 76:       DD[0][0] =  0.53333333333333321;
 77:       DD[0][1] =  0.20000000000000001;
 78:       DD[0][2] = -0.33333333333333331;
 79:       DD[0][3] =   0.0000000000000000;
 80:       DD[0][4] = -0.26666666666666666;
 81:       DD[0][5] = -0.20000000000000001;
 82:       DD[0][6] =  6.66666666666666796E-002;
 83:       DD[0][7] =  6.93889390390722838E-018;
 84:       DD[1][0] =  0.20000000000000001;
 85:       DD[1][1] =  0.53333333333333333;
 86:       DD[1][2] =  7.80625564189563192E-018;
 87:       DD[1][3] =  6.66666666666666935E-002;
 88:       DD[1][4] = -0.20000000000000001;
 89:       DD[1][5] = -0.26666666666666666;
 90:       DD[1][6] = -3.46944695195361419E-018;
 91:       DD[1][7] = -0.33333333333333331;
 92:       DD[2][0] = -0.33333333333333331;
 93:       DD[2][1] =  1.12757025938492461E-017;
 94:       DD[2][2] =  0.53333333333333333;
 95:       DD[2][3] = -0.20000000000000001;
 96:       DD[2][4] =  6.66666666666666935E-002;
 97:       DD[2][5] = -6.93889390390722838E-018;
 98:       DD[2][6] = -0.26666666666666666;
 99:       DD[2][7] =  0.19999999999999998;
100:       DD[3][0] =   0.0000000000000000;
101:       DD[3][1] =  6.66666666666666935E-002;
102:       DD[3][2] = -0.20000000000000001;
103:       DD[3][3] =  0.53333333333333333;
104:       DD[3][4] =  4.33680868994201774E-018;
105:       DD[3][5] = -0.33333333333333331;
106:       DD[3][6] =  0.20000000000000001;
107:       DD[3][7] = -0.26666666666666666;
108:       DD[4][0] = -0.26666666666666666;
109:       DD[4][1] = -0.20000000000000001;
110:       DD[4][2] =  6.66666666666666935E-002;
111:       DD[4][3] =  8.67361737988403547E-019;
112:       DD[4][4] =  0.53333333333333333;
113:       DD[4][5] =  0.19999999999999998;
114:       DD[4][6] = -0.33333333333333331;
115:       DD[4][7] = -3.46944695195361419E-018;
116:       DD[5][0] = -0.20000000000000001;
117:       DD[5][1] = -0.26666666666666666;
118:       DD[5][2] = -1.04083408558608426E-017;
119:       DD[5][3] = -0.33333333333333331;
120:       DD[5][4] =  0.19999999999999998;
121:       DD[5][5] =  0.53333333333333333;
122:       DD[5][6] =  6.93889390390722838E-018;
123:       DD[5][7] =  6.66666666666666519E-002;
124:       DD[6][0] =  6.66666666666666796E-002;
125:       DD[6][1] = -6.93889390390722838E-018;
126:       DD[6][2] = -0.26666666666666666;
127:       DD[6][3] =  0.19999999999999998;
128:       DD[6][4] = -0.33333333333333331;
129:       DD[6][5] =  6.93889390390722838E-018;
130:       DD[6][6] =  0.53333333333333321;
131:       DD[6][7] = -0.20000000000000001;
132:       DD[7][0] =  6.93889390390722838E-018;
133:       DD[7][1] = -0.33333333333333331;
134:       DD[7][2] =  0.19999999999999998;
135:       DD[7][3] = -0.26666666666666666;
136:       DD[7][4] =   0.0000000000000000;
137:       DD[7][5] =  6.66666666666666519E-002;
138:       DD[7][6] = -0.20000000000000001;
139:       DD[7][7] =  0.53333333333333321;

141:     /* BC version of element */
142:     for (i=0; i<8; i++) {
143:       for (j=0; j<8; j++) {
144:         if (i<4 || j < 4) {
145:           if (i==j) DD2[i][j] = .1*DD1[i][j];
146:           else DD2[i][j] = 0.0;
147:         } else DD2[i][j] = DD1[i][j];
148:       }
149:     }
150:   }
151:   {
152:     PetscReal *coords;
153:     PetscMalloc1(m,&coords);
154:     /* forms the element stiffness and coordinates */
155:     for (Ii = Istart/2, ix = 0; Ii < Iend/2; Ii++, ix++) {
156:       j = Ii/(ne+1); i = Ii%(ne+1);
157:       /* coords */
158:       x            = h*(Ii % (ne+1)); y = h*(Ii/(ne+1));
159:       coords[2*ix] = x; coords[2*ix+1] = y;
160:       if (i<ne && j<ne) {
161:         PetscInt jj,ii,idx[4];
162:         /* radius */
163:         PetscReal radius = PetscSqrtReal((x-.5+h/2)*(x-.5+h/2) + (y-.5+h/2)*(y-.5+h/2));
164:         PetscReal alpha  = 1.0;
165:         if (radius < 0.25) alpha = soft_alpha;

167:         idx[0] = Ii; idx[1] = Ii+1; idx[2] = Ii + (ne+1) + 1;  idx[3] = Ii + (ne+1);
168:         for (ii=0; ii<8; ii++) {
169:           for (jj=0;jj<8;jj++) DD[ii][jj] = alpha*DD1[ii][jj];
170:         }
171:         MatSetValuesBlocked(Pmat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
172:         if (j>0) {
173:           MatSetValuesBlocked(Amat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
174:         } else {
175:           /* a BC */
176:           for (ii=0; ii<8; ii++) {
177:             for (jj=0;jj<8;jj++) DD[ii][jj] = alpha*DD2[ii][jj];
178:           }
179:           MatSetValuesBlocked(Amat,4,idx,4,idx,(const PetscScalar*)DD,ADD_VALUES);
180:         }
181:       }
182:       if (j>0) {
183:         PetscScalar v  = h*h;
184:         PetscInt    jj = 2*Ii; /* load in x direction */
185:         VecSetValues(bb,1,&jj,&v,INSERT_VALUES);
186:       }
187:     }
188:     MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
189:     MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
190:     MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
191:     MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);
192:     VecAssemblyBegin(bb);
193:     VecAssemblyEnd(bb);

195:     /* Setup solver */
196:     KSPCreate(PETSC_COMM_WORLD,&ksp);
197:     KSPSetFromOptions(ksp);

199:     /* finish KSP/PC setup */
200:     KSPSetOperators(ksp, Amat, Amat);
201:     if (use_coords) {
202:       PC             pc;
203:       KSPGetPC(ksp, &pc);
204:       PCSetCoordinates(pc, 2, m/2, coords);
205:     }
206:     PetscFree(coords);
207:   }

209:   if (!PETSC_TRUE) {
210:     PetscViewer viewer;
211:     PetscViewerASCIIOpen(comm, "Amat.m", &viewer);
212:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
213:     MatView(Amat,viewer);
214:     PetscViewerPopFormat(viewer);
215:     PetscViewerDestroy(&viewer);
216:   }

218:   /* solve */
219: #if defined(PETSC_USE_LOG)
220:   PetscLogStageRegister("Setup", &stage[0]);
221:   PetscLogStageRegister("Solve", &stage[1]);
222:   PetscLogStagePush(stage[0]);
223: #endif
224:   KSPSetUp(ksp);
225: #if defined(PETSC_USE_LOG)
226:   PetscLogStagePop();
227: #endif

229:   VecSet(xx,.0);

231: #if defined(PETSC_USE_LOG)
232:   PetscLogStagePush(stage[1]);
233: #endif
234:   KSPSolve(ksp, bb, xx);
235: #if defined(PETSC_USE_LOG)
236:   PetscLogStagePop();
237: #endif

239:   KSPGetIterationNumber(ksp,&its);

241:   if (0) {
242:     PetscReal   norm,norm2;
243:     PetscViewer viewer;
244:     Vec         res;

246:     PetscObjectGetComm((PetscObject)bb,&comm);
247:     VecNorm(bb, NORM_2, &norm2);

249:     VecDuplicate(xx, &res);
250:     MatMult(Amat, xx, res);
251:     VecAXPY(bb, -1.0, res);
252:     VecDestroy(&res);
253:     VecNorm(bb, NORM_2, &norm);
254:     PetscPrintf(PETSC_COMM_WORLD,"[%d]%s |b-Ax|/|b|=%e, |b|=%e\n",0,PETSC_FUNCTION_NAME,norm/norm2,norm2);
255:     PetscViewerASCIIOpen(comm, "residual.m", &viewer);
256:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
257:     VecView(bb,viewer);
258:     PetscViewerPopFormat(viewer);
259:     PetscViewerDestroy(&viewer);
260:   }

262:   /* Free work space */
263:   KSPDestroy(&ksp);
264:   VecDestroy(&xx);
265:   VecDestroy(&bb);
266:   MatDestroy(&Amat);
267:   MatDestroy(&Pmat);

269:   PetscFinalize();
270:   return ierr;
271: }



275: /*TEST

277:    test:
278:       nsize: 4
279:       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_type cg -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
280:       output_file: output/ex55_sa.out

282:    test:
283:       suffix: Classical
284:       nsize: 4
285:       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 -mg_levels_esteig_ksp_type cg
286:       output_file: output/ex55_classical.out

288:    test:
289:       suffix: NC
290:       nsize: 4
291:       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_type cg -pc_gamg_esteig_ksp_max_it 10 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.2

293:    test:
294:       suffix: geo
295:       nsize: 4
296:       args: -ne 29 -alpha 1.e-3 -ksp_type cg -pc_type gamg -pc_gamg_type geo -use_coordinates -ksp_monitor_short -mg_levels_esteig_ksp_type cg -ksp_type cg -ksp_norm_type unpreconditioned  -mg_levels_ksp_max_it 3
297:       output_file: output/ex55_0.out
298:       requires: triangle

300:    test:
301:       suffix: hypre
302:       nsize: 4
303:       requires: hypre !complex
304:       args: -ne 29 -alpha 1.e-3 -ksp_type cg -pc_type hypre -pc_hypre_type boomeramg -ksp_monitor_short

306: TEST*/