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

      1: static char help[] = "1D coupled Allen-Cahn and Cahn-Hilliard equation for degenerate 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: ./ex63 -ksp_type fgmres -snes_vi_monitor   -snes_atol 1.e-11 -snes_converged_reason -ksp_converged_reason   -snes_linesearch_monitor -VG 10000000 -draw_fields 1,3,4  -pc_type mg -pc_mg_galerkin -log_summary -dt .000001 -mg_coarse_pc_type svd  -ksp_monitor_true_residual -ksp_rtol 1.e-9 -snes_linesearch basic -T .0020 -P_casc .0005 -snes_monitor_solution -da_refine 10
     13:  */

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

     18: typedef struct {
     19:   PetscReal   dt,T; /* Time step and end time */
     20:   DM          da1,da1_clone,da2;
     21:   Mat         M;    /* Jacobian matrix */
     22:   Mat         M_0;
     23:   Vec         q,wv,cv,wi,ci,eta,cvi,DPsiv,DPsii,DPsieta,Pv,Pi,Piv,logcv,logci,logcvi,Riv;
     24:   Vec         work1,work2,work3,work4;
     25:   PetscScalar Dv,Di,Evf,Eif,A,kBT,kav,kai,kaeta,Rsurf,Rbulk,L,P_casc,VG; /* physics parameters */
     26:   PetscReal   xmin,xmax,ymin,ymax;
     27:   PetscInt    nx;
     28:   PetscBool   voidgrowth;
     29:   PetscBool   degenerate;
     30:   PetscBool   graphics;
     31:   PetscReal   smallnumber;
     32:   PetscReal   initv;
     33:   PetscReal   initeta;
     34: } AppCtx;

     36: PetscErrorCode GetParams(AppCtx*);
     37: PetscErrorCode SetRandomVectors(AppCtx*);
     38: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
     39: PetscErrorCode SetUpMatrices(AppCtx*);
     40: PetscErrorCode UpdateMatrices(AppCtx*);
     41: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
     42: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
     43: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
     44: PetscErrorCode Update_q(AppCtx*);
     45: PetscErrorCode Update_u(Vec,AppCtx*);
     46: PetscErrorCode DPsi(AppCtx*);
     47: PetscErrorCode Llog(Vec,Vec);
     48: PetscErrorCode CheckRedundancy(SNES,IS,IS*,DM);

     52: int main(int argc, char **argv)
     53: {
     55:   Vec            x,r;  /* Solution and residual vectors */
     56:   SNES           snes; /* Nonlinear solver context */
     57:   AppCtx         user; /* Application context */
     58:   Vec            xl,xu; /* Upper and lower bounds on variables */
     59:   Mat            J;
     60:   PetscScalar    t=0.0;
     61:   PetscViewer    view_out, view_p, view_q, view_psi, view_mat;
     62:   PetscReal      bounds[] = {1000.0,-1000.,0.0,1.0,1000.0,-1000.0,0.0,1.0,1000.0,-1000.0};

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

     67:   /* Get physics and time parameters */
     68:   GetParams(&user);
     69:   /* Create a 1D DA with dof = 5; the whole thing */
     70:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC, -4, 5, 1,NULL,&user.da1);
     71:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC, -4, 5, 1,NULL,&user.da1_clone);
     72:   /* Create a 1D DA with dof = 1; for individual componentes */
     73:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC, -4, 1, 1,NULL,&user.da2);

     75:   /* Set Element type (triangular) */
     76:   DMDASetElementType(user.da1,DMDA_ELEMENT_P1);
     77:   DMDASetElementType(user.da2,DMDA_ELEMENT_P1);

     79:   /* Set x and y coordinates */
     80:   DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,NULL,NULL,NULL,NULL);
     81:   DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,NULL,NULL,NULL,NULL);
     82:   /* Get global vector x from DM (da1) and duplicate vectors r,xl,xu */
     83:   DMCreateGlobalVector(user.da1,&x);
     84:   VecDuplicate(x,&r);
     85:   VecDuplicate(x,&xl);
     86:   VecDuplicate(x,&xu);
     87:   VecDuplicate(x,&user.q);

     89:   /* Get global vector user->wv from da2 and duplicate other vectors */
     90:   DMCreateGlobalVector(user.da2,&user.wv);
     91:   VecDuplicate(user.wv,&user.cv);
     92:   VecDuplicate(user.wv,&user.wi);
     93:   VecDuplicate(user.wv,&user.ci);
     94:   VecDuplicate(user.wv,&user.eta);
     95:   VecDuplicate(user.wv,&user.cvi);
     96:   VecDuplicate(user.wv,&user.DPsiv);
     97:   VecDuplicate(user.wv,&user.DPsii);
     98:   VecDuplicate(user.wv,&user.DPsieta);
     99:   VecDuplicate(user.wv,&user.Pv);
    100:   VecDuplicate(user.wv,&user.Pi);
    101:   VecDuplicate(user.wv,&user.Piv);
    102:   VecDuplicate(user.wv,&user.Riv);
    103:   VecDuplicate(user.wv,&user.logcv);
    104:   VecDuplicate(user.wv,&user.logci);
    105:   VecDuplicate(user.wv,&user.logcvi);
    106:   VecDuplicate(user.wv,&user.work1);
    107:   VecDuplicate(user.wv,&user.work2);
    108:   VecDuplicate(user.wv,&user.work3);
    109:   VecDuplicate(user.wv,&user.work4);

    111:   /* Get Jacobian matrix structure from the da for the entire thing, da1 */
    112:   DMCreateMatrix(user.da1,MATAIJ,&user.M);
    113:   /* Get the (usual) mass matrix structure from da2 */
    114:   DMCreateMatrix(user.da2,MATAIJ,&user.M_0);
    115:   SetInitialGuess(x,&user);
    116:   /* Form the jacobian matrix and M_0 */
    117:   SetUpMatrices(&user);
    118:   MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);

    120:   SNESCreate(PETSC_COMM_WORLD,&snes);
    121:   SNESSetDM(snes,user.da1);

    123:   SNESSetFunction(snes,r,FormFunction,(void*)&user);
    124:   SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);

    126:   SetVariableBounds(user.da1,xl,xu);
    127:   SNESVISetVariableBounds(snes,xl,xu);
    128:   SNESSetFromOptions(snes);
    129:   /* SNESVISetRedundancyCheck(snes,(PetscErrorCode (*)(SNES,IS,IS*,void*))CheckRedundancy,user.da1_clone); */
    130:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
    131:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_p",FILE_MODE_WRITE,&view_p);
    132:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
    133:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);
    134:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);
    135:   /* PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds); */

    137:   if (user.graphics) {
    138:     PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds);

    140:     VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
    141:   }
    142:   while (t<user.T) {
    143:     char        filename[PETSC_MAX_PATH_LEN];
    144:     PetscScalar a = 1.0;
    145:     PetscInt    i;
    146:     PetscViewer view;

    149:     SNESSetFunction(snes,r,FormFunction,(void*)&user);
    150:     SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);

    152:     SetRandomVectors(&user);
    153:     DPsi(&user);
    154:     VecView(user.DPsiv,view_psi);
    155:     VecView(user.DPsii,view_psi);
    156:     VecView(user.DPsieta,view_psi);

    158:     VecView(user.Pv,view_p);
    159:     Update_q(&user);
    160:     VecView(user.q,view_q);
    161:     MatView(user.M,view_mat);
    162:     SNESSolve(snes,NULL,x);

    164:     VecView(x,view_out);
    165:     if (user.graphics) {
    166:       VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
    167:     }

    169:     PetscInt its;
    170:     SNESGetIterationNumber(snes,&its);
    171:     PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4g in %d iterations\n",t,its);

    173:     Update_u(x,&user);
    174:     for (i=0; i < (int)(user.T/a) ; i++) {
    175:       if (t/a > i - user.dt/a && t/a < i + user.dt/a) {
    176:         sprintf(filename,"output_%f",t);
    177:         PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&view);
    178:         VecView(user.cv,view);
    179:         VecView(user.ci,view);
    180:         VecView(user.eta,view);
    181:         PetscViewerDestroy(&view);
    182:       }
    183:     }
    184:     UpdateMatrices(&user);
    185:     t    = t + user.dt;

    187:   }

    190:   PetscViewerDestroy(&view_out);
    191:   PetscViewerDestroy(&view_p);
    192:   PetscViewerDestroy(&view_q);
    193:   PetscViewerDestroy(&view_psi);
    194:   PetscViewerDestroy(&view_mat);
    195:   VecDestroy(&x);
    196:   VecDestroy(&r);
    197:   VecDestroy(&xl);
    198:   VecDestroy(&xu);
    199:   VecDestroy(&user.q);
    200:   VecDestroy(&user.wv);
    201:   VecDestroy(&user.cv);
    202:   VecDestroy(&user.wi);
    203:   VecDestroy(&user.ci);
    204:   VecDestroy(&user.eta);
    205:   VecDestroy(&user.cvi);
    206:   VecDestroy(&user.DPsiv);
    207:   VecDestroy(&user.DPsii);
    208:   VecDestroy(&user.DPsieta);
    209:   VecDestroy(&user.Pv);
    210:   VecDestroy(&user.Pi);
    211:   VecDestroy(&user.Piv);
    212:   VecDestroy(&user.Riv);
    213:   VecDestroy(&user.logcv);
    214:   VecDestroy(&user.logci);
    215:   VecDestroy(&user.logcvi);
    216:   VecDestroy(&user.work1);
    217:   VecDestroy(&user.work2);
    218:   VecDestroy(&user.work3);
    219:   VecDestroy(&user.work4);
    220:   MatDestroy(&user.M);
    221:   MatDestroy(&user.M_0);
    222:   DMDestroy(&user.da1);
    223:   DMDestroy(&user.da1_clone);
    224:   DMDestroy(&user.da2);
    225:   SNESDestroy(&snes);
    226:   PetscFinalize();
    227:   return 0;
    228: }

    232: PetscErrorCode Update_u(Vec X,AppCtx *user)
    233: {
    235:   PetscInt       i,n;
    236:   PetscScalar    *xx,*wv_p,*cv_p,*wi_p,*ci_p,*eta_p;
    237:   PetscInt       pv1,pi1,peta1,pv2,pi2,peta2;
    238:   PetscReal      maxv,maxi,maxeta,minv,mini,mineta;

    242:   VecGetLocalSize(user->wv,&n);

    245:   VecMax(user->cv,&pv1,&maxv);
    246:   VecMax(user->ci,&pi1,&maxi);
    247:   VecMax(user->eta,&peta1,&maxeta);
    248:   VecMin(user->cv,&pv2,&minv);
    249:   VecMin(user->ci,&pi2,&mini);
    250:   VecMin(user->eta,&peta2,&mineta);

    253:   VecGetArray(X,&xx);
    254:   VecGetArray(user->wv,&wv_p);
    255:   VecGetArray(user->cv,&cv_p);
    256:   VecGetArray(user->wi,&wi_p);
    257:   VecGetArray(user->ci,&ci_p);
    258:   VecGetArray(user->eta,&eta_p);

    261:   for (i=0; i<n; i++) {
    262:     wv_p[i]  = xx[5*i];
    263:     cv_p[i]  = xx[5*i+1];
    264:     wi_p[i]  = xx[5*i+2];
    265:     ci_p[i]  = xx[5*i+3];
    266:     eta_p[i] = xx[5*i+4];
    267:   }
    268:   VecRestoreArray(X,&xx);
    269:   VecRestoreArray(user->wv,&wv_p);
    270:   VecRestoreArray(user->cv,&cv_p);
    271:   VecRestoreArray(user->wi,&wi_p);
    272:   VecRestoreArray(user->ci,&ci_p);
    273:   VecRestoreArray(user->eta,&eta_p);
    274:   return(0);
    275: }

    279: PetscErrorCode Update_q(AppCtx *user)
    280: {
    282:   PetscScalar    *q_p, *w1, *w2;
    283:   PetscInt       i,n;
    284:   PetscScalar    norm1;

    287:   VecPointwiseMult(user->Riv,user->eta,user->eta);
    288:   VecScale(user->Riv,user->Rsurf);
    289:   VecShift(user->Riv,user->Rbulk);
    290:   VecPointwiseMult(user->Riv,user->ci,user->Riv);
    291:   VecPointwiseMult(user->Riv,user->cv,user->Riv);

    293:   VecCopy(user->Riv,user->work1);
    294:   VecScale(user->work1,user->dt);
    295:   VecAXPY(user->work1,-1.0,user->cv);
    296:   MatMult(user->M_0,user->work1,user->work2);

    298:   VecGetArray(user->q,&q_p);
    299:   VecGetArray(user->work1,&w1);
    300:   VecGetArray(user->work2,&w2);
    301:   VecGetLocalSize(user->wv,&n);
    302:   for (i=0; i<n; i++) q_p[5*i]=w2[i];

    304:   MatMult(user->M_0,user->DPsiv,user->work1);
    305:   for (i=0; i<n; i++) q_p[5*i+1]=w1[i];

    307:   VecCopy(user->Riv,user->work1);
    308:   VecScale(user->work1,user->dt);
    309:   VecAXPY(user->work1,-1.0,user->ci);
    310:   MatMult(user->M_0,user->work1,user->work2);
    311:   for (i=0; i<n; i++) q_p[5*i+2]=w2[i];

    313:   MatMult(user->M_0,user->DPsii,user->work1);
    314:   for (i=0; i<n; i++) q_p[5*i+3]=w1[i];

    316:   VecCopy(user->DPsieta,user->work1);
    317:   VecScale(user->work1,user->L);
    318:   VecAXPY(user->work1,-1.0/user->dt,user->eta);
    319:   MatMult(user->M_0,user->work1,user->work2);
    320:   for (i=0; i<n; i++) q_p[5*i+4]=w2[i];

    322:   VecRestoreArray(user->q,&q_p);
    323:   VecRestoreArray(user->work1,&w1);
    324:   VecRestoreArray(user->work2,&w2);
    325:   return(0);
    326: }

    330: PetscErrorCode DPsi(AppCtx *user)
    331: {
    333:   PetscScalar    Evf=user->Evf,Eif=user->Eif,kBT=user->kBT,A=user->A;
    334:   PetscScalar    *cv_p,*ci_p,*eta_p,*logcv_p,*logci_p,*logcvi_p,*DPsiv_p,*DPsii_p,*DPsieta_p;
    335:   PetscInt       n,i;

    338:   VecGetLocalSize(user->cv,&n);
    339:   VecGetArray(user->cv,&cv_p);
    340:   VecGetArray(user->ci,&ci_p);
    341:   VecGetArray(user->eta,&eta_p);
    342:   VecGetArray(user->logcv,&logcv_p);
    343:   VecGetArray(user->logci,&logci_p);
    344:   VecGetArray(user->logcvi,&logcvi_p);
    345:   VecGetArray(user->DPsiv,&DPsiv_p);
    346:   VecGetArray(user->DPsii,&DPsii_p);
    347:   VecGetArray(user->DPsieta,&DPsieta_p);

    349:   Llog(user->cv,user->logcv);
    350:   Llog(user->ci,user->logci);

    352:   VecCopy(user->cv,user->cvi);
    353:   VecAXPY(user->cvi,1.0,user->ci);
    354:   VecScale(user->cvi,-1.0);
    355:   VecShift(user->cvi,1.0);
    356:   Llog(user->cvi,user->logcvi);

    358:   for (i=0; i<n; i++)
    359:   {
    360:     DPsiv_p[i] = (eta_p[i]-1.0)*(eta_p[i]-1.0)*(Evf + kBT*(logcv_p[i] - logcvi_p[i])) + eta_p[i]*eta_p[i]*2*A*(cv_p[i]-1);
    361:     DPsii_p[i] = (eta_p[i]-1.0)*(eta_p[i]-1.0)*(Eif + kBT*(logci_p[i] - logcvi_p[i])) + eta_p[i]*eta_p[i]*2*A*ci_p[i];

    363:     DPsieta_p[i] = 2.0*(eta_p[i]-1.0)*(Evf*cv_p[i] + Eif*ci_p[i] + kBT*(cv_p[i]* logcv_p[i] + ci_p[i]* logci_p[i] + (1-cv_p[i]-ci_p[i])*logcvi_p[i])) + 2.0*eta_p[i]*A*((cv_p[i]-1.0)*(cv_p[i]-1.0) + ci_p[i]*ci_p[i]);
    364:   }

    366:   VecRestoreArray(user->cv,&cv_p);
    367:   VecRestoreArray(user->ci,&ci_p);
    368:   VecRestoreArray(user->eta,&eta_p);
    369:   VecGetArray(user->logcv,&logcv_p);
    370:   VecGetArray(user->logci,&logci_p);
    371:   VecGetArray(user->logcvi,&logcvi_p);
    372:   VecRestoreArray(user->DPsiv,&DPsiv_p);
    373:   VecRestoreArray(user->DPsii,&DPsii_p);
    374:   VecRestoreArray(user->DPsieta,&DPsieta_p);
    375:   return(0);

    377: }

    382: PetscErrorCode Llog(Vec X, Vec Y)
    383: {
    385:   PetscScalar    *x,*y;
    386:   PetscInt       n,i;

    389:   VecGetArray(X,&x);
    390:   VecGetArray(Y,&y);
    391:   VecGetLocalSize(X,&n);
    392:   for (i=0; i<n; i++) {
    393:     if (x[i] < 1.0e-12) y[i] = log(1.0e-12);
    394:     else y[i] = log(x[i]);
    395:   }
    396:   return(0);
    397: }

    401: PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
    402: {

    405:   PetscInt          n,i,Mda;
    406:   PetscScalar       *xx,*cv_p,*ci_p,*wv_p,*wi_p,*eta_p;
    407:   PetscViewer       view_out;
    408:   /* needed for the void growth case */
    409:   PetscScalar       xmid,cv_v=1.0,cv_m=0.122,ci_v=0.0,ci_m=.00069,eta_v=1.0,eta_m=0.0,h,lambda;
    410:   PetscInt          nele,nen,idx[2];
    411:   const PetscInt    *ele;
    412:   PetscScalar       x[2];
    413:   Vec               coords;
    414:   const PetscScalar *_coords;
    415:   PetscViewer       view;
    416:   PetscScalar       xwidth = user->xmax - user->xmin;

    419:   VecGetLocalSize(X,&n);

    421:   if (user->voidgrowth) {
    422:     DMDAGetInfo(user->da2,NULL,&Mda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
    423:     DMGetCoordinatesLocal(user->da2,&coords);
    424:     VecGetArrayRead(coords,&_coords);

    426:     h      = (user->xmax-user->xmin)/Mda;
    427:     xmid   = (user->xmax + user->xmin)/2.0;
    428:     lambda = 4.0*h;

    430:     DMDAGetElements(user->da2,&nele,&nen,&ele);
    431:     for (i=0; i < nele; i++) {
    432:       idx[0] = ele[2*i]; idx[1] = ele[2*i+1];

    434:       x[0] = _coords[idx[0]];
    435:       x[1] = _coords[idx[1]];

    438:       PetscInt    k;
    439:       PetscScalar vals_cv[2],vals_ci[2],vals_eta[2],s,hhr,r;
    440:       for (k=0; k < 2; k++) {
    441:         s = PetscAbs(x[k] - xmid);
    442:         if (s < xwidth*(5.0/64.0)) {
    443:           vals_cv[k]  = cv_v;
    444:           vals_ci[k]  = ci_v;
    445:           vals_eta[k] = eta_v;
    446:         } else if (s>= xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0)) {
    447:           /* r = (s - xwidth*(6.0/64.0))/(0.5*lambda); */
    448:           r           = (s - xwidth*(6.0/64.0))/(xwidth/64.0);
    449:           hhr         = 0.25*(-r*r*r + 3*r + 2);
    450:           vals_cv[k]  = cv_m + (1.0 - hhr)*(cv_v - cv_m);
    451:           vals_ci[k]  = ci_m + (1.0 - hhr)*(ci_v - ci_m);
    452:           vals_eta[k] = eta_m + (1.0 - hhr)*(eta_v - eta_m);
    453:         } else {
    454:           vals_cv[k]  = cv_m;
    455:           vals_ci[k]  = ci_m;
    456:           vals_eta[k] = eta_m;
    457:         }
    458:       }

    460:       VecSetValuesLocal(user->cv,2,idx,vals_cv,INSERT_VALUES);
    461:       VecSetValuesLocal(user->ci,2,idx,vals_ci,INSERT_VALUES);
    462:       VecSetValuesLocal(user->eta,2,idx,vals_eta,INSERT_VALUES);
    463:     }
    464:     VecAssemblyBegin(user->cv);
    465:     VecAssemblyEnd(user->cv);
    466:     VecAssemblyBegin(user->ci);
    467:     VecAssemblyEnd(user->ci);
    468:     VecAssemblyBegin(user->eta);
    469:     VecAssemblyEnd(user->eta);

    471:     DMDARestoreElements(user->da2,&nele,&nen,&ele);
    472:     VecRestoreArrayRead(coords,&_coords);

    475:   } else {

    477:     VecGetArray(user->cv,&cv_p);
    478:     VecGetArray(user->ci,&ci_p);

    481:     /* VecSet(user->cv,0.122);
    482:      VecSet(user->cv,6.9e-8);
    483:      VecSet(user->ci,6.9e-8); */

    486:     for (i=0; i<n/5; i++) {
    487:       if (i%5 <4) {
    488:         cv_p[i] = 0.0;
    489:         ci_p[i] = 1.0e-2;
    490:       } else {
    491:         cv_p[i] = 1.0e-2;
    492:         ci_p[i] = 0.0;
    493:       }
    494:     }
    495:     VecRestoreArray(user->cv,&cv_p);
    496:     VecRestoreArray(user->ci,&ci_p);

    499:     VecSet(user->eta,0.0);

    501:   }

    503:   DPsi(user);
    504:   VecCopy(user->DPsiv,user->wv);
    505:   VecCopy(user->DPsii,user->wi);

    507:   VecGetArray(X,&xx);
    508:   VecGetArray(user->cv,&cv_p);
    509:   VecGetArray(user->ci,&ci_p);
    510:   VecGetArray(user->eta,&eta_p);
    511:   VecGetArray(user->wv,&wv_p);
    512:   VecGetArray(user->wi,&wi_p);
    513:   for (i=0; i<n/5; i++)
    514:   {
    515:     xx[5*i]   = wv_p[i];
    516:     xx[5*i+1] = cv_p[i];
    517:     xx[5*i+2] = wi_p[i];
    518:     xx[5*i+3] = ci_p[i];
    519:     xx[5*i+4] = eta_p[i];

    523:   }

    525:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);
    526:   VecView(user->wv,view_out);
    527:   VecView(user->cv,view_out);
    528:   VecView(user->wi,view_out);
    529:   VecView(user->ci,view_out);
    530:   VecView(user->eta,view_out);
    531:   PetscViewerDestroy(&view_out);

    533:   VecRestoreArray(X,&xx);
    534:   VecRestoreArray(user->cv,&cv_p);
    535:   VecRestoreArray(user->ci,&ci_p);
    536:   VecRestoreArray(user->wv,&wv_p);
    537:   VecRestoreArray(user->wi,&wi_p);   VecRestoreArray(user->eta,&eta_p);
    538:   return(0);
    539: }

    543: PetscErrorCode SetRandomVectors(AppCtx *user)
    544: {
    546:   PetscInt       i,n,count=0;
    547:   PetscScalar    *w1,*w2,*Pv_p,*eta_p;

    549:   /* static PetscViewer viewer=0; */
    550:   static PetscRandom rand = 0;
    551:   static PetscInt    step = 0;

    554:   if (!rand) {
    555:     PetscRandomCreate(PETSC_COMM_WORLD,&rand);
    556:     PetscRandomSetFromOptions(rand);
    557:   }

    559:   VecSetRandom(user->work1,rand);
    560:   VecSetRandom(user->work2,rand);
    561:   VecGetArray(user->work1,&w1);
    562:   VecGetArray(user->work2,&w2);
    563:   VecGetArray(user->Pv,&Pv_p);
    564:   VecGetArray(user->eta,&eta_p);
    565:   VecGetLocalSize(user->work1,&n);
    566:   for (i=0; i<n; i++) {
    567:     if (eta_p[i]>=0.8 || w1[i]>user->P_casc) Pv_p[i]=0;
    568:     else {
    569:       Pv_p[i]=w2[i]*user->VG;
    570:       count  = count + 1;
    571:     }
    572:   }
    573:   step++;

    575:   VecCopy(user->Pv,user->Pi);
    576:   VecScale(user->Pi,0.9);
    577:   VecPointwiseMult(user->Piv,user->Pi,user->Pv);
    578:   VecRestoreArray(user->work1,&w1);
    579:   VecRestoreArray(user->work2,&w2);
    580:   VecRestoreArray(user->Pv,&Pv_p);
    581:   VecRestoreArray(user->eta,&eta_p);
    582:   printf("Number of radiations count %d out of total mesh points n %d\n",count,n);
    583:   return(0);

    585: }
    588: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
    589: {
    591:   AppCtx         *user=(AppCtx*)ctx;

    594:   MatMultAdd(user->M,X,user->q,F);
    595:   return(0);
    596: }

    600: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
    601: {
    603:   AppCtx         *user=(AppCtx*)ctx;

    606:   *flg = SAME_NONZERO_PATTERN;
    607:   MatCopy(user->M,*J,*flg);
    608:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
    609:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
    610:   return(0);
    611: }
    614: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
    615: {
    617:   PetscScalar    **l,**u;
    618:   PetscInt       xs,xm;
    619:   PetscInt       i;

    622:   DMDAVecGetArrayDOF(da,xl,&l);
    623:   DMDAVecGetArrayDOF(da,xu,&u);

    625:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

    628:   for (i=xs; i < xs+xm; i++) {
    629:     l[i][0] = -SNES_VI_INF;
    630:     l[i][1] = 0.0;
    631:     l[i][2] = -SNES_VI_INF;
    632:     l[i][3] = 0.0;
    633:     l[i][4] = 0.0;
    634:     u[i][0] = SNES_VI_INF;
    635:     u[i][1] = 1.0;
    636:     u[i][2] = SNES_VI_INF;
    637:     u[i][3] = 1.0;
    638:     u[i][4] = 1.0;
    639:   }

    642:   DMDAVecRestoreArrayDOF(da,xl,&l);
    643:   DMDAVecRestoreArrayDOF(da,xu,&u);
    644:   return(0);
    645: }

    649: PetscErrorCode GetParams(AppCtx *user)
    650: {
    652:   PetscBool      flg;

    655:   /* Set default parameters */
    656:   user->xmin     = 0.0; user->xmax = 128.0;
    657:   user->Dv       = 1.0; user->Di   = 1.0;
    658:   user->Evf      = 0.8; user->Eif  = 0.8;
    659:   user->A        = 1.0;
    660:   user->kBT      = 0.11;
    661:   user->kav      = 1.0;    user->kai    = 1.0; user->kaeta = 1.0;
    662:   user->Rsurf    = 10.0;   user->Rbulk  = 0.0;
    663:   user->L        = 10.0;   user->P_casc = 0.05;
    664:   user->T        = 1.0e-2; user->dt     = 1.0e-4;
    665:   user->VG       = 100.0;
    666:   user->initv    = .0001;
    667:   user->initeta  = 0.0;
    668:   user->graphics = PETSC_TRUE;

    670:   user->degenerate = PETSC_FALSE;
    671:   /* void growth */
    672:   user->voidgrowth = PETSC_FALSE;

    674:   PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);
    675:   PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);
    676:   PetscOptionsGetReal(NULL,"-T",&user->T,&flg);
    677:   PetscOptionsGetReal(NULL,"-dt",&user->dt,&flg);
    678:   PetscOptionsBool("-degenerate","Run with degenerate mobility\n","None",user->degenerate,&user->degenerate,&flg);
    679:   PetscOptionsReal("-smallnumber","Small number added to degenerate mobility\n","None",user->smallnumber,&user->smallnumber,&flg);
    680:   PetscOptionsBool("-voidgrowth","Use initial conditions for void growth\n","None",user->voidgrowth,&user->voidgrowth,&flg);
    681:   PetscOptionsBool("-graphics","Contour plot solutions at each timestep\n","None",user->graphics,&user->graphics,&flg);
    682:   return(0);
    683: }

    688: PetscErrorCode SetUpMatrices(AppCtx *user)
    689: {
    691:   PetscInt       nele,nen,i,n;
    692:   const PetscInt *ele;
    693:   PetscScalar    dt=user->dt,h;

    695:   PetscInt    idx[2];
    696:   PetscScalar eM_0[2][2],eM_2[2][2];
    697:   PetscScalar cv_sum, ci_sum;
    698:   Mat         M  =user->M;
    699:   Mat         M_0=user->M_0;
    700:   PetscInt    Mda;
    701:   PetscScalar *cv_p,*ci_p;
    702:   /* newly added */
    703:   Vec cvlocal,cilocal;

    706:   /* new stuff */
    707:   DMGetLocalVector(user->da2,&cvlocal);
    708:   DMGetLocalVector(user->da2,&cilocal);
    709:   DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
    710:   DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
    711:   DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
    712:   DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);

    714:   /* old stuff */
    715:   /*
    716:   VecGetArray(user->cv,&cv_p);
    717:   VecGetArray(user->ci,&ci_p);
    718:    */
    719:   /* new stuff */
    720:   VecGetArray(cvlocal,&cv_p);
    721:   VecGetArray(cilocal,&ci_p);

    723:   MatGetLocalSize(M,&n,NULL);
    724:   DMDAGetInfo(user->da1,NULL,&Mda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

    726:   h = (user->xmax-user->xmin)/Mda;

    728:   eM_0[0][0]=eM_0[1][1]=h/3.0;
    729:   eM_0[0][1]=eM_0[1][0]=h/6.0;

    731:   eM_2[0][0]=eM_2[1][1]=1.0/h;
    732:   eM_2[0][1]=eM_2[1][0]=-1.0/h;

    734:   /* Get local element info */
    735:   DMDAGetElements(user->da1,&nele,&nen,&ele);
    736:   /* for (i=0;i < nele + 1;i++) { */
    737:   for (i=0; i < nele; i++) {

    739:     idx[0] = ele[2*i]; idx[1] = ele[2*i+1];

    741:     if (user->degenerate) {
    742:       cv_sum = (2.0*user->smallnumber + cv_p[idx[0]]+cv_p[idx[1]])*user->Dv/(2.0*user->kBT);
    743:       ci_sum = (2.0*user->smallnumber + ci_p[idx[0]]+ci_p[idx[1]])*user->Di/(2.0*user->kBT);
    744:     } else {
    745:       cv_sum = user->initv*user->Dv/user->kBT;
    746:       ci_sum = user->initv*user->Di/user->kBT;
    747:     }

    750:     PetscInt row,cols[4],r,row_M_0,cols2[2];
    751:     /* PetscScalar vals[4],vals_M_0[1],vals2[2]; */
    752:     PetscScalar vals[4],vals_M_0[2],vals2[2];

    754:     for (r=0; r<2; r++) {
    755:       row_M_0 = idx[r];

    757:       vals_M_0[0]=eM_0[r][0];
    758:       vals_M_0[1]=eM_0[r][1];

    760:       /* vals_M_0[0] = h; */
    761:       MatSetValuesLocal(M_0,1,&row_M_0,2,idx,vals_M_0,ADD_VALUES);
    762:       /* MatSetValuesLocal(M_0,1,&row_M_0,1,&row_M_0,vals_M_0,INSERT_VALUES); */

    764:       row     = 5*idx[r];
    765:       cols[0] = 5*idx[0];     vals[0] = dt*eM_2[r][0]*cv_sum;
    766:       cols[1] = 5*idx[1];     vals[1] = dt*eM_2[r][1]*cv_sum;
    767:       cols[2] = 5*idx[0]+1;   vals[2] = eM_0[r][0];
    768:       cols[3] = 5*idx[1]+1;   vals[3] = eM_0[r][1];

    770:       /* Insert values in matrix M for 1st dof */
    771:       MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);

    773:       row     = 5*idx[r]+1;
    774:       cols[0] = 5*idx[0];     vals[0] = -eM_0[r][0];
    775:       cols[1] = 5*idx[1];     vals[1] = -eM_0[r][1];
    776:       cols[2] = 5*idx[0]+1;   vals[2] = user->kav*eM_2[r][0];
    777:       cols[3] = 5*idx[1]+1;   vals[3] = user->kav*eM_2[r][1];

    779:       /* Insert values in matrix M for 2nd dof */
    780:       MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);

    782:       row     = 5*idx[r]+2;
    783:       cols[0] = 5*idx[0]+2;   vals[0] = dt*eM_2[r][0]*ci_sum;
    784:       cols[1] = 5*idx[1]+2;   vals[1] = dt*eM_2[r][1]*ci_sum;
    785:       cols[2] = 5*idx[0]+3;   vals[2] = eM_0[r][0];
    786:       cols[3] = 5*idx[1]+3;   vals[3] = eM_0[r][1];
    787:       /* Insert values in matrix M for 3nd dof */
    788:       MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);

    790:       row     = 5*idx[r]+3;
    791:       cols[0] = 5*idx[0]+2;   vals[0] = -eM_0[r][0];
    792:       cols[1] = 5*idx[1]+2;   vals[1] = -eM_0[r][1];
    793:       cols[2] = 5*idx[0]+3;   vals[2] = user->kai*eM_2[r][0];
    794:       cols[3] = 5*idx[1]+3;   vals[3] = user->kai*eM_2[r][1];
    795:       /* Insert values in matrix M for 4th dof */
    796:       MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);

    798:       row      = 5*idx[r]+4;
    799:       cols2[0] = 5*idx[0]+4;   vals2[0] = eM_0[r][0]/dt + user->L*user->kaeta*eM_2[r][0];
    800:       cols2[1] = 5*idx[1]+4;   vals2[1] = eM_0[r][1]/dt + user->L*user->kaeta*eM_2[r][1];
    801:       MatSetValuesLocal(M,1,&row,2,cols2,vals2,ADD_VALUES);
    802:     }
    803:   }

    807:   /* new */
    808:   VecRestoreArray(cvlocal,&cv_p);
    809:   VecRestoreArray(cilocal,&ci_p);
    810:   DMRestoreLocalVector(user->da2,&cvlocal);
    811:   DMRestoreLocalVector(user->da2,&cilocal);
    812:   /*
    813:   VecRestoreArray(user->cv,&cv_p);
    814:    VecRestoreArray(user->ci,&ci_p);*/
    815:   MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
    816:   MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);

    818:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
    819:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

    821:   DMDARestoreElements(user->da1,&nele,&nen,&ele);
    822:   return(0);
    823: }

    828: PetscErrorCode UpdateMatrices(AppCtx *user)
    829: {
    831:   PetscInt       nele,nen,i,n,Mda;
    832:   const PetscInt *ele;

    834:   PetscInt    idx[2];
    835:   PetscScalar eM_2[2][2],h;
    836:   Mat         M=user->M;
    837:   PetscScalar *cv_p,*ci_p,cv_sum,ci_sum;
    838:   /* newly added */
    839:   Vec cvlocal,cilocal;

    842:   /*new stuff */
    843:   DMGetLocalVector(user->da2,&cvlocal);
    844:   DMGetLocalVector(user->da2,&cilocal);
    845:   DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
    846:   DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
    847:   DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
    848:   DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
    849:   /* old stuff */
    850:   /*
    851:   VecGetArray(user->cv,&cv_p);
    852:   VecGetArray(user->ci,&ci_p);
    853:    */

    855:   /* new stuff */
    856:   VecGetArray(cvlocal,&cv_p);
    857:   VecGetArray(cilocal,&ci_p);

    859:   /* Create the mass matrix M_0 */
    860:   MatGetLocalSize(M,&n,NULL);
    861:   DMDAGetInfo(user->da1,NULL,&Mda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

    863:   h = (user->xmax-user->xmin)/Mda;

    865:   /* Get local element info */
    866:   DMDAGetElements(user->da1,&nele,&nen,&ele);

    868:   for (i=0; i<nele; i++) {
    869:     idx[0] = ele[2*i]; idx[1] = ele[2*i+1];

    871:     PetscInt    r,row,cols[2];
    872:     PetscScalar vals[2];

    874:     for (r=0; r<2; r++) {
    875:       row     = 5*idx[r];
    876:       cols[0] = 5*idx[0];     vals[0] = 0.0;
    877:       cols[1] = 5*idx[1];     vals[1] = 0.0;

    880:       /* Insert values in matrix M for 1st dof */
    881:       MatSetValuesLocal(M,1,&row,2,cols,vals,INSERT_VALUES);

    883:       row     = 5*idx[r]+2;
    884:       cols[0] = 5*idx[0]+2;   vals[0] = 0.0;
    885:       cols[1] = 5*idx[1]+2;   vals[1] = 0.0;

    887:       /* Insert values in matrix M for 3nd dof */
    888:       MatSetValuesLocal(M,1,&row,2,cols,vals,INSERT_VALUES);
    889:     }
    890:   }

    892:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
    893:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

    895:   eM_2[0][0]=eM_2[1][1]=1.0/h;
    896:   eM_2[0][1]=eM_2[1][0]=-1.0/h;

    898:   for (i=0; i<nele; i++) {
    899:     idx[0] = ele[2*i]; idx[1] = ele[2*i+1];

    901:     PetscInt    row,cols[2],r;
    902:     PetscScalar vals[2];

    904:     for (r=0; r<2; r++) {

    906:       if (user->degenerate) {
    907:         cv_sum = (2.0*user->smallnumber + cv_p[idx[0]] + cv_p[idx[1]])*user->Dv/(2.0*user->kBT);
    908:         ci_sum = (2.0*user->smallnumber + ci_p[idx[0]] + ci_p[idx[1]])*user->Di/(2.0*user->kBT);
    909:       } else {
    910:         cv_sum = user->initv*user->Dv/user->kBT;
    911:         ci_sum = user->initv*user->Di/user->kBT;
    912:       }

    914:       row     = 5*idx[r];
    915:       cols[0] = 5*idx[0];     vals[0] = user->dt*eM_2[r][0]*cv_sum;
    916:       cols[1] = 5*idx[1];     vals[1] = user->dt*eM_2[r][1]*cv_sum;

    918:       /* Insert values in matrix M for 1st dof */
    919:       MatSetValuesLocal(M,1,&row,2,cols,vals,ADD_VALUES);

    921:       row     = 5*idx[r]+2;
    922:       cols[0] = 5*idx[0]+2;   vals[0] = user->dt*eM_2[r][0]*ci_sum;
    923:       cols[1] = 5*idx[1]+2;   vals[1] = user->dt*eM_2[r][1]*ci_sum;

    925:       /* Insert values in matrix M for 3nd dof */
    926:       MatSetValuesLocal(M,1,&row,2,cols,vals,ADD_VALUES);
    927:     }
    928:   }

    930:   /* old stuff
    931:   VecRestoreArray(user->cv,&cv_p);
    932:   VecRestoreArray(user->ci,&ci_p);
    933:    */

    935:   /* new stuff */
    936:   VecRestoreArray(cvlocal,&cv_p);
    937:   VecRestoreArray(cilocal,&ci_p);
    938:   DMRestoreLocalVector(user->da2,&cvlocal);
    939:   DMRestoreLocalVector(user->da2,&cilocal);

    941:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
    942:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
    943:   DMDARestoreElements(user->da1,&nele,&nen,&ele);
    944:   return(0);
    945: }

    950: PetscErrorCode CheckRedundancy(SNES snes, IS act, IS *outact, DM da)
    951: {
    953:   PetscScalar    **uin,**uout;
    954:   Vec            UIN, UOUT;
    955:   PetscInt       xs,xm,*outindex;
    956:   const PetscInt *index;
    957:   PetscInt       k,i,l,n,M,cnt=0;

    960:   DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
    961:   DMGetGlobalVector(da,&UIN);
    962:   VecSet(UIN,0.0);
    963:   DMGetLocalVector(da,&UOUT);
    964:   DMDAVecGetArrayDOF(da,UIN,&uin);
    965:   ISGetIndices(act,&index);
    966:   ISGetLocalSize(act,&n);
    967:   for (k=0; k<n; k++) {
    968:     l        = index[k]%5;
    969:     i        = index[k]/5;
    970:     uin[i][l]=1.0;
    971:   }
    972:   printf("Number of active constraints before applying redundancy %d\n",n);
    973:   ISRestoreIndices(act,&index);
    974:   DMDAVecRestoreArrayDOF(da,UIN,&uin);
    975:   DMGlobalToLocalBegin(da,UIN,INSERT_VALUES,UOUT);
    976:   DMGlobalToLocalEnd(da,UIN,INSERT_VALUES,UOUT);
    977:   DMDAVecGetArrayDOF(da,UOUT,&uout);

    979:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

    981:   for (i=xs; i < xs+xm; i++) {
    982:     if (uout[i-1][1] && uout[i][1] && uout[i+1][1]) uout[i][0] = 1.0;
    983:     if (uout[i-1][3] && uout[i][3] && uout[i+1][3]) uout[i][2] = 1.0;
    984:   }

    986:   for (i=xs; i < xs+xm; i++) {
    987:     for (l=0; l < 5; l++) {
    988:       if (uout[i][l]) cnt++;
    989:     }
    990:   }

    992:   printf("Number of active constraints after applying redundancy %d\n",cnt);

    995:   PetscMalloc(cnt*sizeof(PetscInt),&outindex);
    996:   cnt  = 0;

    998:   for (i=xs; i < xs+xm; i++) {
    999:     for (l=0; l < 5; l++) {
    1000:       if (uout[i][l]) outindex[cnt++] = 5*(i)+l;
    1001:     }
    1002:   }

    1005:   ISCreateGeneral(PETSC_COMM_WORLD,cnt,outindex,PETSC_OWN_POINTER,outact);
    1006:   DMDAVecRestoreArrayDOF(da,UOUT,&uout);
    1007:   DMRestoreGlobalVector(da,&UIN);
    1008:   DMRestoreLocalVector(da,&UOUT);
    1009:   return(0);
    1010: }