Actual source code: ex23.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: -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*,MatStructure*,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,DMDA_BOUNDARY_NONE,DMDA_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:   DMCreateMatrix(user.da,MATAIJ,&user.M);
     73:   DMCreateMatrix(user.da,MATAIJ,&user.S);
     74:   DMCreateMatrix(user.da,MATAIJ,&J);
     75:   /* Form the mass,stiffness matrices and matrix M_0 */
     76:   SetUpMatrices(&user);

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

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

     93:   SetInitialGuess(x,&user);

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

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

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

    108:   TSSetFromOptions(ts);

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

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

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

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

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

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

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

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

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

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

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

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

    200: }

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

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

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

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

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

    224: }

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

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

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

    251:   /* for active set method the matrix does not get changed, so do not need to copy each time,
    252:      if the active set remains the same for several solves the preconditioner does not need to be rebuilt*/
    253:   *flg = SAME_PRECONDITIONER;
    254:   if (!copied) {
    255:     MatCopy(user->S,*J,*flg);
    256:     MatAXPY(*J,a,user->M,*flg);
    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] = -SNES_VI_INF;
    284:       l[j][i][1] = -1.0;
    285:       u[j][i][0] = SNES_VI_INF;
    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: }