Actual source code: ex633d_db.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  1: static char help[] = "3D 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: ./ex633d_db -ksp_type fgmres -snes_vi_monitor -snes_atol 1.e-11 -snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor -VG 10000000  -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_type basic -T .0020 -P_casc .0005 -snes_monitor_solution -da_refine 1
 13: ./ex633d_db -ksp_type fgmres -snes_vi_monitor -snes_atol 1.e-11 -snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor -VG 10000000  -pc_type mg -pc_mg_galerkin -log_summary  -da_refine 1

 15:  */

 17: #include <petscsnes.h>
 18: #include <petscdmda.h>

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

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

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


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


 70:   /* Get physics and time parameters */
 71:   GetParams(&user);

 73:   if (user.periodic) {
 74:     DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX, -3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, 5, 1,NULL,NULL,NULL,&user.da1);
 75:     DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX, -3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, 5, 1,NULL,NULL,NULL,&user.da1_clone);
 76:     DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX, -3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, 1, 1,NULL,NULL,NULL,&user.da2);

 78:   } else {
 79:     /* Create a 1D DA with dof = 5; the whole thing */
 80:     DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX, -3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, 5, 1,NULL,NULL,NULL,&user.da1);
 81:     DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX, -3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, 5, 1,NULL,NULL,NULL,&user.da1_clone);
 82:     /* Create a 1D DA with dof = 1; for individual componentes */
 83:     DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX, -3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, 1, 1,NULL,NULL,NULL,&user.da2);
 84:   }

 86:   /* Set Element type (rectangular) */
 87:   DMDASetElementType(user.da1,DMDA_ELEMENT_Q1);
 88:   DMDASetElementType(user.da2,DMDA_ELEMENT_Q1);

 90:   /* Set x and y coordinates */
 91:   DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,user.ymin,user.ymax, user.zmin,user.zmax);
 92:   DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,user.ymin,user.ymax, user.zmin,user.zmax);

 94:   /* Get global vector x from DM (da1) and duplicate vectors r,xl,xu */
 95:   DMCreateGlobalVector(user.da1,&x);
 96:   VecDuplicate(x,&r);
 97:   VecDuplicate(x,&xl);
 98:   VecDuplicate(x,&xu);
 99:   VecDuplicate(x,&user.q);


102:   /* Get global vector user->wv from da2 and duplicate other vectors */
103:   DMCreateGlobalVector(user.da2,&user.wv);
104:   VecDuplicate(user.wv,&user.cv);
105:   VecDuplicate(user.wv,&user.wi);
106:   VecDuplicate(user.wv,&user.ci);
107:   VecDuplicate(user.wv,&user.eta);
108:   VecDuplicate(user.wv,&user.cvi);
109:   VecDuplicate(user.wv,&user.DPsiv);
110:   VecDuplicate(user.wv,&user.DPsii);
111:   VecDuplicate(user.wv,&user.DPsieta);
112:   VecDuplicate(user.wv,&user.Pv);
113:   VecDuplicate(user.wv,&user.Pi);
114:   VecDuplicate(user.wv,&user.Piv);
115:   VecDuplicate(user.wv,&user.Riv);
116:   VecDuplicate(user.wv,&user.logcv);
117:   VecDuplicate(user.wv,&user.logci);
118:   VecDuplicate(user.wv,&user.logcvi);
119:   VecDuplicate(user.wv,&user.work1);
120:   VecDuplicate(user.wv,&user.work2);
121:   VecDuplicate(user.wv,&user.work3);
122:   VecDuplicate(user.wv,&user.work4);

124:   /* Get Jacobian matrix structure from the da for the entire thing, da1 */
125:   DMSetMatType(user.da1,MATAIJ);
126:   DMCreateMatrix(user.da1,&user.M);
127:   /* Get the (usual) mass matrix structure from da2 */
128:   DMSetMatType(user.da2,MATAIJ);
129:   DMCreateMatrix(user.da2,&user.M_0);

131:   SetInitialGuess(x,&user);

133:   /* Form the jacobian matrix and M_0 */
134:   SetUpMatrices(&user);


137:   MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);

139:   SNESCreate(PETSC_COMM_WORLD,&snes);
140:   SNESSetDM(snes,user.da1);

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

145:   SetVariableBounds(user.da1,xl,xu);
146:   /* SNESVISetRedundancyCheck(snes,(PetscErrorCode (*)(SNES,IS,IS*,void*))CheckRedundancy,user.da1_clone); */
147:   SNESVISetVariableBounds(snes,xl,xu);
148:   SNESSetFromOptions(snes);

150:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
151:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_p",FILE_MODE_WRITE,&view_p);
152:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
153:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);
154:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);


157:   /*
158:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
159:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);

161:   VecView(user.q,view_q);
162:   MatView(user.M,view_mat);

164:   PetscViewerDestroy(&view_q);
165:   PetscViewerDestroy(&view_mat);
166:  */



170:   while (t<user.T) {
171:     char        filename[PETSC_MAX_PATH_LEN];
172:     PetscScalar a = 1.0;
173:     PetscInt    i;
174:     PetscViewer view;


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

180:     SetRandomVectors(&user);
181:     DPsi(&user);
182:     VecView(user.DPsiv,view_psi);
183:     VecView(user.DPsii,view_psi);
184:     VecView(user.DPsieta,view_psi);

186:     VecView(user.Pv,view_p);
187:     Update_q(&user);
188:     VecView(user.q,view_q);

190:     printf("after VecView\n");
191:     MatView(user.M,view_mat);

193:     printf("after MatView\n");

195:     SNESSolve(snes,NULL,x);

197:     VecView(x,view_out);

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

203:     Update_u(x,&user);
204:     /*
205:     for (i=0; i < (int)(user.T/a) ; i++) {
206:       if (t/a > i - user.dt/a && t/a < i + user.dt/a) {
207:         sprintf(filename,"output_%f",t);
208:         PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&view);
209:         VecView(user.cv,view);
210:         VecView(user.ci,view);
211:         VecView(user.eta,view);
212:         PetscViewerDestroy(&view);
213:       }

215:     }
216:     */
217:     UpdateMatrices(&user);
218:     t    = t + user.dt;

220:   }

222:   PetscViewerDestroy(&view_out);
223:   PetscViewerDestroy(&view_p);
224:   PetscViewerDestroy(&view_q);
225:   PetscViewerDestroy(&view_psi);
226:   PetscViewerDestroy(&view_mat);
227:   VecDestroy(&x);
228:   VecDestroy(&r);
229:   VecDestroy(&xl);
230:   VecDestroy(&xu);
231:   VecDestroy(&user.q);
232:   VecDestroy(&user.wv);
233:   VecDestroy(&user.cv);
234:   VecDestroy(&user.wi);
235:   VecDestroy(&user.ci);
236:   VecDestroy(&user.eta);
237:   VecDestroy(&user.cvi);
238:   VecDestroy(&user.DPsiv);
239:   VecDestroy(&user.DPsii);
240:   VecDestroy(&user.DPsieta);
241:   VecDestroy(&user.Pv);
242:   VecDestroy(&user.Pi);
243:   VecDestroy(&user.Piv);
244:   VecDestroy(&user.Riv);
245:   VecDestroy(&user.logcv);
246:   VecDestroy(&user.logci);
247:   VecDestroy(&user.logcvi);
248:   VecDestroy(&user.work1);
249:   VecDestroy(&user.work2);
250:   VecDestroy(&user.work3);
251:   VecDestroy(&user.work4);
252:   MatDestroy(&user.M);
253:   MatDestroy(&user.M_0);
254:   DMDestroy(&user.da1);
255:   DMDestroy(&user.da1_clone);
256:   DMDestroy(&user.da2);
257:   SNESDestroy(&snes);

259:   printf("I am finalize \n");
260:   PetscFinalize();
261:   return 0;
262: }

266: PetscErrorCode Update_u(Vec X,AppCtx *user)
267: {
269:   PetscInt       i,n;
270:   PetscScalar    *xx,*wv_p,*cv_p,*wi_p,*ci_p,*eta_p;

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

275:   VecGetArray(X,&xx);
276:   VecGetArray(user->wv,&wv_p);
277:   VecGetArray(user->cv,&cv_p);
278:   VecGetArray(user->wi,&wi_p);
279:   VecGetArray(user->ci,&ci_p);
280:   VecGetArray(user->eta,&eta_p);


283:   for (i=0; i<n; i++) {
284:     wv_p[i]  = xx[5*i];
285:     cv_p[i]  = xx[5*i+1];
286:     wi_p[i]  = xx[5*i+2];
287:     ci_p[i]  = xx[5*i+3];
288:     eta_p[i] = xx[5*i+4];
289:   }
290:   VecRestoreArray(X,&xx);
291:   VecRestoreArray(user->wv,&wv_p);
292:   VecRestoreArray(user->cv,&cv_p);
293:   VecRestoreArray(user->wi,&wi_p);
294:   VecRestoreArray(user->ci,&ci_p);
295:   VecRestoreArray(user->eta,&eta_p);
296:   return(0);
297: }

301: PetscErrorCode Update_q(AppCtx *user)
302: {
304:   PetscScalar    *q_p, *w1, *w2;
305:   PetscInt       i,n;

308:   VecPointwiseMult(user->Riv,user->eta,user->eta); /* Riv = eta.^2                             */
309:   VecScale(user->Riv,user->Rsurf);                 /* Riv = Rsurf * eta.^2                     */
310:   VecShift(user->Riv,user->Rbulk);                 /* Riv = Rbulk + Rsurf * eta.^2             */
311:   VecPointwiseMult(user->Riv,user->ci,user->Riv);  /* Riv = (Rbulk + Rsurf * eta.^2) * ci      */
312:   VecPointwiseMult(user->Riv,user->cv,user->Riv);  /* Riv = (Rbulk + Rsurf * eta.^2) * ci * cv */

314:   VecCopy(user->Riv,user->work1);
315:   VecScale(user->work1,user->dt);             /* work1 = dt * Riv              */
316:   VecAXPY(user->work1,-1.0,user->cv);         /* work1 = dt * Riv - cv         */
317:   MatMult(user->M_0,user->work1,user->work2); /* work2 = M_0 * (dt * Riv - cv) */

319:   VecGetArray(user->q,&q_p);
320:   VecGetArray(user->work1,&w1);
321:   VecGetArray(user->work2,&w2);
322:   VecGetLocalSize(user->wv,&n);
323:   for (i=0; i<n; i++) q_p[5*i]=w2[i];

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

328:   VecCopy(user->Riv,user->work1);
329:   VecScale(user->work1,user->dt);
330:   VecAXPY(user->work1,-1.0,user->ci);
331:   MatMult(user->M_0,user->work1,user->work2);
332:   for (i=0; i<n; i++) q_p[5*i+2]=w2[i];

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

337:   VecCopy(user->DPsieta,user->work1);
338:   VecScale(user->work1,user->L);
339:   VecAXPY(user->work1,-1.0/user->dt,user->eta);
340:   MatMult(user->M_0,user->work1,user->work2);
341:   for (i=0; i<n; i++) q_p[5*i+4]=w2[i];

343:   VecRestoreArray(user->q,&q_p);
344:   VecRestoreArray(user->work1,&w1);
345:   VecRestoreArray(user->work2,&w2);
346:   return(0);
347: }

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

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

370:   Llog(user->cv,user->logcv);
371:   Llog(user->ci,user->logci);

373:   VecCopy(user->cv,user->cvi);
374:   VecAXPY(user->cvi,1.0,user->ci);
375:   VecScale(user->cvi,-1.0);
376:   VecShift(user->cvi,1.0);
377:   Llog(user->cvi,user->logcvi);

379:   for (i=0; i<n; i++)
380:   {
381:     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);
382:     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];

384:     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]);
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);
396:   return(0);
397: }


402: PetscErrorCode Llog(Vec X, Vec Y)
403: {
405:   PetscScalar    *x,*y;
406:   PetscInt       n,i;

409:   VecGetArray(X,&x);
410:   VecGetArray(Y,&y);
411:   VecGetLocalSize(X,&n);
412:   for (i=0; i<n; i++) {
413:     if (x[i] < 1.0e-12) y[i] = log(1.0e-12);
414:     else y[i] = log(x[i]);
415:   }
416:   return(0);
417: }

421: PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
422: {
423:   /******************** 3D ********************/
425:   PetscInt       n,i,j,Xda,Yda,Zda;
426:   PetscScalar    *xx,*cv_p,*ci_p,*wv_p,*wi_p,*eta_p;
427:   PetscViewer    view_out;

429:   /* needed for the void growth case */
430:   PetscScalar       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;
431:   PetscScalar       xmid,ymid,zmid;
432:   PetscInt          nele,nen,idx[8];
433:   const PetscInt    *ele;
434:   PetscScalar       x[8],y[8],z[8];
435:   Vec               coords;
436:   const PetscScalar *_coords;
437:   PetscViewer       view;
438:   PetscScalar       xwidth = user->xmax - user->xmin, ywidth = user->ymax - user->ymin, zwidth = user->zmax - user->zmin;

441:   VecGetLocalSize(X,&n);

443:   if (user->voidgrowth) {
444:     DMDAGetInfo(user->da2,NULL,&Xda,&Yda,&Zda,NULL,NULL,
445:                        NULL,NULL,NULL,NULL,NULL,NULL,NULL);
446:     DMGetCoordinatesLocal(user->da2,&coords);
447:     VecGetArrayRead(coords,&_coords);

449:     if (user->periodic) h = (user->xmax-user->xmin)/Xda;
450:     else h = (user->xmax-user->xmin)/(Xda-1.0);

452:     xmid   = (user->xmax + user->xmin)/2.0;
453:     ymid   = (user->ymax + user->ymin)/2.0;
454:     zmid   = (user->zmax + user->zmin)/2.0;
455:     lambda = 4.0*h;

457:     DMDAGetElements(user->da2,&nele,&nen,&ele); /* number of local elements, number of element nodes, the local indices of the elements' vertices */
458:     for (i=0; i < nele; i++) {
459:       for (j = 0; j<8; j++) {
460:         idx[j] = ele[8*i + j];
461:         x[j]   = _coords[3*idx[j]];
462:         y[j]   = _coords[3*idx[j]+1];
463:         z[j]   = _coords[3*idx[j]+2];
464:       }

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

488:       VecSetValuesLocal(user->cv,8,idx,vals_cv,INSERT_VALUES);
489:       VecSetValuesLocal(user->ci,8,idx,vals_ci,INSERT_VALUES);
490:       VecSetValuesLocal(user->eta,8,idx,vals_eta,INSERT_VALUES);
491:     }
492:     VecAssemblyBegin(user->cv);
493:     VecAssemblyEnd(user->cv);
494:     VecAssemblyBegin(user->ci);
495:     VecAssemblyEnd(user->ci);
496:     VecAssemblyBegin(user->eta);
497:     VecAssemblyEnd(user->eta);

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

502:   } else {

504:     VecSet(user->cv,6.9e-4);
505:     VecSet(user->ci,6.9e-4);
506:     VecSet(user->eta,0.0);
507:   }

509:   DPsi(user);
510:   VecCopy(user->DPsiv,user->wv);
511:   VecCopy(user->DPsii,user->wi);

513:   VecGetArray(X,&xx);
514:   VecGetArray(user->cv,&cv_p);
515:   VecGetArray(user->ci,&ci_p);
516:   VecGetArray(user->eta,&eta_p);
517:   VecGetArray(user->wv,&wv_p);
518:   VecGetArray(user->wi,&wi_p);
519:   for (i=0; i<n/5; i++) {
520:     xx[5*i]  =wv_p[i];
521:     xx[5*i+1]=cv_p[i];
522:     xx[5*i+2]=wi_p[i];
523:     xx[5*i+3]=ci_p[i];
524:     xx[5*i+4]=eta_p[i];
525:   }

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

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

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

552:   /* static PetscViewer viewer=0; */
553:   static PetscRandom rand = 0;
554:   static PetscInt    step = 0;

557:   if (!rand) {
558:     PetscRandomCreate(PETSC_COMM_WORLD,&rand);
559:     PetscRandomSetFromOptions(rand);
560:   }

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

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

591: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
592: {
594:   AppCtx         *user=(AppCtx*)ctx;

597:   MatMultAdd(user->M,X,user->q,F);
598:   return(0);
599: }

603: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat B,void *ctx)
604: {
606:   AppCtx         *user=(AppCtx*)ctx;

609:   MatCopy(user->M,J,*flg);
610:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
611:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
612:   return(0);
613: }

617: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
618: {
620:   PetscScalar    ****l,****u;
621:   PetscInt       xs,xm,ys,ym,zs,zm;
622:   PetscInt       k,j,i;

625:   DMDAVecGetArrayDOF(da,xl,&l);
626:   DMDAVecGetArrayDOF(da,xu,&u);

628:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);

630:   for (k = zs; k < zs + zm; k++) {
631:     for (j=ys; j<ys+ym; j++) {
632:       for (i=xs; i < xs+xm; i++) {
633:         l[k][j][i][0] = -PETSC_INFINITY;
634:         l[k][j][i][1] = 0.0;
635:         l[k][j][i][2] = -PETSC_INFINITY;
636:         l[k][j][i][3] = 0.0;
637:         l[k][j][i][4] = 0.0;
638:         u[k][j][i][0] = PETSC_INFINITY;
639:         u[k][j][i][1] = 1.0;
640:         u[k][j][i][2] = PETSC_INFINITY;
641:         u[k][j][i][3] = 1.0;
642:         u[k][j][i][4] = 1.0;
643:       }
644:     }
645:   }

647:   DMDAVecRestoreArrayDOF(da,xl,&l);
648:   DMDAVecRestoreArrayDOF(da,xu,&u);
649:   return(0);
650: }

654: PetscErrorCode GetParams(AppCtx *user)
655: {
657:   PetscBool      flg;

660:   /* Set default parameters */
661:   user->xmin    = 0.0; user->xmax = 64.0;
662:   user->ymin    = 0.0; user->ymax = 64.0;
663:   user->zmin    = 0.0; user->zmax = 64.0;
664:   user->Dv      = 1.0; user->Di=1.0;
665:   user->Evf     = 0.8; user->Eif = 0.8;
666:   user->A       = 1.0;
667:   user->kBT     = 0.11;
668:   user->kav     = 1.0; user->kai = 1.0; user->kaeta = 1.0;
669:   user->Rsurf   = 10.0; user->Rbulk = 0.0;
670:   user->L       = 10.0; user->P_casc = 0.05;
671:   user->T       = 1.0e-2;    user->dt = 1.0e-4;
672:   user->VG      = 100.0;
673:   user->initv   = .0001;
674:   user->initeta = 0.0;

676:   user->degenerate = PETSC_FALSE;
677:   /* void growth */
678:   user->voidgrowth = PETSC_FALSE;

680:   PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);
681:   PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);
682:   PetscOptionsGetReal(NULL,"-T",&user->T,&flg);
683:   PetscOptionsGetReal(NULL,"-dt",&user->dt,&flg);
684:   PetscOptionsBool("-degenerate","Run with degenerate mobility\n","None",user->degenerate,&user->degenerate,&flg);
685:   PetscOptionsReal("-smallnumber","Small number added to degenerate mobility\n","None",user->smallnumber,&user->smallnumber,&flg);
686:   PetscOptionsBool("-voidgrowth","Use initial conditions for void growth\n","None",user->voidgrowth,&user->voidgrowth,&flg);
687:   return(0);
688: }


693: PetscErrorCode SetUpMatrices(AppCtx *user)
694: {
696:   PetscInt       nele,nen,i,j,n;
697:   const PetscInt *ele;
698:   PetscScalar    dt=user->dt,hx,hy,hz;

700:   PetscInt    idx[8];
701:   PetscScalar eM_0[8][8],eM_2[8][8];
702:   PetscScalar cv_sum, ci_sum;
703:   PetscScalar tp_cv, tp_ci;

705:   Mat         M  =user->M;
706:   Mat         M_0=user->M_0;
707:   PetscInt    Xda,Yda,Zda;
708:   PetscScalar *cv_p,*ci_p;
709:   Vec         cvlocal,cilocal;

712:   DMGetLocalVector(user->da2,&cvlocal);
713:   DMGetLocalVector(user->da2,&cilocal);
714:   DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
715:   DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
716:   DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
717:   DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);

719:   VecGetArray(cvlocal,&cv_p);
720:   VecGetArray(cilocal,&ci_p);

722:   MatGetLocalSize(M,&n,NULL);
723:   DMDAGetInfo(user->da1,NULL,&Xda,&Yda,&Zda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

725:   if (Xda!=Yda) printf("Currently different Xda and Yda are not supported");
726:   if (user->periodic) {
727:     hx = (user->xmax-user->xmin)/Xda;
728:     hy = (user->ymax-user->ymin)/Yda;
729:     hz = (user->zmax-user->zmin)/Zda;
730:   } else {
731:     hx = (user->xmax-user->xmin)/(Xda-1.0);
732:     hy = (user->ymax-user->ymin)/(Yda-1.0);
733:     hz = (user->zmax-user->zmin)/(Zda-1.0);
734:   }

736:   /** blocks of M_0 **/
737:   for (i = 0; i<8; i++) eM_0[i][i] = 1;

739:   eM_0[0][1] = 0.5; eM_0[0][2] = 0.25; eM_0[0][3] = 0.5; eM_0[0][4] = 0.5; eM_0[0][5] = 0.25; eM_0[0][6] = 0.125; eM_0[0][7] = 0.25;
740:   eM_0[1][2] = 0.5; eM_0[1][3] = 0.25; eM_0[1][4] = 0.25;eM_0[1][5] = 0.5; eM_0[1][6] = 0.25; eM_0[1][7] = 0.125;
741:   eM_0[2][3] = 0.5; eM_0[2][4] = 0.125;eM_0[2][5] = 0.25;eM_0[2][6] = 0.5;eM_0[2][7] = 0.25;
742:   eM_0[3][4] = 0.25;eM_0[3][5] = 0.125;eM_0[3][6] = 0.25;eM_0[3][7] = 0.5;
743:   eM_0[4][5] = 0.5; eM_0[4][6] = 0.25; eM_0[4][7] = 0.5;
744:   eM_0[5][6] = 0.5; eM_0[5][7] = 0.25;
745:   eM_0[6][7] = 0.5;

747:   for (i = 0; i<8; i++) {
748:     for (j =0; j<8; j++) {
749:       if (i>j) eM_0[i][j] = eM_0[j][i];
750:     }
751:   }
752: 
753:   for (i = 0; i<8; i++) {
754:     for (j =0; j<8; j++) {
755:       eM_0[i][j] = hx*hy*hz/27.0 * eM_0[i][j];
756:     }
757:   }

759:   /** blocks of M_2 **/

761:   for (i = 0; i<8; i++) eM_2[i][i] = 12;

763:   eM_2[0][1] = 0; eM_2[0][2] = -3; eM_2[0][3] = 0;  eM_2[0][4] = 0;  eM_2[0][5] = -3; eM_2[0][6] = -3; eM_2[0][7] = -3;
764:   eM_2[1][2] = 0; eM_2[1][3] = -3; eM_2[1][4] = -3; eM_2[1][5] = 0;  eM_2[1][6] = -3; eM_2[1][7] = -3;
765:   eM_2[2][3] = 0; eM_2[2][4] = -3; eM_2[2][5] = -3; eM_2[2][6] = 0;  eM_2[2][7] = -3;
766:   eM_2[3][4] = -3;eM_2[3][5] = -3; eM_2[3][6] = -3; eM_2[3][7] = 0;
767:   eM_2[4][5] = 0; eM_2[4][6] = -3; eM_2[4][7] = 0;
768:   eM_2[5][6] = 0; eM_2[5][7] = -3;
769:   eM_2[6][7] = 0;

771:   for (i = 0; i<8; i++) {
772:     for (j =0; j<8; j++) {
773:       if (i>j) eM_2[i][j] = eM_2[j][i];
774:     }
775:   }

777:   for (i = 0; i<8; i++) {
778:     for (j =0; j<8; j++) {
779:       eM_2[i][j] = hx*hy*hz/36.0 * eM_2[i][j];
780:     }
781:   }

783:   /* Get local element info */
784:   DMDAGetElements(user->da1,&nele,&nen,&ele);
785:   PetscInt    row,cols[16],r,row_M_0,cols3[8];
786:   PetscScalar vals[16],vals_M_0[8],vals3[8];

788:   for (i=0;i < nele;i++) {
789:     for (j=0; j<8; j++) idx[j] = ele[8*i + j];

791:     for (r=0;r<8;r++) {
792:       row_M_0 = idx[r];
793:       for (j=0; j<8; j++) vals_M_0[j]=eM_0[r][j];

795:       MatSetValuesLocal(M_0,1,&row_M_0,8,idx,vals_M_0,ADD_VALUES);

797:       tp_cv = 0.0;
798:       tp_ci = 0.0;
799:       for (j = 0; j < 8; j++) {
800:         tp_cv = tp_cv + cv_p[idx[j]];
801:         tp_ci = tp_ci + ci_p[idx[j]];
802:       }

804:       cv_sum = tp_cv * user->Dv/user->kBT;
805:       ci_sum = tp_ci * user->Di/user->kBT;

807:       row = 5*idx[r];
808:       for (j = 0; j < 8; j++) {
809:         cols[j]     = 5 * idx[j];
810:         vals[j]     = dt*eM_2[r][j]*cv_sum;
811:         cols[j + 8] = 5 * idx[j] + 1;
812:         vals[j + 8] = eM_0[r][j];
813:       }

815:       /* Insert values in matrix M for 1st dof */

817:       MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);

819:       row = 5*idx[r]+1;
820:       for (j = 0; j < 8; j++) {
821:         cols[j]     = 5 * idx[j];
822:         vals[j]     = -eM_0[r][j];
823:         cols[j + 8] = 5 * idx[j] + 1;
824:         vals[j + 8] = user->kav*eM_2[r][j];
825:       }

827:       /* Insert values in matrix M for 2nd dof */
828:       MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);

830:       row = 5*idx[r]+2;
831:       for (j = 0; j < 8; j++) {
832:         cols[j]     = 5 * idx[j] + 2;
833:         vals[j]     = dt*eM_2[r][j]*ci_sum;
834:         cols[j + 8] = 5 * idx[j] + 3;
835:         vals[j + 8] = eM_0[r][j];
836:       }

838:       /* Insert values in matrix M for 3rd dof */
839:       MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);

841:       row = 5*idx[r]+3;
842:       for (j = 0; j < 8; j++) {
843:         cols[j]     = 5 * idx[j] + 2;
844:         vals[j]     = -eM_0[r][j];
845:         cols[j + 8] = 5 * idx[j] + 3;
846:         vals[j + 8] = user->kai*eM_2[r][j];
847:       }

849:       /* Insert values in matrix M for 4th dof */
850:       MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);

852:       row = 5*idx[r]+4;
853:       for (j = 0; j < 8; j++) {
854:         cols3[j] = 5 * idx[j] + 4;
855:         vals3[j] = eM_0[r][j]/dt + user->L*user->kaeta*eM_2[r][j];
856:       }

858:       /* Insert values in matrix M for 5th dof */
859:       MatSetValuesLocal(M,1,&row,8,cols3,vals3,ADD_VALUES);
860:     }
861:   }


864:   VecRestoreArray(cvlocal,&cv_p);
865:   VecRestoreArray(cilocal,&ci_p);

867:   DMRestoreLocalVector(user->da2,&cvlocal);
868:   DMRestoreLocalVector(user->da2,&cilocal);

870:   MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
871:   MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);

873:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
874:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

876:   DMDARestoreElements(user->da1,&nele,&nen,&ele);
877:   return(0);
878: }


883: PetscErrorCode UpdateMatrices(AppCtx *user)
884: {
886:   PetscInt       nele,nen,i,j,n,Xda,Yda,Zda;
887:   const PetscInt *ele;

889:   PetscInt    idx[8];
890:   PetscScalar eM_2[8][8],h;
891:   Mat         M=user->M;
892:   PetscScalar *cv_p,*ci_p,cv_sum,ci_sum;
893:   /* newly added */
894:   Vec cvlocal,cilocal;

897:   DMGetLocalVector(user->da2,&cvlocal);
898:   DMGetLocalVector(user->da2,&cilocal);
899:   DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
900:   DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
901:   DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
902:   DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);

904:   VecGetArray(cvlocal,&cv_p);
905:   VecGetArray(cilocal,&ci_p);

907:   /* Create the mass matrix M_0 */
908:   MatGetLocalSize(M,&n,NULL);
909:   DMDAGetInfo(user->da1,NULL,&Xda,&Yda,&Zda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

911:   if (Xda!=Yda) printf("Currently different Xda and Yda are not supported");
912:   if (user->periodic) h = (user->xmax-user->xmin)/Xda;
913:   else h = (user->xmax-user->xmin)/(Xda-1.0);

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

918:   for (i=0; i<nele; i++) {
919:     for (j=0; j<8; j++) idx[j] = ele[8*i + j];

921:     PetscInt    r,row,cols[8];
922:     PetscScalar vals[8];

924:     for (r=0; r<8; r++) {
925:       row = 5*idx[r];
926:       for (j = 0; j < 8; j++) {
927:         cols[j] = 5*idx[j];
928:         vals[j] = 0.0;
929:       }

931:       /* Insert values in matrix M for 1st dof */
932:       MatSetValuesLocal(M,1,&row,8,cols,vals,INSERT_VALUES);

934:       row = 5*idx[r]+2;
935:       for (j = 0; j < 8; j++) {
936:         cols[j] = 5*idx[j]+2;
937:         vals[j] = 0.0;
938:       }

940:       /* Insert values in matrix M for 3nd dof */
941:       MatSetValuesLocal(M,1,&row,8,cols,vals,INSERT_VALUES);
942:     }
943:   }

945:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
946:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

948:   /** blocks of M_2 **/

950:   for (i = 0; i<8; i++) eM_2[i][i] = 12;

952:   eM_2[0][1] = 0; eM_2[0][2] = -3; eM_2[0][3] = 0;  eM_2[0][4] = 0;  eM_2[0][5] = -3; eM_2[0][6] = -3; eM_2[0][7] = -3;
953:   eM_2[1][2] = 0; eM_2[1][3] = -3; eM_2[1][4] = -3; eM_2[1][5] = 0;  eM_2[1][6] = -3; eM_2[1][7] = -3;
954:   eM_2[2][3] = 0; eM_2[2][4] = -3; eM_2[2][5] = -3; eM_2[2][6] = 0;  eM_2[2][7] = -3;
955:   eM_2[3][4] = -3;eM_2[3][5] = -3; eM_2[3][6] = -3; eM_2[3][7] = 0;
956:   eM_2[4][5] = 0; eM_2[4][6] = -3; eM_2[4][7] = 0;
957:   eM_2[5][6] = 0; eM_2[5][7] = -3;
958:   eM_2[6][7] = 0;

960:   for (i = 0; i<8; i++) {
961:     for (j =0; j<8; j++) {
962:       if (i>j) eM_2[i][j] = eM_2[j][i];
963:     }
964:   }

966:   for (i = 0; i<8; i++) {
967:     for (j =0; j<8; j++) {
968:       eM_2[i][j] = h*h*h/36.0 * eM_2[i][j];
969:     }
970:   }

972:   for (i=0;i<nele;i++) {
973:     for (j=0; j<8; j++) idx[j] = ele[8*i + j];

975:     PetscInt    row,cols[8],r;
976:     PetscScalar vals[8];

978:     for (r=0;r<8;r++) {

980:       if (user->degenerate) {
981:         cv_sum = (2.0*user->smallnumber + cv_p[idx[0]] + cv_p[idx[1]])*user->Dv/(2.0*user->kBT);
982:         ci_sum = (2.0*user->smallnumber + ci_p[idx[0]] + ci_p[idx[1]])*user->Di/(2.0*user->kBT);
983:       } else {
984:         cv_sum = user->initv*user->Dv/user->kBT;
985:         ci_sum = user->initv*user->Di/user->kBT;
986:       }

988:       row = 5*idx[r];
989:       for (j=0; j<8; j++) {
990:         cols[j] = 5*idx[j];
991:         vals[j] = user->dt*eM_2[r][j]*cv_sum;
992:       }

994:       /* Insert values in matrix M for 1st dof */
995:       MatSetValuesLocal(M,1,&row,8,cols,vals,ADD_VALUES);

997:       row = 5*idx[r]+2;
998:       for (j=0; j<8; j++) {
999:         cols[j] = 5*idx[j]+2;
1000:         vals[j] = user->dt*eM_2[r][j]*ci_sum;
1001:       }

1003:       /* Insert values in matrix M for 3nd dof */
1004:       MatSetValuesLocal(M,1,&row,8,cols,vals,ADD_VALUES);
1005:     }
1006:   }

1008:   VecRestoreArray(cvlocal,&cv_p);
1009:   VecRestoreArray(cilocal,&ci_p);
1010:   DMRestoreLocalVector(user->da2,&cvlocal);
1011:   DMRestoreLocalVector(user->da2,&cilocal);

1013:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1014:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
1015:   DMDARestoreElements(user->da1,&nele,&nen,&ele);
1016:   return(0);
1017: }