Actual source code: ex65.c
petsc-3.3-p7 2013-05-11
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);
76:
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,PETSC_NULL,PETSC_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,PETSC_NULL,PETSC_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,PETSC_NULL,PETSC_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,PETSC_NULL,PETSC_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,PETSC_NULL,PETSC_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,PETSC_NULL,PETSC_NULL,&user.da2);
88: }
89: /* Set Element type (triangular) */
90: DMDASetElementType(user.da1,DMDA_ELEMENT_P1);
91: DMDASetElementType(user.da2,DMDA_ELEMENT_P1);
92:
93: /* Set x and y coordinates */
94: DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,user.ymin,user.ymax,PETSC_NULL,PETSC_NULL);
95: DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,user.ymin,user.ymax,PETSC_NULL,PETSC_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);
102:
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);
133:
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); */
152:
153:
154: if (user.graphics) {
155: //PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds);
157: VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
158: }
159:
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,PETSC_NULL,x);
196: //VecView(x,view_out);
198:
199: if (user.graphics) {
200: VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
201: }
202:
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;
221:
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;
274:
276: VecGetLocalSize(user->wv,&n);
277:
278: VecGetArray(X,&xx);
279: VecGetArray(user->wv,&wv_p);
280: VecGetArray(user->cv,&cv_p);
281: VecGetArray(user->eta,&eta_p);
282:
283:
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:
294: return(0);
295: }
299: PetscErrorCode Update_q(AppCtx *user)
300: {
302: PetscScalar *q_p, *w1, *w2;
303: PetscInt n,i;
306:
307: VecGetArray(user->q,&q_p);
308: VecGetArray(user->work1,&w1);
309: VecGetArray(user->work2,&w2);
310: VecGetLocalSize(user->wv,&n);
313: VecSet(user->work1,0.0);
314: if (user->radiation) {
315: VecAXPY(user->work1,-1.0,user->Pv);
316: }
317: if (user->domain) {
318: VecCopy(user->cv,user->work3);
319: VecShift(user->work3,-1.0*user->cv_eq);
320: VecCopy(user->Phi2D_V,user->work4);
321: VecScale(user->work4,-1.0);
322: VecShift(user->work4,1.0);
323: VecPointwiseMult(user->work4,user->work4,user->work3);
324: VecScale(user->work4,user->Svr);
325: // Parameter tuning: user->grain
326: // 5000.0 worked for refinement 2 and 5
327: //10000.0 worked for refinement 2, 3, 4, 5
328: VecAXPY(user->work1,user->grain,user->work4);
329: }
330: VecScale(user->work1,user->dt);
331: VecAXPY(user->work1,-1.0,user->cv);
332: MatMult(user->M_0,user->work1,user->work2);
333: for (i=0;i<n;i++) {
334: q_p[3*i]=w2[i];
335: }
337: MatMult(user->M_0,user->DPsiv,user->work1);
338: for (i=0;i<n;i++) {
339: q_p[3*i+1]=w1[i];
340: }
342: VecCopy(user->DPsieta,user->work1);
343: VecScale(user->work1,user->L*user->dt);
344: VecAXPY(user->work1,-1.0,user->eta);
345: if (user->radiation) {
346: VecAXPY(user->work1,-1.0,user->Piv);
347: }
348: MatMult(user->M_0,user->work1,user->work2);
349: for (i=0;i<n;i++) {
350: q_p[3*i+2]=w2[i];
351: }
353: VecRestoreArray(user->q,&q_p);
354: VecRestoreArray(user->work1,&w1);
355: VecRestoreArray(user->work2,&w2);
356:
357: return(0);
358: }
362: PetscErrorCode DPsi(AppCtx* user)
363: {
364: PetscErrorCode ierr;
365: PetscScalar Evf=user->Evf,A=user->A,B=user->B,cv0=user->cv0;
366: PetscScalar *cv_p,*eta_p,*logcv_p,*logcv2_p,*DPsiv_p,*DPsieta_p;
367: PetscInt n,i;
371: VecGetLocalSize(user->cv,&n);
372: VecGetArray(user->cv,&cv_p);
373: VecGetArray(user->eta,&eta_p);
374: VecGetArray(user->logcv,&logcv_p);
375: VecGetArray(user->logcv2,&logcv2_p);
376: VecGetArray(user->DPsiv,&DPsiv_p);
377: VecGetArray(user->DPsieta,&DPsieta_p);
379: Llog(user->cv,user->logcv);
381: VecCopy(user->cv,user->work1);
382: VecScale(user->work1,-1.0);
383: VecShift(user->work1,1.0);
384: Llog(user->work1,user->logcv2);
386: for (i=0;i<n;i++)
387: {
388: 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];
389:
390: 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];
391:
392: }
394: VecRestoreArray(user->cv,&cv_p);
395: VecRestoreArray(user->eta,&eta_p);
396: VecGetArray(user->logcv,&logcv_p);
397: VecGetArray(user->logcv2,&logcv2_p);
398: VecRestoreArray(user->DPsiv,&DPsiv_p);
399: VecRestoreArray(user->DPsieta,&DPsieta_p);
402: return(0);
404: }
409: PetscErrorCode Llog(Vec X, Vec Y)
410: {
411: PetscErrorCode ierr;
412: PetscScalar *x,*y;
413: PetscInt n,i;
416:
417: VecGetArray(X,&x);
418: VecGetArray(Y,&y);
419: VecGetLocalSize(X,&n);
420: for (i=0;i<n;i++) {
421: if (x[i] < 1.0e-12) {
422: y[i] = log(1.0e-12);
423: }
424: else {
425: y[i] = log(x[i]);
426: }
427: }
428:
429: return(0);
430: }
434: PetscErrorCode SetInitialGuess(Vec X,AppCtx* user)
435: {
436: PetscErrorCode ierr;
437: PetscInt n,i,Mda,Nda;
438: PetscScalar *xx,*cv_p,*wv_p,*eta_p;
439: //PetscViewer view_out;
440: /* needed for the void growth case */
441: PetscScalar xmid,ymid,cv_v=1.0,cv_m=user->Sv_scalar*user->cv0,eta_v=1.0,eta_m=0.0,h,lambda;
442: PetscInt nele,nen,idx[3];
443: const PetscInt *ele;
444: PetscScalar x[3],y[3];
445: Vec coords;
446: const PetscScalar *_coords;
447: PetscScalar xwidth = user->xmax - user->xmin, ywidth = user->ymax - user->ymin;
451: VecGetLocalSize(X,&n);
452:
453: if (!user->radiation) {
454: DMDAGetInfo(user->da2,PETSC_NULL,&Mda,&Nda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
455: DMDAGetGhostedCoordinates(user->da2,&coords);
456: VecGetArrayRead(coords,&_coords);
458: if (user->periodic) {
459: h = (user->xmax-user->xmin)/Mda;
460: } else {
461: h = (user->xmax-user->xmin)/(Mda-1.0);
462: }
463:
464: xmid = (user->xmax + user->xmin)/2.0;
465: ymid = (user->ymax + user->ymin)/2.0;
466: lambda = 4.0*h;
467:
468: DMDAGetElements(user->da2,&nele,&nen,&ele);
469: for (i=0;i < nele; i++) {
470: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
471:
472: x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
473: x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
474: x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
475:
476: PetscInt k;
477: PetscScalar vals_DDcv[3],vals_cv[3],vals_eta[3],s,hhr,r;
478: for (k=0; k < 3 ; k++) {
479: //s = PetscAbs(x[k] - xmid);
480: s = sqrt((x[k] - xmid)*(x[k] - xmid) + (y[k] - ymid)*(y[k] - ymid));
481: if (s <= xwidth*(5.0/64.0)) {
482: vals_cv[k] = cv_v;
483: vals_eta[k] = eta_v;
484: vals_DDcv[k] = 0.0;
485: } else if (s> xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0) ) {
486: //r = (s - xwidth*(6.0/64.0) )/(0.5*lambda);
487: r = (s - xwidth*(6.0/64.0) )/(xwidth/64.0);
488: hhr = 0.25*(-r*r*r + 3*r + 2);
489: vals_cv[k] = cv_m + (1.0 - hhr)*(cv_v - cv_m);
490: vals_eta[k] = eta_m + (1.0 - hhr)*(eta_v - eta_m);
491: vals_DDcv[k] = (cv_v - cv_m)*r*6.0/(lambda*lambda);
492: } else {
493: vals_cv[k] = cv_m;
494: vals_eta[k] = eta_m;
495: vals_DDcv[k] = 0.0;
496: }
497: }
498:
499: VecSetValuesLocal(user->cv,3,idx,vals_cv,INSERT_VALUES);
500: VecSetValuesLocal(user->eta,3,idx,vals_eta,INSERT_VALUES);
501: VecSetValuesLocal(user->work2,3,idx,vals_DDcv,INSERT_VALUES);
502:
503: }
504: VecAssemblyBegin(user->cv);
505: VecAssemblyEnd(user->cv);
506: VecAssemblyBegin(user->eta);
507: VecAssemblyEnd(user->eta);
508: VecAssemblyBegin(user->work2);
509: VecAssemblyEnd(user->work2);
511: DMDARestoreElements(user->da2,&nele,&nen,&ele);
512: VecRestoreArrayRead(coords,&_coords);
513: } else {
514: //VecSet(user->cv,user->initv);
515: //VecSet(user->ci,user->initv);
516: VecSet(user->cv,.05);
517: VecSet(user->eta,user->initeta);
519: }
520: DPsi(user);
521: VecCopy(user->DPsiv,user->wv);
522: VecAXPY(user->wv,-2.0*user->kav,user->work2);
524: VecGetArray(X,&xx);
525: VecGetArray(user->wv,&wv_p);
526: VecGetArray(user->cv,&cv_p);
527: VecGetArray(user->eta,&eta_p);
529: for (i=0;i<n/3;i++)
530: {
531: xx[3*i]=wv_p[i];
532: xx[3*i+1]=cv_p[i];
533: xx[3*i+2]=eta_p[i];
534: }
536: /*
537: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);
538: VecView(user->wv,view_out);
539: VecView(user->cv,view_out);
540: VecView(user->eta,view_out);
541: PetscViewerDestroy(&view_out);
542: */
544: VecRestoreArray(X,&xx);
545: VecRestoreArray(user->wv,&wv_p);
546: VecRestoreArray(user->cv,&cv_p);
547: VecRestoreArray(user->eta,&eta_p);
548: return(0);
549: }
551: typedef struct {
552: PetscReal dt,x,y,strength;
553: } RandomValues;
558: PetscErrorCode SetRandomVectors(AppCtx* user,PetscReal t)
559: {
560: PetscErrorCode ierr;
561: static RandomValues *randomvalues = 0;
562: static PetscInt randindex = 0,n; /* indicates how far into the randomvalues we have currently used */
563: static PetscReal randtime = 0; /* indicates time of last radiation event */
564: PetscInt i,j,M,N,cnt = 0;
565: PetscInt xs,ys,xm,ym;
568: if (!randomvalues) {
569: PetscViewer viewer;
570: char filename[PETSC_MAX_PATH_LEN];
571: PetscBool flg;
572: PetscInt seed;
574: PetscOptionsGetInt(PETSC_NULL,"-random_seed",&seed,&flg);
575: if (flg) {
576: sprintf(filename,"ex61.random.%d",(int)seed);
577: } else {
578: PetscStrcpy(filename,"ex61.random");
579: }
580: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
581: PetscViewerBinaryRead(viewer,&n,1,PETSC_INT);
582: PetscMalloc(n*sizeof(RandomValues),&randomvalues);
583: PetscViewerBinaryRead(viewer,randomvalues,4*n,PETSC_DOUBLE);
584: for (i=0; i<n; i++) randomvalues[i].dt = randomvalues[i].dt*user->dtevent;
585: PetscViewerDestroy(&viewer);
586: }
587: user->maxevents = PetscMin(user->maxevents,n);
589: VecSet(user->Pv,0.0);
590: DMDAGetInfo(user->da1,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
591: DMDAGetGhostCorners(user->da1,&xs,&ys,0,&xm,&ym,0);
592: while (user->maxevents > randindex && randtime + randomvalues[randindex].dt < t + user->dt) { /* radiation event has occured since last time step */
593: i = ((PetscInt) (randomvalues[randindex].x*M)) - xs;
594: j = ((PetscInt) (randomvalues[randindex].y*N)) - ys;
595: if (i >= 0 && i < xm && j >= 0 && j < ym) { /* point is on this process */
597: /* need to make sure eta at the given point is not great than .8 */
598: VecSetValueLocal(user->Pv,i + xm*(j), randomvalues[randindex].strength*user->VG,INSERT_VALUES);
599: }
600: randtime += randomvalues[randindex++].dt;
601: cnt++;
602: }
603: PetscPrintf(PETSC_COMM_WORLD,"Number of radiation events %d\n",cnt);
604: VecAssemblyBegin(user->Pv);
605: VecAssemblyEnd(user->Pv);
607: VecCopy(user->Pv,user->Pi);
608: VecScale(user->Pi,0.9);
609: VecPointwiseMult(user->Piv,user->Pi,user->Pv);
610: return(0);
611:
612: }
617: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ctx)
618: {
620: AppCtx *user=(AppCtx*)ctx;
621:
623: MatMultAdd(user->M,X,user->q,F);
624: return(0);
625: }
629: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
630: {
632: AppCtx *user=(AppCtx*)ctx;
633:
635: *flg = SAME_NONZERO_PATTERN;
636: MatCopy(user->M,*J,*flg);
637: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
638: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
639: return(0);
640: }
643: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
644: {
646: PetscScalar ***l,***u;
647: PetscInt xs,xm,ys,ym;
648: PetscInt i,j;
649:
651: DMDAVecGetArrayDOF(da,xl,&l);
652: DMDAVecGetArrayDOF(da,xu,&u);
653:
654: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
655:
656: for (j=ys; j < ys+ym; j++) {
657: for(i=xs; i < xs+xm;i++) {
658: l[j][i][0] = -SNES_VI_INF;
659: l[j][i][1] = 0.0;
660: l[j][i][2] = 0.0;
661: u[j][i][0] = SNES_VI_INF;
662: u[j][i][1] = 1.0;
663: u[j][i][2] = 1.0;
664: }
665: }
666:
667: DMDAVecRestoreArrayDOF(da,xl,&l);
668: DMDAVecRestoreArrayDOF(da,xu,&u);
669: return(0);
670: }
674: PetscErrorCode GetParams(AppCtx* user)
675: {
677: PetscBool flg;
678:
680:
681: /* Set default parameters */
682: user->xmin = 0.0; user->xmax = 128.0;
683: user->ymin = 0.0; user->ymax = 128.0;
684: user->Mv = 1.0;
685: user->L = 1.0;
686: user->kaeta = 1.0;
687: user->kav = 0.5;
688: user->Evf = 9.09;
689: user->A = 9.09;
690: user->B = 9.09;
691: user->cv0 = 1.13e-4;
692: user->Sv_scalar = 500.0;
693: user->dt = 1.0e-5;
694: user->T = 1.0e-2;
695: user->initv = .00069;
696: user->initeta = .0;
697: user->graphics = PETSC_FALSE;
698: user->periodic = PETSC_FALSE;
699: user->lumpedmass = PETSC_FALSE;
700: user->radiation = PETSC_FALSE;
701: /* multidomain modeling */
702: user->domain = 0;
703: user->grain = 5000.0;
704: user->Svr = 0.5;
705: user->Sir = 0.5;
706: user->cv_eq = 6.9e-4;
707: user->ci_eq = 6.9e-4;
708: user->VG = 100.0;
709: user->maxevents = 10;
710:
711: user->dtevent = user->dt;
712: PetscOptionsReal("-dtevent","Average time between random events\n","None",user->dtevent,&user->dtevent,&flg);
713:
714:
715: PetscOptionsGetReal(PETSC_NULL,"-xmin",&user->xmin,&flg);
716: PetscOptionsGetReal(PETSC_NULL,"-xmax",&user->xmax,&flg);
717: PetscOptionsGetReal(PETSC_NULL,"-T",&user->T,&flg);
718: PetscOptionsGetReal(PETSC_NULL,"-dt",&user->dt,&flg);
719: PetscOptionsBool("-graphics","Contour plot solutions at each timestep\n","None",user->graphics,&user->graphics,&flg);
720: PetscOptionsBool("-periodic","Use periodic boundary conditions\n","None",user->periodic,&user->periodic,&flg);
721: PetscOptionsBool("-radiation","Allow radiation\n","None",user->radiation,&user->radiation,&flg);
722: PetscOptionsBool("-lumpedmass","Use lumped mass matrix\n","None",user->lumpedmass,&user->lumpedmass,&flg);
723: PetscOptionsInt("-domain","Number of domains (0=one domain, 1=two domains, 2=multidomain\n","None",user->domain,&user->domain,&flg);
724: PetscOptionsReal("-VG","Maximum increase in vacancy (or interstitial) concentration due to a cascade event","None",user->VG,&user->VG,&flg);
725: PetscOptionsReal("-grain","Some bogus factor that controls the strength of the grain boundaries, makes the picture more plausible","None",user->grain,&user->grain,&flg);
726: PetscOptionsInt("-maxevents","Maximum random events allowed\n","None",user->maxevents,&user->maxevents,&flg);
727: PetscOptionsInt("-da_refine","da refine \n","None",user->darefine,&user->darefine,&flg);
728: PetscOptionsInt("-da_grid_x","da grid x\n","None",user->dagrid,&user->dagrid,&flg);
730: return(0);
731: }
736: PetscErrorCode SetUpMatrices(AppCtx* user)
737: {
738: PetscErrorCode ierr;
739: PetscInt nele,nen,i,n;
740: const PetscInt *ele;
741: PetscScalar dt=user->dt,h;
742:
743: PetscInt idx[3];
744: PetscScalar eM_0[3][3],eM_2[3][3];
745: Mat M=user->M;
746: Mat M_0=user->M_0;
747: PetscInt Mda,Nda;
749:
753: MatGetLocalSize(M,&n,PETSC_NULL);
754: DMDAGetInfo(user->da1,PETSC_NULL,&Mda,&Nda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
756: if (Mda!=Nda) {
757: printf("Currently different Mda and Nda are not supported");
758: }
759: if (user->periodic) {
760: h = (user->xmax-user->xmin)/Mda;
761: } else {
762: h = (user->xmax-user->xmin)/(Mda-1.0);
763: }
764: if (user->lumpedmass) {
765: eM_0[0][0] = eM_0[1][1] = eM_0[2][2] = h*h/6.0;
766: 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;
767: } else {
768: eM_0[0][0] = eM_0[1][1] = eM_0[2][2] = h*h/12.0;
769: 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;
770: }
771: eM_2[0][0] = 1.0;
772: eM_2[1][1] = eM_2[2][2] = 0.5;
773: eM_2[0][1] = eM_2[0][2] = eM_2[1][0]= eM_2[2][0] = -0.5;
774: eM_2[1][2] = eM_2[2][1] = 0.0;
777: /* Get local element info */
778: DMDAGetElements(user->da1,&nele,&nen,&ele);
779: for(i=0;i < nele;i++) {
780:
781: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
783: PetscInt row,cols[6],r,row_M_0,cols3[3];
784: PetscScalar vals[6],vals_M_0[3],vals3[3];
786: for(r=0;r<3;r++) {
787: row_M_0 = idx[r];
788: vals_M_0[0]=eM_0[r][0];
789: vals_M_0[1]=eM_0[r][1];
790: vals_M_0[2]=eM_0[r][2];
792: MatSetValuesLocal(M_0,1,&row_M_0,3,idx,vals_M_0,ADD_VALUES);
794: row = 3*idx[r];
795: cols[0] = 3*idx[0]; vals[0] = dt*eM_2[r][0]*user->Mv;
796: cols[1] = 3*idx[1]; vals[1] = dt*eM_2[r][1]*user->Mv;
797: cols[2] = 3*idx[2]; vals[2] = dt*eM_2[r][2]*user->Mv;
798: cols[3] = 3*idx[0]+1; vals[3] = eM_0[r][0];
799: cols[4] = 3*idx[1]+1; vals[4] = eM_0[r][1];
800: cols[5] = 3*idx[2]+1; vals[5] = eM_0[r][2];
801:
802: /* Insert values in matrix M for 1st dof */
803: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
804:
805: row = 3*idx[r]+1;
806: cols[0] = 3*idx[0]; vals[0] = -eM_0[r][0];
807: cols[1] = 3*idx[1]; vals[1] = -eM_0[r][1];
808: cols[2] = 3*idx[2]; vals[2] = -eM_0[r][2];
809: cols[3] = 3*idx[0]+1; vals[3] = 2.0*user->kav*eM_2[r][0];
810: cols[4] = 3*idx[1]+1; vals[4] = 2.0*user->kav*eM_2[r][1];
811: cols[5] = 3*idx[2]+1; vals[5] = 2.0*user->kav*eM_2[r][2];
813: /* Insert values in matrix M for 2nd dof */
814: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
815:
816: row = 3*idx[r]+2;
817: cols3[0] = 3*idx[0]+2; vals3[0] = eM_0[r][0] + user->dt*2.0*user->L*user->kaeta*eM_2[r][0];
818: cols3[1] = 3*idx[1]+2; vals3[1] = eM_0[r][1] + user->dt*2.0*user->L*user->kaeta*eM_2[r][1];
819: cols3[2] = 3*idx[2]+2; vals3[2] = eM_0[r][2] + user->dt*2.0*user->L*user->kaeta*eM_2[r][2];
820: MatSetValuesLocal(M,1,&row,3,cols3,vals3,ADD_VALUES);
821: }
822: }
824: MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
825: MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);
826:
827: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
828: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
829:
830: DMDARestoreElements(user->da1,&nele,&nen,&ele);
831:
832:
833: return(0);
834: }
838: PetscErrorCode CheckRedundancy(SNES snes, IS act, IS *outact, DM da)
839: {
841: PetscScalar **uin,**uout;
842: Vec UIN, UOUT;
843: PetscInt xs,xm,*outindex;
844: const PetscInt *index;
845: PetscInt k,i,l,n,M,cnt=0;
846:
848: DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
849: DMGetGlobalVector(da,&UIN);
850: VecSet(UIN,0.0);
851: DMGetLocalVector(da,&UOUT);
852: DMDAVecGetArrayDOF(da,UIN,&uin);
853: ISGetIndices(act,&index);
854: ISGetLocalSize(act,&n);
855: for (k=0;k<n;k++){
856: l = index[k]%5;
857: i = index[k]/5;
858: uin[i][l]=1.0;
859: }
860: printf("Number of active constraints before applying redundancy %d\n",n);
861: ISRestoreIndices(act,&index);
862: DMDAVecRestoreArrayDOF(da,UIN,&uin);
863: DMGlobalToLocalBegin(da,UIN,INSERT_VALUES,UOUT);
864: DMGlobalToLocalEnd(da,UIN,INSERT_VALUES,UOUT);
865: DMDAVecGetArrayDOF(da,UOUT,&uout);
866:
867: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
869: for(i=xs; i < xs+xm;i++) {
870: if (uout[i-1][1] && uout[i][1] && uout[i+1][1])
871: uout[i][0] = 1.0;
872: if (uout[i-1][3] && uout[i][3] && uout[i+1][3])
873: uout[i][2] = 1.0;
874: }
876: for(i=xs; i < xs+xm;i++) {
877: for(l=0;l<5;l++) {
878: if (uout[i][l])
879: cnt++;
880: }
881: }
883: printf("Number of active constraints after applying redundancy %d\n",cnt);
884:
886: PetscMalloc(cnt*sizeof(PetscInt),&outindex);
887: cnt = 0;
888:
889: for(i=xs; i < xs+xm;i++) {
890: for(l=0;l<5;l++) {
891: if (uout[i][l])
892: outindex[cnt++] = 5*(i)+l;
893: }
894: }
895:
896:
897: ISCreateGeneral(PETSC_COMM_WORLD,cnt,outindex,PETSC_OWN_POINTER,outact);
898: DMDAVecRestoreArrayDOF(da,UOUT,&uout);
899: DMRestoreGlobalVector(da,&UIN);
900: DMRestoreLocalVector(da,&UOUT);
901: return(0);
902: }
906: PetscErrorCode Phi(AppCtx* user)
907: {
908: PetscErrorCode ierr;
909: PetscScalar xmid, xqu, lambda, h,x[3],y[3];
910: Vec coords;
911: const PetscScalar *_coords;
912: PetscInt nele,nen,i,idx[3],Mda,Nda;
913: const PetscInt *ele;
914: //PetscViewer view;
918: DMDAGetInfo(user->da1,PETSC_NULL,&Mda,&Nda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
919: DMDAGetGhostedCoordinates(user->da2,&coords);
920: VecGetArrayRead(coords,&_coords);
922: h = (user->xmax - user->xmin)/Mda;
923: xmid = (user->xmin + user->xmax)/2.0;
924: xqu = (user->xmin + xmid)/2.0;
925: lambda = 4.0*h;
928: DMDAGetElements(user->da2,&nele,&nen,&ele);
929: for (i=0;i < nele; i++) {
930: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
931: //printf("idx[0]=%d,idx[1]=%d,idx[2]=%d\n",idx[0],idx[1],idx[2]);
933: x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
934: x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
935: x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
937: //printf("x[0]=%f,x[1]=%f,x[2]=%f\n",x[0],x[1],x[2]);
938: //printf("y[0]=%f,y[1]=%f,y[2]=%f\n",y[0],y[1],y[2]);
939:
940: PetscScalar vals1[3],vals2[3],dist1,dist2,s1,r,hhr,xc1,xc2;
941: PetscInt k;
943: xc1 = user->xmin;
944: xc2 = xmid;
946: //VecSet(user->phi1,0.0);
947: for (k=0;k < 3; k++) {
948: if (x[k]-xqu > 0) {
949: s1 = (x[k] - xqu);
950: } else {
951: s1 = -(x[k] - xqu);
952: }
953: if (x[k] - xc1 > 0) {
954: dist1 = (x[k] - xc1);
955: } else {
956: dist1 = -(x[k] - xc1);
957: }
958: if (x[k] - xc2 > 0) {
959: dist2 = (x[k] - xc2);
960: } else {
961: dist2 = -(x[k] - xc2);
962: }
963: if (dist1 <= 0.5*lambda) {
964: r = (x[k]-xc1)/(0.5*lambda);
965: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
966: vals1[k] = hhr;
967: }
968: else if (dist2 <= 0.5*lambda) {
969: r = (x[k]-xc2)/(0.5*lambda);
970: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
971: vals1[k] = 1.0 - hhr;
972: }
973: else if (s1 <= xqu - 2.0*h) {
974: vals1[k] = 1.0;
975: }
976:
977: //else if ( abs(x[k]-(user->xmax-h)) < 0.1*h ) {
978: else if ( (user->xmax-h)-x[k] < 0.1*h ) {
979: vals1[k] = .15625;
980: }
981: else {
982: vals1[k] = 0.0;
983: }
984: }
985:
986: VecSetValuesLocal(user->phi1,3,idx,vals1,INSERT_VALUES);
988: xc1 = xmid;
989: xc2 = user->xmax;
991: //VecSet(user->phi2,0.0);
992: for (k=0;k < 3; k++) {
993: /*
994: s1 = abs(x[k] - (xqu+xmid));
995: dist1 = abs(x[k] - xc1);
996: dist2 = abs(x[k] - xc2);
997: */
998: if (x[k]-(xqu + xmid) > 0) {
999: s1 = (x[k] - (xqu + xmid));
1000: } else {
1001: s1 = -(x[k] - (xqu + xmid));
1002: }
1003: if (x[k] - xc1 > 0) {
1004: dist1 = (x[k] - xc1);
1005: } else {
1006: dist1 = -(x[k] - xc1);
1007: }
1008: if (x[k] - xc2 > 0) {
1009: dist2 = (x[k] - xc2);
1010: } else {
1011: dist2 = -(x[k] - xc2);
1012: }
1013:
1014: if (dist1 <= 0.5*lambda) {
1015: r = (x[k]-xc1)/(0.5*lambda);
1016: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
1017: vals2[k] = hhr;
1018: }
1019: else if (dist2 <= 0.5*lambda) {
1020: r = -(x[k]-xc2)/(0.5*lambda);
1021: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
1022: vals2[k] = hhr;
1023: }
1024: else if (s1 <= xqu - 2.0*h) {
1025: vals2[k] = 1.0;
1026: }
1027:
1028: else if ( x[k]-(user->xmin) < 0.1*h ) {
1029: vals2[k] = 0.5;
1030: }
1031:
1032:
1033: else if ( (x[k]-(user->xmin+h)) < 0.1*h ) {
1034: vals2[k] = .15625;
1035: }
1036:
1037: else {
1038: vals2[k] = 0.0;
1039: }
1040:
1041: }
1043: VecSetValuesLocal(user->phi2,3,idx,vals2,INSERT_VALUES);
1044: /*
1045: for (k=0;k < 3; k++) {
1046: vals_sum[k] = vals1[k]*vals1[k] + vals2[k]*vals2[k];
1047: }
1048: */
1049: //VecSetValuesLocal(user->Phi2D_V,3,idx,vals_sum,INSERT_VALUES);
1050:
1051: }
1052:
1053: VecAssemblyBegin(user->phi1);
1054: VecAssemblyEnd(user->phi1);
1055: VecAssemblyBegin(user->phi2);
1056: VecAssemblyEnd(user->phi2);
1058: /*
1059: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_phi",FILE_MODE_WRITE,&view);
1060: VecView(user->phi1,view);
1061: VecView(user->phi2,view);
1062: */
1063:
1064: //VecView(user->phi1,0);
1065: //VecView(user->phi2,0);
1066:
1067: VecPointwiseMult(user->phi1,user->phi1,user->phi1);
1068: VecPointwiseMult(user->phi2,user->phi2,user->phi2);
1069: /*
1070: VecView(user->phi1,view);
1071: VecView(user->phi2,view);
1072: */
1074: VecCopy(user->phi1,user->Phi2D_V);
1075: VecAXPY(user->Phi2D_V,1.0,user->phi2);
1076: //VecView(user->Phi2D_V,0);
1078: //VecView(user->Phi2D_V,view);
1079: //PetscViewerDestroy(&view);
1080:
1081: return(0);
1082:
1083: }
1087: PetscErrorCode Phi_read(AppCtx* user)
1088: {
1089: PetscErrorCode ierr;
1090: PetscReal *values;
1091: PetscViewer viewer;
1092: PetscInt power;
1095:
1096: power = user->darefine + (PetscInt)(PetscLogScalar(user->dagrid)/PetscLogScalar(2.0));
1097: switch (power) {
1098: case 6:
1099: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi4",FILE_MODE_READ,&viewer);
1100: VecLoad(user->Phi2D_V,viewer);
1101: PetscViewerDestroy(&viewer);
1102: break;
1103: case 7:
1104: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi3",FILE_MODE_READ,&viewer);
1105: VecLoad(user->Phi2D_V,viewer);
1106: PetscViewerDestroy(&viewer);
1107: break;
1108: case 8:
1109: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi2",FILE_MODE_READ,&viewer);
1110: VecLoad(user->Phi2D_V,viewer);
1111: PetscViewerDestroy(&viewer);
1112: break;
1113: case 9:
1114: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi1",FILE_MODE_READ,&viewer);
1115: VecLoad(user->Phi2D_V,viewer);
1116: PetscViewerDestroy(&viewer);
1117: break;
1118: case 10:
1119: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi",FILE_MODE_READ,&viewer);
1120: VecLoad(user->Phi2D_V,viewer);
1121: PetscViewerDestroy(&viewer);
1122: break;
1123: case 11:
1124: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim1",FILE_MODE_READ,&viewer);
1125: VecLoad(user->Phi2D_V,viewer);
1126: PetscViewerDestroy(&viewer);
1127: break;
1128: case 12:
1129: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim2",FILE_MODE_READ,&viewer);
1130: VecLoad(user->Phi2D_V,viewer);
1131: PetscViewerDestroy(&viewer);
1132: break;
1133: case 13:
1134: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim3",FILE_MODE_READ,&viewer);
1135: VecLoad(user->Phi2D_V,viewer);
1136: PetscViewerDestroy(&viewer);
1137: break;
1138: case 14:
1139: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim4",FILE_MODE_READ,&viewer);
1140: VecLoad(user->Phi2D_V,viewer);
1141: PetscViewerDestroy(&viewer);
1142: break;
1143: }
1144: return(0);
1145: }