Actual source code: ex61.c
1: static char help[] = "2D coupled Allen-Cahn and Cahn-Hilliard equation for constant mobility and triangular elements. Use periodic boundary condidtions.\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: ./ex61 -ksp_type gmres -snes_vi_monitor -snes_atol 1.e-11 -da_refine 3 -T 0.1 -ksp_monitor_true_residual -pc_type lu -pc_factor_mat_solver_package superlu -snes_converged_reason -ksp_converged_reason -ksp_rtol 1.e-9 -snes_ls_monitor -VG 10 -draw_fields 1,3,4 -snes_ls basic
14: ./ex61 -ksp_type gmres -snes_vi_monitor -snes_atol 1.e-11 -da_refine 4 -T 0.1 -ksp_monitor_true_residual -pc_type sor -snes_converged_reason -ksp_converged_reason -ksp_rtol 1.e-9 -snes_ls_monitor -VG 10 -draw_fields 1,3,4 -snes_ls basic
16: ./ex61 -ksp_type fgmres -snes_vi_monitor -snes_atol 1.e-11 -da_refine 5 -T 0.1 -ksp_monitor_true_residual -snes_converged_reason -ksp_converged_reason -ksp_rtol 1.e-9 -snes_ls_monitor -VG 10 -draw_fields 1,3,4 -snes_ls basic -pc_type mg -pc_mg_galerkin
18: ./ex61 -ksp_type fgmres -snes_vi_monitor -snes_atol 1.e-11 -da_refine 5 -snes_converged_reason -ksp_converged_reason -snes_ls_monitor -VG 1 -draw_fields 1,3,4 -pc_type mg -pc_mg_galerkin -log_summary -dt .0000000000001 -mg_coarse_pc_type svd -ksp_monitor_true_residual -ksp_rtol 1.e-9
20: Movie version
21: ./ex61 -ksp_type fgmres -snes_vi_monitor -snes_atol 1.e-11 -da_refine 6 -snes_converged_reason -ksp_converged_reason -snes_ls_monitor -VG 10 -draw_fields 1,3,4 -pc_type mg -pc_mg_galerkin -log_summary -dt .000001 -mg_coarse_pc_type redundant -mg_coarse_redundant_pc_type svd -ksp_monitor_true_residual -ksp_rtol 1.e-9 -snes_ls basic -T .0020
23: */
25: /*
26: Possible additions to the code. At each iteration count the number of solution elements that are at the upper bound and stop the program if large
28: Is the solution at time 0 nonsense? Looks totally different from first time step. Why does cubic line search at beginning screw it up?
30: Add command-line option for constant or degenerate mobility
31: Add command-line option for graphics at each time step
33: Check time-step business; what should it be? How to check that it is good?
34: Make random right hand side forcing function proportional to time step so smaller time steps don't mean more radiation
35: How does the multigrid linear solver work now?
36: What happens when running with degenerate mobility
38:
39: */
41: #include petscsnes.h
42: #include petscdmda.h
44: typedef struct{
45: PetscReal dt,T; /* Time step and end time */
46: PetscReal dtevent; /* time scale of radiation events, roughly one event per dtevent */
47: PetscInt maxevents; /* once this number of events is reached no more events are generated */
48: PetscReal initv; /* initial value of phase variables */
49: PetscReal initeta;
50: PetscBool degenerate; /* use degenerate mobility */
51: PetscReal smallnumber;
52: PetscBool graphics;
53: PetscBool twodomain;
54: PetscBool radiation;
55: PetscBool voidgrowth; /* use initial conditions for void growth */
56: DM da1,da2;
57: Mat M; /* Jacobian matrix */
58: Mat M_0;
59: Vec q,wv,cv,wi,ci,eta,cvi,DPsiv,DPsii,DPsieta,Pv,Pi,Piv,logcv,logci,logcvi,Riv;
60: Vec phi1,phi2,Phi2D_V,Sv,Si; /* for twodomain modeling */
61: Vec work1,work2,work3,work4;
62: PetscScalar Dv,Di,Evf,Eif,A,kBT,kav,kai,kaeta,Rsurf,Rbulk,L,VG; /* physics parameters */
63: PetscScalar Svr,Sir,cv_eq,ci_eq; /* for twodomain modeling */
64: PetscReal asmallnumber; /* gets added to degenerate mobility */
65: PetscReal xmin,xmax,ymin,ymax;
66: PetscInt Mda, Nda;
67: PetscViewer graphicsfile; /* output of solution at each times step */
68: }AppCtx;
70: PetscErrorCode GetParams(AppCtx*);
71: PetscErrorCode SetRandomVectors(AppCtx*,PetscReal);
72: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
73: PetscErrorCode SetUpMatrices(AppCtx*);
74: PetscErrorCode UpdateMatrices(AppCtx*);
75: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
76: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
77: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
78: PetscErrorCode Update_q(AppCtx*);
79: PetscErrorCode Update_u(Vec,AppCtx*);
80: PetscErrorCode DPsi(AppCtx*);
81: PetscErrorCode LaplacianFiniteDifference(AppCtx*);
82: PetscErrorCode Llog(Vec,Vec);
83: PetscErrorCode Phi(AppCtx*);
87: int main(int argc, char **argv)
88: {
89: PetscErrorCode ierr;
90: Vec x,r; /* Solution and residual vectors */
91: SNES snes; /* Nonlinear solver context */
92: AppCtx user; /* Application context */
93: Vec xl,xu; /* Upper and lower bounds on variables */
94: Mat J;
95: PetscScalar t=0.0,normq;
96: /* PetscViewer view_out, view_q, view_psi, view_mat;*/
97: /* PetscViewer view_rand;*/
98: IS inactiveconstraints;
99: PetscInt ninactiveconstraints,N;
100: SNESConvergedReason reason;
101: /*PetscViewer view_out, view_cv,view_eta,view_vtk_cv,view_vtk_eta;*/
102: char cv_filename[80],eta_filename[80];
103: /*PetscReal bounds[] = {1000.0,-1000.,0.0,1.0,1000.0,-1000.0,0.0,1.0,1000.0,-1000.0}; */
105: PetscInitialize(&argc,&argv, (char*)0, help);
106:
107: /* Get physics and time parameters */
108: GetParams(&user);
109: /* Create a 1D DA with dof = 5; the whole thing */
110: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX, -3,-3,PETSC_DECIDE,PETSC_DECIDE, 5, 1,PETSC_NULL,PETSC_NULL,&user.da1);
111:
112: /* Create a 1D DA with dof = 1; for individual componentes */
113: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX, -3,-3,PETSC_DECIDE,PETSC_DECIDE, 1, 1,PETSC_NULL,PETSC_NULL,&user.da2);
116: /* Set Element type (triangular) */
117: DMDASetElementType(user.da1,DMDA_ELEMENT_P1);
118: DMDASetElementType(user.da2,DMDA_ELEMENT_P1);
119:
120: /* Set x and y coordinates */
121: DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,user.ymin,user.ymax,PETSC_NULL,PETSC_NULL);
122: DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,user.ymin,user.ymax,PETSC_NULL,PETSC_NULL);
125: /* Get global vector x from DM (da1) and duplicate vectors r,xl,xu */
126: DMCreateGlobalVector(user.da1,&x);
127: VecGetSize(x,&N);
128: VecDuplicate(x,&r);
129: VecDuplicate(x,&xl);
130: VecDuplicate(x,&xu);
131: VecDuplicate(x,&user.q);
132:
133: /* Get global vector user->wv from da2 and duplicate other vectors */
134: DMCreateGlobalVector(user.da2,&user.wv);
135: VecDuplicate(user.wv,&user.cv);
136: VecDuplicate(user.wv,&user.wi);
137: VecDuplicate(user.wv,&user.ci);
138: VecDuplicate(user.wv,&user.eta);
139: VecDuplicate(user.wv,&user.cvi);
140: VecDuplicate(user.wv,&user.DPsiv);
141: VecDuplicate(user.wv,&user.DPsii);
142: VecDuplicate(user.wv,&user.DPsieta);
143: VecDuplicate(user.wv,&user.Pv);
144: VecDuplicate(user.wv,&user.Pi);
145: VecDuplicate(user.wv,&user.Piv);
146: VecDuplicate(user.wv,&user.logcv);
147: VecDuplicate(user.wv,&user.logci);
148: VecDuplicate(user.wv,&user.logcvi);
149: VecDuplicate(user.wv,&user.work1);
150: VecDuplicate(user.wv,&user.work2);
151: VecDuplicate(user.wv,&user.Riv);
152: VecDuplicate(user.wv,&user.phi1);
153: VecDuplicate(user.wv,&user.phi2);
154: VecDuplicate(user.wv,&user.Phi2D_V);
155: VecDuplicate(user.wv,&user.Sv);
156: VecDuplicate(user.wv,&user.Si);
157: VecDuplicate(user.wv,&user.work3);
158: VecDuplicate(user.wv,&user.work4);
159: /* for twodomain modeling */
160: VecDuplicate(user.wv,&user.phi1);
161: VecDuplicate(user.wv,&user.phi2);
162: VecDuplicate(user.wv,&user.Phi2D_V);
163: VecDuplicate(user.wv,&user.Sv);
164: VecDuplicate(user.wv,&user.Si);
165: VecDuplicate(user.wv,&user.work3);
166: VecDuplicate(user.wv,&user.work4);
169: /* Get Jacobian matrix structure from the da for the entire thing, da1 */
170: DMGetMatrix(user.da1,MATAIJ,&user.M);
171: /* Get the (usual) mass matrix structure from da2 */
172: DMGetMatrix(user.da2,MATAIJ,&user.M_0);
173: SetInitialGuess(x,&user);
174: /* twodomain modeling */
175: if (user.twodomain) {
176: Phi(&user);
177: }
179: /* Form the jacobian matrix and M_0 */
180: SetUpMatrices(&user);
181: MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
182:
183: SNESCreate(PETSC_COMM_WORLD,&snes);
184: SNESSetDM(snes,user.da1);
186: SNESSetFunction(snes,r,FormFunction,(void*)&user);
187: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
188:
190: SNESSetType(snes,SNESVI);
191: SNESSetTolerances(snes,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,100,PETSC_DEFAULT);
193: SNESSetFromOptions(snes);
194: SetVariableBounds(user.da1,xl,xu);
195: SNESVISetVariableBounds(snes,xl,xu);
196:
197: /* PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_rand",FILE_MODE_WRITE,&view_rand);
198: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat2",FILE_MODE_WRITE,&view_mat);
199: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
200: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
201: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);*/
202:
203: /*PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
204:
205: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_cv",FILE_MODE_WRITE,&view_cv);
206: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_eta",FILE_MODE_WRITE,&view_eta);*/
207:
208: /* PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds); */
209: if (user.graphics) {
210: VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
211: }
212: if (user.graphicsfile) {
213: DMView(user.da1,user.graphicsfile);
214: VecView(x,user.graphicsfile);
215: }
216: while (t<user.T) {
217: SNESSetFunction(snes,r,FormFunction,(void*)&user);
218: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
220: SetRandomVectors(&user,t);
221: /* VecView(user.Pv,view_rand);
222: VecView(user.Pi,view_rand);
223: VecView(user.Piv,view_rand);*/
225: DPsi(&user);
226: /* VecView(user.DPsiv,view_psi);
227: VecView(user.DPsii,view_psi);
228: VecView(user.DPsieta,view_psi);*/
230: Update_q(&user);
232: /* VecView(user.q,view_q);*/
233: /* MatView(user.M,view_mat);*/
235:
236:
237: sprintf(cv_filename,"file_cv_%f.vtk",t);
238: sprintf(eta_filename,"file_eta_%f.vtk",t);
239: /* PetscViewerASCIIOpen(PETSC_COMM_WORLD,cv_filename,&view_vtk_cv);
240: PetscViewerASCIIOpen(PETSC_COMM_WORLD,eta_filename,&view_vtk_eta);
241: PetscViewerSetFormat(view_vtk_cv, PETSC_VIEWER_ASCII_VTK);
242: PetscViewerSetFormat(view_vtk_eta, PETSC_VIEWER_ASCII_VTK);
243: DMView(user.da2,view_vtk_cv);
244: DMView(user.da2,view_vtk_eta);
245: VecView(user.cv,view_cv);
246: VecView(user.eta,view_eta);
247: VecView(user.cv,view_vtk_cv);
248: VecView(user.eta,view_vtk_eta);
249: PetscViewerDestroy(&view_vtk_cv);
250: PetscViewerDestroy(&view_vtk_eta);*/
252:
253: VecNorm(user.q,NORM_2,&normq);
254: printf("2-norm of q = %14.12f\n",normq);
255: SNESSolve(snes,PETSC_NULL,x);
256: SNESGetConvergedReason(snes,&reason);
257: if (reason < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Nonlinear solver failed");
258: SNESVIGetInactiveSet(snes,&inactiveconstraints);
259: ISGetSize(inactiveconstraints,&ninactiveconstraints);
260: /* if (ninactiveconstraints < .90*N) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP,"To many active constraints, model has become non-physical"); */
262: /* VecView(x,view_out);*/
263: if (user.graphics) {
264: VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
265: }
266: /* VecView(x,PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD));*/
267: PetscInt its;
268: SNESGetIterationNumber(snes,&its);
269: PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %g in %d iterations\n",t,its);
271: Update_u(x,&user);
272: UpdateMatrices(&user);
273: t = t + user.dt;
274: if (user.graphicsfile) {
275: VecView(x,user.graphicsfile);
276: }
277: }
278:
279: /* PetscViewerDestroy(&view_rand);
280: PetscViewerDestroy(&view_mat);
281: PetscViewerDestroy(&view_q);
282: PetscViewerDestroy(&view_out);
283: PetscViewerDestroy(&view_psi);
284: PetscViewerDestroy(&view_out);
285:
286: PetscViewerDestroy(&view_cv);
287: PetscViewerDestroy(&view_eta);*/
288:
289: if (user.graphicsfile) {
290: PetscViewerDestroy(&user.graphicsfile);
291: }
292:
293: VecDestroy(&x);
294: VecDestroy(&r);
295: VecDestroy(&xl);
296: VecDestroy(&xu);
297: VecDestroy(&user.q);
298: VecDestroy(&user.wv);
299: VecDestroy(&user.cv);
300: VecDestroy(&user.wi);
301: VecDestroy(&user.ci);
302: VecDestroy(&user.eta);
303: VecDestroy(&user.cvi);
304: VecDestroy(&user.DPsiv);
305: VecDestroy(&user.DPsii);
306: VecDestroy(&user.DPsieta);
307: VecDestroy(&user.Pv);
308: VecDestroy(&user.Pi);
309: VecDestroy(&user.Piv);
310: VecDestroy(&user.logcv);
311: VecDestroy(&user.logci);
312: VecDestroy(&user.logcvi);
313: VecDestroy(&user.work1);
314: VecDestroy(&user.work2);
315: VecDestroy(&user.Riv);
316: MatDestroy(&user.M);
317: MatDestroy(&user.M_0);
318: DMDestroy(&user.da1);
319: DMDestroy(&user.da2);
321:
322: PetscFinalize();
323: return 0;
324: }
328: PetscErrorCode Update_u(Vec X,AppCtx *user)
329: {
331: PetscInt i,n;
332: PetscScalar *xx,*wv_p,*cv_p,*wi_p,*ci_p,*eta_p;
333:
335: VecGetLocalSize(user->wv,&n);
336: VecGetArray(X,&xx);
337: VecGetArray(user->wv,&wv_p);
338: VecGetArray(user->cv,&cv_p);
339: VecGetArray(user->wi,&wi_p);
340: VecGetArray(user->ci,&ci_p);
341: VecGetArray(user->eta,&eta_p);
344: for(i=0;i<n;i++) {
345: wv_p[i] = xx[5*i];
346: cv_p[i] = xx[5*i+1];
347: wi_p[i] = xx[5*i+2];
348: ci_p[i] = xx[5*i+3];
349: eta_p[i] = xx[5*i+4];
350: }
351: VecRestoreArray(X,&xx);
352: VecRestoreArray(user->wv,&wv_p);
353: VecRestoreArray(user->cv,&cv_p);
354: VecRestoreArray(user->wi,&wi_p);
355: VecRestoreArray(user->ci,&ci_p);
356: VecRestoreArray(user->eta,&eta_p);
357:
358: return(0);
359: }
363: PetscErrorCode Update_q(AppCtx *user)
364: {
366: PetscScalar *q_p,*w1,*w2,max1;
367: PetscInt i,n;
369:
371:
372: VecPointwiseMult(user->Riv,user->eta,user->eta);
373: VecScale(user->Riv,user->Rsurf);
374: VecShift(user->Riv,user->Rbulk);
375: VecPointwiseMult(user->Riv,user->ci,user->Riv);
376: VecPointwiseMult(user->Riv,user->cv,user->Riv);
377:
378: VecCopy(user->Riv,user->work1);
379: VecScale(user->work1,user->dt);
380: VecAXPY(user->work1,-1.0,user->cv);
381: MatMult(user->M_0,user->work1,user->work2);
382:
383: VecGetArray(user->q,&q_p);
384: VecGetArray(user->work1,&w1);
385: VecGetArray(user->work2,&w2);
386: VecGetLocalSize(user->wv,&n);
387: for (i=0;i<n;i++) {
388: q_p[5*i]=w2[i];
389: }
391: MatMult(user->M_0,user->DPsiv,user->work1);
392: for (i=0;i<n;i++) {
393: q_p[5*i+1]=w1[i];
394: }
396: VecCopy(user->Riv,user->work1);
397: VecScale(user->work1,user->dt);
398: VecAXPY(user->work1,-1.0,user->ci);
399: MatMult(user->M_0,user->work1,user->work2);
400: for (i=0;i<n;i++) {
401: q_p[5*i+2]=w2[i];
402: }
404: MatMult(user->M_0,user->DPsii,user->work1);
405: for (i=0;i<n;i++) {
406: q_p[5*i+3]=w1[i];
407: }
409: VecCopy(user->DPsieta,user->work1);
410: VecScale(user->work1,user->L);
411: VecAXPY(user->work1,-1.0/user->dt,user->eta);
412: MatMult(user->M_0,user->work1,user->work2);
413: for (i=0;i<n;i++) {
414: q_p[5*i+4]=w2[i];
415: }
416:
417: VecRestoreArray(user->q,&q_p);
418: VecRestoreArray(user->work1,&w1);
419: VecRestoreArray(user->work2,&w2);
420:
421: return(0);
422: }
426: PetscErrorCode DPsi(AppCtx* user)
427: {
428: PetscErrorCode ierr;
429: PetscScalar Evf=user->Evf,Eif=user->Eif,kBT=user->kBT,A=user->A;
430: PetscScalar *cv_p,*ci_p,*eta_p,*logcv_p,*logci_p,*logcvi_p,*DPsiv_p,*DPsii_p,*DPsieta_p;
431: PetscInt n,i;
435: VecGetLocalSize(user->cv,&n);
436: VecGetArray(user->cv,&cv_p);
437: VecGetArray(user->ci,&ci_p);
438: VecGetArray(user->eta,&eta_p);
439: VecGetArray(user->logcv,&logcv_p);
440: VecGetArray(user->logci,&logci_p);
441: VecGetArray(user->logcvi,&logcvi_p);
442: VecGetArray(user->DPsiv,&DPsiv_p);
443: VecGetArray(user->DPsii,&DPsii_p);
444: VecGetArray(user->DPsieta,&DPsieta_p);
445:
446: Llog(user->cv,user->logcv);
447:
448: Llog(user->ci,user->logci);
449:
451: VecCopy(user->cv,user->cvi);
452: VecAXPY(user->cvi,1.0,user->ci);
453: VecScale(user->cvi,-1.0);
454: VecShift(user->cvi,1.0);
455: Llog(user->cvi,user->logcvi);
456:
457: for (i=0;i<n;i++)
458: {
459: 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);
461: 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] ;
462:
463: 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]);
464:
465:
466: }
468: VecRestoreArray(user->cv,&cv_p);
469: VecRestoreArray(user->ci,&ci_p);
470: VecRestoreArray(user->eta,&eta_p);
471: VecRestoreArray(user->logcv,&logcv_p);
472: VecRestoreArray(user->logci,&logci_p);
473: VecRestoreArray(user->logcvi,&logcvi_p);
474: VecRestoreArray(user->DPsiv,&DPsiv_p);
475: VecRestoreArray(user->DPsii,&DPsii_p);
476: VecRestoreArray(user->DPsieta,&DPsieta_p);
479: return(0);
481: }
486: PetscErrorCode Llog(Vec X, Vec Y)
487: {
488: PetscErrorCode ierr;
489: PetscScalar *x,*y;
490: PetscInt n,i;
493:
494: VecGetArray(X,&x);
495: VecGetArray(Y,&y);
496: VecGetLocalSize(X,&n);
497: for (i=0;i<n;i++) {
498: if (x[i] < 1.0e-12) {
499: y[i] = log(1.0e-12);
500: }
501: else {
502: y[i] = log(x[i]);
503: }
504: }
505:
506: return(0);
507: }
512: PetscErrorCode SetInitialGuess(Vec X,AppCtx* user)
513: {
514: PetscErrorCode ierr;
515: PetscInt n,i,Mda,Nda;
516: PetscScalar *xx,*cv_p,*ci_p,*wv_p,*wi_p,*eta;
518: /* needed for the void growth case */
519: PetscScalar xmid,ymid,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;
520: PetscInt nele,nen,idx[3];
521: const PetscInt *ele;
522: PetscScalar x[3],y[3];
523: Vec coords;
524: const PetscScalar *_coords;
525: PetscViewer view;
526: PetscScalar xwidth = user->xmax - user->xmin;
530: VecGetLocalSize(X,&n);
532: if (user->voidgrowth) {
533: DMDAGetInfo(user->da2,PETSC_NULL,&Mda,&Nda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
534: DMDAGetGhostedCoordinates(user->da2,&coords);
535: VecGetArrayRead(coords,&_coords);
537: h = (user->xmax-user->xmin)/Mda;
538: xmid = (user->xmax + user->xmin)/2.0;
539: ymid = (user->ymax + user->ymin)/2.0;
540: lambda = 4.0*h;
542: DMDAGetElements(user->da2,&nele,&nen,&ele);
543: for (i=0;i < nele; i++) {
544: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
545:
546: x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
547: x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
548: x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
549:
550: PetscInt k;
551: PetscScalar vals_cv[3],vals_ci[3],vals_eta[3],s,hhr,r;
552: for (k=0; k < 3 ; k++) {
553: s = sqrt((x[k] - xmid)*(x[k] - xmid) + (y[k] - ymid)*(y[k] - ymid));
554: if (s < xwidth*(5.0/64.0)) {
555: vals_cv[k] = cv_v;
556: vals_ci[k] = ci_v;
557: vals_eta[k] = eta_v;
558: } else if (s>= xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0) ) {
559: //r = (s - xwidth*(6.0/64.0) )/(0.5*lambda);
560: r = (s - xwidth*(6.0/64.0) )/(xwidth/64.0);
561: hhr = 0.25*(-r*r*r + 3*r + 2);
562: vals_cv[k] = cv_m + (1.0 - hhr)*(cv_v - cv_m);
563: vals_ci[k] = ci_m + (1.0 - hhr)*(ci_v - ci_m);
564: vals_eta[k] = eta_m + (1.0 - hhr)*(eta_v - eta_m);
565: } else {
566: vals_cv[k] = cv_m;
567: vals_ci[k] = ci_m;
568: vals_eta[k] = eta_m;
569: }
570: }
571: VecSetValuesLocal(user->cv,3,idx,vals_cv,INSERT_VALUES);
572: VecSetValuesLocal(user->ci,3,idx,vals_ci,INSERT_VALUES);
573: VecSetValuesLocal(user->eta,3,idx,vals_eta,INSERT_VALUES);
574: }
575: VecAssemblyBegin(user->cv);
576: VecAssemblyEnd(user->cv);
577: VecAssemblyBegin(user->ci);
578: VecAssemblyEnd(user->ci);
579: VecAssemblyBegin(user->eta);
580: VecAssemblyEnd(user->eta);
582: DMDARestoreElements(user->da2,&nele,&nen,&ele);
583: VecRestoreArrayRead(coords,&_coords);
585: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view);
586: VecView(user->cv,view);
587: VecView(user->ci,view);
588: VecView(user->eta,view);
589: PetscViewerDestroy(&view);
590: }
591: else {
592: VecSet(user->cv,user->initv);
593: VecSet(user->ci,user->initv);
594: VecSet(user->eta,user->initeta);
595: }
597: DPsi(user);
598: VecCopy(user->DPsiv,user->wv);
599: VecCopy(user->DPsii,user->wi);
601: VecGetArray(X,&xx);
602: VecGetArray(user->cv,&cv_p);
603: VecGetArray(user->ci,&ci_p);
604: VecGetArray(user->wv,&wv_p);
605: VecGetArray(user->wi,&wi_p);
606: VecGetArray(user->eta,&eta);
607: for (i=0;i<n/5;i++)
608: {
609: xx[5*i]=wv_p[i];
610: xx[5*i+1]=cv_p[i];
611: xx[5*i+2]=wi_p[i];
612: xx[5*i+3]=ci_p[i];
613: xx[5*i+4]=eta[i];
614: }
616: /* VecView(user->wv,view);
617: VecView(user->cv,view);
618: VecView(user->wi,view);
619: VecView(user->ci,view);
620: PetscViewerDestroy(&view);*/
622: VecRestoreArray(X,&xx);
623: VecRestoreArray(user->cv,&cv_p);
624: VecRestoreArray(user->ci,&ci_p);
625: VecRestoreArray(user->wv,&wv_p);
626: VecRestoreArray(user->wi,&wi_p);
627: VecRestoreArray(user->eta,&eta);
628:
629: return(0);
630: }
632: typedef struct {
633: PetscReal dt,x,y,strength;
634: } RandomValues;
639: PetscErrorCode SetRandomVectors(AppCtx* user,PetscReal t)
640: {
641: PetscErrorCode ierr;
642: static RandomValues *randomvalues = 0;
643: static PetscInt randindex = 0,n; /* indicates how far into the randomvalues we have currently used */
644: static PetscReal randtime = 0; /* indicates time of last radiation event */
645: PetscInt i,j,M,N,cnt = 0;
646: PetscInt xs,ys,xm,ym;
649: if (!randomvalues) {
650: PetscViewer viewer;
651: char filename[PETSC_MAX_PATH_LEN];
652: PetscBool flg;
653: PetscInt seed;
655: PetscOptionsGetInt(PETSC_NULL,"-random_seed",&seed,&flg);
656: if (flg) {
657: sprintf(filename,"ex61.random.%d",(int)seed);
658: } else {
659: PetscStrcpy(filename,"ex61.random");
660: }
661: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
662: PetscViewerBinaryRead(viewer,&n,1,PETSC_INT);
663: PetscMalloc(n*sizeof(RandomValues),&randomvalues);
664: PetscViewerBinaryRead(viewer,randomvalues,4*n,PETSC_DOUBLE);
665: for (i=0; i<n; i++) randomvalues[i].dt = randomvalues[i].dt*user->dtevent;
666: PetscViewerDestroy(&viewer);
667: }
668: user->maxevents = PetscMin(user->maxevents,n);
670: VecSet(user->Pv,0.0);
671: DMDAGetInfo(user->da1,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
672: DMDAGetGhostCorners(user->da1,&xs,&ys,0,&xm,&ym,0);
673: while (user->maxevents > randindex && randtime + randomvalues[randindex].dt < t + user->dt) { /* radiation event has occured since last time step */
674: i = ((PetscInt) (randomvalues[randindex].x*M)) - xs;
675: j = ((PetscInt) (randomvalues[randindex].y*N)) - ys;
676: if (i >= 0 && i < xm && j >= 0 && j < ym) { /* point is on this process */
677: /* need to make sure eta at the given point is not great than .8 */
678: VecSetValueLocal(user->Pv,i + 1 + xm*(j + 1), randomvalues[randindex].strength*user->VG,INSERT_VALUES);
679: }
680: randtime += randomvalues[randindex++].dt;
681: cnt++;
682: }
683: PetscPrintf(PETSC_COMM_WORLD,"Number of radiation events %d\n",cnt);
684: VecAssemblyBegin(user->Pv);
685: VecAssemblyEnd(user->Pv);
687: VecCopy(user->Pv,user->Pi);
688: VecScale(user->Pi,0.9);
689: VecPointwiseMult(user->Piv,user->Pi,user->Pv);
690: return(0);
691:
692: }
696: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ctx)
697: {
699: AppCtx *user=(AppCtx*)ctx;
700:
702: MatMultAdd(user->M,X,user->q,F);
703: return(0);
704: }
708: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
709: {
711: AppCtx *user=(AppCtx*)ctx;
712:
714: *flg = SAME_NONZERO_PATTERN;
715: MatCopy(user->M,*J,*flg);
716: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
717: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
718: return(0);
719: }
722: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
723: {
725: PetscScalar ***l,***u;
726: PetscInt xs,xm,ys,ym;
727: PetscInt j,i;
728:
730: DMDAVecGetArrayDOF(da,xl,&l);
731: DMDAVecGetArrayDOF(da,xu,&u);
732:
733: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
734:
735: for (j=ys; j<ys+ym; j++) {
736: for(i=xs; i < xs+xm;i++) {
737: l[j][i][0] = -SNES_VI_INF;
738: l[j][i][1] = 0.0;
739: l[j][i][2] = -SNES_VI_INF;
740: l[j][i][3] = 0.0;
741: l[j][i][4] = 0.0;
742: u[j][i][0] = SNES_VI_INF;
743: u[j][i][1] = 1.0;
744: u[j][i][2] = SNES_VI_INF;
745: u[j][i][3] = 1.0;
746: u[j][i][4] = 1.0;
747: }
748: }
749:
750:
751: DMDAVecRestoreArrayDOF(da,xl,&l);
752: DMDAVecRestoreArrayDOF(da,xu,&u);
753: return(0);
754: }
758: PetscErrorCode GetParams(AppCtx* user)
759: {
761: PetscBool flg,graphicsfile = PETSC_FALSE;
762:
764:
765: /* Set default parameters */
766: user->xmin = 0.0; user->xmax = 128.0;
767: user->ymin = 0.0; user->ymax = 128.0;
768: user->Dv = 1.0;
769: user->Di = 4.0;
770: user->Evf = 0.8;
771: user->Eif = 1.2;
772: user->A = 1.0;
773: user->kBT = 0.11;
774: user->kav = 1.0;
775: user->kai = 1.0;
776: user->kaeta = 1.0;
777: user->Rsurf = 10.0;
778: user->Rbulk = 1.0;
779: user->VG = 100.0;
780: user->L = 10.0;
782: user->T = 1.0e-2;
783: user->dt = 1.0e-4;
784: user->initv = .00069;
785: user->initeta = 0.0;
786: user->degenerate = PETSC_FALSE;
787: user->maxevents = 1000;
788: user->graphics = PETSC_TRUE;
790: /* twodomain modeling */
791: user->twodomain = PETSC_FALSE;
792: user->Svr = 0.5;
793: user->Sir = 0.5;
794: user->cv_eq = 6.9e-4;
795: user->ci_eq = 6.9e-4;
796: /* void growth */
797: user->voidgrowth = PETSC_FALSE;
799: user->radiation = PETSC_FALSE;
801: /* degenerate mobility */
802: user->smallnumber = 1.0e-3;
803: PetscOptionsBegin(PETSC_COMM_WORLD,PETSC_NULL,"Coupled Cahn-Hillard/Allen-Cahn Equations","Phasefield");
804: PetscOptionsReal("-Dv","Vacancy Diffusivity\n","None",user->Dv,&user->Dv,&flg);
805: PetscOptionsReal("-Di","Interstitial Diffusivity\n","None",user->Di,&user->Di,&flg);
806: PetscOptionsReal("-Evf","Vacancy Formation Energy\n","None",user->Evf,&user->Evf,&flg);
807: PetscOptionsReal("-Eif","Interstitial Formation energy\n","None",user->Eif,&user->Eif,&flg);
808: PetscOptionsReal("-A","???","None",user->A,&user->A,&flg);
809: PetscOptionsReal("-kBT","Boltzmann's Constant times the Absolute Temperature","None",user->kBT,&user->kBT,&flg);
810: PetscOptionsReal("-kav","???","None",user->kav,&user->kav,&flg);
811: PetscOptionsReal("-kai","???","None",user->kai,&user->kai,&flg);
812: PetscOptionsReal("-kaeta","???","None",user->kaeta,&user->kaeta,&flg);
813: PetscOptionsReal("-Rsurf","???","None",user->Rsurf,&user->Rsurf,&flg);
814: PetscOptionsReal("-Rbulk","???","None",user->Rbulk,&user->Rbulk,&flg);
815: PetscOptionsReal("-VG","Maximum increase in vacancy (or interstitial) concentration due to a cascade event","None",user->VG,&user->VG,&flg);
816: PetscOptionsReal("-L","???","None",user->L,&user->L,&flg);
818: PetscOptionsReal("-initv","Initial solution of Cv and Ci","None",user->initv,&user->initv,&flg);
819: PetscOptionsReal("-initeta","Initial solution of Eta","None",user->initeta,&user->initeta,&flg);
820: PetscOptionsBool("-degenerate","Run with degenerate mobility\n","None",user->degenerate,&user->degenerate,&flg);
821: PetscOptionsReal("-smallnumber","Small number added to degenerate mobility\n","None",user->smallnumber,&user->smallnumber,&flg);
822: PetscOptionsBool("-twodomain","Run two domain model\n","None",user->twodomain,&user->twodomain,&flg);
823: PetscOptionsBool("-voidgrowth","Use initial conditions for void growth\n","None",user->voidgrowth,&user->voidgrowth,&flg);
824: PetscOptionsBool("-radiation","Use initial conditions for void growth\n","None",user->radiation,&user->radiation,&flg);
825: PetscOptionsReal("-xmin","Lower X coordinate of domain\n","None",user->xmin,&user->xmin,&flg);
826: PetscOptionsReal("-xmax","Upper X coordinate of domain\n","None",user->xmax,&user->xmax,&flg);
827: PetscOptionsReal("-T","Total runtime\n","None",user->T,&user->T,&flg);
828: PetscOptionsReal("-dt","Time step\n","None",user->dt,&user->dt,&flg);
829: user->dtevent = user->dt;
830: PetscOptionsReal("-dtevent","Average time between random events\n","None",user->dtevent,&user->dtevent,&flg);
831: PetscOptionsInt("-maxevents","Maximum random events allowed\n","None",user->maxevents,&user->maxevents,&flg);
833: PetscOptionsBool("-graphics","Contour plot solutions at each timestep\n","None",user->graphics,&user->graphics,&flg);
834: PetscOptionsBool("-graphicsfile","Save solution at each timestep\n","None",graphicsfile,&graphicsfile,&flg);
835: if (graphicsfile) {
836: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"ex61.data",FILE_MODE_WRITE,&user->graphicsfile);
837: }
838: PetscOptionsEnd();
839: return(0);
840: }
845: PetscErrorCode SetUpMatrices(AppCtx* user)
846: {
847: PetscErrorCode ierr;
848: PetscInt nele,nen,i,n;
849: const PetscInt *ele;
850: PetscScalar dt=user->dt,hx,hy;
851:
852: PetscInt idx[3];
853: PetscScalar eM_0[3][3],eM_2_even[3][3],eM_2_odd[3][3];
854: PetscScalar cv_sum, ci_sum;
855: Mat M=user->M;
856: Mat M_0=user->M_0;
857: PetscInt Mda=user->Mda, Nda=user->Nda;
858: PetscScalar *cv_p,*ci_p;
859: /* newly added */
860: Vec cvlocal,cilocal;
863:
864: /* MatSetOption(M,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
865: MatSetOption(M_0,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);*/
867: /* new stuff */
868: DMGetLocalVector(user->da2,&cvlocal);
869: DMGetLocalVector(user->da2,&cilocal);
870: DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
871: DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
872: DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
873: DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
874: /* old stuff */
875: /*
876: VecGetArray(user->cv,&cv_p);
877: VecGetArray(user->ci,&ci_p);
878: */
879: /* new stuff */
880: VecGetArray(cvlocal,&cv_p);
881: VecGetArray(cilocal,&ci_p);
883: MatGetLocalSize(M,&n,PETSC_NULL);
884: DMDAGetInfo(user->da1,PETSC_NULL,&Mda,&Nda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
885: hx = (user->xmax-user->xmin)/Mda;
886: hy = (user->ymax-user->ymin)/Nda;
888: eM_0[0][0]=eM_0[1][1]=eM_0[2][2]=hx*hy/12.0;
889: eM_0[0][1]=eM_0[0][2]=eM_0[1][0]=eM_0[1][2]=eM_0[2][0]=eM_0[2][1]=hx*hy/24.0;
890:
891: eM_2_odd[0][0] = 1.0;
892: eM_2_odd[1][1] = eM_2_odd[2][2] = 0.5;
893: eM_2_odd[0][1] = eM_2_odd[0][2] = eM_2_odd[1][0]= eM_2_odd[2][0] = -0.5;
894: eM_2_odd[1][2] = eM_2_odd[2][1] = 0.0;
896: eM_2_even[0][0] = 1.0;
897: eM_2_even[1][1] = eM_2_even[2][2] = 0.5;
898: eM_2_even[0][1] = eM_2_even[0][2] = eM_2_even[1][0]= eM_2_even[2][0] = -0.5;
899: eM_2_even[1][2] = eM_2_even[2][1] = 0.0;
901: /* eM_2_even[1][1] = 1.0;
902: eM_2_even[0][0] = eM_2_even[2][2] = 0.5;
903: eM_2_even[0][1] = eM_2_even[1][0] = eM_2_even[1][2] = eM_2_even[2][1] = -0.5;
904: eM_2_even[0][2] = eM_2_even[2][0] = 0.0;
905: */
907: // for(k=0;k < Mda*Nda*2;k++) {
908: DMDAGetElements(user->da1,&nele,&nen,&ele);
909: for (i=0; i < nele; i++) {
910: /*
911: idx[0] = connect[k*3];
912: idx[1] = connect[k*3+1];
913: idx[2] = connect[k*3+2];
914: */
915: idx[0] = ele[3*i];
916: idx[1] = ele[3*i+1];
917: idx[2] = ele[3*i+2];
919: PetscInt row,cols[6],r,row_M_0,cols3[3];
920: PetscScalar vals[6],vals_M_0[3],vals3[3];
921:
922: for(r=0;r<3;r++) {
923: //row_M_0 = connect[k*3+r];
924: row_M_0 = idx[r];
926: vals_M_0[0]=eM_0[r][0];
927: vals_M_0[1]=eM_0[r][1];
928: vals_M_0[2]=eM_0[r][2];
929:
930:
931: MatSetValuesLocal(M_0,1,&row_M_0,3,idx,vals_M_0,ADD_VALUES);
932:
933: if (user->degenerate) {
934: cv_sum = (cv_p[idx[0]] + cv_p[idx[1]] + cv_p[idx[2]])*user->Dv/(3.0*user->kBT);
935: ci_sum = (ci_p[idx[0]] + ci_p[idx[1]] + ci_p[idx[2]])*user->Di/(3.0*user->kBT);
936: } else {
937: cv_sum = user->initv*user->Dv/user->kBT;
938: ci_sum = user->initv*user->Di/user->kBT;
939: }
940:
942:
943: row = 5*idx[r];
944: cols[0] = 5*idx[0]; vals[0] = dt*eM_2_odd[r][0]*cv_sum;
945: cols[1] = 5*idx[1]; vals[1] = dt*eM_2_odd[r][1]*cv_sum;
946: cols[2] = 5*idx[2]; vals[2] = dt*eM_2_odd[r][2]*cv_sum;
947: cols[3] = 5*idx[0]+1; vals[3] = eM_0[r][0];
948: cols[4] = 5*idx[1]+1; vals[4] = eM_0[r][1];
949: cols[5] = 5*idx[2]+1; vals[5] = eM_0[r][2];
951:
952: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
953:
954:
955: row = 5*idx[r]+1;
956: cols[0] = 5*idx[0]; vals[0] = -1.0*eM_0[r][0];
957: cols[1] = 5*idx[1]; vals[1] = -1.0*eM_0[r][1];
958: cols[2] = 5*idx[2]; vals[2] = -1.0*eM_0[r][2];
959: cols[3] = 5*idx[0]+1; vals[3] = user->kav*eM_2_odd[r][0];
960: cols[4] = 5*idx[1]+1; vals[4] = user->kav*eM_2_odd[r][1];
961: cols[5] = 5*idx[2]+1; vals[5] = user->kav*eM_2_odd[r][2];
962:
963: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
964:
965: row = 5*idx[r]+2;
966: cols[0] = 5*idx[0]+2; vals[0] = dt*eM_2_odd[r][0]*ci_sum;
967: cols[1] = 5*idx[1]+2; vals[1] = dt*eM_2_odd[r][1]*ci_sum;
968: cols[2] = 5*idx[2]+2; vals[2] = dt*eM_2_odd[r][2]*ci_sum;
969: cols[3] = 5*idx[0]+3; vals[3] = eM_0[r][0];
970: cols[4] = 5*idx[1]+3; vals[4] = eM_0[r][1];
971: cols[5] = 5*idx[2]+3; vals[5] = eM_0[r][2];
972:
973: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
974:
975:
976: row = 5*idx[r]+3;
977: cols[0] = 5*idx[0]+2; vals[0] = -1.0*eM_0[r][0];
978: cols[1] = 5*idx[1]+2; vals[1] = -1.0*eM_0[r][1];
979: cols[2] = 5*idx[2]+2; vals[2] = -1.0*eM_0[r][2];
980: cols[3] = 5*idx[0]+3; vals[3] = user->kai*eM_2_odd[r][0];
981: cols[4] = 5*idx[1]+3; vals[4] = user->kai*eM_2_odd[r][1];
982: cols[5] = 5*idx[2]+3; vals[5] = user->kai*eM_2_odd[r][2];
983:
984: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
985:
986:
987: row = 5*idx[r]+4;
988: /*
989: cols3[0] = 5*idx[0]+4; vals3[0] = eM_0[r][0]/dt + user->L*user->kaeta*dt*eM_2_odd[r][0];
990: cols3[1] = 5*idx[1]+4; vals3[1] = eM_0[r][1]/dt + user->L*user->kaeta*dt*eM_2_odd[r][1];
991: cols3[2] = 5*idx[2]+4; vals3[2] = eM_0[r][2]/dt + user->L*user->kaeta*dt*eM_2_odd[r][2];
992: */
993: cols3[0] = 5*idx[0]+4; vals3[0] = (eM_0[r][0]/dt + user->L*user->kaeta*eM_2_odd[r][0]);
994: cols3[1] = 5*idx[1]+4; vals3[1] = (eM_0[r][1]/dt + user->L*user->kaeta*eM_2_odd[r][1]);
995: cols3[2] = 5*idx[2]+4; vals3[2] = (eM_0[r][2]/dt + user->L*user->kaeta*eM_2_odd[r][2]);
996:
997: MatSetValuesLocal(M,1,&row,3,cols3,vals3,ADD_VALUES);
999: }
1000: }
1002: /* new */
1003: VecRestoreArray(cvlocal,&cv_p);
1004: VecRestoreArray(cilocal,&ci_p);
1005: DMRestoreLocalVector(user->da2,&cvlocal);
1006: DMRestoreLocalVector(user->da2,&cilocal);
1007: /* old */
1008: /*
1009: VecRestoreArray(user->cv,&cv_p);
1010: VecRestoreArray(user->ci,&ci_p);
1011: */
1012: MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
1013: MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);
1014:
1015: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1016: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
1017:
1018: DMDARestoreElements(user->da1,&nele,&nen,&ele);
1019:
1021: return(0);
1022: }
1027: PetscErrorCode UpdateMatrices(AppCtx* user)
1028: {
1029: PetscErrorCode ierr;
1030: PetscInt i,n,Mda,Nda,nele,nen;
1031: const PetscInt *ele;
1032:
1033: PetscInt idx[3];
1034: PetscScalar eM_2_odd[3][3],eM_2_even[3][3],h,dt=user->dt;
1035: Mat M=user->M;
1036: PetscScalar *cv_p,*ci_p,cv_sum,ci_sum;
1037: /* newly added */
1038: Vec cvlocal,cilocal;
1041:
1042:
1043: MatGetLocalSize(M,&n,PETSC_NULL);
1044:
1045: /* new stuff */
1046: DMGetLocalVector(user->da2,&cvlocal);
1047: DMGetLocalVector(user->da2,&cilocal);
1048: DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
1049: DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
1050: DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
1051: DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
1052: /* new stuff */
1053: VecGetArray(cvlocal,&cv_p);
1054: VecGetArray(cilocal,&ci_p);
1055:
1056: /* old stuff */
1057: /*
1058: VecGetArray(user->cv,&cv_p);
1059: VecGetArray(user->ci,&ci_p);
1060: */
1061: DMDAGetInfo(user->da1,PETSC_NULL,&Mda,&Nda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
1063:
1064:
1065: h = (user->xmax-user->xmin)/Mda;
1067: DMDAGetElements(user->da1,&nele,&nen,&ele);
1069: for(i=0; i < nele; i++) {
1070: /*
1071: idx[0] = connect[k*3];
1072: idx[1] = connect[k*3+1];
1073: idx[2] = connect[k*3+2];
1074: */
1075: idx[0] = ele[3*i];
1076: idx[1] = ele[3*i+1];
1077: idx[2] = ele[3*i+2];
1079: PetscInt r,row,cols[3];
1080: PetscScalar vals[3];
1081: for (r=0;r<3;r++) {
1082: row = 5*idx[r];
1083: cols[0] = 5*idx[0]; vals[0] = 0.0;
1084: cols[1] = 5*idx[1]; vals[1] = 0.0;
1085: cols[2] = 5*idx[2]; vals[2] = 0.0;
1087: /* Insert values in matrix M for 1st dof */
1088: MatSetValuesLocal(M,1,&row,3,cols,vals,INSERT_VALUES);
1089:
1090: row = 5*idx[r]+2;
1091: cols[0] = 5*idx[0]+2; vals[0] = 0.0;
1092: cols[1] = 5*idx[1]+2; vals[1] = 0.0;
1093: cols[2] = 5*idx[2]+2; vals[2] = 0.0;
1095: /* Insert values in matrix M for 3nd dof */
1096: MatSetValuesLocal(M,1,&row,3,cols,vals,INSERT_VALUES);
1097: }
1098: }
1099:
1100: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1101: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
1103: eM_2_odd[0][0] = 1.0;
1104: eM_2_odd[1][1] = eM_2_odd[2][2] = 0.5;
1105: eM_2_odd[0][1] = eM_2_odd[0][2] = eM_2_odd[1][0]= eM_2_odd[2][0] = -0.5;
1106: eM_2_odd[1][2] = eM_2_odd[2][1] = 0.0;
1108: eM_2_even[0][0] = 1.0;
1109: eM_2_even[1][1] = eM_2_even[2][2] = 0.5;
1110: eM_2_even[0][1] = eM_2_even[0][2] = eM_2_even[1][0]= eM_2_even[2][0] = -0.5;
1111: eM_2_even[1][2] = eM_2_even[2][1] = 0.0;
1113: /*
1114: eM_2_even[1][1] = 1.0;
1115: eM_2_even[0][0] = eM_2_even[2][2] = 0.5;
1116: eM_2_even[0][1] = eM_2_even[1][0] = eM_2_even[1][2] = eM_2_even[2][1] = -0.5;
1117: eM_2_even[0][2] = eM_2_even[2][0] = 0.0;
1118: */
1119:
1120: /* Get local element info */
1121: //for(k=0;k < Mda*Nda*2;k++) {
1122: for (i=0; i < nele; i++) {
1123: /*
1124: idx[0] = connect[k*3];
1125: idx[1] = connect[k*3+1];
1126: idx[2] = connect[k*3+2];
1127: */
1128: idx[0] = ele[3*i];
1129: idx[1] = ele[3*i+1];
1130: idx[2] = ele[3*i+2];
1132: PetscInt row,cols[3],r;
1133: PetscScalar vals[3];
1134:
1135: for(r=0;r<3;r++) {
1136:
1137: if (user->degenerate) {
1138: printf("smallnumber = %14.12f\n",user->smallnumber);
1139: cv_sum = (user->smallnumber + cv_p[idx[0]] + cv_p[idx[1]] + cv_p[idx[2]])*user->Dv/(3.0*user->kBT);
1140: ci_sum = (user->smallnumber + ci_p[idx[0]] + ci_p[idx[1]] + ci_p[idx[2]])*user->Di/(3.0*user->kBT);
1141: } else {
1142: cv_sum = user->initv*user->Dv/(user->kBT);
1143: ci_sum = user->initv*user->Di/user->kBT;
1144: }
1146:
1147:
1148: row = 5*idx[r];
1149: cols[0] = 5*idx[0]; vals[0] = dt*eM_2_odd[r][0]*cv_sum;
1150: cols[1] = 5*idx[1]; vals[1] = dt*eM_2_odd[r][1]*cv_sum;
1151: cols[2] = 5*idx[2]; vals[2] = dt*eM_2_odd[r][2]*cv_sum;
1152:
1153: /* Insert values in matrix M for 1st dof */
1154: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
1155:
1156:
1157: row = 5*idx[r]+2;
1158: cols[0] = 5*idx[0]+2; vals[0] = dt*eM_2_odd[r][0]*ci_sum;
1159: cols[1] = 5*idx[1]+2; vals[1] = dt*eM_2_odd[r][1]*ci_sum;
1160: cols[2] = 5*idx[2]+2; vals[2] = dt*eM_2_odd[r][2]*ci_sum;
1161:
1162: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
1164:
1165:
1166: }
1167:
1168: }
1170: /* new stuff */
1171: VecRestoreArray(cvlocal,&cv_p);
1172: VecRestoreArray(cilocal,&ci_p);
1173: DMRestoreLocalVector(user->da2,&cvlocal);
1174: DMRestoreLocalVector(user->da2,&cilocal);
1175: /* old stuff */
1176: /*
1177: VecRestoreArray(user->cv,&cv_p);
1178: VecRestoreArray(user->ci,&ci_p);
1179: */
1181: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1182: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
1185: return(0);
1186: }
1192: PetscErrorCode Phi(AppCtx* user)
1193: {
1194: PetscErrorCode ierr;
1195: PetscScalar xmid, xqu, lambda, h,x[3],y[3];
1196: Vec coords;
1197: const PetscScalar *_coords;
1198: PetscInt nele,nen,i,idx[3],Mda,Nda;
1199: const PetscInt *ele;
1200: PetscViewer view;
1204: DMDAGetInfo(user->da1,PETSC_NULL,&Mda,&Nda,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
1205: DMDAGetGhostedCoordinates(user->da2,&coords);
1206: VecGetArrayRead(coords,&_coords);
1208: h = (user->xmax - user->xmin)/Mda;
1209: xmid = (user->xmin + user->xmax)/2.0;
1210: xqu = (user->xmin + xmid)/2.0;
1211: lambda = 4.0*h;
1214: DMDAGetElements(user->da2,&nele,&nen,&ele);
1215: for (i=0;i < nele; i++) {
1216: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
1217: //printf("idx[0]=%d,idx[1]=%d,idx[2]=%d\n",idx[0],idx[1],idx[2]);
1219: x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
1220: x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
1221: x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
1223: //printf("x[0]=%f,x[1]=%f,x[2]=%f\n",x[0],x[1],x[2]);
1224: //printf("y[0]=%f,y[1]=%f,y[2]=%f\n",y[0],y[1],y[2]);
1225:
1226: PetscScalar vals1[3],vals2[3],dist1,dist2,s1,r,hhr,xc1,xc2;
1227: PetscInt k;
1229: xc1 = user->xmin;
1230: xc2 = xmid;
1232: //VecSet(user->phi1,0.0);
1233: for (k=0;k < 3; k++) {
1234: if (x[k]-xqu > 0) {
1235: s1 = (x[k] - xqu);
1236: } else {
1237: s1 = -(x[k] - xqu);
1238: }
1239: if (x[k] - xc1 > 0) {
1240: dist1 = (x[k] - xc1);
1241: } else {
1242: dist1 = -(x[k] - xc1);
1243: }
1244: if (x[k] - xc2 > 0) {
1245: dist2 = (x[k] - xc2);
1246: } else {
1247: dist2 = -(x[k] - xc2);
1248: }
1249: if (dist1 <= 0.5*lambda) {
1250: r = (x[k]-xc1)/(0.5*lambda);
1251: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
1252: vals1[k] = hhr;
1253: }
1254: else if (dist2 <= 0.5*lambda) {
1255: r = (x[k]-xc2)/(0.5*lambda);
1256: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
1257: vals1[k] = 1.0 - hhr;
1258: }
1259: else if (s1 <= xqu - 2.0*h) {
1260: vals1[k] = 1.0;
1261: }
1262:
1263: //else if ( abs(x[k]-(user->xmax-h)) < 0.1*h ) {
1264: else if ( (user->xmax-h)-x[k] < 0.1*h ) {
1265: vals1[k] = .15625;
1266: }
1267: else {
1268: vals1[k] = 0.0;
1269: }
1270: }
1271:
1272: VecSetValuesLocal(user->phi1,3,idx,vals1,INSERT_VALUES);
1274: xc1 = xmid;
1275: xc2 = user->xmax;
1277: //VecSet(user->phi2,0.0);
1278: for (k=0;k < 3; k++) {
1279: /*
1280: s1 = abs(x[k] - (xqu+xmid));
1281: dist1 = abs(x[k] - xc1);
1282: dist2 = abs(x[k] - xc2);
1283: */
1284: if (x[k]-(xqu + xmid) > 0) {
1285: s1 = (x[k] - (xqu + xmid));
1286: } else {
1287: s1 = -(x[k] - (xqu + xmid));
1288: }
1289: if (x[k] - xc1 > 0) {
1290: dist1 = (x[k] - xc1);
1291: } else {
1292: dist1 = -(x[k] - xc1);
1293: }
1294: if (x[k] - xc2 > 0) {
1295: dist2 = (x[k] - xc2);
1296: } else {
1297: dist2 = -(x[k] - xc2);
1298: }
1299:
1300: if (dist1 <= 0.5*lambda) {
1301: r = (x[k]-xc1)/(0.5*lambda);
1302: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
1303: vals2[k] = hhr;
1304: }
1305: else if (dist2 <= 0.5*lambda) {
1306: r = -(x[k]-xc2)/(0.5*lambda);
1307: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
1308: vals2[k] = hhr;
1309: }
1310: else if (s1 <= xqu - 2.0*h) {
1311: vals2[k] = 1.0;
1312: }
1313:
1314: else if ( x[k]-(user->xmin) < 0.1*h ) {
1315: vals2[k] = 0.5;
1316: }
1317:
1318:
1319: else if ( (x[k]-(user->xmin+h)) < 0.1*h ) {
1320: vals2[k] = .15625;
1321: }
1322:
1323: else {
1324: vals2[k] = 0.0;
1325: }
1326:
1327: }
1329: VecSetValuesLocal(user->phi2,3,idx,vals2,INSERT_VALUES);
1330: /*
1331: for (k=0;k < 3; k++) {
1332: vals_sum[k] = vals1[k]*vals1[k] + vals2[k]*vals2[k];
1333: }
1334: */
1335: //VecSetValuesLocal(user->Phi2D_V,3,idx,vals_sum,INSERT_VALUES);
1336:
1337: }
1338:
1339: VecAssemblyBegin(user->phi1);
1340: VecAssemblyEnd(user->phi1);
1341: VecAssemblyBegin(user->phi2);
1342: VecAssemblyEnd(user->phi2);
1343: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_phi",FILE_MODE_WRITE,&view);
1345: VecView(user->phi1,view);
1346: VecView(user->phi2,view);
1348:
1349: //VecView(user->phi1,0);
1350: //VecView(user->phi2,0);
1351:
1352: VecPointwiseMult(user->phi1,user->phi1,user->phi1);
1353: VecPointwiseMult(user->phi2,user->phi2,user->phi2);
1354: VecView(user->phi1,view);
1355: VecView(user->phi2,view);
1357: VecCopy(user->phi1,user->Phi2D_V);
1358: VecAXPY(user->Phi2D_V,1.0,user->phi2);
1359: //VecView(user->Phi2D_V,0);
1361: VecView(user->Phi2D_V,view);
1362: PetscViewerDestroy(&view);
1363: // VecNorm(user->Phi2D_V,NORM_INFINITY,&max1);
1364: //VecMin(user->Phi2D_V,&loc1,&min1);
1365: //printf("norm phi = %f, min phi = %f\n",max1,min1);
1367: return(0);
1368:
1369: }