Actual source code: ex54.c
petsc-3.4.2 2013-07-02

      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: {
     42:   PetscErrorCode             ierr;
     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_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,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);

     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:   DMCreateMatrix(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:   /* Set the boundary conditions */
     86:   SetVariableBounds(user.da,xl,xu);
     87:   SNESVISetVariableBounds(snes,xl,xu);
     88:   SNESSetFromOptions(snes);

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

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

    125: PetscErrorCode Update_q(Vec q,Vec u,Mat M_0,AppCtx *user)
    126: {
    128:   PetscScalar    *q_arr,*w_arr;
    129:   PetscInt       i,n;

    132:   PetscLogEventBegin(event_update_q,0,0,0,0);
    133:   MatMult(M_0,u,user->work1);
    134:   VecScale(user->work1,-1.0);
    135:   VecGetLocalSize(u,&n);
    136:   VecGetArray(q,&q_arr);
    137:   VecGetArray(user->work1,&w_arr);
    138:   for (i=0; i<n; i++) q_arr[2*i]=q_arr[2*i+1] = w_arr[i];
    139:   VecRestoreArray(q,&q_arr);
    140:   VecRestoreArray(user->work1,&w_arr);
    141:   PetscLogEventEnd(event_update_q,0,0,0,0);
    142:   return(0);
    143: }

    147: PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
    148: {
    150:   PetscScalar    *x,*u;
    151:   PetscInt       n,i;
    152:   Vec            rand;

    155:   /* u = -0.4 + 0.05*rand(N,1)*(rand(N,1) - 0.5) */
    156:   VecDuplicate(user->u,&rand);
    157:   VecSetRandom(rand,NULL);
    158:   VecCopy(rand,user->u);
    159:   VecShift(rand,-0.5);
    160:   VecPointwiseMult(user->u,user->u,rand);
    161:   VecDestroy(&rand);
    162:   VecScale(user->u,0.05);
    163:   VecShift(user->u,-0.4);

    165:   VecGetLocalSize(X,&n);
    166:   VecGetArray(X,&x);
    167:   VecGetArray(user->u,&u);
    168:   /* Set initial guess, only set value for 2nd dof */
    169:   for (i=0; i<n/2; i++) x[2*i+1] = u[i];
    170:   VecRestoreArray(X,&x);
    171:   VecRestoreArray(user->u,&u);
    172:   return(0);
    173: }

    177: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
    178: {
    180:   AppCtx         *user=(AppCtx*)ctx;

    183:   MatMultAdd(user->M,X,user->q,F);
    184:   return(0);
    185: }

    189: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
    190: {
    191:   PetscErrorCode   ierr;
    192:   AppCtx           *user  =(AppCtx*)ctx;
    193:   static PetscBool copied = PETSC_FALSE;

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

    217:   DMDAVecGetArrayDOF(da,xl,&l);
    218:   DMDAVecGetArrayDOF(da,xu,&u);

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

    222:   for (j=ys; j < ys+ym; j++) {
    223:     for (i=xs; i < xs+xm;i++) {
    224:       l[j][i][0] = -SNES_VI_INF;
    225:       l[j][i][1] = -1.0;
    226:       u[j][i][0] = SNES_VI_INF;
    227:       u[j][i][1] = 1.0;
    228:     }
    229:   }

    231:   DMDAVecRestoreArrayDOF(da,xl,&l);
    232:   DMDAVecRestoreArrayDOF(da,xu,&u);
    233:   return(0);
    234: }

    238: PetscErrorCode GetParams(AppCtx *user)
    239: {
    241:   PetscBool      flg;

    244:   /* Set default parameters */
    245:   user->tsmonitor = PETSC_FALSE;
    246:   user->xmin      = 0.0;    user->xmax    = 1.0;
    247:   user->ymin      = 0.0;    user->ymax    = 1.0;
    248:   user->T         = 0.0002; user->dt      = 0.0001;
    249:   user->gamma     = 3.2E-4; user->theta_c = 0;

    251:   PetscOptionsGetBool(NULL,"-ts_monitor",&user->tsmonitor,NULL);
    252:   PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);
    253:   PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);
    254:   PetscOptionsGetReal(NULL,"-ymin",&user->ymin,&flg);
    255:   PetscOptionsGetReal(NULL,"-ymax",&user->ymax,&flg);
    256:   PetscOptionsGetReal(NULL,"-T",&user->T,&flg);
    257:   PetscOptionsGetReal(NULL,"-dt",&user->dt,&flg);
    258:   PetscOptionsGetScalar(NULL,"-gamma",&user->gamma,&flg);
    259:   PetscOptionsGetScalar(NULL,"-theta_c",&user->theta_c,&flg);
    260:   return(0);
    261: }

    263: static void Gausspoints(PetscScalar *xx,PetscScalar *yy,PetscScalar *w,PetscScalar *x,PetscScalar *y)
    264: {

    266:   xx[0] = 2.0/3.0*x[0] + 1.0/6.0*x[1] + 1.0/6.0*x[2];
    267:   xx[1] = 1.0/6.0*x[0] + 2.0/3.0*x[1] + 1.0/6.0*x[2];
    268:   xx[2] = 1.0/6.0*x[0] + 1.0/6.0*x[1] + 2.0/3.0*x[2];

    270:   yy[0] = 2.0/3.0*y[0] + 1.0/6.0*y[1] + 1.0/6.0*y[2];
    271:   yy[1] = 1.0/6.0*y[0] + 2.0/3.0*y[1] + 1.0/6.0*y[2];
    272:   yy[2] = 1.0/6.0*y[0] + 1.0/6.0*y[1] + 2.0/3.0*y[2];

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

    276: }

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

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

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

    290:   /* shape functions */
    291:   phi[0] = pp*(a1 + b1*xx + c1*yy);
    292:   phi[1] = pp*(a2 + b2*xx + c2*yy);
    293:   phi[2] = pp*(a3 + b3*xx + c3*yy);

    295:   /* shape functions derivatives */
    296:   phider[0][0] = pp*b1; phider[0][1] = pp*c1;
    297:   phider[1][0] = pp*b2; phider[1][1] = pp*c2;
    298:   phider[2][0] = pp*b3; phider[2][1] = pp*c3;

    300: }

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

    326:   /* Get ghosted coordinates */
    327:   DMGetCoordinatesLocal(user->da,&coords);
    328:   VecGetArrayRead(coords,&_coords);

    330:   /* Get local element info */
    331:   DMDAGetElements(user->da,&nele,&nen,&ele);
    332:   for (i=0; i < nele; i++) {
    333:     idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
    334:     x[0]   = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
    335:     x[1]   = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
    336:     x[2]   = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];

    338:     PetscMemzero(xx,3*sizeof(PetscScalar));
    339:     PetscMemzero(yy,3*sizeof(PetscScalar));
    340:     Gausspoints(xx,yy,&w,x,y);

    342:     eM_0[0][0]=eM_0[0][1]=eM_0[0][2]=0.0;
    343:     eM_0[1][0]=eM_0[1][1]=eM_0[1][2]=0.0;
    344:     eM_0[2][0]=eM_0[2][1]=eM_0[2][2]=0.0;
    345:     eM_2[0][0]=eM_2[0][1]=eM_2[0][2]=0.0;
    346:     eM_2[1][0]=eM_2[1][1]=eM_2[1][2]=0.0;
    347:     eM_2[2][0]=eM_2[2][1]=eM_2[2][2]=0.0;

    350:     for (m=0; m<3; m++) {
    351:       PetscMemzero(phi,3*sizeof(PetscScalar));
    352:       phider[0][0]=phider[0][1]=0.0;
    353:       phider[1][0]=phider[1][1]=0.0;
    354:       phider[2][0]=phider[2][1]=0.0;

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

    358:       for (j=0; j<3; j++) {
    359:         for (k=0; k<3; k++) {
    360:           eM_0[k][j] += phi[j]*phi[k]*w;
    361:           eM_2[k][j] += phider[j][0]*phider[k][0]*w + phider[j][1]*phider[k][1]*w;
    362:         }
    363:       }
    364:     }

    366:     for (r=0; r<3; r++) {
    367:       row     = 2*idx[r];
    368:       cols[0] = 2*idx[0];     vals[0] = dt*eM_2[r][0];
    369:       cols[1] = 2*idx[0]+1;   vals[1] = eM_0[r][0];
    370:       cols[2] = 2*idx[1];     vals[2] = dt*eM_2[r][1];
    371:       cols[3] = 2*idx[1]+1;   vals[3] = eM_0[r][1];
    372:       cols[4] = 2*idx[2];     vals[4] = dt*eM_2[r][2];
    373:       cols[5] = 2*idx[2]+1;   vals[5] = eM_0[r][2];

    375:       /* Insert values in matrix M for 1st dof */
    376:       MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
    377:       row     = 2*idx[r]+1;
    378:       cols[0] = 2*idx[0];     vals[0] = -eM_0[r][0];
    379:       cols[1] = 2*idx[0]+1;   vals[1] = gamma*eM_2[r][0]-theta_c*eM_0[r][0];
    380:       cols[2] = 2*idx[1];     vals[2] = -eM_0[r][1];
    381:       cols[3] = 2*idx[1]+1;   vals[3] = gamma*eM_2[r][1]-theta_c*eM_0[r][1];
    382:       cols[4] = 2*idx[2];     vals[4] = -eM_0[r][2];
    383:       cols[5] = 2*idx[2]+1;   vals[5] = gamma*eM_2[r][2]-theta_c*eM_2[r][2];

    385:       /* Insert values in matrix M for 2nd dof */
    386:       MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
    387:     }
    388:   }

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

    393:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
    394:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

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

    398:   MatGetLocalSize(M,&n,NULL);
    399:   MatGetOwnershipRange(M,&rstart,NULL);
    400:   ISCreateStride(PETSC_COMM_WORLD,n/2,rstart,2,&isrow);
    401:   ISCreateStride(PETSC_COMM_WORLD,n/2,rstart+1,2,&iscol);

    403:   /* Extract M_0 from M */
    404:   MatGetSubMatrix(M,isrow,iscol,MAT_INITIAL_MATRIX,&user->M_0);
    405:   VecCreate(PETSC_COMM_WORLD,&user->u);
    406:   VecSetSizes(user->u,n/2,PETSC_DECIDE);
    407:   VecSetFromOptions(user->u);
    408:   VecDuplicate(user->u,&user->work1);
    409:   ISDestroy(&isrow);
    410:   ISDestroy(&iscol);
    411:   return(0);
    412: }