Actual source code: ex23.c
petsc-3.3-p7 2013-05-11
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: -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\
10: -implicit <0,1> treat theta_c*M_0*u explicitly/implicitly\n\n";
12: /*
13: 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
14: */
16: #include petscts.h
17: #include petscdmda.h
19: typedef struct{
20: PetscReal dt,T; /* Time step and end time */
21: DM da;
22: Mat M; /* Mass matrix */
23: Mat S; /* stiffness matrix */
24: Mat M_0;
25: Vec q,u,work1;
26: PetscScalar gamma,theta_c; /* physics parameters */
27: PetscReal xmin,xmax,ymin,ymax;
28: PetscBool tsmonitor;
29: PetscInt implicit; /* Evaluate theta_c*Mo*u impliicitly or explicitly */
30: }AppCtx;
32: PetscErrorCode GetParams(AppCtx*);
33: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
34: PetscErrorCode SetUpMatrices(AppCtx*);
35: PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
36: PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat*,Mat*,MatStructure*,void*);
37: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
38: PetscErrorCode Update_q(TS);
39: PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);
43: int main(int argc, char **argv)
44: {
46: Vec x,r; /* Solution and residual vectors */
47: TS ts; /* Timestepping object */
48: AppCtx user; /* Application context */
49: Vec xl,xu; /* Upper and lower bounds on variables */
50: Mat J;
52: PetscInitialize(&argc,&argv, (char*)0, help);
54: /* Get physics and time parameters */
55: GetParams(&user);
56: /* Create a 2D DA with dof = 2 */
57: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,2,1,PETSC_NULL,PETSC_NULL,&user.da);
58: /* Set Element type (triangular) */
59: DMDASetElementType(user.da,DMDA_ELEMENT_P1);
61: /* Set x and y coordinates */
62: DMDASetUniformCoordinates(user.da,user.xmin,user.xmax,user.ymin,user.ymax,0.0,1.0);
63:
64: /* Get global vector x from DM and duplicate vectors r,xl,xu */
65: DMCreateGlobalVector(user.da,&x);
66: VecDuplicate(x,&r);
67: VecDuplicate(x,&xl);
68: VecDuplicate(x,&xu);
69: if(!user.implicit) {
70: VecDuplicate(x,&user.q);
71: }
73: /* Get mass,stiffness, and jacobian matrix structure from the da */
74: DMCreateMatrix(user.da,MATAIJ,&user.M);
75: DMCreateMatrix(user.da,MATAIJ,&user.S);
76: DMCreateMatrix(user.da,MATAIJ,&J);
77: /* Form the mass,stiffness matrices and matrix M_0 */
78: SetUpMatrices(&user);
80: /* Create timestepping solver context */
81: TSCreate(PETSC_COMM_WORLD,&ts);
82: TSSetProblemType(ts,TS_NONLINEAR);
85: /* Set Function evaluation and jacobian evaluation routines */
86: TSSetIFunction(ts,r,FormIFunction,(void*)&user);
87: TSSetIJacobian(ts,J,J,FormIJacobian,(void*)&user);
88: TSSetType(ts,TSBEULER);
89: // TSThetaSetTheta(ts,1.0);/* Explicit Euler */
90: TSSetDuration(ts,PETSC_DEFAULT,user.T);
91: /* Set lower and upper bound vectors */
92: SetVariableBounds(user.da,xl,xu);
93: TSVISetVariableBounds(ts,xl,xu);
95: SetInitialGuess(x,&user);
97: if(!user.implicit) {
98: TSSetApplicationContext(ts,&user);
99: TSSetSolution(ts,x);
100: Update_q(ts);
101: TSSetPostStep(ts,Update_q);
102: }
104: if(user.tsmonitor) {
105: TSMonitorSet(ts,Monitor,&user,PETSC_NULL);
106: }
108: TSSetInitialTimeStep(ts,0.0,user.dt);
110: TSSetFromOptions(ts);
113: /* Run the timestepping solver */
114: TSSolve(ts,x,PETSC_NULL);
116: VecDestroy(&x);
117: VecDestroy(&r);
118: VecDestroy(&xl);
119: VecDestroy(&xu);
120: if(!user.implicit) {
121: VecDestroy(&user.q);
122: VecDestroy(&user.u);
123: VecDestroy(&user.work1);
124: MatDestroy(&user.M_0);
125: }
126: MatDestroy(&user.M);
127: MatDestroy(&user.S);
128: MatDestroy(&J);
129: DMDestroy(&user.da);
130: TSDestroy(&ts);
131: PetscFinalize();
132: return 0;
133: }
137: PetscErrorCode Update_q(TS ts)
138: {
140: AppCtx *user;
141: Vec x;
142: PetscScalar *q_arr,*w_arr;
143: PetscInt i,n;
144: PetscScalar scale;
145:
147: TSGetApplicationContext(ts,&user);
148: TSGetSolution(ts,&x);
149: VecStrideGather(x,1,user->u,INSERT_VALUES);
150: MatMult(user->M_0,user->u,user->work1);
151: scale = -(1.0 - user->implicit)*user->theta_c;
152: VecScale(user->work1,scale);
153: VecGetLocalSize(user->u,&n);
154: VecGetArray(user->q,&q_arr);
155: VecGetArray(user->work1,&w_arr);
156: for(i=0;i<n;i++) {
157: q_arr[2*i+1] = w_arr[i];
158: }
159: VecRestoreArray(user->q,&q_arr);
160: VecRestoreArray(user->work1,&w_arr);
161: return(0);
162: }
166: PetscErrorCode SetInitialGuess(Vec X,AppCtx* user)
167: {
169: PetscScalar *x;
170: PetscInt n,i;
171: PetscRandom rand;
172: PetscScalar value;
175: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
176: PetscRandomSetFromOptions(rand);
177:
178: VecGetLocalSize(X,&n);
179: VecGetArray(X,&x);
180: /* Set initial guess, only set value for 2nd dof */
181: for(i=0;i<n/2;i++) {
182: PetscRandomGetValue(rand,&value);
183: x[2*i+1] = -0.4 + 0.05*value*(value - 0.5);
184: }
185: VecRestoreArray(X,&x);
187: PetscRandomDestroy(&rand);
189: return(0);
190: }
192: static void Gausspoints(PetscScalar *xx,PetscScalar *yy,PetscScalar *w,PetscScalar *x,PetscScalar *y)
193: {
195: xx[0] = 2.0/3.0*x[0] + 1.0/6.0*x[1] + 1.0/6.0*x[2];
196: xx[1] = 1.0/6.0*x[0] + 2.0/3.0*x[1] + 1.0/6.0*x[2];
197: xx[2] = 1.0/6.0*x[0] + 1.0/6.0*x[1] + 2.0/3.0*x[2];
199: yy[0] = 2.0/3.0*y[0] + 1.0/6.0*y[1] + 1.0/6.0*y[2];
200: yy[1] = 1.0/6.0*y[0] + 2.0/3.0*y[1] + 1.0/6.0*y[2];
201: yy[2] = 1.0/6.0*y[0] + 1.0/6.0*y[1] + 2.0/3.0*y[2];
203: *w = PetscAbsScalar(x[0]*(y[2]-y[1]) + x[2]*(y[1]-y[0]) + x[1]*(y[0]-y[2]))/6.0;
205: }
207: static void ShapefunctionsT3(PetscScalar *phi,PetscScalar phider[][2],PetscScalar xx,PetscScalar yy,PetscScalar *x,PetscScalar *y)
208: {
209: PetscScalar area,a1,a2,a3,b1,b2,b3,c1,c2,c3,pp;
211: /* Area of the triangle */
212: area = 1.0/2.0*PetscAbsScalar(x[0]*(y[2]-y[1]) + x[2]*(y[1]-y[0]) + x[1]*(y[0]-y[2]));
214: 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];
215: b1 = y[1]-y[2]; b2 = y[2]-y[0]; b3 = y[0]-y[1];
216: c1 = x[2]-x[1]; c2 = x[0]-x[2]; c3 = x[1]-x[0];
217: pp = 1.0/(2.0*area);
219: /* shape functions */
220: phi[0] = pp*(a1 + b1*xx + c1*yy);
221: phi[1] = pp*(a2 + b2*xx + c2*yy);
222: phi[2] = pp*(a3 + b3*xx + c3*yy);
224: /* shape functions derivatives */
225: phider[0][0] = pp*b1; phider[0][1] = pp*c1;
226: phider[1][0] = pp*b2; phider[1][1] = pp*c2;
227: phider[2][0] = pp*b3; phider[2][1] = pp*c3;
229: }
233: PetscErrorCode FormIFunction(TS ts,PetscReal t, Vec X,Vec Xdot,Vec F,void* ctx)
234: {
236: AppCtx *user=(AppCtx*)ctx;
239: MatMult(user->M,Xdot,F);
240: MatMultAdd(user->S,X,F,F);
241: if(!user->implicit) {
242: VecAXPY(F,1.0,user->q);
243: }
244: return(0);
245: }
249: PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat *J,Mat *B,MatStructure *flg,void *ctx)
250: {
251: PetscErrorCode ierr;
252: AppCtx *user=(AppCtx*)ctx;
253: static PetscBool copied = PETSC_FALSE;
256: /* for active set method the matrix does not get changed, so do not need to copy each time,
257: if the active set remains the same for several solves the preconditioner does not need to be rebuilt*/
258: *flg = SAME_PRECONDITIONER;
259: if (!copied) {
260: MatCopy(user->S,*J,*flg);
261: MatAXPY(*J,a,user->M,*flg);
262: copied = PETSC_TRUE;
263: }
264: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
265: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
266: // MatView(*J,0);
267: // SETERRQ(PETSC_COMM_WORLD,1,"Stopped here\n");
268: return(0);
269: }
273: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
274: {
276: PetscScalar ***l,***u;
277: PetscInt xs,xm,ys,ym;
278: PetscInt j,i;
281: DMDAVecGetArrayDOF(da,xl,&l);
282: DMDAVecGetArrayDOF(da,xu,&u);
284: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
286: for(j=ys; j < ys+ym; j++) {
287: for(i=xs; i < xs+xm;i++) {
288: l[j][i][0] = -SNES_VI_INF;
289: l[j][i][1] = -1.0;
290: u[j][i][0] = SNES_VI_INF;
291: u[j][i][1] = 1.0;
292: }
293: }
294:
295: DMDAVecRestoreArrayDOF(da,xl,&l);
296: DMDAVecRestoreArrayDOF(da,xu,&u);
297: return(0);
298: }
302: PetscErrorCode GetParams(AppCtx* user)
303: {
305: PetscBool flg;
306:
309: /* Set default parameters */
310: user->tsmonitor = PETSC_FALSE;
311: user->xmin = 0.0; user->xmax = 1.0;
312: user->ymin = 0.0; user->ymax = 1.0;
313: user->T = 0.2; user->dt = 0.0001;
314: user->gamma = 3.2E-4; user->theta_c = 1;
315: user->implicit = 0;
317: PetscOptionsGetBool(PETSC_NULL,"-ts_monitor",&user->tsmonitor,PETSC_NULL);
318: PetscOptionsGetReal(PETSC_NULL,"-xmin",&user->xmin,&flg);
319: PetscOptionsGetReal(PETSC_NULL,"-xmax",&user->xmax,&flg);
320: PetscOptionsGetReal(PETSC_NULL,"-ymin",&user->ymin,&flg);
321: PetscOptionsGetReal(PETSC_NULL,"-ymax",&user->ymax,&flg);
322: PetscOptionsGetReal(PETSC_NULL,"-T",&user->T,&flg);
323: PetscOptionsGetReal(PETSC_NULL,"-dt",&user->dt,&flg);
324: PetscOptionsGetScalar(PETSC_NULL,"-gamma",&user->gamma,&flg);
325: PetscOptionsGetScalar(PETSC_NULL,"-theta_c",&user->theta_c,&flg);
326: PetscOptionsGetInt(PETSC_NULL,"-implicit",&user->implicit,&flg);
328: return(0);
329: }
333: PetscErrorCode SetUpMatrices(AppCtx* user)
334: {
335: PetscErrorCode ierr;
336: PetscInt nele,nen,i;
337: const PetscInt *ele;
338: Vec coords;
339: const PetscScalar *_coords;
340: PetscScalar x[3],y[3],xx[3],yy[3],w;
341: PetscInt idx[3];
342: PetscScalar phi[3],phider[3][2];
343: PetscScalar eM_0[3][3],eM_2[3][3];
344: Mat M=user->M;
345: Mat S=user->S;
346: PetscScalar gamma=user->gamma,theta_c=user->theta_c;
347: PetscInt implicit = user->implicit;
350: /* Get ghosted coordinates */
351: DMDAGetGhostedCoordinates(user->da,&coords);
352: VecGetArrayRead(coords,&_coords);
354: /* Get local element info */
355: DMDAGetElements(user->da,&nele,&nen,&ele);
356: for(i=0;i < nele;i++) {
357: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
358: x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
359: x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
360: x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
361:
362: PetscMemzero(xx,3*sizeof(PetscScalar));
363: PetscMemzero(yy,3*sizeof(PetscScalar));
364: Gausspoints(xx,yy,&w,x,y);
365:
366: eM_0[0][0]=eM_0[0][1]=eM_0[0][2]=0.0;
367: eM_0[1][0]=eM_0[1][1]=eM_0[1][2]=0.0;
368: eM_0[2][0]=eM_0[2][1]=eM_0[2][2]=0.0;
369: eM_2[0][0]=eM_2[0][1]=eM_2[0][2]=0.0;
370: eM_2[1][0]=eM_2[1][1]=eM_2[1][2]=0.0;
371: eM_2[2][0]=eM_2[2][1]=eM_2[2][2]=0.0;
373: PetscInt m;
374: for(m=0;m<3;m++) {
375: PetscMemzero(phi,3*sizeof(PetscScalar));
376: phider[0][0]=phider[0][1]=0.0;
377: phider[1][0]=phider[1][1]=0.0;
378: phider[2][0]=phider[2][1]=0.0;
379:
380: ShapefunctionsT3(phi,phider,xx[m],yy[m],x,y);
382: PetscInt j,k;
383: for(j=0;j<3;j++) {
384: for(k=0;k<3;k++) {
385: eM_0[k][j] += phi[j]*phi[k]*w;
386: eM_2[k][j] += phider[j][0]*phider[k][0]*w + phider[j][1]*phider[k][1]*w;
387: }
388: }
389: }
390: PetscInt row,cols[6],r;
391: PetscScalar vals[6];
393: for(r=0;r<3;r++) {
394: row = 2*idx[r];
395: /* Insert values in the mass matrix */
396: cols[0] = 2*idx[0]+1; vals[0] = eM_0[r][0];
397: cols[1] = 2*idx[1]+1; vals[1] = eM_0[r][1];
398: cols[2] = 2*idx[2]+1; vals[2] = eM_0[r][2];
400: MatSetValuesLocal(M,1,&row,3,cols,vals,ADD_VALUES);
402: /* Insert values in the stiffness matrix */
403: cols[0] = 2*idx[0]; vals[0] = eM_2[r][0];
404: cols[1] = 2*idx[1]; vals[1] = eM_2[r][1];
405: cols[2] = 2*idx[2]; vals[2] = eM_2[r][2];
407: MatSetValuesLocal(S,1,&row,3,cols,vals,ADD_VALUES);
409: row = 2*idx[r]+1;
410: cols[0] = 2*idx[0]; vals[0] = -eM_0[r][0];
411: cols[1] = 2*idx[0]+1; vals[1] = gamma*eM_2[r][0]-implicit*theta_c*eM_0[r][0];
412: cols[2] = 2*idx[1]; vals[2] = -eM_0[r][1];
413: cols[3] = 2*idx[1]+1; vals[3] = gamma*eM_2[r][1]-implicit*theta_c*eM_0[r][1];
414: cols[4] = 2*idx[2]; vals[4] = -eM_0[r][2];
415: cols[5] = 2*idx[2]+1; vals[5] = gamma*eM_2[r][2]-implicit*theta_c*eM_2[r][2];
417: /* Insert values in matrix M for 2nd dof */
418: MatSetValuesLocal(S,1,&row,6,cols,vals,ADD_VALUES);
419: }
420: }
422: DMDARestoreElements(user->da,&nele,&nen,&ele);
423: VecRestoreArrayRead(coords,&_coords);
425: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
426: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
428: MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY);
429: MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY);
431: if(!implicit) {
432: /* Create ISs to extract matrix M_0 from M */
433: PetscInt n,rstart;
434: IS isrow,iscol;
436: MatGetLocalSize(M,&n,PETSC_NULL);
437: MatGetOwnershipRange(M,&rstart,PETSC_NULL);
438: ISCreateStride(PETSC_COMM_WORLD,n/2,rstart,2,&isrow);
439: ISCreateStride(PETSC_COMM_WORLD,n/2,rstart+1,2,&iscol);
440:
441: /* Extract M_0 from M */
442: MatGetSubMatrix(M,isrow,iscol,MAT_INITIAL_MATRIX,&user->M_0);
443:
444: VecCreate(PETSC_COMM_WORLD,&user->u);
445: VecSetSizes(user->u,n/2,PETSC_DECIDE);
446: VecSetFromOptions(user->u);
447: VecDuplicate(user->u,&user->work1);
449: ISDestroy(&isrow);
450: ISDestroy(&iscol);
451: }
453: return(0);
454: }
458: PetscErrorCode Monitor(TS ts,PetscInt steps,PetscReal time,Vec x,void* mctx)
459: {
463: PetscPrintf(PETSC_COMM_WORLD,"Solution vector at t = %5.4f\n",time,steps);
464: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
466: return(0);
467: }