Actual source code: ex633d_db.c
petsc-3.4.5 2014-06-29
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*,MatStructure*,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,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_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,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_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,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_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,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_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,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_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,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_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: DMCreateMatrix(user.da1,MATAIJ,&user.M);
126: /* Get the (usual) mass matrix structure from da2 */
127: DMCreateMatrix(user.da2,MATAIJ,&user.M_0);
129: SetInitialGuess(x,&user);
131: /* Form the jacobian matrix and M_0 */
132: SetUpMatrices(&user);
135: MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
137: SNESCreate(PETSC_COMM_WORLD,&snes);
138: SNESSetDM(snes,user.da1);
140: SNESSetFunction(snes,r,FormFunction,(void*)&user);
141: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
143: SetVariableBounds(user.da1,xl,xu);
144: /* SNESVISetRedundancyCheck(snes,(PetscErrorCode (*)(SNES,IS,IS*,void*))CheckRedundancy,user.da1_clone); */
145: SNESVISetVariableBounds(snes,xl,xu);
146: SNESSetFromOptions(snes);
148: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
149: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_p",FILE_MODE_WRITE,&view_p);
150: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
151: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);
152: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);
155: /*
156: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
157: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);
159: VecView(user.q,view_q);
160: MatView(user.M,view_mat);
162: PetscViewerDestroy(&view_q);
163: PetscViewerDestroy(&view_mat);
164: */
168: while (t<user.T) {
169: char filename[PETSC_MAX_PATH_LEN];
170: PetscScalar a = 1.0;
171: PetscInt i;
172: PetscViewer view;
175: SNESSetFunction(snes,r,FormFunction,(void*)&user);
176: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
178: SetRandomVectors(&user);
179: DPsi(&user);
180: VecView(user.DPsiv,view_psi);
181: VecView(user.DPsii,view_psi);
182: VecView(user.DPsieta,view_psi);
184: VecView(user.Pv,view_p);
185: Update_q(&user);
186: VecView(user.q,view_q);
188: printf("after VecView\n");
189: MatView(user.M,view_mat);
191: printf("after MatView\n");
193: SNESSolve(snes,NULL,x);
195: VecView(x,view_out);
197: PetscInt its;
198: SNESGetIterationNumber(snes,&its);
199: PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4g in %d iterations\n",t,its);
201: Update_u(x,&user);
202: /*
203: for (i=0; i < (int)(user.T/a) ; i++) {
204: if (t/a > i - user.dt/a && t/a < i + user.dt/a) {
205: sprintf(filename,"output_%f",t);
206: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&view);
207: VecView(user.cv,view);
208: VecView(user.ci,view);
209: VecView(user.eta,view);
210: PetscViewerDestroy(&view);
211: }
213: }
214: */
215: UpdateMatrices(&user);
216: t = t + user.dt;
218: }
220: PetscViewerDestroy(&view_out);
221: PetscViewerDestroy(&view_p);
222: PetscViewerDestroy(&view_q);
223: PetscViewerDestroy(&view_psi);
224: PetscViewerDestroy(&view_mat);
225: VecDestroy(&x);
226: VecDestroy(&r);
227: VecDestroy(&xl);
228: VecDestroy(&xu);
229: VecDestroy(&user.q);
230: VecDestroy(&user.wv);
231: VecDestroy(&user.cv);
232: VecDestroy(&user.wi);
233: VecDestroy(&user.ci);
234: VecDestroy(&user.eta);
235: VecDestroy(&user.cvi);
236: VecDestroy(&user.DPsiv);
237: VecDestroy(&user.DPsii);
238: VecDestroy(&user.DPsieta);
239: VecDestroy(&user.Pv);
240: VecDestroy(&user.Pi);
241: VecDestroy(&user.Piv);
242: VecDestroy(&user.Riv);
243: VecDestroy(&user.logcv);
244: VecDestroy(&user.logci);
245: VecDestroy(&user.logcvi);
246: VecDestroy(&user.work1);
247: VecDestroy(&user.work2);
248: VecDestroy(&user.work3);
249: VecDestroy(&user.work4);
250: MatDestroy(&user.M);
251: MatDestroy(&user.M_0);
252: DMDestroy(&user.da1);
253: DMDestroy(&user.da1_clone);
254: DMDestroy(&user.da2);
255: SNESDestroy(&snes);
257: printf("I am finalize \n");
258: PetscFinalize();
259: return 0;
260: }
264: PetscErrorCode Update_u(Vec X,AppCtx *user)
265: {
267: PetscInt i,n;
268: PetscScalar *xx,*wv_p,*cv_p,*wi_p,*ci_p,*eta_p;
271: VecGetLocalSize(user->wv,&n);
273: VecGetArray(X,&xx);
274: VecGetArray(user->wv,&wv_p);
275: VecGetArray(user->cv,&cv_p);
276: VecGetArray(user->wi,&wi_p);
277: VecGetArray(user->ci,&ci_p);
278: VecGetArray(user->eta,&eta_p);
281: for (i=0; i<n; i++) {
282: wv_p[i] = xx[5*i];
283: cv_p[i] = xx[5*i+1];
284: wi_p[i] = xx[5*i+2];
285: ci_p[i] = xx[5*i+3];
286: eta_p[i] = xx[5*i+4];
287: }
288: VecRestoreArray(X,&xx);
289: VecRestoreArray(user->wv,&wv_p);
290: VecRestoreArray(user->cv,&cv_p);
291: VecRestoreArray(user->wi,&wi_p);
292: VecRestoreArray(user->ci,&ci_p);
293: VecRestoreArray(user->eta,&eta_p);
294: return(0);
295: }
299: PetscErrorCode Update_q(AppCtx *user)
300: {
302: PetscScalar *q_p, *w1, *w2;
303: PetscInt i,n;
306: VecPointwiseMult(user->Riv,user->eta,user->eta); /* Riv = eta.^2 */
307: VecScale(user->Riv,user->Rsurf); /* Riv = Rsurf * eta.^2 */
308: VecShift(user->Riv,user->Rbulk); /* Riv = Rbulk + Rsurf * eta.^2 */
309: VecPointwiseMult(user->Riv,user->ci,user->Riv); /* Riv = (Rbulk + Rsurf * eta.^2) * ci */
310: VecPointwiseMult(user->Riv,user->cv,user->Riv); /* Riv = (Rbulk + Rsurf * eta.^2) * ci * cv */
312: VecCopy(user->Riv,user->work1);
313: VecScale(user->work1,user->dt); /* work1 = dt * Riv */
314: VecAXPY(user->work1,-1.0,user->cv); /* work1 = dt * Riv - cv */
315: MatMult(user->M_0,user->work1,user->work2); /* work2 = M_0 * (dt * Riv - cv) */
317: VecGetArray(user->q,&q_p);
318: VecGetArray(user->work1,&w1);
319: VecGetArray(user->work2,&w2);
320: VecGetLocalSize(user->wv,&n);
321: for (i=0; i<n; i++) q_p[5*i]=w2[i];
323: MatMult(user->M_0,user->DPsiv,user->work1);
324: for (i=0; i<n; i++) q_p[5*i+1]=w1[i];
326: VecCopy(user->Riv,user->work1);
327: VecScale(user->work1,user->dt);
328: VecAXPY(user->work1,-1.0,user->ci);
329: MatMult(user->M_0,user->work1,user->work2);
330: for (i=0; i<n; i++) q_p[5*i+2]=w2[i];
332: MatMult(user->M_0,user->DPsii,user->work1);
333: for (i=0; i<n; i++) q_p[5*i+3]=w1[i];
335: VecCopy(user->DPsieta,user->work1);
336: VecScale(user->work1,user->L);
337: VecAXPY(user->work1,-1.0/user->dt,user->eta);
338: MatMult(user->M_0,user->work1,user->work2);
339: for (i=0; i<n; i++) q_p[5*i+4]=w2[i];
341: VecRestoreArray(user->q,&q_p);
342: VecRestoreArray(user->work1,&w1);
343: VecRestoreArray(user->work2,&w2);
344: return(0);
345: }
349: PetscErrorCode DPsi(AppCtx *user)
350: {
352: PetscScalar Evf=user->Evf,Eif=user->Eif,kBT=user->kBT,A=user->A;
353: PetscScalar *cv_p,*ci_p,*eta_p,*logcv_p,*logci_p,*logcvi_p,*DPsiv_p,*DPsii_p,*DPsieta_p;
354: PetscInt n,i;
357: VecGetLocalSize(user->cv,&n);
358: VecGetArray(user->cv,&cv_p);
359: VecGetArray(user->ci,&ci_p);
360: VecGetArray(user->eta,&eta_p);
361: VecGetArray(user->logcv,&logcv_p);
362: VecGetArray(user->logci,&logci_p);
363: VecGetArray(user->logcvi,&logcvi_p);
364: VecGetArray(user->DPsiv,&DPsiv_p);
365: VecGetArray(user->DPsii,&DPsii_p);
366: VecGetArray(user->DPsieta,&DPsieta_p);
368: Llog(user->cv,user->logcv);
369: Llog(user->ci,user->logci);
371: VecCopy(user->cv,user->cvi);
372: VecAXPY(user->cvi,1.0,user->ci);
373: VecScale(user->cvi,-1.0);
374: VecShift(user->cvi,1.0);
375: Llog(user->cvi,user->logcvi);
377: for (i=0; i<n; i++)
378: {
379: 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);
380: 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];
382: 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: }
385: VecRestoreArray(user->cv,&cv_p);
386: VecRestoreArray(user->ci,&ci_p);
387: VecRestoreArray(user->eta,&eta_p);
388: VecGetArray(user->logcv,&logcv_p);
389: VecGetArray(user->logci,&logci_p);
390: VecGetArray(user->logcvi,&logcvi_p);
391: VecRestoreArray(user->DPsiv,&DPsiv_p);
392: VecRestoreArray(user->DPsii,&DPsii_p);
393: VecRestoreArray(user->DPsieta,&DPsieta_p);
394: return(0);
395: }
400: PetscErrorCode Llog(Vec X, Vec Y)
401: {
403: PetscScalar *x,*y;
404: PetscInt n,i;
407: VecGetArray(X,&x);
408: VecGetArray(Y,&y);
409: VecGetLocalSize(X,&n);
410: for (i=0; i<n; i++) {
411: if (x[i] < 1.0e-12) y[i] = log(1.0e-12);
412: else y[i] = log(x[i]);
413: }
414: return(0);
415: }
419: PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
420: {
421: /******************** 3D ********************/
423: PetscInt n,i,j,Xda,Yda,Zda;
424: PetscScalar *xx,*cv_p,*ci_p,*wv_p,*wi_p,*eta_p;
425: PetscViewer view_out;
427: /* needed for the void growth case */
428: 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;
429: PetscScalar xmid,ymid,zmid;
430: PetscInt nele,nen,idx[8];
431: const PetscInt *ele;
432: PetscScalar x[8],y[8],z[8];
433: Vec coords;
434: const PetscScalar *_coords;
435: PetscViewer view;
436: PetscScalar xwidth = user->xmax - user->xmin, ywidth = user->ymax - user->ymin, zwidth = user->zmax - user->zmin;
439: VecGetLocalSize(X,&n);
441: if (user->voidgrowth) {
442: DMDAGetInfo(user->da2,NULL,&Xda,&Yda,&Zda,NULL,NULL,
443: NULL,NULL,NULL,NULL,NULL,NULL,NULL);
444: DMGetCoordinatesLocal(user->da2,&coords);
445: VecGetArrayRead(coords,&_coords);
447: if (user->periodic) h = (user->xmax-user->xmin)/Xda;
448: else h = (user->xmax-user->xmin)/(Xda-1.0);
450: xmid = (user->xmax + user->xmin)/2.0;
451: ymid = (user->ymax + user->ymin)/2.0;
452: zmid = (user->zmax + user->zmin)/2.0;
453: lambda = 4.0*h;
455: DMDAGetElements(user->da2,&nele,&nen,&ele); /* number of local elements, number of element nodes, the local indices of the elements' vertices */
456: for (i=0; i < nele; i++) {
457: for (j = 0; j<8; j++) {
458: idx[j] = ele[8*i + j];
459: x[j] = _coords[3*idx[j]];
460: y[j] = _coords[3*idx[j]+1];
461: z[j] = _coords[3*idx[j]+2];
462: }
464: PetscInt k;
465: PetscScalar vals_cv[8],vals_ci[8],vals_eta[8],s,hhr,r;
466: for (k=0; k < 8; k++) {
467: s = PetscSqrtScalar((x[k] - xmid)*(x[k] - xmid) + (y[k] - ymid)*(y[k] - ymid) + (z[k] - zmid)*(z[k] - zmid));
468: if (s < xwidth*(5.0/64.0)) {
469: vals_cv[k] = cv_v;
470: vals_ci[k] = ci_v;
471: vals_eta[k] = eta_v;
472: } else if (s>= xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0)) {
473: /*r = (s - xwidth*(6.0/64.0))/(0.5*lambda);*/
474: r = (s - xwidth*(6.0/64.0))/(xwidth/64.0);
475: hhr = 0.25*(-r*r*r + 3*r + 2);
476: vals_cv[k] = cv_m + (1.0 - hhr)*(cv_v - cv_m);
477: vals_ci[k] = ci_m + (1.0 - hhr)*(ci_v - ci_m);
478: vals_eta[k] = eta_m + (1.0 - hhr)*(eta_v - eta_m);
479: } else {
480: vals_cv[k] = cv_m;
481: vals_ci[k] = ci_m;
482: vals_eta[k] = eta_m;
483: }
484: }
486: VecSetValuesLocal(user->cv,8,idx,vals_cv,INSERT_VALUES);
487: VecSetValuesLocal(user->ci,8,idx,vals_ci,INSERT_VALUES);
488: VecSetValuesLocal(user->eta,8,idx,vals_eta,INSERT_VALUES);
489: }
490: VecAssemblyBegin(user->cv);
491: VecAssemblyEnd(user->cv);
492: VecAssemblyBegin(user->ci);
493: VecAssemblyEnd(user->ci);
494: VecAssemblyBegin(user->eta);
495: VecAssemblyEnd(user->eta);
497: DMDARestoreElements(user->da2,&nele,&nen,&ele);
498: VecRestoreArrayRead(coords,&_coords);
500: } else {
502: VecSet(user->cv,6.9e-4);
503: VecSet(user->ci,6.9e-4);
504: VecSet(user->eta,0.0);
505: }
507: DPsi(user);
508: VecCopy(user->DPsiv,user->wv);
509: VecCopy(user->DPsii,user->wi);
511: VecGetArray(X,&xx);
512: VecGetArray(user->cv,&cv_p);
513: VecGetArray(user->ci,&ci_p);
514: VecGetArray(user->eta,&eta_p);
515: VecGetArray(user->wv,&wv_p);
516: VecGetArray(user->wi,&wi_p);
517: for (i=0; i<n/5; i++) {
518: xx[5*i] =wv_p[i];
519: xx[5*i+1]=cv_p[i];
520: xx[5*i+2]=wi_p[i];
521: xx[5*i+3]=ci_p[i];
522: xx[5*i+4]=eta_p[i];
523: }
525: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);
526: VecView(user->wv,view_out);
527: VecView(user->cv,view_out);
528: VecView(user->wi,view_out);
529: VecView(user->ci,view_out);
530: VecView(user->eta,view_out);
531: PetscViewerDestroy(&view_out);
533: VecRestoreArray(X,&xx);
534: VecRestoreArray(user->cv,&cv_p);
535: VecRestoreArray(user->ci,&ci_p);
536: VecRestoreArray(user->wv,&wv_p);
537: VecRestoreArray(user->wi,&wi_p);
538: VecRestoreArray(user->eta,&eta_p);
539: return(0);
540: }
544: PetscErrorCode SetRandomVectors(AppCtx *user)
545: {
547: PetscInt i,n,count=0;
548: PetscScalar *w1,*w2,*Pv_p,*eta_p;
550: /* static PetscViewer viewer=0; */
551: static PetscRandom rand = 0;
552: static PetscInt step = 0;
555: if (!rand) {
556: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
557: PetscRandomSetFromOptions(rand);
558: }
560: VecSetRandom(user->work1,rand);
561: VecSetRandom(user->work2,rand);
562: VecGetArray(user->work1,&w1);
563: VecGetArray(user->work2,&w2);
564: VecGetArray(user->Pv,&Pv_p);
565: VecGetArray(user->eta,&eta_p);
566: VecGetLocalSize(user->work1,&n);
567: for (i=0; i<n; i++) {
568: if (eta_p[i]>=0.8 || w1[i]>user->P_casc) Pv_p[i]=0;
569: else {
570: Pv_p[i] = w2[i]*user->VG;
571: count = count + 1;
572: }
573: }
574: step++;
576: VecCopy(user->Pv,user->Pi);
577: VecScale(user->Pi,0.9);
578: VecPointwiseMult(user->Piv,user->Pi,user->Pv);
579: VecRestoreArray(user->work1,&w1);
580: VecRestoreArray(user->work2,&w2);
581: VecRestoreArray(user->Pv,&Pv_p);
582: VecRestoreArray(user->eta,&eta_p);
583: printf("Number of radiations count %d out of total mesh points n %d\n",count,n);
584: return(0);
585: }
589: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
590: {
592: AppCtx *user=(AppCtx*)ctx;
595: MatMultAdd(user->M,X,user->q,F);
596: return(0);
597: }
601: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
602: {
604: AppCtx *user=(AppCtx*)ctx;
607: *flg = SAME_NONZERO_PATTERN;
608: MatCopy(user->M,*J,*flg);
609: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
610: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
611: return(0);
612: }
616: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
617: {
619: PetscScalar ****l,****u;
620: PetscInt xs,xm,ys,ym,zs,zm;
621: PetscInt k,j,i;
624: DMDAVecGetArrayDOF(da,xl,&l);
625: DMDAVecGetArrayDOF(da,xu,&u);
627: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
629: for (k = zs; k < zs + zm; k++) {
630: for (j=ys; j<ys+ym; j++) {
631: for (i=xs; i < xs+xm; i++) {
632: l[k][j][i][0] = -SNES_VI_INF;
633: l[k][j][i][1] = 0.0;
634: l[k][j][i][2] = -SNES_VI_INF;
635: l[k][j][i][3] = 0.0;
636: l[k][j][i][4] = 0.0;
637: u[k][j][i][0] = SNES_VI_INF;
638: u[k][j][i][1] = 1.0;
639: u[k][j][i][2] = SNES_VI_INF;
640: u[k][j][i][3] = 1.0;
641: u[k][j][i][4] = 1.0;
642: }
643: }
644: }
646: DMDAVecRestoreArrayDOF(da,xl,&l);
647: DMDAVecRestoreArrayDOF(da,xu,&u);
648: return(0);
649: }
653: PetscErrorCode GetParams(AppCtx *user)
654: {
656: PetscBool flg;
659: /* Set default parameters */
660: user->xmin = 0.0; user->xmax = 64.0;
661: user->ymin = 0.0; user->ymax = 64.0;
662: user->zmin = 0.0; user->zmax = 64.0;
663: user->Dv = 1.0; user->Di=1.0;
664: user->Evf = 0.8; user->Eif = 0.8;
665: user->A = 1.0;
666: user->kBT = 0.11;
667: user->kav = 1.0; user->kai = 1.0; user->kaeta = 1.0;
668: user->Rsurf = 10.0; user->Rbulk = 0.0;
669: user->L = 10.0; user->P_casc = 0.05;
670: user->T = 1.0e-2; user->dt = 1.0e-4;
671: user->VG = 100.0;
672: user->initv = .0001;
673: user->initeta = 0.0;
675: user->degenerate = PETSC_FALSE;
676: /* void growth */
677: user->voidgrowth = PETSC_FALSE;
679: PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);
680: PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);
681: PetscOptionsGetReal(NULL,"-T",&user->T,&flg);
682: PetscOptionsGetReal(NULL,"-dt",&user->dt,&flg);
683: PetscOptionsBool("-degenerate","Run with degenerate mobility\n","None",user->degenerate,&user->degenerate,&flg);
684: PetscOptionsReal("-smallnumber","Small number added to degenerate mobility\n","None",user->smallnumber,&user->smallnumber,&flg);
685: PetscOptionsBool("-voidgrowth","Use initial conditions for void growth\n","None",user->voidgrowth,&user->voidgrowth,&flg);
686: return(0);
687: }
692: PetscErrorCode SetUpMatrices(AppCtx *user)
693: {
695: PetscInt nele,nen,i,j,n;
696: const PetscInt *ele;
697: PetscScalar dt=user->dt,hx,hy,hz;
699: PetscInt idx[8];
700: PetscScalar eM_0[8][8],eM_2[8][8];
701: PetscScalar cv_sum, ci_sum;
702: PetscScalar tp_cv, tp_ci;
704: Mat M =user->M;
705: Mat M_0=user->M_0;
706: PetscInt Xda,Yda,Zda;
707: PetscScalar *cv_p,*ci_p;
708: Vec cvlocal,cilocal;
711: DMGetLocalVector(user->da2,&cvlocal);
712: DMGetLocalVector(user->da2,&cilocal);
713: DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
714: DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
715: DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
716: DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
718: VecGetArray(cvlocal,&cv_p);
719: VecGetArray(cilocal,&ci_p);
721: MatGetLocalSize(M,&n,NULL);
722: DMDAGetInfo(user->da1,NULL,&Xda,&Yda,&Zda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
724: if (Xda!=Yda) printf("Currently different Xda and Yda are not supported");
725: if (user->periodic) {
726: hx = (user->xmax-user->xmin)/Xda;
727: hy = (user->ymax-user->ymin)/Yda;
728: hz = (user->zmax-user->zmin)/Zda;
729: } else {
730: hx = (user->xmax-user->xmin)/(Xda-1.0);
731: hy = (user->ymax-user->ymin)/(Yda-1.0);
732: hz = (user->zmax-user->zmin)/(Zda-1.0);
733: }
735: /** blocks of M_0 **/
736: for (i = 0; i<8; i++) eM_0[i][i] = 1;
738: 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;
739: 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;
740: 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;
741: eM_0[3][4] = 0.25;eM_0[3][5] = 0.125;eM_0[3][6] = 0.25;eM_0[3][7] = 0.5;
742: eM_0[4][5] = 0.5; eM_0[4][6] = 0.25; eM_0[4][7] = 0.5;
743: eM_0[5][6] = 0.5; eM_0[5][7] = 0.25;
744: eM_0[6][7] = 0.5;
746: for (i = 0; i<8; i++) {
747: for (j =0; j<8; j++) {
748: if (i>j) eM_0[i][j] = eM_0[j][i];
749: }
750: }
751:
752: for (i = 0; i<8; i++) {
753: for (j =0; j<8; j++) {
754: eM_0[i][j] = hx*hy*hz/27.0 * eM_0[i][j];
755: }
756: }
758: /** blocks of M_2 **/
760: for (i = 0; i<8; i++) eM_2[i][i] = 12;
762: 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;
763: 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;
764: 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;
765: eM_2[3][4] = -3;eM_2[3][5] = -3; eM_2[3][6] = -3; eM_2[3][7] = 0;
766: eM_2[4][5] = 0; eM_2[4][6] = -3; eM_2[4][7] = 0;
767: eM_2[5][6] = 0; eM_2[5][7] = -3;
768: eM_2[6][7] = 0;
770: for (i = 0; i<8; i++) {
771: for (j =0; j<8; j++) {
772: if (i>j) eM_2[i][j] = eM_2[j][i];
773: }
774: }
776: for (i = 0; i<8; i++) {
777: for (j =0; j<8; j++) {
778: eM_2[i][j] = hx*hy*hz/36.0 * eM_2[i][j];
779: }
780: }
782: /* Get local element info */
783: DMDAGetElements(user->da1,&nele,&nen,&ele);
784: PetscInt row,cols[16],r,row_M_0,cols3[8];
785: PetscScalar vals[16],vals_M_0[8],vals3[8];
787: for (i=0;i < nele;i++) {
788: for (j=0; j<8; j++) idx[j] = ele[8*i + j];
790: for (r=0;r<8;r++) {
791: row_M_0 = idx[r];
792: for (j=0; j<8; j++) vals_M_0[j]=eM_0[r][j];
794: MatSetValuesLocal(M_0,1,&row_M_0,8,idx,vals_M_0,ADD_VALUES);
796: tp_cv = 0.0;
797: tp_ci = 0.0;
798: for (j = 0; j < 8; j++) {
799: tp_cv = tp_cv + cv_p[idx[j]];
800: tp_ci = tp_ci + ci_p[idx[j]];
801: }
803: cv_sum = tp_cv * user->Dv/user->kBT;
804: ci_sum = tp_ci * user->Di/user->kBT;
806: row = 5*idx[r];
807: for (j = 0; j < 8; j++) {
808: cols[j] = 5 * idx[j];
809: vals[j] = dt*eM_2[r][j]*cv_sum;
810: cols[j + 8] = 5 * idx[j] + 1;
811: vals[j + 8] = eM_0[r][j];
812: }
814: /* Insert values in matrix M for 1st dof */
816: MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);
818: row = 5*idx[r]+1;
819: for (j = 0; j < 8; j++) {
820: cols[j] = 5 * idx[j];
821: vals[j] = -eM_0[r][j];
822: cols[j + 8] = 5 * idx[j] + 1;
823: vals[j + 8] = user->kav*eM_2[r][j];
824: }
826: /* Insert values in matrix M for 2nd dof */
827: MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);
829: row = 5*idx[r]+2;
830: for (j = 0; j < 8; j++) {
831: cols[j] = 5 * idx[j] + 2;
832: vals[j] = dt*eM_2[r][j]*ci_sum;
833: cols[j + 8] = 5 * idx[j] + 3;
834: vals[j + 8] = eM_0[r][j];
835: }
837: /* Insert values in matrix M for 3rd dof */
838: MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);
840: row = 5*idx[r]+3;
841: for (j = 0; j < 8; j++) {
842: cols[j] = 5 * idx[j] + 2;
843: vals[j] = -eM_0[r][j];
844: cols[j + 8] = 5 * idx[j] + 3;
845: vals[j + 8] = user->kai*eM_2[r][j];
846: }
848: /* Insert values in matrix M for 4th dof */
849: MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);
851: row = 5*idx[r]+4;
852: for (j = 0; j < 8; j++) {
853: cols3[j] = 5 * idx[j] + 4;
854: vals3[j] = eM_0[r][j]/dt + user->L*user->kaeta*eM_2[r][j];
855: }
857: /* Insert values in matrix M for 5th dof */
858: MatSetValuesLocal(M,1,&row,8,cols3,vals3,ADD_VALUES);
859: }
860: }
863: VecRestoreArray(cvlocal,&cv_p);
864: VecRestoreArray(cilocal,&ci_p);
866: DMRestoreLocalVector(user->da2,&cvlocal);
867: DMRestoreLocalVector(user->da2,&cilocal);
869: MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
870: MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);
872: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
873: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
875: DMDARestoreElements(user->da1,&nele,&nen,&ele);
876: return(0);
877: }
882: PetscErrorCode UpdateMatrices(AppCtx *user)
883: {
885: PetscInt nele,nen,i,j,n,Xda,Yda,Zda;
886: const PetscInt *ele;
888: PetscInt idx[8];
889: PetscScalar eM_2[8][8],h;
890: Mat M=user->M;
891: PetscScalar *cv_p,*ci_p,cv_sum,ci_sum;
892: /* newly added */
893: Vec cvlocal,cilocal;
896: DMGetLocalVector(user->da2,&cvlocal);
897: DMGetLocalVector(user->da2,&cilocal);
898: DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
899: DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
900: DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
901: DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
903: VecGetArray(cvlocal,&cv_p);
904: VecGetArray(cilocal,&ci_p);
906: /* Create the mass matrix M_0 */
907: MatGetLocalSize(M,&n,NULL);
908: DMDAGetInfo(user->da1,NULL,&Xda,&Yda,&Zda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
910: if (Xda!=Yda) printf("Currently different Xda and Yda are not supported");
911: if (user->periodic) h = (user->xmax-user->xmin)/Xda;
912: else h = (user->xmax-user->xmin)/(Xda-1.0);
914: /* Get local element info */
915: DMDAGetElements(user->da1,&nele,&nen,&ele);
917: for (i=0; i<nele; i++) {
918: for (j=0; j<8; j++) idx[j] = ele[8*i + j];
920: PetscInt r,row,cols[8];
921: PetscScalar vals[8];
923: for (r=0; r<8; r++) {
924: row = 5*idx[r];
925: for (j = 0; j < 8; j++) {
926: cols[j] = 5*idx[j];
927: vals[j] = 0.0;
928: }
930: /* Insert values in matrix M for 1st dof */
931: MatSetValuesLocal(M,1,&row,8,cols,vals,INSERT_VALUES);
933: row = 5*idx[r]+2;
934: for (j = 0; j < 8; j++) {
935: cols[j] = 5*idx[j]+2;
936: vals[j] = 0.0;
937: }
939: /* Insert values in matrix M for 3nd dof */
940: MatSetValuesLocal(M,1,&row,8,cols,vals,INSERT_VALUES);
941: }
942: }
944: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
945: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
947: /** blocks of M_2 **/
949: for (i = 0; i<8; i++) eM_2[i][i] = 12;
951: 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;
952: 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;
953: 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;
954: eM_2[3][4] = -3;eM_2[3][5] = -3; eM_2[3][6] = -3; eM_2[3][7] = 0;
955: eM_2[4][5] = 0; eM_2[4][6] = -3; eM_2[4][7] = 0;
956: eM_2[5][6] = 0; eM_2[5][7] = -3;
957: eM_2[6][7] = 0;
959: for (i = 0; i<8; i++) {
960: for (j =0; j<8; j++) {
961: if (i>j) eM_2[i][j] = eM_2[j][i];
962: }
963: }
965: for (i = 0; i<8; i++) {
966: for (j =0; j<8; j++) {
967: eM_2[i][j] = h*h*h/36.0 * eM_2[i][j];
968: }
969: }
971: for (i=0;i<nele;i++) {
972: for (j=0; j<8; j++) idx[j] = ele[8*i + j];
974: PetscInt row,cols[8],r;
975: PetscScalar vals[8];
977: for (r=0;r<8;r++) {
979: if (user->degenerate) {
980: cv_sum = (2.0*user->smallnumber + cv_p[idx[0]] + cv_p[idx[1]])*user->Dv/(2.0*user->kBT);
981: ci_sum = (2.0*user->smallnumber + ci_p[idx[0]] + ci_p[idx[1]])*user->Di/(2.0*user->kBT);
982: } else {
983: cv_sum = user->initv*user->Dv/user->kBT;
984: ci_sum = user->initv*user->Di/user->kBT;
985: }
987: row = 5*idx[r];
988: for (j=0; j<8; j++) {
989: cols[j] = 5*idx[j];
990: vals[j] = user->dt*eM_2[r][j]*cv_sum;
991: }
993: /* Insert values in matrix M for 1st dof */
994: MatSetValuesLocal(M,1,&row,8,cols,vals,ADD_VALUES);
996: row = 5*idx[r]+2;
997: for (j=0; j<8; j++) {
998: cols[j] = 5*idx[j]+2;
999: vals[j] = user->dt*eM_2[r][j]*ci_sum;
1000: }
1002: /* Insert values in matrix M for 3nd dof */
1003: MatSetValuesLocal(M,1,&row,8,cols,vals,ADD_VALUES);
1004: }
1005: }
1007: VecRestoreArray(cvlocal,&cv_p);
1008: VecRestoreArray(cilocal,&ci_p);
1009: DMRestoreLocalVector(user->da2,&cvlocal);
1010: DMRestoreLocalVector(user->da2,&cilocal);
1012: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1013: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
1014: DMDARestoreElements(user->da1,&nele,&nen,&ele);
1015: return(0);
1016: }