Actual source code: ex60.c
petsc-3.5.4 2015-05-23
1: static char help[] = "2D coupled Allen-Cahn and Cahn-Hilliard equation for constant 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: ./ex60 -ksp_type fgmres -pc_type mg -snes_vi_monitor -snes_atol 1.e-11 -da_grid_x 72 -da_grid_y 72 -ksp_rtol 1.e-8 -T 0.1 -VG 100 -pc_type lu -ksp_monitor_true_residual -pc_factor_mat_solver_package superlu -snes_converged_reason -ksp_converged_reason -pc_type sor -ksp_rtol 1.e-9 -snes_linesearch_monitor -VG 10 -draw_fields 1,3,4 -snes_monitor_solution
14: */
16: /*
17: 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
19: Add command-line option for constant or degenerate mobility
20: Add command-line option for graphics at each time step
22: Check time-step business; what should it be? How to check that it is good?
23: Make random right hand side forcing function proportional to time step so smaller time steps don't mean more radiation
24: How does the multigrid linear solver work now?
25: What happens when running with degenerate mobility
28: */
30: #include <petscsnes.h>
31: #include <petscdmda.h>
33: typedef struct {
34: PetscReal dt,T; /* Time step and end time */
35: DM da1,da2;
36: Mat M; /* Jacobian matrix */
37: Mat M_0;
38: Vec q,wv,cv,wi,ci,eta,cvi,DPsiv,DPsii,DPsieta,Pv,Pi,Piv,logcv,logci,logcvi,Rr,Riv;
39: Vec work1,work2,work3,work4;
40: PetscScalar Dv,Di,Evf,Eif,A,kBT,kav,kai,kaeta,Rsurf,Rbulk,L,P_casc,VG; /* physics parameters */
41: PetscReal xmin,xmax,ymin,ymax;
42: PetscInt Mda, Nda;
43: } AppCtx;
45: PetscErrorCode GetParams(AppCtx*);
46: PetscErrorCode SetRandomVectors(AppCtx*);
47: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
48: PetscErrorCode SetUpMatrices(AppCtx*);
49: PetscErrorCode UpdateMatrices(AppCtx*);
50: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
51: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
52: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
53: PetscErrorCode Update_q(AppCtx*);
54: PetscErrorCode Update_u(Vec,AppCtx*);
55: PetscErrorCode DPsi(AppCtx*);
56: PetscErrorCode LaplacianFiniteDifference(AppCtx*);
57: PetscErrorCode Llog(Vec,Vec);
60: int main(int argc, char **argv)
61: {
63: Vec x,r; /* Solution and residual vectors */
64: SNES snes; /* Nonlinear solver context */
65: AppCtx user; /* Application context */
66: Vec xl,xu; /* Upper and lower bounds on variables */
67: Mat J;
68: PetscScalar t=0.0;
69: PetscViewer view_out, view_q, view_psi, view_mat;
70: PetscViewer view_rand;
71: IS inactiveconstraints;
72: PetscInt ninactiveconstraints,N;
74: PetscInitialize(&argc,&argv, (char*)0, help);
76: /* Get physics and time parameters */
77: GetParams(&user);
78: /* Create a 1D DA with dof = 5; the whole thing */
79: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX, -3,-3,PETSC_DECIDE,PETSC_DECIDE, 5, 1,NULL,NULL,&user.da1);
81: /* Create a 1D DA with dof = 1; for individual componentes */
82: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX, -3,-3,PETSC_DECIDE,PETSC_DECIDE, 1, 1,NULL,NULL,&user.da2);
85: /* Set Element type (triangular) */
86: DMDASetElementType(user.da1,DMDA_ELEMENT_P1);
87: DMDASetElementType(user.da2,DMDA_ELEMENT_P1);
89: /* Set x and y coordinates */
90: DMDASetUniformCoordinates(user.da1,user.xmin,user.xmax,user.ymin,user.ymax,NULL,NULL);
91: DMDASetUniformCoordinates(user.da2,user.xmin,user.xmax,user.ymin,user.ymax,NULL,NULL);
92: /* Get global vector x from DM (da1) and duplicate vectors r,xl,xu */
93: DMCreateGlobalVector(user.da1,&x);
94: VecGetSize(x,&N);
95: VecDuplicate(x,&r);
96: VecDuplicate(x,&xl);
97: VecDuplicate(x,&xu);
98: VecDuplicate(x,&user.q);
100: /* Get global vector user->wv from da2 and duplicate other vectors */
101: DMCreateGlobalVector(user.da2,&user.wv);
102: VecDuplicate(user.wv,&user.cv);
103: VecDuplicate(user.wv,&user.wi);
104: VecDuplicate(user.wv,&user.ci);
105: VecDuplicate(user.wv,&user.eta);
106: VecDuplicate(user.wv,&user.cvi);
107: VecDuplicate(user.wv,&user.DPsiv);
108: VecDuplicate(user.wv,&user.DPsii);
109: VecDuplicate(user.wv,&user.DPsieta);
110: VecDuplicate(user.wv,&user.Pv);
111: VecDuplicate(user.wv,&user.Pi);
112: VecDuplicate(user.wv,&user.Piv);
113: VecDuplicate(user.wv,&user.logcv);
114: VecDuplicate(user.wv,&user.logci);
115: VecDuplicate(user.wv,&user.logcvi);
116: VecDuplicate(user.wv,&user.work1);
117: VecDuplicate(user.wv,&user.work2);
118: VecDuplicate(user.wv,&user.Rr);
119: VecDuplicate(user.wv,&user.Riv);
122: /* Get Jacobian matrix structure from the da for the entire thing, da1 */
123: DMSetMatType(user.da1,MATAIJ);
124: DMCreateMatrix(user.da1,&user.M);
125: /* Get the (usual) mass matrix structure from da2 */
126: DMSetMatType(user.da2,MATAIJ);
127: DMCreateMatrix(user.da2,&user.M_0);
128: SetInitialGuess(x,&user);
129: /* Form the jacobian matrix and M_0 */
130: SetUpMatrices(&user);
131: MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
133: SNESCreate(PETSC_COMM_WORLD,&snes);
134: SNESSetDM(snes,user.da1);
136: SNESSetFunction(snes,r,FormFunction,(void*)&user);
137: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
139: SetVariableBounds(user.da1,xl,xu);
140: SNESVISetVariableBounds(snes,xl,xu);
141: SNESSetFromOptions(snes);
143: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_rand",FILE_MODE_WRITE,&view_rand);
144: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_mat",FILE_MODE_WRITE,&view_mat);
145: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_q",FILE_MODE_WRITE,&view_q);
146: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_out",FILE_MODE_WRITE,&view_out);
147: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_psi",FILE_MODE_WRITE,&view_psi);
149: while (t<user.T) {
150: SNESSetFunction(snes,r,FormFunction,(void*)&user);
151: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
153: SetRandomVectors(&user);
154: /* VecView(user.Pv,view_rand);
155: VecView(user.Pi,view_rand);
156: VecView(user.Piv,view_rand);*/
158: DPsi(&user);
159: /* VecView(user.DPsiv,view_psi);
160: VecView(user.DPsii,view_psi);
161: VecView(user.DPsieta,view_psi);*/
163: Update_q(&user);
164: /* VecView(user.q,view_q);
165: MatView(user.M,view_mat);*/
166: SNESSolve(snes,NULL,x);
167: SNESVIGetInactiveSet(snes,&inactiveconstraints);
168: ISGetSize(inactiveconstraints,&ninactiveconstraints);
169: /* if (ninactiveconstraints < .90*N) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP,"To many active constraints, model has become non-physical"); */
171: /* VecView(x,view_out);*/
172: VecView(x,PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD));
173: PetscInt its;
174: SNESGetIterationNumber(snes,&its);
175: PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4f in %d iterations\n",t,its);
177: Update_u(x,&user);
178: UpdateMatrices(&user);
179: t = t + user.dt;
180: }
182: PetscViewerDestroy(&view_rand);
183: PetscViewerDestroy(&view_mat);
184: PetscViewerDestroy(&view_q);
185: PetscViewerDestroy(&view_out);
186: PetscViewerDestroy(&view_psi);
188: VecDestroy(&x);
189: VecDestroy(&r);
190: VecDestroy(&xl);
191: VecDestroy(&xu);
192: VecDestroy(&user.q);
193: VecDestroy(&user.wv);
194: VecDestroy(&user.cv);
195: VecDestroy(&user.wi);
196: VecDestroy(&user.ci);
197: VecDestroy(&user.eta);
198: VecDestroy(&user.cvi);
199: VecDestroy(&user.DPsiv);
200: VecDestroy(&user.DPsii);
201: VecDestroy(&user.DPsieta);
202: VecDestroy(&user.Pv);
203: VecDestroy(&user.Pi);
204: VecDestroy(&user.Piv);
205: VecDestroy(&user.logcv);
206: VecDestroy(&user.logci);
207: VecDestroy(&user.logcvi);
208: VecDestroy(&user.work1);
209: VecDestroy(&user.work2);
210: VecDestroy(&user.Rr);
211: VecDestroy(&user.Riv);
212: MatDestroy(&user.M);
213: MatDestroy(&user.M_0);
214: DMDestroy(&user.da1);
215: DMDestroy(&user.da2);
217: PetscFinalize();
218: return 0;
219: }
223: PetscErrorCode Update_u(Vec X,AppCtx *user)
224: {
226: PetscInt i,n;
227: PetscScalar *xx,*wv_p,*cv_p,*wi_p,*ci_p,*eta_p;
230: VecGetLocalSize(user->wv,&n);
231: VecGetArray(X,&xx);
232: VecGetArray(user->wv,&wv_p);
233: VecGetArray(user->cv,&cv_p);
234: VecGetArray(user->wi,&wi_p);
235: VecGetArray(user->ci,&ci_p);
236: VecGetArray(user->eta,&eta_p);
239: for (i=0; i<n; i++) {
240: wv_p[i] = xx[5*i];
241: cv_p[i] = xx[5*i+1];
242: wi_p[i] = xx[5*i+2];
243: ci_p[i] = xx[5*i+3];
244: eta_p[i] = xx[5*i+4];
245: }
246: VecRestoreArray(X,&xx);
247: VecRestoreArray(user->wv,&wv_p);
248: VecRestoreArray(user->cv,&cv_p);
249: VecRestoreArray(user->wi,&wi_p);
250: VecRestoreArray(user->ci,&ci_p);
251: VecRestoreArray(user->eta,&eta_p);
252: return(0);
253: }
257: PetscErrorCode Update_q(AppCtx *user)
258: {
260: PetscScalar *q_p,*w1,*w2;
261: PetscInt i,n;
264: VecPointwiseMult(user->Rr,user->eta,user->eta);
265: VecScale(user->Rr,user->Rsurf);
266: VecShift(user->Rr,user->Rbulk);
267: VecPointwiseMult(user->Riv,user->cv,user->ci);
268: VecPointwiseMult(user->Riv,user->Rr,user->Riv);
270: VecGetArray(user->q,&q_p);
271: VecGetArray(user->work1,&w1);
272: VecGetArray(user->work2,&w2);
274: VecCopy(user->cv,user->work1);
275: VecAXPY(user->work1,1.0,user->Pv);
276: VecScale(user->work1,-1.0);
277: MatMult(user->M_0,user->work1,user->work2);
278: VecGetLocalSize(user->work1,&n);
280: for (i=0; i<n; i++) q_p[5*i]=w2[i];
282: MatMult(user->M_0,user->DPsiv,user->work1);
283: for (i=0; i<n; i++) q_p[5*i+1]=w1[i];
285: VecCopy(user->ci,user->work1);
286: VecAXPY(user->work1,1.0,user->Pi);
287: VecScale(user->work1,-1.0);
288: MatMult(user->M_0,user->work1,user->work2);
289: for (i=0; i<n; i++) q_p[5*i+2]=w2[i];
291: MatMult(user->M_0,user->DPsii,user->work1);
292: for (i=0; i<n; i++) q_p[5*i+3]=w1[i];
294: VecCopy(user->eta,user->work1);
295: VecScale(user->work1,-1.0/user->dt);
296: VecAXPY(user->work1,user->L,user->DPsieta);
297: VecAXPY(user->work1,-1.0,user->Piv);
298: MatMult(user->M_0,user->work1,user->work2);
299: for (i=0; i<n; i++) q_p[5*i+4]=w2[i];
301: VecRestoreArray(user->q,&q_p);
302: VecRestoreArray(user->work1,&w1);
303: VecRestoreArray(user->work2,&w2);
304: return(0);
305: }
309: PetscErrorCode DPsi(AppCtx *user)
310: {
312: PetscScalar Evf=user->Evf,Eif=user->Eif,kBT=user->kBT,A=user->A;
313: PetscScalar *cv_p,*ci_p,*eta_p,*logcv_p,*logci_p,*logcvi_p,*DPsiv_p,*DPsii_p,*DPsieta_p;
314: PetscInt n,i;
317: VecGetLocalSize(user->cv,&n);
318: VecGetArray(user->cv,&cv_p);
319: VecGetArray(user->ci,&ci_p);
320: VecGetArray(user->eta,&eta_p);
321: VecGetArray(user->logcv,&logcv_p);
322: VecGetArray(user->logci,&logci_p);
323: VecGetArray(user->logcvi,&logcvi_p);
324: VecGetArray(user->DPsiv,&DPsiv_p);
325: VecGetArray(user->DPsii,&DPsii_p);
326: VecGetArray(user->DPsieta,&DPsieta_p);
328: Llog(user->cv,user->logcv);
330: Llog(user->ci,user->logci);
333: VecCopy(user->cv,user->cvi);
334: VecAXPY(user->cvi,1.0,user->ci);
335: VecScale(user->cvi,-1.0);
336: VecShift(user->cvi,1.0);
337: Llog(user->cvi,user->logcvi);
339: for (i=0; i<n; i++)
340: {
341: 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);
342: 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];
344: 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]);
345: }
347: VecRestoreArray(user->cv,&cv_p);
348: VecRestoreArray(user->ci,&ci_p);
349: VecRestoreArray(user->eta,&eta_p);
350: VecRestoreArray(user->logcv,&logcv_p);
351: VecRestoreArray(user->logci,&logci_p);
352: VecRestoreArray(user->logcvi,&logcvi_p);
353: VecRestoreArray(user->DPsiv,&DPsiv_p);
354: VecRestoreArray(user->DPsii,&DPsii_p);
355: VecRestoreArray(user->DPsieta,&DPsieta_p);
356: return(0);
357: }
362: PetscErrorCode Llog(Vec X, Vec Y)
363: {
365: PetscScalar *x,*y;
366: PetscInt n,i;
369: VecGetArray(X,&x);
370: VecGetArray(Y,&y);
371: VecGetLocalSize(X,&n);
372: for (i=0; i<n; i++) {
373: if (x[i] < 1.0e-12) y[i] = log(1.0e-12);
374: else y[i] = log(x[i]);
375: }
376: return(0);
377: }
382: PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
383: {
385: PetscInt n,i;
386: PetscScalar *xx,*cv_p,*ci_p,*wv_p,*wi_p;
387: PetscViewer view;
388: PetscScalar initv = .00069;
391: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"file_initial",FILE_MODE_WRITE,&view);
392: VecGetLocalSize(X,&n);
394: VecSet(user->cv,initv);
395: VecSet(user->ci,initv);
396: VecSet(user->eta,0.0);
398: DPsi(user);
399: VecCopy(user->DPsiv,user->wv);
400: VecCopy(user->DPsii,user->wi);
402: VecGetArray(X,&xx);
403: VecGetArray(user->cv,&cv_p);
404: VecGetArray(user->ci,&ci_p);
405: VecGetArray(user->wv,&wv_p);
406: VecGetArray(user->wi,&wi_p);
407: for (i=0; i<n/5; i++) {
408: xx[5*i] =wv_p[i];
409: xx[5*i+1]=cv_p[i];
410: xx[5*i+2]=wi_p[i];
411: xx[5*i+3]=ci_p[i];
412: xx[5*i+4]=0.0;
413: }
415: VecView(user->wv,view);
416: VecView(user->cv,view);
417: VecView(user->wi,view);
418: VecView(user->ci,view);
419: PetscViewerDestroy(&view);
421: VecRestoreArray(X,&xx);
422: VecRestoreArray(user->cv,&cv_p);
423: VecRestoreArray(user->ci,&ci_p);
424: VecRestoreArray(user->wv,&wv_p);
425: VecRestoreArray(user->wi,&wi_p);
426: return(0);
427: }
431: PetscErrorCode SetRandomVectors(AppCtx *user)
432: {
434: PetscInt i,n,count=0;
435: PetscScalar *w1,*w2,*Pv_p,*eta_p;
438: VecSetRandom(user->work1,NULL);
439: VecSetRandom(user->work2,NULL);
440: VecGetArray(user->work1,&w1);
441: VecGetArray(user->work2,&w2);
442: VecGetArray(user->Pv,&Pv_p);
443: VecGetArray(user->eta,&eta_p);
444: VecGetLocalSize(user->work1,&n);
445: for (i=0; i<n; i++) {
447: if (eta_p[i]>=0.8 || w1[i]>user->P_casc) Pv_p[i]=0;
448: else {
449: Pv_p[i]=w2[i]*user->VG;
450: count =count+1;
451: }
453: }
455: VecCopy(user->Pv,user->Pi);
456: VecScale(user->Pi,0.9);
457: VecPointwiseMult(user->Piv,user->Pi,user->Pv);
458: VecRestoreArray(user->work1,&w1);
459: VecRestoreArray(user->work2,&w2);
460: VecRestoreArray(user->Pv,&Pv_p);
461: VecRestoreArray(user->eta,&eta_p);
462: return(0);
464: }
468: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
469: {
471: AppCtx *user=(AppCtx*)ctx;
474: MatMultAdd(user->M,X,user->q,F);
475: return(0);
476: }
480: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat B,void *ctx)
481: {
483: AppCtx *user=(AppCtx*)ctx;
486: MatCopy(user->M,J,*flg);
487: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
488: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
489: return(0);
490: }
493: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
494: {
496: PetscScalar ***l,***u;
497: PetscInt xs,xm,ys,ym;
498: PetscInt j,i;
501: DMDAVecGetArrayDOF(da,xl,&l);
502: DMDAVecGetArrayDOF(da,xu,&u);
504: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
506: for (j=ys; j<ys+ym; j++) {
507: for (i=xs; i < xs+xm; i++) {
508: l[j][i][0] = -PETSC_INFINITY;
509: l[j][i][1] = 0.0;
510: l[j][i][2] = -PETSC_INFINITY;
511: l[j][i][3] = 0.0;
512: l[j][i][4] = 0.0;
513: u[j][i][0] = PETSC_INFINITY;
514: u[j][i][1] = 1.0;
515: u[j][i][2] = PETSC_INFINITY;
516: u[j][i][3] = 1.0;
517: u[j][i][4] = 1.0;
518: }
519: }
522: DMDAVecRestoreArrayDOF(da,xl,&l);
523: DMDAVecRestoreArrayDOF(da,xu,&u);
524: return(0);
525: }
529: PetscErrorCode GetParams(AppCtx *user)
530: {
532: PetscBool flg;
535: /* Set default parameters */
536: user->xmin = 0.0; user->xmax = 1.0;
537: user->ymin = 0.0; user->ymax = 1.0;
538: user->Dv = 1.0; user->Di=4.0;
539: user->Evf = 0.8; user->Eif = 1.2;
540: user->A = 1.0;
541: user->kBT = 0.11;
542: user->kav = 1.0; user->kai = 1.0; user->kaeta = 1.0;
543: user->Rsurf = 10.0; user->Rbulk = 1.0;
544: user->L = 10.0; user->P_casc = 0.05;
545: user->T = 1.0e-2; user->dt = 1.0e-4;
546: user->VG = 100.0;
548: PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);
549: PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);
550: PetscOptionsGetReal(NULL,"-T",&user->T,&flg);
551: PetscOptionsGetReal(NULL,"-dt",&user->dt,&flg);
552: PetscOptionsGetReal(NULL,"-P_casc",&user->P_casc,&flg);
553: PetscOptionsGetReal(NULL,"-VG",&user->VG,&flg);
554: return(0);
555: }
560: PetscErrorCode SetUpMatrices(AppCtx *user)
561: {
563: PetscInt nele,nen,i,j,n;
564: const PetscInt *ele;
565: PetscScalar dt=user->dt,hx,hy;
567: PetscInt idx[3],*nodes, *connect, k;
568: PetscScalar eM_0[3][3],eM_2_even[3][3],eM_2_odd[3][3];
569: PetscScalar cv_sum, ci_sum;
570: Mat M =user->M;
571: Mat M_0=user->M_0;
572: PetscInt Mda=user->Mda, Nda=user->Nda, ld, rd, ru, lu;
573: PetscScalar *cv_p,*ci_p;
576: /* MatSetOption(M,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
577: MatSetOption(M_0,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);*/
579: /* Create the mass matrix M_0 */
580: VecGetArray(user->cv,&cv_p);
581: VecGetArray(user->ci,&ci_p);
582: MatGetLocalSize(M,&n,NULL);
583: DMDAGetInfo(user->da1,NULL,&Mda,&Nda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
585: PetscMalloc1((Mda+1)*(Nda+1),&nodes);
586: PetscMalloc1(Mda*Nda*2*3,&connect);
587: hx = (user->xmax-user->xmin)/Mda;
588: hy = (user->ymax-user->ymin)/Nda;
589: for (j=0; j < Nda; j++) {
590: for (i=0; i < Mda; i++) nodes[j*(Mda+1)+i] = j*Mda+i;
591: nodes[j*(Mda+1)+Mda] = j*Mda;
592: }
594: for (i=0; i < Mda; i++) nodes[Nda*(Mda+1)+i] = i;
596: nodes[Nda*(Mda+1)+Mda] = 0;
598: k = 0;
599: for (j=0; j<Nda; j++) {
600: for (i=0; i<Mda; i++) {
602: /* ld = nodes[j][i]; */
603: ld = nodes[j*(Mda+1)+i];
604: /* rd = nodes[j+1][i]; */
605: rd = nodes[(j+1)*(Mda+1)+i];
606: /* ru = nodes[j+1][i+1]; */
607: ru = nodes[(j+1)*(Mda+1)+i+1];
608: /* lu = nodes[j][i+1]; */
609: lu = nodes[j*(Mda+1)+i+1];
611: /* connect[k][0]=ld; */
612: connect[k*6]=ld;
613: /* connect[k][1]=lu; */
614: connect[k*6+1]=lu;
615: /* connect[k][2]=ru; */
616: connect[k*6+2]=rd;
617: connect[k*6+3]=lu;
618: connect[k*6+4]=ru;
619: connect[k*6+5]=rd;
621: k = k+1;
622: }
623: }
626: eM_0[0][0]=eM_0[1][1]=eM_0[2][2]=hx*hy/12.0;
627: 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;
629: eM_2_odd[0][0] = 1.0;
630: eM_2_odd[1][1] = eM_2_odd[2][2] = 0.5;
631: eM_2_odd[0][1] = eM_2_odd[0][2] = eM_2_odd[1][0]= eM_2_odd[2][0] = -0.5;
632: eM_2_odd[1][2] = eM_2_odd[2][1] = 0.0;
634: eM_2_even[1][1] = 1.0;
635: eM_2_even[0][0] = eM_2_even[2][2] = 0.5;
636: eM_2_even[0][1] = eM_2_even[1][0] = eM_2_even[1][2] = eM_2_even[2][1] = -0.5;
637: eM_2_even[0][2] = eM_2_even[2][0] = 0.0;
640: for (k=0; k < Mda*Nda*2; k++) {
641: idx[0] = connect[k*3];
642: idx[1] = connect[k*3+1];
643: idx[2] = connect[k*3+2];
645: PetscInt row,cols[6],r,row_M_0,cols3[3];
646: PetscScalar vals[6],vals_M_0[3],vals3[3];
648: for (r=0; r<3; r++) {
649: row_M_0 = connect[k*3+r];
651: vals_M_0[0]=eM_0[r][0];
652: vals_M_0[1]=eM_0[r][1];
653: vals_M_0[2]=eM_0[r][2];
656: MatSetValues(M_0,1,&row_M_0,3,idx,vals_M_0,ADD_VALUES);
658: /* cv_sum = (cv_p[idx[0]] + cv_p[idx[1]] + cv_p[idx[2]])*user->Dv/(3.0*user->kBT); */
659: /* ci_sum = (ci_p[idx[0]] + ci_p[idx[1]] + ci_p[idx[2]])*user->Di/(3.0*user->kBT); */
660: cv_sum = .0000069*user->Dv/user->kBT;
661: ci_sum = .0000069*user->Di/user->kBT;
663: if (k%2 == 0) {
664: row = 5*idx[r];
665: cols[0] = 5*idx[0]; vals[0] = dt*eM_2_odd[r][0]*cv_sum;
666: cols[1] = 5*idx[1]; vals[1] = dt*eM_2_odd[r][1]*cv_sum;
667: cols[2] = 5*idx[2]; vals[2] = dt*eM_2_odd[r][2]*cv_sum;
668: cols[3] = 5*idx[0]+1; vals[3] = eM_0[r][0];
669: cols[4] = 5*idx[1]+1; vals[4] = eM_0[r][1];
670: cols[5] = 5*idx[2]+1; vals[5] = eM_0[r][2];
672: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
674: row = 5*idx[r]+1;
675: cols[0] = 5*idx[0]; vals[0] = -1.0*eM_0[r][0];
676: cols[1] = 5*idx[1]; vals[1] = -1.0*eM_0[r][1];
677: cols[2] = 5*idx[2]; vals[2] = -1.0*eM_0[r][2];
678: cols[3] = 5*idx[0]+1; vals[3] = user->kav*eM_2_odd[r][0];
679: cols[4] = 5*idx[1]+1; vals[4] = user->kav*eM_2_odd[r][1];
680: cols[5] = 5*idx[2]+1; vals[5] = user->kav*eM_2_odd[r][2];
682: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
684: row = 5*idx[r]+2;
685: cols[0] = 5*idx[0]+2; vals[0] = dt*eM_2_odd[r][0]*ci_sum;
686: cols[1] = 5*idx[1]+2; vals[1] = dt*eM_2_odd[r][1]*ci_sum;
687: cols[2] = 5*idx[2]+2; vals[2] = dt*eM_2_odd[r][2]*ci_sum;
688: cols[3] = 5*idx[0]+3; vals[3] = eM_0[r][0];
689: cols[4] = 5*idx[1]+3; vals[4] = eM_0[r][1];
690: cols[5] = 5*idx[2]+3; vals[5] = eM_0[r][2];
692: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
694: row = 5*idx[r]+3;
695: cols[0] = 5*idx[0]+2; vals[0] = -1.0*eM_0[r][0];
696: cols[1] = 5*idx[1]+2; vals[1] = -1.0*eM_0[r][1];
697: cols[2] = 5*idx[2]+2; vals[2] = -1.0*eM_0[r][2];
698: cols[3] = 5*idx[0]+3; vals[3] = user->kai*eM_2_odd[r][0];
699: cols[4] = 5*idx[1]+3; vals[4] = user->kai*eM_2_odd[r][1];
700: cols[5] = 5*idx[2]+3; vals[5] = user->kai*eM_2_odd[r][2];
702: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
704: row = 5*idx[r]+4;
705: /*
706: cols3[0] = 5*idx[0]+4; vals3[0] = eM_0[r][0]/dt + user->L*user->kaeta*dt*eM_2_odd[r][0];
707: cols3[1] = 5*idx[1]+4; vals3[1] = eM_0[r][1]/dt + user->L*user->kaeta*dt*eM_2_odd[r][1];
708: cols3[2] = 5*idx[2]+4; vals3[2] = eM_0[r][2]/dt + user->L*user->kaeta*dt*eM_2_odd[r][2];
709: */
710: cols3[0] = 5*idx[0]+4; vals3[0] = eM_0[r][0]/dt + user->L*user->kaeta*eM_2_odd[r][0];
711: cols3[1] = 5*idx[1]+4; vals3[1] = eM_0[r][1]/dt + user->L*user->kaeta*eM_2_odd[r][1];
712: cols3[2] = 5*idx[2]+4; vals3[2] = eM_0[r][2]/dt + user->L*user->kaeta*eM_2_odd[r][2];
714: MatSetValuesLocal(M,1,&row,3,cols3,vals3,ADD_VALUES);
717: } else {
720: row = 5*idx[r];
721: cols[0] = 5*idx[0]; vals[0] = dt*eM_2_even[r][0]*cv_sum;
722: cols[1] = 5*idx[1]; vals[1] = dt*eM_2_even[r][1]*cv_sum;
723: cols[2] = 5*idx[2]; vals[2] = dt*eM_2_even[r][2]*cv_sum;
724: cols[3] = 5*idx[0]+1; vals[3] = eM_0[r][0];
725: cols[4] = 5*idx[1]+1; vals[4] = eM_0[r][1];
726: cols[5] = 5*idx[2]+1; vals[5] = eM_0[r][2];
728: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
730: row = 5*idx[r]+1;
731: cols[0] = 5*idx[0]; vals[0] = -1.0*eM_0[r][0];
732: cols[1] = 5*idx[1]; vals[1] = -1.0*eM_0[r][1];
733: cols[2] = 5*idx[2]; vals[2] = -1.0*eM_0[r][2];
734: cols[3] = 5*idx[0]+1; vals[3] = user->kav*eM_2_even[r][0];
735: cols[4] = 5*idx[1]+1; vals[4] = user->kav*eM_2_even[r][1];
736: cols[5] = 5*idx[2]+1; vals[5] = user->kav*eM_2_even[r][2];
738: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
740: row = 5*idx[r]+2;
741: cols[0] = 5*idx[0]+2; vals[0] = dt*eM_2_even[r][0]*ci_sum;
742: cols[1] = 5*idx[1]+2; vals[1] = dt*eM_2_even[r][1]*ci_sum;
743: cols[2] = 5*idx[2]+2; vals[2] = dt*eM_2_even[r][2]*ci_sum;
744: cols[3] = 5*idx[0]+3; vals[3] = eM_0[r][0];
745: cols[4] = 5*idx[1]+3; vals[4] = eM_0[r][1];
746: cols[5] = 5*idx[2]+3; vals[5] = eM_0[r][2];
748: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
750: row = 5*idx[r]+3;
751: cols[0] = 5*idx[0]+2; vals[0] = -1.0*eM_0[r][0];
752: cols[1] = 5*idx[1]+2; vals[1] = -1.0*eM_0[r][1];
753: cols[2] = 5*idx[2]+2; vals[2] = -1.0*eM_0[r][2];
754: cols[3] = 5*idx[0]+3; vals[3] = user->kai*eM_2_even[r][0];
755: cols[4] = 5*idx[1]+3; vals[4] = user->kai*eM_2_even[r][1];
756: cols[5] = 5*idx[2]+3; vals[5] = user->kai*eM_2_even[r][2];
758: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
760: row = 5*idx[r]+4;
761: /*
762: cols3[0] = 5*idx[0]+4; vals3[0] = eM_0[r][0]/dt + user->L*user->kaeta*dt*eM_2_even[r][0];
763: cols3[1] = 5*idx[1]+4; vals3[1] = eM_0[r][1]/dt + user->L*user->kaeta*dt*eM_2_even[r][1];
764: cols3[2] = 5*idx[2]+4; vals3[2] = eM_0[r][2]/dt + user->L*user->kaeta*dt*eM_2_even[r][2];
765: */
766: cols3[0] = 5*idx[0]+4; vals3[0] = eM_0[r][0]/dt + user->L*user->kaeta*eM_2_even[r][0];
767: cols3[1] = 5*idx[1]+4; vals3[1] = eM_0[r][1]/dt + user->L*user->kaeta*eM_2_even[r][1];
768: cols3[2] = 5*idx[2]+4; vals3[2] = eM_0[r][2]/dt + user->L*user->kaeta*eM_2_even[r][2];
770: MatSetValuesLocal(M,1,&row,3,cols3,vals3,ADD_VALUES);
772: }
774: }
775: }
777: PetscFree(nodes);
778: PetscFree(connect);
780: VecRestoreArray(user->cv,&cv_p);
781: VecRestoreArray(user->ci,&ci_p);
782: MatAssemblyBegin(M_0,MAT_FINAL_ASSEMBLY);
783: MatAssemblyEnd(M_0,MAT_FINAL_ASSEMBLY);
785: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
786: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
788: DMDARestoreElements(user->da1,&nele,&nen,&ele);
789: return(0);
790: }
795: PetscErrorCode UpdateMatrices(AppCtx *user)
796: {
798: PetscInt i,j,n,Mda,Nda;
800: PetscInt idx[3],*nodes,*connect,k;
801: PetscInt ld,rd,lu,ru;
802: PetscScalar eM_2_odd[3][3],eM_2_even[3][3],h,dt=user->dt;
803: Mat M=user->M;
804: PetscScalar *cv_p,*ci_p,cv_sum,ci_sum;
807: /* Create the mass matrix M_0 */
808: MatGetLocalSize(M,&n,NULL);
809: VecGetArray(user->cv,&cv_p);
810: VecGetArray(user->ci,&ci_p);
811: DMDAGetInfo(user->da1,NULL,&Mda,&Nda,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
813: PetscMalloc1((Mda+1)*(Nda+1),&nodes);
814: PetscMalloc1(Mda*Nda*2*3,&connect);
816: h = 100.0/Mda;
818: for (j=0; j < Nda; j++) {
819: for (i=0; i < Mda; i++) nodes[j*(Mda+1)+i] = j*Mda+i;
820: nodes[j*(Mda+1)+Mda] = j*Mda;
821: }
822: for (i=0; i < Mda; i++) nodes[Nda*(Mda+1)+i]=i;
823: nodes[Nda*(Mda+1)+Mda]=0;
826: k = 0;
827: for (j=0; j<Nda; j++) {
828: for (i=0; i<Mda; i++) {
829: ld = nodes[j*(Mda+1)+i];
830: rd = nodes[(j+1)*(Mda+1)+i];
831: ru = nodes[(j+1)*(Mda+1)+i+1];
832: lu = nodes[j*(Mda+1)+i+1];
834: connect[k*6] = ld;
835: connect[k*6+1] = lu;
836: connect[k*6+2] = rd;
837: connect[k*6+3] = lu;
838: connect[k*6+4] = ru;
839: connect[k*6+5] = rd;
841: k = k+1;
842: }
843: }
845: for (k=0; k < Mda*Nda*2; k++) {
846: idx[0] = connect[k*3];
847: idx[1] = connect[k*3+1];
848: idx[2] = connect[k*3+2];
850: PetscInt r,row,cols[3];
851: PetscScalar vals[3];
852: for (r=0; r<3; r++) {
853: row = 5*idx[r];
854: cols[0] = 5*idx[0]; vals[0] = 0.0;
855: cols[1] = 5*idx[1]; vals[1] = 0.0;
856: cols[2] = 5*idx[2]; vals[2] = 0.0;
858: /* Insert values in matrix M for 1st dof */
859: MatSetValuesLocal(M,1,&row,3,cols,vals,INSERT_VALUES);
861: row = 5*idx[r]+2;
862: cols[0] = 5*idx[0]+2; vals[0] = 0.0;
863: cols[1] = 5*idx[1]+2; vals[1] = 0.0;
864: cols[2] = 5*idx[2]+2; vals[2] = 0.0;
866: /* Insert values in matrix M for 3nd dof */
867: MatSetValuesLocal(M,1,&row,3,cols,vals,INSERT_VALUES);
868: }
869: }
871: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
872: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
874: eM_2_odd[0][0] = 1.0;
875: eM_2_odd[1][1] = eM_2_odd[2][2] = 0.5;
876: eM_2_odd[0][1] = eM_2_odd[0][2] = eM_2_odd[1][0]= eM_2_odd[2][0] = -0.5;
877: eM_2_odd[1][2] = eM_2_odd[2][1] = 0.0;
879: eM_2_even[1][1] = 1.0;
880: eM_2_even[0][0] = eM_2_even[2][2] = 0.5;
881: eM_2_even[0][1] = eM_2_even[1][0] = eM_2_even[1][2] = eM_2_even[2][1] = -0.5;
882: eM_2_even[0][2] = eM_2_even[2][0] = 0.0;
885: /* Get local element info */
886: for (k=0; k < Mda*Nda*2; k++) {
887: idx[0] = connect[k*3];
888: idx[1] = connect[k*3+1];
889: idx[2] = connect[k*3+2];
891: PetscInt row,cols[3],r;
892: PetscScalar vals[3];
894: for (r=0; r<3; r++) {
896: /* cv_sum = (1.0e-3+cv_p[idx[0]] + cv_p[idx[1]] + cv_p[idx[2]])*user->Dv/(3.0*user->kBT); */
897: /* ci_sum = (1.0e-3+ci_p[idx[0]] + ci_p[idx[1]] + ci_p[idx[2]])*user->Di/(3.0*user->kBT); */
898: cv_sum = .0000069*user->Dv/(user->kBT);
899: ci_sum = .0000069*user->Di/user->kBT;
901: if (k%2 == 0) { /* odd triangle */
903: row = 5*idx[r];
904: cols[0] = 5*idx[0]; vals[0] = dt*eM_2_odd[r][0]*cv_sum;
905: cols[1] = 5*idx[1]; vals[1] = dt*eM_2_odd[r][1]*cv_sum;
906: cols[2] = 5*idx[2]; vals[2] = dt*eM_2_odd[r][2]*cv_sum;
908: /* Insert values in matrix M for 1st dof */
909: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
911: row = 5*idx[r]+2;
912: cols[0] = 5*idx[0]+2; vals[0] = dt*eM_2_odd[r][0]*ci_sum;
913: cols[1] = 5*idx[1]+2; vals[1] = dt*eM_2_odd[r][1]*ci_sum;
914: cols[2] = 5*idx[2]+2; vals[2] = dt*eM_2_odd[r][2]*ci_sum;
916: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
918: } else {
919: row = 5*idx[r];
920: cols[0] = 5*idx[0]; vals[0] = dt*eM_2_even[r][0]*cv_sum;
921: cols[1] = 5*idx[1]; vals[1] = dt*eM_2_even[r][1]*cv_sum;
922: cols[2] = 5*idx[2]; vals[2] = dt*eM_2_even[r][2]*cv_sum;
924: /* Insert values in matrix M for 1st dof */
925: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
927: row = 5*idx[r]+2;
928: cols[0] = 5*idx[0]+2; vals[0] = dt*eM_2_even[r][0]*ci_sum;
929: cols[1] = 5*idx[1]+2; vals[1] = dt*eM_2_even[r][1]*ci_sum;
930: cols[2] = 5*idx[2]+2; vals[2] = dt*eM_2_even[r][2]*ci_sum;
931: /* Insert values in matrix M for 3nd dof */
932: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
933: }
934: }
935: }
937: PetscFree(nodes);
938: PetscFree(connect);
939: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
940: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
941: VecRestoreArray(user->cv,&cv_p);
942: VecRestoreArray(user->ci,&ci_p);
943: return(0);
944: }