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: DMDA^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 "petscdmda.h" so that we can use distributed arrays (DMDAs).
34: Include "petscsnes.h" so that we can use SNES solvers.
35: Include "petscpcmg.h" to control the multigrid solvers.
36: Note that these automatically include:
37: petscsys.h - base PETSc routines petscvec.h - vectors
38: 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 <petscdmda.h>
45: #include <petscpcmg.h>
46: #include <petscdmmg.h>
48: #ifdef HAVE_DA_HDF
49: PetscInt DMDAVecHDFOutput(DM,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: PetscReal fnorm;
80: PetscBool 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: PetscBool draw_contours; /* flag - 1 indicates drawing contours */
98: PetscBool second_order;
99: PetscBool PetscPreLoading;
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: DM 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: PetscPreLoadBegin(PETSC_TRUE,"SetUp");
144: param.PetscPreLoading = PetscPreLoading;
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: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_PERIODIC, DMDA_BOUNDARY_PERIODIC, DMDA_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: DMDASetBlockFills(da,dfill,ofill);
160: DMMGSetDM(dmmg,(DM)da);
161: DMDestroy(&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: DMDASetFieldName(DMMGGetDM(dmmg), 0, "phi");
184: DMDASetFieldName(DMMGGetDM(dmmg), 1, "psi");
186: DMDASetFieldName(DMMGGetDM(dmmg), 2, "U");
188: DMDASetFieldName(DMMGGetDM(dmmg), 3, "F");
190: /*======================================================================*/
191: /* Initialize stuff related to time stepping */
192: /*======================================================================*/
194: DMDAGetInfo(DMMGGetDM(dmmg),0,&mx,&my,0,0,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: PetscBool flg;
224: PetscOptionsHasName(PETSC_NULL, "-socket_viewer", &flg);
225: if (flg && !PetscPreLoading) {
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: DMMGSetFromOptions(dmmg);
256: DMMGSetSNESLocali(dmmg,FormFunctionLocali,0,admf_FormFunctionLocali);
257: DMMGSetSNESLocalib(dmmg,FormFunctionLocali4,0,admfb_FormFunctionLocali4);
258:
259: /* attach nullspace to each level of the preconditioner */
260: DMMGSetNullSpace(dmmg,PETSC_FALSE,1,CreateNullSpace);
261:
262: PetscPrintf(PETSC_COMM_WORLD, "finish setupNull!");
264: if (PetscPreLoading) {
265: PetscPrintf(PETSC_COMM_WORLD, "# viscosity = %G, resistivity = %G, "
266: "skin_depth # = %G, larmor_radius # = %G\n",
267: param.nu, param.eta, param.d_e, param.rho_s);
268: DMDAGetInfo(DMMGGetDM(dmmg),0,&m,&n,0,0,0,0,0,0,0,0,0,0);
269: PetscPrintf(PETSC_COMM_WORLD,"Problem size %D by %D\n",m,n);
270: 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);
271: }
273:
274:
275: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
276: Solve the nonlinear system
277: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
278:
279: PetscPreLoadStage("Solve");
281: if (param.draw_contours) {
282: VecView(((AppCtx*)DMMGGetUser(dmmg,param.mglevels-1))->Xold,PETSC_VIEWER_DRAW_WORLD);
283: }
285: Update(dmmg);
287: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
288: Free work space. All PETSc objects should be destroyed when they
289: are no longer needed.
290: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
291:
292: for (i=0; i<param.mglevels; i++) {
293: VecDestroy(&user[i].Xoldold);
294: VecDestroy(&user[i].Xold);
295: }
296: PetscFree(user);
297: DMMGDestroy(dmmg);
299: PetscPreLoadEnd();
300:
301: PetscFinalize();
302: return 0;
303: }
305: /* ------------------------------------------------------------------- */
308: /* ------------------------------------------------------------------- */
309: PetscErrorCode Gnuplot(DM da, Vec X, double mtime)
310: {
311: PetscInt i,j,xs,ys,xm,ym;
312: PetscInt xints,xinte,yints,yinte;
314: Field **x;
315: FILE *f;
316: char fname[PETSC_MAX_PATH_LEN];
317: PetscMPIInt rank;
319: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
321: sprintf(fname, "out-%g-%d.dat", mtime, rank);
323: f = fopen(fname, "w");
325: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
327: DMDAVecGetArray(da,X,&x);
329: xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;
331: for (j=yints; j<yinte; j++) {
332: for (i=xints; i<xinte; i++) {
333: PetscFPrintf(PETSC_COMM_WORLD, f,
334: "%D %D %G %G %G %G %G %G\n",
335: i, j, 0.0, 0.0,
336: PetscAbsScalar(x[j][i].U), PetscAbsScalar(x[j][i].F),
337: PetscAbsScalar(x[j][i].phi), PetscAbsScalar(x[j][i].psi));
338: }
339: PetscFPrintf(PETSC_COMM_WORLD,f, "\n");
340: }
341: DMDAVecRestoreArray(da,X,&x);
342: fclose(f);
343: return 0;
344: }
346: /* ------------------------------------------------------------------- */
349: /* ------------------------------------------------------------------- */
350: PetscErrorCode Initialize(DMMG *dmmg)
351: {
352: AppCtx *appCtx = (AppCtx*)DMMGGetUser(dmmg,0);
353: Parameters *param = appCtx->param;
354: DM da;
356: PetscInt i,j,mx,my,xs,ys,xm,ym;
357: PetscReal two = 2.0,one = 1.0;
358: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx,dhxdhy; /* hxhy */
359: PetscReal d_e,rho_s,de2,xx,yy;
360: Field **x, **localx;
361: Vec localX;
362: PetscBool flg;
365: PetscOptionsHasName(0,"-restart",&flg);
366: if (flg) {
367: PetscViewer viewer;
368: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"binaryoutput",FILE_MODE_READ,&viewer);
369: VecLoad(dmmg[param->mglevels-1]->x,viewer);
370: PetscViewerDestroy(&viewer);
371: return(0);
372: }
374: d_e = param->d_e;
375: rho_s = param->rho_s;
376: de2 = sqr(param->d_e);
378: da = (dmmg[param->mglevels-1]->dm);
379: DMDAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0,0,0);
381: dhx = mx/lx; dhy = my/ly;
382: hx = one/dhx; hy = one/dhy;
383: hxdhy = hx*dhy; hydhx = hy*dhx;
384: /* hxhy = hx*hy; */ dhxdhy = dhx*dhy;
386: /*
387: Get local grid boundaries (for 2-dimensional DMDA):
388: xs, ys - starting grid indices (no ghost points)
389: xm, ym - widths of local grid (no ghost points)
390: */
391: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
393: DMGetLocalVector(da,&localX);
394: /*
395: Get a pointer to vector data.
396: - For default PETSc vectors, VecGetArray() returns a pointer to
397: the data array. Otherwise, the routine is implementation dependent.
398: - You MUST call VecRestoreArray() when you no longer need access to
399: the array.
400: */
401: DMDAVecGetArray(da,dmmg[param->mglevels-1]->x,&x);
402: DMDAVecGetArray(da,localX,&localx);
404: /*
405: Compute initial solution over the locally owned part of the grid
406: */
407: #if defined (PETSC_HAVE_ERF)
408: {
409: PetscReal eps = lx/ly;
410: PetscReal pert = 1e-4;
411: PetscReal k = 1.*eps;
412: PetscReal gam;
414: if (d_e < rho_s) d_e = rho_s;
415: gam = k * d_e;
416: for (j=ys-1; j<ys+ym+1; j++) {
417: yy = j * hy;
418: for (i=xs-1; i<xs+xm+1; i++) {
419: xx = i * hx;
420: if (xx < -3.1416926/2) {
421: localx[j][i].phi = pert * gam / k * erf((xx + 3.1415926) / (sqrt(2.0) * d_e)) * (-sin(k*yy));
422: } else if (xx < 3.1415926/2) {
423: localx[j][i].phi = - pert * gam / k * erf(xx / (sqrt(2.0) * d_e)) * (-sin(k*yy));
424: } else if (xx < 3*3.1415926/2){
425: localx[j][i].phi = pert * gam / k * erf((xx - 3.1415926) / (sqrt(2.0) * d_e)) * (-sin(k*yy));
426: } else {
427: localx[j][i].phi = - pert * gam / k * erf((xx - 2.*3.1415926) / (sqrt(2.0) * d_e)) * (-sin(k*yy));
428: }
429: #ifdef EQ
430: localx[j][i].psi = 0.;
431: #else
432: localx[j][i].psi = cos(xx);
433: #endif
434: }
435: }
437: for (j=ys; j<ys+ym; j++) {
438: for (i=xs; i<xs+xm; i++) {
439: x[j][i].U = Lapl(localx,phi,i,j);
440: x[j][i].F = localx[j][i].psi - de2 * Lapl(localx,psi,i,j);
441: x[j][i].phi = localx[j][i].phi;
442: x[j][i].psi = localx[j][i].psi;
443: }
444: }
445: }
446: #else
447: SETERRQ(PETSC_COMM_SELF,1,"erf() not available on this machine");
448: #endif
450: /*
451: Restore vector
452: */
453: DMDAVecRestoreArray(da,dmmg[param->mglevels-1]->x,&x);
454: DMDAVecRestoreArray(da,localX,&localx);
455: DMRestoreLocalVector(da,&localX);
457: return(0);
458: }
462: PetscErrorCode ComputeMaxima(DM da, Vec X, PetscReal t)
463: {
465: PetscInt i,j,mx,my,xs,ys,xm,ym;
466: PetscInt xints,xinte,yints,yinte;
467: Field **x;
468: double norm[4] = {0,0,0,0};
469: double gnorm[4];
470: MPI_Comm comm;
474: DMDAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0,0,0);
476: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
478: xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;
480: DMDAVecGetArray(da, X, &x);
482: for (j=yints; j<yinte; j++) {
483: for (i=xints; i<xinte; i++) {
484: norm[0] = PetscMax(norm[0],PetscAbsScalar(x[j][i].U));
485: norm[1] = PetscMax(norm[1],PetscAbsScalar(x[j][i].F));
486: norm[2] = PetscMax(norm[2],PetscAbsScalar(x[j][i].phi));
487: norm[3] = PetscMax(norm[3],PetscAbsScalar(x[j][i].psi));
488: }
489: }
491: DMDAVecRestoreArray(da,X,&x);
493: PetscObjectGetComm((PetscObject)da, &comm);
494: MPI_Allreduce(norm, gnorm, 4, MPI_DOUBLE, MPI_MAX, comm);
496: PetscFPrintf(PETSC_COMM_WORLD, stderr,"%G\t%G\t%G\t%G\t%G\n",t, gnorm[0], gnorm[1], gnorm[2], gnorm[3]);
498: return(0);
499: }
503: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
504: {
505: AppCtx *user = (AppCtx*)ptr;
506: TstepCtx *tsCtx = user->tsCtx;
507: Parameters *param = user->param;
509: PetscInt xints,xinte,yints,yinte,i,j;
510: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
511: PassiveReal de2,rhos2,nu,eta,dde2;
512: PassiveReal two = 2.0,one = 1.0,p5 = 0.5;
513: PassiveReal F_eq_x,By_eq;
514: PetscScalar xx;
515: PetscScalar vx,vy,avx,avy,vxp,vxm,vyp,vym;
516: PetscScalar Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
519: de2 = sqr(param->d_e);
520: rhos2 = sqr(param->rho_s);
521: nu = param->nu;
522: eta = param->eta;
523: dde2 = one/de2;
525: /*
526: Define mesh intervals ratios for uniform grid.
527: [Note: FD formulae below are normalized by multiplying through by
528: local volume element to obtain coefficients O(1) in two dimensions.]
529: */
530: dhx = info->mx/lx; dhy = info->my/ly;
531: hx = one/dhx; hy = one/dhy;
532: hxdhy = hx*dhy; hydhx = hy*dhx;
533: hxhy = hx*hy; dhxdhy = dhx*dhy;
535: xints = info->xs; xinte = info->xs+info->xm;
536: yints = info->ys; yinte = info->ys+info->ym;
538: /* Compute over the interior points */
539: for (j=yints; j<yinte; j++) {
540: for (i=xints; i<xinte; i++) {
541: #ifdef EQ
542: xx = i * hx;
543: F_eq_x = - (1. + de2) * sin(PetscAbsScalar(xx));
544: By_eq = sin(PetscAbsScalar(xx));
545: #else
546: F_eq_x = 0.;
547: By_eq = 0.;
548: #endif
550: /*
551: * convective coefficients for upwinding
552: */
554: vx = - D_y(x,phi,i,j);
555: vy = D_x(x,phi,i,j);
556: avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
557: avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
558: #ifndef UPWINDING
559: vxp = vxm = p5*vx;
560: vyp = vym = p5*vy;
561: #endif
563: Bx = D_y(x,psi,i,j);
564: By = - D_x(x,psi,i,j) + By_eq;
565: aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
566: aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
567: #ifndef UPWINDING
568: Bxp = Bxm = p5*Bx;
569: Byp = Bym = p5*By;
570: #endif
572: /* Lap(phi) - U */
573: f[j][i].phi = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
575: /* psi - d_e^2 * Lap(psi) - F */
576: f[j][i].psi = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
578: /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
579: - nu Lap(U) */
580: f[j][i].U =
581: ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
582: vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
583: (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
584: Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
585: nu * Lapl(x,U,i,j)) * hxhy;
586:
587: /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
588: - eta * Lap(psi) */
589: f[j][i].F =
590: ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
591: vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
592: (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
593: Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
594: eta * Lapl(x,psi,i,j)) * hxhy;
595: }
596: }
598: /* Add time step contribution */
599: if (tsCtx->ires) {
600: if ((param->second_order) && (tsCtx->itstep > 0)){
601: AddTSTermLocal2(info,x,f,user);
602: } else {
603: AddTSTermLocal(info,x,f,user);
604: }
605: }
607: /*
608: Flop count (multiply-adds are counted as 2 operations)
609: */
610: /* PetscLogFlops(84.0*info->ym*info->xm); FIXME */
611: return(0);
612: }
614: /*---------------------------------------------------------------------*/
617: PetscErrorCode Update(DMMG *dmmg)
618: /*---------------------------------------------------------------------*/
619: {
620: AppCtx *user = (AppCtx *) DMMGGetUser(dmmg,0);
621: TstepCtx *tsCtx = user->tsCtx;
622: Parameters *param = user->param;
623: SNES snes;
624: PetscErrorCode ierr;
625: PetscInt its,lits,i;
626: PetscInt max_steps;
627: PetscInt nfailsCum = 0,nfails = 0;
628: static PetscInt ic_out;
629: PetscBool ts_monitor = (tsCtx->ts_monitor && !param->PetscPreLoading) ? PETSC_TRUE : PETSC_FALSE;
633: if (param->PetscPreLoading)
634: max_steps = 1;
635: else
636: max_steps = tsCtx->max_steps;
637:
638: Initialize(dmmg);
640: snes = DMMGGetSNES(dmmg);
642: for (tsCtx->itstep = 0; tsCtx->itstep < max_steps; tsCtx->itstep++) {
643:
644: PetscPrintf(PETSC_COMM_WORLD, "time step=%d!\n",tsCtx->itstep);
645: if ((param->second_order) && (tsCtx->itstep > 0))
646: {
647: for (i=param->mglevels-1; i>=0 ;i--)
648: {
649: VecCopy(((AppCtx*)DMMGGetUser(dmmg,i))->Xold,((AppCtx*)DMMGGetUser(dmmg,i))->Xoldold);
650: }
651: }
653: for (i=param->mglevels-1; i>0 ;i--) {
654: MatRestrict(dmmg[i]->R, dmmg[i]->x, dmmg[i-1]->x);
655: VecPointwiseMult(dmmg[i-1]->x,dmmg[i-1]->x,dmmg[i]->Rscale);
656: VecCopy(dmmg[i]->x, ((AppCtx*)DMMGGetUser(dmmg,i))->Xold);
657: }
658: VecCopy(dmmg[0]->x, ((AppCtx*)DMMGGetUser(dmmg,0))->Xold);
660: DMMGSolve(dmmg);
664: if (tsCtx->itstep == 665000)
665: {
666: KSP ksp;
667: PC pc;
668: Mat mat, pmat;
669: MatStructure flag;
670: PetscViewer viewer;
671: char file[PETSC_MAX_PATH_LEN];
673: PetscStrcpy(file, "matrix");
675: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_WRITE, &viewer);
677: SNESGetKSP(snes, &ksp);
679: KSPGetPC(ksp, &pc);
681: PCGetOperators(pc, &mat, &pmat, &flag);
683: MatView(mat, viewer);
685: PetscViewerDestroy(&viewer);
686: SETERRQ(PETSC_COMM_SELF,1,"Done saving Jacobian");
687: }
690: tsCtx->t += tsCtx->dt;
692: /* save restart solution if requested at a particular time, then exit */
693: if (tsCtx->dump_time > 0.0 && tsCtx->t >= tsCtx->dump_time) {
694: Vec v = DMMGGetx(dmmg);
695: VecView(v,PETSC_VIEWER_BINARY_WORLD);
696: SETERRQ1(PETSC_COMM_SELF,1,"Saved solution at time %G",tsCtx->t);
697: }
699: if (ts_monitor)
700: {
701: SNESGetIterationNumber(snes, &its);
702: SNESGetLinearSolveIterations(snes, &lits);
703: SNESGetNonlinearStepFailures(snes, &nfails);
704: SNESGetFunctionNorm(snes, &tsCtx->fnorm);
706: nfailsCum += nfails;
707: if (nfailsCum >= 2) SETERRQ(PETSC_COMM_SELF,1, "unable to find a newton step");
709: PetscPrintf(PETSC_COMM_WORLD,
710: "time step = %D, time = %G, number of nonlinear steps = %D, "
711: "number of linear steps = %D, norm of the function = %G\n",
712: tsCtx->itstep + 1, tsCtx->t, its, lits, PetscAbsScalar(tsCtx->fnorm));
714: /* send solution over to MATLAB, to be visualized (using ex29.m) */
715: if (!param->PetscPreLoading && tsCtx->socketviewer)
716: {
717: Vec v;
718: SNESGetSolution(snes, &v);
719: #if defined(PETSC_USE_SOCKET_VIEWER)
720: VecView(v, tsCtx->socketviewer);
721: #endif
722: }
723: }
725: if (!param->PetscPreLoading) {
726: if (param->draw_contours) {
727: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
728: }
730: if (1 && ts_monitor) {
731: /* compute maxima */
732: ComputeMaxima( dmmg[param->mglevels-1]->dm,dmmg[param->mglevels-1]->x,tsCtx->t);
733: /* output */
734: if (ic_out++ == (int)(tsCtx->dt_out / tsCtx->dt + .5)) {
735: #ifdef HAVE_DA_HDF
736: char fname[PETSC_MAX_PATH_LEN];
738: sprintf(fname, "out-%G.hdf", tsCtx->t);
739: DMDAVecHDFOutput(DMMGGetDM(dmmg), DMMGGetx(dmmg), fname);
740: #else
741: /*
742: Gnuplot(DMMGGetDM(dmmg), DMMGGetx(dmmg), tsCtx->t);
743: */
744: #endif
745: ic_out = 1;
746: }
747: }
748: }
749: } /* End of time step loop */
750:
751: if (!param->PetscPreLoading){
752: SNESGetFunctionNorm(snes,&tsCtx->fnorm);
753: PetscPrintf(PETSC_COMM_WORLD, "timesteps %D fnorm = %G\n",
754: tsCtx->itstep, PetscAbsScalar(tsCtx->fnorm));
755: }
757: if (param->PetscPreLoading) {
758: tsCtx->fnorm_ini = 0.0;
759: }
760:
761: return(0);
762: }
764: /*---------------------------------------------------------------------*/
767: PetscErrorCode AddTSTermLocal(DMDALocalInfo* info,Field **x,Field **f,AppCtx *user)
768: /*---------------------------------------------------------------------*/
769: {
770: TstepCtx *tsCtx = user->tsCtx;
771: DM da = info->da;
773: PetscInt i,j;
774: PetscInt xints,xinte,yints,yinte;
775: PassiveReal hx,hy,dhx,dhy,hxhy;
776: PassiveReal one = 1.0;
777: PassiveScalar dtinv;
778: PassiveField **xold;
782: xints = info->xs; xinte = info->xs+info->xm;
783: yints = info->ys; yinte = info->ys+info->ym;
785: dhx = info->mx/lx; dhy = info->my/ly;
786: hx = one/dhx; hy = one/dhy;
787: hxhy = hx*hy;
789: dtinv = hxhy/(tsCtx->dt);
791: DMDAVecGetArray(da,user->Xold,&xold);
792: for (j=yints; j<yinte; j++) {
793: for (i=xints; i<xinte; i++) {
794: f[j][i].U += dtinv*(x[j][i].U-xold[j][i].U);
795: f[j][i].F += dtinv*(x[j][i].F-xold[j][i].F);
796: }
797: }
798: DMDAVecRestoreArray(da,user->Xold,&xold);
800: return(0);
801: }
803: /*---------------------------------------------------------------------*/
806: PetscErrorCode AddTSTermLocal2(DMDALocalInfo* info,Field **x,Field **f,AppCtx *user)
807: /*---------------------------------------------------------------------*/
808: {
809: TstepCtx *tsCtx = user->tsCtx;
810: DM da = info->da;
812: PetscInt i,j,xints,xinte,yints,yinte;
813: PassiveReal hx,hy,dhx,dhy,hxhy;
814: PassiveReal one = 1.0, onep5 = 1.5, two = 2.0, p5 = 0.5;
815: PassiveScalar dtinv;
816: PassiveField **xoldold,**xold;
820: xints = info->xs; xinte = info->xs+info->xm;
821: yints = info->ys; yinte = info->ys+info->ym;
823: dhx = info->mx/lx; dhy = info->my/ly;
824: hx = one/dhx; hy = one/dhy;
825: hxhy = hx*hy;
827: dtinv = hxhy/(tsCtx->dt);
829: DMDAVecGetArray(da,user->Xoldold,&xoldold);
830: DMDAVecGetArray(da,user->Xold,&xold);
831: for (j=yints; j<yinte; j++) {
832: for (i=xints; i<xinte; i++) {
833: f[j][i].U += dtinv * (onep5 * x[j][i].U - two * xold[j][i].U +
834: p5 * xoldold[j][i].U);
835: f[j][i].F += dtinv * (onep5 * x[j][i].F - two * xold[j][i].F +
836: p5 * xoldold[j][i].F);
837: }
838: }
839: DMDAVecRestoreArray(da,user->Xoldold,&xoldold);
840: DMDAVecRestoreArray(da,user->Xold,&xold);
842: return(0);
843: }
845: /* Creates the null space of the Jacobian for a particular level */
848: PetscErrorCode CreateNullSpace(DMMG dmmg,Vec *nulls)
849: {
851: PetscInt i,N,rstart,rend;
852: PetscScalar scale,*vx;
855: VecGetSize(nulls[0],&N);
856: scale = 2.0/sqrt((PetscReal)N);
857: VecGetArray(nulls[0],&vx);
858: VecGetOwnershipRange(nulls[0],&rstart,&rend);
859: for (i=rstart; i<rend; i++) {
860: if (!(i % 4)) vx[i-rstart] = scale;
861: else vx[i-rstart] = 0.0;
862: }
863: VecRestoreArray(nulls[0],&vx);
864: return(0);
865: }
867: /*
868: This is an experimental function and can be safely ignored.
869: */
870: PetscErrorCode FormFunctionLocali(DMDALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
871: {
872: AppCtx *user = (AppCtx*)ptr;
873: TstepCtx *tsCtx = user->tsCtx;
874: Parameters *param = user->param;
876: PetscInt i,j,c;
877: /* PetscInt xints,xinte,yints,yinte; */
878: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
879: PassiveReal de2,rhos2,nu,eta,dde2;
880: PassiveReal two = 2.0,one = 1.0,p5 = 0.5;
881: PassiveReal F_eq_x,By_eq,dtinv;
882: PetscScalar xx;
883: PetscScalar vx,vy,avx,avy,vxp,vxm,vyp,vym;
884: PetscScalar Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
885: PassiveField **xold;
888: de2 = sqr(param->d_e);
889: rhos2 = sqr(param->rho_s);
890: nu = param->nu;
891: eta = param->eta;
892: dde2 = one/de2;
894: /*
895: Define mesh intervals ratios for uniform grid.
896: [Note: FD formulae below are normalized by multiplying through by
897: local volume element to obtain coefficients O(1) in two dimensions.]
898: */
899: dhx = info->mx/lx; dhy = info->my/ly;
900: hx = one/dhx; hy = one/dhy;
901: hxdhy = hx*dhy; hydhx = hy*dhx;
902: hxhy = hx*hy; dhxdhy = dhx*dhy;
904: /*
905: xints = info->xs; xinte = info->xs+info->xm;
906: yints = info->ys; yinte = info->ys+info->ym;
907: */
909: i = st->i; j = st->j; c = st->c;
911: #ifdef EQ
912: xx = i * hx;
913: F_eq_x = - (1. + de2) * sin(PetscAbsScalar(xx));
914: By_eq = sin(PetscAbsScalar(xx));
915: #else
916: F_eq_x = 0.;
917: By_eq = 0.;
918: #endif
920: /*
921: * convective coefficients for upwinding
922: */
924: vx = - D_y(x,phi,i,j);
925: vy = D_x(x,phi,i,j);
926: avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
927: avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
928: #ifndef UPWINDING
929: vxp = vxm = p5*vx;
930: vyp = vym = p5*vy;
931: #endif
933: Bx = D_y(x,psi,i,j);
934: By = - D_x(x,psi,i,j) + By_eq;
935: aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
936: aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
937: #ifndef UPWINDING
938: Bxp = Bxm = p5*Bx;
939: Byp = Bym = p5*By;
940: #endif
942: DMDAVecGetArray(info->da,user->Xold,&xold);
943: dtinv = hxhy/(tsCtx->dt);
944: switch(c) {
946: case 0:
947: /* Lap(phi) - U */
948: *f = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
949: break;
951: case 1:
952: /* psi - d_e^2 * Lap(psi) - F */
953: *f = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
954: break;
956: case 2:
957: /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
958: - nu Lap(U) */
959: *f =
960: ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
961: vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
962: (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
963: Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
964: nu * Lapl(x,U,i,j)) * hxhy;
965: *f += dtinv*(x[j][i].U-xold[j][i].U);
966: break;
968: case 3:
969: /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
970: - eta * Lap(psi) */
971: *f =
972: ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
973: vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
974: (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
975: Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
976: eta * Lapl(x,psi,i,j)) * hxhy;
977: *f += dtinv*(x[j][i].F-xold[j][i].F);
978: break;
979: }
980: DMDAVecRestoreArray(info->da,user->Xold,&xold);
983: /*
984: Flop count (multiply-adds are counted as 2 operations)
985: */
986: /* PetscLogFlops(84.0*info->ym*info->xm); FIXME */
987: return(0);
988: }
989: /*
990: This is an experimental function and can be safely ignored.
991: */
992: PetscErrorCode FormFunctionLocali4(DMDALocalInfo *info,MatStencil *st,Field **x,PetscScalar *ff,void *ptr)
993: {
994: Field *f = (Field *)ff;
995: AppCtx *user = (AppCtx*)ptr;
996: TstepCtx *tsCtx = user->tsCtx;
997: Parameters *param = user->param;
999: PetscInt i,j;
1000: /* PetscInt xints,xinte,yints,yinte; */
1001: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
1002: PassiveReal de2,rhos2,nu,eta,dde2;
1003: PassiveReal two = 2.0,one = 1.0,p5 = 0.5;
1004: PassiveReal F_eq_x,By_eq,dtinv;
1005: PetscScalar xx;
1006: PetscScalar vx,vy,avx,avy,vxp,vxm,vyp,vym;
1007: PetscScalar Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
1008: PassiveField **xold;
1011: de2 = sqr(param->d_e);
1012: rhos2 = sqr(param->rho_s);
1013: nu = param->nu;
1014: eta = param->eta;
1015: dde2 = one/de2;
1017: /*
1018: Define mesh intervals ratios for uniform grid.
1019: [Note: FD formulae below are normalized by multiplying through by
1020: local volume element to obtain coefficients O(1) in two dimensions.]
1021: */
1022: dhx = info->mx/lx; dhy = info->my/ly;
1023: hx = one/dhx; hy = one/dhy;
1024: hxdhy = hx*dhy; hydhx = hy*dhx;
1025: hxhy = hx*hy; dhxdhy = dhx*dhy;
1027: /*
1028: xints = info->xs; xinte = info->xs+info->xm;
1029: yints = info->ys; yinte = info->ys+info->ym;
1030: */
1032: i = st->i; j = st->j;
1034: #ifdef EQ
1035: xx = i * hx;
1036: F_eq_x = - (1. + de2) * sin(PetscAbsScalar(xx));
1037: By_eq = sin(PetscAbsScalar(xx));
1038: #else
1039: F_eq_x = 0.;
1040: By_eq = 0.;
1041: #endif
1043: /*
1044: * convective coefficients for upwinding
1045: */
1047: vx = - D_y(x,phi,i,j);
1048: vy = D_x(x,phi,i,j);
1049: avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
1050: avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
1051: #ifndef UPWINDING
1052: vxp = vxm = p5*vx;
1053: vyp = vym = p5*vy;
1054: #endif
1056: Bx = D_y(x,psi,i,j);
1057: By = - D_x(x,psi,i,j) + By_eq;
1058: aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
1059: aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
1060: #ifndef UPWINDING
1061: Bxp = Bxm = p5*Bx;
1062: Byp = Bym = p5*By;
1063: #endif
1065: DMDAVecGetArray(info->da,user->Xold,&xold);
1066: dtinv = hxhy/(tsCtx->dt);
1067:
1069:
1070: /* Lap(phi) - U */
1071: f->phi = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
1072:
1074:
1075: /* psi - d_e^2 * Lap(psi) - F */
1076: f->psi = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
1077:
1079:
1080: /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
1081: - nu Lap(U) */
1082: f->U =
1083: ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
1084: vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
1085: (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
1086: Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
1087: nu * Lapl(x,U,i,j)) * hxhy;
1088: f->U += dtinv*(x[j][i].U-xold[j][i].U);
1089:
1091:
1092: /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
1093: - eta * Lap(psi) */
1094: f->F =
1095: ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
1096: vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
1097: (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
1098: Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
1099: eta * Lapl(x,psi,i,j)) * hxhy;
1100: f->F += dtinv*(x[j][i].F-xold[j][i].F);
1101:
1102:
1103: DMDAVecRestoreArray(info->da,user->Xold,&xold);
1106: /*
1107: Flop count (multiply-adds are counted as 2 operations)
1108: */
1109: /* PetscLogFlops(84.0*info->ym*info->xm); FIXME */
1110: return(0);
1111: }