Actual source code: ex61.c
petsc-3.3-p7 2013-05-11
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_linesearch_monitor -VG 10 -draw_fields 1,3,4 -snes_linesearch_type 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_linesearch_monitor -VG 10 -draw_fields 1,3,4 -snes_linesearch_type 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_linesearch_monitor -VG 10 -draw_fields 1,3,4 -snes_linesearch_type 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_linesearch_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_linesearch_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_linesearch_type 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: PetscInt domain;
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*);
84: PetscErrorCode Phi_read(AppCtx*);
88: int main(int argc, char **argv)
89: {
90: PetscErrorCode ierr;
91: Vec x,r; /* Solution and residual vectors */
92: SNES snes; /* Nonlinear solver context */
93: AppCtx user; /* Application context */
94: Vec xl,xu; /* Upper and lower bounds on variables */
95: Mat J;
96: PetscScalar t=0.0,normq;
97: /* PetscViewer view_out, view_q, view_psi, view_mat;*/
98: /* PetscViewer view_rand;*/
99: IS inactiveconstraints;
100: PetscInt ninactiveconstraints,N;
101: SNESConvergedReason reason;
102: /*PetscViewer view_out, view_cv,view_eta,view_vtk_cv,view_vtk_eta;*/
103: char cv_filename[80],eta_filename[80];
104: /*PetscReal bounds[] = {1000.0,-1000.,0.0,1.0,1000.0,-1000.0,0.0,1.0,1000.0,-1000.0}; */
106: PetscInitialize(&argc,&argv, (char*)0, help);
107:
108: /* Get physics and time parameters */
109: GetParams(&user);
110: /* Create a 1D DA with dof = 5; the whole thing */
111: 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);
112:
113: /* Create a 1D DA with dof = 1; for individual componentes */
114: 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);
117: /* Set Element type (triangular) */
118: DMDASetElementType(user.da1,DMDA_ELEMENT_P1);
119: DMDASetElementType(user.da2,DMDA_ELEMENT_P1);
120:
121: /* Set x and y coordinates */
122: DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,user.ymin,user.ymax,PETSC_NULL,PETSC_NULL);
123: DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,user.ymin,user.ymax,PETSC_NULL,PETSC_NULL);
126: /* Get global vector x from DM (da1) and duplicate vectors r,xl,xu */
127: DMCreateGlobalVector(user.da1,&x);
128: VecGetSize(x,&N);
129: VecDuplicate(x,&r);
130: VecDuplicate(x,&xl);
131: VecDuplicate(x,&xu);
132: VecDuplicate(x,&user.q);
133:
134: /* Get global vector user->wv from da2 and duplicate other vectors */
135: DMCreateGlobalVector(user.da2,&user.wv);
136: VecDuplicate(user.wv,&user.cv);
137: VecDuplicate(user.wv,&user.wi);
138: VecDuplicate(user.wv,&user.ci);
139: VecDuplicate(user.wv,&user.eta);
140: VecDuplicate(user.wv,&user.cvi);
141: VecDuplicate(user.wv,&user.DPsiv);
142: VecDuplicate(user.wv,&user.DPsii);
143: VecDuplicate(user.wv,&user.DPsieta);
144: VecDuplicate(user.wv,&user.Pv);
145: VecDuplicate(user.wv,&user.Pi);
146: VecDuplicate(user.wv,&user.Piv);
147: VecDuplicate(user.wv,&user.logcv);
148: VecDuplicate(user.wv,&user.logci);
149: VecDuplicate(user.wv,&user.logcvi);
150: VecDuplicate(user.wv,&user.work1);
151: VecDuplicate(user.wv,&user.work2);
152: VecDuplicate(user.wv,&user.Riv);
153: VecDuplicate(user.wv,&user.phi1);
154: VecDuplicate(user.wv,&user.phi2);
155: VecDuplicate(user.wv,&user.Phi2D_V);
156: VecDuplicate(user.wv,&user.Sv);
157: VecDuplicate(user.wv,&user.Si);
158: VecDuplicate(user.wv,&user.work3);
159: VecDuplicate(user.wv,&user.work4);
160: /* for twodomain modeling */
161: VecDuplicate(user.wv,&user.phi1);
162: VecDuplicate(user.wv,&user.phi2);
163: VecDuplicate(user.wv,&user.Phi2D_V);
164: VecDuplicate(user.wv,&user.Sv);
165: VecDuplicate(user.wv,&user.Si);
166: VecDuplicate(user.wv,&user.work3);
167: VecDuplicate(user.wv,&user.work4);
170: /* Get Jacobian matrix structure from the da for the entire thing, da1 */
171: DMCreateMatrix(user.da1,MATAIJ,&user.M);
172: /* Get the (usual) mass matrix structure from da2 */
173: DMCreateMatrix(user.da2,MATAIJ,&user.M_0);
174: SetInitialGuess(x,&user);
175: /* twodomain modeling */
176: if (user.domain) {
177: switch (user.domain) {
178: case 1:
179: Phi(&user);
180: break;
181: case 2:
182: Phi_read(&user);
183: break ;
184: }
185: }
187: /* Form the jacobian matrix and M_0 */
188: SetUpMatrices(&user);
189: MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
190:
191: SNESCreate(PETSC_COMM_WORLD,&snes);
192: SNESSetDM(snes,user.da1);
194: SNESSetFunction(snes,r,FormFunction,(void*)&user);
195: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
196:
198: SetVariableBounds(user.da1,xl,xu);
199: SNESVISetVariableBounds(snes,xl,xu);
200: SNESSetTolerances(snes,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,100,PETSC_DEFAULT);
201: SNESSetFromOptions(snes);
202:
203: /* PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_rand",FILE_MODE_WRITE,&view_rand);
204: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat2",FILE_MODE_WRITE,&view_mat);
205: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
206: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
207: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);*/
208:
209: /*PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
210:
211: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_cv",FILE_MODE_WRITE,&view_cv);
212: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_eta",FILE_MODE_WRITE,&view_eta);*/
213:
214: /* PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),5,bounds); */
215: if (user.graphics) {
216: VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
217: }
218: /*
219: if (user.graphicsfile) {
220: DMView(user.da1,user.graphicsfile);
221: VecView(x,user.graphicsfile);
222: }
223: */
224: while (t<user.T) {
225: SNESSetFunction(snes,r,FormFunction,(void*)&user);
226: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
228: SetRandomVectors(&user,t);
229: /* VecView(user.Pv,view_rand);
230: VecView(user.Pi,view_rand);
231: VecView(user.Piv,view_rand);*/
233: DPsi(&user);
234: /* VecView(user.DPsiv,view_psi);
235: VecView(user.DPsii,view_psi);
236: VecView(user.DPsieta,view_psi);*/
238: Update_q(&user);
240: /* VecView(user.q,view_q);*/
241: /* MatView(user.M,view_mat);*/
243:
244:
245: sprintf(cv_filename,"file_cv_%f.vtk",t);
246: sprintf(eta_filename,"file_eta_%f.vtk",t);
247: /* PetscViewerASCIIOpen(PETSC_COMM_WORLD,cv_filename,&view_vtk_cv);
248: PetscViewerASCIIOpen(PETSC_COMM_WORLD,eta_filename,&view_vtk_eta);
249: PetscViewerSetFormat(view_vtk_cv, PETSC_VIEWER_ASCII_VTK);
250: PetscViewerSetFormat(view_vtk_eta, PETSC_VIEWER_ASCII_VTK);
251: DMView(user.da2,view_vtk_cv);
252: DMView(user.da2,view_vtk_eta);
253: VecView(user.cv,view_cv);
254: VecView(user.eta,view_eta);
255: VecView(user.cv,view_vtk_cv);
256: VecView(user.eta,view_vtk_eta);
257: PetscViewerDestroy(&view_vtk_cv);
258: PetscViewerDestroy(&view_vtk_eta);*/
260:
261: VecNorm(user.q,NORM_2,&normq);
262: printf("2-norm of q = %14.12f\n",normq);
263: SNESSolve(snes,PETSC_NULL,x);
264: SNESGetConvergedReason(snes,&reason);
265: if (reason < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Nonlinear solver failed");
266: SNESVIGetInactiveSet(snes,&inactiveconstraints);
267: ISGetSize(inactiveconstraints,&ninactiveconstraints);
268: /* if (ninactiveconstraints < .90*N) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP,"To many active constraints, model has become non-physical"); */
270: /* VecView(x,view_out);*/
271: if (user.graphics) {
272: VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
273: }
274: /* VecView(x,PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD));*/
275: PetscInt its;
276: SNESGetIterationNumber(snes,&its);
277: PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %g in %d iterations\n",t,its);
279: Update_u(x,&user);
280: UpdateMatrices(&user);
281: t = t + user.dt;
282: /*
283: if (user.graphicsfile) {
284: VecView(x,user.graphicsfile);
285: }
286: */
287: }
288:
289: /* PetscViewerDestroy(&view_rand);
290: PetscViewerDestroy(&view_mat);
291: PetscViewerDestroy(&view_q);
292: PetscViewerDestroy(&view_out);
293: PetscViewerDestroy(&view_psi);
294: PetscViewerDestroy(&view_out);
295:
296: PetscViewerDestroy(&view_cv);
297: PetscViewerDestroy(&view_eta);*/
298:
299: if (user.graphicsfile) {
300: PetscViewerDestroy(&user.graphicsfile);
301: }
302:
303: VecDestroy(&x);
304: VecDestroy(&r);
305: VecDestroy(&xl);
306: VecDestroy(&xu);
307: VecDestroy(&user.q);
308: VecDestroy(&user.wv);
309: VecDestroy(&user.cv);
310: VecDestroy(&user.wi);
311: VecDestroy(&user.ci);
312: VecDestroy(&user.eta);
313: VecDestroy(&user.cvi);
314: VecDestroy(&user.DPsiv);
315: VecDestroy(&user.DPsii);
316: VecDestroy(&user.DPsieta);
317: VecDestroy(&user.Pv);
318: VecDestroy(&user.Pi);
319: VecDestroy(&user.Piv);
320: VecDestroy(&user.logcv);
321: VecDestroy(&user.logci);
322: VecDestroy(&user.logcvi);
323: VecDestroy(&user.work1);
324: VecDestroy(&user.work2);
325: VecDestroy(&user.Riv);
326: MatDestroy(&user.M);
327: MatDestroy(&user.M_0);
328: DMDestroy(&user.da1);
329: DMDestroy(&user.da2);
331:
332: PetscFinalize();
333: return 0;
334: }
338: PetscErrorCode Update_u(Vec X,AppCtx *user)
339: {
341: PetscInt i,n;
342: PetscScalar *xx,*wv_p,*cv_p,*wi_p,*ci_p,*eta_p;
343:
345: VecGetLocalSize(user->wv,&n);
346: VecGetArray(X,&xx);
347: VecGetArray(user->wv,&wv_p);
348: VecGetArray(user->cv,&cv_p);
349: VecGetArray(user->wi,&wi_p);
350: VecGetArray(user->ci,&ci_p);
351: VecGetArray(user->eta,&eta_p);
354: for(i=0;i<n;i++) {
355: wv_p[i] = xx[5*i];
356: cv_p[i] = xx[5*i+1];
357: wi_p[i] = xx[5*i+2];
358: ci_p[i] = xx[5*i+3];
359: eta_p[i] = xx[5*i+4];
360: }
361: VecRestoreArray(X,&xx);
362: VecRestoreArray(user->wv,&wv_p);
363: VecRestoreArray(user->cv,&cv_p);
364: VecRestoreArray(user->wi,&wi_p);
365: VecRestoreArray(user->ci,&ci_p);
366: VecRestoreArray(user->eta,&eta_p);
367:
368: return(0);
369: }
373: PetscErrorCode Update_q(AppCtx *user)
374: {
376: PetscScalar *q_p,*w1,*w2,max1;
377: PetscInt i,n;
379:
381:
382: VecPointwiseMult(user->Riv,user->eta,user->eta);
383: VecScale(user->Riv,user->Rsurf);
384: VecShift(user->Riv,user->Rbulk);
385: VecPointwiseMult(user->Riv,user->ci,user->Riv);
386: VecPointwiseMult(user->Riv,user->cv,user->Riv);
387:
388: VecCopy(user->Riv,user->work1);
389: if (user->radiation) {
390: VecAXPY(user->work1,-1.0,user->Pv);
391: }
392: if (user->domain) {
393: VecCopy(user->cv,user->work3);
394: VecShift(user->work3,-1.0*user->cv_eq);
395: VecCopy(user->Phi2D_V,user->work4);
396: VecScale(user->work4,-1.0);
397: VecShift(user->work4,1.0);
398: VecPointwiseMult(user->work4,user->work4,user->work3);
399: VecScale(user->work4,user->Svr);
400: VecAXPY(user->work1,300.0,user->work4);
401: }
402: VecScale(user->work1,user->dt);
403: VecAXPY(user->work1,-1.0,user->cv);
404: MatMult(user->M_0,user->work1,user->work2);
405:
406: VecGetArray(user->q,&q_p);
407: VecGetArray(user->work1,&w1);
408: VecGetArray(user->work2,&w2);
409: VecGetLocalSize(user->wv,&n);
410: for (i=0;i<n;i++) {
411: q_p[5*i]=w2[i];
412: }
414: MatMult(user->M_0,user->DPsiv,user->work1);
415: for (i=0;i<n;i++) {
416: q_p[5*i+1]=w1[i];
417: }
419: VecCopy(user->Riv,user->work1);
420: if (user->radiation) {
421: VecAXPY(user->work1,-1.0,user->Pi);
422: }
423: if (user->domain) {
424: VecCopy(user->ci,user->work3);
425: VecShift(user->work3,-1.0*user->ci_eq);
426: VecCopy(user->Phi2D_V,user->work4);
427: VecScale(user->work4,-1.0);
428: VecShift(user->work4,1.0);
429: VecPointwiseMult(user->work4,user->work4,user->work3);
430: VecScale(user->work4,user->Sir);
431: VecAXPY(user->work1,300.0,user->work4);
432: }
433: VecScale(user->work1,user->dt);
434: VecAXPY(user->work1,-1.0,user->ci);
435: MatMult(user->M_0,user->work1,user->work2);
436: for (i=0;i<n;i++) {
437: q_p[5*i+2]=w2[i];
438: }
440: MatMult(user->M_0,user->DPsii,user->work1);
441: for (i=0;i<n;i++) {
442: q_p[5*i+3]=w1[i];
443: }
445: VecCopy(user->DPsieta,user->work1);
446: VecScale(user->work1,user->L);
447: VecAXPY(user->work1,-1.0/user->dt,user->eta);
448: if (user->radiation) {
449: VecAXPY(user->work1,-1.0,user->Piv);
450: }
451: MatMult(user->M_0,user->work1,user->work2);
452: for (i=0;i<n;i++) {
453: q_p[5*i+4]=w2[i];
454: }
455:
456: VecRestoreArray(user->q,&q_p);
457: VecRestoreArray(user->work1,&w1);
458: VecRestoreArray(user->work2,&w2);
459:
460: return(0);
461: }
465: PetscErrorCode DPsi(AppCtx* user)
466: {
467: PetscErrorCode ierr;
468: PetscScalar Evf=user->Evf,Eif=user->Eif,kBT=user->kBT,A=user->A;
469: PetscScalar *cv_p,*ci_p,*eta_p,*logcv_p,*logci_p,*logcvi_p,*DPsiv_p,*DPsii_p,*DPsieta_p;
470: PetscInt n,i;
474: VecGetLocalSize(user->cv,&n);
475: VecGetArray(user->cv,&cv_p);
476: VecGetArray(user->ci,&ci_p);
477: VecGetArray(user->eta,&eta_p);
478: VecGetArray(user->logcv,&logcv_p);
479: VecGetArray(user->logci,&logci_p);
480: VecGetArray(user->logcvi,&logcvi_p);
481: VecGetArray(user->DPsiv,&DPsiv_p);
482: VecGetArray(user->DPsii,&DPsii_p);
483: VecGetArray(user->DPsieta,&DPsieta_p);
484:
485: Llog(user->cv,user->logcv);
486:
487: Llog(user->ci,user->logci);
488:
490: VecCopy(user->cv,user->cvi);
491: VecAXPY(user->cvi,1.0,user->ci);
492: VecScale(user->cvi,-1.0);
493: VecShift(user->cvi,1.0);
494: Llog(user->cvi,user->logcvi);
495:
496: for (i=0;i<n;i++)
497: {
498: 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);
500: 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] ;
501:
502: 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]);
503:
504:
505: }
507: VecRestoreArray(user->cv,&cv_p);
508: VecRestoreArray(user->ci,&ci_p);
509: VecRestoreArray(user->eta,&eta_p);
510: VecRestoreArray(user->logcv,&logcv_p);
511: VecRestoreArray(user->logci,&logci_p);
512: VecRestoreArray(user->logcvi,&logcvi_p);
513: VecRestoreArray(user->DPsiv,&DPsiv_p);
514: VecRestoreArray(user->DPsii,&DPsii_p);
515: VecRestoreArray(user->DPsieta,&DPsieta_p);
518: return(0);
520: }
525: PetscErrorCode Llog(Vec X, Vec Y)
526: {
527: PetscErrorCode ierr;
528: PetscScalar *x,*y;
529: PetscInt n,i;
532:
533: VecGetArray(X,&x);
534: VecGetArray(Y,&y);
535: VecGetLocalSize(X,&n);
536: for (i=0;i<n;i++) {
537: if (x[i] < 1.0e-12) {
538: y[i] = log(1.0e-12);
539: }
540: else {
541: y[i] = log(x[i]);
542: }
543: }
544:
545: return(0);
546: }
551: PetscErrorCode SetInitialGuess(Vec X,AppCtx* user)
552: {
553: PetscErrorCode ierr;
554: PetscInt n,i,Mda,Nda;
555: PetscScalar *xx,*cv_p,*ci_p,*wv_p,*wi_p,*eta;
557: /* needed for the void growth case */
558: 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;
559: PetscInt nele,nen,idx[3];
560: const PetscInt *ele;
561: PetscScalar x[3],y[3];
562: Vec coords;
563: const PetscScalar *_coords;
564: PetscViewer view;
565: PetscScalar xwidth = user->xmax - user->xmin;
569: VecGetLocalSize(X,&n);
571: if (user->voidgrowth) {
572: 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);
573: DMDAGetGhostedCoordinates(user->da2,&coords);
574: VecGetArrayRead(coords,&_coords);
576: h = (user->xmax-user->xmin)/Mda;
577: xmid = (user->xmax + user->xmin)/2.0;
578: ymid = (user->ymax + user->ymin)/2.0;
579: lambda = 4.0*h;
581: DMDAGetElements(user->da2,&nele,&nen,&ele);
582: for (i=0;i < nele; i++) {
583: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
584:
585: x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
586: x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
587: x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
588:
589: PetscInt k;
590: PetscScalar vals_cv[3],vals_ci[3],vals_eta[3],s,hhr,r;
591: for (k=0; k < 3 ; k++) {
592: s = sqrt((x[k] - xmid)*(x[k] - xmid) + (y[k] - ymid)*(y[k] - ymid));
593: if (s < xwidth*(5.0/64.0)) {
594: vals_cv[k] = cv_v;
595: vals_ci[k] = ci_v;
596: vals_eta[k] = eta_v;
597: } else if (s>= xwidth*(5.0/64.0) && s<= xwidth*(7.0/64.0) ) {
598: //r = (s - xwidth*(6.0/64.0) )/(0.5*lambda);
599: r = (s - xwidth*(6.0/64.0) )/(xwidth/64.0);
600: hhr = 0.25*(-r*r*r + 3*r + 2);
601: vals_cv[k] = cv_m + (1.0 - hhr)*(cv_v - cv_m);
602: vals_ci[k] = ci_m + (1.0 - hhr)*(ci_v - ci_m);
603: vals_eta[k] = eta_m + (1.0 - hhr)*(eta_v - eta_m);
604: } else {
605: vals_cv[k] = cv_m;
606: vals_ci[k] = ci_m;
607: vals_eta[k] = eta_m;
608: }
609: }
610: VecSetValuesLocal(user->cv,3,idx,vals_cv,INSERT_VALUES);
611: VecSetValuesLocal(user->ci,3,idx,vals_ci,INSERT_VALUES);
612: VecSetValuesLocal(user->eta,3,idx,vals_eta,INSERT_VALUES);
613: }
614: VecAssemblyBegin(user->cv);
615: VecAssemblyEnd(user->cv);
616: VecAssemblyBegin(user->ci);
617: VecAssemblyEnd(user->ci);
618: VecAssemblyBegin(user->eta);
619: VecAssemblyEnd(user->eta);
621: DMDARestoreElements(user->da2,&nele,&nen,&ele);
622: VecRestoreArrayRead(coords,&_coords);
624: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view);
625: VecView(user->cv,view);
626: VecView(user->ci,view);
627: VecView(user->eta,view);
628: PetscViewerDestroy(&view);
629: }
630: else {
631: //VecSet(user->cv,user->initv);
632: //VecSet(user->ci,user->initv);
633: VecSet(user->cv,.05);
634: VecSet(user->ci,.05);
635: VecSet(user->eta,user->initeta);
636: }
638: DPsi(user);
639: VecCopy(user->DPsiv,user->wv);
640: VecCopy(user->DPsii,user->wi);
642: VecGetArray(X,&xx);
643: VecGetArray(user->cv,&cv_p);
644: VecGetArray(user->ci,&ci_p);
645: VecGetArray(user->wv,&wv_p);
646: VecGetArray(user->wi,&wi_p);
647: VecGetArray(user->eta,&eta);
648: for (i=0;i<n/5;i++)
649: {
650: xx[5*i]=wv_p[i];
651: xx[5*i+1]=cv_p[i];
652: xx[5*i+2]=wi_p[i];
653: xx[5*i+3]=ci_p[i];
654: xx[5*i+4]=eta[i];
655: }
657: /* VecView(user->wv,view);
658: VecView(user->cv,view);
659: VecView(user->wi,view);
660: VecView(user->ci,view);
661: PetscViewerDestroy(&view);*/
663: VecRestoreArray(X,&xx);
664: VecRestoreArray(user->cv,&cv_p);
665: VecRestoreArray(user->ci,&ci_p);
666: VecRestoreArray(user->wv,&wv_p);
667: VecRestoreArray(user->wi,&wi_p);
668: VecRestoreArray(user->eta,&eta);
669:
670: return(0);
671: }
673: typedef struct {
674: PetscReal dt,x,y,strength;
675: } RandomValues;
680: PetscErrorCode SetRandomVectors(AppCtx* user,PetscReal t)
681: {
682: PetscErrorCode ierr;
683: static RandomValues *randomvalues = 0;
684: static PetscInt randindex = 0,n; /* indicates how far into the randomvalues we have currently used */
685: static PetscReal randtime = 0; /* indicates time of last radiation event */
686: PetscInt i,j,M,N,cnt = 0;
687: PetscInt xs,ys,xm,ym;
690: if (!randomvalues) {
691: PetscViewer viewer;
692: char filename[PETSC_MAX_PATH_LEN];
693: PetscBool flg;
694: PetscInt seed;
696: PetscOptionsGetInt(PETSC_NULL,"-random_seed",&seed,&flg);
697: if (flg) {
698: sprintf(filename,"ex61.random.%d",(int)seed);
699: } else {
700: PetscStrcpy(filename,"ex61.random");
701: }
702: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
703: PetscViewerBinaryRead(viewer,&n,1,PETSC_INT);
704: PetscMalloc(n*sizeof(RandomValues),&randomvalues);
705: PetscViewerBinaryRead(viewer,randomvalues,4*n,PETSC_DOUBLE);
706: for (i=0; i<n; i++) randomvalues[i].dt = randomvalues[i].dt*user->dtevent;
707: PetscViewerDestroy(&viewer);
708: }
709: user->maxevents = PetscMin(user->maxevents,n);
711: VecSet(user->Pv,0.0);
712: DMDAGetInfo(user->da1,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
713: DMDAGetGhostCorners(user->da1,&xs,&ys,0,&xm,&ym,0);
714: while (user->maxevents > randindex && randtime + randomvalues[randindex].dt < t + user->dt) { /* radiation event has occured since last time step */
715: i = ((PetscInt) (randomvalues[randindex].x*M)) - xs;
716: j = ((PetscInt) (randomvalues[randindex].y*N)) - ys;
717: if (i >= 0 && i < xm && j >= 0 && j < ym) { /* point is on this process */
718: /* need to make sure eta at the given point is not great than .8 */
719: VecSetValueLocal(user->Pv,i + 1 + xm*(j + 1), randomvalues[randindex].strength*user->VG,INSERT_VALUES);
720: }
721: randtime += randomvalues[randindex++].dt;
722: cnt++;
723: }
724: PetscPrintf(PETSC_COMM_WORLD,"Number of radiation events %d\n",cnt);
725: VecAssemblyBegin(user->Pv);
726: VecAssemblyEnd(user->Pv);
728: VecCopy(user->Pv,user->Pi);
729: VecScale(user->Pi,0.9);
730: VecPointwiseMult(user->Piv,user->Pi,user->Pv);
731: return(0);
732:
733: }
737: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ctx)
738: {
740: AppCtx *user=(AppCtx*)ctx;
741:
743: MatMultAdd(user->M,X,user->q,F);
744: return(0);
745: }
749: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
750: {
752: AppCtx *user=(AppCtx*)ctx;
753:
755: *flg = SAME_NONZERO_PATTERN;
756: MatCopy(user->M,*J,*flg);
757: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
758: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
759: return(0);
760: }
763: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
764: {
766: PetscScalar ***l,***u;
767: PetscInt xs,xm,ys,ym;
768: PetscInt j,i;
769:
771: DMDAVecGetArrayDOF(da,xl,&l);
772: DMDAVecGetArrayDOF(da,xu,&u);
773:
774: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
775:
776: for (j=ys; j<ys+ym; j++) {
777: for(i=xs; i < xs+xm;i++) {
778: l[j][i][0] = -SNES_VI_INF;
779: l[j][i][1] = 0.0;
780: l[j][i][2] = -SNES_VI_INF;
781: l[j][i][3] = 0.0;
782: l[j][i][4] = 0.0;
783: u[j][i][0] = SNES_VI_INF;
784: u[j][i][1] = 1.0;
785: u[j][i][2] = SNES_VI_INF;
786: u[j][i][3] = 1.0;
787: u[j][i][4] = 1.0;
788: }
789: }
790:
791:
792: DMDAVecRestoreArrayDOF(da,xl,&l);
793: DMDAVecRestoreArrayDOF(da,xu,&u);
794: return(0);
795: }
799: PetscErrorCode GetParams(AppCtx* user)
800: {
802: PetscBool flg,graphicsfile = PETSC_FALSE;
803:
805:
806: /* Set default parameters */
807: user->xmin = 0.0; user->xmax = 128.0;
808: user->ymin = 0.0; user->ymax = 128.0;
809: user->Dv = 1.0;
810: user->Di = 4.0;
811: user->Evf = 0.8;
812: user->Eif = 1.2;
813: user->A = 1.0;
814: user->kBT = 0.11;
815: user->kav = 1.0;
816: user->kai = 1.0;
817: user->kaeta = 1.0;
818: user->Rsurf = 10.0;
819: user->Rbulk = 1.0;
820: user->VG = 100.0;
821: user->L = 10.0;
823: user->T = 1.0e-2;
824: user->dt = 1.0e-4;
825: user->initv = .00069;
826: user->initeta = 0.0;
827: user->degenerate = PETSC_FALSE;
828: user->maxevents = 10;
829: user->graphics = PETSC_TRUE;
831: /* multidomain modeling */
832: user->domain = 1;
833: user->Svr = 0.5;
834: user->Sir = 0.5;
835: user->cv_eq = 6.9e-4;
836: user->ci_eq = 6.9e-4;
837: /* void growth */
838: user->voidgrowth = PETSC_FALSE;
840: user->radiation = PETSC_FALSE;
842: /* degenerate mobility */
843: user->smallnumber = 1.0e-3;
844: PetscOptionsBegin(PETSC_COMM_WORLD,PETSC_NULL,"Coupled Cahn-Hillard/Allen-Cahn Equations","Phasefield");
845:
846: PetscOptionsInt("-domain","Number of domains (0=one domain, 1=two domains, 2=multidomain\n","None",user->domain,&user->domain,&flg);
848: PetscOptionsReal("-Dv","Vacancy Diffusivity\n","None",user->Dv,&user->Dv,&flg);
849: PetscOptionsReal("-Di","Interstitial Diffusivity\n","None",user->Di,&user->Di,&flg);
850: PetscOptionsReal("-Evf","Vacancy Formation Energy\n","None",user->Evf,&user->Evf,&flg);
851: PetscOptionsReal("-Eif","Interstitial Formation energy\n","None",user->Eif,&user->Eif,&flg);
852: PetscOptionsReal("-A","???","None",user->A,&user->A,&flg);
853: PetscOptionsReal("-kBT","Boltzmann's Constant times the Absolute Temperature","None",user->kBT,&user->kBT,&flg);
854: PetscOptionsReal("-kav","???","None",user->kav,&user->kav,&flg);
855: PetscOptionsReal("-kai","???","None",user->kai,&user->kai,&flg);
856: PetscOptionsReal("-kaeta","???","None",user->kaeta,&user->kaeta,&flg);
857: PetscOptionsReal("-Rsurf","???","None",user->Rsurf,&user->Rsurf,&flg);
858: PetscOptionsReal("-Rbulk","???","None",user->Rbulk,&user->Rbulk,&flg);
859: PetscOptionsReal("-VG","Maximum increase in vacancy (or interstitial) concentration due to a cascade event","None",user->VG,&user->VG,&flg);
860: PetscOptionsReal("-L","???","None",user->L,&user->L,&flg);
862: PetscOptionsReal("-initv","Initial solution of Cv and Ci","None",user->initv,&user->initv,&flg);
863: PetscOptionsReal("-initeta","Initial solution of Eta","None",user->initeta,&user->initeta,&flg);
864: PetscOptionsBool("-degenerate","Run with degenerate mobility\n","None",user->degenerate,&user->degenerate,&flg);
865: PetscOptionsReal("-smallnumber","Small number added to degenerate mobility\n","None",user->smallnumber,&user->smallnumber,&flg);
867: PetscOptionsBool("-voidgrowth","Use initial conditions for void growth\n","None",user->voidgrowth,&user->voidgrowth,&flg);
868: PetscOptionsBool("-radiation","Use initial conditions for void growth\n","None",user->radiation,&user->radiation,&flg);
869: PetscOptionsReal("-xmin","Lower X coordinate of domain\n","None",user->xmin,&user->xmin,&flg);
870: PetscOptionsReal("-xmax","Upper X coordinate of domain\n","None",user->xmax,&user->xmax,&flg);
871: PetscOptionsReal("-T","Total runtime\n","None",user->T,&user->T,&flg);
872: PetscOptionsReal("-dt","Time step\n","None",user->dt,&user->dt,&flg);
873: user->dtevent = user->dt;
874: PetscOptionsReal("-dtevent","Average time between random events\n","None",user->dtevent,&user->dtevent,&flg);
875: PetscOptionsInt("-maxevents","Maximum random events allowed\n","None",user->maxevents,&user->maxevents,&flg);
877: PetscOptionsBool("-graphics","Contour plot solutions at each timestep\n","None",user->graphics,&user->graphics,&flg);
878: PetscOptionsBool("-graphicsfile","Save solution at each timestep\n","None",graphicsfile,&graphicsfile,&flg);
879: if (graphicsfile) {
880: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"ex61.data",FILE_MODE_WRITE,&user->graphicsfile);
881: }
882: PetscOptionsEnd();
883: return(0);
884: }
889: PetscErrorCode SetUpMatrices(AppCtx* user)
890: {
891: PetscErrorCode ierr;
892: PetscInt nele,nen,i,n;
893: const PetscInt *ele;
894: PetscScalar dt=user->dt,hx,hy;
895:
896: PetscInt idx[3];
897: PetscScalar eM_0[3][3],eM_2_even[3][3],eM_2_odd[3][3];
898: PetscScalar cv_sum, ci_sum;
899: Mat M=user->M;
900: Mat M_0=user->M_0;
901: PetscInt Mda=user->Mda, Nda=user->Nda;
902: PetscScalar *cv_p,*ci_p;
903: /* newly added */
904: Vec cvlocal,cilocal;
907:
908: /* MatSetOption(M,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
909: MatSetOption(M_0,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);*/
911: /* new stuff */
912: DMGetLocalVector(user->da2,&cvlocal);
913: DMGetLocalVector(user->da2,&cilocal);
914: DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
915: DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
916: DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
917: DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
918: /* old stuff */
919: /*
920: VecGetArray(user->cv,&cv_p);
921: VecGetArray(user->ci,&ci_p);
922: */
923: /* new stuff */
924: VecGetArray(cvlocal,&cv_p);
925: VecGetArray(cilocal,&ci_p);
927: MatGetLocalSize(M,&n,PETSC_NULL);
928: 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);
929: hx = (user->xmax-user->xmin)/Mda;
930: hy = (user->ymax-user->ymin)/Nda;
932: eM_0[0][0]=eM_0[1][1]=eM_0[2][2]=hx*hy/12.0;
933: 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;
934:
935: eM_2_odd[0][0] = 1.0;
936: eM_2_odd[1][1] = eM_2_odd[2][2] = 0.5;
937: eM_2_odd[0][1] = eM_2_odd[0][2] = eM_2_odd[1][0]= eM_2_odd[2][0] = -0.5;
938: eM_2_odd[1][2] = eM_2_odd[2][1] = 0.0;
940: eM_2_even[0][0] = 1.0;
941: eM_2_even[1][1] = eM_2_even[2][2] = 0.5;
942: eM_2_even[0][1] = eM_2_even[0][2] = eM_2_even[1][0]= eM_2_even[2][0] = -0.5;
943: eM_2_even[1][2] = eM_2_even[2][1] = 0.0;
945: /* eM_2_even[1][1] = 1.0;
946: eM_2_even[0][0] = eM_2_even[2][2] = 0.5;
947: eM_2_even[0][1] = eM_2_even[1][0] = eM_2_even[1][2] = eM_2_even[2][1] = -0.5;
948: eM_2_even[0][2] = eM_2_even[2][0] = 0.0;
949: */
951: // for(k=0;k < Mda*Nda*2;k++) {
952: DMDAGetElements(user->da1,&nele,&nen,&ele);
953: for (i=0; i < nele; i++) {
954: /*
955: idx[0] = connect[k*3];
956: idx[1] = connect[k*3+1];
957: idx[2] = connect[k*3+2];
958: */
959: idx[0] = ele[3*i];
960: idx[1] = ele[3*i+1];
961: idx[2] = ele[3*i+2];
963: PetscInt row,cols[6],r,row_M_0,cols3[3];
964: PetscScalar vals[6],vals_M_0[3],vals3[3];
965:
966: for(r=0;r<3;r++) {
967: //row_M_0 = connect[k*3+r];
968: row_M_0 = idx[r];
970: vals_M_0[0]=eM_0[r][0];
971: vals_M_0[1]=eM_0[r][1];
972: vals_M_0[2]=eM_0[r][2];
973:
974:
975: MatSetValuesLocal(M_0,1,&row_M_0,3,idx,vals_M_0,ADD_VALUES);
976:
977: if (user->degenerate) {
978: cv_sum = (cv_p[idx[0]] + cv_p[idx[1]] + cv_p[idx[2]])*user->Dv/(3.0*user->kBT);
979: ci_sum = (ci_p[idx[0]] + ci_p[idx[1]] + ci_p[idx[2]])*user->Di/(3.0*user->kBT);
980: } else {
981: cv_sum = user->initv*user->Dv/user->kBT;
982: ci_sum = user->initv*user->Di/user->kBT;
983: }
984:
986:
987: row = 5*idx[r];
988: cols[0] = 5*idx[0]; vals[0] = dt*eM_2_odd[r][0]*cv_sum;
989: cols[1] = 5*idx[1]; vals[1] = dt*eM_2_odd[r][1]*cv_sum;
990: cols[2] = 5*idx[2]; vals[2] = dt*eM_2_odd[r][2]*cv_sum;
991: cols[3] = 5*idx[0]+1; vals[3] = eM_0[r][0];
992: cols[4] = 5*idx[1]+1; vals[4] = eM_0[r][1];
993: cols[5] = 5*idx[2]+1; vals[5] = eM_0[r][2];
995:
996: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
997:
998:
999: row = 5*idx[r]+1;
1000: cols[0] = 5*idx[0]; vals[0] = -1.0*eM_0[r][0];
1001: cols[1] = 5*idx[1]; vals[1] = -1.0*eM_0[r][1];
1002: cols[2] = 5*idx[2]; vals[2] = -1.0*eM_0[r][2];
1003: cols[3] = 5*idx[0]+1; vals[3] = user->kav*eM_2_odd[r][0];
1004: cols[4] = 5*idx[1]+1; vals[4] = user->kav*eM_2_odd[r][1];
1005: cols[5] = 5*idx[2]+1; vals[5] = user->kav*eM_2_odd[r][2];
1006:
1007: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
1008:
1009: row = 5*idx[r]+2;
1010: cols[0] = 5*idx[0]+2; vals[0] = dt*eM_2_odd[r][0]*ci_sum;
1011: cols[1] = 5*idx[1]+2; vals[1] = dt*eM_2_odd[r][1]*ci_sum;
1012: cols[2] = 5*idx[2]+2; vals[2] = dt*eM_2_odd[r][2]*ci_sum;
1013: cols[3] = 5*idx[0]+3; vals[3] = eM_0[r][0];
1014: cols[4] = 5*idx[1]+3; vals[4] = eM_0[r][1];
1015: cols[5] = 5*idx[2]+3; vals[5] = eM_0[r][2];
1016:
1017: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
1018:
1019:
1020: row = 5*idx[r]+3;
1021: cols[0] = 5*idx[0]+2; vals[0] = -1.0*eM_0[r][0];
1022: cols[1] = 5*idx[1]+2; vals[1] = -1.0*eM_0[r][1];
1023: cols[2] = 5*idx[2]+2; vals[2] = -1.0*eM_0[r][2];
1024: cols[3] = 5*idx[0]+3; vals[3] = user->kai*eM_2_odd[r][0];
1025: cols[4] = 5*idx[1]+3; vals[4] = user->kai*eM_2_odd[r][1];
1026: cols[5] = 5*idx[2]+3; vals[5] = user->kai*eM_2_odd[r][2];
1027:
1028: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
1029:
1030:
1031: row = 5*idx[r]+4;
1032: /*
1033: cols3[0] = 5*idx[0]+4; vals3[0] = eM_0[r][0]/dt + user->L*user->kaeta*dt*eM_2_odd[r][0];
1034: cols3[1] = 5*idx[1]+4; vals3[1] = eM_0[r][1]/dt + user->L*user->kaeta*dt*eM_2_odd[r][1];
1035: cols3[2] = 5*idx[2]+4; vals3[2] = eM_0[r][2]/dt + user->L*user->kaeta*dt*eM_2_odd[r][2];
1036: */
1037: cols3[0] = 5*idx[0]+4; vals3[0] = (eM_0[r][0]/dt + user->L*user->kaeta*eM_2_odd[r][0]);
1038: cols3[1] = 5*idx[1]+4; vals3[1] = (eM_0[r][1]/dt + user->L*user->kaeta*eM_2_odd[r][1]);
1039: cols3[2] = 5*idx[2]+4; vals3[2] = (eM_0[r][2]/dt + user->L*user->kaeta*eM_2_odd[r][2]);
1040:
1041: MatSetValuesLocal(M,1,&row,3,cols3,vals3,ADD_VALUES);
1043: }
1044: }
1046: /* new */
1047: VecRestoreArray(cvlocal,&cv_p);
1048: VecRestoreArray(cilocal,&ci_p);
1049: DMRestoreLocalVector(user->da2,&cvlocal);
1050: DMRestoreLocalVector(user->da2,&cilocal);
1051: /* old */
1052: /*
1053: VecRestoreArray(user->cv,&cv_p);
1054: VecRestoreArray(user->ci,&ci_p);
1055: */
1056: MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
1057: MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);
1058:
1059: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1060: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
1061:
1062: DMDARestoreElements(user->da1,&nele,&nen,&ele);
1063:
1065: return(0);
1066: }
1071: PetscErrorCode UpdateMatrices(AppCtx* user)
1072: {
1073: PetscErrorCode ierr;
1074: PetscInt i,n,Mda,Nda,nele,nen;
1075: const PetscInt *ele;
1076:
1077: PetscInt idx[3];
1078: PetscScalar eM_2_odd[3][3],eM_2_even[3][3],h,dt=user->dt;
1079: Mat M=user->M;
1080: PetscScalar *cv_p,*ci_p,cv_sum,ci_sum;
1081: /* newly added */
1082: Vec cvlocal,cilocal;
1085:
1086:
1087: MatGetLocalSize(M,&n,PETSC_NULL);
1088:
1089: /* new stuff */
1090: DMGetLocalVector(user->da2,&cvlocal);
1091: DMGetLocalVector(user->da2,&cilocal);
1092: DMGlobalToLocalBegin(user->da2,user->cv,INSERT_VALUES,cvlocal);
1093: DMGlobalToLocalEnd(user->da2,user->cv,INSERT_VALUES,cvlocal);
1094: DMGlobalToLocalBegin(user->da2,user->ci,INSERT_VALUES,cilocal);
1095: DMGlobalToLocalEnd(user->da2,user->ci,INSERT_VALUES,cilocal);
1096: /* new stuff */
1097: VecGetArray(cvlocal,&cv_p);
1098: VecGetArray(cilocal,&ci_p);
1099:
1100: /* old stuff */
1101: /*
1102: VecGetArray(user->cv,&cv_p);
1103: VecGetArray(user->ci,&ci_p);
1104: */
1105: 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);
1107:
1108:
1109: h = (user->xmax-user->xmin)/Mda;
1111: DMDAGetElements(user->da1,&nele,&nen,&ele);
1113: for(i=0; i < nele; i++) {
1114: /*
1115: idx[0] = connect[k*3];
1116: idx[1] = connect[k*3+1];
1117: idx[2] = connect[k*3+2];
1118: */
1119: idx[0] = ele[3*i];
1120: idx[1] = ele[3*i+1];
1121: idx[2] = ele[3*i+2];
1123: PetscInt r,row,cols[3];
1124: PetscScalar vals[3];
1125: for (r=0;r<3;r++) {
1126: row = 5*idx[r];
1127: cols[0] = 5*idx[0]; vals[0] = 0.0;
1128: cols[1] = 5*idx[1]; vals[1] = 0.0;
1129: cols[2] = 5*idx[2]; vals[2] = 0.0;
1131: /* Insert values in matrix M for 1st dof */
1132: MatSetValuesLocal(M,1,&row,3,cols,vals,INSERT_VALUES);
1133:
1134: row = 5*idx[r]+2;
1135: cols[0] = 5*idx[0]+2; vals[0] = 0.0;
1136: cols[1] = 5*idx[1]+2; vals[1] = 0.0;
1137: cols[2] = 5*idx[2]+2; vals[2] = 0.0;
1139: /* Insert values in matrix M for 3nd dof */
1140: MatSetValuesLocal(M,1,&row,3,cols,vals,INSERT_VALUES);
1141: }
1142: }
1143:
1144: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1145: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
1147: eM_2_odd[0][0] = 1.0;
1148: eM_2_odd[1][1] = eM_2_odd[2][2] = 0.5;
1149: eM_2_odd[0][1] = eM_2_odd[0][2] = eM_2_odd[1][0]= eM_2_odd[2][0] = -0.5;
1150: eM_2_odd[1][2] = eM_2_odd[2][1] = 0.0;
1152: eM_2_even[0][0] = 1.0;
1153: eM_2_even[1][1] = eM_2_even[2][2] = 0.5;
1154: eM_2_even[0][1] = eM_2_even[0][2] = eM_2_even[1][0]= eM_2_even[2][0] = -0.5;
1155: eM_2_even[1][2] = eM_2_even[2][1] = 0.0;
1157: /*
1158: eM_2_even[1][1] = 1.0;
1159: eM_2_even[0][0] = eM_2_even[2][2] = 0.5;
1160: eM_2_even[0][1] = eM_2_even[1][0] = eM_2_even[1][2] = eM_2_even[2][1] = -0.5;
1161: eM_2_even[0][2] = eM_2_even[2][0] = 0.0;
1162: */
1163:
1164: /* Get local element info */
1165: //for(k=0;k < Mda*Nda*2;k++) {
1166: for (i=0; i < nele; i++) {
1167: /*
1168: idx[0] = connect[k*3];
1169: idx[1] = connect[k*3+1];
1170: idx[2] = connect[k*3+2];
1171: */
1172: idx[0] = ele[3*i];
1173: idx[1] = ele[3*i+1];
1174: idx[2] = ele[3*i+2];
1176: PetscInt row,cols[3],r;
1177: PetscScalar vals[3];
1178:
1179: for(r=0;r<3;r++) {
1180:
1181: if (user->degenerate) {
1182: printf("smallnumber = %14.12f\n",user->smallnumber);
1183: cv_sum = (user->smallnumber + cv_p[idx[0]] + cv_p[idx[1]] + cv_p[idx[2]])*user->Dv/(3.0*user->kBT);
1184: ci_sum = (user->smallnumber + ci_p[idx[0]] + ci_p[idx[1]] + ci_p[idx[2]])*user->Di/(3.0*user->kBT);
1185: } else {
1186: cv_sum = user->initv*user->Dv/(user->kBT);
1187: ci_sum = user->initv*user->Di/user->kBT;
1188: }
1190:
1191:
1192: row = 5*idx[r];
1193: cols[0] = 5*idx[0]; vals[0] = dt*eM_2_odd[r][0]*cv_sum;
1194: cols[1] = 5*idx[1]; vals[1] = dt*eM_2_odd[r][1]*cv_sum;
1195: cols[2] = 5*idx[2]; vals[2] = dt*eM_2_odd[r][2]*cv_sum;
1196:
1197: /* Insert values in matrix M for 1st dof */
1198: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
1199:
1200:
1201: row = 5*idx[r]+2;
1202: cols[0] = 5*idx[0]+2; vals[0] = dt*eM_2_odd[r][0]*ci_sum;
1203: cols[1] = 5*idx[1]+2; vals[1] = dt*eM_2_odd[r][1]*ci_sum;
1204: cols[2] = 5*idx[2]+2; vals[2] = dt*eM_2_odd[r][2]*ci_sum;
1205:
1206: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
1208:
1209:
1210: }
1211:
1212: }
1214: /* new stuff */
1215: VecRestoreArray(cvlocal,&cv_p);
1216: VecRestoreArray(cilocal,&ci_p);
1217: DMRestoreLocalVector(user->da2,&cvlocal);
1218: DMRestoreLocalVector(user->da2,&cilocal);
1219: /* old stuff */
1220: /*
1221: VecRestoreArray(user->cv,&cv_p);
1222: VecRestoreArray(user->ci,&ci_p);
1223: */
1225: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
1226: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
1229: return(0);
1230: }
1236: PetscErrorCode Phi(AppCtx* user)
1237: {
1238: PetscErrorCode ierr;
1239: PetscScalar xmid, xqu, lambda, h,x[3],y[3];
1240: Vec coords;
1241: const PetscScalar *_coords;
1242: PetscInt nele,nen,i,idx[3],Mda,Nda;
1243: const PetscInt *ele;
1244: PetscViewer view;
1248: 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);
1249: DMDAGetGhostedCoordinates(user->da2,&coords);
1250: VecGetArrayRead(coords,&_coords);
1252: h = (user->xmax - user->xmin)/Mda;
1253: xmid = (user->xmin + user->xmax)/2.0;
1254: xqu = (user->xmin + xmid)/2.0;
1255: lambda = 4.0*h;
1258: DMDAGetElements(user->da2,&nele,&nen,&ele);
1259: for (i=0;i < nele; i++) {
1260: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
1261: //printf("idx[0]=%d,idx[1]=%d,idx[2]=%d\n",idx[0],idx[1],idx[2]);
1263: x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
1264: x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
1265: x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
1267: //printf("x[0]=%f,x[1]=%f,x[2]=%f\n",x[0],x[1],x[2]);
1268: //printf("y[0]=%f,y[1]=%f,y[2]=%f\n",y[0],y[1],y[2]);
1269:
1270: PetscScalar vals1[3],vals2[3],dist1,dist2,s1,r,hhr,xc1,xc2;
1271: PetscInt k;
1273: xc1 = user->xmin;
1274: xc2 = xmid;
1276: //VecSet(user->phi1,0.0);
1277: for (k=0;k < 3; k++) {
1278: if (x[k]-xqu > 0) {
1279: s1 = (x[k] - xqu);
1280: } else {
1281: s1 = -(x[k] - xqu);
1282: }
1283: if (x[k] - xc1 > 0) {
1284: dist1 = (x[k] - xc1);
1285: } else {
1286: dist1 = -(x[k] - xc1);
1287: }
1288: if (x[k] - xc2 > 0) {
1289: dist2 = (x[k] - xc2);
1290: } else {
1291: dist2 = -(x[k] - xc2);
1292: }
1293: if (dist1 <= 0.5*lambda) {
1294: r = (x[k]-xc1)/(0.5*lambda);
1295: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
1296: vals1[k] = hhr;
1297: }
1298: else if (dist2 <= 0.5*lambda) {
1299: r = (x[k]-xc2)/(0.5*lambda);
1300: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
1301: vals1[k] = 1.0 - hhr;
1302: }
1303: else if (s1 <= xqu - 2.0*h) {
1304: vals1[k] = 1.0;
1305: }
1306:
1307: //else if ( abs(x[k]-(user->xmax-h)) < 0.1*h ) {
1308: else if ( (user->xmax-h)-x[k] < 0.1*h ) {
1309: vals1[k] = .15625;
1310: }
1311: else {
1312: vals1[k] = 0.0;
1313: }
1314: }
1315:
1316: VecSetValuesLocal(user->phi1,3,idx,vals1,INSERT_VALUES);
1318: xc1 = xmid;
1319: xc2 = user->xmax;
1321: //VecSet(user->phi2,0.0);
1322: for (k=0;k < 3; k++) {
1323: /*
1324: s1 = abs(x[k] - (xqu+xmid));
1325: dist1 = abs(x[k] - xc1);
1326: dist2 = abs(x[k] - xc2);
1327: */
1328: if (x[k]-(xqu + xmid) > 0) {
1329: s1 = (x[k] - (xqu + xmid));
1330: } else {
1331: s1 = -(x[k] - (xqu + xmid));
1332: }
1333: if (x[k] - xc1 > 0) {
1334: dist1 = (x[k] - xc1);
1335: } else {
1336: dist1 = -(x[k] - xc1);
1337: }
1338: if (x[k] - xc2 > 0) {
1339: dist2 = (x[k] - xc2);
1340: } else {
1341: dist2 = -(x[k] - xc2);
1342: }
1343:
1344: if (dist1 <= 0.5*lambda) {
1345: r = (x[k]-xc1)/(0.5*lambda);
1346: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
1347: vals2[k] = hhr;
1348: }
1349: else if (dist2 <= 0.5*lambda) {
1350: r = -(x[k]-xc2)/(0.5*lambda);
1351: hhr = 0.25*(-r*r*r + 3.0*r + 2.0);
1352: vals2[k] = hhr;
1353: }
1354: else if (s1 <= xqu - 2.0*h) {
1355: vals2[k] = 1.0;
1356: }
1357:
1358: else if ( x[k]-(user->xmin) < 0.1*h ) {
1359: vals2[k] = 0.5;
1360: }
1361:
1362:
1363: else if ( (x[k]-(user->xmin+h)) < 0.1*h ) {
1364: vals2[k] = .15625;
1365: }
1366:
1367: else {
1368: vals2[k] = 0.0;
1369: }
1370:
1371: }
1373: VecSetValuesLocal(user->phi2,3,idx,vals2,INSERT_VALUES);
1374: /*
1375: for (k=0;k < 3; k++) {
1376: vals_sum[k] = vals1[k]*vals1[k] + vals2[k]*vals2[k];
1377: }
1378: */
1379: //VecSetValuesLocal(user->Phi2D_V,3,idx,vals_sum,INSERT_VALUES);
1380:
1381: }
1382:
1383: VecAssemblyBegin(user->phi1);
1384: VecAssemblyEnd(user->phi1);
1385: VecAssemblyBegin(user->phi2);
1386: VecAssemblyEnd(user->phi2);
1387: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_phi",FILE_MODE_WRITE,&view);
1389: VecView(user->phi1,view);
1390: VecView(user->phi2,view);
1392:
1393: //VecView(user->phi1,0);
1394: //VecView(user->phi2,0);
1395:
1396: VecPointwiseMult(user->phi1,user->phi1,user->phi1);
1397: VecPointwiseMult(user->phi2,user->phi2,user->phi2);
1398: VecView(user->phi1,view);
1399: VecView(user->phi2,view);
1401: VecCopy(user->phi1,user->Phi2D_V);
1402: VecAXPY(user->Phi2D_V,1.0,user->phi2);
1403: //VecView(user->Phi2D_V,0);
1405: VecView(user->Phi2D_V,view);
1406: PetscViewerDestroy(&view);
1407: // VecNorm(user->Phi2D_V,NORM_INFINITY,&max1);
1408: //VecMin(user->Phi2D_V,&loc1,&min1);
1409: //printf("norm phi = %f, min phi = %f\n",max1,min1);
1411: return(0);
1412:
1413: }
1417: PetscErrorCode Phi_read(AppCtx* user)
1418: {
1419: PetscErrorCode ierr;
1420: PetscReal *values;
1421: PetscViewer viewer;
1424:
1425: VecGetArray(user->Phi2D_V,&values);
1426: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi3",FILE_MODE_READ,&viewer);
1427: PetscViewerBinaryRead(viewer,values,16384,PETSC_DOUBLE);
1428: PetscViewerDestroy(&viewer);
1429: return(0);
1430: }