Actual source code: ex29.c
1: /*$Id: $*/
3: /*#define HAVE_DA_HDF*/
5: /* solve the equations for the perturbed magnetic field only */
6: #define EQ
8: /* turning this on causes instability?!? */
9: /* #define UPWINDING */
11: static char help[] = "Hall MHD with in two dimensions with time stepping and multigrid.\n\n\
12: -options_file ex29.options\n\
13: other PETSc options\n\
14: -resistivity <eta>\n\
15: -viscosity <nu>\n\
16: -skin_depth <d_e>\n\
17: -larmor_radius <rho_s>\n\
18: -contours : draw contour plots of solution\n\n";
20: /*T
21: Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
22: Concepts: DA^using distributed arrays;
23: Concepts: multicomponent
24: Processors: n
25: T*/
27: /* ------------------------------------------------------------------------
29: We thank Kai Germaschewski for contributing the FormFunctionLocal()
31: ------------------------------------------------------------------------- */
33: /*
34: Include "petscda.h" so that we can use distributed arrays (DAs).
35: Include "petscsnes.h" so that we can use SNES solvers.
36: Include "petscmg.h" to control the multigrid solvers.
37: Note that these automatically include:
38: petsc.h - base PETSc routines petscvec.h - vectors
39: petscsys.h - system routines petscmat.h - matrices
40: petscis.h - index sets petscksp.h - Krylov subspace methods
41: petscviewer.h - viewers petscpc.h - preconditioners
42: petscsles.h - linear solvers
43: */
44: #include petscsnes.h
45: #include petscda.h
46: #include petscmg.h
48: #ifdef HAVE_DA_HDF
49: int DAVecHDFOutput(DA da, Vec X, char *fname);
50: #endif
52: #define D_x(x,m,i,j) (p5 * (x[(j)][(i)+1].m - x[(j)][(i)-1].m) * dhx)
53: #define D_xm(x,m,i,j) ((x[(j)][(i)].m - x[(j)][(i)-1].m) * dhx)
54: #define D_xp(x,m,i,j) ((x[(j)][(i+1)].m - x[(j)][(i)].m) * dhx)
55: #define D_y(x,m,i,j) (p5 * (x[(j)+1][(i)].m - x[(j)-1][(i)].m) * dhy)
56: #define D_ym(x,m,i,j) ((x[(j)][(i)].m - x[(j)-1][(i)].m) * dhy)
57: #define D_yp(x,m,i,j) ((x[(j)+1][(i)].m - x[(j)][(i)].m) * dhy)
58: #define D_xx(x,m,i,j) ((x[(j)][(i)+1].m - two*x[(j)][(i)].m + x[(j)][(i)-1].m) * hydhx * dhxdhy)
59: #define D_yy(x,m,i,j) ((x[(j)+1][(i)].m - two*x[(j)][(i)].m + x[(j)-1][(i)].m) * hxdhy * dhxdhy)
60: #define Lapl(x,m,i,j) (D_xx(x,m,i,j) + D_yy(x,m,i,j))
61: #define lx (2.*M_PI)
62: #define ly (4.*M_PI)
63: #define sqr(a) ((a)*(a))
65: /*
66: User-defined routines and data structures
67: */
69: typedef struct {
70: PassiveScalar fnorm_ini,dt_ini;
71: PassiveScalar fnorm,dt,dt_out;
72: PassiveScalar ptime;
73: PassiveScalar max_time;
74: PassiveScalar fnorm_ratio;
75: int ires,itstep;
76: int max_steps,print_freq;
77: PassiveScalar t;
79: PetscTruth ts_monitor; /* print information about each time step */
80: PetscReal dump_time; /* time to dump solution to a file */
81: PetscViewer socketviewer; /* socket to dump solution at each timestep for visualization */
82: } TstepCtx;
84: typedef struct {
85: PetscScalar phi,psi,U,F;
86: } Field;
88: typedef struct {
89: PassiveScalar phi,psi,U,F;
90: } PassiveField;
92: typedef struct {
93: int mglevels;
94: int cycles; /* number of time steps for integration */
95: PassiveReal nu,eta,d_e,rho_s; /* physical parameters */
96: PetscTruth draw_contours; /* flag - 1 indicates drawing contours */
97: PetscTruth second_order;
98: PetscTruth PreLoading;
99: } Parameter;
101: typedef struct {
102: Vec Xoldold,Xold,func;
103: TstepCtx *tsCtx;
104: Parameter *param;
105: } AppCtx;
107: extern int FormFunctionLocal(DALocalInfo*,Field**,Field**,void*);
108: extern int Update(DMMG *);
109: extern int Initialize(DMMG *);
110: extern int AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,AppCtx *user);
111: extern int AddTSTermLocal2(DALocalInfo* info,Field **x,Field **f,AppCtx *user);
112: extern int Gnuplot(DA da, Vec X, double time);
113: extern int AttachNullSpace(PC,Vec);
114: extern int FormFunctionLocali(DALocalInfo*,MatStencil*,Field**,PetscScalar*,void*);
118: int main(int argc,char **argv)
119: {
120: DMMG *dmmg; /* multilevel grid structure */
121: AppCtx *user; /* user-defined work context (one for each level) */
122: TstepCtx tsCtx; /* time-step parameters (one total) */
123: Parameter param; /* physical parameters (one total) */
124: int i,ierr,m,n,mx,my;
125: MPI_Comm comm;
126: DA da;
127: PetscTruth flg;
128: PetscReal dt_ratio;
130: int dfill[16] = {1,0,1,0,
131: 0,1,0,1,
132: 0,0,1,1,
133: 0,1,1,1};
134: int ofill[16] = {1,0,0,0,
135: 0,1,0,0,
136: 1,1,1,1,
137: 1,1,1,1};
139: PetscInitialize(&argc,&argv,(char *)0,help);
140: comm = PETSC_COMM_WORLD;
143: PreLoadBegin(PETSC_TRUE,"SetUp");
145: param.PreLoading = PreLoading;
146: DMMGCreate(comm,1,&user,&dmmg);
147: param.mglevels = DMMGGetLevels(dmmg);
149: /*
150: Create distributed array multigrid object (DMMG) to manage
151: parallel grid and vectors for principal unknowns (x) and
152: governing residuals (f)
153: */
154: DACreate2d(comm, DA_XYPERIODIC, DA_STENCIL_STAR, -5, -5,
155: PETSC_DECIDE, PETSC_DECIDE, 4, 1, 0, 0, &da);
156:
158: /* overwrite the default sparsity pattern toone specific for
159: this code's nonzero structure */
160: DASetBlockFills(da,dfill,ofill);
162: DMMGSetDM(dmmg,(DM)da);
163: DADestroy(da);
165: /* default physical parameters */
166: param.nu = 0;
167: param.eta = 1e-3;
168: param.d_e = 0.2;
169: param.rho_s = 0;
171: PetscOptionsGetReal(PETSC_NULL, "-viscosity", ¶m.nu,
172: PETSC_NULL);
173:
175: PetscOptionsGetReal(PETSC_NULL, "-resistivity", ¶m.eta,
176: PETSC_NULL);
177:
179: PetscOptionsGetReal(PETSC_NULL, "-skin_depth", ¶m.d_e,
180: PETSC_NULL);
181:
183: PetscOptionsGetReal(PETSC_NULL, "-larmor_radius", ¶m.rho_s,
184: PETSC_NULL);
185:
187: PetscOptionsHasName(PETSC_NULL, "-contours", ¶m.draw_contours);
188:
190: PetscOptionsHasName(PETSC_NULL, "-second_order", ¶m.second_order);
191:
193: DASetFieldName(DMMGGetDA(dmmg), 0, "phi");
194:
196: DASetFieldName(DMMGGetDA(dmmg), 1, "psi");
197:
199: DASetFieldName(DMMGGetDA(dmmg), 2, "U");
200:
202: DASetFieldName(DMMGGetDA(dmmg), 3, "F");
203:
205: /*======================================================================*/
206: /* Initialize stuff related to time stepping */
207: /*======================================================================*/
209: DAGetInfo(DMMGGetDA(dmmg),0,&mx,&my,0,0,0,0,0,0,0,0);
211: tsCtx.fnorm_ini = 0;
212: tsCtx.max_steps = 1000000;
213: tsCtx.max_time = 1.0e+12;
214: /* use for dt = min(dx,dy); multiplied by dt_ratio below */
215: tsCtx.dt = PetscMin(lx/mx,ly/my);
216: tsCtx.fnorm_ratio = 1.0e+10;
217: tsCtx.t = 0;
218: tsCtx.dt_out = 10;
219: tsCtx.print_freq = tsCtx.max_steps;
220: tsCtx.ts_monitor = PETSC_FALSE;
221: tsCtx.dump_time = -1.0;
223: PetscOptionsGetInt(PETSC_NULL, "-max_st", &tsCtx.max_steps,
224: PETSC_NULL);
225:
227: PetscOptionsGetReal(PETSC_NULL, "-ts_rtol", &tsCtx.fnorm_ratio,
228: PETSC_NULL);
229:
231: PetscOptionsGetInt(PETSC_NULL, "-print_freq", &tsCtx.print_freq,
232: PETSC_NULL);
233:
235: dt_ratio = 1.0;
236: PetscOptionsGetReal(PETSC_NULL, "-dt_ratio", &dt_ratio,
237: PETSC_NULL);
238:
239: tsCtx.dt *= dt_ratio;
242: PetscOptionsHasName(PETSC_NULL, "-ts_monitor", &tsCtx.ts_monitor);
243:
245: PetscOptionsGetReal(PETSC_NULL, "-dump", &tsCtx.dump_time,
246: PETSC_NULL);
247:
249: tsCtx.socketviewer = 0;
250: PetscOptionsHasName(PETSC_NULL, "-socket_viewer", &flg);
251:
252: if (flg && !PreLoading) {
253: tsCtx.ts_monitor = PETSC_TRUE;
254: PetscViewerSocketOpen(PETSC_COMM_WORLD,0,PETSC_DECIDE,&tsCtx.socketviewer);
255: }
256:
258: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
259: Create user context, set problem data, create vector data structures.
260: Also, compute the initial guess.
261: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
262: /* create application context for each level */
263: PetscMalloc(param.mglevels*sizeof(AppCtx),&user);
264: for (i=0; i<param.mglevels; i++) {
265: /* create work vectors to hold previous time-step solution and
266: function value */
267: VecDuplicate(dmmg[i]->x, &user[i].Xoldold);
268: VecDuplicate(dmmg[i]->x, &user[i].Xold);
269: VecDuplicate(dmmg[i]->x, &user[i].func);
270: user[i].tsCtx = &tsCtx;
271: user[i].param = ¶m;
272: dmmg[i]->user = &user[i];
273: }
274:
275: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
276: Create nonlinear solver context
277:
278: Process adiC(20): AddTSTermLocal AddTSTermLocal2 FormFunctionLocal FormFunctionLocali
279: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
280: DMMGSetSNESLocal(dmmg, FormFunctionLocal, 0,
281: ad_FormFunctionLocal, admf_FormFunctionLocal);
282:
283: DMMGSetSNESLocali(dmmg,FormFunctionLocali,ad_FormFunctionLocali,admf_FormFunctionLocali);
285: /* attach nullspace to each level of the preconditioner */
286: {
287: SLES subsles,sles;
288: PC pc,subpc;
289: PetscTruth mg;
291: SNESGetSLES(DMMGGetSNES(dmmg),&sles);
292: SLESGetPC(sles,&pc);
293: AttachNullSpace(pc,DMMGGetx(dmmg));
294: PetscTypeCompare((PetscObject)pc,PCMG,&mg);
295: if (mg) {
296: for (i=0; i<param.mglevels; i++) {
297: MGGetSmoother(pc,i,&subsles);
298: SLESGetPC(subsles,&subpc);
299: AttachNullSpace(subpc,dmmg[i]->x);
300: }
301: }
302: }
303:
304: if (PreLoading) {
305: PetscPrintf(comm, "# viscosity = %g, resistivity = %g, "
306: "skin_depth # = %g, larmor_radius # = %g\n",
307: param.nu, param.eta, param.d_e, param.rho_s);
308:
309: DAGetInfo(DMMGGetDA(dmmg),0,&m,&n,0,0,0,0,0,0,0,0);
310: PetscPrintf(comm,"Problem size %d by %d\n",m,n);
311: PetscPrintf(comm,"dx %g dy %g dt %g ratio dt/min(dx,dy) %g\n",lx/mx,ly/my,tsCtx.dt,dt_ratio);
312: }
314:
315:
316: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
317: Solve the nonlinear system
318: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
319:
320: PreLoadStage("Solve");
322: if (param.draw_contours) {
323: VecView(((AppCtx*)dmmg[param.mglevels-1]->user)->Xold,
324: PETSC_VIEWER_DRAW_WORLD);
325:
326: }
328: Update(dmmg);
329:
331: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
332: Free work space. All PETSc objects should be destroyed when they
333: are no longer needed.
334: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
335:
336: for (i=0; i<param.mglevels; i++) {
337: VecDestroy(user[i].Xoldold);
338: VecDestroy(user[i].Xold);
339: VecDestroy(user[i].func);
340: }
341: PetscFree(user);
342: DMMGDestroy(dmmg);
344: PreLoadEnd();
345:
346: PetscFinalize();
347:
349: return 0;
350: }
352: /* ------------------------------------------------------------------- */
355: /* ------------------------------------------------------------------- */
356: int Gnuplot(DA da, Vec X, double time)
357: {
358: int i,j,xs,ys,xm,ym;
359: int xints,xinte,yints,yinte;
360: int ierr;
361: Field **x;
362: FILE *f;
363: char fname[100];
364: int cpu;
366: MPI_Comm_rank(PETSC_COMM_WORLD,&cpu);
367:
369: sprintf(fname, "out-%g-%d.dat", time, cpu);
371: f = fopen(fname, "w");
373: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
374:
376: DAVecGetArray(da,X,&x);
377:
379: xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;
381: for (j=yints; j<yinte; j++) {
382: for (i=xints; i<xinte; i++) {
383: PetscFPrintf(PETSC_COMM_WORLD, f,
384: "%d %d %g %g %g %g %g %g\n",
385: i, j, 0.0, 0.0,
386: x[j][i].U, x[j][i].F, x[j][i].phi, x[j][i].psi);
387:
388: }
389: PetscFPrintf(PETSC_COMM_WORLD,f, "\n");
390:
391: }
393: DAVecRestoreArray(da,X,&x);
394:
396: fclose(f);
398: return 0;
399: }
401: /* ------------------------------------------------------------------- */
404: /* ------------------------------------------------------------------- */
405: int Initialize(DMMG *dmmg)
406: {
407: AppCtx *appCtx = (AppCtx*)dmmg[0]->user;
408: Parameter *param = appCtx->param;
409: DA da;
410: int i,j,mx,my,ierr,xs,ys,xm,ym;
411: PetscReal two = 2.0,one = 1.0;
412: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
413: PetscReal d_e,rho_s,de2,xx,yy;
414: Field **x, **localx;
415: Vec localX;
416: PetscTruth flg;
419: PetscOptionsHasName(0,"-restart",&flg);
420: if (flg) {
421: PetscViewer viewer;
422: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"binaryoutput",PETSC_BINARY_RDONLY,&viewer);
423: VecLoadIntoVector(viewer,dmmg[param->mglevels-1]->x);
424: PetscViewerDestroy(viewer);
425: return(0);
426: }
428: d_e = param->d_e;
429: rho_s = param->rho_s;
430: de2 = sqr(param->d_e);
432: da = (DA)(dmmg[param->mglevels-1]->dm);
433: DAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0);
435: dhx = mx/lx; dhy = my/ly;
436: hx = one/dhx; hy = one/dhy;
437: hxdhy = hx*dhy; hydhx = hy*dhx;
438: hxhy = hx*hy; dhxdhy = dhx*dhy;
440: /*
441: Get local grid boundaries (for 2-dimensional DA):
442: xs, ys - starting grid indices (no ghost points)
443: xm, ym - widths of local grid (no ghost points)
444: */
445: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
447: DAGetLocalVector(da,&localX);
448: /*
449: Get a pointer to vector data.
450: - For default PETSc vectors, VecGetArray() returns a pointer to
451: the data array. Otherwise, the routine is implementation dependent.
452: - You MUST call VecRestoreArray() when you no longer need access to
453: the array.
454: */
455: DAVecGetArray(da,dmmg[param->mglevels-1]->x,&x);
456: DAVecGetArray(da,localX,&localx);
458: /*
459: Compute initial solution over the locally owned part of the grid
460: */
461: {
462: PetscReal eps = lx/ly;
463: PetscReal pert = 1e-4;
464: PetscReal k = 1.*eps;
465: PetscReal gam;
467: if (d_e < rho_s) d_e = rho_s;
468: gam = k * d_e;
470: for (j=ys-1; j<ys+ym+1; j++) {
471: yy = j * hy;
472: for (i=xs-1; i<xs+xm+1; i++) {
473: xx = i * hx;
475: if (xx < -M_PI/2) {
476: localx[j][i].phi = pert * gam / k * erf((xx + M_PI) / (sqrt(2) * d_e)) * (-sin(k*yy));
477: } else if (xx < M_PI/2) {
478: localx[j][i].phi = - pert * gam / k * erf(xx / (sqrt(2) * d_e)) * (-sin(k*yy));
479: } else if (xx < 3*M_PI/2){
480: localx[j][i].phi = pert * gam / k * erf((xx - M_PI) / (sqrt(2) * d_e)) * (-sin(k*yy));
481: } else {
482: localx[j][i].phi = - pert * gam / k * erf((xx - 2.*M_PI) / (sqrt(2) * d_e)) * (-sin(k*yy));
483: }
484: #ifdef EQ
485: localx[j][i].psi = 0.;
486: #else
487: localx[j][i].psi = cos(xx);
488: #endif
489: }
490: }
491: for (j=ys; j<ys+ym; j++) {
492: for (i=xs; i<xs+xm; i++) {
493: x[j][i].U = Lapl(localx,phi,i,j);
494: x[j][i].F = localx[j][i].psi - de2 * Lapl(localx,psi,i,j);
495: x[j][i].phi = localx[j][i].phi;
496: x[j][i].psi = localx[j][i].psi;
497: }
498: }
499: }
500: /*
501: Restore vector
502: */
503: DAVecRestoreArray(da,dmmg[param->mglevels-1]->x,&x);
504:
505: DAVecRestoreArray(da,localX,&localx);
506:
507: DARestoreLocalVector(da,&localX);
508:
510: return(0);
511: }
515: int ComputeMaxima(DA da, Vec X, PetscReal t)
516: {
517: int ierr,i,j,mx,my,xs,ys,xm,ym;
518: int xints,xinte,yints,yinte;
519: Field **x;
520: double norm[4] = {0,0,0,0};
521: double gnorm[4];
522: MPI_Comm comm;
526: DAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0);
527:
529: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
530:
532: xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;
534: DAVecGetArray(da, X, &x);
535:
537: for (j=yints; j<yinte; j++) {
538: for (i=xints; i<xinte; i++) {
539: norm[0] = PetscMax(norm[0],x[j][i].U);
540: norm[1] = PetscMax(norm[1],x[j][i].F);
541: norm[2] = PetscMax(norm[2],x[j][i].phi);
542: norm[3] = PetscMax(norm[3],x[j][i].psi);
543: }
544: }
546: DAVecRestoreArray(da,X,&x);
547:
549: PetscObjectGetComm((PetscObject)da, &comm);
550:
551: MPI_Allreduce(norm, gnorm, 4, MPI_DOUBLE, MPI_MAX, comm);
552:
554: PetscFPrintf(PETSC_COMM_WORLD, stderr,
555: "%g\t%g\t%g\t%g\t%g\n",
556: t, gnorm[0], gnorm[1], gnorm[2], gnorm[3]);
557:
559: return(0);
560: }
564: int FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
565: {
566: AppCtx *user = (AppCtx*)ptr;
567: TstepCtx *tsCtx = user->tsCtx;
568: Parameter *param = user->param;
569: int ierr,i,j;
570: int xints,xinte,yints,yinte;
571: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
572: PassiveReal de2,rhos2,nu,eta,dde2;
573: PassiveReal two = 2.0,one = 1.0,p5 = 0.5;
574: PassiveReal F_eq_x,By_eq;
576: PetscScalar xx;
577: PetscScalar vx,vy,avx,avy,vxp,vxm,vyp,vym;
578: PetscScalar Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
581: de2 = sqr(user->param->d_e);
582: rhos2 = sqr(user->param->rho_s);
583: nu = user->param->nu;
584: eta = user->param->eta;
585: dde2 = one/de2;
587: /*
588: Define mesh intervals ratios for uniform grid.
589: [Note: FD formulae below are normalized by multiplying through by
590: local volume element to obtain coefficients O(1) in two dimensions.]
591: */
592: dhx = info->mx/lx; dhy = info->my/ly;
593: hx = one/dhx; hy = one/dhy;
594: hxdhy = hx*dhy; hydhx = hy*dhx;
595: hxhy = hx*hy; dhxdhy = dhx*dhy;
597: xints = info->xs; xinte = info->xs+info->xm;
598: yints = info->ys; yinte = info->ys+info->ym;
600: /* Compute over the interior points */
601: for (j=yints; j<yinte; j++) {
602: for (i=xints; i<xinte; i++) {
603: #ifdef EQ
604: xx = i * hx;
605: F_eq_x = - (1. + de2) * sin(xx);
606: By_eq = sin(xx);
607: #else
608: F_eq_x = 0.;
609: By_eq = 0.;
610: #endif
612: /*
613: * convective coefficients for upwinding
614: */
616: vx = - D_y(x,phi,i,j);
617: vy = D_x(x,phi,i,j);
618: avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
619: avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
620: #ifndef UPWINDING
621: vxp = vxm = p5*vx;
622: vyp = vym = p5*vy;
623: #endif
625: Bx = D_y(x,psi,i,j);
626: By = - D_x(x,psi,i,j) + By_eq;
627: aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
628: aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
629: #ifndef UPWINDING
630: Bxp = Bxm = p5*Bx;
631: Byp = Bym = p5*By;
632: #endif
634: /* Lap(phi) - U */
635: f[j][i].phi = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
637: /* psi - d_e^2 * Lap(psi) - F */
638: f[j][i].psi = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
640: /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
641: - nu Lap(U) */
642: f[j][i].U =
643: ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
644: vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
645: (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
646: Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
647: nu * Lapl(x,U,i,j)) * hxhy;
648:
649: /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
650: - eta * Lap(psi) */
651: f[j][i].F =
652: ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
653: vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
654: (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
655: Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
656: eta * Lapl(x,psi,i,j)) * hxhy;
657: }
658: }
660: /* Add time step contribution */
661: if (tsCtx->ires) {
662: if ((param->second_order) && (tsCtx->itstep > 0))
663: {
664: AddTSTermLocal2(info,x,f,user);
665:
666: }
667: else
668: {
669: AddTSTermLocal(info,x,f,user);
670:
671: }
672: }
674: /*
675: Flop count (multiply-adds are counted as 2 operations)
676: */
677: /* PetscLogFlops(84*info->ym*info->xm); FIXME */
678: return(0);
679: }
681: /*---------------------------------------------------------------------*/
684: int Update(DMMG *dmmg)
685: /*---------------------------------------------------------------------*/
686: {
687: AppCtx *user = (AppCtx *) ((dmmg[0])->user);
688: TstepCtx *tsCtx = user->tsCtx;
689: Parameter *param = user->param;
690: SNES snes;
691: int ierr,its,lits,i;
692: int max_steps;
693: int nfailsCum = 0,nfails = 0;
694: static int ic_out;
695: PetscTruth ts_monitor = (tsCtx->ts_monitor && !user->param->PreLoading) ? PETSC_TRUE : PETSC_FALSE;
699: if (user->param->PreLoading)
700: max_steps = 1;
701: else
702: max_steps = tsCtx->max_steps;
703:
704: Initialize(dmmg);
705:
707: for (tsCtx->itstep = 0; tsCtx->itstep < max_steps; tsCtx->itstep++) {
709: if ((param->second_order) && (tsCtx->itstep > 0))
710: {
711: for (i=param->mglevels-1; i>=0 ;i--)
712: {
713: VecCopy(((AppCtx*)dmmg[i]->user)->Xold,
714: ((AppCtx*)dmmg[i]->user)->Xoldold);
715:
716: }
717: }
719: for (i=param->mglevels-1; i>0 ;i--) {
720: MatRestrict(dmmg[i]->R, dmmg[i]->x, dmmg[i-1]->x);
721:
722: VecPointwiseMult(dmmg[i]->Rscale,dmmg[i-1]->x,dmmg[i-1]->x);
723:
724: VecCopy(dmmg[i]->x, ((AppCtx*)dmmg[i]->user)->Xold);
725:
726: }
727: VecCopy(dmmg[0]->x, ((AppCtx*)dmmg[0]->user)->Xold);
728:
730: DMMGSolve(dmmg);
731:
733: snes = DMMGGetSNES(dmmg);
736: if (tsCtx->itstep == 665000)
737: {
738: SLES sles;
739: PC pc;
740: Mat mat, pmat;
741: MatStructure flag;
742: PetscViewer viewer;
743: char file[128];
745: strcpy(file, "matrix");
747: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file,
748: PETSC_BINARY_CREATE, &viewer);
749:
751: SNESGetSLES(snes, &sles);
752:
754: SLESGetPC(sles, &pc);
755:
757: PCGetOperators(pc, &mat, &pmat, &flag);
758:
760: MatView(mat, viewer);
761:
763: PetscViewerDestroy(viewer);
764:
765: SETERRQ(1,"Done saving Jacobian");
766: }
769: tsCtx->t += tsCtx->dt;
771: /* save restart solution if requested at a particular time, then exit */
772: if (tsCtx->dump_time > 0.0 && tsCtx->t >= tsCtx->dump_time) {
773: Vec v = DMMGGetx(dmmg);
774: VecView(v,PETSC_VIEWER_BINARY_WORLD);
775: SETERRQ1(1,"Saved solution at time %g",tsCtx->t);
776: }
778: if (ts_monitor) {
779: SNESGetIterationNumber(snes,&its);
780: SNESGetNumberLinearIterations(snes,&lits);
781: SNESGetNumberUnsuccessfulSteps(snes,&nfails);
782: SNESGetFunctionNorm(snes,&tsCtx->fnorm);
783: nfailsCum += nfails;
784: if (nfailsCum >= 2) SETERRQ(1,"Unable to find a Newton Step");
785: PetscPrintf(PETSC_COMM_WORLD,"Time Step %d time = %g Newton steps %d linear steps %d fnorm %g\n",
786: tsCtx->itstep, tsCtx->t, its, lits, tsCtx->fnorm);
787:
789: /* send solution over to Matlab, to be visualized (using ex29.m) */
790: if (!user->param->PreLoading && tsCtx->socketviewer){
791: Vec v;
792: SNESGetSolution(snes,&v);
793: VecView(v,tsCtx->socketviewer);
794: }
795: }
797: if (!param->PreLoading) {
798: if (param->draw_contours) {
799: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
800:
801: }
803: if (0 && ts_monitor) {
804: /* compute maxima */
805: ComputeMaxima((DA) dmmg[param->mglevels-1]->dm,
806: dmmg[param->mglevels-1]->x,
807: tsCtx->t);
808: /* output */
809: if (ic_out++ == (int)(tsCtx->dt_out / tsCtx->dt + .5)) {
810: #ifdef HAVE_DA_HDF
811: char fname[128];
813: sprintf(fname, "out-%g.hdf", tsCtx->t);
814: DAVecHDFOutput(DMMGGetDA(dmmg), DMMGGetx(dmmg), fname);
815:
816: #else
817: /*
818: Gnuplot(DMMGGetDA(dmmg), DMMGGetx(dmmg), tsCtx->t);
819:
820: */
821: #endif
822: ic_out = 1;
823: }
824: }
825: }
826: } /* End of time step loop */
827:
828: if (!user->param->PreLoading){
829: SNESGetFunctionNorm(snes,&tsCtx->fnorm);
830: PetscPrintf(PETSC_COMM_WORLD, "timesteps %d fnorm = %g\n",
831: tsCtx->itstep, tsCtx->fnorm);
832:
833: }
835: if (user->param->PreLoading) {
836: tsCtx->fnorm_ini = 0.0;
837: }
838:
839: return(0);
840: }
842: /*---------------------------------------------------------------------*/
845: int AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,AppCtx *user)
846: /*---------------------------------------------------------------------*/
847: {
848: TstepCtx *tsCtx = user->tsCtx;
849: DA da = info->da;
850: int ierr,i,j;
851: int xints,xinte,yints,yinte;
852: PassiveReal hx,hy,dhx,dhy,hxhy;
853: PassiveReal one = 1.0;
854: PassiveScalar dtinv;
855: PassiveField **xold;
859: xints = info->xs; xinte = info->xs+info->xm;
860: yints = info->ys; yinte = info->ys+info->ym;
862: dhx = info->mx/lx; dhy = info->my/ly;
863: hx = one/dhx; hy = one/dhy;
864: hxhy = hx*hy;
866: dtinv = hxhy/(tsCtx->dt);
868: DAVecGetArray(da,user->Xold,&xold);
869: for (j=yints; j<yinte; j++) {
870: for (i=xints; i<xinte; i++) {
871: f[j][i].U += dtinv*(x[j][i].U-xold[j][i].U);
872: f[j][i].F += dtinv*(x[j][i].F-xold[j][i].F);
873: }
874: }
875: DAVecRestoreArray(da,user->Xold,&xold);
877: return(0);
878: }
880: /*---------------------------------------------------------------------*/
883: int AddTSTermLocal2(DALocalInfo* info,Field **x,Field **f,AppCtx *user)
884: /*---------------------------------------------------------------------*/
885: {
886: TstepCtx *tsCtx = user->tsCtx;
887: DA da = info->da;
888: int ierr,i,j;
889: int xints,xinte,yints,yinte;
890: PassiveReal hx,hy,dhx,dhy,hxhy;
891: PassiveReal one = 1.0, onep5 = 1.5, two = 2.0, p5 = 0.5;
892: PassiveScalar dtinv;
893: PassiveField **xoldold,**xold;
897: xints = info->xs; xinte = info->xs+info->xm;
898: yints = info->ys; yinte = info->ys+info->ym;
900: dhx = info->mx/lx; dhy = info->my/ly;
901: hx = one/dhx; hy = one/dhy;
902: hxhy = hx*hy;
904: dtinv = hxhy/(tsCtx->dt);
906: DAVecGetArray(da,user->Xoldold,&xoldold);
907: DAVecGetArray(da,user->Xold,&xold);
908: for (j=yints; j<yinte; j++) {
909: for (i=xints; i<xinte; i++) {
910: f[j][i].U += dtinv * (onep5 * x[j][i].U - two * xold[j][i].U +
911: p5 * xoldold[j][i].U);
912: f[j][i].F += dtinv * (onep5 * x[j][i].F - two * xold[j][i].F +
913: p5 * xoldold[j][i].F);
914: }
915: }
916: DAVecRestoreArray(da,user->Xoldold,&xoldold);
917: DAVecRestoreArray(da,user->Xold,&xold);
919: return(0);
920: }
922: /* -------------------------------------------------------------------------*/
925: int AttachNullSpace(PC pc,Vec model)
926: {
927: int i,ierr,rstart,rend,N;
928: MatNullSpace sp;
929: Vec v,vs[1];
930: PetscScalar *vx,scale;
933: VecDuplicate(model,&v);
934: VecGetSize(model,&N);
935: scale = 2.0/sqrt(N);
936: VecGetArray(v,&vx);
937: VecGetOwnershipRange(v,&rstart,&rend);
938: for (i=rstart; i<rend; i++) {
939: if (!(i % 4)) vx[i-rstart] = scale;
940: else vx[i-rstart] = 0.0;
941: }
942: VecRestoreArray(v,&vx);
943: vs[0] = v;
944: MatNullSpaceCreate(PETSC_COMM_WORLD,0,1,vs,&sp);
945: VecDestroy(v);
946: PCNullSpaceAttach(pc,sp);
947: MatNullSpaceDestroy(sp);
948: return(0);
949: }
951: /*
952: This is an experimental function and can be safely ignored.
953: */
954: int FormFunctionLocali(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
955: {
956: AppCtx *user = (AppCtx*)ptr;
957: TstepCtx *tsCtx = user->tsCtx;
958: int ierr,i,j,c;
959: int xints,xinte,yints,yinte;
960: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
961: PassiveReal de2,rhos2,nu,eta,dde2;
962: PassiveReal two = 2.0,one = 1.0,p5 = 0.5;
963: PassiveReal F_eq_x,By_eq,dtinv;
965: PetscScalar xx;
966: PetscScalar vx,vy,avx,avy,vxp,vxm,vyp,vym;
967: PetscScalar Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
968: PassiveField **xold;
971: de2 = sqr(user->param->d_e);
972: rhos2 = sqr(user->param->rho_s);
973: nu = user->param->nu;
974: eta = user->param->eta;
975: dde2 = one/de2;
977: /*
978: Define mesh intervals ratios for uniform grid.
979: [Note: FD formulae below are normalized by multiplying through by
980: local volume element to obtain coefficients O(1) in two dimensions.]
981: */
982: dhx = info->mx/lx; dhy = info->my/ly;
983: hx = one/dhx; hy = one/dhy;
984: hxdhy = hx*dhy; hydhx = hy*dhx;
985: hxhy = hx*hy; dhxdhy = dhx*dhy;
987: xints = info->xs; xinte = info->xs+info->xm;
988: yints = info->ys; yinte = info->ys+info->ym;
991: i = st->i; j = st->j; c = st->c;
993: #ifdef EQ
994: xx = i * hx;
995: F_eq_x = - (1. + de2) * sin(xx);
996: By_eq = sin(xx);
997: #else
998: F_eq_x = 0.;
999: By_eq = 0.;
1000: #endif
1002: /*
1003: * convective coefficients for upwinding
1004: */
1006: vx = - D_y(x,phi,i,j);
1007: vy = D_x(x,phi,i,j);
1008: avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
1009: avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
1010: #ifndef UPWINDING
1011: vxp = vxm = p5*vx;
1012: vyp = vym = p5*vy;
1013: #endif
1015: Bx = D_y(x,psi,i,j);
1016: By = - D_x(x,psi,i,j) + By_eq;
1017: aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
1018: aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
1019: #ifndef UPWINDING
1020: Bxp = Bxm = p5*Bx;
1021: Byp = Bym = p5*By;
1022: #endif
1024: DAVecGetArray(info->da,user->Xold,&xold);
1025: dtinv = hxhy/(tsCtx->dt);
1026: switch(c) {
1028: case 0:
1029: /* Lap(phi) - U */
1030: *f = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
1031: break;
1033: case 1:
1034: /* psi - d_e^2 * Lap(psi) - F */
1035: *f = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
1036: break;
1038: case 2:
1039: /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
1040: - nu Lap(U) */
1041: *f =
1042: ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
1043: vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
1044: (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
1045: Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
1046: nu * Lapl(x,U,i,j)) * hxhy;
1047: *f += dtinv*(x[j][i].U-xold[j][i].U);
1048: break;
1050: case 3:
1051: /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
1052: - eta * Lap(psi) */
1053: *f =
1054: ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
1055: vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
1056: (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
1057: Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
1058: eta * Lapl(x,psi,i,j)) * hxhy;
1059: *f += dtinv*(x[j][i].F-xold[j][i].F);
1060: break;
1061: }
1062: DAVecRestoreArray(info->da,user->Xold,&xold);
1065: /*
1066: Flop count (multiply-adds are counted as 2 operations)
1067: */
1068: /* PetscLogFlops(84*info->ym*info->xm); FIXME */
1069: return(0);
1070: }