Actual source code: ex63.c
petsc-3.3-p7 2013-05-11
1: static char help[] = "1D coupled Allen-Cahn and Cahn-Hilliard equation for degenerate mobility and triangular elements.\n\
2: Runtime options include:\n\
3: -xmin <xmin>\n\
4: -xmax <xmax>\n\
5: -ymin <ymin>\n\
6: -T <T>, where <T> is the end time for the time domain simulation\n\
7: -dt <dt>,where <dt> is the step size for the numerical integration\n\
8: -gamma <gamma>\n\
9: -theta_c <theta_c>\n\n";
11: /*
12: ./ex63 -ksp_type fgmres -snes_vi_monitor -snes_atol 1.e-11 -snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor -VG 10000000 -draw_fields 1,3,4 -pc_type mg -pc_mg_galerkin -log_summary -dt .000001 -mg_coarse_pc_type svd -ksp_monitor_true_residual -ksp_rtol 1.e-9 -snes_linesearch basic -T .0020 -P_casc .0005 -snes_monitor_solution -da_refine 10
13: */
15: #include petscsnes.h
16: #include petscdmda.h
18: typedef struct{
19: PetscReal dt,T; /* Time step and end time */
20: DM da1,da1_clone,da2;
21: Mat M; /* Jacobian matrix */
22: Mat M_0;
23: Vec q,wv,cv,wi,ci,eta,cvi,DPsiv,DPsii,DPsieta,Pv,Pi,Piv,logcv,logci,logcvi,Riv;
24: Vec work1,work2,work3,work4;
25: PetscScalar Dv,Di,Evf,Eif,A,kBT,kav,kai,kaeta,Rsurf,Rbulk,L,P_casc,VG; /* physics parameters */
26: PetscReal xmin,xmax,ymin,ymax;
27: PetscInt nx;
28: PetscBool voidgrowth;
29: PetscBool degenerate;
30: PetscBool graphics;
31: PetscReal smallnumber;
32: PetscReal initv;
33: PetscReal initeta;
34: }AppCtx;
36: PetscErrorCode GetParams(AppCtx*);
37: PetscErrorCode SetRandomVectors(AppCtx*);
38: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
39: PetscErrorCode SetUpMatrices(AppCtx*);
40: PetscErrorCode UpdateMatrices(AppCtx*);
41: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
42: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
43: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
44: PetscErrorCode Update_q(AppCtx*);
45: PetscErrorCode Update_u(Vec,AppCtx*);
46: PetscErrorCode DPsi(AppCtx*);
47: PetscErrorCode Llog(Vec,Vec);
48: PetscErrorCode CheckRedundancy(SNES,IS,IS*,DM);
52: int main(int argc, char **argv)
53: {
55: Vec x,r; /* olution and residual vectors */
56: SNES snes; /* Nonlinear solver context */
57: AppCtx user; /* Application context */
58: Vec xl,xu; /* Upper and lower bounds on variables */
59: Mat J;
60: PetscScalar t=0.0;
61: PetscViewer view_out, view_p, view_q, view_psi, view_mat;
62: PetscReal bounds[] = {1000.0,-1000.,0.0,1.0,1000.0,-1000.0,0.0,1.0,1000.0,-1000.0};
65: PetscInitialize(&argc,&argv, (char*)0, help);
66:
67: /* Get physics and time parameters */
68: GetParams(&user);
69: /* Create a 1D DA with dof = 5; the whole thing */
70: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC, -4, 5, 1,PETSC_NULL,&user.da1);
71: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC, -4, 5, 1,PETSC_NULL,&user.da1_clone);
72: /* Create a 1D DA with dof = 1; for individual componentes */
73: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC, -4, 1, 1,PETSC_NULL,&user.da2);
75: /* Set Element type (triangular) */
76: DMDASetElementType(user.da1,DMDA_ELEMENT_P1);
77: DMDASetElementType(user.da2,DMDA_ELEMENT_P1);
78:
79: /* Set x and y coordinates */
80: DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
81: DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
82: /* Get global vector x from DM (da1) and duplicate vectors r,xl,xu */
83: DMCreateGlobalVector(user.da1,&x);
84: VecDuplicate(x,&r);
85: VecDuplicate(x,&xl);
86: VecDuplicate(x,&xu);
87: VecDuplicate(x,&user.q);
88:
89: /* Get global vector user->wv from da2 and duplicate other vectors */
90: DMCreateGlobalVector(user.da2,&user.wv);
91: VecDuplicate(user.wv,&user.cv);
92: VecDuplicate(user.wv,&user.wi);
93: VecDuplicate(user.wv,&user.ci);
94: VecDuplicate(user.wv,&user.eta);
95: VecDuplicate(user.wv,&user.cvi);
96: VecDuplicate(user.wv,&user.DPsiv);
97: VecDuplicate(user.wv,&user.DPsii);
98: VecDuplicate(user.wv,&user.DPsieta);
99: VecDuplicate(user.wv,&user.Pv);
100: VecDuplicate(user.wv,&user.Pi);
101: VecDuplicate(user.wv,&user.Piv);
102: VecDuplicate(user.wv,&user.Riv);
103: VecDuplicate(user.wv,&user.logcv);
104: VecDuplicate(user.wv,&user.logci);
105: VecDuplicate(user.wv,&user.logcvi);
106: VecDuplicate(user.wv,&user.work1);
107: VecDuplicate(user.wv,&user.work2);
108: VecDuplicate(user.wv,&user.work3);
109: VecDuplicate(user.wv,&user.work4);
111: /* Get Jacobian matrix structure from the da for the entire thing, da1 */
112: DMCreateMatrix(user.da1,MATAIJ,&user.M);
113: /* Get the (usual) mass matrix structure from da2 */
114: DMCreateMatrix(user.da2,MATAIJ,&user.M_0);
115: SetInitialGuess(x,&user);
116: /* Form the jacobian matrix and M_0 */
117: SetUpMatrices(&user);
118: MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
119:
120: SNESCreate(PETSC_COMM_WORLD,&snes);
121: SNESSetDM(snes,user.da1);
123: SNESSetFunction(snes,r,FormFunction,(void*)&user);
124: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
126: SetVariableBounds(user.da1,xl,xu);
127: SNESVISetVariableBounds(snes,xl,xu);
128: SNESSetFromOptions(snes);
129: //SNESVISetRedundancyCheck(snes,(PetscErrorCode (*)(SNES,IS,IS*,void*))CheckRedundancy,user.da1_clone);
130: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
131: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_p",FILE_MODE_WRITE,&view_p);
132: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
133: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);
134: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);
135: /* PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds); */
137: if (user.graphics) {
138: PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds);
140: VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
141: }
142: while(t<user.T) {
143: char filename[PETSC_MAX_PATH_LEN];
144: PetscScalar a = 1.0;
145: PetscInt i;
146: PetscViewer view;
147:
149: SNESSetFunction(snes,r,FormFunction,(void*)&user);
150: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
152: SetRandomVectors(&user);
153: DPsi(&user);
154: VecView(user.DPsiv,view_psi);
155: VecView(user.DPsii,view_psi);
156: VecView(user.DPsieta,view_psi);
158: VecView(user.Pv,view_p);
159: Update_q(&user);
160: VecView(user.q,view_q);
161: MatView(user.M,view_mat);
162: SNESSolve(snes,PETSC_NULL,x);
163:
164: VecView(x,view_out);
165: if (user.graphics) {
166: VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
167: }
168:
169: PetscInt its;
170: SNESGetIterationNumber(snes,&its);
171: PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4g in %d iterations\n",t,its);
173: Update_u(x,&user);
174: for (i=0; i < (int)(user.T/a) ; i++) {
175: if (t/a > i - user.dt/a && t/a < i + user.dt/a) {
176: sprintf(filename,"output_%f",t);
177: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&view);
178: VecView(user.cv,view);
179: VecView(user.ci,view);
180: VecView(user.eta,view);
181: PetscViewerDestroy(&view);
182: }
183:
184: }
185: UpdateMatrices(&user);
186: t = t + user.dt;
187:
188: }
190:
191: PetscViewerDestroy(&view_out);
192: PetscViewerDestroy(&view_p);
193: PetscViewerDestroy(&view_q);
194: PetscViewerDestroy(&view_psi);
195: PetscViewerDestroy(&view_mat);
196: VecDestroy(&x);
197: VecDestroy(&r);
198: VecDestroy(&xl);
199: VecDestroy(&xu);
200: VecDestroy(&user.q);
201: VecDestroy(&user.wv);
202: VecDestroy(&user.cv);
203: VecDestroy(&user.wi);
204: VecDestroy(&user.ci);
205: VecDestroy(&user.eta);
206: VecDestroy(&user.cvi);
207: VecDestroy(&user.DPsiv);
208: VecDestroy(&user.DPsii);
209: VecDestroy(&user.DPsieta);
210: VecDestroy(&user.Pv);
211: VecDestroy(&user.Pi);
212: VecDestroy(&user.Piv);
213: VecDestroy(&user.Riv);
214: VecDestroy(&user.logcv);
215: VecDestroy(&user.logci);
216: VecDestroy(&user.logcvi);
217: VecDestroy(&user.work1);
218: VecDestroy(&user.work2);
219: VecDestroy(&user.work3);
220: VecDestroy(&user.work4);
221: MatDestroy(&user.M);
222: MatDestroy(&user.M_0);
223: DMDestroy(&user.da1);
224: DMDestroy(&user.da1_clone);
225: DMDestroy(&user.da2);
226: SNESDestroy(&snes);
227: PetscFinalize();
228: return 0;
229: }
233: PetscErrorCode Update_u(Vec X,AppCtx *user)
234: {
236: PetscInt i,n;
237: PetscScalar *xx,*wv_p,*cv_p,*wi_p,*ci_p,*eta_p;
238: PetscInt pv1,pi1,peta1,pv2,pi2,peta2;
239: PetscReal maxv,maxi,maxeta,minv,mini,mineta;
241:
243: VecGetLocalSize(user->wv,&n);
245:
246: VecMax(user->cv,&pv1,&maxv);
247: VecMax(user->ci,&pi1,&maxi);
248: VecMax(user->eta,&peta1,&maxeta);
249: VecMin(user->cv,&pv2,&minv);
250: VecMin(user->ci,&pi2,&mini);
251: VecMin(user->eta,&peta2,&mineta);
252:
253:
254: VecGetArray(X,&xx);
255: VecGetArray(user->wv,&wv_p);
256: VecGetArray(user->cv,&cv_p);
257: VecGetArray(user->wi,&wi_p);
258: VecGetArray(user->ci,&ci_p);
259: VecGetArray(user->eta,&eta_p);
260:
261:
262: for(i=0;i<n;i++) {
263: wv_p[i] = xx[5*i];
264: cv_p[i] = xx[5*i+1];
265: wi_p[i] = xx[5*i+2];
266: ci_p[i] = xx[5*i+3];
267: eta_p[i] = xx[5*i+4];
268: }
269: VecRestoreArray(X,&xx);
270: VecRestoreArray(user->wv,&wv_p);
271: VecRestoreArray(user->cv,&cv_p);
272: VecRestoreArray(user->wi,&wi_p);
273: VecRestoreArray(user->ci,&ci_p);
274: VecRestoreArray(user->eta,&eta_p);
275:
276: return(0);
277: }
281: PetscErrorCode Update_q(AppCtx *user)
282: {
284: PetscScalar *q_p, *w1, *w2;
285: PetscInt i,n;
286: PetscScalar norm1;
289:
290: VecPointwiseMult(user->Riv,user->eta,user->eta);
291: VecScale(user->Riv,user->Rsurf);
292: VecShift(user->Riv,user->Rbulk);
293: VecPointwiseMult(user->Riv,user->ci,user->Riv);
294: VecPointwiseMult(user->Riv,user->cv,user->Riv);
295:
296: VecCopy(user->Riv,user->work1);
297: VecScale(user->work1,user->dt);
298: VecAXPY(user->work1,-1.0,user->cv);
299: MatMult(user->M_0,user->work1,user->work2);
300:
301: VecGetArray(user->q,&q_p);
302: VecGetArray(user->work1,&w1);
303: VecGetArray(user->work2,&w2);
304: VecGetLocalSize(user->wv,&n);
305: for (i=0;i<n;i++) {
306: q_p[5*i]=w2[i];
307: }
309: MatMult(user->M_0,user->DPsiv,user->work1);
310: for (i=0;i<n;i++) {
311: q_p[5*i+1]=w1[i];
312: }
314: VecCopy(user->Riv,user->work1);
315: VecScale(user->work1,user->dt);
316: VecAXPY(user->work1,-1.0,user->ci);
317: MatMult(user->M_0,user->work1,user->work2);
318: for (i=0;i<n;i++) {
319: q_p[5*i+2]=w2[i];
320: }
322: MatMult(user->M_0,user->DPsii,user->work1);
323: for (i=0;i<n;i++) {
324: q_p[5*i+3]=w1[i];
325: }
327: VecCopy(user->DPsieta,user->work1);
328: VecScale(user->work1,user->L);
329: VecAXPY(user->work1,-1.0/user->dt,user->eta);
330: MatMult(user->M_0,user->work1,user->work2);
331: for (i=0;i<n;i++) {
332: q_p[5*i+4]=w2[i];
333: }
334:
335: VecRestoreArray(user->q,&q_p);
336: VecRestoreArray(user->work1,&w1);
337: VecRestoreArray(user->work2,&w2);
340:
341: return(0);
342: }
346: PetscErrorCode DPsi(AppCtx* user)
347: {
348: PetscErrorCode ierr;
349: PetscScalar Evf=user->Evf,Eif=user->Eif,kBT=user->kBT,A=user->A;
350: PetscScalar *cv_p,*ci_p,*eta_p,*logcv_p,*logci_p,*logcvi_p,*DPsiv_p,*DPsii_p,*DPsieta_p;
351: PetscInt n,i;
355: VecGetLocalSize(user->cv,&n);
356: VecGetArray(user->cv,&cv_p);
357: VecGetArray(user->ci,&ci_p);
358: VecGetArray(user->eta,&eta_p);
359: VecGetArray(user->logcv,&logcv_p);
360: VecGetArray(user->logci,&logci_p);
361: VecGetArray(user->logcvi,&logcvi_p);
362: VecGetArray(user->DPsiv,&DPsiv_p);
363: VecGetArray(user->DPsii,&DPsii_p);
364: VecGetArray(user->DPsieta,&DPsieta_p);
366: Llog(user->cv,user->logcv);
367: Llog(user->ci,user->logci);
369: VecCopy(user->cv,user->cvi);
370: VecAXPY(user->cvi,1.0,user->ci);
371: VecScale(user->cvi,-1.0);
372: VecShift(user->cvi,1.0);
373: Llog(user->cvi,user->logcvi);
375: for (i=0;i<n;i++)
376: {
377: DPsiv_p[i] = (eta_p[i]-1.0)*(eta_p[i]-1.0)*( Evf + kBT*(logcv_p[i] - logcvi_p[i])) + eta_p[i]*eta_p[i]*2*A*(cv_p[i]-1);
379: DPsii_p[i] = (eta_p[i]-1.0)*(eta_p[i]-1.0)*( Eif + kBT*(logci_p[i] - logcvi_p[i])) + eta_p[i]*eta_p[i]*2*A*ci_p[i] ;
380:
381: DPsieta_p[i] = 2.0*(eta_p[i]-1.0)*( Evf*cv_p[i] + Eif*ci_p[i] + kBT*( cv_p[i]* logcv_p[i] + ci_p[i]* logci_p[i] + (1-cv_p[i]-ci_p[i])*logcvi_p[i] ) ) + 2.0*eta_p[i]*A*( (cv_p[i]-1.0)*(cv_p[i]-1.0) + ci_p[i]*ci_p[i]);
383:
384: }
386: VecRestoreArray(user->cv,&cv_p);
387: VecRestoreArray(user->ci,&ci_p);
388: VecRestoreArray(user->eta,&eta_p);
389: VecGetArray(user->logcv,&logcv_p);
390: VecGetArray(user->logci,&logci_p);
391: VecGetArray(user->logcvi,&logcvi_p);
392: VecRestoreArray(user->DPsiv,&DPsiv_p);
393: VecRestoreArray(user->DPsii,&DPsii_p);
394: VecRestoreArray(user->DPsieta,&DPsieta_p);
397: return(0);
399: }
404: PetscErrorCode Llog(Vec X, Vec Y)
405: {
406: PetscErrorCode ierr;
407: PetscScalar *x,*y;
408: PetscInt n,i;
411:
412: VecGetArray(X,&x);
413: VecGetArray(Y,&y);
414: VecGetLocalSize(X,&n);
415: for (i=0;i<n;i++) {
416: if (x[i] < 1.0e-12) {
417: y[i] = log(1.0e-12);
418: }
419: else {
420: y[i] = log(x[i]);
421: }
422: }
423:
424: return(0);
425: }
429: PetscErrorCode SetInitialGuess(Vec X,AppCtx* user)
430: {
431: PetscErrorCode ierr;
434: PetscInt n,i,Mda;
435: PetscScalar *xx,*cv_p,*ci_p,*wv_p,*wi_p,*eta_p;
436: PetscViewer view_out;
437: /* needed for the void growth case */
438: PetscScalar xmid,cv_v=1.0,cv_m=0.122,ci_v=0.0,ci_m=.00069,eta_v=1.0,eta_m=0.0,h,lambda;
439: PetscInt nele,nen,idx[2];
440: const PetscInt *ele;
441: PetscScalar x[2];
442: Vec coords;
443: const PetscScalar *_coords;
444: PetscViewer view;
445: PetscScalar xwidth = user->xmax - user->xmin;
449: VecGetLocalSize(X,&n);
450:
451: if (user->voidgrowth) {
452: DMDAGetInfo(user->da2,PETSC_NULL,&Mda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
453: DMDAGetGhostedCoordinates(user->da2,&coords);
454: VecGetArrayRead(coords,&_coords);
456: h = (user->xmax-user->xmin)/Mda;
457: xmid = (user->xmax + user->xmin)/2.0;
458: lambda = 4.0*h;
460: DMDAGetElements(user->da2,&nele,&nen,&ele);
461: for (i=0;i < nele; i++) {
462: idx[0] = ele[2*i]; idx[1] = ele[2*i+1];
463:
464: x[0] = _coords[idx[0]];
465: x[1] = _coords[idx[1]];
467:
468: PetscInt k;
469: PetscScalar vals_cv[2],vals_ci[2],vals_eta[2],s,hhr,r;
470: for (k=0; k < 2 ; k++) {
471: s = PetscAbs(x[k] - xmid);
472: if (s < xwidth*(5.0/64.0)) {
473: vals_cv[k] = cv_v;
474: vals_ci[k] = ci_v;
475: vals_eta[k] = eta_v;
476: } else if (s>= xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0) ) {
477: //r = (s - xwidth*(6.0/64.0) )/(0.5*lambda);
478: r = (s - xwidth*(6.0/64.0) )/(xwidth/64.0);
479: hhr = 0.25*(-r*r*r + 3*r + 2);
480: vals_cv[k] = cv_m + (1.0 - hhr)*(cv_v - cv_m);
481: vals_ci[k] = ci_m + (1.0 - hhr)*(ci_v - ci_m);
482: vals_eta[k] = eta_m + (1.0 - hhr)*(eta_v - eta_m);
483: } else {
484: vals_cv[k] = cv_m;
485: vals_ci[k] = ci_m;
486: vals_eta[k] = eta_m;
487: }
488: }
490: VecSetValuesLocal(user->cv,2,idx,vals_cv,INSERT_VALUES);
491: VecSetValuesLocal(user->ci,2,idx,vals_ci,INSERT_VALUES);
492: VecSetValuesLocal(user->eta,2,idx,vals_eta,INSERT_VALUES);
493: }
494: VecAssemblyBegin(user->cv);
495: VecAssemblyEnd(user->cv);
496: VecAssemblyBegin(user->ci);
497: VecAssemblyEnd(user->ci);
498: VecAssemblyBegin(user->eta);
499: VecAssemblyEnd(user->eta);
501: DMDARestoreElements(user->da2,&nele,&nen,&ele);
502: VecRestoreArrayRead(coords,&_coords);
504:
505: } else {
507: VecGetArray(user->cv,&cv_p);
508: VecGetArray(user->ci,&ci_p);
510:
511: /* VecSet(user->cv,0.122);
512: VecSet(user->cv,6.9e-8);
513: VecSet(user->ci,6.9e-8); */
514:
515:
516: for (i=0;i<n/5;i++)
517: {
518: if (i%5 <4 ) {
519: cv_p[i] = 0.0;
520: ci_p[i] = 1.0e-2;
521: }
522: else {
523: cv_p[i] = 1.0e-2;
524: ci_p[i] = 0.0;
525: }
526: }
527: VecRestoreArray(user->cv,&cv_p);
528: VecRestoreArray(user->ci,&ci_p);
531: VecSet(user->eta,0.0);
533: }
535: DPsi(user);
536: VecCopy(user->DPsiv,user->wv);
537: VecCopy(user->DPsii,user->wi);
539: VecGetArray(X,&xx);
540: VecGetArray(user->cv,&cv_p);
541: VecGetArray(user->ci,&ci_p);
542: VecGetArray(user->eta,&eta_p);
543: VecGetArray(user->wv,&wv_p);
544: VecGetArray(user->wi,&wi_p);
545: for (i=0;i<n/5;i++)
546: {
547: xx[5*i]=wv_p[i];
548: xx[5*i+1]=cv_p[i];
549: xx[5*i+2]=wi_p[i];
550: xx[5*i+3]=ci_p[i];
551: xx[5*i+4]=eta_p[i];
553:
554:
555: }
557: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);
558: VecView(user->wv,view_out);
559: VecView(user->cv,view_out);
560: VecView(user->wi,view_out);
561: VecView(user->ci,view_out);
562: VecView(user->eta,view_out);
563: PetscViewerDestroy(&view_out);
565: VecRestoreArray(X,&xx);
566: VecRestoreArray(user->cv,&cv_p);
567: VecRestoreArray(user->ci,&ci_p);
568: VecRestoreArray(user->wv,&wv_p);
569: VecRestoreArray(user->wi,&wi_p); VecRestoreArray(user->eta,&eta_p);
570: return(0);
571: }
575: PetscErrorCode SetRandomVectors(AppCtx* user)
576: {
578: PetscInt i,n,count=0;
579: PetscScalar *w1,*w2,*Pv_p,*eta_p;
581: /* static PetscViewer viewer=0; */
582: static PetscRandom rand = 0;
583: static PetscInt step = 0;
586: if (!rand) {
587: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
588: PetscRandomSetFromOptions(rand);
589: }
590:
591: VecSetRandom(user->work1,rand);
592: VecSetRandom(user->work2,rand);
593: VecGetArray(user->work1,&w1);
594: VecGetArray(user->work2,&w2);
595: VecGetArray(user->Pv,&Pv_p);
596: VecGetArray(user->eta,&eta_p);
597: VecGetLocalSize(user->work1,&n);
598: for (i=0;i<n;i++) {
599: if (eta_p[i]>=0.8 || w1[i]>user->P_casc){
600: Pv_p[i]=0;
601: }
602: else
603: {
604: Pv_p[i]=w2[i]*user->VG;
605: count = count + 1;
606: }
607: }
608: step++;
610: VecCopy(user->Pv,user->Pi);
611: VecScale(user->Pi,0.9);
612: VecPointwiseMult(user->Piv,user->Pi,user->Pv);
613: VecRestoreArray(user->work1,&w1);
614: VecRestoreArray(user->work2,&w2);
615: VecRestoreArray(user->Pv,&Pv_p);
616: VecRestoreArray(user->eta,&eta_p);
617: printf("Number of radiations count %d out of total mesh points n %d\n",count,n);
619: return(0);
620:
621: }
624: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ctx)
625: {
627: AppCtx *user=(AppCtx*)ctx;
628:
630: MatMultAdd(user->M,X,user->q,F);
631: return(0);
632: }
636: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
637: {
639: AppCtx *user=(AppCtx*)ctx;
640:
642: *flg = SAME_NONZERO_PATTERN;
643: MatCopy(user->M,*J,*flg);
644: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
645: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
646: return(0);
647: }
650: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
651: {
653: PetscScalar **l,**u;
654: PetscInt xs,xm;
655: PetscInt i;
656:
658: DMDAVecGetArrayDOF(da,xl,&l);
659: DMDAVecGetArrayDOF(da,xu,&u);
660:
661: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
662:
663:
664: for(i=xs; i < xs+xm;i++) {
665: l[i][0] = -SNES_VI_INF;
666: l[i][1] = 0.0;
667: l[i][2] = -SNES_VI_INF;
668: l[i][3] = 0.0;
669: l[i][4] = 0.0;
670: u[i][0] = SNES_VI_INF;
671: u[i][1] = 1.0;
672: u[i][2] = SNES_VI_INF;
673: u[i][3] = 1.0;
674: u[i][4] = 1.0;
675: }
676:
677:
678: DMDAVecRestoreArrayDOF(da,xl,&l);
679: DMDAVecRestoreArrayDOF(da,xu,&u);
680: return(0);
681: }
685: PetscErrorCode GetParams(AppCtx* user)
686: {
688: PetscBool flg;
689:
691:
692: /* Set default parameters */
693: user->xmin = 0.0; user->xmax = 128.0;
694: user->Dv = 1.0; user->Di=1.0;
695: user->Evf = 0.8; user->Eif = 0.8;
696: user->A = 1.0;
697: user->kBT = 0.11;
698: user->kav = 1.0; user->kai = 1.0; user->kaeta = 1.0;
699: user->Rsurf = 10.0; user->Rbulk = 0.0;
700: user->L = 10.0; user->P_casc = 0.05;
701: user->T = 1.0e-2; user->dt = 1.0e-4;
702: user->VG = 100.0;
703: user->initv = .0001;
704: user->initeta = 0.0;
705: user->graphics = PETSC_TRUE;
706:
707: user->degenerate = PETSC_FALSE;
708: /* void growth */
709: user->voidgrowth = PETSC_FALSE;
711: PetscOptionsGetReal(PETSC_NULL,"-xmin",&user->xmin,&flg);
712: PetscOptionsGetReal(PETSC_NULL,"-xmax",&user->xmax,&flg);
713: PetscOptionsGetReal(PETSC_NULL,"-T",&user->T,&flg);
714: PetscOptionsGetReal(PETSC_NULL,"-dt",&user->dt,&flg);
715: PetscOptionsBool("-degenerate","Run with degenerate mobility\n","None",user->degenerate,&user->degenerate,&flg);
716: PetscOptionsReal("-smallnumber","Small number added to degenerate mobility\n","None",user->smallnumber,&user->smallnumber,&flg);
717: PetscOptionsBool("-voidgrowth","Use initial conditions for void growth\n","None",user->voidgrowth,&user->voidgrowth,&flg);
718: PetscOptionsBool("-graphics","Contour plot solutions at each timestep\n","None",user->graphics,&user->graphics,&flg);
719:
720: return(0);
721: }
726: PetscErrorCode SetUpMatrices(AppCtx* user)
727: {
728: PetscErrorCode ierr;
729: PetscInt nele,nen,i,n;
730: const PetscInt *ele;
731: PetscScalar dt=user->dt,h;
732:
733: PetscInt idx[2];
734: PetscScalar eM_0[2][2],eM_2[2][2];
735: PetscScalar cv_sum, ci_sum;
736: Mat M=user->M;
737: Mat M_0=user->M_0;
738: PetscInt Mda;
739: PetscScalar *cv_p,*ci_p;
740: /* newly added */
741: Vec cvlocal,cilocal;
742:
745: /* new stuff */
746: DMGetLocalVector(user->da2,&cvlocal);
747: DMGetLocalVector(user->da2,&cilocal);
748: DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
749: DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
750: DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
751: DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
753: /* old stuff */
754: /*
755: VecGetArray(user->cv,&cv_p);
756: VecGetArray(user->ci,&ci_p);
757: */
758: /* new stuff */
759: VecGetArray(cvlocal,&cv_p);
760: VecGetArray(cilocal,&ci_p);
762: MatGetLocalSize(M,&n,PETSC_NULL);
763: DMDAGetInfo(user->da1,PETSC_NULL,&Mda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
765: h = (user->xmax-user->xmin)/Mda;
766:
767: eM_0[0][0]=eM_0[1][1]=h/3.0;
768: eM_0[0][1]=eM_0[1][0]=h/6.0;
770: eM_2[0][0]=eM_2[1][1]=1.0/h;
771: eM_2[0][1]=eM_2[1][0]=-1.0/h;
773: /* Get local element info */
774: DMDAGetElements(user->da1,&nele,&nen,&ele);
775: //for(i=0;i < nele + 1;i++) {
776: for(i=0;i < nele;i++) {
778: idx[0] = ele[2*i]; idx[1] = ele[2*i+1];
780: if (user->degenerate) {
781: cv_sum = (2.0*user->smallnumber + cv_p[idx[0]]+cv_p[idx[1]])*user->Dv/(2.0*user->kBT);
782: ci_sum = (2.0*user->smallnumber + ci_p[idx[0]]+ci_p[idx[1]])*user->Di/(2.0*user->kBT);
783: } else {
784: cv_sum = user->initv*user->Dv/user->kBT;
785: ci_sum = user->initv*user->Di/user->kBT;
786: }
789: PetscInt row,cols[4],r,row_M_0,cols2[2];
790: //PetscScalar vals[4],vals_M_0[1],vals2[2];
791: PetscScalar vals[4],vals_M_0[2],vals2[2];
792:
793: for(r=0;r<2;r++) {
794: row_M_0 = idx[r];
795:
796: vals_M_0[0]=eM_0[r][0];
797: vals_M_0[1]=eM_0[r][1];
798:
799: //vals_M_0[0] = h;
800: MatSetValuesLocal(M_0,1,&row_M_0,2,idx,vals_M_0,ADD_VALUES);
801: //MatSetValuesLocal(M_0,1,&row_M_0,1,&row_M_0,vals_M_0,INSERT_VALUES);
802:
803: row = 5*idx[r];
804: cols[0] = 5*idx[0]; vals[0] = dt*eM_2[r][0]*cv_sum;
805: cols[1] = 5*idx[1]; vals[1] = dt*eM_2[r][1]*cv_sum;
806: cols[2] = 5*idx[0]+1; vals[2] = eM_0[r][0];
807: cols[3] = 5*idx[1]+1; vals[3] = eM_0[r][1];
808:
809: /* Insert values in matrix M for 1st dof */
810: MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);
811:
812: row = 5*idx[r]+1;
813: cols[0] = 5*idx[0]; vals[0] = -eM_0[r][0];
814: cols[1] = 5*idx[1]; vals[1] = -eM_0[r][1];
815: cols[2] = 5*idx[0]+1; vals[2] = user->kav*eM_2[r][0];
816: cols[3] = 5*idx[1]+1; vals[3] = user->kav*eM_2[r][1];
818: /* Insert values in matrix M for 2nd dof */
819: MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);
820:
821: row = 5*idx[r]+2;
822: cols[0] = 5*idx[0]+2; vals[0] = dt*eM_2[r][0]*ci_sum;
823: cols[1] = 5*idx[1]+2; vals[1] = dt*eM_2[r][1]*ci_sum;
824: cols[2] = 5*idx[0]+3; vals[2] = eM_0[r][0];
825: cols[3] = 5*idx[1]+3; vals[3] = eM_0[r][1];
826: /* Insert values in matrix M for 3nd dof */
827: MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);
828:
829: row = 5*idx[r]+3;
830: cols[0] = 5*idx[0]+2; vals[0] = -eM_0[r][0];
831: cols[1] = 5*idx[1]+2; vals[1] = -eM_0[r][1];
832: cols[2] = 5*idx[0]+3; vals[2] = user->kai*eM_2[r][0];
833: cols[3] = 5*idx[1]+3; vals[3] = user->kai*eM_2[r][1];
834: /* Insert values in matrix M for 4th dof */
835: MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);
837: row = 5*idx[r]+4;
838: cols2[0] = 5*idx[0]+4; vals2[0] = eM_0[r][0]/dt + user->L*user->kaeta*eM_2[r][0];
839: cols2[1] = 5*idx[1]+4; vals2[1] = eM_0[r][1]/dt + user->L*user->kaeta*eM_2[r][1];
840: MatSetValuesLocal(M,1,&row,2,cols2,vals2,ADD_VALUES);
841: }
842: }
844:
845:
846: /* new */
847: VecRestoreArray(cvlocal,&cv_p);
848: VecRestoreArray(cilocal,&ci_p);
849: DMRestoreLocalVector(user->da2,&cvlocal);
850: DMRestoreLocalVector(user->da2,&cilocal);
851: /*
852: VecRestoreArray(user->cv,&cv_p);
853: VecRestoreArray(user->ci,&ci_p);*/
854: MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
855: MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);
856:
857: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
858: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
859:
860: DMDARestoreElements(user->da1,&nele,&nen,&ele);
861:
862:
863: return(0);
864: }
869: PetscErrorCode UpdateMatrices(AppCtx* user)
870: {
871: PetscErrorCode ierr;
872: PetscInt nele,nen,i,n,Mda;
873: const PetscInt *ele;
874:
875: PetscInt idx[2];
876: PetscScalar eM_2[2][2],h;
877: Mat M=user->M;
878: PetscScalar *cv_p,*ci_p,cv_sum,ci_sum;
879: /* newly added */
880: Vec cvlocal,cilocal;
883:
884: /*new stuff */
885: DMGetLocalVector(user->da2,&cvlocal);
886: DMGetLocalVector(user->da2,&cilocal);
887: DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
888: DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
889: DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
890: DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
891: /* old stuff */
892: /*
893: VecGetArray(user->cv,&cv_p);
894: VecGetArray(user->ci,&ci_p);
895: */
897: /* new stuff */
898: VecGetArray(cvlocal,&cv_p);
899: VecGetArray(cilocal,&ci_p);
901: /* Create the mass matrix M_0 */
902: MatGetLocalSize(M,&n,PETSC_NULL);
903: DMDAGetInfo(user->da1,PETSC_NULL,&Mda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
905: h = (user->xmax-user->xmin)/Mda;
907: /* Get local element info */
908: DMDAGetElements(user->da1,&nele,&nen,&ele);
909:
910: for(i=0;i<nele;i++) {
911: idx[0] = ele[2*i]; idx[1] = ele[2*i+1];
912:
913: PetscInt r,row,cols[2];
914: PetscScalar vals[2];
916: for (r=0;r<2;r++) {
917: row = 5*idx[r];
918: cols[0] = 5*idx[0]; vals[0] = 0.0;
919: cols[1] = 5*idx[1]; vals[1] = 0.0;
920:
922: /* Insert values in matrix M for 1st dof */
923: MatSetValuesLocal(M,1,&row,2,cols,vals,INSERT_VALUES);
924:
925: row = 5*idx[r]+2;
926: cols[0] = 5*idx[0]+2; vals[0] = 0.0;
927: cols[1] = 5*idx[1]+2; vals[1] = 0.0;
929: /* Insert values in matrix M for 3nd dof */
930: MatSetValuesLocal(M,1,&row,2,cols,vals,INSERT_VALUES);
931: }
932: }
933:
934: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
935: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
937: eM_2[0][0]=eM_2[1][1]=1.0/h;
938: eM_2[0][1]=eM_2[1][0]=-1.0/h;
940: for(i=0;i<nele;i++) {
941: idx[0] = ele[2*i]; idx[1] = ele[2*i+1];
942:
943: PetscInt row,cols[2],r;
944: PetscScalar vals[2];
946: for(r=0;r<2;r++) {
948: if (user->degenerate) {
949: cv_sum = (2.0*user->smallnumber + cv_p[idx[0]] + cv_p[idx[1]])*user->Dv/(2.0*user->kBT);
950: ci_sum = (2.0*user->smallnumber + ci_p[idx[0]] + ci_p[idx[1]])*user->Di/(2.0*user->kBT);
951: } else {
952: cv_sum = user->initv*user->Dv/user->kBT;
953: ci_sum = user->initv*user->Di/user->kBT;
954: }
956: row = 5*idx[r];
957: cols[0] = 5*idx[0]; vals[0] = user->dt*eM_2[r][0]*cv_sum;
958: cols[1] = 5*idx[1]; vals[1] = user->dt*eM_2[r][1]*cv_sum;
960: /* Insert values in matrix M for 1st dof */
961: MatSetValuesLocal(M,1,&row,2,cols,vals,ADD_VALUES);
963: row = 5*idx[r]+2;
964: cols[0] = 5*idx[0]+2; vals[0] = user->dt*eM_2[r][0]*ci_sum;
965: cols[1] = 5*idx[1]+2; vals[1] = user->dt*eM_2[r][1]*ci_sum;
967: /* Insert values in matrix M for 3nd dof */
968: MatSetValuesLocal(M,1,&row,2,cols,vals,ADD_VALUES);
969: }
970: }
971:
972: /* old stuff
973: VecRestoreArray(user->cv,&cv_p);
974: VecRestoreArray(user->ci,&ci_p);
975: */
977: /* new stuff */
978: VecRestoreArray(cvlocal,&cv_p);
979: VecRestoreArray(cilocal,&ci_p);
980: DMRestoreLocalVector(user->da2,&cvlocal);
981: DMRestoreLocalVector(user->da2,&cilocal);
982:
983: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
984: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
985: DMDARestoreElements(user->da1,&nele,&nen,&ele);
987: return(0);
988: }
993: PetscErrorCode CheckRedundancy(SNES snes, IS act, IS *outact, DM da)
994: {
996: PetscScalar **uin,**uout;
997: Vec UIN, UOUT;
998: PetscInt xs,xm,*outindex;
999: const PetscInt *index;
1000: PetscInt k,i,l,n,M,cnt=0;
1001:
1003: DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1004: DMGetGlobalVector(da,&UIN);
1005: VecSet(UIN,0.0);
1006: DMGetLocalVector(da,&UOUT);
1007: DMDAVecGetArrayDOF(da,UIN,&uin);
1008: ISGetIndices(act,&index);
1009: ISGetLocalSize(act,&n);
1010: for (k=0;k<n;k++){
1011: l = index[k]%5;
1012: i = index[k]/5;
1013: uin[i][l]=1.0;
1014: }
1015: printf("Number of active constraints before applying redundancy %d\n",n);
1016: ISRestoreIndices(act,&index);
1017: DMDAVecRestoreArrayDOF(da,UIN,&uin);
1018: DMGlobalToLocalBegin(da,UIN,INSERT_VALUES,UOUT);
1019: DMGlobalToLocalEnd(da,UIN,INSERT_VALUES,UOUT);
1020: DMDAVecGetArrayDOF(da,UOUT,&uout);
1021:
1022: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
1024: for(i=xs; i < xs+xm;i++) {
1025: if (uout[i-1][1] && uout[i][1] && uout[i+1][1])
1026: uout[i][0] = 1.0;
1027: if (uout[i-1][3] && uout[i][3] && uout[i+1][3])
1028: uout[i][2] = 1.0;
1029: }
1031: for(i=xs; i < xs+xm;i++) {
1032: for(l=0;l<5;l++) {
1033: if (uout[i][l])
1034: cnt++;
1035: }
1036: }
1038: printf("Number of active constraints after applying redundancy %d\n",cnt);
1039:
1041: PetscMalloc(cnt*sizeof(PetscInt),&outindex);
1042: cnt = 0;
1043:
1044: for(i=xs; i < xs+xm;i++) {
1045: for(l=0;l<5;l++) {
1046: if (uout[i][l])
1047: outindex[cnt++] = 5*(i)+l;
1048: }
1049: }
1050:
1051:
1052: ISCreateGeneral(PETSC_COMM_WORLD,cnt,outindex,PETSC_OWN_POINTER,outact);
1053: DMDAVecRestoreArrayDOF(da,UOUT,&uout);
1054: DMRestoreGlobalVector(da,&UIN);
1055: DMRestoreLocalVector(da,&UOUT);
1056: return(0);
1057: }