Actual source code: ex633d_db.c
petsc-3.5.4 2015-05-23
1: static char help[] = "3D 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: ./ex633d_db -ksp_type fgmres -snes_vi_monitor -snes_atol 1.e-11 -snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor -VG 10000000 -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_type basic -T .0020 -P_casc .0005 -snes_monitor_solution -da_refine 1
13: ./ex633d_db -ksp_type fgmres -snes_vi_monitor -snes_atol 1.e-11 -snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor -VG 10000000 -pc_type mg -pc_mg_galerkin -log_summary -da_refine 1
15: */
17: #include <petscsnes.h>
18: #include <petscdmda.h>
20: typedef struct {
21: PetscReal dt,T; /* Time step and end time */
22: DM da1,da1_clone,da2;
23: Mat M; /* Jacobian matrix */
24: Mat M_0;
25: Vec q,wv,cv,wi,ci,eta,cvi,DPsiv,DPsii,DPsieta,Pv,Pi,Piv,logcv,logci,logcvi,Riv;
26: Vec work1,work2,work3,work4;
27: PetscScalar Dv,Di,Evf,Eif,A,kBT,kav,kai,kaeta,Rsurf,Rbulk,L,P_casc,VG; /* physics parameters */
28: PetscReal xmin,xmax,ymin,ymax,zmin,zmax;
29: PetscInt nx;
30: PetscBool voidgrowth;
31: PetscBool degenerate;
32: PetscBool periodic;
33: PetscReal smallnumber;
34: PetscReal initv;
35: PetscReal initeta;
36: } AppCtx;
38: PetscErrorCode GetParams(AppCtx*);
39: PetscErrorCode SetRandomVectors(AppCtx*);
40: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
41: PetscErrorCode SetUpMatrices(AppCtx*);
42: PetscErrorCode UpdateMatrices(AppCtx*);
43: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
44: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
45: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
46: PetscErrorCode Update_q(AppCtx*);
47: PetscErrorCode Update_u(Vec,AppCtx*);
48: PetscErrorCode DPsi(AppCtx*);
49: PetscErrorCode Llog(Vec,Vec);
50: PetscErrorCode CheckRedundancy(SNES,IS,IS*,DM);
54: int main(int argc, char **argv)
55: {
57: Vec x,r; /* olution and residual vectors */
58: SNES snes; /* Nonlinear solver context */
59: AppCtx user; /* Application context */
60: Vec xl,xu; /* Upper and lower bounds on variables */
61: Mat J;
62: PetscScalar t=0.0;
63: PetscViewer view_out, view_p, view_q, view_psi, view_mat;
64: PetscReal bounds[] = {1000.0,-1000.,0.0,1.0,1000.0,-1000.0,0.0,1.0,1000.0,-1000.0};
67: PetscInitialize(&argc,&argv, (char*)0, help);
70: /* Get physics and time parameters */
71: GetParams(&user);
73: if (user.periodic) {
74: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX, -3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, 5, 1,NULL,NULL,NULL,&user.da1);
75: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX, -3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, 5, 1,NULL,NULL,NULL,&user.da1_clone);
76: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX, -3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, 1, 1,NULL,NULL,NULL,&user.da2);
78: } else {
79: /* Create a 1D DA with dof = 5; the whole thing */
80: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX, -3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, 5, 1,NULL,NULL,NULL,&user.da1);
81: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX, -3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, 5, 1,NULL,NULL,NULL,&user.da1_clone);
82: /* Create a 1D DA with dof = 1; for individual componentes */
83: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX, -3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, 1, 1,NULL,NULL,NULL,&user.da2);
84: }
86: /* Set Element type (rectangular) */
87: DMDASetElementType(user.da1,DMDA_ELEMENT_Q1);
88: DMDASetElementType(user.da2,DMDA_ELEMENT_Q1);
90: /* Set x and y coordinates */
91: DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,user.ymin,user.ymax, user.zmin,user.zmax);
92: DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,user.ymin,user.ymax, user.zmin,user.zmax);
94: /* Get global vector x from DM (da1) and duplicate vectors r,xl,xu */
95: DMCreateGlobalVector(user.da1,&x);
96: VecDuplicate(x,&r);
97: VecDuplicate(x,&xl);
98: VecDuplicate(x,&xu);
99: VecDuplicate(x,&user.q);
102: /* Get global vector user->wv from da2 and duplicate other vectors */
103: DMCreateGlobalVector(user.da2,&user.wv);
104: VecDuplicate(user.wv,&user.cv);
105: VecDuplicate(user.wv,&user.wi);
106: VecDuplicate(user.wv,&user.ci);
107: VecDuplicate(user.wv,&user.eta);
108: VecDuplicate(user.wv,&user.cvi);
109: VecDuplicate(user.wv,&user.DPsiv);
110: VecDuplicate(user.wv,&user.DPsii);
111: VecDuplicate(user.wv,&user.DPsieta);
112: VecDuplicate(user.wv,&user.Pv);
113: VecDuplicate(user.wv,&user.Pi);
114: VecDuplicate(user.wv,&user.Piv);
115: VecDuplicate(user.wv,&user.Riv);
116: VecDuplicate(user.wv,&user.logcv);
117: VecDuplicate(user.wv,&user.logci);
118: VecDuplicate(user.wv,&user.logcvi);
119: VecDuplicate(user.wv,&user.work1);
120: VecDuplicate(user.wv,&user.work2);
121: VecDuplicate(user.wv,&user.work3);
122: VecDuplicate(user.wv,&user.work4);
124: /* Get Jacobian matrix structure from the da for the entire thing, da1 */
125: DMSetMatType(user.da1,MATAIJ);
126: DMCreateMatrix(user.da1,&user.M);
127: /* Get the (usual) mass matrix structure from da2 */
128: DMSetMatType(user.da2,MATAIJ);
129: DMCreateMatrix(user.da2,&user.M_0);
131: SetInitialGuess(x,&user);
133: /* Form the jacobian matrix and M_0 */
134: SetUpMatrices(&user);
137: MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
139: SNESCreate(PETSC_COMM_WORLD,&snes);
140: SNESSetDM(snes,user.da1);
142: SNESSetFunction(snes,r,FormFunction,(void*)&user);
143: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
145: SetVariableBounds(user.da1,xl,xu);
146: /* SNESVISetRedundancyCheck(snes,(PetscErrorCode (*)(SNES,IS,IS*,void*))CheckRedundancy,user.da1_clone); */
147: SNESVISetVariableBounds(snes,xl,xu);
148: SNESSetFromOptions(snes);
150: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
151: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_p",FILE_MODE_WRITE,&view_p);
152: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
153: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);
154: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);
157: /*
158: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
159: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);
161: VecView(user.q,view_q);
162: MatView(user.M,view_mat);
164: PetscViewerDestroy(&view_q);
165: PetscViewerDestroy(&view_mat);
166: */
170: while (t<user.T) {
171: char filename[PETSC_MAX_PATH_LEN];
172: PetscScalar a = 1.0;
173: PetscInt i;
174: PetscViewer view;
177: SNESSetFunction(snes,r,FormFunction,(void*)&user);
178: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
180: SetRandomVectors(&user);
181: DPsi(&user);
182: VecView(user.DPsiv,view_psi);
183: VecView(user.DPsii,view_psi);
184: VecView(user.DPsieta,view_psi);
186: VecView(user.Pv,view_p);
187: Update_q(&user);
188: VecView(user.q,view_q);
190: printf("after VecView\n");
191: MatView(user.M,view_mat);
193: printf("after MatView\n");
195: SNESSolve(snes,NULL,x);
197: VecView(x,view_out);
199: PetscInt its;
200: SNESGetIterationNumber(snes,&its);
201: PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4g in %d iterations\n",t,its);
203: Update_u(x,&user);
204: /*
205: for (i=0; i < (int)(user.T/a) ; i++) {
206: if (t/a > i - user.dt/a && t/a < i + user.dt/a) {
207: sprintf(filename,"output_%f",t);
208: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&view);
209: VecView(user.cv,view);
210: VecView(user.ci,view);
211: VecView(user.eta,view);
212: PetscViewerDestroy(&view);
213: }
215: }
216: */
217: UpdateMatrices(&user);
218: t = t + user.dt;
220: }
222: PetscViewerDestroy(&view_out);
223: PetscViewerDestroy(&view_p);
224: PetscViewerDestroy(&view_q);
225: PetscViewerDestroy(&view_psi);
226: PetscViewerDestroy(&view_mat);
227: VecDestroy(&x);
228: VecDestroy(&r);
229: VecDestroy(&xl);
230: VecDestroy(&xu);
231: VecDestroy(&user.q);
232: VecDestroy(&user.wv);
233: VecDestroy(&user.cv);
234: VecDestroy(&user.wi);
235: VecDestroy(&user.ci);
236: VecDestroy(&user.eta);
237: VecDestroy(&user.cvi);
238: VecDestroy(&user.DPsiv);
239: VecDestroy(&user.DPsii);
240: VecDestroy(&user.DPsieta);
241: VecDestroy(&user.Pv);
242: VecDestroy(&user.Pi);
243: VecDestroy(&user.Piv);
244: VecDestroy(&user.Riv);
245: VecDestroy(&user.logcv);
246: VecDestroy(&user.logci);
247: VecDestroy(&user.logcvi);
248: VecDestroy(&user.work1);
249: VecDestroy(&user.work2);
250: VecDestroy(&user.work3);
251: VecDestroy(&user.work4);
252: MatDestroy(&user.M);
253: MatDestroy(&user.M_0);
254: DMDestroy(&user.da1);
255: DMDestroy(&user.da1_clone);
256: DMDestroy(&user.da2);
257: SNESDestroy(&snes);
259: printf("I am finalize \n");
260: PetscFinalize();
261: return 0;
262: }
266: PetscErrorCode Update_u(Vec X,AppCtx *user)
267: {
269: PetscInt i,n;
270: PetscScalar *xx,*wv_p,*cv_p,*wi_p,*ci_p,*eta_p;
273: VecGetLocalSize(user->wv,&n);
275: VecGetArray(X,&xx);
276: VecGetArray(user->wv,&wv_p);
277: VecGetArray(user->cv,&cv_p);
278: VecGetArray(user->wi,&wi_p);
279: VecGetArray(user->ci,&ci_p);
280: VecGetArray(user->eta,&eta_p);
283: for (i=0; i<n; i++) {
284: wv_p[i] = xx[5*i];
285: cv_p[i] = xx[5*i+1];
286: wi_p[i] = xx[5*i+2];
287: ci_p[i] = xx[5*i+3];
288: eta_p[i] = xx[5*i+4];
289: }
290: VecRestoreArray(X,&xx);
291: VecRestoreArray(user->wv,&wv_p);
292: VecRestoreArray(user->cv,&cv_p);
293: VecRestoreArray(user->wi,&wi_p);
294: VecRestoreArray(user->ci,&ci_p);
295: VecRestoreArray(user->eta,&eta_p);
296: return(0);
297: }
301: PetscErrorCode Update_q(AppCtx *user)
302: {
304: PetscScalar *q_p, *w1, *w2;
305: PetscInt i,n;
308: VecPointwiseMult(user->Riv,user->eta,user->eta); /* Riv = eta.^2 */
309: VecScale(user->Riv,user->Rsurf); /* Riv = Rsurf * eta.^2 */
310: VecShift(user->Riv,user->Rbulk); /* Riv = Rbulk + Rsurf * eta.^2 */
311: VecPointwiseMult(user->Riv,user->ci,user->Riv); /* Riv = (Rbulk + Rsurf * eta.^2) * ci */
312: VecPointwiseMult(user->Riv,user->cv,user->Riv); /* Riv = (Rbulk + Rsurf * eta.^2) * ci * cv */
314: VecCopy(user->Riv,user->work1);
315: VecScale(user->work1,user->dt); /* work1 = dt * Riv */
316: VecAXPY(user->work1,-1.0,user->cv); /* work1 = dt * Riv - cv */
317: MatMult(user->M_0,user->work1,user->work2); /* work2 = M_0 * (dt * Riv - cv) */
319: VecGetArray(user->q,&q_p);
320: VecGetArray(user->work1,&w1);
321: VecGetArray(user->work2,&w2);
322: VecGetLocalSize(user->wv,&n);
323: for (i=0; i<n; i++) q_p[5*i]=w2[i];
325: MatMult(user->M_0,user->DPsiv,user->work1);
326: for (i=0; i<n; i++) q_p[5*i+1]=w1[i];
328: VecCopy(user->Riv,user->work1);
329: VecScale(user->work1,user->dt);
330: VecAXPY(user->work1,-1.0,user->ci);
331: MatMult(user->M_0,user->work1,user->work2);
332: for (i=0; i<n; i++) q_p[5*i+2]=w2[i];
334: MatMult(user->M_0,user->DPsii,user->work1);
335: for (i=0; i<n; i++) q_p[5*i+3]=w1[i];
337: VecCopy(user->DPsieta,user->work1);
338: VecScale(user->work1,user->L);
339: VecAXPY(user->work1,-1.0/user->dt,user->eta);
340: MatMult(user->M_0,user->work1,user->work2);
341: for (i=0; i<n; i++) q_p[5*i+4]=w2[i];
343: VecRestoreArray(user->q,&q_p);
344: VecRestoreArray(user->work1,&w1);
345: VecRestoreArray(user->work2,&w2);
346: return(0);
347: }
351: PetscErrorCode DPsi(AppCtx *user)
352: {
354: PetscScalar Evf=user->Evf,Eif=user->Eif,kBT=user->kBT,A=user->A;
355: PetscScalar *cv_p,*ci_p,*eta_p,*logcv_p,*logci_p,*logcvi_p,*DPsiv_p,*DPsii_p,*DPsieta_p;
356: PetscInt n,i;
359: VecGetLocalSize(user->cv,&n);
360: VecGetArray(user->cv,&cv_p);
361: VecGetArray(user->ci,&ci_p);
362: VecGetArray(user->eta,&eta_p);
363: VecGetArray(user->logcv,&logcv_p);
364: VecGetArray(user->logci,&logci_p);
365: VecGetArray(user->logcvi,&logcvi_p);
366: VecGetArray(user->DPsiv,&DPsiv_p);
367: VecGetArray(user->DPsii,&DPsii_p);
368: VecGetArray(user->DPsieta,&DPsieta_p);
370: Llog(user->cv,user->logcv);
371: Llog(user->ci,user->logci);
373: VecCopy(user->cv,user->cvi);
374: VecAXPY(user->cvi,1.0,user->ci);
375: VecScale(user->cvi,-1.0);
376: VecShift(user->cvi,1.0);
377: Llog(user->cvi,user->logcvi);
379: for (i=0; i<n; i++)
380: {
381: 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);
382: 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];
384: 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]);
385: }
387: VecRestoreArray(user->cv,&cv_p);
388: VecRestoreArray(user->ci,&ci_p);
389: VecRestoreArray(user->eta,&eta_p);
390: VecGetArray(user->logcv,&logcv_p);
391: VecGetArray(user->logci,&logci_p);
392: VecGetArray(user->logcvi,&logcvi_p);
393: VecRestoreArray(user->DPsiv,&DPsiv_p);
394: VecRestoreArray(user->DPsii,&DPsii_p);
395: VecRestoreArray(user->DPsieta,&DPsieta_p);
396: return(0);
397: }
402: PetscErrorCode Llog(Vec X, Vec Y)
403: {
405: PetscScalar *x,*y;
406: PetscInt n,i;
409: VecGetArray(X,&x);
410: VecGetArray(Y,&y);
411: VecGetLocalSize(X,&n);
412: for (i=0; i<n; i++) {
413: if (x[i] < 1.0e-12) y[i] = log(1.0e-12);
414: else y[i] = log(x[i]);
415: }
416: return(0);
417: }
421: PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
422: {
423: /******************** 3D ********************/
425: PetscInt n,i,j,Xda,Yda,Zda;
426: PetscScalar *xx,*cv_p,*ci_p,*wv_p,*wi_p,*eta_p;
427: PetscViewer view_out;
429: /* needed for the void growth case */
430: PetscScalar 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;
431: PetscScalar xmid,ymid,zmid;
432: PetscInt nele,nen,idx[8];
433: const PetscInt *ele;
434: PetscScalar x[8],y[8],z[8];
435: Vec coords;
436: const PetscScalar *_coords;
437: PetscViewer view;
438: PetscScalar xwidth = user->xmax - user->xmin, ywidth = user->ymax - user->ymin, zwidth = user->zmax - user->zmin;
441: VecGetLocalSize(X,&n);
443: if (user->voidgrowth) {
444: DMDAGetInfo(user->da2,NULL,&Xda,&Yda,&Zda,NULL,NULL,
445: NULL,NULL,NULL,NULL,NULL,NULL,NULL);
446: DMGetCoordinatesLocal(user->da2,&coords);
447: VecGetArrayRead(coords,&_coords);
449: if (user->periodic) h = (user->xmax-user->xmin)/Xda;
450: else h = (user->xmax-user->xmin)/(Xda-1.0);
452: xmid = (user->xmax + user->xmin)/2.0;
453: ymid = (user->ymax + user->ymin)/2.0;
454: zmid = (user->zmax + user->zmin)/2.0;
455: lambda = 4.0*h;
457: DMDAGetElements(user->da2,&nele,&nen,&ele); /* number of local elements, number of element nodes, the local indices of the elements' vertices */
458: for (i=0; i < nele; i++) {
459: for (j = 0; j<8; j++) {
460: idx[j] = ele[8*i + j];
461: x[j] = _coords[3*idx[j]];
462: y[j] = _coords[3*idx[j]+1];
463: z[j] = _coords[3*idx[j]+2];
464: }
466: PetscInt k;
467: PetscScalar vals_cv[8],vals_ci[8],vals_eta[8],s,hhr,r;
468: for (k=0; k < 8; k++) {
469: s = PetscSqrtScalar((x[k] - xmid)*(x[k] - xmid) + (y[k] - ymid)*(y[k] - ymid) + (z[k] - zmid)*(z[k] - zmid));
470: if (s < xwidth*(5.0/64.0)) {
471: vals_cv[k] = cv_v;
472: vals_ci[k] = ci_v;
473: vals_eta[k] = eta_v;
474: } else if (s>= xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0)) {
475: /*r = (s - xwidth*(6.0/64.0))/(0.5*lambda);*/
476: r = (s - xwidth*(6.0/64.0))/(xwidth/64.0);
477: hhr = 0.25*(-r*r*r + 3*r + 2);
478: vals_cv[k] = cv_m + (1.0 - hhr)*(cv_v - cv_m);
479: vals_ci[k] = ci_m + (1.0 - hhr)*(ci_v - ci_m);
480: vals_eta[k] = eta_m + (1.0 - hhr)*(eta_v - eta_m);
481: } else {
482: vals_cv[k] = cv_m;
483: vals_ci[k] = ci_m;
484: vals_eta[k] = eta_m;
485: }
486: }
488: VecSetValuesLocal(user->cv,8,idx,vals_cv,INSERT_VALUES);
489: VecSetValuesLocal(user->ci,8,idx,vals_ci,INSERT_VALUES);
490: VecSetValuesLocal(user->eta,8,idx,vals_eta,INSERT_VALUES);
491: }
492: VecAssemblyBegin(user->cv);
493: VecAssemblyEnd(user->cv);
494: VecAssemblyBegin(user->ci);
495: VecAssemblyEnd(user->ci);
496: VecAssemblyBegin(user->eta);
497: VecAssemblyEnd(user->eta);
499: DMDARestoreElements(user->da2,&nele,&nen,&ele);
500: VecRestoreArrayRead(coords,&_coords);
502: } else {
504: VecSet(user->cv,6.9e-4);
505: VecSet(user->ci,6.9e-4);
506: VecSet(user->eta,0.0);
507: }
509: DPsi(user);
510: VecCopy(user->DPsiv,user->wv);
511: VecCopy(user->DPsii,user->wi);
513: VecGetArray(X,&xx);
514: VecGetArray(user->cv,&cv_p);
515: VecGetArray(user->ci,&ci_p);
516: VecGetArray(user->eta,&eta_p);
517: VecGetArray(user->wv,&wv_p);
518: VecGetArray(user->wi,&wi_p);
519: for (i=0; i<n/5; i++) {
520: xx[5*i] =wv_p[i];
521: xx[5*i+1]=cv_p[i];
522: xx[5*i+2]=wi_p[i];
523: xx[5*i+3]=ci_p[i];
524: 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);
540: VecRestoreArray(user->eta,&eta_p);
541: return(0);
542: }
546: PetscErrorCode SetRandomVectors(AppCtx *user)
547: {
549: PetscInt i,n,count=0;
550: PetscScalar *w1,*w2,*Pv_p,*eta_p;
552: /* static PetscViewer viewer=0; */
553: static PetscRandom rand = 0;
554: static PetscInt step = 0;
557: if (!rand) {
558: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
559: PetscRandomSetFromOptions(rand);
560: }
562: VecSetRandom(user->work1,rand);
563: VecSetRandom(user->work2,rand);
564: VecGetArray(user->work1,&w1);
565: VecGetArray(user->work2,&w2);
566: VecGetArray(user->Pv,&Pv_p);
567: VecGetArray(user->eta,&eta_p);
568: VecGetLocalSize(user->work1,&n);
569: for (i=0; i<n; i++) {
570: if (eta_p[i]>=0.8 || w1[i]>user->P_casc) Pv_p[i]=0;
571: else {
572: Pv_p[i] = w2[i]*user->VG;
573: count = count + 1;
574: }
575: }
576: step++;
578: VecCopy(user->Pv,user->Pi);
579: VecScale(user->Pi,0.9);
580: VecPointwiseMult(user->Piv,user->Pi,user->Pv);
581: VecRestoreArray(user->work1,&w1);
582: VecRestoreArray(user->work2,&w2);
583: VecRestoreArray(user->Pv,&Pv_p);
584: VecRestoreArray(user->eta,&eta_p);
585: printf("Number of radiations count %d out of total mesh points n %d\n",count,n);
586: return(0);
587: }
591: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
592: {
594: AppCtx *user=(AppCtx*)ctx;
597: MatMultAdd(user->M,X,user->q,F);
598: return(0);
599: }
603: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat B,void *ctx)
604: {
606: AppCtx *user=(AppCtx*)ctx;
609: MatCopy(user->M,J,*flg);
610: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
611: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
612: return(0);
613: }
617: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
618: {
620: PetscScalar ****l,****u;
621: PetscInt xs,xm,ys,ym,zs,zm;
622: PetscInt k,j,i;
625: DMDAVecGetArrayDOF(da,xl,&l);
626: DMDAVecGetArrayDOF(da,xu,&u);
628: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
630: for (k = zs; k < zs + zm; k++) {
631: for (j=ys; j<ys+ym; j++) {
632: for (i=xs; i < xs+xm; i++) {
633: l[k][j][i][0] = -PETSC_INFINITY;
634: l[k][j][i][1] = 0.0;
635: l[k][j][i][2] = -PETSC_INFINITY;
636: l[k][j][i][3] = 0.0;
637: l[k][j][i][4] = 0.0;
638: u[k][j][i][0] = PETSC_INFINITY;
639: u[k][j][i][1] = 1.0;
640: u[k][j][i][2] = PETSC_INFINITY;
641: u[k][j][i][3] = 1.0;
642: u[k][j][i][4] = 1.0;
643: }
644: }
645: }
647: DMDAVecRestoreArrayDOF(da,xl,&l);
648: DMDAVecRestoreArrayDOF(da,xu,&u);
649: return(0);
650: }
654: PetscErrorCode GetParams(AppCtx *user)
655: {
657: PetscBool flg;
660: /* Set default parameters */
661: user->xmin = 0.0; user->xmax = 64.0;
662: user->ymin = 0.0; user->ymax = 64.0;
663: user->zmin = 0.0; user->zmax = 64.0;
664: user->Dv = 1.0; user->Di=1.0;
665: user->Evf = 0.8; user->Eif = 0.8;
666: user->A = 1.0;
667: user->kBT = 0.11;
668: user->kav = 1.0; user->kai = 1.0; user->kaeta = 1.0;
669: user->Rsurf = 10.0; user->Rbulk = 0.0;
670: user->L = 10.0; user->P_casc = 0.05;
671: user->T = 1.0e-2; user->dt = 1.0e-4;
672: user->VG = 100.0;
673: user->initv = .0001;
674: user->initeta = 0.0;
676: user->degenerate = PETSC_FALSE;
677: /* void growth */
678: user->voidgrowth = PETSC_FALSE;
680: PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);
681: PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);
682: PetscOptionsGetReal(NULL,"-T",&user->T,&flg);
683: PetscOptionsGetReal(NULL,"-dt",&user->dt,&flg);
684: PetscOptionsBool("-degenerate","Run with degenerate mobility\n","None",user->degenerate,&user->degenerate,&flg);
685: PetscOptionsReal("-smallnumber","Small number added to degenerate mobility\n","None",user->smallnumber,&user->smallnumber,&flg);
686: PetscOptionsBool("-voidgrowth","Use initial conditions for void growth\n","None",user->voidgrowth,&user->voidgrowth,&flg);
687: return(0);
688: }
693: PetscErrorCode SetUpMatrices(AppCtx *user)
694: {
696: PetscInt nele,nen,i,j,n;
697: const PetscInt *ele;
698: PetscScalar dt=user->dt,hx,hy,hz;
700: PetscInt idx[8];
701: PetscScalar eM_0[8][8],eM_2[8][8];
702: PetscScalar cv_sum, ci_sum;
703: PetscScalar tp_cv, tp_ci;
705: Mat M =user->M;
706: Mat M_0=user->M_0;
707: PetscInt Xda,Yda,Zda;
708: PetscScalar *cv_p,*ci_p;
709: Vec cvlocal,cilocal;
712: DMGetLocalVector(user->da2,&cvlocal);
713: DMGetLocalVector(user->da2,&cilocal);
714: DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
715: DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
716: DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
717: DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
719: VecGetArray(cvlocal,&cv_p);
720: VecGetArray(cilocal,&ci_p);
722: MatGetLocalSize(M,&n,NULL);
723: DMDAGetInfo(user->da1,NULL,&Xda,&Yda,&Zda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
725: if (Xda!=Yda) printf("Currently different Xda and Yda are not supported");
726: if (user->periodic) {
727: hx = (user->xmax-user->xmin)/Xda;
728: hy = (user->ymax-user->ymin)/Yda;
729: hz = (user->zmax-user->zmin)/Zda;
730: } else {
731: hx = (user->xmax-user->xmin)/(Xda-1.0);
732: hy = (user->ymax-user->ymin)/(Yda-1.0);
733: hz = (user->zmax-user->zmin)/(Zda-1.0);
734: }
736: /** blocks of M_0 **/
737: for (i = 0; i<8; i++) eM_0[i][i] = 1;
739: eM_0[0][1] = 0.5; eM_0[0][2] = 0.25; eM_0[0][3] = 0.5; eM_0[0][4] = 0.5; eM_0[0][5] = 0.25; eM_0[0][6] = 0.125; eM_0[0][7] = 0.25;
740: eM_0[1][2] = 0.5; eM_0[1][3] = 0.25; eM_0[1][4] = 0.25;eM_0[1][5] = 0.5; eM_0[1][6] = 0.25; eM_0[1][7] = 0.125;
741: eM_0[2][3] = 0.5; eM_0[2][4] = 0.125;eM_0[2][5] = 0.25;eM_0[2][6] = 0.5;eM_0[2][7] = 0.25;
742: eM_0[3][4] = 0.25;eM_0[3][5] = 0.125;eM_0[3][6] = 0.25;eM_0[3][7] = 0.5;
743: eM_0[4][5] = 0.5; eM_0[4][6] = 0.25; eM_0[4][7] = 0.5;
744: eM_0[5][6] = 0.5; eM_0[5][7] = 0.25;
745: eM_0[6][7] = 0.5;
747: for (i = 0; i<8; i++) {
748: for (j =0; j<8; j++) {
749: if (i>j) eM_0[i][j] = eM_0[j][i];
750: }
751: }
752:
753: for (i = 0; i<8; i++) {
754: for (j =0; j<8; j++) {
755: eM_0[i][j] = hx*hy*hz/27.0 * eM_0[i][j];
756: }
757: }
759: /** blocks of M_2 **/
761: for (i = 0; i<8; i++) eM_2[i][i] = 12;
763: eM_2[0][1] = 0; eM_2[0][2] = -3; eM_2[0][3] = 0; eM_2[0][4] = 0; eM_2[0][5] = -3; eM_2[0][6] = -3; eM_2[0][7] = -3;
764: eM_2[1][2] = 0; eM_2[1][3] = -3; eM_2[1][4] = -3; eM_2[1][5] = 0; eM_2[1][6] = -3; eM_2[1][7] = -3;
765: eM_2[2][3] = 0; eM_2[2][4] = -3; eM_2[2][5] = -3; eM_2[2][6] = 0; eM_2[2][7] = -3;
766: eM_2[3][4] = -3;eM_2[3][5] = -3; eM_2[3][6] = -3; eM_2[3][7] = 0;
767: eM_2[4][5] = 0; eM_2[4][6] = -3; eM_2[4][7] = 0;
768: eM_2[5][6] = 0; eM_2[5][7] = -3;
769: eM_2[6][7] = 0;
771: for (i = 0; i<8; i++) {
772: for (j =0; j<8; j++) {
773: if (i>j) eM_2[i][j] = eM_2[j][i];
774: }
775: }
777: for (i = 0; i<8; i++) {
778: for (j =0; j<8; j++) {
779: eM_2[i][j] = hx*hy*hz/36.0 * eM_2[i][j];
780: }
781: }
783: /* Get local element info */
784: DMDAGetElements(user->da1,&nele,&nen,&ele);
785: PetscInt row,cols[16],r,row_M_0,cols3[8];
786: PetscScalar vals[16],vals_M_0[8],vals3[8];
788: for (i=0;i < nele;i++) {
789: for (j=0; j<8; j++) idx[j] = ele[8*i + j];
791: for (r=0;r<8;r++) {
792: row_M_0 = idx[r];
793: for (j=0; j<8; j++) vals_M_0[j]=eM_0[r][j];
795: MatSetValuesLocal(M_0,1,&row_M_0,8,idx,vals_M_0,ADD_VALUES);
797: tp_cv = 0.0;
798: tp_ci = 0.0;
799: for (j = 0; j < 8; j++) {
800: tp_cv = tp_cv + cv_p[idx[j]];
801: tp_ci = tp_ci + ci_p[idx[j]];
802: }
804: cv_sum = tp_cv * user->Dv/user->kBT;
805: ci_sum = tp_ci * user->Di/user->kBT;
807: row = 5*idx[r];
808: for (j = 0; j < 8; j++) {
809: cols[j] = 5 * idx[j];
810: vals[j] = dt*eM_2[r][j]*cv_sum;
811: cols[j + 8] = 5 * idx[j] + 1;
812: vals[j + 8] = eM_0[r][j];
813: }
815: /* Insert values in matrix M for 1st dof */
817: MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);
819: row = 5*idx[r]+1;
820: for (j = 0; j < 8; j++) {
821: cols[j] = 5 * idx[j];
822: vals[j] = -eM_0[r][j];
823: cols[j + 8] = 5 * idx[j] + 1;
824: vals[j + 8] = user->kav*eM_2[r][j];
825: }
827: /* Insert values in matrix M for 2nd dof */
828: MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);
830: row = 5*idx[r]+2;
831: for (j = 0; j < 8; j++) {
832: cols[j] = 5 * idx[j] + 2;
833: vals[j] = dt*eM_2[r][j]*ci_sum;
834: cols[j + 8] = 5 * idx[j] + 3;
835: vals[j + 8] = eM_0[r][j];
836: }
838: /* Insert values in matrix M for 3rd dof */
839: MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);
841: row = 5*idx[r]+3;
842: for (j = 0; j < 8; j++) {
843: cols[j] = 5 * idx[j] + 2;
844: vals[j] = -eM_0[r][j];
845: cols[j + 8] = 5 * idx[j] + 3;
846: vals[j + 8] = user->kai*eM_2[r][j];
847: }
849: /* Insert values in matrix M for 4th dof */
850: MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);
852: row = 5*idx[r]+4;
853: for (j = 0; j < 8; j++) {
854: cols3[j] = 5 * idx[j] + 4;
855: vals3[j] = eM_0[r][j]/dt + user->L*user->kaeta*eM_2[r][j];
856: }
858: /* Insert values in matrix M for 5th dof */
859: MatSetValuesLocal(M,1,&row,8,cols3,vals3,ADD_VALUES);
860: }
861: }
864: VecRestoreArray(cvlocal,&cv_p);
865: VecRestoreArray(cilocal,&ci_p);
867: DMRestoreLocalVector(user->da2,&cvlocal);
868: DMRestoreLocalVector(user->da2,&cilocal);
870: MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
871: MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);
873: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
874: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
876: DMDARestoreElements(user->da1,&nele,&nen,&ele);
877: return(0);
878: }
883: PetscErrorCode UpdateMatrices(AppCtx *user)
884: {
886: PetscInt nele,nen,i,j,n,Xda,Yda,Zda;
887: const PetscInt *ele;
889: PetscInt idx[8];
890: PetscScalar eM_2[8][8],h;
891: Mat M=user->M;
892: PetscScalar *cv_p,*ci_p,cv_sum,ci_sum;
893: /* newly added */
894: Vec cvlocal,cilocal;
897: DMGetLocalVector(user->da2,&cvlocal);
898: DMGetLocalVector(user->da2,&cilocal);
899: DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
900: DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
901: DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
902: DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
904: VecGetArray(cvlocal,&cv_p);
905: VecGetArray(cilocal,&ci_p);
907: /* Create the mass matrix M_0 */
908: MatGetLocalSize(M,&n,NULL);
909: DMDAGetInfo(user->da1,NULL,&Xda,&Yda,&Zda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
911: if (Xda!=Yda) printf("Currently different Xda and Yda are not supported");
912: if (user->periodic) h = (user->xmax-user->xmin)/Xda;
913: else h = (user->xmax-user->xmin)/(Xda-1.0);
915: /* Get local element info */
916: DMDAGetElements(user->da1,&nele,&nen,&ele);
918: for (i=0; i<nele; i++) {
919: for (j=0; j<8; j++) idx[j] = ele[8*i + j];
921: PetscInt r,row,cols[8];
922: PetscScalar vals[8];
924: for (r=0; r<8; r++) {
925: row = 5*idx[r];
926: for (j = 0; j < 8; j++) {
927: cols[j] = 5*idx[j];
928: vals[j] = 0.0;
929: }
931: /* Insert values in matrix M for 1st dof */
932: MatSetValuesLocal(M,1,&row,8,cols,vals,INSERT_VALUES);
934: row = 5*idx[r]+2;
935: for (j = 0; j < 8; j++) {
936: cols[j] = 5*idx[j]+2;
937: vals[j] = 0.0;
938: }
940: /* Insert values in matrix M for 3nd dof */
941: MatSetValuesLocal(M,1,&row,8,cols,vals,INSERT_VALUES);
942: }
943: }
945: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
946: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
948: /** blocks of M_2 **/
950: for (i = 0; i<8; i++) eM_2[i][i] = 12;
952: eM_2[0][1] = 0; eM_2[0][2] = -3; eM_2[0][3] = 0; eM_2[0][4] = 0; eM_2[0][5] = -3; eM_2[0][6] = -3; eM_2[0][7] = -3;
953: eM_2[1][2] = 0; eM_2[1][3] = -3; eM_2[1][4] = -3; eM_2[1][5] = 0; eM_2[1][6] = -3; eM_2[1][7] = -3;
954: eM_2[2][3] = 0; eM_2[2][4] = -3; eM_2[2][5] = -3; eM_2[2][6] = 0; eM_2[2][7] = -3;
955: eM_2[3][4] = -3;eM_2[3][5] = -3; eM_2[3][6] = -3; eM_2[3][7] = 0;
956: eM_2[4][5] = 0; eM_2[4][6] = -3; eM_2[4][7] = 0;
957: eM_2[5][6] = 0; eM_2[5][7] = -3;
958: eM_2[6][7] = 0;
960: for (i = 0; i<8; i++) {
961: for (j =0; j<8; j++) {
962: if (i>j) eM_2[i][j] = eM_2[j][i];
963: }
964: }
966: for (i = 0; i<8; i++) {
967: for (j =0; j<8; j++) {
968: eM_2[i][j] = h*h*h/36.0 * eM_2[i][j];
969: }
970: }
972: for (i=0;i<nele;i++) {
973: for (j=0; j<8; j++) idx[j] = ele[8*i + j];
975: PetscInt row,cols[8],r;
976: PetscScalar vals[8];
978: for (r=0;r<8;r++) {
980: if (user->degenerate) {
981: cv_sum = (2.0*user->smallnumber + cv_p[idx[0]] + cv_p[idx[1]])*user->Dv/(2.0*user->kBT);
982: ci_sum = (2.0*user->smallnumber + ci_p[idx[0]] + ci_p[idx[1]])*user->Di/(2.0*user->kBT);
983: } else {
984: cv_sum = user->initv*user->Dv/user->kBT;
985: ci_sum = user->initv*user->Di/user->kBT;
986: }
988: row = 5*idx[r];
989: for (j=0; j<8; j++) {
990: cols[j] = 5*idx[j];
991: vals[j] = user->dt*eM_2[r][j]*cv_sum;
992: }
994: /* Insert values in matrix M for 1st dof */
995: MatSetValuesLocal(M,1,&row,8,cols,vals,ADD_VALUES);
997: row = 5*idx[r]+2;
998: for (j=0; j<8; j++) {
999: cols[j] = 5*idx[j]+2;
1000: vals[j] = user->dt*eM_2[r][j]*ci_sum;
1001: }
1003: /* Insert values in matrix M for 3nd dof */
1004: MatSetValuesLocal(M,1,&row,8,cols,vals,ADD_VALUES);
1005: }
1006: }
1008: VecRestoreArray(cvlocal,&cv_p);
1009: VecRestoreArray(cilocal,&ci_p);
1010: DMRestoreLocalVector(user->da2,&cvlocal);
1011: DMRestoreLocalVector(user->da2,&cilocal);
1013: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1014: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
1015: DMDARestoreElements(user->da1,&nele,&nen,&ele);
1016: return(0);
1017: }