Actual source code: ex633D_DB.c
petsc-3.3-p7 2013-05-11
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);
68:
70: // Get physics and time parameters
71: GetParams(&user);
72:
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,PETSC_NULL,PETSC_NULL,PETSC_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,PETSC_NULL,PETSC_NULL,PETSC_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,PETSC_NULL,PETSC_NULL,PETSC_NULL,&user.da2);
77:
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,PETSC_NULL,PETSC_NULL,PETSC_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,PETSC_NULL,PETSC_NULL,PETSC_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,PETSC_NULL,PETSC_NULL,PETSC_NULL,&user.da2);
84: }
85:
86: // Set Element type (rectangular)
87: DMDASetElementType(user.da1,DMDA_ELEMENT_Q1);
88: DMDASetElementType(user.da2,DMDA_ELEMENT_Q1);
89:
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);
93:
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);
100:
101:
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);
123:
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);
128:
129: SetInitialGuess(x,&user);
130:
131: // Form the jacobian matrix and M_0
132: SetUpMatrices(&user);
133:
134:
135: MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
136:
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);
147:
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);
153:
154:
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);
158:
159: VecView(user.q,view_q);
160: MatView(user.M,view_mat);
161:
162: PetscViewerDestroy(&view_q);
163: PetscViewerDestroy(&view_mat);
164: */
165:
166:
167:
168: while(t<user.T) {
169: char filename[PETSC_MAX_PATH_LEN];
170: PetscScalar a = 1.0;
171: PetscInt i;
172: PetscViewer view;
173:
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);
187:
188: printf("after VecView\n");
189: MatView(user.M,view_mat);
190:
191: printf("after MatView\n");
193: SNESSolve(snes,PETSC_NULL,x);
194:
195: VecView(x,view_out);
196:
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: }
212:
213: }
214: */
215: UpdateMatrices(&user);
216: t = t + user.dt;
217:
218: }
219:
221:
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);
258:
259:
260:
261: printf("I am finalize \n");
262: PetscFinalize();
263: return 0;
264: }
268: PetscErrorCode Update_u(Vec X,AppCtx *user)
269: {
271: PetscInt i,n;
272: PetscScalar *xx,*wv_p,*cv_p,*wi_p,*ci_p,*eta_p;
273:
274:
276: VecGetLocalSize(user->wv,&n);
278: VecGetArray(X,&xx);
279: VecGetArray(user->wv,&wv_p);
280: VecGetArray(user->cv,&cv_p);
281: VecGetArray(user->wi,&wi_p);
282: VecGetArray(user->ci,&ci_p);
283: VecGetArray(user->eta,&eta_p);
284:
285:
286: for(i=0;i<n;i++) {
287: wv_p[i] = xx[5*i];
288: cv_p[i] = xx[5*i+1];
289: wi_p[i] = xx[5*i+2];
290: ci_p[i] = xx[5*i+3];
291: eta_p[i] = xx[5*i+4];
292: }
293: VecRestoreArray(X,&xx);
294: VecRestoreArray(user->wv,&wv_p);
295: VecRestoreArray(user->cv,&cv_p);
296: VecRestoreArray(user->wi,&wi_p);
297: VecRestoreArray(user->ci,&ci_p);
298: VecRestoreArray(user->eta,&eta_p);
299:
300: return(0);
301: }
305: PetscErrorCode Update_q(AppCtx *user)
306: {
308: PetscScalar *q_p, *w1, *w2;
309: PetscInt i,n;
312:
313: VecPointwiseMult(user->Riv,user->eta,user->eta); //Riv = eta.^2
314: VecScale(user->Riv,user->Rsurf); // Riv = Rsurf * eta.^2
315: VecShift(user->Riv,user->Rbulk); // Riv = Rbulk + Rsurf * eta.^2
316: VecPointwiseMult(user->Riv,user->ci,user->Riv); // Riv = (Rbulk + Rsurf * eta.^2) * ci
317: VecPointwiseMult(user->Riv,user->cv,user->Riv);// Riv = (Rbulk + Rsurf * eta.^2) * ci * cv
318:
319: VecCopy(user->Riv,user->work1);
320: VecScale(user->work1,user->dt); // work1 = dt * Riv
321: VecAXPY(user->work1,-1.0,user->cv);// work1 = dt * Riv - cv
322: MatMult(user->M_0,user->work1,user->work2); // work2 = M_0 * (dt * Riv - cv)
323:
324: VecGetArray(user->q,&q_p);
325: VecGetArray(user->work1,&w1);
326: VecGetArray(user->work2,&w2);
327: VecGetLocalSize(user->wv,&n);
328: for (i=0;i<n;i++) {
329: q_p[5*i]=w2[i];
330: }
332: MatMult(user->M_0,user->DPsiv,user->work1);
333: for (i=0;i<n;i++) {
334: q_p[5*i+1]=w1[i];
335: }
337: VecCopy(user->Riv,user->work1);
338: VecScale(user->work1,user->dt);
339: VecAXPY(user->work1,-1.0,user->ci);
340: MatMult(user->M_0,user->work1,user->work2);
341: for (i=0;i<n;i++) {
342: q_p[5*i+2]=w2[i];
343: }
345: MatMult(user->M_0,user->DPsii,user->work1);
346: for (i=0;i<n;i++) {
347: q_p[5*i+3]=w1[i];
348: }
350: VecCopy(user->DPsieta,user->work1);
351: VecScale(user->work1,user->L);
352: VecAXPY(user->work1,-1.0/user->dt,user->eta);
353: MatMult(user->M_0,user->work1,user->work2);
354: for (i=0;i<n;i++) {
355: q_p[5*i+4]=w2[i];
356: }
357:
358: VecRestoreArray(user->q,&q_p);
359: VecRestoreArray(user->work1,&w1);
360: VecRestoreArray(user->work2,&w2);
363:
364: return(0);
365: }
369: PetscErrorCode DPsi(AppCtx* user)
370: {
371: PetscErrorCode ierr;
372: PetscScalar Evf=user->Evf,Eif=user->Eif,kBT=user->kBT,A=user->A;
373: PetscScalar *cv_p,*ci_p,*eta_p,*logcv_p,*logci_p,*logcvi_p,*DPsiv_p,*DPsii_p,*DPsieta_p;
374: PetscInt n,i;
378: VecGetLocalSize(user->cv,&n);
379: VecGetArray(user->cv,&cv_p);
380: VecGetArray(user->ci,&ci_p);
381: VecGetArray(user->eta,&eta_p);
382: VecGetArray(user->logcv,&logcv_p);
383: VecGetArray(user->logci,&logci_p);
384: VecGetArray(user->logcvi,&logcvi_p);
385: VecGetArray(user->DPsiv,&DPsiv_p);
386: VecGetArray(user->DPsii,&DPsii_p);
387: VecGetArray(user->DPsieta,&DPsieta_p);
389: Llog(user->cv,user->logcv);
390: Llog(user->ci,user->logci);
392: VecCopy(user->cv,user->cvi);
393: VecAXPY(user->cvi,1.0,user->ci);
394: VecScale(user->cvi,-1.0);
395: VecShift(user->cvi,1.0);
396: Llog(user->cvi,user->logcvi);
398: for (i=0;i<n;i++)
399: {
400: 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);
402: 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] ;
403:
404: 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]);
406:
407: }
409: VecRestoreArray(user->cv,&cv_p);
410: VecRestoreArray(user->ci,&ci_p);
411: VecRestoreArray(user->eta,&eta_p);
412: VecGetArray(user->logcv,&logcv_p);
413: VecGetArray(user->logci,&logci_p);
414: VecGetArray(user->logcvi,&logcvi_p);
415: VecRestoreArray(user->DPsiv,&DPsiv_p);
416: VecRestoreArray(user->DPsii,&DPsii_p);
417: VecRestoreArray(user->DPsieta,&DPsieta_p);
420: return(0);
422: }
427: PetscErrorCode Llog(Vec X, Vec Y)
428: {
429: PetscErrorCode ierr;
430: PetscScalar *x,*y;
431: PetscInt n,i;
434:
435: VecGetArray(X,&x);
436: VecGetArray(Y,&y);
437: VecGetLocalSize(X,&n);
438: for (i=0;i<n;i++) {
439: if (x[i] < 1.0e-12) {
440: y[i] = log(1.0e-12);
441: }
442: else {
443: y[i] = log(x[i]);
444: }
445: }
446:
447: return(0);
448: }
452: PetscErrorCode SetInitialGuess(Vec X,AppCtx* user)
453: {
454: /////////////////////////// 3D //////////////////////////////
455: PetscErrorCode ierr;
456: PetscInt n,i,j,Xda,Yda,Zda;
457: PetscScalar *xx,*cv_p,*ci_p,*wv_p,*wi_p,*eta_p;
458: PetscViewer view_out;
459: // needed for the void growth case
460: 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;
461: PetscScalar xmid,ymid,zmid;
462: PetscInt nele,nen,idx[8];
463: const PetscInt *ele;
464: PetscScalar x[8],y[8],z[8];
465: Vec coords;
466: const PetscScalar *_coords;
467: PetscViewer view;
468: PetscScalar xwidth = user->xmax - user->xmin, ywidth = user->ymax - user->ymin, zwidth = user->zmax - user->zmin;
469:
471:
472: VecGetLocalSize(X,&n);
473:
474: if (user->voidgrowth) {
475: DMDAGetInfo(user->da2,PETSC_NULL,&Xda,&Yda,&Zda,PETSC_NULL,PETSC_NULL,
476: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
477: DMDAGetGhostedCoordinates(user->da2,&coords);
478: VecGetArrayRead(coords,&_coords);
479:
480: if (user->periodic) {
481: h = (user->xmax-user->xmin)/Xda;
482: } else {
483: h = (user->xmax-user->xmin)/(Xda-1.0);
484: }
485:
486: xmid = (user->xmax + user->xmin)/2.0;
487: ymid = (user->ymax + user->ymin)/2.0;
488: zmid = (user->zmax + user->zmin)/2.0;
489: lambda = 4.0*h;
490:
491: DMDAGetElements(user->da2,&nele,&nen,&ele); // number of local elements, number of element nodes, the local indices of the elements' vertices
492: for (i=0;i < nele; i++)
493: {
494: for (j = 0; j<8; j++)
495: {
496: idx[j] = ele[8*i + j];
497: x[j] = _coords[3*idx[j]];
498: y[j] = _coords[3*idx[j]+1];
499: z[j] = _coords[3*idx[j]+2];
500: }
501:
502: PetscInt k;
503: PetscScalar vals_cv[8],vals_ci[8],vals_eta[8],s,hhr,r;
504: for (k=0; k < 8 ; k++)
505: {
506: s = sqrt((x[k] - xmid)*(x[k] - xmid) + (y[k] - ymid)*(y[k] - ymid) + (z[k] - zmid)*(z[k] - zmid));
507: if (s < xwidth*(5.0/64.0))
508: {
509: vals_cv[k] = cv_v;
510: vals_ci[k] = ci_v;
511: vals_eta[k] = eta_v;
512: }
513: else if (s>= xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0) )
514: {
515: //r = (s - xwidth*(6.0/64.0) )/(0.5*lambda);
516: r = (s - xwidth*(6.0/64.0) )/(xwidth/64.0);
517: hhr = 0.25*(-r*r*r + 3*r + 2);
518: vals_cv[k] = cv_m + (1.0 - hhr)*(cv_v - cv_m);
519: vals_ci[k] = ci_m + (1.0 - hhr)*(ci_v - ci_m);
520: vals_eta[k] = eta_m + (1.0 - hhr)*(eta_v - eta_m);
521: } else
522: {
523: vals_cv[k] = cv_m;
524: vals_ci[k] = ci_m;
525: vals_eta[k] = eta_m;
526: }
527: }
528:
529: VecSetValuesLocal(user->cv,8,idx,vals_cv,INSERT_VALUES);
530: VecSetValuesLocal(user->ci,8,idx,vals_ci,INSERT_VALUES);
531: VecSetValuesLocal(user->eta,8,idx,vals_eta,INSERT_VALUES);
532: }
533: VecAssemblyBegin(user->cv);
534: VecAssemblyEnd(user->cv);
535: VecAssemblyBegin(user->ci);
536: VecAssemblyEnd(user->ci);
537: VecAssemblyBegin(user->eta);
538: VecAssemblyEnd(user->eta);
539:
540: DMDARestoreElements(user->da2,&nele,&nen,&ele);
541: VecRestoreArrayRead(coords,&_coords);
542:
543: }else {
544:
545: VecSet(user->cv,6.9e-4);
546: VecSet(user->ci,6.9e-4);
547: VecSet(user->eta,0.0);
548: }
550: DPsi(user);
551: VecCopy(user->DPsiv,user->wv);
552: VecCopy(user->DPsii,user->wi);
554: VecGetArray(X,&xx);
555: VecGetArray(user->cv,&cv_p);
556: VecGetArray(user->ci,&ci_p);
557: VecGetArray(user->eta,&eta_p);
558: VecGetArray(user->wv,&wv_p);
559: VecGetArray(user->wi,&wi_p);
560: for (i=0;i<n/5;i++)
561: {
562: xx[5*i]=wv_p[i];
563: xx[5*i+1]=cv_p[i];
564: xx[5*i+2]=wi_p[i];
565: xx[5*i+3]=ci_p[i];
566: xx[5*i+4]=eta_p[i];
567: }
569: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view_out);
570: VecView(user->wv,view_out);
571: VecView(user->cv,view_out);
572: VecView(user->wi,view_out);
573: VecView(user->ci,view_out);
574: VecView(user->eta,view_out);
575: PetscViewerDestroy(&view_out);
577: VecRestoreArray(X,&xx);
578: VecRestoreArray(user->cv,&cv_p);
579: VecRestoreArray(user->ci,&ci_p);
580: VecRestoreArray(user->wv,&wv_p);
581: VecRestoreArray(user->wi,&wi_p);
582: VecRestoreArray(user->eta,&eta_p);
583: return(0);
584: }
588: PetscErrorCode SetRandomVectors(AppCtx* user)
589: {
591: PetscInt i,n,count=0;
592: PetscScalar *w1,*w2,*Pv_p,*eta_p;
594: /* static PetscViewer viewer=0; */
595: static PetscRandom rand = 0;
596: static PetscInt step = 0;
599: if (!rand) {
600: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
601: PetscRandomSetFromOptions(rand);
602: }
603:
604: VecSetRandom(user->work1,rand);
605: VecSetRandom(user->work2,rand);
606: VecGetArray(user->work1,&w1);
607: VecGetArray(user->work2,&w2);
608: VecGetArray(user->Pv,&Pv_p);
609: VecGetArray(user->eta,&eta_p);
610: VecGetLocalSize(user->work1,&n);
611: for (i=0;i<n;i++) {
612: if (eta_p[i]>=0.8 || w1[i]>user->P_casc){
613: Pv_p[i]=0;
614: }
615: else
616: {
617: Pv_p[i]=w2[i]*user->VG;
618: count = count + 1;
619: }
620: }
621: step++;
623: VecCopy(user->Pv,user->Pi);
624: VecScale(user->Pi,0.9);
625: VecPointwiseMult(user->Piv,user->Pi,user->Pv);
626: VecRestoreArray(user->work1,&w1);
627: VecRestoreArray(user->work2,&w2);
628: VecRestoreArray(user->Pv,&Pv_p);
629: VecRestoreArray(user->eta,&eta_p);
630: printf("Number of radiations count %d out of total mesh points n %d\n",count,n);
632: return(0);
633:
634: }
637: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ctx)
638: {
640: AppCtx *user=(AppCtx*)ctx;
641:
643: MatMultAdd(user->M,X,user->q,F);
644: return(0);
645: }
649: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
650: {
652: AppCtx *user=(AppCtx*)ctx;
653:
655: *flg = SAME_NONZERO_PATTERN;
656: MatCopy(user->M,*J,*flg);
657: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
658: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
659: return(0);
660: }
663: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
664: {
666: PetscScalar ****l,****u;
667: PetscInt xs,xm,ys,ym,zs,zm;
668: PetscInt k,j,i;
669:
671: DMDAVecGetArrayDOF(da,xl,&l);
672: DMDAVecGetArrayDOF(da,xu,&u);
673:
674: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
675:
676: for (k = zs; k < zs + zm; k++)
677: {
678: for (j=ys; j<ys+ym; j++)
679: {
680: for(i=xs; i < xs+xm;i++)
681: {
682: l[k][j][i][0] = -SNES_VI_INF;
683: l[k][j][i][1] = 0.0;
684: l[k][j][i][2] = -SNES_VI_INF;
685: l[k][j][i][3] = 0.0;
686: l[k][j][i][4] = 0.0;
687: u[k][j][i][0] = SNES_VI_INF;
688: u[k][j][i][1] = 1.0;
689: u[k][j][i][2] = SNES_VI_INF;
690: u[k][j][i][3] = 1.0;
691: u[k][j][i][4] = 1.0;
692: }
693: }
694: }
695:
696:
697: DMDAVecRestoreArrayDOF(da,xl,&l);
698: DMDAVecRestoreArrayDOF(da,xu,&u);
699: return(0);
700: }
704: PetscErrorCode GetParams(AppCtx* user)
705: {
707: PetscBool flg;
708:
710:
711: /* Set default parameters */
712: user->xmin = 0.0; user->xmax = 64.0;
713: user->ymin = 0.0; user->ymax = 64.0;
714: user->zmin = 0.0; user->zmax = 64.0;
715: user->Dv = 1.0; user->Di=1.0;
716: user->Evf = 0.8; user->Eif = 0.8;
717: user->A = 1.0;
718: user->kBT = 0.11;
719: user->kav = 1.0; user->kai = 1.0; user->kaeta = 1.0;
720: user->Rsurf = 10.0; user->Rbulk = 0.0;
721: user->L = 10.0; user->P_casc = 0.05;
722: user->T = 1.0e-2; user->dt = 1.0e-4;
723: user->VG = 100.0;
724: user->initv = .0001;
725: user->initeta = 0.0;
726:
727: user->degenerate = PETSC_FALSE;
728: /* void growth */
729: user->voidgrowth = PETSC_FALSE;
731: PetscOptionsGetReal(PETSC_NULL,"-xmin",&user->xmin,&flg);
732: PetscOptionsGetReal(PETSC_NULL,"-xmax",&user->xmax,&flg);
733: PetscOptionsGetReal(PETSC_NULL,"-T",&user->T,&flg);
734: PetscOptionsGetReal(PETSC_NULL,"-dt",&user->dt,&flg);
735: PetscOptionsBool("-degenerate","Run with degenerate mobility\n","None",user->degenerate,&user->degenerate,&flg);
736: PetscOptionsReal("-smallnumber","Small number added to degenerate mobility\n","None",user->smallnumber,&user->smallnumber,&flg);
737: PetscOptionsBool("-voidgrowth","Use initial conditions for void growth\n","None",user->voidgrowth,&user->voidgrowth,&flg);
738:
739: return(0);
740: }
745: PetscErrorCode SetUpMatrices(AppCtx* user)
746: {
747: PetscErrorCode ierr;
748: PetscInt nele,nen,i,j,n;
749: const PetscInt *ele;
750: PetscScalar dt=user->dt,hx,hy,hz;
751:
752: PetscInt idx[8];
753: PetscScalar eM_0[8][8],eM_2[8][8];
754: PetscScalar cv_sum, ci_sum;
755: PetscScalar tp_cv, tp_ci;
756:
757: Mat M=user->M;
758: Mat M_0=user->M_0;
759: PetscInt Xda,Yda,Zda;
760: PetscScalar *cv_p,*ci_p;
761: Vec cvlocal,cilocal;
762:
765: DMGetLocalVector(user->da2,&cvlocal);
766: DMGetLocalVector(user->da2,&cilocal);
767: DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
768: DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
769: DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
770: DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
772: VecGetArray(cvlocal,&cv_p);
773: VecGetArray(cilocal,&ci_p);
775: MatGetLocalSize(M,&n,PETSC_NULL);
776: DMDAGetInfo(user->da1,PETSC_NULL,&Xda,&Yda,&Zda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
778: if (Xda!=Yda) {
779: printf("Currently different Xda and Yda are not supported");
780: }
781: if (user->periodic) {
782: hx = (user->xmax-user->xmin)/Xda;
783: hy = (user->ymax-user->ymin)/Yda;
784: hz = (user->zmax-user->zmin)/Zda;
785: } else {
786: hx = (user->xmax-user->xmin)/(Xda-1.0);
787: hy = (user->ymax-user->ymin)/(Yda-1.0);
788: hz = (user->zmax-user->zmin)/(Zda-1.0);
789: }
790:
791:
792:
793: /** blocks of M_0 **/
794:
795:
796: for (i = 0; i<8; i++)
797: {
798: eM_0[i][i] = 1;
799: }
800: 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;
801: 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;
802: 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;
803: eM_0[3][4] = 0.25;eM_0[3][5] = 0.125;eM_0[3][6] = 0.25;eM_0[3][7] = 0.5;
804: eM_0[4][5] = 0.5; eM_0[4][6] = 0.25; eM_0[4][7] = 0.5;
805: eM_0[5][6] = 0.5; eM_0[5][7] = 0.25;
806: eM_0[6][7] = 0.5;
807:
808: for (i = 0; i<8; i++)
809: {
810: for (j =0; j<8; j++)
811: {
812: if (i>j)
813: {
814: eM_0[i][j] = eM_0[j][i];
815: }
816: }
817: }
818:
819: for (i = 0; i<8; i++)
820: {
821: for (j =0; j<8; j++)
822: {
823: eM_0[i][j] = hx*hy*hz/27.0 * eM_0[i][j];
824: }
825: }
826:
827: /** blocks of M_2 **/
828:
829: for (i = 0; i<8; i++)
830: {
831: eM_2[i][i] = 12;
832: }
833:
834: 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;
835: 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;
836: 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;
837: eM_2[3][4] = -3;eM_2[3][5] = -3; eM_2[3][6] = -3; eM_2[3][7] = 0;
838: eM_2[4][5] = 0; eM_2[4][6] = -3; eM_2[4][7] = 0;
839: eM_2[5][6] = 0; eM_2[5][7] = -3;
840: eM_2[6][7] = 0;
841:
842: for (i = 0; i<8; i++)
843: {
844: for (j =0; j<8; j++)
845: {
846: if (i>j)
847: {
848: eM_2[i][j] = eM_2[j][i];
849: }
850: }
851: }
852:
853:
854: for (i = 0; i<8; i++)
855: {
856: for (j =0; j<8; j++)
857: {
858: eM_2[i][j] = hx*hy*hz/36.0 * eM_2[i][j];
859: }
860: }
861:
862:
863: /* Get local element info */
864: DMDAGetElements(user->da1,&nele,&nen,&ele);
865: PetscInt row,cols[16],r,row_M_0,cols3[8];
866: PetscScalar vals[16],vals_M_0[8],vals3[8];
867:
868:
869: for(i=0;i < nele;i++)
870: {
871: for (j=0; j<8; j++)
872: {
873: idx[j] = ele[8*i + j];
874: }
875:
876: for(r=0;r<8;r++)
877: {
878: row_M_0 = idx[r];
879: for (j=0; j<8; j++)
880: {
881: vals_M_0[j]=eM_0[r][j];
882: }
883:
884:
885: MatSetValuesLocal(M_0,1,&row_M_0,8,idx,vals_M_0,ADD_VALUES);
886:
887:
888: tp_cv = 0.0;
889: tp_ci = 0.0;
890: for (j = 0; j < 8; j++)
891: {
892: tp_cv = tp_cv + cv_p[idx[j]];
893: tp_ci = tp_ci + ci_p[idx[j]];
894: }
895:
896: cv_sum = tp_cv * user->Dv/user->kBT;
897: ci_sum = tp_ci * user->Di/user->kBT;
898:
899: row = 5*idx[r];
900: for(j = 0; j < 8; j++)
901: {
902: cols[j] = 5 * idx[j];
903: vals[j] = dt*eM_2[r][j]*cv_sum;
904: cols[j + 8] = 5 * idx[j] + 1;
905: vals[j + 8] = eM_0[r][j];
906: }
907:
908: // Insert values in matrix M for 1st dof
909:
910: MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);
911:
912:
913: row = 5*idx[r]+1;
914: for(j = 0; j < 8; j++)
915: {
916: cols[j] = 5 * idx[j];
917: vals[j] = -eM_0[r][j];
918: cols[j + 8] = 5 * idx[j] + 1;
919: vals[j + 8] = user->kav*eM_2[r][j];
920: }
921:
922: // Insert values in matrix M for 2nd dof
923: MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);
924:
925: row = 5*idx[r]+2;
926: for(j = 0; j < 8; j++)
927: {
928: cols[j] = 5 * idx[j] + 2;
929: vals[j] = dt*eM_2[r][j]*ci_sum;
930: cols[j + 8] = 5 * idx[j] + 3;
931: vals[j + 8] = eM_0[r][j];
932: }
933:
934: // Insert values in matrix M for 3rd dof
935: MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);
936:
937:
938: row = 5*idx[r]+3;
939: for(j = 0; j < 8; j++)
940: {
941: cols[j] = 5 * idx[j] + 2;
942: vals[j] = -eM_0[r][j];
943: cols[j + 8] = 5 * idx[j] + 3;
944: vals[j + 8] = user->kai*eM_2[r][j];
945: }
946:
947: // Insert values in matrix M for 4th dof
948: MatSetValuesLocal(M,1,&row,16,cols,vals,ADD_VALUES);
949:
950: row = 5*idx[r]+4;
951: for(j = 0; j < 8; j++)
952: {
953: cols3[j] = 5 * idx[j] + 4;
954: vals3[j] = eM_0[r][j]/dt + user->L*user->kaeta*eM_2[r][j];
955: }
956:
957: // Insert values in matrix M for 5th dof
958: MatSetValuesLocal(M,1,&row,8,cols3,vals3,ADD_VALUES);
959: }
960: } //
962:
963: VecRestoreArray(cvlocal,&cv_p);
964: VecRestoreArray(cilocal,&ci_p);
965:
966: DMRestoreLocalVector(user->da2,&cvlocal);
967: DMRestoreLocalVector(user->da2,&cilocal);
968:
969: MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
970: MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);
971:
972: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
973: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
974:
975: DMDARestoreElements(user->da1,&nele,&nen,&ele);
976:
977:
978: return(0);
979: }
984: PetscErrorCode UpdateMatrices(AppCtx* user)
985: {
986: PetscErrorCode ierr;
987: PetscInt nele,nen,i,j,n,Xda,Yda,Zda;
988: const PetscInt *ele;
989:
990: PetscInt idx[8];
991: PetscScalar eM_2[8][8],h;
992: Mat M=user->M;
993: PetscScalar *cv_p,*ci_p,cv_sum,ci_sum;
994: /* newly added */
995: Vec cvlocal,cilocal;
998:
999: DMGetLocalVector(user->da2,&cvlocal);
1000: DMGetLocalVector(user->da2,&cilocal);
1001: DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
1002: DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
1003: DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
1004: DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
1005:
1006: VecGetArray(cvlocal,&cv_p);
1007: VecGetArray(cilocal,&ci_p);
1009: /* Create the mass matrix M_0 */
1010: MatGetLocalSize(M,&n,PETSC_NULL);
1011: DMDAGetInfo(user->da1,PETSC_NULL,&Xda,&Yda,&Zda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
1013: if (Xda!=Yda) {
1014: printf("Currently different Xda and Yda are not supported");
1015: }
1016: if (user->periodic) {
1017: h = (user->xmax-user->xmin)/Xda;
1018: } else {
1019: h = (user->xmax-user->xmin)/(Xda-1.0);
1020: }
1021:
1023: /* Get local element info */
1024: DMDAGetElements(user->da1,&nele,&nen,&ele);
1025:
1026: for(i=0;i<nele;i++) {
1027: for (j=0; j<8; j++)
1028: {
1029: idx[j] = ele[8*i + j];
1030: }
1031:
1032: PetscInt r,row,cols[8];
1033: PetscScalar vals[8];
1035: for (r=0;r<8;r++) {
1036: row = 5*idx[r];
1037: for(j = 0; j < 8; j++)
1038: {
1039: cols[j] = 5*idx[j]; vals[j] = 0.0;
1040:
1041: }
1042:
1043: /* Insert values in matrix M for 1st dof */
1044: MatSetValuesLocal(M,1,&row,8,cols,vals,INSERT_VALUES);
1045:
1046: row = 5*idx[r]+2;
1047: for(j = 0; j < 8; j++)
1048: {
1049: cols[j] = 5*idx[j]+2; vals[j] = 0.0;
1050: }
1051:
1053: /* Insert values in matrix M for 3nd dof */
1054: MatSetValuesLocal(M,1,&row,8,cols,vals,INSERT_VALUES);
1055: }
1056: }
1057:
1058: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1059: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
1061: /** blocks of M_2 **/
1062:
1063: for (i = 0; i<8; i++)
1064: {
1065: eM_2[i][i] = 12;
1066: }
1067:
1068: 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;
1069: 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;
1070: 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;
1071: eM_2[3][4] = -3;eM_2[3][5] = -3; eM_2[3][6] = -3; eM_2[3][7] = 0;
1072: eM_2[4][5] = 0; eM_2[4][6] = -3; eM_2[4][7] = 0;
1073: eM_2[5][6] = 0; eM_2[5][7] = -3;
1074: eM_2[6][7] = 0;
1075:
1076: for (i = 0; i<8; i++)
1077: {
1078: for (j =0; j<8; j++)
1079: {
1080: if (i>j)
1081: {
1082: eM_2[i][j] = eM_2[j][i];
1083: }
1084: }
1085: }
1086:
1087:
1088: for (i = 0; i<8; i++)
1089: {
1090: for (j =0; j<8; j++)
1091: {
1092: eM_2[i][j] = h*h*h/36.0 * eM_2[i][j];
1093: }
1094: }
1095:
1096: for(i=0;i<nele;i++) {
1097: for (j=0; j<8; j++)
1098: {
1099: idx[j] = ele[8*i + j];
1100: }
1101:
1102: PetscInt row,cols[8],r;
1103: PetscScalar vals[8];
1105: for(r=0;r<8;r++) {
1107: if (user->degenerate) {
1108: cv_sum = (2.0*user->smallnumber + cv_p[idx[0]] + cv_p[idx[1]])*user->Dv/(2.0*user->kBT);
1109: ci_sum = (2.0*user->smallnumber + ci_p[idx[0]] + ci_p[idx[1]])*user->Di/(2.0*user->kBT);
1110: } else {
1111: cv_sum = user->initv*user->Dv/user->kBT;
1112: ci_sum = user->initv*user->Di/user->kBT;
1113: }
1115: row = 5*idx[r];
1116: for (j=0; j<8; j++)
1117: {
1118: cols[j] = 5*idx[j]; vals[j] = user->dt*eM_2[r][j]*cv_sum;
1119: }
1120: /* Insert values in matrix M for 1st dof */
1121: MatSetValuesLocal(M,1,&row,8,cols,vals,ADD_VALUES);
1123: row = 5*idx[r]+2;
1124: for (j=0; j<8; j++)
1125: {
1126: cols[j] = 5*idx[j]+2; vals[j] = user->dt*eM_2[r][j]*ci_sum;
1127: }
1128:
1129: /* Insert values in matrix M for 3nd dof */
1130: MatSetValuesLocal(M,1,&row,8,cols,vals,ADD_VALUES);
1131: }
1132: }
1133:
1134:
1135: VecRestoreArray(cvlocal,&cv_p);
1136: VecRestoreArray(cilocal,&ci_p);
1137: DMRestoreLocalVector(user->da2,&cvlocal);
1138: DMRestoreLocalVector(user->da2,&cilocal);
1139:
1140: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1141: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
1142: DMDARestoreElements(user->da1,&nele,&nen,&ele);
1144: return(0);
1145: }