Actual source code: ex54.c

  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: -T <T>, where <T> is the end time for the time domain simulation\n\
  7: -dt <dt>,where <dt> is the step size for the numerical integration\n\
  8: -gamma <gamma>\n\
  9: -theta_c <theta_c>\n\n";

 11: /*
 12:     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
 13:  */

 15:  #include petscsnes.h
 16:  #include petscdmda.h

 18: typedef struct{
 19:   PetscReal   dt,T; /* Time step and end time */
 20:   DM          da;
 21:   Mat         M;    /* Jacobian 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: }AppCtx;

 29: PetscErrorCode GetParams(AppCtx*);
 30: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
 31: PetscErrorCode SetUpMatrices(AppCtx*);
 32: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 33: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 34: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
 35: PetscErrorCode Update_q(Vec,Vec,Mat,AppCtx*);
 36: PetscLogEvent event_update_q;

 40: int main(int argc, char **argv)
 41: {
 43:   Vec            x,r;  /* Solution and residual vectors */
 44:   SNES           snes; /* Nonlinear solver context */
 45:   AppCtx         user; /* Application context */
 46:   Vec            xl,xu; /* Upper and lower bounds on variables */
 47:   Mat            J;
 48:   PetscReal      t=0.0;
 49:   PETSC_UNUSED PetscLogStage  stage_timestep;
 50:   PetscInt       its;

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

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

 61:   /* Set x and y coordinates */
 62:   DMDASetUniformCoordinates(user.da,user.xmin,user.xmax,user.ymin,user.ymax,0.0,1.0);
 63: 
 64:   /* Get global vector x from DM and duplicate vectors r,xl,xu */
 65:   DMCreateGlobalVector(user.da,&x);
 66:   VecDuplicate(x,&r);
 67:   VecDuplicate(x,&xl);
 68:   VecDuplicate(x,&xu);
 69:   VecDuplicate(x,&user.q);

 71:   /* Get Jacobian matrix structure from the da */
 72:   DMGetMatrix(user.da,MATAIJ,&user.M);
 73:   /* Form the jacobian matrix and M_0 */
 74:   SetUpMatrices(&user);
 75:   MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);

 77:   /* Create nonlinear solver context */
 78:   SNESCreate(PETSC_COMM_WORLD,&snes);
 79:   SNESSetDM(snes,user.da);

 81:   /* Set Function evaluation and jacobian evaluation routines */
 82:   SNESSetFunction(snes,r,FormFunction,(void*)&user);
 83:   SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);

 85:   SNESSetType(snes,SNESVI);
 86:   SNESSetFromOptions(snes);
 87:   /* Set the boundary conditions */
 88:   SetVariableBounds(user.da,xl,xu);
 89:   SNESVISetVariableBounds(snes,xl,xu);

 91:   SetInitialGuess(x,&user);
 92:   PetscLogStageRegister("Time stepping",&stage_timestep);
 93:   PetscLogEventRegister("Update q",MAT_CLASSID,&event_update_q);
 94:   PetscLogStagePush(stage_timestep);
 95:   /* Begin time loop */
 96:   while(t < user.T) {
 97:     Update_q(user.q,user.u,user.M_0,&user);
 98:     SNESSolve(snes,PETSC_NULL,x);
 99:     SNESGetIterationNumber(snes,&its);
100:     if (user.tsmonitor) {
101:       PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4f in %d iterations\n",t,its);
102:     }
103:     VecStrideGather(x,1,user.u,INSERT_VALUES);
104:     t = t + user.dt;
105:   }
106:   PetscLogStagePop();

108:   VecDestroy(&x);
109:   VecDestroy(&r);
110:   VecDestroy(&xl);
111:   VecDestroy(&xu);
112:   VecDestroy(&user.q);
113:   VecDestroy(&user.u);
114:   VecDestroy(&user.work1);
115:   MatDestroy(&user.M);
116:   MatDestroy(&user.M_0);
117:   MatDestroy(&J);
118:   DMDestroy(&user.da);
119:   SNESDestroy(&snes);
120:   PetscFinalize();
121:   return 0;
122: }

126: PetscErrorCode Update_q(Vec q,Vec u,Mat M_0,AppCtx *user)
127: {
129:   PetscScalar    *q_arr,*w_arr;
130:   PetscInt       i,n;
131: 
133:   PetscLogEventBegin(event_update_q,0,0,0,0);
134:   MatMult(M_0,u,user->work1);
135:   VecScale(user->work1,-1.0);
136:   VecGetLocalSize(u,&n);
137:   VecGetArray(q,&q_arr);
138:   VecGetArray(user->work1,&w_arr);
139:   for(i=0;i<n;i++) {
140:     q_arr[2*i]=q_arr[2*i+1] = w_arr[i];
141:   }
142:   VecRestoreArray(q,&q_arr);
143:   VecRestoreArray(user->work1,&w_arr);
144:   PetscLogEventEnd(event_update_q,0,0,0,0);
145:   return(0);
146: }

150: PetscErrorCode SetInitialGuess(Vec X,AppCtx* user)
151: {
153:   PetscScalar    *x,*u;
154:   PetscInt        n,i;
155:   Vec             rand;

158:   /* u = -0.4 + 0.05*rand(N,1)*(rand(N,1) - 0.5) */
159:   VecDuplicate(user->u,&rand);
160:   VecSetRandom(rand,PETSC_NULL);
161:   VecCopy(rand,user->u);
162:   VecShift(rand,-0.5);
163:   VecPointwiseMult(user->u,user->u,rand);
164:   VecDestroy(&rand);
165:   VecScale(user->u,0.05);
166:   VecShift(user->u,-0.4);
167: 
168:   VecGetLocalSize(X,&n);
169:   VecGetArray(X,&x);
170:   VecGetArray(user->u,&u);
171:   /* Set initial guess, only set value for 2nd dof */
172:   for(i=0;i<n/2;i++) {
173:     x[2*i+1] = u[i];
174:   }
175:   VecRestoreArray(X,&x);
176:   VecRestoreArray(user->u,&u);

178:   return(0);
179: }

183: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ctx)
184: {
186:   AppCtx         *user=(AppCtx*)ctx;

189:   MatMultAdd(user->M,X,user->q,F);
190:   return(0);
191: }

195: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
196: {
197:   PetscErrorCode   ierr;
198:   AppCtx           *user=(AppCtx*)ctx;
199:   static PetscBool copied = PETSC_FALSE;

202:   /* for active set method the matrix does not get changed, so do not need to copy each time, 
203:      if the active set remains the same for several solves the preconditioner does not need to be rebuilt*/
204:   *flg = SAME_PRECONDITIONER;
205:   if (!copied) {
206:     MatCopy(user->M,*J,*flg);
207:     copied = PETSC_TRUE;
208:   }
209:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
210:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
211:   return(0);
212: }
215: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
216: {
218:   PetscScalar    ***l,***u;
219:   PetscInt       xs,xm,ys,ym;
220:   PetscInt       j,i;

223:   DMDAVecGetArrayDOF(da,xl,&l);
224:   DMDAVecGetArrayDOF(da,xu,&u);

226:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

228:   for(j=ys; j < ys+ym; j++) {
229:     for(i=xs; i < xs+xm;i++) {
230:       l[j][i][0] = -SNES_VI_INF;
231:       l[j][i][1] = -1.0;
232:       u[j][i][0] = SNES_VI_INF;
233:       u[j][i][1] = 1.0;
234:     }
235:   }
236: 
237:   DMDAVecRestoreArrayDOF(da,xl,&l);
238:   DMDAVecRestoreArrayDOF(da,xu,&u);
239:   return(0);
240: }

244: PetscErrorCode GetParams(AppCtx* user)
245: {
247:   PetscBool      flg;
248: 

251:   /* Set default parameters */
252:   user->tsmonitor = PETSC_FALSE;
253:   user->xmin = 0.0; user->xmax = 1.0;
254:   user->ymin = 0.0; user->ymax = 1.0;
255:   user->T = 0.0002;    user->dt = 0.0001;
256:   user->gamma = 3.2E-4; user->theta_c = 0;

258:   PetscOptionsGetBool(PETSC_NULL,"-ts_monitor",&user->tsmonitor,PETSC_NULL);
259:   PetscOptionsGetReal(PETSC_NULL,"-xmin",&user->xmin,&flg);
260:   PetscOptionsGetReal(PETSC_NULL,"-xmax",&user->xmax,&flg);
261:   PetscOptionsGetReal(PETSC_NULL,"-ymin",&user->ymin,&flg);
262:   PetscOptionsGetReal(PETSC_NULL,"-ymax",&user->ymax,&flg);
263:   PetscOptionsGetReal(PETSC_NULL,"-T",&user->T,&flg);
264:   PetscOptionsGetReal(PETSC_NULL,"-dt",&user->dt,&flg);
265:   PetscOptionsGetScalar(PETSC_NULL,"-gamma",&user->gamma,&flg);
266:   PetscOptionsGetScalar(PETSC_NULL,"-theta_c",&user->theta_c,&flg);

268:   return(0);
269: }
270: 
271: static void Gausspoints(PetscScalar *xx,PetscScalar *yy,PetscScalar *w,PetscScalar *x,PetscScalar *y)
272: {

274:   xx[0] = 2.0/3.0*x[0] + 1.0/6.0*x[1] + 1.0/6.0*x[2];
275:   xx[1] = 1.0/6.0*x[0] + 2.0/3.0*x[1] + 1.0/6.0*x[2];
276:   xx[2] = 1.0/6.0*x[0] + 1.0/6.0*x[1] + 2.0/3.0*x[2];

278:   yy[0] = 2.0/3.0*y[0] + 1.0/6.0*y[1] + 1.0/6.0*y[2];
279:   yy[1] = 1.0/6.0*y[0] + 2.0/3.0*y[1] + 1.0/6.0*y[2];
280:   yy[2] = 1.0/6.0*y[0] + 1.0/6.0*y[1] + 2.0/3.0*y[2];

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

284: }

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

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

293:   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];
294:   b1 = y[1]-y[2]; b2 = y[2]-y[0]; b3 = y[0]-y[1];
295:   c1 = x[2]-x[1]; c2 = x[0]-x[2]; c3 = x[1]-x[0];
296:   pp = 1.0/(2.0*area);

298:   /* shape functions */
299:   phi[0] = pp*(a1 + b1*xx + c1*yy);
300:   phi[1] = pp*(a2 + b2*xx + c2*yy);
301:   phi[2] = pp*(a3 + b3*xx + c3*yy);

303:   /* shape functions derivatives */
304:   phider[0][0] = pp*b1; phider[0][1] = pp*c1;
305:   phider[1][0] = pp*b2; phider[1][1] = pp*c2;
306:   phider[2][0] = pp*b3; phider[2][1] = pp*c3;

308: }

312: PetscErrorCode SetUpMatrices(AppCtx* user)
313: {
314:   PetscErrorCode    ierr;
315:   PetscInt          nele,nen,i;
316:   const PetscInt    *ele;
317:   PetscScalar       dt=user->dt;
318:   Vec               coords;
319:   const PetscScalar *_coords;
320:   PetscScalar       x[3],y[3],xx[3],yy[3],w;
321:   PetscInt          idx[3];
322:   PetscScalar       phi[3],phider[3][2];
323:   PetscScalar       eM_0[3][3],eM_2[3][3];
324:   Mat               M=user->M;
325:   PetscScalar       gamma=user->gamma,theta_c=user->theta_c;
326:   PetscInt          m;
327:   PetscInt          j,k;
328:   PetscInt          row,cols[6],r;
329:   PetscScalar       vals[6];
330:   PetscInt          n,rstart;
331:   IS                isrow,iscol;

334:   /* Get ghosted coordinates */
335:   DMDAGetGhostedCoordinates(user->da,&coords);
336:   VecGetArrayRead(coords,&_coords);

338:   /* Get local element info */
339:   DMDAGetElements(user->da,&nele,&nen,&ele);
340:   for(i=0;i < nele;i++) {
341:     idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
342:     x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
343:     x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
344:     x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
345: 
346:     PetscMemzero(xx,3*sizeof(PetscScalar));
347:     PetscMemzero(yy,3*sizeof(PetscScalar));
348:     Gausspoints(xx,yy,&w,x,y);
349: 
350:     eM_0[0][0]=eM_0[0][1]=eM_0[0][2]=0.0;
351:     eM_0[1][0]=eM_0[1][1]=eM_0[1][2]=0.0;
352:     eM_0[2][0]=eM_0[2][1]=eM_0[2][2]=0.0;
353:     eM_2[0][0]=eM_2[0][1]=eM_2[0][2]=0.0;
354:     eM_2[1][0]=eM_2[1][1]=eM_2[1][2]=0.0;
355:     eM_2[2][0]=eM_2[2][1]=eM_2[2][2]=0.0;


358:     for(m=0;m<3;m++) {
359:       PetscMemzero(phi,3*sizeof(PetscScalar));
360:       phider[0][0]=phider[0][1]=0.0;
361:       phider[1][0]=phider[1][1]=0.0;
362:       phider[2][0]=phider[2][1]=0.0;
363: 
364:       ShapefunctionsT3(phi,phider,xx[m],yy[m],x,y);

366:       for(j=0;j<3;j++) {
367:         for(k=0;k<3;k++) {
368:           eM_0[k][j] += phi[j]*phi[k]*w;
369:           eM_2[k][j] += phider[j][0]*phider[k][0]*w + phider[j][1]*phider[k][1]*w;
370:         }
371:       }
372:     }

374:     for(r=0;r<3;r++) {
375:       row = 2*idx[r];
376:       cols[0] = 2*idx[0];     vals[0] = dt*eM_2[r][0];
377:       cols[1] = 2*idx[0]+1;   vals[1] = eM_0[r][0];
378:       cols[2] = 2*idx[1];     vals[2] = dt*eM_2[r][1];
379:       cols[3] = 2*idx[1]+1;   vals[3] = eM_0[r][1];
380:       cols[4] = 2*idx[2];     vals[4] = dt*eM_2[r][2];
381:       cols[5] = 2*idx[2]+1;   vals[5] = eM_0[r][2];

383:       /* Insert values in matrix M for 1st dof */
384:       MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
385:       row = 2*idx[r]+1;
386:       cols[0] = 2*idx[0];     vals[0] = -eM_0[r][0];
387:       cols[1] = 2*idx[0]+1;   vals[1] = gamma*eM_2[r][0]-theta_c*eM_0[r][0];
388:       cols[2] = 2*idx[1];     vals[2] = -eM_0[r][1];
389:       cols[3] = 2*idx[1]+1;   vals[3] = gamma*eM_2[r][1]-theta_c*eM_0[r][1];
390:       cols[4] = 2*idx[2];     vals[4] = -eM_0[r][2];
391:       cols[5] = 2*idx[2]+1;   vals[5] = gamma*eM_2[r][2]-theta_c*eM_2[r][2];

393:       /* Insert values in matrix M for 2nd dof */
394:       MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
395:     }
396:   }

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

401:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
402:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

404:   /* Create ISs to extract matrix M_0 from M */

406:   MatGetLocalSize(M,&n,PETSC_NULL);
407:   MatGetOwnershipRange(M,&rstart,PETSC_NULL);
408:   ISCreateStride(PETSC_COMM_WORLD,n/2,rstart,2,&isrow);
409:   ISCreateStride(PETSC_COMM_WORLD,n/2,rstart+1,2,&iscol);

411:   /* Extract M_0 from M */
412:   MatGetSubMatrix(M,isrow,iscol,MAT_INITIAL_MATRIX,&user->M_0);
413:   VecCreate(PETSC_COMM_WORLD,&user->u);
414:   VecSetSizes(user->u,n/2,PETSC_DECIDE);
415:   VecSetFromOptions(user->u);
416:   VecDuplicate(user->u,&user->work1);
417:   ISDestroy(&isrow);
418:   ISDestroy(&iscol);

420:   return(0);
421: }