Actual source code: ex29.c
2: /*#define HAVE_DA_HDF*/
4: /* solve the equations for the perturbed magnetic field only */
5: #define EQ
7: /* turning this on causes instability?!? */
8: /* #define UPWINDING */
10: static char help[] = "Hall MHD with in two dimensions with time stepping and multigrid.\n\n\
11: -options_file ex29.options\n\
12: other PETSc options\n\
13: -resistivity <eta>\n\
14: -viscosity <nu>\n\
15: -skin_depth <d_e>\n\
16: -larmor_radius <rho_s>\n\
17: -contours : draw contour plots of solution\n\n";
19: /*T
20: Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
21: Concepts: DA^using distributed arrays;
22: Concepts: multicomponent
23: Processors: n
24: T*/
26: /* ------------------------------------------------------------------------
28: We thank Kai Germaschewski for contributing the FormFunctionLocal()
30: ------------------------------------------------------------------------- */
32: /*
33: Include "petscda.h" so that we can use distributed arrays (DAs).
34: Include "petscsnes.h" so that we can use SNES solvers.
35: Include "petscmg.h" to control the multigrid solvers.
36: Note that these automatically include:
37: petsc.h - base PETSc routines petscvec.h - vectors
38: petscsys.h - system routines petscmat.h - matrices
39: petscis.h - index sets petscksp.h - Krylov subspace methods
40: petscviewer.h - viewers petscpc.h - preconditioners
41: petscksp.h - linear solvers
42: */
43: #include petscsnes.h
44: #include petscda.h
45: #include petscmg.h
46: #include petscdmmg.h
48: #ifdef HAVE_DA_HDF
49: PetscInt DAVecHDFOutput(DA,Vec,char*);
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.*3.1415926)
62: #define ly (4.*3.1415926)
63: #define sqr(a) ((a)*(a))
65: /*
66: User-defined routines and data structures
67: */
69: typedef struct {
70: PetscReal fnorm_ini,dt_ini;
71: PetscReal dt,dt_out;
72: PetscReal ptime;
73: PetscReal max_time;
74: PetscReal fnorm_ratio;
75: PetscInt ires,itstep;
76: PetscInt max_steps,print_freq;
77: PetscReal t;
78: PetscScalar fnorm;
80: PetscTruth ts_monitor; /* print information about each time step */
81: PetscReal dump_time; /* time to dump solution to a file */
82: PetscViewer socketviewer; /* socket to dump solution at each timestep for visualization */
83: } TstepCtx;
85: typedef struct {
86: PetscScalar phi,psi,U,F;
87: } Field;
89: typedef struct {
90: PassiveScalar phi,psi,U,F;
91: } PassiveField;
93: typedef struct {
94: PetscInt mglevels;
95: PetscInt cycles; /* number of time steps for integration */
96: PassiveReal nu,eta,d_e,rho_s; /* physical parameters */
97: PetscTruth draw_contours; /* flag - 1 indicates drawing contours */
98: PetscTruth second_order;
99: PetscTruth PreLoading;
100: } Parameters;
102: typedef struct {
103: Vec Xoldold,Xold;
104: TstepCtx *tsCtx;
105: Parameters *param;
106: } AppCtx;
120: int main(int argc,char **argv)
121: {
123: DMMG *dmmg; /* multilevel grid structure */
124: AppCtx *user; /* user-defined work context (one for each level) */
125: TstepCtx tsCtx; /* time-step parameters (one total) */
126: Parameters param; /* physical parameters (one total) */
127: PetscInt i,m,n,mx,my;
128: DA da;
129: PetscReal dt_ratio;
130: PetscInt dfill[16] = {1,0,1,0,
131: 0,1,0,1,
132: 0,0,1,1,
133: 0,1,1,1};
134: PetscInt 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);
142: PreLoadBegin(PETSC_TRUE,"SetUp");
144: param.PreLoading = PreLoading;
145: DMMGCreate(PETSC_COMM_WORLD,1,&user,&dmmg);
146: param.mglevels = DMMGGetLevels(dmmg);
148: /*
149: Create distributed array multigrid object (DMMG) to manage
150: parallel grid and vectors for principal unknowns (x) and
151: governing residuals (f)
152: */
153: DACreate2d(PETSC_COMM_WORLD, DA_XYPERIODIC, DA_STENCIL_STAR, -5, -5,
154: PETSC_DECIDE, PETSC_DECIDE, 4, 1, 0, 0, &da);
156: /* overwrite the default sparsity pattern toone specific for
157: this code's nonzero structure */
158: DASetBlockFills(da,dfill,ofill);
160: DMMGSetDM(dmmg,(DM)da);
161: DADestroy(da);
163: /* default physical parameters */
164: param.nu = 0;
165: param.eta = 1e-3;
166: param.d_e = 0.2;
167: param.rho_s = 0;
168: param.second_order = PETSC_FALSE;
170: PetscOptionsGetReal(PETSC_NULL, "-viscosity", ¶m.nu,PETSC_NULL);
172: PetscOptionsGetReal(PETSC_NULL, "-resistivity", ¶m.eta,PETSC_NULL);
174: PetscOptionsGetReal(PETSC_NULL, "-skin_depth", ¶m.d_e,PETSC_NULL);
176: PetscOptionsGetReal(PETSC_NULL, "-larmor_radius", ¶m.rho_s,PETSC_NULL);
178: PetscOptionsHasName(PETSC_NULL, "-contours", ¶m.draw_contours);
180: PetscOptionsHasName(PETSC_NULL, "-second_order", ¶m.second_order);
182: DASetFieldName(DMMGGetDA(dmmg), 0, "phi");
184: DASetFieldName(DMMGGetDA(dmmg), 1, "psi");
186: DASetFieldName(DMMGGetDA(dmmg), 2, "U");
188: DASetFieldName(DMMGGetDA(dmmg), 3, "F");
190: /*======================================================================*/
191: /* Initialize stuff related to time stepping */
192: /*======================================================================*/
194: DAGetInfo(DMMGGetDA(dmmg),0,&mx,&my,0,0,0,0,0,0,0,0);
196: tsCtx.fnorm_ini = 0;
197: tsCtx.max_steps = 1000000;
198: tsCtx.max_time = 1.0e+12;
199: /* use for dt = min(dx,dy); multiplied by dt_ratio below */
200: tsCtx.dt = PetscMin(lx/mx,ly/my);
201: tsCtx.fnorm_ratio = 1.0e+10;
202: tsCtx.t = 0;
203: tsCtx.dt_out = 10;
204: tsCtx.print_freq = tsCtx.max_steps;
205: tsCtx.ts_monitor = PETSC_FALSE;
206: tsCtx.dump_time = -1.0;
208: PetscOptionsGetInt(PETSC_NULL, "-max_st", &tsCtx.max_steps,PETSC_NULL);
209: PetscOptionsGetReal(PETSC_NULL, "-ts_rtol", &tsCtx.fnorm_ratio,PETSC_NULL);
210: PetscOptionsGetInt(PETSC_NULL, "-print_freq", &tsCtx.print_freq,PETSC_NULL);
212: dt_ratio = 1.0;
213: PetscOptionsGetReal(PETSC_NULL, "-dt_ratio", &dt_ratio,PETSC_NULL);
214: tsCtx.dt *= dt_ratio;
216: PetscOptionsHasName(PETSC_NULL, "-ts_monitor", &tsCtx.ts_monitor);
217: PetscOptionsGetReal(PETSC_NULL, "-dump", &tsCtx.dump_time,PETSC_NULL);
220: tsCtx.socketviewer = 0;
221: #if defined(PETSC_USE_SOCKET_VIEWER)
222: {
223: PetscTruth flg;
224: PetscOptionsHasName(PETSC_NULL, "-socket_viewer", &flg);
225: if (flg && !PreLoading) {
226: tsCtx.ts_monitor = PETSC_TRUE;
227: PetscViewerSocketOpen(PETSC_COMM_WORLD,0,PETSC_DECIDE,&tsCtx.socketviewer);
228: }
229: }
230: #endif
232: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
233: Create user context, set problem data, create vector data structures.
234: Also, compute the initial guess.
235: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
236: /* create application context for each level */
237: PetscMalloc(param.mglevels*sizeof(AppCtx),&user);
238: for (i=0; i<param.mglevels; i++) {
239: /* create work vectors to hold previous time-step solution and
240: function value */
241: VecDuplicate(dmmg[i]->x, &user[i].Xoldold);
242: VecDuplicate(dmmg[i]->x, &user[i].Xold);
243: user[i].tsCtx = &tsCtx;
244: user[i].param = ¶m;
245: DMMGSetUser(dmmg,i,&user[i]);
246: }
247:
248: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
249: Create nonlinear solver context
250:
251: Process adiC(20): AddTSTermLocal AddTSTermLocal2 FormFunctionLocal FormFunctionLocali
252: Process blockadiC(4): FormFunctionLocali4
253: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
254: DMMGSetSNESLocal(dmmg, FormFunctionLocal, 0,ad_FormFunctionLocal, admf_FormFunctionLocal);
255: DMMGSetSNESLocali(dmmg,FormFunctionLocali,0,admf_FormFunctionLocali);
256: DMMGSetSNESLocalib(dmmg,FormFunctionLocali4,0,admfb_FormFunctionLocali4);
257:
258: /* attach nullspace to each level of the preconditioner */
259: DMMGSetNullSpace(dmmg,PETSC_FALSE,1,CreateNullSpace);
260:
261: PetscPrintf(PETSC_COMM_WORLD, "finish setupNull!");
263: if (PreLoading) {
264: PetscPrintf(PETSC_COMM_WORLD, "# viscosity = %G, resistivity = %G, "
265: "skin_depth # = %G, larmor_radius # = %G\n",
266: param.nu, param.eta, param.d_e, param.rho_s);
267: DAGetInfo(DMMGGetDA(dmmg),0,&m,&n,0,0,0,0,0,0,0,0);
268: PetscPrintf(PETSC_COMM_WORLD,"Problem size %D by %D\n",m,n);
269: PetscPrintf(PETSC_COMM_WORLD,"dx %G dy %G dt %G ratio dt/min(dx,dy) %G\n",lx/mx,ly/my,tsCtx.dt,dt_ratio);
270: }
272:
273:
274: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
275: Solve the nonlinear system
276: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
277:
278: PreLoadStage("Solve");
280: if (param.draw_contours) {
281: VecView(((AppCtx*)DMMGGetUser(dmmg,param.mglevels-1))->Xold,PETSC_VIEWER_DRAW_WORLD);
282: }
284: Update(dmmg);
286: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
287: Free work space. All PETSc objects should be destroyed when they
288: are no longer needed.
289: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
290:
291: for (i=0; i<param.mglevels; i++) {
292: VecDestroy(user[i].Xoldold);
293: VecDestroy(user[i].Xold);
294: }
295: PetscFree(user);
296: DMMGDestroy(dmmg);
298: PreLoadEnd();
299:
300: PetscFinalize();
301: return 0;
302: }
304: /* ------------------------------------------------------------------- */
307: /* ------------------------------------------------------------------- */
308: PetscErrorCode Gnuplot(DA da, Vec X, double mtime)
309: {
310: PetscInt i,j,xs,ys,xm,ym;
311: PetscInt xints,xinte,yints,yinte;
313: Field **x;
314: FILE *f;
315: char fname[PETSC_MAX_PATH_LEN];
316: PetscMPIInt rank;
318: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
320: sprintf(fname, "out-%g-%d.dat", mtime, rank);
322: f = fopen(fname, "w");
324: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
326: DAVecGetArray(da,X,&x);
328: xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;
330: for (j=yints; j<yinte; j++) {
331: for (i=xints; i<xinte; i++) {
332: PetscFPrintf(PETSC_COMM_WORLD, f,
333: "%D %D %G %G %G %G %G %G\n",
334: i, j, 0.0, 0.0,
335: PetscAbsScalar(x[j][i].U), PetscAbsScalar(x[j][i].F),
336: PetscAbsScalar(x[j][i].phi), PetscAbsScalar(x[j][i].psi));
337: }
338: PetscFPrintf(PETSC_COMM_WORLD,f, "\n");
339: }
340: DAVecRestoreArray(da,X,&x);
341: fclose(f);
342: return 0;
343: }
345: /* ------------------------------------------------------------------- */
348: /* ------------------------------------------------------------------- */
349: PetscErrorCode Initialize(DMMG *dmmg)
350: {
351: AppCtx *appCtx = (AppCtx*)DMMGGetUser(dmmg,0);
352: Parameters *param = appCtx->param;
353: DA da;
355: PetscInt i,j,mx,my,xs,ys,xm,ym;
356: PetscReal two = 2.0,one = 1.0;
357: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
358: PetscReal d_e,rho_s,de2,xx,yy;
359: Field **x, **localx;
360: Vec localX;
361: PetscTruth flg;
364: PetscOptionsHasName(0,"-restart",&flg);
365: if (flg) {
366: PetscViewer viewer;
367: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"binaryoutput",FILE_MODE_READ,&viewer);
368: VecLoadIntoVector(viewer,dmmg[param->mglevels-1]->x);
369: PetscViewerDestroy(viewer);
370: return(0);
371: }
373: d_e = param->d_e;
374: rho_s = param->rho_s;
375: de2 = sqr(param->d_e);
377: da = (DA)(dmmg[param->mglevels-1]->dm);
378: DAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0);
380: dhx = mx/lx; dhy = my/ly;
381: hx = one/dhx; hy = one/dhy;
382: hxdhy = hx*dhy; hydhx = hy*dhx;
383: hxhy = hx*hy; dhxdhy = dhx*dhy;
385: /*
386: Get local grid boundaries (for 2-dimensional DA):
387: xs, ys - starting grid indices (no ghost points)
388: xm, ym - widths of local grid (no ghost points)
389: */
390: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
392: DAGetLocalVector(da,&localX);
393: /*
394: Get a pointer to vector data.
395: - For default PETSc vectors, VecGetArray() returns a pointer to
396: the data array. Otherwise, the routine is implementation dependent.
397: - You MUST call VecRestoreArray() when you no longer need access to
398: the array.
399: */
400: DAVecGetArray(da,dmmg[param->mglevels-1]->x,&x);
401: DAVecGetArray(da,localX,&localx);
403: /*
404: Compute initial solution over the locally owned part of the grid
405: */
406: #if defined (PETSC_HAVE_ERF)
407: {
408: PetscReal eps = lx/ly;
409: PetscReal pert = 1e-4;
410: PetscReal k = 1.*eps;
411: PetscReal gam;
413: if (d_e < rho_s) d_e = rho_s;
414: gam = k * d_e;
415: for (j=ys-1; j<ys+ym+1; j++) {
416: yy = j * hy;
417: for (i=xs-1; i<xs+xm+1; i++) {
418: xx = i * hx;
419: if (xx < -3.1416926/2) {
420: localx[j][i].phi = pert * gam / k * erf((xx + 3.1415926) / (sqrt(2.0) * d_e)) * (-sin(k*yy));
421: } else if (xx < 3.1415926/2) {
422: localx[j][i].phi = - pert * gam / k * erf(xx / (sqrt(2.0) * d_e)) * (-sin(k*yy));
423: } else if (xx < 3*3.1415926/2){
424: localx[j][i].phi = pert * gam / k * erf((xx - 3.1415926) / (sqrt(2.0) * d_e)) * (-sin(k*yy));
425: } else {
426: localx[j][i].phi = - pert * gam / k * erf((xx - 2.*3.1415926) / (sqrt(2.0) * d_e)) * (-sin(k*yy));
427: }
428: #ifdef EQ
429: localx[j][i].psi = 0.;
430: #else
431: localx[j][i].psi = cos(xx);
432: #endif
433: }
434: }
436: for (j=ys; j<ys+ym; j++) {
437: for (i=xs; i<xs+xm; i++) {
438: x[j][i].U = Lapl(localx,phi,i,j);
439: x[j][i].F = localx[j][i].psi - de2 * Lapl(localx,psi,i,j);
440: x[j][i].phi = localx[j][i].phi;
441: x[j][i].psi = localx[j][i].psi;
442: }
443: }
444: }
445: #else
446: SETERRQ(1,"erf() not available on this machine");
447: #endif
449: /*
450: Restore vector
451: */
452: DAVecRestoreArray(da,dmmg[param->mglevels-1]->x,&x);
453:
454: DAVecRestoreArray(da,localX,&localx);
455:
456: DARestoreLocalVector(da,&localX);
457:
459: return(0);
460: }
464: PetscErrorCode ComputeMaxima(DA da, Vec X, PetscReal t)
465: {
467: PetscInt i,j,mx,my,xs,ys,xm,ym;
468: PetscInt xints,xinte,yints,yinte;
469: Field **x;
470: double norm[4] = {0,0,0,0};
471: double gnorm[4];
472: MPI_Comm comm;
476: DAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0);
478: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
479:
481: xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;
483: DAVecGetArray(da, X, &x);
485: for (j=yints; j<yinte; j++) {
486: for (i=xints; i<xinte; i++) {
487: norm[0] = PetscMax(norm[0],PetscAbsScalar(x[j][i].U));
488: norm[1] = PetscMax(norm[1],PetscAbsScalar(x[j][i].F));
489: norm[2] = PetscMax(norm[2],PetscAbsScalar(x[j][i].phi));
490: norm[3] = PetscMax(norm[3],PetscAbsScalar(x[j][i].psi));
491: }
492: }
494: DAVecRestoreArray(da,X,&x);
496: PetscObjectGetComm((PetscObject)da, &comm);
497: MPI_Allreduce(norm, gnorm, 4, MPI_DOUBLE, MPI_MAX, comm);
499: PetscFPrintf(PETSC_COMM_WORLD, stderr,"%G\t%G\t%G\t%G\t%G\n",t, gnorm[0], gnorm[1], gnorm[2], gnorm[3]);
501: return(0);
502: }
506: PetscErrorCode FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
507: {
508: AppCtx *user = (AppCtx*)ptr;
509: TstepCtx *tsCtx = user->tsCtx;
510: Parameters *param = user->param;
512: PetscInt xints,xinte,yints,yinte,i,j;
513: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
514: PassiveReal de2,rhos2,nu,eta,dde2;
515: PassiveReal two = 2.0,one = 1.0,p5 = 0.5;
516: PassiveReal F_eq_x,By_eq;
517: PetscScalar xx;
518: PetscScalar vx,vy,avx,avy,vxp,vxm,vyp,vym;
519: PetscScalar Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
522: de2 = sqr(param->d_e);
523: rhos2 = sqr(param->rho_s);
524: nu = param->nu;
525: eta = param->eta;
526: dde2 = one/de2;
528: /*
529: Define mesh intervals ratios for uniform grid.
530: [Note: FD formulae below are normalized by multiplying through by
531: local volume element to obtain coefficients O(1) in two dimensions.]
532: */
533: dhx = info->mx/lx; dhy = info->my/ly;
534: hx = one/dhx; hy = one/dhy;
535: hxdhy = hx*dhy; hydhx = hy*dhx;
536: hxhy = hx*hy; dhxdhy = dhx*dhy;
538: xints = info->xs; xinte = info->xs+info->xm;
539: yints = info->ys; yinte = info->ys+info->ym;
541: /* Compute over the interior points */
542: for (j=yints; j<yinte; j++) {
543: for (i=xints; i<xinte; i++) {
544: #ifdef EQ
545: xx = i * hx;
546: F_eq_x = - (1. + de2) * sin(PetscAbsScalar(xx));
547: By_eq = sin(PetscAbsScalar(xx));
548: #else
549: F_eq_x = 0.;
550: By_eq = 0.;
551: #endif
553: /*
554: * convective coefficients for upwinding
555: */
557: vx = - D_y(x,phi,i,j);
558: vy = D_x(x,phi,i,j);
559: avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
560: avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
561: #ifndef UPWINDING
562: vxp = vxm = p5*vx;
563: vyp = vym = p5*vy;
564: #endif
566: Bx = D_y(x,psi,i,j);
567: By = - D_x(x,psi,i,j) + By_eq;
568: aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
569: aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
570: #ifndef UPWINDING
571: Bxp = Bxm = p5*Bx;
572: Byp = Bym = p5*By;
573: #endif
575: /* Lap(phi) - U */
576: f[j][i].phi = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
578: /* psi - d_e^2 * Lap(psi) - F */
579: f[j][i].psi = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
581: /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
582: - nu Lap(U) */
583: f[j][i].U =
584: ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
585: vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
586: (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
587: Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
588: nu * Lapl(x,U,i,j)) * hxhy;
589:
590: /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
591: - eta * Lap(psi) */
592: f[j][i].F =
593: ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
594: vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
595: (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
596: Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
597: eta * Lapl(x,psi,i,j)) * hxhy;
598: }
599: }
601: /* Add time step contribution */
602: if (tsCtx->ires) {
603: if ((param->second_order) && (tsCtx->itstep > 0)){
604: AddTSTermLocal2(info,x,f,user);
605: } else {
606: AddTSTermLocal(info,x,f,user);
607: }
608: }
610: /*
611: Flop count (multiply-adds are counted as 2 operations)
612: */
613: /* PetscLogFlops(84*info->ym*info->xm); FIXME */
614: return(0);
615: }
617: /*---------------------------------------------------------------------*/
620: PetscErrorCode Update(DMMG *dmmg)
621: /*---------------------------------------------------------------------*/
622: {
623: AppCtx *user = (AppCtx *) DMMGGetUser(dmmg,0);
624: TstepCtx *tsCtx = user->tsCtx;
625: Parameters *param = user->param;
626: SNES snes;
627: PetscErrorCode ierr;
628: PetscInt its,lits,i;
629: PetscInt max_steps;
630: PetscInt nfailsCum = 0,nfails = 0;
631: static PetscInt ic_out;
632: PetscTruth ts_monitor = (tsCtx->ts_monitor && !param->PreLoading) ? PETSC_TRUE : PETSC_FALSE;
636: if (param->PreLoading)
637: max_steps = 1;
638: else
639: max_steps = tsCtx->max_steps;
640:
641: Initialize(dmmg);
643: snes = DMMGGetSNES(dmmg);
645: for (tsCtx->itstep = 0; tsCtx->itstep < max_steps; tsCtx->itstep++) {
646:
647: PetscPrintf(PETSC_COMM_WORLD, "time step=%d!\n",tsCtx->itstep);
648: if ((param->second_order) && (tsCtx->itstep > 0))
649: {
650: for (i=param->mglevels-1; i>=0 ;i--)
651: {
652: VecCopy(((AppCtx*)DMMGGetUser(dmmg,i))->Xold,((AppCtx*)DMMGGetUser(dmmg,i))->Xoldold);
653: }
654: }
656: for (i=param->mglevels-1; i>0 ;i--) {
657: MatRestrict(dmmg[i]->R, dmmg[i]->x, dmmg[i-1]->x);
658: VecPointwiseMult(dmmg[i-1]->x,dmmg[i-1]->x,dmmg[i]->Rscale);
659: VecCopy(dmmg[i]->x, ((AppCtx*)DMMGGetUser(dmmg,i))->Xold);
660: }
661: VecCopy(dmmg[0]->x, ((AppCtx*)DMMGGetUser(dmmg,0))->Xold);
663: DMMGSolve(dmmg);
667: if (tsCtx->itstep == 665000)
668: {
669: KSP ksp;
670: PC pc;
671: Mat mat, pmat;
672: MatStructure flag;
673: PetscViewer viewer;
674: char file[PETSC_MAX_PATH_LEN];
676: PetscStrcpy(file, "matrix");
678: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_WRITE, &viewer);
680: SNESGetKSP(snes, &ksp);
682: KSPGetPC(ksp, &pc);
684: PCGetOperators(pc, &mat, &pmat, &flag);
686: MatView(mat, viewer);
688: PetscViewerDestroy(viewer);
689: SETERRQ(1,"Done saving Jacobian");
690: }
693: tsCtx->t += tsCtx->dt;
695: /* save restart solution if requested at a particular time, then exit */
696: if (tsCtx->dump_time > 0.0 && tsCtx->t >= tsCtx->dump_time) {
697: Vec v = DMMGGetx(dmmg);
698: VecView(v,PETSC_VIEWER_BINARY_WORLD);
699: SETERRQ1(1,"Saved solution at time %G",tsCtx->t);
700: }
702: if (ts_monitor)
703: {
704: SNESGetIterationNumber(snes, &its);
705: SNESGetNumberLinearIterations(snes, &lits);
706: SNESGetNumberUnsuccessfulSteps(snes, &nfails);
707: SNESGetFunctionNorm(snes, &tsCtx->fnorm);
709: nfailsCum += nfails;
710: if (nfailsCum >= 2)
711: SETERRQ(1, "unable to find a newton step");
713: PetscPrintf(PETSC_COMM_WORLD,
714: "time step = %D, time = %G, number of nonlinear steps = %D, "
715: "number of linear steps = %D, norm of the function = %G\n",
716: tsCtx->itstep + 1, tsCtx->t, its, lits, PetscAbsScalar(tsCtx->fnorm));
718: /* send solution over to Matlab, to be visualized (using ex29.m) */
719: if (!param->PreLoading && tsCtx->socketviewer)
720: {
721: Vec v;
722: SNESGetSolution(snes, &v);
723: #if defined(PETSC_USE_SOCKET_VIEWER)
724: VecView(v, tsCtx->socketviewer);
725: #endif
726: }
727: }
729: if (!param->PreLoading) {
730: if (param->draw_contours) {
731: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
732: }
734: if (1 && ts_monitor) {
735: /* compute maxima */
736: ComputeMaxima((DA) dmmg[param->mglevels-1]->dm,dmmg[param->mglevels-1]->x,tsCtx->t);
737: /* output */
738: if (ic_out++ == (int)(tsCtx->dt_out / tsCtx->dt + .5)) {
739: #ifdef HAVE_DA_HDF
740: char fname[PETSC_MAX_PATH_LEN];
742: sprintf(fname, "out-%G.hdf", tsCtx->t);
743: DAVecHDFOutput(DMMGGetDA(dmmg), DMMGGetx(dmmg), fname);
744: #else
745: /*
746: Gnuplot(DMMGGetDA(dmmg), DMMGGetx(dmmg), tsCtx->t);
747:
748: */
749: #endif
750: ic_out = 1;
751: }
752: }
753: }
754: } /* End of time step loop */
755:
756: if (!param->PreLoading){
757: SNESGetFunctionNorm(snes,&tsCtx->fnorm);
758: PetscPrintf(PETSC_COMM_WORLD, "timesteps %D fnorm = %G\n",
759: tsCtx->itstep, PetscAbsScalar(tsCtx->fnorm));
760: }
762: if (param->PreLoading) {
763: tsCtx->fnorm_ini = 0.0;
764: }
765:
766: return(0);
767: }
769: /*---------------------------------------------------------------------*/
772: PetscErrorCode AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,AppCtx *user)
773: /*---------------------------------------------------------------------*/
774: {
775: TstepCtx *tsCtx = user->tsCtx;
776: DA da = info->da;
778: PetscInt i,j;
779: PetscInt xints,xinte,yints,yinte;
780: PassiveReal hx,hy,dhx,dhy,hxhy;
781: PassiveReal one = 1.0;
782: PassiveScalar dtinv;
783: PassiveField **xold;
787: xints = info->xs; xinte = info->xs+info->xm;
788: yints = info->ys; yinte = info->ys+info->ym;
790: dhx = info->mx/lx; dhy = info->my/ly;
791: hx = one/dhx; hy = one/dhy;
792: hxhy = hx*hy;
794: dtinv = hxhy/(tsCtx->dt);
796: DAVecGetArray(da,user->Xold,&xold);
797: for (j=yints; j<yinte; j++) {
798: for (i=xints; i<xinte; i++) {
799: f[j][i].U += dtinv*(x[j][i].U-xold[j][i].U);
800: f[j][i].F += dtinv*(x[j][i].F-xold[j][i].F);
801: }
802: }
803: DAVecRestoreArray(da,user->Xold,&xold);
805: return(0);
806: }
808: /*---------------------------------------------------------------------*/
811: PetscErrorCode AddTSTermLocal2(DALocalInfo* info,Field **x,Field **f,AppCtx *user)
812: /*---------------------------------------------------------------------*/
813: {
814: TstepCtx *tsCtx = user->tsCtx;
815: DA da = info->da;
817: PetscInt i,j,xints,xinte,yints,yinte;
818: PassiveReal hx,hy,dhx,dhy,hxhy;
819: PassiveReal one = 1.0, onep5 = 1.5, two = 2.0, p5 = 0.5;
820: PassiveScalar dtinv;
821: PassiveField **xoldold,**xold;
825: xints = info->xs; xinte = info->xs+info->xm;
826: yints = info->ys; yinte = info->ys+info->ym;
828: dhx = info->mx/lx; dhy = info->my/ly;
829: hx = one/dhx; hy = one/dhy;
830: hxhy = hx*hy;
832: dtinv = hxhy/(tsCtx->dt);
834: DAVecGetArray(da,user->Xoldold,&xoldold);
835: DAVecGetArray(da,user->Xold,&xold);
836: for (j=yints; j<yinte; j++) {
837: for (i=xints; i<xinte; i++) {
838: f[j][i].U += dtinv * (onep5 * x[j][i].U - two * xold[j][i].U +
839: p5 * xoldold[j][i].U);
840: f[j][i].F += dtinv * (onep5 * x[j][i].F - two * xold[j][i].F +
841: p5 * xoldold[j][i].F);
842: }
843: }
844: DAVecRestoreArray(da,user->Xoldold,&xoldold);
845: DAVecRestoreArray(da,user->Xold,&xold);
847: return(0);
848: }
850: /* Creates the null space of the Jacobian for a particular level */
853: PetscErrorCode CreateNullSpace(DMMG dmmg,Vec *nulls)
854: {
856: PetscInt i,N,rstart,rend;
857: PetscScalar scale,*vx;
860: VecGetSize(nulls[0],&N);
861: scale = 2.0/sqrt((PetscReal)N);
862: VecGetArray(nulls[0],&vx);
863: VecGetOwnershipRange(nulls[0],&rstart,&rend);
864: for (i=rstart; i<rend; i++) {
865: if (!(i % 4)) vx[i-rstart] = scale;
866: else vx[i-rstart] = 0.0;
867: }
868: VecRestoreArray(nulls[0],&vx);
869: return(0);
870: }
872: /*
873: This is an experimental function and can be safely ignored.
874: */
875: PetscErrorCode FormFunctionLocali(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
876: {
877: AppCtx *user = (AppCtx*)ptr;
878: TstepCtx *tsCtx = user->tsCtx;
879: Parameters *param = user->param;
881: PetscInt i,j,c;
882: PetscInt xints,xinte,yints,yinte;
883: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
884: PassiveReal de2,rhos2,nu,eta,dde2;
885: PassiveReal two = 2.0,one = 1.0,p5 = 0.5;
886: PassiveReal F_eq_x,By_eq,dtinv;
887: PetscScalar xx;
888: PetscScalar vx,vy,avx,avy,vxp,vxm,vyp,vym;
889: PetscScalar Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
890: PassiveField **xold;
893: de2 = sqr(param->d_e);
894: rhos2 = sqr(param->rho_s);
895: nu = param->nu;
896: eta = param->eta;
897: dde2 = one/de2;
899: /*
900: Define mesh intervals ratios for uniform grid.
901: [Note: FD formulae below are normalized by multiplying through by
902: local volume element to obtain coefficients O(1) in two dimensions.]
903: */
904: dhx = info->mx/lx; dhy = info->my/ly;
905: hx = one/dhx; hy = one/dhy;
906: hxdhy = hx*dhy; hydhx = hy*dhx;
907: hxhy = hx*hy; dhxdhy = dhx*dhy;
909: xints = info->xs; xinte = info->xs+info->xm;
910: yints = info->ys; yinte = info->ys+info->ym;
913: i = st->i; j = st->j; c = st->c;
915: #ifdef EQ
916: xx = i * hx;
917: F_eq_x = - (1. + de2) * sin(PetscAbsScalar(xx));
918: By_eq = sin(PetscAbsScalar(xx));
919: #else
920: F_eq_x = 0.;
921: By_eq = 0.;
922: #endif
924: /*
925: * convective coefficients for upwinding
926: */
928: vx = - D_y(x,phi,i,j);
929: vy = D_x(x,phi,i,j);
930: avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
931: avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
932: #ifndef UPWINDING
933: vxp = vxm = p5*vx;
934: vyp = vym = p5*vy;
935: #endif
937: Bx = D_y(x,psi,i,j);
938: By = - D_x(x,psi,i,j) + By_eq;
939: aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
940: aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
941: #ifndef UPWINDING
942: Bxp = Bxm = p5*Bx;
943: Byp = Bym = p5*By;
944: #endif
946: DAVecGetArray(info->da,user->Xold,&xold);
947: dtinv = hxhy/(tsCtx->dt);
948: switch(c) {
950: case 0:
951: /* Lap(phi) - U */
952: *f = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
953: break;
955: case 1:
956: /* psi - d_e^2 * Lap(psi) - F */
957: *f = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
958: break;
960: case 2:
961: /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
962: - nu Lap(U) */
963: *f =
964: ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
965: vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
966: (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
967: Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
968: nu * Lapl(x,U,i,j)) * hxhy;
969: *f += dtinv*(x[j][i].U-xold[j][i].U);
970: break;
972: case 3:
973: /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
974: - eta * Lap(psi) */
975: *f =
976: ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
977: vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
978: (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
979: Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
980: eta * Lapl(x,psi,i,j)) * hxhy;
981: *f += dtinv*(x[j][i].F-xold[j][i].F);
982: break;
983: }
984: DAVecRestoreArray(info->da,user->Xold,&xold);
987: /*
988: Flop count (multiply-adds are counted as 2 operations)
989: */
990: /* PetscLogFlops(84*info->ym*info->xm); FIXME */
991: return(0);
992: }
993: /*
994: This is an experimental function and can be safely ignored.
995: */
996: PetscErrorCode FormFunctionLocali4(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *ff,void *ptr)
997: {
998: Field *f = (Field *)ff;
999: AppCtx *user = (AppCtx*)ptr;
1000: TstepCtx *tsCtx = user->tsCtx;
1001: Parameters *param = user->param;
1003: PetscInt i,j;
1004: PetscInt xints,xinte,yints,yinte;
1005: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
1006: PassiveReal de2,rhos2,nu,eta,dde2;
1007: PassiveReal two = 2.0,one = 1.0,p5 = 0.5;
1008: PassiveReal F_eq_x,By_eq,dtinv;
1009: PetscScalar xx;
1010: PetscScalar vx,vy,avx,avy,vxp,vxm,vyp,vym;
1011: PetscScalar Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
1012: PassiveField **xold;
1015: de2 = sqr(param->d_e);
1016: rhos2 = sqr(param->rho_s);
1017: nu = param->nu;
1018: eta = param->eta;
1019: dde2 = one/de2;
1021: /*
1022: Define mesh intervals ratios for uniform grid.
1023: [Note: FD formulae below are normalized by multiplying through by
1024: local volume element to obtain coefficients O(1) in two dimensions.]
1025: */
1026: dhx = info->mx/lx; dhy = info->my/ly;
1027: hx = one/dhx; hy = one/dhy;
1028: hxdhy = hx*dhy; hydhx = hy*dhx;
1029: hxhy = hx*hy; dhxdhy = dhx*dhy;
1031: xints = info->xs; xinte = info->xs+info->xm;
1032: yints = info->ys; yinte = info->ys+info->ym;
1035: i = st->i; j = st->j;
1037: #ifdef EQ
1038: xx = i * hx;
1039: F_eq_x = - (1. + de2) * sin(PetscAbsScalar(xx));
1040: By_eq = sin(PetscAbsScalar(xx));
1041: #else
1042: F_eq_x = 0.;
1043: By_eq = 0.;
1044: #endif
1046: /*
1047: * convective coefficients for upwinding
1048: */
1050: vx = - D_y(x,phi,i,j);
1051: vy = D_x(x,phi,i,j);
1052: avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
1053: avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
1054: #ifndef UPWINDING
1055: vxp = vxm = p5*vx;
1056: vyp = vym = p5*vy;
1057: #endif
1059: Bx = D_y(x,psi,i,j);
1060: By = - D_x(x,psi,i,j) + By_eq;
1061: aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
1062: aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
1063: #ifndef UPWINDING
1064: Bxp = Bxm = p5*Bx;
1065: Byp = Bym = p5*By;
1066: #endif
1068: DAVecGetArray(info->da,user->Xold,&xold);
1069: dtinv = hxhy/(tsCtx->dt);
1070:
1072:
1073: /* Lap(phi) - U */
1074: f->phi = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
1075:
1077:
1078: /* psi - d_e^2 * Lap(psi) - F */
1079: f->psi = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
1080:
1082:
1083: /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
1084: - nu Lap(U) */
1085: f->U =
1086: ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
1087: vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
1088: (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
1089: Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
1090: nu * Lapl(x,U,i,j)) * hxhy;
1091: f->U += dtinv*(x[j][i].U-xold[j][i].U);
1092:
1094:
1095: /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
1096: - eta * Lap(psi) */
1097: f->F =
1098: ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
1099: vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
1100: (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
1101: Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
1102: eta * Lapl(x,psi,i,j)) * hxhy;
1103: f->F += dtinv*(x[j][i].F-xold[j][i].F);
1104:
1105:
1106: DAVecRestoreArray(info->da,user->Xold,&xold);
1109: /*
1110: Flop count (multiply-adds are counted as 2 operations)
1111: */
1112: /* PetscLogFlops(84*info->ym*info->xm); FIXME */
1113: return(0);
1114: }