Actual source code: ex65.c

petsc-3.4.4 2014-03-13
  1: static char help[] = "2D coupled Allen-Cahn and Cahn-Hilliard equation for constant mobility. Only c_v and eta are considered.\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: ./ex65 -ksp_type fgmres  -snes_atol 1.e-13  -da_refine 6  -VG 10   -pc_type mg -pc_mg_galerkin -log_summary -dt .000001 -mg_coarse_pc_type redundant -mg_coarse_redundant_pc_type svd  -ksp_rtol 1.e-13 -snes_linesearch_type basic -T .0020  -voidgrowth -da_refine 5 -draw_fileds 0,1,2 -dt .0001 -T 1 -da_grid_x 4 -da_grid_y 4 -periodic -snes_rtol 1.e-13 -ksp_atol 1.e-13  -snes_vi_ignore_function_sign -domain 1
 13:  */

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

 18: typedef struct {
 19:   PetscReal   dt,T; /* Time step and end time */
 20:   PetscReal   dtevent;  /* time scale of radiation events, roughly one event per dtevent */
 21:   PetscInt    maxevents; /* once this number of events is reached no more events are generated */
 22:   PetscReal   initv;
 23:   PetscReal   initeta;
 24:   DM          da1,da1_clone,da2;
 25:   Mat         M;    /* Jacobian matrix */
 26:   Mat         M_0;
 27:   Vec         q,wv,cv,eta,DPsiv,DPsieta,logcv,logcv2,Pv,Pi,Piv;
 28:   Vec         phi1,phi2,Phi2D_V,Sv,Si; /* for twodomain modeling */
 29:   Vec         work1,work2;
 30:   PetscScalar Mv,L,kaeta,kav,Evf,A,B,cv0,Sv_scalar,VG;
 31:   PetscScalar Svr,Sir,cv_eq,ci_eq; /* for twodomain modeling */
 32:   Vec         work3,work4; /* for twodomain modeling*/
 33:   PetscReal   xmin,xmax,ymin,ymax;
 34:   PetscInt    nx;
 35:   PetscBool   graphics;
 36:   PetscBool   periodic;
 37:   PetscBool   lumpedmass;
 38:   PetscBool   radiation; /* Either radiation or void growth */
 39:   PetscInt    domain;
 40:   PetscReal   grain; /* some bogus factor that controls the "strength" of the grain boundaries */
 41:   PetscInt    darefine;
 42:   PetscInt    dagrid;
 43: } AppCtx;

 45: PetscErrorCode GetParams(AppCtx*);
 46: PetscErrorCode SetRandomVectors(AppCtx*,PetscReal);
 47: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
 48: PetscErrorCode SetUpMatrices(AppCtx*);
 49: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 50: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 51: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
 52: PetscErrorCode Update_q(AppCtx*);
 53: PetscErrorCode Update_u(Vec,AppCtx*);
 54: PetscErrorCode DPsi(AppCtx*);
 55: PetscErrorCode Llog(Vec,Vec);
 56: PetscErrorCode CheckRedundancy(SNES,IS,IS*,DM);
 57: PetscErrorCode Phi(AppCtx*);
 58: PetscErrorCode Phi_read(AppCtx*);

 62: int main(int argc, char **argv)
 63: {
 65:   Vec            x,r;  /* solution and residual vectors */
 66:   SNES           snes; /* Nonlinear solver context */
 67:   AppCtx         user; /* Application context */
 68:   Vec            xl,xu; /* Upper and lower bounds on variables */
 69:   Mat            J;
 70:   PetscScalar    t=0.0;
 71:   /* PetscViewer    view_out, view_p, view_q, view_psi, view_mat; */
 72:   PetscReal      bounds[] = {1000.0,-1000.,0.0,1.0,1000.0,-1000.0,0.0,1.0,1000.0,-1000.0};


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

 77:   /* Get physics and time parameters */
 78:   GetParams(&user);

 80:   if (user.periodic) {
 81:     DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE, 3, 1,NULL,NULL,&user.da1);
 82:     DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE, 3, 1,NULL,NULL,&user.da1_clone);
 83:     DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE, 1, 1,NULL,NULL,&user.da2);
 84:   } else {
 85:     DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE, 3, 1,NULL,NULL,&user.da1);
 86:     DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE, 3, 1,NULL,NULL,&user.da1_clone);
 87:     DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE, 1, 1,NULL,NULL,&user.da2);
 88:   }
 89:   /* Set Element type (triangular) */
 90:   DMDASetElementType(user.da1,DMDA_ELEMENT_P1);
 91:   DMDASetElementType(user.da2,DMDA_ELEMENT_P1);

 93:   /* Set x and y coordinates */
 94:   DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,user.ymin,user.ymax,NULL,NULL);
 95:   DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,user.ymin,user.ymax,NULL,NULL);
 96:   /* Get global vector x from DM (da1) and duplicate vectors r,xl,xu */
 97:   DMCreateGlobalVector(user.da1,&x);
 98:   VecDuplicate(x,&r);
 99:   VecDuplicate(x,&xl);
100:   VecDuplicate(x,&xu);
101:   VecDuplicate(x,&user.q);

103:   /* Get global vector user->wv from da2 and duplicate other vectors */
104:   DMCreateGlobalVector(user.da2,&user.wv);
105:   VecDuplicate(user.wv,&user.cv);
106:   VecDuplicate(user.wv,&user.eta);
107:   VecDuplicate(user.wv,&user.DPsiv);
108:   VecDuplicate(user.wv,&user.DPsieta);
109:   VecDuplicate(user.wv,&user.Pv);
110:   VecDuplicate(user.wv,&user.Pi);
111:   VecDuplicate(user.wv,&user.Piv);
112:   VecDuplicate(user.wv,&user.logcv);
113:   VecDuplicate(user.wv,&user.logcv2);
114:   VecDuplicate(user.wv,&user.work1);
115:   VecDuplicate(user.wv,&user.work2);
116:   /* for twodomain modeling */
117:   VecDuplicate(user.wv,&user.phi1);
118:   VecDuplicate(user.wv,&user.phi2);
119:   VecDuplicate(user.wv,&user.Phi2D_V);
120:   VecDuplicate(user.wv,&user.Sv);
121:   VecDuplicate(user.wv,&user.Si);
122:   VecDuplicate(user.wv,&user.work3);
123:   VecDuplicate(user.wv,&user.work4);

125:   /* Get Jacobian matrix structure from the da for the entire thing, da1 */
126:   DMCreateMatrix(user.da1,MATAIJ,&user.M);
127:   /* Get the (usual) mass matrix structure from da2 */
128:   DMCreateMatrix(user.da2,MATAIJ,&user.M_0);
129:   /* Form the jacobian matrix and M_0 */
130:   SetUpMatrices(&user);
131:   SetInitialGuess(x,&user);
132:   MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);

134:   SNESCreate(PETSC_COMM_WORLD,&snes);
135:   SNESSetDM(snes,user.da1);

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

140:   SetVariableBounds(user.da1,xl,xu);
141:   SNESVISetVariableBounds(snes,xl,xu);
142:   SNESSetFromOptions(snes);
143:   /* SNESVISetRedundancyCheck(snes,(PetscErrorCode (*)(SNES,IS,IS*,void*))CheckRedundancy,user.da1_clone); */
144:   /*
145:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
146:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_p",FILE_MODE_WRITE,&view_p);
147:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
148:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);
149:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);
150:    */
151:   /* PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds); */


154:   if (user.graphics) {
155:     /*PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds);*/

157:     VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
158:   }


161:   /* multidomain modeling */
162:   if (user.domain) {
163:     switch (user.domain) {
164:     case 1:
165:       Phi(&user);
166:       break;
167:     case 2:
168:       Phi_read(&user);
169:       break;
170:     }
171:   }

173:   while (t<user.T) {

175:     char        filename[PETSC_MAX_PATH_LEN];
176:     PetscScalar a = 1.0;
177:     PetscInt    i;
178:     /*PetscViewer  view;*/


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

184:     if (user.radiation) {
185:       SetRandomVectors(&user,t);
186:     }
187:     DPsi(&user);
188:     /*VecView(user.DPsiv,view_psi);*/
189:     /*VecView(user.DPsieta,view_psi);*/

191:     Update_q(&user);
192:     /*VecView(user.q,view_q);*/
193:     /*MatView(user.M,view_mat);*/

195:     SNESSolve(snes,NULL,x);
196:     /*VecView(x,view_out);*/


199:     if (user.graphics) {
200:       VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
201:     }

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

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

220:     t = t + user.dt;

222:   }

224:   /*
225:   PetscViewerDestroy(&view_out);
226:   PetscViewerDestroy(&view_p);
227:   PetscViewerDestroy(&view_q);
228:   PetscViewerDestroy(&view_psi);
229:   PetscViewerDestroy(&view_mat);
230:    */
231:   VecDestroy(&x);
232:   VecDestroy(&r);
233:   VecDestroy(&xl);
234:   VecDestroy(&xu);
235:   VecDestroy(&user.q);
236:   VecDestroy(&user.wv);
237:   VecDestroy(&user.cv);
238:   VecDestroy(&user.eta);
239:   VecDestroy(&user.DPsiv);
240:   VecDestroy(&user.DPsieta);
241:   VecDestroy(&user.Pv);
242:   VecDestroy(&user.Pi);
243:   VecDestroy(&user.Piv);
244:   VecDestroy(&user.logcv);
245:   VecDestroy(&user.logcv2);
246:   VecDestroy(&user.work1);
247:   VecDestroy(&user.work2);
248:   /* for twodomain modeling */
249:   VecDestroy(&user.phi1);
250:   VecDestroy(&user.phi2);
251:   VecDestroy(&user.Phi2D_V);
252:   VecDestroy(&user.Sv);
253:   VecDestroy(&user.Si);
254:   VecDestroy(&user.work3);
255:   VecDestroy(&user.work4);

257:   MatDestroy(&user.M);
258:   MatDestroy(&user.M_0);
259:   DMDestroy(&user.da1);
260:   DMDestroy(&user.da1_clone);
261:   DMDestroy(&user.da2);
262:   SNESDestroy(&snes);
263:   PetscFinalize();
264:   return 0;
265: }

269: PetscErrorCode Update_u(Vec X,AppCtx *user)
270: {
272:   PetscInt       i,n;
273:   PetscScalar    *xx,*wv_p,*cv_p,*eta_p;

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

278:   VecGetArray(X,&xx);
279:   VecGetArray(user->wv,&wv_p);
280:   VecGetArray(user->cv,&cv_p);
281:   VecGetArray(user->eta,&eta_p);


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

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

305:   VecGetArray(user->q,&q_p);
306:   VecGetArray(user->work1,&w1);
307:   VecGetArray(user->work2,&w2);
308:   VecGetLocalSize(user->wv,&n);


311:   VecSet(user->work1,0.0);
312:   if (user->radiation) {
313:     VecAXPY(user->work1,-1.0,user->Pv);
314:   }
315:   if (user->domain) {
316:     VecCopy(user->cv,user->work3);
317:     VecShift(user->work3,-1.0*user->cv_eq);
318:     VecCopy(user->Phi2D_V,user->work4);
319:     VecScale(user->work4,-1.0);
320:     VecShift(user->work4,1.0);
321:     VecPointwiseMult(user->work4,user->work4,user->work3);
322:     VecScale(user->work4,user->Svr);
323:     /* Parameter tuning: user->grain */
324:     /* 5000.0 worked for refinement 2 and 5 */
325:     /*10000.0 worked for refinement 2, 3, 4, 5 */
326:     VecAXPY(user->work1,user->grain,user->work4);
327:   }
328:   VecScale(user->work1,user->dt);
329:   VecAXPY(user->work1,-1.0,user->cv);
330:   MatMult(user->M_0,user->work1,user->work2);
331:   for (i=0; i<n; i++) q_p[3*i]=w2[i];

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

336:   VecCopy(user->DPsieta,user->work1);
337:   VecScale(user->work1,user->L*user->dt);
338:   VecAXPY(user->work1,-1.0,user->eta);
339:   if (user->radiation) {
340:     VecAXPY(user->work1,-1.0,user->Piv);
341:   }
342:   MatMult(user->M_0,user->work1,user->work2);
343:   for (i=0; i<n; i++) q_p[3*i+2]=w2[i];

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

353: PetscErrorCode DPsi(AppCtx *user)
354: {
356:   PetscScalar    Evf=user->Evf,A=user->A,B=user->B,cv0=user->cv0;
357:   PetscScalar    *cv_p,*eta_p,*logcv_p,*logcv2_p,*DPsiv_p,*DPsieta_p;
358:   PetscInt       n,i;

361:   VecGetLocalSize(user->cv,&n);
362:   VecGetArray(user->cv,&cv_p);
363:   VecGetArray(user->eta,&eta_p);
364:   VecGetArray(user->logcv,&logcv_p);
365:   VecGetArray(user->logcv2,&logcv2_p);
366:   VecGetArray(user->DPsiv,&DPsiv_p);
367:   VecGetArray(user->DPsieta,&DPsieta_p);

369:   Llog(user->cv,user->logcv);

371:   VecCopy(user->cv,user->work1);
372:   VecScale(user->work1,-1.0);
373:   VecShift(user->work1,1.0);
374:   Llog(user->work1,user->logcv2);

376:   for (i=0; i<n; i++)
377:   {
378:     DPsiv_p[i] = (eta_p[i]-1.0)*(eta_p[i]-1.0)*(eta_p[i]+1.0)*(eta_p[i]+1.0)*(Evf + logcv_p[i] - logcv2_p[i]) - 2.0*A*(cv_p[i] - cv0)*eta_p[i]*(eta_p[i]+2.0)*(eta_p[i]-1.0)*(eta_p[i]-1.0) + 2.0*B*(cv_p[i] - 1.0)*eta_p[i]*eta_p[i];

380:     DPsieta_p[i] = 4.0*eta_p[i]*(eta_p[i]-1.0)*(eta_p[i]+1.0)*(Evf*cv_p[i] + cv_p[i]*logcv_p[i] + (1.0-cv_p[i])*logcv2_p[i]) - A*(cv_p[i] - cv0)*(cv_p[i] - cv0)*(4.0*eta_p[i]*eta_p[i]*eta_p[i] - 6.0*eta_p[i] + 2.0) + 2.0*B*(cv_p[i]-1.0)*(cv_p[i]-1.0)*eta_p[i];

382:   }

384:   VecRestoreArray(user->cv,&cv_p);
385:   VecRestoreArray(user->eta,&eta_p);
386:   VecGetArray(user->logcv,&logcv_p);
387:   VecGetArray(user->logcv2,&logcv2_p);
388:   VecRestoreArray(user->DPsiv,&DPsiv_p);
389:   VecRestoreArray(user->DPsieta,&DPsieta_p);
390:   return(0);
391: }


396: PetscErrorCode Llog(Vec X, Vec Y)
397: {
399:   PetscScalar    *x,*y;
400:   PetscInt       n,i;

403:   VecGetArray(X,&x);
404:   VecGetArray(Y,&y);
405:   VecGetLocalSize(X,&n);
406:   for (i=0; i<n; i++) {
407:     if (x[i] < 1.0e-12) y[i] = log(1.0e-12);
408:     else y[i] = log(x[i]);
409:   }
410:   return(0);
411: }

415: PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
416: {
417:   PetscErrorCode    ierr;
418:   PetscInt          n,i,Mda,Nda;
419:   PetscScalar       *xx,*cv_p,*wv_p,*eta_p;
420:   /*PetscViewer      view_out;*/
421:   /* needed for the void growth case */
422:   PetscScalar       xmid,ymid,cv_v=1.0,cv_m=user->Sv_scalar*user->cv0,eta_v=1.0,eta_m=0.0,h,lambda;
423:   PetscInt          nele,nen,idx[3];
424:   const PetscInt    *ele;
425:   PetscScalar       x[3],y[3];
426:   Vec               coords;
427:   const PetscScalar *_coords;
428:   PetscScalar       xwidth = user->xmax - user->xmin, ywidth = user->ymax - user->ymin;

431:   VecGetLocalSize(X,&n);

433:   if (!user->radiation) {
434:     DMDAGetInfo(user->da2,NULL,&Mda,&Nda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
435:     DMGetCoordinatesLocal(user->da2,&coords);
436:     VecGetArrayRead(coords,&_coords);

438:     if (user->periodic) h = (user->xmax-user->xmin)/Mda;
439:     else h = (user->xmax-user->xmin)/(Mda-1.0);

441:     xmid   = (user->xmax + user->xmin)/2.0;
442:     ymid   = (user->ymax + user->ymin)/2.0;
443:     lambda = 4.0*h;

445:     DMDAGetElements(user->da2,&nele,&nen,&ele);
446:     for (i=0; i < nele; i++) {
447:       idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];

449:       x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
450:       x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
451:       x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];

453:       PetscInt    k;
454:       PetscScalar vals_DDcv[3],vals_cv[3],vals_eta[3],s,hhr,r;
455:       for (k=0; k < 3; k++) {
456:         /*s = PetscAbs(x[k] - xmid);*/
457:         s = PetscSqrtScalar((x[k] - xmid)*(x[k] - xmid) + (y[k] - ymid)*(y[k] - ymid));
458:         if (s <= xwidth*(5.0/64.0)) {
459:           vals_cv[k]   = cv_v;
460:           vals_eta[k]  = eta_v;
461:           vals_DDcv[k] = 0.0;
462:         } else if (s> xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0)) {
463:           /*r = (s - xwidth*(6.0/64.0))/(0.5*lambda);*/
464:           r            = (s - xwidth*(6.0/64.0))/(xwidth/64.0);
465:           hhr          = 0.25*(-r*r*r + 3*r + 2);
466:           vals_cv[k]   = cv_m + (1.0 - hhr)*(cv_v - cv_m);
467:           vals_eta[k]  = eta_m + (1.0 - hhr)*(eta_v - eta_m);
468:           vals_DDcv[k] = (cv_v - cv_m)*r*6.0/(lambda*lambda);
469:         } else {
470:           vals_cv[k]   = cv_m;
471:           vals_eta[k]  = eta_m;
472:           vals_DDcv[k] = 0.0;
473:         }
474:       }

476:       VecSetValuesLocal(user->cv,3,idx,vals_cv,INSERT_VALUES);
477:       VecSetValuesLocal(user->eta,3,idx,vals_eta,INSERT_VALUES);
478:       VecSetValuesLocal(user->work2,3,idx,vals_DDcv,INSERT_VALUES);

480:     }
481:     VecAssemblyBegin(user->cv);
482:     VecAssemblyEnd(user->cv);
483:     VecAssemblyBegin(user->eta);
484:     VecAssemblyEnd(user->eta);
485:     VecAssemblyBegin(user->work2);
486:     VecAssemblyEnd(user->work2);

488:     DMDARestoreElements(user->da2,&nele,&nen,&ele);
489:     VecRestoreArrayRead(coords,&_coords);
490:   } else {
491:     /*VecSet(user->cv,user->initv);*/
492:     /*VecSet(user->ci,user->initv);*/
493:     VecSet(user->cv,.05);
494:     VecSet(user->eta,user->initeta);

496:   }
497:   DPsi(user);
498:   VecCopy(user->DPsiv,user->wv);
499:   VecAXPY(user->wv,-2.0*user->kav,user->work2);

501:   VecGetArray(X,&xx);
502:   VecGetArray(user->wv,&wv_p);
503:   VecGetArray(user->cv,&cv_p);
504:   VecGetArray(user->eta,&eta_p);

506:   for (i=0; i<n/3; i++)
507:   {
508:     xx[3*i]  =wv_p[i];
509:     xx[3*i+1]=cv_p[i];
510:     xx[3*i+2]=eta_p[i];
511:   }

513:   /*
514:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);
515:   VecView(user->wv,view_out);
516:   VecView(user->cv,view_out);
517:   VecView(user->eta,view_out);
518:   PetscViewerDestroy(&view_out);
519:    */

521:   VecRestoreArray(X,&xx);
522:   VecRestoreArray(user->wv,&wv_p);
523:   VecRestoreArray(user->cv,&cv_p);
524:   VecRestoreArray(user->eta,&eta_p);
525:   return(0);
526: }

528: typedef struct {
529:   PetscReal dt,x,y,strength;
530: } RandomValues;


535: PetscErrorCode SetRandomVectors(AppCtx *user,PetscReal t)
536: {
537:   PetscErrorCode      ierr;
538:   static RandomValues *randomvalues = 0;
539:   static PetscInt     randindex     = 0,n; /* indicates how far into the randomvalues we have currently used */
540:   static PetscReal    randtime      = 0; /* indicates time of last radiation event */
541:   PetscInt            i,j,M,N,cnt = 0;
542:   PetscInt            xs,ys,xm,ym;

545:   if (!randomvalues) {
546:     PetscViewer viewer;
547:     char        filename[PETSC_MAX_PATH_LEN];
548:     PetscBool   flg;
549:     PetscInt    seed;

551:     PetscOptionsGetInt(NULL,"-random_seed",&seed,&flg);
552:     if (flg) {
553:       sprintf(filename,"ex61.random.%d",(int)seed);
554:     } else {
555:       PetscStrcpy(filename,"ex61.random");
556:     }
557:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
558:     PetscViewerBinaryRead(viewer,&n,1,PETSC_INT);
559:     PetscMalloc(n*sizeof(RandomValues),&randomvalues);
560:     PetscViewerBinaryRead(viewer,randomvalues,4*n,PETSC_DOUBLE);
561:     for (i=0; i<n; i++) randomvalues[i].dt = randomvalues[i].dt*user->dtevent;
562:     PetscViewerDestroy(&viewer);
563:   }
564:   user->maxevents = PetscMin(user->maxevents,n);

566:   VecSet(user->Pv,0.0);
567:   DMDAGetInfo(user->da1,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
568:   DMDAGetGhostCorners(user->da1,&xs,&ys,0,&xm,&ym,0);
569:   while (user->maxevents > randindex && randtime + randomvalues[randindex].dt < t + user->dt) {  /* radiation event has occured since last time step */
570:     i = ((PetscInt) (randomvalues[randindex].x*M)) - xs;
571:     j = ((PetscInt) (randomvalues[randindex].y*N)) - ys;
572:     if (i >= 0 && i < xm && j >= 0 && j < ym) { /* point is on this process */

574:       /* need to make sure eta at the given point is not great than .8 */
575:       VecSetValueLocal(user->Pv,i  + xm*(j), randomvalues[randindex].strength*user->VG,INSERT_VALUES);
576:     }
577:     randtime += randomvalues[randindex++].dt;
578:     cnt++;
579:   }
580:   PetscPrintf(PETSC_COMM_WORLD,"Number of radiation events %d\n",cnt);
581:   VecAssemblyBegin(user->Pv);
582:   VecAssemblyEnd(user->Pv);

584:   VecCopy(user->Pv,user->Pi);
585:   VecScale(user->Pi,0.9);
586:   VecPointwiseMult(user->Piv,user->Pi,user->Pv);
587:   return(0);
588: }


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

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

605: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
606: {
608:   AppCtx         *user=(AppCtx*)ctx;

611:   *flg = SAME_NONZERO_PATTERN;
612:   MatCopy(user->M,*J,*flg);
613:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
614:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
615:   return(0);
616: }
619: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
620: {
622:   PetscScalar    ***l,***u;
623:   PetscInt       xs,xm,ys,ym;
624:   PetscInt       i,j;

627:   DMDAVecGetArrayDOF(da,xl,&l);
628:   DMDAVecGetArrayDOF(da,xu,&u);

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

632:   for (j=ys; j < ys+ym; j++) {
633:     for (i=xs; i < xs+xm; i++) {
634:       l[j][i][0] = -SNES_VI_INF;
635:       l[j][i][1] = 0.0;
636:       l[j][i][2] = 0.0;
637:       u[j][i][0] = SNES_VI_INF;
638:       u[j][i][1] = 1.0;
639:       u[j][i][2] = 1.0;
640:     }
641:   }

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

650: PetscErrorCode GetParams(AppCtx *user)
651: {
653:   PetscBool      flg;

656:   /* Set default parameters */
657:   user->xmin       = 0.0; user->xmax = 128.0;
658:   user->ymin       = 0.0; user->ymax = 128.0;
659:   user->Mv         = 1.0;
660:   user->L          = 1.0;
661:   user->kaeta      = 1.0;
662:   user->kav        = 0.5;
663:   user->Evf        = 9.09;
664:   user->A          = 9.09;
665:   user->B          = 9.09;
666:   user->cv0        = 1.13e-4;
667:   user->Sv_scalar  = 500.0;
668:   user->dt         = 1.0e-5;
669:   user->T          = 1.0e-2;
670:   user->initv      = .00069;
671:   user->initeta    = .0;
672:   user->graphics   = PETSC_FALSE;
673:   user->periodic   = PETSC_FALSE;
674:   user->lumpedmass = PETSC_FALSE;
675:   user->radiation  = PETSC_FALSE;
676:   /* multidomain modeling */
677:   user->domain    = 0;
678:   user->grain     = 5000.0;
679:   user->Svr       = 0.5;
680:   user->Sir       = 0.5;
681:   user->cv_eq     = 6.9e-4;
682:   user->ci_eq     = 6.9e-4;
683:   user->VG        = 100.0;
684:   user->maxevents = 10;

686:   user->dtevent = user->dt;
687:   PetscOptionsReal("-dtevent","Average time between random events\n","None",user->dtevent,&user->dtevent,&flg);


690:   PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);
691:   PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);
692:   PetscOptionsGetReal(NULL,"-T",&user->T,&flg);
693:   PetscOptionsGetReal(NULL,"-dt",&user->dt,&flg);
694:   PetscOptionsBool("-graphics","Contour plot solutions at each timestep\n","None",user->graphics,&user->graphics,&flg);
695:   PetscOptionsBool("-periodic","Use periodic boundary conditions\n","None",user->periodic,&user->periodic,&flg);
696:   PetscOptionsBool("-radiation","Allow radiation\n","None",user->radiation,&user->radiation,&flg);
697:   PetscOptionsBool("-lumpedmass","Use lumped mass matrix\n","None",user->lumpedmass,&user->lumpedmass,&flg);
698:   PetscOptionsInt("-domain","Number of domains (0=one domain, 1=two domains, 2=multidomain\n","None",user->domain,&user->domain,&flg);
699:   PetscOptionsReal("-VG","Maximum increase in vacancy (or interstitial) concentration due to a cascade event","None",user->VG,&user->VG,&flg);
700:   PetscOptionsReal("-grain","Some bogus factor that controls the strength of the grain boundaries, makes the picture more plausible","None",user->grain,&user->grain,&flg);
701:   PetscOptionsInt("-maxevents","Maximum random events allowed\n","None",user->maxevents,&user->maxevents,&flg);
702:   PetscOptionsInt("-da_refine","da refine \n","None",user->darefine,&user->darefine,&flg);
703:   PetscOptionsInt("-da_grid_x","da grid x\n","None",user->dagrid,&user->dagrid,&flg);
704:   return(0);
705: }


710: PetscErrorCode SetUpMatrices(AppCtx *user)
711: {
713:   PetscInt       nele,nen,i,n;
714:   const PetscInt *ele;
715:   PetscScalar    dt=user->dt,h;

717:   PetscInt    idx[3];
718:   PetscScalar eM_0[3][3],eM_2[3][3];
719:   Mat         M  =user->M;
720:   Mat         M_0=user->M_0;
721:   PetscInt    Mda,Nda;


725:   MatGetLocalSize(M,&n,NULL);
726:   DMDAGetInfo(user->da1,NULL,&Mda,&Nda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

728:   if (Mda!=Nda) printf("Currently different Mda and Nda are not supported");
729:   if (user->periodic) h = (user->xmax-user->xmin)/Mda;
730:   else h = (user->xmax-user->xmin)/(Mda-1.0);

732:   if (user->lumpedmass) {
733:     eM_0[0][0] = eM_0[1][1] = eM_0[2][2] = h*h/6.0;
734:     eM_0[0][1] = eM_0[1][0] = eM_0[0][2] = eM_0[2][0] = eM_0[1][2] = eM_0[2][1] = 0.0;
735:   } else {
736:     eM_0[0][0] = eM_0[1][1] = eM_0[2][2] = h*h/12.0;
737:     eM_0[0][1] = eM_0[0][2] = eM_0[1][0] = eM_0[1][2] = eM_0[2][0] = eM_0[2][1] = h*h/24.0;
738:   }
739:   eM_2[0][0] = 1.0;
740:   eM_2[1][1] = eM_2[2][2] = 0.5;
741:   eM_2[0][1] = eM_2[0][2] = eM_2[1][0]= eM_2[2][0] = -0.5;
742:   eM_2[1][2] = eM_2[2][1] = 0.0;


745:   /* Get local element info */
746:   DMDAGetElements(user->da1,&nele,&nen,&ele);
747:   for (i=0; i < nele; i++) {

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

751:     PetscInt    row,cols[6],r,row_M_0,cols3[3];
752:     PetscScalar vals[6],vals_M_0[3],vals3[3];

754:     for (r=0; r<3; r++) {
755:       row_M_0    = idx[r];
756:       vals_M_0[0]=eM_0[r][0];
757:       vals_M_0[1]=eM_0[r][1];
758:       vals_M_0[2]=eM_0[r][2];

760:       MatSetValuesLocal(M_0,1,&row_M_0,3,idx,vals_M_0,ADD_VALUES);

762:       row     = 3*idx[r];
763:       cols[0] = 3*idx[0];     vals[0] = dt*eM_2[r][0]*user->Mv;
764:       cols[1] = 3*idx[1];     vals[1] = dt*eM_2[r][1]*user->Mv;
765:       cols[2] = 3*idx[2];     vals[2] = dt*eM_2[r][2]*user->Mv;
766:       cols[3] = 3*idx[0]+1;   vals[3] = eM_0[r][0];
767:       cols[4] = 3*idx[1]+1;   vals[4] = eM_0[r][1];
768:       cols[5] = 3*idx[2]+1;   vals[5] = eM_0[r][2];

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

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

781:       /* Insert values in matrix M for 2nd dof */
782:       MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);

784:       row      = 3*idx[r]+2;
785:       cols3[0] = 3*idx[0]+2;   vals3[0] = eM_0[r][0] + user->dt*2.0*user->L*user->kaeta*eM_2[r][0];
786:       cols3[1] = 3*idx[1]+2;   vals3[1] = eM_0[r][1] + user->dt*2.0*user->L*user->kaeta*eM_2[r][1];
787:       cols3[2] = 3*idx[2]+2;   vals3[2] = eM_0[r][2] + user->dt*2.0*user->L*user->kaeta*eM_2[r][2];

789:       MatSetValuesLocal(M,1,&row,3,cols3,vals3,ADD_VALUES);
790:     }
791:   }

793:   MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
794:   MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);

796:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
797:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

799:   DMDARestoreElements(user->da1,&nele,&nen,&ele);
800:   return(0);
801: }

805: PetscErrorCode CheckRedundancy(SNES snes, IS act, IS *outact, DM da)
806: {
808:   PetscScalar    **uin,**uout;
809:   Vec            UIN, UOUT;
810:   PetscInt       xs,xm,*outindex;
811:   const PetscInt *index;
812:   PetscInt       k,i,l,n,M,cnt=0;

815:   DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
816:   DMGetGlobalVector(da,&UIN);
817:   VecSet(UIN,0.0);
818:   DMGetLocalVector(da,&UOUT);
819:   DMDAVecGetArrayDOF(da,UIN,&uin);
820:   ISGetIndices(act,&index);
821:   ISGetLocalSize(act,&n);
822:   for (k=0; k<n; k++) {
823:     l = index[k]%5;
824:     i = index[k]/5;

826:     uin[i][l]=1.0;
827:   }
828:   printf("Number of active constraints before applying redundancy %d\n",n);
829:   ISRestoreIndices(act,&index);
830:   DMDAVecRestoreArrayDOF(da,UIN,&uin);
831:   DMGlobalToLocalBegin(da,UIN,INSERT_VALUES,UOUT);
832:   DMGlobalToLocalEnd(da,UIN,INSERT_VALUES,UOUT);
833:   DMDAVecGetArrayDOF(da,UOUT,&uout);

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

837:   for (i=xs; i < xs+xm;i++) {
838:     if (uout[i-1][1] && uout[i][1] && uout[i+1][1]) uout[i][0] = 1.0;
839:     if (uout[i-1][3] && uout[i][3] && uout[i+1][3]) uout[i][2] = 1.0;
840:   }

842:   for (i=xs; i < xs+xm;i++) {
843:     for (l=0;l<5;l++) {
844:       if (uout[i][l]) cnt++;
845:     }
846:   }

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


851:   PetscMalloc(cnt*sizeof(PetscInt),&outindex);
852:   cnt  = 0;

854:   for (i=xs; i < xs+xm;i++) {
855:     for (l=0;l<5;l++) {
856:       if (uout[i][l]) outindex[cnt++] = 5*(i)+l;
857:     }
858:   }


861:   ISCreateGeneral(PETSC_COMM_WORLD,cnt,outindex,PETSC_OWN_POINTER,outact);
862:   DMDAVecRestoreArrayDOF(da,UOUT,&uout);
863:   DMRestoreGlobalVector(da,&UIN);
864:   DMRestoreLocalVector(da,&UOUT);
865:   return(0);
866: }

870: PetscErrorCode Phi(AppCtx *user)
871: {
872:   PetscErrorCode    ierr;
873:   PetscScalar       xmid, xqu, lambda, h,x[3],y[3];
874:   Vec               coords;
875:   const PetscScalar *_coords;
876:   PetscInt          nele,nen,i,idx[3],Mda,Nda;
877:   const PetscInt    *ele;
878:   /*PetscViewer        view;*/

881:   DMDAGetInfo(user->da1,NULL,&Mda,&Nda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
882:   DMGetCoordinatesLocal(user->da2,&coords);
883:   VecGetArrayRead(coords,&_coords);

885:   h      = (user->xmax - user->xmin)/Mda;
886:   xmid   = (user->xmin + user->xmax)/2.0;
887:   xqu    = (user->xmin + xmid)/2.0;
888:   lambda = 4.0*h;


891:   DMDAGetElements(user->da2,&nele,&nen,&ele);
892:   for (i=0; i < nele; i++) {
893:     idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
894:     /*printf("idx[0]=%d,idx[1]=%d,idx[2]=%d\n",idx[0],idx[1],idx[2]);*/

896:     x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
897:     x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
898:     x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];

900:     /*printf("x[0]=%f,x[1]=%f,x[2]=%f\n",x[0],x[1],x[2]);*/
901:     /*printf("y[0]=%f,y[1]=%f,y[2]=%f\n",y[0],y[1],y[2]);*/

903:     PetscScalar vals1[3],vals2[3],dist1,dist2,s1,r,hhr,xc1,xc2;
904:     PetscInt    k;

906:     xc1 = user->xmin;
907:     xc2 = xmid;

909:     /*VecSet(user->phi1,0.0);*/
910:     for (k=0;k < 3; k++) {
911:       if (x[k]-xqu > 0) s1 = (x[k] - xqu);
912:       else s1 = -(x[k] - xqu);

914:       if (x[k] - xc1 > 0) dist1 = (x[k] - xc1);
915:       else dist1 = -(x[k] - xc1);

917:       if (x[k] - xc2 > 0) dist2 = (x[k] - xc2);
918:       else dist2 = -(x[k] - xc2);

920:       if (dist1 <= 0.5*lambda) {
921:         r = (x[k]-xc1)/(0.5*lambda);
922:         hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
923:         vals1[k] = hhr;
924:       } else if (dist2 <= 0.5*lambda) {
925:         r = (x[k]-xc2)/(0.5*lambda);
926:         hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
927:         vals1[k] = 1.0 - hhr;
928:       } else if (s1 <= xqu - 2.0*h) {
929:         vals1[k] = 1.0;
930:       } else if ((user->xmax-h)-x[k] < 0.1*h) {
931:       /*else if (abs(x[k]-(user->xmax-h)) < 0.1*h) {*/
932:         vals1[k] = .15625;
933:       } else {
934:         vals1[k] = 0.0;
935:       }
936:     }

938:     VecSetValuesLocal(user->phi1,3,idx,vals1,INSERT_VALUES);

940:     xc1 = xmid;
941:     xc2 = user->xmax;

943:     /*VecSet(user->phi2,0.0);*/
944:     for (k=0;k < 3; k++) {
945:       /*
946:       s1 = abs(x[k] - (xqu+xmid));
947:       dist1 = abs(x[k] - xc1);
948:       dist2 = abs(x[k] - xc2);
949:        */
950:       if (x[k]-(xqu + xmid) > 0) s1 = (x[k] - (xqu + xmid));
951:       else s1 = -(x[k] - (xqu + xmid));

953:       if (x[k] - xc1 > 0) dist1 = (x[k] - xc1);
954:       else dist1 = -(x[k] - xc1);

956:       if (x[k] - xc2 > 0) dist2 = (x[k] - xc2);
957:       else dist2 = -(x[k] - xc2);

959:       if (dist1 <= 0.5*lambda) {
960:         r = (x[k]-xc1)/(0.5*lambda);
961:         hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
962:         vals2[k] = hhr;
963:       } else if (dist2 <= 0.5*lambda) {
964:         r = -(x[k]-xc2)/(0.5*lambda);
965:         hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
966:         vals2[k] = hhr;
967:       } else if (s1 <= xqu - 2.0*h) {
968:         vals2[k] = 1.0;
969:       } else if (x[k]-(user->xmin) < 0.1*h) {
970:         vals2[k] = 0.5;
971:       } else if ((x[k]-(user->xmin+h)) < 0.1*h) {
972:         vals2[k] = .15625;
973:       } else {
974:         vals2[k] = 0.0;
975:       }

977:     }

979:     VecSetValuesLocal(user->phi2,3,idx,vals2,INSERT_VALUES);
980:     /*
981:     for (k=0;k < 3; k++) {
982:       vals_sum[k] = vals1[k]*vals1[k] + vals2[k]*vals2[k];
983:     }
984:      */
985:     /*VecSetValuesLocal(user->Phi2D_V,3,idx,vals_sum,INSERT_VALUES);*/

987:   }

989:   VecAssemblyBegin(user->phi1);
990:   VecAssemblyEnd(user->phi1);
991:   VecAssemblyBegin(user->phi2);
992:   VecAssemblyEnd(user->phi2);

994:   /*
995:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_phi",FILE_MODE_WRITE,&view);
996:   VecView(user->phi1,view);
997:   VecView(user->phi2,view);
998:    */

1000:   /*VecView(user->phi1,0);*/
1001:   /*VecView(user->phi2,0);*/

1003:   VecPointwiseMult(user->phi1,user->phi1,user->phi1);
1004:   VecPointwiseMult(user->phi2,user->phi2,user->phi2);
1005:   /*
1006:   VecView(user->phi1,view);
1007:   VecView(user->phi2,view);
1008:    */

1010:   VecCopy(user->phi1,user->Phi2D_V);
1011:   VecAXPY(user->Phi2D_V,1.0,user->phi2);
1012:   /*VecView(user->Phi2D_V,0);*/

1014:   /*VecView(user->Phi2D_V,view);*/
1015:   /*PetscViewerDestroy(&view);*/
1016:   return(0);
1017: }

1021: PetscErrorCode Phi_read(AppCtx *user)
1022: {
1024:   PetscReal      *values;
1025:   PetscViewer    viewer;
1026:   PetscInt       power;

1029:   power = user->darefine + (PetscInt)(PetscLogScalar(user->dagrid)/PetscLogScalar(2.0));
1030:   switch (power) {
1031:   case 6:
1032:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi4",FILE_MODE_READ,&viewer);
1033:     VecLoad(user->Phi2D_V,viewer);
1034:     PetscViewerDestroy(&viewer);
1035:     break;
1036:   case 7:
1037:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi3",FILE_MODE_READ,&viewer);
1038:     VecLoad(user->Phi2D_V,viewer);
1039:     PetscViewerDestroy(&viewer);
1040:     break;
1041:   case 8:
1042:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi2",FILE_MODE_READ,&viewer);
1043:     VecLoad(user->Phi2D_V,viewer);
1044:     PetscViewerDestroy(&viewer);
1045:     break;
1046:   case 9:
1047:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi1",FILE_MODE_READ,&viewer);
1048:     VecLoad(user->Phi2D_V,viewer);
1049:     PetscViewerDestroy(&viewer);
1050:     break;
1051:   case 10:
1052:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi",FILE_MODE_READ,&viewer);
1053:     VecLoad(user->Phi2D_V,viewer);
1054:     PetscViewerDestroy(&viewer);
1055:     break;
1056:   case 11:
1057:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim1",FILE_MODE_READ,&viewer);
1058:     VecLoad(user->Phi2D_V,viewer);
1059:     PetscViewerDestroy(&viewer);
1060:     break;
1061:   case 12:
1062:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim2",FILE_MODE_READ,&viewer);
1063:     VecLoad(user->Phi2D_V,viewer);
1064:     PetscViewerDestroy(&viewer);
1065:     break;
1066:   case 13:
1067:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim3",FILE_MODE_READ,&viewer);
1068:     VecLoad(user->Phi2D_V,viewer);
1069:     PetscViewerDestroy(&viewer);
1070:     break;
1071:   case 14:
1072:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim4",FILE_MODE_READ,&viewer);
1073:     VecLoad(user->Phi2D_V,viewer);
1074:     PetscViewerDestroy(&viewer);
1075:     break;
1076:   }
1077:   return(0);
1078: }