Actual source code: ex54.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\n";
11: /*
12: 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
13: */
15: #include petscsnes.h
16: #include petscdmda.h
18: typedef struct{
19: PetscReal dt,T; /* Time step and end time */
20: DM da;
21: Mat M; /* Jacobian 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: }AppCtx;
29: PetscErrorCode GetParams(AppCtx*);
30: PetscErrorCode SetVariableBounds(DM,Vec,Vec);
31: PetscErrorCode SetUpMatrices(AppCtx*);
32: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
33: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
34: PetscErrorCode SetInitialGuess(Vec,AppCtx*);
35: PetscErrorCode Update_q(Vec,Vec,Mat,AppCtx*);
36: PetscLogEvent event_update_q;
40: int main(int argc, char **argv)
41: {
43: Vec x,r; /* Solution and residual vectors */
44: SNES snes; /* Nonlinear solver context */
45: AppCtx user; /* Application context */
46: Vec xl,xu; /* Upper and lower bounds on variables */
47: Mat J;
48: PetscReal t=0.0;
49: PETSC_UNUSED PetscLogStage stage_timestep;
50: PetscInt its;
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: VecDuplicate(x,&user.q);
71: /* Get Jacobian matrix structure from the da */
72: DMCreateMatrix(user.da,MATAIJ,&user.M);
73: /* Form the jacobian matrix and M_0 */
74: SetUpMatrices(&user);
75: MatDuplicate(user.M,MAT_DO_NOT_COPY_VALUES,&J);
77: /* Create nonlinear solver context */
78: SNESCreate(PETSC_COMM_WORLD,&snes);
79: SNESSetDM(snes,user.da);
81: /* Set Function evaluation and jacobian evaluation routines */
82: SNESSetFunction(snes,r,FormFunction,(void*)&user);
83: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
85: /* Set the boundary conditions */
86: SetVariableBounds(user.da,xl,xu);
87: SNESVISetVariableBounds(snes,xl,xu);
88: SNESSetFromOptions(snes);
90: SetInitialGuess(x,&user);
91: PetscLogStageRegister("Time stepping",&stage_timestep);
92: PetscLogEventRegister("Update q",MAT_CLASSID,&event_update_q);
93: PetscLogStagePush(stage_timestep);
94: /* Begin time loop */
95: while(t < user.T) {
96: Update_q(user.q,user.u,user.M_0,&user);
97: SNESSolve(snes,PETSC_NULL,x);
98: SNESGetIterationNumber(snes,&its);
99: if (user.tsmonitor) {
100: PetscPrintf(PETSC_COMM_WORLD,"SNESVI solver converged at t = %5.4f in %d iterations\n",t,its);
101: }
102: VecStrideGather(x,1,user.u,INSERT_VALUES);
103: t = t + user.dt;
104: }
105: PetscLogStagePop();
107: VecDestroy(&x);
108: VecDestroy(&r);
109: VecDestroy(&xl);
110: VecDestroy(&xu);
111: VecDestroy(&user.q);
112: VecDestroy(&user.u);
113: VecDestroy(&user.work1);
114: MatDestroy(&user.M);
115: MatDestroy(&user.M_0);
116: MatDestroy(&J);
117: DMDestroy(&user.da);
118: SNESDestroy(&snes);
119: PetscFinalize();
120: return 0;
121: }
125: PetscErrorCode Update_q(Vec q,Vec u,Mat M_0,AppCtx *user)
126: {
128: PetscScalar *q_arr,*w_arr;
129: PetscInt i,n;
130:
132: PetscLogEventBegin(event_update_q,0,0,0,0);
133: MatMult(M_0,u,user->work1);
134: VecScale(user->work1,-1.0);
135: VecGetLocalSize(u,&n);
136: VecGetArray(q,&q_arr);
137: VecGetArray(user->work1,&w_arr);
138: for(i=0;i<n;i++) {
139: q_arr[2*i]=q_arr[2*i+1] = w_arr[i];
140: }
141: VecRestoreArray(q,&q_arr);
142: VecRestoreArray(user->work1,&w_arr);
143: PetscLogEventEnd(event_update_q,0,0,0,0);
144: return(0);
145: }
149: PetscErrorCode SetInitialGuess(Vec X,AppCtx* user)
150: {
152: PetscScalar *x,*u;
153: PetscInt n,i;
154: Vec rand;
157: /* u = -0.4 + 0.05*rand(N,1)*(rand(N,1) - 0.5) */
158: VecDuplicate(user->u,&rand);
159: VecSetRandom(rand,PETSC_NULL);
160: VecCopy(rand,user->u);
161: VecShift(rand,-0.5);
162: VecPointwiseMult(user->u,user->u,rand);
163: VecDestroy(&rand);
164: VecScale(user->u,0.05);
165: VecShift(user->u,-0.4);
166:
167: VecGetLocalSize(X,&n);
168: VecGetArray(X,&x);
169: VecGetArray(user->u,&u);
170: /* Set initial guess, only set value for 2nd dof */
171: for(i=0;i<n/2;i++) {
172: x[2*i+1] = u[i];
173: }
174: VecRestoreArray(X,&x);
175: VecRestoreArray(user->u,&u);
177: return(0);
178: }
182: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ctx)
183: {
185: AppCtx *user=(AppCtx*)ctx;
188: MatMultAdd(user->M,X,user->q,F);
189: return(0);
190: }
194: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ctx)
195: {
196: PetscErrorCode ierr;
197: AppCtx *user=(AppCtx*)ctx;
198: static PetscBool copied = PETSC_FALSE;
201: /* for active set method the matrix does not get changed, so do not need to copy each time,
202: if the active set remains the same for several solves the preconditioner does not need to be rebuilt*/
203: *flg = SAME_PRECONDITIONER;
204: if (!copied) {
205: MatCopy(user->M,*J,*flg);
206: copied = PETSC_TRUE;
207: }
208: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
209: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
210: return(0);
211: }
214: PetscErrorCode SetVariableBounds(DM da,Vec xl,Vec xu)
215: {
217: PetscScalar ***l,***u;
218: PetscInt xs,xm,ys,ym;
219: PetscInt j,i;
222: DMDAVecGetArrayDOF(da,xl,&l);
223: DMDAVecGetArrayDOF(da,xu,&u);
225: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
227: for(j=ys; j < ys+ym; j++) {
228: for(i=xs; i < xs+xm;i++) {
229: l[j][i][0] = -SNES_VI_INF;
230: l[j][i][1] = -1.0;
231: u[j][i][0] = SNES_VI_INF;
232: u[j][i][1] = 1.0;
233: }
234: }
235:
236: DMDAVecRestoreArrayDOF(da,xl,&l);
237: DMDAVecRestoreArrayDOF(da,xu,&u);
238: return(0);
239: }
243: PetscErrorCode GetParams(AppCtx* user)
244: {
246: PetscBool flg;
247:
250: /* Set default parameters */
251: user->tsmonitor = PETSC_FALSE;
252: user->xmin = 0.0; user->xmax = 1.0;
253: user->ymin = 0.0; user->ymax = 1.0;
254: user->T = 0.0002; user->dt = 0.0001;
255: user->gamma = 3.2E-4; user->theta_c = 0;
257: PetscOptionsGetBool(PETSC_NULL,"-ts_monitor",&user->tsmonitor,PETSC_NULL);
258: PetscOptionsGetReal(PETSC_NULL,"-xmin",&user->xmin,&flg);
259: PetscOptionsGetReal(PETSC_NULL,"-xmax",&user->xmax,&flg);
260: PetscOptionsGetReal(PETSC_NULL,"-ymin",&user->ymin,&flg);
261: PetscOptionsGetReal(PETSC_NULL,"-ymax",&user->ymax,&flg);
262: PetscOptionsGetReal(PETSC_NULL,"-T",&user->T,&flg);
263: PetscOptionsGetReal(PETSC_NULL,"-dt",&user->dt,&flg);
264: PetscOptionsGetScalar(PETSC_NULL,"-gamma",&user->gamma,&flg);
265: PetscOptionsGetScalar(PETSC_NULL,"-theta_c",&user->theta_c,&flg);
267: return(0);
268: }
269:
270: static void Gausspoints(PetscScalar *xx,PetscScalar *yy,PetscScalar *w,PetscScalar *x,PetscScalar *y)
271: {
273: xx[0] = 2.0/3.0*x[0] + 1.0/6.0*x[1] + 1.0/6.0*x[2];
274: xx[1] = 1.0/6.0*x[0] + 2.0/3.0*x[1] + 1.0/6.0*x[2];
275: xx[2] = 1.0/6.0*x[0] + 1.0/6.0*x[1] + 2.0/3.0*x[2];
277: yy[0] = 2.0/3.0*y[0] + 1.0/6.0*y[1] + 1.0/6.0*y[2];
278: yy[1] = 1.0/6.0*y[0] + 2.0/3.0*y[1] + 1.0/6.0*y[2];
279: yy[2] = 1.0/6.0*y[0] + 1.0/6.0*y[1] + 2.0/3.0*y[2];
281: *w = PetscAbsScalar(x[0]*(y[2]-y[1]) + x[2]*(y[1]-y[0]) + x[1]*(y[0]-y[2]))/6.0;
283: }
285: static void ShapefunctionsT3(PetscScalar *phi,PetscScalar phider[][2],PetscScalar xx,PetscScalar yy,PetscScalar *x,PetscScalar *y)
286: {
287: PetscScalar area,a1,a2,a3,b1,b2,b3,c1,c2,c3,pp;
289: /* Area of the triangle */
290: area = 1.0/2.0*PetscAbsScalar(x[0]*(y[2]-y[1]) + x[2]*(y[1]-y[0]) + x[1]*(y[0]-y[2]));
292: 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];
293: b1 = y[1]-y[2]; b2 = y[2]-y[0]; b3 = y[0]-y[1];
294: c1 = x[2]-x[1]; c2 = x[0]-x[2]; c3 = x[1]-x[0];
295: pp = 1.0/(2.0*area);
297: /* shape functions */
298: phi[0] = pp*(a1 + b1*xx + c1*yy);
299: phi[1] = pp*(a2 + b2*xx + c2*yy);
300: phi[2] = pp*(a3 + b3*xx + c3*yy);
302: /* shape functions derivatives */
303: phider[0][0] = pp*b1; phider[0][1] = pp*c1;
304: phider[1][0] = pp*b2; phider[1][1] = pp*c2;
305: phider[2][0] = pp*b3; phider[2][1] = pp*c3;
307: }
311: PetscErrorCode SetUpMatrices(AppCtx* user)
312: {
313: PetscErrorCode ierr;
314: PetscInt nele,nen,i;
315: const PetscInt *ele;
316: PetscScalar dt=user->dt;
317: Vec coords;
318: const PetscScalar *_coords;
319: PetscScalar x[3],y[3],xx[3],yy[3],w;
320: PetscInt idx[3];
321: PetscScalar phi[3],phider[3][2];
322: PetscScalar eM_0[3][3],eM_2[3][3];
323: Mat M=user->M;
324: PetscScalar gamma=user->gamma,theta_c=user->theta_c;
325: PetscInt m;
326: PetscInt j,k;
327: PetscInt row,cols[6],r;
328: PetscScalar vals[6];
329: PetscInt n,rstart;
330: IS isrow,iscol;
333: /* Get ghosted coordinates */
334: DMDAGetGhostedCoordinates(user->da,&coords);
335: VecGetArrayRead(coords,&_coords);
337: /* Get local element info */
338: DMDAGetElements(user->da,&nele,&nen,&ele);
339: for(i=0;i < nele;i++) {
340: idx[0] = ele[3*i]; idx[1] = ele[3*i+1]; idx[2] = ele[3*i+2];
341: x[0] = _coords[2*idx[0]]; y[0] = _coords[2*idx[0]+1];
342: x[1] = _coords[2*idx[1]]; y[1] = _coords[2*idx[1]+1];
343: x[2] = _coords[2*idx[2]]; y[2] = _coords[2*idx[2]+1];
344:
345: PetscMemzero(xx,3*sizeof(PetscScalar));
346: PetscMemzero(yy,3*sizeof(PetscScalar));
347: Gausspoints(xx,yy,&w,x,y);
348:
349: eM_0[0][0]=eM_0[0][1]=eM_0[0][2]=0.0;
350: eM_0[1][0]=eM_0[1][1]=eM_0[1][2]=0.0;
351: eM_0[2][0]=eM_0[2][1]=eM_0[2][2]=0.0;
352: eM_2[0][0]=eM_2[0][1]=eM_2[0][2]=0.0;
353: eM_2[1][0]=eM_2[1][1]=eM_2[1][2]=0.0;
354: eM_2[2][0]=eM_2[2][1]=eM_2[2][2]=0.0;
357: for(m=0;m<3;m++) {
358: PetscMemzero(phi,3*sizeof(PetscScalar));
359: phider[0][0]=phider[0][1]=0.0;
360: phider[1][0]=phider[1][1]=0.0;
361: phider[2][0]=phider[2][1]=0.0;
362:
363: ShapefunctionsT3(phi,phider,xx[m],yy[m],x,y);
365: for(j=0;j<3;j++) {
366: for(k=0;k<3;k++) {
367: eM_0[k][j] += phi[j]*phi[k]*w;
368: eM_2[k][j] += phider[j][0]*phider[k][0]*w + phider[j][1]*phider[k][1]*w;
369: }
370: }
371: }
373: for(r=0;r<3;r++) {
374: row = 2*idx[r];
375: cols[0] = 2*idx[0]; vals[0] = dt*eM_2[r][0];
376: cols[1] = 2*idx[0]+1; vals[1] = eM_0[r][0];
377: cols[2] = 2*idx[1]; vals[2] = dt*eM_2[r][1];
378: cols[3] = 2*idx[1]+1; vals[3] = eM_0[r][1];
379: cols[4] = 2*idx[2]; vals[4] = dt*eM_2[r][2];
380: cols[5] = 2*idx[2]+1; vals[5] = eM_0[r][2];
382: /* Insert values in matrix M for 1st dof */
383: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
384: row = 2*idx[r]+1;
385: cols[0] = 2*idx[0]; vals[0] = -eM_0[r][0];
386: cols[1] = 2*idx[0]+1; vals[1] = gamma*eM_2[r][0]-theta_c*eM_0[r][0];
387: cols[2] = 2*idx[1]; vals[2] = -eM_0[r][1];
388: cols[3] = 2*idx[1]+1; vals[3] = gamma*eM_2[r][1]-theta_c*eM_0[r][1];
389: cols[4] = 2*idx[2]; vals[4] = -eM_0[r][2];
390: cols[5] = 2*idx[2]+1; vals[5] = gamma*eM_2[r][2]-theta_c*eM_2[r][2];
392: /* Insert values in matrix M for 2nd dof */
393: MatSetValuesLocal(M,1,&row,6,cols,vals,ADD_VALUES);
394: }
395: }
397: DMDARestoreElements(user->da,&nele,&nen,&ele);
398: VecRestoreArrayRead(coords,&_coords);
400: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
401: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
403: /* Create ISs to extract matrix M_0 from M */
405: MatGetLocalSize(M,&n,PETSC_NULL);
406: MatGetOwnershipRange(M,&rstart,PETSC_NULL);
407: ISCreateStride(PETSC_COMM_WORLD,n/2,rstart,2,&isrow);
408: ISCreateStride(PETSC_COMM_WORLD,n/2,rstart+1,2,&iscol);
410: /* Extract M_0 from M */
411: MatGetSubMatrix(M,isrow,iscol,MAT_INITIAL_MATRIX,&user->M_0);
412: VecCreate(PETSC_COMM_WORLD,&user->u);
413: VecSetSizes(user->u,n/2,PETSC_DECIDE);
414: VecSetFromOptions(user->u);
415: VecDuplicate(user->u,&user->work1);
416: ISDestroy(&isrow);
417: ISDestroy(&iscol);
419: return(0);
420: }