Actual source code: ex23.c
petsc-3.4.5 2014-06-29
1: static char help[] = "Cahn-Hilliard-2d problem 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: -gamma <gamma>\n\
7: -theta_c <theta_c>\n\
8: -implicit <0,1> treat theta_c*M_0*u explicitly/implicitly\n\n";
10: /*
11: Run with for example: -pc_type mg -pc_mg_galerkin -T .01 -da_grid_x 65 -da_grid_y 65 -pc_mg_levels 4 -ksp_type fgmres -snes_atol 1.e-14 -mat_no_inode
12: */
14: #include petscts.h
15: #include petscdmda.h
17: typedef struct {
18: PetscReal dt,T; /* Time step and end time */
19: DM da;
20: Mat M; /* Mass matrix */
21: Mat S; /* stiffness matrix */
22: Mat M_0;
23: Vec q,u,work1;
24: PetscScalar gamma,theta_c; /* physics parameters */
25: PetscReal xmin,xmax,ymin,ymax;
26: PetscBool tsmonitor;
27: PetscInt implicit; /* Evaluate theta_c*Mo*u impliicitly or explicitly */
28: } AppCtx;
30: PetscErrorCode GetParams(AppCtx*);
31: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
32: PetscErrorCode SetUpMatrices(AppCtx*);
33: PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
34: PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat*,Mat*,MatStructure*,void*);
35: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
36: PetscErrorCode Update_q(TS);
37: PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);
41: int main(int argc, char **argv)
42: {
44: Vec x,r; /* Solution and residual vectors */
45: TS ts; /* Timestepping object */
46: AppCtx user; /* Application context */
47: Vec xl,xu; /* Upper and lower bounds on variables */
48: Mat J;
50: PetscInitialize(&argc,&argv, (char*)0, help);
52: /* Get physics and time parameters */
53: GetParams(&user);
54: /* Create a 2D DA with dof = 2 */
55: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&user.da);
56: /* Set Element type (triangular) */
57: DMDASetElementType(user.da,DMDA_ELEMENT_P1);
59: /* Set x and y coordinates */
60: DMDASetUniformCoordinates(user.da,user.xmin,user.xmax,user.ymin,user.ymax,0.0,1.0);
62: /* Get global vector x from DM and duplicate vectors r,xl,xu */
63: DMCreateGlobalVector(user.da,&x);
64: VecDuplicate(x,&r);
65: VecDuplicate(x,&xl);
66: VecDuplicate(x,&xu);
67: if (!user.implicit) {
68: VecDuplicate(x,&user.q);
69: }
71: /* Get mass,stiffness, and jacobian matrix structure from the da */
72: DMCreateMatrix(user.da,MATAIJ,&user.M);
73: DMCreateMatrix(user.da,MATAIJ,&user.S);
74: DMCreateMatrix(user.da,MATAIJ,&J);
75: /* Form the mass,stiffness matrices and matrix M_0 */
76: SetUpMatrices(&user);
78: /* Create timestepping solver context */
79: TSCreate(PETSC_COMM_WORLD,&ts);
80: TSSetProblemType(ts,TS_NONLINEAR);
83: /* Set Function evaluation and jacobian evaluation routines */
84: TSSetIFunction(ts,r,FormIFunction,(void*)&user);
85: TSSetIJacobian(ts,J,J,FormIJacobian,(void*)&user);
86: TSSetType(ts,TSBEULER);
87: /* TSThetaSetTheta(ts,1.0);*/ /* Explicit Euler */
88: TSSetDuration(ts,PETSC_DEFAULT,user.T);
89: /* Set lower and upper bound vectors */
90: SetVariableBounds(user.da,xl,xu);
91: TSVISetVariableBounds(ts,xl,xu);
93: SetInitialGuess(x,&user);
95: if (!user.implicit) {
96: TSSetApplicationContext(ts,&user);
97: TSSetSolution(ts,x);
98: Update_q(ts);
99: TSSetPostStep(ts,Update_q);
100: }
102: if (user.tsmonitor) {
103: TSMonitorSet(ts,Monitor,&user,NULL);
104: }
106: TSSetInitialTimeStep(ts,0.0,user.dt);
108: TSSetFromOptions(ts);
111: /* Run the timestepping solver */
112: TSSolve(ts,x);
114: VecDestroy(&x);
115: VecDestroy(&r);
116: VecDestroy(&xl);
117: VecDestroy(&xu);
118: if (!user.implicit) {
119: VecDestroy(&user.q);
120: VecDestroy(&user.u);
121: VecDestroy(&user.work1);
122: MatDestroy(&user.M_0);
123: }
124: MatDestroy(&user.M);
125: MatDestroy(&user.S);
126: MatDestroy(&J);
127: DMDestroy(&user.da);
128: TSDestroy(&ts);
129: PetscFinalize();
130: return 0;
131: }
135: PetscErrorCode Update_q(TS ts)
136: {
138: AppCtx *user;
139: Vec x;
140: PetscScalar *q_arr,*w_arr;
141: PetscInt i,n;
142: PetscScalar scale;
145: TSGetApplicationContext(ts,&user);
146: TSGetSolution(ts,&x);
147: VecStrideGather(x,1,user->u,INSERT_VALUES);
148: MatMult(user->M_0,user->u,user->work1);
149: scale = -(1.0 - user->implicit)*user->theta_c;
150: VecScale(user->work1,scale);
151: VecGetLocalSize(user->u,&n);
152: VecGetArray(user->q,&q_arr);
153: VecGetArray(user->work1,&w_arr);
154: for (i=0; i<n; i++) q_arr[2*i+1] = w_arr[i];
155: VecRestoreArray(user->q,&q_arr);
156: VecRestoreArray(user->work1,&w_arr);
157: return(0);
158: }
162: PetscErrorCode SetInitialGuess(Vec X,AppCtx *user)
163: {
165: PetscScalar *x;
166: PetscInt n,i;
167: PetscRandom rand;
168: PetscScalar value;
171: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
172: PetscRandomSetFromOptions(rand);
174: VecGetLocalSize(X,&n);
175: VecGetArray(X,&x);
176: /* Set initial guess, only set value for 2nd dof */
177: for (i=0; i<n/2; i++) {
178: PetscRandomGetValue(rand,&value);
179: x[2*i+1] = -0.4 + 0.05*value*(value - 0.5);
180: }
181: VecRestoreArray(X,&x);
183: PetscRandomDestroy(&rand);
184: return(0);
185: }
187: static void Gausspoints(PetscScalar *xx,PetscScalar *yy,PetscScalar *w,PetscScalar *x,PetscScalar *y)
188: {
190: xx[0] = 2.0/3.0*x[0] + 1.0/6.0*x[1] + 1.0/6.0*x[2];
191: xx[1] = 1.0/6.0*x[0] + 2.0/3.0*x[1] + 1.0/6.0*x[2];
192: xx[2] = 1.0/6.0*x[0] + 1.0/6.0*x[1] + 2.0/3.0*x[2];
194: yy[0] = 2.0/3.0*y[0] + 1.0/6.0*y[1] + 1.0/6.0*y[2];
195: yy[1] = 1.0/6.0*y[0] + 2.0/3.0*y[1] + 1.0/6.0*y[2];
196: yy[2] = 1.0/6.0*y[0] + 1.0/6.0*y[1] + 2.0/3.0*y[2];
198: *w = PetscAbsScalar(x[0]*(y[2]-y[1]) + x[2]*(y[1]-y[0]) + x[1]*(y[0]-y[2]))/6.0;
200: }
202: static void ShapefunctionsT3(PetscScalar *phi,PetscScalar phider[][2],PetscScalar xx,PetscScalar yy,PetscScalar *x,PetscScalar *y)
203: {
204: PetscScalar area,a1,a2,a3,b1,b2,b3,c1,c2,c3,pp;
206: /* Area of the triangle */
207: area = 1.0/2.0*PetscAbsScalar(x[0]*(y[2]-y[1]) + x[2]*(y[1]-y[0]) + x[1]*(y[0]-y[2]));
209: a1 = x[1]*y[2]-x[2]*y[1]; a2 = x[2]*y[0]-x[0]*y[2]; a3 = x[0]*y[1]-x[1]*y[0];
210: b1 = y[1]-y[2]; b2 = y[2]-y[0]; b3 = y[0]-y[1];
211: c1 = x[2]-x[1]; c2 = x[0]-x[2]; c3 = x[1]-x[0];
212: pp = 1.0/(2.0*area);
214: /* shape functions */
215: phi[0] = pp*(a1 + b1*xx + c1*yy);
216: phi[1] = pp*(a2 + b2*xx + c2*yy);
217: phi[2] = pp*(a3 + b3*xx + c3*yy);
219: /* shape functions derivatives */
220: phider[0][0] = pp*b1; phider[0][1] = pp*c1;
221: phider[1][0] = pp*b2; phider[1][1] = pp*c2;
222: phider[2][0] = pp*b3; phider[2][1] = pp*c3;
224: }
228: PetscErrorCode FormIFunction(TS ts,PetscReal t, Vec X,Vec Xdot,Vec F,void *ctx)
229: {
231: AppCtx *user=(AppCtx*)ctx;
234: MatMult(user->M,Xdot,F);
235: MatMultAdd(user->S,X,F,F);
236: if (!user->implicit) {
237: VecAXPY(F,1.0,user->q);
238: }
239: return(0);
240: }
244: PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat *J,Mat *B,MatStructure *flg,void *ctx)
245: {
246: PetscErrorCode ierr;
247: AppCtx *user =(AppCtx*)ctx;
248: static PetscBool copied = PETSC_FALSE;
251: /* for active set method the matrix does not get changed, so do not need to copy each time,
252: if the active set remains the same for several solves the preconditioner does not need to be rebuilt*/
253: *flg = SAME_PRECONDITIONER;
254: if (!copied) {
255: MatCopy(user->S,*J,*flg);
256: MatAXPY(*J,a,user->M,*flg);
257: copied = PETSC_TRUE;
258: }
259: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
260: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
261: /* MatView(*J,0); */
262: /* SETERRQ(PETSC_COMM_WORLD,1,"Stopped here\n"); */
263: return(0);
264: }
268: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
269: {
271: PetscScalar ***l,***u;
272: PetscInt xs,xm,ys,ym;
273: PetscInt j,i;
276: DMDAVecGetArrayDOF(da,xl,&l);
277: DMDAVecGetArrayDOF(da,xu,&u);
279: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
281: for (j=ys; j < ys+ym; j++) {
282: for (i=xs; i < xs+xm; i++) {
283: l[j][i][0] = -SNES_VI_INF;
284: l[j][i][1] = -1.0;
285: u[j][i][0] = SNES_VI_INF;
286: u[j][i][1] = 1.0;
287: }
288: }
290: DMDAVecRestoreArrayDOF(da,xl,&l);
291: DMDAVecRestoreArrayDOF(da,xu,&u);
292: return(0);
293: }
297: PetscErrorCode GetParams(AppCtx *user)
298: {
300: PetscBool flg;
303: /* Set default parameters */
304: user->tsmonitor = PETSC_FALSE;
305: user->xmin = 0.0; user->xmax = 1.0;
306: user->ymin = 0.0; user->ymax = 1.0;
307: user->T = 0.2; user->dt = 0.0001;
308: user->gamma = 3.2E-4; user->theta_c = 1;
309: user->implicit = 0;
311: PetscOptionsGetBool(NULL,"-monitor",&user->tsmonitor,NULL);
312: PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);
313: PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);
314: PetscOptionsGetReal(NULL,"-ymin",&user->ymin,&flg);
315: PetscOptionsGetReal(NULL,"-ymax",&user->ymax,&flg);
316: PetscOptionsGetScalar(NULL,"-gamma",&user->gamma,&flg);
317: PetscOptionsGetScalar(NULL,"-theta_c",&user->theta_c,&flg);
318: PetscOptionsGetInt(NULL,"-implicit",&user->implicit,&flg);
319: return(0);
320: }
324: PetscErrorCode SetUpMatrices(AppCtx *user)
325: {
326: PetscErrorCode ierr;
327: PetscInt nele,nen,i;
328: const PetscInt *ele;
329: Vec coords;
330: const PetscScalar *_coords;
331: PetscScalar x[3],y[3],xx[3],yy[3],w;
332: PetscInt idx[3];
333: PetscScalar phi[3],phider[3][2];
334: PetscScalar eM_0[3][3],eM_2[3][3];
335: Mat M = user->M;
336: Mat S = user->S;
337: PetscScalar gamma = user->gamma,theta_c=user->theta_c;
338: PetscInt implicit = user->implicit;
341: /* Get ghosted coordinates */
342: DMGetCoordinatesLocal(user->da,&coords);
343: VecGetArrayRead(coords,&_coords);
345: /* Get local element info */
346: DMDAGetElements(user->da,&nele,&nen,&ele);
347: for (i=0; i < nele; i++) {
348: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
349: x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
350: x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
351: x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
353: PetscMemzero(xx,3*sizeof(PetscScalar));
354: PetscMemzero(yy,3*sizeof(PetscScalar));
355: Gausspoints(xx,yy,&w,x,y);
357: eM_0[0][0]=eM_0[0][1]=eM_0[0][2]=0.0;
358: eM_0[1][0]=eM_0[1][1]=eM_0[1][2]=0.0;
359: eM_0[2][0]=eM_0[2][1]=eM_0[2][2]=0.0;
360: eM_2[0][0]=eM_2[0][1]=eM_2[0][2]=0.0;
361: eM_2[1][0]=eM_2[1][1]=eM_2[1][2]=0.0;
362: eM_2[2][0]=eM_2[2][1]=eM_2[2][2]=0.0;
364: PetscInt m;
365: for (m=0; m<3; m++) {
366: PetscMemzero(phi,3*sizeof(PetscScalar));
367: phider[0][0]=phider[0][1]=0.0;
368: phider[1][0]=phider[1][1]=0.0;
369: phider[2][0]=phider[2][1]=0.0;
371: ShapefunctionsT3(phi,phider,xx[m],yy[m],x,y);
373: PetscInt j,k;
374: for (j=0;j<3;j++) {
375: for (k=0;k<3;k++) {
376: eM_0[k][j] += phi[j]*phi[k]*w;
377: eM_2[k][j] += phider[j][0]*phider[k][0]*w + phider[j][1]*phider[k][1]*w;
378: }
379: }
380: }
381: PetscInt row,cols[6],r;
382: PetscScalar vals[6];
384: for (r=0;r<3;r++) {
385: row = 2*idx[r];
386: /* Insert values in the mass matrix */
387: cols[0] = 2*idx[0]+1; vals[0] = eM_0[r][0];
388: cols[1] = 2*idx[1]+1; vals[1] = eM_0[r][1];
389: cols[2] = 2*idx[2]+1; vals[2] = eM_0[r][2];
391: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
393: /* Insert values in the stiffness matrix */
394: cols[0] = 2*idx[0]; vals[0] = eM_2[r][0];
395: cols[1] = 2*idx[1]; vals[1] = eM_2[r][1];
396: cols[2] = 2*idx[2]; vals[2] = eM_2[r][2];
398: MatSetValuesLocal(S,1,&row,3,cols,vals,ADD_VALUES);
400: row = 2*idx[r]+1;
401: cols[0] = 2*idx[0]; vals[0] = -eM_0[r][0];
402: cols[1] = 2*idx[0]+1; vals[1] = gamma*eM_2[r][0]-implicit*theta_c*eM_0[r][0];
403: cols[2] = 2*idx[1]; vals[2] = -eM_0[r][1];
404: cols[3] = 2*idx[1]+1; vals[3] = gamma*eM_2[r][1]-implicit*theta_c*eM_0[r][1];
405: cols[4] = 2*idx[2]; vals[4] = -eM_0[r][2];
406: cols[5] = 2*idx[2]+1; vals[5] = gamma*eM_2[r][2]-implicit*theta_c*eM_2[r][2];
408: /* Insert values in matrix M for 2nd dof */
409: MatSetValuesLocal(S,1,&row,6,cols,vals,ADD_VALUES);
410: }
411: }
413: DMDARestoreElements(user->da,&nele,&nen,&ele);
414: VecRestoreArrayRead(coords,&_coords);
416: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
417: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
419: MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY);
420: MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY);
422: if (!implicit) {
423: /* Create ISs to extract matrix M_0 from M */
424: PetscInt n,rstart;
425: IS isrow,iscol;
427: MatGetLocalSize(M,&n,NULL);
428: MatGetOwnershipRange(M,&rstart,NULL);
429: ISCreateStride(PETSC_COMM_WORLD,n/2,rstart,2,&isrow);
430: ISCreateStride(PETSC_COMM_WORLD,n/2,rstart+1,2,&iscol);
432: /* Extract M_0 from M */
433: MatGetSubMatrix(M,isrow,iscol,MAT_INITIAL_MATRIX,&user->M_0);
435: VecCreate(PETSC_COMM_WORLD,&user->u);
436: VecSetSizes(user->u,n/2,PETSC_DECIDE);
437: VecSetFromOptions(user->u);
438: VecDuplicate(user->u,&user->work1);
440: ISDestroy(&isrow);
441: ISDestroy(&iscol);
442: }
443: return(0);
444: }
448: PetscErrorCode Monitor(TS ts,PetscInt steps,PetscReal time,Vec x,void *mctx)
449: {
453: PetscPrintf(PETSC_COMM_WORLD,"Solution vector at t = %5.4f\n",time,steps);
454: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
455: return(0);
456: }