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: petscksp.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: KSP subksp,ksp;
288: PC pc,subpc;
289: PetscTruth mg;
290: KSP ksp;
292: SNESGetKSP(DMMGGetSNES(dmmg),&ksp);
293: KSPGetPC(ksp,&pc);
294: AttachNullSpace(pc,DMMGGetx(dmmg));
295: PetscTypeCompare((PetscObject)pc,PCMG,&mg);
296: if (mg) {
297: for (i=0; i<param.mglevels; i++) {
298: MGGetSmoother(pc,i,&subksp);
299: KSPGetPC(subksp,&subpc);
300: AttachNullSpace(subpc,dmmg[i]->x);
301: }
302: }
303: }
304:
305: if (PreLoading) {
306: PetscPrintf(comm, "# viscosity = %g, resistivity = %g, "
307: "skin_depth # = %g, larmor_radius # = %g\n",
308: param.nu, param.eta, param.d_e, param.rho_s);
309:
310: DAGetInfo(DMMGGetDA(dmmg),0,&m,&n,0,0,0,0,0,0,0,0);
311: PetscPrintf(comm,"Problem size %d by %d\n",m,n);
312: PetscPrintf(comm,"dx %g dy %g dt %g ratio dt/min(dx,dy) %g\n",lx/mx,ly/my,tsCtx.dt,dt_ratio);
313: }
315:
316:
317: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
318: Solve the nonlinear system
319: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
320:
321: PreLoadStage("Solve");
323: if (param.draw_contours) {
324: VecView(((AppCtx*)dmmg[param.mglevels-1]->user)->Xold,
325: PETSC_VIEWER_DRAW_WORLD);
326:
327: }
329: Update(dmmg);
330:
332: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
333: Free work space. All PETSc objects should be destroyed when they
334: are no longer needed.
335: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
336:
337: for (i=0; i<param.mglevels; i++) {
338: VecDestroy(user[i].Xoldold);
339: VecDestroy(user[i].Xold);
340: VecDestroy(user[i].func);
341: }
342: PetscFree(user);
343: DMMGDestroy(dmmg);
345: PreLoadEnd();
346:
347: PetscFinalize();
348:
350: return 0;
351: }
353: /* ------------------------------------------------------------------- */
356: /* ------------------------------------------------------------------- */
357: int Gnuplot(DA da, Vec X, double time)
358: {
359: int i,j,xs,ys,xm,ym;
360: int xints,xinte,yints,yinte;
361: int ierr;
362: Field **x;
363: FILE *f;
364: char fname[100];
365: int cpu;
367: MPI_Comm_rank(PETSC_COMM_WORLD,&cpu);
368:
370: sprintf(fname, "out-%g-%d.dat", time, cpu);
372: f = fopen(fname, "w");
374: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
375:
377: DAVecGetArray(da,X,&x);
378:
380: xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;
382: for (j=yints; j<yinte; j++) {
383: for (i=xints; i<xinte; i++) {
384: PetscFPrintf(PETSC_COMM_WORLD, f,
385: "%d %d %g %g %g %g %g %g\n",
386: i, j, 0.0, 0.0,
387: x[j][i].U, x[j][i].F, x[j][i].phi, x[j][i].psi);
388:
389: }
390: PetscFPrintf(PETSC_COMM_WORLD,f, "\n");
391:
392: }
394: DAVecRestoreArray(da,X,&x);
395:
397: fclose(f);
399: return 0;
400: }
402: /* ------------------------------------------------------------------- */
405: /* ------------------------------------------------------------------- */
406: int Initialize(DMMG *dmmg)
407: {
408: AppCtx *appCtx = (AppCtx*)dmmg[0]->user;
409: Parameter *param = appCtx->param;
410: DA da;
411: int i,j,mx,my,ierr,xs,ys,xm,ym;
412: PetscReal two = 2.0,one = 1.0;
413: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
414: PetscReal d_e,rho_s,de2,xx,yy;
415: Field **x, **localx;
416: Vec localX;
417: PetscTruth flg;
420: PetscOptionsHasName(0,"-restart",&flg);
421: if (flg) {
422: PetscViewer viewer;
423: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"binaryoutput",PETSC_FILE_RDONLY,&viewer);
424: VecLoadIntoVector(viewer,dmmg[param->mglevels-1]->x);
425: PetscViewerDestroy(viewer);
426: return(0);
427: }
429: d_e = param->d_e;
430: rho_s = param->rho_s;
431: de2 = sqr(param->d_e);
433: da = (DA)(dmmg[param->mglevels-1]->dm);
434: DAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0);
436: dhx = mx/lx; dhy = my/ly;
437: hx = one/dhx; hy = one/dhy;
438: hxdhy = hx*dhy; hydhx = hy*dhx;
439: hxhy = hx*hy; dhxdhy = dhx*dhy;
441: /*
442: Get local grid boundaries (for 2-dimensional DA):
443: xs, ys - starting grid indices (no ghost points)
444: xm, ym - widths of local grid (no ghost points)
445: */
446: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
448: DAGetLocalVector(da,&localX);
449: /*
450: Get a pointer to vector data.
451: - For default PETSc vectors, VecGetArray() returns a pointer to
452: the data array. Otherwise, the routine is implementation dependent.
453: - You MUST call VecRestoreArray() when you no longer need access to
454: the array.
455: */
456: DAVecGetArray(da,dmmg[param->mglevels-1]->x,&x);
457: DAVecGetArray(da,localX,&localx);
459: /*
460: Compute initial solution over the locally owned part of the grid
461: */
462: {
463: PetscReal eps = lx/ly;
464: PetscReal pert = 1e-4;
465: PetscReal k = 1.*eps;
466: PetscReal gam;
468: if (d_e < rho_s) d_e = rho_s;
469: gam = k * d_e;
471: for (j=ys-1; j<ys+ym+1; j++) {
472: yy = j * hy;
473: for (i=xs-1; i<xs+xm+1; i++) {
474: xx = i * hx;
476: if (xx < -M_PI/2) {
477: localx[j][i].phi = pert * gam / k * erf((xx + M_PI) / (sqrt(2) * d_e)) * (-sin(k*yy));
478: } else if (xx < M_PI/2) {
479: localx[j][i].phi = - pert * gam / k * erf(xx / (sqrt(2) * d_e)) * (-sin(k*yy));
480: } else if (xx < 3*M_PI/2){
481: localx[j][i].phi = pert * gam / k * erf((xx - M_PI) / (sqrt(2) * d_e)) * (-sin(k*yy));
482: } else {
483: localx[j][i].phi = - pert * gam / k * erf((xx - 2.*M_PI) / (sqrt(2) * d_e)) * (-sin(k*yy));
484: }
485: #ifdef EQ
486: localx[j][i].psi = 0.;
487: #else
488: localx[j][i].psi = cos(xx);
489: #endif
490: }
491: }
492: for (j=ys; j<ys+ym; j++) {
493: for (i=xs; i<xs+xm; i++) {
494: x[j][i].U = Lapl(localx,phi,i,j);
495: x[j][i].F = localx[j][i].psi - de2 * Lapl(localx,psi,i,j);
496: x[j][i].phi = localx[j][i].phi;
497: x[j][i].psi = localx[j][i].psi;
498: }
499: }
500: }
501: /*
502: Restore vector
503: */
504: DAVecRestoreArray(da,dmmg[param->mglevels-1]->x,&x);
505:
506: DAVecRestoreArray(da,localX,&localx);
507:
508: DARestoreLocalVector(da,&localX);
509:
511: return(0);
512: }
516: int ComputeMaxima(DA da, Vec X, PetscReal t)
517: {
518: int ierr,i,j,mx,my,xs,ys,xm,ym;
519: int xints,xinte,yints,yinte;
520: Field **x;
521: double norm[4] = {0,0,0,0};
522: double gnorm[4];
523: MPI_Comm comm;
527: DAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0);
528:
530: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
531:
533: xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;
535: DAVecGetArray(da, X, &x);
536:
538: for (j=yints; j<yinte; j++) {
539: for (i=xints; i<xinte; i++) {
540: norm[0] = PetscMax(norm[0],x[j][i].U);
541: norm[1] = PetscMax(norm[1],x[j][i].F);
542: norm[2] = PetscMax(norm[2],x[j][i].phi);
543: norm[3] = PetscMax(norm[3],x[j][i].psi);
544: }
545: }
547: DAVecRestoreArray(da,X,&x);
548:
550: PetscObjectGetComm((PetscObject)da, &comm);
551:
552: MPI_Allreduce(norm, gnorm, 4, MPI_DOUBLE, MPI_MAX, comm);
553:
555: PetscFPrintf(PETSC_COMM_WORLD, stderr,
556: "%g\t%g\t%g\t%g\t%g\n",
557: t, gnorm[0], gnorm[1], gnorm[2], gnorm[3]);
558:
560: return(0);
561: }
565: int FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
566: {
567: AppCtx *user = (AppCtx*)ptr;
568: TstepCtx *tsCtx = user->tsCtx;
569: Parameter *param = user->param;
570: int ierr,i,j;
571: int xints,xinte,yints,yinte;
572: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
573: PassiveReal de2,rhos2,nu,eta,dde2;
574: PassiveReal two = 2.0,one = 1.0,p5 = 0.5;
575: PassiveReal F_eq_x,By_eq;
577: PetscScalar xx;
578: PetscScalar vx,vy,avx,avy,vxp,vxm,vyp,vym;
579: PetscScalar Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
582: de2 = sqr(user->param->d_e);
583: rhos2 = sqr(user->param->rho_s);
584: nu = user->param->nu;
585: eta = user->param->eta;
586: dde2 = one/de2;
588: /*
589: Define mesh intervals ratios for uniform grid.
590: [Note: FD formulae below are normalized by multiplying through by
591: local volume element to obtain coefficients O(1) in two dimensions.]
592: */
593: dhx = info->mx/lx; dhy = info->my/ly;
594: hx = one/dhx; hy = one/dhy;
595: hxdhy = hx*dhy; hydhx = hy*dhx;
596: hxhy = hx*hy; dhxdhy = dhx*dhy;
598: xints = info->xs; xinte = info->xs+info->xm;
599: yints = info->ys; yinte = info->ys+info->ym;
601: /* Compute over the interior points */
602: for (j=yints; j<yinte; j++) {
603: for (i=xints; i<xinte; i++) {
604: #ifdef EQ
605: xx = i * hx;
606: F_eq_x = - (1. + de2) * sin(xx);
607: By_eq = sin(xx);
608: #else
609: F_eq_x = 0.;
610: By_eq = 0.;
611: #endif
613: /*
614: * convective coefficients for upwinding
615: */
617: vx = - D_y(x,phi,i,j);
618: vy = D_x(x,phi,i,j);
619: avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
620: avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
621: #ifndef UPWINDING
622: vxp = vxm = p5*vx;
623: vyp = vym = p5*vy;
624: #endif
626: Bx = D_y(x,psi,i,j);
627: By = - D_x(x,psi,i,j) + By_eq;
628: aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
629: aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
630: #ifndef UPWINDING
631: Bxp = Bxm = p5*Bx;
632: Byp = Bym = p5*By;
633: #endif
635: /* Lap(phi) - U */
636: f[j][i].phi = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
638: /* psi - d_e^2 * Lap(psi) - F */
639: f[j][i].psi = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
641: /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
642: - nu Lap(U) */
643: f[j][i].U =
644: ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
645: vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
646: (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
647: Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
648: nu * Lapl(x,U,i,j)) * hxhy;
649:
650: /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
651: - eta * Lap(psi) */
652: f[j][i].F =
653: ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
654: vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
655: (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
656: Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
657: eta * Lapl(x,psi,i,j)) * hxhy;
658: }
659: }
661: /* Add time step contribution */
662: if (tsCtx->ires) {
663: if ((param->second_order) && (tsCtx->itstep > 0))
664: {
665: AddTSTermLocal2(info,x,f,user);
666:
667: }
668: else
669: {
670: AddTSTermLocal(info,x,f,user);
671:
672: }
673: }
675: /*
676: Flop count (multiply-adds are counted as 2 operations)
677: */
678: /* PetscLogFlops(84*info->ym*info->xm); FIXME */
679: return(0);
680: }
682: /*---------------------------------------------------------------------*/
685: int Update(DMMG *dmmg)
686: /*---------------------------------------------------------------------*/
687: {
688: AppCtx *user = (AppCtx *) ((dmmg[0])->user);
689: TstepCtx *tsCtx = user->tsCtx;
690: Parameter *param = user->param;
691: SNES snes;
692: int ierr,its,lits,i;
693: int max_steps;
694: int nfailsCum = 0,nfails = 0;
695: static int ic_out;
696: PetscTruth ts_monitor = (tsCtx->ts_monitor && !user->param->PreLoading) ? PETSC_TRUE : PETSC_FALSE;
700: if (user->param->PreLoading)
701: max_steps = 1;
702: else
703: max_steps = tsCtx->max_steps;
704:
705: Initialize(dmmg);
706:
708: for (tsCtx->itstep = 0; tsCtx->itstep < max_steps; tsCtx->itstep++) {
710: if ((param->second_order) && (tsCtx->itstep > 0))
711: {
712: for (i=param->mglevels-1; i>=0 ;i--)
713: {
714: VecCopy(((AppCtx*)dmmg[i]->user)->Xold,
715: ((AppCtx*)dmmg[i]->user)->Xoldold);
716:
717: }
718: }
720: for (i=param->mglevels-1; i>0 ;i--) {
721: MatRestrict(dmmg[i]->R, dmmg[i]->x, dmmg[i-1]->x);
722:
723: VecPointwiseMult(dmmg[i]->Rscale,dmmg[i-1]->x,dmmg[i-1]->x);
724:
725: VecCopy(dmmg[i]->x, ((AppCtx*)dmmg[i]->user)->Xold);
726:
727: }
728: VecCopy(dmmg[0]->x, ((AppCtx*)dmmg[0]->user)->Xold);
729:
731: DMMGSolve(dmmg);
732:
734: snes = DMMGGetSNES(dmmg);
737: if (tsCtx->itstep == 665000)
738: {
739: KSP ksp;
740: PC pc;
741: Mat mat, pmat;
742: MatStructure flag;
743: PetscViewer viewer;
744: char file[128];
745: KSP ksp;
747: PetscStrcpy(file, "matrix");
749: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file,
750: PETSC_FILE_CREATE, &viewer);
751:
753: SNESGetKSP(snes, &ksp);
754:
756: KSPGetPC(ksp, &pc);
757:
759: PCGetOperators(pc, &mat, &pmat, &flag);
760:
762: MatView(mat, viewer);
763:
765: PetscViewerDestroy(viewer);
766:
767: SETERRQ(1,"Done saving Jacobian");
768: }
771: tsCtx->t += tsCtx->dt;
773: /* save restart solution if requested at a particular time, then exit */
774: if (tsCtx->dump_time > 0.0 && tsCtx->t >= tsCtx->dump_time) {
775: Vec v = DMMGGetx(dmmg);
776: VecView(v,PETSC_VIEWER_BINARY_WORLD);
777: SETERRQ1(1,"Saved solution at time %g",tsCtx->t);
778: }
780: if (ts_monitor)
781: {
782: SNESGetIterationNumber(snes, &its);
783: SNESGetNumberLinearIterations(snes, &lits);
784: SNESGetNumberUnsuccessfulSteps(snes, &nfails);
785: SNESGetFunctionNorm(snes, &tsCtx->fnorm);
787: nfailsCum += nfails;
788: if (nfailsCum >= 2)
789: SETERRQ(1, "unable to find a newton step");
791: PetscPrintf(PETSC_COMM_WORLD,
792: "time step = %d, time = %g, number of nonlinear steps = %d, "
793: "number of linear steps = %d, norm of the function = %g\n",
794: tsCtx->itstep + 1, tsCtx->t, its, lits, tsCtx->fnorm);
795:
797: /* send solution over to Matlab, to be visualized (using ex29.m) */
798: if (!user->param->PreLoading && tsCtx->socketviewer)
799: {
800: Vec v;
801: SNESGetSolution(snes, &v);
802: VecView(v, tsCtx->socketviewer);
803: }
804: }
806: if (!param->PreLoading) {
807: if (param->draw_contours) {
808: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
809:
810: }
812: if (1 && ts_monitor) {
813: /* compute maxima */
814: ComputeMaxima((DA) dmmg[param->mglevels-1]->dm,
815: dmmg[param->mglevels-1]->x,
816: tsCtx->t);
817: /* output */
818: if (ic_out++ == (int)(tsCtx->dt_out / tsCtx->dt + .5)) {
819: #ifdef HAVE_DA_HDF
820: char fname[128];
822: sprintf(fname, "out-%g.hdf", tsCtx->t);
823: DAVecHDFOutput(DMMGGetDA(dmmg), DMMGGetx(dmmg), fname);
824:
825: #else
826: /*
827: Gnuplot(DMMGGetDA(dmmg), DMMGGetx(dmmg), tsCtx->t);
828:
829: */
830: #endif
831: ic_out = 1;
832: }
833: }
834: }
835: } /* End of time step loop */
836:
837: if (!user->param->PreLoading){
838: SNESGetFunctionNorm(snes,&tsCtx->fnorm);
839: PetscPrintf(PETSC_COMM_WORLD, "timesteps %d fnorm = %g\n",
840: tsCtx->itstep, tsCtx->fnorm);
841:
842: }
844: if (user->param->PreLoading) {
845: tsCtx->fnorm_ini = 0.0;
846: }
847:
848: return(0);
849: }
851: /*---------------------------------------------------------------------*/
854: int AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,AppCtx *user)
855: /*---------------------------------------------------------------------*/
856: {
857: TstepCtx *tsCtx = user->tsCtx;
858: DA da = info->da;
859: int ierr,i,j;
860: int xints,xinte,yints,yinte;
861: PassiveReal hx,hy,dhx,dhy,hxhy;
862: PassiveReal one = 1.0;
863: PassiveScalar dtinv;
864: PassiveField **xold;
868: xints = info->xs; xinte = info->xs+info->xm;
869: yints = info->ys; yinte = info->ys+info->ym;
871: dhx = info->mx/lx; dhy = info->my/ly;
872: hx = one/dhx; hy = one/dhy;
873: hxhy = hx*hy;
875: dtinv = hxhy/(tsCtx->dt);
877: DAVecGetArray(da,user->Xold,&xold);
878: for (j=yints; j<yinte; j++) {
879: for (i=xints; i<xinte; i++) {
880: f[j][i].U += dtinv*(x[j][i].U-xold[j][i].U);
881: f[j][i].F += dtinv*(x[j][i].F-xold[j][i].F);
882: }
883: }
884: DAVecRestoreArray(da,user->Xold,&xold);
886: return(0);
887: }
889: /*---------------------------------------------------------------------*/
892: int AddTSTermLocal2(DALocalInfo* info,Field **x,Field **f,AppCtx *user)
893: /*---------------------------------------------------------------------*/
894: {
895: TstepCtx *tsCtx = user->tsCtx;
896: DA da = info->da;
897: int ierr,i,j;
898: int xints,xinte,yints,yinte;
899: PassiveReal hx,hy,dhx,dhy,hxhy;
900: PassiveReal one = 1.0, onep5 = 1.5, two = 2.0, p5 = 0.5;
901: PassiveScalar dtinv;
902: PassiveField **xoldold,**xold;
906: xints = info->xs; xinte = info->xs+info->xm;
907: yints = info->ys; yinte = info->ys+info->ym;
909: dhx = info->mx/lx; dhy = info->my/ly;
910: hx = one/dhx; hy = one/dhy;
911: hxhy = hx*hy;
913: dtinv = hxhy/(tsCtx->dt);
915: DAVecGetArray(da,user->Xoldold,&xoldold);
916: DAVecGetArray(da,user->Xold,&xold);
917: for (j=yints; j<yinte; j++) {
918: for (i=xints; i<xinte; i++) {
919: f[j][i].U += dtinv * (onep5 * x[j][i].U - two * xold[j][i].U +
920: p5 * xoldold[j][i].U);
921: f[j][i].F += dtinv * (onep5 * x[j][i].F - two * xold[j][i].F +
922: p5 * xoldold[j][i].F);
923: }
924: }
925: DAVecRestoreArray(da,user->Xoldold,&xoldold);
926: DAVecRestoreArray(da,user->Xold,&xold);
928: return(0);
929: }
931: /* -------------------------------------------------------------------------*/
934: int AttachNullSpace(PC pc,Vec model)
935: {
936: int i,ierr,rstart,rend,N;
937: MatNullSpace sp;
938: Vec v,vs[1];
939: PetscScalar *vx,scale;
942: VecDuplicate(model,&v);
943: VecGetSize(model,&N);
944: scale = 2.0/sqrt(N);
945: VecGetArray(v,&vx);
946: VecGetOwnershipRange(v,&rstart,&rend);
947: for (i=rstart; i<rend; i++) {
948: if (!(i % 4)) vx[i-rstart] = scale;
949: else vx[i-rstart] = 0.0;
950: }
951: VecRestoreArray(v,&vx);
952: vs[0] = v;
953: MatNullSpaceCreate(PETSC_COMM_WORLD,0,1,vs,&sp);
954: VecDestroy(v);
955: PCNullSpaceAttach(pc,sp);
956: MatNullSpaceDestroy(sp);
957: return(0);
958: }
960: /*
961: This is an experimental function and can be safely ignored.
962: */
963: int FormFunctionLocali(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
964: {
965: AppCtx *user = (AppCtx*)ptr;
966: TstepCtx *tsCtx = user->tsCtx;
967: int ierr,i,j,c;
968: int xints,xinte,yints,yinte;
969: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
970: PassiveReal de2,rhos2,nu,eta,dde2;
971: PassiveReal two = 2.0,one = 1.0,p5 = 0.5;
972: PassiveReal F_eq_x,By_eq,dtinv;
974: PetscScalar xx;
975: PetscScalar vx,vy,avx,avy,vxp,vxm,vyp,vym;
976: PetscScalar Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
977: PassiveField **xold;
980: de2 = sqr(user->param->d_e);
981: rhos2 = sqr(user->param->rho_s);
982: nu = user->param->nu;
983: eta = user->param->eta;
984: dde2 = one/de2;
986: /*
987: Define mesh intervals ratios for uniform grid.
988: [Note: FD formulae below are normalized by multiplying through by
989: local volume element to obtain coefficients O(1) in two dimensions.]
990: */
991: dhx = info->mx/lx; dhy = info->my/ly;
992: hx = one/dhx; hy = one/dhy;
993: hxdhy = hx*dhy; hydhx = hy*dhx;
994: hxhy = hx*hy; dhxdhy = dhx*dhy;
996: xints = info->xs; xinte = info->xs+info->xm;
997: yints = info->ys; yinte = info->ys+info->ym;
1000: i = st->i; j = st->j; c = st->c;
1002: #ifdef EQ
1003: xx = i * hx;
1004: F_eq_x = - (1. + de2) * sin(xx);
1005: By_eq = sin(xx);
1006: #else
1007: F_eq_x = 0.;
1008: By_eq = 0.;
1009: #endif
1011: /*
1012: * convective coefficients for upwinding
1013: */
1015: vx = - D_y(x,phi,i,j);
1016: vy = D_x(x,phi,i,j);
1017: avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
1018: avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
1019: #ifndef UPWINDING
1020: vxp = vxm = p5*vx;
1021: vyp = vym = p5*vy;
1022: #endif
1024: Bx = D_y(x,psi,i,j);
1025: By = - D_x(x,psi,i,j) + By_eq;
1026: aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
1027: aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
1028: #ifndef UPWINDING
1029: Bxp = Bxm = p5*Bx;
1030: Byp = Bym = p5*By;
1031: #endif
1033: DAVecGetArray(info->da,user->Xold,&xold);
1034: dtinv = hxhy/(tsCtx->dt);
1035: switch(c) {
1037: case 0:
1038: /* Lap(phi) - U */
1039: *f = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
1040: break;
1042: case 1:
1043: /* psi - d_e^2 * Lap(psi) - F */
1044: *f = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
1045: break;
1047: case 2:
1048: /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
1049: - nu Lap(U) */
1050: *f =
1051: ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
1052: vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
1053: (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
1054: Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
1055: nu * Lapl(x,U,i,j)) * hxhy;
1056: *f += dtinv*(x[j][i].U-xold[j][i].U);
1057: break;
1059: case 3:
1060: /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
1061: - eta * Lap(psi) */
1062: *f =
1063: ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
1064: vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
1065: (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
1066: Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
1067: eta * Lapl(x,psi,i,j)) * hxhy;
1068: *f += dtinv*(x[j][i].F-xold[j][i].F);
1069: break;
1070: }
1071: DAVecRestoreArray(info->da,user->Xold,&xold);
1074: /*
1075: Flop count (multiply-adds are counted as 2 operations)
1076: */
1077: /* PetscLogFlops(84*info->ym*info->xm); FIXME */
1078: return(0);
1079: }