Actual source code: ex23.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: static char help[] = "Cahn-Hilliard-2d problem for constant mobility and triangular elements.\n\
  2: Runtime options include:\n\
  3: -xmin <xmin>\n\
  4: -xmax <xmax>\n\
  5: -ymin <ymin>\n\
  6: -gamma <gamma>\n\
  7: -theta_c <theta_c>\n\
  8: -implicit <0,1> treat theta_c*M_0*u explicitly/implicitly\n\n";

 10: /*
 11:     Run with for example: -pc_type mg -pc_mg_galerkin -T .01 -da_grid_x 65 -da_grid_y 65 -pc_mg_levels 4 -ksp_type fgmres -snes_atol 1.e-14 -mat_no_inode
 12:  */

 14: #include <petscts.h>
 15: #include <petscdmda.h>

 17: typedef struct {
 18:   PetscReal   dt,T; /* Time step and end time */
 19:   DM          da;
 20:   Mat         M; /* Mass matrix */
 21:   Mat         S; /* stiffness matrix */
 22:   Mat         M_0;
 23:   Vec         q,u,work1;
 24:   PetscScalar gamma,theta_c; /* physics parameters */
 25:   PetscReal   xmin,xmax,ymin,ymax;
 26:   PetscBool   tsmonitor;
 27:   PetscInt    implicit; /* Evaluate theta_c*Mo*u impliicitly or explicitly */
 28: } AppCtx;

 30: PetscErrorCode GetParams(AppCtx*);
 31: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
 32: PetscErrorCode SetUpMatrices(AppCtx*);
 33: PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
 34: PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
 35: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
 36: PetscErrorCode Update_q(TS);
 37: PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);

 41: int main(int argc, char **argv)
 42: {
 44:   Vec            x,r;  /* Solution and residual vectors */
 45:   TS             ts;   /* Timestepping object */
 46:   AppCtx         user; /* Application context */
 47:   Vec            xl,xu; /* Upper and lower bounds on variables */
 48:   Mat            J;

 50:   PetscInitialize(&argc,&argv, (char*)0, help);

 52:   /* Get physics and time parameters */
 53:   GetParams(&user);
 54:   /* Create a 2D DA with dof = 2 */
 55:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&user.da);
 56:   /* Set Element type (triangular) */
 57:   DMDASetElementType(user.da,DMDA_ELEMENT_P1);

 59:   /* Set x and y coordinates */
 60:   DMDASetUniformCoordinates(user.da,user.xmin,user.xmax,user.ymin,user.ymax,0.0,1.0);

 62:   /* Get global vector x from DM and duplicate vectors r,xl,xu */
 63:   DMCreateGlobalVector(user.da,&x);
 64:   VecDuplicate(x,&r);
 65:   VecDuplicate(x,&xl);
 66:   VecDuplicate(x,&xu);
 67:   if (!user.implicit) {
 68:     VecDuplicate(x,&user.q);
 69:   }

 71:   /* Get mass,stiffness, and jacobian matrix structure from the da */
 72:   DMSetMatType(user.da,MATAIJ);
 73:   DMCreateMatrix(user.da,&user.M);
 74:   DMCreateMatrix(user.da,&user.S);
 75:   DMCreateMatrix(user.da,&J);
 76:   /* Form the mass,stiffness matrices and matrix M_0 */
 77:   SetUpMatrices(&user);

 79:   /* Create timestepping solver context */
 80:   TSCreate(PETSC_COMM_WORLD,&ts);
 81:   TSSetProblemType(ts,TS_NONLINEAR);


 84:   /* Set Function evaluation and jacobian evaluation routines */
 85:   TSSetIFunction(ts,r,FormIFunction,(void*)&user);
 86:   TSSetIJacobian(ts,J,J,FormIJacobian,(void*)&user);
 87:   TSSetType(ts,TSBEULER);
 88:   /*  TSThetaSetTheta(ts,1.0);*/ /* Explicit Euler */
 89:   TSSetDuration(ts,PETSC_DEFAULT,user.T);
 90:   /* Set lower and upper bound vectors */
 91:   SetVariableBounds(user.da,xl,xu);
 92:   TSVISetVariableBounds(ts,xl,xu);

 94:   SetInitialGuess(x,&user);

 96:   if (!user.implicit) {
 97:     TSSetApplicationContext(ts,&user);
 98:     TSSetSolution(ts,x);
 99:     Update_q(ts);
100:     TSSetPostStep(ts,Update_q);
101:   }

103:   if (user.tsmonitor) {
104:     TSMonitorSet(ts,Monitor,&user,NULL);
105:   }

107:   TSSetInitialTimeStep(ts,0.0,user.dt);

109:   TSSetFromOptions(ts);


112:   /* Run the timestepping solver */
113:   TSSolve(ts,x);

115:   VecDestroy(&x);
116:   VecDestroy(&r);
117:   VecDestroy(&xl);
118:   VecDestroy(&xu);
119:   if (!user.implicit) {
120:     VecDestroy(&user.q);
121:     VecDestroy(&user.u);
122:     VecDestroy(&user.work1);
123:     MatDestroy(&user.M_0);
124:   }
125:   MatDestroy(&user.M);
126:   MatDestroy(&user.S);
127:   MatDestroy(&J);
128:   DMDestroy(&user.da);
129:   TSDestroy(&ts);
130:   PetscFinalize();
131:   return 0;
132: }

136: PetscErrorCode Update_q(TS ts)
137: {
139:   AppCtx         *user;
140:   Vec            x;
141:   PetscScalar    *q_arr,*w_arr;
142:   PetscInt       i,n;
143:   PetscScalar    scale;

146:   TSGetApplicationContext(ts,&user);
147:   TSGetSolution(ts,&x);
148:   VecStrideGather(x,1,user->u,INSERT_VALUES);
149:   MatMult(user->M_0,user->u,user->work1);
150:   scale = -(1.0 - user->implicit)*user->theta_c;
151:   VecScale(user->work1,scale);
152:   VecGetLocalSize(user->u,&n);
153:   VecGetArray(user->q,&q_arr);
154:   VecGetArray(user->work1,&w_arr);
155:   for (i=0; i<n; i++) q_arr[2*i+1] = w_arr[i];
156:   VecRestoreArray(user->q,&q_arr);
157:   VecRestoreArray(user->work1,&w_arr);
158:   return(0);
159: }

163: PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
164: {
166:   PetscScalar    *x;
167:   PetscInt       n,i;
168:   PetscRandom    rand;
169:   PetscScalar    value;

172:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
173:   PetscRandomSetFromOptions(rand);

175:   VecGetLocalSize(X,&n);
176:   VecGetArray(X,&x);
177:   /* Set initial guess, only set value for 2nd dof */
178:   for (i=0; i<n/2; i++) {
179:     PetscRandomGetValue(rand,&value);
180:     x[2*i+1] = -0.4 + 0.05*value*(value - 0.5);
181:   }
182:   VecRestoreArray(X,&x);

184:   PetscRandomDestroy(&rand);
185:   return(0);
186: }

188: static void Gausspoints(PetscScalar *xx,PetscScalar *yy,PetscScalar *w,PetscScalar *x,PetscScalar *y)
189: {

191:   xx[0] = 2.0/3.0*x[0] + 1.0/6.0*x[1] + 1.0/6.0*x[2];
192:   xx[1] = 1.0/6.0*x[0] + 2.0/3.0*x[1] + 1.0/6.0*x[2];
193:   xx[2] = 1.0/6.0*x[0] + 1.0/6.0*x[1] + 2.0/3.0*x[2];

195:   yy[0] = 2.0/3.0*y[0] + 1.0/6.0*y[1] + 1.0/6.0*y[2];
196:   yy[1] = 1.0/6.0*y[0] + 2.0/3.0*y[1] + 1.0/6.0*y[2];
197:   yy[2] = 1.0/6.0*y[0] + 1.0/6.0*y[1] + 2.0/3.0*y[2];

199:   *w = PetscAbsScalar(x[0]*(y[2]-y[1]) + x[2]*(y[1]-y[0]) + x[1]*(y[0]-y[2]))/6.0;

201: }

203: static void ShapefunctionsT3(PetscScalar *phi,PetscScalar phider[][2],PetscScalar xx,PetscScalar yy,PetscScalar *x,PetscScalar *y)
204: {
205:   PetscScalar area,a1,a2,a3,b1,b2,b3,c1,c2,c3,pp;

207:   /* Area of the triangle */
208:   area = 1.0/2.0*PetscAbsScalar(x[0]*(y[2]-y[1]) + x[2]*(y[1]-y[0]) + x[1]*(y[0]-y[2]));

210:   a1 = x[1]*y[2]-x[2]*y[1]; a2 = x[2]*y[0]-x[0]*y[2]; a3 = x[0]*y[1]-x[1]*y[0];
211:   b1 = y[1]-y[2]; b2 = y[2]-y[0]; b3 = y[0]-y[1];
212:   c1 = x[2]-x[1]; c2 = x[0]-x[2]; c3 = x[1]-x[0];
213:   pp = 1.0/(2.0*area);

215:   /* shape functions */
216:   phi[0] = pp*(a1 + b1*xx + c1*yy);
217:   phi[1] = pp*(a2 + b2*xx + c2*yy);
218:   phi[2] = pp*(a3 + b3*xx + c3*yy);

220:   /* shape functions derivatives */
221:   phider[0][0] = pp*b1; phider[0][1] = pp*c1;
222:   phider[1][0] = pp*b2; phider[1][1] = pp*c2;
223:   phider[2][0] = pp*b3; phider[2][1] = pp*c3;

225: }

229: PetscErrorCode FormIFunction(TS ts,PetscReal t, Vec X,Vec Xdot,Vec F,void *ctx)
230: {
232:   AppCtx         *user=(AppCtx*)ctx;

235:   MatMult(user->M,Xdot,F);
236:   MatMultAdd(user->S,X,F,F);
237:   if (!user->implicit) {
238:     VecAXPY(F,1.0,user->q);
239:   }
240:   return(0);
241: }

245: PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J,Mat B,void *ctx)
246: {
247:   PetscErrorCode   ierr;
248:   AppCtx           *user  =(AppCtx*)ctx;
249:   static PetscBool copied = PETSC_FALSE;

252:   /* for active set method the matrix does not get changed, so do not need to copy each time,
253:      if the active set remains the same for several solves the preconditioner does not need to be rebuilt*/
254:   if (!copied) {
255:     MatCopy(user->S,J,DIFFERENT_NONZERO_PATTERN);
256:     MatAXPY(J,a,user->M,DIFFERENT_NONZERO_PATTERN);
257:     copied = PETSC_TRUE;
258:   }
259:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
260:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
261:   /*  MatView(*J,0); */
262:   /* SETERRQ(PETSC_COMM_WORLD,1,"Stopped here\n"); */
263:   return(0);
264: }

268: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
269: {
271:   PetscScalar    ***l,***u;
272:   PetscInt       xs,xm,ys,ym;
273:   PetscInt       j,i;

276:   DMDAVecGetArrayDOF(da,xl,&l);
277:   DMDAVecGetArrayDOF(da,xu,&u);

279:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

281:   for (j=ys; j < ys+ym; j++) {
282:     for (i=xs; i < xs+xm; i++) {
283:       l[j][i][0] = -PETSC_INFINITY;
284:       l[j][i][1] = -1.0;
285:       u[j][i][0] = PETSC_INFINITY;
286:       u[j][i][1] = 1.0;
287:     }
288:   }

290:   DMDAVecRestoreArrayDOF(da,xl,&l);
291:   DMDAVecRestoreArrayDOF(da,xu,&u);
292:   return(0);
293: }

297: PetscErrorCode GetParams(AppCtx *user)
298: {
300:   PetscBool      flg;

303:   /* Set default parameters */
304:   user->tsmonitor = PETSC_FALSE;
305:   user->xmin      = 0.0; user->xmax = 1.0;
306:   user->ymin      = 0.0; user->ymax = 1.0;
307:   user->T         = 0.2;    user->dt = 0.0001;
308:   user->gamma     = 3.2E-4; user->theta_c = 1;
309:   user->implicit  = 0;

311:   PetscOptionsGetBool(NULL,"-monitor",&user->tsmonitor,NULL);
312:   PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);
313:   PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);
314:   PetscOptionsGetReal(NULL,"-ymin",&user->ymin,&flg);
315:   PetscOptionsGetReal(NULL,"-ymax",&user->ymax,&flg);
316:   PetscOptionsGetScalar(NULL,"-gamma",&user->gamma,&flg);
317:   PetscOptionsGetScalar(NULL,"-theta_c",&user->theta_c,&flg);
318:   PetscOptionsGetInt(NULL,"-implicit",&user->implicit,&flg);
319:   return(0);
320: }

324: PetscErrorCode SetUpMatrices(AppCtx *user)
325: {
326:   PetscErrorCode    ierr;
327:   PetscInt          nele,nen,i;
328:   const PetscInt    *ele;
329:   Vec               coords;
330:   const PetscScalar *_coords;
331:   PetscScalar       x[3],y[3],xx[3],yy[3],w;
332:   PetscInt          idx[3];
333:   PetscScalar       phi[3],phider[3][2];
334:   PetscScalar       eM_0[3][3],eM_2[3][3];
335:   Mat               M        = user->M;
336:   Mat               S        = user->S;
337:   PetscScalar       gamma    = user->gamma,theta_c=user->theta_c;
338:   PetscInt          implicit = user->implicit;

341:   /* Get ghosted coordinates */
342:   DMGetCoordinatesLocal(user->da,&coords);
343:   VecGetArrayRead(coords,&_coords);

345:   /* Get local element info */
346:   DMDAGetElements(user->da,&nele,&nen,&ele);
347:   for (i=0; i < nele; i++) {
348:     idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
349:     x[0]   = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
350:     x[1]   = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
351:     x[2]   = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];

353:     PetscMemzero(xx,3*sizeof(PetscScalar));
354:     PetscMemzero(yy,3*sizeof(PetscScalar));
355:     Gausspoints(xx,yy,&w,x,y);

357:     eM_0[0][0]=eM_0[0][1]=eM_0[0][2]=0.0;
358:     eM_0[1][0]=eM_0[1][1]=eM_0[1][2]=0.0;
359:     eM_0[2][0]=eM_0[2][1]=eM_0[2][2]=0.0;
360:     eM_2[0][0]=eM_2[0][1]=eM_2[0][2]=0.0;
361:     eM_2[1][0]=eM_2[1][1]=eM_2[1][2]=0.0;
362:     eM_2[2][0]=eM_2[2][1]=eM_2[2][2]=0.0;

364:     PetscInt m;
365:     for (m=0; m<3; m++) {
366:       PetscMemzero(phi,3*sizeof(PetscScalar));
367:       phider[0][0]=phider[0][1]=0.0;
368:       phider[1][0]=phider[1][1]=0.0;
369:       phider[2][0]=phider[2][1]=0.0;

371:       ShapefunctionsT3(phi,phider,xx[m],yy[m],x,y);

373:       PetscInt j,k;
374:       for (j=0;j<3;j++) {
375:         for (k=0;k<3;k++) {
376:           eM_0[k][j] += phi[j]*phi[k]*w;
377:           eM_2[k][j] += phider[j][0]*phider[k][0]*w + phider[j][1]*phider[k][1]*w;
378:         }
379:       }
380:     }
381:     PetscInt    row,cols[6],r;
382:     PetscScalar vals[6];

384:     for (r=0;r<3;r++) {
385:       row = 2*idx[r];
386:       /* Insert values in the mass matrix */
387:       cols[0] = 2*idx[0]+1;     vals[0] = eM_0[r][0];
388:       cols[1] = 2*idx[1]+1;     vals[1] = eM_0[r][1];
389:       cols[2] = 2*idx[2]+1;     vals[2] = eM_0[r][2];

391:       MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);

393:       /* Insert values in the stiffness matrix */
394:       cols[0] = 2*idx[0];   vals[0] = eM_2[r][0];
395:       cols[1] = 2*idx[1];   vals[1] = eM_2[r][1];
396:       cols[2] = 2*idx[2];   vals[2] = eM_2[r][2];

398:       MatSetValuesLocal(S,1,&row,3,cols,vals,ADD_VALUES);

400:       row     = 2*idx[r]+1;
401:       cols[0] = 2*idx[0];     vals[0] = -eM_0[r][0];
402:       cols[1] = 2*idx[0]+1;   vals[1] = gamma*eM_2[r][0]-implicit*theta_c*eM_0[r][0];
403:       cols[2] = 2*idx[1];     vals[2] = -eM_0[r][1];
404:       cols[3] = 2*idx[1]+1;   vals[3] = gamma*eM_2[r][1]-implicit*theta_c*eM_0[r][1];
405:       cols[4] = 2*idx[2];     vals[4] = -eM_0[r][2];
406:       cols[5] = 2*idx[2]+1;   vals[5] = gamma*eM_2[r][2]-implicit*theta_c*eM_2[r][2];

408:       /* Insert values in matrix M for 2nd dof */
409:       MatSetValuesLocal(S,1,&row,6,cols,vals,ADD_VALUES);
410:     }
411:   }

413:   DMDARestoreElements(user->da,&nele,&nen,&ele);
414:   VecRestoreArrayRead(coords,&_coords);

416:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
417:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

419:   MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY);
420:   MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY);

422:   if (!implicit) {
423:     /* Create ISs to extract matrix M_0 from M */
424:     PetscInt n,rstart;
425:     IS       isrow,iscol;

427:     MatGetLocalSize(M,&n,NULL);
428:     MatGetOwnershipRange(M,&rstart,NULL);
429:     ISCreateStride(PETSC_COMM_WORLD,n/2,rstart,2,&isrow);
430:     ISCreateStride(PETSC_COMM_WORLD,n/2,rstart+1,2,&iscol);

432:     /* Extract M_0 from M */
433:     MatGetSubMatrix(M,isrow,iscol,MAT_INITIAL_MATRIX,&user->M_0);

435:     VecCreate(PETSC_COMM_WORLD,&user->u);
436:     VecSetSizes(user->u,n/2,PETSC_DECIDE);
437:     VecSetFromOptions(user->u);
438:     VecDuplicate(user->u,&user->work1);

440:     ISDestroy(&isrow);
441:     ISDestroy(&iscol);
442:   }
443:   return(0);
444: }

448: PetscErrorCode Monitor(TS ts,PetscInt steps,PetscReal time,Vec x,void *mctx)
449: {

453:   PetscPrintf(PETSC_COMM_WORLD,"Solution vector at t = %5.4f\n",time,steps);
454:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
455:   return(0);
456: }