Actual source code: ex65.c
petsc-3.5.4 2015-05-23
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,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,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE, 3, 1,NULL,NULL,&user.da1);
82: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE, 3, 1,NULL,NULL,&user.da1_clone);
83: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE, 1, 1,NULL,NULL,&user.da2);
84: } else {
85: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE, 3, 1,NULL,NULL,&user.da1);
86: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE, 3, 1,NULL,NULL,&user.da1_clone);
87: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_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: DMSetMatType(user.da1,MATAIJ);
127: DMCreateMatrix(user.da1,&user.M);
128: /* Get the (usual) mass matrix structure from da2 */
129: DMSetMatType(user.da2,MATAIJ);
130: DMCreateMatrix(user.da2,&user.M_0);
131: /* Form the jacobian matrix and M_0 */
132: SetUpMatrices(&user);
133: SetInitialGuess(x,&user);
134: MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
136: SNESCreate(PETSC_COMM_WORLD,&snes);
137: SNESSetDM(snes,user.da1);
139: SNESSetFunction(snes,r,FormFunction,(void*)&user);
140: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
142: SetVariableBounds(user.da1,xl,xu);
143: SNESVISetVariableBounds(snes,xl,xu);
144: SNESSetFromOptions(snes);
145: /* SNESVISetRedundancyCheck(snes,(PetscErrorCode (*)(SNES,IS,IS*,void*))CheckRedundancy,user.da1_clone); */
146: /*
147: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
148: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_p",FILE_MODE_WRITE,&view_p);
149: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
150: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);
151: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);
152: */
153: /* PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds); */
156: if (user.graphics) {
157: /*PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds);*/
159: VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
160: }
163: /* multidomain modeling */
164: if (user.domain) {
165: switch (user.domain) {
166: case 1:
167: Phi(&user);
168: break;
169: case 2:
170: Phi_read(&user);
171: break;
172: }
173: }
175: while (t<user.T) {
177: char filename[PETSC_MAX_PATH_LEN];
178: PetscScalar a = 1.0;
179: PetscInt i;
180: /*PetscViewer view;*/
183: SNESSetFunction(snes,r,FormFunction,(void*)&user);
184: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
186: if (user.radiation) {
187: SetRandomVectors(&user,t);
188: }
189: DPsi(&user);
190: /*VecView(user.DPsiv,view_psi);*/
191: /*VecView(user.DPsieta,view_psi);*/
193: Update_q(&user);
194: /*VecView(user.q,view_q);*/
195: /*MatView(user.M,view_mat);*/
197: SNESSolve(snes,NULL,x);
198: /*VecView(x,view_out);*/
201: if (user.graphics) {
202: VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
203: }
205: PetscInt its;
206: SNESGetIterationNumber(snes,&its);
207: PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4g in %d iterations\n",t,its);
209: Update_u(x,&user);
210: /*
211: for (i=0; i < (int)(user.T/a) ; i++) {
212: if (t/a > i - user.dt/a && t/a < i + user.dt/a) {
213: sprintf(filename,"output_%f",t);
214: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&view);
215: VecView(user.cv,view);
216: VecView(user.eta,view);
217: PetscViewerDestroy(&view);
218: }
219: }
220: */
222: t = t + user.dt;
224: }
226: /*
227: PetscViewerDestroy(&view_out);
228: PetscViewerDestroy(&view_p);
229: PetscViewerDestroy(&view_q);
230: PetscViewerDestroy(&view_psi);
231: PetscViewerDestroy(&view_mat);
232: */
233: VecDestroy(&x);
234: VecDestroy(&r);
235: VecDestroy(&xl);
236: VecDestroy(&xu);
237: VecDestroy(&user.q);
238: VecDestroy(&user.wv);
239: VecDestroy(&user.cv);
240: VecDestroy(&user.eta);
241: VecDestroy(&user.DPsiv);
242: VecDestroy(&user.DPsieta);
243: VecDestroy(&user.Pv);
244: VecDestroy(&user.Pi);
245: VecDestroy(&user.Piv);
246: VecDestroy(&user.logcv);
247: VecDestroy(&user.logcv2);
248: VecDestroy(&user.work1);
249: VecDestroy(&user.work2);
250: /* for twodomain modeling */
251: VecDestroy(&user.phi1);
252: VecDestroy(&user.phi2);
253: VecDestroy(&user.Phi2D_V);
254: VecDestroy(&user.Sv);
255: VecDestroy(&user.Si);
256: VecDestroy(&user.work3);
257: VecDestroy(&user.work4);
259: MatDestroy(&user.M);
260: MatDestroy(&user.M_0);
261: DMDestroy(&user.da1);
262: DMDestroy(&user.da1_clone);
263: DMDestroy(&user.da2);
264: SNESDestroy(&snes);
265: PetscFinalize();
266: return 0;
267: }
271: PetscErrorCode Update_u(Vec X,AppCtx *user)
272: {
274: PetscInt i,n;
275: PetscScalar *xx,*wv_p,*cv_p,*eta_p;
278: VecGetLocalSize(user->wv,&n);
280: VecGetArray(X,&xx);
281: VecGetArray(user->wv,&wv_p);
282: VecGetArray(user->cv,&cv_p);
283: VecGetArray(user->eta,&eta_p);
286: for (i=0; i<n; i++) {
287: wv_p[i] = xx[3*i];
288: cv_p[i] = xx[3*i+1];
289: eta_p[i] = xx[3*i+2];
290: }
291: VecRestoreArray(X,&xx);
292: VecRestoreArray(user->wv,&wv_p);
293: VecRestoreArray(user->cv,&cv_p);
294: VecRestoreArray(user->eta,&eta_p);
295: return(0);
296: }
300: PetscErrorCode Update_q(AppCtx *user)
301: {
303: PetscScalar *q_p, *w1, *w2;
304: PetscInt n,i;
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++) q_p[3*i]=w2[i];
335: MatMult(user->M_0,user->DPsiv,user->work1);
336: for (i=0; i<n; i++) q_p[3*i+1]=w1[i];
338: VecCopy(user->DPsieta,user->work1);
339: VecScale(user->work1,user->L*user->dt);
340: VecAXPY(user->work1,-1.0,user->eta);
341: if (user->radiation) {
342: VecAXPY(user->work1,-1.0,user->Piv);
343: }
344: MatMult(user->M_0,user->work1,user->work2);
345: for (i=0; i<n; i++) q_p[3*i+2]=w2[i];
347: VecRestoreArray(user->q,&q_p);
348: VecRestoreArray(user->work1,&w1);
349: VecRestoreArray(user->work2,&w2);
350: return(0);
351: }
355: PetscErrorCode DPsi(AppCtx *user)
356: {
358: PetscScalar Evf=user->Evf,A=user->A,B=user->B,cv0=user->cv0;
359: PetscScalar *cv_p,*eta_p,*logcv_p,*logcv2_p,*DPsiv_p,*DPsieta_p;
360: PetscInt n,i;
363: VecGetLocalSize(user->cv,&n);
364: VecGetArray(user->cv,&cv_p);
365: VecGetArray(user->eta,&eta_p);
366: VecGetArray(user->logcv,&logcv_p);
367: VecGetArray(user->logcv2,&logcv2_p);
368: VecGetArray(user->DPsiv,&DPsiv_p);
369: VecGetArray(user->DPsieta,&DPsieta_p);
371: Llog(user->cv,user->logcv);
373: VecCopy(user->cv,user->work1);
374: VecScale(user->work1,-1.0);
375: VecShift(user->work1,1.0);
376: Llog(user->work1,user->logcv2);
378: for (i=0; i<n; i++)
379: {
380: 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];
382: 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];
384: }
386: VecRestoreArray(user->cv,&cv_p);
387: VecRestoreArray(user->eta,&eta_p);
388: VecGetArray(user->logcv,&logcv_p);
389: VecGetArray(user->logcv2,&logcv2_p);
390: VecRestoreArray(user->DPsiv,&DPsiv_p);
391: VecRestoreArray(user->DPsieta,&DPsieta_p);
392: return(0);
393: }
398: PetscErrorCode Llog(Vec X, Vec Y)
399: {
401: PetscScalar *x,*y;
402: PetscInt n,i;
405: VecGetArray(X,&x);
406: VecGetArray(Y,&y);
407: VecGetLocalSize(X,&n);
408: for (i=0; i<n; i++) {
409: if (x[i] < 1.0e-12) y[i] = log(1.0e-12);
410: else y[i] = log(x[i]);
411: }
412: return(0);
413: }
417: PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
418: {
419: PetscErrorCode ierr;
420: PetscInt n,i,Mda,Nda;
421: PetscScalar *xx,*cv_p,*wv_p,*eta_p;
422: /*PetscViewer view_out;*/
423: /* needed for the void growth case */
424: PetscScalar xmid,ymid,cv_v=1.0,cv_m=user->Sv_scalar*user->cv0,eta_v=1.0,eta_m=0.0,h,lambda;
425: PetscInt nele,nen,idx[3];
426: const PetscInt *ele;
427: PetscScalar x[3],y[3];
428: Vec coords;
429: const PetscScalar *_coords;
430: PetscScalar xwidth = user->xmax - user->xmin, ywidth = user->ymax - user->ymin;
433: VecGetLocalSize(X,&n);
435: if (!user->radiation) {
436: DMDAGetInfo(user->da2,NULL,&Mda,&Nda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
437: DMGetCoordinatesLocal(user->da2,&coords);
438: VecGetArrayRead(coords,&_coords);
440: if (user->periodic) h = (user->xmax-user->xmin)/Mda;
441: else h = (user->xmax-user->xmin)/(Mda-1.0);
443: xmid = (user->xmax + user->xmin)/2.0;
444: ymid = (user->ymax + user->ymin)/2.0;
445: lambda = 4.0*h;
447: DMDAGetElements(user->da2,&nele,&nen,&ele);
448: for (i=0; i < nele; i++) {
449: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
451: x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
452: x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
453: x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
455: PetscInt k;
456: PetscScalar vals_DDcv[3],vals_cv[3],vals_eta[3],s,hhr,r;
457: for (k=0; k < 3; k++) {
458: /*s = PetscAbs(x[k] - xmid);*/
459: s = PetscSqrtScalar((x[k] - xmid)*(x[k] - xmid) + (y[k] - ymid)*(y[k] - ymid));
460: if (s <= xwidth*(5.0/64.0)) {
461: vals_cv[k] = cv_v;
462: vals_eta[k] = eta_v;
463: vals_DDcv[k] = 0.0;
464: } else if (s> xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0)) {
465: /*r = (s - xwidth*(6.0/64.0))/(0.5*lambda);*/
466: r = (s - xwidth*(6.0/64.0))/(xwidth/64.0);
467: hhr = 0.25*(-r*r*r + 3*r + 2);
468: vals_cv[k] = cv_m + (1.0 - hhr)*(cv_v - cv_m);
469: vals_eta[k] = eta_m + (1.0 - hhr)*(eta_v - eta_m);
470: vals_DDcv[k] = (cv_v - cv_m)*r*6.0/(lambda*lambda);
471: } else {
472: vals_cv[k] = cv_m;
473: vals_eta[k] = eta_m;
474: vals_DDcv[k] = 0.0;
475: }
476: }
478: VecSetValuesLocal(user->cv,3,idx,vals_cv,INSERT_VALUES);
479: VecSetValuesLocal(user->eta,3,idx,vals_eta,INSERT_VALUES);
480: VecSetValuesLocal(user->work2,3,idx,vals_DDcv,INSERT_VALUES);
482: }
483: VecAssemblyBegin(user->cv);
484: VecAssemblyEnd(user->cv);
485: VecAssemblyBegin(user->eta);
486: VecAssemblyEnd(user->eta);
487: VecAssemblyBegin(user->work2);
488: VecAssemblyEnd(user->work2);
490: DMDARestoreElements(user->da2,&nele,&nen,&ele);
491: VecRestoreArrayRead(coords,&_coords);
492: } else {
493: /*VecSet(user->cv,user->initv);*/
494: /*VecSet(user->ci,user->initv);*/
495: VecSet(user->cv,.05);
496: VecSet(user->eta,user->initeta);
498: }
499: DPsi(user);
500: VecCopy(user->DPsiv,user->wv);
501: VecAXPY(user->wv,-2.0*user->kav,user->work2);
503: VecGetArray(X,&xx);
504: VecGetArray(user->wv,&wv_p);
505: VecGetArray(user->cv,&cv_p);
506: VecGetArray(user->eta,&eta_p);
508: for (i=0; i<n/3; i++)
509: {
510: xx[3*i] =wv_p[i];
511: xx[3*i+1]=cv_p[i];
512: xx[3*i+2]=eta_p[i];
513: }
515: /*
516: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);
517: VecView(user->wv,view_out);
518: VecView(user->cv,view_out);
519: VecView(user->eta,view_out);
520: PetscViewerDestroy(&view_out);
521: */
523: VecRestoreArray(X,&xx);
524: VecRestoreArray(user->wv,&wv_p);
525: VecRestoreArray(user->cv,&cv_p);
526: VecRestoreArray(user->eta,&eta_p);
527: return(0);
528: }
530: typedef struct {
531: PetscReal dt,x,y,strength;
532: } RandomValues;
537: PetscErrorCode SetRandomVectors(AppCtx *user,PetscReal t)
538: {
539: PetscErrorCode ierr;
540: static RandomValues *randomvalues = 0;
541: static PetscInt randindex = 0,n; /* indicates how far into the randomvalues we have currently used */
542: static PetscReal randtime = 0; /* indicates time of last radiation event */
543: PetscInt i,j,M,N,cnt = 0;
544: PetscInt xs,ys,xm,ym;
547: if (!randomvalues) {
548: PetscViewer viewer;
549: char filename[PETSC_MAX_PATH_LEN];
550: PetscBool flg;
551: PetscInt seed;
553: PetscOptionsGetInt(NULL,"-random_seed",&seed,&flg);
554: if (flg) {
555: sprintf(filename,"ex61.random.%d",(int)seed);
556: } else {
557: PetscStrcpy(filename,"ex61.random");
558: }
559: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
560: PetscViewerBinaryRead(viewer,&n,1,PETSC_INT);
561: PetscMalloc1(n,&randomvalues);
562: PetscViewerBinaryRead(viewer,randomvalues,4*n,PETSC_DOUBLE);
563: for (i=0; i<n; i++) randomvalues[i].dt = randomvalues[i].dt*user->dtevent;
564: PetscViewerDestroy(&viewer);
565: }
566: user->maxevents = PetscMin(user->maxevents,n);
568: VecSet(user->Pv,0.0);
569: DMDAGetInfo(user->da1,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
570: DMDAGetGhostCorners(user->da1,&xs,&ys,0,&xm,&ym,0);
571: while (user->maxevents > randindex && randtime + randomvalues[randindex].dt < t + user->dt) { /* radiation event has occured since last time step */
572: i = ((PetscInt) (randomvalues[randindex].x*M)) - xs;
573: j = ((PetscInt) (randomvalues[randindex].y*N)) - ys;
574: if (i >= 0 && i < xm && j >= 0 && j < ym) { /* point is on this process */
576: /* need to make sure eta at the given point is not great than .8 */
577: VecSetValueLocal(user->Pv,i + xm*(j), randomvalues[randindex].strength*user->VG,INSERT_VALUES);
578: }
579: randtime += randomvalues[randindex++].dt;
580: cnt++;
581: }
582: PetscPrintf(PETSC_COMM_WORLD,"Number of radiation events %d\n",cnt);
583: VecAssemblyBegin(user->Pv);
584: VecAssemblyEnd(user->Pv);
586: VecCopy(user->Pv,user->Pi);
587: VecScale(user->Pi,0.9);
588: VecPointwiseMult(user->Piv,user->Pi,user->Pv);
589: return(0);
590: }
595: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
596: {
598: AppCtx *user=(AppCtx*)ctx;
601: MatMultAdd(user->M,X,user->q,F);
602: return(0);
603: }
607: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat B,void *ctx)
608: {
610: AppCtx *user=(AppCtx*)ctx;
613: MatCopy(user->M,J,*flg);
614: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
615: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
616: return(0);
617: }
620: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
621: {
623: PetscScalar ***l,***u;
624: PetscInt xs,xm,ys,ym;
625: PetscInt i,j;
628: DMDAVecGetArrayDOF(da,xl,&l);
629: DMDAVecGetArrayDOF(da,xu,&u);
631: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
633: for (j=ys; j < ys+ym; j++) {
634: for (i=xs; i < xs+xm; i++) {
635: l[j][i][0] = -PETSC_INFINITY;
636: l[j][i][1] = 0.0;
637: l[j][i][2] = 0.0;
638: u[j][i][0] = PETSC_INFINITY;
639: u[j][i][1] = 1.0;
640: u[j][i][2] = 1.0;
641: }
642: }
644: DMDAVecRestoreArrayDOF(da,xl,&l);
645: DMDAVecRestoreArrayDOF(da,xu,&u);
646: return(0);
647: }
651: PetscErrorCode GetParams(AppCtx *user)
652: {
654: PetscBool flg;
657: /* Set default parameters */
658: user->xmin = 0.0; user->xmax = 128.0;
659: user->ymin = 0.0; user->ymax = 128.0;
660: user->Mv = 1.0;
661: user->L = 1.0;
662: user->kaeta = 1.0;
663: user->kav = 0.5;
664: user->Evf = 9.09;
665: user->A = 9.09;
666: user->B = 9.09;
667: user->cv0 = 1.13e-4;
668: user->Sv_scalar = 500.0;
669: user->dt = 1.0e-5;
670: user->T = 1.0e-2;
671: user->initv = .00069;
672: user->initeta = .0;
673: user->graphics = PETSC_FALSE;
674: user->periodic = PETSC_FALSE;
675: user->lumpedmass = PETSC_FALSE;
676: user->radiation = PETSC_FALSE;
677: /* multidomain modeling */
678: user->domain = 0;
679: user->grain = 5000.0;
680: user->Svr = 0.5;
681: user->Sir = 0.5;
682: user->cv_eq = 6.9e-4;
683: user->ci_eq = 6.9e-4;
684: user->VG = 100.0;
685: user->maxevents = 10;
687: user->dtevent = user->dt;
688: PetscOptionsReal("-dtevent","Average time between random events\n","None",user->dtevent,&user->dtevent,&flg);
691: PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);
692: PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);
693: PetscOptionsGetReal(NULL,"-T",&user->T,&flg);
694: PetscOptionsGetReal(NULL,"-dt",&user->dt,&flg);
695: PetscOptionsBool("-graphics","Contour plot solutions at each timestep\n","None",user->graphics,&user->graphics,&flg);
696: PetscOptionsBool("-periodic","Use periodic boundary conditions\n","None",user->periodic,&user->periodic,&flg);
697: PetscOptionsBool("-radiation","Allow radiation\n","None",user->radiation,&user->radiation,&flg);
698: PetscOptionsBool("-lumpedmass","Use lumped mass matrix\n","None",user->lumpedmass,&user->lumpedmass,&flg);
699: PetscOptionsInt("-domain","Number of domains (0=one domain, 1=two domains, 2=multidomain\n","None",user->domain,&user->domain,&flg);
700: PetscOptionsReal("-VG","Maximum increase in vacancy (or interstitial) concentration due to a cascade event","None",user->VG,&user->VG,&flg);
701: PetscOptionsReal("-grain","Some bogus factor that controls the strength of the grain boundaries, makes the picture more plausible","None",user->grain,&user->grain,&flg);
702: PetscOptionsInt("-maxevents","Maximum random events allowed\n","None",user->maxevents,&user->maxevents,&flg);
703: PetscOptionsInt("-da_refine","da refine \n","None",user->darefine,&user->darefine,&flg);
704: PetscOptionsInt("-da_grid_x","da grid x\n","None",user->dagrid,&user->dagrid,&flg);
705: return(0);
706: }
711: PetscErrorCode SetUpMatrices(AppCtx *user)
712: {
714: PetscInt nele,nen,i,n;
715: const PetscInt *ele;
716: PetscScalar dt=user->dt,h;
718: PetscInt idx[3];
719: PetscScalar eM_0[3][3],eM_2[3][3];
720: Mat M =user->M;
721: Mat M_0=user->M_0;
722: PetscInt Mda,Nda;
726: MatGetLocalSize(M,&n,NULL);
727: DMDAGetInfo(user->da1,NULL,&Mda,&Nda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
729: if (Mda!=Nda) printf("Currently different Mda and Nda are not supported");
730: if (user->periodic) h = (user->xmax-user->xmin)/Mda;
731: else h = (user->xmax-user->xmin)/(Mda-1.0);
733: if (user->lumpedmass) {
734: eM_0[0][0] = eM_0[1][1] = eM_0[2][2] = h*h/6.0;
735: 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;
736: } else {
737: eM_0[0][0] = eM_0[1][1] = eM_0[2][2] = h*h/12.0;
738: 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;
739: }
740: eM_2[0][0] = 1.0;
741: eM_2[1][1] = eM_2[2][2] = 0.5;
742: eM_2[0][1] = eM_2[0][2] = eM_2[1][0]= eM_2[2][0] = -0.5;
743: eM_2[1][2] = eM_2[2][1] = 0.0;
746: /* Get local element info */
747: DMDAGetElements(user->da1,&nele,&nen,&ele);
748: for (i=0; i < nele; i++) {
750: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
752: PetscInt row,cols[6],r,row_M_0,cols3[3];
753: PetscScalar vals[6],vals_M_0[3],vals3[3];
755: for (r=0; r<3; r++) {
756: row_M_0 = idx[r];
757: vals_M_0[0]=eM_0[r][0];
758: vals_M_0[1]=eM_0[r][1];
759: vals_M_0[2]=eM_0[r][2];
761: MatSetValuesLocal(M_0,1,&row_M_0,3,idx,vals_M_0,ADD_VALUES);
763: row = 3*idx[r];
764: cols[0] = 3*idx[0]; vals[0] = dt*eM_2[r][0]*user->Mv;
765: cols[1] = 3*idx[1]; vals[1] = dt*eM_2[r][1]*user->Mv;
766: cols[2] = 3*idx[2]; vals[2] = dt*eM_2[r][2]*user->Mv;
767: cols[3] = 3*idx[0]+1; vals[3] = eM_0[r][0];
768: cols[4] = 3*idx[1]+1; vals[4] = eM_0[r][1];
769: cols[5] = 3*idx[2]+1; vals[5] = eM_0[r][2];
771: /* Insert values in matrix M for 1st dof */
772: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
774: row = 3*idx[r]+1;
775: cols[0] = 3*idx[0]; vals[0] = -eM_0[r][0];
776: cols[1] = 3*idx[1]; vals[1] = -eM_0[r][1];
777: cols[2] = 3*idx[2]; vals[2] = -eM_0[r][2];
778: cols[3] = 3*idx[0]+1; vals[3] = 2.0*user->kav*eM_2[r][0];
779: cols[4] = 3*idx[1]+1; vals[4] = 2.0*user->kav*eM_2[r][1];
780: cols[5] = 3*idx[2]+1; vals[5] = 2.0*user->kav*eM_2[r][2];
782: /* Insert values in matrix M for 2nd dof */
783: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
785: row = 3*idx[r]+2;
786: cols3[0] = 3*idx[0]+2; vals3[0] = eM_0[r][0] + user->dt*2.0*user->L*user->kaeta*eM_2[r][0];
787: cols3[1] = 3*idx[1]+2; vals3[1] = eM_0[r][1] + user->dt*2.0*user->L*user->kaeta*eM_2[r][1];
788: cols3[2] = 3*idx[2]+2; vals3[2] = eM_0[r][2] + user->dt*2.0*user->L*user->kaeta*eM_2[r][2];
790: MatSetValuesLocal(M,1,&row,3,cols3,vals3,ADD_VALUES);
791: }
792: }
794: MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
795: MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);
797: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
798: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
800: DMDARestoreElements(user->da1,&nele,&nen,&ele);
801: return(0);
802: }
806: PetscErrorCode CheckRedundancy(SNES snes, IS act, IS *outact, DM da)
807: {
809: PetscScalar **uin,**uout;
810: Vec UIN, UOUT;
811: PetscInt xs,xm,*outindex;
812: const PetscInt *index;
813: PetscInt k,i,l,n,M,cnt=0;
816: DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
817: DMGetGlobalVector(da,&UIN);
818: VecSet(UIN,0.0);
819: DMGetLocalVector(da,&UOUT);
820: DMDAVecGetArrayDOF(da,UIN,&uin);
821: ISGetIndices(act,&index);
822: ISGetLocalSize(act,&n);
823: for (k=0; k<n; k++) {
824: l = index[k]%5;
825: i = index[k]/5;
827: uin[i][l]=1.0;
828: }
829: printf("Number of active constraints before applying redundancy %d\n",n);
830: ISRestoreIndices(act,&index);
831: DMDAVecRestoreArrayDOF(da,UIN,&uin);
832: DMGlobalToLocalBegin(da,UIN,INSERT_VALUES,UOUT);
833: DMGlobalToLocalEnd(da,UIN,INSERT_VALUES,UOUT);
834: DMDAVecGetArrayDOF(da,UOUT,&uout);
836: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
838: for (i=xs; i < xs+xm;i++) {
839: if (uout[i-1][1] && uout[i][1] && uout[i+1][1]) uout[i][0] = 1.0;
840: if (uout[i-1][3] && uout[i][3] && uout[i+1][3]) uout[i][2] = 1.0;
841: }
843: for (i=xs; i < xs+xm;i++) {
844: for (l=0;l<5;l++) {
845: if (uout[i][l]) cnt++;
846: }
847: }
849: printf("Number of active constraints after applying redundancy %d\n",cnt);
852: PetscMalloc1(cnt,&outindex);
853: cnt = 0;
855: for (i=xs; i < xs+xm;i++) {
856: for (l=0;l<5;l++) {
857: if (uout[i][l]) outindex[cnt++] = 5*(i)+l;
858: }
859: }
862: ISCreateGeneral(PETSC_COMM_WORLD,cnt,outindex,PETSC_OWN_POINTER,outact);
863: DMDAVecRestoreArrayDOF(da,UOUT,&uout);
864: DMRestoreGlobalVector(da,&UIN);
865: DMRestoreLocalVector(da,&UOUT);
866: return(0);
867: }
871: PetscErrorCode Phi(AppCtx *user)
872: {
873: PetscErrorCode ierr;
874: PetscScalar xmid, xqu, lambda, h,x[3],y[3];
875: Vec coords;
876: const PetscScalar *_coords;
877: PetscInt nele,nen,i,idx[3],Mda,Nda;
878: const PetscInt *ele;
879: /*PetscViewer view;*/
882: DMDAGetInfo(user->da1,NULL,&Mda,&Nda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
883: DMGetCoordinatesLocal(user->da2,&coords);
884: VecGetArrayRead(coords,&_coords);
886: h = (user->xmax - user->xmin)/Mda;
887: xmid = (user->xmin + user->xmax)/2.0;
888: xqu = (user->xmin + xmid)/2.0;
889: lambda = 4.0*h;
892: DMDAGetElements(user->da2,&nele,&nen,&ele);
893: for (i=0; i < nele; i++) {
894: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
895: /*printf("idx[0]=%d,idx[1]=%d,idx[2]=%d\n",idx[0],idx[1],idx[2]);*/
897: x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
898: x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
899: x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
901: /*printf("x[0]=%f,x[1]=%f,x[2]=%f\n",x[0],x[1],x[2]);*/
902: /*printf("y[0]=%f,y[1]=%f,y[2]=%f\n",y[0],y[1],y[2]);*/
904: PetscScalar vals1[3],vals2[3],dist1,dist2,s1,r,hhr,xc1,xc2;
905: PetscInt k;
907: xc1 = user->xmin;
908: xc2 = xmid;
910: /*VecSet(user->phi1,0.0);*/
911: for (k=0;k < 3; k++) {
912: if (x[k]-xqu > 0) s1 = (x[k] - xqu);
913: else s1 = -(x[k] - xqu);
915: if (x[k] - xc1 > 0) dist1 = (x[k] - xc1);
916: else dist1 = -(x[k] - xc1);
918: if (x[k] - xc2 > 0) dist2 = (x[k] - xc2);
919: else dist2 = -(x[k] - xc2);
921: if (dist1 <= 0.5*lambda) {
922: r = (x[k]-xc1)/(0.5*lambda);
923: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
924: vals1[k] = hhr;
925: } else if (dist2 <= 0.5*lambda) {
926: r = (x[k]-xc2)/(0.5*lambda);
927: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
928: vals1[k] = 1.0 - hhr;
929: } else if (s1 <= xqu - 2.0*h) {
930: vals1[k] = 1.0;
931: } else if ((user->xmax-h)-x[k] < 0.1*h) {
932: /*else if (abs(x[k]-(user->xmax-h)) < 0.1*h) {*/
933: vals1[k] = .15625;
934: } else {
935: vals1[k] = 0.0;
936: }
937: }
939: VecSetValuesLocal(user->phi1,3,idx,vals1,INSERT_VALUES);
941: xc1 = xmid;
942: xc2 = user->xmax;
944: /*VecSet(user->phi2,0.0);*/
945: for (k=0;k < 3; k++) {
946: /*
947: s1 = abs(x[k] - (xqu+xmid));
948: dist1 = abs(x[k] - xc1);
949: dist2 = abs(x[k] - xc2);
950: */
951: if (x[k]-(xqu + xmid) > 0) s1 = (x[k] - (xqu + xmid));
952: else s1 = -(x[k] - (xqu + xmid));
954: if (x[k] - xc1 > 0) dist1 = (x[k] - xc1);
955: else dist1 = -(x[k] - xc1);
957: if (x[k] - xc2 > 0) dist2 = (x[k] - xc2);
958: else dist2 = -(x[k] - xc2);
960: if (dist1 <= 0.5*lambda) {
961: r = (x[k]-xc1)/(0.5*lambda);
962: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
963: vals2[k] = hhr;
964: } else if (dist2 <= 0.5*lambda) {
965: r = -(x[k]-xc2)/(0.5*lambda);
966: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
967: vals2[k] = hhr;
968: } else if (s1 <= xqu - 2.0*h) {
969: vals2[k] = 1.0;
970: } else if (x[k]-(user->xmin) < 0.1*h) {
971: vals2[k] = 0.5;
972: } else if ((x[k]-(user->xmin+h)) < 0.1*h) {
973: vals2[k] = .15625;
974: } else {
975: vals2[k] = 0.0;
976: }
978: }
980: VecSetValuesLocal(user->phi2,3,idx,vals2,INSERT_VALUES);
981: /*
982: for (k=0;k < 3; k++) {
983: vals_sum[k] = vals1[k]*vals1[k] + vals2[k]*vals2[k];
984: }
985: */
986: /*VecSetValuesLocal(user->Phi2D_V,3,idx,vals_sum,INSERT_VALUES);*/
988: }
990: VecAssemblyBegin(user->phi1);
991: VecAssemblyEnd(user->phi1);
992: VecAssemblyBegin(user->phi2);
993: VecAssemblyEnd(user->phi2);
995: /*
996: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_phi",FILE_MODE_WRITE,&view);
997: VecView(user->phi1,view);
998: VecView(user->phi2,view);
999: */
1001: /*VecView(user->phi1,0);*/
1002: /*VecView(user->phi2,0);*/
1004: VecPointwiseMult(user->phi1,user->phi1,user->phi1);
1005: VecPointwiseMult(user->phi2,user->phi2,user->phi2);
1006: /*
1007: VecView(user->phi1,view);
1008: VecView(user->phi2,view);
1009: */
1011: VecCopy(user->phi1,user->Phi2D_V);
1012: VecAXPY(user->Phi2D_V,1.0,user->phi2);
1013: /*VecView(user->Phi2D_V,0);*/
1015: /*VecView(user->Phi2D_V,view);*/
1016: /*PetscViewerDestroy(&view);*/
1017: return(0);
1018: }
1022: PetscErrorCode Phi_read(AppCtx *user)
1023: {
1025: PetscReal *values;
1026: PetscViewer viewer;
1027: PetscInt power;
1030: power = user->darefine + (PetscInt)(PetscLogScalar(user->dagrid)/PetscLogScalar(2.0));
1031: switch (power) {
1032: case 6:
1033: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi4",FILE_MODE_READ,&viewer);
1034: VecLoad(user->Phi2D_V,viewer);
1035: PetscViewerDestroy(&viewer);
1036: break;
1037: case 7:
1038: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi3",FILE_MODE_READ,&viewer);
1039: VecLoad(user->Phi2D_V,viewer);
1040: PetscViewerDestroy(&viewer);
1041: break;
1042: case 8:
1043: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi2",FILE_MODE_READ,&viewer);
1044: VecLoad(user->Phi2D_V,viewer);
1045: PetscViewerDestroy(&viewer);
1046: break;
1047: case 9:
1048: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi1",FILE_MODE_READ,&viewer);
1049: VecLoad(user->Phi2D_V,viewer);
1050: PetscViewerDestroy(&viewer);
1051: break;
1052: case 10:
1053: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi",FILE_MODE_READ,&viewer);
1054: VecLoad(user->Phi2D_V,viewer);
1055: PetscViewerDestroy(&viewer);
1056: break;
1057: case 11:
1058: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim1",FILE_MODE_READ,&viewer);
1059: VecLoad(user->Phi2D_V,viewer);
1060: PetscViewerDestroy(&viewer);
1061: break;
1062: case 12:
1063: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim2",FILE_MODE_READ,&viewer);
1064: VecLoad(user->Phi2D_V,viewer);
1065: PetscViewerDestroy(&viewer);
1066: break;
1067: case 13:
1068: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim3",FILE_MODE_READ,&viewer);
1069: VecLoad(user->Phi2D_V,viewer);
1070: PetscViewerDestroy(&viewer);
1071: break;
1072: case 14:
1073: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim4",FILE_MODE_READ,&viewer);
1074: VecLoad(user->Phi2D_V,viewer);
1075: PetscViewerDestroy(&viewer);
1076: break;
1077: }
1078: return(0);
1079: }