Actual source code: ex63.c

  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_ls_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_ls 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;  /* olution 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);
 66: 
 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,PETSC_NULL,&user.da1);
 71:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC, -4, 5, 1,PETSC_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,PETSC_NULL,&user.da2);

 75:   /* Set Element type (triangular) */
 76:   DMDASetElementType(user.da1,DMDA_ELEMENT_P1);
 77:   DMDASetElementType(user.da2,DMDA_ELEMENT_P1);
 78: 
 79:   /* Set x and y coordinates */
 80:   DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
 81:   DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_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);
 88: 
 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:   DMGetMatrix(user.da1,MATAIJ,&user.M);
113:   /* Get the (usual) mass matrix structure from da2 */
114:   DMGetMatrix(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);
119: 
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:   SNESSetType(snes,SNESVI);
127:   SNESSetFromOptions(snes);
128:   SetVariableBounds(user.da1,xl,xu);
129:   //SNESVISetRedundancyCheck(snes,(PetscErrorCode (*)(SNES,IS,IS*,void*))CheckRedundancy,user.da1_clone);
130:   SNESVISetVariableBounds(snes,xl,xu);
131:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
132:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_p",FILE_MODE_WRITE,&view_p);
133:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
134:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);
135:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);
136:   /* PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds); */

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

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

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

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

159:     VecView(user.Pv,view_p);
160:     Update_q(&user);
161:     VecView(user.q,view_q);
162:     MatView(user.M,view_mat);
163:     SNESSolve(snes,PETSC_NULL,x);
164: 
165:     VecView(x,view_out);
166:     if (user.graphics) {
167:       VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
168:     }
169: 
170:     PetscInt its;
171:     SNESGetIterationNumber(snes,&its);
172:     PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4g in %d iterations\n",t,its);

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

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

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

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

246: 
247:  VecMax(user->cv,&pv1,&maxv);
248:   VecMax(user->ci,&pi1,&maxi);
249:   VecMax(user->eta,&peta1,&maxeta);
250:   VecMin(user->cv,&pv2,&minv);
251:   VecMin(user->ci,&pi2,&mini);
252:   VecMin(user->eta,&peta2,&mineta);
253: 
254: 
255:   VecGetArray(X,&xx);
256:   VecGetArray(user->wv,&wv_p);
257:   VecGetArray(user->cv,&cv_p);
258:   VecGetArray(user->wi,&wi_p);
259:   VecGetArray(user->ci,&ci_p);
260:   VecGetArray(user->eta,&eta_p);
261: 
262: 
263:   for(i=0;i<n;i++) {
264:     wv_p[i] = xx[5*i];
265:     cv_p[i] = xx[5*i+1];
266:     wi_p[i] = xx[5*i+2];
267:     ci_p[i] = xx[5*i+3];
268:     eta_p[i] = xx[5*i+4];
269:   }
270:   VecRestoreArray(X,&xx);
271:   VecRestoreArray(user->wv,&wv_p);
272:   VecRestoreArray(user->cv,&cv_p);
273:   VecRestoreArray(user->wi,&wi_p);
274:   VecRestoreArray(user->ci,&ci_p);
275:   VecRestoreArray(user->eta,&eta_p);
276: 
277:   return(0);
278: }

282: PetscErrorCode Update_q(AppCtx *user)
283: {
285:   PetscScalar    *q_p, *w1, *w2;
286:   PetscInt       i,n;
287:   PetscScalar    norm1;

290: 
291:   VecPointwiseMult(user->Riv,user->eta,user->eta);
292:   VecScale(user->Riv,user->Rsurf);
293:   VecShift(user->Riv,user->Rbulk);
294:   VecPointwiseMult(user->Riv,user->ci,user->Riv);
295:   VecPointwiseMult(user->Riv,user->cv,user->Riv);
296: 
297:   VecCopy(user->Riv,user->work1);
298:   VecScale(user->work1,user->dt);
299:   VecAXPY(user->work1,-1.0,user->cv);
300:   MatMult(user->M_0,user->work1,user->work2);
301: 
302:   VecGetArray(user->q,&q_p);
303:   VecGetArray(user->work1,&w1);
304:   VecGetArray(user->work2,&w2);
305:   VecGetLocalSize(user->wv,&n);
306:   for (i=0;i<n;i++) {
307:     q_p[5*i]=w2[i];
308:   }

310:   MatMult(user->M_0,user->DPsiv,user->work1);
311:   for (i=0;i<n;i++) {
312:     q_p[5*i+1]=w1[i];
313:   }

315:   VecCopy(user->Riv,user->work1);
316:   VecScale(user->work1,user->dt);
317:   VecAXPY(user->work1,-1.0,user->ci);
318:   MatMult(user->M_0,user->work1,user->work2);
319:   for (i=0;i<n;i++) {
320:     q_p[5*i+2]=w2[i];
321:   }

323:   MatMult(user->M_0,user->DPsii,user->work1);
324:   for (i=0;i<n;i++) {
325:     q_p[5*i+3]=w1[i];
326:   }

328:   VecCopy(user->DPsieta,user->work1);
329:   VecScale(user->work1,user->L);
330:   VecAXPY(user->work1,-1.0/user->dt,user->eta);
331:   MatMult(user->M_0,user->work1,user->work2);
332:   for (i=0;i<n;i++) {
333:     q_p[5*i+4]=w2[i];
334:   }
335: 
336:   VecRestoreArray(user->q,&q_p);
337:   VecRestoreArray(user->work1,&w1);
338:   VecRestoreArray(user->work2,&w2);


341: 
342:   return(0);
343: }

347: PetscErrorCode DPsi(AppCtx* user)
348: {
349:   PetscErrorCode  ierr;
350:   PetscScalar     Evf=user->Evf,Eif=user->Eif,kBT=user->kBT,A=user->A;
351:   PetscScalar     *cv_p,*ci_p,*eta_p,*logcv_p,*logci_p,*logcvi_p,*DPsiv_p,*DPsii_p,*DPsieta_p;
352:   PetscInt        n,i;


356:   VecGetLocalSize(user->cv,&n);
357:   VecGetArray(user->cv,&cv_p);
358:   VecGetArray(user->ci,&ci_p);
359:   VecGetArray(user->eta,&eta_p);
360:   VecGetArray(user->logcv,&logcv_p);
361:   VecGetArray(user->logci,&logci_p);
362:   VecGetArray(user->logcvi,&logcvi_p);
363:   VecGetArray(user->DPsiv,&DPsiv_p);
364:   VecGetArray(user->DPsii,&DPsii_p);
365:   VecGetArray(user->DPsieta,&DPsieta_p);

367:   Llog(user->cv,user->logcv);
368:   Llog(user->ci,user->logci);

370:   VecCopy(user->cv,user->cvi);
371:   VecAXPY(user->cvi,1.0,user->ci);
372:   VecScale(user->cvi,-1.0);
373:   VecShift(user->cvi,1.0);
374:   Llog(user->cvi,user->logcvi);

376:   for (i=0;i<n;i++)
377:   {
378:     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);

380:     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] ;
381: 
382:     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]);

384: 
385:   }

387:   VecRestoreArray(user->cv,&cv_p);
388:   VecRestoreArray(user->ci,&ci_p);
389:   VecRestoreArray(user->eta,&eta_p);
390:   VecGetArray(user->logcv,&logcv_p);
391:   VecGetArray(user->logci,&logci_p);
392:   VecGetArray(user->logcvi,&logcvi_p);
393:   VecRestoreArray(user->DPsiv,&DPsiv_p);
394:   VecRestoreArray(user->DPsii,&DPsii_p);
395:   VecRestoreArray(user->DPsieta,&DPsieta_p);


398:   return(0);

400: }


405: PetscErrorCode Llog(Vec X, Vec Y)
406: {
407:   PetscErrorCode    ierr;
408:   PetscScalar       *x,*y;
409:   PetscInt          n,i;

412: 
413:   VecGetArray(X,&x);
414:   VecGetArray(Y,&y);
415:   VecGetLocalSize(X,&n);
416:   for (i=0;i<n;i++) {
417:     if (x[i] < 1.0e-12) {
418:       y[i] = log(1.0e-12);
419:     }
420:     else {
421:       y[i] = log(x[i]);
422:     }
423:   }
424: 
425:   return(0);
426: }

430: PetscErrorCode SetInitialGuess(Vec X,AppCtx* user)
431: {
432:   PetscErrorCode    ierr;


435:   PetscInt         n,i,Mda;
436:   PetscScalar           *xx,*cv_p,*ci_p,*wv_p,*wi_p,*eta_p;
437:   PetscViewer      view_out;
438:   /* needed for the void growth case */
439:   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;
440:   PetscInt          nele,nen,idx[2];
441:   const PetscInt    *ele;
442:   PetscScalar       x[2];
443:   Vec               coords;
444:   const PetscScalar *_coords;
445:   PetscViewer       view;
446:   PetscScalar       xwidth = user->xmax - user->xmin;


450:   VecGetLocalSize(X,&n);
451: 
452:   if (user->voidgrowth) {
453:     DMDAGetInfo(user->da2,PETSC_NULL,&Mda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
454:     DMDAGetGhostedCoordinates(user->da2,&coords);
455:     VecGetArrayRead(coords,&_coords);

457:     h = (user->xmax-user->xmin)/Mda;
458:     xmid = (user->xmax + user->xmin)/2.0;
459:     lambda = 4.0*h;

461:     DMDAGetElements(user->da2,&nele,&nen,&ele);
462:     for (i=0;i < nele; i++) {
463:       idx[0] = ele[2*i]; idx[1] = ele[2*i+1];
464: 
465:       x[0] = _coords[idx[0]];
466:       x[1] = _coords[idx[1]];

468: 
469:       PetscInt k;
470:       PetscScalar vals_cv[2],vals_ci[2],vals_eta[2],s,hhr,r;
471:       for (k=0; k < 2 ; k++) {
472:         s = PetscAbs(x[k] - xmid);
473:         if (s < xwidth*(5.0/64.0)) {
474:           vals_cv[k] = cv_v;
475:           vals_ci[k] = ci_v;
476:           vals_eta[k] = eta_v;
477:         } else if (s>= xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0) ) {
478:           //r = (s - xwidth*(6.0/64.0) )/(0.5*lambda);
479:           r = (s - xwidth*(6.0/64.0) )/(xwidth/64.0);
480:           hhr = 0.25*(-r*r*r + 3*r + 2);
481:           vals_cv[k] = cv_m + (1.0 - hhr)*(cv_v - cv_m);
482:           vals_ci[k] = ci_m + (1.0 - hhr)*(ci_v - ci_m);
483:           vals_eta[k] = eta_m + (1.0 - hhr)*(eta_v - eta_m);
484:         } else {
485:           vals_cv[k] = cv_m;
486:           vals_ci[k] = ci_m;
487:           vals_eta[k] = eta_m;
488:         }
489:       }

491:       VecSetValuesLocal(user->cv,2,idx,vals_cv,INSERT_VALUES);
492:       VecSetValuesLocal(user->ci,2,idx,vals_ci,INSERT_VALUES);
493:       VecSetValuesLocal(user->eta,2,idx,vals_eta,INSERT_VALUES);
494:     }
495:     VecAssemblyBegin(user->cv);
496:     VecAssemblyEnd(user->cv);
497:     VecAssemblyBegin(user->ci);
498:     VecAssemblyEnd(user->ci);
499:     VecAssemblyBegin(user->eta);
500:     VecAssemblyEnd(user->eta);

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

505: 
506:   } else {

508:     VecGetArray(user->cv,&cv_p);
509:     VecGetArray(user->ci,&ci_p);

511: 
512:     /* VecSet(user->cv,0.122);
513:      VecSet(user->cv,6.9e-8);
514:      VecSet(user->ci,6.9e-8); */
515: 
516: 
517:     for (i=0;i<n/5;i++)
518:     {
519:       if (i%5 <4 ) {
520:         cv_p[i] = 0.0;
521:         ci_p[i] = 1.0e-2;
522:       }
523:       else {
524:         cv_p[i] = 1.0e-2;
525:         ci_p[i] = 0.0;
526:       }
527:     }
528:     VecRestoreArray(user->cv,&cv_p);
529:     VecRestoreArray(user->ci,&ci_p);


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

534:   }

536:   DPsi(user);
537:   VecCopy(user->DPsiv,user->wv);
538:   VecCopy(user->DPsii,user->wi);

540:   VecGetArray(X,&xx);
541:   VecGetArray(user->cv,&cv_p);
542:   VecGetArray(user->ci,&ci_p);
543:   VecGetArray(user->eta,&eta_p);
544:   VecGetArray(user->wv,&wv_p);
545:   VecGetArray(user->wi,&wi_p);
546:   for (i=0;i<n/5;i++)
547:   {
548:     xx[5*i]=wv_p[i];
549:     xx[5*i+1]=cv_p[i];
550:     xx[5*i+2]=wi_p[i];
551:     xx[5*i+3]=ci_p[i];
552:     xx[5*i+4]=eta_p[i];

554: 
555: 
556:   }

558:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);
559:   VecView(user->wv,view_out);
560:   VecView(user->cv,view_out);
561:   VecView(user->wi,view_out);
562:   VecView(user->ci,view_out);
563:   VecView(user->eta,view_out);
564:   PetscViewerDestroy(&view_out);

566:   VecRestoreArray(X,&xx);
567:   VecRestoreArray(user->cv,&cv_p);
568:   VecRestoreArray(user->ci,&ci_p);
569:   VecRestoreArray(user->wv,&wv_p);
570:   VecRestoreArray(user->wi,&wi_p);   VecRestoreArray(user->eta,&eta_p);
571:   return(0);
572: }

576: PetscErrorCode SetRandomVectors(AppCtx* user)
577: {
579:   PetscInt       i,n,count=0;
580:   PetscScalar    *w1,*w2,*Pv_p,*eta_p;

582:   /* static PetscViewer viewer=0; */
583:   static PetscRandom rand = 0;
584:   static PetscInt    step = 0;

587:   if (!rand) {
588:     PetscRandomCreate(PETSC_COMM_WORLD,&rand);
589:     PetscRandomSetFromOptions(rand);
590:   }
591: 
592:   VecSetRandom(user->work1,rand);
593:   VecSetRandom(user->work2,rand);
594:   VecGetArray(user->work1,&w1);
595:   VecGetArray(user->work2,&w2);
596:   VecGetArray(user->Pv,&Pv_p);
597:   VecGetArray(user->eta,&eta_p);
598:   VecGetLocalSize(user->work1,&n);
599:   for (i=0;i<n;i++) {
600:     if (eta_p[i]>=0.8 || w1[i]>user->P_casc){
601:       Pv_p[i]=0;
602:     }
603:     else
604:     {
605:       Pv_p[i]=w2[i]*user->VG;
606:       count = count + 1;
607:     }
608:   }
609:   step++;

611:   VecCopy(user->Pv,user->Pi);
612:   VecScale(user->Pi,0.9);
613:   VecPointwiseMult(user->Piv,user->Pi,user->Pv);
614:   VecRestoreArray(user->work1,&w1);
615:   VecRestoreArray(user->work2,&w2);
616:   VecRestoreArray(user->Pv,&Pv_p);
617:   VecRestoreArray(user->eta,&eta_p);
618:   printf("Number of radiations count %d out of total mesh points n %d\n",count,n);

620:   return(0);
621: 
622: }
625: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ctx)
626: {
628:   AppCtx         *user=(AppCtx*)ctx;
629: 
631:   MatMultAdd(user->M,X,user->q,F);
632:   return(0);
633: }

637: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
638: {
640:   AppCtx         *user=(AppCtx*)ctx;
641: 
643:   *flg = SAME_NONZERO_PATTERN;
644:   MatCopy(user->M,*J,*flg);
645:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
646:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
647:   return(0);
648: }
651: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
652: {
654:   PetscScalar    **l,**u;
655:   PetscInt       xs,xm;
656:   PetscInt       i;
657: 
659:   DMDAVecGetArrayDOF(da,xl,&l);
660:   DMDAVecGetArrayDOF(da,xu,&u);
661: 
662:   DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
663: 
664: 
665:   for(i=xs; i < xs+xm;i++) {
666:     l[i][0] = -SNES_VI_INF;
667:     l[i][1] = 0.0;
668:     l[i][2] = -SNES_VI_INF;
669:     l[i][3] = 0.0;
670:     l[i][4] = 0.0;
671:     u[i][0] = SNES_VI_INF;
672:     u[i][1] = 1.0;
673:     u[i][2] = SNES_VI_INF;
674:     u[i][3] = 1.0;
675:     u[i][4] = 1.0;
676:   }
677: 
678: 
679:   DMDAVecRestoreArrayDOF(da,xl,&l);
680:   DMDAVecRestoreArrayDOF(da,xu,&u);
681:   return(0);
682: }

686: PetscErrorCode GetParams(AppCtx* user)
687: {
689:   PetscBool      flg;
690: 
692: 
693:   /* Set default parameters */
694:   user->xmin = 0.0; user->xmax = 128.0;
695:   user->Dv = 1.0; user->Di=1.0;
696:   user->Evf = 0.8; user->Eif = 0.8;
697:   user->A = 1.0;
698:   user->kBT = 0.11;
699:   user->kav = 1.0; user->kai = 1.0; user->kaeta = 1.0;
700:   user->Rsurf = 10.0; user->Rbulk = 0.0;
701:   user->L = 10.0; user->P_casc = 0.05;
702:   user->T = 1.0e-2;    user->dt = 1.0e-4;
703:   user->VG = 100.0;
704:   user->initv = .0001;
705:   user->initeta = 0.0;
706:   user->graphics = PETSC_TRUE;
707: 
708:   user->degenerate = PETSC_FALSE;
709:   /* void growth */
710:   user->voidgrowth = PETSC_FALSE;

712:   PetscOptionsGetReal(PETSC_NULL,"-xmin",&user->xmin,&flg);
713:   PetscOptionsGetReal(PETSC_NULL,"-xmax",&user->xmax,&flg);
714:   PetscOptionsGetReal(PETSC_NULL,"-T",&user->T,&flg);
715:   PetscOptionsGetReal(PETSC_NULL,"-dt",&user->dt,&flg);
716:   PetscOptionsBool("-degenerate","Run with degenerate mobility\n","None",user->degenerate,&user->degenerate,&flg);
717:   PetscOptionsReal("-smallnumber","Small number added to degenerate mobility\n","None",user->smallnumber,&user->smallnumber,&flg);
718:  PetscOptionsBool("-voidgrowth","Use initial conditions for void growth\n","None",user->voidgrowth,&user->voidgrowth,&flg);
719: PetscOptionsBool("-graphics","Contour plot solutions at each timestep\n","None",user->graphics,&user->graphics,&flg);
720: 
721:   return(0);
722:  }


727: PetscErrorCode SetUpMatrices(AppCtx* user)
728: {
729:   PetscErrorCode    ierr;
730:   PetscInt          nele,nen,i,n;
731:   const PetscInt    *ele;
732:   PetscScalar       dt=user->dt,h;
733: 
734:   PetscInt          idx[2];
735:   PetscScalar       eM_0[2][2],eM_2[2][2];
736:   PetscScalar       cv_sum, ci_sum;
737:   Mat               M=user->M;
738:   Mat               M_0=user->M_0;
739:   PetscInt          Mda;
740:   PetscScalar       *cv_p,*ci_p;
741:   /* newly added */
742:   Vec               cvlocal,cilocal;
743: 

746:   /* new stuff */
747:   DMGetLocalVector(user->da2,&cvlocal);
748:   DMGetLocalVector(user->da2,&cilocal);
749:   DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
750:   DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
751:   DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
752:   DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);

754:   /* old stuff */
755:   /*
756:   VecGetArray(user->cv,&cv_p);
757:   VecGetArray(user->ci,&ci_p);
758:    */
759:   /* new stuff */
760:   VecGetArray(cvlocal,&cv_p);
761:   VecGetArray(cilocal,&ci_p);

763:   MatGetLocalSize(M,&n,PETSC_NULL);
764:   DMDAGetInfo(user->da1,PETSC_NULL,&Mda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);

766:   h = (user->xmax-user->xmin)/Mda;
767: 
768:   eM_0[0][0]=eM_0[1][1]=h/3.0;
769:   eM_0[0][1]=eM_0[1][0]=h/6.0;

771:   eM_2[0][0]=eM_2[1][1]=1.0/h;
772:   eM_2[0][1]=eM_2[1][0]=-1.0/h;

774:   /* Get local element info */
775:   DMDAGetElements(user->da1,&nele,&nen,&ele);
776:   //for(i=0;i < nele + 1;i++) {
777:   for(i=0;i < nele;i++) {

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

781:       if (user->degenerate) {
782:         cv_sum = (2.0*user->smallnumber + cv_p[idx[0]]+cv_p[idx[1]])*user->Dv/(2.0*user->kBT);
783:         ci_sum = (2.0*user->smallnumber + ci_p[idx[0]]+ci_p[idx[1]])*user->Di/(2.0*user->kBT);
784:       } else {
785:         cv_sum = user->initv*user->Dv/user->kBT;
786:         ci_sum = user->initv*user->Di/user->kBT;
787:       }


790:     PetscInt    row,cols[4],r,row_M_0,cols2[2];
791:     //PetscScalar vals[4],vals_M_0[1],vals2[2];
792:     PetscScalar vals[4],vals_M_0[2],vals2[2];
793: 
794:     for(r=0;r<2;r++) {
795:       row_M_0 = idx[r];
796: 
797:       vals_M_0[0]=eM_0[r][0];
798:       vals_M_0[1]=eM_0[r][1];
799: 
800:       //vals_M_0[0] = h;
801:       MatSetValuesLocal(M_0,1,&row_M_0,2,idx,vals_M_0,ADD_VALUES);
802:       //MatSetValuesLocal(M_0,1,&row_M_0,1,&row_M_0,vals_M_0,INSERT_VALUES);
803: 
804:       row = 5*idx[r];
805:       cols[0] = 5*idx[0];     vals[0] = dt*eM_2[r][0]*cv_sum;
806:       cols[1] = 5*idx[1];     vals[1] = dt*eM_2[r][1]*cv_sum;
807:       cols[2] = 5*idx[0]+1;   vals[2] = eM_0[r][0];
808:       cols[3] = 5*idx[1]+1;   vals[3] = eM_0[r][1];
809: 
810:       /* Insert values in matrix M for 1st dof */
811:       MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);
812: 
813:       row = 5*idx[r]+1;
814:       cols[0] = 5*idx[0];     vals[0] = -eM_0[r][0];
815:       cols[1] = 5*idx[1];     vals[1] = -eM_0[r][1];
816:       cols[2] = 5*idx[0]+1;   vals[2] = user->kav*eM_2[r][0];
817:       cols[3] = 5*idx[1]+1;   vals[3] = user->kav*eM_2[r][1];

819:       /* Insert values in matrix M for 2nd dof */
820:       MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);
821: 
822:       row = 5*idx[r]+2;
823:       cols[0] = 5*idx[0]+2;   vals[0] = dt*eM_2[r][0]*ci_sum;
824:       cols[1] = 5*idx[1]+2;   vals[1] = dt*eM_2[r][1]*ci_sum;
825:       cols[2] = 5*idx[0]+3;   vals[2] = eM_0[r][0];
826:       cols[3] = 5*idx[1]+3;   vals[3] = eM_0[r][1];
827:       /* Insert values in matrix M for 3nd dof */
828:       MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);
829: 
830:       row = 5*idx[r]+3;
831:       cols[0] = 5*idx[0]+2;   vals[0] = -eM_0[r][0];
832:       cols[1] = 5*idx[1]+2;   vals[1] = -eM_0[r][1];
833:       cols[2] = 5*idx[0]+3;   vals[2] = user->kai*eM_2[r][0];
834:       cols[3] = 5*idx[1]+3;   vals[3] = user->kai*eM_2[r][1];
835:       /* Insert values in matrix M for 4th dof */
836:       MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);

838:       row = 5*idx[r]+4;
839:       cols2[0] = 5*idx[0]+4;   vals2[0] = eM_0[r][0]/dt + user->L*user->kaeta*eM_2[r][0];
840:       cols2[1] = 5*idx[1]+4;   vals2[1] = eM_0[r][1]/dt + user->L*user->kaeta*eM_2[r][1];
841:       MatSetValuesLocal(M,1,&row,2,cols2,vals2,ADD_VALUES);
842:     }
843:   }

845: 
846: 
847:   /* new */
848:   VecRestoreArray(cvlocal,&cv_p);
849:   VecRestoreArray(cilocal,&ci_p);
850:   DMRestoreLocalVector(user->da2,&cvlocal);
851:   DMRestoreLocalVector(user->da2,&cilocal);
852:   /*
853:   VecRestoreArray(user->cv,&cv_p);
854:    VecRestoreArray(user->ci,&ci_p);*/
855:   MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
856:   MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);
857: 
858:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
859:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
860: 
861:   DMDARestoreElements(user->da1,&nele,&nen,&ele);
862: 
863: 
864:   return(0);
865: }


870: PetscErrorCode UpdateMatrices(AppCtx* user)
871: {
872:   PetscErrorCode    ierr;
873:   PetscInt          nele,nen,i,n,Mda;
874:   const PetscInt    *ele;
875: 
876:   PetscInt          idx[2];
877:   PetscScalar       eM_2[2][2],h;
878:   Mat               M=user->M;
879:   PetscScalar       *cv_p,*ci_p,cv_sum,ci_sum;
880:   /* newly added */
881:   Vec               cvlocal,cilocal;

884: 
885:   /*new stuff */
886:   DMGetLocalVector(user->da2,&cvlocal);
887:   DMGetLocalVector(user->da2,&cilocal);
888:   DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
889:   DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
890:   DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
891:   DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
892:   /* old stuff */
893:   /*
894:   VecGetArray(user->cv,&cv_p);
895:   VecGetArray(user->ci,&ci_p);
896:    */

898:   /* new stuff */
899:   VecGetArray(cvlocal,&cv_p);
900:   VecGetArray(cilocal,&ci_p);

902:   /* Create the mass matrix M_0 */
903:   MatGetLocalSize(M,&n,PETSC_NULL);
904:   DMDAGetInfo(user->da1,PETSC_NULL,&Mda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);

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

908:   /* Get local element info */
909:   DMDAGetElements(user->da1,&nele,&nen,&ele);
910: 
911:   for(i=0;i<nele;i++) {
912:     idx[0] = ele[2*i]; idx[1] = ele[2*i+1];
913: 
914:     PetscInt     r,row,cols[2];
915:     PetscScalar  vals[2];

917:     for (r=0;r<2;r++) {
918:       row = 5*idx[r];
919:       cols[0] = 5*idx[0];     vals[0] = 0.0;
920:       cols[1] = 5*idx[1];     vals[1] = 0.0;
921: 

923:       /* Insert values in matrix M for 1st dof */
924:       MatSetValuesLocal(M,1,&row,2,cols,vals,INSERT_VALUES);
925: 
926:       row = 5*idx[r]+2;
927:       cols[0] = 5*idx[0]+2;   vals[0] = 0.0;
928:       cols[1] = 5*idx[1]+2;   vals[1] = 0.0;

930:       /* Insert values in matrix M for 3nd dof */
931:       MatSetValuesLocal(M,1,&row,2,cols,vals,INSERT_VALUES);
932:     }
933:   }
934: 
935:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
936:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

938:   eM_2[0][0]=eM_2[1][1]=1.0/h;
939:   eM_2[0][1]=eM_2[1][0]=-1.0/h;

941:   for(i=0;i<nele;i++) {
942:     idx[0] = ele[2*i]; idx[1] = ele[2*i+1];
943: 
944:     PetscInt    row,cols[2],r;
945:     PetscScalar vals[2];

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

949:       if (user->degenerate) {
950:         cv_sum = (2.0*user->smallnumber + cv_p[idx[0]] + cv_p[idx[1]])*user->Dv/(2.0*user->kBT);
951:         ci_sum = (2.0*user->smallnumber + ci_p[idx[0]] + ci_p[idx[1]])*user->Di/(2.0*user->kBT);
952:       } else {
953:         cv_sum = user->initv*user->Dv/user->kBT;
954:         ci_sum = user->initv*user->Di/user->kBT;
955:       }

957:       row = 5*idx[r];
958:       cols[0] = 5*idx[0];     vals[0] = user->dt*eM_2[r][0]*cv_sum;
959:       cols[1] = 5*idx[1];     vals[1] = user->dt*eM_2[r][1]*cv_sum;

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

964:       row = 5*idx[r]+2;
965:       cols[0] = 5*idx[0]+2;   vals[0] = user->dt*eM_2[r][0]*ci_sum;
966:       cols[1] = 5*idx[1]+2;   vals[1] = user->dt*eM_2[r][1]*ci_sum;

968:       /* Insert values in matrix M for 3nd dof */
969:       MatSetValuesLocal(M,1,&row,2,cols,vals,ADD_VALUES);
970:     }
971:   }
972: 
973:   /* old stuff
974:   VecRestoreArray(user->cv,&cv_p);
975:   VecRestoreArray(user->ci,&ci_p);
976:    */

978:   /* new stuff */
979:   VecRestoreArray(cvlocal,&cv_p);
980:   VecRestoreArray(cilocal,&ci_p);
981:   DMRestoreLocalVector(user->da2,&cvlocal);
982:   DMRestoreLocalVector(user->da2,&cilocal);
983: 
984:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
985:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
986:   DMDARestoreElements(user->da1,&nele,&nen,&ele);

988:   return(0);
989: }


994: PetscErrorCode CheckRedundancy(SNES snes, IS act, IS *outact, DM da)
995: {
997:   PetscScalar    **uin,**uout;
998:   Vec            UIN, UOUT;
999:   PetscInt       xs,xm,*outindex;
1000:   const PetscInt *index;
1001:   PetscInt       k,i,l,n,M,cnt=0;
1002: 
1004:   DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1005:   DMGetGlobalVector(da,&UIN);
1006:   VecSet(UIN,0.0);
1007:   DMGetLocalVector(da,&UOUT);
1008:   DMDAVecGetArrayDOF(da,UIN,&uin);
1009:   ISGetIndices(act,&index);
1010:   ISGetLocalSize(act,&n);
1011:   for (k=0;k<n;k++){
1012:     l = index[k]%5;
1013:     i = index[k]/5;
1014:     uin[i][l]=1.0;
1015:   }
1016:   printf("Number of active constraints before applying redundancy %d\n",n);
1017:   ISRestoreIndices(act,&index);
1018:   DMDAVecRestoreArrayDOF(da,UIN,&uin);
1019:   DMGlobalToLocalBegin(da,UIN,INSERT_VALUES,UOUT);
1020:   DMGlobalToLocalEnd(da,UIN,INSERT_VALUES,UOUT);
1021:   DMDAVecGetArrayDOF(da,UOUT,&uout);
1022: 
1023:   DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);

1025:   for(i=xs; i < xs+xm;i++) {
1026:     if (uout[i-1][1] && uout[i][1] && uout[i+1][1])
1027:              uout[i][0] = 1.0;
1028:     if (uout[i-1][3] && uout[i][3] && uout[i+1][3])
1029:              uout[i][2] = 1.0;
1030:   }

1032:   for(i=xs; i < xs+xm;i++) {
1033:     for(l=0;l<5;l++) {
1034:       if (uout[i][l])
1035:         cnt++;
1036:     }
1037:   }

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

1042:   PetscMalloc(cnt*sizeof(PetscInt),&outindex);
1043:   cnt = 0;
1044: 
1045:   for(i=xs; i < xs+xm;i++) {
1046:     for(l=0;l<5;l++) {
1047:       if (uout[i][l])
1048:         outindex[cnt++] = 5*(i)+l;
1049:     }
1050:   }
1051: 
1052: 
1053:   ISCreateGeneral(PETSC_COMM_WORLD,cnt,outindex,PETSC_OWN_POINTER,outact);
1054:   DMDAVecRestoreArrayDOF(da,UOUT,&uout);
1055:   DMRestoreGlobalVector(da,&UIN);
1056:   DMRestoreLocalVector(da,&UOUT);
1057:   return(0);
1058: }