Actual source code: biharmonic.c
2: static char help[] = "Solves biharmonic equation in 1d.\n";
4: /*
5: Solves the equation
7: u_t = - kappa \Delta \Delta u
8: Periodic boundary conditions
10: Evolve the biharmonic heat equation:
11: ---------------
12: ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -mymonitor
14: Evolve with the restriction that -1 <= u <= 1; i.e. as a variational inequality
15: ---------------
16: ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -mymonitor
18: u_t = kappa \Delta \Delta u + 6.*u*(u_x)^2 + (3*u^2 - 12) \Delta u
19: -1 <= u <= 1
20: Periodic boundary conditions
22: Evolve the Cahn-Hillard equations: double well Initial hump shrinks then grows
23: ---------------
24: ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 6 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -ts_monitor_draw_solution --mymonitor
26: Initial hump neither shrinks nor grows when degenerate (otherwise similar solution)
28: ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 6 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -degenerate -ts_monitor_draw_solution --mymonitor
30: ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 6 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -snes_vi_ignore_function_sign -ts_monitor_draw_solution --mymonitor
32: Evolve the Cahn-Hillard equations: double obstacle
33: ---------------
34: ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -energy 2 -snes_linesearch_monitor -ts_monitor_draw_solution --mymonitor
36: Evolve the Cahn-Hillard equations: logarithmic + double well (never shrinks and then grows)
37: ---------------
38: ./biharmonic -ts_monitor -snes_monitor -pc_type lu --snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -kappa .0001 -ts_dt 5.96046e-06 -cahn-hillard -energy 3 -snes_linesearch_monitor -theta .00000001 -ts_monitor_draw_solution --ts_max_time 1. -mymonitor
40: ./biharmonic -ts_monitor -snes_monitor -pc_type lu --snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -kappa .0001 -ts_dt 5.96046e-06 -cahn-hillard -energy 3 -snes_linesearch_monitor -theta .00000001 -ts_monitor_draw_solution --ts_max_time 1. -degenerate -mymonitor
42: Evolve the Cahn-Hillard equations: logarithmic + double obstacle (never shrinks, never grows)
43: ---------------
44: ./biharmonic -ts_monitor -snes_monitor -pc_type lu --snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -energy 4 -snes_linesearch_monitor -theta .00000001 -ts_monitor_draw_solution --mymonitor
46: */
47: #include <petscdm.h>
48: #include <petscdmda.h>
49: #include <petscts.h>
50: #include <petscdraw.h>
52: extern PetscErrorCode FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DM,Vec),MyMonitor(TS,PetscInt,PetscReal,Vec,void*),MyDestroy(void**),FormJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
53: typedef struct {PetscBool cahnhillard;PetscBool degenerate;PetscReal kappa;PetscInt energy;PetscReal tol;PetscReal theta,theta_c;PetscInt truncation;PetscBool netforce; PetscDrawViewPorts *ports;} UserCtx;
55: int main(int argc,char **argv)
56: {
57: TS ts; /* nonlinear solver */
58: Vec x,r; /* solution, residual vectors */
59: Mat J; /* Jacobian matrix */
60: PetscInt steps,Mx;
62: DM da;
63: PetscReal dt;
64: PetscBool mymonitor;
65: UserCtx ctx;
67: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: Initialize program
69: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
71: ctx.kappa = 1.0;
72: PetscOptionsGetReal(NULL,NULL,"-kappa",&ctx.kappa,NULL);
73: ctx.degenerate = PETSC_FALSE;
74: PetscOptionsGetBool(NULL,NULL,"-degenerate",&ctx.degenerate,NULL);
75: ctx.cahnhillard = PETSC_FALSE;
76: PetscOptionsGetBool(NULL,NULL,"-cahn-hillard",&ctx.cahnhillard,NULL);
77: ctx.netforce = PETSC_FALSE;
78: PetscOptionsGetBool(NULL,NULL,"-netforce",&ctx.netforce,NULL);
79: ctx.energy = 1;
80: PetscOptionsGetInt(NULL,NULL,"-energy",&ctx.energy,NULL);
81: ctx.tol = 1.0e-8;
82: PetscOptionsGetReal(NULL,NULL,"-tol",&ctx.tol,NULL);
83: ctx.theta = .001;
84: ctx.theta_c = 1.0;
85: PetscOptionsGetReal(NULL,NULL,"-theta",&ctx.theta,NULL);
86: PetscOptionsGetReal(NULL,NULL,"-theta_c",&ctx.theta_c,NULL);
87: ctx.truncation = 1;
88: PetscOptionsGetInt(NULL,NULL,"-truncation",&ctx.truncation,NULL);
89: PetscOptionsHasName(NULL,NULL,"-mymonitor",&mymonitor);
91: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: Create distributed array (DMDA) to manage parallel grid and vectors
93: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10,1,2,NULL,&da);
95: DMSetFromOptions(da);
96: DMSetUp(da);
97: DMDASetFieldName(da,0,"Biharmonic heat equation: u");
98: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,0,0,0,0,0,0);
99: dt = 1.0/(10.*ctx.kappa*Mx*Mx*Mx*Mx);
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Extract global vectors from DMDA; then duplicate for remaining
103: vectors that are the same types
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: DMCreateGlobalVector(da,&x);
106: VecDuplicate(x,&r);
108: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: Create timestepping solver context
110: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111: TSCreate(PETSC_COMM_WORLD,&ts);
112: TSSetDM(ts,da);
113: TSSetProblemType(ts,TS_NONLINEAR);
114: TSSetRHSFunction(ts,NULL,FormFunction,&ctx);
115: DMSetMatType(da,MATAIJ);
116: DMCreateMatrix(da,&J);
117: TSSetRHSJacobian(ts,J,J,FormJacobian,&ctx);
118: TSSetMaxTime(ts,.02);
119: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_INTERPOLATE);
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Create matrix data structure; set Jacobian evaluation routine
124: Set Jacobian matrix data structure and default Jacobian evaluation
125: routine. User can override with:
126: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
127: (unless user explicitly sets preconditioner)
128: -snes_mf_operator : form preconditioning matrix as set by the user,
129: but use matrix-free approx for Jacobian-vector
130: products within Newton-Krylov method
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: Customize nonlinear solver
135: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136: TSSetType(ts,TSCN);
138: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: Set initial conditions
140: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: FormInitialSolution(da,x);
142: TSSetTimeStep(ts,dt);
143: TSSetSolution(ts,x);
145: if (mymonitor) {
146: ctx.ports = NULL;
147: TSMonitorSet(ts,MyMonitor,&ctx,MyDestroy);
148: }
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: Set runtime options
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: TSSetFromOptions(ts);
155: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: Solve nonlinear system
157: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158: TSSolve(ts,x);
159: TSGetStepNumber(ts,&steps);
160: VecView(x,PETSC_VIEWER_BINARY_WORLD);
162: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: Free work space. All PETSc objects should be destroyed when they
164: are no longer needed.
165: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166: MatDestroy(&J);
167: VecDestroy(&x);
168: VecDestroy(&r);
169: TSDestroy(&ts);
170: DMDestroy(&da);
172: PetscFinalize();
173: return ierr;
174: }
175: /* ------------------------------------------------------------------- */
176: /*
177: FormFunction - Evaluates nonlinear function, F(x).
179: Input Parameters:
180: . ts - the TS context
181: . X - input vector
182: . ptr - optional user-defined context, as set by SNESSetFunction()
184: Output Parameter:
185: . F - function vector
186: */
187: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
188: {
189: DM da;
191: PetscInt i,Mx,xs,xm;
192: PetscReal hx,sx;
193: PetscScalar *x,*f,c,r,l;
194: Vec localX;
195: UserCtx *ctx = (UserCtx*)ptr;
196: PetscReal tol = ctx->tol, theta=ctx->theta,theta_c=ctx->theta_c,a,b; /* a and b are used in the cubic truncation of the log function */
199: TSGetDM(ts,&da);
200: DMGetLocalVector(da,&localX);
201: DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
203: hx = 1.0/(PetscReal)Mx; sx = 1.0/(hx*hx);
205: /*
206: Scatter ghost points to local vector,using the 2-step process
207: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
208: By placing code between these two statements, computations can be
209: done while messages are in transition.
210: */
211: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
212: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
214: /*
215: Get pointers to vector data
216: */
217: DMDAVecGetArrayRead(da,localX,&x);
218: DMDAVecGetArray(da,F,&f);
220: /*
221: Get local grid boundaries
222: */
223: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
225: /*
226: Compute function over the locally owned part of the grid
227: */
228: for (i=xs; i<xs+xm; i++) {
229: if (ctx->degenerate) {
230: c = (1. - x[i]*x[i])*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
231: r = (1. - x[i+1]*x[i+1])*(x[i] + x[i+2] - 2.0*x[i+1])*sx;
232: l = (1. - x[i-1]*x[i-1])*(x[i-2] + x[i] - 2.0*x[i-1])*sx;
233: } else {
234: c = (x[i-1] + x[i+1] - 2.0*x[i])*sx;
235: r = (x[i] + x[i+2] - 2.0*x[i+1])*sx;
236: l = (x[i-2] + x[i] - 2.0*x[i-1])*sx;
237: }
238: f[i] = -ctx->kappa*(l + r - 2.0*c)*sx;
239: if (ctx->cahnhillard) {
240: switch (ctx->energy) {
241: case 1: /* double well */
242: f[i] += 6.*.25*x[i]*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (3.*x[i]*x[i] - 1.)*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
243: break;
244: case 2: /* double obstacle */
245: f[i] += -(x[i-1] + x[i+1] - 2.0*x[i])*sx;
246: break;
247: case 3: /* logarithmic + double well */
248: f[i] += 6.*.25*x[i]*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (3.*x[i]*x[i] - 1.)*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
249: if (ctx->truncation==2) { /* log function with approximated with a quadratic polynomial outside -1.0+2*tol, 1.0-2*tol */
250: if (PetscRealPart(x[i]) < -1.0 + 2.0*tol) f[i] += (.25*theta/(tol-tol*tol))*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
251: else if (PetscRealPart(x[i]) > 1.0 - 2.0*tol) f[i] += (.25*theta/(tol-tol*tol))*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
252: else f[i] += 2.0*theta*x[i]/((1.0-x[i]*x[i])*(1.0-x[i]*x[i]))*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (theta/(1.0-x[i]*x[i]))*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
253: } else { /* log function is approximated with a cubic polynomial outside -1.0+2*tol, 1.0-2*tol */
254: a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol));
255: b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol);
256: if (PetscRealPart(x[i]) < -1.0 + 2.0*tol) f[i] += -1.0*a*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (-1.0*a*x[i] + b)*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
257: else if (PetscRealPart(x[i]) > 1.0 - 2.0*tol) f[i] += 1.0*a*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + ( a*x[i] + b)*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
258: else f[i] += 2.0*theta*x[i]/((1.0-x[i]*x[i])*(1.0-x[i]*x[i]))*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (theta/(1.0-x[i]*x[i]))*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
259: }
260: break;
261: case 4: /* logarithmic + double obstacle */
262: f[i] += -theta_c*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
263: if (ctx->truncation==2) { /* quadratic */
264: if (PetscRealPart(x[i]) < -1.0 + 2.0*tol) f[i] += (.25*theta/(tol-tol*tol))*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
265: else if (PetscRealPart(x[i]) > 1.0 - 2.0*tol) f[i] += (.25*theta/(tol-tol*tol))*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
266: else f[i] += 2.0*theta*x[i]/((1.0-x[i]*x[i])*(1.0-x[i]*x[i]))*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (theta/(1.0-x[i]*x[i]))*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
267: } else { /* cubic */
268: a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol));
269: b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol);
270: if (PetscRealPart(x[i]) < -1.0 + 2.0*tol) f[i] += -1.0*a*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (-1.0*a*x[i] + b)*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
271: else if (PetscRealPart(x[i]) > 1.0 - 2.0*tol) f[i] += 1.0*a*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + ( a*x[i] + b)*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
272: else f[i] += 2.0*theta*x[i]/((1.0-x[i]*x[i])*(1.0-x[i]*x[i]))*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (theta/(1.0-x[i]*x[i]))*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
273: }
274: break;
275: }
276: }
278: }
280: /*
281: Restore vectors
282: */
283: DMDAVecRestoreArrayRead(da,localX,&x);
284: DMDAVecRestoreArray(da,F,&f);
285: DMRestoreLocalVector(da,&localX);
286: return(0);
287: }
289: /* ------------------------------------------------------------------- */
290: /*
291: FormJacobian - Evaluates nonlinear function's Jacobian
293: */
294: PetscErrorCode FormJacobian(TS ts,PetscReal ftime,Vec X,Mat A,Mat B,void *ptr)
295: {
296: DM da;
298: PetscInt i,Mx,xs,xm;
299: MatStencil row,cols[5];
300: PetscReal hx,sx;
301: PetscScalar *x,vals[5];
302: Vec localX;
303: UserCtx *ctx = (UserCtx*)ptr;
306: TSGetDM(ts,&da);
307: DMGetLocalVector(da,&localX);
308: DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
310: hx = 1.0/(PetscReal)Mx; sx = 1.0/(hx*hx);
312: /*
313: Scatter ghost points to local vector,using the 2-step process
314: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
315: By placing code between these two statements, computations can be
316: done while messages are in transition.
317: */
318: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
319: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
321: /*
322: Get pointers to vector data
323: */
324: DMDAVecGetArrayRead(da,localX,&x);
326: /*
327: Get local grid boundaries
328: */
329: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
331: /*
332: Compute function over the locally owned part of the grid
333: */
334: for (i=xs; i<xs+xm; i++) {
335: row.i = i;
336: if (ctx->degenerate) {
337: /*PetscScalar c,r,l;
338: c = (1. - x[i]*x[i])*(x[i-1] + x[i+1] - 2.0*x[i])*sx;
339: r = (1. - x[i+1]*x[i+1])*(x[i] + x[i+2] - 2.0*x[i+1])*sx;
340: l = (1. - x[i-1]*x[i-1])*(x[i-2] + x[i] - 2.0*x[i-1])*sx; */
341: } else {
342: cols[0].i = i - 2; vals[0] = -ctx->kappa*sx*sx;
343: cols[1].i = i - 1; vals[1] = 4.0*ctx->kappa*sx*sx;
344: cols[2].i = i ; vals[2] = -6.0*ctx->kappa*sx*sx;
345: cols[3].i = i + 1; vals[3] = 4.0*ctx->kappa*sx*sx;
346: cols[4].i = i + 2; vals[4] = -ctx->kappa*sx*sx;
347: }
348: MatSetValuesStencil(B,1,&row,5,cols,vals,INSERT_VALUES);
350: if (ctx->cahnhillard) {
351: switch (ctx->energy) {
352: case 1: /* double well */
353: /* f[i] += 6.*.25*x[i]*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (3.*x[i]*x[i] - 1.)*(x[i-1] + x[i+1] - 2.0*x[i])*sx; */
354: break;
355: case 2: /* double obstacle */
356: /* f[i] += -(x[i-1] + x[i+1] - 2.0*x[i])*sx; */
357: break;
358: case 3: /* logarithmic + double well */
359: break;
360: case 4: /* logarithmic + double obstacle */
361: break;
362: }
363: }
365: }
367: /*
368: Restore vectors
369: */
370: DMDAVecRestoreArrayRead(da,localX,&x);
371: DMRestoreLocalVector(da,&localX);
372: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
373: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
374: if (A != B) {
375: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
376: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
377: }
378: return(0);
379: }
380: /* ------------------------------------------------------------------- */
381: PetscErrorCode FormInitialSolution(DM da,Vec U)
382: {
383: PetscErrorCode ierr;
384: PetscInt i,xs,xm,Mx,N,scale;
385: PetscScalar *u;
386: PetscReal r,hx,x;
387: const PetscScalar *f;
388: Vec finesolution;
389: PetscViewer viewer;
392: DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
394: hx = 1.0/(PetscReal)Mx;
396: /*
397: Get pointers to vector data
398: */
399: DMDAVecGetArray(da,U,&u);
401: /*
402: Get local grid boundaries
403: */
404: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
406: /*
407: Seee heat.c for how to generate InitialSolution.heat
408: */
409: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"InitialSolution.heat",FILE_MODE_READ,&viewer);
410: VecCreate(PETSC_COMM_WORLD,&finesolution);
411: VecLoad(finesolution,viewer);
412: PetscViewerDestroy(&viewer);
413: VecGetSize(finesolution,&N);
414: scale = N/Mx;
415: VecGetArrayRead(finesolution,&f);
417: /*
418: Compute function over the locally owned part of the grid
419: */
420: for (i=xs; i<xs+xm; i++) {
421: x = i*hx;
422: r = PetscSqrtReal((x-.5)*(x-.5));
423: if (r < .125) u[i] = 1.0;
424: else u[i] = -.5;
426: /* With the initial condition above the method is first order in space */
427: /* this is a smooth initial condition so the method becomes second order in space */
428: /*u[i] = PetscSinScalar(2*PETSC_PI*x); */
429: u[i] = f[scale*i];
430: }
431: VecRestoreArrayRead(finesolution,&f);
432: VecDestroy(&finesolution);
434: /*
435: Restore vectors
436: */
437: DMDAVecRestoreArray(da,U,&u);
438: return(0);
439: }
441: /*
442: This routine is not parallel
443: */
444: PetscErrorCode MyMonitor(TS ts,PetscInt step,PetscReal time,Vec U,void *ptr)
445: {
446: UserCtx *ctx = (UserCtx*)ptr;
447: PetscDrawLG lg;
449: PetscScalar *u,l,r,c;
450: PetscInt Mx,i,xs,xm,cnt;
451: PetscReal x,y,hx,pause,sx,len,max,xx[4],yy[4],xx_netforce,yy_netforce,yup,ydown,y2,len2;
452: PetscDraw draw;
453: Vec localU;
454: DM da;
455: int colors[] = {PETSC_DRAW_YELLOW,PETSC_DRAW_RED,PETSC_DRAW_BLUE,PETSC_DRAW_PLUM,PETSC_DRAW_BLACK};
456: /*
457: const char *const legend[3][3] = {{"-kappa (\\grad u,\\grad u)","(1 - u^2)^2"},{"-kappa (\\grad u,\\grad u)","(1 - u^2)"},{"-kappa (\\grad u,\\grad u)","logarithmic"}};
458: */
459: PetscDrawAxis axis;
460: PetscDrawViewPorts *ports;
461: PetscReal tol = ctx->tol, theta=ctx->theta,theta_c=ctx->theta_c,a,b; /* a and b are used in the cubic truncation of the log function */
462: PetscReal vbounds[] = {-1.1,1.1};
465: PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,vbounds);
466: PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),800,600);
467: TSGetDM(ts,&da);
468: DMGetLocalVector(da,&localU);
469: DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
470: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
471: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
472: hx = 1.0/(PetscReal)Mx; sx = 1.0/(hx*hx);
473: DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
474: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
475: DMDAVecGetArrayRead(da,localU,&u);
477: PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,&lg);
478: PetscDrawLGGetDraw(lg,&draw);
479: PetscDrawCheckResizedWindow(draw);
480: if (!ctx->ports) {
481: PetscDrawViewPortsCreateRect(draw,1,3,&ctx->ports);
482: }
483: ports = ctx->ports;
484: PetscDrawLGGetAxis(lg,&axis);
485: PetscDrawLGReset(lg);
487: xx[0] = 0.0; xx[1] = 1.0; cnt = 2;
488: PetscOptionsGetRealArray(NULL,NULL,"-zoom",xx,&cnt,NULL);
489: xs = xx[0]/hx; xm = (xx[1] - xx[0])/hx;
491: /*
492: Plot the energies
493: */
494: PetscDrawLGSetDimension(lg,1 + (ctx->cahnhillard ? 1 : 0) + (ctx->energy == 3));
495: PetscDrawLGSetColors(lg,colors+1);
496: PetscDrawViewPortsSet(ports,2);
497: x = hx*xs;
498: for (i=xs; i<xs+xm; i++) {
499: xx[0] = xx[1] = xx[2] = x;
500: if (ctx->degenerate) yy[0] = PetscRealPart(.25*(1. - u[i]*u[i])*ctx->kappa*(u[i-1] - u[i+1])*(u[i-1] - u[i+1])*sx);
501: else yy[0] = PetscRealPart(.25*ctx->kappa*(u[i-1] - u[i+1])*(u[i-1] - u[i+1])*sx);
503: if (ctx->cahnhillard) {
504: switch (ctx->energy) {
505: case 1: /* double well */
506: yy[1] = .25*PetscRealPart((1. - u[i]*u[i])*(1. - u[i]*u[i]));
507: break;
508: case 2: /* double obstacle */
509: yy[1] = .5*PetscRealPart(1. - u[i]*u[i]);
510: break;
511: case 3: /* logarithm + double well */
512: yy[1] = .25*PetscRealPart((1. - u[i]*u[i])*(1. - u[i]*u[i]));
513: if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) yy[2] = .5*theta*(2.0*tol*PetscLogReal(tol) + PetscRealPart(1.0-u[i])*PetscLogReal(PetscRealPart(1.-u[i])/2.0));
514: else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) yy[2] = .5*theta*(PetscRealPart(1.0+u[i])*PetscLogReal(PetscRealPart(1.0+u[i])/2.0) + 2.0*tol*PetscLogReal(tol));
515: else yy[2] = .5*theta*(PetscRealPart(1.0+u[i])*PetscLogReal(PetscRealPart(1.0+u[i])/2.0) + PetscRealPart(1.0-u[i])*PetscLogReal(PetscRealPart(1.0-u[i])/2.0));
516: break;
517: case 4: /* logarithm + double obstacle */
518: yy[1] = .5*theta_c*PetscRealPart(1.0-u[i]*u[i]);
519: if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) yy[2] = .5*theta*(2.0*tol*PetscLogReal(tol) + PetscRealPart(1.0-u[i])*PetscLogReal(PetscRealPart(1.-u[i])/2.0));
520: else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) yy[2] = .5*theta*(PetscRealPart(1.0+u[i])*PetscLogReal(PetscRealPart(1.0+u[i])/2.0) + 2.0*tol*PetscLogReal(tol));
521: else yy[2] = .5*theta*(PetscRealPart(1.0+u[i])*PetscLogReal(PetscRealPart(1.0+u[i])/2.0) + PetscRealPart(1.0-u[i])*PetscLogReal(PetscRealPart(1.0-u[i])/2.0));
522: break;
523: default:
524: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"It will always be one of the values");
525: }
526: }
527: PetscDrawLGAddPoint(lg,xx,yy);
528: x += hx;
529: }
530: PetscDrawGetPause(draw,&pause);
531: PetscDrawSetPause(draw,0.0);
532: PetscDrawAxisSetLabels(axis,"Energy","","");
533: /* PetscDrawLGSetLegend(lg,legend[ctx->energy-1]); */
534: PetscDrawLGDraw(lg);
536: /*
537: Plot the forces
538: */
539: PetscDrawLGSetDimension(lg,0 + (ctx->cahnhillard ? 2 : 0) + (ctx->energy == 3));
540: PetscDrawLGSetColors(lg,colors+1);
541: PetscDrawViewPortsSet(ports,1);
542: PetscDrawLGReset(lg);
543: x = xs*hx;
544: max = 0.;
545: for (i=xs; i<xs+xm; i++) {
546: xx[0] = xx[1] = xx[2] = xx[3] = x;
547: xx_netforce = x;
548: if (ctx->degenerate) {
549: c = (1. - u[i]*u[i])*(u[i-1] + u[i+1] - 2.0*u[i])*sx;
550: r = (1. - u[i+1]*u[i+1])*(u[i] + u[i+2] - 2.0*u[i+1])*sx;
551: l = (1. - u[i-1]*u[i-1])*(u[i-2] + u[i] - 2.0*u[i-1])*sx;
552: } else {
553: c = (u[i-1] + u[i+1] - 2.0*u[i])*sx;
554: r = (u[i] + u[i+2] - 2.0*u[i+1])*sx;
555: l = (u[i-2] + u[i] - 2.0*u[i-1])*sx;
556: }
557: yy[0] = PetscRealPart(-ctx->kappa*(l + r - 2.0*c)*sx);
558: yy_netforce = yy[0];
559: max = PetscMax(max,PetscAbs(yy[0]));
560: if (ctx->cahnhillard) {
561: switch (ctx->energy) {
562: case 1: /* double well */
563: yy[1] = PetscRealPart(6.*.25*u[i]*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (3.*u[i]*u[i] - 1.)*(u[i-1] + u[i+1] - 2.0*u[i])*sx);
564: break;
565: case 2: /* double obstacle */
566: yy[1] = -PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx;
567: break;
568: case 3: /* logarithmic + double well */
569: yy[1] = PetscRealPart(6.*.25*u[i]*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (3.*u[i]*u[i] - 1.)*(u[i-1] + u[i+1] - 2.0*u[i])*sx);
570: if (ctx->truncation==2) { /* quadratic */
571: if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) yy[2] = (.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx;
572: else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) yy[2] = (.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx;
573: else yy[2] = PetscRealPart(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx);
574: } else { /* cubic */
575: a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol));
576: b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol);
577: if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) yy[2] = PetscRealPart(-1.0*a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (-1.0*a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx);
578: else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) yy[2] = PetscRealPart(1.0*a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + ( a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx);
579: else yy[2] = PetscRealPart(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx);
580: }
581: break;
582: case 4: /* logarithmic + double obstacle */
583: yy[1] = theta_c*PetscRealPart(-(u[i-1] + u[i+1] - 2.0*u[i]))*sx;
584: if (ctx->truncation==2) {
585: if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) yy[2] = (.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx;
586: else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) yy[2] = (.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx;
587: else yy[2] = PetscRealPart(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx);
588: }
589: else {
590: a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol));
591: b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol);
592: if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) yy[2] = PetscRealPart(-1.0*a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (-1.0*a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx);
593: else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) yy[2] = PetscRealPart(1.0*a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + ( a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx);
594: else yy[2] = PetscRealPart(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx);
595: }
596: break;
597: default:
598: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"It will always be one of the values");
599: }
600: if (ctx->energy < 3) {
601: max = PetscMax(max,PetscAbs(yy[1]));
602: yy[2] = yy[0]+yy[1];
603: yy_netforce = yy[2];
604: } else {
605: max = PetscMax(max,PetscAbs(yy[1]+yy[2]));
606: yy[3] = yy[0]+yy[1]+yy[2];
607: yy_netforce = yy[3];
608: }
609: }
610: if (ctx->netforce) {
611: PetscDrawLGAddPoint(lg,&xx_netforce,&yy_netforce);
612: } else {
613: PetscDrawLGAddPoint(lg,xx,yy);
614: }
615: x += hx;
616: /*if (max > 7200150000.0) */
617: /* printf("max very big when i = %d\n",i); */
618: }
619: PetscDrawAxisSetLabels(axis,"Right hand side","","");
620: PetscDrawLGSetLegend(lg,NULL);
621: PetscDrawLGDraw(lg);
623: /*
624: Plot the solution
625: */
626: PetscDrawLGSetDimension(lg,1);
627: PetscDrawViewPortsSet(ports,0);
628: PetscDrawLGReset(lg);
629: x = hx*xs;
630: PetscDrawLGSetLimits(lg,x,x+(xm-1)*hx,-1.1,1.1);
631: PetscDrawLGSetColors(lg,colors);
632: for (i=xs; i<xs+xm; i++) {
633: xx[0] = x;
634: yy[0] = PetscRealPart(u[i]);
635: PetscDrawLGAddPoint(lg,xx,yy);
636: x += hx;
637: }
638: PetscDrawAxisSetLabels(axis,"Solution","","");
639: PetscDrawLGDraw(lg);
641: /*
642: Print the forces as arrows on the solution
643: */
644: x = hx*xs;
645: cnt = xm/60;
646: cnt = (!cnt) ? 1 : cnt;
648: for (i=xs; i<xs+xm; i += cnt) {
649: y = yup = ydown = PetscRealPart(u[i]);
650: c = (u[i-1] + u[i+1] - 2.0*u[i])*sx;
651: r = (u[i] + u[i+2] - 2.0*u[i+1])*sx;
652: l = (u[i-2] + u[i] - 2.0*u[i-1])*sx;
653: len = -.5*PetscRealPart(ctx->kappa*(l + r - 2.0*c)*sx)/max;
654: PetscDrawArrow(draw,x,y,x,y+len,PETSC_DRAW_RED);
655: if (ctx->cahnhillard) {
656: if (len < 0.) ydown += len;
657: else yup += len;
659: switch (ctx->energy) {
660: case 1: /* double well */
661: len = .5*PetscRealPart(6.*.25*u[i]*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (3.*u[i]*u[i] - 1.)*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max;
662: break;
663: case 2: /* double obstacle */
664: len = -.5*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx/max;
665: break;
666: case 3: /* logarithmic + double well */
667: len = .5*PetscRealPart(6.*.25*u[i]*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (3.*u[i]*u[i] - 1.)*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max;
668: if (len < 0.) ydown += len;
669: else yup += len;
671: if (ctx->truncation==2) { /* quadratic */
672: if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) len2 = .5*(.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx/max;
673: else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) len2 = .5*(.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx/max;
674: else len2 = PetscRealPart(.5*(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max);
675: } else { /* cubic */
676: a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol));
677: b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol);
678: if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) len2 = PetscRealPart(.5*(-1.0*a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (-1.0*a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max);
679: else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) len2 = PetscRealPart(.5*(a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + ( a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max);
680: else len2 = PetscRealPart(.5*(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max);
681: }
682: y2 = len < 0 ? ydown : yup;
683: PetscDrawArrow(draw,x,y2,x,y2+len2,PETSC_DRAW_PLUM);
684: break;
685: case 4: /* logarithmic + double obstacle */
686: len = -.5*theta_c*PetscRealPart(-(u[i-1] + u[i+1] - 2.0*u[i])*sx/max);
687: if (len < 0.) ydown += len;
688: else yup += len;
690: if (ctx->truncation==2) { /* quadratic */
691: if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) len2 = .5*(.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx/max;
692: else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) len2 = .5*(.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx/max;
693: else len2 = PetscRealPart(.5*(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max);
694: } else { /* cubic */
695: a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol));
696: b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol);
697: if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) len2 = .5*PetscRealPart(-1.0*a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (-1.0*a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max;
698: else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) len2 = .5*PetscRealPart(a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + ( a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max;
699: else len2 = .5*PetscRealPart(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max;
700: }
701: y2 = len < 0 ? ydown : yup;
702: PetscDrawArrow(draw,x,y2,x,y2+len2,PETSC_DRAW_PLUM);
703: break;
704: }
705: PetscDrawArrow(draw,x,y,x,y+len,PETSC_DRAW_BLUE);
706: }
707: x += cnt*hx;
708: }
709: DMDAVecRestoreArrayRead(da,localU,&x);
710: DMRestoreLocalVector(da,&localU);
711: PetscDrawStringSetSize(draw,.2,.2);
712: PetscDrawFlush(draw);
713: PetscDrawSetPause(draw,pause);
714: PetscDrawPause(draw);
715: return(0);
716: }
718: PetscErrorCode MyDestroy(void **ptr)
719: {
720: UserCtx *ctx = *(UserCtx**)ptr;
724: PetscDrawViewPortsDestroy(ctx->ports);
725: return(0);
726: }
728: /*TEST
730: test:
731: TODO: currently requires initial condition file generated by heat
733: TEST*/