Actual source code: ex63.c
petsc-3.5.4 2015-05-23
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,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; /* Solution 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);
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,DM_BOUNDARY_PERIODIC, -4, 5, 1,NULL,&user.da1);
71: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC, -4, 5, 1,NULL,&user.da1_clone);
72: /* Create a 1D DA with dof = 1; for individual componentes */
73: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC, -4, 1, 1,NULL,&user.da2);
75: /* Set Element type (triangular) */
76: DMDASetElementType(user.da1,DMDA_ELEMENT_P1);
77: DMDASetElementType(user.da2,DMDA_ELEMENT_P1);
79: /* Set x and y coordinates */
80: DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,NULL,NULL,NULL,NULL);
81: DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,NULL,NULL,NULL,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);
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: DMSetMatType(user.da1,MATAIJ);
113: DMCreateMatrix(user.da1,&user.M);
114: /* Get the (usual) mass matrix structure from da2 */
115: DMSetMatType(user.da2,MATAIJ);
116: DMCreateMatrix(user.da2,&user.M_0);
117: SetInitialGuess(x,&user);
118: /* Form the jacobian matrix and M_0 */
119: SetUpMatrices(&user);
120: MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
122: SNESCreate(PETSC_COMM_WORLD,&snes);
123: SNESSetDM(snes,user.da1);
125: SNESSetFunction(snes,r,FormFunction,(void*)&user);
126: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
128: SetVariableBounds(user.da1,xl,xu);
129: SNESVISetVariableBounds(snes,xl,xu);
130: SNESSetFromOptions(snes);
131: /* SNESVISetRedundancyCheck(snes,(PetscErrorCode (*)(SNES,IS,IS*,void*))CheckRedundancy,user.da1_clone); */
132: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
133: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_p",FILE_MODE_WRITE,&view_p);
134: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
135: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);
136: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);
137: /* PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds); */
139: if (user.graphics) {
140: PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds);
142: VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
143: }
144: while (t<user.T) {
145: char filename[PETSC_MAX_PATH_LEN];
146: PetscScalar a = 1.0;
147: PetscInt i;
148: PetscViewer view;
151: SNESSetFunction(snes,r,FormFunction,(void*)&user);
152: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
154: SetRandomVectors(&user);
155: DPsi(&user);
156: VecView(user.DPsiv,view_psi);
157: VecView(user.DPsii,view_psi);
158: VecView(user.DPsieta,view_psi);
160: VecView(user.Pv,view_p);
161: Update_q(&user);
162: VecView(user.q,view_q);
163: MatView(user.M,view_mat);
164: SNESSolve(snes,NULL,x);
166: VecView(x,view_out);
167: if (user.graphics) {
168: VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
169: }
171: PetscInt its;
172: SNESGetIterationNumber(snes,&its);
173: PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4g in %d iterations\n",t,its);
175: Update_u(x,&user);
176: for (i=0; i < (int)(user.T/a) ; i++) {
177: if (t/a > i - user.dt/a && t/a < i + user.dt/a) {
178: sprintf(filename,"output_%f",t);
179: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&view);
180: VecView(user.cv,view);
181: VecView(user.ci,view);
182: VecView(user.eta,view);
183: PetscViewerDestroy(&view);
184: }
185: }
186: UpdateMatrices(&user);
187: t = t + user.dt;
189: }
192: PetscViewerDestroy(&view_out);
193: PetscViewerDestroy(&view_p);
194: PetscViewerDestroy(&view_q);
195: PetscViewerDestroy(&view_psi);
196: PetscViewerDestroy(&view_mat);
197: VecDestroy(&x);
198: VecDestroy(&r);
199: VecDestroy(&xl);
200: VecDestroy(&xu);
201: VecDestroy(&user.q);
202: VecDestroy(&user.wv);
203: VecDestroy(&user.cv);
204: VecDestroy(&user.wi);
205: VecDestroy(&user.ci);
206: VecDestroy(&user.eta);
207: VecDestroy(&user.cvi);
208: VecDestroy(&user.DPsiv);
209: VecDestroy(&user.DPsii);
210: VecDestroy(&user.DPsieta);
211: VecDestroy(&user.Pv);
212: VecDestroy(&user.Pi);
213: VecDestroy(&user.Piv);
214: VecDestroy(&user.Riv);
215: VecDestroy(&user.logcv);
216: VecDestroy(&user.logci);
217: VecDestroy(&user.logcvi);
218: VecDestroy(&user.work1);
219: VecDestroy(&user.work2);
220: VecDestroy(&user.work3);
221: VecDestroy(&user.work4);
222: MatDestroy(&user.M);
223: MatDestroy(&user.M_0);
224: DMDestroy(&user.da1);
225: DMDestroy(&user.da1_clone);
226: DMDestroy(&user.da2);
227: SNESDestroy(&snes);
228: PetscFinalize();
229: return 0;
230: }
234: PetscErrorCode Update_u(Vec X,AppCtx *user)
235: {
237: PetscInt i,n;
238: PetscScalar *xx,*wv_p,*cv_p,*wi_p,*ci_p,*eta_p;
239: PetscInt pv1,pi1,peta1,pv2,pi2,peta2;
240: PetscReal maxv,maxi,maxeta,minv,mini,mineta;
244: VecGetLocalSize(user->wv,&n);
247: VecMax(user->cv,&pv1,&maxv);
248: VecMax(user->ci,&pi1,&maxi);
249: VecMax(user->eta,&peta1,&maxeta);
250: VecMin(user->cv,&pv2,&minv);
251: VecMin(user->ci,&pi2,&mini);
252: VecMin(user->eta,&peta2,&mineta);
255: VecGetArray(X,&xx);
256: VecGetArray(user->wv,&wv_p);
257: VecGetArray(user->cv,&cv_p);
258: VecGetArray(user->wi,&wi_p);
259: VecGetArray(user->ci,&ci_p);
260: VecGetArray(user->eta,&eta_p);
263: for (i=0; i<n; i++) {
264: wv_p[i] = xx[5*i];
265: cv_p[i] = xx[5*i+1];
266: wi_p[i] = xx[5*i+2];
267: ci_p[i] = xx[5*i+3];
268: eta_p[i] = xx[5*i+4];
269: }
270: VecRestoreArray(X,&xx);
271: VecRestoreArray(user->wv,&wv_p);
272: VecRestoreArray(user->cv,&cv_p);
273: VecRestoreArray(user->wi,&wi_p);
274: VecRestoreArray(user->ci,&ci_p);
275: VecRestoreArray(user->eta,&eta_p);
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: VecPointwiseMult(user->Riv,user->eta,user->eta);
290: VecScale(user->Riv,user->Rsurf);
291: VecShift(user->Riv,user->Rbulk);
292: VecPointwiseMult(user->Riv,user->ci,user->Riv);
293: VecPointwiseMult(user->Riv,user->cv,user->Riv);
295: VecCopy(user->Riv,user->work1);
296: VecScale(user->work1,user->dt);
297: VecAXPY(user->work1,-1.0,user->cv);
298: MatMult(user->M_0,user->work1,user->work2);
300: VecGetArray(user->q,&q_p);
301: VecGetArray(user->work1,&w1);
302: VecGetArray(user->work2,&w2);
303: VecGetLocalSize(user->wv,&n);
304: for (i=0; i<n; i++) q_p[5*i]=w2[i];
306: MatMult(user->M_0,user->DPsiv,user->work1);
307: for (i=0; i<n; i++) q_p[5*i+1]=w1[i];
309: VecCopy(user->Riv,user->work1);
310: VecScale(user->work1,user->dt);
311: VecAXPY(user->work1,-1.0,user->ci);
312: MatMult(user->M_0,user->work1,user->work2);
313: for (i=0; i<n; i++) q_p[5*i+2]=w2[i];
315: MatMult(user->M_0,user->DPsii,user->work1);
316: for (i=0; i<n; i++) q_p[5*i+3]=w1[i];
318: VecCopy(user->DPsieta,user->work1);
319: VecScale(user->work1,user->L);
320: VecAXPY(user->work1,-1.0/user->dt,user->eta);
321: MatMult(user->M_0,user->work1,user->work2);
322: for (i=0; i<n; i++) q_p[5*i+4]=w2[i];
324: VecRestoreArray(user->q,&q_p);
325: VecRestoreArray(user->work1,&w1);
326: VecRestoreArray(user->work2,&w2);
327: return(0);
328: }
332: PetscErrorCode DPsi(AppCtx *user)
333: {
335: PetscScalar Evf=user->Evf,Eif=user->Eif,kBT=user->kBT,A=user->A;
336: PetscScalar *cv_p,*ci_p,*eta_p,*logcv_p,*logci_p,*logcvi_p,*DPsiv_p,*DPsii_p,*DPsieta_p;
337: PetscInt n,i;
340: VecGetLocalSize(user->cv,&n);
341: VecGetArray(user->cv,&cv_p);
342: VecGetArray(user->ci,&ci_p);
343: VecGetArray(user->eta,&eta_p);
344: VecGetArray(user->logcv,&logcv_p);
345: VecGetArray(user->logci,&logci_p);
346: VecGetArray(user->logcvi,&logcvi_p);
347: VecGetArray(user->DPsiv,&DPsiv_p);
348: VecGetArray(user->DPsii,&DPsii_p);
349: VecGetArray(user->DPsieta,&DPsieta_p);
351: Llog(user->cv,user->logcv);
352: Llog(user->ci,user->logci);
354: VecCopy(user->cv,user->cvi);
355: VecAXPY(user->cvi,1.0,user->ci);
356: VecScale(user->cvi,-1.0);
357: VecShift(user->cvi,1.0);
358: Llog(user->cvi,user->logcvi);
360: for (i=0; i<n; i++)
361: {
362: 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);
363: 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];
365: 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]);
366: }
368: VecRestoreArray(user->cv,&cv_p);
369: VecRestoreArray(user->ci,&ci_p);
370: VecRestoreArray(user->eta,&eta_p);
371: VecGetArray(user->logcv,&logcv_p);
372: VecGetArray(user->logci,&logci_p);
373: VecGetArray(user->logcvi,&logcvi_p);
374: VecRestoreArray(user->DPsiv,&DPsiv_p);
375: VecRestoreArray(user->DPsii,&DPsii_p);
376: VecRestoreArray(user->DPsieta,&DPsieta_p);
377: return(0);
379: }
384: PetscErrorCode Llog(Vec X, Vec Y)
385: {
387: PetscScalar *x,*y;
388: PetscInt n,i;
391: VecGetArray(X,&x);
392: VecGetArray(Y,&y);
393: VecGetLocalSize(X,&n);
394: for (i=0; i<n; i++) {
395: if (x[i] < 1.0e-12) y[i] = log(1.0e-12);
396: else y[i] = log(x[i]);
397: }
398: return(0);
399: }
403: PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
404: {
407: PetscInt n,i,Mda;
408: PetscScalar *xx,*cv_p,*ci_p,*wv_p,*wi_p,*eta_p;
409: PetscViewer view_out;
410: /* needed for the void growth case */
411: 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;
412: PetscInt nele,nen,idx[2];
413: const PetscInt *ele;
414: PetscScalar x[2];
415: Vec coords;
416: const PetscScalar *_coords;
417: PetscViewer view;
418: PetscScalar xwidth = user->xmax - user->xmin;
421: VecGetLocalSize(X,&n);
423: if (user->voidgrowth) {
424: DMDAGetInfo(user->da2,NULL,&Mda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
425: DMGetCoordinatesLocal(user->da2,&coords);
426: VecGetArrayRead(coords,&_coords);
428: h = (user->xmax-user->xmin)/Mda;
429: xmid = (user->xmax + user->xmin)/2.0;
430: lambda = 4.0*h;
432: DMDAGetElements(user->da2,&nele,&nen,&ele);
433: for (i=0; i < nele; i++) {
434: idx[0] = ele[2*i]; idx[1] = ele[2*i+1];
436: x[0] = _coords[idx[0]];
437: x[1] = _coords[idx[1]];
440: PetscInt k;
441: PetscScalar vals_cv[2],vals_ci[2],vals_eta[2],s,hhr,r;
442: for (k=0; k < 2; k++) {
443: s = PetscAbs(x[k] - xmid);
444: if (s < xwidth*(5.0/64.0)) {
445: vals_cv[k] = cv_v;
446: vals_ci[k] = ci_v;
447: vals_eta[k] = eta_v;
448: } else if (s>= xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0)) {
449: /* r = (s - xwidth*(6.0/64.0))/(0.5*lambda); */
450: r = (s - xwidth*(6.0/64.0))/(xwidth/64.0);
451: hhr = 0.25*(-r*r*r + 3*r + 2);
452: vals_cv[k] = cv_m + (1.0 - hhr)*(cv_v - cv_m);
453: vals_ci[k] = ci_m + (1.0 - hhr)*(ci_v - ci_m);
454: vals_eta[k] = eta_m + (1.0 - hhr)*(eta_v - eta_m);
455: } else {
456: vals_cv[k] = cv_m;
457: vals_ci[k] = ci_m;
458: vals_eta[k] = eta_m;
459: }
460: }
462: VecSetValuesLocal(user->cv,2,idx,vals_cv,INSERT_VALUES);
463: VecSetValuesLocal(user->ci,2,idx,vals_ci,INSERT_VALUES);
464: VecSetValuesLocal(user->eta,2,idx,vals_eta,INSERT_VALUES);
465: }
466: VecAssemblyBegin(user->cv);
467: VecAssemblyEnd(user->cv);
468: VecAssemblyBegin(user->ci);
469: VecAssemblyEnd(user->ci);
470: VecAssemblyBegin(user->eta);
471: VecAssemblyEnd(user->eta);
473: DMDARestoreElements(user->da2,&nele,&nen,&ele);
474: VecRestoreArrayRead(coords,&_coords);
477: } else {
479: VecGetArray(user->cv,&cv_p);
480: VecGetArray(user->ci,&ci_p);
483: /* VecSet(user->cv,0.122);
484: VecSet(user->cv,6.9e-8);
485: VecSet(user->ci,6.9e-8); */
488: for (i=0; i<n/5; i++) {
489: if (i%5 <4) {
490: cv_p[i] = 0.0;
491: ci_p[i] = 1.0e-2;
492: } else {
493: cv_p[i] = 1.0e-2;
494: ci_p[i] = 0.0;
495: }
496: }
497: VecRestoreArray(user->cv,&cv_p);
498: VecRestoreArray(user->ci,&ci_p);
501: VecSet(user->eta,0.0);
503: }
505: DPsi(user);
506: VecCopy(user->DPsiv,user->wv);
507: VecCopy(user->DPsii,user->wi);
509: VecGetArray(X,&xx);
510: VecGetArray(user->cv,&cv_p);
511: VecGetArray(user->ci,&ci_p);
512: VecGetArray(user->eta,&eta_p);
513: VecGetArray(user->wv,&wv_p);
514: VecGetArray(user->wi,&wi_p);
515: for (i=0; i<n/5; i++)
516: {
517: xx[5*i] = wv_p[i];
518: xx[5*i+1] = cv_p[i];
519: xx[5*i+2] = wi_p[i];
520: xx[5*i+3] = ci_p[i];
521: xx[5*i+4] = eta_p[i];
525: }
527: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);
528: VecView(user->wv,view_out);
529: VecView(user->cv,view_out);
530: VecView(user->wi,view_out);
531: VecView(user->ci,view_out);
532: VecView(user->eta,view_out);
533: PetscViewerDestroy(&view_out);
535: VecRestoreArray(X,&xx);
536: VecRestoreArray(user->cv,&cv_p);
537: VecRestoreArray(user->ci,&ci_p);
538: VecRestoreArray(user->wv,&wv_p);
539: VecRestoreArray(user->wi,&wi_p); VecRestoreArray(user->eta,&eta_p);
540: return(0);
541: }
545: PetscErrorCode SetRandomVectors(AppCtx *user)
546: {
548: PetscInt i,n,count=0;
549: PetscScalar *w1,*w2,*Pv_p,*eta_p;
551: /* static PetscViewer viewer=0; */
552: static PetscRandom rand = 0;
553: static PetscInt step = 0;
556: if (!rand) {
557: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
558: PetscRandomSetFromOptions(rand);
559: }
561: VecSetRandom(user->work1,rand);
562: VecSetRandom(user->work2,rand);
563: VecGetArray(user->work1,&w1);
564: VecGetArray(user->work2,&w2);
565: VecGetArray(user->Pv,&Pv_p);
566: VecGetArray(user->eta,&eta_p);
567: VecGetLocalSize(user->work1,&n);
568: for (i=0; i<n; i++) {
569: if (eta_p[i]>=0.8 || w1[i]>user->P_casc) Pv_p[i]=0;
570: else {
571: Pv_p[i]=w2[i]*user->VG;
572: count = count + 1;
573: }
574: }
575: step++;
577: VecCopy(user->Pv,user->Pi);
578: VecScale(user->Pi,0.9);
579: VecPointwiseMult(user->Piv,user->Pi,user->Pv);
580: VecRestoreArray(user->work1,&w1);
581: VecRestoreArray(user->work2,&w2);
582: VecRestoreArray(user->Pv,&Pv_p);
583: VecRestoreArray(user->eta,&eta_p);
584: printf("Number of radiations count %d out of total mesh points n %d\n",count,n);
585: return(0);
587: }
590: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
591: {
593: AppCtx *user=(AppCtx*)ctx;
596: MatMultAdd(user->M,X,user->q,F);
597: return(0);
598: }
602: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat B,void *ctx)
603: {
605: AppCtx *user=(AppCtx*)ctx;
608: MatCopy(user->M,J,*flg);
609: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
610: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
611: return(0);
612: }
615: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
616: {
618: PetscScalar **l,**u;
619: PetscInt xs,xm;
620: PetscInt i;
623: DMDAVecGetArrayDOF(da,xl,&l);
624: DMDAVecGetArrayDOF(da,xu,&u);
626: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
629: for (i=xs; i < xs+xm; i++) {
630: l[i][0] = -PETSC_INFINITY;
631: l[i][1] = 0.0;
632: l[i][2] = -PETSC_INFINITY;
633: l[i][3] = 0.0;
634: l[i][4] = 0.0;
635: u[i][0] = PETSC_INFINITY;
636: u[i][1] = 1.0;
637: u[i][2] = PETSC_INFINITY;
638: u[i][3] = 1.0;
639: u[i][4] = 1.0;
640: }
643: DMDAVecRestoreArrayDOF(da,xl,&l);
644: DMDAVecRestoreArrayDOF(da,xu,&u);
645: return(0);
646: }
650: PetscErrorCode GetParams(AppCtx *user)
651: {
653: PetscBool flg;
656: /* Set default parameters */
657: user->xmin = 0.0; user->xmax = 128.0;
658: user->Dv = 1.0; user->Di = 1.0;
659: user->Evf = 0.8; user->Eif = 0.8;
660: user->A = 1.0;
661: user->kBT = 0.11;
662: user->kav = 1.0; user->kai = 1.0; user->kaeta = 1.0;
663: user->Rsurf = 10.0; user->Rbulk = 0.0;
664: user->L = 10.0; user->P_casc = 0.05;
665: user->T = 1.0e-2; user->dt = 1.0e-4;
666: user->VG = 100.0;
667: user->initv = .0001;
668: user->initeta = 0.0;
669: user->graphics = PETSC_TRUE;
671: user->degenerate = PETSC_FALSE;
672: /* void growth */
673: user->voidgrowth = PETSC_FALSE;
675: PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);
676: PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);
677: PetscOptionsGetReal(NULL,"-T",&user->T,&flg);
678: PetscOptionsGetReal(NULL,"-dt",&user->dt,&flg);
679: PetscOptionsBool("-degenerate","Run with degenerate mobility\n","None",user->degenerate,&user->degenerate,&flg);
680: PetscOptionsReal("-smallnumber","Small number added to degenerate mobility\n","None",user->smallnumber,&user->smallnumber,&flg);
681: PetscOptionsBool("-voidgrowth","Use initial conditions for void growth\n","None",user->voidgrowth,&user->voidgrowth,&flg);
682: PetscOptionsBool("-graphics","Contour plot solutions at each timestep\n","None",user->graphics,&user->graphics,&flg);
683: return(0);
684: }
689: PetscErrorCode SetUpMatrices(AppCtx *user)
690: {
692: PetscInt nele,nen,i,n;
693: const PetscInt *ele;
694: PetscScalar dt=user->dt,h;
696: PetscInt idx[2];
697: PetscScalar eM_0[2][2],eM_2[2][2];
698: PetscScalar cv_sum, ci_sum;
699: Mat M =user->M;
700: Mat M_0=user->M_0;
701: PetscInt Mda;
702: PetscScalar *cv_p,*ci_p;
703: /* newly added */
704: Vec cvlocal,cilocal;
707: /* new stuff */
708: DMGetLocalVector(user->da2,&cvlocal);
709: DMGetLocalVector(user->da2,&cilocal);
710: DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
711: DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
712: DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
713: DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
715: /* old stuff */
716: /*
717: VecGetArray(user->cv,&cv_p);
718: VecGetArray(user->ci,&ci_p);
719: */
720: /* new stuff */
721: VecGetArray(cvlocal,&cv_p);
722: VecGetArray(cilocal,&ci_p);
724: MatGetLocalSize(M,&n,NULL);
725: DMDAGetInfo(user->da1,NULL,&Mda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
727: h = (user->xmax-user->xmin)/Mda;
729: eM_0[0][0]=eM_0[1][1]=h/3.0;
730: eM_0[0][1]=eM_0[1][0]=h/6.0;
732: eM_2[0][0]=eM_2[1][1]=1.0/h;
733: eM_2[0][1]=eM_2[1][0]=-1.0/h;
735: /* Get local element info */
736: DMDAGetElements(user->da1,&nele,&nen,&ele);
737: /* for (i=0;i < nele + 1;i++) { */
738: for (i=0; i < nele; i++) {
740: idx[0] = ele[2*i]; idx[1] = ele[2*i+1];
742: if (user->degenerate) {
743: cv_sum = (2.0*user->smallnumber + cv_p[idx[0]]+cv_p[idx[1]])*user->Dv/(2.0*user->kBT);
744: ci_sum = (2.0*user->smallnumber + ci_p[idx[0]]+ci_p[idx[1]])*user->Di/(2.0*user->kBT);
745: } else {
746: cv_sum = user->initv*user->Dv/user->kBT;
747: ci_sum = user->initv*user->Di/user->kBT;
748: }
751: PetscInt row,cols[4],r,row_M_0,cols2[2];
752: /* PetscScalar vals[4],vals_M_0[1],vals2[2]; */
753: PetscScalar vals[4],vals_M_0[2],vals2[2];
755: for (r=0; r<2; r++) {
756: row_M_0 = idx[r];
758: vals_M_0[0]=eM_0[r][0];
759: vals_M_0[1]=eM_0[r][1];
761: /* vals_M_0[0] = h; */
762: MatSetValuesLocal(M_0,1,&row_M_0,2,idx,vals_M_0,ADD_VALUES);
763: /* MatSetValuesLocal(M_0,1,&row_M_0,1,&row_M_0,vals_M_0,INSERT_VALUES); */
765: row = 5*idx[r];
766: cols[0] = 5*idx[0]; vals[0] = dt*eM_2[r][0]*cv_sum;
767: cols[1] = 5*idx[1]; vals[1] = dt*eM_2[r][1]*cv_sum;
768: cols[2] = 5*idx[0]+1; vals[2] = eM_0[r][0];
769: cols[3] = 5*idx[1]+1; vals[3] = eM_0[r][1];
771: /* Insert values in matrix M for 1st dof */
772: MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);
774: row = 5*idx[r]+1;
775: cols[0] = 5*idx[0]; vals[0] = -eM_0[r][0];
776: cols[1] = 5*idx[1]; vals[1] = -eM_0[r][1];
777: cols[2] = 5*idx[0]+1; vals[2] = user->kav*eM_2[r][0];
778: cols[3] = 5*idx[1]+1; vals[3] = user->kav*eM_2[r][1];
780: /* Insert values in matrix M for 2nd dof */
781: MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);
783: row = 5*idx[r]+2;
784: cols[0] = 5*idx[0]+2; vals[0] = dt*eM_2[r][0]*ci_sum;
785: cols[1] = 5*idx[1]+2; vals[1] = dt*eM_2[r][1]*ci_sum;
786: cols[2] = 5*idx[0]+3; vals[2] = eM_0[r][0];
787: cols[3] = 5*idx[1]+3; vals[3] = eM_0[r][1];
788: /* Insert values in matrix M for 3nd dof */
789: MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);
791: row = 5*idx[r]+3;
792: cols[0] = 5*idx[0]+2; vals[0] = -eM_0[r][0];
793: cols[1] = 5*idx[1]+2; vals[1] = -eM_0[r][1];
794: cols[2] = 5*idx[0]+3; vals[2] = user->kai*eM_2[r][0];
795: cols[3] = 5*idx[1]+3; vals[3] = user->kai*eM_2[r][1];
796: /* Insert values in matrix M for 4th dof */
797: MatSetValuesLocal(M,1,&row,4,cols,vals,ADD_VALUES);
799: row = 5*idx[r]+4;
800: cols2[0] = 5*idx[0]+4; vals2[0] = eM_0[r][0]/dt + user->L*user->kaeta*eM_2[r][0];
801: cols2[1] = 5*idx[1]+4; vals2[1] = eM_0[r][1]/dt + user->L*user->kaeta*eM_2[r][1];
802: MatSetValuesLocal(M,1,&row,2,cols2,vals2,ADD_VALUES);
803: }
804: }
808: /* new */
809: VecRestoreArray(cvlocal,&cv_p);
810: VecRestoreArray(cilocal,&ci_p);
811: DMRestoreLocalVector(user->da2,&cvlocal);
812: DMRestoreLocalVector(user->da2,&cilocal);
813: /*
814: VecRestoreArray(user->cv,&cv_p);
815: VecRestoreArray(user->ci,&ci_p);*/
816: MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
817: MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);
819: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
820: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
822: DMDARestoreElements(user->da1,&nele,&nen,&ele);
823: return(0);
824: }
829: PetscErrorCode UpdateMatrices(AppCtx *user)
830: {
832: PetscInt nele,nen,i,n,Mda;
833: const PetscInt *ele;
835: PetscInt idx[2];
836: PetscScalar eM_2[2][2],h;
837: Mat M=user->M;
838: PetscScalar *cv_p,*ci_p,cv_sum,ci_sum;
839: /* newly added */
840: Vec cvlocal,cilocal;
843: /*new stuff */
844: DMGetLocalVector(user->da2,&cvlocal);
845: DMGetLocalVector(user->da2,&cilocal);
846: DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
847: DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
848: DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
849: DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
850: /* old stuff */
851: /*
852: VecGetArray(user->cv,&cv_p);
853: VecGetArray(user->ci,&ci_p);
854: */
856: /* new stuff */
857: VecGetArray(cvlocal,&cv_p);
858: VecGetArray(cilocal,&ci_p);
860: /* Create the mass matrix M_0 */
861: MatGetLocalSize(M,&n,NULL);
862: DMDAGetInfo(user->da1,NULL,&Mda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
864: h = (user->xmax-user->xmin)/Mda;
866: /* Get local element info */
867: DMDAGetElements(user->da1,&nele,&nen,&ele);
869: for (i=0; i<nele; i++) {
870: idx[0] = ele[2*i]; idx[1] = ele[2*i+1];
872: PetscInt r,row,cols[2];
873: PetscScalar vals[2];
875: for (r=0; r<2; r++) {
876: row = 5*idx[r];
877: cols[0] = 5*idx[0]; vals[0] = 0.0;
878: cols[1] = 5*idx[1]; vals[1] = 0.0;
881: /* Insert values in matrix M for 1st dof */
882: MatSetValuesLocal(M,1,&row,2,cols,vals,INSERT_VALUES);
884: row = 5*idx[r]+2;
885: cols[0] = 5*idx[0]+2; vals[0] = 0.0;
886: cols[1] = 5*idx[1]+2; vals[1] = 0.0;
888: /* Insert values in matrix M for 3nd dof */
889: MatSetValuesLocal(M,1,&row,2,cols,vals,INSERT_VALUES);
890: }
891: }
893: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
894: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
896: eM_2[0][0]=eM_2[1][1]=1.0/h;
897: eM_2[0][1]=eM_2[1][0]=-1.0/h;
899: for (i=0; i<nele; i++) {
900: idx[0] = ele[2*i]; idx[1] = ele[2*i+1];
902: PetscInt row,cols[2],r;
903: PetscScalar vals[2];
905: for (r=0; r<2; r++) {
907: if (user->degenerate) {
908: cv_sum = (2.0*user->smallnumber + cv_p[idx[0]] + cv_p[idx[1]])*user->Dv/(2.0*user->kBT);
909: ci_sum = (2.0*user->smallnumber + ci_p[idx[0]] + ci_p[idx[1]])*user->Di/(2.0*user->kBT);
910: } else {
911: cv_sum = user->initv*user->Dv/user->kBT;
912: ci_sum = user->initv*user->Di/user->kBT;
913: }
915: row = 5*idx[r];
916: cols[0] = 5*idx[0]; vals[0] = user->dt*eM_2[r][0]*cv_sum;
917: cols[1] = 5*idx[1]; vals[1] = user->dt*eM_2[r][1]*cv_sum;
919: /* Insert values in matrix M for 1st dof */
920: MatSetValuesLocal(M,1,&row,2,cols,vals,ADD_VALUES);
922: row = 5*idx[r]+2;
923: cols[0] = 5*idx[0]+2; vals[0] = user->dt*eM_2[r][0]*ci_sum;
924: cols[1] = 5*idx[1]+2; vals[1] = user->dt*eM_2[r][1]*ci_sum;
926: /* Insert values in matrix M for 3nd dof */
927: MatSetValuesLocal(M,1,&row,2,cols,vals,ADD_VALUES);
928: }
929: }
931: /* old stuff
932: VecRestoreArray(user->cv,&cv_p);
933: VecRestoreArray(user->ci,&ci_p);
934: */
936: /* new stuff */
937: VecRestoreArray(cvlocal,&cv_p);
938: VecRestoreArray(cilocal,&ci_p);
939: DMRestoreLocalVector(user->da2,&cvlocal);
940: DMRestoreLocalVector(user->da2,&cilocal);
942: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
943: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
944: DMDARestoreElements(user->da1,&nele,&nen,&ele);
945: return(0);
946: }
951: PetscErrorCode CheckRedundancy(SNES snes, IS act, IS *outact, DM da)
952: {
954: PetscScalar **uin,**uout;
955: Vec UIN, UOUT;
956: PetscInt xs,xm,*outindex;
957: const PetscInt *index;
958: PetscInt k,i,l,n,M,cnt=0;
961: DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
962: DMGetGlobalVector(da,&UIN);
963: VecSet(UIN,0.0);
964: DMGetLocalVector(da,&UOUT);
965: DMDAVecGetArrayDOF(da,UIN,&uin);
966: ISGetIndices(act,&index);
967: ISGetLocalSize(act,&n);
968: for (k=0; k<n; k++) {
969: l = index[k]%5;
970: i = index[k]/5;
971: uin[i][l]=1.0;
972: }
973: printf("Number of active constraints before applying redundancy %d\n",n);
974: ISRestoreIndices(act,&index);
975: DMDAVecRestoreArrayDOF(da,UIN,&uin);
976: DMGlobalToLocalBegin(da,UIN,INSERT_VALUES,UOUT);
977: DMGlobalToLocalEnd(da,UIN,INSERT_VALUES,UOUT);
978: DMDAVecGetArrayDOF(da,UOUT,&uout);
980: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
982: for (i=xs; i < xs+xm; i++) {
983: if (uout[i-1][1] && uout[i][1] && uout[i+1][1]) uout[i][0] = 1.0;
984: if (uout[i-1][3] && uout[i][3] && uout[i+1][3]) uout[i][2] = 1.0;
985: }
987: for (i=xs; i < xs+xm; i++) {
988: for (l=0; l < 5; l++) {
989: if (uout[i][l]) cnt++;
990: }
991: }
993: printf("Number of active constraints after applying redundancy %d\n",cnt);
996: PetscMalloc1(cnt,&outindex);
997: cnt = 0;
999: for (i=xs; i < xs+xm; i++) {
1000: for (l=0; l < 5; l++) {
1001: if (uout[i][l]) outindex[cnt++] = 5*(i)+l;
1002: }
1003: }
1006: ISCreateGeneral(PETSC_COMM_WORLD,cnt,outindex,PETSC_OWN_POINTER,outact);
1007: DMDAVecRestoreArrayDOF(da,UOUT,&uout);
1008: DMRestoreGlobalVector(da,&UIN);
1009: DMRestoreLocalVector(da,&UOUT);
1010: return(0);
1011: }